Why Simulate Power for Gray’s Test?
A sample-size calculation for Gray’s test must reflect more than the number of participants. Power also depends on how often each cause occurs, how the cause-specific hazards differ between groups, how many participants are censored, and how long they are followed. When those features do not fit a convenient approximation, simulation lets the design reproduce the planned study and estimate how often the test rejects under a specified alternative.
The result is conditional on the scenario supplied to the simulation. It is not a guarantee that a trial of that size will achieve the estimated power: the true hazards and censoring patterns may differ from the planning values. The earlier material on sensitivity analysis is useful for examining how that uncertainty affects a design; here the focus is constructing and using the simulation itself.
Specify the Data-Generating Model First
Suppose participants enter one of two groups, indexed by \(g\), and may experience the event of interest or one or more competing events. For participant \(i\) in group \(g\), let \(T_{ik}\) be the latent time to event type \(k\), with constant cause-specific hazard \(\lambda_{gk}\). Under the exponential model, the \(T_{ik}\) are generated independently as exponential random variables with rates \(\lambda_{gk}\). The observed event time is the earliest latent event time; its cause is the corresponding cause. This construction assumes the specified exponential event-time model and independent latent event times given group.
Let \(C_i\) denote a latent time to loss to follow-up. In the simple model below, loss times are independent exponential random variables with group-specific rate \(\gamma_g\), independent of the event times given group. The observed time is the minimum of the event time and the censoring time, and the outcome records whether an event occurred and which cause it was. If the study also has administrative censoring at a time \(\tau_i\), use the minimum of all three times: event, loss to follow-up, and administrative censoring.
These assumptions are modeling choices, not universal features of competing-risks data. The simulation should use the planned allocation, accrual, follow-up, administrative stopping rule, censoring mechanism, stratification, and test settings. For example, a uniform accrual period of length \(A\) followed by \(F\) additional years of follow-up gives a participant entering at time \(U_i\sim\text{Uniform}(0,A)\) administrative follow-up \(A+F-U_i\). If a different calendar schedule is planned, simulate that schedule instead.
Generate One Simulated Data Set
For exponential times, inverse-transform generation is direct. If \(U\) is uniformly distributed on \((0,1)\), then \(-\log(U)/\lambda\) has an exponential distribution with rate \(\lambda\). Generate one such time for each cause and, where appropriate, a loss-to-follow-up time. Take the minimum together with administrative censoring to obtain the record used by the test. Independent continuous times have probability zero of tying; software should nevertheless use a documented tie rule for ties created by rounding or discrete recording.
Worked Example: Generating One Participant’s Outcome
Consider a hypothetical participant in group 1 with cause-specific hazards \(\lambda_{11}=0.18\) and \(\lambda_{12}=0.12\) per year, and loss-to-follow-up rate \(\gamma_1=0.04\) per year. Suppose the generated uniform values for the two event times and loss time are \(0.70\), \(0.62\), and \(0.74\), respectively. There is no earlier administrative censoring in this example.
The generated times are \(\displaystyle T_{i1}=-\log(0.70)/0.18=1.981527\ldots\) years, \(\displaystyle T_{i2}=-\log(0.62)/0.12=3.983631\ldots\) years, and \(\displaystyle C_i=-\log(0.74)/0.04=7.527627\ldots\) years.
Rounded to three decimal places, these are \(1.982\), \(3.984\), and \(7.528\) years. The earliest time is \(T_{i1}\), so the simulated record is an event of cause 1 at \(1.982\) years. The competing event and loss-to-follow-up times are latent and are not recorded as later outcomes for this participant.
After generating all participants, apply Gray’s test to the simulated records exactly as planned for the study. For a two-group comparison of one cause-specific cumulative incidence function, the usual test statistic is compared with a chi-square distribution with one degree of freedom under the null hypothesis. At significance level \(\alpha=0.05\), reject when the test’s p-value is below \(0.05\), equivalently when the statistic exceeds approximately \(3.841\). If the design has more groups, strata, or a different prespecified testing rule, the simulation must use that rule rather than silently substituting this two-group example.
The simulated rejection indicator is 1 when the test rejects and 0 otherwise. Repeat the entire data-generation and testing process independently \(M\) times. If \(R\) of those simulated data sets yield rejection, the estimated power is \(\hat{\pi}=R/M\).
Monte Carlo Precision and Sample-Size Selection
Simulation adds Monte Carlo error: a different set of random draws would usually produce a slightly different \(\hat{\pi}\), even with the same sample size and model. Because each rejection indicator is a Bernoulli observation with success probability equal to the simulated power \(\pi\), the estimated Monte Carlo standard error is \(\displaystyle SE(\hat{\pi})=\sqrt{\hat{\pi}(1-\hat{\pi})/M}\). A rough 95% Monte Carlo interval is \(\hat{\pi}\pm1.96SE(\hat{\pi})\). It describes simulation uncertainty conditional on the assumed model; it does not include uncertainty about whether the planning hazards or censoring assumptions are correct.
Worked Example: Quantifying Monte Carlo Error
Suppose a hypothetical simulation uses \(M=2{,}000\) independent replicates for a candidate design. Gray’s test rejects in \(R=1{,}268\) replicates. The estimated power is \(\hat{\pi}=1{,}268/2{,}000=0.634\), or 63.4%.
Its estimated Monte Carlo standard error is \(\displaystyle SE(\hat{\pi})=\sqrt{0.634(1-0.634)/2{,}000}=\sqrt{0.000116022}\approx0.01077\). The rough 95% Monte Carlo interval is \(\displaystyle 0.634\pm1.96(0.01077)=(0.6129,0.6551)\).
Conditional on the simulation model and this candidate design, the estimated rejection probability is 63.4%, with Monte Carlo uncertainty of about 1.08 percentage points in standard-error units. The interval is not a confidence interval for the true clinical effect and does not account for model misspecification. More replicates would reduce this Monte Carlo uncertainty.
A useful planning approach is to evaluate a prespecified grid of candidate sample sizes, using the same allocation and other design rules at each size. Choose a target power in advance, such as 80% or 90%, and identify the smallest candidate that meets the target with adequate Monte Carlo precision. If a candidate is close to the target, increase \(M\) or simulate additional candidates near that point rather than treating a noisy estimate as definitive.
Worked Example: Choosing Among Candidate Sample Sizes
Consider a hypothetical two-group trial with equal allocation, a target cause of interest, cause-specific hazards \(\lambda_{11}=0.20\) and \(\lambda_{21}=0.30\) per year for that cause, and a competing-cause hazard of \(0.15\) per year in both groups. Assume independent loss to follow-up at rate \(0.05\) per year and administrative follow-up of three years for every participant. The planned analysis is a two-group Gray test at \(\alpha=0.05\). These values define a hypothetical planning scenario, not observations from a real study.
For orientation, with constant event hazards and no censoring, the cause-1 cumulative incidence by time \(t\) in a group with cause hazards \(\lambda_1,\lambda_2\) is \(\displaystyle F_1(t)=\int_0^t \lambda_1\exp\{-(\lambda_1+\lambda_2)u\}\,du =\frac{\lambda_1}{\lambda_1+\lambda_2}\left[1-\exp\{-(\lambda_1+\lambda_2)t\}\right]\). At three years, this gives approximately \(0.371\) in group 1 and \(0.494\) in group 2. These values describe the event process before censoring; the simulation still generates loss and applies the actual follow-up rule before running the test.
Suppose an illustrative simulation with \(M=2{,}000\) replicates at each candidate produces the following rejection counts. The counts and resulting powers are hypothetical outputs used to show the calculation.
| Participants per group | Rejections | Estimated power |
|---|---|---|
| 100 | 840 of 2,000 | 0.420 |
| 150 | 1,160 of 2,000 | 0.580 |
| 200 | 1,440 of 2,000 | 0.720 |
| 250 | 1,640 of 2,000 | 0.820 |
For the 250-per-group candidate, \(\hat{\pi}=1{,}640/2{,}000=0.820\), and \(\displaystyle SE(\hat{\pi})=\sqrt{0.820(0.180)/2{,}000}\approx0.00859\). Its rough 95% Monte Carlo interval is \(0.820\pm1.96(0.00859)=(0.8032,0.8368)\). In this illustrative grid, 250 participants per group is the smallest evaluated candidate with estimated power above 80%; the next smaller candidate has estimated power 72%. This supports selecting 250 per group under the stated assumptions, while recognizing that the calculation has not established performance under other hazards or censoring patterns.
Implementation Choices, Diagnostics, and Alternatives
A simulation is only as relevant as its correspondence to the proposed study. Before interpreting power, check that the simulated data and analysis match the design. In particular, verify that group allocation, accrual, calendar end, competing causes, loss to follow-up, administrative censoring, and any stratification are represented. Use the same target cause, test options, and significance level in every replicate. A practical diagnostic is to summarize simulated event counts and censoring by group and cause and compare them with the planning expectations. Implausible event counts, unexpectedly heavy censoring, or a sample-size curve that behaves erratically warrant investigation before the estimate is used.
The key assumptions deserve particular scrutiny. Exponential event times impose constant cause-specific hazards over time; that assumption may be unsuitable if the planned group differences vary with time. Independent loss to follow-up may be questionable if participants leave because of evolving health status. Administrative censoring must reflect the actual follow-up plan, not an arbitrary common duration when participants enter over a period. Simulation does not repair a misspecified model or an inappropriate test; it calculates consequences of the assumptions it is given.
Simulation is especially useful when unequal allocation, accrual, censoring, or a complicated hazard pattern makes a simple approximation difficult. An analytical approximation can be faster and can clarify how design quantities influence power, but may rely on simplifying assumptions. Simulation can incorporate more realistic rules, at the cost of Monte Carlo error and dependence on the specified model. The log-rank test is not an interchangeable shortcut when the design question concerns a cause-specific cumulative incidence function: the test should match the estimand and planned analysis, as explained in the earlier material on Gray’s test versus log-rank testing.
Do not use one power estimate to imply that the assumed treatment effect is known, or that a calculated sample size protects against every plausible data-generating process. The alternative hazards are planning inputs. Report them alongside the censoring and follow-up assumptions, candidate sample size, number of replicates, estimated power, and Monte Carlo precision. Sensitivity analyses over defensible alternative inputs help show whether the design conclusion is robust.
State the alternative cause-specific hazards, censoring, allocation, follow-up, test, and target power before running the calculation.
Generate competing event times and censoring for each participant, then record the earliest observed time and event type.
Run Gray’s test on each replicate using the prespecified analysis and count rejections at the chosen \(\alpha\).
Estimate power and its Monte Carlo error at each sample size, then choose a candidate that meets the target without relying on an imprecise boundary result.
Check Your Understanding
Answer each question using the definitions and procedures above. Brief solutions follow the questions.
- If each latent cause time and the loss time are generated independently, what determines the observed event type and time?
- A simulation has \(M=1{,}000\) replicates and rejects \(720\) times. Calculate estimated power and its Monte Carlo standard error.
- Why should a simulation with a three-year follow-up rule not simply censor every participant at three years if accrual continues for one year before the study closes?
- What does a Monte Carlo interval quantify, and what important uncertainty does it not quantify?
Solutions. 1. The observed record is determined by the minimum of the latent cause times and all applicable censoring times. If the earliest time is an event time, that cause is recorded; if a censoring time is earliest, the observation is censored. 2. The estimated power is \(720/1{,}000=0.720\). Its standard error is \(\sqrt{0.720(0.280)/1{,}000}=\sqrt{0.0002016}\approx0.01420\), or 1.42 percentage points. 3. Participants entering later have less calendar time before study closure. Assigning everyone three years would overstate follow-up for later entrants and fail to reproduce the planned administrative censoring. 4. A Monte Carlo interval describes random error from using a finite number of simulation replicates, conditional on the specified model. It does not quantify uncertainty about whether the hazards, censoring assumptions, or other planning inputs are correct.