> Natasha <natasha.lloyd <at> gmail.com> writes:
For future reference, you may have better luck *either* on
the r-sig-mixed-models help list (for advice about GLMMs) *or*
on r-sig-ecology (for advice about ecological studies) [it's not
considered polite to cross-post: try looking at the archives of
both lists to decide which is more appropriate].
> Background information:
> Our study design: We used mark-recapture trapping to obtain density
> estimates and demography of prairie dogs. The location of our study sites
> are nested, we have 6 trapping plots within 3 colonies of prairie dogs.
> Within each colony 3 plots are treatment (given supplemental food) and 3
are
> control sites. Therefore, in total we have 18 study plots (9 control and 9
> treatment).
> Our study spans over two years and within each year we trapped prairie dogs
> twice (Spring 2008, Summer 2008, Spring 2009, Summer 2009) for a total of
> four trapping sessions. Therefore we also have temporal repeated measures
in
> our sampling design. From the mark-recapture data we estimated density
using
> Program DENSITY for each study plot and for each session (a total of 4
> density estimates for each 18 plots over time). We also did vegetation
> sampling to obtain estimates of the natural food available in each plot,
> which gave us total biomass and % edible vegetation as covariates.
>
> We are using the lme4 package in program R (2.11.1) for analysis:
> We fitted a GLMM with Laplace restricted maximum likelihood estimation.
> Since our response variable density is a count with overdispersion we used
a
> corrected poisson distribution with a log link function.
(I'm not sure "corrected" is a technically meaningful term:
do Zuur et al use this?)
> And in order to
> account for the nested and repeated measures in our design we want to
> include Plot nested in Colony as random effects. The explanatory variables
> include Treatment, Session, Total biomass, and/or %Edible food for the
fixed
> effects. (Colony, Plot, Treatment, and Session are all factor variables).
>
> And an example of the R code that we started with is:
>
> PDmodel_1 <- lmer(Density~1+Treatment+Session+Biomass+(1|Colony/Plot),
> family=quasipoisson, data=DensityPD)
(You don't need the "1" in this formula: the intercept is
implicitly
included.)
> My main question is:
> Does this accurately take into account our nested and repeated measures
> design? Have we accounted for any possible pseudoreplication with this type
> of analysis? Do I need to also look at type AR-1 error structures (or can
I
> do that in a glmm?) Or is there a simpler way you might suggest? I can
also
> look at the change in density between time intervals instead density at
time
> intervals if that makes anything easier to analyze? I can also split
> Session into season and year if that is helpful too.
This looks reasonable. You're likely to have trouble estimating
the among-colony variance, because there are only 3 colonies (try
estimating a variance from 3 samples sometime) -- I wouldn't be
surprised if the estimate of the among-colony variance is zero, unless
the colonies are very different. You might alternatively try
Colony as a fixed effect and Colony:Plot as the random effect
(or equivalently make sure your plots have 18 distinct labels rather
than being numbered 1..6 in each plot).
It's a bit tricky to incorporate autoregressive structures in a GLMM --
if your densities are high (means > 5) you can try this in MASS::glmmPQL.
However, I suspect you're unlikely to be able to resolve AR structure
in this data set -- if you're worried, take the residuals and examine
the correlation between r(t) and r(t+1) [you'll have 54 points -- one for
each plot/session combination, not counting the final session]. If that
correlation is small/non-significant, you're off the hook.
Quasipoisson fits are a little bit problematic in lme4 -- two alternatives
are (1) fitting a lognormal-Poisson model by adding a random effect with
a distinct level for each observation (2) fitting a negative binomial
model with glmm.admb (glmm.admb only allows a single random effect,
but you could do that if you make Colony fixed as suggested above).
Also keep in mind the rule of thumb that you need at least 10-20 observations
per parameter you want to estimate -- with 72 observations, you're pushing
it a little (intercept (1) + treatment (1) + session (3) + biomass (1) +
random effects (2)) (it's hard to know exactly how to count random effects,
but this is a crude rule of thumb). I certainly wouldn't try to make
your model any *more* complicated ...
good luck