1 들어가며 — overview 의 § 4.1~4.2 를 정전 깊이로
Klein 시리즈 Ch.4 의 첫 번째 deep-dive. Overview 편에서 KM·NA 의 골격을 다뤘다면, 본 편은 그 골격의 모든 뼈와 인대 를 풀어낸다.
| 편 | 주제 |
|---|---|
| Ch.4 Overview | 7 개 절 조망 (KM·NA·CI·CB·평균·중앙·LT·CIF) |
| § 4.1~4.2 (본 편) | KM·NA estimator 의 정의·5 가지 유도·분산·예제·실무 주의점 |
| § 4.3 (예정) | Pointwise CI — linear · log · arcsine 변환 |
| § 4.4 (예정) | Confidence Bands — EP · Hall-Wellner |
| § 4.5 (예정) | Mean · Median (Brookmeyer-Crowley) |
| § 4.6 (예정) | Left-Truncation |
| § 4.7 (예정) | Competing Risks (1-KM · CIF · CP) |
“우측 censoring 데이터의 비모수 생존 추정은 두 단계로 이뤄진다 — (1) 사건 시점 직전 위험집합 \(Y_i\) 중 사건 수 \(d_i\) 의 비율 \(d_i/Y_i\) 를 사건 시점에서의 hazard 로 추정하고, (2) 이 비율을 곱하면 생존함수 KM, 합하면 누적위험 NA 가 된다. 두 추정량은 5 가지 다른 유도를 통해 동일하게 도출되며 — 이 다중 유도가 KM 의 robustness 의 근거 — 점근적으로 동치이지만 소표본에서 NA 가 약간 안정적이다.”
이 단순한 비율 \(d_i/Y_i\) 가 Ch.5~13 의 모든 비모수·반모수·모수 추정의 building block 이다.
2 § 4.1 — Introduction
2.1 왜 비모수 추정이 필요한가
censoring 이 없는 데이터에서는 empirical survival 함수
\[ \widehat{S}_{\text{emp}}(t) = \frac{\#\{i : X_i > t\}}{n} \]
가 자명한 비모수 추정량이다 — 분자는 \(t\) 까지 살아남은 사람 수, 분모는 전체.
문제: censoring 시 분자가 모호. 환자 \(i\) 의 censoring time \(T_i = 8\) 인데 \(t = 10\) 의 생존 여부를 묻는다면? 8 시점까지는 살았지만, 10 시점에서는 죽었는지 살았는지 모른다. “죽었다” 로 처리하면 생존을 과소 추정, “살았다” 로 처리하면 과대 추정.
→ censoring 정보를 likelihood 에 정확히 반영하는 추정량이 필요. KM 이 그 답.
만약 모수 모형 (Exponential, Weibull 등) 을 가정하면 Ch.3 의 likelihood (식 3.5.1) 로 MLE 가 가능하지만 — 분포 가정이 틀리면 결론이 무너진다. 비모수 추정은 분포 가정 없이 데이터의 자연스러운 step function 형태로 \(S\) 를 추정한다.
2.2 표기
Ch.4 전체에서 사용되는 표기를 정리한다.
| 기호 | 정의 | 비고 |
|---|---|---|
| \(X_i\) | 개체 \(i\) 의 (잠재적) 사건 시점 | 관측 안 될 수 있음 |
| \(C_{r,i}\) | 개체 \(i\) 의 (잠재적) censoring 시점 | 관측 안 될 수 있음 |
| \(T_i = \min(X_i, C_{r,i})\) | 관측된 시간 | study time |
| \(\delta_i = I(X_i \leq C_{r,i})\) | 사건 indicator | \(\delta_i = 1\) 사건, \(\delta_i = 0\) censored |
| \(D\) | 서로 다른 사건 시점의 수 | \(D \leq n\) |
| \(t_1 < t_2 < \cdots < t_D\) | 서로 다른 사건 시점 (정렬) | 오름차순 |
| \(d_i\) | \(t_i\) 에서의 사건 수 | tie 가능 (\(d_i \geq 1\)) |
| \(Y_i\) | \(t_i\) 직전 위험집합 크기 | “just prior to \(t_i\)” |
\(Y_i\) 는 시점 \(t_i\) 직전에 여전히 추적 중이고, 아직 사건을 경험하지 않은 개체의 수:
\[ Y_i = \#\{j : T_j \geq t_i\} \]
직관: “\(t_i\) 라는 사건이 일어날 수 있는 후보자 수”. \(t_i\) 시점의 hazard 를 추정할 때 분모가 된다.
\(Y_i\) 는 \(t_i\) 직전 의 사람 수, 즉 \(T_j \geq t_i\) 인 사람 수 (강한 부등호 \(>\) 가 아니라 \(\geq\)).
이유: \(t_i\) 에 사망한 사람도 “\(t_i\) 직전 살아있던 사람” 에 포함되어야 한다. \(d_i\) 는 \(Y_i\) 의 부분집합 (\(d_i \leq Y_i\)).
\(d_i / Y_i\) 의 의미는 따라서:
“\(t_i\) 직전까지 살아있던 사람 중에서, \(t_i\) 에 사건을 겪은 비율”
= \(P(T = t_i \mid T \geq t_i) \approx h(t_i) \cdot dt\)
2.3 핵심 가정 — 독립 (비정보적) censoring
본 chapter 의 모든 추정량은 다음 가정에 의존한다:
\[ X_i \perp C_{r,i} \]
또는 약화된 형태: \(C_{r,i}\) 의 분포가 \(X_i\) 의 분포 모수와 무관 (non-informative censoring, Ch.3.2).
가정 위반 시 발생하는 일:
- 예후 나쁜 환자만 선택적 dropout (예: 부작용으로 약 끊은 환자) → KM 이 사망률을 과소 추정 (살아남는 비율이 인위적으로 높음).
- 좋은 환자 (회복) 만 추적 종료 → KM 이 사망률을 과대 추정.
Klein & Moeschberger (1984) 는 이런 informative censoring 시 추정량이 어떤 다른 함수를 추정하게 되는지를 보여준다.
2.4 어떤 censoring 형태에 적용 가능한가
본 § 4.2 의 추정량은 Ch.3.2 의 다음 형태에 모두 적용된다:
| Censoring 형태 | KM·NA 적용 |
|---|---|
| Type I (사전 종료 시점) | O |
| Type II (첫 r 사건까지) | O |
| Progressive (다단계 sacrifice) | O |
| Generalized Type I (개체별 입학 시점) | O |
| Random (개체별 censoring time) | O |
| Left censoring | X (Ch.5) |
| Interval censoring | X (Ch.5 Turnbull) |
| Left truncation | O (식 4.6.1 의 \(Y_i\) 재정의) |
| Right truncation | X (Ch.5 reverse-time KM) |
핵심 — 우측 censoring 의 모든 형태가 동일한 KM·NA 식으로 처리됨. censoring 형태별로 다른 추정량이 필요하지 않다.
3 § 4.2 — Kaplan-Meier (Product-Limit) Estimator
3.1 정의
\[ \widehat{S}(t) = \begin{cases} 1 & \text{if } t < t_1 \\ \displaystyle\prod_{t_i \leq t} \left(1 - \frac{d_i}{Y_i}\right) & \text{if } t_1 \leq t \end{cases} \]
Kaplan & Meier (1958) 가 도입했고, “product-limit estimator” 라는 이름은 곱 형태 + 이산 분포의 극한이라는 도출 방식에서 유래.
3.2 왜 곱인가 — 사슬 분해 직관
\(t_k\) 까지의 생존을 조건부 확률의 사슬로 분해:
\[ \begin{aligned} S(t_k) &= P(X > t_k) \\ &= \underbrace{P(X > t_k \mid X > t_{k-1})}_{\text{step } k} \cdot \underbrace{P(X > t_{k-1} \mid X > t_{k-2})}_{\text{step } k-1} \cdots \underbrace{P(X > t_1)}_{\text{step } 1} \end{aligned} \]
(이 분해는 이산 분포에서만 정확. 연속 분포에서는 \(S(t_k) = \exp(-H(t_k))\) 의 적분 형태가 자연. 그러나 KM 은 데이터를 이산 분포로 본다 — 점프가 사건 시점에만 있음).
각 조건부 확률을 추정:
\[ \widehat{P}(X > t_i \mid X > t_{i-1}) = \frac{Y_i - d_i}{Y_i} = 1 - \frac{d_i}{Y_i} \]
- 분모 \(Y_i\) : \(t_{i-1}\) 까지 살아있던 후보 수 (이 중 일부는 \(t_i\) 직전에 censoring 으로 사라짐).
- 분자 \(Y_i - d_i\) : \(t_i\) 에 사건 안 겪고 넘어간 사람 수.
이를 \(i = 1, \ldots, k\) 에 대해 곱하면 KM. 따라서 KM 의 곱은 “조건부 생존확률의 사슬” 의 자연스러운 표현이다.
KM 은 step function 으로, 사건 시점 \(t_i\) 에서만 점프한다. 점프 크기:
\[ \widehat{S}(t_i^-) - \widehat{S}(t_i) = \widehat{S}(t_i^-) \cdot \frac{d_i}{Y_i} \]
즉 점프의 절대 크기는 이전까지의 KM 값 \(\widehat{S}(t_i^-)\) 와 현재 비율 \(d_i/Y_i\) 모두에 의존한다.
\(Y_i\) 자체는 \(t_i\) 직전까지의 censoring 패턴에 의해 결정된다. 따라서 동일한 \(d_i\) 라도 censoring 이 많이 발생한 시기 (Y_i 가 작음) 의 사건은 KM 을 더 크게 낮춘다.
→ 한 번의 사건이 곡선에 미치는 영향은 “그 시점에서 위험집합이 얼마나 작았는가” 에 따라 결정. 후반부의 적은 표본에서 발생한 사건이 KM 을 급격히 떨어뜨리는 이유.
3.3 Greenwood 분산 공식
\[ \widehat{V}[\widehat{S}(t)] = \widehat{S}(t)^2 \sum_{t_i \leq t} \frac{d_i}{Y_i (Y_i - d_i)} \]
표준오차: \(\sqrt{\widehat{V}[\widehat{S}(t)]}\).
3.4 Greenwood 도출 — 단계별 (delta method)
Step 1 — log 변환:
\[ \log \widehat{S}(t) = \sum_{t_i \leq t} \log\left(1 - \frac{d_i}{Y_i}\right) \]
각 항이 독립이라고 가정 (실제로는 근사적으로 — counting process 이론에서 정확히 보임).
Step 2 — 각 항의 분산:
\(d_i \sim \text{Binomial}(Y_i, h_i)\) 라고 보면 (위험집합 \(Y_i\) 명 중 \(h_i\) 의 확률로 사건 발생):
\[ \text{Var}(d_i) = Y_i h_i (1 - h_i) \approx Y_i \cdot \frac{d_i}{Y_i} \cdot \frac{Y_i - d_i}{Y_i} = \frac{d_i (Y_i - d_i)}{Y_i} \]
Step 3 — log-항의 분산 (delta method):
\(f(d) = \log(1 - d/Y_i)\), \(f'(d) = -1/(Y_i - d)\). delta method:
\[ \text{Var}\left[\log\left(1 - \frac{d_i}{Y_i}\right)\right] \approx \left(\frac{1}{Y_i - d_i}\right)^2 \cdot \frac{d_i (Y_i - d_i)}{Y_i} = \frac{d_i}{Y_i (Y_i - d_i)} \]
Step 4 — 합과 변환:
\[ \text{Var}[\log \widehat{S}(t)] \approx \sum_{t_i \leq t} \frac{d_i}{Y_i (Y_i - d_i)} \]
Step 5 — \(\widehat{S}\) 분산 (delta method 다시): \(g(\log S) = e^{\log S} = S\), \(g' = S\).
\[ \text{Var}[\widehat{S}(t)] \approx \widehat{S}(t)^2 \cdot \text{Var}[\log \widehat{S}(t)] = \widehat{S}(t)^2 \sum_{t_i \leq t} \frac{d_i}{Y_i (Y_i - d_i)} \]
이것이 Greenwood. → delta method 두 번 + Binomial 분산 = Greenwood. 각 단계가 자연스럽다.
3.5 Nelson-Aalen Estimator
\[ \widetilde{H}(t) = \begin{cases} 0 & \text{if } t \leq t_1 \\ \displaystyle\sum_{t_i \leq t} \frac{d_i}{Y_i} & \text{if } t_1 \leq t \end{cases} \]
Nelson (1972) 가 신뢰성 공학에서 처음 제안하고, Aalen (1978b) 가 counting process 이론으로 재발견. NA 로부터 survival 추정도 가능: \(\widetilde{S}(t) = \exp[-\widetilde{H}(t)]\).
분산 (식 4.2.4): \[ \sigma_H^2(t) = \sum_{t_i \leq t} \frac{d_i}{Y_i^2} \]
3.6 NA 분산 도출
NA 는 합. Counting process 이론에서 \(d_i \sim \text{Poisson}\) 으로 근사하면 \(\text{Var}(d_i) \approx d_i\) (Poisson 의 평균 = 분산). 따라서
\[ \text{Var}\left[\frac{d_i}{Y_i}\right] = \frac{\text{Var}(d_i)}{Y_i^2} \approx \frac{d_i}{Y_i^2} \]
각 항이 독립이라고 가정하면 (counting process martingale 에서 정확히 보임):
\[ \text{Var}[\widetilde{H}(t)] \approx \sum_{t_i \leq t} \frac{d_i}{Y_i^2} \]
→ NA 분산은 Greenwood 보다 훨씬 단순. 이유: NA 는 합 (delta method 한 번 불필요), Poisson 근사가 자연. 이것이 소표본에서 NA 가 KM 보다 안정적인 이유 중 하나.
3.7 KM 의 5 가지 유도 — Why 다섯 길이 모두 KM 으로 모이는가
KM 은 다섯 가지 다른 출발점에서 동일한 추정량으로 도출된다. 이 다중 유도가 KM 의 robustness 의 수학적 근거.
3.7.1 유도 1 — Reduced-Sample (조건부 사슬 분해)
이산 분포 \(S(t)\) 의 점프가 사건 시점에만 있다고 가정 (정보가 없는 censoring 시점에서 점프할 이유 없음).
\[ S(t_i) = \prod_{j \leq i} \frac{S(t_j)}{S(t_{j-1})} = \prod_{j \leq i} P(X > t_j \mid X \geq t_j) \]
(이산 분포에서 \(S(t_{j-1}) = P(X > t_{j-1}) = P(X \geq t_j)\)).
각 조건부 확률을 \(\widehat{P}(X > t_j \mid X \geq t_j) = (Y_j - d_j)/Y_j\) 로 추정. 곱하면 KM. → § 들어가며 의 사슬 분해와 같은 도출이지만 더 명시적.
3.7.2 유도 2 — Redistribute-to-the-Right (Efron 1967)
핵심 아이디어: censoring 이 없으면 각 관측에 mass \(1/n\) 을 주어 empirical CDF 를 구한다. censoring 이 있으면 그 mass 를 “그 이후 시점들에 균등 분배” 한다.
알고리즘:
- 시작: 각 관측에 mass \(1/n\).
- 가장 작은 censoring 시점부터 시작.
- 그 mass 를 더 큰 모든 관측 (사건 + censoring) 에 균등 재분배.
- 다음 censoring 시점에서 반복.
예제 — 10 명 데이터: 3, 4, 5+, 6, 6+, 8+, 11, 14, 15, 16+ (+ = censoring).
| 데이터 | Step 0 | Step 1 (5+) | Step 2 (8+) | Step 3 (16+) | \(\widehat{S}\) |
|---|---|---|---|---|---|
| 3 | 0.100 | 0.100 | 0.100 | 0.100 | 0.900 |
| 4 | 0.100 | 0.100 | 0.100 | 0.100 | 0.800 |
| 5+ | 0.100 | 0.000 | — | — | 0.800 |
| 6 | 0.100 | 0.114 | 0.114 | 0.114 | 0.686 |
| 6+ | 0.100 | 0.114 | 0.000 | — | 0.686 |
| 8+ | 0.100 | 0.114 | 0.000 | — | 0.686 |
| 11 | 0.100 | 0.114 | 0.137 | 0.171 | 0.515 |
| 14 | 0.100 | 0.114 | 0.137 | 0.171 | 0.343 |
| 15 | 0.100 | 0.114 | 0.137 | 0.171 | 0.171 |
| 16+ | 0.100 | 0.114 | 0.137 | 0.171 | 0.000* |
(*Efron 의 tail 처리 — 마지막 mass 를 그 점에 두면 곡선이 0 으로 떨어짐).
Step 1 (5+ 처리): 5+ 의 mass 0.100 을 더 큰 7 개 관측 (6, 6+, 8+, 11, 14, 15, 16+) 에 균등 분배: 각각 \(0.100/7 = 0.0143\) 추가. → 새로운 mass \(0.100 + 0.0143 = 0.1143 \approx 0.114\).
Step 2 (6+ 와 8+ 처리): 6+ 의 mass 와 8+ 의 mass 를 더 큰 5 개 관측 (11, 14, 15, 16+ 그리고 8+ 의 경우 자기 이후) 에 분배. 결과 0.137.
Step 3 (16+ 처리): Efron 안에서는 mass 를 그 자리에 두므로 \(\widehat{S}(16+) = 0\). Gill 안에서는 \(\widehat{S}(t > 16) = \widehat{S}(15) = 0.171\) 유지.
→ Efron (1967) 은 이 알고리즘의 결과가 KM 식 4.2.1 과 동일함을 증명. 직관적으로 매우 깔끔한 유도.
3.7.3 유도 3 — Self-Consistency (Efron 1967)
추정량 \(\widehat{S}\) 가 다음 fixed-point 방정식을 만족하면 self-consistent:
\[ \widehat{S}(t) = \frac{1}{n}\left[\sum_{T_i > t} \phi(T_i) + \sum_{\delta_i = 0, T_i \leq t} \frac{\widehat{S}(t)}{\widehat{S}(T_i)}\right] \]
여기서 \(\phi(y) = I(y > t)\).
Censoring 이 없다면 \(\widehat{S}(t) = \frac{1}{n} \sum_i I(X_i > t)\) — empirical survival 함수.
Censoring 이 있을 때:
- \(T_i > t\) 인 관측 (사건 또는 censored): 100% 가 \(t\) 보다 큼이 확실 → 1 로 카운트.
- \(T_i \leq t\) 이고 \(\delta_i = 1\) (사건): 0% 가 \(t\) 보다 큼 → 0 으로 카운트.
- \(T_i \leq t\) 이고 \(\delta_i = 0\) (censored, 즉 \(T_i\) 시점에 살아있었지만 그 이후 추적 끊김): 이 사람의 진짜 사건 시점 \(X_i > t\) 일 확률을 추정해야 함. 그 확률이 바로 \(S(t)/S(T_i)\).
- 이유: \(X_i > T_i\) 임은 알고 (censoring 시점 이후도 살아있었음을 의미하므로 \(X_i > T_i\)), \(X_i > t\) 일 조건부 확률 = \(P(X > t \mid X > T_i) = S(t)/S(T_i)\).
- 모르는 \(S\) 를 \(\widehat{S}\) 로 대체 → fixed-point 방정식.
→ “추정량이 자기 자신을 일관되게 사용” 한다는 의미에서 self-consistent.
EM 알고리즘과의 관계: 식 4.2.7 은 EM 의 원형. censored 관측의 missing event time 을 conditional expectation 으로 대체하는 E-step + 새 추정량으로 업데이트하는 M-step 의 반복. Efron (1967) 은 이 반복이 KM 으로 수렴함을 증명. → KM 은 EM 의 자명한 비모수 적용.
3.7.4 유도 4 — Counting Process (Aalen 1975, Andersen et al. 1993)
Ch.3.6 에서 다룬 counting process 이론. \(N(t)\) = 사건 카운트, \(Y(t)\) = at-risk indicator. KM 의 product integral 표현:
\[ \widehat{S}(t) = \prod_{u \leq t} \left(1 - \frac{dN(u)}{Y(u)}\right) \]
이산 점프에서는 \(dN(u) = d_i\), \(Y(u) = Y_i\), 곱은 식 4.2.1 과 일치. Andersen-Borgan-Gill-Keiding (1993) 은 martingale 이론으로 KM 의 점근 정규성·일관성·신뢰대를 통일적으로 도출.
연속 분포에서 \(S(t) = \exp(-\int_0^t h(u) du)\). Hazard 의 무한소 기여 \(h(u) du\) 가 누적된 형태.
이산 점프에서는 \(h\) 가 atom 으로 변환: \(h(u) du \to dN(u)/Y(u)\). 무한소 기여의 곱 \(\prod_{u \leq t}(1 - h(u)du) \to \exp(-\int h du)\) (Riemann product → exponential).
→ KM 은 hazard 의 product integral. counting process language 에서는 이게 가장 자연스러운 표현.
3.7.5 유도 5 — Nonparametric MLE (NPMLE)
분포 가족: 점프가 사건 시점 \(t_1 < \cdots < t_D\) 에만 있는 모든 이산 분포.
각 시점에서의 hazard \(h_i = P(X = t_i \mid X \geq t_i)\) 를 자유 모수로. likelihood (Ch.3 식 3.5.1):
\[ L(h_1, \ldots, h_D) \propto \prod_{i=1}^D h_i^{d_i} (1 - h_i)^{Y_i - d_i} \]
이는 각 \(h_i\) 에 대한 독립 Binomial likelihood. 각 항을 최대화:
\[ \widehat{h}_i = \frac{d_i}{Y_i} \]
\(S(t) = \prod_{t_i \leq t}(1 - h_i)\) 의 MLE (invariance):
\[ \widehat{S}(t) = \prod_{t_i \leq t}\left(1 - \frac{d_i}{Y_i}\right) \]
→ KM. Nelson-Aalen 도 비슷하게 NPMLE (Klein Theoretical Note 5).
왜 5 가지 유도가 모두 같은 결과인가: 다섯 가지 모두 “censoring 의 정보를 likelihood 에 정확히 반영” 하는 다른 표현일 뿐. 본질적으로 KM 의 robustness 가 매우 깊다는 의미.
3.8 KM ↔︎ NA 의 관계 — Taylor 1차
\(\log(1 - x) = -x - x^2/2 - x^3/3 - \cdots\). \(x = d_i/Y_i\) 가 작으면 (즉 위험집합이 충분히 큼):
\[ -\log(1 - d_i/Y_i) \approx d_i/Y_i \]
따라서
\[ -\log \widehat{S}(t) = -\sum_{t_i \leq t} \log\left(1 - \frac{d_i}{Y_i}\right) \approx \sum_{t_i \leq t} \frac{d_i}{Y_i} = \widetilde{H}(t) \]
→ NA 는 \(-\log\) KM 의 Taylor 1차 근사 (Klein Theoretical Note 7).
따라서 \(\widetilde{S}(t) = \exp[-\widetilde{H}(t)] \approx \exp[\log \widehat{S}(t)] = \widehat{S}(t)\).
차이의 부호: \(-\log(1-x) > x\) for \(x \in (0, 1)\) → \(\widetilde{H} > -\log \widehat{S}\) → \(\widetilde{S} > \widehat{S}\). 즉 NA-기반 survival 이 KM 보다 약간 큼.
큰 위험집합 \(Y_i\) 에서는 차이가 무시할 수 있고, 작은 \(Y_i\) 에서는 차이가 두드러진다.
3.9 KM vs NA 비교표
| 측면 | Kaplan-Meier \(\widehat{S}\) | Nelson-Aalen \(\widetilde{S} = e^{-\widetilde{H}}\) |
|---|---|---|
| 자연성 | 생존함수 직접 | 누적위험 → 변환 |
| 점프 처리 | \(1 - d_i/Y_i\) (multiplicative) | \(-d_i/Y_i\) (additive log) |
| 분산 식 | Greenwood (delta method 두 번) | Aalen (Poisson 근사 한 번) |
| 분산 형태 | \(\widehat{S}^2 \sum d_i / [Y_i(Y_i - d_i)]\) | \(\sum d_i / Y_i^2\) |
| 소표본 | 음의 편향 | \(\widehat{S}\) 보다 큼 |
| 점근적 | 동치 | 동치 |
| 모형 식별 | 약함 | \(\widetilde{H}\) vs \(t\) 직선이면 Exp 분포 |
| Hazard 추정 | 어려움 (점프 크기 비균일) | 점프 크기 \(d_i/Y_i\) 가 \(h\) 의 거친 추정 |
- 생존곡선 보고: KM (의학·임상시험에서 표준).
- 분포 식별 plot: NA — \(\widetilde{H}(t)\) vs \(t\) 가 직선이면 Exponential, \(\log \widetilde{H}\) vs \(\log t\) 가 직선이면 Weibull (Ch.12).
- Hazard 모양 진단: NA 의 점프 크기 \(d_i/Y_i\) 를 시간으로 평활 (kernel smoothing — Ch.6).
- 소표본 (n < 30 사건): NA 가 약간 더 안정적 (Klein 1991).
4 Klein Example 4.1 — 6-MP 21 명 손풀이
4.1 데이터 설명
Freireich (1963) 의 6-MP 군 21 명 (Klein § 1.2 의 Leukemia 임상시험 처치군). 단위 = 주.
관측 데이터 (\(+\) = censoring): 6, 6, 6, 6+, 7, 9+, 10, 10+, 11+, 13, 16, 17+, 19+, 20+, 22, 23, 25+, 32+, 32+, 34+, 35+.
- 사건 (재발) 9 건, censoring 12 건.
- 서로 다른 사건 시점 \(D = 7\): \(t_1=6, t_2=7, t_3=10, t_4=13, t_5=16, t_6=22, t_7=23\).
- \(t_1 = 6\) 에서 tie 3 건.
4.2 KM 계산표 (Klein Table 4.1A·4.1B 재구성)
| \(i\) | \(t_i\) | \(d_i\) | \(Y_i\) | \(1 - d_i/Y_i\) | \(\widehat{S}(t_i)\) | \(\dfrac{d_i}{Y_i(Y_i-d_i)}\) | \(\sum\) | \(\widehat{V}[\widehat{S}]\) | \(\text{SE}\) |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 6 | 3 | 21 | 0.857 | 0.857 | \(\dfrac{3}{21 \cdot 18} = 0.00794\) | 0.00794 | \(0.857^2 \cdot 0.00794 = 0.00583\) | 0.076 |
| 2 | 7 | 1 | 17 | 0.941 | 0.807 | \(\dfrac{1}{17 \cdot 16} = 0.00368\) | 0.01162 | \(0.807^2 \cdot 0.01162 = 0.00756\) | 0.087 |
| 3 | 10 | 1 | 15 | 0.933 | 0.753 | \(\dfrac{1}{15 \cdot 14} = 0.00476\) | 0.01638 | \(0.00929\) | 0.096 |
| 4 | 13 | 1 | 12 | 0.917 | 0.690 | \(\dfrac{1}{12 \cdot 11} = 0.00758\) | 0.02396 | \(0.01140\) | 0.107 |
| 5 | 16 | 1 | 11 | 0.909 | 0.628 | \(\dfrac{1}{11 \cdot 10} = 0.00909\) | 0.03305 | \(0.01303\) | 0.114 |
| 6 | 22 | 1 | 7 | 0.857 | 0.538 | \(\dfrac{1}{7 \cdot 6} = 0.02381\) | 0.05687 | \(0.01643\) | 0.128 |
| 7 | 23 | 1 | 6 | 0.833 | 0.448 | \(\dfrac{1}{6 \cdot 5} = 0.03333\) | 0.09020 | \(0.01807\) | 0.135 |
\(i = 1\) 행 (시점 \(t_1 = 6\) 주):
- \(Y_1 = 21\) — 시작 시점에서 21 명 모두 위험집합.
- \(d_1 = 3\) — 6 주에 3 명 재발.
- \(1 - d_1/Y_1 = 1 - 3/21 = 18/21 = 0.857\).
- \(\widehat{S}(6) = 1 \cdot 0.857 = 0.857\).
- Greenwood term: \(3/(21 \cdot 18) = 0.00794\).
- \(\widehat{V}[\widehat{S}(6)] = 0.857^2 \cdot 0.00794 = 0.00583\).
- SE \(= \sqrt{0.00583} = 0.076\).
\(i = 2\) 행 (시점 \(t_2 = 7\) 주):
- \(Y_2 = 17\) — 21 명 중 사라진 사람은: 6 주 사망 3 명 + 6+ censoring 1 명 = 4 명. 남은 17 명.
- \(d_2 = 1\) — 7 주에 1 명 재발.
- \(\widehat{S}(7) = \widehat{S}(6) \cdot (1 - 1/17) = 0.857 \cdot 0.941 = 0.807\).
- Greenwood: 누적 0.01162.
\(i = 3\) 행 (시점 \(t_3 = 10\) 주):
- \(Y_3 = 15\) — 7 주 후 censoring (9+) 1 명 → 17 - 1 - 1 = 15.
- 그렇다. \(\widehat{S}(10) = 0.807 \cdot 14/15 = 0.753\).
\(i = 4\) 행 (시점 \(t_4 = 13\) 주):
- \(Y_4 = 12\) — 10 주 후 censoring (10+, 11+) 2 명 → 15 - 1 - 2 = 12.
\(i = 5\) 행 (시점 \(t_5 = 16\) 주):
- \(Y_5 = 11\) — 13 주 후 censoring 0 명 → 12 - 1 = 11.
\(i = 6\) 행 (시점 \(t_6 = 22\) 주):
- \(Y_6 = 7\) — 16 주 후 censoring (17+, 19+, 20+) 3 명 → 11 - 1 - 3 = 7. \(Y\) 의 갑작스러운 감소 — censoring 이 많이 발생.
\(i = 7\) 행 (시점 \(t_7 = 23\) 주):
- \(Y_7 = 6\) — 22 주 후 censoring 0 명 → 7 - 1 = 6.
\(t > 35\) 영역: 마지막 관측 35+ 가 censoring 이므로 KM 정의 범위 외 (Practical Note 2 의 tail 처리 필요).
4.3 NA 계산표 (Klein Table 4.2 재구성)
| \(i\) | \(t_i\) | \(d_i\) | \(Y_i\) | \(d_i/Y_i\) | \(\widetilde{H}(t_i)\) | \(d_i/Y_i^2\) | \(\sigma_H^2\) | \(\sigma_H\) |
|---|---|---|---|---|---|---|---|---|
| 1 | 6 | 3 | 21 | 0.1428 | 0.1428 | \(3/441 = 0.0068\) | 0.0068 | 0.0825 |
| 2 | 7 | 1 | 17 | 0.0588 | 0.2017 | \(1/289 = 0.0035\) | 0.0103 | 0.1015 |
| 3 | 10 | 1 | 15 | 0.0667 | 0.2683 | \(1/225 = 0.0044\) | 0.0147 | 0.1212 |
| 4 | 13 | 1 | 12 | 0.0833 | 0.3517 | \(1/144 = 0.0069\) | 0.0217 | 0.1473 |
| 5 | 16 | 1 | 11 | 0.0909 | 0.4426 | \(1/121 = 0.0083\) | 0.0299 | 0.1729 |
| 6 | 22 | 1 | 7 | 0.1429 | 0.5854 | \(1/49 = 0.0204\) | 0.0503 | 0.2243 |
| 7 | 23 | 1 | 6 | 0.1667 | 0.7521 | \(1/36 = 0.0278\) | 0.0781 | 0.2795 |
\(\widetilde{H}\) 의 값: * 0 → 0.143 → 0.202 → 0.268 → 0.352 → 0.443 → 0.585 → 0.752.
차이 (점프 크기): 0.143, 0.059, 0.067, 0.083, 0.091, 0.143, 0.167.
- 시점 6 의 점프 0.143 은 tie 3 건 때문 (이 한 시점에서만 사건 3 건).
- 그 외 점프는 0.06~0.17 범위, 점진적으로 증가.
- \(\widetilde{H}(t)\) vs \(t\) plot 의 기울기가 거의 일정 → constant hazard 의심 → Exponential 분포 후보.
기울기로부터 hazard 의 거친 추정: \(\widetilde{H}(23) - \widetilde{H}(6) = 0.752 - 0.143 = 0.609\). \(23 - 6 = 17\) 주. → \(\widehat{h} \approx 0.609/17 = 0.036\) /week.
(Ch.3 Ex 3.4 의 Exponential MLE \(\widehat{\lambda} = r/S_T = 9/359 = 0.025\) /week 와 같은 자릿수, 완벽히 일치하지는 않음 — KM 은 사건 시점 사이의 censoring 정보만 사용, MLE 는 모든 person-time 사용.)
4.4 KM ↔︎ NA 비교 (Klein Figure 4.1A·4.1B)
| \(t\) | KM \(\widehat{S}\) | NA \(\widetilde{S} = e^{-\widetilde{H}}\) | 차이 |
|---|---|---|---|
| 6 | 0.857 | \(e^{-0.143} = 0.867\) | 0.010 |
| 10 | 0.753 | \(e^{-0.268} = 0.765\) | 0.012 |
| 16 | 0.628 | \(e^{-0.443} = 0.642\) | 0.014 |
| 22 | 0.538 | \(e^{-0.585} = 0.557\) | 0.019 |
| 23 | 0.448 | \(e^{-0.752} = 0.471\) | 0.023 |
→ NA 추정 survival 이 KM 보다 일관되게 크다 (Taylor 1차의 부호와 일치). 차이는 시간이 갈수록 커지는데, 이는 위험집합이 작아져 \(d_i/Y_i\) 가 커지기 때문 (Taylor 근사의 오차가 커짐).
5 Klein Example 4.2 — BMT 137 명
5.1 데이터 설명
Copelan (1991) 의 137 명 골수이식 (BMT) 데이터 (Klein § 1.3). 3 개 disease group:
| Group | \(n\) | 가장 긴 관측 (일) | 사건 정의 |
|---|---|---|---|
| ALL (Acute Lymphoblastic Leukemia) | 38 | 2081 | 재발 또는 사망 (whichever first) |
| AML low risk | 54 | 2569 | 동일 |
| AML high risk | 45 | 2640 | 동일 |
Disease-free survival = “재발도 사망도 안 함”. \(\delta_3 = \max(\delta_1, \delta_2)\) (\(\delta_1\) 재발, \(\delta_2\) 사망).
5.2 ALL 군 24 행 KM·NA 계산 (Klein Table 4.3 발췌)
| \(t_i\) | \(d_i\) | \(Y_i\) | \(\widehat{S}(t_i)\) | \(\sqrt{\widehat{V}[\widehat{S}]}\) | \(\widetilde{H}(t_i)\) | \(\sigma_H\) |
|---|---|---|---|---|---|---|
| 1 | 1 | 38 | 0.9737 | 0.0260 | 0.0263 | 0.0263 |
| 55 | 1 | 37 | 0.9474 | 0.0362 | 0.0533 | 0.0377 |
| 74 | 1 | 36 | 0.9211 | 0.0437 | 0.0811 | 0.0468 |
| 86 | 1 | 35 | 0.8947 | 0.0498 | 0.1097 | 0.0549 |
| 104 | 1 | 34 | 0.8684 | 0.0548 | 0.1391 | 0.0623 |
| 107 | 1 | 33 | 0.8421 | 0.0592 | 0.1694 | 0.0692 |
| 109 | 1 | 32 | 0.8158 | 0.0629 | 0.2007 | 0.0760 |
| 110 | 1 | 31 | 0.7895 | 0.0661 | 0.2329 | 0.0825 |
| 122 | 2 | 30 | 0.7368 | 0.0714 | 0.2996 | 0.0950 |
| 129 | 1 | 28 | 0.7105 | 0.0736 | 0.3353 | 0.1015 |
| 172 | 1 | 27 | 0.6842 | 0.0754 | 0.3723 | 0.1081 |
| 192 | 1 | 26 | 0.6579 | 0.0770 | 0.4108 | 0.1147 |
| 194 | 1 | 25 | 0.6316 | 0.0783 | 0.4508 | 0.1215 |
| 230 | 1 | 23 | 0.6041 | 0.0795 | 0.4943 | 0.1290 |
| 276 | 1 | 22 | 0.5767 | 0.0805 | 0.5397 | 0.1368 |
| 332 | 1 | 21 | 0.5492 | 0.0812 | 0.5873 | 0.1449 |
| 383 | 1 | 20 | 0.5217 | 0.0817 | 0.6373 | 0.1532 |
| 418 | 1 | 19 | 0.4943 | 0.0819 | 0.6900 | 0.1620 |
| 466 | 1 | 18 | 0.4668 | 0.0818 | 0.7455 | 0.1713 |
| 487 | 1 | 17 | 0.4394 | 0.0815 | 0.8044 | 0.1811 |
| 526 | 1 | 16 | 0.4119 | 0.0809 | 0.8669 | 0.1916 |
| 609 | 1 | 14 | 0.3825 | 0.0803 | 0.9383 | 0.2045 |
| 662 | 1 | 13 | 0.3531 | 0.0793 | 1.0152 | 0.2185 |
| 2081 | 0 | 1 | 0.3531 | 0.0793 | 1.0152 | 0.2185 |
5.3 3 년 (1095 일) DFS 비교
마지막 사건 시점 662 일을 지나 1095 일까지는 사건이 없으므로 \(\widehat{S}\) 가 0.3531 그대로 유지. 따라서
| Group | \(\widehat{S}(1095)\) | SE | 95% linear CI |
|---|---|---|---|
| ALL | 0.3531 | 0.0793 | (0.198, 0.508) |
| AML low risk | 0.5470 | 0.0691 | (0.412, 0.682) |
| AML high risk | 0.2444 | 0.0641 | (0.119, 0.370) |
- AML low risk 가 가장 좋은 예후 — 3 년 DFS 약 55%.
- AML high risk 가 가장 나쁨 — 3 년 DFS 약 24%.
- ALL 은 중간 — 3 년 DFS 약 35%.
CI 가 겹치는지 확인: ALL 과 AML high 의 CI 가 약간 겹침 (0.198~0.508 vs 0.119~0.370). AML low 와 다른 두 군은 명확히 분리. → 공식적인 비교는 log-rank 검정 (Ch.7) 으로 확정 필요.
5.4 NA 누적위험의 직선성 (Klein Figure 4.3)
ALL 군의 \(\widetilde{H}\): * \(t = 100\): 0.110. * \(t = 700\): 1.015.
기울기 \(\approx 0.0014\) /day = 0.043 /month → “ALL 환자의 hazard 는 첫 2 년간 약 4.3% /month 정도 일정”.
AML low 의 \(\widetilde{H}\) 기울기 \(\approx 0.020\) /month, AML high 의 기울기 \(\approx 0.060\) /month. 모두 거의 직선 → constant hazard 가정 (Exponential) 이 첫 2 년간 합리적.
6 Practical Notes (8 개)
KM·NA 의 모든 결과는 censoring 의 비정보성에 의존. 실무 점검:
- 임상시험 dropout: 부작용으로 약 끊은 환자가 나쁜 예후를 의미하면 비정보성 위반.
- 의료기록 retrospective: 의료시설 변경으로 추적 끊김은 보통 비정보성 (랜덤).
- 품질 관리: 시험 종료 시점이 부품 상태와 무관하면 비정보성.
위반 의심 시 sensitivity analysis: dropout 환자를 모두 사건 발생 / 모두 사건 안 발생 으로 가정하고 KM 양 끝을 비교 (Klein 2003 Ch.4).
마지막 관측이 censoring 일 때 \(\widehat{S}(t)\) 를 \(t_{\max}\) 너머로 정의하는 방법:
| 방법 | 가정 | \(\widehat{S}(t > t_{\max})\) | 편향 |
|---|---|---|---|
| Efron (1967) | 마지막 censored = 즉시 사건 | 0 | 음 (under) |
| Gill (1980) | 마지막 censored = \(\infty\) | \(\widehat{S}(t_{\max})\) (constant) | 양 (over) |
| Brown-Hollander-Kowar (BHK 1974) | exponential tail | \(\exp\{t \ln[\widehat{S}(t_{\max})]/t_{\max}\}\) | 모수 가정 추가 |
Klein (1991) 의 small-sample 실험 결과: Gill 안이 표준. 이유: Efron 은 “마지막 사람이 즉시 죽음” 이라는 비현실적 가정, BHK 는 exponential 모수 가정 추가.
BMT ALL 군 BHK tail: \(t_{\max} = 2081\) 에서 \(\widehat{S} = 0.3531\) → \(\widehat{S}(t > 2081) = \exp(-0.0005 t)\).
Moeschberger-Klein (1985) 는 BHK 를 Weibull 분포로 확장 (더 유연한 tail).
Greenwood 외에 또 다른 분산 추정량 (Aalen & Johansen 1978):
\[ \widetilde{V}[\widehat{S}(t)] = \widehat{S}(t)^2 \sum_{t_i \leq t} \frac{d_i}{Y_i^2} \]
(Greenwood 의 \(Y_i(Y_i - d_i)\) 가 \(Y_i^2\) 로 대체).
Klein (1991) 비교:
- 두 추정량 모두 소표본에서 진짜 분산을 과소 추정.
- Greenwood 가 평균적으로 진짜 값에 더 가까움.
- 단 \(Y_1\) 이 매우 작으면 Aalen-Johansen 이 더 작은 분산.
관행: Greenwood 를 default 로. R survival::survfit 의 error="greenwood".
NA 의 Aalen 분산 (식 4.2.4) 외 대안:
\[ \widehat{V}[\widetilde{H}(t)] = \sum_{t_i \leq t} \frac{(Y_i - d_i) d_i}{Y_i^3} \]
이 추정량은 진짜 분산을 과소 추정, 식 4.2.4 는 과대 추정.
Klein (1991): 식 4.2.4 가 편향이 더 작음 → 4.2.4 권장.
직관: 식 4.2.6 은 finite-population correction 적용 (Binomial 의 hypergeometric 보정). 그러나 점근적으로 4.2.4 가 더 robust.
\(\widetilde{S} = \exp(-\widetilde{H})\) 의 분산은 식 4.2.2 또는 4.2.5 에서 \(\widehat{S}\) 를 \(\widetilde{S}\) 로 대체. Delta method:
\[ \text{Var}[\widetilde{S}] = \widetilde{S}^2 \cdot \text{Var}[\widetilde{H}] \]
(왜냐하면 \(\frac{d \exp(-H)}{d H} = -\exp(-H) = -\widetilde{S}\)).
모든 \(\delta_i = 1\) 이면 \(Y_i = n - i + 1\), \(d_i = 1\) (각 시점에 1 건). KM:
\[ \widehat{S}(t) = \prod_{t_i \leq t}\left(1 - \frac{1}{n-i+1}\right) = \frac{n - i}{n} \]
이는 empirical survival \(1 - i/n\). → censoring 없는 한계에서 KM 이 empirical CDF 로 환원. 이 일관성이 KM 의 자명한 성질.
| 패키지 | KM | NA |
|---|---|---|
R survival |
survfit(Surv(time, status) ~ 1) |
survfit(..., type="fh") (Fleming-Harrington) |
R KMsurv |
KM 계산 헬퍼 | — |
Python lifelines |
KaplanMeierFitter |
NelsonAalenFitter |
Python scikit-survival |
kaplan_meier_estimator |
nelson_aalen_estimator |
| SAS | PROC LIFETEST |
PROC LIFETEST METHOD=NA |
| SPSS | Survival > Kaplan-Meier | (수동) |
| Stata | sts list, survival |
sts list, na |
R survival 의 default 분산은 Greenwood. survfit(..., error="tsiatis") 는 Tsiatis 분산 (Aalen-Johansen 변형).
7 Theoretical Notes (9 개)
이산 분포 + 사건 시점에만 점프 가정. 조건부 확률 사슬 분해:
\[ S(t_i) = \prod_{j \leq i} P(X > t_j \mid X \geq t_j) \]
\(\widehat{P}(X > t_j \mid X \geq t_j) = (Y_j - d_j)/Y_j\). 곱하면 KM. 출발점이 가장 자명한 유도.
각 censoring mass 를 우측 관측에 균등 분배. 알고리즘이 KM 으로 수렴함을 Efron (1967) 가 증명.
의의: KM 을 “censoring 정보를 어떻게 분배할지” 라는 직관적 알고리즘으로 도출 — 비통계학자도 이해 가능.
식 4.2.7 의 fixed-point 방정식. 임의 초기값에서 시작하여 반복 → KM 으로 수렴.
EM 알고리즘과의 관계: censored 관측의 missing event time 의 conditional expectation 을 사용. → 비모수 EM 의 원형. Dempster-Laird-Rubin (1977) 의 EM 일반화 이전에 이미 KM 의 도출에 등장.
현대적 의의: interval censoring 의 Turnbull (Ch.5), 좌절단의 Andersen et al., 일반 missing data 모형의 EM — 모두 self-consistency 의 변형.
Andersen-Borgan-Gill-Keiding (1993) 의 통일 framework (Ch.3.6 참조). KM = \(\prod(1 - dN/Y)\) product integral, NA = \(\int J/Y \, dN\) stochastic integral.
의의: 점근 정규성·신뢰대·log-rank 등 모든 추론을 martingale CLT 로 일관되게 도출. Cox model (Ch.8), Aalen additive model (Ch.10), Schoenfeld residuals (Ch.11) 의 수학적 토대.
regularity 조건 하에서 KM·NA 는 nonparametric MLE.
증명 스케치 (위의 유도 5):
- 점프가 사건 시점에만 있는 분포족.
- likelihood 가 각 \(h_i\) 의 Binomial 로 분해.
- \(\widehat{h}_i = d_i/Y_i\) → \(\widehat{S}, \widetilde{H}\).
왜 중요한가: NPMLE 라는 사실이 KM 의 robustness 의 깊은 이유. 분포 가족을 어떻게 잡든 (사건 시점에 점프 가능한 모든 이산 분포), likelihood 의 최대값은 KM.
regularity 조건 하에서:
- \(\widehat{S}(t) \to S(t)\) (in probability) as \(n \to \infty\).
- \(\widetilde{S}(t) \to S(t)\) (in probability).
- \(\widehat{S}(t) - \widetilde{S}(t) = O_p(1/n)\) (점근적 동치).
→ 큰 표본에서 KM 과 NA 중 어느 것을 써도 같은 결론.
\(\log(1-x) = -x - x^2/2 - \cdots\):
\[ -\log \widehat{S}(t) = \sum -\log(1 - d_i/Y_i) = \sum (d_i/Y_i + (d_i/Y_i)^2/2 + \cdots) \]
1차 항 = \(\widetilde{H}\), 2차 이상은 \(O(1/Y_i)\). \(Y_i \to \infty\) 면 차이 사라짐.
소표본 함의: 작은 \(Y_i\) 에서 \(-\log \widehat{S} \neq \widetilde{H}\), 즉 KM 과 NA 가 실질적으로 다름. 후반부 (위험집합 작음) 에서 둘의 차이 두드러짐.
- Guerts (1987): KM 의 소표본 성질 — 음의 편향, 표본 크기 증가에 따른 수렴 속도.
- Wellner (1985): KM 의 점근적 분포 + 신뢰대.
- Klein (1991): KM·NA 분산 추정량의 비교 — Greenwood 가 평균적으로 가장 가까움, 식 4.2.4 가 NA 분산으로 권장.
실무 함의: \(n < 30\) 사건에서는 KM·NA 의 차이를 의식적으로 보고 (둘 다 plot), 분산은 Greenwood 사용.
regularity 조건 하에서:
\[ \sqrt{n}(\widehat{S}(t) - S(t)) \xrightarrow{d} W^*(t) \]
where \(W^*\) is a mean-zero Gaussian process with covariance \(\sigma^2(s, t) = S(s)S(t) \int_0^{\min(s,t)} dF/[(1-F)^2 G^{-1}]\).
의의:
- 점별 정규성: 고정 \(t\) 에서 \(\widehat{S}(t) \approx N(S(t), \widehat{V}[\widehat{S}(t)])\) — § 4.3 의 신뢰구간의 기반.
- 신뢰대: process 자체의 sup 분포로 § 4.4 의 confidence band 도출.
- Brownian bridge 표현: 표준화 후 Brownian bridge 로 수렴 — Hall-Wellner band 의 기반.
→ § 4.3·4.4 의 모든 신뢰구간/대 가 본 정리 위에서 작동.
8 응용 분야
| 분야 | 활용 | 구체적 예시 |
|---|---|---|
| 임상시험 | 일차 종점 KM | 신약 vs 위약의 OS·PFS 비교 |
| 종양학 | 5 년 생존률 | KM 5-년 OS = 65% (95% CI 58~72) |
| 신뢰성 공학 | 제품 수명 | NA 기울기로 Weibull/Exp 식별 |
| HIV 역학 | 잠복기 분포 | reverse-time KM (Ch.5) 이전의 우측 censoring |
| 보험·금융 | 부도 시점 분포 | 채무자 default 곡선 |
| 인구 통계 | Cohort 사망률 | 좌절단 보정 KM (§ 4.6) |
| 재발 모형 | 첫 재발 시점 | KM 의 사건 정의를 재발로 |
9 코드 예시
9.1 Step 1 — 순수 NumPy 로 KM·NA 직접 구현
import numpy as np
# Klein Example 4.1: 6-MP 21 명
times = np.array([6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16,
17, 19, 20, 22, 23, 25, 32, 32, 34, 35])
events = np.array([1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1,
0, 0, 0, 1, 1, 0, 0, 0, 0, 0])
def km_na(times, events):
"""Kaplan-Meier + Nelson-Aalen + Greenwood + Aalen variance"""
order = np.argsort(times)
t_sorted = times[order]
e_sorted = events[order]
unique_event_times = np.unique(t_sorted[e_sorted == 1])
S = 1.0 # KM
H = 0.0 # NA
var_S_term = 0.0 # Greenwood 누적합
var_H = 0.0 # Aalen 누적합
rows = []
for ti in unique_event_times:
Y = int(np.sum(t_sorted >= ti)) # 위험집합
d = int(np.sum((t_sorted == ti) & (e_sorted == 1))) # 사건 수
# KM (식 4.2.1)
S *= (1 - d / Y)
# NA (식 4.2.3)
H += d / Y
# Greenwood (식 4.2.2): 누적 후 S^2 곱
var_S_term += d / (Y * (Y - d)) if (Y - d) > 0 else 0
var_S = (S ** 2) * var_S_term
# Aalen (식 4.2.4)
var_H += d / (Y ** 2)
rows.append((ti, Y, d, S, np.sqrt(var_S),
H, np.sqrt(var_H)))
return rows
print(f"{'t':>3} {'Y':>3} {'d':>2} {'S':>6} {'SE_S':>6} "
f"{'H':>6} {'SE_H':>6}")
for ti, Y, d, S, se_S, H, se_H in km_na(times, events):
print(f"{ti:>3} {Y:>3} {d:>2} {S:.3f} {se_S:.3f} "
f"{H:.3f} {se_H:.3f}")출력 (Klein Table 4.1B + Table 4.2 와 일치):
t Y d S SE_S H SE_H
6 21 3 0.857 0.076 0.143 0.082
7 17 1 0.807 0.087 0.202 0.102
10 15 1 0.753 0.096 0.268 0.121
13 12 1 0.690 0.107 0.352 0.147
16 11 1 0.628 0.114 0.443 0.173
22 7 1 0.538 0.128 0.585 0.224
23 6 1 0.448 0.135 0.752 0.2809.2 Step 2 — Self-consistency 알고리즘 직접 구현
식 4.2.7 의 fixed-point 반복으로 KM 도출.
def self_consistent_km(times, events, max_iter=100, tol=1e-8):
"""Efron (1967) self-consistency 알고리즘"""
n = len(times)
grid = np.unique(times)
# 초기값: empirical (censored 무시)
S = np.array([(times > t).mean() for t in grid])
for iteration in range(max_iter):
S_new = np.zeros_like(S)
for k, t in enumerate(grid):
term1 = np.sum(times > t) / n
mask = (events == 0) & (times <= t)
term2 = 0.0
for i in np.where(mask)[0]:
idx_Ti = np.searchsorted(grid, times[i])
if S[idx_Ti] > 0:
term2 += S[k] / S[idx_Ti] / n
S_new[k] = term1 + term2
if np.max(np.abs(S_new - S)) < tol:
print(f"수렴: iteration {iteration}")
break
S = S_new
return grid, S
grid, S_self = self_consistent_km(times, events)
# KM 결과와 동일함을 확인9.3 Step 3 — lifelines 로 실무 코드
import pandas as pd
from lifelines import KaplanMeierFitter, NelsonAalenFitter
df = pd.DataFrame({"t": times, "event": events})
kmf = KaplanMeierFitter().fit(df["t"], df["event"], label="6-MP")
naf = NelsonAalenFitter().fit(df["t"], df["event"], label="6-MP")
# 핵심 출력
print(kmf.survival_function_) # KM 곡선
print(kmf.confidence_interval_) # 점별 95% CI (default)
print(naf.cumulative_hazard_) # NA 곡선
# 두 추정량 비교
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
kmf.plot_survival_function(ax=axes[0], label="KM")
S_from_NA = np.exp(-naf.cumulative_hazard_)
S_from_NA.plot(ax=axes[0], label="exp(-NA)", linestyle="--")
axes[0].set_title("KM vs exp(-NA)")
axes[0].legend()
naf.plot_cumulative_hazard(ax=axes[1])
H_from_KM = -np.log(kmf.survival_function_)
H_from_KM.plot(ax=axes[1], label="-log(KM)", linestyle="--")
axes[1].set_title("NA vs -log(KM)")
axes[1].legend()9.4 Step 4 — R survival 패키지
library(survival)
library(KMsurv)
# 6-MP 데이터
sixmp <- data.frame(
time = c(6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16,
17, 19, 20, 22, 23, 25, 32, 32, 34, 35),
status = c(1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1,
0, 0, 0, 1, 1, 0, 0, 0, 0, 0)
)
# KM (default Greenwood variance)
km_fit <- survfit(Surv(time, status) ~ 1, data = sixmp)
summary(km_fit)
# Time N Events Survival StdErr Lower 95% CI Upper 95% CI
# 6 21 3 0.857 0.0764 0.7194 1.0000
# 7 17 1 0.807 0.0869 0.6531 0.9971
# 10 15 1 0.753 0.0963 0.5859 0.9669
# ...
# NA (Fleming-Harrington type = NA-based survival)
na_fit <- survfit(Surv(time, status) ~ 1, data = sixmp, type = "fh")
summary(na_fit)
# Aalen-Johansen variance (식 4.2.5)
km_aj <- survfit(Surv(time, status) ~ 1, data = sixmp,
error = "tsiatis")
summary(km_aj)
# 시각화
plot(km_fit, conf.int = TRUE,
xlab = "Weeks", ylab = "Survival Probability",
main = "Kaplan-Meier — 6-MP Group (Klein Ex 4.1)")9.5 Step 5 — BMT 137 명 그룹별 비교
from lifelines.datasets import load_bmt
df = load_bmt() # patients with bone marrow transplants
# disease group 별 KM
fig, ax = plt.subplots(figsize=(10, 6))
for group, label in zip([1, 2, 3], ["ALL", "AML low", "AML high"]):
mask = df["group"] == group
kmf = KaplanMeierFitter().fit(df.loc[mask, "T"],
df.loc[mask, "event"],
label=label)
kmf.plot_survival_function(ax=ax, ci_show=True)
ax.set_title("BMT — Disease-Free Survival by Group")
ax.set_xlabel("Days Post Transplant")
ax.set_ylabel("DFS Probability")
# 3 년 (1095 일) 시점 출력
for group, label in zip([1, 2, 3], ["ALL", "AML low", "AML high"]):
mask = df["group"] == group
kmf = KaplanMeierFitter().fit(df.loc[mask, "T"],
df.loc[mask, "event"])
print(f"{label}: S(1095) = {kmf.predict(1095):.4f}")10 핵심 takeaway
\(d_i/Y_i\) 가 모든 비모수 추정의 building block — 이 한 가지 비율의 곱이 KM, 합이 NA, 적분이 RMST, 역변환이 median, 재정의가 left-truncated KM, 분리가 cause-specific hazard. 도구 7 개가 한 가지 양의 변형이다.
5 가지 유도가 모두 같은 결과 — reduced-sample, redistribute-to-the-right, self-consistency, counting process, NPMLE 가 모두 KM 으로 모인다. 이 다중성이 KM 의 robustness 의 수학적 깊이.
Greenwood 분산은 delta method 두 번 — log 변환 후 분산 합산, 다시 exp 변환. Aalen 분산은 Poisson 근사 한 번. 두 도출이 서로 다른 가정에서 출발하지만 점근적으로 동치.
KM 과 NA 는 점근적으로 동치이지만 소표본에서 다르다 — NA = -log KM 의 Taylor 1차. 작은 \(Y_i\) 에서 \(\widetilde{S} > \widehat{S}\). Hazard 모양 진단은 NA, 생존곡선 보고는 KM.
Tail 처리는 가정에 의존 — Efron (immediate death), Gill (immortal), BHK (exponential decay). Klein (1991) 의 small-sample 비교 결과 Gill 권장. 분석 보고에 어떤 안을 썼는지 명시.
독립 censoring 가정 위반은 추정 대상을 바꾼다 — informative dropout 시 KM 이 다른 함수를 추정. sensitivity analysis (worst-case bounds) 로 점검 필수.
11 관련 주제
선행 지식
- Ch.4 Overview — 7 절 조망
- Ch.2 — 4 함수의 동등성 + 9 parametric models
- § 3.6 — Counting Process (Aalen 1975)
- § 3.5 — Likelihood master 식 (KM·NA 의 likelihood 출발점)
후속 주제
- § 4.3 — Pointwise Confidence Intervals (linear · log · arcsine)
- § 4.4 — Confidence Bands (EP · Hall-Wellner)
- § 4.5 — Mean (RMST) · Median (Brookmeyer-Crowley)
- § 4.6 — Left-Truncation 의 \(Y_i\) 재정의
- Ch.6 — Kernel Hazard Smoothing (NA 점프의 평활)
- Ch.7 — Log-rank, Wilcoxon, Tarone-Ware (KM 의 두 곡선 비교)
- Ch.12 — Parametric Model Identification (NA plot 으로 분포 식별)
관련 개념