통계적 함정: 다중 검정, p-hacking, 심슨의 역설, 상관 ≠ 인과

import numpy as np
from scipy import stats

np.random.seed(42)

n_tests = 20
n_samples = 100
significant_count = 0
n_simulations = 10000

for _ in range(n_simulations):
    p_values = []
    for _ in range(n_tests):
        # 두 그룹 모두 같은 분포에서 추출 (효과 없음)
        group_a = np.random.normal(0, 1, n_samples)
        group_b = np.random.normal(0, 1, n_samples)
        _, p = stats.ttest_ind(group_a, group_b)
        p_values.append(p)
    # 20개 중 하나라도 p < 0.05이면 "발견"
    if min(p_values) < 0.05:
        significant_count += 1

print(f"20개 검정 시 거짓 발견 비율: {significant_count/n_simulations:.1%}")
# 20개 검정 시 거짓 발견 비율: 64.0%

import numpy as np
from scipy import stats

np.random.seed(42)

n_tests = 20
n_samples = 100

# 20개 검정 중 처음 3개만 진짜 효과 있음
p_values = []
for i in range(n_tests):
    group_a = np.random.normal(0, 1, n_samples)
    if i < 3:
        group_b = np.random.normal(0.5, 1, n_samples)  # 효과 크기 0.5
    else:
        group_b = np.random.normal(0, 1, n_samples)     # 효과 없음
    _, p = stats.ttest_ind(group_a, group_b)
    p_values.append(p)

p_values = np.array(p_values)

# Bonferroni 보정
bonferroni_sig = p_values < (0.05 / n_tests)

# Benjamini-Hochberg FDR 보정
sorted_idx = np.argsort(p_values)
sorted_p = p_values[sorted_idx]
bh_threshold = np.arange(1, n_tests + 1) / n_tests * 0.05
bh_reject = sorted_p <= bh_threshold
# BH는 마지막으로 기각되는 지점까지 모두 기각
if bh_reject.any():
    max_reject_idx = np.max(np.where(bh_reject))
    bh_reject_final = np.zeros(n_tests, dtype=bool)
    bh_reject_final[sorted_idx[:max_reject_idx + 1]] = True
else:
    bh_reject_final = np.zeros(n_tests, dtype=bool)

print("검정 | p-value  | Bonferroni | BH-FDR | 진짜 효과")
print("-" * 55)
for i in range(n_tests):
    true_effect = "O" if i < 3 else "X"
    bonf = "기각" if bonferroni_sig[i] else "  -"
    bh = "기각" if bh_reject_final[i] else "  -"
    print(f" {i+1:2d}  | {p_values[i]:.4f}  |    {bonf}   |  {bh}  |    {true_effect}")

import numpy as np
from scipy import stats

np.random.seed(123)

# 효과가 전혀 없는 데이터
x = np.random.normal(50, 10, 200)
y = np.random.normal(50, 10, 200)

print("=== p-hacking 시뮬레이션 ===\n")

# 1단계: 전체 데이터로 t-검정 → 유의하지 않음
t, p = stats.ttest_ind(x, y)
print(f"1) 전체 데이터 (n=200 vs 200): p = {p:.4f}")

# 2단계: 상위 50%만 선택 → 하위 그룹에서 유의!
x_sub = x[x > np.median(x)]
y_sub = y[y > np.median(y)]
t, p = stats.ttest_ind(x_sub, y_sub)
print(f"2) 상위 50%만 (n={len(x_sub)} vs {len(y_sub)}): p = {p:.4f}")

# 3단계: 이상치 제거 (평균±1.5σ 바깥 제거)
x_mean, x_std = np.mean(x), np.std(x)
y_mean, y_std = np.mean(y), np.std(y)
x_clean = x[np.abs(x - x_mean) < 1.5 * x_std]
y_clean = y[np.abs(y - y_mean) < 1.5 * y_std]
t, p = stats.ttest_ind(x_clean, y_clean)
print(f"3) 이상치 제거 (n={len(x_clean)} vs {len(y_clean)}): p = {p:.4f}")

# 4단계: 여러 분할 시도 → 일부 구간에서 유의한 결과
print("\n4) 구간별 분석:")
best_p = 1.0
best_desc = ""
for low, high in [(0, 50), (25, 75), (50, 100), (30, 70), (40, 90), (10, 60),
                  (20, 80), (35, 85), (15, 65), (45, 95)]:
    mask_x = (x >= np.percentile(x, low)) & (x <= np.percentile(x, high))
    mask_y = (y >= np.percentile(y, low)) & (y <= np.percentile(y, high))
    if mask_x.sum() > 5 and mask_y.sum() > 5:
        t, p_sub = stats.ttest_ind(x[mask_x], y[mask_y])
        if p_sub < best_p:
            best_p = p_sub
            best_desc = f"x[{low}%-{high}%] vs y[{low}%-{high}%]"

print(f"   10개 구간 중 최소 p-value: {best_p:.4f} ({best_desc})")

# 5단계: 변수 변환 시도
transforms = {
    "log(x+100)": (np.log(x + 100), np.log(y + 100)),
    "sqrt(|x|)": (np.sqrt(np.abs(x)), np.sqrt(np.abs(y))),
    "x^2": (x**2, y**2),
    "1/x": (1/(x + 0.01), 1/(y + 0.01)),
}
print("\n5) 변수 변환 시도:")
for name, (xt, yt) in transforms.items():
    t, p_tr = stats.ttest_ind(xt, yt)
    marker = " ← 유의!" if p_tr < 0.05 else ""
    print(f"   {name}: p = {p_tr:.4f}{marker}")

import numpy as np
import pandas as pd

np.random.seed(42)

# UC Berkeley 스타일 시뮬레이션 데이터
data = {
    '학과': ['A', 'A', 'B', 'B', 'C', 'C', 'D', 'D'],
    '성별': ['남', '여', '남', '여', '남', '여', '남', '여'],
    '지원자': [825, 108, 560, 25, 325, 593, 417, 375],
    '합격자': [512, 89, 353, 17, 120, 202, 138, 131],
}
df = pd.DataFrame(data)
df['합격률'] = df['합격자'] / df['지원자']

print("=== 학과별 합격률 ===")
for dept in ['A', 'B', 'C', 'D']:
    dept_data = df[df['학과'] == dept]
    male = dept_data[dept_data['성별'] == '남'].iloc[0]
    female = dept_data[dept_data['성별'] == '여'].iloc[0]
    print(f"학과 {dept}: 남 {male['합격률']:.1%} ({male['지원자']}명) | "
          f"여 {female['합격률']:.1%} ({female['지원자']}명)")

# 전체 합산
male_total = df[df['성별'] == '남']
female_total = df[df['성별'] == '여']
male_rate = male_total['합격자'].sum() / male_total['지원자'].sum()
female_rate = female_total['합격자'].sum() / female_total['지원자'].sum()
print(f"\n=== 전체 합산 ===")
print(f"남성 합격률: {male_rate:.1%} ({male_total['지원자'].sum()}명)")
print(f"여성 합격률: {female_rate:.1%} ({female_total['지원자'].sum()}명)")
print(f"\n→ 전체: 남성이 높다. 학과별: 여성이 같거나 높다 → 심슨의 역설!")

import numpy as np
from scipy import stats

np.random.seed(42)
n = 50

# 교란 변수 Z가 X, Y 모두에 영향
z = np.random.normal(0, 1, n)  # 교란 변수 (예: 기온)
x = 2 * z + np.random.normal(0, 1, n)  # 아이스크림 판매
y = 3 * z + np.random.normal(0, 1, n)  # 범죄율

# X와 Y의 상관
r_xy, p_xy = stats.pearsonr(x, y)
print(f"X-Y 상관: r = {r_xy:.3f}, p = {p_xy:.4f}")

# Z를 통제한 편상관(Partial Correlation)
# X ~ Z의 잔차, Y ~ Z의 잔차의 상관
slope_xz = np.polyfit(z, x, 1)
slope_yz = np.polyfit(z, y, 1)
resid_x = x - np.polyval(slope_xz, z)
resid_y = y - np.polyval(slope_yz, z)
r_partial, p_partial = stats.pearsonr(resid_x, resid_y)
print(f"Z 통제 후 편상관: r = {r_partial:.3f}, p = {p_partial:.4f}")
print(f"\n→ 교란 변수 Z를 통제하면 X-Y 상관이 사라진다!")

import numpy as np

# 기저율 무시 시뮬레이션
prevalence = 0.001     # 유병률 0.1%
sensitivity = 0.99     # 민감도 (진양성률)
specificity = 0.99     # 특이도 (진음성률)

# 베이즈 정리
p_positive = sensitivity * prevalence + (1 - specificity) * (1 - prevalence)
p_disease_given_positive = (sensitivity * prevalence) / p_positive

print("=== 기저율 무시 예제 ===\n")
print(f"유병률: {prevalence:.1%}")
print(f"민감도: {sensitivity:.0%}")
print(f"특이도: {specificity:.0%}")
print(f"\n양성 판정 시 실제 감염 확률: {p_disease_given_positive:.1%}")
print(f"양성 판정 시 거짓 양성 확률: {1 - p_disease_given_positive:.1%}")

# 100만 명 검사 시뮬레이션
population = 1_000_000
actual_sick = int(population * prevalence)
actual_healthy = population - actual_sick

true_positive = int(actual_sick * sensitivity)
false_positive = int(actual_healthy * (1 - specificity))

print(f"\n=== 100만 명 검사 결과 ===")
print(f"실제 환자: {actual_sick:,}명")
print(f"진양성 (정확한 양성): {true_positive:,}명")
print(f"거짓 양성 (건강한데 양성): {false_positive:,}명")
print(f"전체 양성: {true_positive + false_positive:,}명")
print(f"양성 중 실제 환자 비율: {true_positive/(true_positive+false_positive):.1%}")

import numpy as np
from scipy import stats

np.random.seed(42)

print("=== 표본 크기와 p-value의 관계 ===\n")
print(f"{'표본 크기':>10} | {'평균 차이':>10} | {'p-value':>12} | {'Cohen d':>10} | 판정")
print("-" * 70)

# 진짜 평균 차이가 0.02 (체감 불가능한 미미한 차이, Cohen's d ≈ 0.02)
for n in [100, 1000, 10000, 100000, 1000000]:
    group_a = np.random.normal(100, 1, n)
    group_b = np.random.normal(100.02, 1, n)

    t, p = stats.ttest_ind(group_a, group_b)
    mean_diff = np.mean(group_b) - np.mean(group_a)
    pooled_std = np.sqrt((np.std(group_a)**2 + np.std(group_b)**2) / 2)
    cohens_d = mean_diff / pooled_std

    sig = "유의 *" if p < 0.05 else "비유의"
    print(f"{n:>10,} | {mean_diff:>10.4f} | {p:>12.6f} | {cohens_d:>10.4f} | {sig}")

print(f"\n→ Cohen's d ≈ 0.02 (사실상 무의미)인데, n이 커지면 p-value만 작아진다!")

import numpy as np
from scipy import stats

np.random.seed(42)

n_groups = 20
group_means_x = np.random.uniform(30, 70, n_groups)
group_means_y = 0.8 * group_means_x + np.random.normal(0, 5, n_groups)

# 집단 수준 상관
r_group, _ = stats.pearsonr(group_means_x, group_means_y)
print(f"집단 수준 상관 (n={n_groups} 그룹): r = {r_group:.3f}")

# 개인 수준 데이터 생성 — 집단 내에서는 상관이 약하거나 없음
all_x, all_y = [], []
for gx, gy in zip(group_means_x, group_means_y):
    n_ind = 50
    ind_x = np.random.normal(gx, 15, n_ind)
    ind_y = np.random.normal(gy, 15, n_ind)  # 개인 수준에서 x와 y는 독립
    all_x.extend(ind_x)
    all_y.extend(ind_y)

r_individual, _ = stats.pearsonr(all_x, all_y)
print(f"개인 수준 상관 (n={n_groups * 50} 명): r = {r_individual:.3f}")
print(f"\n→ 집단 수준에서 강한 상관이 개인 수준에서는 약해진다!")