이표본 및 다집단 검정¶

from scipy import stats

# 독립인 두 집단
group_a = [12, 15, 18, 22, 25]
group_b = [8, 10, 14, 19, 21, 24]

# 순위합 검정. 동점 보정을 하지 않으므로 동점이 많으면 mannwhitneyu 를 쓴다.
stat, p = stats.ranksums(group_a, group_b, alternative="two-sided")
print(f"Z = {stat:.4f}, p = {p:.2%}")
# Z = 0.7303, p = 46.52%

Z = 0.7303, p = 46.52%

import itertools
import numpy as np
from scipy import stats

m, n = len(group_a), len(group_b)
N = m + n
ranks = stats.rankdata(group_a + group_b)
W = ranks[:m].sum()
E = m * (N + 1) / 2
V = m * n * (N + 1) / 12
print(f"A 의 순위 {list(ranks[:m])}  ->  W = {W}")
print(f"E[W] = {E},  Var(W) = {V},  Z = {(W - E) / np.sqrt(V):.6f}")

# 귀무분포를 전수 열거한다. C(11, 5) = 462 가지뿐이다.
Ws = np.array([sum(c) for c in itertools.combinations(range(1, N + 1), m)])
print(f"\n열거한 가지 수 {len(Ws)},  평균 {Ws.mean()},  분산 {Ws.var()}")
hit = (np.abs(Ws - E) >= abs(W - E)).sum()
print(f"|W - {E:.0f}| >= {abs(W - E):.0f} 인 것 {hit} 개")
print(f"정확      양측 p = {hit / len(Ws):.6f}")
print(f"정규근사  양측 p = {2 * stats.norm.sf(abs(W - E) / np.sqrt(V)):.6f}")
print(f"연속성보정 양측 p = {2 * stats.norm.sf((abs(W - E) - 0.5) / np.sqrt(V)):.6f}")
print(f"mannwhitneyu(method='exact'): {stats.mannwhitneyu(group_a, group_b, method='exact')}")
print(f"U = W - m(m+1)/2 = {W - m * (m + 1) / 2}")

A 의 순위 [3.0, 5.0, 6.0, 9.0, 11.0]  ->  W = 34.0
E[W] = 30.0,  Var(W) = 30.0,  Z = 0.730297

열거한 가지 수 462,  평균 30.0,  분산 30.0
|W - 30| >= 4 인 것 248 개
정확      양측 p = 0.536797
정규근사  양측 p = 0.465209
연속성보정 양측 p = 0.522817
mannwhitneyu(method='exact'): MannwhitneyuResult(statistic=19.0, pvalue=0.5367965367965368)
U = W - m(m+1)/2 = 19.0

from scipy import stats
import numpy as np

paired_data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]
])

# 대응표본에는 부호순위검정을 쓴다. 짝을 무시하고 독립표본 검정을 쓰면
# 개체 간 변동이 그대로 남아 검정력을 크게 잃는다.
stat, p = stats.wilcoxon(
    paired_data[:, 0], paired_data[:, 1],
    alternative="two-sided", method="approx", zero_method="pratt"
)
print(f"W = {stat}, p = {p:.4f}")
# W = 11.0, p = 0.0086

W = 11.0, p = 0.0086

from collections import Counter

d = paired_data[:, 0] - paired_data[:, 1]
n, n0 = len(d), int((d == 0).sum())
r = stats.rankdata(np.abs(d))          # pratt: 0 도 함께 순위를 매긴다
print(f"차이  {list(d)}")
print(f"순위  {list(r)}   (0 이 {n0} 개)")
W_plus, W_minus = r[d > 0].sum(), r[d < 0].sum()
zero_rank_sum = r[d == 0].sum()
print(f"W+ = {W_plus},  W- = {W_minus},  합 = {W_plus + W_minus}"
      f"   (n(n+1)/2 - 0 의 순위합 = {n * (n + 1) / 2} - {zero_rank_sum})")

E = n * (n + 1) / 4 - n0 * (n0 + 1) / 4
tie = sum(c ** 3 - c for c in Counter(r[d != 0]).values())
V24 = n * (n + 1) * (2 * n + 1) - n0 * (n0 + 1) * (2 * n0 + 1) - tie / 2
V = V24 / 24
z = (W_minus - E) / np.sqrt(V)
print(f"\nE[W] = {E},  24Var = {n * (n + 1) * (2 * n + 1)} - "
      f"{n0 * (n0 + 1) * (2 * n0 + 1)} - {tie / 2} = {V24}")
print(f"Var = {V},  sigma = {np.sqrt(V):.6f}")
print(f"z = {z:.6f},  양측 p = {2 * stats.norm.cdf(z):.12f}")
print(f"scipy (pratt)  = {p:.12f}")

w = stats.wilcoxon(paired_data[:, 0], paired_data[:, 1],
                   method="approx", zero_method="wilcox")
print(f"scipy (wilcox) = W {w.statistic},  p = {w.pvalue:.12f}")

차이  [17, -2, 6, -3, 14, 0, 6, 10, 12, 8, 9, 2, 0, 3, 0]
순위  [15.0, 4.5, 8.5, 6.5, 14.0, 2.0, 8.5, 12.0, 13.0, 10.0, 11.0, 4.5, 2.0, 6.5, 2.0]   (0 이 3 개)
W+ = 103.0,  W- = 11.0,  합 = 114.0   (n(n+1)/2 - 0 의 순위합 = 120.0 - 6.0)

E[W] = 57.0,  24Var = 7440 - 84 - 9.0 = 7347.0
Var = 306.125,  sigma = 17.496428
z = -2.629108,  양측 p = 0.008560915879
scipy (pratt)  = 0.008560915879
scipy (wilcox) = W 5.0,  p = 0.007579201614

from scipy import stats

data0 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
data1 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]

# 크기가 다른 두 집단이어도 상관없다. 순위만 쓰기 때문이다.
stat, p = stats.mannwhitneyu(data0, data1)
print(f"U = {stat}, p = {p:.2%}")
# U = 49.0, p = 0.53%

U = 49.0, p = 0.53%

import numpy as np
from collections import Counter

m, n = len(data0), len(data1)
N = m + n
combined = data0 + data1
W = stats.rankdata(combined)[:m].sum()
U = W - m * (m + 1) / 2
mu = m * n / 2
counts = Counter(combined)
tie = sum(t ** 3 - t for t in counts.values())
print(f"W = {W},  U = W - m(m+1)/2 = {U},  E[U] = {mu}")
print(f"동점 묶음 {dict((k, v) for k, v in sorted(counts.items()) if v > 1)}"
      f"  ->  sum(t^3-t) = {tie}")

V_tie = m * n / 12 * ((N + 1) - tie / (N * (N - 1)))
V_plain = m * n * (N + 1) / 12
print(f"Var(U) 동점보정 = {V_tie:.6f} (sigma {np.sqrt(V_tie):.6f}),"
      f"  보정없음 = {V_plain} (sigma {np.sqrt(V_plain):.6f})")
print(f"동점항 / (N+1) = {tie / (N * (N - 1)) / (N + 1):.6%}")

for label, V, cc in (("동점 + 연속성", V_tie, 0.5),
                     ("동점만        ", V_tie, 0.0),
                     ("연속성만      ", V_plain, 0.5)):
    z = (U - mu + cc) / np.sqrt(V)
    print(f"  {label}: z = {z:9.6f},  p = {2 * stats.norm.cdf(z):.7f}")
print(f"  scipy 기본값             : p = {p:.7f}")
print(f"  scipy use_continuity=False: "
      f"p = {stats.mannwhitneyu(data0, data1, use_continuity=False).pvalue:.7f}")

W = 185.0,  U = W - m(m+1)/2 = 49.0,  E[U] = 120.0
동점 묶음 {14: 2, 31: 3, 43: 2}  ->  sum(t^3-t) = 36
Var(U) 동점보정 = 639.225806 (sigma 25.282915),  보정없음 = 640.0 (sigma 25.298221)
동점항 / (N+1) = 0.120968%
  동점 + 연속성: z = -2.788444,  p = 0.0052962
  동점만        : z = -2.808220,  p = 0.0049816
  연속성만      : z = -2.786757,  p = 0.0053238
  scipy 기본값             : p = 0.0052962
  scipy use_continuity=False: p = 0.0049816

from scipy import stats

data0 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
data1 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]
data2 = [0, 3, 9, 22, 23, 25, 25, 33, 34, 34, 40, 45, 46, 48, 62, 67, 84]

# 집단이 셋이 되면 Kruskal-Wallis 로 간다. 유의하다는 결론이 나면 어느 쌍이
# 다른지는 따로 사후비교를 해야 한다.
stat, p = stats.kruskal(data0, data1, data2)
print(f"H = {stat:.4f}, p = {p:.2%}")
# H = 7.6480, p = 2.18%

H = 7.6480, p = 2.18%

import numpy as np
from collections import Counter

groups = [data0, data1, data2]
sizes = [len(g) for g in groups]
allv = sum(groups, [])
N = len(allv)
ranks = stats.rankdata(allv)
edges = np.cumsum([0] + sizes)
R = [ranks[edges[i]:edges[i + 1]].sum() for i in range(3)]
print(f"n = {sizes},  N = {N}")
print(f"R = {R},  합 = {sum(R)}  (N(N+1)/2 = {N * (N + 1) / 2})")
print(f"평균순위 = {[round(R[i] / sizes[i], 4) for i in range(3)]}"
      f"   (전체 {(N + 1) / 2})")

H = 12 / (N * (N + 1)) * sum(R[i] ** 2 / sizes[i] for i in range(3)) - 3 * (N + 1)
counts = Counter(allv)
tie = sum(t ** 3 - t for t in counts.values())
C = 1 - tie / (N ** 3 - N)
print(f"\n본문 식의 H          = {H:.6f}")
print(f"동점 묶음 {dict((k, v) for k, v in sorted(counts.items()) if v > 1)}")
print(f"sum(t^3-t) = {tie},  C = 1 - {tie}/{N ** 3 - N} = {C:.6f}")
print(f"동점 보정한 H = H/C  = {H / C:.6f}")
print(f"scipy 가 준 H        = {stat:.6f}")
print(f"p: 보정 전 {stats.chi2.sf(H, 2):.6f}  ->  보정 후 {stats.chi2.sf(H / C, 2):.6f}"
      f"  (scipy {p:.6f})")

n = [16, 15, 17],  N = 48
R = [295.0, 483.0, 398.0],  합 = 1176.0  (N(N+1)/2 = 1176.0)
평균순위 = [18.4375, 32.2, 23.4118]   (전체 24.5)

본문 식의 H          = 7.640535
동점 묶음 {14: 2, 22: 2, 25: 3, 31: 3, 33: 2, 34: 3, 40: 2, 43: 2, 48: 2}
sum(t^3-t) = 108,  C = 1 - 108/110544 = 0.999023
동점 보정한 H = H/C  = 7.648007
scipy 가 준 H        = 7.648007
p: 보정 전 0.021922  ->  보정 후 0.021840  (scipy 0.021840)

from scipy import stats

data0 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
data1 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]
data2 = [0, 3, 9, 22, 23, 25, 25, 33, 34, 34, 40, 45, 46, 48, 62, 67, 84]

# 중앙값 검정은 전체 중앙값을 기준으로 위아래 개수만 센다. 정보를 많이
# 버리므로 위 Kruskal-Wallis 보다 p-값이 크게 나온다. 같은 자료에서
# 한쪽은 유의하고 한쪽은 그렇지 않은 것이 그 차이다.
result = stats.median_test(data0, data1, data2)
print(f"Grand median = {result.median}")
print(f"Contingency table:\n{result.table}")
print(f"p-value = {result.pvalue:.4f}")
# Grand median = 34.0
# Contingency table:
# [[ 5 10  7]
#  [11  5 10]]
# p-value = 0.1261

Grand median = 34.0
Contingency table:
[[ 5 10  7]
 [11  5 10]]
p-value = 0.1261

import numpy as np

allv = np.array(data0 + data1 + data2)
s = np.sort(allv)
N = len(allv)
print(f"N = {N},  24번째 = {s[23]},  25번째 = {s[24]},  전체 중앙값 = {np.median(allv)}")
print(f"> 34 인 것 {(allv > 34).sum()} 개,  < 34 인 것 {(allv < 34).sum()} 개,"
      f"  = 34 인 것 {(allv == 34).sum()} 개")

T = np.array([[5, 10, 7], [11, 5, 10]])
row, col = T.sum(1), T.sum(0)
E = np.outer(row, col) / T.sum()
print(f"\n행합 {list(row)},  열합 {list(col)}")
print("기대도수\n", np.round(E, 4))
chi2 = ((T - E) ** 2 / E).sum()
print(f"칸별 기여 {np.round((T - E) ** 2 / E, 6).tolist()}")
print(f"chi2 = {chi2:.6f},  p = {stats.chi2.sf(chi2, 2):.6f}"
      f"  (chi2_2 의 5% 임계값 {stats.chi2.ppf(0.95, 2):.4f})")

print("\nties 인수를 바꾸어 본다")
for tie in ("below", "above", "ignore"):
    r = stats.median_test(data0, data1, data2, ties=tie)
    print(f"  ties={tie:7s} {r.table.tolist()}  "
          f"chi2 = {r.statistic:.6f},  p = {r.pvalue:.6f}")
for i, g in enumerate([data0, data1, data2]):
    print(f"  집단 {i} 에 34 가 {np.sum(np.array(g) == 34)} 개")

N = 48,  24번째 = 34,  25번째 = 34,  전체 중앙값 = 34.0
> 34 인 것 22 개,  < 34 인 것 23 개,  = 34 인 것 3 개

행합 [22, 26],  열합 [16, 15, 17]
기대도수
 [[7.3333 6.875  7.7917]
 [8.6667 8.125  9.2083]]
칸별 기여 [[0.742424, 1.420455, 0.080437], [0.628205, 1.201923, 0.068062]]
chi2 = 4.141506,  p = 0.126091  (chi2_2 의 5% 임계값 5.9915)

ties 인수를 바꾸어 본다
  ties=below   [[5, 10, 7], [11, 5, 10]]  chi2 = 4.141506,  p = 0.126091
  ties=above   [[5, 11, 9], [11, 4, 8]]  chi2 = 5.501708,  p = 0.063873
  ties=ignore  [[5, 10, 7], [11, 4, 8]]  chi2 = 4.868277,  p = 0.087673
  집단 0 에 34 가 0 개
  집단 1 에 34 가 1 개
  집단 2 에 34 가 2 개

from scipy import stats
a = [72, 78, 81, 85, 90]; b = [68, 74, 77, 83, 88, 92]
print(stats.ranksums(a, b))
# RanksumsResult(statistic=0.18257, pvalue=0.85513)
print(stats.mannwhitneyu(a, b, method='exact'))
# MannwhitneyuResult(statistic=16.0, pvalue=0.93074)

RanksumsResult(statistic=0.18257418583505536, pvalue=0.8551321405847059)
MannwhitneyuResult(statistic=16.0, pvalue=0.9307359307359306)

from scipy import stats
t1 = [14, 18, 22, 25]; t2 = [19, 23, 27, 30, 35]; t3 = [10, 15, 20]
print(stats.kruskal(t1, t2, t3))
# KruskalResult(statistic=5.42564, pvalue=0.06637)

KruskalResult(statistic=5.425641025641035, pvalue=0.06634940320987275)

import numpy as np
from scipy import stats

data0 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
data1 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]

combined = np.array(data0 + data1)
ranks = stats.rankdata(combined)
W = ranks[:16].sum()
print(f"W = {W}")                        # W = 185.0
print(f"W - m(m+1)/2 = {W - 16*17/2}")   # 49.0

U, p = stats.mannwhitneyu(data0, data1)
print(f"U = {U}")                        # U = 49.0

W = 185.0
W - m(m+1)/2 = 49.0
U = 49.0

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
x = rng.normal(0, 1, 12); y = rng.normal(0.5, 1, 15)
H = stats.kruskal(x, y).statistic
Z = stats.ranksums(x, y).statistic
print(H, Z ** 2, abs(H - Z ** 2) < 1e-10)
# True

0.7714285714285722 0.7714285714285715 True