Kwangmin Kim - Klein § 10.3–10.4 — Lin-Ying Additive Hazards Model & Exercises

1 § 10.3 Lin-Ying 모형의 동기 — 왜 상수 계수인가

1.1 Aalen의 강점과 약점

§ 10.2에서 살펴본 Aalen 모형은 회귀 계수를 시간의 함수 \(\beta_k(t)\) 로 둔다. 이는 유연성의 원천이지만 동시에 부담의 원천이기도 하다.

측면	Aalen \(\beta_k(t)\)	Lin-Ying \(\alpha_k\) (상수)
유연성	효과의 시간 변화 포착 가능	시간 평균 효과만 추정
추정 대상	함수 \(B_k(t) = \int \beta_k(u) du\)	스칼라 \(\alpha_k\)
결과 표현	도표 (시각적)	수치 (HR-like)
추론	가중 카이제곱 함수 검정	단순 t-검정 / 카이제곱
닫힌 형식	시점별 OLS 누적	단일 행렬 연산
효율성	진실이 시간 변동 시 우수	진실이 상수 시 우수

실무에서는 종종 효과가 시간에 거의 일정하다고 가정해도 무방한 경우가 많다. 그렇다면 불필요한 비모수 자유도를 제거하고 더 강력한 추론을 얻는 것이 합리적이다. 이 동기에서 Lin & Ying(1994, 1995, 1997)이 제안한 모형이 다음이다.

1.2 핵심 모형 — 식 (10.3.1)

개체 \(j\) 의 조건부 위험률을 다음과 같이 모형화한다:

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

핵심 차이:

\(\alpha_0(t)\): 임의의 비모수 기저 함수 — 형태에 제약 없음
\(\alpha_k\) (\(k \geq 1\)): 상수 — 시간에 무관

이는 Cox 모형의 가산 버전이라 볼 수 있다:

	Cox PH	Lin-Ying
모형	\(h_0(t) \cdot \exp(\beta^\top Z)\)	\(\alpha_0(t) + \alpha^\top Z\)
기저	\(h_0(t)\) 비모수	\(\alpha_0(t)\) 비모수
효과	상수 \(\beta\) (곱셈)	상수 \(\alpha\) (덧셈)
추정	편우도 (반복)	추정방정식 (닫힌 형식)
분산	Fisher 정보 역행렬	샌드위치

직관적 의미

Lin-Ying에서 \(\alpha_k\) 는 “공변량 \(Z_k\) 가 1단위 증가할 때 위험률이 항상 더해지는 양”이다.

예: \(\alpha_1 = 0.005\) → “흡연(Z=1)은 비흡연(Z=0)에 비해 단위 시간당 0.005만큼의 위험을 추가로 발생시킨다.” 이 추가량은 시점에 무관하다.

Cox의 HR이 “위험이 몇 배인가”라면, Lin-Ying의 \(\alpha\) 는 “위험이 얼마나 더해지는가”를 직접 답한다. 공중보건 의사결정(“연간 추가 사망자 수”)에서 후자가 더 직접적이다.

1.3 실무적 장점: 닫힌 형식

본 절에서는 모든 공변량이 시점 0에 고정된 경우를 다룬다. 이 경우 반복 알고리즘 없이 단일 행렬 공식으로 \(\hat{\alpha}\) 가 얻어진다. 이는 Cox 모형(Newton-Raphson 반복 필요)과 비교했을 때 매우 큰 계산상 이점이다. 시간의존 공변량을 다루는 일반 형태는 Lin & Ying(1994, 1997)을 참조한다.

2 추정량 구성 — 세 행렬 A, B, C

Lin-Ying 추정량은 세 가지 핵심 행렬로부터 도출된다. 각 행렬의 역할을 직관적으로 이해하면 식 자체가 자연스럽게 보인다.

2.1 위험 집합 평균 — 식 (10.3.2)

먼저 시점 \(t\) 에서 위험에 노출된 개체들의 공변량 평균을 정의한다:

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

분자는 위험 집합의 공변량 합, 분모는 위험 집합 크기 \(Y_\cdot(t)\). 직관: “이 시점에 살아 있는 사람들의 평균 공변량값”. 이는 Cox 모형의 부분우도에서도 등장하는 핵심 양으로, 위험 집합 내 공변량 분포를 요약한다.

2.2 행렬 A — 공변량의 시간 가중 분산 (식 10.3.3)

\[ \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)] \tag{식 10.3.3} \]

여기서 시점은 \(0 = T_0 \leq T_1 \leq T_2 \leq \cdots \leq T_n\) 으로 정렬되어 있고, \((T_j - T_{j-1})\) 는 인접 사건 사이의 간격이다.

A 행렬의 직관

내부의 \([Z_i - \bar{Z}(T_j)]^\top [Z_i - \bar{Z}(T_j)]\) 는 시점 \(T_j\) 에서 개체 \(i\) 의 공변량이 위험 집합 평균으로부터 얼마나 떨어졌는지를 나타내는 \(p \times p\) 행렬이다. 이를:

모든 위험 집합 시점 \(T_j\) 에 대해 (시간 가중 \(T_j - T_{j-1}\) 곱하기)
모든 개체 \(i\) 에 대해

합산한 것이 \(\mathbf{A}\) 이다. 의미적으로:

통상 OLS의 \(\sum (X_i - \bar{X})^2\) (총제곱합) 의 생존 분석 버전
“공변량의 시간 가중 분산” 또는 정보 행렬(Fisher information의 가산 모형 대응물)

\(\mathbf{A}\) 가 클수록 공변량이 풍부한 정보를 제공한다 → 추정의 효율이 높다.

2.3 벡터 B — 사건 발생자의 공변량 편차 (식 10.3.4)

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

이는 사건이 발생한 개체에 한해서(따라서 \(\delta_i = 1\) 인 항만 살아남음) 공변량의 위험집합 평균 대비 편차를 합한 것이다.

B 벡터의 직관

만약 공변량 효과가 없다면 (\(H_0: \alpha = 0\)), 사건은 위험 집합 내에서 무작위로 발생하므로 사건 발생자의 평균 공변량은 위험 집합 평균과 같아야 한다. 즉:

\[ E\left[\sum_i \delta_i (Z_i - \bar{Z}(T_i))\right] = 0 \quad \text{under } H_0 \]

따라서 \(\mathbf{B}\) 의 절댓값이 클수록 “사건 발생자의 공변량이 비사건자보다 체계적으로 다르다” → 공변량 효과의 증거. 이는 스코어 함수의 역할을 한다.

비교: Cox 모형의 스코어 함수 \(\sum_i \delta_i (Z_i - \tilde{Z}(T_i))\) 도 정확히 같은 구조를 갖는다. 단지 가중 평균 \(\tilde{Z}\) 가 위험률 가중치를 사용하는 점만 다르다.

2.4 행렬 C — 사건 발생자의 공변량 편차 제곱합 (식 10.3.5)

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

\(\mathbf{B}\) 와 같은 합산 범위(사건 발생자만)이지만, 편차의 제곱(외적)을 합한다. 이는 \(\mathbf{B}\) 의 표본 분산을 추정하는 역할이다 → 로버스트 분산 구성 요소.

2.5 추정량과 분산 — 식 (10.3.6), (10.3.7)

이제 본 추정량:

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

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

식의 구조 — 회귀와 샌드위치

식 (10.3.6)을 OLS 정규방정식 \(\hat{\beta} = (X^\top X)^{-1} X^\top y\) 와 비교하라:

가산 모형	OLS 대응
\(\mathbf{A}\)	\(X^\top X\) (정보)
\(\mathbf{B}\)	\(X^\top y\) (스코어)
\(\hat{\alpha} = A^{-1} B\)	\(\hat{\beta} = (X^\top X)^{-1} X^\top y\)

식 (10.3.7)은 샌드위치 분산(sandwich variance) 형태이다:

\[ A^{-1} C A^{-1} \quad \leftrightarrow \quad (X^\top X)^{-1} \cdot \widehat{\text{cov}}(\text{score}) \cdot (X^\top X)^{-1} \]

이는 모형이 정확히 맞지 않더라도 작동하는 로버스트 분산이다. 모형이 정확히 맞을 때 \(\mathbf{C} \approx \mathbf{A}\) 가 되어 분산이 \(\mathbf{A}^{-1}\) 로 단순화된다 (Fisher 정보의 역행렬).

3 가설 검정

3.1 단일 모수 — 식 (10.3.8)

\(H_j: \alpha_j = 0\) 을 검정하려면 표준 z-통계량을 사용한다:

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

3.2 동시 검정 — 식 (10.3.9)

부분집합 \(J \subset \{1, \ldots, p\}\) 에 대해 \(\alpha_J = 0\) 을 동시 검정:

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

이는 다변량 Wald 검정의 표준 형태이다. Cox 모형의 Wald 검정과 구조적으로 동일하다.

4 두 표본 문제: 명시적 풀이

공변량이 단일 이항 변수 \(Z = 0/1\) 인 경우 Lin-Ying 식이 매우 간단해진다.

4.1 \(\bar{Z}(t)\) 의 의미

\[ \bar{Z}(t) = \frac{N_1(t)}{N_1(t) + N_2(t)} \]

여기서 \(N_k(t)\) 는 시점 \(t\) 의 그룹 \(k\) 위험 집합 크기. 이는 “시점 \(t\) 에 살아 있는 사람이 그룹 1에 속할 확률”이다. 두 그룹 크기가 같으면 0.5, 그룹 1이 빨리 죽으면 점차 작아진다.

4.2 유방암 예제 (Example 10.1)

면역과산화효소(immunoperoxidase) 음성 vs 양성 유방암 환자 비교:

모형: \(h(t \mid Z) = \alpha_0(t) + \alpha Z\) (\(Z = 1\) if 음성)
결과: \(\hat{\alpha} = 0.00803\), SE = 0.00471
검정: \(Z = 1.704\), p = 0.088

해석

\(\hat{\alpha} = 0.00803\) 은 “음성 환자가 양성 환자에 비해 단위 시간당 약 0.008건의 추가 사망 위험을 갖는다”는 뜻이다 (시간 단위에 따라 해석).

p = 0.088 → 5% 수준에서 유의하지 않으나 경계적. Aalen 모형의 결과(p = 0.088)와 거의 동일한 결론이다.

이는 우연이 아니다. § 10.2 마지막에서 보았듯, \(\hat{B}_1(t)\) 도표가 거의 직선이면 효과가 시간에 걸쳐 일정 → 상수 모형이 적합 → 두 모형 결과가 일치한다.

4.3 로그순위와의 비교 (재방문)

같은 데이터에서 로그순위 검정은 \(Z = -2.344\), p = 0.019였다. Lin-Ying은 p = 0.088. 차이의 원인은 § 10.2에서와 동일하다:

로그순위: \(H_0\) 가정 하 분산 추정 (작은 분산)
Lin-Ying: 일반 모형 하 분산 추정 (\(C\) 행렬, 샌드위치)

따라서 로그순위가 더 강력하지만(power), 모형 misspecification에 덜 robust하다.

5 다중 공변량: 후두암 예제 (Example 10.2)

90명 남성 후두암 환자, 4개 공변량:

\(Z_1\): Stage II 지시 변수
\(Z_2\): Stage III 지시 변수
\(Z_3\): Stage IV 지시 변수
\(Z_4\): 진단 시 연령 (평균 64.11세 중심화)

5.1 추정 결과

\[ \hat{\alpha} = (0.01325, \; 0.07654, \; 0.37529, \; 0.00256)^\top \]

분산 행렬 (대각선만):

\[ \widehat{V}_{11} = 1.90 \times 10^{-3}, \quad \widehat{V}_{22} = 1.95 \times 10^{-3}, \quad \widehat{V}_{33} = 1.97 \times 10^{-2}, \quad \widehat{V}_{44} = 4.74 \times 10^{-6} \]

5.2 ANOVA 표

효과	\(\hat{\alpha}\)	SE	\(\chi^2\)	df	p-값
\(Z_1\): Stage II	0.01325	0.0435	0.09	1	0.761
\(Z_2\): Stage III	0.07654	0.0442	3.00	1	0.083
\(Z_3\): Stage IV	0.37529	0.1403	7.16	1	0.008
\(Z_4\): 연령	0.00256	0.0022	1.38	1	0.240

해석의 강력함

Stage IV 환자는 Stage I 대비 연간 약 0.375건의 초과 사망 위험을 갖는다 (시간 단위 = 년 가정). 이는 다음과 같이 직접 사용 가능하다:

100명의 Stage IV 환자에서 연간 약 38명 추가 사망 예상
Stage I 환자의 기저 사망률 위에 단순히 더해서 계산

Cox 모형의 HR이 “위험이 몇 배”라고 알려준다면, Lin-Ying의 \(\hat{\alpha}\) 는 “위험이 얼마나 더해지는가”를 직접 알려준다. 공중보건 자원 배분이나 신뢰구간 기반 의사결정에서 후자가 압도적으로 직관적이다.

5.3 동시 검정: 병기 효과

병기 효과 (\(\alpha_1 = \alpha_2 = \alpha_3 = 0\))의 동시 검정:

\[ \chi^2 = \hat{\alpha}_q^\top \widehat{V}_q^{-1} \hat{\alpha}_q = 9.84 \]

자유도 3의 카이제곱 분포 → p = 0.020. 병기 효과는 유의하다.

Aalen 모형의 동시 검정 결과 (\(\chi^2 = 9.18\), p = 0.027)와 매우 가깝다. 두 모형이 거의 같은 결론을 내린다는 것은 효과가 시간에 걸쳐 거의 일정함을 시사한다.

6 기저 위험 함수의 추정

6.1 Practical Note 1 — 식 (10.3.11) 단순화

\(\hat{\alpha}\) 를 구한 후, 누적 기저 위험 \(A_0(t) = \int_0^t \alpha_0(u) du\) 을 추정할 수 있다. 모든 공변량이 고정 시점이면:

\[ \hat{A}_0(t) = \underbrace{\sum_{T_i \leq t} \frac{\delta_i}{Y_\cdot(T_i)}}_{\text{Nelson-Aalen}} - \underbrace{\sum_{j=1}^{k} \sum_{i=j}^{n} \hat{\alpha}^\top Z_i \frac{T_j - T_{j-1}}{Y_\cdot(T_j)} - \sum_{i=k+1}^{n} \hat{\alpha}^\top Z_i \frac{t - T_k}{Y_\cdot(T_k)}}_{\text{공변량 보정}} \]

(여기서 \(T_k \leq t \leq T_{k+1}\), \(Y_\cdot(t)\) = 총 위험 집합 크기)

직관

첫째 항은 공변량 무시 시 통상의 Nelson-Aalen 추정량이다.

둘째 항은 공변량이 위험에 기여한 부분을 빼는 보정이다. 모든 개체의 공변량이 \(\hat{\alpha}\) 만큼 위험률에 기여하므로, 그것을 누적해서 차감해야 순수한 기저 위험이 남는다.

이는 계단 함수가 아니라 piecewise linear 함수이다. Nelson-Aalen 부분은 사건 시점에서 점프하지만, 보정 항은 시간 간격에 걸쳐 연속적으로 누적된다.

7 계수과정 기반 유도 — 왜 식 (10.3.6)인가

7.1 강도 모형

각 개체의 계수과정 \(N_i(t)\) 를 마팅게일로 분해:

\[ N_i(t) = M_i(t) + \int_0^t Y_i(u) [\alpha_0(u) + \alpha^\top Z_i(u)] du \tag{식 10.3.10} \]

\(\alpha\) 가 안다고 가정하면, Nelson-Aalen 유사 추정량으로 \(A_0(t)\) 를 추정:

\[ \hat{A}_0(t) = \int_0^t \frac{\sum_i \{dN_i(u) - Y_i(u) \alpha^\top Z_i(u) du\}}{\sum_i Y_i(u)} \tag{식 10.3.11} \]

분자는 “관찰된 사건 - 공변량으로 설명되는 사건”, 분모는 위험 집합 크기. 즉, 공변량 효과를 빼고 남은 위험을 평균낸 것이다.

7.2 스코어 방정식의 모방

Cox 모형의 부분 스코어 방정식 (Theoretical Note 2 in § 8.2):

\[ \mathbf{U}(\beta) = \sum_i \int_0^\infty Z_i(t) [dN_i(t) - Y_i(t) \exp(\beta^\top Z_i(t)) d\hat{\Lambda}_0(t, \beta)] \]

Lin-Ying은 이 형태를 가산 모형에 대해 그대로 모방한다:

\[ \mathbf{U}(\alpha) = \sum_i \int_0^\infty Z_i(t) [dN_i(t) - Y_i(t) d\hat{A}_0(t) - Y_i(t) \alpha^\top Z_i(t) dt] \tag{식 10.3.12} \]

\(\hat{A}_0(t)\) 를 식 (10.3.11)로 대입하면:

\[ \mathbf{U}(\alpha) = \sum_i \int_0^\infty [Z_i(t) - \bar{Z}(t)] [dN_i(t) - Y_i(t) \alpha^\top Z_i(t) dt] \tag{식 10.3.13} \]

핵심 — 평균으로의 환원

식 (10.3.12) → (10.3.13) 의 변환에서 결정적인 단계는 \(\hat{A}_0\) 의 대입을 통해 \(Z_i\) 가 위험집합 평균 \(\bar{Z}\) 로부터의 편차 \(Z_i - \bar{Z}\) 로 바뀌는 것이다.

이 형태가 자연스러운 이유는 명확하다 — 공변량 효과가 없다면 \(\bar{Z}\) 는 위험집합 평균이고 \(dN_i\) 는 위험집합 내 무작위로 분포하므로 식 (10.3.13)의 기댓값이 0이 된다.

\(\hat{\alpha}\) 는 이 식을 0으로 만드는 값으로 정의되며, 모든 공변량이 고정 시점인 경우 식 (10.3.6)의 닫힌 형식을 얻는다.

8 § 10.4 연습문제 풀이

§ 10.4에는 세 개의 연습문제(10.1, 10.2, 10.3)가 있다. 모두 가산 위험 모형의 실전 적용을 다룬다.

8.1 Exercise 10.1 — 림프종 이식 (Aalen)

§ 1.10의 림프종 환자 데이터: Allogenic vs Autologous 이식 (호지킨 림프종 vs 비호지킨 림프종).

(a) 이식 유형 효과만: 단일 공변량 모형 \(h(t \mid Z) = \beta_0(t) + \beta_1(t) Z\) (\(Z = 1\) if Auto). \(\hat{B}_1(t)\) 도표를 그려 누적 초과 위험의 시간 패턴을 본다.

(b) 검정 (가중치 = 위험 집합 크기): 식 (10.2.11)-(10.2.13) 적용. \(\chi^2\) 통계량으로 \(H_0: \beta_1(t) \equiv 0\) 검정.

(c) 교호작용 모형: \(Z_1\) = 질환 (NHL=1), \(Z_2\) = 이식 (Allo=1), \(Z_3 = Z_1 \times Z_2\). 교호작용 검정은 \(H_0: \beta_3(t) \equiv 0\) 의 카이제곱 검정으로 수행한다.

풀이 전략

이식 유형이 질환에 따라 다른 효과를 갖는지가 핵심 질문이다. \(Z_3\) 의 누적 회귀 함수 \(\hat{B}_3(t)\) 가:

평평한 0 근방 → 교호작용 없음
단조 증가/감소 → 시간 일정 교호작용
곡률 → 시간 변동 교호작용

Aalen 모형은 이 시간 패턴을 직접 시각화하므로 NHL 환자에서만 Allo가 우수한지, 아니면 효과가 시간에 따라 변하는지 등을 도표에서 읽을 수 있다.

8.2 Exercise 10.2 — 화상 환자 포도상구균 감염

§ 1.6의 화상 환자 데이터: 두 가지 청결 방법(routine vs chlorhexidine)과 화상 면적의 영향.

(a) Aalen 적합: 두 공변량 모형 — \(Z_1\) = chlorhexidine 지시, \(Z_2\) = 면적 (평균 중심화). \(\hat{B}_0(t), \hat{B}_1(t), \hat{B}_2(t)\) 도표 작성.

(b) 생존 함수 추정: 면적 50%인 개체에 대해 두 청결 방법별 \(\hat{S}(t \mid Z)\) 를 그린다:

\[ \hat{S}(t \mid Z) = \exp\left[-\hat{B}_0(t) - \hat{B}_1(t) Z_1 - \hat{B}_2(t) Z_2\right] \]

이 식은 Practical Note 6의 누적 위험 보정 공식에서 직접 도출된다.

(c) 검정: \(H_0: \beta_1(t) \equiv 0\), 가중치 = 위험 집합 크기.

(d) 커널 평활: biweight 커널, 대역폭 10일로 \(\hat{\beta}_1(t)\) 평활 추정. 도표가:

일정 음수 → chlorhexidine이 일관되게 위험을 낮춤
초기 큰 음수 → 평수 → 초기에만 효과
음→양 변환 → 시간에 따라 효과가 역전 (이상 패턴)

이런 시간 패턴 정보는 Cox 모형으로는 절대 얻을 수 없다.

8.3 Exercise 10.3 — 림프종 이식 (Lin-Ying)

10.1과 같은 데이터를 Lin-Ying 모형으로 분석.

(a) 이식 유형 효과: \(\hat{\alpha} = A^{-1} B\) 계산. 닫힌 형식이므로 한 번의 행렬 연산.

(b) 교호작용 검정: \(\chi^2 = \hat{\alpha}_3^2 / \widehat{V}_{33}\) 으로 \(H_0: \alpha_3 = 0\) 검정.

Aalen vs Lin-Ying 비교

10.1과 10.3은 같은 데이터·같은 모형 구조에 두 가지 추정 방식을 적용하라는 의도이다. 결과를 비교했을 때:

결과가 비슷 → 효과가 시간에 거의 일정 → Lin-Ying이 더 효율적
결과가 다름 → 효과가 시간에 따라 변동 → Aalen 도표가 추가 정보를 줌

이 비교 자체가 모형 진단의 도구로 활용된다.

9 코드 예시

9.1 Step 1: 순수 Python으로 식 (10.3.6) 구현

import numpy as np

# 작은 데이터 예시
np.random.seed(42)
n = 50
T = np.random.exponential(2.0, n)
delta = np.random.binomial(1, 0.7, n)
Z = np.column_stack([
    np.random.binomial(1, 0.4, n),  # 공변량 1
    np.random.normal(0, 1, n),      # 공변량 2
])
p = Z.shape[1]

# 시간 정렬
order = np.argsort(T)
T, delta, Z = T[order], delta[order], Z[order]

# 식 (10.3.2): 위험집합 평균 Z_bar(T_j)
def Z_bar_at(t):
    Y = (T >= t).astype(float)
    if Y.sum() == 0:
        return np.zeros(p)
    return (Z * Y[:, None]).sum(axis=0) / Y.sum()

# 식 (10.3.3): A 행렬 — 공변량의 시간 가중 분산
A = np.zeros((p, p))
T_full = np.concatenate(([0.0], T))
for j in range(1, n + 1):
    dt = T_full[j] - T_full[j - 1]
    Zb = Z_bar_at(T_full[j])
    for i in range(j - 1, n):  # i >= j-1 이면 위험집합에 있음
        if T[i] >= T_full[j]:
            diff = (Z[i] - Zb).reshape(-1, 1)
            A += dt * (diff @ diff.T)

# 식 (10.3.4): B 벡터 — 사건 발생자의 공변량 편차 합
B = np.zeros(p)
for i in range(n):
    if delta[i] == 1:
        B += Z[i] - Z_bar_at(T[i])

# 식 (10.3.5): C 행렬 — 사건 발생자의 편차 외적 합
C = np.zeros((p, p))
for i in range(n):
    if delta[i] == 1:
        diff = (Z[i] - Z_bar_at(T[i])).reshape(-1, 1)
        C += diff @ diff.T

# 식 (10.3.6), (10.3.7): 추정량 + 샌드위치 분산
A_inv = np.linalg.inv(A)
alpha_hat = A_inv @ B
V_hat = A_inv @ C @ A_inv

# z-검정 (식 10.3.8)
SE = np.sqrt(np.diag(V_hat))
z_stat = alpha_hat / SE
p_values = 2 * (1 - 0.5 * (1 + np.sign(z_stat) * (1 - 2 * np.exp(-0.717 * np.abs(z_stat) - 0.416 * z_stat**2))))

print("Lin-Ying 가산 모형 결과:")
for k in range(p):
    print(f"  α_{k+1} = {alpha_hat[k]:7.4f}  SE = {SE[k]:.4f}  z = {z_stat[k]:6.3f}")

9.2 Step 2: lifelines 비교 — Aalen vs Lin-Ying

import pandas as pd
import numpy as np
from lifelines import AalenAdditiveFitter

# 합성 데이터
np.random.seed(0)
n = 200
df = pd.DataFrame({
    'T': np.random.exponential(2.5, n),
    'E': np.random.binomial(1, 0.7, n),
    'treatment': np.random.binomial(1, 0.5, n),
    'age': np.random.normal(0, 5, n),
})

# Aalen 모형 (시간 변동 계수)
aaf = AalenAdditiveFitter(coef_penalizer=0.5, fit_intercept=True)
aaf.fit(df, duration_col='T', event_col='E')

# 누적 회귀 함수 도표
ax = aaf.plot(columns=['treatment', 'age'])
# 직선 패턴이면 Lin-Ying 상수 모형이 적합

# Lin-Ying은 lifelines에 직접 구현되어 있지 않으나
# coef_penalizer를 충분히 높여 상수에 가깝게 만들 수 있음
# 또는 R의 timereg::aalen 함수에서 const() 옵션 사용

# 진단: 누적 회귀 함수의 직선성
import matplotlib.pyplot as plt
B_treatment = aaf.cumulative_hazards_['treatment']
times = B_treatment.index

# 선형 회귀 적합도
from numpy.polynomial import polynomial as P
slope, intercept = np.polyfit(times, B_treatment.values, 1)
predicted = slope * times + intercept
r_squared = 1 - np.sum((B_treatment.values - predicted)**2) / np.sum((B_treatment.values - B_treatment.mean())**2)

print(f"treatment의 B(t) 직선 적합도 R² = {r_squared:.3f}")
# R²이 1에 가까우면 Lin-Ying 모형이 적합

9.3 Step 3: R로 Lin-Ying 직접 적합

# R: timereg 패키지의 aalen() 함수에서 const() 사용
library(timereg)

# data: T, E, treatment, age 컬럼
fit <- aalen(Surv(T, E) ~ const(treatment) + const(age), data = df)
summary(fit)
# const() 안의 변수는 시간 일정 효과 (Lin-Ying 스타일)
# const() 밖의 변수는 시간 변동 효과 (Aalen 스타일)

# 상수 가정 검정
test_stationary <- summary(fit)$pval.testBeq0
# H_0: 효과가 시간 일정인가?

10 핵심 요약

§ 10.3-10.4 한 줄 요약

Lin-Ying 가산 위험 모형은 Aalen 모형의 시간 변동 계수를 상수로 제한하여, 세 행렬 \(\mathbf{A}\) (시간 가중 정보), \(\mathbf{B}\) (스코어), \(\mathbf{C}\) (로버스트 분산 구성)로부터 닫힌 형식 추정량 \(\hat{\alpha} = \mathbf{A}^{-1} \mathbf{B}^\top\) 와 샌드위치 분산 \(\widehat{V} = \mathbf{A}^{-1} \mathbf{C} \mathbf{A}^{-1}\) 를 얻는다. 효과가 시간에 거의 일정한 상황에서 Aalen보다 효율적이며, Cox 모형과 같은 단일 추정 + Wald 검정 워크플로를 가산 척도에서 제공한다.

핵심 식	역할
식 (10.3.1)	Lin-Ying 가산 모형 정의
식 (10.3.2)	위험집합 공변량 평균 \(\bar{Z}(t)\)
식 (10.3.3)	\(\mathbf{A}\) — 정보 행렬 대응물
식 (10.3.4)	\(\mathbf{B}\) — 스코어 벡터
식 (10.3.5)	\(\mathbf{C}\) — 로버스트 분산 핵
식 (10.3.6)	\(\hat{\alpha} = A^{-1} B^\top\)
식 (10.3.7)	샌드위치 분산
식 (10.3.13)	평균 편차 스코어 방정식

11 관련 주제

선행 지식

후속 주제

Ch.11 — 회귀 진단

관련 개념

샌드위치 분산 추정량 — 모형 robust 분산
현대 ML 생존 분석 확장