操作変数法のベイズモデリング―『統計的因果推論』(下) 25章より

インベンス・ルービン『統計的因果推論』(下)の25章では操作変数法のベイズモデリングを扱っています。この記事ではその分析をStanで行います。

データと推定したいもの

データはインフルエンザワクチン接種データです。最初の10行は以下のようになっています。
https://github.com/MatsuuraKentaro/reproducing_references/blob/main/Imbens_and_Rubin_2015/chap25/input/flu.csv#L1-L10

• treatment.assigned：操作変数$Z_i$。患者$i$の主治医に対して「インフルエンザのリスクがあると考えられる患者にワクチンを接種するように促す」手紙が送られたか。ここでは無作為に選ばれた医師たちに手紙を送りました。
• treatment.received：処置変数$W_i$。患者$i$がインフルエンザの予防接種をうけたか。
• outcome：結果変数$Y_i$。患者$i$がインフルエンザに関連した病気で入院したか。
• age、copd、heart.disease：共変量ベクトル$X_i$。順に年齢、慢性閉塞性肺疾患の二値指標、過去の心臓疾患に関する指標。

女性だけに絞って分析を行うことにします。操作変数法の文脈では、各ユニットが以下の4タイプのグループ($G_i$)のいずれかに属すると仮定します。

• complier (co)：操作変数を順守する。操作変数が0なら処置も0、1なら処置も1。
• never taker (nt)：常に処置を拒否する。操作変数が何であれ、処置は0。
• always taker (at)：常に処置を受け入れる。操作変数が何であれ、処置は1。
• defier (df)：操作変数に反抗する。操作変数が0なら処置は1、1なら処置は0。

dfが存在しないと仮定し、以下で定義される、「coに対する平均処置効果($\tau_{\text{late}}$)」を推定したいとします（$N_{\text{co}}$はcoの人数です）。

$$\tau_{\text{late}}=\frac{1}{N_{\text{co}}} \sum_{i:G_i=\text{co}} \left(Y_i(1)-Y_i(0)\right)$$

モデル

各ユニットがco/nt/atのいずれになるかについて三値ロジスティック回帰モデルを仮定します。結果変数が0/1のいずれかになるかについてロジスティック回帰モデルを仮定します。難しいのは、例えば$\left(Z_i,W_i\right)=(0,0)$の各ユニットがcoかntか観察できない点です。そのため、データにあてはめる(尤度を計算する)場合には、coだったときのロジスティック回帰モデルとntだったときのロジスティック回帰モデルの混合分布になります。一方、仮定より、$\left(Z_i,W_i\right)=(0,1)$のユニットはatで確定、$\left(Z_i,W_i\right)=(1,0)$のユニットはntで確定です。$\left(Z_i,W_i\right)=(1,1)$の各ユニットについてもcoとatの混合分布とします。

事前分布

人工的なユニットが観察されたと仮定された事前分布を使用します。ここでは、$\left(Z,G,X,Y\right)$ のすべての組み合わせが1回ずつ測定されたとし、その組み合わせ数で割って1ユニット分の重みとします。そのあとで事前分布のユニット数相当のハイパーパラメータを掛けて重みを調節します。データとして観測されないco/nt/atのグループ情報($G_i$)も観測されたと仮定することで三値ロジスティック回帰の係数の推定が安定になります。

Stanコード

Stanコードは以下になりました。
https://github.com/MatsuuraKentaro/reproducing_references/blob/main/Imbens_and_Rubin_2015/chap25/model/model.stan

• 2-7行目：$\left(Z_i,W_i\right)=(0,0)$や$\left(Z_i,W_i\right)=(1,1)$に属するユニットについて、1人分の尤度を計算する関数を定義しています。混合分布なのでlog_sum_exp関数を使います。log_sum_exp関数については詳しく知りたい読者は『StanとRでベイズ統計モデリング』の11章を読んでください。カテゴリカル分布のグループの添字は1がco、2がnt、3がatを表します。よってこの関数は「co」と「ntかatのどちらか一方であるg」が混ざっていてYを生成する分布を表します。
• 9-14行目：上の関数と同じ引数で、あるユニットがcoかどうかの二値を返すベルヌーイ分布を定義しています。coに対する平均処置効果を求めるために、$(0,0)$の中のユニットについてこの分布を使って乱数を生成し、coとなるかどうかを判定します。$(1,1)$の中のユニットについてについても使います。そのため、はじめに$\left(Z_i,W_i\right)=(0,0)$の各ユニットがcoである確率（あるいは$\left(Z_i,W_i\right)=(1,1)$の各ユニットがcoである確率）を求めています。これがちょっと難しいですが、以下の式（の対数を経由して）から求めています。

$$\begin{aligned} &p\left(G_i=\text{co}\middle| Y_i=y, W_i=0, Z_i=0, X_i=x, \theta\right) \\ &= \frac{p\left(G_i=\text{co}, Y_i=y \middle| W_i=0, Z_i=0, X_i=x, \theta\right)}{p\left(Y_i=y \middle| W_i=0, Z_i=0, X_i=x, \theta\right)} \quad (\because 積の公式) \\ &= \frac{p\left(G_i=\text{co}, Y_i=y \middle| W_i=0, Z_i=0, X_i=x, \theta\right)}{\sum_{g} p\left(G_i=g, Y_i=y \middle| W_i=0, Z_i=0, X_i=x, \theta\right)} \quad (\because 周辺化の逆) \\ &= \frac{p\left(\text{co}\middle| x,\gamma\right)\cdot p\left(y\middle| x,\beta_{\text{co,c}}\right)}{p\left(\text{co}\middle| x,\gamma\right)\cdot p\left(y\middle| x,\beta_{\text{co,c}}\right)+p\left(\text{nt}\middle| x,\gamma\right)\cdot p\left(y\middle| x,\beta_{\text{nt}}\right)} \\ \end{aligned}$$

• 19-22行目：$\left(Z_i,W_i\right)=(0,0),(0,1),(1,0),(1,1)$のユニット数をそれぞれ定義しています。今回のモデリングだとここを分割してモデリングしたほうが明快と感じたのでそうしています。
• 23-26行目：$\left(Z_i,W_i\right)=(0,0),(0,1),(1,0),(1,1)$のユニットの添字の配列を定義しています。これがあると53-56行目や76行目や89行目のようにすっきり書くことができます。
• 36-43行目：$Y$に対するロジスティック回帰モデルの係数です。各グループで係数が異なるという柔軟な仮定にしています。
• 47-48行目：三値ロジスティック回帰モデルの係数です。識別性のため、coに相当する係数を0に固定しています。
• 58-71行目：事前分布です。例えば60行目では、「coが観察された」ことに相当します。文献ではここは1ユニット分の尤度ですが、個人的にはここは12人が各グループに4回ずつ観測されたと考えた方が分かりやすいと思うので、4倍したほうがよいと思っています。63～66行目では、「（coが観察されたもとで）coで$(Z,Y)=(0,0),(0,1),(1,0),(1,1)$の4つのデータが観察された」ことに相当します。共変量もすべて使います。71行目では、すべての組み合わせが$12\times N$通りありますので、それで割っています。事前分布のユニット数に相当するハイパーパラメータを30.0に設定しています。
• 93行目：このgenerated quantitiesブロックでは、$\left(Z_i,W_i\right)=(0,0),(1,1)$の各ユニットについて、coである確率を計算してその確率に応じた離散分布よりグループを生成し、coとなった場合には、観測されていない潜在的結果変数を生成します。coのユニット数もカウントしておきます。最後に$\tau_{\text{late}}$を計算しています。

実行するRコード

https://github.com/MatsuuraKentaro/reproducing_references/blob/main/Imbens_and_Rubin_2015/chap25/run-model.R
こちらは特に変わったところはありません。

推定結果

$\tau_{\text{late}}$の分布、全ユニット中coの割合の分布

各係数の事後平均と95%区間

tau_late: -0.15 (-0.36, 0.00)
ITT_W: 0.10 (0.06, 0.13)
ITT_Y: -0.01 (-0.03, 0.00)
a_co_A0,intercept: -3.10 (-6.29, -1.02)
b_co_A0,age: 0.17 (-1.24, 1.75)
b_co_A0,copd: 2.28 (-4.15, 8.94)
b_co_A0,heart: 1.78 (-1.43, 5.59)
a_co_A1,intercept: -3.10 (-6.29, -1.02)
b_co_A1,age: 0.17 (-1.24, 1.75)
b_co_A1,copd: 2.28 (-4.15, 8.94)
b_co_A1,heart: 1.78 (-1.43, 5.59)
a_nt,intercept: -3.26 (-3.75, -2.82)
b_nt,age: -0.17 (-0.35, 0.02)
b_nt,copd: 0.60 (0.08, 1.10)
b_nt,heart: 0.86 (0.33, 1.42)
a_at,intercept: -2.59 (-3.37, -1.91)
b_at,age: 0.02 (-0.32, 0.37)
b_at,copd: -0.08 (-1.02, 0.79)
b_at,heart: 0.52 (-0.27, 1.37)
a_g nt,intercept: 1.93 (1.38, 2.79)
a_g at,intercept: 0.33 (-0.35, 1.33)
b_g nt,age: -0.00 (-0.32, 0.25)
b_g nt,copd: 1.85 (-0.18, 6.66)
b_g nt,heart: -0.03 (-0.98, 0.86)
b_g at,age: 0.21 (-0.16, 0.53)
b_g at,copd: 2.42 (0.24, 7.24)
b_g at,heart: 0.11 (-0.98, 1.15)

tau_lateの事後平均は一致していますが、95%区間は書籍より狭くなっており、書籍の4番目の分析と近い幅になっています。係数の方は微妙にずれてはいますが、だいたい同じような傾向となっていました。

まとめ

操作変数法のモデルベースアプローチとStanによるベイズモデリングはかなり相性が良いと感じました。Enjoy!