Functionality

Here is an overview of the main functions. You can read their documentation and see examples with ?function_name.

The shapr package implements Kernel SHAP estimation of dependence-aware Shapley values with eight different Monte Carlo-based approaches for estimating the conditional distributions of the data, namely "empirical", "gaussian", "copula", "ctree", "vaeac", "categorical", "timeseries", and "independence". shapr has also implemented two regression-based approaches "regression_separate" and "regression_surrogate". See Estimation approaches and plotting functionality below for examples. It is also possible to combine the different approaches, see the combined approach.

The package allows for parallelized computation through the futurepackage, see Parallelization for details.

The level of detail in the output can be controlled through the verbose argument. In addition, progress updates on the process of estimating the v(S)’s (and training the "vaeac" model) is available through the progressr package, supporting progress updates also for parallelized computation. See Verbosity and progress updates for details.

Moreover, the default behavior is to estimate the Shapley values iteratively, with increasing number of feature coalitions being added, and to stop estimation as the estimated Shapley values has achieved a certain level of stability. More information about this is provided in Iterative estimation The above, combined with batch computation of the v(S) values, enables fast and accurate estimation of the Shapley values in a memory friendly manner.

The package also provides functionality for computing Shapley values for groups of features, and custom function explanation, see Advanced usage. Finally, explanation of multiple output time series forecasting models are discussed in Explaining forecasting models.

Default behavior of explain

Below we provide brief descriptions of the most important parts of the default behavior of the explain function.

By default explain always compute feature-wise Shapley values. Groups of features can be explained by providing the feature groups through the group argument.

When there are five or less features (or feature groups), iterative estimation is by default disabled. The reason for this is that it is usually faster to estimate the Shapley values for all possible coalitions (v(S)), than to estimate the uncertainty of the Shapley values, and potentially stop estimation earlier. While iterative estimation is the default starting from six features, it is mainly when there are more than ten features, that it is most beneficial, and can save a lot of computation time. The reason for this is that the number of possible coalitions grows exponentially. These defaults can be overridden by setting the iterative argument to TRUE or FALSE. When using the iterative argument, the estimation for an observation is stopped when all Shapley value standard deviations are below t times the range of the Shapley values. The t value controls the convergence tolerance, defaults to 0.02, and can be set through the iterative_args$convergence_tol argument, see iterative estimation for more details.

Since the iterativeness default changes based on the number of features (or feature groups), the default is also to have no upper bound on the number of coalitions considered. This can be controlled through the max_n_coalitions argument.

The Kernel SHAP Method

Assume a predictive model f(x) for a response value y with features x ∈ ℝ^M, trained on a training set, and that we want to explain the predictions for new sets of data. This may be done using ideas from cooperative game theory, letting a single prediction take the place of the game being played and the features the place of the players. Letting N denote the set of all M players, and S ⊆ N be a subset of |S| players, the “contribution” function v(S) describes the total expected sum of payoffs the members of S can obtain by cooperation. The Shapley value (Shapley (1953)) is one way to distribute the total gains to the players, assuming that they all collaborate. The amount that player i gets is then

$$\phi_i(v) = \phi_i = \sum_{S \subseteq N \setminus\{i\}} \frac{|S| ! (M-| S| - 1)!}{M!}(v(S\cup \{i\})-v(S)),$$

that is, a weighted mean over all subsets S of players not containing player i. Lundberg and Lee (2017) define the contribution function for a certain subset S of these features x_S as v(S) = E[f(x)|x_S], the expected output of the predictive model conditional on the feature values of the subset. Lundberg and Lee (2017) names this type of Shapley values SHAP (SHapley Additive exPlanation) values. Since the conditional expectations can be written as

the conditional distributions p(x_S̄|x_S = x_S^*) are needed to compute the contributions. The Kernel SHAP method of Lundberg and Lee (2017) assumes feature independence, so that p(x_S̄|x_S = x_S^*) = p(x_S̄). If samples x_S̄^k, k = 1, …, K, from p(x_S̄|x_S = x_S^*) are available, the conditional expectation above can be approximated by

In Kernel SHAP, x_S̄^k, k = 1, …, K are sampled from the S̄-part of the training data, independently of x_S. This is motivated by using the training set as the empirical distribution of x_S̄, and assuming that x_S̄ is independent of x_S = x_S^*. Due to the independence assumption, if the features in a given model are highly dependent, the Kernel SHAP method may give a completely wrong answer. This can be avoided by estimating the conditional distribution p(x_S̄|x_S = x_S^*) directly and generating samples from this distribution. With this small change, the contributions and Shapley values may then be approximated as in the ordinary Kernel SHAP framework. Aas, Jullum, and Løland (2021) propose three different approaches for estimating the conditional probabilities which are implemented: empirical, gaussian and copula. The package also implements the ctree method from Redelmeier, Jullum, and Aas (2020), the vaeac method from Olsen et al. (2022) and a categorical for cateogorical data. The original independence approach of Lundberg and Lee (2017) is also available. The methods may also be combined, such that e.g. one method is used when conditioning on a small number of features, while another method is used otherwise. The shapr package also supports directly estimating the contribution function using regression. We briefly introduce the regression-based methods below, but we refer to the separate regression vignette (Shapley value explanations using the regression paradigm) and Olsen et al. (2024) for an in-depth explanation of the regression paradigm.

Multivariate Gaussian Distribution Approach

The first approach arises from the assumption that the feature vector x stems from a multivariate Gaussian distribution with some mean vector μ and covariance matrix Σ. Under this assumption, the conditional distribution $p(\boldsymbol{x}_{\bar{\mathcal{S}}} |\boldsymbol{x}_{\mathcal{S}}=\boldsymbol{x}_{\mathcal{S}}^*)$ is also multivariate Gaussian
$\text{N}_{|\bar{\mathcal{S}}|}(\boldsymbol{\mu}_{\bar{\mathcal{S}}|\mathcal{S}},\boldsymbol{\Sigma}_{\bar{\mathcal{S}}|\mathcal{S}})$, with analytical expressions for the conditional mean vector $\boldsymbol{\mu}_{\bar{\mathcal{S}}|\mathcal{S}}$ and covariance matrix $\boldsymbol{\Sigma}_{\bar{\mathcal{S}}|\mathcal{S}}$, see Aas, Jullum, and Løland (2021) for details. Hence, instead of sampling from the marginal empirical distribution of $\boldsymbol{x}_{\bar{\mathcal{S}}}$ approximated by the training data, we can sample from the Gaussian conditional distribution, which is fitted using the training data. Using the resulting samples $\boldsymbol{x}_{\bar{\mathcal{S}}}^k, k=1,\ldots,K$, the conditional expectations be approximated as in the Kernel SHAP.

Gaussian Copula Approach

If the features are far from multivariate Gaussian, an alternative approach is to instead represent the marginals by their empirical distributions, and model the dependence structure by a Gaussian copula. Assuming a Gaussian copula, we may convert the marginals of the training data to Gaussian features using their empirical distributions, and then fit a multivariate Gaussian distribution to these.

To produce samples from the conditional distribution $p(\boldsymbol{x}_{\bar{\mathcal{S}}} |\boldsymbol{x}_{\mathcal{S}}=\boldsymbol{x}_{\mathcal{S}}^*)$, we convert the marginals of x_𝒮 to Gaussians, sample from the conditional Gaussian distribution as above, and convert the marginals of the samples back to the original distribution. Those samples are then used to approximate the sample from the resulting multivariate Gaussian conditional distribution. While other copulas may be used, the Gaussian copula has the benefit that we may use the analytical expressions for the conditionals $\boldsymbol{\mu}_{\bar{\mathcal{S}}|\mathcal{S}}$ and $\boldsymbol{\Sigma}_{\bar{\mathcal{S}}|\mathcal{S}}$. Finally, we may convert the marginals back to their original distribution, and use the resulting samples to approximate the conditional expectations as in Kernel SHAP.

Empirical Conditional Distribution Approach

If both the dependence structure and the marginal distributions of x are very far from the Gaussian, neither of the two aforementioned methods will work very well. Few methods exists for the non-parametric estimation of conditional densities, and the classic kernel estimator (Rosenblatt (1956)) for non-parametric density estimation suffers greatly from the curse of dimensionality and does not provide a way to generate samples from the estimated distribution. For such situations, Aas, Jullum, and Løland (2021) propose an empirical conditional approach to sample approximately from $p(\boldsymbol{x}_{\bar{\mathcal{S}}}|\boldsymbol{x}_{\mathcal{S}}^*)$. The idea is to compute weights w_𝒮(x^*, xⁱ), i = 1, ..., n_train for all training instances based on their Mahalanobis distances (in the S subset only) to the instance x^* to be explained. Instead of sampling from this weighted (conditional) empirical distribution, Aas, Jullum, and Løland (2021) suggests a more efficient variant, using only the K instances with the largest weights:

$$v_{\text{condKerSHAP}}(\mathcal{S}) = \frac{\sum_{k=1}^K w_{\mathcal{S}}(\boldsymbol{x}^*, \boldsymbol{x}^{[k]}) f(\boldsymbol{x}_{\bar{\mathcal{S}}}^{[k]}, \boldsymbol{x}_{\mathcal{S}}^*)}{\sum_{k=1}^K w_{\mathcal{S}}(\boldsymbol{x}^*,\boldsymbol{x}^{[k]})},$$

The number of samples K to be used in the approximate prediction can for instance be chosen such that the K largest weights accounts for a fraction η, for example 0.9, of the total weight. If K exceeds a certain limit, for instance 5, 000, it might be set to that limit. A bandwidth parameter σ used to scale the weights, must also be specified. This choice may be viewed as a bias-variance trade-off. A small σ puts most of the weight to a few of the closest training observations and thereby gives low bias, but high variance. When σ → ∞, this method converges to the original Kernel SHAP assuming feature independence. Typically, when the features are highly dependent, a small σ is typically needed such that the bias does not dominate. Aas, Jullum, and Løland (2021) show that a proper criterion for selecting σ is a small-sample-size corrected version of the AIC known as AICc. As calculation of it is computationally intensive, an approximate version of the selection criterion is also suggested. Details on this is found in Aas, Jullum, and Løland (2021).

Conditional Inference Tree Approach

The previous three methods can only handle numerical data. This means that if the data contains categorical/discrete/ordinal features, the features first have to be one-hot encoded. When the number of levels/features is large, this is not feasible. An approach that handles mixed (i.e numerical, categorical, discrete, ordinal) features and both univariate and multivariate responses is conditional inference trees (Hothorn, Hornik, and Zeileis (2006)).

Conditional inference trees is a special tree fitting procedure that relies on hypothesis tests to choose both the splitting feature and the splitting point. The tree fitting procedure is sequential: first a splitting feature is chosen (the feature that is least independent of the response), and then a splitting point is chosen for this feature. This decreases the chance of being biased towards features with many splits (Hothorn, Hornik, and Zeileis (2006)).

We use conditional inference trees (ctree) to model the conditional distribution, $p(\boldsymbol{x}_{\bar{\mathcal{S}}}|\boldsymbol{x}_{\mathcal{S}}^*)$, found in the Shapley methodology. First, we fit a different conditional inference tree to each conditional distribution. Once a tree is fit for given dependent features, the end node of x_𝒮^* is found. Then, we sample from this end node and use the resulting samples, $\boldsymbol{x}_{\bar{\mathcal{S}}}^k, k=1,\ldots,K$, when approximating the conditional expectations as in Kernel SHAP. See Redelmeier, Jullum, and Aas (2020) for more details.

The conditional inference trees are fit using the party or partykit packages (Hothorn and Zeileis (2015)).

Variational AutoEncoder with Arbitrary Conditioning (vaeac) Approach

Another approach that supports mixed features is the Variational AutoEncoder with Arbitrary Conditioning (Olsen et al. (2022)), abbreviated to vaeac. The vaeac is an extension of the regular variational autoencoder (Kingma and Welling (2014)), but instead of giving a probabilistic representation of the distribution p(x) it gives a probabilistic representation of the conditional distribution $p(\boldsymbol{x}_{\bar{\mathcal{S}}} \mid \boldsymbol{x}_{\mathcal{S}})$, for all possible feature subsets 𝒮 ⊆ ℳ simultaneously, where ℳ is the set of all features. That is, only a single vaeac model is needed to model all conditional distributions.

The vaeac consists of three neural networks: a full encoder, a masked encoder, and a decoder. The encoders map the full and masked/conditional input representations, i.e., x and x_𝒮, respectively, to latent probabilistic representations. Sampled instances from this latent probabilistic representations are sent to the decoder, which maps them back to the feature space and provides a samplable probabilistic representation for the unconditioned features $\boldsymbol{x}_{\bar{\mathcal{S}}}$. The full encoder is only used during the training phase of the vaeac model to guide the training process of the masked encoder, as the former relies on the full input sample x, which is not accessible in the deployment phase (when we generate the Monte Carlo samples), as we only have access to x_𝒮. The networks are trained by minimizing a variational lower bound. See Section 3 in Olsen et al. (2022) for an in-depth introduction to the vaeac methodology. We use the vaeac model at the epoch which obtains the lowest validation IWAE score to generate the Monte Carlo samples used in the Shapley value computations.

We fit the vaeac model using the R-package torch(Falbel and Luraschi (2023)). The main parameters are the the number of layers in the networks (vaeac.depth), the width of the layers (vaeac.width), the number of dimensions in the latent space (vaeac.latent_dim), the activation function between the layers in the networks (vaeac.activation_function), the learning rate in the ADAM optimizer (vaeac.lr), the number of vaeac models to initiate to remedy poorly initiated model parameter values (vaeac.n_vaeacs_initialize), and the number of learning epochs (vaeac.epochs). See ?shapr::setup_approach.vaeac for a more detailed description of the parameters.

There are additional extra parameters which can be set by including a named list in the call to the explain() function. For example, we can the change the batch size to 32 by including vaeac.extra_parameters = list(vaeac.batch_size = 32) as an argument in the explain() function. See ?shapr::vaeac_get_extra_para_default for a description of the possible extra parameters to the vaeac approach. Note that the main parameters are entered as arguments directly in explain() function, while the extra parameters are specified through a named list called vaeac.extra_parameters.

Categorical Approach

When the features are all categorical, we can estimate the conditional expectations using basic statistical formulas. For example, if we have three features, x₁, x₂, x₃ with three levels each (indicated as 1, 2, 3), and we are provided with a table of counts indicating how many times each combination of feature values occurs, we can estimate the marginal and conditional probabilities as follows. Marginal probabilities are estimated by dividing the number of times a given feature (or features) takes on a certain value in the data set with the total number of observations in the data set. Condititional probabilities (for example, P(X₁ = 1|X₂ = 1)) are estimated by first subsetting the data set to reflect the conditioning (i.e., extracting all rows where X₂ = 1), and then dividing the number of times the feature on the left hand side of | takes the given value in this subset by the total number of observations in this subset. Once the marginal and conditional probabilities are estimated for all combinations of feature values, each conditional expectation can be calculated. For example, the expected value of X₁ given X₂ = 1 and X₃ = 2 is $$E(X_1|X_2, X_3) = \sum_{x}x P(X_1 = x | X_2=1, X_3=2) = \sum_{x} x \frac{P(X_1 = x, X_2 = 1, X_3 = 2)}{P(X_2=1, X_3=2)}.$$.

Separate and Surrogate Regression Approaches

Another paradigm for estimating the contribution function is the regression paradigm. In contrast to the methods above, which belong to the Monte Carlo paradigm, the regression based methods use regression models to estimate the contribution function v(S) = E[f(x)|x_S = x_S^*] directly. The separate regression method class fits a separate regression model for each coalition S, while the surrogate regression method class fits a single regression model to simultaneously predict the contribution function for all coalitions. We refer to Olsen et al. (2024) for when one should use the different paradigms, method classes, and methods.

In a separate vignette (Shapley value explanations using the regression paradigm), we elaborate and demonstrate the regression paradigm. We describe how to specify the regression model, enable automatic cross-validation of the model’s hyperparameters, and apply pre-processing steps to the data before fitting the regression models. Olsen et al. (2024) divides the regression paradigm into the separate and surrogate regression method classes. In the separate vignette, we briefly introduce the two method classes. For an in-depth explanation, we refer the reader to Sections 3.5 and 3.6 in Olsen et al. (2024).

MSEv evaluation criterion

We can use the MSE_v criterion proposed by Frye et al. (2021), and later used by, e.g., Olsen et al. (2022) and Olsen et al. (2024), to evaluate and rank the approaches/methods. The MSE_v is given by

where v̂_q is the estimated contribution function using method q and N_𝒮 = |𝒫^*(ℳ)| = 2^M − 2, i.e., we have removed the empty (𝒮 = ∅) and the grand combinations (𝒮 = ℳ) as they are method independent. Meaning that these two combinations do not influence the ranking of the methods as the methods are not used to compute the contribution function for them.

The motivation behind the MSE_v criterion is that 𝔼_𝒮𝔼_x(v_true(𝒮, x) − v̂_q(𝒮, x))² can be decomposed as

see Appendix A in Covert, Lundberg, and Lee (2020). The first term on the right-hand side of the equation above can be estimated by MSE_v, while the second term is a fixed (unknown) constant not influenced by the approach . Thus, a low value of MSE_v indicates that the estimated contribution function v̂_q is closer to the true counterpart v_true than a high value.

In shapr, we allow for weighting the combinations in the MSE_v evaluation criterion either uniformly or by using the corresponding Shapley kernel weights (or the sampling frequencies when sampling of combinations is used). This is determined by the logical parameter MSEv_uniform_comb_weights in the explain() function, and the default is to do uniform weighting, that is, MSEv_uniform_comb_weights = TRUE.

# We use more explicands here for more stable confidence intervals
ind_x_explain_many <- 1:25
x_train <- data[-ind_x_explain_many, ..x_var]
y_train <- data[-ind_x_explain_many, get(y_var)]
x_explain <- data[ind_x_explain_many, ..x_var]

# Fitting a basic xgboost model to the training data
model <- xgboost::xgboost(
  data = as.matrix(x_train),
  label = y_train,
  nround = 20,
  verbose = FALSE
)

# Specifying the phi_0, i.e. the expected prediction without any features
p0 <- mean(y_train)

# Independence approach
explanation_independence <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "independence",
  phi0 = p0,
  n_MC_samples = 1e2,
  MSEv_uniform_comb_weights = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:35 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: independence
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 25
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9321f2ff20.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Empirical approach
explanation_empirical <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0,
  n_MC_samples = 1e2,
  MSEv_uniform_comb_weights = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:36 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 25
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9372002e18.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Gaussian 1e1 approach
explanation_gaussian_1e1 <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "gaussian",
  phi0 = p0,
  n_MC_samples = 1e1,
  MSEv_uniform_comb_weights = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:43 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 25
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9371245336.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Gaussian 1e2 approach
explanation_gaussian_1e2 <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "gaussian",
  phi0 = p0,
  n_MC_samples = 1e2,
  MSEv_uniform_comb_weights = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:43 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 25
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c93155a4322.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Combined approach
explanation_combined <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = c("gaussian", "empirical", "independence"),
  phi0 = p0,
  n_MC_samples = 1e2,
  MSEv_uniform_comb_weights = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:43 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian, empirical, and independence
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 25
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c937ba6bac.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Create a list of explanations with names
explanation_list_named <- list(
  "Ind." = explanation_independence,
  "Emp." = explanation_empirical,
  "Gaus. 1e1" = explanation_gaussian_1e1,
  "Gaus. 1e2" = explanation_gaussian_1e2,
  "Combined" = explanation_combined
)

# Create the MSEv plots with approximate 95% confidence intervals
MSEv_plots <- plot_MSEv_eval_crit(explanation_list_named,
  plot_type = c("overall", "comb", "explicand"),
  CI_level = 0.95
)

# 5 plots are made
names(MSEv_plots)
#> [1] "MSEv_explicand_bar"        "MSEv_explicand_line_point" "MSEv_coalition_bar"       
#> [4] "MSEv_coalition_line_point" "MSEv_bar"

# The main plot of the overall MSEv averaged over both the combinations and observations
MSEv_plots$MSEv_bar

# The MSEv averaged over only the explicands for each combinations
MSEv_plots$MSEv_combination_bar
#> NULL

# The MSEv averaged over only the combinations for each observation/explicand
MSEv_plots$MSEv_explicand_bar

# To see which coalition S each of the `id_combination` corresponds to,
# i.e., which features that are conditions on.
explanation_list_named[[1]]$MSEv$MSEv_combination[, c("id_combination", "features")]
#> NULL

# We can specify which test observations or combinations to plot
plot_MSEv_eval_crit(explanation_list_named,
  plot_type = "explicand",
  index_x_explain = c(1, 3:4, 6),
  CI_level = 0.95
)$MSEv_explicand_bar

plot_MSEv_eval_crit(explanation_list_named,
  plot_type = "comb",
  id_coalition = c(3, 4, 9, 13:15),
  CI_level = 0.95
)$MSEv_combination_bar
#> NULL

bar_text_n_decimals <- 1
plot_MSEv_eval_crit(explanation_list_named) +
  ggplot2::scale_x_discrete(limits = rev(levels(MSEv_plots$MSEv_bar$data$Method))) +
  ggplot2::coord_flip() +
  ggplot2::scale_fill_brewer(palette = "Paired") +
  ggplot2::theme_minimal() + # This must be set before other theme calls
  ggplot2::theme(
    plot.title = ggplot2::element_text(size = 10),
    legend.position = "bottom"
  ) +
  ggplot2::geom_text(
    ggplot2::aes(label = sprintf(
      paste("%.", sprintf("%d", bar_text_n_decimals), "f", sep = ""),
      round(MSEv, bar_text_n_decimals)
    )),
    vjust = -0.35, # This number might need altering for different plots sizes
    hjust = 1.1, # This number might need altering for different plots sizes
    color = "black",
    position = ggplot2::position_dodge(0.9),
    size = 4
  )

Batch computation

The computational complexity of Shapley value based explanations grows fast in the number of features, as the number of conditional expectations one needs to estimate in the Shapley formula grows exponentially. As outlined above, the estimating of each of these conditional expectations is also computationally expensive, typically requiring estimation of a conditional probability distribution, followed by Monte Carlo integration. These computations are not only heavy for the CPU, they also require a lot of memory (RAM), which typically is a limited resource. By doing the most resource hungry computations (the computation of v(S)) in sequential batches with different feature subsets S, the memory usage can be significantly reduces. The user can control the number of batches by setting the two arguments extra_computation_args$max_batch_size (defaults to 10) and extra_computation_args$min_n_batches (defaults to 10).

Parallelized computation

In addition to reducing the memory consumption, the batch computing allows the computations within each batch to be performed in parallel. The parallelization in shapr::explain() is handled by the future_apply package which builds on the future environment. These packages work on all OS, allows the user to decide the parallelization backend (mutliple R procesess or forking), works directly with hpc clusters, and also supports progress updates for the parallelized task via the associated progressr package (see Verbosity and progress updates).

Note that, since it takes some time to duplicate data into different processes/machines when running in parallel, it is not always preferrable to run shapr::explain() in parallel, at least not with many parallel sessions/workers. Parallelization also increases the memory consumption proportionally, so you may want to limit the number of workers for that reason too. Below is a basic example of a parallelization with two workers.

library(future)
future::plan(multisession, workers = 2)

explanation_par <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 32, 
#> and is therefore set to 2^n_features = 32.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:51 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 5
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9348628d8d.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 32 of 32 coalitions.

future::plan(sequential) # To return to non-parallel computation

Combined approach

In addition to letting the user select one of the five aforementioned approaches for estimating the conditional distribution of the data (i.e. approach equals either "gaussian", "copula", "empirical", "ctree", "vaeac", "categorical") or "timeseries", the package allows the user to combine the given approaches. The 'regression_surrogate' and 'regression_separate approaches are not supported for the combined approach. To simplify the usage, the flexibility is restricted such that the same approach is used when conditioning on the same number of features. This is also in line Aas, Jullum, and Løland (2021, sec. 3.4).

This can be done by setting approach equal to a character vector, where the length of the vector is one less than the number of features in the model. Consider a situation where you have trained a model that consists of 10 features, and you would like to use the "empirical" approach when you condition on 1-3 features, the "copula" approach when you condition on 4-5 features, and the "gaussian" approach when conditioning on 6 or more features. This can be applied by simply passing approach = c(rep("empirical", 3), rep("copula", 2), rep("gaussian", 4)), i.e. approach[i] determines which method to use when conditioning on i features. Conditioning on all features needs no approach as that is given by the complete prediction itself, and should thus not be part of the vector.

The code below exemplifies this approach for a case where there are four features, using "empirical", "copula" and "gaussian" when conditioning on respectively 1, 2 and 3 features.

library(xgboost)
library(data.table)

data("airquality")
data <- data.table::as.data.table(airquality)
data <- data[complete.cases(data), ]

x_var <- c("Solar.R", "Wind", "Temp", "Month")
y_var <- "Ozone"

ind_x_explain <- 1:6
x_train <- data[-ind_x_explain, ..x_var]
y_train <- data[-ind_x_explain, get(y_var)]
x_explain <- data[ind_x_explain, ..x_var]

# Set seed for reproducibility
set.seed(123)

# Fitting a basic xgboost model to the training data
model <- xgboost::xgboost(
  data = as.matrix(x_train),
  label = y_train,
  nround = 20,
  verbose = FALSE
)

# Specifying the phi_0, i.e. the expected prediction without any features
p0 <- mean(y_train)


# Use the combined approach
explanation_combined <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = c("empirical", "copula", "gaussian"),
  phi0 = p0
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:54 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: empirical, copula, and gaussian
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c936a965bd5.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.
# Plot the resulting explanations for observations 1 and 6, excluding
# the no-covariate effect
plot(explanation_combined, bar_plot_phi0 = FALSE, index_x_explain = c(1, 6))

As a second example using "ctree" to condition on 1 and 2 features, and "empirical" when conditioning on 3 features:

# Use the combined approach
explanation_combined <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = c("ctree", "ctree", "empirical"),
  phi0 = p0
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:57 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: ctree, ctree, and empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c93207aeb18.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

Explain groups of features

In some cases, especially when the number of features is very large, it may be more appropriate to explain predictions in terms of groups of features instead of single features, see (Jullum, Redelmeier, and Aas (2021)) for intuition and real world examples. Explaining prediction in terms of groups of features is very easy using shapr:

# Define the feature groups
group_list <- list(
  A = c("Temp", "Month"),
  B = c("Wind", "Solar.R")
)

# Use the empirical approach
explanation_group <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0,
  group = group_list
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_groups = 4, 
#> and is therefore set to 2^n_groups = 4.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:44:59 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of group-wise Shapley values: 2
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c93241851b9.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 4 of 4 coalitions.
# Prints the group-wise explanations
explanation_group
#>    explain_id  none      A        B
#>         <int> <num>  <num>    <num>
#> 1:          1 43.09 -29.25  16.0731
#> 2:          2 43.09 -15.17  -7.8373
#> 3:          3 43.09 -13.07 -10.8778
#> 4:          4 43.09 -17.47   0.6653
#> 5:          5 43.09 -28.27   3.5289
#> 6:          6 43.09 -20.59  -3.3793
# Plots the group-wise explanations
plot(explanation_group, bar_plot_phi0 = TRUE, index_x_explain = c(1, 6))

Explain custom models

shapr currently natively supports explanation of predictions from models fitted with the following functions:

• stats::lm
• stats::glm
• ranger::ranger
• mgcv::gam
• xgboost::xgboost/xgboost::xgb.train
• workflows::workflow

Any continuous response regression model or binary classification model of these model classes, can be explained with the package directly as exemplified above, while we give an example for the workflows::workflow in the tidymodels/workflows section. Moreover, essentially any feature dependent prediction model can be explained by the package by specifying two (or one) simple additional functions for your model.

The first function is predict_model, taking the model and data (as a matrix or data.frame/data.table) as input and outputting the corresponding prediction as a numeric vector. The second (optional, but highly recommended) function is get_model_specs, taking the model as input and outputting a list with the following elements: labels (vector with the feature names to compute Shapley values for), classes (a named vector with the labels as names and the class type as elements), factor_levels (a named list with the labels as names and vectors with the factor levels as elements (NULL if the feature is not a factor)). The get_model_specs function is used to check that the format of the data passed to explain have the correct format in terms of the necessary feature columns being available and having the correct class/attributes. It is highly recommended to do such checks in order to ensure correct usage of explain. If, for some reason, such checking is not desirable, one does not have to provide the get_model_specs function. This will, however, throw a warning that all feature consistency checking against the model is disabled.

Once the above functions are created, you can explain predictions from this model as before by passing the functions through the input arguments predict_model and get_model_specs of explain().

These functions can be made general enough to handle all supported model types of that class, or they can be made minimal, possibly only allowing explanation of the specific version of the model class at hand. Below we give examples of both full support versions of these functions and a minimal version which skips the get_model_specs function. We do this for the gbm model class from the gbm package, fitted to the same airquality data set as used above.

library(gbm)

formula_gbm <- as.formula(paste0(y_var, "~", paste0(x_var, collapse = "+")))
# Fitting a gbm model
set.seed(825)
model_gbm <- gbm::gbm(
  formula_gbm,
  data = cbind(x_train, Ozone = y_train),
  distribution = "gaussian"
)

#### Full feature versions of the three required model functions ####
MY_predict_model <- function(x, newdata) {
  if (!requireNamespace("gbm", quietly = TRUE)) {
    stop("The gbm package is required for predicting train models")
  }
  model_type <- ifelse(
    x$distribution$name %in% c("bernoulli", "adaboost"),
    "classification",
    "regression"
  )
  if (model_type == "classification") {
    predict(x, as.data.frame(newdata), type = "response", n.trees = x$n.trees)
  } else {
    predict(x, as.data.frame(newdata), n.trees = x$n.trees)
  }
}
MY_get_model_specs <- function(x) {
  feature_specs <- list()
  feature_specs$labels <- labels(x$Terms)
  m <- length(feature_specs$labels)
  feature_specs$classes <- attr(x$Terms, "dataClasses")[-1]
  feature_specs$factor_levels <- setNames(vector("list", m), feature_specs$labels)
  feature_specs$factor_levels[feature_specs$classes == "factor"] <- NA # model object doesn't contain factor levels info
  return(feature_specs)
}

# Compute the Shapley values
set.seed(123)
p0 <- mean(y_train)
explanation_custom <- explain(
  model = model_gbm,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0,
  predict_model = MY_predict_model,
  get_model_specs = MY_get_model_specs
)
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:01 ────────────────────────────────────
#> • Model class: <gbm>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c933748f819.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Plot results
plot(explanation_custom, index_x_explain = c(1, 6))

#### Minimal version of the custom model setup ####
# Note: Working only for this exact version of the model class
# Avoiding to define get_model_specs skips all feature
# consistency checking between your data and model
MY_MINIMAL_predict_model <- function(x, newdata) {
  predict(x, as.data.frame(newdata), n.trees = x$n.trees)
}

# Compute the Shapley values
set.seed(123)
explanation_custom_minimal <- explain(
  model = model_gbm,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0,
  predict_model = MY_MINIMAL_predict_model
)
#> Note: You passed a model to explain() which is not natively supported, and did not supply a 'get_model_specs' function to explain().
#> Consistency checks between model and data is therefore disabled.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:07 ────────────────────────────────────
#> • Model class: <gbm>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9344b0ff8e.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Plot results
plot(explanation_custom_minimal, index_x_explain = c(1, 6))

Tidymodels and workflows

In this section, we demonstrate how to use shapr to explain tidymodels models fitted using workflows. In the example above, we directly used the xgboost package to fit the xgboost model. However, we can also fit the xgboost model using the tidymodels package. These fits will be identical as tidymodels calls xgboost internally. which we demonstrate in the example below. Note that we can replace xgboost (i.e., parsnip::boost_tree) with any other fitted tidymodels in the workflows procedure outlined below.

# Fitting a basic xgboost model to the training data using tidymodels
set.seed(123) # Set the same seed as above
all_var <- c(y_var, x_var)
train <- data[-ind_x_explain, ..all_var]

# Fitting the `tidymodels` model using `workflows`
model_tidymodels <- parsnip::fit(
  workflows::add_recipe(
    workflows::add_model(
      workflows::workflow(),
      parsnip::boost_tree(trees = 20, engine = "xgboost", mode = "regression")
    ),
    recipes::recipe(Ozone ~ ., data = train)
  ),
  data = train
)

# # We can also specify the same model using pipes `%>%` by (if pipes are installed/loaded)
# model_tidymodels <-
#   workflows::workflow() %>%
#   workflows::add_model(parsnip::boost_tree(trees = 20, engine = "xgboost", mode = "regression")) %>%
#   workflows::add_recipe(recipes::recipe(Ozone ~ ., data = train)) %>%
#   parsnip::fit(data = train)

# See that the output of the two models are identical
all.equal(predict(model_tidymodels, x_train)$.pred, predict(model, as.matrix(x_train)))
#> [1] TRUE

# Create the Shapley values for the tidymodels version
explanation_tidymodels <- explain(
  model = model_tidymodels,
  x_explain = x_explain,
  x_train = x_train,
  approach = "empirical",
  phi0 = p0
)
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:13 ────────────────────────────────────
#> • Model class: <workflow>
#> • Approach: empirical
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c93476395f2.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# See that the Shapley value explanations are identical too
all.equal(explanation$shapley_values_est, explanation_tidymodels$shapley_values_est)
#> [1] TRUE

The parameters of the vaeac approach

The vaeac approach is a very flexible method that supports mixed data. The main parameters are the the number of layers in the networks (vaeac.depth), the width of the layers (vaeac.width), the number of dimensions in the latent space (vaeac.latent_dim), the activation function between the layers in the networks (vaeac.activation_function), the learning rate in the ADAM optimizer (vaeac.lr), the number of vaeac models to initiate to remedy poorly initiated model parameter values (vaeac.n_vaeacs_initialize), and the number of learning epochs (vaeac.epochs). Call ?shapr::setup_approach.vaeac for a more detailed description of the parameters.

There are additional extra parameters which can be set by including a named list in the call to the explain() function. For example, we can the change the batch size to 32 by including vaeac.extra_parameters = list(vaeac.batch_size = 32) as a parameter in the call the explain() function. See ?shapr::vaeac_get_extra_para_default for a description of the possible extra parameters to the vaeac approach. The main parameters are directly entered to the explain() function, while the extra parameters are included in a named list called vaeac.extra_parameters.

explanation_vaeac <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "vaeac",
  phi0 = p0,
  n_MC_samples = 100,
  vaeac.width = 16,
  vaeac.depth = 2,
  vaeac.epochs = 3,
  vaeac.n_vaeacs_initialize = 2
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:19 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: vaeac
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c932349a0b1.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

Can look at the training and validation error for the trained vaeac model and see that vaeac.epochs = 3 is likely to few epochs as it still seems like the vaeac model is learning.

# Look at the training and validation errors.
vaeac_plot_eval_crit(list("Vaeac 3 epochs" = explanation_vaeac), plot_type = "method")

explanation_vaeac_early_stop <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "vaeac",
  phi0 = p0,
  n_MC_samples = 100,
  vaeac.width = 16,
  vaeac.depth = 2,
  vaeac.epochs = 1000, # Set it to a large number
  vaeac.n_vaeacs_initialize = 2,
  vaeac.extra_parameters = list(vaeac.epochs_early_stopping = 2)
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 16, 
#> and is therefore set to 2^n_features = 16.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:34 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: vaeac
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 4
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c9361ebec61.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 16 of 16 coalitions.

# Look at the training and validation errors.
vaeac_plot_eval_crit(
  list("Vaeac 3 epochs" = explanation_vaeac, "Vaeac early stopping" = explanation_vaeac_early_stop),
  plot_type = "method"
)

plot_MSEv_eval_crit(list("Vaeac 3 epochs" = explanation_vaeac, "Vaeac early stopping" = explanation_vaeac_early_stop))

Continued computation

In this section, we demonstrate how to continue to improve estimation accuracy with additional coalition samples, from a previous Shapley value computation based on shapr::explain() with the iterative estimation procedure. This can be done either by passing an existing object of class shapr, or by passing a string with the path to the intermediately saved results. The latter is found at SHAPR_OBJ$saving_path, defaults to a temporary folder, and is updated after each iteration. This can be particularly handy for long-running computations.

# First we run the computation with the iterative estimation procedure for a limited number of coalition samples
library(xgboost)
library(data.table)

data("airquality")
data <- data.table::as.data.table(airquality)
data <- data[complete.cases(data), ]

x_var <- c("Solar.R", "Wind", "Temp", "Month", "Day")
y_var <- "Ozone"

ind_x_explain <- 1:6
x_train <- data[-ind_x_explain, ..x_var]
y_train <- data[-ind_x_explain, get(y_var)]
x_explain <- data[ind_x_explain, ..x_var]

# Set seed for reproducibility
set.seed(123)

# Fitting a basic xgboost model to the training data
model <- xgboost::xgboost(
  data = as.matrix(x_train),
  label = y_train,
  nround = 20,
  verbose = FALSE
)

# Specifying the phi_0, i.e. the expected prediction without any features
p0 <- mean(y_train)

# Initial explanation computation
ex_init <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "gaussian",
  phi0 = p0,
  max_n_coalitions = 20,
  iterative = TRUE
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:45:57 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian
#> • Iterative estimation: TRUE
#> • Number of feature-wise Shapley values: 5
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c93667c79.rds'
#> 
#> ── iterative computation started ──
#> 
#> ── Iteration 1 ───────────────────────────────────────────────────────────────────────────
#> ℹ Using 5 of 32 coalitions, 5 new.
#> 
#> ── Iteration 2 ───────────────────────────────────────────────────────────────────────────
#> ℹ Using 10 of 32 coalitions, 4 new.
#> 
#> ── Iteration 3 ───────────────────────────────────────────────────────────────────────────
#> ℹ Using 12 of 32 coalitions, 2 new.
#> 
#> ── Iteration 4 ───────────────────────────────────────────────────────────────────────────
#> ℹ Using 16 of 32 coalitions, 4 new.
#> 
#> ── Iteration 5 ───────────────────────────────────────────────────────────────────────────
#> ℹ Using 18 of 32 coalitions, 2 new.

# Using the ex_init object to continue the computation with 5 more coalition samples
ex_further <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "gaussian",
  phi0 = p0,
  max_n_coalitions = 25,
  iterative_args = list(convergence_tol = 0.005), # Decrease the convergence threshold
  prev_shapr_object = ex_init
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:46:01 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 5
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c937cd2ac71.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 24 of 32 coalitions.

print(ex_further$saving_path)
#> [1] "/tmp/RtmpCcS4cy/shapr_obj_1b8c937cd2ac71.rds"

# Using the ex_init object to continue the computation for the remaining coalition samples
# but this time using the path to the saved intermediate estimation object
ex_even_further <- explain(
  model = model,
  x_explain = x_explain,
  x_train = x_train,
  approach = "gaussian",
  phi0 = p0,
  max_n_coalitions = NULL,
  prev_shapr_object = ex_further$saving_path
)
#> Note: Feature classes extracted from the model contains NA.
#> Assuming feature classes from the data are correct.
#> Success with message:
#> max_n_coalitions is NULL or larger than or 2^n_features = 32, 
#> and is therefore set to 2^n_features = 32.
#> 
#> ── Starting `shapr::explain()` at 2024-11-21 20:46:02 ────────────────────────────────────
#> • Model class: <xgb.Booster>
#> • Approach: gaussian
#> • Iterative estimation: FALSE
#> • Number of feature-wise Shapley values: 5
#> • Number of observations to explain: 6
#> • Computations (temporary) saved at: '/tmp/RtmpCcS4cy/shapr_obj_1b8c936aa37ec8.rds'
#> 
#> ── Main computation started ──
#> 
#> ℹ Using 32 of 32 coalitions.