ランダムフォレスト特徴量の重要度（C++の実装例つき）

はじめに

今回の記事は、alfredplplさんの以下の記事とだいぶかぶっています…。図つきで、とても分かりやすい記事なので、お勧めです。こちらをはじめに読んだほうが良いと思います。<a href="http://alfredplpl.hatenablog.com/entry/2013/12/24/225420">Random Forestで計算できる特徴量の重要度 - なにメモ（なにかしらのメモ帳）</a>alfredplpl.hatenablog.com
alfredplplさんの記事は分類を例にあげてましたので、こちらは回帰を例にあげて説明します。基本的な中身は一緒ですので、比較しながら読んでもらえればと思います。

Rのランダムフォレストの重要度

Rのランダムフォレストには２種類の重要度(importance)があります。この記事で説明するのは、IncMSEとIncNodePurityです。

アプローチ	回帰(regression)	分類(classification)
シャッフルしたOOBデータで予測して求める	IncMSE	MeanDecreaseAccuracy
決定木のノードの純度(Purity)の増分から求める	IncNodePurity	MeanDecreaseGini

Rでの求め方

(以下の例はTopCoder機械学習マッチ"ChildStuntedness5"でのデータを使用しました。要ログイン→TopCoder)
(1)randomForestの引数にimportance=TRUEをセットする。（注：以下はただの例ですので、他の引数とかは適宜変更してください。）

forest <- randomForest(GENIQ~., data=data.train, ntree=200, importance=TRUE)

(2) プロットしたいならvarImpPlot、しないならimportanceを呼ぶ。

ret1 = varImpPlot(forest)
ret2 = importance(forest)

結果は以下のようになります。

f:id:shindannin:20150424021604p:plain

IncMSEとIncNodePurityは別なので、重要度の値はもちろんのこと、上記のように順位が異なってくる場合もあります。
上記の方法ではなく、importance(forest)で重要度を出力すると、IncNodePurityは標準誤差で割られた値になります。*1

IncMSE

C++での実装例

C++の実装例と実行結果をideoneにアップしました。#define INCMSE (1)で使用するしないを切り替えられます。#if INCMSEで囲まれている部分を見れば、どう計算してるか分かると思います。もし、ランダムフォレストそのものの部分が分からない場合は、ランダムフォレストのつくりかた（C++の実装例つき） - じじいのプログラミングを参照してください。

アルゴリズム

入力変数X, 目標変数をyとします。データサンプルは(X1,y1),(X2,y2),...(Xi,yi),...(Xn,yn)のn個あるとします。
またXはベクトルでXi = {xi1, xi2, xi3, ..., xim}とm個の特徴量を持つとします。

以下の処理は、ランダムフォレストの決定木１本１本に対して行います。

ランダムフォレストで、決定木を作るときに訓練に使われなかったデータサンプル（ブートストラップサンプリングで外れたデータ）を使います。これをOOBデータと呼びます（Out-of-bag(OOB) data、OOB examplesとも）。
まずOOBエラーを求めます。OOBデータのXを決定木に入力すると、予測した結果y'を求められます。y'と目標変数yがどれだけズレてるかがOOBエラーです。回帰のときは、MSE(mean square error、平均二乗誤差)をOOBエラーとします。つまりΣ(yi-yi')^2 / kとなります（kはOOBデータのデータ個数）。
m個の特徴量それぞれに対して重要度を求めます。単純にm回ループして、ループ変数を特徴量fとして、以下のことをします。
1. OOBデータの特徴量fだけをシャッフルします。他の特徴量とyはシャッフルしないので、気を付けてください。
  - X1=( 1.0, 2.1, 3.2) y1=1.2
  - X2=( 0.5, 1.6,-2.1) y2=2.1
  - X3=(-1.2, 1.5, 2.0) y3=3.5
  - ↓シャッフル
  - X1=( 1.0, 1.6, 3.2) y1=1.2
  - X2=( 0.5, 1.5,-2.1) y2=2.1
  - X3=(-1.2, 2.1, 2.0) y3=3.5
2. シャッフルしたデータを使って、(2)と同じようにOOBエラーを求めます。特徴量fだけシャッフルしたので、特徴量fが重要なのであれば、重要なデータがシャッフルでおかしくなったわけなので予測が悪くなるはずです。予測が悪化した値、つまり、(シャッフルしたデータのOOBエラー) - (元のOOBエラー) を重要度IncMSEとします。

最後に、全ての決定木のIncMSEの総和を求めて、最後に木の本数で割って平均を求めます。

ポイント

計算時間は、かかります。訓練の計算は1回だけで済みますが、予測は、（OOBデータ個数）＊（木の本数）＊（特徴量の数）＊（特徴量1つあたりのシャッフル予測の回数 nPerm）回だけ行うことになります。予測は訓練に比べたら計算時間はずっと短いですが、回数が増えればもちろん時間はかかります。時間がかかりすぎの場合は、上記の4つのどれかを減らす必要があります。
重要度として、それなりに信頼できます。シャッフルする部分以外は、通常のOOBデータを使った予測と変わらないので。
マイナスの値もとりえます。予測に影響をほぼ及ぼさない特徴量だと、シャッフルして、たまたま逆に予測精度があがることもあり得ます。もちろん、重要度は最悪という意味です。
Rでの変数名は、%IncMSEと%がついていますが、全部足して100%になるわけではありません。

R Projectによる、C言語での実装例

本体はC言語で書かれています。CRAN - Package randomForestからrandomForest_4.6-10.tar.gzをダウンロードしましょう。 src/regrf.c の262-295行と327-338行のあたり、if(varimp){} のスコープ内のコードを読むと分かりやすいと思います。

IncNodePurity

C++での実装例

先程のコードの#if INCNODEPURITYで囲まれている部分で、求めています。

アルゴリズム

ランダムフォレストのノードの分け方を知ってるなら、何も新しい知識はいりません。
決定木を作る際に、各ノードを子2つに分けるベストな分け方は、（親のMSE）-（子2つのMSE）が最大をとるときです（（親のMSE）は固定なので、（子2つのMSE）の最小化でもいいです）。数式で書けば以下のようになります。

Aを集合として、i番目の目的変数をYiとして、親の集合Pを、子の集合Lと集合Rに2つに分割する。
要素数 |A|
総和 S[A] = {Σ[i∈A]Yi}
平均 E[A] = S[A]/|A|
と表す。MSE（平均2乗誤差）の変化量は以下のようになる。
  {Σ[i∈P](Yi-E[P])^2} - {Σ[i∈L](Yi-E[L])^2} - {Σ[i∈R](Yi-E[R])^2}

これを選ばれた特徴量別に足していくだけです。最後に木の数で割ります。

ポイント

計算時間は速いです。決定木を作る際に求まるので、訓練をすればいいだけで、予測の時間はかかりません。これはIncNodePurityの最大のメリットです。
重要度としては、あまり信頼できません。確かに分割の際に使われる特徴量は重要とは言えますが、あまり直接的ではないので。ただ、例えば、特徴量が1000個あるときに下位50個を削りたいと言った、ざっくりとした特徴量抽出には使いやすいです。（実際、1位を取ったマッチToxCast Challenge | EPA TopCoderでは、これを使いました。）
決定木の数が増えても、値は増えません（木の数で割って平均をとってるので）。
サンプルデータの数に比例して、値は大きくなります。
経験上、決定木の数が少ない場合でも、かなり値は安定することが多いです。最初は多くの決定木を使って、次に少ない決定木を使って、結果があまり変わらなかったら、それ以降は少ない本数で試すというアプローチが、時間を節約にもなって良いと思います。

R Projectによる、C言語での実装例

先ほどダウンロードした、randomForest_4.6-10.tar.gzを参照してください。

src/regrf.c の変数tginiがIncNodePurityになります。
src/regTree.cで、tginiを計算していますが、変数crit→critvar→critmax→decsplit→tginiというふうにデータが渡されています。247行目 crit = (suml * suml / npopl) + (sumr * sumr / npopr) - critParent; が重要です。平均2乗誤差を使う場合は式展開すると、以下のように簡単になり、これで合ってます。


  {Σ[i∈P](Yi-E[P])^2}                          - {Σ[i∈L](Yi-E[L])^2}                          - {Σ[i∈R](Yi-E[R])^2}
= {Σ[i∈P]Yi^2} - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2 - {Σ[i∈L]Yi^2} + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2 - {Σ[i∈R]Yi^2} + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... 2乗を展開した
=                - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2                  + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2                  + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... {Σ[i∈P]Yi^2} = {Σ[i∈L]Yi^2} + {Σ[i∈R]Yi^2} で消した。
=                - 2E[P]S[P]         + |P|E[P]^2                  + 2E[L]S[L]         - |L|E[L]^2                  + 2E[R]S[R]         - |R|E[R]^2　... 総和Sであらわした
=                - 2S[P]^2/|P|       + S[P]^2/|P|                 + 2S[L]^2/|L|       - S[L]^2/|L|                 + 2S[R]^2/|R|       - S[R]^2/|R| ... E[A] = S[A]/|A|で平均Eを消した
=                -  S[P]^2/|P|                                    +  S[L]^2/|L|                                    +  S[R]^2/|R|                    ... 足し算

重要度をどう使うか？

特徴量選択

無意味な特徴量を削ることで、予測精度を上げることができます。また、計算時間も速くなり、多くの木を作れるので、それでも精度を上げることができます。

以下に、自分が行う特徴量選択の方法をまとめてみました。今回の2つは(3)(4)にあたります。

手法	選択性能	計算時間
(1) 特徴量を削りクロスバリデーションで予測するのを、全ての特徴量で試す	良い*2	とても遅い
(2) 特徴量を削りOOBデータを使って予測するのを、全ての特徴量で試す	良い	遅い
(3) IncMSEで重要度の低い特徴量を削ってみる	普通	普通
(4) IncNodePurityで重要度の低い特徴量を削ってみる	悪い	短い

かけられる時間やデータ量にもよるのですが、(4)→(3)→(2)→(1)の順に試すのが良いと、自分は思います。基本は(2)になると思います
(1)(2)(3)については、MSEではなく、精度の評価関数（scoring function）に置き換えることもできます。例えば、yの大きい順を予測したい順位当ての問題のときは、順位差の絶対値の総和をMSEの変わりに使うといったこともできるでしょう。
ただ「重要度を低い特徴量を削る」と簡単に書きましたが、貪欲に１番重要度が低い特徴量を削っていけば良いというものでもなく、とても難しいです*3。特徴量間に強い関係があれば、特徴例をけずるたびに、重要度も大きく変わってきます。良い方法があったら（というか研究されていると思いますが）、ぜひ教えてもらえるとうれしいです。一応、自分のアイディアレベルのネタとしては、
- (2)(3)であれば、複数の特徴量をシャッフルする手はあるでしょう。ただ、（OOBデータ個数）＊（木の本数）＊（特徴量の数 ^ シャッフル個数）＊（特徴量1つあたりのシャッフル予測の回数 nPerm）回となるので、特徴量が多いときは、シャッフル個数2個にするのすら遅すぎて難しいかもしれません。
- 探索枝仮り・焼きなまし法・ビームサーチといった最適化手法で削っていく方法もあるでしょう。
- 特徴量を１つ削ったとき、他の重要度がどう変化するかで、どの特徴量間に強い関係があるかというグラフのようなものを作れば、最適化手法を使っていく上で役に立つかもしれません。例えば、「特徴量はx1,x2は片方を削ると、もう片方の重要度が大きく下がるので、たぶんペアで意味がある特徴量だ。焼きなまし法を試すときは、この２つはまとめて追加したり削除したりするようにしよう。」といったことはできそうな気がします。

ランダムフォレストのパラメータ調整がひどいのに気付くきっかけ

最初にあげた図では、SUBJIDという特徴量は、重要度IncMSEで最小、重要度IncNodePurityで最大になっています（SUBJIDは、本当に無意味なIDで、実際はIncMSEが正しかったです）。
このように、IncMSEの大きさの順位とIncNodePurityの大きさの順位が異なる場合、経験上ですが、予測精度も悪くなる傾向にあります（）。恐らく原因としては、

元の特徴量が少なすぎる
木を深いところまで作りすぎている（nodeSizeが大きい or maxLevelが大きい）

のあたりでしょう。
ideoneのコードの645行目

	rf->train(trainingFeatures, trainingAnswers, 5, 999, 2, 5, 0.66f, numTrees);

のnodeSizeとmaxLevelのパラメータを変化させ、IncMSE・IncNodePurity・totalErrorの変化をみてみると良いと思います。

以上です。ご感想・ご意見お待ちしてます。

*1:今回は詳細は述べませんが、詳しく知りたい方は、RのソースコードをimportanceSD, impSDで検索したり、importance.Rからソースを追ってみてください。SDという名前で標準偏差に見えますが、なぜか標準誤差なので注意してください（ネットでも多数の間違いが見られます…）

*2:「とても良い」の可能性もありますが、OOBデータ vs クロスバリデーションの議論は奥深いようなので踏み込みません

*3:あくまで自分のイメージですが、「特徴量の追加削除」を近傍、「重要度」を関数の傾きと置き換えると、近傍にいっただけで関数の傾きが大きく変化したりする関数を最適化しないといけないようなものなので、やばいかんじしかしません…