L'optimiseur lme4 par défaut nécessite de nombreuses itérations pour les données de grande dimension

TL; DR: l' lme4optimisation semble être linéaire dans le nombre de paramètres du modèle par défaut, et est beaucoup plus lente qu'un glmmodèle équivalent avec des variables factices pour les groupes. Puis-je faire quelque chose pour l'accélérer?

J'essaie d'adapter un modèle logit hiérarchique assez grand (~ 50k lignes, 100 colonnes, 50 groupes). L'ajustement d'un modèle logit normal aux données (avec des variables factices pour le groupe) fonctionne bien, mais le modèle hiérarchique semble se coincer: la première phase d'optimisation se termine bien, mais la seconde passe par beaucoup d'itérations sans rien changer et sans s'arrêter .

EDIT: Je soupçonne que le problème est principalement dû au fait que j'ai beaucoup de paramètres, car lorsque j'essaie de définir maxfnune valeur inférieure, cela donne un avertissement:

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

Cependant, les estimations des paramètres ne changent pas du tout au cours de l'optimisation, donc je ne sais toujours pas quoi faire. Lorsque j'ai essayé de définir maxfnles commandes de l'optimiseur (malgré l'avertissement), il semblait se bloquer après avoir terminé l'optimisation.

Voici un code qui reproduit le problème pour les données aléatoires:

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)

Cela produit:

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>

J'ai essayé de définir ncold'autres valeurs, et il semble que le nombre d'itérations effectuées soit (environ) 40 par colonne. Évidemment, cela devient une énorme douleur lorsque j'ajoute plus de colonnes. Y a-t-il des ajustements à apporter à l'algorithme d'optimisation qui réduiront la dépendance au nombre de colonnes?

Vous pouvez essayer de changer l'optimiseur. Voir le commentaire de Ben Bolker dans ce numéro de github . L'implémentation nlopt de bobyqa est généralement beaucoup plus rapide que la valeur par défaut (au moins chaque fois que je l'essaye).

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

Consultez également cette réponse pour plus d'options et ce fil de discussion des modèles mixtes R-sig (qui semble plus pertinent pour votre problème).

Edit: Je vous ai donné des informations obsolètes liées à nloptr. Dans lme4 1.1-7et vers le haut, nloptrest automatiquement importé (voir ?nloptwrap). Tout ce que vous avez à faire est d'ajouter

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

à votre appel.