Klein § 8.7~8.8 — Model Building Strategy · Survival Function Estimation

1 들어가며 — Ch.8 네 번째 deep-dive

편	주제
Ch.8 Overview	9 절 조망
§ 8.1~8.2	Cox 모형 + Coding
§ 8.3~8.4	Partial Likelihood + Ties
§ 8.5~8.6	Local Tests + Discretizing
§ 8.7~8.8 (본 편)	Model Building + Survival Estimation
§ 8.9	Exercises

§ 8.7~8.8 의 한 줄 요약

“§ 8.7 — 모형 구축의 두 시나리오. (1) 가설 주도: 주효과 고정 + 교란 순차 탐색. (2) 탐색적: 전진 선택 + AIC \(= -2\log L + 2p\). § 8.8 — Breslow 추정량 \(\widehat{H}_0(t) = \sum d_i / W(t_i; b)\) (NA 의 일반화) → \(\widehat{S}(t|Z_0) = \widehat{S}_0(t)^{\exp(b'Z_0)}\) (Lehmann alternative 의 직접). 분산 = baseline 불확실 \(Q_1\) + \(\beta\) 추정 불확실 \(Q_2\). Klein Example 8.2: 60세 Stage I 5년 \(\widehat{S} = 0.703\), Stage IV \(\widehat{S} = 0.147\).”

왜 이 두 절을 묶는가

§ 8.7 은 “어떤 공변량을 모형에 넣을 것인가”를 다루고, § 8.8 은 “결정된 모형에서 개별 환자의 예후를 어떻게 예측하는가”를 다룬다. 모형 구축 → 예측이라는 자연스러운 파이프라인의 마지막 두 단계이다. § 8.3~8.6 에서 편우도 추정과 검정 도구를 확보했으므로, 이제 그 도구를 조합하여 실전 모형을 만들고 환자별 생존곡선을 산출하는 것이 본 편의 목표이다.

2 § 8.7 — Model Building Using the Proportional Hazards Model

2.1 두 가지 회귀 시나리오

Cox 모형을 적합하는 상황은 크게 두 가지로 나뉜다. 이 구분이 중요한 이유는 모형 구축 전략이 완전히 달라지기 때문이다.

시나리오 구분

시나리오 1 — 가설 주도(Hypothesis-Driven)

특정 가설이 사전에 존재한다 (예: “치료군 간 생존 차이가 있는가?”).
관심 변수(treatment, risk group)를 모형에 고정하고, 나머지 변수는 교란(confounder) 탐색 용도로만 추가한다.
목표: 교란 보정 후 주효과의 유의성을 평가.

시나리오 2 — 탐색적(Exploratory)

사전 가설 없이 다수의 설명변수 중 생존과 관련된 변수를 찾는다.
전진 선택(forward selection), 후진 제거(backward elimination), 단계적 선택(stepwise)을 사용한다.
목표: 예측 모형 구축 또는 향후 연구를 위한 가설 생성(hypothesis generating).

두 시나리오의 핵심 차이를 직관적으로 비유하면: 시나리오 1은 법정에서 “피고(주효과)가 유죄인가”를 심리할 때 방해 요소(교란)를 제거하는 것이고, 시나리오 2는 수사 단계에서 “용의자 전체 중 누가 의심스러운가”를 좁혀가는 것이다. 법정에서는 피고를 바꾸지 않지만, 수사에서는 가장 의심스러운 인물을 순서대로 추적한다.

2.2 시나리오 1: 가설 주도 교란 탐색

절차는 다음과 같다:

주효과 변수만으로 global test (§ 8.3~8.4)를 수행하여 비보정(unadjusted) 관계를 확인한다.
주효과를 모형에 고정한 채, 나머지 변수 각각을 하나씩 추가하여 local test (§ 8.5)를 수행한다.
가장 유의한 교란 변수를 모형에 추가한다.
주효과 + 추가된 교란을 고정한 채, 남은 변수에 대해 local test를 반복한다.
더 이상 유의한 교란이 없으면 중단한다.

이 절차의 핵심은 주효과 변수가 절대 모형에서 빠지지 않는다는 것이다. 교란 변수가 생존과 무관하면 교란도 아니므로 추가할 필요가 없다.

2.3 시나리오 2: 탐색적 전진 선택

각 설명변수를 단독으로 global test 하여 단변량 관계를 파악한다.
가장 유의한 변수를 모형에 넣고, 나머지 각각에 대해 local test를 수행한다.
가장 유의한 변수를 추가하고 반복한다.
p-value 기준으로 유의하지 않으면 중단한다.

시나리오 1 과 달리 어떤 변수도 사전에 고정되지 않는다. 순수하게 데이터가 변수 선택을 주도한다.

2.4 AIC (Akaike Information Criterion)

정의: AIC

\[ \text{AIC} = -2 \log L + kp \]

\(L\): 편우도(partial likelihood) 값.
\(p\): 모형 내 회귀 모수의 수.
\(k\): 벌칙 상수 (표준: \(k = 2\)).

AIC 의 직관은 “적합도(fit)”와 “복잡도(complexity)”의 균형이다. \(-2\log L\) 은 모형이 데이터를 얼마나 잘 설명하는지를 나타내며 작을수록 좋다. \(kp\) 는 모수를 추가할 때마다 부과하는 벌칙이다. 변수를 추가하면 \(-2\log L\) 은 항상 감소하지만, 벌칙 \(kp\) 가 증가하므로 어느 시점에서 AIC 가 증가로 반전한다. 이 반전 지점이 “불필요한 변수”의 시작이다.

AIC 는 최소제곱 회귀의 adjusted \(R^2\) 와 유사한 역할을 한다. 둘 다 모수 수로 적합도를 보정하여 과적합을 방지하지만, AIC 는 정보이론에 기반하고 adjusted \(R^2\) 는 분산 보정에 기반한다는 차이가 있다.

p-value 접근과 AIC 접근은 같은 데이터에서 다른 최종 모형을 줄 수 있다. 일반적으로 AIC 가 p-value 보다 더 많은 변수를 포함하는 경향이 있다. 이는 AIC 의 벌칙 \(k = 2\) 가 대략 \(\alpha \approx 0.16\) 의 유의수준에 대응하기 때문이다 (Klein & Moeschberger, 2003, Ch.8).

2.5 Klein Example 8.5: BMT Disease-Free Survival (가설 주도)

데이터 배경

§ 1.3 에서 소개된 골수이식(BMT) 데이터이다. 무작위 배정이 아니므로 관찰 연구의 편향을 보정해야 한다.

주효과: 3 risk group — ALL, AML low-risk, AML high-risk.
- \(Z_1 = 1\) (AML low-risk), \(Z_2 = 1\) (AML high-risk). ALL 이 기저(baseline).
교란 후보: 환자 요인 (\(Z_3\): 대기시간, \(Z_4\): FAB class, \(Z_5\): MTX 여부), 환자-공여자 결합 요인 (성별 \(Z_6, Z_7, Z_8\); CMV \(Z_9, Z_{10}, Z_{11}\); 나이 \(Z_{12}, Z_{13}, Z_{14}\)).

Step 1 — 비보정 global test:

Risk group 단독 모형의 Wald \(\chi^2 = 13.01\) (df=2, p=0.001). AIC = 737.29.

직관: 보정 없이도 risk group 간 생존 차이가 강하게 존재한다. 그러나 관찰 연구이므로 이 차이가 치료 효과인지 환자 특성 차이인지 구분할 수 없다. 교란 탐색이 필수이다.

Step 2 — 첫 번째 교란 탐색 (주효과 \(Z_1, Z_2\) 고정):

교란 후보	df	Wald \(\chi^2\)	p	AIC
대기시간 (\(Z_3\))	1	1.18	0.277	737.95
FAB class (\(Z_4\))	1	8.08	0.004	731.02
MTX (\(Z_5\))	1	2.03	0.155	737.35
성별 (\(Z_6, Z_7, Z_8\))	3	1.91	0.591	741.44
CMV (\(Z_9, Z_{10}, Z_{11}\))	3	0.19	0.980	743.10
나이 (\(Z_{12}, Z_{13}, Z_{14}\))	3	11.98	0.007	733.18

FAB class 가 p=0.004 로 가장 유의하고 AIC 도 가장 낮다. 모형에 추가.

직관: FAB class (M4/M5 vs 기타)는 AML 환자의 예후에 직접 영향을 미치는 생물학적 특성이다. 이것이 risk group 과 상관되어 있다면 FAB 을 보정하지 않으면 risk group 효과가 과대 또는 과소 추정된다.

Step 3 — 두 번째 교란 탐색 (\(Z_1, Z_2, Z_4\) 고정):

교란 후보	df	Wald \(\chi^2\)	p	AIC
대기시간 (\(Z_3\))	1	1.18	0.277	731.68
MTX (\(Z_5\))	1	2.05	0.152	731.06
성별 (\(Z_6, Z_7, Z_8\))	3	0.92	0.820	736.11
CMV (\(Z_9, Z_{10}, Z_{11}\))	3	0.02	0.999	737.00
나이 (\(Z_{12}, Z_{13}, Z_{14}\))	3	13.05	0.004	725.98

나이 factor 가 p=0.004 로 유의. 모형에 추가.

직관: 공여자-수혜자 나이 interaction (\(Z_{14} = Z_{12} \times Z_{13}\))이 유의하다는 것은 공여자가 젊고 수혜자가 나이 들수록(또는 그 반대) 위험이 비선형적으로 증가한다는 뜻이다. 이 교호작용을 무시하면 나이의 효과가 평균화되어 과소 추정된다.

Step 4 — 세 번째 교란 탐색 (\(Z_1, Z_2, Z_4, Z_{12}, Z_{13}, Z_{14}\) 고정):

교란 후보	df	Wald \(\chi^2\)	p	AIC
대기시간 (\(Z_3\))	1	0.46	0.495	727.48
MTX (\(Z_5\))	1	1.44	0.229	726.58
성별 (\(Z_6, Z_7, Z_8\))	3	1.37	0.713	730.61
CMV (\(Z_9, Z_{10}, Z_{11}\))	3	0.58	0.902	731.42

모든 local test 가 비유의, AIC 도 이전 모형(725.98)보다 높다. 중단.

최종 ANOVA Table (Table 8.9):

변수	\(b\)	SE(\(b\))	Wald \(\chi^2\)	p
\(Z_1\) (AML low-risk)	\(-1.091\)	0.354	9.48	0.002
\(Z_2\) (AML high-risk)	\(-0.404\)	0.363	1.24	0.265
\(Z_4\) (FAB M4/M5)	0.837	0.279	9.03	0.003
\(Z_{12}\) (donor age \(-28\))	0.004	0.018	0.05	0.831
\(Z_{13}\) (recipient age \(-28\))	0.007	0.020	0.12	0.728
\(Z_{14}\) (\(Z_{12} \times Z_{13}\))	0.003	0.001	11.01	0.001

해석:

\(Z_1\) (AML low-risk vs ALL): \(\text{RR} = e^{-1.091} = 0.34\). AML low-risk 이 ALL 보다 위험이 66% 낮다 (p=0.002). FAB 과 나이를 보정한 후에도 유의하다.
\(Z_2\) (AML high-risk vs ALL): \(\text{RR} = e^{-0.404} = 0.67\). 비유의 (p=0.265). 보정 후 AML high-risk 와 ALL 의 차이가 사라진다.
\(Z_4\) (FAB M4/M5): \(\text{RR} = e^{0.837} = 2.31\). FAB M4/M5 분류 환자는 위험이 2.3배 높다.
\(Z_{14}\): 나이 interaction 이 p=0.001 로 강력 유의. 공여자와 수혜자의 나이 차이가 예후에 비선형적 영향을 미친다.

primary hypothesis (risk group 전체)의 local test: p=0.003. 교란 보정 후에도 risk group 간 생존 차이가 유의하게 존재한다.

2.6 Klein Example 8.6: Weaning Time (탐색적)

데이터 배경

§ 1.14 의 모유수유 이유(weaning) 시간 데이터이다. 동점(ties)이 많으므로 discrete likelihood 를 사용한다.

후보 변수 7 개: 인종(2 df), 빈곤여부(1 df), 흡연(1 df), 음주(1 df), 나이(1 df), 교육(2 df), 산전관리 부재(1 df).

Step 1 — 단변량 global test:

Factor	df	Wald \(\chi^2\)	p	AIC
인종	2	8.03	0.018	5481.67
빈곤	1	0.71	0.399	5486.69
흡연	1	10.05	0.002	5477.61
음주	1	2.01	0.157	5485.48
나이	1	0.15	0.698	5487.26
교육	2	6.95	0.031	5482.36
산전관리 부재	1	0.16	0.687	5487.25

흡연이 p=0.002, AIC 최소. 모형에 추가.

직관: 흡연 모와 이유 시간의 관계 — 흡연 모는 비흡연 모보다 이유를 더 빨리 하는 경향이 있다. \(b = 0.308\) 이므로 \(\text{RR} = e^{0.308} = 1.36\). 흡연 모의 이유 위험이 36% 높다(= 이유 시간이 더 짧다).

Step 2 — local test (흡연 보정):

Factor	df	Wald \(\chi^2\)	p	AIC
인종	2	12.38	0.002	5469.71
빈곤	1	1.42	0.234	5478.17
음주	1	1.04	0.307	5478.59
나이	1	0.01	0.954	5479.61
교육	2	3.87	0.145	5477.71
산전관리 부재	1	0.02	0.888	5479.59

인종이 p=0.002 로 유의. 모형에 추가.

주목할 점: 단변량에서 교육(p=0.031)이 인종(p=0.018)보다 덜 유의했는데, 흡연 보정 후 인종(p=0.002)이 교육(p=0.145)보다 훨씬 유의해진다. 이는 흡연과 교육이 교란되어 있기 때문이다. 흡연을 보정하면 교육의 “독립적” 효과가 크게 줄어든다.

Step 3 — local test (흡연 + 인종 보정):

Factor	df	Wald \(\chi^2\)	p	AIC
빈곤	1	2.99	0.084	5468.65
음주	1	1.16	0.281	5470.58
나이	1	0.19	0.660	5471.51
교육	2	2.08	0.353	5471.60
산전관리 부재	1	0.03	0.854	5471.67

p-value 접근 (\(\alpha = 0.05\)): 모든 변수가 비유의. 최종 모형 = 흡연 + 인종 (2-factor).

변수	\(b\)	SE(\(b\))	Wald \(\chi^2\)	p
흡연	0.308	0.081	14.34	\(<0.001\)
인종-Black	0.156	0.111	1.98	0.159
인종-Other	0.350	0.102	11.75	\(<0.001\)

AIC 접근: 빈곤의 AIC(5468.65)가 현재 모형의 AIC(5469.71)보다 낮으므로 빈곤을 추가. 그 이후 더 이상 AIC 가 감소하지 않아 중단. 최종 모형 = 흡연 + 인종 + 빈곤 (3-factor).

변수	\(b\)	SE(\(b\))	Wald \(\chi^2\)	p
흡연	0.328	0.082	15.96	\(<0.001\)
인종-Black	0.184	0.112	2.70	0.100
인종-Other	0.374	0.103	13.18	\(<0.001\)
빈곤	\(-0.163\)	0.094	2.99	0.084

빈곤의 p=0.084 는 \(\alpha = 0.05\) 에서 비유의하지만, AIC 기준에서는 모형에 포함된다. 빈곤 계수 \(b = -0.163\) 이므로 \(\text{RR} = e^{-0.163} = 0.85\). 빈곤가구 모의 이유 위험이 15% 낮다(= 이유 시간이 더 길다). 이는 경제적 이유로 모유수유를 더 오래 유지하는 경향을 반영할 수 있다.

2.7 Practical Notes

Wald, Score, LR 중 어떤 통계량을 사용해도 무방하다. 고차원 모형에서는 Score 통계량이 효율적이다 — 각 step 에서 고차원 모형을 적합하지 않고도 통계량을 계산할 수 있기 때문이다.
전진 선택의 대안으로 후진 제거(backward elimination) 가 있다. 전체 모형에서 시작하여 가장 비유의한 변수를 순차적으로 제거한다. 단계적 선택(stepwise) 은 전진과 후진을 결합한다.
AIC 의 벌칙 상수 \(k\) 를 키우면 보수적(변수를 적게 포함)이 된다. \(k = \log n\) 이면 BIC(Bayesian Information Criterion)이며, AIC 보다 더 간결한 모형을 선호한다.

3 § 8.8 — Estimation of the Survival Function

3.1 왜 생존함수 추정이 필요한가

Cox 모형의 편우도 추정은 \(\beta\) 만 추정하고, 기저 위험함수 \(h_0(t)\) 는 추정하지 않는다. 그 덕분에 \(h_0\) 의 형태를 가정하지 않는 유연성을 얻었다. 그러나 실무에서는 “이 환자의 5년 생존확률은 얼마인가?”라는 질문에 답해야 한다. 이를 위해서는 \(\beta\) 뿐 아니라 \(h_0(t)\) 의 추정이 필요하다.

직관적으로: 위험비(hazard ratio)는 “A 가 B 보다 2배 위험하다”는 상대적 정보만 준다. 생존확률 추정은 “A 의 5년 생존확률은 70%이다”라는 절대적 정보를 준다. 임상에서 환자에게 예후를 설명하려면 절대적 수치가 필요하다.

3.2 Breslow 추정량: 기저 누적위험

정의: Breslow 누적 기저 위험 추정량

\(t_1 < t_2 < \cdots < t_L\) 을 구별 사건시간, \(d_i\) 를 시간 \(t_i\) 에서의 사건 수라 하자. 가중 위험집합 합을 다음과 같이 정의한다:

\[ W(t_i; b) = \sum_{j \in R(t_i)} \exp(b'Z_j) \tag{8.8.1} \]

기저 누적위험의 Breslow 추정량은

\[ \widehat{H}_0(t) = \sum_{t_i \leq t} \frac{d_i}{W(t_i; b)} \tag{8.8.2} \]

기저 생존함수는

\[ \widehat{S}_0(t) = \exp\left[-\widehat{H}_0(t)\right] \tag{8.8.3} \]

식 8.8.2 의 직관은 Nelson-Aalen 추정량의 자연스러운 일반화이다. 공변량이 없으면 \(W(t_i; b) = \sum_{j \in R(t_i)} \exp(0) = n_i\) (위험집합의 크기)이므로 \(\widehat{H}_0(t) = \sum d_i / n_i\) 로 정확히 Nelson-Aalen 이 된다.

공변량이 있을 때의 차이: \(W(t_i; b)\) 는 “가중” 위험집합 크기이다. 위험이 높은 사람(\(\exp(b'Z_j)\) 가 큰 사람)이 위험집합에 있으면 \(W\) 가 커지고, 기저 위험 기여분 \(d_i / W\) 가 작아진다. 직관적으로, 사건이 발생했지만 위험집합에 고위험 환자가 많으면 “기저 수준에서는 그리 위험하지 않았다”고 해석하는 것이다.

3.3 공변량 보정 생존함수

정의: 공변량 보정 생존 추정

공변량 벡터 \(Z = Z_0\) 인 개인의 생존함수 추정은

\[ \widehat{S}(t | Z = Z_0) = \widehat{S}_0(t)^{\exp(b'Z_0)} \tag{8.8.4} \]

이 식은 Ch.2 에서 소개한 Lehmann alternative, 즉 \(S(t|Z) = S_0(t)^{\exp(\beta'Z)}\) 의 직접적 추정 형태이다.

식 8.8.4 의 직관: “기저 생존곡선을 공변량에 따라 위아래로 조절한다.”

\(\exp(b'Z_0) = 1\) (기저 수준): \(\widehat{S}(t|Z_0) = \widehat{S}_0(t)\). 기저 곡선 그대로.
\(\exp(b'Z_0) > 1\) (고위험): 기저 생존을 1 보다 큰 지수로 올리므로 \(\widehat{S}\) 가 더 빠르게 감소. 예: \(\widehat{S}_0 = 0.8\) 이고 \(\exp(b'Z_0) = 2\) 이면 \(\widehat{S} = 0.8^2 = 0.64\).
\(\exp(b'Z_0) < 1\) (저위험): \(\widehat{S}\) 가 기저보다 덜 감소. 예: \(\widehat{S}_0 = 0.8\) 이고 \(\exp(b'Z_0) = 0.5\) 이면 \(\widehat{S} = 0.8^{0.5} = 0.894\).

\(Z = 0\) 인 기저 개인과 \(Z = Z_0\) 인 개인의 생존곡선은 로그-로그 변환하면 평행하다:

\[ \log[-\log \widehat{S}(t|Z_0)] = \log[-\log \widehat{S}_0(t)] + b'Z_0 \]

이것이 바로 PH 가정의 시각적 표현이다.

3.4 분산 추정

정의: 생존 추정량의 분산

\[ \widehat{V}\left[\widehat{S}(t|Z_0)\right] = \left[\widehat{S}(t|Z_0)\right]^2 \left[Q_1(t) + Q_2(t; Z_0)\right] \tag{8.8.5} \]

여기서

\[ Q_1(t) = \sum_{t_i \leq t} \frac{d_i}{W(t_i; b)^2} \tag{8.8.6} \]

\[ Q_2(t; Z_0) = Q_3(t; Z_0)' \, \widehat{V}(b) \, Q_3(t; Z_0) \tag{8.8.7} \]

\(Q_3\) 의 \(k\)-번째 원소는

\[ Q_3(t; Z_0)_k = \sum_{t_i \leq t} \left[\frac{W^{(k)}(t_i; b)}{W(t_i; b)} - Z_{0k}\right] \frac{d_i}{W(t_i; b)}, \quad k = 1, \ldots, p \tag{8.8.8} \]

이며 \(W^{(k)}(t_i; b) = \sum_{j \in R(t_i)} Z_{jk} \exp(b'Z_j)\) 이다.

분산의 두 항은 서로 다른 불확실성 원천을 포착한다:

\(Q_1\) (baseline 불확실성): \(\beta\) 를 참값으로 알고 있더라도 기저 위험의 추정에서 발생하는 불확실성이다. Nelson-Aalen 분산의 일반화이다. \(\sum d_i / W^2\) 형태로, 사건 수가 적거나 위험집합이 작은 시점에서 크다.
\(Q_2\) (\(\beta\) 추정 불확실성): \(\beta\) 를 추정했기 때문에 추가로 발생하는 불확실성이다. \(Q_3\) 는 “위험집합의 평균 공변량”과 “\(Z_0\)” 의 차이에 비례한다. 예측하려는 \(Z_0\) 가 위험집합의 평균 공변량에서 멀수록 \(Q_2\) 가 커진다. 직관적으로, 데이터에서 드문 공변량 조합에 대한 예측은 불확실하다.

실무적 함의: 고위험 환자(극단적 \(Z_0\))의 생존 추정은 \(Q_2\) 가 크기 때문에 넓은 신뢰구간을 갖는다. 이는 “데이터가 부족한 영역에서의 외삽은 불확실하다”는 상식과 일치한다.

3.5 Klein Example 8.2 (continued): 후두암 60세 남성

Table 8.1 의 Cox 모형 (\(Z_1, Z_2, Z_3\): Stage II/III/IV dummy, \(Z_4\): age) 결과를 사용한다. 60세 남성에 대한 Stage 별 생존 추정:

Stage I: \(\widehat{S}(t) = \widehat{S}_0(t)^{\exp(0.0189 \times 60)}\)
Stage II: \(\widehat{S}(t) = \widehat{S}_0(t)^{\exp(0.0189 \times 60 + 0.1386)}\)
Stage III: \(\widehat{S}(t) = \widehat{S}_0(t)^{\exp(0.0189 \times 60 + 0.6383)}\)
Stage IV: \(\widehat{S}(t) = \widehat{S}_0(t)^{\exp(0.0189 \times 60 + 1.6931)}\)

5년 생존확률 추정:

Stage	\(\exp(b'Z_0)\)	\(\widehat{S}(5)\)	SE	95% CI (log-transform)
I	\(e^{1.134} = 3.11\)	0.703	0.074	(0.532, 0.822)
II	\(e^{1.273} = 3.57\)	0.667	0.106	(0.418, 0.829)
III	\(e^{1.772} = 5.88\)	0.513	0.095	(0.317, 0.679)
IV	\(e^{2.827} = 16.89\)	0.147	0.100	(0.022, 0.383)

해석:

Stage I 과 II 의 5년 생존확률은 비슷하다 (0.703 vs 0.667). 이는 § 8.5 에서 Stage I vs II 의 hazard ratio 가 비유의(RR=1.15, p=0.76)했던 결과와 일관된다.
Stage IV 의 5년 생존확률은 14.7%로 매우 낮다. 신뢰구간도 (0.022, 0.383)으로 넓어 불확실성이 크다.
Stage IV 의 넓은 신뢰구간은 두 가지 원인이다: (1) Stage IV 환자 수가 적어 \(Q_1\) 이 큼, (2) \(\exp(b'Z_0) = 16.89\) 로 극단적이어서 \(Q_2\) 도 큼.

로그 변환 신뢰구간을 사용하는 이유

생존확률은 \([0, 1]\) 범위이므로, 선형 신뢰구간 \(\widehat{S} \pm z_{\alpha/2} \cdot \text{SE}\) 는 범위를 벗어날 수 있다. 로그 변환을 적용하면

\[ \exp\left[\log \widehat{S}(t) \pm z_{\alpha/2} \cdot \frac{\text{SE}(\widehat{S})}{\widehat{S}(t)}\right] \]

이 항상 \((0, 1)\) 에 속한다. 보완 로그-로그(complementary log-log) 변환도 비슷한 효과를 준다. Ch.4 § 4.3 에서와 동일한 원리이다.

3.6 대안 추정량 비교

Breslow 추정량 외에 3 가지 대안이 있다:

4 가지 생존 추정량

추정량 1: Breslow (식 8.8.3)

\[ \widehat{S}_1(t|0) = \exp\left[-\widehat{H}_0(t)\right] = \exp\left[-\sum_{t_i \leq t} \frac{d_i}{W(t_i; b)}\right] \]

NA 의 일반화. 연속적 지수 변환.

추정량 2: Breslow product-limit

\[ \widehat{S}_2(t|0) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{W(t_i; b)}\right) \]

KM 의 일반화. 각 사건시간에서 곱극한.

추정량 3: Kalbfleisch-Prentice (KP)

\[ \widetilde{H}_0(t) = \sum_{t_i \leq t} \left[1 - \left(1 - \frac{\delta_i \exp(b'Z_i)}{W(t_i; b)}\right)^{\exp(-b'Z_i)}\right] \]

개별 관측치별 기여. 동점 시 수치해법 필요. SAS 의 기본 추정량.

추정량 4: 보정 product-limit

\[ \widehat{S}_4(t|Z_0) = \prod_{t_i \leq t} \left[1 - \frac{d_i \exp(b'Z_0)}{W(t_i; b)}\right] \]

기저가 아닌 \(Z_0\) 에 대해 직접 곱극한.

Andersen & Klein (1996)의 Monte Carlo 비교 결과:

추정량	편향(bias)	MSE	음수 가능	권장
\(\widehat{S}_1\) (Breslow exp)	중간	중간	아니오	권장
\(\widehat{S}_2\) (Breslow PL)	최소	최소	예	권장
\(\widehat{S}_3\) (KP)	큼	큼	아니오	비권장
\(\widehat{S}_4\) (보정 PL)	최소	최소	예	권장

\(\widehat{S}_2\) 와 \(\widehat{S}_4\) 가 가장 우수하지만 음수가 될 수 있다. 이는 위험집합이 작은 꼬리(right tail)에서 발생하며, 극단적 공변량 값에 대한 예측을 시도할 때 나타난다. 음수는 “이 공변량 조합에 대해서는 이 시점 이후 예측이 부적절하다”는 경고 신호이다.

\(\widehat{S}_3\) (KP)은 SAS 의 기본 추정량이지만 Monte Carlo 에서 가장 성능이 나쁘다. R 의 survival::survfit 은 Breslow 추정량을 기본으로 사용한다.

신뢰구간: 네 추정량 모두 동일한 점근 분산(식 8.8.5)을 갖는다. 로그 변환 신뢰구간이 가장 우수하고, arcsine-square-root 변환이 차선이다. 선형 신뢰구간은 권장하지 않는다 (Andersen & Klein, 1996).

3.7 코드 예시

3.7.1 Step 1: 순수 Python 구현 (Breslow 추정 원리)

import numpy as np

times = np.array([1, 2, 3, 4, 5, 6, 7, 8])
events = np.array([1, 1, 1, 0, 1, 0, 1, 0])
stage = np.array([1, 2, 1, 3, 2, 1, 4, 3])

Z = np.zeros((len(times), 3))
Z[stage == 2, 0] = 1
Z[stage == 3, 1] = 1
Z[stage == 4, 2] = 1

beta = np.array([0.1386, 0.6383, 1.6931])

exp_bz = np.exp(Z @ beta)

event_times = np.sort(np.unique(times[events == 1]))
H0 = 0.0
S0_values = {}

for t in event_times:
    d_i = np.sum((times == t) & (events == 1))
    at_risk = times >= t
    W_i = np.sum(exp_bz[at_risk])
    H0 += d_i / W_i
    S0_values[t] = np.exp(-H0)

print("Breslow baseline cumulative hazard and survival:")
for t, s in S0_values.items():
    print(f"  t={t}: H0={-np.log(s):.4f}, S0={s:.4f}")

age = 60
age_coef = 0.0189
for stage_label, z_vec in [("I", [0, 0, 0]),
                            ("II", [0.1386, 0, 0]),
                            ("III", [0, 0.6383, 0]),
                            ("IV", [0, 0, 1.6931])]:
    lp = sum(z_vec) + age_coef * age
    S_final = list(S0_values.values())[-1] ** np.exp(lp)
    print(f"  Stage {stage_label}: exp(b'Z0)={np.exp(lp):.2f}, "
          f"S(t_max)={S_final:.4f}")

3.7.2 Step 2: lifelines 구현 (실무 활용)

import pandas as pd
from lifelines import CoxPHFitter

df = pd.DataFrame({
    'T': [1, 2, 3, 4, 5, 6, 7, 8, 10, 12],
    'E': [1, 1, 1, 0, 1, 0, 1, 0, 1, 0],
    'stage_II':  [0, 1, 0, 0, 1, 0, 0, 0, 1, 0],
    'stage_III': [0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
    'stage_IV':  [0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
    'age': [55, 60, 65, 70, 58, 62, 72, 68, 63, 59],
})

cph = CoxPHFitter()
cph.fit(df, duration_col='T', event_col='E')
cph.print_summary()

cph.plot_partial_effects_on_outcome(
    covariates='stage_IV',
    values=[0, 1],
    cmap='coolwarm'
)

new_patient = pd.DataFrame({
    'stage_II': [0, 1, 0, 0],
    'stage_III': [0, 0, 1, 0],
    'stage_IV': [0, 0, 0, 1],
    'age': [60, 60, 60, 60],
})

surv_func = cph.predict_survival_function(new_patient)
print(surv_func)

3.7.3 R 코드 (survival 패키지)

library(survival)

fit <- coxph(Surv(time, status) ~ stage + age, data = larynx)
summary(fit)

new_data <- data.frame(
  stage = factor(c("I", "II", "III", "IV")),
  age = rep(60, 4)
)

surv_est <- survfit(fit, newdata = new_data)
summary(surv_est, times = 5)
plot(surv_est, col = 1:4, lwd = 2,
     xlab = "Years", ylab = "Survival Probability",
     main = "Estimated Survival: 60-year-old by Stage")
legend("bottomleft", levels(new_data$stage),
       col = 1:4, lwd = 2)

4 관련 주제

선행 지식

§ 8.1~8.2 — Cox 모형 · Coding
§ 8.3~8.4 — Partial Likelihood · Ties
§ 8.5~8.6 — Local Tests · Discretizing
Ch.4 — 비모수 추정 (KM · NA) — Breslow 추정량은 NA 의 일반화

후속 주제

§ 8.9 Exercises (예정) — Ch.8 연습문제 풀이
Ch.9 — Cox Refinements — 시간의존 공변량, 층화 모형
Ch.11 — Regression Diagnostics — Cox-Snell · 마팅게일 · Schoenfeld 잔차
Ch.12 — Parametric Regression — AFT 모형, Weibull 회귀

관련 개념

GLM Ch.12 — Models for Survival Data — GLM 관점의 생존 모형
AIC 와 모형 선택 — 정보이론 기반 모형 비교