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

Lecture on Tokyo Medical University titled "theory of generalized random forest and its application to causal inference."

Transcript

1. 一般化ランダムフォレストの理論と 統計的因果推論への応用 中村　知繁 January 27, 2024 順天堂大学健康データサイエンス学部
2. Table of contents 1. イントロダクション 2. 一般化ランダムフォレスト (GRF) 3. GRF 推定量に対する漸近理論 4. Tree の分割ルールの構成 5. 一般化ランダムフォレストの 統計的因果推論への応用 6. サンプルサイズが小さい場合の推論について 7. 今後の展望 1/52
3. イントロダクション
4. Breiman [2001] で提案されたランダムフォレスト（オリジナル） Breiman [2001] で提案されたオリジナ ルのランダムフォレストは、次の手順 で学習するものを指す。 • サイズ n の訓練データから、サイ ズ s (< N) のブートストラップサ ンプルを B 回とる。 • ブートストラップサンプル b = 1, 2, ...B に対して、CART を当 てはめる。 • 新たな入力 x に対して、B 本の Tree それぞれの結果を平均し、 random forest の推定量として返す このとき、Tree 間の相関をさげるた め、 CART の当てはめにおいては、 分割 に使用する変数を確率的に選択する。 Fig. 1: Breiman’s RF 2/52
5. ランダムフォレストの理論発展の概略 • ランダムフォレストは、Breiman [2001] において提案されて以来、回帰や判別 の様々なタスクへ応用されてきた。 • Breiman [2001, 2004] においては、条件付き期待値に対するランダムフォレスト 予測値の汎化誤差上界が評価されたが、推定量の漸近論を構成するには至って いない。 • 理論発展において重要な一歩は、Lin and Joen [2006] と、Biau and Devroye [2010] における適応的最近傍法と、ランダムフォレストの関係性を繋げ、汎化 誤差の下界を評価したこと。 • Wager and Walther [2014], Mentch and Hooker [2014], Scornet et al. [2015] では、 adaptive CART の一致性や、それぞれ別のアプローチで random forest の一致性 が示された。 • Mentch and Hooker [2016] がランダムフォレストの漸近正規性を証明（U 統計 量） 、Wager and Athey [2018] もバイアス評価、一致性、漸近正規性、漸近分散 の計算法を提案。 3/52
6. CART やランダムフォレストの条件付き平均以外への応用 randomforests の研究は、提案されて間もないころから、分位点推定や生存時間推 定へ応用されてきた歴史がある。 • Quantile regression forest [Meinshausen, 2006] • Bagging survival trees [Hothorn, 2004] • Random survival forests [Ishwaran et al., 2008] この背景には、CART の生存時間への応用研究がある。 • CART を生存時間へ応用するために、リスク関数や Goodness-of-fit を設計、ま たアルゴリズムの提案を行った [LeBlance and Crowley, 1992,1993] CART を用いたモデル分割アプローチ • Model-based recursive partitioning [Zeileis, Hothorn and Hornik, 2008] : Generalized M-fluctuation test [Zeileis and Hornik, 2007] を用いてモデルを特定 の共変量で分割する方が当てはまりが良くなるかを検定することで、モデルを 再帰的に分割する 4/52
7. 近年のランダムフォレストの研究 • 条件付き平均の推定（一般的なランダムフォレスト）: — Breiman (2001), Biau et al. (2008), Biau (2012), Wager and Walther (2015), Scornet (2016), Davis and Nielsen (2020) • 統計的因果推論: — Athey and Imbens (2016), Wager and Athey (2018) • 条件付き分位点推定: — Meinshausen (2006), Shiraishi, Nakamura and Shibuki (2024) • 生存時間解析: — Hothorn et al.(2006), Ishwaran et al. (2008), Zhu and Kosorok (2012), Steingrimsson et al. (2016), Steingrimsson et al. (2019). • 局所推定方程式（一般化ランダムフォレスト）: — Athey, Tibshirani and Wager [2019], Friedberg et al. [2020], Oprescu et al. [2019] • Causal survival forest: — Cui et al. [2023] 5/52
8. 今日の発表内容について 今日の発表では、次の 2 つについて発表を行います。 1. 一般化ランダムフォレストの理論的な理解 ランダムフォレストは、 （同質な個体に対する）重みづけカーネルを推定するた めの手法であるという視点の変化が非常に重要で、これが理解できるとあとは 直感的にカーネル重みづけ回帰と話が結び付きます。 2. causal survival forest 1. の理解に基づけば、ターゲットとなるパラメータに対する局所推定方程式が 立式できればよいということになります。その例の 1 つとして、Cui et al. [2023] の causal survival forest を扱います。 6/52
9. 一般化ランダムフォレスト (GRF)
10. GRF model Definition 1. (GRF model) 分布 F に従う確率変数の組 (X, O) ∈ X × O が次の局所推定方程式を満たすと仮 定する。 E ψθ(x),ν(x) (O) | X = x = 0 for all x ∈ X (1) ただし、 • θ = (θ(x))x∈X ∈ Θ = θ : X → Rkp : 興味のある（汎関数）パラメータ • ν = (ν(x))x∈X ∈ Θ = ν : X → Rkq : 局外（汎関数）パラメータ • ψ·,· (·) : X × O → Rkψ : スコア関数 (Note) スコア関数 ψ の例 • 条件付き平均: ψθ(x) (y) = y − θ(x) • 条件付き τ ∈ (0, 1) 分位点: ψθ(x) (y) = τ − 1{y≤θ(x)} • 密度関数 f のもとでの局所最尤推定量: ψθ(x) (y) = ∇ log(f θ(x) (y)) 7/52
11. GRF estimator Definition 2. (GRF estimator) 確率変数の組 Dn := {(Xi, Oi)}i=1,...,n がそれぞれ独立同一に分布 F に従うと仮定 する。このとき、式 (1) で定義されるパラメータ θ ∈ Θ に対する GRF 推定量を以 下で定義する。 (ˆ θ(x), ˆ ν(x)) ∈ argmin θ(x),ν(x) n i=1 αGRF i (x)ψθ(x),ν(x) (Oi) for all x ∈ X ただし、 • αGRF i (x) ∈ [0, 1] : ランダムフォレストを用いて計算される重み関数 αGRF i (x) = 1 B B b=1 αGRF bi (x), αGRF bi (x) = 1{Xi∈Lb(x)} |Lb(x)| • B : フォレストを構成する木の本数 • Lb(x) : テスト点 x ∈ X を含む b 番目の木の葉（または終端ノード） • |Lb(x)| : 終端ノード Lb(x) に属する訓練データの数 8/52
12. Forest Weights の構成の図解 Fig. 2: 上段では木 b = 1, 2, 3, ... において点 x と同じ終端ノードを共有する標本 i の重みが αbi(x) = 1 となる点を黒く塗っている。下段の図では各 i について αbi(x) の和の大きさが点 の大きさである。(Athey et al., [2019] より引用) 9/52
13. αi(x) とカーネル重みづけ関数のつながり n 組の観測 (Xi, Yi), i = 1, 2, ..., n が観測された場合に、θ(x) = E[Yi | Xi = x] を CART で推定する問題を例に、この 2 つのつながりを説明します。 • 回帰樹の推定量は、テスト点 x が与えられたもとで、x に距離的に近い訓練デ ータを用いて推定量を計算している。 ⇒「x に近い訓練データに重みづけるような局所カーネルを用いて、カーネル 重みづけ平均を計算している」 • 回帰樹の学習アルゴリズムである再帰的分割では、Tree の θ(x) への当てはまり が良くなる（2 乗誤差を最適化する）ように分割が行われる。 ⇒「θ(x) の勾配が大きい方向の葉の長さは短くなり、θ(x) の勾配が大きい方向 の葉の長さは長くなる。 つまり、αGRF bi (x) を足し合わせることで生成される GRF-weights αGRF i (x) は、点 x におけるターゲット θ(x) の変化が緩やかな方向へ長く、一方で大きく変化する方 向に短い重みをかけるような関数となる。 10/52
14. 例: αGRF i (x) の例 0 2 4 6 8 10 12 14 16 18 20 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 X1 X2 Fig. 3: 図: forest-based kernel による重みづけ 11/52
15. GRF 推定量に対する漸近理論
16. イントロダクション 一般化ランダムフォレスト (GRF) GRF 推定量に対する漸近理論 Tree の分割ルールの構成 一般化ランダムフォレストの 統計的因果推論への応用 サンプルサイズが小さい場合の推論について 今後の展望 12/52
17. 漸近理論の構成 • 以上の議論から、直観的にランダムフォレストの重みづけ関数は、カーネル重 みづけ関数と同様に議論できそう。 • スコア関数と、random forest を構成する tree に制約を付けることで、 GRF-weights の挙動を制御すれば、GRF 推定量は漸近的に良い性質を持つ。 議論における注意。 • 推定方程式の解である θ(x) について、特徴空間 X を再帰的に分割する方法を 議論していないが、推定量が漸近的に正規性を持つために特に必要な条件では ない。 • 有限の n に対する推定精度においては、再帰的分割の損失関数の設計は極めて 重要な役割を持つ（理論のあとに議論を展開する） •（明らかになっていないこととして）漸近有効性と損失関数の設計については、 現状十分な議論ができるほど明らかにはなっていない。 13/52
18. Tree の再帰的分割に対する制約（Specification） Specification • ω-Regular: 親ノード P を 2 つの排反なノード C1 , C2 に分割する際、親ノード に含まれる観測の 100ω% 以上が、子ノード C1 , C2 に含まれる。すなわち、 nP , nC1 , nC2 をノード P, C1 , C2 のサンプルサイズとするとき、以下が成り立つ。 min(nC1 , nC2 ) ≥ ω × nP • Random Split: 任意の分割において、すべての特徴量が選択される確率が π > 0 で抑えることができる。すなわち、 P(ξ = j) ≥ π > 0 for all j ∈ {1, . . . , p} • PNN (Potential Nearest Neighbor) k-set: 解析者が設定した k ∈ N に対して、N を木の終端ノードの集合とするとき、任意のノード L ∈ N が、k ≤ |L| ≤ 2k − 1 を満たす（葉に入っているサンプル数が k 以上 2k − 1 以下になるまで再帰的分 割を行う） 。 • Subsample size: ランダムフォレストを構成する木は、サイズ n の観測データ からサイズ s のサブサンプルをとったもので学習する。このとき、β ∈ (0, 1) に 対してサブサンプルサイズ s が、s = nβ を満たすようにとる。 14/52
19. Double sample tree algorithm Forest を構成する Tree の学習は Double sample tree algorithm [Athey and Imbens, 2016] を用いて行う。 Training data (𝐷𝐷𝑛𝑛 ) Test data Subsample (𝐷𝐷𝑠𝑠 (1)) 𝐽𝐽-sample (𝐷𝐷𝑠𝑠 𝑏𝑏)𝐽𝐽 Subsample (𝐷𝐷𝑠𝑠 (𝑏𝑏)) Subsample (𝐷𝐷𝑠𝑠 (𝐵𝐵)) ・・・ ・・・ 𝐼𝐼-sample (𝐷𝐷𝑠𝑠 𝑏𝑏)𝐼𝐼 (i) Randomly choose subsample with size s = 𝑠𝑠(𝑛𝑛) without replacement (ii) Randomly divide into 𝐽𝐽-sample and 𝐼𝐼- sample with size 𝑠𝑠/2 (iii) Create regression tree based on 𝐽𝐽-sample (iv) Create partitioning of space 𝜒𝜒 based on (iii) (v) Categorize 𝐼𝐼-sample falling in each rectangle (leaf) based on (iv) Fig. 4: Illustration of double sample tree algorithm 15/52
20. スコア関数に対する仮定 スコア関数に対しては、以下の 6 つの仮定をおく。ここで、スコア関数 ψθ,ν (Oi) に対して、スコア関数の条件付き期待値を Mθ,ν (x) := E[ψθ,ν (Oi) | Xi = x] で定義 する。 Assumptions (A1-A3) • Lipschitz x-signal : 任意の (θ, ν) に対して、Mθ,ν (x) は x ∈ X についてリプシッツ連続である. • Smooth identification : Mθ,ν は任意の (θ, ν) において 2 階連続微分可能であり、一様有界な 2 次導関数 を持つ。また、任意の x ∈ X において V(x) := V θ(x),ν(x) (x) は正則である。ただ し、Vθ,ν (x) := ∂/∂(θ, ν)Mθ,ν (x) |θ(x),ν(x) である。 • Lipschitz (θ, ν)-variogram: スコア関数 ψθ,ν (Oi) は (θ, ν) に対して連続な共分散関数を持つ。 sup x∈X ∥Var [ψθ,ν (Oi) − ψθ′,ν′ (Oi)|X = x]∥F ≤ L θ ν − θ′ ν′ 2 16/52
21. スコア関数に対する仮定（続き） Assumptions (A4-A6) • Regularity of ψ: スコア関数 ψθ,ν が次の形で分解可能である。ψθ,ν = λ(θ, ν; Oi) + ζθ,ν (g(Oi)) た だし、λ(θ, ν; Oi) は (θ, ν) についてリプシッツ連続な関数、g : {Oi} → R は Oi を 1 次元へ射影する関数であり、ζθ,ν : R → R は単調かつ有界な関数の任意の 族である。 • Existence of solutions: n i=1 αi = 1 を満たす任意の重み αi に対して、局所重みづけ推定方程式の解 (ˆ θ, ˆ ν) が存在し、ある定数 C ≥ 0 に対して ∥ n i=1 αiψˆ θ,ˆ ν (Oi)∥2 ≤ C max αi に 解が存在する。 • Convexity: スコア関数 ψθ,ν (Oi) が (θ, ν) に関する凸関数の負の劣勾配であり、スコア関数 の期待値 Mθ,ν (Oi) が (θ, ν) に関する狭義凸関数の負の勾配である。 17/52
22. GRF 推定量に対する漸近正規性 Athey et al. [2019] 定理 5 仮定 (A1)-(A6) および、Specification のもとで、Tree の学習がサブサンプルサイズ s = nβ で βmin := 1 − 1 + π−1 log(ω−1) / log (1 − ω)−1 −1 < β < 1 を満たすとする。 さらに、ρ∗ i (x) := −ξTV(x)−1ψθ(x),ν(x) (Oi) が、Var[ρ∗ i (x) | Xi = x] > 0 を満たすと する。 このとき、σ2 n = polylog(n/s)−1 を満たす列 {σn} が存在して、 (ˆ θn(x) − θ(x))/σn(x) d −→ N(0, 1) が成り立つ。ただし、polylog(n/s)−1 は 0 より大きい値をとる関数で、サンプリ ング比 n/s = n1−β の逆対数について高々多項式の速さで増加する関数である。 18/52
23. デルタ法による信頼区間の構成 (1) Athey et al. [2019] の定理 5 の結果から、ˆ σn(x)/σn(x) p → 1 を満たす ˆ σn(x) を ˆ θ(x) から構成できれば、Slutsky の補題から、 lim n→∞ E θ(x) ∈ (ˆ θn(x) ± Φ−1 1 − α 2 ˆ σn(x)) = α を示すことができる。 ここで、Athey et al.[2019] 定理 5 の証明から Var[˜ θ∗(x)] = Var[θ(x) + n i=1 αi(x)ρ∗ i (x)] (2) = ξTV(x)−1Hn(x; θ(x), ν(x))V(x)−T ξ (3) に対して、Var[˜ θ∗(x)]/ˆ σ2 n(x) → 1 が成り立つことがわかる。ただし、 Hn(x; θ, ν) = Var n i=1 αi(x)ψθ.ν (Oi) (4) である。 19/52
24. デルタ法による信頼区間の構成 (2) V(x) に対する推定量は、スコア関数の勾配の条件付き期待値なので推定すること は可能なので、Hn(x; θ(x), ν(x)) に対する推定量のみ考える。ここでは、Efron [1982] の Bootstrap に関する議論から、次の Half-sampling 推定量 ˆ HHS n (x) = n ⌊n/2⌋ −1 {H:|H|=⌊n/2⌋} ΨH ˆ θ(x), ˆ ν(x) − Ψ ˆ θ(x), ˆ ν(x) に対して、次の定理が成り立つ。 Athey et al. [2019] 定理 6 定理 5 の条件もとで、ˆ HHS n (x) は、∥ˆ HHS n (x) − Hn(x; θ, ν)∥F/∥Hn(x; θ, ν)∥F p → 0 が 成り立つ。さらに、∥ˆ V(x) − V(x)∥F p → 0 を満たす V(x) に対する一致推定量 ˆ V(x) が与えられたもとで、 ˆ σ2 n = ξT ˆ Vn(x)−1 ˆ HHS n (x)ˆ V(x)−T ξ は、ˆ σ2 n/σ2 n → 1 を満たす。 20/52
25. デルタ法による信頼区間の構成 (3) • Half-sampling estimator とは、観測されたデータのうち半分のデータを用いて 計算される統計量であるがすべての Half サンプルに対して和を取ることは実質 的に不可能なので、これに対して Athey et al. [2019] では Sexton and Laake [2009] では近似的な推定量を用いている • ランダムフォレストにおいては bagging を行った際にサンプル (Xi, Oi) を含ま ない tree は約 37% 存在するため、それらを用いて簡易的に random forest を構 成することができる。この考え方を応用することで、新たにランダムフォレス トを計算せずに GRF に対する分散推定量を得ることができる。 • この場合の推定量は ˆ HHS n (x) に対して不偏ではないので、バイアスの補正が必 要になる。 • バイアス補正の方法については、Sexton and Laake [2009] の jacknife after bagging が参考になる。 21/52
26. Tree の分割ルールの構成
27. イントロダクション 一般化ランダムフォレスト (GRF) GRF 推定量に対する漸近理論 Tree の分割ルールの構成 一般化ランダムフォレストの 統計的因果推論への応用 サンプルサイズが小さい場合の推論について 今後の展望 22/52
28. 漸近正規性と再帰的分割の分割ルールの関係性について • Athey et al. [2019] の結果は、Specification を満たすような再帰的分割アルゴリ ズムすべてに対して漸近正規性が成り立つことを主張している。 • 分割ルールによる漸近的な推定量の有効性については現状では結果が得られて いない。 有限標本の下での GRF 推定量のバイアスや分散は損失関数に強く依存する。その ため、 • Tree がパラメータ θ(x) をうまく予測するように損失関数を定義し、再帰的分割 を行うような Tree を設計する。 ⇒ 点 x に対して、パラメータ θ(x) の値が同じような値を持つ観測 i の、αbi(x) が 1 になるような近傍を見つける。 αi(x) が、うまく θ(x) がほとんど同じようなサンプルに重みを与えるようにする ためには、木モデルの θ(x) 推定量 ˆ θ(x) と真のパラメータ θ(x) が小さくなるよう な分割を考えればよい。 23/52
29. パラメータ θ(x) に対する Tree の予測誤差 具体的に、P を排反な 2 つのノード C1 , C2 に分割することを考える。親ノード P ⊆ X における推定方程式の解を (ˆ θP , ˆ νP ) とする。 (ˆ θP , ˆ νP ) ∈ argmin θ,ν    Xi∈P ψθ,ν (Oi) 2    (5) 同様にノード C1 , C2 における推定方程式の解をそれぞれ (ˆ θC1 , ˆ νC1 ), (ˆ θC2 , ˆ νC2 ) とす る。このとき、ノード C1 , C2 での推定方程式の解 ˆ θC1 , ˆ θC2 の、新たなテスト点 X に対する 2 乗推定誤差は以下で定義される。 err(C1 , C2 ) = 2 j=1 Pr[X ∈ Cj | X ∈ P] E (ˆ θCj − θ(X))2 | X ∈ Cj (6) X = (X(1), X(2), ..., X(p)) のように p 個の特徴量で構成されていると仮定すると、 最適な分割は以下を満たすような変数 X(k) と閾値 c ∈ R の組 (k, c) として推定で きる。 (k, c) = argmin k,c err(C1 , C2 ) where C1 = {Xi | X(k) i ≤ c, Xi ∈ P} C2 = {Xi | X(k) i > c, Xi ∈ P} (7) 24/52
30. パラメータ θ(x) に対する Tree の予測誤差の問題点 この方法の問題点は 2 つある。 • 1 つ目は、err(C1 , C2 ) には真のパラメータ τ(X) が含まれており、ground-truth が観測できないため、データからこの基準を最適化することができない。 • 2 つ目の問題点はノード P に含まれている訓練データの数を nP とするとき、 最適な分割を探すためには、変数の数 p と、閾値の候補数 nP の積の回数 p × nP 回推定方程式を計算し ˆ θC1 , ˆ θC2 を計算する必要があること 25/52
31. Ground-truth が観測できない問題の解消 Athey et al.[2019] Proposition 1. ノード P の半径が r よりも小さく、分割後のノード内のサンプルサイズが nC1 , nC2 ≫ r−2 を満たすとき、 err(C1 , C2 ) = K(P) − E[∆(C1 , C2 )] + o(r2) where ∆(C1 , C2 ) := nC1 nC2 n2 P ˆ θC1 − ˆ θC2 2 (8) が漸近的に成り立つ。 この結果から、漸近的には θ(X) の ground-truth を含む err(C1 , C2 ) の最小化問題 が、E[∆(C1 , C2 )] の最大化問題と同値であるとみなして、分割ルールを作ればよい ことがわかる。 26/52
32. p × nP 回推定方程式を計算する問題の解消 (1) ノードにおける ˆ θC1 , ˆ θC2 の近似計算を行うことで解消する。ここでは、簡単のた め nuisance parameter ν(x) はないものとして、ψθ,ν (Oi) = ψθ (x) の場合を考える。 いま、ψˆ θ が θ について微分可能なら、ノード Cj における推定方程式を ˆ θP 周りで テイラー展開することによって、 0 = {Xi∈Cj} ψˆ θCj (Oi) = {Xi∈Cj} ψˆ θP (Oi) + (ˆ θCj − ˆ θP ) {Xi∈Cj} ∇θ ψˆ θP (Oi) を得る。ここで、 AP,j = 1 |{i : Xi ∈ Cj}| {Xi∈Cj} ∇θ ψθ (Oi) θ=ˆ θP とおけば、ノード Cj における推定方程式の解は、次のようにして計算することが できることがわかる。 ˆ θCj = ˆ θP − 1 |{i : Xi ∈ Cj}| {Xi∈Cj} A−1 P,j ψˆ θP (Oi) (9) 27/52
33. p × nP 回推定方程式を計算する問題の解消 (2) ここで、スコア関数の微分について、ノード P の半径が十分に小さければ AP,j が ノード P における和で近似できると仮定する。 AP,j ≈ AP = 1 |{i : Xi ∈ P}| {Xi∈P} ∇θ ψθ (Oi) θ=ˆ θP これは例えば因果推論において、ノード P とノード Cj における {Wi − ˆ e(Xi)}2 の 平均がほぼ等しいことを意味しており、平均の意味で処置を受ける確率がほぼ一 定であることを示している。 この仮定のもとで、ˆ θCj は以下のようにして i に依存する項が ψˆ θP (Oi) のみになる。 ˆ θCj = ˆ θP − 1 |{i : Xi ∈ Cj}| {Xi∈Cj} A−1 P ψˆ θP (Oi) (10) よって、ρi = −A−1 P ψˆ θP (Oi) とおく。この ρi を pseudo-outcome をいう [Athey et al., 2019]. 28/52
34. p × nP 回推定方程式を計算する問題の解消 (3) いま、E[∆(C1 , C2 )] に対する推定量として、 ˜ ∆(C1 , C2 ) = E[∆(C1 , C2 )] = 2 j=1 1 |{i : Xi ∈ Cj}|   {i:Xi∈Cj} ρi   2 を構成する。 Athey et al. [2019] Proposition 2. Athey et al. [2019] Proposition 1. の仮定の下で、|AP − ∇ E[ψˆ θP (Oi) | Xi ∈ P]| p → 0 が成り立つならば、 ˜ ∆(C1 , C2 ) = ∆(C1 , C2 ) + oP(max{r2, 1/nC1 , 1/nC2 }) が漸近的に成立する。 この結果は、疑似結果変数 ˜ ∆(C1 , C2 ) を最大化する分割 (k, c) の組を探すことが、 err(C1 , C2 ) 最小化する (k, c) を探すことが漸近的に同等であることを意味してい る。つまり、各ノードで ρi を 1 回だけ計算すればよく、あとは ρi に対して Breiman [1984] の CART アルゴリズムを用いて木を成長させることができる。 29/52
35. （追加）フォレストの重み α(x) に関する理論的・経験的話題として • ランダムフォレストにおける重み αi(x) の形状は、Tree が変数に沿った axis-aligned な分割を生成するため、軸に沿った方向に伸びやすい。 • 例えば X1 , X2 の方向は長くなるが、X1 + X2 方向のような方向で重みのかかる サンプルは少なくなる。 • ランダムフォレストによる重みの推定がうまく働く理由は、Tree の当てはめが 「簡易的な変数選択」になっており、ターゲット θ(x) に関係しない変数が、分割 に用いられにくいためである。 30/52
36. GRF 推定量まとめ ここまでの議論で、次の内容を紹介した。 • GRF 推定量の漸近正規性のための条件 • GRF 推定量に対する分散の計算法と信頼区間の構成 • ターゲットパラメータ θ(x) の異質性が最大になるような再帰的アルゴリズムの 分割ルールの構成（Gradient を用いて疑似結果変数を構成し分割を行う Tree のことを gradient tree という） ここまでの内容で、質問はありませんか？ 31/52
37. 一般化ランダムフォレストの 統計的因果推論への応用
38. 生存時間に対する因果効果の推定 生存時間に対する条件付き因果効果 分布 F からの i.i.d. な変数の組 {Xi, Ti, Ci, Wi} ∈ X × R+ × R+ × {0, 1} を考える。 このとき、Ti = Ti(Wi) を満たす潜在結果変数 {Ti(0), Ti(1)} に対して条件付き因 果効果を以下で定義する。 τ(x) = E[y(Ti(1)) − y(Ti(0)) | Xi = x] (11) ここで、y(·) は生存時間に対する適当な変数変換である。 • y(T) = T とすると、式 (11) は一般的な条件付き因果効果となる。 • ある決められた時刻 h を生存時間の上界として設定し、y(T) = T ∧ h とすると restricted mean survival time に対する条件付き因果効果となる。 • ある時刻 h を設定し y(T) = 1{T ≥ h} とすると、h 時間以上生存する確率に対 する条件付き因果効果となる。 32/52
39. 打ち切りありの観測データ 実際の生存時間解析では、打ち切りを含むような生存時間が観測される。 • Ci を打ち切り時間、打ち切られた生存時間を Ui = Ti ∧ Ci と、非打ち切り指示 関数を ∆i = 1{Ti ≤ Ci} とする。 • データには最大追跡時間が存在するため、結果変数に対する変換 y(·) に対して、 以下の仮定を置く。 仮定 1: Finite horizon 結果変数に対する変換 y(·) は、maximum horizon 0 < h < ∞ のもとで、すべての t ≥ h に対して y(t) = y(h) を満たす。 定義 1 :Effective non-censoring indicator 仮定 1 のもとで、Effective non-censoring indicator を以下で定義する。 ∆h i = 1{(Ti ∧ h) ≤ Ci} 以上より、観測データは {Xi, Ui, ∆h i , Wi} i = 1, 2, ..., n が得られたもとで、いくつ かの仮定の下で条件付き生存因果効果 τ(x) の推定を考える。 33/52
40. 識別可能性（潜在結果変数と処置に対する仮定） 処置効果 τ(x) を識別可能にするための条件として、統計的因果推論における一般 的な仮定を置く。 。 Assumptions (A1-A3) • Potential outcomes : 潜在結果変数 {Ti(0), Ti(1)} は、Ti = Ti(Wi) を almost surely で満たす。 • Ignorability : 処置の割り付けは共変量を条件づけた下でランダム割り付けとみなせる。すな わち {Ti(0), Ti(1)} ⊥ ⊥ Wi | Xi が成り立つ。 • Overlap: 傾向スコア e(x) = Pr[Wi = 1 | Xi = x] は一様に 0 よりも大きく 1 よりも小さ い。すなわち、ある 0 < ηe ≤ 1 2 に対して、 ηe ≤ e(x) ≤ 1 − ηe を満たす。 34/52
41. 識別可能性（生存時間と打ち切り変数に対する仮定） 処置効果 τ(x) を識別可能にするための条件として、生存時間と打ち切りに対して も一般的な仮定を置く。これらの仮定は [Fleming and Harrington, 2011] で与えら れているような生存時間解析の一般的なものである。 Assumptions (A4-A5) • Ignorable censoring : 打ち切りは処置と共変量が与えられたもとで生存時間とは独立である。すな わち、 Ti ⊥ ⊥ Ci | Xi, Wi • Positivity : 処置と共変量が与えらえれたもとで、打ち切りが maximum horizon h よりも前 に起こるかどうかが確率的である。すなわち、ある ηC が存在して、 Pr[Ci < h | Xi, Wi] ≤ 1 − ηC を満たす。 35/52
42. Athey et al., [2019] の causal forest (1) Athey et al., [2019] の causal forest は、前半に述べた「random forest のカーネル 重みづけ的視点」に、次の Robinson [1988] の部分線形モデルの結果を合わせたも のである。 部分線形モデル 任意の x ∈ X に対して処置効果が一定である場合（すなわち、τ(x) = τ） 、仮定 (2)-(4) のもとで、局外パラメータが十分な収束速度 [Chernozhukov et al., 2018] を もって推定できるのであれば、以下の方程式の解として得られる ˆ τ は 1/ √ n の収 束レートを持つ。 n i=1 ψ(c) ˆ τ (Xi, y(Ti), Wi; ˆ e, ˆ m) (12) ψ(c) ˆ τ (Xi, y(Ti), Wi; ˆ e, ˆ m) = [Wi − ˆ e(Xi)][y(Ti) − ˆ m(Xi) − τ × (Wi − ˆ e(Xi))] (13) ここで、e(x) = Pr[Wi = 1 | Xi = x], m(x) = E[y(Ti) | Xi = x] である。 また、ˆ e(Xi), ˆ m(Xi) はこれらのパラメータに対する cross-fitting 推定量である。 36/52
43. Athey et al., [2019] の causal forest (2) Athey et al., [2019] においては、仮定 (A1-A3) のもとで、τ(x) が局所推定方程式 E ψ(c) τ(x) (Xi, y(Ti), Wi; e, m) | Xi = x = 0 を満たすことを用いて、ランダムフォレストの重み αGRF i (x) = 1 B B b=1 1{Xi ∈ Nb(x)} |Nb| を推定し（Nb は点 x を含む tree b の終端ノード） ˆ τ = argmin τ n i=1 αi(x)ψ(c) τ(x) (Xi, y(Ti), Wi; ˆ e, ˆ m) 2 を解くことで、条件付き因果効果 τ(x) に対する因果効果の推定量を構成している。 右側打ち切りを含む生存時間に対しては、打ち切りで重みづけたスコアを構成す る必要がある。 37/52
44. Cui et al., [2023] survival causal forest のスコア関数 Survival causal forest に対するスコア関数は、Steingrimsson et al., [2019] の結果 を用いて、Cui et al., [2023] では以下のように構成している。 ψτ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ) = ∆h i ψc τ (Xi, y(Ui), Wi) SC Wi (Ui ∧ h | Xi) (14) + (1 − ∆h i ) E ψ(c) τ (Xi, y(Ti), Wi) | (Ti ∧ h) > (Ui ∧ h), Wi, Xi SC Wi (Ui ∧ h | Xi) (15) − Ui∧h 0 λC Wi (s | Xi) SC Wi (s | Xi) E ψ(c) τ (Xi, y(Ti), Wi) | (Ti ∧ h) > s, Wi, Xi ds (16) ここで、 • SC w(s | x) = Pr[Ci ≥ s | Wi = w, Xi = x] : 打ち切り過程に対する条件付き生存関 数で、Pr[∆h i = 1 | Xi, Wi, Ti] = SC Wi (Ti ∧ h | Xi) である。 • λC w(s | x) = − d ds log SC w(s | x) : 条件付きハザード関数。 38/52
45. Cui et al., [2023] survival causal forest のスコア関数 これに対して、推定量をプラグインしたスコア関数 ψτ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ; ˆ e, ˆ m, ˆ λC w, ˆ SC w, ˆ Qw) = (Wi − ˆ e(Xi)) × ˆ QWi (Ui ∧ h | Xi) + ∆h i [y(Ui) − ˆ QWi (Ui ∧ h | Xi)] − ˆ m(Xi) − τ(Wi − ˆ e(Xi)) ˆ SC Wi (Ui ∧ h | Xi) − Ui∧h 0 ˆ λC Wi (s | Xi) ˆ SC Wi (s | Xi) ˆ QWi (Ui ∧ h | Xi) − ˆ m(Xi) − τ(Wi − ˆ e(Xi)) ds ここで、 • Qw(s | x) = E[y(Ti) | Wi = w, Xi = x, Ti ∧ h > s] : y(·) で変換された生存時間に 対する条件付き期待値. • ˆ λC w(s | x), ˆ SC w(s | x), ˆ Qw(s | x) は、それぞれ cross-fit によって得られる推定量． • ˆ Qw(s | x) は、推定された生存関数の和によって計算でき、ˆ λw(s | x) は、 − log(ˆ SC w(s | x) の前進差分によって推定できる。 39/52
46. τ(x) = τ の下での推定量の性質 スコア関数 ψτ (·) に対する generalized random forest が漸近正規性を持つために は、部分線形モデル同様に処置効果が一定の場合（τ(x) = τ）に、以下の推定量の 性質を確かめる必要がある。 ˆ τ = argmin τ n i=1 ψτ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ; ˆ e, ˆ m, ˆ λC w, ˆ SC w, ˆ Qw) Cui et al.,[2023] では、 Chernozhukov et al., [2018] の意味でこのスコア関数がネイマ ン直交性を持ち、nuisance component が n1/4 の収束オーダーを持つ時に、ˆ τ が τ に対して √ n の収束レートを持つことと、漸近正規性を有することを示している。 40/52
47. 漸近正規性のための仮定 (1) Assumptions (B1-B5) w = 0, 1 に対して、次の条件が成り立つ。 • 傾向スコア: sup x∈X |ˆ e(x) − e(x)| = op(1) • 条件付き平均: sup x∈X | ˆ m(x) − m(x)| = op(1) • 条件付き生存関数: sup x∈X,s≤h |ˆ SC w(s | x) − SC w(s | x)| = op(1) • 条件付き生存時間: sup x∈X,s≤h |ˆ QC w(s | x) − QC w(s | x)| = op(1) • 条件付きハザード関数と生存関数の比: sup x∈X,s≤h ˆ λC w(s | x) ˆ SC w(s | x) − λC w(s | x) SC w(s | x) = op(1) これらは nuisance parameter の cross-fitting 推定量が一様一致性を持つことの仮 定である。 41/52
48. 漸近正規性のための仮定 (2) Assumptions (C1-C5) w = 0, 1 に対して、次の条件が成り立つ。 • 傾向スコア: sup x∈X E[|ˆ e(x) − e(x)|2] = o(b2 n) • 条件付き平均: sup x∈X E[| ˆ m(x) − m(x)|2] = o(c2 n) • 条件付き生存関数: sup x∈X,s≤h E[|ˆ SC w(s | x) − SC w(s | x)|2] = o(c2 n) • 条件付き生存時間: sup x∈X,s≤h E[|ˆ QC w(s | x) − QC w(s | x)|2] = o(c2 n) • 条件付きハザード関数と生存関数の比: sup x∈X,s≤h E ˆ λC w(s | x) ˆ SC w(s | x) − λC w(s | x) SC w(s | x) 2 = o(d2 n) 推定量をスコア関数を用いて漸近展開した場合に、これらの項または、これらの 項の積によって展開される。よって推定量の漸近的な性質は、これらの項の収束 レートによって決定される。 この直感的な理解は Chernozhukov et al., [2018] のイントロダクションの例がわか りやすいのでここで参考にあげておく。 42/52
49. τ(x) = τ の下での推定量の漸近正規性 Cui et al., [2023] Proposition 1. 条件付き処置効果が x ∈ X に対して τ(x) = τ（定数）であると仮定する。このと き ˜ τ を以下で定義する。 ˜ τ = argmin τ n i=1 ψτ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ; e, m, λC w, SC w, Qw) このとき、仮定 (A1-A5)(B1-B5) および (C1-C5) のもとで以下の収束レートを満たす。 ˜ τ − ˆ τ = op(max{(cn + dn)cn, bn cn, b2 n}) さらに、nuisance パラメータに対する推定量 ˆ e, ˆ m, ˆ λC w/ˆ SC w, ˆ SC w, ˆ Qw が cross-fitting によって推定され、RMSE が n−1/4 の収束オーダーを持つとき、ˆ τ は 1/ √ n の収束 オーダーを持ち、漸近正規性を満たす。 43/52
50. causal survival forest の構成 以上の議論をまとめれば、スコア関数 ψτ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ; ˆ e, ˆ m, ˆ λC w, ˆ SC w, ˆ Qw) を用いて、generalized random forest に対してフォレストの重み αi(x) を計算して、 n i=1 αi(x)ψˆ τ (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ; ˆ e, ˆ m, ˆ λC w, ˆ SC w, ˆ Qw) = 0 を解けば、生存時間に対する異質な処置効果 τ(x) に対する推定量 ˆ τ(x) が構成で きる方針が立つ。 また、gradient-tree に対する疑似結果変数 ρi は、先ほどの議論を用いれば ρi = ψi ˆ τP ×   1 |j : Xj ∈ P| ∑ j:Xj∈P (Wj − ˆ e(Xj))2   1 ˆ SC Wj (Uj ∧ h | Xj) − ∫ Uj∧h 0 ˆ λC Wj (Uj ∧ h | Xj) ˆ SC Wj (Uj ∧ h | Xj) ds     −1 と計算で得きるので、gradient tree で計算可能である。 44/52
51. causal survival forest 推定量の漸近正規性のための仮定 (1) ここで、スコア関数 ψτ (·) が、Athey et al.,[2019] において GRF 推定量が漸近性を 持つように仮定を導入する。 Assumptions (A7) 処置効果の関数 τ(x) が x についてリプシッツ連続であり、nuisance parameter e(x), m(x) も x についてリプシッツ連続であるとする。さらに、 Qw(s | x), λC w(s | x), SC w(s | x) は、w = 0, 1 とすべての s ≤ h において x についてリ プシッツ連続である。 45/52
52. causal survival forest 推定量の漸近正規性のための仮定 (2) さらに先ほどは関数全体に対する平均の議論だったのに対して、各点での収束レ ートについて以下の仮定が必要となる。 Assumptions (A8) w = 0, 1 に対して、次の条件が成り立つ。 • 傾向スコア: E[sup x∈X |ˆ e(x) − e(x)|2] = o(b2 n) • 条件付き平均: E[sup x∈X | ˆ m(x) − m(x)|2] = o(c2 n) • 条件付き生存関数: sup s≤h E[sup x∈X |ˆ SC w(s | x) − SC w(s | x)|2] = o(c2 n) • 条件付き生存時間: sup s≤h E[sup x∈X |ˆ QC w(s | x) − QC w(s | x)|2] = o(c2 n) • 条件付きハザード関数と生存関数の比: sup s≤h E sup x∈X ˆ λC w(s | x) ˆ SC w(s | x) − λC w(s | x) SC w(s | x) 2 = o(d2 n) 46/52
53. causal survival forest 推定量の漸近正規性 Cui el al., 2023 定理 3 仮定 (A1-A8) のもとで、Tree が Athey et al.,[2019] の Specification 1. を満たし、サ ブサンプルサイズ s = nβ, βmin < β < 1 を満たすとする。 このとき、op(max{(cn + dn)cn, bn cn, b2 n}) が、polylog(n/s)−1/2(s/n)1/2 よりも速 いレートで収束するならば、列 {σn} が存在して、任意の x ∈ X に対して、 ˆ τ(x) − τ(x) σn(x) d → N(0, 1) を満たす。ここで、σ2 n = polylog(n/s)−1/2(s/n)1/2 である。 47/52
54. nuisance component の収束レートに関する議論 • b2 n （傾向スコア）: Biau [2012], Wager and Walther[2015] によって、シグナルの次元が 0.54p よりも 小さい場合、b2 n は n−2/(p+2) よりも早い収束レートであることが知られている • c2 n （条件付き平均・条件付き生存関数）: Cui et al., [2022] から、survival forest によって c2 n = n−1/(p+2) が達成できるこ とが知られている • d2 n （ハザード関数と生存関数の比）: Sun et al.[2019] から d2 n = n1/2+κ の収束レートが得られる。ここで κ は特徴空 間の次元 p に依存する。 これらの結果はあくまでのノンパラメトリックな推定手法を利用した場合であり、 パラメトリックモデルによって nuisance component の推定がうまくできるので あれば、推定精度の改善が見込める。 48/52
55. サンプルサイズが小さい場合の推論について
56. 条件付き因果効果の線形射影に対する推論 (1) GRF を現実の問題に適用する際に、議論になるのはサンプルサイズそれほど大き くない（N = 300, 500）の場合である。特に処置関数 τ(x) が複雑な形状をしてい る場合には、τ(x) のバイアスが大きくなり、信頼区間も広くなる。特に、GRF に よる CATE の推定はブレが大きく、問題によっては大きなバイアスを生じることも 少なくないと報告されている。 そこで、[Semenova and Chernozhukov, 2021] では、τ(x) に対する低次元の線形射 影を考えることで、妥当な inference を行う方法を提案している。すなわち、共変 量の集合 Ai ⊂ Xi が与えられたもとで、τ(·) に対する最良線形射影推定量は {β∗ 0 , β∗} = argmin β0,β E [{τ(Xi) − β0 − Aiβ}] で与えらえる。 49/52
57. 条件付き因果効果の線形射影に対する推論 (2) この手法を実装するためには、まず最初に Doubly robust score [Robins et al., 1994] を構成する。 ˆ Γi = ˆ τ(Xi) + ψˆ τ(Xi) (Xi, y(Ui), Ui ∧ h, Wi, ∆h i ) ˆ e(X)(1 − ˆ e(Xi)) ここで、ˆ τ(Xi) は GRF で推定した条件付き処置効果である。 Semenova and Chernozhukov [2021] では、Γi を Ai で回帰した場合に、nuisance components が適当な収束レートを持てば β, β0 に対して妥当な推論が可能である という結論を導いている（ただし、完全なノンパラメトリック推定だと、このレ ートが満たされない） 。 理論的には Semenova and Chernozhukov [2021] の条件を満たさないが、最良線形 推定量を考えることで、処置効果に対する推論が役に立つ可能性は大いにある。 50/52
58. 今後の展望
59. まとめと今後の展望 • 今回の発表では、一般化ランダムフォレストとその理論の応用として causal survival forest について説明した。 • 現在、機械学習手法を用いた因果効果の推定は様々な手法が提案されている。 — Orthogonal Machine Learning: [Foster and Syrgkanis ,2019] — Bayesian causal forest [Hahn, Murray, and Carvalho, 2017]; Tree-based BAFTC [Sun and Song, 2023] — Meta-learners: S-, T-, X- Leaners[Kunzel, 2019]; R-learner[Nie and Wager, 2021]; DR learner[Kennedy, 2020] — Targeted learning [Luedtke and van der Laan, 2016] — causal forests [Wager and Athey, 2018]; generalized random forests [Athey et al., 2019]; causal survival forests [Cui et al., 2023]; Orthogonal random forest [Oprescu et al., 2019] • 解釈性に関する論文も数多く提案されている。 — Sobol-MDA [Benard et al.,2021]; SHAFF [Benard et al.,2022] — Algorithmic variable importance [Williamson et al., 2021]. — Variable importance for causal forest [Benard et al., 2023] — Target learning on VIM [Li, Hubbard and van der Laan, 2023] 機会があれば、BART や VIM の最近の展開もお話します。 51/52
60. 参考文献
61. References (in part) • Athey and Imbens. (2016). ”Recursive partitioning for heterogeneous causal effects”. PNAS, 113, 27, 7353-7360. • Athey et al. (2019). ”Generalized random forests”. The Annals of Statistics, 47, 2, 1148-1178 • Cui et al., (2023). ”Estimating heterogeneous treatment effects with right-censored data via causal survival forests”, JRSS B, 85, 2, 179–211, • Bénard et al. (2022). ”Mean decrease accuracy for random forests: inconsistency, and a practical solution via the Sobol-MDA”. Biometrika, 109,4,881–900. • Chernozhukov et al. (2018). ”Double/debiased machine learning for treatment and structural parameters”. The Econometrics Journal, 21,1,1–68. • Hahn et al. (2020). ”Bayesian Regression Tree Models for Causal Inference: Regularization, Confounding, and Heterogeneous Effects”. Bayesian Anal. 15, 3, 965-1056. • Nie and Wager. (2021). ”Quasi-oracle estimation of heterogeneous treatment effects”. Biometrika, 108, 2, 299–319, 2021. • Semenova and Chernozhukov. (2023). ”Debiased machine learning of conditional average treatment effects and other causal functions”. The Econometrics Journal, 24, 2, 264–289. • Wager and Athey. (2018). ”Estimation and Inference of Heterogeneous Treatment Effects using Random Forests”. JASA, 113, 523, 1228-1242. 52/52
62. ご清聴ありがとうございました。 素敵な 1 日をお過ごしください。 52/52