Klein § 6.3~6.4 — Excess Mortality · Bayesian NPMLE

§ 6.3 Excess Mortality 두 모형 — Multiplicative h_j(t) = β(t)θ_j(t) (식 6.3.1) + B̂(t) = ∑ d_i/Q(t_i) (식 6.3.2) + Q(t) = ∑ θ_j(t)Y_j(t) — SMR (Breslow 1975) 의 시간 변동 일반화 / Additive h_j(t) = α(t) + θ_j(t) (식 6.3.4) + Â(t) = H̃(t) - Θ(t) (식 6.3.6) + Θ(t) = ∑ ∫ θ_j(u) Y_j(u)/Y(u) du (식 6.3.5) + corrected survival S^C(t) = Ŝ(t)/S*(t) / Klein Example 6.3 Iowa psychiatric 26 명 (Klein § 1.15) + 1959 Iowa mortality table — β ≈ 25 첫 2 년 (정신질환자 일반인 25 배 사망률) / § 6.4 Bayesian NPMLE — squared-error loss → posterior mean / Dirichlet process prior (Ferguson 1973) for S(t) — α(t,∞) = c S_0(t), 식 6.4.1 closed form S̃_D / Beta process prior (Hjort 1990) for H(t) — c(t) weight, 식 6.4.2 closed form S̃_B / Gibbs sampler (Gelfand-Smith 1990) — 식 6.4.3 Dirichlet density + 식 6.4.4 posterior mean — interval cens, double cens 등 일반 sampling scheme 처리 / Klein Example 6.4 6-MP 21 명, S_0(t) = e^{-0.1t}, c = 5 — 두 prior + Gibbs 결과 비교

Klein & Moeschberger Ch.6 의 § 6.3 (Excess Mortality) + § 6.4 (Bayesian NPMLE) 를 deep-dive. § 6.1~6.2 가 hazard 의 평활을 다뤘다면, 본 편은 (1) 표준 인구와 비교 (excess mortality) 와 (2) prior 정보 결합 (Bayesian) 의 두 정교화. § 6.3 — Excess Mortality: 환자 그룹의 mortality 를 표준 인구와 비교. 두 모형이 같은 데이터에 다른 질문에 답. Multiplicative (relative mortality) — h_j(t) = β(t) · θ_j(t) (식 6.3.1) 에서 β(t) 가 환자/표준 ratio. β(t) > 1 면 빠른 사망, β(t) = 1 표준 동일, β(t) < 1 더 좋음. 추정량 B̂(t) = ∑_{t_i ≤ t} d_i/Q(t_i) (식 6.3.2) where Q(t) = ∑_j θ_j(t) Y_j(t) — 모든 위험 환자의 표준 인구 hazard 가중 합. 분산 식 6.3.3. SMR (Standardized Mortality Ratio, Breslow 1975) 의 시간-변동 일반화 — β(t) = β_0 일정 가정 시 SMR = ∑ d_i / E where E = ∫ Q(u) du. Additive — h_j(t) = α(t) + θ_j(t) (식 6.3.4) 에서 α(t) 가 표준 외 추가 위험. 추정량 Â(t) = H̃(t) - Θ(t) (식 6.3.6) where Θ(t) = ∑_j ∫ θ_j(u) Y_j(u)/Y(u) du (식 6.3.5) — 표준 인구 가정 시 기대 cumulative hazard. 분산 식 6.3.7 = NA 분산. Corrected survival S^C(t) = Ŝ(t)/S(t) where S(t) = exp(-Θ(t)) — cancer relative survival 의 표준. Klein Example 6.3 26 명 Iowa psychiatric (Klein § 1.15) — 1959-1961 Iowa State life table 기반. λ_S(a) = -ln S(a) - (-ln S(a+1)) 로 표준 hazard 계산. Q(t), Θ(t), B̂(t), Â(t) 모두 손풀이. 첫 2 년 β(t) ≈ 20-30 → 정신질환자 일반인 대비 20-30 배 사망률. 30 년 후 cumulative excess Â(30) = 0.36 (95% CI 0.04~0.68) → 100 명당 36 명 추가 사망. § 6.4 — Bayesian NPMLE: prior \(\pi(S)\) + likelihood → posterior. Squared-error loss → Bayes estimator = posterior mean. Dirichlet process prior (Ferguson 1973) for S(t): α(t,∞) = c S_0(t) where S_0 = prior guess + c = prior strength. Prior mean = S_0(t), variance = S_0(1-S_0)/(c+1) — c+1 명의 가짜 표본. Posterior 도 Dirichlet (conjugate). Bayes estimator 식 6.4.1 의 closed form S̃_D(t). Beta process prior (Hjort 1990) for H(t): H_0(t) prior guess + c(t) 시간별 weight. Sample path 가 Dirichlet 보다 매끄러움. Closed form Bayes estimator 식 6.4.2 S̃_B(t). Gibbs sampler (Gelfand-Smith 1990): Dirichlet 와 beta process prior 외 다른 sampling scheme (interval cens, double cens) 또는 regression 의 일반화. 핵심 idea: censored 관측의 사건 시점을 잠재변수로 처리, multinomial sampling 으로 분배 → Dirichlet posterior 갱신 → 반복. 식 6.4.3 의 Dirichlet density + 식 6.4.4 의 posterior mean. Klein Example 6.4 6-MP 21 명 (Klein § 1.2), S_0(t) = e^{-0.1t}, c = 5. Dirichlet S̃_D(t) + beta S̃_B(t) closed form 손풀이. Gibbs 1000~10000 회 반복으로 동일 결과 도출. 두 prior 의 차이 — beta process 가 prior 에 더 가까움 (sample path 매끄러움). c → 0 또는 n → ∞ 에서 KM 으로 수렴 — Bayes 가 KM 의 일반화. Finite sample 에서 prior 가 안정화 효과.

Statistics
Survival Analysis
Klein-Moeschberger
Excess-Mortality
SMR
Bayesian-NPMLE
Dirichlet-Process
Gibbs-Sampler
저자

Kwangmin Kim

공개

2026년 04월 28일

1 들어가며 — Ch.6 의 두 번째·세 번째 정교화

주제
Ch.6 Overview 4 절 조망
§ 6.1~6.2 Kernel Hazard Smoothing
§ 6.3~6.4 (본 편) Excess Mortality + Bayesian NPMLE
§ 6.5 (예정) 10 Exercises
§ 6.3~6.4 의 한 줄 요약

“§ 6.3 의 excess mortality 는 환자 그룹의 mortality 를 표준 인구와 비교 — multiplicative \(\beta(t) = h_j/\theta_j\) 는 ‘몇 배 빠른가’, additive \(\alpha(t) = h_j - \theta_j\) 는 ‘얼마 추가인가’. SMR 의 시간-변동 일반화 + cancer relative survival 의 표준. § 6.4 의 Bayesian NPMLE 는 prior 정보 + KM 데이터를 결합 — Dirichlet process (Ferguson 1973) for \(S\) 와 beta process (Hjort 1990) for \(H\) 두 conjugate prior 의 closed form Bayes estimator + Gibbs sampler 의 일반화. 두 정교화 모두 Ch.4 의 raw KM/NA 위에 추가 layer 를 올리는 도구.”

2 § 6.3 — Excess Mortality

2.1 동기 — 왜 표준 인구와 비교하는가

직관 — 절대 mortality vs 상대 mortality

특정 환자 그룹의 mortality 를 보고할 때:

절대 척도 (Ch.4 KM):

“정신질환자 5 년 생존률 80%”

  • 절대값.
  • 같은 나이·성별의 일반인 대비 어떤지 모름.
  • 의미: “X% 환자 살아남음” — 임상 보고에 직관적이지만 비교 불가.

상대 척도 (§ 6.3 excess mortality):

“정신질환자가 같은 나이·성별 일반인보다 25 배 빠른 사망”

  • 상대값.
  • 표준 인구 대비 비교.
  • 의미: “환자 고유의 위험” — 표준 mortality 효과 제거.

→ 임상 보고에 두 척도 모두 필요. § 6.3 은 두 가지 모형으로 상대 척도 추정.

가정: 표준 인구 hazard \(\theta_j(t)\) 가 환자 \(j\) 의 특성 (나이·성별·인종 등) 에 따라 알려져 있음 — life table 또는 cohort 데이터.

2.2 Multiplicative Model — Relative Mortality

정의: Multiplicative (Relative Mortality) Model (식 6.3.1)

환자 \(j\) 의 hazard rate:

\[ h_j(t) = \beta(t) \cdot \theta_j(t), \quad j = 1, \ldots, n \]

  • \(\beta(t)\): relative mortality function (환자 그룹 / 표준 인구).
  • \(\beta(t) > 1\): 환자 그룹이 표준보다 빠른 사망.
  • \(\beta(t) = 1\): 표준과 동일.
  • \(\beta(t) < 1\): 환자 그룹이 더 좋음 (예: 운동 선수의 심혈관).

Cumulative: \(B(t) = \int_0^t \beta(u) du\).

정의: B̂(t) — Cumulative Relative Mortality (식 6.3.2)

\[ \widehat{B}(t) = \sum_{t_i \leq t} \frac{d_i}{Q(t_i)} \]

여기서 \(Q(t) = \sum_{j=1}^n \theta_j(t) Y_j(t)\).

분산 (식 6.3.3):

\[ \widehat{V}[\widehat{B}(t)] = \sum_{t_i \leq t} \frac{d_i}{Q(t_i)^2} \]

직관 — 식 6.3.2 의 분해
  • 분자 \(d_i\): 시점 \(t_i\) 에서 실제 관측된 사건 수.
  • 분모 \(Q(t_i) = \sum_j \theta_j(t_i) Y_j(t_i)\): “표준 인구 mortality 가정 시 \(t_i\) 시점의 기대 사건 수의 가중치”.
  • 비율 \(d_i / Q(t_i)\): “실제 / 기대” — 시간별 relative mortality 의 거친 추정.

NA (식 4.2.3) 와의 비교:

  • NA: \(\widetilde{H}(t) = \sum d_i / Y_i\) — 분모에 위험집합만.
  • B̂(t): \(\widehat{B}(t) = \sum d_i / Q(t_i)\) — 분모에 표준 인구 hazard 가중을 추가.

→ B̂ 가 NA 의 일반화.

직관 — SMR 의 시간-변동 일반화

Standardized Mortality Ratio (SMR) — Breslow (1975):

\(\beta(t) = \beta_0\) (일정 가정) 의 단일 추정량:

\[ \widehat{\beta}_0 = \frac{\sum_{i=1}^D d_i}{\int_0^{\infty} Q(u) du} = \frac{\text{관측 사건 수}}{\text{기대 사건 수}} \]

  • 분자 = 실제 사건 수 합.
  • 분모 = 표준 인구 가정 시 기대 사건 수 (모든 시점 적분).

SMR 의 한계: 시간 변동 무시. 환자 그룹의 mortality 가 시간에 따라 변할 때 (예: 첫 2 년 매우 높고 그 이후 안정) 단일 SMR 로는 시간 패턴 못 잡음.

B̂(t) 의 일반화: 시간 변동을 누적 곡선으로 표현. 기울기가 시간별 \(\beta(t)\) 의 거친 추정. → SMR 의 자연스러운 일반화.

2.3 Additive Model

정의: Additive Model (식 6.3.4)

\[ h_j(t) = \alpha(t) + \theta_j(t) \]

  • \(\alpha(t)\): excess (additional) hazard — 환자 그룹의 추가 위험.
  • \(\alpha(t) > 0\): 표준 외 추가 위험.
  • \(\alpha(t) = 0\): 표준과 동일.
  • \(\alpha(t) < 0\): 표준보다 안전 (드뭄).

Cumulative: \(A(t) = \int_0^t \alpha(u) du\).

정의: Â(t) — Cumulative Excess Mortality (식 6.3.6)

\[ \widehat{A}(t) = \widetilde{H}(t) - \Theta(t) \]

  • 첫 항: 관측 cumulative hazard (NA).
  • 두 번째 항 \(\Theta(t)\): “표준 인구 가정 시 기대 cumulative hazard” (식 6.3.5):

\[ \Theta(t) = \sum_{j=1}^n \int_0^t \theta_j(u) \frac{Y_j(u)}{Y(u)} du \]

where \(Y(t) = \sum_j Y_j(t)\).

분산 (식 6.3.7):

\[ \widehat{V}[\widehat{A}(t)] = \sum_{t_i \leq t} \frac{d_i}{Y(t_i)^2} \]

(NA 분산과 동일).

직관 — 식 6.3.5 의 의미

\(\Theta(t)\) 는 “표준 인구 hazard 의 시간별 가중 평균 적분”:

  • \(\theta_j(u)\): 환자 \(j\) 의 시점 \(u\) 표준 hazard.
  • \(Y_j(u)/Y(u)\): 시점 \(u\) 에 위험인 환자 중 \(j\) 의 비중.
  • \(\theta_j(u) \cdot Y_j(u)/Y(u)\): \(j\) 의 표준 hazard 기여.
  • \(\sum_j\): 모든 위험 환자의 평균 표준 hazard.
  • 적분 \(\int_0^t\): 시점 \(t\) 까지 누적.

→ “표준 인구 mortality 만 작용하는 가상 세계의 cumulative hazard”.

차이 \(\widehat{A}(t) = \widetilde{H}(t) - \Theta(t)\) = “관측 - 가상 세계” = excess.

→ 100 명당 추가 사망 수의 누적.

2.4 Corrected Survival Curve

정의: Corrected Survival (Cancer Relative Survival)

\[ S^C(t) = \frac{\widehat{S}(t)}{S^*(t)} \]

where:

  • \(\widehat{S}(t)\): 관측 KM.
  • \(S^*(t) = \exp(-\Theta(t))\): 표준 인구 기대 생존곡선.
  • \(S^C(t)\): “환자 그룹의 표준 인구 보정 survival”.
직관 — Cancer relative survival

암 역학의 표준 도구:

  • 관측 KM: “암 환자의 절대 5 년 생존률 60%”.
  • 일반 인구 기대 생존: “같은 나이·성별 일반인의 5 년 생존률 90%”.
  • Relative survival: \(60\% / 90\% = 67\%\) → “암 자체로 인한 사망 외 배제한 환자의 생존률”.

해석:

  • \(S^C(t) = 1\): 환자가 표준 인구와 동일.
  • \(S^C(t) < 1\): 환자가 표준보다 못함.
  • \(S^C(t) > 1\): 환자가 표준보다 좋음 (드뭄).

주의: \(S^C(t)\) 는 monotone non-increasing 가정 안 됨. 작은 \(t\) 에서 \(S^C > 1\) 일 수 있어 해석 주의.

2.5 Multiplicative vs Additive — 어느 모형?

두 모형이 답하는 다른 질문
모형 답하는 질문 비유
Multiplicative “위험이 몇 배?” 1.5 배 빠른 자동차 (모든 도로)
Additive “추가 위험이 얼마?” 시속 50 km 더 빠른 자동차 (모든 도로)

:

  • \(\theta_j(t) = 0.01\) 환자 (저위험) + \(\beta = 2\)\(h_j = 0.02\) (배 차이 적음).

  • \(\theta_j(t) = 0.05\) 환자 (고위험) + \(\beta = 2\)\(h_j = 0.10\) (배 차이 큼).

  • 같은 \(\beta\) 라도 baseline 에 따라 절대 hazard 차이 변동.

  • \(\theta_j(t) = 0.01\) + \(\alpha = 0.01\)\(h_j = 0.02\) (작은 baseline 에 큰 영향).

  • \(\theta_j(t) = 0.05\) + \(\alpha = 0.01\)\(h_j = 0.06\) (큰 baseline 에 작은 영향).

선택 기준:

  • 위험 요인이 비례적 효과 (예: 흡연 → 폐암, 모든 연령 공통 ratio): multiplicative.
  • 위험 요인이 추가적 효과 (예: 직업 노출 → 일정 추가 위험): additive.

→ 두 모형 모두 보고하고 비교 권장. Klein Example 6.3 도 두 모형 모두 추정.

2.6 Klein Example 6.3 — 26 명 Iowa Psychiatric

데이터 (Klein § 1.15)

Woolson (1981) 의 26 명 정신질환자 (Iowa). Entry age + study time + death indicator.

표준 인구: 1959-1961 Iowa State life table (US Dept. of Health and Human Services 1959).

  • \(\lambda_S(a)\): 나이 \(a\) 의 표준 hazard rate.
  • 식: \(\lambda_S(a) = -\ln S(a) - (-\ln S(a+1))\) — 1 년 구간 내 constant hazard 가정.

Klein Table 6.2 의 일부:

나이 \(\lambda_M\) (남) \(\lambda_F\) (여)
35 0.00181 0.00191
50 0.00810 0.00414
65 0.03020 0.01421
75 0.06795 0.04207
손풀이 — Q(t) 와 Θ(t) 계산

각 환자 \(j\) 의 entry age \(a_j\) 와 sex 에서:

\[ \theta_j(t) = \lambda_S(a_j + t) \]

(시점 \(t\) 에 환자 나이는 \(a_j + t\)).

Q(t) (식 6.3.2 분모):

\[ Q(t) = \sum_j \theta_j(t) Y_j(t) = \sum_{j: \text{at risk at } t} \lambda_S(a_j + t) \]

Θ(t) (식 6.3.5 의 적분):

연 단위 가정: 1 년 구간 내 constant hazard. 각 정수 시점 \(t\) 에서:

\[ \Theta(t) = \Theta(t-1) + \frac{\sum_{j: \text{at risk}} \lambda_S(a_j + t - 1)}{Y(t)} \]

(이전 누적 + 현재 구간 평균 표준 hazard).

비-정수 시점은 선형 보간.

결과 (Klein Table 6.4)

26 명 데이터의 첫 10 년 일부:

\(t_i\) \(d_i\) \(Y(t_i)\) \(\widetilde{H}\) \(\Theta\) \(\widehat{A}\) \(\text{SE}\) \(\widehat{S}\) \(S^*\) \(S^C\)
1 2 26 0.0769 0.0021 0.0748 0.0544 0.9231 0.9979 0.9250
2 1 24 0.1186 0.0041 0.1145 0.0685 0.8846 0.9959 0.8882
10 0 23 0.1186 0.0234 0.0952 0.0685 0.8846 0.9769 0.9056

B̂(t) (multiplicative) Klein Figure 6.8:

  • B̂(2) ≈ 25 → 첫 2 년 정신질환자가 일반 인구보다 25 배 빠른 사망.
  • B̂(10) ≈ 50 → 10 년 누적.
  • B̂(30) ≈ 250 → 30 년 누적.

Â(t) (additive) Klein § 6.3 본문:

  • 첫 2 년 기울기 약 0.05 → α(t) ≈ 0.05.
  • 2 < t < 21 기울기 약 0 → α(t) ≈ 0 (안정).
  • t > 21 기울기 약 0.05 → α(t) ≈ 0.05 (재상승).
  • Â(30) ≈ 0.36 (95% CI 0.04~0.68) → 100 명당 36 명 추가 사망 (30 년 누적).

임상 함의: 정신질환자의 mortality 는 일반 인구보다 매우 높음. 첫 2 년 + 후반부에 위험 집중. 사회적 지원 + 의료 모니터링 강화의 통계적 근거.

3 § 6.4 — Bayesian Nonparametric Methods

3.1 동기 — Prior 정보의 활용

직관 — 왜 Bayesian 비모수가 필요한가

표준 KM (Ch.4) 의 한계:

  • Prior 정보 무시 — 이전 연구·전문가 의견·표준 인구 mortality 등 외부 지식 활용 못 함.
  • Finite sample 불안정 — 작은 표본에서 KM 이 불안정 (특히 tail).

Bayesian 접근:

  • Prior \(\pi(S)\): \(S(t)\) 에 대한 사전 분포.
  • Likelihood \(L(\text{data} \mid S)\): 표본 데이터의 가능도.
  • Posterior \(\pi(S \mid \text{data}) \propto L \cdot \pi\).

Squared-error loss:

\[ L(S, \widehat{S}) = \int_0^\infty [\widehat{S}(t) - S(t)]^2 dw(t) \]

→ Bayes estimator = posterior mean (loss minimizer).

→ 큰 표본 (\(n \to \infty\)) 에서 prior 영향 0 → KM 으로 수렴. 작은 표본에서 prior 가 dominant.

두 conjugate prior (계산이 closed form 이어서 표준):

  1. Dirichlet process for \(S(t)\) — Ferguson (1973).
  2. Beta process for \(H(t)\) — Hjort (1990).

3.2 Dirichlet Process Prior — Ferguson (1973)

정의: Dirichlet Process Prior

\(S(t) \sim DP(\alpha)\) where \(\alpha\) 는 measure on \([0, \infty)\).

표준 형태: \(\alpha([t, \infty)) = c \cdot S_0(t)\).

  • \(S_0(t)\): prior guess (investigator 의 best estimate).
  • \(c\): prior strength.

Prior 통계:

\[ E[S(t)] = S_0(t), \quad V[S(t)] = \frac{S_0(t)[1 - S_0(t)]}{c + 1} \]

→ “\(c + 1\) 명의 가짜 표본” 에 해당하는 prior — \(c\) 클수록 강한 prior (분산 작음).

직관 — Dirichlet 의 sample paths

Klein Figure 6.11: \(S_0(t) = e^{-0.1 t}\), \(c = 5\) 의 10 개 sample paths.

  • 각 path 가 nonincreasing function (생존곡선의 기본 성질).
  • \(S(0) = 1\) — 시작 시점 정확히 1.
  • 연속 함수 — 점프 없음.
  • 그러나 약간 거침 — \(c\) 작을수록 거침.

\(c\) 의 의미:

  • \(c \to 0\): 매우 약한 prior, 거의 무관계.
  • \(c = 5\): “5 + 1 = 6 명의 가짜 표본” 가치.
  • \(c \to \infty\): 매우 강한 prior, \(S_0\) 에 고정.

→ 실무에서 \(c\) 는 “내가 prior 에 얼마나 확신하는가” 를 가짜 표본 크기로 표현.

정의: Bayes Estimator (식 6.4.1)

right-censored 데이터 + Dirichlet prior 의 closed form posterior mean:

\[ \widetilde{S}_D(t) = \frac{\alpha(t, \infty) + Y_{i+1}}{\alpha(0, \infty) + n} \prod_{k=1}^i \frac{\alpha(t_k, \infty) + Y_{k+1} + \lambda_k}{\alpha(t_k, \infty) + Y_{k+1}} \]

for \(t_i \leq t < t_{i+1}\), where:

  • \(\alpha(t, \infty) = c \cdot S_0(t)\).
  • \(Y_k\): 시점 \(t_k\) 직전 위험집합.
  • \(\lambda_k\): 시점 \(t_k\) 의 right-censored 수.

→ Closed form! Monte Carlo 불필요.

3.3 Beta Process Prior — Hjort (1990)

정의: Beta Process Prior

대안: \(H(t)\) 에 prior 부여.

작은 구간 \([a_{i-1}, a_i)\) 에서 \(W_i = H(a_i) - H(a_{i-1})\) 가 베타 분포:

\[ W_i \sim \text{Beta}(p_i, q_i) \]

with:

  • \(p_i = c \cdot [H_0(a_i) - H_0(a_{i-1})]\).
  • \(q_i = c \cdot \{1 - [H_0(a_i) - H_0(a_{i-1})]\}\).

Limit 구간 length → 0 → beta process.

정의: Bayes Estimator (식 6.4.2)

\(c\) 가 일정한 경우:

\[ \widetilde{S}_B(t) = \exp\left\{-\sum_{k=1}^i \frac{c[H_0(t_k) - H_0(t_{k-1})]}{c + Y_k} - \frac{c[H_0(t) - H_0(t_i)]}{c + Y_{i+1}}\right\} \prod_{k: t_k \leq t}\left[1 - \frac{c h_0(t_k) + d_k}{c + Y_k}\right]^{\Delta_k} \]

직관 — Dirichlet vs Beta process
측면 Dirichlet Beta Process
Prior on \(S(t)\) \(H(t)\)
도메인 \([0, 1]\) (확률) \([0, \infty)\) (cumulative hazard)
Sample path 약간 거침 더 매끄러움
Closed form 식 6.4.1 식 6.4.2
권장 \(S_0\) 형태 prior \(h_0\) 형태 prior

왜 beta process 가 더 매끄러운가: \(H\) 의 sample path 는 monotone non-decreasing 자연. \(S = e^{-H}\) 변환은 매끄러움 보존.

c → 0 또는 n → ∞: 두 prior 모두 KM 으로 수렴.

3.4 Klein Example 6.4 — 6-MP 손풀이

6-MP 21 명 (Klein § 1.2) — 두 prior + 손풀이

Prior 설정:

  • \(S_0(t) = e^{-0.1 t}\)\(\alpha(t, \infty) = 5 \cdot e^{-0.1t}\).
  • \(c = 5\) (“5 + 1 = 6 명 가짜 표본”).
  • Beta process: \(H_0(t) = 0.1 t\), \(c(t) = 5\).

\(t \in [0, 6)\) 시점 (사건 시점 이전):

Dirichlet (식 6.4.1):

\[ \widetilde{S}_D(t) = \frac{5 e^{-0.1t} + 21}{5 + 21} = \frac{5 e^{-0.1t} + 21}{26} \]

  • \(t = 0\): \(\widetilde{S}_D(0) = (5 + 21)/26 = 1.0\) ✓.
  • \(t = 5\): \(\widetilde{S}_D(5) = (5 \cdot 0.6065 + 21)/26 = 24.03/26 = 0.924\).

Beta process (식 6.4.2):

\[ \widetilde{S}_B(t) = \exp\left[-\frac{5 \cdot 0.1 t}{5 + 21}\right] = \exp\left(-\frac{0.5 t}{26}\right) \]

  • \(t = 5\): \(\widetilde{S}_B(5) = \exp(-2.5/26) = \exp(-0.0962) = 0.908\).

\(t \in [6, 7)\) 시점 (첫 사건 6 주 후):

Dirichlet:

\[ \widetilde{S}_D(t) = \frac{5 e^{-0.1t} + 17}{26} \cdot \frac{5 e^{-0.6} + 18}{5 e^{-0.6} + 17} \]

(첫 항: 사건 후 갱신된 위험집합 \(Y_2 = 17\), 곱 항: 사건 시점에서의 점프 보정).

Beta process:

\[ \widetilde{S}_B(t) = \exp\left\{-\frac{5 \cdot 0.6}{26} - \frac{5(0.1 t - 0.6)}{22}\right\} \cdot \left[1 - \frac{0.5 + 3}{26}\right] \]

(첫 두 항: 누적 hazard, 곱 항: 사건 점프).

직관 — Klein Figure 6.13 의 비교

4 개 곡선 비교:

  • KM (Ch.4): step function, 첫 사건 6 주에 0.857 로 점프, 23 주에 0.448.
  • Prior \(S_0(t) = e^{-0.1t}\): \(t = 6\) 에서 0.549, \(t = 23\) 에서 0.100.
  • Dirichlet \(\widetilde{S}_D\): KM 과 prior 사이의 가중 평균. 사건 시점에 점프 (KM 처럼) + 사이에는 매끄러운 감소 (prior 처럼).
  • Beta process \(\widetilde{S}_B\): prior 에 더 가까움 — sample path 가 매끄러움.

관찰:

  • Prior 가 매우 비관적 (\(t = 23\) 에서 0.10 vs KM 0.45) — prior 가 control 군 mean rate 기반이어서 6-MP 의 이점 반영 안 됨.
  • 하지만 \(c = 5\) 가 작으므로 prior 영향 제한적 — Bayes estimator 가 KM 에 더 가까움.
  • \(c\) 가 컸다면 (\(c = 50\)) prior 쪽으로 더 끌렸을 것.

→ Finite sample 에서 prior 가 안정화 효과. 큰 표본에서는 KM 으로 수렴.

3.5 Gibbs Sampler — 일반 Sampling Scheme 처리

직관 — 왜 Gibbs sampler 가 필요한가

Dirichlet + beta process 의 closed form (식 6.4.1·6.4.2) 은 right-censored 데이터만.

다른 sampling scheme (interval cens, double cens, 좌절단) 또는 regression (Ch.8 의 Cox PH) 에서 closed form 없음 → Monte Carlo (Gibbs sampler) 필요.

Gibbs sampler 의 장점:

  • 임의 sampling scheme + 임의 prior 처리.
  • Posterior 분포의 Monte Carlo sample 추출.
  • 점추정 + 신뢰구간 + posterior probability 모두 계산.

비용: 계산 시간 (1000~10000 회 반복).

정의: Gibbs Sampler (Gelfand-Smith 1990)

right-censored 데이터의 latent variable approach.

핵심 idea:

  • Censored 관측의 진짜 사건 시점을 잠재 변수 \(Z\) 로 처리.
  • Posterior 에서 \(Z\) 와 모수 \(\theta\) 를 번갈아 sample (각각 다른 conditional posterior 에서).
  • 충분히 반복 → posterior sample 추출.

알고리즘:

  • Step 0: 초기값 \(\theta^0 \sim\) prior (식 6.4.3 의 Dirichlet density).
  • Iteration \(i+1\):
    • Latent step: 각 censored \(\lambda_j\) 명에 대해 multinomial sampling — 사건이 어느 미래 구간에 있을지 (식 6.4.x 의 \(\rho_k\) 확률).
    • Parameter step: 갱신된 사건 분포로 새 \(\theta^{i+1} \sim\) Dirichlet posterior.
  • 수렴 후 식 6.4.4 의 posterior mean 으로 Bayes estimator.
Klein Example 6.4 의 Gibbs sampler 결과

24 개 grid intervals (Klein Table 6.5) — 사건 시점 + censoring 시점 모두 포함. 1000 회 Gibbs iteration:

\(j\) \((t_{j-1}, t_j]\) \(d_j\) \(\lambda_j\) \(\alpha_j\) Posterior \(\theta_j\) (SE)
1 \((0, 6^-]\) 0 0 2.256 0.0867 (0)
2 \((6^-, 6]\) 3 1 0.000 0.1154 (0)
3 \((6, 7^-]\) 0 0 0.261 0.0105 (0.0001)
4 \((7^-, 7]\) 1 0 0.000 0.0408 (0.0003)
17 \((22^-, 22]\) 1 0 0.000 0.0678 (0.0012)

관찰:

  • 사건 시점 (\(d_j > 0\)): Bayes 가 KM 의 점프와 유사.
  • Censoring 시점 (\(\lambda_j > 0\)): Bayes 가 prior 분배 보정.
  • SE 가 0 인 시점: 사건이 알려진 시점 (분배 불확실성 없음).
  • SE 가 큰 시점: censoring 분배의 Monte Carlo 변동 반영.

→ Closed form (식 6.4.1) 와 Gibbs (식 6.4.4) 결과 일치 — Gibbs 의 일반화 가능성을 확인.

4 Practical Notes (3 개)

Practical Note 1 — c 의 선택

Prior strength \(c\) 의 실무 가이드:

  • \(c\) 작음 (\(c \leq 1\)): prior 거의 무관, KM 과 거의 동일.
  • \(c\) 적당 (\(c = 5 \sim 10\)): 작은 표본 (\(n = 20 \sim 30\)) 안정화.
  • \(c\) 큼 (\(c \geq 50\)): 매우 강한 prior, 데이터를 압도.

실무 조언:

  • “내가 prior guess 에 얼마나 확신하는가” 를 가짜 표본 크기로 결정.
  • 표본 \(n\) 의 10~20% 가 적절 (prior 영향 명확하지만 데이터에 압도되지 않음).
Practical Note 2 — Excess Mortality 의 임상 보고

관행:

  1. 두 모형 (multiplicative + additive) 모두 추정.
  2. SMR (시간 일정 가정) 또는 B̂(t) (시간 변동) 둘 다 보고.
  3. Cancer 분야는 Â(t) + S^C(t) (relative survival) 표준.

주의:

  • 표준 인구 hazard \(\theta_j(t)\) 의 정확성이 핵심.
  • 환자 그룹과 표준 인구의 적절한 매칭 (나이·성별·인종) 중요.
Practical Note 3 — Bayesian 의 응용

언제 Bayesian 이 유용한가:

  1. 작은 표본 (n < 30): finite sample 안정화.
  2. Phase II 임상시험: 이전 phase I 결과를 prior 로.
  3. 신약 개발: 기존 약 effect 를 prior 로 (informative prior).
  4. 외부 데이터 통합: 전문가 의견 + 문헌 정보 + 표준 인구.

언제 Bayesian 이 부적절한가:

  1. 큰 표본 (\(n > 200\)): prior 영향 무시할 만 → KM 으로 충분.
  2. Prior 정보 없음: \(c \to 0\) 이 자연.
  3. 임상시험의 frequentist 보고 요구 (FDA/EMA 가이드라인).

5 Theoretical Notes

Theoretical Note 1 — c 와 효과적 표본 크기

Dirichlet prior 의 분산 공식:

\[ V[S(t)] = \frac{S_0(t)[1 - S_0(t)]}{c + 1} \]

이 분산은 “\(c + 1\) 명의 표본에서 추정한 비율의 분산” 과 같음:

\[ V[\widehat{p}] = \frac{p(1-p)}{n} \]

→ Prior 가 “효과적으로 \(c + 1\) 명의 가짜 표본” 가치.

Theoretical Note 2 — 점근 동치

regularity 조건 하에서 \(n \to \infty\):

\[ \widetilde{S}_D(t) \to \widehat{S}_{KM}(t), \quad \widetilde{S}_B(t) \to \widehat{S}_{KM}(t) \]

(in probability).

→ Bayes estimator 가 KM 의 일반화 — 작은 표본 + prior 정보 결합, 큰 표본에서 KM 회귀.

6 응용 분야

분야 도구 구체적 예시
종양학 relative survival § 6.3 multiplicative 5 년 cancer 생존률 (일반 인구 보정)
정신건강 mortality § 6.3 multiplicative Iowa psychiatric SMR ≈ 25
직업 노출 mortality § 6.3 additive 화학 노출 → 추가 hazard
HIV 코호트 추적 § 6.3 + § 6.4 표준 인구 + Bayesian smoothing
Phase II 임상시험 § 6.4 Bayesian n = 20 + prior (phase I 기반)
신약 개발 § 6.4 informative prior 기존 약 effect 를 prior 로
의료기기 PMA § 6.4 Bayesian FDA 의 informative prior 가이드

7 코드 예시

7.1 Step 1 — Excess Mortality (R)

library(survival)
library(KMsurv)
library(relsurv)  # relative survival 패키지

# Klein Example 6.3 — Iowa psychiatric
data(psych)

# 표준 인구 mortality table (1959-1961 Iowa)
ratetable <- survexp.us  # 또는 사용자 정의 life table

# Multiplicative — relative survival
fit <- rs.surv(Surv(time, event) ~ 1,
               data = psych,
               ratetable = ratetable,
               method = "ederer1")  # multiplicative

plot(fit, xlab = "Years on Study",
     ylab = "Relative Survival",
     main = "Klein Example 6.3 — Iowa Psychiatric")

# Additive — excess mortality
fit_excess <- rs.add(Surv(time, event) ~ 1,
                      data = psych,
                      ratetable = ratetable)
print(fit_excess)

7.2 Step 2 — Bayesian Dirichlet Process (R)

library(bayesSurv)

# Klein Example 6.4 — 6-MP
data(sixmp)  # 21 명 6-MP 데이터

# Dirichlet process prior + closed form
fit <- bayesBisurvreg(
  Surv(time, status) ~ 1,
  data = sixmp,
  prior = list(
    type = "dirichlet",
    S0 = function(t) exp(-0.1 * t),
    c = 5
  )
)

plot(fit, plot.type = "survival")

7.3 Step 3 — Gibbs Sampler 직접 구현 (Python)

import numpy as np
from scipy.stats import dirichlet

def gibbs_sampler_npmle(times, events, S_0_func, c, n_iter=1000, burn_in=100):
    """
    Klein § 6.4 Gibbs sampler for right-censored Bayesian NPMLE.
    times, events: 데이터.
    S_0_func: prior guess function.
    c: prior strength.
    """
    n = len(times)
    grid = np.unique(np.concatenate([[0], times, [np.inf]]))
    M = len(grid) - 2

    # 식 6.4.3 prior parameters
    alpha = c * np.array([S_0_func(grid[j-1]) - S_0_func(grid[j])
                          for j in range(1, M+2)])

    # d_j, lambda_j 계산
    d = np.zeros(M+1, dtype=int)
    lam = np.zeros(M+1, dtype=int)
    for t, e in zip(times, events):
        j = np.searchsorted(grid, t, side='right') - 1
        if e == 1:
            d[j] += 1
        else:
            lam[j] += 1

    # Initialize theta from prior
    theta = dirichlet.rvs(alpha)[0]

    samples = []
    for it in range(n_iter):
        # Step 1: latent variable (Z) sampling
        Z = np.zeros((M+1, M+1), dtype=int)
        for j in range(M+1):
            if lam[j] > 0:
                rho = theta[j+1:] / theta[j+1:].sum()
                Z[j+1:, j] = np.random.multinomial(lam[j], rho)

        # Step 2: theta sampling from Dirichlet posterior
        R = alpha + d + Z.sum(axis=1)
        theta = dirichlet.rvs(R)[0]

        if it >= burn_in:
            samples.append(theta.copy())

    # 식 6.4.4 — posterior mean
    theta_mean = np.mean(samples, axis=0)
    return grid, theta_mean

# Klein Example 6.4
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])

S_0 = lambda t: np.exp(-0.1 * t)
grid, theta = gibbs_sampler_npmle(times, events, S_0, c=5,
                                    n_iter=10000, burn_in=1000)

# Survival 추정
S_bayes = 1 - np.cumsum(theta)
print("Bayesian S(t):", list(zip(grid[1:], S_bayes)))

7.4 Step 4 — Python lifelines (Bayesian via PyMC3)

import pymc3 as pm
import numpy as np

# Klein Example 6.4 with PyMC3
with pm.Model() as model:
    # Prior: log-hazard ~ Normal
    log_h = pm.Normal('log_h', mu=np.log(0.1), sigma=1.0)
    h = pm.Deterministic('h', pm.math.exp(log_h))

    # Likelihood: Exponential survival
    obs = pm.Exponential('obs',
                          lam=h,
                          observed=times[events == 1])
    cens = pm.Potential('cens',
                         -h * times[events == 0].sum())

    trace = pm.sample(2000, tune=1000)

pm.summary(trace)
# h 의 posterior mean 으로 Bayesian 추정

8 핵심 takeaway

§ 6.3~6.4 의 6 가지 교훈
  1. 두 excess mortality 모형 (§ 6.3) — Multiplicative (\(\beta(t) = h_j/\theta_j\)) 는 “몇 배 빠른가”, additive (\(\alpha(t) = h_j - \theta_j\)) 는 “얼마 추가인가”. 같은 데이터에서 다른 질문에 답.

  2. B̂(t) 가 SMR 의 시간-변동 일반화 — Breslow 1975 의 단일 SMR 이 시간 일정 가정. B̂(t) 의 기울기가 시간별 \(\beta(t)\) — 환자 그룹 mortality 의 시간 패턴 정량화.

  3. Corrected survival \(S^C = \widehat{S}/S^*\) — Cancer relative survival 의 표준. “환자 고유의 mortality” 분리 — 표준 인구 효과 제거.

  4. Iowa psychiatric β ≈ 25 (Klein Example 6.3) — 정신질환자가 일반인보다 첫 2 년에 25 배 빠른 사망. 통계적 도구가 공중보건 의사결정의 근거.

  5. 두 conjugate Bayesian prior (§ 6.4) — Dirichlet process (Ferguson 1973) for \(S\) vs Beta process (Hjort 1990) for \(H\). 두 prior 모두 closed form Bayes estimator (식 6.4.1·6.4.2). \(c \to 0\) 또는 \(n \to \infty\) 에서 KM 으로 수렴.

  6. Gibbs sampler (Gelfand-Smith 1990) 의 일반화 — Closed form 없는 일반 sampling scheme (interval cens, regression) 에서도 Bayesian 추정 가능. Latent variable 처리 + Dirichlet posterior 반복.

9 관련 주제

선행 지식

후속 주제

  • § 6.5 — 10 Exercises 풀이
  • Ch.7 — Hypothesis Testing (log-rank · Wilcoxon)
  • Ch.8 — Cox PH (excess mortality 와의 연결)
  • Ch.10 — Aalen Additive Model (additive 의 회귀 일반화)

관련 개념

  • SMR (Standardized Mortality Ratio, Breslow 1975) — 인구 역학 표준
  • Relative survival (Cancer epidemiology) — § 6.3 multiplicative
  • Dirichlet process (Ferguson 1973) — Bayesian 비모수의 표준
  • Beta process (Hjort 1990) — hazard 의 Bayesian 표준
  • Gibbs sampler (Gelfand-Smith 1990) — Monte Carlo 표준
  • Conjugate prior 일반론 (Bayesian inference)

Subscribe

Enjoy this blog? Get notified of new posts by email: