Worked MCMC Analysis Of A Small Binary-Endpoint Trial

From a Small Trial to a Posterior Summary

A binary endpoint is often summarized by the number of participants who experience an event in each randomized group. With a small sample, however, the observed event proportions can be unstable, and a baseline risk factor may be important for interpreting the treatment contrast. Bayesian logistic regression provides a model for the event probability conditional on treatment and baseline information, while MCMC can approximate the resulting posterior distribution when it has no convenient closed form.

The aim here is to work through one hypothetical trial from data specification to an odds-ratio summary. The focus is the complete analysis chain: define the estimand and model, specify priors, form the posterior, inspect sampling diagnostics, and interpret the resulting uncertainty. The numerical MCMC summaries below are illustrative output for the stated hypothetical data, not results from a real study.

Trial Data and Model

Suppose 24 participants are randomized equally to a new treatment or control. The endpoint is an adverse event by a prespecified follow-up time. A binary baseline covariate identifies participants at higher baseline risk. It was measured before randomization and is coded $x_i=1$ for higher risk and $x_i=0$ otherwise. Treatment assignment is coded $z_i=1$ for treatment and $z_i=0$ for control. The table gives the number of participants and events in each treatment-by-covariate cell.

Group	Baseline risk $x$	Participants	Adverse events
Control	0, lower risk	6	1
Treatment	0, lower risk	6	1
Control	1, higher risk	6	4
Treatment	1, higher risk	6	3

Model: For participant $i$, let $Y_i=1$ indicate an adverse event and $Y_i=0$ otherwise. Assume $Y_i\mid z_i,x_i\sim\operatorname{Bernoulli}(p_i)$, with $$ \operatorname{logit}(p_i)=\log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_Tz_i+\beta_Xx_i. $$ The coefficient $\beta_T$ is the treatment-versus-control log odds ratio conditional on baseline risk. This model includes no treatment-by-risk interaction: it assumes the conditional treatment odds ratio is the same for both covariate groups.

The model relies on independent participant outcomes conditional on their modeled covariates, a correctly coded endpoint and treatment indicator, and an adequate logistic mean model. The common treatment effect across the two risk groups is a substantive assumption, not a consequence of randomization. Randomization supports a causal interpretation of treatment assignment under the trial's design and conduct, but does not by itself guarantee that this particular regression model fits well. With only 24 participants, the model is deliberately small; adding many covariates or interactions would be difficult to support.

Likelihood and Prior Specification

For a given coefficient vector, the probability of each participant's observed outcome follows from the Bernoulli model. Writing $\eta_i=\beta_0+\beta_Tz_i+\beta_Xx_i$ and $p_i=\operatorname{logit}^{-1}(\eta_i)$, the individual log-likelihood contribution can be written in either of two equivalent forms:

y_i\log(p_i)+(1-y_i)\log(1-p_i) = y_i\eta_i-\log\{1+\exp(\eta_i)\}.

The second form is convenient for computation. For the four cells in the table, the log likelihood, omitting binomial coefficients that do not depend on the coefficients, is

\ell(\beta_0,\beta_T,\beta_X) = \sum_{j=1}^{4} \left[y_j\eta_j-n_j\log\{1+\exp(\eta_j)\}\right],

where $n_j$ and $y_j$ are the participant and event counts for cell $j$, and $\eta_j$ is the linear predictor for that cell. Omitting the binomial coefficients does not change the posterior for the coefficients because they are constant with respect to those coefficients.

For this worked analysis, assign independent normal priors with mean zero and standard deviation 1.5 to all three coefficients. These priors are proper and regularize very large coefficient values, a useful feature when event counts are small. They are not “flat” priors: for example, the prior for the treatment odds ratio $\exp(\beta_T)$ places 95% of its probability approximately between $\exp(-2.94)=0.0528$ and $\exp(2.94)=18.9$. The scale should be considered in light of the coding and clinical context, and a real analysis should examine whether conclusions are sensitive to defensible alternative priors.

\beta_0,\beta_T,\beta_X\ \overset{\mathrm{ind}}{\sim}\ N(0,1.5^2), \qquad \pi(\beta_0,\beta_T,\beta_X\mid y) \propto \exp\{\ell(\beta_0,\beta_T,\beta_X)\} \prod_{k\in\{0,T,X\}}\pi(\beta_k).

This posterior expression follows directly from Bayes' rule: posterior density is proportional to likelihood times prior density. The likelihood is not conjugate to the normal coefficient priors, so the posterior is not available as a standard closed-form distribution. MCMC supplies dependent draws whose distribution, when the sampler works as intended, approximates this posterior.

Worked Example: Evaluating the Log Posterior

Evaluate the log posterior up to an additive constant at $(\beta_0,\beta_T,\beta_X)=(-1.10,-0.25,1.20)$. The cell predictors, in the order control/lower risk, treatment/lower risk, control/higher risk, treatment/higher risk, are $-1.10,-1.35,0.10,-0.15$.

The event counts are $1,1,4,3$, and each cell has $n_j=6$. Substituting into the grouped log likelihood gives

\begin{aligned} \ell &=1(-1.10)+1(-1.35)+4(0.10)+3(-0.15)\\ &\quad-6\left[\log(1+e^{-1.10})+\log(1+e^{-1.35}) +\log(1+e^{0.10})+\log(1+e^{-0.15})\right]\\ &\approx -2.50-6(0.2873+0.2305+0.7444+0.6210)\\ &\approx -13.80. \end{aligned}

Ignoring the normalizing constants of the normal priors, their log-density contribution is

-\frac{(-1.10)^2+(-0.25)^2+(1.20)^2}{2(1.5)^2} =-\frac{2.7125}{4.5} \approx -0.6028.

Thus the log posterior up to an additive constant is approximately $-13.80-0.6028=-14.40$. This calculation does not establish that this coefficient vector is likely or that it is a posterior estimate; it shows how the likelihood and prior combine at one point in parameter space.

Sampling and Reading the Output

An HMC implementation can sample the three coefficients jointly. The earlier material on Hamiltonian Monte Carlo for logistic regression and divergent transitions provides the sampler background; the analysis here emphasizes what to check for this particular trial. A practical run uses multiple dispersed initial values, warmup to adapt the sampler, and retained draws from each chain. For illustration, suppose four chains are run and 4,000 post-warmup draws per chain are retained, for 16,000 retained draws in total.

The first requirement is not a particular number of iterations but adequate exploration of the posterior. Inspect trace plots for sustained movement through the same region across chains, rather than chains stuck at distinct values or showing long trends. Rank-normalized $\hat{R}$ close to 1 and adequate bulk and tail effective sample sizes provide complementary evidence; none is a substitute for checking the model and sampler behavior. Divergent transitions in HMC require investigation rather than being dismissed because a summary statistic looks acceptable.

Worked Example: Interpreting Illustrative MCMC Diagnostics

Suppose the four-chain output has rank-normalized $\hat{R}=1.001$ for $\beta_T$, bulk effective sample size 6,200, tail effective sample size 5,400, and no divergent transitions. The trace plots overlap and show repeated movement across the sampled range.

These results are reassuring for the treatment coefficient: the chains appear to have mixed across a common distribution, and the effective sample sizes indicate that the 16,000 correlated draws contain substantially fewer than 16,000 independent-draw equivalents. They do not prove convergence, verify the logistic model, or establish that the prior is appropriate. The absence of divergences is also a sampler diagnostic, not a test of the clinical assumptions.

If the posterior standard deviation of $\beta_T$ is about 0.74, the approximate Monte Carlo standard error of its posterior mean is

SE_{\mathrm{MC}}(\bar{\beta}_T)\approx\frac{0.74}{\sqrt{6200}} \approx 0.0094.

This estimates simulation error in the reported posterior mean, not uncertainty about the treatment effect itself. The latter is represented by the posterior distribution, which is much wider in this example.

Odds-Ratio Summary

For a one-unit increase in $z_i$, holding $x_i$ fixed, the log odds change by $\beta_T$. Therefore the conditional treatment odds ratio is $\exp(\beta_T)$. Transform each retained draw of $\beta_T$ by exponentiation to obtain draws from the posterior distribution of that odds ratio. Because exponentiation is increasing, quantiles transform directly: exponentiating the 2.5th and 97.5th percentiles of $\beta_T$ gives the corresponding equal-tailed 95% credible interval for the odds ratio.

Worked Example: Transforming the Treatment Effect

Suppose the illustrative MCMC output has posterior median $\beta_T=-0.24$, with equal-tailed 95% credible interval $[-1.70,1.18]$. The posterior median odds ratio is

\exp(-0.24)=0.787\approx 0.79,

and the 95% credible interval is

[\exp(-1.70),\exp(1.18)] \approx [0.183,3.25].

In this hypothetical analysis, the posterior median conditional odds of an adverse event under treatment are about 0.79 times those under control at the same baseline-risk level. The 95% credible interval ranges from about 0.18 to 3.25 for that odds ratio, and includes 1. Under the specified model and priors, the interval expresses posterior uncertainty; it is not a claim that the parameter has a 95% chance of lying in this particular interval after the interval has been computed. The wide range shows that these small data do not estimate the treatment contrast precisely. An odds ratio is also not a risk ratio, especially when the endpoint is not rare.

Diagnostics, Alternatives, and Limitations

In addition to sampler diagnostics, check the data and model. Confirm event coding, denominators, treatment allocation, and covariate timing. Compare observed event patterns with model-implied probabilities by treatment and risk group, and investigate whether influential observations or separation-like patterns are driving the fit. A treatment-by-risk interaction may be clinically important, but four cells of six participants each provide little information for estimating it reliably. Any extension should be prespecified or justified and accompanied by suitably cautious uncertainty assessment.

The logistic model uses the baseline covariate to make an adjusted conditional comparison. A simpler beta-binomial analysis of event proportions by arm can be a useful alternative when the aim is an unadjusted group comparison and covariate adjustment is not needed. It does not answer the same conditional question. Frequentist logistic regression is another way to estimate an adjusted odds ratio, but its interval and inferential interpretation differ from the Bayesian posterior summary. Polya-Gamma sampling is an alternative MCMC strategy for Bayesian logistic regression; it changes the computation, not the model or the need to check assumptions and prior sensitivity.

Sparse data can yield broad posterior intervals even when chains mix well. A proper prior prevents an improper posterior in many problematic settings, but does not create information in the trial or guarantee that the treatment estimate is robust. Report the endpoint definition, covariate coding, model, prior scales, sampler settings, diagnostic results, and odds-ratio uncertainty. Consider sensitivity to plausible prior scales, particularly when the posterior is materially influenced by the prior.

Common pitfall: A learner may see $\exp(\beta_T)=0.79$ and conclude that treatment reduces the event risk by 21%. That is not what an odds ratio says: it is a 21% lower odds on the ratio scale, not necessarily a 21% lower probability or risk. The correction is to state the odds contrast and, if clinically useful, obtain model-based event probabilities at specified covariate values and describe their uncertainty separately.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For the model in this tutorial, what does $\exp(\beta_T)$ compare, and what assumption makes that odds ratio common across baseline-risk groups?
Why can the binomial coefficients be omitted from the log posterior calculation in the worked likelihood example?
If the posterior median of $\beta_T$ is $-0.40$, calculate the corresponding odds-ratio median to two significant figures.
What does an effective sample size of 6,200 mean when 16,000 draws were retained, and what does it not mean?
Why is a well-mixing chain insufficient evidence that the fitted analysis is suitable for a clinical conclusion?

Solutions. (1) It compares treatment with control in the odds of an adverse event, conditional on the baseline covariate. The model has no treatment-by-covariate interaction, so the treatment coefficient does not vary with baseline risk. (2) Those coefficients depend only on the observed counts, not on the regression coefficients; multiplying the likelihood by a parameter-independent constant leaves the posterior density after normalization unchanged. (3) $\exp(-0.40)=0.670$, or about 0.67; this is an odds ratio, not a 33% reduction in risk. (4) It is an estimate of the information in the correlated draws in independent-draw equivalents for the relevant posterior summary. It is not 6,200 new trial participants and does not measure clinical uncertainty. (5) Mixing addresses computation, not whether endpoint coding, model assumptions, covariate choices, or prior scales are defensible for the trial question.