コンテンツにスキップ

一般化ランダムフォレストと因果推論への応用

1. ランダムフォレストの概観

1.1 Leo Breiman によるランダムフォレストの提案から、2012年までのランダムフォレストの理論解析の歩み

Breiman (2001)年に提案したランダムフォレスト（Random Forest）は、複数の決定木をランダムに生成し、それらの予測を平均する強力なアンサンブル学習法である。提案当初は、クラス分類問題へのランダムフォレストの応用がメインであったため、決定木がアンサンブルの弱学習器となっていたが、現在では回帰の文脈でも用いられるため、この論文では回帰木と呼ぶことにする。ランダムフォレストの基本的な構成要素は2つのランダム化メカニズムによる。1つ目はデータのサブサンプリングである。各回帰木の訓練において、訓練データセットからランダムに抽出されたサブサンプルを用いる。このサブサンプリングは、ブートストラップサンプリングを用いる場合もあれば、各サンプルは1回しかとらないという方法を用いることもある。もう1つは特徴量のランダムな選択である。回帰木を学習する際に、各ノードの分割で用いることのできる変数をランダムに選択するというメカニズムである（この点に日本語の文献の説明では誤りがある場合がある。各決定木の学習の際に用いる変数をランダムに選択すると書かれているが、この表現は誤解を招くので適切ではない）。2001年に提案されて以来、ランダムフォレストは応用の文脈では高次元データへの適応力や汎用性の高さから様々な応用が行われている。（Belgiu and Dragut,2016, Parmar et al.2018, Hu and Szymczak, 2022, Iranzad and Liu, Tyralis et al.,2019, Biau and Scornet, 2016, Chen and Ishwaran, 2012, Walker et al.,2022）

しかし、木の構造がデータに強く依存し、さらにランダム化要素も含むため、その数学的解析は困難であり、広く使われているにもかかわらず理論的性質の多くがわかっていないと指摘されてきた。例えば、Breiman (2004) においても理論的な解析を試みているが、シンプルなランダムフォレストに対する性質の解析の留まり、一般的なランダムフォレストの理論解析には至っていない。ランダムフォレストの理論研究が大きく進んだのは2010年以降であるが、その理論解析に大きな貢献を果たしているのは、Lin and Jean (2006)の研究である。Lin and Jean (2006) では、回帰木によって生成される木構造によって、点\(x\)における回帰木の推定値の計算が、k-近傍法と見なすことができることを示した。この結果は、Biau and Devroye (2010)によってランダムフォレストの理論解析へと拡張されている。そして、2010年以前の結果は、2014年以降に続くランダムフォレストの理論的な結果の大きな礎となっている。また、同時期にはランダムフォレストを用いて、条件付き平均関数以外を推定することに取り組む研究も発展している。例えば、Meinshausen (2006) では、ランダムフォレストによる条件付き分位点関数の推定法である quantile regression forest を提案し、その一致性に言及されている。また、Ishwaran et al.,2008 では、累積ハザード関数が共変量によってどのように変動するのかを明らかにするための random survival forest を提案している。

また、ランダムフォレストが注目された1つの理由に変数重要度の概念がある。変数重要度は、Breiman (2001)が木構造の特徴を活かして提案したランダムフォレストにおいて、どの特徴量が予測寄与度が大きいかを定量化するための方法である。ランダムフォレストには、Mean Decrease Accuracy (MDA) と Mean Decrease Impurity (MDI) の2つが変数重要度として利用される。MDAは、ランダムフォレストを一度学習させた後で、ランダムフォレストの予測値の\(R^2\)決定係数と、特定の変数をランダムにシャッフルした場合のランダムフォレストの予測値の\(R^{2}\)の決定係数を計算し、その差によって変数の重要度を定量化する。一方で、MDIは各回帰木の学習において変数による分割で不純度が改善する量のノード内サンプルの重み付け和として定量化し、それを全ての回帰木で平均かすることで変数の重要度を計算する。これら2つの変数重要度についても、Breiman (2001) において提案されて以来、予測に関係する変数の発見で用いられてきた。しかしながら、ランダムフォレストの理論基盤が確立していないこともあり、変数重要度に対する理論もまた確立されていなかった。よって、近年までの変数重要度の結果からは統計的な推測の妥当性などは述べることができていない。しかし、2020年から変数重要度に関する理論的な解析結果の多くが示されており、従来経験的にしかわかっていなかった変数重要度の問題点の解消に新たな光が当てられている。

1.2 2012年以降のランダムフォレストの理論解析の歩み：一致性と、漸近正規性の証明

Biau (2012)​ はBreimanの提案したアルゴリズムに近いランダムフォレストモデル（例えば各木でランダムに特徴次元を選ぶ方法など簡略化した設定）について初めて厳密な解析を行いました​。この研究では、「ランダムフォレストは一貫して（consistent）動作し、またスパース性（不要なノイズ特徴が多数存在する状況）に適応する」ことが示されています​。スパース性への適応とは、真に有効な特徴変数がごく一部であっても、その収束レート（予測誤差の減少速度）が有効特徴の次元数にのみ依存し、無関係な特徴の数には依存しないことを意味します​。これは、ランダムフォレストが高次元でも不要な変数を無視し、本質的な変数に集中できることを直感的に裏付ける結果です。

Scornet et al. 2015の貢献

続いて、Erwan Scornet と Gérard Biau らによるさらなる発展があります。Scornet・Biau・Vert (2015)​ では、Breiman (2001) のオリジナルのランダムフォレストアルゴリズムにおける（非ブートストラップサンプリングと、\(m_{try}\) パラメータによるランダムな特徴選択による分割）ランダムフォレスト推定量の一致性が、推定対象となる真の関数が特徴量に対して加法的である場合に、一致性が成り立つことを示しました。

  • 加法性の仮定 $$ Y = \sum_{j=1}^{p}m_{j}(X^{(j)})+\varepsilon \quad \mathrm{where} \quad X = (X^{(1)},...,X^{(p)}) \sim U([0,1]^{p}),\quad \varepsilon \sim N(0,\sigma^2) $$ ここで、Scornet・Biau・Vert (2015)​の一致性の証明に限らず、一致性や漸近正規性の証明では、ランダムフォレストを構成する木のパラメータをどのように制御するかが本質的であり、損失関数をどのように設定しているかは一致性や漸近正規性に対して影響を与えていません。この点には注意が必要です。

Scornet・Biau・Vert (2015)においては、木に対して次の仮定が置かれます。まず、木を構成する際のサブサンプルサイズ\(a_n\)は、観測されたサンプル\(n\) に対して、無限大に発散します\(a_n \rightarrow \infty\). また、木を構成する葉の数\(t_n\)も、\(t_n \rightarrow \infty\)を満たしますが、この発散速度には制約があり、\(t_n (\log a_n)^9/a_n \rightarrow 0\) となることが条件として与えられます。つまり、サンプルの発散速度に対して、木の葉を増やす速度はそれよりも遅い速度が要求されます。

また、\(m_{try} = p\)の状況下において、真の関数が、観測された変数\(p\)個のうち、モデルに含まれる変数が\(s\)個である状況を考えます。 $$ Y = \sum_{j=1}^{s} m_j (X^{(j)}) + \varepsilon $$ このとき、任意の\(m_j\) が 任意の\(Y\)の区間\([a,b]\)において定数関数でないならば、十分高い確率で木の分割変数は\(\{1,2,...,s\}\)から選ばれることを示しました。これは、ランダムフォレストが高次元の観測のもとで有効に動作することを示す根拠の1つと捉えることができます。

Wager et al., 2014-2018の貢献

次に、Stefan Wager が Stanfordの研究グループで取り組んだランダムフォレストの一致性と漸近正規性、および漸近分散の導出、さらにオリジナルのcausal forestの提案までの流れを説明します。まず、ランダムフォレストの一致性についてWager and Walther (2014) では、Scornet et al.,(2015)とは異なる仮定をおいて一致性を示しています。まず、木の葉の数\(k_n\)に対して、 $$ \lim_{n\rightarrow \infty}\frac{\log(n)\max \left(\log(d), \log\log(n)\right)}{k} = 0 $$ の仮定および Sparse signal の仮定をおく。すなわち、\(p\)次元の特徴量\(X\)に対して、有効な次元の集合\(\mathcal{Q} =\{1,2,...,s\}\)次元で、\((Y,X_{\mathcal{Q}}) \perp \!\!\! \perp (X_{-\mathcal{Q}})\) が成り立つ。 さらに、Monotone signal の仮定をおく。この仮定は、\(j \in \mathcal{Q}\)（有効な次元）に対して、\(x_{(-j)} \in [0,1]^{d-1}\) を固定した時に、以下の式を満たすような最小の効果\(\beta > 0\)が存在することである。 

    \left|\mathbb{E} \left[ Y_i \mid (X_i)_{-j} = x_{-j}, (X_i)_{j} > \frac{1}{2} \right] - \mathbb{E}\left[ Y_i \mid (X_i)_{-j} = x_{-j}, (X_i)_{j} \leq \frac{1}{2} \right] \right| \geq \beta
さらに、\(E[Y|X=x]\) に対するLipchitz連続性を仮定したもとで、以下で定義されるGuess-and-Check forestは、一様一致性を持つ。Guess-and-check forest は次のように定義される。

Guess-and-Check Forest

Input

  • $ n $ 個の訓練データ \((X_i, Y_i)\)
  • 最小葉ノードサイズ $ k $
  • バランスパラメータ $ 0 < \alpha < 1/2 $

アルゴリズム

Guess-and-Check Tree は、以下の分割手順を再帰的に適用し、分割が不可能になるまで処理を行う。すなわち、全ての終端ノードの訓練データ数が $ 2k $ 未満になるか、そもそも (12) の条件を満たす分割が存在しない場合まで繰り返す。

  1. ノード $ \nu $ を選択し、そこに少なくとも $ 2k $ 個の訓練データが含まれることを確認する。
  2. 候補となる分割変数 $ j \in {1, \dots, d} $ を一様ランダムに選択する。
  3. 最小二乗誤差の分割点 $ \hat{\theta} $ を選択する。より具体的には、 $$ \hat{\theta} = \arg\max_{\theta} \ell (\theta) $$ ここで、目的関数 $ \ell(\theta) $ は次のように定義される。 $$ \ell (\theta) := \frac{4 N^-(\theta) N^+(\theta)}{(N^-(\theta) + N^+(\theta))^2} \Delta^2(\theta) $$ 分割点 $ \theta $ は、ノード $ \nu $ に属するサンプル $ X_i $ の成分 $ (X_i)j $ のいずれかに対応する値とする。 $$ \alpha \times | { i : X_i \in \nu } |, k \, \leq\, N^-(\theta), \quad N^+(\theta) $$ また、以下の定義を用いる。 $$ \Delta(\theta) = \frac{1}{N^+}\sum{{ i : X_i \in \nu(x), (X_i)j > \theta }} Y_i - \frac{1}{N^-}\sum Y_i $$ ただし、 $$ N^-(\theta) = | { i : X_i \in \nu, (X_i)_j \leq \theta } |, \, N^+(\theta) = | { i : X_i \in \nu, (X_i)_j > \theta } | $$
  4. 分割条件の判定:
    • 変数 $ j $ に対して、すでに成功した分割が存在する場合。
    • または、以下の条件が満たされる場合。 $$ \ell (\hat{\theta}) \geq \left( 2 \times 9M \sqrt{\frac{\log(n) \log(d)}{k \log((1 - \alpha)^{-1})}} \right)^2 $$ この場合、$ j $ 番目の変数に対してノード $ \nu $ を $ \hat{\theta} $ で分割する。そうでない場合は、ノード $ \nu $ はこの時点では分割されない。

Guess-and-Check Forest

Guess-and-Check Forest は、上記のアルゴリズムに従って 独立に生成された \( B \) 本の Guess-and-Check 木の平均 を取ることで構築される。

次は、Random forestの漸近分散の推定法に関する結果である。まず、漸近正規性の結果に先立って、Wager, Hastie and Efron (2014) では、Efron (2014) においてモデル選択によるデータ駆動型のアプローチがもたらす予測誤差の過小評価を補正する方法として提案した Infinitesimal Jacknifeによる方法を、ランダムフォレストに拡張し、ランダムフォレストによる推定量のばらつきを評価する方法を提案した。 \(b=1,2,...,B\)番目の回帰木による推定値を\(\hat{f}_{b}(x)\)とし、 \(N_{ib} \in \{0,1\}\) が \(b\)番目の回帰木（Double Sample Treeの場合には、いずれか一方のグループに含まれる場合に\(N_{bi}=1\)と定義する）に含まれるかどうかを表す指示関数とする。 $$ \hat{V}{IJ}(x) = \frac{n-1}{n} \left(\frac{n}{n-s}\right)^2 \sum^{}^{n}\mathrm{Cov}\left[\hat{f}_{b}^{}(x), N_{ib}^{}\right]^2 $$ これに対して、有限標本化での推定量は $$ \hat{V}{IJ}^{B}(x) = \frac{n-1}{n} \left(\frac{n}{n-s}\right)^2 \sum^{}^{n}\frac{1}{B}\sum_{b=1}^{B}\left(N_{bi}-\frac{s}{n}\right)\left(t_{b}^{}(x)-\bar{t}^{}(x)\right) $$ となる。ただし、この結果はEmpiricalな評価にとどまっており、実際に漸近分散と一致することを示したのは、このあとのWager and Athey (2019)である。漸近正規性の議論に移る前にもう1つの重要な論文は、Athey and Imbens (2016) において指摘された木構造モデルの学習におけるバイアスの問題である。木構造モデルを学習する際には、テストデータへの当てはまりがよくなるように学習を行う必要がある。しかし、一般的なCARTのような学習では訓練データを用いて木構造の学習と木による予測値の両方を学習するため、実はテストデータに対しての当てはまりを最適化しているのではないということを示している。また、この論文では、ランダム化比較試験に対する介入効果の異質性を発見する回帰木としてCausal Treeを提案している。従来のCART的なアプローチでは、共変量\(X\)の分割と分割によって生成された葉の値には相関があるため、相関に依存した因果効果が検出されてしまう。そこで、Honest（誠実性） と呼ばれる概念を導入することで、この問題を解消する。誠実性とは、データを2つに分け、一方で木構造を学習させ、もう一方で予測値の計算と検定を行えば、分割と予測値は独立になるため、効果の検証ができるということである。ここで、導入された誠実性が次に示す漸近正規性を示すための1つの鍵となる。

Wager and Athey (2019) では、次の性質を満たすbagging treeに対して、有限標本化でのバイアスを評価し、漸近正規性を示した。 - サブサンプルを2つに分割し、サブサンプルの片方を木構造の学習に、もう一方を葉の推定値の計算に用いる Double Sample Tree を 弱学習器として採用する - 任意の特徴量\(X^{(j)}, \, j=1,2,...,d\)が、ノードの分割で選択される確率が\(0\)ではない。 - subsample size \(s_n\) は、サンプルサイズ\(n\)、次元\(d\)、および分割を制御するパラメータ\(\alpha \leq 0.2\) に対して、次の関係式を満たす。 $$ s_n \approx n^{\beta}, \qquad 1-\left(1+\frac{d}{\pi}\frac{\log(\alpha^{-1})}{\log\left((1-\alpha)^{-1}\right)}\right)^{-1} < \beta < 1 $$ - \(\mu(x) := E[Y|X]\) が リプシッツ連続 である。

また、この論文では上記に述べた\(\hat{V}_{IJ}\)がランダムフォレストの漸近分散に対する一致推定量であることが示されており、ランダムフォレストによる推定量の分散を与えた点で画期的な論文であった。

これに加えてWager and Athey (2019)では、因果効果を推定する方法として、causal forestを提案している。このcausal forestは現在はほとんど採用されないが、概念だけ紹介しておく。causal forest は、Athey and Imbens (2016) で提案された causal tree を弱学習器として用いるランダムフォレストである。causal tree は、回帰におけるノードの分割基準を、処置効果推定の文脈へと拡張したものである。ただし、単純に拡張したものではなく、テストデータへの当てはまりを最適化するような基準を用いている。 $$ \begin{align} -\widehat{\text{EMSE}}{r} \left( \mathcal{S}^{\text{tr}}, N^{\text{est}}, \Pi \right) &\equiv \frac{1}{N^{\text{tr}}} \sum, \Pi \right) \ &\quad - \left( \frac{1}{N^{\text{tr}}} + \frac{1}{N^{\text{est}}} \right) \cdot \sum_{\ell \in \Pi} \left( \frac{S^2_{\mathcal{S}^{\text{tr}}}^{\text{tr}}} \hat{\tau}^2 \left( X_i, \mathcal{S}^{\text{tr}{\text{treat}}} (\ell)}{p} + \frac{S^2 \right). \end{align}^{\text{tr}}_{\text{control}}} (\ell)}{1 - p} $$ この基準は、第1項が葉の推定値間の分散を最大化させる項であり、異質性を捉えるための項となっている。一方で、第2項は葉内の処置群および対照群の分散の和であり、当てはまりに対する罰則項である。つまり、予測の異質性を高めつつも、葉の中のサンプルが極端に少なくなったりすることで予測全体の分散が大きくなりすぎないように制御していると捉えることができる。

一般化ランダムフォレストの一致性・漸近正規性

一般化ランダムフォレスト (Athey, Tibshirani and Wager, 2020) は、ランダムフォレストが回帰と分類の問題にしか対象としないのに対して、統計モデルなどで定義されるパラメータの推定を対象するランダムフォレストの拡張の1つである。例えば、後にも扱う因果効果は、データから直接推定することはできないため、通常のランダムフォレストでは因果効果を直接的に推定することはできない。これに対して、一般化ランダムフォレストは、推定方程式をランダムフォレストで直接取り扱う方法を提供する。

この枠組みを提供するために、まずはランダムフォレストに対する視点転換を行い、ランダムフォレストをカーネル重み付け関数の推定法として捉える。いくつかの記号の導入を行い、その後でランダムフォレストのカーネル重み付け推定としての側面を説明する。

問題設定として、\(p\)次元の共変量 \(X\) と 1次元の結果変数 \(Y\) を考え、\((X,Y)\)の従う確率分布を\(P_{XY}\)とする。次に、\(P_{XY}\)からの独立同一な観測の組を \(D_{n} = (X_{i}, Y_{i}), i=1,2,...,n\)とする。いま、観測データ\(D_n\)上で学習された木による点\(x\)の予測を\(T(x;D_n)\)とし、木の葉のうち点\(x\)を含むものを\(L(x)\)とする。すると、 $$ T(x;D_n) = \sum_{i=1}^{n} \frac{1{X_i \in L(x) }}{|L(x)|}\cdot Y_i $$ ただし、\(|L(x)|\)は点\(x\)を含む葉に含まれる観測データのサンプルサイズである。この式は、点\(x\)に対する回帰木による予測とは、点\(x\)を含む葉に存在する観測データに対して、\(1/|L(x)|\)の重み付けて和を取ることで得られることを意味している。

同様にして、ランダムフォレストを考えてみる。いま、サイズ \(n\) のサンプルから、サイズ \(s_n\)のサブサンプル を \(B\)回とり、それぞれを \(D_{n}^{(b)}, b=1,2,...,B\) をとする。ランダムフォレストは、サブサンプル \(D_{n}^{(b)}\) 上に、回帰木を学習させ、各回帰木の出力の平均値を テスト点\(x\)の推定値として返す。各回帰木を\(T_{b}\)とし、点\(x\)に対する回帰木\(T_{b}\)による予測を\(T_{b}(x;D_n, \xi)\)する。ここで、\(\xi\) は回帰木を学習させる際のランダムネスに対応する。例えば、splitの候補となる変数の集合をランダムに選ぶことなどに対応する確率変数である。いま、回帰木\(T_b\)において 点\(x\)を含む葉を\(L_{b}(x)\)とする。すると、先ほどと同様の議論から、 $$ \alpha_{bi} = \frac{1{X_i \in L_{b}(x) }}{|L_{b}(x)|}, \qquad T(x;D_n) = \sum_{i=1}^{n} \frac{1{X_i \in L_{b}(x) }}{|L_{b}(x)|}\cdot Y_i = \sum_{i=1}^{n}\alpha_{bi}(x)Y_i. $$ となる。ランダムフォレストによる点\(x\)の予測値 \(RF_{B}(x; D_n)\)は、回帰木の平均であるから、 $$ RF_{B}(x;D_n) = \frac{1}{B}\sum_{b=1}^{B}\sum_{i=1}^{n}\alpha_{bi}(x)Y_i = \sum_{i=1}^{n}\left{\frac{1}{B}\sum_{b=1}^{B}\alpha_{bi}(x)\right}Y_i = \sum_{i=1}^{n}\alpha_{i}(x)Y_i $$ ここで、 $$ \alpha_{i}(x) = \left{\frac{1}{B}\sum_{b=1}^{B}\alpha_{bi}(x)\right} $$ とおいた。\(\alpha_{i}(x)\)は、回帰木それぞれが点\(i\)を含むかどうかを確率変数とみなしたとき、\(B\)個の確率変数の平均と解釈できる。つまり、点\(x\)を同じ葉に含まれやすいサンプル\(i\)に対しては、大きな重みがかかり、同じ葉にほとんど含まれないサンプル\(i\)に対しては小さな重みがかかる。

一方で、\(\hat{\theta}_{B}(x) = RF_{B}(x;D_n)\) とおくと、これは以下の推定方程式の解として解釈することができる。 $$ \hat{\theta}{B}(x) = \arg\min(x)(Y_i - \theta)^2 $$ よって、ランダムフォレストは回帰木のアンサンブルと捉えることもできるが、回帰木によってテスト点}\sum_{i=1}^{n}\alpha_{i\(x\)と観測データの関係の強さを\(\alpha_{i}(x)\)として定量化し、それによって推定方程式を解く方法と捉え直すことができる。この視点の切り替えによって、ランダムフォレストが\(\alpha_{i}(x)\)をカーネル重み付け関数としたカーネル重み付け推定量であることが理解できる。この性質を用いて、ランダムフォレストでパラメータの推定を行うというのが一般化ランダムフォレストである。

一般化ランダムフォレストと局所推定方程式

確率変数の組 \(Z = (X, O)\) と考え、これらが従う分布を\(P\)とおく。いま、\(P\)からの独立同一なサンプルを\(Z_i, (i=1,2,...,n)\) とする。ただし、\(O\)は回帰であれば \(Y\) であるが、因果推論などにおいては処置変数\(T\)と結果変数\(Y\)の組 \(O=\{T, Y\}\) となる場合もあり、問題によって\(O\)は読み替えるものとする。いま、スコア関数 \(\psi(\cdot)\) によって特徴づけられる局所推定方程式 $$ E[\psi_{\theta(x),\nu(x)}(O_i)\mid X_i = x] = 0 $$ を考え、このスコア方程式の解として定義されるパラメータ \(\theta(x), \nu(x)\) に興味があるとする。例えば、\(\theta(x) = E[Y\mid X=x]\) に興味がある場合には、\(\psi(O) = \psi_{\theta(x)}(O) = Y-\theta(x)\)とする。また、統計的因果推論の文脈において、観測される結果変数\(Y\)が、共変量\(X\)と処置変数\(T \in \{0,1\}\)によって、\(Y = \mu(X) + T \cdot \tau(X) + \varepsilon\) とモデル化され、\(E[\varepsilon \mid X] = 0\) および \(\varepsilon \perp\!\!\!\perp T \mid X\) を満たすとする。このとき、\(E[\varepsilon \mid X] = 0\) および、\(E[\varepsilon \cdot T \mid X] = 0\) であるから、以下のようになる。 $$ \psi_{\theta(x),\nu(x)}(O) = \psi_{\mu(x),\tau(x)}(T,Y) = \left(Y - \mu(x) - T\cdot \tau(x)\right)\begin{pmatrix}1 \ T \end{pmatrix} $$

次に、局所推定方程式によって特徴づけられるパラメータの推定について考える。いま、分布\(P\)からの独立で同一なサンプルを\((X_i, O_i), i=1,2,...,n\)とする。このとき、パラメータ \(\theta(x)\) の推定量は、一般的には点\(x\)と\(z\)の距離 \(u = |x-z|\)に依存する、バンド幅を\(h\)とするカーネル重み付け関数 \(K_{h}(u)\) に基づいて重み付け推定方程式を解くことによって構成される。 $$ \hat\theta(x) = \arg\min_{\theta,\nu}\left|\sum_{i=1}^{n}K_{h}(X_i - x)\psi_{\theta,\nu}(O_i)\right|_{2} $$ この重み付け推定方程式の解 \(\hat\theta(x)\) は、適当な正則条件のもとで一致性と漸近正規性を満たすことはよく知られている。

ランダムフォレストによって構成される \(\alpha_{i}(x)\) は、点\(x\) とサンプル\(i\)の関連性の強さの指標であり、これはカーネル重み付け関数 \(K_{h}(X_i - x)\) に概念的に対応づけることができる。よって、直感的には \(\alpha_{i}(x)\) がカーネル重み付け関数と同様の性質を満たすのであれば、ランダムフォレストによって得られる推定方程式の解もまた、一致性と漸近正規性を満たすということがわかる。

仮定と理論的な背景

ここで、理論的結果を述べる上で重要な指摘を行なっておく。現状、Wager and Athey (2019) や、Athey et al. (2020) などで述べられている理論的な結果では、回帰木を構成する際の損失関数に対しての制約はなく、これから述べるような仮定を満たす回帰木を用いていればよい。よって、あるパラメータに対する推定法を構成する上で最適な損失関数の選択などについては、ここでは全く議論されていないことに注意して欲しい。これらの推定における細やかな指針については本節の後半で紹介する。

まずは回帰木の構成に関する仮定をおく。 - (A1) 回帰木の構成はサンプルの順序に依存しない（訓練データの添字には依存しない, i.i.d. であれば関係はない） - (A2) ノードを分割する際には、分割によって生成される子ノードが、分割元のノードに含まれるサンプルの\(0 < \omega < 0.5\) 以上の割合を含む（極端な分割を抑制する） - (A3) ノードの分割において、任意の次元 \(j\) が分割される確率 \(\pi > 0\) 以上である（つまり、どの次元\(j\)も分割される可能性がある）。 - (A4) 任意の回帰木の葉は\(k\)以上、\(2k-1\)以下のサンプルを含む \((k \geq 1)\)。 - (A5) Double sample trees の回帰木の構成方法を利用する。ただし、random forest を構成する際のサブサンプルのサイズ \(s_n\) は、\(s_n \rightarrow \infty\) かつ \(s/n \rightarrow 0\) を満たす。

これらの仮定は、Wager and Athey (2019) で用いられた漸近正規性を満たすランダムフォレスト推定量のための仮定と同じである。ここで、(A3) の実装においては、min(max{Poisson(m),1},p) を grf のパッケージでは用いることで、どの変数にも分割される確率が存在することを保証している。

さらに、局所推定方程式に対して以下の仮定を記述するために、 $$ M_{\theta,\nu}(x) := E[\psi_{\theta,\nu}(O) \mid X=x] $$ とおく。

  • (B1) 固定された\((\theta,\nu)\)に対して、\(M_{\theta,\nu}(x)\) は \(x\) に対して連続である。
  • (B2) \(x\)を固定した時、\(M_{\theta,\nu}\)は、\((\theta,\nu)\) について2階連続微分可能であり、\(V_{\theta,\nu}(x) := \partial/\partial(\theta,\nu)M_{\theta,\nu}(x)\mid_{\theta(x),\nu(x)}\) が、任意の\(x \in \mathcal{X}\) に対して正則である。
  • (B3) \(\gamma\) を以下で定義されるバリオグラムとする $$ \gamma\left(\begin{pmatrix}\theta \ \nu\end{pmatrix}, \begin{pmatrix}\theta' \ \nu'\end{pmatrix}\right) := \sup_{x\in\mathcal{X}}\left{| Var[\psi_{\theta,\nu}(O_i) - \psi_{\theta',\nu'}(O_i)]\mid X=x|_{F}\right} $$ このとき、\((\theta,\nu)\)について、\(\gamma\) はリプシッツ連続である。
  • (B4) \(\psi\)が、\(\psi_{\theta,\nu}(O) = \lambda(\theta,\nu ; O_i) + \zeta_{\theta,\nu}(g(O_i))\) の形状に分解される。ただし、\(\lambda\) は \((\theta,\nu)\)についてリプシッツ連続かつ、\(g : \{O_i\} \rightarrow \mathbb{R}\)、および \(\zeta_{\theta,\nu} : \mathbf{R} \rightarrow \mathbf{R}\) の単調かつ有界な関数である。
  • (B5) \(\sum_{i}\alpha_i = 1\) を満たす重み\(\alpha_i\) に対して、重み付け推定方程式の解\(\hat\theta, \hat\nu\)が、ある定数\(C\)に対して、\(\|\sum_{i=1}^{n}\alpha_{i}\psi_{\hat\theta,\hat\nu}(O_i)\|_{2} \leq C \max\{\alpha_i\}\) を満たす。
  • (B6) スコア関数 \(\psi_{\theta,\nu}(O_i)\) は、凸関数の負の劣導関数となっており、\(M_{\theta,\nu}\)は、狭義凸関数の負の導関数となっている。

これらの仮定は証明において本質的であるが、回帰の問題、分位点推定の問題、2値の処置変数に対する因果パラメータの推定、2値の処置変数と2値の操作変数に対する操作変数法による因果パラメータの推定などのさまざまな設定に対して満たされる。また、これらの仮定の中で大きな特徴なのは、スコア関数の期待値である\(M\)が微分可能であれば、\(\psi\)は微分可能ではなくてよいというところである。これによって、スコア関数として不連続な関数となる分位点回帰についても一般化ランダムフォレストは理論的な保証を与えることができる。

一般化ランダムフォレストの一致性

以上の仮定のもとで構成された回帰木のアンサンブルによって生成される重みを \(\alpha_{i}(x)\)とするとき、ランダムフォレストによる重みづけ推定方程式 $$ (\hat\theta(x),\hat\nu(x)) = \arg\min_{\theta,\nu}\left|\sum_{i=1}^{n}\alpha_{i}(x)\psi_{\theta,\nu}(O_i)\right|{2} $$ の解 \((\hat\theta(x),\hat\nu(x))\) は、\((\theta(x),\nu(x))\) に対する一致推定量となる。この証明の中で最も重要なのは、ランダムフォレストの生成する重みの性質である。Wager and Athey (2018) においては、\(\alpha_{i}(x)\)の局所化が、上記の仮定を満たす回帰木では次のレートで起こることを示した。 $$ E[\sup{|X_i - x|\right) $$ この局所化は、点} : \alpha_{i}(x) > 0}] = O\left(s^{-\frac{\pi}{2}\cdot \frac{\log((1-\omega)^{-1})}{\log(\omega^{-1})}\(x\)を推定する際に、L2距離でどの程度離れたサンプルまで重みがかかるのかを示している。\(\omega\) は分割を制御するパラメータで、\(\omega = 0.5\) の場合はサンプルを\(1:1\)に分割する点のみを分割点として採用することになり、\(\omega\)の値が\(0\)に近づくほどより柔軟なノードの分割が可能になる。このとき、\(0 < \omega < 0.5 \(に対して、\)\frac{\log((1-\omega)^{-1})}{\log(\omega^{-1})}\)は単調増加で、\([0,1]\)の値を取るので、極端な分割を許容するとランダムフォレストの生成する重みの裾が\(x\)の周辺への収束するのが遅くなることがわかる。

しかし、いずれにせよ\(\alpha_{i}(x)\)が\(x\)の周辺に収束するということは、重みづけ推定方程式が\(\sum_{i=1}^{n}\alpha_{i}(x)\psi_{\theta,\nu}(O_i)\) が \(E[\psi_{\theta,\nu}(O)\mid X=x]\) に収束することに対応するため、任意の \(x\in \mathcal{X}\)得られる推定方程式の解も \((\hat\theta,\hat\nu) \rightarrow (\theta(x),\nu(x))\) の確率収束が成り立つことを意味している。これらは直感的な説明であるが、詳細な証明の過程は経験過程を用いた議論となる。

一般化ランダムフォレストの漸近正規性

一般化ランダムフォレストの漸近正規性についても、ここでは証明の直感的なプロセスを紹介する。いま、パラメータ\(\theta(x)\)に対する \(i\)番目のサンプルに対する影響関数を $$ \rho_{i}^{}(x) := - \xi^{T}V(x)^{-1}\psi_{\theta(x),\nu(x)}(O_i) $$ とおく。次に、擬似的なランダムフォレストの推定量を $$ \tilde{\theta}^{}(x) := \theta(x) + \sum_{i=1}^{n}\alpha_{i}(x)\rho_{i}^{*}(x) $$ と構成する。この擬似的なランダムフォレスト推定量は、\(\theta(x) + \rho_{i}^{*}(x)\)を結果変数とする重み\(\alpha_{i}(x)\)による回帰フォレストによる推定量となっている。この\(\tilde{\theta}^{*}(x)\)を用いて次の2つを示す。

  1. \((\tilde{\theta}^{*}(x)-\theta(x))\) が適当な収束レートのもとで漸近的に正規分布に収束すること
  2. \((\tilde{\theta}^{*}(x) - \hat\theta)\) が \((\tilde{\theta}^{*}(x)-\theta(x))\) の収束レートよりも速いレートで\(0\)に確率収束すること

この2つの結果が得られれば、ランダムフォレストによる重みづけ推定方程式によって得られる \(\hat\theta\) が \(\theta(x)\) に対して漸近正規性を持つことを示すことができる。まず、1つ目の結果は infinite-order U統計量の議論となる。これは、Mentch and Hooker (2016) の結果から示すことができる。2つ目の結果が、Athey et al.(2020) の結果の1つで Lemma 4として記されているものである。このとき、収束のレートをコントロールするために必要なのが、ランダムフォレストを構成する際のサブサンプルサイズ \(s_n\)である。Athey et al. (2020) の結果では、\(s_n = n^{\beta}\)の指数部について、以下のレートを要求している $$ \beta_{\min} = 1 - \left(1 + \frac{\pi^{-1}\log(\omega^{-1})}{\log((1-\omega)^{-1})}\right)^{-1} < \beta < 1 $$ このレートは先ほどの重み\(\alpha_{i}(x)\)でみた対数の逆数に似た形状の項が左片にあり、\(\omega\)が\(0\)に近いほど、\(\beta\)の下限が\(1\)に近づくようになっている。このサブサンプルレートのもとで、以下の結果が得られる。

$$ \frac{\hat\theta(x) - \theta(x)}{\sigma_{n}(x)} \rightarrow N(0,1), \quad \sigma_{n}^{2}(x) = (s_n/n)\cdot\mathrm{polylog}(s_n/n)^{-1} $$ ただし、\(\mathrm{polylog}(u) = a_{k}(\log(u))^{k} + \cdots + a_{1}(\log(u))^{1} + a_{0}\)を意味する。よって、この結果から一般化ランダムフォレストに基づく推定量が漸近正規性を持つことがわかる。ここで重要なのは、\(\sigma_{n}^{2}(x)\)の収束速度である。

簡単な例でランダムフォレストの収束レートを考える。いま、共変量を\(1\)次元であると仮定する。すなわち\(\pi = 1\)である。このとき、\(\beta\)の下限\(\beta_{\min}\)は、\(0 < w < 0.5\) に対して、\((1, 0.5]\)の区間での単調減少関数となる。一方で、共変量の次元を\(2,3,4\)と増やしていくと、\((1,0.66]\)、\((1,0.75]\)、\((1,0.8]\)と下限の\(\beta\)の値が大きくなる。ここで、ランダムフォレストの収束レートは\(n^{1-\beta}\)に依存するので、この結果から共変量の次元に対して収束のレートが悪くなることがわかる。

実際、ランダムフォレストは探索する次元の大きさに対して、分割可能な回数がサンプルサイズに依存するため、サンプルサイズが小さい場合には点\(x\)を含む葉が大きくなり、それだけバイアスも大きくなる点には注意が必要である。以上の点から、ランダムフォレストを利用する場合には、次元の大きさとサンプルサイズには十分に注意が必要である。

一般化ランダムフォレストの漸近分散推定

次に、一般化ランダムフォレストに対する漸近分散の推定を考える。Wager and Athey (2018) での議論では、Infinitesimal Jacknifeを用いてランダムフォレスト推定量の漸近分散を推定したが、ここでは一般化ランダムフォレストの推定量が推定方程式の解として得られることを用いてデルタ法を応用することで漸近分散の推定を考える。

まず、上記の漸近正規性の議論から、漸近分散に寄与する項が\((\tilde{\theta}^{*}(x)-\theta(x))\)の部分であるから、 $$ \frac{\mathrm{Var}[\tilde{\theta}^{}(x)]}{\sigma_{n}^{2}(x)} \rightarrow 1 $$ であるから、\(\tilde{\theta}^{*}(x)\) の分散について考えれば良いということになる。いま、 $$ H_{n}(x; \theta,\nu) = \mathrm{Var}\left[ \sum_{i=1}^{n}\alpha_{i}(x)\psi_{\theta,\nu}(O_i) \right] $$ とおくと、\(\tilde{\theta}^{*}(x)\)の定義から、 $$ \mathrm{Var}[\tilde{\theta}^{}(x)] = \xi^{T}V(x)^{-1}H_{n}(x; \theta,\nu)V(x)^{-T}\xi $$ となる。各項の計算について見ていくと、まず\(V(x)^{-1}\)については、\(V(x)\)がスコア関数の条件付き期待値の導関数なので、これを推定して逆行列をとる。例えば、\(\psi_{\theta}(O) = Y - \theta\) のようなシンプルな回帰モデルの場合には、\(V(x) = -1\) となり、\(x\)には依存しない。一方で、因果推論の場合には \(V(x) = - \begin{pmatrix}1 & E[T|X=x] \\ E[T|X=x] & E[T|X=x]\end{pmatrix}\)となる。よって、傾向スコアの推定が必要である（現在の grf では推定方法が異なるので必要はない）。

よって、分散の推定に必要なのは \(H_{n}(x; \theta,\nu)\) に対する推定量の構成ということになる。\(H_{n}(x; \theta,\nu)\)の推定のための基本的な方針は、ブートストラップを用いることである。記号を簡単にするために、 $$ \Psi(\hat\theta(x),\hat\nu(x)) = \sum_{i=1}^{n}\alpha_{i}(x)\psi_{\theta,\nu}(O_i) $$ とおく。次のような Half サンプルに基づく推定量を考える。 $$ \widehat{H}{n}^{HS} = \binom{n}{\lfloor n/2 \rfloor} \sum(\hat\theta(x),\hat\nu(x)) - \Psi(\hat\theta(x),\hat\nu(x)) \right)^2 $$ ここで、}:|\mathcal{H}|= \lfloor n/2 \rfloor} \left(\Psi_{\mathcal{H}\(\Psi_{\mathcal{H}}\) は Half サンプル \(\mathcal{H} \in \{1, 2, ..., n\}\) を回帰木の\(\mathcal{I}\)-サブサンプルが含むような回帰木によって推定される\(\Psi\)によって計算される。ここで、\(\mathcal{I} \subseteq \mathcal{H}\)であることを満たしているようなすべての回帰木を考えれば良いので、サブサンプルの\(\mathcal{J}\)の選び方については \(\{1,2,...,n\}\backslash \mathcal{I}\) の自由度があり、\(\mathcal{I} \subseteq \mathcal{H}\)であればよいから、\(\mathcal{I}\)についてはこれを含む形でより大きな集合をとっても良いので、ランダムフォレストを構成できる。このようなランダムフォレストの構成をすべての\(\mathcal{H}\)の全てのパターンについて考えるというものである。実際、この推定量は\(H_{n}(x; \theta,\nu)\)に対する一致性を持つことは、Athey et al.(2020) の Theorem 6からわかるが、すべての \(\mathcal{H}\) について計算するのは実質的に不可能である。

そこで、Sexton and Laake (2009) で提案された方法を拡張することで、\(\widehat{H}_{n}^{HS}\) に対する推定量を構成する。いま、\(\ell \geq 2\) を little bag size とし、\(B\) を\(\ell\)に乗じる定数とする。いま、\(g = 1, ..., B/\ell\) 回 \(|\mathcal{H}_{g}| = n/2\) となる Halfサンプル \(\mathcal{H}_{g} \in \{1,2,...,n\}\) をとる。 次に、回帰木を構成する際の予測に用いるサンプル \(\mathcal{I}_{b}\) を、$\mathcal{I}{b} \subseteq \mathcal{H} $ となるように \(b=1,2,...,B\) に対して取る。

例えば、\(\ell = 2\), \(B=10\) の場合を考えると、\(\mathcal{H}_{g}, (g=1,2,3,4,5)\) が Halfサンプルとして生成される。次に、ランダムフォレストの木を \(B\) 回生成するときに、それぞれの木を構成する際のDouble Sampleで予測に用いるサンプル $\mathcal{I}{b} = \mathcal{H} $ となるようにするので、\(\mathcal{I}_{1} = \mathcal{H}_{\lceil 1/2 \rceil} = \mathcal{H}_{1}\), \(\mathcal{I}_{2} = \mathcal{H}_{\lceil 2/2 \rceil} = \mathcal{H}_{1}\), \(\mathcal{I}_{3} = \mathcal{H}_{\lceil 3/2 \rceil} = \mathcal{H}_{2}...\) のようになる。

ここで、データが固定されたもとでのサブサンプリングメカニズムについての期待値を\(E_{ss}\)で表すことにすると、以下の形に分割できる。 $$ E_{ss}\left[\left(\frac{1}{\ell}\sum_{b=1}^{\ell}\Psi_{b} - \Psi\right)^{2}\right] = \widehat{H}{n}^{HS} + \frac{1}{\ell-1}E\right] $$ この展開は、左辺の展開において}\left[\frac{1}{\ell}\sum_{b=1}^{\ell}\left(\Psi_{b} - \frac{1}{\ell}\sum_{b=1}^{\ell}\Psi_{b}\right)^{2\(\Psi_{\mathcal{H}}\)を挟んで分解する。また、\(E_{ss}[\Psi_{\mathcal{H}} - \Psi]\)は、すべての Half サンプル\(\mathcal{H}\)について平均を取る操作を表すから、これは\(\widehat{H}_{n}^{HS}\)　そのものであることに注意する。左辺は、little in Bagsによって構成されるフォレストの分散（全分散)が、これが \(\widehat{H}_{n}^{HS}\) と、little in Bags内の分散に分解されるという意味になる。

以上の議論から、十分大きな\(B\)に対しては上記の関係性が成り立つから、漸近分散に対する一致推定量が計算可能となる。この方法による漸近分散の計算は、回帰木が最初に与えた仮定を満たしていれば妥当となるので、以降に説明する様々な一般化ランダムフォレストによる推定に対して漸近分散を与えることができる。

推定方程式に対応した勾配ベースの回帰木

ここまでの議論では、推定方程式によって定義されるパラメータを扱っているにもかかわらず回帰木の構成法に触れていない。これは、回帰木の構成法にランダムフォレストによって生成された重みに基づく推定方程式の解が一致性や漸近正規性を有しており、漸近分散の導出結果を用いてよいということを表している。しかし、一般的にランダムフォレストを構成する際に、回帰木の学習は損失関数の最小化することによって行われる。これは、実際その方が効率的な学習ができるということと、有効な変数以外でのノードの分割を避けることは、予測における実質的な次元削減の効果があることから有効である。ここでは、Athey et al.(2020) で提案された推定方程式によって定義されるパラメータの推定において有効かつ計算効率の良いノードの繰り返し分割法である勾配ベースの回帰木について導入を行う。

推定方程式によって定義されるパラメータは一般的に観測されず推定方程式を解くことによって得られる場合がほとんどである。そのため、予測の対象となる結果変数\(Y\)が観測されたもとでの、条件付き平均 \(\mu(x)=E[Y|X=x]\)の場合とは異なり、損失関数を工夫する必要がある。条件付き平均を予測する問題では、ノードを2つに分割する際には、分割された後のノード内での結果変数の分散の重み付け和が最小になるようにするのが一般的である。勾配ベースの回帰木においてもこの考え方を継承し、ノードを分割した後でパラメータ \(\theta(x)\)のノード内での分散の重み付け和が最小になるように分割を考える。

いま、ノード\(P\)を2つの子ノード\(C_1, C_2\)に分割する場合を考え、ノード\(P, C_1, C_2\)における推定方程式の解を $$ (\hat\theta_{P}, \hat\nu_{P}) \in \argmin_{\theta,\nu}\left{\left|\sum_{{i\in \mathcal{J}:X_i \in P}}\psi_{\theta,\nu}(O_i)\right|{2}\right} $$ $$ (\hat\theta(O_i)\right|}}, \hat\nu_{C_{j}}) \in \argmin_{\theta,\nu}\left{\left|\sum_{{i\in \mathcal{J}:X_i \in C_{j}}}\psi_{\theta,\nu{2}\right}, \qquad j=1,2 $$ とする。このとき、損失関数 \(\mathrm{err}(C_1, C_2)\) は $$ \mathrm{err}(C_1, C_2) = \sum - \theta(X))^2 \mid X\in C_j] $$ とかける。この損失関数を用いて回帰木を学習させるには2つの問題がある。1つ目は}\Pr(X\in C_j \mid X\in P) E[(\hat\theta_{C_{j}\(\theta(X)\)が観測できないことである。もう1つは、候補となる分割を試すたびに推定方程式を解く必要があるということである。特に、大量の回帰木を作成するランダムフォレストの場合、2つ目の問題はクリティカルな課題となる。そこで、漸近的には同値な言い換え (Proposition 1. of Athey et al.(2020))を用いて、 $$ \widetilde{\mathrm{err}}(C_1, C_2) = \frac{n_{C_1}n_{C_2}}{n_{P}^{2}}(\hat\theta_{C_1}-\hat\theta_{C_2})^{2} $$ とする。ここで、\(n_{P}, n_{C_1}, n_{C_2}\)はそれぞれノード\(P, C_1, C_2\)のノードに含まれるサンプルサイズである。この言い換えは、ノード内の重み付け分散の最小化問題を、ノード間のパラメータ推定量の差の大きさに言い換えたものであるから自然である。

次に2つ目の問題であった推定方程式の繰り返し計算の問題を、推定方程式の解を漸近展開することで親ノードの推定方程式の解と影響関数を用いて近似する。具体的には、 $$ \hat\theta_{C_j} \approx \hat\theta_{P} - \frac{1}{|{i:X_i\in C}|}\sum_{{i\in X_i \in C_j }}\xi^{T}A_{P}^{-1}\psi_{\hat\theta_{P},\hat\nu_{P}}(O_i) $$ とする。ここで、 \(A_{P}\)は\(\nabla E[\psi_{\hat\theta_{P},\hat\nu_{P}}(O_i)\mid X\in P]\) に対する一致推定量であればよく、スコア関数が微分可能な場合には、 $$ A_{P} = \frac{1}{|{i:X_i\in P}|}\sum_{{i\in X_i \in P }}\nabla\psi_{\hat\theta_{P},\hat\nu_{P}}(O_i) $$ によって推定する。この計算の中には、分割後のノード\(C_j\)における推定方程式の計算が行われていないことに注意する。よって、 $$ \rho_{i} = -\xi^{T}A_{P}^{-1}\psi_{\hat\theta_{P},\hat\nu_{P}}(O_i) $$ とおいて、 \(\widetilde{\mathrm{err}}(C_1, C_2)\) に代入したものを計算すると、 $$ \widehat{\mathrm{err}}(C_1, C_2) := \left(\sqrt{\frac{n_{C_2}}{n_{C_1}}}\sum_{{i\in X_i \in C_1 }}\rho_{i} - \sqrt{\frac{n_{C_1}}{n_{C_2}}}\sum_{{i\in X_i \in C_2 }}\rho_{i}\right)^{2} $$ となる。よって、一般化ランダムフォレストにおける推定では、近似的な損失関数を定義した回帰木をもとにして学習を行うことで、計算量を減らしつつランダムフォレストの重み \(\alpha_{i}(x)\) を推定する。何度も指摘しておくが、この方式で学習しないと一致性がないということではなく、パラメータ\(\theta(x)\)に対する回帰木を効率的に学習し、\(\alpha_{i}(x)\)がパラメータ\(\theta(x)\)の変動をうまく反映するように構成するためのアプローチの1つであると認識してもらいたい。

一般化ランダムフォレストの応用

一般化ランダムフォレストは、その一般的な理論的枠組みによって様々な回帰の問題を扱うことができる。ここでは、一般化ランダムフォレストで扱える代表的な回帰モデルについて紹介する。これらの方法の大枠では、1) 興味のあるパラメータを規定する局所推定方程式 \(\psi_{\theta,\nu}(O_i)\)を構成し、次に局所推定方程式を用いて Double sample gradient tree （以降、勾配回帰木） を学習させる。パラメータの値を推定したい点\(x\)に対して、学習済み勾配回帰木を用いて推定方程式に観測データ\(X_i\)の重み \(\alpha_{i}(x)\) を計算する。最後に、点\(x\)におけるパラメータの推定量 \(\hat\theta(x), \hat\nu(x)\)を\(\alpha_{i}(x)\)で重み付けた推定方程式を解くことで推定する。 $$ \hat\theta(x),\hat\nu(x) = \operatorname{argmin}{\theta,\nu} \left{\left| \sum \right} $$ よって、}^{n}\alpha_{i}(x)\psi_{\theta,\nu}(O_i)\right|_{2\(\hat\theta(x), \hat\nu(x)\) は正規方程式の解となる。この方針に従うことで、様々な推定方程式の解を求めることができる。

quantile regression forest と、一般化ランダムフォレストによる分位点推定

ここでは、確率変数の組 \((X_i, Y_i) \in \mathbb{R}^{p}\times \mathbb{R}\) (\(i=1,2,...,n\)) が観測されたもとでの、条件付きq-分位点関数 \(\theta_{Y|X}(q) = \inf\{y : F_{Y|X}(y) \geq q\}\) を推定するための方法である。この問題は、スコア関数として $$ \psi_{\theta}(Y_i) = q - \mathbb{1}({Y_i \leq \theta}) $$ を設定したときの解である。条件付き分位点に対するスコア関数自体は微分できないが、この期待値をとった関数が微分可能であればGRFを用いた際の理論的妥当性は保証される。しかしながら、この方法に基づいた条件付き分位点関数の推定量は、GRFによる quantile_regression 関数では、分位点関数が滑らかではない上に、また2つの分位点\(q_1 < q_2\)に対する条件付き分位点関数が交差しないなどの条件が満たされないため、ノンパラメトリックに分位点を推定するというタスクにおいては、qgam (Fasiolo and Wood) や、GP(Boukouvalas et al) などを用いることが望ましいと言える。

一般化ランダムフォレストによる条件付き因果効果の推定

一般化ランダムフォレストが最も注目を集めたのは、条件付き因果効果の推定をノンパラメトリックに行えるという点である。以下では、処置変数が2値の場合を扱うが、処置変数は連続変数でも仮定の拡張によって同様の議論が成り立つため、推定の枠組みは同じである。

確率変数の組 \((X, T, Y(0), Y(1)) \in \mathbb{R}^{p}\times \{0,1\} \times \mathbb{R} \times \mathbb{R}\)が、\(Y(0), Y(1) \perp\!\!\!\perp T \mid X\) を満たし、任意の\(x \in [0,1]^{p}\)に対して\(0 < \Pr(T \mid X = x) < 1\)、確率変数\(Y = TY(1)+ (1-T)Y(0)\)であるとする。確率変数の組 \((X_i, T_i, Y_i), i=1,2,...,n\) が、\((X, T, Y(0), Y(1))\)の従う分布から観測されたとき、\(\tau(x) = E[Y(1)-Y(0)\mid X=x]\)を推定するというのが因果効果推定の一般的な枠組みである。\(\mu(x) = E[Y(0)\mid X=x]\)とおくと、強く無視可能な割り付けから\(\mu(x), \tau(x)\)に対するスコア関数は次のように導かれる。 $$ \psi_{\theta(x),\nu(x)}(O) = \psi_{\mu(x),\tau(x)}(T,Y) = \left(Y - \mu(x) - T\cdot \tau(x)\right)\begin{pmatrix}1 \ T \end{pmatrix} $$ 推定方程式 \(E[\psi_{\mu(x),\tau(x)}(O_i)\mid X_{i}=x] = 0\) に対しては、正規方程式と解くことで解が求められるので、 $$ \theta(x) = \xi^{T}Var(W_i \mid X_i = x)^{-1}Cov[W_i, Y_i \mid X_i = x] $$ となる。よって、ランダムフォレストによる重み \(\alpha_{i}\) を用いて、それぞれの期待値に対するカーネル重み付け推定量を計算することで $$ \hat\theta(x) = \xi^{T}\left(\sum_{i=1}^{n}\alpha_{i}(x)(W_i - \bar{W}{\alpha})(W_i - \bar{W}})^{T}\right)^{-1} \sum_{i=1}^{n}\alpha_{i}(x)(W_i - \bar{W{\alpha})(Y_i - \bar{Y}) $$ ただし、\(\bar{W}_{\alpha} = \sum_{i=1}^{n}\alpha_{i}(x)W_{i}\) および \(\bar{Y}_{\alpha} = \sum_{i=1}^{n}\alpha_{i}(x)Y_{i}\) とかける。

Local Centering (R-Learner by Nie and Wager., 2021)

条件付き因果効果の推定と同様の設定で、異なるアプローチを取ることもできる。Nie and Wager (2021) では、Partial linear model (Robinson 1988) の考え方を用いて、\(\tau(x)\)を推定する方法を提案しており、この考え方を取り入れて\(\tau(x)\)を推定する方法について説明する。いま、\(m(x) = E[Y_{i} \mid X=x]\) および、\(e(x) = E[W_{i} \mid X=x]\)とおくと、 $$ m(x) = E[(1-T)Y(0) + TY(1) \mid X=x] = (1-e(x))\mu(x) + e(x)(\mu(x)+\tau(x)) $$ であるから、 $$ Y - m(X) = (T - e(X)) \tau(X) + \varepsilon $$ が成立する。ただし、\(T \perp\!\!\! \perp \varepsilon \mid X\) および、\(E[\varepsilon \mid X] = 0\) である。

よって、\(m(X)\)と\(e(X)\)が既知であるときは、スコア関数を $$ \psi_{\tau}(Y,T) = (Y - m(X)) - (T - e(X)) \tau(X) $$ とすることで、\(\tau(x)\)を推定できることがわかる。すなわち、正規方程式を考えて\(\tau(x)\)について解くと、任意の\(x\)に対する局所推定方程式の解は、 $$ \begin{align} \theta(x) &= \xi^{T}Var[(W_i - E[W_i \mid X_i = x]) \mid X_i = x]^{-1} \ &\times Cov[(W_i - E[W_i \mid X_i = x]),(Y_i - E[Y_i \mid X_i]) \mid X_i = x] \end{align} $$ となる。よって、条件付き因果効果の推定と同様に、ランダムフォレストの重み \(\alpha_{i}(x)\) で重み付けた empirical version を推定量として用いる。

しかし、\(m(x)\)と\(e(x)\)は一般的には未知であるため、推定する必要がある。ただし、\(m(x)\)と\(e(x)\)をランダムフォレストや、ブースティングなどを用いて推定する場合には、真の関数への収束レートが\(\sqrt{n}\)オーダーではないためバイアスを生じるため、推定にはcross-fitting (Chernozhukov et al., 2018) を用いることが推奨されている。

cross-fitting は、サンプル\(i\)に対する\(m(x)\)と\(e(x)\)の推定値\(\hat{m}(X_i)\)を計算する際に用いるモデル\(\hat{m}\)の推定に、\(i\)サンプルを用いないという方法である。例えば、Kennedy et al., (2020) では、観測データを3つのグループ \(\mathcal{I}, \mathcal{J}, \mathcal{K}\)に分け、\(\mathcal{I}\)に対するスコア関数を計算する際には、\(m\)の推定量\(\hat{m}\)を\(\mathcal{J}\)を用いて構成し、\(e\)の推定量\(\hat{e}\)を\(\mathcal{K}\)を用いて構成し、\(i \in \mathcal{I}\)に属するサンプルに対する\(m(X_i)\)と\(e(X_i)\)の推定量を\(\hat{m}(X_i)\)と\(\hat{e}(X_i)\)として構成する。

この方法については、推定モデルが\(i\)サンプルを用いていなければよいため、ここでは明示的に\(m\)と\(e\)の推定量はそれぞれ\(i\)サンプルを用いていないことを明示するため、\(\hat{m}^{(-i)}\)および、\(\hat{e}^{(-i)}\) と書くことにする。新たにスコア関数として、 $$ \tilde{\psi}_{\tau}(Y,T) = (Y - \hat{m}^{(-i)}(X)) - (T - \hat{e}^{(-i)}(X)) \tau(X) $$ を定義し、このもとで勾配回帰木を実行する。ここで、注意として、\(\hat\tau(x)\)の一致性は、ここでの \(\hat{m}, \hat{e}\) の一致性が直接影響するわけではなく、 \(\alpha_{i}(x)\) で重み付けたあとのスコア方程式の設計に依存することを述べておく。

一般化ランダムフォレストによるInstrumental variableを用いた因果効果の推定

処置変数に対する因果効果の推定において、\(T_i\)と誤差\(\varepsilon\)が条件付き独立ではない場合に、因果効果を推定する手法として操作変数による因果効果の推定がある。独立同一な確率変数の組 \((X_i, Y_i, T_i, Z_i) \in \mathcal{X}\times \mathbb{R} \times \{0,1\}\times \{0,1\}\) (\(i=1,2,...,n\)) が観測されたもとで、以下のモデルを考える。 $$ Y_i = \mu(X_i) + T_i \tau(X_i) + \varepsilon $$ ただし\(T_i\) と \(\varepsilon\) は\(X_i\)が与えられたもとで独立ではなく、任意の\(x\in \mathcal{X}\)において \(Cov[Z_i, W_i \mid X=x] \neq 0\)かつ、\(Z_i \perp\!\!\!\perp \varepsilon \mid X=x\) が成り立つと仮定する。このとき、\(\tau(x)\) は識別可能で、 $$ \tau(x) = \frac{Cov[Y_i, Z_i \mid X_i = x]}{Cov[T_i, Z_i \mid X_i = x]} $$ となる。このとき、\(\tau(x)\) はモデルに基づいたスコア関数 $$ \psi_{\theta(x),\nu(x)}(O) = \psi_{\mu(x),\tau(x)}(T,Z,Y) = \left(Y - \mu(x) - T\cdot \tau(x)\right)\begin{pmatrix}1 \ Z \end{pmatrix} $$ を用いて推定可能である。また、この推定方程式もやはり Local centering を行うことができて、\(E[Y_{i}|X=x], E[T_{i}|X=x], E[Z_{i}|X=x]\) に対する推定量 \(m^{(-i)}(X_i)\), \(e^{(-i)}(X_i)\), \(\eta^{(-i)}(X_i)\)を用いて、 $$ \tau(x) = \frac{Cov[Y_i - m^{(-i)}(X_i), Z_i -\eta^{(-i)}(X_i) \mid X_i = x]}{Cov[T_i - e^{(-i)}(X_i), Z_i - \eta^{(-i)}(X_i) \mid X_i = x]} $$ とできるので、ランダムフォレストによって構成された重み \(\alpha_{i}\) による重み付けで推定量を計算することができる。

Random Survival Forest と 一般化ランダムフォレストによる生存関数推定

Random survival forest (SRF, Ishwaran et al.,2008) は、ランダムフォレスト手法を生存分析の設定に拡張し、制限的な仮定を課すことなく、共変量と生存結果間の複雑な関係を捉えることができる柔軟なノンパラメトリックな手法である。SRFはランダムフォレストを、共変量条件付きの生存関数および累積ハザード関数へ拡張するため、回帰木におけるノード分割において、ログランク統計量を用いて検定統計量が最も大きくなるように分割を行っている。回帰木を生存時間へと拡張した survival trees (生存木) については、Hamad et al.,(2011) が詳しい。また、Zhou and McArdle (2015) では、　生存木 について応用事例などをまとめている。

いま、\(T^{*}\)を生存時間（イベントまでの時間）とし、\(C\)を右側打ち切り時間とする。観測されるイベント発生時刻を\(T\)とすると、 \(\(T = \min(T^{*}, C)\)\) と表すことができる。また、右側打ち切りに対する指示変数を、 \(\(\Delta = \mathbf{1}\{T^{*} \leq C\}\)\) とする。共変量ベクトルを\(X \in \mathbb{R}^{p}\)とする。ただし、共変量ベクトルは時間変化しないものとする。また、生存時間における一般的な仮定である、打ち切りに対する条件付き独立性 \(\(T \perp C | X\)\)を仮定する。いま、\((T^{*},C,X)\) が従う分布\(P\)からの独立同一なサイズ\(n\)の確率変数 \((T^{*}_{i}, C_{i}, X_i)\) に基づいて、\((T_{i}, \Delta_{i}, X_i)\)が観測されたとする。このとき、SRFに基づく生存時間解析では、条件付き生存関数 \(S(t|x) = P(T > t|X = x)\) や、累積ハザード関数 \(\Lambda(t|x) = -\log S(t|x)\) を推定することが目標である。SRFでは、通常のランダムフォレスト同様に、サイズ\(n\)の訓練データからサイズ\(s\)のサブサンプリングを\(B\)回行い、それぞれに対して 生存木 を当てはめた結果を、平均することで条件付き生存関数および、累積ハザード関数を推定する。

SRFを構成する 生存木 では、通常の回帰木とは異なりログランク統計量を用いる。ログランク統計量は、2つの生存曲線を比較し、差があるかどうかを検定するための用いる統計量である。ここでは、親ノードを2つの子ノードに分割するための 生存木 の損失関数について説明する。いま、親ノード\(L\)内のイベント発生時間を\(\{t_1, t_2, \ldots, t_M\}\)とし、各時点\(t_j\)について、

  • \(d_j = \sum_{i \in L}\mathbb{1}\{T_i = t_j , \Delta_{i} = 1\}\)：親ノードにおける時間\(t_j\)でのイベント数
  • \(d_{j,L}\)：左子ノードにおける時間\(t_j\)でのイベント数
  • \(N_j = \sum_{i \in L}\mathbb{1}\{T_i \geq t_j\}\)：親ノードにおける時間\(t_j\)でリスク下にある対象者数
  • \(N_{j,L}\)：左子ノードにおける時間\(t_j\)でリスク下にある対象者数

と定義する。このとき、ログランク統計量は次のように与えられる。

\[L = \frac{\left(\sum_{j=1}^M \left(d_{j,L} - N_{j,L} \cdot \frac{d_j}{N_j}\right)\right)^2}{\sum_{j=1}^M \frac{N_{j,L}}{N_j} \cdot \left(1 - \frac{N_{j,L}}{N_j}\right) \cdot \frac{N_j - d_j}{N_j - 1} \cdot d_j}\]

この統計量は、二つの子ノードにおけるハザード関数が同じであるという帰無仮説の下で、漸近的に自由度1の\(\chi^2\)分布に従い、この値が大きほど生存時間関数に差があることを意味する。よって生存木では、生存時間関数のノード間での異質性を捉えるために、ログランク統計量を最大にする分割を各ノードで選択し、再起的分割を行い木構造を推定する。木構造が推定されたあとで、生存木の葉それぞれにおいてカプラン・マイヤー推定量（KM推定量）と、ネルソン・アーレン推定量（AL推定量）を計算し、条件付き生存関数と、条件付き累積ハザード関数を推定する。生存関数に対するKM推定量は、生存木の各葉内のサンプルに対して、 $$ \hat{S}{\text{KM}}(t) = \prod\right) $$ によって計算される。また、累積ハザード関数に対するAL推定量を、各葉に対して、 $$ \hat{\Lambda}} \left(1 - \frac{d_j}{N_j{\text{NA}}(t) = \sum $$ によって推定し、AL推定量に基づいて生存関数の推定を行う。 $$ \hat{S}} \frac{d_j}{N_j{\text{NA}}(t) = \exp(-\hat{\Lambda}(t)) $$ サバイバルフォレストでは、点}\(x\)における生存関数または累積ハザード関数を、生存木の各葉で推定されるKM推定量または、AL推定量の平均によって推定する。

生存木に用いる分割基準の損失関数は、ログランク統計量に限る必要はなく、2つの生存関数の非類似度を測る指標であればよい。すなわち、損失関数の設計は研究者や実務家が、何を使うかを研究の文脈において決定することが望ましい。例えば、Li et al. (2015) において報告されているように、ログランク検定は比例ハザード性が成り立つ場合において検出力が高いが、生存時間の交差が認められる場合などでは検出力が低下するため、代替的な手法を用いることが推奨されている。この論文では、Qui-Shengの2段階手順や、適応的ネイマン平滑化検定などが、生存関数の交差においてその差を検出できたという報告がなされている。また、Ishwaran et al.(2014) によってGray検定などを用いることで条件付き累積発生関数を推定することで、競合リスクがある場合へ対応した SRF も提案されている。

一方で、一般化ランダムフォレストでは、生存木ではKM推定量や、AL推定量を計算せずに重みを計算し、その重み付けによってカーネル重み付け生存関数の推定を行う。いま、\(T_{b}, b=1,2,...,B\)を、生存木であるとし、点\(x\)を含む葉を\(L_{b}(x)\)とすると、点\(x\)を含むランダムフォレストの重みは、回帰木の場合同様に $$ \alpha_i(x) = \frac{1}{B} \sum_{b=1}^B \frac{\mathbf{1}{X_i \in L_b(x)}}{|L_b(x)|} $$ となる。よって、サバイバルランダムフォレストにおけるKM推定量は\(\alpha_{i}(x)\)による重み付け関数として得られる。 $$ \hat{S}{KM}^{(GRF)}(t|x) = \prod\right) $$ また、AL推定量も同様に得られる。 $$ \hat{S}^{(GRF))}} \left(1 - \frac{\sum_{i=1}^n \alpha_i(x) \cdot \mathbf{1}{T_i = t_j, \Delta_i = 1}}{\sum_{i=1}^n \alpha_i(x) \cdot \mathbf{1}{T_i \geq t_j}{\text{NA}}(t|x) = \exp\left(-\sum\right) $$ GRFによる生存時間の推定量に対する一致性については結果が証明されているわけではなく、今後の課題である。関連する文献をいくつか挙げると、} \frac{\sum_{i=1}^n \alpha_i(x) \cdot \mathbf{1}{T_i = t_j, \Delta_i = 1}}{\sum_{i=1}^n \alpha_i(x) \cdot \mathbf{1}{T_i \geq t_j}Dabrowska (1989) において、カーネル重み付け条件付きKM推定量については一様一致性の結果が得られている。また、Bladt and Furrer (2025)競合リスク下における 条件付きAalen-Johansen推定量に対する一致性と漸近正規性の結果も示されている。

一般化ランダムフォレストによる 条件付き causal survival effectの推定

処置を受けたグループ（処置群）と、処置を受けなかったグループ（対照群）を比較して、その生存関数が個体の特性によってどう異なるのかに着目するのが、Conditional causal survival effect の推定である。ここでは、Cui et al. (2023) を中心に説明を行う。

確率変数の組 \(\{X, T(0), T(1), C, W\} \in \mathcal{X} \times \mathbb{R}_{+} \times \mathbb{R}_{+} \times \mathbb{R}_{+} \times \{0,1\}\) が従う分布を\(P\)とする。 \(X\) は特徴量であり、\(T(0)\) と \(T(1)\) は潜在的な生存時間、 \(C\) は打ち切り時間、 \(W\) はバイナリーの処置変数である。

確率分布\(P\)に独立同一に従うサイズ\(n\)の確率変数の組を \(\{X_{i}, T_{i}(0), T_{i}(0), C_{i}, W_{i}\}_{i=1}^{n}\) とし、ここから観測データ \(\{X_{i}, T_{i}, C_{i}, W_{i}\}_{i=1}^{n}\) が得られるとする。ただし、\(T_i = W_i T_{i}(1) + (1-W_i) T_{i}(0)\) である。Cui et al.,(2023)では、生存時間\(T\)に対する変換 \(y\)を用いて生存時間における条件付き因果効果 \(\tau(x)\)を定義する。 $$ \tau(x) = E[y(T_{i}(1)) - y(T_{i}(0)) \mid X_{i} = x] $$ ここで、Cui et al.,(2023) では、\(y\)に対して 定数 \(0 < h < \infty\) に対して、\(y(t) = y(h), \; (t \geq h)\)を仮定して、制限付き平均生存時間を考える。打ち切りが存在しない理想的な状況においては、標準的な因果推論の仮定のもとで、\(\tau(x)\)の推定は Local centering を行なった causal forest と同様の枠組みで推定が可能であり、スコア関数 \(\psi^{(c)}\)を $$ \psi_{\tau}^{(c)}\bigl(X_i, y(T_i), W_i; \hat{e}, \hat{m}\bigr) = \bigl[W_i - \hat{e}^{(-i)}(X_i)\bigr]\bigl[y(T_i) - \hat{m}^{(-i)}(X_i) - \tau(W_i - \hat{e}(X_i))\bigr] $$ とした一般化ランダムフォレストによって\(\tau(x)\)が推定可能である。ただし、\(\hat{m}^{(-i)}\)および \(\hat{e}^{(-i)}\) は \(e(x) = E(W_i = 1 \mid X = x)\)および、\(m(x) = E[y(T_i)\mid X_i = x]\) に対するcross-fitting推定量である。しかし、生存時間の打ち切りの影響を含んで因果効果を推定する場合には、打ち切りの影響を制御するための拡張が必要である。ここでは、識別可能性を保証するために、一般的な因果推論の仮定のほかに打ち切りに対して2つの仮定を置く。 - 処置変数と共変量を条件づけたもとでの打ち切り変数と生存時間は独立性： \(\(T_i \perp\!\!\!\perp C_i \mid X_i , W_i\)\) - 任意の処置変数と共変量を条件づけたもとでの打ち切り発生は確率的である：ある\(0 < \eta_{C} \leq 1\) \(\(\Pr[C_{i} < h \mid X_i, W_i] \leq 1-\eta_{C}\)\)

これらの仮定のもとで、Cui et al.,(2023) では、\(\tau(x)\)に対する推定量として、打ち切り逆確率重み付け（Inverse Probability censoring weighting）推定量を提案している。打ち切り過程に対する確率に対するモデル \(S_{w}^{C}(s\mid x)\) を $$ S_{w}^{C}(s\mid x) = \Pr[C_i \geq s \mid W_i = w, X_i=x] $$ で定義する。非打ち切りに対応する確率変数を、\(\Delta_{i}^{h} = 1\{(T_{i} \wedge h) \leq C_i\}\)とすると、非打ち切りデータに打ち切り過程の確率の推定値の逆数で重み付けた推定方程式 $$ \sum_{{i : \Delta_{i}^{h}=1}}\frac{\alpha_{i}(x)}{\hat{S}{w}^{C}(s\mid x)}\psi\bigr) = 0 $$ を解くことで、}^{(c)}\bigl(X_i, y(T_i), W_i; \hat{e}, \hat{m\(\tau(x)\)が推定できる。ここで、\(\hat{S}_{w}^{C}(s\mid x)\)は、survival forest を用いて推定することができる。 Cui et al., (2023) では、IPCW推定量を用いる場合の問題点として、\(\hat{S}_{w}^{C}(s\mid x)\) の推定精度が、 \(\tau(x)\) の推定精度に大きな影響を与えることと、非打ち切りデータのみに絞った解析の問題として、観測データを捨てることによる情報量ロスの問題を挙げている。そこで、Tsiatis (2007) において提案された打ち切りに対して2重に頑健な推定量の構成の考え方を用いて、Robins et al.(1994) を拡張した観測データに基づく新たなスコア関数を提案した。

\(U_i = T_i \wedge C_i\; (i=1,2,...,n)\)とし、スコア関数\(\psi_{\tau}(X_i, y(U_i), U_i \wedge h, W_i, \Delta_{i}^{h})\)を以下で定義する。

$$ \begin{aligned} \psi_{\tau}(X_i, y(U_i), U_i \wedge h, W_i, \Delta_i^{h}) &= \frac{\Delta_i^{h}\psi_{\tau}^{(c)}(X_i, y(U_i), W_i)}{S_{W_i}^{C}(U_i \wedge h|X_i)} \ &\quad + \frac{(1 - \Delta_i^{h})\mathbb{E}\left[\psi_{\tau}^{(c)}(X_i, y(T_i), W_i)|T_i \wedge h > U_i \wedge h, W_i, X_i\right]}{S_{W_i}^{C}(U_i \wedge h|X_i)} \ &\quad - \int_{0}^{U_i \wedge h}\frac{\lambda_{W_i}^{C}(s|X_i)}{S_{W_i}^{C}(s|X_i)}\mathbb{E}\left[\psi_{\tau}^{(c)}(X_i, y(T_i), W_i)|T_i \wedge h > s, W_i, X_i\right]ds, \end{aligned} $$ ここで、\(\lambda_{w}^{C}(s\mid x) = -\frac{d}{ds}\log S_{w}^{C}(s\mid x)\) である。第1項は打ち切られなかったデータに対するIPCWに対応するスコア関数である。第2項は、打ち切りデータに対して、\(E[\cdot \mid T_i \wedge h > U_i \wedge h, W_i, X_i]\) の条件は観測データに基づいた条件付きスコア関数の推定量である。これは Robins et al.(1994) でも述べられている Augmented IPW の補完項に対応する。最後の項は、Tsiatis (2006) の理論に基づく補正項で、\(\frac{\lambda_{W_i}^{C}(s|X_i)}{S_{W_i}^{C}(s|X_i)}\) は時間\(s\)における条件付きの打ち切りハザードであり、特定時点での打ち切りが発生する条件付き確率に対応する。

実際の因果効果推定においては、\(\mathbb{E}\left[\psi_{\tau}^{(c)}(X_i, y(T_i), W_i)|T_i \wedge h > U_i \wedge h, W_i, X_i\right]\) および、\(\lambda_{W_i}^{C}(s|X_i) \(、\) S_{W_i}^{C}(s|X_i)\)をデータから推定したものをプラグインしたスコア関数を用いる。ここで、推定量の構成は cross-fitting を用いて行う。 $$ \begin{aligned} &\psi_{\tau}(X_i, y(U_i), U_i \wedge h, W_i, \Delta_i^{h}; \hat{e}, \hat{m}, \hat{\lambda}^{C}{w}, \hat{S}^{C}}, \hat{Q{w}) \ &= \left(\frac{\hat{Q}}(U_i \wedge h|X_i) + \Delta_i^{h}[y(U_i) - \hat{Q{W_i}(U_i \wedge h|X_i)] - \hat{m}(X_i) - \tau(W_i - \hat{e}(X_i))}{\hat{S}^{C}\right. \ &\quad - \left.\int_{0}^{U_i \wedge h}\frac{\hat{\lambda}^{C}}(U_i \wedge h|X_i){W_i}(s|X_i)}{\hat{S}^{C}(X_i)), \end{aligned} $$ ここで、}(s|X_i)}[\hat{Q}_{W_i}(s|X_i) - \hat{m}(X_i) - \tau(W_i - \hat{e}(X_i))]\,ds\right)(W_i - \hat{e\(\hat{Q}_{W_i}(U_i \wedge h|X_i)\)は、\(\mathbb{E}\left[\psi_{\tau}^{(c)}(X_i, y(T_i), W_i)|T_i \wedge h > U_i \wedge h, W_i, X_i\right]\)に対する推定量である。

Cui et al.(2023) では、 $$ \sum_{i=1}^{n}\alpha_{i}(x)\psi_{\tau}(X_i, y(U_i), U_i \wedge h, W_i, \Delta_i^{h}; \hat{e}, \hat{m}, \hat{\lambda}^{C}{w}, \hat{S}^{C}) = 0 $$ の解として推定される条件付き因果効果の推定量 }, \hat{Q}_{w\(\hat\tau(x)\) は、cross-fitting 推定量が \(\hat{e}\),\(\hat{m}\),\(\hat{Q}\),\(\hat{\lambda}_{w}^{C}\),\(\hat{S}_{w}^{C}\) が一様一致性を満たす場合には、\(\hat{\tau}(x)\)が一致性を持つことを示している。また、条件付き因果効果の推定と同様に、cross-fitting推定量の収束速度が適当な正則条件を満たす場合には、causal forest の推定量と同様に漸近正規性が成り立つ。また、漸近分散の推定方法についても一般化ランダムフォレストの枠組みで推定することがで可能である。ただし、Cui et al.(2023) では、推定された\(\hat\tau(x)\)は、しばしば複雑な関数となることと、漸近分散は非常に大きくなるため、解釈を行う際には適当な特徴量の部分集合に対して線形射影した結果に対して解釈を行うことを推奨している。

Causal survival forest を応用した事例としては、Cui et al.,(2023) においては、AIDS臨床試験グループプロトコル175（ACTG175）[Hammer et al., 1996]データへの応用を行なっており、このデータでは4つの治療群へ患者はランダムに割り付けられており、そのうち2つの治療群間の間での条件付き生存因果効果の推定を行っている。また、Sverdrup and Wager (2024) では、National Job Training Partnership Act（JTPA）に基づくランダム化比較試験（RCT）のデータを用いて、職業訓練プログラムへの参加資格が失業期間に与える影響を、打ち切り時間（2年）を考慮してRMSTベースの因果効果を推定している。

一般化ランダムフォレストの推定精度の改善

一般化ランダムフォレストによる解は、漸近的な性能は良いが、実際のデータに当て嵌めを行った場合、ランダムフォレスト同様の問題が生じる場合が多い。ここでは、ランダムフォレスト特有の問題を解消するための方法についてLocal Linear forest について説明する。

Locally Linear Forest

ランダムフォレストは、点\(x\)に対する推定量をその周辺のサンプルに重み付けして推定するという性質上、特徴空間\(\mathcal{X}\)の裾などの境界上において推定精度が低下するという問題がある。Local regression forest (LLF, Friedberg et al.2020) は、この問題に対処するために局所回帰の考え方を用いる。いま、確率変数の組 \((X_1, Y_1), ..., (X_n, Y_n) \in [0,1]^{p} \times \mathbb{R}\) が観測されたと仮定し、\(Y_i = \mu(X_i) + \varepsilon\) を満たすとする。このとき、点 \(x \in \mathcal{X}\) に対する \(\mu(x) = E[Y \mid X=x]\)に対するLLF推定量を、 $$ \begin{pmatrix}\hat\mu(x_0) \ \hat\theta(x_0)\end{pmatrix} = \argmin_{\mu,\theta}\left{\sum_{i=1}^{n}\alpha_{i}(x_0)(Y_i - \mu(x_0) - (X_i - x_0)\theta(x_0))^{2} + \lambda|\theta(x_0)|_{2}^{2} \right} $$ によって定義する。\(\hat\mu(x_0)\) は \(\mu(x_0)\)に対する推定量であり、\(\theta(x_0)\)は、\(X_i - x_0\) による局所的な変動を補正するための項である。また、リッジ罰則 \(\lambda\|\theta(x_0)\|_{2}^{2}\) は局所変動への過剰な当てはまりを抑制する項である。

いま、\(A\)を\(A_{i,i} = \alpha_{i}(x_0)\) である対角行列とし、\(J\)を \((d+1)\times (d+1)\) の対角行列とし、\(J_{1,1} = 0\) かつ \(J_{i,i} = 1 (i=2,3,...,d+1)\) とする（この行列は、リッジ回帰の罰則項を\(\mu(x_0)\)にのみかけないことに対応する）。さらに、\(\Delta\)を切片と\(x_0 = (x_{0,1}, x_{0,2},...,x_{0,p})\)に中心化した特徴量からなる行列で、\(\Delta_{i,1} = 1\) かつ \(\Delta_{i,j+1} = X_{i,j} - x_{0,j}\) (\(j=1,2,...,p\)) とすると、LLF推定量は $$ \begin{pmatrix}\hat\mu(x_0) \ \hat\theta(x_0)\end{pmatrix} = (\Delta^{T}A\Delta + \lambda J)^{-1} \Delta^{T}AY $$ と陽に表すことができる。

次に、LLFにおける分割基準について考える。LLFの目標は、局所的な変動を除いて\(\mu(x)\)を推定することである。そこで、分割においても局所的な変動を除いた上で、\(\mu(x)\)の変動を捉えるような分割を考えるのは自然である。親ノード \(P\) において、サイズ\(n_P\) の観測データ \((X_1,Y_1),...,(X_n,Y_n)\) が含まれるとする。2つのノード\(C_1, C_2\) への分割は、分割後のノードの \(Y\) の平均を、\(\bar{Y_1}, \bar{Y_2}\) とおくと、2乗損失の和 $$ \sum_{i : X_i \in C_1}\frac{#{i : X_i \in C_1}}{n_{P}}(Y_i - \bar{Y}1)^{2} + \sum}\frac{#{i : X_i \in C_2}}{n_{P}}(Y_i - \bar{Y{2})^{2} $$ を最小にするのが一般的である。これに対して、LLFでは局所的な変動を除いた\(\mu(x)\)への当てはまりをよくする。親ノードで、\(Y_i\) に対してリッジ回帰モデルを当てはめた推定値を $$ \hat{Y}} = \hat{\alpha{P} + X $$ とする。ここで、}^{T}\hat\beta_{P}, \quad \hat\beta_{P} = (X_{P}^{T}X_{P}+ \lambda J)^{-1}X_{P}^{T}Y_{P\(X_{P}\) と \(Y_{P}\) はそれぞれ、親ノードに含まれる \(X_{i}\)と\(Y_{i}\)からなる特徴量の行列、および結果変数のベクトルである。LLFにおけるノードの分割では局所変動を除いた斬差 $$ Y-\hat{Y}_{i} $$ に対してCARTによる分割を行う。

ここで、リッジ回帰における罰則\(\lambda\)は、大きい値を用いると\(\mu(x)\)の変動の学習も抑制するため、回帰木の学習時には小さい値を使うことが推奨されている。一方で、推定方程式を解く場合には、回帰木の学習時よりも大きい罰則の値を用いる。

この方法は、条件付き因果効果の推定にも応用することができる。因果効果の推定の場合には、回帰の問題の場合に加えて、因果効果の局所変動も捉えたい。そこで、次のL2罰則付きのリッジ回帰を用いる。 $$ \left{\hat\tau(x_0), \hat\theta_{\tau}(x_0), \hat{a}(x_0), \hat\theta_{a}(x_0) \right} = \argmin_{\tau,\theta}\left{ \sum_{i=1}^{n}\alpha_{i}(x_0)\left( Y_i - \hat{m}^{(-i)}(X_i) - a - (X_i - x_0)\theta_{a} - (\tau + \theta_{\tau}(X_i - x_0))(W_i - \hat{e}^{(-i)}(X_i)) \right)^2 + \lambda_{\tau}|\theta_{\tau}|{2}^{2} + \lambda \right} $$}|\theta_{a}|_{2}^{2

一般化ランダムフォレストのハイパーパラメータ選択

一般的な機械学習と同様に、ランダムフォレストの予測精度は、ハイパーパラメータの選択に依存している。ここでは、ランダムフォレストのパラメータチューニングの方法と、それによる予測性能の改善について説明する。まず、一般的な結果として、Probst et al.,(2019) において、ランダムフォレストのデフォルトのパラメータと、ハイパーパラメータチューニングを行った場合の最小2乗誤差やAUCを比較し、ハイパーパラメータのチューニングがパフォーマンス改善に寄与するかを、SVMやXgBoostなどと比較をすることで、ランダムフォレストの tunability を評価している。その結果、k-近傍法についでランダムフォレストのパラメータはデフォルトと、CV等でチューニングした結果が変わらないという結果が得られている。よって、他のブースティングやNN等と比較して、ハイパーパラメータの選択が予測結果に大きな影響を与えないという点は、評価できる手法であることがわかる。

しかし、ランダムフォレストにおいてパラメータチューニングが全く必要ないというわけではなく、いくつかのパラメータのチューニングは重要であることがわかっている。Scornet (2020)の指摘にもある通り、ランダムフォレストにおいて推定量に影響を与えるのは、min-node-size と subsample size そして、mtry である。

min-node-size は、ランダムフォレストの生成するカーネル \(\alpha_{i}(x)\) がどの範囲に重みを与えるのかに関わっている。また、min-node-size が小さいと、狭い範囲に重みがかかることになるので、データのノイズが大きい場合などは、推定値が不安定になる可能性がある。一方で広すぎると、関数の変動を捉えきれなくなる。

subsample size は、データを2つに分割する回帰木のアルゴリズム（Double sample）を用いている限りは、木構造の学習と予測に使われるサンプルがランダムに分割されるため、異なる回帰木間の相関は小さくなるため、\(n \times [0.9,1.0)\) 程度に大きくとっても原理上は問題にはならない。

mtryは、回帰木1本を構築する際に用いる特徴量の数で、package{ranger} などを見る限り一般的には \(\sqrt{p}\) が用いられる。しかし、一般化ランダムフォレストにおいては mtry = min{p, 20 + $\sqrt{p}$} が用いられる。木における各ノードの分割で使用する変数の選択を \(\min\{1+\mathrm{Poisson}(mtry),p\}\) で行い、ランダム性を持たせているため、他のランダムフォレストとは扱いが異なるが、基本的にデフォルトで問題はない。

また、ランダムフォレストを構成する回帰木の本数 \(B\) は、サンプルと同じオーダー以上が望ましいとされることと、\(B\)は大きい方が回帰木のランダムネスが平均化されるため結果が安定する。Mentch and Hooker (2016) では、\(B=500\)以上を選択しており、この他 package{grf} などでも \(B=2000\)をデフォルトにしていることを指摘しておく。

以上の議論から、一般化ランダムフォレストを用いる場合にはsubsample size がチューニングにおいて最も重要になるため、この変数についてはクロスバリデーションなどを通して、慎重に選択することが望ましい。

最後にチューニングとは直接的に関係はないが、一般化ランダムフォレストにおける Double sampleのメカニズムによって、多くの論文でのパフォーマンスの評価は、\(n=1000\)以上が多く比較的大きな規模のデータセットを対象としていることには注意していただきたい。特に因果推論などの文脈において、\(n=400\)程度で一般化ランダムフォレストを使っているケースが見受けられるが、提案時において\(n\)が小さい状況はあまり想定されているものではないと著者は考えている。

一般化ランダムフォレスト予測の解釈

一般化ランダムフォレストに対する予測の解釈については、実は理論が大きく深まっているわけではなく現在発展途上である。従来のランダムフォレストでは、主に Mead Decrease Impurity (MDI;平均不純度低下)または、Mean Decrease Accuracy (MDA; 平均精度低下) のいずれかを用いて、ランダムフォレストの変数の重要度が測られてきた。

MDIは、各回帰木における分割で低下した損失に、元のノードに含まれるサンプルサイズ（の全体に対する）比を乗じた値が、選択された変数の重要度となる。これを全ての木で平均することでMDIが計算される。一方のMDAは、観測された変数のパーミュテーションを用いて計算される統計量で、ある変数 \(j \in \{1,2,...,p\}\)をモデルから除いた場合の推定制度の低下を、観測データのうち変数\(j\)をパーミュテーションすることで代替的に観測することで、推定誤差の上昇量によって変数\(j\)の重要度を捉える方法である。しかし、一般化ランダムフォレストでの関心は、パラメータ\(\theta(x)\)に対する変数重要度となるため、従来のランダムフォレストの枠組みで推定することが難しい。ここでは、Benard and Josse (2023) による条件付き因果効果に対する変数重要度の推定法を紹介する。

先ほどの条件付き因果効果の例において、\(\tau(x)\)に対する変数重要度について考える。いま、条件付き因果効果が観測変数\(X = (X_1,...,X_p)\)のうち、\(\mathcal{S} \subset \{1,2,...,p\}\)に依存する場合を考える。結果変数\(Y\)は、 $$ Y = \mu(X) + T \cdot \tau(X^{\mathcal{S}}) + \varepsilon $$ を満たすとし、前に指摘したモデルの仮定が成り立つと仮定する。このとき、条件付き因果効果に対する 変数\(j\)の変数重要度 \(I^{(j)}\) を以下のように定義する。 $$ I^{(j)} = \frac{Var[\tau(X^{\mathcal{S}})] - Var[\tau(X^{\mathcal{S}})\mid X^{(-j)}]}{Var[\tau(X^{\mathcal{S}})]} = \frac{E[(\tau(X^{\mathcal{S}}) - E[\tau(X^{\mathcal{S}})\mid X^{(-j)}])^{2}]}{Var[\tau(X^{\mathcal{S}})]} $$ \(I^{(j)}\) は \(j \not\in \mathcal{S}\)の場合には、\(I^{(j)}=0\) となり、変数\(X_{j}\) が\(\tau(X^{\mathcal{S}})\) に対して大きな影響を持つならば、\(I^{(j)}\)が大きくなることがわかる。ここで、定義から \(j \in \mathcal{S}\)ならば、\(0 < I^{(j)} \leq 1\)である。

ここで、\(\tau(X^{\mathcal{S}})\)は Local centering による推定方程式の解であるから、以下の方程式を満たす。 $$ \tau(X^{\mathcal{S}}) \times \operatorname{Var}[(T-e(X))\mid X^{\mathcal{S}}] -\operatorname{Cov}[T-e(X),Y-m(X) \mid X^{\mathcal{S}}] = 0 $$ 一方で、\(X^{(-j)} = (X_1, ..., X_{j-1}, X_{j+1},...,X_{p})\)を条件づけた場合には、以下の関係式が成り立つ。 $$ \begin{align} &E[\tau(X^{\mathcal{S}})\mid X^{(-j)}] \times \operatorname{Var}[(T-e(X))\mid X^{(-j)}] - \operatorname{Cov}[T-e(X),Y-m(X) \mid X^{(-j)}]\ & \phantom{====} + \operatorname{Cov}[\tau(X^{\mathcal{S}}),e(X)(1-e(X)) \mid X^{(-j)}] = 0 \end{align} $$ この結果から、\(E[\tau(X^{\mathcal{S}})\mid X^{(-j)}]\)に対して、以下の関係式が得られる。 $$ E[\tau(X^{\mathcal{S}})\mid X^{(-j)}] = \frac{\operatorname{Cov}[T-e(X),Y-m(X) \mid X^{(-j)}]}{\operatorname{Var}[(T-e(X))\mid X^{(-j)}]} - \frac{\operatorname{Cov}[\tau(X^{\mathcal{S}}),e(X)(1-e(X)) \mid X^{(-j)}]}{\operatorname{Var}[(T-e(X))\mid X^{(-j)}]} $$ この第1項は、\(\tilde{Y}_{i} = Y_{i} - m(X_{i})\) および \(\tilde{T}_{i} = T_{i} - e(X_{i})\) と定義した時に、\((X_{i}^{(-j)}, \tilde{T}_{i}, \tilde{Y}_{i})\) (\(i=1,2,...,n\)) に対して、causal forest を学習させて得られる条件付き因果効果の推定値である。第2項は \(X_{j}\) を条件から外すことで生じるバイアスに対応する。

よって、このバイアス項に対応する推定量が得られれば、causal forest の推定量に対する変数重要度を推定することができる。\((X_{i}, \tilde{T}_{i}, \tilde{Y}_{i}), i=1,2,...,n\) 上で学習された causal forest の点 \(x\) に対する推定量を\(\hat{\tau}(x)\) とする。一方で、特徴量の集合から特徴量\(X_{j}\)を除いたデータ \((X_{i}^{(-j)}, \tilde{T}_{i}, \tilde{Y}_{i}), i=1,2,...,n\) 上で学習された causal forest が生成するフォレストの重みを \(\alpha'_{i}(x^{(-j)})\)とし、条件付き因果効果の推定量を\(\hat{\tau}^{(-j)}(x)\) とする。このとき、\(E[\tau(X^{\mathcal{S}})\mid X^{(-j)}]\)に対するバイアス修正済みの条件付き因果効果の推定量を、 $$ \hat{\theta}^{(-j)}(x) = \hat{\tau}^{(-j)}(x) - \frac{\sum_{i=1}^{n}\alpha'{i}(x^{(-j)})\tilde{T}}^{2}\hat{\tau}(X_i) - \bar{T{\alpha'}^{2}\bar{\tau}}}{\bar{T{\alpha'}^{2} - (\bar{T} $$ で定義する。ただし、})^{2}\(\bar{T}_{\alpha'}^{2} = \sum_{i=1}^{n}\alpha'_{i}(x^{(-j)})\tilde{T}_{i}^{2}\), \(\bar{T}_{\alpha'} = \sum_{i=1}^{n}\alpha'_{i}(x^{(-j)})\tilde{T}_{i}\) および、\(\bar{\tau}_{\alpha'} = \sum_{i=1}^{n}\alpha'_{i}(x^{(-j)}) \hat\tau(X_i)\) である。それぞれ、\(X^{(-j)}\)を除いた causal forest によって得られた重みを用いて構成されている。

変数重要度は、得られた推定量で\(I^{(j)}\)の各項を置き換えることで得られる。いま、観測データの組 \(D= \{(X_i, T_i, Y_i)\}_{i=1}^{n}\) が従う分布\(P\) から、新たにサイズ\(n\)の観測データとは独立に i.i.d. なサンプル \(D' = \{(X'_i, T'_i, Y'_i)\}_{i=1}^{n}\) が得られたする。変数重要度の推定量は、\(D'\)に基づいて以下のように構成する。 $$ \hat{I}^{(j)} = \frac{\sum_{i=1}^{n}[\hat\tau(X'i) - \hat{\theta}^{(-j)}(X'_i)]^{2}}{\sum[\hat\tau(X'}^{ni) - \bar{\tau}]} - \hat{I}^{(0)} $$ ここで、\(I^{(0)}\) は、\(\hat{I}^{(j)}\) に対するバイアスの修正項である。一般的に、ランダムフォレストでは木の構成にランダムネスが存在するため、変数重要度の推定ではその影響を取り除いている。 $$ \hat{I}^{(j)} = \frac{\sum[\hat\tau(X'}^{ni) - \hat{\theta}^{(-j)}(X'_i)]^{2}}{\sum $$}^{n}[\hat\tau(X'_i) - \bar{\tau}]} - \hat{I}^{(0)

Benard and Josse (2023) の結果は、一般化ランダムフォレストによる推定量から、変数重要度を計算する場合には、従来の回帰の文脈では生じないバイアスが生じる可能性がある。バイアスの生じ方は、推定方程式の形状に依存するため、例えば生存データの文脈では再解析を行う必要がある。また、変数重要度の推定精度は、サンプルサイズに依存する。Benard and Josse (2023) では、 \(n=3000\) を用いてシミュレーションを行っている。変数重要度は結果の解釈においては有効であるが、その一方でランダムフォレストの推定精度に依存することには注意が必要である。

また、変数重要度の推定で注意が必要なのは、\(X_{j}\)の変数重要度を評価する際に、相関の強い変数\(X_{k}\)が存在すると、\(X_{j}\)を取り除いても代替する情報を\(X_{k}\)が与えるので、変数重要度が小さくなるという問題がある。これに対する対処としては、ドメイン知識などで変数をグループ化することで、変数単一の重要度ではなく、グループ単位の変数重要度を評価するなどの対処を行うことが望ましい。

一般化ランダムフォレストに対する変数重要度の今後の方針として、特に重要なアプリケーションの1つである生存時間解析への拡張は今後行われていくと予想される。また近年、変数重要度に関する研究は、セミパラメトリック推論からの理論解析、Model-X knockoffs などの概念、Shapley Value などの分野で発展している。一般化ランダムフォレストの推定量を、これらの枠組みの中で解釈する取り組みが現在も行われている。

一般化ランダムフォレストの応用事例

医療分野

  • Suicide Risk Among Hospitalized Versus Discharged Deliberate Self-Harm Patients: Generalized Random Forest Analysis Using a Large Claims Data Set

自傷行為を行って救急部に受診した患者は、自殺リスクが非常に高いため治療（入院）を行うことが推奨されている。この論文では、治療が自殺リスクに与える影響が明らかにすることを目的としている。このデータでは、57,312人の患者データに対して、causal forestを当てはめることで、交絡要因を取り除いて治療の因果効果を推定している。また、E-valueを用いて未観測交絡に対する影響を評価し、結果に対しては慎重であるべきであることを付記している。

経済分野

  • Heterogeneous effects of Medicaid coverage on cardiovascular risk factors: secondary analysis of randomized controlled trial

この論文では、オレゴン健康保険実験のデータを用いて、Medicaid（低所得者向け公的医療保険）の心血管リスク因子（主に収縮期血圧とヘモグロビンA1c）の改善効果における個人差（ヘテロジニアス効果）を機械学習の手法で評価している。12,134人の参加者を対象として、抽選でMedicaid申請可能か（実際に加入したかどうかは、別の変数）どうかを振り分けており、ランダム化比較試験となっている。収縮機血圧と、ヘモグロビンA1cを目的変数として、Medicate加入によりこれらの健康指標が改善するかを確かめている。共変量としては、年齢、性別、人種・民族、教育水準、既往症（高血圧、糖尿病など）および過去の医療費（総医療費や救急部利用費用）が観測されている。このモデルを、抽選を操作変数\(Z\)、Medicaidへの加入を処置 \(T\)として、instrumental variable causal forest を当てはめて、条件付き局所平均処置効果（Conditional LATE）効果を推定している。この分析の結果では、全体としての効果は小さいが、一部のサブグループで改善が見られることが、CLATEの推定によって明らかになったと結果を報告している。

金融分野

  • Predicting Value at Risk for Cryptocurrencies With Generalized Random Forests

この論文は、暗号資産（クリプトカレンシー）のリスク管理において重要な指標であるVaR（Value at Risk）の予測に焦点を当てています。従来の手法（例えば、分位点回帰、GARCH系モデル、CAViaRモデルなど）は、暗号資産特有の極端な変動性やスパイクに対応するには柔軟性が不足していることが指摘されている。そこで、著者ら Generalized Random Forests (GRF)による分位点推定を活用している。その結果、局所的な大きな変動に対応は、GRFによる条件付き分位点の推定が優れていることを報告している。この結果は、Shiraishi, Nakamura and Shibuki (2024)によって報告されたGRFによる時系列分位点推定の一致性の論文でも同様の結果が報告されている。

マーケティング分野

  • GCF: Generalized Causal Forest for Heterogeneous Treatment Effects Estimation in Online Marketplace

この論文では、オンラインマーケットプレイス（特にライドシェアリングサービス）の文脈において、連続的な治療（例：価格や割引）の異質的な介入効果（HTE）を非パラメトリックに推定するための手法として、Generalized Causal Forest を用いることを提案している。ライドシェアリングなどのプラットフォームでは、需要と供給のバランスを取るために、各地域・時間帯で価格や割引といった連続的な介入が重要である。例えば、ある地域で割引を実施することで利用者のリクエストを促進し、供給過多の状況を解消しようとする一方、過度な割引は混雑やサービス品質の低下を招く可能性があるため、データから適切な割引率などを推定したいという需要がある。特に、混雑時・平常時や、時間帯で提供するディスカウントを使い分けたいという需要も存在する。そこで、Colangelo and Lee (2020) によって提案された連続処置変数に対する用量反応関数の DML推定を GRFに組み合わせることで推定を実現している。実際のデータで、学習したアルゴリズムと別のアルゴリズムの比較も行なっており、その結果として提案した手法による方法がA/Bテストの結果、改善をもたらすことを報告している。

因果効果推定のいくつかの重要なトピックス

一般化ランダムフォレストを用いた因果効果の推定においては、モデルの性質や推定の精度に影響を与える重要な課題が複数存在する。ここでは、現在の因果推論研究において特に重要な「共変量の外れ値と重複の仮定」、「結果変数の外れ値とロバスト推定」、および「共形予測」について論じる。

因果推論における特徴量空間の外れ値とoverlap仮定

因果効果の推定において最も重要な仮定の一つが "Overlap Assumption" または "Common Support Assumption" である。この仮定は、任意の共変量 \(X = x\) に対して、処置群と対照群の両方に十分なサンプルが存在することを意味している。 $$ 0 < \Pr(T=1 \mid X=x) < 1 \quad (x \in \mathcal{X}) $$ この仮定が満たされない場合、つまり特定の \(x\) の値に対して処置群または対照群のサンプルが極端に少ないもしくは存在しない場合、その領域での因果効果の推定は外挿に依存することになり、不安定かつ信頼性の低いものとなる。特に、機械学習の文脈における性能比較などでこの仮定が満たされないベンチマークデータなどで比較が行われているケースが存在する。特に、一般化ランダムフォレストによる推定においても、共変量空間の端や\(\mathcal{X}\)の外れ値領域（X-outlier）は問題となる。

因果推論における結果変数の外れ値の問題

機械学習を用いた因果推論においては、多くの場合\(Y\)の分布に対して意識されないことが多い。特に、一般化ランダムフォレストで因果推論を行うケースにおいては、観測データをそのまま訓練データとして用いている場合が多い。一方で、\(Y\)の分布の一部のデータが大きな値を取るなどの分布の形状は、回帰木の学習時の2乗損失の最小化を通して、分割に影響を与える。そのため、\(Y\)のスケーリングは機械学習を用いた因果推論においては本質的な問題となる。医療データにおける医療費の分布や、マーケティングのデータにおける購入頻度の分布などは、ロングテールな分布であり、これらを直接機械学習に入力する場合には問題が生じる。そのためには、結果変数\(Y\)に対して変数変換を施して、機械学習に適した分布形状に変換し、損失関数が外れ値の影響を大きく受けすぎないような工夫をすることが重要となる。例えば、Bayesian additive regression trees (Chipman et al.,2010) においては、最大最小スケーリングを用いることを推奨している。ランダムフォレストにおいても、遺伝子データの解析などではBox-Cox変換が利用されることもあるため、一般化ランダムフォレストにおいても同様に変数変換については今後考えていく必要がある。または、\(Y\) の分布に関する問題については、この他の対処として損失関数をパラメータの2乗損失から、Huber損失や、擬似Huber損失などを採用することにより、より頑健な推定を行うという方針も考えられる。

Conformal Predictions

因果効果の推定において、点推定値だけでなく予測区間を提供することは、不確実性の評価と意思決定の観点から極めて重要である。特に、漸近論に基づかず予測区間が求めるカバレッジを達成することに関心が高まっている。 Gammerman et al.,1998によって提案されたConformal Prediction (CP) は、モデルに依存しない（model-agnostic）予測区間を構築するための手法であり、限られた仮定のもとで有限サンプルにおける正確な被覆確率（coverage probability）の下限を保証することができる。Papadopoulos et al.,2002 によって、CPの効率的な計算方法である Split Conformal Prediction が導入され、その後 Lei et al.,(2013)、 Lei and Wasserman (2014) および Lei et al (2018) によって体系的に整理された。Conformal prediction についての理論的な解説は Angelopoulos et al.(2024) が詳しい。その後、Tibshirani et al.,(2019) によって共変量シフトのもとでのConformal prediction の理論が提案されている。Lei and Candes (2021) は、共変量シフトのもとでのCPの枠組みを因果推論へと拡張し、Causal Conformal Inference を提案している。また、Alaa et al.,(2023) は、2値の処置における Meta-learner に対するCPの枠組みを提案している。Schröder et al.,(2024) では、連続処置に対するCPの枠組みを提案している。平均処置効果に着目したConformal Predictionとしては、Wang et al.,(2024) がある。

因果推論におけるCPの応用は、特に政策決定や臨床意思決定など、不確実性の適切な評価が重要な文脈で価値が高い。例えば、医療介入の効果に関する予測区間を提供することで、患者ごとの治療の不確実性を定量化し、より情報に基づいた意思決定が可能になるなどの利点がある。現在の一般化ランダムフォレストの推定は漸近論に基づいた結果となっているが、今後は conformal prediction等を通した、有限標本下でのカバレッジ保証を行なっていくことも求められている。

まとめ

本レビュー論文では、一般化ランダムフォレストとその因果推論への応用について体系的に概観した。Breiman (2001) によって提案された当初のランダムフォレストは、その予測精度と汎用性の高さから多くの分野で広く利用されてきたが、理論的な解析は困難とされてきた。2010年代以降、Lin and Jean (2006) のカーネル重み付け推定量としての解釈を基礎として、Biau (2012) や Scornet et al. (2015) らによる一致性の証明、Wager and Athey (2019) による漸近正規性の証明と漸近分散の推定方法の確立など、理論的な基盤が徐々に整備されてきた。

Athey, Tibshirani, and Wager (2020) によって提案された一般化ランダムフォレストは、ランダムフォレストをカーネル重み付け関数の推定法として捉え直し、局所推定方程式の枠組みに拡張することで、条件付き平均以外の多様なパラメータの推定を可能にした画期的な手法である。本手法の最大の利点は、条件付き因果効果、分位点推定、操作変数法による因果パラメータの推定など、従来のランダムフォレストでは困難だった統計量の柔軟な推定を理論的保証付きで行える点にある。

特に因果推論の分野では、Wager and Athey (2019) によるcausal forestの提案以降、様々な拡張が行われている。Local centering (Nie and Wager, 2021) によるアプローチは推定の効率性を高め、Local Linear Forest (Friedberg et al., 2020) は境界付近での推定精度の問題を改善し、推定の局所的な変動への適応力を向上させた。これらの技術的進展は、医療、経済、金融、マーケティングなど様々な分野での実証研究を可能にしている。

しかしながら、一般化ランダムフォレストを用いた推定にはいくつかの重要な課題も残されている。特徴空間における重複（overlap）の問題、結果変数の外れ値の扱い、サンプルサイズの小さい状況での適用限界などは、実践上の重要な課題である。また、理論的には共変量の次元が大きくなるにつれて収束のレートが悪化するという「次元の呪い」の問題も完全には解決されていない。

さらに、推定結果の解釈においても課題が残されている。Benard and Josse (2023) による条件付き因果効果に対する変数重要度の推定法の提案など進展はあるものの、一般化ランダムフォレストの枠組みにおける変数重要度の理論は発展途上である。また、予測の不確実性の定量化において、有限標本下でのカバレッジを保証するConformal Predictionの枠組みとの統合も今後の重要な研究方向である。

近年では、因果推論とランダムフォレストの両方の分野で理論的な発展が著しく、一般化ランダムフォレストはその二つの領域を橋渡しする重要な手法となっている。今後は、高次元データへの適応、モデル解釈性の向上、有限標本での性能保証、計算効率の改善など多くの方向性が考えられる。特に、複雑な因果構造を持つデータに対する拡張や、時系列データや空間データなど特殊な構造を持つデータへの応用は重要な研究課題である。

一般化ランダムフォレストは、理論と実践の両面で重要な貢献をもたらした手法であり、統計的機械学習と因果推論の融合を促進している。今後もデータ駆動型の意思決定が重要性を増す様々な分野において、その応用範囲はさらに拡大していくと考えられる。理論的基盤の深化と実践的な適用範囲の拡大という両面からの発展を通じて、一般化ランダムフォレストは因果推論の方法論的ツールキットの中で今後も重要な役割を果たすと考えられる。本レビュー論文が理論研究者と実務家の双方にとって一般化ランダムフォレストについて有益な橋渡しとなることを願っています。