Chap. 10 Pairwise Randomized Experiments
10.7 Model-based analysis of pairwise randomized
experiments
Imbens & Rubin (2015)の10.7節「Model-based analysis of pairwise randomized experiments」におけるベイズ分析の説明に従って、Pythonスクリプトを書いてみた。
対データ（統制群、実験群）は次の確率に従うとされる（Imbens & Rubin, 2015, p. 232）。
ここで、unit iはj番目の対のものである。統制群と実験群の差は、
で表されている。
は共変量である。
各
は次式の確率モデルに従うとされる（Imbens & Rubin, 2015, p. 232）。
以上の確率モデルを表すStanスクリプトを、リスト１に示した。
リスト１ ベイズ分析用Stanスクリプト（ファイル名Pairwise.stan）
data {
int N;
array[N] vector[2] Y01;
array[N] vector[2] X01;
}
parameters {
real vmu; real vg; real vb;
real<lower=0.0>
vsgm_mu2; real<lower=0.0> vsgm_c2;
real<lower=0.0>
vsgm_t2; array[N] real muj;
}
transformed parameters {
real mu; real g; real b;
real sgm_mu2; real sgm_c2;
real sgm_t2;
mu = vmu * 10; g = vg * 10;
b = vb * 10;
sgm_mu2 = vsgm_mu2 + 0.001;
sgm_c2 = vsgm_c2 + 0.001;
sgm_t2 = vsgm_t2 + 0.001;
}
model {
vmu ~ normal(0.0, 10.0);
vg ~ normal(0.0, 10.0); vb ~
normal(0.0, 10.0);
vsgm_mu2 ~
scaled_inv_chi_square(1.0, 1.0);
vsgm_c2 ~
scaled_inv_chi_square(1.0, 1.0);
vsgm_t2 ~
scaled_inv_chi_square(1.0, 1.0);
for (j in 1:N) muj[j] ~
normal(mu, sqrt(sgm_mu2));
for (i in 1:N) {
Y01[i][1] ~ normal(muj[i] + X01[i][1] * b, sqrt(sgm_c2));
Y01[i][2] ~ normal(muj[i] + g + X01[i][2] * b, sqrt(sgm_t2));
}
}
の事前分布は、平均０，標準偏差100の正規分布であるが（Imbens &
Rubin, 2015, p. 233）、PyStan 3の警告を避けるために、スケーリングしている。
、
、
の事前分布は、パラメータが１と１の逆カイ二乗分布を設定している（Imbens & Rubin, 2015, p. 233）。
Imbens & Rubin (2015, p. 232)では、多変量正規分布が用いられているが、PyStan 3（PyStan 3.9.0）ではエラーが頻発したので、1変量正規分布を用いた。共分散が０であるので、１変量正規分布に分解できる。
リスト１を用いてImbens & Rubin (2015)のTable 10.1のデータを分析するPythonスクリプトをリスト２のように用意した。
リスト２ リスト１のStanスクリプト（Pairwise.stan）を用いるPythonスクリプト（PyStan 3対応）
import numpy as np
import stan
import arviz as az
import pandas as pd
import matplotlib.pyplot as plt
#
#
From Table 10.1, not from Table 10.2
#
Y01 = [[54.6, 60.6], [56.5, 55.5],
[75.2, 84.8], [75.6, 101.9],
[55.3,
70.6], [59.3, 78.4], [87.0, 84.2], [73.7, 108.6]]
X01 = [[12.9, 12.0], [15.1, 12.3],
[16.8, 17.2], [15.8, 18.9],
[13.9,
15.3], [14.5, 16.6], [17.0, 16.0], [15.8, 20.1]]
with open('Pairwise.stan') as f_stan:
sm =
stan.build(f_stan.read(),
data = {'N':len(Y01), 'Y01':Y01, 'X01':X01})
fit = sm.sample(num_samples = 20_000)
i_d = az.from_pystan(posterior = fit,
posterior_model = sm )
print('\nSummary =\n',
az.summary(i_d))
リスト２ではImbens & Rubin
(2015)のTable 10.1のデータを分析の入力値としている。Y01がPost-Test Score 
を表し、X01がPre-Test Score 
を表している。
リスト２のスクリプトを実行すると、以下の出力を得る。
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
vmu
-2.229 0.856 -3.927 -0.673 0.019 0.015 2224.0 1302.0 1.00
vg
0.859 0.161 0.559 1.159 0.002 0.002 5029.0 8306.0 1.00
vb
0.587 0.054 0.493 0.695 0.001 0.001 2748.0 988.0 1.00
vsgm_mu2 25.573 21.874 0.197 58.667 0.257 0.182 4482.0 5541.0 1.00
vsgm_c2 13.845 13.780 0.160 33.840 0.191 0.135 1964.0 899.0 1.00
vsgm_t2 5.798 8.825 0.138 19.203 0.125 0.089 567.0 299.0 1.01
muj[0] -19.756 7.205 -33.862 -7.085 0.138 0.099 2729.0 3394.0 1.00
muj[1] -26.751 7.269 -40.686 -13.646 0.142 0.102 2577.0 3442.0 1.00
muj[2] -24.114 8.952 -41.776 -7.667 0.191 0.161 2294.0 805.0 1.00
muj[3] -18.122 9.657 -37.593 -1.096 0.203 0.169 2391.0 818.0 1.00
muj[4] -26.788 7.844 -42.413 -12.447 0.175 0.145 2122.0 787.0 1.00
muj[5] -26.545 8.319 -42.977 -11.017 0.186 0.157 2133.0 774.0 1.00
muj[6] -17.431 9.296 -35.978 -0.633 0.186 0.156 2506.0 825.0 1.00
muj[7] -18.801 9.925 -38.902 -1.297 0.217 0.181 2215.0 809.0 1.00
mu
-22.294 8.560
-39.268 -6.730 0.186 0.153 2224.0 1302.0 1.00
g
8.590 1.608 5.593 11.592 0.022 0.016 5029.0 8306.0 1.00
b
5.874 0.541 4.929 6.949 0.010 0.008 2748.0 988.0 1.00
sgm_mu2 25.574 21.874 0.198 58.668 0.257 0.182 4482.0 5541.0 1.00
sgm_c2 13.846 13.780 0.161 33.841 0.191 0.135 1964.0 899.0 1.00
sgm_t2 5.799 8.825 0.139 19.204 0.125 0.089 567.0 299.0 1.01
（g）と
(b)の平均値は、8.6と5.9であり、Table 10.4に記載されている値と一致する。ただし、MCMCサンプリングの場合、プログラムの実行ごとに値はランダムに変わりうることに注意されたい。
muの値は、Table 10.4における-9.2とはっきりと異なる。
の推定を、次のスクリプトで行ってみる。
import numpy as np
Y01 = [[54.6, 60.6], [56.5, 55.5],
[75.2, 84.8], [75.6, 101.9],
[55.3,
70.6], [59.3, 78.4], [87.0, 84.2], [73.7, 108.6]]
X01 = [[12.9, 12.0], [15.1, 12.3],
[16.8, 17.2], [15.8, 18.9],
[13.9,
15.3], [14.5, 16.6], [17.0, 16.0], [15.8, 20.1]]
b = 5.9
g = 8.6
Res = np.empty((8,2))
for i in range(8):
Res[i][0] = Y01[i][0] -
(X01[i][0] * b)
Res[i][1] = Y01[i][1] - (g +
X01[i][1] * b)
print('Res =')
for v in Res:
print(v, f' mean = {np.mean(v):.3f}')
実行すると以下の出力を得る。
Res =
[-21.51 -18.8 ] mean = -20.155
[-32.59 -25.67] mean = -29.130
[-23.92 -25.28] mean = -24.600
[-17.62 -18.21] mean = -17.915
[-26.71 -28.27] mean = -27.490
[-26.25 -28.14] mean = -27.195
[-13.3 -18.8] mean = -16.050
[-19.52 -18.59] mean = -19.055
リスト２の実行により得られた
mu -22.294
は、上の出力と合うが、Table 10.4に記載されている
は、不適切であると思われる。校正ミスであろう。
Table 10.1に記載されている値は、Table 10.2のものとは異なるものがあり、Table 10.2の値に対して分析を行うと、Table 10.4に記載されているとの値は得られない。校正ミスが散見される。本サイトの構築者は、
First published 2015 (version 13, Novenber 2023)
と印刷されているものを使っている。
なお、Pythonプログラミングの解説を岡本（2019）で行っている。
Imbens, G. W., & Rubin, D. B. (2015), Causal Inference for Statistics, Social, and Biomedical Sciences: An introduction. Cambridge University Press.