🐈

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

2024/02/12に公開

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

データと推定したいもの

データはインフルエンザワクチン接種データです。最初の10行は以下のようになっています。

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コードは以下になりました。

2-4行目： $\left(Z_i,W_i\right)=(0,1)$ や $\left(Z_i,W_i\right)=(1,0)$ に属するユニットについて、全員分の尤度を計算する関数を定義しています。例えば、 $\left(Z_i,W_i\right)=(0,1)$ の場合、ntになる確率がsum(log_p1)で、 $Y$ を観察する確率がbernoulli_lpmf(y | p2)です。
6-11行目： $\left(Z_i,W_i\right)=(0,0)$ や $\left(Z_i,W_i\right)=(1,1)$ に属するユニットについて、全員分の尤度を計算する関数を定義しています。混合分布なのでlog_sum_exp関数を使います。log_sum_exp関数については詳しく知りたい読者は『StanとRでベイズ統計モデリング』の11章を読んでください。
16-19行目： $\left(Z_i,W_i\right)=(0,0),(0,1),(1,0),(1,1)$ のユニット数をそれぞれ定義しています。今回のモデリングだとここを分割してモデリングしたほうが明快と感じたのでそうしています。
33-36行目： $Y$ に対するロジスティック回帰モデルの係数です。各グループで係数が異なるという柔軟な仮定にしています。
37-38行目：三値ロジスティック回帰モデルの係数です。識別性のため、coに相当する係数を0に固定しています。
69-76行目：事前分布です。例えば70行目では、「（atが観察されたもとで）atで $(Z,Y)=(0,0),(0,1),(1,0),(1,1)$ の4つのデータが観察された」ことに相当します。共変量もすべて使います。74行目では、「coとntとatが観察された」ことに相当します。76行目では、すべての組み合わせが $12\times N$ 通りありますので、それで割っています。事前分布のユニット数に相当するハイパーパラメータを30.0に設定しています。
91行目：このgenerated quantitiesブロックでは、 $\left(Z_i,W_i\right)=(0,0),(1,1)$ の各ユニットについて、coである確率を計算してその確率に応じた離散分布よりグループを生成し、coとなった場合には、観測されていない潜在的結果変数を生成します。coのユニット数もカウントしておきます。最後に $\tau_{\text{late}}$ を計算しています。
107-111行目： $\left(Z_i,W_i\right)=(0,0)$ の各ユニットが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}

122-126行目：同様に $\left(Z_i,W_i\right)=(1,1)$ の各ユニットがcoである確率を求めています。

実行するRコード

こちらは特に変わったところはありません。31-44行目では、いったんMAP推定で点推定値を求めて、それをMCMCの初期値に設定しています。

推定結果

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

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

tau_late: -0.16 (-0.37, 0.00)
ITT_W: 0.10 (0.06, 0.13)
ITT_Y: -0.01 (-0.03, 0.00)
beta_co_c,intercept: -3.11 (-6.32, -0.94)
beta_co_c,age: 0.16 (-1.37, 1.77)
beta_co_c,copd: 2.25 (-4.21, 8.95)
beta_co_c,heart: 1.80 (-1.60, 5.66)
beta_co_t,intercept: -4.18 (-7.56, -1.98)
beta_co_t,age: 0.14 (-0.93, 1.48)
beta_co_t,copd: 1.38 (-5.23, 8.35)
beta_co_t,heart: 0.83 (-2.54, 4.53)
beta_nt,intercept: -3.26 (-3.75, -2.82)
beta_nt,age: -0.17 (-0.35, 0.02)
beta_nt,copd: 0.61 (0.05, 1.12)
beta_nt,heart: 0.85 (0.32, 1.40)
beta_at,intercept: -2.61 (-3.36, -1.94)
beta_at,age: 0.02 (-0.31, 0.37)
beta_at,copd: -0.08 (-1.02, 0.81)
beta_at,heart: 0.54 (-0.24, 1.40)
gamma_nt,intercept: 1.94 (1.37, 2.91)
gamma_nt,age: -0.01 (-0.34, 0.26)
gamma_nt,copd: 1.98 (-0.22, 7.49)
gamma_nt,heart: -0.04 (-1.07, 0.89)
gamma_at,intercept: 0.35 (-0.35, 1.48)
gamma_at,age: 0.21 (-0.18, 0.54)
gamma_at,copd: 2.54 (0.19, 8.14)
gamma_at,heart: 0.10 (-1.11, 1.19)

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

まとめ

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