동질성 검정¶

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def compute_expected(observed):
    """행과 열이 무관하다는 가정 아래 기대도수를 구한다.

    무관하다면 결합확률이 주변확률의 곱이 된다. 그 곱에 전체 도수를 곱한
    것이 각 칸의 기대도수다.
    """
    row_sum = observed.sum(axis=1)
    row_pmf = row_sum.reshape((-1, 1)) / row_sum.sum()

    column_sum = observed.sum(axis=0)
    column_pmf = column_sum.reshape((1, -1)) / column_sum.sum()

    joint_pmf = row_pmf * column_pmf
    expected = joint_pmf * row_sum.sum()
    return expected

def main():
    """동질성 검정을 정의대로 계산하고 p-값을 그림으로 보인다."""
    # 행이 나라, 열이 응답 범주다. 나라마다 응답 분포가 같은지를 묻는다.
    observed = np.array([[541, 75, 231], [498, 71, 213],
                         [779, 96, 321], [282, 50, 345], [65, 19, 120]])
    expected = compute_expected(observed)

    # 자유도는 (행-1)(열-1). 행합과 열합이 묶여 있어 자유롭게 움직일 수 있는
    # 칸이 그만큼뿐이다.
    df = (observed.shape[0] - 1) * (observed.shape[1] - 1)

    statistic = np.sum((observed - expected)**2 / expected)
    p_value = stats.chi2(df).sf(statistic)
    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.02%}")

    _, ax = plt.subplots(figsize=(12, 4))

    x = np.linspace(0, statistic, 1000)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='b', linewidth=3)

    x = np.concatenate([[0], x, [statistic], [0]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='b', alpha=0.1)

    x = np.linspace(statistic, 300, 100)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='r', linewidth=3)

    x = np.concatenate([[statistic], x, [20], [statistic]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='r', alpha=0.1)

    xy = (250, 0.01)
    xytext = (250, 0.08)
    arrowprops = dict(color='k', width=0.2, headwidth=8)
    ax.annotate(f'{p_value = :.02%}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.spines['left'].set_position("zero")

    plt.show()

if __name__ == "__main__":
    main()

statistic = 212.94
p_value   = 0.00%

import numpy as np

O = np.array([[541, 75, 231], [498, 71, 213],
              [779, 96, 321], [282, 50, 345], [65, 19, 120]])
R, C, N = O.sum(1), O.sum(0), O.sum()
E = np.outer(R, C) / N
contrib = (O - E) ** 2 / E

print("칸별 기여  (O-E)^2/E")
print(f"{'':8s}{'US':>9s}{'Canada':>9s}{'Mexico':>9s}{'행 합':>9s}")
for lab, row in zip(["별 5개", "별 4개", "별 3개", "별 2개", "별 1개"], contrib):
    print(f"{lab:8s}" + "".join(f"{v:9.2f}" for v in row) + f"{row.sum():9.2f}")
print(f"{'열 합':8s}" + "".join(f"{v:9.2f}" for v in contrib.sum(0))
      + f"{contrib.sum():9.2f}")
print(f"\n최소 기대도수 {E.min():.2f}")
print(f"크래머 V = sqrt(chi2 / (n min(r-1,c-1))) "
      f"= sqrt({contrib.sum():.2f} / (3706 x 2)) = {np.sqrt(contrib.sum() / (N * 2)):.4f}")

칸별 기여  (O-E)^2/E
               US   Canada   Mexico      행 합
별 5개         4.31     0.22     8.93    13.46
별 4개         3.71     0.44     8.35    12.50
별 3개         9.23     0.19    14.53    23.95
별 2개        32.57     0.82    64.42    97.80
별 1개        24.63     0.21    40.39    65.22
열 합         74.45     1.87   136.62   212.94

최소 기대도수 17.12
크래머 V = sqrt(chi2 / (n min(r-1,c-1))) = sqrt(212.94 / (3706 x 2)) = 0.1695

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def main():
    """같은 계산을 scipy 의 chi2_contingency 로 대신한다."""
    observed = np.array([[541, 75, 231], [498, 71, 213],
                         [779, 96, 321], [282, 50, 345], [65, 19, 120]])

    # 기대도수와 자유도까지 함께 돌려준다. 앞의 손계산과 값이 맞아야 한다.
    statistic, p_value, df, expected = stats.chi2_contingency(observed)
    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.02%}")

    _, ax = plt.subplots(figsize=(12, 4))

    x = np.linspace(0, statistic, 1000)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='b', linewidth=3)

    x = np.concatenate([[0], x, [statistic], [0]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='b', alpha=0.1)

    x = np.linspace(statistic, 300, 100)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='r', linewidth=3)

    x = np.concatenate([[statistic], x, [20], [statistic]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='r', alpha=0.1)

    xy = (250, 0.01)
    xytext = (250, 0.08)
    arrowprops = dict(color='k', width=0.2, headwidth=8)
    ax.annotate(f'{p_value = :.02%}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.spines['left'].set_position("zero")

    plt.show()

if __name__ == "__main__":
    main()

statistic = 212.94
p_value   = 0.00%

from scipy.stats import chi2

chi2_obs = 212.93563265
print(f"정확한 p-값          {chi2.sf(chi2_obs, 8):.4e}")
print(f"임계값 chi2_0.05(8)  {chi2.ppf(0.95, 8):.4f}")
print(f"관측값은 임계값의    {chi2_obs / chi2.ppf(0.95, 8):.1f} 배")

정확한 p-값          1.1951e-41
임계값 chi2_0.05(8)  15.5073
관측값은 임계값의    13.7 배

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def compute_expected(observed):
    row_sum = observed.sum(axis=1)
    row_pmf = row_sum.reshape((-1, 1)) / row_sum.sum()
    column_sum = observed.sum(axis=0)
    column_pmf = column_sum.reshape((1, -1)) / column_sum.sum()
    joint_pmf = row_pmf * column_pmf
    expected = joint_pmf * row_sum.sum()
    return expected

def main():
    # 이번에는 세 나라의 분포를 비슷하게 맞춘 자료다. 앞과 같은 절차인데
    # p-값이 크게 나온다. 검정이 무엇에 반응하는지가 이 대비에서 드러난다.
    observed = np.array([[541, 530, 550], [498, 490, 503],
                         [779, 750, 760], [282, 270, 265], [65, 60, 58]])
    expected = compute_expected(observed)
    df = (observed.shape[0] - 1) * (observed.shape[1] - 1)

    statistic = np.sum((observed - expected)**2 / expected)
    p_value = stats.chi2(df).sf(statistic)
    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.02%}")

if __name__ == "__main__":
    main()

statistic = 1.10
p_value   = 99.75%

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def main():
    """열을 둘로 줄여 두 나라만 견준다."""
    observed = np.array([[541, 75], [498, 71], [779, 96], [282, 50], [65, 19]])

    statistic, p_value, df, expected = stats.chi2_contingency(observed)
    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.02%}")
    print("expected")
    print(expected)

if __name__ == "__main__":
    main()

statistic = 11.73
p_value   = 1.95%
expected
[[538.62681745  77.37318255]
 [497.53029079  71.46970921]
 [765.09491115 109.90508885]
 [290.29886914  41.70113086]
 [ 73.44911147  10.55088853]]

import itertools
import numpy as np
from scipy import stats

O = np.array([[541, 75, 231], [498, 71, 213],
              [779, 96, 321], [282, 50, 345], [65, 19, 120]])
names = ["US", "Canada", "Mexico"]

print("쌍마다 열 둘만 남겨 다시 검정한다 (df = 4)")
total = 0.0
for i, j in itertools.combinations(range(3), 2):
    res = stats.chi2_contingency(O[:, [i, j]])
    total += res.statistic
    print(f"  {names[i]:>6s} - {names[j]:<6s}  chi2 = {res.statistic:8.4f}   "
          f"p = {res.pvalue:9.3g}   최소기대 {res.expected_freq.min():6.2f}")
print(f"\n세 쌍의 합 {total:.4f}   전체 검정 212.9356")
print(f"본페로니 임계 p = 0.05/3 = {0.05 / 3:.4f}")

쌍마다 열 둘만 남겨 다시 검정한다 (df = 4)
      US - Canada  chi2 =  11.7255   p =    0.0195   최소기대  10.55
      US - Mexico  chi2 = 210.5657   p =  2.01e-44   최소기대  67.03
  Canada - Mexico  chi2 =  27.2786   p =  1.75e-05   최소기대  28.05

세 쌍의 합 249.5698   전체 검정 212.9356
본페로니 임계 p = 0.05/3 = 0.0167

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def main():
    observed = np.array([[30, 10], [15, 25], [15, 5]])

    statistic, p_value, df, expected = stats.chi2_contingency(observed)
    print(f"{statistic = :.02f}")
    print(f"{p_value   = :.04f}")
    print(f"\nExpected frequencies:")
    print(expected)

    _, ax = plt.subplots(figsize=(12, 4))

    x = np.linspace(0, statistic)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='b', linewidth=3)

    x = np.concatenate([[0], x, [statistic], [0]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='b', alpha=0.1)

    x = np.linspace(statistic, 20, 100)
    y = stats.chi2(df).pdf(x)
    ax.plot(x, y, color='r', linewidth=3)

    x = np.concatenate([[statistic], x, [20], [statistic]])
    y = np.concatenate([[0], y, [0], [0]])
    ax.fill(x, y, color='r', alpha=0.1)

    xy = (15.0, 0.01)
    xytext = (16.5, 0.10)
    arrowprops = dict(color='k', width=0.2, headwidth=8)
    ax.annotate(f'{p_value = :.04f}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.spines['left'].set_position("zero")

    plt.show()

if __name__ == "__main__":
    main()

statistic = 14.06
p_value   = 0.0009

Expected frequencies:
[[24. 16.]
 [24. 16.]
 [12.  8.]]

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

table = np.array([[60, 40], [55, 45], [45, 55]], float)
labels = ["A", "B", "C"]

chi2, p, df, exp = stats.chi2_contingency(table, correction=False)
print(f"전체 동질성:  χ² = {chi2:.4f},  df = {df},  p = {p:.4f}\n")

print("쌍별 비교 (보정 전)")
pvals, pairs = [], []
for i, j in combinations(range(3), 2):
    c, pp, _, _ = stats.chi2_contingency(table[[i, j]], correction=False)
    pvals.append(pp)
    pairs.append(f"{labels[i]}-{labels[j]}")
    print(f"  {labels[i]}-{labels[j]}:  χ² = {c:.4f},  p = {pp:.4f}")

pvals = np.array(pvals)
print()
for method, name in [("bonferroni", "본페로니"), ("holm", "홀름  "),
                     ("fdr_bh", "BH    ")]:
    rej, adj, _, _ = multipletests(pvals, alpha=0.05, method=method)
    detail = "  ".join(f"{pairs[k]} {adj[k]:.4f}{'*' if rej[k] else ' '}"
                       for k in range(3))
    print(f"  {name}: {detail}")

전체 동질성:  χ² = 4.6875,  df = 2,  p = 0.0960

쌍별 비교 (보정 전)
  A-B:  χ² = 0.5115,  p = 0.4745
  A-C:  χ² = 4.5113,  p = 0.0337
  B-C:  χ² = 2.0000,  p = 0.1573

  본페로니: A-B 1.0000   A-C 0.1010   B-C 0.4719
  홀름  : A-B 0.4745   A-C 0.1010   B-C 0.3146
  BH    : A-B 0.4745   A-C 0.1010   B-C 0.2359

def n_per_group_for(p1, p2, power=0.80, alpha=0.05):
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    pbar = (p1 + p2) / 2
    num = (za * np.sqrt(2 * pbar * (1 - pbar))
           + zb * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2)))**2
    return int(np.ceil(num / (p1 - p2)**2))

print(f"A 대 C (0.60 대 0.45) 를 80% 로 탐지: 군당 "
      f"{n_per_group_for(0.60, 0.45)}명")
print(f"A 대 B (0.60 대 0.55) 를 80% 로 탐지: 군당 "
      f"{n_per_group_for(0.60, 0.55)}명")

A 대 C (0.60 대 0.45) 를 80% 로 탐지: 군당 173명
A 대 B (0.60 대 0.55) 를 80% 로 탐지: 군당 1534명

import numpy as np
from scipy import stats

def cochran_armitage(success, total, scores=None):
    """순서형 집단에 대한 선형 추세 검정. (z, 양측 p) 를 돌려준다."""
    success = np.asarray(success, float)
    total = np.asarray(total, float)
    x = (np.arange(len(success), dtype=float) if scores is None
         else np.asarray(scores, float))
    N, S = total.sum(), success.sum()
    p_bar = S / N
    t = (success * x).sum() - p_bar * (total * x).sum()
    var = p_bar * (1 - p_bar) * ((total * x**2).sum()
                                 - (total * x).sum()**2 / N)
    z = t / np.sqrt(var)
    return z, 2 * stats.norm.sf(abs(z))

tot = [50, 50, 50, 50]
for label, succ in [("단조 증가", [10, 14, 18, 24]),
                    ("비단조(뒤섞음)", [10, 24, 14, 18])]:
    table = np.array([[s, t - s] for s, t in zip(succ, tot)], float)
    c, p, df, _ = stats.chi2_contingency(table, correction=False)
    z, pz = cochran_armitage(succ, tot)
    print(f"{label}  성공 {succ} / 각 50")
    print(f"  일반 동질성   χ² = {c:.4f}, df = {df}, p = {p:.4f}")
    print(f"  추세검정      z = {z:.4f}, df = 1, p = {pz:.4f}\n")

단조 증가  성공 [10, 14, 18, 24] / 각 50
  일반 동질성   χ² = 9.6789, df = 3, p = 0.0215
  추세검정      z = 3.0936, df = 1, p = 0.0020

비단조(뒤섞음)  성공 [10, 24, 14, 18] / 각 50
  일반 동질성   χ² = 9.6789, df = 3, p = 0.0215
  추세검정      z = 0.9415, df = 1, p = 0.3464

rng = np.random.default_rng(1357)
M = 5_000
for label, ps in [("단조 증가 (0.20→0.44)", [0.20, 0.28, 0.36, 0.44]),
                  ("V자 (0.35,0.20,0.20,0.35)", [0.35, 0.20, 0.20, 0.35])]:
    a = b = 0
    for _ in range(M):
        s = [rng.binomial(50, q) for q in ps]
        table = np.array([[x, 50 - x] for x in s], float)
        if table.sum(0).min() > 0:
            a += stats.chi2_contingency(table, correction=False)[1] < 0.05
        b += cochran_armitage(s, [50] * 4)[1] < 0.05
    print(f"{label:>26s}:  일반 χ² {a / M:.4f}   추세검정 {b / M:.4f}")

         단조 증가 (0.20→0.44):  일반 χ² 0.6220   추세검정 0.7768
  V자 (0.35,0.20,0.20,0.35):  일반 χ² 0.4856   추세검정 0.0592

import numpy as np
from scipy import stats

rng = np.random.default_rng(2468)
M, m, n_cluster = 5_000, 10, 10      # 군당 군집 10개 × 크기 10 = 100명

print(f"{'ICC':>5s} {'DEFF':>6s} {'보정 없음':>10s} {'DEFF 로 나눔':>13s} "
      f"{'χ²/df 평균':>11s}")
for icc in [0.0, 0.05, 0.10, 0.20]:
    deff = 1 + (m - 1) * icc
    a = b = 0.5 * (1 / icc - 1) if icc > 0 else None
    raw = adj = 0
    phis = []
    for _ in range(M):
        table = np.zeros((3, 2))
        for g in range(3):
            for _ in range(n_cluster):
                p_c = rng.beta(a, b) if icc > 0 else 0.5
                k = rng.binomial(m, p_c)
                table[g, 0] += k
                table[g, 1] += m - k
        if table.sum(0).min() > 0:
            chi2, p, df, _ = stats.chi2_contingency(table, correction=False)
            raw += p < 0.05
            adj += stats.chi2.sf(chi2 / deff, df) < 0.05
            phis.append(chi2 / df)
    print(f"{icc:5.2f} {deff:6.2f} {raw / M:10.4f} {adj / M:13.4f} "
          f"{np.mean(phis):11.4f}")

  ICC   DEFF      보정 없음     DEFF 로 나눔    χ²/df 평균
 0.00   1.00     0.0524        0.0524      1.0084
 0.05   1.45     0.1208        0.0466      1.4339
 0.10   1.90     0.2092        0.0516      1.9139
 0.20   2.80     0.3520        0.0510      2.8554

여러 집단의 분포가 같은지 묻는다
    │
    ├─ ① 표집 설계 확인
    │     ├─ 집단별로 따로 표집 → 동질성 검정
    │     ├─ 한 표본을 두 변수로 분류 → 독립성 검정 (계산은 같다)
    │     └─ 같은 개체를 반복 측정 → McNemar / 코크런 Q
    │
    ├─ ② 독립성 확인
    │     └─ 군집 구조가 있으면 DEFF 보정 (연습문제 9)
    │
    ├─ ③ 집단에 순서가 있는가
    │     ├─ 예, 단조 관계 예상 → 추세검정 (연습문제 8)
    │     └─ 아니오 → 일반 동질성 검정
    │
    ├─ ④ 기대도수 확인 → 필요하면 병합·정확검정
    │
    ├─ ⑤ 옴니버스 검정
    │     └─ 기각 못 하면 여기서 멈춘다
    │
    ├─ ⑥ 사후분석 (연습문제 7)
    │     ├─ 쌍별 비교 + 다중비교 보정
    │     └─ 또는 조정 표준화 잔차
    │
    └─ ⑦ 효과크기 + 신뢰구간 + 원 도수표