-
Notifications
You must be signed in to change notification settings - Fork 0
/
explainer_tutorial_binary_classification.Rmd
377 lines (264 loc) · 13 KB
/
explainer_tutorial_binary_classification.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
---
title: "classification model interpretation using explainer package"
author: "Ramtin Zargari Marandi"
output:
word_document: default
html_document:
bibliography: references.bib
citation_package: natbib
runtime: shiny
---
This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Ctrl+Shift+Enter*.
This is an example on how to use the package "explainer" developed by Ramtin Zargari Marandi (email:[ramtin.zargari.marandi\@regionh.dk](mailto:ramtin.zargari.marandi@regionh.dk){.email})
# Loading a dataset and training a machine learning model
This first code chunk loads a dataset and creates a binary classification task and train a "random forest" model using mlr3 package.
```{r}
Sys.setenv(LANG = "en") # change R language to English!
RNGkind("L'Ecuyer-CMRG") # change to L'Ecuyer-CMRG in case it uses default "Mersenne-Twister"
library("explainer")
# set seed for reproducibility
seed <- 246
set.seed(seed)
# load the BreastCancer data from the mlbench package
data("BreastCancer", package = "mlbench")
# keep the target column as "Class"
target_col <- "Class"
# change the positive class to "malignant"
positive_class <- "malignant"
# keep only the predictor variables and outcome
mydata <- BreastCancer[, -1] # 1 is ID
# remove rows with missing values
mydata <- na.omit(mydata)
# create a vector of sex categories
sex <- sample(c("Male", "Female"), size = nrow(mydata), replace = TRUE)
# create a vector of sex categories
mydata$age <- as.numeric(sample(seq(18,60), size = nrow(mydata), replace = TRUE))
# add a sex column to the mydata data frame (for fairness analysis)
mydata$sex <- factor(sex, levels = c("Male", "Female"), labels = c(1, 0))
# create a classification task
maintask <- mlr3::TaskClassif$new(id = "my_classification_task",
backend = mydata,
target = target_col,
positive = positive_class)
# create a train-test split
set.seed(seed)
splits <- mlr3::partition(maintask)
# add a learner (machine learning model base)
# library("mlr3learners")
# library("mlr3extralearners")
# mlr_learners$get("classif.randomForest")
# here we use random forest for example (you can use any other available model)
# mylrn <- mlr3::lrn("classif.randomForest", predict_type = "prob")
library("mlr3learners")
mylrn <- mlr3::lrn("classif.ranger", predict_type = "prob")
# train the model
mylrn$train(maintask, splits$train)
# make predictions on new data
mylrn$predict(maintask, splits$test)
```
## SHAP analysis to extract feature (variable) impacts on predictions
The following code chunk uses eSHAP_plot function to estimate SHAP values for the test set and create an interactive SHAP plot. This is an enhanced SHAP plot that means it provides additional information such as whether the predictions were correct (TP or TN). The color mapping provides enhanced visual inspection of the SHAP plot.
```{r}
library("magrittr")
library("plotly")
# enhanced SHAP plot
SHAP_output <- eSHAP_plot(task = maintask,
trained_model = mylrn,
splits = splits,
sample.size = 30,
seed = seed,
subset = .8)
# display the SHAP plot
myplot <- SHAP_output[[1]]
myplot
```
The following plot displays SHAP values associate with the predicted probabilities.
```{r}
SHAP_output[[5]]
```
## Visualize model performance by confusion matrix
The following code chunk uses eCM_plot function to visualize the confusion matrix to evaluate model performance for the train and test sets. More information can be found here: <https://en.wikipedia.org/wiki/Confusion_matrix> <https://cran.r-project.org/web/packages/cvms/vignettes/Creating_a_confusion_matrix.html>
```{r}
# enhanced confusion matrix
confusionmatrix_plot <- eCM_plot(task = maintask,
trained_model = mylrn,
splits = splits)
print(confusionmatrix_plot)
```
## Decision curve analysis
The provided code chunk employs the eDecisionCurve function to conduct "decision curve analysis" on the test set within the model. For an in-depth understanding of this methodology, interested readers are encouraged to explore the following authoritative references:
Decision Curve Analysis: https://en.wikipedia.org/wiki/Decision_curve_analysis
"Decision curve analysis: a novel method for evaluating prediction models" by Andrew J. Vickers and Elia B. Elkin. Link: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2577036/
These references offer comprehensive insights into the principles and applications of decision curve analysis, providing a solid foundation for further exploration and understanding of the methodology employed in the presented code.
```{r}
# enhanced decision curve plot
eDecisionCurve(task = maintask,
trained_model = mylrn,
splits = splits,
seed = seed)
```
## Model evaluation (multi-metrics and visual inspection of ROC curves)
By running the next code chunk, you will get the following model evaluation metrics and visualizations:
**AUC (Area Under the Curve):**
AUC quantifies the binary classification model's performance by assessing the area under the ROC curve, which plots sensitivity against 1-specificity across various threshold values. A value of 0.5 suggests random chance performance, while 1 signifies perfect classification.
**BACC (Balanced Accuracy):**
BACC addresses class imbalance by averaging sensitivity and specificity. Ranging from 0 to 1, a score of 0 indicates chance performance, and 1 signifies perfect classification.
**MCC (Matthews Correlation Coefficient):**
MCC evaluates binary classification model quality, considering true positives, true negatives, false positives, and false negatives. Ranging from -1 to 1, -1 represents complete disagreement, 0 implies chance performance, and 1 indicates perfect classification.
**BBRIER (Brier Score):**
BBRIER gauges the accuracy of probabilistic predictions by measuring the mean squared difference between predicted probabilities and true binary outcomes. Values range from 0 to 1, with 0 indicating perfect calibration and 1 indicating poor calibration.
**PPV (Positive Predictive Value):**
PPV, or precision, measures the proportion of true positive predictions out of all positive predictions made by the model.
**NPV (Negative Predictive Value):**
NPV quantifies the proportion of true negative predictions out of all negative predictions made by the model.
**Specificity:**
Specificity calculates the proportion of true negative predictions out of all actual negative cases in a binary classification problem.
**Sensitivity:**
Sensitivity, also known as recall or true positive rate, measures the proportion of true positive predictions out of all actual positive cases in a binary classification problem.
**PRAUC (Precision-Recall Area Under the Curve):**
PRAUC assesses binary classification model performance based on precision and recall, quantifying the area under the precision-recall curve. A PRAUC value of 1 indicates perfect classification performance.
Additionally, the analysis involves the visualization of ROC and Precision-Recall curves for both development and test sets.
```{r}
eROC_plot(task = maintask,
trained_model = mylrn,
splits = splits)
```
## ROC curves with annotated thresholds
By running the next code chunk, you will get ROC and Precision Recall curves for the development and test sets this time with probability threshold information.
```{r}
ePerformance(task = maintask,
trained_model = mylrn,
splits = splits)
```
## loading SHAP results for downstream analysis
Now we can get the outputs from eSHAP_plot function to apply clustering on SHAP values
```{r}
shap_Mean_wide <- SHAP_output[[2]]
shap_Mean_long <- SHAP_output[[3]]
shap <- SHAP_output[[4]]
```
## SHAP values in association with feature values
```{r}
ShapFeaturePlot(shap_Mean_long)
```
# partial dependence of features
Partial dependence plots (PDPs): PDPs can be used to visualize the marginal effect of a single feature on the model prediction.
```{r}
ShapPartialPlot(shap_Mean_long = shap_Mean_long)
```
# extract feature values and predicted probabilities as output of the model to analyze
```{r}
fval_predprob <- reshape2::dcast(shap, sample_num + pred_prob + predcorrectness ~ feature, value.var = "feature.value")
```
# Run a Shiny app to visualize 2-way partial dependence plots
```{r}
library(shiny)
library(ggplot2)
# assuming your data.table is named `fval_predprob`
ui <- fluidPage(
titlePanel("Feature Plot"),
sidebarLayout(
sidebarPanel(
selectInput("x_feature", "X-axis Feature:",
choices = names(fval_predprob)),
selectInput("y_feature", "Y-axis Feature:",
choices = names(fval_predprob)),
width = 2,
actionButton("stop_button", "Stop")
),
mainPanel(
plotOutput("feature_plot")
)
)
)
server <- function(input, output) {
# create a reactive value for the app state
app_state <- reactiveValues(running = TRUE)
output$feature_plot <- renderPlot({
ggplot(fval_predprob, aes(x = .data[[input$x_feature]], y = .data[[input$y_feature]], fill = pred_prob)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "steelblue", limits = c(0, 1)) +
xlab(input$x_feature) +
ylab(input$y_feature) +
labs(shape = "correct prediction") +
egg::theme_article()
})
# observe the stop button
observeEvent(input$stop_button, {
app_state$running <- FALSE
})
# stop the app if the running state is FALSE
observe({
if(!app_state$running) {
stopApp()
}
})
}
shinyApp(ui = ui, server = server)
```
```{r}
str(shap_Mean_long)
```
```{r}
# This is the data table that includes SHAP values in wide format. Each row contains sample_num as sample ID and the SHAP values for each feature.
shap_Mean_wide
```
```{r}
# This is the data table of SHAP values in long format and includes feature name, mean_phi_test: mean shap value for the feature across samples, scaled feature values, shap values for each feature, sample ID, and whether the prediction was correct (i.e. predicted class = actual class)
shap_Mean_long
```
# Patient subgroups determined by SHAP clusters
SHAP clustering is a method to better understand why a model may perform better for some patients than others. Here for example, you can identify patient subgroups that have specific patterns that are different from other subgroups and that explains why the model you developed have perhaps better or worse performance for those patients than the average performance for the whole dataset. You can see the difference in the SHAP plots that if you group all together provide the overall SHAP summary plot. Again here the edges reflect how features may interact with each other in each individual sample (instance).
```{r}
# the number of clusters can be changed
SHAP_plot_clusters <- SHAPclust(task = maintask,
trained_model = mylrn,
splits = splits,
shap_Mean_wide = shap_Mean_wide,
shap_Mean_long = shap_Mean_long,
num_of_clusters = 3,
seed = seed,
subset = .8)
# note that the subset must be the same value as the SHAP analysis done earlier
# display the SHAP cluster plots
SHAP_plot_clusters[[1]]
```
```{r}
# display the confusion matrices corresponding to the SHAP clusters (patient subsets determined by SHAP clusters)
SHAP_plot_clusters[[2]]
```
# Model fairness (sensitivity analysis)
Sometimes we would like to investigate whether our model performs fairly well for different subgroups based on categories of variables such as sex.
```{r}
# you should decide what variables to use to be tested
# here we chose sex from the variables existing in the dataset
Fairness_results <- eFairness(task = maintask,
trained_model = mylrn,
splits = splits,
target_variable = "sex",
var_levels = c("Male", "Female"))
# ROC curves for the subgroups for the development (left) and test (right) sets
Fairness_results[[1]]
# performance in the subgroups for the development set
Fairness_results[[2]]
# performance in the subgroups for the test set
Fairness_results[[3]]
```
## Model parameters
```{r}
# get model parameters
model_params <- mylrn$param_set
print(data.table::as.data.table(model_params))
```
## Report packages that have been used
The current R package has been developed with the following dependencies, ensuring its functionality and compatibility:
For a detailed reference on each package, please consult the respective documentation and citation information available on CRAN (Comprehensive R Archive Network) or the official package repositories.
```{r}
library("report") # to report packages have been used
oursessionreport <- report(sessionInfo())
summary(oursessionreport)
oursessionreport_df <- as.data.frame(oursessionreport)
data.table::fwrite(oursessionreport_df, file = paste0("sessioninfo",seed,".xlsx"))
summary(as.data.frame(oursessionreport))
```