대응 검정과 이표본 검정 중 무엇을 쓸 것인가¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(88)
M, n = 20_000, 20

print(f"{'ρ':>6s} {'올바른 분석':>11s} {'잘못된 분석':>12s} {'상황':>22s}")
# ① 실제로는 대응인데 독립으로 분석
for rho in [-0.3, 0.0, 0.3, 0.7]:
    L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
    z = rng.standard_normal((M, n, 2)) @ L.T
    x, y = z[:, :, 0], z[:, :, 1]            # H0 참
    d = y - x
    t1 = d.mean(1) / (d.std(1, ddof=1) / np.sqrt(n))
    c_pair = np.mean(np.abs(t1) > stats.t.ppf(0.975, n - 1))
    sp = np.sqrt((x.var(1, ddof=1) + y.var(1, ddof=1)) / 2)
    t2 = (y.mean(1) - x.mean(1)) / (sp * np.sqrt(2 / n))
    c_ind = np.mean(np.abs(t2) > stats.t.ppf(0.975, 2 * n - 2))
    print(f"{rho:6.1f} {c_pair:11.4f} {c_ind:12.4f} "
          f"{'대응 자료를 독립으로':>22s}")

# ② 실제로는 독립인데 대응으로 분석 (짝이 임의)
x = rng.standard_normal((M, n))
y = rng.standard_normal((M, n))
d = y - x
t1 = d.mean(1) / (d.std(1, ddof=1) / np.sqrt(n))
print(f"\n독립 자료를 대응으로 분석: 제1종 오류율 "
      f"{np.mean(np.abs(t1) > stats.t.ppf(0.975, n - 1)):.4f}")

     ρ      올바른 분석       잘못된 분석                     상황
  -0.3      0.0488       0.0846            대응 자료를 독립으로
   0.0      0.0482       0.0472            대응 자료를 독립으로
   0.3      0.0512       0.0222            대응 자료를 독립으로
   0.7      0.0497       0.0008            대응 자료를 독립으로

독립 자료를 대응으로 분석: 제1종 오류율 0.0514

# 검정력 비교 (참 차이 0.6, ρ=0)
x = rng.standard_normal((M, n))
y = rng.standard_normal((M, n)) + 0.6
d = y - x
t1 = d.mean(1) / (d.std(1, ddof=1) / np.sqrt(n))
sp = np.sqrt((x.var(1, ddof=1) + y.var(1, ddof=1)) / 2)
t2 = (y.mean(1) - x.mean(1)) / (sp * np.sqrt(2 / n))
print(f"독립 자료(ρ=0), 참 차이 0.6:")
print(f"  올바른 독립 분석 검정력 "
      f"{np.mean(np.abs(t2) > stats.t.ppf(0.975, 2 * n - 2)):.4f}")
print(f"  대응으로 분석한 검정력   "
      f"{np.mean(np.abs(t1) > stats.t.ppf(0.975, n - 1)):.4f}")

독립 자료(ρ=0), 참 차이 0.6:
  올바른 독립 분석 검정력 0.4555
  대응으로 분석한 검정력   0.4395

import numpy as np
from scipy import stats

def simulate(delta, seed=64, M=5_000, n1=15, n2=10, n3=10, rho=0.7):
    rng = np.random.default_rng(seed)
    L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
    hit_cc = hit_ind = hit_comb = 0
    for _ in range(M):
        z = rng.standard_normal((n1, 2)) @ L.T
        xp, yp = z[:, 0], z[:, 1] + delta      # 완전한 쌍
        xo = rng.standard_normal(n2)           # X 만
        yo = rng.standard_normal(n3) + delta   # Y 만

        # ① 완전사례 대응
        hit_cc += stats.ttest_1samp(yp - xp, 0).pvalue < 0.05
        # ② 짝을 무시하고 전부 독립으로
        hit_ind += stats.ttest_ind(np.r_[yp, yo], np.r_[xp, xo],
                                   equal_var=False).pvalue < 0.05
        # ③ 역분산 가중 결합 + 새터스웨이트 자유도
        d = yp - xp
        est_p, v_p, df_p = d.mean(), d.var(ddof=1) / n1, n1 - 1
        a, b = yo.var(ddof=1) / n3, xo.var(ddof=1) / n2
        est_i, v_i = yo.mean() - xo.mean(), a + b
        df_i = v_i**2 / (a**2 / (n3 - 1) + b**2 / (n2 - 1))
        w1, w2 = 1 / v_p, 1 / v_i
        V = 1 / (w1 + w2)
        est = (w1 * est_p + w2 * est_i) * V
        c1, c2 = w1 * V, w2 * V
        df = V**2 / ((c1**2 * v_p)**2 / df_p + (c2**2 * v_i)**2 / df_i)
        hit_comb += 2 * stats.t.sf(abs(est / np.sqrt(V)), df) < 0.05
    return hit_cc / M, hit_ind / M, hit_comb / M

print(f"{'참 차이':>7s} {'완전사례 대응':>13s} {'전부 독립':>11s} {'결합 추정':>11s}")
for delta in [0.0, 0.4, 0.7]:
    cc, ind, comb = simulate(delta)
    print(f"{delta:7.1f} {cc:13.4f} {ind:11.4f} {comb:11.4f}")

   참 차이       완전사례 대응       전부 독립       결합 추정
    0.0        0.0456      0.0086      0.0552
    0.4        0.4714      0.2260      0.5626
    0.7        0.9040      0.7226      0.9522

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
n_per, M = 15, 5_000
tau, period = 1.0, 0.8          # 처리효과 1.0, 기간효과 0.8

bias_naive, bias_cross = [], []
for _ in range(M):
    # AB 군: 기간1=A, 기간2=B / BA 군: 기간1=B, 기간2=A
    sub_ab = rng.normal(0, 1, n_per)
    sub_ba = rng.normal(0, 1, n_per)
    e = lambda: rng.normal(0, 0.5, n_per)
    ab_p1 = sub_ab + 0 + e()
    ab_p2 = sub_ab + tau + period + e()
    ba_p1 = sub_ba + tau + e()
    ba_p2 = sub_ba + 0 + period + e()

    # ① 순진한 분석: AB 군만 보고 "기간2 - 기간1"
    bias_naive.append((ab_p2 - ab_p1).mean() - tau)
    # ② 교차 분석: 두 군의 개체 내 차이를 결합
    d_ab = ab_p2 - ab_p1          # = τ + 기간
    d_ba = ba_p2 - ba_p1          # = -τ + 기간
    bias_cross.append((d_ab.mean() - d_ba.mean()) / 2 - tau)

print(f"참 처리효과 {tau},  기간효과 {period}")
print(f"순진한 분석(AB 군만)  편향 {np.mean(bias_naive):+.4f}")
print(f"교차 분석            편향 {np.mean(bias_cross):+.4f}")

참 처리효과 1.0,  기간효과 0.8
순진한 분석(AB 군만)  편향 +0.7957
교차 분석            편향 -0.0032

import numpy as np
from scipy import stats

def naive_combine_level(delta, seed=101, M=5_000,
                        n_pairs=15, rho=0.7, n_broken=3):
    """②를 그대로 적용했을 때의 기각률."""
    rng = np.random.default_rng(seed)
    L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
    hit_drop = hit_comb = 0
    for _ in range(M):
        z = rng.standard_normal((n_pairs, 2)) @ L.T
        a, b = z[:, 0], z[:, 1] + delta
        broken = rng.permutation(n_pairs)[:n_broken]   # b 가 결측
        keep = np.setdiff1d(np.arange(n_pairs), broken)

        d = (b - a)[keep]                              # ① 쌍 삭제
        hit_drop += stats.ttest_1samp(d, 0).pvalue < 0.05

        # ② 짝 없는 a 를 완전사례의 b 와 대비해 결합 (자료 공유!)
        est_p, v_p, df_p = d.mean(), d.var(ddof=1) / len(d), len(d) - 1
        a_only = a[broken]
        u, v = (b[keep].var(ddof=1) / len(keep),
                a_only.var(ddof=1) / len(a_only))
        est_i, v_i = b[keep].mean() - a_only.mean(), u + v
        df_i = v_i**2 / (u**2 / (len(keep) - 1) + v**2 / (len(a_only) - 1))
        w1, w2 = 1 / v_p, 1 / v_i
        V = 1 / (w1 + w2)
        est = (w1 * est_p + w2 * est_i) * V
        c1, c2 = w1 * V, w2 * V
        df = V**2 / ((c1**2 * v_p)**2 / df_p + (c2**2 * v_i)**2 / df_i)
        hit_comb += 2 * stats.t.sf(abs(est / np.sqrt(V)), df) < 0.05
    return hit_drop / M, hit_comb / M

for delta, name in [(0.0, "수준 (H0 참)"), (0.5, "검정력")]:
    drop, comb = naive_combine_level(delta)
    print(f"{name:14s}  쌍 삭제 {drop:.4f}   결합 {comb:.4f}")

수준 (H0 참)       쌍 삭제 0.0450   결합 0.0810
검정력             쌍 삭제 0.5324   결합 0.6136

def compare_estimators(n_broken, rho, n_pairs=15, delta=0.8,
                       M=20_000, seed=101):
    rng = np.random.default_rng(seed)
    L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
    cc, ml = [], []
    for _ in range(M):
        z = rng.standard_normal((n_pairs, 2)) @ L.T
        a, b = z[:, 0], z[:, 1] + delta
        broken = rng.permutation(n_pairs)[:n_broken]
        keep = np.setdiff1d(np.arange(n_pairs), broken)
        ak, bk = a[keep], b[keep]
        cc.append((bk - ak).mean())                    # 완전사례
        beta = np.cov(ak, bk, ddof=1)[0, 1] / ak.var(ddof=1)
        alpha = bk.mean() - beta * ak.mean()
        ml.append(alpha + (beta - 1) * a.mean())       # 회귀 보정
    cc, ml = np.array(cc), np.array(ml)
    return cc.std(ddof=1), ml.std(ddof=1), ml.var(ddof=1) / cc.var(ddof=1)

print(f"{'결측 쌍':>7s} {'ρ':>5s} {'완전사례 SD':>12s} "
      f"{'ML SD':>9s} {'분산비':>8s}")
for nb in [3, 6, 9]:
    for rho in [0.3, 0.7, 0.9]:
        s_cc, s_ml, ratio = compare_estimators(nb, rho)
        print(f"{nb:7d} {rho:5.1f} {s_cc:12.4f} {s_ml:9.4f} {ratio:8.4f}")

   결측 쌍     ρ      완전사례 SD     ML SD      분산비
      3   0.3       0.3401    0.3305   0.9441
      3   0.7       0.2230    0.2216   0.9876
      3   0.9       0.1289    0.1295   1.0093
      6   0.3       0.3921    0.3745   0.9121
      6   0.7       0.2573    0.2578   1.0034
      6   0.9       0.1489    0.1524   1.0470
      9   0.3       0.4839    0.4667   0.9299
      9   0.7       0.3170    0.3315   1.0938
      9   0.9       0.1831    0.1985   1.1753

import numpy as np
from scipy import stats

rng = np.random.default_rng(55)
M, n = 20_000, 20

c_fixed = c_adaptive = 0
for _ in range(M):
    x = rng.standard_normal(n)
    y = rng.standard_normal(n)             # H0 참, 짝 없음(ρ=0)
    # ① 사전에 대응으로 고정
    d = y - x
    c_fixed += stats.ttest_1samp(d, 0).pvalue < 0.05
    # ② 관측된 상관을 보고 고르기
    r = np.corrcoef(x, y)[0, 1]
    if r > 0.2:
        p = stats.ttest_1samp(d, 0).pvalue
    else:
        p = stats.ttest_ind(y, x).pvalue
    c_adaptive += p < 0.05
print(f"사전 고정(대응)    제1종 오류율 {c_fixed / M:.4f}")
print(f"상관 보고 선택     제1종 오류율 {c_adaptive / M:.4f}")

사전 고정(대응)    제1종 오류율 0.0513
상관 보고 선택     제1종 오류율 0.0607

c_min = 0
for _ in range(M // 2):
    x = rng.standard_normal(n)
    y = rng.standard_normal(n)
    p1 = stats.ttest_1samp(y - x, 0).pvalue
    p2 = stats.ttest_ind(y, x).pvalue
    c_min += min(p1, p2) < 0.05
print(f"두 p-값 중 작은 것 선택: 제1종 오류율 "
      f"{c_min / (M // 2):.4f}")

두 p-값 중 작은 것 선택: 제1종 오류율 0.0585

같은 개체를 두 조건에서 측정할 수 있는가
    │
    ├─ 아니오 ─────────────────────────→ 독립 설계
    │      (완치되는 질환, 비가역적 결과, 파괴검사)
    │
    └─ 예
         │
         ├─ 이월효과를 세척기간으로 제거할 수 있는가
         │      │
         │      ├─ 아니오 ──────────────→ 독립 설계
         │      │      (학습효과, 긴 반감기, 비가역적 변화)
         │      │
         │      └─ 예 ────────────────→ 교차 설계
         │
         └─ 시간이 문제라면
                └─ 짝지을 변수가 있는가 ──→ 짝지은 독립 설계
                                          (쌍둥이, 매칭, 블록)