如何在 R 中的多级模型中测试随机效应

机器算法验证 r 假设检验 混合模式 多层次分析 面板数据
2022-03-13 05:35:23

我一直在阅读 Judith Singer 和 John Willet 的一本好书，名为Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence。该书表明，通过在 2 个层次上建模，我们可以在 1 级和 2 级模型中对个体变化进行建模，以实现系统性的个体间变化差异。

示例的R 代码仅显示如何使用lme()来估计固定效应和随机效应。但是，文本建议我们应该测试方差分量以确定随机效应是否显着。

例如，其中一个代码仅执行以下操作：

library(nlme)

model.a <- lme(alcuse~ 1, alcohol1, random= ~1 |id)
summary(model.a)

Linear mixed-effects model fit by REML
 Data: alcohol1 
       AIC      BIC    logLik
  679.0049 689.5087 -336.5025

Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev:   0.7570578 0.7494974

Fixed effects: alcuse ~ 1 
                Value  Std.Error  DF  t-value p-value
(Intercept) 0.9219549 0.09629638 164 9.574139       0

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8892070 -0.3079143 -0.3029178  0.6110925  2.8562135 

Number of Observations: 246
Number of Groups: 82

但是文本列出了以下内容：

  • 固定效果：0.922*** (s = 0.096) -> 在输出中可用
  • 人内方差：0.562*** (s = 0.062) -> 可以从输出中获得（随机效应残差标准偏差平方）
  • 人与人之间的差异：0.564*** (s = 0.119)

我的工作涉及大量纵向数据分析，所以我真的需要理解这个想法。非常感激您的帮忙。

4个回答

第一点：

如果要测试随机效应的方差是否为 0，则需要小心：用于 LR渐近线不适用于 LRT 或方差分量的受限似然比检验，因为零假设位于参数空间的边缘。使用这种错误的参考分布将产生荒谬的保守测试。χ2anova()

资料来源：Self & Liang (1987)、Crainiceanu & Ruppert (2003)、Greven、Crainiceanu、Küchenhoff (2008)。

对于具有不相关随机效应的线性混合模型中的方差的精确 LR 测试，您可以使用我的包RLRsim。对于广义线性混合模型或具有相关随机效应的模型，如果 LRT 或 RLRT（后者具有更大的功率，请参阅上面的来源）位于临界值附近，我强烈建议使用参数引导来近似正确的 p 值（错误的）标准参考分布的值。代码lme4在帖子末尾。

第二点：

使用标准置信区间进行方差参数估计可能很危险，因为它们的分布通常非常偏斜，因此使用对称区间有点简单化——Douglas Bates 在这里有一些材料。

lme4-models的参数引导代码：

#m0 is the lmer model under the null hypothesis (i.e. the smaller model)
#mA is the lmer model under the alternative
bootstrapAnova <- function(mA, m0, B=1000){
     oneBootstrap <- function(m0, mA){
         d <- drop(simulate(m0))
         m2 <-refit(mA, newresp=d)
         m1 <-refit(m0, newresp=d)
         return(anova(m2,m1)$Chisq[2])
     }  
     nulldist <-  if(!require(multicore)){
         replicate(B, oneBootstrap(m0, mA))
     } else {
         unlist(mclapply(1:B, function(x) oneBootstrap(m0, mA)))
     }   
     ret <- anova(mA, m0)
     ret$"Pr(>Chisq)"[2] <- mean(ret$Chisq[2] < nulldist)
     names(ret)[7] <- "Pr_boot(>Chisq)"
     attr(ret, "heading") <- c(attr(ret, "heading")[1], 
          paste("Parametric bootstrap with", B,"samples."),
          attr(ret, "heading")[-1])
     attr(ret, "nulldist") <- nulldist
     return(ret)
}
#use like this (and increase B if you want reviewers to believe you):
(bRLRT <- bootstrapAnova(mA=<BIG MODEL>, m0=<SMALLER MODEL>))

（我会将其发布为对“chi”先前答案的评论，但在这里看不到这是可能的。）

对包装要非常小心intervals()，它们具有@fabians 上面指出的缺陷——它们依赖于标准卡方分布，适用于方差似然曲线形状的二次近似。较新的包（在 r-forge 上）允许为方差参数创建似然分布，它负责二次近似部分（尽管不是分布部分）。anova()nlmelme4a

此外，关于丢弃非显着方差分量是否是一个好主意的争论相当激烈（没有立即参考，但这已在 r-sig-mixed-models 上进行了讨论）。

该intervals()函数应该为您置信区间，有关更多信息，请参阅。您还可以测试是否可以通过使用从模型中删除任何方差分量（这相当于在两个嵌套模型之间进行 LRT）。100(1−α)help(intervals.lme)anova()

R 中多级分析的最佳资源之一是 John Fox 的“An R and S-PLUS Companion to Applied Regression”文本的网络附录。它提供了从 R NLME 输出计算一些更熟悉的度量的方法和方法的一个很好的概述。该附录可在 CRAN 上获得。

链接在这里：

http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-mixed-models.pdf

其它你可能感兴趣的问题
上一篇声望魔术师悖论 下一篇神经网络可以处理负输入和零输入吗？