IPCW-based time-dependent sensitivity/specificity metrics for survival models

[abichat]

Hi,

At the moment, roc_auc_survival() provides the time-dependent AUC at a chosen horizon t, which relies on time-dependent sensitivity and specificity (estimated with IPCW).

However, these underlying quantities are not directly exposed as user-facing metrics. It would be very useful to have dedicated functions like sens_survival(), spec_survival(), ppv_survival(), npv_survival()... that would return the IPCW-adjusted time-dependent metric.

This may be related to #537.

Thanks for considering this, and for your work on the tidymodels ecosystem!

[topepo]

We do have roc_curve_survival(). Would that be sufficient to obtain the four statistics you are interested in?

[abichat]

Thanks for the suggestion, I hadn’t considered using roc_curve_survival(), but I see two main limitations with this workaround:

• While it’s possible to derive NPV and PPV from sensitivity, specificity, and prevalence, computing the prevalence would require an additional auxiliary step outside the usual tidymodels workflow. At that point, I personally prefer relying directly on timeROC::SeSpPPVNPV, which handles these calculations in one call (but that’s my own preference).

• More importantly, this approach doesn’t allow computing these metrics on the fly when fitting multiple models with e.g. tune_grid(). And without dedicated yardstick functions for time-dependent sens/spec/PPV/NPV, it becomes impossible to use them as objective metrics during tuning.

Thanks again!

[topepo]

Would you like to use these values for plotting or diagnostics, or for optimization purposes? Do you need a yardstick function or another function that can perform the computations?

Keep in mind that the prevalence would be pretty complex to specify (since it is basically a prior that is time-dependent). We have the reverse KM curve, but it is unclear how one would specify or estimate these values.

We could make a yardstick function to compute those statistics en masse that would use roc_curve_survival() under the hood. Would you like to create a PR?

[abichat]

I my use cases, it could be both diagnostics and optimization purposes. I think if it's developed as a yardstick function, it will work for both, right?

To avoid having to explicitly specify a time-dependent prevalence, would it make sense to add a new prediction type that returns whether the event has occurred (0/1) at time t (and that could be integrated with tailor for cutoff choice)?. With such a prediction type (together with IPCW weights), computing time-dependent sens/spec/PPV/NPV essentially reduces to applying the usual binary metrics at time t, which could then be handled by yardstick in the standard way.

Below is a reprex showing how to currently compute this (three methods; the IPCW computation is slightly different in timeROC):

library(tidymodels)
library(survival)
library(censored)
library(timeROC)

veteran <- mutate(veteran, surv = Surv(time, status))

tpred <- 100
cutoff <- 0.5

coxph <- proportional_hazards(mode = "censored regression")
coxph_fitted <- fit(
  coxph,
  surv ~ age + karno,
  data = veteran
)

preds <-
  coxph_fitted %>%
  augment(new_data = veteran, eval_time = tpred) %>%
  unnest(.pred) %>%
  select(surv, time, status, .pred_survival, .weight_censored) %>% 
  mutate(
    risk_t = 1 - .pred_survival,
    pred_event = .pred_survival < cutoff,
    is_case = time <= tpred & status == 1,
    is_ctrl = time > tpred,
    w_case = if_else(is_case, .weight_censored, 0),
    w_ctrl = if_else(is_ctrl, .weight_censored, 0)
  )


pi_t <- sum(preds$w_case) / sum(preds$w_case + preds$w_ctrl)
perf_withprev <-
  preds %>%
  summarise(
    sens = sum(w_case * pred_event) / sum(w_case),
    spec = sum(w_ctrl * (1 - pred_event)) / sum(w_ctrl)
  ) %>%
  mutate(
    ppv = (sens * pi_t) / (sens * pi_t + (1 - spec) * (1 - pi_t)),
    npv = (spec * (1 - pi_t)) / (spec * (1 - pi_t) + (1 - sens) * pi_t)
  )

perf_tidymodels <-
  preds %>%
  mutate(
    pred_event = factor(pred_event, levels = c("TRUE", "FALSE")),
    is_case = factor(is_case, levels = c("TRUE", "FALSE"))
  ) %>%
  summarise(
    sens = sens_vec(
      truth = is_case,
      estimate = pred_event,
      case_weights = .weight_censored
    ),
    spec = spec_vec(
      truth = is_case,
      estimate = pred_event,
      case_weights = .weight_censored
    ),
    ppv = ppv_vec(
      truth = is_case,
      estimate = pred_event,
      case_weights = .weight_censored
    ),
    npv = npv_vec(
      truth = is_case,
      estimate = pred_event,
      case_weights = .weight_censored
    )
  )

res_troc <-
  SeSpPPVNPV(
    cutpoint = cutoff,
    T = preds$time,
    delta = preds$status,
    marker = preds$risk_t,
    cause = 1,
    times = tpred
  )

perf_troc <-
  c(sens = "TP", spec = "FP", ppv = "PPV", npv = "NPV") %>%
  map(~ pluck(res_troc, ., 2)) %>%
  as.data.frame() %>%
  as_tibble()
perf_troc$spec <- 1 - perf_troc$spec


bind_rows(
  withprev = perf_withprev,
  tidymodels = perf_tidymodels,
  troc = perf_troc,
  .id = "method"
) %>% 
  column_to_rownames("method")
#>                 sens      spec       ppv       npv
#> withprev   0.7696859 0.6792453 0.7729439 0.6752077
#> tidymodels 0.7696859 0.6792453 0.7729439 0.6752077
#> troc       0.7665815 0.6792453 0.7665419 0.6792935

^{Created on 2025-11-17 with reprex v2.1.1}

[topepo]

I think that we can do that, with some internal modifications, so long as there is a single probability cutoff value. Otherwise, we would have to modify the structure of metric sets, which has implications for how we summarize resamples, choose, and optimization goals for racing/iterative optimization, and so on.

Is that an acceptable constraint?

[abichat]

As long as the single cutoff value can be selected (and potentially tuned), this works perfectly well for my use case. Thanks!

[topepo]

Selected yes. The interface is a little kludgy if you want to modify the threshold. The CCC example shows how to do this. It does enable you to look at multiple thresholds but you would have to make a new cloned metric for each threshold (and differentiate them by name).

If you mean can I optimize my tning parameters to maximize the sensitivity when the probability threshold is X then yes.

[abichat]

Oh, I see what you meant.

This constraint arises when treating the survival probability at time t as the prediction and computing the metric directly from it. In that setting, the cutoff has to be embedded inside the metric definition.

But if we instead consider a predicted event indicator at time t (which doesn’t exist yet for survival models but does exist for classification—and is exactly what sens(), spec()... operate on), then the cutoff can be selected and tuned externally, for example via tailor::adjust_probability_threshold.
With such a prediction type, one would be able to choose the cutoff that maximizes a given metric (e.g. sensitivity) without having to clone metrics for each threshold.

[topepo]

I see your point about that. A new dynamic prediction type a good idea but I'm not sure that we are able to add one right now ("now" is probably the next 2-4 months).

We do make that assignment in yardstick code using the system defined by Graf. One numerical issue is that there will be a lot of missing outcome designations due to ambiguities (Graf's three cases). The rate of missingness will change with the threshold (which is the reason for the propensity weights).

the cutoff can be selected and tuned externally,

That would be an interesting tailor method.

Looking at a lot of the time-dependent curves, I feel like it is making the best of a bad situation. There's usually a high class imbalance at later time points, and the probability statement we would make is complex due to the missingness and propensity weights.

If we were earlier (or later) in the year, this would make a great summer internship project. Our short- and medium-term path is fairly set right now.

The good news is that there are tools that can let you do this by saving the predictions. Integration will take a bit longer.

[abichat]

I’m not in a hurry, for now I have a (somewhat verbose) workaround. And I understand that creating a new prediction type and the associated metrics can be challenging. Still, I believe this could become a useful feature for other users as well.

Thanks again for taking the time to walk through all the details.