比率（確率）の差の分析
ロジスティック回帰分析
ロジスティック回帰分析による比率（確率）の差の分析を行うためのStan/Pythonスクリプトを書いてみた。ロジスティック回帰分析については、＜岡本安晴「心理学データ分析と測定」勁草書房、第４章＞で説明している。
従属変数が連続量
である場合の重回帰モデルは、次式で与えられる。
これに対して、従属変量が確率
である場合は、その変換値
を重回帰モデル
で表す。ここで、変換関数
はロジスティック回帰モデルでは次式
で与えられ、ロジット変換（logit transformation）と呼ばれている。このとき、
である。ここで、従属変数が次の確率で1あるいは0の値をとるとしている。
ロジスティックモデルによる分析を行うためのStanスクリプトを個人単位で分析するものと条件ごとにまとめられたグループ単位のデータを扱うものの２つを用意した。それぞれ順番に説明する。
個人単位のデータとして表１のものを例として考える。
|
表１．（合否、模試成績、体重）の人工データ。合否は合格を1、不合格を0で表している。 |
|
(1, 790, 67) (0, 626, 57) (1, 869, 53) (1, 688, 54) (0, 687, 53) (0, 614, 48) (1, 763, 55) (0, 744, 58) (0, 629, 59) (0, 660, 62) (0, 680, 62) (1, 651, 59) (0, 771, 59) (0, 787, 72) (0, 641, 60) (1, 706, 62) (1, 821, 60) (1, 840, 59) (0, 644, 60) (0, 595, 48) (1, 753, 49) (1, 697, 72) (1, 743, 56) (0, 593, 56) (1, 828, 72) (1, 726, 55) (0, 713, 67) (1, 728, 58) (1, 748, 51) (0, 521, 56) (1, 770, 71) (1, 751, 77) (1, 752, 59) (0, 615, 48) (1, 956, 65) (1, 698, 62) (0, 668, 50) (1, 781, 54) (1, 583, 51) (1, 768, 72) (0, 705, 66) (1, 769, 62) (1, 744, 67) (1, 826, 62) (1, 864, 57) (0, 646, 70) (1, 896, 75) (1, 879, 56) (0, 664, 46) (0, 594, 73) (1, 696, 68) (0, 630, 65) (0, 729, 60) (1, 930, 54) (1, 816, 74) (0, 645, 62) (1, 777, 68) (1, 809, 47) (1, 761, 51) (0, 740, 60) (1, 816, 50) (1, 660, 53) (0, 711, 64) (1, 681, 61) (1, 650, 71) (1, 913, 68) (1, 616, 63) (0, 743, 50) (0, 687, 53) (0, 564, 63) (1, 620, 59) (1, 724, 61) (0, 650, 67) (0, 631, 69) (0, 748, 60) (0, 612, 58) (0, 522, 48) (0, 699, 61) (0, 561, 68) (1, 879, 59) (1, 863, 60) (1, 761, 55) (1, 803, 52) (1, 791, 56) (1, 786, 63) (1, 738, 57) (0, 707, 63) (1, 680, 49) (1, 695, 48) (1, 873, 67) (1, 580, 53) (1, 799, 57) (0, 614, 46) (0, 607, 55) (1, 783, 54) (1, 728, 47) (1, 734, 55) (0, 667, 79) (0, 721, 57) (1, 734, 78) |
作成したStanスクリプトは以下のようである。
＝＝＝＝＝＝＝＝＝＝＝＝ LgstcRegInd.stan ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
data {
int N;
int p;
int<lower = 0, upper =
1> Y[N];
matrix[N, p] X;
}
parameters{
vector[p] b;
}
transformed parameters {
real prob[N];
real v;
for (i in 1:N){
v =
0.0;
for
(j in 1:p)
v = v + b[j] * X[i, j];
prob[i] = inv_logit(v);
}
}
model {
for (i in 1:N)
Y[i]
~ bernoulli(prob[i]);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のStanスクリプトを用いて表１の形式のデータを分析するPythonスクリプトを以下のように作成した。これらのスクリプトファイルとサンプルデータファイルはまとめて圧縮ファイルPrgIndvdl.zipとした。ファイル名PrgIndvdl.zipをクリックしてダウンロードすることができる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
from pystan import StanModel
import pickle
def calc_stats( d ):
data = sorted(d)
n = len(data)
p025 = data[int(0.025 * n)]
med = data[int(0.5 * n)]
p975 = data[int(0.975 * n)]
sum = 0.0
for v in data:
sum += v
mean = sum / float(n)
n_neg = 0
while data[n_neg] < 0.0:
n_neg += 1
if
n_neg >= n:
break
prop_neg = n_neg / n
n_pos = 0
while data[n - 1 - n_pos]
> 0.0:
n_pos += 1
if
n_pos >= n:
break
prop_pos = n_pos / n
return mean, med, p025,
p975, prop_neg, prop_pos
#
# Prepare
the input data file
#
fn_in = input("Input data file =
")
f_in = open(fn_in, "r")
#
# Set the
contents of the input data file in the object data_f
#
data_f = f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
fn_out = input("Output data file
= ")
f_out = open(fn_out, "w")
f_out.write("Input data file =
" + fn_in + "\n\n")
pos = 0
while True:
if len(data_f[pos]) > 0:
if
data_f[pos][0] == "/":
break
pos += 1
pos += 1
p = int(data_f[pos])
pos += 1
var_name = data_f[pos].split()
print('\nlen(var_name) = ',
len(var_name))
print('\np = ', p, '\nVar_Names = ', var_name)
f_out.write('\np = {0}\nVar_Names =
{1}'.format(p, var_name))
N = 0
Y = []
X = []
pos += 1
while not(data_f[pos][0] == '/'):
temp_str =
data_f[pos].split()
Y.append(int(temp_str[0]))
X.append([])
X[N].append(1)
for i in range(p):
X[N].append(float(temp_str[i+1]))
N += 1
pos += 1
f_out.write('\nData...\n')
for i in range(N):
f_out.write('{0:>5}: {1}'.format(i+1, Y[i]))
for j in range(p):
f_out.write('{0:>10}'.format(X[i][j+1]))
f_out.write('\n')
f_out.write('\n')
ib = []
for i in range(p+1):
ib.append(0.0)
def f_init():
return dict(b = ib)
Data = {'N': N, 'p': (p + 1), 'Y': Y,
'X': X}
fit = pystan.stan(file =
'LgstcRegInd.stan', data = Data,
init = f_init, seed = 999, iter = 10000, warmup = 2000,
pars = ['b'], n_jobs = 1)
print(fit)
b_smpls = fit['b']
# print('\nlen(b_smpls) = ',
len(b_smpls))
# print('\nb_smpls[0] = ', b_smpls[0])
n_smpls = len(b_smpls)
param_b = []
for i in range(p+1):
param_b.append([])
for i in range(n_smpls):
for j in range(p+1):
param_b[j].append(b_smpls[i][j])
mean, med, p025, p975, prop_neg,
prop_pos = calc_stats(param_b[0])
print('\nb[0]<const>:\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean,
med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write('\nb[0]<const>:\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean, med, p025,
p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
for i in range(p):
print('\nb[{0}]<{1}>:'.format(i+1, var_name[i+1]))
f_out.write('\nb[{0}]<{1}>:'.format(i+1,
var_name[i+1]))
mean, med, p025, p975,
prop_neg, prop_pos = calc_stats(param_b[i+1])
print(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean, med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write('\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean, med, p025, p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
f_out.write('\n')
while True:
pos += 1
while True:
if
len(data_f[pos]) > 0:
if data_f[pos][0] == '/':
break
pos
+= 1
pos += 1
ck = int(data_f[pos])
if ck < 0:
break
pos += 1
s_cntrst =
data_f[pos].split()
cntrst = []
for i in range(p):
cntrst.append(float(s_cntrst[i]))
print('\ncontrast = ', end =
'')
for i in range(p):
print(' {0:<.5}<{1}>'.format(cntrst[i],
var_name[i+1]), end = '')
print(' ')
f_out.write('\ncontrast = ')
for i in range(p):
f_out.write('
{0:<.5}<{1}>'.format(cntrst[i], var_name[i+1]))
f_out.write('\n')
c_b = []
for i in range(n_smpls):
v =
0.0
for
j in range(p):
v += cntrst[j] * param_b[j+1][i]
c_b.append(v)
mean, med, p025, p975,
prop_neg, prop_pos = calc_stats(c_b)
print(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean, med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean, med, p025, p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
f_out.close()
print('\n')
print(fn_out, ' was saved.')
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のスクリプトを実行すると、入力データファイル名と出力用テキストファイル名の設定が求められる（図１）。
図１
入力データファイルは、図２の形式で用意する。
・
・
・
図２
スラッシュ/で始まる最初の行までは、無視される（読み飛ばされる）。
このスラッシュ/で始まる最初の行の次の行に、独立変数の数qを書く。Pythonスクリプトではqの値は変数pに読み込まれ、モデル（１）あるいは（２）における確率はStanスクリプトでは配列propで表されている。Stanスクリプトでは、pの値は、Pythonスクリプトでのpの値より１つ増やされているが、これは定数項
に対する独立変数
がStanスクリプトを実行するときに追加されているからである。
その次の行には、変数のラベルを、従属変数、独立変数１、・・・、独立変数qの順に並べる。
このラベルを書いた行の次の行から、1行に1組のデータを、従属変数、独立変数１、・・・、独立変数qの値の順に書いていく。
すべての組のデータを書いた最後は、スラッシュ/で始まる行を置いて、データの終わりを表す。
次のスラッシュ/で始まる行の次の行に書かれている整数値が０以下の値であれば、プログラムの実行終了となる。
正の値が書かれているときは、その次の行に各独立変数に対するウェイトを並べる。例えば、
と書かれているときは、係数の重み付き和
の事後分布が求められる。ここで、回帰モデル
において、入力した独立変数の他に定数項と組み合わせる変数
がプログラム内で自動的に追加されていることに注意。
上の例の場合、ウェイト
1 0
に対しては
の事後分布が求められ、ウェイト
0 1
に対しては、
の事後分布が求められる。
入力データファイル名と出力用ファイル名を設定すると計算が始まる。出力用ファイル名は任意のテキストファイル名である。計算が終了すると、パラメータの事後分布の統計量が表示されてプログラムの実行終了となる（図３）。統計量は、ランダムサンプリングによる推定であるので、サンプリングの条件に依存した変動がある。他の方法（乱数の初期値を変えるなど）による結果と完全には一致しないのが普通である。
図３
プログラムの実行終了後、出力ファイルを開くと、表１のデータの場合、以下のような出力になっている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Input data file = DataI.txt
p = 2
Var_Names = ['Result', 'Score',
'Weight']
Data...
1: 1 790.0 67.0
2: 0 626.0 57.0
3: 1 869.0 53.0
4: 1 688.0 54.0
5: 0 687.0 53.0
・
・
・
95: 1 783.0 54.0
96: 1 728.0 47.0
97: 1 734.0 55.0
98: 0 667.0 79.0
99: 0 721.0 57.0
100: 1 734.0 78.0
b[0]<const>:
-12.988(mean) -12.818(med) 95%CI [-20.048, -6.7665]
Prop(v <= 0.0) =
1.0000 Prop(v
>= 0.0) = 0.0000
b[1]<Score>:
0.020731(mean) 0.020558(med) 95%CI [0.012799, 0.029819]
Prop(v <= 0.0) =
0.0000 Prop(v
>= 0.0) = 1.0000
b[2]<Weight>:
-0.022119(mean) -0.022008(med) 95%CI [-0.087451, 0.042867]
Prop(v <= 0.0) =
0.7508 Prop(v
>= 0.0) = 0.2492
contrast = 1.0<Score> 0.0<Weight>
0.020731(mean) 0.020558(med) 95%CI [0.012799, 0.029819]
Prop(v <= 0.0) =
0.0000 Prop(v
>= 0.0) = 1.0000
contrast = 0.0<Score> 1.0<Weight>
-0.022119(mean) -0.022008(med) 95%CI [-0.087451, 0.042867]
Prop(v <= 0.0) =
0.7508 Prop(v
>= 0.0) = 0.2492
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
まず、読み込んだ入力データファイル名とデータが出力されている。これにより、データが意図通り読み込まれたかどうかの確認ができる。
その後、各パラメータの事後分布の統計量が出力されている。＜Score＞の効果が有意であることが分かる。
パラメータごとの統計量に続いて、コントラストの統計量が出力されている。上の場合は、コントラストのウェイトは単一のパラメータのみを取り出すものなので、先に出力されている各パラメータの統計量と同じである。
上の表１のデータは、個人単位（ケース単位）のものであるが、データが条件ごとの頻度で与えられている場合も同様に分析できる。
|
表２ 数学・算数の授業法と数学・算数の好き嫌いの人数． |
|||
|
|
従来の教授法 |
新しい教授法 |
|
|
小学３年生 |
好き |
20 |
30 |
|
嫌い |
20 |
10 |
|
|
好きの比率 |
0.5 |
0.75 |
|
|
中学２年生 |
好き |
5 |
60 |
|
嫌い |
35 |
200 |
|
|
好きの比率 |
0.125 |
0.231 |
|
例えば、表２のデータは、数学（算数）の２種類の授業法（従来の方法と新しい方法）の効果を比較したもの（人工データ）である。従来の方法による授業を受けた場合、小学３年生では２０人が数学（算数）が好きと答え、２０人が嫌いと答えたが、新しい授業法を受けた生徒は、３０人が好きと答え、１０人が嫌いと答えている。中学２年生の場合は、従来の授業法を受けた場合は５人が数学が好きと答え、３５人が嫌いと答えたが、新しい授業法では６０人が好きと答え２００人が嫌いと答えている。好きと答えた者の比率は、小学３年生では従来の方法の場合は0.5であるのに対して、新しい方法の場合は0.75であり、中学２年生では従来の方法の場合は0.125であり、新しい授業法では0.231である。小学３年生、中学２年生とも、同じ学年同士の比較では、新しい授業法を受けた生徒の方が従来の授業法を受けた生徒より、数学・算数が好きと答えた者の比率が高い。このデータをロジスティック回帰モデルで分析してみる。
好きと答える確率を
とおき、ロジスティック回帰モデルを
とおく。ここで、
とおいた。データがどの条件のものであるかを示し、それらの効果を表すために導入された変数
、
、
はダミー変数と呼ばれている。ダミー変数については、＜岡本安晴「データ分析のための統計学入門」おうふう、pp.
130-133＞を参照されたい。
式（３）の各条件における値を表３に示す。
|
表３ ロジスティック回帰モデル（３）の値。 |
||
|
|
従来の授業法 |
新しい授業法 |
|
小学３年生 |
|
|
|
中学２年生 |
|
|
係数は小学３年生における従来の授業法を受けた場合の「好き」の比率に対する式（３）の値を表し、この条件での値を基準にして授業法の効果および学年の効果を表している。
は授業法の主効果、
は学年の主効果、
は交互作用である。
表3の形式のデータを分析するためのより一般的モデルを以下のように設定する。
式（３）をグループgにおける独立変数がq個のものとして、
とおく。
確率
の表す事象（表２の場合は、「好き」）のグループgにおける度数を
、確率の表す事象ではない事象（表２の場合は、「嫌い」）のグループgにおける度数を
で表すと、尤度は
で与えられる。式（４）に対するStanスクリプトを以下のように用意した。Stanスクリプトでは、グループ数がＮ、独立変数の数q+1がpで表されている。独立変数の数は、定数項に対応して変数
を加えてq+1である。
＝＝＝＝＝＝＝＝＝＝＝＝＝ LgstcRegGrp.stan ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
data {
int N;
int p;
int Y[N, 2];
matrix[N, p] X;
}
parameters{
vector[p] b;
}
transformed parameters {
simplex[2] prob[N];
real v;
for (i in 1:N){
v =
0.0;
for
(j in 1:p)
v = v + b[j] * X[i, j];
prob[i][1] = inv_logit(v);
prob[i][2] = 1.0 - prob[i][1];
}
}
model {
for (i in 1:N)
Y[i]
~multinomial(prob[i]);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のStanスクリプトを利用するPythonスクリプトを以下のように用意した。プログラムスクリプトファイルとサンプルデータファイルはまとめて圧縮ファイルPrgGrpd.zipとした。ファイル名PrgGrpd.zipをクリックしてダウンロードすることができる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
def calc_stats( d ):
data = sorted(d)
n = len(data)
p025 = data[int(0.025 * n)]
med = data[int(0.5 * n)]
p975 = data[int(0.975 * n)]
sum = 0.0
for v in data:
sum
+= v
mean = sum / float(n)
n_neg = 0
while data[n_neg] < 0.0:
n_neg += 1
if
n_neg >= n:
break
prop_neg = n_neg / n
n_pos = 0
while data[n - 1 - n_pos]
> 0.0:
n_pos += 1
if
n_pos >= n:
break
prop_pos = n_pos / n
return mean, med, p025,
p975, prop_neg, prop_pos
#
# Prepare
the input data file
#
fn_in = input("Input data file =
")
f_in = open(fn_in, "r")
#
# Set the
contents of the input data file in the object data_f
#
data_f = f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
fn_out = input("Output data file
= ")
f_out = open(fn_out, "w")
f_out.write("Input data file =
" + fn_in + "\n\n")
pos = 0
while True:
if len(data_f[pos]) > 0:
if
data_f[pos][0] == "/":
break
pos += 1
pos += 1
p = int(data_f[pos])
pos += 1
var_name = data_f[pos].split()
print('\nlen(var_name) = ',
len(var_name))
print('\np = ', p, '\nVar_Names = ',
var_name)
f_out.write('\np = {0}\nVar_Names =
{1}'.format(p, var_name))
N = 0
Y = []
X = []
pos += 1
while not(data_f[pos][0] == '/'):
temp_str =
data_f[pos].split()
Y.append([])
Y[N].append(int(temp_str[0]))
Y[N].append(int(temp_str[1]))
X.append([])
X[N].append(1)
for i in range(p):
X[N].append(float(temp_str[i+2]))
N += 1
pos += 1
f_out.write('\nData...\n')
for i in range(N):
f_out.write('{0:>5}: {1:<7} {2:<7} '.format(i+1, Y[i][0], Y[i][1]))
for j in range(p):
f_out.write('{0:>10}'.format(X[i][j+1]))
f_out.write('\n')
f_out.write('\n')
Data = {'N': N, 'p': (p + 1), 'Y': Y,
'X': X}
fit = pystan.stan(file =
'LgstcRegGrp.stan', data = Data, seed = 999, n_jobs = 1)
print(fit)
b_smpls = fit['b']
n_smpls = len(b_smpls)
param_b = []
for i in range(p+1):
param_b.append([])
for i in range(n_smpls):
for j in range(p+1):
param_b[j].append(b_smpls[i][j])
mean, med, p025, p975, prop_neg,
prop_pos = calc_stats(param_b[0])
print('\nb[0]<const>:\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean,
med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write('\nb[0]<const>:\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean,
med, p025, p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
for i in range(p):
print('\nb[{0}]<{1}>:'.format(i+1, var_name[i+2]))
f_out.write('\nb[{0}]<{1}>:'.format(i+1, var_name[i+2]))
mean, med, p025, p975,
prop_neg, prop_pos = calc_stats(param_b[i+1])
print(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean, med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write('\n {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean, med, p025, p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
f_out.write('\n')
while True:
pos += 1
while True:
if
len(data_f[pos]) > 0:
if data_f[pos][0] == '/':
break
pos
+= 1
pos += 1
ck = int(data_f[pos])
if ck < 0:
break
pos += 1
s_cntrst =
data_f[pos].split()
cntrst = []
for i in range(p):
cntrst.append(float(s_cntrst[i]))
print('\ncontrast = ', end =
'')
for i in range(p):
print('
{0:<.5}<{1}>'.format(cntrst[i], var_name[i+2]), end = '')
print(' ')
f_out.write('\ncontrast = ')
for i in range(p):
f_out.write(' {0:<.5}<{1}>'.format(cntrst[i],
var_name[i+2]))
f_out.write('\n')
c_b = []
for i in range(n_smpls):
v =
0.0
for
j in range(p):
v += cntrst[j] * param_b[j+1][i]
c_b.append(v)
mean, med, p025, p975, prop_neg,
prop_pos = calc_stats(c_b)
print(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]'.
format(mean, med, p025, p975))
print(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}'.format(prop_neg, prop_pos))
f_out.write(' {0:>.5}(mean) {1:>.5}(med) 95%CI [{2:>.5}, {3:>.5}]\n'.
format(mean, med, p025, p975))
f_out.write(' Prop(v <= 0.0) =
{0:<.4f}
Prop(v >= 0.0) = {1:<.4f}\n'.
format(prop_neg, prop_pos))
f_out.close()
print('\n')
print(fn_out, ' was saved.')
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のPythonスクリプトを実行すると入力データファイル名と出力用ファイル名の設定が求められる（図４）。
図４
入力データファイルは、表２のデータの場合、図５のように用意する。
図５
スラッシュ/で始まる最初の行までは無視される（読み飛ばされる）。
このスラッシュ/で始まる最初の行の次の行に独立変数の数qを書く。
その次の行には、変数のラベルを、従属変数（「好き」のデータ数、すなわち確率で生起する事象の度数を表すもの）、従属変数（「嫌い」のデータ数、すなわち確率
で生起する事象の度数を表すもの）,独立変数１、・・・、独立変数qの順に並べる。
このラベルを書いた行の次の行から、1行に1組のデータを、従属変数、従属変数、独立変数１、・・・、独立変数qの値の順に書いていく。
すべての組のデータを書いた最後は、スラッシュ/で始まる行を置いて、データの終わりを表す。
次のスラッシュ/で始まる行の次の行に書かれている整数値が０以下の値であれば、プログラムの実行終了となる。
正の値が書かれているときは、その次の行に各独立変数に対するウェイトを並べる。例えば、
と書かれているときは、係数の重み付き和
の事後分布が求められる。ここで、回帰モデル
において、入力した独立変数の他に定数項と組み合わせる変数がPythonスクリプト内で自動的に追加されていることに注意。Pythonスクリプトでは、次のモデルが用いられている。
上の入力データファイルの例における重み
1 0 1
に対応する重み付きパラメタの和は
を表す。これは、中学２年生における新しい授業法と従来の授業法の差
を表すものである（表３参照）。
入力データファイル名と出力用ファイル名を設定すると計算が始まる。出力用ファイル名は任意のテキストファイル名である。計算が終了すると、Stanのサンプリング結果の統計量などが出力される（図６）。統計量は、ランダムサンプリングによる推定であるので、サンプリングの条件に依存した変動がある。他の方法（乱数の初期値を変えるなど）による結果と完全には一致しないのが普通である。
図６
プログラムの実行終了後、出力ファイルを開くと表２のデータの場合は、図７のような出力なっている。
図７
上で説明した中学２年生における新しい授業法と従来の授業法の差
を表すコントラストの事後分布をみるとぎりぎりで有意差があると解釈できるが、結論を出すにはさらに慎重なデータ収集と分析が必要であると解釈されるものである。