1 도입 — § 13.1-13.2 가 다루는 두 질문
Ch.13 전체의 질문은 “개체 독립 가정이 깨질 때 어떻게 분석하는가” 이다. § 13.1 과 § 13.2 는 그 중 가장 기초적인 두 질문에 답한다.
| 절 | 질문 | 도구 |
|---|---|---|
| § 13.1 | 그룹 내 연관을 어떻게 모형화하는가? | Shared frailty 모형 — 보이지 않는 random effect \(u_i\) |
| § 13.2 | 그 연관이 정말 존재하는가? (모형 적합 전에 검증) | Commenges-Andersen score test — 분포 무관 검정 |
§ 13.3 (gamma frailty 추정) 과 § 13.4 (marginal model) 는 § 13.1 의 모형을 바탕으로 더 정교한 추정·대안을 제공한다. 그러나 그 출발점은 모두 § 13.1 의 모형 정의이다. 그리고 § 13.2 는 § 13.3 의 무거운 EM 알고리즘에 들어가기 전에 “그럴 가치가 있는지” 를 빠르게 점검하는 게이트 역할을 한다.
이 두 절을 깊이 이해하면 Ch.13 의 나머지가 자연스럽게 따라온다.
2 § 13.1 Frailty 모형 — 보이지 않는 그룹 효과의 도입
2.1 동기 — 독립성이 깨지는 4가지 상황
Klein 은 § 13.1 도입에서 다음 네 가지 시나리오를 제시한다.
| 시나리오 | 그룹 단위 | 공유되는 것 |
|---|---|---|
| Litter (한 배 새끼) 비교 | 어미 한 마리 | 유전·자궁환경·초기 영양 |
| 부부 코호트 | 가구 | 식습관·환경 노출·생활 패턴 |
| 같은 환자의 여러 사건 | 한 환자 | 면역계·체질·만성 질환 |
| 다발성 종양 / 양측 신장 | 한 환자의 좌우 장기 | 같은 면역계·유전 배경 |
이 모든 경우에 그룹 내 두 개체의 사건 시점이 공통 요인 을 통해 연관된다. Ch.4-12 에서 구축한 비모수·Cox·모수 모형은 모두 \(T_1, T_2\) 의 독립을 전제했으므로, 이 가정이 깨지면 분산 추정이 잘못되고 회귀 계수가 attenuate 된다.
2.3 다른 표현 — 식 (13.1.2)
식 (13.1.1) 에서 \(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} \]
이 표현이 frailty 모형의 표준 형태 이며 본 장의 모든 후속 절이 이를 기반으로 한다.
두 식은 수학적으로 동치이지만 해석의 출발점이 다르다.
| 표현 | 형태 | 직관 |
|---|---|---|
| 식 (13.1.1) | \(\exp(\sigma w_i + \beta^t Z_{ij})\) | “frailty 는 또 하나의 (보이지 않는) 공변량이다” |
| 식 (13.1.2) | \(u_i \cdot \exp(\beta^t Z_{ij})\) | “frailty 는 위험률에 곱해지는 그룹별 인자다” |
식 (13.1.2) 가 곱셈 형태로 더 직관적이라 Klein 은 § 13.3 부터는 이 표현을 주로 쓴다. 다만 § 13.2 score test 의 가설 (\(H_0: \sigma = 0\)) 은 식 (13.1.1) 의 \(\sigma\) 를 자연스럽게 가리키므로 두 표현을 모두 익숙하게 두는 것이 좋다.
2.4 \(u_i\) 가 만드는 효과 — Frail 한 사람이 먼저 죽는다
식 (13.1.2) 에서 \(u_i > 1\) 인 그룹은 위험률이 부풀려져 사건이 빨리 발생하고, \(u_i < 1\) 인 그룹은 늦게 발생한다.
이 단순한 사실에서 두 가지 흥미로운 결과가 따라온다.
시간이 지날수록 살아남은 사람들은 robust group (\(u_i\) 가 작은) 출신 일 가능성이 높아진다. 가장 frail 한 사람이 먼저 죽기 때문이다.
이는 인구 단위 (population-level) 위험률을 관찰할 때 다음 패턴을 만든다:
- 초기에는 모든 그룹이 섞여 있어 평균 위험률이 비교적 높다
- 시간이 지나면 frail 그룹이 사라지고 robust 그룹만 남아 평균 위험률이 떨어진다
- 따라서 frailty 모형은 시간에 따른 위험률 감소 를 자연스럽게 설명한다
이 selection effect 는 임상에서 자주 관찰된다. 예: 심부전 환자의 1년 후 사망률이 6 개월 시점보다 낮은 경우가 있는데, 이는 위험이 감소해서가 아니라 frail 환자가 이미 사망했기 때문이다. Cox PH 만으로는 이 패턴을 모형화하기 어려우나 frailty 가 자연스럽게 잡아낸다.
Frailty 가 있는데 일반 Cox 모형으로 적합하면 어떻게 되는가? 회귀 계수 \(\beta\) 가 0 쪽으로 끌려가는 경향이 있다 (attenuation).
이유: 관측된 효과 = 진짜 공변량 효과 + frailty 효과의 평균. 그런데 selection effect 로 시간이 지나면 frail 한 사람이 빠지므로, 후기에 관측되는 효과 크기가 작아진다. 단일 \(\beta\) 는 이 시간 평균을 추정하므로 진짜 효과보다 작게 추정된다.
오늘 13-overview 의 Mantel rat 예제 (\(\hat{\beta}_{\text{indep}} = 0.897\) vs \(\hat{\beta}_{\text{frailty}} = 0.904\)) 가 이 패턴을 보여준다. 차이는 작지만 방향은 일관 — frailty 를 고려하면 효과가 더 크게 드러난다.
2.5 Joint Survival Function — 식 (13.1.3) 과 Laplace 변환
데이터의 우도를 쓰려면 그룹 내 개체들의 결합 (joint) 생존함수 가 필요하다. \(u_i\) 가 주어졌을 때 그룹 \(i\) 의 개체들은 conditional independent 이므로:
\[ S(x_{i1}, \ldots, x_{in_i} \mid u_i) = \prod_{j=1}^{n_i} \exp\!\left[-u_i H_0(x_{ij}) e^{\beta^t Z_{ij}}\right] = \exp\!\left[-u_i \sum_{j=1}^{n_i} H_0(x_{ij}) e^{\beta^t Z_{ij}}\right] . \]
여기서 \(u_i\) 에 대해 기댓값을 취해 marginal 결합 생존함수를 얻는다. \(V = \sum_j H_0(x_{ij}) e^{\beta^t Z_{ij}}\) 로 쓰면:
\[ S(x_{i1}, \ldots, x_{in_i}) = E_U[\exp(-U V)] = \mathcal{L}_U(V) , \]
여기서 \(\mathcal{L}_U(\cdot)\) 는 frailty 분포 \(U\) 의 Laplace 변환 이다. 따라서:
\[ S(x_{i1}, \ldots, x_{in_i}) = \mathcal{L}_U\!\left[ \sum_{j=1}^{n_i} H_0(x_{ij}) e^{\beta^t Z_{ij}} \right] . \tag{식 13.1.3} \]
\(E[e^{-Uv}]\) 는 \(U\) 의 분포 \(G\) 에 대한 Laplace 변환의 정의 그 자체이다:
\[ \mathcal{L}_U(v) = E_U[e^{-Uv}] = \int_0^\infty e^{-uv} \, dG(u) . \]
frailty 모형에서 \(\exp(-u_i V)\) 형태의 식은 항상 등장하는데 (조건부 생존함수가 위험률에 곱셈 작용하므로), 거기에 \(u_i\) 적분을 취하면 자동으로 Laplace 변환이 된다. 이는 우연이 아니라 모형 구조의 필연이다.
실용적 함의:
- Frailty 분포 선택 = Laplace 변환 선택. 분포의 PDF 가 복잡해도 Laplace 변환이 폐쇄형이면 우도 계산이 가능하다.
- Gamma 가 가장 많이 쓰이는 이유는 \(\mathcal{L}(v) = (1 + \theta v)^{-1/\theta}\) 의 깔끔한 폐쇄형 때문이다. 이 한 가지 사실이 § 13.3 의 EM 알고리즘 전체를 가능하게 한다.
- 다른 분포의 Laplace 변환:
- Positive stable: \(\mathcal{L}(v) = \exp(-v^\alpha)\), \(0 < \alpha \leq 1\)
- Inverse Gaussian: \(\mathcal{L}(v) = \exp\!\left[\frac{1}{\theta}\left(1 - \sqrt{1 + 2\theta v}\right)\right]\)
- Log-normal: 폐쇄형 없음 → 수치 적분 필요 (느림)
따라서 분포의 “친화도” 는 거의 Laplace 변환의 모양으로 결정된다.
2.6 Frailty 분포의 선택 — 무엇을 쓸 것인가
문헌에 제안된 frailty 분포들:
| 분포 | 제안자 | \(\mathcal{L}(v)\) | 특징 |
|---|---|---|---|
| One-parameter Gamma | Clayton (1978), Clayton & Cuzick (1985) | \((1 + \theta v)^{-1/\theta}\) | EM 가능, Kendall \(\tau\) 폐쇄형 |
| Positive stable | Hougaard (1986a) | \(\exp(-v^\alpha)\) | Marginal 도 PH 유지 (유일) |
| Inverse Gaussian | Hougaard (1986b) | \(\exp[\frac{1-\sqrt{1+2\theta v}}{\theta}]\) | 시간이 지나면 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. Kendall \(\tau = \theta/(\theta+2)\) 폐쇄형으로 임상 해석이 직관적.
Marginal 분포도 PH 를 유지하고 싶음 → Positive stable. 식 (13.1.2) 의 conditional PH 와 marginal PH 가 동시에 성립하는 유일한 분포. 다만 모수 해석이 까다로워 실무 사용은 적다.
시간에 따른 frailty 영향 약화 (생존자 평균화) → Inverse Gaussian. 장기 추적 데이터에서 자연스럽게 selection effect 가 더 강하게 표현된다.
Linear mixed model 친숙 → Log-normal. 수치 적분 부담은 있지만 random effect 의 추가/곱 모형을 자유롭게 섞을 수 있다.
불확실하면 → Gamma 가 가장 무난. EM 알고리즘이 표준화되어 있고 (§ 13.3) 패키지 지원이 가장 좋다.
2.7 Frailty 의 두 가지 사용 — 다변량 vs Overdispersion
식 (13.1.2) 는 두 개의 다른 문제 에 같은 형태로 적용된다.
전형적 셋업: \(G\) 그룹, 그룹 \(i\) 에 \(n_i \geq 2\) 명. 같은 그룹 구성원 사이의 사건 시점 연관을 모형화.
예: - 50 litter × 3 마리 = 150 마리 쥐 (\(G=50, n_i=3\)) - 100 가구 × 2 명 = 200 명 부부 (\(G=100, n_i=2\)) - 50 환자 × 2 신장 = 100 신장 (\(G=50, n_i=2\))
\(\sigma\) (또는 \(\theta\)) 의 추정값이 곧 그룹 내 연관 강도이다.
전형적 셋업: \(G = N\) (개체 수), \(n_i = 1\) (그룹 크기 1). 즉 모든 개체가 자기 자신만의 그룹.
해석: 모형에 포함되지 않은 공변량들의 효과를 random effect 로 흡수. Cox PH 의 잔차 분산이 모형 예상보다 클 때 (overdispersion) 이를 frailty 로 표현.
식의 형태는 동일하지만 \(\sigma = 0\) 검정의 검정력은 매우 낮다 (\(n_i = 1\) 이면 그룹당 정보가 거의 없음). 따라서 이 사용은 주로 모형의 robust 화 의미이지 검정 결과로 결론짓는 분석은 아니다.
이 두 사용이 같은 식으로 표현된다는 점이 frailty 모형의 보편성을 보여준다. 식 (13.1.1) 하나로 두 문제를 다 푼다.
3 § 13.2 Score Test for Association — 분포 무관 검정
3.1 동기 — 무거운 모형 적합 전에 빠른 검증
§ 13.3 의 gamma frailty EM 알고리즘은 계산이 무겁고 (수렴 느림, 정보 행렬 거대), 더구나 frailty 분포 가정이 하나 들어간다. 만약 데이터에 정말 연관이 없다면 (즉 \(\sigma = 0\)) 이 모든 작업은 헛수고다.
Commenges & Andersen (1995) 은 이 문제에 깔끔한 해법을 제안했다: 모형 적합 없이 score test 만으로 \(\sigma = 0\) 을 검정한다.
3.2 가설과 분포 무관성
\[ H_0: \sigma = 0 \quad \text{vs} \quad H_1: \sigma > 0 . \]
핵심 매력 두 가지:
- 일반 Cox PH 모형 (frailty 없음) 만 적합하면 된다. EM, 수렴 검사 등 모두 불필요.
- frailty 분포에 어떤 가정도 하지 않는다. Gamma 든 IG 든 log-normal 이든, \(H_0\) 하에서는 모두 동일한 score 통계량으로 환원된다.
귀무가설 \(\sigma = 0\) 은 “frailty 분산이 0” 즉 “\(u_i \equiv 1\) (degenerate)” 라는 점이다. 이 점에서는 frailty 분포 자체가 사라진다 (분산 0 인 분포는 점 분포이다). 따라서 \(H_0\) 근방의 score 함수는 frailty 분포의 형태에 의존하지 않는다.
수학적으로: 임의의 분포 \(G\) 에서 \(\sigma \to 0\) 극한을 취하면 \(u_i = \exp(\sigma w_i) \to 1\) 이 된다 (확률적 의미). 따라서 score 함수의 1차 편미분 \(\partial \ell / \partial \sigma |_{\sigma=0}\) 은 어떤 \(G\) 든 동일한 형태가 된다.
이 사실 덕분에 검정 통계량이 단순해지고 실무 적용이 쉬워진다.
3.3 검정 통계량 — 식 (13.2.2) 와 (13.2.4)
귀무 모형 (일반 Cox) 을 적합해 \(\hat{\beta}\) 와 Breslow 베이스라인 \(\widehat{H}_0(t)\) 를 얻는다 (Ch.8.8). 마팅게일 잔차를 계산한다 (Ch.11.3):
\[ M_{ij} = \delta_{ij} - \widehat{H}_0(T_{ij}) \exp(\hat{\beta}^t Z_{ij}) . \]
그러면 score 통계량은:
\[ T = \sum_{i=1}^G \left( \sum_{j=1}^{n_i} M_{ij} \right)^2 - D + C , \tag{식 13.2.2} \]
여기서 \(D = \sum_{i,j} \delta_{ij}\) 는 총 사건 수, \(C\) 는 식 (13.2.3) 으로 정의되는 보정 항으로 \(N \to \infty\) 시 0 으로 수렴한다:
\[ C = \sum_{i=1}^G \sum_{j=1}^{n_i} \frac{\delta_{ij}}{[S^{(0)}(T_{ij})]^2} \sum_{b=1}^G \left[\sum_{k=1}^{n_b} Y_{bk}(T_{ij}) e^{\hat{\beta}^t Z_{bk}}\right]^2 . \tag{식 13.2.3} \]
(여기서 \(S^{(0)}(t) = \sum_{i,j} Y_{ij}(t) e^{\hat{\beta}^t Z_{ij}}\).)
식 (13.2.2) 를 다시 쓰면 통계량의 의미가 명확해진다:
\[ T = \underbrace{\sum_{i=1}^G \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\) 의 마팅게일 잔차 곱의 합이다.
- 두 잔차가 같은 방향 (둘 다 양수 또는 둘 다 음수) 이면 곱이 양수이다.
- 그룹 내 잔차들이 체계적으로 같이 움직이면 (= 연관 존재) 이 항이 커진다.
- 독립이면 \(E[M_{ij} M_{ik}] = 0\) 이라 이 항이 0 주변에서 변동한다.
이 항이 score test 의 주된 신호 항 이다. 사실상 “연관이 있다 = 그룹 내 잔차가 같이 움직인다” 의 직접 정량화이다.
(2) Overdispersion: \(\sum M_{ij}^2\) 가 표본 크기 \(N\) 보다 큰지 측정.
- PH 모형이 정확하면 마팅게일 잔차의 표본 모멘트는 \(E[M_{ij}^2] \approx 1\) 에 가까워야 한다 (이론적으로 정확히 1 은 아니지만 단위 분산 근처).
- 이 항이 양수로 크면 individual variability 가 모형 예상보다 크다는 신호.
- 단변량 frailty (그룹 크기 1) 에서는 (1) 이 0 이고 이 항만 작동한다.
(3) 보정 \(C\): 마팅게일 잔차의 추정 오차 보정.
- \(M_{ij}\) 가 진짜 마팅게일이 아니라 추정량이므로 (분모에 \(\hat{\beta}, \widehat{H}_0\) 가 들어감) 약간의 편향이 있다.
- \(C\) 는 이 편향의 1차 근사 보정이며 \(N \to \infty\) 시 0 으로 사라진다.
- 작은 표본에서는 무시할 수 없으나 큰 표본에서는 \(T \approx (1) + (2)\) 로 간소화된다.
따라서 score test 통계량 = “그룹 내 상관 정도” + “overdispersion 정도” + “유한 표본 보정” 이다.
3.4 분산 추정량 \(V\) — 식 (13.2.5)-(13.2.9)
\(T\) 의 점근 분산은 약간 복잡하다. 다음 양들을 차례로 계산한다.
Step 1. 위험집합 정규화 가중치 (식 13.2.5):
\[ p_{ij}(t) = \frac{Y_{ij}(t) \exp(\hat{\beta}^t Z_{ij})}{S^{(0)}(t)}, \quad \bar{p}_i(t) = \sum_{j=1}^{n_i} p_{ij}(t) . \]
직관: \(p_{ij}(t)\) 는 “시점 \(t\) 의 위험집합에서 개체 \((i,j)\) 가 차지하는 가중 비율” 이고, \(\bar{p}_i(t)\) 는 그룹 \(i\) 가 차지하는 비율이다.
Step 2. 각 사건 시점 \(t_k\) 에서 그룹별 양 (식 13.2.6-7):
\[ M_{ij}(t_k) = \delta_{ij} \mathbb{1}\{T_{ij} \leq t_k\} - Y_{ij}(t_k') \widehat{H}_0(t_k') e^{\hat{\beta}^t Z_{ij}} , \]
(시점 \(t_k\) 에 적합한 적절한 형태로; 본문 식 13.2.6 참조). 그룹 합 \(\bar{M}_i(t_k) = \sum_j M_{ij}(t_k)\) 를 정의한다.
Step 3. 사건 시점별 quadratic 양 (식 13.2.7):
\[ Q_i(t_k) = 2\!\left[ \bar{M}_i(t_{k-1}) - \sum_{g=1}^G \bar{M}_g(t_{k-1}) \bar{p}_g(t_k) - \bar{p}_i(t_k) + \sum_{g=1}^G \bar{p}_i(t_k)^2 \right] . \]
Step 4. 공변량 효과 정정 항 (식 13.2.8):
\[ \theta_b = \sum_{i=1}^G \sum_{k=1}^d Q_i(t_k) d_k \left[ \sum_{j=1}^{n_i} p_{ij}(t_k) Z_{ijb} \right] , \quad b = 1, \ldots, p . \]
Step 5. 분산 (식 13.2.9):
\[ V = \sum_{i=1}^G \sum_{k=1}^d Q_i(t_k)^2 \, \bar{p}_i(t_k) d_k + \sum_{b=1}^p \sum_{c=1}^p \theta_b \theta_c \widehat{\sigma}_{bc} , \tag{식 13.2.9} \]
여기서 \(\widehat{\sigma}_{bc}\) 는 \(\hat{\beta}\) 의 추정 공분산 행렬 원소.
식 (13.2.9) 의 \(V\) 는 두 부분으로 구성된다:
첫째 항 \(\sum Q_i^2 \bar{p}_i d_k\): 각 그룹의 변동성을 사건 시점 가중으로 합산. 표본 크기와 그룹 수가 클수록 안정적으로 커진다.
둘째 항 \(\sum_b \sum_c \theta_b \theta_c \widehat{\sigma}_{bc}\): \(\hat{\beta}\) 의 추정 변동성이 \(T\) 의 분산에 미치는 영향. \(\hat{\beta}\) 도 추정량이라 그 자체가 변동을 가진다는 점을 보정.
만약 \(\hat{\beta}\) 를 알려진 진짜 값으로 사용한다면 둘째 항은 사라진다. 그러나 우리는 \(\hat{\beta}\) 도 데이터에서 추정하므로 그 변동성이 \(T\) 의 분산에 propagate 된다 (델타 방법의 정신).
식이 복잡해 보이지만 본질은 단순하다 — score 통계량의 분산을 그룹별 변동 + 회귀 추정 변동으로 분해 한 것이다.
3.5 표준화 통계량과 검정 결정
\[ Z = \frac{T}{\sqrt{V}} \xrightarrow{d} N(0, 1) \quad \text{under } H_0 . \]
표준 정규 분포에서 \(|Z| > 1.96\) 또는 \(Z > 1.645\) (단측, 양의 연관만 검정) 로 5 % 수준 결정.
\(\sigma > 0\) 만 의미가 있다 (\(\sigma\) 는 척도 모수라 음수 불가). 따라서 단측 검정 이 자연스럽다 — \(H_1: \sigma > 0\).
단측이면 임계값 \(z_{0.95} = 1.645\). 양측이면 \(z_{0.975} = 1.96\). 본문 Mantel rat 예제는 양측 p-value (0.184) 를 보고 하는데, 이는 보수적 보고 관행 때문이다.
실무에서는 양측 p-value 를 보고하되, 단측이 정당화되는 맥락 (예: “litter 내 음의 연관은 의학적으로 비현실적”) 이면 단측으로 결론을 좁힐 수 있다.
3.6 절차 요약 — 5 단계 워크플로
1. 일반 Cox 모형 적합 (frailty 없음)
→ β̂, Ĥ₀(t) 획득
↓
2. 마팅게일 잔차 계산
M_ij = δ_ij - Ĥ₀(T_ij) exp(β̂'Z_ij)
↓
3. 통계량 T 계산 (식 13.2.2)
T = Σ_i (Σ_j M_ij)² - D + C
↓
4. 분산 V 계산 (식 13.2.5-9)
↓
5. 표준화 + 결정
Z = T/√V, p = P(|N(0,1)| > |Z|)EM 알고리즘 같은 반복 수렴 없음. 행렬 한 번 계산으로 끝.
3.7 예제 13.1 — Mantel et al. (1977) Litter Rat 데이터
3.7.1 데이터 셋업
50 litter 의 암컷 새끼 쥐를 사용한 종양 유발 연구. 각 litter 의 새끼 3 마리 중:
- Treated 1 마리 (drug 투여, \(Z = 1\))
- Control 2 마리 (placebo, \(Z = 0\))
총 표본 \(N = 150\). 관심 사건: 종양 발생 시점 (개월).
질문 두 가지:
- 약물 효과 가 있는가? (전통적 회귀)
- Litter 효과 가 있는가? (그룹 내 연관)
3.7.2 1단계 — 일반 Cox 모형 (frailty 무시)
단일 공변량 \(Z\) (1 if treated):
- \(\hat{\beta} = 0.8975\)
- \(\widehat{\text{Var}}(\hat{\beta}) = 0.1007\)
- HR \(= e^{0.8975} \approx 2.45\) — 약물군이 placebo 군보다 약 2.45 배 빠르게 종양 발생
같은 어미의 형제들이 연관되어 있다면 SE = \(\sqrt{0.1007} = 0.317\) 은 과소추정 일 수 있다. 형제 효과가 강할수록 유효 표본 크기가 150 보다 작아지고, 진짜 SE 는 더 커야 한다. 이를 검증하기 위해 score test 를 한다.
3.7.3 2단계 — Score Test 적용
마팅게일 잔차 \(M_{ij}\) 계산 후 식 (13.2.2) 와 (13.2.9):
- \(T = 8.91\)
- \(V = 45.03\)
- \(Z = T/\sqrt{V} = 8.91 / \sqrt{45.03} = 1.33\)
- 양측 \(p = 2 \cdot P(N(0,1) > 1.33) = 0.184\)
3.7.4 결론
\(p = 0.184 > 0.05\) → 5 % 수준에서 litter 효과의 증거 없음.
이는 다음을 의미한다:
- 같은 어미 새끼들 사이의 연관이 통계적으로 유의하지 않다.
- 따라서 일반 Cox 의 SE (\(0.317\)) 를 그대로 받아들여도 큰 문제가 없다.
- § 13.3 의 무거운 gamma frailty EM 을 굳이 돌릴 필요가 없다 — 추정해봐야 \(\hat{\sigma}\) 가 0 근처로 나올 것이다.
이 한 번의 score test 로 § 13.3 ~ § 13.4 의 모든 후속 분석 부담을 절약한다.
검정의 비유의는 “탐지하지 못했다” 일 뿐이지 “없다” 의 증명이 아니다. 다음 한계를 유념한다:
- 표본 크기: 50 litter × 3 마리 = 150. 작은 그룹 수 (50) 에 작은 그룹 크기 (3). frailty 분산 \(\theta\) 가 작을수록 검정력이 떨어진다.
- 검열: 종양이 일어나지 않은 (검열) 마리는 마팅게일 잔차가 음수로 나오는데, 그 정보가 약하다.
- 단순 모형: \(Z\) 가 단일 이항이라 다른 보정 변수가 없다.
따라서 score test 가 비유의여도 § 13.3 의 frailty 추정과 § 13.4 의 marginal model SE 를 함께 보는 것이 강건한 분석이다. Mantel rat 예제에서는 세 접근이 모두 일관 (연관 무) 이므로 결론이 단단하지만, 일반적으로 항상 그런 것은 아니다.
3.8 Score Test 의 한계와 권장
1. 분포 가정 없음: 어떤 frailty 분포 (Gamma, IG, log-normal 등) 도 가정하지 않으므로 모형 misspecification 위험이 작다.
2. 단일 적합으로 충분: 일반 Cox 한 번만 적합하면 된다. EM 같은 반복 알고리즘 불필요.
3. 그러나 추정값을 주지 않는다: \(\sigma\) 또는 \(\theta\) 의 점추정·신뢰구간을 얻지는 못한다. 강도가 알고 싶으면 § 13.3.
4. 작은 표본 / 작은 \(\theta\) 에서 검정력 약함: \(\sigma = 0\) 이 모수공간 경계라 점근 이론의 검정력이 약하다. 그룹 수 \(G < 30\) 정도면 신뢰 어려움.
4 분포 무관성을 다시 본다 — 왜 이게 가능한가
이 절의 핵심 신비 중 하나는 “어떤 frailty 분포든 같은 통계량으로 환원” 된다는 점이다. 직관적 이유를 좀 더 정리한다.
\(\sigma\) 가 척도 모수: 식 (13.1.1) 의 \(\exp(\sigma w_i)\) 에서 \(\sigma\) 는 frailty 의 변동성 폭을 결정한다.
- \(\sigma > 0\): \(w_i\) 가 확률적으로 변동 → 그룹마다 다른 \(u_i\)
- \(\sigma = 0\): \(\exp(0 \cdot w_i) = 1\) → 모든 그룹의 \(u_i = 1\) (degenerate)
따라서 \(\sigma = 0\) 인 점에서는 frailty 가 사라지고 모형이 일반 Cox 와 정확히 같다. 그 점에서의 score 함수 \(\partial \ell / \partial \sigma\) 는 frailty 분포의 모양과 무관 — 어차피 분포가 점 분포로 collapse 했으므로.
수학적으로:
\[ \left. \frac{\partial \ell(\sigma, \beta, h_0)}{\partial \sigma} \right|_{\sigma=0} \]
이 식이 분포 \(G\) 의 형태에 의존하지 않는다 (Commenges-Andersen 의 핵심 보조정리). \(G\) 에 의존하는 부분은 모두 \(\sigma\) 의 고차 항으로 들어가 1차에서는 없다.
이 사실 덕분에:
- 검정 통계량이 어떤 \(G\) 든 같은 마팅게일 잔차 형태로 표현된다 (식 13.2.2).
- 분산 \(V\) 도 \(G\) 무관하게 일반 Cox 의 양들로만 계산된다.
- “Gamma 만 검정” 또는 “IG 만 검정” 이 아니라 모든 frailty 모형에 대한 단일 검정 이다.
이는 LRT 와의 본질적 차이이다 — LRT 는 두 개의 명시 모형 (Null vs Alt) 을 비교하므로 Alt 의 분포 가정이 필요하다. Score test 는 Null 만 적합하므로 Alt 의 분포 가정이 불필요.
5 코드 예시
5.1 Step 1 — 마팅게일 잔차 + Score Test 직접 구현 (R)
R 의 survival 패키지 + 직접 식 (13.2.2) 구현. (Python lifelines 는 frailty score test 를 직접 제공하지 않으므로 R 이 우선.)
library(survival)
# Mantel rat 데이터 (예시)
# data: time, status, treatment, litter
fit0 <- coxph(Surv(time, status) ~ treatment, data = rats)
# 1. 마팅게일 잔차 (Ch.11.3)
M <- residuals(fit0, type = "martingale")
D <- sum(rats$status) # 총 사건 수
N <- nrow(rats)
# 2. 그룹별 마팅게일 잔차 합 제곱
group_M_sum <- tapply(M, rats$litter, sum)
T_paired <- sum(group_M_sum^2) - sum(M^2) # = (1) 그룹 내 쌍 상관 항
# 3. C (보정) 계산은 식 (13.2.3) 직접 구현 — 복잡
# 실무에서는 Klein 의 SAS 매크로 또는 R 의 frailty score test 패키지 사용
# 4. 표준화 통계량 (간이)
T_full <- sum(group_M_sum^2) - D # 식 (13.2.2) 의 T - C
# V 는 식 (13.2.5-9) 로 별도 계산 필요
# 본 예에서는 Klein 보고값 사용: V = 45.03
Z_stat <- T_full / sqrt(45.03)
p_val <- 2 * (1 - pnorm(abs(Z_stat)))
cat("Z =", Z_stat, ", p =", p_val, "\n")5.2 Step 2 — Frailty 모형 점검 (R coxph + frailty)
직접 score test 가 어렵다면, 같은 모형의 LRT 로 우회한다 (분포 가정이 추가되지만 결과 비슷):
5.3 Step 3 — Python lifelines 차선 (직접 구현)
lifelines 는 frailty score test 미지원. 그러나 마팅게일 잔차 기반 통계량은 직접 계산 가능:
import numpy as np
import pandas as pd
from lifelines import CoxPHFitter
# data: T, E, treatment, litter
df = pd.DataFrame({...})
# 1. 일반 Cox 적합
cph = CoxPHFitter()
cph.fit(df, duration_col="T", event_col="E", formula="treatment")
# 2. 마팅게일 잔차
M = cph.compute_residuals(df, kind="martingale").values.flatten()
D = df["E"].sum()
N = len(df)
# 3. 그룹별 합 (식 13.2.4 의 (1) 항)
df["M"] = M
group_sum = df.groupby("litter")["M"].sum().values
T_pair = (group_sum ** 2).sum() # Σ (Σ_j M_ij)²
# 4. 통계량 (보정 C 미포함, 큰 표본 근사)
T_stat = T_pair - D
# 5. 분산은 부트스트랩 또는 Klein 식 직접 구현
# 간이 부트스트랩
B = 1000
T_boot = []
for _ in range(B):
idx = np.random.choice(df["litter"].unique(),
size=df["litter"].nunique(), replace=True)
df_b = pd.concat([df[df.litter == i] for i in idx])
cph_b = CoxPHFitter().fit(df_b, "T", "E", formula="treatment")
M_b = cph_b.compute_residuals(df_b, kind="martingale").values.flatten()
df_b["M_b"] = M_b
g_b = df_b.groupby("litter")["M_b"].sum().values
T_boot.append((g_b ** 2).sum() - df_b["E"].sum())
V_boot = np.var(T_boot)
Z_stat = T_stat / np.sqrt(V_boot)
print(f"Z = {Z_stat:.3f}, p = {2*(1-stats.norm.cdf(abs(Z_stat))):.4f}")6 핵심 요약
Shared frailty 모형 \(h_{ij}(t) = h_0(t) u_i \exp(\beta^t Z_{ij})\) (식 13.1.2) 는 그룹 \(i\) 의 모든 구성원이 공유하는 latent random effect \(u_i\) 로 그룹 내 연관을 모형화한다. Commenges-Andersen score test \(T = \sum_i (\sum_j M_{ij})^2 - D + C\) (식 13.2.2) 는 일반 Cox 의 마팅게일 잔차로부터 계산되며, frailty 분포 가정 없이 \(\sigma = 0\) 을 검정한다. 통계량은 (1) 그룹 내 쌍 상관 + (2) overdispersion + (3) 보정 의 세 부분으로 분해되며, 그룹 내 쌍 상관 항이 주된 신호이다.
| 항목 | 핵심 식 |
|---|---|
| Frailty 모형 | \(h_{ij}(t) = h_0(t) \exp(\sigma w_i + \beta^t Z_{ij})\) (식 13.1.1) |
| Frailty 모형 (다른 표현) | \(h_{ij}(t) = h_0(t) u_i \exp(\beta^t Z_{ij})\) (식 13.1.2) |
| Joint survival | \(\mathcal{L}_U[\sum_j H_0(x_{ij}) e^{\beta^t Z_{ij}}]\) (식 13.1.3) |
| Score 통계량 | \(T = \sum_i (\sum_j M_{ij})^2 - D + C\) (식 13.2.2) |
| 분해 형태 | \(T = \sum \sum_{j \neq k} M_{ij} M_{ik} + (\sum M_{ij}^2 - N) + C\) (식 13.2.4) |
| 분산 | \(V\) (식 13.2.9) |
| 결정 | \(Z = T/\sqrt{V} \xrightarrow{d} N(0,1)\) |
Score test 가 유의하면 (또는 비유의여도 강도가 궁금하면) § 13.3 의 gamma frailty EM 추정 으로 진행한다. 거기서:
- \(\hat{\theta}\) 의 점추정 + SE
- Kendall \(\tau = \theta/(\theta+2)\) 로 임상 해석
- 회귀 계수 \(\hat{\beta}\) 의 frailty-adjusted 추정 (attenuation 보정)
이 모든 게 § 13.3 의 EM 알고리즘에서 한 번에 나온다. § 13.2 score test 가 그 무거운 작업의 게이트 역할을 한다는 점을 기억한다.
7 관련 주제
선행 지식
- Ch.13 Overview — Multivariate Survival Analysis — frailty / score test / marginal model 의 통합 그림
- Ch.8 — Cox 비례위험 모형 — 식 (13.1.1) 의 모태, 베이스라인 \(\widehat{H}_0\) 추정 (§ 8.8)
- § 11.3-11.4 — 마팅게일 잔차 — score test 의 입력 (\(M_{ij}\)) 정의
Ch.13 시리즈
- § 13.3-13.4 — Gamma Frailty EM Estimation & Marginal Model — 본 검정 후 강도 추정 (EM) 또는 분산만 보정 (sandwich) 의 두 길
관련 개념
- Ch.10 — 가산 위험 모형 — Cox PH 의 또 다른 일반화 (분포가 아닌 함수 형태로)
- Modern Survival Analysis (RSF / DeepSurv) — 계층적 random effect 의 ML 확장
8 참고 문헌
- Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer. § 13.1-13.2.
- Commenges, D., & Andersen, P. K. (1995). Score test of homogeneity for survival data. Lifetime Data Analysis, 1(2), 145-156. (§ 13.2 의 원전)
- 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)
- Mantel, N., Bohidar, N. R., & Ciminera, J. L. (1977). Mantel-Haenszel analyses of litter-matched time-to-response data, with modifications for recovery of interlitter information. Cancer Research, 37(11), 3863-3868. (Example 13.1 의 데이터 출처)
- Hougaard, P. (2001). Analysis of Multivariate Survival Data. Springer. (Ch.13 주제의 단행본 확장)