08. Interpreting the Black Box with SHAP & survex

Introduction

Flexible machine learning models can be difficult to interpret, so users may want summaries of how features contribute to a specific prediction.

SuperSurv provides an optional SHAP (SHapley Additive exPlanations) workflow through kernelshap. The explanation is computed for the fitted ensemble prediction function and should be interpreted for the selected prediction target and background sample.

This tutorial covers Global feature importance, Local (patient-level) explanations, and Time-Dependent survival analysis using the survex package.

1. Setup and Model Fitting

Let’s train a diverse Super Learner on our metabric dataset.

For a quick demonstration, we use 200 sampled patients and five predictors, with a 50-tree forest when available. We explain four patients relative to a background sample of ten. These small settings illustrate the workflow rather than provide a full-cohort interpretation.

library(SuperSurv)
library(survival)

data("metabric", package = "SuperSurv")
set.seed(42)
metabric <- metabric[sample(seq_len(nrow(metabric)), 200), ]

train_idx <- sample(1:nrow(metabric), 0.7 * nrow(metabric))
train <- metabric[train_idx, ]
test  <- metabric[-train_idx, ]

x_cols <- grep("^x", names(metabric))[1:5]
X_tr <- train[, x_cols, drop = FALSE]
X_te <- test[, x_cols, drop = FALSE]
new.times <- seq(50, 200, by = 25)
X_explain_subset <- X_te[1:4, , drop = FALSE]
X_background_subset <- X_tr[1:10, , drop = FALSE]

my_library <- c("surv.coxph", "surv.weibull")
if (has_rfsrc) {
  rf_library <- create_grid("surv.rfsrc", list(ntree = 50))
  my_library <- c(my_library, rf_library)
}

fit_sl <- SuperSurv(
  time = train$duration,
  event = train$event,
  X = X_tr,
  newdata = X_te,
  new.times = new.times,
  event.library = my_library,
  cens.library = c("surv.coxph"),
  control = list(
    saveFitLibrary = TRUE,
    event.t.grid = seq(0, max(train$duration[train$event == 1]), length.out = 30),
    cens.t.grid = seq(0, max(train$duration[train$event == 0]), length.out = 30)
  ),
  verbose = FALSE,
  selection = "ensemble",
  nFolds = 3
)

2. Global Explanations (Kernel SHAP)

We explain event probability by 100 months, 1 - S(100 | X). An explicit eval_time keeps every learner on the same probability scale. X_explain contains the patients to explain; X_background defines the reference population. All positive ensemble weights are retained, and the stored survival prediction methods apply the same screening and calibration used for ordinary prediction.

# Explain the ensemble event probability at one common horizon.
shap_vals <- explain_kernel(
  model = fit_sl, 
  X_explain = X_explain_subset, 
  X_background = X_background_subset, 
  nsim = 20,
  eval_time = 100
)

plot_global_importance(shap_vals, top_n = 5)

plot_beeswarm(shap_vals, data = X_explain_subset, top_n = 5)

3. Local (Patient-Level) Explanations

In precision medicine, we often need to explain why a specific patient has a high risk score. The Waterfall plot breaks down the exact algorithmic logic for an individual.

# Explain Patient #1 from our test subset
plot_patient_waterfall(shap_vals, patient_index = 1, top_n = 5)

4. Feature Dependence Plots

Does a specific biomarker have a linear or non-linear relationship with mortality risk?

plot_dependence(shap_vals, data = X_explain_subset, feature_name = "x0")

5. Patient Risk Trajectories (Survival Heatmap)

We can visualize patient trajectories over time using a survival heatmap generated directly from SuperSurv predictions.

# Plot the survival trajectories for the four patients being explained.
plot_survival_heatmap(fit_sl, newdata = X_explain_subset, times = new.times)

Interpretation: Patients at the top experience rapid drops in survival probability (high risk), while patients at the bottom maintain high survival probabilities.

6. Time-Dependent Explanations (survex Integration)

Traditional SHAP looks at an overall “risk score,” but survival analysis is fundamentally about time. A feature might be highly predictive of early mortality (Time = 50) but irrelevant for long-term survival (Time = 200).

SuperSurv natively bridges to the survex package to evaluate dynamic feature importance across the entire survival curve \(S(t)\).

library(survex)

# 1. Create the true survival object for the same explanation subset
y_explain <- survival::Surv(
  test$duration[seq_len(nrow(X_explain_subset))],
  test$event[seq_len(nrow(X_explain_subset))]
)

# 2. Build the survex explainer using our custom function
surv_explainer <- explain_survex(
  model = fit_sl, 
  data = X_explain_subset, 
  y = y_explain, 
  times = new.times
)

Once the explainer is created, you have full access to the survex ecosystem. Let’s look at three powerful time-dependent visualizations.

The following specialist analyses are shown as code but are not executed during the package vignette build because their runtime and accepted prediction representations depend on the installed survex release.

# Calculate time-dependent model parts (permutation feature importance)
time_importance <- model_parts(surv_explainer)

# Plot the dynamic importance over time
plot(time_importance)

# Calculate the partial dependence profile for feature 'x0'
pdp_time <- model_profile(surv_explainer, variables = "x0")

plot(pdp_time)

# Explain Patient #1 over time
patient_1_data <- X_explain_subset[1, , drop = FALSE]

# Calculate SurvSHAP(t)
survshap_t <- predict_parts(surv_explainer, new_observation = patient_1_data, type = "survshap")

plot(survshap_t)

# Optional packages can change their accepted prediction representation across
# releases, so keep this complementary illustration failure-tolerant.
survex_perf <- tryCatch(
  model_performance(surv_explainer),
  error = function(e) {
    message("The installed survex version could not evaluate this explainer: ",
            conditionMessage(e))
    NULL
  }
)

# Plot the Brier score and AUC curves
if (!is.null(survex_perf)) {
  plot(survex_perf)
}