[R] Predicting complicated GAMMs on response scale

[Gavin Simpson]

On Mon, 2009-05-18 at 11:48 -0700, William Paterson wrote:
> Hi,
> 
> I am using GAMMs to show a relationship of temperature differential over
> time with a model that looks like this:-
> 
> gamm(Diff~s(DaysPT)+AirToC,method="REML") 
> 
> where DaysPT is time in days since injury and Diff is repeat measures of
> temperature differentials with regards to injury sites compared to
> non-injured sites in individuals over the course of 0-24 days. I use the
> following code to plot this model on the response scale with 95% CIs which
> works fine:-
> 
> g.m<-gamm(Diff~s(DaysPT)+AirToC,method="REML")
> p.d<-data.frame(DaysPT=seq(min(DaysPT),max(DaysPT)))
> p.d$AirToC<-(6.7)
> b<-predict.gam(g.m$gam,p.d,se=TRUE)
> range<-c(min(b$fit-2*b$se.fit),max(b$fit+2*b$se.fit))
> plot(p.d$DaysPT,b$fit,ylim=c(-4,12),xlab="Days post-tagging",ylab="dTmax
> (ºC)",type="l",lab=c(24,4,12),las=1,cex.lab=1.5, cex.axis=1,lwd=2)
> lines(p.d$DaysPT,b$fit+b$se.fit*1.96,lty=2,lwd=1.5)
> lines(p.d$DaysPT,b$fit-b$se.fit*1.96,lty=2,lwd=1.5)
> points(DaysPT,Diff)
> 
> 
> However, when I add a correlation structure and/or a variance structure so
> that the model may look like:- 
> 
> 
> gamm(Diff~s(DaysPT3)+AirToC,correlation=corCAR1(form=~DaysPT|
> Animal),weights=varPower(form=~DaysPT),method="REML")
> 
> 
> I get this message at the point of inputting the line
> "b<-predict.gam(g.m$gam,p.d,se=TRUE)"

Note that p.d doesn't contain Animal. Not sure this is the problem, but
I would have thought you'd need to supply new values of Animal for the
data you wish to predict for in order to get the CAR(1) errors correct.
Is it possible that the model is finding another Animal variable in the
global environment?

I have predicted from several thousand GAMMs containing correlation
structures similar to the way you do above so this does work in general.
If the above change to p.d doesn't work, you'll probably need to speak
to Simon Wood to take this further.

Is mgcv up-to-date? I am using 1.5-5 that was released in the last week
or so.

For example, this dummy example runs without error for me and is similar
to your model

y1 <- arima.sim(list(order = c(1,0,0), ar = 0.5), n = 200, sd = 1)
y2 <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 200, sd = 3)
x1 <- rnorm(200)
x2 <- rnorm(200)
ind <- rep(1:2, each = 200)
d <- data.frame(Y = c(y1,y2), X = c(x1,x2), ind = ind, time = rep(1:200,
times = 2))
require(mgcv)
mod <- gamm(Y ~ s(X), data = d, corr = corCAR1(form = ~ time | ind),
            weights = varPower(form = ~ time))
p.d <- data.frame(X = rep(seq(min(d$X), max(d$X), len = 20), 2),
                  ind = rep(1:2, each = 20),
                  time = rep(1:20, times = 2))
pred <- predict(mod$gam, newdata = p.d, se = TRUE)

Does this work for you? If not, the above would be a reproducible
example (as asked for in the posting guide) and might help Simon track
down the problem if you are running an up-to-date mgcv.

HTH

G

> 
> 
> Error in model.frame(formula, rownames, variables, varnames, extras,
> extranames,  : 
>         variable lengths differ (found for 'DaysPT')
> In addition: Warning messages:
> 1: not all required variables have been supplied in  newdata!
>  in: predict.gam(g.m$gam, p.d, se = TRUE) 
> 2: 'newdata' had 25 rows but variable(s) found have 248 rows 
> 
> 
> Is it possible to predict a more complicated model like this on the response
> scale? How can I incorporate a correlation structure and variance structure
> in a dataframe when using the predict function for GAMMs?
> 
> Any help would be greatly appreciated.
> 
> William Paterson
> 
> 
> 
> 
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%