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 [EMAIL PROTECTED]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
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
