r-如何在线性混合效应模型中找到随机效应的p值

Question

我在R:中运行以下代码行

model = lme(divedepth ~ oarea, random=~1|deployid, data=GDataTimes, method="REML")
summary(model)

我看到了这个结果：

Linear mixed-effects model fit by REML
Data: GDataTimes 
AIC     BIC   logLik
2512718 2512791 -1256352
Random effects:
Formula: ~1 | deployid
(Intercept) Residual
StdDev:    9.426598 63.50004
Fixed effects:  divedepth ~ oarea 
Value Std.Error     DF   t-value p-value
(Intercept) 25.549003  3.171766 225541  8.055135  0.0000
oarea2      12.619669  0.828729 225541 15.227734  0.0000
oarea3       1.095290  0.979873 225541  1.117787  0.2637
oarea4       0.852045  0.492100 225541  1.731447  0.0834
oarea5       2.441955  0.587300 225541  4.157933  0.0000
[snip]
Number of Observations: 225554
Number of Groups: 9 

但是，我找不到随机变量deployID的p值。如何查看此值？

Answer 1

如评论所述，GLMM常见问题解答中有关于随机效应显著性测试的内容。你绝对应该考虑：

• 为什么你真的对p值感兴趣(这不是从来没有感兴趣，但这是一种不寻常的情况(
• 事实上，似然比检验在检验方差参数时是非常保守的(在这种情况下，它给出的p值太大了2倍(

这里有一个例子表明，没有随机效应的lme()拟合和相应的lm()模型具有相称的对数似然性(即，它们以可比较的方式计算(，并且可以与anova():进行比较

加载包和模拟数据(零随机效应方差(

library(lme4)
library(nlme)
set.seed(101)
dd <- data.frame(x = rnorm(120), f = factor(rep(1:3, 40)))
dd$y <- simulate(~ x + (1|f),
newdata = dd,
newparams = list(beta = rep(1, 2),
theta = 0,
sigma = 1))[[1]]

拟合模型(请注意，您不能将拟合REML的模型与没有随机效应的模型进行比较(。

m1 <- lme(y ~ x , random = ~ 1 | f, data = dd, method = "ML")
m0 <- lm(y ~ x, data = dd)

测试：

anova(m1, m0)
##    Model df      AIC      BIC    logLik   Test      L.Ratio p-value
## m1     1  4 328.4261 339.5761 -160.2131                            
## m0     2  3 326.4261 334.7886 -160.2131 1 vs 2 6.622332e-08  0.9998

在这里，测试正确地确定了两个模型是相同的，并给出了1的p值。

如果您使用lme4::lmer而不是lme，那么您还有一些其他更准确(但速度较慢(的选项(用于基于模拟的测试的RLRsim和PBmodcomp包(：请参阅GLMM常见问题解答。