Re: [R] Prediction interval with GAM?

From: Simon Wood <s.wood_at_bath.ac.uk>
Date: Tue, 19 Apr 2011 16:39:18 +0100

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

## 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

______________________________________________
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.
Received on Tue 19 Apr 2011 - 15:42:15 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 19 Apr 2011 - 16:30:31 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive