Markov Chain Monte Carlo (MCMC) for Clinical Statisticians

Introduction

Bayesian statistics provides a natural framework for clinical research because it allows investigators to combine prior information with observed clinical data and obtain a posterior distribution for parameters of interest.

The fundamental Bayesian calculation is:

where:

• \(\theta\) represents the unknown model parameters.
• \(y\) represents the observed data.
• \(p(\theta)\) is the prior distribution.
• \(p(y\mid\theta)\) is the likelihood.
• \(p(\theta\mid y)\) is the posterior distribution.

For simple Bayesian models, the posterior distribution may have a closed-form solution. For many realistic clinical models, however, the posterior distribution cannot be written down in a convenient analytical form.

Examples include:

• Bayesian logistic regression
• Hierarchical clinical trial models
• Longitudinal mixed-effects models
• Pharmacokinetic and pharmacodynamic models
• Survival models with complex random effects
• Bayesian dose-response models
• Meta-analytic models with multiple hierarchical levels
• Models with nonlinear parameters

This is where Markov Chain Monte Carlo (MCMC) becomes important.

What Problem Is MCMC Solving?

Suppose we have a parameter vector:

and observed clinical data \(y\). The posterior distribution is proportional to:

The proportionality statement is important. The right-hand side can often be evaluated numerically even when the normalized posterior distribution cannot be obtained analytically.

The normalizing constant is:

For a one-dimensional parameter this integral may be manageable. For a 10-dimensional parameter vector it becomes substantially more difficult. For a hierarchical model containing hundreds or thousands of latent parameters, direct numerical integration may be impractical.

Posterior Expectations Are Integrals

Suppose we want the posterior mean of some function \(g(\theta)\):

MCMC approximates this integral using an average of simulated draws. If the MCMC algorithm produces draws:

then:

This is the central computational idea behind MCMC.

Why Not Just Simulate Directly From the Posterior?

If we could independently generate:

then ordinary Monte Carlo simulation would be enough.

The difficulty is that the posterior distribution is usually exactly the distribution we do not know how to sample from directly.

MCMC solves this by constructing a stochastic process whose stationary distribution is the desired posterior.

What Is a Markov Chain?

A Markov chain is a sequence of random variables:

with the defining property that the distribution of the next state depends only on the current state, not on the entire history.

Formally:

The chain therefore moves through parameter space one step at a time.

Importantly, the draws generated by a Markov chain are generally not independent.

The Target Distribution

The distribution we want the chain to converge to is called the target distribution. In Bayesian analysis, the target is usually the posterior:

The MCMC algorithm is constructed so that, after sufficient iterations, the distribution of the chain is approximately \(\pi(\theta)\).

Stationarity

A distribution \(\pi\) is stationary for a Markov chain if applying the transition mechanism leaves that distribution unchanged.

If the chain has transition kernel \(K\), stationarity means:

In practical Bayesian computation, the goal is to construct a transition mechanism for which the posterior distribution is stationary.

Once the chain has converged sufficiently toward this stationary distribution, its draws can be used for posterior inference.

The MCMC Workflow

Burn-In and Warm-Up

Suppose the chain begins at:

There is no guarantee that this initial value resembles a draw from the posterior distribution.

Early iterations may therefore be strongly influenced by the starting point.

Historically, analysts often discarded an initial portion of the chain called the burn-in period.

Modern samplers such as Hamiltonian Monte Carlo typically use the term warm-up because these iterations are not merely discarded. They are also used to adapt tuning parameters such as step size and, depending on the algorithm, the geometry of the posterior.

Metropolis-Hastings

The Metropolis-Hastings algorithm is one of the foundational MCMC methods. Suppose the current state is:

We propose a new state:

The proposal is then accepted with probability:

If the proposed value is accepted:

Otherwise:

Why the Acceptance Ratio Works

Suppose the proposal distribution is symmetric:

For example, a normal random-walk proposal:

The proposal terms cancel. The acceptance probability becomes:

If the proposed location has higher posterior density, it is accepted with probability 1. If it has lower posterior density, it may still be accepted.

A Simple Metropolis-Hastings Example

Suppose the target distribution is a standard normal distribution:

Start at:

Suppose the proposal is:

If the proposed value is \(\theta^\star=0.5\), then:

Therefore, the proposal is accepted with probability approximately 88.3%.

A proposal to \(\theta^\star=2\) would have a lower acceptance probability:

Thus, it would be accepted approximately 13.5% of the time.

Random-Walk Behavior

Metropolis-Hastings can behave like a random walk through parameter space. If the proposal variance is too small, the chain moves only tiny distances. If it is too large, most proposals may be rejected.

This tuning problem is one reason more advanced samplers have become popular for complex Bayesian clinical models.

Gibbs Sampling

Gibbs sampling takes advantage of models in which the full conditional distributions are easy to sample from. Suppose:

The algorithm alternates between:

and:

The process continues for many iterations.

Why Gibbs Sampling Is Useful

Consider a Bayesian normal model with conjugate priors. Some conditional distributions may have recognizable forms such as:

• Normal
• Gamma
• Inverse-Gamma
• Beta

If each full conditional can be sampled directly, Gibbs sampling can be very efficient.

However, modern clinical models often contain nonlinear relationships, hierarchical structures, constrained parameters, and strong posterior correlations. In those settings, Gibbs sampling may become inefficient or unavailable.

Metropolis-Hastings vs. Gibbs Sampling

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) takes a fundamentally different approach. Instead of proposing a new parameter value using a local random walk, HMC introduces auxiliary momentum variables and uses gradients of the log posterior to move efficiently through parameter space.

Define a potential energy:

and introduce momentum \(r\) with kinetic energy:

The Hamiltonian is:

The algorithm uses the geometry of this Hamiltonian system to make long, directed moves through the posterior.

Hamiltonian Monte Carlo Intuition

Imagine trying to explore a mountainous landscape. A random-walk algorithm repeatedly takes small steps and asks whether the new location is acceptable. HMC instead uses information about the slope of the landscape to travel efficiently through broad regions of the posterior.

The "slope" is supplied by the gradient:

This is particularly powerful in high-dimensional continuous models.

The Leapfrog Integrator

HMC generally uses a numerical integration scheme called the leapfrog integrator to simulate Hamiltonian dynamics. For step size \(\epsilon\):

Several leapfrog steps are combined to produce a proposed trajectory. The proposal is then subject to an acceptance mechanism that corrects for numerical integration error.

What Is NUTS?

The No-U-Turn Sampler (NUTS) is an adaptive extension of HMC. One challenge with ordinary HMC is selecting the number of leapfrog steps. Too few steps and the chain does not move far enough. Too many steps and computation is wasted.

NUTS automatically determines when a simulated trajectory has traveled far enough by detecting when it begins to turn back toward itself.

This makes NUTS particularly attractive for general-purpose Bayesian inference. Modern Bayesian software such as Stan commonly uses NUTS as its default sampling algorithm for continuous parameter models.

Why HMC Is So Important for Clinical Bayesian Models

Many clinical Bayesian models contain correlated parameters. For example:

• Random intercepts and random slopes
• Between-study heterogeneity
• Treatment-by-covariate interactions
• PK/PD parameters
• Hierarchical treatment effects
• Patient-level latent effects
• Nonlinear dose-response parameters

The posterior may therefore resemble a narrow, highly correlated ridge rather than a spherical distribution.

A random-walk sampler can struggle to move efficiently along such a ridge. Gradient-based samplers are often much better suited to this geometry.

A Clinical Example: Bayesian Logistic Regression

Consider a Phase II clinical study comparing treatment response according to a binary treatment indicator. Let:

• \(Y_i=1\) if patient \(i\) responds.
• \(Y_i=0\) otherwise.
• \(T_i=1\) for experimental treatment.
• \(T_i=0\) for control.

We model:

with:

The treatment odds ratio is:

Suppose we place normal priors on the coefficients:

The posterior is:

There is generally no convenient closed-form expression for the normalized joint posterior. MCMC provides a way to obtain draws from it.

What the MCMC Output Represents

Suppose the sampler produces 4,000 post-warm-up draws of \(\beta\):

Then each draw represents a plausible value of the treatment log odds ratio under the posterior distribution.

We can transform every draw:

The resulting collection:

approximates the posterior distribution of the treatment odds ratio.

Posterior Mean

The posterior mean can be estimated as:

This is simply the average of the retained MCMC draws.

Posterior Median

The posterior median is the 50th percentile of the retained draws.

Unlike a frequentist estimator, the posterior median is directly interpretable as a summary of the posterior probability distribution.

Credible Intervals

A 95% equal-tail credible interval can be estimated using the 2.5th and 97.5th percentiles of the posterior draws:

For example:

The Bayesian interpretation is: conditional on the model, data, and prior, 95% of the posterior probability lies within the reported interval.

Posterior Probability of Clinical Benefit

MCMC also makes probability statements about clinically meaningful thresholds straightforward. Suppose the clinically important odds ratio is 1.5. We can estimate:

directly from the MCMC draws:

mean(OR_draws > 1.5)

For example, if 3,720 of 4,000 retained draws exceed 1.5:

The posterior probability of clinically meaningful benefit would therefore be approximately 93%.

Posterior Predictive Inference

MCMC can also generate predictions for future patients. Suppose:

For each posterior draw \(\theta^{(s)}\), generate a simulated future outcome:

The resulting collection of simulated outcomes approximates the posterior predictive distribution:

This is especially useful for:

• Bayesian predictive probability
• Trial monitoring
• Future response prediction
• Clinical decision analysis
• Bayesian sample size determination
• Posterior predictive checks

Why MCMC Draws Are Correlated

Suppose the chain is:

Adjacent values tend to be similar because the chain moves gradually through parameter space.

This means the draws contain less information than the same number of independent draws.

The degree of dependence is related to the autocorrelation function.

Autocorrelation

For lag \(k\), the autocorrelation can be expressed conceptually as:

High positive autocorrelation means that draws separated by \(k\) iterations remain similar.

Low autocorrelation generally indicates more efficient exploration.

Effective Sample Size

The effective sample size, or ESS, adjusts the nominal number of MCMC draws for autocorrelation. A common approximation is:

where \(S\) is the nominal number of retained draws.

If:

but strong autocorrelation reduces the effective sample size to:

then the simulation contains information roughly comparable to 800 independent draws for that parameter.

Monte Carlo Standard Error

The MCMC estimate of a posterior quantity is itself subject to simulation error. For a posterior mean:

a rough relationship is:

A large effective sample size reduces Monte Carlo error.

This leads to an important distinction:

Convergence

Before using MCMC draws for inference, we need evidence that the chains have adequately explored the target posterior distribution. This is referred to as convergence.

Convergence is not something that can be established by a single diagnostic. A good assessment combines multiple pieces of evidence.

Trace Plots

A trace plot shows the sampled value of a parameter against iteration. A well-behaved trace often resembles a relatively stable, rapidly mixing "hairy caterpillar."

For example:

plot(beta_draws,
     type = "l",
     xlab = "Iteration",
     ylab = "Treatment log odds ratio")

Warning signs include:

• Long trends over time
• Distinct regions separated by long periods
• Very high autocorrelation
• Chains that appear to occupy different regions
• Sudden unexplained changes in behavior

Multiple Chains

Modern Bayesian workflows commonly run several chains from dispersed starting values. Suppose we run four chains:

If all chains converge to the same posterior distribution, their behavior should become broadly similar after warm-up.

The Gelman-Rubin Diagnostic: \(\hat R\)

The potential scale reduction factor, commonly written \(\hat R\), compares within-chain and between-chain variation. Conceptually:

When chains have converged, between-chain and within-chain variation should be similar. Modern implementations generally use a stricter rank-based version of the diagnostic.

A value close to:

is desirable.

Values materially above 1 indicate that additional sampling or model investigation may be necessary.

Effective Sample Size Diagnostics

A parameter can have:

• Excellent convergence but poor effective sample size.
• Large nominal sample size but strong autocorrelation.
• A reasonable mean estimate but insufficient ESS in the tails.

Modern Bayesian software often reports separate bulk and tail effective sample sizes.

This is particularly useful when the scientific conclusion depends on extreme posterior quantiles.

Why Tail ESS Matters Clinically

Suppose the clinical decision depends on:

or a 97.5th percentile. The tails of the posterior may be estimated less precisely than the center. A large bulk ESS does not automatically guarantee accurate tail probabilities.

For decision-making based on posterior tail probabilities, tail ESS and MCSE should therefore be considered.

Divergent Transitions

Divergent transitions are particularly important in HMC-based Bayesian analysis. They indicate that the numerical Hamiltonian trajectory encountered difficult posterior geometry that the sampler could not accurately resolve.

Divergences should not simply be ignored because the estimated posterior mean "looks reasonable."

Common Causes of Difficult Posterior Geometry

Difficult MCMC behavior can arise from:

• Strongly correlated parameters
• Weakly identified parameters
• Hierarchical variance parameters near zero
• Badly scaled predictors
• Weakly informative data
• Separation in logistic regression
• Nonlinear parameter relationships
• Funnel-shaped hierarchical posteriors

Many of these problems are fundamentally statistical or modeling problems, not merely computational problems.

Parameterization Matters

One of the most powerful tools for improving MCMC performance is reparameterization. Consider a hierarchical model:

A centered parameterization represents the latent effects directly. A non-centered parameterization can instead use:

The non-centered representation can dramatically improve sampling when the hierarchical parameters are weakly identified.

Prior Predictive Checking

MCMC is not only used after observing data. Before fitting the model, investigators should consider whether the prior distributions imply plausible clinical outcomes.

For example, if:

in a logistic regression, then the implied odds ratios may span an enormous range.

The prior predictive distribution asks:

MCMC can be used to simulate from this distribution.

If the simulated clinical outcomes are implausible, the prior or model should be reconsidered before analyzing the actual data.

Posterior Predictive Checking

After fitting the model, we can simulate replicated datasets:

We then compare the replicated data with the observed data. For example, in a binary-response study we might compare:

• Overall response rates
• Treatment-group response rates
• Number of responders
• Covariate-specific response patterns

Poor posterior predictive performance may indicate model misspecification even when MCMC convergence diagnostics look excellent.

R Implementation with Stan

For modern clinical Bayesian analyses, a practical workflow is to specify the model in Stan and use an interface such as cmdstanr or rstan.

Consider a simple treatment-response model.

library(cmdstanr)
library(posterior)

stan_code <- '
data {
  int<lower=1> N;
  array[N] int<lower=0,upper=1> y;
  vector[N] treatment;
}

parameters {
  real alpha;
  real beta;
}

model {
  alpha ~ normal(0, 2.5);
  beta  ~ normal(0, 2.5);

  y ~ bernoulli_logit(alpha + beta * treatment);
}

generated quantities {
  real odds_ratio;
  odds_ratio = exp(beta);
}
'

The model contains:

• An intercept \(\alpha\)
• A treatment effect \(\beta\)
• Normal prior distributions
• A Bernoulli likelihood using the logit link
• A generated treatment odds ratio

Prepare the Clinical Data

dat <- list(
  N = nrow(df),
  y = df$response,
  treatment = df$treatment
)

Here:

• response is coded 0/1.
• treatment is coded 0/1.

Compile the Model

mod <- cmdstan_model(
  write_stan_file(stan_code)
)

The Stan program is compiled into a model object that can then be sampled.

Run the MCMC Sampler

fit <- mod$sample(
  data = dat,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 2000,
  seed = 2026
)

This requests:

• 4 independent chains
• 1,000 warm-up iterations per chain
• 2,000 post-warm-up iterations per chain

The nominal number of post-warm-up draws is therefore:

But the effective sample size will generally be smaller than 8,000 because the MCMC draws are correlated.

Inspect the Posterior Summary

fit$summary(
  variables = c("alpha", "beta", "odds_ratio")
)

A typical summary contains quantities such as:

Extract Posterior Draws

draws <- fit$draws()

beta_draws <- as_draws_vector(
  draws,
  variable = "beta"
)

OR_draws <- as_draws_vector(
  draws,
  variable = "odds_ratio"
)

Once the draws are available in R, posterior inference becomes a calculation on the simulated distribution.

Calculate the Posterior Probability of Benefit

Suppose the clinically meaningful treatment effect is:

The posterior probability of benefit can be calculated as:

mean(OR_draws > 1.5)

This is one of the simplest and most powerful ways to use posterior draws.

Calculate a 95% Credible Interval

quantile(
  OR_draws,
  probs = c(0.025, 0.975)
)

This returns the empirical 95% equal-tail credible interval from the MCMC sample.

Calculate the Posterior Probability of Harm

Suppose harm is defined as:

Then:

mean(OR_draws < 0.80)

provides the posterior probability that the odds ratio is below 0.80.

Calculate the Probability of a Clinically Meaningful Effect

A clinical decision may involve both benefit and harm thresholds. Suppose:

is considered clinically meaningful benefit. Then:

prob_benefit <-
  mean(OR_draws > 1.5)

prob_benefit

Similarly:

prob_harm <-
  mean(OR_draws < 0.80)

prob_harm

The posterior distribution can therefore be mapped directly onto a clinical decision framework.

Trace Plot in R

library(bayesplot)

mcmc_trace(
  as.array(fit),
  pars = c("alpha", "beta")
)

The chains should generally appear to overlap and mix without persistent trends or separated regions.

Autocorrelation Plot

mcmc_acf(
  as.array(fit),
  pars = "beta"
)

This shows how strongly draws remain correlated at different lags.

Rank Plot

Rank plots are another useful way to examine whether multiple chains are sampling from similar distributions.

mcmc_rank_hist(
  as.array(fit),
  pars = "beta"
)

Substantial differences between chains can indicate incomplete mixing.

Posterior Density Plot

mcmc_dens(
  as.array(fit),
  pars = "odds_ratio"
)

The resulting density provides a visual representation of the posterior distribution of the treatment effect.

What Does "Good MCMC" Look Like?

A well-behaved analysis typically has:

• Chains that mix well.
• Trace plots without persistent trends.
• \(\hat R\) very close to 1.
• Adequate bulk ESS.
• Adequate tail ESS.
• Small Monte Carlo standard errors relative to posterior uncertainty.
• No problematic divergent transitions.
• No problematic maximum treedepth saturation.
• No pathological numerical warnings.

Increasing the Number of Iterations

Suppose the posterior estimates are noisy. The first response might be to increase:

the number of post-warm-up draws. This can reduce Monte Carlo error.

However, increasing iterations does not necessarily fix:

• Non-identifiability
• Bad parameterization
• Severe posterior geometry
• Incorrect priors
• Model misspecification
• Persistent divergent transitions

Those problems require model or sampler investigation.

Thinning

Traditional MCMC workflows sometimes recommended retaining every \(k\)-th draw to reduce autocorrelation. For example:

thin = 10

Modern workflows generally do not recommend routine thinning. Why? Because if storage is not a problem, retaining correlated draws preserves information that can be used to estimate the effective sample size.

Warm-Up Is Not Part of the Posterior Sample

Suppose each chain uses:

warm-up iterations and:

sampling iterations. Only the post-warm-up draws are normally used for posterior inference.

Thus four chains produce:

retained draws, not 12,000.

What Happens if Chains Do Not Converge?

Suppose:

for a clinically important parameter. The appropriate response is not simply to report the posterior estimate and move on.

Investigate:

• Trace plots
• Chain-specific posterior distributions
• Initialization
• Posterior geometry
• Parameter scaling
• Prior informativeness
• Model identification
• Potential multimodality

Then rerun the analysis after addressing the underlying problem.

Multimodality

Some posterior distributions have multiple modes. For example:

If an MCMC chain becomes trapped in only one mode, the resulting posterior sample can be misleading.

Multiple chains are especially useful in this setting.

Identifiability and MCMC

Some computational problems are really information problems. Suppose two parameters always appear together as a product:

and the data contain little information that separates them. The posterior may contain a long ridge of nearly equivalent parameter combinations.

MCMC may then exhibit strong correlations and slow mixing.

The solution may require:

• More informative data
• More informative but scientifically defensible priors
• Reparameterization
• Removing redundant parameters
• Reconsidering the model structure

MCMC and Bayesian Hierarchical Models

Hierarchical models are particularly important in clinical research. Suppose treatment effects are estimated across \(J\) studies:

where:

• \(\theta_j\) is the treatment effect in study \(j\).
• \(\mu\) is the overall mean treatment effect.
• \(\tau\) is between-study heterogeneity.

A Bayesian hierarchical model estimates the entire posterior distribution:

This can be difficult to evaluate analytically. MCMC allows the joint posterior to be explored.

Partial Pooling

One important feature of hierarchical Bayesian models is partial pooling. The study-specific effects are informed by both:

• The individual study data
• The population-level distribution

MCMC provides draws for:

allowing posterior inference for every study as well as the overall effect.

MCMC in Pharmacokinetic and Pharmacodynamic Models

MCMC is particularly relevant to nonlinear PK/PD modeling. Suppose a model contains parameters:

• Clearance \(CL\)
• Volume \(V\)
• Absorption rate \(k_a\)
• Maximum effect \(E_{\max}\)
• EC50
• Inter-individual variability parameters

The posterior may be:

with substantial nonlinear dependence among parameters.

MCMC can generate joint posterior samples that preserve those dependencies.

Posterior Prediction from PK/PD Models

Suppose a future concentration is:

For every MCMC draw:

The resulting distribution:

automatically propagates parameter uncertainty into the prediction.

Why Plug-In Estimates Can Be Misleading

Suppose two parameters are strongly negatively correlated. If we independently use their posterior means:

we may obtain a different result from the posterior predictive mean:

because generally:

This is especially important for nonlinear PK/PD and exposure-response applications.

MCMC and Clinical Decision Thresholds

Suppose a Bayesian clinical trial has a decision rule:

MCMC makes this calculation straightforward. For each posterior draw, evaluate:

Then:

The Bayesian decision can therefore be implemented directly from the posterior sample.

MCMC and Bayesian Predictive Probability

Suppose a clinical trial has enrolled \(n\) patients and we want the probability that the final study will satisfy a success criterion. For every posterior draw:

1. Draw a future response probability from the posterior draw.
2. Simulate future patient outcomes.
3. Evaluate whether the final success criterion is met.

The proportion of simulated future trials that succeed estimates the Bayesian predictive probability.

This connects MCMC directly to Bayesian trial monitoring and predictive probability.

MCMC Does Not Mean "Bayesian Simulation"

It is useful to distinguish MCMC from ordinary simulation.

MCMC is therefore a particular class of Monte Carlo computation, not a synonym for all Bayesian simulation.

Common MCMC Misconceptions

1. "MCMC gives independent posterior samples." Usually false. The draws are generally autocorrelated.
2. "A large number of iterations guarantees convergence." False. A poorly parameterized model can remain problematic even with many iterations.
3. "R-hat exactly equal to 1 proves the model is correct." False. R-hat evaluates aspects of chain behavior, not scientific model validity.
4. "MCMC is only necessary for complicated Bayesian models." Closed-form Bayesian models do not require MCMC, but MCMC becomes particularly valuable when analytical calculations are unavailable.
5. "The posterior mean is always the best summary." Not necessarily. Skewed or multimodal posteriors may be better described using medians, quantiles, probabilities, and the full posterior distribution.
6. "A credible interval is the same as a confidence interval." They have different probability interpretations.
7. "Divergences can be ignored if the parameter estimates look plausible." They should be investigated because they can indicate problematic posterior geometry.
8. "Increasing the number of iterations fixes all MCMC problems." No. Some problems require reparameterization or model revision.
9. "Thinning is required whenever autocorrelation is present." No. Autocorrelation is expected in MCMC; thinning is generally not a first-line solution.
10. "MCMC validates the statistical model." No. MCMC validates the computational approximation to the specified model; it does not establish that the model is scientifically appropriate.

Metropolis-Hastings vs. HMC

Gibbs vs. HMC

Discrete Parameters and MCMC

Many clinical Bayesian models contain discrete quantities such as:

• Latent class membership
• Mixture component assignments
• Unknown model structures
• Discrete treatment states

Gradient-based samplers such as HMC require continuous differentiable parameters. Consequently, models involving discrete latent variables may require marginalization or alternative sampling strategies.

Posterior Correlation

Suppose two posterior parameters satisfy:

The posterior is highly elongated. A random-walk sampler may move inefficiently because a step in one coordinate often needs to be accompanied by a compensating step in the other.

This is one reason posterior correlation is useful to examine during model diagnostics.

Pair Plots

mcmc_pairs(
  as.array(fit),
  pars = c("alpha", "beta")
)

Pair plots can reveal:

• Strong posterior correlations
• Nonlinear relationships
• Multimodality
• Outlying regions
• Potential problematic geometry

Scaling Predictors

Suppose a continuous covariate is measured in units ranging from 0 to 100,000. The model:

may be computationally less convenient than using a standardized predictor:

Then:

Scaling can improve numerical stability and make priors easier to interpret.

Prior Choice and MCMC

Prior distributions influence both statistical inference and computational behavior. Suppose a hierarchical standard deviation has:

Therefore, prior specification should be considered part of the computational design of a Bayesian model.

How Many MCMC Draws Are Enough?

There is no universal number. The required number depends on:

• Posterior geometry
• Autocorrelation
• Effective sample size
• Precision required for the inferential target
• Tail probabilities of interest
• Computational complexity

Instead of asking: "Did I run 10,000 iterations?" ask: "Is the Monte Carlo error sufficiently small for the quantity I need to report?"

Monte Carlo Error vs. Clinical Uncertainty

Suppose the posterior SD of an odds ratio is 0.45 but the MCSE of its posterior mean is only 0.003. Then:

which indicates that simulation error is small relative to posterior uncertainty.

If instead:

then the MCMC approximation itself may be too imprecise for reliable reporting.

MCMC for a Bayesian Clinical Trial

Suppose a trial estimates a treatment effect \(\theta\) and has a decision threshold \(\delta\). The Bayesian analysis might specify:

After MCMC sampling, the probability is estimated by:

The analysis then becomes:

Why the Posterior Sample Is So Powerful

Once the MCMC draws exist, many quantities can be computed without deriving new analytical formulas. For example:

• Posterior means
• Posterior medians
• Credible intervals
• Odds ratios
• Risk differences
• Relative risks
• Probabilities of clinical benefit
• Probabilities of harm
• Predicted outcomes
• Future trial success probabilities
• Expected utility

Each can be calculated by transforming the posterior draws.

Derived Clinical Quantities

Suppose MCMC produces posterior draws for two treatment response probabilities:

The posterior risk difference draws are:

The relative risk draws are:

The odds ratio draws are:

The posterior probability of a clinically meaningful risk difference can then be calculated directly.

RD <- p_treatment - p_control

mean(RD > 0.10)

quantile(
  RD,
  c(0.025, 0.975)
)

MCMC and Missing Data

Bayesian models can incorporate missing-data mechanisms directly. For example, missing outcomes may be treated as latent quantities under an explicit missing-data model. MCMC can sample these latent values jointly with the parameters.

However, Bayesian computation does not make missing-data assumptions disappear. The validity of the resulting inference still depends on the assumed missing data mechanism.

MCMC and Model Comparison

Posterior simulation can also support model comparison through quantities such as:

• Posterior predictive performance
• Expected log predictive density
• Leave-one-out cross-validation
• Other predictive criteria

Model comparison should focus on the scientific purpose of the model rather than simply choosing the model with the most complicated structure.

Common Clinical Reporting Language

A Bayesian clinical report might state:

"The posterior median treatment effect was 0.18, with a 95% credible interval of 0.04 to 0.32. The posterior probability that the treatment effect exceeded the prespecified clinically meaningful threshold of 0.10 was 0.91."

The corresponding MCMC diagnostics should also be available to support the computational validity of the reported posterior quantities.

What Should Be Included in a Bayesian Statistical Analysis Plan?

For an MCMC-based clinical analysis, the statistical analysis plan should prespecify:

• The likelihood
• The model structure
• Prior distributions
• Clinical decision thresholds
• Posterior quantities of interest
• Decision criteria
• Number of chains
• Warm-up strategy
• Sampling strategy
• Convergence diagnostics
• Effective sample size criteria
• Handling of divergences and other sampler diagnostics
• Posterior predictive checks
• Missing-data assumptions
• Sensitivity analyses

Sensitivity to Priors

A Bayesian clinical analysis should generally examine whether important conclusions are sensitive to reasonable alternative priors. Suppose the primary analysis uses:

A sensitivity analysis might consider a different scientifically defensible prior.

The question is not whether every possible prior gives exactly the same answer. The question is whether the substantive clinical conclusion is robust to reasonable prior assumptions.

MCMC Does Not Eliminate Statistical Judgment

MCMC is a computational tool. It does not decide:

• What the clinically meaningful treatment effect is.
• Which prior is scientifically appropriate.
• Whether the endpoint is clinically relevant.
• Whether missing-data assumptions are plausible.
• Whether a hierarchical structure makes scientific sense.
• Whether the model captures important sources of variation.

The algorithm can produce extremely precise samples from a posterior distribution that is based on poor assumptions.

A Practical MCMC Workflow for Clinical Statisticians

A Compact Diagnostic Checklist

Worked Clinical Interpretation

Suppose a Bayesian clinical analysis produces:

The posterior median suggests a treatment odds ratio of approximately 1.72. The 95% credible interval suggests substantial posterior uncertainty but places most posterior mass above an odds ratio of 1.

However, if the prespecified clinically meaningful threshold is 1.5, the posterior probability of exceeding that threshold is only 74%.

This illustrates why reporting only whether a credible interval excludes 1 can be insufficient. The clinically relevant threshold may be different from the null value.

Clinical Thresholds vs. Statistical Null Values

A statistical null value might be:

But the clinically meaningful threshold could be:

These are different questions:

asks whether the treatment effect is positive.

Whereas:

asks whether the treatment effect is large enough to be clinically meaningful.

MCMC and Frequentist Simulation

MCMC should not be confused with frequentist simulation studies. A frequentist simulation study might generate many hypothetical datasets under a known parameter value and evaluate:

• Type I error
• Power
• Bias
• Coverage
• Sample size

MCMC instead generates samples from a posterior distribution conditional on the observed data and model.

Why MCMC Is a Core Skill for Modern Clinical Statisticians

Bayesian methods are increasingly used for complex clinical models because they can naturally represent:

• Prior information
• Hierarchical structure
• Partial pooling
• Parameter uncertainty
• Predictive uncertainty
• Clinical decision probabilities
• Complex nonlinear relationships

MCMC is the computational machinery that makes many of these models feasible.

A clinical statistician does not necessarily need to derive every MCMC algorithm from first principles. But they should understand what the algorithm is doing and, critically, whether it appears to be doing it correctly.

What a Clinical Statistician Should Be Able to Explain

After learning MCMC, you should be able to explain:

• Why the posterior cannot be evaluated analytically.
• What distribution the sampler is targeting.
• Why the MCMC draws are dependent.
• What warm-up accomplishes.
• How multiple chains are used.
• What \(\hat R\) measures.
• What effective sample size means.
• What Monte Carlo error means.
• Why divergent transitions matter.
• Why parameterization can affect computation.
• How posterior draws become clinical estimates.
• How posterior probabilities are calculated.
• How posterior predictive distributions are generated.
• Why convergence does not establish model validity.

Common MCMC Reporting Mistakes

1. Reporting the posterior without diagnostics. The reader should have evidence that the computational approximation was adequate.
2. Reporting only the number of iterations. Effective sample size and Monte Carlo error are more informative.
3. Ignoring divergences. Problematic HMC diagnostics require investigation.
4. Using one chain only. Multiple chains provide stronger evidence about convergence.
5. Calling posterior draws independent. MCMC samples are generally autocorrelated.
6. Confusing posterior uncertainty with Monte Carlo error. They represent different sources of uncertainty.
7. Failing to check prior sensitivity. Clinical conclusions can sometimes depend materially on prior assumptions.
8. Using posterior means in nonlinear models without preserving joint uncertainty. This can produce misleading predictions.
9. Assuming convergence means the model is correct. MCMC only addresses computational sampling from the specified model.

The Most Important Concept

The most important conceptual point is simple: MCMC converts an analytically difficult posterior distribution into a large simulated sample that can be analyzed empirically.

Instead of trying to calculate:

in closed form, we generate:

whose distribution approximates the posterior.

Then:

and:

The entire framework follows from this idea.

Final Summary

Markov Chain Monte Carlo is a computational framework for obtaining samples from probability distributions that are difficult to sample from directly. In Bayesian clinical statistics, the target distribution is usually the posterior distribution:

MCMC constructs a Markov chain whose stationary distribution is the posterior.

Important MCMC algorithms include:

• Metropolis-Hastings: proposes new values and accepts or rejects them according to an acceptance probability.
• Gibbs sampling: samples sequentially from full conditional distributions.
• Hamiltonian Monte Carlo: uses gradients and simulated Hamiltonian dynamics to explore continuous posterior distributions efficiently.
• NUTS: adaptively selects HMC trajectory lengths and is widely used in modern Bayesian computation.

Once posterior draws are available, they can be used to calculate:

• Posterior means and medians
• Credible intervals
• Posterior probabilities
• Clinical decision probabilities
• Risk differences
• Odds ratios
• Relative risks
• Posterior predictions
• Predictive probabilities
• Expected clinical outcomes

But generating draws is only part of the analysis. A valid MCMC workflow also requires evaluation of:

• Trace plots
• Multiple-chain behavior
• \(\hat R\)
• Bulk effective sample size
• Tail effective sample size
• Monte Carlo standard error
• Divergent transitions
• Posterior geometry
• Posterior predictive performance

References

Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. & Teller, E. (1953). Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21, 1087–1092.
Hastings, W.K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97–109.
Gelfand, A.E. & Smith, A.F.M. (1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association, 85, 398–409.
Gelman, A. & Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472.
Hoffman, M.D. & Gelman, A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, 1593–1623.
Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv:1701.02434.
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. & Rubin, D.B. (2013). Bayesian Data Analysis, 3rd ed. Chapman & Hall/CRC.
McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with Examples in R and Stan, 2nd ed. CRC Press.
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved \(\hat R\) for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.