O otimizador lme4 padrão requer muitas iterações para dados de alta dimensão

12

TL; DR: a lme4otimização parece linear no número de parâmetros do modelo por padrão e é muito mais lenta que um glmmodelo equivalente com variáveis ​​dummy para grupos. Existe algo que eu possa fazer para acelerar isso?


Estou tentando ajustar um modelo de logit hierárquico bastante grande (~ 50k linhas, 100 colunas, 50 grupos). Ajustar um modelo de logit normal aos dados (com variáveis ​​fictícias para o grupo) funciona bem, mas o modelo hierárquico parece estar travando: a primeira fase de otimização é concluída com perfeição, mas a segunda passa por muitas iterações sem que nada mude e sem parar .

EDIT: Suspeito que o problema seja principalmente o fato de eu ter muitos parâmetros, porque quando tento definir maxfnum valor menor, ele emite um aviso:

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

No entanto, as estimativas de parâmetros não estão mudando no decorrer da otimização, então ainda estou confuso sobre o que fazer. Quando tentei definir maxfnos controles do otimizador (apesar do aviso), ele pareceu travar após concluir a otimização.

Aqui está um código que reproduz o problema para dados aleatórios:

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)

Isso gera:

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>

Tentei definir ncoloutros valores e parece que o número de iterações feitas é (aproximadamente) 40 por coluna. Obviamente, isso se torna uma dor enorme à medida que adiciono mais colunas. Existem ajustes que posso fazer no algoritmo de otimização que reduzirá a dependência do número de colunas?

Ben Kuhn
fonte
1
Seria útil conhecer o modelo específico que você está tentando ajustar (especialmente a estrutura de efeitos aleatórios).
Patrick S. Forscher
Infelizmente, o modelo preciso é proprietário. Há um nível de efeitos aleatórios, com tamanhos de grupo que variam entre ~ 100 e 5000. Deixe-me saber se posso fornecer outras informações relevantes sobre o modelo.
Ben Kuhn
OK, adicionei um código que reproduz o problema.
Ben Kuhn
1
Eu não tenho uma resposta completa para você, então deixarei isso como um comentário. Na minha experiência, glmeré bastante lento, especialmente para modelos que possuem uma estrutura complexa de efeitos aleatórios (por exemplo, muitas inclinações aleatórias, efeitos aleatórios cruzados etc.). Minha primeira sugestão seria tentar novamente com uma estrutura simplificada de efeitos aleatórios. No entanto, se você estiver enfrentando esse problema apenas com um modelo de interceptação aleatória, seu problema pode ser simplesmente o número de casos; nesse caso, você precisará tentar algumas ferramentas especializadas em big data.
Patrick S. Forscher
Ele tem o mesmo problema com 2 grupos em vez de 50. Além disso, testando com um número menor de colunas, parece que o número de iterações é aproximadamente linear no número de colunas ... Existem métodos de otimização que farão melhor aqui ?
Ben Kuhn

Respostas:

12

Uma coisa que você pode tentar é mudar o otimizador. Veja o comentário de Ben Bolker nesta edição do github . A implementação nlopt do bobyqa geralmente é muito mais rápida que o padrão (pelo menos sempre que eu tento).

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

Além disso, consulte esta resposta para obter mais opções e este tópico de R-sig-mixed-models (que parece mais relevante para o seu problema).

Editar: forneci algumas informações desatualizadas relacionadas a nloptr. Para lme4 1.1-7cima e para cima, nloptré importado automaticamente (consulte ?nloptwrap). Tudo o que você precisa fazer é adicionar

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

à sua ligação.

alexforrence
fonte
Obrigado! Estou tentando o código nlopt agora. Eu quero saber se há algo que não seja uma implementação otimizador ruim acontecendo, desde que se monte um glm dummified quase equivalente foi muito mais rápido, mas vou ver ...
Ben Kuhn
Bem, era certamente mais rápido, mas parou com um erro: PIRLS step-halvings failed to reduce deviance in pwrssUpdate. Você tem alguma ideia do que pode estar acontecendo aqui? A mensagem de erro não é exatamente transparente ...
Ben Kuhn
Para chutes, você pode tentar definir nAGQ = 0 (consulte o tópico que eu vinculei para mais algumas idéias). Não me lembro o que causa o erro PIRLS, mas vou dar uma olhada.
alexforrence
Muito obrigado! Você poderia me indicar um recurso em que eu pudesse aprender mais sobre os detalhes desses métodos para poder resolver problemas como esse no futuro? A otimização parece muito com magia negra para mim no momento.
Ben Kuhn
2
nAGQ = 0 funcionou para mim no seu exemplo de teste com o bobyqa padrão (executado em ~ 15 segundos) e em 11 segundos com o nloptrbobyqa. Aqui está uma entrevista com John C. Nash (co-autor dos pacotes optime optimx), onde ele faz uma explicação de alto nível sobre otimização. Se você olhar para cima optimxou nloptrem CRAN, seus respectivos manuais de referência vai lhe dizer mais sobre a sintaxe. nloptrtambém possui uma vinheta disponível, que detalha um pouco mais.
precisa saber é o seguinte