This post is heavily based on the new {shapviz} vignette<\/a>.<\/p>\n\n\n\n Besides other features, a model with geographic components contains features like<\/p>\n\n\n\n Like any feature, the effect of a single geographic feature can be described using SHAP dependence<\/strong> plots. However, studying the effect of latitude (or any other location dependent feature) alone is often not very illuminating – simply due to strong interaction effects and correlations with other geographic features.<\/p>\n\n\n\n That’s where the additivity of SHAP<\/strong> values comes into play: The sum of SHAP values of all geographic components represent the total geographic effect, and this sum can be visualized as a heatmap or 3D scatterplot against latitude\/longitude (or any other geographic representation).<\/p>\n\n\n\n For illustration, we will use a beautiful house price dataset containing information on about 14’000 houses sold in 2016 in Miami-Dade County. Some of the columns are as follows:<\/p>\n\n\n\n (Italic features are geographic components.) For more background on this dataset, see Mayer et al [2].<\/a><\/p>\n\n\n\n We will fit an XGBoost model to explain log(price) as a function of lat\/long, size, and quality\/age.<\/p>\n\n\n Let’s first study selected SHAP dependence plots, evaluated on the validation dataset with around 2800 observations. Note that we could as well use the training data for this purpose, but it is a bit large.<\/p>\n\n\n And now the two-dimensional plot of the sum of SHAP values:<\/p>\n\n\n The last plot gives a good impression on price levels, but note:<\/p>\n\n\n\n We will now change above model in two ways, not unlike the model in Mayer et al [2].<\/a><\/p>\n\n\n\nSetting<\/h2>\n\n\n\n
\n
A first example <\/h2>\n\n\n\n
\n
devtools::install_github(\"ModelOriented\/shapviz\", dependencies = TRUE)\nlibrary(xgboost)\nlibrary(ggplot2)\nlibrary(shapviz) # Needs development version 0.9.0 from github\n\nhead(miami)\n\nx_coord <- c(\"LATITUDE\", \"LONGITUDE\")\nx_nongeo <- c(\"TOT_LVG_AREA\", \"LND_SQFOOT\", \"structure_quality\", \"age\")\nx <- c(x_coord, x_nongeo)\n\n# Train\/valid split\nset.seed(1)\nix <- sample(nrow(miami), 0.8 * nrow(miami))\nX_train <- data.matrix(miami[ix, x])\nX_valid <- data.matrix(miami[-ix, x])\ny_train <- log(miami$SALE_PRC[ix])\ny_valid <- log(miami$SALE_PRC[-ix])\n\n# Fit XGBoost model with early stopping\ndtrain <- xgb.DMatrix(X_train, label = y_train)\ndvalid <- xgb.DMatrix(X_valid, label = y_valid)\n\nparams <- list(learning_rate = 0.2, objective = \"reg:squarederror\", max_depth = 5)\n\nfit <- xgb.train(\n params = params, \n data = dtrain, \n watchlist = list(valid = dvalid), \n early_stopping_rounds = 20,\n nrounds = 1000,\n callbacks = list(cb.print.evaluation(period = 100))\n)<\/pre><\/div>\n\n<\/div>
%load_ext lab_black\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.datasets import fetch_openml\n\ndf = fetch_openml(data_id=43093, as_frame=True)\nX, y = df.data, np.log(df.target)\nX.head()\n\n# Data split and model\nfrom sklearn.model_selection import train_test_split\nimport xgboost as xgb\n\nx_coord = [\"LONGITUDE\", \"LATITUDE\"]\nx_nongeo = [\"TOT_LVG_AREA\", \"LND_SQFOOT\", \"structure_quality\", \"age\"]\nx = x_coord + x_nongeo\n\nX_train, X_valid, y_train, y_valid = train_test_split(\n X[x], y, test_size=0.2, random_state=30\n)\n\n# Fit XGBoost model with early stopping\ndtrain = xgb.DMatrix(X_train, label=y_train)\ndvalid = xgb.DMatrix(X_valid, label=y_valid)\n\nparams = dict(learning_rate=0.2, objective=\"reg:squarederror\", max_depth=5)\n\nfit = xgb.train(\n params=params,\n dtrain=dtrain,\n evals=[(dvalid, \"valid\")],\n verbose_eval=100,\n early_stopping_rounds=20,\n num_boost_round=1000,\n)\n<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
SHAP dependence plots<\/h3>\n\n\n\n
sv <- shapviz(fit, X_pred = X_valid)\nsv_dependence(\n sv, \n v = c(\"TOT_LVG_AREA\", \"structure_quality\", \"LONGITUDE\", \"LATITUDE\"), \n alpha = 0.2\n)<\/pre><\/div>\n\n<\/div>
import shap\n\nxgb_explainer = shap.Explainer(fit)\nshap_values = xgb_explainer(X_valid)\n\nv = [\"TOT_LVG_AREA\", \"structure_quality\", \"LONGITUDE\", \"LATITUDE\"]\nshap.plots.scatter(shap_values[:, v], color=shap_values[:, v])<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2023\/06\/image-1024x352.png\")
Total coordindate effect<\/h3>\n\n\n\n
sv_dependence2D(sv, x = \"LONGITUDE\", y = \"LATITUDE\") +\n coord_equal()<\/pre><\/div>\n\n<\/div>
shap_coord = shap_values[:, x_coord]\nplt.scatter(*list(shap_coord.data.T), c=shap_coord.values.sum(axis=1), s=4)\nax = plt.gca()\nax.set_aspect(\"equal\", adjustable=\"box\")\nplt.colorbar()\nplt.title(\"Total location effect\")\nplt.show()<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2023\/06\/image-1.png\")
\n
Two modifications<\/h2>\n\n\n\n
\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2023\/06\/image-2.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2023\/06\/image-3.png\")