1 § 10.1 Introduction — 왜 또 다른 회귀 모형인가
1.1 비례위험 모형의 한계
Ch.8-9에서 다룬 Cox 비례위험(PH) 모형은 다음 형태였다:
\[ h(t \mid Z) = h_0(t) \exp\left(\sum_{k=1}^{p} \beta_k Z_k\right) \]
이 모형의 핵심 가정은 두 가지이다:
- 곱셈적 작용: 공변량은 기저 위험 \(h_0(t)\) 에 곱해져 작용한다
- 상수 계수: \(\beta_k\) 는 시간에 무관하다 (PH 가정)
이 가정이 깨지는 상황을 만나면 어떻게 할 것인가? Ch.9에서 두 가지 대응책을 살펴봤다:
- 시간의존 공변량: \(Z_k(t) = Z_k \times g(t)\) 형태로 시간 의존성을 명시적으로 모형화
- 층화: 비례위험을 만족시키지 못하는 변수를 stratum으로 분리
그러나 이 접근들은 여전히 “곱셈적” 틀 안에 머문다. 만약 효과 자체가 본질적으로 덧셈적이라면, 그리고 그 덧셈량이 시간에 따라 변한다면, 어떻게 해야 하는가?
1.2 가산 위험 모형의 일반 형태
Klein § 10.1은 다음 형태를 제시한다:
\[ h[t \mid Z(t)] = \beta_0(t) + \sum_{k=1}^{p} \beta_k(t) Z_k(t) \tag{식 10.1.1} \]
이 식이 가산 위험 모형의 출발점이며, 두 가지 주요 변형이 파생된다:
| 모형 | 계수 형태 | 추정 원리 | 다루는 절 |
|---|---|---|---|
| Aalen (1989) | \(\beta_k(t)\) — 시간 함수 | 최소제곱 (counting process) | § 10.2 |
| Lin-Ying (1995) | \(\beta_k\) — 상수 | 추정방정식 (score-mimic) | § 10.3 |
본 포스트는 § 10.1의 도입과 § 10.2의 Aalen 모형을 깊이 다룬다. Lin-Ying 모형은 별도 포스트에서 다룬다.
1.3 곱셈 vs 덧셈: 직관적 차이
비례위험과 가산 위험의 차이를 구체적 숫자로 보면 명확하다. 기저 위험률 \(h_0(t) = 0.01\) (월간 1% 사망)이고, 어떤 위험 인자 \(Z = 1\) 이 미치는 효과를 비교해보자.
| 모형 | 효과 척도 | 해석 |
|---|---|---|
| Cox: \(\beta = \ln 2\) | \(\text{HR} = 2\) | “위험이 2배 커진다” → \(h(t) = 0.02\) |
| 가산: \(\beta(t) = 0.005\) | 절대 차이 | “위험이 월 0.5%p 증가한다” → \(h(t) = 0.015\) |
곱셈 모형에서는 \(h_0(t)\) 가 작을 때 절대 증가량도 작다. 반면 가산 모형에서는 기저 위험과 무관하게 일정한 절대 증가량을 갖는다. 공중보건 의사결정(“연간 추가 사망자 수”)에는 가산 모형이 직관적이다.
2 § 10.2 Aalen 비모수 가산 위험 모형
2.1 모형 설정과 데이터 구조
데이터는 통상의 생존 분석 형태를 따른다: \((T_j, \delta_j, Z_j(t))\), \(j = 1, \ldots, n\).
- \(T_j\): 관측 시간 (사건 또는 우중도절단)
- \(\delta_j\): 사건 지시자 (1 = 사건, 0 = 절단)
- \(Z_j(t)\): \(p\)-차원 (시간의존 가능) 공변량 벡터
각 개체에 대해 위험 집합 지시자를 정의한다:
\[ Y_j(t) = \begin{cases} 1 & \text{개체 } j \text{ 가 시점 } t \text{ 에 관찰 중} \\ 0 & \text{그렇지 않음} \end{cases} \]
좌절단이 있다면 \(Y_j(t) = 1\) 은 진입 시점부터 종료 시점까지만이다. 우중도절단만 있다면 \(t \leq T_j\) 일 때 1이다.
2.2 핵심 모형: 식 (10.2.1)
개체 \(j\) 의 조건부 위험률을 다음과 같이 모형화한다:
\[ h[t \mid Z_j(t)] = \beta_0(t) + \sum_{k=1}^{p} \beta_k(t) Z_{jk}(t) \tag{식 10.2.1} \]
핵심 차이점: \(\beta_k(t)\) 가 시간의 함수이다. 이는 두 가지 의미를 갖는다:
- 유연성: 고정 공변량이라도 그 효과가 시간에 따라 변할 수 있다. 예: 수술 후 1년간만 합병증 위험이 큰 인자
- 추정의 어려움: \(\beta_k(t)\) 자체는 무한 차원이므로 직접 추정이 비실용적
2.3 누적 회귀 함수 \(B_k(t)\) 라는 우회로
\(\beta_k(t)\) 를 점별로 추정하는 대신, 적분된 형태를 추정한다:
\[ B_k(t) = \int_0^t \beta_k(u) \, du, \quad k = 0, 1, \ldots, p \tag{식 10.2.2} \]
이 우회는 Ch.4의 Nelson-Aalen 추정량과 동일한 논리이다. 위험률 \(h(t)\) 자체보다 누적 위험 \(H(t) = \int_0^t h(u) du\) 가 추정하기 쉬운 것과 같다. 왜 그런가?
연속 시점 \(t\) 에서 위험률 \(\beta_k(t)\) 를 추정하려면 그 시점 근방의 데이터가 필요하다. 그런데 생존 데이터는 사건 시점에서만 정보를 준다 — 사건 사이의 빈 구간은 정보가 없다. 적분 형태 \(B_k(t)\) 는 사건 시점들의 계단 함수로 추정할 수 있다. 각 사건 시점에서 점프(jump)가 일어나고, 사건 사이에서는 일정하다. 이는 데이터의 자연스러운 구조와 정확히 부합한다.
기울기 \(\beta_k(t) \approx \Delta B_k(t) / \Delta t\) 는 사후적으로 도표에서 읽거나 커널 평활로 얻을 수 있다.
2.4 도표 진단: \(B_k(t)\) 가 말하는 것
추정된 \(\hat{B}_k(t)\) 의 도표 형태는 효과의 시간적 동태를 직접 보여준다:
| 도표 형태 | 의미 | 함의 |
|---|---|---|
| 직선 (일정 기울기) | \(\beta_k(t)\) 가 상수 | Lin-Ying 상수 모형이 적합 |
| 볼록 (기울기 증가) | 효과가 시간에 따라 강해짐 | 잠복기 후 효과 발현 |
| 오목 (기울기 감소) | 효과가 약해짐 | 초기 충격이 점차 소멸 |
| 수평 → 상승 → 수평 | 특정 구간에만 효과 | 변환점(change-point) |
| 거의 수평 | 효과 없음 | \(\beta_k(t) \approx 0\) |
이 진단은 비모수 모형이기 때문에 가능하다. 사용자가 “효과의 시간 형태”를 가정하지 않고 데이터로부터 읽어내는 것이다.
2.5 최소제곱 추정량 — 식 (10.2.3)
추정량 유도의 핵심 도구는 설계행렬 \(\mathbf{X}(t)\) 이다. 이는 \(n \times (p+1)\) 행렬로, \(i\) 번째 행은:
\[ X_i(t) = Y_i(t) \cdot (1, Z_{i1}(t), Z_{i2}(t), \ldots, Z_{ip}(t)) \]
즉:
- 시점 \(t\) 에 위험 집합에 있는 개체: \((1, Z_{i1}(t), \ldots, Z_{ip}(t))\)
- 위험 집합에 없는 개체: \((0, 0, \ldots, 0)\) — 영벡터
또한 \(\mathbf{I}(t)\) 는 \(n \times 1\) 벡터로, 시점 \(t\) 에 사건이 발생한 개체에서만 1, 나머지는 0이다.
이 설정에서 \(B(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) \tag{식 10.2.3} \]
이 식이 무서워 보이지만, 본질은 사건 시점마다 미니 선형회귀를 수행하고 그 결과를 누적하는 것이다.
각 사건 시점 \(T_i\) 에서:
- 공변량 정보: \(\mathbf{X}(T_i)\) — 그 시점에 누가 위험 집합에 있고, 그들의 공변량 값은 무엇인가
- 반응 정보: \(\mathbf{I}(T_i)\) — 그 시점에 누가 사건을 경험했는가 (보통 1명, 동순위 시 여러 명)
- 회귀: \([\mathbf{X}^\top \mathbf{X}]^{-1} \mathbf{X}^\top \mathbf{I}\) — 통상의 OLS 공식 그 자체
각 시점의 OLS 결과는 그 시점의 점프 크기 \(\Delta \hat{B}_k(T_i)\) 를 준다. 이를 시간에 걸쳐 누적하면 \(\hat{B}_k(t)\) 의 계단 함수가 된다.
2.6 분산 추정량 — 식 (10.2.4)
\(\hat{\mathbf{B}}(t)\) 의 분산-공분산 행렬은:
\[ \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{D}}(T_i) \mathbf{X}(T_i) \{[\mathbf{X}^\top(T_i) \mathbf{X}(T_i)]^{-1}\}^\top \tag{식 10.2.4} \]
여기서 \(\mathbf{I}^{\text{D}}(t)\) 는 \(\mathbf{I}(t)\) 를 대각 원소로 갖는 대각 행렬이다.
이 분산은 통상 OLS의 \(\sigma^2 (X^\top X)^{-1}\) 과 다른 형태인데, 사건의 발생 자체가 베르누이 시행이고 그 분산이 평균(=1) × (1−평균≈0)에 비례한다는 점에서 유래한다 (마팅게일 분산 구조). 점근적으로 \(\hat{\mathbf{B}}(t)\) 는 정규 분포를 따른다.
2.7 신뢰구간 — 식 (10.2.10)
지점별(pointwise) 95% 신뢰구간은 단순한 정규 근사:
\[ \hat{B}_k(t) \pm 1.96 \sqrt{\widehat{\text{Var}}[\hat{B}_k(t)]} \tag{식 10.2.10} \]
이는 Ch.4의 Nelson-Aalen 신뢰구간과 정확히 동일한 구조이다. 변환된 척도(예: log-log) CI도 같은 방식으로 구성할 수 있다.
2.8 추정 가능 영역의 한계: \(\tau\)
\(\hat{\mathbf{B}}(t)\) 는 \(\mathbf{X}^\top(T_i) \mathbf{X}(T_i)\) 가 비특이(invertible)인 시점까지만 정의된다.
설계행렬이 특이가 되는 상황:
- 두 표본 문제: 한 그룹의 마지막 환자가 사건 또는 절단을 경험한 시점
- 다군 문제: 어느 한 군의 위험 집합이 비는 시점
- 연속 공변량: 위험 집합 내 공변량 변동이 사라지는 시점
이 한계 시점을 \(\tau\) 라 부르며, 모든 결과는 \(t \leq \tau\) 에서만 유효하다. Klein 본문에서 유방암 예제는 \(\tau = 144\) 개월, 후두암 예제는 \(\tau = 4.4\) 년이다.
3 두 표본 문제로 본 추정량의 정체
가산 모형의 단순한 경우를 직접 풀어보면 추정량의 의미가 분명해진다. 공변량이 하나뿐:
\[ Z_{j1} = \begin{cases} 1 & \text{표본 1} \\ 0 & \text{표본 2} \end{cases} \]
이 경우 \(\mathbf{X}^\top(t) \mathbf{X}(t)\) 와 그 역행렬은 매우 간단한 \(2 \times 2\) 행렬이다 (Klein 본문 식 참조). 이를 식 (10.2.3)에 대입하면 놀라운 결과가 나온다.
3.1 \(\hat{B}_0(t)\) 는 표본 2의 Nelson-Aalen 추정량
\[ \hat{B}_0(t) = \sum_{\substack{T_i \leq t \\ i \in \text{sample 2}}} \frac{d_i}{N_2(T_i)} \tag{식 10.2.5} \]
이는 정확히 표본 2(기저군) 만의 Nelson-Aalen 추정량이다. 즉, \(\beta_0(t)\) 는 “공변량 = 0”인 그룹의 위험률을 의미하므로, 그 누적은 표본 2의 누적 위험과 같다.
3.2 \(\hat{B}_1(t)\) 는 두 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)} \tag{식 10.2.6} \]
이는 두 그룹의 누적 위험률 차이이다. 시점 \(t\) 에서의 기울기 \(\beta_1(t)\) 는 표본 1이 표본 2 대비 갖는 순간 초과 위험률(excess hazard rate)이 된다.
식 (10.2.6)이 말하는 것은 단순하다: “가산 모형의 회귀 계수는 두 그룹의 위험 차이를 누적한 것이다.” 이는 비례위험 모형과 근본적으로 다르다:
- Cox HR: 두 그룹의 위험 비율이 일정하다고 가정 → \(\log[h_1(t)/h_2(t)] = \beta_1\)
- Aalen: 두 그룹의 위험 차이의 시간 패턴을 직접 추적 → \(h_1(t) - h_2(t) = \beta_1(t)\)
따라서 \(\hat{B}_1(t)\) 도표는 “두 그룹 간 누적 초과 위험”의 시각화 그 자체이다. 도표가 점점 평평해지면 “초과 위험이 사라졌다”, 가팔라지면 “초과 위험이 누적되고 있다”로 직접 해석된다.
3.3 분산도 직관적
식 (10.2.7)-(10.2.9)에 따르면:
\[ \widehat{\text{Var}}[\hat{B}_0(t)] = \sum_{\substack{T_i \leq t \\ i \in \text{sample 2}}} \frac{d_i}{N_2(T_i)^2} \]
\[ \widehat{\text{Var}}[\hat{B}_1(t)] = \sum_{\substack{T_i \leq t \\ i \in \text{sample 1}}} \frac{d_i}{N_1(T_i)^2} + \sum_{\substack{T_i \leq t \\ i \in \text{sample 2}}} \frac{d_i}{N_2(T_i)^2} \]
이는 정확히 두 Nelson-Aalen 분산의 합이다. 두 독립 추정량의 차이에 대한 분산이 각 분산의 합인 것과 일치한다 (\(\text{Var}(A - B) = \text{Var}(A) + \text{Var}(B)\), \(A \perp B\)).
4 \(\beta_k(t)\) 의 평활 추정
도표를 보고 효과의 시간 패턴을 시각적으로 판단하는 것을 넘어, \(\beta_k(t)\) 자체의 추정이 필요할 때가 있다. Ch.6의 커널 평활(kernel smoothing) 기법을 그대로 가져온다.
기본 아이디어: \(\hat{B}_k(t)\) 는 사건 시점에서 점프하는 계단 함수이다. 점프의 크기 \(\Delta \hat{B}_k(T_i)\) 를 위험률의 “이산 관측”으로 간주하고, 커널을 씌워 부드러운 함수로 만든다:
\[ \hat{\beta}_k(t) = \int K\left(\frac{t - u}{b}\right) \frac{1}{b} \, d\hat{B}_k(u) = \sum_{T_i} \frac{1}{b} K\left(\frac{t - T_i}{b}\right) \Delta \hat{B}_k(T_i) \]
여기서 \(K\) 는 커널 함수(예: biweight), \(b\) 는 대역폭(bandwidth)이다. 이는 식 (6.2.4)의 \(\Delta \tilde{H}(t)\) 자리에 \(\Delta \hat{B}_k(t)\) 를 대입하는 것과 같다.
후두암 예제(§10.2 본문)에서 biweight 커널, 대역폭 1년을 사용한 평활 결과는 다음을 보여준다:
- Stage II: 전 구간 거의 0 (효과 없음)
- Stage III: 초기 1-2년 양수, 약 2.5년 후 0 근방으로 수렴
- Stage IV: 첫 2년간 약 0.4 사망/년, 이후 감소
이는 단일 위험비로는 절대 포착할 수 없는 시간적 패턴이다.
5 가설 검정: \(\beta_k(t) = 0\) 인가
5.1 검정의 목표
가설은 “공변량 \(k\) 의 효과가 모든 시점에서 0이다”:
\[ H_0: \beta_k(t) = 0, \quad \forall t \leq \tau, \quad k \in K \]
이는 단일 시점에 대한 검정이 아니라 함수 전체에 대한 검정이다. Ch.7에서 본 가중 로그순위 검정의 Aalen 모형 버전이다.
5.2 검정 통계량 — 식 (10.2.11), (10.2.13)
가중 검정 통계량 벡터:
\[ \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) \tag{식 10.2.11} \]
여기서 \(\mathbf{W}(t)\) 는 대각 가중 행렬이다. \(\mathbf{U}\) 는 식 (10.2.3)의 추정량 증분에 가중치를 곱해 합한 형태이며, \(H_0\) 하에서 평균이 0이다.
공분산 행렬 \(\mathbf{V}\) 는 식 (10.2.12)로 주어지고, 동시 검정 통계량은:
\[ \chi^2 = \mathbf{U}_J^\top \mathbf{V}_J^{-1} \mathbf{U}_J \sim \chi^2_{|J|} \tag{식 10.2.13} \]
5.3 가중 함수 선택
| 가중치 | 형태 | 특성 |
|---|---|---|
| Aalen 제안 | \(W(t) = \{\text{diag}[(\mathbf{X}^\top \mathbf{X})^{-1}]\}^{-1}\) | 부분가설마다 다름 → 일관성 문제 |
| 위험 집합 크기 | \(W_j(t) = N(t)\) | 모든 부분가설에 공통 (권장) |
| 단위 가중 | \(W_j(t) = 1\) | 단순하지만 효율 낮음 |
Klein은 공통 가중치를 권장한다. Aalen 가중은 두 표본 문제에서 우아한 결과를 주지만, 다중 공변량에서 부분가설(예: \(\beta_1 = 0\) 만 검정 vs \(\beta_1 = \beta_2 = 0\) 동시 검정)마다 통계량이 일관되지 않을 수 있다.
5.4 로그순위 검정과의 관계
두 표본 문제에서 Aalen 가중치 \(W_{22}(t) = [N_1^{-1}(t) + N_2^{-1}(t)]^{-1} = N_1 N_2 / (N_1 + N_2)\) 를 사용하면:
\[ U_1 = \sum_{i=1}^n \delta_i \left( Z_{i1} - \frac{N_1(T_i)}{N_1(T_i) + N_2(T_i)} \right) = D_1 - \sum_i \frac{\delta_i N_1(T_i)}{N_1(T_i) + N_2(T_i)} \]
이는 정확히 로그순위 통계량의 분자이다 (식 7.3.3). 즉 Aalen 검정의 통계량 자체는 로그순위와 같다.
Aalen 검정과 로그순위 검정의 차이는 분산 추정에 있다:
- 로그순위: \(H_0\) (두 그룹 동일) 가정 하에서 분산 추정 — 더 작은 분산 → 더 강한 검정력
- Aalen: 일반 모형 하에서 분산 추정 (두 그룹 동일성 가정 안 함) — 더 큰 분산 → 더 보수적
따라서 같은 \(U_1\) 에 대해서도 p-값이 다르게 나온다. 유방암 예제(Klein):
| 검정 | \(U_1\) | \(\sqrt{V_{11}}\) | \(Z\) | p-값 |
|---|---|---|---|---|
| 로그순위 | -4.187 | √3.191 | -2.344 | 0.019 |
| Aalen | -4.187 | √6.040 | -1.704 | 0.088 |
같은 데이터에서 결론이 갈린다. Aalen이 더 robust한 분산을 사용하므로 결론을 더 신중하게 받아들여야 할 때 유용하다.
6 마팅게일 기반 유도 — 왜 식 (10.2.3)이 옳은가
6.1 계수과정 모형
각 개체 \(j\) 에 대해 계수과정(counting process) \(N_j(t) = I(T_j \leq t, \delta_j = 1)\) 을 정의한다. 가산 위험 모형 하에서 \(N_j(t)\) 의 강도(intensity) 함수는:
\[ \lambda_j(t) = Y_j(t) \cdot [Y_j(t), Y_j(t) Z_{j1}(t), \ldots, Y_j(t) Z_{jp}(t)] \cdot \boldsymbol{\beta}(t) \]
이는 정확히 설계행렬 \(\mathbf{X}(t)\) 의 \(j\) 행에 \(\boldsymbol{\beta}(t)\) 를 곱한 것이다.
6.2 마팅게일 분해
계수과정의 핵심 정리(Doob-Meyer 분해)에 의해:
\[ \mathbf{M}(t) = \mathbf{N}(t) - \int_0^t \mathbf{X}(u) \boldsymbol{\beta}(u) \, du \]
는 평균 0 마팅게일이다. 이를 미분 형태로 쓰면:
\[ d\mathbf{N}(t) = \mathbf{X}(t) \boldsymbol{\beta}(t) \, dt + d\mathbf{M}(t) \tag{식 10.2.14} \]
이는 선형 회귀 형태이다: 신호 \(\mathbf{X}(t)\boldsymbol{\beta}(t)dt\) + 잡음 \(d\mathbf{M}(t)\).
6.3 형식적 풀이
잡음을 0으로 놓고 \(\mathbf{B}(t) = \int_0^t \boldsymbol{\beta}(u) du\) 를 풀면:
\[ \hat{\mathbf{B}}(t) = \int_0^t \mathbf{X}^-(u) \, d\mathbf{N}(u) \tag{식 10.2.15} \]
여기서 \(\mathbf{X}^-(t)\) 는 \(\mathbf{X}(t)\) 의 일반화 역행렬이다. 일반화 역행렬은 유일하지 않지만, 가장 자연스러운 선택은 OLS의 의사역행렬:
\[ \mathbf{X}^-(t) = [\mathbf{X}^\top(t) \mathbf{X}(t)]^{-1} \mathbf{X}^\top(t) \]
이를 대입하면 정확히 식 (10.2.3)을 얻는다. 즉 Aalen 추정량은 마팅게일 회귀의 최소제곱 해이다. 이는 Nelson-Aalen 추정량이 단변량 위험률 회귀의 해인 것의 자연스러운 다변량 일반화이다.
6.4 분산 유도
\((\hat{\mathbf{B}} - \mathbf{B})(t) = \int_0^t \mathbf{X}^-(u) \, d\mathbf{M}(u)\) 는 마팅게일에 대한 확률적분이다. 이의 예측 변동(predictable variation) 과정은:
\[ \langle \hat{\mathbf{B}} - \mathbf{B} \rangle(t) = \int_0^t \mathbf{X}^-(u) [\text{diag } h(u)] \mathbf{X}^{-\top}(u) \, du \]
여기서 \(h(u)\) 는 각 개체의 위험률 벡터이다. 데이터로 추정하면:
\[ \hat{\Sigma}(t) = \int_0^t \mathbf{X}^-(u) [\text{diag } d\mathbf{N}(u)] \mathbf{X}^{-\top}(u) \]
OLS 의사역행렬을 대입하면 식 (10.2.4)와 일치한다.
7 두 표본 검정 코드 예시
7.1 Step 1: 순수 Python으로 식 (10.2.3) 구현
import numpy as np
# 두 표본 데이터
times_g1 = np.array([3, 5, 7, 12, 15, 18])
cens_g1 = np.array([1, 1, 0, 1, 1, 0])
times_g2 = np.array([2, 6, 8, 10, 14, 20])
cens_g2 = np.array([1, 1, 1, 0, 1, 1])
n1, n2 = len(times_g1), len(times_g2)
n = n1 + n2
# 통합 데이터: (시간, 사건 지시자, 그룹 지시자)
T = np.concatenate([times_g1, times_g2])
D = np.concatenate([cens_g1, cens_g2])
Z = np.concatenate([np.ones(n1), np.zeros(n2)]) # 그룹 1 = 1
# 사건 시점만 정렬
event_times = np.sort(np.unique(T[D == 1]))
# 각 사건 시점에서 식 (10.2.3) 적용
B0_cum, B1_cum = 0.0, 0.0
B0_path, B1_path = [], []
var_B0_cum, var_B1_cum = 0.0, 0.0
for ti in event_times:
# 위험 집합 Y_j(t)
at_risk = (T >= ti).astype(float)
# 설계행렬 X(t): n × 2, 열 = (1, Z)
X = np.column_stack([at_risk, at_risk * Z])
# 사건 발생 벡터 I(t)
I_vec = ((T == ti) & (D == 1)).astype(float)
XtX = X.T @ X
if np.linalg.det(XtX) < 1e-10:
break # tau 도달
XtX_inv = np.linalg.inv(XtX)
delta_B = XtX_inv @ X.T @ I_vec # (2,)
B0_cum += delta_B[0]
B1_cum += delta_B[1]
# 분산 식 (10.2.4)
I_diag = np.diag(I_vec)
V = XtX_inv @ X.T @ I_diag @ X @ XtX_inv.T
var_B0_cum += V[0, 0]
var_B1_cum += V[1, 1]
B0_path.append((ti, B0_cum, var_B0_cum))
B1_path.append((ti, B1_cum, var_B1_cum))
print("시점 | B0(t) | B1(t) | SE(B1)")
for ti, b0, _ in B0_path:
b1 = next(b for t, b, _ in B1_path if t == ti)
sb1 = next(np.sqrt(v) for t, _, v in B1_path if t == ti)
print(f"{ti:4d} | {b0:7.4f} | {b1:7.4f} | {sb1:.4f}")
# B1(tau) 가 양수이면 그룹 1의 누적 초과 위험이 더 크다
print(f"\n최종 B1(tau) = {B1_cum:.4f} ± 1.96 × {np.sqrt(var_B1_cum):.4f}")7.2 Step 2: lifelines로 다중 공변량 Aalen 모형
from lifelines import AalenAdditiveFitter
import pandas as pd
import numpy as np
# 후두암 스타일 합성 데이터
np.random.seed(42)
n = 100
data = pd.DataFrame({
'T': np.random.exponential(2.5, n),
'E': np.random.binomial(1, 0.7, n),
'stage_II': np.random.binomial(1, 0.25, n),
'stage_III': np.random.binomial(1, 0.30, n),
'stage_IV': np.random.binomial(1, 0.20, n),
'age_centered': np.random.normal(0, 10, n),
})
aaf = AalenAdditiveFitter(coef_penalizer=0.5, fit_intercept=True)
aaf.fit(data, duration_col='T', event_col='E')
# 누적 회귀 함수 도표 — 기울기 변화 = 효과의 시간 패턴
ax = aaf.plot(columns=['stage_III', 'stage_IV'])
# 직선이면 Lin-Ying 상수 모형으로 충분
# 곡선이면 Aalen 모형이 추가 정보를 제공
# 평활 hazards (β_k(t))
print(aaf.smoothed_hazards_(bandwidth=0.5).head())8 실무 노트
8.1 음수 추정치 문제
최소제곱 추정은 위험률이 비음이라는 제약을 두지 않는다. 따라서:
- 연속 공변량을 평균 중심화하지 않으면 \(\hat{B}_0(t)\) 가 단조 증가하지 않을 수 있다
- “기저 위험”의 의미가 모호해진다
해결책: 모든 연속 공변량을 평균 중심화한다. 그러면 \(\beta_0(t)\) 는 “평균적 개체”의 위험률이 되어 음수 가능성이 크게 줄어든다.
8.2 Aalen의 특이한 성질
Aalen(1989)이 증명한 흥미로운 성질:
공변량 독립성 하 축소 불변성: 만약 한 공변량 \(Z_k\) 가 나머지 공변량들과 통계적으로 독립이면, 모형에서 \(Z_k\) 를 제거해도 나머지 회귀 함수 \(\beta_l(t)\) (\(l \neq k\))는 동일하게 추정된다. 변하는 것은 \(\beta_0(t)\) 만이다.
이는 Cox 모형에서는 성립하지 않는 성질이다. Cox에서는 어떤 공변량을 빼도 다른 공변량의 계수가 변할 수 있다 (omitted variable bias). 가산 모형의 이 성질은 회귀 결과의 해석을 더 간단하게 만든다.
8.3 측정 오차 강건성
또 다른 성질: \(Z_k\) 가 가법적 정규 측정 오차를 갖는다면 (즉 관측값 = 참값 + \(\epsilon_k\), \(\epsilon_k \sim N(0, \sigma^2)\)):
- 가산 모형: 회귀 함수 \(\beta_k(t)\) 가 일정한 비율로 축소되지만 모형 형태는 보존됨
- Cox 모형: 모형 형태 자체가 변형됨 (지수의 로그 정규 가중)
따라서 측정 오차가 우려되는 상황에서 가산 모형이 해석적으로 더 안정적이다.
8.4 생존 함수 복원
추정된 \(\hat{B}_k(t)\) 들을 사용하면 특정 공변량 값 \(\mathbf{Z}\) 를 가진 개체의 누적 위험과 생존 함수를 추정할 수 있다:
\[ \hat{H}(t \mid \mathbf{Z}) = \hat{B}_0(t) + \sum_{k=1}^p \hat{B}_k(t) Z_k \]
\[ \hat{S}(t \mid \mathbf{Z}) = \exp[-\hat{H}(t \mid \mathbf{Z})] \]
또는 product-limit 형태:
\[ \hat{S}^*(t \mid \mathbf{Z}) = \prod_{T_j \leq t} [1 - \Delta \hat{H}(T_j \mid \mathbf{Z})] \]
두 추정량은 점근적으로 동등하지만, \(\hat{H}(t \mid \mathbf{Z})\) 가 단조 증가하지 않을 수 있으므로 해석에 주의가 필요하다.
9 핵심 요약
Aalen 가산 위험 모형은 위험률을 공변량의 시간 변동 선형 결합으로 모형화하고, 누적 회귀 함수 \(B_k(t) = \int_0^t \beta_k(u) du\) 를 시점별 OLS의 누적으로 추정한다. \(\hat{B}_k(t)\) 도표의 기울기 변화는 효과의 시간적 동태를 직접 시각화하며, 두 표본 문제에서 추정량은 두 Nelson-Aalen 추정량의 차이로 환원된다.
| 핵심 식 | 역할 |
|---|---|
| 식 (10.2.1) | 가산 위험 모형 정의 |
| 식 (10.2.2) | 누적 회귀 함수 \(B_k(t)\) |
| 식 (10.2.3) | 최소제곱 추정량 |
| 식 (10.2.4) | 분산 추정 |
| 식 (10.2.6) | 두 표본: NA 차이 |
| 식 (10.2.13) | 가설 검정 통계량 |
| 식 (10.2.14) | 마팅게일 회귀 분해 |
10 관련 주제
선행 지식
후속 주제
관련 개념
- Ch.6 — 커널 평활 — \(\hat{\beta}_k(t)\) 평활 추정
- § 6.3 — 초과 사망 모형 — Aalen의 특수 사례