[R] gam confidence interval (package mgcv)

[David Winsemius]

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