# The metafor Package

A Meta-Analysis Package for R

### Site Tools

tips:convergence_problems_rma_mv

## Convergence Problems with the rma.mv() Function

Model fitting with the rma.mv() function makes use of iterative methods that attempt to maximize the log likelihood function over the model parameters. The optimization techniques that are used for this purpose sometimes do not converge, in which case the function will inform the user that the optimizer did not achieve convergence. Below, I illustrate this problem with an example and discuss some remedies.

### Example

Consider the following dataset:

library(metafor)
dat <- structure(list(study = c(1L, 2L, 2L, 2L, 2L, 2L, 3L, 4L, 5L, 6L, 6L,
6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 12L,
12L, 12L, 12L, 13L, 14L, 14L, 15L, 15L, 16L, 17L, 17L, 17L, 17L, 18L, 19L,
19L, 20L), id = 1:46, yij = c(0.686, 0.672, 0.742, 0.914, 0.861, 0.676, 0.711,
0.688, 0.355, 0.774, 1.067, 0.869, 1.121, 0.328, 0.559, 0.62, 0.442, 0.299,
0.387, 0.376, 0.44, 0.414, 0.414, 0.39, 0.061, 0.334, -0.015, 0.508, 0.958,
0.256, 0.539, 0.941, 0.57, 0.308, 0.274, 0.834, 0.633, 0.195, 0.504, 0.216,
0.545, 0.069, 0.097, 0.662, 0.633, 0.317 ), vij = c(0.085, 0.007, 0.085,
0.026, 0.076, 0.028, 0.074, 0.024, 0.085, 0.027, 0.053, 0.054, 0.055, 0.045,
0.014, 0.078, 0.095, 0.028, 0.029, 0.084, 0.017, 0.064, 0.084, 0.05, 0.067,
0.059, 0.092, 0.027, 0.082, 0.096, 0.086, 0.049, 0.056, 0.042, 0.016, 0.068,
0.006, 0.085, 0.041, 0.048, 0.073, 0.048, 0.066, 0.009, 0.016, 0.085)), class
= "data.frame", row.names = c(NA, -46L))
dat
   study id   yij   vij
1      1  1 0.686 0.085
2      2  2 0.672 0.007
3      2  3 0.742 0.085
4      2  4 0.914 0.026
5      2  5 0.861 0.076
6      2  6 0.676 0.028
7      3  7 0.711 0.074
8      4  8 0.688 0.024
9      5  9 0.355 0.085
10     6 10 0.774 0.027
11     6 11 1.067 0.053
12     6 12 0.869 0.054
13     .  .     .     .

Variable yij contains the effect sizes and variable vij the corresponding sampling variances. The dataset has a multilevel structure, with some studies providing multiple effect size estimates (e.g., study 2 provides 5 estimates). Leaving aside further details, a sensible model for these data may be a multilevel random-effects model, which we can fit with:

res <- rma.mv(yij, vij, random = ~ 1 | study / id, data=dat)

Unfortunately, this yields the error:

Error : Optimizer (nlminb) did not achieve convergence (convergence = 1).

Hence, the algorithm used for the model fitting (which by default is the nlminb() function, which uses a Quasi-Newton method for the optimization) did not converge with the default settings.

By setting verbose=TRUE, we can obtain information on the progress of the optimization algorithm:

res <- rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, verbose=TRUE)
Iteration 0     ll = -0.6174   sigma2 = 0.0107  0.0107
Iteration 1     ll = -0.6174   sigma2 = 0.0107  0.0107
...
Iteration 7     ll =  0.7454   sigma2 = 0.0224  0.0054
...
Iteration 14    ll =  1.3502   sigma2 = 0.0283  0.0015
...
Iteration 21    ll =  1.4727   sigma2 = 0.0222  0.0008
...
Iteration 56    ll =  1.6225   sigma2 = 0.0252  0.0000
...
Iteration 124   ll =  1.6284   sigma2 = 0.0252  0.0000
Iteration 125   ll =  1.6285   sigma2 = 0.0252  0.0000
...
Iteration 126   ll =  1.6285   sigma2 = 0.0252  0.0000
Iteration 306   ll =  1.6285   sigma2 = 0.0252  0.0000

Starting around iteration 126, there is no further change in the log likelihood when rounded to 4 decimal places (and neither in the two variance components that are estimated), but nlminb() uses a much smaller threshold for determining when convergence has occurred.

Via the control argument, one can adjust the settings of the optimizer that is used. For example, reading help(nlminb) reveals that one can adjust the number of iterations via the iter.max control argument (which by default is 150). So let's try increasing this to something larger:

res <- rma.mv(yij, vij, random = ~ 1 | study / id, data=dat,
verbose=TRUE, control=list(iter.max=1000))

The output (not shown) is the same as before, so this didn't really help with achieving convergence. The problem in this case is the rather small threshold that is used by default for determining whether the change in the log likelihood is small enough to declare that convergence has occurred. As we see above, convergence actually seems to occur, but is not detected by the algorithm. Let's try a slightly less strict threshold by adjusting rel.tol (which by default is 1e-10):

res <- rma.mv(yij, vij, random = ~ 1 | study / id, data=dat,
verbose=TRUE, control=list(rel.tol=1e-8))
res

Now the error disappears and the output shows that convergence was reached after 226 iterations. We can then proceed to examine the results:

Multivariate Meta-Analysis Model (k = 46; method: REML)

Variance Components:

estim    sqrt  nlvls  fixed    factor
sigma^2.1  0.0252  0.1588     20     no     study
sigma^2.2  0.0000  0.0000     46     no  study/id

Test for Heterogeneity:
Q(df = 45) = 68.9392, p-val = 0.0124

Model Results:

estimate      se     zval    pval   ci.lb   ci.ub
0.5245  0.0488  10.7373  <.0001  0.4288  0.6203  ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We see that the second variance component ($\sigma^2_1$), which provides an estimate of the amount of within-study heterogeneity, converged to (essentially) 0. There is nothing inherently wrong with this (it could very well be the case that all of the heterogeneity in these data is due to differences between studies and none due to differences in the effect sizes within studies), but whenever parameter estimates are drifting towards bounds like this, it does make life more difficult for the optimization routines.

The figure below shows the (restricted) log likelihood surface for these data/model as a function of the two variance components. The black line indicates the path that the optimizer takes over the surface, eventually reaching the maximum at the red dot (for an animated version of this figure, click here).

Instead of adjusting the settings of the default optimizer, we can also try switching to a different optimization routine. The rma.mv() function actually provides access to more than a dozen different optimizers, which can be chosen via the optimizer control argument. Below, I fit the same model as above, trying out 12 other optimizers (besides the default one). Where necessary, I adjusted some other settings to achieve convergence with all of the optimizers.

res <- list(
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="nlminb", rel.tol=1e-8)),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="Nelder-Mead")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="BFGS")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="bobyqa")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="nloptr", maxeval=1000)),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="nlm")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="hjk")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="nmk")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="mads")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="ucminf")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="lbfgsb3c")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="subplex")),
rma.mv(yij, vij, random = ~ 1 | study / id, data=dat, control=list(optimizer="BBoptim")))

We can then extract the log likelihood and estimates of the two variance components and put everything into a data frame as follows:

tab <- lapply(res, function(x) data.frame(optimizer=x$control$optimizer, ll=logLik(x),
sigma2.1=x$sigma2[1], sigma2.2=x$sigma2[2]))
tab <- do.call(rbind, tab)
dfround(tab, digits=5)
     optimizer      ll sigma2.1 sigma2.2
1       nlminb 1.62847  0.02522        0
3         BFGS 1.62830  0.02522        0
4       bobyqa 1.62847  0.02522        0
5       nloptr 1.62841  0.02522        0
6          nlm 1.62847  0.02522        0
7          hjk 1.62847  0.02522        0
8          nmk 1.62847  0.02523        0
13     BBoptim 1.62847  0.02522        0
Given that some other optimizers converged without any adjustments to their default settings, one might wonder why nlminb() is the default optimizer in rma.mv() (which did require some tweaking to achieve convergence in this example). The example above was constructed specifically to lead to non-convergence with nlminb() under its default settings. For many other datasets, nlminb() will converge without difficulties, while some of the other optimizers might struggle. There is simply not a single algorithm that works flawlessly in all cases. In my experience, nlminb() is actually a quite reasonable and robust default choice and typically works without problems.