Pengoptimal lme4 default membutuhkan banyak iterasi untuk data dimensi tinggi

TL; DR: lme4optimasi tampaknya linier dalam jumlah parameter model secara default, dan jauh lebih lambat daripada glmmodel yang setara dengan variabel dummy untuk grup. Apakah ada yang bisa saya lakukan untuk mempercepatnya?

Saya mencoba menyesuaikan model logit hierarkis yang cukup besar (~ baris 50k, 100 kolom, 50 grup). Memasukkan model logit normal ke data (dengan variabel dummy untuk grup) berfungsi dengan baik, tetapi model hierarkis tampaknya macet: fase optimisasi pertama selesai dengan baik, tetapi yang kedua melewati banyak iterasi tanpa perubahan apa pun dan tanpa berhenti .

EDIT: Saya menduga masalahnya terutama karena saya memiliki begitu banyak parameter, karena ketika saya mencoba mengatur maxfnke nilai yang lebih rendah itu memberikan peringatan:

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

Namun, perkiraan parameter tidak berubah sama sekali selama optimasi, jadi saya masih bingung tentang apa yang harus dilakukan. Ketika saya mencoba mengatur maxfnkontrol pengoptimal (meskipun ada peringatan), sepertinya hang setelah menyelesaikan optimasi.

Berikut beberapa kode yang mereproduksi masalah untuk data acak:

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)

Output ini:

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>

Saya mencoba pengaturan ncolke nilai-nilai lain, dan tampaknya jumlah iterasi yang dilakukan adalah (kurang-lebih) 40 per kolom. Jelas, ini menjadi sangat menyakitkan ketika saya menambahkan lebih banyak kolom. Apakah ada tweak yang bisa saya buat untuk algoritma optimasi yang akan mengurangi ketergantungan pada jumlah kolom?

Satu hal yang bisa Anda coba adalah mengubah pengoptimal. Lihat komentar Ben Bolker di masalah github ini . Implementasi nlopt dari bobyqa biasanya jauh lebih cepat daripada default (setidaknya setiap kali saya mencobanya).

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

Juga, lihat jawaban ini untuk lebih banyak opsi dan utas ini dari R-sig-mixed-models (yang terlihat lebih relevan dengan masalah Anda).

Sunting: Saya memberi Anda beberapa informasi kedaluwarsa terkait nloptr. Masuk lme4 1.1-7dan naik, nloptrsecara otomatis diimpor (lihat ?nloptwrap). Yang harus Anda lakukan adalah menambahkan

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

untuk panggilan Anda.