Simulation Diagnostics for Bias and Coverage

Why Simulate the Stopped Procedure?

A group sequential estimator is not just a formula applied to a fixed sample size. The observed sample size depends on interim results, so the estimate and interval are evaluated in a distribution shaped by the stopping rule. A simulation diagnostic makes that full procedure explicit: generate trial data, apply the rule, calculate the estimate and interval at the stage actually reached, and then summarize performance over many trials.

The quantities of interest should be specified before simulation. Stage-specific bias describes the average estimation error among trials that stop at a particular stage. Overall bias averages estimation error over all trials under the design. Coverage and interval length can likewise be summarized conditional on stopping stage or over the entire design. These summaries answer different questions; one should not substitute for another.

Model, Stopping Rule, and Performance Measures

Consider a hypothetical two-look trial. Individual outcomes \(X_i\) are independent and identically distributed as \(N(\theta,\sigma^2)\), with known \(\sigma=1\) outcome unit. The parameter \(\theta\) is the population mean in outcome units. The first analysis occurs after \(n_1=50\) observations and the final analysis after \(n_2=100\). At look 1, stop for efficacy if \(Z_1\geq 2.5\); otherwise continue to look 2 and stop there. For \(k=1,2\), let \(\bar X_k\) be the mean of the first \(n_k\) observations and \(Z_k=\sqrt{n_k}\bar X_k/\sigma\).

This setup assumes independent observations, the stated normal outcome model, and a prespecified rule applied to the interim statistic. It has no futility stopping, missing data, or change in the analysis plan. A simulation intended to represent another design must instead reproduce its actual data-generating assumptions and stopping rule.

Definition: Let \(T\) be the stopping stage, so \(T=1\) when \(Z_1\geq 2.5\), and \(T=2\) otherwise. For an estimator \(\hat\theta_T\), stage-specific bias at stage \(k\) is \(E(\hat\theta_T-\theta\mid T=k)\). Overall bias is \(E(\hat\theta_T-\theta)\), averaging over the stopping stages. For a reported interval \(C_T=[L_T,U_T]\), stage-specific coverage is \(P(\theta\in C_T\mid T=k)\), while overall coverage is \(P(\theta\in C_T)\). Interval length is \(U_T-L_T\).

For a concrete diagnostic, use the stopped sample mean \(\hat\theta_T=\bar X_T\) and the usual known-variance 95% interval \(C_T=\bar X_T\pm1.96\sigma/\sqrt{n_T}\). This interval is a candidate to evaluate, not an assumption that its coverage remains 95% after stopping. A stopping-aware interval can be assessed with the same simulation by replacing this formula with the specified interval procedure.

Simulation Design and a Useful Identity

In each replicate, generate 50 independent observations for the first look and 50 new independent observations for the second increment. Calculate \(Z_1\), apply the boundary, and generate or reveal the second increment only if the rule continues. At the stopping stage, record the estimate, interval endpoints, length, stage, and whether the interval contains the simulation value of \(\theta\). The simulation value of \(\theta\) is known because it is specified when generating data.

The decomposition below explains why stage-specific summaries and overall summaries must both be reported. It follows by partitioning the expectation according to the mutually exclusive stopping stages:

E(\hat\theta_T-\theta) =\sum_{k=1}^{2}P(T=k)E(\hat\theta_T-\theta\mid T=k).

The same partition applies to coverage probability and mean interval length. For example, overall coverage is a probability-weighted average of conditional coverage, not the unweighted average of the two stage-specific coverage estimates. The weights are the probabilities of reaching each stage under the chosen value of \(\theta\).

Specify the scenario.
Set \(\theta\), \(\sigma\), the look sizes, the stopping boundary, the estimator, and the interval procedure before generating replicates.

Generate and stop.
Generate outcomes under the model, calculate the interim statistic, and apply the prespecified boundary exactly. Continue only when the rule says to continue.

Record replicate-level results.
For each replicate save its stage, estimation error, coverage indicator, and interval length. Retaining these records allows checks and alternative summaries.

Summarize with Monte Carlo uncertainty.
Compute stage-specific summaries among replicates at that stage, then compute overall summaries across all replicates. Report Monte Carlo standard errors and stage counts.

Worked Examples

Worked Example: Applying the Stopping Rule to One Replicate

Suppose \(\theta=0.224\) outcome units. In one hypothetical replicate, the sum of the first 50 observations is 18.0, so \(\bar X_1=18.0/50=0.360\). The interim statistic is

Z_1=\frac{\sqrt{50}(0.360)}{1}=2.546.

Since \(2.546\geq2.5\), the replicate stops at stage 1. Its estimate is \(0.360\) outcome units. The stated candidate interval is

0.360\pm1.96\frac{1}{\sqrt{50}} =0.360\pm0.2772 =[0.0828,\ 0.6372].

Its length is \(0.6372-0.0828=0.5544\) outcome units, and it covers the generating value \(0.224\). This single covered interval contributes one to the coverage count, but it does not establish the coverage probability or the estimator's bias.

Worked Example: Computing the Overall Bias From Stage Results

Suppose a hypothetical simulation of \(100{,}000\) replicates gives stage 1 probability \(0.180\), stage-specific bias \(0.204\) outcome units at stage 1, and stage-specific bias \(-0.0450\) outcome units at stage 2. The implied overall bias is

0.180(0.204)+0.820(-0.0450) =0.03672-0.03690 =-0.000180\text{ outcome units}.

The two conditional biases are substantial relative to the nearly zero overall bias: early stopping selects unusually large interim estimates, while the stage-2 group consists of trials that did not cross the early efficacy boundary. The overall result does not mean there is no stage-specific selection effect. These are hypothetical simulation summaries, and their uncertainty must also be assessed.

Worked Example: Combining Coverage and Its Monte Carlo Error

Suppose the same hypothetical run has \(18{,}000\) stage-1 replicates with estimated coverage \(0.840\), and \(82{,}000\) stage-2 replicates with estimated coverage \(0.946\). The overall coverage estimate is

\widehat C =\frac{18{,}000(0.840)+82{,}000(0.946)}{100{,}000} =0.92692.

Thus, the estimated overall coverage is \(92.692\%\), below the nominal 95% for this candidate interval. The approximate Monte Carlo standard error for the overall proportion is

\sqrt{\frac{0.92692(1-0.92692)}{100{,}000}} =0.000823.

This is about \(0.0823\) percentage points. Conditional Monte Carlo standard errors are different because their denominators are the stage-specific counts: for stage 1, \(\sqrt{0.840(0.160)/18{,}000}=0.00273\), or about \(0.273\) percentage points. The stage-1 estimate is less precise because it is based on fewer replicates.

Deriving the Simulation Summaries

Let replicate \(r\) produce stopping stage \(T_r\), estimate \(\hat\theta_r\), interval \(C_r=[L_r,U_r]\), and coverage indicator \(I_r=1\{\theta\in C_r\}\), for \(r=1,\ldots,M\). If \(M_k\) replicates stop at stage \(k\), the natural conditional estimates are

\widehat{\operatorname{Bias}}_k =\frac{1}{M_k}\sum_{r:T_r=k}(\hat\theta_r-\theta),\qquad \widehat C_k=\frac{1}{M_k}\sum_{r:T_r=k}I_r,\qquad \bar L_k=\frac{1}{M_k}\sum_{r:T_r=k}(U_r-L_r).

Overall estimates use all \(M\) replicates: average the estimation errors for overall bias, the indicators for overall coverage, and the lengths for overall mean length. Equivalently, the overall summaries can be formed by weighting stage-specific summaries by \(M_k/M\). For bias, the identity follows directly because every replicate belongs to exactly one stage; the same argument applies to the coverage indicators and lengths.

The Monte Carlo standard error for a conditional bias estimate is \(s_{e,k}/\sqrt{M_k}\), where \(s_{e,k}\) is the sample standard deviation of replicate-level errors among stage-\(k\) replicates. For an estimated coverage proportion \(\widehat C_k\), use \(\sqrt{\widehat C_k(1-\widehat C_k)/M_k}\). The latter is the usual standard error for a proportion based on independent simulation replicates. These quantify simulation noise, not uncertainty from having observed only one real trial.

Diagnostics, Alternatives, and Limitations

Several checks help establish that the simulation is evaluating the intended procedure rather than a coding error:

Check stage frequencies. The counts must sum to \(M\). Compare the early-stop frequency across repeated runs or with its model-based value when it is straightforward to calculate.
Verify boundary logic. Confirm that \(Z_1=2.5\) is treated according to the stated greater-than-or-equal rule, and that stage-2 summaries include only trials that continued.
Inspect replicate records. Confirm that stage, sample size, estimate, interval endpoints, coverage, and length agree for a small number of hand-checked paths.
Monitor Monte Carlo precision. A small stage count can make conditional estimates noisy even when the total number of replicates is large. Increase \(M\) or report the resulting uncertainty rather than interpreting unstable digits.
Repeat with independent random-number streams. Large discrepancies between runs can signal inadequate simulation precision, an implementation problem, or both.

A tempting alternative is to evaluate the fixed-sample estimator or interval at \(n_1\) and \(n_2\) without selecting on the interim statistic. That comparison may help explain selection effects, but it is not a diagnostic of the stopped design: it omits the event that determines which estimate and interval are reported. A stopping-aware analytic calculation can provide a useful benchmark where available. A bootstrap assessment must also reproduce the stopping mechanism; resampling only the observed stopped data does not automatically reproduce the sampling law under the design.

Simulation results depend on the data-generating model, parameter values, and exact analysis rule. Normal outcomes with known variance do not establish performance for skewed outcomes, estimated variances, missing data, or a different boundary. Coverage is a probability over repeated trials at a fixed \(\theta\); it is not the probability that the parameter lies in a particular observed interval. Finally, an overall average can conceal poor performance at a rarely reached stage, so report stage counts and conditional summaries alongside overall results.

Common pitfall: Reporting a stage-1 conditional bias as the overall bias, or averaging stage-specific coverage estimates equally, ignores how often each stage occurs. Use the stage probabilities (or simulation proportions \(M_k/M\)) as weights; also retain the conditional results so that weighting does not hide stage-specific behavior.

Check Your Understanding

Answer each question using the definitions and simulation logic above.

If stage-specific biases are \(0.12\) and \(-0.02\) outcome units, and the stage probabilities are \(0.25\) and \(0.75\), what is overall bias?
Why is the stage-1 conditional coverage estimate not the same quantity as overall coverage?
A simulation reports \(M=50{,}000\), \(M_1=7{,}000\), and \(M_2=42{,}000\). Identify the diagnostic problem.
What is the approximate Monte Carlo standard error of an overall estimated coverage of \(0.90\) based on \(10{,}000\) independent replicates?
Why does a simulation that analyzes every replicate at the final sample size fail to assess the specified early-stopping design?

Solutions. 1. The overall bias is \(0.25(0.12)+0.75(-0.02)=0.015\) outcome units. 2. Stage-1 coverage conditions on stopping early, while overall coverage averages over both stopping stages in proportion to their probabilities. 3. The stage counts sum to \(49{,}000\), not \(50{,}000\); \(1{,}000\) replicates are unaccounted for, suggesting a recording or stage-classification error. 4. The standard error is \(\sqrt{0.90(0.10)/10{,}000}=0.00300\), or \(0.300\) percentage points. 5. It discards the specified rule that determines whether stage-1 results are reported, so it simulates a different sampling and reporting procedure.