As Frank Harrell says here and other places, it's better to compare two predictive models (Cox proportional hazards in this case) wrt discrimination (C-index) using the paired U-statistic (Hmisc::rcorrp.cens) rather than compare two C-indices. But this test can only really be done in independent validation data, since any overfitting makes the test non-informative (i.e., overfitting will necessarily improve discrimination in the same data). So this suggests splitting the data into two, training and testing.

But bootstrap is also advocated over data splitting, to avoid the big reduction in sample size from splitting and to get better precision in estimates.

So does it make sense to perform the paired test in bootstrap, and then report its bootstrapped quantiles for Somers' D? Any gotchas here?

[Edit: sample R code]




n <- 1000
d <- data.frame(time=rexp(n),
   event=sample(0:1, n, replace=TRUE),

cutoff <- 3
B <- 500
S.orig <- with(d, Surv(time, event))
n <- nrow(d)

# bootstrap with testing on samples not selected in current rep
res.boot <- foreach(b=1:B, .combine="c") %dopar% {
   cat("rep", b, "\n")
   s <- sample(n, replace=TRUE)
   ns <- (1:n)[!1:n %in% s]
   d.boot <- d[s, ]
   S.boot <- S.orig[s, ]
   fit1 <- cph(Surv(time, event) ~ Pred1, data=d.boot,
      x=TRUE, y=TRUE, surv=TRUE, time.inc=cutoff)
   fit2 <- cph(Surv(time, event) ~ Pred1 + Pred2, data=d.boot,
      x=TRUE, y=TRUE, surv=TRUE, time.inc=cutoff)

   pr1 <- survest(fit1, newdata=d[ns,], times=cutoff)
   pr2 <- survest(fit2, newdata=d[ns,], times=cutoff)

   rcorrp.cens(pr1$surv, pr2$surv, S.orig[ns, ])["Dxy"]

quantile(res.boot, prob=c(0.025, 0.975))

Nice work. You can save a bit of time by using cph(..., data=original.data, subset=s) instead of creating a new data frame that is a subset.

The gold standard test is the likelihood ratio $\chi^2$ test, and the next best test is the Wald test. If you want to use the rcorrp.cens "more concordant" test, it can likely still work if applied the naive way if the two models have the same amount of overfitting. To do it the right way in general, along the lines you have done, it may be better to use the optimism bootstrap as the rms package validate functions do. You estimate the drop-off in the special $D_{xy}$ statistic from rcorrp.cens when comparing performance of fits from bootstrap samples to the performance of the bootstrap-derived model on the original dataset.

This would give you the overfitting-corrected $D_{xy}$ but not the standard error needed for a hypothesis test.

The more ad hoc non-bootstrap approach of ignoring overfitting can possibly work because what is being assessed is relative concordance not absolute concordance. One might think that the amount of overfitting of the two models is the same one of the following are true:

  1. The two models were both pre-specified and estimating the same number of parameters.
  2. The two models were both developed from strategies that used the same number of candidate parameters.
  3. The heuristic shrinkage coefficients $\hat{\gamma} = \frac{LR - p}{LR}$ [where $LR$ is the likelihood ratio $\chi^2$ statistic for a model and $p$ is the model's effective degrees of freedom] are equal. $\hat{\gamma}$ estimates the calibration curve slope. The effective d.f. will be closer to the number of candidate parameters than to the number of parameters in the final model if stepwise selection was used.
  • $\begingroup$ Any how would you know if the two models had same amount of overfitting, assess via the optimism bootstrap for each one individually? $\endgroup$ – purple51 Oct 10 '14 at 5:31
  • $\begingroup$ See additional information in the edited answer. $\endgroup$ – Frank Harrell Oct 10 '14 at 12:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.