Ensemble methods
Introduction
Today, we will learn how to use different ensemble methods in R, recap on how to evaluate the performance of the methods, and learn how we can substantively interpret the model output.
In this practical we will work with the ILPD (Indian Liver Patient Dataset) from the UCI Machine Learning Repository (you can find the data here). This data set contains data on 414 liver disease patients, and 165 non-patients. In general, medical researchers have two distinct goals when doing research: (1) to be able to classify people in their waiting room as either patients or non-patients, and (2) get insight into the factors that are associated with the disease. In this practical, we will look at both aspects.
In this practical, we will use the tidyverse, magrittr, psych, GGally and caret libraries.
library(tidyverse)
library(magrittr)
library(psych)
library(caret)
library(gbm)
library(xgboost)
library(data.table)
library(ggforce)First, we specify a seed and load the training data. We will use this data to make inferences and to train a prediction model.
set.seed(45)
df <- readRDS("data/train_disease.RDS")1. Get an impression of the data by looking at the structure of the data and creating some descriptive statistics.
head(df)tail(df)df %>%
select(-c(Gender, Disease)) %>%
describeBy(df$Disease, fast = TRUE)##
## Descriptive statistics by group
## group: Healthy
## vars n mean sd min max range se
## Age 1 360 45.81 15.60 10.0 90.0 80.0 0.82
## Total_Bilirubin 2 360 4.05 7.15 0.4 75.0 74.6 0.38
## Direct_Bilirubin 3 360 1.86 3.16 0.1 19.7 19.6 0.17
## Alkaline_Phosphotase 4 360 309.47 250.37 63.0 2110.0 2047.0 13.20
## Alamine_Aminotransferase 5 360 97.19 200.17 12.0 1680.0 1668.0 10.55
## Aspartate_Aminotransferase 6 360 133.66 324.59 11.0 4929.0 4918.0 17.11
## Total_Protiens 7 360 6.47 1.11 2.7 9.6 6.9 0.06
## Albumin 8 360 3.07 0.80 0.9 5.5 4.6 0.04
## Ratio_Albumin_Globulin 9 360 0.92 0.32 0.3 2.8 2.5 0.02
## ------------------------------------------------------------
## group: Disease
## vars n mean sd min max range se
## Age 1 140 41.51 16.98 4.00 85.0 81.00 1.43
## Total_Bilirubin 2 140 1.14 1.01 0.50 7.3 6.80 0.09
## Direct_Bilirubin 3 140 0.40 0.53 0.10 3.6 3.50 0.04
## Alkaline_Phosphotase 4 140 220.56 149.76 90.00 1580.0 1490.00 12.66
## Alamine_Aminotransferase 5 140 33.31 24.50 10.00 181.0 171.00 2.07
## Aspartate_Aminotransferase 6 140 39.79 34.79 10.00 285.0 275.00 2.94
## Total_Protiens 7 140 6.53 1.07 3.70 9.2 5.50 0.09
## Albumin 8 140 3.34 0.80 1.40 5.0 3.60 0.07
## Ratio_Albumin_Globulin 9 140 1.04 0.30 0.37 1.9 1.53 0.03# It becomes directly clear that there are substantial differences between
# the diseased and non-diseased in the data.2. To further explore the data we work with, create some interesting data visualizations that show whether there are interesting patterns in the data.
Hint: Think about adding a color aesthetic for the variable Disease.
df %>%
select(-Gender) %>%
pivot_longer(where(is.numeric)) %>%
ggplot(aes(x = value, col = Disease, fill = Disease)) +
geom_boxplot(alpha = 0.8) +
facet_wrap(~name, scales = "free") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme_minimal()
df %>%
select(-Gender) %>%
pivot_longer(where(is.numeric)) %>%
ggplot(aes(x = value, col = Disease, fill = Disease)) +
geom_density(alpha = 0.8) +
facet_wrap(~name, scales = "free") +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
theme_minimal()
prop.table(table(df$Gender, df$Disease), margin = 1) %>%
as.data.frame %>%
select(Gender = Var1, Disease = Var2, `Relative Frequency` = Freq) %>%
ggplot(aes(y = `Relative Frequency`, x = Gender, col = Disease, fill = Disease)) +
geom_histogram(alpha = 0.8, stat = "identity", position = "dodge") +
scale_fill_brewer(palette = "Paired") +
scale_color_brewer(palette = "Paired") +
theme_minimal()## Warning: Ignoring unknown parameters: binwidth, bins, pad
## There are some differences between distributions for the two
## Disease categories, but the differences do not seem to be dramatic.
## Additionally, there are relatively more women with the liver
## disease than men.3. Shortly reflect on the difference between bagging, random forests, and boosting.
## Bagging: fit a regression tree to N bootstrap samples of the training data
## take the average of all classification trees to base predictions on
## Note: out-of-bag data can serve as internal validation set.
## Random forest: Similarly to bagging, classification trees are trained on
## a bootstrap sample of the data. However, the decision trees
## are trained using less than all features in the data.
## Boosting: We build a decision tree sequentially. Given the current
## we fit a (small) tree on the residuals of the current model,
## rather than on the outcome YWe are going to apply different machine learning models using the caret package.
4. Apply bagging to the training data, to predict the outcome Disease, using the caret library.
Note. We first specify the internal validation settings, like so:
cvcontrol <- trainControl(method = "repeatedcv",
number = 10,
allowParallel = TRUE)These settings can be inserted within the train function from the caret package. Make sure to use the treebag method, to specify cvcontrol as the trControl argument and to set importance = TRUE.
bag_train <- train(Disease ~ .,
data = df,
method = 'treebag',
trControl = cvcontrol,
importance = TRUE)5. Interpret the variable importance measure using the varImp function on the trained model object.
bag_train %>%
varImp %>%
plot
6. Create training set predictions based on the bagged model, and use the confusionMatrix() function from the caret package to assess it’s performance.`
Hint: You will have to create predictions based on the trained model for the training data, and evaluate these against the observed values of the training data.
confusionMatrix(predict(bag_train, type = "raw"),
df$Disease)## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Disease
## Healthy 360 0
## Disease 0 140
##
## Accuracy : 1
## 95% CI : (0.9926, 1)
## No Information Rate : 0.72
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.00
## Specificity : 1.00
## Pos Pred Value : 1.00
## Neg Pred Value : 1.00
## Prevalence : 0.72
## Detection Rate : 0.72
## Detection Prevalence : 0.72
## Balanced Accuracy : 1.00
##
## 'Positive' Class : Healthy
## ## We have achieved a perfect training set performance. However, this shows
## nothing more than that we have been able to train the model. We need to
## evaluate our 7. Now ask for the output of the bagged model. Explain why the under both approaches differ.
bag_train## Bagged CART
##
## 500 samples
## 10 predictor
## 2 classes: 'Healthy', 'Disease'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 1 times)
## Summary of sample sizes: 450, 450, 450, 450, 450, 450, ...
## Resampling results:
##
## Accuracy Kappa
## 0.678 0.1311671We will now follow the same approach, but rather than bagging, we will train a random forest on the training data.
8. Fit a random forest to the training data to predict the outcome Disease, using the caret library.
Note. Use the same cvcontrol settings as in the previous model.
rf_train <- train(Disease ~ .,
data = df,
method = 'rf',
trControl = cvcontrol,
importance = TRUE)9. Again, interpret the variable importance measure using the varImp function on the trained model object. Do you draw the same conclusions as under the bagged model?
rf_train %>%
varImp %>%
plot
## The random forest model indicates that other variables are more important,
## as compared to the bagged model.10. Output the model output from the random forest. Are we doing better than with the bagged model?
rf_train## Random Forest
##
## 500 samples
## 10 predictor
## 2 classes: 'Healthy', 'Disease'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 1 times)
## Summary of sample sizes: 450, 450, 450, 450, 450, 450, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.688 0.1179249
## 6 0.692 0.1649948
## 10 0.692 0.1661658
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 6.## Yes, the most accurate model indicates that we do just slightly better than
## with the bagged model. However, this might well be due to chance.11. Now, fit a boosting model using the caret library to predict disease status.`
Hint: Use gradient boosting (the gbm method in caret).
gbm_train <- train(Disease ~ .,
data = df,
method = "gbm",
verbose = F,
trControl = cvcontrol)12. Again, interpret the variable importance measure. You will have to call for summary() on the model object you just created. Compare the output to the previously obtained variable importance measures.
summary(gbm_train)
13. Output the model output from our gradient boosting procedure. Are we doing better than with the bagged and random forest model?
gbm_train## Stochastic Gradient Boosting
##
## 500 samples
## 10 predictor
## 2 classes: 'Healthy', 'Disease'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 1 times)
## Summary of sample sizes: 450, 450, 450, 450, 450, 450, ...
## Resampling results across tuning parameters:
##
## interaction.depth n.trees Accuracy Kappa
## 1 50 0.712 0.09569625
## 1 100 0.712 0.16033979
## 1 150 0.702 0.15411887
## 2 50 0.694 0.09651620
## 2 100 0.684 0.12744126
## 2 150 0.676 0.12138676
## 3 50 0.688 0.12456452
## 3 100 0.698 0.15242249
## 3 150 0.694 0.13992459
##
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
##
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were n.trees = 50, interaction.depth =
## 1, shrinkage = 0.1 and n.minobsinnode = 10.# Yes, our best model is doing slightly better then the previous two models.
# However, this might still be random variation.For now, we will continue with extreme gradient boosting, although we will use a difference procedure.
We will use xgboost to train a binary classification model, and create some visualizations to obtain additional insight in our model. We will create the visualizations using SHAP (SHapley Additive exPlanations) values, which are a measure of importance of the variables in the model. In fact, SHAP values indicate the influence of each input variable on the predicted probability for each person. Essentially, these give an indication of the difference between the predicted probability with and without that variable, for each person’s score.
14. Download the file shap.R from this Github repository.
Note. There are multiple ways to this, of which the simplest is to run the following code.
library(devtools)
source_url("https://github.com/pablo14/shap-values/blob/master/shap.R?raw=TRUE")Additionally, you could simply go to the file shap.R and copy-and-paste the code into the current repository. However, you could also fork and clone the repository, to make adjustments to the functions that are already created.
15. Specify your model as follows, and use it to create predictions on the training data.
train_x <- model.matrix(Disease ~ ., df)[,-1]
train_y <- as.numeric(df$Disease) - 1
xgboost_train <- xgboost(data = train_x,
label = train_y,
max.depth = 10,
eta = 1,
nthread = 4,
nrounds = 4,
objective = "binary:logistic",
verbose = 2)
pred <- tibble(Disease = predict(xgboost_train, newdata = train_x)) %>%
mutate(Disease = factor(ifelse(Disease < 0.5, 1, 2),
labels = c("Healthy", "Disease")))
table(pred$Disease, df$Disease)16. First, calculate the SHAP rank scores for all variables in the data, and create a variable importance plot using these values. Interpret the plot.
shap_results <- shap.score.rank(xgboost_train,
X_train = train_x,
shap_approx = F)## make SHAP score by decreasing ordervar_importance(shap_results)17. Plot the SHAP values for every individual for every feature and interpret them.
shap_long <- shap.prep(shap = shap_results,
X_train = train_x)
plot.shap.summary(shap_long)
xgb.plot.shap(train_x, features = colnames(train_x), model = xgboost_train, n_col = 3)
## The first plot shows, for example, that those with a high value for
## Direct_Bilirubin have a lower probability of being diseased. Also,
## Those with a higher age have a lower probability of being diseased,
## while those with a higher Albumin have a higher probability of being diseased.
## The second set of plots displays the marginal relationships of the SHAP values with the predictors. This conveys the same information, but in greater detail. The interpretability may be a bit tricky for the inexperienced data analyst. 18. Verify which of the models you created in this practical performs best on the test data.
test <- readRDS("data/test_disease.RDS")
bag_test <- predict(bag_train, newdata = test)
rf_test <- predict(rf_train, newdata = test)
gbm_test <- predict(gbm_train, newdata = test)
xgb_test <- predict(xgboost_train, newdata = model.matrix(Disease ~ ., test)[,-1]) %>%
factor(x = ifelse(. < 0.5, 1, 2), levels = c(1,2), labels = c("Healthy", "Disease"))
list(bag_test,
rf_test,
gbm_test,
xgb_test) %>%
map(~ confusionMatrix(.x, test$Disease))## [[1]]
## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Disease
## Healthy 48 16
## Disease 6 9
##
## Accuracy : 0.7215
## 95% CI : (0.6093, 0.8165)
## No Information Rate : 0.6835
## P-Value [Acc > NIR] : 0.27598
##
## Kappa : 0.2788
##
## Mcnemar's Test P-Value : 0.05501
##
## Sensitivity : 0.8889
## Specificity : 0.3600
## Pos Pred Value : 0.7500
## Neg Pred Value : 0.6000
## Prevalence : 0.6835
## Detection Rate : 0.6076
## Detection Prevalence : 0.8101
## Balanced Accuracy : 0.6244
##
## 'Positive' Class : Healthy
##
##
## [[2]]
## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Disease
## Healthy 51 18
## Disease 3 7
##
## Accuracy : 0.7342
## 95% CI : (0.6228, 0.8273)
## No Information Rate : 0.6835
## P-Value [Acc > NIR] : 0.19982
##
## Kappa : 0.2675
##
## Mcnemar's Test P-Value : 0.00225
##
## Sensitivity : 0.9444
## Specificity : 0.2800
## Pos Pred Value : 0.7391
## Neg Pred Value : 0.7000
## Prevalence : 0.6835
## Detection Rate : 0.6456
## Detection Prevalence : 0.8734
## Balanced Accuracy : 0.6122
##
## 'Positive' Class : Healthy
##
##
## [[3]]
## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Disease
## Healthy 53 21
## Disease 1 4
##
## Accuracy : 0.7215
## 95% CI : (0.6093, 0.8165)
## No Information Rate : 0.6835
## P-Value [Acc > NIR] : 0.276
##
## Kappa : 0.1802
##
## Mcnemar's Test P-Value : 5.104e-05
##
## Sensitivity : 0.9815
## Specificity : 0.1600
## Pos Pred Value : 0.7162
## Neg Pred Value : 0.8000
## Prevalence : 0.6835
## Detection Rate : 0.6709
## Detection Prevalence : 0.9367
## Balanced Accuracy : 0.5707
##
## 'Positive' Class : Healthy
##
##
## [[4]]
## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Disease
## Healthy 44 14
## Disease 10 11
##
## Accuracy : 0.6962
## 95% CI : (0.5825, 0.7947)
## No Information Rate : 0.6835
## P-Value [Acc > NIR] : 0.4577
##
## Kappa : 0.2663
##
## Mcnemar's Test P-Value : 0.5403
##
## Sensitivity : 0.8148
## Specificity : 0.4400
## Pos Pred Value : 0.7586
## Neg Pred Value : 0.5238
## Prevalence : 0.6835
## Detection Rate : 0.5570
## Detection Prevalence : 0.7342
## Balanced Accuracy : 0.6274
##
## 'Positive' Class : Healthy
## Hand-in
When you have finished the practical,
enclose all files of the project
09_ensemble.Rproj(i.e. all.Rand/or.Rmdfiles including the one with your answers, and the.Rprojfile) in a zip file, andhand in the zip by PR from your fork here. Try do so as soon as possible.