Domyślny optymalizator Lme4 wymaga dużej liczby iteracji dla danych wielowymiarowych

TL; DR: lme4optymalizacja wydaje się domyślnie liniowa pod względem liczby parametrów modelu i jest znacznie wolniejsza niż równoważny glmmodel ze zmiennymi fikcyjnymi dla grup. Czy mogę coś przyspieszyć?

Próbuję dopasować dość duży hierarchiczny model logit (~ 50 000 wierszy, 100 kolumn, 50 grup). Dopasowywanie normalnego modelu logu do danych (ze zmiennymi fikcyjnymi dla grupy) działa dobrze, ale wydaje się, że model hierarchiczny utknął: pierwsza faza optymalizacji kończy się dobrze, ale druga przechodzi wiele iteracji bez zmiany i bez zatrzymywania .

EDYCJA: Podejrzewam, że problem polega głównie na tym, że mam tak wiele parametrów, ponieważ gdy próbuję ustawić maxfnniższą wartość, pojawia się ostrzeżenie:

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

Jednak szacunki parametrów wcale się nie zmieniają w trakcie optymalizacji, więc nadal nie jestem pewien, co robić. Kiedy próbowałem ustawić maxfnelementy sterujące optymalizatora (pomimo ostrzeżenia), wydawało się, że zawiesiło się ono po zakończeniu optymalizacji.

Oto kod, który odtwarza problem losowych danych:

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)

To daje:

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>

Próbowałem ustawić ncolinne wartości i wydaje się, że liczba wykonanych iteracji wynosi (około) 40 na kolumnę. Oczywiście staje się to ogromnym bólem, gdy dodam więcej kolumn. Czy mogę ulepszyć algorytm optymalizacji, który zmniejszy zależność od liczby kolumn?

Jedną z rzeczy, których możesz spróbować, jest zmiana optymalizatora. Zobacz komentarz Bena Bolkera w tym numerze github . Implementacja nlopt bobyqa jest zwykle znacznie szybsza niż domyślna (przynajmniej za każdym razem, gdy próbuję).

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

Zobacz także tę odpowiedź, aby uzyskać więcej opcji i ten wątek z modeli R-sig-mixed-models (który wydaje się bardziej odpowiedni dla Twojego problemu).

Edycja: Dałem ci trochę nieaktualnych informacji związanych z nloptr. W lme4 1.1-7górę i w górę nloptrjest automatycznie importowany (patrz ?nloptwrap). Wszystko, co musisz zrobić, to dodać

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

na twój telefon.