Computing SHAP values for rule
ensembles
Marjolein Fokkema
Introduction
SHAP (SHapley Additive exPlanations) values were introduced by Lundberg & Lee (2017) and quantify the
contribution of individual predictors to a model’s predictions. SHAP
values quantify these contributions for specific observations, and
therefore provide measures of local importance Molnar (2025). Although exact SHAP values can be
computationally expensive to obtain, in Spadaccini, Fokkema, & Wiel (2025) we
derived a simplified expression for computing exact marginal SHAP values
for prediction rule ensembles. This functionality is implemented in
function shap. This vignette illustrates its use, and how
the results can be plotted using package shapviz.
Example: Fitting a rule ensemble
We fit a PRE to predict Ozone (daily ozone readings)
with the airquality data:
airq <- na.omit(airquality)
set.seed(42)
airq.ens <- pre(Ozone ~ ., data = airq)
airq.ens ## equivalent to print(airq.ens)
##
## Final ensemble with cv error within 1se of minimum:
##
## lambda = 3.543968
## number of terms = 12
## mean cv error (se) = 352.395 (99.13754)
##
## cv error type : Mean-Squared Error
##
## rule coefficient description
## (Intercept) 68.48270407 1
## rule191 -10.97368180 Wind > 5.7 & Temp <= 87
## rule173 -10.90385519 Wind > 5.7 & Temp <= 82
## rule42 -8.79715538 Wind > 6.3 & Temp <= 84
## rule204 7.16114781 Wind <= 10.3 & Solar.R > 148
## rule10 -4.68646145 Temp <= 84 & Temp <= 77
## rule192 -3.34460038 Wind > 5.7 & Temp <= 87 & Day <= 23
## rule51 -2.27864287 Wind > 5.7 & Temp <= 84
## rule93 2.18465676 Temp > 77 & Wind <= 8.6
## rule74 -1.36479545 Wind > 6.9 & Temp <= 84
## rule28 -1.15326093 Temp <= 84 & Wind > 7.4
## rule25 -0.70818400 Wind > 6.3 & Temp <= 82
## rule166 -0.04751152 Wind > 6.9 & Temp <= 82
Twelve rules were selected, involving predictors Temp,
Wind, Solar.R and Day. No linear
terms were selected.
By default, the 1SE rule is used to select the final ensemble by the
print method, which can be overridden by specifying the
penalty.par.val argument.
Computing marginal SHAP values
We compute SHAP values using function shap (it uses the
same default for the penalty.par.val argument, which can be
overridden):
airq.shaps <- shap(airq.ens, newdata = airq)
The newdata argument specifies for which observations
SHAP values should be computed. This can be one or more new (or test)
observations, but here we requested SHAP values for the full training
dataset. Not specifying the newdata argument would have
yielded the exact same result.
The result is a list of two elements, marginal and
interactions:
## List of 2
## $ marginal : num [1:111, 1:5] 1.839 -3.935 0.677 0.677 1.839 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:111] "1" "2" "3" "4" ...
## .. ..$ : chr [1:5] "Solar.R" "Wind" "Temp" "Month" ...
## $ interactions: NULL
Marginal SHAP values are computed for each predictor, by default. The
marginal SHAP values are returned as an \(N
\times p\) matrix, with SHAP values for each observation (rows)
and predictor (columns). Interaction SHAP values will be computed when
interactions = TRUE is specified.
Visualizing marginal SHAP values with shapviz
We use R package shapviz to visualize the SHAP values.
We need to construct a shapviz object from the SHAP values,
predictor variable values and a baseline:
library("shapviz")
baseline <- mean(predict(airq.ens, newdata = airq))
sv <- shapviz(object = airq.shaps$marginal,
X = airq,
baseline = baseline)
With the default reference data used by function shap
(see also argument reference_data for changing it), the
baseline is the mean prediction over the the training data. Furthermore,
argument X specifies the original predictor values
corresponding to the rows of the SHAP matrix. Note that categorical
predictors (if present) should be supplied in their original factor
representation.
From this shapviz object, we can obtain e.g. beeswarm
plots:
sv_importance(sv, kind = "bee")
The predictors are ordered by the variation in their SHAP values.
More variation corresponds to greater magnitude, because a variable that
does not contribute to predictions would have SHAP values of 0.
Temp thus appears most important, followed by
Wind, then by Solar.R, then by
Day which has a very minor effect and finally
Month which has no effect on predictions.
The contributions of Temp and Solar.R
appear monotonically increasing, the effect of Wind
monotonically decreasing. If observations with similar feature values
(and therefore similar colors) have substantially different SHAP values,
this suggests interactions with other predictors, which does not seem to
be the case here.
Month does not affect predictions, but its points are
colored mostly yellow because the observations in airq are
ordered by Month. The apparent color pattern is therefore
an artifact of the observation order, which is also used for
plotting.
Often, waterfall plots are used for explaining individual
predictions, e.g.:
sv_waterfall(sv, row_id = 1)
The waterfall starts at the baseline prediction and successively adds
the SHAP contributions to arrive at the model prediction. For the first
observation in the dataset, Solar.R and Temp
have small positive effects and Day has a small negative
effect. We observe zero contribution for Month, because it
does not appear in any of the terms of the rule ensemble.
Other functions from package shapviz, like
sv_importance and sv_dependence can also be
applied to sv.
Computing interaction SHAP values
We can obtain SHAP interaction values by specifying
interactions = TRUE:
airq.shaps <- shap(airq.ens, newdata = airq, interactions = TRUE)
str(airq.shaps)
## List of 2
## $ marginal : num [1:111, 1:5] 1.839 -3.935 0.677 0.677 1.839 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:111] "1" "2" "3" "4" ...
## .. ..$ : chr [1:5] "Solar.R" "Wind" "Temp" "Month" ...
## $ interactions: num [1:111, 1:5, 1:5] 1.35 -3.03 1.35 1.35 1.35 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:111] "1" "2" "3" "4" ...
## .. ..$ : chr [1:5] "Solar.R" "Wind" "Temp" "Month" ...
## .. ..$ : chr [1:5] "Solar.R" "Wind" "Temp" "Month" ...
Interaction SHAP values are saved in a three-dimensional \(N \times p \times p\) array. That is, for
each of the \(N\) observations, a \(p \times p\) matrix of SHAP values is
returned. The diagonal entries are predictor variables’ main-effect SHAP
values, and the off-diagonal entries represent their interaction SHAP
values. E.g., for the first observation the SHAP interaction matrix
looks as follows:
airq.shaps$interactions[1, , ]
## Solar.R Wind Temp Month Day
## Solar.R 1.3548117 0.48386134 0.00000000 0 0.00000000
## Wind 0.4838613 1.43358614 -1.77812376 0 -0.03766442
## Temp 0.0000000 -1.77812376 -8.98161032 0 -0.06779595
## Month 0.0000000 0.00000000 0.00000000 0 0.00000000
## Day 0.0000000 -0.03766442 -0.06779595 0 -0.56747724
The diagonal entries tend to have stronger magnitude than the
off-diagonal entries, indicating that the model’s effects for this
observation are mostly main effects, while interactions are minor. There
is one exception to this rule: The interaction SHAP for
Wind and Temp, which has stronger magnitude
than the main effect of Wind (but weaker than the main
effect of Temp) for this observation.
Understanding (interaction) SHAP values
Interaction SHAP values decompose the marginal SHAP values into
contributions due to main effects and to interactions. Because the
matrix of interaction SHAP values is symmetric, the row sums and column
sums are identical, and they are also identical to the marginal SHAP
values. E.g. for the first observation we have:
rowSums(airq.shaps$interactions[1, , ])
## Solar.R Wind Temp Month Day
## 1.8386731 0.1016593 -10.8275300 0.0000000 -0.6729376
colSums(airq.shaps$interactions[1, , ])
## Solar.R Wind Temp Month Day
## 1.8386731 0.1016593 -10.8275300 0.0000000 -0.6729376
## Solar.R Wind Temp Month Day
## 1.8386731 0.1016593 -10.8275300 0.0000000 -0.6729376
SHAP values represent the predictor variables’ contribution to the
predicted value, so by definition they sum to the difference between the
model’s predicted value and the baseline prediction. More formally, a
model’s prediction \(f(x)\) is given
by:
\[f(x) = E[f(X)] + \sum_j
\phi_j(x),\] where \(f\) is the
predictive model, \(x\) is a vector of
predictor variable values, \(E[f(X)]\)
is the expected (or mean) model prediction over the reference
distribution and \(\phi_j(x)\) is the
SHAP value for predictor \(j\).
The SHAP values \(\phi_j\) can be
further decomposed into a main-effect contribution (diagonal elements)
and interaction contributions with all other predictors (the
corresponding off-diagonal elements). More formally,
\[ \phi_j(x) = \phi_{jj}(x) + \sum_{k\neq
j}\phi_{jk}(x), \]
where \(\phi_{jj}\) is the main
effect (diagonal entries of the SHAP interaction matrix) and \(\phi_{jk}\) are pairwise interactions
(off-diagonal entries of the SHAP interaction matrix).
Visualizing interaction SHAP values
From the computed (interaction) SHAP values, we again create a
shapviz object:
baseline <- mean(predict(airq.ens, newdata = airq))
sv_int <- shapviz(object = airq.shaps$marginal,
X = airq,
baseline = baseline,
S_inter = airq.shaps$interactions)
sv_interaction(sv_int)
The plots on the diagonal show greatest variation in SHAP values,
indicating again that main effects appear strongest in the model.
Evidence for interaction effects would be found in off-diagonal plots
showing variation in SHAP values.
We again observe that Month does not contribute to
predictions, indicated by all SHAP values being zero.
Solar.R and Day contribute only somewhat.
Wind and Temp show the strongest effects, and
there might be some interaction between them. Note that the two plots
for this interaction are identical, save for the coloring, which
indicates the values of Wind in the 2nd row, 1st column,
and the values of Temp in the 1st row, 2nd column. However,
the direction and form of this interaction are not easy to discern from
these plots alone.
Session info
In case you obtained different results, the results above were
obtained using the following:
## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=C LC_CTYPE=Dutch_Netherlands.utf8
## [3] LC_MONETARY=Dutch_Netherlands.utf8 LC_NUMERIC=C
## [5] LC_TIME=Dutch_Netherlands.utf8
##
## time zone: Europe/Amsterdam
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] shapviz_0.10.3 pROC_1.19.0.1 caret_7.0-1 lattice_0.22-9 ggplot2_4.0.3
## [6] pre_1.1.0 mice_3.19.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.6 deldir_2.0-4 rlang_1.2.0
## [4] magrittr_2.0.4 otel_0.2.0 compiler_4.5.3
## [7] vctrs_0.7.2 reshape2_1.4.5 stringr_1.6.0
## [10] pkgconfig_2.0.3 shape_1.4.6.1 fastmap_1.2.0
## [13] backports_1.5.0 inum_1.0-5 labeling_0.4.3
## [16] rmarkdown_2.31 prodlim_2026.03.11 ggfittext_0.10.4
## [19] nloptr_2.2.1 MatrixModels_0.5-4 purrr_1.2.1
## [22] xfun_0.57 glmnet_5.0 jomo_2.7-6
## [25] cachem_1.1.0 shades_1.5.0 jsonlite_2.0.0
## [28] recipes_1.3.2 gggenes_0.7.0 pan_1.9
## [31] broom_1.0.12 parallel_4.5.3 R6_2.6.1
## [34] bslib_0.10.0 stringi_1.8.7 RColorBrewer_1.1-3
## [37] parallelly_1.47.0 boot_1.3-32 rpart_4.1.24
## [40] lubridate_1.9.5 jquerylib_0.1.4 xgboost_3.2.1.1
## [43] Rcpp_1.1.1 iterators_1.0.14 knitr_1.51
## [46] future.apply_1.20.2 Matrix_1.7-4 splines_4.5.3
## [49] nnet_7.3-20 timechange_0.4.0 tidyselect_1.2.1
## [52] rstudioapi_0.18.0 yaml_2.3.12 partykit_1.2-27
## [55] timeDate_4052.112 codetools_0.2-20 listenv_0.10.1
## [58] tibble_3.3.1 plyr_1.8.9 withr_3.0.3
## [61] S7_0.2.2 evaluate_1.0.5 future_1.70.0
## [64] survival_3.8-6 pillar_1.11.1 foreach_1.5.2
## [67] stats4_4.5.3 reformulas_0.4.4 generics_0.1.4
## [70] scales_1.4.0 minqa_1.2.8 plotmo_3.7.0
## [73] globals_0.19.1 class_7.3-23 glue_1.8.0
## [76] tools_4.5.3 interp_1.1-6 data.table_1.18.4
## [79] lme4_2.0-1 ModelMetrics_1.2.2.2 gower_1.0.2
## [82] mvtnorm_1.3-6 grid_4.5.3 plotrix_3.8-14
## [85] tidyr_1.3.2 rbibutils_2.4.1 libcoin_1.0-12
## [88] ipred_0.9-15 nlme_3.1-168 earth_5.3.6
## [91] Formula_1.2-5 cli_3.6.6 viridisLite_0.4.3
## [94] lava_1.9.0 dplyr_1.2.0 gtable_0.3.6
## [97] sass_0.4.10 digest_0.6.39 farver_2.1.2
## [100] htmltools_0.5.9 lifecycle_1.0.5 hardhat_1.4.3
## [103] mitml_0.4-5 MASS_7.3-65
References
Lundberg, S. M., & Lee, S.-I. (2017). A unified approach to
interpreting model predictions.
Advances in Neural Information
Processing Systems,
30. Retrieved from
https://proceedings.neurips.cc/paper_files/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf
Molnar, C. (2025).
Interpretable machine learning: A guide for
making black box models explainable (3rd ed.). Retrieved from
https://christophm.github.io/interpretable-ml-book/
Spadaccini, G., Fokkema, M., & Wiel, M. A. van de. (2025). Discovery
and inference beyond linearity for epidemiological data by integrating
bayesian regression, tree ensembles and shapley values.
arXiv
Preprint arXiv:2505.00571.
https://doi.org/10.48550/arXiv.2505.00571