1

将 16 个处理 (4*2*2) 的因子组合重复 3 次，并布置在条带分割块中。处理包括八个整地处理（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 函数。我设法使用 aov 获得了正确的结果 function result ~ Burning * Weeding * Site.prep + Error(Rep/Burning*Site.prep)。使用 lmer 我使用了这个函数 result ~ Burning*Site.prep*Weeding+(1|Rep/(Burning:Site.prep))，只给了我部分正确的结果。Burning、Site.prep 和 Burning:Site.prep 的 SS 值和 F 值与 Genstat 结果有偏差（尽管不是太多），但 Weeding 和 Weeding 相互作用给出了与 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)
4

1 回答 1

0

应用@cuttlefish44 提到的站点中的三向裂区因子方差分析示例，会导致：

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）。我不确定这是否是建模方法差异的结果：多层次方法同时考虑了内部和效应之间的影响。

于 2017-07-30T20:04:47.277 回答