Type II SS 在混合模型中应该如何计算？

机器算法验证 r 方差分析 混合模式 平方和
2022-04-01 07:57:37

lmerTest我有一个数据集（和相应的混合模型），当使用 I 型（顺序，注意它是最后一个）和 II 型（使用和car包）进行测试时，它在双向交互之一中获得了非常不同的 p 值。

顺序和car包“同意”并且lmerTest包是不同的。

这些方法在做什么，为什么它们不同，哪些是“正确的”？

（请注意，这里有很多关于 Type I/II/III 的问题，但我没有看到回答这个关于 Type II 的特定问题。另外，当我使用 R 并且对 R 中使用的特定实现感兴趣时，更大的问题是统计问题，而不是编码问题。）

library(lme4)
library(lmerTest)
library(car)

d <- data.frame(
  Y = c(6, 4, 5, 7, 6, 8, 9, 0, 10, 9, 10, 8, 7, 7, 6, 5, 6, 7, 8, 
        7, 9, 10, 10, 6, 5, 8, 10, 4, 6, 7, 7, 10, 6, 10, 10, 7, 6, 3, 
        6, 7, 5, 7, 8, 9, 11, 10, 7, 8, 5, 6, 6, 7, 8, 5, 4, 8, 8, 8, 
        7, 9, 5, 4, 4, 6, 5, 7, 8, 3, 9, 8, 8, 7, 5, 7, 4, 5, 5, 4, 6, 
        8, 8, 6, 7, 8, 4, 2, 4, 4, 3, 6, 7, 8, 9, 8, 7, 6, 4, 2, 3, 3, 
        3, 0, 1, 4, 6, 4, 2, 0, 9, 9, 10, 8, 10, 9, 7, 9, 6, 5, 9, 7, 
        6, 2, 2, 3, 5, 7, 8, 9, 9, 8, 7, 8, 1, 0, 3, 2, 2, 5, 3, 5, 8, 
        6, 6, 6, 8, 6, 2, 3, 4, 5, 7, 6, 7, 8, 3, 3),
  A = rep(c(2, 5, 4, 4, 2, 3, 3, 5, 9, 3, 7, 2, 11), each=12),
  B = factor(rep(1:2, each=6)),
  C = factor(1:6),
  ID = factor(rep(1:13, each=12))
)
options(contrasts=c("contr.sum","contr.poly")) # for Type III car::Anova

这是我的模型，我注意到 B:C 交互的差异。

m1 <- lmer(Y ~ A * B * C + (1|ID), data=d)

I 型方差分析（来自lmerTest），除了 A:B:C 之外，B:C 相互作用最后，给出 B:C 的 p=0.04。（但是，顺序似乎并不重要。）

anova(m1, type=1, ddf="Kenward-Roger")
#> Analysis of Variance Table of type I  with  Kenward-Roger 
#> approximation for degrees of freedom
#>        Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A       8.794   8.794     1    11   2.908   0.11621    
#> B     136.641 136.641     1   121  45.178 6.269e-10 ***
#> C      14.128   2.826     5   121   0.934   0.46141    
#> A:B     0.381   0.381     1   121   0.126   0.72330    
#> A:C    45.874   9.175     5   121   3.034   0.01291 *  
#> B:C    34.821   6.964     5   121   2.303   0.04882 *  
#> A:B:C  21.860   4.372     5   121   1.446   0.21285

包装中的II 型 Anovacar给出了相同的结果。

Anova(m1, type=2, test="F")
#> Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
#> 
#> Response: Y
#>             F Df Df.res    Pr(>F)    
#> A      2.9075  1     11   0.11621    
#> B     45.1784  1    121 6.269e-10 ***
#> C      0.9343  5    121   0.46141    
#> A:B    0.1259  1    121   0.72330    
#> A:C    3.0335  5    121   0.01291 *  
#> B:C    2.3026  5    121   0.04882 *  
#> A:B:C  1.4456  5    121   0.21285

来自 的 II 型 Anova lmerTest，对于 B:C 给出 p=0.88，并且在 和 中都lmerTest与III 型相同car。

anova(m1, type=2, ddf="Kenward-Roger")
#> Analysis of Variance Table of type II  with  Kenward-Roger 
#> approximation for degrees of freedom
#>       Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A      8.794   8.794     1    11  2.9075 0.1162113    
#> B     41.500  41.500     1   121 13.7214 0.0003209 ***
#> C     50.879  10.176     5   121  3.3645 0.0070072 ** 
#> A:B    0.381   0.381     1   121  0.1259 0.7233027    
#> A:C   45.874   9.175     5   121  3.0335 0.0129149 *  
#> B:C    5.174   1.035     5   121  0.3422 0.8864025    
#> A:B:C 21.860   4.372     5   121  1.4456 0.2128521   

Anova(m1, type=3, test="F")
#> Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
#> 
#> Response: Y
#>                   F Df Df.res    Pr(>F)    
#> (Intercept) 88.2434  1     11 1.376e-06 ***
#> A            2.9075  1     11 0.1162113    
#> B           13.7214  1    121 0.0003209 ***
#> C            3.3645  5    121 0.0070072 ** 
#> A:B          0.1259  1    121 0.7233027    
#> A:C          3.0335  5    121 0.0129149 *  
#> B:C          0.3422  5    121 0.8864025    
#> A:B:C        1.4456  5    121 0.2128521    

anova(m1, type=3, ddf="Kenward-Roger")
#> Analysis of Variance Table of type III  with  Kenward-Roger 
#> approximation for degrees of freedom
#>       Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A      8.794   8.794     1    11  2.9075 0.1162113    
#> B     41.500  41.500     1   121 13.7214 0.0003209 ***
#> C     50.879  10.176     5   121  3.3645 0.0070072 ** 
#> A:B    0.381   0.381     1   121  0.1259 0.7233027    
#> A:C   45.874   9.175     5   121  3.0335 0.0129149 *  
#> B:C    5.174   1.035     5   121  0.3422 0.8864025    
#> A:B:C 21.860   4.372     5   121  1.4456 0.2128521
2个回答

II 类测试car::Anova与 lmerTest 的anova方法之间的区别在于如何处理连续解释变量。详细信息部分的第一个步骤help(Anova)描述了如何Anova处理它们：

“type-II”和“type-III”的名称是从 SAS 借用的，但这里使用的定义与 SAS 使用的定义并不完全对应。Type-II测试是根据边缘性原则计算的，每个term都先测试，除了忽略term的高阶亲属；所谓的 III 型测试违反了边际性，将模型中的每个项都测试在所有其他项之后。II 类检验的定义对应于 SAS 为方差分析模型生成的检验，其中所有预测变量都是因子，但不是更普遍的（即，当有定量预测变量时）。为 III 型测试制定模型时要非常小心，否则测试的假设将毫无意义。

因此，虽然 lmerTest 实现了本文中描述的 Type II 的 SAS 定义，car::Anova但做了一些不同的事情。此处描述了 SAS 定义，这意味着对一个术语的测试对于所有不包含它的术语都是边缘化的。包含的定义在SAS 文档中给出为：

给定两个效应F1和F2，如果满足以下两个条件，F1则称为包含在：F2

  • 两种效应都涉及相同的连续变量（如果有）。
  • F2有更多的 CLASS [R 中的因素] 变量F1，如果F1有 CLASS 变量，它们都出现在F2

所以在你的模型B:C中不包含，A:B:C因为A是连续的，B:C因此 II 型假设是边际的A:B:C。car::Anova在这种情况下不区分因子和协变量，因此如果A是一个因子car::Anova并且 lmerTest-anova同意：

d$A2 <- cut(d$A, quantile(d$A), include.lowest = TRUE)
m2 <- lmer(Y ~ A2 * B * C + (1|ID), data=d)
anova(m2, type=2, ddf="Ken") 
Type II Analysis of Variance Table with Kenward-Roger's method
        Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
A2      17.210   5.737     3     9  1.7035   0.23533    
B      136.641 136.641     1    99 40.5756 5.987e-09 ***
C       14.128   2.826     5    99  0.8391   0.52512    
A2:B    11.137   3.712     3    99  1.1024   0.35193    
A2:C    38.816   2.588    15    99  0.7684   0.70879    
B:C     34.821   6.964     5    99  2.0680   0.07575 .  
A2:B:C  50.735   3.382    15    99  1.0044   0.45709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova(m2, type=2, test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: Y
             F Df Df.res    Pr(>F)    
A2      1.7035  3      9   0.23533    
B      40.5756  1     99 5.987e-09 ***
C       0.8391  5     99   0.52512    
A2:B    1.1024  3     99   0.35193    
A2:C    0.7684 15     99   0.70879    
B:C     2.0680  5     99   0.07575 .  
A2:B:C  1.0044 15     99   0.45709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

该lmerTest::show_tests()函数可用于准确查看模型系数的哪个线性函数构成了测试假设。例如，我们有

show_tests(anova(m1, type=2))$`B:C`
      (Intercept) A B1 C1 C2 C3 C4 C5 A:B1 A:C1 A:C2 A:C3 A:C4 A:C5 B1:C1 B1:C2 B1:C3
B1:C1           0 0  0  0  0  0  0  0    0    0    0    0    0    0     1     0     0
B1:C2           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     1     0
B1:C3           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     1
B1:C4           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     0
B1:C5           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     0
      B1:C4 B1:C5 A:B1:C1 A:B1:C2 A:B1:C3 A:B1:C4 A:B1:C5
B1:C1     0     0       0       0       0       0       0
B1:C2     0     0       0       0       0       0       0
B1:C3     0     0       0       0       0       0       0
B1:C4     1     0       0       0       0       0       0
B1:C5     0     1       0       0       0       0       0

所以这种对比对于模型中的所有其他项来说都是微不足道的。因为car::AnovaType II 测试B:C等价于 Type I 测试，我们可以看到这种对比可以表示为：

show_tests(anova(m1, type=1), fractions = TRUE)$`B:C`
      (Intercept) A     B1    C1    C2    C3    C4    C5    A:B1  A:C1  A:C2  A:C3 
B1:C1     0           0     0     0     0     0     0     0     0     0     0     0
B1:C2     0           0     0     0     0     0     0     0     0     0     0     0
B1:C3     0           0     0     0     0     0     0     0     0     0     0     0
B1:C4     0           0     0     0     0     0     0     0     0     0     0     0
B1:C5     0           0     0     0     0     0     0     0     0     0     0     0
      A:C4  A:C5  B1:C1 B1:C2 B1:C3 B1:C4 B1:C5 A:B1:C1 A:B1:C2 A:B1:C3 A:B1:C4 A:B1:C5
B1:C1     0     0     1   1/2   1/2   1/2   1/2 60/13   30/13   30/13   30/13   30/13  
B1:C2     0     0     0     1   1/3   1/3   1/3     0   60/13   20/13   20/13   20/13  
B1:C3     0     0     0     0     1   1/4   1/4     0       0   60/13   15/13   15/13  
B1:C4     0     0     0     0     0     1   1/5     0       0       0   60/13   12/13  
B1:C5     0     0     0     0     0     0     1     0       0       0       0   60/13

这对于 3 向交互项来说并不重要。

不同之处在于A不是以零为中心。尚不确定为什么这很重要，或者是否应该（似乎不应该...）；其他更全面地探讨这一点的答案是最受欢迎的。

d$Ac <- d$A - mean(d$A)
m2 <- lmer(Y ~ Ac*B*C + (1|ID), data=d)

anova(m2, type=1)
#> Type I Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285    

anova(m2, type=2)
#> Type II Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285    

anova(m2, type=3)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285
其它你可能感兴趣的问题
上一篇我是否需要丢弃 90% 的实验以使样本独立？ 下一篇在现代深度神经网络中使用 MSE 损失有什么意义？