# How to calculate SHAP for a factor in a linear model?

Shapley additive explanation (SHAP) are used to explain the prediction of a model $$Y = f(X_1,...,X_p) + \varepsilon$$. If we observe $$x_1,...,x_p$$ and predict $$y$$, then for each $$i$$ the contribution of $$x_i$$ to the prediction $$y$$ is defined by

where $$S$$ is a subset of indices from the set $$F=\{1,...,p\}$$ and $$X_S = \{X_i, i\in S\}$$. The weights $$w$$ sum up to 1 and have an explicit formula (not important here I guess). The paper "A Unified Approach to Interpreting Model Predictions" of Lundberg and Lee is the reference here.

In the case of a linear model $$f(X_1,...,X_p) = \sum \beta_i X_i$$ where no correlations among the covariates are assumed, then the previous formula simplifies to $$\phi_i = \beta_i (x_i - \mathbb{E}(X_i))$$

How do I translate that into categorical variables? Assume we have X_1 a categorical variable with two levels "A" and "B". Assume that we observe "A". In principle, if I base the calculation on the last displayed equation, then the contribution has to be $$\phi_1 = \beta_1(1-m/n)$$ where $$m/n$$ is the proportion of appearance of "A" in the training set (I estimate the expectation from the training set). I used R packages "shapr" and "iml", and none of them gives that result. Any idea?

• In a blog post, I have demonstrated how to achieve it. As far as I know, the SHAP values will exactly correspond to the fitted coefficients, for additive predictors: lorentzen.ch/index.php/2022/12/21/… Commented Dec 31, 2022 at 14:29

## 2 Answers

It is actually true, but R packages "iml" and "shapr" don't have exact implementations of SHAP for linear models. I checked the calculation with the python package "shap" of the authors of the paper "A Unified Approach to Interpreting Model Predictions" of Lundberg and Lee. They have a separate function "LinearExplainer" The SHAP value corresponds actually to $$\phi_1 (1-m/n)$$ with one extra detail. If you have one categorical variable with two levels A and B, then its hot-encoding is 0-1. The SHAP value for level A corresponds to $$\phi_1 (1-m0/n)$$ where $$m0/n$$ is the proportion of zeros and the SHAP value for level B corresponds to $$\phi_1 (1-(n-m0)/n)$$. For more complicated scenarios, we have to use a hot encoding then finally for each categorical variable sum up everything.

In principle, if I base the calculation on the last displayed equation, then the contribution has to be $$\phi_1 = \beta_1(1-m/n)$$ where $$m/n$$ is the proportion of appearance of "A" in the training set (I estimate the expectation from the training set). I used R packages "shapr" and "iml", and none of them gives that result. Any idea?

The iml package doesn't make an exact computation but uses an estimate based on a Monte Carlo simulation.

The code below demonstrates this. If we sample $$10^5$$ times, then estimate is better in comparison to the standard $$10^2$$ times, and the estimated value $$-1.848816$$ is very close to the computed value $$-1.85$$.

library(iml)
library(datasets)
set.seed(1)

### some dataset with a factor to be used as example
data = datasets::ToothGrowth

### the linear model
mod_lm <- lm(len ~ ., data = data)

### predicting and computing phi effect contributions with shapley from the iml package
mod_pr <- Predictor$$new(mod_lm, data = data[,-1]) # ,-1 removes the len column which is no predictor x.interest = data[1,-1] shapley <- Shapley$$new(mod_pr, x.interest = x.interest, sample.size = 10^5)

### results
shapley
# Head of results:
# feature       phi   phi.var feature.value
# 1    supp -1.848816  3.422533       supp=VC
# 2    dose -6.489162 37.185843      dose=0.5
(1-mean(data['supp'] == 'VC'))*mod_lm\$coefficients['suppVC']
# suppVC
# -1.85
#