[R] autocorrelation in gams

[Gavin Simpson]

On Thu, 2008-08-14 at 16:12 +0100, Abigail McQuatters-Gollop wrote:
> Hi,
> 
> I am looking at the effects of two explanatory variables on chlorophyll.
> The data are an annual time-series (so are autocorrelated) and the
> relationships are non-linear. I want to account for autocorrelation in
> my model. 
> 
>  
> 
> The model I am trying to use is this:
> 
>  
> 
> Library(mgcv)
> 
>  
> 
> gam1 <-gam(Chl~s(wintersecchi)+s(SST),family=gaussian,
> na.action=na.omit, correlation=corAR1(form =~ Year)) 
> 

There is no correlation argument in mgcv::gam you need gamm(). gam() has
a ... argument which I suspect is silently mopping up your correlation
argument so that no error/warning is raised.

Note that gamm() uses lme from the nlme package (in the Gaussian case)
and works that function very hard (see Wood 2006 GAM book). In my
experience with this function separating trend from the correlation is
quite difficult when also estimating the degree of smoothness and one
has to work hard with the options.

As an alternative you might also take a look at the paper by Ferguson et
al (2007):

http://www3.interscience.wiley.com/journal/119392062/abstract?CRETRY=1&SRETRY=0

Which has R code using the sm package to do a very similar thing.

HTH

G
 
> 
> the result I get is this: 
> 
>  
> 
> Family: gaussian 
> 
> Link function: identity 
> 
>  
> 
> Formula:
> 
> CPRChl ~ s(wintersecchi) + s(SST)
> 
>  
> 
> Parametric coefficients:
> 
>             Estimate Std. Error t value Pr(>|t|)    
> 
> (Intercept)  3.57000    0.05061   70.54   <2e-16 ***
> 
> ---
> 
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
>  
> 
> Approximate significance of smooth terms:
> 
>                   edf Est.rank     F p-value   
> 
> s(wintersecchi) 2.445        5 4.672 0.00887 **
> 
> s(SST)          2.408        5 4.301 0.01237 * 
> 
> ---
> 
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
>  
> 
> R-sq.(adj) =  0.676   Deviance explained = 75.4%
> 
> GCV score = 0.074563   Scale est. = 0.053781  n = 21
> 
>  
> 
> The result look good - significant, with a lot of deviance explained,
> but I am not convinced the model is actually accounting for the
> autocorrelation (based on the formula in the results). How can I tell? 
> 
>  
> 
> Many thanks,
> 
>  
> 
> 
> 
> 
> 
> Dr Abigail McQuatters-Gollop
> 
> Sir Alister Hardy Foundation for Ocean Science (SAHFOS)
> 
> The Laboratory
> 
> Citadel Hill
> 
> Plymouth UK PL1 2PB
> 
> tel: +44 1752 633233
> 
>  
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.