r中的anova分区和比较(正交单df)

Question

我想在anova(固定或混合模型)中进行单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年有3个品种的苜蓿和4个第三次扦插日期.实验单元分为6个区块,每个区块分为4个区块.将苜蓿(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

因此anova输出看起来像:

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

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

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

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

要记住三件事：

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