카플란-마이어 생존곡선과 로그순위 검정¶

import numpy as np

def kaplan_meier(times, censored):
    """카플란-마이어 생존함수 추정값을 구한다.

    사건이 일어난 시점마다 "그 직전까지 살아 있던 사람 중 그 시점을
    넘긴 비율"을 곱해 나간다. 중도절단된 사람은 절단 시점까지만
    위험집합에 남아 있다가 조용히 빠진다 — 이것이 중도절단 자료를
    버리지 않고 쓰는 방법이다.

    매개변수
    --------
    times    : 관측된 시각(사건 또는 절단)
    censored : 1 이면 중도절단, 0 이면 사건이 관측됨

    돌려주는 값
    ----------
    t_plot, s_plot : 계단그림에 쓸 시각과 생존확률
    """
    order = np.argsort(times)
    times = times[order]
    censored = censored[order]

    event_times = times[censored == 0]
    unique_events = np.unique(event_times)

    s = 1.0
    t_list = [0.0]
    s_list = [1.0]

    for t_j in unique_events:
        # 위험집합: 그 시점에 아직 사건도 절단도 겪지 않은 사람 수
        n_at_risk = np.sum(times >= t_j)
        d_j = np.sum((times == t_j) & (censored == 0))
        s *= (n_at_risk - d_j) / n_at_risk
        t_list.append(t_j)
        s_list.append(s)

    t_list.append(times.max())
    s_list.append(s_list[-1])

    return np.array(t_list), np.array(s_list)

# --- (1) 의 표와 그린우드 표준오차를 확인한다 ---
times = np.array([3., 5., 7., 7., 10., 12., 15., 18.])
censored = np.array([0, 1, 0, 0, 1, 0, 1, 0])

t, s = kaplan_meier(times, censored)
print("t_plot =", t)
print("s_plot =", np.round(s, 6))

acc, S = 0.0, 1.0
for t_j in np.unique(times[censored == 0]):
    n = int(np.sum(times >= t_j))
    d = int(np.sum((times == t_j) & (censored == 0)))
    S *= (n - d) / n
    if n > d:
        acc += d / (n * (n - d))
    print(f"  t={t_j:5.1f}  n={n}  d={d}  S={S:.6f}  "
          f"그린우드합={acc:.6f}  SE={S * np.sqrt(acc):.6f}")

# --- (2) 절단을 무시한 평균과 곡선 아래 넓이 ---
print(f"절단을 사건으로 보고 낸 평균 = {times.mean():.4f}")
area = np.sum(s[:-2] * np.diff(t[:-1]))
print(f"곡선 아래 넓이 (18 까지의 제한평균생존시간) = {area:.4f}")

t_plot = [ 0.  3.  7. 12. 18. 18.]
s_plot = [1.       0.875    0.583333 0.388889 0.       0.      ]
  t=  3.0  n=8  d=1  S=0.875000  그린우드합=0.017857  SE=0.116927
  t=  7.0  n=6  d=2  S=0.583333  그린우드합=0.101190  SE=0.185561
  t= 12.0  n=3  d=1  S=0.388889  그린우드합=0.267857  SE=0.201269
  t= 18.0  n=1  d=1  S=0.000000  그린우드합=0.267857  SE=0.000000
절단을 사건으로 보고 낸 평균 = 9.6250
곡선 아래 넓이 (18 까지의 제한평균생존시간) = 11.7500

from scipy import stats

def logrank_test(times_1, censored_1, times_2, censored_2):
    """이표본 로그순위 검정.

    사건 시점마다 2x2 분할표를 만들어 관측 사건 수와 기대 사건 수를
    비교한다. 그 차이를 모든 시점에 걸쳐 누적한 것이 통계량이다.
    두 생존곡선이 같다는 귀무가설 아래에서 자유도 1 인 카이제곱을 따른다.

    돌려주는 값
    ----------
    chi2, p_value
    """
    event_1 = times_1[censored_1 == 0]
    event_2 = times_2[censored_2 == 0]
    all_event_times = np.unique(np.concatenate([event_1, event_2]))

    O1, E1, V = 0.0, 0.0, 0.0

    for t_j in all_event_times:
        r1 = np.sum(times_1 >= t_j)
        r2 = np.sum(times_2 >= t_j)
        r  = r1 + r2

        d1 = np.sum(event_1 == t_j)
        d2 = np.sum(event_2 == t_j)
        d  = d1 + d2

        # 두 집단의 생존이 같다면, 그 시점의 사건은 위험집합 크기에
        # 비례해 나뉘어야 한다. 그것이 기대 사건 수 e1 이다.
        e1 = r1 * d / r if r > 0 else 0
        v  = r1 * r2 * d * (r - d) / (r**2 * (r - 1)) if r > 1 else 0

        O1 += d1
        E1 += e1
        V  += v

    chi2 = (O1 - E1)**2 / V if V > 0 else 0
    p_value = stats.chi2(1).sf(chi2)
    return chi2, p_value

# --- (1) 의 합산과 순위 불변성을 확인한다 ---
t1 = np.array([5., 8., 12., 15.]); c1 = np.array([0, 0, 1, 0])
t2 = np.array([2., 4., 6., 10.]);  c2 = np.array([0, 0, 1, 0])

chi2, p = logrank_test(t1, c1, t2, c2)
print(f"chi2 = {chi2:.6f},  p = {p:.4f}")

for name, g in [("log(1+t)", np.log1p), ("t^3", lambda x: x ** 3),
                ("sqrt(t)", np.sqrt)]:
    c, pp = logrank_test(g(t1), c1, g(t2), c2)
    print(f"  {name:9s}: chi2 = {c:.6f}  p = {pp:.4f}")

chi2 = 1.183373,  p = 0.2767
  log(1+t) : chi2 = 1.183373  p = 0.2767
  t^3      : chi2 = 1.183373  p = 0.2767
  sqrt(t)  : chi2 = 1.183373  p = 0.2767

import matplotlib.pyplot as plt

np.random.seed(0)

# 1집단: 사건이 늦게 일어난다(평균 20). 20%는 중도절단된다.
# 중도절단은 "그 시점까지는 살아 있었다"는 정보를 주므로 버리지 않는다.
n1 = 40
times_1 = np.random.exponential(scale=20, size=n1)
censored_1 = (np.random.rand(n1) < 0.2).astype(int)

# 2집단: 사건이 빨리 일어난다(평균 12).
n2 = 40
times_2 = np.random.exponential(scale=12, size=n2)
censored_2 = (np.random.rand(n2) < 0.2).astype(int)

t1, s1 = kaplan_meier(times_1, censored_1)
t2, s2 = kaplan_meier(times_2, censored_2)

# 그림으로 확인
fig, ax = plt.subplots(figsize=(10, 5))
# where="post" 가 계단을 오른쪽으로 뻗게 한다. 생존함수는 사건이 일어난
# 그 순간에 떨어지고 다음 사건까지 평평하므로, 이 설정이라야 맞다.
ax.step(t1, s1, where="post", linewidth=2, label="Group 1 (slow)")
ax.step(t2, s2, where="post", linewidth=2, label="Group 2 (fast)")
ax.set_xlabel("Time")
ax.set_ylabel("Survival Probability")
ax.set_title("Kaplan-Meier Survival Curves")
ax.set_ylim(-0.02, 1.05)
ax.legend()
plt.tight_layout()
plt.show()

# --- (1) 의 예측과 맞추어 본다 ---
def km_median(t, s):
    idx = np.where(s <= 0.5)[0]
    return t[idx[0]] if len(idx) else np.nan

q1, q2 = censored_1.mean(), censored_2.mean()
print(f"실제 절단 비율: 1집단 {q1:.3f}, 2집단 {q2:.3f}  (설계값 0.200)")
print(f"KM 중앙값          : 1집단 {km_median(t1, s1):6.3f}, 2집단 {km_median(t2, s2):6.3f}")
print(f"참 중앙값 th*ln2    : 1집단 {20 * np.log(2):6.3f}, 2집단 {12 * np.log(2):6.3f}")
print(f"이름표 절단이 뜻하는 값: 1집단 {20 * np.log(2) / (1 - q1):6.3f}, "
      f"2집단 {12 * np.log(2) / (1 - q2):6.3f}")

chi2, p = logrank_test(times_1, censored_1, times_2, censored_2)
print(f"로그순위: chi2 = {chi2:.4f},  p = {p:.5f}")

실제 절단 비율: 1집단 0.325, 2집단 0.175  (설계값 0.200)
KM 중앙값          : 1집단 20.763, 2집단 10.446
참 중앙값 th*ln2    : 1집단 13.863, 2집단  8.318
이름표 절단이 뜻하는 값: 1집단 20.538, 2집단 10.082
로그순위: chi2 = 12.6889,  p = 0.00037