傾向スコアのモデルに含める共変量選択のアプローチ

傾向スコアのモデルに含める共変量をどうするか？という話題は、2010年ごろから因果推論の話題になるたびに議論されてきました。今回は、それらを整理し、共変量を選択する2つのアプローチについて紹介します。

Transcript

1. 傾向スコアのモデルに含める共変量選択のアプローチ Outcome-adaptive lasso と PACS の紹介 中村　知繁 December 1, 2022 慶應義塾大学
2. Table of contents 1. Introduction 2. Outcome-adaptive Lasso (Shortreed and Ertefaie, 2017) 3. Propensity score adapted covariate selection (Zhou and Jia, 2021) 4. おわりに 1
3. Intro
4. 因果推論における変数選択問題 • 潜在結果変数を用いた因果推論の枠組みは、Rosenbaum and Rubin (1983) において傾向スコアが提案されて以来、多様な分野に浸透した。 • 傾向スコアを用いた解析は、医学・疫学・社会科学・スポーツなど様々 な分野がある。 傾向スコアを用いた解析において、気になる問題の 1 つとして「傾向スコア」 と「結果変数モデル」に含める変数の問題があり、これらは 2 つに大別さ れる。 (1) そもそも、強く無視可能な割り付けを成り立たせるうえで十分な共変量 が観測されているか？ （未観測の交絡の問題） (2) 十分な数の調整に用いる共変量は観測されているが、モデルにどのよう な変数を含めれば推定量の性質が良くなるかわからない（変数選択の問 題） 。 今回は、(2) に関わる話題を紹介する。 2
5. 因果推論における変数選択問題 モデルに対して必要以上の変数を含めることの影響について、よく知られて いるのは、Brookhart et al. (2006), De Luna et al. (2011) および Patrick et at. (2011) の結果である。 De Luna et al. (2011) and Patrick et al. (2011) 結果変数には影響を与えず、処置のみに影響を与える変数を傾向スコアの モデルに含めると、平均処置効果の推定量の分散が大きくなり、バイアス の原因となることがある。　 （⇒）処置変数にのみ影響を与える変数はモデルに含めない方が良い。 Brookhart et al. (2006) 結果変数のみに影響を与える変数を傾向スコアのモデルに含めると、平均 処置効果の推定量の分散が小さくなる。　 （⇒）交絡変数ではないが結果変数のみに影響を与える変数は、傾向スコアの モデルに含めるほうが良い。 3
6. 因果推論における変数選択問題 これらの事実は知られているが、これらの結果を利用した解析手法について は、あまり知られていないのが現状である。 • 本発表では因果推論における変数選択の問題へのアプローチを 2 つ紹介 する。 • Outcome-adaptive lasso by Shortreed and Ertefaie (2017) • Propensity Score Adapted Covariate selection (PACS) by Zhou and Jia (2021) これらのほかにも、 • penalized credible regions (Wilson and Reigh, 2014) • sufficient dimension reduction を用いた次元削減の方法 (Ma et al., 2019) などが提案されているが、ここでの紹介にとどめる。 4
7. Outcome-adaptive Lasso
8. Adaptive Lasso (Zhou, 2006) • outcome-adaptive lasso (OAL) は、傾向スコアに含める変数の選択に、 adaptive lasso (Zhou, 2006) の手法を応用した手法である • まず、Adaptive lasso について導入し、議論を因果推論における傾向ス コアに含める変数選択へと拡張する。 議論においては、以下の記号を用いる。 • Y ∈ R : 結果変数 • X = (X1, ..., Xd) ∈ Rd : 観測された共変量 • A ∈ {0, 1} : 処置変数 • Ya=1, Ya=0 ∈ R : 潜在結果変数 • π(x) = P(A = 1|X = x) : 傾向スコア ただし、E[Xj] = 0 および Var[Xj] = 1 に調整されているとする。 5
9. Adaptive Lasso (Zhou, 2006) ここでは、Adaptive Lasso について説明する。 • サイズ n のデータ (Xi, Yi), i = 1, 2, ..., n が、観測されたとする。 • 観測された共変量のうち、結果変数に対して影響を与える変数の個数が d0 < d であるとする。 • β でパラメトライズされた対数尤度関数を ℓn(β; Y, X) とする。 このとき、Adaptive Lasso 推定量 ˆ β(AL) は以下で定義される。 Adaptive Lasso 推定量 ˆ β(AL) = arg min β    ℓn(β; Y, X) + λn n j=1 ˆ ωj|βj|    (1) ここで、ˆ ωj = |˜ βj|−γ (γ > 0)、λn は正則化パラメータである。 また、˜ β は最尤推定量である。˜ β = arg minβ ℓn(β; Y, X). 6
10. Adaptive Lasso (Zhou, 2006) Adaptive lasso 推定量は、次の仮定 (1),(2) の下で oracle property を満たす。 (1) Yi = XT i β∗ + ε. ただし、ε1, ..., εn は i.i.d. で、平均 0 と、有限な分散 σ2 < ∞ を持つ。また、β∗ = (β∗ 1 , ..., β∗ p ). (2) XTX/n → C, ただし、X は観測データ行列であり、C はある正定値行列。 次に、AL で推定された 0 ではない回帰係数の集合と、0 ではない真のパラメ ータの集合について以下で定義する。 • An = {j : ˆ β(AL) ̸= 0} • A = {j : β∗ j ̸= 0} Adaptive Lasso 推定量の oracle property λn/ √ n → 0 および、λn n(γ−1)/2 → ∞ を満たすと仮定する。このとき、 adaptive lasso 推定量は以下の性質を満たす。 • 変数選択の一致性: limn Pr(A∗ n = A) = 1 • 漸近正規性: √ n(ˆ βA (AL) − β∗ A ) d −→ N(0, σ2C11) 7
11. Adaptive Lasso (Zhou, 2006) ここで C11 は、行列 C の集合 A に対応する行列である。また、βA もパラメー タベクトルの A に対応する要素のみからなるベクトルである。 • adaptive lasso を使いたい場合は、{glmnet} のパッケージを使用すると よいです。 • penalty の引数を OLS 推定量の逆数の γ 乗などとすればよいです。 • この記事が参考になります。 • https://bit.ly/3Uo2QhB 8
12. 因果推論のための記号と仮定 ここからは、adaptive lasso を傾向スコアのモデル選択に応用する outcome adaptive lasso について説明する。まず、以下の仮定を置く。 Assumption : Consistency Y = AYa=1 + (1 − A)Ya=0 Assumption : Unconfoundedness A ⊥ ⊥ Ya|X Assumption : Posititvity 0 < π(x) < 1 今回は、平均処置効果の推定を考える。 Definition: Average Treatment Effect（ATE） θATE = E[Ya=1 − Ya=0] 9
13. 因果推論における変数選択の考え方 観測された共変量は、4 つの背反な添字集合 C, P, I, S ⊂ {1, 2, ..., d} でグル ープに分けることができる。 • X C : 潜在結果変数 Ya と、処置変 数 A の両方に影響を与える共変量 • X P : 潜在結果変数 Ya に影響を与 えるが、処置変数 A には影響を与 えない共変量 • X I : 潜在結果変数 Ya に影響を与 えないが、処置変数 A には影響を 与える共変量 • X S : 潜在結果変数 Ya と、 処置変数 A の両方に影響を与えない共変量 Variables Ya A X C ◦ ◦ X P ◦ × X I × ◦ X S × × このとき、De Luna et al. (2011) の結果から、傾向スコアのモデルは X C , X P を 含み、X I , X S をモデルに含まないようにモデル選択をすることを目指す。 10
14. Outcome-adaptive Lasso • 傾向スコアに対するモデルとして、ロジスティック回帰モデルを仮定す る。 • De Luna et al (2011) の結果から、j ∈ C ∪ P のみを説明変数変数として含 む以下のようなモデルを推定したい。 logit{π(X; ˆ α)} = j∈C∪P ˆ αj Xj + j∈I∪S 0 · Xj (2) • α∗ をモデル (2) の ˆ α の母集団パラメータとする。 • モデル (2) を Reduced PS model と呼ぶ。 Remark 真の傾向スコアのモデルにおける X P の回帰係数は 0 であるが、モデルに含 めることによって、有限標本バイアスを小さくし、推定される傾向スコアと 結果変数の相関が大きくすることで因果効果の推定量の推定分散を小さくし ている。 11
15. Outcome-adaptive Lasso 変数として、j ∈ C ∪ P を選択するために、L1 正則化を用いると次のような問 題が起こる。 • すべての変数をモデルに含めたロジスティック回帰モデルに、L1-正則化 を用いると、X P に対する回帰係数は shrinkage されて 0 に推定されやす くなり、一方で、有効な変数は、j ∈ C ∪ I であるから、これらが推定さ れる。 • そのため、L1 正則化を直接的に用いて推定を行った場合には、因果効果 の推定に有効な変数集合が推定できていない可能性がある。 j ∈ C ∪ P を選択するためには、従来通りの L1 正則化では難しいが、 outcome-adaptive lasso はここに Y の情報を用いた ωj を構築することで、 j ∈ C ∪ P を推定する方法である。 12
16. Outcome-adaptive Lasso いま結果変数に対するモデルを、 Y = ηA + Xβ + ε とし、データへ当てはめて推定された回帰係数を (˜ β, ˜ η) とする。 (˜ β, ˜ η) = arg min β,η ℓn(β, η; Y, X, A). このとき、outcome-adaptive lasso (OAL) は以下で定義される。 Outcome adaptive lasso ˆ α(OAL) = arg min α n i=1 Ai log π(Xi; α) + (1 − Ai) log{1 − π(Xi; α)} +λn d j=1 ˆ ωj|αj|    (3) ここで、ˆ ωj = |˜ βj|−γ（γ > 1）である。 13
17. Outcome-adaptive Lasso (Theoretical aspects) 傾向スコアの推定において、モデルに含めたい変数の添え字集合を A = C ∪ P、取り除きたい変数の添え字集合を Ac = I ∪ S とする。 このとき、一般性を失わずに X = (X A , X Ac ) と並び替えることができ、 A = {1, 2, ..., d0}、Ac = {d0 + 1, d0 + 2, ..., d} とする。ここで、 d0 < d = |C| + |P| + |I| + |S| である。 ここで、すべての変数を predictor として用いたロジスティック回帰モデル のフィッシャー情報行列を I(α∗) = I11 I12 I21 I22 とする。ここで、I11 は d0 × d0 行列であり、Reduced PS モデルの変数に対応 する。 14
18. Outcome-adaptive Lasso (Theoretical aspects) Theorem.1 (Shortreed and Ertefaie, 2017) γ > 1、λn/ √ n → 0 および λn nγ/2−1 → ∞ を仮定する。このとき、 Outcome-adaptive Lasso は適当な正則条件のもとで次の 2 つの性質を満 たす。 1 Consistency in variable selection: lim n P{ˆ αj(OAL) = 0 | j ∈ I ∪ S} = 1 (4) 2 Asymptotic Normality : √ n(ˆ α(OAL) − α∗ A) d −→ N(0, I−1 11 ) (5) Outcome-Adaptive Lasso における罰則は |αj|/|βj|γ という形をしているため、 そのため、j ∈ P である変数に対しては、˜ βj = OP(1) であるから、罰則が小さ くなり、j ∈ I, S である変数に対しては罰則が大きくなる。 また、もともとの adaptive lasso (Zou, 2006) では sparsity の性質を得るため に必要な条件は、λn nγ/2−1/2 → ∞ であったが、Outcome-adaptive lasso で は λn nγ/2−1 → ∞ と条件が強くなっている 15
19. Selecting λn 罰則 λn を選択する際には、グループ間の weighted absolute mean difference （wAMD）を最小にするように選択する (Shortreed and Ertefaie, 2017)。 wAMD(λn) = d j=1 ˜ βj n i=1 ˆ τλn i Ai Xi n i=1 ˆ τλn i Ai − n i=1 ˆ τλn i {1 − Ai}Xi n i=1 ˆ τλn i {1 − Ai} (6) ここで、ˆ τλn i は、λn のもとで Outcome-adaptive lasso によって推定された傾 向スコア ˆ πλn i (·) から作られる inverse weights である。 ˆ τλn i = Ai ˆ πλn i {Xi, ˆ α(OAL)} + 1 − Ai 1 − ˆ πλn i {Xi, ˆ α(OAL)} (7) 改めて、 ˜ βj は、Y に対する X と A の罰則なし条件付き回帰モデルにおける、 X の回帰係数である。すなわち、wAMD は結果変数 Y に対して影響を与える 共変量 X A を強く釣り合わせるような λn を選択するための基準である。 16
20. Propensity score adapted covariate selection
21. Propensity score adapted covariate selection (PACS) PACS は、Outcome-adaptive lasso とアプローチとしては非常に近い手法で ある。Outcome-adaptive lasso では、結果変数 Y に対する線形モデル Y = ηA + Xβ + ε をデータに対して当てはめ、得られた β の推定量 ˜ β を傾向スコアを推定する 際の adaptive lasso の重み ˆ ωj = |˜ βj|−γ として用いることで傾向スコアの変数 を選択するという手法であった。 一方の PACS では、次の 3 段階で変数の選択を行う 1. データから傾向スコアを推定する。( √ n-consistent な推定) 2. 処置群と対照群それぞれの結果変数に対して、傾向スコアによる逆確率 で重みづけた Adaptive-Lasso を当てはめる。 3. adaptive-lasso の回帰係数のいずれかが 0 ではない、変数を propensity score モデルの推定に含めて、PS を推定しなおす。 次のスライドでは、具体的な計算手順について示す。 17
22. PACS - Step1(Propensity score estimation) 傾向スコアの推定においては、推定された傾向スコアが真の傾向スコアに対 して、 √ n 一致性を満たすように推定する。ここでは、ロジスティック回帰 モデルを仮定する。またこのステップにおいては、ロジスティック回帰モデ ルにすべての候補となる変数 X = (X A , X A ) を含める。傾向スコアの推定にお いては、推定された傾向スコアが真の傾向スコアに対して、 √ n 一致性を満 たすように推定する。ここでは、ロジスティック回帰モデルを仮定する。ま たこのステップにおいては、ロジスティック回帰モデルにすべての候補とな る変数 X = (X A , X A ) を含める。 log π(X) 1 − π(X) = j∈A αj Xj + j∈Ac αj Xj (8) ここで、α の最尤推定量を ˆ α とする。 ˆ α = arg max α {Di(XT i α) − log(1 + exp XT i α)} また、ˆ α を代入して得られる傾向スコアの推定量を ˆ π(X) = ˆ π(X, ˆ α) とする。 18
23. PACS - Step1(Propensity score estimation) 傾向スコアの推定においては、推定された傾向スコアが真の傾向スコアに対 して、 √ n 一致性を満たすように推定する。ここでは、ロジスティック回帰 モデルを仮定する。またこのステップにおいては、ロジスティック回帰モデ ルにすべての候補となる変数 X = (X A , X Ac ) を含める。 log π(X) 1 − π(X) = j∈A αj Xj + j∈Ac αj Xj (9) ここで、α の最尤推定量を ˆ α とする。 ˆ α = arg max α {Di(XT i α) − log(1 + exp XT i α)} また、ˆ α を代入して得られる傾向スコアの推定量を ˆ π(X) = ˆ π(X, ˆ α) とする。 19
24. PACS - Step2(β estimation) ステップ 2 では、処置群と対照群に対して、X での Adaptive Lasso 回帰モデ ルを当てはめる。ここでは、Y に対するモデルを正しく特定していることは 仮定しない。まず、処置群において罰則なしの回帰モデルを当てはめる。 (˜ βtreat , ˜ ηtreat ) = arg min β,η i∈{i:Ai=1} 1 ˆ π(Xi) (Yi − η − βTXi)2 次に、X に対する回帰係数 ˜ βtreat を用いて、Adaptive-lasso を当てはめる。 (ˆ βtreat PACS, ˆ ηtreat PACS) = arg min β,η i∈{i:Ai=1} 1 ˆ π(Xi) (Yi − η − βTXi)2 + λn p j=1 ˆ ωtreat j |βj| ここで、λn > 0, ˆ ωtreat j = |˜ βtreat |−γ (γ > 0) である。 同様の計算を対照群においても行う。その場合、傾向スコアの逆数の重みは (1 − ˆ π(Xi)−1 ) とすることに注意する。対照群における adaptive-lasso の回帰 係数の推定量を、ˆ βtreat PACS に対応させて、ˆ βctrl PACS とする。 20
25. PACS - Step3 (variable selection) 最後に、ˆ βtreat PACS および ˆ βctrl PACS を用いて、次のルールに基づいて傾向スコアのモ デルに含める変数を選択する。 変数の選択ルール (Zhou and Jia, 2021) j = 1, 2, ..., p それぞれに対して、ˆ βtreat PACS ̸= 0 かつ ˆ βctrl PACS ̸= 0 の場合に、Xj を傾 向スコアのモデルに加える。 Outcome-adaptive Lasso vs PACS Outcome-adaptive lasso では、結果変数に対する回帰モデルの回帰係数の 推定量で重み付けを行った adaptive-lasso で傾向スコアを推定した。これ は結果変数に対して有効な変数で重み付けをした上で、傾向スコアを Lasso で推定しているため、傾向スコアに対してのみ強く有効な変数の影 響は標本が有限なもとでは残りやすい。 一方で、PACS では結果変数に対して強い影響を持つ変数を Adaptive lasso で推定しているため、傾向スコアに影響を与える変数は adaptive-lasso に よって shrinkage されやすくなる。 21
26. Oracle properties of PACS 次に、PACS による結果の Oracle property について述べる。いま、 ˆ π(x) P −→ π(x) より、結果変数に対する線形モデルが正しいかどうかに関わら ず ˆ βtreat PACS および ˆ βctrl PACS は、それぞれ適当な ˆ βtreat∗ および ˆ βctrl∗ に収束する。 ここで、先程同様に confounder と結果変数にのみ影響を与える変数の添字 集合を A で表し、|A| = p0 < p とする。このとき、次の条件を仮定する。 Assumption: Instrumental variable X A ⊥ ⊥ X A⌋ Assumption: Exclusion restrection Ya=1 ⊥ ⊥ X Ac and Ya=0 ⊥ ⊥ X Ac Linear association condition for Ya=1 すべての j = 1, 2, ..., p0 に対して、Cov(X A , X A )−1Cov(X A , Ya=1)j ̸= 0 Linear association condition for Ya=0 すべての j = 1, 2, ..., p0 に対して、Cov(X A , X A )−1Cov(X A , Ya=0)j ̸= 0 22
27. Oracle properties of PACS Theorem.1 (Zhou and Jia, 2021) γ > 0 に対して、λn/ √ n → 0 および λn nγ/2−1/2 → ∞ を仮定する。このと き、適当な正則条件のもとで次の性質を満たす。 (1) Linear association condition for Ya=1 のもとで、次の結果が成り立つ。 1 limn P{ˆ βtreat PACS,j ̸= 0 | ∀j ∈ A} = 1 2 limn P{ˆ βtreat PACS,j = 0 | ∀j ∈ Ac } = 1 3 √ n(ˆ βtreat PACS,j − βtreat∗) is asymptotically normal. (2) Linear association condition for Ya=0 のもとでも、ˆ βctrl PACS,j に対して、同様 の結果が成り立つ。 (3) さらに、上記 2 つの条件のもとで、 1 limn P{ˆ βtreat PACS,j ˆ βctrl PACS,j ̸= 0 | ∀j ∈ A} = 1 2 limn P{ˆ βtreat PACS,j ˆ βctrl PACS,j = 0 | ∀j ∈ Ac } = 1 が成り立つ。 23
28. PACS についての考察 Theorem.1 (Zhou and Jia, 2021) によって、結果変数に対するモデルを正確に 特定することなく、変数の選択が可能であることがわかる。ただし、傾向ス コアに対しては、 √ n 一致性を満たすような推定量が必要である。これは、 どのパラメトリックな変数選択法でも変わらないので、仕方ないかなという 印象。 こちらも、package 等の実装はないが、R や python で adaptive-lasso を用い ることができれば、問題なく実行可能です。 24
29. おわりに
30. まとめ • 今回の発表では、outcome-adaptive lasso と、propensity score adapted covariate selection を紹介した。 • 傾向スコアを用いた因果推論の問題は、多くの問題が「モデル誤特定」 に焦点が当たっているが、今回の話題は含めるべき変数は何かという視 点からアプローチを行っている。 • 傾向スコアを用いた因果推論を行う際に、少しでもこれらの点を意識で きるようにしていけると良さそうですね。 25
31. ありがとうございました 25