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

12

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?

Ben Kuhn
źródło
1
Przydałaby się znajomość konkretnego modelu, który próbujesz dopasować (zwłaszcza struktura efektów losowych).
Patrick S. Forscher,
Niestety dokładny model jest zastrzeżony. Istnieje jeden poziom efektów losowych, z wielkościami grup od ~ 100 do 5000. Daj mi znać, czy mogę podać inne istotne informacje na temat modelu.
Ben Kuhn,
OK, dodałem kod, który odtwarza problem.
Ben Kuhn,
1
Nie mam dla ciebie pełnej odpowiedzi, więc zostawię to jako komentarz. Z mojego doświadczenia glmerwynika , że jest dość wolny, szczególnie w przypadku modeli, które mają złożoną strukturę efektów losowych (np. Wiele losowych nachyleń, skrzyżowane efekty losowe itp.). Moją pierwszą sugestią byłoby spróbować ponownie z uproszczoną strukturą efektów losowych. Jeśli jednak występuje ten problem tylko z przypadkowym modelem przechwytywania, problemem może być po prostu liczba przypadków, w którym to przypadku należy wypróbować narzędzia wyspecjalizowane dla dużych zbiorów danych.
Patrick S. Forscher,
Ma ten sam problem z 2 grupami zamiast 50. Ponadto, testując z mniejszą liczbą kolumn, wydaje się, że liczba iteracji jest mniej więcej liniowa w liczbie kolumn ... Czy istnieją metody optymalizacji, które lepiej tutaj ?
Ben Kuhn

Odpowiedzi:

12

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.

alexforrence
źródło
Dziękuję Ci! Próbuję teraz kodu nlopt. Zastanawiam się, czy dzieje się coś innego niż zła implementacja optymalizatora, ponieważ dopasowanie prawie równoważnego zmumifikowanego glm było o wiele szybsze, ale zobaczę ...
Ben Kuhn
Dobrze, że to na pewno szybciej, ale zatrzymał się z błędem: PIRLS step-halvings failed to reduce deviance in pwrssUpdate. Czy masz pojęcie, co się tutaj dzieje? Komunikat o błędzie nie jest dokładnie przezroczysty ...
Ben Kuhn
W przypadku kopnięć możesz spróbować ustawić nAGQ = 0 (zobacz wątek, który podłączyłem, aby uzyskać więcej pomysłów). Nie pamiętam, co powoduje błąd PIRLS, ale rozejrzę się.
alexforrence
Dzięki wielkie! Czy możesz wskazać mi zasób, w którym mógłbym dowiedzieć się więcej o szczegółach tych metod, abym mógł samodzielnie rozwiązać takie problemy w przyszłości? Optymalizacja wydaje mi się w tej chwili bardzo podobna do czarnej magii.
Ben Kuhn
2
nAGQ = 0 działało dla mnie na twoim przykładzie testowym z domyślnym bobyqa (działało w ~ 15 sekund), a w 11 sekund z nloptrbobyqa. Oto wywiad z Johnem C. Nashem (współautorem pakietów optimi optimx), w którym wyjaśnia on na wysokim poziomie optymalizację. Jeśli spojrzysz w górę optimxlub nloptrna CRAN, ich odpowiednie podręczniki zawierają więcej informacji na temat składni. nloptrdostępna jest również winieta, która jest nieco bardziej szczegółowa.
alexforrence