[R] lme, lmer, convergence

[Ben Bolker]

Erin McMullen Jonaitis <jonaitis <at> wisc.edu> writes:

>  Hello, all, I am running some simulations to estimate power for a
> complicated epidemiological study, and am using lme and lmer to get
> these estimates.  I have to run a few thousand iterations, and once
> in a great while, an iteration will create fake data such that the
> model won't converge.  I see from Google searches that this is not
> an uncommon situation.
 
> My question: is there a way to extract the convergence value from an
> lme or lmer model?  It prints on the screen, but how can I get hold
> of it to evaluate it with some function?  What I'd like to do is
> build a failsafe into my program so that if a particular model in an
> iteration doesn't converge, I call a redo on that iteration.  This
> way the program will keep running and not stop in a fit of pique in
> the middle of my long simulation.
 
> If I can't do this, my fallback will be to try setting lmeControl
> options such that even "bad" models return parameter estimates etc
> -- once or twice in 10,000 iterations should not ruin things too
> badly -- but I'd like to try it the cleaner way first.

  There's a somewhat hack-ish solution, which is to use options(warn=2)
to 'upgrade' warnings to errors, and then use try() or tryCatch() to
catch them.

  More fancily, I used code that looked something like this to save
warnings as I went along (sorry about the <<- ) in a recent simulation
study.  You could also check w$message to do different things in the
case of different warnings.

withCallingHandlers(tryCatch(fun(n=nvec[j],tau=tauvec[i],...),
                error = function(e) {
                  warn[k,i,j] <<- paste("ERROR:",e$message)
              NA_ans}),
               warning = function(w) {
                  warn[k,i,j] <<- w$message
                  invokeRestart("muffleWarning")
             })

In the slightly longer run, we are working on getting the development
(lme4Eigen) version of lme4 to save convergence warnings in an
accessible slot.  I don't know if I would hold my breath for this
to be back-ported to nlme, though ...

  Some of these discussions might be better suited for r-sig-mixed-models
<at> r-project.org


  Ben Bolker