12
$\begingroup$

TL;DR: lme4 optimization appears to be linear in the number of model parameters by default, and is way slower than an equivalent glm model with dummy variables for groups. Is there anything I can do to speed it up?


I'm trying to fit a fairly large hierarchical logit model (~50k rows, 100 columns, 50 groups). Fitting a normal logit model to the data (with dummy variables for group) works fine, but the hierarchical model appears to be getting stuck: the first optimization phase completes fine, but the second goes through a lot of iterations without anything changing and without stopping.

EDIT: I suspect the problem is mainly that I have so many parameters, because when I try to set maxfn to a lower value it gives a warning:

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

However, the parameter estimates aren't changing at all over the course of the optimization, so I'm still confused about what to do. When I tried to set maxfn in the optimizer controls (despite the warning), it seemed to hang after finishing the optimization.

Here's some code that reproduces the problem for random data:

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)

This outputs:

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>

I tried setting ncol to other values, and it appears that the number of iterations done is (approximately) 40 per column. Obviously, this becomes a huge pain as I add more columns. Are there tweaks I can make to the optimization algorithm that will reduce the dependence on the number of columns?

$\endgroup$
  • 1
    $\begingroup$ It would be helpful to know the specific model that you're trying to fit (especially the random effects structure). $\endgroup$ – Patrick S. Forscher Jan 9 '15 at 20:36
  • $\begingroup$ Unfortunately the precise model is proprietary. There's one level of random effects, with group sizes ranging between ~100 and 5000. Let me know if I can provide any other relevant info about the model. $\endgroup$ – Ben Kuhn Jan 9 '15 at 21:16
  • $\begingroup$ OK, I've added some code that reproduces the problem. $\endgroup$ – Ben Kuhn Jan 9 '15 at 22:13
  • 1
    $\begingroup$ I don't have a full answer for you, so I'll leave this as a comment. In my experience, glmer is quite slow, especially for models that have a complex random effects structure (e.g., many random slopes, crossed random effects, etc.). My first suggestion would be to try again with a simplified random effects structure. However, if you're experiencing this problem with a random intercepts model only, your problem may simply be the number of cases, in which case you'll need to try some tools specialized for big data. $\endgroup$ – Patrick S. Forscher Jan 9 '15 at 22:28
  • $\begingroup$ It has the same problem with 2 groups instead of 50. Also, testing with a smaller number of columns, it seems as though the number of iterations is roughly linear in the number of columns... Are there optimization methods that will do better here? $\endgroup$ – Ben Kuhn Jan 9 '15 at 23:47
12
+50
$\begingroup$

One thing you could try is to change the optimizer. See Ben Bolker's comment at this github issue. The nlopt implementation of bobyqa is usually much faster than the default (at least whenever I try it).

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

Also, see this answer for more options and this thread from R-sig-mixed-models (which looks more relevant to your issue).

Edit: I gave you some out-of-date info related to nloptr. In lme4 1.1-7 and up, nloptr is automatically imported (see ?nloptwrap). All you have to do is add

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

to your call.

$\endgroup$
  • $\begingroup$ Thank you! I'm trying the nlopt code right now. I do wonder if there's something other than a bad optimizer implementation going on, since fitting an almost-equivalent dummified glm was so much faster, but I'll see... $\endgroup$ – Ben Kuhn Jan 15 '15 at 21:55
  • $\begingroup$ Well, it was certainly faster, but it stopped with an error: PIRLS step-halvings failed to reduce deviance in pwrssUpdate. Do you have any idea what might be going on here? The error message isn't exactly transparent... $\endgroup$ – Ben Kuhn Jan 15 '15 at 22:00
  • $\begingroup$ For kicks, you could try setting nAGQ = 0 (see the thread I linked for a few more ideas). I don't remember what causes the PIRLS error, but I'll look around. $\endgroup$ – alexforrence Jan 15 '15 at 22:03
  • $\begingroup$ Thanks so much! Could you point me to a resource where I could learn more about the details of these methods so that I could solve problems like this myself in the future? Optimization feels very much like black magic to me at the moment. $\endgroup$ – Ben Kuhn Jan 15 '15 at 22:11
  • 2
    $\begingroup$ nAGQ = 0 worked for me on your test example with the default bobyqa (ran in ~15 sec), and in 11 sec with the nloptr bobyqa. Here's an interview with John C. Nash (co-author of the optim and optimx packages) where he does a high-level explanation of optimization. If you look up optimx or nloptr on CRAN, their respective reference manuals will tell you more about the syntax. nloptr also has a vignette available, which goes a little further into detail. $\endgroup$ – alexforrence Jan 16 '15 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.