I have data from an experiment with the following design:
- 11 blocks within a forest along an elevational gradient
- Within each block, three treatment plots (control, low and high treatment) with unique elevations
- From each plot we assessed the height growth of three trees
- Sampling was performed annually during 10 years
Here is an illustration of the design:
And here some data:
year block treatment tree elevation growth
1 2001 b1 C a 10.3 5.34
2 2001 b1 C b 10.3 7.12
3 2001 b1 C c 10.3 4.62
4 2001 b1 L a 12.2 4.33
5 2001 b1 L b 12.2 5.01
6 2001 b1 L c 12.2 6.35
7 2001 b1 H a 14.4 7.11
8 2001 b1 H b 14.4 6.54
9 2001 b1 H c 14.4 4.24
10 2001 b2 C a 1.5 5.98
11 2001 b2 C b 1.5 4.45
12 2001 b2 C c 1.5 5.64
13 2001 b2 L a 3.1 5.96
14 2001 b2 L b 3.1 5.44
.
.
.
Goal: I want to estimate the effects of treatment, elevation and year and their interaction. To do so, I'm trying to fit a mixed model to the data.
My approach: I specified the following model using nlme in R:
M1 <- lme(growth ~ treatment * elevation * year, random = ~1|block/treatment/tree)
Problem: I have limited statistical knowledge. I'm afraid that the random effect of block confounds the fixed effect of elevation, as each block covers a distinct range of elevations. Could that be the case?
Is this a valid way of specifying the random effects given the design?
Best Answer
Here's what I would do. I expanded your dataset a bit so we have something to work with:
Now I am not super familiar with
lme()
so I am usinglmer()
from thelme4
package. Given your design, I would removeyear
as a fixed effect. It has 10 levels, and since it is a growth experiment, it will probably be significant because trees grow over time. Making sense of a model output with 3 main effects and their interactions whenyear
has 10 levels might be a bit too much. That reduces your model totreatment
andblock
(I wouldn't useelevation
unless these small differences inelevation
within the plots were purposely chosen).As for the random part, I would include
year
(to account for multiple measurements on the same trees) andblock:treatment
, which represents the plots, because the 3 trees in the plots may not be spatially independent(?!), which should be accounted for.To get p-values for the fixed effects you could use the
mixed()
function from theafex
package:However, note that this may not be the only way to analyze your dataset. Depending on the actual questions you have, you can build multiple models and compare them. To generally get a better idea about fixed and random effects, you could start reading here.