Using lme to analyse a complete randomized block design with repeated measures: Is my model correct?

Question

I've got a completely randomized block design with three treatments and four replications. Biodiversity was measured in four successive years.

I figured that a mixed model with repeated measures as random terms should be appropriate to analyse this design.

My hypothesis is that considering all years, biodiversity is different between the treatments.

This is my analysis:

library(nlme)
library(multcomp)
# Made-up random dataset
mydata <- data.frame(
  Treatment=rep(c("Control", "Irrigation", "Fertailization"), 16), 
  Block=rep(1:4, 12), 
  Year=rep(2000:2003, 12), 
  Value=runif(48, 0.5, 1.5)
)
# Model Treatment is a fixed effect, Year is a random effect
fit <- lme(Value ~ Treatment,  random = ~1|Year, data = mydata)
# Post-hoc comparison
summary(glht(fit,linfct=mcp(Treatment="Tukey")))

My questions:

1. Is my model correct?

2. Is the post-hoc comparison appropriate?

3. How could I include the bock effect?

---

If I understood you right, "Year" should be nested in "Block" - so the correct model would be coded like this:

fit <- lme(Value ~ Treatment,  random = ~1|Block/Year, data = mydata)

There seems to be a linear temporal trend in the data. How could I account for this in the model?

Answer 1

This is the model I might start with:

fit <- lme(Value ~ Treatment * Year,  random = ~1|Block, data = mydata)

I would include the year as a fixed effect, since a temporal trend of biodiversity usually can be expected and it would also be of interest. However, this is guesswork, because I don't know the background of the experiment nor the actual data. Wether the temporal effect can be assumed to be linear and wether you need the interaction, you would have to judge from your data. Block is clearly a random effect here and is needed to account for repeated measures.

Usually with this experimental setup you have treatment plots within your blocks, which often stay the same over the whole measurement period. Then it could be necessary to account for that, too:

fit <- lme(Value ~ Treatment * Year,  random = ~1|Block/Plot, data = mydata)

You did not mention, what kind of variable Value is. You might need a generalized linear model (look a the family parameter of lme) or need to transform your dependend.