\[ % 定義 \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \]
今回の分析で使うパッケージを読み込む。
決定木(decision tree)は特徴量の値に応じて条件分岐を繰り返して木のような構造をとるアルゴリズムである。
以下の図は\(x_1,x_2\)という特徴量を使って目的変数\(y\)の分類を決定木によって行うときのイメージである。 左側の図が木構造のモデルの様子を表しており、まず予測対象のレコードが「\(x_1<a\)」であるかどうかによって分岐し、もしそうであれば\(x_1<a\)がYesであったほうの枝において更に「\(x_2<b\)」であるかどうかで分岐するようにして、そのレコードがどちらのクラスに属するかの分類を行う。 右側の図はこの条件分岐によって生み出される識別境界(\(x_1, x_2\)の軸からなる空間の分割)を示している。
なお、木の分岐の結節点や終端をノード(node)と呼び、木の末端にあるノードは終端ノード(terminal node)や葉ノード(leaf node)と呼ぶ。各葉ノードは、右側の図の分割された各部分空間に対応している。
線形回帰の場合、予測モデルの学習とはパラメータ\(\beta\)を推定することであったが、決定木では最適な分割を推定して木を分岐させていくことが学習にあたる。以下でその仕組みを概説する。
目的変数\(y\)がカテゴリカル変数である場合、決定木の予測値は葉ノード内の訓練データのうち、最も割合の多いクラスを使う(上の図のように)。
目的変数\(y\)が量的変数である場合、葉ノード内の学習データの\(y\)の平均値を予測値とする。
決定木では、他の機械学習アルゴリズムが誤差関数で予測精度の良し悪しを評価していたのと同様に、分割の良し悪しについての指標である不純度(impurity)という関数を用いて予測精度が最もよい分割を探索する。
目的変数\(y\)が\(K\)個のクラスからなるカテゴリカル変数である場合、不純度にはジニ係数(Gini index)
\[ I(t) = \sum_{k=1}^K p_{tk} (1 - p_{tk}), \hspace{1em} (k=1,2, ...,K) \]
が使われることが多い。ここで\(t\)はノード、\(k\)は目的変数のクラスで、\(p_{tk}\)はノード\(t\)の領域に属する訓練データのサンプルサイズ\(N(t)\)のうち目的変数がクラス\(k\)であるデータ数\(N_k(t)\)の割合
\[ p_{kt} = \frac{N_k(t)}{N(t)} \]
である。 ベルヌーイ分布の分散は\(p(1-p)\)であるため、ジニ係数は数式の形としては各クラスごとの出現率の分散を計算して総和をとったような形となっている。
分類すべきクラスが2つ(\(K=2\))である場合のジニ係数は次の図のような形になる。
例えば\(p_{tk}=0.5\)のとき、すなわち2つのクラスに属するデータが半々の比率で混じり合っている場合には、その葉ノードのジニ係数は
\[ \begin{align} I(t) &= 0.5 \times (1 - 0.5) + 0.5 \times (1 - 0.5)\\ &= (0.25 + 0.25) = 0.5 \end{align} \]
となる。グラフからわかるように、これは2クラス分類の状況下での最大の不純度である。
逆に、同クラスに属するデータのみをうまく切り分けることに成功し、1つめのクラスに属するデータのみが葉ノード\(t\)に含まれる場合は\(p(k=1|t)=1.0\)であるから
\[ I(t) = 1 \times (1 - 1) + 0 \times (1 - 0) = 0 \]
となる。
回帰(目的変数が量的変数)の場合、不純度\(I(t)\)にはノード\(t\)内の学習データでの予測値\(\hat{y}(t)\)と実測値\(y_i\)の平均2乗誤差
\[ I(t) = \frac{1}{N(t)} \sum_{i \in t} \big( y_i - \hat{y}(t) \big)^2 \]
を使用する。なお、予測値\(\hat{y}\)がそのノード\(t\)内での平均値であるため、これはノード内での\(y\)の分散である。
不純度を使って、冒頭の例における「\(x_1<a\)」のような分岐の条件を探索していく。
ノード\(t\)から左右に分割される子ノード達を\(t_{L},t_{R}\)と表記する。不純度\(I(t_L),I(t_R)\)と子ノードの領域に属するサンプルサイズ\(N(t_L),N(t_R)\)を用いて算出される、分割後の不純度の重み付き和
\[ \left\{\frac{N(t_L)}{N(t)} I(t_{L}) + \frac{N(t_R)}{N(t)} I(t_{R})\right\} \]
と分割前の不純度\(I(t)\)との差
\[ I(t) - \left\{\frac{N(t_L)}{N(t)} I(t_{L}) + \frac{N(t_R)}{N(t)} I(t_{R})\right\} \]
を情報利得(information gain)と呼び、これを最大化する(分割後の不純度の和を最も下げる)ような分割の条件を探索する。
より具体的には、以下のような手順で学習が行われていく。
{carData}
パッケージに含まれるタイタニック号の乗客データTitanicSurvival
を使う。
## survived sex age passengerClass
## Allen, Miss. Elisabeth Walton yes female 29.0000 1st
## Allison, Master. Hudson Trevor yes male 0.9167 1st
## Allison, Miss. Helen Loraine no female 2.0000 1st
## Allison, Mr. Hudson Joshua Crei no male 30.0000 1st
## Allison, Mrs. Hudson J C (Bessi no female 25.0000 1st
## Anderson, Mr. Harry yes male 48.0000 1st
このデータセットは1912年のタイタニック号の沈没事故の乗客の生死に関するデータで、次の変数が含まれている
survived
:生存したかどうかsex
:性別age
:年齢(1歳に満たない幼児は小数)。263の欠損値を含む。passengerClass
:船室の等級このデータセットには欠損値が含まれているため、まず欠損値を除去する
そしてデータを学習用・テスト用に分割する。
# ID列を追加
df <- titanic %>% rownames_to_column("ID")
# 80%を学習用データにランダムサンプリング
set.seed(1) # 乱数の種を固定
train <- df %>% sample_frac(size = 0.8)
# 学習用データに使っていないIDの行をテスト用データに
test <- anti_join(df,
train,
by = "ID")
# ID列は予測に使わないため削除しておく
train <- train %>% dplyr::select(-ID)
test <- test %>% dplyr::select(-ID)
決定木の実行には{rpart}
パッケージを使う。決定木の作図用に{partykit}
も使う。
rpart()
関数を用いて決定木を実行できる。lm()
関数と同様に、formula
の引数にモデルの構造すなわち被特徴量・特徴量の関係を記述しdata
の引数にデータのオブジェクトを入れる。
# rpartを実施し、treeに結果を保存
# lmと同様に「y ~ x1 + x2」の形で変数を指定できる。"." は「残りのすべての変数」の意味
tree <- rpart::rpart(survived ~ . , data = train)
ここでは表示しないが、結果はsummary()
やprint()
で表示できる。
{partykit}
パッケージを使えば決定木を描くことができる。
まず性別で分岐し、男性の場合は年齢が9.5歳未満だと生存する割合が高くなるが、9.5歳以上だと生存する割合は低くなっている。
女性の場合、船室の等級が1stや2ndであれば生存する割合が極めて高いが、3rdの場合は年齢が高いほど生存する割合は低くなっていることがわかる。
決定木はこのようにデータの解釈・説明を行うことができる。そのため線形回帰などと並んでマーケティングの現場などビジネスにおいても使われる。
予測値はpredict()
で算出できる。type = "class"
とすると目的変数のクラスがfactor型で出力され、type = "prob"
だと各クラスの予測確率値が出力される。
# 予測値を算出
y_pred <- predict(tree, test, type = "class")
# 混同行列: 実際の値と予測値を比べる
table(test$survived, y_pred)
## y_pred
## no yes
## no 120 3
## yes 41 45
## [1] 0.7894737
回帰問題の場合で、特徴量\(x\)を横軸に、目的変数\(y\)を縦軸にとった場合の散布図と、決定木と線形回帰の予測値を線としてプロットしたものが次の図である。
このような非線形データに対しては線形回帰よりも決定木のほうが当てはまりがよい。
これまで扱ってきた機械学習アルゴリズム(線形回帰、ロジスティック回帰、リッジ回帰など)は線形モデルと呼ばれる。 線形モデルは回帰問題においては特徴量と目的変数が線形の関係にあるようなデータを仮定し、それ以外のデータでは予測誤差が大きくなる。
次の図は特徴量\(x1, x2\)をそれぞれ横軸と縦軸にとり目的変数\(y\in\{0, 1\}\)で色を塗り分けた散布図に、線形モデル(ロジスティック回帰)と決定木の予測値がどちらのクラスに決定されるかの境目(決定境界 decision boundary)の線を書き込んだものである。
2次元平面において2つの集合を1つの線で分断できることを線形分離可能(linearly separable)という。 線形モデルは線形の決定境界を生むため、線形分離可能なデータの分類問題とは相性が良い。
しかし、線形モデルは線形分離不可能問題ではうまく分類することができない。 決定木は非線形モデルであるため、線形分離不可能なデータでもうまく扱うことができる。
決定木は異なる学習データで学習し直した時に学習結果が大きく変化しやすい特性があり、単一の決定木での予測精度はあまり高くない。
しかし、複数の予測モデルを統合して予測値を出すアンサンブル学習(ensemble learning)という学習方法を用いることで、予測モデルの予測精度を向上させることができることが知られている。
例えば、60%の確率で正しい分類ができる予測モデルがあるとする。もしこの予測モデルを単体で使うのであればあまり信用できないモデルである。しかし、もしそのようなモデルが10個あり互いの予測の相関がなければ、平均的には10個のうち6個のモデルは正しい分類を行うと考えられるため、それらで多数決を行えば良い精度の予測モデルが得られると考えられる。
以下の図はこの考えをもとにシミュレーションを行ったものである。「単体」の予測モデルは試行回数を増やすにつれて0.6に収束していくが、「10個の多数決」のほうは約0.85の累積正解率に収束しているのが分かる。
ランダムフォレスト(random forest)は上記の例のような考え方に基づき「互いの相関が低い複数の決定木を生成して多数決をとる」というアルゴリズムである(Breiman, 2001)。 ランダムフォレストの学習は以下のような流れで行われる。
なお、このようにブートストラップによって多数の予測モデルを構築した後にそれらの予測結果を集約(aggregating)するアンサンブル学習の方法はバギング(bagging; bootstrap aggregatingの略)と呼ばれる。
{randomForest}
パッケージを使う。
set.seed(0)
# randomForestを実施し、titanic_RFに保存。
titanic_RF <- randomForest(survived ~ .,
data = train,
mtry = 3, # ランダムに選んで使用する特徴量の数
nodesize = 5, # 葉ノードに含まれる学習データの最小数
sampsize = 100, # ブートストラップサンプルのサンプルサイズ
ntree = 2000) # ブートストラップの反復回数であり、作成する決定木の数
結果の概略はprint
によって表示できる。
##
## Call:
## randomForest(formula = survived ~ ., data = train, mtry = 3, nodesize = 5, sampsize = 100, ntree = 2000)
## Type of random forest: classification
## Number of trees: 2000
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 20.31%
## Confusion matrix:
## no yes class.error
## no 448 48 0.09677419
## yes 122 219 0.35777126
なお、ここでのOOB estimate of error rate
とはランダムフォレストの予測誤差の推定値であり、ブートストラップサンプリングを行った際に使われなかった学習データが生じることを利用して予測精度を測っている(OOBはout-of-bagの意味)。
titanic_RF$err.rate
には木の数ごとのOOB誤り率が含まれているため、ntree
に指定すべき最適な木の数を考える際の参考になる。
テストデータを用いて予測精度を測る。
## y_pred
## no yes
## no 115 8
## yes 32 54
## [1] 0.8086124
また、ランダムフォレストでは、$importance
でジニ係数の減少に基づいて算出される特徴量重要度(feature
importance)を確認することができる。これは分割時の不純度の減少量の観点から、予測モデルにおいて各特徴量がどの程度予測に寄与しているのかの目安を示している。
## MeanDecreaseGini
## sex 13.764060
## age 15.434422
## passengerClass 6.568634
ランダムフォレストとはまた別のアンサンブル学習の方法の一つにブースティング(boosting)がある。 ブースティングでは複数の予測モデルを用いて逐次的に学習を行い、\(m\)回目の学習ではその1つ前の\(m-1\)回目に学習した予測モデルで予測誤差が多かったサンプルを重視する、というようなアプローチで学習していく。
複数の関数に重みをかけて足し合わせた関数
\[ f(x) = f_0(x; \theta_0) + \beta_1 f_1(x; \theta_1) + \cdots + \beta_M f_M(x; \theta_M) \]
の形で予測モデルを構築することを考える。ここで\(\theta_0, \dots, \theta_M\)は関数を形づくるパラメータ(例えば線形回帰の重みや決定木の分岐の閾値)である。
このモデルは\(\beta_1, \beta_2, \dots, \beta_M\)と\(\theta_0, \theta_1, \dots, \theta_M\)のパラメータを推定する必要がある。 今回はすべてのパラメータを一度に学習するのではなく、\(\beta_m f_m(x; \theta_m)\)を一つずつ学習していく方法を考える。具体的には次のように行う。
前向き段階的加法モデリング(forward stagewise additive modeling)
これを前向き段階的加法モデリング(forward stagewise additive modeling)という。ブースティングはこの方法でアンサンブル学習を行う。
誤差関数が二乗誤差\(L(y, f(x)) = (y - f(x))^2\)の場合、
\[ \begin{align} L(y_i, f_{m-1}(x_i) + \beta f(x_i; \theta)) &= (y_i - f_{m-1}(x_i) - \beta f(x_i; \theta))^2 \\ &= (\text{residual}_{i,m-1} - \beta f(x_i; \theta))^2 \end{align} \]
となり、\(m-1\)回目のモデルの残差\(\text{residual}_{i,m-1} = y_i - f_{m-1}(x_i)\)を近似するように\(m\)回目のモデル\(\beta f(x_i; \theta)\)を学習させていると捉えることができる。残差が大きければそれだけ訓練中に重視されるため「間違えた箇所を重点的に学習する手法」とも捉えることができる。
数理最適化において関数の最小化問題
\[ \min_x f(x) \]
を解く方法のひとつに勾配降下法(gradient descent method)あるいは最急降下法(steepest descent method)と呼ばれるものがある。これは目的関数の1次の導関数を用いて
\[ x_{m} = x_{m-1} - \alpha_{m-1} \frac{\partial f(x_{m-1})}{\partial x_{m-1}} \]
という値の更新を何度も繰り返して最適化を行っていく。ここで\(\alpha\)は学習率と呼ばれるパラメータで、値の更新量が多すぎると最適解を通り過ぎてしまうことがあるので小さめの値を乗じて更新幅を抑えるために用いられる。
(ちなみに名前の由来である勾配(gradient)とは偏微分のベクトル\(\nabla f(\boldsymbol{w}) = \left( \frac{\partial f(\boldsymbol{w})}{\partial w_1}, \dots, \frac{\partial f(\boldsymbol{w})}{\partial w_n} \right)\)のことである)
最終的に\(M\)回反復して得た最適解\(x^*\)は
\[ x^* = x_0 - \alpha_1 \frac{\partial f(x_{1})}{\partial x_{1}} - \alpha_2 \frac{\partial f(x_{2})}{\partial x_{2}} - \cdots - \alpha_M \frac{\partial f(x_{M})}{\partial x_{M}} \]
となり、ブースティングにより得られる予測モデル
\[ f(x) = f_0(x; \theta_0) + \beta_1 f_1(x; \theta_1) + \cdots + \beta_M f_M(x; \theta_M) \]
と同様に重み付き和の形になる。
ブースティングは勾配降下法を機械学習で行っていると捉えることができる。
機械学習においては予測値\(f(x)\)と実測値\(y\)の誤差の最小化問題
\[ \min_{f(x)} L(y, f(x)) \]
を解きたいため、勾配は誤差関数の予測モデルによる微分\(\frac{ \partial L(y, f(x)) }{ \partial f(x) }\)によって得られる。
二乗誤差\(L(y, f(x)) = \frac{1}{2} (y - f(x))^2\)の場合、負の勾配は残差である
\[ - \frac{ \partial L(y, f(x)) }{ \partial f(x) } = y - f(x) = \text{residual} \]
前向き段階的加法モデルの節で「二乗誤差の場合は残差を近似するように学習している」と述べた。
\[ L(y_i, f_{m-1}(x_i) + \beta f(x_i; \theta)) = (\text{residual}_{i,m-1} - \beta f(x_i; \theta))^2 \]
これにより、学習されるモデル\(\beta f(x; \theta)\)は負の勾配を学習するようになり、最終的にそれらの和となるモデルは勾配降下法を解いた状態を近似することになる。
勾配ブースティング決定木(gradient boosting decision tree: GBDT)は上記のようなブースティングによる加法モデルの構築を決定木\(T(x; \theta)\)を用いて行うアルゴリズムである。
\[ f(x) = T_0(x; \theta_0) + \beta_1 T_1(x; \theta_1) + \cdots + \beta_M T_M(x; \theta_M) \]
勾配ブースティング決定木は表形式データ(tabular data)に対して高い性能を発揮することが多く、頻繁に利用される。 有名な勾配ブースティング決定木のアルゴリズムにXGBoostとLightGBMがあり、Rからも利用可能なのでそれらを紹介していく。
XGBoost(Chen & Guestrin, 2016)は勾配ブースティング決定木を拡張したアルゴリズムおよびパッケージである。 内部のアルゴリズムは上記のような古典的な勾配ブースティングからかなり改善されているが、基本的な考え方は同じなので本講義では割愛し、Rでの使い方の解説を重視していく。
XGBoostはfactor型を使うことができないため、データ型を数値のみに変更する。factor型をas.numeric()
で数値型にするとカテゴリに応じた整数に変換されるので、それを用いる。
(なお、決定木ベースのアルゴリズムでは回帰分析のような線形モデルとは異なり、カテゴリカル変数をダミー変数にせずにカテゴリに応じた整数にするだけでも扱うことができる。決定木の内部で3 < x & x <= 4
のような分岐を作れば「x
= カテゴリ4」を表現できるためである。)
# factor型の変数たちを数値に変えていく
y_train <- 1 * (train$survived == "yes")
X_train <- train %>% mutate(
sex = as.numeric(sex),
passengerClass = as.numeric(passengerClass)
) %>% select(-survived)
xgboost()
関数を呼び出すことで学習を行うことができる。
# xgboostの学習
xgb <- xgboost(
data = as.matrix(X_train), # 特徴量(matrix型にする必要があるので注意)
label = y_train, # 目的変数(教師ラベル)
objective = "binary:logistic", # 目的関数(誤差関数) binary = 二値分類
max.depth = 4, # 決定木の深さ。大きすぎると過学習を招く
nrounds = 100, # 反復回数
learning_rate = 0.1 # 学習率
)
XGBoostには様々なハイパーパラメータが存在し、モデルの挙動を細かく設定していくことができる。
objective
: 目的関数(誤差関数)
"reg:squarederror"
と指定するnrounds
: 反復回数。使用する決定木の数でもある。learning_rate
:
学習率。追加される決定木にどれだけの重みを乗じるか。
リッジ回帰やLASSOのように正則化の設定を行うこともできる。
lambda
: L2正則化alpha
: L1正則化ランダムフォレストのようにブートストラップサンプリングを行って決定木を構築することもできる。これは次の2つのハイパーパラメータを1より小さい値にすればよい。
subsample
: 訓練データのサンプリングの割合colsample_bytree
:
決定木の構築に使う特徴量のサンプリングの割合# testもfactor型を数値に変える
y_test <- 1 * (test$survived == "yes")
X_test <- test %>% mutate(
sex = as.numeric(sex),
passengerClass = as.numeric(passengerClass)
) %>% select(-survived)
# 予測値を算出する
p_pred <- predict(xgb, as.matrix(X_test)) # predictの返り値はy=1の確率値P(y=1)なので注意
y_pred <- 1 * (p_pred > mean(p_pred)) # 平均以上なら1(生存)とすると変換
# 混同行列
table(y_test, y_pred)
## y_pred
## y_test 0 1
## 0 107 16
## 1 26 60
## [1] 0.7990431
LightGBM (Ke et al. 2017) は計算が非常に高速で大規模データとの相性が良い勾配ブースティング決定木のアルゴリズムである。 XGBoostなど従来の勾配ブースティング決定木に比べて、
などの特徴を持っている1。
通常の決定木は分岐点を探索する時間計算量(計算回数)が多い。例えばある特徴量にユニークな値が\(m\)個あった場合、分岐点の候補はそれらの間となるため\(m-1\)個の候補点があり、それぞれの点において分岐するかどうかを判定していく必要がある。
LightGBMはHistogram-based Treeを採用している。これはヒストグラムのように連続変数をいくつかのグループ(bin)に分割して扱っていくタイプの決定木である。最大のbinの数もハイパーパラメータの一つであり、デフォルトでは256であるため、その場合は分岐の候補点は255個に絞られる。
この処置は\(m-1\)回の探索を厳密に行う方法に比べると大幅な近似となってしまうため解の厳密性を落とすものの、サンプル数が非常に多くなり\(m\)が数百万となっているようなデータセットにおいてはこの工夫による計算量の削減が著しいものとなる。
また近似しているといってもLightGBMがXGBoostに比べて性能が劣るとは限らず、経験的には多くのタスクにおいてXGBoostと同程度の性能を発揮することが知られている。
通常の決定木はx <= 10
のような大小比較で分岐を行うが、LightGBMではx == 10
のような一致の判定や複数カテゴリの集合に含まれるかどうか(x in (1, 3, 6)
)の判定を行うことができる。
これによりカテゴリカル変数をより効率よく扱うことができるようになっている。
d_train <- lgb.Dataset(
data = as.matrix(X_train),
label = y_train,
categorical_feature = c(1L, 3L) # カテゴリカル変数のインデックス番号を指定(1Lは1番目の意味。1Lはinteger型になる)
)
lgb <- lgb.train(
params = list(
objective = "binary", # 目的関数
max.depth = 4, # 決定木の深さ。大きすぎると過学習を招く
learning_rate = 0.1 # 学習率
),
data = d_train,
nrounds = 10, # 反復回数
)
# 予測値を算出する
p_pred <- predict(lgb, as.matrix(X_test)) # predictの返り値はy=1の確率値P(y=1)なので注意
y_pred <- 1 * (p_pred > mean(p_pred)) # 平均以上なら1(生存)とすると変換
# 混同行列
table(y_test, y_pred)
## y_pred
## y_test 0 1
## 0 106 17
## 1 26 60
## [1] 0.7942584
Random Forests
XGBoost
LightGBM
書籍
ウェブサイト & スライド
厳密にはLightGBMの論文で新規に提案された技術はGOSS(勾配情報に基づく効率的なサブサンプリング)とEFB(ダミー変数のようにゼロが多くスパースなデータを効率良く扱う技術)である。しかしLightGBMのパッケージにおいてGOSSはデフォルトでは無効化されており、EFBはスパースではない通常のデータではあまり恩恵がないため、実用上のLightGBMの目立った特徴は列挙した2点であると判断した。↩︎