Integrate the two types of weights

Author: mayer79

NEWS.md

@@ -3,9 +3,12 @@
 ### Major bug fix
 
 `kernelshap()` used a wrong weighting logic, leading to values slightly off. This has
-  been fixed with the help of Prof Mario Wuethrich of ETHZ and 
-  [Ian Covert's wonderful Github repo](https://github.com/iancovert/shapley-regression).
-  Now, exact Kernel SHAP returns identical values as exact permutation SHAP.
+been fixed with the help of Prof Mario Wuethrich of ETHZ and Ian Covert and his
+[wonderful Github repo](https://github.com/iancovert/shapley-regression).
+Now, exact Kernel SHAP returns identical values as exact permutation SHAP.
+All variants of `kernelshap()` had been affected by this (exact, sampling, hybrid).
+For models with interactions up to two, the bug had no consequences - which is why
+it went unnoticed.
 
 Fixed in [#168](https://github.com/ModelOriented/kernelshap/pull/168).
 

R/utils_kernelshap.R

@@ -107,7 +107,7 @@ solver <- function(A, b, constraint) {
 # to Kernel SHAP weights -> (m x p) matrix.
 # The argument S can be used to restrict the range of sum(z).
 sample_Z <- function(p, m, feature_names, S = 1:(p - 1L)) {
-  probs <- kernel_weights_per_coalition_size(p, S = S)
+  probs <- kernel_weights(p, per_coalition_size = TRUE, S = S)
   N <- S[sample.int(length(S), m, replace = TRUE, prob = probs)]
 
   # Then, conditional on that number, set random positions of z to 1
@@ -158,7 +158,7 @@ input_sampling <- function(p, m, deg, feature_names) {
 input_exact <- function(p, feature_names) {
   Z <- exact_Z(p, feature_names = feature_names)
   Z <- Z[2L:(nrow(Z) - 1L), , drop = FALSE]
-  kw <- kernel_weights(p) # Kernel weights for all subsets
+  kw <- kernel_weights(p, per_coalition_size = FALSE) # Kernel weights for all subsets
   w <- kw[rowSums(Z)] # Corresponding weight for each row in Z
   w <- w / sum(w)
   list(Z = Z, w = w, A = crossprod(Z, w * Z))
@@ -203,7 +203,7 @@ input_partly_exact <- function(p, deg, feature_names) {
     stop("p must be >=2*deg")
   }
 
-  kw <- kernel_weights(p)
+  kw <- kernel_weights(p, per_coalition_size = FALSE)
 
   Z <- vector("list", deg)
   for (k in seq_len(deg)) {
@@ -216,15 +216,16 @@ input_partly_exact <- function(p, deg, feature_names) {
   list(Z = Z, w = w, A = crossprod(Z, w * Z))
 }
 
-# Kernel weights
-kernel_weights <- function(p, S = seq_len(p - 1L)) {
-  probs <- (p - 1L) / (choose(p, S) * S * (p - S))
-  return(probs / sum(probs))
-}
-
-# Kernel weights per coalition size
-kernel_weights_per_coalition_size <- function(p, S = seq_len(p - 1L)) {
-  probs <- 1 / (S * (p - S))
+# Kernel weight distribution
+#
+# `per_coalition_size = TRUE` is required, e.g., when one wants to sample random masks
+# according to the Kernel SHAP distribution: Pick a coalition size as per
+# these weights, then randomly place "on" positions. `FALSE` refer to weights
+# if all masks has been calculated and one wants to calculate their weights based
+# on the number of "on" positions.
+kernel_weights <- function(p, per_coalition_size, S = seq_len(p - 1L)) {
+  const <- if (per_coalition_size) 1 else choose(p, S)
+  probs <- (p - 1) / (const * S * (p - S)) # could drop the numerator
   return(probs / sum(probs))
 }
 
@@ -234,7 +235,7 @@ prop_exact <- function(p, deg) {
   if (deg == 0) {
     return(0)
   }
-  w <- kernel_weights_per_coalition_size(p)
+  w <- kernel_weights(p, per_coalition_size = TRUE)
   w_total <- 2 * sum(w[seq_len(deg)]) - w[deg] * (p == 2 * deg)
   return(w_total)
 }

tests/testthat/test-kernelshap-utils.R

@@ -1,11 +1,13 @@
 test_that("sum of kernel weights is 1", {
   for (p in 2:10) {
-    expect_equal(sum(kernel_weights(p)), 1.0)
+    expect_equal(sum(kernel_weights(p, per_coalition_size = FALSE)), 1.0)
+    expect_equal(sum(kernel_weights(p, per_coalition_size = TRUE)), 1.0)
   }
 })
 
 test_that("Sum of kernel weights is 1, even for subset of domain", {
-  expect_equal(sum(kernel_weights(10L, S = 2:5)), 1.0)
+  expect_equal(sum(kernel_weights(10L, S = 2:5, per_coalition_size = FALSE)), 1.0)
+  expect_equal(sum(kernel_weights(10L, S = 2:5, per_coalition_size = TRUE)), 1.0)
 })
 
 p <- 10L
@@ -105,6 +107,17 @@ test_that("hybrid weights sum to 1 for different p and degree 2", {
   }
 })
 
+test_that("sampling input A is comparable from exact input", {
+  set.seed(1)
+
+  for (p in 2:6) {
+    feature_names <- LETTERS[1:p]
+    pa <- input_exact(p, feature_names)
+    sa <- input_sampling(p, m = 100000L, deg = 0, feature_names = feature_names)
+    expect_true(all(abs(pa$A - sa$A) < 0.01))
+  }
+})
+
 test_that("partly_exact_Z(p, k) fails for bad p or k", {
   expect_error(partly_exact_Z(0L, k = 1L, feature_names = LETTERS[1:p]))
   expect_error(partly_exact_Z(5L, k = 3L, feature_names = LETTERS[1:p]))