p-해킹 시연¶

import numpy as np
from scipy import stats

np.random.seed(42)

n_experiments = 10_000
n_per_group = 30
pvals = np.zeros(n_experiments)

for i in range(n_experiments):
    a = np.random.normal(0, 1, n_per_group)
    b = np.random.normal(0, 1, n_per_group)
    _, pvals[i] = stats.ttest_ind(a, b)

# 두 집단 모두 같은 N(0,1)에서 뽑았으니 H0가 참인 상황이다.
# 그런데도 5%는 기각된다. 그것이 alpha의 정의다.
false_pos_rate = np.mean(pvals < 0.05)
print(f"False positive rate: {false_pos_rate:.4f}  (expected: 0.05)")

# (1) 의 예측과 견준다.
count = int(np.sum(pvals < 0.05))
mean_y = n_experiments * 0.05
sd_y = np.sqrt(n_experiments * 0.05 * 0.95)
print(f"유의한 개수 {count}  예측 {mean_y:.0f} +- {sd_y:.3f}"
      f"   거리 {(count - mean_y) / sd_y:+.3f} SD")

# p-값이 정말 균등한가.
print(f"\n사분위수  관측 {np.percentile(pvals, [25, 50, 75]).round(4)}"
      f"   이론 [0.25 0.5  0.75]")
print(f"평균 {pvals.mean():.4f} (이론 0.5)"
      f"   표준편차 {pvals.std(ddof=1):.4f} (이론 {1 / np.sqrt(12):.4f})")
ks = stats.kstest(pvals, "uniform")
print(f"균등성 KS 검정: D = {ks.statistic:.4f}, p = {ks.pvalue:.4f}")

False positive rate: 0.0508  (expected: 0.05)
유의한 개수 508  예측 500 +- 21.794   거리 +0.367 SD

사분위수  관측 [0.2527 0.5023 0.7477]   이론 [0.25 0.5  0.75]
평균 0.4996 (이론 0.5)   표준편차 0.2878 (이론 0.2887)
균등성 KS 검정: D = 0.0058, p = 0.8898

n_outcomes = 20
min_pvals = np.zeros(1000)

for i in range(1000):
    ps = []
    for _ in range(n_outcomes):
        a = np.random.normal(0, 1, 30)
        b = np.random.normal(0, 1, 30)
        _, p = stats.ttest_ind(a, b)
        ps.append(p)
    min_pvals[i] = min(ps)

# 20개를 검정하고 그중 **가장 작은** p-값만 보고하는 상황을 흉내 낸다.
phack_rate = np.mean(min_pvals < 0.05)
print(f"Cherry-pick rate: {phack_rate:.4f}  (theoretical: 0.6415)")

# (1) 의 예측과 견준다.
k = n_outcomes
theory = 1 - 0.95 ** k
cnt = int(np.sum(min_pvals < 0.05))
sd = np.sqrt(1000 * theory * (1 - theory))
print(f"유의한 횟수 {cnt}  예측 {1000 * theory:.1f} +- {sd:.3f}"
      f"   거리 {(cnt - 1000 * theory) / sd:+.3f} SD")
print(f"min p 의 평균   {min_pvals.mean():.6f}  이론 1/(k+1) = {1 / (k + 1):.6f}")
print(f"min p 의 중앙값 {np.median(min_pvals):.6f}"
      f"  이론 1-2^(-1/k) = {1 - 2 ** (-1 / k):.6f}")

Cherry-pick rate: 0.6580  (theoretical: 0.6415)
유의한 횟수 658  예측 641.5 +- 15.165   거리 +1.087 SD
min p 의 평균   0.045322  이론 1/(k+1) = 0.047619
min p 의 중앙값 0.033638  이론 1-2^(-1/k) = 0.034064

n_experiments = 1000
n_max = 200
check_interval = 10

stopped_pvals = []
for _ in range(n_experiments):
    a, b = [], []
    for n in range(check_interval, n_max + 1, check_interval):
        a.extend(np.random.normal(0, 1, check_interval).tolist())
        b.extend(np.random.normal(0, 1, check_interval).tolist())
        _, p = stats.ttest_ind(a, b)
        if p < 0.05:
            stopped_pvals.append(p)
            break
    else:
        stopped_pvals.append(p)

# for-else 구문이다. break 없이 반복이 끝나면 else가 실행된다.
# 즉 20번을 다 엿봐도 유의하지 않았던 경우에는 마지막 p-값을 기록한다.
stop_rate = np.mean(np.array(stopped_pvals) < 0.05)
print(f"Optional stopping rate: {stop_rate:.4f}")

# 되풀이를 10만 번으로 늘려 기준값을 다시 잡는다.
# 반복문 대신 누적합으로 모든 엿보기의 t 통계량을 한꺼번에 만든다.
def stopping_rate(B, n_looks, step=10, alpha=0.05, seed=0):
    """step 명씩 n_looks 번 엿볼 때 '한 번이라도 기각'할 비율."""
    rng = np.random.default_rng(seed)
    n_max = step * n_looks
    a = rng.normal(0, 1, (B, n_max))
    b = rng.normal(0, 1, (B, n_max))
    ca, cb = np.cumsum(a, 1), np.cumsum(b, 1)
    ca2, cb2 = np.cumsum(a * a, 1), np.cumsum(b * b, 1)
    rejected = np.zeros(B, dtype=bool)
    for n in range(step, n_max + 1, step):
        i = n - 1
        ma, mb = ca[:, i] / n, cb[:, i] / n
        va = (ca2[:, i] - n * ma ** 2) / (n - 1)
        vb = (cb2[:, i] - n * mb ** 2) / (n - 1)
        t = (ma - mb) / np.sqrt((va + vb) / 2 * 2 / n)
        rejected |= 2 * stats.t.sf(np.abs(t), 2 * n - 2) < alpha
    return rejected.mean()

B = 100_000
print(f"\n{'엿본 횟수':>8} {'거짓 양성률':>11} {'SE':>8} {'독립이라면':>10}")
for L in (1, 2, 5, 10, 20):
    r = stopping_rate(B, L)
    print(f"{L:>8} {r:>11.4f} {np.sqrt(r * (1 - r) / B):>8.4f}"
          f" {1 - 0.95 ** L:>10.4f}")

Optional stopping rate: 0.2380

   엿본 횟수      거짓 양성률       SE      독립이라면
       1      0.0503   0.0007     0.0500
       2      0.0841   0.0009     0.0975
       5      0.1426   0.0011     0.2262
      10      0.1944   0.0013     0.4013
      20      0.2492   0.0014     0.6415

stopped = []
for _ in range(1000):
    a, b = [], []
    for n in range(2, 201):  # need at least 2 per group
        a.append(np.random.normal(0, 1))
        b.append(np.random.normal(0, 1))
        if len(a) >= 2:
            _, p = stats.ttest_ind(a, b)
            if p < 0.05:
                stopped.append(p)
                break
    else:
        stopped.append(p)
rate = np.mean(np.array(stopped) < 0.05)
print(f"Rate with every-observation peeking: {rate:.4f}")

Rate with every-observation peeking: 0.4350

import numpy as np
from scipy import stats

rng = np.random.default_rng(2024)
n, M = 30, 2_000

def study(rng):
    g = np.repeat([0, 1], n)
    y1 = rng.normal(0, 1, 2 * n)
    y2 = rng.normal(0, 1, 2 * n)
    cov = rng.normal(0, 1, 2 * n)
    ps = []
    for y in (y1, y2, (y1 + y2) / 2):
        ps.append(stats.ttest_ind(y[g == 0], y[g == 1]).pvalue)
        X = np.column_stack([np.ones(2 * n), g, cov])
        b, *_ = np.linalg.lstsq(X, y, rcond=None)
        r = y - X @ b
        df = 2 * n - 3
        se = np.sqrt((r @ r / df) * np.linalg.inv(X.T @ X)[1, 1])
        ps.append(2 * stats.t.sf(abs(b[1] / se), df))
    return min(ps)

cnt = sum(study(rng) < 0.05 for _ in range(M))
print(f"6가지 분석 중 최소 p < 0.05 인 비율: {cnt / M:.4f}")

6가지 분석 중 최소 p < 0.05 인 비율: 0.1345

import numpy as np
from scipy import stats

rng = np.random.default_rng(2024)

def peek(rng, step, nmax=100, nmin=10):
    x = rng.normal(0, 1, nmax)
    for n in range(nmin, nmax + 1, step):
        s = x[:n]
        t = s.mean() / (s.std(ddof=1) / np.sqrt(n))
        if 2 * stats.t.sf(abs(t), n - 1) < 0.05:
            return True
    return False

for step in [1, 5, 10, 50]:
    c = sum(peek(rng, step) for _ in range(5_000))
    print(f"{step:3d}개마다 확인: 거짓양성률 {c / 5_000:.4f}  "
          f"(확인 횟수 {len(range(10, 101, step))})")

  1개마다 확인: 거짓양성률 0.2894  (확인 횟수 91)
  5개마다 확인: 거짓양성률 0.2290  (확인 횟수 19)
 10개마다 확인: 거짓양성률 0.2014  (확인 횟수 10)
 50개마다 확인: 거짓양성률 0.0926  (확인 횟수 2)

rng = np.random.default_rng(7)
M = 200_000
for K in [1, 2, 3, 5, 10]:
    inc = rng.standard_normal((M, K)) / np.sqrt(K)
    z = np.cumsum(inc, axis=1) / np.sqrt(np.arange(1, K + 1) / K)
    c = np.percentile(np.abs(z).max(1), 95)
    print(f"확인 {K:2d}회: 포콕 임계값 {c:.3f}  각 시점의 명목 α "
          f"{2 * stats.norm.sf(c):.4f}")

확인  1회: 포콕 임계값 1.955  각 시점의 명목 α 0.0506
확인  2회: 포콕 임계값 2.179  각 시점의 명목 α 0.0293
확인  3회: 포콕 임계값 2.290  각 시점의 명목 α 0.0220
확인  5회: 포콕 임계값 2.416  각 시점의 명목 α 0.0157
확인 10회: 포콕 임계값 2.558  각 시점의 명목 α 0.0105

import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
M = 200_000
za = stats.norm.ppf(0.975)
bins = np.array([0.01, 0.02, 0.03, 0.04, 0.05])

print(f"{'상황':>16s} " + " ".join(f"{f'~{b:.2f}':>8s}" for b in bins))

# ① 효과 없음
p = rng.random(M)
sig = p[p < 0.05]
h = np.histogram(sig, bins=np.r_[0, bins])[0] / len(sig)
print(f"{'효과 없음':>16s} " + " ".join(f"{v:8.3f}" for v in h))

# ② 효과 있음(검정력 0.5, 0.8)
for pw in [0.50, 0.80]:
    nc = za + stats.norm.ppf(pw)
    p = 2 * stats.norm.sf(np.abs(rng.normal(nc, 1, M)))
    sig = p[p < 0.05]
    h = np.histogram(sig, bins=np.r_[0, bins])[0] / len(sig)
    print(f"{f'효과 있음(검정력 {pw:.0%})':>16s} "
          + " ".join(f"{v:8.3f}" for v in h))

              상황    ~0.01    ~0.02    ~0.03    ~0.04    ~0.05
           효과 없음    0.202    0.199    0.196    0.205    0.198
  효과 있음(검정력 50%)    0.537    0.176    0.121    0.092    0.075
  효과 있음(검정력 80%)    0.737    0.116    0.066    0.046    0.034

import numpy as np
from scipy import stats

rng = np.random.default_rng(303)
M, n, k = 20_000, 60, 10        # 결과변수 10개, 총 n명

def one_trial(rng, split=True, effect=0.0):
    g = np.repeat([0, 1], n // 2)
    y = rng.normal(0, 1, (n, k))
    y[g == 1, 0] += effect                      # 첫 변수에만 참 효과
    if not split:                               # 전체로 고르고 전체로 검정
        ps = [stats.ttest_ind(y[g == 0, j], y[g == 1, j]).pvalue
              for j in range(k)]
        return min(ps) < 0.05
    h = n // 2
    idx = rng.permutation(n)
    a, b = idx[:h], idx[h:]                     # 탐색 / 확증
    ps = [stats.ttest_ind(y[np.intersect1d(a, np.where(g == 0)[0]), j],
                          y[np.intersect1d(a, np.where(g == 1)[0]), j]).pvalue
          for j in range(k)]
    j = int(np.argmin(ps))                      # 탐색에서 고른 변수
    return stats.ttest_ind(y[np.intersect1d(b, np.where(g == 0)[0]), j],
                           y[np.intersect1d(b, np.where(g == 1)[0]), j]).pvalue < 0.05

for eff, name in [(0.0, "효과 없음(제1종 오류)"), (0.8, "첫 변수에 효과(검정력)")]:
    a = np.mean([one_trial(rng, False, eff) for _ in range(M // 10)])
    b = np.mean([one_trial(rng, True, eff) for _ in range(M // 10)])
    print(f"{name:24s} 분할 없음 {a:.4f}   분할 {b:.4f}")

효과 없음(제1종 오류)            분할 없음 0.4040   분할 0.0460
첫 변수에 효과(검정력)            분할 없음 0.9125   분할 0.3450