동질성 검정 (scipy)¶

import numpy as np
from scipy import stats

# 모집단 3개(행), 범주 4개(열)
# 동질성 검정과 독립성 검정은 **계산이 완전히 같다**. 다른 것은 표집 설계뿐이다.
# 여기서는 모집단마다 표본을 따로 뽑았으므로 행 합계가 설계자에 의해 고정되어 있다.
observed = np.array([
    [25, 30, 20, 25],   # 모집단 1
    [18, 22, 35, 25],   # 모집단 2
    [30, 25, 15, 30],   # 모집단 3
], dtype=float)

chi2, p, df, expected = stats.chi2_contingency(observed, correction=False)

print("=== Chi-square Test of Homogeneity (scipy) ===")
print(f"chi2 = {chi2:.4f}, df = {df}, p = {p:.6f}")
print()
print("Expected counts under H0 (same proportions):")
print(expected)

=== Chi-square Test of Homogeneity (scipy) ===
chi2 = 14.1697, df = 6, p = 0.027796

Expected counts under H0 (same proportions):
[[24.33333333 25.66666667 23.33333333 26.66666667]
 [24.33333333 25.66666667 23.33333333 26.66666667]
 [24.33333333 25.66666667 23.33333333 26.66666667]]

import numpy as np

O = np.array([[25, 30, 20, 25],
              [18, 22, 35, 25],
              [30, 25, 15, 30]], dtype=float)
R, C, n = O.sum(1), O.sum(0), O.sum()
E = np.outer(R, C) / n
contrib = (O - E) ** 2 / E
pearson = (O - E) / np.sqrt(E)
adj = (O - E) / np.sqrt(E * (1 - R[:, None] / n) * (1 - C[None, :] / n))

for name, M in [("칸별 기여  (O-E)^2/E", contrib),
                ("피어슨 잔차  (O-E)/sqrt(E)", pearson),
                ("조정 잔차", adj)]:
    print(name)
    for i, row in enumerate(M):
        print(f"  모집단 {i + 1}  " + " ".join(f"{v:+8.4f}" for v in row))
print("\n열별 기여 합  " + " ".join(f"{v:8.4f}" for v in contrib.sum(0)))
print(f"전체 chi2 = {contrib.sum():.4f}   가장 큰 칸의 몫 "
      f"{contrib.max() / contrib.sum():.3f}   범주 3 열의 몫 "
      f"{contrib.sum(0)[2] / contrib.sum():.3f}")
print(f"크래머 V = {np.sqrt(contrib.sum() / (n * 2)):.4f},   최소 기대도수 {E.min():.4f}")

칸별 기여  (O-E)^2/E
  모집단 1   +0.0183  +0.7316  +0.4762  +0.1042
  모집단 2   +1.6484  +0.5238  +5.8333  +0.1042
  모집단 3   +1.3196  +0.0173  +2.9762  +0.4167
피어슨 잔차  (O-E)/sqrt(E)
  모집단 1   +0.1351  +0.8553  -0.6901  -0.3227
  모집단 2   -1.2839  -0.7237  +2.4152  -0.3227
  모집단 3   +1.1488  -0.1316  -1.7252  +0.6455
조정 잔차
  모집단 1   +0.1903  +1.2150  -0.9652  -0.4616
  모집단 2   -1.8077  -1.0281  +3.3783  -0.4616
  모집단 3   +1.6174  -0.1869  -2.4131  +0.9232

열별 기여 합    2.9863   1.2727   9.2857   0.6250
전체 chi2 = 14.1697   가장 큰 칸의 몫 0.412   범주 3 열의 몫 0.655
크래머 V = 0.1537,   최소 기대도수 23.3333

import numpy as np
from scipy import stats

observed = np.array([[30,20],[28,22],[35,15],[25,25],[32,18]], dtype=float)
chi2, p, df, expected = stats.chi2_contingency(observed, correction=False)
print(f"chi2 = {chi2:.4f}, p = {p:.4f}, df = {df}")

chi2 = 4.8333, p = 0.3048, df = 4

import numpy as np
from scipy import stats

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

n = obs.sum()
row_tot, col_tot = obs.sum(1), obs.sum(0)
print(f"행 합 {row_tot.tolist()},  열 합 {col_tot.tolist()},  n = {n:.0f}")

exp = np.outer(row_tot, col_tot) / n          # E_ij = R_i C_j / n
print(f"\n기대도수\n{np.round(exp, 4)}")

chi2 = ((obs - exp)**2 / exp).sum()
df = (obs.shape[0] - 1) * (obs.shape[1] - 1)
print(f"\n손계산  χ² = {chi2:.4f},  df = {df},  "
      f"p = {stats.chi2.sf(chi2, df):.6f}")
c, p, d, e = stats.chi2_contingency(obs, correction=False)
print(f"scipy   χ² = {c:.4f},  df = {d},  p = {p:.6f}")
print(f"일치? {np.isclose(chi2, c)} / {np.isclose(exp, e).all()}")

행 합 [100.0, 100.0, 100.0],  열 합 [73.0, 77.0, 70.0, 80.0],  n = 300

기대도수
[[24.3333 25.6667 23.3333 26.6667]
 [24.3333 25.6667 23.3333 26.6667]
 [24.3333 25.6667 23.3333 26.6667]]

손계산  χ² = 14.1697,  df = 6,  p = 0.027796
scipy   χ² = 14.1697,  df = 6,  p = 0.027796
일치? True / True

unequal = np.array([[25, 30, 20, 25],      # n=100
                    [36, 44, 70, 50],      # n=200
                    [15, 12,  8, 15]],     # n=50
                   dtype=float)
_, _, _, e2 = stats.chi2_contingency(unequal, correction=False)
print("표본크기가 다를 때의 기대도수")
print(np.round(e2, 3))
print(f"행 합 {unequal.sum(1).tolist()}")
print(f"각 행의 기대도수 합 {np.round(e2.sum(1), 3).tolist()}  ← 행 합과 일치")

표본크기가 다를 때의 기대도수
[[21.714 24.571 28.    25.714]
 [43.429 49.143 56.    51.429]
 [10.857 12.286 14.    12.857]]
행 합 [100.0, 200.0, 50.0]
각 행의 기대도수 합 [100.0, 200.0, 50.0]  ← 행 합과 일치

import numpy as np
from scipy import stats

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

a = b = 0
for _ in range(M):
    # ① 독립성 표집: 전체 n 을 한 번에 다항으로
    t1 = rng.multinomial(n, P.ravel()).reshape(3, 4).astype(float)
    if (t1.sum(0) > 0).all() and (t1.sum(1) > 0).all():
        a += stats.chi2_contingency(t1, correction=False)[1] < 0.05

    # ② 동질성 표집: 행마다 100 씩 따로
    t2 = np.array([rng.multinomial(100, p_col) for _ in range(3)], float)
    if (t2.sum(0) > 0).all():
        b += stats.chi2_contingency(t2, correction=False)[1] < 0.05

print(f"독립성 표집(총 n 만 고정):  기각률 {a / M:.4f}")
print(f"동질성 표집(행 합 고정):    기각률 {b / M:.4f}")

독립성 표집(총 n 만 고정):  기각률 0.0521
동질성 표집(행 합 고정):    기각률 0.0508

print("동질성 설계에서 추정 가능한 것과 불가능한 것")
obs = np.array([[25, 30, 20, 25], [18, 22, 35, 25], [30, 25, 15, 30]], float)
print(f"  집단별 범주 비율(가능):\n{np.round(obs / obs.sum(1, keepdims=True), 4)}")
print(f"  전체 집단 비율(불가능): {np.round(obs.sum(1) / obs.sum(), 4).tolist()}")
print("    ← 이 값은 연구자가 1:1:1 로 정한 것이지 모집단의 성질이 아니다")

동질성 설계에서 추정 가능한 것과 불가능한 것
  집단별 범주 비율(가능):
[[0.25 0.3  0.2  0.25]
 [0.18 0.22 0.35 0.25]
 [0.3  0.25 0.15 0.3 ]]
  전체 집단 비율(불가능): [0.3333, 0.3333, 0.3333]
    ← 이 값은 연구자가 1:1:1 로 정한 것이지 모집단의 성질이 아니다

import numpy as np
from scipy import stats

obs = np.array([[25, 30, 20, 25],
                [18, 22, 35, 25],
                [30, 25, 15, 30]], dtype=float)
k = obs.shape[0]

z_single = stats.norm.ppf(0.975)
z_bonf = stats.norm.ppf(1 - 0.025 / k)                 # 본페로니
z_sidak = stats.norm.ppf((1 + 0.95**(1 / k)) / 2)      # 시닥
print(f"임계값:  단일 {z_single:.4f}   본페로니 {z_bonf:.4f}   "
      f"시닥 {z_sidak:.4f}\n")

print("범주 3 의 집단별 비율")
for i in range(k):
    n_i = obs[i].sum()
    p_hat = obs[i, 2] / n_i
    se = np.sqrt(p_hat * (1 - p_hat) / n_i)
    print(f"  모집단{i + 1}: p̂ = {p_hat:.4f}   "
          f"단일 ({p_hat - z_single * se:.4f}, {p_hat + z_single * se:.4f})   "
          f"본페로니 ({p_hat - z_bonf * se:.4f}, {p_hat + z_bonf * se:.4f})")

임계값:  단일 1.9600   본페로니 2.3940   시닥 2.3877

범주 3 의 집단별 비율
  모집단1: p̂ = 0.2000   단일 (0.1216, 0.2784)   본페로니 (0.1042, 0.2958)
  모집단2: p̂ = 0.3500   단일 (0.2565, 0.4435)   본페로니 (0.2358, 0.4642)
  모집단3: p̂ = 0.1500   단일 (0.0800, 0.2200)   본페로니 (0.0645, 0.2355)

from itertools import combinations
for i, j in combinations(range(k), 2):
    n_i, n_j = obs[i].sum(), obs[j].sum()
    p_i, p_j = obs[i, 2] / n_i, obs[j, 2] / n_j
    se = np.sqrt(p_i * (1 - p_i) / n_i + p_j * (1 - p_j) / n_j)
    lo, hi = (p_i - p_j) - z_bonf * se, (p_i - p_j) + z_bonf * se
    mark = "*" if not (lo <= 0 <= hi) else " "
    print(f"  모집단{i + 1} - 모집단{j + 1}: 차이 {p_i - p_j:+.4f}   "
          f"본페로니 95% CI ({lo:+.4f}, {hi:+.4f}){mark}")

  모집단1 - 모집단2: 차이 -0.1500   본페로니 95% CI (-0.2990, -0.0010)*
  모집단1 - 모집단3: 차이 +0.0500   본페로니 95% CI (-0.0784, +0.1784)
  모집단2 - 모집단3: 차이 +0.2000   본페로니 95% CI (+0.0574, +0.3426)*

import numpy as np
from scipy import stats

def max_table(r, c, n):
    """완전 연관인 표: 작은 쪽의 각 수준이 큰 쪽의 서로 다른 묶음에 대응."""
    obs = np.zeros((r, c))
    if r <= c:
        for i, g in enumerate(np.array_split(np.arange(c), r)):
            obs[i, g] = n / (r * len(g))
    else:
        for j, g in enumerate(np.array_split(np.arange(r), c)):
            obs[g, j] = n / (c * len(g))
    return obs

def max_check(r, c, n):
    q = min(r, c)
    chi2, _, _, _ = stats.chi2_contingency(max_table(r, c, n),
                                           correction=False)
    return chi2, n * (q - 1), np.sqrt(chi2 / (n * (q - 1)))

print(f"{'표':>7s} {'n':>5s} {'관측 χ²':>10s} {'n(q-1)':>9s} {'V':>7s}")
for r, c, n in [(2, 2, 100), (3, 3, 300), (4, 4, 400),
                (2, 4, 200), (4, 2, 200), (3, 5, 300)]:
    got, bound, v = max_check(r, c, n)
    print(f"{f'{r}×{c}':>7s} {n:5d} {got:10.4f} {bound:9.0f} {v:7.4f}")

      표     n      관측 χ²    n(q-1)       V
    2×2   100   100.0000       100  1.0000
    3×3   300   600.0000       600  1.0000
    4×4   400  1200.0000      1200  1.0000
    2×4   200   200.0000       200  1.0000
    4×2   200   200.0000       200  1.0000
    3×5   300   600.0000       600  1.0000

import numpy as np
from scipy import stats

def homogeneity_test(obs, labels=None, alpha=0.05):
    """동질성 검정과 함께 진단 정보를 한 번에 돌려준다."""
    obs = np.asarray(obs, float)
    r, c = obs.shape
    n = obs.sum()
    chi2, p, df, exp = stats.chi2_contingency(obs, correction=False)
    v = np.sqrt(chi2 / (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))

    print(f"χ² = {chi2:.4f},  df = {df},  p = {p:.6f}")
    print(f"크라메르 V = {v:.4f}   (독립일 때 기댓값 ≈ "
          f"{np.sqrt(df / (n * (min(r, c) - 1))):.4f})")
    print(f"E_min = {exp.min():.3f}"
          + ("   ⚠ 5 미만" if exp.min() < 5 else ""))
    print(f"행 비율\n{np.round(obs / obs.sum(1, keepdims=True), 4)}")
    if p < alpha:
        print(f"조정 잔차 (|·| 최대 {np.abs(adj).max():.3f})\n{np.round(adj, 3)}")
    return {"chi2": chi2, "df": df, "p": p, "V": v, "adj": adj}

_ = homogeneity_test([[25, 30, 20, 25],
                      [18, 22, 35, 25],
                      [30, 25, 15, 30]])

χ² = 14.1697,  df = 6,  p = 0.027796
크라메르 V = 0.1537   (독립일 때 기댓값 ≈ 0.1000)
E_min = 23.333
행 비율
[[0.25 0.3  0.2  0.25]
 [0.18 0.22 0.35 0.25]
 [0.3  0.25 0.15 0.3 ]]
조정 잔차 (|·| 최대 3.378)
[[ 0.19   1.215 -0.965 -0.462]
 [-1.808 -1.028  3.378 -0.462]
 [ 1.617 -0.187 -2.413  0.923]]