Klein § 13.3-13.4 — Gamma Frailty EM Estimation & Marginal Model

1 도입 — 두 절의 핵심 대비

§ 13.2 score test 가 연관의 존재 만 검정했다면, § 13.3 과 § 13.4 는 연관이 있을 때 어떻게 추정할 것인가 를 다룬다. 두 절은 같은 문제에 정반대 철학으로 접근한다.

절	철학	출력
§ 13.3 Frailty	“그룹 효과를 모형에 직접 추가한다.” Conditional PH 가정.	\(\hat{\beta}\) (frailty-adjusted) + \(\hat{\theta}\) (연관 강도) + Kendall \(\tau\)
§ 13.4 Marginal	“그룹 효과는 모형에 안 넣고 분산만 보정한다.” Marginal PH 가정.	\(\hat{\beta}\) (independence working) + robust SE

직관 — 두 철학의 차이

Frailty 접근: “보이지 않는 그룹 효과 \(u_i\) 가 진짜 있다고 모형에 적자. 그러면 잔차에 \(u_i\) 의 흔적이 사라지고 \(\hat{\beta}\) 도 정확해진다.”

• 장점: 연관 강도 \(\theta\) 자체가 추정값으로 나옴. 임상 해석 (Kendall \(\tau\)) 가능.
• 단점: frailty 분포 (Gamma 등) 가정. EM 알고리즘 무거움.

Marginal 접근: “\(u_i\) 의 정체는 모르겠지만, 그게 그룹 내 연관을 만드는 건 안다. \(\hat{\beta}\) 점추정은 일반 Cox 가 일치성 있게 주니까 받아들이고, 분산만 보정하자.”

• 장점: 분포 가정 없음. 일반 Cox 한 번 + 샌드위치 공식.
• 단점: 연관 강도는 모름. \(\hat{\beta}\) 의 의미가 다름 (conditional 이 아닌 marginal effect).

언제 어느 쪽: 그룹 효과 자체가 관심 (예: 가족력의 강도) → § 13.3. 회귀 계수만 관심, 연관은 nuisance → § 13.4. 둘 다 적용해 결과 일관성 보면 가장 강건.

이 한 줄 요약을 머리에 두고 두 절의 디테일을 본다.

2 § 13.3 Gamma Frailty Estimation

2.1 동기 — Score Test 후 무엇이 더 필요한가

§ 13.2 의 score test 는 “연관이 있냐 없냐” 만 답한다. 그러나 다음 질문에는 답하지 못한다:

• 연관의 강도 가 어느 정도인가? (Kendall \(\tau = 0.1\) 인지 \(0.5\) 인지)
• 연관을 고려한 회귀 계수 \(\hat{\beta}\) 는 일반 Cox 와 얼마나 다른가? (attenuation 정도)
• 베이스라인 위험률 \(\widehat{H}_0(t)\) 가 frailty 보정 후 어떻게 바뀌는가?

이 셋을 한꺼번에 답하는 것이 § 13.3 의 EM 알고리즘이다.

2.2 Gamma 분포의 선택 — 식 (13.3.1)

다양한 frailty 분포 후보 중 (§ 13.1 참조) gamma 가 EM 알고리즘이 가장 자연스럽게 작동 하므로 본 절에서는 gamma 만 다룬다. \(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} \]

“One-parameter” gamma 의 의미

표준 gamma 분포는 형상 모수 \(k\) 와 척도 모수 \(\lambda\) 의 두 모수를 가진다. 그러나 frailty 모형에서는 평균 = 1 제약 (식별성 위해) 을 두므로 한 모수만 자유롭다.

\(E[U] = k \cdot \lambda = 1\) 제약 하에 \(\theta = \text{Var}(U)\) 만 추정한다. 따라서:

• \(k = 1/\theta\) (형상)
• \(\lambda = \theta\) (척도)
• \(\theta \to 0\): \(U\) 가 1 주변에 집중 → 거의 frailty 없음 → 일반 Cox PH 로 환원
• \(\theta \to \infty\): \(U\) 의 분산이 매우 큼 → 그룹별 위험률 차이 매우 큼 → 강한 연관

이 단일 모수 \(\theta\) 가 연관 강도를 완전히 결정한다.

2.3 Joint Survival 의 폐쇄형

식 (13.1.3) 에 gamma 의 Laplace 변환 \(\mathcal{L}(v) = (1 + \theta v)^{-1/\theta}\) 을 대입하면:

\[ 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} . \]

직관 — 폐쇄형의 가치

다른 frailty 분포 (log-normal 등) 는 결합 생존함수가 적분으로만 표현되어 우도 계산에 수치 적분이 필요하다 (느림). Gamma 는 이 적분이 폐쇄형 \((1 + \theta v)^{-1/\theta}\) 으로 풀려 우도가 깔끔하다.

이 한 가지 사실이 § 13.3 의 EM 알고리즘 전체를 가능하게 한다. Conjugate 구조 + 폐쇄형 marginal = E-step·M-step 모두 명시적 식.

2.4 Kendall’s \(\tau\) — 임상 해석의 도구

연관 강도를 임상 친화적으로 표현하는 양:

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

\(\theta\)	Kendall \(\tau\)	해석
0	0	독립
0.5	0.20	약한 양의 연관
2	0.50	중간 연관
8	0.80	강한 연관
\(\infty\)	1	완전 연관

왜 \(\theta\) 자체보다 \(\tau\) 를 보고하는가

\(\theta\) 는 척도 모수라 임상 의미가 모호 — “\(\theta = 0.5\)” 가 강한지 약한지 즉답하기 어렵다. 반면 Kendall \(\tau\) 는 0 (독립) ~ 1 (완전 일치) 의 정규화된 척도라 직관적이다.

또한 \(\tau\) 는 분포 무관 측도이다 (rank-based). 다른 frailty 분포로 적합해도 \(\tau\) 의 추정값은 비슷한 범위로 나오므로 분포 선택의 부담이 줄어든다.

실무 보고에서는 \(\hat{\theta}\) 점추정과 \(\hat{\tau}\) 둘 다 제시하는 것이 표준이다.

2.5 Marginal Log-Likelihood — 식 (13.3.2)

데이터로부터 우도를 쓸 때 \(u_i\) 는 latent 이므로 적분으로 없앤 marginal log-likelihood 를 사용한다. \(D_i = \sum_j \delta_{ij}\) 를 그룹 \(i\) 의 사건 수로 두면:

\[ \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=1}^{n_i} H_0(T_{ij}) e^{\beta^t Z_{ij}} \right] \\ &+ \sum_{j=1}^{n_i} \delta_{ij} \big[ \beta^t Z_{ij} + \ln h_0(T_{ij}) \big] \Big\} . \end{aligned} \tag{식 13.3.2} \]

식 (13.3.2) 의 4 부분 분해

이 우도가 복잡해 보이지만 4 부분으로 나누어 읽으면 명확하다:

(1) \(D_i \ln \theta - \ln \Gamma(1/\theta) + \ln \Gamma(1/\theta + D_i)\): gamma 분포 정규화 상수에서 온 항. \(\theta\) 만 의존, \(\beta\) 와 \(H_0\) 무관.

(2) \(-(1/\theta + D_i) \ln[1 + \theta \sum H_0(T_{ij}) e^{\beta^t Z_{ij}}]\): gamma Laplace 변환의 로그. frailty 모형의 핵심 항 — 그룹 내 위험률 합이 어떻게 marginal 우도에 들어가는지 보여준다.

(3) \(\sum \delta_{ij} \beta^t Z_{ij}\): 사건자의 공변량 효과. 일반 Cox 와 동일.

(4) \(\sum \delta_{ij} \ln h_0(T_{ij})\): 사건자의 베이스라인 기여.

핵심 관찰: \(\theta = 0\) 으로 보내면 (2) → \(-\sum_i \sum_j H_0(T_{ij}) e^{\beta^t Z_{ij}}\) 가 되어 일반 Cox 의 expected event count 항으로 환원된다. (1) 도 사라진다 (gamma 가 점 분포로 collapse).

2.6 모수적 베이스라인의 경우 — 직접 MLE

만약 \(h_0(t)\) 를 모수적 (예: Weibull \(h_0(t) = \alpha \lambda t^{\alpha-1}\)) 으로 두면, 식 (13.3.2) 의 모든 양이 \((\theta, \beta, \lambda, \alpha)\) 의 닫힌 함수가 된다. 그러면:

\[ (\hat{\theta}, \hat{\beta}, \hat{\lambda}, \hat{\alpha}) = \arg\max_{\theta, \beta, \lambda, \alpha} L(\theta, \beta, \lambda, \alpha) . \]

수치 최적화 (Newton-Raphson, BFGS 등) 로 직접 풀면 끝. 분산은 정보 행렬의 역수.

이 접근은 작은 표본에서 장점이 있다 (모수 가정이 정보를 보강) 그러나 베이스라인 분포 misspecification 의 위험을 떠안는다. 본 장은 비모수 베이스라인에 집중한다.

2.7 비모수 베이스라인 — EM 알고리즘이 등장하는 이유

베이스라인 \(h_0\) 를 비모수로 두면 식 (13.3.2) 가 무한 차원이라 직접 최대화 불가능. EM 알고리즘이 이 문제를 푸는 표준 도구 이다.

핵심 아이디어 — \(u_i\) 가 관측 가능하다면 우도가 단순해진다. 그래서 \(u_i\) 의 conditional 분포를 추정해 가며 우도를 최대화하는 반복 절차를 짠다.

만약 \(u_i\) 를 안다고 가정하면 완전 우도 (full likelihood):

\[ L_{\text{FULL}}(\theta, \beta, H_0; u) = L_1(\theta; u) + L_2(\beta, H_0; u) , \]

여기서

\[ L_1(\theta; u) = -G\left[\frac{1}{\theta}\ln\theta + \ln\Gamma\!\left(\frac{1}{\theta}\right)\right] + \sum_{i=1}^G \left[\frac{1}{\theta} + D_i - 1\right]\ln u_i - \frac{u_i}{\theta} , \tag{식 13.3.3} \]

\[ L_2(\beta, H_0; u) = \sum_{i=1}^G \sum_{j=1}^{n_i} \delta_{ij}\big[\beta^t Z_{ij} + \ln h_0(T_{ij})\big] - u_i H_0(T_{ij}) e^{\beta^t Z_{ij}} . \tag{식 13.3.4} \]

식 분리의 의미

\(L_{\text{FULL}}\) 이 \(L_1(\theta; u) + L_2(\beta, H_0; u)\) 로 깔끔하게 분리된다는 점이 EM 알고리즘 설계의 결정적 장점이다.

• \(L_1\): \(\theta\) 와 \(u_i\) 만 포함. \(\beta, H_0\) 무관 → \(\theta\) 추정 시 \(\beta, H_0\) 고려 안 해도 됨.
• \(L_2\): \(\beta, H_0\) 와 \(u_i\) 만 포함. \(\theta\) 가 사라짐 → \(\beta, H_0\) 추정이 일반 Cox + offset 형태가 됨 (offset = \(\ln u_i\)).

이 분리 덕분에 M-step 이 두 단계로 나뉘어 각각 표준적 방법 (gamma MLE + Cox partial likelihood) 으로 처리 가능하다.

2.8 EM 의 E-step — 식 (13.3.5), 베이지언 갱신

\(u_i\) 의 conditional 분포를 데이터·현재 모수 추정값으로 결정. Gamma 가 conjugate prior 이므로 conditional 분포도 gamma 가 된다:

\[ u_i \mid \text{data}, (\hat{\theta}, \hat{\beta}, \widehat{H}_0) \sim \text{Gamma}(A_i, C_i) , \]

여기서

\[ A_i = \frac{1}{\hat{\theta}} + D_i , \quad C_i = \frac{1}{\hat{\theta}} + \sum_{j=1}^{n_i} \widehat{H}_0(T_{ij}) e^{\hat{\beta}^t Z_{ij}} . \]

따라서 conditional 평균과 로그 평균:

\[ E[u_i \mid \text{data}] = \frac{A_i}{C_i}, \quad E[\ln u_i \mid \text{data}] = \psi(A_i) - \ln C_i , \tag{식 13.3.5} \]

여기서 \(\psi\) 는 digamma 함수.

직관 — \(\hat{u}_i\) 가 무엇을 측정하는가

\(\hat{u}_i = A_i / C_i\) 의 분자와 분모를 보자:

분자 \(A_i = 1/\theta + D_i\): - \(D_i\) 는 그룹 \(i\) 의 관측 사건 수 — 이 그룹에서 실제로 죽은 사람 수. - \(1/\theta\) 는 prior 정보 (베이지언 prior shape).

분모 \(C_i = 1/\theta + \sum H_0(T_{ij}) e^{\beta^t Z_{ij}}\): - \(\sum H_0(T_{ij}) e^{\beta^t Z_{ij}}\) 는 그룹 \(i\) 의 누적 위험 노출 — 모형이 예측한 기대 사건 수. - \(1/\theta\) 는 prior 정보 (베이지언 prior rate).

따라서:

\[ \hat{u}_i = \frac{1/\theta + D_i}{1/\theta + \text{기대 사건 수}} \approx \frac{\text{관측 사건 수}}{\text{기대 사건 수}} \quad (\theta \text{ 가 작을 때}) . \]

해석: \(\hat{u}_i\) 는 “이 그룹이 모형 예측보다 얼마나 더/덜 죽었나” 의 비율.

• 관측 > 기대 → \(\hat{u}_i > 1\) → frail group (위험률 부풀려짐)
• 관측 < 기대 → \(\hat{u}_i < 1\) → robust group (위험률 줄어듦)
• 관측 = 기대 → \(\hat{u}_i = 1\) → 평균 그룹

\(1/\theta\) 가 prior 로 들어간 형태는 shrinkage 를 만든다 — 작은 그룹 (\(D_i\) 작음) 은 prior 쪽으로 끌려가 \(\hat{u}_i \approx 1\). 큰 그룹은 데이터 쪽으로 끌려가 \(\hat{u}_i\) 가 더 극단으로 간다. 이는 random effect 의 자연스러운 정규화이다.

2.9 EM 의 M-step — 식 (13.3.6), (13.3.7)

E-step 의 \(\hat{u}_i\) 를 known 값으로 대입한 expected \(L_2\) 를 최대화:

\[ E[L_2(\beta, H_0) \mid \text{data}] = \sum_{i=1}^G \sum_{j=1}^{n_i} \delta_{ij}\big[\beta^t Z_{ij} + \ln h_0(T_{ij})\big] - \frac{A_i}{C_i} H_0(T_{ij}) e^{\beta^t Z_{ij}} . \tag{식 13.3.6} \]

이를 최대화하는 \(\beta\) 와 \(H_0\) 를 찾는다. 핵심 통찰:

M-step 은 사실상 일반 Cox + frailty offset 이다

식 (13.3.6) 을 다시 쓰면:

\[ \sum_{i,j} \delta_{ij}\big[\beta^t Z_{ij} + \ln h_0(T_{ij})\big] - \hat{u}_i H_0(T_{ij}) e^{\beta^t Z_{ij}} . \]

이 형태는 각 개체 \(j\) 의 위험률에 \(\hat{u}_i\) 가 곱해진 일반 Cox 우도 와 동치이다. 즉:

\[ \tilde{h}_{ij}(t) = \hat{u}_i \cdot h_0(t) \cdot e^{\beta^t Z_{ij}} = h_0(t) e^{\beta^t Z_{ij} + \ln \hat{u}_i} . \]

따라서 \(\ln \hat{u}_i\) 를 알려진 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} \]

여기서 \(t_{(k)}\) 는 \(k\) 번째 사건 시점, \(m_{(k)}\) 는 그 시점의 사건 수, \(R(t_{(k)})\) 는 위험집합, \(S_{(k)}\) 는 사건자 공변량 합.

식 (13.3.7) 은 Breslow partial likelihood 와 정확히 같은 형태에 위험집합 가중치만 \(\hat{u}_b\) 가 곱해진 것이다. 즉 frail group 의 개체는 위험집합 합에 더 큰 가중치로 들어간다. Standard Cox 코드를 약간 수정해 구현 가능.

베이스라인 점프도 update:

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

이는 Breslow estimator 에 frailty 가중치만 추가된 형태.

\(\theta\) 의 update 는 expected \(L_1\) 로:

\[ L_4(\theta) = E[L_1(\theta) \mid \text{data}] = -G\!\left[\frac{1}{\theta}\ln\theta + \ln\Gamma\!\left(\frac{1}{\theta}\right)\right] + \sum_{i=1}^G \left[\frac{1}{\theta} + D_i - 1\right]\big[\psi(A_i) - \ln C_i\big] - \frac{A_i}{\theta C_i} . \]

이는 1차원 함수라 1D Newton-Raphson 또는 Brent’s method 로 빠르게 풀린다.

2.10 전체 EM 알고리즘 — 단계별

Step 0 (초기화):
    β = 0, θ = 1, H_0 = Nelson-Aalen (frailty 무시)

Step 1 (E-step): 식 (13.3.5)
    A_i = 1/θ + D_i
    C_i = 1/θ + Σ_j Ĥ_0(T_ij) exp(β'Z_ij)
    û_i = A_i / C_i
    E[ln u_i] = ψ(A_i) - ln C_i

Step 2 (M-step):
    Step 2a: û_i 를 offset 으로 두고 식 (13.3.7) Cox partial
             likelihood 최대화 → β update
    Step 2b: 식 (13.3.8) 로 ĥ_k0 update
    Step 2c: 식 L_4(θ) 1D 최대화 → θ update

Step 3: |θ_new - θ_old|, ||β_new - β_old|| < ε 까지 Step 1~2 반복

EM 의 수렴은 느릴 수 있다

EM 알고리즘은 단조 우도 증가는 보장하지만 수렴 속도가 매우 느린 경우 가 있다. 특히:

• \(\theta\) 가 0 근처에 있을 때 (frailty 거의 없음): \(\hat{u}_i\) 가 모두 1 근처라 정보가 적어 \(\theta\) 추정이 크게 변하지 않음.
• 그룹 수가 작을 때 (\(G < 30\)): 정보 부족.
• 우도의 곡률이 평탄할 때.

수십~수백 회 iteration 이 필요할 수 있고, 수렴 진단은 우도 변화 외에도 모수 변화로 함께 봐야 한다.

2.11 Nielsen et al. (1992) Profile EM — 가속

표준 EM 의 느린 수렴 문제를 해결하기 위해 Nielsen 등은 다음 절차를 제안했다:

1. \(\theta\) 의 격자 (예: \(\theta = 0.05, 0.1, 0.2, 0.5, 1, 2, 5\)) 를 미리 잡는다.
2. 각 격자 점 \(\theta\) 마다 그것을 고정하고 \(\beta_\theta\) 만 EM 으로 추정 (M-step 만 반복).
3. profile log-likelihood \(L(\theta, \hat{\beta}_\theta)\) 를 격자 점들에서 평가.
4. profile likelihood 를 최대화하는 \(\theta\) 를 최종 추정값으로.

왜 profile 방식이 더 빠른가

표준 EM 은 매 iteration 에서 \(\theta\) 와 \(\beta\) 를 동시에 update 하는데, 이 둘이 강하게 상관되어 있어 (서로 영향) 진동·느린 수렴 발생.

Profile EM 은 \(\theta\) 를 외부 grid 로 분리 — 각 \(\theta\) 에서 \(\beta\) 만 풀면 conditional 문제라 빠르게 수렴. 격자가 충분히 촘촘하면 (\(\theta\) 의 정확한 최대값 위치를 1D 탐색 가능) 표준 EM 보다 5~10 배 빠르다.

두 방식 모두 같은 MLE 에 수렴한다 (Nielsen et al. 1992).

2.12 분산 추정 — Andersen et al. (1997)

EM 수렴 후 \((\hat{\theta}, \hat{\beta}, \widehat{H}_0)\) 의 분산은 식 (13.3.2) 의 marginal log-likelihood 의 관측 정보 행렬 의 역수이다. 정보 행렬은 \(D + p + 1\) 차 정방 (사건 시점 수 \(D\) + 회귀 계수 \(p\) + frailty 모수 1).

행렬의 원소들 (식 본문 참조) 은 베이스라인 점프 \(h_{k0}\), 회귀 계수 \(\beta_v\), frailty \(\theta\) 의 모든 짝에 대한 2차 편미분 \(-\partial^2 L / \partial \cdot \partial \cdot\) 으로 채운다.

비표준 분산 — 완전한 점근 정당화 부족

이 방식은 비표준이다 — frailty 의 marginal 우도가 conventional 정칙 조건을 일부 위반 (\(\theta = 0\) 이 모수공간 경계). 그러나 Andersen et al. (1997) 의 경험적·시뮬레이션 결과 (Morsing 1994) 가 충분히 정당화한다.

흥미로운 점: 같은 방식을 일반 Cox 모형에 적용 하면 정보 행렬의 역수가 정확히 표준 Cox 분산 공식이 나온다. 즉 frailty 가 0 으로 가는 극한에서 표준 이론과 일치.

작은 표본 (\(G < 30\)) 에서는 분산 추정이 신뢰 어려울 수 있어 부트스트랩으로 보강하는 것을 권장.

2.13 Example 13.1 (계속) — Mantel Litter Rat 의 Gamma Frailty 적합

§ 13.2 의 Mantel rat 데이터 (50 litter × 3 마리, 처치/대조) 에 EM 적용.

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

2.13.1 결과 해석

세 가지 관찰

1. Frailty 검정 (Wald): \(\hat{\theta} / SE(\hat{\theta}) = 0.472 / 0.462 = 1.02\), \(p \approx 0.31\) → 유의 X.

§ 13.2 의 score test (\(p = 0.184\)) 와 일관. 두 검정 모두 litter 효과를 통계적으로 탐지하지 못함.

2. \(\hat{\beta}\) 의 변화 (attenuation 보정): 독립 모형 \(0.897\) → frailty 모형 \(0.904\). frailty 를 고려하면 약물 효과가 약간 더 강하게 추정.

차이가 작은 (0.007) 이유는 \(\hat{\theta}\) 자체가 작아 (0.472) 그룹 내 연관이 약하기 때문. \(\theta\) 가 더 컸으면 \(\hat{\beta}\) 의 차이도 컸을 것.

3. Kendall’s \(\tau\): \(\hat{\tau} = 0.472 / (0.472 + 2) = 0.19\), SE(\(\hat{\tau}\)) \(\approx 0.43\).

점추정은 약한 양의 연관을 시사하지만 SE 가 커서 유의 X. 표본 크기 (50 litter) 가 frailty 추정에 충분치 않은 한계 — 연관 강도를 정밀하게 추정하려면 더 많은 그룹이 필요하다 (보통 \(G \geq 100\) 권장).

학습 포인트 — Score test vs Wald 의 일관성

같은 데이터에서:

• § 13.2 score test: \(p = 0.184\)
• § 13.3 Wald (frailty 적합 후): \(p = 0.31\)

두 검정의 p-value 가 다른 이유는 검정 방식 차이 (score 는 null hypothesis 만 적합, Wald 는 alternative 적합). 그러나 결론 (유의 X) 은 일관.

작은 표본에서는 score test 가 일반적으로 검정력이 약간 더 강하다 (\(p\) 가 더 작게 나온다). 그러나 강한 연관 (큰 \(\theta\)) 이 있을 때는 Wald 가 더 강해진다.

3 § 13.4 Marginal Model — 분산만 보정하는 길

3.1 동기 — Frailty 의 한 가지 한계

§ 13.3 의 frailty 모형은 conditional PH 를 가정한다:

\[ h_{ij}(t \mid u_i) = h_0(t) u_i e^{\beta^t Z_{ij}} . \]

즉 “그룹 효과 \(u_i\) 가 주어진 조건부에서” PH 가 성립. 그런데 우리가 실제로 관측하는 것은 marginal (집단 평균) 분포이고, marginal hazard 는 일반적으로 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 모형에서 추정한 \(\hat{\beta}\) 는 conditional 효과 (“같은 frailty 그룹 내에서 \(Z = 1\) vs \(Z = 0\) 의 위험비”) 라는 점이다. marginal hazard ratio (“전체 모집단에서 \(Z = 1\) vs \(Z = 0\) 의 위험비”) 와 다르다.

3.2 Lee, Wei, Amato (1992) 의 출발점

만약 우리의 진정한 관심이 marginal effect 라면 어떻게 해야 하는가? Lee et al. 는 다른 접근을 제안했다:

“조건부 모형은 만들지 말자. 각 개체의 marginal hazard 만 PH 로 가정하고, 그룹 내 연관은 분산 보정으로 처리하자.”

식 (13.4.1):

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

이는 일반 Cox PH 와 형태가 똑같다. 단지 그룹 내 연관 \(E[T_{ij}, T_{ik}] \neq 0\) (\(j \neq k\)) 을 허용할 뿐.

Frailty vs Marginal — 핵심 대비 다시

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

두 접근 모두 \(\hat{\beta}\) 의 점추정은 비슷 한 경우가 많지만, frailty 의 \(\hat{\beta}\) 와 marginal 의 \(\hat{\beta}\) 가 같은 양 을 추정하지 않는다는 점은 명심.

중간 강도의 연관 (\(\theta < 1\)) 에서는 두 추정값이 거의 같아 실무 차이가 작다. 강한 연관 (\(\theta > 2\)) 에서는 차이가 클 수 있어 어느 쪽을 보고할지 결정해야 한다.

3.3 Independence Working Model — 일반 Cox 그대로 적합

추정 절차의 첫 단계는 놀랍도록 단순:

1. 모든 개체를 독립으로 가정 하고 일반 Cox partial likelihood 최대화.
2. → \(\hat{\beta}\) 의 점추정 획득.

Lee et al. (1992) 가 증명한 핵심: 이 \(\hat{\beta}\) 는 marginal 모형 (식 13.4.1) 이 옳다면 일치 추정량 (consistent). 즉 그룹 내 연관 여부와 무관하게 점추정은 정확.

직관 — 왜 점추정은 영향 받지 않는가

부분우도의 score 함수 \(U(\beta) = \sum_i \sum_j \int [Z_{ij} - \bar{Z}(t)] dN_{ij}(t)\) 는 그룹 내 연관 여부와 무관하게 평균 0 인 마팅게일 (marginal model 이 옳을 때). 따라서 \(U(\hat{\beta}) = 0\) 의 해 \(\hat{\beta}\) 는 일치 추정량.

그러나 분산 은 다르다. Score 함수의 분산이 그룹 내 연관 때문에 부풀려진다. 이를 보정하지 않으면 SE 가 잘못된다.

3.4 Score 잔차 — 식 (13.4.2)

분산 보정에 필요한 핵심 양은 score 잔차 (Ch.11.6 식 11.6.1 참조). \(j\) 번째 그룹의 \(i\) 번째 개체의 \(k\) 번째 공변량에 대한 score 잔차:

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

Score 잔차의 두 항 다시 보기

(Ch.11.6 의 직관을 재현)

첫째 항 \(\delta_{ij}[Z_{ijk} - \bar{Z}_k(T_{ij})]\): 사건이 발생한 경우만 (\(\delta_{ij}=1\)) 활성화. 자기 자신의 사건 시점에서 위험집합 평균과의 편차. 이는 § 11.4 의 Schoenfeld 부분 잔차와 정확히 같다.

둘째 항 \(-\sum_{t_b \leq T_{ij}}[Z_{ijk} - \bar{Z}_k(t_b)] d\widehat{H}_0(t_b)\): 자신이 위험집합에 있던 모든 사건 시점 에서의 편차의 누적 합. 다른 사건들에 자신이 기여한 정보의 양.

종합: score 잔차 = “자기 사건 시점의 편차” - “다른 사건들에 기여한 정보”.

이는 부분우도 score function 의 \(j\) 번째 개체의 기여분이며, 합하면 정확히 \(U(\hat{\beta}) = 0\) 이 된다.

3.5 Sandwich Variance — 식 (13.4.3)

그룹 내 연관을 잡아내는 핵심 양 \(C\) 는 score 잔차의 그룹 내 합 외적:

\[ 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} \]

여기서 안쪽 두 합 \(\sum_j \sum_m\) 이 같은 그룹 \(i\) 내의 두 개체 \(j, m\) 의 score 잔차 곱을 모은다. 따라서:

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

각 그룹의 score 잔차 합의 outer product 를 더한 것.

\(C\) 가 그룹 내 연관을 어떻게 잡아내는가

독립 인 경우 (\(u_i\) 효과 없음): \(S_{ij}\) 와 \(S_{im}\) (\(j \neq m\)) 이 독립이라 \(E[S_{ij} S_{im}] = E[S_{ij}] E[S_{im}] = 0\). 따라서 \(C \approx \sum_i \sum_j S_{ij}^2\) 정도 (대각만 살아남음). 이 값은 일반 정보 행렬 \(\widehat{V}\) 와 거의 같음.

연관 이 있을 때 (\(u_i > 1\) 인 그룹): 그 그룹 내 \(S_{ij}\) 들이 같은 방향으로 치우침 → \(\sum_j S_{ij}\) 가 0 이 아닌 큰 값이 됨 → outer product 가 큰 값 → \(C\) 가 부풀려짐.

Sandwich variance:

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

• \(\widehat{V}\): 일반 정보 행렬 역수 (독립 working model 의 SE 제곱)
• 두 빵 (\(\widehat{V}\)) 사이에 속재료 (\(C\)) 를 끼운 형태 → “샌드위치”
• 그룹 내 연관이 있으면 $C > $ 일반 분산 → \(\widetilde{V} > \widehat{V}\) → SE 가 보정으로 커짐

Sandwich Form 의 보편성

\(\widehat{V} C \widehat{V}\) 형태는 통계학 전반에서 등장한다:

• White (1980): heteroskedasticity-robust SE for OLS
• GEE (Liang & Zeger, 1986): marginal model for longitudinal data
• GMM (Hansen, 1982): generalized method of moments
• Pseudo-likelihood: misspecified model 의 robust 분산

핵심은 모두 같다: 모형이 정확하지 않을 때도 작동하는 robust 분산. 식 (13.4.3) 은 이 패턴을 Cox 모형에 적용한 것이다.

실용적 함의: \(\widehat{V} C \widehat{V}\) 가 “표준” 이라는 점에서 통계 패키지 (R coxph 의 cluster(), Python lifelines 의 cluster_col 등) 에서 직접 지원한다.

3.6 검정 — Wald, Score, LRT 모두 가능

\(\widetilde{V}\) 를 SE 로 사용해 표준 Wald 검정:

\[ Z = \frac{\hat{\beta}_k}{\sqrt{\widetilde{V}_{kk}}} \xrightarrow{d} N(0, 1) . \]

다중 가설 동시 검정:

\[ \chi^2 = (\hat{\beta}_J)^t \widetilde{V}_J^{-1} \hat{\beta}_J \sim \chi^2_{|J|} . \]

LRT 도 가능하지만 일반적으로 Wald 가 더 자주 사용된다 (계산 단순).

3.7 Example 13.1 (계속) — Mantel Litter Rat 의 Marginal 적합

같은 50 litter 데이터:

• 일반 Cox: \(\hat{\beta} = 0.897\), \(\widehat{V} = 0.1007\)
• 식 (13.4.3) 으로 계산: \(C = 8.893\)
• \(\widetilde{V} = 0.1007^2 \times 8.893 = 0.0902\)
• robust SE \(= \sqrt{0.0902} = 0.3000\)

3.7.1 95% 신뢰구간 비교 (HR 척도)

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

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

이 데이터에서 샌드위치 분산 (0.0902) 이 naive 분산 (0.1007) 보다 작다. 보정 후 CI 가 좁아진 셈.

일반적인 직관 (“연관이 있으면 분산이 부풀려진다”) 과 반대 결과. 가능한 해석:

약한 음의 연관 가능성: Litter 내 새끼들 사이에 약한 음의 상관 (한 마리가 약하면 다른 마리가 강한 경쟁 효과 등) 이 있다면 score 잔차의 그룹 합이 평균 0 에 가까워 \(C\) 가 일반 분산보다 작아질 수 있다.

우연 변동: \(G = 50\) 의 작은 표본에서 \(C\) 자체의 변동성이 크다. 이 결과를 강하게 해석하지 말 것.

이 예제는 score test (\(p = 0.184\)), frailty Wald (\(p = 0.31\)), marginal robust CI 모두 “litter 효과 강한 증거 없음” 으로 일관된 결론을 준다. 어느 하나의 결과만 두고 결론짓지 말고 세 접근을 모두 비교하는 것이 강건 분석.

3.8 Wei-Lin-Weissfeld (WLW) 확장 — 그룹별 다른 베이스라인

Wei, Lin, Weissfeld (1989) 는 식 (13.4.1) 을 더 일반화 — 그룹별로 다른 베이스라인 \(h_{0i}(t)\) 를 허용:

\[ h_{ij}(t \mid Z_{ij}) = h_{0i}(t) \exp(\beta^t Z_{ij}) , \quad i = 1, \ldots, G . \]

언제 WLW 확장이 필요한가

Recurrent event 분석이 대표적 예 — 한 환자의 첫 사건과 두 번째 사건이 다른 베이스라인을 따른다 (첫 사건 후 환자 상태 변화). 또는 multi-state 모형에서 transition 마다 다른 베이스라인.

이 경우 식 (13.4.1) 의 단일 \(h_0(t)\) 가정은 너무 강하다. WLW 는 각 transition (또는 각 사건 발생 순서) 마다 별도 stratified Cox 를 적합하고 모든 stratum 에 공통인 \(\beta\) 를 추정한다.

추정·분산 구조는 식 (13.4.1) 과 같은 샌드위치 패턴이고, 단지 stratification 이 추가될 뿐이다. 자세한 절차는 Lin (1993) 참조.

4 통합 비교 — 두 접근의 결과 일관성

같은 데이터에 두 접근을 모두 적용해 결과를 비교하는 것이 권장 분석. Mantel rat 데이터의 세 접근 종합:

접근	\(\hat{\beta}\)	SE	\(p\)-value	연관 강도 추정
일반 Cox (독립)	0.897	0.317	0.005	—
Score test (§ 13.2)	—	—	0.184 (단측 0.092)	정성적: “약함”
Gamma frailty (§ 13.3)	0.904	0.323	0.005	\(\hat{\theta} = 0.472\), \(\hat{\tau} = 0.19\)
Marginal model (§ 13.4)	0.897	0.300 (robust)	0.003	— (분산만 보정)

통찰 — 일관된 결론에서의 강건성

세 접근이 모두 “약물은 종양 위험을 약 2.4 배 증가시킴 (\(\hat{\beta} \approx 0.9\), HR \(\approx 2.5\), \(p < 0.01\))” 을 일관되게 보고. 이 일관성이 결론의 강건성을 보장한다.

만약 세 접근이 갈라진다면 (예: Cox \(p = 0.05\), marginal \(p = 0.20\)) 그룹 내 연관이 결과에 영향을 미치고 있다는 신호이며, 더 깊은 진단 (frailty 분포 점검, score residual 도표 등) 이 필요하다.

실무 권고: 세 접근의 점추정·SE·p-value 를 표로 함께 보고. 어느 하나만으로 결론짓지 말 것.

5 코드 예시

5.1 Step 1 — Python: lifelines 의 cluster_col (Marginal model)

from lifelines import CoxPHFitter
import pandas as pd

# data: T, E, treatment, litter
df = pd.DataFrame({...})

# § 13.4 Marginal model — robust SE (식 13.4.3 자동 계산)
cph_robust = CoxPHFitter()
cph_robust.fit(df, duration_col="T", event_col="E",
               formula="treatment", cluster_col="litter")
print(cph_robust.summary)
# robust SE 출력. cluster_col 없으면 naive SE.

# 비교: cluster_col 없는 일반 Cox
cph_naive = CoxPHFitter()
cph_naive.fit(df, duration_col="T", event_col="E", formula="treatment")
print("Naive SE:", cph_naive.standard_errors_)
print("Robust SE:", cph_robust.standard_errors_)

lifelines 는 frailty (§ 13.3 EM) 를 직접 지원하지 않으므로 R 사용을 권장.

5.2 Step 2 — R: coxph 의 frailty / cluster

library(survival)

# § 13.3 Gamma frailty (Therneau)
fit_frailty <- coxph(Surv(T, E) ~ treatment + frailty(litter, dist = "gamma"),
                     data = rats)
summary(fit_frailty)
# theta 추정값 + LRT p-value 자동 출력

# § 13.4 Marginal model (cluster)
fit_marginal <- coxph(Surv(T, E) ~ treatment + cluster(litter),
                      data = rats)
summary(fit_marginal)
# robust SE 자동 출력 (식 13.4.3)

# 세 접근 비교 표
fit_naive <- coxph(Surv(T, E) ~ treatment, data = rats)
results <- data.frame(
  Method = c("Naive Cox", "Frailty", "Marginal"),
  Beta = c(coef(fit_naive), coef(fit_frailty)[1], coef(fit_marginal)[1]),
  SE = c(sqrt(vcov(fit_naive)),
         sqrt(vcov(fit_frailty)[1,1]),
         sqrt(vcov(fit_marginal)[1,1]))
)
print(results)

5.3 Step 3 — coxme 로 log-normal frailty (Gamma 가 아닌 경우)

library(coxme)

# Log-normal frailty (random intercept by litter)
fit_lognormal <- coxme(Surv(T, E) ~ treatment + (1 | litter), data = rats)
print(fit_lognormal)
# random effect variance 추정

# Gamma vs log-normal 비교
AIC(fit_frailty)
AIC(fit_lognormal)

5.4 Step 4 — 직접 EM 알고리즘 구현 (학습용 minimal)

em_gamma_frailty <- function(T, E, Z, group, max_iter = 100, tol = 1e-5) {
  # 초기화
  beta <- 0
  theta <- 1
  G <- length(unique(group))
  D_i <- tapply(E, group, sum)  # 그룹별 사건 수

  for (iter in 1:max_iter) {
    # 일반 Cox 적합 (offset = log(u_hat))
    # ... (생략, lifelines/coxph 의 offset 옵션 사용)

    # E-step: 식 (13.3.5)
    H0_T <- ...  # Breslow 베이스라인
    risk_sum <- tapply(H0_T * exp(beta * Z), group, sum)
    A_i <- 1/theta + D_i
    C_i <- 1/theta + risk_sum
    u_hat <- A_i / C_i

    # M-step: 식 (13.3.7) Cox + offset
    offset <- log(u_hat[group])
    fit <- coxph(Surv(T, E) ~ Z + offset(offset))
    beta_new <- coef(fit)

    # theta update: 1D Newton on L_4
    theta_new <- ...  # optimize() 사용

    if (abs(theta_new - theta) < tol && abs(beta_new - beta) < tol) break
    beta <- beta_new
    theta <- theta_new
  }
  list(beta = beta, theta = theta, u_hat = u_hat, iters = iter)
}

(완전 구현은 길어 생략. 학습 목적은 R coxph(frailty=...) 가 어떤 일을 하는지 이해하는 것.)

6 핵심 요약

§ 13.3-13.4 한 줄 요약

§ 13.3 Gamma frailty 는 conditional PH 가정 하에서 EM 알고리즘으로 \((\hat{\theta}, \hat{\beta}, \widehat{H}_0)\) 를 동시 추정한다. E-step 의 \(u_i \mid \text{data} \sim \text{Gamma}(A_i, C_i)\) (conjugate posterior, 식 13.3.5) 와 M-step 의 Cox partial likelihood + frailty offset (식 13.3.7) 이 핵심 메커니즘이다. § 13.4 Marginal model 은 conditional 가정을 버리고 marginal PH 만 가정 (식 13.4.1), 일반 Cox 점추정에 score 잔차 기반 샌드위치 분산 \(\widetilde{V} = \widehat{V} C \widehat{V}\) (식 13.4.3) 을 결합한다. 두 접근은 같은 데이터의 다른 측면을 보는 도구이며, 함께 보면 결론이 단단해진다.

절	핵심 식	출력
13.3.1	Gamma frailty 분포	\(\theta = \text{Var}(U)\)
13.3.2	Marginal log-likelihood	\(L(\theta, \beta, H_0)\)
13.3.5	E-step (\(u_i \mid \text{data}\))	\(\hat{u}_i = A_i/C_i\)
13.3.7	M-step (Cox + offset)	\(\hat{\beta}\), \(\widehat{H}_0\)
Kendall	\(\tau = \theta/(\theta+2)\)	임상 해석
13.4.1	Marginal PH 모형	\(\hat{\beta}\) (working)
13.4.2	Score 잔차 (= 식 11.6.1)	\(S_{ijk}\)
13.4.3	\(C\) 행렬	그룹 내 연관 잡아냄
Sandwich	\(\widetilde{V} = \widehat{V} C \widehat{V}\)	robust SE

Mantel Rat 의 세 접근 결과 일관

접근	결과
Score test (§ 13.2)	\(T/\sqrt{V} = 1.33\), \(p = 0.184\) — 연관 무
Gamma frailty (§ 13.3)	\(\hat{\theta} = 0.472\), Wald \(p = 0.31\), \(\hat{\tau} = 0.19\) — 연관 무
Marginal (§ 13.4)	robust SE \(= 0.300\) < naive \(0.317\) — 연관 무

세 접근 모두 “litter 효과 강한 증거 없음 + 약물 효과 HR ≈ 2.5 (\(p < 0.01\))” 일관 보고.

7 관련 주제

선행 지식

• Ch.13 Overview — Multivariate Survival Analysis — 세 접근 (score test / frailty / marginal) 의 통합 그림
• § 13.1-13.2 — Frailty Models & Score Test — 본 포스트의 직접 선행
• Ch.8 — Cox 비례위험 모형 — 식 (13.4.1) 의 모태, Breslow 베이스라인
• § 11.5-11.6 — Deviance + dfbeta — score 잔차 (식 11.6.1) 의 정의와 직관

Ch.13 시리즈

• § 13.5 — 연습문제 풀이 — Batchelor-Hackett 화상 환자 allograft + McGilchrist-Aisbett kidney catheter 두 데이터에 § 13.1-13.4 도구 종합 적용

관련 개념

• Ch.10 — 가산 위험 모형 — Cox PH 의 또 다른 일반화
• GEE for Survival — Marginal model 의 더 일반화 (Liang & Zeger)
• 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.3-13.4.
• 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 의 원전)
• Clayton, D. G., & Cuzick, J. (1985). Multivariate generalizations of the proportional hazards model. JRSS A, 148(2), 82-117. (Gamma frailty 추정의 초기 work)
• Nielsen, G. G., Gill, R. D., Andersen, P. K., & Sørensen, T. I. A. (1992). A counting process approach to maximum likelihood estimation in frailty models. Scand. J. Statist., 19(1), 25-43. (EM 알고리즘 + profile EM)
• Klein, J. P. (1992). Semiparametric estimation of random effects using the Cox model based on the EM algorithm. Biometrics, 48(3), 795-806.
• Andersen, P. K., Klein, J. P., Knudsen, K. M., & y Palacios, R. T. (1997). Estimation of variance in Cox’s regression model with shared gamma frailties. Biometrics, 53(4), 1475-1484. (분산 추정 정당화)
• 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 확장)
• Lin, D. Y. (1993). MULCOX: A computer program for the Cox regression analysis of multiple failure time data. Computer Methods and Programs in Biomedicine, 40(4), 279-293. (FORTRAN 구현)
• 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 주제의 단행본 확장)