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

[Gavin Simpson]

On Thu, 2013-07-18 at 14:34 +0000, Henkelman, Jonathan wrote:
> Gavin,
> 
> Welcome to the prairies...

Cheers

<snip />

> 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.

Exactly; I've had the same problem modelling smooth trends with residual
correlations.

> 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?  

Sorry I missed one important aspect of approach 2) (from my earlier
email). What I suggested was fit the random effects model without the
autocorrelation structure, then add the appropriate correlation
structure in in a second run after looking at the ACF/PACF to determine
what order ARMA to use. What I should have said is that you will have to
set the parameters of the ARMA model explicitly in `corARMA()` (see
`?corARMA` for details on how to do that.

To actually do this then you'll have to estimate an ARMA for the
normalised residuals from the random effects fit. You'll have to do this
*within* the Plots. That will give many models (`Plot` models) and
you'll need to choose common parameter values from the set of fits -
`corARMA()` fits the same structure within each level of Plots. So the
resulting values you plug into `corARMA()` will be say an average of the
coefficients derived from the individual Plot-level ARMA fits. Note this
won't be perfect, but as Zuur et al note, the idea is not to model the
covariance structure perfectly but to get something that is reasonable.

Beyond that, you're into Bayesian modelling to have fine control, or do
as Brian Cade suggests and put everything in the fixed effects if you
have enough DFs.

Fitting these models (LMEs with correlation structures) can get very
hard; it is nice when they work without having to impart much effort,
but very often the optimisation heads off down the wrong path and the
only way to correct that is to estimate somethings outside lme() and
then plug that back into lme() as known.

> I think question 2) is answered given the citation above -- it is
> valid conceptually.  Perhaps nlme IS the wrong package, or there is a
> bug in nlme since the information _seems_ to be there -- can anyone
> comment on this.

This isn't a bug, it is just a feature (infelicity) of the optimisation
algorithm that is choosing parameters. It is trying to choose values for
the variance of the random effects and the parameters for the ARMA part.
It obviously got a better fit by saying that the residual variance was
better fitted by the AR (or ARMA) than by putting some in the random
effects. Hence by estimating the ARMA coefficients after fitting the
random effect model, and then using these as fixed values in `corARMA()`
might allow the model to be fitted.

HTH

G

>   But I do still wonder about suppostion 3: Is it valid to use a
> two-step approach: run a basic autocorrelation model to remove the
> autocorrelation effect (normalized residuals shown in Figs 1,2), then
> run a mixed intercept model on the residuals to deal with the nested
> structure (and also the heterogeneity between plots).  This runs,
> gives reasonable results that agree with the graphs, but is it
> statistically valid?
> 
> Thanks,
> Jonathan
> ________________________________________
> From: Gavin Simpson [gavin.simpson at ucl.ac.uk]
> Sent: July 17, 2013 10:06 PM
> To: Henkelman, Jonathan
> Cc: r-sig-ecology at r-project.org; Johnstone, Jill
> Subject: Re: [R-sig-eco] Mixed effect (intercept) model with AR1 autocorrelation structure
> 
> On Wed, 2013-07-17 at 16:10 +0000, Henkelman, Jonathan wrote:
> > Perhaps I should clarify.  There is a time-series trend -- daily
> > temperature fluctuates randomly throughout the summer.  But there is
> > not a clear long-term signal.  I have modelled the time-series effect
> > using a gam to see if that can adequately compensate for the effect.
> > However, I believe this is a fundamentally flawed approach:
> >
> > 1) We are not interested in modelling the time-series; it merely is a
> > way of estimating the temperature response in plot.  That is, I don't
> > want to ask the question, can I predict the temperature given the
> > treatment and the date, but rather, is there a treatment effect.  My
> > inference ability is badly reduced if a model a gam.
> 
> It is *irrelevant* whether you want to model the time series. You *must*
> model it either via the fixed effects or in the covariance matrix of the
> residuals.
> 
> You have two options:
> 
> 1) include a time effect (either linear or spline) plus possibly a
> simple time series model for the residuals (AR(1) would be a start).
> Your null would include the time effect thus you can assess the effect
> of Treatment relative to this null.
> 
> 2) proceed as you did but use a much more complex ARMA model for the
> residuals. AR(1) is clearly inappropriate with such a high \phi
> estimate.
> 
> Which is most useful will depend on the data; it sounds as if you don't
> have a deterministic trend as you seem to be describing a stochastic
> trend. However, I'd be very surprised if you didn't have a cyclic
> temperature effect with temperatures varying smoothly through the year
> in a seasonal signal/cycle, with noise superimposed. But I am just
> guessing.
> 
> Perhaps start with the mixed effects and **no** correlation structure if
> you don't want to do 1). Extract the normalised residuals and look at
> the ACF and the PACF of these. That will help identify what order ARMA
> model might be needed. Then add the correlation structure with
> corARMA().
> 
> I didn't notice you were at USask; I move to the UofR a few months ago.
> If it would help to discuss offlist, my UofR contact details are in the
> footer.
> 
> HTH
> 
> G
> 
> > 2) The current analysis is for a single season.  In a few years we
> > will be re-running this analysis of 5 years of data.  I do not expect
> > the random fluctuation in seasonal temperature will be the same each
> > year.  Hence, while this analysis sort of works now, it won't in the
> > future.  However, it seems reasonable to model the autocorrelation
> > effect within the time-series as constant through time.
> >
> > 3) When I look at the process of temperature I can say, yes today is
> > more likely to be similar to yesterday than the day before.  There is
> > autocorrelation and random fluctuation, hence it makes sense to model
> > it this way.  For the record, as simple AR1 model better account for
> > the seasonal fluctuation than a gam, and my ARMA(2,0) model does an
> > even better job.
> >
> > Hope this helps, J
> >
> 
> --
> 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
> 
> 
> 

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