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?
 A: 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.
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?

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
#

