효과크기와 Cramér의 V¶

import numpy as np
from scipy import stats

def cramers_v(observed):
    """분할표의 Cramér V를 계산한다.

    V는 카이제곱을 표본크기로 나눠 정규화한 것이라 0과 1 사이에 놓인다.
    카이제곱 통계량 자체는 n에 비례해 커지므로 연구 사이 비교에 쓸 수 없지만
    V는 쓸 수 있다는 것이 요점이다.
    """
    chi2, p_value, df, expected = stats.chi2_contingency(observed)
    n = observed.sum()
    # q는 행 수와 열 수 중 작은 것에서 1을 뺀 값이다.
    # 카이제곱의 최댓값이 n*q라서 이것으로 나누면 상한이 1이 된다.
    q = min(observed.shape) - 1
    v = np.sqrt(chi2 / (n * q))
    return v, chi2, p_value

# 보기: 성별과 손잡이
observed = np.array([[934, 1070], [113, 92], [20, 8]])
v, chi2, p_value = cramers_v(observed)

print(f"Chi-square statistic: {chi2:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Cramér's V: {v:.4f}")

# (1) 비율로 적은 항등식
n = observed.sum()
P = observed / n
pr, pc = P.sum(1, keepdims=True), P.sum(0, keepdims=True)
print(f"\n항등식  n(S-1) = {n * ((P**2 / (pr * pc)).sum() - 1):.6f}"
      f"   정의대로 = {chi2:.6f}")

# (2) 상한. 완전히 예측되는 표에서 등호가 성립한다.
perfect = np.array([[50, 0], [0, 50], [30, 0]])
c2p = stats.chi2_contingency(perfect)[0]
print(f"완전 예측 표: chi2 = {c2p:.1f}  n(q-1) = {perfect.sum() * 1}"
      f"  V = {np.sqrt(c2p / perfect.sum()):.1f}")

# (3) 유의 문턱을 V 로, 그리고 n 으로
crit = stats.chi2.ppf(0.95, 2)
print(f"\nn = {n:.0f} 에서 유의해지는 V 문턱 {np.sqrt(crit / n):.6f}")
print(f"V = {v:.6f} 가 유의해지는 n   {crit / v**2:.2f}")
print(f"H0 아래 V 의 대략적 크기 sqrt(df/n) = {np.sqrt(2 / n):.6f}")
for m in [300, 1136, 2237, 20000]:
    x = m * v**2
    print(f"  n = {m:6d}:  chi2 = {x:8.3f}   p = {stats.chi2(2).sf(x):.3g}"
          f"   V = {v:.4f}")

Chi-square statistic: 11.8061
p-value: 0.0027
Cramér's V: 0.0726

항등식  n(S-1) = 11.806135   정의대로 = 11.806135
완전 예측 표: chi2 = 130.0  n(q-1) = 130  V = 1.0

n = 2237 에서 유의해지는 V 문턱 0.051753
V = 0.072648 가 유의해지는 n   1135.25
H0 아래 V 의 대략적 크기 sqrt(df/n) = 0.029901
  n =    300:  chi2 =    1.583   p = 0.453   V = 0.0726
  n =   1136:  chi2 =    5.995   p = 0.0499   V = 0.0726
  n =   2237:  chi2 =   11.806   p = 0.00273   V = 0.0726
  n =  20000:  chi2 =  105.553   p = 1.2e-23   V = 0.0726

import numpy as np
from scipy import stats

rng = np.random.default_rng(909)

def cramers_v(obs):
    chi2 = stats.chi2_contingency(obs, correction=False)[0]
    n, q = obs.sum(), min(obs.shape) - 1
    return np.sqrt(chi2 / (n * q))

def cramers_v_corrected(obs):
    """베르그스마(2013)의 편향보정."""
    chi2 = stats.chi2_contingency(obs, correction=False)[0]
    n, (r, c) = obs.sum(), obs.shape
    phi2 = max(0.0, chi2 / n - (r - 1) * (c - 1) / (n - 1))
    r_t = r - (r - 1)**2 / (n - 1)
    c_t = c - (c - 1)**2 / (n - 1)
    return np.sqrt(phi2 / min(r_t - 1, c_t - 1))

print(f"{'표':>7s} {'n':>5s} {'V 평균':>9s} {'V 중앙값':>10s} {'보정 V 평균':>12s}")
layouts = [((2, 2), np.array([.5, .5]), np.array([.5, .5])),
           ((3, 3), np.ones(3) / 3, np.ones(3) / 3),
           ((4, 5), np.ones(4) / 4, np.ones(5) / 5)]
for shape, p_row, p_col in layouts:
    prob = np.outer(p_row, p_col).ravel()      # 완전 독립: 참 V = 0
    for n in [30, 100, 500]:
        vs, vc = [], []
        for _ in range(2_000):
            obs = rng.multinomial(n, prob).reshape(shape).astype(float)
            if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
                continue
            vs.append(cramers_v(obs))
            vc.append(cramers_v_corrected(obs))
        print(f"{str(shape):>7s} {n:5d} {np.mean(vs):9.4f} "
              f"{np.median(vs):10.4f} {np.mean(vc):12.4f}")

      표     n      V 평균      V 중앙값      보정 V 평균
 (2, 2)    30    0.1457     0.1336       0.0624
 (2, 2)   100    0.0814     0.0705       0.0351
 (2, 2)   500    0.0360     0.0313       0.0153
 (3, 3)    30    0.2487     0.2424       0.0810
 (3, 3)   100    0.1355     0.1323       0.0449
 (3, 3)   500    0.0595     0.0580       0.0188
 (4, 5)    30    0.3639     0.3617       0.0887
 (4, 5)   100    0.1975     0.1971       0.0495
 (4, 5)   500    0.0878     0.0874       0.0217

import numpy as np
from scipy import stats
from scipy.optimize import brentq

def ncp_ci(chi2, df, alpha=0.05, big=1e6):
    """비중심 모수 λ 의 신뢰구간.
    sf(chi2; df, λ) 는 λ 에 대해 증가하고 cdf 는 감소한다."""
    lo = (brentq(lambda l: stats.ncx2.sf(chi2, df, l) - alpha / 2, 0.0, big)
          if stats.ncx2.sf(chi2, df, 0.0) < alpha / 2 else 0.0)
    hi = (brentq(lambda l: stats.ncx2.cdf(chi2, df, l) - alpha / 2, 0.0, big)
          if stats.ncx2.cdf(chi2, df, 0.0) > alpha / 2 else 0.0)
    return lo, hi

def v_with_ci(obs, alpha=0.05):
    chi2, p, df, _ = stats.chi2_contingency(obs, correction=False)
    n, q = obs.sum(), min(obs.shape) - 1
    lo, hi = ncp_ci(chi2, df, alpha)
    return (np.sqrt(chi2 / (n * q)),
            np.sqrt(lo / (n * q)), np.sqrt(hi / (n * q)), chi2, p, df)

tables = [np.array([[934, 1070], [113, 92], [20, 8]], float),
          np.array([[12, 8], [6, 14]], float),
          np.array([[40, 10], [10, 40]], float)]
for obs in tables:
    v, lo, hi, chi2, p, df = v_with_ci(obs)
    print(f"{obs.astype(int).tolist()}")
    print(f"   n={obs.sum():.0f}  χ²={chi2:.4f}  df={df}  p={p:.4f}")
    print(f"   V = {v:.4f}   95% CI ({lo:.4f}, {hi:.4f})")

[[934, 1070], [113, 92], [20, 8]]
   n=2237  χ²=11.8061  df=2  p=0.0027
   V = 0.0726   95% CI (0.0261, 0.1116)
[[12, 8], [6, 14]]
   n=40  χ²=3.6364  df=1  p=0.0565
   V = 0.3015   95% CI (0.0000, 0.6114)
[[40, 10], [10, 40]]
   n=100  χ²=36.0000  df=1  p=0.0000
   V = 0.6000   95% CI (0.4040, 0.7960)

rng = np.random.default_rng(55)
P = np.array([[.35, .15], [.15, .35]])
p_row, p_col = P.sum(1), P.sum(0)
v_true = np.sqrt((((P - np.outer(p_row, p_col))**2)
                  / np.outer(p_row, p_col)).sum())
print(f"참 V = {v_true:.4f}")
for n in [50, 200, 1000]:
    cov = ok = 0
    for _ in range(2_000):
        obs = rng.multinomial(n, P.ravel()).reshape(2, 2).astype(float)
        if (obs.sum(0) == 0).any() or (obs.sum(1) == 0).any():
            continue
        ok += 1
        _, lo, hi, *_ = v_with_ci(obs)
        cov += lo <= v_true <= hi
    print(f"  n={n:4d}: 포함률 {cov / ok:.4f}")

참 V = 0.4000
  n=  50: 포함률 0.9560
  n= 200: 포함률 0.9680
  n=1000: 포함률 0.9695

import numpy as np
from scipy import stats

def report(obs, label):
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    n = obs.sum()
    v = np.sqrt(chi2 / (n * (min(obs.shape) - 1)))
    p_row, p_col = obs.sum(1) / n, obs.sum(0) / n
    # 조정 표준화 잔차: H0 아래 근사적으로 N(0,1)
    adj = (obs - exp) / np.sqrt(exp * np.outer(1 - p_row, 1 - p_col))
    print(f"{label}   n={n:.0f}  χ²={chi2:.3f}  df={df}  "
          f"p={p:.4f}  V={v:.4f}")
    print(f"  조정 표준화 잔차\n{np.round(adj, 2)}")
    print(f"  |잔차| > 2 인 칸 {int((np.abs(adj) > 2).sum())}개, "
          f"최대 {np.abs(adj).max():.2f}\n")

concentrated = np.array([[102, 45, 45], [45, 45, 45], [45, 45, 45]], float)
diffuse = np.array([[64, 43, 43], [43, 64, 43], [43, 43, 64]], float)
report(concentrated, "집중형 (한 칸만 큼)")
report(diffuse, "분산형 (대각선 전체)")

집중형 (한 칸만 큼)   n=462  χ²=18.097  df=4  p=0.0012  V=0.1399
  조정 표준화 잔차
[[ 4.25 -2.3  -2.3 ]
 [-2.3   1.25  1.25]
 [-2.3   1.25  1.25]]
  |잔차| > 2 인 칸 5개, 최대 4.25

분산형 (대각선 전체)   n=450  χ²=17.640  df=4  p=0.0015  V=0.1400
  조정 표준화 잔차
[[ 2.97 -1.48 -1.48]
 [-1.48  2.97 -1.48]
 [-1.48 -1.48  2.97]]
  |잔차| > 2 인 칸 3개, 최대 2.97

import numpy as np
from scipy import stats

def measures(obs, label):
    obs = obs.astype(float)
    n = obs.sum()
    chi2, p, df, _ = stats.chi2_contingency(obs, correction=False)
    r, c = obs.shape
    v = np.sqrt(chi2 / (n * (min(r, c) - 1)))
    w = np.sqrt(chi2 / n)                                  # 코헨 w
    cc = np.sqrt(chi2 / (chi2 + n))                        # 분할계수
    t = np.sqrt((chi2 / n) / np.sqrt((r - 1) * (c - 1)))   # 추프로프 T

    P = obs / n
    p_row, p_col = P.sum(1), P.sum(0)
    h_row = -(p_row * np.log(p_row)).sum()
    h_col = -(p_col * np.log(p_col)).sum()
    h_joint = -(P[P > 0] * np.log(P[P > 0])).sum()
    u = 2 * (h_row + h_col - h_joint) / (h_row + h_col)     # 불확실성 계수

    print(f"{label}   n={n:.0f}  χ²={chi2:.3f}  p={p:.4f}")
    print(f"  V {v:.4f} | 코헨 w {w:.4f} | 분할계수 C {cc:.4f} | "
          f"추프로프 T {t:.4f} | 불확실성 U {u:.4f}")
    if obs.shape == (2, 2):
        a, b, c_, d = obs.ravel()
        print(f"  오즈비 {(a * d) / (b * c_):.4f} | "
              f"위험비 {(a / (a + b)) / (c_ / (c_ + d)):.4f} | "
              f"위험차 {a / (a + b) - c_ / (c_ + d):+.4f}")

measures(np.array([[40, 10], [10, 40]]), "강한 연관 (2×2)")
measures(np.array([[55, 45], [45, 55]]), "약한 연관 (2×2)")
measures(np.array([[934, 1070], [113, 92], [20, 8]]), "보기 1 (3×2)")

강한 연관 (2×2)   n=100  χ²=36.000  p=0.0000
  V 0.6000 | 코헨 w 0.6000 | 분할계수 C 0.5145 | 추프로프 T 0.6000 | 불확실성 U 0.2781
  오즈비 16.0000 | 위험비 4.0000 | 위험차 +0.6000
약한 연관 (2×2)   n=200  χ²=2.000  p=0.1573
  V 0.1000 | 코헨 w 0.1000 | 분할계수 C 0.0995 | 추프로프 T 0.1000 | 불확실성 U 0.0072
  오즈비 1.4938 | 위험비 1.2222 | 위험차 +0.1000
보기 1 (3×2)   n=2237  χ²=11.806  p=0.0027
  V 0.0726 | 코헨 w 0.0726 | 분할계수 C 0.0725 | 추프로프 T 0.0611 | 불확실성 U 0.0050

import numpy as np
from scipy import stats

def n_for_v(v, r, c, power=0.80, alpha=0.05):
    df, q = (r - 1) * (c - 1), min(r, c) - 1
    crit = stats.chi2.ppf(1 - alpha, df)
    for n in range(10, 200_000):
        if stats.ncx2.sf(crit, df, n * v * v * q) >= power:
            return n

header = " ".join(f"{'V=' + str(v):>8s}" for v in [0.1, 0.2, 0.3, 0.5])
print(f"{'표 크기':>8s} {'df':>4s} {header}")
for r, c in [(2, 2), (2, 3), (3, 3), (4, 5), (5, 5)]:
    row = " ".join(f"{n_for_v(v, r, c):8d}" for v in [0.1, 0.2, 0.3, 0.5])
    print(f"{f'{r}×{c}':>8s} {(r - 1) * (c - 1):4d} {row}")

print("\n보기 1 (3×2, V=0.0726, n=2237) 에 대해")
for v, label in [(0.05, "무시할 만함"), (0.10, "작음"), (0.30, "중간")]:
    print(f"  참 V={v} ({label}) 를 80% 로 탐지하려면 n = {n_for_v(v, 3, 2)}")

    표 크기   df    V=0.1    V=0.2    V=0.3    V=0.5
     2×2    1      785      197       88       32
     2×3    2      964      241      108       39
     3×3    4      597      150       67       24
     4×5   12      578      145       65       24
     5×5   16      482      121       54       20

보기 1 (3×2, V=0.0726, n=2237) 에 대해
  참 V=0.05 (무시할 만함) 를 80% 로 탐지하려면 n = 3854
  참 V=0.1 (작음) 를 80% 로 탐지하려면 n = 964
  참 V=0.3 (중간) 를 80% 로 탐지하려면 n = 108