Cochran의 Q 검정 (관련된 k개의 결과)¶

import numpy as np
from scipy import stats

def cochran_q(data):
    """Cochran의 Q 검정.

    data : (피험자 수, 조건 수) 모양의 0/1 행렬.
           행 하나가 피험자, 열 하나가 조건이다.
    돌려주는 값은 (Q 통계량, chi2(k-1) 근사 p-값).
    """
    data = np.asarray(data, dtype=float)
    n, k = data.shape

    T_j = data.sum(axis=0)          # 조건별 성공 수
    L_i = data.sum(axis=1)          # 피험자별 성공 수
    grand_T = T_j.sum()

    # 분자는 조건 사이의 변동, 분모는 피험자 안의 변동을 잰다.
    # L_i가 0이거나 k인 피험자, 즉 전부 성공하거나 전부 실패한 피험자는
    # 분모에 아무 기여도 하지 않는다. 조건 간 비교에 정보가 없기 때문이다.
    numerator = (k - 1) * (k * np.sum(T_j**2) - grand_T**2)
    denominator = k * grand_T - np.sum(L_i**2)

    Q = numerator / denominator
    p_value = stats.chi2(k - 1).sf(Q)
    return Q, p_value

# 피험자 12명이 과제 3개를 수행한 결과 (성공=1, 실패=0)
tasks = np.array([
    [0, 1, 0],
    [1, 1, 0],
    [1, 1, 1],
    [0, 0, 0],
    [1, 0, 0],
    [0, 1, 1],
    [0, 0, 0],
    [1, 1, 0],
    [0, 1, 0],
    [0, 1, 0],
    [0, 1, 0],
    [0, 1, 0],
])

Q, p = cochran_q(tasks)

print(f"Cochran's Q = {Q:.4f}")
print(f"p-value     = {p:.4f}")

if p < 0.05:
    print("Reject H0: success rates differ across tasks.")
else:
    print("Fail to reject H0: no significant difference.")

# (1) 두 꼴이 같은지, 기여하는 사람이 몇 명인지
T_j, L_i = tasks.sum(axis=0), tasks.sum(axis=1)
k = tasks.shape[1]
print(f"\n분자  k(k-1)Σ(T_j-T̄)² = {k * (k - 1) * ((T_j - T_j.mean())**2).sum():.0f}")
print(f"분모  ΣL_i(k-L_i)      = {(L_i * (k - L_i)).sum():.0f}"
      f"   기여하는 사람 {np.sum((L_i > 0) & (L_i < k))}명")

# (2) 자유도 2 의 닫힌 꼴
print(f"exp(-Q/2) = {np.exp(-Q / 2):.7f}")

# (3) 쌍별 McNemar 와 본페로니
print()
for a, b_ in [(0, 1), (0, 2), (1, 2)]:
    b = int(((tasks[:, a] == 1) & (tasks[:, b_] == 0)).sum())
    c = int(((tasks[:, a] == 0) & (tasks[:, b_] == 1)).sum())
    chi = (b - c) ** 2 / (b + c)
    p_chi = stats.chi2(1).sf(chi)
    p_exact = min(1.0, 2 * stats.binom(b + c, 0.5).cdf(min(b, c)))
    print(f"과제 {a+1} 대 {b_+1}:  b={b} c={c}  chi2={chi:.4f}  "
          f"p={p_chi:.4f}  정확 p={p_exact:.4f}"
          f"{'  ← 유의 (0.0167)' if p_chi < 0.05 / 3 else ''}")

Cochran's Q = 8.6667
p-value     = 0.0131

분자  k(k-1)Σ(T_j-T̄)² = 156
분모  ΣL_i(k-L_i)      = 18   기여하는 사람 9명
exp(-Q/2) = 0.0131237

과제 1 대 2:  b=1 c=6  chi2=3.5714  p=0.0588  정확 p=0.1250
과제 1 대 3:  b=3 c=1  chi2=1.0000  p=0.3173  정확 p=0.6250
과제 2 대 3:  b=7 c=0  chi2=7.0000  p=0.0082  정확 p=0.0156  ← 유의 (0.0167)

import numpy as np
from scipy import stats

def cochran_q(data):
    data = np.asarray(data, float)
    n, k = data.shape
    T_j, L_i = data.sum(0), data.sum(1)
    grand_T = T_j.sum()
    Q = ((k - 1) * (k * np.sum(T_j**2) - grand_T**2)
         / (k * grand_T - np.sum(L_i**2)))
    return Q, stats.chi2(k - 1).sf(Q)

tasks = np.array([[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 0, 0],
                  [1, 0, 0], [0, 1, 1], [0, 0, 0], [1, 1, 0],
                  [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0]], float)
n, k = tasks.shape
L = tasks.sum(1)

Q, p = cochran_q(tasks)
print(f"전체 {n}명:  Q = {Q:.4f},  p = {p:.4f}")
print(f"  열 합 T_j = {tasks.sum(0).tolist()}")
print(f"  행 합 L_i 의 분포 = {np.bincount(L.astype(int)).tolist()}"
      f"   (L=0 이 {int((L == 0).sum())}명, L=k 가 {int((L == k).sum())}명)")

informative = tasks[(L > 0) & (L < k)]
Q2, p2 = cochran_q(informative)
print(f"\n정보 있는 {len(informative)}명만:  Q = {Q2:.4f},  p = {p2:.4f}")
print(f"  같은가? {np.isclose(Q, Q2)}")

전체 12명:  Q = 8.6667,  p = 0.0131
  열 합 T_j = [4.0, 9.0, 2.0]
  행 합 L_i 의 분포 = [2, 6, 3, 1]   (L=0 이 2명, L=k 가 1명)

정보 있는 9명만:  Q = 8.6667,  p = 0.0131
  같은가? True

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

tasks = np.array([[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 0, 0],
                  [1, 0, 0], [0, 1, 1], [0, 0, 0], [1, 1, 0],
                  [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0]], float)
print(f"과제별 성공률 {(tasks.mean(0)).round(4).tolist()}\n")

print("쌍별 McNemar (정확 이항)")
pvals, labels = [], []
for i, j in combinations(range(3), 2):
    b = int(((tasks[:, i] == 1) & (tasks[:, j] == 0)).sum())
    c = int(((tasks[:, i] == 0) & (tasks[:, j] == 1)).sum())
    p = min(1.0, 2 * stats.binom.cdf(min(b, c), b + c, 0.5)) \
        if b + c > 0 else 1.0
    pvals.append(p)
    labels.append(f"과제{i + 1}-과제{j + 1}")
    print(f"  과제{i + 1}-과제{j + 1}:  b = {b},  c = {c},  "
          f"불일치 {b + c}쌍,  p = {p:.4f}")

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

과제별 성공률 [0.3333, 0.75, 0.1667]

쌍별 McNemar (정확 이항)
  과제1-과제2:  b = 1,  c = 6,  불일치 7쌍,  p = 0.1250
  과제1-과제3:  b = 3,  c = 1,  불일치 4쌍,  p = 0.6250
  과제2-과제3:  b = 7,  c = 0,  불일치 7쌍,  p = 0.0156
  본페로니: 과제1-과제2 0.3750    과제1-과제3 1.0000    과제2-과제3 0.0469*
  홀름  : 과제1-과제2 0.2500    과제1-과제3 0.6250    과제2-과제3 0.0469*

import numpy as np
from scipy import stats

rng = np.random.default_rng(159)

def cochran_q_stat(data):
    data = np.asarray(data, float)
    n, k = data.shape
    T_j, L_i = data.sum(0), data.sum(1)
    grand_T = T_j.sum()
    den = k * grand_T - np.sum(L_i**2)
    if den == 0:
        return 0.0
    return (k - 1) * (k * np.sum(T_j**2) - grand_T**2) / den

def perm_pvalue(data, B, rng):
    """행 안에서 조건 이름표를 섞는 순열검정."""
    q_obs = cochran_q_stat(data)
    cnt = 1
    for _ in range(B):
        shuffled = np.array([rng.permutation(row) for row in data])
        cnt += cochran_q_stat(shuffled) >= q_obs - 1e-9
    return cnt / (B + 1)

tasks = np.array([[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 0, 0],
                  [1, 0, 0], [0, 1, 1], [0, 0, 0], [1, 1, 0],
                  [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0]], float)
q = cochran_q_stat(tasks)
print(f"보기 2:  Q = {q:.4f}")
print(f"  카이제곱 근사 p = {stats.chi2(2).sf(q):.4f}")
print(f"  순열 p          = {perm_pvalue(tasks, 9_999, rng):.4f}\n")

M = 3_000
print(f"{'n':>4s} {'k':>3s} {'p':>5s} {'카이제곱 근사':>14s} {'순열':>9s}")
for n, k, p_true in [(12, 3, 0.3), (20, 3, 0.3), (40, 3, 0.3), (12, 4, 0.5)]:
    a = b = 0
    for _ in range(M):
        d = (rng.random((n, k)) < p_true).astype(float)
        a += stats.chi2(k - 1).sf(cochran_q_stat(d)) < 0.05
        b += perm_pvalue(d, 299, rng) < 0.05
    print(f"{n:4d} {k:3d} {p_true:5.1f} {a / M:14.4f} {b / M:9.4f}")

보기 2:  Q = 8.6667
  카이제곱 근사 p = 0.0131
  순열 p          = 0.0152

   n   k     p        카이제곱 근사        순열
  12   3   0.3         0.0557    0.0353
  20   3   0.3         0.0433    0.0310
  40   3   0.3         0.0543    0.0457
  12   4   0.5         0.0457    0.0370

import numpy as np
from scipy import stats

rng = np.random.default_rng(753)
M = 3_000

def cochran_q_stat(data):
    data = np.asarray(data, float)
    n, k = data.shape
    T_j, L_i = data.sum(0), data.sum(1)
    grand_T = T_j.sum()
    den = k * grand_T - np.sum(L_i**2)
    return 0.0 if den == 0 else \
        (k - 1) * (k * np.sum(T_j**2) - grand_T**2) / den

print(f"{'n':>5s} {'조건별 성공률':>22s} {'정보 있는 피험자':>16s} {'검정력':>8s}")
settings = [(20, [0.3, 0.5, 0.7]), (40, [0.3, 0.5, 0.7]),
            (80, [0.3, 0.5, 0.7]),
            (40, [0.4, 0.5, 0.6]), (40, [0.05, 0.10, 0.15]),
            (40, [0.80, 0.85, 0.90])]
for n, ps in settings:
    ps = np.array(ps)
    rej, info = 0, []
    for _ in range(M):
        d = (rng.random((n, len(ps))) < ps).astype(float)
        L = d.sum(1)
        info.append(((L > 0) & (L < len(ps))).sum())
        rej += stats.chi2(len(ps) - 1).sf(cochran_q_stat(d)) < 0.05
    print(f"{n:5d} {str(ps.tolist()):>22s} {np.mean(info):16.1f} "
          f"{rej / M:8.4f}")

    n                조건별 성공률        정보 있는 피험자      검정력
   20        [0.3, 0.5, 0.7]             15.8   0.6077
   40        [0.3, 0.5, 0.7]             31.6   0.9183
   80        [0.3, 0.5, 0.7]             63.0   0.9983
   40        [0.4, 0.5, 0.6]             30.4   0.3520
   40      [0.05, 0.1, 0.15]             10.9   0.2203
   40       [0.8, 0.85, 0.9]             15.4   0.1810

같은 대상을 여러 조건에서 측정했다
    │
    ├─ 결과가 이진인가
    │      │
    │      ├─ 조건이 2개 ────→ McNemar
    │      └─ 조건이 k≥3 ───→ 코크런 Q
    │
    ├─ 결과가 k개 범주(순서 없음)
    │      └─ 조건이 2개 ────→ 보커 대칭성 / 주변 동질성
    │
    ├─ 결과가 순서형
    │      └─→ 대응 순위 검정, 순서형 혼합모형
    │
    └─ 공변량·결측·불균형이 있다
           └─→ 혼합효과 로지스틱 / GEE / 조건부 로지스틱