Debeer, D., Stroble, C. (2020) Conditional permulation importance revisited. BMC Bioinformatics, 21:307.

　仕事の都合で読んだ。ランダム・フォレストにおいて変数重要性を評価する手法のひとつに、「他の変数で条件づけたパーミュテーション重要性」というのがあるんだけど、その算出方法についての論文。
　Rにおける既存の実装についていくつか疑問点があったのでなんとなく目を通して見たら、なんと著者らはRパッケージpartyの中の人であり(途中で気が付いた)、しかも、途中でなんだかとんでもないことを言い出した…

背景
　ランダムフォレスト(RF)の提案者Breimanさんたちは、RFにおける変数重要性指標として、Mean Decrease in Impurity (MDI)とMean Decrease in Accuracyの二種類を考えていた。後者を本論文はpermutation重要性(PI)と呼ぶ。
　その後いろんな提案があった。大きく二種類に分けられる。

• 木の構造に基づく指標。Epifani(2017, BMC Bioinformatics)のIntervention in Prediction Measure, Ishwaranら(2010, JASA)のMinimal Depthなど。
• 関心ある予測子にノイズをいれることによる予測の正確性の変化に基づく指標。BreimanのPIがこれ。ほかにIshwaranら(2007, ElectronJStat), Stroble et al(2008, BMC Bioinformatics)[←Rのparty::varimp(condtional=T)ね] など。

周辺的重要性と部分的重要性
　そもそも重要性指標というのは、正確には何を測ってんのかいまいちはっきりしない。これはRFベースの重要性指標に限ったことではなく、そもそも線形回帰においてさえそうである。Hoffman(1960, Psych.Bull.), Pratt(1987), Thomas et al.(1998 Soc.Indic.Res.)とか、彼らを批判するBring(1996, Am.Stat.), Darlington(1968, Psych.Bull.), Ward(1962, Psych.Bull) とか読め。というかGromping(2015, Wiley Interdiscip Rev Comput Stat)を読め。
　読んでないお前らのために要約してやると[←意訳]、ふたつの極端な考え方があるわけよ。

• 右端はmarginal重要性。他の予測子を考慮しない予測子のインパクトとして解釈できる。線形回帰でいえば相関係数である。
• 左端はpartial重要性ないしconditional重要性。他のすべての予測子を考慮した予測子のインパクトとして解釈できる。線形回帰でいえばたとえば半偏相関係数がそれだ[部分相関係数のことね。\(y, x_1, x_2\)の3変数ならば \(r[y, (x_1|x_2)]\)のこと]。

すべての重要性指標はこの軸上のどこかにいる。重要性指標にとってどこが理想点なのかについて合意はない。これはそれぞれが自分のリサーチ・クエスチョンに適した場所を選んでいるということであり、合意がないからといって重要性指標自体がダメということにはならない。Breimanさんも講演で両方の視点が大事だと云っている[←へー]。

RFにおける変数重要性
　線形回帰とちがって、RFにはアウトカムについての統計的モデルがないし、重要性指標を真のデータ生成メカニズムの特性と捉えることもできない。単に、フィットしたRFにおける予測子の貢献を示すのみである。
　そういうわけで、RFにおけるmarginal重要性とpartial重要性とはなんなのかも、実のところいまいちよくわからんのだが、Stroble et al.(2008)はこの概念をRFに持ち込み、オリジナルのPIはmarginal重要性なのでpartial重要性の指標(CPIと呼ぶ)も作ろうと主張したわけである。いまではCPIは広く使われている。[←へー。まあ俺もR Journalで読んでかなり説得されたよね]

　本研究はCPIの改善を提案する。実は最初はRのpartyパッケージを改善したりしてたんだけど、もっといろいろ思いついちゃったもんで、新たにpermimpパッケージというのを作りました。
　[あっ！この論文の第二著者はStroblさんその人だ！ いま気が付いた…]

オリジナルのpermutation重要性
　背後にある考え方はこうである。\(Y\)と\(X_k\)がなんらかの依存的な構造をもっているとき、\(X_k\)にノイズを足してやれば、その構造は壊れるはずである。そうすることでRFによる予測の正確性が下がるなら、\(X_k\)は予測において重要だったのだということになる。
　木\(t\)におけるOOB標本\(\beta^{(t)}\)に基づく予測誤差を\(R^{(t)}\)とする。分類木なら$$ R^{(t)} = \sum_{i \in \beta^{(t)}} \frac{I(\hat{y}^{(t)}_i \neq y_i)}{|\beta^{(t)}|}$$ 回帰木なら $$ R^{(t)} = \sum_{i \in \beta^{(t)}} \frac{(\hat{y}^{(t)}_i – y_i)^2}{|\beta^{(t)}|}$$としておく。で、OOB標本において\(X_k\)の値をpermuteしたあとの予測誤差を\(R^{(t)}_{(k)}\)とする。木ごとのPIを$$ PI^{(t)}_{(k)} = R^{(t)}_{(k)} – R^{(t)}$$ とし、森のPIは木のPIの平均とする。
　さて、\(X_k\)以外の予測子を\(Z_{(-k)}\)と書くとして、PIは\(X_k\)と\(Y, Z_{(-k)}\)が独立な時に0になるが、大きな値になったとしてそれがどっちの依存性のせいかはわからない。実際、すべての\(X_k\)がそれぞれ\(Y\)と独立であっても、\(X_k\)と\(Z_{(-k)}\)が独立でないときにはPIは大きな正の値になりうる。[ちょ、ちょっと待って、それってどういう状況よ… Nicodemus, et al.(2010, BMC Bioinformatics)というのが引用されている]

条件つきpermutation重要性
　そうではなくて、\(Z_{(-k)}\)で条件付けたもとで\(X_k\)と\(Y\)が独立である時に0となる指標をつくろう。というのがCPIである。具体的には、予測子の空間を\(Z_{(-k)}\)でパーティションに区切って、各パーティションのなかで\(X_k\)の値をpermuteする。
　ここで問題になるのは、筋が通っていてかつ計算的にみても実現可能なパーティションをどう決めるかである。Stroblらの提案は、木の成長アルゴリズムによって作られたパーティションを使う、である。それならもう作られてるし、OOBには依存してないし、カテゴリカルでない予測子についても分け方が明確である。
　party::varimp(object, condtional = TRUE)での実装は次の通り。

1. それぞれの\(X_k\)について、\(Z_{(-k)}\)に含めるべき予測子を決める。それは\(X_k\)と連関している変数に限られるべきである。なので、個々の予測子について\(X_k\)との連関性の検定のp値を求め、1-pがユーザが定義した閾値を超えている奴を選ぶ。選ばれた予測子のセットを\(Z^{(s)}_{(-k)}\)とする。
   [覚え書きとしてメモしておくと、これはparty:::create_cond_list()でやっている。目的変数を\(X_k\)としたctree()を改めてコールしているようだ。ここでの予測子の分岐点の定義は元の木とは異なることになるわけで、コードを読んでてちょっとびっくりした点なんだけど、まあ考え方としては理解できる]
2. それぞれの木\(t\)における\(Z^{(s)}_{(-k)}\)のすべての分岐点を結合したグリッドを作り、予測子の空間を分割する。このグリッドにおいて予測子の空間はそれぞれの分岐によって完全に二分される、というところがポイント(木による分割の場合は、ある分岐点は前の分岐点によって限定された下位空間を分割するだけである)。で、分割された空間のなかで\(X_k\)のOOB値をpermuteし、\(R^{(t)}_{(k)}\)を求める。
   [そうそう、ここもvarimp()を使っててびびった点であった。たとえば\(X_k\)が\(X_1\)で\(Z_{(-k)}\)が\(X_2\)だけだとして、どっちも量的変数で、元の木ではまず\(X_1\)の高低で分割が生じ、高群では\(X_2\)が3以上かどうかで分割が起き低群では\(X_2\)が2以上かどうかで分割が起きているとする。結局「2以上3未満」という分割は元の木には起きていない。しかしvarimp()はOOB標本を、\(X_2\)が2未満、2以上3未満、3以上の3つに分割することになる。\(Z_{(-k)}\)に複数の変数があったらそのすべての組み合わせに分割するわけだから、結果として元の木からみるとありえないくらいの滅茶苦茶たくさんの分割が生まれ、その分割の中にはデータ点がひとつしかないことも多く、結果として\(X_k\)をpermutationしてもOOB予測誤差はそんなに大きくならないのだ。直感として、それってどうなの、って思うんですよね。実際の木においてはそこまで厳しく条件づけられてはいないわけだから]

4つの問題
さて、CPIが広く用いられるようになった結果、次の4つの問題があきらかになった。

1. 計算速度。party::varimp()は遅いという評判である。よく考えてみると実装に結構な無駄がある。たとえばあるパーティションのケースがすべて同じ終端に落ちる場合、そのパーティションでpermutationする意味はない。というわけで、partyパッケージを改善してかなり速くなり、さらにpermimpパッケージではもっと速くなっている。
   [いやいやいや先生… permimpパッケージは知らないが、最新のparty::varimp()についていえば、そもそもRコードとしてかなり下手な書き方してますよ… forループで行列要素をいじるのをやめてベクトル処理を徹底するだけで劇的に速くなりましたけど…?]
2. 線形の関連性という制限。\(Z^{(s)}_{(-k)}\)を決める際、たとえば両方が連続変数だったら線形の関連性しかみていない。permimpパッケージでは非線形な関連性があるときも選ぶようにした。
   [なるほどね？ でもそれってvarimp()というよりctree()の限界ってことですよね？]
3. 標本サイズへの依存。\(Z_{(-k)}\)を決める際、変数の型によって異なる連関指標を用いることになり、効果量を比較できないので、仕方なくp値を使ったわけだけど、(1)考えてみると木における分岐点のどっち側に落ちるかだけが問題なのであってローデータなんてどうでもいいはずだし、木はIn-bagデータを使って成長させるのに検定には全データを使っているのもおかしい。(2)p値は標本サイズに依存する。木を成長させていく上ではおかしくないが、\(Z_{(-k)}\)を決める際には、標本サイズが大きいほど貪欲な変数選択が起きることになっておかしい。
   [あっそうだ！これは気が付かなかった。permutationの条件付けの際にやたらたくさんの共変量が選択される理由がわかったよ…]
4. 不安定性。シミュレーションしてみるとparty::varimp()で求めたCPIは標本抽出に伴ってものすごく変動する。おそらくは\(Z_{(-k)}\)の選択が不安定だからだろう。木ごとじゃなくて森全体について決めちゃってるのがいけなかったね。
   [ちょっとぉーー！！！ いまさらそれはないんじゃないの？！ えええええ？！
   CPIの不安定性については前に指摘している人がいたけど、あなたたち全力で反論してたじゃないの… うわぁ、ひっでえなあ… ]

permimpの実装
　そんなわけで、permimpパッケージにおける実装についてご紹介する次第です。
　大きな違いは\(Z^{(s)}_{(-k)}\)の決め方。

• 森について決めるんじゃなくて、木ごとに決めます。\(Z^{(t)(s)}_{(-k)}\)と書くことにします。さらにいうと、その木のIn-Bagデータだけ使って決めます。また、閾値はもっと厳しくします。
• ローデータを使って決めるのをやめます。個々の予測子についてその木の分割点で分割してからカイ二乗検定で決めます。非線形性も捉えることができるし、標本サイズへの依存性も減ります。

[うーん、条件づける共変量を木ごとに決めるというのは納得である。いっぽう、共変量を問答無用でカテゴリ化してから共変量選択するというのは、よく考えてみると、それはそれでちょっと変な気がしますね。いま、木というブラックボックスの挙動を調べたい人がいる。そのために、予測子のなかに\(X_1\)と\(X_2, X_3, \ldots\)があるときに\((X_1|X_2, X_3, \ldots)\)が予測に貢献しているかをシミュレーションで調べたい。そこで、\(X_2, X_3, \ldots\)のなかからシミュレーションにおいて条件づけるべき変数を選択したい。よし、\(X_1\)と\(X_2\)との関係、\(X_1\)と\(X_3\)との関係… を調べていき、強い奴だけ条件づけてシミュレーションしよう。と思っていたら、そこに密告者が現れて、兄貴、あの木は\(X_1, X_2, X_3, \ldots\)をこんな分割点でカテゴリ化した\(X’_1, X’_2, X’_3, \ldots\)を使ってますぜ、とタレコミしてくれた。オーケー、そんなら\(X’_1\)と\(X’_2\)との関係、\(X’_1\)と\(X’_3\)との関係… を調べてけばいっか、と割り切れるのか？ という話である。実際には、\(X’_1\)と\(X’_2\)の間に関係がなくても\(X_1\)と\(X_2\)の間には関係があり得るし、本当に知りたいのは\((X_1|X_2, X_3, \ldots)\)の貢献なんだから、タレコミは無視すべきなのではないかという気がするんだけど？ うーん、まあいいか。共変量選択という問題に正解はないのかもしれない。
　それよりも、予測子の空間を\(Z^{(t)(s)}_{(-k)}\)のすべての組み合わせで細分化してpermutationするという方針は変わらないってことね。ほんとはそこに違和感を感じるんだけどな…]

　なお、permimpパッケージはpartyパッケージの森(ctreeに基づく)だけでなく、randomForestパッケージの森(impurity reduction原理に基づく)にも対応するようにしました。

閾値の解釈
　permimpにおける閾値は、CPIがmarginal-partial軸状のどこに位置するかを決めるパラメータだとみることができる。閾値を1にすればunconditionalなpermutationになる。閾値を0にすると、\(X_k\)以外のすべての予測子が\(Z^{(t)(s)}_{(-k)}=Z^{(t)}_{(-k)}\)に含まれるようになり、パーティションは極端に小さくなり、CPIの値はゼロに近づく(あまり意味がなくなる)。partyでは閾値のデフォルトを0.2としていたが今にして思えば小さすぎた。
　結局どうしても、閾値の影響はサンプルサイズで変わってきてしまう。サンプルサイズが大きいときは閾値を上げないとpartial寄りになる。

シミュレーション
　[ちゃんと読まずに飛ばしたけど、閾値は0.8以下であれば動かしてもCPIの全体的な大きさが変わるだけでパターンは変わらない、つまりパーティションを細分化しすぎることによる弊害はなさそう、とかなんとか書いてあった。ふーん]

考察
　閾値は1を含めていくつか試してみることをお勧めする。permimpのデフォルトは0.95にした。0.8以下を試す必要はなさそう。サンプルサイズが大きいときは閾値を高めにするようにね。

　今後の課題:

• 欠損データの扱い。当面は多重代入をお勧めする。
• シミュレーションの拡張。回帰問題でやったけど、分類問題だとどうなるか、とか。
• mtryなどのハイパーパラメータの効き方。考え方としては、ハイパーパラメータはあくまで予測の正確性という観点から最適化すればよくて、PI/CPIの算出はそのあとの話なのだけれど。
• 木の本数(ntree)。RFに本来必要な本数よりも多めにしたほうがよさそう。
• パーティションが細かすぎる問題。現実にはあまり問題にならないとは思うけど。
• パーティションのための分割点を木から得るとき、ある深さまでしかみない、という手があるんじゃない？って査読者からいわれたけど、仰る通りでございますね。

　そもそもの問題として、「変数重要性」とは正確にはなんなのか、いまだ合意がない。本論文ではpartial/marginalを区別した。手元のリサーチ・クエスチョンに照らしてどちらの視点が重要かを考えるように。

云々。
——————–
いやー、もうね… party::varimp()の実装についてめっちゃカジュアルに自己批判してて、もうびっくりですよ。小松政夫の往年のギャグ「悪いね、悪いね、ワリーネ・ディートリッヒ」を思い出してしまった。いや、正直であることは素晴らしいですけれども…