Recognizing When Trial Models Need MCMC

Why Posterior Computation Can Stop Being Algebraic

Some trial models have a special pairing of likelihood and prior: multiplying them gives a posterior in a familiar distributional family. For example, with a binomial response model and a beta prior on the response probability, the posterior update is available in closed form. In more complex models, treatment effects enter through a nonlinear link, or parameters are coupled in the likelihood. The resulting posterior may not belong to a standard distributional family, even when the prior itself is simple.

The practical consequence is not that the posterior cannot be defined. Its density is still proportional to the likelihood times the prior. Rather, its normalizing constant and summaries may lack a convenient analytic expression. Numerical integration or numerical sampling can then be considered. MCMC is one broad sampling approach for such targets, but it is not automatically necessary: low-dimensional numerical integration or a suitable approximation may be feasible.

Definition: For parameter vector $\theta$, data $y$, likelihood $L(\theta\mid y)$, and prior density $\pi(\theta)$, the posterior kernel is $$ \pi(\theta\mid y)\ \propto\ L(\theta\mid y)\pi(\theta). $$ It is the posterior density up to a multiplicative constant that does not depend on $\theta$. A model is called nonconjugate for a specified prior when this product does not have a posterior distribution in the same standard family as that prior.

Nonconjugacy is a statement about a likelihood-prior pairing, not a permanent label attached to a likelihood alone. A binomial likelihood is conjugate to a beta prior on its probability parameter. But if probabilities are linked to regression coefficients, a prior specified on those coefficients need not produce a standard conjugate posterior. The parameterization and prior matter.

Logistic Treatment Effects and Their Posterior Kernels

Consider a parallel-group trial with binary outcomes. Let $Y_i=1$ indicate response for participant $i$, and let $x_i=0$ for control and $x_i=1$ for treatment. Assume the outcomes are conditionally independent given the model parameters and $$ Y_i\mid \alpha,\beta \sim \operatorname{Bernoulli}(p_i),\qquad \operatorname{logit}(p_i)=\alpha+\beta x_i. $$ Here $\alpha$ is the control-group log odds of response and $\beta$ is the treatment log odds ratio. The inverse-logit relationship is $p_i=\exp(\alpha+\beta x_i)/[1+\exp(\alpha+\beta x_i)]$.

For independent normal priors $\alpha\sim N(m_\alpha,s_\alpha^2)$ and $\beta\sim N(m_\beta,s_\beta^2)$, the posterior kernel is the product of the Bernoulli likelihood contributions and the two prior kernels. For observed binary outcomes, it can be written as follows.

\pi(\alpha,\beta\mid y)\ \propto\ \prod_{i=1}^{n}p_i^{y_i}(1-p_i)^{1-y_i} \exp\!\left[-\frac{(\alpha-m_\alpha)^2}{2s_\alpha^2} -\frac{(\beta-m_\beta)^2}{2s_\beta^2}\right].

This is generally not a bivariate normal density. In a normal model with a normal prior, the likelihood and prior kernels combine into a quadratic expression in the parameters. Here, the likelihood includes denominators $1+\exp(\alpha+\beta x_i)$. Those terms prevent the product from reducing to a normal kernel. A useful log-kernel, which avoids multiplying many small likelihood values, is

\log \pi(\alpha,\beta\mid y) = \sum_{i=1}^{n}\left[y_i(\alpha+\beta x_i) -\log\{1+\exp(\alpha+\beta x_i)\}\right] -\frac{(\alpha-m_\alpha)^2}{2s_\alpha^2} -\frac{(\beta-m_\beta)^2}{2s_\beta^2} +C,

where $C$ does not depend on $\alpha$ or $\beta$. The transformation from coefficients to response probabilities also couples interpretation: $p_0=\operatorname{logit}^{-1}(\alpha)$, while $p_1=\operatorname{logit}^{-1}(\alpha+\beta)$. The treatment risk difference is $p_1-p_0$, which depends on both coefficients, not on $\beta$ alone.

Worked Example: A Two-Arm Logistic Trial

Suppose, hypothetically, that 10 participants per arm have binary outcomes: 4 of 10 control participants respond and 7 of 10 treatment participants respond. Assume independent outcomes within arms conditional on their response probabilities, with $\operatorname{logit}(p_0)=\alpha$ and $\operatorname{logit}(p_1)=\alpha+\beta$. Use independent $N(0,2.5^2)$ priors for $\alpha$ and $\beta$. The data are compatible with this model; the calculation below evaluates the kernel at one parameter point, rather than estimating the posterior summary.

Ignoring binomial coefficients that do not depend on the parameters, the likelihood is

L(\alpha,\beta)\ \propto\ p_0^4(1-p_0)^6p_1^7(1-p_1)^3.

Evaluate at $\alpha=0,\ \beta=\log 2$. Then $p_0=1/(1+1)=1/2$, and $p_1=2/(1+2)=2/3$. The likelihood kernel is

(1/2)^{10}(2/3)^7(1/3)^3 =\frac{1}{1024}\frac{128}{59049} =\frac{1}{472392} \approx 0.000002117.

The product of the two normal prior kernels, with constants omitted, is

\exp\!\left[-\frac{0^2+(\log 2)^2}{2(2.5)^2}\right] =\exp[-(\log 2)^2/12.5] \approx 0.9623.

Thus the posterior kernel at this point is approximately $0.000002117(0.9623)=0.000002037$, using the unrounded factors for the multiplication. This value is not a posterior probability; a kernel value has no probability interpretation on its own. To obtain a posterior density, the kernel must be normalized over the parameter space.

A Short Justification for the Nonconjugacy

The logistic log-likelihood has a form that makes the obstruction visible. Write $\eta_i=\alpha+\beta x_i$. Since $p_i=\exp(\eta_i)/[1+\exp(\eta_i)]$, each Bernoulli contribution satisfies

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

Summing gives the log-likelihood. A normal prior adds a quadratic function of $\alpha$ and $\beta$. For a normal posterior kernel, the full log-kernel would need to be quadratic in those parameters. The term $\log(1+\exp(\eta_i))$ is not quadratic in $\eta_i$, so this sum is generally not quadratic either. Therefore the product is not, in general, a standard normal posterior kernel. This establishes why the usual normal-prior/normal-posterior update does not apply; it does not rule out numerical evaluation.

The same reasoning applies to more elaborate logistic regressions. Adding covariates changes each $\eta_i$, but the log-sum-exponential terms remain. A normal prior on a coefficient vector is still a reasonable prior specification in many settings, but it does not make the posterior normal.

Worked Example: Individual Covariates in a Logistic Model

Consider four hypothetical participants with covariate values $x=(0,0,1,2)$ and outcomes $y=(1,0,1,0)$. Assume conditional independence and the logistic model $\operatorname{logit}(p_i)=\alpha+\beta x_i$, with independent $N(0,2.5^2)$ priors. At $\alpha=0,\ \beta=\log 2$, the response probabilities for $x=0,1,2$ are $1/2,\ 2/3,\ 4/5$, respectively.

The likelihood kernel for the four records is

(1/2)(1-1/2)(2/3)(1-4/5) =\frac{1}{2}\frac{1}{2}\frac{2}{3}\frac{1}{5} =\frac{1}{30} \approx 0.03333.

The prior kernel at this point is $\exp[-(\log 2)^2/12.5]\approx0.9623$. The posterior kernel is therefore approximately $0.03333(0.9623)=0.03208$, using unrounded values. Each participant has a different probability determined by their covariate value; the likelihood is not one binomial count with one shared response probability. The product still defines a valid posterior kernel, but it does not yield a standard conjugate update for the regression coefficients.

Nonlinearity Beyond Logistic Regression

Unknown nonlinear parameters can also disrupt a conjugate update in trial models with time-to-event outcomes. For example, suppose event times follow a Weibull model with shape $k>0$ and rate-like scale parameter $\lambda>0$, using survival function $S(t)=\exp(-\lambda t^k)$. For an observed time $t_i$ and event indicator $\delta_i$, where $\delta_i=1$ denotes an observed event and $\delta_i=0$ right censoring, the likelihood contribution is the event density when $\delta_i=1$ and the survival probability when $\delta_i=0$. Assume independent observations and noninformative censoring conditional on the model parameters.

Let $D=\sum_i\delta_i$. The likelihood kernel is

L(\lambda,k)\ \propto\ \lambda^D k^D \left(\prod_{i:\delta_i=1}t_i^{k-1}\right) \exp\!\left(-\lambda\sum_i t_i^k\right).

If $k$ were fixed and $\lambda$ had a gamma prior, the factors involving $\lambda$ have a gamma form. When $k$ is unknown, however, both $\sum_i t_i^k$ and the event-time powers depend nonlinearly on $k$. A gamma prior for $\lambda$ alone does not provide a standard joint conjugate posterior for $(\lambda,k)$.

Worked Example: An Unknown Weibull Shape

Suppose, hypothetically, that three participants have observed times $1,2,3$ years, with events at years 1 and 3 and right censoring at year 2. Thus $D=2$. Use independent gamma priors with shape 2 and rate 1 for $\lambda$ and $k$, and evaluate the posterior kernel at $\lambda=0.2$ per year and $k=1$. The gamma priors are proper; the censoring assumption is that censoring is noninformative conditional on the parameters.

At $k=1$, $\sum_i t_i^k=1+2+3=6$, and the event-time powers are $1^{0}3^{0}=1$. Hence the likelihood kernel is

(0.2)^2(1)^2(1)\exp[-0.2(6)] =0.04\exp(-1.2) \approx0.01205.

A gamma density with shape 2 and rate 1 has kernel $z^{2-1}\exp(-z)$. The joint prior kernel at this point is $(0.2)e^{-0.2}(1)e^{-1}=0.2e^{-1.2}\approx0.06024$. Therefore the joint posterior kernel is $0.04e^{-1.2}\times0.2e^{-1.2}=0.008e^{-2.4}\approx0.0007256$. The pointwise calculation is straightforward, but it does not provide a closed-form joint posterior distribution: away from $k=1$, the powers $t_i^k$ and the event-time factors vary nonlinearly with $k$.

From a Kernel to a Computational Choice

Once a posterior kernel is available, the next question is whether it can be normalized or summarized adequately by analytic work, deterministic numerical integration, or numerical sampling. MCMC is useful when the posterior is too complicated for a direct update and numerical integration over all parameters is impractical. The kernel, rather than the unknown normalizing constant, is the essential target for many sampling algorithms. A stable implementation usually evaluates the log-kernel, because products of many small probabilities can underflow numerically.

A useful log-kernel for independent observations and a normal coefficient prior can be calculated by summing the individual log-likelihood contributions and adding the prior log-density, omitting constants independent of the coefficients. Before choosing a computational method, check the following:

Model and data alignment: verify outcome coding, treatment indicators, covariate scaling, and any offsets. A coding error changes the target posterior, not merely the numerical method.
Parameter support: respect restrictions such as $k>0$ and $\lambda>0$. A kernel evaluated outside its parameter space is not a valid posterior target.
Identifiability: determine whether the data can distinguish the coefficients. Sparse outcomes, strongly correlated covariates, or complete separation can produce broad, highly correlated, or extreme posterior regions.
Posterior propriety: check that the posterior is a proper probability distribution. An improper posterior cannot be repaired by running a sampling algorithm longer.
Numerical stability: calculate sums of log contributions and use numerically stable forms for the logistic log-likelihood when predictors are very large in magnitude.

These checks are not optional details. If a chain later appears not to explore the target well, a coding error, weak identification, or an inappropriate parameterization may be involved. Sampling diagnostics can reveal computational trouble, but they do not establish that the model is clinically appropriate or that the posterior is proper.

How the Main Alternatives Differ

Approach	When it is attractive	Main qualification
Conjugate analytic update	The likelihood-prior product belongs to a known family, allowing direct posterior summaries.	Not generally available for logistic coefficients with normal priors or for joint nonlinear survival parameters.
Deterministic numerical integration	The parameter dimension is small and accurate integration is feasible.	Computational cost can rise rapidly with dimension and complex posterior geometry.
Analytic approximation, such as a Laplace approximation	A compact approximation near a well-behaved posterior mode is sufficient for the purpose.	May be inaccurate for skewed, heavy-tailed, multimodal, or boundary-concentrated posteriors.
MCMC	The posterior kernel is available but direct computation is difficult, especially with several parameters.	Requires careful implementation and assessment of whether the generated draws adequately represent the target posterior.

These are alternatives, not a hierarchy in which one method is always best. The model dimension, posterior shape, desired accuracy, and inferential target all matter. A simple conjugate model may also be embedded in a larger nonconjugate model; then some parameters may have familiar conditional distributions even though the joint posterior is not conjugate.

Common Mistakes and Practical Limits

Wrong reasoning: “The prior is normal, so the posterior coefficients are normal.” Correction: a normal prior yields a normal posterior only under likelihood and parameter conditions that produce a quadratic log-kernel. Logistic likelihood terms generally do not.
Wrong reasoning: “Because the posterior cannot be written in closed form, it is undefined.” Correction: the likelihood-times-prior kernel defines the posterior whenever it can be normalized; the difficulty is computation, not definition.
Wrong reasoning: “Any model without a conjugate update requires MCMC.” Correction: quadrature, approximation, or other numerical methods may be suitable, especially in low dimensions.
Wrong reasoning: “A long MCMC run will fix separation or a misspecified likelihood.” Correction: computation cannot replace a defensible model, suitable priors, or an assessment of identifiability and posterior propriety.

Finally, a posterior kernel evaluated at a point is not itself a posterior probability, and its numerical magnitude depends on parameterization and omitted constants. Use the kernel to define and compute the target; interpret posterior summaries only after the distribution has been properly represented.

Check Your Understanding

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

In the two-arm logistic model, what are the control response probability and treatment response probability when $\alpha=0$ and $\beta=\log 2$?
Why does a normal prior on logistic regression coefficients generally not give a normal posterior?
In the Weibull example, which parameter has a gamma-shaped conditional likelihood-prior pairing when the shape $k$ is fixed?
Name one reason MCMC might be considered and one problem MCMC cannot fix.

Solutions. (1) $p_0=\operatorname{logit}^{-1}(0)=1/2$, while $p_1=\operatorname{logit}^{-1}(\log 2)=2/3$. (2) The logistic log-likelihood contains terms $\log(1+\exp(\eta_i))$, which are not quadratic in the coefficients; adding a quadratic normal-prior log-kernel therefore does not generally produce a normal kernel. (3) With $k$ fixed, the likelihood terms in $\lambda$ are proportional to $\lambda^D\exp[-\lambda\sum_i t_i^k]$; multiplying by a gamma prior gives a gamma-shaped kernel for $\lambda$. (4) MCMC may be considered when the posterior kernel is available but direct integration is impractical. It cannot fix a misspecified model, an improper posterior, or a coding error.