Package 'fastshap'

Fast approximate Shapley values

Description

Compute fast (approximate) Shapley values for a set of features using the Monte Carlo algorithm described in Strumbelj and Igor (2014). An efficient algorithm for tree-based models, commonly referred to as Tree SHAP, is also supported for lightgbm and xgboost models; see Lundberg et. al. (2020) for details.

Usage

explain(object, ...)

## Default S3 method:
explain(
  object,
  feature_names = NULL,
  X = NULL,
  nsim = 1,
  pred_wrapper = NULL,
  newdata = NULL,
  adjust = FALSE,
  baseline = NULL,
  shap_only = TRUE,
  parallel = FALSE,
  raw = FALSE,
  seed = NULL,
  exact = FALSE,
  batch_size = NULL,
  method = c("strumbelj", "permutation"),
  ...
)

## S3 method for class 'lm'
explain(
  object,
  feature_names = NULL,
  X,
  nsim = 1,
  pred_wrapper,
  newdata = NULL,
  adjust = FALSE,
  exact = FALSE,
  baseline = NULL,
  shap_only = TRUE,
  parallel = FALSE,
  ...
)

## S3 method for class 'xgb.Booster'
explain(
  object,
  feature_names = NULL,
  X = NULL,
  nsim = 1,
  pred_wrapper,
  newdata = NULL,
  adjust = FALSE,
  exact = FALSE,
  baseline = NULL,
  shap_only = TRUE,
  parallel = FALSE,
  ...
)

## S3 method for class 'lgb.Booster'
explain(
  object,
  feature_names = NULL,
  X = NULL,
  nsim = 1,
  pred_wrapper,
  newdata = NULL,
  adjust = FALSE,
  exact = FALSE,
  baseline = NULL,
  shap_only = TRUE,
  parallel = FALSE,
  ...
)

Arguments

object	A fitted model object (e.g., a ranger::ranger(), xgboost::xgboost(), or earth::earth() object, to name a few).
...	Additional optional arguments to be passed on to foreach::foreach() whenever parallel = TRUE (e.g., .packages for loading packages on the workers); ignored otherwise.
feature_names	Character string giving the names of the predictor variables (i.e., features) of interest. If NULL (default) they will be taken from the column names of X.
X	A matrix-like R object (e.g., a data frame or matrix) containing ONLY the feature columns from the training data (or suitable background data set). NOTE: This argument is required whenever exact = FALSE.
nsim	The number of Monte Carlo repetitions to use for estimating each Shapley value (only used when exact = FALSE). Default is 1. NOTE: To obtain the most accurate results, nsim should be set as large as feasibly possible.
pred_wrapper	Prediction function that requires two arguments, object and newdata. NOTE: This argument is required whenever exact = FALSE. The output of this function should be determined according to: Regression A numeric vector of predicted outcomes. Binary classification A vector of predicted class probabilities for the reference class. Multiclass classification A vector of predicted class probabilities for the reference class.
newdata	A matrix-like R object (e.g., a data frame or matrix) containing ONLY the feature columns for the observation(s) of interest; that is, the observation(s) you want to compute explanations for. Default is NULL which will produce approximate Shapley values for all the rows in X (i.e., the training data).
adjust	Logical indicating whether or not to adjust the sum of the estimated Shapley values to satisfy the local accuracy property; that is, to equal the difference between the model's prediction for that sample and the average prediction over all the training data (i.e., X). Default is FALSE and setting to TRUE requires nsim > 1.
baseline	Numeric baseline to use when adjusting the computed Shapley values to achieve local accuracy. Adjusted Shapley values for a single prediction (fx) will sum to the difference fx - baseline. Defaults to NULL, which corresponds to the average predictions computed from X, and zero otherwise (i.e., no additional predictions will be computed and the baseline attribute of the output will be set to zero).
shap_only	Logical indicating whether or not to include additional output useful for plotting (i.e., newdata and the baseline value.). This is convenient, for example, when using shapviz::shapviz() for plotting. Default is TRUE.
parallel	Logical indicating whether or not to compute the approximate Shapley values in parallel across features; default is FALSE. NOTE: setting parallel = TRUE requires the foreach package (in Suggests) and an appropriate (i.e., system-specific) parallel backend (e.g., via the doParallel package). For reproducible parallel results, use the doRNG package to register the backend. Note that the internal batching of predictions (see batch_size) means models whose predict() methods are already multithreaded (e.g., ranger or xgboost) often gain little from parallel = TRUE. The fastshap namespace is always added to .packages automatically (socket-based clusters, the default on Windows, start with no packages attached), so any additional packages needed by pred_wrapper should be supplied via .packages as usual. If fastshap (or a dependency) is installed somewhere other than the default library path (e.g., a development checkout, renv/pak project library, or during ⁠R CMD check⁠, which installs the package under test into an ephemeral library), fresh workers won't see it either; propagate the current session's library path to the cluster before registering it, e.g. parallel::clusterCall(cl, function(lp) .libPaths(lp), .libPaths()).
raw	Logical indicating whether or not to return the raw per-simulation Shapley values from each Monte Carlo replication. If TRUE, a 3-D array of dimensions ⁠n x p x nsim⁠ is returned, where n is the number of observations, p is the number of features, and nsim is the number of Monte Carlo replications; for example, apply(result, 1:2, sd) computes standard errors for each (observation, feature) pair. Only supported when adjust = FALSE. Default is FALSE.
seed	Integer specifying a random seed for reproducibility; passed to base::set.seed(). Default is NULL (no seed). NOTE: the Monte Carlo loop was restructured (vectorized) in version 0.2.0, so seeded results differ from those produced by fastshap (<= 0.1.5).
exact	Logical indicating whether to compute exact Shapley values. Currently only supported for stats::lm(), xgboost::xgboost(), and lightgbm::lightgbm() objects (binary/regression only — multiclass is not yet supported). Passing exact = TRUE for any other model type issues a warning and falls back to the Monte Carlo approximation. Note that exact = TRUE for stats::lm() returns explanations for each of the stats::terms() in the model. Default is FALSE.
batch_size	Optional positive integer giving the maximum number of rows to pass to pred_wrapper() per call. By default (NULL), all nsim Monte Carlo replications for a feature are stacked and evaluated in a single pair of prediction calls, which is fastest but requires 2 * nrow(newdata) * nsim * ncol(X) values of working memory per feature; set batch_size to bound the size of each prediction call instead. All of the randomness is drawn up front, so for a fixed seed the results are identical regardless of batch_size.
method	Character string specifying the Monte Carlo estimator to use (only relevant when exact = FALSE). "strumbelj" (the default) is the Strumbelj and Kononenko (2014) estimator described above: independent random coalitions per feature, 2 * nsim predictions per feature (2 * nsim * p total for p features), and supports feature_names subsetting (requesting fewer features costs proportionally less) and parallel = TRUE. "permutation" walks a single random permutation through all features per replication, flipping one feature at a time from background to foreground and reading off each feature's contribution as the consecutive prediction difference; this costs about nsim * p predictions for all features combined (roughly half of "strumbelj"'s cost for computing every feature), and each replication's contributions sum exactly to f(x) - f(background) (verifiable via raw = TRUE), so adjust = TRUE is inapplicable (and an error) for this method. Because the telescoping property requires walking every feature regardless of what is requested, feature_names only subsets the output under "permutation" — it does not reduce the amount of work performed, unlike "strumbelj". parallel = TRUE is not yet supported for "permutation".

Value

If shap_only = TRUE (the default), a matrix is returned with one column for each feature specified in feature_names (if feature_names = NULL, the default, there will be one column for each feature in X) and one row for each observation in newdata (if newdata = NULL, the default, there will be one row for each observation in X). Additionally, the returned matrix will have an attribute called "baseline" containing the baseline value. If shap_only = FALSE, then a list is returned with three components:

• shapley_values - a matrix of Shapley values (as described above);

• feature_values - the corresponding feature values (for plotting with shapviz::shapviz());

• baseline - the corresponding baseline value (for plotting with shapviz::shapviz()).

Note

Setting exact = TRUE with a linear model (i.e., an stats::lm() or stats::glm() object) assumes that the input features are independent. Also, setting adjust = TRUE is experimental and we follow the same approach as in shap.

References

Strumbelj, E., and Igor K. (2014). Explaining prediction models and individual predictions with feature contributions. Knowledge and information systems, 41(3), 647-665.

Lundberg, S. M., Erion, G., Chen, H., DeGrave, A., Prutkin, J. M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., and Lee, Su-In (2020). From local explanations to global understanding with explainable AI for trees. Nature Machine Intelligence, 2(1), 2522–5839.

See Also

You can find more examples (with larger and more realistic data sets) on the fastshap GitHub repository: https://github.com/bgreenwell/fastshap.

Examples

#
# A projection pursuit regression (PPR) example
#

# Load the sample data; see ?datasets::mtcars for details
data(mtcars)

# Fit a projection pursuit regression model
fit <- ppr(mpg ~ ., data = mtcars, nterms = 5)

# Prediction wrapper
pfun <- function(object, newdata) {  # needs to return a numeric vector
  predict(object, newdata = newdata)
}

# Compute approximate Shapley values using 10 Monte Carlo simulations
set.seed(101)  # for reproducibility
shap <- explain(fit, X = subset(mtcars, select = -mpg), nsim = 10,
                pred_wrapper = pfun)
head(shap)

---

Main, total, and interaction effect decomposition

Description

For each feature, decompose its contribution to a prediction into a main effect (the feature's effect in isolation), a total effect (the feature's effect together with every interaction it participates in), and the gap between them (an interaction diagnostic). Unlike explain(), which estimates Shapley values by averaging over randomly sampled coalitions, explain_effects() uses only the two coalitions with the highest Shapley weight — the empty coalition and the full coalition — so by default (nsim = NULL) it is deterministic given the background data, with no Monte Carlo simulation at all.

Usage

explain_effects(
  object,
  feature_names = NULL,
  X = NULL,
  pred_wrapper = NULL,
  newdata = NULL,
  nsim = NULL,
  batch_size = NULL,
  seed = NULL,
  ...
)

## S3 method for class 'explain_effects'
print(x, ...)

Arguments

object	A fitted model object (e.g., a ranger::ranger(), xgboost::xgboost(), or earth::earth() object, to name a few).
feature_names	Character string giving the names of the predictor variables (i.e., features) of interest. If NULL (default) they will be taken from the column names of X.
X	A matrix-like R object (e.g., a data frame or matrix) containing ONLY the feature columns from the training data (or suitable background data set). NOTE: This argument is required whenever exact = FALSE.
pred_wrapper	Prediction function that requires two arguments, object and newdata. NOTE: This argument is required whenever exact = FALSE. The output of this function should be determined according to: Regression A numeric vector of predicted outcomes. Binary classification A vector of predicted class probabilities for the reference class. Multiclass classification A vector of predicted class probabilities for the reference class.
newdata	A matrix-like R object (e.g., a data frame or matrix) containing ONLY the feature columns for the observation(s) of interest; that is, the observation(s) you want to compute explanations for. Default is NULL which will produce approximate Shapley values for all the rows in X (i.e., the training data).
nsim	Optional positive integer giving the number of background rows to sample per observation. Default (NULL) averages over every row of X (deterministic, exact given the background data, but requires 2 * nrow(newdata) * length(feature_names) * nrow(X) predictions); set nsim to average over a random sample of that many background rows instead (requires seed for reproducibility).
batch_size	Optional positive integer giving the maximum number of rows to pass to pred_wrapper() per call. By default (NULL), all nsim Monte Carlo replications for a feature are stacked and evaluated in a single pair of prediction calls, which is fastest but requires 2 * nrow(newdata) * nsim * ncol(X) values of working memory per feature; set batch_size to bound the size of each prediction call instead. All of the randomness is drawn up front, so for a fixed seed the results are identical regardless of batch_size.
seed	Integer specifying a random seed for reproducibility; passed to base::set.seed(). Default is NULL (no seed). NOTE: the Monte Carlo loop was restructured (vectorized) in version 0.2.0, so seeded results differ from those produced by fastshap (<= 0.1.5).
...	Additional optional arguments to be passed on to foreach::foreach() whenever parallel = TRUE (e.g., .packages for loading packages on the workers); ignored otherwise.
x	An object of class "explain_effects".

Details

For feature $j$ and observation $x$, averaging over background draws $w$:

main

⁠E_w[f(x_j, w_{-j})] - E_w[f(w)]⁠. Since no other feature of x is present, this is free of interaction effects — it is exactly the centered partial dependence of feature $j$ evaluated at $x_j$.

total

⁠f(x) - E_w[f(w_j, x_{-j})]⁠. This captures the main effect of feature $j$ together with every interaction it participates in, in the context of $x$.

interaction

total - main. Zero exactly when feature $j$ enters the model additively (no interactions involving $j$) at $x$; for a two-way interaction the gap splits equally between the two participating features.

shapley_values

(main + total) / 2. This equals the exact (background-marginal) Shapley value for any model with at most pairwise interactions; for higher-order interactions it is a cheap but biased approximation.

As with explain(), this uses marginal (interventional) sampling of the background, so the usual caveat applies: with correlated features, the constructed hybrid observations can be unrealistic combinations of feature values.

Value

An object of class "explain_effects": a list with components main, total, and interaction (each an nrow(newdata) x length(feature_names) matrix), shapley_values (same shape), and baseline (the average prediction over X).

References

Strumbelj, E., and Igor K. (2014). Explaining prediction models and individual predictions with feature contributions. Knowledge and information systems, 41(3), 647-665.

Examples

data(mtcars)
fit <- ppr(mpg ~ ., data = mtcars, nterms = 5)
pfun <- function(object, newdata) predict(object, newdata = newdata)
eff <- explain_effects(fit, X = subset(mtcars, select = -mpg),
                       pred_wrapper = pfun)
eff

---

Friedman benchmark data

Description

Simulate data from the Friedman 1 benchmark problem. These data were originally described in Friedman (1991) and Breiman (1996). For details, see sklearn.datasets.make_friedman1.

Usage

gen_friedman(n_samples = 100, n_features = 10, sigma = 0.1, seed = NULL)

Arguments

n_samples	Integer specifying the number of samples (i.e., rows) to generate. Default is 100.
n_features	Integer specifying the number of features to generate. Default is 10.
sigma	Numeric specifying the standard deviation of the noise.
seed	Integer specifying the random seed. If NULL (the default) the results will be different each time the function is run.

Note

This function is mostly used for internal testing.

References

Breiman, Leo (1996) Bagging predictors. Machine Learning 24, pages 123-140.

Friedman, Jerome H. (1991) Multivariate adaptive regression splines. The Annals of Statistics 19 (1), pages 1-67.

Examples

gen_friedman()

---

Survival of Titanic passengers

Description

A data set containing the survival outcome, passenger class, age, sex, and the number of family members for a large number of passengers aboard the ill-fated Titanic.

Usage

titanic

Format

A data frame with 1309 observations on the following 6 variables:

survived

binary with levels "yes" for survived and "no" otherwise;

pclass

integer giving the corresponding passenger (i.e., ticket) class with values 1–3;

age

the age in years of the corresponding passenger (with 263 missing values);

sex

factor giving the sex of each passenger with levels "male" and "female";

sibsp

integer giving the number of siblings/spouses aboard for each passenger (ranges from 0–8);

parch

integer giving the number of parents/children aboard for each passenger (ranges from 0–9).

Note

As mentioned in the column description, age contains 263 NAs (or missing values). For a complete version (or versions) of the data set, see titanic_mice.

Source

https://hbiostat.org/data/.

---

Survival of Titanic passengers

Description

The titanic data set contains 263 missing values (i.e., NA's) in the age column. This version of the data contains imputed values for the age column using multivariate imputation by chained equations via the mice package. Consequently, this is a list containing 11 imputed versions of the observations containd in the titanic data frame; each completed data sets has the same dimension and column structure as titanic.

Usage

titanic_mice

Format

An object of class mild (inherits from list) of length 21.

Source

Greenwell, Brandon M. (2022). Tree-Based Methods for Statistical Learning in R. CRC Press.

---