Package {shapley}

Type: Package
Title: Weighted Mean SHAP and CI for Robust Feature Assessment in ML Grid
Version: 1.0.0
Depends: R (≥ 4.1.0)
Description: This R package introduces Weighted Mean SHapley Additive exPlanations (WMSHAP), a framework for estimating and combining SHAP values across a grid of fine-tuned machine learning models or the base learners of stacked ensembles. Rather than relying on a single best-performing model, WMSHAP combines SHAP values from multiple models, weighting their contributions according to predictive performance. Model weights can be calculated using R-squared for regression and classification models or, for binary classifiers, the area under the precision-recall curve (AUCPR), the area under the receiver operating characteristic curve (AUC), or the F2 score. The package also provides weighted confidence intervals for mean SHAP values, allowing feature importance to be evaluated across multiple models while accounting for variability in their contributions. This approach particularly improves accuracy if feature importance assessment at both local (individual) and global (feature) levels and improves the stability of feature importance rankins (see the GitHub repository for references and details). By integrating SHAP values across models, WMSHAP reduces reliance on any single potentially overfitted or biased model, providing a more comprehensive and generalizable assessment of feature importance. The package further implements hypothesis tests to assess the statistical significance of individual feature contributions, as well as comparative tests to evaluate differences in SHAP contributions between features. It also provides criteria for automatic feature selection across model grids and stacked ensembles, eliminating the need to specify an arbitrary number of top features. This approach supports a more representative assessment of feature importance by considering evidence from multiple models rather than relying exclusively on the best-performing model. The package also provides visualizations of SHAP values at individual, feature, and factor levels, including feature-selection plots based on weighted mean WMSHAP ratios. Supported SHAP backends include 'h2o', 'treeshap', 'shapr', and 'kernelshap', as well as native 'gpboost' explanations of the boosted-tree component.
License: MIT + file LICENSE
Encoding: UTF-8
Imports: ggplot2 (≥ 3.4.2), h2o (≥ 3.34.0.0), treeshap (≥ 0.4.0), kernelshap (≥ 0.9.1), pander (≥ 0.6.5), parallel, stats, utils
Suggests: gpboost (≥ 1.7.4), shapr (≥ 1.1.0), ranger
URL: https://github.com/haghish/shapley
BugReports: https://github.com/haghish/shapley/issues
Config/roxygen2/version: 8.1.0
NeedsCompilation: no
Packaged: 2026-11-05 23:47:34 UTC; haghish
Author: E. F. Haghish [aut, cre, cph]
Maintainer: E. F. Haghish <haghish@hotmail.com>
Repository: CRAN
Date/Publication: 2026-11-06 10:10:02 UTC

Weighted permutation test for paired difference of means

Description

Performs a weighted permutation test for the null hypothesis that the weighted mean of (var1 - var2) is zero.

Usage

feature.test(var1, var2, weights, n = 2000)

Arguments

var1

A numeric vector.

var2

A numeric vector of the same length as var1.

weights

A numeric vector of non-negative weights of the same length as var1 and var2.

n

Integer. Number of permutations (default 2000).

Value

A list with:

mean_wmshap_diff

Observed weighted mean difference (var1 - var2).

p_value

Monte Carlo permutation p-value.

Examples

## Not run: 
var1 <- rnorm(100)
var2 <- rnorm(100)
weights <- runif(100)
result <- shapley:::feature.test(var1, var2, weights)
result$mean_wmshap_diff
result$p_value

## End(Not run)

Extract model IDs from H2O AutoML or Grid objects

Description

Extracts model IDs from a "H2OAutoML" object (via the leaderboard) or from a "H2OGrid" object.

Usage

h2o.get_ids(h2oboard)

Arguments

h2oboard

An object inheriting from "H2OAutoML" or "H2OGrid".

Value

A character vector of trained model IDs.

Author(s)

E. F. Haghish

Examples


## Not run: 
library(h2o)
h2o.init(ignore_config = TRUE, nthreads = 2, bind_to_localhost = FALSE, insecure = TRUE)
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)
y <- "CAPSULE"
prostate[,y] <- as.factor(prostate[,y])  #convert to factor for classification
aml <- h2o.automl(y = y, training_frame = prostate, max_runtime_secs = 30)

# get the model IDs
ids <- h2o.get_ids(aml)

## End(Not run)

Normalize a vector based on specified minimum and maximum values

Description

This function normalizes a vector based on specified minimum and maximum values. If the minimum and maximum values are not specified, the function will use the minimum and maximum values of the vector (ignoring missing values).

Usage

normalize(x, min = NULL, max = NULL)

Arguments

x

numeric vector

min

minimum value

max

maximum value

Value

A numeric vector of the same length as x

Author(s)

E. F. Haghish

Examples

## Not run: 
# the function is not exported
normalize(c(0, 5, 10))
normalize(c(1, 1, 1))
normalize(c(NA, 2, 3))

## End(Not run)

Weighted Mean SHAP (WMSHAP) and Confidence Interval for Multiple Models

Description

Computes Weighted Mean SHAP ratios (WMSHAP) and confidence intervals to assess feature importance across a set of near-optimal models. Instead of reporting SHAP contributions for a single model, this function summarizes feature importance across multiple models and weights each model by a chosen performance metric. Models trained by the h2o machine learning platform use the native H2O procedure. Other models can be explained using treeshap, kernelshap, shapr, or the native gpboost SHAP method.

Usage

shapley(
  models,
  newdata,
  background_data = NULL,
  pred_fun = NULL,
  model_performance = NULL,
  performance_metric = "r2",
  backend = "auto",
  sample_size = nrow(newdata),
  bootstrap = 0L,
  n_core = -1L,
  cutoff = 0,
  plot = TRUE,
  top_n_features = NULL,
  target_class = NULL,
  shapr_method = "empirical",
  shapr_monte_carlo = 1000L,
  ...
)

Arguments

models

A model, which can be a fitted R model, a named list of models, an H2O AutoML object, H2O grid object, autoEnsemble object, a character vector of H2O model IDs. a wide range of R models are supported.

newdata

A data frame (for R models) or an H2OFrame (for H2O models), SHAP and WMSHAP contributions are computed on these observations.

background_data

Data frame containing reference observations. Internally passed as data to treeshap::unify(), bg_X to kernelshap::kernelshap(), and x_train to shapr::explain(). This argument is not needed for models trained with 'h2o' or 'gpboost' packages. background_data can be the training set with predictor columns only.

pred_fun

Optional prediction function of the form function(model, newdata). kernelshap package accepts multiple outputs; shapr package requires one numeric prediction per observation, selected using target_class when necessary. This argument is ignored for models trained with 'h2o' package, or when the 'treeshap' package is defined as the backend.

model_performance

Named numeric vector containing one model performance estimate per non-H2O model. Use comparable cross-validation estimates, such as AUCPR, for all models. Names should correspond to model list names (see the examples).

performance_metric

Character. Performance metric used to weight models. Options are "r2" (regression), "aucpr", "auc", and "f2" (classification metrics).

backend

Character. Specifies the backend used to compute SHAP values. the default is "auto", which detects the appropriate backend based on the models' class (assuming they are trained with a single learner). Alternative available options are "kernelshap" (universal), "h2o" (only for h2o models), "gpboost" (only for gpboost models), "treeshap" (for supported tree-based models), and "shapr". Each backend uses a different method to estimate SHAP values. The gpboost and shapr packages are optional dependencies that must be installed separately. Native GPBoost SHAP values explain only the boosted-tree component, excluding Gaussian-process and random effects. Note that these backends compute SHAP values in different ways.

sample_size

Integer. Number of rows in newdata used for SHAP assessment. Defaults to number of rows in the newdata. Decreasing this can increase the computation speed, but decreases accuracy and stability of WMSHAP. Note that for 'shapr' package, this argument defines of Monte Carlo samples passed to shapr::explain().

bootstrap

Integer. Number of observation-level bootstrap replications used to estimate sampling variability in WMSHAP and Between-Model Variability (BMV). The default is 0L, which disables bootstrapping. For each bootstrap sample, the function recalculates the model-specific SHAP importance ratios, their performance-weighted mean, and between-model standard deviation. The root mean square of the bootstrap between-model standard deviations is used to construct a descriptive BMV interval. The function also estimates bootstrap standard errors and percentile 95% confidence intervals for WMSHAP and the between-model SD. Note that bootstrap is only reasonable for global and factor (domain) level WMSHAP, and is not useful for individual-level (local) WMSHAP. Moreover, the current implementation of bootstrap is not suitable for hierarchical data (i.e., gpboost models).

n_core

Integer. Number of parallel R workers for computing SHAP contributions across fitted models. The default -1 uses all detected physical CPU cores, limited to the number of models. Use 1 for sequential computation or a positive integer for a fixed number of workers. Each worker explains one model at a time. This argument does not control internal parallelism in the SHAP backends.

cutoff

Numeric. WMSHAP mean cutoff used to identify and select important features. The recommended default is 1 divided by the number of predictors in the model.

plot

Logical. If TRUE, plots the WMSHAP summary using shapley.plot().

top_n_features

Integer or NULL. If specified, select the top N features. All features are used to compute global importance ratios; this argument controls feature selection afterwards.

target_class

Name or index of the prediction class to explain for models returning more than one prediction per observation. The second class is used by default in binary classification.

shapr_method

Character vector. Method used by shapr::explain() to estimate feature contributions. Available methods are "empirical", "gaussian", "copula", "ctree", "arf", "vaeac", "categorical", "independence", "timeseries", "regression_separate", and "regression_surrogate". The default is "empirical".

shapr_monte_carlo

Integer. Number of Monte Carlo samples used to estimate SHAP values when backend = "shapr". Larger values may improve the precision of SHAP estimates but increase computation time. The default is 1000L. This argument is ignored when other SHAP backends are used.

...

arguments that are used internally for compatibility with older versions

Details

If the near-optimal set contains models requiring different non-H2O backends, kernelshap is used for all included models to maintain a common calculation method. The shapr backend can be selected explicitly to calculate dependence-aware SHAP values; it requires installing the optional shapr package.

The function works as follows:

  1. For each model, SHAP contributions are computed on newdata.

  2. For each model, feature-level absolute SHAP contributions are aggregated and converted to a ratio (share of total absolute SHAP across features).

  3. Models are weighted by a normalized performance metric (e.g., "r2" for regression or "auc" / "aucpr" for classification).

  4. The weighted mean SHAP ratio (WMSHAP) is computed for each feature, along with an interval across models.

Value

An object of class "shapley" (a named list) containing:

ids

Character vector of model IDs originally supplied or extracted.

included_models

Character vector of included model IDs.

ignored_models

Data frame of excluded models and their performance.

weights

Legacy rescaled model weights (sum equals the number of models).

model_weights

Normalized model-performance weights (sum equals 1).

results

Data frame of row-level SHAP contributions across models.

summaryShaps

Data frame of global WMSHAP ratios and intervals as well as bootstrap summary statistics, if applied.

selectedFeatures

Character vector of selected important features.

feature_importance

List of feature importance values by model.

contributionPlot

Plot data for local SHAP contributions.

plot

A ggplot object when plot = TRUE.

localShaps

Individual-level WMSHAP and between-model intervals for all supported models, including H2O.

model_backends

Selected SHAP backend for each included non-H2O model.

explanation_scope

For native GPBoost SHAP, "gpboost_tree_component". This excludes random effects and Gaussian-process contributions.

Note

SHAP values from different backends are only directly comparable if they explain the same prediction output and use compatible feature-dependence assumptions. The intervals summarize between-model variation, not population uncertainty. The H2O implementation directly retrieves local SHAP contributions using h2o::h2o.predict_contributions(), without creating a SHAP plot per model. Internal performance weights sum to 1; legacy output weights sums to the number of included models for compatibility, while model_weights contains the normalized weights. Sampling uncertainty from shapr is different and is not included in these intervals. For GPBoost mixed-effects models, native SHAP contributions describe the boosted-tree component only. They do not explain group-level random effects, Gaussian-process effects, or the entire response prediction. Model weights may be based on full-model cross-validation performance, but the explained component remains the boosted-tree function.

Author(s)

E. F. Haghish

Examples

## Not run: 

# ============================================================
# An Example Using h2o Package
# ============================================================
library(shapley)
library(h2o)
h2o.init(ignore_config = TRUE, nthreads = 2,
         bind_to_localhost = FALSE, insecure = TRUE)
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)
y <- "CAPSULE"
prostate[, y] <- as.factor(prostate[, y])
grid <- h2o.grid(algorithm = "gbm", y = y,
                 training_frame = prostate,
                 hyper_params = list(ntrees = seq(1, 50)),
                 grid_id = "ensemble_grid", nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

result <- shapley(models = grid,
                  newdata = prostate,
                  performance_metric = "aucpr",
                  bootstrap = 100)

# ============================================================
# An Example Using ranger Package
# ============================================================
library(tidymodels)
library(ranger)
library(shapley)

# 1. prepare the dataset for training and testing
prostate <- read.csv(system.file("extdata", "prostate.csv", package = "h2o"))
prostate$ID <- NULL
prostate$CAPSULE <- factor(prostate$CAPSULE, levels = c(1, 0))

set.seed(2026)
split <- initial_split(prostate, prop = 0.80, strata = CAPSULE)
train <- training(split)
test <- testing(split)
predictors <- setdiff(names(train), "CAPSULE")

# 2. define a grid for randomforest configurations.
rf <- rand_forest(
  trees = tune(),
  mtry = tune(),
  min_n = tune()
) |>
  set_engine("ranger", num.threads = 1L, seed = 2026L) |>
  set_mode("classification")

wf <- workflow() |>
  add_model(rf) |>
  add_formula(CAPSULE ~ .)

grid <- tidyr::crossing(
  trees = c(100L, 200L, 300L, 400L, 500L),
  mtry = c(2L, 3L, 4L, 5L, 6L, 7L),
  min_n = c(2L, 4L, 6L, 8L)
)

# 3. evaluate the configurations using 10-fold CV and AUCPR as evaluation metrics.
folds <- vfold_cv(train, v = 10, strata = CAPSULE)

cv <- tune_grid(
  wf,
  resamples = folds,
  grid = grid,
  metrics = metric_set(pr_auc),
  control = control_grid(verbose = TRUE, save_pred = FALSE,
                         save_workflow = FALSE, event_level = "first")
)

# mean AUCPR across CV folds
scores <- collect_metrics(cv) |>
  dplyr::filter(.metric == "pr_auc", is.finite(mean)) |>
  dplyr::arrange(dplyr::desc(mean))

# 4: select the near-optimal configurations. for example,
# retain top models with at least 98\
near <- scores |>
  dplyr::filter(mean >= 0.98 * max(scores$mean))

# 5: fit each retained configuration on the complete training data.
models <- lapply(seq_len(nrow(near)), function(i) {
  pars <- near[i, c("trees", "mtry", "min_n")]
  fitted <- finalize_workflow(wf, pars) |>
    fit(data = train)
  workflows::extract_fit_engine(fitted)
})

names(models) <- paste0("model", seq_along(models))
cv_aucpr <- stats::setNames(near$mean, names(models))

# computing WMSHAP using treeshap package
# ------------------------------------------------------------
result_tree <- shapley(
  models = models,
  newdata = test[, predictors],
  background_data = train[, predictors],
  model_performance = cv_aucpr,
  performance_metric = "aucpr",
  backend = "treeshap",
  target_class = "1",
  bootstrap = 100,
  plot = TRUE
)

# computing WMSHAP using kernelshap package
# ------------------------------------------------------------
pred_fun <- function(model, X) {
  requireNamespace("ranger", quietly = TRUE)
  as.numeric(stats::predict(model, data = as.data.frame(X), num.threads = 1)$predictions[, "1"])
}

result_kernel <- shapley(
  models = models,
  newdata = test[, predictors],
  background_data = train[, predictors],
  model_performance = cv_aucpr,
  performance_metric = "aucpr",
  backend = "kernelshap",
  pred_fun = pred_fun,
  target_class = "1",
  bootstrap = 100,
  n_core = 4,
  plot = FALSE
)

# omputing WMSHAP using shapr package
# ------------------------------------------------------------
result_shapr <- shapley(
  models = models,
  newdata = test[, predictors],
  background_data = train[, predictors],
  model_performance = cv_aucpr,
  performance_metric = "aucpr",
  backend = "shapr",
  pred_fun = pred_fun,
  target_class = "1",
  shapr_method = "empirical",
  bootstrap = 100,
  shapr_monte_carlo = 1000,
  plot = TRUE
)

# WMSHAP global results (feature-level)
# ------------------------------------------------------------
print(result_tree$summaryShaps)
print(result_kernel$summaryShaps)
print(result_shapr$summaryShaps)

# WMSHAP local results (individual-level)
# ------------------------------------------------------------
print(head(result_tree$localShaps))
print(head(result_kernel$localShaps))
print(head(result_shapr$localShaps))


## End(Not run)

Compute and plot weighted mean SHAP contributions at group level (factors or domains)

Description

Aggregates SHAP contributions across user-defined domains (groups of features), computes weighted means and descriptive between-model intervals, and returns a plot plus summary tables.

Usage

shapley.domain(
  shapley,
  domains,
  plot = TRUE,
  print = FALSE,
  colorcode = NULL,
  xlab = "Domains"
)

Arguments

shapley

Object of class "shapley", as returned by the 'shapley' function

domains

Named list of character vectors. Each element name is a domain name; each element value is a character vector of feature names assigned to that domain.

plot

Logical. If TRUE, a bar plot of domain WMSHAP contributions is created.

print

Logical. If TRUE, prints the domain WMSHAP summary table.

colorcode

Character vector for specifying the color names for each domain in the plot.

xlab

Character. Specify the ggplot 'xlab' label in the plot (default is "Domains")

Author(s)

E. F. Haghish

Independent weighted bootstrap test for the difference between two domain-level WMSHAP means

Description

Compares the weighted mean contribution of a specific domain between two "shapley" objects. The observed effect is defined as the difference between the weighted mean domain contributions in the two shapley objects.

Usage

shapley.domain.compare.test(
  shapley1,
  shapley2,
  domains,
  domain,
  n = 2000,
  conf.level = 0.95,
  report = TRUE,
  seed = NULL
)

Arguments

shapley1

Object of class "shapley".

shapley2

Object of class "shapley".

domains

Named list of domains passed to shapley.domain(). The list must contain at least 2 domains.

domain

Character string naming the domain to be compared between the two "shapley" objects.

n

Integer, number of bootstrap resamples used for the confidence interval and p-value calculation. Must be at least 100.

conf.level

Numeric scalar between 0 and 1 giving the confidence level for the bootstrap confidence interval.

report

logical; if TRUE (default) the results are explained.

seed

Optional integer seed for reproducibility.

Details

The function carries out an independent weighted bootstrap comparison of domain-level weighted mean SHAP (WMSHAP) values. For the selected domain, model-specific domain contribution ratios are extracted from each "shapley" object. The model weights stored in each object are then normalized to sum to 1 and used to compute the weighted mean domain contribution in each sample.

The observed difference is defined as the weighted mean domain contribution in shapley1 minus the weighted mean domain contribution in shapley2. A nonparametric weighted bootstrap is then used to estimate the confidence interval of this difference. Within each sample, bootstrap resamples are drawn with replacement from the model-specific domain contributions, using the normalized model weights as sampling probabilities. The bootstrap distribution of the difference between the two resampled means is used to construct the confidence interval.

To test the null hypothesis that the two weighted mean domain contributions are equal, the function generates a pooled weighted bootstrap null distribution. The domain contributions from both samples are pooled, the corresponding weights are combined and renormalized, and bootstrap samples of the original sample sizes are repeatedly drawn from this pooled weighted distribution. The two-sided p-value is computed as the proportion of null bootstrap differences that are at least as extreme as the observed difference, using the finite-sample correction (b + 1) / (n + 1).

For reporting in APA style, the following values are usually sufficient: the domain name, the weighted mean domain contribution in each sample, the difference between the two weighted means, the confidence interval, and the p-value. A typical report can be written as follows: For the [domain] domain, the weighted mean contribution was M_w = x.xx in Sample 1 and M_w = x.xx in Sample 2, with a difference of \DeltaWMSHAP = x.xx, 95% CI [x.xx, x.xx], p = .xxx, based on n weighted bootstrap resamples. Because this function uses bootstrap inference rather than a t test, no t statistic or degrees of freedom are reported.

Value

A list containing the selected domain, the weighted mean domain contribution in the first sample (mean_wmshap_1), the weighted mean domain contribution in the second sample (mean_wmshap_2), the difference between the two weighted means (mean_wmshap_diff), the lower and upper confidence interval bounds (ci_lower, ci_upper), and the two-sided bootstrap p_value.

Author(s)

E. F. Haghish

Weighted permutation test for difference between two domains

Description

Computes domain-level contribution ratios (via shapley.domain()) and tests whether two domains differ using a weighted paired permutation test across models.

Usage

shapley.domain.test(shapley, domains, n = 2000)

Arguments

shapley

Object of class "shapley", as returned by the 'shapley' function

domains

A named list of length 2. Each element is a character vector of feature names defining a domain; the two element names are the domain labels to be compared.

n

Integer, number of permutations (default 2000)

Value

A list with mean_wmshap_diff (observed weighted mean difference) and p_value.

Author(s)

E. F. Haghish

Examples


## Not run: 
# load the required libraries for building the base-learners and the ensemble models
library(h2o)            #shapley supports h2o models
library(autoEnsemble)   #autoEnsemble models, particularly useful under severe class imbalance
library(shapley)

# initiate the h2o server
h2o.init(ignore_config = TRUE, nthreads = 2, bind_to_localhost = FALSE, insecure = TRUE)

# upload data to h2o cloud
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)

### H2O provides 2 types of grid search for tuning the models, which are
### AutoML and Grid. Below, I demonstrate how weighted mean shapley values
### can be computed for both types.

set.seed(10)

#######################################################
### PREPARE AutoML Grid (takes a couple of minutes)
#######################################################
# run AutoML to tune various models (GBM) for 60 seconds
y <- "CAPSULE"
prostate[,y] <- as.factor(prostate[,y])  #convert to factor for classification
aml <- h2o.automl(y = y, training_frame = prostate, max_runtime_secs = 120,
                 include_algos=c("GBM"),

                 # this setting ensures the models are comparable for building a meta learner
                 seed = 2023, nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

### call 'shapley' function to compute the weighted mean and weighted confidence intervals
### of SHAP values across all trained models.
### Note that the 'newdata' should be the testing dataset!
result <- shapley(models = aml, newdata = prostate, plot = TRUE)

#######################################################
### Significance testing of contributions of two domains (or latent factors)
#######################################################
domains = list(Demographic = c("RACE", "AGE"),
               Cancer = c("VOL", "PSA", "GLEASON"))
shapley.domain.test(result, domains = domains, n=5000)

## End(Not run)

Weighted permutation test for WMSHAP difference between two features

Description

Performs a weighted paired permutation test to assess whether two features have different contributions (e.g., weighted mean SHAP, referred to as WMSHAP) across models in a shapley object.

Usage

shapley.feature.test(shapley, features, n = 2000)

Arguments

shapley

object of class "shapley", as returned by the 'shapley' function

features

Character vector of length 2 giving the names of the two features to compare.

n

Integer. Number of permutations (default 2000).

Value

A list with mean_wmshap_diff (observed weighted mean difference) and p_value.

Author(s)

E. F. Haghish

Examples


## Not run: 
# load the required libraries for building the base-learners and the ensemble models
library(h2o)            #shapley supports h2o models
library(autoEnsemble)   #autoEnsemble models, particularly useful under severe class imbalance
library(shapley)

# initiate the h2o server
h2o.init(ignore_config = TRUE, nthreads = 2, bind_to_localhost = FALSE, insecure = TRUE)

# upload data to h2o cloud
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)

### H2O provides 2 types of grid search for tuning the models, which are
### AutoML and Grid. Below, I demonstrate how weighted mean shapley values
### can be computed for both types.

set.seed(10)

#######################################################
### PREPARE AutoML Grid (takes a couple of minutes)
#######################################################
# run AutoML to tune various models (GBM) for 60 seconds
y <- "CAPSULE"
prostate[,y] <- as.factor(prostate[,y])  #convert to factor for classification
aml <- h2o.automl(y = y, training_frame = prostate, max_runtime_secs = 120,
                 include_algos=c("GBM"),

                 # this setting ensures the models are comparable for building a meta learner
                 seed = 2023, nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

### call 'shapley' function to compute the weighted mean and weighted confidence intervals
### of SHAP values across all trained models.
### Note that the 'newdata' should be the testing dataset!
result <- shapley(models = aml, newdata = prostate, plot = TRUE)

#######################################################
### Significance testing of contributions of two features
#######################################################

shapley.feature.test(result, features = c("GLEASON", "PSA"), n = 5000)

## End(Not run)

Plot Global WMSHAP Values

Description

Plot global feature importance using Weighted Mean SHAP ratios or plot individual-level weighted mean SHAP contributions. The bootstrap percentile confidence intervals are distinct from descriptive between-model and BMV intervals. The WMSHAP contribution plot uses density-based jitter and the classic plot styling of the original H2O implementation of the shapley package.

Usage

shapley.plot(
  shapley,
  plot = "bar",
  method = "mean",
  cutoff = 0.01,
  top_n_features = NULL,
  features = NULL,
  legendstyle = "continuous",
  scale_colour_gradient = NULL,
  labels = NULL,
  interval = "auto"
)

Arguments

shapley

An object returned by shapley().

plot

Plot type: "bar" or "wmshap".

method

Feature selection method: "mean", "lowerCI", "bootstrap_lowerCI", or "bmv_lower".

cutoff

Minimum global WMSHAP ratio or lower confidence interval.

top_n_features

Optional number of features to plot.

features

Optional character vector of features to plot.

legendstyle

Legend type used for the contribution plot.

scale_colour_gradient

Optional colour vector for the contribution plot.

labels

Optional feature labels.

interval

Interval type to display: "auto", "bmv", "bootstrap", "between_model", "both", or "all".

Value

A ggplot object.

Plot Individual-Level WMSHAP Values

Description

Computes and plots model-performance-weighted SHAP contributions and between-model intervals for an individual observation.

Usage

shapley.row.plot(
  shapley,
  row_index,
  top_n_features = NULL,
  features = NULL,
  nonzeroCI = FALSE,
  plot = TRUE,
  print = FALSE
)

Arguments

shapley

An object returned by shapley().

row_index

Index of the observation to explain.

top_n_features

Optional number of important predictors to plot.

features

Optional predictor names.

nonzeroCI

If TRUE, show only contributions whose intervals do not contain zero.

plot

Logical. If TRUE, print the plot.

print

Logical. If TRUE, print the results table.

Value

A list containing the summary table and a ggplot object.

Create WMSHAP Summary Table Based on the Given Criterion

Description

#' Generates a summary table of weighted mean SHAP ratios (WMSHAP) and confidence intervals for each feature based on a weighted SHAP analysis. The function filters the SHAP summary table (from a shapley object) by selecting features that meet or exceed a specified cutoff using a selection method (default "mean").

The output is sorted by WMSHAP and formatted as either a markdown table (via pander) or a data frame.

Usage

shapley.table(
  shapley,
  method = "mean",
  cutoff = 0.01,
  round = 3,
  exclude_features = NULL,
  dict = NULL,
  markdown.table = TRUE,
  split.tables = 120,
  split.cells = 50,
  interval = "auto"
)

Arguments

shapley

A shapley object returned by shapley that contains a data frame summaryShaps.

method

Character. The column name in summaryShaps used for feature selection. Default is "mean", which selects important features which have weighted mean shap ratio (WMSHAP) higher than the specified cutoff. Other alternative is "lowerCI", which selects features which their lower bound of confidence interval is higher than the cutoff.

cutoff

Numeric. The threshold cutoff for the selection method; only features with a value in the method column greater than or equal to this value are retained. Default is 0.01.

round

Integer. The number of decimal places to round the SHAP mean and confidence interval values. Default is 3.

exclude_features

Character vector. A vector of feature names to be excluded from the summary table. Default is NULL.

dict

A data frame containing at least two columns named "name" and "description". If provided, the function uses this dictionary to add human-readable feature descriptions. Default is NULL.

markdown.table

Logical. If TRUE, the output is formatted as a markdown table using the pander package; otherwise, a data frame is returned. Default is TRUE.

split.tables

Integer. Controls table splitting in pander(). Default is 120.

split.cells

Integer. Controls cell splitting in pander(). Default is 50.

interval

Interval type to display: "auto", "bmv", "bootstrap", and "between_model".

Value

If markdown.table = TRUE, returns a markdown table (invisibly) showing two columns: "Description" and "WMSHAP". If markdown.table = FALSE, returns a data frame with these columns.

Author(s)

E. F. Haghish

Examples

## Not run: 
# load the required libraries for building the base-learners and the ensemble models
library(h2o)            #shapley supports h2o models
library(shapley)

# initiate the h2o server
h2o.init(ignore_config = TRUE, nthreads = 2, bind_to_localhost = FALSE, insecure = TRUE)

# upload data to h2o cloud
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)

set.seed(10)

### H2O provides 2 types of grid search for tuning the models, which are
### AutoML and Grid. Below, I demonstrate how weighted mean shapley values
### can be computed for both types.

#######################################################
### PREPARE AutoML Grid (takes a couple of minutes)
#######################################################
# run AutoML to tune various models (GBM) for 60 seconds
y <- "CAPSULE"
prostate[,y] <- as.factor(prostate[,y])  #convert to factor for classification
aml <- h2o.automl(y = y, training_frame = prostate, max_runtime_secs = 120,
                 include_algos=c("GBM"),

                 # this setting ensures the models are comparable for building a meta learner
                 seed = 2023, nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

### call 'shapley' function to compute the weighted mean and weighted confidence intervals
### of SHAP values across all trained models.
### Note that the 'newdata' should be the testing dataset!
result <- shapley(models = aml, newdata = prostate, performance_metric = "aucpr", plot = TRUE)

#######################################################
### PREPARE H2O Grid (takes a couple of minutes)
#######################################################
# make sure equal number of "nfolds" is specified for different grids
grid <- h2o.grid(algorithm = "gbm", y = y, training_frame = prostate,
                 hyper_params = list(ntrees = seq(1,50,1)),
                 grid_id = "ensemble_grid",

                 # this setting ensures the models are comparable for building a meta learner
                 seed = 2023, fold_assignment = "Modulo", nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

result2 <- shapley(models = grid, newdata = prostate, performance_metric = "aucpr", plot = TRUE)

# get the output as a Markdown table:
md_table <- shapley.table(shapley = result2,
                          method = "mean",
                          cutoff = 0.01,
                          round = 3,
                          markdown.table = TRUE)
head(md_table)

## End(Not run)

Flag and rank features by WMSHAP cutoffs

Description

This function applies different criteria simultaniously to identify the most important features in a model. The criteria include: 1) minimum limit of lower weighted confidence intervals of SHAP values relative to the feature with highest SHAP value. 2) minimum limit of percentage of weighted mean SHAP values relative to over all SHAP values of all features. These are specified with two different cutoff values.

Usage

shapley.top(shapley, mean = 0.01, lowerCI = 0.01)

Arguments

shapley

object of class 'shapley', as returned by the 'shapley' function

mean

Numeric. specifying the cutoff of weighted mean SHAP ratio (WMSHAP). The default is 0.01. Lower values will be more generous in defining "importance", while higher values are more restrictive. However, these default values are not generalizable to all situations and algorithms.

lowerCI

numeric. Specifying the limit of lower bound of 95% WMSHAP The default is 0.01. Lower values will be more generous in defining "importance", while higher values are more restrictive. However, these default values are not generalizable to all situations and algorithms.

Value

data.frame of selected features

Author(s)

E. F. Haghish

Examples


## Not run: 
# load the required libraries for building the base-learners and the ensemble models
library(h2o)            #shapley supports h2o models
library(shapley)

# initiate the h2o server
h2o.init(ignore_config = TRUE, nthreads = 2, bind_to_localhost = FALSE, insecure = TRUE)

# upload data to h2o cloud
prostate_path <- system.file("extdata", "prostate.csv", package = "h2o")
prostate <- h2o.importFile(path = prostate_path, header = TRUE)

### H2O provides 2 types of grid search for tuning the models, which are
### AutoML and Grid. Below, I demonstrate how weighted mean shapley values
### can be computed for both types.

set.seed(10)

#######################################################
### PREPARE AutoML Grid (takes a couple of minutes)
#######################################################
# run AutoML to tune various models (GBM) for 60 seconds
y <- "CAPSULE"
prostate[,y] <- as.factor(prostate[,y])  #convert to factor for classification
aml <- h2o.automl(y = y, training_frame = prostate, max_runtime_secs = 120,
                 include_algos=c("GBM"),

                 # this setting ensures the models are comparable for building a meta learner
                 seed = 2023, nfolds = 10,
                 keep_cross_validation_predictions = TRUE)

### call 'shapley' function to compute the weighted mean and weighted confidence intervals
### of SHAP values across all trained models.
### Note that the 'newdata' should be the testing dataset!
result <- shapley(models = aml, newdata = prostate, plot = TRUE)

#######################################################
### Select top features
#######################################################
shapley.top(result, mean = 0.005, lowerCI = 0.01)

## End(Not run)