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
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
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.83183 0.03885 47.152 < 2e-16 ***
treatmentB 0.44513 0.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 at 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.
>