02 December 2016

In last week’s post I explored whether machine learning models can be applied to predict flu deaths from the 2013 outbreak of influenza A H7N9 in China. There, I compared random forests, elastic-net regularized generalized linear models, k-nearest neighbors, penalized discriminant analysis, stabilized linear discriminant analysis, nearest shrunken centroids, single C5.0 tree and partial least squares.


Extreme gradient boosting

Extreme gradient boosting (XGBoost) is a faster and improved implementation of gradient boosting for supervised learning and has recently been very successfully applied in Kaggle competitions. Because I’ve heard XGBoost’s praise being sung everywhere lately, I wanted to get my feet wet with it too. So this week I want to compare the prediction success of gradient boosting with the same dataset. Additionally, I want to test the influence of different preprocessing methods on the outcome.


“XGBoost uses a more regularized model formalization to control over-fitting, which gives it better performance.” Tianqi Chen, developer of xgboost

XGBoost is a tree ensemble model, which means the sum of predictions from a set of classification and regression trees (CART). In that, XGBoost is similar to Random Forests but it uses a different approach to model training.


Starting with the same test and training data (partitioned into validation test and validation train subsets) from last week’s post, I am training extreme gradient boosting models as implemented in the xgboost and caret packages with different preprocessing settings.

Out of the different implementations and variations of gradient boosting algorithms, caret performed best on PCA-preprocessed data in the validation set. These paramteres were then used to predict the outcome in the test set and compare it to last week’s predictions.

Compared to last week, there is much less uncertainty in the predictions from XGBoost. Overall, I would say that this algorithm is superior to the others I have used before.


xgboost

Extreme gradient boosting is implemented in the xgboost package.

# install the stable/pre-compiled version from CRAN

install.packages('xgboost')

# or install from weekly updated drat repo

install.packages("drat", repos="https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("xgboost", repos="http://dmlc.ml/drat/", type="source")

XGBoost supports only numbers, so the outcome classes have to be converted into integers and both training and test data have to be in numeric matrix format.

library(xgboost)
matrix_train <- apply(val_train_X, 2, function(x) as.numeric(as.character(x)))
outcome_death_train <- ifelse(val_train_data$outcome == "Death", 1, 0)

matrix_test <- apply(val_test_X, 2, function(x) as.numeric(as.character(x)))
outcome_death_test <- ifelse(val_test_data$outcome == "Death", 1, 0)

xgb_train_matrix <- xgb.DMatrix(data = as.matrix(matrix_train), label = outcome_death_train)
xgb_test_matrix <- xgb.DMatrix(data = as.matrix(matrix_test), label = outcome_death_test)

watchlist <- list(train = xgb_train_matrix, test = xgb_test_matrix)
label <- getinfo(xgb_test_matrix, "label")

I am using cross validation to evaluate the error rate.

param <- list("objective" = "binary:logistic")

xgb.cv(param = param, 
       data = xgb_train_matrix, 
       nfold = 3,
       label = getinfo(xgb_train_matrix, "label"),
       nrounds = 5)
## [1]  train-error:0.133950+0.022131   test-error:0.341131+0.129919 
## [2]  train-error:0.124941+0.033432   test-error:0.358674+0.105191 
## [3]  train-error:0.115932+0.045501   test-error:0.340156+0.092598 
## [4]  train-error:0.124941+0.033432   test-error:0.340156+0.092598 
## [5]  train-error:0.089616+0.051301   test-error:0.394737+0.098465


Training with gbtree

gbtree is the default booster for xgb.train.

bst_1 <- xgb.train(data = xgb_train_matrix, 
                   label = getinfo(xgb_train_matrix, "label"),
                   max.depth = 2, 
                   eta = 1, 
                   nthread = 4, 
                   nround = 50, # number of trees used for model building
                   watchlist = watchlist, 
                   objective = "binary:logistic")
## [1]  train-error:0.250000    test-error:0.347826 
## [2]  train-error:0.178571    test-error:0.304348 
## [3]  train-error:0.142857    test-error:0.347826 
## [4]  train-error:0.089286    test-error:0.260870 
## [5]  train-error:0.053571    test-error:0.347826 
## [6]  train-error:0.053571    test-error:0.391304 
## [7]  train-error:0.035714    test-error:0.347826 
## [8]  train-error:0.000000    test-error:0.391304 
## [9]  train-error:0.017857    test-error:0.391304 
## [10] train-error:0.000000    test-error:0.391304 
## [11] train-error:0.000000    test-error:0.391304 
## [12] train-error:0.000000    test-error:0.347826 
## [13] train-error:0.000000    test-error:0.347826 
## [14] train-error:0.000000    test-error:0.347826 
## [15] train-error:0.000000    test-error:0.347826 
## [16] train-error:0.000000    test-error:0.347826 
## [17] train-error:0.000000    test-error:0.347826 
## [18] train-error:0.000000    test-error:0.347826 
## [19] train-error:0.000000    test-error:0.347826 
## [20] train-error:0.000000    test-error:0.347826 
## [21] train-error:0.000000    test-error:0.347826 
## [22] train-error:0.000000    test-error:0.347826 
## [23] train-error:0.000000    test-error:0.347826 
## [24] train-error:0.000000    test-error:0.347826 
## [25] train-error:0.000000    test-error:0.347826 
## [26] train-error:0.000000    test-error:0.347826 
## [27] train-error:0.000000    test-error:0.347826 
## [28] train-error:0.000000    test-error:0.347826 
## [29] train-error:0.000000    test-error:0.347826 
## [30] train-error:0.000000    test-error:0.347826 
## [31] train-error:0.000000    test-error:0.347826 
## [32] train-error:0.000000    test-error:0.347826 
## [33] train-error:0.000000    test-error:0.347826 
## [34] train-error:0.000000    test-error:0.347826 
## [35] train-error:0.000000    test-error:0.347826 
## [36] train-error:0.000000    test-error:0.347826 
## [37] train-error:0.000000    test-error:0.347826 
## [38] train-error:0.000000    test-error:0.347826 
## [39] train-error:0.000000    test-error:0.347826 
## [40] train-error:0.000000    test-error:0.347826 
## [41] train-error:0.000000    test-error:0.347826 
## [42] train-error:0.000000    test-error:0.347826 
## [43] train-error:0.000000    test-error:0.347826 
## [44] train-error:0.000000    test-error:0.347826 
## [45] train-error:0.000000    test-error:0.347826 
## [46] train-error:0.000000    test-error:0.347826 
## [47] train-error:0.000000    test-error:0.347826 
## [48] train-error:0.000000    test-error:0.347826 
## [49] train-error:0.000000    test-error:0.347826 
## [50] train-error:0.000000    test-error:0.347826

Each feature is grouped by importance with k-means clustering. Gain is the improvement in accuracy that the addition of a feature brings to the branches it is on.

features = colnames(matrix_train)
importance_matrix_1 <- xgb.importance(features, model = bst_1)
print(importance_matrix_1)
##                   Feature       Gain      Cover  Frequency
## 1:                    age 0.54677689 0.45726557 0.43333333
## 2:  days_onset_to_outcome 0.14753698 0.11686537 0.11111111
## 3:            early_onset 0.06724531 0.07365416 0.07777778
## 4:          early_outcome 0.06679880 0.08869426 0.07777778
## 5:      province_Shanghai 0.05325295 0.04317030 0.04444444
## 6: days_onset_to_hospital 0.03883631 0.11300882 0.13333333
## 7:         province_other 0.03502971 0.03977657 0.04444444
## 8:               gender_f 0.02610959 0.02484170 0.01111111
## 9:               hospital 0.01841345 0.04272324 0.06666667
xgb.ggplot.importance(importance_matrix_1) +
  theme_minimal()

pred_1 <- predict(bst_1, xgb_test_matrix)

result_1 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction_p_death = round(pred_1, digits = 2),
                       prediction = as.integer(pred_1 > 0.5),
                       prediction_eval = ifelse(as.integer(pred_1 > 0.5) != label, "wrong", "correct"))
result_1
##     case_ID outcome label prediction_p_death prediction prediction_eval
## 1  case_123   Death     1               0.01          0           wrong
## 2  case_127 Recover     0               0.00          0         correct
## 3  case_128 Recover     0               0.96          1           wrong
## 4   case_14 Recover     0               0.76          1           wrong
## 5   case_19   Death     1               0.95          1         correct
## 6    case_2   Death     1               0.85          1         correct
## 7   case_20 Recover     0               0.98          1           wrong
## 8   case_21 Recover     0               0.06          0         correct
## 9   case_34   Death     1               1.00          1         correct
## 10  case_37 Recover     0               0.00          0         correct
## 11   case_5 Recover     0               0.03          0         correct
## 12  case_51 Recover     0               0.27          0         correct
## 13  case_55 Recover     0               0.11          0         correct
## 14   case_6   Death     1               0.26          0           wrong
## 15  case_61   Death     1               1.00          1         correct
## 16  case_65 Recover     0               0.00          0         correct
## 17  case_74 Recover     0               0.92          1           wrong
## 18  case_78   Death     1               0.03          0           wrong
## 19  case_79 Recover     0               0.00          0         correct
## 20   case_8   Death     1               0.06          0           wrong
## 21  case_87   Death     1               0.64          1         correct
## 22  case_91 Recover     0               0.00          0         correct
## 23  case_94 Recover     0               0.01          0         correct
err <- as.numeric(sum(as.integer(pred_1 > 0.5) != label))/length(label)
print(paste("test-error =", round(err, digits = 2)))
## [1] "test-error = 0.35"


Training with gblinear

bst_2 <- xgb.train(data = xgb_train_matrix, 
                   booster = "gblinear", 
                   label = getinfo(xgb_train_matrix, "label"),
                   max.depth = 2, 
                   eta = 1, 
                   nthread = 4, 
                   nround = 50, # number of trees used for model building
                   watchlist = watchlist, 
                   objective = "binary:logistic")
## [1]  train-error:0.339286    test-error:0.434783 
## [2]  train-error:0.303571    test-error:0.391304 
## [3]  train-error:0.285714    test-error:0.347826 
## [4]  train-error:0.303571    test-error:0.347826 
## [5]  train-error:0.285714    test-error:0.347826 
## [6]  train-error:0.285714    test-error:0.391304 
## [7]  train-error:0.285714    test-error:0.391304 
## [8]  train-error:0.267857    test-error:0.391304 
## [9]  train-error:0.250000    test-error:0.391304 
## [10] train-error:0.250000    test-error:0.478261 
## [11] train-error:0.250000    test-error:0.434783 
## [12] train-error:0.232143    test-error:0.434783 
## [13] train-error:0.232143    test-error:0.434783 
## [14] train-error:0.232143    test-error:0.434783 
## [15] train-error:0.232143    test-error:0.434783 
## [16] train-error:0.232143    test-error:0.434783 
## [17] train-error:0.232143    test-error:0.434783 
## [18] train-error:0.232143    test-error:0.434783 
## [19] train-error:0.232143    test-error:0.434783 
## [20] train-error:0.232143    test-error:0.434783 
## [21] train-error:0.214286    test-error:0.434783 
## [22] train-error:0.214286    test-error:0.434783 
## [23] train-error:0.214286    test-error:0.434783 
## [24] train-error:0.214286    test-error:0.434783 
## [25] train-error:0.214286    test-error:0.434783 
## [26] train-error:0.214286    test-error:0.434783 
## [27] train-error:0.214286    test-error:0.434783 
## [28] train-error:0.214286    test-error:0.434783 
## [29] train-error:0.214286    test-error:0.434783 
## [30] train-error:0.232143    test-error:0.434783 
## [31] train-error:0.214286    test-error:0.434783 
## [32] train-error:0.214286    test-error:0.434783 
## [33] train-error:0.214286    test-error:0.434783 
## [34] train-error:0.214286    test-error:0.434783 
## [35] train-error:0.214286    test-error:0.434783 
## [36] train-error:0.214286    test-error:0.434783 
## [37] train-error:0.214286    test-error:0.434783 
## [38] train-error:0.214286    test-error:0.434783 
## [39] train-error:0.214286    test-error:0.434783 
## [40] train-error:0.214286    test-error:0.434783 
## [41] train-error:0.214286    test-error:0.434783 
## [42] train-error:0.214286    test-error:0.434783 
## [43] train-error:0.214286    test-error:0.434783 
## [44] train-error:0.214286    test-error:0.434783 
## [45] train-error:0.214286    test-error:0.434783 
## [46] train-error:0.214286    test-error:0.434783 
## [47] train-error:0.214286    test-error:0.434783 
## [48] train-error:0.214286    test-error:0.434783 
## [49] train-error:0.214286    test-error:0.434783 
## [50] train-error:0.214286    test-error:0.434783

Each feature is grouped by importance with k-means clustering. Gain is the improvement in accuracy that the addition of a feature brings to the branches it is on.

features = colnames(matrix_train)
importance_matrix_2 <- xgb.importance(features, model = bst_2)
print(importance_matrix_2)
##                    Feature     Weight
##  1:                    age  0.0439012
##  2:               hospital -0.3131160
##  3:               gender_f  0.8067830
##  4:       province_Jiangsu -0.8807380
##  5:      province_Shanghai -0.4065330
##  6:      province_Zhejiang -0.1284590
##  7:         province_other  0.4510330
##  8:  days_onset_to_outcome  0.0209952
##  9: days_onset_to_hospital -0.0856504
## 10:            early_onset  0.6452190
## 11:          early_outcome  2.5134300
xgb.ggplot.importance(importance_matrix_2) +
  theme_minimal()

pred_2 <- predict(bst_2, xgb_test_matrix)

result_2 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction_p_death = round(pred_2, digits = 2),
                       prediction = as.integer(pred_2 > 0.5),
                       prediction_eval = ifelse(as.integer(pred_2 > 0.5) != label, "wrong", "correct"))
result_2
##     case_ID outcome label prediction_p_death prediction prediction_eval
## 1  case_123   Death     1               0.10          0           wrong
## 2  case_127 Recover     0               0.03          0         correct
## 3  case_128 Recover     0               0.13          0         correct
## 4   case_14 Recover     0               0.15          0         correct
## 5   case_19   Death     1               0.28          0           wrong
## 6    case_2   Death     1               0.28          0           wrong
## 7   case_20 Recover     0               0.80          1           wrong
## 8   case_21 Recover     0               0.79          1           wrong
## 9   case_34   Death     1               0.74          1         correct
## 10  case_37 Recover     0               0.06          0         correct
## 11   case_5 Recover     0               0.14          0         correct
## 12  case_51 Recover     0               0.77          1           wrong
## 13  case_55 Recover     0               0.23          0         correct
## 14   case_6   Death     1               0.47          0           wrong
## 15  case_61   Death     1               0.52          1         correct
## 16  case_65 Recover     0               0.03          0         correct
## 17  case_74 Recover     0               0.70          1           wrong
## 18  case_78   Death     1               0.22          0           wrong
## 19  case_79 Recover     0               0.06          0         correct
## 20   case_8   Death     1               0.36          0           wrong
## 21  case_87   Death     1               0.65          1         correct
## 22  case_91 Recover     0               0.05          0         correct
## 23  case_94 Recover     0               0.07          0         correct
err <- as.numeric(sum(as.integer(pred_2 > 0.5) != label))/length(label)
print(paste("test-error =", round(err, digits = 2)))
## [1] "test-error = 0.43"


caret

Extreme gradient boosting is also implemented in the caret package. Caret also provides options for preprocessing, of which I will compare a few.


No preprocessing

library(caret)

set.seed(27)
model_xgb_null <-train(outcome ~ .,
                 data=val_train_data,
                 method="xgbTree",
                 preProcess = NULL,
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
confusionMatrix(predict(model_xgb_null, val_test_data[, -1]), val_test_data$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Death Recover
##    Death       5       4
##    Recover     4      10
##                                           
##                Accuracy : 0.6522          
##                  95% CI : (0.4273, 0.8362)
##     No Information Rate : 0.6087          
##     P-Value [Acc > NIR] : 0.4216          
##                                           
##                   Kappa : 0.2698          
##  Mcnemar's Test P-Value : 1.0000          
##                                           
##             Sensitivity : 0.5556          
##             Specificity : 0.7143          
##          Pos Pred Value : 0.5556          
##          Neg Pred Value : 0.7143          
##              Prevalence : 0.3913          
##          Detection Rate : 0.2174          
##    Detection Prevalence : 0.3913          
##       Balanced Accuracy : 0.6349          
##                                           
##        'Positive' Class : Death           
## 


Scaling and centering

With this method the column variables are centered (subtracting the column mean from each value in a column) and standardized (dividing by the column standard deviation).

set.seed(27)
model_xgb_sc_cen <-train(outcome ~ .,
                 data=val_train_data,
                 method="xgbTree",
                 preProcess = c("scale", "center"),
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
confusionMatrix(predict(model_xgb_sc_cen, val_test_data[, -1]), val_test_data$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Death Recover
##    Death       5       4
##    Recover     4      10
##                                           
##                Accuracy : 0.6522          
##                  95% CI : (0.4273, 0.8362)
##     No Information Rate : 0.6087          
##     P-Value [Acc > NIR] : 0.4216          
##                                           
##                   Kappa : 0.2698          
##  Mcnemar's Test P-Value : 1.0000          
##                                           
##             Sensitivity : 0.5556          
##             Specificity : 0.7143          
##          Pos Pred Value : 0.5556          
##          Neg Pred Value : 0.7143          
##              Prevalence : 0.3913          
##          Detection Rate : 0.2174          
##    Detection Prevalence : 0.3913          
##       Balanced Accuracy : 0.6349          
##                                           
##        'Positive' Class : Death           
## 

Scaling and centering did not improve the predictions.

pred_3 <- predict(model_xgb_sc_cen, val_test_data[, -1])
pred_3b <- round(predict(model_xgb_sc_cen, val_test_data[, -1], type="prob"), digits = 2)

result_3 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction = pred_3,
                       pred_3b)
result_3$prediction_eval <- ifelse(result_3$prediction != result_3$outcome, "wrong", "correct")
result_3
##     case_ID outcome label prediction Death Recover prediction_eval
## 1  case_123   Death     1    Recover  0.04    0.96           wrong
## 2  case_127 Recover     0    Recover  0.04    0.96         correct
## 3  case_128 Recover     0      Death  0.81    0.19           wrong
## 4   case_14 Recover     0      Death  0.70    0.30           wrong
## 5   case_19   Death     1      Death  0.67    0.33         correct
## 6    case_2   Death     1      Death  0.66    0.34         correct
## 7   case_20 Recover     0      Death  0.76    0.24           wrong
## 8   case_21 Recover     0    Recover  0.17    0.83         correct
## 9   case_34   Death     1      Death  0.94    0.06         correct
## 10  case_37 Recover     0    Recover  0.02    0.98         correct
## 11   case_5 Recover     0    Recover  0.09    0.91         correct
## 12  case_51 Recover     0    Recover  0.41    0.59         correct
## 13  case_55 Recover     0    Recover  0.12    0.88         correct
## 14   case_6   Death     1    Recover  0.30    0.70           wrong
## 15  case_61   Death     1      Death  0.99    0.01         correct
## 16  case_65 Recover     0    Recover  0.01    0.99         correct
## 17  case_74 Recover     0      Death  0.85    0.15           wrong
## 18  case_78   Death     1    Recover  0.21    0.79           wrong
## 19  case_79 Recover     0    Recover  0.01    0.99         correct
## 20   case_8   Death     1    Recover  0.20    0.80           wrong
## 21  case_87   Death     1      Death  0.65    0.35         correct
## 22  case_91 Recover     0    Recover  0.01    0.99         correct
## 23  case_94 Recover     0    Recover  0.07    0.93         correct


Box-Cox transformation

The Box-Cox power transformation is used to normalize data.

set.seed(27)
model_xgb_BoxCox <-train(outcome ~ .,
                 data=val_train_data,
                 method="xgbTree",
                 preProcess = "BoxCox",
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
confusionMatrix(predict(model_xgb_BoxCox, val_test_data[, -1]), val_test_data$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Death Recover
##    Death       5       4
##    Recover     4      10
##                                           
##                Accuracy : 0.6522          
##                  95% CI : (0.4273, 0.8362)
##     No Information Rate : 0.6087          
##     P-Value [Acc > NIR] : 0.4216          
##                                           
##                   Kappa : 0.2698          
##  Mcnemar's Test P-Value : 1.0000          
##                                           
##             Sensitivity : 0.5556          
##             Specificity : 0.7143          
##          Pos Pred Value : 0.5556          
##          Neg Pred Value : 0.7143          
##              Prevalence : 0.3913          
##          Detection Rate : 0.2174          
##    Detection Prevalence : 0.3913          
##       Balanced Accuracy : 0.6349          
##                                           
##        'Positive' Class : Death           
## 

Box-Cox transformation did not improve the predictions either.

pred_4 <- predict(model_xgb_BoxCox, val_test_data[, -1])
pred_4b <- round(predict(model_xgb_BoxCox, val_test_data[, -1], type="prob"), digits = 2)

result_4 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction = pred_4,
                       pred_4b)
result_4$prediction_eval <- ifelse(result_4$prediction != result_4$outcome, "wrong", "correct")
result_4
##     case_ID outcome label prediction Death Recover prediction_eval
## 1  case_123   Death     1    Recover  0.01    0.99           wrong
## 2  case_127 Recover     0    Recover  0.01    0.99         correct
## 3  case_128 Recover     0      Death  0.88    0.12           wrong
## 4   case_14 Recover     0      Death  0.81    0.19           wrong
## 5   case_19   Death     1      Death  0.88    0.12         correct
## 6    case_2   Death     1      Death  0.79    0.21         correct
## 7   case_20 Recover     0      Death  0.84    0.16           wrong
## 8   case_21 Recover     0    Recover  0.07    0.93         correct
## 9   case_34   Death     1      Death  0.97    0.03         correct
## 10  case_37 Recover     0    Recover  0.01    0.99         correct
## 11   case_5 Recover     0    Recover  0.06    0.94         correct
## 12  case_51 Recover     0    Recover  0.27    0.73         correct
## 13  case_55 Recover     0    Recover  0.09    0.91         correct
## 14   case_6   Death     1    Recover  0.30    0.70           wrong
## 15  case_61   Death     1      Death  1.00    0.00         correct
## 16  case_65 Recover     0    Recover  0.00    1.00         correct
## 17  case_74 Recover     0      Death  0.90    0.10           wrong
## 18  case_78   Death     1    Recover  0.11    0.89           wrong
## 19  case_79 Recover     0    Recover  0.00    1.00         correct
## 20   case_8   Death     1    Recover  0.13    0.87           wrong
## 21  case_87   Death     1      Death  0.61    0.39         correct
## 22  case_91 Recover     0    Recover  0.00    1.00         correct
## 23  case_94 Recover     0    Recover  0.04    0.96         correct


Principal Component Analysis (PCA)

PCA is used for dimensionality reduction. When applied as a preprocessing method the number of features are reduced by using the eigenvectors of the covariance matrix.

set.seed(27)
model_xgb_pca <-train(outcome ~ .,
                 data=val_train_data,
                 method="xgbTree",
                 preProcess = "pca",
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
confusionMatrix(predict(model_xgb_pca, val_test_data[, -1]), val_test_data$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Death Recover
##    Death       4       2
##    Recover     5      12
##                                           
##                Accuracy : 0.6957          
##                  95% CI : (0.4708, 0.8679)
##     No Information Rate : 0.6087          
##     P-Value [Acc > NIR] : 0.2644          
##                                           
##                   Kappa : 0.3207          
##  Mcnemar's Test P-Value : 0.4497          
##                                           
##             Sensitivity : 0.4444          
##             Specificity : 0.8571          
##          Pos Pred Value : 0.6667          
##          Neg Pred Value : 0.7059          
##              Prevalence : 0.3913          
##          Detection Rate : 0.1739          
##    Detection Prevalence : 0.2609          
##       Balanced Accuracy : 0.6508          
##                                           
##        'Positive' Class : Death           
## 

PCA did improve the predictions slightly.

pred_5 <- predict(model_xgb_pca, val_test_data[, -1])
pred_5b <- round(predict(model_xgb_pca, val_test_data[, -1], type="prob"), digits = 2)

result_5 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction = pred_5,
                       pred_5b)
result_5$prediction_eval <- ifelse(result_5$prediction != result_5$outcome, "wrong", "correct")
result_5
##     case_ID outcome label prediction Death Recover prediction_eval
## 1  case_123   Death     1    Recover  0.02    0.98           wrong
## 2  case_127 Recover     0    Recover  0.29    0.71         correct
## 3  case_128 Recover     0    Recover  0.01    0.99         correct
## 4   case_14 Recover     0    Recover  0.44    0.56         correct
## 5   case_19   Death     1    Recover  0.07    0.93           wrong
## 6    case_2   Death     1    Recover  0.26    0.74           wrong
## 7   case_20 Recover     0    Recover  0.16    0.84         correct
## 8   case_21 Recover     0      Death  0.57    0.43           wrong
## 9   case_34   Death     1      Death  0.93    0.07         correct
## 10  case_37 Recover     0    Recover  0.03    0.97         correct
## 11   case_5 Recover     0    Recover  0.11    0.89         correct
## 12  case_51 Recover     0    Recover  0.19    0.81         correct
## 13  case_55 Recover     0    Recover  0.16    0.84         correct
## 14   case_6   Death     1      Death  0.55    0.45         correct
## 15  case_61   Death     1    Recover  0.22    0.78           wrong
## 16  case_65 Recover     0    Recover  0.08    0.92         correct
## 17  case_74 Recover     0      Death  0.67    0.33           wrong
## 18  case_78   Death     1      Death  0.87    0.13         correct
## 19  case_79 Recover     0    Recover  0.16    0.84         correct
## 20   case_8   Death     1    Recover  0.22    0.78           wrong
## 21  case_87   Death     1      Death  0.66    0.34         correct
## 22  case_91 Recover     0    Recover  0.01    0.99         correct
## 23  case_94 Recover     0    Recover  0.14    0.86         correct


Median imputation

set.seed(27)
model_xgb_medianImpute <-train(outcome ~ .,
                 data=val_train_data,
                 method="xgbTree",
                 preProcess = "medianImpute",
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
confusionMatrix(predict(model_xgb_medianImpute, val_test_data[, -1]), val_test_data$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Death Recover
##    Death       5       4
##    Recover     4      10
##                                           
##                Accuracy : 0.6522          
##                  95% CI : (0.4273, 0.8362)
##     No Information Rate : 0.6087          
##     P-Value [Acc > NIR] : 0.4216          
##                                           
##                   Kappa : 0.2698          
##  Mcnemar's Test P-Value : 1.0000          
##                                           
##             Sensitivity : 0.5556          
##             Specificity : 0.7143          
##          Pos Pred Value : 0.5556          
##          Neg Pred Value : 0.7143          
##              Prevalence : 0.3913          
##          Detection Rate : 0.2174          
##    Detection Prevalence : 0.3913          
##       Balanced Accuracy : 0.6349          
##                                           
##        'Positive' Class : Death           
## 

Median imputation did not improve the predictions either.

pred_6 <- predict(model_xgb_medianImpute, val_test_data[, -1])
pred_6b <- round(predict(model_xgb_medianImpute, val_test_data[, -1], type="prob"), digits = 2)

result_6 <- data.frame(case_ID = rownames(val_test_data),
                       outcome = val_test_data$outcome, 
                       label = label, 
                       prediction = pred_6,
                       pred_6b)
result_6$prediction_eval <- ifelse(result_6$prediction != result_6$outcome, "wrong", "correct")
result_6
##     case_ID outcome label prediction Death Recover prediction_eval
## 1  case_123   Death     1    Recover  0.04    0.96           wrong
## 2  case_127 Recover     0    Recover  0.04    0.96         correct
## 3  case_128 Recover     0      Death  0.81    0.19           wrong
## 4   case_14 Recover     0      Death  0.70    0.30           wrong
## 5   case_19   Death     1      Death  0.67    0.33         correct
## 6    case_2   Death     1      Death  0.66    0.34         correct
## 7   case_20 Recover     0      Death  0.76    0.24           wrong
## 8   case_21 Recover     0    Recover  0.17    0.83         correct
## 9   case_34   Death     1      Death  0.94    0.06         correct
## 10  case_37 Recover     0    Recover  0.02    0.98         correct
## 11   case_5 Recover     0    Recover  0.09    0.91         correct
## 12  case_51 Recover     0    Recover  0.41    0.59         correct
## 13  case_55 Recover     0    Recover  0.12    0.88         correct
## 14   case_6   Death     1    Recover  0.30    0.70           wrong
## 15  case_61   Death     1      Death  0.99    0.01         correct
## 16  case_65 Recover     0    Recover  0.01    0.99         correct
## 17  case_74 Recover     0      Death  0.85    0.15           wrong
## 18  case_78   Death     1    Recover  0.21    0.79           wrong
## 19  case_79 Recover     0    Recover  0.01    0.99         correct
## 20   case_8   Death     1    Recover  0.20    0.80           wrong
## 21  case_87   Death     1      Death  0.65    0.35         correct
## 22  case_91 Recover     0    Recover  0.01    0.99         correct
## 23  case_94 Recover     0    Recover  0.07    0.93         correct


Comparison of extreme gradient boosting models

Combining results

library(dplyr)
result <- left_join(result_1[, c(1, 2, 6)], result_2[, c(1, 6)], by = "case_ID")
result <- left_join(result, result_3[, c(1, 7)], by = "case_ID")
result <- left_join(result, result_4[, c(1, 7)], by = "case_ID")
result <- left_join(result, result_5[, c(1, 7)], by = "case_ID")
result <- left_join(result, result_6[, c(1, 7)], by = "case_ID")
colnames(result)[-c(1:2)] <- c("pred_xgboost_gbtree", "pred_xgboost_gblinear", "model_xgb_sc_cen", "model_xgb_BoxCox", "pred_xgbTree_pca", "model_xgb_medianImpute")


What’s the rate of correctly predicted cases in the validation data?

round(sum(result$pred_xgboost_gbtree == "correct")/nrow(result), digits = 2)
## [1] 0.65
round(sum(result$pred_xgboost_gblinear == "correct")/nrow(result), digits = 2)
## [1] 0.57
round(sum(result$model_xgb_sc_cen == "correct")/nrow(result), digits = 2)
## [1] 0.65
round(sum(result$model_xgb_BoxCox == "correct")/nrow(result), digits = 2)
## [1] 0.65
round(sum(result$pred_xgbTree_pca == "correct")/nrow(result), digits = 2)
## [1] 0.7
round(sum(result$model_xgb_medianImpute == "correct")/nrow(result), digits = 2)
## [1] 0.65

PCA preprocessing combined with caret’s implementation of extreme gradient boosting had the highest number of accurately predicted cases in the validation data set.


Predicting unknown output

Because PCA preprocessing and caret’s implementation of extreme gradient boosting had the highest prediction accuracy, this model will be used to predict the unknown cases from the original dataset.

Uncertainty of prediction is accounted for with a cutoff below 0.7.

set.seed(27)
model_xgb_pca <-train(outcome ~ .,
                 data = train_data,
                 method = "xgbTree",
                 preProcess = "pca",
                 trControl = trainControl(method = "repeatedcv", number = 5, repeats = 10, verboseIter = FALSE))
pred <- predict(model_xgb_pca, test_data)
predb <- round(predict(model_xgb_pca, test_data, type="prob"), digits = 2)

result <- data.frame(case_ID = rownames(test_data),
                       prediction = pred,
                       predb)
result$predicted_outcome <- ifelse(result$Death > 0.7, "Death",
                                   ifelse(result$Recover > 0.7, "Recover", "uncertain"))
result
##     case_ID prediction Death Recover predicted_outcome
## 1  case_100    Recover  0.08    0.92           Recover
## 2  case_101    Recover  0.11    0.89           Recover
## 3  case_102    Recover  0.04    0.96           Recover
## 4  case_103    Recover  0.01    0.99           Recover
## 5  case_104    Recover  0.19    0.81           Recover
## 6  case_105    Recover  0.01    0.99           Recover
## 7  case_108      Death  0.99    0.01             Death
## 8  case_109    Recover  0.30    0.70         uncertain
## 9  case_110    Recover  0.16    0.84           Recover
## 10 case_112      Death  0.54    0.46         uncertain
## 11 case_113    Recover  0.02    0.98           Recover
## 12 case_114    Recover  0.30    0.70         uncertain
## 13 case_115    Recover  0.02    0.98           Recover
## 14 case_118      Death  0.65    0.35         uncertain
## 15 case_120      Death  0.89    0.11             Death
## 16 case_122    Recover  0.00    1.00           Recover
## 17 case_126    Recover  0.04    0.96           Recover
## 18 case_130    Recover  0.42    0.58         uncertain
## 19 case_132    Recover  0.47    0.53         uncertain
## 20 case_136      Death  0.86    0.14             Death
## 21  case_15      Death  0.99    0.01             Death
## 22  case_16      Death  0.91    0.09             Death
## 23  case_22      Death  0.80    0.20             Death
## 24  case_28      Death  0.91    0.09             Death
## 25  case_31      Death  1.00    0.00             Death
## 26  case_32      Death  0.93    0.07             Death
## 27  case_38      Death  0.87    0.13             Death
## 28  case_39      Death  0.88    0.12             Death
## 29   case_4      Death  0.99    0.01             Death
## 30  case_40      Death  1.00    0.00             Death
## 31  case_41      Death  0.95    0.05             Death
## 32  case_42    Recover  0.49    0.51         uncertain
## 33  case_47      Death  0.91    0.09             Death
## 34  case_48      Death  0.95    0.05             Death
## 35  case_52      Death  0.62    0.38         uncertain
## 36  case_54    Recover  0.43    0.57         uncertain
## 37  case_56      Death  0.81    0.19             Death
## 38  case_62      Death  0.58    0.42         uncertain
## 39  case_63    Recover  0.44    0.56         uncertain
## 40  case_66    Recover  0.18    0.82           Recover
## 41  case_67    Recover  0.02    0.98           Recover
## 42  case_68    Recover  0.27    0.73           Recover
## 43  case_69    Recover  0.35    0.65         uncertain
## 44  case_70    Recover  0.02    0.98           Recover
## 45  case_71    Recover  0.17    0.83           Recover
## 46  case_80    Recover  0.05    0.95           Recover
## 47  case_84    Recover  0.02    0.98           Recover
## 48  case_85      Death  0.99    0.01             Death
## 49  case_86    Recover  0.02    0.98           Recover
## 50  case_88    Recover  0.02    0.98           Recover
## 51   case_9      Death  0.98    0.02             Death
## 52  case_90    Recover  0.08    0.92           Recover
## 53  case_92    Recover  0.04    0.96           Recover
## 54  case_93    Recover  0.43    0.57         uncertain
## 55  case_95    Recover  0.11    0.89           Recover
## 56  case_96      Death  0.68    0.32         uncertain
## 57  case_99    Recover  0.04    0.96           Recover


Comparison with predicted outcome from last week’s analyses


Plotting predicted outcome

Results from last week’s analysis come from the following part of the analysis:

results <- data.frame(randomForest = predict(model_rf, newdata = test_data, type="prob"),
glmnet = predict(model_glmnet, newdata = test_data, type="prob"),
kknn = predict(model_kknn, newdata = test_data, type="prob"),
pda = predict(model_pda, newdata = test_data, type="prob"),
slda = predict(model_slda, newdata = test_data, type="prob"),
pam = predict(model_pam, newdata = test_data, type="prob"),
C5.0Tree = predict(model_C5.0Tree, newdata = test_data, type="prob"),
pls = predict(model_pls, newdata = test_data, type="prob"))

results$sum_Death <- rowSums(results[, grep("Death", colnames(results))])
results$sum_Recover <- rowSums(results[, grep("Recover", colnames(results))])
results$log2_ratio <- log2(results$sum_Recover/results$sum_Death)
results$predicted_outcome <- ifelse(results$log2_ratio > 1.5, "Recover", ifelse(results$log2_ratio < -1.5, "Death", "uncertain"))
results[, -c(1:16)]

results <- results_last_week

Combining the table with predicted outcomes from this and last week with the original data.

result <- merge(result[, c(1, 5)], results_last_week[, ncol(results_last_week), drop = FALSE], by.x = "case_ID", by.y = "row.names")
colnames(result)[2:3] <- c("predicted_outcome_xgboost", "predicted_outcome_last_week")

results_combined <- merge(result, fluH7N9.china.2013[which(fluH7N9.china.2013$case.ID %in% result$case_ID), ], 
                          by.x = "case_ID", by.y = "case.ID")
results_combined <- results_combined[, -c(6,7)]

For plotting with ggplot2, the dataframe needs to be gathered.

library(tidyr)
results_combined_gather <- results_combined %>%
  gather(group_dates, date, date.of.onset:date.of.hospitalisation)

results_combined_gather$group_dates <- factor(results_combined_gather$group_dates, levels = c("date.of.onset", "date.of.hospitalisation"))

results_combined_gather$group_dates <- mapvalues(results_combined_gather$group_dates, from = c("date.of.onset", "date.of.hospitalisation"), 
                                             to = c("Date of onset", "Date of hospitalisation"))

results_combined_gather$gender <- mapvalues(results_combined_gather$gender, from = c("f", "m"), 
                                             to = c("Female", "Male"))
levels(results_combined_gather$gender) <- c(levels(results_combined_gather$gender), "unknown")
results_combined_gather$gender[is.na(results_combined_gather$gender)] <- "unknown"

results_combined_gather <- results_combined_gather %>%
  gather(group_pred, prediction, predicted_outcome_xgboost:predicted_outcome_last_week)

results_combined_gather$group_pred <- mapvalues(results_combined_gather$group_pred, from = c("predicted_outcome_xgboost", "predicted_outcome_last_week"), 
                                             to = c("Predicted outcome from XGBoost", "Predicted outcome from last week"))

Setting a custom theme for plotting:

library(ggplot2)
my_theme <- function(base_size = 12, base_family = "sans"){
  theme_minimal(base_size = base_size, base_family = base_family) +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.5),
    axis.title = element_text(size = 14),
    panel.grid.major = element_line(color = "grey"),
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = "aliceblue"),
    strip.background = element_rect(fill = "lightgrey", color = "grey", size = 1),
    strip.text = element_text(face = "bold", size = 12, color = "black"),
    legend.position = "bottom",
    legend.justification = "top", 
    legend.box = "horizontal",
    legend.box.background = element_rect(colour = "grey50"),
    legend.background = element_blank(),
    panel.border = element_rect(color = "grey", fill = NA, size = 0.5),
    panel.spacing = unit(1, "lines")
  )
}
results_combined_gather$province <- mapvalues(results_combined_gather$province, 
                                                from = c("Anhui", "Beijing", "Fujian", "Guangdong", "Hebei", "Henan", "Hunan", "Jiangxi", "Shandong", "Taiwan"), 
                                                to = rep("Other", 10))

levels(results_combined_gather$gender) <- c(levels(results_combined_gather$gender), "unknown")
results_combined_gather$gender[is.na(results_combined_gather$gender)] <- "unknown"

results_combined_gather$province <- factor(results_combined_gather$province, levels = c("Jiangsu",  "Shanghai", "Zhejiang", "Other"))
ggplot(data = subset(results_combined_gather, group_dates == "Date of onset"), aes(x = date, y = as.numeric(age), fill = prediction)) +
  stat_density2d(aes(alpha = ..level..), geom = "polygon") +
  geom_jitter(aes(color = prediction, shape = gender), size = 2) +
  geom_rug(aes(color = prediction)) +
  labs(
    fill = "Predicted outcome",
    color = "Predicted outcome",
    alpha = "Level",
    shape = "Gender",
    x = "Date of onset in 2013",
    y = "Age",
    title = "2013 Influenza A H7N9 cases in China",
    subtitle = "Predicted outcome of cases with unknown outcome",
    caption = ""
  ) +
  facet_grid(group_pred ~ province) +
  my_theme() +
  scale_shape_manual(values = c(15, 16, 17)) +
  scale_color_brewer(palette="Set1", na.value = "grey50") +
  scale_fill_brewer(palette="Set1")


There is much less uncertainty in the XGBoost data, even tough I used slightly different methods for classifying uncertainty: In last week’s analysis I based uncertainty on the ratio of combined prediction values from all analyses, this week uncertainty is based on the prediction value from one analysis.


For a comparison of feature selection methods see here.


If you are interested in more machine learning posts, check out the category listing for machine_learning.


sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
## 
## locale:
## [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                           LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidyr_0.6.0     plyr_1.8.4      xgboost_0.6-0   caret_6.0-73    ggplot2_2.2.0   lattice_0.20-34 mice_2.25       Rcpp_0.12.8     dplyr_0.5.0     outbreaks_1.0.0
## 
## loaded via a namespace (and not attached):
##  [1] Ckmeans.1d.dp_3.4.6-4 RColorBrewer_1.1-2    compiler_3.3.2        nloptr_1.0.4          class_7.3-14          iterators_1.0.8       tools_3.3.2           rpart_4.1-10          digest_0.6.10         lme4_1.1-12           evaluate_0.10         tibble_1.2            gtable_0.2.0          nlme_3.1-128          mgcv_1.8-16           Matrix_1.2-7.1        foreach_1.4.3         DBI_0.5-1             parallel_3.3.2        yaml_2.1.14           SparseM_1.74          e1071_1.6-7           stringr_1.1.0         knitr_1.15            MatrixModels_0.4-1    stats4_3.3.2          grid_3.3.2            nnet_7.3-12           data.table_1.9.6      R6_2.2.0              survival_2.40-1       rmarkdown_1.1         minqa_1.2.4           reshape2_1.4.2        car_2.1-3             magrittr_1.5          scales_0.4.1          codetools_0.2-15      ModelMetrics_1.1.0    htmltools_0.3.5       MASS_7.3-45           splines_3.3.2         assertthat_0.1        pbkrtest_0.4-6        colorspace_1.3-0      labeling_0.3          quantreg_5.29         stringi_1.1.2         lazyeval_0.2.0        munsell_0.4.3         chron_2.3-47


blog comments powered by Disqus