lmerTest::anova 中的自由度是否正确？它们与 RM-ANOVA 非常不同

机器算法验证方差分析混合模式重复测量 lme4-nlme 自由程度

2022-03-19 13:48:24

我正在分析 R 中的反应时间实验的结果。

我运行了重复测量方差分析（1 个具有 2 个水平的受试者内因素和 1 个具有 2 个水平的受试者间因素）。我运行了一个类似的线性混合模型，我想用 ANOVA 表的形式总结 lmer 结果lmerTest::anova。

不要误会我的意思：我没想到会得到相同的结果，但是我不确定lmerTest::anova结果的自由度。~~在我看来，它更像是反映了一个在主题级别上没有聚合的方差分析。~~

我知道在混合效应模型中计算自由度很棘手，但在更新的主题（包）lmerTest::anova中被提及为一种可能的解决方案。?pvalueslme4

这个计算正确吗？结果lmerTest::anova是否正确反映了指定的模型？

更新：我使个体差异更大。的自由度lmerTest::anova与简单的方差分析有很大不同，但我仍然不确定，为什么它们对于主体内因素/交互作用如此之大。

# mini example with ANT dataset from ez package
library(ez); library(lme4); library(lmerTest)

# repeated measures ANOVA with ez package
data(ANT)
ANT.2 <- subset(ANT, !error)
# update: make individual differences larger
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
  within = .(direction), between = .(group))
anova.ez

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)

# simple ANOVA on all available data
m <- lm(rt ~ group * direction, data = ANT.2)
anova(m)

上面代码的结果[更新]：

方差分析

$方差分析

           Effect DFn DFd         F          p p<.05          ges
2           group   1  18 2.6854464 0.11862957       0.1294475137
3       direction   1  18 0.9160571 0.35119193       0.0001690471
4 group:direction   1  18 4.9169156 0.03970473     * 0.0009066868

lmerTest::anova（模型）

Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
                Df Sum Sq Mean Sq F value Denom Pr(>F)
group            1  13293   13293  2.6830    18 0.1188
direction        1   1946    1946  0.3935  5169 0.5305
group:direction  1  11563   11563  2.3321  5169 0.1268

方差分析（米）

Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1  1791568 1791568 242.3094 <2e-16 ***
direction          1      728     728   0.0985 0.7537    
group:direction    1    12024   12024   1.6262 0.2023    
Residuals       5187 38351225    7394                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2个回答

我认为这lmerTest是正确的，ezanova在这种情况下是错误的。

结果lmerTest符合我的直觉/理解
lmerTest（Satterthwaite 和 Kenward-Roger）中的两种不同的计算同意
他们也同意nlme::lme
当我运行它时，ezanova会给出一个警告，我不完全理解，但不应忽视它......

重新运行示例：

library(ez); library(lmerTest); library(nlme)
data(ANT)
ANT.2 <- subset(ANT, !error)
set.seed(101)  ## for reproducibility
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

弄清楚实验设计

with(ANT.2,table(subnum,group,direction))

所以看起来个体subnum（个人

(anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
    within = .(direction), between = .(group)))
## $ANOVA
##            Effect DFn DFd         F          p p<.05          ges
## 2           group   1  18 2.4290721 0.13651174       0.1183150147
## 3       direction   1  18 0.9160571 0.35119193       0.0002852171
## 4 group:direction   1  18 4.9169156 0.03970473     * 0.0015289914

在这里，我得到Warning: collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate. 分母 DF 看起来有点古怪（都等于 18）：我认为方向和组：方向应该更大，可以独立测试（但如果添加(direction|subnum)到模型中会更小）？

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
##                 Df  Sum Sq Mean Sq F value Denom Pr(>F)
## group            1 12065.7 12065.7  2.4310    18 0.1364
## direction        1  1952.2  1952.2  0.3948  5169 0.5298
## group:direction  1 11552.2 11552.2  2.3299  5169 0.1270

这里的Df列指的是分子df，Denom（倒数第二个）给出了估计的分母df；他们同意古典直觉。更重要的是，对于 F 值，我们也得到了不同的答案……

我们还可以用 Kenward-Roger 仔细检查（非常慢，因为它涉及多次重新拟合模型）

lmerTest::anova(model,ddf="Kenward-Roger")

结果是相同的。

对于这个例子lme（来自nlme包）实际上很好地猜测了适当的分母 df（F 值和 p 值略有不同）：

model3 <- lme(rt ~ group * direction, random=~1|subnum, data = ANT.2)
anova(model3)[-1,]
##                 numDF denDF   F-value p-value
## group               1    18 2.4334314  0.1362
## direction           1  5169 0.3937316  0.5304
## group:direction     1  5169 2.3298847  0.1270

direction如果我适合和subnumdf之间的交互direction并且group:direction更小（我原以为他们会是 18 岁，但也许我弄错了）：

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
lmerTest::anova(model2)
##                 Df  Sum Sq Mean Sq F value   Denom Pr(>F)
## group            1 20334.7 20334.7  2.4302  17.995 0.1364
## direction        1  1804.3  1804.3  0.3649 124.784 0.5469
## group:direction  1 10616.6 10616.6  2.1418 124.784 0.1459

我大体上同意 Ben 的分析，但让我补充几点意见和一点直觉。

一、总体结果：

使用 Satterthwaite 方法的 lmerTest 结果是正确的
Kenward-Roger 方法也是正确的，并且同意 Satterthwaite

Ben 概述了subnum嵌套在groupwhile中direction 并与 .group:direction交叉的设计subnum。这意味着自然误差项（即所谓的“封闭误差层”）group是subnum，而其他项（包括subnum）的封闭误差层是残差。

这种结构可以用所谓的因子结构图表示：

names <- c(expression("[I]"[5169]^{5191}),
           expression("[subnum]"[18]^{20}), expression(grp:dir[1]^{4}),
           expression(dir[1]^{2}), expression(grp[1]^{2}), expression(0[1]^{1}))
x <- c(2, 4, 4, 6, 6, 8)
y <- c(5, 7, 5, 3, 7, 5)
plot(NA, NA, xlim=c(2, 8), ylim=c(2, 8), type="n", axes=F, xlab="", ylab="")
text(x, y, names) # Add text according to ’names’ vector
# Define coordinates for start (x0, y0) and end (x1, y1) of arrows:
x0 <- c(1.8, 1.8, 4.2, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 7, 5, 5, 3, 7)
x1 <- c(2.7, 2.7, 5, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 7, 7, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length=0.1)

这里随机项用括号括起来，0代表总体平均值（或截距），[I]代表误差项，上标数字是水平数，下标数字是假设平衡设计的自由度数。该图表明是的自然误差项（封闭误差层）group和subnum的分子 dfsubnum等于的分母 dfgroup是 18: 20 减去 1 dfgroup和 1 df 的整体平均值。此处第 2 章提供了对因子结构图的更全面介绍：https ://02429.compute.dtu.dk/eBook 。

如果数据完全平衡，我们将能够从提供的 SSQ 分解构建 F 检验anova.lm。由于数据集非常接近平衡，我们可以获得如下近似 F 检验：

ANT.2 <- subset(ANT, !error)
set.seed(101)
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
fm <- lm(rt ~ group * direction + subnum, data=ANT.2)
(an <- anova(fm))
Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1   994365  994365 200.5461 <2e-16 ***
direction          1     1568    1568   0.3163 0.5739    
subnum            18  7576606  420923  84.8927 <2e-16 ***
group:direction    1    11561   11561   2.3316 0.1268    
Residuals       5169 25629383    4958                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这里所有的F和p值都是假设所有项都具有作为其封闭误差层的残差，并且对于除“组”之外的所有项都是如此。组的“平衡正确” F检验是：

F_group <- an["group", "Mean Sq"] / an["subnum", "Mean Sq"]
c(Fvalue=F_group, pvalue=pf(F_group, 1, 18, lower.tail = FALSE))
   Fvalue    pvalue 
2.3623466 0.1416875

我们在F值分母中使用subnumMS 而不是ResidualsMS 。

请注意，这些值与 Satterthwaite 结果非常匹配：

model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
anova(model, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

剩余的差异是由于数据不完全平衡。

OP 与比较anova.lm，anova.lmerModLmerTest这没关系，但要比较 like 和 like，我们必须使用相同的对比。在这种情况下，它们之间存在差异anova.lm，anova.lmerModLmerTest因为它们默认分别生成 I 类和 III 类测试，对于这个数据集，I 类和 III 类对比之间存在（小）差异：

show_tests(anova(model, type=1))$group
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1    0.005202759                     0.5013477

show_tests(anova(model, type=3))$group # type=3 is default
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1              0                           0.5

如果数据集完全平衡，则 I 型对比将与 III 型对比相同（不受观察到的样本数量的影响）。

最后一点是，Kenward-Roger 方法的“缓慢”不是由于模型重新拟合，而是因为它涉及使用观测值/残差的边际方差-协方差矩阵（在这种情况下为 5191x5191）进行计算，这不是萨特思韦特方法的例子。

关于模型2

至于模型2，情况变得更加复杂，我认为用另一个模型开始讨论会更容易，其中我已经包含了和之间的“经典”subnum交互direction：

model3 <- lmer(rt ~ group * direction + (1 | subnum) +
                 (1 | subnum:direction), data = ANT.2)
VarCorr(model3)
 Groups           Name        Std.Dev.  
 subnum:direction (Intercept) 1.7008e-06
 subnum           (Intercept) 4.0100e+01
 Residual                     7.0415e+01

因为与交互作用相关的方差基本上为零（在存在subnum随机主效应的情况下），交互项对分母自由度、F值和p值的计算没有影响：

anova(model3, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

但是，如果我们删除所有相关的subnum:directionSSQ会退回到subnumsubnumsubnum:direction

model4 <- lmer(rt ~ group * direction +
                 (1 | subnum:direction), data = ANT.2)

现在自然误差项group,direction和group:direction是 subnum:direction和nlevels(with(ANT.2, subnum:direction))= 40 和四个参数，这些项的分母自由度应该约为 36：

anova(model4, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
group           24004.5 24004.5     1 35.994  4.8325 0.03444 *
direction          50.6    50.6     1 35.994  0.0102 0.92020  
group:direction   273.4   273.4     1 35.994  0.0551 0.81583  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这些F检验也可以用“平衡正确”的F检验来近似：

an4 <- anova(lm(rt ~ group*direction + subnum:direction, data=ANT.2))
an4[1:3, "F value"] <- an4[1:3, "Mean Sq"] / an4[4, "Mean Sq"]
an4[1:3, "Pr(>F)"] <- pf(an4[1:3, "F value"], 1, 36, lower.tail = FALSE)
an4
Analysis of Variance Table

Response: rt
                   Df   Sum Sq Mean Sq F value Pr(>F)    
group               1   994365  994365  4.6976 0.0369 *  
direction           1     1568    1568  0.0074 0.9319    
group:direction     1    10795   10795  0.0510 0.8226    
direction:subnum   36  7620271  211674 42.6137 <2e-16 ***
Residuals        5151 25586484    4967                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

现在转向模型2：

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)

该模型描述了一个相当复杂的具有 2x2 方差-协方差矩阵的随机效应协方差结构。默认参数化不容易处理，我们最好对模型进行重新参数化：

model2 <- lmer(rt ~ group * direction + (0 + direction | subnum), data = ANT.2)

如果我们比较model2，model4它们有同样多的随机效应；每个 2 个subnum，即总共 2*20=40。虽然model4为所有 40 个随机效应规定了单个方差参数，但model2规定每subnum对随机效应具有双变量正态分布，具有 2x2 方差-协方差矩阵，其参数由下式给出

VarCorr(model2)
 Groups   Name           Std.Dev. Corr 
 subnum   directionleft  38.880        
          directionright 41.324   1.000
 Residual                70.405

这表明过度拟合，但让我们将其留到另一天。这里的重点是是model4的特例model2 ，model也是的特例model2。松散（和直观）地说，(direction | subnum)包含或捕捉与主要效果subnum 以及交互相关的变化direction:subnum。就随机效应而言，我们可以将这两种效应或结构视为分别捕获行和逐列之间的变化：

head(ranef(model2)$subnum)
  directionleft directionright
1    -25.453576     -27.053697
2     16.446105      17.479977
3    -47.828568     -50.835277
4     -1.980433      -2.104932
5      5.647213       6.002221
6     41.493591      44.102056

在这种情况下，这些随机效应估计以及方差参数估计都表明我们实际上只有一个随机主效应subnum（行间变化）。这一切导致的是 Satterthwaite 分母的自由度

anova(model2, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
group           12059.8 12059.8     1  17.998  2.4329 0.1362
direction        1803.6  1803.6     1 125.135  0.3638 0.5475
group:direction 10616.6 10616.6     1 125.136  2.1418 0.1458

是这些主效应和交互结构之间的折衷：DenDF 组保持在 18（subnum设计嵌套），但direction和 group:directionDenDF 是 36 ( model4) 和 5169 ( model) 之间的折衷。

我认为这里没有任何迹象表明 Satterthwaite 近似值（或其在lmerTest中的实现）是错误的。

Kenward-Roger 方法的等效表给出

anova(model2, type=1, ddf="Ken")
Type I Analysis of Variance Table with Kenward-Roger's method
                 Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
group           12059.8 12059.8     1 18.000  2.4329 0.1362
direction        1803.2  1803.2     1 17.987  0.3638 0.5539
group:direction 10614.7 10614.7     1 17.987  2.1414 0.1606

KR 和 Satterthwaite 可能不同并不奇怪，但出于所有实际目的，p值的差异很小。我上面的分析表明DenDFfordirection和group:direction不应小于 ~36 并且可能大于考虑到我们基本上只有当前的随机主效应direction，所以如果有的话，我认为这表明 KR 方法变得DenDF太低了在这种情况下。但请记住，数据并不真正支持(group | direction)结构，因此比较有点人为 - 如果模型真的得到支持会更有趣。

其它你可能感兴趣的问题

上一篇评估回归模型下一篇为什么在标准偏差公式中对样本计数“N”取平方根？