Files
shap--shap/shap/explainers/_linear.py
T
2026-07-13 13:22:52 +08:00

516 lines
22 KiB
Python

import warnings
from collections.abc import Callable
from typing import Any, Literal
import numpy as np
import numpy.typing as npt
import pandas as pd
from scipy.sparse import issparse
from tqdm.auto import tqdm
from .. import links, maskers
from ..utils._exceptions import (
DimensionError,
InvalidFeaturePerturbationError,
InvalidModelError,
)
from ._explainer import Explainer
class LinearExplainer(Explainer):
"""Computes SHAP values for a linear model, optionally accounting for inter-feature correlations.
This computes the SHAP values for a linear model and can account for the
correlations among the input features. Assuming features are independent
leads to interventional SHAP values which for a linear model are ``coef[i] *
(x[i] - X.mean(0)[i])`` for the ith feature. If instead we account for
correlations, then we prevent any problems arising from collinearity and
share credit among correlated features. Accounting for correlations can be
computationally challenging, but ``LinearExplainer`` uses sampling to
estimate a transform that can then be applied to explain any prediction of
the model.
Parameters
----------
model : (coef, intercept) or sklearn.linear_model.*
User supplied linear model either as either a parameter pair or sklearn object.
masker : function, numpy.array, pandas.DataFrame, tuple of (mean, cov), shap.maskers.Masker
A callable Python object used to "mask" out hidden features of the form
``masker(binary_mask, x)``. It takes a single input sample and a binary
mask and returns a matrix of masked samples. These masked samples are
evaluated using the model function and the outputs are then averaged.
As a shortcut for the standard masking using by SHAP you can pass a
background data matrix instead of a function and that matrix will be
used for masking.
You can also provide a tuple of ``(mean, covariance)``, or pass in a
masker meant for tabular data (i.e., :class:`.maskers.Independent`,
:class:`.maskers.Impute`, or :class:`.maskers.Partition`) directly.
data : (mean, cov), numpy.array, pandas.DataFrame, iml.DenseData or scipy.csr_matrix
The background dataset to use for computing conditional expectations.
Note that only the mean and covariance of the dataset are used. This
means passing a raw data matrix is just a convenient alternative to
passing the mean and covariance directly.
nsamples : int
Number of samples to use when estimating the transformation matrix used
to account for feature correlations.
feature_perturbation : None (default), "interventional" or "correlation_dependent"
DEPRECATED: this option is now deprecated in favor of using the appropriate
tabular masker and will be removed in a future release.
There are two ways we might want to compute SHAP values, either the full
conditional SHAP values or the interventional SHAP values.
- For interventional SHAP values we break any dependence structure between
features in the model and so uncover how the model would behave if we
intervened and changed some of the inputs. This approach is used by the
Independent and Partition maskers.
- For the full conditional SHAP values we respect the correlations among the
input features, so if the model depends on one input but that input is
correlated with another input, then both get some credit for the model's
behavior. This approach is used by the Impute masker.
The interventional option stays "true to the model" meaning it will only
give credit to features that are actually used by the model, while the
correlation option stays "true to the data" in the sense that it only
considers how the model would behave when respecting the correlations in
the input data. For sparse case only interventional option is supported.
Examples
--------
See `Linear explainer examples <https://shap.readthedocs.io/en/latest/api_examples/explainers/LinearExplainer.html>`_
"""
feature_perturbation: Literal["interventional", "correlation_dependent"]
nsamples: int
coef: Any
intercept: Any
mean: npt.NDArray[np.floating[Any]]
cov: npt.NDArray[np.floating[Any]]
expected_value: float | npt.NDArray[np.floating[Any]]
M: int
valid_inds: npt.NDArray[np.intp]
avg_proj: npt.NDArray[np.floating[Any]]
mean_transformed: npt.NDArray[np.floating[Any]]
x_transform: npt.NDArray[np.floating[Any]]
def __init__(
self,
model: Any,
masker: Any,
link: Callable[[Any], Any] = links.identity,
nsamples: int = 1000,
feature_perturbation: None | Literal["interventional", "correlation_dependent"] = None,
**kwargs: Any,
) -> None:
if "feature_dependence" in kwargs:
emsg = "The option feature_dependence has been renamed to feature_perturbation!"
raise ValueError(emsg)
if feature_perturbation is not None: # pragma: no cover
wmsg = (
"The feature_perturbation option is now deprecated in favor of using the appropriate "
"masker (maskers.Independent, maskers.Partition or maskers.Impute)."
)
warnings.warn(wmsg, FutureWarning)
else:
feature_perturbation = "interventional"
if feature_perturbation not in ("interventional", "correlation_dependent"):
emsg = "feature_perturbation must be one of 'interventional' or 'correlation_dependent'"
raise InvalidFeaturePerturbationError(emsg)
self.feature_perturbation = feature_perturbation
# wrap the incoming masker object as a shap.Masker object before calling
# parent class constructor, which does the same but without respecting
# the user-provided feature_perturbation choice
if isinstance(masker, pd.DataFrame) or (
(isinstance(masker, np.ndarray) or issparse(masker)) and len(masker.shape) == 2
):
if self.feature_perturbation == "correlation_dependent":
masker = maskers.Impute(masker)
else:
masker = maskers.Independent(masker)
elif issubclass(type(masker), tuple) and len(masker) == 2:
if self.feature_perturbation == "correlation_dependent":
masker = maskers.Impute({"mean": masker[0], "cov": masker[1]}, method="linear")
else:
masker = maskers.Independent({"mean": masker[0], "cov": masker[1]})
super().__init__(model, masker, link=link, **kwargs)
self.nsamples = nsamples
# extract what we need from the given model object
self.coef, self.intercept = LinearExplainer._parse_model(model)
# extract the data
if issubclass(type(self.masker), (maskers.Independent, maskers.Partition)):
self.feature_perturbation = "interventional"
elif issubclass(type(self.masker), maskers.Impute):
self.feature_perturbation = "correlation_dependent"
else:
raise NotImplementedError(
"The Linear explainer only supports the Independent, Partition, and Impute maskers right now!"
)
data = getattr(self.masker, "data", None)
# convert DataFrame's to numpy arrays
if isinstance(data, pd.DataFrame):
data = data.values
# get the mean and covariance of the model
if getattr(self.masker, "mean", None) is not None:
self.mean = self.masker.mean
self.cov = self.masker.cov
elif isinstance(data, dict) and len(data) == 2:
self.mean = data["mean"]
if isinstance(self.mean, pd.Series):
self.mean = self.mean.values
self.cov = data["cov"]
if isinstance(self.cov, pd.DataFrame):
self.cov = self.cov.values
elif isinstance(data, tuple) and len(data) == 2:
self.mean = data[0]
if isinstance(self.mean, pd.Series):
self.mean = self.mean.values
self.cov = data[1]
if isinstance(self.cov, pd.DataFrame):
self.cov = self.cov.values
elif data is None:
raise ValueError("A background data distribution must be provided!")
else:
if issparse(data):
self.mean = np.array(np.mean(data, 0))[0]
if self.feature_perturbation != "interventional":
raise NotImplementedError(
"Only feature_perturbation = 'interventional' is supported for sparse data"
)
else:
self.mean = np.array(np.mean(data, 0)).flatten() # assumes it is an array
if self.feature_perturbation == "correlation_dependent":
data_shape = np.shape(data)
if len(data_shape) > 1:
n_samples = data_shape[0]
n_features = data_shape[1]
if n_samples <= n_features:
warnings.warn(
f"The number of samples ({n_samples}) is less than or equal to "
f"the number of features ({n_features}). This will produce a "
"singular covariance matrix and may result in unreliable SHAP "
"values when using feature_perturbation='correlation_dependent'.",
UserWarning,
stacklevel=2,
)
self.cov = np.cov(data, rowvar=False)
# print(self.coef, self.mean.flatten(), self.intercept)
# Note: mean can be numpy.matrixlib.defmatrix.matrix or numpy.matrix type depending on numpy version
if issparse(self.mean) or str(type(self.mean)).endswith("matrix'>"):
# accept both sparse and dense coef
# if not issparse(self.coef):
# self.coef = np.asmatrix(self.coef)
self.expected_value = np.dot(self.coef, self.mean) + self.intercept
# unwrap the matrix form
if len(self.expected_value) == 1: # type: ignore[arg-type]
self.expected_value = self.expected_value[0, 0] # type: ignore[index]
else:
self.expected_value = np.array(self.expected_value)[0]
else:
self.expected_value = np.dot(self.coef, self.mean) + self.intercept
self.M = len(self.mean)
# if needed, estimate the transform matrices
if self.feature_perturbation == "correlation_dependent":
self.valid_inds = np.where(np.diag(self.cov) > 1e-8)[0]
self.mean = self.mean[self.valid_inds]
self.cov = self.cov[:, self.valid_inds][self.valid_inds, :]
self.coef = self.coef[self.valid_inds]
# group perfectly redundant variables together
self.avg_proj, sum_proj = duplicate_components(self.cov)
self.cov = np.matmul(np.matmul(self.avg_proj, self.cov), self.avg_proj.T)
self.mean = np.matmul(self.avg_proj, self.mean)
self.coef = np.matmul(sum_proj, self.coef)
# if we still have some multi-collinearity present then we just add regularization...
e, _ = np.linalg.eig(self.cov)
if e.min() < 1e-7:
self.cov = self.cov + np.eye(self.cov.shape[0]) * 1e-6 # type: ignore[assignment]
mean_transform, x_transform = self._estimate_transforms(nsamples)
self.mean_transformed = np.matmul(mean_transform, self.mean)
self.x_transform = x_transform
elif self.feature_perturbation == "interventional":
if nsamples != 1000:
warnings.warn("Setting nsamples has no effect when feature_perturbation = 'interventional'!")
else:
raise InvalidFeaturePerturbationError(
"Unknown type of feature_perturbation provided: " + self.feature_perturbation
)
def _estimate_transforms(
self, nsamples: int
) -> tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]:
"""Uses block matrix inversion identities to quickly estimate transforms.
After a bit of matrix math we can isolate a transform matrix (# features x # features)
that is independent of any sample we are explaining. It is the result of averaging over
all feature permutations, but we just use a fixed number of samples to estimate the value.
TODO: Do a brute force enumeration when # feature subsets is less than nsamples. This could
happen through a recursive method that uses the same block matrix inversion as below.
"""
M = len(self.coef)
mean_transform = np.zeros((M, M))
x_transform = np.zeros((M, M))
inds = np.arange(M, dtype=int)
for _ in tqdm(range(nsamples), "Estimating transforms"):
np.random.shuffle(inds)
cov_inv_SiSi = np.zeros((0, 0))
cov_Si = np.zeros((M, 0))
for j in range(M):
i = inds[j]
# use the last Si as the new S
cov_S = cov_Si
cov_inv_SS = cov_inv_SiSi
# get the new cov_Si
cov_Si = self.cov[:, inds[: j + 1]] # type: ignore[assignment]
# compute the new cov_inv_SiSi from cov_inv_SS
d = cov_Si[i, :-1].T
t = np.matmul(cov_inv_SS, d)
Z = self.cov[i, i]
u = Z - np.matmul(t.T, d)
cov_inv_SiSi = np.zeros((j + 1, j + 1))
if j > 0:
cov_inv_SiSi[:-1, :-1] = cov_inv_SS + np.outer(t, t) / u
cov_inv_SiSi[:-1, -1] = cov_inv_SiSi[-1, :-1] = -t / u
cov_inv_SiSi[-1, -1] = 1 / u
# + coef @ (Q(bar(Sui)) - Q(bar(S)))
mean_transform[i, i] += self.coef[i]
# + coef @ R(Sui)
coef_R_Si = np.matmul(self.coef[inds[j + 1 :]], np.matmul(cov_Si, cov_inv_SiSi)[inds[j + 1 :]])
mean_transform[i, inds[: j + 1]] += coef_R_Si
# - coef @ R(S)
coef_R_S = np.matmul(self.coef[inds[j:]], np.matmul(cov_S, cov_inv_SS)[inds[j:]])
mean_transform[i, inds[:j]] -= coef_R_S
# - coef @ (Q(Sui) - Q(S))
x_transform[i, i] += self.coef[i]
# + coef @ R(Sui)
x_transform[i, inds[: j + 1]] += coef_R_Si
# - coef @ R(S)
x_transform[i, inds[:j]] -= coef_R_S
mean_transform /= nsamples
x_transform /= nsamples
return mean_transform, x_transform
@staticmethod
def _parse_model(model: Any) -> tuple[Any, Any]:
"""Attempt to pull out the coefficients and intercept from the given model object."""
# raw coefficients
if isinstance(model, tuple) and len(model) == 2:
coef = model[0]
intercept = model[1]
# sklearn style model
elif hasattr(model, "coef_") and hasattr(model, "intercept_"):
# work around for multi-class with a single class
if len(model.coef_.shape) > 1 and model.coef_.shape[0] == 1:
coef = model.coef_[0]
try:
intercept = model.intercept_[0]
except TypeError:
intercept = model.intercept_
else:
coef = model.coef_
intercept = model.intercept_
else:
raise InvalidModelError("An unknown model type was passed: " + str(type(model)))
return coef, intercept
@staticmethod
def supports_model_with_masker(model: Any, masker: Any) -> bool:
"""Determines if we can parse the given model."""
if not isinstance(masker, (maskers.Independent, maskers.Partition, maskers.Impute)):
return False
try:
LinearExplainer._parse_model(model)
except Exception:
return False
return True
def explain_row(
self,
*row_args: Any,
max_evals: int | Literal["auto"],
main_effects: bool,
error_bounds: bool,
outputs: Any,
silent: bool,
**kwargs: Any,
) -> dict[str, Any]:
"""Explains a single row and returns the tuple (row_values, row_expected_values, row_mask_shapes)."""
assert len(row_args) == 1, "Only single-argument functions are supported by the Linear explainer!"
X = row_args[0]
if len(X.shape) == 1:
X = X.reshape(1, -1)
# convert dataframes
if isinstance(X, (pd.Series, pd.DataFrame)):
X = X.values
if len(X.shape) not in (1, 2):
raise DimensionError(f"Instance must have 1 or 2 dimensions! Not: {len(X.shape)}")
if self.feature_perturbation == "correlation_dependent":
if issparse(X):
raise InvalidFeaturePerturbationError(
"Only feature_perturbation = 'interventional' is supported for sparse data"
)
phi = (
np.matmul(np.matmul(X[:, self.valid_inds], self.avg_proj.T), self.x_transform.T) - self.mean_transformed
)
phi = np.matmul(phi, self.avg_proj)
full_phi = np.zeros((phi.shape[0], self.M))
full_phi[:, self.valid_inds] = phi
phi = full_phi
elif self.feature_perturbation == "interventional":
if issparse(X):
phi = np.array(np.multiply(X - self.mean, self.coef))
# if len(self.coef.shape) == 1:
# return np.array(np.multiply(X - self.mean, self.coef))
# else:
# return [np.array(np.multiply(X - self.mean, self.coef[i])) for i in range(self.coef.shape[0])]
else:
phi = np.array(X - self.mean) * self.coef
# if len(self.coef.shape) == 1:
# phi = np.array(X - self.mean) * self.coef
# return np.array(X - self.mean) * self.coef
# else:
# return [np.array(X - self.mean) * self.coef[i] for i in range(self.coef.shape[0])]
return {
"values": phi.T,
"expected_values": self.expected_value,
"mask_shapes": (X.shape[1:],),
"main_effects": phi.T,
"clustering": None,
}
def shap_values(self, X: npt.NDArray[np.floating[Any]] | pd.DataFrame | pd.Series) -> npt.NDArray[np.floating[Any]]:
"""Estimate the SHAP values for a set of samples.
Parameters
----------
X : numpy.array, pandas.DataFrame or scipy.csr_matrix
A matrix of samples (# samples x # features) on which to explain the model's output.
Returns
-------
np.array
Estimated SHAP values, usually of shape ``(# samples x # features)``.
Each row sums to the difference between the model output for that
sample and the expected value of the model output (which is stored
as the ``expected_value`` attribute of the explainer).
The shape of the returned array depends on the number of model outputs:
* one output: array of shape ``(#num_samples, *X.shape[1:])``.
* multiple outputs: array of shape ``(#num_samples, *X.shape[1:],
#num_outputs)``.
.. versionchanged:: 0.45.0
Return type for models with multiple outputs changed from list to np.ndarray.
"""
# convert dataframes
if isinstance(X, (pd.Series, pd.DataFrame)):
X = X.values
# assert isinstance(X, np.ndarray), "Unknown instance type: " + str(type(X))
if len(X.shape) not in (1, 2):
raise DimensionError(f"Instance must have 1 or 2 dimensions! Not: {len(X.shape)}")
if self.feature_perturbation == "correlation_dependent":
if issparse(X):
raise InvalidFeaturePerturbationError(
"Only feature_perturbation = 'interventional' is supported for sparse data"
)
phi = (
np.matmul(np.matmul(X[:, self.valid_inds], self.avg_proj.T), self.x_transform.T) - self.mean_transformed
)
phi = np.matmul(phi, self.avg_proj)
full_phi = np.zeros((phi.shape[0], self.M))
full_phi[:, self.valid_inds] = phi
return full_phi
elif self.feature_perturbation == "interventional":
if issparse(X):
if len(self.coef.shape) == 1:
return np.array(np.multiply(X - self.mean, self.coef))
else:
return np.stack(
[np.array(np.multiply(X - self.mean, self.coef[i])) for i in range(self.coef.shape[0])], axis=-1
)
else:
if len(self.coef.shape) == 1:
return np.array(X - self.mean) * self.coef
else:
return np.stack(
[np.array(X - self.mean) * self.coef[i] for i in range(self.coef.shape[0])], axis=-1
)
def duplicate_components(
C: npt.NDArray[np.floating[Any]],
) -> tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]:
D = np.diag(1 / np.sqrt(np.diag(C)))
C = np.matmul(np.matmul(D, C), D)
components = -np.ones(C.shape[0], dtype=int)
count = -1
for i in range(C.shape[0]):
found_group = False
for j in range(C.shape[0]):
if components[j] < 0 and np.abs(2 * C[i, j] - C[i, i] - C[j, j]) < 1e-8:
if not found_group:
count += 1
found_group = True
components[j] = count
proj = np.zeros((len(np.unique(components)), C.shape[0]))
proj[0, 0] = 1
for i in range(1, C.shape[0]):
proj[components[i], i] = 1
return (proj.T / proj.sum(1)).T, proj