r中的anova分区和比较（正交单df）

Question

我想在方差分析（固定或混合模型）中进行单 df 正交对比。这里只是一个例子：

require(nlme)
data (Alfalfa)
  Variety: a factor with levels Cossack, Ladak, and Ranger
  Date : a factor with levels None S1 S20 O7
  Block: a factor with levels 1 2 3 4 5 6
  Yield : a numeric vector

这些数据在 Snedecor 和 Cochran (1980) 中作为示例进行了描述 的裂区设计。实验中使用的处理结构 是一个 3×4 全因子，有三个品种的紫花苜蓿和四个 1943年第三次切割的日期。安排了实验单元 分成六个街区，每个街区又分为四个地块。苜蓿的品种 （Cossac、Ladak 和 Ranger）被随机分配到区块和 第三次切割的日期（无，S1-9 月 1 日，S20-9 月 20 日， 和 O7—10 月 7 日）被随机分配到这些地块。 每个区块都使用了所有四个日期。

model<-with (Alfalfa, aov(Yield~Variety*Date +Error(Block/Date/Variety)))

    > summary(model)

Error: Block
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  5   4.15    0.83

Error: Block:Date
          Df Sum Sq Mean Sq F value   Pr(>F)
Date       3 1.9625  0.6542   17.84 3.29e-05 ***
Residuals 15 0.5501  0.0367
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Block:Date:Variety
             Df Sum Sq Mean Sq F value Pr(>F)
Variety       2 0.1780 0.08901   1.719  0.192
Variety:Date  6 0.2106 0.03509   0.678  0.668
Residuals    40 2.0708 0.05177

我想进行一些比较（组内的正交对比），例如对于日期，两个对比：

   (a) S1 vs others (S20 O7)
   (b) S20 vs 07,

对于品种因素二对比：

  (c)  Cossack vs others (Ladak and Ranger)
   (d) Ladak vs Ranger

因此方差分析输出如下所示：

Error: Block
              Df Sum Sq Mean Sq F value Pr(>F)
Residuals  5   4.15    0.83

Error: Block:Date
          Df Sum Sq Mean Sq F value   Pr(>F)
Date       3 1.9625  0.6542   17.84 3.29e-05 ***
       (a) S1 vs others    ?        ?
       (b)  S20 vs 07       ?        ?
Residuals 15 0.5501  0.0367
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Block:Date:Variety
             Df Sum Sq Mean Sq F value Pr(>F)
Variety       2 0.1780 0.08901   1.719  0.192
     (c)  Cossack vs others ?        ?    ?
     (d)  Ladak vs Ranger    ?       ?     ?
Variety:Date  6 0.2106 0.03509   0.678  0.668
Residuals    40 2.0708 0.05177

我该如何执行此操作？ ....................

Answer 1

首先，为什么要使用方差分析？您可以使用nlme 包中的lme，除了aov 为您提供的假设检验之外，您还可以获得效应大小和效应方向的可解释估计值。无论如何，我想到了两种方法：

• 手动指定变量的对比，如here 所述。
• 安装multcomp 包并使用glht。

glht 对于预测变量中的多元模型有点固执己见。不过，长话短说，如果您要创建一个对角矩阵cm0，其尺寸和暗名与模型的vcov 相同（假设它是一个lme 拟合，称为model0），那么summary(glht(model0,linfct=cm0))应该给出与summary(model0)$tTable 相同的估计、SE 和测试统计数据（但不正确 p 值）。现在，如果你弄乱了来自cm0 的行的线性组合并创建具有与cm0 相同的列数但这些线性组合作为行的新矩阵，你最终会找出创建矩阵的模式为您提供每个单元格的截距估计值（对照predict(model0,level=0) 检查）。现在，另一个具有该矩阵各行之间差异的矩阵将为您提供相应的组间差异。可以使用相同的方法但将数值设置为 1 而不是 0 来获取每个单元格的斜率估计值。然后可以使用这些斜率估计之间的差异来获得组间斜率差异。

记住三件事：

• 正如我所说，除了lm、（可能没有尝试过）aov 和某些生存模型之外，p 值将是错误的。这是因为glht 默认采用z 分布而不是t 分布（lm 除外）。要获得正确的 p 值，请使用测试统计 glht 计算并手动执行 2*pt(abs(STAT),df=DF,lower=F) 以获得双尾 p -value，其中 STAT 是 glht 返回的测试统计信息，DF 是 df 来自 summary(model0)$tTable 中相应类型的默认对比。
• 您的对比可能不再检验独立假设，如果还没有，则需要进行多次检验校正。通过 p.adjust 运行 p 值。
• 这是我自己从教授和同事的大量手稿中提炼出来的，以及对相关主题的 Crossvalidated 和 Stackoverflow 的大量阅读。我可能在很多方面都错了，如果我错了，希望有更博学的人能纠正我们俩。