Forest and funnel plots

The forest plot: what each element encodes

A forest plot looks like a pile of dots and lines, but every visual element carries specific information:

Box size is the element most often ignored and the one carrying the most information. It tells you who is actually driving the pooled result. A common situation: one very large trial’s box is big enough to take most of the weight, while a dozen small studies together contribute little — at which point saying “we pooled fifteen studies” is somewhat misleading, because what you are really looking at is close to the result of that one large trial.

library(metafor)
library(metadat)

data(dat.bcg, package = "metadat")

# Turn the 2x2 tables into a log risk ratio and its variance for each trial
bcg <- escalc(measure = "RR",
              ai = tpos, bi = tneg, ci = cpos, di = cneg,
              data = dat.bcg, slab = paste(author, year, sep = ", "))

# Random-effects model (REML)
res <- rma(yi, vi, data = bcg, method = "REML")
res                                   # pooled estimate, I2, tau2, Q test

forest(res, atransf = exp,            # draw on the RR scale
       at = log(c(0.05, 0.25, 1, 4)),
       header = c("Trial", "Risk ratio [95% CI]"))

import numpy as np
from statsmodels.stats.meta_analysis import combine_effects

# The 2x2 counts of metadat::dat.bcg, 13 trials, carried here so this
# block runs on its own.
tpos = np.array([4, 6, 3, 62, 33, 180, 8, 505, 29, 17, 186, 5, 27])
tneg = np.array([119, 300, 228, 13536, 5036, 1361, 2537, 87886, 7470, 1699,
                 50448, 2493, 16886])
cpos = np.array([11, 29, 11, 248, 47, 372, 10, 499, 45, 65, 141, 3, 29])
cneg = np.array([128, 274, 209, 12619, 5761, 1079, 619, 87892, 7232, 1600,
                 27197, 2338, 17825])
rr        = (tpos / (tpos + tneg)) / (cpos / (cpos + cneg))
se_log_rr = np.sqrt(1 / tpos - 1 / (tpos + tneg) + 1 / cpos - 1 / (cpos + cneg))

# Input is the log RR and its variance (computed just above)
eff = np.log(rr)            # each study's log risk ratio
var = se_log_rr ** 2        # the matching variance

res = combine_effects(eff, var, method_re="dl")   # DerSimonian-Laird
print(res.summary_frame())
print("I2 =", res.i2)

The pooled result is RR 0.49 (95% CI 0.34–0.70) across 13 studies. Every box and every horizontal line in that plot is one study’s log RR together with its variance — where those two numbers come from, and why it is the variance that decides how big the box is, is on effect measures and their variances.

Fixed effect and random effects ask different questions

The diamond above came out of a random-effects model, and which model drew it changes what it means: a fixed-effect diamond estimates one common true effect, a random-effects diamond estimates the mean of a distribution of study-specific effects. Reading a forest plot therefore starts with finding out which of the two produced the diamond — and with noticing that random effects also spreads the weights more evenly, so the boxes are not the sizes a fixed-effect plot would have drawn.

Which assumption is defensible for a given set of trials, what the choice costs in interval width, and which τ² estimator was used are the subject of fixed-effect and random-effects models; this page does not repeat them.

I²: it measures a proportion, not a magnitude

Most forest plots print I² beneath the diamond — here 92.2% — usually alongside the Q test (here Q = 152.2, df = 12, p = < 0.001). Both are quoted far more often than they are read correctly, and each has one standard misreading you need before the plot makes sense.

I² is the percentage of the total variation that is between-study variation — a ratio, so a large I² says nothing about how large the differences between studies actually are. For that, look at τ² (the between-study variance, 0.313 here) or at a prediction interval. And the Q test has very little power when the studies are few, so failing to reject its null is not evidence that the studies agree.

What each statistic does mean, the two estimators of I² that do not give the same number, and a simulation of how little power that Q test really has are all on heterogeneity: I², τ² and prediction intervals.

And the reason the BCG dataset is a teaching classic is that its heterogeneity has an explicable source: the latitude of the trial site. BCG works markedly better at higher latitudes (generally attributed to exposure to environmental non-tuberculous mycobacteria). That is the right way to handle heterogeneity — not to look for a way of suppressing it, but to find where it comes from and explain it with meta-regression or subgroup analysis.

Funnel plots and asymmetry

A funnel plot puts each study’s effect size on the x-axis against its precision on the y-axis (usually the standard error, inverted so that the precise studies sit at the top).

The logic is: in the absence of bias, large precise studies cluster near the true value (the neck of the funnel) and small imprecise studies scatter symmetrically on either side (the mouth), giving an inverted funnel shape.

Egger’s test formalises the eyeballing: regress each study’s effect estimate on its own standard error and ask whether the slope differs from zero. A symmetric funnel puts that slope at zero — the estimate does not depend on how precise the study was; a slope away from zero says the smaller studies sit systematically to one side. (In this parameterisation the intercept is the other useful quantity: it is the pooled estimate extrapolated to a study of infinite precision, which is where the Egger-adjusted estimate mentioned below comes from. Egger’s original write-up regresses the standard normal deviate on precision instead, and there it is the intercept that carries the test — same test, different algebra, which is why sources disagree about which coefficient to read.)

Pooling this dataset gives OR 0.46 (95% CI 0.31–0.70), with an Egger test p < 0.0001. Its I² is 73.6%, so between-study variation is one of the possible sources that has to be considered when explaining this asymmetry; I² is a relative proportion of variance and cannot on its own settle either the absolute clinical size of the heterogeneity or the cause of the asymmetry.

The intercept of that same regression — the estimate extrapolated to a study of infinite precision — is OR 1.105 on this dataset, on the other side of the null. It is an extrapolation with no test attached to it, useful only for showing which way the data point.

This is one of the most famous episodes in the history of meta-analysis: a set of small trials consistently showed a significant benefit of intravenous magnesium in myocardial infarction, and pooling them produced a large apparent effect; the very large trial that followed (ISIS-4, with over fifty thousand participants) detected no protective effect — in this dataset its point estimate even falls to the right of 1. The funnel plot’s asymmetry was already a warning sign before the large trial was published.

Egger’s is not the only asymmetry test

Egger’s test is the one people cite, but it is neither the only one nor the right one for every outcome type. In practice you meet three tests and one sensitivity analysis. All four are aimed at the same thing — does the reported effect track how large the study was — and differ in what stands in for “large” and in how the regression is weighted.

Here is what they say about the magnesium trials (16 of them).

Egger’s known bias for binary outcomes is why Peters exists

Look again at the list of small-study effects above: the fourth entry is the mathematical dependence between the effect measure and its variance. For an OR or an RR that is not an abstract warning. The point estimate of the log OR and its standard error are computed from the same cells of the same 2x2 table, so a trial that happens to observe a larger effect also gets a larger estimated standard error. Regressing the effect on the standard error therefore has a non-zero slope built into it even when nothing is biased at all — most visibly when events are rare, the arms are unbalanced, or the true effect is large. Egger’s test inflates its type I error rate accordingly.

Peters replaces the predictor with a quantity that does not contain the observed effect: one over the total sample size. The weights change too, to an effective sample size built from the marginal totals of events and non-events rather than from the inverse variance of the effect measure. On the same data Egger gives p < 0.0001 and Peters gives p = 0.0014 — the same direction, but the Peters p-value is more than an order of magnitude larger, and part of that gap is the artefactual correlation just described.

The three tests disagree on the same dataset

On these 16 trials, Egger and Peters both flag asymmetry and the rank correlation test does not (p = 0.4503). That is not a contradiction, it is a difference in power: the rank test gives up the linearity assumption and pays for it by needing more studies to detect the same degree of asymmetry. A non-significant rank test cannot be used to argue against the other two — it simply did not detect anything.

Run the same four methods on the BCG dataset from the first half of this page and everything flips:

On the BCG side all three tests detected no funnel plot asymmetry. That sentence cannot be extended: it is not the same as “there is no publication bias”, nor as “the funnel is symmetric”. There are only 13 studies here, sitting right on the Cochrane threshold below which asymmetry tests should not be run, and at that number all three tests have low power to begin with — low enough that “nothing detected” is close to the default outcome rather than information. Note also that the BCG dataset has an I² of 92.2%: extreme heterogeneity with no detectable asymmetry. Heterogeneity and funnel asymmetry are different things, and a high value of one does not imply a high value of the other.

Contour-enhanced funnels: which region the gap falls in

An ordinary funnel plot tells you where the gap is but not whether studies sitting in that gap would have been statistically significant — which is precisely the piece of information that decides whether the gap looks like publication bias. A contour-enhanced funnel paints the significance regions onto the background.

Centring the contours on the null is the whole point of the plot. Drawn around the pooled estimate instead — the default in several packages — the shading looks almost identical but answers a completely different question: how far each study sits from the pooled value, rather than whether each study was significant on its own.

There is only one reading rule, and it splits the five mechanisms into two piles:

• The gap is in the non-significant region — consistent with publication bias. The studies that go unpublished are exactly the small ones that failed to reach significance.
• The gap is in a significant region — publication bias does not explain it, because nobody suppresses a significant result. Look instead at heterogeneity, the methodological quality of small studies, the mathematical dependence between effect and variance, or plain chance.

Here is how this dataset falls out:

Of the 14 trials left of the null, 6 reach nominal significance; on the right there are only 2, and neither is significant. What is missing is the small trial pointing towards the control arm and failing to reach significance — the white region. That is the shape publication bias leaves behind, and it is why this dataset should have been unsettling before ISIS-4 was published.

Trim-and-fill imputes studies; it does not find them

The algorithm has three steps: trim off the extreme studies responsible for the asymmetry and re-estimate the pooled effect from what is left; mirror the trimmed studies across that new pooled value to fill in the other side; then re-pool using the completed set.

On this dataset the algorithm imputes 7 studies, all on the right, taking the count from 16 to 23. The pooled OR moves from 0.465 to 0.676, 95% CI 0.443–1.031 — the adjusted interval crosses the null. A protective effect that looked settled no longer reaches statistical significance once the funnel is filled.

Read that carefully. It is not evidence that intravenous magnesium does not work, for three reasons:

1. The 7 imputed points are not recovered studies. They have no patients, no event counts and no authors; they are coordinates the algorithm computed from a symmetry assumption. What the genuinely missing studies looked like, if they exist, nobody knows.
2. It assumes the asymmetry comes from missing studies. If the asymmetry is really heterogeneity — this dataset has an I² of 73.6% — then the imputed studies are compensating for a gap that is not there and pushing the estimate in a direction nothing supports.
3. It is sensitive to how τ² is estimated. Change the random-effects estimator and both the number imputed and the distance moved change with it.

So the output of trim-and-fill should be read as a conditional: if the asymmetry really is missing studies, the pooled estimate moves at most this far. That is one bound in a sensitivity analysis, not a corrected answer.

Applied to the BCG dataset the same procedure imputes only 1 study and moves the RR from 0.489 to 0.518 (95% CI 0.365–0.736), an interval that still excludes the null. When a funnel is not conspicuously asymmetric there is little for this algorithm to do — which is it behaving normally, not it confirming an absence of bias.

Run it yourself

library(metafor)
library(metadat)

data(dat.egger2001, package = "metadat")
mag <- escalc(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i,
              data = dat.egger2001, slab = study)
res <- rma(yi, vi, data = mag, method = "REML")

# 1. Egger: effect on standard error. The slope carries the test;
#    the intercept is the limit estimate at infinite precision
regtest(res, model = "lm")

# 2. Begg-Mazumdar: Kendall rank correlation of standardised effects
#    against their variances
ranktest(res)

# 3. Peters: predictor becomes 1/total sample size, weights become the
#    marginal event and non-event totals. metafor's
#    regtest(predictor = "ninv") swaps only the predictor and keeps 1/vi
#    weights, which is not Peters, so write out the published regression
a  <- dat.egger2001$ai
b  <- dat.egger2001$n1i - a
cc <- dat.egger2001$ci
d  <- dat.egger2001$n2i - cc
ni <- a + b + cc + d
wi <- 1 / (1 / (a + cc) + 1 / (b + d))
summary(lm(mag$yi ~ I(1 / ni), weights = wi))   # read the second row: slope

# 4. Trim-and-fill: not a test, a sensitivity analysis
tf <- trimfill(res)
tf                                    # k0 = how many imputed, side = which side
predict(tf, transf = exp)             # pooled estimate and interval after filling

# Contour-enhanced funnel. refline = 0 is the load-bearing argument: the
# contours must be centred on the null, because centred on the pooled
# estimate the same picture answers a different question
funnel(res, atransf = exp, refline = 0,
       level = c(90, 95, 99),
       shade = c("white", "gray88", "gray72"),
       legend = TRUE)

funnel(tf, atransf = exp)             # imputed studies get a different symbol

import numpy as np
import statsmodels.api as sm
from scipy.stats import kendalltau

# The magnesium trials of metadat::dat.egger2001, 16 studies. The
# continuity correction goes ONLY to the study with a zero cell, matching
# escalc's default add = 1/2, to = "only0". Adding half to every cell puts
# 14 of the 16 yi values out of step with R.
ai  = np.array([1, 9, 2, 1, 10, 1, 1, 0, 6, 1, 2, 5, 4, 90, 4, 2216])
n1i = np.array([40, 135, 200, 48, 150, 59, 25, 22, 76, 27, 89, 23, 130,
                1159, 107, 29011])
ci  = np.array([2, 23, 7, 1, 8, 9, 3, 1, 11, 7, 12, 13, 8, 118, 17, 2103])
n2i = np.array([36, 135, 200, 46, 148, 56, 23, 21, 75, 27, 80, 33, 122,
                1157, 108, 29039])
a, b = ai.astype(float), (n1i - ai).astype(float)
c, d = ci.astype(float), (n2i - ci).astype(float)
zero = (a == 0) | (b == 0) | (c == 0) | (d == 0)
a[zero] += 0.5; b[zero] += 0.5; c[zero] += 0.5; d[zero] += 0.5
yi = np.log((a * d) / (b * c))
vi = 1 / a + 1 / b + 1 / c + 1 / d

# yi = each study's log OR, vi = its variance, both computed just above
sei = np.sqrt(vi)

# 1. Egger: weighted least squares with 1/vi weights, testing the slope
egger = sm.WLS(yi, sm.add_constant(sei), weights=1 / vi).fit()
print(egger.tvalues[1], egger.pvalues[1])

# 2. Begg-Mazumdar: Kendall tau. metafor standardises the effects first,
#    so ranktest() will not give exactly this number
print(kendalltau(yi, vi))

# 3. Peters: predictor 1/n, weights from the marginal event and
#    non-event totals
ni = a + b + c + d
wi = 1 / (1 / (a + c) + 1 / (b + d))
peters = sm.WLS(yi, sm.add_constant(1 / ni), weights=wi).fit()
print(peters.tvalues[1], peters.pvalues[1])

# 4. trim-and-fill and contour-enhanced funnels: no maintained Python
#    implementation, so these two steps go back to R

Common misuses

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B7-04-forest-funnel.R

Watch next

▶How to Interpret a Forest Plot

▶醫學統計 EP18 統合分析：加權整合多項研究、看懂森林圖

▶統合分析的異質性（Heterogeneity）

▶統合分析異質性和小樣本研究偏誤

Report a problem on this page