R help - Feb 2005 - prediction, gam, mgcv

I fitted a GAM model with Poisson distribution
using the function gam() in the mgcv package.

My model is of the form:
mod<-gam(y~s(x0)+s(x1)+s(x2),family=poisson).


To extract estimates at a specified set of covariate
values I used the gam `predict' method. 

But I want to get
estimate and standard error of the difference of two fitted values.

Can someone explain what should I do?

Thank you
gabriele

> My model is of the form:
> mod<-gam(y~s(x0)+s(x1)+s(x2),family=poisson).
> But I want to get
> estimate and standard error of the difference of two fitted values.
> 
The following code shows you how to do this (a) for differences on the 
scale of the linear predictor, and (b) for differences on the scale of the 
response. (Could be made more efficient at the cost of being a bit less 
readable)

best,
Simon


## simulating data
n <- 300
x <- runif(n)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
f <- exp(f/5)
y <- rpois(f,f)

## fit gam
mod <- gam(y~s(x),family=poisson)

## creat prediction data frame
pd <- data.frame(x=c(.5,.6))

## create matrix, Xp, mapping parameters to linear predictor
Xp <- predict(mod,pd,type="lpmatrix")

### First working on scale of linear predictor, repeatedly using
### fact that if X=CY where C is a matrix of coefficients and
### X and Y are random vectors, then cov(X) = C cov(Y) C'

## extract model coefficients and their covariance matrix
b <- coef(mod)
Vb <- mod$Vp

## obtain predictions and their covariance matrix
pred <- Xp%*%b
Vpred <- Xp%*%Vb%*%t(Xp)

## find difference and associated standard deviation
diff <- c(1,-1)
diff.pred <- t(diff)%*%pred
diff.sd <- sqrt(t(diff)%*%Vpred%*%diff)

diff.pred;diff.sd

### Now working on response scale, by simulation

library(MASS)
## simulate 1000 rep. param. vectors from posterior distn...
br <- mvrnorm(n=1000,b,Vb)
diff.pred <- rep(0,1000) 
for (i in 1:1000) {     ## for each rep
  pred <- Xp%*%br[i,]   ## get l.p. predictions
  ## and hence calculate the required response scale difference
  diff.pred[i] <- exp(pred[1])-exp(pred[2])
}
## diff.pred now contains 1000 rep. differences, so can easily estimated 
## their standard deviation....
sd(diff.pred) 
## and here is the point estimate ... 
pred <- Xp%*%b
exp(pred[1])-exp(pred[2])


_____________________________________________________________________
> Simon Wood simon at stats.gla.ac.uk        www.stats.gla.ac.uk/~simon/
>>  Department of Statistics, University of Glasgow, Glasgow, G12 8QQ
>>>   Direct telephone: (0)141 330 4530          Fax: (0)141 330 4814