split block in R lmer

Question

A factorial combination of 16 treatments (4*2*2) was replicated three times and laid out in a strip-split block. Treatments consisted of eight site preparations (4*2) applied as whole plot treatments and two levels of weeding(weeding/no-weeding) were applied randomly to subplots. The analysis was run in Genstat giving the following results:

Variate: result

Source of variation d.f.    s.s.    m.s.    v.r.    F pr.

Rep stratum            2     35.735  17.868 
Rep.Burning stratum
Burning                1     0.003   0.003   0.00    0.972
Residual               2     3.933   1.966   1.53    
Rep.Site_prep stratum
Site_prep              3     7.981   2.660   0.45    0.727
Residual               6     35.477  5.913   4.61    
Rep.Burning.Site_prep stratum
Burning.Site_prep      3     2.395   0.798   0.62    0.626
Residual               6     7.691   1.282   0.60   
Rep.Burning.Site_prep.*Units* stratum
Weeding                1     13.113  13.113  6.13    0.025
Burning.Weeding        1     0.486   0.486   0.23    0.640
Site_prep.Weeding      3     17.703  5.901   2.76    0.076
Burning.Site_prep.Weed.3     3.425   1.142   0.53    0.666
Residual               16    34.248  2.141       
Total                  47    162.190

I want to repeat these results in R. I used both the base::aov function and the lmerTest::lmer function. I managed to get the correct results with aov using function result ~ Burning * Weeding * Site.prep + Error(Rep/Burning*Site.prep). With lmer I used the function result ~ Burning*Site.prep*Weeding+(1|Rep/(Burning:Site.prep)) giving me only partially correct results. The SS values and the F-values for Burning, Site.prep and Burning:Site.prep deviated (although not too much)from the Genstat results, but the Weeding and Weeding interactions gave the same SS and F-valus as the Genstat output. I would like to know how I should specify the lmer model to reproduce the Genstat and aov results. Data and code below:

    x <- structure(list(
  Rep = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", "2", "3"
 ), class = "factor"),Burning = structure(c(1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L), .Label = c("Burn", 
"No-burn"), class = "factor"), Site.prep = structure(c(4L, 4L,4L, 4L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 
 2L, 2L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L), 
 .Label = c("Chop_Pit", "Chop_Rip", "Pit", "Rip"), class = "factor"), Weeding = structure(c(1L, 
2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L),
.Label = c("Weedfree",  "Weedy"), class = "factor"), 
Dbh14 = c(27.4, 28.4083333333333, 27.7066666666667, 27.3461538461538, 28.6, 28.3333333333333, 27.0909090909091, 
27.8076923076923, 27.1833333333333, 27.5461538461538, 24.3076923076923, 
29.3461538461538, 27.4, 25.1, 26.61, 28.0461538461538, 27.71, 
25.2533333333333, 25.3833333333333, 24.2307692307692, 24.2533333333333, 
24.95, 24.34375, 26.9909090909091, 24.775, 25.9076923076923, 
25.1666666666667, 25.9933333333333, 27.0466666666667, 30.5625, 
27.36, 25.2636363636364, 29.6846153846154, 27.7, 28.3071428571429, 
29.4857142857143, 27.025, 30.1, 31.2454545454545, 24.2888888888889, 
28.4875, 29.23, 30, 28.5, 29.3615384615385, 27.45, 28.8153846153846, 
29.1866666666667)), .Names = c("Rep", "Burning", "Site.prep", 
"Weeding", "result"), class = "data.frame", row.names = c(NA, -48L))

model1 <- aov(result ~ Burning* Weeding*Site.prep+ Error(Rep/Burning*Site.prep), data=x)
summary(model1)

model2 <- lmer(result ~ Burning*Site.prep*Weeding+(1|Rep/(Burning:Site.prep)),data=x)
anova(model2)

[This site](http://www.uni-kiel.de/psychologie/rexrepos/posts/anovaMixed.html) will give you the answer. — cuttlefish44, Sep 11 '16 at 09:59

score 0 · Answer 1 · answered Jul 30 '17 at 20:04

Applying the three-way split-plot-factorial ANOVA example from the site mentioned by @cuttlefish44, leads to:

library(lme4)
library(nlme)
m1 <- aov(result ~ Weeding*Burning*Site.prep + Error(Rep/Burning*Site.prep), data=x)
m2 <- lmer(result ~ Weeding*Burning*Site.prep + (1|Rep) + (1|Burning:Rep) +
               (1|Site.prep:Rep), data=x)
m3 <- anova(lme(result ~ Weeding*Burning*Site.prep,
          random=list(Rep=pdBlocked(list(~1, pdIdent(~Burning-1), pdIdent(~Site.prep-1)))),
          method="ML", data=x))
summary(m1)
anova(m2)
m3

Except for Site.prep, the results match. Moreover, the results between lmer() and lme() are pretty similar (also for Site.prep). I'm not sure whether this is the result of differences in modelling approaches: the multi-level approach takes both within and between effects into account.

Just a set of parentheses needed to make `m1` match exactly: `aov(result ~ Weeding*Burning*Site.prep + Error(Rep/(Burning*Site.prep)), data=x)` — Aaron left Stack Overflow, Aug 09 '19 at 01:05

split block in R lmer

1 Answers1