[R-sig-ME] Understanding and analysing split split design using lmer()

[Peter Claussen]

> On Sep 9, 2019, at 1:35 PM, Rabin KC <kc000001 using umn.edu> wrote:
> 
> Hello R mixed models community,
> 
> …

> model1<aov(cover_coverage~rep+crop*tillage*cover+Error(rep/crop/tillage),data=data)
> 
> For my  analysis using lmer, and lme,  I have the following code:
> 
> model2<-lmer(biomass~crop*tillage*cover+ (1|rep),data=data)
> 
> model3<-lme(fixed=biomass~crop*tillage*cover, data=data, random=~1|rep)
> 
> Now my confusion and questions are:
> 1. I have used crop, tillage, and cover as fixed effects. And only the rep
> (block) as random. Is the random effect properly assigned for this design?

Simply, no. 

When you have three levels of independent randomization, then your design model has three random effects

Consider the expected output from agricolae. There should be an error term for crop, as randomized within whole blocks, designated Ea. This is the crop x replicate interaction term, and would also be the residual error if you were to analyze these data as a simple RCB of 2 treatments and 4 replicates (averaging biomass over whole plots). The F-test you would use to test the significance of crop would be MS(crop) / Ea, while the F-test for replicate effects would be MS(rep)/Ea.

There will be a second error strata, representing the independent randomization of tillage within whole plots (replicate:crop). That should be labelled Eb. This would also be the residual error of a simple split-plot experiment with crop as whole plot and tillage as subplot treatments (averaging biomass over subplots) and would be represented as rep:crop:tillage. This is the error term to test tillage and crop-tillage interaction.

The final error strata is the independent randomization of cover within crop:tillage subplots, and this will work out to be equivalent to residual error. So, to duplicate the decomposition of the aov from agricolae, we might write the linear model as

lm(biomass ~ crop*tillage*cover + rep +  rep:crop + rep:crop:tillage, data=data)
or, more briefly, as
lm(biomass ~ crop*tillage*cover + rep/crop/tillage,data=data)

That won’t give the correct F-tests, though, so to get aov to perform appropriate tests, we should write

aov(biomass ~ crop*tillage*cover + rep + Error(rep:crop/tillage), data=data)

(not what you’ve written, since we test rep against rep:crop)

This leads to an lmer model of
lmer(biomass ~ crop*tillage*cover+(1|rep/crop/tillage), data=data)

> 
> 2. I have read a lot about crossed and nested design. Crawley and Oehlert
> give particular examples about it, but I have trouble understanding if this
> design has nested factors or crossed factors.

crop, tillage and cover are all crossed factors, since each possible combination of crop, tillage and cover are included, and you might be interested in the interactions among these factors. Since this is a split-split-plot experiment, the experimental units (rep:crop = whole plot, rep:crop:tillage=subplot, rep:crop:tillage:cover = sub subplot) are nested, and not all combinations of treatments are subject to the same experimental errors.

That’s where the analysis starts to become tricky. If you want to compare treatments that are all applied within same random effects, that is,

Crop A : Tillage A : Cover A vs Crop A : Tillage A : Cover B, 

then the error terms for mean comparisons are simple. Comparisons across strata, like

Crop A : Tillage A : Cover A vs Crop B : Tillage A : Cover A

require estimates of the variances for each level of error, and with only a few observations (Ea (crop:rep) should have only 3 d.f.) you might not get estimates from lmer, and error terms derived from aov might include negative variances. AOV of the linear model (lm) is useful here, in that if any of the F-ratios of the random effects (rep, rep:crop or rep:crop:tillage) are less than 1, then you should expect lmer to produce a variance estimate of effectively 0. I would need to write out the expected means squares for a split-split-plot, though, to justify this statement.

By your email address and the subject matter, I’m guessing you’re in agronomy at the U of M. If so, I can suggest a friendly neighborhood statistician to consult.

Cheers,

Peter Claussen 
(work email Peter using gdmdata.com <mailto:Peter using gdmdata.com>)


	[[alternative HTML version deleted]]