Kwangmin Kim - Klein Ch.10 — Additive Hazards Regression Models

1 도입: 왜 가산 위험 모형인가

Cox 비례위험(PH) 모형은 공변량 효과를 상대적 승수로 표현한다:

\[ h(t \mid Z) = h_0(t) \exp(\beta^\top Z) \]

이 구조에서 \(\exp(\beta_k)\) 는 해당 공변량의 위험비(hazard ratio)로, “기저 위험이 몇 배 커지는가”를 알려준다. 그러나 공중보건 현장에서는 “이 위험 인자로 인해 단위 시간당 몇 명이 더 사망하는가?”라는 절대적 변화가 의사결정에 더 유용하다. 가산(additive) 위험 모형은 이 요구에 직접 답한다.

\[ h[t \mid Z(t)] = \beta_0(t) + \sum_{k=1}^{p} \beta_k(t) Z_k(t) \]

직관적으로, 각 \(\beta_k(t)\) 는 해당 공변량이 1단위 증가할 때 위험률에 더해지는 양이다. 비례위험 모형이 “배율”을 추정한다면, 가산 모형은 “덧셈량”을 추정한다.

비교 항목	Cox PH 모형	가산 위험 모형
공변량 작용 방식	곱셈(multiplicative)	덧셈(additive)
효과 척도	위험비 \(e^\beta\)	위험률 차이 \(\beta_k\)
공중보건 해석	“몇 배 위험”	“연간 몇 명 초과 사망”
계수의 시간 의존성	기본적으로 상수	Aalen: 시간 함수 / Lin-Ying: 상수
추정 원리	편우도(partial likelihood)	최소제곱 / 추정방정식

가산 모형의 핵심 장점은 효과의 시간적 변화를 자연스럽게 탐색할 수 있다는 것이다. \(B_k(t) = \int_0^t \beta_k(u) du\) 도표의 기울기 변화가 곧 효과의 시간적 변동을 시각적으로 보여준다. 이는 비례위험 가정이 의심될 때 대안이자 진단 도구로 활용된다.

단, 가산 모형에는 고유한 한계도 있다. 추정된 위험률이 음수가 될 수 있어 해석에 주의가 필요하며, 연속 공변량은 평균 중심화(centering)를 하지 않으면 기저 위험 추정치가 음수가 될 가능성이 커진다.

2 Aalen 비모수 가산 위험 모형

2.1 모형 구조

Aalen(1989) 모형은 회귀 계수 자체가 시간의 함수인 완전 비모수 모형이다:

\[ h[t \mid Z_j(t)] = \beta_0(t) + \sum_{k=1}^{p} \beta_k(t) Z_{jk}(t) \]

여기서 \(\beta_k(t)\) 는 시점 \(t\) 에서 공변량 \(Z_k\) 의 효과이다. 이 모형은 고정 공변량에 대해서도 효과가 시간에 따라 변할 수 있도록 허용한다. 예를 들어, 수술 직후 1년간은 특정 위험 인자의 효과가 크지만 이후 감소하는 패턴을 포착할 수 있다.

2.2 누적 회귀 함수 \(B_k(t)\)

\(\beta_k(t)\) 를 직접 추정하기는 어렵다. 대신 누적 회귀 함수(cumulative regression function)를 추정한다:

\[ B_k(t) = \int_0^t \beta_k(u) \, du, \quad k = 0, 1, \ldots, p \]

이는 Nelson-Aalen 추정량이 \(H(t) = \int_0^t h(u) du\) 를 추정하는 것과 동일한 논리이다. \(B_k(t)\) 의 기울기가 곧 \(\beta_k(t)\) 이므로:

직선: 효과가 시간에 걸쳐 일정하다 (Lin-Ying 모형이 적합)
기울기 증가: 효과가 점점 커진다
기울기 감소 / 수평: 효과가 약해지거나 사라진다

2.3 최소제곱 추정

\(n \times (p+1)\) 설계행렬 \(\mathbf{X}(t)\) 를 정의한다. \(i\) 번째 행은:

\[ X_i(t) = Y_i(t) \cdot (1, Z_{i1}(t), \ldots, Z_{ip}(t)) \]

여기서 \(Y_i(t)\) 는 위험 집합 지시 함수이다 (시점 \(t\) 에 관찰 중이면 1, 아니면 0). 즉, 관찰이 끝난 개체는 해당 행이 영벡터가 된다. 이 설계를 사용하면, 각 사건 시점에서의 “누가 죽었는가”를 \(n \times 1\) 벡터 \(\mathbf{I}(t)\) 로 나타낼 수 있다.

누적 회귀 함수의 최소제곱 추정량은:

\[ \hat{\mathbf{B}}(t) = \sum_{T_i \leq t} [\mathbf{X}^\top(T_i) \mathbf{X}(T_i)]^{-1} \mathbf{X}^\top(T_i) \mathbf{I}(T_i) \]

직관은 간단하다. 각 사건 시점 \(T_i\) 마다, “이 시점에 누가 위험에 노출되어 있었는가(설계행렬)”와 “누가 사건을 경험했는가(반응벡터)”를 사용해 시점별 최소제곱 해를 구하고, 이를 누적한다. 이는 선형 회귀의 \(\hat{\beta} = (X^\top X)^{-1} X^\top y\) 를 시점별로 반복 적용하는 것과 같다.

분산-공분산 행렬은:

\[ \widehat{\text{Var}}(\hat{\mathbf{B}}(t)) = \sum_{T_i \leq t} [\mathbf{X}^\top(T_i) \mathbf{X}(T_i)]^{-1} \mathbf{X}^\top(T_i) \mathbf{I}^{\text{diag}}(T_i) \mathbf{X}(T_i) \{[\mathbf{X}^\top(T_i) \mathbf{X}(T_i)]^{-1}\}^\top \]

여기서 \(\mathbf{I}^{\text{diag}}(t)\) 는 \(\mathbf{I}(t)\) 를 대각 행렬로 만든 것이다.

추정량의 존재 조건

\(\hat{\mathbf{B}}(t)\) 는 \(\mathbf{X}^\top(T_i) \mathbf{X}(T_i)\) 가 비특이(nonsingular)인 시점까지만 정의된다. 즉, 모든 공변량 수준에 최소 1명이 위험 집합에 남아 있어야 한다. 예를 들어 두 표본 문제에서 한 그룹의 마지막 환자가 사건을 경험하면, 그 이후 시점에서는 추정이 불가능하다.

2.4 두 표본 문제에서의 해석

공변량이 \(Z_1 = 1\) (표본 1) 또는 \(0\) (표본 2)인 경우:

\[ \hat{B}_0(t) = \sum_{\substack{T_i \leq t \\ i \in \text{sample 2}}} \frac{d_i}{N_2(T_i)} \]

이는 표본 2의 Nelson-Aalen 추정량 그 자체이다. 또한:

\[ \hat{B}_1(t) = \sum_{\substack{T_i \leq t \\ i \in \text{sample 1}}} \frac{d_i}{N_1(T_i)} - \sum_{\substack{T_i \leq t \\ i \in \text{sample 2}}} \frac{d_i}{N_2(T_i)} \]

이는 두 Nelson-Aalen 추정량의 차이이다. 즉 \(\hat{B}_1(t)\) 의 기울기는 표본 1이 표본 2에 비해 갖는 초과 위험률(excess hazard rate)을 나타낸다.

2.5 가설 검정

“공변량 \(k\) 의 효과가 전 시간에 걸쳐 0이다” (\(H_0: \beta_k(t) = 0, \forall t \leq \tau\))를 검정한다. 가중 검정통계량:

\[ \mathbf{U} = \sum_{T_i} \mathbf{W}(T_i) [\mathbf{X}^\top(T_i) \mathbf{X}(T_i)]^{-1} \mathbf{X}^\top(T_i) \mathbf{I}(T_i) \]

여기서 \(\mathbf{W}(t)\) 는 대각 가중 행렬이다. 가중치로는 위험 집합 크기 \(W_j(t) = n(t)\) 가 일반적으로 사용된다. 동시 검정 통계량:

\[ \chi^2 = \mathbf{U}_J^\top \mathbf{V}_J^{-1} \mathbf{U}_J \sim \chi^2_{|J|} \]

Aalen 가중 \(W(t) = \{\text{diag}[(\mathbf{X}^\top \mathbf{X})^{-1}]\}^{-1}\) 을 사용하면 두 표본 문제에서 \(U_1\) 이 로그순위 검정의 분자와 동일해진다. 그러나 분산 추정이 다르므로 p-값은 일반적으로 다르다. 로그순위 검정은 \(H_0\) 하에서 분산을 추정하고, Aalen 검정은 일반 모형에서 분산을 추정한다.

3 Lin-Ying 반모수 가산 위험 모형

3.1 모형 구조

Lin & Ying(1994, 1995) 모형은 Aalen 모형의 특수 사례로, 회귀 계수를 시간에 무관한 상수로 제한한다:

\[ h(t \mid Z_j(t)) = \alpha_0(t) + \sum_{k=1}^{p} \alpha_k Z_{jk}(t) \]

여기서 \(\alpha_0(t)\) 는 임의의 기저 위험 함수이고, \(\alpha_k\) 는 추정할 상수 모수이다. 이 모형은 Cox 모형과 구조적으로 대응한다:

	Cox PH	Lin-Ying 가산
모형	\(h_0(t) \exp(\beta^\top Z)\)	\(\alpha_0(t) + \alpha^\top Z\)
기저 함수	\(h_0(t)\) 임의	\(\alpha_0(t)\) 임의
공변량 효과	상수 \(\beta\) (곱셈)	상수 \(\alpha\) (덧셈)
추정 원리	편우도 최대화	추정방정식 풀이
해의 형태	반복법(Newton-Raphson)	닫힌 형식(closed-form)

Lin-Ying 모형의 가장 큰 실무적 장점은 추정량에 명시적 공식이 존재한다는 것이다. Cox 모형처럼 반복 수렴을 기다릴 필요가 없다.

3.2 추정량 유도

시점 \(t\) 에서 공변량의 위험 집합 평균:

\[ \bar{Z}(t) = \frac{\sum_{i=1}^{n} Z_i Y_i(t)}{\sum_{i=1}^{n} Y_i(t)} \]

이 평균으로부터의 편차를 누적하여 세 행렬을 구성한다:

\[ \mathbf{A} = \sum_{i=1}^{n} \sum_{j=1}^{i} (T_j - T_{j-1}) [Z_i - \bar{Z}(T_j)]^\top [Z_i - \bar{Z}(T_j)] \]

\[ \mathbf{B}^\top = \sum_{i=1}^{n} \delta_i [Z_i - \bar{Z}(T_i)] \]

\[ \mathbf{C} = \sum_{i=1}^{n} \delta_i [Z_i - \bar{Z}(T_i)]^\top [Z_i - \bar{Z}(T_i)] \]

직관적으로:

\(\mathbf{A}\) 는 공변량의 시간 가중 분산 (정보 행렬 역할)
\(\mathbf{B}\) 는 사건 발생자의 공변량 편차 합 (스코어 역할)
\(\mathbf{C}\) 는 사건 발생자의 공변량 편차 제곱합 (분산 추정 역할)

추정량과 분산:

\[ \hat{\alpha} = \mathbf{A}^{-1} \mathbf{B}^\top \]

\[ \widehat{\text{Var}}(\hat{\alpha}) = \mathbf{A}^{-1} \mathbf{C} \mathbf{A}^{-1} \]

이 형태는 선형 회귀의 샌드위치 분산 추정량과 구조적으로 동일하다. 모형이 정확히 맞으면 \(\mathbf{C} \approx \mathbf{A}\) 가 되어 분산이 \(\mathbf{A}^{-1}\) 로 단순화되지만, 일반적으로 샌드위치 형태가 로버스트하다.

3.3 가설 검정

개별 검정은:

\[ Z = \frac{\hat{\alpha}_j}{\sqrt{\hat{V}_{jj}}} \xrightarrow{d} N(0,1) \]

동시 검정은:

\[ \chi^2 = \hat{\alpha}_J^\top \hat{\mathbf{V}}_J^{-1} \hat{\alpha}_J \sim \chi^2_{|J|} \]

3.4 기저 누적 위험 함수 추정

\(\hat{\alpha}\) 를 구한 후, 기저 누적 위험 \(A_0(t) = \int_0^t \alpha_0(u) du\) 를 추정할 수 있다:

\[ \hat{A}_0(t) = \sum_{T_i \leq t} \frac{\delta_i}{Y_\cdot(T_i)} - \sum_{j=1}^{k} \sum_{i=j}^{n} \hat{\alpha}^\top Z_i \frac{T_j - T_{j-1}}{Y_\cdot(T_j)} \]

첫째 항은 통상적인 Nelson-Aalen 추정량이고, 둘째 항은 공변량 효과를 보정한 부분이다.

4 예제: 유방암 두 표본 비교

면역과산화효소(immunoperoxidase) 양성 vs 음성 유방암 환자의 사망 시간 비교 (Klein, Example 10.1).

Aalen 모형 결과:

\(\hat{B}_0(t)\): 양성(baseline) 환자의 누적 위험. 20-90개월 구간에서 기울기가 거의 일정하여 \(\beta_0(t) \approx 0.009\) 사망/월
\(\hat{B}_1(t)\): 음성 환자의 누적 초과 위험. 기울기가 양수이므로 음성 환자의 위험이 더 높다

Aalen 검정: \(U_1 = -4.187\), \(\sigma_{11} = 6.04\) → \(Z = -1.704\), p = 0.088 (유의하지 않음)

Lin-Ying 모형 결과: \(\hat{\alpha} = 0.00803\), SE = 0.00471, \(Z = 1.704\), p = 0.088

두 모형이 거의 동일한 결론을 내린다. 이는 \(B_1(t)\) 도표가 거의 직선이어서, 효과가 시간에 걸쳐 일정하기 때문이다. 로그순위 검정(p = 0.019)과 차이가 나는 이유는 분산 추정 방식의 차이 때문이다: 로그순위는 \(H_0\) 가정 하에서 분산을 구하고, 가산 모형 검정은 일반 모형에서 구한다.

5 예제: 후두암 다중 공변량

90명 남성 후두암 환자, 4가지 병기(Stage I-IV) + 진단 시 연령 (평균 64.11세에 중심화) (Klein, Example 10.2).

Aalen ANOVA 결과 (가중치 = 위험 집합 크기):

공변량	\(\chi^2\)	df	p-값
\(Z_1\): Stage II	0.14	1	0.711
\(Z_2\): Stage III	3.16	1	0.075
\(Z_3\): Stage IV	6.61	1	0.010
\(Z_4\): 연령	0.28	1	0.591

Lin-Ying 결과:

공변량	\(\hat{\alpha}\)	SE	p-값
Stage II	0.01325	0.0435	0.761
Stage III	0.07654	0.0442	0.083
Stage IV	0.37529	0.1403	0.008
연령	0.00256	0.0022	0.240

해석: Stage IV 환자는 Stage I 환자에 비해 연간 약 0.38건의 초과 사망 위험을 갖는다. 이것이 가산 모형의 힘이다 — Cox 모형은 “위험이 \(e^{0.38}\) 배”라고 말하지만, 가산 모형은 직접적으로 “연간 0.38건 더 사망한다”고 말해준다. Stage II는 Stage I과 유의한 차이가 없으며, Stage III는 경계적이다.

병기 효과의 동시 검정: Aalen \(\chi^2 = 9.18\) (p = 0.027), Lin-Ying \(\chi^2 = 9.84\) (p = 0.020).

커널 평활(kernel smoothing) 추정에 따르면, Stage IV의 초과 위험은 진단 후 첫 2년간 약 0.4 사망/년 수준이며, Stage III는 약 2.5년 후 초과 위험이 소멸된다. 이 시간적 패턴은 Cox 모형의 단일 위험비로는 포착할 수 없다.

6 계수과정 이론 기반 유도

6.1 Aalen 추정량의 마팅게일 유도

계수과정 \(N_j(t)\) 의 분해:

\[ d\mathbf{N}(t) = \mathbf{X}(t) \boldsymbol{\beta}(t) dt + d\mathbf{M}(t) \]

여기서 \(\mathbf{M}(t)\) 는 마팅게일 벡터(통계적 잡음)이다. 잡음을 0으로 놓고 일반화 역행렬 \(\mathbf{X}^-(t)\) 를 적용하면:

\[ \hat{\mathbf{B}}(t) = \int_0^t \mathbf{X}^-(u) \, d\mathbf{N}(u) = \sum_{T_j \leq t} \delta_j \mathbf{X}^-(T_j) \]

이는 Nelson-Aalen 추정량의 다변량 일반화이다. 일반화 역행렬로 \(\mathbf{X}^-(t) = [\mathbf{X}^\top(t)\mathbf{X}(t)]^{-1}\mathbf{X}^\top(t)\) 를 사용하면 앞서의 최소제곱 추정량을 얻는다.

6.2 Lin-Ying 추정량의 프로파일 구성

기저 위험 \(\alpha_0(t)\) 가 알려져 있다면, Nelson-Aalen 유사 추정:

\[ \hat{A}_0(t) = \int_0^t \frac{\sum_{i=1}^{n} \{dN_i(u) - Y_i(u) \alpha^\top Z_i(u) du\}}{\sum_{i=1}^{n} Y_i(u)} \]

이를 스코어 방정식에 대입하면(프로파일 구성):

\[ \mathbf{U}(\alpha) = \sum_{i=1}^{n} \int_0^\infty [Z_i(t) - \bar{Z}(t)] [dN_i(t) - Y_i(t) \alpha^\top Z_i(t) dt] \]

\(\mathbf{U}(\alpha) = 0\) 을 풀면 닫힌 형식의 \(\hat{\alpha}\) 를 얻는다. 이 과정은 Cox 모형에서 Breslow 추정을 프로파일 우도에 대입하는 것과 정확히 대응한다.

7 실무적 고려사항

7.1 가산 vs 비례: 언제 어떤 모형을 쓸 것인가

상황	권장 모형
상대적 효과(HR)가 관심	Cox PH
절대적 초과 위험이 관심	가산 위험
PH 가정이 의심됨	Aalen (진단 도구로 활용)
효과의 시간 패턴 탐색	Aalen \(B_k(t)\) 도표
간결한 요약 통계 필요	Lin-Ying
음수 위험률 문제 회피	Cox PH (항상 비음)

7.2 음수 위험률 문제

가산 모형에서는 추정된 위험률 \(\hat{h}(t \mid Z) = \hat{\beta}_0(t) + \sum \hat{\beta}_k(t) Z_k\) 가 음수가 될 수 있다. 이는 모형의 한계이지 오류가 아니다. 대응 방법:

연속 공변량을 평균 중심화하면 기저 위험이 “평균적 개체”에 대한 것이 되어 음수 가능성이 줄어든다
생존 함수 추정 시 누적 위험의 단조성이 보장되지 않으므로 주의가 필요하다

7.3 Aalen 모형의 특이한 성질

공변량 독립성 하 축소: 한 공변량이 나머지와 독립이면, 그 공변량을 제거해도 나머지 회귀 함수는 동일하다 (기저 함수만 변함). Cox 모형에서는 이 성질이 성립하지 않는다
측정 오차 강건성: 공변량에 가법적 정규 오차가 있으면, 선형 모형은 보존되고 회귀 함수는 상수배만큼 축소된다. Cox 모형에서는 모형 형태 자체가 변한다
초과 사망 모형과의 관계: \(Z_j(t) = \theta_j(t)\) (알려진 참조 위험)으로 놓으면, 6.3절의 초과 사망 모형이 특수 사례가 된다

8 코드 예시

8.1 Step 1: 순수 Python으로 Aalen 추정 원리 구현

import numpy as np

# 두 표본 문제: 각 그룹의 Nelson-Aalen 차이 = B_1(t)
times_g1 = np.array([3, 5, 7, 12, 15])   # 그룹 1 사건 시점
times_g2 = np.array([2, 6, 8, 10, 20])   # 그룹 2 사건 시점
cens_g1 = np.array([1, 1, 0, 1, 1])      # 사건 지시자
cens_g2 = np.array([1, 1, 1, 0, 1])

n1_total, n2_total = len(times_g1), len(times_g2)

# 모든 사건 시점 정렬
all_times = np.sort(np.unique(np.concatenate([
    times_g1[cens_g1 == 1], times_g2[cens_g2 == 1]
])))

# 누적 회귀 함수 B_1(t) = NA(그룹1) - NA(그룹2) 계산
B1 = 0.0
B1_values = []
for t in all_times:
    n1_at_risk = np.sum(times_g1 >= t)
    n2_at_risk = np.sum(times_g2 >= t)
    d1 = np.sum((times_g1 == t) & (cens_g1 == 1))
    d2 = np.sum((times_g2 == t) & (cens_g2 == 1))
    if n1_at_risk > 0:
        B1 += d1 / n1_at_risk
    if n2_at_risk > 0:
        B1 -= d2 / n2_at_risk
    B1_values.append(B1)

print("시점별 누적 초과 위험 B_1(t):")
for t, b in zip(all_times, B1_values):
    print(f"  t={t:3d}: B_1 = {b:.4f}")
# B_1(t) > 0 이면 그룹 1의 위험이 그룹 2보다 높다

8.2 Step 2: lifelines를 이용한 Aalen 가산 모형 적합

from lifelines import AalenAdditiveFitter
import pandas as pd

# 후두암 데이터 예시 (단순화)
data = pd.DataFrame({
    'T': [0.6, 1.3, 2.4, 3.2, 0.4, 1.5, 3.0, 4.1, 0.2, 0.9],
    'E': [1, 1, 0, 1, 1, 1, 1, 0, 1, 1],
    'stage_III': [0, 0, 0, 0, 1, 1, 1, 1, 0, 0],
    'stage_IV':  [0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
})

aaf = AalenAdditiveFitter(coef_penalizer=0.1)
aaf.fit(data, duration_col='T', event_col='E')

# 누적 회귀 함수 도표 — 기울기가 효과의 시간적 변화를 보여준다
aaf.plot()

# 기울기가 직선이면 Lin-Ying 상수 계수 모형이 적합함을 시사
print(aaf.summary)

9 관련 주제

선행 지식

Ch.10 시리즈

후속 주제

관련 개념

초과 사망 모형 — 가산 모형의 특수 사례
현대 ML 생존 분석 확장 — 딥러닝 기반 가산 위험