[R] Comparing GAMMs

[Gavin Simpson]

On Wed, 2009-07-08 at 11:24 -0400, Paul Simonin wrote:
> Greetings!
>   I am looking for advice regarding the best way to compare GAMMs. I 
> know other model outputs return enough information for R's AIC, ANOVA, 
> etc. commands to function, but this is not the case with GAMM unless one 
> specifies the gam or lme portion. I know these parts of the gamm contain 
> items that will facilitate comparisons between gamms. Is it correct to 
> simply use these values for this purpose? For example, the lme portion 
> of the gamm returns a log liklihood value that could be used to 
> calculate information criteria. However, I am wondering whether entire 
> gamms be compared using this, or only the lme part.
>   Maybe my thinking about the lme and gam portions of gamms is 
> incorrect? If this appears to be the case, let me know! In general, if 
> someone could clarify my understanding in any way it would be much 
> appreciated.
> Thank you very much!
> Sincerely,
> Paul Simonin

Hi Paul,

Are your GAMMs Guassian (i.e. AMM) or non-Gaussian? If they are
Gaussian, then

anova(mod1$lme, mod2$lme)

gives an approximate LRT for the two models. That will also yield AIC
and BIC which might also be used for inference. Your AMM in this case is
just a linear mixed model and these usual forms of inference apply, with
the caveat that the hypothesis testing is approximate. You end up using
both the $lme and the $gam components for various aspects of model
inspection, interrogation etc, but for hypothesis testing, the lme bit
is sufficient. You can also use things like intervals(mod1$lme) to look
at confidence on the smoothing parameters. See Simon Wood's book [1]
section 6.7 for more details, and preceding sections on how the
smoothers can be formulated as a mixed model.

If your GAMMS are generalised then I'm not sure what the best approach
for comparison or hypothesis testing might be - especially as this is an
ongoing research topic for GLMMs, and also because of the method by
which GAMMs are fitted in mgcv. Simon Wood says as much in his 2006
monograph [1, page 318, section 6.6.2]. The non-Gaussian case uses
glmmPQL from package MASS, and this doesn't return a likelihood and
hence no AIC (in the same way that quasi families in glm() fits don't
return likelihoods).

So having said that, if you do have a likelihood, then you must be
fitting AMM via gamm() and the first half of my reply would seem most
appropriate.

[1] Wood, S.N. (2006) Generalized Additive Models; an Introduction with
R. Chapman & Hall/CRC.

HTH

G

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%