Mil*_*los 2 r linear-regression mixed-models
我有六个固定因子:A, B, C, D, EandF和一个随机因子R。我想使用语言 R 测试线性项、纯二次项和双向交互。因此,我构建了完整的线性混合模型并尝试使用以下方法测试其项drop1:
full.model <- lmer(Z ~ A + B + C + D + E + F
+ I(A^2) + I(B^2) + I(C^2) + I(D^2) + I(E^2) + I(F^2)
+ A:B + A:C + A:D + A:E + A:F
+ B:C + B:D + B:E + B:F
+ C:D + C:E + C:F
+ D:E + D:F
+ E:F
+ (1 | R), data=mydata, REML=FALSE)
drop1(full.model, test="Chisq")
Run Code Online (Sandbox Code Playgroud)
似乎drop1完全忽略了线性项:
Single term deletions
Model:
Z ~ A + B + C + D + E + F + I(A^2) + I(B^2) + I(C^2) + I(D^2) +
I(E^2) + I(F^2) + A:B + A:C + A:D + A:E + A:F + B:C + B:D +
B:E + B:F + C:D + C:E + C:F + D:E + D:F + E:F + (1 | R)
Df AIC LRT Pr(Chi)
<none> 127177
I(A^2) 1 127610 434.81 < 2.2e-16 ***
I(B^2) 1 127378 203.36 < 2.2e-16 ***
I(C^2) 1 129208 2032.42 < 2.2e-16 ***
I(D^2) 1 127294 119.09 < 2.2e-16 ***
I(E^2) 1 127724 548.84 < 2.2e-16 ***
I(F^2) 1 127197 21.99 2.747e-06 ***
A:B 1 127295 120.24 < 2.2e-16 ***
A:C 1 127177 1.75 0.185467
A:D 1 127240 64.99 7.542e-16 ***
A:E 1 127223 48.30 3.655e-12 ***
A:F 1 127242 66.69 3.171e-16 ***
B:C 1 127180 5.36 0.020621 *
B:D 1 127202 27.12 1.909e-07 ***
B:E 1 127300 125.28 < 2.2e-16 ***
B:F 1 127192 16.60 4.625e-05 ***
C:D 1 127181 5.96 0.014638 *
C:E 1 127298 122.89 < 2.2e-16 ***
C:F 1 127176 0.77 0.380564
D:E 1 127223 47.76 4.813e-12 ***
D:F 1 127182 6.99 0.008191 **
E:F 1 127376 201.26 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Run Code Online (Sandbox Code Playgroud)
如果我从模型中排除交互:
full.model <- lmer(Z ~ A + B + C + D + E + F
+ I(A^2) + I(B^2) + I(C^2) + I(D^2) + I(E^2) + I(F^2)
+ (1 | R), data=mydata, REML=FALSE)
drop1(full.model, test="Chisq")
Run Code Online (Sandbox Code Playgroud)
然后测试线性项:
Single term deletions
Model:
Z ~ A + B + C + D + E + F + I(A^2) + I(B^2) + I(C^2) + I(D^2) +
I(E^2) + I(F^2) + (1 | R)
Df AIC LRT Pr(Chi)
<none> 127998
A 1 130130 2133.9 < 2.2e-16 ***
B 1 130177 2181.0 < 2.2e-16 ***
C 1 133464 5467.6 < 2.2e-16 ***
D 1 129484 1487.9 < 2.2e-16 ***
E 1 130571 2575.0 < 2.2e-16 ***
F 1 128009 12.7 0.0003731 ***
I(A^2) 1 128418 422.2 < 2.2e-16 ***
I(B^2) 1 128193 197.4 < 2.2e-16 ***
I(C^2) 1 129971 1975.1 < 2.2e-16 ***
I(D^2) 1 128112 115.6 < 2.2e-16 ***
I(E^2) 1 128529 533.0 < 2.2e-16 ***
I(F^2) 1 128017 21.3 3.838e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Run Code Online (Sandbox Code Playgroud)
因为这是drop1工作方式(它不是特定于混合模型 - 您会发现这种行为也适用于拟合的常规线性模型lm)。来自?drop1:
在考虑要添加或删除的项时,要考虑层次结构:二阶交互中包含的所有主效应必须保留,依此类推。
我在这篇 CrossValidated 帖子中详细讨论了这一点
统计上棘手的部分是,在还包含更高级别交互的模型中测试较低级别的交互(取决于您与谁交谈)要么(i)难以正确执行,要么(ii)只是很愚蠢(对于后一种情况) ,请参阅 Bill Venables 的“线性模型注释”的第 5 部分)。这方面的准则是边缘性原则。至少,低阶项的含义敏感地取决于模型中的对比度是如何编码的(例如处理与中点/总和为零)。我的默认规则是,如果您不确定自己是否完全理解这可能是一个问题的原因,则不应违反边缘性原则。
但是,正如 Venables 在链接文章中实际描述的那样,如果您愿意,您可以让 R 违反边缘性(第 15 页):
令我高兴的是,我看到因子项之间的边际约束在默认情况下得到尊重,并且学生不会沿着逻辑上滑溜溜的“类型 III 平方和”路径走下去。我们讨论了为什么没有显示主效应,它是一个有用的教程点。
具有讽刺意味的是,当然,如果人们了解它们的真正含义以及如何获得它们,那么 III 型平方和一直可用。如果调用
drop1包含任何式作为第二个参数,对应于所有非截距项的模型矩阵的部分被省略逐一从模型中,给某种测试的用于主效应...如果您使用了带有零和列的对比矩阵,它们将是唯一的,它们就是臭名昭著的“III 型平方和”。但是,如果您使用
contr.treatment对比,那么列的总和不为零,您就会变得无意义。这种对在这种情况下应该是任意的事情的敏感性应该足以提醒任何人注意正在做一些愚蠢的事情。
换句话说,使用scope = . ~ .将迫使drop1忽略边缘性。这样做的风险由您自己承担 - 当您遵循此程序时,您绝对应该能够向自己解释您实际测试的内容......
例如:
set.seed(101)
dd <- expand.grid(A=1:10,B=1:10,g=factor(1:10))
dd$y <- rnorm(1000)
library(lme4)
m1 <- lmer(y~A*B+(1|g),data=dd)
drop1(m1,scope=.~.)
## Single term deletions
##
## Model:
## y ~ A * B + (1 | g)
## Df AIC
## <none> 2761.9
## A 1 2761.7
## B 1 2762.4
## A:B 1 2763.1
Run Code Online (Sandbox Code Playgroud)