如何在线性混合效应模型中找到随机效应的 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

但是，我找不到随机变量的 p 值：deployID。我怎样才能看到这个值？

Answer 1

如评论中所述，GLMM FAQ 中有关于随机效应显着性检验的内容。你绝对应该考虑：

为什么您真的对 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 常见问题解答。

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

How do I find the p-value for my random effect in my linear mixed effect model?

r

lme4

p-value