[R] prediction, gam, mgcv
Simon Wood
simon at stats.gla.ac.uk
Mon Feb 28 15:45:40 CET 2005
> 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
More information about the R-help
mailing list