はじめに
このブログ記事では、永田靖著『サンプルサイズの決め方』（朝倉書店）の7章を参考に、2つの母平均の差の検定、検出力の計算方法、適切なサンプルサイズの設計方法について解説します。これら統計学の基本的な概念をPythonを用いて実装しながら確認していきます。
2つの母平均の差の検定
2つの母平均の差の検定をする手順は書籍では以下のように紹介されています。
-
手順1 仮説の設定
- 帰無仮説
H₀: μ_1 = μ_2 μ₀ - 対立仮説:
- (1)
H₁: μ_1 ≠ μ_2 - (2)
H₁: μ_1 > μ_2 - (3)
H₁: μ_1 < μ_2
- (1)
- 帰無仮説
-
手順2 有意水準の設定通常は α = 0.05 を使用します。
- 有意水準
\alpha \alpha = 0.05
- 有意水準
-
手順3 棄却域の設定
- 対立仮説(1)の場合:
R: |t₀| ≥ t(\phi, α) - 対立仮説(2)の場合:
R: t₀ ≥ t(\phi, 2α) - 対立仮説(3)の場合:
R: t₀ ≤ -t(\phi, 2α)
- 対立仮説(1)の場合:
ここで、
-
手順4 検定統計量の計算
- データ
x_{11}, x_{12}, \ldots, x_{1n_1} x_{21}, x_{22}, \ldots, x_{2n_2} t₀ t₀ = \frac{(\bar{x}_1 - \bar{x}_2)}{\sqrt{V(1/n_1+1/n_2)}} V = \frac{S_1+S_2}{(n_1+n_2-2)} S_1 = \Sigma^{n_1}_{i=1}(x_{1i}-\bar{x}_1)^2 = \Sigma{x}^2_{1i} - \frac{(\Sigma x_{1i})^2}{n_1} S_2 = \Sigma^{n_2}_{i=1}(x_{2i}-\bar{x}_2)^2 = \Sigma{x}^2_{2i} - \frac{(\Sigma x_{2i})^2}{n_2} \phi = (n_1-1)+(n_2-1)=n_1+n_2-2
- データ
-
手順5 判定
- 計算された
t₀ R H₀
- 計算された
この検定手順に則って検定を行う関数をpythonで実装しました。
import numpy as np
from scipy import stats
def two_sample_t_test(data1, data2, alpha=0.05, test_type='two-sided'):
"""
パラメータ:
data1, data2 : array
2つの標本データ
alpha : float, optional
有意水準 (デフォルト: 0.05)
test_type : {'two-sided', 'greater', 'less'}, optional
検定の種類 (デフォルト: 'two-sided')
戻り値:
dict
検定結果を含む辞書
"""
# データの準備
x1 = np.array(data1)
x2 = np.array(data2)
n1 = len(x1)
n2 = len(x2)
# 検定統計量の計算
mean1 = np.mean(x1)
mean2 = np.mean(x2)
s1 = np.sum((x1 - mean1)**2)
s2 = np.sum((x2 - mean2)**2)
v = (s1 + s2) / (n1 + n2 - 2)
t0 = (mean1 - mean2) / np.sqrt(v * (1/n1 + 1/n2))
df = n1 + n2 - 2
# 棄却域の設定
if test_type == 'two-sided':
t_crit = stats.t.ppf(1 - alpha/2, df)
p_value = 2 * (1 - stats.t.cdf(abs(t0), df))
elif test_type == 'greater':
t_crit = stats.t.ppf(1 - alpha, df)
p_value = 1 - stats.t.cdf(t0, df)
elif test_type == 'less':
t_crit = -stats.t.ppf(1 - alpha, df)
p_value = stats.t.cdf(t0, df)
else:
raise ValueError("test_type must be 'two-sided', 'greater', or 'less'")
# 判定
if test_type == 'two-sided':
reject = abs(t0) >= t_crit
elif test_type == 'greater':
reject = t0 >= t_crit
else: # 'less'
reject = t0 <= t_crit
return {
't_statistic': t0,
'critical_value': t_crit,
'reject_null': reject,
}
# 問題7.1
data1 = [6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3]
data2 = [5.3, 6.2, 5.9, 7.3, 8.4, 7.3, 6.9, 7.6, 8.5, 8.1, 6.7, 7.7]
result = two_sample_t_test(data1, data2, alpha=0.05, test_type='two-sided')
print(f"7.1: t統計量 = {result['t_statistic']:.4f}")
print(f" 臨界値 = {result['critical_value']:.4f}")
print(f" 帰無仮説棄却 = {result['reject_null']}")
結果
7.1: t統計量 = -1.7857
臨界値 = 2.0860
帰無仮説棄却 = False
# 問題7.2
data1 = [10.8, 11.2, 9.7, 9.9, 12.0, 9.6, 10.5, 10.7, 10.1]
data2 = [10.2, 10.1, 9.9, 8.2, 10.2, 9.4, 10.4, 10.0]
result = two_sample_t_test(data1, data2, alpha=0.05, test_type='greater')
print(f"7.2: t統計量 = {result['t_statistic']:.4f}")
print(f" 臨界値 = {result['critical_value']:.4f}")
print(f" 帰無仮説棄却 = {result['reject_null']}")
結果
7.2: t統計量 = 1.9274
臨界値 = 1.7531
帰無仮説棄却 = True
# 問題7.3
data1 = [21, 19, 16, 19, 22, 18, 20, 21]
data2 = [19, 22, 21, 22, 25, 19, 24, 23, 19, 22]
result = two_sample_t_test(data1, data2, alpha=0.05, test_type='less')
print(f"7.3: t統計量 = {result['t_statistic']:.4f}")
print(f" 臨界値 = {result['critical_value']:.4f}")
print(f" 帰無仮説棄却 = {result['reject_null']}")
結果
7.3: t統計量 = -2.1732
臨界値 = -1.7459
帰無仮説棄却 = True
検出力の計算
検出力を計算する手順は書籍では以下のように紹介されています。
となります。しかし、非心
import numpy as np
from scipy import stats
def calculate_power(n1, n2, mu1, mu2, sigma, alpha=0.05, test_type='two-sided'):
"""
パラメータ:
n1, n2 (int): 各グループのサンプルサイズ
mu1, mu2 (float): 各グループの母平均
sigma (float): 共通の母標準偏差
alpha (float): 有意水準（デフォルトは0.05）
test_type (str): 検定の種類（'two-sided', 'greater', 'less'）
戻り値:
float: 計算された検出力
"""
# 自由度
phi = n1 + n2 - 2
# 非心パラメータ
delta = (mu1 - mu2) / sigma
lambda_param = np.sqrt(n1 * n2 / (n1 + n2)) * delta
# 臨界値の計算
if test_type == 'two-sided':
cv = stats.t.ppf(1 - alpha/2, phi)
elif test_type in ['greater', 'less']:
cv = stats.t.ppf(1 - alpha, phi)
else:
raise ValueError("無効な検定タイプです。'two-sided', 'greater', 'less'のいずれかを指定してください。")
# 非心t分布の正規分布による近似
def noncentral_t_approx(w):
return stats.norm.cdf((w * (1 - 1/(4*phi)) - lambda_param) / np.sqrt(1 + w**2/(2*phi)))
# 検出力の計算
if test_type == 'two-sided':
power = 1 - noncentral_t_approx(cv) + noncentral_t_approx(-cv)
elif test_type == 'greater':
power = 1 - noncentral_t_approx(cv)
else: # 'less'
power = noncentral_t_approx(-cv)
return power
# 問題7.4
power_7_4_1 = calculate_power(10, 12, 0.5, 0, 1, alpha=0.05, test_type='two-sided')
power_7_4_2 = calculate_power(20, 22, 0.5, 0, 1, alpha=0.05, test_type='two-sided')
print(f"7.4: (1) 検出力 = {power_7_4_1:.4f}")
print(f" (2) 検出力 = {power_7_4_2:.4f}")
結果
7.4: (1) 検出力 = 0.1995
(2) 検出力 = 0.3520
サンプルサイズを増やすことで検出力が0.1995から0.3520に上昇しました。
# 問題7.5
power_7_5_1 = calculate_power(9, 8, 1.0, 0, 1, alpha=0.05, test_type='greater')
power_7_5_2 = calculate_power(9, 8, 1.5, 0, 1, alpha=0.05, test_type='greater')
print(f"7.5: (1) 検出力 = {power_7_5_1:.4f}")
print(f" (2) 検出力 = {power_7_5_2:.4f}")
結果
7.5: (1) 検出力 = 0.6249
(2) 検出力 = 0.9029
# 問題7.6
power_7_6_1 = calculate_power(8, 10, -1.0, 0, 1, alpha=0.05, test_type='less')
power_7_6_2 = calculate_power(15, 17, -1.0, 0, 1, alpha=0.05, test_type='less')
print(f"7.6: (1) 検出力 = {power_7_6_1:.4f}")
print(f" (2) 検出力 = {power_7_6_2:.4f}")
結果
7.6: (1) 検出力 = 0.6451
(2) 検出力 = 0.8672
同じ
サンプルサイズの設計方法
サンプルサイズの設計方法は書籍では以下のように紹介されています。
上記の式通りにサンプルサイズを計算する関数をpythonで実装しました。
この関数を用いて練習問題を解いていきます。
import math
import numpy as np
from scipy import stats
def calculate_sample_size(alpha, beta, delta, test_type='two-sided'):
"""
パラメータ:
alpha (float): 有意水準
beta (float): 第二種の誤りの確率 (1 - 検出力)
delta (float): 効果量
test_type (str): 検定の種類（'two-sided', 'greater', 'less'）
戻り値:
int: 必要なサンプルサイズ（切り上げ）
"""
if test_type == 'two-sided':
z_alpha = stats.norm.ppf(1 - alpha/2)
else:
z_alpha = stats.norm.ppf(1 - alpha)
z_beta = stats.norm.ppf(beta)
n = 2 * ((z_alpha - z_beta) / delta) ** 2 + (z_alpha ** 2) / 4
return math.ceil(n)
def find_optimal_sample_size(alpha, target_power, delta, test_type='two-sided'):
"""
パラメータ:
alpha (float): 有意水準
target_power (float): 目標とする検出力
delta (float): 効果量
test_type (str): 検定の種類（'two-sided', 'greater', 'less'）
戻り値:
tuple: (最適なサンプルサイズ, 実際の検出力)
"""
beta = 1 - target_power
n = calculate_sample_size(alpha, beta, delta, test_type)
while calculate_power(n, n, delta, 0, 1, alpha, test_type) < target_power:
n += 1
while calculate_power(n-1, n-1, delta, 0, 1, alpha, test_type) >= target_power:
n -= 1
actual_power = calculate_power(n, n, delta, 0, 1, alpha, test_type)
return n, actual_power
# 問題7.7
print("問題7.7:")
n_7_7_1, power_7_7_1 = find_optimal_sample_size(0.05, 0.90, 0.8, 'two-sided')
print(f"(1) サンプルサイズ: {n_7_7_1}, 実際の検出力: {power_7_7_1:.4f}")
n_7_7_2, power_7_7_2 = find_optimal_sample_size(0.05, 0.95, 0.8, 'two-sided')
print(f"(2) サンプルサイズ: {n_7_7_2}, 実際の検出力: {power_7_7_2:.4f}")
結果
問題7.7:
(1) サンプルサイズ: 34, 実際の検出力: 0.9015
(2) サンプルサイズ: 42, 実際の検出力: 0.9518
検出力を0.90から0.95に上げるために、サンプルサイズを34から42に増やす必要がありました。高い検出力を得るには、より大きなサンプルサイズが必要となることがわかります。
# 問題7.8
print("\n問題7.8:")
n_7_8_1, power_7_8_1 = find_optimal_sample_size(0.05, 0.90, 1.0, 'greater')
print(f"(1) サンプルサイズ: {n_7_8_1}, 実際の検出力: {power_7_8_1:.4f}")
n_7_8_2, power_7_8_2 = find_optimal_sample_size(0.05, 0.90, 0.5, 'greater')
print(f"(2) サンプルサイズ: {n_7_8_2}, 実際の検出力: {power_7_8_2:.4f}")
結果
問題7.8:
(1) サンプルサイズ: 18, 実際の検出力: 0.9023
(2) サンプルサイズ: 70, 実際の検出力: 0.9030
# 問題7.9
print("\n問題7.9:")
n_7_9_1, power_7_9_1 = find_optimal_sample_size(0.05, 0.90, -1.2, 'less')
print(f"(1) サンプルサイズ: {n_7_9_1}, 実際の検出力: {power_7_9_1:.4f}")
n_7_9_2, power_7_9_2 = find_optimal_sample_size(0.05, 0.95, -1.5, 'less')
print(f"(2) サンプルサイズ: {n_7_9_2}, 実際の検出力: {power_7_9_2:.4f}")
結果
問題7.9:
(1) サンプルサイズ: 13, 実際の検出力: 0.9077
(2) サンプルサイズ: 11, 実際の検出力: 0.9600
最後に
本ブログでは、永田靖著『サンプルサイズの決め方』の第7章「2つの母平均の差の検定」の内容を、Pythonを用いて実践的に解説しました。ここで学んだ主要なポイントは以下の通りです：
- 2つの母平均の差の検定の基本概念と手順
- 検出力の計算方法とその重要性
- 適切なサンプルサイズの設計方法
次回のブログでは、『サンプルサイズの決め方』の第8章の内容を取り上げる予定です。このブログを通して少しでもサンプルサイズの設計方法について理解が進めば幸いです。
Discussion