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

12

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?

Ben Kuhn
la source
1
Il serait utile de connaître le modèle spécifique que vous essayez d'adapter (en particulier la structure des effets aléatoires).
Patrick S.Forscher
Malheureusement, le modèle précis est propriétaire. Il existe un niveau d'effets aléatoires, avec des tailles de groupe comprises entre ~ 100 et 5000. Faites-moi savoir si je peux fournir d'autres informations pertinentes sur le modèle.
Ben Kuhn
OK, j'ai ajouté du code qui reproduit le problème.
Ben Kuhn
1
Je n'ai pas de réponse complète pour vous, alors je vais laisser cela comme un commentaire. D'après mon expérience, glmerest assez lent, en particulier pour les modèles qui ont une structure complexe d'effets aléatoires (par exemple, de nombreuses pentes aléatoires, des effets aléatoires croisés, etc.). Ma première suggestion serait d'essayer à nouveau avec une structure d'effets aléatoires simplifiée. Cependant, si vous rencontrez ce problème avec un modèle d'interceptions aléatoires uniquement, votre problème peut simplement être le nombre de cas, auquel cas vous devrez essayer des outils spécialisés pour les mégadonnées.
Patrick S.Forscher
Il a le même problème avec 2 groupes au lieu de 50. Aussi, en testant avec un plus petit nombre de colonnes, il semble que le nombre d'itérations soit à peu près linéaire dans le nombre de colonnes ... Existe-t-il des méthodes d'optimisation qui feront mieux ici ?
Ben Kuhn

Réponses:

12

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.

alexforrence
la source
Je vous remercie! J'essaie le code nlopt en ce moment. Je me demande s'il y a autre chose qu'une mauvaise implémentation d'optimiseur en cours, car l'ajustement d'un glm factice presque équivalent était tellement plus rapide, mais je vais voir ...
Ben Kuhn
Eh bien, il était certainement plus rapide, mais il est arrêté avec une erreur: PIRLS step-halvings failed to reduce deviance in pwrssUpdate. Avez-vous une idée de ce qui pourrait se passer ici? Le message d'erreur n'est pas exactement transparent ...
Ben Kuhn
Pour les coups de pied, vous pouvez essayer de définir nAGQ = 0 (voir le fil que j'ai lié pour quelques idées supplémentaires). Je ne me souviens pas des causes de l'erreur PIRLS, mais je vais regarder autour de moi.
alexforrence
Merci beaucoup! Pourriez-vous m'indiquer une ressource où je pourrais en savoir plus sur les détails de ces méthodes afin de pouvoir résoudre moi-même des problèmes comme celui-ci à l'avenir? L'optimisation me ressemble beaucoup à la magie noire en ce moment.
Ben Kuhn
2
nAGQ = 0 a fonctionné pour moi sur votre exemple de test avec le bobyqa par défaut (exécuté en ~ 15 sec), et en 11 sec avec le nloptrbobyqa. Voici une interview de John C. Nash (co-auteur des packages optimet optimx) où il fait une explication de haut niveau de l'optimisation. Si vous recherchez optimxou nloptrsur CRAN, leurs manuels de référence respectifs vous en diront plus sur la syntaxe. nloptra également une vignette disponible, qui va un peu plus loin dans les détails.
alexforrence