默认 lme4 优化器需要对高维数据进行大量迭代

机器算法验证 r 混合模式 优化 lme4-nlme 数字
2022-03-25 04:47:19

TL；DR：lme4默认情况下，优化在模型参数的数量上似乎是线性的，并且比具有组虚拟变量的等效模型要慢得多。glm我有什么办法可以加快速度吗？

我正在尝试拟合一个相当大的分层 logit 模型（约 50k 行，100 列，50 个组）。将正常的 logit 模型拟合到数据（使用组的虚拟变量）工作正常，但分层模型似乎卡住了：第一个优化阶段完成得很好，但第二个优化阶段经历了很多迭代，没有任何改变，也没有停止.

编辑：我怀疑问题主要是我有这么多参数，因为当我尝试设置maxfn为较低的值时，它会发出警告：

Warning message:
In commonArgs(par, fn, control, environment()) :
  maxfun < 10 * length(par)^2 is not recommended.

但是，参数估计在优化过程中根本没有改变，所以我仍然对该怎么做感到困惑。当我尝试maxfn在优化器控件中设置时（尽管有警告），它似乎在完成优化后挂起。

这是一些重现随机数据问题的代码：

library(lme4)

set.seed(1)

SIZE <- 50000
NGRP <- 50
NCOL <- 100

test.case <- data.frame(i=1:SIZE)
test.case[["grouping"]] <- sample(NGRP, size=SIZE, replace=TRUE, prob=1/(1:NGRP))
test.case[["y"]] <- sample(c(0, 1), size=SIZE, replace=TRUE, prob=c(0.05, 0.95))

test.formula = y ~ (1 | grouping)

for (i in 1:NCOL) {
    colname <- paste("col", i, sep="")
    test.case[[colname]] <- runif(SIZE)
    test.formula <- update.formula(test.formula, as.formula(paste(". ~ . +", colname)))
}

print(test.formula)

test.model <- glmer(test.formula, data=test.case, family='binomial', verbose=TRUE)

这输出：

start par. =  1 fn =  19900.78 
At return
eval:  15 fn:      19769.402 par:  0.00000
(NM) 20: f = 19769.4 at           0     <other numbers>
(NM) 40: f = 19769.4 at           0     <other numbers>

我尝试设置ncol为其他值，似乎完成的迭代次数是（大约）每列 40 次。显然，当我添加更多列时，这变得非常痛苦。我可以对优化算法进行调整以减少对列数的依赖吗？

1个回答

您可以尝试的一件事是更改优化器。请参阅 Ben Bolker 在此 github 问题上的评论。bobyqa 的 nlopt 实现通常比默认的快得多（至少每当我尝试它时）。

library(nloptr)
defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1e-6,maxeval=1e5)
nloptwrap2 <- function(fn,par,lower,upper,control=list(),...) {
    for (n in names(defaultControl)) 
      if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
    res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...)
    with(res,list(par=solution,
                  fval=objective,
                  feval=iterations,
                  conv=if (status>0) 0 else status,
                  message=message))
}

system.time(test.model <- glmer(test.formula, data=test.case, 
family='binomial', verbose=TRUE))

system.time(test.model2 <- update(test.model,
control=glmerControl(optimizer="nloptwrap2"))

另外，请参阅此答案以获取更多选项以及来自 R-sig-mixed-models 的此线程（看起来与您的问题更相关）。

编辑： 我给了你一些与nloptr. Inlme4 1.1-7及以上，nloptr是自动导入的（请参阅 参考资料?nloptwrap）。您所要做的就是添加

control = [g]lmerControl(optimizer = "nloptwrap") # +g if fitting with glmer

接您的电话。

其它你可能感兴趣的问题
上一篇Python中的QQ图 下一篇为排名数据（Spearman 相关性）绘制回归线是否“可以”？