출판 편향과 파일서랍 문제¶

"""승자의 저주: 유의한 연구만 출판되면 문헌의 효과크기가 부풀려진다."""

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n_studies = 2000
n_per_group = 25          # 집단당 25명 — 작고 검정력이 낮은 연구
true_effect = 0.3         # 코헨의 d. 실재하지만 크지 않은 효과

# === 독립적인 연구 2000개를 돌린다 ===
# 중요한 전제: 2000개 모두 참 효과 0.3인 같은 현상을 연구한다.
# 즉 연구자들의 실력도, 현상도 모두 같다. 다른 것은 우연뿐이다.
effects, pvals = [], []
for _ in range(n_studies):
    control = rng.normal(0, 1, n_per_group)
    treated = rng.normal(true_effect, 1, n_per_group)
    t_stat, p = stats.ttest_ind(treated, control)

    # 코헨의 d = (평균 차이) / (합동표준편차)
    # 단위에 의존하지 않는 표준화된 효과크기라 연구끼리 비교할 수 있다.
    pooled_sd = np.sqrt((control.var(ddof=1) + treated.var(ddof=1)) / 2)
    effects.append((treated.mean() - control.mean()) / pooled_sd)
    pvals.append(p)

effects, pvals = np.array(effects), np.array(pvals)

# === 출판 필터 ===
# p < 0.05 인 연구만 학술지에 실린다고 두자.
# 검정력이 낮으므로, 통과하려면 우연히 효과가 크게 나와야 한다.
# 다시 말해 **필터가 효과크기가 큰 표본을 골라낸다.** 이것이 승자의 저주다.
published = pvals < 0.05

# === 서랍 속에 남은 것과 문헌에 보이는 것을 비교한다 ===
print(f"True effect                : {true_effect:.2f}")
print(f"All studies,      mean d   : {effects.mean():.2f}  (n = {n_studies})")
print(f"Published only,   mean d   : {effects[published].mean():.2f}  "
      f"(n = {published.sum()})")
print(f"Inflation                  : "
      f"{100 * (effects[published].mean() / true_effect - 1):.0f}%")
print(f"Studies in the file drawer : {(~published).sum()}")

# --- (1)(2) 의 이론값과 맞추어 본다 ---
from scipy.integrate import quad

df = 2 * n_per_group - 2
ncp = true_effect * np.sqrt(n_per_group / 2)
t_c = stats.t.ppf(0.975, df)
power = stats.nct.sf(t_c, df, ncp) + stats.nct.cdf(-t_c, df, ncp)
k = np.sqrt(2 / n_per_group)

f = lambda t: stats.nct.pdf(t, df, ncp)
num = quad(lambda t: t * f(t), t_c, 50)[0] + quad(lambda t: t * f(t), -50, -t_c)[0]
den = quad(f, t_c, 50)[0] + quad(f, -50, -t_c)[0]

print(f"이론 검정력 = {power:.4f}  (관측 {published.mean():.4f})")
print(f"출판 문턱 |d| > {k * t_c:.4f}")
print(f"이론 E[d | 유의] = {k * num / den:.4f}  (관측 {effects[published].mean():.4f})")
print(f"관측 전체 평균 d = {effects.mean():.4f},  "
      f"SE = {effects.std(ddof=1) / np.sqrt(n_studies):.4f}")

True effect                : 0.30
All studies,      mean d   : 0.32  (n = 2000)
Published only,   mean d   : 0.72  (n = 381)
Inflation                  : 141%
Studies in the file drawer : 1619
이론 검정력 = 0.1800  (관측 0.1905)
출판 문턱 |d| > 0.5687
이론 E[d | 유의] = 0.7232  (관측 0.7231)
관측 전체 평균 d = 0.3178,  SE = 0.0064

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

# 축이름과 제목이 한글이므로 한글 글꼴을 지정한다.
plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

rng = np.random.default_rng(0)

def run_studies(K=400, d=0.3):
    """표본 크기가 제각각인 연구 K 편을 돌린다."""
    ns = rng.integers(10, 300, K)
    eff = np.empty(K); se = np.empty(K); pv = np.empty(K)
    for i, n in enumerate(ns):
        c = rng.normal(0, 1, n)
        t = rng.normal(d, 1, n)
        sd = np.sqrt((c.var(ddof=1) + t.var(ddof=1)) / 2)
        eff[i] = (t.mean() - c.mean()) / sd
        se[i] = np.sqrt(2 / n)               # 코헨 d 의 근사 표준오차
        pv[i] = stats.ttest_ind(t, c)[1]
    return eff, se, pv

eff, se, pv = run_studies()
published = pv < 0.05

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True, sharey=True)
for ax, mask, title in [(axes[0], np.ones(len(eff), bool), "수행된 모든 연구 — 대칭"),
                        (axes[1], published, "유의한 연구만 — 왼쪽 아래가 비었다")]:
    ax.scatter(eff[mask], 1 / se[mask], s=14, alpha=0.6)
    ax.axvline(0.3, color="red", ls="--", lw=1.5, label="참 효과 0.3")
    ax.set_xlabel("효과크기 (코헨의 d)")
    ax.set_title(title, fontsize=10)
    ax.legend(fontsize=8)
axes[0].set_ylabel("정밀도 (1 / SE)")
fig.tight_layout()
plt.show()

for label, mask in [("전체", np.ones(len(eff), bool)), ("출판된 것만", published)]:
    pooled = np.average(eff[mask], weights=1 / se[mask] ** 2)
    print(f"{label:>10}: {mask.sum():>3}편   통합 효과 {pooled:.4f}")

        전체: 400편   통합 효과 0.3066
    출판된 것만: 275편   통합 효과 0.3375

import numpy as np
import statsmodels.api as sm
from scipy import stats

rng = np.random.default_rng(0)

def run_studies(K=400, d=0.3):
    ns = rng.integers(10, 300, K)
    eff = np.empty(K); se = np.empty(K); pv = np.empty(K)
    for i, n in enumerate(ns):
        c = rng.normal(0, 1, n)
        t = rng.normal(d, 1, n)
        sd = np.sqrt((c.var(ddof=1) + t.var(ddof=1)) / 2)
        eff[i] = (t.mean() - c.mean()) / sd
        se[i] = np.sqrt(2 / n)
        pv[i] = stats.ttest_ind(t, c)[1]
    return eff, se, pv

eff, se, pv = run_studies()

print(f"{'문헌':>12}{'연구 수':>9}{'통합 효과':>11}{'Egger 절편':>12}{'p 값':>12}")
for label, mask in [("전체", np.ones(len(eff), bool)), ("유의한 것만", pv < 0.05)]:
    y = eff[mask] / se[mask]
    X = sm.add_constant(1 / se[mask])
    fit = sm.OLS(y, X).fit()
    pooled = np.average(eff[mask], weights=1 / se[mask] ** 2)
    print(f"{label:>12}{mask.sum():>9}{pooled:>11.4f}"
          f"{fit.params[0]:>+12.3f}{fit.pvalues[0]:>12.2e}")

          문헌     연구 수      통합 효과    Egger 절편         p 값
          전체      400     0.3066      +0.026    8.73e-01
      유의한 것만      275     0.3375      +1.925    7.22e-16

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

print(f"{'상황':>14}{'유의 건수':>11}   p 구간별 비율 [0-.01, .01-.02, .02-.03, .03-.04, .04-.05]")
for d, label in [(0.0, "참 효과 없음"), (0.5, "참 효과 0.5")]:
    ps = np.array([stats.ttest_ind(rng.normal(d, 1, 25), rng.normal(0, 1, 25))[1]
                   for _ in range(20_000)])
    sig = ps[ps < 0.05]
    counts = np.histogram(sig, bins=[0, .01, .02, .03, .04, .05])[0]
    print(f"{label:>14}{len(sig):>11}   {np.round(counts / len(sig), 3)}")

            상황      유의 건수   p 구간별 비율 [0-.01, .01-.02, .02-.03, .03-.04, .04-.05]
       참 효과 없음       1026   [0.202 0.189 0.199 0.19  0.22 ]
      참 효과 0.5       8230   [0.473 0.193 0.133 0.109 0.092]

import numpy as np
from scipy import stats

alpha = 0.05
print(f"{'사전 오즈 R':>12}{'검정력 0.8':>12}{'검정력 0.5':>12}{'검정력 0.2':>12}")
for R, desc in [(1.0, "확증 시험"), (0.2, "탐색적 연구"),
                (0.1, "후보 유전자"), (0.02, "무작위 탐색")]:
    row = [pw * R / (pw * R + alpha) for pw in (0.8, 0.5, 0.2)]
    print(f"{R:>12.2f}" + "".join(f"{v:>12.3f}" for v in row) + f"   {desc}")

# 모의로 확인
rng = np.random.default_rng(0)
R, power, N = 0.1, 0.2, 400_000
true = rng.random(N) < R / (1 + R)
sig = np.where(true, rng.random(N) < power, rng.random(N) < alpha)
print(f"\n모의 확인 (R=0.1, 검정력 0.2): 유의 결과 중 참인 비율 "
      f"{true[sig].mean():.4f}   이론 {power * R / (power * R + alpha):.4f}")

     사전 오즈 R     검정력 0.8     검정력 0.5     검정력 0.2
        1.00       0.941       0.909       0.800   확증 시험
        0.20       0.762       0.667       0.444   탐색적 연구
        0.10       0.615       0.500       0.286   후보 유전자
        0.02       0.242       0.167       0.074   무작위 탐색

모의 확인 (R=0.1, 검정력 0.2): 유의 결과 중 참인 비율 0.2830   이론 0.2857

def ppv(R, power, u, alpha=0.05):
    beta = 1 - power
    return ((1 - beta) * R + u * beta * R) / (
        R + alpha - beta * R + u - u * alpha + u * beta * R)

print(f"{'편향 u':>8}{'PPV':>10}   (R = 0.1, 검정력 0.2)")
for u in (0.0, 0.1, 0.3, 0.5):
    print(f"{u:>8.1f}{ppv(0.1, 0.2, u):>10.3f}")

    편향 u       PPV   (R = 0.1, 검정력 0.2)
     0.0     0.286
     0.1     0.162
     0.3     0.116
     0.5     0.103