为什么 lme 和 aov 对于 R 中的重复测量 ANOVA 返回不同的结果？

机器算法验证 r 方差分析 混合模式 重复测量 lme4-nlme
2022-01-17 16:45:08

我正在尝试从使用ez包转移到lme重复测量方差分析（因为我希望我能够在 with 上使用自定义对比lme）。

按照这篇博文的建议，我能够同时使用aov（如ez请求时一样）和lme. 但是，虽然在那篇文章中给出的示例中， F值确实在aov和之间完全一致lme（我检查了它，它们确实如此），但我的数据并非如此。尽管F值相似，但它们并不相同。

aov返回 f 值 1.3399，lme返回 1.36264。我愿意接受aov结果作为“正确”结果，因为这也是 SPSS 返回的结果（这对我的领域/主管很重要）。

问题：

  1. 如果有人能解释为什么存在这种差异以及我如何使用它lme来提供可信的结果，那就太好了。（如果它给出“正确”的结果，我也愿意使用lmer而不是lme这种类型的东西。但是，到目前为止我还没有使用它。）

  2. 解决这个问题后，我想进行对比分析。特别是我会对合并前两个因子水平（即c("MP", "MT")）的对比感兴趣，并将其与第三个因子水平（即 ）进行比较"AC"。此外，测试因子的第三与第四水平（即，"AC"与"DA"）。

数据：

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

和代码：

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))
2个回答

它们是不同的，因为 lme 模型迫使 的方差分量id大于零。查看所有项的原始方差分析表，我们看到 id 的均方误差小于残差的均方误差。 

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131

当我们计算方差分量时，这意味着由于 id 引起的方差将为负数。我对预期均方记忆的记忆不稳定，但计算类似于

(0.15052-0.16131)/3 = -0.003597.

这听起来很奇怪，但可能会发生。这意味着每个 id 的平均值比您预期的更接近，因为模型中的剩余变化量。

相反，使用 lme 会强制此方差大于零。

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

这报告标准偏差，平方以获得9.553e-10id 方差和0.1586164残差方差的方差收益。

现在，您应该知道，aov仅当您认为所有重复测量对之间的相关性相同时，才适合使用重复测量；这称为复合对称。（从技术上讲，球形度是必需的，但现在这已经足够了。）使用lmeover的一个原因aov是它可以处理不同类型的相关结构。

在这个特定的数据集中，这种相关性的估计是负的；这有助于解释 id 的均方误差如何小于残差平方误差。负相关意味着如果一个人的第一次测量值低于平均水平，平均而言，他们的第二次测量值将高于平均水平，这使得个体的总平均值比我们在零相关或正相关时预期的变化小。

使用lme随机效应等效于拟合复合对称模型，其中该相关性被迫为非负；我们可以使用以下方法拟合允许相关性为负的模型gls：

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

aov此 ANOVA 表与拟合表和拟合表一致lm。

好的，那又怎样？好吧，如果您认为id观测值之间的方差和相关性应该是非负的，那么lme拟合实际上比使用aov或拟合更合适，lm因为它对残差方差的估计略好。但是，如果您认为观察结果之间的相关性可能是负的，或者aov更好。 lmgls

您可能还对进一步探索相关结构感兴趣；要查看一般的相关结构，您可以执行类似的操作

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

这里我只将输出限制为相关结构。值 1 到 4 代表 的四个级别factor；我们看到因子 1 和因子 4 具有相当强的负相关性：

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

在这些模型之间进行选择的一种方法是使用似然比检验；这表明随机效应模型和一般相关结构模型在统计上没有显着差异；当这种情况发生时，通常首选更简单的模型。

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965

aov()lm()通过使用最小二乘拟合模型，lme通过最大似然拟合。估计线性模型参数的方式的差异可能是 f 值差异（非常小的）的原因。

在实践中，（例如假设检验）这些估计是相同的，所以我不明白如何认为一个比另一个“更可信”。它们来自不同的模型拟合范例。

对于对比，您需要为您的因子设置对比矩阵。Venebles 和 Ripley 在第 4 版的第 143 页、第 146 页和第 293-294 页中展示了如何做到这一点。

其它你可能感兴趣的问题
上一篇观测级别的马氏距离分布 下一篇“基准”是什么意思（在统计的背景下）？