Re: [R] gam confidence interval (package mgcv)
not sure if I'm missing something here, but since you are using a log link, isn't the ratio you are looking for given by the `treatmentB' parameter in the summary (independent of X) > summary(gfit) [snip] Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.831830.03885 47.152 < 2e-16 *** treatmentB 0.445130.05567 7.996 1.4e-09 *** --- [snip] Let mu = E(y), and b be a binary indicator for treatment B. You want mu|b=1/mu|b=0 log(mu|b=1) = intercept + treatmentB + s(X) log(mu|b=0) = intercept + s(X) => log(mu|b=1) - log(mu|b=0) = treatmentB so mu|b=1/mu|b = exp(treatmentB) So you can get the required interval by finding and interval for treatment B and exponentiating... tB <- coef(gfit)[2] se.tB <- sqrt(vcov(gfit)[2,2]) exp(c(tB - 2*se.tB,tB+2*se.tB)) On 06/28/2011 03:45 AM, Remko Duursma wrote: Dear R-helpers, I am trying to construct a confidence interval on a prediction of a gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant but I am not able to apply that to this, different, problem. Any help is appreciated! Basically I have a function Y = f(X) for two different treatments A and B. I am interested in the treatment ratios : Y(treatment = B) / Y(treatment = A) as a function of X, including a confidence interval for this treatment ratio (because we are testing this ratio against some value, across the range of X). The X values that Y is measured at differs between the treatments, but the ranges are similar. # Reproducible problem: X1<- runif(20, 0.5, 4) X2<- runif(20, 0.5, 4) Y1<- 20*exp(-0.5*X1) + rnorm(20) Y2<- 30*exp(-0.5*X2) + rnorm(20) # Look at data: plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5)) points(X2, Y2, pch=19, col="red") # Full dataset dfr<- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A",20),rep("B",20))) # Fit gam # I use a gamma family here although it is not necessary: in the real problem it is, though. gfit<- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log)) # I am interested in the relationship: # Y(treatment =="B") / Y(treatment=="A") as a function of X, with a confidence interval! Do I just do a bootstrap here? Or is there a more appropriate method? Thanks a lot for any help. greetings, Remko - Remko Duursma Research Lecturer Hawkesbury Institute for the Environment University of Western Sydney Hawkesbury Campus, Richmond Mobile: +61 (0)422 096908 www.remkoduursma.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] gam confidence interval (package mgcv)
But that just gives me the prediction of Y for treatment A or B, not the ratio. As I stated: # I am interested in the relationship: # Y(treatment =="B") / Y(treatment=="A") as a function of X, with a confidence interval! I can get the SE for either of them using predict.gam without a problem, but I don't know how to get the CI for the ratio! thanks, remko - Remko Duursma Research Lecturer Hawkesbury Institute for the Environment University of Western Sydney Hawkesbury Campus, Richmond Mobile: +61 (0)422 096908 www.remkoduursma.com On Tue, Jun 28, 2011 at 2:09 PM, David Winsemius wrote: > > On Jun 27, 2011, at 10:45 PM, Remko Duursma wrote: > >> Dear R-helpers, >> >> I am trying to construct a confidence interval on a prediction of a >> gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant >> but I am not able to apply that to this, different, problem. >> >> Any help is appreciated! >> >> Basically I have a function Y = f(X) for two different treatments A >> and B. I am interested in the treatment ratios : Y(treatment = B) / >> Y(treatment = A) as a function of X, including a confidence interval >> for this treatment ratio (because we are testing this ratio against >> some value, across the range of X). >> >> The X values that Y is measured at differs between the treatments, but >> the ranges are similar. >> >> >> # Reproducible problem: >> X1 <- runif(20, 0.5, 4) >> X2 <- runif(20, 0.5, 4) >> >> Y1 <- 20*exp(-0.5*X1) + rnorm(20) >> Y2 <- 30*exp(-0.5*X2) + rnorm(20) >> >> # Look at data: >> plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5)) >> points(X2, Y2, pch=19, col="red") >> >> # Full dataset >> dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2), >> treatment=c(rep("A",20),rep("B",20))) >> >> # Fit gam >> # I use a gamma family here although it is not necessary: in the real >> problem it is, though. >> gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log)) >> >> # I am interested in the relationship: >> # Y(treatment =="B") / Y(treatment=="A") as a function of X, > > Can't you use predict.gam? > > plot(predict(gfit, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2), > treatment=c(rep("A",37),rep("B",37) ) ) )[1:37] ) > lines(predict(gfit3, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2), > treatment=c(rep("A",37),rep("B",37) ) ) )[-(1:37)]) > >> with a confidence interval! > > There is an se.fit argument to predict.gam(). > > > -- > > David Winsemius, MD > West Hartford, CT > > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] gam confidence interval (package mgcv)
On Jun 27, 2011, at 10:45 PM, Remko Duursma wrote: Dear R-helpers, I am trying to construct a confidence interval on a prediction of a gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant but I am not able to apply that to this, different, problem. Any help is appreciated! Basically I have a function Y = f(X) for two different treatments A and B. I am interested in the treatment ratios : Y(treatment = B) / Y(treatment = A) as a function of X, including a confidence interval for this treatment ratio (because we are testing this ratio against some value, across the range of X). The X values that Y is measured at differs between the treatments, but the ranges are similar. # Reproducible problem: X1 <- runif(20, 0.5, 4) X2 <- runif(20, 0.5, 4) Y1 <- 20*exp(-0.5*X1) + rnorm(20) Y2 <- 30*exp(-0.5*X2) + rnorm(20) # Look at data: plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5)) points(X2, Y2, pch=19, col="red") # Full dataset dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A", 20),rep("B",20))) # Fit gam # I use a gamma family here although it is not necessary: in the real problem it is, though. gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log)) # I am interested in the relationship: # Y(treatment =="B") / Y(treatment=="A") as a function of X, Can't you use predict.gam? plot(predict(gfit, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2), treatment=c(rep("A",37),rep("B",37) ) ) )[1:37] ) lines(predict(gfit3, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2), treatment=c(rep("A",37),rep("B",37) ) ) )[-(1:37)]) with a confidence interval! There is an se.fit argument to predict.gam(). -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] gam confidence interval (package mgcv)
Dear R-helpers, I am trying to construct a confidence interval on a prediction of a gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant but I am not able to apply that to this, different, problem. Any help is appreciated! Basically I have a function Y = f(X) for two different treatments A and B. I am interested in the treatment ratios : Y(treatment = B) / Y(treatment = A) as a function of X, including a confidence interval for this treatment ratio (because we are testing this ratio against some value, across the range of X). The X values that Y is measured at differs between the treatments, but the ranges are similar. # Reproducible problem: X1 <- runif(20, 0.5, 4) X2 <- runif(20, 0.5, 4) Y1 <- 20*exp(-0.5*X1) + rnorm(20) Y2 <- 30*exp(-0.5*X2) + rnorm(20) # Look at data: plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5)) points(X2, Y2, pch=19, col="red") # Full dataset dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A",20),rep("B",20))) # Fit gam # I use a gamma family here although it is not necessary: in the real problem it is, though. gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log)) # I am interested in the relationship: # Y(treatment =="B") / Y(treatment=="A") as a function of X, with a confidence interval! Do I just do a bootstrap here? Or is there a more appropriate method? Thanks a lot for any help. greetings, Remko - Remko Duursma Research Lecturer Hawkesbury Institute for the Environment University of Western Sydney Hawkesbury Campus, Richmond Mobile: +61 (0)422 096908 www.remkoduursma.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.