我正在处理具有高斯误差分布的多级模型,该模型具有约 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)