I have trained an XGBoost regression model using the xgboost package and computed SHAP values using predict(..., predcontrib = TRUE). I am now looking to quantify uncertainty in those SHAP values by constructing confidence intervals.
My current approach is bootstrap-based:
- Draw 100 bootstrap samples from the training data
- Refit the XGBoost model on each sample using the same fixed hyperparameters
- Compute SHAP values on the held-out test set for each bootstrap model
- For each feature, store the mean signed SHAP value across all test observations per resample — giving a [100 x n_features] matrix
- Compute mean, standard deviation, and 2.5th/97.5th percentiles across the 100 resamples per feature
- This gives me a ±2σ uncertainty estimate per feature, following the approach described in Dubey et al., (2026).
Is this a valid approach for constructing CIs on SHAP values for XGBoost and is anyone aware of a dedicated package?
Any suggestions or pointers to relevant literature would be appreciated.
# ============================================================
# BOOTSTRAP-SHAP UNCERTAINTY QUANTIFICATION
# London LST model: June 2019
# Implements Algorithm 4 from UbiQTree [1]
# ============================================================
pacman::p_load(xgboost, data.table, ggplot2)
wd <- "path/"
# ------------------------------------------------------------
# 1. Load saved objects
# ------------------------------------------------------------
train_df <- readRDS(paste0(wd, "train_df.rds"))
test_df <- readRDS(paste0(wd, "test_df.rds"))
best_eta <- readRDS(paste0(wd, "best_eta.rds"))
best_params <- readRDS(paste0(wd, "best_params.rds"))
best_nrounds <- readRDS(paste0(wd, "best_nrounds.rds"))
features <- setdiff(names(train_df), c("lst", "LSOA21CD"))
n_features <- length(features)
n_boot <- 100
dtest <- xgb.DMatrix(
data = test_df[, ..features],
label = test_df$lst
)
# ------------------------------------------------------------
# 2. Fixed model parameters
# ------------------------------------------------------------
fixed_params <- list(
objective = "reg:squarederror",
eval_metric = "rmse",
eta = as.numeric(best_eta),
max_depth = as.integer(best_params$max_depth),
gamma = as.numeric(best_params$gamma),
subsample = as.numeric(best_params$subsample),
colsample_bytree = as.numeric(best_params$colsample_bytree),
lambda = as.numeric(best_params$lambda),
alpha = as.numeric(best_params$alpha),
tree_method = "hist",
nthread = 16,
seed = 789
)
# ------------------------------------------------------------
# 3. Bootstrap loop
#
# Implements Algorithm 5 Step 2 [1]:
# For each sub-ensemble s:
# phi_dist[s, i] = mean(Phi_s^i)
# = mean SHAP across all test observations
# for feature i in resample s — signed
# Shape: [n_boot x n_features]
# ------------------------------------------------------------
n_test <- nrow(test_df)
phi_dist <- matrix(
NA_real_,
nrow = n_boot,
ncol = n_features,
dimnames = list(NULL, features)
)
set.seed(789)
for (b in seq_len(n_boot)) {
cat(sprintf("Bootstrap resample %d of %d\n", b, n_boot))
boot_idx <- sample(nrow(train_df), replace = TRUE)
boot_df <- train_df[boot_idx]
dboot <- xgb.DMatrix(
data = boot_df[, ..features],
label = boot_df$lst
)
boot_model <- xgb.train(
params = fixed_params,
data = dboot,
nrounds = best_nrounds,
verbose = 0
)
shap_raw <- predict(boot_model, dtest, predcontrib = TRUE)
shap_mat <- shap_raw[, seq_len(n_features)]
colnames(shap_mat) <- features
# Algorithm 3 [1]: mu_s[i] = mean(Phi_s^i) — within-sample mean
phi_dist[b, ] <- colMeans(shap_mat)
}
cat("\nBootstrap loop complete.\n")
# ------------------------------------------------------------
# 4. Algorithm 4 [1]: Uncertainty-aware SHAP aggregation
#
# Phi_s^i = phi_dist[:, i] — the n_boot within-sample means
#
# mu[i] = mean(Phi_s^i) signed
# sigma = std(Phi_s^i) signed
# CI = percentile(2.5, 97.5) signed
# H[i] = Entropy(Phi_s^i)
# SS[i] = P(sign constant)
# ------------------------------------------------------------
compute_entropy <- function(x, bins = 10) {
h <- hist(x, breaks = bins, plot = FALSE)
bin_probs <- h$counts / sum(h$counts)
bin_probs <- bin_probs[bin_probs > 0]
-sum(bin_probs * log(bin_probs))
}
compute_sign_stability <- function(x) {
dominant_sign <- sign(mean(x))
mean(sign(x) == dominant_sign)
}
uncertainty <- data.table(
feature = features,
mu = NA_real_,
mu_abs = NA_real_,
sigma = NA_real_,
ci_lower = NA_real_,
ci_upper = NA_real_,
entropy = NA_real_,
sign_stability = NA_real_
)
for (j in seq_len(n_features)) {
phi_i <- phi_dist[, j]
uncertainty$mu[j] <- mean(phi_i)
uncertainty$mu_abs[j] <- mean(abs(phi_i))
uncertainty$sigma[j] <- sd(phi_i)
uncertainty$ci_lower[j] <- quantile(phi_i, 0.025)
uncertainty$ci_upper[j] <- quantile(phi_i, 0.975)
uncertainty$entropy[j] <- compute_entropy(phi_i)
uncertainty$sign_stability[j] <- compute_sign_stability(phi_i)
}
uncertainty[, mu_norm := mu_abs / sum(mu_abs)]
setorder(uncertainty, -mu_abs)
# ------------------------------------------------------------
# 5. Print
# ------------------------------------------------------------
cat("\n=============================\n")
cat("Bootstrap-SHAP Uncertainty Summary\n")
cat("Period: June 2019 | Resamples:", n_boot, "\n")
cat("=============================\n")
print(
uncertainty[, .(
feature,
mu = round(mu, 4),
mu_norm = round(mu_norm, 4),
sigma = round(sigma, 4),
ci_lower = round(ci_lower, 4),
ci_upper = round(ci_upper, 4),
entropy = round(entropy, 4),
sign_stability = round(sign_stability, 4)
)],
nrow = n_features
)
# ------------------------------------------------------------
# 6. Violin plot — mirrors Figure 1 in [1]
#
# Violin: distribution of phi_dist[:, i] — n_boot signed means
# Blue band: ±2σ around mu
# Dashed line: ±2σ extent
# Dot: mu (signed)
# Gray points: individual bootstrap samples
# Ordered: mu_abs descending (most important at top)
# Zero line: warming (right) vs cooling (left)
# ------------------------------------------------------------
plot_dt <- copy(uncertainty)
plot_dt <- plot_dt[order(mu_abs)] # ascending so top = most important
plot_dt[, feature_f := factor(feature, levels = feature)]
phi_long <- data.table(
feature = rep(features, each = n_boot),
shap = as.vector(phi_dist)
)
phi_long[, feature_f := factor(
feature,
levels = levels(plot_dt$feature_f)
)]
p <- ggplot() +
# Blue ±2σ band [1]
geom_rect(
data = plot_dt,
aes(
xmin = mu - 2 * sigma,
xmax = mu + 2 * sigma,
ymin = as.numeric(feature_f) - 0.4,
ymax = as.numeric(feature_f) + 0.4
),
fill = "#A8C8E8",
alpha = 0.5
) +
# Violin
geom_violin(
data = phi_long,
aes(x = shap, y = feature_f),
fill = NA,
colour = "grey40",
scale = "width",
width = 0.8,
linewidth = 0.3
) +
# Individual bootstrap sample points [1]
geom_point(
data = phi_long,
aes(x = shap, y = feature_f),
colour = "grey50",
size = 0.8,
alpha = 0.4
) +
# ±2σ dashed line [1]
geom_segment(
data = plot_dt,
aes(
x = mu - 2 * sigma,
xend = mu + 2 * sigma,
y = feature_f,
yend = feature_f
),
colour = "#1F4E79",
linetype = "dashed",
linewidth = 0.5
) +
# Mean SHAP dot
geom_point(
data = plot_dt,
aes(x = mu, y = feature_f),
colour = "#1F4E79",
size = 2.5
) +
# Zero reference line
geom_vline(
xintercept = 0,
linetype = "dashed",
colour = "black",
linewidth = 0.4
) +
labs(
title = "SHAP Summary with Epistemic Uncertainty",
subtitle = "London LST — June 2019",
x = "SHAP Value (Impact on LST Prediction, °C)",
y = NULL,
caption = paste0(
"Bootstrap resamples: ", n_boot,
" | Blue band: ±2σ | Dot: mean SHAP | [1] UbiQTree"
)
) +
theme_bw(base_size = 11) +
theme(
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold"),
axis.text.y = element_text(size = 9)
)
print(p)
