如何获得混合模型(lme4)中分类因子的“整体”p值和效应大小?

机器算法验证 r 假设检验 方差分析 混合模式 lme4-nlme
2022-02-09 23:59:29

我想获得一个独立分类变量(具有多个级别)的 p 值和效果大小 - 这是“整体”而不是单独针对每个级别,就像lme4R 中的正常输出一样。就像人们在运行 ANOVA 时报告的内容。

我怎样才能得到这个?

3个回答

您提到的两个概念(线性混合模型的 p 值和效应大小)都存在固有问题。关于效果大小,引用 Doug Bates,原作者lme4

假设一个人想要定义一个度量,我认为可以提出一个论点来处理来自线性混合模型的惩罚残差平方和,就像我们考虑来自线性模型的残差平方和一样。或者可以只使用没有惩罚的残差平方和,或者可以从给定的一组项中获得的最小残差平方和,这对应于一个无限精度矩阵。我不知道,真的。这取决于您要描述的内容。R2

有关更多信息,您可以查看此线程此线程此消息基本上,问题在于没有一致同意的方法来包含和分解模型中随机效应的方差。但是,使用了一些标准。如果您查看为/由 r-sig-mixed-models 邮件列表设置的 Wiki,则列出了几种方法。

建议的方法之一着眼于拟合值和观察值之间的相关性。这可以按照Jarrett Byrnes 在其中一个线程中的建议在R中实现:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

例如,假设我们估计以下线性混合模型:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

我们可以使用上面定义的函数计算效果大小:

r2.corr.mer(fm1)
# [1] 0.0160841

Ronghui Xu的一篇论文中推荐了一个类似的替代方案,称为,并且可以在R中简单地计算:Ω02

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

关于 p 值,这是一个更具争议性的问题(至少在R /lme4社区中)。请参阅此处此处此处问题中的讨论。再次参考 Wiki 页面,有几种方法可以测试线性混合模型中效果的假设。从“最差到最好”列出(根据我认为包括 Doug Bates 和 Ben Bolker 在内的Wiki 页面的作者,他们在这里做出了很多贡献):

  • Wald Z 检验
  • 对于可以计算 df 的平衡嵌套 LMM:Wald t-tests
  • 似然比检验,通过设置模型以便可以隔离/删除参数(通过anovadrop1),或通过计算似然分布
  • MCMC 或参数引导置信区间

他们推荐马尔可夫链蒙特卡洛抽样方法,并列出了从伪贝叶斯方法和完全贝叶斯方法实现这一点的多种可能性,如下所列。

伪贝叶斯:

  • 事后抽样,通常(1)假设平坦的先验和(2)从 MLE 开始,可能使用近似方差-协方差估计来选择候选分布
  • 通过mcmcsamp(如果适用于您的问题:即具有简单随机效应的 LMM — 不是 GLMM 或复杂随机效应)
    通过pvals.fnclanguageR包中,一个包装器mcmcsamp
  • 在 AD 模型生成器中,可能通过glmmADMB包(使用mcmc=TRUE选项)或R2admb包(在 AD 模型生成器中编写您自己的模型定义),或在R之外
  • 通过包中的sim函数arm(仅模拟 beta(固定效应)系数的后验

完全贝叶斯方法:

  • 通过MCMCglmm
  • 使用glmmBUGS(WinBUGS 包装器/ R接口)
  • 使用 JAGS/WinBUGS/OpenBUGS 等通过rjags///r2jagsR2WinBUGSBRugs

为了说明这可能看起来像什么,下面是MCMCglmm使用包的估计值MCMCglmm,您将看到它产生与上述模型相似的结果,并且具有某种贝叶斯 p 值:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

我希望这会有所帮助。我认为对于从线性混合模型开始并尝试在R中估计它们的人的最佳建议是阅读Wiki常见问题解答,其中大部分信息都是从中提取的。它是各种混合效果主题的绝佳资源,从基础到高级,从建模到绘图。

关于计算显着性 ( p ) 值,Luke (2016) Evaluating显着性在 R 中的线性混合效应模型报告说,最佳方法是自由度的 Kenward-Roger 或 Satterthwaite 近似(在 R 中可用,包如lmerTestafex)。

抽象的

混合效应模型在实验数据分析中的使用越来越频繁。然而,在 R 中的 lme4 包中,用于评估这些模型中固定效应显着性的标准(即,获得 p 值)有些模糊。这有充分的理由,但由于在许多情况下需要使用这些模型的研究人员报告 p 值,因此需要一些方法来评估模型输出的重要性。本文报告的模拟结果表明,评估显着性的两种最常用方法,使用似然比检验和将 z 分布应用于模型输出 (t-as-z) 中的 Wald t 值,在某种程度上是反保守的,特别是对于较小的样本量。其他评估显着性的方法,这些模拟的结果表明,当使用 REML 拟合模型并使用 Kenward-Roger 或 Satterthwaite 近似值导出 p 值时,类型 1 错误率最接近于 0.05,因为这些近似值都产生了可接受的类型 1 错误率,即使对于较小的样品。

(重点补充)

我用这个lmerTest包。这很方便地包括对我的 MLM 分析输出中的 p 值的估计anova(),但由于此处其他帖子中给出的原因,它没有给出效果大小。