必要なパッケージを読み込む。
::p_load(tidyverse,
pacman
ISLR,
rsample,
gmodels,
rpart,
rattle,
rpart.plot,
C50,
randomForest, ROCR)
次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
<- "Yu Gothic"
my_font else {
} <- "Japan1"
my_font
}else if (capabilities("aqua")) {
} <- "HiraginoSans-W3"
my_font else {
} <- "IPAexGothic"
my_font
}
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
を作る。
<- Carseats %>%
cs 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)
<- initial_split(cs)
cs_split <- training(cs_split)
cs_train <- testing(cs_split) cs_test
ジニ不純度(ジニ係数)に基づく決定木による分類は、rpartパッケージのrpart()
で実行できる。 学習データで学習する。
<- rpart(High ~ .,
tree1 data = cs_train,
method = "class")
学習結果は、rattle パッケージの fancyRpartPlot()
で可視化することができる。
fancyRpartPlot(tree1, caption = NULL)
この図から、どのような木が育てられたかを把握することができる。たとえば、最初の分岐では、ShelveLoc
の値が Bad
または Medium
のものが左下の枝に、それ以外(値が Good
)が右下の枝に分けられている。
学習したモデルによって、検証データの分類を行おう。 分類結果としてカテゴリ(今回の例では “Yes” または “No”)を得たいときは type = "class"
を、カテゴリに属する確率を得たいときは type = "prob"
を指定する。
<- predict(tree1,
pred1 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 |
## -------------|-----------|-----------|-----------|
##
##
この結果から、分類精度を評価するためのさまざまな指標を計算することができる(今回は省略するが、各自で計算して確かめること)。
決定木による分類で、どの説明変数が相対的に重要だったかを確認する。
$variable.importance tree1
## 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()
を使う。
<- prune(tree1, cp = 0.13) tree1_prune
剪定結果を可視化する。
fancyRpartPlot(tree1_prune, caption = NULL)
木が刈り込まれ、シンプルなモデルになった。ここで学習したモデルによって、検証データの分類を行おう。
<- predict(tree1_prune, type = "class", newdata = cs_test) pred1_prune
この結果を利用して混同行列を作る。
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()
で実行できる。 まず、学習データで学習する。
<- C5.0(High ~ .,
tree2 data = cs_train)
学習結果を可視化する。
plot(tree2)
ちょっとつらい……
学習したモデルを使って検証データの分類を行う。
<- predict(tree2, type = "class", newdata = cs_test) pred2
この結果を使って混同行列を作る(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
を指定して実行してみよう。
<- C5.0Control(minCases = 10)
tree_control <- C5.0(High ~ .,
tree2_prune data = cs_train,
control = tree_control)
結果を可視化する。
plot(tree2_prune)
学習したモデルを使って検証データの分類を行う。
<- predict(tree2_prune, type = "class", newdata = cs_test) pred2_prune
この結果を使って混同行列を作る。
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
で指定する。
ntree = 5
で学習を実行する。
<- randomForest(High ~ .,
rf5 data = cs_train,
ntree = 5)
学習したモデルを利用して、検証データの分類を行う。
<- predict(rf5, newdata = cs_test) pred_rf5
混同行列を作る。
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 |
## -------------|-----------|-----------|-----------|
##
##
ntree = 30
で学習を実行する。
<- randomForest(High ~ .,
rf30 data = cs_train,
ntree = 30)
学習したモデルを利用して、検証データの分類を行う。
<- predict(rf30, newdata = cs_test) pred_rf30
混同行列を作成する。
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の決定木
<- predict(tree1, type = "prob", newdata = cs_test)[, 2] %>%
roc_tree1 prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## C5.0 の決定木
<- predict(tree2, type = "prob", newdata = cs_test)[, 2] %>%
roc_tree2 prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## 木の数が5のランダムフォレスト
<- predict(rf5, type = "prob", newdata = cs_test)[, 2] %>%
roc_rf5 prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
## 木の数が30のランダムフォレスト
<- predict(rf30, type = "prob", newdata = cs_test)[, 2] %>%
roc_rf30 prediction(labels = cs_test$High) %>%
performance(measure = "tpr", x.measure = "fpr")
<- list(
fprs 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]]
)
<- list(
tprs 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]]
)
<- sapply(fprs, length)
roc_lengths
<- tibble(
df_rocs 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)"))
)
<- ggplot(df_rocs, aes(x = fpr, y = tpr, color = model)) +
ROCs geom_line() +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率(1 - 特異度)",
y = "真陽性率(感度)") +
scale_color_brewer(palette = "Set2",
name = "") +
coord_fixed()
plot(ROCs)
木の数が30のランダムフォレストの性能が最も良さそうだ。
AUCを計算してみよう。
<- predict(tree1, type = "prob", newdata = cs_test)[, 2] %>%
auc1 prediction(labels = cs_test$High) %>%
performance(measure = "auc")
<- predict(tree2, type = "prob", newdata = cs_test)[, 2] %>%
auc2 prediction(labels = cs_test$High) %>%
performance(measure = "auc")
<- predict(rf5, type = "prob", newdata = cs_test)[, 2] %>%
auc3 prediction(labels = cs_test$High) %>%
performance(measure = "auc")
<- predict(rf30, type = "prob", newdata = cs_test)[, 2] %>%
auc4 prediction(labels = cs_test$High) %>%
performance(measure = "auc")
cbind(auc1@y.values[[1]],
@y.values[[1]],
auc2@y.values[[1]],
auc3@y.values[[1]]) auc4
## [,1] [,2] [,3] [,4]
## [1,] 0.7706981 0.8273133 0.8112825 0.8780438
木の数が30のランダムフォレストのAUCが0.9を上回っており、性能が良いことがわかる。 反対に、rpart()
によるジニ不純度に基づく決定木の分類はAUCが0.8未満で、それほど良い分類ではなさそうだ。