「数値シミュレーションで読み解く統計のしくみ」をPythonで写経 ～ 第６章6.4「いろいろな検定におけるサンプルサイズ設計の実践」
第６章「適切な検定のためのサンプルサイズ設計」
書籍の著者 小杉考司 先生、紀ノ定保礼 先生、清水裕士 先生
この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第６章「適切な検定のためのサンプルサイズ設計」の6.4節「いろいろな検定におけるサンプルサイズ設計の実践」の Python写経活動 を取り扱います。
第６章はシミュレーションを通じて、調査・実験の事前段階でサンプルサイズを設計する方法を学びます。
今回は３つの検定のサンプルサイズ設計を学びます。
R の基本的なプログラム記法はざっくり Python の記法と近いです。
コードの文字の細部をなぞって、R と Python の両方に接近してみましょう！
ではテキストを開いてシミュレーションの旅に出発です🚀
はじめに
テキスト「数値シミュレーションで読み解く統計のしくみ」のご紹介
このシリーズは書籍「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」（技術評論社、「テキスト」と呼びます）の Python写経です。
テキストは、2023年9月に発売され、副題「Rでためしてわかる心理統計」のとおり、統計処理に定評のある R の具体的なコードを用いて、心理統計の理解に役立つ数値シミュレーションを実践する素晴らしい書籍です。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」初版第１刷、著者 小杉考司・紀ノ定保礼・清水裕士、技術評論社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
6.4 いろいろな検定におけるサンプルサイズ設計の実践
この節では以下の３つのサンプルサイズ設計を実践します。
・相関係数
・一元配置分散分析
・反復測定分散分析
Pythonで必要サンプルサイズを計算する際には、ぜひ、pingouin ライブラリをご利用ください！
記事では「アディショナルタイム」で pingouin を用いたサンプルサイズ設計のコードを掲載しています。
この記事は Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
概ね確率分布の特性値算出には scipy.stats を用い、乱数生成には numpy.random.generator を用いるようにしています。
主に利用するライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
# 確率・統計
import scipy.stats as stats
from statsmodels.stats import power # サンプルサイズ
import pingouin as pg # 統計便利ツール
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
6.4.1 相関係数のサンプルサイズ設計
相関係数の検定統計量である$${t}$$値は、帰無仮説の相関係数と母相関係数のズレを反映した、非心ｔ分布に従います。
テキストから、母相関係数$${\rho}$$における非心度の数式をお借りします。
$$
\lambda = \cfrac{\rho}{\sqrt{1 - \rho^2}} \sqrt{n}
$$
テキストに合わせて、相関係数のタイプⅡエラー確率計算を関数化します。
上記の非心度$${\lambda}$$の数式は、関数内の「非心度の計算」で使用します。
scipy.stats の nct で非心ｔ分布を扱います。
### 245ページ 相関係数のタイプⅡエラー確率を算出する関数
def t2e_cor(alpha, rho, n):
# 自由度の計算
df = n - 2
# 非心度の計算
lambda_ = rho / np.sqrt(1 - rho**2) * np.sqrt(n)
# 臨界値の計算
cv = stats.t.ppf(q=1 - alpha / 2, df=df)
# タイプⅡエラー確率の計算
type2_error = stats.nct.cdf(x=cv, df=df, nc=lambda_)
return type2_error
■ 検出したい相関係数の大きさを指定してサンプルサイズをシミュレーション
検出したい相関係数を$${0.3}$$、有意水準$${\alpha=5\%}$$、タイプⅡエラー確率の閾値$${\beta=20\%}$$の条件で、必要サンプルサイズを計算します。
シミュレーション回数の上限を1000回としています。
### 245ページ 相関係数ののサンプルサイズ設計
# 検出したい母相関係数の大きさを0.3に見積もった場合、有意水準α=0.05,
# タイプⅡエラー確率β=0.20（検出力=0.80）の下で必要サンプルサイズは82
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
rho = 0.3 # 検出したい母相関係数の大きさ
max_iter = 1000 # シミュレーション回数の最大値
## シミュレーション
# タイプⅡエラー確率が閾値β以下になるまでサンプルサイズを増やす（上限1000）
for n in range(5, max_iter+1):
# サンプルサイズnのときのタイプⅡエラー確率を算出
type2_error = t2e_cor(alpha, rho, n)
# タイプⅡエラー確率が閾値β以下になったらシミュレーションを終了
if type2_error <= beta:
break
## 結果
# 必要サンプルサイズ
print(f'必要サンプルサイズ: {n}')
# そのときのタイプⅡエラー確率
print(f'タイプⅡエラー確率: {type2_error:.7f}')【実行結果】
テキストと同じ値になりました。
テキストの言葉をお借りしますと「検出したい相関係数の大きさを 0.3 に見積もる場合、有意水準$${\alpha=0.05}$$、タイプⅡエラー確率の閾値$${\beta=0.20}$$の下で、必要サンプルサイズは$${82}$$となりました」、です。
アディショナルタイム１
テキストから少々離れて、ここからはアディショナルタイムに突入です！
Pythonのライブラリを用いて相関係数のサンプルサイズ計算を行います。
■ アディショナルタイム①「pingouinのpower_corr」
pingouin ライブラリを用いて必要サンプルサイズを計算します。
### 別解 pingouinのpower_corrで必要サンプルサイズを算出
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
rho = 0.3 # 検出したい母相関係数の大きさ
## 必要サンプルサイズの算出 3を調整
n = np.ceil(pg.power_corr(r=rho, power=1 - beta, alpha=alpha,
alternative='two-sided') - 3)
n【実行結果】
テキストの計算結果と同じになりました！
サンプルサイズ$${n=82}$$のときのタイプⅡエラー確率を計算します。
### 別解 続き pingouinのpower_corrで必要サンプルサイズのタイプⅡエラー確率を算出
1 - pg.power_corr(r=rho, n=n+3, alpha=alpha, alternative='two-sided')【実行結果】
テキストの計算結果とほぼほぼ同じになりました！
■ アディショナルタイム②「statsmodelsのNormalIndPower」
statsmodels ライブラリを用いて必要サンプルサイズを計算します。
### 別解 statsmodelsのNormalIndPowerで必要サンプルサイズを算出
# 参考
# https://www.monotalk.xyz/blog/With-python-statsmodel-calculate-the-verification-power-and-the-Sample-Size/
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
rho = 0.3 # 検出したい母相関係数の大きさ
## 必要サンプルサイズの算出
n = np.ceil(
power.NormalIndPower().solve_power(
effect_size=np.arctanh(rho), alpha=alpha, ratio=0, power=1-beta,
alternative='two-sided'))
n【実行結果】
テキストの計算結果と同じになりました！
サンプルサイズ$${n=82}$$のときのタイプⅡエラー確率を計算します。
### 別解 statsmodelsのNormalIndPowerで必要サンプルサイズのタイプⅡエラー確率を算出
1 - power.NormalIndPower().power(effect_size=np.arctanh(rho), nobs1=n,
alpha=alpha, ratio=0, alternative='two-sided')【実行結果】
テキストの計算結果と概ね同じになりました！（若干、値は大きいかもです）
6.4.2 一元配置分散分析のサンプルサイズ設計
テキストによると、分散分析は帰無仮説が正しくないとき、検定統計量は非心$${F}$$分布に従います。
テキストの$${F}$$分布と非心$${F}$$分布のチャートを描画します。
### 246ページ 図6.8 F分布（破線）と非心F分布（実線）
## 設定と準備
f_sq = 0.1 # 効果量
k = 4 # 群の数
ng = 30 # 群ごとのサンプルサイズ
lambda_ = f_sq * ng * k # 非心度
## 描画処理
# 描画領域の設定
plt.figure(figsize=(6, 3))
# x軸の値
line_x = np.linspace(0, 10, 1001)
# F分布の確率密度関数の描画 scipy.statsのfを使用
plt.plot(line_x, stats.f.pdf(line_x, dfn=k - 1, dfd=(ng - 1)*k),
color='tab:orange', ls='--', label='帰無分布')
# 非心F分布の確率密度関数の描画 scipy.statsのncfを使用
plt.plot(line_x, stats.ncf.pdf(line_x, dfn=k - 1, dfd=(ng - 1)*k, nc=lambda_),
color='tab:blue', label='非心F分布')
# x軸ラベル、y軸ラベル、タイトル、凡例の設定
plt.xlabel('x')
plt.ylabel('Density')
plt.title('F分布と非心F分布')
plt.legend();【実行結果】
非心$${F}$$分布の非心度は、分散分析の母効果量$${\eta^2}$$（相関比の二乗）を用いて計算できるそうです。
$$
\lambda = \cfrac{\eta^2}{1 - \eta^2} \ n
$$
テキストに合わせて、一元配置分散分析のタイプⅡエラー確率計算を関数化します。
上記の非心度$${\lambda}$$の数式は、関数内の「非心度の計算」で使用します。
サンプルサイズ$${n}$$は、群の数$${k}$$と群ごとのサンプルサイズ$${n_g}$$を用いて、$${n_g k}$$で求められます。
scipy.stats の ncf で非心$${F}$$分布を扱います。
### 247ページ 一元配置分散分析のタイプⅡエラー確率を算出する関数
# alpha:有意水準, eta_sq:母効果量（相関比の二乗）, k:群の数, ng:群ごとのサンプルサイズ
def t2e_1way(alpha, eta_sq, k, ng):
# 自由度の計算
df1 = k - 1 # 群間の自由度
df2 = (ng - 1) * k # 群内の自由度
# 非心度の計算
lambda_ = eta_sq / (1 - eta_sq) * ng * k
# 臨界値の計算
cv = stats.f.ppf(q=1 - alpha, dfn=df1, dfd=df2)
# タイプⅡエラー確率の計算
type2_error = stats.ncf.cdf(x=cv, dfn=df1, dfd=df2, nc=lambda_)
return type2_error
■ 一元配置分散分析のサンプルサイズをシミュレーション
3群で検出したい効果量を$${0.06}$$、有意水準$${\alpha=5\%}$$、タイプⅡエラー確率の閾値$${\beta=20\%}$$の条件で、必要サンプルサイズを計算します。
シミュレーション回数の上限を10000回としています。
### 247～248ページ 一元配置分散分析のサンプルサイズ設計
# ３群で検出したい効果量をη²=0.06で設定したとき、必要サンプルサイズは
# 有意水準α=0.05, タイプⅡエラー確率β=0.20（検出力=0.80）の下で各群52, 全体で156
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
k = 3 # 群の数
eta_sq = 0.06 # 検出したいの効果量（相関比の二乗）の大きさ
iter = 10000 # シミュレーション回数の最大値
## シミュレーション
# タイプⅡエラー確率が閾値β以下になるまでサンプルサイズを増やす（上限10000）
for ng in range(5, iter + 1):
# サンプルサイズnのときのタイプⅡエラー確率を算出
type2_error = t2e_1way(alpha, eta_sq, k, ng)
# タイプⅡエラー確率が閾値β以下になったらシミュレーションを終了
if type2_error <= beta:
break
## 結果
# 必要サンプルサイズ（各群）
print(f'各群の必要サンプルサイズ: {ng}')
# 必要サンプルサイズ（全体）
print(f'全体の必要サンプルサイズ: {ng * k}')
# そのときのタイプⅡエラー確率
print(f'タイプⅡエラー確率: {type2_error:.7f}')【実行結果】
テキストと同じ値になりました。
テキストの言葉をお借りしますと「３群で検出したい効果量を$${0.06}$$で設定したとき、有意水準$${\alpha=0.05}$$、タイプⅡエラー確率の閾値$${\beta=0.20}$$の下で、必要サンプルサイズは各群$${52}$$人、全体で$${156}$$でした」、です。
アディショナルタイム２
テキストから少々離れて、ここからはアディショナルタイムに突入です！
Pythonのライブラリを用いて一元配置分散分析のサンプルサイズ計算を行います。
■ アディショナルタイム①「pingouinのpower_anova」
pingouin ライブラリを用いて各群の必要サンプルサイズを計算します。
### 別解 pingouinのpower_anovaで必要サンプルサイズ(各群)を算出
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
k = 3 # 群の数
eta_sq = 0.06 # 検出したいの効果量（相関比の二乗）の大きさ
# 必要サンプルサイズの算出
n = np.ceil(pg.power_anova(eta_squared=eta_sq, k=k, power=1-beta, alpha=alpha))
n【実行結果】
各群の必要サンプルサイズはテキストの計算結果と同じになりました！
各群のサンプルサイズ$${n_g=52}$$のときのタイプⅡエラー確率を計算します。
### 別解 続き pingouinのpower_anovaで必要サンプルサイズのタイプⅡエラー確率を算出
1 - pg.power_anova(eta_squared=eta_sq, k=k, n=n, alpha=alpha)【実行結果】
テキストの計算結果と同じになりました！
■ アディショナルタイム②「statsmodelsのFTestAnovaPower」
statsmodels ライブラリを用いて必要サンプルサイズを計算します。
効果量 effect_size の設定が上手くいっているか心配です…
### 別解 statsmodelsのFTestAnovaPowerで必要サンプルサイズ(各群)を算出
# 効果量に設定する値は、η²ではなく、平均を標準偏差でわった値とのこと
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
k = 3 # 群の数
eta_sq = 0.06 # 検出したいの効果量（相関比の二乗）の大きさ
# 必要サンプルサイズの算出
n = np.ceil(power.FTestAnovaPower().solve_power(effect_size=0.445, alpha=alpha,
power=1-beta, k_groups=k))
n【実行結果】
各群の必要サンプルサイズはテキストの計算結果と同じになりました！
各群のサンプルサイズ$${n_g=52}$$のときのタイプⅡエラー確率を計算します。
### 別解 statsmodelsのFTestAnovaPowerで必要サンプルサイズのタイプⅡエラー確率を算出
1 - power.FTestAnovaPower().solve_power(effect_size=0.445, nobs=n, alpha=alpha,
k_groups=k)【実行結果】
テキストの計算結果と概ね同じになりました！（若干、値は大きいかもです）
6.4.3 反復測定分散分析
反復測定分散分析は、比較する群間で回答する人が同じ「参加者内計画」の分散分析です（対応がある分散分析）。
反復測定分散分析の検定統計量・$${F}$$値も非心$${F}$$分布に従います。
反復測定分散分析の場合の非心$${F}$$分布の非心度は、反復測定数$${m}$$、測定間の相関係数$${\rho}$$、母効果量$${\eta^2}$$（相関比の二乗）を用いて計算できるそうです。
$$
\lambda = \cfrac{\eta^2}{1 - \eta^2} \ n \ \cfrac{m}{1 - \rho}
$$
テキストに合わせて、反復測定分散分析のタイプⅡエラー確率計算を関数化します。
上記の非心度$${\lambda}$$の数式は、関数内の「非心度の計算」で使用します。
自由度補正項$${\varepsilon}$$に注意、とのこと。
### 249ページ 反復測定分散分析のタイプⅡエラー確率を算出する関数
# alpha:有意水準, eta_sq:母効果量（相関比の二乗）, m:反復測定数, rho:測定間の相関係数
# epsilon:自由度の補正項, n:サンプルサイズ
def t2e_repeated(alpha, eta_sq, m, rho, epsilon, n):
# 自由度の計算
df1 = m - 1 # 要因の自由度
df2 = (n - 1) * (m - 1) # 誤差項の自由度
# 非心度の計算
lambda_ = eta_sq / (1 - eta_sq) * n * m / (1 - rho)
# 臨界値の計算
cv = stats.f.ppf(q=1 - alpha, dfn=df1 * epsilon, dfd=df2 * epsilon)
# タイプⅡエラー確率の計算
type2_error = stats.ncf.cdf(
x=cv, dfn=df1 * epsilon, dfd=df2 * epsilon, nc=lambda_ * epsilon)
return type2_error
■ 反復測定分散分析のサンプルサイズをシミュレーション
４群で検出したい効果量を$${0.06}$$、有意水準$${\alpha=5\%}$$、タイプⅡエラー確率の閾値$${\beta=20\%}$$の条件で、必要サンプルサイズを計算します。
シミュレーション回数の上限を10000回としています。
### 249～250ページ 反復測定分散分析のサンプルサイズ設計
# 検出したい効果量をη²=0.06で設定したとき、必要サンプルサイズは
# 有意水準α=0.05, タイプⅡエラー確率β=0.20（検出力=0.80）の下で27
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
m = 4 # 測定数
eta_sq = 0.06 # 検出したいの効果量（相関比の二乗）の大きさ
rho = 0.5 # 測定間の相関係数
epsilon = 0.8 # 自由度補正項
iter = 10000 # シミュレーション回数の最大値
## シミュレーション
# タイプⅡエラー確率が閾値β以下になるまでサンプルサイズを増やす（上限10000）
for n in range(5, iter + 1):
# サンプルサイズnのときのタイプⅡエラー確率を算出
type2_error = t2e_repeated(alpha, eta_sq, m, rho, epsilon, n)
# タイプⅡエラー確率が閾値β以下になったらシミュレーションを終了
if type2_error <= beta:
break
## 結果
# 必要サンプルサイズ
print(f'必要サンプルサイズ: {n}')
# そのときのタイプⅡエラー確率
print(f'タイプⅡエラー確率: {type2_error:.7f}')【実行結果】
テキストと同じ値になりました。
アディショナルタイム３
テキストから少々離れて、ここからはアディショナルタイムに突入です！
Pythonのライブラリを用いて反復測定分散分析のサンプルサイズ計算を行います。
■ アディショナルタイム「pingouinのpower_rm_anova」
pingouin ライブラリを用いて各群の必要サンプルサイズを計算します。
### 別解 pingouinのpower_rm_anovaで必要サンプルサイズ(各群)を算出
## 設定と準備
alpha = 0.05 # 有意水準
beta = 0.20 # タイプⅡエラー確率の閾値
m = 4 # 測定数
eta_sq = 0.06 # 検出したいの効果量（相関比の二乗）の大きさ
rho = 0.5 # 測定間の相関係数
epsilon = 0.8 # 自由度補正項
n = np.ceil(pg.power_rm_anova(eta_squared=eta_sq, m=m, power=1-beta,
alpha=alpha, corr=rho, epsilon=epsilon))
n【実行結果】
各群の必要サンプルサイズはテキストの計算結果と同じになりました！
各群のサンプルサイズ$${n=27}$$のときのタイプⅡエラー確率を計算します。
### 別解 続き pingouinのpower_rm_anovaで必要サンプルサイズのタイプⅡエラー確率を算出
1- pg.power_rm_anova(eta_squared=eta_sq, m=m, n=n, alpha=alpha, corr=rho,
epsilon=epsilon)【実行結果】
テキストの計算結果と同じになりました！
statsmodels ライブラリに関しては、反復測定分散分析の必要サンプルサイズを計算の情報を見つけられませんでした・・・
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。