Validating Recursive Calculations by Simulation

Why Validate a Recursive Calculation by Simulation?

A recursive calculation propagates probability through successive analyses, retaining only paths that have not crossed a stopping boundary. It can compute stagewise rejection probabilities and power efficiently, but implementation details such as transition distributions, boundary conventions, and stopping rules can be mishandled. Simulation offers a separate check: generate complete trial paths from the assumed joint distribution, apply the same stopping rule, and compare the observed crossing frequencies with the recursive probabilities.

This is a validation check, not a replacement for the recursive calculation. Simulation has random error, while numerical integration has approximation error. A useful comparison therefore reports both, rather than treating a difference between two estimates as proof that either one is wrong.

Specify the Model and the Stopping Rule

For a standard canonical group-sequential model, let \(t_j\) be the information fraction at analysis \(j\), with \(0<t_1<\cdots<t_K=1\). Write \(Z_j\) for the standardized test statistic at analysis \(j\). Under the null hypothesis, the joint distribution is multivariate normal, with mean zero, unit marginal variances, and

Under a simple alternative with standardized drift \(\theta\), \(Z_j\) has mean \(\theta\sqrt{t_j}\) and variance 1, with the same correlation structure. This model presumes the canonical independent-increment information process and correctly specified information fractions. A simulation check only validates calculations relative to the model simulated; it cannot establish that the model describes a real trial.

For example, let \(a_j\) be an upper efficacy boundary and \(b_j\) a lower futility boundary, with \(b_j<a_j\). At each analysis, apply the rule in the specified order: if the trial has not already stopped, declare efficacy if \(Z_j\ge a_j\); otherwise stop for futility if \(Z_j\le b_j\), when futility stopping is part of the design. Paths that cross neither boundary continue. A lower futility boundary is nonbinding if crossing it does not formally require stopping; the simulation must reflect that distinction. For two-sided efficacy, rejection occurs on either designated efficacy crossing. The simulated rule must match the rule used in the recursive calculation exactly.

One convenient way to generate paths is to use independent standard normal variables \(\varepsilon_j\). Set \(Z_0=0\) and \(t_0=0\), then generate

The first term carries forward the previous statistic, the second adds the alternative drift, and the last represents the new independent information. Under \(H_0\), set \(\theta=0\). This recursion generates the specified normal means, variances, and correlations; it does not mean that the \(Z_j\) values at different analyses are independent.

Estimate Stagewise Probabilities and Their Uncertainty

For \(N\) independent simulated paths, let \(C_j\) be the number that first cross an efficacy boundary at analysis \(j\). The estimated stagewise probability is \(\hat q_j=C_j/N\). The total estimated rejection probability is \(\hat q=\sum_j C_j/N\) when the stagewise events are mutually exclusive, as they are under first-crossing classification. Each simulated path either rejects or does not, so the total crossing count is binomial with success probability \(q\), and its Monte Carlo standard error is

The same formula applies to a single stagewise event using \(\hat q_j\). For the full vector of stagewise counts, the categories are multinomial: their estimates are not independent. In particular, a path cannot first reject at two different stages. Comparing each stage separately is useful, but should not be mistaken for a set of independent checks.

A rough 95% Monte Carlo interval for a probability is \(\hat q\pm1.96SE_{\mathrm{MC}}(\hat q)\), when the count is large enough for a normal approximation to be reasonable. It describes simulation uncertainty under the assumed model, not uncertainty about an unknown treatment effect or the probability that the design is correct. For rare crossings or small counts, use a binomial interval with better small-sample behavior rather than relying on this approximation.

Check Power Under an Alternative

To validate power, repeat the simulation under the specified alternative drift and count paths that reject according to the design. A path that stops for futility without crossing an efficacy boundary is not a rejection. If the design has two-sided efficacy boundaries, either specified efficacy crossing contributes to power. Thus, power is not necessarily the probability of crossing an upper boundary alone.

Separate Monte Carlo Error From Integration Error

The recursive result may itself depend on grid spacing, interpolation, truncation, or quadrature tolerances. Those numerical choices create integration error. Monte Carlo uncertainty instead arises because a finite number of random paths estimates a probability. The two sources should be assessed separately:

Diagnostics, Alternatives, and Limitations

A useful simulation check examines more than a single total. Compare the first-crossing frequencies by analysis, confirm that the counts sum to the total number of simulated paths when nonrejection outcomes are also classified, and inspect whether the estimates stabilize as \(N\) increases. Under \(H_0\), also check that the simulated marginal means are near zero and marginal variances near one; under the alternative, compare the means with \(\theta\sqrt{t_j}\). These checks can expose incorrect drift, information fractions, or transition variances before interpreting crossing probabilities.

Large disagreement relative to Monte Carlo uncertainty warrants investigation, but random fluctuation remains possible. Check for a reversed boundary inequality, an incorrect update of the conditional mean, accidental independent draws of the \(Z_j\), or a mismatch in whether futility is binding. Repeated simulation batches can reveal instability; increasing \(N\) reduces Monte Carlo standard errors at rate \(1/\sqrt{N}\), so reducing them by half requires about four times as many paths.

Simulation complements, rather than supersedes, recursive integration or direct multivariate normal probability calculations. Recursive calculations can provide precise results efficiently but require numerical choices to be controlled. Direct multivariate normal calculations are another numerical benchmark, as covered in the earlier material on direct multivariate normal probabilities. Simulation is often easy to implement and can validate the stopping logic, but may be inefficient for very small probabilities and never removes model or coding errors automatically. All three approaches can share a mistaken boundary specification, so verify inputs and event definitions independently.