Bayesian Hierarchical Models in Clinical Research

Introduction

Clinical research frequently produces data that are naturally organized into groups, levels, or clusters.

Patients may be treated at different clinical centers. Patients may belong to different geographic regions. Multiple clinical trials may evaluate related treatments. Repeated observations may be collected from the same patient. Subgroups may have related but non-identical treatment effects.

These structures create statistical dependence and information-sharing opportunities that are not adequately represented by treating every observation as completely independent.

A hierarchical model, also called a multilevel model, explicitly represents these levels of variation.

A Bayesian hierarchical model goes one step further by placing probability distributions on parameters at different levels of the hierarchy.

Why Hierarchical Models Matter in Clinical Research

Consider a randomized clinical trial conducted at 20 hospitals.

The treatment effect may differ somewhat from hospital to hospital because of differences in:

• Patient characteristics
• Clinical practice
• Patient recruitment
• Background therapy
• Investigator experience
• Measurement procedures
• Chance variation

Let:

A model that assumes one identical treatment effect for every center is a complete-pooling approach.

At the opposite extreme, estimating every center effect completely independently is a no-pooling approach.

A hierarchical model provides a middle ground.

Each center has its own treatment effect \(\theta_j\), but those effects are modeled as arising from a common population distribution with mean \(\mu\) and between-center standard deviation \(\tau\).

The Three Pooling Strategies

The distinction between complete pooling, no pooling, and partial pooling is fundamental to understanding hierarchical models.

Complete Pooling

Suppose there are \(J\) clinical centers. A complete-pooling model assumes:

There is only one treatment-effect parameter.

This is statistically efficient if the treatment effect really is identical across centers.

However, it does not allow genuine between-center heterogeneity.

No Pooling

At the other extreme, each center can receive its own independent parameter:

Each effect is estimated using only the observations from that center.

This allows maximum heterogeneity but can produce unstable estimates when some centers have small sample sizes.

Partial Pooling

A hierarchical model allows center-specific effects while connecting them through a common distribution:

Here:

• \(\theta_j\) is the effect for center \(j\).
• \(\mu\) is the overall mean effect.
• \(\tau\) is the between-center standard deviation.

The center-specific effects are neither completely identical nor completely independent.

They are allowed to differ, but the model learns how much variation among them is supported by the data.

The Hierarchical Structure

A simple hierarchical model has two levels.

At the patient level:

where \(Y_{ij}\) is the outcome for patient \(i\) in group \(j\).

At the group level:

The population-level parameters receive their own priors:

The resulting model contains several sources of uncertainty simultaneously.

Exchangeability

Hierarchical Bayesian models are closely connected to the concept of exchangeability.

Suppose treatment effects from several studies are represented by:

If the studies are considered exchangeable, their labels do not determine their prior treatment-effect distributions.

A simple representation is:

This does not mean that all studies have the same treatment effect.

It means that before observing the study-specific data, the effects are considered to arise from the same population distribution.

Between-Group Heterogeneity

The parameter \(\tau\) controls the amount of between-group heterogeneity.

If:

the group-specific effects are concentrated near \(\mu\).

The groups behave approximately like a common-effect population.

If \(\tau\) is large, the model permits substantial differences between groups.

Partial Pooling and Shrinkage

One of the most important consequences of hierarchical modeling is shrinkage.

Suppose center \(j\) has an observed treatment effect that is unusually large. A hierarchical model does not necessarily accept that estimate at face value.

Instead, it asks: How much evidence does this center provide for having a genuinely unusual effect, given both its own data and the distribution of effects across all centers?

An estimate based on a small sample will generally be pulled more strongly toward the population mean.

An estimate based on a large sample will generally be pulled less.

A Simple Normal-Normal Example

Suppose study \(j\) produces an approximately normal estimate:

Assume the hierarchical prior:

If \(\mu\) and \(\tau\) are treated as known, the posterior distribution is also normal.

The posterior mean is:

where:

The weight \(w_j\) determines how strongly the study-specific estimate influences the posterior.

Understanding the Shrinkage Weight

Consider the two extremes.

If the study is very precise:

then:

and the posterior estimate approaches the observed estimate.

If the study is very imprecise:

then:

and the posterior estimate approaches the population mean.

Worked Shrinkage Example

Suppose the overall treatment-effect distribution has:

Consider a small study with:

and standard error:

The shrinkage weight is:

Therefore the posterior mean is:

The raw estimate was 0.50, but the hierarchical estimate is 0.26.

The reason is that the study provides relatively weak information about its true effect.

What If the Study Is Large?

Suppose the same observed effect occurs in a much larger study with:

Then:

The posterior mean becomes:

The estimate is still shrunk, but much less strongly.

This illustrates why hierarchical models do not simply force small studies to look like large studies.

Instead, the amount of borrowing depends on the information contained in the individual group's data and the estimated heterogeneity across groups.

Bayesian Hierarchical Models for Multi-Center Trials

Suppose a randomized clinical trial includes \(J\) centers. Let:

A simple continuous-outcome model could be:

where:

Here \(T_{ij}\) indicates treatment assignment.

The center-specific treatment effects can then be modeled hierarchically:

The parameter \(\beta\) represents the population-average treatment effect, while \(\tau_\beta\) represents heterogeneity in treatment effects across centers.

Random Intercepts

Hierarchical models do not require treatment effects to vary by group.

A common model allows centers to have different baseline outcomes while assuming a common treatment effect.

with:

The center-specific intercepts account for baseline differences among centers.

This is commonly called a random-intercept model.

Random Slopes

If the treatment effect itself varies across centers, a random slope can be introduced:

with:

The covariance matrix \(\Sigma\) can represent variation in baseline outcomes, variation in treatment effects, and correlation between the two.

Random Intercept vs. Random Slope

Binary Outcomes

Hierarchical models are especially useful when clinical endpoints are binary.

Suppose:

A hierarchical logistic model can be written as:

where:

This allows center-specific baseline response probabilities.

Hierarchical Treatment Effects for Binary Outcomes

To allow treatment effects to vary:

and:

The treatment-specific odds ratio for center \(j\) is:

The hierarchical model therefore produces a posterior distribution for each center's treatment effect as well as the overall distribution of treatment effects.

Hierarchical Models for Continuous Outcomes

For a continuous endpoint such as change from baseline:

A simple model is:

with:

This can represent center-specific baseline levels while estimating a common treatment effect.

If treatment effects vary:

The same conceptual framework applies.

Hierarchical Models for Count Outcomes

Clinical research also frequently involves counts, such as:

• Number of exacerbations
• Number of hospitalizations
• Number of adverse events
• Number of disease episodes

A Poisson hierarchical model could be:

with a hierarchical center effect:

For overdispersed outcomes, a negative-binomial likelihood may be more appropriate.

Hierarchical Models for Time-to-Event Outcomes

Hierarchical structures can also be incorporated into survival models.

For example, a proportional hazards model could be written:

where the baseline hazard can vary by center.

Alternatively, treatment effects can vary hierarchically:

The resulting model can estimate both an overall treatment effect and between-center heterogeneity.

Hierarchical Models for Repeated Measures

Repeated-measures data have a natural hierarchical structure because observations are nested within patients.

For patient \(i\) at time \(t\):

A random intercept model might use:

This allows patients to have different baseline values.

A random slope model could additionally allow individual trajectories to vary:

with:

This is particularly useful when patients have different baseline disease severity and different rates of change over time.

Hierarchical Models for Subgroup Effects

Clinical trials often evaluate treatment effects in subgroups such as:

• Age categories
• Sex
• Disease severity
• Biomarker-defined populations
• Geographic regions
• Prior treatment status

A separate analysis for every subgroup can produce unstable estimates, especially when some subgroups are small.

A hierarchical model can allow subgroup-specific effects:

where \(k\) indexes subgroups.

The model estimates each subgroup while allowing the subgroups to share information.

Why Hierarchical Subgroup Models Can Be Preferable

Suppose a trial has ten subgroups. If each subgroup is analyzed independently, a small subgroup may produce an extreme estimate simply because its sample size is limited.

A hierarchical model provides regularization.

A subgroup with little information is pulled toward the overall effect, whereas a subgroup with strong evidence for a genuinely different effect can remain distinct.

Historical Borrowing

One of the most important clinical applications of Bayesian hierarchical models is borrowing information from historical or external data.

Suppose a current clinical trial has a control group with a relatively small sample size. Historical trials may contain substantial information about the control response rate.

Instead of treating the historical response rate as either completely irrelevant or perfectly interchangeable with the current control group, a hierarchical model can allow the historical and current effects to be related but not identical.

A Hierarchical Historical-Control Model

Suppose study \(j\) has control response probability \(p_j\). A convenient transformed scale is the logit:

The study-specific effects can be modeled as:

The current study is treated as another member of the distribution, but the model estimates the amount of between-study heterogeneity.

If the historical studies are highly compatible with the current study, more borrowing occurs.

If they are inconsistent, the estimated heterogeneity can reduce the amount of borrowing.

Borrowing Is Controlled by Heterogeneity

The central parameter is again:

Small \(\tau\) implies that the studies are similar.

Large \(\tau\) implies that the studies may differ substantially.

Thus the hierarchical model can adapt the amount of borrowing according to the estimated compatibility among studies.

Complete Pooling vs. Hierarchical Historical Borrowing

Meta-Analysis as a Hierarchical Model

Bayesian hierarchical models provide a natural framework for meta-analysis.

Suppose study \(k\) estimates a treatment effect \(\hat\theta_k\) with standard error \(s_k\). A simple likelihood is:

The study-specific effects follow:

The overall effect is \(\mu\), while \(\tau\) represents between-study heterogeneity.

This is the Bayesian analogue of a random-effects meta-analysis.

Hierarchical Models and Clinical Trial Estimands

A sophisticated hierarchical model does not replace the need to define the clinical estimand.

Before specifying the model, investigators should determine what treatment effect is actually being estimated.

For example:

• Population-average treatment effect
• Center-specific treatment effect
• Subgroup treatment effect
• Average treatment effect across a distribution of populations
• Probability that treatment benefit exceeds a clinical threshold

A hierarchical model can estimate several of these quantities, but they are not automatically interchangeable.

Population-Average vs. Group-Specific Effects

Suppose:

The parameter \(\mu\) describes the mean of the distribution of group-specific effects.

It does not necessarily equal every group's treatment effect.

For a randomly selected group from the modeled population:

But individual groups may have:

This distinction becomes particularly important when interpreting subgroup effects.

Conditional vs. Marginal Effects

In nonlinear models such as logistic regression, the population-average (marginal) effect is not necessarily equal to the average of the conditional group-specific effects after transformation.

For example, if:

then \(e^\beta\) is a conditional odds ratio on the model's log-odds scale.

A marginal population-level treatment effect may require integrating over the distribution of the random effects.

Therefore, analysts should distinguish between:

• Conditional effects
• Population-average effects
• Subject-specific effects
• Cluster-specific effects

Bayesian Posterior Distribution in a Hierarchical Model

Let the observed data be \(D\), group-specific parameters be \(\theta=(\theta_1,\ldots,\theta_K)\), and hyperparameters be \(\phi=(\mu,\tau)\). The posterior is:

The likelihood describes the patient-level observations.

The hierarchical distribution describes how group-specific parameters relate to the population-level parameters.

The hyperprior describes uncertainty about the population-level parameters.

Hyperparameters and Hyperpriors

Parameters that govern the distribution of other parameters are often called hyperparameters.

In:

\(\mu\) and \(\tau\) are hyperparameters for the distribution of \(\theta_j\).

They can themselves have prior distributions. For example:

The exact scale should be chosen based on the scale and plausible range of the clinical parameter.

Why Priors on Heterogeneity Matter

The between-group standard deviation \(\tau\) is often difficult to estimate, especially when there are relatively few groups.

Consequently, the prior on \(\tau\) can have substantial influence.

An overly permissive prior may imply implausibly large heterogeneity.

An overly restrictive prior may force the model toward near-complete pooling.

Prior Predictive Checking

Before fitting a hierarchical model, investigators can simulate outcomes from the prior predictive distribution.

Conceptually:

The resulting simulated datasets should be examined for plausibility.

If the model frequently generates impossible or clinically absurd outcomes, the prior specification should be reconsidered.

Posterior Predictive Checking

After fitting the model, the same idea can be applied using the posterior.

Draw parameters from:

and generate replicated datasets:

The replicated outcomes can be compared with the observed data.

Useful comparisons may include:

• Overall event rates
• Group-specific event rates
• Distribution of outcomes
• Extreme observations
• Between-center variability
• Within-patient variability
• Correlation structures

Posterior Prediction for a New Center

One advantage of hierarchical models is that they naturally support prediction for a new group.

Suppose a future center has not yet contributed data. Its treatment effect can be predicted from the population distribution:

When uncertainty in \(\mu\) and \(\tau\) is included:

This produces a posterior predictive distribution for a future center.

Clinical Interpretation of a New-Center Prediction

Suppose a trial estimates:

This could be interpreted as: Given the hierarchical model and observed trial data, the posterior probability that a future exchangeable center has a positive treatment effect is 87%.

This is different from:

because the latter concerns the mean treatment effect across the modeled population of centers.

Hierarchical Models for Rare Diseases

Rare disease trials often have small sample sizes.

This creates an especially challenging combination:

• Few patients
• Few events
• Multiple clinical centers
• Potentially substantial patient heterogeneity
• Limited opportunities for independent replication

Hierarchical models can potentially improve estimation by allowing related populations to share information.

Examples include:

• Multiple rare disease trials
• Related disease subtypes
• Historical controls
• Multiple endpoints
• Age-defined populations

However, the validity of the borrowing depends critically on the plausibility of the exchangeability assumptions.

Borrowing Across Related Indications

Suppose a therapy is studied in several related diseases. Let:

A hierarchical model could specify:

This permits the indications to have related effects without requiring them to be identical.

If one indication has sparse data, its estimate may borrow information from the other indications.

Dynamic Borrowing

Hierarchical models are sometimes described as providing dynamic borrowing.

The amount of borrowing depends on how compatible the data are with the assumed common distribution.

When groups appear similar, the model can borrow more strongly.

When groups differ substantially, the estimated heterogeneity increases and borrowing becomes weaker.

This is one of the most attractive features of hierarchical modeling compared with simply combining all observations.

But Dynamic Borrowing Is Not Magic

The model can only adapt to differences that are represented by the hierarchical structure.

If there is a systematic difference between groups that is not modeled, the hierarchical model can still borrow inappropriately.

For example, historical studies may differ systematically because of:

• Changes in standard of care
• Different eligibility criteria
• Different diagnostic definitions
• Different follow-up periods
• Changes in supportive therapy
• Different endpoint assessment procedures

A hierarchical distribution does not automatically correct these problems.

Adding Covariates to a Hierarchical Model

Hierarchical models can incorporate patient-level and group-level covariates.

Suppose treatment effect varies with a patient-level biomarker \(X_{ij}\):

The interaction parameter \(\delta\) describes how treatment effect changes with the biomarker.

The center-specific intercept may still be modeled hierarchically:

The treatment effect itself could also vary hierarchically.

Cross-Level Interactions

Hierarchical models can also represent interactions between patient-level and group-level characteristics.

For example, suppose \(X_{ij}\) is a patient biomarker and \(Z_j\) is a center-level characteristic. A model could include:

The interaction term \(\eta\) allows the association between the patient-level biomarker and outcome to depend on the center-level characteristic.

These structures can become complex quickly, so prespecification and careful model diagnostics are important.

Multilevel Models vs. Fixed-Effects Models

A center can be included in a model using fixed effects or random effects.

A fixed-effects approach might include center indicators:

A hierarchical model instead treats center effects as draws from a population distribution.

These approaches answer somewhat different questions and rely on different assumptions.

Hierarchical Bayesian Models and Multiple Comparisons

Hierarchical models are sometimes described as a solution to multiplicity. That statement requires care.

Suppose 20 subgroup treatment effects are estimated. A hierarchical model can regularize extreme estimates by shrinking them toward the population distribution.

This can reduce overinterpretation of noisy subgroup estimates.

However, shrinkage is not identical to a conventional frequentist multiplicity adjustment.

The model addresses the problem through a probability distribution and regularization structure rather than automatically controlling a frequentist family-wise error rate.

Hierarchical Models and Multiplicity in Subgroup Analysis

Suppose there are \(K\) subgroup effects:

A hierarchical prior:

encourages estimates to remain compatible with the overall distribution unless there is sufficient subgroup-specific information to support divergence.

This can be particularly useful for exploratory subgroup analyses.

However, confirmatory subgroup claims may require additional prespecification and decision criteria.

Non-Centered Parameterization

Hierarchical Bayesian models can create challenging posterior geometries, particularly when the number of groups is small or the between-group standard deviation is near zero.

A centered parameterization is:

where:

This representation is often called a non-centered parameterization.

It can improve sampling efficiency in some hierarchical models.

The appropriate parameterization depends on the data and model.

MCMC Estimation

Many hierarchical Bayesian models do not have closed-form posterior distributions.

They are commonly estimated using Markov chain Monte Carlo methods such as:

• Hamiltonian Monte Carlo
• No-U-Turn Sampling
• Other MCMC algorithms

The goal is to obtain draws from the posterior distribution:

Posterior summaries can then be calculated from these draws.

Convergence Diagnostics

A hierarchical Bayesian analysis should assess whether the computational algorithm has adequately explored the posterior distribution.

Common diagnostics include:

• Trace plots
• Effective sample size
• \(\hat R\) convergence diagnostics
• Monte Carlo standard error
• Divergent transitions for HMC-based methods
• Maximum treedepth warnings
• Energy diagnostics where appropriate

Why Divergent Transitions Matter

In Hamiltonian Monte Carlo, divergent transitions can indicate that the algorithm is having difficulty exploring parts of the posterior geometry.

They are particularly common in hierarchical models with difficult variance-component structures.

Possible remedies can include:

• Improving parameterization
• Using more appropriate priors
• Rescaling predictors
• Increasing the HMC target acceptance rate
• Reconsidering model complexity

Simply ignoring divergent transitions is not good practice.

Centering and Scaling Predictors

Hierarchical regression models often benefit from sensible scaling of continuous covariates.

Suppose age is entered as:

A centered variable might be:

This can make the intercept easier to interpret and can improve computational geometry.

Scaling can be especially useful when predictors have very different magnitudes.

Bayesian Hierarchical Model in R

A simple hierarchical normal model can be fit using a Bayesian modeling package such as brms.

Suppose the data contain:

• y: continuous outcome
• treatment: treatment indicator
• center: clinical center

A model with center-specific intercepts and a common treatment effect can be specified as:

library(brms)

fit <- brm(
  y ~ treatment + (1 | center),
  data = dat,
  family = gaussian(),
  chains = 4,
  iter = 4000,
  seed = 1234
)

The expression:

(1 | center)

specifies a varying intercept by center.

A Random-Slope Model in R

If the treatment effect is expected to vary across centers, the model can include a varying treatment slope:

fit <- brm(
  y ~ treatment + (1 + treatment | center),
  data = dat,
  family = gaussian(),
  chains = 4,
  iter = 4000,
  seed = 1234
)

The model now estimates:

• A population-average intercept
• A population-average treatment effect
• Center-specific intercept deviations
• Center-specific treatment-effect deviations
• Between-center variability
• Correlation between center intercepts and treatment effects

Inspecting the Posterior

Posterior summaries can be obtained with:

summary(fit)

Posterior draws can be extracted with:

draws <- as_draws_df(fit)

A posterior distribution for the overall treatment effect can then be summarized using the model's posterior draws.

Posterior Probability of Clinical Benefit

Suppose the clinically important treatment effect is \(\delta\). The Bayesian question may be:

For example, if the treatment effect is stored as b_treatment:

mean(
  draws$b_treatment > 0.20
)

This calculates the proportion of posterior draws in which the treatment effect exceeds 0.20.

If the result is 0.94, the posterior probability is approximately 94%.

Extracting Center-Specific Effects

For a varying-slope model, posterior draws for center-specific effects can be examined using:

coef(fit)$center

This allows investigators to examine the posterior treatment effect for each center.

The key advantage is that these estimates are not independent estimates based only on each center's data.

They are informed by both center-specific observations and the population-level hierarchical distribution.

Posterior Predictive Checking in R

Posterior predictive checks can be performed using:

pp_check(fit)

Additional checks can examine specific aspects of the outcome distribution.

For example:

pp_check(
  fit,
  type = "dens_overlay"
)

The goal is not to make the replicated data look identical to the observed data.

The goal is to determine whether the model reproduces important features of the observed data.

Comparing Complete and Hierarchical Models

Suppose there are two competing models.

A complete-pooling model:

fit_pool <- brm(
  y ~ treatment,
  data = dat,
  family = gaussian()
)

and a hierarchical model:

fit_hier <- brm(
  y ~ treatment + (1 + treatment | center),
  data = dat,
  family = gaussian()
)

The models make different assumptions about center-level variation.

The hierarchical model should not automatically be preferred simply because it is more flexible.

Model comparison and predictive diagnostics should be considered alongside scientific plausibility.

Model Complexity

A hierarchical model can become very complex.

For example, a clinical trial might have:

• Patients nested within centers
• Centers nested within regions
• Repeated measurements within patients
• Multiple treatment arms
• Biomarker-defined subgroups

A model could theoretically represent all of these levels.

However, complexity should be justified by the scientific question and the available data.

Three-Level Clinical Data

Consider patients nested within centers, with centers nested within regions.

A three-level model could be:

where:

• \(i\) indexes patients
• \(j\) indexes centers
• \(k\) indexes regions
• \(u_k\) is the regional effect
• \(v_{jk}\) is the center-within-region effect

The variance structure can be:

This separates regional and center-level heterogeneity.

Crossed Hierarchical Structures

Not every clinical dataset has a simple nested structure.

For example, patients may be associated with:

• Clinical centers
• Investigators
• Manufacturing lots

where these classifications are not strictly nested.

Hierarchical models can represent multiple random effects when appropriate.

For example:

The scientific interpretation of each source of variation should be carefully considered.

Hierarchical Models for Multiple Treatment Arms

Suppose a trial evaluates several experimental doses. Let:

The dose-specific effects could be modeled jointly:

This can provide regularization across doses.

Alternatively, a dose-response model may be more scientifically appropriate. For example:

A hierarchical model can then allow parameters such as \(E_{\max}\) or \(ED_{50}\) to vary across patients or studies.

Hierarchical Models in PK/PD Research

Hierarchical Bayesian modeling is particularly natural in pharmacometric applications.

For example, individual pharmacokinetic parameters may be modeled as:

and:

Patient-level parameters vary around population-level parameters.

Covariates can then explain some of the between-patient variability.

This same hierarchical logic underlies population PK and many PK/PD models.

Borrowing Information Across Patients

Hierarchical models can also improve estimation of individual-level parameters.

Suppose patient \(i\) has limited observations. A hierarchical population model can combine:

• The patient's own observations
• The population distribution
• Covariate information

The resulting individual estimate can be more stable than an estimate based only on the patient's observations.

This is conceptually similar to shrinkage in population PK.

Hierarchical Models and Missing Data

Bayesian hierarchical models can naturally incorporate uncertainty about missing outcomes through the joint probability model.

For example:

However, simply using a Bayesian model does not make missing-data assumptions disappear.

The analyst still needs to consider whether the assumed missingness mechanism is scientifically plausible.

Sensitivity analyses may be required under alternative assumptions.

Missing at Random Is Still an Assumption

A hierarchical model may assume:

where \(R\) represents the missingness indicator.

This is related to a missing-at-random assumption.

If missingness depends on unobserved outcomes even after conditioning on observed variables, a sensitivity analysis may be necessary.

Hierarchical Models and Informative Cluster Size

An additional complication arises when cluster sizes differ substantially.

For example, some clinical centers may enroll hundreds of patients while others enroll only a handful.

The interpretation of the population-level effect can depend on whether the scientific target weights:

• Patients equally
• Centers equally
• Regions equally

The estimand should therefore be clearly defined before interpreting a hierarchical model.

Prior Sensitivity in Hierarchical Models

Prior sensitivity is especially important when there are few groups.

Suppose only five centers contribute data. The between-center standard deviation \(\tau\) may be poorly identified.

Different plausible priors for \(\tau\) can therefore produce materially different amounts of shrinkage.

A responsible analysis should evaluate multiple scientifically defensible choices.

An Example of Prior Sensitivity

Suppose the primary prior is:

A sensitivity analysis might consider:

and:

The exact choices should be based on the clinical scale of the treatment effect.

The purpose is not to select the prior that produces the preferred conclusion.

The purpose is to determine whether the scientific conclusion is robust to reasonable prior assumptions.

Robust Hierarchical Models

A normal hierarchical distribution may not always adequately represent between-group heterogeneity.

If occasional groups genuinely have extreme effects, a heavy-tailed distribution may be more appropriate.

For example:

A Student-\(t\) distribution can permit more extreme deviations than a normal distribution.

This can be useful when the assumption that all group effects are tightly concentrated around a normal distribution is questionable.

Mixture Models and Non-Exchangeability

Sometimes the population may contain genuinely different classes of groups.

For example, centers might belong to two qualitatively different clinical practice environments.

A single exchangeable distribution may then be inappropriate.

Possible alternatives include:

• Including relevant covariates
• Using separate hierarchical populations
• Mixture models
• Non-exchangeable historical-data models

The important point is that borrowing should reflect scientific similarity, not merely statistical convenience.

Bayesian Hierarchical Models and Trial Design

Hierarchical Bayesian models can be incorporated into clinical trial design before data are collected.

Simulation can evaluate:

• Bias
• Mean squared error
• Coverage
• Posterior probability behavior
• Borrowing behavior
• Type I error under relevant null scenarios
• Power or assurance
• Decision probabilities

This is particularly important when the trial will use adaptive decisions or external information.

Operating Characteristics

Suppose a Bayesian trial declares success when:

The operating characteristics of this rule should be evaluated under a range of true treatment effects.

For example:

Simulation can determine how often the trial declares success under each scenario.

Borrowing and Type I Error

When historical information is borrowed, investigators may need to evaluate frequentist operating characteristics even when the primary analysis is Bayesian.

For example, simulation can evaluate the probability of declaring treatment benefit when the current treatment is actually ineffective.

This can be particularly important for regulatory discussions.

A Bayesian analysis can therefore be evaluated using both Bayesian quantities and frequentist operating characteristics.

A Simple Simulation Framework

A conceptual simulation study can follow this structure:

Common Clinical Applications

Advantages of Bayesian Hierarchical Models

Bayesian hierarchical models offer several important advantages.

• Partial pooling across related groups
• Stable estimates for small groups
• Direct modeling of between-group heterogeneity
• Formal incorporation of external information
• Natural posterior prediction for new groups
• Flexible multilevel structures
• Direct probability statements about clinically meaningful quantities
• Natural integration with decision analysis
• Useful framework for complex clinical datasets

Limitations of Bayesian Hierarchical Models

The flexibility of hierarchical models also creates important challenges.

• Models can be computationally demanding
• Few groups can make variance components weakly identified
• Prior specification can materially affect results
• Exchangeability assumptions may be questionable
• Excessive borrowing can produce misleading estimates
• Complex models can be difficult to validate
• MCMC diagnostics may be required
• Interpretation of conditional and marginal effects can be complicated
• Regulatory expectations may require additional operating-characteristic analyses

Common Misconception: "Hierarchical Means All Groups Are the Same"

This is incorrect.

The model:

explicitly allows groups to differ.

The parameter \(\tau\) measures how much they differ.

If \(\tau\) is large, the model supports substantial heterogeneity.

Hierarchical modeling therefore does not eliminate heterogeneity. It estimates it.

Common Misconception: "Shrinkage Means the Model Is Biased"

Shrinkage can feel counterintuitive because an estimate is intentionally moved away from the raw group-specific estimate.

But the purpose is to reduce the influence of noisy estimates.

For a small group, an extreme observed estimate may have high variance. A hierarchical model recognizes that uncertainty.

The resulting posterior estimate may therefore have better predictive performance than the unpooled estimate.

Common Misconception: "More Borrowing Is Always Better"

Borrowing is beneficial only when the groups are sufficiently related.

If two studies are fundamentally different, forcing them to share information can introduce bias.

The key question is: Would it be scientifically reasonable to believe that these groups could have arisen from the same underlying population of effects?

Common Misconception: "A Large Heterogeneity Parameter Means the Model Failed"

Not necessarily.

A large posterior estimate of \(\tau\) may be an important scientific result.

It may indicate that:

• Treatment effects genuinely vary across populations
• Clinical practice modifies treatment response
• Patient characteristics are important effect modifiers
• The populations may not be exchangeable

The correct response may be to investigate the source of heterogeneity rather than simply force a common-effect model.

Common Misconception: "Hierarchical Models Solve Small Sample Problems"

Hierarchical modeling can stabilize estimates, but it cannot manufacture information that does not exist.

With very little information, posterior results can still be strongly influenced by the prior.

The model should therefore communicate uncertainty rather than create a false impression of precision.

Common Mistakes

1. Ignoring the hierarchy in clustered data. Patients within the same center or repeated observations within the same patient are generally not independent.
2. Using complete pooling when meaningful heterogeneity is expected. Forcing all groups to share one effect can conceal important variation.
3. Using no pooling for many small groups. Independent estimates can become extremely unstable.
4. Assuming hierarchical models automatically produce unbiased estimates. They remain dependent on the likelihood, prior, exchangeability, and model structure.
5. Ignoring the prior on the heterogeneity parameter. The prior for \(\tau\) can materially affect shrinkage when the number of groups is small.
6. Borrowing from historically incompatible populations. Differences in standard of care, eligibility, endpoints, and follow-up can invalidate exchangeability assumptions.
7. Interpreting shrinkage as proof that groups are identical. Shrinkage reflects uncertainty and information sharing, not equality.
8. Ignoring divergent MCMC transitions. Computational warnings can indicate that the posterior has not been sampled reliably.
9. Relying only on convergence diagnostics. A converged model can still be scientifically inappropriate.
10. Failing to perform prior sensitivity analyses. Small numbers of groups can make hierarchical variance parameters weakly identified.
11. Confusing conditional and marginal treatment effects. This is especially important for nonlinear models such as logistic regression.
12. Adding random effects without a scientific reason. Complexity should correspond to the data-generating structure and estimand.

A Practical Workflow for Bayesian Hierarchical Modeling

How to Interpret a Hierarchical Bayesian Result

Suppose a clinical trial estimates an overall treatment effect:

with a 95% credible interval:

The posterior probability of exceeding a clinically meaningful threshold of 0.10 is:

The hierarchical model also estimates:

This suggests that treatment effects vary across centers, but the estimated between-center variation is substantially smaller than the overall treatment effect.

The correct interpretation should distinguish the population-level effect from the center-specific effects.

What Should Be Reported?

A clinical Bayesian hierarchical analysis should generally report:

• Primary estimand
• Hierarchy and clustering structure
• Patient-level likelihood
• Group-level distributions
• Population-level priors
• Prior for heterogeneity parameters
• Rationale for exchangeability assumptions
• Posterior estimation method
• Convergence diagnostics
• Posterior predictive checks
• Population-level estimates
• Group-specific estimates when clinically relevant
• Credible intervals
• Clinically meaningful posterior probabilities
• Prior sensitivity analyses
• Model sensitivity analyses
• Borrowing sensitivity when external information is used

Reporting Historical Borrowing

If historical data are used, the analysis should additionally describe:

• Source of historical data
• Eligibility criteria
• Endpoint definitions
• Follow-up differences
• Standard-of-care differences
• Exchangeability assumptions
• Hierarchical borrowing model
• Prior for between-study heterogeneity
• Sensitivity to alternative borrowing assumptions
• Operating characteristics under incompatibility scenarios

Hierarchical Models vs. Simple Pooling

Hierarchical Models vs. No Pooling

The Most Important Mathematical Idea

The fundamental hierarchical model can be written:

combined with:

and:

The resulting posterior is:

This equation captures the entire hierarchy.

The data inform the group effects. The group effects inform the population distribution. The population distribution controls the degree of information sharing.

A Mental Model for Remembering Hierarchical Bayes

Think of three levels.

The key idea is that information flows in both directions.

Individual observations inform their group. Groups collectively inform the population distribution. The population distribution regularizes individual groups.

Bottom Line

Bayesian hierarchical models provide a powerful framework for clinical research data in which observations are naturally organized into related groups.

The central structure is:

This structure allows the analyst to estimate group-specific effects while borrowing information across groups.

The resulting partial pooling often produces more stable estimates than either complete pooling or no pooling, particularly when some groups have limited information.

In clinical research, hierarchical Bayesian models can be useful for:

• Multi-center clinical trials
• Subgroup analyses
• Historical-data borrowing
• Rare disease studies
• Meta-analysis
• Repeated-measures studies
• Population PK/PD modeling
• Multi-region development programs
• Multiple treatment arms
• Longitudinal clinical outcomes

But hierarchical modeling is not simply a statistical technique for making small datasets appear more precise.

The credibility of the analysis depends on whether the hierarchy is scientifically reasonable, whether the exchangeability assumptions are defensible, whether the prior distributions are appropriate, whether the amount of borrowing is reasonable, and whether the computational implementation has been adequately diagnosed.

References

Gelman, A. & Hill, J. (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. & Rubin, D.B. (2013). Bayesian Data Analysis, 3rd ed. CRC Press.
McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with Examples in R and Stan, 2nd ed. CRC Press.
Bürkner, P.-C. (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80(1), 1–28.
Spiegelhalter, D.J., Abrams, K.R. & Myles, J.P. (2004). Bayesian Approaches to Clinical Trials and Health-Care Evaluation. Wiley.
Berry, D.A. (2006). Bayesian clinical trials. Nature Reviews Drug Discovery, 5, 27–36.
Neuenschwander, B., Branson, M. & Gsponer, T. (2008). Critical aspects of the Bayesian approach to phase I cancer trials. Statistics in Medicine, 27, 2420–2439.
Ibrahim, J.G., Chen, M.-H. & Sinha, D. (2001). Bayesian Survival Analysis. Springer.
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515–534.
Piironen, J. & Vehtari, A. (2017). Comparison of Bayesian predictive methods for model selection. Statistics and Computing, 27, 711–735.