実験群と統制群の平均値の差の分析例
共変量の分布に重なりが無い場合
実験群と統制群の平均値の差の分析において、傾向スコア（propensity score）による共変量の分布の調整が説明されている（Imbens & Rubin, 2015）。しかし、2群における共変量の分布に重なりが無いときは、傾向スコアによる調整はできない。共変量の分布に重なりが無いデータの例として、図１のものを考える。
図１
Unitはデータのユニット番号、Yは差を調べる変量、Wは条件（0：統制群、1：実験群）、Xは共変量である。
共変量Xのヒストグラムは、図２のようである。
図２
統制群（W＝0）における共変量Xの分布と実験群（W＝1）における共変量Xの分布に重なりがないことが分かる。
統制群と実験群のXとYの散布図を図３に示す。
図３
統制群（W＝0、青丸）と実験群（W＝1、橙色丸）においてYの値の平均値に差が無いようであるが、Yの値は共変量Xの値が増加すると統制群および実験群において減少している。したがって、共変量Xの全域においてYはXの増加とともに減少している可能性がある。Yと共変量Xとの関係を考えて統制群と実験群の平均値の差を検討する必要がある。
以下では、まず、統制群と実験群の比較を直接行い、続いて共変量とYとの関係を調べる。最後に、共変量の影響を考慮した統制群と実験群の比較を行う。
なお、Jupyter Notebook用のファイルも用意した。Jupyter Notebook用では最小2乗法により回帰分析を行っている。NotebookのファイルとデータファイルをまとめてScript_JN.zipとした。ダウンロード展開すれば自由に用いることができる。また、Jupyter Notebookでのhtmlファイル出力mvregression.htmlも用意した。クリックすれば表示される。
統制群と実験群の比較を、モデル（１）により行う。
ここで、YとWは図１のデータにおけるYとWであり、比較対象の変量Yと条件を表す変量（W＝0：統制群、W＝1：実験群）である。Eは残差（誤差）項である。Wの係数ｔは統制群に対して実験群の平均値がどれだけ大きいかを表し、aは統制群におけるＹの平均値である。Ｅは平均０の正規分布に従う確率変数である。
モデル(1)によるデータ（図１）のベイズ分析を行うスクリプトをリスト１のように用意した。
リスト１のスクリプトを実行すると、MCMCサンプルのトレース図が図W1のように表示される。
図W1
パラメータｔの事後分布は、図Ｗ２のようである。
図W2
中央値が0.02、mad sdが11.21であるので、事後分布の中ほどに０があり、ｔの値は０と見なせるようである。
パラメータaとｔ、それぞれの事後分布の中央値を値とする回帰直線
のグラフを図Ｗ３に示す。
図W3
回帰直線はほぼ水平である。
以上より、統制群と実験群の測定値Ｙの平均値には差はないと言えそうであるが、共変量が無視されている。次に、共変量Ｘと測定値Ｙとの関係を調べた後、共変量も考慮して、統制群と実験群の差を分析する。
リスト１ モデル(1)によるベイズ分析スクリプト（PyStan 3対応）
import pandas as pd
import numpy as np
import stan
import arviz as az
import matplotlib.pyplot as plt
infl_nm = 'YWX_1.xlsx'
with pd.ExcelFile(infl_nm) as
f_xlsx:
data_xlsx =
pd.read_excel(f_xlsx)
print(data_xlsx.keys())
w_data = data_xlsx.values
print(w_data)
Y = np.array(data_xlsx['Y'],
dtype='float')
W = np.array(data_xlsx['W'],
dtype='float')
X = np.array(data_xlsx['X'],
dtype='float')
for v in zip(Y, W, X):
print(v)
stan_code = """
data {
int
N;
array[N] real Y;
array[N] real W;
}
parameters {
real
va;
real
vt;
real<lower=0.0> vsgm;
}
transformed parameters {
real
a;
real
t;
real
sgm;
a =
va * 1000.0;
t =
vt * 1000.0;
sgm
= 0.0001 + vsgm * 100.0;
}
model {
va ~
normal(0, 1);
vt ~
normal(0, 1);
vsgm
~ exponential(1.0);
for
(i in 1:N) {
Y[i] ~ normal(a + t*W[i], sgm);
}
}
"""
Data = {'N':len(Y), 'Y':Y, 'W':W}
sm = stan.build(stan_code, data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
az.plot_trace(infdata, var_names=['a',
't', 'sgm'])
plt.tight_layout()
plt.show()
dframe = fit.to_frame()
a_med = np.median(dframe['a'])
t_med = np.median(dframe['t'])
t_mad_sd = 1.483 *
np.median(np.abs(dframe['t'] - t_med))
az.plot_kde(dframe['t'])
plt.title(f'Med.={t_med:.2f}, mad sd ={t_mad_sd:.2f}', fontsize=18)
plt.xlabel('t', fontsize=18)
plt.show()
plt.scatter(W, Y)
plt.plot([W[0], W[9]], [a_med +
t_med*W[0], a_med + t_med*W[9]],
label=f'Y=a+tW')
plt.plot([], color='#ffffff',
label=f'a={a_med:.2f}(med.)\nt={t_med:.2f}(med.)')
plt.xticks([0, 1], [0, 1], fontsize=
16)
plt.xlabel('W: 0/control, 1/treated',
fontsize=16)
plt.ylabel('Y', fontsize=18)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
共変量Xと測定値Yとの関係を調べるために、次の単回帰モデルを設定する。
ここで、Eは残差（誤差）である。
モデル（2）により図１のデータのベイズ分析を行うために、リスト２のスクリプトを作成した。
リスト２を実行すると、図X1のトレース図を得る。
図X1
測定値Yと共変量Xの関係を表すパラメータｂの事後分布は図X2のようである。
図X2
中央値が−1.96で、mad sdが1.67であるので、共変量Xが増加すると観測値Yは減少していると言えそうである。
モデル（2）のパラメータaとbの中央値をパラメータの値とするときの回帰直線
のグラフを描くと図X3に示すものとなる。
図X3
共変量Xと測定値Yとの間の関係は無視できないと思われる。測定値Yの統制群と実験群の差を共変量Xも考慮して分析する。
リスト２ モデル（２）による単回帰分析（PyStan 3対応）
import pandas as pd
import numpy as np
import stan
import arviz as az
import matplotlib.pyplot as plt
infl_nm = 'YWX_1.xlsx' # input('Input
data file (*xlsx) = ')
with pd.ExcelFile(infl_nm) as
f_xlsx:
data_xlsx =
pd.read_excel(f_xlsx)
print(data_xlsx.keys())
w_data = data_xlsx.values
print(w_data)
Y = np.array(data_xlsx['Y'],
dtype='float')
W = np.array(data_xlsx['W'],
dtype='float')
X = np.array(data_xlsx['X'],
dtype='float')
for v in zip(Y, W, X):
print(v)
stan_code = """
data {
int
N;
array[N] real Y;
array[N] real X;
}
parameters {
real
va;
real
vb;
real<lower=0.0> vsgm;
}
transformed parameters {
real
a;
real
b;
real
sgm;
a =
va * 100.0;
b =
vb * 100.0;
sgm
= 0.0001 + vsgm * 100.0;
}
model {
va ~
normal(0, 1);
vb ~
normal(0, 1);
vsgm
~ exponential(1.0);
for
(i in 1:N) {
Y[i] ~ normal(a + b*X[i], sgm);
}
}
"""
Data = {'N':len(Y), 'Y':Y, 'X':X}
sm = stan.build(stan_code, data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
az.plot_trace(infdata, var_names=['a',
'b', 'sgm'])
plt.tight_layout()
plt.show()
dframe = fit.to_frame()
a_med = np.median(dframe['a'])
b_med = np.median(dframe['b'])
b_mad_sd = 1.483 *
np.median(np.abs(dframe['b'] - b_med))
az.plot_kde(dframe['b'])
plt.title(f'Med.={b_med:.2f}, mad_sd ={b_mad_sd:.2f}', fontsize= 18)
plt.xlabel('b', fontsize=16)
plt.tight_layout()
plt.show()
plt.scatter(X, Y)
plt.plot([X[0], X[9]], [a_med +
b_med*X[0], a_med + b_med*X[9]],
label ='Y=a+bX')
plt.plot([],
color="#ffffff", label=f'a={a_med:.2f}\nb={b_med:.2f}')
plt.xlabel('X', fontsize=16)
plt.ylabel('Y', fontsize=16)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
共変量Xを考慮した統制群（W＝０）と実験群（W＝１）の測定値Yの差を次の重回帰モデルで表す。
ここで、Eは残差（誤差）を表す。
モデル（3）によるベイズ分析を行うスクリプトをリスト３のように用意した。
リスト３を実行すると、MCMCサンプリング後、サンプリングのトレース図が表示される（図ＷＸ１）。
図WX１
パラメータｔの事後分布は、図ＷＸ２のようである。
図WX２
統制群と実験群の平均値の差を表すパラメータｔの事後分布の中央値が50.21、mad sdが0.78であるので、実験群の測定値Yは統制群より大きいと言えそうである。
パラメータｂの事後分布は、図ＷＸ３のようである。
図WX３
パラメータｂは共変量の影響を表すが、事後分布の中央値が−10.00、mad sdが0.13であるので、共変量Ｘが大きくなると測定値Ｙは減少すると言えそうである。
モデル（3）の回帰モデルをパラメータの事後分布の中央値をパラメータ値として、W＝０とW＝１の場合に分けて描くと図WX４のようになる。
図WX４
W＝０の場合の回帰直線
が赤の直線で描かれており、W＝１の場合の回帰直線
が橙色の直線で描かれている。
直線（W0）は、統制群のデータ（緑色の小円）によく合っており、直線（W1）は実験群のデータ（青色の小円）によく合っている。
赤色の星印は、統制群の回帰モデル（W0）において、実験群の共変量に対応する値を統制群（W＝０）であった場合の欠測値をimputeしたものである。同じ共変量に対する赤色の星印と青色の小円との縦方向の差が実験群の統制群に対する差（パラメータｔ）である。
リスト３ モデル（3）によるベイズ分析（PyStn 3対応）
import pandas as pd
import numpy as np
import stan
import arviz as az
import matplotlib.pyplot as plt
infl_nm = 'YWX_1.xlsx' # input('Input
data file (*xlsx) = ')
with pd.ExcelFile(infl_nm) as
f_xlsx:
data_xlsx = pd.read_excel(f_xlsx)
Y = np.array(data_xlsx['Y'],
dtype='float')
W = np.array(data_xlsx['W'],
dtype='float')
X = np.array(data_xlsx['X'],
dtype='float')
for v in zip(Y, W, X):
print(v)
IX = np.concatenate([np.ones((len(Y),
1)), np.reshape(W, (len(Y),1)),
np.reshape(X, (len(Y), 1))], axis=1)
beta = np.linalg.inv(IX.T @ IX) @ IX.T
@ np.reshape(Y, (len(Y), 1))
print(beta)
stan_code = """
data {
int
N;
array[N] real Y;
array[N] real W;
array[N] real X;
real
t_a;
real
t_t;
real
t_b;
}
parameters {
real
a;
real
t;
real
b;
real<lower=0.0> vsgm;
}
transformed parameters {
real
sgm;
sgm
= 0.0001 + vsgm;
}
model {
a ~
normal(t_a, 10);
t ~
normal(t_t, 10);
b ~
normal(t_b, 10);
vsgm
~ exponential(1.0);
for
(i in 1:N) {
Y[i] ~ normal(a + t*W[i] + b*X[i], sgm);
}
}
"""
Data = {'N':len(Y), 'Y':Y, 'W':W,
'X':X, 't_a':beta[0][0],
't_t':beta[1][0], 't_b':beta[2][0]}
sm = stan.build(stan_code, data=Data)
fit = sm.sample(num_samples=10000)
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary = az.summary(infdata, var_names=['a',
't', 'b', 'sgm'])
print('\nSummary...\n', summary)
az.plot_trace(infdata, var_names=['a',
't', 'b', 'sgm'])
plt.tight_layout()
plt.show()
dframe = fit.to_frame()
a_med = np.median(dframe['a'])
t_med = np.median(dframe['t'])
b_med = np.median(dframe['b'])
t_mad_sd = 1.483 *
np.median(np.abs(dframe['t'] - t_med))
b_mad_sd = 1.483 *
np.median(np.abs(dframe['b'] - b_med))
az.plot_kde(dframe['t'])
plt.title(f'med.={t_med:.2f}, mad sd ={t_mad_sd:.2f}', fontsize=18)
plt.xlabel('t', fontsize=18)
plt.show()
az.plot_kde(dframe['b'])
plt.title(f'med.={b_med:.2f}, mad sd ={b_mad_sd:.2f}', fontsize=18)
plt.xlabel('b', fontsize=16)
plt.show()
plt.scatter(X[:5], Y[:5], color='g',
label='W=0')
plt.scatter(X[5:], Y[5:], color='b',
label='W=1')
plt.plot([X[0], X[9]], [a_med +
b_med*X[0], a_med + b_med*X[9]],
color='magenta', label=f'W=0,
b={b_med:.2f}')
plt.plot([X[0], X[9]], [a_med + t_med
+ b_med*X[0],
a_med + t_med + b_med*X[9]],
color='orange', label=f'W=1, t={t_med:.2f}')
Y_imp = [a_med + b_med*v for v in
X[5:]]
plt.scatter(X[5:], Y_imp, marker='*',
color='r', label='W=0(imputed)')
plt.xlabel('X', fontsize=16)
plt.ylabel('Y', fontsize=16)
plt.legend(fontsize=12)
plt.show()
Imbens G. W. and Rubin, D. B. (2015). Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction. Cambridge University Press.