ModelStu001

ModelStu001

library("Numero")
library("DALEX")

## Welcome to DALEX (version: 1.0).
## Find examples and detailed introduction at: https://pbiecek.github.io/ema/
## Additional features will be available after installation of: ggpubr.
## Use 'install_dependencies()' to get all suggested dependencies

library("ranger")
library(data.table)
setwd("C:/Users/s-das/Syncplicity/MyProjects_IMP/MY_Papers_V2/TRB 2021/EScotter_BayesianRule/")
it01 <- fread("IT_aadtMaster.csv")

## Warning in require_bit64_if_needed(ans): Some columns are type 'integer64'
## but package bit64 is not installed. Those columns will print as strange
## looking floating point data. There is no need to reload the data. Simply
## install.packages('bit64') to obtain the integer64 print method and print the
## data again.

mn= subset(it01, State=="MN")
dim(mn)

## [1] 11498    86

mn01 <- mn[, c("FC_RU", "Default_AADT", "HU", "Pop" , "WAC", "RAC", "Agg_Inc", "Agg_Veh", "Empl" )]
mn02= na.omit(mn01)
mn02$FC_RU= as.factor(mn02$FC_RU)

model_HR_ranger <- ranger(FC_RU~.,  data = mn02, probability = TRUE, num.trees = 50)
model_HR_ranger

## Ranger result
## 
## Call:
##  ranger(FC_RU ~ ., data = mn02, probability = TRUE, num.trees = 50) 
## 
## Type:                             Probability estimation 
## Number of trees:                  50 
## Sample size:                      10780 
## Number of independent variables:  8 
## Mtry:                             2 
## Target node size:                 10 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error (Brier s.):  0.1744989

explain_HR_ranger <- explain(model_HR_ranger,
                             data = mn02[,-1],
                             y = mn02$FC_RU,
                             label = "Ranger Multilabel Classification",
                             colorize = FALSE)

## Preparation of a new explainer is initiated
##   -> model label       :  Ranger Multilabel Classification 
##   -> data              :  10780  rows  8  cols 
##   -> target variable   :  10780  values 
##   -> target variable   :  Please note that 'y' is a factor.  (  WARNING  )
##   -> target variable   :  Consider changing the 'y' to a logical or numerical vector.
##   -> target variable   :  Otherwise I will not be able to calculate residuals or loss function.
##   -> model_info        :  package ranger , ver. 0.12.1 , task classification (  default  ) 
##   -> predict function  :  yhat.ranger  will be used (  default  )
##   -> predicted values  :  predict function returns multiple columns:  3  (  WARNING  ) some of functionalities may not work 
##   -> residual function :  difference between y and yhat (  default  )

## Warning in Ops.factor(y, predict_function(model, data)): '-' not meaningful for
## factors

##   -> residuals         :  numerical, min =  NA , mean =  NA , max =  NA  
##   A new explainer has been created!

y_HR <- as.numeric(mn02$FC_RU)
explain_HR_ranger_new_y <- explain(model_HR_ranger,
                                   data = mn02[,-1],
                                   y = y_HR,
                                   label = "Ranger Multilabel Classification",
                                   colorize = FALSE)

## Preparation of a new explainer is initiated
##   -> model label       :  Ranger Multilabel Classification 
##   -> data              :  10780  rows  8  cols 
##   -> target variable   :  10780  values 
##   -> model_info        :  package ranger , ver. 0.12.1 , task classification (  default  ) 
##   -> predict function  :  yhat.ranger  will be used (  default  )
##   -> predicted values  :  predict function returns multiple columns:  3  (  WARNING  ) some of functionalities may not work 
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  0.002222222 , mean =  1.727706 , max =  3  
##   A new explainer has been created!

vi <- variable_importance(explain_HR_ranger_new_y, loss_function = loss_cross_entropy)
plot(vi)

ve_p <- variable_profile(explain_HR_ranger, variables = "Default_AADT", type = "partial")

## Warning in if (class(new_observation_ext) != "data.frame") {: the condition has
## length > 1 and only the first element will be used

ve_p$color = "_label_"
plot(ve_p)

bd <- variable_attribution(explain_HR_ranger, mn02[2,], type = "break_down")
plot(bd)

shap <- variable_attribution(explain_HR_ranger, mn02[2,], type = "shap")
plot(shap)