Monte Carlo Checks for Stagewise P-Values

Why Check a Stagewise P-Value by Simulation?

A stagewise p-value is a null probability defined by the observed stopped outcome and a specified ordering of stopped outcomes. Numerical integration can evaluate that probability accurately, but its implementation may contain errors: a stopping condition can be reversed, a path can be assigned to the wrong look, or the ordering event can be coded incorrectly. Monte Carlo simulation provides an independent check by generating paths under the null, applying the stopping rule, and counting the simulated outcomes that are at least as extreme as the observed outcome under the same ordering.

The simulation is a validation device, not a new definition of the p-value. It must preserve the null model, look schedule, stopping rule, and ordering used in the calculation being checked. If any of these differ, the two results need not agree, even if both computations are internally correct.

Definition: Let \(E_{\mathrm{obs}}\) be the event that a null-generated stopped path is at least as extreme as the observed stopped outcome under the specified ordering. The stagewise p-value is \(p_{\mathrm{stage}}=P_{H_0}(E_{\mathrm{obs}})\). For \(M\) independently simulated null paths, let \(K\) be the number in \(E_{\mathrm{obs}}\). The Monte Carlo estimate is \(\hat p_{\mathrm{stage}}=K/M\).

Null Paths and the Stopping Rule

For the examples, suppose the null model is the canonical Gaussian group-sequential model with two analyses at information levels \(I_1<I_2\). Under \(H_0\), the standardized statistics \(Z_1\) and \(Z_2\) have mean zero, unit variances, and correlation \(\rho=\sqrt{I_1/I_2}\). Equivalently, they can be generated from independent increments of a standard Brownian motion \(W(t)\):

Z_j=\frac{W(I_j)}{\sqrt{I_j}},\qquad \operatorname{Corr}(Z_1,Z_2)=\sqrt{\frac{I_1}{I_2}}.

In a simulation, set \(I_0=0\), generate independent increments \(\Delta W_j\sim N(0,I_j-I_{j-1})\), form \(W(I_j)=\sum_{r=1}^{j}\Delta W_r\), and standardize at each look. This construction gives the required dependence; generating independent standard normal values for \(Z_1\) and \(Z_2\) would not.

Consider a rule that stops for efficacy if \(Z_1\geq 2\), and otherwise continues to the second look. The simulated outcome is therefore either an early-stopped outcome recorded at look 1 or a final outcome recorded at look 2. To check a stagewise p-value, the code must first apply this rule to every path and then evaluate the ordering event for that path. It must not treat every simulated \(Z_2\) as observed: paths stopped at look 1 have no second-look outcome in the stopped experiment.

Generate a null path.
Draw the independent Gaussian information increments and calculate the correlated \(Z\)-statistics.

Apply the stopping rule.
Record the first look at which the path stops, or the final look if it continues.

Apply the ordering.
Set an indicator to 1 if the stopped outcome is at least as extreme as the observed outcome, and 0 otherwise.

Average the indicators.
The proportion of indicators equal to 1 estimates the null probability defining the stagewise p-value.

Why the Simulation Proportion Estimates the Target

For simulated path \(m\), define \(Y_m=1\) if its stopped outcome belongs to \(E_{\mathrm{obs}}\), and \(Y_m=0\) otherwise. Under correct simulation, each \(Y_m\) is Bernoulli with success probability \(p_{\mathrm{stage}}\). Independent paths make the \(Y_m\)'s independent, so

E(\hat p_{\mathrm{stage}})= E\left(\frac{1}{M}\sum_{m=1}^{M}Y_m\right) =p_{\mathrm{stage}}, \qquad \operatorname{Var}(\hat p_{\mathrm{stage}}) =\frac{p_{\mathrm{stage}}(1-p_{\mathrm{stage}})}{M}.

Thus the simulation proportion is unbiased, and its Monte Carlo standard error is estimated by

SE_{\mathrm{MC}}= \sqrt{\frac{\hat p_{\mathrm{stage}}(1-\hat p_{\mathrm{stage}})}{M}}.

This standard error describes variation from using a finite number of simulated paths, conditional on the assumed null model and the event being evaluated. It does not measure uncertainty about whether the model, ordering, or stopping rule is appropriate. Nor is it the standard error of a treatment-effect estimate.

Worked Example: Checking an Early-Stopped Outcome's Ordering Event

Suppose \(I_1/I_2=0.5\), so \(\rho=\sqrt{0.5}\approx0.7071\), and the efficacy rule is to stop at look 1 when \(Z_1\geq2\). The observed outcome is at look 2, and the specified stagewise ordering says that outcomes at least as extreme are \(E_{\mathrm{obs}}=\{Z_1\geq2.1\}\cup\{Z_1<2,\ Z_2\geq2.1\}\). The first event counts only sufficiently large early stops; the second counts sufficiently large final statistics among paths that continued.

Under the stated Gaussian null model, \(Z_2\mid Z_1=x\sim N(\rho x,1-\rho^2)\). Therefore the event probability can be evaluated by one-dimensional integration:

p_{\mathrm{stage}} =P(Z_1\geq2.1)+ \int_{-\infty}^{2}\phi(x) \left[1-\Phi\left(\frac{2.1-\rho x}{\sqrt{1-\rho^2}}\right)\right]\,dx.

The first term is \(P(Z_1\geq2.1)\approx0.01786\). Numerical integration gives approximately \(0.01144\) for the second term, hence \(p_{\mathrm{stage}}\approx0.02930\).

Now simulate \(M=100{,}000\) null paths using the independent-increment construction, apply the look-1 stopping rule, and count paths in the displayed event. Suppose \(K=2{,}931\) paths meet it. Then

\hat p_{\mathrm{stage}}=\frac{2{,}931}{100{,}000}=0.02931,\qquad SE_{\mathrm{MC}}\approx \sqrt{\frac{0.02931(1-0.02931)}{100{,}000}} =0.000533.

The simulation estimate differs from the integrated value by \(0.00001\), far less than one Monte Carlo standard error. This is good evidence that the simulation and integration evaluate the same event under the stated model; it is not proof that either implementation is free of all errors.

A Second Check: A Different Stagewise Event

Changing the observed outcome or the ordering changes \(E_{\mathrm{obs}}\), even when the null model and stopping rule stay fixed. The simulation must use the event for the particular p-value under review, rather than reusing a tail event from another analysis.

Worked Example: Checking a Final-Look Tail Probability

Use the same two-look null model with \(\rho=\sqrt{0.5}\) and the rule \(Z_1\geq2\) to stop early. Suppose the observed final outcome and its ordering give \(E_{\mathrm{obs}}=\{Z_1\geq2\}\cup\{Z_1<2,\ Z_2\geq1.75\}\). The first term includes all early efficacy stops as at least as extreme under this specified ordering; the second is the qualifying final-look tail among paths that continued.

Conditioning on \(Z_1=x\) gives the corresponding integration formula:

p_{\mathrm{stage}} =P(Z_1\geq2)+ \int_{-\infty}^{2}\phi(x) \left[1-\Phi\left(\frac{1.75-\rho x}{\sqrt{1-\rho^2}}\right)\right]\,dx.

The first term is \(P(Z_1\geq2)\approx0.02275\), and numerical integration gives approximately \(0.02965\) for the second. Thus \(p_{\mathrm{stage}}\approx0.05240\). The equivalent decomposition \(P(Z_2\geq1.75)+P(Z_1\geq2,\ Z_2<1.75)\) is a useful independent check: its terms are approximately \(0.04006\) and \(0.0123\), respectively, giving the same total to the stated precision.

For \(M=100{,}000\) simulated paths, suppose \(K=5{,}238\) fall in the event. Then \(\hat p_{\mathrm{stage}}=5{,}238/100{,}000=0.05238\), and \(SE_{\mathrm{MC}}\approx\sqrt{0.05238(1-0.05238)/100{,}000}=0.000705\). The difference from \(0.05240\) is \(0.00002\), consistent with ordinary Monte Carlo variation.

Interpreting Agreement and Diagnosing Disagreement

An empirical proportion will almost never equal a numerical integration result exactly. The relevant question is whether the discrepancy is plausible given the simulation's Monte Carlo variability. As a rough diagnostic, a difference of a few estimated Monte Carlo standard errors may arise by chance; a persistent, much larger difference warrants investigation. This comparison is a diagnostic, not a universal pass-or-fail test. In particular, a very large \(M\) can make a small discrepancy statistically noticeable without making it practically consequential.

A useful check is to repeat the simulation with independent random-number streams or a larger \(M\). If the discrepancy shrinks in the expected manner as simulation uncertainty decreases, that supports the implementation. If it remains, inspect the event and path generation rather than simply increasing \(M\). Common causes include:

Generating independent \(Z\)-statistics instead of the required correlated statistics.
Using the wrong information ratio or reversing an inequality at a stopping boundary.
Counting outcomes from paths that should have stopped at an earlier look.
Applying a different ordering, or treating ties differently, from the numerical calculation.
Comparing an unconditional probability with a probability conditional on reaching a particular look.

With continuous Gaussian statistics, exact ties at a fixed boundary have probability zero. In a discrete-statistic setting, ties may have positive probability, so the tie convention must be applied consistently in both calculations. A Monte Carlo check cannot settle which ordering or tie convention is scientifically intended; it can only check the calculation for the convention supplied.

Worked Example: Choosing a Simulation Size for a Useful Check

Suppose a numerical calculation gives a stagewise p-value near \(0.0293\), and the goal is a Monte Carlo estimate with an approximate 95% margin of error no larger than \(0.0004\). Using \(p(1-p)\approx0.03(0.97)=0.0291\), the normal-approximation planning calculation is

M\approx \frac{p(1-p)(1.96)^2}{(0.0004)^2} = \frac{0.0291(1.96)^2}{0.00000016} \approx698{,}700.

Thus about \(699{,}000\) independent paths are a reasonable planning target for that approximate precision, under the assumed model. If a run instead uses \(M=1{,}000{,}000\) and records \(K=29{,}310\) qualifying paths, then \(\hat p_{\mathrm{stage}}=0.02931\) and \(SE_{\mathrm{MC}}\approx\sqrt{0.02931(0.97069)/1{,}000{,}000}=0.000169\). An approximate 95% Monte Carlo margin is \(1.96(0.000169)\approx0.000331\). The estimate is compatible with an integration value near \(0.02930\).

The calculation is approximate: it uses a normal margin and an assumed probability near the target. It helps choose a simulation budget, but it does not replace checking the observed hit count or validating the path-generation code.

What Monte Carlo Can and Cannot Establish

Numerical integration and Monte Carlo are complementary. Integration can be efficient and highly precise for a low-dimensional Gaussian probability, while simulation is flexible for complicated stopping rules and orderings. Because their numerical errors arise differently, agreement is reassuring. However, agreement does not establish that the null model is correct, that the chosen ordering is the desired one, or that the event corresponds to the intended stagewise p-value. Both methods can share the same mistaken event definition.

A simulation also becomes inefficient for very small tail probabilities: with limited \(M\), few qualifying paths may be observed, making the proportion noisy or even zero. A zero count means no event was seen in that finite run; it does not show that the true probability is zero. The appropriate interpretation depends on the simulation size and the uncertainty in the estimate. This check is most useful when enough paths are generated to make Monte Carlo variability small relative to the discrepancy one wants to detect.

Key takeaway: A valid Monte Carlo check simulates the null dependence, applies the actual stopping rule, and counts the event defined by the same stagewise ordering. Report \(K\), \(M\), \(\hat p_{\mathrm{stage}}\), and \(SE_{\mathrm{MC}}\), then judge agreement in light of that simulation uncertainty.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

Under the two-look canonical Gaussian null model, what correlation should be used if \(I_1/I_2=0.36\)? Why is independent generation of \(Z_1\) and \(Z_2\) incorrect?
A simulation has \(M=50{,}000\) paths, of which \(K=1{,}500\) meet the specified event. Calculate the estimated p-value and its Monte Carlo standard error.
Why must the stopping rule be applied before evaluating the ordering event on a simulated path?
A simulated estimate differs from the numerical value by four estimated Monte Carlo standard errors. Give two sensible checks before concluding that the integration is wrong.

Solutions: (1) The correlation is \(\sqrt{0.36}=0.6\). Independent generation would impose correlation zero and therefore simulate a different null joint distribution. (2) \(\hat p_{\mathrm{stage}}=1{,}500/50{,}000=0.0300\), and \(SE_{\mathrm{MC}}\approx\sqrt{0.03(0.97)/50{,}000}=0.000763\). (3) The ordering applies to the stopped outcome. A path that stops early has no later observed statistic, so treating it as though it continued can assign the wrong outcome and event membership. (4) Check that the simulated dependence and information ratio are correct, and verify the stopping inequalities and ordering event in both implementations. Also check whether the discrepancy persists across independent runs or a larger simulation; four standard errors merits investigation but does not by itself identify which calculation is wrong.