Pairwise Meta-Analysis: The Foundation of NMA

Author

Xiaoge Zhang, PhD

Published

September 18, 2026

Pooling direct evidence for one comparison, from the likelihood up

1 Executive Summary: One Comparison, Many Trials

A pairwise meta-analysis pools the results of several trials that all compare the same two treatments, say a new drug B against standard of care A. Each trial on its own is a noisy estimate of the treatment effect; pooled, they give one estimate with a narrower interval.

It is also the building block of every method in this library. A Standard Bayesian NMA is several pairwise meta-analyses estimated together, tied by the consistency assumption so that evidence on A vs B and A vs C also informs B vs C. Everything an NMA adds sits on top of the two questions answered here:

How should trials of different sizes be weighted?
What if the trials do not all estimate the same effect?

The worked example is the one NICE DSU TSD 2 uses for the same purpose, 22 trials of beta-blockers after myocardial infarction, rebuilt here from the published counts. Its posterior estimates and fit statistics are reproduced independently in Stan and set side by side with the published ones.

2 Theoretical Framework: From a 2×2 Table to a Pooled Estimate

2.1 The Effect Size in a Single Trial

For a binary outcome, where each patient either has the event of interest or does not, each trial reduces to a 2×2 table:

	Events	Non-events
Treatment B	\(a\)	\(b\)
Treatment A	\(c\)	\(d\)

Index the trials by \(i = 1, \dots, k\) (\(k = 22\) in the worked example below). In trial \(i\), arm B enrols \(n_B\) patients and arm A enrols \(n_A\). The event counts are random variables, \(a \sim \text{Binomial}(n_B, p_B)\) and \(c \sim \text{Binomial}(n_A, p_A)\), where \(p_B\) and \(p_A\) are the true event rates; the non-events are \(b = n_B - a\) and \(d = n_A - c\). The subscript \(i\) is dropped from these quantities to keep the formulas readable.

The odds ratio and its logarithm are

\[OR_i = \frac{a / b}{c / d} = \frac{ad}{bc}, \qquad y_i = \ln OR_i\]

Meta-analysis works on the log scale because the OR is bounded below by 0 and asymmetric around 1, whereas \(\ln OR\) is symmetric around 0 and close to Normal in reasonably sized trials.

The variance of \(y_i\) is

\[v_i = \frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}\]

Each arm contributes the reciprocals of its own two counts. Larger trials have larger counts, hence smaller \(v_i\).

Where the variance formula comes from

Step 1: the event rate in one arm. Take arm B, and write \(n = n_B\) and \(p = p_B\) for short. Since \(a \sim \text{Binomial}(n, p)\), the observed rate \(\hat p = a / n\) has variance

\[\text{Var}(\hat p) = \frac{p(1-p)}{n} \tag{1}\]

Step 2: from the rate to the log odds. The log OR is built not from \(\hat p\) but from its log odds, \(g(\hat p) = \ln\frac{\hat p}{1 - \hat p} = \ln\frac{a}{b}\). The delta method handles this change of scale. Near \(p\), the curve \(g\) is close to its tangent line:

\[g(\hat p) \approx g(p) + g'(p)\,(\hat p - p) \tag{2}\]

On the right of (2), \(g(p)\) and \(g'(p)\) are fixed numbers; only \(\hat p\) is random. Multiplying a random quantity by a constant multiplies its variance by the constant squared, and adding a constant leaves it unchanged, so

\[\text{Var}\big(g(\hat p)\big) \approx g'(p)^2 \, \text{Var}(\hat p) \tag{3}\]

Step 3: the slope of the log odds. Writing \(g(p) = \ln p - \ln(1 - p)\) and differentiating:

\[g'(p) = \frac{1}{p} + \frac{1}{1-p} = \frac{1}{p(1-p)} \tag{4}\]

Step 4: combine. Substituting (4) and (1) into (3):

\[\text{Var}\left(\ln\frac{a}{b}\right) \approx \frac{1}{p^2(1-p)^2} \cdot \frac{p(1-p)}{n} = \frac{1}{np(1-p)} \tag{5}\]

Step 5: back to the counts. Because \((1 - p) + p = 1\), the right-hand side of (5) splits into two fractions:

\[\frac{1}{np(1-p)} = \frac{(1 - p) + p}{np(1-p)} = \frac{1}{np} + \frac{1}{n(1-p)} \tag{6}\]

In (6), \(np\) and \(n(1-p)\) are the expected numbers of events and non-events. Replacing them with the observed counts \(a\) and \(b\) gives the variance for arm B:

\[\text{Var}\left(\ln\frac{a}{b}\right) \approx \frac{1}{a} + \frac{1}{b} \tag{7}\]

Steps 1 to 5 apply unchanged to arm A, with its counts \(c\) and \(d\):

\[\text{Var}\left(\ln\frac{c}{d}\right) \approx \frac{1}{c} + \frac{1}{d} \tag{8}\]

Step 6: two arms. The log OR is a difference of the two log odds, \(y_i = \ln\frac{a}{b} - \ln\frac{c}{d}\). For independent quantities, \(\text{Var}(X - Y) = \text{Var}(X) + \text{Var}(Y)\), and the two arms enrol different patients, so adding (7) and (8):

\[v_i = \text{Var}(y_i) \approx \frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d} \tag{9}\]

which is the formula above.

What it implies. The smallest cell dominates (9). A trial of 1,000 patients per arm with only 5 events on A still has \(v_i > 1/5\), so a large trial with a rare event can be less precise than a small trial with a common one. When a cell is zero, \(v_i\) is infinite, which is why a continuity correction, typically adding 0.5 to every cell, is applied. That correction is itself an assumption, and one that can move results when events are sparse.

2.2 Fixed Effect: Inverse-Variance Weighting

The fixed-effect model assumes every trial estimates one common effect \(\theta\), and differs from it only through sampling error:

\[y_i \sim N(\theta, v_i), \qquad i = 1, \dots, k\]

treating each \(v_i\) as known. Strictly, the variance of \(y_i\) depends on the true event rates \(p_B\) and \(p_A\), which are unknown; \(v_i\) replaces them with the observed rates, so it is an estimate. Treating it as a known constant is an approximation that works well with large counts, but with few events it makes both the weights and the pooled estimate unreliable, because \(y_i\) and \(v_i\) are computed from the same counts.

The log-likelihood is

\[\ell(\theta) = -\frac{1}{2} \sum_{i=1}^{k} \frac{(y_i - \theta)^2}{v_i} + \text{const}\]

Setting its derivative to zero gives the maximum likelihood estimate:

\[\sum_{i=1}^{k} \frac{y_i - \theta}{v_i} = 0 \quad \Rightarrow \quad \hat\theta_{FE} = \frac{\sum w_i y_i}{\sum w_i}, \qquad w_i = \frac{1}{v_i}\]

So the pooled estimate is a weighted average, and each trial’s weight is the inverse of its variance: precise trials count for more. The second derivative of \(\ell(\theta)\) is \(-\sum w_i\), which gives

\[\text{Var}(\hat\theta_{FE}) = \frac{1}{\sum w_i}\]

This is the same curvature argument as in The Core of Bayesian Estimation: every trial adds its weight to the curvature of the log-likelihood, so the interval narrows as trials accumulate.

2.3 Random Effects: When Trials Estimate Different Effects

Trials differ in populations, doses and settings, so the true effect may differ between them. The random-effects model lets each trial have its own true effect \(\theta_i\), drawn from a common distribution:

\[y_i \sim N(\theta_i, v_i), \qquad \theta_i \sim N(\mu, \tau^2)\]

Here \(\mu\) is the average effect across settings and \(\tau^2\) the between-study variance. \(v_i\) is the within-trial variance computed from \(a, b, c, d\): it describes how \(y_i\) varies around trial \(i\)’s own true effect \(\theta_i\). A single trial has only one \(\theta_i\), so it cannot show how true effects differ between trials; \(\tau^2\) is estimated from how much the trials scatter beyond what their own \(v_i\) can explain, as described in Section 2.4.2.

Integrating out \(\theta_i\) gives

\[y_i \sim N(\mu, \; v_i + \tau^2)\]

which is the fixed-effect model with every variance increased by \(\tau^2\). The same algebra therefore applies, with weights \(w_i^* = 1 / (v_i + \tau^2)\):

\[\hat\mu_{RE} = \frac{\sum w_i^* y_i}{\sum w_i^*}, \qquad \text{Var}(\hat\mu_{RE}) = \frac{1}{\sum w_i^*}\]

Adding the same \(\tau^2\) to every trial matters most for the trials whose own \(v_i\) is small. Large trials lose weight, small trials gain it, and the weights become more equal.

2.4 Measuring Heterogeneity: Q, τ² and I²

The random-effects model needs a value for \(\tau^2\), and the choice between the two models depends on how much the trials disagree. Three quantities describe that disagreement, and each answers a different question:

Quantity	Question it answers	Scale
Cochran’s \(Q\)	Do the trials disagree more than chance alone would explain?	A test statistic
\(\tau^2\)	How far apart are the true effects?	Variance of the log OR
\(I^2\)	What share of the observed scatter is real rather than chance?	A proportion, 0 to 100%

Note 1: Choosing between the two models

In order of what should carry the most weight:

Clinical judgement, first. Do the trials plausibly estimate one common effect: similar populations, doses, comparators, outcome definitions and follow-up? If they clearly do not, the fixed-effect assumption is hard to defend and the random-effects model is the honest default. A fixed-effect analysis then needs a positive argument, not just a favourable statistic.
The size of the heterogeneity, not its significance. Read \(\tau\) and the prediction interval (Section 2.5), and ask whether the spread they describe is large enough to change a decision. A statistically significant \(Q\) with a prediction interval that stays well away from no effect is a different situation from a non-significant \(Q\) with a prediction interval that crosses it.
Model fit. Compare the two models on residual deviance and DIC: lower is better, and the DIC already penalises the extra parameter the random-effects model brings. NICE DSU TSD 2 recommends consulting both and sets no numerical threshold for how large a difference counts; the figure of 3 to 5 points often quoted is a rule of thumb from the wider DIC literature, not from TSD 2. When the two models are close, TSD 2’s advice elsewhere is to decide on the data-generating process and the interpretability of the results rather than on the statistic.
How many trials there are. With few trials \(\tau\) is badly estimated, and a DerSimonian–Laird random-effects interval can be too narrow (Section 2.4.2). The remedies are the Hartung–Knapp–Sidik–Jonkman adjustment, REML, or a Bayesian model with an informative prior on \(\tau\). Switching to a fixed-effect model is not a remedy: it assumes the problem away.

Three things not to do: choose by whether \(Q\) is significant, choose by an \(I^2\) threshold alone, or choose the model that happens to give the narrower interval. NICE DSU TSD 2 sets out the fit-based comparison; the Standard Bayesian NMA guide lists the same judgement in a network setting.

2.4.1 Cochran’s Q: More Scatter Than Chance?

\(Q\) adds up how far each trial sits from the fixed-effect estimate:

\[Q = \sum_{i=1}^{k} w_i (y_i - \hat\theta_{FE})^2 = \sum_{i=1}^{k} \left( \frac{y_i - \hat\theta_{FE}}{\sqrt{v_i}} \right)^2\]

The second form shows what each term is: the distance of trial \(i\) from the pooled estimate, measured in that trial’s own standard errors, then squared. A small trial far from the pooled estimate adds little, because its standard error is large; a large trial the same distance away adds a lot.

If every trial estimates one common effect, each standardised distance behaves like a draw from \(N(0, 1)\), and the sum of \(k\) squared standard Normal draws follows a \(\chi^2\) distribution. One degree of freedom is lost because \(\hat\theta_{FE}\) is itself estimated from the same trials, so

\[Q \sim \chi^2_{k-1}, \qquad E(Q) = k - 1 \quad \text{under homogeneity}\]

That gives a test: a \(Q\) well above \(k - 1\) is unlikely if the trials share one effect. The test has two weaknesses that pull in opposite directions:

Few trials, low power. With five trials, \(Q\) has 4 degrees of freedom and real heterogeneity often fails to reach significance. A threshold of \(p < 0.10\) rather than \(0.05\) is commonly used for this reason, and a non-significant \(Q\) is not evidence that the trials agree.
Many large trials, too much power. With precise trials, differences too small to matter clinically still produce a significant \(Q\).

Either way, \(Q\) says whether heterogeneity is detectable, not how large it is.

2.4.2 τ²: How Far Apart Are the True Effects?

\(\tau^2\) is the variance of the true log ORs across settings, and \(\tau\) their standard deviation. It is the quantity with a direct clinical reading: under the random-effects model, about 95% of true effects lie within \(\mu \pm 1.96\,\tau\) on the log scale. On the OR scale that is a multiplicative range. With \(\tau = 0.3\), for example, true ORs span a factor of \(e^{1.96 \times 0.3} \approx 1.8\) either side of the average OR: an average OR of 2 would be consistent with true ORs from about 1.1 to 3.6 in different settings.

The DerSimonian–Laird estimator obtains \(\tau^2\) from \(Q\). Under the random-effects model, the expected value of \(Q\) rises in a straight line with \(\tau^2\):

\[E(Q) = (k - 1) + \tau^2 \left( \sum w_i - \frac{\sum w_i^2}{\sum w_i} \right)\]

The first term is the scatter chance alone produces; the second is the extra scatter from real differences. Setting the observed \(Q\) equal to its expected value and solving gives

\[\hat\tau^2_{DL} = \max\left(0, \; \frac{Q - (k - 1)}{\sum w_i - \sum w_i^2 / \sum w_i}\right)\]

This is a method-of-moments estimator: it matches one observed statistic to its expectation, with no likelihood involved. Two consequences follow:

Truncation at zero. By chance, \(Q\) can fall below \(k - 1\), which would give a negative variance. The estimate is then set to 0. An estimate of \(\hat\tau^2 = 0\) therefore means the scatter did not exceed chance, not that the true effects are identical.
One number, treated as known. The random-effects weights and interval use \(\hat\tau^2\) as if it were the true value. With few trials the estimate is very imprecise, and that imprecision is not carried into the pooled interval. Restricted maximum likelihood (REML) and Paule–Mandel are common alternative estimators; the Bayesian model in this guide instead gives \(\tau\) a posterior distribution.

Deriving E(Q) and the DerSimonian–Laird estimator

Write \(W = \sum w_i\) and \(\hat\theta = \sum w_i y_i / W\). Under the random-effects model, \(y_i = \mu + (\theta_i - \mu) + \varepsilon_i\) with \(\text{Var}(y_i) = v_i + \tau^2\), and the trials are independent.

Step 1: split Q. Adding and subtracting \(\mu\) inside the square, the cross term cancels because \(\sum w_i (y_i - \hat\theta) = 0\):

\[Q = \sum w_i (y_i - \mu)^2 - W (\hat\theta - \mu)^2\]

Step 2: the first term. Each \((y_i - \mu)^2\) has expectation \(v_i + \tau^2\), and \(w_i v_i = 1\):

\[E\left[\sum w_i (y_i - \mu)^2\right] = \sum w_i (v_i + \tau^2) = k + \tau^2 W\]

Step 3: the second term. \(\hat\theta\) is a weighted average of independent \(y_i\), so

\[\text{Var}(\hat\theta) = \frac{\sum w_i^2 (v_i + \tau^2)}{W^2} = \frac{W + \tau^2 \sum w_i^2}{W^2}\]

using \(w_i^2 v_i = w_i\). Since \(\hat\theta\) is unbiased for \(\mu\), \(E[W(\hat\theta - \mu)^2] = W\,\text{Var}(\hat\theta) = 1 + \tau^2 \sum w_i^2 / W\).

Step 4: combine.

\[E(Q) = (k + \tau^2 W) - \left(1 + \tau^2 \frac{\sum w_i^2}{W}\right) = (k - 1) + \tau^2 \left( W - \frac{\sum w_i^2}{W} \right)\]

With \(\tau^2 = 0\) this is the \(k - 1\) of the homogeneity test. Replacing \(E(Q)\) with the observed \(Q\) and solving for \(\tau^2\) gives the DerSimonian–Laird estimator.

2.4.3 I²: What Share of the Scatter Is Real?

\[I^2 = \max\left(0, \; \frac{Q - (k - 1)}{Q}\right)\]

\(Q\) is the total scatter, \(k - 1\) is the part chance would produce, and the difference is the part attributed to real differences. \(I^2\) is that difference as a share of the total. The Cochrane Handbook gives rough bands, deliberately overlapping: 0–40% might not be important, 30–60% may be moderate, 50–90% may be substantial, and 75–100% considerable.

The most common misreading is to treat \(I^2\) as a measure of how different the true effects are. It is not, because it compares \(\tau^2\) with the precision of the trials. The same \(I^2\) can be written as

\[I^2 = \frac{\hat\tau^2}{\hat\tau^2 + s^2}\]

where \(s^2\) is a typical within-trial variance. The same true spread therefore gives very different values of \(I^2\) depending on trial size. With \(\tau^2 = 0.09\) (\(\tau = 0.3\)):

Trials	Typical within-trial variance \(s^2\)	\(I^2\)
Small	0.20	\(0.09 / 0.29 \approx 31\%\)
Large	0.02	\(0.09 / 0.11 \approx 82\%\)

Large trials leave little room for chance, so almost all the scatter they show counts as real. Two further cautions apply:

I² is itself uncertain. With few trials its confidence interval commonly spans most of the 0–100% range, so a single value such as 59% should not be read to the nearest percentage point.
I² does not tell you what to expect in a new setting. For that, \(\tau\) and the prediction interval below are the relevant quantities.

Why I² equals τ² / (τ² + s²)

Write \(c = W - \sum w_i^2 / W\), the denominator of the DerSimonian–Laird estimator, so that \(\hat\tau^2 = (Q - (k - 1)) / c\). Higgins and Thompson (2002) define the typical within-trial variance as

\[s^2 = \frac{(k - 1)\,W}{W^2 - \sum w_i^2} = \frac{k - 1}{c}\]

Then, whenever \(Q > k - 1\),

\[\frac{\hat\tau^2}{\hat\tau^2 + s^2} = \frac{(Q - (k - 1)) / c}{(Q - (k - 1)) / c + (k - 1) / c} = \frac{Q - (k - 1)}{Q} = I^2\]

So \(I^2\) and \(\hat\tau^2_{DL}\) carry the same information from \(Q\); \(I^2\) simply rescales it by the precision of the trials.

2.5 The Prediction Interval

The confidence interval for \(\mu\) describes uncertainty about the average effect. A decision about a new population needs the range of effects to expect in one new setting, which also includes the spread \(\tau^2\):

\[\hat\mu_{RE} \pm 1.96 \sqrt{\hat\tau^2 + \text{Var}(\hat\mu_{RE})}\]

With substantial heterogeneity, the prediction interval can include no effect even when the confidence interval does not.

2.6 The Bayesian Version: The Same Model, Estimated Differently

The models above reduce each trial to two numbers, \(y_i\) and \(v_i\), and treat \(y_i\) as Normal with known variance. The Bayesian model used for NMA skips that reduction and works directly with the number of responders and patients in each arm, through a Binomial likelihood:

\[r_{iA} \sim \text{Binomial}(n_{iA}, p_{iA}), \qquad r_{iB} \sim \text{Binomial}(n_{iB}, p_{iB})\]

\[\text{logit}(p_{iA}) = \mu_i, \qquad \text{logit}(p_{iB}) = \mu_i + \delta_i, \qquad \delta_i \sim N(d, \tau^2)\]

This is the model of the Standard Bayesian NMA guide with a single comparison, and its unknowns are the ones already met above under different names: \(d\) is the average effect (\(\mu\) in the random-effects model), \(\delta_i\) is trial \(i\)’s own true effect (\(\theta_i\)), and \(\tau\) is the spread between trials. Only the baselines \(\mu_i\) are new, because reducing each trial to \(y_i\) had already removed them by taking a difference within the trial. Note that \(\mu_i\) here is a baseline log odds, not the average effect \(\mu\) of the random-effects model; the clash of letters is the convention used in the NMA literature.

What changes is not the model but what is treated as unknown. The frequentist approach treats \(\mu\) and \(\tau\) as fixed constants to be estimated, and \(\theta_i\) as a random draw that is never itself reported. The Bayesian approach treats all of them as unknown quantities with a joint posterior distribution:

\[ \begin{aligned} p(d, \tau, \delta, \mu \mid \text{counts}) \;\propto\; & \prod_{i=1}^{k} P(r_{iA}, r_{iB} \mid \mu_i, \delta_i) && \text{likelihood} \\[4pt] \times \; & \prod_{i=1}^{k} N(\delta_i \mid d, \tau^2) && \text{between-trial distribution} \\[4pt] \times \; & \; p(d)\; p(\tau) \prod_{i=1}^{k} p(\mu_i) && \text{priors} \end{aligned} \tag{10}\]

There is nothing to solve in (10); MCMC maps it by sampling, exactly as in The Core of Bayesian Estimation.

Counting the unknowns: with \(k\) trials there are \(2k + 2\) of them, a baseline and a true effect for each trial plus \(d\) and \(\tau\). The baselines are nuisance parameters, each identified by its own trial’s control arm and never pooled; the frequentist model removes them by taking a difference within each trial. The \(\delta_i\) are not free either: the between-trial distribution holds them together, which is what makes \(d\) and \(\tau\) estimable and what stops each \(\delta_i\) from simply reproducing that trial’s own log OR.

Two things follow, and both show up in the results below.

With τ held fixed, the two approaches give the same answer

Take the Normal approximation of the earlier sections, \(y_i \mid \theta_i \sim N(\theta_i, v_i)\) and \(\theta_i \sim N(\mu, \tau^2)\), so that marginally \(y_i \sim N(\mu, v_i + \tau^2)\). Hold \(\tau^2\) fixed and put a flat prior on \(\mu\), \(p(\mu) \propto 1\). The posterior is then proportional to the likelihood:

\[p(\mu \mid y) \;\propto\; \prod_{i=1}^{k} \exp\left(-\frac{(y_i - \mu)^2}{2(v_i + \tau^2)}\right) = \exp\left(-\frac{1}{2}\sum_{i=1}^{k} w_i^* (y_i - \mu)^2\right) \tag{11}\]

with \(w_i^* = 1/(v_i + \tau^2)\) as before. Expanding the sum in (11) and completing the square in \(\mu\) gives a Normal posterior:

\[\mu \mid y \;\sim\; N\left( \frac{\sum w_i^* y_i}{\sum w_i^*}, \; \frac{1}{\sum w_i^*} \right) \tag{12}\]

The posterior mean in (12) is exactly the random-effects estimate \(\hat\mu_{RE}\), and the posterior variance is exactly \(\text{Var}(\hat\mu_{RE})\). So the DerSimonian–Laird interval is the credible interval of this Bayesian model under a flat prior with \(\tau\) fixed at \(\hat\tau\).

The Bayesian model in this guide differs from that special case in two ways: it gives \(\tau\) a prior and integrates over it rather than fixing it, and it uses the Binomial counts rather than the Normal approximation. The first difference is what widens the interval in the results below.

Shrinkage: each trial’s own effect is revised

The posterior also reports each \(\delta_i\), which the frequentist model never estimates. Using the same Normal approximation, with \(\mu\) and \(\tau^2\) treated as known, the posterior for one trial’s true effect combines two sources: its own result \(y_i\), with precision \(1/v_i\), and the distribution the trials come from, with precision \(1/\tau^2\). Weighting each by its precision,

\[E(\theta_i \mid y_i) = B_i \, y_i + (1 - B_i)\, \mu, \qquad B_i = \frac{\tau^2}{\tau^2 + v_i} \tag{13}\]

Each trial’s estimate is pulled from its own result towards the average of the trials. How far it is pulled depends on the comparison in \(B_i\):

A small trial has a large \(v_i\), so \(B_i\) is small and its result is pulled a long way towards \(\mu\). An extreme result from a small trial is largely discounted.
A large trial has a small \(v_i\), so \(B_i\) is close to 1 and its result stands nearly as observed.
If \(\tau^2 \to 0\), every \(B_i \to 0\) and all trials collapse onto \(\mu\), which is the fixed-effect model.
If \(\tau^2\) is large, \(B_i \to 1\) and the trials are left alone, as if each were analysed separately.

In an NMA this is the mechanism behind a question EAGs often ask: how much of the estimate for a treatment supported by one small trial comes from that trial, and how much is borrowed from the rest of the network.

3 The Worked Example: Beta-Blockers After Myocardial Infarction

The example is the one NICE DSU TSD 2 uses throughout: 22 randomised trials of beta-blockers for preventing death after myocardial infarction, originally collated by Carlin (1992) and reproduced in TSD 2 Table 1. It is a pairwise comparison of two treatments with a binary outcome, which is exactly the setting of the models above.

Treatment A: control or placebo
Treatment B: beta-blocker
Outcome: death. The event is a harm, so an OR below 1 favours the beta-blocker.

Using a published example means the results can be checked: TSD 2 reports its own posterior estimates and fit statistics, and this page reproduces them independently in Stan.

View the code

library(knitr)

blocker <- data.frame(
    trial = 1:22,
    r_C = c(3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3, 40, 43, 39),
    n_C = c(39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147, 213, 122, 154, 134, 218, 364, 674),
    r_T = c(3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32, 27, 22),
    n_T = c(38, 114, 69, 1533, 355, 59, 945, 632, 278, 1916, 873, 263, 291, 858, 154, 207, 251, 151, 174, 209, 391, 680)
)

The 22 trials cover 20,290 patients and 1,811 deaths, and they differ in size by almost two orders of magnitude, from 77 patients to 3,837. That spread is what makes the weighting question concrete. The 22 rows of counts are listed in Section 8.

4 Implementation: From Scratch, `metafor` and Stan

4.1 From Scratch, Checked Against `metafor`

Every quantity from the theoretical framework, computed directly and then verified against the standard R package:

View the code

library(metafor)

fmt <- function(x) sprintf("%.2f", x)
ci_text <- function(est, se) sprintf("%.2f (%.2f, %.2f)", exp(est), exp(est - 1.96 * se), exp(est + 1.96 * se))

# Effect size and variance per trial, from the 2x2 cells
resp_B <- blocker$r_T
non_B <- blocker$n_T - blocker$r_T
resp_A <- blocker$r_C
non_A <- blocker$n_C - blocker$r_C

y <- log((resp_B / non_B) / (resp_A / non_A)) # log OR
v <- 1 / resp_B + 1 / non_B + 1 / resp_A + 1 / non_A
k <- length(y)

# Fixed effect
w <- 1 / v
lor_fe <- sum(w * y) / sum(w)
se_fe <- sqrt(1 / sum(w))

# Heterogeneity
Q <- sum(w * (y - lor_fe)^2)
Q_p <- pchisq(Q, df = k - 1, lower.tail = FALSE)
tau2 <- max(0, (Q - (k - 1)) / (sum(w) - sum(w^2) / sum(w))) # DerSimonian-Laird
I2 <- max(0, (Q - (k - 1)) / Q)

# Random effects, and the 95% prediction interval on the OR scale
w_re <- 1 / (v + tau2)
lor_re <- sum(w_re * y) / sum(w_re)
se_re <- sqrt(1 / sum(w_re))
pi_half <- 1.96 * sqrt(tau2 + se_re^2)
pi_lo <- exp(lor_re - pi_half)
pi_hi <- exp(lor_re + pi_half)

w_fe_pct <- 100 * w / sum(w)
w_re_pct <- 100 * w_re / sum(w_re)

# Same analysis through metafor, for verification
es <- escalc(measure = "OR", ai = r_T, n1i = n_T, ci = r_C, n2i = n_C, data = blocker)
fit_fe <- rma(yi, vi, data = es, method = "FE")
fit_re <- rma(yi, vi, data = es, method = "DL")
max_diff <- max(abs(
    c(lor_fe, se_fe, lor_re, se_re, tau2, 100 * I2) -
        c(coef(fit_fe), fit_fe$se, coef(fit_re), fit_re$se, fit_re$tau2, fit_re$I2)
))

kable(
    data.frame(
        c("Fixed effect", "Random effects (DerSimonian-Laird)"),
        c(ci_text(lor_fe, se_fe), ci_text(lor_re, se_re)),
        c("&ndash;", sprintf("%.2f to %.2f", pi_lo, pi_hi))
    ),
    col.names = c("Model", "Pooled OR (95% CI)", "95% prediction interval")
)

Model	Pooled OR (95% CI)	95% prediction interval
Fixed effect	0.77 (0.70, 0.85)	–
Random effects (DerSimonian-Laird)	0.78 (0.70, 0.87)	0.64 to 0.94

Heterogeneity is modest: \(Q\) = 23.26 on 21 degrees of freedom (p = 0.33), \(\hat\tau\) = 0.079 on the log OR scale, and \(I^2\) = 10%. The two models therefore give nearly the same answer, and the weights move only a little: the largest trial (trial 10, with 3,837 patients) carries 18% of the weight under fixed effects and 15% under random effects. metafor reproduces all six quantities above, the largest discrepancy being 6e-14.

The forest plot shows each trial with its 95% CI, the random-effects pooled estimate as a diamond, and the prediction interval as the dotted line through it. Each square is sized by the trial’s weight, and an arrow means the interval runs past the edge of the axis:

View the code

forest(
    fit_re,
    slab = paste("Trial", blocker$trial),
    atransf = exp,
    at = log(c(0.25, 0.5, 1, 2, 4)),
    addpred = TRUE,
    header = c("Trial", "OR [95% CI]")
)

4.2 The Bayesian Models in Stan

TSD 2 fits this example in WinBUGS with a fixed-effect and a random-effects program. The two Stan models below are the same specifications: Binomial likelihood, logit link, one baseline per trial, and TSD 2’s own vague priors, \(N(0, 100^2)\) for the baselines and the pooled effect and \(\text{Uniform}(0, 5)\) for the between-trial SD.

Stan model: blocker_fe.stan (fixed effect)

// Fixed effect model for the Blocker example, following NICE DSU TSD 2 Program 1(d)
data {
  int<lower=1> S;                        // number of trials
  array[S] int<lower=0> n_C;             // patients, control arm
  array[S] int<lower=0> r_C;             // deaths, control arm
  array[S] int<lower=0> n_T;             // patients, beta-blocker arm
  array[S] int<lower=0> r_T;             // deaths, beta-blocker arm
}

parameters {
  vector[S] mu;                          // trial baselines, log odds of death on control
  real d;                                // pooled log OR, beta-blocker vs control
}

model {
  mu ~ normal(0, 100);                   // TSD 2: vague N(0, 100^2)
  d ~ normal(0, 100);

  r_C ~ binomial_logit(n_C, mu);
  r_T ~ binomial_logit(n_T, mu + d);
}

generated quantities {
  real OR = exp(d);
  vector[S] fit_C;                       // expected deaths, control arm
  vector[S] fit_T;                       // expected deaths, treatment arm
  real dev = 0;                          // residual deviance, this draw

  for (i in 1:S) {
    fit_C[i] = n_C[i] * inv_logit(mu[i]);
    fit_T[i] = n_T[i] * inv_logit(mu[i] + d);
    dev += 2 * (r_C[i] * log(r_C[i] / fit_C[i])
                + (n_C[i] - r_C[i]) * log((n_C[i] - r_C[i]) / (n_C[i] - fit_C[i])));
    dev += 2 * (r_T[i] * log(r_T[i] / fit_T[i])
                + (n_T[i] - r_T[i]) * log((n_T[i] - r_T[i]) / (n_T[i] - fit_T[i])));
  }
}

Stan model: blocker_re.stan (random effects)

// Random effects model for the Blocker example, following NICE DSU TSD 2 Program 1(c)
data {
  int<lower=1> S;                        // number of trials
  array[S] int<lower=0> n_C;             // patients, control arm
  array[S] int<lower=0> r_C;             // deaths, control arm
  array[S] int<lower=0> n_T;             // patients, beta-blocker arm
  array[S] int<lower=0> r_T;             // deaths, beta-blocker arm
}

parameters {
  vector[S] mu;                          // trial baselines, log odds of death on control
  real d;                                // mean log OR, beta-blocker vs control
  real<lower=0, upper=5> sigma;          // between-trial SD; TSD 2: Uniform(0, 5)
  vector[S] z;                           // standardised trial deviations
}

transformed parameters {
  vector[S] delta = d + sigma * z;       // trial-specific log OR (non-centred)
}

model {
  mu ~ normal(0, 100);                   // TSD 2: vague N(0, 100^2)
  d ~ normal(0, 100);
  z ~ std_normal();                      // implies delta ~ normal(d, sigma)

  r_C ~ binomial_logit(n_C, mu);
  r_T ~ binomial_logit(n_T, mu + delta);
}

generated quantities {
  real OR = exp(d);
  real OR_new = exp(normal_rng(d, sigma)); // OR in a new trial setting
  vector[S] fit_C;                       // expected deaths, control arm
  vector[S] fit_T;                       // expected deaths, treatment arm
  real dev = 0;                          // residual deviance, this draw

  for (i in 1:S) {
    fit_C[i] = n_C[i] * inv_logit(mu[i]);
    fit_T[i] = n_T[i] * inv_logit(mu[i] + delta[i]);
    dev += 2 * (r_C[i] * log(r_C[i] / fit_C[i])
                + (n_C[i] - r_C[i]) * log((n_C[i] - r_C[i]) / (n_C[i] - fit_C[i])));
    dev += 2 * (r_T[i] * log(r_T[i] / fit_T[i])
                + (n_T[i] - r_T[i]) * log((n_T[i] - r_T[i]) / (n_T[i] - fit_T[i])));
  }
}

Both models also return the residual deviance at every iteration, which is what TSD 2’s fit statistics are built from:

Residual deviance \(\bar{D}_{res}\): the posterior mean of the deviance relative to a saturated model. Each independent data point should contribute about 1, so with 22 two-arm trials a well-fitting model gives roughly 44.
Effective number of parameters \(p_D\): the posterior mean deviance minus the deviance evaluated at the posterior mean of the fitted values.
DIC \(= \bar{D}_{res} + p_D\): fit plus a penalty for complexity, lower being better.

View the code

library(rstan)

stan_data <- list(
    S = nrow(blocker),
    n_C = blocker$n_C, r_C = blocker$r_C,
    n_T = blocker$n_T, r_T = blocker$r_T
)

fit_stan <- function(file) {
    stan(
        file = file, data = stan_data,
        chains = 4, iter = 6000, warmup = 2000,
        seed = 42, control = list(adapt_delta = 0.99), refresh = 0
    )
}

bayes_fe <- fit_stan("blocker_fe.stan")
bayes_re <- fit_stan("blocker_re.stan")

# TSD 2 fit statistics, from the deviance contributions the models return
fit_stats <- function(fit) {
1    D_bar <- mean(as.matrix(fit, pars = "dev"))
    fit_C <- colMeans(as.matrix(fit, pars = "fit_C"))
    fit_T <- colMeans(as.matrix(fit, pars = "fit_T"))
    dev_at <- function(r, n, f) 2 * (r * log(r / f) + (n - r) * log((n - r) / (n - f)))
2    D_hat <- sum(dev_at(blocker$r_C, blocker$n_C, fit_C)) + sum(dev_at(blocker$r_T, blocker$n_T, fit_T))
    c(D_res = D_bar, pD = D_bar - D_hat, DIC = D_bar + (D_bar - D_hat))
}

stats_fe <- fit_stats(bayes_fe)
stats_re <- fit_stats(bayes_re)

pars_fe <- summary(bayes_fe, pars = c("d", "OR"))$summary
pars_re <- summary(bayes_re, pars = c("d", "OR", "sigma", "OR_new"))$summary
posterior_tbl <- rbind(pars_fe, pars_re)[, c("mean", "sd", "2.5%", "50%", "97.5%", "n_eff", "Rhat")]
rownames(posterior_tbl) <- c(
    "Fixed effect: d (log OR)", "Fixed effect: OR",
    "Random effects: d (log OR)", "Random effects: OR",
    "Random effects: sigma (between-trial SD)", "Random effects: OR in a new trial"
)
kable(
    posterior_tbl,
    digits = c(3, 3, 3, 3, 3, 0, 3),
    col.names = c("Mean", "SD", "2.5%", "Median", "97.5%", "ESS", "R-hat")
)

n_div <- function(fit) sum(sapply(get_sampler_params(fit, inc_warmup = FALSE), function(x) sum(x[, "divergent__"])))

1: Posterior mean of the residual deviance, summed over all 44 trial arms
2: The same deviance evaluated once, at the posterior mean of the fitted values

	Mean	SD	2.5%	Median	97.5%	ESS	R-hat
Fixed effect: d (log OR)	-0.262	0.050	-0.360	-0.261	-0.165	11112	1.000
Fixed effect: OR	0.771	0.039	0.698	0.770	0.848	11088	1.000
Random effects: d (log OR)	-0.249	0.066	-0.375	-0.250	-0.116	9340	1.000
Random effects: OR	0.781	0.051	0.687	0.779	0.891	9262	1.000
Random effects: sigma (between-trial SD)	0.135	0.082	0.009	0.128	0.316	3358	1.001
Random effects: OR in a new trial	0.790	0.147	0.545	0.776	1.125	13979	1.000

Neither sampler reported any divergent transitions (0 and 0), the largest R-hat is 1.001, and the smallest effective sample size is 3,358.

4.3 The Same Quantities, Both Ways

The frequentist and Bayesian analyses estimate the same things. Setting them next to each other shows where they agree and where they cannot:

View the code

or_ci <- function(est, lo, hi) sprintf("%.2f (%.2f, %.2f)", exp(est), exp(lo), exp(hi))
post <- rstan::extract(bayes_re, pars = c("sigma", "OR_new"))

kable(
    data.frame(
        c(
            "Pooled OR, fixed effect",
            "Pooled OR, random effects",
            "Between-trial SD ($\\hat\\tau$ or $\\sigma$)",
            "95% prediction interval, OR",
            "Cochran's $Q$, $I^2$"
        ),
        c(
            or_ci(lor_fe, lor_fe - 1.96 * se_fe, lor_fe + 1.96 * se_fe),
            or_ci(lor_re, lor_re - 1.96 * se_re, lor_re + 1.96 * se_re),
            sprintf("%.3f", sqrt(tau2)),
            sprintf("%.2f to %.2f", pi_lo, pi_hi),
            sprintf("%.2f, %.0f%%", Q, 100 * I2)
        ),
        c(
            or_ci(pars_fe["d", "50%"], pars_fe["d", "2.5%"], pars_fe["d", "97.5%"]),
            or_ci(pars_re["d", "50%"], pars_re["d", "2.5%"], pars_re["d", "97.5%"]),
            sprintf("%.3f (%.3f, %.3f)", median(post$sigma), quantile(post$sigma, 0.025), quantile(post$sigma, 0.975)),
            sprintf("%.2f to %.2f", quantile(post$OR_new, 0.025), quantile(post$OR_new, 0.975)),
            "no counterpart"
        )
    ),
    align = "lcc",
    col.names = c("Quantity", "Frequentist (inverse variance, DL)", "Bayesian (Stan)")
)

Quantity	Frequentist (inverse variance, DL)	Bayesian (Stan)
Pooled OR, fixed effect	0.77 (0.70, 0.85)	0.77 (0.70, 0.85)
Pooled OR, random effects	0.78 (0.70, 0.87)	0.78 (0.69, 0.89)
Between-trial SD (\(\hat\tau\) or \(\sigma\))	0.079	0.128 (0.009, 0.316)
95% prediction interval, OR	0.64 to 0.94	0.54 to 1.12
Cochran’s \(Q\), \(I^2\)	23.26, 10%	no counterpart

The pooled effects agree closely, in both models and on both scales. That is not a coincidence: Section 2.6 shows that with \(\tau\) held fixed and a flat prior, the Bayesian posterior mean for the pooled effect is the inverse-variance weighted average, so the two can only differ through the treatment of \(\tau\) and through the Normal approximation the frequentist analysis makes for each trial.

The heterogeneity estimates do not agree, and the quantities that depend on them inherit the difference. Everything that follows from \(\tau\) is affected: the random-effects interval is a little wider in the Bayesian analysis, and the prediction interval noticeably so.

The last row is the one to note in an audit. \(Q\) and \(I^2\) are constructions of the Cochrane tradition and have no counterpart in this model; TSD 2 reports neither for this example, and defers heterogeneity to TSD 3. The two traditions summarise the same data differently: one reports heterogeneity statistics, the other model fit.

5 Results: Reproducing TSD 2

TSD 2 reports its own posterior summaries and fit statistics for this example in its Table 2. Placing them next to the Stan output shows whether the reconstruction holds:

View the code

# Published values, exactly as printed in TSD 2 Table 2
tsd2_fe <- c("-0.26", "0.050", "-0.36", "-0.16", "&ndash;", "46.8", "23.0", "69.8")
tsd2_re <- c("-0.25", "0.066", "-0.38", "-0.12", "0.13", "41.9", "28.1", "70.0")

ours_row <- function(pars, stats, sigma = NA) {
    c(
        sprintf("%.2f", pars["d", "50%"]), sprintf("%.3f", pars["d", "sd"]),
        sprintf("%.2f", pars["d", "2.5%"]), sprintf("%.2f", pars["d", "97.5%"]),
        if (is.na(sigma)) "&ndash;" else sprintf("%.2f", sigma),
        sprintf("%.1f", stats["D_res"]), sprintf("%.1f", stats["pD"]), sprintf("%.1f", stats["DIC"])
    )
}

kable(
    rbind(
        c("Fixed effect, TSD 2 Table 2", tsd2_fe),
        c("Fixed effect, this page", ours_row(pars_fe, stats_fe)),
        c("Random effects, TSD 2 Table 2", tsd2_re),
        c("Random effects, this page", ours_row(pars_re, stats_re, pars_re["sigma", "50%"]))
    ),
    align = "lcccccccc",
    escape = FALSE,
    col.names = c("", "$d$ (log OR)", "SD", "2.5%", "97.5%", "$\\sigma$", "$\\bar{D}_{res}$", "$p_D$", "DIC")
)

	\(d\) (log OR)	SD	2.5%	97.5%	\(\sigma\)	\(\bar{D}_{res}\)	\(p_D\)	DIC
Fixed effect, TSD 2 Table 2	-0.26	0.050	-0.36	-0.16	–	46.8	23.0	69.8
Fixed effect, this page	-0.26	0.050	-0.36	-0.16	–	46.8	23.2	70.0
Random effects, TSD 2 Table 2	-0.25	0.066	-0.38	-0.12	0.13	41.9	28.1	70.0
Random effects, this page	-0.25	0.066	-0.38	-0.12	0.13	41.9	28.3	70.2

Every quantity reproduces to within Monte Carlo error. The pooled log OR agrees to two decimal places in both models, the residual deviances match at 46.81 and 41.89, and \(p_D\) and DIC differ by about 0.2, which is the order of the simulation noise in these statistics. The reconstruction uses a different sampler (Stan rather than WinBUGS) and a different parameterisation of the random effects, so agreement at this level is a check on the model specification, not a repetition of the same arithmetic.

Three readings of this table matter for an HTA submission.

Fixed versus random effects, when the two agree. The random-effects model fits better (\(\bar{D}_{res}\) 41.89 against 46.81, for 44 data points) but uses more effective parameters (28.32 against 23.19), and the DIC ends up essentially tied: 70.22 against 70.00. TSD 2 draws the conclusion that “there is little to choose between the two models and the fixed effect model may be preferred since it is easier to interpret”. This is the situation described in Note 1: no threshold decides it, and the argument moves to interpretability.

What \(p_D\) shows about the parameter count. The fixed-effect model has 23 parameters, 22 baselines and one pooled effect, and \(p_D\) comes out at 23.19: with a vague prior, every parameter counts as one. The random-effects model nominally adds 22 trial-specific effects and \(\sigma\), which would be 46, but \(p_D\) is only 28.32. The 22 trial effects are held together by the between-trial distribution, so they are worth about 5.14 free parameters between them rather than 22. This is the shrinkage of Section 2.6 measured on real data.

What the pooled estimate means, and what it does not. The random-effects OR is 0.78 with a 95% credible interval of 0.69 to 0.89: on average across settings like these trials, beta-blockers reduce mortality. The prediction interval for a new setting is wider, 0.54 to 1.12, and its upper end sits above 1. TSD 2 does not report a prediction interval for this example, and neither do most submissions; with \(\sigma\) estimated at 0.13 the gap between the two intervals is real but moderate here. In an appraisal, the confidence interval answers “what is the average effect”, while the prediction interval answers “what might this treatment do in the population we are appraising”.

The two random-effects analyses also disagree about how much heterogeneity there is: DerSimonian-Laird returns \(\hat\tau\) = 0.08, the Bayesian model a posterior median of 0.13. With 22 trials that scatter little more than chance alone would explain, the data say little about \(\sigma\), so the vague \(\text{Uniform}(0, 5)\) prior keeps the posterior away from zero while the moment estimator is free to fall almost to it. That difference is why the frequentist prediction interval (0.64 to 0.94) is narrower than the Bayesian one, and it is the reason TSD 2 insists on inspecting the posterior for \(\sigma\) against its prior.

6 HTA Audit Points & Interview Q&A

Should I choose fixed or random effects based on I²?

No. I² is a supporting signal, not a switch, and a low I² with few trials is weak evidence of homogeneity, because \(Q\) has little power to detect heterogeneity when \(k\) is small. The judgement starts from whether the trials plausibly estimate one common effect; see Note 1 for the full order.

The beta-blocker trials give an I² of about 10%. Does that settle it?

It says that only about a tenth of the scatter between these trials looks like real difference rather than chance, which is consistent with the small \(\hat\tau\) and with the two models agreeing. It does not by itself justify a fixed-effect analysis: \(I^2\) is a share, not a size, and with 22 trials of very unequal precision the same true spread would read differently (Section 2.4.3). The prediction interval, which still reaches above 1, is the more decision-relevant summary.

Why is DerSimonian–Laird criticised when there are few trials?

It estimates \(\tau^2\) from one moment equation and then treats the estimate as fixed, so the confidence interval ignores the uncertainty in \(\tau^2\). With few trials that uncertainty is large, and the interval is too narrow. Common remedies are the Hartung–Knapp–Sidik–Jonkman adjustment, which rescales the variance and uses a \(t\) distribution, usually widening the interval, restricted maximum likelihood for \(\tau^2\), or a Bayesian model that carries the uncertainty in \(\tau\) through, as above.

When is a pairwise meta-analysis enough, and when do you need an NMA?

A pairwise meta-analysis is enough when the decision problem is one comparison and there are head-to-head trials of it. Once the decision involves several treatments, or the comparator of interest has never been tested directly against the new treatment, the evidence has to be connected through common comparators, which is what an NMA does.

7 Conclusion

Pairwise meta-analysis comes down to a weighted average with two decisions attached: how to weight the trials, and whether they share one effect. The fixed-effect model weights by precision alone; the random-effects model adds between-study variance, which evens out the weights, widens the interval and makes the prediction interval available. The Bayesian version works from the counts and carries the uncertainty in \(\tau\), and it is the same model that a Standard Bayesian NMA extends to a network of comparisons.

On the beta-blocker data all of this reproduces what TSD 2 reports, down to the residual deviance and the effective number of parameters. The two models are hard to separate on fit, which is the ordinary case rather than the exception, and the decision then rests on interpretability and on what the prediction interval says about a new population.

Methodological References: NICE DSU TSD 2 (Dias et al., 2011, updated 2016); Cochrane Handbook for Systematic Reviews of Interventions, Chapter 10; DerSimonian & Laird (1986); Higgins, Thompson & Spiegelhalter (2009) on prediction intervals.

8 Appendix: The Trial Data

Deaths and patients in each arm of the 22 trials, as given in NICE DSU TSD 2 Table 1 (Carlin, 1992).

Trial	Deaths / patients, control	Deaths / patients, beta-blocker	OR
1	3 / 39	3 / 38	1.03
2	14 / 116	7 / 114	0.48
3	11 / 93	5 / 69	0.58
4	127 / 1520	102 / 1533	0.78
5	27 / 365	28 / 355	1.07
6	6 / 52	4 / 59	0.56
7	152 / 939	98 / 945	0.60
8	48 / 471	60 / 632	0.92
9	37 / 282	25 / 278	0.65
10	188 / 1921	138 / 1916	0.72
11	52 / 583	64 / 873	0.81
12	47 / 266	45 / 263	0.96
13	16 / 293	9 / 291	0.55
14	45 / 883	57 / 858	1.33
15	31 / 147	25 / 154	0.73
16	38 / 213	33 / 207	0.87
17	12 / 122	28 / 251	1.15
18	6 / 154	8 / 151	1.38
19	3 / 134	6 / 174	1.56
20	40 / 218	32 / 209	0.80
21	43 / 364	27 / 391	0.55
22	39 / 674	22 / 680	0.54

1 Executive Summary: One Comparison, Many Trials

2 Theoretical Framework: From a 2×2 Table to a Pooled Estimate

2.1 The Effect Size in a Single Trial

2.2 Fixed Effect: Inverse-Variance Weighting

2.3 Random Effects: When Trials Estimate Different Effects

2.4 Measuring Heterogeneity: Q, τ² and I²

2.4.1 Cochran’s Q: More Scatter Than Chance?

2.4.2 τ²: How Far Apart Are the True Effects?

2.4.3 I²: What Share of the Scatter Is Real?

2.5 The Prediction Interval

2.6 The Bayesian Version: The Same Model, Estimated Differently

3 The Worked Example: Beta-Blockers After Myocardial Infarction

4 Implementation: From Scratch, metafor and Stan

4.1 From Scratch, Checked Against metafor

4.2 The Bayesian Models in Stan

4.3 The Same Quantities, Both Ways

5 Results: Reproducing TSD 2

6 HTA Audit Points & Interview Q&A

7 Conclusion

8 Appendix: The Trial Data

4 Implementation: From Scratch, `metafor` and Stan

4.1 From Scratch, Checked Against `metafor`