자료 수집을 설계하라 (통계학 패러다임)¶

• 처리: 새 교수법(예: 거꾸로 교실)을 한 학기 동안 적용.
• 대조: 같은 과목, 같은 학기의 표준 교수법(평소대로).
• 무작위화 단위: 학생(개인 단위 무작위화가 비현실적이면 분반). 무작위 배정은 처리군과 대조군을 능력, 동기, 배경에서 평균적으로 비교 가능하게 만든다(어느 한 번의 배정에서 균형이 잡힌다는 보장은 아니다).
• 주요 결과변수: 기말시험 점수(실험 시작 전에 사전 지정).
• 표본 크기: 탐지하고자 하는 최소 효과 크기에 근거한 검정력 분석으로 결정.
• 눈가림: 채점자는 각 학생이 어느 집단인지 모르도록 눈가림해야 한다.

층화무작위표집은 모집단을 동질적인 층으로 나누고 각 층 안에서 독립적으로 표집한다. 층 내 변동성이 전체 변동성보다 작으므로 각 층의 평균이 더 정밀하게 추정되어 더 효율적이다.

예: 두 동네로 이루어진 도시의 평균 가구 소득 추정. 한쪽은 부유하고(평균 20만 달러, 분산 작음) 다른 쪽은 소박하다(평균 4만 달러, 분산 작음). 단순무작위표본은 우연히 한쪽 동네를 과대 또는 과소 대표할 수 있어 변동성이 커진다. 동네로 층화하여 각각에서 비례 표집하면 양쪽이 모두 대표되도록 보장되어 전체 평균 추정값의 분산이 줄어든다. 층 평균이 크게 다르고 층 내 변동성이 작을수록 이득이 크다.

타당한 추론에는 무작위표본이 선호된다. 편의표본은 가장 싸지만 가장 편향되어 있고, 할당표본은 불균형의 한 원천만 통제하는 절충안이다.

검정력은 대립가설이 참일 때 귀무가설을 올바르게 기각할 확률이다. 즉 \(\text{Power} = 1 - \beta\)이며, 여기서 \(\beta\)는 제2종 오류율(거짓음성 확률)이다.

어떤 연구가 5점 차이를 탐지할 검정력이 80%라는 것은 이런 뜻이다. 참된 처리효과가 정확히 5점이라면, 그 연구가 통계적으로 유의한 결과를 낼(선택한 유의수준에서 \(H_0\)을 기각할) 확률이 80%다. 달리 말하면 실제로 존재하는 5점 효과를 탐지하지 못할 확률이 20%다.

80%라는 기준은 관례적인 것이다. 검정력은 유의수준 \(\alpha\), 효과 크기, 표본 크기 \(n\)의 세 요인에 달려 있다. \(n\)이나 \(\alpha\)를 키우면 검정력이 커지고, 더 작은 효과를 탐지하려면 더 큰 표본이 필요하다.

군당 \(n\)이 같은 이표본 \(t\)-검정에서 (표준편차 단위의) 효과 크기 \(d\)와 유의수준 \(\alpha\)에 대해 검정력은 근사적으로

이다. \(d = 0.3\), \(n = 50\), \(\alpha = 0.05\), \(z_{0.975} = 1.96\)을 넣으면

이므로 비중심 모수는 1.5다. 검정력 \(\approx P(Z > 1.96 - 1.5) + P(Z < -1.96 - 1.5) = P(Z > 0.46) + P(Z < -3.46) \approx 0.323 + 0.0003 \approx 0.32\).

검정력이 32%에 불과하므로 이 연구는 검정력이 심각하게 부족하다. \(d = 0.3\)에서 80% 검정력에 이르려면 군당 대략 \(n \approx 175\)가 필요하다. 검정력이 부족한 연구를 수행하면 제2종 오류(실제 효과를 놓침)와 승자의 저주(유의하게 나온 소수의 결과가 효과 크기를 체계적으로 과대평가함)를 모두 감수하게 된다.

갈래길의 정원(Gelman & Loken, 2013)은 연구자가 내리는 수많은 재량적 선택 — 변수 정의, 제외 기준, 변환, 하위집단 분석 — 을 가리킨다. 자료를 본 뒤에 이런 선택을 하면, 암묵적으로 이루어진 수많은 비교 중에서 하나의 "유의한" 결과를 만들어낼 수 있다. 의도적인 p-해킹이 없더라도 제1종 오류율이 부풀어 오른다.

사전등록은 자료를 보기 전에 이 선택들을 못 박는다. 분석 계획은 주요 가설, 정확한 통계 검정, 모집단, 제외 기준, 그리고 (있다면) 다중비교를 어떻게 처리할지를 지정한다. 그 계획에서 벗어난 것은 무엇이든 확증적이 아니라 탐색적인 것으로 보고해야 한다.

현실적 우려: 엄격한 사전등록은 정당한 사후 발견을 위축시킨다. 등록된 계획이 예상하지 못한 뚜렷한 패턴이 자료에서 나타나면 연구자는 이를 탐색적 결과로 보고해야 하는데, 출판 시스템은 흔히 그런 결과를 덜 우호적으로 다룬다. 해법은 혼합이다. 확증적 가설은 엄격한 유의수준으로 사전등록하고, 탐색적 분석은 탐색적임을 투명하게 밝혀 보고한다. 등록 보고서(자료 수집 전에 동료심사가 설계를 승인하는 방식)를 도입한 학술지는 귀무 결과에 불리한 출판 유인을 없애 출판 편향을 더욱 줄이고 있다.

양측 \(\alpha\), 검정력 \(1-\beta\), 군당 \(n\)인 두 표본 검정에서 필요한 표본 크기는

이다. \(n\)이 \(d^2\)에 반비례한다는 것이 핵심이다. 효과가 절반이면 표본은 네 배가 필요하다.

from statsmodels.stats.power import TTestIndPower

power = TTestIndPower()

print("검정력 0.80, alpha 0.05 양측")
print(f"{'효과크기 d':>11}{'군당 n':>10}{'n x d^2':>10}")
for d in (0.2, 0.3, 0.5, 0.8, 1.0):
    n = power.solve_power(effect_size=d, power=0.8, alpha=0.05)
    print(f"{d:>11.1f}{n:>10.1f}{n * d * d:>10.1f}")

print("\n군당 n=50 일 때의 검정력")
for d in (0.2, 0.3, 0.5, 0.8):
    print(f"  d={d}: {power.power(effect_size=d, nobs1=50, alpha=0.05):.4f}")

출력:

검정력 0.80, alpha 0.05 양측
     효과크기 d      군당 n   n x d^2
        0.2     393.4      15.7
        0.3     175.4      15.8
        0.5      63.8      15.9
        0.8      25.5      16.3
        1.0      16.7      16.7

군당 n=50 일 때의 검정력
  d=0.2: 0.1677
  d=0.3: 0.3178
  d=0.5: 0.6969
  d=0.8: 0.9773

셋째 열 \(n \times d^2\)가 \(15.7\)–\(16.7\)로 거의 일정하다. 공식의 상수 \(2(1.96+0.84)^2 = 15.7\)과 맞는다(작은 \(n\)에서 정규 근사 대신 \(t\) 분포를 쓰기 때문에 조금 커진다).

실무에서 이 관계가 뜻하는 것.

심리학·의학의 실제 효과크기가 대개 \(0.2\)–\(0.4\) 범위라는 점을 생각하면, 군당 \(50\)명 남짓인 연구가 왜 그렇게 재현되지 않는지 이 표가 설명한다. 연습문제 5의 \(d = 0.3\), \(n = 50\)은 검정력이 \(0.32\)에 불과하다. 참 효과가 있어도 세 번 중 두 번은 놓친다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
B, n, tau = 20_000, 20, 0.5
crit_p = stats.t.ppf(0.975, n - 1)
crit_i = stats.t.ppf(0.975, 2 * n - 2)

print(f"{'개체 간 상관 rho':>17}{'대응표본 (20명 x 2회)':>24}{'독립표본 (40명)':>20}")
for rho in (0.0, 0.5, 0.8, 0.95):
    s, t = np.sqrt(rho), np.sqrt(1 - rho)

    # 대응: 같은 개체를 두 번 측정한다
    sub = rng.normal(0, 1, (B, n))
    y0 = s * sub + t * rng.normal(0, 1, (B, n))
    y1 = s * sub + t * rng.normal(0, 1, (B, n)) + tau
    d = y1 - y0
    t_paired = d.mean(1) / (d.std(1, ddof=1) / np.sqrt(n))

    # 독립: 서로 다른 40명을 20명씩 두 군으로
    a = s * rng.normal(0, 1, (B, n)) + t * rng.normal(0, 1, (B, n))
    b = s * rng.normal(0, 1, (B, n)) + t * rng.normal(0, 1, (B, n)) + tau
    t_ind = stats.ttest_ind(b, a, axis=1)[0]

    print(f"{rho:>17.2f}{(np.abs(t_paired) > crit_p).mean():>24.4f}"
          f"{(np.abs(t_ind) > crit_i).mean():>20.4f}")

출력:

      개체 간 상관 rho         대응표본 (20명 x 2회)          독립표본 (40명)
             0.00                  0.3235              0.3337
             0.50                  0.5676              0.3402
             0.80                  0.9187              0.3385
             0.95                  1.0000              0.3371

독립표본의 검정력은 \(\rho\)와 무관하게 \(0.34\)로 일정하다. 대응표본은 \(\rho\)가 커질수록 급격히 좋아진다.

왜 그런가. 대응표본에서 검정하는 것은 차이 \(d_i = y_{1i} - y_{0i}\)이고

이다. \(\rho\)가 클수록 개체 고유의 변동이 차분에서 소거된다. 사람마다 기본 혈압이 크게 다르더라도, 같은 사람의 전후 차이는 훨씬 안정적이다.

\(\rho = 0\)일 때는 오히려 조금 손해다(\(0.32\) 대 \(0.33\)). 짝지어서 얻는 것이 없는데 자유도만 \(2n-2 = 38\)에서 \(n-1 = 19\)로 절반이 되기 때문이다. 짝짓기는 짝지을 이유가 있을 때만 이득이다.

일반화: 층화·블록·공변량 보정이 모두 같은 착상이다.

• 쌍둥이 연구, 전후 비교, 같은 환자의 좌우 눈 비교
• 농업 실험의 구획(밭의 비옥도 차이를 블록으로 흡수)
• 회귀에서 사전 측정값을 공변량으로 넣기(ANCOVA)

셋 다 결과와 상관된 알려진 변동 원천을 설계나 분석에서 제거하는 것이다. 무작위화는 그 변동을 평균적으로 균형 잡아 줄 뿐 없애 주지는 않는다. 짝짓기는 없앤다. \(\square\)

import numpy as np

rng = np.random.default_rng(1)
B, N = 20_000, 40
a_eff, b_eff, ab = 0.6, 0.6, 0.5        # A 효과, B 효과, 상호작용

# 요인설계: 네 조합에 N/4 씩
k = N // 4
y = {(A, Bb): rng.normal(A * a_eff + Bb * b_eff + A * Bb * ab, 1, (B, k))
     for A in (0, 1) for Bb in (0, 1)}
main_A = ((y[(1, 0)].mean(1) + y[(1, 1)].mean(1))
          - (y[(0, 0)].mean(1) + y[(0, 1)].mean(1))) / 2

# OFAT: 기준 / A만 / B만 에 N/3 씩
m = N // 3
base = rng.normal(0, 1, (B, m))
onlyA = rng.normal(a_eff, 1, (B, m))
ofat_A = onlyA.mean(1) - base.mean(1)

print(f"요인설계 (N={N}): A 주효과 {main_A.mean():.4f}"
      f"   표준오차 {main_A.std():.4f}")
print(f"OFAT     (N={N}): A 효과   {ofat_A.mean():.4f}"
      f"   표준오차 {ofat_A.std():.4f}")
print(f"→ 표준오차가 {ofat_A.std() / main_A.std():.2f}배 작다")

출력:

요인설계 (N=40): A 주효과 0.8537   표준오차 0.3182
OFAT     (N=40): A 효과   0.5981   표준오차 0.3896
→ 표준오차가 1.22배 작다

요인설계가 두 가지 면에서 이긴다.

첫째, 정밀도. 표준오차가 \(0.318\) 대 \(0.390\)으로 \(1.22\)배 낫다. 요인설계에서는 모든 관측이 두 요인의 효과 추정에 전부 쓰이기 때문이다. \(A\)의 주효과를 구할 때 \(B=0\)인 절반과 \(B=1\)인 절반을 모두 쓴다. OFAT에서는 각 비교가 자료의 삼분의 이만 쓴다. 요인 수가 늘수록 격차가 커진다.

둘째, 상호작용. OFAT는 상호작용을 원리적으로 볼 수 없다. \(A\)만 켠 조건과 \(B\)만 켠 조건만 있고 둘 다 켠 조건이 없기 때문이다. 요인설계는 \(\mathrm{AB}\) 칸이 있으므로 상호작용 \(0.5\)를 추정할 수 있다.

상호작용을 놓치는 것이 실무에서 치명적인 이유는, OFAT의 결론이 "각각 좋으니 둘 다 하면 더 좋다" 로 이어지기 쉽기 때문이다. 상호작용이 음수라면 이 추론은 틀린다.

주의: 요인설계의 주효과는 다른 요인에 걸쳐 평균한 값이다. 위 출력에서 요인설계의 A 주효과가 \(0.85\)로 나오는데, 이는 \(a_{\text{eff}} + \mathrm{ab}/2 = 0.6 + 0.25\)다. 상호작용이 있으면 "A의 효과"라는 말 자체가 B의 수준에 따라 달라지므로, 주효과만 보고하는 것은 오해를 부른다.

이것이 산업 실험계획법의 출발점이다. 요인이 많아지면 \(2^k\)가 너무 커지므로 부분요인설계로 일부 조합만 돌리고, 고차 상호작용이 작다는 가정 아래 주효과와 저차 상호작용을 추정한다. \(\square\)

착상. 귀무가설이 "처리가 어느 개체에도 아무 효과가 없다"라면, 관측된 각 값은 배정과 무관하게 정해져 있었던 것이다. 그렇다면 실제로 일어날 수 있었던 모든 배정에 대해 검정통계량을 계산해 관측값의 극단성을 재면 된다. 분포 가정이 전혀 필요 없다. 무작위화 절차 자체가 귀무분포를 만들어 준다.

import numpy as np
from itertools import combinations
from scipy import stats

y = np.array([12.1, 11.4, 15.2, 13.8, 9.9, 14.6, 10.2, 13.1, 16.4, 12.7])
treated = np.array([1, 0, 1, 1, 0, 1, 0, 0, 1, 0], bool)
obs = y[treated].mean() - y[~treated].mean()

diffs = []
for c in combinations(range(10), 5):          # 10명 중 5명을 처리군으로
    m = np.zeros(10, bool); m[list(c)] = True
    diffs.append(y[m].mean() - y[~m].mean())
diffs = np.array(diffs)

print(f"관측된 차이 {obs:.4f},  가능한 배정 {len(diffs)}가지")
print(f"무작위화 검정 p = {(np.abs(diffs) >= abs(obs)).mean():.4f}")
print(f"t 검정 p        = {stats.ttest_ind(y[treated], y[~treated])[1]:.4f}")

출력:

관측된 차이 2.9600,  가능한 배정 252가지
무작위화 검정 p = 0.0238
t 검정 p        = 0.0153

두 \(p\) 값이 비슷하다. 자료가 정규에 가까우면 \(t\) 검정이 무작위화 검정의 좋은 근사이기 때문이다. 차이는 가정이 깨질 때 드러난다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
B = 20_000
err_t = err_perm = 0
for _ in range(B):
    z = rng.standard_cauchy(16)               # 꼬리가 아주 두꺼운 자료
    g = np.zeros(16, bool); g[rng.choice(16, 8, replace=False)] = True
    obs = z[g].mean() - z[~g].mean()

    perm = np.empty(199)
    for j in range(199):
        m = np.zeros(16, bool); m[rng.choice(16, 8, replace=False)] = True
        perm[j] = z[m].mean() - z[~m].mean()

    if (np.sum(np.abs(perm) >= abs(obs)) + 1) / 200 <= 0.05:
        err_perm += 1
    if stats.ttest_ind(z[g], z[~g])[1] <= 0.05:
        err_t += 1

print(f"참 효과가 없는 코시 자료에서의 1종 오류율 (명목 0.05)")
print(f"  t 검정        {err_t / B:.4f}")
print(f"  무작위화 검정 {err_perm / B:.4f}")

출력:

참 효과가 없는 코시 자료에서의 1종 오류율 (명목 0.05)
  t 검정        0.0197
  무작위화 검정 0.0500

무작위화 검정은 \(0.0500\)으로 명목 수준을 정확히 지킨다. \(t\) 검정은 \(0.0197\)로 지나치게 보수적이며, 이는 곧 참 효과가 있을 때 놓친다는 뜻이다(검정력 손실).

핵심. 무작위화 검정의 타당성은 자료의 분포가 아니라 연구자가 실제로 동전을 던졌다는 사실에서 나온다. 관찰연구에는 이 검정을 쓸 수 없다. 배정이 무작위가 아니었으므로 귀무분포를 만들 근거가 없기 때문이다.