Revert "Rename things"

Author: mayer79

CRAN-SUBMISSION

@@ -0,0 +1,3 @@
+Version: 0.3.7
+Date: 2023-05-17 06:52:31 UTC
+SHA: b6e4ce87f93a54e5c451cd06315ab810bb29eb8a

DESCRIPTION

@@ -2,9 +2,17 @@ Package: kernelshap
 Title: Kernel SHAP
 Version: 0.3.8
 Authors@R: c(
-    person("Michael", "Mayer", , "[email protected]", role = c("aut", "cre"))
+    person("Michael", "Mayer", , "[email protected]", role = c("aut", "cre")),
+    person("David", "Watson", , "[email protected]", role = "aut"),
+    person("Przemyslaw", "Biecek", , "[email protected]", role = "ctb",
+           comment = c(ORCID = "0000-0001-8423-1823"))
   )
-Description: Implementation of ...  The package plays well together
+Description: Efficient implementation of Kernel SHAP, see Lundberg and Lee
+    (2017) <https://dl.acm.org/doi/10.5555/3295222.3295230>, and Covert
+    and Lee (2021) <http://proceedings.mlr.press/v130/covert21a>.  For
+    models with up to eight features, the results are exact regarding the
+    selected background data.  Otherwise, an almost exact hybrid algorithm
+    involving iterative sampling is used.  The package plays well together
     with meta-learning packages like 'tidymodels', 'caret' or 'mlr3'.
     Visualizations can be done using the R package 'shapviz'.
 License: GPL (>= 2)

R/exact.R

@@ -9,7 +9,19 @@ input_exact <- function(p) {
   Z <- exact_Z(p)
   # Each Kernel weight(j) is divided by the number of vectors z having sum(z) = j
   w <- kernel_weights(p) / choose(p, 1:(p - 1L))
-  list(Z = Z, w = w[rowSums(Z)])
+  list(Z = Z, w = w[rowSums(Z)], A = exact_A(p))
+}
+
+# Calculates exact A. Notice the difference to the off-diagnonals in the Supplement of 
+# Covert and Lee (2021). Credits to David Watson for figuring out the correct formula,
+# see our discussions in https://github.com/ModelOriented/kernelshap/issues/22
+exact_A <- function(p) {
+  S <- 1:(p - 1L)
+  c_pr <- S * (S - 1) / p / (p - 1)
+  off_diag <- sum(kernel_weights(p) * c_pr)
+  A <- matrix(off_diag, nrow = p, ncol = p)
+  diag(A) <- 0.5
+  A
 }
 
 # Creates (2^p-2) x p matrix with all on-off vectors z of length p
@@ -53,10 +65,10 @@ input_partly_exact <- function(p, deg) {
   if (p < 2L * deg) {
     stop("p must be >=2*deg")
   }
-
+  
   kw <- kernel_weights(p)
   Z <- w <- vector("list", deg)
-
+  
   for (k in seq_len(deg)) {
     Z[[k]] <- partly_exact_Z(p, k = k)
     n <- nrow(Z[[k]])
@@ -65,6 +77,40 @@ input_partly_exact <- function(p, deg) {
   }
   w <- unlist(w, recursive = FALSE, use.names = FALSE)
   Z <- do.call(rbind, Z)
-
-  list(Z = Z, w = w)
+  
+  list(Z = Z, w = w, A = crossprod(Z, w * Z))
 }
+
+# Case p = 1 returns exact Shapley values
+case_p1 <- function(n, nms, v0, v1, X, verbose) {
+  txt <- "Exact Shapley values (p = 1)"
+  if (verbose) {
+    message(txt)
+  }
+  S <- v1 - v0[rep(1L, n), , drop = FALSE]
+  SE <- matrix(numeric(n), dimnames = list(NULL, nms))
+  if (ncol(v1) > 1L) {
+    SE <- replicate(ncol(v1), SE, simplify = FALSE)
+    S <- lapply(
+      asplit(S, MARGIN = 2L), function(M) as.matrix(M, dimnames = list(NULL, nms))
+    )
+  } else {
+    colnames(S) <- nms      
+  }
+  out <- list(
+    S = S, 
+    X = X, 
+    baseline = as.vector(v0), 
+    SE = SE, 
+    n_iter = integer(n), 
+    converged = rep(TRUE, n),
+    m = 0L,
+    m_exact = 0L,
+    prop_exact = 1,
+    exact = TRUE,
+    txt = txt,
+    predictions = v1
+  )
+  class(out) <- "kernelshap"
+  out
+}

R/kernelshap.R

@@ -1,39 +1,39 @@
 #' Print Method
+#' 
+#' Prints the first two rows of the matrix (or matrices) of SHAP values. 
 #'
-#' Prints the first two rows of the matrix (or matrices) of SHAP values.
-#'
-#' @param x An object of class "permshap".
+#' @param x An object of class "kernelshap".
 #' @param n Maximum number of rows of SHAP values to print.
 #' @param ... Further arguments passed from other methods.
 #' @returns Invisibly, the input is returned.
 #' @export
 #' @examples
 #' fit <- stats::lm(Sepal.Length ~ ., data = iris)
-#' s <- permshap(fit, iris[1:3, -1], bg_X = iris[-1])
+#' s <- kernelshap(fit, iris[1:3, -1], bg_X = iris[-1])
 #' s
-#' @seealso [permshap()]
-print.permshap <- function(x, n = 2L, ...) {
+#' @seealso [kernelshap()]
+print.kernelshap <- function(x, n = 2L, ...) {
   cat("SHAP values of first", n, "observations:\n")
   print(head_list(getElement(x, "S"), n = n))
   invisible(x)
 }
 
 #' Summary Method
 #'
-#' @param object An object of class "permshap".
-#' @param compact Set to `TRUE` to hide printing the top n SHAP values,
-#'   standard errors and feature values.
-#' @param n Maximum number of rows of SHAP values, standard errors and feature values
+#' @param object An object of class "kernelshap".
+#' @param compact Set to `TRUE` to hide printing the top n SHAP values, 
+#'   standard errors and feature values. 
+#' @param n Maximum number of rows of SHAP values, standard errors and feature values 
 #'   to print.
 #' @param ... Further arguments passed from other methods.
 #' @returns Invisibly, the input is returned.
 #' @export
 #' @examples
 #' fit <- stats::lm(Sepal.Length ~ ., data = iris)
-#' s <- permshap(fit, iris[1:3, -1], bg_X = iris[-1])
+#' s <- kernelshap(fit, iris[1:3, -1], bg_X = iris[-1])
 #' summary(s)
-#' @seealso [permshap()]
-summary.permshap <- function(object, compact = FALSE, n = 2L, ...) {
+#' @seealso [kernelshap()]
+summary.kernelshap <- function(object, compact = FALSE, n = 2L, ...) {
   cat(getElement(object, "txt"))
 
   S <- getElement(object, "S")
@@ -68,19 +68,19 @@ summary.permshap <- function(object, compact = FALSE, n = 2L, ...) {
   invisible(object)
 }
 
-#' Check for permshap
+#' Check for kernelshap
 #'
-#' Is object of class "permshap"?
+#' Is object of class "kernelshap"?
 #'
 #' @param object An R object.
-#' @returns `TRUE` if `object` is of class "permshap", and `FALSE` otherwise.
+#' @returns `TRUE` if `object` is of class "kernelshap", and `FALSE` otherwise.
 #' @export
 #' @examples
 #' fit <- stats::lm(Sepal.Length ~ ., data = iris)
-#' s <- permshap(fit, iris[1:2, -1], bg_X = iris[-1])
-#' is.permshap(s)
-#' is.permshap("a")
-#' @seealso [permshap()]
-is.permshap <- function(object){
-  inherits(object, "permshap")
+#' s <- kernelshap(fit, iris[1:2, -1], bg_X = iris[-1])
+#' is.kernelshap(s)
+#' is.kernelshap("a")
+#' @seealso [kernelshap()]
+is.kernelshap <- function(object){
+  inherits(object, "kernelshap")
 }

R/methods.R

@@ -1,8 +1,8 @@
 # Kernel SHAP algorithm for a single row x
 # If exact, a single call to predict() is necessary.
 # If sampling is involved, we need at least two additional calls to predict().
-permshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact, deg,
-                         paired, m, tol, max_iter, v0, precalc, ...) {
+kernelshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact, deg, 
+                           paired, m, tol, max_iter, v0, precalc, ...) {
   p <- length(feature_names)
 
   # Calculate A_exact and b_exact
@@ -12,28 +12,28 @@ permshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact, de
     Z <- precalc[["Z"]]                                           #  (m_ex x p)
     m_exact <- nrow(Z)
     v0_m_exact <- v0[rep(1L, m_exact), , drop = FALSE]            #  (m_ex x K)
-
+    
     # Most expensive part
     vz <- get_vz(                                                 #  (m_ex x K)
       X = x[rep(1L, times = nrow(bg_X_exact)), , drop = FALSE],   #  (m_ex*n_bg x p)
       bg = bg_X_exact,                                            #  (m_ex*n_bg x p)
       Z = Z,                                                      #  (m_ex x p)
-      object = object,
+      object = object, 
       pred_fun = pred_fun,
       feature_names = feature_names,
-      w = bg_w,
+      w = bg_w, 
       ...
     )
     # Note: w is correctly replicated along columns of (vz - v0_m_exact)
     b_exact <- crossprod(Z, precalc[["w"]] * (vz - v0_m_exact))   #  (p x K)
-
+    
     # Some of the hybrid cases are exact as well
     if (exact || trunc(p / 2) == deg) {
       beta <- solver(A_exact, b_exact, constraint = v1 - v0)      #  (p x K)
-      return(list(beta = beta, sigma = 0 * beta, n_iter = 1L, converged = TRUE))
+      return(list(beta = beta, sigma = 0 * beta, n_iter = 1L, converged = TRUE))  
     }
-  }
-
+  } 
+  
   # Iterative sampling part, always using A_exact and b_exact to fill up the weights
   bg_X_m <- precalc[["bg_X_m"]]                                   #  (m*n_bg x p)
   X <- x[rep(1L, times = nrow(bg_X_m)), , drop = FALSE]           #  (m*n_bg x p)
@@ -48,32 +48,32 @@ permshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact, de
     A_exact <- A_sum
     b_exact <- b_sum
   }
-
+  
   while(!isTRUE(converged) && n_iter < max_iter) {
     n_iter <- n_iter + 1L
     input <- input_sampling(p = p, m = m, deg = deg, paired = paired)
     Z <- input[["Z"]]
-
+      
     # Expensive                                                              #  (m x K)
     vz <- get_vz(
-      X = X,
-      bg = bg_X_m,
-      Z = Z,
-      object = object,
-      pred_fun = pred_fun,
-      feature_names = feature_names,
-      w = bg_w,
+      X = X, 
+      bg = bg_X_m, 
+      Z = Z, 
+      object = object, 
+      pred_fun = pred_fun, 
+      feature_names = feature_names, 
+      w = bg_w, 
       ...
     )
-
+    
     # The sum of weights of A_exact and input[["A"]] is 1, same for b
     A_temp <- A_exact + input[["A"]]                                         #  (p x p)
     b_temp <- b_exact + crossprod(Z, input[["w"]] * (vz - v0_m))             #  (p x K)
     A_sum <- A_sum + A_temp                                                  #  (p x p)
     b_sum <- b_sum + b_temp                                                  #  (p x K)
-
-    # Least-squares with constraint that beta_1 + ... + beta_p = v_1 - v_0.
-    # The additional constraint beta_0 = v_0 is dealt via offset
+    
+    # Least-squares with constraint that beta_1 + ... + beta_p = v_1 - v_0. 
+    # The additional constraint beta_0 = v_0 is dealt via offset   
     est_m[[n_iter]] <- solver(A_temp, b_temp, constraint = v1 - v0)          #  (p x K)
 
     # Covariance calculation would fail in the first iteration
@@ -116,7 +116,7 @@ ginv <- function (X, tol = sqrt(.Machine$double.eps)) {
   } else if (!any(Positive)) {
     array(0, dim(X)[2L:1L])
   } else {
-    Xsvd$v[, Positive, drop = FALSE] %*%
+    Xsvd$v[, Positive, drop = FALSE] %*% 
       ((1 / Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))
   }
 }
@@ -126,11 +126,11 @@ get_vz <- function(X, bg, Z, object, pred_fun, feature_names, w, ...) {
   m <- nrow(Z)
   not_Z <- !Z
   n_bg <- nrow(bg) / m   # because bg was replicated m times
-
+  
   # Replicate not_Z, so that X, bg, not_Z are all of dimension (m*n_bg x p)
   g <- rep(seq_len(m), each = n_bg)
   not_Z <- not_Z[g, , drop = FALSE]
-
+  
   if (is.matrix(X)) {
     # Remember that columns of X and bg are perfectly aligned in this case
     X[not_Z] <- bg[not_Z]
@@ -143,7 +143,7 @@ get_vz <- function(X, bg, Z, object, pred_fun, feature_names, w, ...) {
     }
   }
   preds <- check_pred(pred_fun(object, X, ...), n = nrow(X))
-
+  
   # Aggregate
   if (is.null(w)) {
     return(rowsum(preds, group = g, reorder = FALSE) / n_bg)
@@ -162,15 +162,15 @@ weighted_colMeans <- function(x, w = NULL, ...) {
     if (nrow(x) != length(w)) {
       stop("Weights w not compatible with matrix x")
     }
-    out <- colSums(x * w, ...) / sum(w)
+    out <- colSums(x * w, ...) / sum(w)  
   }
   matrix(out, nrow = 1L)
 }
 
 # Binds list of matrices along new first axis
 abind1 <- function(a) {
   out <- array(
-    dim = c(length(a), dim(a[[1L]])),
+    dim = c(length(a), dim(a[[1L]])), 
     dimnames = c(list(NULL), dimnames(a[[1L]]))
   )
   for (i in seq_along(a)) {
@@ -196,9 +196,9 @@ reorganize_list <- function(alist, nms) {
 # Checks and reshapes predictions to (n x K) matrix
 check_pred <- function(x, n) {
   if (
-    !is.vector(x) &&
-    !is.matrix(x) &&
-    !is.data.frame(x) &&
+    !is.vector(x) && 
+    !is.matrix(x) && 
+    !is.data.frame(x) && 
     !(is.array(x) && length(dim(x)) <= 2L)
   ) {
     stop("Predictions must be a vector, matrix, data.frame, or <=2D array")
@@ -235,7 +235,7 @@ summarize_strategy <- function(p, exact, deg) {
   }
   if (deg == 0L) {
     return("Kernel SHAP values by iterative sampling")
-  }
+  } 
   paste("Kernel SHAP values by the hybrid strategy of degree", deg)
 }