🐈

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

2023/11/18に公開

概要

the Student/Teacher Achivement Ratio experiment、略してSTARプロジェクトを扱います。解析単位はクラス(もしくは教師)になります。学校ごとに、標準規模クラス(教師1人につき生徒22～25人)が2クラスと小規模クラス(教師1人につき生徒13-17人)が2～4クラスあります。教師がどっちのクラスに割り付けられるかと生徒がどっちのクラスに割り付けられるかはランダムで決まります。クラスを小規模にすることによる因果効果が推定したいです。表記は一部改変します。

データの説明

データはRのAERパッケージに生徒レベルのデータが含まれています。しかし、ここからの前処理の記述がざっくりすぎて、9.2節のクラスレベルのデータが再現できないです。そこで、生データからの再現は諦めて、前処理後のデータである表9.1のデータを手動で入力し、縦持ち形式(long format)で保存しました。

schoolID: 学校のID
classID: クラスのID
W: 割り付け。0が対照群で標準規模クラス、1が処置群で小規模のクラス。
Y: 基準化された数学の点数のクラス平均。

モデル

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

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

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

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

\begin{pmatrix} \mu_{c}(j) \\ \mu_{t}(j) \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \gamma_{c} \\ \gamma_{t} \end{pmatrix} , \Sigma\right)

$\gamma_c$ , $\gamma_t$ , $\sigma^2$ , $\Sigma$ に事前分布を設定します。 $\gamma_c$ と $\gamma_t$ の事前分布には無情報事前分布を設定します。書籍では $\sigma^2$ には標準逆カイ二乗分布を、 $\Sigma$ には二種類の逆ウィシャート分布を事前分布に設定していますが、Wikipediaの記述やRやStanのパッケージとは異なる関数形で与えられているようで特定できませんでした。そこで厳密な書籍の再現は諦めて、標準偏差に対して無情報事前分布を設定します。分散共分散行列は「無情報事前分布に従う標準偏差」と「一様分布に従う相関係数」から構築します。二変量から変数がもっと増えたらLKJ相関分布を使って構築する方がよいです(『StanとRでベイズ統計モデリング』の10.2.4節参照)。また、弱情報事前分布を使う場合にはStan開発者らによるPrior Choice Recommendationsにしたがって事前分布を設定するのが良いと思います。

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

階層モデルを使うときは、6行目のi2jのような個人の添字iを所属する学校の添字jに変換する変数を用意することで、33,35,49,51行目のようにすっきり書くことができるのでした。
21-23行目で分散共分散行列を構築しています。
27-29行目で各学校の平均が二変量正規分布に従っています。[x0,x1]は行ベクトルを作成します。'は転置を表します。
残りの記述は前回の記事のモデル1とだいたい同じです。

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

結果は、書籍より因果効果の分散が少し小さめに推定されています。

概要

データの説明

モデル

Discussion