python statsmodel の power module で、検定力 と、Sample Size の計算ができます。
ドキュメントを見る限り、R の pwr と同じようなことができるようなので、実際に動かして Sample Size の計算をしてみました。
実施した内容を記載します。

pwr と、statsmodel のモジュールの対応表

モジュールの対応関係は以下の通りです。
tt_ind_solve_power、zt_ind_solve_power は、それぞれ、TTestIndPowerと、GofChisquarePower のショートカットです。
基本的にできることは、pwr のほうが多く、他にstatsmodel には、片側t検定のためのTTestPower、tt_solve_power があります。
半分近くは、NormalIndPower を駆使して、pwr と同様の値を計算できるようで、試している記事がありました。
ただ、計算精度の問題でRと結果が異なってきたりする場合があるようなので、statsmodel が想定している使い方の範囲で使うのがよいかと思います。

説明 pwrのモジュール名 statsmodelのモジュール名
平均値に関するt検定(1群,2群,対応あり) pwr.t.test statsmodels.stats.power.TTestIndPower,statsmodels.stats.power.tt_ind_solve_power
サンプルサイズの異なる2群の平均値に関するt検定 pwr.t2n.test statsmodels.stats.power.TTestIndPower
正規分布の平均値の検定 pwr.norm.test -
比率の検定(1標本) pwr.p.test statsmodels.stats.power.NormalIndPower
２つの比率の差の検定(サンプル数が等しい場合) pwr.2p.test statsmodels.stats.power.NormalIndPower
２つの比率の差の検定(サンプル数が異なる場合) pwr.2p2n.test statsmodels.stats.power.NormalIndPower
カイ二乗検定 pwr.chisq.test statsmodels.stats.power.GofChisquarePower,statsmodels.stats.power.zt_ind_solve_power
相関係数の検定 pwr.r.test statsmodels.stats.power.NormalIndPower
一元配置分散分析 pwr.anova.test statsmodels.stats.power.FTestAnovaPower
一般線形モデル(分散分析・回帰分析) pwr.f2.test statsmodels.stats.power.FTestPower

参考

各モジュールの実行結果

Sample Size と、検定力 を計算して出力します。

平均値に関するt検定(1群,2群,対応あり) statsmodels.stats.power.tt_ind_solve_power

  • Sample Size を計算 
from statsmodels.stats.power import tt_ind_solve_power
# 平均値の差
delta = 0.5
# 標準偏差
sd    = 1
# 有意水準
sig   = 0.05
# 検定力
power = 0.8
# 効果量 (平均値の差 / 標準偏差 )
effect_size = delta / sd
# 両側検定
alternative='two-sided'
n = tt_ind_solve_power(effect_size=effect_size,  alpha=sig, power=power, alternative=alternative)
print(n)

63.765611775409525
  • 検定力 を計算
    TTestIndPower().power() を使用して、検定力を計算します。 
# 平均値の差
delta = 0.5
# 標準偏差
sd    = 1
# 有意水準
sig   = 0.05
# サンプルサイズ   
n = 10
# 効果量 あまりよくわかっていない。 平均値の差 / 標準偏差     
effect_size = delta / sd
# 両側検定
alternative='two-sided'
import statsmodels.stats.power as smp
n = smp.TTestIndPower().power(effect_size, n, sig,  alternative=alternative)
print(n)

0.185095655379

サンプルサイズの異なる2群の平均値に関するt検定 statsmodels.stats.power.TTestIndPower

tt_ind_solve_power、TTestIndPower().solve_power() を使います。
サンプルサイズが異なることは、パラメータ ratio を使って表現します。

  • Sample Size を計算 
from statsmodels.stats.power import tt_ind_solve_power
# 有意水準
sig   = 0.05
# 検定力
power = 0.8
# 効果量
effect_size = 0.6
# レシオ  サンプルサイズは、2:3
ratio = 60./90
# 片側検定
alternative ='larger'
n = tt_ind_solve_power(effect_size=effect_size,  alpha=sig, power=power, alternative=alternative, ratio=ratio)
print("alternative=larger: ", n)

# 両側検定
alternative='two-sided'
n = tt_ind_solve_power(effect_size=effect_size,  alpha=sig, power=power, alternative=alternative, ratio=ratio)
print("alternative=two-sided: ", n)

alternative=larger:  43.76924760999283
alternative=two-sided:  55.68220224349466
  • 検定力 を計算
    TTestIndPower().power を使用して検定力を計算します。 
# サンプルサイズ1
n1 = 90
# サンプルサイズ2
n2 = 60
# レシオ      
ratio = n2/n1
# 効果量      
effect_size = 0.6
# 両側検定か、片側検定か、'two-sided' (default), 'larger', 'smaller'　が指定可能    
alternative ='larger'
# 有意水準
sig   = 0.05
import statsmodels.stats.power as smp
p= smp.TTestIndPower().power(effect_size=effect_size, nobs1=n1, ratio=ratio, alpha=sig, alternative=alternative)
# alternative 'larger' 指定だと、片側検定になる
print("alternative=larger:", p)
p= smp.TTestIndPower().power(effect_size=effect_size, nobs1=n1, ratio=ratio, alpha=sig)
# 指定なしは、two-sided 両側検定
print("alternative=two-sided:", p)

alternative=larger: 0.973726135488
alternative=two-sided: 0.947015360249

比率の検定(1標本) statsmodels.stats.power.NormalIndPower

  • Sample Size を計算 
 from statsmodels.stats.proportion import proportion_effectsize as proportion_effectsize
proportion_effectsize(0.5, 0.4)
import statsmodels.stats.power as smp
p=0.8
smp.NormalIndPower().solve_power(effect_size=0.2013579207903309, power=p, alpha=0.05, ratio=0)

193.58387311494354
  • 検定力 を計算
    proportion_effectsize で、比率を計算し、計算結果を、NormalIndPower().power のeffect_sizeとして設定します。 
 from statsmodels.stats.proportion import proportion_effectsize as proportion_effectsize
proportion_effectsize(0.5, 0.4)
import statsmodels.stats.power as smp
smp.NormalIndPower().power(effect_size=0.2013579207903309, nobs1=60, alpha=0.05, ratio=0)

0.34470140912721525

２つの比率の差の検定(サンプル数が等しい場合) statsmodels.stats.power.NormalIndPower

  • Sample Size を計算 
import statsmodels.stats.power as smp
# 効果量
effect_size=0.3   
# 検出力
power = 0.8
# 有意水準
sig = 0.05
p=smp.NormalIndPower().solve_power(effect_size=effect_size, power=power, alpha=sig, alternative='two-sided')
print(p)

174.41912242947043
  • 検定力を計算
import statsmodels.stats.power as smp
# 効果量
effect_size=0.3   
# サンプルサイズ
n = 80
# 有意水準
sig = 0.05
p=smp.NormalIndPower().power(effect_size=effect_size, nobs1=n, alpha=sig, alternative='two-sided')
print(p)

0.475100870573

２つの比率の差の検定(サンプル数が異なる場合) statsmodels.stats.power.NormalIndPower

サンプリング時の最適なサンプルサイズをRパッケージ{pwr}で求める - 六本木で働くデータサイエンティストのブログ の母不良率（○○率）を求めたい場合 と同様の値が計算可能か試してみます。

  • Sample Size を計算
import statsmodels.stats.power as smp
# 有意水準
sig   = 0.05
# 検定力
power = 0.9
# サンプルサイズ   
n = 2.5e8
# 効果量
effect_size = 0.2
# 両側検定
alternative='two-sided'
result = smp.NormalIndPower().solve_power(effect_size, nobs1=n, alpha=sig, power=power, ratio=None, alternative=alternative)
print(result * n)

1706.6078223432073

上記、一致しません。おそらく計算精度の問題で、うまく算出できていないようです。
サンプルサイズをn = 2.5e4あたりにしたところ、Rの、pwr.2p2n.test の結果と一致するようになりました。

  • 検定力を計算 
import statsmodels.stats.power as smp
# 有意水準
sig   = 0.05
# サンプルサイズ   
n = 2.5e8
# 効果量
effect_size = 0.2
# 両側検定
alternative='two-sided'
result = smp.NormalIndPower().power(effect_size, n, sig, ratio=262.6855/n, alternative=alternative)
print(result)

0.899999717107

検定力については、四捨五入するとRの計算結果と、ほぼ一致します。

カイ二乗検定 statsmodels.stats.power.GofChisquarePower

カイ二乗検定の、サンプル数、検定力の計算には、statsmodels.stats.power.GofChisquarePowerを使用します。

  • Sample Size を計算 
import statsmodels.stats.power as smp
# サンプル数の計算、n_binsは自由度、 + 1する
smp.GofChisquarePower().solve_power(0.1, n_bins=(5-1)*(6-1) + 1, alpha=0.05, power=0.8)

2096.0784652617726
  • 検定力を計算 
import statsmodels.stats.power as smp
# 検定力の計算、n_bins は自由度に + 1 する
smp.GofChisquarePower().power(0.346, nobs=140, n_bins=(2-1)*(3-1) + 1, alpha=0.01)

0.88540533954389722

相関係数の検定 statsmodels.stats.power.NormalIndPower

相関係数の検定には、statsmodels.stats.power.NormalIndPower　と、numpy.arctanh(x) を使います。
numpy.arctanh(x) を使うのは母相関の検定と推定 あたりの内容が理由かと思われます。
参考記事に記載がありますが、Rの計算結果とは完全に一致せず、だいたい一致する結果が得られます。

  • Sample Size を計算 
import statsmodels.stats.power as smp
import numpy as np
# サンプル数の調整のため、何故か3を足す。理由はわからず。
smp.NormalIndPower().solve_power(np.arctanh(0.3), alpha=0.05, ratio=0, power= 0.8, alternative='two-sided') + 3

84.92761044616469
  • 検定力を計算 
import statsmodels.stats.power as smp
import numpy as np
# サンプル数の調整のため、3を引く。理由はわからず。
smp.NormalIndPower().power(np.arctanh(0.3), nobs1=50-3, alpha=0.05, ratio=0, alternative='two-sided')

0.56436763896354003

一元配置分散分析 statsmodels.stats.power.FTestAnovaPower

FTestAnovaPower を使って、計算できます。

  • Sample Size を計算 
import statsmodels.stats.power as smp
# 効果量
effect_size=0.28
# 効果量
alpha=0.05
# グループの数
k_groups=4
#  検出力
power=0.8
smp.FTestAnovaPower().solve_power(effect_size=effect_size, alpha=alpha, power=power, k_groups=k_groups)

143.03088865487553
  • 検定力を計算 
import statsmodels.stats.power as smp
# 効果量
effect_size=0.28
# サンプル数
nobs=80
# 有意水準
alpha=0.05
# グループの数
k_groups=4
smp.FTestAnovaPower().power(effect_size, nobs, alpha,  k_groups=k_groups)

0.51498349798649112

一般線形モデル(分散分析・回帰分析) statsmodels.stats.power.FTestPower

FTestPower を使って計算できます。

  • Sample Size の計算 
import statsmodels.stats.power as smp
#  効果量
effect_size=np.sqrt(0.1/(1-0.1))
# 有意水準
alpha=0.05
# 検出力
power = 0.8
# 分母の自由度
df_denom = 5
smp.FTestPower().solve_power(effect_size=effect_size, alpha=alpha, power=power, df_denom=df_denom)

115.10364322775642
  • 検定力を計算 
import statsmodels.stats.power as smp
#  効果量
effect_size=np.sqrt(0.1/(1-0.1))
# 分母の自由度
df_denom = 5
# 分子の自由度
df_num = 89
# 有意水準
alpha=0.05
smp.FTestPower().power(effect_size=effect_size, df_num=df_num, df_denom=df_denom, alpha=alpha)

0.67358904832076494

感想

statsmodel で、Sample Size の計算を実施しました。
以下、感想をまとめます。

  • pwrと同レベルのことは実施できないが、NormalIndPower を駆使するとそれなりにできる。
  • NormalIndPower を駆使する場合、計算誤差が発生する場合があるので、注意。
  • joepy: Statistical Power in Statsmodels が非常に参考になった。
  • 実業務で使用する際は、TTestIndPower、GofChisquarePower、FTestAnovaPower、FTestPower を使うと思う。

以上です。

コメント

カテゴリー

最近の投稿