効果量と検定力

５　効果量と検定力

効果量δおよび効果量と検定力の関係を調べるシミュレーションプログラム例を用意した。サンプルプログラムファイルは圧縮ファイルchap5_files.zipに用意したので、ダウンロード・展開すれば自由に利用できる。

サンプルプログラムは、Jupyter Notebook用にまとめて圧縮ファイルchap5_prgs.zipとした。ダウンロード展開して自由に利用されたい。

効果量δと検定力

効果量δと検定力の関係を調べるシミュレーションを本書のアルゴリズム5.1.1によって行うPythonスクリプト例をリスト１のように作成した。リスト１のスクリプトを実行すると、計算結果が図１のように表示される。

図１　効果量δと検定力

図１はそれぞれデータ数MおよびNの２条件での平均値の差の検定を行う場合の検定力（Power）を調べたものである。効果量が中ほどとされるの条件では、検定力が0.5より大きく、すなわち1／2以上の確率で有意差が得られるためには、データ数MおよびNは20より大きくなければならないことが分かる。

リスト１　効果量δと検定力を関係を調べるPythonスクリプト例（s5_effect_size_delta.py）

# アルゴリズム5.1.1によるシミュレーション

import matplotlib.pyplot as plt

import numpy as np

import scipy.stats as ss

NSamples = 10_000

L_MNs = [5, 10, 20, 50, 100, 200, 500]

L_deltas = [0.2, 0.5, 0.8]

v_alpha = 0.05

powers = np.empty((len(L_deltas), len(L_MNs)))

for i_delta, v_delta in enumerate(L_deltas):

for j_M, M in enumerate(L_MNs):

N = M

df = M + N - 2

t_a = ss.t.ppf(0.95, df) # 限界値

count, i = 0, 1 # 初期化

while True:

v1 = ss.norm.rvs(loc = v_delta * ((M * N) / (M + N)) ** 0.5, scale = 1)

v2 = ss.chi2.rvs(df)

t_v = v1 / ((v2 / df) ** 0.5)

if t_a < t_v: # 有意？

count += 1

i += 1

if i > NSamples: # 設定条件でのシミュレーション終了？

break

print('i_delta =', i_delta, ' delta =', v_delta, ' j_M =', j_M, ' M = N =', M)

print('Power =', count / NSamples)

powers[i_delta][j_M] = count / NSamples # 検定力のシミュレーション推定値

print('Powers =\n', powers)

# シミュレーション結果の描画

v_log_x = np.log(L_MNs)

L_log_x = ['{}'.format(v) for v in L_MNs]

for j in range(len(L_deltas)):

plt.plot(v_log_x, powers[j], label = '$\delta$ = {}'.format(L_deltas[j]))

plt.xticks(v_log_x, L_log_x)

plt.xlabel('M = N', fontsize = 16)

plt.yticks([0, 0.5, 0.95, 1])

plt.ylabel('Power', fontsize = 16)

plt.hlines(1.0, np.log(5), np.log(500), lw = 1)

plt.hlines(0.95, np.log(5), np.log(500), lw = 1)

plt.hlines(0.5, np.log(5), np.log(500), lw = 1)

plt.hlines(0, np.log(5), np.log(500), lw = 1)

plt.legend(loc = 'lower right')

plt.title('Effect Size $\delta$, Sample Size M, and Power', fontsize = 18)

plt.tight_layout()

plt.show()

効果量とと検定力

効果量と検定力の関係を、本書のアルゴリズム5.2.1によるシミュレーションで調べるPythonスクリプトをリスト２のように作成した。リスト２を実行すると、計算結果が図２のように表示される。

図２　効果量と検定力の関係

効果量の大きさが中ほどの値とされるのときは、検定力が0.5以上、すなわち1／2以上の確率で有意な結果を得るためには、各水準内のデータ数Nは20以上でなければならないことが判る。

リスト２　効果量と検定力の関係を調べるシミュレーションスクリプト（s5_effect_size_f.py）

# アルゴリズム5.2.1によるシミュレーション

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats as ss

NSamples = 10_000

L_Ns = [5, 10, 20, 50, 100, 200, 500]

L_f = [0.1, 0.25, 0.40]

v_alpha = 0.05 # 有意水準

a = 3 # ３水準

powers = np.empty((len(L_f), len(L_Ns)))

for i_f, v_f in enumerate(L_f):

for j_N, N in enumerate(L_Ns):

df_A = a - 1

df_S_A = a * (N -1)

lmbd = a * N * (v_f ** 2) # lambdaの値

F_a = ss.f.ppf(0.95, df_A, df_S_A) # 限界値

count, i = 0, 1 # 初期値設定

while True:

v1 = ss.ncx2.rvs(df_A, lmbd)

v2 = ss.chi2.rvs(df_S_A)

F_v = (v1/df_A) / (v2/df_S_A)

if F_a < F_v: # 有意？

count += 1

i += 1

if i > NSamples: # 設定条件でのシミュレーション終了？

break

power = count / NSamples # 検定力

powers[i_f][j_N] = power

print('i_f =', i_f, ' j_N =', j_N)

print('power = ', power)

for i in range(len(L_f)):

print(powers[i])

# シミュレーション結果のグラフ描画

v_log_x = np.log(L_Ns)

L_log_x = ['{}'.format(v) for v in L_Ns]

for j in range(len(L_f)):

plt.plot(v_log_x, powers[j], label = 'f = {}'.format(L_f[j]))