독립성 검정 - 타이타닉 생존과 성별¶

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd

URL = ("https://raw.githubusercontent.com/datasciencedojo/"
       "datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv")
df = pd.read_csv(URL, index_col="PassengerId")

# Survived 와 Sex 에는 결측이 없다. 891명 전수를 쓴다.
print(f"Survived 결측 {df['Survived'].isna().sum()}개, "
      f"Sex 결측 {df['Sex'].isna().sum()}개")

tab = pd.crosstab(df["Sex"], df["Survived"])
print("\n관측도수")
print(pd.crosstab(df["Sex"], df["Survived"], margins=True).to_string())

# 기대도수를 정의대로 만든다. scipy 를 쓰지 않는다.
n = tab.values.sum()
row = tab.sum(axis=1)        # 314, 577
col = tab.sum(axis=0)        # 549, 342
exp = pd.DataFrame(np.outer(row, col) / n,
                   index=tab.index, columns=tab.columns)

print("\n기대도수 (H0: 독립)")
print(exp.round(2).to_string())

# 한 칸만 손으로 확인해 본다.
print(f"\n손계산 E[여성, 생존] = 314 * 342 / 891 = "
      f"{row['female'] * col[1] / n:.2f}")

print("\n관측 - 기대")
print((tab - exp).round(2).to_string())

Survived 결측 0개, Sex 결측 0개

관측도수
Survived    0    1  All
Sex
female     81  233  314
male      468  109  577
All       549  342  891

기대도수 (H0: 독립)
Survived       0       1
Sex
female    193.47  120.53
male      355.53  221.47

손계산 E[여성, 생존] = 314 * 342 / 891 = 120.53

관측 - 기대
Survived       0       1
Sex
female   -112.47  112.47
male      112.47 -112.47

# 칸마다 (관측-기대)^2 / 기대 를 구해 더한다. 그것이 전부다.
cells = (tab - exp) ** 2 / exp
print("칸별 기여도")
print(cells.round(2).to_string())

chi2_stat = cells.values.sum()
print(f"\nchi2 = {chi2_stat:.4f}")

# 기여도가 큰 칸부터 본다. 어디서 독립이 깨졌는지 알려 준다.
flat = cells.stack().sort_values(ascending=False)
print("\n기여도가 큰 칸부터")
for (s, surv), v in flat.items():
    obs_v = tab.loc[s, surv]
    exp_v = exp.loc[s, surv]
    mark = "관측이 많다" if obs_v > exp_v else "관측이 적다"
    print(f"  {s:6s} 생존={surv}  기여 {v:6.2f}  "
          f"({obs_v:3d} 대 {exp_v:6.2f}, {mark})")

칸별 기여도
Survived      0       1
Sex
female    65.39  104.96
male      35.58   57.12

chi2 = 263.0506

기여도가 큰 칸부터
  female 생존=1  기여 104.96  (233 대 120.53, 관측이 많다)
  female 생존=0  기여  65.39  ( 81 대 193.47, 관측이 적다)
  male   생존=1  기여  57.12  (109 대 221.47, 관측이 적다)
  male   생존=0  기여  35.58  (468 대 355.53, 관측이 많다)

from scipy.stats import chi2_contingency
from scipy.stats import chi2 as chi2_dist

# scipy 의 기본값은 correction=True 이며, 2x2 표에만 예이츠 보정을 건다.
# 2x2 를 다룰 때 이 기본값을 모르면 손계산과 답이 어긋난다.
for flag in [False, True]:
    res = chi2_contingency(tab, correction=flag)
    label = "보정 없음" if not flag else "보정 있음(scipy 기본값)"
    print(f"{label:22s} chi2 = {res.statistic:8.4f}   "
          f"p = {res.pvalue:.4e}")

print(f"\n정의대로 구한 값        chi2 = {chi2_stat:8.4f}")
print(f"임계값 chi2_0.05(1)     {chi2_dist.ppf(0.95, 1):.4f}")
print(f"p값을 직접 구하면       {chi2_dist.sf(chi2_stat, 1):.4e}")

# 보정은 각 칸의 이탈을 0.5 만큼 줄인 뒤 계산한다.
yates = ((tab - exp).abs() - 0.5) ** 2 / exp
print(f"\n손으로 건 예이츠 보정   {yates.values.sum():.4f}")

보정 없음                  chi2 = 263.0506   p = 3.7117e-59
보정 있음(scipy 기본값)       chi2 = 260.7170   p = 1.1974e-58

정의대로 구한 값        chi2 = 263.0506
임계값 chi2_0.05(1)     3.8415
p값을 직접 구하면       3.7117e-59

손으로 건 예이츠 보정   260.7170

import numpy as np

rng = np.random.default_rng(0)

n, n_f, n_s = 891, 314, 342          # 전체, 여성, 생존자

def chi2_2x2(a):
    """여성 생존자 수 a 하나로 2x2 표가 결정된다 (주변합이 고정되어 있으므로)."""
    b = n_f - a                      # 여성 사망
    c = n_s - a                      # 남성 생존
    d = n - n_f - c                  # 남성 사망
    return n * (a * d - b * c) ** 2 / (n_f * (n - n_f) * n_s * (n - n_s))

obs = chi2_2x2(233)
print(f"관측 chi2 = {obs:.4f}  (여성 생존자 233명)")

# 생존자 342명을 성별과 무관하게 다시 고른다. 이것이 H0 가 참인 세계다.
B = 10000
is_female = np.zeros(n, dtype=bool)
is_female[:n_f] = True

a_perm = np.empty(B, dtype=int)
for i in range(B):
    survivors = rng.permutation(n)[:n_s]
    a_perm[i] = is_female[survivors].sum()

stats = chi2_2x2(a_perm)

print(f"\n재표본 {B}회")
print(f"  여성 생존자 수  관측 233명,  재표본 {a_perm.min()}~{a_perm.max()}명")
print(f"  chi2 최댓값     {stats.max():.4f}")
print(f"  chi2 평균 {stats.mean():.4f}, 분산 {stats.var():.4f}"
      f"   (chi2(1) 이론값 1, 2)")
print(f"  관측값 이상     {(stats >= obs).sum()}회")

관측 chi2 = 263.0506  (여성 생존자 233명)

재표본 10000회
  여성 생존자 수  관측 233명,  재표본 97~146명
  chi2 최댓값     13.4943
  chi2 평균 1.0140, 분산 1.9319   (chi2(1) 이론값 1, 2)
  관측값 이상     0회

import numpy as np
from scipy.stats import hypergeom

n, R1, R2, C1, C2 = 891, 314, 577, 549, 342     # 전체, 여성, 남성, 사망, 생존

A = np.arange(0, R1 + 1)             # 여성 생존자 수가 취할 수 있는 모든 값
pmf = hypergeom.pmf(A, n, C2, R1)    # 전체 n, 표시된 것 C2, 뽑는 수 R1

mu = R1 * C2 / n
var = R1 * R2 * C1 * C2 / (n ** 2 * (n - 1))
const = n ** 3 / (R1 * R2 * C1 * C2)
chis = const * (A - mu) ** 2         # 보기 2 의 닫힌 꼴을 A 로 적은 것

m1 = (pmf * chis).sum()                      # 평균
m2 = (pmf * (chis - m1) ** 2).sum()          # 분산
m4 = (pmf * (chis - m1) ** 4).sum()          # 네 번째 중심적률

print(f"확률의 총합      {pmf.sum():.12f}")
print(f"E[A] = {mu:.4f},  sd(A) = {np.sqrt(var):.4f}")
print(f"Var(A) * const = {var * const:.7f}   이론 n/(n-1) = {n / (n - 1):.7f}")
print(f"정확한 평균 {m1:.7f}   모의 1.0140   오차 {np.sqrt(m2 / 10000):.4f}"
      f"   차이/오차 {abs(1.0140 - m1) / np.sqrt(m2 / 10000):.2f}")
print(f"정확한 분산 {m2:.7f}   모의 1.9319   오차 "
      f"{np.sqrt((m4 - m2 ** 2) / 10000):.4f}   차이/오차 "
      f"{abs(1.9319 - m2) / np.sqrt((m4 - m2 ** 2) / 10000):.2f}")
print(f"E[Z^4] = {m2 / (n / (n - 1)) ** 2 + 1:.4f}   (정규분포라면 3)")

확률의 총합      1.000000000000
E[A] = 120.5253,  sd(A) = 6.9387
Var(A) * const = 1.0011236   이론 n/(n-1) = 1.0011236
정확한 평균 1.0011236   모의 1.0140   오차 0.0141   차이/오차 0.91
정확한 분산 2.0009729   모의 1.9319   오차 0.0747   차이/오차 0.93
E[Z^4] = 2.9965   (정규분포라면 3)

import numpy as np
from scipy.stats import chi2_contingency

# 타이타닉의 비율(여성 74.20%, 남성 18.89%, 여성 비율 35.24%)을 그대로 두고
# 전체 인원만 바꾼다. 효과의 크기는 조금도 변하지 않는다.
f_rate, m_rate, f_share = 0.7420, 0.1889, 0.3524

print(f"{'n':>7s}{'chi2':>11s}{'p':>12s}{'phi':>9s}"
      f"{'생존율 차':>11s}{'유의':>7s}")
for n_tot in [20, 50, 100, 300, 891, 5000]:
    nf = round(n_tot * f_share); nm = n_tot - nf
    fs = round(nf * f_rate);     ms = round(nm * m_rate)
    T = np.array([[nf - fs, fs], [nm - ms, ms]])
    res = chi2_contingency(T, correction=False)
    phi = np.sqrt(res.statistic / n_tot)
    diff = fs / nf - ms / nm
    print(f"{n_tot:>7d}{res.statistic:>11.3f}{res.pvalue:>12.2e}"
          f"{phi:>9.4f}{diff:>11.4f}"
          f"{'예' if res.pvalue < 0.05 else '아니오':>7s}")

      n       chi2           p      phi      생존율 차     유의
     20      6.282    1.22e-02   0.5604     0.5604      예
     50     13.981    1.85e-04   0.5288     0.5347      예
    100     30.092    4.12e-08   0.5486     0.5582      예
    300     88.890    4.17e-21   0.5443     0.5546      예
    891    263.051    3.71e-59   0.5434     0.5531      예
   5000   1474.237    0.00e+00   0.5430     0.5528      예

import numpy as np
from scipy.stats import chi2 as chi2_dist

crit = chi2_dist.ppf(0.95, 1)
print(f"chi2_0.05(1) = {crit:.4f}")
for phi, label in [(0.5430, "타이타닉"), (0.30, "중간"),
                   (0.10, "약함"), (0.0205, "40% 대 38%")]:
    print(f"  phi = {phi:.4f} ({label:10s})  n >= {crit / phi ** 2:9.1f}")

# 40% 대 38% 짜리 효과를 큰 표본으로 재면
phi = 0.020502
print()
for n_tot in [9000, 10000, 100000]:
    chi2_v = n_tot * phi ** 2
    print(f"  n = {n_tot:>6d}   chi2 = {chi2_v:7.3f}   "
          f"p = {chi2_dist.sf(chi2_v, 1):.3e}")

chi2_0.05(1) = 3.8415
  phi = 0.5430 (타이타닉      )  n >=      13.0
  phi = 0.3000 (중간        )  n >=      42.7
  phi = 0.1000 (약함        )  n >=     384.1
  phi = 0.0205 (40% 대 38% )  n >=    9140.9

  n =   9000   chi2 =   3.783   p = 5.178e-02
  n =  10000   chi2 =   4.203   p = 4.034e-02
  n = 100000   chi2 =  42.033   p = 8.974e-11

import numpy as np

a = tab.loc["female", 1]; b = tab.loc["female", 0]   # 여성 생존 / 사망
c = tab.loc["male",   1]; d = tab.loc["male",   0]   # 남성 생존 / 사망
n1, n2 = a + b, c + d
p1, p2 = a / n1, c / n2

print(f"여성 생존율 {p1:.4f} ({a}/{n1})")
print(f"남성 생존율 {p2:.4f} ({c}/{n2})")

print(f"\n비율 차이  {p1 - p2:.4f}   ({(p1 - p2) * 100:.2f}%포인트)")
print(f"위험비     {p1 / p2:.4f}배")
print(f"오즈비     {(a * d) / (b * c):.4f}")

# 파이 계수는 카이제곱을 n 으로 나눈 것의 제곱근이다.
# 2x2 에서는 Cramer 의 V 와 같은 값이 된다.
phi = np.sqrt(chi2_stat / n)
print(f"\nphi = sqrt(chi2/n) = sqrt({chi2_stat:.4f}/{n}) = {phi:.4f}")

# 주변 비율이 다르면 phi 가 1에 닿지 못한다. 도달 가능한 최댓값을 구해 둔다.
p_s = col[1] / n                 # 생존 비율
q_f = row["female"] / n          # 여성 비율
phi_max = (min(p_s, q_f) - p_s * q_f) / np.sqrt(
    p_s * (1 - p_s) * q_f * (1 - q_f))
print(f"\n생존 비율 {p_s:.4f}, 여성 비율 {q_f:.4f}")
print(f"phi_max = {phi_max:.4f},  도달률 = {phi / phi_max:.4f}")

여성 생존율 0.7420 (233/314)
남성 생존율 0.1889 (109/577)

비율 차이  0.5531   (55.31%포인트)
위험비     3.9280배
오즈비     12.3507

phi = sqrt(chi2/n) = sqrt(263.0506/891) = 0.5434

생존 비율 0.3838, 여성 비율 0.3524
phi_max = 0.9347,  도달률 = 0.5813

import numpy as np
from scipy.stats import norm

z = norm.ppf(0.975)

# (1) 비율 차이 -- 원 척도에서 대칭 구간을 쓴다.
se_d = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
diff = p1 - p2
print(f"비율 차이 {diff:.4f},  SE {se_d:.4f}")
print(f"  95% CI [{diff - z * se_d:.4f}, {diff + z * se_d:.4f}]")

# (2) 위험비, (3) 오즈비 -- 로그 척도에서 구간을 만든 뒤 되돌린다.
# 두 비율 모두 표본분포가 오른쪽으로 치우쳐 있어 원 척도에서는 대칭이 아니다.
se_lr = np.sqrt(1 / a - 1 / n1 + 1 / c - 1 / n2)
lr = np.log(p1 / p2)
print(f"\n위험비 {np.exp(lr):.4f}")
print(f"  95% CI [{np.exp(lr - z * se_lr):.4f}, "
      f"{np.exp(lr + z * se_lr):.4f}]")

se_lo = np.sqrt(1 / a + 1 / b + 1 / c + 1 / d)
lo = np.log((a * d) / (b * c))
print(f"\n오즈비 {np.exp(lo):.4f},  log(OR) {lo:.4f},  SE {se_lo:.4f}")
print(f"  95% CI [{np.exp(lo - z * se_lo):.4f}, "
      f"{np.exp(lo + z * se_lo):.4f}]")
print(f"  로그 척도에서 [{lo - z * se_lo:.4f}, {lo + z * se_lo:.4f}] "
      f"를 지수로 되돌린 것이다")

비율 차이 0.5531,  SE 0.0296
  95% CI [0.4951, 0.6111]

위험비 3.9280
  95% CI [3.2770, 4.7084]

오즈비 12.3507,  log(OR) 2.5137,  SE 0.1672
  95% CI [8.9000, 17.1393]
  로그 척도에서 [2.1860, 2.8414] 를 지수로 되돌린 것이다

from scipy.stats import chi2 as chi2_dist

chi2_stat = 263.0506
for alpha in [0.05, 0.01, 0.001]:
    crit = chi2_dist.ppf(1 - alpha, 1)
    print(f"alpha={alpha:<6} 임계값 {crit:7.4f}   "
          f"관측 {chi2_stat:.2f}   "
          f"{'기각' if chi2_stat > crit else '기각 못 함'}")

alpha=0.05   임계값  3.8415   관측 263.05   기각
alpha=0.01   임계값  6.6349   관측 263.05   기각
alpha=0.001  임계값 10.8276   관측 263.05   기각

import numpy as np
import pandas as pd
from scipy.stats import chi2_contingency

# 타이타닉의 비율을 유지한 채 표본만 줄인다.
f_rate, m_rate, f_share = 0.7420, 0.1889, 0.3524

print(f"{'n':>6s}{'최소 기대도수':>14s}{'보정 없음':>11s}"
      f"{'보정 있음':>11s}{'차이(%)':>10s}")
for n_tot in [20, 40, 100, 300, 891]:
    nf = round(n_tot * f_share); nm = n_tot - nf
    fs = round(nf * f_rate);     ms = round(nm * m_rate)
    T = np.array([[nf - fs, fs], [nm - ms, ms]])
    r0 = chi2_contingency(T, correction=False)
    r1 = chi2_contingency(T, correction=True)
    print(f"{n_tot:>6d}{r0.expected_freq.min():>14.2f}"
          f"{r0.statistic:>11.4f}{r1.statistic:>11.4f}"
          f"{100 * (1 - r1.statistic / r0.statistic):>10.2f}")

     n       최소 기대도수      보정 없음      보정 있음     차이(%)
    20          2.45     6.2819     4.0599     35.37
    40          5.25    10.5788     8.4689     19.94
   100         13.30    30.0920    27.7692      7.72
   300         40.99    88.8899    86.5669      2.61
   891        120.53   263.0506   260.7170      0.89

import numpy as np
from scipy.stats import chi2_contingency

a, b = tab.loc["female", 0], tab.loc["female", 1]
c, d = tab.loc["male",   0], tab.loc["male",   1]
n_ = a + b + c + d

denom = (a + b) * (c + d) * (a + c) * (b + d)
phi_f  = (a * d - b * c) / np.sqrt(denom)
chi2_f = n_ * (a * d - b * c) ** 2 / denom

print(f"공식으로 구한 phi   {phi_f:+.6f}")
print(f"공식으로 구한 chi2  {chi2_f:.6f}")
print(f"scipy chi2          "
      f"{chi2_contingency(tab, correction=False).statistic:.6f}")
print(f"n * phi^2           {n_ * phi_f ** 2:.6f}")

공식으로 구한 phi   -0.543351
공식으로 구한 chi2  263.050574
scipy chi2          263.050574
n * phi^2           263.050574

import numpy as np

def phi_max(p, q):
    return (min(p, q) - p * q) / np.sqrt(p * (1 - p) * q * (1 - q))

print(f"{'p':>6s}{'q':>6s}{'phi_max':>10s}")
for pp, qq in [(0.50, 0.50), (0.38, 0.35), (0.50, 0.35),
               (0.50, 0.10), (0.10, 0.90), (0.05, 0.50)]:
    print(f"{pp:>6.2f}{qq:>6.2f}{phi_max(pp, qq):>10.4f}")

print(f"\n타이타닉 p=0.3838, q=0.3524 -> "
      f"phi_max = {phi_max(0.3838, 0.3524):.4f}")
print(f"관측 |phi| = 0.5434,  도달률 = "
      f"{0.5434 / phi_max(0.3838, 0.3524):.4f}")

# 연관의 '세기'(오즈비)는 같은데 주변 비율만 다른 두 표
print("\n오즈비는 같고 주변 비율만 다를 때")
print(f"{'표':>16s}{'phi':>9s}{'OR':>9s}")
for lab, T in [("균형 (50/50)", np.array([[40, 10], [10, 40]])),
               ("한쪽이 드묾",  np.array([[80, 20], [ 2,  8]]))]:
    a, b, c, d = T[0, 0], T[0, 1], T[1, 0], T[1, 1]
    phi = (a * d - b * c) / np.sqrt(
        (a + b) * (c + d) * (a + c) * (b + d))
    print(f"{lab:>16s}{phi:>9.4f}{a * d / (b * c):>9.4f}")

     p     q   phi_max
  0.50  0.50    1.0000
  0.38  0.35    0.9373
  0.50  0.35    0.7338
  0.50  0.10    0.3333
  0.10  0.90    0.1111
  0.05  0.50    0.2294

타이타닉 p=0.3838, q=0.3524 -> phi_max = 0.9347
관측 |phi| = 0.5434,  도달률 = 0.5814

오즈비는 같고 주변 비율만 다를 때
               표      phi       OR
      균형 (50/50)   0.6000  16.0000
          한쪽이 드묾   0.3960  16.0000

import numpy as np
import pandas as pd
from scipy.stats import chi2_contingency

print("전체 (등급 무시)")
fr = 233 / 314; mr = 109 / 577
print(f"  여성 {fr:.4f}  남성 {mr:.4f}  "
      f"차이 {fr - mr:+.4f}  OR 12.3507  chi2 263.05")

print("\n등급별")
print(f"{'등급':>5s}{'n':>6s}{'여성 생존율':>12s}{'남성 생존율':>12s}"
      f"{'차이':>9s}{'오즈비':>10s}{'chi2':>10s}{'p':>12s}")
for pc in [1, 2, 3]:
    s = df[df["Pclass"] == pc]
    t = pd.crosstab(s["Sex"], s["Survived"])
    f_ = t.loc["female", 1] / t.loc["female"].sum()
    m_ = t.loc["male",   1] / t.loc["male"].sum()
    orr = (t.loc["female", 1] * t.loc["male", 0]) / \
          (t.loc["female", 0] * t.loc["male", 1])
    r = chi2_contingency(t, correction=False)
    print(f"{pc:>5d}{len(s):>6d}{f_:>12.4f}{m_:>12.4f}"
          f"{f_ - m_:>+9.4f}{orr:>10.4f}"
          f"{r.statistic:>10.4f}{r.pvalue:>12.2e}")

print("\n등급 자체의 효과 (성별 무시)")
print(f"{'등급':>5s}{'생존율':>10s}{'여성 비율':>11s}")
for pc in [1, 2, 3]:
    s = df[df["Pclass"] == pc]
    print(f"{pc:>5d}{s['Survived'].mean():>10.4f}"
          f"{(s['Sex'] == 'female').mean():>11.4f}")

전체 (등급 무시)
  여성 0.7420  남성 0.1889  차이 +0.5531  OR 12.3507  chi2 263.05

등급별
   등급     n      여성 생존율      남성 생존율       차이       오즈비      chi2           p
    1   216      0.9681      0.3689  +0.5992   51.9037   81.7530    1.54e-19
    2   184      0.9211      0.1574  +0.7636   62.4510  104.3632    1.68e-24
    3   491      0.5000      0.1354  +0.3646    6.3830   73.6556    9.30e-18

등급 자체의 효과 (성별 무시)
   등급       생존율      여성 비율
    1    0.6296     0.4352
    2    0.4728     0.4130
    3    0.2424     0.2933

3등실:  생존율 0.2424 (가장 낮음),  여성 비율 0.2933 (가장 낮음)
1등실:  생존율 0.6296 (가장 높음),  여성 비율 0.4352 (가장 높음)

from scipy.stats import fisher_exact, chi2_contingency

# 행·열 순서를 [생존, 사망] x [여성, 남성] 으로 맞춰 오즈비 방향을 정한다.
T = [[233, 81],       # 여성 생존, 여성 사망
     [109, 468]]      # 남성 생존, 남성 사망
res_f = fisher_exact(T)
res_c = chi2_contingency(T, correction=False)

print(f"피셔 정확검정   OR = {res_f[0]:.4f}   p = {res_f[1]:.4e}")
print(f"카이제곱 검정   chi2 = {res_c.statistic:.4f}   "
      f"p = {res_c.pvalue:.4e}")
print(f"\n최소 기대도수 {res_c.expected_freq.min():.2f}")

피셔 정확검정   OR = 12.3507   p = 6.4639e-60
카이제곱 검정   chi2 = 263.0506   p = 3.7117e-59

최소 기대도수 120.53

import numpy as np, pandas as pd
from scipy import stats

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
t = pd.read_csv(URL)
ct = pd.crosstab(t.Sex, t.Survived)
O = ct.values.astype(float)
chi2, p, dof, E = stats.chi2_contingency(O, correction=False)

r = (O - E) / np.sqrt(E)                       # 피어슨 잔차
n = O.sum()
rowp, colp = O.sum(1) / n, O.sum(0) / n
adj = (O - E) / np.sqrt(E * np.outer(1 - rowp, 1 - colp))   # 조정잔차

print(f"{'칸':>16}{'관측':>8}{'기대':>10}{'피어슨잔차':>12}{'조정잔차':>11}{'기여%':>9}")
for i, s in enumerate(ct.index):
    for j, c in enumerate(ct.columns):
        lab = f"{s},{'생존' if c == 1 else '사망'}"
        print(f"{lab:>16}{O[i,j]:>8.0f}{E[i,j]:>10.2f}"
              f"{r[i,j]:>12.3f}{adj[i,j]:>11.3f}{r[i,j]**2/chi2*100:>9.1f}")

               칸      관측        기대       피어슨잔차       조정잔차      기여%
       female,사망      81    193.47      -8.086    -16.219     24.9
       female,생존     233    120.53      10.245     16.219     39.9
         male,사망     468    355.53       5.965     16.219     13.5
         male,생존     109    221.47      -7.558    -16.219     21.7

import numpy as np, pandas as pd
from scipy import stats

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
t = pd.read_csv(URL)

num = den = 0.0
numz = denz = 0.0
print(f"{'등급':>5}{'여생':>6}{'여사':>6}{'남생':>6}{'남사':>6}{'n':>6}{'층별 OR':>10}")
for pc in (1, 2, 3):
    s = t[t.Pclass == pc]
    c = pd.crosstab(s.Sex, s.Survived)
    a, b = c.loc['female', 1], c.loc['female', 0]
    cc, d = c.loc['male', 1],   c.loc['male', 0]
    n = a + b + cc + d
    num += a * d / n
    den += b * cc / n
    r1, c1 = a + b, a + cc                       # 행합, 열합
    numz += a - r1 * c1 / n
    denz += r1 * (n - r1) * c1 * (n - c1) / (n ** 2 * (n - 1))
    print(f"{pc:>5}{a:>6}{b:>6}{cc:>6}{d:>6}{n:>6}{(a*d)/(b*cc):>10.4f}")

ct = pd.crosstab(t.Sex, t.Survived)
crude = (ct.loc['female',1] * ct.loc['male',0]) / (ct.loc['female',0] * ct.loc['male',1])
mh = numz ** 2 / denz
print(f"\n전체(조) OR   = {crude:.4f}")
print(f"맨텔-헨첼 OR  = {num/den:.4f}")
print(f"맨텔-헨첼 chi2 = {mh:.4f}   p = {stats.chi2.sf(mh, 1):.3e}")

   등급    여생    여사    남생    남사     n     층별 OR
    1    91     3    45    77   216   51.9037
    2    70     6    17    91   184   62.4510
    3    72    72    47   300   491    6.3830

전체(조) OR   = 12.3507
맨텔-헨첼 OR  = 13.7586
맨텔-헨첼 chi2 = 250.4473   p = 2.075e-56

import numpy as np, pandas as pd
from scipy import stats

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
t = pd.read_csv(URL)
ct3 = pd.crosstab(t.Pclass, t.Survived)
c3, p3, d3, E3 = stats.chi2_contingency(ct3.values, correction=False)
print(ct3)

V = np.sqrt(c3 / (len(t) * min(ct3.shape[0] - 1, ct3.shape[1] - 1)))
print(f"\nchi2={c3:.4f}  dof={d3}  p={p3:.3e}   Cramer V={V:.4f}")

r3 = (ct3.values - E3) / np.sqrt(E3)
print("\n표준화잔차:")
for i, s in enumerate(ct3.index):
    print(f"  {s}등석: 사망 {r3[i,0]:+.3f}  생존 {r3[i,1]:+.3f}")

Survived    0    1
Pclass
1          80  136
2          97   87
3         372  119

chi2=102.8890  dof=2  p=4.549e-23   Cramer V=0.3398

표준화잔차:
  1등석: 사망 -4.602  생존 +5.831
  2등석: 사망 -1.538  생존 +1.948
  3등석: 사망 +3.994  생존 -5.060

import numpy as np
from scipy import stats

crit = stats.chi2.ppf(0.95, 1)
print(f"{'n':>6}{'w=0.1':>10}{'w=0.3':>10}{'w=0.5':>10}")
for nn in (50, 100, 200, 500, 891):
    row = [stats.ncx2.sf(crit, 1, nn * w ** 2) for w in (0.1, 0.3, 0.5)]
    print(f"{nn:>6}{row[0]:>10.4f}{row[1]:>10.4f}{row[2]:>10.4f}")

print()
for w in (0.1, 0.3, 0.5):
    nn = 1
    while stats.ncx2.sf(crit, 1, nn * w ** 2) < 0.8:
        nn += 1
    print(f"  검정력 0.8 에 필요한 n (w={w}): {nn}")

     n     w=0.1     w=0.3     w=0.5
    50    0.1090    0.5641    0.9424
   100    0.1701    0.8508    0.9988
   200    0.2930    0.9888    1.0000
   500    0.6088    1.0000    1.0000
   891    0.8473    1.0000    1.0000

  검정력 0.8 에 필요한 n (w=0.1): 785
  검정력 0.8 에 필요한 n (w=0.3): 88
  검정력 0.8 에 필요한 n (w=0.5): 32