[R-sig-ME] Generalised Linear Mixed Model with Moving Average

[Ben Bolker]

Mark Bilton <mcbilton at ...> writes:

> I am having a few problems with a glmer I am trying to create in R.

 [snip]
 
> My data to be modelled is the abundance of a plant species in a fixed quadrat
> (and for the purposes of this post, just 1 species).
> I have a spatial hierachical design. At a site, I have 15 'plots', each plot
> contains 5 'quadrats'.

  A single site, or multiple sites?
 
> I have a 'treatment' of 3 climates (control, +water, -water), and
> so for the 15
> plots 5 of them were assigned to each treatment (therefore 5 is my 'real'
> replication per treatment).
> 
> The data for each species was also recorded for 8 'years'.
> The first year no treatment was applied to any plots.
> 
> My objective is to find if there is a change in species abundance across time
> (slope) for the different climate treatments.
> 
> Firstly, I will only talk about the full model 
> (and not the comparisons for full
> significance testing).
> 
> Year is a continuous variable.
> Plot, Treatment, and Quadrat are factors.

  Do you really need to keep quadrats separate?  If you take the
mean of quadrats to get a single value for each plot, are your
data adequately/approximately normally distributed?  This would let
you do a LMM instead of a GLMM, which would simplify your life.
In fact, then you might not need a LMM at all.

> 
> The model I have then is:
> 
> Model1<-glmer(Species~Year*Treatment+(1|Plot), REML = FALSE, family=poisson,
> na.action = na.omit)

  REML=FALSE is redundant/ignored in glmer, as explained in a recent
thread on this list.

> QUESTION 1:
> I do not include 'Quadrat' as a random variable 
> (as stated in Crawley R book for
> the lowest level of the hierarchy). However, I am not sure if this is correct,
> as I do get slightly different results when included or not. 
> So, for the spatial
> nesting should I use ((1|Plot/Quadrat) ??

  If you include quadrat, then you are essentially allowing for
extra-Poisson variation in the samples, or equivalently using
a lognormal-Poisson rather than a Poisson distribution.  If your
model is feasible and there is a non-negligible amount of variability
estimated at the among-quadrat/within-plot level, you probably
do want to retain quadrat as a random variable (or lump and hope
for normality, see above).

> QUESTION 2:
> Do I need to include 'Year' in the random variables as well ?
> I have tried  (1|Plot) + (Year|Plot) (doesn't converge) and just (Year|Plot)
> (doesn't converge) and (as.factor(Year)|Plot) (doesn't converge). The lack of
> converging may be due to a large amount of zero's in my dataset, 
> but not sure if
> this is the case or not.

  If you include (Year|Plot), you are allowing for
the possibility that there is random variation in slopes (as well
as intercepts) across plots.  as.factor(Year)|Plot allows for random
(not necessarily linear) year-to-year variation that varies across plots.

(1|Plot) and (Year|Plot) should be the same, because the latter includes
an implicit intercept term.
> 
> QUESTION 3:
> I could include a covariate for the level of the species in the first year.
> Although, I am not sure this is necessary as I am really interested 
> in slope (ie
> the interaction between Year and Treatment), and a covariate should not alter
> this. Do you agree with my statement, or do you think I should include it ?
> If I do include it, is this a fixed or random variable. 
> My feeling was random, but a
> colleague disagreed.

  It's certainly not impossible that there could be a relationship
between overall level and slope. Have you tried looking at the data?
Fixed makes more sense to me.
 
> QUESTION 4:
> I also want to test some moving average models using this framework. Is it
> possible to do this with glmer ?? I tried coding a simple function myself for
> different lag distances, but the averages then are not integers, which the
> poisson link function does not like.

  glmer does not allow so-called 'R-side' correlation structures;
nor does lmer.  You can do this in lme, and possibly handle a GLMM
with such correlation structures in glmmPQL (in the MASS package),
but glmmPQL has some disadvantages.

  I would suggest you look at the residuals from an otherwise
full model to see if there is any hint of temporal autocorrelation ...
this is a pretty small/short data set for detecting temporal autocorrelation.
> 

  I'm going to let someone else have question 5.