TL은, DR은 : lme4최적화 기본적 모델 파라미터의 수에 선형으로 나타날 것이다 방법 등가보다 느리게 glm그룹 더미 변수 모델. 속도를 높이기 위해 할 수있는 일이 있습니까?

상당히 큰 계층 적 로짓 모델 (~ 50k 행, 100 열, 50 그룹)을 맞추려고합니다. 일반 로짓 모델을 데이터에 그룹화하면 (그룹에 더미 변수가 있음) 제대로 작동하지만 계층 적 모델이 멈춘 것처럼 보입니다. 첫 번째 최적화 단계는 정상적으로 완료되지만 두 번째는 변경하지 않고 멈추지 않고 많은 반복을 거칩니다. .

편집 : 문제는 주로 너무 많은 매개 변수가 있다고 생각합니다 maxfn. 낮은 값 으로 설정하려고 하면 경고가 표시 되기 때문입니다 .

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

그러나 매개 변수 추정치는 최적화 과정에서 전혀 변경되지 않으므로 수행해야 할 작업에 대해 여전히 혼란 스럽습니다. maxfn옵티 마이저 컨트롤에서 설정 을 시도했을 때 (경고에도 불구하고) 최적화가 완료된 후 멈추는 것처럼 보였습니다.

무작위 데이터에 대한 문제를 재현하는 코드는 다음과 같습니다.

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)

이 결과는 다음과 같습니다.

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>

ncol다른 값으로 설정 을 시도했는데 반복 횟수가 열 당 40 정도 인 것으로 보입니다. 분명히 더 많은 열을 추가함에 따라 이것은 큰 고통이됩니다. 열 수에 대한 의존성을 줄이는 최적화 알고리즘을 조정할 수 있습니까?

답변

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

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