사례 연구: 절단점 하나가 결론을 바꾼다¶

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")

# 나이가 없는 177명은 이 절의 분석에서 아예 쓸 수 없다.
# 앞 절 연습문제 3에서 보았듯 이들은 무작위로 빠진 것이 아니다.
d = df.dropna(subset=["Age"]).copy()
print(f"전체 {len(df)}명 중 나이가 있는 {len(d)}명으로 분석한다")
print(f"나이 범위 {d['Age'].min():.2f} ~ {d['Age'].max():.2f}세\n")

CUT = 21                                   # ← 이 한 줄이 이 절의 주제다
d["Age_Group"] = np.where(d["Age"] < CUT, f"<{CUT}", f">={CUT}")
ORDER = [f"<{CUT}", f">={CUT}"]

tab = pd.crosstab(d["Age_Group"], d["Survived"]).reindex(ORDER)
print("도수")
print(tab.to_string())

pct = pd.crosstab(d["Age_Group"], d["Survived"],
                  normalize="index").reindex(ORDER) * 100
print("\n집단 안에서의 생존율 (%)")
print(pct.round(2).to_string())

lo = d["Age"] < CUT
r_lo, r_hi = d.loc[lo, "Survived"].mean(), d.loc[~lo, "Survived"].mean()
print(f"\n{CUT}세 미만 {lo.sum():3d}명  생존율 {r_lo:.4f}")
print(f"{CUT}세 이상 {(~lo).sum():3d}명  생존율 {r_hi:.4f}")
print(f"차이 {r_lo - r_hi:+.4f}")

print(f"\n상관계수(파이) {np.corrcoef(d['Survived'], lo.astype(int))[0, 1]:+.4f}")

a, b = tab.iloc[0, 1], tab.iloc[0, 0]      # 어린 쪽 생존/사망
c, e = tab.iloc[1, 1], tab.iloc[1, 0]      # 나이 든 쪽 생존/사망
print(f"오즈비 {(a * e) / (b * c):.4f}")

# --- 파이 계수를 생존율 차이로부터 다시 만들어 본다 ---
# phi = sqrt(pi(1-pi)) / sqrt(ybar(1-ybar)) * (r_lo - r_hi),  pi = 어린 쪽 비율
n = len(d)
pi = lo.sum() / n
ybar = d["Survived"].mean()
print(f"\npi = {lo.sum()}/{n} = {pi:.6f},  ybar = {ybar:.6f}")
print(f"  sqrt(pi(1-pi))       = {np.sqrt(pi * (1 - pi)):.6f}")
print(f"  sqrt(ybar(1-ybar))   = {np.sqrt(ybar * (1 - ybar)):.6f}")
print(f"  차이 Delta           = {r_lo - r_hi:.6f}")
print(f"  식이 주는 phi        = {np.sqrt(pi * (1 - pi)) / np.sqrt(ybar * (1 - ybar)) * (r_lo - r_hi):.6f}")
print(f"  numpy 가 준 phi      = {np.corrcoef(d['Survived'], lo.astype(int))[0, 1]:.6f}")
print(f"  오즈비를 유리수로     = {a * e}/{b * c} = {(a * e) / (b * c):.6f}")

전체 891명 중 나이가 있는 714명으로 분석한다
나이 범위 0.42 ~ 80.00세

도수
Survived     0    1
Age_Group
<21         98   82
>=21       326  208

집단 안에서의 생존율 (%)
Survived       0      1
Age_Group
<21        54.44  45.56
>=21       61.05  38.95

21세 미만 180명  생존율 0.4556
21세 이상 534명  생존율 0.3895
차이 +0.0660

상관계수(파이) +0.0584
오즈비 1.3114

pi = 180/714 = 0.252101,  ybar = 0.406162
  sqrt(pi(1-pi))       = 0.434219
  sqrt(ybar(1-ybar))   = 0.491116
  차이 Delta           = 0.066042
  식이 주는 phi        = 0.058391
  numpy 가 준 phi      = 0.058391
  오즈비를 유리수로     = 26732/20384 = 1.311421

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import seaborn as sns

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
# 글꼴을 바꾸면 마이너스 기호가 깨지므로 unicode_minus 도 함께 꺼 준다.
plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

KO = [f"{CUT}세 미만", f"{CUT}세 이상"]
DIED, LIVED = "#90A4AE", "#1565C0"

fig, ax = plt.subplots(2, 2, figsize=(13, 9))

# (1) 도수 -- 두 집단의 크기가 크게 다르다 (180 대 534)
sns.countplot(data=d, x="Age_Group", hue="Survived", order=ORDER,
              ax=ax[0, 0], palette=[DIED, LIVED])
ax[0, 0].set_title(f"(1) 나이 집단별 도수 (절단점 {CUT}세)")
ax[0, 0].set_xlabel("나이 집단"); ax[0, 0].set_ylabel("사람 수")
ax[0, 0].set_xticks([0, 1]); ax[0, 0].set_xticklabels(KO)
ax[0, 0].legend(["사망", "생존"], title="생존 여부")

# (2) 생존율 -- 막대 두 개의 높이가 비슷하다
rate = d.groupby("Age_Group")["Survived"].mean().reindex(ORDER) * 100
rate.plot(kind="bar", ax=ax[0, 1], color=["#E65100", "#1565C0"],
          edgecolor="black", width=0.7)
ax[0, 1].set_title(f"(2) 나이 집단별 생존율 (%) (절단점 {CUT}세)")
ax[0, 1].set_xlabel("나이 집단"); ax[0, 1].set_ylabel("생존율 (%)")
ax[0, 1].set_xticks([0, 1]); ax[0, 1].set_xticklabels(KO)
ax[0, 1].tick_params(axis="x", rotation=0)
ax[0, 1].grid(True, alpha=0.3, axis="y"); ax[0, 1].set_ylim(0, 60)
for i, v in enumerate(rate):
    ax[0, 1].text(i, v + 1.5, f"{v:.1f}%", ha="center", fontweight="bold")

# (3) 누적 막대
tab.plot(kind="bar", stacked=True, ax=ax[1, 0],
         color=["#D32F2F", "#33691E"], edgecolor="black", width=0.7)
ax[1, 0].set_title("(3) 누적 막대 (도수)")
ax[1, 0].set_xlabel("나이 집단"); ax[1, 0].set_ylabel("사람 수")
ax[1, 0].set_xticks([0, 1]); ax[1, 0].set_xticklabels(KO)
ax[1, 0].tick_params(axis="x", rotation=0)
ax[1, 0].legend(["사망", "생존"], title="생존 여부")

# (4) 비율 열지도 -- 색 범위를 0~1 로 고정한다
prop = pd.crosstab(d["Age_Group"], d["Survived"],
                   normalize="index").reindex(ORDER)
sns.heatmap(prop, annot=True, fmt=".2%", cmap="Blues", ax=ax[1, 1],
            vmin=0, vmax=1, cbar_kws={"label": "비율"},
            linewidths=2, linecolor="black")
ax[1, 1].set_title("(4) 비율 열지도")
ax[1, 1].set_xlabel("생존 여부"); ax[1, 1].set_ylabel("나이 집단")
ax[1, 1].set_xticklabels(["사망", "생존"])
ax[1, 1].set_yticklabels(KO, rotation=0)

fig.suptitle(f"타이타닉: 나이 집단과 생존, 절단점 {CUT}세", y=1.00)
fig.tight_layout()
fig.savefig("titanic_age21_four.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# --- 세로축을 자르면 막대 높이의 비가 어떻게 바뀌는가 ---
for lohi in [(0, 60), (38, 46)]:
    l, h = lohi
    h1 = (rate.iloc[0] - l) / (h - l)
    h2 = (rate.iloc[1] - l) / (h - l)
    print(f"세로축 [{l}, {h}]: 그려지는 높이 {h1:.4f} 대 {h2:.4f},  비 {h1 / h2:.2f}")
print(f"생존율 자체의 비 = {rate.iloc[0] / rate.iloc[1]:.4f}")
print(f"\n두 집단의 크기 {int(tab.iloc[0].sum())}명 대 {int(tab.iloc[1].sum())}명 "
      f"({tab.iloc[1].sum() / tab.iloc[0].sum():.2f} 배)")

세로축 [0, 60]: 그려지는 높이 0.7593 대 0.6492,  비 1.17
세로축 [38, 46]: 그려지는 높이 0.9444 대 0.1189,  비 7.94
생존율 자체의 비 = 1.1696

두 집단의 크기 180명 대 534명 (2.97 배)

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
PALETTE = ["#E65100", "#33691E"]          # 사망 / 생존

# 생존/사망 각각의 나이 분포를 겹쳐 그린다.
sns.histplot(data=d, x="Age", hue="Survived", kde=True, ax=ax[0],
             palette=PALETTE, bins=30)
ax[0].axvline(CUT, color="#D32F2F", ls="--", lw=2)
ax[0].text(CUT + 1, ax[0].get_ylim()[1] * 0.92, f"절단점 {CUT}세",
           color="#D32F2F", fontweight="bold")
ax[0].set_title("생존 여부별 나이 분포")
ax[0].set_xlabel("나이 (세)"); ax[0].set_ylabel("사람 수")
ax[0].legend(["생존", "사망"], title="생존 여부")

sns.boxplot(data=d, x="Survived", y="Age", ax=ax[1], hue="Survived",
            palette=PALETTE, legend=False)
ax[1].axhline(CUT, color="#D32F2F", ls="--", lw=2)
ax[1].text(1.35, CUT + 1.5, f"절단점 {CUT}세", color="#D32F2F",
           fontweight="bold")
ax[1].set_title("생존 여부별 나이 상자그림")
ax[1].set_xlabel("생존 여부"); ax[1].set_ylabel("나이 (세)")
ax[1].set_xticks([0, 1])
ax[1].set_xticklabels(["사망", "생존"])

fig.tight_layout()
fig.savefig("titanic_age_dist.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# 그림이 말하는 것을 수치로도 확인한다.
print(f"{'':14s}{'n':>6s}{'평균':>8s}{'중앙값':>8s}{'Q1':>7s}{'Q3':>7s}")
for lab, v in [("사망", d.loc[d["Survived"] == 0, "Age"]),
               ("생존", d.loc[d["Survived"] == 1, "Age"])]:
    print(f"{lab:14s}{len(v):>6d}{v.mean():>8.2f}{v.median():>8.2f}"
          f"{v.quantile(.25):>7.2f}{v.quantile(.75):>7.2f}")

# --- 상자그림의 수염과 울타리: 어린아이는 이상치로도 찍히지 않는다 ---
print("\n상자그림이 그리는 것 (울타리 = Q1-1.5IQR, Q3+1.5IQR)")
for lab, v in [("사망", d.loc[d["Survived"] == 0, "Age"]),
               ("생존", d.loc[d["Survived"] == 1, "Age"])]:
    q1, q3 = v.quantile(.25), v.quantile(.75)
    iqr = q3 - q1
    print(f"  {lab}: 최솟값 {v.min():5.2f},  아래 울타리 {q1 - 1.5 * iqr:6.2f}"
          f"  ->  아래쪽 이상치 {int((v < q1 - 1.5 * iqr).sum())}명")

# --- 평균 차이가 어디서 오는가 ---
m0 = d["Survived"] == 0
dm = d.loc[~m0, "Age"].mean() - d.loc[m0, "Age"].mean()
sp = np.sqrt(((m0.sum() - 1) * d.loc[m0, "Age"].var()
              + ((~m0).sum() - 1) * d.loc[~m0, "Age"].var()) / (len(d) - 2))
print(f"\n평균 차이 (생존 - 사망) = {dm:+.4f}세,  합동 표준편차 {sp:.4f},  표준화 차이 {dm / sp:+.4f}")
big = d[d["Age"] >= 10]
m0b = big["Survived"] == 0
print(f"  10세 미만 {int((d['Age'] < 10).sum())}명을 빼면 평균 차이 = "
      f"{big.loc[~m0b, 'Age'].mean() - big.loc[m0b, 'Age'].mean():+.4f}세")
print(f"  10세 미만의 생존율 {d.loc[d['Age'] < 10, 'Survived'].mean():.4f}"
      f"  (전체 {d['Survived'].mean():.4f})")
print(f"  생존자 중 10세 미만 비율 {(d.loc[~m0, 'Age'] < 10).mean():.4f},"
      f"  사망자 중 {(d.loc[m0, 'Age'] < 10).mean():.4f}")

# --- 결측 177명을 버리는 대신 메우면 답이 바뀌는가 ---
def cut21(age, sur, tag):
    m = age < CUT
    a, b = sur[m].sum(), m.sum() - sur[m].sum()
    c, e = sur[~m].sum(), (~m).sum() - sur[~m].sum()
    print(f"  {tag:26s} n={len(age):3d}  n<21={m.sum():3d}  r<={sur[m].mean():.4f}"
          f"  r>={sur[~m].mean():.4f}  차 {sur[m].mean() - sur[~m].mean():+.4f}"
          f"  OR {(a * e) / (b * c):.4f}")

print("\n결측 177명을 어떻게 다루는가")
cut21(d["Age"].to_numpy(), d["Survived"].to_numpy(), "(ㄱ) 버린다")
med = d["Age"].median()
cut21(df["Age"].fillna(med).to_numpy(), df["Survived"].to_numpy(),
      f"(ㄴ) 중앙값 {med:.1f} 로 메운다")
by = df.groupby(["Sex", "Pclass"])["Age"].transform(lambda s: s.fillna(s.median()))
cut21(by.to_numpy(), df["Survived"].to_numpy(), "(ㄷ) 성별·등급 중앙값")
print(f"  성별·등급 중앙값은 모두 21세 이상이라 메운 177명이 전원 '21세 이상' 으로 간다")
print(f"  중앙값으로 메우면 28.0세인 사람이 {int((d['Age'] == med).sum())}명에서 "
      f"{int((df['Age'].fillna(med) == med).sum())}명으로 늘어난다")

                   n      평균     중앙값     Q1     Q3
사망               424   30.63   28.00  21.00  39.00
생존               290   28.34   28.00  19.00  36.00

상자그림이 그리는 것 (울타리 = Q1-1.5IQR, Q3+1.5IQR)
  사망: 최솟값  1.00,  아래 울타리  -6.00  ->  아래쪽 이상치 0명
  생존: 최솟값  0.42,  아래 울타리  -6.50  ->  아래쪽 이상치 0명

평균 차이 (생존 - 사망) = -2.2825세,  합동 표준편차 14.4933,  표준화 차이 -0.1575
  10세 미만 62명을 빼면 평균 차이 = -0.0680세
  10세 미만의 생존율 0.6129  (전체 0.4062)
  생존자 중 10세 미만 비율 0.1310,  사망자 중 0.0566

결측 177명을 어떻게 다루는가
  (ㄱ) 버린다                    n=714  n<21=180  r<=0.4556  r>=0.3895  차 +0.0660  OR 1.3114
  (ㄴ) 중앙값 28.0 로 메운다         n=891  n<21=180  r<=0.4556  r>=0.3657  차 +0.0899  OR 1.4514
  (ㄷ) 성별·등급 중앙값              n=891  n<21=180  r<=0.4556  r>=0.3657  차 +0.0899  OR 1.4514
  성별·등급 중앙값은 모두 21세 이상이라 메운 177명이 전원 '21세 이상' 으로 간다
  중앙값으로 메우면 28.0세인 사람이 25명에서 202명으로 늘어난다

age = d["Age"].to_numpy()
sur = d["Survived"].to_numpy()

# 양쪽에 최소 10명은 남는 절단점만 고려한다.
CUTS = [c for c in range(1, 80)
        if (age < c).sum() >= 10 and (age >= c).sum() >= 10]
print(f"검사 대상 절단점 {len(CUTS)}개 ({CUTS[0]}세 ~ {CUTS[-1]}세)\n")

print(f"{'절단점':>7s}{'n<c':>6s}{'n>=c':>6s}{'생존율<':>9s}"
      f"{'생존율>=':>10s}{'차이':>9s}{'오즈비':>9s}")
gaps = []
for c in CUTS:
    m = age < c
    lo_r, hi_r = sur[m].mean(), sur[~m].mean()
    gaps.append((c, lo_r - hi_r))
    if c % 5 == 0:                           # 5세 간격만 출력한다
        odds = (lo_r / (1 - lo_r)) / (hi_r / (1 - hi_r))
        print(f"{c:>7d}{m.sum():>6d}{(~m).sum():>6d}{lo_r:>9.4f}"
              f"{hi_r:>10.4f}{lo_r - hi_r:>+9.4f}{odds:>9.4f}")

wide = max(gaps, key=lambda z: z[1])
narrow = min(gaps, key=lambda z: abs(z[1]))
print(f"\n차이가 가장 큰 절단점:   {wide[0]}세, {wide[1]:+.4f}")
print(f"차이가 가장 작은 절단점: {narrow[0]}세, {narrow[1]:+.4f}")
print(f"차이가 10%포인트를 넘는 절단점 "
      f"{sum(1 for c, g in gaps if abs(g) > 0.10)}개 / {len(CUTS)}개")

# --- Delta 를 ybar 와 pi 만으로 다시 쓴다 ---
# Delta(c) = (r_<c - ybar) / (1 - pi(c)),  pi(c) = (age<c) 의 비율
# phi(c)   = sqrt(pi(1-pi)/(ybar(1-ybar))) * Delta(c)
ybar = sur.mean()
print(f"\n{'c':>4}{'pi':>9}{'r<c':>9}{'Delta':>10}{'(r<c-y)/(1-pi)':>17}{'phi':>9}{'numpy phi':>11}")
for c in (2, 5, 21, 30, 65):
    m = age < c
    pi = m.mean()
    dl = sur[m].mean() - sur[~m].mean()
    ph = np.sqrt(pi * (1 - pi) / (ybar * (1 - ybar))) * dl
    print(f"{c:>4}{pi:>9.4f}{sur[m].mean():>9.4f}{dl:>+10.4f}"
          f"{(sur[m].mean() - ybar) / (1 - pi):>+17.4f}{ph:>+9.4f}"
          f"{np.corrcoef(sur, m.astype(int))[0, 1]:>+11.4f}")

# --- 차이가 가장 큰 절단점과 연관이 가장 센 절단점은 다르다 ---
rows = []
for c in CUTS:
    m = age < c
    pi = m.mean()
    dl = sur[m].mean() - sur[~m].mean()
    rows.append((c, pi, dl, np.sqrt(pi * (1 - pi) / (ybar * (1 - ybar))) * dl))
bd = max(rows, key=lambda z: abs(z[2]))
bp = max(rows, key=lambda z: abs(z[3]))
print(f"\n|Delta| 가 가장 큰 절단점: {bd[0]}세  pi={bd[1]:.4f}  Delta={bd[2]:+.4f}  phi={bd[3]:+.4f}")
print(f"|phi|   가 가장 큰 절단점: {bp[0]}세  pi={bp[1]:.4f}  Delta={bp[2]:+.4f}  phi={bp[3]:+.4f}")
print("|phi| 상위 다섯:", [(c, round(p, 4)) for c, _, _, p in sorted(rows, key=lambda z: -abs(z[3]))[:5]])
print("|Delta| 상위 다섯:", [(c, round(g, 4)) for c, _, g, _ in sorted(rows, key=lambda z: -abs(z[2]))[:5]])

검사 대상 절단점 64개 (2세 ~ 65세)

    절단점   n<c  n>=c     생존율<     생존율>=       차이      오즈비
      5    40   674   0.6750    0.3902  +0.2848   3.2457
     10    62   652   0.6129    0.3865  +0.2264   2.5132
     15    78   636   0.5769    0.3852  +0.1917   2.1763
     20   164   550   0.4817    0.3836  +0.0981   1.4932
     25   278   436   0.4245    0.3945  +0.0300   1.1320
     30   384   330   0.4062    0.4061  +0.0002   1.0008
     35   479   235   0.4092    0.4000  +0.0092   1.0389
     40   551   163   0.4156    0.3742  +0.0414   1.1892
     45   599   115   0.4124    0.3739  +0.0384   1.1749
     50   640    74   0.4109    0.3649  +0.0461   1.2144
     55   672    42   0.4122    0.3095  +0.1027   1.5644
     60   688    26   0.4113    0.2692  +0.1421   1.8966
     65   703    11   0.4111    0.0909  +0.3202   6.9807

차이가 가장 큰 절단점:   2세, +0.4600
차이가 가장 작은 절단점: 30세, +0.0002
차이가 10%포인트를 넘는 절단점 29개 / 64개

   c       pi      r<c     Delta   (r<c-y)/(1-pi)      phi  numpy phi
   2   0.0196   0.8571   +0.4600          +0.4600  +0.1299    +0.1299
   5   0.0560   0.6750   +0.2848          +0.2848  +0.1334    +0.1334
  21   0.2521   0.4556   +0.0660          +0.0660  +0.0584    +0.0584
  30   0.5378   0.4062   +0.0002          +0.0002  +0.0002    +0.0002
  65   0.9846   0.4111   +0.3202          +0.3202  +0.0803    +0.0803

|Delta| 가 가장 큰 절단점: 2세  pi=0.0196  Delta=+0.4600  phi=+0.1299
|phi|   가 가장 큰 절단점: 7세  pi=0.0658  Delta=+0.3168  phi=+0.1600
|phi| 상위 다섯: [(7, 0.16), (6, 0.1557), (8, 0.153), (9, 0.1517), (16, 0.136)]
|Delta| 상위 다섯: [(2, 0.46), (64, 0.3353), (65, 0.3202), (6, 0.318), (7, 0.3168)]

fig, ax = plt.subplots(1, 2, figsize=(12, 4.8))

dif = [sur[age < c].mean() - sur[age >= c].mean() for c in CUTS]

# (1) 절단점에 따른 생존율 차이
ax[0].plot(CUTS, dif, "o-", ms=4, color="#33691E")
ax[0].axhline(0, color="#37474F", lw=1)
ax[0].axvline(21, color="#E65100", ls=":", lw=2, label="절단점 21세")
ax[0].set_xlabel("나이 절단점 $c$ (세)")
# mathtext 는 \ge 를 모른다. \geq 로 써야 한다.
ax[0].set_ylabel("생존율 차이  $r(<c) - r(\\geq c)$")
ax[0].set_title("절단점에 따른 생존율 차이")
ax[0].legend(); ax[0].grid(alpha=0.25)

# (2) 5년 구간별 생존율 -- 자르지 않고 본 실제 모양
edges = np.arange(0, 85, 5)
mid, rt, cnt = [], [], []
for i in range(len(edges) - 1):
    m = (age >= edges[i]) & (age < edges[i + 1])
    if m.sum() >= 5:
        mid.append((edges[i] + edges[i + 1]) / 2)
        rt.append(sur[m].mean()); cnt.append(m.sum())
ax[1].plot(mid, rt, "o-", color="#6A1B9A", ms=6)
for x_, y_, n_ in zip(mid, rt, cnt):        # 각 점 위에 표본 크기를 적는다
    ax[1].annotate(f"{n_}", (x_, y_), textcoords="offset points",
                   xytext=(0, 7), ha="center", fontsize=7, color="#90A4AE")
ax[1].axhline(sur.mean(), color="#37474F", ls="--", lw=1,
              label=f"전체 생존율 {sur.mean():.3f}")
ax[1].axvline(21, color="#E65100", ls=":", lw=2, label="절단점 21세")
ax[1].set_xlabel("나이 (5년 구간, 점 위의 수는 구간별 인원)")
ax[1].set_ylabel("생존율")
ax[1].set_title("관계가 단조가 아니다")
ax[1].legend(); ax[1].grid(alpha=0.25)

fig.suptitle("자유로운 모수 하나, 여러 개의 결론", y=1.02)
fig.tight_layout()
fig.savefig("titanic_cutoff_sweep.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# 10년 단위로 묶어 숫자로도 본다.
g = d.groupby(pd.cut(d["Age"], [0, 10, 20, 30, 40, 50, 60, 81],
                     right=False))["Survived"].agg(["count", "sum", "mean"])
print(g.round(4).to_string())

# --- 21세 미만의 생존율을 두 조각의 가중평균으로 분해한다 ---
m1 = age < 10
m2 = (age >= 10) & (age < 21)
w1 = m1.sum() / (m1 | m2).sum()
print(f"\n10세 미만   n={m1.sum():3d}  생존 {int(sur[m1].sum()):3d}  r={sur[m1].mean():.6f}")
print(f"10~21세     n={m2.sum():3d}  생존 {int(sur[m2].sum()):3d}  r={sur[m2].mean():.6f}")
print(f"  가중평균 {w1:.6f}*{sur[m1].mean():.6f} + {1 - w1:.6f}*{sur[m2].mean():.6f}"
      f" = {w1 * sur[m1].mean() + (1 - w1) * sur[m2].mean():.6f}")
print(f"  실제 21세 미만 생존율                                      = {sur[m1 | m2].mean():.6f}")
print(f"  전체 생존율 ybar = {sur.mean():.6f}")

          count  sum    mean
Age
[0, 10)      62   38  0.6129
[10, 20)    102   41  0.4020
[20, 30)    220   77  0.3500
[30, 40)    167   73  0.4371
[40, 50)     89   34  0.3820
[50, 60)     48   20  0.4167
[60, 81)     26    7  0.2692

10세 미만   n= 62  생존  38  r=0.612903
10~21세     n=118  생존  44  r=0.372881
  가중평균 0.344444*0.612903 + 0.655556*0.372881 = 0.455556
  실제 21세 미만 생존율                                      = 0.455556
  전체 생존율 ybar = 0.406162

21세 미만 = [0~10세의 0.613]  +  [10~20세의 0.402]  -> 섞여서 0.456
21세 이상 = 20대부터 60대까지 대체로 0.35~0.44   -> 0.390

두 평균이 비슷해진다 -> 차이 +0.066

age = d["Age"].to_numpy()
sur = d["Survived"].to_numpy()
CUTS = [c for c in range(1, 80)
        if (age < c).sum() >= 10 and (age >= c).sum() >= 10]

gaps = [(c, sur[age < c].mean() - sur[age >= c].mean()) for c in CUTS]
gaps.sort(key=lambda z: abs(z[1]))
print("생존율 차이가 0에 가장 가까운 절단점 다섯")
print(f"{'절단점':>7s}{'n<c':>6s}{'생존율<':>9s}{'생존율>=':>10s}{'차이':>11s}")
for c, g in gaps[:5]:
    print(f"{c:>7d}{(age < c).sum():>6d}{sur[age < c].mean():>9.4f}"
          f"{sur[age >= c].mean():>10.4f}{g:>+11.6f}")

# 부호가 바뀌는 곳을 찾는다.
seq = [(c, sur[age < c].mean() - sur[age >= c].mean()) for c in CUTS]
print("\n차이의 부호가 바뀌는 구간")
for i in range(1, len(seq)):
    if np.sign(seq[i][1]) != np.sign(seq[i - 1][1]):
        print(f"  {seq[i-1][0]}세 ({seq[i-1][1]:+.6f}) "
              f"-> {seq[i][0]}세 ({seq[i][1]:+.6f})")

생존율 차이가 0에 가장 가까운 절단점 다섯
    절단점   n<c     생존율<     생존율>=         차이
     30   384   0.4062    0.4061  +0.000189
     29   364   0.4066    0.4057  +0.000879
     32   428   0.4065    0.4056  +0.000948
     27   319   0.4075    0.4051  +0.002460
     31   411   0.4039    0.4092  -0.005348

차이의 부호가 바뀌는 구간
  30세 (+0.000189) -> 31세 (-0.005348)
  31세 (-0.005348) -> 32세 (+0.000948)
  47세 (+0.021336) -> 48세 (-0.010930)
  48세 (-0.010930) -> 49세 (+0.021017)

import numpy as np
import pandas as pd

print("성별로 나눈 나이대별 생존율")
bins = [0, 10, 20, 30, 40, 50, 81]
d["AgeBin"] = pd.cut(d["Age"], bins, right=False)
t = d.pivot_table(index="AgeBin", columns="Sex",
                  values="Survived", aggfunc=["count", "mean"],
                  observed=True)
print(t.round(4).to_string())

성별로 나눈 나이대별 생존율
          count         mean
Sex      female male  female    male
AgeBin
[0, 10)      30   32  0.6333  0.5938
[10, 20)     45   57  0.7556  0.1228
[20, 30)     72  148  0.7222  0.1689
[30, 40)     60  107  0.8333  0.2150
[40, 50)     32   57  0.6875  0.2105
[50, 81)     22   52  0.9091  0.1346

남자아이(10세 미만) 생존율 0.594
남자   (10세 이상) 생존율 0.123 ~ 0.215

-> 남성에게는 '어린이' 지위가 결정적이었다

# 두 열이 각각 어느 쪽으로 움직이는지 첫 구간과 끝 구간으로 견준다.
rates = t["mean"]
for sex in ["female", "male"]:
    col = rates[sex]
    print(f"{sex:7s} 10세 미만 {col.iloc[0]:.4f}  ->  "
          f"50세 이상 {col.iloc[-1]:.4f}   "
          f"({col.iloc[-1] - col.iloc[0]:+.4f})")
    print(f"{'':7s} 최소 {col.min():.4f}, 최대 {col.max():.4f}")

female  10세 미만 0.6333  ->  50세 이상 0.9091   (+0.2758)
        최소 0.6333, 최대 0.9091
male    10세 미만 0.5938  ->  50세 이상 0.1346   (-0.4591)
        최소 0.1228, 최대 0.5938

전체에서 본 나이 효과 (21세 절단점에서 차이 +0.066, 아주 약함)
  = 여성의 상승 추세  와  남성의 하강 추세  를 뭉뚱그린 결과

-> 방향이 반대인 두 흐름이 서로 상쇄되어 약해진 것이다

import numpy as np

a = d.loc[d["Survived"] == 0, "Age"].to_numpy()
b = d.loc[d["Survived"] == 1, "Age"].to_numpy()

print(f"{'분위수':>8s}{'사망':>9s}{'생존':>9s}{'차이':>9s}")
for q in [0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95]:
    qa, qb = np.quantile(a, q), np.quantile(b, q)
    print(f"{q:>8.2f}{qa:>9.2f}{qb:>9.2f}{qb - qa:>+9.2f}")

print(f"\n{'평균':>8s}{a.mean():>9.2f}{b.mean():>9.2f}{b.mean() - a.mean():>+9.2f}")
print(f"{'표준편차':>8s}{a.std(ddof=1):>9.2f}{b.std(ddof=1):>9.2f}")

# 차이가 어느 구간에 몰려 있는지 센다.
print(f"\n10세 미만 비율   사망 {(a < 10).mean():.4f}   생존 {(b < 10).mean():.4f}")
print(f"20~50세 비율    사망 {((a >= 20) & (a < 50)).mean():.4f}   "
      f"생존 {((b >= 20) & (b < 50)).mean():.4f}")

     분위수       사망       생존       차이
    0.05     9.00     2.45    -6.55
    0.10    16.00     5.00   -11.00
    0.25    21.00    19.00    -2.00
    0.50    28.00    28.00    +0.00
    0.75    39.00    36.00    -3.00
    0.90    50.00    49.00    -1.00
    0.95    58.00    54.00    -4.00

      평균    30.63    28.34    -2.28
    표준편차    14.17    14.95

10세 미만 비율   사망 0.0566   생존 0.1310
20~50세 비율    사망 0.6887   생존 0.6345

import numpy as np
import pandas as pd

print("두 절단점 (c1, c2) 로 세 집단을 만든다")
print(f"{'c1':>4s}{'c2':>4s}{'n1':>5s}{'n2':>5s}{'n3':>5s}"
      f"{'r1':>8s}{'r2':>8s}{'r3':>8s}{'최대차':>8s}")
for c1 in [5, 10, 15]:
    for c2 in [20, 30, 40, 50, 60]:
        g = np.digitize(age, [c1, c2])
        ns = [(g == k).sum() for k in range(3)]
        if min(ns) < 10:
            continue
        rs = [sur[g == k].mean() for k in range(3)]
        if c1 == 10 or c2 == 60:
            print(f"{c1:>4d}{c2:>4d}{ns[0]:>5d}{ns[1]:>5d}{ns[2]:>5d}"
                  f"{rs[0]:>8.4f}{rs[1]:>8.4f}{rs[2]:>8.4f}"
                  f"{max(rs) - min(rs):>8.4f}")

# 단일 절단점 10세로 만든 두 집단과 견준다.
m = age < 10
print(f"\n단일 절단점 10세: {m.sum()}명 {sur[m].mean():.4f}  "
      f"vs  {(~m).sum()}명 {sur[~m].mean():.4f}  "
      f"(차이 {sur[m].mean() - sur[~m].mean():+.4f})")

두 절단점 (c1, c2) 로 세 집단을 만든다
  c1  c2   n1   n2   n3      r1      r2      r3     최대차
   5  60   40  648   26  0.6750  0.3951  0.2692  0.4058
  10  20   62  102  550  0.6129  0.4020  0.3836  0.2293
  10  30   62  322  330  0.6129  0.3665  0.4061  0.2464
  10  40   62  489  163  0.6129  0.3906  0.3742  0.2387
  10  50   62  578   74  0.6129  0.3893  0.3649  0.2480
  10  60   62  626   26  0.6129  0.3914  0.2692  0.3437
  15  60   78  610   26  0.5769  0.3902  0.2692  0.3077

단일 절단점 10세: 62명 0.6129  vs  652명 0.3865  (차이 +0.2264)

import numpy as np
import pandas as pd

full = df.copy()
miss = full["Age"].isna()
print(f"결측 {miss.sum()}명의 생존율 {full.loc[miss, 'Survived'].mean():.4f}")
print(f"관측 {(~miss).sum()}명의 생존율 {full.loc[~miss, 'Survived'].mean():.4f}\n")

# 관측된 자료에서의 21세 분할
obs = full[~miss]
lo_n = (obs["Age"] < 21).sum(); lo_s = obs.loc[obs["Age"] < 21, "Survived"].sum()
hi_n = (obs["Age"] >= 21).sum(); hi_s = obs.loc[obs["Age"] >= 21, "Survived"].sum()
m_n = miss.sum(); m_s = full.loc[miss, "Survived"].sum()
print(f"관측:  <21  {lo_s}/{lo_n},  >=21  {hi_s}/{hi_n}")
print(f"결측:  생존 {m_s}/{m_n}\n")

print("결측자를 전부 한쪽에 몰아넣으면")
print(f"{'가정':>26s}{'rate<21':>10s}{'rate>=21':>11s}{'차이':>9s}"
      f"{'오즈비':>9s}")
for lab, (a, n1, c, n2) in {
    "관측 자료만 (실제 분석)": (lo_s, lo_n, hi_s, hi_n),
    "전부 21세 미만": (lo_s + m_s, lo_n + m_n, hi_s, hi_n),
    "전부 21세 이상": (lo_s, lo_n, hi_s + m_s, hi_n + m_n),
}.items():
    p1, p2 = a / n1, c / n2
    odds = (p1 / (1 - p1)) / (p2 / (1 - p2))
    print(f"{lab:>26s}{p1:>10.4f}{p2:>11.4f}{p1 - p2:>+9.4f}"
          f"{odds:>9.4f}")

결측 177명의 생존율 0.2938
관측 714명의 생존율 0.4062

관측:  <21  82/180,  >=21  208/534
결측:  생존 52/177

결측자를 전부 한쪽에 몰아넣으면
                        가정   rate<21   rate>=21       차이      오즈비
            관측 자료만 (실제 분석)    0.4556     0.3895  +0.0660   1.3114
                 전부 21세 미만    0.3754     0.3895  -0.0142   0.9418
                 전부 21세 이상    0.4556     0.3657  +0.0899   1.4514

# 결측자의 다른 특성으로 나이를 짐작해 본다.
print("\n결측 여부에 따른 다른 변수")
for col, lab in [("Pclass", "평균 객실등급"), ("SibSp", "평균 형제자매/배우자"),
                 ("Parch", "평균 부모/자녀"), ("Fare", "평균 운임")]:
    print(f"  {lab:18s} 관측 {full.loc[~miss, col].mean():7.3f}"
          f"   결측 {full.loc[miss, col].mean():7.3f}")

결측 여부에 따른 다른 변수
  평균 객실등급            관측   2.237   결측   2.599
  평균 형제자매/배우자        관측   0.513   결측   0.565
  평균 부모/자녀           관측   0.431   결측   0.181
  평균 운임              관측  34.695   결측  22.159

연속형 변수를 잘라야 하는가?
    │
    ├─ 자료 밖에서 정해진 절단점이 있는가?
    │     (진단 기준, 법적 연령, 규정된 등급)
    │     └─ 예 -> 그 값을 쓴다. 출처를 밝힌다.
    │
    ├─ 해석의 편의를 위해서인가?
    │     └─ 회귀계수도 충분히 해석 가능하다. 자르지 않는다.
    │
    ├─ 관계가 비선형이라서인가?
    │     └─ 스플라인 · 다항항 · 구간 더미를 쓴다.
    │
    └─ 자료를 보고 좋은 절단점을 찾았는가?
          └─ 절대 그대로 보고하지 않는다. (1종 오류 0.52)
             순열검정으로 보정하거나, 독립 표본에서 확인한다.

나이와 생존 (나이 기록 714/891명)

  나이를 자르지 않고 5년 구간별 생존율로 제시한다.
    0~10세 0.613, 20대 0.350, 30대 0.437, 60세 이상 0.269
    -> 비단조. 절단점 하나로는 요약할 수 없다.

  참고로 21세 기준 이분화는 45.6% 대 39.0% (차이 +6.6%포인트,
  오즈비 1.31)이나, 이 값은 절단점 21의 선택에 전적으로
  의존한다. 2세에서는 +46%포인트, 30세에서는 +0.02%포인트다.
  절단점을 해석의 근거로 삼지 않았다.

  결측 177명에 대한 민감도 분석 결과 차이의 범위는
  -0.014 ~ +0.090 으로 부호가 바뀔 수 있다.

import numpy as np, pandas as pd

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
d = pd.read_csv(URL).dropna(subset=["Age"])
sur = d["Survived"].to_numpy()

for name in ("Fare", "Age"):
    v = d[name].to_numpy()
    hi = 101 if name == "Fare" else 80
    cuts = [c for c in range(1, hi)
            if (v < c).sum() >= 10 and (v >= c).sum() >= 10]
    gaps = [(c, sur[v < c].mean() - sur[v >= c].mean()) for c in cuts]
    g = [x for _, x in gaps]
    neg = sum(1 for x in g if x < 0)
    print(f"[{name}] 절단점 {len(cuts)}개 ({cuts[0]}~{cuts[-1]})")
    print(f"   차이 범위 [{min(g):+.4f}, {max(g):+.4f}]")
    print(f"   음수인 절단점 {neg}/{len(g)}개\n")

[Fare] 절단점 94개 (7~100)
   차이 범위 [-0.4167, -0.2349]
   음수인 절단점 94/94개

[Age] 절단점 64개 (2~65)
   차이 범위 [-0.0109, +0.4600]
   음수인 절단점 2/64개

import pandas as pd

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
d = pd.read_csv(URL).dropna(subset=["Age"])

for name, binning in [("동일 폭", pd.cut(d.Age, bins=5)),
                      ("분위수", pd.qcut(d.Age, q=5))]:
    g = d.groupby(binning, observed=True).Survived.agg(["mean", "size"])
    print(f"[{name}]")
    for idx, row in g.iterrows():
        print(f"  {str(idx):>18}  생존율 {row['mean']:.4f}  n={int(row['size']):>4}")
    print()

[동일 폭]
      (0.34, 16.336]  생존율 0.5500  n= 100
    (16.336, 32.252]  생존율 0.3699  n= 346
    (32.252, 48.168]  생존율 0.4043  n= 188
    (48.168, 64.084]  생존율 0.4348  n=  69
      (64.084, 80.0]  생존율 0.0909  n=  11

[분위수]
       (0.419, 19.0]  생존율 0.4817  n= 164
        (19.0, 25.0]  생존율 0.3285  n= 137
        (25.0, 31.8]  생존율 0.3937  n= 127
        (31.8, 41.0]  생존율 0.4375  n= 144
        (41.0, 80.0]  생존율 0.3732  n= 142

import numpy as np, pandas as pd

URL = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
d = pd.read_csv(URL).dropna(subset=["Age"])
age, sur = d["Age"].to_numpy(), d["Survived"].to_numpy()
CUTS = [c for c in range(1, 80)
        if (age < c).sum() >= 10 and (age >= c).sum() >= 10]

rng = np.random.default_rng(0)
best = []
for _ in range(500):
    i = rng.choice(len(age), int(0.8 * len(age)), replace=False)
    a, s = age[i], sur[i]
    gg = [(c, s[a < c].mean() - s[a >= c].mean()) for c in CUTS
          if (a < c).sum() >= 10 and (a >= c).sum() >= 10]
    best.append(max(gg, key=lambda z: z[1])[0])

best = np.array(best)
u, cnt = np.unique(best, return_counts=True)
for j in np.argsort(-cnt)[:6]:
    print(f"  {u[j]:>3}세 : {cnt[j]:>4}회 ({cnt[j]/5:.1f}%)")
print(f"\n  범위 {best.min()}~{best.max()}세,  서로 다른 값 {len(u)}개")

print("\n최적 절단점 주변의 집단 크기")
for c in (2, 3, 64):
    m = age < c
    print(f"  c={c:>3}: n(<c)={m.sum():>4} 생존율 {sur[m].mean():.4f} | "
          f"n(>=c)={(~m).sum():>4} 생존율 {sur[~m].mean():.4f}")