Supplementary Material

Supporting information for paper: ‘When Are Scoring Rules Proper? Bridging Theory and Practice in Survival Model Evaluation’

Load libraries:

library(survival)
library(dplyr)
library(tibble)
library(ggplot2)
library(patchwork)
library(akima)
library(purrr)
library(tidyr)
library(data.table)
library(DT)

Weibull distribution variability

To illustrate the diversity of the Weibull distributions used in our simulation design, we plot both the density f(t) and survival S(t) functions for 10 randomly sampled shape \alpha and scale \sigma parameter pairs. These distributions represent either the true survival, predicted survival, or censoring times across simulations. The wide parameter range induces substantial variability, ensuring a rich family of distributions and increasing the sensitivity of our evaluation to detect violations of properness. By construction, all data-generating mechanisms assume Y independent of C, satisfying the random censoring assumption.

The figure below corresponds to Figure 1 in the main manuscript.

set.seed(0L)

shape = runif(9, 0.5, 5)
scale = runif(9, 0.5, 5)

time_grid = seq(0, 6, 0.1)

params = tibble(
  id = seq_along(shape),
  shape = shape,
  scale = scale,
  group = sprintf("(%0.2f, %0.2f)", shape, scale)
)

grid = tidyr::expand_grid(params, t = time_grid)

dw = grid |>
  mutate(value = dweibull(t, shape, scale), type = "Density f(t)")

pw = grid |>
  mutate(value = pweibull(t, shape, scale, lower.tail = FALSE), type = "Survival S(t)")

plot_df = bind_rows(dw, pw) |>
  mutate(type = factor(type, levels = c("Density f(t)", "Survival S(t)")))

ggplot(plot_df, aes(x = t, y = value, color = group, group = group)) +
  geom_line(linewidth = 0.8, alpha = 0.9) +
  facet_wrap(~ type, ncol = 1, scales = "free_y") +
  scale_color_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 0:6) +
  labs(
    x = expression("Time"),
    y = NULL,
    color = expression("(" * alpha * "," * sigma * ")"),
    title = "Weibull Density and Survival Curves"
  ) +
  theme_bw(base_size = 16) +
  theme(
    legend.position = "right",
    strip.text = element_text(face = "bold")
  )

Weibull density (top) and survival (bottom) curves for randomly generated shape, \alpha, and scale, \sigma parameters. Each curve corresponds to a different parameter pair, illustrating the broad variability across simulated distributions.

#ggsave("weibull_plots.png", device = png, dpi = 600, width = 8, height = 5)

Table 1

Note

properness_test.R is the core file related to the assessment of empirical violations of properness.

Example execution: Rscript properness_test.R 10000 1000 10000 FALSE 42.

The above arguments correspond to (in that order): K = 10\,000 simulations, m = 1\,000 distributions per simulation and n = 10\,000 samples drawn from each distribution. FALSE means that S_C(t) is used directly from the true Weibull distribution and not estimated, 42 is just a seed number for reproducibility. There is a hard-coded argument in the script, num_cores = 60, to allow parallelization of the simulations across multiple CPU cores. Change as needed for your system.

To run for different sample sizes n (sequentially), use run_tests.sh.

To get the results displayed, we executed: ./run_tests.sh 10 25 50 100 250 500 1000 2500 5000 10000 (this will take a while to run, depending on your system).

Merged simulation results for all sample sizes are available here.

Results from 10 randomly sampled simulations:

res = readRDS("results/res_sims10000_distrs1000_0.rds")

res[sample(1:nrow(res), 10), ] |>
  DT::datatable(
    rownames = FALSE,
    options = list(dom = "t", searching = FALSE, lengthChange = FALSE)
  ) |>
  formatRound(columns = 3:40, digits = c(rep(3,8), rep(5,30)))

Each row corresponds to one simulation, and the quantites are averaged across m = 1000 draws from Weibull survival distributions at a fixed sample size n \in \{10, \dots, 10\,000\}.

The output columns are:

• sim: simulation index (k in the paper)
• n: sample size
• n_events: number of uncensored events observed in the simulated dataset
• max_t: maximum observed follow-up time
• m_true, m_pred: true and predicted mean survival times from the generating Weibull models
• m_true_est, m_pred_est: estimated mean survival times using the trapezoidal rule on the observed time point support
• S_Y_q10, S_Y_med, S_Y_q90: values of the true event-time survival function S_Y(t) evaluated at the 10th percentile, median, and 90th percentile of the observed time points.
• S_C_q10, S_C_med, S_C_q90: Values of the censoring survival function S_C(t) evaluated at the 10th percentile, median, and 90th percentile of the observed time points.
• S_Y_tail, S_C_tail: tail survival probability at max_t for event and censoring distributions
• epsilon: \hat\varepsilon = \frac{\sum_{d=1}^{1000} \left[S_Y^d(t_{\\max}^d) \times S_C^d(t_{\max}^d)\right]}{1000}
• {surv|pred|cens|}_{scale|shape} the scale and shape of the Weibull distributions for true, predicted and censoring times respectively
• prop_cens: proportion of censoring in each simulation
• tv_dist: the total variation distance between the true and the predicted Weibull distribution (closer to 0 means more similar distributions)
• {score}_diff, {score}_sd: the mean and standard deviation of the score difference D between true and predicted distributions across draws

Scoring rules analyzed:

• The SBS at the 10th, 50th and 90th percentiles of observed times
• The ISBS, using 50 equidistant points between the 5th and 80th percentile of observed times
• The RCLL
• An experimental weighted version of RCLL, which we called wRCLL (or RCLL*), see weighted_RCLL.R for more details. This is not a proper scoring rule in general (at least in theory), and was included in initial versions of this work as a possible (even more proper) extension of RCLL.

Important

We compute 95% confidence intervals (CIs) for the mean score difference \bar{D} (true - predicted) using a t-distribution. A violation is marked statistically significant if:

• The mean score difference exceeds a positive threshold: \bar{D} > 0.0001
• The CI excludes zero (i.e., \text{CI}_\text{lower} > 0)

res = res |>
  select(!matches("shape|scale|prop_cens|tv_dist|n_events|max_t|m_true|m_pred"))

measures = c("SBS_q10", "SBS_median", "SBS_q90", "ISBS", "RCLL", "wRCLL")
n_distrs = 1000 # `m` value from paper's experiment
tcrit = qt(0.975, df = n_distrs - 1) # 95% t-test
threshold = 1e-4

data = measures |>
  lapply(function(m) {
    mean_diff_col = paste0(m, "_diff")
    mean_sd_col = paste0(m, "_sd")

    res |>
      select(sim, n, !!mean_diff_col, !!mean_sd_col, starts_with("shift_"),
             epsilon, starts_with("S_")) |>
      select(-ends_with("tail")) |>
      mutate(
        se = .data[[mean_sd_col]] / sqrt(n_distrs),
        CI_lower = .data[[mean_diff_col]] - tcrit * se,
        CI_upper = .data[[mean_diff_col]] + tcrit * se,
        signif_violation = .data[[mean_diff_col]] > threshold & CI_lower > 0,
        # Keep relevant SBS columns
        shift_q10 = if (m == "SBS_q10") shift_q10 else NA_real_,
        S_Y_q10   = if (m == "SBS_q10") S_Y_q10 else NA_real_,
        S_C_q10   = if (m == "SBS_q10") S_C_q10 else NA_real_,
        shift_med = if (m == "SBS_median") shift_med else NA_real_,
        S_Y_med   = if (m == "SBS_median") S_Y_med else NA_real_,
        S_C_med   = if (m == "SBS_median") S_C_med else NA_real_,
        shift_q90 = if (m == "SBS_q90") shift_q90 else NA_real_,
        S_Y_q90   = if (m == "SBS_q90") S_Y_q90 else NA_real_,
        S_C_q90   = if (m == "SBS_q90") S_C_q90 else NA_real_
      ) |>
      select(!se) |>
      mutate(metric = !!m) |>
      relocate(metric, .after = 1) |>
      rename(diff = !!mean_diff_col,
             sd   = !!mean_sd_col)
  }) |>
  bind_rows()

data$metric = factor(
  data$metric,
  levels = measures,
  labels = c("SBS (Q10)", "SBS (Median)", "SBS (Q90)", "ISBS", "RCLL", "RCLL*")
)

For example, the results for the first simulation and n = 250 across all scoring rules show no violations of properness:

data |>
  filter(sim == 1, n == 250) |>
  select(-starts_with("shift"), -epsilon, -starts_with("S_")) |>
  DT::datatable(
    rownames = FALSE,
    options = list(dom = "t", searching = FALSE, lengthChange = FALSE)
  ) |>
  formatRound(columns = 4:7, digits = 5)

We summarize violations across simulations (sim), by computing for each score & sample size:

• Number of significant violations
• Violation rate
• Average, median, min and max score difference among simulations where violations occurred

all_stats =
  data |>
  group_by(n, metric) |>
  summarize(
    n_violations = sum(signif_violation),
    violation_rate = mean(signif_violation),
    diff_mean = if (any(signif_violation)) mean(diff[signif_violation]) else NA_real_,
    diff_median = if (any(signif_violation)) median(diff[signif_violation]) else NA_real_,
    diff_min = if (any(signif_violation)) min(diff[signif_violation]) else NA_real_,
    diff_max = if (any(signif_violation)) max(diff[signif_violation]) else NA_real_,
    .groups = "drop"
  )

all_stats |>
  arrange(metric) |>
  DT::datatable(
    rownames = FALSE,
    options = list(searching = FALSE),
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: center; font-size:150%",
      "Table 1: Empirical violations of properness"
    )
  ) |>
  formatRound(columns = 4:8, digits = 5)

TipSummary

• RCLL and RCLL* show practically no violations across any simulation or sample size.
• SBS shows time- and sample-size-dependent violations, mostly at smaller n, especially for early evaluation times (\tau^* = q_{0.1}). All differences are small (typically < 0.01).
• ISBS showed minor violations only at n \le 50; none for n > 100.

SBS Improperness: Bias and Tail Effects

Here we use the results from the previous section. The theoretical shift (i.e. bias) of the global minimum x^* of the expected SBS from the true survival probability at \tau^* is (for \varepsilon \ge 0):

x^\star - S_Y(\tau^*) = - \frac{\varepsilon (1-S_Y(\tau^*))}{S_C(\tau^*) - \varepsilon} \approx - \varepsilon \cdot \frac{1 - S_Y(\tau^*)}{S_C(\tau^*)}

Filter SBS simulation data for plotting:

plot_data = data |>
  filter(grepl("^SBS", metric)) |>
  pivot_longer(
    cols = starts_with("shift_"),
    names_to = "shift_type",
    values_to = "shift",
    values_drop_na = TRUE  # only keep the relevant shift per metric
  ) |>
  mutate(
    S_Y = case_when(
      metric == "SBS (Q10)" ~ S_Y_q10,
      metric == "SBS (Median)" ~ S_Y_med,
      metric == "SBS (Q90)" ~ S_Y_q90
    ),
    S_C = case_when(
      metric == "SBS (Q10)" ~ S_C_q10,
      metric == "SBS (Median)" ~ S_C_med,
      metric == "SBS (Q90)" ~ S_C_q90
    )
  ) |>
  select(-S_Y_q10, -S_Y_med, -S_Y_q90, -S_C_q10, -S_C_med, -S_C_q90) |>
  mutate(
    metric = recode(
      metric,
      "SBS (Q10)"   = "SBS (Early, q10)",
      "SBS (Median)"= "SBS (Median, q50)",
      "SBS (Q90)"   = "SBS (Late, q90)"
    )
  )

• The \hat \varepsilon vs bias plot corresponds to Figure 2 in the main manuscript.
• The Bias surface over (S_Y,S_C) (n = 50) corresponds to Figure 3 in the main manuscript.

\hat \varepsilon vs bias

ann_df = data.frame(
  metric = unique(plot_data$metric)[1],
  x = 0.01,
  y = -0.07,
  xend = 0.01,
  yend = -0.175,
  label = "larger bias\n(improperness)"
)

# shift vs epsilon
set.seed(42)
plot_data |>
  slice_sample(prop = 0.2) |> # 20% of: 3 SBS metrics x 10 sample sizes (n) x 10000 simulations per n
  ggplot(aes(x = epsilon, y = shift, color = factor(n))) +
    geom_point(alpha = 0.5, size = 0.1) +
    facet_wrap(~metric, scales = "fixed") +
    scale_color_discrete(
      name = "Sample size (n)",
    ) +
    scale_x_continuous(
      limits = c(0, 0.11),
      breaks = c(0, 0.025, 0.05, 0.075, 0.1)
    ) +
    labs(
      x = "Tail product estimate ε̂ ", #expression(hat(epsilon) * " "),
      y = expression(x^"*" ~ "-" ~ S[Y](tau^ ~ symbol("\052"))), # shift from optimal prediction
      color = "Sample size (n)"
    ) +
    theme_bw(base_size = 14, base_family = "Arial") +
    #theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
    guides(color = guide_legend(override.aes = list(size = 4))) +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "black") +
    geom_segment(
      data = ann_df,
      aes(x = x, y = y, xend = xend, yend = yend),
      inherit.aes = FALSE,linewidth = 1,
      arrow = arrow(length = unit(0.25, "cm"))
    ) +
    geom_text(
      data = ann_df,
      aes(x = x + 0.02, y = y - 0.05, label = label),
      inherit.aes = FALSE,
      angle = 90,
      size = 6,
      fontface = "bold"
    )

Relationship between survival tail product estimate \hat{\varepsilon} = \frac{\sum_{d=1}^{1000} \left[S_Y^d(t_{\max}^d) \times S_C^d(t_{\max}^d)\right]}{1000} and theoretical shift from the global minimum, using a sample of data from the Weibull simulation experiments

#ggsave(filename = "epsilon_vs_bias.png", dpi = 600, height = 5, width = 10, units = "in")

S_Y(\tau^*) vs bias

set.seed(42)
plot_data |>
  slice_sample(prop = 0.2) |> # 3 SBS metrics x 10 sample sizes (n) x 10000 simulations per n
  ggplot(aes(x = S_Y, y = shift, color = factor(n))) +
  geom_point(alpha = 0.5, size = 0.3) +
  facet_wrap(~metric, scales = "fixed") +
  scale_color_discrete(name = "Sample size (n)") +
  scale_x_continuous(
    limits = c(0, 1),
    breaks = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  labs(
    x = expression(S[Y](tau^ ~ symbol("\052"))),
    y = expression(x^"*" ~ "-" ~ S[Y](tau^ ~ symbol("\052"))),
  ) +
  theme_bw(base_size = 14) +
  theme(panel.spacing = unit(1, "lines")) +
  guides(color = guide_legend(override.aes = list(size = 4))) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "black") +
  geom_segment(
    data = ann_df,
    aes(x = x, y = y, xend = xend, yend = yend),
    inherit.aes = FALSE,linewidth = 1,
    arrow = arrow(length = unit(0.25, "cm"))
  ) +
  geom_text(
    data = ann_df,
    aes(x = x + 0.17, y = y - 0.05, label = label),
    inherit.aes = FALSE,
    angle = 90,
    size = 6,
    fontface = "bold"
  )

Relationship between true survival S_Y(\tau^*) and empirical shift from the global minimum, using a sample of data from the Weibull simulation experiments

S_C(\tau^*) vs bias

set.seed(42)
plot_data |>
  slice_sample(prop = 0.2) |> # 3 SBS metrics x 10 sample sizes (n) x 10000 simulations per n
  ggplot(aes(x = S_C, y = shift, color = factor(n))) +
  geom_point(alpha = 0.5, size = 0.3) +
  facet_wrap(~metric, scales = "fixed") +
  scale_color_discrete(name = "Sample size (n)") +
  scale_x_continuous(
    limits = c(0, 1),
    breaks = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  labs(
    x = expression(S[C](tau^ ~ symbol("\052"))),
    y = expression(x^"*" ~ "-" ~ S[Y](tau^ ~ symbol("\052"))),
  ) +
  theme_bw(base_size = 14) +
  theme(panel.spacing = unit(1, "lines")) +
  guides(color = guide_legend(override.aes = list(size = 4))) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "black") +
  geom_segment(
    data = ann_df,
    aes(x = x, y = y, xend = xend, yend = yend),
    inherit.aes = FALSE,linewidth = 1,
    arrow = arrow(length = unit(0.25, "cm"))
  ) +
  geom_text(
    data = ann_df,
    aes(x = x + 0.17, y = y - 0.05, label = label),
    inherit.aes = FALSE,
    angle = 90,
    size = 6,
    fontface = "bold"
  )

Relationship between censoring survival S_C(\tau^*) and empirical shift from the global minimum, using a sample of data from the Weibull simulation experiments

Bias surface over (S_Y,S_C) (n = 50)

base_data = plot_data |>
  dplyr::filter(n == 50)

interp_shift = function(df, nx = 150, ny = 150) {
  if (nrow(df) == 0) return(NULL)

  interp_res = with(df,
    akima::interp(x = S_C, y = S_Y, z = shift, nx = nx, ny = ny, duplicate = "mean")
  )

  expand.grid(S_C = interp_res$x, S_Y = interp_res$y) |>
    mutate(shift = as.vector(interp_res$z),
           metric = unique(df$metric))
}

interp_list = base_data |>
  dplyr::group_split(metric, .keep = TRUE) |>
  purrr::map(.f = interp_shift)

surface_data = bind_rows(interp_list) |>
  dplyr::filter(!is.na(shift))

range_shift = range(surface_data$shift, na.rm = TRUE)

min_shift = range_shift[1] # most negative
max_shift = range_shift[2] # closest to zero

ggplot(surface_data, aes(x = S_C, y = S_Y, fill = shift)) +
  geom_tile() +
  facet_wrap(~ metric) +
  scale_fill_viridis_c(
    direction = -1,
    limits = c(min_shift, 0),
    name = expression(x^"*" ~ "-" ~ S[Y](tau^ ~ symbol("\052")))
  ) +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(
    x = expression(S[C](tau^ ~ symbol("\052"))),
    y = expression(S[Y](tau^ ~ symbol("\052")))
  ) +
  theme_bw(base_size = 14) +
  theme(panel.spacing = unit(1, "lines"))

Interpolated surface of empirical shift x^\star - S_Y(\tau^*) as a function of S_Y(\tau^*) and S_C(\tau^*) for n = 50 (with \hat{\varepsilon} ~ 0.02). Panels correspond to early, median, and late evaluation times.

#ggsave(filename = "bias_surface_n50.png", dpi = 600, height = 5, width = 10, units = "in")

Bias surface over (S_Y,S_C) (n = 100)

base_data2 = plot_data |>
  dplyr::filter(n == 100)

interp_list2 = base_data2 |>
  dplyr::group_split(metric, .keep = TRUE) |>
  purrr::map(.f = interp_shift)

surface_data2 = bind_rows(interp_list2) |>
  dplyr::filter(!is.na(shift))

range_shift2 = range(surface_data2$shift, na.rm = TRUE)

min_shift = range_shift2[1] # most negative
max_shift = range_shift2[2] # closest to zero

ggplot(surface_data2, aes(x = S_C, y = S_Y, fill = shift)) +
  geom_tile() +
  facet_wrap(~ metric) +
  scale_fill_viridis_c(
    direction = -1,
    limits = c(min_shift, 0),
    name = expression(x^"*" ~ "-" ~ S[Y](tau^ ~ symbol("\052")))
  ) +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(
    x = expression(S[C](tau^ ~ symbol("\052"))),
    y = expression(S[Y](tau^ ~ symbol("\052")))
  ) +
  theme_bw(base_size = 14) +
  theme(panel.spacing = unit(1, "lines"))

Interpolated surface of empirical shift x^\star - S_Y(\tau^*) as a function of S_Y(\tau^*) and S_C(\tau^*) for n = 100 (with \hat{\varepsilon} ~ 0.01). Panels correspond to early, median, and late evaluation times.

#ggsave(filename = "bias_surface_n100.png", dpi = 600, height = 5, width = 10, units = "in")

plot_data |>
  group_by(n) |>
  summarise(
    mean_hat_epsilon = mean(epsilon)
  ) |>
  DT::datatable(
    rownames = FALSE,
    width = "52%",
    options = list(dom = "t", searching = FALSE, lengthChange = FALSE),
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: center; font-size:150%",
      "Empirical tail product estimate vs sample size"
    )
  ) |>
  formatRound(columns = 2, digits = 4)

SBS Improperness: administrative censoring

Note

Here we investigate the impact of administritive censoring on SBS properness. Administritive censoring is applied at the 80% percentile of true event times. SBS is calculated at the 10th, 50th and 90th percentile of observed times (after application of the administrative censoring).

Use properness_test_admin_cens.R to run the experiment.

Execute: Rscript properness_test_admin_cens.R 10000 1000 10000 FALSE 42 to run for n = 10\,000.

Simulation results are available here.

res = readRDS("results/res_sims10000_distrs1000_n10000_0_admin_cens.rds") |>
  select(!matches("shape|scale|prop_cens|tv_dist|n_events|max_t|m_true|m_pred"))

measures = c("SBS_q10", "SBS_median", "SBS_q90", "ISBS")
n_distrs = 1000 # `m` value from paper's experiment
tcrit = qt(0.975, df = n_distrs - 1) # 95% t-test
threshold = 1e-4

data = measures |>
  lapply(function(m) {
    mean_diff_col = paste0(m, "_diff")
    mean_sd_col = paste0(m, "_sd")

    res |>
      select(sim, n, !!mean_diff_col, !!mean_sd_col, starts_with("shift_"),
             epsilon, starts_with("S_")) |>
      select(-ends_with("tail")) |>
      mutate(
        se = .data[[mean_sd_col]] / sqrt(n_distrs),
        CI_lower = .data[[mean_diff_col]] - tcrit * se,
        CI_upper = .data[[mean_diff_col]] + tcrit * se,
        signif_violation = .data[[mean_diff_col]] > threshold & CI_lower > 0,
        # Keep relevant SBS columns
        shift_q10 = if (m == "SBS_q10") shift_q10 else NA_real_,
        S_Y_q10   = if (m == "SBS_q10") S_Y_q10 else NA_real_,
        S_C_q10   = if (m == "SBS_q10") S_C_q10 else NA_real_,
        shift_med = if (m == "SBS_median") shift_med else NA_real_,
        S_Y_med   = if (m == "SBS_median") S_Y_med else NA_real_,
        S_C_med   = if (m == "SBS_median") S_C_med else NA_real_,
        shift_q90 = if (m == "SBS_q90") shift_q90 else NA_real_,
        S_Y_q90   = if (m == "SBS_q90") S_Y_q90 else NA_real_,
        S_C_q90   = if (m == "SBS_q90") S_C_q90 else NA_real_
      ) |>
      select(!se) |>
      mutate(metric = !!m) |>
      relocate(metric, .after = 1) |>
      rename(diff = !!mean_diff_col,
             sd   = !!mean_sd_col)
  }) |>
  bind_rows()

data$metric = factor(
  data$metric,
  levels = measures,
  labels = c("SBS (Early, q10)", "SBS (Median, q50)", "SBS (Late, q90)", "ISBS")
)

admin_all_stats =
  data |>
  group_by(n, metric) |>
  summarize(
    n_violations = sum(signif_violation),
    violation_rate = mean(signif_violation),
    mean_hat_epsilon = mean(epsilon),
    diff_mean = if (any(signif_violation)) mean(diff[signif_violation]) else NA_real_,
    diff_median = if (any(signif_violation)) median(diff[signif_violation]) else NA_real_,
    diff_min = if (any(signif_violation)) min(diff[signif_violation]) else NA_real_,
    diff_max = if (any(signif_violation)) max(diff[signif_violation]) else NA_real_,
    .groups = "drop"
  )

admin_all_stats |>
  arrange(metric) |>
  DT::datatable(
    rownames = FALSE,
    options = list(dom = "t", searching = FALSE),
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: center; font-size:142%",
      "SBS: Empirical violations of properness for large sample size + administrative censoring")
  ) |>
  formatRound(columns = 4:9, digits = 5)

Table C1

Note

The only difference with the previous experiment is that S_C(t) is now being estimated using the marginal Kaplan-Meier model via survival::survfit() instead of using the true Weibull censoring distribution (see helper.R). For SBS/ISBS we use constant interpolation of the censoring survival distribution S_C(t). For RCLL* we use linear interpolation of S_C(t) to mitigate density estimation issues, i.e. for f_C(t).

Use properness_test.R to run the experiments. To run for different sample sizes n, use run_tests.sh, changing estimate_cens to TRUE.

To get the results displayed, we executed the following command: ./run_tests.sh 10 25 50 100 250 500 1000 2500 5000 10000 (this will take a while to run, depending on your system).

Merged simulation results for all sample sizes are available here.

Load results (same output columns as in the previous section):

res = readRDS("results/res_sims10000_distrs1000_1.rds") |>
  select(!matches("shape|scale|prop_cens|tv_dist|n_events|max_t|m_true|m_pred")) # remove columns

As before, we compute 95% confidence intervals for the score differences across m = 1000 draws per simulation:

measures = c("SBS_q10", "SBS_median", "SBS_q90", "ISBS", "RCLL", "wRCLL")
n_distrs = 1000 # `m` value from paper's experiment
tcrit = qt(0.975, df = n_distrs - 1) # 95% t-test
threshold = 1e-4

data = measures |>
  lapply(function(m) {
    mean_diff_col = paste0(m, "_diff")
    mean_sd_col = paste0(m, "_sd")

    res |>
      select(sim, n, !!mean_diff_col, !!mean_sd_col, starts_with("shift_"),
             epsilon, starts_with("S_")) |>
      select(-ends_with("tail")) |>
      mutate(
        se = .data[[mean_sd_col]] / sqrt(n_distrs),
        CI_lower = .data[[mean_diff_col]] - tcrit * se,
        CI_upper = .data[[mean_diff_col]] + tcrit * se,
        signif_violation = .data[[mean_diff_col]] > threshold & CI_lower > 0,
        # Keep relevant SBS columns, just in case
        shift_q10 = if (m == "SBS_q10") shift_q10 else NA_real_,
        S_Y_q10   = if (m == "SBS_q10") S_Y_q10 else NA_real_,
        S_C_q10   = if (m == "SBS_q10") S_C_q10 else NA_real_,
        shift_med = if (m == "SBS_median") shift_med else NA_real_,
        S_Y_med   = if (m == "SBS_median") S_Y_med else NA_real_,
        S_C_med   = if (m == "SBS_median") S_C_med else NA_real_,
        shift_q90 = if (m == "SBS_q90") shift_q90 else NA_real_,
        S_Y_q90   = if (m == "SBS_q90") S_Y_q90 else NA_real_,
        S_C_q90   = if (m == "SBS_q90") S_C_q90 else NA_real_
      ) |>
      select(!se) |>
      mutate(metric = !!m) |>
      relocate(metric, .after = 1) |>
      rename(diff = !!mean_diff_col,
             sd   = !!mean_sd_col)
  }) |>
  bind_rows()

data$metric = factor(
  data$metric,
  levels = measures,
  labels = c("SBS (Q10)", "SBS (Median)", "SBS (Q90)", "ISBS", "RCLL", "RCLL*")
)

Lastly, we summarize the significant violations across sample sizes and scoring rules:

all_stats =
  data |>
  group_by(n, metric) |>
  summarize(
    n_violations = sum(signif_violation),
    violation_rate = mean(signif_violation),
    diff_mean = if (any(signif_violation)) mean(diff[signif_violation]) else NA_real_,
    diff_median = if (any(signif_violation)) median(diff[signif_violation]) else NA_real_,
    diff_min = if (any(signif_violation)) min(diff[signif_violation]) else NA_real_,
    diff_max = if (any(signif_violation)) max(diff[signif_violation]) else NA_real_,
    .groups = "drop"
  )

all_stats |>
  arrange(metric) |>
  DT::datatable(
    rownames = FALSE,
    options = list(searching = FALSE),
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: center; font-size:150%",
      "Table C1: Empirical violations of properness using G(t)")
  ) |>
  DT::formatRound(columns = 4:8, digits = 5)

TipSummary

Estimating S_C(t) via Kaplan-Meier produced nearly identical results to using the true censoring distribution, with no violations for RCLL/RCLL* and only small-sample effects for SBS and ISBS.

Sensitivity to Model Misspecification

To evaluate how effectively each scoring rule discriminates between increasingly misspecified models, we designed a simulation study based on a non-proportional hazards data-generating process (DGP). The DGP is a log-normal accelerated failure time (AFT) model with covariate-dependent location and group-specific scale parameters, producing crossing survival curves. Further details on the DGP, generated tasks, models, and benchmark design are provided in the manuscript.

Note

The simulation and benchmarking pipeline is implemented in the simulation_benchmark directory. Key scripts include:

• simulate.R: Defines the DGP and can generate data sets with different censoring rates.
• train.R: Generates two tasks, one with a low censoring rate and one with a high censoring rate, and saves them as tasks.rds. Then trains all models on the generated tasks. The trained models are saved as trained_low.rds, trained_high.rds, and separate files for the large Random Survival Forest objects (trained_low_rsf.rds, trained_high_rsf.rds).
• execute_dense.R: Executes a benchmark where the specified models are evaluated on test sets of varying sizes using a dense prediction time grid (unique event times from the training task). Results are saved in results/sim_bm_dense.rds.
• RCLL_investigation.R: Investigates the sensitivity of the RCLL to the resolution of the prediction time grid. Results are saved in results/rcll_sim.rds.

DGP and censoring scenarios

The code below generates Kaplan–Meier curves for both tasks with additional summary statistics, stratified by the binary treatment variable. See Figure 4 in the paper. Below, we also display the estimated censoring survival function S_C(t) via Kaplan–Meier alongside the true exponential censoring distribution for each DGP.

tasks = readRDS("simulation_benchmark/tasks.rds")

plot_km_task = function(task, title, show_censor = FALSE, extend_to_last_time = TRUE) {
  d = task$data()

  # summary quantities
  cens_total = mean(d$status == 0)
  cens_x20 = mean(d$status == 0 & d$x2 == 0)
  cens_x21 = mean(d$status == 0 & d$x2 == 1)
  n_events = sum(d$status == 1)
  admin_cens = mean(d$time == 10)

  subtitle =
    paste0(
      sprintf("Total censoring: %.1f%% (x2 = 0: %.1f%%, x2 = 1: %.1f%%)\n",
        100 * cens_total, 100 * cens_x20, 100 * cens_x21
      ),
      sprintf("Number of events: %d\n", n_events),
      sprintf("Administrative censoring: %.1f%%", 100*admin_cens)
    )

  # Kaplan-Meier fit
  fit = survival::survfit(
    survival::Surv(time, status) ~ factor(x2),
    data = d
  )

  s = summary(fit)
  df = data.frame(
    time = s$time,
    surv = s$surv,
    lower = s$lower,
    upper = s$upper,
    strata = s$strata,
    n_censor = s$n.censor
  )

  # ensure S(0) = 1 for each strata
  by_strata = split(df, df$strata)
  s0_rows = lapply(by_strata, function(dd) {
    if (any(dd$time == 0)) {
      NULL
    } else {
      data.frame(
        time = 0,
        surv = 1,
        lower = 1,
        upper = 1,
        strata = dd$strata[1],
        n_censor = 0
      )
    }
  })
  s0_rows = do.call(rbind, s0_rows)
  if (!is.null(s0_rows)) {
    df = rbind(s0_rows, df)
  }

  if (extend_to_last_time) {
    max_time = max(d$time)
    by_strata = split(df, df$strata)
    ext = lapply(by_strata, function(dd) {
      last_row = dd[nrow(dd), , drop = FALSE]
      if (max_time > last_row$time) {
        last_row$time = max_time
        last_row
      } else {
        NULL
      }
    })
    ext = do.call(rbind, ext)
    if (!is.null(ext)) {
      df = rbind(df, ext)
    }
  }

  censor_df = df[df$n_censor > 0, , drop = FALSE]

  p =
    ggplot(df, aes(x = time, y = surv, color = strata, linetype = strata)) +
    geom_ribbon(
      aes(ymin = lower, ymax = upper, group = strata),
      fill = "grey80",
      alpha = 0.4,
      color = NA
    ) +
    geom_step() +
    scale_color_manual(
      values = c("black", "red"),
      labels = c(
        expression(x[2] == 0 ~ paste("(", sigma, "  = ", 0.5, ")")),
        expression(x[2] == 1 ~ paste("(", sigma, "  = ", 1.5, ")"))
      )
    ) +
    scale_linetype_manual(
      values = c(1, 1),
      labels = c(
        expression(x[2] == 0 ~ paste("(", sigma, "  = ", 0.5, ")")),
        expression(x[2] == 1 ~ paste("(", sigma, "  = ", 1.5, ")"))
      )
    ) +
    labs(
      title = title,
      subtitle = subtitle,
      x = "Time",
      y = "Survival probability",
      color = NULL,
      linetype = NULL
    ) +
    theme_bw(base_size = 14) +
    theme(
      legend.position = c(0.98, 0.98),
      legend.justification = c(1, 1),
      legend.background = element_rect(fill = "white", color = "grey80"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 16),
      legend.key.size = grid::unit(1, "lines")
    )

  if (show_censor && nrow(censor_df) > 0) {
    p = p + geom_point(
      data = censor_df,
      aes(x = time, y = surv, color = strata),
      shape = 3,
      stroke = 0.9,
      size = 2,
      show.legend = FALSE
    )
  }

  p
}

plot_cens_surv = function(tasks, task_id) {
  surv_obj = tasks[[task_id]]$truth()
  time = surv_obj[, 1]
  status = surv_obj[, 2] # 1 for event, 0 for censored
  surv_cens = Surv(time, 1 - status) # now "event" means "censored"

  # Fit KM to the censoring survival
  km_cens_fit = survfit(surv_cens ~ 1)
  km_cens_df = data.frame(
    time = c(0, km_cens_fit$time),
    surv = c(1, km_cens_fit$surv)
  )

  # True lambda for the censoring distribution
  lambda = ifelse(task_id == "low", 0.075, 0.45)

  # Create the comparison plot
  ggplot() +
    # KM for censoring: black
    geom_step(
      data = km_cens_df,
      aes(x = time, y = surv, color = "Kaplan-Meier"),
      linewidth = 1,
      linetype = "solid"
    ) +
    # True exponential survival for censoring: blue
    stat_function(
      aes(color = "True"),
      fun = function(t) exp(-lambda * t),
      linewidth = 1.2
    ) +
    scale_color_manual(
      name = "",
      values = c("Kaplan-Meier" = "black", "True" = "blue")
    ) +
    labs(
      title = bquote("" * S[C](t) * " Estimation (True " * lambda[C] * " = " * .(lambda) * ")"),
      x = "Time",
      y = "Censoring Survival Probability"
    ) +
    theme_bw(base_size = 14)
}

Low censoring DGP

p_low  = plot_km_task(tasks$low, title = "Low censoring task", show_censor = TRUE)
p_low

plot_cens_surv(tasks, "low")

High censoring DGP

p_high = plot_km_task(tasks$high, title = "High censoring task", show_censor = TRUE)
p_high

plot_cens_surv(tasks, "high")

Large sample size sensitivity

We evaluate how effectively each scoring rule discriminates between correctly specified and increasingly misspecified models under a large test set. For each task, we generate 100 Monte Carlo test sets of size n = 1000. For each replicate, we compute the excess loss,

\Delta L = L(\hat S) - L(S_Y),

i.e. the difference in empirical mean loss between a candidate model and the true data-generating process. We then plot it against the mean integrated absolute error (MIAE) between the predicted and true survival functions, which serves as a measure of misspecification.

The panels below correspond to Figure 5 in the main manuscript and show the relationship between excess loss and MIAE for four commonly used evaluation metrics:

• RCLL: right-censored log-loss (density estimated via linear interpolation of the survival function)
• SBS: Survival Brier Score evaluated at t = 5
• ISBS: Integrated Survival Brier Score, integrated up to the time corresponding to the 80\% proportion of censoring in each test set
• Harrell’s C-index: a discrimination metric (not a proper scoring rule), included for comparison purposes

res = readRDS("results/sim_bm_dense.rds")

# Plot score differences vs. MIAE/MISE for multiple measures and one test size value
# Score diffs are (predicted - true)
plot_score_diff = function(
  res,
  y = "miae", # "mise" or "miae"
  task_id_val = "low", # "low" or "high"
  n_test_val = 1000, # test set size to display
  drop_learners = c("CoxPH"), # remove undesired learners
  rcll_truth = "interp", # "true" => (true S, true f), "interp" => (true S, interpolated f from linear S)
  meas = c("rcll", "cindex", "sbs", "isbs"), # measures to display
  legend_position = "right" # "right", "left", "top", "bottom", "none"
) {
  # Subset to the chosen task
  dt = res[task_id == task_id_val & n_test == n_test_val]
  # Remove unwanted learners
  dt = dt[!learner_id %in% drop_learners]

  # Rename learner_ids
  label_map = c(
    "LogNorm_int_shape_x2"   = "Oracle",
    "LogLog_int_shape_x2"    = "LLogis",
    "Weib_int_shape_x2"      = "Weibull",
    "LogNorm_noint_shape_x2" = "LogNorm_scale",
    "LogNorm_int_noshape"    = "LogNorm_int",
    "LogNorm_noint_noshape"  = "LogNorm",
    "CoxPH_int"              = "Cox_int",
    "RSF"                    = "RSF",
    "KM"                     = "KM"
  )

  dt[, learner_id := ifelse(
    learner_id %in% names(label_map),
    label_map[learner_id],
    learner_id
  )]

  # --- Extract true model references per replicate ---
  true_dt = dt[learner_id == "true",
    .(rsmp_id,
      rcll_true_full = true_rcll,
      rcll_true_interp = rcll,
      true_cindex = cindex,
      true_sbs = sbs,
      true_isbs = isbs)]

  dt = merge(dt, true_dt, by = "rsmp_id")

  # --- Compute differences ---
  # RCLL: compare to either full true Likelihood or f interpolated from linear S
  if (rcll_truth == "interp") {
    dt[, diff_rcll := rcll - rcll_true_interp]
  } else if (rcll_truth == "true") {
    # comparing with the true model's RCLL, makes the rest of the
    # metric comparisons a bit less comparable (as we don't have the true ISBS or SBS)
    dt[, diff_rcll := rcll - rcll_true_full]
  }

  # C‑index: error = 1 - cindex
  # error_pred - error_true = (1 - cindex_pred) - (1 - cindex_true) = cindex_true - cindex_pred
  dt[, diff_cindex := true_cindex - cindex]

  # Brier score
  dt[, diff_sbs := sbs - true_sbs]

  # Integrated Brier score
  dt[, diff_isbs := isbs - true_isbs]

  # Remove true model
  dt = dt[learner_id != "true"]

  # --- Reshape to long format ---
  measures = c("rcll", "cindex", "sbs", "isbs")
  measures = intersect(measures, meas)

  long_dt = melt(dt,
                 id.vars = c("learner_id", "mise", "miae", "rsmp_id"),
                 measure.vars = paste0("diff_", measures),
                 variable.name = "measure",
                 value.name = "diff")
  long_dt[, measure := sub("diff_", "", measure)]

  # Learner colours
  cols = c(
    Oracle        = "#dfce0cff",
    LogNorm_scale = "#2171B5",
    LogNorm_int   = "#6BAED6",
    LogNorm       = "#C6DBEF",
    LLogis        = "#FFA500",
    Weibull       = "#756BB1",
    Cox_int       = "#238B45",
    RSF           = "#D62728",
    KM            = "#7F7F7F"
  )

  # Keep only colours that appear in the data
  present_learners = unique(long_dt$learner_id)
  cols = cols[names(cols) %in% present_learners]

  # Label measures nicely
  measure_labels = c(
    rcll   = if (rcll_truth == "interp") "RCLL" else "RCLL (difference to true likelihood)",
    sbs    = "SBS (t = 5)",
    isbs   = "ISBS (τ* = q90)",
    cindex = "1 − C‑index (error)"
  )
  measure_labels = measure_labels[names(measure_labels) %in% measures]
  long_dt[, measure_label := factor(measure,
                                    levels = names(measure_labels),
                                    labels = measure_labels)]

  # Order learners for the legend
  learner_order = c(
    "Oracle",
    "LLogis",
    "Weibull",
    "LogNorm_scale",
    "LogNorm_int",
    "LogNorm",
    "Cox_int",
    "RSF",
    "KM"
  )
  learner_order = learner_order[learner_order %in% present_learners]
  long_dt[, learner_id := factor(learner_id, levels = learner_order)]

  # Create the plot
  ggplot(long_dt, aes(x = .data[[y]], y = diff, color = learner_id)) +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = cols) +
    facet_wrap(~ measure_label, scales = "free_y") +
    labs(
      title = sprintf(
        "%s censoring task (Test set size = %d)", # unique event times as grid by default
        paste0(toupper(substring(task_id_val, 1, 1)), substring(task_id_val, 2)),
        n_test_val
      ),
      x = ifelse(y == "mise",
                 "MISE (distance to true S)",
                 "MIAE (distance to true S)"),
      #y = "Score difference (predicted − true model)",
      y = expression(L[pred] - ~L[true]),
      color = "Model"
    ) +
    theme_bw(base_size = 12, base_family = "Arial") +
    theme(
      strip.background = element_rect(fill = "grey95"),
      strip.text = element_text(face = "bold"),
      legend.position = if (legend_position == "none") "none" else legend_position,
      legend.title = element_text(size = 14),
      legend.text = element_text(size = 12)
    ) +
    guides(color = guide_legend(override.aes = list(size = 5)))
}

Low censoring DGP

plot_score_diff(res, task_id_val = "low", meas = c("rcll", "sbs", "isbs", "cindex"),
                n_test_val = 1000, legend_position = "bottom")

High censoring DGP

plot_score_diff(res, task_id_val = "high", meas = c("rcll", "sbs", "isbs", "cindex"),
                n_test_val = 1000, legend_position = "bottom")

TipSummary

• RCLL shows the strongest discrimination. Excess loss increases almost monotonically with MIAE, producing well-separated model clusters with low variability. All misspecified models have positive excess loss, consistent with empirical properness.
• SBS (at t = 5) provides poor discrimination. Models exhibit substantial overlap, weak correspondence with MIAE, and high variability, particularly under heavy censoring, making SBS unsuitable for reliably ranking models.
• ISBS largely preserves the global ordering induced by MIAE but compresses models into a few broad clusters with similar excess loss, limiting discrimination among similarly performing models. Occasional negative \Delta L values arise from finite-sample variability, especially under high censoring.
• Harrell’s C-index is not aligned with model misspecification: under high censoring it can assign better scores to misspecified models than to the true data-generating process, illustrating its lack of properness.

Effect of Test Set Size

To investigate the finite-sample behaviour of the scoring rules, we repeat the simulation benchmark while varying the test set size. As before, performance is assessed using the excess loss relative to the true data-generating process, but now across increasingly smaller test sets. This experiment illustrates how sample size affects both the discriminatory power of the scoring rules and the frequency of empirical properness violations.

The figures below correspond to Figure 6 (low censoring) and Figure C1 (high censoring) in the manuscript.

# Plot score differences vs. MIAE/MISE for one measure and multiple test set sizes
# Score diffs are (predicted - true)
plot_score_diff_by_ntest = function(
  res,
  measure_id = "rcll",# "rcll", "cindex", "sbs", "isbs"
  y = "miae", # "mise" or "miae"
  task_id_val = "low", # "low" or "high"
  n_test_vals = c(10, 25, 50, 100, 250, 500), # test sizes to facet
  drop_learners = c("CoxPH"),
  rcll_truth = "interp", # "true" or "interp"
  legend_position = "right" # "right", "left", "top", "bottom", "none"
) {
  # subset to the chosen task and test sizes
  dt = res[task_id == task_id_val & n_test %in% n_test_vals]
  # Remove unwanted learners
  dt = dt[!learner_id %in% drop_learners]
  # Rename learner_ids
  label_map = c(
    "LogNorm_int_shape_x2"   = "Oracle",
    "LogLog_int_shape_x2"    = "LLogis",
    "Weib_int_shape_x2"      = "Weibull",
    "LogNorm_noint_shape_x2" = "LogNorm_scale",
    "LogNorm_int_noshape"    = "LogNorm_int",
    "LogNorm_noint_noshape"  = "LogNorm",
    "CoxPH_int"              = "Cox_int",
    "RSF"                    = "RSF",
    "KM"                     = "KM"
  )

  dt[, learner_id := ifelse(
    learner_id %in% names(label_map),
    label_map[learner_id],
    learner_id
  )]

  # True references
  true_dt = dt[learner_id == "true",
    .(rsmp_id,
      n_test,
      rcll_true_full = true_rcll,
      rcll_true_interp = rcll,
      true_cindex = cindex,
      true_sbs = sbs,
      true_isbs = isbs)
  ]

  dt = merge(dt, true_dt, by = c("rsmp_id", "n_test"))

  # Differences
  if (measure_id == "rcll") {
    if (rcll_truth == "interp") {
      dt[, diff := rcll - rcll_true_interp]
    } else {
      dt[, diff := rcll - rcll_true_full]
    }
  } else if (measure_id == "cindex") {
    dt[, diff := true_cindex - cindex]
  } else if (measure_id == "sbs") {
    dt[, diff := sbs - true_sbs]
  } else if (measure_id == "isbs") {
    dt[, diff := isbs - true_isbs]
  }

  # Remove true model
  dt = dt[learner_id != "true"]

  # Ordering
  dt[, n_test := factor(n_test, levels = sort(n_test_vals))]

  learner_order = c(
    "Oracle",
    "LLogis",
    "Weibull",
    "LogNorm_scale",
    "LogNorm_int",
    "LogNorm",
    "Cox_int",
    "RSF",
    "KM"
  )
  present_learners = unique(dt$learner_id)
  learner_order = learner_order[learner_order %in% present_learners]
  dt[, learner_id := factor(learner_id, levels = learner_order)]

  cols = c(
    Oracle        = "#dfce0cff",
    LLogis        = "#FFA500",
    Weibull       = "#756BB1",
    LogNorm_scale = "#2171B5",
    LogNorm_int   = "#6BAED6",
    LogNorm       = "#C6DBEF",
    Cox_int       = "#238B45",
    RSF           = "#D62728",
    KM            = "#7F7F7F"
  )
  cols = cols[names(cols) %in% present_learners]

  measure_labels = c(
    rcll   = if (rcll_truth == "interp") "RCLL" else "RCLL (difference to true likelihood)",
    sbs    = "SBS (t = 5)",
    isbs   = "ISBS (τ* = q90)",
    cindex = "1 − C‑index (error)"
  )

  ggplot(dt, aes(x = .data[[y]], y = diff, color = learner_id)) +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = cols) +
    facet_grid(
      . ~ n_test,
      labeller = labeller(
        n_test = function(x) paste0("n['test']==", x),
        .default = label_parsed
      )
    ) +
    labs(
      title = sprintf(
        "%s censoring task — %s",
        paste0(toupper(substring(task_id_val, 1, 1)), substring(task_id_val, 2)),
        measure_labels[[measure_id]]
      ),
      x = ifelse(y == "mise", "MISE (distance to true S)", "MIAE (distance to true S)"),
      y = expression(L[pred] - ~L[true]),
      color = "Model"
    ) +
    theme_bw(base_size = 14, base_family = "Arial") +
    theme(
      strip.background = element_rect(fill = "grey95"),
      strip.text = element_text(size = 16, face = "bold"),
      legend.position = if (legend_position == "none") "none" else legend_position,
      legend.title = element_text(size = 14),
      legend.text = element_text(size = 12)
    ) +
    guides(color = guide_legend(override.aes = list(size = 5)))
}

Low censoring DGP

p1 = plot_score_diff_by_ntest(res, measure_id = "rcll", task_id_val = "low", legend_position = "none")
p2 = plot_score_diff_by_ntest(res, measure_id = "isbs", task_id_val = "low", legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5, alpha = 1), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

High censoring DGP

p1 = plot_score_diff_by_ntest(res, measure_id = "rcll", task_id_val = "high", legend_position = "none")
p2 = plot_score_diff_by_ntest(res, measure_id = "isbs", task_id_val = "high", legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5, alpha = 1), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

TipSummary

• Discriminatory power decreases as the test set becomes smaller. For both RCLL and ISBS, excess loss becomes increasingly variable, leading to greater overlap between models and weaker correspondence with model misspecification.

• RCLL generally retains better discrimination at small sample sizes. Although model separation deteriorates as n_{\mathrm{test}} decreases, RCLL can still distinguish strongly misspecified models more effectively than ISBS.

• Empirical properness is increasingly affected by finite-sample variability. Negative excess losses become more frequent as the test set size decreases, particularly for ISBS. For test sets of at least 100 observations, such violations are uncommon for both scoring rules.

• Heavy censoring further amplifies finite-sample effects. Under high censoring, variability increases and model discrimination deteriorates more rapidly, especially when evaluating with the ISBS.

Impact of Prediction Grid Resolution on RCLL

Unlike the Brier score, the RCLL requires both the survival function and event density. For machine learning models that predict survival probabilities only at a discrete set of time points, the density must therefore be approximated by interpolating the survival curve. To assess the impact of this approximation, we progressively reduce the prediction time grid by uniformly subsampling the unique event times, from 100% down to 2% of the original grid.

The figures below show how the resulting approximation (grid size B) affects the excess RCLL for different test set sizes. The results for n = 250 correspond to Figure 7, Table 3 and Table C2 in the manuscript.

# plots score differences (relative to TRUE likelihood) vs. MIAE/MISE for RCLL
# on multiple prediction grid resolutions (proportions of unique event times)
plot_rcll_diff_by_grid = function(
  res,
  y = "miae", # "mise" or "miae"
  task_id_val = "low", # "low" or "high"
  n_test_val = 500, # test set size to display
  prop_vals = c(0.02, 0.05, 0.1, 0.2, 0.5, 0.8, 1), # unique event time proportions to display
  drop_learners = c("CoxPH"),
  legend_position = "right" # "right", "left", "top", "bottom", "none"
) {
  # subset to the chosen task and test size
  dt = res[task_id == task_id_val & prop %in% prop_vals & n_test == n_test_val]
  # Remove unwanted learners
  dt = dt[!learner_id %in% drop_learners]
  # Rename learner_ids
  label_map = c(
    "LogNorm_int_shape_x2"   = "Oracle",
    "LogLog_int_shape_x2"    = "LLogis",
    "Weib_int_shape_x2"      = "Weibull",
    "LogNorm_noint_shape_x2" = "LogNorm_scale",
    "LogNorm_int_noshape"    = "LogNorm_int",
    "LogNorm_noint_noshape"  = "LogNorm",
    "CoxPH_int"              = "Cox_int",
    "RSF"                    = "RSF",
    "KM"                     = "KM"
  )

  dt[, learner_id := ifelse(
    learner_id %in% names(label_map),
    label_map[learner_id],
    learner_id
  )]

  # Extract true references (per replicate + prop)
  true_dt = dt[learner_id == "true",
    .(rsmp_id,
      prop,
      rcll_true_full = true_rcll, # true S + true f
      rcll_true_interp = rcll # true S + interpolated f
    )
  ]

  dt = merge(dt, true_dt, by = c("rsmp_id", "prop"))

  # Compute differences relative to TRUE likelihood
  dt[, diff := rcll - rcll_true_full]

  # Rename true model with interpolated density
  dt[learner_id == "true", learner_id := "Truth (true S, interp f)"]

  # Factor ordering
  # Create facet label: B = n_times (percentage)
  dt[, facet_label := sprintf("B == %d~'('*%d*'%%'*')'", n_times, round(prop * 100))]
  dt[, facet_label := factor(facet_label, levels = dt[order(prop), unique(facet_label)])]

  learner_order = c(
    "Truth (true S, interp f)",
    "Oracle",
    "LLogis",
    "Weibull",
    "LogNorm_scale",
    "LogNorm_int",
    "LogNorm",
    "Cox_int",
    "RSF",
    "KM"
  )
  present_learners = unique(dt$learner_id)
  learner_order = learner_order[learner_order %in% present_learners]
  dt[, learner_id := factor(learner_id, levels = learner_order)]

  # Colors
  cols = c(
    `Truth (true S, interp f)` = "#000000",
    Oracle        = "#dfce0cff",
    LLogis        = "#FFA500",
    Weibull       = "#756BB1",
    LogNorm_scale = "#2171B5",
    LogNorm_int   = "#6BAED6",
    LogNorm       = "#C6DBEF",
    Cox_int       = "#238B45",
    RSF           = "#D62728",
    KM            = "#7F7F7F"
  )
  cols = cols[names(cols) %in% present_learners]

  # Plot
  ggplot(dt, aes(x = .data[[y]], y = diff, color = learner_id)) +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = cols) +
    facet_grid(
      . ~ facet_label,
      labeller = label_parsed
    ) +
    labs(
      title = sprintf(
        "%s censoring task (Test set size = %d)",
        paste0(toupper(substring(task_id_val, 1, 1)), substring(task_id_val, 2)),
        n_test_val
      ),
      x = ifelse(y == "mise", "MISE (distance to true S)", "MIAE (distance to true S)"),
      y = expression(L[pred] - L[true]),
      color = "Model"
    ) +
    theme_bw(base_size = 14, base_family = "Arial") +
    theme(
      strip.background = element_rect(fill = "grey95"),
      strip.text = element_text(size = 16, face = "bold"),
      legend.position = if (legend_position == "none") "none" else legend_position,
      legend.title = element_text(size = 14),
      legend.text = element_text(size = 12)
    ) +
    guides(color = guide_legend(override.aes = list(size = 5, alpha = 1)))
}

# load the results of the RCLL grid resolution investigation
res = readRDS("results/rcll_sim.rds")

# duplicate the truth model
truth_exact = res[
  learner_id == "true",
  .(
    learner_id = "Truth (exact)",
    rcll = true_rcll,
    task_id,
    rsmp_id,
    prop,
    n_times,
    n_test
  )
]

truth_interp = res[
  learner_id == "true",
  .(
    learner_id = "Truth (interp.)",
    rcll,
    task_id,
    rsmp_id,
    prop,
    n_times,
    n_test
  )
]

others = res[
  !learner_id %in% c("true", "CoxPH"),
  .(
    learner_id,
    rcll,
    task_id,
    rsmp_id,
    prop,
    n_times,
    n_test
  )
]

res_tbl = rbindlist(list(truth_interp, truth_exact, others))

# rename learners (and order according to misspecification)
label_map = c(
  "Truth (exact)"    = "Truth (true S, true f)",
  "Truth (interp.)"  = "Truth (true S, interp f)",
  "LogNorm_int_shape_x2"   = "Oracle",
  "LogLog_int_shape_x2"    = "LLogis",
  "Weib_int_shape_x2"      = "Weibull",
  "LogNorm_noint_shape_x2" = "LogNorm_scale",
  "LogNorm_int_noshape"    = "LogNorm_int",
  "LogNorm_noint_noshape"  = "LogNorm",
  "CoxPH_int"              = "Cox_int",
  "RSF"                    = "RSF",
  "KM"                     = "KM"
)

res_tbl[, learner_id := factor(
  label_map[learner_id],
  levels = c(
    "Truth (true S, true f)",
    "Truth (true S, interp f)",
    "Oracle",
    "LLogis",
    "Weibull",
    "LogNorm_scale",
    "LogNorm_int",
    "LogNorm",
    "Cox_int",
    "KM",
    "RSF"
  )
)]

# easy way to print summarized absolute RCLL values
plot_rcll_table = function(
  res_tbl,
  task = c("low", "high"),
  n = 250,
  caption = NULL,
  digits = 3
) {
  tbl = res_tbl[
    task_id == task & n_test == n,
    .(
      mean = mean(rcll),
      sd = sd(rcll)
    ),
    by = .(learner_id, n_times)
  ]

  fmt = paste0("%.", digits, "f \u00B1 %.", digits, "f")
  tbl[, value := sprintf(fmt, mean, sd)]

  # to long format
  tbl = dcast(
    tbl,
    learner_id ~ n_times,
    value.var = "value"
  )

  DT::datatable(
    tbl,
    rownames = FALSE,
    options = list(
      pageLength = nrow(tbl),
      dom = "t"
    ),
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: center; font-size:150%",
      sprintf("RCLL (mean \u00B1 SD) for %s censoring task, n = %d.", task, n)
    )
  )
}

n = 50

p1 = plot_rcll_diff_by_grid(res, task_id_val = "low", n_test_val = 50, legend_position = "none")
p2 = plot_rcll_diff_by_grid(res, task_id_val = "high", n_test_val = 50, legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

plot_rcll_table(res_tbl, task = "low", n = 50)

plot_rcll_table(res_tbl, task = "high", n = 50)

n = 100

p1 = plot_rcll_diff_by_grid(res, task_id_val = "low", n_test_val = 100, legend_position = "none")
p2 = plot_rcll_diff_by_grid(res, task_id_val = "high", n_test_val = 100, legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

plot_rcll_table(res_tbl, task = "low", n = 100)

plot_rcll_table(res_tbl, task = "high", n = 100)

n = 250

p1 = plot_rcll_diff_by_grid(res, task_id_val = "low", n_test_val = 250, legend_position = "none")
p2 = plot_rcll_diff_by_grid(res, task_id_val = "high", n_test_val = 250, legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

plot_rcll_table(res_tbl, task = "low", n = 250)

plot_rcll_table(res_tbl, task = "high", n = 250)

n = 500

p1 = plot_rcll_diff_by_grid(res, task_id_val = "low", n_test_val = 500, legend_position = "none")
p2 = plot_rcll_diff_by_grid(res, task_id_val = "high", n_test_val = 500, legend_position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5), nrow = 1, byrow = TRUE)) +
  theme(legend.direction = "horizontal")
(p1 / p2)

plot_rcll_table(res_tbl, task = "low", n = 500)

plot_rcll_table(res_tbl, task = "high", n = 500)

TipSummary

• Prediction grid resolution has little influence on model ranking. Across all grid sizes and test set sizes, the ordering of the fitted models remains largely unchanged in both censoring scenarios, indicating that interpolation mainly affects the absolute scale of the RCLL rather than its ability to rank models.

• The absolute scale of the RCLL depends on the density of the prediction grid. The general trend is that for the true model, coarse prediction grids inflate the interpolated RCLL above the exact likelihood, whereas all fitted models exhibit systematically smaller RCLL values, making them appear artificially closer to the truth. While moderate grid coarsening has little effect on model ranking, sufficiently dense prediction grids are needed to minimize interpolation bias.

• The absolute scale of the RCLL depends strongly on the censoring rate. All models attain substantially smaller RCLL values under high censoring compared to under low censoring.

Benchmark (Real-World Data)

Note

This benchmark complements the simulation study by comparing how different evaluation metrics assess and rank survival learners across eight heterogeneous, low-dimensional, real-world datasets.

To reproduce this benchmark, see real_data_benchmark.R.

Results are available in results/real_data_bm.rds.

Datasets

All the datasets are available from mlr3proba as mlr3 tasks.

library(mlr3proba)
library(mlr3misc)
# silence mlr3 logging messages
lgr::get_logger("mlr3")$set_threshold(0)

res = readRDS("results/real_data_bm.rds")

task_ids = unique(res$task_id)
tasks = mlr3::tsks(task_ids)

metrics = c("C-index", "D-calib", "RCLL*", "RCLL", "ISBS")
res_long = melt(
  res,
  id.vars = c("task_id", "learner_id", "iteration"),
  measure.vars = metrics,
  variable.name = "metric",
  value.name = "value"
)

# remove RCLL* from the boxplots
res_long = res_long[metric != "RCLL*"]

# truncate D-calibration to 50 to avoid extreme outliers in the boxplots
res_long[
  metric == "D-calib",
  value := pmin(value, 50)
]

# Reorder metric factor to the desired sequence
res_long[, metric := factor(metric, levels = metrics)]

task_info = lapply(tasks, function(task) {
  id = task$id
  n = task$nrow
  n_events = length(task$unique_event_times())
  p = task$n_features
  feature_dt = task$col_info[id %nin% c("..row_id", task$target_names)]
  n_factor = sum(feature_dt$type == "factor")
  n_numeric = sum(feature_dt$type %in% c("numeric", "integer"))
  cens_rate = task$cens_prop()
  admin_cens_rate = task$admin_cens_prop()
  prop_haz = ifelse(task$prop_haz() < 0.05, FALSE, TRUE)

  list(
    id = id,
    n = n,
    n_events = n_events,
    p = p,
    n_factor = n_factor,
    n_numeric = n_numeric,
    censoring_rate = cens_rate,
    admin_censoring_rate = admin_cens_rate,
    prop_haz = prop_haz
  )
}) |> rbindlist()

In the following table:

• id is the dataset identifier
• n is the number of observations
• n_events is the number of unique event times
• p is the number of features
• n_factor is the number of factor features
• n_numeric is the number of numeric features
• censoring_rate is the proportion of censored observations in the dataset
• admin_censoring_rate is the proportion of observations censored due to administrative censoring
• prop_haz is a boolean derived from the p-value of the global proportional hazards test computed with survival::cox.zph() at the 5% significance level: TRUE indicates no evidence against proportional hazards (p \ge 0.05), whereas FALSE indicates a violation (p < 0.05).

task_info |>
  arrange(censoring_rate) |>
  DT::datatable(
    rownames = FALSE,
    filter = "top",
    options = list(dom = "t", searching = TRUE)
  ) |>
  DT::formatRound(columns = c(7, 8), digits = 2) |>
  DT::formatStyle(
    columns = "censoring_rate",
    backgroundColor = DT::styleInterval(
      c(0.25, 0.5, 0.75),
      c("#e5f5e0", "#a1d99b", "#e7785f", "#f8555b")
    )
  )

Models

We used the following models:

• Kaplan-Meier (KM; non-parametric)
• Cox proportional hazards model (Cox; semi-parametric)
• Random survival forests (RSF; non-parametric)
• Log-normal accelerated failure time model (LogNorm; parametric)
• Weibull accelerated failure time model (Weibull; parametric)

Model performance is evaluated using three survival scoring rules (RCLL, RCLL*, and ISBS), Harrell’s C-index (discrimination), and D-calibration (calibration). All models are assessed using 10 repetitions of 5-fold cross-validation without hyperparameter tuning.

Performance per-dataset

Predictive performance can vary considerably across datasets due to differences in sample size, censoring proportion, covariate distributions and the suitability of the underlying model assumptions. To provide a detailed overview, the following figures show the distribution of each evaluation metric across all repeated cross-validation resamples, separately for each dataset. D-calibration is truncated at 50 to avoid extreme outliers in the boxplots. RCLL* is excluded from the boxplots as it is an experimental measure in this analysis and has not been used in practice. Datasets are ordered by increasing censoring rate, and the title of each panel includes the dataset identifier, sample size, number of features, and whether the proportional hazards assumption is satisfied. Within each dataset, learners are ordered by their median RCLL to facilitate comparison across the different evaluation metrics.

learner_ids = sort(unique(res_long$learner_id))
cols = scales::hue_pal()(length(learner_ids))
names(cols) = learner_ids

plot_task = function(task) {
  dat = res_long[task_id == task]

  # Order learners by median RCLL
  rcll_medians = dat[metric == "RCLL", .(med = median(value)), by = learner_id]
  learner_order = rcll_medians[order(med), learner_id]
  dat[, learner_id := factor(learner_id, levels = learner_order)]

  p = dat |>
    ggplot(aes(
      x = value,
      y = learner_id,
      colour = learner_id,
      fill = learner_id
    )) +
    geom_boxplot(alpha = 0.25, outlier.shape = NA) +
    geom_jitter(width = 0, height = 0.15, alpha = 0.25, size = 0.7) +
    facet_wrap(~metric, scales = "free_x", ncol = 2) +
    # ---- Add the manual scales ----
    scale_colour_manual(values = cols) +
    scale_fill_manual(values = cols) +
    labs(
      x = NULL,
      y = NULL,
      title = task
    ) +
    theme_bw(base_size = 14, base_family = "Arial") +
    theme(
      legend.position = "none",
      # strip.background = element_blank(),
      strip.background = element_rect(fill = "grey95"),
      strip.text = element_text(face = "bold"),
      panel.grid.major.y = element_blank(),
      plot.title = element_text(face = "bold")
    )

  p
}

# order datasets by censoring rate
task_ids = task_info[order(censoring_rate)]$id

cat("::: {.panel-tabset}\n\n")

mgus (n=176, p=7, PH)

veteran (n=137, p=6, PH violated)

lung (n=168, p=7, PH)

whas (n=481, p=9, PH violated)

gbsg (n=686, p=8, PH violated)

grace (n=1000, p=6, PH violated)

gbcs (n=686, p=8, PH violated)

rats (n=300, p=3, PH)

cat(":::\n")

For easier numerical comparison, the table below reports the mean performance of every learner on each dataset, averaged over the repeated cross-validation resamples.

res_mean = res[,
  lapply(.SD, mean, na.rm = TRUE), by = .(task_id, learner_id), .SDcols = metrics
]

res_mean |>
  DT::datatable(
    rownames = FALSE,
    filter = "top",
    options = list(searching = TRUE, pageLength = 10)
  ) |>
  DT::formatRound(columns = c("RCLL", "RCLL*", "C-index", "ISBS"), digits = 3) |>
  DT::formatSignif(columns = "D-calib", digits = 3)

Rank-correlation between evaluation metrics

Since the evaluation metrics operate on different numerical scales and some (such as RCLL) are affected by dataset characteristics like censoring, we compare learners using within-dataset rankings rather than absolute scores. For every dataset and resampling iteration, model performances are converted into ranks (1 = best). Harrell’s C-index is converted to an error metric (1 - C-index) so that all metrics are oriented in the same direction (smaller values indicate better performance).

We compute Spearman rank correlations and visualize the empirical agreement between evaluation metrics as a correlation heatmap.

library(ggcorrplot)

# convert C-index to 1 - C
res = res[, `C-index-error` := 1 - `C-index`]

metrics = c("C-index-error", "D-calib", "RCLL*", "RCLL", "ISBS")

# Compute ranks per task for each metric
rank_cols = paste0("rank_", metrics)
res_ranks = copy(res)

# calculate ranks for each metric, grouped by task_id and iteration
for (i in seq_along(metrics)) {
  m = metrics[i]
  res_ranks[,
    paste0("rank_", m) := frank(get(m), ties.method = "min"),
    by = .(task_id, iteration)
  ]
}

cor_all = cor(res_ranks[, ..rank_cols], method = "pearson")

ggcorrplot(
  cor_all,
  method = "square",
  type = "lower",
  hc.order = TRUE,
  lab = TRUE,
  lab_size = 4,
  outline.color = "white"
)

TipSummary

The real-world benchmark confirms that the choice of evaluation metric can substantially influence the ranking of survival learners.

• RCLL provides the greatest discrimination between competing models. Across nearly all datasets, RCLL produces the clearest separation between learners. With ISBS, models are bit less distinguishable, and the C-index often fails to differentiate between models with noticeably different predictive distributions. ISBS shows moderate agreement with both RCLL (\rho = 0.43) and the C-index (\rho = 0.64).

• RCLL values depends strongly on the censoring rate. As censoring increases, absolute RCLL values become smaller (and vice versa), making comparisons across datasets less meaningful. This behaviour follows directly from the RCLL definition: L_{RCLL} = - (\delta \log f(t) + (1 - \delta) \log S(t)). As the censoring rate increases, fewer observations contribute through the density term and more contribute through the survival term. For all learners and datasets considered here, the density term is consistently larger than the survival term on average, resulting in smaller absolute RCLL values with increasing censoring. Absolute RCLL scores are therefore most informative for comparing models within the same dataset, although learner rankings remain meaningful across datasets.

• D-calibration captures complementary information. Its near-zero correlation with the remaining metrics indicates that calibration reflects a different aspect of predictive performance than discrimination or probabilistic accuracy.

• Closely related scoring rules need not rank models similarly. Despite being a variant of RCLL, RCLL* exhibits only weak agreement with RCLL (\rho = 0.16), highlighting that seemingly similar evaluation metrics can induce substantially different learner rankings in practice.

Overall, these findings complement the simulation study by demonstrating that the empirical behaviour of survival evaluation metrics varies considerably on real datasets. Even metrics that are theoretically related—or theoretically proper—may emphasize different aspects of predictive performance and therefore lead to different conclusions about which model performs best.

R session info

utils::sessionInfo()

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_DK.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_DK.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_DK.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Oslo
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggcorrplot_0.1.4.1 DT_0.33            data.table_1.18.4  tidyr_1.3.1       
 [5] purrr_1.0.2        akima_0.6-3.6      patchwork_1.3.2    ggplot2_4.0.2     
 [9] tibble_3.3.0       dplyr_1.1.4        survival_3.7-0    

loaded via a namespace (and not attached):
 [1] gtable_0.3.6         xfun_0.57            bslib_0.10.0        
 [4] htmlwidgets_1.6.4    lattice_0.22-6       vctrs_0.7.1         
 [7] tools_4.4.2          crosstalk_1.2.1      generics_0.1.3      
[10] parallel_4.4.2       pkgconfig_2.0.3      Matrix_1.7-1        
[13] checkmate_2.3.4      RColorBrewer_1.1-3   mlr3pipelines_0.10.0
[16] S7_0.2.1             uuid_1.2-2           lifecycle_1.0.5     
[19] stringr_1.6.0        set6_0.2.6           compiler_4.4.2      
[22] farver_2.1.2         codetools_0.2-20     htmltools_0.5.9     
[25] sass_0.4.10          yaml_2.3.12          pillar_1.11.0       
[28] crayon_1.5.3         jquerylib_0.1.4      cachem_1.1.0        
[31] parallelly_1.47.0    tidyselect_1.2.1     digest_0.6.39       
[34] stringi_1.8.7        future_1.70.0        reshape2_1.4.4      
[37] listenv_0.10.1       labeling_0.4.3       splines_4.4.2       
[40] fastmap_1.2.0        grid_4.4.2           cli_3.6.6           
[43] magrittr_2.0.4       paradox_1.0.1        withr_3.0.3         
[46] scales_1.4.0         backports_1.5.1      sp_2.1-4            
[49] ooplah_0.2.0         rmarkdown_2.31       globals_0.19.1      
[52] distr6_1.8.4         evaluate_1.0.5       knitr_1.51          
[55] dictionar6_0.1.3     mlr3misc_0.22.0      rlang_1.2.0         
[58] Rcpp_1.1.1           glue_1.8.0           param6_0.2.4        
[61] palmerpenguins_0.1.1 mlr3proba_0.8.9      jsonlite_2.0.0      
[64] plyr_1.8.9           lgr_0.5.2            R6_2.6.1            
[67] mlr3_1.7.1