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:


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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.