情報科学3

モデルアンサンブル

準備

必要なパッケージを読み込む。

pacman::p_load(tidyverse,
               ISLR,
               rsample,
               gmodels,
               rpart,
               ipred,
               irr, 
               randomForest,
               adabag)

決定木 の実習でも用いた、 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-11-12)
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")

テストデータの類を行う。

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 |        43 |        14 |        57 | 
##              |     0.430 |     0.140 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        19 |        24 |        43 | 
##              |     0.190 |     0.240 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        62 |        38 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

この学習モデルの正解率は、

0.43 + 0.24

## [1] 0.67

である。

バギング

バギングを行う。仕組みを理解するために、自分でコードを書いて実装する。

まず、ブートストラップ標本を作る方法を確認しよう。 標本サイズ\(N\)の訓練データから、標本サイズ\(N\)のブートストラップ標本を抽出する。

N <- nrow(cs_train)
bs1_index <- sample(1:N, 
                    size = N, 
                    replace = TRUE)
bs1 <- cs_train[bs1_index, ]

中身を確認してみよう（結果は省略）。

glimpse(bs1)

ブートストラップ標本を繰り返し抽出し、そのそれぞれで決定木による分類学習を行う関数を作る。 nbagg 引数で、ブートストラップの回数を指定する。

mybag <- function(training, test, nbagg = 25) {
  pred <- list()
  N <- nrow(training)
  for (i in 1:nbagg) {
    bs_ind <- sample(1:N, 
                     size = N, 
                     replace = TRUE)
    bs <- training[bs_ind, ]
    tree <- rpart(High ~ ., 
                  data = bs,
                  method = "class")
    pred[[i]] <- predict(tree, 
                         type = "class", 
                         newdata = test)
  }
  return(pred)
}

この関数を使ってみよう。

bagging1 <- mybag(training = cs_train,
                  test = cs_test,
                  nbagg = 25)

この結果から1つのブートストラップ標本に基づく学習結果を取り出して混同行列を作ってみよう。

たとえば、5番目の結果は次のように取り出せる。

CrossTable(x = cs_test$High, 
           y = bagging1[[5]],
           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 |        45 |        12 |        57 | 
##              |     0.450 |     0.120 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        22 |        21 |        43 | 
##              |     0.220 |     0.210 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        67 |        33 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

同様に、12番目の結果を取り出してみよう。

CrossTable(x = cs_test$High, 
           y = bagging1[[12]],
           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 |        41 |        16 |        57 | 
##              |     0.410 |     0.160 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        16 |        27 |        43 | 
##              |     0.160 |     0.270 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        57 |        43 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

学習結果にばらつきがあることが確認できる。

これらの学習結果を組み合わせよう。ここでは、多数決によって分類を決定する。すなわち、25の学習モデルのうち、13個以上で “Yes” になったものを “Yes” に、それ以外を “No” に分類しよう。

そのために、まず結果をtibble に変換し、各行にIDを割り当てる

names(bagging1) <- paste0("model", 1:25)
bagging1_df <- bagging1 %>% 
  as_tibble() %>% 
  mutate(id = 1:n())

横長データを縦長に変換する。

bagging1_long <- bagging1_df %>% 
  pivot_longer(cols = model1:model25,
               values_to = "pred",
               names_to = "model",
               names_prefix = "model")

ID（テストデータの観測点）ごとに、多数決を取り、分類結果を決める。

bagging1_res <- bagging1_long %>% 
  group_by(id) %>% 
  summarize(pred = ifelse(sum(pred == "Yes") > 12, 1, 0),
            .groups = "drop") %>% 
  mutate(pred = factor(pred,
                       levels = 0:1,
                       labels = c("No", "Yes")))

バギングの結果を混同行列にする。

CrossTable(x = cs_test$High, 
           y = bagging1_res$pred,
           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 |        46 |        11 |        57 | 
##              |     0.460 |     0.110 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        18 |        25 |        43 | 
##              |     0.180 |     0.250 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        64 |        36 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

この学習モデルの正解率は、

0.46 + 0.25

## [1] 0.71

であり、バギングを使わない場合と比べて0.04ポイント性能が向上している。

バギング用いない場合と用いた場合でコーエンの\(\kappa\)を比較してみよう。

df_kappa <- tibble(
  observed = cs_test$High,
  pred1 = pred1,
  pred2 = bagging1_res$pred,
)
cbind(kappa2(df_kappa[, 1:2])$value,
      kappa2(df_kappa[,  -2])$value)

##          [,1]      [,2]
## [1,] 0.317053 0.3963364

\(\kappa\)で測っても、バギングによって性能が上がったことがわかる。

バギングを行うための関数として、ipredパッケージのbagging() がある。adabaggパッケージにも同じ名前の関数があるので、パッケージを明示的に指定する。

bagging2 <- ipred::bagging(High ~ .,
                           data = cs_train,
                           nbagg = 25)

テストデータの分類は、他の関数と同様に predict()で行える。

pred_bagg2 <- predict(bagging2,
                      newdata = cs_test)

このバギングの結果を混同行列にする。

CrossTable(x = cs_test$High, 
           y = pred_bagg2,
           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 |         9 |        57 | 
##              |     0.480 |     0.090 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        15 |        28 |        43 | 
##              |     0.150 |     0.280 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        63 |        37 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

上で行ったバギングの結果に近い結果が得られた。

ランダムフォレスト

木の数を50とするランダムフォレストで分類学習を実行する。randomForestパッケージのrandomForest() を使う。

rf1 <- randomForest(High ~ .,
                    data = cs_train,
                    num.trees = 50)

学習したモデルを利用して、テストデータの分類を行う。

pred_rf1 <- predict(rf1,
                    newdata = cs_test) 

混同行列を作成する。

CrossTable(x = cs_test$High, 
           y = pred_rf1,
           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 |         8 |        57 | 
##              |     0.490 |     0.080 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        16 |        27 |        43 | 
##              |     0.160 |     0.270 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        65 |        35 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

この学習モデルの正解率は、

0.49 + 0.27

## [1] 0.76

である。

コーエンの \(\kappa\) は、

df_kappa$rf1 <- pred_rf1
kappa2(df_kappa[,  c(1, 4)])$value

## [1] 0.4989562

なので、1つの決定木による分類よりは性能が良い。

（ランダムフォレストについては、第10回の実習資料 も参照。）

ブースティング

adabagパッケージのboosting() を使ってブースティングを実行しよう。

adab <- boosting(High ~ .,
                 data = cs_train)

テストデータの分類は、他の関数と同様に predict()で行える。 ただし、この結果には分類さらたクラス（カテゴリ）以外の情報を含まれる。分類の結果自体は class に保存されている。

pred_adab <- predict(adab, 
                     newdata = cs_test)

分類結果を混同行列にする。

CrossTable(x = cs_test$High, 
           y = pred_adab$class,
           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 |         9 |        57 | 
##              |     0.480 |     0.090 |           | 
## -------------|-----------|-----------|-----------|
##          Yes |        13 |        30 |        43 | 
##              |     0.130 |     0.300 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |        61 |        39 |       100 | 
## -------------|-----------|-----------|-----------|
## 
## 

この学習モデルの正解率は、

0.48 + 0.30

## [1] 0.78

である。

コーエンの \(\kappa\) は、

df_kappa$adab <- pred_adab$class
kappa2(df_kappa[,  c(1, 5)])$value

## [1] 0.5460173

なので、今回の例では最も性能が高いことがわかる。

授業の内容に戻る