Der standardmäßige lme4-Optimierer erfordert viele Iterationen für hochdimensionale Daten

12

TL; DR: Die lme4Optimierung scheint standardmäßig in Bezug auf die Anzahl der Modellparameter linear zu sein und ist viel langsamer als ein äquivalentes glmModell mit Dummy-Variablen für Gruppen. Kann ich irgendetwas tun, um es zu beschleunigen?


Ich versuche, ein ziemlich großes hierarchisches Logit-Modell (~ 50.000 Zeilen, 100 Spalten, 50 Gruppen) anzupassen. Das Anpassen eines normalen Logit-Modells an die Daten (mit Dummy-Variablen für die Gruppe) funktioniert einwandfrei, aber das hierarchische Modell scheint hängen zu bleiben: Die erste Optimierungsphase wird einwandfrei abgeschlossen, während die zweite viele Iterationen durchläuft, ohne dass sich etwas ändert und ohne anzuhalten .

EDIT: Ich vermute, dass das Problem hauptsächlich darin besteht, dass ich so viele Parameter habe, denn wenn ich versuche maxfn, einen niedrigeren Wert einzustellen , wird eine Warnung ausgegeben:

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

Die Parameterschätzungen ändern sich jedoch im Verlauf der Optimierung überhaupt nicht, sodass ich immer noch verwirrt bin, was zu tun ist. Als ich versuchte, maxfndie Steuerelemente des Optimierers einzustellen (trotz der Warnung), schien es nach Abschluss der Optimierung zu hängen.

Hier ist ein Code, der das Problem für zufällige Daten reproduziert:

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)

Dies gibt aus:

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>

Ich habe versucht ncol, andere Werte festzulegen, und es scheint, dass die Anzahl der durchgeführten Iterationen (ungefähr) 40 pro Spalte beträgt. Offensichtlich wird dies ein großer Schmerz, wenn ich mehr Spalten hinzufüge. Gibt es Optimierungen am Optimierungsalgorithmus, die die Abhängigkeit von der Anzahl der Spalten verringern?

Ben Kuhn
quelle
1
Es wäre hilfreich, das spezifische Modell zu kennen, das Sie anpassen möchten (insbesondere die Zufallseffektstruktur).
Patrick S. Forscher
Leider ist das genaue Modell proprietär. Es gibt eine Ebene von zufälligen Effekten mit Gruppengrößen zwischen ~ 100 und 5000. Lassen Sie mich wissen, ob ich weitere relevante Informationen zum Modell bereitstellen kann.
Ben Kuhn
OK, ich habe Code hinzugefügt, der das Problem reproduziert.
Ben Kuhn
1
Ich habe keine vollständige Antwort für Sie, daher werde ich dies als Kommentar hinterlassen. Nach meiner Erfahrung glmerist es ziemlich langsam, insbesondere für Modelle, die eine komplexe zufällige Effektstruktur aufweisen (z. B. viele zufällige Steigungen, gekreuzte zufällige Effekte usw.). Mein erster Vorschlag wäre, es noch einmal mit einer vereinfachten Zufallseffektstruktur zu versuchen. Wenn dieses Problem jedoch nur bei einem zufälligen Abfangmodell auftritt, liegt das Problem möglicherweise in der Anzahl der Fälle. In diesem Fall müssen Sie einige Tools ausprobieren, die auf Big Data spezialisiert sind.
Patrick S. Forscher
Es hat das gleiche Problem mit 2 statt 50 Gruppen. Außerdem scheint es, als ob die Anzahl der Iterationen bei einer geringeren Anzahl von Spalten in etwa linear ist. Gibt es Optimierungsmethoden, die hier besser funktionieren? ?
Ben Kuhn,

Antworten:

12

Sie können versuchen, das Optimierungsprogramm zu ändern. Siehe Ben Bolkers Kommentar zu dieser Github-Ausgabe . Die nlopt-Implementierung von bobyqa ist normalerweise viel schneller als die Standardimplementierung (zumindest, wenn ich es versuche).

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

In dieser Antwort finden Sie auch weitere Optionen und diesen Thread von R-sig-mixed-models (der für Ihr Problem relevanter aussieht).

Bearbeiten: Ich gab Ihnen einige veraltete Informationen im Zusammenhang mit nloptr. In lme4 1.1-7und up nloptrwird automatisch importiert (siehe ?nloptwrap). Alles was Sie tun müssen, ist hinzuzufügen

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

zu Ihrem Anruf.

alexforrence
quelle
Vielen Dank! Ich versuche gerade den nlopt-Code. Ich frage mich, ob etwas anderes als eine schlechte Optimierer-Implementierung im Gange ist, da die Anpassung eines fast äquivalenten dummifizierten glm so viel schneller war, aber ich werde sehen ...
Ben Kuhn
Nun, es war auf jeden Fall schneller, aber es hielt mit einem Fehler: PIRLS step-halvings failed to reduce deviance in pwrssUpdate. Haben Sie eine Idee, was hier los sein könnte? Die Fehlermeldung ist nicht gerade transparent ...
Ben Kuhn
Für Kicks können Sie versuchen, nAGQ = 0 zu setzen (siehe den von mir verlinkten Thread für ein paar weitere Ideen). Ich erinnere mich nicht, was den PIRLS-Fehler verursacht, aber ich werde mich umsehen.
alexforrence
Vielen Dank! Könnten Sie mich auf eine Ressource verweisen, in der ich mehr über die Details dieser Methoden erfahren kann, damit ich in Zukunft Probleme wie diese selbst lösen kann? Die Optimierung fühlt sich im Moment sehr nach schwarzer Magie an.
Ben Kuhn
2
nAGQ = 0 arbeitete für mich an Ihrem Testbeispiel mit dem Standard-Bobyqa (lief in ~ 15 Sekunden) und in 11 Sekunden mit dem nloptrBobyqa. Hier ist ein Interview mit John C. Nash (Co-Autor der optimund optimx-Pakete), in dem er eine allgemeine Erklärung zur Optimierung gibt. Wenn Sie nachschlagen optimxoder nloptrCRAN verwenden, erfahren Sie in den entsprechenden Referenzhandbüchern mehr über die Syntax. nloptrEs gibt auch eine Vignette, die etwas detaillierter ist.
alexforrence