Dropout

データサイエンスについて勉強したことを書いていきます。機械学習、解釈性、因果推論など。

XGBoostの論文を読んだのでGBDTについてまとめた

はじめに

今更ですが、XGboostの論文を読んだので、2章GBDT部分のまとめ記事を書こうと思います。*1
この記事を書くにあたって、できるだけ数式の解釈を書くように心がけました。数式の意味をひとつひとつ追っていくことは、実際にXGBoost(またはLightGBMやCatBoostなどのGBDT実装)を使う際にも役立つと考えています。たとえばハイパーパラメータがどこに効いているかを理解することでチューニングを効率化したり、モデルを理解することでよりモデルに合った特徴量のエンジニアリングができるのではないかと思います。

また、この記事に限りませんが、記述に間違いや不十分な点などあればご指摘頂ければ嬉しいです。

XGBoost論文

目的関数の設定


一般的な状況として、サンプルサイズが Iで特徴量の数が Mのデータ \mathcal{D} = \left\{ (\mathbf{x}_i, y_i) \right\}(i \in\mathcal{I} = \{1, \dots, I\}, \; \mathbf{x}_i \in \mathbb{R}^M, \; y_i \in\mathbb{R})に対する予測モデルを構築することを想定しましょう。*2
今回はツリーをアンサンブルした予測モデルを構築します。
 \mathcal{K}  =\{1, \dots, K\}のツリーを加法的に組み合わせた予測モデルは以下のように定式化できます。*3


\begin{align}
\hat{y}_{i} &= \phi\left(\mathbf{x}_{i}\right)=\sum_{k\in\mathcal{K}} f_{k}\left(\mathbf{x}_{i}\right),\\
\text{where}\quad f_{k} \in \mathcal{F} &= \left\{f(\mathbf{x})=w_{q(\mathbf{x})}\right\}\left(q : \mathbb{R}^{m} \rightarrow \mathcal{T},\; \mathcal{T} = \{1, \dots, T\}, \; w \in \mathbb{R}^{T}\right)
\end{align}

ここで、 f_kはひとつひとつのツリーを表しています。ツリー f(\mathbf{x})は特徴量 \mathbf{x} が与えられると、それを q(\mathbf{x})に従って各ノード t = 1, \dots, Tに紐づけ、それぞれのノードに対応する予測値 w_{q(\mathbf{x})}を返します。そして、ひとつひとつのツリーの予測値を足し合わせることで、最終的な予測結果 \hat{y}_iとします。



では、具体的にツリーをどうやって作っていくかを決めるために、最小化したい目的関数 \mathcal{L}(\phi)を設定します。

\begin{align}
\mathcal{L}(\phi) &= \sum_{i\in \mathcal{I}} l\left(y_{i}, \hat{y}_{i}\right)+\sum_{k\in \mathcal{K}} \Omega\left(f_{k}\right), \\
\text{where}\quad\Omega(f) &= \gamma T+\frac{1}{2} \lambda\|w\|^{2} = \gamma T+\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2}
\end{align}

ここで、 l(y_{i}, \hat{y}_{i})は損失関数で、たとえば二乗誤差になります。ただし、単に二乗誤差を最小化するのではなく、過適合を回避して汎化性能を上げるために正則化 \Omega(f)が追加されています。なお、 \gamma \lambda はハイパーパラメータであり、交差検証などで最適な値を探索する必要があります。

  •  \Omega(f)の第一項 \gamma Tはツリーのノードの数に応じてペナルティが課されるようになっています。ハイパーパラメータ \gammaを大きくするとよりノード数少ないツリーが好まれるようになります。ツリーの大きさに制限をかけることで過適合を回避することが目的です。
  •  \Omega(f)の第二項 \frac{1}{2}\lambda\|w\|^{2}は各ノードが返す値の大きさに対してペナルティがかかることを意味しています。ハイパーパラメータ \lambdaを大きくすると、(絶対値で見て)より小さい wが好まれるようになります。 wが小さいということは最終的な出力を決める \sum_{k\in\mathcal{K}} f_{k}部分で足し合わされる値が小さくなるので、過適合を避けることに繋がります。

勾配ブースティング

さて、目的関数 \mathcal{L}(\phi)を最小化するような K個のツリー構築したいわけですが、 K個ツリーを同時に構築して最適化するのではなく、 k個目のツリーを作る際には、 k-1個目までに構築したツリーを所与として、目的関数を最小化するようなツリーを作ることにしましょう。*4


\begin{align}
\min_{f_k}\;\mathcal{L}^{(k)}=\sum_{i\in\mathcal{I}} l\left(y_{i}, \hat{y}_i^{(k - 1)}+f_{k}\left(\mathbf{x}_{i}\right)\right)+\Omega\left(f_{k}\right)
\end{align}

このステップで作成する k個目のツリーを合わせた予測値は \hat{y}_i^{(k)} = \hat{y}_{i}^{(k - 1)}+f_{k}\left(\mathbf{x}_{i}\right)であり、 k-1個目までのツリーではうまく予測できていない部分に対してフィットするように新しいツリーを構築すると解釈できます。このように残差に対してフィットするツリーを逐次的に作成していく手法をブースティングと呼びます。



さて、損失関数 \sum_{i\in\mathcal{I}} l\left(y_{i}, \hat{y}_{i}^{(k - 1)}+f_{k}\left(\mathbf{x}_{i}\right)\right)を直接最適化するのではなく、その2階近似を最適化することにしましょう。後にわかるように、2次近似によって解析的に解を求めることができます。 f_k = 0の周りで2階のテイラー展開を行うと、目的関数 \mathcal{L}^{(k)}は以下で近似できます。

\begin{align}
\mathcal{L}^{(k)} &\approx \sum_{i\in\mathcal{I}}\left[l\left(y_{i}, \hat{y}^{(k - 1)}\right)+g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right] + \Omega\left(f_{k}\right),\\
\text{where} \quad g_i &= \frac{\partial }{\partial \hat{y}^{(k - 1)}}l\left(y_{i}, \hat{y}^{(k - 1)}\right),\\
h_i &= \frac{\partial^2 }{\partial \left(\hat{y}^{(k - 1)}\right)^2}l\left(y_{i}, \hat{y}^{(k - 1)}\right)
\end{align}

ここで、 g_i h_iはそれぞれ損失関数の1階と2階の勾配情報になります。勾配情報を使ったブースティングなので勾配ブースティングと呼ばれています。*5
今回 f_kを動かすことで目的関数を最小化するので、 f_kと関係ない第一項は無視できます。

\begin{align}
\tilde{\mathcal{L}}^{(k)} &=\sum_{i\in\mathcal{I}}\left[g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right]+\gamma T+\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2} \\
&= \sum_{t \in \mathcal{T}}\left[\sum_{i \in \mathcal{I}_{t}} g_{i} f_{k}\left(\mathbf{x}_{i}\right)+\frac{1}{2} \sum_{i \in \mathcal{I}_{t}} h_{i} f_{k}^{2}\left(\mathbf{x}_{i}\right)\right]+\gamma T +\frac{1}{2} \lambda \sum_{t \in \mathcal{T}} w_{t}^{2}\\
&=\sum_{t \in \mathcal{T}}\left[\left(\sum_{i \in \mathcal{I}_{t}} g_{i}\right) w_{t}+\frac{1}{2}\left(\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}\right) w_{t}^{2}\right]+\gamma T
\end{align}

  • 1行目の式では、 f_kと関係ない第一項を取り除き、 \Omega(f_k)の中身を書き下しました。
  • 2行目への変換ですが、全ての i\in\mathcal{I}について足し合わせている部分を、まずノード tに所属する部分 i \in \mathcal{I}_t (\mathcal{I}_t = \{i | q(\mathbf{x}_i) = t \}) を足し合わせてから、全てのノード t \in \mathcal{T}について足し合わせるように分解しています。
  • 3行目への変換では、同じノードに所属する f_k(\mathbf{x}_i)は全て w_tを返すというツリーの性質を利用しています。また、 w^2_tの共通部分をくくっています。

さて、 \tilde{\mathcal{L}}^{(k)} w_tに関しての2次式なので、解析的に解くことができます。

\begin{align}
w_{t}^{*}=-\frac{\sum_{i \in \mathcal{I}_{t}} g_{i}}{\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}}
\end{align}

以上で、 k個目のツリーに関して、各ノードが返すべき値 w_t^*が解析的に求まりました。*6
この式からもハイパーパラメータ \lambdaを大きくすると w^*_tが(絶対値で見て)小さくなることが見て取れます。この w^*_tを元の目的関数に代入してあげることで

\begin{align}
\tilde{\mathcal{L}}^{(k)}(q)=-\frac{1}{2} \sum_{t\in\mathcal{T}} \frac{\left(\sum_{i \in \mathcal{I}_{t}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_{t}} h_{i}}+\gamma T
\end{align}

を得ます。あとはツリーの構造 q、言い換えれば特徴量の分割ルールを決める必要があります。たとえば、一番シンプルなケースとして、全く分割を行わない場合( \mathcal{I})と一度だけ分割を行う場合( \mathcal{I}_L, \mathcal{I}_Rに分割)を比較しましょう。分割による目的関数の値の減少分は

\begin{align}
\mathcal{L}_{\text{split}}&= -\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}} h_{i}}+\gamma - \left(-\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}_L} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_L} h_{i}}-\frac{1}{2} \frac{\left(\sum_{i \in \mathcal{I}_R} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_R} h_{i}}+2\gamma \right)\\
&= \frac{1}{2}\left(\frac{\left(\sum_{i \in \mathcal{I}_L} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_L} h_{i}}+\frac{\left(\sum_{i \in \mathcal{I}_R} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}_R} h_{i}} - \frac{\left(\sum_{i \in \mathcal{I}} g_{i}\right)^2}{\lambda + \sum_{i \in \mathcal{I}} h_{i}}\right) - \gamma
\end{align}


であり、これがプラスなら分割を行い、マイナスなら分割を行わないということになります。上式からも、ハイパーパラメータ \gamma を大きくするとより分割が行われなくなることが見て取れます。

ところで、そもそもどの特徴量のどの値で分割するべきかなのでしょうか?一番ナイーブな考え方は、全ての変数に対して全ての分割点を考慮して、一番目的関数の値を減少させるような分割を選ぶというものがあります。ただし、この方法は膨大な計算量が必要になるため、XGBoostでは近似手法が提供されており、3章に記述されています。さらに、4章移行では並列計算や比較実験などが記されています。


まとめ

XGboostの論文を読んだので、自身の理解を深めるために2章GBDT部分のまとめ記事を書きました。
今までなんとなく使っていたハイパーパラメータが具体的にどの部分に効いているのか学ぶことができて、とても有意義だったと感じています。

*1:なお、元論文のノーテーションがおかしい/統一的でないように感じたので、一部表記を変更しています。

*2:元論文では |\mathcal{D}| = Nですが、 i Nが混じるとややこしい気もするので Iにしました。

*3:元論文の q : \mathbb{R}^{m} \rightarrow Tとなっているのですが、 qはインデックス 1, \dots, Tを返す関数なので、タイポかと思われます。

*4:元論文の添字 tがツリーの数 Tとややこしいので kのまま進めることにしました。

*5:たぶん

*6:論文にはこれがimpuityみたいなものと書かれているのですが、不勉強で理解できませんでした。

tidymodelsによるtidyな機械学習(その2:Cross Varidation)

はじめに

本記事ではtidymodelsを用いたCross Validationとハイパーパラメータのチューニングについて紹介したいと思います。 なお、tidymodelsの基本的な操作方法については以下の記事をご覧下さい。

dropout009.hatenablog.com

前処理

まずは前回の記事と同様、訓練/テストデータの分割と前処理を行います。なお、例によってデータはdiamondsを用います。

# パッケージ
library(tidyverse)
library(tidymodels)
set.seed(42)

# 分割
df_split = initial_split(diamonds,  p = 0.8)

df_train = training(df_split)
df_test  = testing(df_split)

# 前処理レシピ
rec = recipe(price ~ ., data = df_train) %>% 
  step_log(price) %>% 
  step_ordinalscore(all_nominal())

Cross Validation

さて、ここからが新しい内容になります。

df_cv = vfold_cv(df_train, v = 10)

> df_cv
#  10-fold cross-validation 
# A tibble: 10 x 2
   splits               id    
   <list>               <chr> 
 1 <split [38.8K/4.3K]> Fold01
 2 <split [38.8K/4.3K]> Fold02
 3 <split [38.8K/4.3K]> Fold03
 4 <split [38.8K/4.3K]> Fold04
 5 <split [38.8K/4.3K]> Fold05
 6 <split [38.8K/4.3K]> Fold06
 7 <split [38.8K/4.3K]> Fold07
 8 <split [38.8K/4.3K]> Fold08
 9 <split [38.8K/4.3K]> Fold09
10 <split [38.8K/4.3K]> Fold10

tidymodels配下のrsampleパケージにはCross Validationを行うための関数vfold_cv()があるので、これを使います。

> df_cv$splits[[1]]
<38837/4316/43153>

splitsを確認します。全43153サンプルを38837の訓練データと4316の検証データに分割していることを示しています。

> df_cv$splits[[1]] %>% analysis()
# A tibble: 38,837 x 10
   carat cut       color clarity depth table price     x     y     z
   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
 1 0.23  Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
 2 0.21  Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
 3 0.23  Good      E     VS1      56.9    65   327  4.05  4.07  2.31
 4 0.290 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
 5 0.31  Good      J     SI2      63.3    58   335  4.34  4.35  2.75
 6 0.24  Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48
 7 0.26  Very Good H     SI1      61.9    55   337  4.07  4.11  2.53
 8 0.22  Fair      E     VS2      65.1    61   337  3.87  3.78  2.49
 9 0.23  Very Good H     VS1      59.4    61   338  4     4.05  2.39
10 0.3   Good      J     SI1      64      55   339  4.25  4.28  2.73

> df_cv$splits[[1]] %>% assessment()
# A tibble: 4,316 x 10
   carat cut       color clarity depth table price     x     y     z
   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
 1  0.24 Very Good I     VVS1     62.3    57   336  3.95  3.98  2.47
 2  0.3  Very Good J     VS2      62.2    57   357  4.28  4.3   2.67
 3  0.23 Very Good F     VS1      59.8    57   402  4.04  4.06  2.42
 4  0.31 Good      H     SI1      64      54   402  4.29  4.31  2.75
 5  0.32 Good      H     SI2      63.1    56   403  4.34  4.37  2.75
 6  0.22 Premium   E     VS2      61.6    58   404  3.93  3.89  2.41
 7  0.3  Very Good I     SI1      63      57   405  4.28  4.32  2.71
 8  0.24 Premium   E     VVS1     60.7    58   553  4.01  4.03  2.44
 9  0.26 Very Good D     VVS2     62.4    54   554  4.08  4.13  2.56
10  0.86 Fair      E     SI2      55.1    69  2757  6.45  6.33  3.52
# ... with 4,306 more rows

訓練データにはanalysis()で、検証データにはassesment()でアクセスできます。

CVデータが準備できたので、学習と予測を行いましょう。まずは適当なRandom Forestモデルを用意します。

rf = rand_forest(mode = "regression",
                 trees = 50,
                 min_n = 10,
                 mtry = 3) %>%
  set_engine("ranger", num.threads = parallel::detectCores(), seed = 42)

CVデータに前処理レシピを適用するにはmap()prepper()を使います。

> df_cv %>% 
+   mutate(recipes = map(splits, prepper, recipe = rec))
#  10-fold cross-validation 
# A tibble: 10 x 3
   splits               id     recipes     
 * <list>               <chr>  <list>      
 1 <split [38.8K/4.3K]> Fold01 <S3: recipe>
 2 <split [38.8K/4.3K]> Fold02 <S3: recipe>
 3 <split [38.8K/4.3K]> Fold03 <S3: recipe>
 4 <split [38.8K/4.3K]> Fold04 <S3: recipe>
 5 <split [38.8K/4.3K]> Fold05 <S3: recipe>
 6 <split [38.8K/4.3K]> Fold06 <S3: recipe>
 7 <split [38.8K/4.3K]> Fold07 <S3: recipe>
 8 <split [38.8K/4.3K]> Fold08 <S3: recipe>
 9 <split [38.8K/4.3K]> Fold09 <S3: recipe>
10 <split [38.8K/4.3K]> Fold10 <S3: recipe>

prepper()prep()前のレシピを渡すことで、CV訓練データを用いてprep()したレシピを返してくれます。

次に前処理済みのCVデータでモデルを学習します。

> df_cv %>% 
+   mutate(recipes = map(splits, prepper, recipe = rec),
+          fitted = map(recipes, ~ fit(rf, price ~ ., data = juice(.))))
#  10-fold cross-validation 
# A tibble: 10 x 4
   splits               id     recipes      fitted  
 * <list>               <chr>  <list>       <list>  
 1 <split [38.8K/4.3K]> Fold01 <S3: recipe> <fit[+]>
 2 <split [38.8K/4.3K]> Fold02 <S3: recipe> <fit[+]>
 3 <split [38.8K/4.3K]> Fold03 <S3: recipe> <fit[+]>
 4 <split [38.8K/4.3K]> Fold04 <S3: recipe> <fit[+]>
 5 <split [38.8K/4.3K]> Fold05 <S3: recipe> <fit[+]>
 6 <split [38.8K/4.3K]> Fold06 <S3: recipe> <fit[+]>
 7 <split [38.8K/4.3K]> Fold07 <S3: recipe> <fit[+]>
 8 <split [38.8K/4.3K]> Fold08 <S3: recipe> <fit[+]>
 9 <split [38.8K/4.3K]> Fold09 <S3: recipe> <fit[+]>
10 <split [38.8K/4.3K]> Fold10 <S3: recipe> <fit[+]>

学習モデルを用いて検証用データに対して予測を行います。

# ラッパーを作っておく
pred_wrapper = function(split_obj, rec_obj, model_obj, ...) {
  df_pred = bake(rec_obj, assessment(split_obj)) %>%
    bind_cols(predict(model_obj, .)) %>% 
    select(price, .pred)
  
  return(df_pred)
}

> df_cv %>% 
+   mutate(recipes = map(splits, prepper, recipe = rec),
+          fitted = map(recipes, ~ fit(rf, price ~ ., data = juice(.))),
+          predicted = pmap(list(splits, recipes, fitted), pred_wrapper))
#  10-fold cross-validation 
# A tibble: 10 x 5
   splits               id     recipes      fitted   predicted           
 * <list>               <chr>  <list>       <list>   <list>              
 1 <split [38.8K/4.3K]> Fold01 <S3: recipe> <fit[+]> <tibble [4,316 × 2]>
 2 <split [38.8K/4.3K]> Fold02 <S3: recipe> <fit[+]> <tibble [4,316 × 2]>
 3 <split [38.8K/4.3K]> Fold03 <S3: recipe> <fit[+]> <tibble [4,316 × 2]>
 4 <split [38.8K/4.3K]> Fold04 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
 5 <split [38.8K/4.3K]> Fold05 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
 6 <split [38.8K/4.3K]> Fold06 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
 7 <split [38.8K/4.3K]> Fold07 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
 8 <split [38.8K/4.3K]> Fold08 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
 9 <split [38.8K/4.3K]> Fold09 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>
10 <split [38.8K/4.3K]> Fold10 <S3: recipe> <fit[+]> <tibble [4,315 × 2]>

pmap()は入力列が複数の場合のmap()の拡張です。

予測結果を用いて検証用データで精度を計算して完了です。

> df_cv %>% 
+   mutate(recipes = map(splits, prepper, recipe = rec),
+          fitted = map(recipes, ~ fit(rf, price ~ ., data = juice(.))),
+          predicted = pmap(list(splits, recipes, fitted), pred_wrapper),
+          evaluated = map(predicted, metrics, price, .pred))
#  10-fold cross-validation 
# A tibble: 10 x 6
   splits               id     recipes      fitted   predicted            evaluated       
 * <list>               <chr>  <list>       <list>   <list>               <list>          
 1 <split [38.8K/4.3K]> Fold01 <S3: recipe> <fit[+]> <tibble [4,316 × 2]> <tibble [3 × 3]>
 2 <split [38.8K/4.3K]> Fold02 <S3: recipe> <fit[+]> <tibble [4,316 × 2]> <tibble [3 × 3]>
 3 <split [38.8K/4.3K]> Fold03 <S3: recipe> <fit[+]> <tibble [4,316 × 2]> <tibble [3 × 3]>
 4 <split [38.8K/4.3K]> Fold04 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
 5 <split [38.8K/4.3K]> Fold05 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
 6 <split [38.8K/4.3K]> Fold06 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
 7 <split [38.8K/4.3K]> Fold07 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
 8 <split [38.8K/4.3K]> Fold08 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
 9 <split [38.8K/4.3K]> Fold09 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>
10 <split [38.8K/4.3K]> Fold10 <S3: recipe> <fit[+]> <tibble [4,315 × 2]> <tibble [3 × 3]>

ハイパーパラメータのサーチ

Cross Validationで精度を計算する方法がわかったので、これを用いてハイパーパラメータのサーチを行います。 今回はmin_nのみグリッドサーチすることにします。

grid_min_n = c(5, 10, 15)
df_result_cv = tibble() # 結果をいれる
for (n in grid_min_n) {
  
  rf = rand_forest(mode = "regression",
                   trees = 50,
                   min_n = n,
                   mtry = 3) %>%
    set_engine("ranger", num.threads = parallel::detectCores(), seed = 42)
  
  tmp_result = df_cv %>% 
    mutate(recipes = map(splits, prepper, recipe = rec),
           fitted = map(recipes, ~ fit(rf, price ~ ., data = juice(.))),
           predicted = pmap(list(splits, recipes, fitted), pred_wrapper),
           evaluated = map(predicted, metrics, price, .pred))
  
  df_result_cv = df_result_cv %>% 
    bind_rows(tmp_result %>% 
                select(id, evaluated) %>% 
                mutate(min_n = n))
}

> df_result_cv %>% 
+   unnest() %>% 
+   group_by(min_n, .metric) %>% 
+   summarise(mean(.estimate))
# A tibble: 9 x 3
# Groups:   min_n [?]
  min_n .metric `mean(.estimate)`
  <dbl> <chr>               <dbl>
1     5 mae                0.0671
2     5 rmse               0.0932
3     5 rsq                0.992 
4    10 mae                0.0679
5    10 rmse               0.0940
6    10 rsq                0.991 
7    15 mae                0.0686
8    15 rmse               0.0947
9    15 rsq                0.991 

どうやらmin_n = 5の精度が最も良さそうです。 このパラメータを用いて訓練データ全体でモデルを再学習し、予測精度を確認します。 parsnipで作ったモデルはupdate()でパラメータを更新できます。

rf_best = update(rf, min_n = 5)

rec_preped = rec %>% 
  prep(df_train)

fitted = rf_best %>% 
  fit(price ~ ., data = juice(rec_preped))


df_test_baked = bake(rec_preped, df_test)


> df_test_baked %>% 
+   bind_cols(predict(fitted, df_test_baked)) %>% 
+   metrics(price, .pred)
# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard      0.0905
2 rsq     standard      0.992 
3 mae     standard      0.0665

予測結果と実際の値もプロットして確認します。

df_test_baked %>% 
  bind_cols(predict(fitted, df_test_baked)) %>% 
  ggplot(aes(.pred, price)) +
  geom_hex() +
  geom_abline(slope = 1, intercept = 0) + 
  coord_fixed() +
  scale_x_continuous(breaks = seq(6, 10, 1), limits = c(5.8, 10.2)) + 
  scale_y_continuous(breaks = seq(6, 10, 1), limits = c(5.8, 10.2)) + 
  scale_fill_viridis_c() + 
  labs(x = "Prediction", y = "Truth") + 
  theme_minimal(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black", size =12))

f:id:dropout009:20190109213420p:plain

大きく外したり系統的なバイアスも見られないので、うまく予測できていると言えるのではないかと思います。

まとめ

tidymodels配下のパッケージ、特にpurrrrsampleを用いることでtidyにCross Validationを行うことができました。

参考

変数重要度とPartial Dependence Plotで機械学習モデルを解釈する

はじめに

RF/GBDT/NNなどの機械学習モデルは古典的な線形回帰モデルよりも高い予測精度が得られる一方で、インプットとアウトプットの関係がよくわからないという解釈性の問題を抱えています。 この予測精度と解釈性のトレードオフでどちらに重点を置くかは解くべきタスクによって変わってくると思いますが、私が仕事で行うデータ分析はクライアントの意思決定に繋げる必要があり、解釈性に重きを置いていることが多いです。

とはいえ機械学習モデルの高い予測精度は惜しく、悩ましかったのですが、学習アルゴリズムによらずモデルに解釈性を与えられる手法が注目され始めました。 本記事では変数重要度とPDP/ICE Plot (Partial Dependence/Individual Conditional Expectation)を用いて所謂ブラックボックスモデルを解釈する方法を紹介します。

※この記事で紹介する変数重要度やPD/ICEを含んだ、機械学習の解釈手法に関する本を書きました!

モデルの学習

まずはパッケージを読み込みます。

library(tidyverse)
library(tidymodels)

library(pdp) # partial dependence plot
library(vip) # variable importance plot

例によってdiamondsデータを使用し、Rondom Forestでダイヤの価格を予測するモデルを作ります。 tidymodelsの使い方は以前の記事をご覧下さい。

set.seed(42)

df_split = initial_split(diamonds,  p = 0.9)

df_train = training(df_split)
df_test  = testing(df_split)


rec = recipe(price ~ carat + clarity + color + cut, data = df_train) %>% 
  step_log(price) %>% 
  step_ordinalscore(all_nominal()) %>% 
  prep(df_train)

df_train_baked = juice(rec)
df_test_baked = bake(rec, df_test)

fitted = rand_forest(mode = "regression", 
                     trees = 100,
                     min_n = 5,
                     mtry = 2) %>% 
  set_engine("ranger", num.threads = parallel::detectCores(), seed = 42) %>% 
  fit(price ~ ., data = juice(rec))

> fitted
parsnip model object

Ranger result

Call:
 ranger::ranger(formula = formula, data = data, mtry = ~2, num.trees = ~100,      min.node.size = ~5, num.threads = ~parallel::detectCores(),      seed = ~42, verbose = FALSE) 

Type:                             Regression 
Number of trees:                  100 
Sample size:                      48547 
Number of independent variables:  4 
Mtry:                             2 
Target node size:                 5 
Variable importance mode:         none 
Splitrule:                        variance 
OOB prediction error (MSE):       0.01084834 
R squared (OOB):                  0.9894755 

変数重要度

モデルの解釈において最も重要なのは、どの変数がアウトプットに強く影響し、どの変数は影響しないのかを特定すること、つまり変数重要度を見ることだと思います。 たとえば特に重要な変数を見定めて施策を打つことでアウトプットをより効率よく改善できる可能性が高まりますし*1、ありえない変数の重要度が高く出ていればデータやモデルが間違っていることに気づけるかもしれません。

tree系のアルゴリズムだとSplit(分割の回数)やGain(分割したときの誤差の減少量)などで変数重要度を定義することもできますが、ここではアルゴリズムに依存しない変数重要度であるPermutationベースのものを紹介します。 Permutationベースの変数重要度のコンセプトは極めて直感的で、「ある変数の情報を壊した時にモデルの予測精度がすごく落ちるならそれは大事な変数だ」というものです。

f:id:dropout009:20190107002315p:plain

たとえば今回のモデルでcaratの変数重要度を見たいとします。

  1. 訓練データを使って学習済みモデルを用意し、テストデータに対する予測精度を出します。
  2. テストデータのcarat列の値をシャッフルしますし、シャッフル済みテストデータを使って同様に予測を行います。
  3. 2つの予測精度を比較し、予測精度の減少度合いを変数重要度とします。

carat列のシャッフルによって、もしcaratが重要な変数ならモデルは的はずれな予測をするようになりますし、逆に重要でないなら特に影響は出ないはずで、これはある意味においての変数の重要度を捉えらていると思います。

また、この手法はあくまで予測の際に用いるデータをシャッフルしていて、学習をやり直しているわけではないため、計算が軽いことも利点です。

さて、実際にPermutationベースの変数重要度を計算してみましょう。 パッケージはvipを使います。変数重要度が計算できるパッケージは他にもいくつかありますが、vipはPermutationベースの変数重要度が計算でき、次に見るPartial Dependence Plotで利用するパッケージpdpと同じシンタックスが使えます(作者が同じなので)

# objectとnewdataを受けて予測結果をベクトルで返す関数を作っておく必要がある
pred_wrapper = function(object, newdata) {
  return(predict(object, newdata) %>% pull(.pred))
}


fitted %>% 
  vip(method = "permute",  # 変数重要度の計算方法
      pred_fun = pred_wrapper, 
      target = df_test_baked %>% pull(price),
      train = df_test_baked %>% select(-price),
      metric = "rsquared", # 予測精度の評価関数) 

f:id:dropout009:20190106230934p:plain

一番手っ取り早い方法はvip()関数を使った可視化です(グラフを自分で細かくいじりたい場合はvi()での計算結果を自分で可視化します)。 caratの変数重要度がその他の変数と比べてダントツであることが見て取れます。

Partial Dependence Plot

各変数の重要度がわかったら、次に行うべきは重要な変数とアウトカムの関係を見ることだと思います。 ただ、一般にブラックボックスモデルにおいてインプットとアウトカムの関係は非常に複雑で、可視化することは困難です。

そこで、複雑な関係を要約する手法の一つにPartial Dependence Plot(PDP)があります。PDPは興味のある変数以外の影響を周辺化して消してしまうことで、インプットとアウトカムの関係を単純化しようというものです。

特徴量を  x とした学習済みモデル  \hat{f}(x) があるとします。 x を興味のある変数  x_s とその他の変数 x_c に分割し、以下のpartial dependence function

 \begin{align} \hat{f_{s}}(x_s) = E_c\left[ \hat{f}(x_s, x_c) \right] = \int \hat{f}(x_s, x_c) p(x_c) d x_c \end{align}

を定義し、これを

\begin{align} \bar{f_{s}}(x_s) = \frac{1}{N}\sum_i\hat{f}\left(x_s, x_c^{(i)}\right) \end{align}

で推定します。具体的には以下のようなことをやっています。

f:id:dropout009:20190107002309p:plain

なお、PDPは平均をとりますが、平均を取らずに各 iについて個別にプロットしたものをIndividual Conditional Expectation (ICE) Plotと呼びます。 特に変数間の交互作用がある場合にPDPでは見失ってしまう関係を確認することができます。

それでは、特に重要だった変数caratについて実際にPDPを作成してみましょう。パッケージはpdpを使います。

ice_carat = fitted %>% 
  partial(pred.var = "carat",
          pred.fun = pred_wrapper,
          train = df_test_baked, 
          type = "regression")

ice_carat %>% 
  autoplot()

f:id:dropout009:20190107011226p:plain

partial()でPD/ICE Plotの元となるデータを計算できます。細部にこだわらなければautoplot()で可視化するのが一番手っ取り早いです。 赤い線がPDP、黒い線がICEになります。caratlog(price)に与える影響が徐々に逓減していく様子が見て取れます。 事前に関数形の仮定を置かずにインプットとアウトカムの非線形な関係を柔軟に捉えることができるのがRF+PDPの強力な点かと思います。

まとめ

アルゴリズムによらずモデルを解釈する手法としては他にもALE、SHAP value、LIMEなどがあり、次回以降に紹介できればと思っています。

参考

*1:変数間の関係が相関関係なのか因果関係なのかという問題はありますが

tidymodelsによるtidyな機械学習(その1:データ分割と前処理から学習と性能評価まで)

目次

※この記事をベースにした2019年12月7日に行われたJapan.R 2019での発表資料は以下になります。 tidymodelsによるtidyな機械学習 - Speaker Deck

はじめに

本記事ではtidymodelsを用いたtidyな機械学習フローを紹介したいと思います。

tidyverseはデータハンドリングと可視化のためのメタパッケージでしたが、tidymodelstydyverseにフィットするやり方で統計モデリング/機械学習をするためのメタパッケージになります。

tidymodels配下のパッケージは量が多く使い所が限られているパッケージも多いため、一度に全ては紹介できません。 ですので、今回は典型的な

  1. 訓練データとテストデータの分割
  2. 特徴量エンジニアリング
  3. モデルの学習
  4. モデルの精度評価

という機械学習フローを想定し、各パッケージの機能と連携をコンパクトにまとめたいと思います。

各パッケージの用途は以下になります。

  • rsample:訓練データ/テストデータの分割
  • recipe:特徴量エンジニアリング
  • parsnip:異なるパッケージのアルゴリズムを同じシンタックスで扱えるようにするラッパー
  • yardstic:モデルの精度評価

tidyな機械学習フロー

tidyversetidymodelsを読み込みます。

library(tidyverse)
library(tidymodels)

訓練データとテストデータの分割

まず第一に大元のデータを訓練データとテストデータに分割しましょう。

データはdiamondsを使って、ダイヤの価格を予測することにします。

set.seed(42)

df_split = initial_split(diamonds,  p = 0.8)

df_train = training(df_split)
df_test  = testing(df_split)

データ分割の役割はrsampleが担っています。

initial_split()で分割するデータと、訓練データが占める割合を指定します。

training()で訓練データを、testing()でテストデータを抽出できます。

特徴量エンジニアリング

次に特徴量エンジニアリングのためのレシピを作成します。作成したレシピを実際にデータに適用することで、処理済みのデータが作成されます。

どのような特徴量エンジニアリングが必要かは機械学習モデルに依存しますが、今回はRandom Forestを用いるための最低限の前処理を行います。

rec = recipe(price ~ carat + color + cut + clarity + depth, data = df_train) %>% 
  step_log(price) %>% 
  step_ordinalscore(all_nominal())

> rec
Data Recipe

Inputs:

      role #variables
   outcome          1
 predictor          9

Operations:

Log transformation on price
Scoring for all_nominal()

特徴量エンジニアリングのための関数はrecipeにまとめられています。

まずはrecipe()で用いる特徴量とデータを指定し、そこに処理をパイプで重ねていきます。

関数は全てstep_*()で統一されていて、今回用いたstep_log()は指定した変数の対数をとる関数で、step_ordinalscore()は順序のあるカテゴリカル変数を順序(1, 2, 3..)に変換します。 (tree系のモデルではなく線形モデルを用いる場合はstep_dummy()でダミー変数にできます。 *1

all_nominal()は全てのカテゴリカル変数を指定しています。他にもたとえばall_numeric()は全ての連続変数を指定できますし、starts_with()などのtidyverseではお馴染みの関数で指定することもできます。詳細は?selectionsをご確認下さい。

レシピができたので実際にデータに適用しましょう。

rec_preped = rec %>% 
  prep(df_train)

df_train_baked = rec_preped %>% 
  juice() 

> df_train_baked
# A tibble: 43,153 x 10
   carat   cut color clarity depth table     x     y     z price
   <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 0.23      5     2       2  61.5    55  3.95  3.98  2.43  5.79
 2 0.21      4     2       3  59.8    61  3.89  3.84  2.31  5.79
 3 0.23      2     2       5  56.9    65  4.05  4.07  2.31  5.79
 4 0.290     4     6       4  62.4    58  4.2   4.23  2.63  5.81
 5 0.31      2     7       2  63.3    58  4.34  4.35  2.75  5.81
 6 0.24      3     7       6  62.8    57  3.94  3.96  2.48  5.82
 7 0.24      3     6       7  62.3    57  3.95  3.98  2.47  5.82
 8 0.26      3     5       3  61.9    55  4.07  4.11  2.53  5.82
 9 0.22      1     2       4  65.1    61  3.87  3.78  2.49  5.82
10 0.23      3     5       5  59.4    61  4     4.05  2.39  5.82
# ... with 43,143 more rows

prep()に訓練データを指定し、訓練データの変換を行います。

さらにjuice()で変換されたデータを抽出することができます。

priceに対数がとられてカテゴリカル変数が(cut, color, clarity)が順序になっているのがわかるかと思います。

次に同じレシピを用いてテストデータも変換します。

df_test_beked = rec_preped %>% 
  bake(df_test)

> df_test_baked
# A tibble: 10,787 x 10
   carat   cut color clarity depth table price     x     y     z
   <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1  0.22     4     3       3  60.4    61  5.83  3.88  3.84  2.33
 2  0.3      3     7       3  62.7    59  5.86  4.21  4.27  2.66
 3  0.23     3     2       4  63.8    55  5.86  3.85  3.92  2.48
 4  0.23     3     1       4  60.5    61  5.88  3.96  3.97  2.4 
 5  0.23     3     2       5  59.5    58  6.00  4.01  4.06  2.4 
 6  0.23     2     3       5  58.2    59  6.00  4.06  4.08  2.37
 7  0.33     5     6       2  61.8    55  6.00  4.49  4.51  2.78
 8  0.33     5     7       3  61.1    56  6.00  4.49  4.55  2.76
 9  0.3      3     6       3  62.6    57  6.00  4.25  4.28  2.67
10  0.3      4     1       3  62.6    59  6.31  4.23  4.27  2.66
# ... with 10,777 more rows

レシピを新しいデータに適用する場合はbake()を使います*2

モデルの学習

recipeで訓練データ/テストデータをモデルに投入できる形に変換できたので、次はparsnipで学習モデルを定義しましょう。

Rで機械学習を行う際、学習アルゴリズムごとにパッケージが分かれていて、それ故シンタックスが統一されておらず、各々のシンタックスを覚えるのが大変でした(たとえば入力はformulaなのかmatrixなのかなど)。

これを統一的なシンタックスで扱えるようにするメタパッケージとして、たとえばcaretがありましたが、最近よりtidyverseに馴染むパッケージとしてparsnipが開発されました*3。 ちなみに開発者はcaretと同じくMax Kuhnです。

rf = rand_forest(mode = "regression", 
                 trees = 50,
                 min_n = 10,
                 mtry = 3) %>% 
  set_engine("ranger", num.threads = parallel::detectCores(), seed = 42)

> rf
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 3
  trees = 50
  min_n = 10

Engine-Specific Arguments:
  num.threads = parallel::detectCores()
  seed = 42

Computational engine: ranger 

rand_forest()は学習アルゴリズムとしてRandom Forestを使用します。今回は回帰問題なのでmode=regressionとしています。 rand_forest()で指定できるハイパーパラメータは3つあり、それぞれ

  • trees:作成する決定木の数。
  • min_n:ノードに最低限含まれるサンプル数。
  • mtry:各ノードで分割に用いる特徴量の数。

となります。今回は適当に指定しました。

さらにset_engine()でRandom Forestに用いるパッケージを指定でき、ここではranger を指定しています。 set_engine()内でranger固有のパラメータを指定でき、コア数とシードを指定しました。

モデルが定義できたので、学習を行います。これはfit()にモデル式と入力データを指定するだけで完結します。

fitted = rf %>% 
  fit(price ~ ., data = df_train_baked)

> fitted
parsnip model object

Ranger result

Call:
 ranger::ranger(formula = formula, data = data, mtry = ~3, num.trees = ~50,      min.node.size = ~10, num.threads = ~parallel::detectCores(),      seed = ~42, verbose = FALSE) 

Type:                             Regression 
Number of trees:                  50 
Sample size:                      43153 
Number of independent variables:  9 
Mtry:                             3 
Target node size:                 10 
Variable importance mode:         none 
Splitrule:                        variance 
OOB prediction error (MSE):       0.009112244 
R squared (OOB):                  0.9911622 

モデルの精度評価

学習が済んだので、最後に予測精度を評価しましょう。 テストデータへの予測はpredict() で行います。

df_result = df_test_baked %>% 
  select(price) %>% 
  bind_cols(predict(fitted, df_test_baked))
  
> df_result
# A tibble: 10,787 x 2
   price .pred
   <dbl> <dbl>
 1  5.83  5.95
 2  5.86  5.98
 3  5.86  5.97
 4  5.88  5.97
 5  6.00  6.00
 6  6.00  5.99
 7  6.00  6.14
 8  6.00  6.14
 9  6.00  6.03
10  6.31  6.33
# ... with 10,777 more rows

予測結果に対してyardsticの評価指標関数を適用することでモデルの性能が評価できます。 yardsticには様々な評価指標関数が含まれており、個別で指定することもできますが、matrics()は回帰モデルの結果に対してRMSE、MAE、R2をまとめて計算してくれます。

> df_result %>% 
+   metrics(price, .pred)
# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard      0.0910
2 rsq     standard      0.992 
3 mae     standard      0.0671

まとめ

本記事ではtidymodelsによる

  1. rsampleで訓練データ/テストデータの分割し、
  2. recipeで特徴量エンジニアリングを行い、
  3. parsnipでモデルを定義・学習し、
  4. yardsticで性能を評価

する、tidyな機械学習フローを紹介しました。

各パッケージにはまだまだ多くの機能があり、今回紹介できたのは最低限の機能のみです。

次回はtidymodelsで交差検証を行い、ハイパーパラメータをチューニングする方法をご紹介できればと思っています。

参考文献

*1: recipeの他のstep関数やselectorは以下をご確認下さい。Function reference • recipes

*2:今回tree系のモデルなので標準化はやっていませんが、たとえば標準化を行う場合は訓練データの平均と分散を用いてテストデータの標準化が行われます。

*3:parsnipが対応しているパッケージは以下を御覧ください:List of Models • parsnip

Ridge回帰の使い所を考える

はじめに

本記事では回帰係数の推定方法としてのRidge回帰(L2正則化)の使い所について考えたいと思います。
結論から言うと、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方ができそうです。

OLS推定量とRidge推定量

簡単のために、以下の線形回帰モデルを考えます。

\begin{align}
Y = X\beta + \epsilon
\end{align}

ここで、 X N \times Kの行列で、 \beta K\times 1の列ベクトルとします。  Nはサンプルサイズ、  Kは特徴量の数になります。
このとき、OLS推定量は二乗誤差を最小にするような bで与えられます。
\begin{align}
\widehat{\beta}_{\mathrm{OLS}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) \\
&= \left(X^\top X \right)^{-1}X^\top Y
\end{align}

標準的な仮定のもとで、OLS推定量は推定量 \widehat{\beta}_{\mathrm{OLS}} は真の値 \betaの一致推定量となることが知られています。*1

これに対して、 b のとれる値に制限をかけることで過学習を抑えようとするのがRidge回帰になります。
具体的にはOLSの目的関数に正則化\lambda||b||^2を付け加えたものを最小化します。\lambda正則化の強さを表すハイパーパラメータになります。

\begin{align}
\widehat{\beta}_{\mathrm{Ridge}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) + \lambda||b||^2 \\
&= \left(X^\top X + \lambda I_K \right)^{-1}X^\top Y
\end{align}

ここで   I_K K\times K単位行列です。
 \widehat{\beta}_{\mathrm{OLS}} \widehat{\beta}_{\mathrm{Ridge}} を比較すれば明らかですが、両者の違いは  \lambda I_Kのみです。
なので、  \lambda = 0のときOLS推定量とRidge推定量は完全に一致します。反対に  \lambdaを大きくするとRidge推定量は0方向に引っ張られます。
つまり、Ridge回帰はパラメータが極端な値を取らないよう安定させる(過学習を抑える)というリターンを得るために、(多少の)バイアスがかかるという犠牲を払うモデルであると言えます。
これは機械学習の文脈でバイアスとバリアンスのトレードオフとしてよく知られています。

シミュレーション

それでは実際にOLSとRidge回帰の挙動の違いについてシミュレーションで確認してみましょう。
まずは必要なパッケージを読み込みます。*2

import pandas as pd
import numpy as np

from sklearn.linear_model import LinearRegression, Ridge,  RidgeCV
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_spd_matrix #分散共分散行列を作る

from plotnine import * # 可視化用。


ちょっと極端な方が結果がクリアに見えると思うので、サンプルサイズは100しかないんだけど80個の特徴量について回帰係数を推定したい、というシチュエーションを想定することにします。
下記のコードではランダムにデータを生成してOLSとRidge回帰で回帰係数を推定する、ということを繰り返しています。

np.random.seed(42)
I = 1000
N = 100
K = 80

beta = np.random.uniform(-1, 1, K)

betas = np.zeros((I, 2, K))
for i in range(I):
    X = np.random.multivariate_normal(np.zeros(K), make_spd_matrix(K), N)
    X = StandardScaler().fit_transform(X)
    U = np.random.normal(0, 1, N)
    Y = X @ beta + U
    
    if i == 0: #一回目のループでRidgeのハイパーパラメータを決める
        ridge_alpha = RidgeCV(alphas=np.random.uniform(0.1, 20, 100), cv=5).fit(X, Y).alpha_
        print(f'ridge: {ridge_alpha}')
    
    beta_ols = LinearRegression().fit(X, Y).coef_
    
    beta_ridge = Ridge(alpha=ridge_alpha).fit(X, Y).coef_
    
    betas[i, :, :] = np.concatenate([beta_ols[:, np.newaxis], beta_ridge[:, np.newaxis]], axis=1).T
    
    if i % 200 == 0:
        print(i)

# 可視化のためにシミュレーションの結果をデータフレームとして持ち直す
df = None
for i in range(I):
    tmp = (pd.DataFrame(data = betas[i, :, :].T, columns=['OLS', 'Ridge'])
           .assign(Truth = beta, k = lambda df: df.index, i = i))
    
    if df is None:
        df = tmp
    else:
        df = df.append(tmp)
        
df = df.melt(id_vars=['i', 'k', 'Truth'], value_vars=['OLS', 'Ridge'], value_name='coefficient', var_name='method')


シミュレーションの結果を可視化してみましょう。

g = (ggplot(aes('factor(k)', 'coefficient', fill='method', color='method'), df[df.k.between(1, 10)]) +
     geom_hline(yintercept=0, size=1, color='gray') +
     geom_violin(alpha=0.5) +
     geom_point(aes('factor(k)', 'Truth'), size=4, color='black', fill='black') +
     scale_fill_brewer(name='Method', type='qual',palette='Dark2') +
     scale_color_brewer(name = 'Method', type='qual',palette='Dark2') + 
     lims(y=(-3, 3)) + 
     labs(x='beta_hat(k)', y='Coefficient') + 
     theme_bw() +
     theme(figure_size=(15, 7)))

print(g)

f:id:dropout009:20190102204239p:plain



図は初めの10個の回帰係数を取り出して、シミュレーション1000回分の回帰係数の分布を確認しています。
縦軸が回帰係数の値、横軸が回帰係数の番号です。
真ん中の黒い点が真の回帰係数の値で、緑がOLSで求めた回帰係数の分布、オレンジがRidge回帰での分布になります。


図からは各手法で求めた回帰係数について、以下の性質が見て取れます。

  • OLSの場合はバイアスがなく真の値を中心に分布しているが、その分散は大きく、見当はずれの値を取ることもある
  • Ridge回帰の場合はゼロ方向にバイアスがかかっている。その一方で回帰係数は非常に安定しており、真の値に近い部分に密度高く分布しており、極端な値を取ることもない。

Ridge回帰の使いどころ

数式とシミュレーションの結果から、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方が考えられそうです。*3
特に、特徴量のインパクトを見える際により保守的な推定をしたい場合には有効と言えそうです。Ridge回帰は極端な回帰係数をとらず、またRidge回帰のバイアスは0方向にかかるので、Ridge回帰ではより控え目なインパクトが推定されると言えるかと思います。

*1:https://www.ssc.wisc.edu/~bhansen/econometrics/

*2:ちなみにplotnineはRの代表的な可視化パッケージであるggplot2と同じ記法が使えるpython用のパッケージです。若干制限もありますが、ggplot2に慣れたRユーザーにとっては使いやすいかと思います。

*3:サンプルサイズが十分な状況ではOLSで一致推定を行うのが妥当かと思います。

purrrとbroomで複数の回帰モデルを効率的に管理する

私は探索的にデータを見てく段階では、可視化に加えて複数の回帰モデルを作成して比較をする、ということをよくやっています。 モデルの数が少ない場合は個別にモデルを作成してsummary()で見ていく事もできますが、モデルの数が増えるにつれてそのやり方では管理が難しくなってきます。 そこで、本記事では、purrrmap()broomtidy(), glance()を用いて複数の回帰モデルを効率的に扱う方法を紹介したいと思います。

まずはライブラリを読み込みます。tidyverseはおなじみのデータハンドリングと可視化のためのパッケージ群です。tidymodelsモデリングをtidyなやり方で統一的に扱えるようにするパッケージ群になります。今回はbroomのみ用いますが、後日他のパッケージの紹介記事も書ければと思っています。

library(tidyverse)
library(tidymodels)

# 可視化用
theme_scatter = theme_minimal() + 
    theme(panel.grid.minor = element_blank(), 
          axis.text  = element_text(color = "black")) 

theme_minimal2 = theme_scatter + 
    theme(panel.grid.major.x = element_blank()) 

データセットはこれまたおなじみのdiamondを使います。 詳細は?diamondsを見て頂ければと思いますが、ダイヤの重さ(carat)や透明度(clarity)とその値段(price)などが入ったデータセットになります。 今回はこのデータセットを用いて、ダイヤの属性がダイヤの値段にどんな影響を与えるのかを探索することにします。

可視化によるデータ分析

まずは非説明変数である価格の分布を見てみましょう。

df = diamonds #面倒なので名前をdfに

df %>% 
    ggplot(aes(price)) +
    geom_histogram(fill = "#56B4E9", color = "white") +
    theme_minimal2

f:id:dropout009:20190102144514p:plain

ものすごく歪んだ分布をしています。 歪んだ分布への簡便な対応策として、今回は対数を取ることにします。

df %>% 
    ggplot(aes(log(price))) +
    geom_histogram(fill = "#56B4E9", color = "white") +
    theme_minimal2

f:id:dropout009:20190102144521p:plain

きれいな正規分布とまではいきませんが先程よりは中心によった分布になりました。価格は対数を取って分析することにします。

次にダイヤの属性と値段の関係を見ていきましょう。 たとえばダイヤが重ければ重いほど(大きければ大きいほど)値段は高くなりそうです。

df %>% 
    sample_frac(0.1) %>% #データが多いので減らす
    ggplot(aes(carat, log(price))) +
    geom_point(color = "#0072B2", alpha = 0.5) +
    theme_scatter

f:id:dropout009:20190102144526p:plain

実際に可視化してみると、どうやらこの仮説は正しそうです。 ただ、caratlog(price)の関係は非線形に見えます。 caratlog(price)に与えるインパクトはcaratが大きくなるにつれて逓減するという関係を反映するために、caratにも対数をとってみましょう。

df %>% 
    sample_frac(0.1) %>% #データが多いので減らす
    ggplot(aes(log(carat), log(price))) +
    geom_point(color = "#0072B2", alpha = 0.5) +
    theme_scatter

f:id:dropout009:20190102144533p:plain

両変数に対数をとることで線形の関係が構築できました! これなら線形回帰でうまくモデリングできそうです。

同様に透明度と価格の関係も見てみましょう。 直感的には透明度が高ければ高いほど価格が高くなる傾向がありそうに思えます。

df %>% 
    ggplot(aes(clarity, log(price))) +
    geom_boxplot(fill = "#56B4E9") +
    theme_minimal2

f:id:dropout009:20190102144540p:plain 直感に反する結果となりました。上のグラフは右に行くほど透明度のランクが高いのですが、透明度が高いほどダイヤが安くなっています。

これは一体どういうことでしょうか? 透明度とカラット数の関係をみることでこの謎が解けます。

df %>% 
    ggplot(aes(clarity, log(carat))) +
    geom_boxplot(fill = "#56B4E9") +
    theme_minimal2

f:id:dropout009:20190102144545p:plain 透明度が高い場合にはカラット数が小さくなる傾向が見て取れます。 僕はダイヤモンドに詳しくないのですが、大きくて透明なダイヤを作るのは難しいということのようです。*1

つまり、ここでは「ダイヤのカラット数と価格には正の相関がある一方で、カラット数と透明度には負の相関があり、結果として透明度と価格に負の相関があるように見えてしまう」という典型的な交絡の問題が起きています。 この場合は、重回帰分析でカラット数の影響を取り除くことで透明度と価格の正しい関係を見ることができます。 今回のように、単純に一つ一つの説明変数と被説明変数の関係を見ているだけでは間違った結論を下していまう可能性があるので注意が必要です。

モデリングによるデータ分析

ここまでの可視化で得られた仮説を回帰モデルで確かめてみましょう。

回帰モデルを作る前に、カテゴリカル変数がordere=TRUEになっているとlm()の挙動が面倒なので、factor型の変数は全てordered=FALSEにしておきます。

df_input = df %>% 
    mutate_if(is.ordered, factor, ordered = FALSE)

mutate_if()は複数の列に同じ処理を行う際に便利な関数です。 第一引数で条件を指定して(is.ordered)、条件に当てはまった列にのみ第二引数の関数を適用します(factor)。関数のオプションは後ろにくっつけて指定すればOKです(ordered = FALSE)。

複数のモデルを当てはめる場合

今回は

  1. 価格を透明度のみで説明するモデル
  2. 透明度とカラット数で説明するモデル
  3. 上と同じだが、カラット数に対数をとったモデル

の3つのモデルを作成してみます。

formulas = c(log(price) ~ clarity,
             log(price) ~ clarity + carat,
             log(price) ~ clarity + log(carat)) %>% 
    enframe("model_no", "formula")

formulas 

# A tibble: 3 x 2
  model_no formula      
     <int> <list>       
1        1 <S3: formula>
2        2 <S3: formula>
3        3 <S3: formula>

enframe()はベクトルをデータフレームにしてくれる関数です。データフレームとして持っていたほうが後々管理がしやすいと思うので、変換しておきました。

上のデータフレームにmap()を適用することで複数の回帰モデルを一気に作成できます。

df_result = formulas %>% 
    mutate(model = map(formula, lm, data = df_input), #(1)
           tidied = map(model, tidy), #(2)
           glanced = map(model, glance)) #(3)

df_result

# A tibble: 3 x 5
  model_no formula       model    tidied           glanced          
     <int> <list>        <list>   <list>           <list>           
1        1 <S3: formula> <S3: lm> <tibble [8 × 5]> <tibble [1 × 11]>
2        2 <S3: formula> <S3: lm> <tibble [9 × 5]> <tibble [1 × 11]>
3        3 <S3: formula> <S3: lm> <tibble [9 × 5]> <tibble [1 × 11]>

(1)ではmap()を用いて「formula列の一つ一つの値を引数としてlm()を実行する」ということをやっています。lm()はデータも指定する必要がありますが、map()mutate_if()同様後ろにくっつけて指定できます。 さらに、(2), (3)では推定されたモデルに対して、それぞれtidy()glance()を適用しています。

分析結果をデータフレームとして持つことで、どのモデルがどの結果に対応するかを間違うことなく効率的に管理することができます。

tidy()は回帰モデルの係数をtidyなデータフレームとして持ってきてくれる関数です。データフレームの中のデータフレームを取り出すためにunnest()を使います。

df_coef = df_result %>% 
    select(model_no, tidied) %>% 
    unnest() %>% 
    mutate_if(is.double, round, digits=2) 

df_coef

# A tibble: 26 x 6
   model_no term        estimate std.error statistic p.value
      <int> <chr>          <dbl>     <dbl>     <dbl>   <dbl>
 1        1 (Intercept)    8.03       0.04    221.         0
 2        1 claritySI2     0.14       0.04      3.69       0
 3        1 claritySI1    -0.18       0.04     -4.81       0
 4        1 clarityVS2    -0.26       0.04     -7.08       0
 5        1 clarityVS1    -0.3        0.04     -7.99       0
 6        1 clarityVVS2   -0.5        0.04    -12.8        0
 7        1 clarityVVS1   -0.7        0.04    -17.7        0
 8        1 clarityIF     -0.62       0.04    -14.4        0
 9        2 (Intercept)    5.36       0.01    372.         0
10        2 claritySI2     0.570      0.01     40.1        0
# ... with 16 more rows

3つのモデルを比較してみましょう。 モデルを横に並べるためにspread()を使います。 そのままだと回帰係数がアルファベット順になってしまうので、並び順を維持するためにfct_inorder()を使っています(出てきた順に並ぶ)。

df_coef %>% 
    mutate(term = fct_inorder(term)) %>% 
    select(model_no, term, estimate) %>% 
    spread(model_no, estimate)

# A tibble: 10 x 4
   term           `1`    `2`   `3`
   <fct>        <dbl>  <dbl> <dbl>
 1 (Intercept)   8.03  5.36   7.77
 2 claritySI2    0.14  0.570  0.48
 3 claritySI1   -0.18  0.72   0.62
 4 clarityVS2   -0.26  0.82   0.78
 5 clarityVS1   -0.3   0.86   0.82
 6 clarityVVS2  -0.5   0.93   0.98
 7 clarityVVS1  -0.7   0.92   1.03
 8 clarityIF    -0.62  1      1.11
 9 carat        NA     2.08  NA   
10 log(carat)   NA    NA      1.81

モデル1は透明度のみのモデルであり、カラット数の影響を取り除いていないので、透明度と価格には負の関係があるように見えています。 その一方で、モデル2はカラット数を変数として加えることで影響を取り除いた上での透明度と価格の影響を見ています。モデル2の係数を見ると透明度が高いほど価格が高くなる傾向が見て取れ、こちらのほうがより尤もらしい結果であると言えそうです。

モデル3はカラット数に対数をとったモデルになります。 モデル2もモデル3も透明度と価格の関係には大きな変化はなさそうです。

glance()はモデルの性能をtidyなデータフレームとしてまとめてくれる関数です。

df_result %>% 
    select(model_no, glanced) %>% 
    unnest() %>% 
    mutate_if(is.double, round, digits=2) 

# A tibble: 3 x 12
  model_no r.squared adj.r.squared sigma statistic p.value    df  logLik
     <int>     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>
1        1      0.05          0.05  0.99      415.       0     8 -75907.
2        2      0.87          0.87  0.37    43737.       0     9 -23024.
3        3      0.97          0.97  0.19   187918.       0     9  13378.
# ... with 4 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
#   df.residual <int>

様々な指標がまとめて出力されますが、今回は自由度調整済み決定係数(adj.r.squared)を見ることにしましょう。 透明度のみモデル(モデル1)と比べると、カラット数を加えることでモデルの説明力は大幅に上昇しています(モデル2)。 カラット数に対数をとるとさらに説明力が改善されており(モデル3)、3つのモデルの中ではモデル3がベストであると言えそうです。

サブサンプルに分けて分析する場合

先程は複数のモデルを同じデータに当てはめましたが、同じモデルを複数のデータにまとめてフィットさせることもできます。 たとえば透明度によって価格とカラット数の関係が異なるのかを調べてみることにしましょう。 透明度の低いダイヤは大きくてもあまり価値は上がらないが、透明度が高い場合はサイズが大きくなると価値が大きく上昇する、といった関係があるかもしれません。 これは透明度でデータをサブサンプルに分割して、各データに対して回帰モデルを当てはめることで確認できます。

df_nested = df %>% 
    group_by(clarity) %>% #nest(-clarity)でも同じ
    nest() %>% 
    arrange(clarity)

df_nested

# A tibble: 8 x 2
  clarity data                 
  <ord>   <list>               
1 I1      <tibble [741 × 9]>   
2 SI2     <tibble [9,194 × 9]> 
3 SI1     <tibble [13,065 × 9]>
4 VS2     <tibble [12,258 × 9]>
5 VS1     <tibble [8,171 × 9]> 
6 VVS2    <tibble [5,066 × 9]> 
7 VVS1    <tibble [3,655 × 9]> 
8 IF      <tibble [1,790 × 9]> 

nest()はグループごとにデータフレームを分割してdata列に格納する関数です。 data列にmap()を適用することで各サブサンプルに対して同じモデルをまとめてフィットさせることができます。

df_result2 = df_nested %>% 
    mutate(model = map(data, ~lm(log(price) ~ log(carat), data = .)),
           tidied = map(model, tidy))

モデルが作成できたので、回帰係数を見てみましょう。

df_result2 %>% 
    select(clarity, tidied) %>% 
    unnest() %>% 
    filter(term != "(Intercept)") %>% #切片は無視
    mutate_if(is.double, round, digits=2) 

# A tibble: 8 x 6
  clarity term       estimate std.error statistic p.value
  <ord>   <chr>         <dbl>     <dbl>     <dbl>   <dbl>
1 I1      log(carat)     1.53      0.02      84.4       0
2 SI2     log(carat)     1.79      0        557.        0
3 SI1     log(carat)     1.82      0        676.        0
4 VS2     log(carat)     1.78      0        563.        0
5 VS1     log(carat)     1.83      0        473.        0
6 VVS2    log(carat)     1.85      0.01     338.        0
7 VVS1    log(carat)     1.83      0.01     245.        0
8 IF      log(carat)     1.87      0.01     169.        0

I1だけは係数が小さくなっており、その他はほとんど同じ係数です。 どうやらと透明度が最低ランクの場合は、他のランクと比べてカラット数が価格に与える影響は小さいようです。

散布図も確認すると、確かにI1のみ傾きが少し緩やかになっているような気がします。

df %>% 
    sample_frac(0.3) %>% #データが多いので減らす
    ggplot(aes(log(carat), log(price))) +
    geom_point(color = "gray", alpha = 0.5) +
    facet_wrap(~clarity) + 
    geom_smooth(method = "lm", se = F, color = "#0072B2")  +
    theme_scatter

f:id:dropout009:20190102144550p:plain

まとめ

このように、探索的にデータを見ていく際には、purrrbroomを使うと複数の回帰モデルを効率的に比較することができ、とても便利です。