Klein Ch.13 — Multivariate Survival Analysis

1 도입 — 독립성 가정이 깨지는 순간

Klein 교재의 Ch.4 (Kaplan-Meier) 부터 Ch.12 (모수 회귀) 까지 모든 추정·검정은 한 가지 공통 가정을 깔고 있었다: 개체들의 생존시간이 서로 독립이다. 이 가정 덕분에 부분우도 (Cox), Nelson-Aalen, MLE 모두 합곱 형태로 표현되어 깔끔한 점근 이론을 얻었다.

그런데 다음 상황을 보자.

• Litter (한 배 새끼) 효과: 같은 어미에서 태어난 새끼 쥐들은 유전·자궁환경·초기 영양을 공유한다. 한 형제의 종양 발생이 다른 형제의 위험과 무관하다고 가정할 근거가 없다.
• 부부 코호트: 같은 가정에서 살아온 부부는 식습관·흡연·환경 노출을 공유한다. 한 배우자의 사망 위험과 다른 배우자의 위험은 공통 요인을 통해 연관된다.
• 재발 사건: 한 환자에서 첫 번째 심근경색 시점과 두 번째 심근경색 시점은 같은 사람의 사건이므로 강한 연관이 있다.
• 다발성 종양: 한 환자의 좌우 신장 종양 발생 시점은 같은 면역계·유전 배경을 공유한다.

이 모든 경우에 그룹 내 (within-group) 연관이 존재한다. 이를 무시하고 독립 모형을 적용하면:

• 분산 추정이 잘못된다 (보통 과소추정): 정보 행렬 기반 SE 가 실제 변동성보다 작아 신뢰구간이 좁아지고 검정이 과민해진다.
• 점추정도 편향될 수 있다 (frailty 의 경우): 관측되지 않은 그룹 효과를 무시하면 회귀 계수가 0 쪽으로 끌려간다 (attenuation).

직관 — 독립 가정이 깨지면 정보가 부풀려진다

100 마리 쥐가 50 litter 에 두 마리씩 들어 있다고 하자. 형제가 강하게 연관되어 있으면, 두 마리의 정보는 한 마리분에 가깝다 (둘은 거의 같은 운명을 따른다). 그런데 독립 가정 모형은 두 마리를 완전히 다른 정보로 세어 표본 크기를 100 으로 사용한다. 이는 유효 표본 크기 (effective sample size) 를 과대평가한 것이고, 그 결과 SE 가 너무 작게 나온다.

극단 사례: 50 쌍이 모두 일란성 쌍둥이라면 실질 표본은 50 이지 100 이 아니다. 분산이 두 배쯤 부풀려져야 한다 — 이것이 본 장에서 다루는 “분산 보정” 의 본질이다.

Ch.13 은 이 문제에 세 가지 접근을 제시한다.

접근	절	핵심 아이디어
Score test (Commenges-Andersen 1995)	§ 13.2	연관 자체의 존재를 분포 가정 없이 검정 (\(H_0: \sigma = 0\))
Frailty 모형 (Clayton 1978, Hougaard 1986)	§ 13.3	보이지 않는 random effect \(u_i\) 로 그룹 내 연관을 명시적 모형화
Marginal 모형 (Lee et al. 1992)	§ 13.4	각 개체에 marginal PH 만 가정, 분산은 샌드위치로 보정

세 접근은 목적이 다르다. 검정만 필요하면 Score test, 연관 강도까지 추정하려면 Frailty, 회귀 계수와 그 robust SE 만 원하면 Marginal. 같은 데이터 (Mantel litter rat 50 그룹) 에 세 접근을 모두 적용해 결과를 비교하는 것이 본 장의 흐름이다.

2 § 13.1 Frailty 모형 — 보이지 않는 그룹 효과

2.1 모형의 기본 형태

Klein 은 비례위험 모형의 공유 frailty (shared frailty) 확장으로 출발한다. \(i\) 번째 그룹의 \(j\) 번째 개체의 위험률을, 그룹 공유 random effect \(w_i\) 와 함께 다음과 같이 모형화한다:

\[ h_{ij}(t) = h_0(t) \exp(\sigma w_i + \beta^t Z_{ij}), \quad i = 1, \ldots, G, \quad j = 1, \ldots, n_i . \tag{식 13.1.1} \]

여기서:

• \(h_0(t)\): 임의의 베이스라인 위험률
• \(Z_{ij}\): 개체 공변량
• \(w_i\): frailty — 그룹 \(i\) 가 공유하는 보이지 않는 random effect (평균 0, 분산 1)
• \(\sigma\): frailty 의 척도 — \(\sigma = 0\) 이면 일반 Cox 모형

다른 표현으로 \(u_i = \exp(\sigma w_i)\) 를 써서:

\[ h_{ij}(t) = h_0(t) \, u_i \exp(\beta^t Z_{ij}), \quad u_i \sim \text{i.i.d. } G(\cdot), \quad E[u_i] = 1 . \tag{식 13.1.2} \]

이 표현이 더 직관적이다. \(u_i\) 는 그룹별 위험률 곱셈 인자이다.

직관 — frailty \(u_i\) 가 의미하는 것

\(u_i\) 는 그룹 \(i\) 의 “체질적 취약성” 이라고 보면 된다.

• \(u_i = 1\) (평균): 평균적인 그룹. 모든 개체가 베이스라인 위험률을 따른다.
• \(u_i > 1\) (frail group): 위험률이 부풀려진다. 예: \(u_i = 2\) 이면 이 그룹의 모든 형제는 평균 그룹의 두 배 빠르게 사건 경험.
• \(u_i < 1\) (robust group): 위험률이 줄어든다. 이 그룹은 평균보다 오래 산다.

핵심은 \(u_i\) 를 우리가 관측할 수 없다는 점이다. 데이터에서 우리가 보는 것은 \(T_{ij}, \delta_{ij}, Z_{ij}\) 뿐이고, \(u_i\) 는 latent variable 로 적분해서 없애야 한다.

또한 “frail 한 사람이 먼저 죽는다” 라는 자연스러운 결과가 따라온다. 시간이 지날수록 살아남은 사람들은 robust group 출신일 가능성이 높아진다 — 이것이 frailty 모형이 만들어내는 선택 효과 (selection effect) 이며, 인구 단위 위험률이 시간에 따라 떨어지는 패턴을 설명한다.

2.2 Joint Survival Function 과 Laplace 변환

\(u_i\) 가 관측되지 않으므로 그룹 내 개체들의 결합 (joint) 생존함수 는 \(u_i\) 에 대해 적분된 형태로 주어진다. \(u_i\) 의 분포 함수를 \(G\) 라 하면:

\[ S(x_{i1}, \ldots, x_{in_i}) = E_U\!\left[ \exp\!\left(-U \sum_{j=1}^{n_i} H_0(x_{ij}) e^{\beta^t Z_{ij}}\right)\right] . \]

이를 Laplace 변환 \(\mathcal{L}(v) = E[e^{-Uv}]\) 으로 표현하면:

\[ S(x_{i1}, \ldots, x_{in_i}) = \mathcal{L}\!\left[ \sum_{j=1}^{n_i} H_0(x_{ij}) e^{\beta^t Z_{ij}} \right] . \tag{식 13.1.3} \]

왜 Laplace 변환이 등장하는가

\(u_i\) 가 주어졌을 때 그룹 내 개체들은 conditional independent. 따라서:

\[ S(x_{i1}, \ldots \mid u_i) = \prod_j \exp[-u_i H_0(x_{ij}) e^{\beta^t Z_{ij}}] = \exp\!\left[-u_i \sum_j H_0(x_{ij}) e^{\beta^t Z_{ij}}\right] . \]

여기서 \(u_i\) 에 대한 기댓값을 취하면 \(E[e^{-u_i v}]\) 형태가 나오는데, 이것이 정확히 \(u_i\) 의 Laplace 변환이다. 따라서 frailty 분포의 Laplace 변환이 곧 marginal joint survival 을 결정한다.

이 사실 덕분에 frailty 분포 선택 = Laplace 변환 선택 = joint survival 형태 선택이라는 1대1 대응이 성립한다. 분포의 모양이 직접 보이지 않더라도 Laplace 변환의 폐쇄형이 깔끔하면 그것으로 충분하다.

2.3 Frailty 분포 후보 — 어느 것을 쓸 것인가

문헌에 제안된 frailty 분포들은 다음과 같다.

분포	제안자	특징
One-parameter Gamma	Clayton (1978), Clayton & Cuzick (1985)	Laplace 폐쇄형, EM 알고리즘 가능 (§ 13.3)
Positive stable	Hougaard (1986a)	Marginal 도 PH 를 유지하는 유일한 분포
Inverse Gaussian	Hougaard (1986b)	시간에 따라 frailty 가 약해짐
Log-normal	McGilchrist & Aisbett (1991)	linear mixed model 과 호환
Power variance	Aalen (1988, 1992)	Gamma·stable·IG 를 포괄하는 가족
Uniform	Lee & Klein (1988)	유한 지지
Threshold	Lindley & Singpurwalla (1986)	임계값 모형

Gamma 가 가장 많이 쓰이는 이유

세 가지 실용적 이점이 있다.

1. Laplace 폐쇄형: \(\mathcal{L}(v) = (1 + \theta v)^{-1/\theta}\) — 기초 함수로 표현되어 우도가 깔끔하다.
2. Kendall’s \(\tau\) 폐쇄형: 그룹 내 연관도가 \(\tau = \theta / (\theta + 2)\) 로 한 번에 해석된다.
3. EM 알고리즘: \(u_i \mid \text{data}\) 가 다시 Gamma 가 되어 (conjugate) E-step 이 명시적이다.

다만 marginal 분포가 PH 를 따르지 않는다는 단점이 있다. Marginal PH 가 필요하면 positive stable 만 가능하지만, 모수 해석이 어려워 실무에서는 거의 사용되지 않는다.

2.4 Frailty 의 두 가지 용도

식 (13.1.2) 는 두 가지 다른 문제에 같은 형태로 적용된다.

용도 1 — 그룹 내 연관 (다변량 생존): \(G\) 그룹, 각 그룹에 \(n_i \geq 2\) 명. Litter, 부부, 다발 종양 등.

용도 2 — Overdispersion 보정 (단변량 생존): \(G = N\) (개체 수), \(n_i = 1\) (그룹 크기 1). 측정되지 않은 공변량들의 효과를 random effect 로 흡수.

같은 식이 두 문제를 다 푼다

다변량 frailty: “관측되지 않은 그룹 공유 효과가 그룹 내 연관을 만든다.”

단변량 overdispersion: “관측되지 않은 개인별 공변량이 모형의 잔차 분산을 부풀린다.”

두 경우 모두 latent random effect \(u_i\) 를 도입해 보이지 않는 이질성을 흡수하는 것이 핵심이다. 식 (13.1.1) 의 형태는 동일하다.

다만 단변량의 경우 \(\sigma = 0\) 검정의 검정력이 매우 낮다. \(G = N\) 이면 그룹당 정보가 너무 적어 연관을 식별하기 어렵다. 따라서 overdispersion 보정은 본질적으로 모형의 robust 화 의미이지 검정 결과로 결론짓는 분석은 아니다.

3 § 13.2 Score Test — 연관이 존재하는가

3.1 검정의 목표

Frailty 모형의 모든 추정 절차 (§ 13.3) 는 계산이 무겁다. 그러나 연관이 존재하는지 여부만 알고 싶다면 모형 적합 없이 score test 만 수행하면 된다. Commenges & Andersen (1995) 의 검정은 다음 가설을 다룬다:

\[ H_0: \sigma = 0 \quad \text{vs} \quad H_1: \sigma > 0 . \]

핵심 매력: frailty 분포에 대한 어떤 가정도 필요 없다. Gamma 든 IG 든 log-normal 이든, \(H_0\) 하에서는 모두 동일한 score 통계량으로 환원된다.

3.2 검정 통계량의 구성

귀무 모형 (frailty 없는 일반 Cox) 을 적합해 마팅게일 잔차 \(M_{ij} = \delta_{ij} - \widehat{H}_0(T_{ij}) \exp(\widehat{\beta}^t Z_{ij})\) 를 얻는다 (Ch.11.3). 그러면 score 통계량은:

\[ T = \sum_{i=1}^G \left( \sum_{j=1}^{n_i} M_{ij} \right)^2 - D + C , \tag{식 13.2.2} \]

여기서 \(D\) 는 총 사건 수, \(C\) 는 보정 항 (식 13.2.3, \(N \to \infty\) 시 0 으로 수렴).

이를 풀어쓰면:

\[ T = \underbrace{\sum_{i} \sum_{j \neq k} M_{ij} M_{ik}}_{\text{(1) 그룹 내 쌍 상관}} \;+\; \underbrace{\left(\sum_{i,j} M_{ij}^2 - N\right)}_{\text{(2) overdispersion}} \;+\; \underbrace{C}_{\text{(3) 보정}} . \tag{식 13.2.4} \]

직관 — 세 항이 무엇을 잡아내는가

(1) 그룹 내 쌍 상관: 같은 그룹 내 두 개체 \(j, k\) 의 마팅게일 잔차 곱의 합. 두 잔차가 같은 방향 (둘 다 양수 또는 둘 다 음수) 이면 곱이 양수. 그룹 내 잔차가 체계적으로 같이 움직이면 (= 연관 존재) 이 항이 커진다. 주된 신호 항이다.

(2) Overdispersion: 개별 잔차의 제곱합 \(\sum M_{ij}^2\) 가 표본 크기 \(N\) 보다 큰지. PH 모형이 정확하면 \(E[M_{ij}^2] \approx 1\) 이므로 \(\sum M_{ij}^2 \approx N\). 이 항이 양수로 크면 individual variability 가 모형 예상보다 크다는 신호.

(3) 보정 \(C\): 마팅게일 잔차의 추정 오차로 인한 편향을 1차 근사로 보정. \(N \to \infty\) 에서는 무시할 수준.

표준화된 통계량: \(T / \sqrt{V} \xrightarrow{d} N(0, 1)\) 이며, \(V\) 는 식 (13.2.5)–(13.2.9) 의 복잡한 분산 추정량이다.

3.3 예제 13.1 — Mantel et al. (1977) Litter Rat 데이터

50 litter 의 암컷 새끼 쥐를 사용한 종양 유발 연구. 각 litter 에서 1 마리는 약물 투여 (treated), 2 마리는 placebo (control). 관심 사건: 종양 발생 시점.

먼저 일반 Cox 모형 (litter 무시):

• 단일 공변량 \(Z = 1\) (treated), \(Z = 0\) (control)
• \(\hat{\beta} = 0.8975\), \(\widehat{\text{Var}}(\hat{\beta}) = 0.1007\)
• 약물이 종양 발생을 빠르게 함 (HR \(\approx e^{0.9} \approx 2.5\))

Score test 결과:

• \(T = 8.91\)
• \(V = 45.03\)
• \(T/\sqrt{V} = 8.91 / \sqrt{45.03} = 1.33\)
• \(p = 0.184\)

해석

\(p = 0.184 > 0.05\) → 5 % 수준에서 litter 효과의 증거 없음. 같은 어미 새끼들 사이의 연관이 통계적으로 유의하지 않다.

이 결과는 우리에게 약물 효과 (\(\hat{\beta} = 0.90\)) 의 SE 를 그대로 사용해도 큰 문제가 없다고 알려준다. 그러나 검정의 유의 X 는 “연관이 없다” 의 증명이 아니라 “탐지하지 못했다” 일 뿐이다. 표본 크기가 작거나 (50 litter × 3 = 150) frailty 분산이 작으면 검정력 부족으로 통과한다. 따라서 § 13.3, § 13.4 에서 frailty 모형과 marginal 모형도 함께 적합해 보는 것이 권장된다.

4 § 13.3 Gamma Frailty Estimation — 연관 강도 추정

4.1 Gamma Frailty 의 폐쇄형

\(u_i\) 가 평균 1, 분산 \(\theta\) 의 one-parameter Gamma 분포를 따른다고 가정한다:

\[ g(u) = \frac{u^{1/\theta - 1} \exp(-u/\theta)}{\Gamma(1/\theta) \, \theta^{1/\theta}} . \tag{식 13.3.1} \]

이 분포의 Laplace 변환 \(\mathcal{L}(v) = (1 + \theta v)^{-1/\theta}\) 을 식 (13.1.3) 에 대입하면 joint survival 의 폐쇄형이 나온다:

\[ S(x_{i1}, \ldots, x_{in_i}) = \left[ 1 + \theta \sum_{j=1}^{n_i} H_0(x_{ij}) e^{\beta^t Z_{ij}} \right]^{-1/\theta} . \]

직관 — \(\theta\) 의 의미

• \(\theta = 0\): \(u_i \equiv 1\) (분산 0) → 일반 Cox PH 모형, 그룹 내 독립.
• \(\theta\) 가 크다: \(u_i\) 의 분산이 크다 → 그룹별 위험률 차이가 크다 → 그룹 내 연관이 강하다.

연관 강도는 Kendall’s \(\tau\) 로 표현하면 직관적이다:

\[ \tau = \frac{\theta}{\theta + 2} \in [0, 1) . \]

• \(\theta = 0\): \(\tau = 0\) (독립)
• \(\theta = 2\): \(\tau = 0.5\) (중간 강도)
• \(\theta \to \infty\): \(\tau \to 1\) (완전 연관)

Kendall’s \(\tau\) 는 두 개체의 사건 시점이 같은 순서로 일어날 (concordant) 확률이 다른 순서일 (discordant) 확률보다 얼마나 큰가를 측정한다. 따라서 \(\theta\) 자체보다 \(\tau\) 로 보고하면 임상 독자에게 친숙하다.

4.2 관측 우도 — Marginal Likelihood

\(u_i\) 가 관측되지 않으므로 데이터로 우도를 쓰려면 \(u_i\) 를 적분해서 없앤 marginal log-likelihood 를 사용한다:

\[ \begin{aligned} L(\theta, \beta, H_0) = \sum_{i=1}^G \Big\{ &D_i \ln \theta - \ln \Gamma(1/\theta) + \ln \Gamma(1/\theta + D_i) \\ &- (1/\theta + D_i) \ln\!\left[ 1 + \theta \sum_j H_0(T_{ij}) e^{\beta^t Z_{ij}} \right] \\ &+ \sum_j \delta_{ij} \big[ \beta^t Z_{ij} + \ln h_0(T_{ij}) \big] \Big\} , \end{aligned} \tag{식 13.3.2} \]

\(D_i = \sum_j \delta_{ij}\) 는 그룹 \(i\) 의 사건 수.

베이스라인 \(h_0\) 를 모수적 (예: Weibull) 으로 두면 식 (13.3.2) 를 직접 최대화하면 끝. 비모수로 두면 \(H_0\) 가 무한 차원이라 이 그대로는 풀 수 없다 — 여기서 EM 알고리즘이 등장한다.

4.3 EM 알고리즘 — 비모수 베이스라인

핵심 아이디어: \(u_i\) 를 관측할 수 있다고 가정하면 우도가 단순 Cox + offset 형태로 환원된다 (\(L_{\text{FULL}} = L_1(\theta) + L_2(\beta, H_0)\)). 그러나 \(u_i\) 는 latent. EM 은 다음을 반복한다.

E-step — \(u_i\) 를 추측하기

현재 모수 추정값과 데이터를 조건으로 한 \(u_i\) 의 conditional 분포는 다시 Gamma 가 된다 (Gamma 가 conjugate prior 인 덕분):

\[ u_i \mid \text{data} \sim \text{Gamma}(A_i, C_i), \quad A_i = 1/\theta + D_i, \quad C_i = 1/\theta + \sum_j \widehat{H}_0(T_{ij}) e^{\widehat{\beta}^t Z_{ij}} . \]

따라서 \(E[u_i \mid \text{data}] = A_i / C_i\), \(E[\ln u_i \mid \text{data}] = \psi(A_i) - \ln C_i\) (\(\psi\) 는 digamma 함수).

직관: 그룹 \(i\) 의 사건 수 \(D_i\) 가 많으면 (= 그룹이 자주 죽으면) → \(A_i\) 가 커짐 → \(\hat{u}_i\) 추정값 증가. 그룹의 누적 위험 노출 (\(\sum \widehat{H}_0 e^{\widehat{\beta}^t Z}\)) 이 크면 → \(C_i\) 증가 → \(\hat{u}_i\) 감소. “위험에 많이 노출되었는데도 적게 죽으면 robust group, 노출이 적은데 많이 죽으면 frail group” 의 베이지언 갱신이다.

M-step — Cox 부분우도 + frailty offset

E-step 의 \(\hat{u}_i\) 를 known offset 으로 사용한 Cox partial likelihood 를 최대화:

\[ L_3(\beta) = \sum_{k=1}^D \left\{ S_{(k)} - m_{(k)} \ln\!\left[ \sum_{b \in R(t_{(k)})} \hat{u}_b \exp(\beta^t Z_b) \right] \right\} , \tag{식 13.3.7} \]

이는 일반 Breslow partial likelihood 에서 위험 score \(\exp(\beta^t Z_b)\) 자리에 \(\hat{u}_b \exp(\beta^t Z_b)\) 가 들어간 형태. frail group 의 개체는 위험집합 합에 더 큰 가중치로 들어간다.

베이스라인 점프도 update:

\[ \widehat{h}_{k0} = \frac{m_{(k)}}{\sum_{b \in R(t_{(k)})} \hat{u}_b \exp(\hat{\beta}^t Z_b)} . \]

\(\theta\) 는 \(L_4 = E[L_1(\theta) \mid \text{data}]\) 를 최대화해 update.

EM 은 수렴이 느릴 수 있다. Nielsen et al. (1992) 의 profile EM 은 \(\theta\) 를 격자로 잡고 각 \(\theta\) 에 대해 \(\beta_\theta\) 만 EM 으로 구한 뒤 \(L(\theta, \hat{\beta}_\theta)\) 를 비교하는 방식으로 더 빠르다.

4.4 예제 13.1 (계속) — Litter Rat Frailty 적합

같은 50 litter 데이터에 Gamma frailty Cox 모형 적합:

모형	\(\hat{\beta}\) (treatment)	SE(\(\hat{\beta}\))	\(\hat{\theta}\)	SE(\(\hat{\theta}\))
Cox + Frailty	0.904	0.323	0.472	0.462
Cox (독립)	0.897	0.317	—	—

결과 해석

Frailty 검정: \(\hat{\theta} / SE(\hat{\theta}) = 0.472 / 0.462 = 1.02\), \(p \approx 0.31\) → 유의 X. § 13.2 score test (\(p = 0.184\)) 와 일관된 결론.

\(\hat{\beta}\) 의 변화: 독립 모형 (\(\hat{\beta} = 0.897\)) → Frailty 모형 (\(\hat{\beta} = 0.904\)). Frailty 를 고려하면 회귀 계수가 0 에서 멀어졌다 (절댓값 증가). 이는 frailty 모형의 일반적 성질이다 — 관측되지 않은 그룹 효과를 무시하면 회귀 계수가 0 쪽으로 attenuate 되는 경향이 있고, frailty 를 모형화하면 그 효과가 풀려난다.

Kendall’s \(\tau\): \(\hat{\tau} = 0.472 / (0.472 + 2) = 0.19\), SE(\(\hat{\tau}\)) \(\approx 0.43\). 점추정은 약한 양의 연관을 시사하지만 SE 가 커서 유의하지 않음. 표본 크기 (50 litter) 가 frailty 추정에 충분치 않은 한계.

Frailty 분산 추정의 어려움

\(\theta\) 의 SE 가 거의 추정값만 한 크기 (\(\hat{\theta} = 0.47\), SE \(= 0.46\)) 인 점이 흔하다. 이유:

• \(\theta = 0\) 이 모수공간 경계에 있어 표준 점근 이론이 약하게 성립한다 (Wald 검정의 검정력 낮음).
• 그룹 수 \(G\) 가 frailty 추정의 effective sample size 인데, 50 정도는 작다.

따라서 frailty 검정은 score test (§ 13.2) 가 우선이다. EM 적합 후 Wald 검정만으로 결론짓지 말 것.

5 § 13.4 Marginal Model — 분산만 보정하기

5.1 동기 — Frailty 모형의 한 가지 한계

Frailty 모형은 conditional PH 를 가정한다: \(h_{ij}(t \mid u_i) = h_0(t) u_i e^{\beta^t Z_{ij}}\). 그런데 우리가 실제로 관측하는 것은 marginal 분포다. Marginal 위험률은 일반적으로 PH 가 아니다.

예를 들어 Gamma frailty 의 marginal hazard 는:

\[ h_{ij}^{\text{marg}}(t) = \frac{h_0(t) e^{\beta^t Z_{ij}}}{1 + \theta H_0(t) e^{\beta^t Z_{ij}}} , \]

분모에 \(H_0(t)\) 가 들어가 시간에 따라 변하므로 PH 가 아니다. 따라서 frailty 모형에서 추정한 \(\beta\) 는 conditional 효과이며 marginal hazard ratio 와 다르다.

5.2 Lee et al. (1992) 의 접근

Lee, Wei, Amato (1992) 는 다른 길을 택했다: conditional 모형은 만들지 않는다. 각 개체의 marginal hazard 만 PH 로 가정한다.

\[ h_{ij}(t \mid Z_{ij}) = h_0(t) \exp(\beta^t Z_{ij}) . \tag{식 13.4.1} \]

그룹 내 연관은 허용하지만 그 형태는 모형화하지 않는다.

직관 — Frailty vs Marginal 의 핵심 차이

측면	Frailty (§ 13.3)	Marginal (§ 13.4)
PH 가 성립하는 단위	그룹 (conditional on \(u_i\))	개체 (marginal)
회귀 계수 \(\beta\) 의 의미	“같은 그룹 내” 효과	“전체 모집단 평균” 효과
연관 강도 추정	\(\hat{\theta}\), \(\hat{\tau}\) 명시	추정 안 함 (그냥 분산만 부풀림)
추정 복잡도	EM 또는 직접 MLE	일반 Cox 적합 + 샌드위치
분포 가정	Gamma (또는 다른 frailty 분포)	없음

언제 어느 쪽인가: 그룹 효과 자체가 관심사 (예: 가족력의 강도) → frailty. 회귀 계수가 관심사이고 연관은 nuisance → marginal.

5.3 Independence Working Model + 샌드위치 분산

Marginal 모형의 추정 절차는 단순하다.

1. 모든 개체를 독립으로 가정 하고 일반 Cox partial likelihood 를 최대화 → \(\hat{\beta}\).
2. Lee et al. (1992) 가 증명: 이 \(\hat{\beta}\) 는 marginal 모형 (식 13.4.1) 이 정확하면 일치 추정량 (consistent).
3. 그러나 분산 추정은 잘못된다 — 정보 행렬 기반 SE 는 그룹 내 상관을 무시하므로 과소추정.
4. 샌드위치 분산 (sandwich variance) 으로 보정.

샌드위치 분산의 공식:

\[ \widetilde{V} = \widehat{V} \, C \, \widehat{V} , \]

• \(\widehat{V}\): 독립 working model 의 정보 행렬 역수 (= 일반 Cox 의 SE 제곱)
• \(C\): score 잔차의 그룹 내 합의 외적

Score 잔차 (Ch.11.6 식 11.6.1) 를 이용해:

\[ S_{ijk} = \delta_{ij}[Z_{ijk} - \bar{Z}_k(T_{ij})] - \sum_{t_b \leq T_{ij}} [Z_{ijk} - \bar{Z}_k(t_b)] \, d\widehat{H}_0(t_b) , \tag{식 13.4.2} \]

\[ C_{b, k} = \sum_{i=1}^G \sum_{j=1}^{n_i} \sum_{m=1}^{n_i} S_{ijb} S_{imk} . \tag{식 13.4.3} \]

직관 — 샌드위치가 그룹 내 연관을 어떻게 잡아내는가

\(C\) 의 핵심은 \(\sum_{j} \sum_{m} S_{ijb} S_{imk}\) — 같은 그룹 \(i\) 내 두 개체의 score 잔차 곱을 합한다 (\(j \neq m\) 항 포함).

• 그룹 내 score 잔차가 같은 방향 이면 곱이 양수 → \(C\) 증가 → 샌드위치 분산 증가.
• 독립이면 \(E[S_{ij} S_{im}] = 0\) (\(j \neq m\)) → \(C\) 가 일반 정보 행렬에 가까워져 보정 효과 사라짐.

\(\widehat{V} C \widehat{V}\) 형태가 “샌드위치”인 이유: 빵 두 개 (\(\widehat{V}\)) 사이에 속재료 (\(C\)) 가 들어 있는 구조. White (1980) 의 heteroskedasticity-robust SE, GEE, GMM 의 분산 추정량이 모두 같은 형태이다. 모형을 정확히 알지 않아도 작동하는 로버스트 분산 의 표준 패턴이다.

5.4 예제 13.1 (계속) — Litter Rat Marginal 적합

같은 데이터:

• \(\hat{\beta} = 0.897\) (독립 Cox 와 동일)
• \(\widehat{V} = 0.1007\), naive SE \(= \sqrt{0.1007} = 0.3173\)
• \(C = 8.893\)
• \(\widetilde{V} = (0.1007)^2 \times 8.893 = 0.0902\), robust SE \(= \sqrt{0.0902} = 0.3000\)

95 % 신뢰구간 (HR):

분산 추정	95 % CI for HR
Naive (독립 가정)	\(\exp(0.897 \pm 1.96 \times 0.3173) = (1.32, 4.57)\)
샌드위치 (보정)	\(\exp(0.897 \pm 1.96 \times 0.3000) = (1.36, 4.42)\)

흥미로운 결과 — 보정 후 CI 가 더 좁아졌다

이 데이터에서는 샌드위치 분산이 naive 분산보다 작다. 즉 그룹 내 연관이 양이 아니라 약한 음의 방향일 가능성을 시사한다 (또는 우연 변동).

일반적으로 양의 연관이 있으면 샌드위치 분산 > naive 분산 (CI 가 더 넓어짐). 음의 연관이 있으면 (예: litter 내 경쟁 효과 — 한 마리가 약하면 다른 마리가 강하다) 샌드위치 분산이 더 작아질 수 있다. Litter 내 종양 발생에서는 약한 음의 효과가 그럴듯한 시나리오 (한정된 자원의 경쟁 또는 특정 약물 분배 차이).

이 예제는 score test (p=0.184), frailty Wald (p=0.31), marginal CI 비교 모두 “강한 그룹 효과 없음” 으로 일관된 결론을 준다.

5.5 Wei-Lin-Weissfeld 확장

Wei, Lin, Weissfeld (1989, WLW) 는 그룹별로 다른 베이스라인 \(h_{0i}(t)\) 를 허용하는 marginal 모형을 제안했다 (예: 한 환자의 첫 사건 시간과 두 번째 사건 시간이 다른 베이스라인을 따름). 식 (13.4.1) 의 \(h_0(t)\) 가 \(h_{0i}(t)\) 로 바뀌고, 추정·분산은 같은 샌드위치 구조 (Lin 1993).

6 세 접근의 비교 — 같은 데이터, 세 결과

접근	출력	가정	사용 시점
Score test (§ 13.2)	\(T/\sqrt{V}\), p-값	분포 무관, \(H_0\) 만 검정	“연관이 있는가?” 만 알고 싶을 때
Gamma Frailty (§ 13.3)	\(\hat{\beta}\), \(\hat{\theta}\), \(\hat{\tau}\)	Gamma frailty + conditional PH	연관 강도 자체가 관심사
Marginal (§ 13.4)	\(\hat{\beta}\), robust SE	Marginal PH (분포 무관)	회귀 계수만 관심, 연관은 nuisance

실무 권고 흐름

1. 먼저 score test 로 연관 존재 여부 검정. 비계산 부담이 작아 모든 다변량 데이터에 우선 적용.
2. 연관이 유의하지 않으면 일반 Cox 그대로 사용해도 SE 가 크게 잘못되지는 않음 (다만 marginal sandwich 로 한 번 더 확인하면 안전).
3. 연관이 유의하면:
   • 연관 강도가 관심사 → Gamma frailty 적합. Kendall’s \(\tau\) 보고.
   • 연관은 nuisance, 회귀 계수만 → Marginal model 의 robust SE 사용.
4. 민감도 분석: 가능하면 두 접근을 모두 적용해 \(\hat{\beta}\) 와 SE 가 일관되는지 확인. 큰 차이가 나면 모형 가정 (conditional vs marginal PH) 의 적합성 재검토.

7 코드 예시 — Python lifelines / R survival

7.1 Step 1 — Python: Frailty 모형은 lifelines 가 직접 지원 X, 차선 접근

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

# Litter rat 데이터: T, E, treatment, litter_id
df = pd.DataFrame(...)  # 50 litters × 3 rats

# 1. 독립 Cox (베이스라인)
cph = CoxPHFitter()
cph.fit(df, duration_col="T", event_col="E", formula="treatment")
print("Naive SE:", cph.standard_errors_)

# 2. Marginal model — robust 샌드위치 SE
cph_robust = CoxPHFitter()
cph_robust.fit(df, duration_col="T", event_col="E",
               formula="treatment", cluster_col="litter_id")
print("Robust SE:", cph_robust.standard_errors_)
# cluster_col 옵션이 식 (13.4.3) 의 그룹 내 합 외적을 자동 계산

7.2 Step 2 — R: 모든 접근을 직접 지원

library(survival)
library(coxme)

# 1. 독립 Cox
fit_indep <- coxph(Surv(T, E) ~ treatment, data = rats)

# 2. Marginal model — robust SE (식 13.4.3 자동 계산)
fit_marg <- coxph(Surv(T, E) ~ treatment + cluster(litter), data = rats)
summary(fit_marg)
# robust SE 가 출력됨

# 3. Gamma Frailty (Therneau)
fit_frail <- coxph(Surv(T, E) ~ treatment + frailty(litter, dist = "gamma"),
                   data = rats)
summary(fit_frail)
# theta 추정값과 LRT p-value 출력

# 4. coxme 로 log-normal frailty (대안)
fit_lognormal <- coxme(Surv(T, E) ~ treatment + (1 | litter), data = rats)
print(fit_lognormal)

7.3 Step 3 — Score test 직접 구현

# Commenges-Andersen score test 식 (13.2.2)
fit0 <- coxph(Surv(T, E) ~ treatment, data = rats)
M <- residuals(fit0, type = "martingale")
D <- sum(rats$E)
N <- nrow(rats)

# 그룹별 마팅게일 잔차 합 제곱
group_M_sum <- tapply(M, rats$litter, sum)
T_stat <- sum(group_M_sum^2) - D + C  # C 는 식 (13.2.3) 으로 별도 계산

# V 는 식 (13.2.5)–(13.2.9) 로 계산 — 복잡
# 실무에서는 SAS 매크로 또는 ggcoxzph 등 패키지 사용 권장
Z <- T_stat / sqrt(V)
p_value <- 2 * (1 - pnorm(abs(Z)))

8 Ch.13 시리즈

본 장의 세부 절은 후속 포스트에서 다룬다.

• § 13.1-13.2 — Shared Frailty Models & Score Test for Association — frailty 모형 정의, Laplace 변환 유도, Commenges-Andersen score test 의 통계량 분해와 Mantel rat 예제
• § 13.3-13.4 — Gamma Frailty EM Estimation & Marginal Model — EM 알고리즘 상세 (E-step conjugate posterior + M-step Cox + offset), Lee et al. marginal model 과 샌드위치 분산 \(\widetilde{V} = \widehat{V} C \widehat{V}\)
• § 13.5 — 연습문제 풀이 — Skin allograft (Batchelor-Hackett) · Kidney catheter (McGilchrist-Aisbett) 두 데이터에 score test · gamma frailty · marginal model 종합 적용 + 세 데이터 그룹 효과 강도 비교
• Appendix A — Numerical Techniques for Maximization — Cox·Weibull MLE·frailty EM 등 본문의 거의 모든 추정에 등장하는 수치 최적화: Bisection·Secant·Newton-Raphson·Steepest Ascent·Marquardt 의 직관과 trade-off

9 핵심 요약

Ch.13 한 줄 요약

다변량 생존 데이터에서는 shared frailty 모형 (\(h_{ij}(t) = h_0(t) u_i e^{\beta^t Z_{ij}}\), 식 13.1.1) 으로 그룹 내 연관을 명시적으로 모형화하거나, marginal 모형 (식 13.4.1) + 샌드위치 분산으로 회귀 계수만 robust 하게 추정한다. 연관 자체의 존재 여부는 분포 무관한 Commenges-Andersen score test (식 13.2.2) 로 검정한다.

절	핵심 도구	핵심 식	출력
13.1	Shared frailty 모형	식 13.1.1, 13.1.3 (Laplace)	모형 정의
13.2	Commenges-Andersen score test	식 13.2.2-9	\(T/\sqrt{V} \sim N(0,1)\), p-값
13.3	Gamma frailty + EM	식 13.3.1-7	\(\hat{\beta}\), \(\hat{\theta}\), \(\hat{\tau} = \theta/(\theta+2)\)
13.4	Lee et al. marginal model	식 13.4.1-3	\(\hat{\beta}\), 샌드위치 SE

Litter rat 예제 (Mantel et al. 1977) 의 세 결과 일치:

접근	결과
Score test	\(T/\sqrt{V} = 1.33\), \(p = 0.184\) — 연관 무
Gamma frailty	\(\hat{\theta} = 0.472\), Wald \(p = 0.31\), \(\hat{\tau} = 0.19\) — 연관 무
Marginal	robust SE \(= 0.300\) < naive \(0.317\) — 연관 무
공통 결론	Litter 내 종양 발생 연관은 통계적으로 유의하지 않음

Ch.13 의 위치 — 책 전체에서

Klein 교재의 마지막 본문 장. Ch.4-7 (비모수), Ch.8-9 (Cox), Ch.10 (가산), Ch.11-12 (진단·모수) 모두 개체 독립을 전제했다. Ch.13 은 그 가정을 마지막으로 풀어내며 클러스터 데이터 (longitudinal, family, recurrent event) 분석의 입문을 제공한다.

후속 학습 방향: - Recurrent event: Andersen-Gill, PWP (Prentice-Williams-Peterson), WLW 차이 - GEE for survival: Marginal model 의 더 일반화된 형태 - Joint frailty: 두 사건 (예: 입원과 사망) 의 joint analysis - Cure model: 일부 환자가 사건을 절대 경험하지 않는 모형 - Modern: DeepSurv 의 계층적 확장 — neural-net frailty

10 관련 주제

선행 지식

• Ch.8 — Cox 비례위험 모형 — 식 (13.1.1) 의 모태
• Ch.11 — 회귀 진단 — 마팅게일 잔차 (§ 13.2 score test 재료) 와 score 잔차 (§ 13.4 식 13.4.2 재료)
• § 11.3-11.4 — 마팅게일 + PH 검정
• § 11.5-11.6 — Deviance + dfbeta

관련 개념

• Ch.10 — 가산 위험 모형 — PH 의 또 다른 일반화
• Ch.12 — 모수적 회귀 모형 — frailty 의 모수적 베이스라인 옵션
• Modern Survival Analysis (RSF / DeepSurv) — 계층적 확장의 ML 버전

11 참고 문헌

• Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer. Chapter 13.
• Clayton, D. G. (1978). A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika, 65(1), 141-151. (Gamma frailty 의 원전)
• Hougaard, P. (1986a). Survival models for heterogeneous populations derived from stable distributions. Biometrika, 73(2), 387-396. (Positive stable frailty)
• Commenges, D., & Andersen, P. K. (1995). Score test of homogeneity for survival data. Lifetime Data Analysis, 1(2), 145-156. (§ 13.2 score test 의 원전)
• Lee, E. W., Wei, L. J., & Amato, D. A. (1992). Cox-type regression analysis for large numbers of small groups of correlated failure time observations. In Survival Analysis: State of the Art (pp. 237-247). Springer. (§ 13.4 marginal 모형)
• Wei, L. J., Lin, D. Y., & Weissfeld, L. (1989). Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. JASA, 84(408), 1065-1073. (WLW 확장)
• Therneau, T. M., & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer. (R coxph 의 frailty 와 cluster 옵션 매뉴얼)
• Hougaard, P. (2001). Analysis of Multivariate Survival Data. Springer. (Ch.13 주제의 단행본 확장)