I am having troubles coming up with the proper statistical model and R code for a forestry trial we did several years ago. We looked at whether wood piles left over in clear cuts were beneficial for biodiversity. I include a schematic of the design:enter image description here

We looked at wood piles left over in forest stands 2,4,and 6 years post harvest. Within each stand we examined 8 wood piles, that we sampled using 3 traps at different locations, 1 in the center of the wood pile, 1 at the edge of the wood pile and 1 on the ground between wood piles. I understand that this is not the optimal design for an experiment, but in an operational forestry landscape it is often the best that can be done.

What I am trying to analyze is the catch of insects (y), with special attention to the Harvest age*Trap Location interaction. I want to stress that it is not a repeated measures experiment, as each of the time periods post harvest were sampled in the same year on different forest stands. I am interested in using generalized linear mixed effects modelling in R (package lme4). I am certain that the underlying distribution of the data is a negative binomial (which is very common distribution for insect data). My model summary is as follows:

Total number of samples = 144, therefore total df = 143;

Harvest Age (A) - Fixed effect - N=3 - df=2;
Trap Location (L) - Fixed effect - N=3 - df=2;
A*L - Fixed effect - (A-1)(L-1) - df=4;

I know that both Stands (S; N=2, df=1) and wood piles (P; N=8, df=7) are both Random effects. Forest stands could be thought of as blocks (replicates of Harvest Age)? I think wood piles are the experimental unit, and are nested within the ALS interaction, and would probably be the total error term in the model. I am stuck trying to figure out how to account for the variation among stands and wood piles. Any help with the model and R-code would be greatly appreciated.


1 Answer 1


Ok, I think the design is closest to a Split-Split plot design. ANOVA table below.

enter image description here

However, I have been trying to code different models into a simple ANOVA in r to check the model degrees of freedom, but cannot get the correct error terms:

aovtest = aov(abund ~ yr + loc*yr + Error(yr*blk + loc*pil + loc*yr*pil), data=abund)  

I have also tried other variations on this theme but cannot get the right mode. Anyone out there have a suggestion?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.