Gregorutti, B., Michel, B., Saint-Pierre, P. (2015) Grouped variable importance with random forests and application to multiple functional data analysis. Computational Statistics & Data Analysis, 90, 15-35.
仕事の都合で読んだ奴。ランダムフォレストとかで、個々の変数についてのpermutation重要性じゃなくて、分析者が定義したなにかしらの変数グループについてpermutation重要性を求めるという論文。
第一著者は、おそらくこの論文の提案手法を実装したであろうRパッケージRFgrooveを公開しているが、開発したきりメンテしておらず、CRANからは最近削除されている模様。
28頁あるけど本文は19頁だ、なんとかなるさ！と思って読み始めたけど…
1. イントロダクション
[メモ省略]
2. 変数グループの重要性指標
\(Y\)を確率変数、\(X^\top = (X_1, \ldots, X_p)\)を確率変数ベクトルとする。回帰関数を\(f(x) = E[Y|X=x]\)とする。
まずはBreiman(2011)が考えたpermutation重要性について。
\(X_j\)の変数重要性を以下とする。\(X_j\)だけを他のすべてから独立にreplicateした確率ベクトル\(X_{(j)}\)をつくり、$$ I(X_j) \equiv E[(Y – f(X_{(j)}))^2] – E[(Y – f(X))^2]$$ [えーと、つまり具体的には、検証用のデータ行列の\(X_j\)の列だけ縦にシャッフルするとMSEがどれだけ上がるかってことだと思う。そうですよね？]
さて、ご提案する変数グループ重要性は以下のとおりです。\(X_j\)から\(k\)個を選んだサブベクトル\(X_j = (X_{j_1}, X_{j_2}, \ldots, X_{j_k})^\top\)について、他のすべてから独立にreplicateして元の\(X_j\)に戻したやつを\(X_{(J)}\) として、$$ I(X_J) \equiv E[(Y – f(X_{(J)}))^2] – E[(Y – f(X))^2]$$ [はいはい。まあここまでは誰でも思いつく話だよな。他のすべてから独立にreplicateするというのは、\(k\)列についてそれぞれ縦にシャッフルするという話ではなくて、\(k\)列を取り出して同時に縦にシャッフルするという話であろう]
\(X\)のうち\(X_J\)に選ばれなかった残りの部分を\(X_{\bar{J}}\)とする。以下の加法回帰モデルを考えよう。$$ Y = f(X) + \epsilon = f_J(X_J) + f_{\bar{J}}(X_{\bar{J}}) + \epsilon$$ ただし\(E[\epsilon|X] = 0\)で、\(E[\epsilon^2|X]\)は有限とする。このとき以下が成り立ちます。証明は付録を見よ。$$ I(X_J) = 2 Var[f_J(X_J)]$$ [えーと、いま誰かがなにか回帰モデルみたいなのを作ったとしますね。それをよくみたら、\(k\)個の変数を使う前半部分と残りの変数を使う後半部分との足し算でできていた、としますね(ふつうの重回帰ならば常にそうだ)。前半部分が生み出す分散がありますやんか。それを二倍しますやんか。その値は「もとのモデルで\(Y\)を予測するとき、前半部分の\(k\)個の変数をこっそり縦にシャッフルしたらMSEがどんだけ上がっちゃうか」になりますよ。という話であろう。
そ、それはあれですよね、元のモデルのパラメータの値をMSEを最小化するようにして決めていればそうなる、ってことですよね？ もしパラメータの値をめっちゃ適当に決めてたら、モデルの前半部分は分散を生むけど、説明変数を縦にシャッフルしてもMSEは上がらないじゃないですか]
ここからさらに以下が成り立ちます。いま$$ f_J(x_J) = \sum_{j \in J} f_j(x_j)$$ であるとします。このとき、
- もし\( (X(j))_{j \in J}\)が独立ならば、$$ I(X_j) = 2 \sum_{j \in J} Var(f_j(X_j)) = \sum_{j \in J} I(X_j)$$
- もし任意の \(j \in J\)について\(f_j\)が\(f_j(x_j) = a_j x_j\)のような線形関数ならば、\(\alpha_J = (\alpha_j)_{j \in J}\)として $$ I(X_J) = 2 \alpha^\top_J Cov(X_j) \alpha_J$$
[よし、これも日常語に訳してみよう。えーと、さっきの前半部分が、さらによくみると、\(k\)個の変数それぞれを使う部分の足し算になっていたとしますね(ふつうの重回帰ならば常にそうだ)。
ひとつめ: もし\(k\)個の変数が互いに独立なら、\(k\)個の変数のグループとしての重要性、つまり「もとのモデルで\(Y\)を予測するときそれらを縦にシャッフルしたらMSEがどんだけ上がるか」は、\(k\)個それぞれの変数の重要性の合計だ。うん、これはわかるね。
ふたつめ: \(k\)個の変数が互いに独立だという話は忘れてくれ。かわりに、前半部分はまじでふつうの重回帰モデルだと思ってくれ。このときには、\(k\)個の変数の共分散行列を、前半部分の偏回帰係数ベクトルで挟んでスカラーに変えて二倍したのが、\(k\)個の変数のグループとしての重要性になる。うん、これもわかる。だってそれって前半部分が生む分散を二倍した奴だもんね]
ついでに、変数グループ重要性のリスケールしたバージョンをご提案します。$$ I_{nor}(X_J) \equiv \frac{1}{|J|} I(X_J)$$ [\(|J|\)ってのはつまり\(k\)だと思うんだけど、なぜそう書かないんだろう。ほんとにもうね、数学ができる人って性格悪いよね]
ここからはランダムフォレストの話。
おさらいします。\( (X, Y) \)の\(n\)個のiidなreplicateがあるとして、これを訓練データ\(D_n\)としましょう。ここから\(M\)個のブートストラップ標本\(D^1_n, \ldots, D^M_n\)をつくり、それぞれに基づいて回帰木による推定量 \(\hat{f}_1, \ldots, \hat{f}_M\)をつくります(ふうつのCARTと違って、各ノードの分割ルールに使える変数セットはランダムに決めます)。でもって、これを平均します。
\(m\)番目のブートストラップ標本について、それに対応するOOB標本を \( \bar{D}^m_n \equiv D_n \backslash D^m_n\)とします。\(\hat{f}_m\)のリスクのOOB推定量は$$ \hat{R}(\hat{f}_m, \bar{D}^m_n) = \frac{1}{|\bar{D}^m_n|} \sum_{i:(X_i, Y_i) \in \bar{D}^m_n} (Y_i – \hat{f}_m (X_i))^2$$ \(X_j\)のpermutation重要性指標は次の通り。それぞれのOOB標本\(\bar{D}^m_n\)について\(X_k\)をランダムにpermuteしたのを\(\bar{D}^{mj}_n\)として、$$ \hat{I}(X_j) = \frac{1}{M} \sum_{m=1}^M \left[ \hat{R}(\hat{m}, \bar{D}^{mj}_n) – \hat{R}(\hat{m}, \bar{D}^{m}_n) \right]$$
では、変数グループのpermutation重要性のご提案。
それぞれの\( \bar{D}^m_n\)において変数グループ\(X_J\)をランダムにpermuteしたのを\(\bar{D}^{mJ}_n\)とします。つまり、\(X_j\)の経験的な同時分布は変えず、\(X_J\)とその他の変数とのつながりだけ切るわけです。でもって、$$ \hat{I}(X_J) = \frac{1}{M} \sum_{m=1}^M \left[ \hat{R}(\hat{m}, \bar{D}^{mJ}_n) – \hat{R}(\hat{m}, \bar{D}^{m}_n) \right]$$ [はいはい。これも誰でも思いつく話である]
3. 変数グループの重要性を用いたmultiple functional data analysis
[この節から関数データ分析の話になる。あれですよね、なんかこう時系列とかあるときに、それにむりくり関数をあてはめちゃって、そのパラメータをデータとみよう、というような話ですよね… 前にセミナーで聴いた程度で、ほぼ未知の世界だ。のみならず、冒頭からヒルベルト空間とかDaubechiesウェーブレット基底とかが乱舞するいかつい世界で、１パラグラフも読まずに断念した]
4. 事例: aviation safetyの変数選択
[なんかこうソコハカとなく面倒になってきちゃったので、この章も読まずにとばした]
5. 結論
[メモ省略]
———–
ちゃんと読んでないのにアレですけど、変数グループのpermutation重要性の定義は、うん、まあそういう話だろうな、という内容であった。要は、注目している変数グループの列をまとめて縦方向にシャッフルしてモデルに放り込んだとき、MSEだかなんだかがどう変わるかみてみましょうという話だと思います。
つまらない感想だが、ランダムフォレストにおけるOOBベースのpermutation重要性の算出はもうパッケージに実装されちゃっているから、変数グループについてpermutationするのはかえって面倒くさそうだなあ。RのrandomForest::randomForest() なり party::cforest() なりが返すオブジェクトから、個々の木にとって訓練データのどの行がOOBかを調べるのってどうやんのさ？ party::varimp()のコードを読めばわかるんだろうけど… もう面倒くさいからk-folds CVでやり直したほうが早いような気がしてきた。