동질성 잔차 열지도¶

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

observed = np.array([
    [25, 30, 20, 25],
    [18, 22, 35, 25],
    [30, 25, 15, 30],
], dtype=float)

row_tot = observed.sum(axis=1, keepdims=True)
col_tot = observed.sum(axis=0, keepdims=True)
tot = observed.sum()
expected = (row_tot @ col_tot) / tot

# Pearson 표준화 잔차. 제곱해서 모두 더하면 카이제곱 통계량이 된다.
# 분모가 sqrt(E)인 것은 H0 아래에서 각 칸 도수의 분산이 근사적으로 E이기 때문이다.
resid = (observed - expected) / np.sqrt(expected)

# 칸마다 z-검정을 하는 셈이라 다중검정 문제가 생긴다. 그래서 아래에서 보정한다.
z = resid.ravel()
pvals = 2 * (1 - stats.norm.cdf(np.abs(z)))
reject, pvals_bonf, _, _ = multipletests(pvals, method="bonferroni")
pvals_bonf = pvals_bonf.reshape(observed.shape)
reject = reject.reshape(observed.shape)

print("Standardized residuals:")
print(resid)
print()
print("Bonferroni-adjusted per-cell p-values:")
print(pvals_bonf)

# (2) 전체 검정
chi2 = (resid ** 2).sum()
df = (observed.shape[0] - 1) * (observed.shape[1] - 1)
print(f"\nchi2 = {chi2:.4f}, df = {df}, p = {stats.chi2(df).sf(chi2):.5f}")

# (3) 표준화 잔차의 분산은 1 이 아니다
var_factor = (1 - row_tot / tot) * (1 - col_tot / tot)
print(f"Var(R_ij) 범위  {var_factor.min():.4f} ~ {var_factor.max():.4f}")
print(f"Var(R_ij) 의 합 {var_factor.sum():.4f}   (r-1)(c-1) = {df}")

# (4) 조정 잔차
adj = resid / np.sqrt(var_factor)
adj_bonf = np.minimum(1.0, 12 * 2 * stats.norm.sf(np.abs(adj)))
print(f"\n조정 잔차\n{np.round(adj, 4)}")
print(f"최대 |조정 잔차| {np.abs(adj).max():.4f}  "
      f"본페로니 임계값 {stats.norm.ppf(1 - 0.025 / 12):.4f}")
print(f"그 칸의 보정 p-값  표준화 {pvals_bonf.min():.4f}  ->  조정 {adj_bonf.min():.4f}")

Standardized residuals:
[[ 0.13514748  0.8553372  -0.69006556 -0.32274861]
 [-1.28390102 -0.72374686  2.41522946 -0.32274861]
 [ 1.14875354 -0.13159034 -1.7251639   0.64549722]]

Bonferroni-adjusted per-cell p-values:
[[1.        1.        1.        1.       ]
 [1.        1.        0.1887036 1.       ]
 [1.        1.        1.        1.       ]]

chi2 = 14.1697, df = 6, p = 0.02780
Var(R_ij) 범위  0.4889 ~ 0.5111
Var(R_ij) 의 합 6.0000   (r-1)(c-1) = 6

조정 잔차
[[ 0.1903  1.215  -0.9652 -0.4616]
 [-1.8077 -1.0281  3.3783 -0.4616]
 [ 1.6174 -0.1869 -2.4131  0.9232]]
최대 |조정 잔차| 3.3783  본페로니 임계값 2.8653
그 칸의 보정 p-값  표준화 0.1887  ->  조정 0.0088

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 4))
im = ax.imshow(resid, aspect="auto")
ax.set_title("Standardized residuals heatmap")
ax.set_xlabel("Category")
ax.set_ylabel("Population")
plt.colorbar(im, ax=ax, shrink=0.8)

# 본페로니 보정 뒤에도 유의한 칸에 표시를 남긴다
for i in range(observed.shape[0]):
    for j in range(observed.shape[1]):
        # 보정 후 유의한 칸에만 별표를 붙인다. 이 보기에서는 하나도 없다.
        mark = "*" if reject[i, j] else ""
        ax.text(j, i, f"{resid[i, j]:.2f}{mark}",
                ha="center", va="center", fontsize=10)

plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

rng = np.random.default_rng(3690)
M, n = 20_000, 300
p_row = np.array([1 / 3, 1 / 3, 1 / 3])
p_col = np.array([0.25, 0.26, 0.23, 0.26])
P = np.outer(p_row, p_col)                   # H0: 독립이 참

std_r, adj_r = [], []
for _ in range(M):
    obs = rng.multinomial(n, P.ravel()).reshape(3, 4).astype(float)
    if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
        continue
    exp = np.outer(obs.sum(1), obs.sum(0)) / n
    rp, cp = obs.sum(1) / n, obs.sum(0) / n
    std_r.append(((obs - exp) / np.sqrt(exp))[1, 2])
    adj_r.append(((obs - exp)
                  / np.sqrt(exp * np.outer(1 - rp, 1 - cp)))[1, 2])

std_r, adj_r = np.array(std_r), np.array(adj_r)
print(f"표준화 잔차 R    평균 {std_r.mean():+.4f}  표준편차 "
      f"{std_r.std(ddof=1):.4f}  |·|>1.96 비율 {np.mean(np.abs(std_r) > 1.96):.4f}")
print(f"조정 잔차   R^adj 평균 {adj_r.mean():+.4f}  표준편차 "
      f"{adj_r.std(ddof=1):.4f}  |·|>1.96 비율 {np.mean(np.abs(adj_r) > 1.96):.4f}")
print(f"\n이론값: R 의 표준편차 = √((1-R_i/n)(1-C_j/n)) = "
      f"√((1-1/3)(1-0.23)) = {np.sqrt((1 - 1 / 3) * (1 - 0.23)):.4f}")

표준화 잔차 R    평균 +0.0033  표준편차 0.7161  |·|>1.96 비율 0.0062
조정 잔차   R^adj 평균 +0.0047  표준편차 0.9991  |·|>1.96 비율 0.0494

이론값: R 의 표준편차 = √((1-R_i/n)(1-C_j/n)) = √((1-1/3)(1-0.23)) = 0.7165

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

obs = np.array([[25, 30, 20, 25],
                [18, 22, 35, 25],
                [30, 25, 15, 30]], dtype=float)
chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
n = obs.sum()
rp, cp = obs.sum(1) / n, obs.sum(0) / n

std_r = (obs - exp) / np.sqrt(exp)
adj_r = (obs - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))

print(f"전체 검정:  χ² = {chi2:.4f},  df = {df},  p = {p:.6f}")
print(f"검산  ΣR² = {np.sum(std_r**2):.4f}\n")
print("표준화 잔차 R\n", np.round(std_r, 4))
print("\n조정 잔차 R^adj\n", np.round(adj_r, 4))
print(f"\n|R| 최대 {np.abs(std_r).max():.4f},  "
      f"|R^adj| 최대 {np.abs(adj_r).max():.4f}\n")

for arr, name in [(std_r, "표준화"), (adj_r, "조정  ")]:
    pv = 2 * stats.norm.sf(np.abs(arr)).ravel()
    for method, label in [("bonferroni", "본페로니"),
                          ("holm", "홀름   "),
                          ("fdr_bh", "BH     ")]:
        rej, adj_p, _, _ = multipletests(pv, alpha=0.05, method=method)
        print(f"  {name} 잔차 + {label}: 최소 조정 p = {adj_p.min():.4f}, "
              f"유의한 칸 {int(rej.sum())}개")

전체 검정:  χ² = 14.1697,  df = 6,  p = 0.027796
검산  ΣR² = 14.1697

표준화 잔차 R
 [[ 0.1351  0.8553 -0.6901 -0.3227]
 [-1.2839 -0.7237  2.4152 -0.3227]
 [ 1.1488 -0.1316 -1.7252  0.6455]]

조정 잔차 R^adj
 [[ 0.1903  1.215  -0.9652 -0.4616]
 [-1.8077 -1.0281  3.3783 -0.4616]
 [ 1.6174 -0.1869 -2.4131  0.9232]]

|R| 최대 2.4152,  |R^adj| 최대 3.3783

  표준화 잔차 + 본페로니: 최소 조정 p = 0.1887, 유의한 칸 0개
  표준화 잔차 + 홀름   : 최소 조정 p = 0.1887, 유의한 칸 0개
  표준화 잔차 + BH     : 최소 조정 p = 0.1887, 유의한 칸 0개
  조정   잔차 + 본페로니: 최소 조정 p = 0.0088, 유의한 칸 1개
  조정   잔차 + 홀름   : 최소 조정 p = 0.0088, 유의한 칸 1개
  조정   잔차 + BH     : 최소 조정 p = 0.0088, 유의한 칸 1개

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

rng = np.random.default_rng(1470)
M, n = 3_000, 400

def cell_pvalues(obs):
    """조정 잔차에서 칸별 양측 p 값을 만든다."""
    total = obs.sum()
    exp = np.outer(obs.sum(1), obs.sum(0)) / total
    rp, cp = obs.sum(1) / total, obs.sum(0) / total
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))
    return 2 * stats.norm.sf(np.abs(adj)).ravel()

scenarios = {
    "H0 (독립)": np.outer([1 / 3] * 3, [0.25] * 4),
    "한 칸만 이탈": None,       # 아래에서 만든다
    "여러 칸 이탈": None,
}
base = np.outer([1 / 3] * 3, [0.25] * 4)
p1 = base.copy(); p1[1, 2] += 0.04; p1[1, 0] -= 0.04
p2 = base.copy()
p2[0, 0] += 0.02; p2[0, 1] -= 0.02
p2[1, 2] += 0.02; p2[1, 3] -= 0.02
p2[2, 1] += 0.02; p2[2, 0] -= 0.02
scenarios["한 칸만 이탈"] = p1
scenarios["여러 칸 이탈"] = p2

print(f"{'상황':>14s} {'보정':>8s} {'적어도 하나 기각':>16s} {'평균 기각 칸':>13s}")
for label, P in scenarios.items():
    for method, name in [("bonferroni", "본페로니"), ("holm", "홀름"),
                         ("fdr_bh", "BH"), (None, "보정 없음")]:
        any_rej = tot_rej = 0
        for _ in range(M):
            obs = rng.multinomial(n, P.ravel()).reshape(3, 4).astype(float)
            if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
                continue
            pv = cell_pvalues(obs)
            rej = pv < 0.05 if method is None else \
                multipletests(pv, alpha=0.05, method=method)[0]
            any_rej += rej.any()
            tot_rej += rej.sum()
        print(f"{label:>14s} {name:>8s} {any_rej / M:16.4f} {tot_rej / M:13.4f}")

            상황       보정        적어도 하나 기각       평균 기각 칸
       H0 (독립)     본페로니           0.0453        0.0500
       H0 (독립)       홀름           0.0453        0.0537
       H0 (독립)       BH           0.0510        0.0723
       H0 (독립)    보정 없음           0.3620        0.5923
       한 칸만 이탈     본페로니           0.6363        1.0953
       한 칸만 이탈       홀름           0.6517        1.1263
       한 칸만 이탈       BH           0.6613        1.4780
       한 칸만 이탈    보정 없음           0.9533        2.8807
       여러 칸 이탈     본페로니           0.5123        0.8883
       여러 칸 이탈       홀름           0.5193        0.9187
       여러 칸 이탈       BH           0.5447        1.4087
       여러 칸 이탈    보정 없음           0.9250        3.0647

import numpy as np
from scipy import stats

def show(obs, label):
    obs = np.asarray(obs, float)
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    n = obs.sum()
    rp, cp = obs.sum(1) / n, obs.sum(0) / n
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - rp, 1 - cp))
    contrib = (obs - exp)**2 / exp
    print(f"{label}   n={n:.0f}  χ²={chi2:.4f}  df={df}  p={p:.4f}")
    print(f"  관측\n{obs.astype(int)}")
    print(f"  기대\n{np.round(exp, 2)}")
    print(f"  조정 잔차\n{np.round(adj, 3)}")
    print(f"  χ² 기여율(%)\n{np.round(contrib / chi2 * 100, 1)}\n")

show([[980, 20], [960, 40]], "① 주변합이 한쪽에 쏠림")
show([[500, 500], [400, 600]], "② 주변합이 균형")

① 주변합이 한쪽에 쏠림   n=2000  χ²=6.8729  df=1  p=0.0088
  관측
[[980  20]
 [960  40]]
  기대
[[970.  30.]
 [970.  30.]]
  조정 잔차
[[ 2.622 -2.622]
 [-2.622  2.622]]
  χ² 기여율(%)
[[ 1.5 48.5]
 [ 1.5 48.5]]

② 주변합이 균형   n=2000  χ²=20.2020  df=1  p=0.0000
  관측
[[500 500]
 [400 600]]
  기대
[[450. 550.]
 [450. 550.]]
  조정 잔차
[[ 4.495 -4.495]
 [-4.495  4.495]]
  χ² 기여율(%)
[[27.5 22.5]
 [27.5 22.5]]

obs = np.array([[980, 20], [960, 40]], float)
p1, p2 = obs[0, 1] / obs[0].sum(), obs[1, 1] / obs[1].sum()
print(f"둘째 열 비율: {p1:.4f} 대 {p2:.4f}")
print(f"  위험차 {p2 - p1:+.4f}   위험비 {p2 / p1:.4f}")

둘째 열 비율: 0.0200 대 0.0400
  위험차 +0.0200   위험비 2.0000

① 옴니버스 검정
    χ² 이 유의하지 않으면 → 여기서 멈춘다
          ↓ 유의함
② 표의 크기를 본다
    2×2  → 잔차 분석 불필요 (자유도 1)
    그 외 → 계속
          ↓
③ 조정 표준화 잔차를 계산한다
    R^adj = (O-E) / √[E(1-R_i/n)(1-C_j/n)]
    ※ 단순 표준화 잔차를 쓰면 이탈을 놓친다
          ↓
④ 다중비교 보정 (칸이 r×c 개)
    홀름을 기본, 탐색이면 BH
          ↓
⑤ 살아남은 칸을 원 도수·비율과 함께 보고
          ↓
⑥ 그림으로 전체 패턴 확인 (모자이크 등)