{"id":1734,"date":"2024-06-28T16:00:41","date_gmt":"2024-06-28T14:00:41","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1734"},"modified":"2024-07-04T17:38:36","modified_gmt":"2024-07-04T15:38:36","slug":"shap-values-of-additive-models","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/06\/28\/shap-values-of-additive-models\/","title":{"rendered":"SHAP Values of Additive Models"},"content":{"rendered":"\n

Within only a few years, SHAP (Shapley additive explanations) has emerged as the number 1 way to investigate black-box models. The basic idea is to decompose model predictions into additive contributions of the features in a fair way. Studying decompositions of many predictions allows to derive global properties of the model.<\/p>\n\n\n\n

What happens if we apply SHAP algorithms to additive models? Why would this ever make sense?<\/strong><\/p>\n\n\n\n

In the spirit of our “Lost In Translation” series, we provide both high-quality Python and R code.<\/p>\n\n\n\n

The models<\/h2>\n\n\n\n

Let’s build the models using a dataset with three highly correlated covariates and a (deterministic) response.<\/p>\n\n\n

\n\t\t\t
\n\t\t\t\t
\n\t\t\t\t
R<\/div>\n\t\t\t<\/div>
\n\t\t\t\t
Python<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
\n\n
library(lightgbm)\nlibrary(kernelshap)\nlibrary(shapviz)\n\n#===================================================================\n# Make small data\n#===================================================================\n\nmake_data <- function(n = 100) {\n  x1 <- seq(0.01, 1, length = n)\n  data.frame(\n    x1 = x1,\n    x2 = log(x1),\n    x3 = x1 > 0.7\n  ) |>\n    transform(y = 1 + 0.2 * x1 + 0.5 * x2 + x3 + sin(2 * pi * x1))\n}\ndf <- make_data()\nhead(df)\ncor(df) |>\n  round(2)\n\n#      x1   x2   x3    y\n# x1 1.00 0.90 0.80 0.46\n# x2 0.90 1.00 0.58 0.58\n# x3 0.80 0.58 1.00 0.51\n# y  0.46 0.58 0.51 1.00\n\n#===================================================================\n# Additive linear model and additive boosted trees\n#===================================================================\n\n# Linear regression\nfit_lm <- lm(y ~ poly(x1, 3) + poly(x2, 3) + x3, data = df)\nsummary(fit_lm)\n\n# Boosted trees\nxvars <- setdiff(colnames(df), \"y\")\nX <- data.matrix(df[xvars])\n\nparams <- list(\n  learning_rate = 0.05,\n  objective = \"mse\",\n  max_depth = 1,\n  colsample_bynode = 0.7\n)\n\nfit_lgb <- lgb.train(\n  params = params,\n  data = lgb.Dataset(X, label = df$y),\n  nrounds = 300\n)<\/pre><\/div>\n\n<\/div>
\n\n
import numpy as np\nimport lightgbm as lgb\nimport shap\nfrom sklearn.preprocessing import PolynomialFeatures\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.linear_model import LinearRegression\n\n#===================================================================\n# Make small data\n#===================================================================\n\ndef make_data(n=100):\n    x1 = np.linspace(0.01, 1, n)\n    x2 = np.log(x1)\n    x3 = x1 > 0.7\n    X = np.column_stack((x1, x2, x3))\n\n    y = 1 + 0.2 * x1 + 0.5 * x2 + x3 + np.sin(2 * np.pi * x1)\n    \n    return X, y\n\nX, y = make_data()\n\n#===================================================================\n# Additive linear model and additive boosted trees\n#===================================================================\n\n# Linear model with polynomial terms\npoly = PolynomialFeatures(degree=3, include_bias=False)\n\npreprocessor =  ColumnTransformer(\n    transformers=[\n        (\"poly0\", poly, [0]),\n        (\"poly1\", poly, [1]),\n        (\"other\", \"passthrough\", [2]),\n    ]\n)\n\nmodel_lm = Pipeline(\n    steps=[\n        (\"preprocessor\", preprocessor),\n        (\"lm\", LinearRegression()),\n    ]\n)\n_ = model_lm.fit(X, y)\n\n# Boosted trees with single-split trees\nparams = dict(\n    learning_rate=0.05,\n    objective=\"mse\",\n    max_depth=1,\n    colsample_bynode=0.7,\n)\n\nmodel_lgb = lgb.train(\n    params=params,\n    train_set=lgb.Dataset(X, label=y),\n    num_boost_round=300,\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
SHAP<\/h2>\n\n\n\n
For both models, we use exact permutation SHAP and exact Kernel SHAP. Furthermore, the linear model is analyzed with “additive SHAP”, and the tree-based model with TreeSHAP.<\/p>\n\n\n\n
Do the algorithms provide the same?<\/p>\n\n\n
\n\t\t\t
\n\t\t\t\t
\n\t\t\t\t
R<\/div>\n\t\t\t<\/div>
\n\t\t\t\t
Python<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
\n\n
system.time({  # 1s\n  shap_lm <- list(\n    add = shapviz(additive_shap(fit_lm, df)),\n    kern = kernelshap(fit_lm, X = df[xvars], bg_X = df),\n    perm = permshap(fit_lm, X = df[xvars], bg_X = df)\n  )\n\n  shap_lgb <- list(\n    tree = shapviz(fit_lgb, X),\n    kern = kernelshap(fit_lgb, X = X, bg_X = X),\n    perm = permshap(fit_lgb, X = X, bg_X = X)\n  )\n})\n\n# Consistent SHAP values for linear regression\nall.equal(shap_lm$add$S, shap_lm$perm$S)\nall.equal(shap_lm$kern$S, shap_lm$perm$S)\n\n# Consistent SHAP values for boosted trees\nall.equal(shap_lgb$lgb_tree$S, shap_lgb$lgb_perm$S)\nall.equal(shap_lgb$lgb_kern$S, shap_lgb$lgb_perm$S)\n\n# Linear coefficient of x3 equals slope of SHAP values\ntail(coef(fit_lm), 1)                # 1.112096\ndiff(range(shap_lm$kern$S[, \"x3\"]))  # 1.112096\n\nsv_dependence(shap_lm$add, xvars)sv_dependence(shap_lm$add, xvars, color_var = NULL)<\/pre><\/div>\n\n<\/div>
\n\n
shap_lm = {\n    \"add\": shap.Explainer(model_lm.predict, masker=X, algorithm=\"additive\")(X),\n    \"perm\": shap.Explainer(model_lm.predict, masker=X, algorithm=\"exact\")(X),\n    \"kern\": shap.KernelExplainer(model_lm.predict, data=X).shap_values(X),\n}\n\nshap_lgb = {\n    \"tree\": shap.Explainer(model_lgb)(X),\n    \"perm\": shap.Explainer(model_lgb.predict, masker=X, algorithm=\"exact\")(X),\n    \"kern\": shap.KernelExplainer(model_lgb.predict, data=X).shap_values(X),\n}\n\n# Consistency for additive linear regression\neps = 1e-12\nassert np.abs(shap_lm[\"add\"].values - shap_lm[\"perm\"].values).max() < eps\nassert np.abs(shap_lm[\"perm\"].values - shap_lm[\"kern\"]).max() < eps\n\n# Consistency for additive boosted trees\nassert np.abs(shap_lgb[\"tree\"].values - shap_lgb[\"perm\"].values).max() < eps\nassert np.abs(shap_lgb[\"perm\"].values - shap_lgb[\"kern\"]).max() < eps\n\n# Linear effect of last feature in the fitted model\nmodel_lm.named_steps[\"lm\"].coef_[-1]  # 1.112096\n\n# Linear effect of last feature derived from SHAP values (ignore the sign)\nshap_lm[\"perm\"][:, 2].values.ptp()    # 1.112096\n\nshap.plots.scatter(shap_lm[\"add\"])<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
\"\"
SHAP dependence plot of the additive linear model and the additive explainer (Python).<\/figcaption><\/figure>\n\n\n\n
Yes – the three algorithms within model provide the same SHAP values. Furthermore, the SHAP values reconstruct the additive components of the features. <\/p>\n\n\n\n
Didactically, this is very helpful when introducing SHAP as a method: Pick a white-box and a black-box model and compare their SHAP dependence plots. For the white-box model, you simply see the additive components, while the dependence plots of the black-box model show scatter due to interactions.<\/p>\n\n\n\n
Remark: The exact equivalence between algorithms is lost, when<\/strong><\/p>\n\n\n\n