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

[Peter Claussen]

Rabin,

OK, finally getting back to this.

A couple cases to consider that might be leading to a singular fit. The first should be, I think, the most obvious - you simply do not have enough samples from a random variable to obtain an estimate. If we can assume that your experimental replicates and whole plots sample over comparable spatial heterogeneity, then combining trials will give you more samples from those populations, thus perhaps a better estimate for whole pilot variance terms.

The other is that your analysis model is incorrect for these experiments. 

Remember that one of the assumptions of AOV is that residual errors are independent and identically distributed. This assumption allows us to compute variances from expected mean squares. I would suggest https://dl.sciencesocieties.org/publications/aj/abstracts/81/4/AJ0810040665 <https://dl.sciencesocieties.org/publications/aj/abstracts/81/4/AJ0810040665> for an example of how to write up the expected mean squares. 

We can compute variance components from mean squares, based on their expectations, but these usually dependent on an estimate of residual variance, and when residuals are not i.i.d., I think we get poor estimates of the other variance terms.

To explore this, I considered the data you posted and fit the model

model1 <-lm(biomass ~ crop*tillage*cover + rep/crop/tillage,data=data)
par(mfrow=c(2,2))
plot(model1)

We might need to accept that the residuals are not i.i.d, i particular, that the errors variances are not comparable between corn and soybeans. In lmer, we might attempt to address this with

model1.lmer <- lmer(biomass ~ crop*tillage*cover+(1|rep/crop/tillage), data=data)
model2.lmer <- lmer(biomass ~ crop*tillage*cover+(crop | rep/crop/tillage), data=data)
plot(model1.lmer)
plot(model2.lmer)

There appears to be some improvement, at least with soybeans. 

You might consider what this suggests about your experimental design and randomization model. The primary plots may have been exchangeable units when either corn or soybean were seeded, but over the course of the season, may no have been exchangeable with regard to cover crop - that is, errors in application or germination may be different when seeding into corn vs soybeans. Perhaps there is an interaction with tillage?

Another issue you might be facing is the assumption of independence. Consider these:

data$wholeplot <- data$crop:data$rep
data$subplot <- data$wholeplot:data$tillage
data$residuals <- residuals(model1)
plot(residuals ~ rep,data=data)
plot(residuals ~ wholeplot,data=data)

There may be a spatial component that introduces correlations among the residuals, making the estimate of residual unreliable. You would need to consider your trial map. It might be possible to model a fixed effect spatial trend and continue with lmer; capturing residual correlations in two dimensions (i.e. AR1 x AR1), as far as I know, is only possible in R using AS-REML.


You most certainly will want to combine experiments, but for this post I’m limiting myself to a simpler discussion of computational options for mixed models in R. 

FWIW,  I presented some options for mixed model analysis of agricultural trials at one of the ASA/CSSA/SSSA meetings - see https://gdmdata.com/media/documents/handouts/2016ASA_RMixedModelsPresentation.pdf <https://gdmdata.com/media/documents/handouts/2016ASA_RMixedModelsPresentation.pdf> . I also gave a presentation on our current approach to spatial modeling - the presentation (https://gdmdata.com/media/documents/handouts/2017ASA_BeyondRCBD_Spatial.pdf <https://gdmdata.com/media/documents/handouts/2017ASA_BeyondRCBD_Spatial.pdf>) doesn’t include R code but there is supplementary material with literate documentation at https://gdmdata.com/media/documents/handouts/2017ASA_BeyondRCBD.zip <https://gdmdata.com/media/documents/handouts/2017ASA_BeyondRCBD.zip> . I left autocorrelated errors out of this presentation because I was unable to find an acceptable solution in lme. 

Cheers




> On Sep 9, 2019, at 6:05 PM, Rabin KC <kc000001 using umn.edu> wrote:
> 
> Dear Peter, 
> 
> Your reply has made my day and very clearly answered the questions I was asking myself for a few weeks. Thank you for correcting my mistakes as well. 
> 
> The lmer() does produce a variance of 0 for random effect (crop: rep). I now know the term is called boundary (singular) fit. What would be the correct thing to do in this case? I do have this experiment conducted for 2 years in 2 different locations, which might take care of the precious degrees of freedom, but again, will add new variables making the model more complicated. 
> 
> I ran an Anova II (wald test) and was thinking of doing post hoc analysis on significant terms. (If this would take care of singularity). Any advice on that would be greatly appreciated.
> 
> And yes, I am in Agronomy in U of M. I will be writing to you in your work email about the friendly statistician you mentioned!!!!
> 
> 
> Best regards, 
> Rabin
> 
> On Mon, Sep 9, 2019 at 4:43 PM Peter Claussen <dakotajudo using mac.com <mailto:dakotajudo using mac.com>> wrote:
> 
> 
>> On Sep 9, 2019, at 1:35 PM, Rabin KC <kc000001 using umn.edu <mailto: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]]