[R] gamm() and predict()

[Wakefield, Ewan]

Dear All,
I have a query relating to the use of the ‘predict’ and ‘gamm’ functions.  I am dealing with large (approx. 5000) sets of presence/absence data, which I am trying to model as a function of different of environmental covariates.  Ideally my models should include individual and colony as random factors. I have been trying to fit binomial models using the gamm function to achieve this. For the sake of simplicity I have adapted a some of the example code from ?gamm to illustrate some problems I have been having predicting values using this approach.
### Begin example ###
library(mgcv)
## Generate some example data
set.seed(0)
n <- 400
sig <- 2
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sig)
y <- f + e
## Change the response to binary
y <-round(y/max(y))
## Add a factor to the linear predictor, to be modelled as random
fac <- rep(1:4,n/4)
f <- f + fac*3
fac<-as.factor(fac)
## Fit an additive model
mod<-gamm(y~s(x0)+s(x1)+s(x2)+s(x3),family=binomial,random=list(fac=~1))
## Generate some new example data
new.dat = data.frame(x0 = runif(n, 0, 1), x1 = runif(n, 0, 1), x2 = runif(n, 0, 1),x3 = runif(n, 0, 1), fac = fac)
## Predict response values using the original data and the gam part of the model…
predict(mod$gam, type = "response")
## Predict response values using the new data and the gam part of the model…
predict(mod$gam, type = "response", new.dat)
## Predict response values using the original data and the glmm part of the model…
predict(mod$lme, level = 0, type = "response")
## Predict response values using the new data and the gam part of the model…
predict(mod$lme, level = 0, type = "response", new.dat)
# This produces the error message ‘Error in eval(expr, envir, enclos) : object "fixed" not found’
## End example ###
My questions are as follows:
1. I presume predict(mod.$gam) produces population level predictions. Is this correct?
2. Is it possible to extract standard errors using predict(mod.$gam) or is there a more suitable approach to estimating confidence in prections made with gamms?
3. It seems that predict(mod.$lme) results in predictions at the level of random factors. Furthermore, these appear to be on the scale of the linear predictor regardless of how level is specified (see ?glmmPQL). Is this correct?
4. The code predict(mod$lme, new.dat) produces an error message, seemingly indicating the the fixed effects are missing from my new data frame (see example). Am I doing something wrong here?
5. Is it possible to predict both population and random factor level predictions using new data with gamm objects?
I have read all the relevant help files including those associated with glmmPQL and also Simon Woods book and I am still a bit confused so any help would be gratefully received. I am using R 2.6.1 with Windows XP pro, mgcv version 1.3-29.
Thanks,
Ewan Wakefield
British Antarctic Survey
High Cross
Madingley Road
Cambridge
UK