[R-sig-ME] nested random effects with temporal correlation

[Thierry Onkelinx]

Dear Charlotte,

The reviewer is correct, at least from a conceptual point of view.
Pair and transect are both random effects and transect is nested in
pair. Temporal auto correlation is plausible.

However from a practical point of view, you don't have enough pairs to
use it as a random effects (see
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#should-i-treat-factor-xxx-as-fixed-or-random).
The number of transects is a tiny bit larger. Therefore I would use
only the transects as random intercepts. Any effect at the pairs
levels will be absorbed by the transects.

Your glmmPQL model estimates the **residual** auto correlation
**within** the most detailed random effect level. The residuals from
year t are correlated to the residuals from year t-1 **within** the
same transect. Residuals **between** different transects are assumed
to be independent. Given that your model contains the treatment - year
interaction and year is used as a factor, then much of the temporal
pattern will be already explained by the fixed effects. Hence the
residual auto correlation might be overkill.

Another option is to model Year as a random effect with temporal auto
correlation between the years. You can do that with the INLA package.

Best regards,

Thierry



ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given body
of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////




2018-06-04 13:38 GMT+02:00 Charlotte Reemts <creemts using tnc.org>:
> Mixed modelers,
> I am analyzing data from a fire experiment. We counted the number of point intercepts of a grass along 3 pairs of transects (so n=6). One transect in each pair was burned. We collected data in 2006 (before burning) and for three years after burning. We calculated frequency of the grass by dividing the intercepts by the total number of contacts of all species, so I am using the binomial family for analysis. We want to know whether the frequency of the grass decreased in the burned transects compared to the unburned transects and how that frequency changes with time since burning (I am using glht contrasts, not shown, to do this).
>
> A reviewer would like me to include nested random effects and account for temporal correlation, since we repeatedly sampled the same transects. I created the glmmPQL model below, which converges, but I am concerned that it is a very complex model for a small dataset. Do you have any suggestions for 1) the best way to simplify the model and 2) the best way to justify that simplification to the reviewer?
>
> Thanks,
> Charlotte
>
> krdata2<-read.table(header=T, text= "
> Pair
>
> Treatment
>
> Year
>
> Transect
>
> totalcontacts
>
> freq
>
> 1
>
> burned
>
> 2006
>
> T1
>
> 190
>
> 0.778947
>
> 1
>
> burned
>
> 2007
>
> T1
>
> 231
>
> 0.337662
>
> 1
>
> burned
>
> 2008
>
> T1
>
> 250
>
> 0.508
>
> 1
>
> burned
>
> 2009
>
> T1
>
> 148
>
> 0.52027
>
> 1
>
> unburned
>
> 2006
>
> C1
>
> 188
>
> 0.946809
>
> 1
>
> unburned
>
> 2007
>
> C1
>
> 210
>
> 0.92381
>
> 1
>
> unburned
>
> 2008
>
> C1
>
> 214
>
> 0.878505
>
> 1
>
> unburned
>
> 2009
>
> C1
>
> 162
>
> 0.962963
>
> 2
>
> burned
>
> 2006
>
> T2
>
> 196
>
> 0.80102
>
> 2
>
> burned
>
> 2007
>
> T2
>
> 270
>
> 0.414815
>
> 2
>
> burned
>
> 2008
>
> T2
>
> 266
>
> 0.56015
>
> 2
>
> burned
>
> 2009
>
> T2
>
> 210
>
> 0.847619
>
> 2
>
> unburned
>
> 2006
>
> C2
>
> 193
>
> 0.782383
>
> 2
>
> unburned
>
> 2007
>
> C2
>
> 211
>
> 0.85782
>
> 2
>
> unburned
>
> 2008
>
> C2
>
> 194
>
> 0.938144
>
> 2
>
> unburned
>
> 2009
>
> C2
>
> 198
>
> 0.959596
>
> 3
>
> burned
>
> 2006
>
> T3
>
> 193
>
> 0.632124
>
> 3
>
> burned
>
> 2007
>
> T3
>
> 275
>
> 0.192727
>
> 3
>
> burned
>
> 2008
>
> T3
>
> 222
>
> 0.405405
>
> 3
>
> burned
>
> 2009
>
> T3
>
> 176
>
> 0.642045
>
> 3
>
> unburned
>
> 2006
>
> C3
>
> 198
>
> 0.747475
>
> 3
>
> unburned
>
> 2007
>
> C3
>
> 207
>
> 0.758454
>
> 3
>
> unburned
>
> 2008
>
> C3
>
> 207
>
> 0.772947
>
> 3
>
> unburned
>
> 2009
>
> C3
>
> 143
>
> 0.944056
>
> ")
>
> krdata2$Year<-as.factor(krdata2$Year)
> aus.kr.pql2<-glmmPQL(freq ~ Treatment*Year, random = ~1|Pair/Treatment,
>                     correlation=corCAR1(form = ~Year|Pair/Treatment),
>                     family=binomial(link="logit"), weights=totalcontacts, data=krdata2)?
>
>
>
> ___________________________________________
> Charlotte Reemts, M.S.
> Research and Monitoring Ecologist
> creemts using tnc.org
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models