在分析纵向数据集时使用 aov() 和 lme() 有什么区别？

机器算法验证 r 混合模式 重复测量 面板数据
2022-02-17 10:30:35

谁能告诉我使用aov()和lme()分析纵向数据之间的区别以及如何解释这两种方法的结果？

下面，我使用aov()和分析相同的数据集，lme()并得到 2 个不同的结果。使用aov()，我在治疗交互作用的时间上得到了显着的结果，但拟合线性混合模型，治疗交互作用的时间是微不足道的。

> UOP.kg.aov <- aov(UOP.kg~time*treat+Error(id), raw3.42)
> summary(UOP.kg.aov)

Error: id
          Df  Sum Sq Mean Sq F value Pr(>F)
treat      1   0.142  0.1421  0.0377 0.8471
Residuals 39 147.129  3.7725               

Error: Within
            Df  Sum Sq Mean Sq  F value  Pr(>F)    
time         1 194.087 194.087 534.3542 < 2e-16 ***
time:treat   1   2.077   2.077   5.7197 0.01792 *  
Residuals  162  58.841   0.363                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> UOP.kg.lme <- lme(UOP.kg~time*treat, random=list(id=pdDiag(~time)), 
                    na.action=na.omit, raw3.42)
> summary(UOP.kg.lme)
Linear mixed-effects model fit by REML
 Data: raw3.42 
       AIC      BIC    logLik
  225.7806 248.9037 -105.8903

Random effects:
 Formula: ~time | id
 Structure: Diagonal
        (Intercept)      time  Residual
StdDev:   0.6817425 0.5121545 0.1780466

Fixed effects: UOP.kg ~ time + treat + time:treat 
                 Value Std.Error  DF   t-value p-value
(Intercept)  0.5901420 0.1480515 162  3.986059  0.0001
time         0.8623864 0.1104533 162  7.807701  0.0000
treat       -0.2144487 0.2174843  39 -0.986042  0.3302
time:treat   0.1979578 0.1622534 162  1.220053  0.2242
 Correlation: 
           (Intr) time   treat 
time       -0.023              
treat      -0.681  0.016       
time:treat  0.016 -0.681 -0.023

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.198315285 -0.384858426  0.002705899  0.404637305  2.049705655 

Number of Observations: 205
Number of Groups: 41
3个回答

根据您的描述，您似乎有一个具有单一处理因素的重复测量模型。由于我无权访问数据集 ( raw3.42)，因此我将使用nlme包中的 Orthodont 数据来说明此处发生的情况。数据结构是相同的（对两个不同的组进行重复测量——在这种情况下是男性和女性）。

如果您运行以下代码：

library(nlme)
data(Orthodont)

res <- lme(distance ~ age*Sex, random = ~ 1 | Subject, data = Orthodont)
anova(res)

你会得到以下结果：

            numDF denDF  F-value p-value
(Intercept)     1    79 4123.156  <.0001
age             1    79  122.450  <.0001
Sex             1    25    9.292  0.0054
age:Sex         1    79    6.303  0.0141

如果你运行：

res <- aov(distance ~ age*Sex + Error(Subject), data = Orthodont)
summary(res)

你会得到：

Error: Subject
          Df Sum Sq Mean Sq F value   Pr(>F)   
Sex        1 140.46 140.465  9.2921 0.005375 **
Residuals 25 377.91  15.117                    

Error: Within
          Df  Sum Sq Mean Sq  F value  Pr(>F)    
age        1 235.356 235.356 122.4502 < 2e-16 ***
age:Sex    1  12.114  12.114   6.3027 0.01410 *  
Residuals 79 151.842   1.922

请注意，F 检验完全相同。

对于lme()，您使用list(id=pdDiag(~time))，它不仅为模型添加了随机截距，还添加了随机斜率。此外，通过使用pdDiag，您将随机截距和斜率之间的相关性设置为零。这是与您指定的模型不同的模型aov()，因此您会得到不同的结果。

我只想补充一点，您可能希望安装该car软件包并使用Anova()该软件包提供的软件包，而不是anova()因为 foraov()和lm()objects，vanillaanova()使用顺序平方和，这对于不相等的样本大小会给出错误的结果，而lme()它使用任何一种类型-I 或类型 III 的平方和取决于type论点，但类型 III 的平方和违反了边际性——即它对待交互作用与主效应没有区别。

R-help 列表对 I 型和 III 型平方和没有什么好说的，但这些是唯一的选择！去搞清楚。

编辑：实际上，如果有一个重要的交互项，那么 type-II 似乎是无效的，而且似乎任何人都可以做的最好的事情就是在有交互时使用 type-III。我通过回答我自己的一个问题得到了提示，而这些问题又将我引向了这篇文章。

在我看来，您每次对每个 id 都有多个度量值。您需要为 aov 汇总这些信息，因为它不公平地增加了该分析的能力。我并不是说进行汇总会使结果相同，但它应该使它们更加相似。

dat.agg <- aggregate(UOP.kg ~ time + treat + id, raw3.42, mean)

然后像用 dat.agg 替换数据之前一样运行您的 aov 模型。

另外，我相信 anova(lme) 更适合您比较结果。效应的方向和大小与模型方差与误差的比率不同。

（顺便说一句，如果你对聚合数据进行 lme 分析，你不应该这样做，并检查 anova(lme) 你会得到与 aov 几乎相同的结果）

其它你可能感兴趣的问题
上一篇R中的梯度下降与lm（）函数？ 下一篇R中的面向对象编程教程