Covariate adjustment does not automatically buy power in a randomised survival trial: a plasmode resampling study of the NCCTG adjuvant colon cancer trial

Abstract

Prespecified adjustment for prognostic baseline covariates is routinely recommended to improve the efficiency of randomised trials. For the Cox model the argument is indirect: adjustment inflates the standard error of the treatment coefficient, but non-collapsibility inflates the conditional estimand by more, so the test is expected to gain power. That reasoning has been validated almost entirely by parametric simulation, in which covariates are balanced in expectation by construction. We tested it in a real trial's own resampling world. Using the 594 complete-case death records of the observation and levamisole-plus-fluorouracil arms of the NCCTG/Mayo adjuvant colon cancer trial, we built three plasmode data-generating mechanisms by bootstrapping the observed records: as-randomised, risk-balanced, and a permutation null. Over 4,000 replicates per cell at 150 and 300 patients per arm we compared the log-rank test, unadjusted Cox, Cox adjusted for eight prespecified covariates, Cox with within-replicate backward elimination, and restricted mean survival time (RMST) differences at 5 and 6.5 years. Adjustment lost power under the as-randomised mechanism, by 5.08 percentage points (Monte Carlo standard error, MCSE, 0.62) at 150 per arm and 2.70 points (0.43) at 300, because the trial's realised covariate imbalance cancelled the non-collapsibility inflation (estimand ratio 0.998) while the standard-error cost (ratio 1.035) remained. Resampling with baseline risk balanced restored the textbook gain (+9.55 points, MCSE 0.56). The adjusted Wald test was also anti-conservative (6.35% and 6.45% at alpha = 0.05). Trialists should treat the efficiency gain from Cox adjustment as a population property, not a guarantee for their trial.

Keywords: covariate adjustment; Cox proportional hazards model; non-collapsibility; plasmode simulation; restricted mean survival time; statistical power

1. Introduction

Randomisation makes a crude treatment comparison unbiased, but it does not make it efficient. Adjusting the primary analysis for strongly prognostic baseline covariates is therefore standard advice, and the empirical case for it in real trials is documented [1]. For linear models the case is simple: the estimand is unchanged and the residual variance falls. For non-linear models it is not. Robinson and Jewell showed for logistic regression that adjustment increases the standard error of the log odds ratio, yet increases the magnitude of the estimand by more, so that the test still gains [2]. Ford, Norrie and Ahmadi gave the survival analogue under the heading of model inconsistency: if conditional hazards are proportional, marginal hazards generally are not, so the adjusted and unadjusted Cox models cannot both be correctly specified and the marginal log hazard ratio is attenuated toward the null relative to the conditional one [3]. The unadjusted fit remains interpretable as a consistent estimator of a well-defined limiting parameter under this misspecification [4].

The resulting power claim has been quantified almost entirely by parametric Monte Carlo simulation. Hernández, Eijkemans and Steyerberg asked directly how much prespecified adjustment increases power for time-to-event outcomes [5], and Momal et al. report sample-size reductions from 3.1% to 85.7% depending on prognostic strength and cumulative incidence [6]. In every such design the covariates are balanced across arms in expectation by construction, and the realised imbalance in a replicate is a fresh random draw.

A real trial is not like that: its realised imbalance is a fixed feature of the dataset. Adjustment then does two things at once, inflating the estimand through non-collapsibility and correcting the crude estimate for that particular imbalance. The two can cancel. Whether the promised gain survives in a real trial is therefore an empirical question a parametric simulation cannot answer.

This paper answers it for one trial and explains the answer rather than merely reporting it. Our contributions are: (i) a plasmode design whose data-generating mechanism is a real randomised trial, with every estimand's reference value computed exactly inside the resampling population; (ii) a decomposition of the adjusted-versus-unadjusted power difference into an estimand ratio and a standard-error ratio, whose relative size fixes the sign of the difference; (iii) a paired evaluation, with discordance-based Monte Carlo standard errors, of five analysis strategies including within-replicate covariate selection and RMST contrasts.

2. Methods

Let $T_i$ denote the survival time of patient $i$, $C_i$ the censoring time, $Y_i = \min(T_i, C_i)$, $\delta_i = \mathbb{1}\{T_i \le C_i\}$, $A_i \in \{0,1\}$ the randomised arm and $Z_i \in \mathbb{R}^{p}$ a vector of baseline covariates. Two Cox models are considered,

whose treatment coefficients are the marginal and conditional log hazard ratios. Both are estimated by maximising the Breslow partial likelihood, with $\hat\beta$ solving $U(\beta) = \sum_i \delta_i \{x_i - \bar{x}(Y_i;\beta)\} = 0$ where $\bar{x}(t;\beta) = \sum_{j: Y_j \ge t} e^{x_j^{\top}\beta} x_j / \sum_{j: Y_j \ge t} e^{x_j^{\top}\beta}$, and model-based standard errors taken from the inverse observed information. Under non-collapsibility $|\beta^{\mathrm{C}}| > |\beta^{\mathrm{M}}|$ when $Z$ is prognostic and balanced [3]; under misspecification $\hat\beta^{\mathrm{M}}$ is still consistent for a well-defined functional [4].

Write $\kappa = |\beta^{\mathrm{C}}| / |\beta^{\mathrm{M}}|$ for the estimand ratio and $\rho = \mathrm{SE}(\hat\beta^{\mathrm{C}}) / \mathrm{SE}(\hat\beta^{\mathrm{M}})$ for the standard-error ratio. Since the Wald statistic is $|\hat\beta| / \mathrm{SE}(\hat\beta)$, adjustment increases power if and only if $\kappa > \rho$. Both quantities are estimable by resampling; the study measures them directly.

Six analyses were applied to every replicate: M1, the log-rank test with the hypergeometric variance; M2, the Wald test of $\hat\beta^{\mathrm{M}}$; M3, the Wald test of $\hat\beta^{\mathrm{C}}$ with all eight prespecified covariates; M4, the same after within-replicate backward elimination, dropping at each pass every covariate with Wald $p > 0.10$ and refitting (treatment always retained, at most four passes); and M5, M6, the difference in restricted mean survival time $\mathrm{RMST}(\tau) = \int_0^{\tau} \hat S(u)\,\mathrm{d}u$ from arm-specific Kaplan-Meier curves with the Greenwood-type variance $\sum_{t_i \le \tau} \{\int_{t_i}^{\tau} \hat S(u)\mathrm{d}u\}^2 d_i / \{n_i (n_i - d_i)\}$, at $\tau_1 = 5$ and $\tau_2 = 6.5$ years [7]. All tests are two-sided at $\alpha = 0.05$; all intervals are 95% Wald intervals. The Cox and log-rank routines were written in NumPy for speed and validated against statsmodels: coefficients and standard errors agreed with PHReg (Breslow ties) to $8.7\times10^{-13}$, and the log-rank statistic with survdiff to $2.1\times10^{-14}$.

3. Simulation Study

3.1 Aims

To determine whether prespecified Cox adjustment increases power over the unadjusted log-rank test in the resampling population of a real randomised trial; to decompose any difference into $\kappa$ and $\rho$; and to assess type I error, bias and coverage for the same strategies, within-replicate covariate selection and RMST contrasts.

3.2 Data-generating mechanisms

The data-generating mechanism is a real trial, not a parametric generator. We used the colon dataset of the R survival package, mirrored by the Rdatasets archive (https://vincentarelbundock.github.io/Rdatasets/csv/survival/colon.csv, retrieved 2026-09-08T15:47:32Z, sha256 6f3472a6...), which holds the patient records of the NCCTG/Mayo adjuvant chemotherapy trial in stage B/C colon cancer reported by Moertel et al. [8]. The file has 1,858 rows in long form, two per patient (etype = 1 recurrence, etype = 2 death), so 929 patients are represented. time is days from randomisation to event or censoring and status the 0/1 indicator; rx is the randomised arm (Obs, Lev, Lev+5FU); age is in years, nodes is the number of positive lymph nodes, differ tumour differentiation (1 well to 3 poor), extent depth of invasion (1 submucosa to 4 contiguous structures), and sex, obstruct, perfor, adhere, surg, node4 are 0/1 clinical indicators.

Filters, applied in this order: death records only (etype = 2), arms Obs and Lev+5FU only (619 patients), and complete cases on nodes and differ (25 patients dropped, 4.0%). The analysed set is 594 patients, 289 on Lev+5FU and 305 on observation, with 281 deaths (52.7% censored), median observed time to death or censoring 2,003.5 days and maximum 3,309 days (9.1 years). The prespecified covariate vector is $Z = \{$age, sex, obstruct, adhere, surg, $\log(\text{nodes}+1)$, differ, extent$\}$; perfor was excluded a priori as too rare for stable estimation, and node4 is a dichotomisation of nodes and so redundant.

Three mechanisms resample these 594 records with master seed 20260908 (Table 1). D1 (as-randomised) draws $n$ patients with replacement within each arm, preserving the real treatment effect, covariate joint distribution, censoring pattern and realised imbalance. D2 (risk-balanced) does the same but stratified on tertiles of the baseline prognostic index $\hat\gamma^{\top} Z$ from a covariate-only Cox fit, allocating both arms to the pooled tertile proportions. D3 (permutation null) draws $2n$ patients with replacement from the pooled sample and splits them at random into two arms of $n$, giving an exactly null treatment effect with the real covariate-outcome relationship intact. Per-arm sizes were $n = 150$ and $n = 300$; the latter reproduces the size of the trial itself.

D2 needs a caution that governs how the whole D1-to-D2 contrast should be read. Forcing both arms to the pooled tertile proportions changes each arm's covariate distribution, not merely the difference between the two, so D2 is a different resampling population from D1 and not D1 with chance imbalance switched off. The shift in the marginal estimand from $-0.3854$ under D1 to $-0.3243$ under D2 (Table 1) is a consequence of that change of population. Balancing is also only first-order: the index is one-dimensional and estimated from the same 594 records, so residual within-stratum imbalance, and imbalance in directions orthogonal to the index, survive.

Following the caution of Shaw et al. that a plasmode reference value must belong to the resampling population [9], each estimand's reference value was computed exactly, by fitting the same estimator to the 594 records with case weights equal to the per-record probability mass of that population, rather than by a large draw.

Table 1. The analysed trial sample, the three resampling populations and their exact estimand reference values.

Notes. Per-arm sizes $n$ = 150 and 300 in every mechanism; 4,000 replicates per cell, 24,000 in total, master seed 20260908. Reference values are exact population parameters, obtained by weighted Breslow Cox and weighted Kaplan-Meier fits to the 594 records, so they carry no Monte Carlo error. D1 and D2 are different populations, so their reference values are not directly comparable. The prognostic index uses the covariate-only Cox coefficients age 0.0023, sex -0.0991, obstruct 0.1094, adhere 0.2564, surg 0.2929, log(nodes+1) 0.7805, differ 0.1195, extent 0.3591; its mean is 0.072 lower (better) in the Lev+5FU arm, which is the realised chance imbalance. In the full sample 251 of the 281 deaths occur within 5 years; 178 and 153 patients are at risk at 5 years, 85 and 64 at 6.5 years.

3.3 Estimands

The marginal log hazard ratio $\beta^{\mathrm{M}}$ (M1, M2), the conditional log hazard ratio $\beta^{\mathrm{C}}$ (M3, approximately M4), and the RMST difference at $\tau_1$ and $\tau_2$ (M5, M6). These are different quantities: rejection rates are comparable across them, estimation performance is not.

3.4 Performance measures

Empirical rejection rate at $\alpha = 0.05$ (power under D1 and D2, type I error under D3); bias, empirical standard error, RMSE and 95% interval coverage against the exact reference; mean model-based standard error; mean number of covariates retained by M4; and non-convergence counts. Following Morris, White and Crowther [10], all methods were applied to the same replicate, so head-to-head differences are reported with a discordance-based Monte Carlo standard error, $\mathrm{MCSE} = \{n_{01} + n_{10} - (n_{10} - n_{01})^2 / B\}^{1/2} / B$, where $n_{10}$ and $n_{01}$ count replicates in which only one method rejects. With $B = 4{,}000$ the MCSE of a rejection rate is at most 0.79 points and that of a paired difference typically 0.4 to 0.8, which resolves differences of 2 points. No replicate failed to converge.

3.5 Results

Table 2. Empirical rejection rate (%) at $\alpha$ = 0.05, with Monte Carlo standard errors, from 4,000 replicates per cell.

Notes. D1 and D2 rows are power, D3 rows type I error. Rates are not comparable between D1 and D2, which are different populations. M4 retained on average 3.03 of 8 covariates at $n$ = 150 and 3.87 at $n$ = 300 under D1.

Table 3. Paired difference in rejection rate against the log-rank test (percentage points), with discordance-based Monte Carlo standard errors.

Notes. Bold entries are the study's primary comparison. Negative values favour the log-rank test. Differences are within-mechanism and paired on the replicate. Under D1 at $n$ = 150 there were 411 replicates in which only the log-rank test rejected against 208 in which only the adjusted Cox test did.

The primary hypothesis is refuted on both clauses. Under D1 the adjusted Cox test lost 5.08 points of power at 150 per arm and 2.70 at 300, each more than four MCSEs on the wrong side of zero. Under D3 it rejected 6.35% and 6.45% of the time against a nominal 5%, while the log-rank and unadjusted Cox tests were correctly sized (5.03% to 5.20%): the resampling machinery is calibrated, and the inflation sits in the adjusted model.

Table 4 records the mechanism of that inflation, though not its cause. Under D3 the adjusted estimator's mean model-based standard error is 6.9% below its empirical standard error (0.1724 against 0.1852 at $n$ = 150), against 1.9% for the unadjusted estimator, and coverage follows at 93.65% versus 94.98%. That discrepancy is the inflation restated rather than an explanation of it, and this study does not separate the two candidate causes: finite-sample bias of the observed information when eight nuisance parameters are estimated alongside the treatment effect, and the tie burden created by drawing 150 or 300 patients with replacement from 594 records, handled throughout by the Breslow approximation. All estimators are close to unbiased against their exact reference values (Tables 4 and 5); the largest absolute bias in a log hazard ratio anywhere in the study is 0.0068, for M4 under D2 at $n$ = 150, which is 1.8% of the estimand.

Table 4. Estimation performance for the treatment log hazard ratio: unadjusted (M2), adjusted (M3) and selected (M4) Cox, 4,000 replicates per cell.

Notes. Coverage is of the exact reference value by the 95% Wald interval; its MCSE is at most 0.36 points. M4's estimand changes from replicate to replicate, so its bias and coverage are measured against the full adjusted conditional estimand and are approximate. Under D2 the unadjusted estimator's empirical SE falls below its model SE, hence coverage above 95%: the stratified resample removes between-arm variation in baseline risk that the unstratified variance formula still charges for.

Table 5. Estimation performance for the restricted mean survival time difference in days, $\tau$ = 5 years (M5) and 6.5 years (M6), 4,000 replicates per cell.

Notes. All quantities are in days. The largest absolute bias is 2.20 days on a reference of 149.65. At $\tau$ = 5 years no replicate had either arm's follow-up end before the restriction time, so no Kaplan-Meier curve required extrapolation; the equivalent count was not recorded at 6.5 years. The over-coverage under D2 has the same origin as in Table 4.

Table 6. The two-ingredient decomposition: adjustment gains power if and only if the estimand ratio exceeds the standard-error ratio.

Notes. $\hat\kappa$ is the ratio of replicate-averaged adjusted to unadjusted estimates, $\hat\rho$ the ratio of their mean model-based standard errors; the exact population $\kappa$ is 0.998 (D1) and 1.175 (D2), and is undefined under the null. The sign of $\hat\kappa / \hat\rho - 1$ matches the sign of the observed difference in all four alternative cells; the apparent D3 "gain" of 1.3 points is size inflation, not power. Because D2 is a different population from D1 rather than D1 with imbalance removed, the rise in $\hat\kappa$ between the two rows should be read as the estimand ratio of a balanced population, not as a within-population effect of balancing.

Table 6 carries the mechanism. The standard-error cost of adjustment is real but small and nearly constant, $\hat\rho$ = 1.02 to 1.04 in every cell. The estimand ratio is what moves: in the trial as randomised there is no inflation at all ($\kappa$ = 0.998 exactly), because the Lev+5FU arm happened to receive patients with a mean prognostic index 0.072 lower than the observation arm, and correcting for that advantage pushes the adjusted estimate back toward the null by almost exactly as much as non-collapsibility pushes it away. In the balanced population D2, $\kappa$ is 1.175 and adjustment wins by 8 to 10 points. The theory of Ford et al. [3] is thus vindicated as a statement about populations and refuted as a guarantee about any one trial.

Two secondary findings. First, within-replicate backward elimination, retaining about 3 to 4 of the 8 covariates, did not inflate type I error beyond what full adjustment already costs (6.05% and 6.10% against 6.35% and 6.45%), and under D1 it recovered part of the lost power (60.13% against 58.35% at $n$ = 150) by discarding covariates that add variance without prognostic content; the common expectation that data-driven selection inflates size was not borne out here. Second, RMST contrasts were markedly less powerful than the log-rank test on this trial's empirical distribution, by 20.85 points at $\tau$ = 5 years and 12.13 at $\tau$ = 6.5 years ($n$ = 150, D1), even though 251 of the 281 deaths occur within five years; their type I error was nominal (4.80% to 5.30%). Because the deficit roughly halves as $\tau$ is extended, this is a cost of restriction rather than a defect of the estimand, consistent with Tian et al. [7] that the log-rank test is more efficient when hazards are approximately proportional.

4. Illustrative Example: the trial itself

This section is not an independent application: it analyses the same 594 patients that define the resampling populations, so it is the finite-sample realisation of D1 rather than an out-of-sample illustration. Its purpose is to show that the mechanism is visible in the single conventional analysis a trialist would report. Table 7 gives that analysis: crude and adjusted hazard ratios are indistinguishable (0.680 both times) while the adjusted standard error is 1.4% larger, so the adjusted analysis is the weaker of the two, exactly as the resampling study predicts for this population.

Table 7. Conventional analysis of the 594 analysed patients of the NCCTG adjuvant colon cancer trial (real data, overall survival).

Notes. Lev+5FU versus observation; 289 versus 305 patients, 281 deaths, 52.7% censored, median observed time to death or censoring 2,003.5 days. Adjustment covariates as in Section 3.2. Confidence intervals are Wald intervals on the log scale. The estimated hazard reduction of 32% is the trial's own overall-survival result [8].

5. Discussion

Findings. In the resampling world of a real randomised trial, prespecified Cox adjustment for a strongly prognostic covariate set reduced power by 5.08 points at 150 patients per arm and 2.70 points at 300, and produced a Wald test that rejected 6.4% of the time under a genuine null. In a balanced resampling population built from the same records, the same adjustment gained 8 to 10 points. Both results follow from one inequality, $\kappa > \rho$: the price of adjustment is a 2 to 4% larger standard error, and whether it is worth paying depends on whether non-collapsibility inflation survives the trial's realised imbalance.

Recommendations. Three follow from the scenarios studied. (i) The efficiency gain from Cox adjustment should be presented as an expectation over randomisations, not a property of a given trial: where the realised imbalance favours the experimental arm, adjustment can cost power, as it did here. (ii) The size inflation we measured for the adjusted Cox Wald test, 6.4% against a nominal 5%, is a caution about this design and not an established property of the test in a real trial of 150 to 300 patients per arm. It was observed under a with-replacement bootstrap whose replicates duplicate patients and therefore carry heavier ties than a fresh sample of the same size; all fits used the Breslow approximation, and Efron's approximation, the score test and the likelihood ratio test were not evaluated. It is a reason to check the size of an adjusted Wald test in one's own setting, not a general finding. (iii) Under approximately proportional hazards, an RMST contrast should be expected to be substantially less powerful than the log-rank test unless $\tau$ is set close to the end of follow-up.

Limitations. This is one trial, one endpoint and one covariate set; the sign of the D1 result depends on the realised imbalance of this particular randomisation, which these data cannot replicate. The bootstrap resamples with replacement, so replicates contain duplicated patients and more tied event times than a fresh sample would, and the resulting tie burden is confounded with every small-sample effect reported here. D2 removes only first-order imbalance, in tertiles of a one-dimensional index estimated from the same records, and changes the resampling population as Section 3.2 describes; it is an idealisation, not a re-randomisation. Twenty-five of 619 patients (4.0%) were dropped for missing nodes or differ with no imputation, so the resampling population is the complete-case population. Because the selected model's estimand varies across replicates, the bias and coverage of M4 are approximate. The treatment effect is fixed at the observed magnitude (HR 0.68) and the covariate set at eight variables, so the study yields single points rather than a curve, and the numerical size of the differences should not be extrapolated to trials with different effect sizes, event rates or prognostic strength. The RMST conclusions depend on the two restriction times chosen; at 6.5 years only 85 and 64 patients remain at risk in the full sample.

Extensions. Natural next steps are to repeat the design across a library of real trials, so that the distribution of $\kappa$ over randomisations is characterised rather than a single draw from it; to sweep effect size and prognostic strength within one trial's records; and to compare Efron's tie handling and robust or small-sample variance estimators for the adjusted Wald test, which would settle the cause of the size inflation left open here.

Data Availability and Software

The data are public. The colon dataset of the R survival package was retrieved from the Rdatasets archive at https://vincentarelbundock.github.io/Rdatasets/csv/survival/colon.csv on 2026-09-08 (sha256 6f3472a6..., 1,858 rows, 84,108 bytes), redistributed under LGPL (>= 2) with the package's COPYRIGHTS file recording "Copyright 2000 Mayo Foundation for Medical Education and Research"; the underlying trial is reported by Moertel et al. [8]. The analysis is one self-contained Python 3 script using NumPy, SciPy and statsmodels, with master seed 20260908; it loads the real file, applies the Section 3.2 filters and runs 24,000 replicates in under a minute. No data were simulated: every replicate is a bootstrap draw from the observed patient records.

References

1. Kahan BC, Jairath V, Doré CJ, Morris TP. The risks and rewards of covariate adjustment in randomized trials: an assessment of 12 outcomes from 8 studies. Trials. 2014;15:139. doi:10.1186/1745-6215-15-139
2. Robinson LD, Jewell NP. Some surprising results about covariate adjustment in logistic regression models. Int Stat Rev. 1991;59(2):227-240. doi:10.2307/1403444
3. Ford I, Norrie J, Ahmadi S. Model inconsistency, illustrated by the Cox proportional hazards model. Stat Med. 1995;14(8):735-746. doi:10.1002/sim.4780140804
4. Lin DY, Wei LJ. The robust inference for the Cox proportional hazards model. J Am Stat Assoc. 1989;84(408):1074-1078. doi:10.1080/01621459.1989.10478874
5. Hernández AV, Eijkemans MJC, Steyerberg EW. Randomized controlled trials with time-to-event outcomes: how much does prespecified covariate adjustment increase power? Ann Epidemiol. 2006;16(1):41-48. doi:10.1016/j.annepidem.2005.09.007
6. Momal R, Li H, Trichelair P, Blum MGB, Balazard F. More efficient and inclusive time-to-event trials with covariate adjustment: a simulation study. Trials. 2023;24:380. doi:10.1186/s13063-023-07375-0
7. Tian L, Fu H, Ruberg SJ, Uno H, Wei LJ. Efficiency of two sample tests via the restricted mean survival time for analyzing event time observations. Biometrics. 2018;74(2):694-702. doi:10.1111/biom.12770
8. Moertel CG, Fleming TR, Macdonald JS, et al. Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma. N Engl J Med. 1990;322(6):352-358. doi:10.1056/NEJM199002083220602
9. Shaw PA, Gruber S, Williamson BD, Desai R, Shortreed SM, Krakauer C, Nelson JC, van der Laan MJ. A cautionary note for plasmode simulation studies in the setting of causal inference. Stat Med. 2026;45(18-19). doi:10.1002/sim.70676
10. Morris TP, White IR, Crowther MJ. Using simulation studies to evaluate statistical methods. Stat Med. 2019;38(11):2074-2102. doi:10.1002/sim.8086