Exact Enumeration and Simulation for Predictive Probabilities

Why Compare Enumeration and Simulation?

An interim predictive probability is a probability about what may happen among the patients still to be observed, conditional on the responses available now. With a binary response and a fixed number of future patients, the possible future response counts form a finite set. This makes exact calculation possible: evaluate the predictive probability at each count and add the probabilities for counts that meet the terminal success rule.

A second approach is Monte Carlo simulation: repeatedly draw a possible future outcome under the posterior predictive distribution and record whether it meets the rule. Simulation is useful when the event is complicated or exact summation is cumbersome, but its estimate has random error. That error matters when an interim decision depends on whether the predictive probability lies above or below a boundary.

Exact Enumeration Under the Beta-Binomial Model

Suppose the response probability \(p\) has a \(\mathrm{Beta}(a,b)\) prior. At an interim analysis, \(x\) responses have been observed among \(n\) patients. Under the model, patient responses are conditionally independent Bernoulli observations given \(p\). The posterior distribution is \(\mathrm{Beta}(A,B)\), where \(A=a+x\) and \(B=b+n-x\). If \(m\) patients remain, let \(Y\) be their total number of responses. The posterior predictive probability of \(Y=y\) is beta-binomial:

Here \(\mathrm{B}(\cdot,\cdot)\) is the beta function. The distribution accounts for uncertainty in \(p\): it averages the future binomial probability over the posterior distribution of \(p\). It is not the same as plugging a single estimate of \(p\) into a binomial distribution.

In some designs, terminal success is defined through a posterior criterion rather than a response-count boundary. For a fixed final sample size and a specified terminal rule, one can first identify which possible final response counts satisfy that rule, then sum their predictive probabilities. The calculation below uses a response-count boundary so that the event is explicit.

Why the Predictive Probability Has This Form

Conditional on a particular value of \(p\), the future count has a binomial distribution. Integrating that conditional probability over the posterior gives the predictive probability. In particular, the beta density is proportional to \(p^{A-1}(1-p)^{B-1}\), so multiplying it by the binomial likelihood for \(Y=y\) gives powers \(p^{A+y-1}\) and \((1-p)^{B+m-y-1}\). Integrating these powers from zero to one yields the beta function in the numerator:

The enumeration is exact apart from numerical rounding in evaluating the beta functions and summing the terms. As a check, the probabilities for \(y=0,\ldots,m\) must be nonnegative and sum to one, up to numerical precision.

Monte Carlo Estimation and Its Error

To simulate from the same predictive distribution, repeat the following steps independently \(R\) times: draw \(p^{(j)}\) from the interim posterior \(\mathrm{Beta}(A,B)\), then draw \(Y^{(j)}\) from \(\mathrm{Binomial}(m,p^{(j)})\). Define \(I_j=1\) if the simulated outcome meets the terminal success rule and \(I_j=0\) otherwise. The Monte Carlo estimate is the proportion of simulated outcomes meeting the rule.

Each \(I_j\) is a Bernoulli random variable with success probability \(q\), provided each replicate is generated from the specified posterior predictive distribution and the replicates are independent. Thus \(E(\widehat q_{\mathrm{MC}})=q\), and

The standard error describes simulation noise in estimating \(q\); it is not the posterior standard deviation of the treatment response probability \(p\), nor a measure of clinical uncertainty in the trial outcome. A rough normal-approximation interval for the Monte Carlo estimate is \(\widehat q_{\mathrm{MC}}\pm 1.96\widehat{SE}_{\mathrm{MC}}\). This approximation can be unreliable when \(R\) is small or the estimated probability is near zero or one; in those settings, use a binomial interval designed for the simulation successes and failures.

Choosing the Number of Replicates

For a fixed number \(R\) of independent replicates, the variance \(q(1-q)\) is at most \(1/4\), since \(q(1-q)\leq1/4\) for \(0\leq q\leq1\). This gives a conservative standard-error bound of \(1/(2\sqrt{R})\). A simulation size should be chosen for a stated precision target, not merely because a round number of replicates is customary.

Diagnostics, Alternatives, and Limitations

Exact enumeration is generally preferable for a single binary endpoint when the number of future patients is modest: it has no Monte Carlo sampling error and is straightforward to reproduce. Monte Carlo becomes attractive when the future event involves many components, complex dependence, or calculations that make a direct sum unwieldy. It can also provide a practical check on an exact implementation. Neither method repairs a misspecified prior, likelihood, or terminal success rule; both calculate a probability under the model supplied.

Useful diagnostics for simulation include checking that the posterior draws and future outcomes follow the intended distributions, verifying that the simulated event indicator corresponds exactly to the terminal rule, and comparing estimates across independent batches. Large differences between batches suggest that \(R\) may be inadequate, although similar batch estimates do not prove adequate precision. For rare events, a run with zero simulated successes does not establish that \(q=0\); it may simply mean that too few replicates were generated to observe the event.

For enumeration, check that all qualifying values of \(y\) are included and no others are included, particularly when the boundary is an integer count. Confirm that the predictive probabilities sum to one. If \(m\) is large, direct calculations of beta functions or binomial coefficients can overflow or underflow in ordinary arithmetic; stable implementations use log-gamma or log-probability calculations. Such numerical precautions improve computation but do not change the mathematical event being evaluated.