독립성 검정 템플릿 함수¶

import numpy as np
from scipy import stats

def chi2_independence(observed: np.ndarray, correction: bool = False):
    """카이제곱 독립성 검정.

    correction의 기본값을 **False**로 두었다는 점에 주의하라.
    scipy의 기본값은 True이며, 2x2 표에 Yates 연속성 보정을 자동으로 적용한다.
    모르고 쓰면 작은 2x2 표에서 통계량이 조용히 줄어들어 결론이 달라질 수 있다.
    돌려주는 값은 (통계량, p-값, 자유도, 기대도수).
    """
    return stats.chi2_contingency(observed, correction=correction)


# 기본값의 차이를 눈으로 확인한다.
tab = np.array([[10, 5], [3, 12]], dtype=float)
print("correction=False:", round(chi2_independence(tab, False)[0], 4))
print("correction=True :", round(chi2_independence(tab, True)[0], 4))

# (1)(2) 닫힌 꼴과 맞는지
(a, b), (c, d) = tab
n = tab.sum()
D = (a + b) * (c + d) * (a + c) * (b + d)
delta = abs(a * d - b * c)
print(f"\n닫힌 꼴  n(ad-bc)^2/D          = {n * delta**2 / D:.6f}")
print(f"예이츠   n(|ad-bc|-n/2)^2/D    = {n * (delta - n / 2)**2 / D:.6f}")
print(f"차       (n^2/D)(|ad-bc|-n/4)  = {n**2 / D * (delta - n / 4):.6f}")
phi = (a * d - b * c) / np.sqrt(D)
print(f"phi = (ad-bc)/sqrt(D) = {phi:.6f}   n*phi^2 = {n * phi**2:.6f}")

# (3) 자유도 1 의 p-값은 표준정규로 적힌다
for name, x in [("보정 없음", chi2_independence(tab, False)[0]),
                ("예이츠   ", chi2_independence(tab, True)[0])]:
    print(f"{name}  chi2 = {x:.6f}  p(chi2) = {stats.chi2(1).sf(x):.7f}"
          f"  2*Phi(-sqrt) = {2 * stats.norm.sf(np.sqrt(x)):.7f}")
print(f"1% 임계값 chi2_(1, 0.01) = {stats.chi2.ppf(0.99, 1):.6f}")

correction=False: 6.6516
correction=True : 4.8869

닫힌 꼴  n(ad-bc)^2/D          = 6.651584
예이츠   n(|ad-bc|-n/2)^2/D    = 4.886878
차       (n^2/D)(|ad-bc|-n/4)  = 1.764706
phi = (ad-bc)/sqrt(D) = 0.470871   n*phi^2 = 6.651584
보정 없음  chi2 = 6.651584  p(chi2) = 0.0099068  2*Phi(-sqrt) = 0.0099068
예이츠     chi2 = 4.886878  p(chi2) = 0.0270616  2*Phi(-sqrt) = 0.0270616
1% 임계값 chi2_(1, 0.01) = 6.634897

observed = np.array([[30, 20, 10],
                     [12, 25, 18]], dtype=float)

# 2x2 표가 아니면 연속성 보정은 뜻이 없다. 그래서 correction=False 로 둔다.
chi2, p, df, exp = chi2_independence(observed, correction=False)
print(f"chi2 = {chi2:.3f}, p = {p:.4f}, df = {df}")
print("expected:\n", exp)

# (1) 2x3 에서는 correction 인자가 무시된다
print(f"\ncorrection=True 로 주면 {chi2_independence(observed, True)[0]:.6f}"
      f"   (False 와 같다)")

# (2) 2 x c 의 닫힌 꼴
n = observed.sum()
R, C = observed.sum(1), observed.sum(0)
e = observed[0] - R[0] * C / n
closed = n**2 / (R[0] * R[1]) * np.sum(e**2 / C)
print(f"e_j = {np.round(e, 5)}   합 {e.sum():.1e}")
print(f"닫힌 꼴 {closed:.5f}   정의대로 {chi2:.5f}")
print(f"p = exp(-chi2/2) = {np.exp(-chi2 / 2):.7f}   sf = {p:.7f}")

# (3) 조정 잔차
adj = (observed - exp) / np.sqrt(exp * np.outer(1 - R / n, 1 - C / n))
print(f"\n조정 잔차\n{np.round(adj, 4)}")
print(f"본페로니 임계값 {stats.norm.ppf(1 - 0.025 / 6):.4f}")
print(f"칸별 기여 (%)\n{np.round((observed - exp)**2 / exp / chi2 * 100, 1)}")

chi2 = 10.358, p = 0.0056, df = 2
expected:
 [[21.91304348 23.47826087 14.60869565]
 [20.08695652 21.52173913 13.39130435]]

correction=True 로 주면 10.357744   (False 와 같다)
e_j = [ 8.08696 -3.47826 -4.6087 ]   합 -1.8e-15
닫힌 꼴 10.35774   정의대로 10.35774
p = exp(-chi2/2) = 0.0056344   sf = 0.0056344

조정 잔차
[[ 3.1354 -1.3304 -2.0046]
 [-3.1354  1.3304  2.0046]]
본페로니 임계값 2.6383
칸별 기여 (%)
[[28.8  5.  14. ]
 [31.4  5.4 15.3]]

chi2, p, df, exp = chi2_independence(
    np.array([[50, 50], [50, 50]], dtype=float)
)
print(f"chi2 = {chi2}, p = {p}, df = {df}")

chi2 = 0.0, p = 1.0, df = 1

observed = np.array([[30, 20, 10],
                     [12, 25, 18]], dtype=float)
_, _, _, expected = chi2_independence(observed)
# 표준화 잔차. 제곱해서 모두 더하면 카이제곱 통계량이 된다.
residuals = (observed - expected) / np.sqrt(expected)
print(residuals)
print("제곱합 =", round((residuals**2).sum(), 3))

[[ 1.72756246 -0.71784254 -1.20579175]
 [-1.80438014  0.74976208  1.25940841]]
제곱합 = 10.358

table = np.array([[45, 30, 25],
                  [35, 40, 25]], dtype=float)
chi2, p, df, exp = chi2_independence(table, correction=False)
print(f"chi2 = {chi2:.4f}, p = {p:.4f}, df = {df}")

chi2 = 2.6786, p = 0.2620, df = 2

import numpy as np
from scipy import stats

tab = np.array([[10, 5], [3, 12]], float)
chi2, p, df, exp = stats.chi2_contingency(tab, correction=False)
n = tab.sum()
rp, cp = tab.sum(1) / n, tab.sum(0) / n

resid = (tab - exp) / np.sqrt(exp)
adj = (tab - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))
contrib = (tab - exp)**2 / exp

print(f"χ² = {chi2:.4f},  df = {df},  p = {p:.4f}")
print(f"기대도수\n{np.round(exp, 4)}")
print(f"\n표준화 잔차\n{np.round(resid, 4)}")
print(f"조정 잔차\n{np.round(adj, 4)}")
print(f"χ² 기여율(%)\n{np.round(contrib / chi2 * 100, 2)}")
print(f"\n검산  ΣR² = {np.sum(resid**2):.4f}   (= χ²)")

χ² = 6.6516,  df = 1,  p = 0.0099
기대도수
[[6.5 8.5]
 [6.5 8.5]]

표준화 잔차
[[ 1.3728 -1.2005]
 [-1.3728  1.2005]]
조정 잔차
[[ 2.5791 -2.5791]
 [-2.5791  2.5791]]
χ² 기여율(%)
[[28.33 21.67]
 [28.33 21.67]]

검산  ΣR² = 6.6516   (= χ²)

a, b, c, d = tab.ravel()
p1, p2 = a / (a + b), c / (c + d)
print(f"행별 비율: {p1:.4f} 대 {p2:.4f}   차이 {p1 - p2:+.4f}")
print(f"오즈비 {(a * d) / (b * c):.4f}")
odds, p_f = stats.fisher_exact(tab)
print(f"피셔 정확검정 p = {p_f:.4f}   (카이제곱 p = {p:.4f})")

행별 비율: 0.6667 대 0.2000   차이 +0.4667
오즈비 8.0000
피셔 정확검정 p = 0.0253   (카이제곱 p = 0.0099)

import numpy as np
from scipy import stats

def chi2_report(observed, alpha=0.05, correction=False, labels=None):
    """독립성 검정과 진단 정보를 한 번에 제시한다.

    correction 의 기본값은 False 다 (scipy 와 반대).
    """
    obs = np.asarray(observed, float)
    if obs.ndim != 2:
        raise ValueError("2차원 분할표가 필요하다")
    if obs.sum(0).min() == 0 or obs.sum(1).min() == 0:
        raise ValueError("합이 0 인 행 또는 열이 있다")

    r, c = obs.shape
    n = obs.sum()
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=correction)
    v = np.sqrt(chi2 / (n * (min(r, c) - 1)))
    v_null = np.sqrt(df / (n * (min(r, c) - 1)))
    rp, cp = obs.sum(1) / n, obs.sum(0) / n
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))
    z_bonf = stats.norm.ppf(1 - alpha / 2 / (r * c))

    print(f"{r}×{c} 표,  n = {n:.0f}")
    print(f"  χ² = {chi2:.4f},  df = {df},  p = {p:.6f}"
          f"   → {'기각' if p < alpha else '기각 못 함'}")
    print(f"  크라메르 V = {v:.4f}   (독립일 때 기댓값 {v_null:.4f}, "
          f"비 {v / v_null:.2f}배)")
    print(f"  E_min = {exp.min():.3f}"
          + ("   ⚠ 5 미만 — 정확검정을 고려하라" if exp.min() < 5 else ""))
    if r == 2 and c == 2:
        a, b, cc, d = obs.ravel()
        if b * cc > 0:
            print(f"  오즈비 = {(a * d) / (b * cc):.4f}"
                  f"   (2×2 이므로 잔차 분석은 생략)")
    elif p < alpha:
        flagged = np.abs(adj) > z_bonf
        print(f"  조정 잔차 (본페로니 임계값 {z_bonf:.3f}, "
              f"유의한 칸 {int(flagged.sum())}개)")
        print(np.round(adj, 3))
    return {"chi2": chi2, "df": df, "p": p, "V": v,
            "expected": exp, "adj_resid": adj}

cases = {
    "처치 × 중증도": [[20, 15, 5], [10, 20, 10]],
    "성별 × 손잡이": [[934, 1070], [113, 92], [20, 8]],
    "작은 2×2": [[10, 5], [3, 12]],
}
for name, T in cases.items():
    print(f"[{name}]")
    chi2_report(T)
    print()

[처치 × 중증도]
2×3 표,  n = 80
  χ² = 5.7143,  df = 2,  p = 0.057433   → 기각 못 함
  크라메르 V = 0.2673   (독립일 때 기댓값 0.1581, 비 1.69배)
  E_min = 7.500

[성별 × 손잡이]
3×2 표,  n = 2237
  χ² = 11.8061,  df = 2,  p = 0.002731   → 기각
  크라메르 V = 0.0726   (독립일 때 기댓값 0.0299, 비 2.43배)
  E_min = 13.355
  조정 잔차 (본페로니 임계값 2.638, 유의한 칸 2개)
[[-3.03   3.03 ]
 [ 2.233 -2.233]
 [ 2.53  -2.53 ]]

[작은 2×2]
2×2 표,  n = 30
  χ² = 6.6516,  df = 1,  p = 0.009907   → 기각
  크라메르 V = 0.4709   (독립일 때 기댓값 0.1826, 비 2.58배)
  E_min = 6.500
  오즈비 = 8.0000   (2×2 이므로 잔차 분석은 생략)

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

def batch_chi2(tables, alpha=0.05, method="holm"):
    """여러 분할표를 한 번에 검정하고 다중비교를 보정한다."""
    rows = []
    for name, T in tables.items():
        T = np.asarray(T, float)
        chi2, p, df, exp = stats.chi2_contingency(T, correction=False)
        n = T.sum()
        v = np.sqrt(chi2 / (n * (min(T.shape) - 1)))
        rows.append([name, n, chi2, df, p, v, exp.min()])

    pvals = np.array([r[4] for r in rows])
    rej, adj, _, _ = multipletests(pvals, alpha=alpha, method=method)

    print(f"{'표':>14s} {'n':>6s} {'χ²':>9s} {'df':>3s} {'p':>10s} "
          f"{'보정 p':>10s} {'V':>7s} {'E_min':>7s}")
    for row, a, r in zip(rows, adj, rej):
        name, n, chi2, df, p, v, emin = row
        flag = "*" if r else " "
        warn = " ⚠" if emin < 5 else ""
        print(f"{name:>14s} {n:6.0f} {chi2:9.4f} {df:3d} {p:10.6f} "
              f"{a:10.6f}{flag} {v:7.4f} {emin:7.2f}{warn}")
    return rows, adj, rej

tables = {
    "결제 × 요일": [[30, 50, 20], [40, 60, 30]],
    "성별 × 손잡이": [[934, 1070], [113, 92], [20, 8]],
    "처치 × 중증도": [[20, 15, 5], [10, 20, 10]],
    "지역 × 선호": [[30, 20], [25, 25], [20, 30], [35, 15], [15, 35]],
}
_ = batch_chi2(tables)

             표      n        χ²  df          p       보정 p       V   E_min
       결제 × 요일    230    0.4320   2   0.805748   0.805748   0.0433   21.74
      성별 × 손잡이   2237   11.8061   2   0.002731   0.008193*  0.0726   13.36
      처치 × 중증도     80    5.7143   2   0.057433   0.114865   0.2673    7.50
       지역 × 선호    250   20.0000   4   0.000499   0.001998*  0.2828   25.00

import numpy as np
from scipy import stats

def text_mosaic(observed, row_labels=None, col_labels=None, width=40):
    """행별 비율을 막대로 그리고 조정 잔차를 기호로 표시한다."""
    obs = np.asarray(observed, float)
    r, c = obs.shape
    n = obs.sum()
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    rp, cp = obs.sum(1) / n, obs.sum(0) / n
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))
    z_bonf = stats.norm.ppf(1 - 0.025 / (r * c))

    row_labels = row_labels or [f"행{i + 1}" for i in range(r)]
    col_labels = col_labels or [f"열{j + 1}" for j in range(c)]
    marks = "".join(f"{lab:>10s}" for lab in col_labels)
    print(f"{'':>8s}{marks}{'  n':>6s}")

    blocks = "░▒▓█"
    for i in range(r):
        props = obs[i] / obs[i].sum()
        cells = ""
        for j in range(c):
            if abs(adj[i, j]) > z_bonf:
                sym = "+" if adj[i, j] > 0 else "-"
            elif abs(adj[i, j]) > 1.96:
                sym = "." 
            else:
                sym = " "
            bar = blocks[min(3, int(props[j] * 4))]
            cells += f"{bar * 3}{props[j] * 100:4.0f}%{sym:>2s}"
        print(f"{row_labels[i]:>8s}{cells}{obs[i].sum():6.0f}")

    print(f"\nχ² = {chi2:.4f},  df = {df},  p = {p:.4f},  "
          f"V = {np.sqrt(chi2 / (n * (min(r, c) - 1))):.4f}")
    print(f"기호:  + 유의하게 많음   - 유의하게 적음   "
          f". |조정잔차|>1.96 (보정 전)")

text_mosaic([[934, 1070], [113, 92], [20, 8]],
            ["오른손", "왼손", "양손"], ["남", "여"])

                 남         여     n
     오른손▒▒▒  47% -▓▓▓  53% +  2004
      왼손▓▓▓  55% .▒▒▒  45% .   205
      양손▓▓▓  71% .▒▒▒  29% .    28

χ² = 11.8061,  df = 2,  p = 0.0027,  V = 0.0726
기호:  + 유의하게 많음   - 유의하게 적음   . |조정잔차|>1.96 (보정 전)

def bad_style(obs):
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    return chi2, p, df, exp                    # 위치로만 구분되는 튜플

def good_style(obs):
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    n, q = obs.sum(), min(obs.shape) - 1
    return {"chi2": chi2, "df": df, "p": p,
            "V": np.sqrt(chi2 / (n * q)),
            "expected": exp, "E_min": exp.min()}

import numpy as np
from scipy import stats
out = good_style(np.array([[10, 5], [3, 12]], float))
print({k: (round(v, 4) if np.isscalar(v) else "…") for k, v in out.items()})

{'chi2': 6.6516, 'df': 1, 'p': 0.0099, 'V': 0.4709, 'expected': '…', 'E_min': 6.5}