情報科学3
決定木（分類木）
矢内 勇生
2021-10-31
準備
必要なパッケージを読み込む。
pacman::p_load(tidyverse,
ISLR,
rsample,
gmodels,
rpart,
rattle,
rpart.plot,
C50,
randomForest,
ROCR)次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
my_font <- "Yu Gothic"
} else {
my_font <- "Japan1"
}
} else if (capabilities("aqua")) {
my_font <- "HiraginoSans-W3"
} else {
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))データの準備
ISLRパッケージに含まれる Carseats データ (James, Witten, Hastie, Tibshirani 2018) を使って実習を行う。 このデータは、チャイルドシート (child car seats) の売上に関するデータである。
data(Carseats)このデータセットには Sales（売上）という変数がある。これは量的変数だが、分類の練習をするために Sales > 8 のものを売上が「高い」、それ以外を売上が「低い」ものとして扱う二値変数（ダミー変数）High を作る。
cs <- Carseats %>%
mutate(High = factor(ifelse(Sales > 8, 1, 0),
levels = 0:1,
labels = c("No", "Yes"))) %>%
select(!Sales) %>%
mutate_if(is.double, scale) # 量的変数を標準化するデータを訓練データと検証データに分割する。
set.seed(2021-10-31)
cs_split <- initial_split(cs)
cs_train <- training(cs_split)
cs_test <- testing(cs_split)木による分類
ジニ不純度に基づく分類
ジニ不純度（ジニ係数）に基づく決定木による分類は、rpartパッケージのrpart() で実行できる。 学習データで学習する。
tree1 <- rpart(High ~ .,
data = cs_train,
method = "class")学習結果は、rattle パッケージの fancyRpartPlot() で可視化することができる。
fancyRpartPlot(tree1, caption = NULL)
この図から、どのような木が育てられたかを把握することができる。たとえば、最初の分岐では、ShelveLoc の値が Bad または Medium のものが左下の枝に、それ以外（値が Good）が右下の枝に分けられている。
学習したモデルによって、検証データの分類を行おう。 分類結果としてカテゴリ（今回の例では “Yes” または “No”）を得たいときは type = "class" を、カテゴリに属する確率を得たいときは type = "prob" を指定する。
pred1 <- predict(tree1,
type = "class",
newdata = cs_test)この結果を利用して混同行列を作る。
CrossTable(x = cs_test$High,
y = pred1,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 48 | 8 | 56 |
## | 0.480 | 0.080 | |
## -------------|-----------|-----------|-----------|
## Yes | 18 | 26 | 44 |
## | 0.180 | 0.260 | |
## -------------|-----------|-----------|-----------|
## Column Total | 66 | 34 | 100 |
## -------------|-----------|-----------|-----------|
##
## この結果から、分類精度を評価するためのさまざまな指標を計算することができる（今回は省略するが、各自で計算して確かめること）。
決定木による分類で、どの説明変数が相対的に重要だったかを確認する。
tree1$variable.importance## ShelveLoc Price CompPrice Advertising US Age
## 27.124401 25.030883 13.568839 11.571902 6.351470 5.811129
## Income Population Education
## 5.489632 3.234744 2.446623チャイルドシートの売上が高いかどうかを分類するには、ShelveLoc （売り場の棚の位置）とPrice （価格）が重要であることがわかる。
rpart() は、木を育てる際に内部で交差検証（10分割交差検証、詳細は今後の授業で説明する）を行っている。 その結果は、printcp() で確認できる。
printcp(tree1)##
## Classification tree:
## rpart(formula = High ~ ., data = cs_train, method = "class")
##
## Variables actually used in tree construction:
## [1] Advertising CompPrice Income Price ShelveLoc US
##
## Root node error: 120/300 = 0.4
##
## n= 300
##
## CP nsplit rel error xerror xstd
## 1 0.308333 0 1.00000 1.00000 0.070711
## 2 0.054167 1 0.69167 0.69167 0.064569
## 3 0.041667 4 0.51667 0.70833 0.065041
## 4 0.025000 6 0.43333 0.67500 0.064080
## 5 0.016667 8 0.38333 0.69167 0.064569
## 6 0.010000 9 0.36667 0.64167 0.063046この結果の一番下に示されている表の読み方をごく簡単に説明する。まず、CP は complexity parameter（複雑度パラメタ）で、（名前とは裏腹に）値が小さいほどモデルが複雑になる。言い換えると、CP の値が小さいほど木が大きく育つ。rpart() のデフォルトでは、cp = 0.01 が採用されている。
nsplit は木が分岐する回数を示す。CP の値が小さいほど、分岐の回数が増え、木が大きく育っていることがわかる。上で可視化した木を見ると、9個の分岐点があることが確認できる。
rel error は、分岐がまったくない場合と比べて、誤分類がどの程度あるかを示している。nsplitが0のときのrel error は必ず1になる。木が育つほど、この値は小さくなる。xerror は、cross-validation error（交差検証誤差）で、交差検証による誤分類の割合が、分岐がまったくない木と比べてどれくらいの割合になるかを示している。nsplitが0のときのxerror も必ず1になる。xstd は xerror の標準偏差である。
これらの値を利用して、木を剪定するかどうかを判断する。慣例として、xerrorの最小値にその標準偏差1つ分を加えた値未満のxerror をもつ木の中で最小のものを選ぶ。この例では、xerror の最小値は0.64167で、その標準偏差は0.063046 である。\(0.64167 + 0.063046 \approx 0.70\) なので、xerror が0.70未満の木のなかで最も小さな木がほしい。この例では、そのような木の CP は0.054167 である。 これは、plotcp()で視覚的に確認できる。
plotcp(tree1)
この図を見ると、cp が0.13 のときに xerror (X-val Relative Error) が0.7 を下回っている。そこで、cp = 0.13 として木の剪定 (pruning) を行う。prune() を使う。
tree1_prune <- prune(tree1, cp = 0.13)剪定結果を可視化する。
fancyRpartPlot(tree1_prune, caption = NULL)
木が刈り込まれ、シンプルなモデルになった。ここで学習したモデルによって、検証データの分類を行おう。
pred1_prune <- predict(tree1_prune, type = "class", newdata = cs_test)この結果を利用して混同行列を作る。
CrossTable(x = cs_test$High,
y = pred1_prune,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 53 | 3 | 56 |
## | 0.530 | 0.030 | |
## -------------|-----------|-----------|-----------|
## Yes | 31 | 13 | 44 |
## | 0.310 | 0.130 | |
## -------------|-----------|-----------|-----------|
## Column Total | 84 | 16 | 100 |
## -------------|-----------|-----------|-----------|
##
## 剪定前と比べると、偽陽性が減った一方で、偽陰性が増えている。
エントロピーに基づく分類
エントロピー（情報量基準）に基づく分類は、C50パッケージのC5.0()で実行できる。 まず、学習データで学習する。
tree2 <- C5.0(High ~ .,
data = cs_train)学習結果を可視化する。
plot(tree2)
ちょっとつらい……
学習したモデルを使って検証データの分類を行う。
pred2 <- predict(tree2, type = "class", newdata = cs_test)この結果を使って混同行列を作る（summary(tree2) でも混同行列が表示される）。
CrossTable(x = cs_test$High,
y = pred2,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 50 | 6 | 56 |
## | 0.500 | 0.060 | |
## -------------|-----------|-----------|-----------|
## Yes | 17 | 27 | 44 |
## | 0.170 | 0.270 | |
## -------------|-----------|-----------|-----------|
## Column Total | 67 | 33 | 100 |
## -------------|-----------|-----------|-----------|
##
## この結果から、分類精度を評価するためのさまざまな指標を計算することができる（今回は省略するが、各自で計算して確かめること）。
木が成長し過ぎているかもしれないので、少し剪定する。 分岐先のケース数が10以上になるように minCases = 10 を指定して実行してみよう。
tree_control <- C5.0Control(minCases = 10)
tree2_prune <- C5.0(High ~ .,
data = cs_train,
control = tree_control)結果を可視化する。
plot(tree2_prune)
学習したモデルを使って検証データの分類を行う。
pred2_prune <- predict(tree2_prune, type = "class", newdata = cs_test)この結果を使って混同行列を作る。
CrossTable(x = cs_test$High,
y = pred2_prune,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 49 | 7 | 56 |
## | 0.490 | 0.070 | |
## -------------|-----------|-----------|-----------|
## Yes | 18 | 26 | 44 |
## | 0.180 | 0.260 | |
## -------------|-----------|-----------|-----------|
## Column Total | 67 | 33 | 100 |
## -------------|-----------|-----------|-----------|
##
## 木を単純化した結果、偽陽性も偽陰性も増えてしまった。
ランダムフォレスト
続いて、ランダムフォレストで分類を行ってみよう。randomForest パッケージの randomeForest() で分類することができる。木の数は ntree で指定する。
木の数が5のとき
ntree = 5 で学習を実行する。
rf5 <- randomForest(High ~ .,
data = cs_train,
ntree = 5)学習したモデルを利用して、検証データの分類を行う。
pred_rf5 <- predict(rf5, newdata = cs_test)混同行列を作る。
CrossTable(x = cs_test$High,
y = pred_rf5,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 49 | 7 | 56 |
## | 0.490 | 0.070 | |
## -------------|-----------|-----------|-----------|
## Yes | 21 | 23 | 44 |
## | 0.210 | 0.230 | |
## -------------|-----------|-----------|-----------|
## Column Total | 70 | 30 | 100 |
## -------------|-----------|-----------|-----------|
##
## 木の数が30のとき
ntree = 30 で学習を実行する。
rf30 <- randomForest(High ~ .,
data = cs_train,
ntree = 30)学習したモデルを利用して、検証データの分類を行う。
pred_rf30 <- predict(rf30, newdata = cs_test)混同行列を作成する。
CrossTable(x = cs_test$High,
y = pred_rf30,
dnn = c("実際「高い」？", "分類「高い」？"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 100
##
##
## | 分類「高い」？
## 実際「高い」？ | No | Yes | Row Total |
## -------------|-----------|-----------|-----------|
## No | 53 | 3 | 56 |
## | 0.530 | 0.030 | |
## -------------|-----------|-----------|-----------|
## Yes | 18 | 26 | 44 |
## | 0.180 | 0.260 | |
## -------------|-----------|-----------|-----------|
## Column Total | 71 | 29 | 100 |
## -------------|-----------|-----------|-----------|
##
## 分類の比較
ROC曲線を描いて比較する。
## rpartの決定木
roc_tree1 <- predict(tree1, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## C5.0 の決定木
roc_tree2 <- predict(tree2, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## 木の数が5のランダムフォレスト
roc_rf5 <- predict(rf5, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## 木の数が30のランダムフォレスト
roc_rf30 <- predict(rf30, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
fprs <- list(
fpr_1 = roc_tree1@x.values[[1]],
fpr_2 = roc_tree2@x.values[[1]],
fpr_3 = roc_rf5@x.values[[1]],
fpr_4 = roc_rf30@x.values[[1]]
)
tprs <- list(
tpr_tree1 = roc_tree1@y.values[[1]],
tpr_2 = roc_tree2@y.values[[1]],
tpr_3 = roc_rf5@y.values[[1]],
tpr_4 = roc_rf30@y.values[[1]]
)
roc_lengths <- sapply(fprs, length)
df_rocs <- tibble(
fpr = unlist(fprs),
tpr = unlist(tprs),
model = factor(rep(1:4, times = roc_lengths),
levels = 1:4,
labels = c("rpart", "C5.0",
"RF (n = 5)", "RF (n = 30)"))
)
ROCs <- ggplot(df_rocs, aes(x = fpr, y = tpr, color = model)) +
geom_line() +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率（1 - 特異度）",
y = "真陽性率（感度）") +
scale_color_brewer(palette = "Set2",
name = "") +
coord_fixed()
plot(ROCs)
木の数が30のランダムフォレストの性能が最も良さそうだ。
AUCを計算してみよう。
auc1 <- predict(tree1, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "auc")
auc2 <- predict(tree2, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "auc")
auc3 <- predict(rf5, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "auc")
auc4 <- predict(rf30, type = "prob", newdata = cs_test)[, 2] %>%
prediction(labels = cs_test$High) %>%
performance(measure = "auc")
cbind(auc1@y.values[[1]],
auc2@y.values[[1]],
auc3@y.values[[1]],
auc4@y.values[[1]])## [,1] [,2] [,3] [,4]
## [1,] 0.7706981 0.8273133 0.8112825 0.8780438木の数が30のランダムフォレストのAUCが0.9を上回っており、性能が良いことがわかる。 反対に、rpart()によるジニ不純度に基づく決定木の分類はAUCが0.8未満で、それほど良い分類ではなさそうだ。