Implementing a Reusable Boundary-Integration Routine

Why Make Boundary Calculation a Routine?

A group-sequential calculation often needs to be repeated with different boundary shapes, information schedules, or alternative effects. A reusable routine reduces the chance that each repetition silently changes the probability model or counts a crossing incorrectly. Its essential task is to carry forward only the probability mass that has not crossed a stopping boundary, then integrate that mass to obtain crossing probabilities at each look.

The earlier material on survivor-density updates, numerical error budgets, and calibration establishes the calculation’s ingredients. Here the emphasis is the interface and control flow that make those ingredients a reliable computational procedure. The routine described below calculates probabilities for deterministic boundaries under a canonical normal model; it does not assign p-values or order stopped outcomes.

Model and assumptions: There are \(K\) planned analyses at fixed, strictly increasing information times \(0<t_1<\cdots<t_K\). The standardized statistic at look \(i\) is \(Z_i\). Under standardized drift \(\theta\), the vector has means \(E(Z_i)=\theta\sqrt{t_i}\), variances 1, and correlations \(\operatorname{Corr}(Z_i,Z_j)=\sqrt{t_i/t_j}\) for \(i<j\). Equivalently, the standardized process has independent increments on the information scale. Boundaries are deterministic, and the stated normal model and information times are treated as correctly specified.

At look \(i\), let \(l_i\) and \(u_i\) be the lower and upper boundaries, with \(l_i<u_i\). A value below \(l_i\) or above \(u_i\) crosses a boundary; a value strictly between them survives to the next look. A boundary may represent efficacy, futility, or another prespecified stopping action. The routine must therefore receive an action label for each side at each look. It must not assume that every crossing is an efficacy rejection.

Routine Contract: Inputs and Outputs

A useful interface separates design inputs, the probability model, and numerical controls. It should reject invalid settings before integration starts. At minimum, the routine needs:

Design: the number of looks \(K\), information times \(t_1,\ldots,t_K\), and lower and upper boundaries at each look.
Actions: labels identifying which boundary crossings count as efficacy, futility, or other stopping events.
Probability settings: the drift \(\theta\) for the requested calculation; use \(\theta=0\) for a null calculation.
Calibration settings, when requested: a target cumulative null efficacy probability, a boundary family indexed by one scalar parameter, and an initial bracket for that parameter.
Numerical controls: integration tolerances, a maximum refinement or evaluation limit, and a rule for reporting failure to meet tolerance.

The output should include, for every look, lower- and upper-crossing probabilities, the mass surviving that look, and cumulative probabilities for the requested stopping actions. A power calculation returns the cumulative probability of crossing an efficacy boundary under the specified alternative drift. A calibration calculation returns the selected boundary parameter and the achieved null probability, not just the parameter itself. Useful diagnostics include total probability accounted for, the estimated integration error, the calibration residual, and a status indicating whether all tolerances were met.

These outputs make an important distinction visible. A stage-specific crossing probability is the probability of stopping at that particular look; a cumulative crossing probability through look \(i\) is the sum of the relevant stage-specific probabilities through that look. The latter is the quantity to compare with a cumulative error target.

The Recursive Calculation

Write \(\mu_i=\theta\sqrt{t_i}\). Under the canonical model, conditional on \(Z_{i-1}=x\), the next statistic is normal with conditional mean \(\mu_i+\rho_i(x-\mu_{i-1})\) and conditional variance \(1-\rho_i^2\), where \(\rho_i=\sqrt{t_{i-1}/t_i}\). Thus the transition density is

k_i(y\mid x)= \frac{1}{s_i}\phi\!\left( \frac{y-\mu_i-\rho_i(x-\mu_{i-1})}{s_i} \right), \qquad s_i=\sqrt{1-\rho_i^2}.

Here \(\phi\) is the standard normal density. Let \(f_i(y)\) denote the unconditional density of reaching look \(i\) with value \(y\) without having crossed any earlier boundary. At the first look, and then at each subsequent look, the survivor density is

f_1(y)=\phi(y-\mu_1)\,\mathbf{1}\{l_1<y<u_1\}, \qquad f_i(y)=\mathbf{1}\{l_i<y<u_i\} \int_{l_{i-1}}^{u_{i-1}} f_{i-1}(x)k_i(y\mid x)\,dx.

The indicator removes values that cross at the current look. Integrating before this removal gives the mass arriving at the look; integrating after it gives the surviving mass. In particular, the lower- and upper-crossing probabilities at look \(i\) are obtained by integrating the arriving density over \((-\infty,l_i]\) and \([u_i,\infty)\), respectively. At the first look, the arriving density is \(\phi(y-\mu_1)\); at later looks, it is the integral in the survivor update before multiplication by the indicator.

Justification: The independent-increment normal model implies the stated conditional normal distribution for \(Z_i\) given \(Z_{i-1}\). The density of surviving to look \(i-1\) is \(f_{i-1}\). Multiplying by the conditional transition density and integrating over all surviving \(x\) therefore gives the density arriving at look \(i\). Restricting \(y\) to the current continuation interval gives \(f_i\). Since crossing and survival partition the arriving probability mass, their probabilities plus the survivor mass must equal the mass arriving at that look, apart from numerical integration error.

A two-sided design can count either side as efficacy, or assign different actions to the two sides. The routine should compute the side-specific probabilities first and then sum only the crossings labelled as efficacy when it reports power. Under the null, the same calculation supplies the cumulative null efficacy probability used for type I error calibration.

Algorithm and Calibration Logic

Validate the request.
Check that information times are positive and strictly increasing, boundaries satisfy \(l_i<u_i\), all action labels are recognized, and numerical tolerances are positive. Check that the requested calibration target lies between zero and one.

Set the drift and transition parameters.
For a null error calculation set \(\theta=0\); for power, use the specified alternative drift. Compute \(\mu_i=\theta\sqrt{t_i}\), \(\rho_i\), and \(s_i\) for the required looks.

Integrate arriving mass and classify crossings.
At each look, integrate the arriving density below \(l_i\), inside \((l_i,u_i)\), and above \(u_i\). Record the two crossing probabilities and use only the interior density for the next update.

Aggregate the requested probability.
Sum the stage-specific probabilities whose action labels match the requested event. For power, sum efficacy crossings under the alternative drift; for null calibration, sum null efficacy crossings.

Calibrate if requested, then report diagnostics.
Repeat the probability calculation as the calibration parameter changes. Stop when the target residual and integration error meet their tolerances, or report failure. Return boundaries, achieved probabilities, residuals, and probability-mass checks.

Calibration requires care. Suppose a family of boundaries is indexed by a scalar \(c\), and the cumulative null efficacy probability \(A(c)\) is continuous and monotone over a bracket \([c_{\min},c_{\max}]\). A bracketing root solver can then seek \(A(c)-\alpha=0\). For a conventional upper boundary, increasing \(c\) reduces the probability of crossing, so \(A(c)\) is nonincreasing. The routine should verify that the target is bracketed; it should not silently return an endpoint when the requested target cannot be attained. If the boundary family is not monotone, this calibration method is not justified without additional analysis.

Worked Examples

Worked Example: A Two-Look Survivor Update

Consider a hypothetical null calculation with \(t_1=0.5\), \(t_2=1\), and upper boundaries \(u_1=u_2=0\). There is no lower stopping boundary, so take \(l_1=l_2=-\infty\). Both upper crossings are efficacy events. The model assumptions hold by construction for this illustration: the times increase, the statistics follow the canonical joint normal model, and the boundaries are fixed.

At the first look, \(Z_1\sim N(0,1)\), so the probability of an efficacy crossing is \(P(Z_1\geq 0)=0.5\), and the survivor mass is \(P(Z_1<0)=0.5\). The correlation between looks is \(\rho_2=\sqrt{0.5/1}=0.7071\), and the conditional standard deviation is \(s_2=\sqrt{1-0.5}=0.7071\).

The mass surviving both looks is \(P(Z_1<0,Z_2<0)\). For two standard normal variables with correlation \(\rho\), this probability is \(1/4+\arcsin(\rho)/(2\pi)\). Substituting \(\rho=0.7071\), for which \(\arcsin(\rho)=\pi/4\), gives \(1/4+1/8=0.375\). Consequently, the probability of crossing at look 2 among all hypothetical trials is \(0.5-0.375=0.125\). Cumulative null efficacy crossing is \(0.5+0.125=0.625\), and the remaining survivor mass is \(0.375\). The masses sum to \(1.000\), as required. This deliberately liberal design is a useful check of stage-specific and cumulative bookkeeping, not a recommended error-controlled design.

Worked Example: Calibrating a One-Look Boundary

Suppose a hypothetical one-look, one-sided test uses \(t_1=1\) and an upper efficacy boundary \(c\). The target cumulative null efficacy probability is \(\alpha=0.025\). There is one analysis, so this is a special case of the routine with no recursive update. Under the null, \(Z_1\sim N(0,1)\), and the calibration equation is \(P(Z_1\geq c)=0.025\).

Equivalently, \(\Phi(c)=0.975\), where \(\Phi\) is the standard normal cumulative distribution function. Thus \(c=\Phi^{-1}(0.975)=1.960\) to three significant figures. Substitution gives \(P(Z_1\geq1.960)\approx0.0250\), so the achieved null probability matches the target to the displayed precision. A reusable routine can obtain this value by root-finding on \(A(c)-0.025\), while returning the achieved probability and residual as well as \(c\).

Worked Example: Evaluating Power at the Calibrated Boundary

Use the one-look design from the preceding example, but suppose the hypothetical alternative drift is \(\theta=0.4\). The model specifies \(Z_1\sim N(\theta\sqrt{t_1},1)=N(0.4,1)\), with the same fixed boundary \(c=1.960\). Because the only stopping event is efficacy, power is the probability of crossing that boundary under this alternative.

Standardizing under the alternative gives \(P(Z_1\geq1.960)=1-\Phi(1.960-0.400)=1-\Phi(1.560)\). Using \(\Phi(1.560)\approx0.9406\), power is approximately \(0.0594\). The routine must use the alternative drift for this calculation; reusing the null drift would return \(0.0250\), the null crossing probability, rather than power. Both probabilities refer to proportions of hypothetical repeated trials, and neither is a probability that a hypothesis is true.

Diagnostics, Alternatives, and Limitations

At each look, the arriving mass should equal lower crossing plus upper crossing plus surviving mass, within numerical tolerance. Across the full design, all mutually exclusive stopping probabilities plus final survivor mass should sum to one. A deficit suggests probability-mass leakage; a total above one suggests duplicated mass or an integration error. The error budget should include integration at every look and, during calibration, the uncertainty in the final probability difference from the target.

Other warning signs include a calibration residual that stops improving, a target outside the probabilities at the calibration bracket endpoints, and results that change materially when integration tolerances are tightened. A successful solver status alone is not validation: repeat calculations at finer numerical settings and, where appropriate, compare with direct multivariate normal calculations or simulation. A simulation check has its own Monte Carlo uncertainty and should be treated as an independent diagnostic rather than an exact benchmark.

Recursive integration is especially useful when many alternative drifts or boundary candidates must be evaluated, because each calculation reuses the same stagewise structure. Direct multivariate normal probability calculations provide an important comparison for manageable numbers of looks, but can be less convenient as the number of looks grows. Simulation can handle flexible stopping rules, but introduces random error and may be inefficient for small tail probabilities. No method removes the need to check that its probability model matches the design.

The routine is only as reliable as its assumptions and inputs. Incorrect information times, a misspecified canonical normal model, or boundaries that change in response to unmodelled data can invalidate its results. Fine grids or tight tolerances can also make integration expensive. Finally, a calibration result applies only to the boundary family, null model, and stopping-action labels used in that calibration; it does not automatically establish control for a different design.

A Common Implementation Mistake

Common pitfall: Treating a stage-specific crossing probability as the cumulative error rate. For example, in the two-look illustration, the second-look crossing probability is \(0.125\), but the cumulative null efficacy probability is \(0.5+0.125=0.625\). Calibration to a cumulative target must use the sum of efficacy crossings across all looks. The correction is to store stage-specific probabilities separately, aggregate them according to their action labels, and check that crossing plus survivor mass is conserved at every update.

A related mistake is to integrate the next-stage transition from all previous values, including values that already crossed. The recursion integrates only over the previous continuation interval. Otherwise, stopped probability mass is allowed to re-enter the calculation, inflating later crossing probabilities.

Check Your Understanding

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

For \(t_1=0.25\) and \(t_2=1\), what are \(\rho_2\) and the conditional standard deviation \(s_2\)?
Why does the survivor-density update integrate only over \(l_{i-1}<x<u_{i-1}\)?
In the two-look illustration, what is the cumulative null efficacy probability after the second look?
For a one-look upper boundary \(c=1.645\) under a standard normal null, what is the approximate null crossing probability?
Name two diagnostics that should be reported or checked for a numerical calculation.

Solutions. (1) \(\rho_2=\sqrt{0.25/1}=0.5\), and \(s_2=\sqrt{1-0.5^2}=\sqrt{0.75}=0.8660\). (2) Values outside the continuation interval crossed a boundary at the previous look, so they have already stopped and cannot contribute to a later crossing probability. (3) It is \(0.5+0.125=0.625\), including efficacy crossings at both looks. (4) \(P(Z\geq1.645)=1-\Phi(1.645)\approx0.0500\). (5) Examples include the probability-mass conservation check, the integration-error estimate, the calibration residual, and stability under tighter numerical tolerances.