🐈

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

2023/11/25に公開

データサイエンス

Stan

因果推論

tech

概要

子供の読解能力を改善することを目的として、教育テレビ番組「The Electric Company」を見せます。各学校において2クラスのペアが選ばれて、片方のクラスには標準的な授業(対照群)を、もう片方のクラスには「The Electric Company」を見せる授業(処置群)をランダムに割り付けます。「The Electric Company」を見せることによる因果効果を推定したいです。表記は一部改変します。

データの説明

生データは見つかりませんでした。そこで、生データからの再現は諦めて、前処理後のデータである表10.1のデータを手動で入力して保存しました。

G: ペアのID
W: 割り付け。0が対照群、1が処置群。
X: 共変量。実験前の点数。
Y: 結果変数。実験後の点数。

モデル

学校(ペア)が8校あるため、8個の層(ブロック)が存在すると考えます。そして潜在的結果変数に学校差がある階層モデルを使って欠測値補完します。このように階層モデルを簡単に組み込めるところがRubin流のベイズモデルベースの因果効果推定のメリットと思います。

ある学校 $j$ における潜在的結果変数が以下に従うと考えます。標準偏差はすべての学校で共通と仮定します：

\begin{pmatrix} Y_{i}(0) \\ Y_{i}(1) \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu(j) + \vec{X}_i\cdot\vec{\beta} \\ \mu(j) + \gamma + \vec{X}_i\cdot\vec{\beta} \end{pmatrix} , \begin{pmatrix} \sigma_{c}^2 & 0 \\ 0 & \sigma_{t}^2 \end{pmatrix} \right)

学校ごとの平均は独立に以下の正規分布に従うと仮定します：

\mu(j) \sim \mathcal{N} \left(\mu_{\text{all}},\sigma_{\mu}^2\right)

$\mu_{\text{all}}$ , $\sigma_{\mu}^2$ , $\sigma_c$ , $\sigma_t$ , $\gamma$ , $\vec{\beta}$ に事前分布を設定します。書籍では分散パラメータに対しては2パラメータ持つ逆カイ二乗分布を使っているのですが、それがわからんのよ～。そこで厳密な書籍の再現は諦めて、標準偏差に対して無情報事前分布を設定します。

コードは以下の通りです。Stanでは事前分布を特に記載しなければ、無情報事前分布(非常に幅の広い一様分布)が設定されることを利用しています。

階層モデルを使うときは、7行目のi2jのようなクラスの添字iを学校の添字jに変換する変数を用意することで、26,28,41,43行目のようにすっきり書くことができるのでした。
22行目で各学校の平均が正規分布に従っています。
残りの記述は以前の記事のモデル1とだいたい同じです。

実行するRコードはこちらです。

4行目：共変量の数が1個であるため、 $N \times 1$ の行列にしています。 $\mu(j)$ があるので切片項は不要な点に注意です。

結果は書籍とほぼ一致しました。

概要

データの説明

モデル

Discussion