I am seeking advice on structuring the random-effects component of a mixed model in R. The data that I am modeling come from an ecological experiment with a repeatedly-measured split-plot design. At each of 5 sites (site = A, B, C, D, E) experimenters set up 2 treatments (x = treatment, control) in separate plots. A response variable (y) was measured every quarter (quarter = Q1, Q2, Q3, Q4) through time for 8 years (year = 1–8). There is no replication of treatments within sites (thus each site contains 1 treatment plot and 1 control plot). Hence this is a split-plot design (site / treatment) with repeated measures.
Some fake data:
dat <- expand.grid('year' = 2001:2008, 'quarter' = paste("Q", 1:4, sep = ""), 'site' = LETTERS[1:5], 'treatment' = c('control', 'manipulate')) dat$yr.qt = dat$year + as.numeric(ifelse(dat$quarter == "Q1", 0.125 + 0, ifelse(dat$quarter == "Q2", 0.125 + 0.25, ifelse(dat$quarter == "Q3", 0.125 + 0.5, 0.125 + 0.75)))) dat$time.step <- ((dat$yr.qt - min(dat$yr.qt)) / 0.25) + 1 dat$y <- sample(seq(0, 10, by = 0.01), size = nrow(dat), replace = TRUE) # Structure of the data xtabs(~ site + time.step, data = dat)
The focal question is, "What is the effect of the treatment on the response?". A secondary question is, "Does the treatment effect change over time?"
Assume for the moment that sites were selected randomly from a range of potential sites and that the experiments are not concerned with site-specific trends. Also assume that we only want to control for seasonality ('quarter') as a covariate.
After reading Pinheiro and Bates (2000) and Zuur et al. (2009), I have come to two potential ways to structure the random effects of a mixed model using nlme (I have preferred using nlme over lme4 for the moment so that I can explore structuring the variance to deal with heterogeneity or serial autocorrelation). These models are:
library('nlme') mod1 <- lme(y ~ treatment + quarter + time.step + treatment:time.step, random = ~ 1 | site, data = dat) mod2 <- lme(y ~ treatment + quarter + time.step + treatment:time.step, random = ~ 1 | site/treatment, data = dat)
Considering that I have no replication of treatments within sites, which is the more appropriate model structure? Examination of the ANOVA tables for 'mod1' vs. 'mod2' shows 309 vs. 4 denominator degrees of freedom, respectively. Is 'mod1' pseudo-replicating the data? Or, is 'mod2' overdetermined?
Many thanks for your feedback!