Kernel SHAP<\/h2>\n\n\n\n
SHAP is one of the most used model interpretation technique in Machine Learning. It decomposes predictions into additive contributions of the features in a fair way. For tree-based methods, the fast TreeSHAP algorithm exists. For general models, one has to resort to computationally expensive Monte-Carlo sampling or the faster Kernel SHAP algorithm. Kernel SHAP uses a regression trick to get the SHAP values of an observation with a comparably small number of calls to the predict function of the model. Still, it is much slower than TreeSHAP.<\/p>\n\n\n\n
Two good references for Kernel SHAP:<\/p>\n\n\n\n
- Scott M. Lundberg and Su-In Lee. A Unified Approach to Interpreting Model Predictions. Advances in Neural Information Processing Systems 30, 2017.<\/li>
- Ian Covert and Su-In Lee. Improving KernelSHAP: Practical Shapley Value Estimation Using Linear Regression. Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, PMLR 130:3457-3465, 2021.<\/li><\/ol>\n\n\n\n
In our last post<\/a>, we introduced our new “kernelshap” package in R. Since then, the package has been substantially improved, also by the big help of David Watson<\/a>: <\/p>\n\n\n\n
- The package now supports multi-dimensional predictions.<\/li>
- It received a massive speed-up<\/li>
- Additionally, parallel computing can be activated for even faster calculations.<\/li>
- The interface has become more intuitive.<\/li>
- If the number of features is small (up to ten or eleven), it can provide exact Kernel SHAP values just like the reference Python implementation<\/a>.<\/li>
- For a larger number of features, it now uses partly-exact (“hybrid”) calculations, very similar to the logic in the Python implementation.<\/li><\/ol>\n\n\n\n
With those changes, the R implementation is about to meet the Python version at eye level.<\/p>\n\n\n\n
Example with four features <\/h2>\n\n\n\n
In the following, we use the diamonds data to fit a linear regression with<\/p>\n\n\n\n
- log(price) as response<\/li>
- log(carat) as numeric feature<\/li>
- clarity, color and cut as categorical features (internally dummy encoded)<\/li>
- interactions between log(carat) and the other three “C” variables. Note that the interactions are very weak<\/li><\/ul>\n\n\n\n
Then, we calculate SHAP decompositions for about 1000 diamonds (every 53th diamond), using 120 diamonds as background dataset. In this case, both R and Python will use exact calculations based on m=2^4 – 2 = 14 possible binary on-off vectors (a value of 1 representing a feature value picked from the original observation, a value of 0 a value picked from the background data).<\/p>\n\n\n
\n\t\t\t\n\t\t\t\t\n\t\t\t\tR<\/div>\n\t\t\t<\/div>\n\t\t\t\tPython<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t\n\nlibrary(ggplot2)\nlibrary(kernelshap)\n\n# Turn ordinal factors into unordered\nord <- c(\"clarity\", \"color\", \"cut\")\ndiamonds[, ord] <- lapply(diamonds[ord], factor, ordered = FALSE)\n\n# Fit model\nfit <- lm(log(price) ~ log(carat) * (clarity + color + cut), data = diamonds)\n\n# Subset of 120 diamonds used as background data\nbg_X <- diamonds[seq(1, nrow(diamonds), 450), ]\n\n# Subset of 1018 diamonds to explain\nX_small <- diamonds[seq(1, nrow(diamonds), 53), c(\"carat\", ord)]\n\n# Exact KernelSHAP (5 seconds)\nsystem.time(\n ks <- kernelshap(fit, X_small, bg_X = bg_X) \n)\nks\n\n# SHAP values of first 2 observations:\n# carat clarity color cut\n# [1,] -2.050074 -0.28048747 0.1281222 0.01587382\n# [2,] -2.085838 0.04050415 0.1283010 0.03731644\n\n# Using parallel backend\nlibrary(\"doFuture\")\n\nregisterDoFuture()\nplan(multisession, workers = 2) # Windows\n# plan(multicore, workers = 2) # Linux, macOS, Solaris\n\n# 3 seconds on second call\nsystem.time(\n ks3 <- kernelshap(fit, X_small, bg_X = bg_X, parallel = TRUE) \n)\n\n# Visualization\nlibrary(shapviz)\n\nsv <- shapviz(ks)\nsv_importance(sv, \"bee\")<\/pre><\/div>\n\n<\/div>
\n\nimport numpy as np\nimport pandas as pd\nfrom plotnine.data import diamonds\nfrom statsmodels.formula.api import ols\nfrom shap import KernelExplainer\n\n# Turn categoricals into integers because, inconveniently, kernel SHAP\n# requires numpy array as input\nord = [\"clarity\", \"color\", \"cut\"]\nx = [\"carat\"] + ord\ndiamonds[ord] = diamonds[ord].apply(lambda x: x.cat.codes)\nX = diamonds[x].to_numpy()\n\n# Fit model with interactions and dummy variables\nfit = ols(\n \"np.log(price) ~ np.log(carat) * (C(clarity) + C(cut) + C(color))\", \n data=diamonds\n).fit()\n\n# Background data (120 rows)\nbg_X = X[0:len(X):450]\n\n# Define subset of 1018 diamonds to explain\nX_small = X[0:len(X):53]\n\n# Calculate KernelSHAP values\nks = KernelExplainer(\n model=lambda X: fit.predict(pd.DataFrame(X, columns=x)), \n data = bg_X\n)\nsv = ks.shap_values(X_small) # 74 seconds\nsv[0:2]\n\n# array([[-2.05007406, -0.28048747, 0.12812216, 0.01587382],\n# [-2.0858379 , 0.04050415, 0.12830103, 0.03731644]])<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2022\/09\/imp.jpeg\")
SHAP summary plot (R model)<\/figcaption><\/figure>\n\n\n\n The results match, hurray!<\/strong><\/p>\n\n\n\n
Example with nine features<\/h2>\n\n\n\n
The computation effort of running exact<\/strong> Kernel SHAP explodes with the number of features. For nine features, the number of relevant on-off vectors is 2^9 – 2 = 510, i.e. about 36 times larger than with four features.<\/p>\n\n\n\n
We now modify above example, adding five additional features to the model. Note that the model structure is completely non-sensical. We just use it to get a feeling about what impact a 36 times larger workload has.<\/em><\/p>\n\n\n\n
Besides exact calculations, we use an almost exact hybrid approach for both R and Python, using 126 on-off vectors (p*(p+1) for the exact part and 4p for the sampling part, where p is the number of features), resulting in a significant speed-up both in R and Python.<\/p>\n\n\n
\n\t\t\t\n\t\t\t\t\n\t\t\t\tR<\/div>\n\t\t\t<\/div>\n\t\t\t\tPython<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t\n\nfit <- lm(\n log(price) ~ log(carat) * (clarity + color + cut) + x + y + z + table + depth, \n data = diamonds\n)\n\n# Subset of 1018 diamonds to explain\nX_small <- diamonds[seq(1, nrow(diamonds), 53), setdiff(names(diamonds), \"price\")]\n\n# Exact Kernel SHAP: 61 seconds\nsystem.time(\n ks <- kernelshap(fit, X_small, bg_X = bg_X, exact = TRUE) \n)\nks\n# carat cut color clarity depth table x y z\n# [1,] -1.842799 0.01424231 0.1266108 -0.27033874 -0.0007084443 0.0017787647 -0.1720782 0.001330275 -0.006445693\n# [2,] -1.876709 0.03856957 0.1266546 0.03932912 -0.0004202636 -0.0004871776 -0.1739880 0.001397792 -0.006560624\n\n# Default, using an almost exact hybrid algorithm: 17 seconds\nsystem.time(\n ks <- kernelshap(fit, X_small, bg_X = bg_X, parallel = TRUE) \n)\n# carat cut color clarity depth table x y z\n# [1,] -1.842799 0.01424231 0.1266108 -0.27033874 -0.0007084443 0.0017787647 -0.1720782 0.001330275 -0.006445693\n# [2,] -1.876709 0.03856957 0.1266546 0.03932912 -0.0004202636 -0.0004871776 -0.1739880 0.001397792 -0.006560624<\/pre><\/div>\n\n<\/div>
\n\nx = [\"carat\"] + ord + [\"table\", \"depth\", \"x\", \"y\", \"z\"]\nX = diamonds[x].to_numpy()\n\n# Fit model with interactions and dummy variables\nfit = ols(\n \"np.log(price) ~ np.log(carat) * (C(clarity) + C(cut) + C(color)) + table + depth + x + y + z\", \n data=diamonds\n).fit()\n\n# Background data (120 rows)\nbg_X = X[0:len(X):450]\n\n# Define subset of 1018 diamonds to explain\nX_small = X[0:len(X):53]\n\n# Calculate KernelSHAP values: 12 minutes\nks = KernelExplainer(\n model=lambda X: fit.predict(pd.DataFrame(X, columns=x)), \n data = bg_X\n)\nsv = ks.shap_values(X_small)\nsv[0:2]\n# array([[-1.84279897e+00, -2.70338744e-01, 1.26610769e-01,\n# 1.42423108e-02, 1.77876470e-03, -7.08444295e-04,\n# -1.72078182e-01, 1.33027467e-03, -6.44569296e-03],\n# [-1.87670887e+00, 3.93291219e-02, 1.26654599e-01,\n# 3.85695742e-02, -4.87177593e-04, -4.20263565e-04,\n# -1.73988040e-01, 1.39779179e-03, -6.56062359e-03]])\n\n# Now, using a hybrid between exact and sampling: 5 minutes\nsv = ks.shap_values(X_small, nsamples=126)\nsv[0:2]\n# array([[-1.84279897e+00, -2.70338744e-01, 1.26610769e-01,\n# 1.42423108e-02, 1.77876470e-03, -7.08444295e-04,\n# -1.72078182e-01, 1.33027467e-03, -6.44569296e-03],\n# [-1.87670887e+00, 3.93291219e-02, 1.26654599e-01,\n# 3.85695742e-02, -4.87177593e-04, -4.20263565e-04,\n# -1.73988040e-01, 1.39779179e-03, -6.56062359e-03]])<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
Again, the results are essentially the same between R and Python, but also between the hybrid algorithm and the exact algorithm. This is interesting, because the hybrid algorithm is significantly faster than the exact one.<\/p>\n\n\n\n
Wrap-Up<\/h2>\n\n\n\n
- R is catching up with Python’s superb “shap” package.<\/li>
- For two non-trivial linear regressions with interactions, the “kernelshap” package in R provides the same output as Python.<\/li>
- The hybrid between exact and sampling KernelSHAP (as implemented in Python and R) offers a very good trade-off between speed and accuracy.<\/li>
kernelshap()<\/code>in R is fast!<\/li><\/ul>\n\n\n\nThe Python and R codes can be found here:<\/p>\n\n\n\n
- R co<\/a>de<\/a><\/li>
- Python code<\/a><\/li><\/ul>\n\n\n\n
The examples were run on a Windows notebook with an Intel i7-8650U 4 core CPU.<\/p>\n","protected":false},"excerpt":{"rendered":"
“R <-> Python” continued… Kernel SHAP<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[16,17,9],"tags":[10,6,5],"class_list":["post-952","post","type-post","status-publish","format-standard","hentry","category-machine-learning","category-programming","category-statistics","tag-lost-in-translation","tag-python","tag-r"],"featured_image_src":null,"author_info":{"display_name":"Michael Mayer","author_link":"https:\/\/lorentzen.ch\/index.php\/author\/michael\/"},"_links":{"self":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/952","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/comments?post=952"}],"version-history":[{"count":8,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/952\/revisions"}],"predecessor-version":[{"id":963,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/952\/revisions\/963"}],"wp:attachment":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/media?parent=952"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/categories?post=952"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/tags?post=952"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
- Python code<\/a><\/li><\/ul>\n\n\n\n
- R co<\/a>de<\/a><\/li>
- For a larger number of features, it now uses partly-exact (“hybrid”) calculations, very similar to the logic in the Python implementation.<\/li><\/ol>\n\n\n\n