Up

効果量と最適な有意水準

対立仮説で表されている条件差は効果量として標準化されるが、与えられた効果量に対応する最適な有意水準は期待効用の観点から決めることができる。独立な２つのサンプルの平均値の検定、および被験者間１要因分散分析における場合について計算を試みた。まず、期待効用について説明する。その後、独立な２つのサンプルの平均値の差の効果量、および被験者間１要因分散分析の効果量について具体的に取り上げて、最適な有意水準の検討を行う。最適な有意水準は、効果量に応じて変わる。

期待効用

記号を以下のように導入する。

: 帰無仮説が正しい確率

第一種の誤りの確率、すなわち有意水準

帰無仮説が正しいときに、効果有りと判断することの効用

帰無仮説が正しいときに、効果無しと判断することの効用

: 対立仮説が正しい確率

第二種の誤りの確率、すなわち: 検定力

対立仮説が正しいときに、効果有りと判断することの効用

対立仮説が正しいときに、効果無しと判断することの効用

検定統計量の確率分布は、図１のように模擬的に表すことができる。

図１

このとき、期待効用は有意水準の関数として、次式で表すことができる。

式（１）で与えられる期待効用を、以下の効果量に対して求めてみる。

独立な２つのサンプルの平均値の場合

独立な２つのサンプルの平均値の差の効果量は、次式で与えられる。

ここで、とは２つのサンプルの母集団の平均値を表し、は標準偏差を表す。

独立な２つのサンプルのサイズがとであるとき、帰無仮説「平均値が等しい」のもとでは、ｔ値は自由度のｔ分布に従い、対立仮説の下では、非心度

の非心ｔ分布に従う。

上記のモデルに対するスクリプトをリスト１に示す。

リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。

= RESTART: C:\*****\two_meansfiles\two_sample_means.py

Input data file (*.csv) = sampled08.xlsx

入力データファイルsampled08.xlsxは、図A1に示すように用意されている。

図A1

図A1のデータは、表A1に示す形式で用意されている。

表A1
	パラメータ値	確率	有意な結果の効用	非有意な結果の効用
帰無仮説		帰無仮説が正しい確率	帰無仮説が正しいときの有意な結果の効用	帰無仮説が正しいときの非有意な結果の効用
対立仮説		対立仮設が正しい確率	対立仮説が正しいときの有意な結果の効用	対立仮説が正しいときの非有意な結果の効用
効果量	効果量の値
M	個数
N	個数

図Ａ1のデータに対してリスト１のスクリプトを実行すると、図Ａ2の結果を得る。

図A2

左側のグラフにおいて青の実線が期待効用のグラフである。グラフの最大値の位置（最適な有意水準）が赤丸で示されている。最適な有意水準（opt-α）として2.4％が得られている。

図A1のパラメータ値における効果量0.8の値を0.5に下げたものを図A3のように用意した。

図A3

図A3のファイルを入力データファイルとしてリスト１のスクリプトを実行すると、図A4の結果を得た。

図A4

最適な有意水準がopt-α=10.8%になっている。

有意水準が大きくなると、差がないのに差があると判断する（False Alarm: FA）確率が高くなる。FAの確率を下げるために、FAの効用をからにさらに負の値として絶対値の大きいものを設定してみた（図A5）。

図A5

図A5のデータに対してリスト１のスクリプトを実行すると、図A6の結果を得る。

図A6

最適な有意水準が2.8％になっている。

なお、用いたファイルは圧縮ファイルtwo_meansfiles.zipにまとめた。

被験者間１要因分散分析の場合

被験者間1要因分散分析における効果量は、次式で与えられる。

ここで、は水準間の効果の分散であり、は水準内の分散を表す。

要因の数を、各水準内のデータ数をとするとき、効果量がである対立仮説の下では、検定統計量は非心度の非心F分布に従う。ここで、

である。自由度は、、である。

帰無仮説の下では、自由度は対立仮説の非心F分布の自由度と同じであるが、非心度は０のF分布に従う。

上のモデルに基づいて期待効用を算出するスクリプトとしてリスト２に示すものを作成した。

リスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。

= RESTART: C:\***** \anovafiles\anova_effect_size.py

Input data file (*.csv) = sample.xlsx

入力データファイルは、Excelファイルであるが、上の例では図B1に示すファイルが設定されている。

図B1

図B1のExcelファイルは、表B1の形式に従って用意されている。

表B1
	パラメータ値	確率	有意な結果の効用	非有意な結果の効用
帰無仮説		帰無仮説が正しい確率	帰無仮説が正しいときの有意な結果の効用	帰無仮説が正しいときの非有意な結果の効用
対立仮説		対立仮設が正しい確率	対立仮説が正しいときの有意な結果の効用	対立仮説が正しいときの非有意な結果の効用
効果量	効果量の値
a	水準数
N	水準内の個体数

入力データファイル名を設定してEnterキーを押すと、計算が始まり、図B1のファイルの場合は、図B2のグラフが表示される。

図B2

左のグラフの青の実線が期待効用のグラフである。が最適な有意水準であることが示されている。

効果量をに小さくしたファイルを図B3のように用意した。

図B3

図B3のファイルを入力データファイルとしてリスト２のスクリプトを実行すると、図B4の結果を得る。

図B4

最適な有意水準としてが示されている。有意水準が通常用いられている値より大きいので、帰無仮説が正しいときに有意な結果を得る（False Alarm: FA）ことに対する負の効用をより強調した図B5のデータを用意した。

図B5

図B5の入力データに対してリスト２のスクリプトを実行すると、図B6の結果を得る。

最適な有意水準がになっている。

なお、用いたファイルは圧縮ファイルanovafiles.zipにまとめた。

参考文献

Hogg, R. V., McKean, J. W., and Craig, A. T. (2005) Introduction to mathematical statistics, 6th ed. Pearson Prentice Hall.

Johnson, N. L., Kotz, S., and Balakrishnan, N. (1995) Continuous univariate distributions, Vol.2, 2nd ed. John wiley & Sons, Inc.

Kirk, R. E. (1995) Experimental design: Procedures for the behavioral sciences, 3rd ed. Brooks/Cole Publishing Company

Myers, J. L., and Well, A. D. (2003) Research design and statistical analysis, 2nd ed. Lawrence Erlbaum Associates, Publishers

岡本安晴（2009）データ分析のための統計学入門――統計学の考え方――．おうふう（第５章　効果量と検定力）

芝祐順祐・渡部洋・石塚智一（1984）　統計用語辞典、新曜社

Winer, B. J., Brown, D. R. and Michels, K. M. (1991) Statistical principles in experimental design, 3rd ed. McGraw-Hill, Inc.

リスト１　独立な２つのサンプルの場合

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats as ss

import scipy.special as sspec

class ExpectedUtil():

"""

期待効用値算出クラス、t分布対応

"""

def __init__(self, p_H0, p_H1, U0sig, U0nsig, U1sig, U1nsig,

df, nc):

self.p_H0 = p_H0

self.p_H1 = p_H1

self.U0sig = U0sig

self.U0nsig = U0nsig

self.U1sig = U1sig

self.U1nsig = U1nsig

self.df = df

self.nc = nc

def calc_criterion(self, a):

c = sspec.stdtrit(self.df, 1.0 - a)

return c

def expectedU(self, p):

a = p

c = self.calc_criterion(a)

b = sspec.nctdtr(self.df, self.nc, c)

exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \

self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)

return exp_util

def exp_uh0(self, p):

"""

帰無仮説の下での期待効用

"""

a = p

exp_util = a * self.U0sig + (1 - a) * self.U0nsig

return exp_util

def exp_uh1(self, p):

"""

対立仮設の下での期待効用

"""

a = p

c = self.calc_criterion(a)

b = sspec.nctdtr(self.df, self.nc, c)

exp_util = (1 - b) * self.U1sig + b * self.U1nsig

return exp_util

def calc_beta(self, a):

c = sspec.stdtrit(self.df, 1.0 - a)

return sspec.nctdtr(self.df, self.nc, c)

pd_df = pd.read_excel(input('Input data file (*.csv) = '))

print(pd_df)

vals = pd_df.values

print(vals)

e_s = vals[2][1] # standardized effect size delta

M = vals[3][1]

N = vals[4][1]

df = M + N - 2

nc = e_s * (((M * N) / (M + N)) ** 0.5)

p_H0 = vals[0][2] / (vals[0][2] + vals[1][2])

p_H1 = vals[1][2] / (vals[0][2] + vals[1][2])

U0sig = vals[0][3]

U0nsig = vals[0][4]

U1sig = vals[1][3]

U1nsig = vals[1][4]

print('e_s =', e_s)

print('M =', M)

print('N =', N)

print('df =', df)

print('nc =', nc)

print('p_H0 =', p_H0)

print('p_H1 =', p_H1)

print('U0sig =', U0sig)

print('U0nsig =', U0nsig)

print('U1sig =', U1sig)

print('U1nsig =', U1nsig)

plt.figure(figsize=(12,5))

plt.subplot(1,3,(1,2)) # 左２区画に期待効用曲線を描画

e_u = ExpectedUtil(p_H0, p_H1, U0sig, U0nsig, U1sig, U1nsig,

df, nc)

ps = np.linspace(1.0e-9, 1.0 - 1.0e-9, 10000)

expected_utils = e_u.expectedU(ps)

uh0 = e_u.exp_uh0(ps)

uh1 = e_u.exp_uh1(ps)

ck_idx = np.argsort(-expected_utils) # 最大値を与えるインデックス

opt_a = ps[ck_idx[0]] # 最大期待効用を与える有意水準

print('opt_a =\n', opt_a)

plt.plot(ps, expected_utils, lw=3, label=f'Exp.U')

plt.plot(ps, uh0, lw=3, ls='--', label='H0')

plt.plot(ps, uh1, lw=3, ls='-.', label='H1')

plt.plot([opt_a], [expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,

label=fr'opt-$\alpha$={opt_a:.3f}')

plt.xlabel(r'$\alpha$', fontsize= 18)

plt.ylabel('Expect.Utility', fontsize=18)

plt.legend(handlelength=5)

plt.title(fr'$\delta$={e_s}, nc={nc:.3f}, df={df}' + '\n' +

fr'U0sig={U0sig}, U0nsig={U0nsig}, p0={p_H0:.3f}' + '\n' +

fr'U1sig={U1sig}, U1nsig={U1nsig}, p1={p_H1:.3f}', fontsize=16)

plt.tight_layout()

plt.subplot(1, 3, 3) # αに対するβと検出力のグラフ

alphas = np.linspace(1.0e-9, 1.0 - 1.0e-9, 10_000)

betas = e_u.calc_beta(alphas)

pwrs = 1.0 - betas

plt.plot(alphas, betas, label=r'$\beta$')

plt.plot(alphas, pwrs, label='power')

plt.plot(opt_a, e_u.calc_beta(opt_a), 'ro', markersize=12,

label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}' + '\n' +

fr'$\alpha$={opt_a:.3f}')

plt.xlabel(r'$\alpha$', fontsize=18)

plt.ylabel(r'$\beta$, power', fontsize=14)

plt.legend()

plt.tight_layout()

plt.show()

リスト２　被験者間１要因分散分析の場合

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats as ss

import scipy.special as sspec

class ExpectedUtil():

"""

期待効用値算出クラス、F分布対応

"""

def __init__(self, p_H0, p_H1, U0sig, U0nsig, U1sig, U1nsig,

df1, df2, nc):

self.p_H0 = p_H0

self.p_H1 = p_H1

self.U0sig = U0sig

self.U0nsig = U0nsig

self.U1sig = U1sig

self.U1nsig = U1nsig

self.df1 = df1

self.df2 = df2

self.nc = nc

def calc_criterion(self, a):

c = sspec.fdtri(self.df1, self.df2, 1.0 - a)

return c

def expectedU(self, p):

a = p

c = self.calc_criterion(a)

b = sspec.ncfdtr(self.df1, self.df2, self.nc, c)

exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \

self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)

return exp_util

def e_uH0(self, p):

a = p

#c = self.calc_criterion(a)

#b = sspec.ncfdtr(self.df1, self.df2, self.nc, c)

exp_util = a * self.U0sig + (1 - a) * self.U0nsig

return exp_util

def e_uH1(self, p):

a = p

c = self.calc_criterion(a)

b = sspec.ncfdtr(self.df1, self.df2, self.nc, c)

exp_util = (1 - b) * self.U1sig + b * self.U1nsig

return exp_util

def calc_beta(self, a):

c = sspec.fdtri(self.df1, self.df2, 1.0 - a)

return sspec.ncfdtr(self.df1, self.df2, self.nc, c)

pd_df = pd.read_excel(input('Input data file (*.csv) = '))

print(pd_df)

vals = pd_df.values

print(vals)

e_s = vals[2][1] # f: effect size

a = vals[3][1]

N = vals[4][1]

df1 = a - 1

df2 = a * (N - 1)

nc = a * N * (e_s**2)

p_H0 = vals[0][2] / (vals[0][2] + vals[1][2])

p_H1 = vals[1][2] / (vals[0][2] + vals[1][2])

U0sig = vals[0][3]

U0nsig = vals[0][4]

U1sig = vals[1][3]

U1nsig = vals[1][4]

print('e_s =', e_s)

print('a =', a)

print('N =', N)

print('nc =', nc)

print('p_H0 =', p_H0)

print('p_H1 =', p_H1)

print('U0sig =', U0sig)

print('U0nsig =', U0nsig)

print('U1sig =', U1sig)

print('U1nsig =', U1nsig)

plt.figure(figsize=(12,5))

plt.subplot(1,3,(1,2)) # 左２区画に期待効用曲線を描画

e_u = ExpectedUtil(p_H0, p_H1, U0sig, U0nsig, U1sig, U1nsig,

df1, df2, nc)

ps = np.linspace(1.0e-9, 1.0 - 1.0e-9, 10_000)

expected_utils = e_u.expectedU(ps)

e_uh0 = e_u.e_uH0(ps)

e_uh1 = e_u.e_uH1(ps)

ck_idx = np.argsort(-expected_utils) # 最大値を与えるインデックス

opt_a = ps[ck_idx[0]] # 最大期待効用を与える有意水準

print('opt_a =\n', opt_a)

plt.plot(ps, expected_utils, label=f'Exp.U')

plt.plot(ps, e_uh0, ls='--', label='e_uh0')

plt.plot(ps, e_uh1, ls='-.', label='e_uh1')

plt.plot([opt_a], [expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,

label=fr'opt-$\alpha$={opt_a:.3f}')

plt.xlabel(r'$\alpha$', fontsize= 18)

plt.ylabel('Expect.Utility', fontsize=18)

plt.legend(handlelength=5)

plt.title(fr'$f$={e_s}, nc={nc:.3f}, df1={df1}, df2={df2}' + '\n' +

fr'U0sig={U0sig}, U0nsig={U0nsig}, p0={p_H0:.3f}' + '\n' +

fr'U1sig={U1sig}, U1nsig={U1nsig}, p1={p_H1:.3f}', fontsize=16)

plt.tight_layout()

plt.subplot(1, 3, 3) # αに対するβと検出力のグラフ

alphas = np.linspace(1.0e-9, 1.0 - 1.0e-9, 10_000)

betas = e_u.calc_beta(alphas)

pwrs = 1.0 - betas

plt.plot(alphas, betas, label=r'$\beta$')

plt.plot(alphas, pwrs, label='power')

plt.plot(opt_a, e_u.calc_beta(opt_a), 'ro', markersize=12,

label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}' + '\n' +

fr'$\alpha$={opt_a:.3f}')

plt.xlabel(r'$\alpha$', fontsize=18)

plt.ylabel(r'$\beta$, power', fontsize=14)

plt.legend()

plt.tight_layout()

plt.show()

Home