Data for Titanic survival

Let’s see an example for DALEX package for classification models for the survival problem for Titanic dataset. Here we are using a dataset titanic avaliable in teh DALEX package. Note that this data was copied from the stablelearner package.

library("DALEX")
head(titanic)
#>   gender age class    embarked       country  fare sibsp parch survived
#> 1   male  42   3rd Southampton United States  7.11     0     0       no
#> 2   male  13   3rd Southampton United States 20.05     0     2       no
#> 3   male  16   3rd Southampton United States 20.05     1     1       no
#> 4 female  39   3rd Southampton       England 20.05     1     1      yes
#> 5 female  16   3rd Southampton        Norway  7.13     0     0      yes
#> 6   male  25   3rd Southampton United States  7.13     0     0      yes

Model for Titanic survival

Ok, not it’s time to create a model. Let’s use the Random Forest model.

# prepare model
library("randomForest")
titanic <- na.omit(titanic)
model_titanic_rf <- randomForest(survived == "yes" ~ gender + age + class + embarked +
                                   fare + sibsp + parch,  data = titanic)
model_titanic_rf
#> 
#> Call:
#>  randomForest(formula = survived == "yes" ~ gender + age + class +      embarked + fare + sibsp + parch, data = titanic) 
#>                Type of random forest: regression
#>                      Number of trees: 500
#> No. of variables tried at each split: 2
#> 
#>           Mean of squared residuals: 0.1427006
#>                     % Var explained: 34.89

Explainer for Titanic survival

The third step (it’s optional but useful) is to create a DALEX explainer for random forest model.

library("DALEX")
explain_titanic_rf <- explain(model_titanic_rf, 
                      data = titanic[,-9],
                      y = titanic$survived == "yes", 
                      label = "Random Forest v7")

Variable importance plots

Use the variable_importance() explainer to present importance of particular features. Note that type = "difference" normalizes dropouts, and now they all start in 0.

vi_rf <- variable_importance(explain_titanic_rf)
head(vi_rf)
#>       variable dropout_loss            label
#> 1 _full_model_     112.0961 Random Forest v7
#> 2      country     112.0961 Random Forest v7
#> 3        sibsp     119.4988 Random Forest v7
#> 4        parch     119.8038 Random Forest v7
#> 5     embarked     121.7538 Random Forest v7
#> 6         fare     138.1498 Random Forest v7
plot(vi_rf)

Variable effects

As we see the most important feature is Sex. Next three importnat features are Pclass, Age and Fare. Let’s see the link between model response and these features.

Such univariate relation can be calculated with variable_response().

Age

Kids 5 years old and younger have much higher survival probability.

vr_age  <- variable_response(explain_titanic_rf, variable =  "age")
head(vr_age)
#>           x            y var type            label
#> 1 0.1666667 0.526169.... age  pdp Random Forest v7
#> 2 1.6433333 0.559907.... age  pdp Random Forest v7
#> 3 3.1200000 0.575562.... age  pdp Random Forest v7
#> 4 4.5966667 0.543471.... age  pdp Random Forest v7
#> 5 6.0733333 0.519516.... age  pdp Random Forest v7
#> 6 7.5500000 0.520232.... age  pdp Random Forest v7
plot(vr_age, use_facets = TRUE)

Passanger class

Passangers in the first class have much higher survival probability.

vr_class  <- variable_response(explain_titanic_rf, variable =  "class")
plot(vr_class)

Fare

Very cheap tickets are linked with lower chances.

vr_fare  <- variable_response(explain_titanic_rf, variable =  "fare")
plot(vr_fare, use_facets = TRUE)

Embarked

Passangers that embarked from C have highest survival.

vr_embarked  <- variable_response(explain_titanic_rf, variable =  "embarked")
plot(vr_embarked)

Instance level explanations

Let’s see break down explanation for model predictions for 8 years old male from 1st class that embarked from port C.

new_passanger <- data.frame(
  class = factor("1st", levels = c("1st", "2nd", "3rd", "deck crew", "engineering crew", "restaurant staff", "victualling crew")),
  gender = factor("male", levels = c("female", "male")),
  age = 8,
  sibsp = 0,
  parch = 0,
  fare = 72,
  embarked = factor("Southampton", levels = c("Belfast", "Cherbourg", "Queenstown", "Southampton"))
)

sp_rf <- single_prediction(explain_titanic_rf, new_passanger)
plot(sp_rf)

It looks like the most important feature for this passenger is age and sex. After all his odds for survival are higher than for the average passenger. Mainly because of the young age and despite of being a male.

More models

Let’s train more models for survival.

Logistic regression

library("rms")
model_titanic_lmr <- lrm(survived == "yes" ~ class + gender + rcs(age) + sibsp +
                   parch + fare + embarked, titanic)
explain_titanic_lmr <- explain(model_titanic_lmr, data = titanic, 
                       y = titanic$survived == "yes", 
                       predict_function = function(m,x) predict(m, x, type="fitted"),
                       label = "Logistic regression")

Generalized Boosted Models (GBM)

library("gbm")
model_titanic_gbm <- gbm(survived == "yes" ~ class + gender + age + sibsp +
                     parch + fare + embarked, data = titanic, n.trees = 15000)
#> Distribution not specified, assuming bernoulli ...
explain_titanic_gbm <- explain(model_titanic_gbm, data = titanic, 
                       y = titanic$survived == "yes", 
                       predict_function = function(m,x) predict(m, x, n.trees = 15000, type = "response"),
                       label = "Generalized Boosted Models")

Support Vector Machines (SVM)

library("e1071")
model_titanic_svm <- svm(survived == "yes" ~ class + gender + age + sibsp +
                     parch + fare + embarked, data = titanic, 
             type = "C-classification", probability = TRUE)
explain_titanic_svm <- explain(model_titanic_svm, data = titanic, 
                       y = titanic$survived == "yes", 
                       label = "Support Vector Machines")

k-Nearest Neighbours (kNN)

library("caret")
model_titanic_knn <- knn3(survived == "yes" ~ class + gender + age + sibsp +
                     parch + fare + embarked, data = titanic, k = 5)
explain_titanic_knn <- explain(model_titanic_knn, data = titanic, 
                       y = titanic$survived == "yes", 
                       predict_function = function(m,x) predict(m, x)[,2],
                       label = "k-Nearest Neighbours")

Variable performance

vi_rf <- variable_importance(explain_titanic_rf)
vi_lmr <- variable_importance(explain_titanic_lmr)
vi_gbm <- variable_importance(explain_titanic_gbm)
vi_svm <- variable_importance(explain_titanic_svm)
vi_knn <- variable_importance(explain_titanic_knn)

plot(vi_rf, vi_lmr, vi_gbm, vi_svm, vi_knn, bar_width = 4)

Single variable

vr_age_rf  <- variable_response(explain_titanic_rf, variable =  "age")
vr_age_lmr  <- variable_response(explain_titanic_lmr, variable =  "age")
vr_age_gbm  <- variable_response(explain_titanic_gbm, variable =  "age")
vr_age_svm  <- variable_response(explain_titanic_svm, variable =  "age")
vr_age_knn  <- variable_response(explain_titanic_knn, variable =  "age")
plot(vr_age_rf, vr_age_lmr, vr_age_gbm, vr_age_svm, vr_age_knn)

plot(vr_age_rf, vr_age_lmr, vr_age_gbm, vr_age_svm, vr_age_knn, use_facets = TRUE)

Instance level explanations

sp_rf <- single_prediction(explain_titanic_rf, new_passanger)
sp_lmr <- single_prediction(explain_titanic_lmr, new_passanger)
sp_gbm <- single_prediction(explain_titanic_gbm, new_passanger)
sp_svm <- single_prediction(explain_titanic_svm, new_passanger)
sp_knn <- single_prediction(explain_titanic_knn, new_passanger)
plot(sp_rf, sp_lmr, sp_gbm, sp_svm, sp_knn)

Session info

#> R version 3.5.3 (2019-03-11)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.2 LTS
#> 
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] caret_6.0-81        e1071_1.7-0.1       gbm_2.1.5          
#>  [4] rms_5.1-3           SparseM_1.77        Hmisc_4.2-0        
#>  [7] ggplot2_3.1.1       Formula_1.2-3       survival_2.43-3    
#> [10] lattice_0.20-38     randomForest_4.6-14 DALEX_0.4          
#> 
#> loaded via a namespace (and not attached):
#>   [1] TH.data_1.0-10      colorspace_1.4-1    deldir_0.1-16      
#>   [4] class_7.3-15        rprojroot_1.3-2     htmlTable_1.13.1   
#>   [7] base64enc_0.1-3     fs_1.2.7            rstudioapi_0.10    
#>  [10] proxy_0.4-23        roxygen2_6.1.1      ggpubr_0.2         
#>  [13] MatrixModels_0.4-1  lubridate_1.7.4     prodlim_2018.04.18 
#>  [16] mvtnorm_1.0-10      xml2_1.2.0          codetools_0.2-16   
#>  [19] splines_3.5.3       knitr_1.22          cluster_2.0.7-1    
#>  [22] shiny_1.3.2         compiler_3.5.3      backports_1.1.4    
#>  [25] assertthat_0.2.1    Matrix_1.2-17       lazyeval_0.2.2     
#>  [28] later_0.8.0         acepack_1.4.1       htmltools_0.3.6    
#>  [31] quantreg_5.38       tools_3.5.3         coda_0.19-2        
#>  [34] gtable_0.3.0        agricolae_1.3-0     glue_1.3.1         
#>  [37] reshape2_1.4.3      dplyr_0.8.0.1       gmodels_2.18.1     
#>  [40] Rcpp_1.0.1          pkgdown_1.3.0       spdep_1.0-2        
#>  [43] gdata_2.18.0        nlme_3.1-137        iterators_1.0.10   
#>  [46] timeDate_3043.102   gower_0.1.2         xfun_0.6           
#>  [49] stringr_1.4.0       mime_0.6            miniUI_0.1.1.1     
#>  [52] breakDown_0.1.6     gtools_3.8.1        polspline_1.1.14   
#>  [55] zoo_1.8-5           LearnBayes_2.15.1   MASS_7.3-51.1      
#>  [58] scales_1.0.0        ipred_0.9-8         promises_1.0.1     
#>  [61] sandwich_2.5-0      expm_0.999-3        RColorBrewer_1.1-2 
#>  [64] yaml_2.2.0          memoise_1.1.0       gridExtra_2.3      
#>  [67] rpart_4.1-13        latticeExtra_0.6-28 stringi_1.4.3      
#>  [70] highr_0.8           klaR_0.6-14         AlgDesign_1.1-7.3  
#>  [73] desc_1.2.0          foreach_1.4.4       checkmate_1.9.1    
#>  [76] boot_1.3-20         lava_1.6.5          spData_0.3.0       
#>  [79] rlang_0.3.4         pkgconfig_2.0.2     commonmark_1.7     
#>  [82] evaluate_0.13       purrr_0.3.2         sf_0.7-3           
#>  [85] recipes_0.1.4       htmlwidgets_1.3     labeling_0.3       
#>  [88] cowplot_0.9.4       tidyselect_0.2.5    factorMerger_0.3.6 
#>  [91] plyr_1.8.4          magrittr_1.5        R6_2.4.0           
#>  [94] generics_0.0.2      multcomp_1.4-10     combinat_0.0-8     
#>  [97] DBI_1.0.0           pillar_1.3.1        foreign_0.8-70     
#> [100] withr_2.1.2         units_0.6-2         sp_1.3-1           
#> [103] nnet_7.3-12         tibble_2.1.1        crayon_1.3.4       
#> [106] questionr_0.7.0     rmarkdown_1.12      grid_3.5.3         
#> [109] data.table_1.12.2   ModelMetrics_1.2.2  digest_0.6.18      
#> [112] classInt_0.3-1      pdp_0.7.0           xtable_1.8-4       
#> [115] httpuv_1.5.1        stats4_3.5.3        munsell_0.5.0