Rで決定木分析(rpartによるCARTとrangerによるランダムフォレスト)


準備

決定木(decision tree)分析をする際、まず目的変数の種類とアルゴリズムを決定する。

アルゴリズム

  • CART
  • CHAID
  • ID3 / C4.5 / C5.0

目的変数の型

目的変数の型によって扱いが変わる

  • 質的変数(2値変数):分類木→目的変数が0/1, T/Fの場合はas.factor()でfactor型にデータ変換しておく
  • 量的変数:回帰木
  • survivalオブジェクト
  • (生起を表す2カラム)

CARTはすべて対応、C4.5/C5.0は質的変数のみ

ここではCARTアルゴリズムでツリーモデルを生成するrpartと、ランダムフォレストrangerを中心に説明する。

1回限りの決定木

実行

require(rpart)
purchase.dt %>% 
  filter(!is.na(age_range)) %>%
  mutate(
    age_range = ifelse(age_range == '(0,30]', '30>=', '30<'),
    purchase = as.factor(purchase > 0)
  ) %>%
  dplyr::select(-ends_with('_id'), -ends_with('_date')) %>%
  rpart(
    formula = purchase ~ .,
    data = .,
    method = 'class',
    parms = list(split='information'),
    control = rpart.control(minsplit=10, cp=0.001)
  ) -> purchase.rpart
  • 投入するデータテーブルの変数を絞っておくと変数指定が楽
  • プロットしたときのラベルをわかりやすく変換しておくといい
  • 分類木の場合は目的変数をfactor型に変換しておく
  • 重要なパラメータ
    • method(通常は目的変数の型によって自動で最適なものが選択される)
      • ‘class’で分類木(目的変数がfactor型)
      • ‘poisson’で生起(目的変数が2カラムの生起データ)
      • ‘exp’で生存(目的変数がsurvicalオブジェクト)
      • ‘anova’で回帰木(目的変数が上記のいずれでもない)
    • parms: method = 'class'の場合、以下の指標に基づいて分割。method = 'anova'の場合は指定しない
      • parms = list(split='gini')でジニ係数を使う(デフォルト)
      • parms = list(split='information')でエントロピーを使う
    • rpart.control
      • minsplitは1ノードのサイズの下限
      • cpは小さいほど細かく分岐する。あとで粗くできるので最初は細かく分けておくといい

見る

summary(purchase.rpart)

変数の重要度を確認

purchase.rpart$variable.importance

チューニング(cpを調整)

cpはツリーモデルの複雑さを表すパラメータ。値が小さいものほどモデルが細かくなる。

cpを見る

printcp(purchase.rpart)
plotcp(purchase.rpart)

cpを調整(cpの小さいツリーモデルからcpの大きいツリーモデルへ)

purchase.rpart <- prune(purchase.rpart, cp=0.01)

プロット(2通りのライブラリで)

ビルトインのplot

par(xpd = TRUE)
plot(purchase.rpart, compress = TRUE)
text(purchase.rpart, use.n = TRUE)

(決定木のプロット例)

ビルトインのplot

ライブラリ{rpart.plot}

require(rpart.plot)
rpart.plot(purchase.rpart)

(決定木のプロット例)

rpart.plot

ライブラリ{partykit}

require(partykit)
plot(as.party(purchase.rpart))
plot(as.party(purchase.rpart), gp = gpar(fontsize = 7))

(決定木のプロット例)

partyをplot

どれを使ってもいいが、使いやすいものを確保しておくといい。

分類されたノードを元データに紐づける

purchase.dt[!is.na(age_range), node := purchase.rpart$where]

予測

predict(purchase.rpart, test.dt)

よくある処理をまとめて

# 学習用データの抽出
purchase.dt[,customer_id] %>% unique %>% sample(size = 3000) -> customer_id.train
purchase.dt[customer_id %in% customer_id.train,] -> purchase.dt_train
purchase.dt[!customer_id %in% customer_id.train,] -> purchase.dt_test
# 決定木の実行
purchase.dt_train %>% 
  filter( ... ) %>%
    : 
  rpart( ... ) -> purchase.rpart
# 確認
printcp(purchase.rpart)
purchase.rpart %>% prune(cp = 0.05) %>% summary %>% as.party %>% plot(gp = gpar(fontsize = 7))
# ノードの紐づけ
purchase.dt_train[, node := purchase.rpart$where]
# 検証
predict(purchase.rpart, purchase.dt_test)

ランダムフォレスト(random forest)

rangerパッケージが便利。
以下の目的変数ごとのツリーモデルをサポートしている。

  • 質的変数(2値変数):分類木→目的変数が0/1, T/Fの場合はas.factor()でfactor型に変換しておく
  • 量的変数:回帰木
  • survivalオブジェクト

実行

require(ranger)
purchase.dt %>% 
  filter(!is.na(age_range)) %>%
  mutate(
    age_range = ifelse(age_range == '(0,30]', '30>=', '30<'),
    purchase = as.factor(purchase > 0)
  ) %>%
  dplyr::select(-ends_with('_id'), -ends_with('_date')) %>%
  ranger(
    formula = purchase ~ .,
    data = .,
    num.trees = 1000,
    mtry = 5,
    write.forest = TRUE,
    importance = 'impurity',
    na.action = na.omit,
    probability = FALSE
  ) -> purchase.ranger

パラメータ

ランダムフォレストのハイパーパラメータは

  • num.trees: 試す決定木の数
  • mtry: モデルに採用する変数の数

mtryをグリッドサーチするならcaretでmethod='ranger'を指定する(後述)。

その他の主なパラメータ

  • min.node.sizeでノードサイズの下限を指定できる
  • importanceを与えると変数の重要度を返す。回帰木と分類木では’impurity’を(ジニ係数)、生存では’permutation’を指定する(デフォルトでは重要度を返さない)。
  • 目的変数が質的変数の時、probability = TRUEで確率を返す。デフォルトのFALSEではT/Fの応答(ただしfactor)を返す。確率は2列の行列で、logical型の目的変数をfactor型にしてranger()をかけている場合、2列目がTRUEとなる(FALSE, TRUEの順)

なおデフォルトでは

  • 分類木ではGini係数に基づいて、
  • 回帰木では分散に基づいて、
  • 生存モデルではログランクに基づいて

分割する。

戻り値

戻り値のrangerオブジェクトはリストで、よく使う属性が

  • predictionsが予測結果
  • variable.importanceが変数の重要度

なお結果のprediction.error

  • 分類木では誤分類の割合
  • 回帰木では平均二乗誤差(MSE)
  • 生存モデルではc-index

が使われる。

モデルの評価

partial dependence plot

partial dependence plot(部分従属プロット)を描くにはedarfパッケージを使う。ranger以外にもrandomForest, RandomForest, rfsrcのランダムフォレストオブジェクトに対応している。

require(edarf)
pd <- partial_dependence(customer.ranger, vars = c('n_purchase','age'), data = as.data.frame(customer_train.dt))
plot_pd(pd)

partial_dependence()の引数dataはdata.tableではダメで、data.frameでなければならない

このパッケージedarfはランダムフォレストの診断に便利。
https://qiita.com/nakamichi/items/bed7a2f180ea9ce86d94

精度

分類木の場合

# 予測
## probability = TRUEでranger()を実行している場合、2列目がTRUE
purchase.predicted <- predict(customer.ranger, customer_test.dt)$predictions[,2]
## probability = FALSEでranger()を実行している場合
purchase.predicted <- predict(customer.ranger, customer_test.dt)$predictions
# 正解データ
purchase.label <- customer_test.dt$response # logical型

confusion matrix(混同行列)

table(purchase.predicted, purchase.label)

ROC曲線
ROC曲線を描くにはROCRパッケージを使う。

require(ROCR)
pred <- prediction(purchase.predicted, purchase.label)
performance(pred, "tpr", "fpr") %>% plot

AUC

unlist(performance(pred,'auc')@y.values)

ハイパーパラメータをグリッドサーチ

purchase.dt %>% 
  filter(!is.na(age_range)) %>%
  mutate(
    age_range = ifelse(age_range == '(0,30]', '30>=', '30<'),
    purchase = as.factor(purchase > 0)
  ) %>%
  dplyr::select(-ends_with('_id'), -ends_with('_date')) %>%
  train(
    formula = purchase ~ .,
    data = .,
    method = 'ranger',
    metric = 'Accuracy',
    write.forest = TRUE, importance = 'impurity',
    num.trees = 100,
    na.action = na.omit,
    tuneGrid = expand.grid(mtry = 3:7),
    trControl = trainControl(method = 'cv', number = 5, allowParallel = TRUE)
  ) -> purchase.ranger.train
  • na.actionはレコードにNAを含む場合の扱いを指定する。デフォルトはna.failで処理に失敗する。NAを含む行を除外して実行する場合はna.omitを指定する。

分類木では(ロジスティック回帰のようなその他のクラス分類においても)metricに’Accuracy’, ‘Kappa’, ‘ROC’を指定できる。ROCを使う場合、trControlの引数にclassProbs = TRUEを加える必要がある。

参考 – rpartに対してtrainすると

purchase.dt %>% 
  filter(!is.na(age_range)) %>%
  mutate(
    age_range = ifelse(age_range == '(0,30]', '30>=', '30<'),
    purchase = as.factor(purchase > 0)
  ) %>%
  dplyr::select(-ends_with('_id'), -ends_with('_date')) %>%
  train(
    formula = purchase ~ .,
    data = .,
    na.action =  na.omit, # or na.fail
    method = 'rpart',
    parms = list(split='information'),
    control = rpart.control(minsplit = 10, cp = 0.001),
    metric = 'Accuracy',
    trControl = trainControl(method = 'repeatedcv', number = 10, repeats = 10, allowParallel = TRUE)
#  metric = 'ROC',
#  trControl = trainControl(method = 'repeatedcv', number = 10, repeats = 10, allowParallel = TRUE, classProbs = TRUE),
) -> a.rpart.train

cpの値を探索することになる。

R関連記事