5

我想在方差分析（固定或混合模型）中进行单 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全因子，三个苜蓿品种和4个1943年第三次扦插日期。实验单元分为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

因此，方差分析输出看起来像：

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

我该怎么做？.....................

4

1 回答 1

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正如我所说的，对于除, （可能没有尝试过）aov和某些生存模型之外的模型，p 值将是错误的。这是因为默认情况下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 的大量阅读的提炼。我可能在很多方面都错了，如果我错了，希望有更博学的人能纠正我们俩。
于 2013-07-25T20:19:52.703 回答