🐈

『インベンス・ルービン統計的因果推論（上）』の8.10節をStanで再現する

2023/11/11に公開

データサイエンス

Stan

因果推論

tech

概要

『インベンス・ルービン統計的因果推論（上）』はランダム割り付けの場合の因果効果推定を扱っています。ランダム割り付けの場合だけで約300ページ！この本の特徴として「DAGがいっさい登場しない」「潜在的結果変数の欠測値補完の問題を通して因果効果を推定する」「割り付けメカニズムを確率変数と見なす」などがあります。この記事では、8章「完全無作為化実験に対するモデルベースの推論」の問題をStanで実装します。表記は一部改変しています。

データの説明

具体的なデータの説明をします。ここでは職業訓練プログラム (Natinal Supported Workプログラム, NSWプログラム) によって、どれほど収入が上がったかの効果を推定します。プログラムへの割り付けはランダムです(完全無作為化)。データはRのMatchingパッケージの中に含まれていますし、Web上 (直接リンク注意)から入手することも可能です。

結果変数：訓練プログラム終了後の1978年の収入(re78)
使用する共変量：年齢(age)、教育年数(education)、結婚しているかどうか(married)、高校中退かどうか(nodegree)、人種(black)、1974年の収入(re74)、1974年の収入がゼロより大きいかどうか(re74pos)、1975年の収入(re75)、1975年の収入がゼロより大きいかどうか(re75pos)

ここでは作図を省略しますが、結果変数のヒストグラムはゼロ過剰で右に尾を引いた分布となっています。

求める因果効果

この本の方針としまして、 $i$ さんが潜在的結果変数 $Y_{i}(0)$ 、 $Y_{i}(1)$ を持っていると考えます。 $i$ さんが対照群( $0$ )に割り付けられたときには $Y_{i}(0)$ が観測され、処置群( $1$ )に割り付けられたときには $Y_{i}(1)$ が観測されると考えます。個人レベルでは必ずどちらか欠測になるため、個人レベルの因果効果を求めることは(非常に強い仮定を置かない限りは)できません。そこで集団レベルの因果効果を考えることにします。集団レベルの因果効果とは、「その集団全員が処置を受けなかった場合の結果変数の平均」と「その集団全員が処置を受けた場合の結果変数の平均」を比較することで求められる効果になります。こちらは比較的緩い仮定を置くことで求めることができます。これらの仮定の詳細については書籍に譲り、ここでは説明しません。

モデル

8.10節ではモデルを4種類試しています。

モデル1: 共変量なし、潜在的結果変数が相関あり、分散が共通
モデル2: 共変量なし、潜在的結果変数が独立、分散が別々
モデル3: 共変量あり、潜在的結果変数が独立、分散が別々
モデル4: 共変量あり、潜在的結果変数が独立、分散が別々、ゼロ過剰対数正規分布を使う

以下で一つずつ詳しく見ていきます。

モデル1

潜在的結果変数が以下の二変量正規分布に従うと仮定します。

\begin{pmatrix} Y_{i}(0) \\ Y_{i}(1) \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu_{c} \\ \mu_{t} \end{pmatrix} , \begin{pmatrix} \sigma^2 & \sigma^2 \\ \sigma^2 & \sigma^2 \end{pmatrix} \right)

ここで、未知パラメータは $\mu_{c}$ , $\mu_{t}$ , $\sigma^2$ です。添え字の $c$ はcontrolの、 $t$ はtreatmentの頭文字と思います。 $\mu_{c}$ と $\mu_{t}$ の事前分布には平均がゼロで分散が $100^2$ の無情報事前分布を設定します。 $\sigma^2$ の事前分布には、パラメータ $1$ と $0.01$ の逆ガンマ分布を設定するとのことです。

分散共分散行列の相関係数が1に決め打ちされているのは理由があります。6章などで、相関が1のときに、推定された因果効果の分散が最も大きくなる(保守的な評価になる)ことが示されているからです。

このモデル式にしたがって欠測となった潜在的結果変数を補完して、集団レベルの因果効果を推定します。コードは以下の通りです。

5行目で分散共分散行列の相関係数をデータとして与えます。44-45行目と47-48行目で観測されなかった結果変数の欠測を補完しています。片方が観測された二変量正規分布なので、条件付き分布から生成している点に注意です。分からない人は「二変量正規分布　条件付き分布」とかで検索してみてください。36行目が超母集団に対する平均因果効果に相当します。37行目が持っているデータに対応する集団(有限標本; finite sample)に対する平均因果効果です。38-40行目は、因果効果の分位点ではなくて、対照群の分位点と処置群の分位点の差を見ます(詳しくは書籍を参照してください)。

実行するRコードは以下の通りです。

5行目で相関係数を与えています。Rでは標準偏差0の正規分布でも正常に動きますが、Stanではエラーになりますので、小さい値を引いておきます。

結果は書籍と完全一致です。

モデル2

潜在的結果変数が以下の二変量正規分布に従うと仮定します。

\begin{pmatrix} Y_{i}(0) \\ Y_{i}(1) \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu_{c} \\ \mu_{t} \end{pmatrix} , \begin{pmatrix} \sigma_{c}^2 & 0 \\ 0 & \sigma_{t}^2 \end{pmatrix} \right)

$\sigma_{c}^2$ と $\sigma_{t}^2$ の事前分布には、モデル1と同様にパラメータ $1$ と $0.01$ の逆ガンマ分布を設定するとのことです。

コードはモデル1とほぼ同じなので省略します。結果は0.01だけずれているぐらいで書籍とほぼ一致です。

モデル3

潜在的結果変数が以下の二変量正規分布に従うと仮定します。

\begin{pmatrix} Y_{i}(0) \\ Y_{i}(1) \end{pmatrix} \bigg| \vec{X}_i \sim \mathcal{N} \left( \begin{pmatrix} \vec{X_i}\cdot\vec{\beta_{c}} \\ \vec{X_i}\cdot\vec{\beta_{t}} \end{pmatrix} , \begin{pmatrix} \sigma_{c}^2 & 0 \\ 0 & \sigma_{t}^2 \end{pmatrix} \right)

$\vec{\beta}_{c}$ と $\vec{\beta}_{t}$ の事前分布には平均がゼロで分散が $100^2$ の無情報事前分布を設定します。 $\sigma_{c}^2$ と $\sigma_{t}^2$ の事前分布はモデル2と同じです。

コードは以下の通りです。

6行目に共変量が追加されました。22-23行目で共変量の係数をベクトルで定義しています。29行目、31行目において、正規分布の平均が共変量と係数の線形結合になっています。

結果は0.01～0.02だけずれているぐらいで書籍とほぼ一致です。

モデル4

上記のモデル1～3では、観測した結果変数においてゼロが多いのを特にケアしないで正規分布をあてはめていましたが、ここではゼロ過剰対数正規分布をあてはめます(書籍では「ゼロ過剰対数正規分布」とは書いていませんが、そのようにとらえた方が見通しが良いと思っています)。すなわち、ゼロより大きくなるかがベルヌーイ分布で決まり、いったんゼロより大きいとなれば値は対数正規分布に従うという分布です：

\operatorname{Pr}\left(Y_i(0)>0 \mid \vec{X_i}\right)=\frac{1}{1+\exp \left(-\vec{X_i}\cdot\vec{\gamma_c}\right)}

Y_i(0) \mid Y_i(0)>0, \vec{X_i} \sim \operatorname{logNormal}\left(\vec{X_i}\cdot\vec{\beta_c}, \sigma_c^2\right)

処置群も同様にゼロ過剰対数正規分布をあてはめます。

コードは以下の通りです。

2-8行目でZero-Inflated Log Normal(ZILN)分布の尤度を定義しています。10-17行目でZILN分布の乱数生成を定義しています。

実行するRコードは以下の通りです。

結果は…微妙に一致しません！この一致しない問題は過去にStanのForumでも取り上げられています。ただ、共変量の回帰係数はほぼ一致します。公式の実装のコードが公開されていないので詳細は分かりません。みなさんは本や論文を書いたら実装のコードは必ず公開しましょうね。

概要

データの説明

求める因果効果

モデル

モデル1

モデル2

モデル3

モデル4

Discussion