# Mixed model random-effect structure for repeated measures, split plot design in R (nlme or lme4)

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?