Prediction interval with GAM?
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