[R] Multiple interaction terms in GAMM-model

[Gavin Simpson]

On Thu, 2013-07-25 at 04:56 -0700, Jeroen wrote:
> Dear all,
> 
> I am trying to correlate a variable tau1 to a set of other variables (x1,
> x2, x3, x4), taking into account an interaction with time ('doy') and place
> ('region'), and taking into account dependency of data in time per object
> ID. My dataset looks like:

<snip />

> Since the data is dependent in time per objectID, I use a GAMM model with an
> autocorrelation function. Since each variable x1, x2, etc. is dependent on
> time and place, I should incorporate this as well. 
> 
> Therefore I am wondering if the following gamm-model is correct for my
> situation:
> 
> model <- gamm( tau1 ~ te( x1, by= doy ) + te( x1, by= factor( region ) ) +
> ... + te( x4, by= doy ) + te( x4, by= factor( region ) ) + factor( region ),
> correlation= corAR1(form= ~ doy|objectID ), na.action= na.omit ).
> 
> Does anyone know if this is ok? 

The model looks a little odd - the way you've written it you don't need
`te()` smooths. For the interaction with doy, I would try

te(x1, doy, bs = c("cr", "cc"))

assuming x1 and doy are continuous. This specifies cubic splines and
cyclic cubic splines respectively for x1 and doy. I use a cyclic spline
for doy as we might not expect Jan 1st and Dec 31st to be much
different. By separating the interaction between x1 and doy and the one
for x1 and region you are saying that the way the effect of x1 varies
through the year does not change between regions. Is that reasonable?

I think you need something for the objectID - a random effect or a fixed
effect will account for differences in mean values of taul between
objects. A random effect seems most useful otherwise you'll eat into
your degrees of freedom, but you have plenty of those.

The correlation part assumes that the residuals follow an AR(1) applied
within the objects, with each AR(1) having the same coefficient.
Autocorrelation decreases exponentially with lag/separation. It is a
reasonable start. Fitting will be much quicker without this. Could you
try fitting without and then extract the normalised residuals to check
for residual correlation?

I hope you have a lot of memory and a lot of spare time available; my
experience of fitting similarly sized (though not as complex - I had one
*very* long series) was that gamm() used large amounts of RAM (I had
16GB and most was used for a single model) and took many days to fit.
You may find the random effect for object fights with the AR(1) so you
could end up with an odd fit if it fits at all.

You are going to want to turn on verbosity when fitting so you can see
how the optimisation is proceeding. Something like

require(nlme)
## need a control object
ctrl <- lmeControl(msVerbose = TRUE,
                   maxIter = 400,
                   msMaxIter = 400,
                   niterEM = 0, ## this is VERY important for gamm()!
                   tolerance = 1e-8,
                   msTol = 1e-8,
                   msMaxEval = 400)

then add `control = ctrl` to the `gamm()` call.

I would approach the expecting the model to fail to fit so be prepared
to simplify it etc.

Good luck,

G

> Or should I use a model which also includes terms like " te( x1 ) + ... +
> te( x4 )".
> And is the correlation function correct?
> 
> Thanks so much!!
> 
> Jeroen
> 
> 
> 
> --
> View this message in context: http://r.789695.n4.nabble.com/Multiple-interaction-terms-in-GAMM-model-tp4672297.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Gavin Simpson, PhD                          [t] +1 306 337 8863
Adjunct Professor, Department of Biology    [f] +1 306 337 2410
Institute of Environmental Change & Society [e] gavin.simpson at uregina.ca
523 Research and Innovation Centre          [tw] @ucfagls
University of Regina
Regina, SK S4S 0A2, Canada