【R】{randomForestSRC}パッケージの基本事項備忘用メモ

なんか日本語情報が全然なかったので、自分用メモ
別にrandom forestの計算ロジック自体が変わるわけではないので*1、情報がなくても不思議じゃないんだけど。

ちょいちょい触ってみて、{randomForest}パッケージより使いやすいと思ったので必要最小限に書き残す

参考URLなど
普通に使う
パラメータチューニング
分析後の可視化いろいろ
交互作用をみつける
感想

参考URLなど

基本的には公式チートシートみればよくて、それが分からなかった関数一覧のページ見ればだいたい事足りる（んだけど一応自分用にメモする）

普通に使う

とりあえず著名なデータセットであるボストン住宅価格データとアヤメ(iris)データ*3を使う

data(Boston , package = "MASS")
data( iris)

実行のコードはあまり{randomForest}とかわらない。
だいたい引数でいじるのは、ntree（生成する木の数）, mtry（選ぶ変数の数。回帰問題ならデフォ値 $p/3$ でそれ以外はデフォ値 $\sqrt{p}$ ）, nodesize（ノードに含まれるべき最小ケース数）あたりくらい。

分類問題

rf_cls0 <- rfsrc( Species~ . , data = iris ,  nodesize =  4 )
print( rf_cls0 )

（summary()ではなく）print()したときに割といい感じの欲しい情報をまとめて出してくれるので、おぼえとくとよい

回帰問題

rf_reg0 <- rfsrc( medv ~ . , data =  Boston, nodesize =5 )
print( rf_reg0)

生存問題

（いつかつかうときが来たらここに追記する）

パラメータチューニング

randomForestSRC::tune()を使うと、mtryとntreeに関して最適な値をみつけてきてくれる

## 回帰問題 ##
# mtryとnodesizeについて最適化
ptune_reg0 <- randomForestSRC::tune( medv~. , data = Boston)
# nodesizeを固定にする
ptune_reg0_2 <- randomForestSRC::tune( medv~. , data = Boston,nodesizeTry = 5) 

#結果
ptune_reg0$optimal
# nodesize     mtry 
# 1            　 13 

## 分類問題 ##
ptune_cls0 <- randomForestSRC::tune( Species~.,data =iris)
#結果
ptune_cls0$optimal
#　nodesize     mtry 
#    3      　　 3

分析後の可視化いろいろ

予測精度を評価する

MSEとか一致率を求める
$predictedでIn-Bagに対する予測値を、$predicted.oob はOut-Of-Bagに対する予測値を取り出せる

ただ計算するだけなのもアレなので、伝統的手法とくらべてみる

library( tidyverse)

##分類：Multinomial Logitと比較
library(mlogit)
iris_md <-  mlogit.data(  shape = "wide" , data = iris , choice = "Species" )
iris_ml0 <- mlogit( Species~ 0 | Sepal.Length + Sepal.Width+ Petal.Length+ Petal.Width ,data= iris_md)
pred_lbl_ml <- iris_ml0$probabilities %>% 
  apply( 1 , function(x){
    return( which( x == max(x)))
  }) %>% unlist


pred_lbl <- rf_cls0$predicted.oob %>% 
  apply( 1 , function(x){
    return( which( x == max(x)))
  }) %>% unlist

err_rate_ml <-  1 -  ((as.numeric( iris$Species) ==pred_lbl_ml ) %>% mean) # 0.01333333 (Error rate)
err_rate_rf <- 1 -  ((as.numeric( iris$Species) ==pred_lbl ) %>% mean) #0.04 ( Error rate)

##回帰：普通にOLSする場合との比較
ols0_bos <- lm( medv~ . , data = Boston) # R^2 = 0.7406
mse_ols0 <-  mean( (Boston$medv-  ols0_bos$fitted.values ) **2) #21.89483(MSE)
mse_rf0 <-  mean( (Boston$medv -  rf_reg0$predicted.oob) **2) #11.22804(MSE)

まぁ、造作もないですが。
なんかirisデータの分類に関しては普通に多項ロジットやったほうが正確だ....

変数重要度

○回帰問題

vimp_reg0 <-  vimp(rf_reg0)
plot( vimp_reg0)

おなじみの図

○分類問題

vimp_cls0 <- vimp( rf_cls0)
plot( vimp_cls0)
vimp_cls0$importance #重要度一覧

おなじみの図,part2

Partial Dependence Plot

各変数に対するpartial dependence plot（わかりやすい解説は「変数重要度とPartial Dependence Plotで機械学習モデルを解釈する - Dropout」とか参照）は、plot.variable()で引数partial=Tを指定することで得られる

plot.variable( rf_reg0,partial = T) #回帰
plot.variable( rf_cls0,partial = T) #分類

交互作用をみつける

結構便利な機能。アウトカムYに対して、変数 $X_1, X_2, \ldots, X_J$ の中から交互作用の大きい組み合わせをみつける手助けをする
method引数に"vimp"（変数重要度の増減への寄与基準）か”maxsubtree"（分割に各変数が現れる深さ基準）かを選べる。
vimpのほうが直感的で好みである。その説明（URL）をみよう

method="vimp"

This invokes a joint-VIMP approach.
Two variables are paired and their paired VIMP calculated (refered to as 'Paired' importance). 
The VIMP for each separate variable is also calculated. The sum of these two values is refered to as 'Additive' importance.
A large positive or negative difference between 'Paired' and 'Additive' indicates 
an association worth pursuing if the univariate VIMP for each of the paired-variables is reasonably large.
See Ishwaran (2007) for more details.

とある。まぁ詳細は元論文をよまないと分からんけど、
二変数の組み合わせから算出できる"paired importance"と単なる二変数のimportanceの和（＝"additive" importance）の値を比べて、そこに差があってなおかつ元の各変数のimportanceも高ければ、そこには分析しがいのある交互作用があるよ、ということらしい。
※プラスはなんとなくわかるけど、マイナスってどういうこと？ってのがよくわからない。。
※paired importance基本的に総当たりで計算しているので、変数の数が多すぎるときは引数nvarで絞ると良い

とりあえずやってみる

### 交互作用をみつける ###
int_reg0_1 <-  find.interaction(rf_reg0,method= "vimp")
int_reg0_1 %>% head()

以下のようなデータフレームがアウトプットとなっている

Differenceの絶対値でソートして、top5をみてみる

dfInt_reg0 <- 
int_reg0_1 %>% as.data.frame() %>% 
  mutate( abs_diff = abs(Difference)) %>% 
  arrange( desc( abs_diff))
head( dfInt_reg0,5)

difference上位の交互作用項がほんとに有意を、OLSであらためてたしかめてみる

ols1_bos <- lm( medv ~. + lstat:rm + rm:ptratio, data = Boston)
summary( ols1_bos)

OLSでも交互作用項が有意になっていることが確認できる

感想

わりと使いやすそうだし、何よりわかりやすめなので、今後random forestを使うときは{randomForestSRC}も良い選択肢だな、と思いました。
PDPの図示とか交互作用の探索とかもパッとできるので、プレ分析用のツールとしてはよい気がする

www.youtube.com

Enjoy!!

*1:あと、今はRよりPython使うひとのほうが多いだろうし

*2:へぇ..そんなのあるんだ感。普通に二回やればいいんじゃないかって思うけど、Y1とY2の両方を説明するXをみつける、みたいなのが求められる問題も世界のどこかにはあるんだろう。時間あるときに数式のとこ読んでみたい

*3:なんかこのデータを使うこと倫理的妥当性について喧々諤々が最近あった気がするけど、気にしない

論理の流刑地

地獄の底を、爆笑しながら闊歩する