library(mlbench)
data(PimaIndiansDiabetes)
df <- PimaIndiansDiabetes
tar.idx <- 9xgboost
パッケージの概要
xgboostパッケージは、勾配ブースティング決定木（Gradient Boosting Decision Tree, GBDT）を高速に実装したパッケージとなります。分類問題、回帰問題のどちらの問題にも対応可です。
使用例（分類問題）：PimaIndiansDiabetesデータの分類
PimaIndiansDiabetesデータを用いて、種々の項目からその患者が糖尿病患者であるかを判定するGBDTモデルを構築します。
データの確認
データの中身の確認です。（head）
head(df) pregnant glucose pressure triceps insulin mass pedigree age diabetes
1 6 148 72 35 0 33.6 0.627 50 pos
2 1 85 66 29 0 26.6 0.351 31 neg
3 8 183 64 0 0 23.3 0.672 32 pos
4 1 89 66 23 94 28.1 0.167 21 neg
5 0 137 40 35 168 43.1 2.288 33 pos
6 5 116 74 0 0 25.6 0.201 30 negデータの中身の確認です。（str）
str(df)'data.frame': 768 obs. of 9 variables:
$ pregnant: num 6 1 8 1 0 5 3 10 2 8 ...
$ glucose : num 148 85 183 89 137 116 78 115 197 125 ...
$ pressure: num 72 66 64 66 40 74 50 0 70 96 ...
$ triceps : num 35 29 0 23 35 0 32 0 45 0 ...
$ insulin : num 0 0 0 94 168 0 88 0 543 0 ...
$ mass : num 33.6 26.6 23.3 28.1 43.1 25.6 31 35.3 30.5 0 ...
$ pedigree: num 0.627 0.351 0.672 0.167 2.288 ...
$ age : num 50 31 32 21 33 30 26 29 53 54 ...
$ diabetes: Factor w/ 2 levels "neg","pos": 2 1 2 1 2 1 2 1 2 2 ...各変数の意味は以下の通りです。
| 変数名 | データ型 | 概要 |
|---|---|---|
| pregnant | num | 妊娠回数 |
| glucose | num | 血漿グルコース濃度（負荷試験） |
| pressure | num | 拡張期血圧 [mm Hg] |
| triceps | num | 上腕三頭筋皮膚襞の厚さ [mm ] |
| insulin | num | 2時間血清インスリン [mu U/ml] |
| mass | num | BMI |
| pedigree | num | 糖尿病血統機能 |
| age | num | 年齢 |
| diabetes | Factor | クラス変数（糖尿病の検査） |
summaryを表示します。
summary(df) pregnant glucose pressure triceps
Min. : 0.000 Min. : 0.0 Min. : 0.00 Min. : 0.00
1st Qu.: 1.000 1st Qu.: 99.0 1st Qu.: 62.00 1st Qu.: 0.00
Median : 3.000 Median :117.0 Median : 72.00 Median :23.00
Mean : 3.845 Mean :120.9 Mean : 69.11 Mean :20.54
3rd Qu.: 6.000 3rd Qu.:140.2 3rd Qu.: 80.00 3rd Qu.:32.00
Max. :17.000 Max. :199.0 Max. :122.00 Max. :99.00
insulin mass pedigree age diabetes
Min. : 0.0 Min. : 0.00 Min. :0.0780 Min. :21.00 neg:500
1st Qu.: 0.0 1st Qu.:27.30 1st Qu.:0.2437 1st Qu.:24.00 pos:268
Median : 30.5 Median :32.00 Median :0.3725 Median :29.00
Mean : 79.8 Mean :31.99 Mean :0.4719 Mean :33.24
3rd Qu.:127.2 3rd Qu.:36.60 3rd Qu.:0.6262 3rd Qu.:41.00
Max. :846.0 Max. :67.10 Max. :2.4200 Max. :81.00 その他のPimaIndiansDiabetesデータのEDA（探索的データ解析）については他のパッケージのコード活用例を参照ください。
各種手法計算用の関数を定義します。
create_df_res <- function(cmat){
TN <- cmat[1,1]
FP <- cmat[1,2]
FN <- cmat[2,1]
TP <- cmat[2,2]
# 正解率：どれだけ正しく分類できたかの割合
acc <- round((TP + TN) / (TP + TN + FN + FP), 3)
# 適合率：陽性と判定されたものがどれだけ正しく陽性であるかの割合
prec <- round(TP / (TP + FP), 3)
# 再現率（真陽性率）：実際に陽性のものをどれだけ正しく陽性と判定できたかの割合
rec <- round(TP / (TP + FN), 3)
# F値：適合率と再現率の調和平均（両者はトレードオフの関係）
Fval <- round(2 * rec * prec / (rec + prec), 3)
# 真陰性率：実際に陰性のものをどれだけ正しく陰性と判定できたかの割合
TNRat <- round(TN / (TN + FP), 3)
df_rat <- data.frame(
Item = c("正解率", "適合率", "再現率（真陽性率）", "F値", "真陰性率"),
Rate = c(acc, prec, rec, Fval, TNRat)
)
df_rat
}data.frameの出力整形用の関数を定義します。
library(knitr)
library(kableExtra)
show_df <- function(df, cap=""){
if(is_html_output()){
kable(df, format = "html", caption = cap) %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "left")
}else{
kable(df, format = "pipe", caption = cap)
}
}PimaIndiansDiabetesデータを、モデル生成のための訓練データとモデル評価のためのテストデータに分割します。データ割合は、訓練データ7割・テストデータ3割とします。確認のためにデータサイズを出力します。
# 再現性のためにシードを設定
set.seed(123)
# データの分割
sample_indices <- sample(1:nrow(df), 0.7 * nrow(df))
train_data <- df[sample_indices, ]
test_data <- df[-sample_indices, ]
# データサイズの確認
c(nrow(df), nrow(train_data), nrow(test_data))[1] 768 537 231モデルの構築
訓練データを用いてGBDTモデルを構築します。
library(xgboost)データを準備します。データはxgboostで使用されるxgb.DMatrix型とします。全てのデータは数値でなくてはならない点にご注意ください。因子型（Factor）もNGとなります。目的変数は1:陰性・2:陽性ですが、扱いやすくするため-1して0:陰性・1:陽性として扱うことにします。
train_X <- as.matrix(train_data[-tar.idx])
train_y <- as.numeric(train_data[[tar.idx]]) - 1
dtrain <- xgb.DMatrix(data = train_X, label = train_y)モデルを設定します。目的関数はparamsのobjectiveとして設定します。今回は二値の分類問題であるので”binary:logistic”を使用します。
model <- xgb.train(
data = dtrain,
nrounds = 100,
params = list(objective = "binary:logistic", eval_metric = "logloss")
)テストデータでの結果を確認します。こちらもまずはデータを準備します。
test_X <- as.matrix(test_data[-tar.idx])
res_act <- as.numeric(test_data[[tar.idx]]) - 1予測結果を変数に格納します。結果は数値（小数）になっているため、分類（0 or 1）に変更します。
res_pred <- predict(model, as.matrix(test_X))
res_pred <- ifelse(res_pred > 0.5, 1, 0)実際の結果と予測結果を比較できる混同行列を出力します。
(conf_mat <- table(res_act, res_pred)) res_pred
res_act 0 1
0 124 26
1 32 49評価指標を確認します。他のコード活用例との比較として、tuneRangerを用いてチューニングしたmlrの分類木において正解率0.719・F値0.545でしたので、高いモデル性能であることが分かります。
df_res <- create_df_res (conf_mat)
show_df(df_res)| Item | Rate |
|---|---|
| 正解率 | 0.749 |
| 適合率 | 0.653 |
| 再現率（真陽性率） | 0.605 |
| F値 | 0.628 |
| 真陰性率 | 0.827 |
ハイパーパラメータのチューニング
xgboostにおいてGBDT構築に使用される主なハイパーパラメータは以下の通りです。
| パラメータ名 | 概要 |
|---|---|
| eta | 学習率（η） |
| nrounds | 作成する木の本数 |
| min_child_weight | 葉に必要な最小重み |
| max_depth | 木の高さ |
| gamma | 木の分割時に必要な最小改善量 |
このうち、max_depthの違いによるモデルの差を見てゆきます。まずはmax_depthを引数として訓練、性能検証を一度にやってくれる関数を作成します。
train_test.max_depth <- function(max_depth_){
model <- xgb.train(
data = dtrain,
nrounds = 100,
params = list(max_depth = max_depth_, objective = "binary:logistic", eval_metric = "logloss")
)
res_pred <- predict(model, as.matrix(test_X))
res_pred <- ifelse(res_pred > 0.5, 1, 0)
conf_mat <- table(res_act, res_pred)
df_res <- create_df_res (conf_mat)
df_res
}max_depthを4, 6, 8, 10と変えたときのモデル性能を比較します。max_depthが8のときに正解率やF値が最も高くなっています。
n_loop <- 4
n_row <- 5
df_list <- vector("list", n_loop)
for (i in seq_len(n_loop)) {
depth <- 4 + 2 * (i - 1)
df_res <- train_test.max_depth(depth)
df_list[[i]] <- df_res
colnames(df_list[[i]])[2] <- paste0("depth = ", depth)
}
df_all <- Reduce(
function(x, y) merge(x, y, by = "Item", all = TRUE, sort = FALSE),
df_list
)
show_df(df_all)| Item | depth = 4 | depth = 6 | depth = 8 | depth = 10 |
|---|---|---|---|---|
| 正解率 | 0.736 | 0.749 | 0.762 | 0.749 |
| 適合率 | 0.632 | 0.653 | 0.659 | 0.642 |
| 再現率（真陽性率） | 0.593 | 0.605 | 0.667 | 0.642 |
| F値 | 0.612 | 0.628 | 0.663 | 0.642 |
| 真陰性率 | 0.813 | 0.827 | 0.813 | 0.807 |