勉強の記録

機械学習、情報処理について勉強した事柄など

Hernán MA, Robins JM (2020). Causal Inference: What If.のPart IIをざっくりとまとめる①

因果推論

Hernan先生のCausal Inferenceを読み進めています．英語で読んでいくとなかなか全体像がつかみづらかったりするので，適宜まとめながら読んでいました．Part IIについては全体の流れは概ね理解できた気がするのでせっかくなので記事にしておきます．Oct 26の版に準拠して該当するページや章を記載していますが，今後の改訂によって多少ずれるかもしれません．

Part Iは無限のpopulationを仮定した理想的な状況での因果推論の定性的な議論です．Part IIは実問題に応用できるようどのようにモデル化して解析するかの手法が紹介されており，Part IIIは時間依存性介入や時間依存性交絡についてのSectionという3部構成になっています．

Part Iのcausal graph(DAG)やconfounding，counterfactual outcomeなどの説明は割愛します．

Notation

$A$: treatment

$L$: confounders

$Y$: outcome

$Y^{a}$: conterfactual outcome

TP: Technical Point

FP: Fine Point

この本の因果推論における前提 (Ch13.5)

  • identifiability conditions: exchangeability, positivity, consistency (Chapter 3)
  • all variables are correctly measured. (Chapter 9 measurement error)
  • all models are correctly specified (Chapter 11). model misspecification.

これらの条件を確認するのが，ingestigator's most important taskであるとのこと．この本で説明されている手法の多くは上記の条件を前提としています．

また，本文中にはrandomでないcensoringに対応する手法も記載がありますがconfoundingへの対処と本質的にあまり差がない気がするのでこの記事では省略します．

conditional exchangeabilityを仮定する推定方法

因果推論における最大の関心領域であり，この仮定を前提とした手法が多いです．これは$Y^{a} ⊥ A|L$という条件であり，$L$が一致している人々の中では実際に治療された群も治療されなかった群も潜在的なアウトカムが等しいということです．この条件のもとでは$E[Y|A,L] = E[Y^{a}|L]$が成立するので，$L$が等しい人々のなかでは単純に治療群と非治療群のアウトカムを比較すれば$association=causation$となって交絡なく治療効果を推定することができます．実際には，$E[Y|A,L] = E[Y^{a}|L]$はmean exchangeabilityと呼ばれ，独立ではなくてもこちらさえ成立していれば以下に扱うようなaverage causal effectの推論には十分なようです（Technical Point 2.1参照）．

なお，Part IのCh.8あたりに詳しく説明されていますが，$A→Y$の因果効果推定に際して調整すべき交絡因子(confounding factors)は一意に定まるわけではありません．例えば，$A←L1→L2→Y$というcausal graphを考慮したとき$L={L1}$も，$L={L2}$も，$L={L1, L2}$もいずれであっても，調整することで十分に交絡に対処できます(Ch.8)．

この本では事前知識などにより構築したcausal graphのもとで調整すべき$L$を検討するという立場を取っており，データからcausal graphを構築する因果探索的なアプローチはこの本のscope外です．以下に，conditional exchangeabilityを仮定する推定方法を挙げていきます．本文ではsubgroup上のATEや，ある共変量によるeffect modificationを推定する方法も登場しますが，この記事ではまずaverage causal effect (ATE)を推定する方法を取り上げます．

Nonparametric estimators of the conditional mean(Ch.11.3)

まずひとつ目はこれ．Part Iでも取り上げられました．すべての$A$と$L$が離散的で、かつサンプル数が無数にある（というやや非現実的な）とき、すべての$A, L$の組み合わせに対して十分なサンプルがあるのでその各セルに対してconditional mean: $E[Y|A, L]$を直接計算することで$E[Y^{a}|L]$が求められます．そして後で取り上げるstratification(=g-formula)の考え方を使って，全サンプルの$L=l$に対して$E[Y|A=a, L=l]$を計算して平均を取ることでaverage causal effect (ATE)を推定することができます．

モデル化の必要性

しかしA, Lに連続量を含む場合やLが高次元な場合など，$E[Y|A, L]$を$A, L$の組み合わせごとに求めていてはうまく計算できないことが多く，何らかのモデルを仮定してsmoothingする必要があります（言い換えると，Lが少し違うサンプルの情報を借りてきて，$E[Y^{a}|A,L]$を推定しなければいけない）．

モデル化の方法としては，この本での主題となるg-methodsとstratification-based methodsに大きく2つに分けることができます(p95)．以下，一つずつ見ていきます．

g-methods

g-methods①: IP weighting

IPWによるpseudo-population上で$E[Y|A] = \theta_0 + \theta_1A$を回帰することで，marginal structural mean model: $E[Ya] = \beta_0 + \beta_1a$を推定します．marginal structural mean modelは個々のサンプルについてのYの回帰式ではなく，Yの平均構造についての回帰式だと思います．p157に，もし知りたいパラメータが母集団全体のATEであれば，marginal strucural mean modelは共変量の項を含まないとあります．（なおmarginal structural mean modelはAが2値であればAについてはsatured modelであり，本質的にはmodelではなくsmoothingの効果はない．Aがcontinuousな場合は何らかのsmoothingがなされる）．

Aについてsaturedなmodelでは関係ありませんが，Aについてもモデル化するような場合は，stabilized IP weighingと言って，$\frac{1}{\pi[a|l]}$のかわりに$\frac{f(A=a)}{\pi[a|l]}$を用いた方がvalidな推定ができると書かれていました(Ch.12.3)．

IP weighing法におけるモデル化の核は，$E[A|L] =\pi (L)$のモデル化であり，logistic regressionなどが用いられることが多いです．Propensity scoreを用いる手法全てに共通することですが，$L$にoutcomeには弱く影響しtreatmentに強く影響するような（言ってしまうと操作変数に近いような）変数を含んでいると推定が不安定になることが欠点です．

また，$A$がcontinuousな場合は$\pi(L)$は何らかの確率密度分布probability density function (pdf)を出力する必要があり，かつこの誤設定に敏感なのでかなり慎重にならなければならないと述べられています．

(p.156) Unfortunately, pdfs are generally hard to estimate correctly, which is why using IP weighting for continuous treatments will often be dangerous. One should be careful when using IP weighting for continuous treatments because the effect estimates may be exquisitely sensitive to the choice of the model for the conditional density f(A|L).

IPW自体はたくさんのblogでも紹介されているので深入りしませんが，一応TP3.1の式変形を追っておきます．

$$ \begin{eqnarray} E[Y^{a}] &=& \sum_{l}E[Y^{a}|L=l]Pr[L=l]\\ &=&\sum_{l} {E[Y|A=a, L=l] Pr[L=l]}\\ &=&\sum_{l} \frac{1}{\pi [a|l]}{E[Y|A=a, L=l] \pi[a|l]Pr[L=l]}\\ &=&\sum_{l} \frac{1}{\pi [a|l]} I(A=a)Y Pr[L=l]\\ &=& E\left[\frac{I(A=a)Y}{\pi[a|l]}\right]\\ &=& \frac{\sum \frac{I(A=a)Y}{\pi[a|l]}} {\sum \frac{1}{\pi[a|l]}} \end{eqnarray} $$

$E[Y|A=a, L=l] \pi[a|l] = I(A=a)Y$の部分の変形がややトリッキーですが，$A=a$の時の期待値に$A=a$となる割合をかけているので，この部分は，

$$ I(A=a)Y = \left\{ \begin{array}{} Y & (A=a)\\ 0 & (A\neq a) \end{array} \right. $$ の期待値に等しくなります．また，最後は${\sum \frac{1}{\pi[a|l]}}$がtotal populationのNに等しくなるためです．

g-methods②: standardization and the parametric g-formula (Ch.13)

こちらはIPWより見かける機会は少ないかもしれません． $$ ATE = \sum_{l} E[Y|A=a, L=l] Pr[L=l] $$ IP weighingではサンプル毎に重みをかけてconterfactualを推定しましたが，L=lの層ごとに考えて層の重みでsummationしています．summationはLでの期待値とみなせるので，L=lの層内でYの期待値(conditional mean outcome)を取った後にLについて期待値を取ったdouble expectationと解釈できます．

実はIP weightingの式変形の途中で同じ式が出現しており，結局は同じものを求めていますが，モデル化する部分が違うのです．

上式のconditional mean outcomeにそのestimateを代入（plug-in）したものを，plug-in g-formulaと呼びます．例えば，$\hat{E}[Y|A=a, L=l]$として，$A, L$を変数としたoutcome regression model（例：$A$，$L$の一次項，二次項，一次の交互作用項などを入れた線形モデル）を考えることができます．

なお，$Pr[L=l]$は実際には推定する必要はなく，サンプルの全ての$l$について$\hat{E}[Y|A=a, L=l]$を計算して平均すればATEを求めることができます（empirical distribution）．

また，後のstratification and outcome regressionでも，outcome regression model: $E[Y|A, L]$が出てきます．standardization (g-formula)のポイントは，$Pr[L=l]$をかけて和を取ることで母集団全体についての平均（＝期待値）を取っているという部分だと思います．

g-methods③: doubly robust estimators（TP13.2）

doubly robust estimatorにもいろいろあるようで，FP13.2とTP13.2で少し違うモデルが取り上げられていますが，ここではTP13.2のものを載せておきます．

$$ \pi(L) = Pr[A=1|L=l]\\ b(L) = E[Y|A=1,L] = E[\frac{AY}{\pi(L)}] $$

として，それぞれの推定式を$\hat{\pi}(L)$，$\hat{b}(L)$と表します．ここで，doubly robust estimatorは， f:id:tmitani-tky:20191104194625p:plain となります（Texclipではうまく表示できるがhatenaだとうまく表示できず...）．

証明は追えていないのですが，この時のbiasはIPWのerror$\frac{1}{\pi(l)}-\frac{1}{\hat{\pi}(l)}$と，outcome modelのerror$b(l)-\hat{b}(l)$との積になるため，doubly robustと呼ばれます．

g-methods④: structural nested models and g-estimation (Ch.14)

これも星野本などではあまり取り上げられていない手法だと思います．ただ，時間依存性介入のような複雑な状況でなければover killであるとも述べられており，Part IIIへの準備という側面が強いのかもしれません．p180には"In fact, structural nested models of any type have rarely been used"とあります．

これはconditional exchangeabilityを逆手に取ったような手法です．conditional exchangeabilityのもとでは$A$と$Y^{a}$は独立であり，Lに加えて$Y^{a}$を条件づけたところでA=aの分布は変わらない，つまり， $$ Pr[A=1|Y^{a=0},L]=Pr[A=1|L] $$ が成り立っているはずで， $$ logit Pr[A=1|Y^{a=0},L] = \alpha_{0} + \alpha_{1}Y^{a=0}+\alpha_{2}L $$ を推定すると$\alpha_{1}=0$になるはずだ，というロジックです．そこで$Y^{a=0}$の候補を虱潰しにサーチして，$ \alpha_{1}$が0になればそれが正しい$Y^{a=0}$だろうというわけです．

かといって任意の$Y^{a=0}$を調べることは計算量が爆発しますし，解も無限に出てきてしまいそうです．そこで，causal effectの式を少し変形してstructural model: $$ E[Y^{a}-Y^{a=0}|A=a,L]=\beta_{1}a $$ を考えます．このmodelはstructural nested mean modelとも呼ばるそうですが，"nested"はtime-varying treatmentでのみnestedらしいのであまり気にしないでよさそうです(p173)．

仮想的な$Y^{a=0}$として， $$ H(\psi^{\dagger}) = Y -\psi^{\dagger}A $$ を考えます．$\beta_{1}$としてありえそうな値の範囲として，（例えば教科書の例では）-20から+20まで0.01刻みで計算しています．その値を$\psi^{\dagger}$に代入した$H(\psi^{\dagger})$を使ってlogistic regressionまで解き，$ \alpha_{1}=0$にもっとも近い結果が得られる$\psi^{\dagger}$が，$\beta_{1}$の不偏推定になるということのようです．

なお，ここでは単一パラメータのモデルのみ取り上げましたが，複数のパラメータを持つモデルでは制約が一つだと複数の候補がでてきてしまいます．そこでパラメータ数と同じだけの制約を設けるために，$E[Y^{a}-Y^{a=0}|A=a,L]=\beta_{1}a + \beta_{2}aL$であれば，$logit Pr[A=1|Y^{a=0},L] = \alpha_{0} + \alpha_{1}Y^{a=0}+\alpha_{1}Y^{a=0}L + \alpha_{3}L$など項が対応したlogstic modelを考えてあげると良いです．

また，パラメータが増えると計算量がexponentialに増えてしまいますが，実際にはしらみつぶしに計算せずにclosed formで解くことができるそうなので，気になる方はTechnical Point(14.2)を見てみてください（よく分かりませんでした...）．

長くなってきたので分けます．stratification-based methodsは次の記事で．