大型数据集的 LMER 对象预测均值的置信区间

如何解决大型数据集的 LMER 对象预测均值的置信区间

我想在大型数据集（约 40k 行）上获得线性混合效应模型的预测均值的置信区间 (CI)，该数据集本身是更大数据集的子集。然后使用该 CI 来估计另一个使用平均值及其相关 CI 作为输入数据的计算的不确定性。

我设法为完整数据集创建了预测估计和区间，但预测区间与 CI 不同且大得多。除了引导（处理这么多数据需要太多时间）之外，我找不到一种方法可以让我估计 CI - 要么是因为它会抛出错误，要么是因为它只提供计算预测区间。

我最近才搬到 LME，因此我可能已经监督了一些明显的方法。

以下是我目前所做的更详细的工作：

输入数据是保密的，因此很遗憾我不能分享任何摘录。

但一般来说，我们有一个表示事件概率的因变量 (y) 和 2 个分类变量（c1 和 c2）和两个带有一些加权因子 (w1) 的连续变量（x1 和 x2）。数据集中的某些值丢失。数据第一行的摘录可能类似于以下示例：

c1	c2	x1	x2	w1	y
伦敦	小	1	10	不适用	不适用
伦敦	小	1	20	不适用	不适用
伦敦	大	2	10	0.2	0.1
巴黎	小	1	10	0.2	0.23
巴黎	大	2	10	0.3	0.3

根据这些输入数据，我将按照以下形式拟合 LMER 模型：

lmer1 <- lme4::lmer( y ~ x1 * poly(x2,5) + ((x1 * poly(x2,5)) | c1),data = df,weights = w1,control = lme4::lmerControl(check.conv.singular = lme4::.makeCC(action = "ignore",tol = 1e-3)))

这会运行几分钟并返回几个警告：

警告信息：1：在 optwrap(optimizer,devfun,getStart(start,rho$pp),lower = rho$lower,: nloptwrap 的收敛代码 5: NLOPT_MAXEVAL_REACHED：优化停止，因为 maxeval（上）已达成。

2: 在 checkConv(attr(opt,“derivs”),opt$par,ctrl = control$checkConv,: 无法评估缩放梯度

3: 在 checkConv(attr(opt,: 模型未能收敛：退化 Hessian 11个负特征值

我增加了 MAXEVAL 参数，但这仍然无助于摆脱警告，我发现尽管有这些警告，模型仍然适合。因此，我开始应用不同的方法来预测整个数据集的均值和均值的相关 CI。

预测间隔

我开始为完整数据集创建预测区间：

predictions <- merTools::predictInterval(lmer1,newdata = df,which = "full",n.sims = 1000,include.resid.var = FALSE,level=0.95,stat="mean")

但是，如上所述，预测区间与 CI 不同（另见 https://datascienceplus.com/prediction-interval-the-wider-sister-of-confidence-interval/）。

我发现一般的预测函数可以选择将区间设置为“预测”或“置信度”，但这个选项不存在于来自 LMER 对象的预测。而且我找不到另一种从预测间隔切换到 CI 的可能性——尽管我相信绘制的数据应该足以做到这一点。

约束

然后我看到有一个叫做“confint”的函数，但是在运行这个函数时我得到以下错误：

predicition_ci = lme4::confint.merMod(lmer1)

正在计算配置文件置信区间...

zeta(shiftpar,start = opt[seqpar1][-w]) 中的错误：分析检测到新的、较低的偏差

此外：警告消息：

1: 在 commonArgs(par,fn,control,environment()) 中：maxfun

2: 在 optwrap(optimizer,x@theta,lower = x@lower,calc.derivs = TRUE,: bobyqa 的收敛代码 1: bobyqa -- 超过最大函数评估次数

我找到了这个线程 (Error when estimating CI for GLMM using confint())，它说我需要通过设置不同的配置文件来减少“devtol”参数。但是这样做会导致相同的错误：

lmer1_devtol = profile(lmer1,devtol = 1e-7)

zeta(shiftpar,: bobyqa 的收敛代码 1: bobyqa -- 超过最大函数评估次数

add_ci

我找到了函数“add_ci”，但这又导致了另一个错误：

predictions_ci = ciTools::add_ci(df,lmer1,alpha = 0.05)

levelfun(r,n,allow.new.levels = allow.new.levels) 中的错误：新在新数据中检测到的水平

然后我像预测函数的描述一样将新的“allow.new.levels”参数设置为TRUE，但这个参数似乎没有被执行：

predictions_ci = ciTools::add_ci(df,alpha = 0.05,allow.new.levels = TRUE)

levelfun(r,allow.new.levels = allow.new.levels) 中的错误：新在新数据中检测到的水平

诊断

我找到了一种计算 sleepstudy 数据的 CI 间隔的方法，该方法使用带有 diag 的矩阵转换。

Designmat <- model.matrix(as.formula("y ~ x1 * poly(x2,5)")[-2],df)
   predvar <- diag(Designmat %*% vcov(lmer1) %*% t(Designmat)) 
    
#With new data
newdat = df
newdat$pred <- predict(lmer1,newdat,allow.new.levels = TRUE)
    
 Designmat <- model.matrix(formula(lmer1)[-2],newdat)

但是 diag 方法不适用于如此大的数据集。

bootMer

如前所述，使用 bootMer 增强置信区间对这个数据子集花费了太多时间（我在 1 天前开始它并且它仍在运行）。我尝试对 sleepstudy 样本数据进行一些并行处理，但这无法显着提高速度，因此我认为它对我的大型数据集也有同样的效果。

merBoot <- bootMer(lmer1,predict,nsim = 1000,re.form = NA)

其他

我已经通读了所有这些帖子（以及更多），但没有一个可以帮助我在合理的时间内为我的案例获得 CI。但也许我监督了一些事情。

大型数据集的 LMER 对象预测均值的置信区间

如何解决大型数据集的 LMER 对象预测均值的置信区间

相关推荐