効果量と検出力 — サンプルサイズ設計の実践

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# (a) 効果量の直感的理解
ax = axes[0][0]
x = np.linspace(-4, 6, 500)
d_values = [0.2, 0.5, 0.8, 1.2]

for d in d_values:
    # 2つの正規分布
    pdf1 = stats.norm.pdf(x, 0, 1)
    pdf2 = stats.norm.pdf(x, d, 1)

    # 重なりの計算
    overlap = np.sum(np.minimum(pdf1, pdf2)) * (x[1] - x[0])

    ax.fill_between(x, pdf1, alpha=0.1, color='blue')
    ax.fill_between(x, pdf2, alpha=0.1, color='red')

ax.set_xlim(-3, 4)
# 効果量ごとの分布を個別にプロット
for i, d in enumerate([0.2, 0.5, 0.8]):
    offset = i * 0.03
    ax.annotate(f'd={d}', xy=(d/2, 0.35 - offset),
                fontsize=10, ha='center')

d_example = 0.8
pdf1 = stats.norm.pdf(x, 0, 1)
pdf2 = stats.norm.pdf(x, d_example, 1)
ax.plot(x, pdf1, 'b-', linewidth=2, label='Group 1')
ax.plot(x, pdf2, 'r-', linewidth=2, label=f'Group 2 (d={d_example})')
ax.fill_between(x, np.minimum(pdf1, pdf2), alpha=0.3, color='purple',
                label='Overlap region')
ax.set_xlabel('Value', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f"Cohen's d = {d_example}: Visual Meaning", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) 検出力曲線
ax = axes[0][1]
n_range = np.arange(5, 201)
d_values = [0.2, 0.3, 0.5, 0.8, 1.0]
alpha = 0.05

for d in d_values:
    powers = []
    for n in n_range:
        # 非心度パラメータ
        ncp = d * np.sqrt(n / 2)
        # 棄却点
        t_crit = stats.t.ppf(1 - alpha/2, 2*n - 2)
        # 検出力
        power = (1 - stats.nct.cdf(t_crit, 2*n-2, ncp)
                 + stats.nct.cdf(-t_crit, 2*n-2, ncp))
        powers.append(power)
    ax.plot(n_range, powers, linewidth=2, label=f'd = {d}')

ax.axhline(0.8, color='red', linewidth=1.5, linestyle='--', alpha=0.5,
           label='Power = 0.80')
ax.set_xlabel('Sample size per group $n$', fontsize=12)
ax.set_ylabel('Power', fontsize=12)
ax.set_title(f'Power Curves ($\\alpha = {alpha}$, two-sided)', fontsize=13)
ax.legend(fontsize=9, ncol=2)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)

# (c) 必要サンプルサイズの計算
ax = axes[1][0]
d_range = np.linspace(0.1, 1.5, 100)
target_powers = [0.70, 0.80, 0.90, 0.95]

for target_power in target_powers:
    z_beta = stats.norm.ppf(target_power)
    z_alpha = stats.norm.ppf(1 - alpha/2)
    n_required = 2 * (z_alpha + z_beta)**2 / d_range**2
    ax.plot(d_range, n_required, linewidth=2,
            label=f'Power = {target_power:.0%}')

ax.set_xlabel("Cohen's d", fontsize=12)
ax.set_ylabel('Required $n$ per group', fontsize=12)
ax.set_title('Required Sample Size', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 500)
ax.set_xlim(0.1, 1.5)

# (d) p値と効果量の関係
ax = axes[1][1]
np.random.seed(42)
n_sim = 1000

for n in [20, 50, 200]:
    d_sim_range = np.linspace(0, 1.5, 30)
    median_p = []
    for d in d_sim_range:
        p_vals = []
        for _ in range(200):
            g1 = np.random.normal(0, 1, n)
            g2 = np.random.normal(d, 1, n)
            _, p = stats.ttest_ind(g1, g2)
            p_vals.append(p)
        median_p.append(np.median(p_vals))

    ax.semilogy(d_sim_range, median_p, linewidth=2, label=f'n = {n}')

ax.axhline(0.05, color='red', linewidth=1.5, linestyle='--', label='$\\alpha = 0.05$')
ax.set_xlabel("True Cohen's d", fontsize=12)
ax.set_ylabel('Median p-value', fontsize=12)
ax.set_title('p-value vs Effect Size', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_ylim(1e-10, 1)

plt.tight_layout()
plt.savefig('effect_size_power.png', dpi=150, bbox_inches='tight')
plt.show()

# サンプルサイズの表
print("\n必要サンプルサイズ（各群、α=0.05、両側、Power=0.80）:")
z_b = stats.norm.ppf(0.80)
z_a = stats.norm.ppf(0.975)
for d in [0.2, 0.3, 0.5, 0.8, 1.0]:
    n_req = int(np.ceil(2 * (z_a + z_b)**2 / d**2))
    print(f"  d = {d}: n = {n_req} per group, total = {2*n_req}")