Klein § 6.1~6.2 — Kernel Hazard Smoothing

1 들어가며 — Ch.6 첫 번째 deep-dive

편	주제
Ch.6 Overview	4 절 조망
§ 6.1~6.2 (본 편)	Kernel Hazard Smoothing 정전 풀이
§ 6.3 (예정)	Excess Mortality (multiplicative + additive)
§ 6.4 (예정)	Bayesian NPMLE (Dirichlet + beta process)
§ 6.5 (예정)	10 Exercises

§ 6.1~6.2 의 한 줄 요약

“Ch.4 의 NA 추정량 \(\widetilde{H}(t)\) 는 step function — 점프 크기 \(\Delta\widetilde{H}(t_i) = d_i/Y_i\) 가 hazard 의 거친 추정. 그러나 점프만으로 hazard 모양 파악 어려움. § 6.2 의 kernel smoothing 은 점프들의 시간 가중 평균 \(\widehat{h}(t) = (1/b) \sum K((t-t_i)/b) \Delta\widetilde{H}(t_i)\) 으로 매끄러운 hazard 추정. 두 핵심 선택 — bandwidth \(b\) (bias-variance trade-off) 와 kernel \(K()\) (Epanechnikov 가 MSE optimal). Boundary effect 는 Gasser-Müller 1979 의 asymmetric kernel 로 해결, optimal bandwidth 는 Ramlau-Hansen 1983 의 cross-validation 으로 선택.”

2 § 6.1 — Introduction

2.1 Ch.4 NA 의 한계

직관 — NA 가 hazard rate 를 직접 추정 못 하는 이유

Ch.4 § 4.2 의 NA 추정량:

\[ \widetilde{H}(t) = \sum_{t_i \leq t} \frac{d_i}{Y_i} \]

연속 시간 hazard 의 정의:

\[ h(t) = \frac{dH(t)}{dt} \]

→ NA 의 미분이 hazard. 그러나 NA 는 step function — 사건 시점에서만 점프, 사이는 상수.

• 점프 시점 직후 \(\widetilde{H}'(t_i^+) = \infty\).
• 점프 사이 \(\widetilde{H}'(t) = 0\).

→ “hazard 가 사건 시점에 무한대, 사이는 0” 이라는 무의미한 결론.

점프 크기 \(\Delta\widetilde{H}(t_i) = d_i/Y_i\): \(t_i\) 시점 hazard 의 거친 추정 (Ch.4 § 4.2 의 hazard rate \(h(t_i) \cdot dt\)).

함정 — 점프만 보고 hazard 모양 파악 불가

NA plot 의 점프 패턴이 hazard 모양과 직결되지 않음:

• Constant hazard (Exponential): 점프가 일정하지 않을 수 있음 (분모 \(Y_i\) 가 변하므로).
• Monotone increasing hazard (Weibull α>1): 점프가 시간에 따라 커지나, 작은 \(Y_i\) 효과로 혼란.
• U-shape hazard (bathtub): 두 봉우리 vs 중앙의 차이가 점프로 명확히 안 보임.

→ “점프들의 시간 평균” 이 진짜 hazard 추정. 시간 윈도우 평활 이 필요.

2.2 Ch.6 의 3 정교화

본 편은 § 6.2 의 kernel smoothing 만 다룬다. § 6.3 (excess mortality) + § 6.4 (Bayesian) 는 추후 deep-dive.

3 § 6.2 — Kernel-Smoothed Hazard Estimator

3.1 Kernel Smoothing 식

정의: Kernel-smoothed hazard estimator (식 6.2.4)

시점 \(t\) 에서 (\(b \leq t \leq t_D - b\) 영역):

\[ \widehat{h}(t) = \frac{1}{b} \sum_{i=1}^D K\left(\frac{t - t_i}{b}\right) \Delta\widetilde{H}(t_i) \]

분산 (식 6.2.5):

\[ \sigma^2[\widehat{h}(t)] = \frac{1}{b^2} \sum_{i=1}^D K\left(\frac{t - t_i}{b}\right)^2 \Delta\widehat{V}[\widetilde{H}(t_i)] \]

여기서:

• \(b > 0\): bandwidth — “얼마나 넓은 시간 윈도우를 평균할지”.
• \(K(\cdot)\): kernel function — \([-1, +1]\) 에서 정의, \(\int K = 1\), \(K \geq 0\).
• \(\Delta\widetilde{H}(t_i) = d_i/Y_i\) (NA 점프).
• \(\Delta\widehat{V}[\widetilde{H}(t_i)] = d_i/Y_i^2\) (Aalen 분산 점프).

직관 — 식의 분해

식 6.2.4 의 각 인자 의미:

• \(\Delta\widetilde{H}(t_i) = d_i/Y_i\): \(t_i\) 시점 hazard 의 거친 추정 (점프 크기).
• \(K((t - t_i)/b)\): \(t_i\) 가 \(t\) 에 가까울수록 큰 가중, \(|t - t_i| > b\) 면 0.
• \(1/b\): 정규화 — bandwidth 가 클수록 더 많은 점이 가중되므로 분모 보정.
• 합 \(\sum\): 모든 사건 시점의 가중 평균.

비유 — 비 오는 날의 빗소리:

• 빗방울 하나하나 = 사건 시점의 점프 (point event).
• “초당 빗방울 수” = hazard rate (count per unit time).
• 우리가 듣는 비의 강도 = kernel 가중 평균.
• Bandwidth = 귀가 평균하는 시간 윈도우 (1 초 vs 10 초).

→ NA 점프들을 “시간 윈도우 내 평균” 으로 연결한 것이 kernel-smoothed hazard.

3.2 3 종 Kernel — 식 6.2.1~6.2.3

정의: 3 가지 표준 kernel

Uniform (식 6.2.1):

\[ K(x) = \frac{1}{2}, \quad -1 \leq x \leq 1 \]

Epanechnikov (식 6.2.2):

\[ K(x) = \frac{3}{4}(1 - x^2), \quad -1 \leq x \leq 1 \]

Biweight (식 6.2.3):

\[ K(x) = \frac{15}{16}(1 - x^2)^2, \quad -1 \leq x \leq 1 \]

세 kernel 모두 \(\int_{-1}^1 K(x) dx = 1\) 만족.

직관 — 3 kernel 의 시각적 차이

각 kernel 의 모양 (Klein Figure 6.3 에서 시각화):

Kernel	\(K(0)\) 중심 가중	\(K(\pm 1)\) 양 끝	모양
Uniform	0.5	0.5	사각형 (cliff)
Epanechnikov	0.75	0	포물선 (smooth)
Biweight	0.9375	0	더 뾰족한 포물선

중심 대 양 끝 비:

• Uniform: 1 (모든 점 동등).
• Epanechnikov: \(\infty\) (양 끝 0).
• Biweight: \(\infty\), 중심에 더 강한 집중.

MSE 비교 (Hodges-Lehmann 1956 + Epanechnikov 1969):

• Epanechnikov 가 MSE optimal (asymptotic mean squared error 최소).
• Biweight 는 약간 더 큰 MSE 이지만 결과가 더 매끄러움.
• Uniform 은 가장 큰 MSE — 이론적으로 비추.

실무 권장:

• Epanechnikov 가 표준 (MSE 최적).
• Biweight 는 시각적으로 매끄러운 plot 필요 시.
• Uniform 은 단순 평균이 명확히 필요할 때만.

3.3 손풀이 — Klein Example 6.1 t=150 (Symmetric Kernel)

BMT ALL 38 명, \(t = 150\), Epanechnikov + \(b = 100\) days

\(t = 150\) 은 \(b \leq t \leq t_D - b\) 범위 (\(100 \leq 150 \leq 562\)) 안. Symmetric kernel 사용.

식 6.2.4 직접 계산. Klein Table 6.1 의 일부:

\(t_i\)	\(\Delta\widetilde{H}\)	\((150-t_i)/100\)	\(K(\cdot)\) Epanechnikov	기여 \(K \cdot \Delta\widetilde{H}\)
55	0.0270	0.95	\(0.75(1-0.95^2) = 0.0731\)	0.00197
74	0.0278	0.76	\(0.75(1-0.76^2) = 0.3168\)	0.00881
86	0.0286	0.64	\(0.75(1-0.64^2) = 0.4428\)	0.01266
104	0.0294	0.46	\(0.75(1-0.46^2) = 0.5913\)	0.01738
122	0.0667	0.28	\(0.75(1-0.28^2) = 0.6912\)	0.04610
129	0.0357	0.21	\(0.75(1-0.21^2) = 0.7169\)	0.02559
172	0.0370	-0.22	\(0.75(1-0.22^2) = 0.7137\)	0.02641
192	0.0385	-0.42	\(0.75(1-0.42^2) = 0.6177\)	0.02378
194	0.0400	-0.44	\(0.75(1-0.44^2) = 0.6048\)	0.02419
230	0.0435	-0.80	\(0.75(1-0.80^2) = 0.2700\)	0.01175

(다른 \(t_i\) 는 \(|t - t_i|/b > 1\) 이라 \(K = 0\) — 즉 \(t_i < 50\) 또는 \(t_i > 250\)).

합 \(\approx 0.2566\).

\[ \widehat{h}(150) = \frac{0.2566}{100} = 0.00257 \text{ /day} \]

연 환산: \(0.00257 \times 365 = 0.94\) → “연 94% hazard rate”.

임상 함의: BMT 후 150 일 시점 ALL 환자의 시간당 hazard 가 매우 높음. 절반 이상의 환자가 1 년 안에 사건 (재발 또는 사망) 경험.

표준오차 (식 6.2.5): \(\sigma(\widehat{h}(150)) = 0.00073\) → 95% CI 약 \((0.0011, 0.0040)\) /day.

3.4 Boundary Effect — Asymmetric Kernel (식 6.2.6~6.2.8)

함정 — 시작·끝 영역에서 symmetric kernel 작동 안 함

대칭 kernel 은 \(-1 \leq x \leq 1\) 모두 데이터 가능 가정.

시작 영역 (\(t < b\)):

• \(t - b < 0\) 영역에 사건 데이터 없음 (시간 음수 안 됨).
• Kernel 의 일부가 “데이터 없는 영역” 으로 빠짐.
• 적분이 1 보다 작아짐 → 추정량 음의 편향.

끝 영역 (\(t > t_D - b\)):

• \(t + b > t_D\) 영역에 사건 데이터 없음.
• 동일한 문제.

→ Boundary 에서 단순 truncation 은 부정확.

정의: Gasser-Müller (1979) Asymmetric Kernel

\(q = t/b\) (boundary 까지의 비율) 정의.

Uniform asymmetric (식 6.2.6):

\[ K_q(x) = \frac{4(1+q^3)}{(1+q)^4} + \frac{6(1-q)}{(1+q)^3} x, \quad -1 \leq x \leq q \]

Epanechnikov asymmetric (식 6.2.7):

\[ K_q(x) = K(x) (\alpha_E + \beta_E x), \quad -1 \leq x \leq q \]

with

\[ \alpha_E = \frac{64(2 - 4q + 6q^2 - 3q^3)}{(1+q)^4 (19 - 18q + 3q^2)}, \quad \beta_E = \frac{240(1-q)^2}{(1+q)^4 (19 - 18q + 3q^2)} \]

Biweight asymmetric (식 6.2.8): 유사한 \(\alpha_{BW} + \beta_{BW} x\) 형태.

오른쪽 tail 처리: \(q = (t_D - t)/b\) 대입 + \(K_q(x)\) 의 \(x\) 를 \(-x\) 로 치환.

직관 — boundary kernel 의 의미

Asymmetric kernel 의 효과:

• “데이터 없는 영역” 이 잘려나가는 만큼, 남은 영역의 가중을 키워 적분 = 1 회복.
• \(q = 1\) (\(t = b\)) 면 symmetric kernel 와 일치.
• \(q < 1\) 면 한쪽이 잘려 비대칭 모양.

Klein Example 6.1 의 \(t = 50\) (boundary):

\(q = 50/100 = 0.5\). 식 6.2.7 의 \(\alpha_E, \beta_E\):

\[ \alpha_E = \frac{64(2 - 2 + 1.5 - 0.375)}{(1.5)^4 (19 - 9 + 0.75)} = \frac{64 \cdot 1.125}{5.0625 \cdot 10.75} = \frac{72}{54.42} = 1.323 \]

\[ \beta_E = \frac{240 \cdot 0.25}{5.0625 \cdot 10.75} = \frac{60}{54.42} = 1.102 \]

→ \(K_{0.5}(-0.05) = 0.75 \cdot [1.323 + 1.102 \cdot (-0.05)] \cdot (1 - 0.05^2) = 0.75 \cdot 1.268 \cdot 0.9975 = 0.9485\).

식 6.2.4 적용: \(\widehat{h}(50) = 0.0015\) /day.

왜 작아지는가: \(t = 50\) 이전 (즉 \(t_i < 0\)) 의 데이터 없으므로 가중 분배가 한쪽으로 치우쳐 추정값이 작음.

3.5 Bandwidth Bias-Variance Trade-off

bandwidth 선택의 trade-off

작은 \(b\):

• 좁은 윈도우 → 적은 점들 평균 → 분산 큼.
• 시간 변동 민감 → bias 작음.
• Plot: 거칠고 noise 많음.

큰 \(b\):

• 넓은 윈도우 → 많은 점들 평균 → 분산 작음.
• 시간 변동 평탄화 → bias 큼.
• Plot: 매끄럽지만 진짜 모양 놓침.

→ MSE = bias² + variance 최소화하는 optimal bandwidth 가 존재.

직관 — Klein Figure 6.4 (kidney transplant)

같은 데이터에 4 가지 bandwidth (0.5, 1.0, 1.5, 2.0 yr) Epanechnikov:

• \(b = 0.5\): 매우 거침 (peaks 와 valleys 많음).
• \(b = 1.0\): 적당히 매끄러움.
• \(b = 1.5\): 매우 매끄러움 (1 ~ 3 년의 변동 거의 안 보임).
• \(b = 2.0\): 거의 평탄 (정보 많이 손실).

→ “어느 정도 매끄러운 게 좋은가” 는 진짜 hazard 의 모양에 의존. 데이터 기반 cross-validation 이 객관적 선택.

3.6 MISE 와 Bias 분해

정의: Mean Integrated Squared Error (MISE)

\[ \text{MISE}(b) = E\int_{\tau_L}^{\tau_U} [\widehat{h}(u) - h(u)]^2 du \]

\(\widehat{h}\) 의 추정 대상이 \(h^*(t) = (1/b) \int K((t-u)/b) h(u) du\) (smoothed version) 임을 인식하면 다음 분해:

\[ \text{MISE} = \underbrace{\int [h^*(u) - h(u)]^2 du}_{\text{bias}^2} + \underbrace{\int E[\widehat{h}(u) - h^*(u)]^2 du}_{\text{variance}} \]

직관 — Bias 의 점근 형태 (Theoretical Note 2)

큰 표본에서:

\[ \text{Bias}[\widehat{h}(t)] \approx \frac{1}{2} b^2 h''(t) k^* \]

where \(k^* = \int_{-1}^1 s^2 K(s) ds\) (kernel 의 second moment).

• \(b^2\): bandwidth 가 클수록 bias 큼 (제곱으로).
• \(h''(t)\): hazard 의 곡률 — 직선이면 bias 0, 곡선이면 bias 큼.
• \(k^*\): kernel 별 상수 (Epanechnikov \(k^* = 0.2\), biweight \(k^* = 1/7\)).

핵심: hazard 가 진짜 변동이 큰 시점 (변곡점) 에서 bias 가 크다. → kernel smoothing 은 곡률이 큰 구간을 평탄화.

3.7 Cross-Validation Optimal Bandwidth (Ramlau-Hansen 1983)

정의: Cross-Validation Bandwidth Selection

최적 \(b\) 는 다음 함수 \(g(b)\) 를 최소화:

\[ g(b) = \int_{\tau_L}^{\tau_U} \widehat{h}^2(u) du - \frac{2}{b} \sum_{i \neq j} K\left(\frac{t_i - t_j}{b}\right) \Delta\widetilde{H}(t_i) \Delta\widetilde{H}(t_j) \]

수치 계산:

• 첫 항: \(\widehat{h}\) 를 grid 점 (\(u_1 < \cdots < u_M\)) 에서 계산 후 trapezoid rule 적분:

\[ \sum_{i=1}^{M-1} \frac{u_{i+1} - u_i}{2} [\widehat{h}^2(u_i) + \widehat{h}^2(u_{i+1})] \]

• 두 번째 항: 모든 pair \((i, j)\), \(i \neq j\) 의 합산.

직관 — Cross-validation 의 의미

핵심 아이디어 (Ramlau-Hansen 1983a·b):

“MISE 의 둘째·셋째 항을 데이터로 추정 — 첫 항은 직접 계산, 둘째는 leave-one-out 형식으로 추정”.

• MISE 첫 항 \(\int h^{*2}\): \(\int \widehat{h}^2\) 으로 추정 (직접).
• MISE 둘째 항 \(-2 \int h^* h\): \(-2/b \sum_{i \neq j} K \Delta\widetilde{H}_i \Delta\widetilde{H}_j\) 으로 추정 (cross-validation).
• MISE 셋째 항 \(\int h^2\): \(b\) 와 무관 → 무시.

→ \(g(b)\) 의 최소값을 주는 \(b^*\) 가 데이터 기반 optimal bandwidth.

3.8 Klein Example 6.2 — Kidney Transplant

Cross-validation 결과 (Klein Figure 6.5)

OSU 863 명 kidney transplant 의 hazard 추정 (\(\tau_L = 0\), \(\tau_U = 6\) yr):

Kernel	Optimal \(b\)	\(g(b)\) at minimum
Uniform	0.17 yr	작음
Epanechnikov	0.20 yr	더 작음 (best)
Biweight	0.23 yr	약간 큼

관찰:

• Epanechnikov 가 가장 작은 \(g(b)\) — MSE 관점 최적.
• Biweight 의 \(b\) 가 더 큼 — biweight 가 자체 smoothness 가 크므로 더 넓은 윈도우 필요.
• Uniform 의 \(b\) 가 작음 — uniform 이 거칠어 작은 윈도우라도 OK.

Klein Figure 6.6: biweight + \(b = 0.23\) yr 의 hazard plot + 95% pointwise CI. 첫 1 년에 hazard 정점 (0.06~0.08), 그 이후 안정.

3.9 95% Pointwise Confidence Interval

정의: Log-transformed CI for \(\widehat{h}(t)\)

\[ \widehat{h}(t) \exp\left[\pm \frac{Z_{1-\alpha/2} \sigma(\widehat{h}(t))}{\widehat{h}(t)}\right] \]

(Ch.4 § 4.3 의 cumulative hazard CI 와 유사한 형태).

함정 — Practical Note 1: 추정 대상은 \(h^*\) 가 아닌 진짜 \(h\)

CI 가 추정하는 것은 smoothed hazard \(h^*(t)\):

\[ h^*(t) = \frac{1}{b} \int K\left(\frac{t-u}{b}\right) h(u) du \]

진짜 \(h(t)\) 가 아닌 그것의 평활 버전. CI 는 \(h^*\) 에 대한 것이므로 \(h\) 와의 bias 가 있음 (Theoretical Note 2 의 \(0.5 b^2 h''(t) k^*\)).

→ “곡률이 큰 시점” 에서 CI 해석 주의. 작은 \(b\) 일수록 \(h^* \approx h\) 이지만 분산 큼.

4 Practical Notes (3 개)

Practical Note 2 — Left-Truncation 도 적용 가능

식 6.2.4 의 모든 도구는 \(\widetilde{H}\) 와 그 분산이 있으면 적용 가능. 따라서:

• Right-censored 데이터 (Ch.4 § 4.2): 표준 NA 사용.
• Left-truncated 데이터 (Ch.4 § 4.6): \(Y_i\) 재정의 후 NA 그대로 적용.

→ Ch.5 의 다양한 sampling scheme 에서도 (Turnbull NPMLE 의 출력 이후) kernel smoothing 가능.

Practical Note 3 — Ramlau-Hansen 1983 source

• Smoothed estimator 의 첫 도입: Ramlau-Hansen (1983a, 1983b).
• Large-sample 점근 성질: Andersen et al. (1993), Ch.IV.2.
• Smoothing 일반론: Izenman (1991) survey.

R muhaz 패키지가 표준 구현. Python 에는 직접 구현 또는 scikit-survival 의 일부 기능.

5 Theoretical Notes (2 개)

Theoretical Note 1 — MISE 의 Bias-Variance 분해

\(\text{MISE}(b) = \text{bias}^2 + \text{variance}\).

• Bias: \(\int [h^*(u) - h(u)]^2 du\) — kernel 가 진짜 hazard 를 평탄화한 양.
• Variance: \(\int E[\widehat{h}(u) - h^*(u)]^2 du\) — 표본 추정의 변동.

작은 \(b\) → bias 작고 variance 큼. 큰 \(b\) → bias 크고 variance 작음.

→ Optimal \(b\) 가 두 항을 균형.

Theoretical Note 2 — Bias 의 점근 형태

큰 \(n\) 에서:

\[ \text{Bias}[\widehat{h}(t)] \approx \frac{1}{2} b^2 h''(t) k^* \]

where \(k^* = \int_{-1}^1 s^2 K(s) ds\).

각 kernel 의 \(k^*\) 값:

• Uniform: \(k^* = 1/3\).
• Epanechnikov: \(k^* = 1/5\).
• Biweight: \(k^* = 1/7\).

→ Biweight 가 가장 작은 bias coefficient — smoother kernel 일수록 같은 \(b\) 에서 bias 작음.

6 응용 분야

분야	Hazard 모양	Kernel + b
BMT post-transplant	Hump (peak 100~150 days)	Epanechnikov, \(b\) ≈ 100 days
Kidney transplant	Hump (peak 1 yr) + plateau	Biweight, \(b\) = 0.23 yr
만성 질환 (당뇨)	Constant or slowly increasing	Uniform 또는 Epanechnikov, 큰 \(b\)
인구 mortality	Bathtub (infant + senescence)	작은 \(b\) 로 변동 보존
임상시험 secondary 분석	Variable	Cross-validation 권장

7 코드 예시

7.1 Step 1 — R muhaz 패키지

library(muhaz)
library(KMsurv)
data(bmt)

# Klein Example 6.1 — BMT ALL hazard with Epanechnikov b=100
all_data <- subset(bmt, group == 1)

fit <- muhaz(times = all_data$t1,
             delta = all_data$d1,
             max.time = 730,        # 2 yr
             bw.method = "fixed",
             bw.smooth = 100,        # b = 100 days
             kern = "epanechnikov",
             b.cor = "left")          # boundary correction (left tail)

plot(fit, xlab = "Days", ylab = "h(t)",
     main = "BMT ALL — Epanechnikov b=100")

# t = 150 의 hazard 추정
h_150 <- approx(fit$est.grid, fit$haz.est, xout = 150)$y
print(paste("h(150) =", round(h_150, 5)))  # ≈ 0.00257

7.2 Step 2 — Cross-validation Optimal Bandwidth

# Klein Example 6.2 — kidney transplant CV
library(muhaz)
library(KMsurv)
data(kidtran)

fit_cv <- muhaz(times = kidtran$time,
                delta = kidtran$delta,
                max.time = 6 * 365,    # 6 yr in days
                bw.method = "global",   # MISE-based CV
                kern = "epanechnikov")

print(fit_cv$bw.glob)   # CV optimal b ≈ 0.20 yr (in study units)

plot(fit_cv, conf.level = 0.95,        # log-transformed CI
     main = "Kidney Transplant Hazard (Optimal b)")

7.3 Step 3 — Python 직접 구현 (Symmetric Kernel)

import numpy as np
from scipy.stats import norm

def epanechnikov(x):
    """Epanechnikov kernel (식 6.2.2)"""
    return np.where(np.abs(x) <= 1, 0.75 * (1 - x**2), 0.0)

def kernel_smoothed_hazard(times, events, t_eval, b, kernel=epanechnikov):
    """식 6.2.4 + 6.2.5 (symmetric kernel only)"""
    sorted_idx = np.argsort(times)
    t_sorted = times[sorted_idx]
    e_sorted = events[sorted_idx]
    event_times = np.unique(t_sorted[e_sorted == 1])

    # NA 점프 + Aalen variance 점프
    delta_H = np.zeros_like(event_times, dtype=float)
    delta_V = np.zeros_like(event_times, dtype=float)

    for i, t_i in enumerate(event_times):
        Y = np.sum(t_sorted >= t_i)
        d = np.sum((t_sorted == t_i) & (e_sorted == 1))
        delta_H[i] = d / Y
        delta_V[i] = d / Y ** 2

    # 식 6.2.4 + 6.2.5
    h_hat = np.zeros_like(t_eval, dtype=float)
    sigma_h = np.zeros_like(t_eval, dtype=float)

    for k, t in enumerate(t_eval):
        K_vals = kernel((t - event_times) / b)
        h_hat[k] = np.sum(K_vals * delta_H) / b
        sigma_h[k] = np.sqrt(np.sum(K_vals**2 * delta_V) / b**2)

    return h_hat, sigma_h


# Klein Example 6.1
times = np.array([1, 55, 74, 86, 104, 107, 109, 110, 122, 122,
                  129, 172, 192, 194, 230, 276, 332, 383, 418,
                  466, 487, 526, 609, 662, 2081])
events = np.concatenate([np.ones(23), [0, 0]])  # 23 events + 2 censored

t_eval = np.array([150])
h, se = kernel_smoothed_hazard(times, events, t_eval, b=100)
print(f"h(150) = {h[0]:.5f} (SE {se[0]:.5f})")  # ≈ 0.00257 (0.00073)

7.4 Step 4 — Bandwidth Cross-Validation

def g_function(b, times, events, kernel=epanechnikov,
                tau_L=0, tau_U=730, n_grid=100):
    """식 (Klein § 6.2 Cross-validation g(b))"""
    # 첫 항 — trapezoid rule of ∫ ĥ²
    grid = np.linspace(tau_L, tau_U, n_grid)
    h_grid, _ = kernel_smoothed_hazard(times, events, grid, b, kernel)
    integral_1 = np.trapz(h_grid**2, grid)

    # 두 번째 항 — cross-validation
    sorted_idx = np.argsort(times)
    e_sorted = events[sorted_idx]
    t_sorted = times[sorted_idx]
    event_times = np.unique(t_sorted[e_sorted == 1])
    delta_H = []
    for t_i in event_times:
        Y = np.sum(t_sorted >= t_i)
        d = np.sum((t_sorted == t_i) & (e_sorted == 1))
        delta_H.append(d / Y)
    delta_H = np.array(delta_H)

    sum_2 = 0
    for i in range(len(event_times)):
        for j in range(len(event_times)):
            if i != j:
                sum_2 += kernel((event_times[i] - event_times[j]) / b) \
                         * delta_H[i] * delta_H[j]

    return integral_1 - 2 * sum_2 / b


# Optimal b 탐색
from scipy.optimize import minimize_scalar
result = minimize_scalar(g_function, bounds=(10, 300),
                          method='bounded',
                          args=(times, events))
print(f"Optimal b = {result.x:.1f} days")

8 핵심 takeaway

§ 6.1~6.2 의 5 가지 교훈

1. NA 의 점프 = 거친 hazard 추정 — \(\Delta\widetilde{H}(t_i) = d_i/Y_i\) 가 \(h(t_i) \cdot dt\) 의 추정. 그러나 점프만으로 hazard 모양 파악 어려워 평활 필요.

2. Kernel-smoothed estimator (식 6.2.4) — \(\widehat{h}(t) = (1/b) \sum K \Delta\widetilde{H}\) — NA 점프들의 시간 가중 평균. Bandwidth \(b\) + kernel \(K()\) 가 두 핵심 선택.

3. Epanechnikov kernel 이 MSE optimal — Hodges-Lehmann 1956 + Epanechnikov 1969. Biweight 는 smoother, uniform 은 거침. 실무에서 Epanechnikov 표준.

4. Boundary effect → Gasser-Müller asymmetric kernel — 시작·끝 영역에서 symmetric kernel 의 적분 ≠ 1. 식 6.2.6~6.2.8 의 \(q = t/b\) 의존 modified kernel 로 보정.

5. Cross-validation (Ramlau-Hansen 1983) 으로 optimal \(b\) — MISE 의 두 항 (bias² + variance) 의 trade-off. \(g(b) = \int \widehat{h}^2 - 2/b \sum K \Delta\widetilde{H} \Delta\widetilde{H}\) 최소화. Bias 의 점근 \(\approx 0.5 b^2 h''(t) k^*\).

9 관련 주제

선행 지식

• Ch.6 Overview
• § 4.1~4.2 — KM·NA + Aalen 분산
• § 4.3~4.4 — log/arcsine 변환 + Brownian bridge

후속 주제

• § 6.3 — Excess Mortality (multiplicative + additive)
• § 6.4 — Bayesian NPMLE (Dirichlet + beta process)
• § 6.5 — 10 Exercises

관련 개념

• Density estimation 의 일반론 (Silverman 1986)
• Variance stabilizing 변환 (§ 4.3 의 arcsine)
• Cross-validation (machine learning 일반)
• Hodges-Lehmann (1956) optimal kernel
• Epanechnikov (1969) kernel 도입