Solved – Split-split-plot design and lme

experiment-designlme4-nlmersplit-plot

I’m working on a data set in order to evaluate the impact of drying on sediment microbial activities. The objective is to determine if the impact of drying varies among sediment types and/or depth within the sediment.

The experimental design is as follows:
The first factor Sediment corresponds to three types of sediment (coded Sed1, Sed2, Sed3).
For each type of Sediment, sampling was carried out on three sites (3 sites for Sed1, 3 sites for Sed2, 3 sites for Sed3). Site is coded : Site1, Site2, …, Site9.
The next factor is Hydrology: within each site, sampling is carried out in a dry plot and in a wet plot (coded Dry /Wet).
Within each of the previous plot, sampling is carried out at two Depths (D1, D2) in triplicate.

There are a total of n = 108 samples = 3 Sediment * 3 Sites * 2 Hydrology * 2 Depths * 3 Replicates.

I use the lme function in R (lnme package) as follows :

Sediment<-as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))
Site<-as.factor(rep(c("Site1","Site2","Site3","Site4","Site5","Site6","Site7","Site8","Site9"),each=12))
Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),9))
Depth<-as.factor(rep(rep(c("D1","D2"),each=3),18))
Variable<-rnorm(108)

mydata<-data.frame(Sediment,Site,Hydrology,Depth,Variable)

mod1<-lme(Variable~Sediment*Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth)

I found an example of a comparable split-split-plot design and its analysis in :
http://www3.imperial.ac.uk/portal/pls/portallive/docs/1/1171923.PDF

Could someone confirm that this is the right way to analyze these data?
Do you think that the random structure is correctly specified according to my experimental design?

Best Answer

This comes rather late, but I think your analysis is generally correct, with 3 comments:

  1. Make sure it is ok to treat depth as random rather than fixed. I guess this depends on your definition of depth. Is it just 'topsoil' and 'subsoil' or some control levels like A1, B2, C2, etc....
  2. A colleague of mine always recommend people creating another variable so it doesn't confuse people when the fixed and random effects have the same name. Something like lme(Variable~Sediment_ef*Hydrology_ef*Depth_ef, data=mydata, random=~1|Site/Hydrology/Depth), even if X_ef and X are identical columns.
  3. This is obviously the full model, where you might (or might not) want to reduce it to achieve parsimony.