example.R

library(mgcv)
n<-200
sig <- 2
dat <- gamSim(1,n=n,scale=sig)

b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)

newd <- data.frame(x0=(0:90)/90,x1=(0:90)/90,x2=(0:90)/90,x3=(0:90)/90)
pred <- predict.gam(b,newd)

#########################################################
## now get variance of sum of predictions using lpmatrix
#########################################################

Xp <- predict(b,newd,type="lpmatrix")

##################################################################
## The following shows how to use use an "lpmatrix" as a lookup
## table for approximate prediction. The idea is to create
## approximate prediction matrix rows by appropriate linear
## interpolation of an existing prediction matrix. The additivity
## of a GAM makes this possible.
## There is no reason to ever do this in R, but the following
## code provides a useful template for predicting from a fitted
## gam *outside* R: all that is needed is the coefficient vector
## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or
## higher order interpolation for higher accuracy.
###################################################################

xn <- c(.341,.122,.476,.981) ## want prediction at these values
ncols <- length(colnames(Xp))
print( Xp )
x0 <- 1         ## intercept column
dx <- 1/90      ## covariate spacing in `newd'
for (j in 0:2) { ## loop through smooth terms
    cols <- 1+j*9 +1:9      ## relevant cols of Xp
    i <- floor(xn[j+1]*90)  ## find relevant rows of Xp
    print( paste0( 'i', i ) )
    w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
    ## find approx. predict matrix row portion, by interpolation
    x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
}
dim(x0)<-c(1,ncols)
fv <- x0%*%coef(b) + xn[4]    ## evaluate and add offset
## compare to normal prediction
print( predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
                                    x2=xn[3],x3=xn[4]),se=FALSE) )
print( fv )
# print( coef(b) )