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: highs...@highstat.com
URL: www.highstat.com
URL: www.brodgar.com

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Reply via email to