How to visually assess the convergence of a mixed-effects model by plotting various optimizers

Question

To assess whether convergence warnings render the results invalid, or on the contrary, the results can be deemed valid in spite of the warnings, Bates et al. (2023) suggest refitting models affected by convergence warnings with a variety of optimizers. The authors argue that, if the different optimizers produce practically-equivalent results, the results are valid. The allFit function from the ‘lme4’ package allows the refitting of models using a number of optimizers. To use the seven optimizers listed above, two extra packages must be installed: ‘dfoptim’ and ‘optimx’ (see lme4 manual). The output from allFit() contains several statistics on the fixed and the random effects fitted by each optimizer.

library(lme4)
#> Loading required package: Matrix
library(dfoptim)
library(optimx)


# Create data using code by Ben Bolker from 
# https://stackoverflow.com/a/38296264/7050882

set.seed(101)
spin = runif(600, 1, 24)
reg = runif(600, 1, 15)
ID = rep(c("1","2","3","4","5", "6", "7", "8", "9", "10"))
day = rep(1:30, each = 10)
testdata <- data.frame(spin, reg, ID, day)
testdata$fatigue <- testdata$spin * testdata$reg/10 * rnorm(30, mean=3, sd=2)

# Model
fit = lmer(fatigue ~ spin * reg + (1|ID),
           data = testdata, REML = TRUE)

# Refit model using all available algorithms
multi_fit = allFit(fit)
#> bobyqa : [OK]
#> Nelder_Mead : [OK]
#> nlminbwrap : [OK]
#> nmkbw : [OK]
#> optimx.L-BFGS-B : [OK]
#> nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
#> nloptwrap.NLOPT_LN_BOBYQA : [OK]

# Show results 
summary(multi_fit)$fixef
#>                               (Intercept)      spin       reg  spin:reg
#> bobyqa                          -2.975678 0.5926561 0.1437204 0.1834016
#> Nelder_Mead                     -2.975675 0.5926559 0.1437202 0.1834016
#> nlminbwrap                      -2.975677 0.5926560 0.1437203 0.1834016
#> nmkbw                           -2.975678 0.5926561 0.1437204 0.1834016
#> optimx.L-BFGS-B                 -2.975680 0.5926562 0.1437205 0.1834016
#> nloptwrap.NLOPT_LN_NELDERMEAD   -2.975666 0.5926552 0.1437196 0.1834017
#> nloptwrap.NLOPT_LN_BOBYQA       -2.975678 0.5926561 0.1437204 0.1834016

^{Created on 2023-06-24 with reprex v2.0.2}

Numbers are hard for the common people, however. We are far better at vision. So, could there be a way to visualise the output from allFit(), and look at the parameters for a set of predictors across various optimizers (e.g., bobyqa, Nelder-Mead, etc.)?

Works for me without explicitly loading dfoptim and optimx, though. — jay.sf, Jun 25 '23 at 11:10

score 2 · Answer 1 · answered Jun 25 '23 at 11:34

Just sapply fixef over the multi_fit object and plot it with Map.

pframe <- as.data.frame(t(sapply(multi_fit, fixef)))

layout(matrix(c(1:5, 5), 3, 2, byrow=TRUE), heights=c(2, 2, 1.5))
par(mar=c(4, 4, 3, 3))

Map(plot, as.list(pframe), ylab=colnames(pframe), col=list(seq_len(nrow(pframe))), pch=3)
plot.new(); legend('top', ncol=2, pch=3, legend=rownames(pframe), col=seq_len(nrow(pframe)), cex=.9)

Pablo Bernabeu · Accepted Answer · 2023-06-27T22:17:05.303

Yes, there is a way, using a custom function called plot.fixef.allFit, described at https://pablobernabeu.github.io/2023/a-new-function-to-plot-convergence-diagnostics-from-lme4-allfit.

library(lme4)
#> Loading required package: Matrix
library(dfoptim)
library(optimx)

# Create data using code by Ben Bolker from 
# https://stackoverflow.com/a/38296264/7050882

set.seed(101)
spin = runif(600, 1, 24)
reg = runif(600, 1, 15)
ID = rep(c("1","2","3","4","5", "6", "7", "8", "9", "10"))
day = rep(1:30, each = 10)
testdata <- data.frame(spin, reg, ID, day)
testdata$fatigue <- testdata$spin * testdata$reg/10 * rnorm(30, mean=3, sd=2)

# Model
fit = lmer(fatigue ~ spin * reg + (1|ID),
           data = testdata, REML = TRUE)

# Refit model using all available algorithms
multi_fit = allFit(fit)
#> bobyqa : [OK]
#> Nelder_Mead : [OK]
#> nlminbwrap : [OK]
#> nmkbw : [OK]
#> optimx.L-BFGS-B : [OK]
#> nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
#> nloptwrap.NLOPT_LN_BOBYQA : [OK]

# Show results 
summary(multi_fit)$fixef
#>                               (Intercept)      spin       reg  spin:reg
#> bobyqa                          -2.975678 0.5926561 0.1437204 0.1834016
#> Nelder_Mead                     -2.975675 0.5926559 0.1437202 0.1834016
#> nlminbwrap                      -2.975677 0.5926560 0.1437203 0.1834016
#> nmkbw                           -2.975678 0.5926561 0.1437204 0.1834016
#> optimx.L-BFGS-B                 -2.975680 0.5926562 0.1437205 0.1834016
#> nloptwrap.NLOPT_LN_NELDERMEAD   -2.975666 0.5926552 0.1437196 0.1834017
#> nloptwrap.NLOPT_LN_BOBYQA       -2.975678 0.5926561 0.1437204 0.1834016

# Read in function from GitHub
source('https://raw.githubusercontent.com/pablobernabeu/plot.fixef.allFit/main/plot.fixef.allFit.R')

plot.fixef.allFit(multi_fit, 
                  
                  select_predictors = c('spin', 'reg', 'spin:reg'), 
                  
                  # Increase padding at top and bottom of Y axis
                  multiply_y_axis_limits = 1.3,
                  
                  y_title = 'Fixed effect (*b*)')
#> Loading required package: dplyr
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
#> Loading required package: reshape2
#> Loading required package: stringr
#> Loading required package: scales
#> Loading required package: ggplot2
#> Loading required package: ggtext
#> Loading required package: patchwork


# Alternative using plot-specific Y axes and other modified settings

plot.fixef.allFit(multi_fit, 
                  
                  select_predictors = c('spin', 'spin:reg'), 
                  
                  # Use plot-specific Y axis limits
                  shared_y_axis_limits = FALSE,
                  
                  decimal_places = 7, 
                  
                  # Move up Y axis title
                  y_title_hjust = 4.5,
                  
                  y_title = 'Fixed effect (*b*)')

^{Created on 2023-06-26 with reprex v2.0.2}

How to visually assess the convergence of a mixed-effects model by plotting various optimizers

2 Answers2