The metafor Package

The Mantel-Haenszel method is an approach for fitting meta-analytic equal-effects models when dealing with studies providing data in the form of 2x2 tables or in the form of event counts (i.e., person-time data) for two groups (Mantel & Haenszel, 1959). The method is particularly advantageous when aggregating a large number of studies with small sample sizes (the so-called sparse data or increasing strata case).

In Cochrane Reviews, one often finds results that are based on a 'random-effects Mantel-Haenszel method', which is a bit confusing, given that the Mantel-Haenszel method is not based on a random-effects model. When looking at the details (see also the relevant section from the Cochrane Handbook), it turns out that this method is really just the standard standard random-effects model, but when estimating the amount of heterogeneity (i.e., $\tau^2$) using the DerSimonian-Laird (DL) estimator, the method uses the Q-statistic from the Mantel-Haenszel method, as opposed to using the one from a standard equal-effects model. The difference is typically negligible. Here is a small example to illustrate this.

We will use a small dataset from a Cochrane review on the use of antibiotics for the treatment of sore throat (Spinks et al., 2021). In Analysis 1.4.1, we see the results from three studies where antipyretics were also administered:

dat <- read.table(header=TRUE, text = "
study year treat.events treat.total placebo.events placebo.total
'Brumfitt' 1957 21 42 26 40
'Chapple' 1956 40 135 37 65
'De Meyere' 1992 18 82 59 91
")
dat

      study year treat.events treat.total placebo.events placebo.total
1  Brumfitt 1957           21          42             26            40
2   Chapple 1956           40         135             37            65
3 De Meyere 1992           18          82             59            91

The event counts are the number of participants who still had a sore throat on day 3 after either receiving antibiotics or placebo. Let's compute (log) risk ratios for the three studies:

library(metafor)
dat <- escalc(measure="RR", ai=treat.events,   n1i=treat.total,
                            ci=placebo.events, n2i=placebo.total,
                            data=dat)
dat

      study year treat.events treat.total placebo.events placebo.total      yi     vi
1  Brumfitt 1957           21          42             26            40 -0.2624 0.0373
2   Chapple 1956           40         135             37            65 -0.6529 0.0292
3 De Meyere 1992           18          82             59            91 -1.0830 0.0493

Now we can fit a standard random-effects model to the log risk ratios with:

res <- rma(yi, vi, data=dat, method="DL")
res

Random-Effects Model (k = 3; tau^2 estimator: DL)
 
tau^2 (estimated amount of total heterogeneity): 0.1095 (SE = 0.1482)
tau (square root of estimated tau^2 value):      0.3309
I^2 (total heterogeneity / total variability):   74.38%
H^2 (total variability / sampling variability):  3.90
 
Test for Heterogeneity:
Q(df = 2) = 7.8053, p-val = 0.0202
 
Model Results:
 
estimate      se     zval    pval    ci.lb    ci.ub
 -0.6552  0.2219  -2.9534  0.0031  -1.0901  -0.2204  **
 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note that the DL estimator of $\tau^2$ was used above, since this will also be the method to estimate the amount of heterogeneity when we apply the random-effects Mantel-Haenszel method. As a reminder, the DL estimator is computed with $$\hat{\tau}^2 = \frac{Q - (k-1)}{\sum w_i - \frac{\sum w_i^2}{\sum w_i}},$$ where $Q$ is the value of the Q-statistic ($Q$ = 7.8053 in the output above), $k$ is the number studies (3 in this example), and $w_i$ is the inverse sampling variance of the $i$th study. Let's see if we can compute the estimate of $\tau^2$ manually using this equation:

wi <- 1 / dat$vi
tau2 <- (res$QE - (res$k-1)) / (sum(wi) - sum(wi^2) / sum(wi))
tau2

0.1095074

This is indeed the estimate of $\tau^2$ given in the output above.

For easier interpretation of the results, we can back-transform the pooled estimate and confidence interval bounds via exponentiation:

predict(res, transf=exp, digits=2)

 pred ci.lb ci.ub pi.lb pi.ub
 0.52  0.34  0.80  0.24  1.13

Sidenote: With just three studies, fitting a random-effects model is quite 'ambitious', but we will ignore this issue for the purposes of this illustration.

To apply the random-effects Mantel-Haenszel method, we first use the Mantel-Haenszel method with:

res <- rma.mh(measure="RR", ai=treat.events,   n1i=treat.total,
                            ci=placebo.events, n2i=placebo.total,
                            data=dat, digits=2)
res

Equal-Effects Model (k = 3)
 
I^2 (total heterogeneity / total variability):  75.78%
H^2 (total variability / sampling variability): 4.13
 
Test for Heterogeneity:
Q(df = 2) = 8.26, p-val = 0.02
 
Model Results (log scale):
 
estimate    se   zval  pval  ci.lb  ci.ub
   -0.71  0.11  -6.25  <.01  -0.93  -0.48
 
Model Results (RR scale):
 
estimate  ci.lb  ci.ub
    0.49   0.40   0.62

These results are based on an equal-effects model, so we are not particularly interested in them in this example. However, we need the value of the Q-statistic from this model, so we can compute the DL estimator as described earlier:

wi <- 1 / dat$vi
tau2 <- (res$QE - (res$k-1)) / (sum(wi) - sum(wi^2) / sum(wi))
tau2

0.1180208

The value is slightly different compared to the one we obtained above. Now we can use this estimate of $\tau^2$ for the purposes of fitting a random-effects model with:

res <- rma(yi, vi, data=dat, tau2=tau2)
res

Random-Effects Model (k = 3; user-specified tau^2 value)
 
tau^2 (specified amount of total heterogeneity): 0.1180 (SE = 0.1563)
tau (square root of specified tau^2 value):      0.3435
I^2 (total heterogeneity / total variability):   75.78%
H^2 (total variability / sampling variability):  4.13
 
Test for Heterogeneity:
Q(df = 2) = 7.8053, p-val = 0.0202
 
Model Results:
 
estimate      se     zval    pval    ci.lb    ci.ub
 -0.6558  0.2282  -2.8740  0.0041  -1.1030  -0.2086  **
 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The back-transformed results are:

predict(res, transf=exp, digits=2)

 pred ci.lb ci.ub pi.lb pi.ub
 0.52  0.33  0.81  0.23  1.16

The results are essentially the same as what we obtained with the standard random-effects model. Of course one can construct examples where the difference is more noticeable.

Personally, I think the name of the method is a bit of a misnomer, as it has very little to do with the Mantel-Haenszel method (which is a method based on an equal/fixed-effects model). In a way, it is yet just another method for estimating $\tau^2$ (based on the DL estimator, but using a slightly different Q-statistic) and I am not aware of any simulation studies that have ever examined this method in detail. Certainly, the advantageous properties of the Mantel-Haenszel method when dealing with sparse data won't apply to the random-effects Mantel-Haenszel method.

Mantel, N., & Haenszel, W. (1959). Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute, 22(4), 719–748.

Spinks, A., Glasziou, P. P., & Mar, C. B. D. (2021). Antibiotics for treatment of sore throat in children and adults. Cochrane Database of Systematic Reviews, 12, CD000023.