{"id":772,"date":"2022-04-25T09:54:41","date_gmt":"2022-04-25T07:54:41","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=772"},"modified":"2022-04-25T09:54:42","modified_gmt":"2022-04-25T07:54:42","slug":"let-the-flashlight-shine-with-plotly","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2022\/04\/25\/let-the-flashlight-shine-with-plotly\/","title":{"rendered":"Let the flashlight shine with plotly"},"content":{"rendered":"\n

There are different R packages devoted to model agnostic interpretability, DALEX<\/a> and iml<\/a> being among the best known. In 2019, I added flashlight<\/a> <\/p>\n\n\n\n

\"logo.png\"<\/a><\/figure>\n\n\n\n

for a couple of reasons:<\/p>\n\n\n\n

  1. Its explainers work with case weights.<\/li>
  2. Multiple explainers can be combined to a multi-explainer.<\/li>
  3. Stratified calculation is possible.<\/li><\/ol>\n\n\n\n

    Since almost all plots in flashlight are constructed with ggplot<\/code>, it is super easy to turn them into interactive plotly objects: just add a simple ggplotly()<\/code> to the end of the call. <\/p>\n\n\n\n

    However… it is not straightforward to show interactive plots in a blog! Thus, we show only screenshots of the resulting plots here and refer to the complete HTML report here: https:\/\/mayer79.github.io\/flashlight_plotly\/flashlight_plotly.html<\/a><\/p>\n\n\n\n

    We will use a sweet dataset with more than 20\u2019000 houses to model house prices by a set of derived features such as the logarithmic living area. The location will be represented by the postal code.<\/p>\n\n\n\n

    Data preparation<\/h3>\n\n\n\n

    We first load the data and prepare some of the columns for modeling. Furthermore, we specify the set of features and the 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><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    library(dplyr)\nlibrary(flashlight)\nlibrary(plotly)\nlibrary(ranger)\nlibrary(lme4)\nlibrary(moderndive)\nlibrary(splitTools)\nlibrary(MetricsWeighted)\n\nset.seed(4933)\n\ndata(\"house_prices\")\n\nprep <- house_prices %>% \n  mutate(\n    log_price = log(price),\n    log_sqft_living = log(sqft_living),\n    log_sqft_lot = log(sqft_lot),\n    log_sqft_basement = log1p(sqft_basement),\n    year = as.numeric(format(date, '%Y')),\n    age = year - yr_built\n  )\n\nx <- c(\n  \"year\", \"age\", \"log_sqft_living\", \"log_sqft_lot\", \n  \"bedrooms\", \"bathrooms\", \"log_sqft_basement\", \n  \"condition\", \"waterfront\", \"zipcode\"\n)\n\ny <- \"log_price\"\n\nhead(prep[c(y, x)])\n\n## # A tibble: 6 x 11\n##   log_price  year   age log_sqft_living log_sqft_lot bedrooms bathrooms\n##       <dbl> <dbl> <dbl>           <dbl>        <dbl>    <int>     <dbl>\n## 1      12.3  2014    59            7.07         8.64        3      1   \n## 2      13.2  2014    63            7.85         8.89        3      2.25\n## 3      12.1  2015    82            6.65         9.21        2      1   \n## 4      13.3  2014    49            7.58         8.52        4      3   \n## 5      13.1  2015    28            7.43         9.00        3      2   \n## 6      14.0  2014    13            8.60        11.5         4      4.5 \n## # ... with 4 more variables: log_sqft_basement <dbl>, condition <fct>,\n## #   waterfront <lgl>, zipcode <fct><\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    Train \/ test split<\/h3>\n\n\n\n
    Then, we split the dataset into 80% training and 20% test rows, stratified on the (binned) response log_price<\/code>.<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    idx <- partition(prep[[y]], c(train = 0.8, test = 0.2), type = \"stratified\")\n\ntrain <- prep[idx$train, ]\ntest <- prep[idx$test, ]<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    Models<\/h3>\n\n\n\n
    We fit two models:<\/p>\n\n\n\n
    1. A linear mixed model with random postal code effect.<\/li>
    2. A random forest with 500 trees.<\/li><\/ol>\n\n\n
      \n\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t
      R<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      # Mixed-effects model\nfit_lmer <- lmer(\n  update(reformulate(x, \"log_price\"), . ~ . - zipcode + (1 | zipcode)),\n  data = train\n)\n\n# Random forest\nfit_rf <- ranger(\n  reformulate(x, \"log_price\"),\n  always.split.variables = \"zipcode\",\n  data = train\n)\ncat(\"R-squared OOB:\", fit_rf$r.squared)\n## R-squared OOB: 0.8463311<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      Model inspection<\/h3>\n\n\n\n
      Now, we are ready to inspect our two models regarding performance, variable importance, and effects.<\/p>\n\n\n\n
      Set up explainers<\/h4>\n\n\n\n
      First, we pack all model dependent information into flashlights (the explainer objects) and combine them to a multiflashlight. As evaluation dataset, we pass the test data. This ensures that interpretability tools using the response (e.g., performance measures and permutation importance) are not being biased by overfitting.<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      fl_lmer <- flashlight(model = fit_lmer, label = \"LMER\")\nfl_rf <- flashlight(\n  model = fit_rf,\n  label = \"RF\",\n  predict_function = function(mod, X) predict(mod, X)$predictions\n)\nfls <- multiflashlight(\n  list(fl_lmer, fl_rf),\n  y = \"log_price\",\n  data = test,\n  metrics = list(RMSE = rmse, `R-squared` = r_squared)\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      Model performance<\/h4>\n\n\n\n
      Let\u2019s evaluate model RMSE and R-squared on the hold-out dataset. Here, the mixed-effects model performs a tiny little bit better than the random forest:<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      (light_performance(fls) %>%\n  plot(fill = \"darkred\") +\n    labs(title = \"Model performance\", x = element_blank())) %>%\n  ggplotly()<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      \"\"
      Model performance (png)<\/figcaption><\/figure>\n\n\n\n
      Permutation importance<\/h4>\n\n\n\n
      Next, we inspect the variable strength based on permutation importance. It shows by how much the RMSE is being increased when shuffling a variable before prediction. The results are quite similar between the two models.<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      (light_importance(fls, v = x) %>%\n    plot(fill = \"darkred\") +\n    labs(title = \"Permutation importance\", y = \"Drop in RMSE\")) %>%\n  ggplotly()<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      \"\"
      Variable importance (png)<\/figcaption><\/figure>\n\n\n\n
      ICE plot<\/h4>\n\n\n\n
      To get an impression of the effect of the living area, we select 200 observations and profile their predictions with increasing (log) living area, keeping everything else fixed (Ceteris Paribus). These ICE (individual conditional expectation) plots are vertically centered in order to highlight potential interaction effects. If all curves coincide, there are no interaction effects and we can say that the effect of the feature is modelled in an additive way (no surprise for the additive linear mixed-effects model).<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      (light_ice(fls, v = \"log_sqft_living\", n_max = 200, center = \"middle\") %>%\n    plot(alpha = 0.05, color = \"darkred\") +\n    labs(title = \"Centered ICE plot\", y = \"log_price (shifted)\")) %>%\n  ggplotly()<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      \"\"<\/figure>\n\n\n\n
      Partial dependence plots<\/h4>\n\n\n\n
      Averaging many uncentered ICE curves provides the famous partial dependence plot, introduced in Friedman\u2019s seminal paper on gradient boosting machines (2001).<\/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><\/div>\n\t\t\t<\/div>\n\t\t\t
      \n\n
      (light_profile(fls, v = \"log_sqft_living\", n_bins = 21) %>%\n    plot(rotate_x = FALSE) +\n    labs(title = \"Partial dependence plot\", y = y) +\n    scale_colour_viridis_d(begin = 0.2, end = 0.8)) %>%\n  ggplotly()<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
      \"\"
      Partial dependence plots (png)<\/figcaption><\/figure>\n\n\n\n
      Multiple effects visualized together<\/h4>\n\n\n\n
      The last figure extends the partial dependence plot with three additional curves, all evaluated on the hold-out dataset:<\/p>\n\n\n\n