Transcribing 'Understanding Statistical Mechanisms through Numerical Simulation' in Python - Chapter 6, Section 6.3 'Practical Sample Size Design'

Chapter 6: 'Sample Size Design for Appropriate Hypothesis Testing'

Authors: Dr. Koji Kosugi, Dr. Yasunori Kinoshita, Dr. Yuji Shimizu

This article covers the Python transcription activity for Section 6.3, 'Practical Sample Size Design,' in Chapter 6, 'Sample Size Design for Appropriate Hypothesis Testing,' of the textbook 'Understanding Statistical Mechanisms through Numerical Simulation'.

In Chapter 6, we learn how to design sample sizes for the pre-study phase of surveys and experiments through simulation.
This time, we will learn about sample size design for two t-tests.
We will get friendly with the non-central t-distribution!

The basic programming syntax of R is roughly similar to that of Python.
Let's trace the details of the code characters and get closer to both R and Python!
Now, let's open the textbook and set off on a simulation journey🚀

Introduction

Introduction to the textbook 'Understanding Statistical Mechanisms through Numerical Simulation'

This series is a Python transcription of the book 'Understanding Statistical Mechanisms through Numerical Simulation: Learning Psychological Statistics by Trying with R' (Gijutsu-Hyohron Co., Ltd., referred to as the 'textbook').

The textbook was released in September 2023, and as the subtitle 'Learning Psychological Statistics by Trying with R' suggests, it is a wonderful book that practices numerical simulations useful for understanding psychological statistics using concrete code in R, which has a strong reputation for statistical processing.

Citation Notation

This article cites text and code published in the book listed in the source, and modifies the published text and code as appropriate.
[Source]
'Understanding Statistical Mechanisms through Numerical Simulation: Learning Psychological Statistics by Trying with R', First Edition, First Printing, Authors: Koji Kosugi, Yasunori Kinoshita, Yuji Shimizu, Gijutsu-Hyohron Co., Ltd.

The illustrations in this article are borrowed from 'Cute Free Material Collection Irasutoya'.
Thank you!

6.3 Practical Sample Size Design

We will practice sample size design for a one-sample t-test and an independent two-sample t-test.

When calculating the required sample size in Python, please use the pingouin library!
In the article, I have included code for sample size design using pingouin in the 'Additional Time' section.

This article writes Python code in Jupyter Notebook format (extension .ipynb).
I generally use scipy.stats for calculating characteristic values of probability distributions and numpy.random.generator for random number generation.

Import the libraries mainly used.

### インポート

# 数値計算
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'

Define the drawing function for Type I error probability and Type II error probability.

### タイプⅠエラー確率とタイプⅡエラー確率の描画関数

def plot_statistical_error(lambda_, df, alpha, i, n=None):

    # 描画領域の設定
    plt.figure(figsize=(6, 3))
    # t分布の描画
    x = np.linspace(-3, 7, 1001)
    plt.plot(x, stats.t.pdf(x, df=df), color='tab:orange', linestyle='--')

    # 帰無分布の描画
    # 臨界値から8までの値を200区切りで用意
    xx = np.linspace(stats.t.ppf(q=1 - alpha / 2, df=df), 7, i)
    yy = stats.t.pdf(xx, df=df) # xxに対応したt分布の密度を得る
    # タイプⅠエラー確率を色付け
    plt.fill_between(xx, yy, color='tab:orange', ec='gray', alpha=0.2,
                    label='タイプⅠエラー確率')

    # 非心分布の描画
    plt.plot(x, stats.nct.pdf(x, df=df, nc=lambda_), color='tab:blue')
    # -3から臨界値までの値を200区切りで用意
    xx = np.linspace(-3, stats.t.ppf(q=1 - alpha / 2, df=df), i)
    yy = stats.nct.pdf(xx, df=df, nc=lambda_) # xxに対応した非心t分布の密度を得る
    # タイプⅡエラー確率を色付け
    plt.fill_between(xx, yy, color='tab:blue', ec='gray', alpha=0.2,
                    label='タイプⅡエラー確率')

    # タイトルの設定
    if n is None:
        title = f'非心度={lambda_}のときの\nタイプⅠエラー確率とタイプⅡエラー確率'
    else:
        title = f'n={n}のときのタイプⅡエラー確率'
    plt.title(title)
    # x軸ラベル、y軸ラベル、凡例の設定
    plt.xlabel('x')
    plt.ylabel('Density')
    plt.legend()
    plt.show()

6.3.1 Sample Size Design for One-Sample t-test

We will work on sample size design for a one-sample t-test!
Following the text, we will calculate the Type II error probability using the non-central t-distribution, and design the sample size by adjusting the non-centrality parameter so that it is less than the pre-set Type II error probability threshold β.

We will proceed by following the five steps listed in the text.

■ Type II error probability for sample size n=5
We will calculate the Type II error probability with a significance level of α=5%, a Type II error probability threshold of β=20%, and an effect size of δ0=0.5.

### 239ページ 1標本のt検定のサンプルサイズ設計

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値（この値未満に抑えたい）
delta = 0.5    # 見積もった効果量

# 1.最初に定めたサンプルサイズ
n = 5
# 2.非心度λの計算
lambda_ = delta * np.sqrt(n)
# 3.1.自由度の計算
df = n -1
# 3.2.臨界値の計算: t分布のパーセント点
cv = stats.t.ppf(q=1 - alpha/2, df=df)
# 4.タイプⅡエラー確率の計算: 非心t分布の累積分布関数
type2_error = stats.nct.cdf(x=cv, df=df, nc=lambda_)
# タイプⅡエラー確率の出力
type2_error

[Execution Result]

The Type II error probability is 86%, which is significantly higher than β=20%.
Let's visualize it to check the Type I error probability and Type II error probability.

# 240ページ 図6.6
# 帰無分布（点線）の臨界値が非心分布（青い線）のエラー確率の上端になる

plot_statistical_error(lambda_, df, alpha, 200, n)

[Execution Result]

In accordance with the text, we will functionalize the Type II error probability calculation for the one-sample t-test.

### 240ページ 1標本のt検定のタイプⅡエラー確率算出関数

def t2e_ttest(alpha, delta, n):

    # 2.非心度λの計算
    lambda_ = delta * np.sqrt(n)
    # 3.1.自由度の計算
    df = n -1
    # 3.2.臨界値の計算: t分布のパーセント点
    cv = stats.t.ppf(q=1 - alpha/2, df=df)
    # 4.タイプⅡエラー確率の計算: 非心t分布の累積分布関数
    type2_error = stats.nct.cdf(x=cv, df=df, nc=lambda_)
    
    # 戻り値: タイプⅡエラー確率、描画のため非心度と自由度も返す
    return type2_error, lambda_, df

■ Type II error probability for sample sizes n=10, 15, 20
We will gradually increase the sample size to observe the changes in the Type II error probability.
The value of the Type II error probability gradually decreases.

This is the case for n=10.

### 240ページ 1標本のt検定のサンプルサイズ設計 n=10

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値（この値未満に抑えたい）
delta = 0.5    # 見積もった効果量

# n=10の場合
n = 10
type2_error, lambda_, df = t2e_ttest(alpha, delta, n)
type2_error

[Execution Result]

# 描画
plot_statistical_error(lambda_, df, alpha, 200, n)

[Execution Result]

This is the case for n=15.

### 241ページ n=15の場合
n = 15
type2_error, lambda_, df = t2e_ttest(alpha, delta, n)
type2_error

[Execution Result]

# 描画
plot_statistical_error(lambda_, df, alpha, 200, n)

[Execution Result]

This is the case for n=20.

### 241ページ n=20の場合
n = 20
type2_error, lambda_, df = t2e_ttest(alpha, delta, n)
type2_error

[Execution Result]

# 描画
plot_statistical_error(lambda_, df, alpha, 200, n)

[Execution Result]

■ Increase sample size using a for-loop until the threshold β = 20% is reached
We will increase the sample size using a for-loop (up to 10,000) to find the sample size where the Type II error probability becomes less than or equal to the threshold β = 20%.

### 241ページ 1標本のt分布のサンプルサイズ・シミュレーション

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値（この値未満に抑えたい）
delta = 0.5    # 見積もった効果量
iter = 10000   # シミュレーション回数

## シミュレーション
# タイプⅡエラー確率が閾値β以下になるまでサンプルサイズを増やす（上限10000）
for n in range(5, iter+1):
    # サンプルサイズnのときのタイプⅡエラー確率を算出
    type2_error, lambda_, df = t2e_ttest(alpha, delta, n)
    # タイプⅡエラー確率が閾値β以下になったらシミュレーションを終了
    if type2_error <= beta:
        break

## 結果
# 条件を満たすn
print(f'必要サンプルサイズ: {n}')
# タイプⅡエラー確率
print(f'タイプⅡエラー確率: {type2_error:.7f}')

[Execution Result]
When the sample size is n=34, the Type II error probability becomes 19.2%!

Let's visualize the Type I error probability and Type II error probability for a sample size of n=34.

# 描画
plot_statistical_error(lambda_, df, alpha, 200, n)

[Execution Results]

Additional Time 1

Stepping away from the text for a moment, we are now entering additional time!
We will perform sample size calculations using Python libraries.

■ Additional Time 1: 'pingouin's power_ttest'
We will calculate the required sample size using the pingouin library.

### 別解 pingouinのpower_ttestで必要サンプルサイズを算出

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値
delta = 0.5    # 見積もった効果量

## 算出
n = np.ceil(pg.power_ttest(d=delta, power=1 - beta, alpha=alpha,
                           contrast='one-sample', alternative='two-sided'))
n

[Execution Results]
The results match the calculations in the text!

We will calculate the Type II error probability when the sample size is $${n=34}$$.

### 別解 続き pingouinのpower_ttestで必要サンプルサイズのタイプⅡエラー確率を算出

1 - pg.power_ttest(d=delta, n=n, alpha=alpha, contrast='one-sample',
                   alternative='two-sided')

[Execution Results]
The results are almost identical to the calculations in the text!

■ Additional Time 2: 'statsmodels's tt_solve_power'
We will calculate the required sample size using the statsmodels library.

### 別解 statsmodelsのtt_solve_powerで必要サンプルサイズを算出

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値
delta = 0.5    # 見積もった効果量

## 算出
n = np.ceil(power.tt_solve_power(effect_size=delta, alpha=alpha, power=1-beta,
                                 alternative='two-sided'))
n

[Execution Results]
The results match the calculations in the text!

We will calculate the Type II error probability when the sample size is $${n=34}$$.

### 別解 続き statsmodelsのtt_solve_powerで必要サンプルサイズのタイプⅡエラー確率算出
1- power.tt_solve_power(effect_size=delta, nobs=n, alpha=alpha,
                        alternative='two-sided')

[Execution Results]
The results are almost identical to the calculations in the text!

6.3.2 Independent t-test

Next, we will work on sample size design for an independent two-sample t-test!
We will borrow the formula from the text for the non-centrality parameter in the case of an independent two-sample t-test.

Following the text, we will functionalize the Type II error probability calculation for the independent two-sample t-test.

### 243ページ 対応のない2標本のt検定のタイプⅡエラー確率算出関数の定義

def t2e_ttest_ind(alpha, delta, n1, n2):
    
    # 自由度の計算
    df = n1 + n2 - 2
    # 非心度の計算
    lambda_ = delta * np.sqrt((n1 * n2) / (n1 + n2))
    # 臨界値の計算
    cv = stats.t.ppf(q=1 - alpha / 2, df=df)
    # タイプⅡエラー確率の計算
    type2_error = stats.nct.cdf(x=cv, df=df, nc=lambda_)

    return type2_error

■ Increase the sample size using a for-loop until the threshold $${\beta = 20\%}$$ is reached
We will increase the sample size using a for-loop (up to 10,000) to search for the sample size where the Type II error probability is below the threshold $${\beta = 20\%}$$.
For determining the sample size of the two groups, we will adopt the method of setting the ratio of $${n_2}$$ to $${n_1}$$, following the text.
In this simulation, the ratio is 1.

### 243ページ 対応のない2標本のt検定のサンプルサイズ・シミュレーション

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値（この値未満に抑えたい）
delta = 0.5    # 見積もった効果量
ratio = 1      # n1に対するn2の大きさを表す比率
iter = 10000   # シミュレーション回数

## シミュレーション
# タイプⅡエラー確率が閾値β以下になるまでサンプルサイズを増やす（上限10000）
for n1 in range(5, iter+1):
    # n2の算出
    n2 = np.ceil(n1 * ratio)  # np.ceil()は正の無限大への切り上げ
    # サンプルサイズn1, n2のときのタイプⅡエラー確率を算出
    type2_error = t2e_ttest_ind(alpha, delta, n1, n2)
    # タイプⅡエラー確率が閾値β以下になったらシミュレーションを終了
    if type2_error <= beta:
        break

## 結果
# 必要サンプルサイズ
print(f'必要サンプルサイズ: {n1 + n2:.0f}')
# タイプⅡエラー確率
print(f'タイプⅡエラー確率: {type2_error:.7f}')

[Execution Results]
The required sample size and Type II error probability match those in the text.
As the text says, 'We can see that a surprisingly large sample of 128 people is required'.

Additional Time 2

Moving away from the text once again, we are now entering additional time!
We will use Python libraries to perform sample size calculations for an independent two-sample t-test.

■ Additional Time 1: 'pingouin's power_ttest'
We will calculate the required sample size using the pingouin library.
The difference from the one-sample case is that we provide 'two-samples' to the contrast argument.

### 別解 pingouinのpower_ttestで必要サンプルサイズを算出
#       サンプルサイズが異なるケースでは必要サンプルサイズを算出できない感じ

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値
delta = 0.5    # 見積もった効果量

## 必要サンプルサイズの算出
n = np.ceil(pg.power_ttest(d=delta, power=1 - beta, alpha=alpha,
                           contrast='two-samples', alternative='two-sided'))
n * 2

【Execution Result】
It matches the calculation result in the text!

We will calculate the Type II error probability when the sample size is n=128.

### 別解 続き pingouinのpower_ttestで必要サンプルサイズのタイプⅡエラー確率
1 - pg.power_ttest(d=delta, n=n, alpha=alpha, contrast='two-samples',
                   alternative='two-sided')

【Execution Result】
It is almost exactly the same as the calculation result in the text!

■ Additional Time 2: 'statsmodels's tt_ind_solve_power'
We will calculate the required sample size using the statsmodels library.
The difference from the one-sample case is that we use the tt_ind_solve_power function.

### 別解 statsmodelsのtt_ind_solve_powerで必要サンプルサイズを算出

## 設定と準備
alpha = 0.05   # 有意水準
beta = 0.20    # タイプⅡエラー確率の閾値（この値未満に抑えたい）
delta = 0.5    # 見積もった効果量
ratio = 1      # n1に対するn2の大きさを表す比率

## 算出
n = np.ceil(power.tt_ind_solve_power(effect_size=delta, alpha=alpha, ratio=ratio,
                                     power=1-beta, alternative='two-sided'))
n * 2

【Execution Result】
It matches the calculation result in the text!

We will calculate the Type II error probability when the sample size is n=128.

### 別解 続き statsmodelsのtt_ind_solve_powerで
#            必要サンプルサイズのタイプⅡエラー確率を算出
1- power.tt_ind_solve_power(effect_size=delta, nobs1=n, alpha=alpha, ratio=ratio,
                            alternative='two-sided')

【Execution Result】
It is almost exactly the same as the calculation result in the text!

That concludes this copy-coding session.

Series Articles

Next Article

Previous Article

Table of Contents

Blog Introduction

I am writing seven series of articles on note.
Please feel free to take a look!

1. Relaxed Statistics

This is a blog that roughly explores probability and statistics using the Statistical Test Grade 2 problem collection as a guide.
Feel free to read it like a casual conversation. Please do take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT Version.
Sample code for Python and EXCEL is also available.

2. Experiments! Fun Bayesian Modeling 1 & 2 with PyMC Ver. 5

I will write and analyze the Bayesian models used in psychological research from the books 'Fun Bayesian Modeling' and 'Fun Bayesian Modeling 2' using PyMC Ver. 5.
Like this book, many Bayesian models are written in R + Stan.
I will explore the possibilities of PyMC and strive to make Bayesian modeling easy to practice.
Since these are familiar and easy-to-visualize themes, please try running them in PyMC and let's enjoy it together!

3. Experiments! Bayesian Modeling from Iwanami Data Science 1 with PyMC Ver. 5

I will write and analyze the Bayesian models by four Bayesians from the book 'Experiments! Iwanami Data Science Vol. 1' using PyMC Ver. 5.
This book is a great resource for learning the basics of Bayesian programming.
I feel like I've become friends with Bayesian methods by happily running PyMC models.
Everyone, please try running them in PyMC, and let's play and learn together!

4. Fun Copying: Bayesian, Python, etc.

I will blog about the results of my 'book copying activities' for Bayesian, Python, and others.
I am mainly working on translations into Python.
I hope this serves as sample code for fellow learners who are also copying code. 🍀

5. Introduction to Time Series Analysis for Psychology with R and Stan, using Python and PyMC Ver. 5

I will practice the time series analysis from the book 'Introduction to Time Series Analysis for Psychology with R and Stan' using Python and PyMC Ver. 5.
This book is packed with themes on time series analysis!
I realized the depth of time series analysis.
I will enjoy learning time series analysis with my favorite language, Python.

6. Writing about Data Science-like things

I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
There are many articles related to statistics and data science books.
Series on 'Statistics', 'Python', 'Mathematics and Python', and 'R' have been created.

7. Python Machine Learning Programming Practice Log

I wrote articles about my various thoughts while studying the book 'Python Machine Learning Programming: PyTorch & scikit-learn Edition'.
This book is a textbook for scikit-learn and PyTorch.
Please feel free to try it out if you like.

Thank you for reading until the end.