机器算法验证 - 对于在 lmer+lmerTest 中具有交互作用的混合模型，summary() 和 anova() 的结果冲突 - 吾爱随笔录

我最近遇到了一个问题，我认为包中的anova()函数如何lmerTest计算其 F 统计量和来自混合效应模型的固定效应的相应 P 值。首先让我说我知道围绕从混合效应模型计算 P 值的争议（原因在这里讨论）。尽管如此，许多人仍然想要 P 值，因此已经开发了许多方法来适应这种情况（参见此处）。在这里，我想展示一种常用方法的结果——即包中的anova函数lmerTest——并希望有人知道为什么结果不太有意义。

首先是我的数据。由于它的大小，我不得不链接到它。请注意，生物量列已标准化（平均值 = 0，sd = 1），因此为负值。这不会改变输出。下载并指定工作目录后，可以按如下方式读取文件：

dat <- read.csv("StackOverflow_Data.csv", header = T)

下面是我使用lmer函数的模型lme4。在这个模型中，我将植物生物量作为响应变量和三个因素——A、B 和 C——每个都有两个水平作为预测变量。植物基因型和空间块作为随机效应包括在内。

model <- lmer(Biomass ~ A + B + C + 
            A:B + A:C + 
            B:C + A:B:C +
            (1 | Genotype) + (1 | Block) , 
          data = dat, REML = T)

总结上述模型，summary(model)我们得到：

Linear mixed model fit by maximum likelihood t-tests use Satterthwaite approximations

  to degrees of freedom [lmerMod]
Formula: Biomass ~ A + B + C + A:B + A:C + B:C + A:B:C + (1 | Genotype) +  
    (1 | Block)
   Data: dat

 AIC      BIC   logLik deviance df.resid 
  1059.7   1111.0   -518.8   1037.7      776 
Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.04330 -0.63914  0.00315  0.69108  2.82368 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.07509  0.2740  
 Block    (Intercept) 0.01037  0.1018  
 Residual             0.19038  0.4363  
Number of obs: 787, groups:  Genotype, 50; Block, 6

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           2.27699    0.08162  47.50000  27.897  < 2e-16 ***
AYes                 -0.02308    0.09958  99.30000  -0.232  0.81719    
BReduced             -0.11036    0.06232 733.00000  -1.771  0.07700 .  
CSupp                -0.02152    0.06243 733.70000  -0.345  0.73039    
AYes:BReduced         0.25113    0.08838 733.70000   2.841  0.00462 ** 
AYes:CSupp            0.02179    0.08854 734.50000   0.246  0.80567    
BReduced:CSupp        0.19436    0.08838 733.10000   2.199  0.02817 *  
AYes:BReduced:CSupp  -0.21746    0.12507 734.20000  -1.739  0.08251 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) AYes   BRedcd CSupp  AYs:BR AYs:CS BRd:CS
AYes        -0.607                                          
BReduced    -0.379  0.311                                   
CSupp       -0.379  0.311  0.498                            
AYes:BRedcd  0.269 -0.444 -0.706 -0.354                     
AYes:CSupp   0.268 -0.444 -0.352 -0.708  0.503              
BRedcd:CSpp  0.267 -0.219 -0.706 -0.705  0.498  0.500       
AYs:BRdc:CS -0.190  0.315  0.500  0.502 -0.709 -0.709 -0.708

上面的摘要使用该lmerTest包使用 Satterthwaites 对分母自由度的近似值从 t 统计量计算 P 值。由此我们看到，在 p = 0.05 的水平上，A:B和交互作用都显着。B:C从理论上讲，这些结果至少在质量上应该与包anova()中的函数产生的结果一致，该函数以lmerTest相同的方式计算 P 值。然而事实并非如此；这是来自的输出anova(model, type = 3)。注意 SS 的类型 3 测试

Analysis of Variance Table of type III  with  Satterthwaite 
approximation for degrees of freedom
       Sum Sq Mean Sq NumDF  DenDF F.value  Pr(>F)  
A     0.09492 0.09492     1  49.87  0.4986 0.48342  
B     0.66040 0.66040     1 732.66  3.4688 0.06294 .
C     0.20207 0.20207     1 733.90  1.0614 0.30324  
A:B   0.99470 0.99470     1 732.56  5.2247 0.02255 *
A:C   0.36903 0.36903     1 733.66  1.9383 0.16427  
B:C   0.35867 0.35867     1 733.20  1.8839 0.17031  
A:B:C 0.57552 0.57552     1 734.23  3.0230 0.08251 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这些结果明显不同。B:C交互作用不再显着，交互作用的 P 值要高A:B得多。两种模型都应该以相似的方式计算 P 值，因此很难想象它们会有如此不同。

为什么它们不同？

这是原始问题的一部分，但它可能会产生误导，请参阅下面的答案。

实际上，看起来该anova(model, type = 3)功能实际上是在使用 2 型 SS，我们可以通过运行来验证anova(model, type = 2)。

Analysis of Variance Table of type II  with  Satterthwaite 
approximation for degrees of freedom
       Sum Sq Mean Sq NumDF  DenDF F.value  Pr(>F)  
A     0.09526 0.09526     1  49.87  0.5004 0.48263  
B     0.65996 0.65996     1 732.66  3.4665 0.06302 .
C     0.19639 0.19639     1 733.91  1.0315 0.31013  
A:B   0.99282 0.99282     1 732.56  5.2148 0.02268 *
A:C   0.37018 0.37018     1 733.65  1.9444 0.16362  
B:C   0.35523 0.35523     1 733.20  1.8659 0.17237  
A:B:C 0.57552 0.57552     1 734.23  3.0230 0.08251 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结果非常相似，考虑到模型中存在交互，情况不应如此。为了证明它lmerTest::anova()实际上使用的是 2 型 SS，而不是它在输出中显示的 3 型 SS，我们可以使用包中的Anova()函数car。Anova(model, type = 2, test.statistic = 'F')产生：

Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: Biomass
           F Df Df.res  Pr(>F)  
A     0.4857  1  48.28 0.48917  
B     3.4537  1 726.63 0.06351 .
C     1.0337  1 727.77 0.30962  
A:B   5.1456  1 726.54 0.02360 *
A:C   1.9302  1 727.55 0.16517  
B:C   1.8776  1 727.12 0.17103  
A:B:C 2.9915  1 728.06 0.08413 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

请注意，对于我的数据，使用 Kenward-Roger ddf 并没有太大改变结果。很清楚的是，2 型 SS 来自于Car封装的生产结果，类似于 3 型 SS 来自于lmerTest包装的结果。这表明该lmerTest包实际上是计算类型 2 SS。我很难弄清楚为什么会出现这种情况，除非从lmerTest包中计算 P 值时出现问题。我错过了什么吗？

欢迎任何建议或想法。非常感谢！

编辑：2016 年 12 月 6 日，上午 11:40

一些人表示这个问题是从这里重复的。但是我不明白这是怎么回事。那篇文章旨在了解为什么aov()并lme()产生不同的 F 统计量，结果证明这与如何从不同的函数计算方差分量有关。在这里，我只运行一个模型，使用并lmer试图理解为什么会产生不同的 P 值，尽管它们应该以类似的方式计算。似乎使用的是 2 型 SS 而不是报告的 3 型 SS，这仅在存在交互的情况下才重要，另一篇文章未包含在任何列出的模型中。lmerTest::anova(model)summary(model)lmerTest::anova()