Some of the new features of SHAPforxgboost<\/strong> <\/p>\n\n\n\n An interesting alternative to calculate and plot SHAP values for different tree-based models is the treeshap<\/a> package by Szymon Maksymiuk et al. Keep an eye on this one – it is actively being developed!<\/p>\n\n\n\n A couple of years ago, the concept of Shapely values from game theory from the 1950ies was discovered e.g. by Scott Lundberg as an interesting approach to explain predictions of ML models. <\/p>\n\n\n\n The basic idea is to decompose a prediction in a fair way into additive contributions of features. Repeating the process for many predictions provides a brilliant way to investigate the model as a whole. <\/p>\n\n\n\n The main resource on the topic is Scott Lundberg’s site<\/a>. Besides this, I’d recomment to go through these two fantastic blog posts, even if you already know what SHAP values are:<\/p>\n\n\n\n As an example, we will try to model log house prices of 20’000 sold houses in Kings County. The dataset is available e.g. on OpenML.org<\/a> under ID 42092.<\/p>\n\n\n\n We start by downloading the data and preparing it for modelling.<\/p>\n\n\n Next, we fit a manually tuned XGBoost model to the data.<\/p>\n\n\n The resulting model consists of about 600 trees and reaches a validation RMSE of 0.16. This means that about 2\/3 of the predictions are within 16% of the observed price, using the empirical rule<\/a>.<\/p>\n\n\n\n ML models are rarely of any use without interpreting its results, so let’s use SHAP to peak into the model.<\/p>\n\n\n\n The analysis includes a first plot with SHAP importances. Then, with decreasing importance, dependence plots are shown to get an impression on the effects of each feature.<\/p>\n\n\n Some of the plots are shown below. The code actually produces all plots, see the corresponding html<\/a> output on github.<\/p>\n\n\n\n Let’s try out the Note: LightGBM Version 3.2.1 on CRAN is not working properly under Windows. This will be fixed in the next release of LightGBM. As a temporary solution, you need to build it from the current master branch<\/a><\/em>.<\/p>\n\n\n Early stopping on the validation data selects about 900 trees as being optimal and results in a validation RMSE of also 0.16.<\/p>\n\n\n\n We use exactly the same short snippet to analyze the model by SHAP.<\/p>\n\n\n Again, we only show some of the output and refer to the html<\/a> of the corresponding rmarkdown. Overall, the model seems to be very similar to the one obtained by XGBoost.<\/p>\n\n\n\nshap.plot.dependence<\/code>() has received the option to select the heuristically strongest interacting feature on the color scale, see last section for details.<\/li>shap.plot.dependence<\/code>() now allows jitter and alpha transparency.<\/li>shap.importance<\/code>() returns SHAP importances without plotting them.<\/li>![\"logo.png\"\/](\"https:\/\/github.com\/liuyanguu\/SHAPforxgboost\/blob\/master\/man\/figures\/logo.png?raw=true\"){kind=logo}
<\/figure><\/div>\n\n\n\nWhat is SHAP?<\/h2>\n\n\n\n
Illustration<\/h2>\n\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/image-1024x224.png\")
Fetch and prepare data<\/h3>\n\n\n\n
library(farff)\nlibrary(OpenML)\nlibrary(dplyr)\nlibrary(xgboost)\nlibrary(ggplot2)\nlibrary(SHAPforxgboost)\n\n# Load King Country house prices dataset on OpenML\n# ID 42092, https:\/\/www.openml.org\/d\/42092\ndf <- getOMLDataSet(data.id = 42092)$data\nhead(df)\n\n# Prepare\ndf <- df %>%\n mutate(\n log_price = log(price),\n log_sqft_lot = log(sqft_lot),\n year = as.numeric(substr(date, 1, 4)),\n building_age = year - yr_built,\n zipcode = as.integer(as.character(zipcode))\n )\n\n# Define response and features\ny <- \"log_price\"\nx <- c(\"grade\", \"year\", \"building_age\", \"sqft_living\",\n \"log_sqft_lot\", \"bedrooms\", \"bathrooms\", \"floors\", \"zipcode\",\n \"lat\", \"long\", \"condition\", \"waterfront\")\n\n# random split\nset.seed(83454)\nix <- sample(nrow(df), 0.8 * nrow(df))<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
Fit XGBoost model<\/h3>\n\n\n\n
dtrain <- xgb.DMatrix(data.matrix(df[ix, x]),\n label = df[ix, y])\ndvalid <- xgb.DMatrix(data.matrix(df[-ix, x]),\n label = df[-ix, y])\n\nparams <- list(\n objective = \"reg:squarederror\",\n learning_rate = 0.05,\n subsample = 0.9,\n colsample_bynode = 1,\n reg_lambda = 2,\n max_depth = 5\n)\n\nfit_xgb <- xgb.train(\n params,\n data = dtrain,\n watchlist = list(valid = dvalid),\n early_stopping_rounds = 20,\n print_every_n = 100,\n nrounds = 10000 # early stopping\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
Compact SHAP analysis<\/h4>\n\n\n\n
# Step 1: Select some observations\nX <- data.matrix(df[sample(nrow(df), 1000), x])\n\n# Step 2: Crunch SHAP values\nshap <- shap.prep(fit_xgb, X_train = X)\n\n# Step 3: SHAP importance\nshap.plot.summary(shap)\n\n# Step 4: Loop over dependence plots in decreasing importance\nfor (v in shap.importance(shap, names_only = TRUE)) {\n p <- shap.plot.dependence(shap, v, color_feature = \"auto\", \n alpha = 0.5, jitter_width = 0.1) +\n ggtitle(v)\n print(p)\n}<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/xgbshap-1-1024x731.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/xgbshap-3-1-1024x731.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/xgbshap-8-1024x731.png\")
Same workflow for LightGBM<\/h3>\n\n\n\n
SHAPforxgboost<\/code> package with LightGBM. <\/p>\n\n\n\nlibrary(lightgbm)\n\ndtrain <- lgb.Dataset(data.matrix(df[ix, x]),\n label = df[ix, y])\ndvalid <- lgb.Dataset(data.matrix(df[-ix, x]),\n label = df[-ix, y])\n\nparams <- list(\n objective = \"regression\",\n learning_rate = 0.05,\n subsample = 0.9,\n reg_lambda = 2,\n num_leaves = 15\n)\n\nfit_lgb <- lgb.train(\n params,\n data = dtrain,\n valids = list(valid = dvalid),\n early_stopping_rounds = 20,\n eval_freq = 100,\n eval = \"rmse\",\n nrounds = 10000\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
SHAP analysis<\/h4>\n\n\n\n
X <- data.matrix(df[sample(nrow(df), 1000), x])\nshap <- shap.prep(fit_lgb, X_train = X)\nshap.plot.summary(shap)\n\nfor (v in shap.importance(shap, names_only = TRUE)) {\n p <- shap.plot.dependence(shap, v, color_feature = \"auto\", \n alpha = 0.5, jitter_width = 0.1) +\n ggtitle(v)\n print(p)\n}<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/lgbfit-1-1024x731.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/lgbfit-3-1024x731.png\")
How does the dependence plot selects the color variable?<\/h2>\n\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/interaction-2-1024x731.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/interaction-4-1024x731.png\")