So I've done a lot of reading and chatting to people and I have a solution.
My experimental design is a split plot design, which is quite different from a nested or hierarchical design. I was originally confusing the terms. As Robert correctly states in his answer, what is needed is a mixed effects model. Thus:
Fixed effects: Year, Treatment1, Treatment2
Random effects: Year, Block, Treatment1
The model is specified thus:
mod<- lmer(Richness~Treatment1*Treatment2*Year+(1|Block/Treatment1)+(1|Year),data=dat,poisson)
The fixed effects are the terms specified in the brackets. Since none of these are continuous (the effect of Year doesn't necessarily increase each year in a linear fashion so I have classed it as a categorical fixed effect), they are specified 1|fixed effect, where 1 represents the intercept.
If Block were actually a continuous fixed effect (obviously hypothetical!) then the fixed effects might be specified +(Block|Treatment1)+(1|Year).
The model can then be simplified as appropriate.
Several things to note:
1) When specified as a random effect, Year is listed separately from Block and Treatment1, since it doesn't have an intuitive "level" at which to be nested between them (Year isn't any different at any plot size of the experiment: for every block, plot and subplot Year is the same.
2) Treatment 2 does not need to be specified as a random effect since it represents the highest level of replication in the experiment and therefore will not be psuedoreplicated
3) In mixed effects models it is possible to specify an error distribution if errors are not normal. I have specified poisson here, since my response data are counts - this improved the distribution of the model residuals.
Best Answer
Assuming that each subject is only in one group, you have a nested design. Conceptually, it makes more sense to treat group as a fixed effect. As you have only three groups, it wouldn't make much sense statistically speaking anyways because you'd be asking the software to estimate a variance for
group
assuming a normal distribution from only three observations.Your sample size is quite low which might be a problem. Nevertheless, I suggest starting with a simple random intercept model:
This model fits a global intercept which is simply the intercept for the reference group, deviations from that intercept for the remaining groups, a single slope for the effect of
day
and a random intercept forsubject
. Hence, this model assumes that each group has the same trajectory over time, the same slope but different intercepts. To allow each group their own slope, you could fit the following model with an interaction betweenday
andgroup
:Lastly, you could also allow for random slopes:
All these models assume a linear relationship between the
outcome
andday
which might be unrealistic. If you suspect nonlinear relationships, you could easily accomodate those by using polynomials or (my recommendation) restricted cubic splines (aka natural splines) with, say, 3 knots. Finally, have a look at this post which goes into more details about such longitudinal models.