[R] Prediction interval with GAM?

[Simon Wood]

 > Is it possible to estimate prediction interval using GAM?  I looked > 
through ?gam,

- The easiest and most general way is by posterior simulation. Here's an 
example...

## Prediction interval example for Gamma GAM

library(mgcv)

## simulate some data...
f <- function(x) (0.2 * x^11 * (10 * (1 - x))^6 + 10 *
             (10 * x)^3 * (1 - x)^10)/2
x <- runif(200)
fx <- f(x)
Ey <- exp(fx);scale <- .5 ## mean and GLM scale parameter
## Note that `shape' and `scale' in `rgamma' are almost
## opposite terminology to that used with GLM/GAM...
set.seed(8)
y <- rgamma(Ey*0,shape=1/scale,scale=Ey*scale)

## fit smooth model to x, y data...
b <- gam(y~s(x,k=20),family=Gamma(link=log),method="REML")

## extract parameter estiamtes and cov matrix...
beta <- coef(b);Vb <- vcov(b)

## simulate replicate beta vectors from posterior...
Cv <- chol(Vb)
n.rep=10000;nb <- length(beta)
br <- t(Cv) %*% matrix(rnorm(n.rep*nb),nb,n.rep) + beta

## turn these into replicate linear predictors...
xp <- 0:200/200
Xp <- predict(b,newdata=data.frame(x=xp),type="lpmatrix")
lp <- Xp%*%br
fv <- exp(lp) ## ... finally, replicate expected value vectors

## now simulate from Gamma deviates with mean as in fv
## and estimated scale...

yr <- matrix(rgamma(fv*0,shape=1/b$scale,scale=fv*scale),nrow(fv),ncol(fv))

plot(rep(xp,n.rep),yr,pch=".") ## plotting replicates
points(x,y,pch=19,cex=.5) ## and original data

## compute 95% prediction interval...
PI <- apply(yr,1,quantile,prob=c(.025,0.975))
## and plot it...
lines(xp,PI[1,],col=2,lwd=2);lines(xp,PI[2,],col=2,lwd=2)

## Confidence interval for comparison...
pred <- predict(b,newdata=data.frame(x=xp),se=TRUE)
lines(xp,exp(pred$fit),col=3,lwd=2)
u.ci <- exp(pred$fit + 2*pred$se.fit)
l.ci <- exp(pred$fit - 2*pred$se.fit)
lines(xp,u.ci,col=3,lwd=2);lines(xp,l.ci,col=3,lwd=2)





On 19/04/11 10:04, yosuke kimura wrote:
> Hello,
>
> Is it possible to estimate prediction interval using GAM?  I looked through
> ?gam, ?predict.gam etc and the mgcv.pdf Simon Wood.  I found it can
> calculate confidence interval but not clear if I can get it to calculate
> prediction interval.  I read "Inference for GAMs is difficult and somewhat
> contentious." in Kuhnert and Venable An Introduction to R, and wondering why
> and if that is the reason that current version of mgcv does not provide
> prediction interval.
>
> I am thinking about GAM because ggplot2 uses it to fit the data and it looks
> working well.  There are many other models and wondering what is can be a
> good model to gives prediction interval.  I am trying to estimate perimeter
> of buildings from projected area of buildings.  I have both measurement for
> some of data but I know only area for the rest.
>
> Thank you for your help,
>
> Yosuke
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283