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 at tnc.org
nested random effects with temporal correlation
2 messages · Charlotte Reemts, Thierry Onkelinx
2 days later
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 at 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 at 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 at tnc.org
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models