为什么 R 中 lmer 模型的参数的 Wald 和自举置信区间之间存在很大差异?

机器算法验证 r lme4-nlme
2022-04-17 09:59:58

我正在处理具有高斯误差分布的多级模型,该模型具有约 21,000 个观测值和 5000 个聚类。该模型具有简单的形式:

lmer(y ~ a + b + a:b + (1|z), weights=b)

并且具有与变量成比例的权重b

我遇到的行为发生在计算固定效应参数的置信区间时:

confint(mod, method="Wald")
confint(mod, method="profile") 
confint(mod1, method="boot", nsim=1000, parm="beta_")

自举的结果给出的置信区间比 Wald 结果宽约 3 倍。配置文件结果会引发许多警告,例如:

1:在 profile.merMod(object, which = parm, signames = oldNames, ...) 中:(拦截)的非单调配置文件

6:在 confint.thpr(pp, level = level, zeta = zeta) 中:样条拟合不良(截距):回退到线性插值

我搜索了许多比较这些方法的旧线程,我确实希望这些方法的结果会有所不同。然而,以前的帖子(以及我自己的经验)表明,这些方法通常会产生相差不大的结果。我正在寻找一些关于为什么会发生这种情况的直觉,以及我认为自举结果可能更现实是否正确(我想这意味着我也应该对报告的 SE 持怀疑态度summary()?)。

对于没有提供可重现的示例,我深表歉意,因为我无法共享原始数据,并且我尝试模拟相同问题只会导致 Wald 和引导式 CI 相似的情况。

编辑:简介图:

在此处输入图像描述

编辑#2:为了回应 Ben Bolker 在下面的评论,以下代码(从 Internet 上的其他地方窃取)模拟混合模型的数据,并演示 confint(..., method="profile") 在包含权重时失败以及来自不同方法的 CI 的一些差异。

library(mvtnorm)
set.seed(2345)

N <- 150
unit.df <- data.frame(unit = c(1:N), a = rnorm(N))
unit.df <-  within(unit.df, {
  E.alpha.given.a <-  1 - 0.15 * a
  E.beta.given.a <-  3 + 0.3 * a
})

q = 0.2
r = 0.9
s = 0.5
cov.matrix <- matrix(c(q^2, r * q * s, r * q * s, s^2), nrow = 2,
                     byrow = TRUE)
random.effects <- rmvnorm(N, mean = c(0, 0), sigma = cov.matrix)
unit.df$alpha <- unit.df$E.alpha.given.a + random.effects[, 1]
unit.df$beta <- unit.df$E.beta.given.a + random.effects[, 2]

J <- 300
M = J * N  #Total number of observations
x.grid = seq(-4, 4, by = 8/J)[0:30]
within.unit.df <-  data.frame(unit = sort(rep(c(1:N), J)), j = rep(c(1:J),N), x =rep(x.grid, N))
flat.df = merge(unit.df, within.unit.df)

flat.df <-  within(flat.df, y <-  alpha + x * beta + 0.75 * rnorm(n = M))
simple.df <-  flat.df[, c("unit", "a", "x", "y")]
simple.df$wht <- rpois(n = dim(simple.df)[1], lambda = 5)+1

my.lmer <-  lmer(y ~ x + a + x * a + wht + (1 | unit), data = simple.df, weights = wht)
summary(my.lmer)

confint(my.lmer, method="Wald")
confint(my.lmer, method="profile")
confint(my.lmer, method="boot", nsim=100)

1个回答

基本上,Wald 统计量并不好,对于混合模型,您不应该相信它。profile与使用andboot.ci方法相比,它使用了对实际可能性的更粗略的近似。如果 R(以及 SAS 和 JMP 和...)今天已经编写好了,他们就不会费心实施 Wald 统计了。这就是该summary.merMod方法有意从固定效应系数输出中省略按照今天的标准,profile/bootstrap 的计算强度最多只有几分钟,但在过去,这需要几周时间。因此,预计分析师将进行大量测试和变量转换方法,以便 Wald stat 可能具有良好的属性。p

编辑:下面是 2010 年我、大卫·达尔和道格拉斯·贝茨之间的对话片段,当时我试图建议.pxtable

您的 lme4 包的用户希望在 lme4 的 mer 对象上使用 xtable。这意味着定义一个函数“xtable.mer”。他建议以下实施。遗憾的是我对 lme4 不是很熟悉。你有什么建议吗?

我很欣赏亚当的建议和他提供的实施。遗憾的是,至少可以说,我认为实施会引起争议,而且我不希望成为后果的接受者。关于固定效应参数测试的 p 值,lme4 存在一个长期存在的问题。对于线性混合模型,人们普遍认为您可以计算 t 统计量(此处标记为“z 值”)并通过确定自由度的近似数量的简单权宜之计将其转换为 p 值. 事实上,SAS PROC MIXED 提供了几种(我相信是 6 种)不同且不兼容的方法来确定此类自由度和相应的 p 值。这些给出不同答案的事实并没有

实际上,这种统计量的分布不是学生 T。它比这要复杂得多,我提倡计算置信区间或检验假设的其他方法。在广义线性混合模型的情况下,我确实从标准正态分布计算 p 值,不是因为 GLMM 的近似值优于 LMM,而是因为它更差。

我正在为 lme4 上的 Springer 写一本书(章节草稿可在http://lme4.R-forge.R-project.org/book/获得),其中我描述了使用似然比检验进行假设检验和基于分析的技术用于生成参数置信区间的 LRT 统计量。那本书中的示例基于使用不同模型表示的包的开发版本。实现不完整,这就是为什么我没有将它作为 lme4 发布,但现在我需要集中精力写作,因为这本书将用于下周开始的研讨会。