El optimizador predeterminado de lme4 requiere muchas iteraciones para datos de alta dimensión

TL; DR: la lme4optimización parece ser lineal en la cantidad de parámetros del modelo por defecto, y es mucho más lenta que un glmmodelo equivalente con variables ficticias para grupos. ¿Hay algo que pueda hacer para acelerarlo?

Estoy tratando de ajustar un modelo logit jerárquico bastante grande (~ 50k filas, 100 columnas, 50 grupos). El ajuste de un modelo logit normal a los datos (con variables ficticias para el grupo) funciona bien, pero el modelo jerárquico parece estar estancado: la primera fase de optimización se completa bien, pero la segunda pasa por muchas iteraciones sin que nada cambie y sin detenerse .

EDITAR: sospecho que el problema es principalmente porque tengo tantos parámetros, porque cuando trato de establecer maxfnun valor más bajo, me da una advertencia:

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

Sin embargo, las estimaciones de los parámetros no cambian en absoluto a lo largo de la optimización, por lo que todavía estoy confundido sobre qué hacer. Cuando intenté configurar maxfnlos controles del optimizador (a pesar de la advertencia), pareció bloquearse después de finalizar la optimización.

Aquí hay un código que reproduce el problema de los datos aleatorios:

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)

Esto produce:

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>

Intenté establecer ncolotros valores, y parece que el número de iteraciones realizadas es (aproximadamente) 40 por columna. Obviamente, esto se convierte en un gran dolor a medida que agrego más columnas. ¿Hay ajustes que pueda hacer al algoritmo de optimización que reducirán la dependencia en el número de columnas?

Una cosa que puedes probar es cambiar el optimizador. Vea el comentario de Ben Bolker en este tema de github . La implementación nlopt de bobyqa suele ser mucho más rápida que la predeterminada (al menos cada vez que lo intento).

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"))

Además, consulte esta respuesta para obtener más opciones y este hilo de R-sig-mixed-models (que parece más relevante para su problema).

Editar: te di información desactualizada relacionada con nloptr. En lme4 1.1-7y hacia arriba, nloptrse importa automáticamente (ver ?nloptwrap). Todo lo que tienes que hacer es agregar

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

a su llamada