{"id":1391,"date":"2024-01-07T10:37:24","date_gmt":"2024-01-07T09:37:24","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1391"},"modified":"2024-01-19T09:33:51","modified_gmt":"2024-01-19T08:33:51","slug":"explain-that-tidymodels-blackbox","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/01\/07\/explain-that-tidymodels-blackbox\/","title":{"rendered":"Explain that tidymodels blackbox!"},"content":{"rendered":"\n

Let’s explain a {tidymodels} random forest by classic explainability methods (permutation importance, partial dependence plots (PDP), Friedman’s H statistics), and also fancy SHAP.<\/p>\n\n\n\n

Disclaimer: {hstats}, {kernelshap} and {shapviz} are three of my own packages.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Diabetes data<\/h2>\n\n\n\n

We will use the diabetes prediction dataset of Kaggle<\/a> to model diabetes (yes\/no) as a function of six demographic features (age, gender, BMI, hypertension, heart disease, and smoking history). It has 100k rows.<\/p>\n\n\n\n

Note: The data additionally contains the typical diabetes indicators HbA1c level and blood glucose level, but we wont use them to avoid potential causality issues, and to gain insights also for people that do not know these values.<\/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
# https:\/\/www.kaggle.com\/datasets\/iammustafatz\/diabetes-prediction-dataset\n\nlibrary(tidyverse)\nlibrary(tidymodels)\nlibrary(hstats)\nlibrary(kernelshap)\nlibrary(shapviz)\nlibrary(patchwork)\n\ndf0 <- read.csv(\"diabetes_prediction_dataset.csv\")  # from above Kaggle link\ndim(df0)  # 100000 9\nhead(df0)\n# gender age hypertension heart_disease smoking_history   bmi HbA1c_level blood_glucose_level diabetes\n# Female  80            0             1           never 25.19         6.6                 140        0\n# Female  54            0             0         No Info 27.32         6.6                  80        0\n#   Male  28            0             0           never 27.32         5.7                 158        0\n# Female  36            0             0         current 23.45         5.0                 155        0\n#   Male  76            1             1         current 20.14         4.8                 155        0\n# Female  20            0             0           never 27.32         6.6                  85        0\n\nsummary(df0)\nanyNA(df0)  # FALSE\ntable(df0$smoking_history, useNA = \"ifany\")\n\n# DATA PREPARATION\n\n# Note: tidymodels needs a factor response for classification\ndf1 <- df0 |>\n  transform(\n    y = factor(diabetes, levels = 0:1, labels = c(\"No\", \"Yes\")),\n    female = (gender == \"Female\") * 1,\n    smoking_history = factor(\n      smoking_history, \n      levels = c(\"No Info\", \"never\", \"former\", \"not current\", \"current\", \"ever\")\n    ),\n    bmi = pmin(bmi, 50)\n  )\n\n# UNIVARIATE ANALYSIS\n\nggplot(df1, aes(diabetes)) +\n  geom_bar(fill = \"chartreuse4\")\n\ndf1  |>  \n  select(age, bmi, HbA1c_level, blood_glucose_level) |> \n  pivot_longer(everything()) |> \n  ggplot(aes(value)) +\n  geom_histogram(fill = \"chartreuse4\", bins = 19) +\n  facet_wrap(~ name, scale = \"free_x\")\n\nggplot(df1, aes(smoking_history)) +\n  geom_bar(fill = \"chartreuse4\")\n\ndf1 |> \n  select(heart_disease, hypertension, female) |>\n  pivot_longer(everything()) |> \n  ggplot(aes(name, value)) +\n  stat_summary(fun = mean, geom = \"bar\", fill = \"chartreuse4\") +\n  xlab(element_blank())\n<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
\"\"
“yes” proportion of binary variables (including the response)<\/figcaption><\/figure>\n\n\n\n
\"\"
Distribution of numeric variables<\/figcaption><\/figure>\n\n\n\n
\"\"
Distribution of smoking_history<\/figcaption><\/figure>\n\n\n\n
Modeling<\/h2>\n\n\n\n
Let’s fit a random forest via tidymodels with {ranger} backend. <\/p>\n\n\n\n
We add a predict function pf()<\/code> that outputs only the probability of the “Yes” class.<\/p>\n\n\n\n
set.seed(1)\nix <- initial_split(df1, strata = diabetes, prop = 0.8)\ntrain <- training(ix)\ntest <- testing(ix)\n\nxvars <- c("age", "bmi", "smoking_history", "heart_disease", "hypertension", "female")\n\nrf_spec <- rand_forest(trees = 500) |> \n  set_mode("classification") |> \n  set_engine("ranger", num.threads = NULL, seed = 49)\n\nrf_wf <- workflow() |> \n  add_model(rf_spec) |>\n  add_formula(reformulate(xvars, "y"))\n\nmodel <- rf_wf |> \n    fit(train)\n\n# predict() gives No\/Yes columns\npredict(model, head(test), type = "prob")\n# .pred_No .pred_Yes\n#    0.981    0.0185\n\n# We need to extract only the "Yes" probabilities\npf <- function(m, X) {\n  predict(m, X, type = "prob")$.pred_Yes\n}\npf(model, head(test))  # 0.01854290 ...\n<\/pre><\/div>\n\n\n\n
Classic explanation methods<\/h2>\n\n\n\n
# 4 times repeated permutation importance wrt test logloss\nimp <- perm_importance(\n  model, X = test, y = "diabetes", v = xvars, pred_fun = pf, loss = "logloss"\n)\nplot(imp) +\n  xlab("Increase in test logloss")\n\n# Partial dependence of age\npartial_dep(model, v = "age", train, pred_fun = pf) |> \n  plot()\n\n# All PDP in one patchwork\np <- lapply(xvars, function(x) plot(partial_dep(model, v = x, X = train, pred_fun = pf)))\nwrap_plots(p) &\n  ylim(0, 0.23) &\n  ylab("Probability")\n\n# Friedman's H stats\nsystem.time( # 20 s\n  H <- hstats(model, train[xvars], approx = TRUE, pred_fun = pf)\n)\nH  # 15% of prediction variability comes from interactions\nplot(H)\n\n# Stratified PDP of strongest interaction\npartial_dep(model, "age", BY = "bmi", X = train, pred_fun = pf) |> \n  plot(show_points = FALSE)<\/pre><\/div>\n\n\n\n
Feature importance<\/h3>\n\n\n\n
Permutation importance<\/em> measures by how much the average test loss (in our case log loss) increases when a feature is shuffled before calculating the losses. We repeat the process four times and also show standard errors.<\/p>\n\n\n\n
\"\"
Permutation importance: Age and BMI are the two main risk factors.<\/figcaption><\/figure>\n\n\n\n
Main effects<\/h3>\n\n\n\n
Main effects are estimated by PDP. They show how the average prediction changes with a feature, keeping every other feature fixed. Using a fixed vertical axis helps to grasp the strenght of the effect.<\/p>\n\n\n\n
\"\"
PDPs: The diabetes risk tends to increase with age, high (and very low) BMI, presence of heart disease\/hypertension, and it is a bit lower for females and non-smoker.<\/figcaption><\/figure>\n\n\n\n
Interaction strength<\/h3>\n\n\n\n
Interaction strength can be measured by Friedman’s H statistics, see the earlier blog post<\/a>. A specific interaction can then be visualized by a stratified PDP.<\/p>\n\n\n\n
\"\"
Friedman’s H statistics: Left: BMI and age are the two features with clearly strongest interactions. Right: Their pairwise interaction explains about 10% of their joint effect variability.<\/figcaption><\/figure>\n\n\n\n
\"\"
Stratified PDP: The strong interaction between age and BMI is clearly visible. A high BMI makes the age effect on diabetes stronger.<\/figcaption><\/figure>\n\n\n\n
SHAP<\/h2>\n\n\n\n
What insights does a SHAP analysis bring? <\/p>\n\n\n\n
We will crunch slow exact permutation SHAP values via kernelshap::permshap()<\/code>. If we had more features, we could switch to <\/p>\n\n\n\n