Sensitivity and influence analysis

Two different questions: sensitivity and influence

One of the most frequent methodological requests in a reviewer’s report is please add a sensitivity analysis. There are two different questions underneath that sentence, and they need different pictures:

The first is sensitivity analysis, covered in general on its own page; the second is influence analysis. This page is about the second, because it is the only place on this site where the included studies are examined one at a time. Leave-one-out belongs to both: dropping a data point is also changing an inclusion criterion.

This page uses metadat::dat.bcg, the 13 controlled trials of BCG vaccination against tuberculosis. It is the same dataset and the same pooled analysis as forest and funnel plots: random effects (REML), pooled RR 0.489 (95% CI 0.344–0.696, p < 0.001), I² 92.2%, τ² 0.3132, Q = 152.2 (df = 12, p < 0.001). The script behind this page pins that pooled value with an assertion, so the two pages cannot drift apart.

Leave-one-out: drop one trial at a time

The most direct of the three. Remove each of the 13 trials in turn, refit the random-effects model, and see where the pooled estimate goes.

Three things are worth reading off it:

1. The pooled estimate moves between RR 0.4517 and 0.5334, against 0.4894 for the full set; the largest single relative change is 9.0% (omitting Hart & Sutherland, 1977).
2. All 13 leave-one-out intervals still exclude 1. No single trial is holding the conclusion up.
3. The heterogeneity does not go away whichever trial is dropped. I² is lowest at 87.0% (omitting TPT Madras, 1980) and highest at 93.2% (omitting Aronson, 1948), against 92.2% for the full set — a range far smaller than the quantity itself. It is not the doing of any one trial, which is the concrete form here of the argument on the heterogeneity page.

One technical trap while you are here: everything leave1out() returns is on the log scale, untransformed. Exponentiating the point estimate and forgetting the interval limits produces a table whose intervals do not contain their own point estimates — and such a table does not look odd enough to make anyone stop.

Per-study influence diagnostics

Leave-one-out only shows how far the pooled value moves. The regression-diagnostics machinery behind influence() breaks the same question apart, giving each study its own set of numbers:

Being an outlier and being influential are two different things, and that is the distinction the whole set exists to make. A small, badly off study can be a plain outlier and still barely touch the pooled value, because its weight is tiny. A huge, entirely unremarkable trial can be nowhere near an outlier and still dominate the result, because its weight is large. Cook’s D is the usual summary precisely because it takes in both.

Nothing was flagged, and one trial still moved the pooled estimate

metafor’s influence() prints a column of asterisks marking the studies it considers influential. The rule is a union of four conditions; any one of them is enough:

Across these 13 trials the number of studies flagged is 0 — not one.

And that is not a lenient rule failing to fire: the largest Cook’s D reaches 50.4% of its threshold, the largest absolute DFFITS 60.1%, the largest hat 44.2%. Nothing comes close to any of the lines.

The Baujat plot: two questions in one picture

The Baujat plot compresses those diagnostics into one scatter plot whose two axes are the two questions:

• Horizontal = that study’s contribution to the heterogeneity, the squared standardised residual resid²/(τ² + vᵢ). Further right means it disagrees with the others. Across these 13 trials the horizontal values sum to 12.1.
• Vertical = the squared shift in the pooled estimate when the study is dropped, divided by the leave-one-out variance. Higher means the pooled estimate moves when it goes.

What each quadrant means:

On these data the upper-right quadrant holds 6 trials: Ferguson & Simes, 1949; Hart & Sutherland, 1977; Vandiviere et al, 1973; TPT Madras, 1980; Rosenthal et al, 1961; Comstock et al, 1976. The upper left holds only Frimodt-Moller et al, 1973 and the lower right only Comstock & Webster, 1969; the remaining 5 sit in the lower left.

The most influential trial: not broken, just sitting at one end of a moderator

All three tools point at the same study: Hart & Sutherland, 1977. It has the largest Cook’s D (0.229), the largest leave-one-out shift (9.0%), and the highest position in the upper-right quadrant of the Baujat plot.

What kind of trial is it? Every cell below comes straight out of dat.bcg, with no added assumptions:

This is not a problematic study. It was randomised, it is among the larger trials, it has plenty of events and a narrow interval. It is influential because its effect is a great deal stronger than the pooled one — its own RR is 0.237 against a pooled 0.489 — and because it carries enough weight for that difference to be heard.

So why is its effect so much stronger? The heterogeneity page has already given the answer: latitude. BCG performs markedly better at high latitudes, usually attributed to partial cross-immunity from exposure to environmental non-tuberculous mycobacteria nearer the equator. And Hart & Sutherland, 1977 sits at 52 degrees — the 2nd furthest from the equator in the set, and the largest of all the trials above 40 degrees.

Put latitude into the model and this becomes visible immediately:

That last row is the point of the whole page. In the intercept-only model the studentized residual of Hart & Sutherland, 1977 is -1.451, which reads as off-trend. Once latitude is in the model it shrinks to -0.664. In other words, as soon as the model knows the trial was run at 52 degrees, the log RR it observed (-1.442) is not far at all from what the model predicts (-1.262).

It was never an outlier. It is sitting at one end of a moderator.

What to do once you have found an influential study

An order you can actually follow:

1. First rule out a data error. Events, denominators, and direction — swapping the treatment and control arms is the classic one — and whether the effect measure is on the scale you think. This is the only step where finding a problem means fixing it immediately, and what gets fixed is the data, not the inclusion criteria.
2. Work out which dimension it differs on. Population, intensity of the intervention, what the control arm received, length of follow-up, how the outcome was defined, era, location, risk of bias. This step is done by reading the paper, not by computing.
3. Put that dimension into the model. With enough studies, meta-regression or a pre-specified subgroup analysis (at least ten studies per covariate — see the heterogeneity page); with too few, write it into the discussion as a possible source of heterogeneity.
4. Exclude only for a reason independent of the result. It is at high risk of bias is a reason. It moves the pooled estimate is not. Exclusion criteria belong in the protocol; if the decision was made after seeing the data, report both the included and the excluded analysis and say which one is primary.
5. Write the process down. Sensitivity analyses showed the result to be robust carries no information on its own. What to report is which perturbations were made, what each produced, and whether any of them changed the conclusion.

How to read the report

1. Was any influence analysis done at all? Without one, a reader cannot tell whether the pooled value rests on one or two studies. Box sizes on the forest plot are a crude substitute, but they show weight, not movement.
2. Does the leave-one-out report intervals, or only point estimates? A row of point estimates does not show whether the conclusion would flip.
3. Is the flag being treated as the only criterion? The data on this page are the counterexample: nothing flagged, and one trial still moving the estimate by close to a tenth.
4. Were studies excluded, and was the reason independent of the result? Pre-specified and post hoc exclusions are methodologically different things.
5. Is the post-exclusion result presented as primary or as a sensitivity analysis? Those two positions mean very different things.
6. Is there an explanation? Identifying the most influential trial and writing only the conclusion was unchanged after exclusion throws away the most informative observation in the analysis.

Run it yourself

library(metafor)
library(metadat)

data(dat.bcg, package = "metadat")
d <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg,
            data = dat.bcg, slab = paste(author, year, sep = ", "))
res <- rma(yi, vi, data = d, method = "REML")

# -- 1. Leave-one-out -------------------------------------------------------
# Everything comes back on the log scale: the CI limits need exp() as well.
# Exponentiating only the estimate gives a table whose intervals do not contain
# their own point estimates, and that table does not look odd enough to notice.
l1 <- leave1out(res)
data.frame(omitted = res$slab,
           rr = exp(l1$estimate),
           lb = exp(l1$ci.lb), ub = exp(l1$ci.ub),
           I2 = l1$I2)

# Percent change against the all-studies pooled value -- the column to look at
100 * (exp(l1$estimate) - exp(as.numeric(res$beta))) / exp(as.numeric(res$beta))

# -- 2. Per-study influence diagnostics -------------------------------------
infl <- influence(res)
infl                 # rstudent, dffits, cook.d, cov.r, tau2.del, QE.del, hat, weight
plot(infl)           # one small panel per statistic, cut-offs as dashed lines
infl$is.infl         # the asterisk column as a logical vector; all FALSE here

# The cut-offs are not magic numbers. p is the number of model parameters,
# m the degrees of freedom for Cook's D.
p <- infl$p; k <- res$k
3 * sqrt(p / (k - p))       # DFFITS
qchisq(0.5, df = infl$m)    # Cook's D
3 * p / k                   # hat

# -- 3. Baujat plot ---------------------------------------------------------
# It always draws, and returns the coordinates invisibly. For the coordinates
# alone, open a null device first with pdf(NULL).
bj <- baujat(res)
head(bj[order(-bj$y), ])    # highest on y = the ones that move the pooled value

# -- 4. Why the most influential trial is influential ------------------------
top <- which.max(infl$inf$cook.d)
dat.bcg[top, ]              # its covariates: latitude, allocation, 2x2 table

# Put the moderator in and its residual shrinks back towards zero
mr <- rma(yi, vi, mods = ~ ablat, data = d, method = "REML")
infl$inf$rstudent[top]              # intercept-only model
influence(mr)$inf$rstudent[top]     # model with latitude

# What dropping it actually costs: the slope and R^2 of the latitude model
rma(yi, vi, mods = ~ ablat, data = d, method = "REML", subset = -top)

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])
log_rr = np.log((tpos / (tpos + tneg)) / (cpos / (cpos + cneg)))
var    = 1 / tpos - 1 / (tpos + tneg) + 1 / cpos - 1 / (cpos + cneg)

# log_rr and var are computed just above (escalc)
k = len(log_rr)
full = combine_effects(log_rr, var, method_re="dl")
mu_full = full.mean_effect_re

# Leave-one-out has to be written as a loop; there is no ready-made function
for i in range(k):
    keep = np.arange(k) != i
    r = combine_effects(log_rr[keep], var[keep], method_re="dl")
    pct = 100 * (np.exp(r.mean_effect_re) - np.exp(mu_full)) / np.exp(mu_full)
    print(i, np.exp(r.mean_effect_re), round(pct, 2), r.i2)

# Cook's D, DFBETAS, hat and the Baujat coordinates have no counterpart in
# statsmodels. Reproducing them means implementing metafor's formulae yourself,
# refitted tau^2 and all. In practice, use R for this page.

Common misuses

Regenerating every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B7-07-influence-diagnostics.R

Watch next

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

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

Report a problem on this page