[R-sig-eco] Mixed effect (intercept) model with AR1 autocorrelation, structure

[Highland Statistics Ltd]

> Gavin,
>
> Welcome to the prairies...
>
> I think this thread is getting off track... a few comments on the autocorrelation model, then bringing it back to the question I was trying to ask:
>
> I agree completely that we *must* model the time-series because our readings are not independent and hence we are violating the assumptions of a standard model.  There is no clear deterministic trend to the temperature readings: we are only modelling summer temperature and the single arch of the season is lost in weekly stochasticity.  There is also no significant linear trend (p-value = 0.3).  As I state in my first paragraph I have modelled the summer trend with a gam as an alternative approach (highly significant with 9 df).  However, even the basic AR1 model gives cleaner residuals.  Also, as I mention in bullet 3, I have an ARMA(2,0) model working that does a better job of dealing with the autocorrelation (as shown by ACF); I will expand this model once the conceptual problems I presented are sorted out.
>
> The problems I have with modelling the summer temperatures as a deterministic effect (fixed effect, gam) are outlined in my previous post: it limits our inference (as per Zuur 327, although he is discussing a slightly different context); and whatever gam I come up with this year won't be consistent with next year, whereas an autocorrelation model will be.  I also think it is more ecologically valid. I can see there could be different points of view here however and debate is healthy.
>
> Coming back to what I see as the problem at hand...  I _have_ several models that adequately deal with the autocorrelation.  Normalized residuals for the AR1 model plotted through time are shown in Figure 2 (see end of paragraph)  However, residuals plotted against Plot# are shown in Figure 1.  I don't have (too much) of a problem with Figure 2 (although variance seems to be decreasing with time)... I do have a big problem with figure 1!  This is what I wanted to address in my new model -- a random intercept model to account for the different means (and variances) of the residuals in each plot.
>
> [My plots don't include in the posting so I've put them on the web:
> http://homepage.usask.ca/~jmh842/Figure%201.pdf
> http://homepage.usask.ca/~jmh842/Figure%202.pdf
> ]
>
> However, this model doesn't work! -- it runs, but gives no random effect (see original post for code and results...).  The residuals of this new model show an identical pattern to Figure 1.  Hence something is going wrong in the computation element within nlme::lme.  This is the question I wanted to address with this post.
>
> I was up examining Zuur last night and found a passing reference to this problem on page 160 "It may be an option to extend the model with the AR1 correlation, with a random intercept on nest [a case which would correspond exactly to my scenario]. Such a model allows for the compound correlation between all observations from the same nest, and temporal correlation between observations from the same nest _and_ night.  *But the danger is that the random intercept and auto-correlation will fight with each other for the same information.* " [* are mine...] I think this is exactly what is happening: the two correlation specifications are fighting with each other and the mixed intercept is loosing.
>
> I have had many suggestions about other ways to model this dataset.  I have been and will continue to pursue these ideas (I'd also like to build a simple Baysian model as Zuur does in chapter 23 to deal with a similar problem).  However, I have not had anyone comment directly on my two main questions:
> 1) is this analysis doable or is nlme perhaps the wrong package?; and
> 2) Is this a reasonable question to ask?

If you want to go this route then I suggest that you follow Chapter 5 in 
'Beginner's Guide to GLM & GLMM' by Zuur, Hilbe, Ieno. Is shows you how 
to include an AR 1 correlation within a Poisson GLM. And Chapter 9 in 
Zuur et al (2012) shows how to program your own smoother. You can easily 
do all this stuff in JAGS.

As to smoothers and auto-correlation..I suggest you fix either the 
amount of smoothing of the smoother(s), or the auto-correlation 
parameter to a fixed value. Or perhaps both. This is the key to getting 
it to run!

Alain






-- 

Dr. Alain F. Zuur
First author of:

1. Analysing Ecological Data (2007).
Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p.
URL: www.springer.com/0-387-45967-7


2. Mixed effects models and extensions in ecology with R. (2009).
Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer.
http://www.springer.com/life+sci/ecology/book/978-0-387-87457-9


3. A Beginner's Guide to R (2009).
Zuur, AF, Ieno, EN, Meesters, EHWG. Springer
http://www.springer.com/statistics/computational/book/978-0-387-93836-3


4. Zero Inflated Models and Generalized Linear Mixed Models with R. (2012) Zuur, Saveliev, Ieno.
http://www.highstat.com/book4.htm

Other books: http://www.highstat.com/books.htm


Statistical consultancy, courses, data analysis and software
Highland Statistics Ltd.
6 Laverock road
UK - AB41 6FN Newburgh
Tel: 0044 1358 788177
Email: highstat at highstat.com
URL: www.highstat.com
URL: www.brodgar.com