What are the effects of continuity correction in binomial meta-analysis in the metafor package?

Question

The metafor package in R can be used to fit the random effects model with the exact binomial likelihood as described in Stijnen et al. (2010). The metafor website also shows how to reproduce the examples from that paper (link). One of the benefits of using the binomial likelihood instead of normal approximations is that studies with zero events can be properly modelled instead of adding a continuity correction of 0.5 to those studies.

The documentation for the appropriate rma.glmm function says that the continutity correction is added by default, trough the add argument. Here is the problem: Trying to do the analysis without the continuity correction gives an error. From the metafor webpage we can do the analysis from the Stijnen et al (2010) paper that works fine:

library(metafor) # version 1.9-7
dat <- get(data(dat.nielweise2007))
rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat) # add=1/2 by default.

If instead the continuity correction is removed an error occurs:

rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat, add=0) # add=0 means no correction.

Error in model.frame.default(formula = yi ~ X - 1, drop.unused.levels = TRUE) : 
  variable lengths differ (found for 'X')
In addition: Warning messages:
1: In escalc.default(measure = measure, xi = xi, mi = mi, add = add,  :
  Some yi and/or vi values equal to +-Inf. Recoded to NAs.
2: In rma.glmm(measure = "PLO", xi = ci, ni = n2i, data = dat, add = 0) :
  Some yi/vi values are NA.

I am able to reproduce the analysis by using the lme4 package without the correction and get the same estimates:

library(lme4) # version 1.1-7

# nAGQ=7 is what the metafor package uses 
mm <- glmer(cbind(ci, n2i-ci) ~ (1|study), family=binomial(link='logit'), data=dat, nAGQ=7)
summary(mm)

Why does it not work to remove the continuity correction in the rma.glmm function? I also tried to add a smaller correction (add=0.0000001) and the estimates for tau, I^2 and H^2 changed a bit, but the numbers under the model results heading did not. What is going on?

References:

Stijnen, T., Hamza, T. H., & Ozdemir, P. (2010). Random effects meta-analysis of event outcome in the framework of the generalized linear mixed model with applications in sparse data. Statistics in Medicine, 29(29), 3046–3067.

Answer 1

There are two separate issues here, one being a slight misunderstanding of how continuity corrections are used in the metafor package when fitting so-called "binomial-normal models" with the rma.glmm() function, the other being a potential bug report.

Let me start with the first issue, namely how continuity corrections come into play when using the rma.glmm() function. Actually, let me quote from the function documentation (help(rma.glmm)):

The various models do not require the calculation of the observed outcomes of the individual studies (e.g., the observed odds ratios of the k studies) and directly make use of the table/event counts. Zero cells/events are not a problem (except in extreme cases, such as when one of the two outcomes never occurs or when there are no events in any of the studies). Therefore, it is unnecessary to add some constant to the cell/event counts when there are zero cells/events.

However, for plotting and various other functions, it is necessary to calculate the observed outcomes for the k studies. Here, zero cells/events can be problematic, so adding a constant value to the cell/event counts ensures that all k values can be calculated. The add and to arguments are used to specify what value should be added to the cell/event counts and under what circumstances when calculating the observed outcomes. The documentation of the escalc function explains how the add and to arguments work. Note that drop00 is set to TRUE by default, since studies where ai=ci=0 or bi=di=0 or studies where x1i=x2i=0 are uninformative about the size of the effect.

Some parts of this text (e.g., the last sentence) are only relevant when analyzing two-group data (e.g., 2x2 table data), but the point is hopefully clear: No continuity corrections are applied when actually fitting the various types of models that are implemented in the rma.glmm() function. Continuity corrections only come into play when calculating the observed outcomes (e.g., the observed log odds or log odds ratios), for example for plotting purposes.

So, no matter what value you specify for add, you should always get the same model results. The following illustrates that this is indeed the case:

library(metafor)
dat <- get(data(dat.nielweise2007))
res1 <- rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat) # note: 'add=1/2' is the default
res2 <- rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat, add=1)
res3 <- rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat, add=.1)
res1
res2
res3

But forest plots of the individual outcomes (with the model-based summary at the bottom) using forest() shows how the add argument is affecting the display of the data:

forest(res1)
forest(res2)
forest(res3)

Only study 15 had 0 events, so focus on that one.

I hope this clarifies what is going on.

Now the second part is about what happens when you set add=0. Note that doing so would imply that the log odds (and corresponding sampling variance) for study 15 cannot be calculated. Compare:

escalc(measure="PLO", xi=ci, ni=n2i, data=dat, subset=14:16) # note: 'add=1/2' is again the default

yields:

  study  author year ai n1i ci n2i      yi     vi
1    14    Leon 2004  6 187 11 180 -2.7320 0.0968
2    15   Yucel 2004  0 118  0 105 -5.3519 2.0095
3    16 Moretti 2005  0 252  1 262 -5.5645 1.0038

And:

escalc(measure="PLO", xi=ci, ni=n2i, data=dat, subset=14:16, add=0)

yields:

  study  author year ai n1i ci n2i      yi     vi
1    14    Leon 2004  6 187 11 180 -2.7320 0.0968
2    15   Yucel 2004  0 118  0 105      NA     NA
3    16 Moretti 2005  0 252  1 262 -5.5645 1.0038
Warning message:
In escalc.default(measure = "PLO", xi = ci, ni = n2i, data = dat,  :
  Some yi and/or vi values equal to +-Inf. Recoded to NAs.

So, the same thing will happen if you use add=0 in rma.glmm(), so if you were to plot the results, for example with the forest() function, you would not have an observed outcome for study 15. Now this may be exactly what you want -- or you may still want to show that there actually was a study 15, but not have an observed outcome computed for it based on some continuity correction. You can find ways of handling such cases/situations described here:

http://www.metafor-project.org/doku.php/tips:handling_missing_data

However, actually running:

rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat, add=0)

indeed yields:

Error in model.frame.default(formula = yi ~ X - 1, drop.unused.levels = TRUE) : 
  variable lengths differ (found for 'X')
In addition: Warning messages:
1: In escalc.default(measure = measure, xi = xi, mi = mi, add = add,  :
  Some yi and/or vi values equal to +-Inf. Recoded to NAs.
2: In rma.glmm(measure = "PLO", xi = ci, ni = n2i, data = dat, add = 0) :
  Some yi/vi values are NA.

So, oops. That's a bug. I'll take a look right away and hopefully will have a fix soon.

Update: This has been fixed for the next version of the package (1.9-8) to be released in the future. A development version of the metafor package that already incorporates the fix can be downloaded as described here:

http://www.metafor-project.org/doku.php/installation#development_version

Also, one clarification. Whether you use add=0 (or some other value) will not have an influence on the results, except that the $I^2$ and $H^2$ values can change, since these statistics involve the sampling variances of the individual studies, and how those are computed is indeed influenced by what you specify for the add argument. In fact, when you set add=0, then the sampling variances for studies with no events cannot be computed, so the $I^2$ and $H^2$ statistics are then based on only the remaining studies for which the sampling variances can actually be computed.