p-해킹 시연 - 타이타닉 절단점 탐색¶

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy.stats import chi2_contingency

URL = ("https://raw.githubusercontent.com/datasciencedojo/"
       "datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv")
df = pd.read_csv(URL, index_col="PassengerId")
d = df.dropna(subset=["Age"]).copy()

age = d["Age"].to_numpy()
sur = d["Survived"].to_numpy()

# 양쪽에 최소 10명은 남는 절단점만 고려한다.
CUTS = [c for c in range(1, 80)
        if (age < c).sum() >= 10 and (age >= c).sum() >= 10]
print(f"검사 대상 절단점 {len(CUTS)}개 ({CUTS[0]}세 ~ {CUTS[-1]}세)\n")

rows = []
print(f"{'절단점':>7s}{'n<c':>6s}{'차이':>9s}{'chi2':>9s}{'p':>11s}")
for c in CUTS:
    m = age < c
    res = chi2_contingency(pd.crosstab(m, sur), correction=False)
    rows.append((c, res.statistic, res.pvalue))
    if c % 5 == 0:
        print(f"{c:>7d}{m.sum():>6d}"
              f"{sur[m].mean() - sur[~m].mean():>+9.4f}"
              f"{res.statistic:>9.3f}{res.pvalue:>11.2e}")

best = max(rows, key=lambda z: z[1])
sig = [c for c, s, p in rows if p < 0.05]
print(f"\n가장 강한 절단점 {best[0]}세: chi2 = {best[1]:.4f}, "
      f"p = {best[2]:.3e}")
print(f"p < 0.05 인 절단점 {len(sig)}개 / {len(CUTS)}개")
print(f"  {sig}")

검사 대상 절단점 64개 (2세 ~ 65세)

    절단점   n<c       차이     chi2          p
      5    40  +0.2848   12.697   3.66e-04
     10    62  +0.2264   12.032   5.23e-04
     15    78  +0.1917   10.586   1.14e-03
     20   164  +0.0981    5.038   2.48e-02
     25   278  +0.0300    0.632   4.27e-01
     30   384  +0.0002    0.000   9.96e-01
     35   479  +0.0092    0.055   8.14e-01
     40   551  +0.0414    0.893   3.45e-01
     45   599  +0.0384    0.591   4.42e-01
     50   640  +0.0461    0.584   4.45e-01
     55   672  +0.1027    1.728   1.89e-01
     60   688  +0.1421    2.098   1.48e-01
     65   703  +0.3202    4.603   3.19e-02

가장 강한 절단점 7세: chi2 = 18.2719, p = 1.915e-05
p < 0.05 인 절단점 21개 / 64개
  [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 64, 65]

from scipy.stats import chi2 as chi2dist

n = len(sur)
pbar = sur.mean()
print(f"전체 생존율 p-bar = {pbar:.6f},  N = {n}")
print(f"{'c':>4}{'n1':>5}{'n2':>5}{'차이':>11}{'공식':>11}{'표의 chi2':>11}{'차':>10}")
for c in (5, 7, 10, 20, 30, 40, 50, 60, 65):
    m = age < c
    n1, n2 = int(m.sum()), int((~m).sum())
    diff = sur[m].mean() - sur[~m].mean()
    f = diff ** 2 / (pbar * (1 - pbar) * (1 / n1 + 1 / n2))
    tab = rows[CUTS.index(c)][1]
    print(f"{c:>4}{n1:>5}{n2:>5}{diff:>+11.6f}{f:>11.5f}{tab:>11.5f}"
          f"{abs(f - tab):>10.1e}")

# 64 개 전부에서 확인
diffs = np.array([sur[age < c].mean() - sur[age >= c].mean() for c in CUTS])
n1s = np.array([(age < c).sum() for c in CUTS])
form = diffs ** 2 / (pbar * (1 - pbar) * (1 / n1s + 1 / (n - n1s)))
tabs = np.array([s for _, s, _ in rows])
print(f"\n절단점 64개 전체에서 |공식 - 표| 의 최대 = {np.abs(form - tabs).max():.3e}")

# 같은 차이가 절단점에 따라 얼마나 다른 chi2 를 주는가
print("\n차이를 0.30 으로 고정하고 n1 만 바꾸면")
print(f"{'n1':>6}{'chi2':>10}{'p':>12}")
for n1 in (20, 40, 100, 357, 600, 700):
    f = 0.30 ** 2 / (pbar * (1 - pbar) * (1 / n1 + 1 / (n - n1)))
    print(f"{n1:>6}{f:>10.4f}{chi2dist.sf(f, 1):>12.3e}")

전체 생존율 p-bar = 0.406162,  N = 714
   c   n1   n2         차이         공식    표의 chi2         차
   5   40  674  +0.284792   12.69728   12.69728   7.1e-15
   7   47  667  +0.316820   18.27192   18.27192   0.0e+00
  10   62  652  +0.226400   12.03170   12.03170   1.8e-15
  20  164  550  +0.098071    5.03758    5.03758   4.4e-15
  30  384  330  +0.000189    0.00003    0.00003   2.4e-17
  40  551  163  +0.041375    0.89278    0.89278   7.8e-16
  50  640   74  +0.046073    0.58376    0.58376   0.0e+00
  60  688   26  +0.142106    2.09760    2.09760   8.9e-16
  65  703   11  +0.320186    4.60349    4.60349   8.9e-16

절단점 64개 전체에서 |공식 - 표| 의 최대 = 8.882e-15

차이를 0.30 으로 고정하고 n1 만 바꾸면
    n1      chi2           p
    20    7.2538   7.075e-03
    40   14.0895   1.743e-04
   100   32.0882   1.473e-08
   357   66.6060   3.316e-16
   600   35.7465   2.247e-09
   700    5.1216   2.363e-02

from scipy.stats import chi2 as chi2dist

n = len(sur)
# 절단점마다의 지시벡터를 미리 쌓아 두면 한 번에 계산할 수 있다.
M = np.array([age < c for c in CUTS])           # (64, 714)

def chi2_all(y):
    """모든 절단점의 2x2 카이제곱을 한꺼번에 (닫힌 꼴)."""
    tot = y.sum()
    n1 = M.sum(axis=1)
    s1 = M @ y
    a, b = s1, n1 - s1                          # 어린 쪽 생존/사망
    c, e = tot - s1, (n - n1) - (tot - s1)      # 나이 든 쪽 생존/사망
    den = (a + b) * (c + e) * (a + c) * (b + e)
    out = np.zeros(len(CUTS))
    ok = den > 0
    out[ok] = n * (a[ok] * e[ok] - b[ok] * c[ok]) ** 2 / den[ok]
    return out

obs = chi2_all(sur).max()
print(f"관측된 최대 카이제곱 {obs:.4f}")
print(f"  이것을 단일 검정으로 읽으면 p = {chi2dist.sf(obs, 1):.4e}\n")

# 생존 여부를 무작위로 섞으면 나이와의 관계가 사라진다 (H0 가 참인 세계).
# 그 세계에서도 '최대 카이제곱'이 얼마나 커지는지 본다.
rng = np.random.default_rng(0)
B = 20_000
y = sur.copy()
null = np.empty(B)
for i in range(B):
    rng.shuffle(y)
    null[i] = chi2_all(y).max()

crit = chi2dist.ppf(0.95, 1)
print(f"귀무가설에서 '최대 카이제곱'의 분포 ({B:,}회)")
print(f"  평균 {null.mean():.4f}")
print(f"  95분위 {np.quantile(null, 0.95):.4f}   "
      f"(자유도 1 카이제곱의 95분위는 {crit:.4f})")
print(f"  99분위 {np.quantile(null, 0.99):.4f}")

print(f"\n나이와 생존이 완전히 무관한데도")
print(f"  '어떤 절단점에서든 p<0.05' 가 나올 확률 = "
      f"{np.mean(null > crit):.4f}")

print(f"\n탐색을 반영한 올바른 p = "
      f"{(np.sum(null >= obs) + 1) / (B + 1):.4f}")

관측된 최대 카이제곱 18.2719
  이것을 단일 검정으로 읽으면 p = 1.9151e-05

귀무가설에서 '최대 카이제곱'의 분포 (20,000회)
  평균 4.4309
  95분위 9.0294   (자유도 1 카이제곱의 95분위는 3.8415)
  99분위 12.1888

나이와 생존이 완전히 무관한데도
  '어떤 절단점에서든 p<0.05' 가 나올 확률 = 0.5194

탐색을 반영한 올바른 p = 0.0004

from scipy import stats

p_single = chi2dist.sf(obs, 1)
k = len(CUTS)
print(f"단일 p-값                     {p_single:.6e}")
print(f"본페로니 (k = {k})             {min(1, k * p_single):.6e}")
print(f"시데크   1-(1-p)^{k}          {1 - (1 - p_single) ** k:.6e}")
k_eff = np.log(1 - np.mean(null > crit)) / np.log(0.95)
print(f"유효 후보 수 {k_eff:.1f} 개로 보정     "
      f"{1 - (1 - p_single) ** k_eff:.6e}")
p_perm = (np.sum(null >= obs) + 1) / (B + 1)
print(f"순열이 준 p                   {p_perm:.6e}"
      f"   (초과 횟수 {int(np.sum(null >= obs))} / {B:,})")

# 20,000 회로 이만한 꼬리를 잴 수 있는가
hits = int(np.sum(null >= obs))
lo, hi = stats.beta.ppf([0.025, 0.975], hits + 0.5, B - hits + 0.5)
print(f"\n몬테카를로 정밀도")
print(f"  분해능 1/(B+1) = {1 / (B + 1):.2e}")
print(f"  초과 횟수 {hits} 에 대한 95% 구간 = [{lo:.2e}, {hi:.2e}]")
print(f"  표준오차 약 {np.sqrt(p_perm * (1 - p_perm) / B):.2e}")

단일 p-값                     1.915098e-05
본페로니 (k = 64)             1.225663e-03
시데크   1-(1-p)^64          1.224924e-03
유효 후보 수 14.3 개로 보정     2.735741e-04
순열이 준 p                   3.999800e-04   (초과 횟수 7 / 20,000)

몬테카를로 정밀도
  분해능 1/(B+1) = 5.00e-05
  초과 횟수 7 에 대한 95% 구간 = [1.57e-04, 6.87e-04]
  표준오차 약 1.41e-04

rng = np.random.default_rng(0)
B = 8_000
crit = chi2dist.ppf(0.95, 1)

print(f"명목 유의수준 0.05, 임계값 {crit:.4f}")
print(f"{'후보 수':>7s}{'실제 1종 오류':>14s}{'팽창':>8s}"
      f"{'올바른 임계값':>14s}")
for k in [1, 2, 3, 5, 10, 20, 64]:
    # 후보를 전체 범위에 고르게 k 개 배치한다.
    idx = [0] if k == 1 else \
          [int(round(i * (len(CUTS) - 1) / (k - 1))) for i in range(k)]
    sub = M[sorted(set(idx))]

    y = sur.copy()
    null_k = np.empty(B)
    for i in range(B):
        rng.shuffle(y)
        tot = y.sum()
        n1 = sub.sum(axis=1); s1 = sub @ y
        a, b = s1, n1 - s1
        c, e = tot - s1, (n - n1) - (tot - s1)
        den = (a + b) * (c + e) * (a + c) * (b + e)
        v = np.zeros(len(sub)); ok = den > 0
        v[ok] = n * (a[ok] * e[ok] - b[ok] * c[ok]) ** 2 / den[ok]
        null_k[i] = v.max()

    fwer = np.mean(null_k > crit)
    print(f"{len(sub):>7d}{fwer:>14.4f}{fwer / 0.05:>8.1f}배"
          f"{np.quantile(null_k, 0.95):>14.4f}")

명목 유의수준 0.05, 임계값 3.8415
   후보 수      실제 1종 오류      팽창       올바른 임계값
      1        0.0521     1.0배        4.1047
      2        0.1111     2.2배        4.7761
      3        0.1476     3.0배        5.6209
      5        0.2129     4.3배        6.4261
     10        0.3399     6.8배        7.5406
     20        0.4324     8.6배        8.0108
     64        0.5211    10.4배        9.1297

# 개별 절단점의 실제 수준과 양 끝 두 후보의 상관
rng = np.random.default_rng(0)
B2 = 40_000
y = sur.copy()
mat = np.empty((B2, len(CUTS)))
for i in range(B2):
    rng.shuffle(y)
    mat[i] = chi2_all(y)

lev = (mat > crit).mean(axis=0)
print(f"개별 절단점에서 카이제곱 근사의 실제 수준 ({B2:,}회)")
print(f"  최소 {lev.min():.4f} ({CUTS[int(lev.argmin())]}세), "
      f"최대 {lev.max():.4f} ({CUTS[int(lev.argmax())]}세), "
      f"평균 {lev.mean():.4f}")
print(f"  k=2 가 쓰는 두 후보: {CUTS[0]}세 {lev[0]:.4f}, "
      f"{CUTS[-1]}세 {lev[-1]:.4f}")
print(f"  두 통계량의 상관 = {np.corrcoef(mat[:, 0], mat[:, -1])[0, 1]:+.4f}")
mx2 = mat[:, [0, len(CUTS) - 1]].max(axis=1)
print(f"  두 후보만 쓸 때 FWER = {(mx2 > crit).mean():.4f}")
print(f"  독립 예측 1-(1-{lev[0]:.4f})(1-{lev[-1]:.4f}) = "
      f"{1 - (1 - lev[0]) * (1 - lev[-1]):.4f}")

# 유효 후보 수
print(f"\n{'후보 수':>7}{'FWER':>9}{'1-0.95^k':>10}{'유효 후보 수':>13}")
for kk, fw in zip([1, 2, 3, 5, 10, 20, 64],
                  [0.0521, 0.1111, 0.1476, 0.2129, 0.3399, 0.4324, 0.5211]):
    print(f"{kk:>7}{fw:>9.4f}{1 - 0.95 ** kk:>10.4f}"
          f"{np.log(1 - fw) / np.log(0.95):>13.2f}")

개별 절단점에서 카이제곱 근사의 실제 수준 (40,000회)
  최소 0.0376 (6세), 최대 0.0673 (58세), 평균 0.0507
  k=2 가 쓰는 두 후보: 2세 0.0538, 65세 0.0576
  두 통계량의 상관 = -0.0016
  두 후보만 쓸 때 FWER = 0.1083
  독립 예측 1-(1-0.0538)(1-0.0576) = 0.1083

   후보 수     FWER  1-0.95^k      유효 후보 수
      1   0.0521    0.0500         1.04
      2   0.1111    0.0975         2.30
      3   0.1476    0.1426         3.11
      5   0.2129    0.2262         4.67
     10   0.3399    0.4013         8.10
     20   0.4324    0.6415        11.04
     64   0.5211    0.9625        14.35

import statsmodels.api as sm

y = d["Survived"].to_numpy()
models = {
    "절편만 (나이 무시)": np.ones((len(y), 1)),
    "21세 이분화": sm.add_constant((age < 21).astype(float)),
    "7세 이분화 (탐색으로 고름)": sm.add_constant((age < 7).astype(float)),
    "나이 연속 (선형)": sm.add_constant(age),
    "나이 + 나이 제곱": sm.add_constant(np.column_stack([age, age ** 2])),
    "10년 구간 더미": sm.add_constant(
        pd.get_dummies(np.digitize(age, [10, 20, 30, 40, 50, 60]),
                       drop_first=True).to_numpy(dtype=float)),
}
print(f"{'모형':>26s}{'모수':>6s}{'로그가능도':>12s}{'AIC':>10s}")
for lab, X in models.items():
    m = sm.Logit(y, X).fit(disp=0)
    print(f"{lab:>26s}{X.shape[1]:>6d}{m.llf:>12.3f}{m.aic:>10.3f}")

# 연속형 나이의 효과가 정말 있는지 직접 검정한다.
m = sm.Logit(y, sm.add_constant(age)).fit(disp=0)
print(f"\n로지스틱 회귀 (나이 연속)")
print(f"  나이 계수 {m.params[1]:+.6f}  (한 살 많아질수록 로그오즈 변화)")
print(f"  p = {m.pvalues[1]:.4f}")
print(f"  오즈비 (10살 차이) {np.exp(m.params[1] * 10):.4f}")

                        모형    모수       로그가능도       AIC
               절편만 (나이 무시)     1    -482.258   966.516
                   21세 이분화     2    -481.049   966.099
          7세 이분화 (탐색으로 고름)     2    -473.249   950.499
                나이 연속 (선형)     2    -480.114   964.228
                나이 + 나이 제곱     3    -478.906   963.812
                 10년 구간 더미     7    -473.913   961.827

로지스틱 회귀 (나이 연속)
  나이 계수 -0.010963  (한 살 많아질수록 로그오즈 변화)
  p = 0.0397
  오즈비 (10살 차이) 0.8962

print(f"모수 하나를 더 써서 AIC 가 낮아지는 조건: LR > 2, "
      f"곧 p < {chi2dist.sf(2, 1):.6f}")
print(f"모수 여섯을 더 쓰는 경우:              LR > 12, "
      f"곧 p < {chi2dist.sf(12, 6):.6f}")

fits = {lab: sm.Logit(y, X).fit(disp=0) for lab, X in models.items()}
print(f"\n{'비교':>34}{'LR':>8}{'df':>4}{'p':>9}{'dAIC':>8}")
pairs = [("절편만 (나이 무시)", "21세 이분화"),
         ("절편만 (나이 무시)", "나이 연속 (선형)"),
         ("나이 연속 (선형)", "나이 + 나이 제곱"),
         ("절편만 (나이 무시)", "10년 구간 더미")]
for a, b in pairs:
    ma, mb = fits[a], fits[b]
    lr = 2 * (mb.llf - ma.llf)
    dfd = models[b].shape[1] - models[a].shape[1]
    print(f"{a[:12] + ' -> ' + b[:14]:>34}{lr:>8.4f}{dfd:>4}"
          f"{chi2dist.sf(lr, dfd):>9.4f}{mb.aic - ma.aic:>8.4f}")
    assert abs((mb.aic - ma.aic) - (2 * dfd - lr)) < 1e-9
print("\n모든 쌍에서 dAIC = 2*df - LR 이 성립한다 (단언문 통과)")

모수 하나를 더 써서 AIC 가 낮아지는 조건: LR > 2, 곧 p < 0.157299
모수 여섯을 더 쓰는 경우:              LR > 12, 곧 p < 0.061969

                                비교      LR  df        p    dAIC
            절편만 (나이 무시) -> 21세 이분화  2.4173   1   0.1200 -0.4173
         절편만 (나이 무시) -> 나이 연속 (선형)  4.2876   1   0.0384 -2.2876
          나이 연속 (선형) -> 나이 + 나이 제곱  2.4164   1   0.1201 -0.4164
          절편만 (나이 무시) -> 10년 구간 더미 16.6893   6   0.0105 -4.6893

모든 쌍에서 dAIC = 2*df - LR 이 성립한다 (단언문 통과)

rng = np.random.default_rng(1)
B = 2_000

print("참 효과가 선형인 자료에서의 검정력 (n=300, 명목 0.05)")
print(f"{'beta':>6s}{'연속 그대로':>11s}{'중앙값 분할':>12s}{'손실':>8s}")
for beta in [0.2, 0.3, 0.4, 0.5]:
    a = b = 0
    for _ in range(B):
        x = rng.standard_normal(300)
        yy = (rng.random(300) < 1 / (1 + np.exp(-beta * x))).astype(int)
        a += sm.Logit(yy, sm.add_constant(x)).fit(disp=0).pvalues[1] < 0.05
        xd = (x > np.median(x)).astype(float)
        b += sm.Logit(yy, sm.add_constant(xd)).fit(disp=0).pvalues[1] < 0.05
    print(f"{beta:>6.1f}{a / B:>11.4f}{b / B:>12.4f}"
          f"{1 - (b / B) / (a / B):>8.1%}")

참 효과가 선형인 자료에서의 검정력 (n=300, 명목 0.05)
  beta     연속 그대로      중앙값 분할      손실
   0.2     0.3905      0.2840   27.3%
   0.3     0.7180      0.5280   26.5%
   0.4     0.9140      0.7610   16.7%
   0.5     0.9840      0.9120    7.3%

from scipy.stats import norm

rho = np.sqrt(2 / np.pi)
print(f"E[X 1{{X>0}}] = phi(0) = {1 / np.sqrt(2 * np.pi):.6f}")
print(f"corr(X, D)  = {1 / np.sqrt(2 * np.pi)} / (1 * 0.5) = {rho:.6f}"
      f"  (= sqrt(2/pi))")
print(f"설명하는 분산의 비 rho^2 = {rho**2:.6f}")
print(f"필요한 표본크기의 배율 1/rho^2 = pi/2 = {np.pi / 2:.6f}")

print(f"\n{'beta':>6}{'예측(연속)':>11}{'모의':>8}{'예측(분할)':>12}{'모의':>8}"
      f"{'예측 손실':>10}{'모의 손실':>10}")
sim_c = {0.2: 0.3905, 0.3: 0.7180, 0.4: 0.9140, 0.5: 0.9840}
sim_d = {0.2: 0.2840, 0.3: 0.5280, 0.4: 0.7610, 0.5: 0.9120}
z95 = 1.959964
for beta in (0.2, 0.3, 0.4, 0.5):
    lam = 0.5 * beta * np.sqrt(300)          # 0.5 = sqrt(p q) at p = 1/2
    pc = norm.sf(z95 - lam) + norm.cdf(-z95 - lam)
    pd_ = norm.sf(z95 - rho * lam) + norm.cdf(-z95 - rho * lam)
    print(f"{beta:>6.1f}{pc:>11.4f}{sim_c[beta]:>8.4f}{pd_:>12.4f}"
          f"{sim_d[beta]:>8.4f}{1 - pd_ / pc:>10.1%}"
          f"{1 - sim_d[beta] / sim_c[beta]:>10.1%}")

E[X 1{X>0}] = phi(0) = 0.398942
corr(X, D)  = 0.3989422804014327 / (1 * 0.5) = 0.797885  (= sqrt(2/pi))
설명하는 분산의 비 rho^2 = 0.636620
필요한 표본크기의 배율 1/rho^2 = pi/2 = 1.570796

  beta     예측(연속)      모의      예측(분할)      모의     예측 손실     모의 손실
   0.2     0.4100  0.3905      0.2821  0.2840     31.2%     27.3%
   0.3     0.7383  0.7180      0.5450  0.5280     26.2%     26.5%
   0.4     0.9337  0.9140      0.7893  0.7610     15.5%     16.7%
   0.5     0.9911  0.9840      0.9325  0.9120      5.9%      7.3%

# 보기 2 의 순열에서 절단점별 통계량을 모두 저장한다.
rng = np.random.default_rng(0)
B = 5_000
y = sur.copy()
mat = np.empty((B, len(CUTS)))
for i in range(B):
    rng.shuffle(y)
    mat[i] = chi2_all(y)

C = np.corrcoef(mat.T)
adj = [C[i, i + 1] for i in range(len(CUTS) - 1)]
print(f"이웃한 절단점 통계량의 상관: 평균 {np.mean(adj):.4f}, "
      f"최소 {min(adj):.4f}")
print(f"양 끝({CUTS[0]}세와 {CUTS[-1]}세)의 상관: {C[0, -1]:+.4f}")

print(f"\n독립이라고 가정하면 1 - 0.95^{len(CUTS)} = "
      f"{1 - 0.95 ** len(CUTS):.4f}")
print(f"실제 (보기 2)                      = 0.5194")

# 상관이 높으면 '서로 다른 검정'의 유효 개수가 줄어든다.
eff = np.log(1 - 0.5194) / np.log(0.95)
print(f"\n0.5194 를 내는 '유효 후보 수' = {eff:.1f}개")
print(f"  실제 후보 {len(CUTS)}개의 {eff / len(CUTS):.1%} 에 해당한다")

이웃한 절단점 통계량의 상관: 평균 0.8796, 최소 0.5626
양 끝(2세와 65세)의 상관: +0.0081

독립이라고 가정하면 1 - 0.95^64 = 0.9625
실제 (보기 2)                      = 0.5194

0.5194 를 내는 '유효 후보 수' = 14.3개
  실제 후보 64개의 22.3% 에 해당한다

from scipy.stats import chi2 as chi2dist

k = len(CUTS)
bonf = chi2dist.ppf(1 - 0.05 / k, 1)
perm = 9.0294                      # 보기 2 의 순열 95분위

print(f"후보 {k}개")
print(f"  보정 없음      임계값 {chi2dist.ppf(0.95, 1):.4f}")
print(f"  Bonferroni    임계값 {bonf:.4f}  (alpha/{k} = {0.05 / k:.6f})")
print(f"  순열검정       임계값 {perm:.4f}")

# 각 임계값에서 실제 1종 오류가 얼마인지 순열분포로 확인한다.
rng = np.random.default_rng(0)
B = 20_000
y = sur.copy()
null = np.empty(B)
for i in range(B):
    rng.shuffle(y)
    null[i] = chi2_all(y).max()

print(f"\n{'방법':>14s}{'임계값':>10s}{'실제 1종 오류':>14s}")
for lab, cv in [("보정 없음", chi2dist.ppf(0.95, 1)),
                ("Bonferroni", bonf),
                ("순열검정", np.quantile(null, 0.95))]:
    print(f"{lab:>14s}{cv:>10.4f}{np.mean(null > cv):>14.4f}")

후보 64개
  보정 없음      임계값 3.8415
  Bonferroni    임계값 11.2853  (alpha/64 = 0.000781)
  순열검정       임계값 9.0294

            방법       임계값      실제 1종 오류
         보정 없음    3.8415        0.5194
    Bonferroni   11.2853        0.0169
          순열검정    9.0294        0.0500

import statsmodels.api as sm
from scipy.stats import chi2 as chi2dist

y = d["Survived"].to_numpy()

# (1) 이차항이 필요한가 -- 우도비 검정
m1 = sm.Logit(y, sm.add_constant(age)).fit(disp=0)
m2 = sm.Logit(y, sm.add_constant(
    np.column_stack([age, age ** 2]))).fit(disp=0)
lr = 2 * (m2.llf - m1.llf)
print(f"(1) 선형 대 이차")
print(f"    우도비 = {lr:.4f}, df=1, p = {chi2dist.sf(lr, 1):.4f}")
print(f"    2차 계수 {m2.params[2]:+.8f}  (양수면 U자)")

# (2) 10년 구간 더미가 선형보다 나은가
D = pd.get_dummies(np.digitize(age, [10, 20, 30, 40, 50, 60]),
                   drop_first=True).to_numpy(dtype=float)
m3 = sm.Logit(y, sm.add_constant(D)).fit(disp=0)
m0 = sm.Logit(y, np.ones((len(y), 1))).fit(disp=0)
lr2 = 2 * (m3.llf - m0.llf)
print(f"\n(2) 절편만 대 10년 구간 더미")
print(f"    우도비 = {lr2:.4f}, df={D.shape[1]}, "
      f"p = {chi2dist.sf(lr2, D.shape[1]):.4f}")
lr3 = 2 * (m3.llf - m1.llf)
print(f"    선형 대 구간 더미: 우도비 = {lr3:.4f}, df={D.shape[1] - 1}, "
      f"p = {chi2dist.sf(lr3, D.shape[1] - 1):.4f}")

# (3) 10세 미만을 따로 떼면 나머지는 평평한가
older = age >= 10
m4 = sm.Logit(y[older], sm.add_constant(age[older])).fit(disp=0)
print(f"\n(3) 10세 이상만 (n={older.sum()})")
print(f"    나이 계수 {m4.params[1]:+.6f}, p = {m4.pvalues[1]:.4f}")
m5 = sm.Logit(y[~older], sm.add_constant(age[~older])).fit(disp=0)
print(f"    10세 미만만 (n={(~older).sum()})")
print(f"    나이 계수 {m5.params[1]:+.6f}, p = {m5.pvalues[1]:.4f}")

(1) 선형 대 이차
    우도비 = 2.4164, df=1, p = 0.1201
    2차 계수 +0.00039646  (양수면 U자)

(2) 절편만 대 10년 구간 더미
    우도비 = 16.6893, df=6, p = 0.0105
    선형 대 구간 더미: 우도비 = 12.4017, df=5, p = 0.0297

(3) 10세 이상만 (n=652)
    나이 계수 -0.000421, p = 0.9470
    10세 미만만 (n=62)
    나이 계수 -0.201384, p = 0.0376

전체 714명:      나이 계수 -0.0110,  p = 0.0397  (유의)
10세 이상 652명: 나이 계수 -0.0004,  p = 0.9470  (효과 없음)
10세 미만  62명: 나이 계수 -0.2014,  p = 0.0376  (강한 효과)

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(2)
B = 2_000
BETA = 0.3

def power(n_, dichotomize):
    hit = 0
    for _ in range(B):
        x = rng.standard_normal(n_)
        yy = (rng.random(n_) < 1 / (1 + np.exp(-BETA * x))).astype(int)
        if yy.sum() in (0, n_):
            continue
        v = (x > np.median(x)).astype(float) if dichotomize else x
        hit += sm.Logit(yy, sm.add_constant(v)).fit(disp=0).pvalues[1] < 0.05
    return hit / B

base = power(300, False)
print(f"기준: 연속형, n=300 -> 검정력 {base:.4f}\n")
print(f"{'n':>6s}{'이분화 검정력':>13s}{'기준 대비':>10s}")
for n_ in [300, 400, 450, 500, 600]:
    p_ = power(n_, True)
    print(f"{n_:>6d}{p_:>13.4f}{p_ / base:>10.4f}")

기준: 연속형, n=300 -> 검정력 0.7130

     n      이분화 검정력     기준 대비
   300       0.5375    0.7539
   400       0.6490    0.9102
   450       0.7265    1.0189
   500       0.7585    1.0638
   600       0.8115    1.1381

import numpy as np, pandas as pd
from scipy import stats

URL = ("https://raw.githubusercontent.com/datasciencedojo/"
       "datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv")
d = pd.read_csv(URL).dropna(subset=["Age"])
age = d["Age"].to_numpy()
sur = d["Survived"].to_numpy().astype(float)
CUTS = np.array([c for c in range(1, 80)
                 if (age < c).sum() >= 10 and (age >= c).sum() >= 10])

def chi_all(a, s, cuts, minn=10):
    """모든 절단점의 chi2 를 2x2 닫힌 공식으로 한 번에 계산한다."""
    M = (a[:, None] < cuts[None, :])
    n = len(a); n1 = M.sum(0); n0 = n - n1
    s1 = (M * s[:, None]).sum(0); s0 = s.sum() - s1
    A, B, C, D = s1, n1 - s1, s0, n0 - s0
    den = (A + B) * (C + D) * (A + C) * (B + D)
    with np.errstate(divide='ignore', invalid='ignore'):
        chi = np.where(den > 0, n * (A * D - B * C) ** 2 / den, 0.0)
    return np.where((n1 >= minn) & (n0 >= minn), chi, np.nan)

pv = stats.chi2.sf(chi_all(age, sur, CUTS), 1)
print(f"관측: p<0.05 {int((pv<0.05).sum())}개, "
      f"p<0.01 {int((pv<0.01).sum())}개, 중앙값 {np.median(pv):.4f}")

rng = np.random.default_rng(0)
c5, c1, med = [], [], []
for _ in range(2000):
    p2 = stats.chi2.sf(chi_all(age, rng.permutation(sur), CUTS), 1)
    c5.append((p2 < 0.05).sum()); c1.append((p2 < 0.01).sum())
    med.append(np.median(p2))
print(f"귀무: p<0.05 {np.mean(c5):.2f}개, "
      f"p<0.01 {np.mean(c1):.2f}개, 중앙값 {np.mean(med):.4f}")
print(f"균등분포 기대: {0.05*64:.2f}개, {0.01*64:.2f}개, 0.5")

관측: p<0.05 21개, p<0.01 17개, 중앙값 0.2306
귀무: p<0.05 3.13개, p<0.01 0.54개, 중앙값 0.4895
균등분포 기대: 3.20개, 0.64개, 0.5

import numpy as np, pandas as pd

URL = ("https://raw.githubusercontent.com/datasciencedojo/"
       "datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv")
d = pd.read_csv(URL).dropna(subset=["Age"])
age = d["Age"].to_numpy()
sur = d["Survived"].to_numpy().astype(float)
CUTS = np.array([c for c in range(1, 80)
                 if (age < c).sum() >= 10 and (age >= c).sum() >= 10])

def scan(a, s, cuts, minn=5):
    M = (a[:, None] < cuts[None, :])
    n = len(a); n1 = M.sum(0); n0 = n - n1
    s1 = (M * s[:, None]).sum(0); s0 = s.sum() - s1
    A, B, C, D = s1, n1 - s1, s0, n0 - s0
    den = (A + B) * (C + D) * (A + C) * (B + D)
    with np.errstate(divide='ignore', invalid='ignore'):
        chi = np.where(den > 0, n * (A * D - B * C) ** 2 / den, 0.0)
    chi = np.where((n1 >= minn) & (n0 >= minn), chi, -1.0)
    return chi, n1, n0, s1, s0

rng = np.random.default_rng(0)
naive, split = [], []
for _ in range(2000):
    i = rng.permutation(len(age)); h = len(age) // 2
    A_, B_ = i[:h], i[h:]
    chi, n1, n0, s1, s0 = scan(age[A_], sur[A_], CUTS)
    if chi.max() < 0:
        continue
    j = int(np.argmax(chi))
    naive.append(s1[j] / n1[j] - s0[j] / n0[j])       # 고른 곳에서 재기
    mB = age[B_] < CUTS[j]
    if mB.sum() >= 5 and (~mB).sum() >= 5:
        split.append(sur[B_][mB].mean() - sur[B_][~mB].mean())   # 다른 쪽에서 재기

print(f"  같은 절반에서 고르고 재기      : {np.mean(naive):+.4f}")
print(f"  절반에서 고르고 나머지에서 재기: {np.mean(split):+.4f}")
print(f"  과대추정 폭                    : {np.mean(naive)-np.mean(split):+.4f}")

  같은 절반에서 고르고 재기      : +0.3506
  절반에서 고르고 나머지에서 재기: +0.2520
  과대추정 폭                    : +0.0986

import numpy as np, pandas as pd
from scipy import stats
# (연습문제 5의 chi_all, CUTS, age, sur 를 그대로 쓴다)

pv = stats.chi2.sf(chi_all(age, sur, CUTS), 1)
m, a = len(pv), 0.05

bon = a / m
o = np.argsort(pv)
passed = pv[o] <= np.arange(1, m + 1) / m * a
k = np.max(np.where(passed)[0]) + 1 if passed.any() else 0

print(f"검정 {m}개, alpha={a}")
print(f"  보정 없이 p<0.05 : {int((pv<0.05).sum())}개")
print(f"  Bonferroni (p < {bon:.6f}) : {int((pv<bon).sum())}개")
print(f"  BH (FDR<=0.05)   : {k}개,  최대 기각 p = {pv[o][k-1]:.6f}")
print(f"\n  BH 로 기각된 절단점: {sorted(np.array(CUTS)[o][:k])}")

검정 64개, alpha=0.05
  보정 없이 p<0.05 : 21개
  Bonferroni (p < 0.000781) : 8개
  BH (FDR<=0.05)   : 17개,  최대 기각 p = 0.009148

  BH 로 기각된 절단점: [2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]

import numpy as np

M = np.column_stack([(age < c).astype(float) for c in CUTS])
R = np.corrcoef(M, rowvar=False)
ev = np.linalg.eigvalsh(R)
m = len(CUTS)
Meff = 1 + (m - 1) * (1 - np.var(ev, ddof=1) / m)

print(f"  명목 검정 수 m        = {m}")
print(f"  상관구조 기반 M_eff   = {Meff:.2f}")
print(f"  오류율에서 역산 m_eff = {np.log(1-0.52)/np.log(1-0.05):.2f}")

  명목 검정 수 m        = 64
  상관구조 기반 M_eff   = 49.98
  오류율에서 역산 m_eff = 14.31

import numpy as np, pandas as pd
from scipy import stats

URL = ("https://raw.githubusercontent.com/datasciencedojo/"
       "datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv")
dd = pd.read_csv(URL).dropna(subset=["Age"]).copy()
dd["FamSize"] = dd.SibSp + dd.Parch
sur2 = dd.Survived.to_numpy().astype(float)
VARS = {k: dd[k].to_numpy().astype(float)
        for k in ("Age", "Fare", "SibSp", "Parch", "FamSize")}

def cand(v):
    return np.array([c for c in np.unique(v)
                     if (v < c).sum() >= 10 and (v >= c).sum() >= 10])

tot = 0
print(f"{'변수':>9}{'절단점 수':>10}{'유의(p<.05)':>13}{'최소 p':>12}")
for k, v in VARS.items():
    cs = cand(v)
    p = stats.chi2.sf(chi_all(v, sur2, cs), 1)
    tot += len(cs)
    print(f"{k:>9}{len(cs):>10}{int(np.nansum(p<0.05)):>13}{np.nanmin(p):>12.2e}")
print(f"\n전체 후보 수 {tot}개")

rng = np.random.default_rng(0)
hits = 0
for _ in range(1000):
    sp = rng.permutation(sur2)
    if any(np.nanmin(stats.chi2.sf(chi_all(v, sp, cand(v)), 1)) < 0.05
           for v in VARS.values()):
        hits += 1
print(f"귀무에서 '어느 변수든 유의한 절단점 하나' 발견 확률: {hits/1000:.3f}")

       변수     절단점 수    유의(p<.05)        최소 p
      Age        76           22    1.92e-05
     Fare       213          213    7.90e-17
    SibSp         4            3    6.20e-03
    Parch         4            3    1.11e-05
  FamSize         6            4    1.60e-07

전체 후보 수 303개
귀무에서 '어느 변수든 유의한 절단점 하나' 발견 확률: 0.874

1. 자를 이유가 자료 밖에 있는가?
   - 법정 연령, 의학적 기준, 정책 경계 -> 자른다 (절단점이 사전에 정해짐)
   - "해 보니 이게 제일 잘 나와서"   -> 자르지 않는다
2. 자르지 않고 답할 수 있는가?
   - 로지스틱 회귀, 스플라인, 평활 -> 대개 가능하다
3. 그래도 잘라야 한다면 몇 개를 시도할 것인가?
   - 시도 횟수를 미리 적어 둔다

나이와 생존의 관계 (n = 714, 나이 결측 177명 제외)

  탐색: 2~65세의 정수 절단점 64개에서 카이제곱 검정
  최대: 7세에서 chi2 = 18.27
  순열 임계값(10,000회, alpha=0.05): 9.03  -> 유의
  단, 절단점 7세는 자료에서 고른 것이므로
    이 위치 자체를 결론으로 삼지 않는다.
  재표본 검사에서 최적 절단점은 2~64세로 흔들린다.

  효과크기: 표본을 반으로 나눠 추정한 생존율 차이 0.25
    (같은 자료에서 고르고 재면 0.35 로 과대추정)

  연속형 분석: 나이를 그대로 쓴 로지스틱 회귀에서
    나이 효과는 비단조이며, 이분화는 이 구조를 지운다.