Kwangmin Kim - Klein § 11.3–11.4 — Martingale Residuals & Graphical PH Checks

1 도입 — 함수 형태와 PH 가정

이 두 절은 Cox 모형 진단의 가장 중요한 두 가지 측면을 다룬다:

측면	질문	진단 도구 (절)
함수 형태	$Z$ 의 효과가 선형인가, 비선형인가, 이산화해야 하는가?	마팅게일 잔차 (§ 11.3)
PH 가정	위험비가 시간에 무관하게 일정한가?	log-cum / Andersen / Arjas / score (§ 11.4)

함수 형태가 잘못된 모형으로 PH 가정을 검정하는 것은 무의미하므로, 함수 형태 → PH 검정의 순서가 권고된다.

§ 11.2와의 연결

§ 11.2 Cox-Snell 잔차는 전반적 적합도를 보여주지만 위반의 종류를 알려주지 않는다. § 11.3-11.4는 이 한계를 정확히 보완한다:

함수 형태 위반 → 마팅게일 도표가 선형이 아님
PH 위반 → log-cumulative 도표가 평행하지 않음, score process가 한계선을 벗어남

2 § 11.3 — 마팅게일 잔차로 함수 형태 결정

2.1 정의 — 식 (11.3.1)

마팅게일 잔차의 일반 형태 (시간의존 공변량 포함):

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

여기서:

$N_j(t) = \mathbb{1}\{T_j \leq t, \delta_j = 1\}$: 시점 $t$ 까지 사건 발생 지시 과정
$Y_j(t) = \mathbb{1}\{T_j \geq t\}$: $t$ 직전 위험집합 포함 지시
$\hat{H}_0(t)$: Breslow 베이스라인 누적 위험

2.2 단순화 — 식 (11.3.2)

고정 공변량 + 우중도절단 환경에서:

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

여기서 $r_j$ 는 § 11.2의 Cox-Snell 잔차이다. 즉 마팅게일 = 사건 지시자 − Cox-Snell.

“관측 − 기대” 구조

식 (11.3.2)는 통상의 잔차 개념과 정확히 같다:

\[ \hat{M}_j = \underbrace{\delta_j}_{\text{관측 사건 수 (0 또는 1)}} - \underbrace{r_j}_{\text{모형 예측 기대 사건 수}} \]

예시: 어떤 환자의 Cox-Snell 잔차 $r_j = 0.3$ 이라면 모형은 그 환자가 약 0.3건의 사건을 경험했어야 한다고 예측한다 (실제 한계 [0, 약간 초과]).

사건 발생 ($\delta_j = 1$): $\hat{M}_j = 1 - 0.3 = 0.7$ — 모형이 위험을 약간 과소예측
절단 ($\delta_j = 0$): $\hat{M}_j = 0 - 0.3 = -0.3$ — 사건이 안 일어났는데 모형은 일어났을 거라 예측

마팅게일 잔차의 부호와 크기는 모형 예측과 관측의 괴리를 직접 보여준다.

2.3 핵심 성질

성질	의미
$\sum_j \hat{M}_j = 0$	평균 0
$\hat{M}_j \in (-\infty, 1]$	비대칭 분포
큰 표본에서 무상관 표본처럼 행동	잔차 도표 해석 가능

마팅게일 잔차의 비대칭성은 사건 발생 시 최댓값 $+1$ 만 가능하지만 절단 시 큰 음수가 가능하기 때문이다. 이 비대칭성으로 인해 마팅게일은 이상치 탐지에는 부적합하며 (deviance 변환으로 해결, § 11.5), 함수 형태 결정에는 적합하다.

2.4 함수 형태 결정 — 핵심 절차

미지의 공변량 $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)$ 이 우리가 찾는 미지의 함수. $Z^*$ 는 함수 형태가 이미 알려진 다른 공변량들이며, $Z_1$ 은 $Z^*$ 와 독립이라고 가정한다.

함수 형태 진단 4단계

Step 1: $Z_1$ 을 제외하고 $Z^*$ 만으로 Cox 모형 적합

Step 2: 마팅게일 잔차 $\hat{M}_j$ 계산 (식 11.3.2)

Step 3: $(\hat{M}_j, Z_{1j})$ 산점도 작성

Step 4: LOWESS 또는 spline 평활 곡선 적합

평활 곡선의 형태가 곧 $f(Z_1)$ 의 형태이다.

2.5 왜 LOWESS 곡선 = 함수 형태인가

이론적 근거 (Fleming-Harrington 1991, Therneau et al. 1990):

식 (11.3.3) 모형이 옳다고 가정하고 $g(Z_1) = \exp[f(Z_1)]$ 로 두면, 마팅게일의 조건부 기댓값은 근사적으로:

\[ E[M_j(t) \mid Z_{1j}] \approx \left[1 - \frac{g^*}{g(Z_{1j})}\right] E[N_j(t) \mid Z_{1j}] \]

여기서 $g^*$ 는 위험집합 분포에서의 $g(Z_1)$ 의 평균. Taylor 전개로 단순화하면:

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

여기서 $w = $ 사건 수 / 표본 크기. 핵심은 상수 $-\log(g^*) \cdot w$ 를 빼면 평활 곡선이 정확히 $f(Z_{1j})$ 에 비례한다는 것이다.

LOWESS 곡선 형태	권장 함수 형태 $f$
직선	$\beta Z$ — 변환 불필요
로그 모양	$\beta \log Z$
U-shape / 포물선	$\beta_1 Z + \beta_2 Z^2$
단조 변환점 (S-curve)	이산화 $\mathbb{1}\{Z \geq Z_0\}$
평탄 (수평)	$Z_1$ 효과 없음 — 모형에서 제거

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

데이터: 비-Hodgkin 림프종 / Hodgkin 환자의 동종 vs 자가 골수이식 (§ 1.10).

목표: 이식 대기 시간 $Z_1$ (월) 의 함수 형태 결정.

Step 1: $Z_1$ 을 빼고 4개 공변량(이식 종류, 질환, 교호작용, Karnofsky)으로 Cox 적합.

Step 2-4: 마팅게일 잔차 vs $Z_1$ 산점도 + LOWESS.

결과 (Figure 11.4):

잔차 ↑
     │  ●   ●  ●  ●
   0 │──●─●─●─●─────────────
     │              ●  ●
     │                ●●●
  -1 │                    ●●●●●  ───
     │
     └─────┬──────┬──────┬─────→ Z₁ (월)
           50    100    150

세 구간이 식별된다:

0-50개월: 평탄 (효과 없음)
50-100개월: 단조 감소
100개월 이상: 다시 평탄 (변동 큼)

이 패턴은 단순 선형이나 로그가 아닌 변환점이 있는 이산화를 시사한다.

최적 컷포인트: § 8.4의 기법으로 $Z_0 = 84$ 개월 발견. 이 컷포인트의 scaled score 통계량은 0.697 (유의하지 않음) — 즉 컷포인트의 위치가 데이터 의존적으로 우연 발견된 것이 아니다.

최종 모형:

\[ h(t) = h_0(t) \exp(\beta^{*\top} Z^* + \beta_5 \mathbb{1}\{Z_1 \geq 84\}) \]

함수 형태 진단의 실전 가치

이 예제의 통찰:

연속 공변량을 자동으로 선형으로 모형화하지 마라. 마팅게일 잔차 도표를 먼저 보면 비선형성이 자주 발견된다
LOWESS의 곡선 형태 → 변환의 형태. 단조 감소 + 평탄 패턴 → 이산화 + 단일 컷포인트
컷포인트는 § 8.4 기법으로 통계적으로 결정. 시각적 추정만으로 결정하지 말 것

3 § 11.4 — PH 가정의 그래프 검정

비례위험(PH) 가정은 Cox 모형의 핵심이며, 위반 시 단일 $\beta$ 추정값은 시간 평균된 효과만 반영하여 의미가 모호해진다. § 11.4는 PH 가정을 검토하는 다섯 가지 그래프 도구를 제시한다.

3.1 다섯 도구 한눈에

도구	그리는 대상	PH 하 기대 형태
log-cumulative	$\log \hat{H}_{g0}(t)$ vs $t$	평행한 곡선들
차이 도표	$\log \hat{H}_{g0}(t) - \log \hat{H}_{10}(t)$ vs $t$	수평선
Andersen plot	$\hat{H}_{g0}(t)$ vs $\hat{H}_{10}(t)$	원점 통과 직선
Arjas plot	$N_g(t)$ vs $\text{TOT}_g(t)$	45도선
Score (Schoenfeld)	표준화된 $W_k(t)$ vs $t$	묶인 브라운 다리

3.2 도구 1-2 — Log-cumulative 및 차이 도표

3.2.1 핵심 원리

연속 공변량 $Z_1$ 을 $K$ 층으로 이산화 (또는 이산 공변량의 $K$ 수준)하고, 각 층에서 PH 모형을 층화 적합. 층 $g$ 의 베이스라인 누적 위험을 $\hat{H}_{g0}(t)$ 라 하자.

PH 가정 하에서:

\[ H_{g0}(t) = \exp(\gamma_g) H_{10}(t) \]

양변에 로그:

\[ \log H_{g0}(t) - \log H_{10}(t) = \gamma_g \quad (\text{시간에 무관한 상수}) \]

두 도표의 직관

Log-cumulative 도표: 모든 층의 $\log \hat{H}_{g0}(t)$ 를 같은 시간축에 그린다. PH 하에서는 곡선들이 수직 거리가 일정 — 즉 평행해야 한다.

차이 도표 (권장): $\log \hat{H}_{g0}(t) - \log \hat{H}_{10}(t)$ 를 $t$ 에 대해 그린다. PH 하에서는 수평선이어야 한다.

차이 도표가 더 우수한 이유: 사람의 눈은 “두 곡선이 평행한가”보다 “곡선이 수평인가”를 훨씬 정확히 판단한다.

3.3 도구 3 — Andersen Plot

Andersen (1982)이 제안한 도표.

PH 하에서 $H_{g0}(t) = e^{\gamma_g} H_{10}(t)$ 이므로, $\hat{H}_{g0}(t)$ vs $\hat{H}_{10}(t)$ 의 산점도는 원점을 지나는 기울기 $e^{\gamma_g}$ 직선이어야 한다.

도표 형태 (Gill-Schumacher 1987)	의미
직선 (원점)	PH 성립, 기울기 $\approx e^{\gamma_g}$
볼록 (convex)	$h_{g0}(t) / h_{10}(t)$ 가 시간에 따라 증가
오목 (concave)	$h_{g0}(t) / h_{10}(t)$ 가 시간에 따라 감소
조각선형	위험비가 조각상수

“상대 추세 함수”의 추정

이론적으로 Andersen plot은 상대 추세 함수 $y = H_g^{-1}[H_1(t)]$ 의 추정이다. PH 하에서 이 함수는 직선이며, Dabrowska et al. (1989)는 이 도표 기반의 형식 검정도 제안했다.

3.4 도구 4 — Arjas Plot

Arjas (1988)가 제안한 도표.

3.4.1 정의 — 식 (11.4.1), (11.4.2)

기존 Cox 모형 (공변량 $Z^*$, $Z_1$ 제외)에서 적합된 누적 위험 $\hat{H}(t \mid Z_j^*)$ 를 사용. 각 사건 시점 $t_i$ 에서, $Z_1$ 의 각 층 $g$ 별로:

\[ \text{TOT}_g(t_i) = \sum_{Z_{1j} = g} \hat{H}\left(\min(t_i, T_j) \mid Z_j^*\right) \tag{식 11.4.1} \]

\[ N_g(t_i) = \sum_{Z_{1j} = g} \delta_j \, \mathbb{1}\{T_j \leq t_i\} \tag{식 11.4.2} \]

도표: $N_g(t_i)$ vs $\text{TOT}_g(t_i)$.

Arjas plot의 직관 — 마팅게일 관점

$Z_1$ 이 모형에 필요 없다면, $Z^*$ 만의 모형의 마팅게일 잔차 $M_j(t) = N_j(t) - \int_0^t Y_j(u) e^{\beta^\top Z_j^*(u)} dH_0(u)$ 는 평균 0이다. 즉:

\[ E[N_j(t)] = E\left[\int_0^t Y_j(u) e^{\beta^\top Z_j^*(u)} dH_0(u)\right] \]

각 층 $g$ 에서 합산하면 $\sum_{Z_{1j} = g} E[N_j(t)] = \text{TOT}_g(t)$ 의 추정량과 같아야 한다. 따라서 $N_g$ vs $\text{TOT}_g$ 도표는 45도선 이어야 한다.

3.4.2 도표 해석

Arjas 도표 형태	의미
모든 곡선 45도선	$Z_1$ 효과 없음, 모형에서 제거 가능
직선이지만 기울기 $\neq 1$	$Z_1$ 효과 있고 PH 성립 → 모형에 추가
비선형 (볼록/오목)	$Z_1$ 효과 있고 PH 위반 → 층화 또는 시간의존

3.4.3 등가 표현 — Total Time on Test

식 (11.2.2)의 Cox-Snell 잔차로 보면, Arjas plot은 각 층별 Total Time on Test 도표이다 (Barlow & Campo, 1975).

3.5 도구 5 — Score (Schoenfeld) 잔차

연속 공변량을 이산화 없이 다루는 가장 강력한 도구.

3.5.1 사전 정의

위험집합에서의 공변량 평균:

\[ \bar{Z}_k(t) = \frac{\sum_{j=1}^n Y_j(t) Z_{jk}(t) \exp[b^\top Z_j(t)]}{\sum_{j=1}^n Y_j(t) \exp[b^\top Z_j(t)]} \tag{식 11.4.3} \]

마팅게일 잔차 (시점 $t$ 까지):

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

3.5.2 Score 잔차 — 식 (11.4.5)

$k$ 번째 공변량과 $j$ 번째 개체의 score 잔차:

\[ S_{jk}(t) = \int_0^t [Z_{jk}(u) - \bar{Z}_k(u)] \, d\hat{M}_j(u) \tag{식 11.4.5} \]

전체 score process:

\[ U_k(t) = \sum_{j=1}^n S_{jk}(t) \tag{식 11.4.6} \]

3.5.3 고정 공변량 단순화 — 식 (11.4.7)

모든 공변량이 시점 0에 고정된 경우:

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

Schoenfeld 잔차의 정체

식 (11.4.7)의 각 사건 시점 $T_j$ 에서의 항:

\[ \boxed{Z_{jk} - \bar{Z}_k(T_j)} \]

이것이 Schoenfeld (1982) 부분 잔차이다. 의미:

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

PH 가정 하에서 사건이 무작위로 발생한다면 $\bar{Z}_k(T_j)$ 가 정확히 $Z_{jk}$ 의 기댓값이며, Schoenfeld 잔차의 평균은 0이다. 시간 트렌드가 보이면 PH 위반 — 이것이 Grambsch-Therneau (1994) 검정의 핵심이다.

3.5.4 Score Process 도표

표준화된 process:

\[ W_k(t) = U_k(t) \times \text{SE}(b_k) \]

이론 (Therneau et al. 1990): PH 하에서 $W_k(t)$ 는 묶인 브라운 다리(tied-down Brownian bridge) 로 약수렴 (단, $\text{Cov}(b_k, b_{k'}) = 0, k \neq k'$).

성질:

$W_k(0) = 0$
$W_k(\infty) = 0$ (부분우도 정규방정식)
사이에서 무작위 진동
5% 검정 임계값: $\sup_t |W_k(t)| > 1.3581$ 이면 PH 기각 (브라운 다리 supremum 분포)

3.5.5 장단점

항목	내용
장점	연속 공변량 이산화 불필요 + 단일 모형 적합으로 모든 공변량 동시 검정
단점	검출력이 Andersen/Arjas 대비 충분히 비교되지 않음
권고	세 도표를 모두 사용하여 일관된 결론 확보

3.6 Klein Examples — BMT 데이터

3.6.1 Example 11.3 — 동종 vs 자가 이식 (단일 이항 공변량)

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

차이 도표 (Figure 11.6): 시간에 따라 변동 → 초기 auto 우월, 후기 allo 우월. 단일 HR로는 포착 불가능.

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

Arjas plot (Figure 11.14): 두 곡선 모두 45도선에서 비선형 이탈 → 비례위험 모형 부적합.

Score process (Figure 11.17): $\pm 1.3581$ 한계선을 벗어남 → 5% 수준에서 PH 기각.

해석: 초기에는 auto가 유리(이식 관련 사망률 낮음), 후기에는 allo가 유리(재발률 낮음). § 7.6의 Renyi 통계량 사용 권고.

3.6.2 Example 11.1 — MTX (8공변량 모형)

Andersen + Arjas (Figure 11.13, 11.15): 두 도표 모두 비선형 이탈 → MTX의 PH 위반 확인.

Score process (Figure 11.18): 한계선 약간 벗어남 → 비례위험 의심.

대기 시간 변수 (Figure 11.16, 11.19): 4개 층 모두 45도선 / 한계선 내 → PH 성립 + 효과 없음 → 모형에서 제거 가능.

PH 검정의 5가지 실무 원칙 (Klein 권고)

1. 그래프적 검정 우선: 형식 검정은 큰 표본에서 거의 항상 기각된다 (PH는 근사일 뿐). 도표로 위반의 형태와 정도를 판단하라.

2. 차이 도표 > log-cumulative 단순 도표: 평행성보다 수평성이 시각적으로 판별 쉽다.

3. Arjas plot의 장점: 기본 Cox 모형 재적합 없이 층화 방법을 변경 가능.

4. 우측 꼬리는 신중히: 베이스라인 추정의 변동이 가장 크다. 큰 시간에서의 이탈은 과해석 금지.

5. 세 도표 종합: log-cumulative + Andersen + Arjas + score 중 적어도 3개를 사용. 한 도표만으로 결론 내리지 말 것.

4 종합 — 진단 워크플로

§ 11.3-11.4를 실전에서 어떻게 결합하는가:

1. 기본 Cox 모형 후보 작성
   ↓
2. 각 연속 공변량 Z_k에 대해:
   a) Z_k를 빼고 적합
   b) 마팅게일 잔차 vs Z_k LOWESS
   c) 형태에 따라 변환 (선형 / 로그 / 이산화 / 다항)
   ↓
3. 변환된 함수 형태로 모형 재적합
   ↓
4. 각 공변량의 PH 가정 검토:
   a) log-cumulative + 차이 도표 (이산 공변량)
   b) Andersen plot (보충)
   c) Arjas plot (모형 재적합 없이 층화 변경)
   d) Score process (연속 공변량 자연 처리)
   ↓
5. PH 위반 발견 시:
   - 시간의존 공변량 추가 (§ 9.2)
   - 위반 공변량으로 층화 (§ 8.6)
   - 가산 모형으로 전환 (Ch.10)

두 절의 핵심 통찰

§ 11.3 마팅게일 잔차는 올바른 함수 형태를 찾고, § 11.4 그래프 도구는 PH 가정의 타당성을 검증한다. 두 단계를 거친 모형만이 § 11.5-11.6의 이상치/영향력 분석을 의미 있게 진행할 수 있다.

5 코드 예시

5.1 Step 1 — 마팅게일 잔차로 함수 형태 진단

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi
import statsmodels.nonparametric.smoothers_lowess as sl

df = load_rossi()

# Z_1 = age 의 함수 형태 진단 — age를 제외하고 적합
df_no_age = df.drop(columns=["age"])
cph_no_age = CoxPHFitter()
cph_no_age.fit(df_no_age, duration_col="week", event_col="arrest")

# 마팅게일 잔차 (식 11.3.2)
mart = cph_no_age.compute_residuals(df_no_age, kind="martingale").iloc[:, 0]

# LOWESS 평활
ages = df["age"].values
smoothed = sl.lowess(mart.values, ages, frac=0.5, return_sorted=True)

plt.figure(figsize=(8, 5))
plt.scatter(ages, mart, alpha=0.3, s=20, label="잔차")
plt.plot(smoothed[:, 0], smoothed[:, 1], "r-", lw=2, label="LOWESS")
plt.axhline(0, color="gray", ls="--", alpha=0.5)
plt.xlabel("Age (years)"); plt.ylabel(r"Martingale residual $\hat{M}_j$")
plt.title(r"함수 형태 진단 — age vs $\hat{M}_j$")
plt.legend(); plt.grid(alpha=0.3)
print("LOWESS 형태로 함수 변환 결정:")
print("  직선 → 선형 그대로 / 곡선 → 변환 / 단조 변환점 → 이산화")

5.2 Step 2 — Schoenfeld 잔차 + Score Process로 PH 검정

from lifelines.statistics import proportional_hazard_test

cph = CoxPHFitter()
cph.fit(df, duration_col="week", event_col="arrest")

# Grambsch-Therneau PH 검정 (Schoenfeld 잔차 기반)
ph_test = proportional_hazard_test(cph, df, time_transform="rank")
print(ph_test)
# 각 공변량의 PH 가정 p-value 출력

# 시각적 검정 — Schoenfeld 잔차 vs 시간
sch = cph.compute_residuals(df, kind="schoenfeld")
times = df.loc[df["arrest"] == 1, "week"].sort_values().values

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
for i, col in enumerate(sch.columns[:6]):
    ax = axes[i // 3, i % 3]
    sch_smoothed = sl.lowess(sch[col].values, times, frac=0.7)
    ax.scatter(times, sch[col], alpha=0.3, s=15)
    ax.plot(sch_smoothed[:, 0], sch_smoothed[:, 1], "r-", lw=2)
    ax.axhline(0, color="gray", ls="--", alpha=0.5)
    ax.set_title(f"Schoenfeld: {col}")
    ax.set_xlabel("week"); ax.set_ylabel("residual")
plt.tight_layout()

5.3 Step 3 — log-cumulative + Andersen plot (이항 공변량)

from lifelines import KaplanMeierFitter, NelsonAalenFitter

# fin = 0 vs fin = 1 두 그룹의 PH 검정
group0 = df[df["fin"] == 0]; group1 = df[df["fin"] == 1]

naf0 = NelsonAalenFitter().fit(group0["week"], group0["arrest"])
naf1 = NelsonAalenFitter().fit(group1["week"], group1["arrest"])

# log-cumulative 도표
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

axes[0].plot(naf0.timeline, np.log(naf0.cumulative_hazard_.values + 1e-10),
             label="fin=0")
axes[0].plot(naf1.timeline, np.log(naf1.cumulative_hazard_.values + 1e-10),
             label="fin=1", ls="--")
axes[0].set_title("log-cumulative hazard")
axes[0].set_xlabel("time"); axes[0].set_ylabel(r"$\log \hat{H}_{g0}(t)$")
axes[0].legend()

# 차이 도표 — 공통 시간 그리드 필요
common_t = np.intersect1d(naf0.timeline, naf1.timeline)
H0_common = np.interp(common_t, naf0.timeline,
                       naf0.cumulative_hazard_.values.flatten())
H1_common = np.interp(common_t, naf1.timeline,
                       naf1.cumulative_hazard_.values.flatten())
diff = np.log(H1_common + 1e-10) - np.log(H0_common + 1e-10)
axes[1].plot(common_t, diff)
axes[1].axhline(diff[len(diff)//2], color="r", ls="--", label="평균")
axes[1].set_title("차이 도표 (PH 하 수평선)"); axes[1].legend()

# Andersen plot
axes[2].plot(H0_common, H1_common)
axes[2].plot([0, H0_common.max()], [0, H0_common.max()],
             "r--", label="PH 기준선")
axes[2].set_title("Andersen plot")
axes[2].set_xlabel(r"$\hat{H}_{0}(t)$"); axes[2].set_ylabel(r"$\hat{H}_{1}(t)$")
axes[2].legend()
plt.tight_layout()

5.4 Step 4 — R survival 패키지 (참고)

library(survival)

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

# § 11.3 함수 형태 — null 모형 마팅게일 잔차
null_fit <- coxph(Surv(time, status) ~ 1, data = lung)
mart <- residuals(null_fit, type = "martingale")
plot(lung$age, mart, xlab = "Age", ylab = "Martingale residual")
lines(lowess(lung$age, mart), col = "red", lwd = 2)
abline(h = 0, col = "gray", lty = 2)

# § 11.4 PH 검정 — Grambsch-Therneau (Schoenfeld 잔차)
ph_test <- cox.zph(fit)
print(ph_test)         # 검정 결과
plot(ph_test)          # 표준화된 Schoenfeld 잔차 도표

# § 11.4 log-cumulative 도표 (이산 공변량 sex)
fit_strata <- coxph(Surv(time, status) ~ age + ph.ecog + strata(sex), data = lung)
plot(survfit(fit_strata), fun = "cloglog",
     xlab = "log(time)", ylab = "log(-log S(t))")

6 핵심 요약

§ 11.3–11.4 한 줄 요약

마팅게일 잔차 $\hat{M}_j = \delta_j - r_j$ 의 LOWESS 평활은 공변량의 함수 형태를 알려주며, log-cumulative · Andersen · Arjas · score(Schoenfeld) 다섯 그래프는 PH 가정을 검정한다. 함수 형태 결정 → PH 검정 순서가 권고된다.

절	도구	핵심 식	주 용도
11.3	Martingale 잔차	$\hat{M}_j = \delta_j - r_j$ (식 11.3.2)	함수 형태 LOWESS 진단
11.4	log-cumulative	$\log \hat{H}_{g0}(t)$ vs $t$	평행성 검토
11.4	차이 도표	$\log \hat{H}_{g0} - \log \hat{H}_{10}$	수평성 검토
11.4	Andersen plot	$\hat{H}_{g0}(t)$ vs $\hat{H}_{10}(t)$	직선성 검토
11.4	Arjas plot	$N_g(t)$ vs $\text{TOT}_g(t)$ (식 11.4.1-2)	45도선 검토
11.4	Score / Schoenfeld	$U_k(t) = \sum (Z_{jk} - \bar{Z}_k)$ (식 11.4.7)	연속 공변량 PH

7 관련 주제

선행 지식

Ch.11 시리즈

관련 개념

Ch.10 — 가산 위험 모형 — PH 위반 시 대안
Modern Survival Analysis (RSF / DeepSurv)