> 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