Klein Ch.11 — Regression Diagnostics for the Cox Model

1 도입 — 왜 회귀 진단인가

Ch.8-9에서는 Cox 비례위험(PH) 모형의 추정과 검정을 다루었지만, 모형이 정말 옳은가라는 질문은 거의 다루지 않았다. 통상의 선형 회귀에서 잔차 도표가 모형 진단의 표준 도구인 것처럼, Cox 모형에서도 잔차에 기반한 진단 절차가 필수이다.

그러나 Cox 모형의 잔차는 선형 회귀처럼 단순하지 않다. 반응 변수가 우중도절단된 시간이고, 베이스라인 위험 $h_0(t)$ 이 비모수이며, 적합 자체가 부분우도 기반이기 때문이다. 이 복잡성을 다루기 위해 여러 종류의 잔차가 제안되었으며, 각 잔차는 모형의 서로 다른 측면을 진단한다.

1.1 진단할 4가지 핵심 측면

Klein § 11.1에 따르면 Cox 모형 진단은 다음 네 가지 질문에 답해야 한다:

공변량의 함수 형태: $Z$ 의 효과를 $\exp(\beta Z)$ 로 모형화하는 게 맞는가, 아니면 $\exp(\beta \log Z)$, $\exp(\beta Z^2)$, 또는 이산화된 형태가 더 적합한가?
PH 가정의 타당성: 위험비가 정말 시간에 무관하게 일정한가?
개별 관측치의 적합도: 모형이 예측한 것보다 너무 빨리/늦게 사망한 환자(이상치)가 있는가?
영향력(leverage): 한 관측치가 추정 결과를 과도하게 좌우하는가?

각 질문에 대응하는 잔차가 따로 있다는 점이 Cox 진단의 핵심이다.

1.2 다섯 잔차 한눈에

잔차	기호	주 용도	해당 절
Cox-Snell	$r_j$	전반적 적합도 (단위 지수 검증)	§ 11.2
Martingale	$\hat{M}_j$	공변량 함수 형태 결정	§ 11.3
Schoenfeld / Score	$S_{jk}$	PH 가정 검토	§ 11.4
Deviance	$D_j$	이상치(outlier) 탐지	§ 11.5
dfbeta (∆)	$\Delta_{jk}$	영향력 관측치 식별	§ 11.6

잔차 간 관계 — 마팅게일이 중심에 있다

Cox 진단의 중심 개념은 마팅게일 잔차이다. 다른 잔차들은 마팅게일 잔차의 변형 또는 도출물이다:

Cox-Snell $r_j$ = $\delta_j - \hat{M}_j$ (사건 지시자에서 마팅게일을 뺀 것)
Score $S_{jk}$ = 마팅게일 잔차에 공변량 편차를 곱한 적분
Deviance $D_j$ = 마팅게일을 정규형으로 변환한 것
dfbeta $\Delta_{jk}$ = score 잔차를 Fisher 정보로 표준화한 것

따라서 마팅게일을 이해하면 Ch.11 전체가 한 줄기로 연결된다.

2 § 11.2 Cox-Snell 잔차 — 전반적 적합도

2.1 정의 — 식 (11.2.1)

가장 먼저 제안된 Cox 모형의 잔차이다 (Cox & Snell, 1968). Cox PH 모형 $h(t \mid Z_j) = h_0(t) \exp(\beta^\top Z_j)$ 를 적합한 후, 각 개체의 잔차를 다음과 같이 정의한다:

\[ r_j = \hat{H}_0(T_j) \exp\left(\sum_{k=1}^p Z_{jk} b_k\right), \quad j = 1, \ldots, n \tag{식 11.2.1} \]

여기서 $\hat{H}_0(t)$ 는 § 8.8의 Breslow 기저 누적 위험 추정량이다.

2.2 핵심 원리 — 확률 적분 변환

왜 Cox-Snell 잔차가 적합도를 알려주는가

확률론의 표준 정리: 임의 확률변수 $X$ 에 대해 $U = F(X) \sim \text{Uniform}(0, 1)$ 이다. 이를 위험 함수로 바꾸면 $H(X) \sim \text{Exp}(1)$ (단위 지수분포). 즉:

모형이 정확하면 — $r_j = \hat{H}(T_j \mid Z_j) \approx H(T_j \mid Z_j)$ 이고, 이는 단위 지수분포에서 추출한 표본처럼 보여야 한다.

Cox-Snell 잔차는 우중도절단을 그대로 갖고 있으므로 절단된 단위 지수 표본처럼 보여야 한다.

2.3 진단 절차

Cox 모형 적합 → $b$ 와 $\hat{H}_0(t)$ 획득
각 개체에 대해 $r_j$ 계산
$(r_j, \delta_j)$ 를 새 데이터로 보고 Nelson-Aalen 추정량 $\hat{H}_r(r)$ 계산
$\hat{H}_r(r)$ vs $r$ 도표 를 그림

도표가 원점을 지나는 45도 직선(기울기 1)에 가까우면 모형이 적합하다. 단위 지수분포의 누적 위험이 $H_E(t) = t$ 이기 때문이다.

2.4 시간의존 / 층화 모형으로 확장 — 식 (11.2.2)

시간의존 공변량 또는 층화 베이스라인을 갖는 모형에서는:

\[ r_j = \hat{H}_{0j}(T_j) \exp\left[\sum_{k=1}^p Z_{jk}(T_j) b_k\right] \tag{식 11.2.2} \]

여기서 $\hat{H}_{0j}(t)$ 는 $j$ 번째 개체에 적용되는(층 또는 시간의존 형태의) 베이스라인 누적 위험이다.

2.5 Klein Example 11.1 — 골수이식 데이터

8개 공변량 Cox 모형으로 무병 생존을 분석한 결과, 전반적 도표는 45도선과 비슷했지만 MTX 사용 vs 미사용으로 분리하여 그리니 MTX 그룹이 45도선 위로 벗어났다. 이는 PH 가정이 MTX 변수에 대해 의심스럽다는 신호였다. 시간의존 공변량 $Z_9(t) = Z_8 \ln t$ Wald 검정 결과 p = 0.032로 PH가 깨졌음이 확인되었고, MTX로 층화한 후의 Cox-Snell 도표는 다시 45도선에 잘 맞았다.

Cox-Snell 잔차의 한계

도표가 45도선에서 벗어나면 “모형이 안 맞는다”는 신호이지만, 어디서 어떻게 안 맞는지는 알려주지 않는다. 함수 형태가 잘못된 건지, PH가 깨진 건지, 분포 가정이 틀린 건지 구별하기 어렵다. 따라서 Cox-Snell은 첫 검진이며, 다른 잔차로 후속 진단이 필요하다.

3 § 11.3 마팅게일 잔차 — 공변량 함수 형태

3.1 정의 — 식 (11.3.1), (11.3.2)

마팅게일 잔차의 일반 형태:

\[ \hat{M}_j = N_j(\infty) - \int_0^\infty Y_j(t) \exp[b^\top Z_j(t)] \, d\hat{H}_0(t) \tag{식 11.3.1} \]

고정 공변량과 우중도절단만 있는 경우:

\[ \hat{M}_j = \delta_j - \hat{H}_0(T_j) \exp\left(\sum_{k=1}^p Z_{jk} b_k\right) = \delta_j - r_j \tag{식 11.3.2} \]

마팅게일 잔차의 직관

식 (11.3.2)는 단순하다:

$_j = $ 관측 사건 수($\delta_j$, 0 또는 1) $-$ 모형이 예측한 기대 사건 수($r_j$)

이는 “관측 - 기대” 형태로, 통상의 잔차 개념과 정확히 일치한다. 차이점은 다음과 같다:

최댓값: $+1$ (사건 발생, 기대값 0)
최솟값: $-\infty$ (절단, 기대값이 매우 큼)
합: $\sum \hat{M}_j = 0$ (성질)

비대칭성(매우 음수로 갈 수 있음)이 마팅게일 잔차의 특징이며, 이로 인해 이상치 탐지에는 부적합하고 deviance로 변환해서 사용한다.

3.2 핵심 용도 — 함수 형태 식별

$Z_1$ 의 함수 형태가 불확실한 상황을 가정한다. 진정한 모형은:

\[ H(t \mid Z^*, Z_1) = H_0(t) \exp(\beta^{*\top} Z^*) \exp[f(Z_1)] \tag{식 11.3.3} \]

여기서 $f(Z_1)$ 이 우리가 찾고자 하는 미지의 함수 형태이다.

진단 절차 (Therneau et al., 1990 / Fleming-Harrington, 1991):

$Z_1$ 을 제외한 $Z^*$ 만으로 Cox 모형 적합
각 개체의 마팅게일 잔차 $\hat{M}_j$ 계산
$\hat{M}_j$ vs $Z_{1j}$ 산점도 작성
LOWESS 등 평활 곡선을 적합

이론적으로 (마팅게일 잔차의 조건부 기댓값 전개):

\[ E(\hat{M}_j \mid Z_{1j}) \approx [f(Z_{1j}) - \log(g^*)] \cdot w \]

여기서 $w = $ 사건 수 / 표본 크기. 따라서 평활 곡선의 형태가 곧 $f(Z_1)$ 의 형태이다.

평활 곡선 형태	$f$ 함수 형태
직선	$f(Z) = \beta Z$ — 변환 불필요
로그 곡선	$f(Z) = \beta \log Z$
단조 변환점	이산화 ($Z \geq Z_0$ 지시 변수)
U-shape	$\beta_1 Z + \beta_2 Z^2$
비단조 / 평탄	효과 없음

3.3 Klein Example 11.2 — 림프종 이식 대기 시간

이식 대기 시간 $Z_1$ 의 함수 형태를 결정하기 위해, 다른 4개 공변량(이식 종류, 질환, 교호작용, Karnofsky)으로만 Cox 모형을 적합하고 마팅게일 잔차 vs 대기시간 LOWESS 곡선을 그렸다. 결과:

0-50개월: 잔차 약 0 (효과 없음)
50-100개월: 단조 감소
100개월+: 평탄

이는 변환점이 있는 이산화가 적절함을 시사한다. § 8.4 기법으로 최적 컷포인트를 찾으니 84개월에서 분기되었으며, scaled score 통계량 0.697(유의하지 않음)을 통해 컷포인트 위치를 확정했다.

4 § 11.4 PH 가정의 그래프적 검정

PH 가정은 Cox 모형의 핵심이지만, 이 가정이 깨지면 결과가 심각하게 왜곡된다. § 11.4는 네 가지 그래프적 검정을 제시한다.

4.1 4가지 그래프 도구

4.1.1 (1) Log-cumulative hazard 도표

공변량 $Z_1$ 으로 데이터를 $K$ 층으로 나누고, 각 층에서 베이스라인 누적 위험 $\hat{H}_{g0}(t)$ 를 추정한다. PH 가정 하에서:

\[ H_{g0}(t) = \exp(\gamma_g) H_{10}(t) \quad \Rightarrow \quad \log \hat{H}_{g0}(t) - \log \hat{H}_{10}(t) = \gamma_g \]

도표: $\log \hat{H}_{g0}(t)$ vs $t$ 의 곡선들이 수직 거리가 일정해야 한다 (평행). 또는 차이 $\log \hat{H}_{g0}(t) - \log \hat{H}_{10}(t)$ vs $t$ 도표가 수평선이어야 한다.

4.1.2 (2) Andersen 도표

$\hat{H}_{g0}(t)$ vs $\hat{H}_{10}(t)$ 의 산점도. PH 가정 하에서 원점을 지나는 직선 ($\text{slope} = e^{\gamma_g}$).

도표 형태	의미 (Gill-Schumacher, 1987)
직선 (원점)	PH 성립
볼록 (convex)	$h_g(t) / h_1(t)$ 가 시간에 따라 증가
오목 (concave)	$h_g(t) / h_1(t)$ 가 시간에 따라 감소
조각선형	위험비가 조각상수

4.1.3 (3) Arjas 도표

특정 $g$ 그룹의 누적 위험 추정량을 사건 수에 대해 그린다. PH 가정 하에서 직선이어야 한다.

4.1.4 (4) Score / Schoenfeld 잔차 도표

각 사건 시점에서 공변량의 위험집합 평균 편차를 통해 PH 위반을 진단한다.

Schoenfeld 잔차의 핵심

식 (11.4.7) — 모든 공변량이 시점 0에 고정인 경우의 score process:

\[ U_k(t) = \sum_{\text{deaths} \leq t} \left( Z_{jk} - \bar{Z}_k(T_j) \right) \tag{식 11.4.7} \]

각 사건 시점 $T_j$ 에서의 단일 항 $Z_{jk} - \bar{Z}_k(T_j)$ 가 바로 Schoenfeld 잔차이다. 의미:

“이 사건 시점에 죽은 사람의 공변량이, 그 시점 위험집합 평균과 얼마나 차이나는가”

PH 가정 하에서는 Schoenfeld 잔차의 시간 트렌드가 0이어야 한다. 트렌드가 보이면 PH가 깨진 신호이다 (Grambsch-Therneau 검정의 기반).

표준화된 score process $W_k(t) = U_k(t) \times \text{SE}(b_k)$ 는 PH 하에서 묶인 브라운 다리(tied-down Brownian bridge) 로 수렴한다. 즉 $W_k(0) = W_k(\infty) = 0$ 이고 사이에서는 무작위 진동을 보여야 한다. 절댓값이 어느 시점에서 비정상적으로 크면 PH 위반.

세 도표를 모두 사용하라

Klein은 log-cumulative, Andersen, score 잔차 세 가지를 모두 사용할 것을 권장한다. 이유:

각 도표가 다른 형태의 PH 위반에 민감하다
score 잔차의 PH 위반 검출력은 다른 도표와 비교가 충분히 되어 있지 않다
연속 공변량은 score 잔차로 자연스럽게 다루지만, log-cumulative는 이산화가 필요

세 도표가 일관된 신호를 보내면 결론이 단단해진다.

4.2 Klein Example 11.3 — 골수이식 (Auto vs Allo)

Log-cumulative 도표: 두 곡선이 교차 → PH 명백히 깨짐.

Andersen 도표: 오목 형태 → $h_{\text{allo}}(t) / h_{\text{auto}}(t)$ 가 시간에 따라 감소.

해석: 초기에는 auto가 우월(이식 관련 사망률 낮음), 후기에는 allo가 우월(재발률 낮음). 이는 단일 HR 모형으로는 절대 포착할 수 없는 시간 변동 효과이다.

5 § 11.5 Deviance 잔차 — 이상치 탐지

5.1 마팅게일 잔차의 한계

마팅게일 잔차 $\hat{M}_j \in (-\infty, 1]$ 는 매우 비대칭이다:

사건 발생자: $\hat{M}_j$ 가 양수에 가까움
절단자: $\hat{M}_j$ 가 큰 음수가 될 수 있음

이런 비대칭성 때문에 마팅게일 잔차의 절댓값을 보고 이상치를 판정하면 절단자가 부당하게 이상치로 판정될 수 있다.

5.2 Deviance 잔차 정의 — 식 (11.5.1)

이를 보정하기 위한 정규형 변환:

\[ D_j = \text{sign}[\hat{M}_j] \cdot \left\{-2 \left[ \hat{M}_j + \delta_j \log(\delta_j - \hat{M}_j) \right]\right\}^{1/2} \tag{식 11.5.1} \]

Deviance 변환의 직관

세 가지 효과를 동시에 달성한다:

부호 보존: $\text{sign}[\hat{M}_j]$ 로 양/음 정보 유지
스케일 변환: $\sqrt{\cdot}$ 으로 큰 음수를 줄이고 작은 음수의 가시성 키움
로그 보정: $\delta_j \log(\delta_j - \hat{M}_j)$ 로 사건자/절단자의 분포 차이 보정

결과적으로 $D_j$ 는 근사적으로 표준정규 분포를 따른다. 따라서 $|D_j| > 2$ 또는 $|D_j| > 3$ 인 관측치를 잠재적 이상치로 간주할 수 있다 (통상 회귀와 동일).

$\hat{M}_j = 0$ 이면 $D_j = 0$, $\hat{M}_j$ 가 1에 가까우면 $D_j$ 가 부풀려지고, 큰 음수에서는 $D_j$ 가 축소된다.

5.3 진단 절차

최종 Cox 모형 적합
각 관측치의 risk score $\sum_k b_k Z_{jk}$ 계산
$D_j$ vs risk score 산점도

경 절단(light censoring) 환경에서는 $D_j$ 가 정규 잡음처럼 보여야 한다. $|D_j|$ 가 큰 점이 잠재 이상치이다.

5.4 Klein Example 11.2 (계속)

이식 데이터에서 risk score = $-1.33 Z_1 - 2.31 Z_2 + 2.11 Z_3 - 0.056 Z_4 - 2.10 Z_5$ 로 계산했다. Deviance 도표에서 두 명의 이상치 발견:

Risk score = $-6.40$, $D_j = 2.26$ — 공변량 $(1, 0, 0, 90, 0)$, 1개월 사망
Risk score = $-7.38$, $D_j = 2.78$ — 공변량 $(0, 1, 0, 90, 0)$, 0.93개월 사망

두 환자 모두 위험 인자가 좋은데 너무 일찍 사망한 사례. Karnofsky 90 (좋은 상태), 비-NHL/HL 호의적 조합. 의학적 검토가 필요한 관측치들이다.

6 § 11.6 dfbeta — 영향력 관측치

6.1 핵심 질문

관측치 $j$ 를 빼면 추정값 $b_k$ 가 얼마나 바뀌는가?

가장 직접적인 방법: $j$ 번째 관측치를 빼고 Cox 모형을 다시 적합 → $b_{(j)}$. 차이 $b - b_{(j)}$ 가 영향력의 정의이다. 그러나 $n$ 개 관측치에 대해 이 절차를 반복하면 $n$ 번의 Cox 적합이 필요하다 — 큰 데이터에서는 비현실적이다.

6.2 Score 잔차 기반 근사

다행히도 score 잔차로 좋은 근사를 얻을 수 있다. 식 (11.6.1)에서, 모든 공변량이 시점 0에 고정인 경우 score 잔차는:

\[ S_{jk} = \delta_j[Z_{jk} - \bar{Z}_k(T_j)] - \sum_{t_b \leq T_j} [Z_{jk} - \bar{Z}_k(t_b)] \exp(b^\top Z_j) [\hat{H}_0(t_b) - \hat{H}_0(t_{b-1})] \tag{식 11.6.1} \]

식 (11.6.1)의 두 부분

이 식은 두 항으로 구성된다:

첫째 항 $\delta_j[Z_{jk} - \bar{Z}_k(T_j)]$ = Schoenfeld 부분 잔차 — 사건 발생 시점에서 공변량 편차
둘째 항 = 위험집합 누적 항으로, 개체가 위험집합에 있던 모든 사건 시점에서의 가중 편차의 합

전체적으로 score 잔차 $S_{jk}$ 는 “관측치 $j$ 가 모형 적합에 기여한 정보량”의 추정이다.

6.3 dfbeta 근사

정의:

\[ \Delta_{jk} \approx \mathbf{I}(b)^{-1} (S_{j1}, S_{j2}, \ldots, S_{jp})^\top \]

여기서 $\mathbf{I}(b) = \text{cov}(b)^{-1}$ 은 관측 Fisher 정보이다. $\mathbf{I}(b)^{-1}$ 는 추정 $b$ 의 공분산 행렬이므로 dfbeta는 score 잔차를 표준 오차의 척도로 변환한 것에 해당한다.

이는 단 한 번의 Cox 적합으로 모든 관측치의 영향력을 동시에 평가할 수 있게 해준다.

6.4 진단 절차

최종 Cox 적합 → $b$, $\mathbf{I}(b)^{-1}$
각 관측치의 score 잔차 $S_{jk}$ 계산
dfbeta 근사 $\Delta_{jk}$ 계산
각 공변량 $k$ 별로 $\Delta_{jk}$ vs 관측치 번호 또는 $\Delta_{jk}$ vs $Z_{jk}$ 산점도 작성

큰 절댓값의 $\Delta_{jk}$ 를 보이는 관측치가 영향력 관측치이다.

6.5 Klein Example 11.2 (계속)

이식 데이터에서 dfbeta 분석 결과:

관측치 17: $Z_1$, $Z_2$, $Z_3$ 추정에 가장 큰 영향. NHL/auto 환자(불리한 조합)인데 가장 오래 생존한 outlier. Risk vector $(0, 0, 0, 70, 0)$, 사망 2.633개월
관측치 12: Karnofsky 계수에 가장 큰 영향. 낮은 Karnofsky인데 오래 생존
관측치 3, 11: 좋은 Karnofsky인데 빨리 사망

영향력 관측치의 처리 원칙

영향력이 큰 관측치를 자동으로 제거하면 안 된다. Klein 권고:

통계학자와 임상의의 협력: 영향력 관측치 발견 시 데이터 기록 오류 가능성을 먼저 확인. 의학적 타당성을 임상의와 검토 후 제거 여부 결정. 단순히 “영향력이 크다”는 이유만으로 제거하면 모형 일반화 가능성이 낮아질 수 있다.

작은 표본 크기에서는 모든 관측치가 어느 정도 영향력을 갖는다는 점도 유념해야 한다.

7 종합 진단 워크플로

이 다섯 잔차를 어떻게 순서대로 사용할 것인가? Klein은 다음 절차를 권고한다.

1. 모형 적합 (Cox PH)
   ↓
2. 함수 형태 결정 (§ 11.3)
   - 후보 공변량 빼고 적합 → 마팅게일 잔차 vs 공변량 산점도
   - 평활 곡선의 형태로 함수 형태 / 이산화 결정
   ↓
3. 최종 함수 형태로 재적합
   ↓
4. PH 가정 검토 (§ 11.4)
   - log-cumulative + Andersen + Score(Schoenfeld) 도표
   - PH 위반 시: 시간의존 공변량 / 층화 / 가산 모형(Ch.10)
   ↓
5. 전반적 적합도 (§ 11.2)
   - Cox-Snell 잔차의 단위 지수 도표
   ↓
6. 이상치 탐지 (§ 11.5)
   - Deviance vs risk score 도표
   - |D_j| > 2-3 관측치 검토
   ↓
7. 영향력 검토 (§ 11.6)
   - dfbeta 도표 (각 공변량별)
   - 큰 영향력 관측치는 임상적 검토

진단의 우선순위

1단계 (필수): 함수 형태 + PH 가정 — 모형 구조 자체가 잘못되면 다른 진단이 의미 없다. 2단계 (권장): Cox-Snell 적합도 — 전반적 모형 평가. 3단계 (선택): 이상치 + 영향력 — 데이터 품질 관리.

순서를 지키지 않고 영향력부터 보면, 함수 형태가 잘못된 모형의 영향력을 분석하는 셈이 된다.

8 코드 예시

8.1 Step 1: 순수 Python으로 마팅게일 / Cox-Snell 잔차 구현

import numpy as np

# 작은 합성 데이터 (Cox 모형 적합 후)
T = np.array([1.2, 2.3, 3.0, 4.1, 5.8, 6.5, 7.2, 8.5])
delta = np.array([1, 1, 0, 1, 1, 0, 1, 1])
Z = np.array([0, 1, 0, 1, 0, 1, 0, 1]).reshape(-1, 1)

# Cox 적합 결과 (가정)
b = np.array([0.85])  # 추정 계수

# Breslow 베이스라인 누적 위험 추정
def breslow_H0(T, delta, Z, b):
    n = len(T)
    order = np.argsort(T)
    T_sorted = T[order]
    delta_sorted = delta[order]
    Z_sorted = Z[order]

    risk = np.exp(Z_sorted @ b)
    H0 = np.zeros(n)
    cum = 0.0
    for i in range(n):
        if delta_sorted[i] == 1:
            denom = risk[i:].sum()
            cum += 1.0 / denom
        H0[i] = cum

    # 원래 순서로 복원
    H0_orig = np.empty(n)
    H0_orig[order] = H0
    return H0_orig

H0_at_T = breslow_H0(T, delta, Z, b)

# Cox-Snell 잔차 식 (11.2.1)
risk_score = (Z @ b)
r = H0_at_T * np.exp(risk_score)

# 마팅게일 잔차 식 (11.3.2)
M = delta - r

print("개체 | T   | δ | Cox-Snell r | Martingale M")
for j in range(len(T)):
    print(f"  {j+1}  | {T[j]:.1f} | {delta[j]} | {r[j]:.4f}      | {M[j]:.4f}")

print(f"\nΣM_j = {M.sum():.6f}  (이론적으로 0)")

# Deviance 잔차 식 (11.5.1)
def deviance_residual(M, delta):
    eps = 1e-12
    inside = M + delta * np.log(np.where(delta - M > eps, delta - M, eps))
    return np.sign(M) * np.sqrt(-2 * inside)

D = deviance_residual(M, delta)
print(f"\nDeviance 잔차: {D}")

8.2 Step 2: lifelines로 통합 진단

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi
import matplotlib.pyplot as plt
import pandas as pd

df = load_rossi()  # 실제 생존 데이터
cph = CoxPHFitter()
cph.fit(df, duration_col='week', event_col='arrest')

# § 11.4 PH 가정 검정 — Schoenfeld 잔차 기반 Grambsch-Therneau
ph_test = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=True)
# 각 공변량별 PH 검정 + score 잔차 시간 도표

# § 11.3 마팅게일 잔차로 함수 형태 진단
martingale = cph.compute_residuals(df, kind='martingale')
# Z 변수 vs martingale 산점도 + LOWESS

import statsmodels.nonparametric.smoothers_lowess as sl
for col in ['age', 'fin']:
    smoothed = sl.lowess(martingale.values.ravel(), df[col].values, frac=0.5)
    plt.figure()
    plt.scatter(df[col], martingale, alpha=0.3)
    plt.plot(smoothed[:, 0], smoothed[:, 1], 'r-')
    plt.xlabel(col); plt.ylabel('Martingale residual')
    plt.title(f'Functional form check for {col}')

# § 11.5 Deviance 잔차로 이상치 탐지
deviance = cph.compute_residuals(df, kind='deviance')
risk_score = cph.predict_log_partial_hazard(df)
plt.figure()
plt.scatter(risk_score, deviance, alpha=0.5)
plt.axhline(2, color='r', ls='--')
plt.axhline(-2, color='r', ls='--')
plt.xlabel('Risk score'); plt.ylabel('Deviance residual')

# § 11.6 dfbeta 영향력 분석
delta_betas = cph.compute_residuals(df, kind='dfbeta')
plt.figure(figsize=(10, 4))
for col in delta_betas.columns:
    plt.plot(delta_betas[col], label=col, alpha=0.6)
plt.legend(); plt.xlabel('Observation'); plt.ylabel('dfbeta')

8.3 Step 3: R로 종합 진단 (survival 패키지)

library(survival)

fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)

# § 11.2 Cox-Snell — martingale + delta 로 환산
mart <- residuals(fit, type = "martingale")
cs <- lung$status - mart  # Cox-Snell = delta - martingale
H_cs <- survfit(Surv(cs, lung$status) ~ 1)
plot(H_cs$time, -log(H_cs$surv), type = "l")
abline(0, 1, col = "red")  # 45도 기준선

# § 11.3 함수 형태 — null 모형 마팅게일 잔차
null_fit <- coxph(Surv(time, status) ~ 1, data = lung)
null_mart <- residuals(null_fit, type = "martingale")
plot(lung$age, null_mart); lines(lowess(lung$age, null_mart), col = "red")

# § 11.4 PH 가정 — Grambsch-Therneau 검정
ph <- cox.zph(fit)
print(ph); plot(ph)

# § 11.5 Deviance
dev <- residuals(fit, type = "deviance")
plot(predict(fit), dev)

# § 11.6 dfbeta
dfb <- residuals(fit, type = "dfbeta")
matplot(dfb, type = "l")

9 핵심 요약

Ch.11 한 줄 요약

Cox 모형 진단은 다섯 가지 잔차로 네 가지 측면을 점검한다 — Cox-Snell(전반 적합), 마팅게일(함수 형태), Schoenfeld/score(PH 가정), deviance(이상치), dfbeta(영향력). 모든 잔차의 중심에는 마팅게일 $\hat{M}_j = \delta_j - \hat{H}_0(T_j) e^{b^\top Z_j}$ 가 있으며, 다른 잔차는 이의 변형 또는 도출물이다.

절	잔차	용도	핵심 식
11.2	Cox-Snell	전반 적합도	식 (11.2.1)
11.3	Martingale	함수 형태	식 (11.3.2)
11.4	Score / Schoenfeld	PH 검정	식 (11.4.7)
11.5	Deviance	이상치	식 (11.5.1)
11.6	dfbeta	영향력	식 (11.6.1)

10 관련 주제

선행 지식

후속 주제

관련 개념

Ch.10 — 가산 위험 모형 — PH 위반 시 대안
Schoenfeld-Grambsch-Therneau PH 검정