R lmer中的拆分块

Question

16个处理(4*2*2)的析因组合被重复三次，并在条裂区块中布置。处理包括8个样地准备(4×2)作为整个小区处理和两个水平的除草(除草/不除草)随机应用于子样地。在Genstat中运行分析，得到以下结果：

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    

我想在R中重复这些结果。我使用了base::aov函数和lmerTest::lmer函数。我设法使用函数result ~ Burning * Weeding * Site.prep + Error(Rep/Burning*Site.prep)获得了aov的正确结果。对于lmer，我使用了函数result ~ Burning*Site.prep*Weeding+(1|Rep/(Burning:Site.prep))，只给出了部分正确的结果。燃烧、Site.prep和Burning:Site.prep的SS值和F值与Genstat结果的偏差(尽管不是太大)，但除草和除草交互作用产生的SS和F值与Genstat输出的SS和F值相同。我想知道我应该如何指定lmer模型来重现Genstat和aov结果。数据和代码如下：

    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)

Answer 1

应用@cuttlefish44 44提到的site中的三向分裂图阶乘方差分析示例，可以得出：

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

除了Site.prep之外，结果都是匹配的。此外，lmer()和lme()之间的结果非常相似(对于Site.prep也是如此)。我不确定这是否是建模方法差异的结果:多级方法既考虑了内部影响，也考虑了影响之间的影响。