Understanding SHAP one path at a time

Build a SHAP explanation from the predictions produced while transforming a background observation into the observation we want to explain.
R
machine-learning
explainability
Author

Joshua Kunst

Published

August 14, 2026

SHAP is usually introduced through Shapley values, coalitions and game theory. All of that matters, but I find it easier to first look at the actual rows that a model sees.

The idea in this post is simple: start from one background observation, transform it into the observation we want to explain one variable at a time, and watch how the model prediction changes along the way.

A model and a prediction

Assume we already have a credit-risk model, a training sample, a test sample and a vector with the predictor names. The details of the model are not important for what follows, so the preparation code is folded below.

Show model preparation
library(dplyr)
library(ggbeeswarm)
library(ggplot2)
library(modeldata)
library(tibble)
library(tidyr)
library(xgboost)

predictors <- c(
  "seniority", "time", "age", "expenses", "income",
  "assets", "debt", "amount", "price"
)

credit_data <- modeldata::credit_data |>
  tibble::as_tibble() |>
  dplyr::rename_with(tolower) |>
  tidyr::drop_na(status, dplyr::all_of(predictors)) |>
  dplyr::mutate(status_bad = as.integer(status == "bad")) |>
  dplyr::select(status_bad, dplyr::all_of(predictors))

set.seed(2026)
train_id <- sample.int(nrow(credit_data), floor(0.75 * nrow(credit_data)))

train <- credit_data[train_id, ]
test <- credit_data[-train_id, ]

dtrain <- xgboost::xgb.DMatrix(
  as.matrix(train[, predictors]),
  label = train$status_bad
)

set.seed(2026)
model <- xgboost::xgb.train(
  data = dtrain,
  nrounds = 150,
  params = list(
    objective = "binary:logistic",
    eval_metric = "logloss",
    max_depth = 4,
    eta = 0.05
  ),
  verbose = 0
)

predict_model <- function(model, newdata) {
  stats::predict(
    model,
    xgboost::xgb.DMatrix(as.matrix(newdata[, predictors]))
  )
}

Let’s take one observation from the test sample. This is the prediction we want to explain.

Code
x <- test |>
  dplyr::slice(1) |>
  dplyr::select(dplyr::all_of(predictors))

prediction_x <- predict_model(model, x)

Before doing anything else, let’s look at the actual values of the observation.

Code
x |>
  dplyr::mutate(prediction = prediction_x) |>
  knitr::kable(digits = 4)
seniority time age expenses income assets debt amount price prediction
1 60 36 75 214 3500 0 650 1645 0.1047

The model gives this observation a probability of 10.5%.

Our question is simply: where does this prediction come from?

Start from a background observation

SHAP needs a reference. Let’s take a small sample from the training data as our background and, for the moment, choose only one observation from it.

Code
set.seed(2026)
background <- train |>
  dplyr::slice_sample(n = 50) |>
  dplyr::select(dplyr::all_of(predictors))

z <- background |>
  dplyr::slice(1)

prediction_z <- predict_model(model, z)

One background row is useful for understanding a single path, but it is not a privileged reference. The eventual explanation uses many background observations, each offering a different route toward the same target x.

Code
set.seed(2026)

background_diagram <- tibble::tibble(
  x = runif(18, 0.05, 0.32),
  y = runif(18, 0.08, 0.92),
  xend = 0.82,
  yend = 0.50
)

ggplot(background_diagram) +
  geom_segment(
    aes(x, y, xend = xend, yend = yend),
    color = scales::alpha("#3487d4", 0.23),
    linewidth = 0.55,
    arrow = grid::arrow(length = grid::unit(0.09, "inches"))
  ) +
  geom_point(aes(x, y), color = "#3487d4", size = 2.6) +
  geom_point(aes(xend, yend), color = "#17324d", size = 7) +
  annotate(
    "text", x = 0.18, y = 0.99,
    label = "background observations  z₁, …, zₙ",
    color = "#52677b", size = 4
  ) +
  annotate(
    "text", x = 0.82, y = 0.62,
    label = "target  x",
    color = "#17324d", fontface = "bold", size = 4.5
  ) +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1), clip = "off") +
  theme_void()

This is only a conceptual diagram—the real paths live in predictor space—but it captures the role of the background sample: many plausible starting points are transformed toward the same observation we want to explain.

Now we can put the observation we want to explain, x, next to the background observation, z.

Code
comparison <- dplyr::bind_rows(
  target = x,
  background = z,
  .id = "observation"
)

comparison$prediction <- predict_model(model, comparison)

comparison |>
  knitr::kable(digits = 4)
observation seniority time age expenses income assets debt amount price prediction
target 1 60 36 75 214 3500 0 650 1645 0.1047
background 1 12 61 35 58 18000 0 500 1342 0.0460

We have two rows, two sets of values and two model predictions. Instead of jumping directly from z to x, let’s move from one to the other one variable at a time.

One path from z to x

First choose a random order for the predictors.

Code
set.seed(1)
order_1 <- sample(predictors)
order_1
[1] "price"     "expenses"  "debt"      "seniority" "time"      "assets"   
[7] "age"       "amount"    "income"

Now start from z. Replace the first variable with its value in x, then the second one, then the third one, and continue until every value comes from x.

Code
build_path <- function(x, z, variable_order) {
  current <- z
  states <- list(current)

  for (variable in variable_order) {
    current[[variable]] <- x[[variable]]
    states[[length(states) + 1L]] <- current
  }

  dplyr::bind_rows(states) |>
    dplyr::mutate(
      step = 0:length(variable_order),
      changed = c("background", variable_order),
      .before = 1
    )
}

With nine predictors we get ten rows: the original background observation, one row for each variable replacement, and finally x itself.

Code
path_1 <- build_path(x, z, order_1)
path_1$prediction <- predict_model(model, path_1)

If we display the columns in the same order in which they were introduced, the construction has a staircase-like shape. At each row, one more value has moved from z to x.

Code
path_1 |>
  dplyr::select(step, changed, dplyr::all_of(order_1), prediction) |>
  knitr::kable(digits = 4)
step changed price expenses debt seniority time assets age amount income prediction
0 background 1342 35 0 1 12 18000 61 500 58 0.0460
1 price 1645 35 0 1 12 18000 61 500 58 0.0460
2 expenses 1645 75 0 1 12 18000 61 500 58 0.0582
3 debt 1645 75 0 1 12 18000 61 500 58 0.0582
4 seniority 1645 75 0 1 12 18000 61 500 58 0.0582
5 time 1645 75 0 1 60 18000 61 500 58 0.1939
6 assets 1645 75 0 1 60 3500 61 500 58 0.3100
7 age 1645 75 0 1 60 3500 36 500 58 0.3628
8 amount 1645 75 0 1 60 3500 36 650 58 0.3784
9 income 1645 75 0 1 60 3500 36 650 214 0.1047

The model does not know anything about SHAP here. It simply predicts ten observations.

Because only one variable changes between consecutive rows, the difference between two consecutive predictions can be assigned to the variable that just changed.

Code
contributions_1 <- path_1 |>
  dplyr::transmute(
    variable = changed,
    prediction,
    contribution = prediction - dplyr::lag(prediction)
  ) |>
  dplyr::filter(variable != "background")

contributions_1 |>
  knitr::kable(digits = 4)
variable prediction contribution
price 0.0460 0.0000
expenses 0.0582 0.0122
debt 0.0582 0.0000
seniority 0.0582 0.0000
time 0.1939 0.1357
assets 0.3100 0.1161
age 0.3628 0.0528
amount 0.3784 0.0157
income 0.1047 -0.2737

For this path, if replacing debt changes the prediction from 0.10 to 0.16, then the contribution of debt along this path is simply +0.06.

Nothing more sophisticated has happened yet: we are just taking differences between consecutive predictions.

Same observations, different path

There is an immediate problem: the contribution of a variable can depend on what was changed before it.

Let’s keep exactly the same x and the same z and change only the order of the variables.

Code
set.seed(2)
order_2 <- sample(predictors)
order_2
[1] "income"    "debt"      "assets"    "amount"    "seniority" "price"    
[7] "expenses"  "time"      "age"

Build and predict the second path.

Code
path_2 <- build_path(x, z, order_2)
path_2$prediction <- predict_model(model, path_2)

Again, the first row is exactly z and the last row is exactly x, but the observations in between are different.

Code
path_2 |>
  dplyr::select(step, changed, dplyr::all_of(order_2), prediction) |>
  knitr::kable(digits = 4)
step changed income debt assets amount seniority price expenses time age prediction
0 background 58 0 18000 500 1 1342 35 12 61 0.0460
1 income 214 0 18000 500 1 1342 35 12 61 0.0156
2 debt 214 0 18000 500 1 1342 35 12 61 0.0156
3 assets 214 0 3500 500 1 1342 35 12 61 0.0316
4 amount 214 0 3500 650 1 1342 35 12 61 0.0390
5 seniority 214 0 3500 650 1 1342 35 12 61 0.0390
6 price 214 0 3500 650 1 1645 35 12 61 0.0387
7 expenses 214 0 3500 650 1 1645 75 12 61 0.0429
8 time 214 0 3500 650 1 1645 75 60 61 0.0957
9 age 214 0 3500 650 1 1645 75 60 36 0.1047

And therefore the intermediate changes can also be different.

Code
contributions_2 <- path_2 |>
  dplyr::transmute(
    variable = changed,
    prediction,
    contribution = prediction - dplyr::lag(prediction)
  ) |>
  dplyr::filter(variable != "background")

path_comparison <- dplyr::bind_rows(
  path_1 = contributions_1,
  path_2 = contributions_2,
  .id = "path"
) |>
  dplyr::select(path, variable, contribution) |>
  tidyr::pivot_wider(
    names_from = path,
    values_from = contribution
  ) |>
  dplyr::mutate(difference = path_2 - path_1)

path_comparison |>
  knitr::kable(digits = 4)
variable path_1 path_2 difference
price 0.0000 -0.0003 -0.0004
expenses 0.0122 0.0042 -0.0080
debt 0.0000 0.0000 0.0000
seniority 0.0000 0.0000 0.0000
time 0.1357 0.0527 -0.0829
assets 0.1161 0.0160 -0.1001
age 0.0528 0.0090 -0.0438
amount 0.0157 0.0074 -0.0083
income -0.2737 -0.0304 0.2434

Both paths start at the same prediction, p(z), and finish at the same prediction, p(x). What changes is how the total difference is distributed among the variables.

This is the first key idea behind SHAP:

There is no single privileged order in which the variables should enter the prediction.

For nonlinear models and models with interactions, the context in which a variable is introduced matters.

From one path to many paths

Trying every possible order quickly becomes impossible. With nine variables there are already 9! = 362,880 possible orders.

Instead, we can sample paths.

There is another arbitrary choice in our example: the background observation itself. Why should one particular z define our reference?

So we repeat the same experiment over many observations from the background. For each background observation we draw a random variable order, build the path to x, predict every row and keep the consecutive differences.

Code
one_path <- function(model, x, z, variable_order) {
  path <- build_path(x, z, variable_order)
  path$prediction <- predict_model(model, path)

  path |>
    dplyr::transmute(
      variable = changed,
      contribution = prediction - dplyr::lag(prediction)
    ) |>
    dplyr::filter(variable != "background")
}
Code
set.seed(2026)

trace <- dplyr::bind_rows(lapply(seq_len(nrow(background)), function(i) {
  variable_order <- sample(predictors)

  one_path(
    model = model,
    x = x,
    z = background[i, ],
    variable_order = variable_order
  ) |>
    dplyr::mutate(background_id = i, .before = 1)
}))

glimpse(trace)
Rows: 450
Columns: 3
$ background_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,…
$ variable      <chr> "price", "seniority", "assets", "income", "age", "expens…
$ contribution  <dbl> 4.982203e-05, 0.000000e+00, 3.571071e-02, -5.034730e-02,…

Now every variable has many observed contributions, obtained under different starting points and different contexts.

The SHAP approximation is simply their average.

Code
shap_values <- trace |>
  dplyr::summarise(
    shap = mean(contribution),
    .by = variable
  ) |>
  dplyr::arrange(desc(abs(shap)))

shap_values |>
  knitr::kable(digits = 4)
variable shap
amount -0.0959
income -0.0912
seniority 0.0630
price -0.0358
assets -0.0173
time 0.0144
debt -0.0122
age 0.0095
expenses 0.0072

This is a Monte Carlo approximation of marginal SHAP: instead of enumerating every possible permutation, we sample paths and average the contribution assigned to each variable.

Adding the contributions

There is a useful property hiding in plain sight.

The official SHAP introduction shows the same idea with a waterfall: begin at the expected model output and add the feature contributions until reaching the prediction for one observation.1 Our path construction lets us see where those additions come from. For an R-focused explanation, the Explanatory Model Analysis guide shows how DALEX::predict_parts() decomposes one prediction into contributions from individual variables.2

For any single path, the contributions are consecutive differences:

\[ [p(s_1)-p(s_0)] + [p(s_2)-p(s_1)] + \cdots + [p(x)-p(s_{p-1})]. \]

Everything in the middle cancels, leaving

\[ p(x) - p(z). \]

Averaging over the background observations gives

\[ p(x) = E[p(z)] + \sum_j \phi_j. \]

We can check that directly.

Code
background_mean <- mean(predict_model(model, background))
reconstructed <- background_mean + sum(shap_values$shap)

reconstruction <- tibble::tibble(
  background_mean = background_mean,
  shap_sum = sum(shap_values$shap),
  reconstructed = reconstructed,
  prediction = prediction_x
)

reconstruction |>
  knitr::kable(digits = 6)
background_mean shap_sum reconstructed prediction
0.263 -0.158304 0.104695 0.104695

This is why the familiar SHAP waterfall works: start from the background prediction, add the contribution of each variable, and arrive at the prediction for the observation we wanted to explain.

From SHAP values to a waterfall

Now that we have the pieces, we can draw the same type of waterfall used in the SHAP Explorer app.

The first bar is the mean prediction over the background sample. Each variable then moves the prediction up or down by its SHAP value. Red contributions increase the predicted probability of default, blue contributions decrease it, and the final bar is the prediction for x.

The order of the bars here is only a display choice: the SHAP values have already been averaged across the sampled paths.

Code
library(highcharter)

profile_labels <- c(
  seniority = "Seniority",
  time = "Loan term",
  age = "Age",
  expenses = "Expenses",
  income = "Income",
  assets = "Assets",
  debt = "Debt",
  amount = "Loan amount",
  price = "Price"
)

contribution_points <- lapply(seq_len(nrow(shap_values)), function(i) {
  value <- 100 * shap_values$shap[[i]]
  variable <- shap_values$variable[[i]]

  list(
    name = unname(profile_labels[[variable]]),
    y = unname(value),
    color = if (value >= 0) "#d95f59" else "#4c91d9"
  )
})

waterfall_data <- unname(c(
  list(list(
    name = "Background<br/>mean",
    y = 100 * background_mean,
    color = "#e0e0e0"
  )),
  contribution_points,
  list(list(
    name = "Predicted<br/>PD",
    isSum = TRUE,
    color = "#34495e"
  ))
))

highchart() |>
  hc_add_dependency("modules/waterfall.js") |>
  hc_chart(type = "waterfall") |>
  hc_xAxis(type = "category") |>
  hc_yAxis(title = list(text = "Probability of default (%)")) |>
  hc_legend(enabled = FALSE) |>
  hc_tooltip(pointFormat = "{point.y:.1f}") |>
  hc_plotOptions(series = list(
    borderWidth = 0,
    dataLabels = list(
      enabled = TRUE,
      inside = FALSE,
      useHTML = TRUE,
      style = list(
        color = "#495057",
        fontWeight = "normal",
        textOutline = "none"
      ),
      formatter = JS(paste(
        "function () {",
        "  if (this.point.isSum || this.point.index === 0) {",
        "    const value = Highcharts.numberFormat(this.y, 1) + '%';",
        "    return '<span style=\"font-size: 13px; font-weight: 600\">' + value + '</span>';",
        "  }",
        "  return (this.y >= 0 ? '+' : '') + Highcharts.numberFormat(this.y, 1) + ' pp';",
        "}"
      ))
    )
  )) |>
  hc_add_series(
    name = "PD",
    data = waterfall_data
  )

The chart is just another view of the identity we checked above:

\[ \text{background mean} + \sum_j \phi_j = p(x). \]

The whole idea

Before moving to the summary, it is worth distinguishing two familiar SHAP views. The waterfall above is local: it explains one prediction. A beeswarm is global: every point represents one observation’s contribution for one variable. Horizontal position shows whether that variable moved the prediction down or up; color shows whether the original feature value was relatively low or high. The vertical spreading only prevents overlapping points and reveals the distribution.

Code
set.seed(2026)

global_targets <- test |>
  dplyr::slice_sample(n = 60) |>
  dplyr::select(dplyr::all_of(predictors))

global_trace <- dplyr::bind_rows(lapply(seq_len(nrow(global_targets)), function(target_id) {
  target_x <- global_targets[target_id, , drop = FALSE]
  background_ids <- sample.int(nrow(background), 15, replace = TRUE)

  dplyr::bind_rows(lapply(background_ids, function(background_id) {
    one_path(
      model = model,
      x = target_x,
      z = background[background_id, , drop = FALSE],
      variable_order = sample(predictors)
    )
  })) |>
    dplyr::summarise(
      shap = mean(contribution),
      .by = variable
    ) |>
    dplyr::mutate(target_id = target_id)
}))

feature_values <- global_targets |>
  dplyr::mutate(target_id = dplyr::row_number()) |>
  tidyr::pivot_longer(
    cols = dplyr::all_of(predictors),
    names_to = "variable",
    values_to = "feature_value"
  )

global_shap <- global_trace |>
  dplyr::left_join(feature_values, by = c("target_id", "variable")) |>
  dplyr::mutate(
    scaled_value = scales::rescale(feature_value, to = c(0, 1)),
    .by = variable
  ) |>
  dplyr::mutate(
    variable = reorder(variable, abs(shap), FUN = mean)
  )
Code
ggplot(global_shap, aes(shap, variable, color = scaled_value)) +
  geom_vline(xintercept = 0, color = "#8a96a3", linewidth = 0.5) +
  ggbeeswarm::geom_quasirandom(
    groupOnX = FALSE,
    width = 0.32,
    size = 1.7,
    alpha = 0.82
  ) +
  scale_color_gradientn(
    colors = c("#1687e8", "#7b3fc6", "#f40064"),
    limits = c(0, 1),
    breaks = c(0, 1),
    labels = c("Low", "High")
  ) +
  labs(
    title = "How variables move predictions across observations",
    subtitle = "Each point is one approximate local SHAP value",
    x = "SHAP value (impact on predicted probability)",
    y = NULL,
    color = "Feature value"
  ) +
  theme(
    panel.grid.major.y = element_line(color = "#e4e9ee", linetype = 3),
    legend.position = "right"
  )

The complete calculation can be summarized without game theory notation:

  1. Pick an observation x to explain.
  2. Pick a background observation z.
  3. Randomly order the predictors.
  4. Transform z into x, one variable at a time.
  5. Predict every intermediate row.
  6. Take differences between consecutive predictions.
  7. Repeat for more background observations and variable orders.
  8. Average the differences by variable.

That average is our SHAP approximation.

I like this construction because the explanation emerges from objects we already know how to inspect: rows, predictions and differences. The formal Shapley framework tells us why averaging across different contexts is the right thing to do, but the staircase of observations shows what the model is actually being asked to evaluate.

The helper one_path() exists only to package the repeated unit of work: build one path, predict its intermediate states and return the consecutive differences. Keeping that operation named makes the sampling loop readable without hiding the algorithm.

The same idea is used in my SHAP Explorer, where the readable implementation is intentionally kept next to the optimized version. The optimized code is faster; the slow version is there because it makes the algorithm easier to see.

Try it yourself

The SHAP Explorer below uses the same construction interactively. Change the client profile or the model and watch how the local contributions, the waterfall and the predicted probability move together.

Footnotes

  1. See the official SHAP guide, An introduction to explainable AI with Shapley values.↩︎

  2. See Break-down Plots for Additive Attributions in Explanatory Model Analysis.↩︎