Up

Chap.8 Model-Based Inference for Completely Randomized Experiments

8.5 Simulation methods in the model-based approach

 

Imbens & Rubin (2015)の8.5節では、簡単なシミュレーションが説明されている。このシミュレーションを行うPythonスクリプトをリスト１のように作成した。

このシミュレーションでは、事後分は解析的に与えられたもの（式（8.27））が用いられているが、事後分布をMCMCサンプリングで推定するStanスクリプトも試みた。このサイトの後半で説明する。

まず、「8.4.3 An Analytic Example with Six Units」の説明に従って分析する。分析データは、Imbens & Rubin (2015)のTable 8.2に示されているものである。Unit数は６人で、統制群３人（unit 1, unit 3, unit 5），実験群３人（unit 2, unit 4, unit 6）である。

実験条件への振り分けをWで表すと、

W=[0, 1, 0, 1, 0, 1]

である。

unit iの条件jにおける結果（観測値の場合も欠測値の場合も含める）をYi(j)で表すと、Table 8.2に示されているデータは以下のように書くことができる。

統制群：

　[Y1(0), Y3(0), Y5(0)] = [0, 12.4, 0]

実験群：

　[Y2(1), Y4(1), Y6(1)] = [9.9, 3.6, 24.9]

このとき、観測値は

　[Y1obs, Y2obs, Y3obs, Y4obs, Y5obs, Y6pbs] = [Y1(0), Y2(1), Y3(0), Y4(1), Y5(0), Y6(1)]

であり、欠測値は

　[Y1mis, Y2mis, Y3mis, Y4mis, Y5mis, Y6mis] = [Y1(1), Y2(0), Y3(1), Y4(0), Y5(1), Y6(0)]

である。簡単のために、次の確率モデルが設定される（Imbens & Rubin, 2015, p. 156）。

このとき、パラメータの事後分布は次式で与えられる（Imbens & Rubin, 2015, p. 160）。

これは、リスト１では次のコードで表されている。

 

    mus = nr.default_rng().multivariate_normal([4.1, 12.8],

                                               [[5.8**2, 0], [0, 4.6**2]])

 

欠測値の事後予測は次式で与えられる（Imbens & Rubin, 2015, p. 158）。

これは、リスト１のコードでは次のようになっている。

 

cov = np.diagflat([64,100,64,100,64,100])

Ys_miss = nr.default_rng().multivariate_normal([mus[1], mus[0], mus[1],

                                                    mus[0], mus[1], mus[0]],

                                                    cov)

 

このとき、統制群と実験群の差は次式で表される（Imbens & Rubin, 2015, pp. 146-147）。

実験群と統制群の差が、観測値Yobsiと欠測値Ymisiの差で計算されていることに注意。

 

リスト１のスクリプトを実行すると、シミュレーションの終了後、のヒストグラムが表示される（図１）。

図１

 

図１のフォームを閉じると、プログラムの実行終了となり、端末にの平均値と標準偏差が表示される。

 

mean_tau = 8.693083664754294

STD = 5.232783490133446

 

この値は、解析的に得られる値、8.7と5.2（p. 162, Eq. (8.32)）、に十分に近い

 

 

 

リスト１　シミュレーションスクリプト例

 

import numpy as np

import numpy.random as nr

import matplotlib.pyplot as plt

 

Ys_obs = np.array([0, 9.9, 12.4, 3.6, 0, 24.9])

W = np.array([0, 1, 0, 1, 0, 1])

cov = np.diagflat([64,100,64,100,64,100])

NR = 100_000

taus = []

for r in range(NR):

    mus = nr.default_rng().multivariate_normal([4.1, 12.8],

                                               [[5.8**2, 0], [0, 4.6**2]])

    Ys_miss = nr.default_rng().multivariate_normal([mus[1], mus[0], mus[1],

                                                    mus[0], mus[1], mus[0]],

                                                    cov)

    tau_fs = np.dot((2*W - 1), (Ys_obs - Ys_miss)) / 6

    taus.append(tau_fs)

    if r % 1000 == 0:

        print(np.mean(taus), f'  r = {r:,}/{NR:,}')

 

mean_tau = np.mean(taus)

print('mean_tau =', mean_tau)

u_sd = np.std(taus, ddof = 1)

print('STD =', u_sd)

 

plt.hist(taus)

plt.xlabel('$\\tau_{fs}$', fontsize=18)

plt.show()

 

 

 

事後分布のMCMCサンプリングによる推定

 

パラメータとの事後分布を解析的に求めるのではなく（p. 156）、MCMCサンプリングで求める。このためのStanスクリプトをリストML1のように用意した。

 

 

リストML1　MCMCサンプリングのためのStanスクリプト（mcmc.stan）

 

data {

    int N;

    array[N] real Y0obs;

    array[N] real Y1obs;

}

transformed data {

    array[2*N] real W;

    array[2*N] real Yobs;

    for (i in 1:N) {

        W[2*i-1] = 0;  W[2*i] = 1;

        Yobs[2*i-1] = Y0obs[i];

        Yobs[2*i] = Y1obs[i];

    }

}

parameters {

    real vmu_c;   real vmu_t;

}

transformed parameters {

    real mu_c;  real mu_t;

    //   Rescaling (Imbens & Rubin (2015), Eq.(8.18))

    mu_c = vmu_c * 10;

    mu_t = vmu_t * 10; 

}

model {

    vmu_c ~ normal(0, 10);  vmu_t ~ normal(0, 10);

    //   Imbens & Rubin (2015), Eq.(8.17

    for (i in 1:N) {

        Y0obs[i] ~ normal(mu_c, 10);

        Y1obs[i] ~ normal(mu_t, 8);

    }

}

generated quantities {

    array[2*N] real Ymiss;

    real tau;

    for (i in 1:N) {

        Ymiss[2*i-1] = normal_rng(mu_t, 8);

        Ymiss[2*i] = normal_rng(mu_c, 10);

    }

    //     Imbens & Rubin (2015), Eq.(8.8)

    tau = 0.0;

    for (i in 1:2*N)  tau += (2*W[i] - 1) * (Yobs[i] - Ymiss[i]);

    tau /= (2*N);

}

 

 

 

平均treatment効果は、次式で与えている（Eq. (8.8)）。

データは、統制群のデータを奇数番目に、実験群のデータを偶数番目に置いている（リストML2）。これは、Imbens & Rubin (2015)のTable 8.2に合わせたものである。

 

 

リストML2　リストML1のStanスクリプトを読み込んでベイズ分析を行うプログラム

 

import numpy as np

import stan

import arviz as az

 

Ys_obs = np.array([0, 9.9, 12.4, 3.6, 0, 24.9])

W = np.array([0, 1, 0, 1, 0, 1])

 

idx = np.array([True, False]*3)

Y0obs =Ys_obs[idx]

Y1obs = Ys_obs[idx!=True]

print(Ys_obs)

print(Y0obs)

print(Y1obs)

 

with open('mcmc.stan') as f_stan:

    sm = stan.build(f_stan.read(), data = {'N':len(Y0obs),

                                           'Y0obs':Y0obs, 'Y1obs':Y1obs})

   

fit = sm.sample()

i_d = az.from_pystan(posterior = fit, posterior_model = sm )

print('\nSummary =\n', az.summary(i_d)) 

 

 

次のコード

 

idx = np.array([True, False]*3)

Y0obs =Ys_obs[idx]

Y1obs = Ys_obs[idx!=True]

 

より、データをY0obsと実験群のデータY1obsに分けている。

分けられたデータを用いて、リストMＬ1のStanスクリプトを次のコードで呼び出している。

 

with open('CompRandomExp.stan') as f_stan:

    sm = stan.build(f_stan.read(), data = {'N':len(Y0obs),

                                           'Y0obs':Y0obs, 'Y1obs':Y1obs})

 

リストMＬ2はPyStan 3対応である。

実行すると以下の結果を得る。

 

 

Summary =

             mean      sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

vmu_c      0.414   0.574  -0.634    1.528      0.009    0.008    4290.0    2829.0    1.0

vmu_t      1.286   0.455   0.378    2.082      0.008    0.005    3529.0    2896.0    1.0

mu_c       4.143   5.739  -6.344   15.279      0.088    0.075    4290.0    2829.0    1.0

mu_t      12.857   4.555   3.783   20.821      0.077    0.055    3529.0    2896.0    1.0

Ymiss[0]  13.065   9.108  -3.687   29.979      0.142    0.101    4085.0    4084.0    1.0

Ymiss[1]   4.022  11.416 -16.680   25.947      0.179    0.128    4075.0    3665.0    1.0

Ymiss[2]  12.955   9.192  -4.749   29.920      0.150    0.106    3762.0    3658.0    1.0

Ymiss[3]   3.833  11.510 -17.488   25.499      0.180    0.130    4077.0    3914.0    1.0

Ymiss[4]  12.703   9.227  -3.832   30.913      0.149    0.106    3818.0    3359.0    1.0

Ymiss[5]   4.028  11.542 -16.629   27.251      0.182    0.133    4017.0    3735.0    1.0

tau        8.807   5.061  -0.783   17.979      0.079    0.056    4121.0    3774.0    1.0

 

 

tauの平均が8.807、標準偏差が5.061であるのは、解析的に与えられている値（Eq. (8.32)）、平均8.7、標準偏差5.2に合っているといえる。

ちなみに、MCMCサンプリングはシミュレーションであるので、実行結果は変動する。再度、実行すると以下の結果となった。tauの平均は8.700、標準偏差は5.251である。

 

 

Summary =

             mean      sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

vmu_c      0.427   0.574  -0.659    1.476      0.010    0.008    3087.0    2568.0    1.0

vmu_t      1.285   0.459   0.425    2.149      0.008    0.006    2917.0    2592.0    1.0

mu_c       4.274   5.741  -6.590   14.765      0.103    0.081    3087.0    2568.0    1.0

mu_t      12.852   4.593   4.251   21.488      0.085    0.061    2917.0    2592.0    1.0

Ymiss[0]  12.800   9.171  -5.167   29.314      0.149    0.106    3766.0    3925.0    1.0

Ymiss[1]   4.044  11.701 -18.192   26.213      0.194    0.146    3640.0    3596.0    1.0

Ymiss[2]  13.110   9.189  -4.393   30.041      0.154    0.109    3551.0    3789.0    1.0

Ymiss[3]   4.247  11.444 -17.859   24.954      0.190    0.134    3640.0    4035.0    1.0

Ymiss[4]  12.952   9.319  -4.521   30.216      0.148    0.107    3952.0    3930.0    1.0

Ymiss[5]   4.373  11.500 -16.729   26.191      0.185    0.131    3850.0    3794.0    1.0

tau        8.700   5.251  -0.573   19.027      0.087    0.062    3665.0    3598.0    1.0

 

 

 

なお、Pythonプログラミングの解説を岡本（2019）で行っている。

 

 

引用文献

Imbens, G. W., & Rubin, D. B. (2015), Causal Inference for Statistics, Social, and Biomedical Sciences: An introduction. Cambridge University Press.

 

 

 

Home