P値と事後分布
帰無仮説検定のP値とベイズ分析の事後分布を比較してみた。
ｐ値の使用について、American Statistical Association (ASA)から次の声明が出されている。
The ASA Statement on p-Values: Context, Process, and Purpose
Ronald L. Wasserstein & Nicole A. Lazar, 2016, pp. 129-131.
ASA Statement on Statistical Significance and P-Values
Ronald L. Wassersteina, 2016, pp. 131-133.
これらは
https://doi.org/10.1080/00031305.2016.1154108
でアクセスすることができる。Pp.129-131までが前半、その後にpp. 131-133が掲載されている。
上のASA statementでは、P値はuseful statistical measureであり得るが、誤用と誤解があると論じられている（p. 131）。
Natureの2019年の論文
"Scientists rise up against statistical significance"
V. Amrhein et al., Nature 567, 305-307 (2019)
https://www.nature.com/articles/d41586-019-00857-9
では、ｐ値のみによる仮説の正否・２者択一的決定は問題であり、ｐ値以外の指標を用いた場合も、この分析結果のカテゴリ的決定に伴う問題は起こり得ると指摘されている（p. 306）。
上のASAの論文は、いずれも３ページであり、Natureの論文も３ページである。
ここでは、ｐ値の問題を念頭に置きながら、ｐ値と事後分布のシミュレーションを行ってみた。
正規分布に従う２つの母集団の平均値の差の検定を扱う。母集団の平均値の差が小さく、実質的に差が無いと見なされる場合と、差が十分に大きい場合について考える。使用したプログラムは圧縮ファイルprgs.zipにまとめたので、自由にダウンロードして利用されたい。また、PythonでのStanの使用については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章で説明している。
なお、効果量と検定力の関係をシミュレーションで調べたものをこのウェブサイトに掲載している。
２つの母集団においてその平均値の差が小さく、実質的には差がないと考えられる場合を想定する。シミュレーションでは、母集団の平均値は、103と100であり、それぞれ標準偏差10の正規分布に従っているとする。ここで、平均値の差は、5以下であれば実質的な差はないとして扱われる場合を考える。２つの母集団の平均値の差は103−100＝3であるので、実質的には差は無いと考えられる場合である。ｐ値に及ぼす標本サイズの影響を調べるため、各母集団からの標本サイズとしてN＝5とN＝200の場合についてサンプリングを行い、ｔ検定における両側検定のｐ値を求めた。N＝５およびN＝200の場合とも1000回のシミュレーションである。Pythonスクリプトをリスト１に示す。
リスト１ ｐ値のシミュレーション。母集団の差が小さい場合。
import numpy as np
import numpy.random as nr
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
mu1 = 103.0
mu2 = 100.0
sgm = 10.0
n_trials = 1000
n_samples = [5, 200]
p_vals = []
for s_size in range(2):
p_vals.append([])
n = n_samples[s_size]
for t in range(n_trials):
x1 =
nr.normal(loc = mu1, scale = sgm, size = n)
x2 =
nr.normal(loc = mu2, scale = sgm, size = n)
results = ss.ttest_ind(x1, x2)
p_vals[s_size].append(results.pvalue)
X = []
G = []
for s_size in range(2):
X += p_vals[s_size]
G +=
['{}'.format(n_samples[s_size])] * len(p_vals[s_size])
p_05 = []
for s_size in range(2):
p =
[np.array(p_vals[s_size]) < 0.05]
p_05.append(np.array(p).mean())
sb.boxplot(x = 'G', y = 'X', data =
{'X': X, 'G': G}, order = ['5', '200'])
plt.title(('Prop(p-val < 0.05) =
{0:.3}: N = {1:} \n' +
'Prop(p-val < 0.05) = {2:.3}: N = {3:}').format
(p_05[0], n_samples[0], p_05[1], n_samples[1]), fontsize = 16)
plt.xlabel('N', fontsize = 16)
plt.ylabel('Proportion', fontsize = 14)
plt.show()
上のスクリプトを実行すると図1.1の結果を得る。
図1.1
各母集団からN個のサンプルを採り、ｐ値を求め箱図を描いたものである。サンプル数がN＝５と少ないときは、ｐ値が0.05以下になる比率は１割以下であるが、サンプル数をN＝200に増やすとｐ値が0.05以下になる比率は８割を超えている。ｐ値だけで母集団の平均値に差があるかどうかを決めると、Ｎ＝５のときは１割以下の場合で差があると結論することになり、Ｎ＝200のときは８割以上の場合で差があると結論することになる。
次にベイズ分析における事後分布のシミュレーションを行う。
母集団からのサンプリングの条件は、ｐ値の場合と同じであるが、シミュレーションの回数は100回とした。プログラムをリスト２に示す。
リスト２ 事後分布のシミュレーション。差が小さい場合。
import numpy as np
import numpy.random as nr
import scipy.stats as ss
import pystan
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
StanCode = """
data {
int N;
real x1[N];
real x2[N];
}
parameters {
real mu1;
real mu2;
real<lower = 0.0> sgm;
}
model {
mu1 ~ normal(0.0, 1000.0);
mu2 ~ normal(0.0, 1000.0);
sgm ~ exponential(0.001);
for (i in 1:N){
x1[i] ~ normal(mu1, sgm);
x2[i] ~ normal(mu2, sgm);
}
}
"""
sm = pystan.StanModel(model_code =
StanCode)
with open('sm.pkl', 'wb') as f:
pickle.dump(sm, f)
mu1 = 103.0
mu2 = 100.0
sgm = 10.0
n_trials = 100
s_size = [5, 200]
ddata = dict()
X = []
G = []
for v_size in s_size:
p_eff = []
for t in range(n_trials):
if t
% 10 == 0:
print('Trial {} started.'.format(t))
x1 =
nr.normal(loc = mu1, scale = sgm, size = v_size)
x2 =
nr.normal(loc = mu2, scale = sgm, size = v_size)
Data
= {'N': len(x1), 'x1': x1, 'x2': x2}
sm =
pickle.load(open('sm.pkl', 'rb'))
fit
= sm.sampling(data = Data,n_jobs = 1)
smpl_mu1 = fit['mu1']
smpl_mu2 = fit['mu2']
smpl_mu12 = smpl_mu1 - smpl_mu2
tmp
= np.abs(smpl_mu12) < 5
p_eff.append(tmp.mean())
if t
< 10:
sb.kdeplot(smpl_mu12)
X += p_eff
G += ['{}'.format(v_size)] *
len(p_eff)
plt.title('Posterior
Distributions\n' +
'N = {}'.format(v_size), fontsize = 16)
plt.xlabel('$\mu_1$ -
$\mu_2$', fontsize = 14)
plt.yticks([])
plt.show()
sb.boxplot(x = 'G', y = 'X', data =
{'G':G, 'X':X}, order = ['5', '200'])
sb.stripplot(x = 'G', y = 'X', data =
{'G':G, 'X':X}, order = ['5', '200'], color = 'k')
plt.xlabel('N', fontsize = 14)
plt.ylabel('P', fontsize = 14)
plt.title('Prob(|$\mu_1$ - $\mu_2$|
< 5)\nfor posterior distributions', fontsize = 16)
plt.show()
上のスクリプトを実行すると、まず各母集団に対するサンプル数Ｎ＝５の場合におけるシミュレーションが100回行われ、最初の10回のシミュレーションにおける平均値パラメータの差
の事後分布のグラフが描かれる（図1.2）。母集団平均値103である母集団に対する平均値パラメータが
であり、平均値100である母集団に対する平均値パラメータが
である。平均値の差は５以下は無視できる場合を考えているが、図1.2における事後分布を見ると、５以上の領域にもかなりの分布が見られる。
図1.2
図1.2のウィンドウを閉じると、サンプル数Ｎ＝200の場合のシミュレーションが始まり、100回のシミュレーションが終了すると、図1.3のグラフが表示される。
図1.3
図1.3には10回分のシミュレーションにおける事後分布が表示されているが、かなりの分布が区間［0，5］の範囲にあることが分かる。すなわち、母集団平均のパラメータの差は実質的には無視できる範囲にある。
図1.3のウィンドウを閉じると、図1.4のグラフが表示される。
図1.4
平均値の差は、５以下であれば無視できる場合を想定しているので、事後分布において差が５以下の確率
を求めた。サンプル数Nの値別に、各サンプルに対する事後分布の式（１）の値を求め、箱図と紐図を描いたものが図1.4である。サンプル数N＝５のときは、式（１）の値が小さい、すなわち差が無視できる範囲にある確率が0.5以下のものがかなり認められるが、サンプル数がN＝200の場合は、差が無視できる範囲である確率が0.9以上のものが多くなっている。
２つの母集団においてその平均値に無視できない差がある場合を想定する。シミュレーションにおける母集団の平均値は、115と100であり、それぞれ標準偏差10の正規分布に従っているとする。ここで、平均値の差は、5以下であれば実質的な差はないとして扱われる場合を考える。２つの母集団の平均値の差は115−100＝15であるので、実質的な差は十分にあると考えられる。ｐ値に対する標本サイズの影響を調べるため、各母集団からの標本サイズとしてN＝5とN＝200の場合についてサンプリングを行い、ｔ検定を行い両側検定のｐ値を求めた。サンプリングの回数は、N＝５およびN＝200の場合とも1000回のシミュレーションである。Pythonスクリプトはリスト１において母集団のパラメータを
mu1 = 115.0
mu2 = 100.0
sgm = 10.0
としたものである。このスクリプトも圧縮ファイルprgs.zipに含まれている。
スクリプトを実行すると図2.1の結果を得る。
図2.1
各母集団からN個のサンプルを採り、ｐ値を求めて箱図を描いたものである。サンプル数がN＝５と少ないときでも、ｐ値が0.05以下になる比率は５割ほどある。サンプル数をN＝200に増やすとｐ値は確実に0.05以下になるといえる。
次にベイズ分析における事後分布のシミュレーションを行う。
母集団からのサンプリングの条件は、ｐ値の場合と同じであるが、シミュレーションの回数は100回とした。プログラムはリスト２において母集団パラメータを
mu1 = 115.0
mu2 = 100.0
sgm = 10.0
と設定したものである。プログラムは圧縮ファイルprgs.zipに含まれている。
このスクリプトを実行すると、まず各母集団に対するサンプル数Ｎ＝５の場合におけるシミュレーションが100回行われ、最初の10回のシミュレーションにおける平均値パラメータの差の事後分布のグラフが描かれる（図2.2）。母集団平均値115である母集団に対する平均値パラメータがであり、平均値100である母集団に対する平均値パラメータがである。平均値の差は５以下は無視できる場合を考えているが、図2.2における事後分布を見ると、５以上の領域にかなり分布しているといえる。
図2.2
図2.2のウィンドウを閉じると、サンプル数Ｎ＝200の場合のシミュレーションが始まり、100回のシミュレーションが終了すると、図2.3のグラフが表示される。
図2.3
図2.3には10回分のシミュレーションにおける事後分布が表示されているが、すべての分布が区間［0，5］より右側に離れていることが分かる。すなわち、母集団平均のパラメータの差は実質的に十分に認められるものであるといえる。
図2.3のウィンドウを閉じると、図2.4のグラフが表示される。
図2.4
平均値の差は、５以下であれば無視できる場合を想定しているので、事後分布において差が５以下の確率
を求めた。サンプル数Nの値別に、各サンプルに対する事後分布の式（１）の値を求め、箱図と紐図を描いたものが図2.4である。サンプル数N＝５のときは、式（１）の値は0.1ぐらいで無視出来ないが、サンプル数がN＝200の場合は、式（１）の値は０、すなわち差が無視できる範囲には平均値パラメータ値の差は存在しないといえる。