Hierarchical modeling applied to a meta-analysis

Meta-analysis is an increasingly popular and important process of summarizing and integrating the findings of research studies in a particular area. As a method for combining information from several parallel data sources, meta-analysis is closely connected to hierarchical modeling. In this section we consider a relatively simple application of hierarchical modeling to a meta-analysis in medicine. We consider another meta-analysis problem in the context of a decision problem later.

The data in our medical example are displayed in the first three columns of Table, which summarize mortality after myocardial infarction in 22 clinical trials, each consisting of two groups of heart attack patients randomly allocated to receive or not receive beta-blockers (a family of drugs that affect the central nervous system and can relax the heart muscles). Mortality varies from 3% to 21% across the studies, most of which show a modest, though not “statistically significant,” benefit from the use of beta-blockers. The aim of a meta-analysis is to provide a combined analysis of the studies that indicates the overall strength of the evidence for a beneficial effect of the treatment under study. Before proceeding to a formal meta-analysis, it is important to apply rigorous criteria in determining which studies are included.

---

Defining a parameter for each study

In the beta-blocker example, the meta-analysis involves data in the form of several $2 \times 2$ tables. If clinical trial $j$ involves the use of $n_{0j}$​ subjects in the control group and $n_{1j}$ in the treatment group, giving rise to $y_{0j}$ and $y_{1j}$ deaths in control and treatment groups, respectively, then the usual sampling model involves two independent binomial distributions with probabilities of death $p_{0j}$ and $p_{1j}$, respectively.

Estimands of interest include

• the difference in probabilities: $p_{1j} – p_{0j}$,
• the risk ratio: $\dfrac{p_{1j}}{p_{0j}}$​​,
• and the odds ratio: $\rho_j = \frac{p_{1j}/(1 – p_{1j})}{p_{0j}/(1 – p_{0j})}.$

For reasons of interpretability and approximate normality, we focus on the log odds ratio, and we label it $\theta_j = \log \rho_j.$

---

A normal approximation to the likelihood

We can do a relatively simple Bayesian meta-analysis using the normal-theory results from the earlier sections if we reduce each study to an approximate normal likelihood for $\theta_j$​. One standard approximation is based on empirical logits. For each study $j$, define the estimate

$y_j = \log \frac{y_{1j}}{n_{1j} – y_{1j}} – \log \frac{y_{0j}}{n_{0j} – y_{0j}}$ (5.23)

with approximate sampling variance

$\sigma_j^2 = \frac{1}{y_{1j}} + \frac{1}{n_{1j} – y_{1j}} + \frac{1}{y_{0j}} + \frac{1}{n_{0j} – y_{0j}}$ (5.24)

We reuse the notation $y_j$​ and $\sigma_j^2$​ so that we can plug directly into the hierarchical normal model from the previous section. In practice one often adds 0.5 to each cell to stabilize the logit, but with the generally large sample sizes here this is not crucial.

The estimated log-odds ratios $y_j$​ and their estimated standard errors $\sigma_j$ are shown in the fourth and fifth columns of Table. We will use a hierarchical Bayesian analysis to combine information from the 22 studies and obtain improved estimates of each $\theta_j$​, along with estimates of the mean and variance of the study effects.

---

Goals of inference in meta-analysis

When doing meta-analysis, one should be clear about the estimand. Under the assumption that the studies are broadly comparable, there are three natural perspectives:

1. Complete pooling view: the studies are regarded as identical replications, so all individuals are effectively sampled from one common population; then there is a single common effect.
2. No pooling view: the studies are so different that one study tells us nothing about another; then we estimate each effect separately.
3. Exchangeable view: the studies are not identical but can be regarded as exchangeable draws from some population of study effects; there may be real differences, but not ones we can predict ahead of time in favor of a particular study.

As in earlier sections, this third view is the most flexible and is what the Bayesian hierarchical model formalizes.

Exchangeability itself does not force the population distribution to be normal, but a normal model is convenient and usually adequate—and we can check it.

In this setting, useful estimands are:

• the population mean effect (overall average log-odds reduction),
• the effect in any of the observed studies,
• the effect in a future, comparable study.

---

A hierarchical normal model

Since each study provides an estimate $y_j$​ of $\theta_j$ with (approximately) known sampling variance $\sigma_j^2$​, we can write the first level as $y_j \mid \theta_j, \sigma_j \sim N(\theta_j, \sigma_j^2), \qquad j = 1, \dots, J$,

where here $J = 22$. At the second level we place an exchangeable normal prior on the study-specific effects:

$\theta_j \mid \mu, \tau \sim N(\mu, \tau^2)$

where $\mu$ is the population mean log-odds ratio and $\tau$ is the between-study standard deviation of the log-odds ratios. Finally, we assign a (locally) uniform prior for $\mu$ and likewise a diffuse (locally uniform) prior for $\tau$, just as in the SAT coaching example; the data from 22 trials are informative enough about $\mu$.

So the model is the same as before, just with “schools” replaced by “trials” and “effects on SAT-V” replaced by “log odds ratios of mortality.”

---

Results of the analysis and comparison to simpler methods

Applying the same computational strategy as in the previous section (simulate $\tau$, then $\mu$, then the $\theta_j$​’s), we obtain a marginal posterior for $\tau$ that peaks at a nonzero value, but still puts nontrivial mass near 0. In fact, the posterior density at $\tau = 0$ is only about 25% lower than at the mode—so the data say “the trials are fairly homogeneous, but maybe not exactly identical.”

Posterior quantiles for the 22 study effects $\theta_j$​ (on the log-odds scale) are shown in Table (last columns in the original text). Because the posterior for $\tau$ is concentrated on values smaller than or comparable to the sampling standard errors $\sigma_j$​, we get considerable shrinkage toward the overall mean, with the biggest shrinkage for the noisiest studies (e.g. studies 1, 6, 18).

Above Table summarizes the posterior for the hyperparameters and for a future study:

• for the mean: $\mu \text{ (median)} \approx -0.25$ with 95% posterior interval $[-0.37, -0.11]$;
• for the between-study s.d.: $\tau \text{ (median)} \approx 0.13,$ with most mass below about 0.3;
• for a predicted effect in a new study, $\tilde{\theta}_j$​ has a wider interval, reflecting both estimation uncertainty and genuine study-to-study variation.

Because $\tau$ is not literally 0, the hierarchical Bayesian intervals are wider than the intervals you’d get by complete pooling (i.e. fixing $\tau = 0$). For example, converting $\mu$’s 95% interval $[-0.37, -0.11]$ to the odds-ratio scale by exponentiating gives about $[0.69, 0.90]$. Complete pooling would have given about $[0.70, 0.85]$, which was criticized as “unusually narrow.” The hierarchical model explains why: complete pooling is pretending all 22 trials are identical—i.e. it is fixing $\tau \equiv 0$—whereas the data actually allow $\tau$ to be as high as 0.3.

Finally, note that if the real interest is “what might happen in a population like these, but not exactly the same as any of the 22?”, then the relevant quantity is the predictive study effect $\tilde{\theta}_j$​, not just the mean $\mu$. That predictive distribution is wider still, and in this case there is a bit more than 10% posterior probability that $\tilde{\theta}_j > 0,$ i.e. that in a new, exchangeable study the treatment could even increase mortality.