Klein § 11.1–11.2 — Introduction & Cox-Snell Residuals

Cox 모형 진단의 4가지 측면 + Cox-Snell 잔차로 전반적 적합도 검증하기

Cox 비례위험 모형 진단의 첫걸음. § 11.1에서 회귀 진단의 4가지 측면(함수 형태, PH 가정, 개별 정확도, 영향력)을 정리하고, § 11.2에서 가장 처음 제안된 잔차인 Cox-Snell 잔차로 모형의 전반적 적합도를 검증하는 방법을 다룬다. 확률 적분 변환 원리, Nelson-Aalen 점검 도표, BMT 데이터로의 실전 적용까지 포함. (Klein & Moeschberger, 2003, § 11.1–11.2)

Statistics
Survival Analysis
저자

Kwangmin Kim

공개

2026년 04월 29일

1 도입 — 왜 회귀 진단이 필요한가

Ch.8-9에서는 Cox 비례위험(PH) 모형의 추정 절차와 가설 검정을 다루었다. 그러나 그 과정에서 한 가지 결정적 가정을 거의 검토하지 않았다: 모형 자체가 옳은가?

선형 회귀에서는 잔차 도표(잔차 vs 적합값, Q-Q 플롯, 영향력 통계 등)가 거의 자동화된 진단 도구로 사용된다. Cox 모형에서도 비슷한 절차가 필요하지만, 잔차의 정의 자체가 단순하지 않다. 그 이유는:

  1. 반응 변수가 우중도절단된 시간 — 모든 개체가 사건을 경험하지 않음
  2. 베이스라인 위험 \(h_0(t)\) 가 비모수 — 잔차에서 차감해야 할 “예측값”이 비모수 추정치
  3. 부분우도(partial likelihood) 기반 추정 — 가능도가 표준 회귀의 잔차 제곱합 형태가 아님

이런 복잡성 때문에 Cox 모형용 잔차는 여러 종류가 제안되었으며, 각 잔차는 모형의 서로 다른 측면을 진단한다.

2 § 11.1 — 진단의 4가지 측면

Klein은 Cox 모형 진단을 다음 네 가지 질문으로 정리한다.

2.1 측면 1 — 공변량의 함수 형태 (Functional Form)

공변량 \(Z\) 의 효과를 어떤 함수로 모형화해야 하는가?

연속 공변량 \(Z\) 의 영향이 다음 중 어떤 형태로 가장 잘 표현되는가:

형태 모형 의미
선형 \(h_0(t) \exp(\beta Z)\) \(Z\) 가 1 단위 증가하면 위험이 \(e^\beta\) 배
로그 \(h_0(t) \exp(\beta \log Z)\) \(Z\) 가 두 배 되면 위험이 \(2^\beta\) 배
이산화 \(h_0(t) \exp(\beta \mathbb{1}\{Z \geq Z_0\})\) 임계점 \(Z_0\) 기준 이분화
다항 \(h_0(t) \exp(\beta_1 Z + \beta_2 Z^2)\) U-shape 가능
왜 함수 형태가 중요한가

선형 가정 \(\exp(\beta Z)\) 는 “단위 증가 효과가 어디서나 동일”을 의미한다. 그러나 의학 데이터에서는:

  • 나이: 50세 이전 효과 약함, 50세 이후 급증 → 이산화 또는 비선형
  • 약물 농도: 너무 낮으면 효과 없음, 너무 높으면 부작용 → U-shape
  • 백혈구 수: 정상 범위 효과 약함, 정상 이탈 시 위험 → 이산화

함수 형태가 잘못되면 효과가 통째로 사라지거나 반대로 측정된다.

이 측면은 § 11.3 마팅게일 잔차로 진단한다.

2.2 측면 2 — PH 가정의 타당성

위험비가 정말 시간에 무관하게 일정한가?

Cox 모형의 핵심은 다음 식이다:

\[ \frac{h(t \mid Z = 1)}{h(t \mid Z = 0)} = \exp(\beta) \quad \text{(시간 } t \text{ 에 의존하지 않음)} \]

이 가정이 깨지면(예: 초기에는 \(Z = 1\) 이 위험하지만 후기에는 \(Z = 0\) 이 위험) 단일 \(\beta\) 추정값은 의미를 잃는다. § 9.2에서 시간의존 공변량 \(Z(t) = Z \cdot \log t\) 로 검정하는 방법을 다뤘지만, 그래프적 검정이 추가 통찰을 준다.

이 측면은 § 11.4 score / Schoenfeld 잔차 + Andersen / Arjas / log-cumulative 도표로 진단한다.

2.3 측면 3 — 개별 관측치의 적합도 (Outliers)

모형이 예측한 것보다 너무 빨리 또는 너무 늦게 사망한 환자가 있는가?

이는 통상 회귀의 이상치 탐지와 같은 개념이다. 위험 점수가 낮은(좋은 예후) 환자가 빨리 사망하면 그 관측치는 모형에 잘 맞지 않는다.

이 측면은 § 11.5 deviance 잔차로 진단한다.

2.4 측면 4 — 영향력 (Influence / Leverage)

한 관측치가 추정 결과를 과도하게 좌우하는가?

작은 표본에서는 단 한 개의 극단 관측치가 \(b\) 추정값을 크게 흔들 수 있다. 이런 관측치를 식별하면 모형의 안정성을 평가할 수 있다.

이 측면은 § 11.6 dfbeta 잔차로 진단한다.

2.5 다섯 잔차 한눈에

잔차 기호 진단 측면 해당 절
Cox-Snell \(r_j\) 전반적 적합도 § 11.2
Martingale \(\hat{M}_j\) 함수 형태 § 11.3
Schoenfeld / Score \(S_{jk}\) PH 가정 § 11.4
Deviance \(D_j\) 이상치 § 11.5
dfbeta \(\Delta_{jk}\) 영향력 § 11.6
잔차의 위계 — 마팅게일이 중심

다섯 잔차는 독립적 발명이 아니라 마팅게일 잔차의 변형이다:

                  Martingale 잔차 M̂_j = δ_j − r_j
                          │
        ┌─────────────────┼─────────────────┐
        ▼                 ▼                 ▼
  Cox-Snell r_j     Score S_jk        Deviance D_j
  (= δ_j − M̂_j)    (편차 적분)       (정규형 변환)
                         │
                         ▼
                  dfbeta Δ_jk
                  (= I⁻¹ S_jk)

따라서 Cox-Snell과 마팅게일은 사실상 같은 정보의 두 표현이다. § 11.2에서 보듯, \(r_j = \delta_j - \hat{M}_j\) 가 정확히 성립한다.

3 § 11.2 — Cox-Snell 잔차

3.1 역사적 맥락

Cox-Snell 잔차는 Cox & Snell (1968) 이 일반화된 잔차 개념으로 처음 제안한 것이다. 흥미롭게도 1968년 시점에는 아직 마팅게일 이론이 생존 분석에 도입되지 않았으며, 이 잔차는 통상 회귀의 잔차 개념을 일반화한 것에서 출발했다.

3.2 정의 — 식 (11.2.1)

데이터 \((T_j, \delta_j, Z_j)\), \(j = 1, \ldots, n\) 에 Cox 모형을 적합하여 \(b = (b_1, \ldots, b_p)^\top\) 와 Breslow 베이스라인 누적 위험 \(\hat{H}_0(t)\) 를 얻었다고 하자.

Cox-Snell 잔차:

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

이는 곧 개체 \(j\) 에 대한 추정된 누적 위험:

\[ r_j = \hat{H}(T_j \mid Z_j) \]

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

왜 Cox-Snell 잔차가 적합도를 측정하는가

확률론의 기초 정리부터 시작한다.

확률 적분 변환 (Probability Integral Transform):

임의의 연속 확률변수 \(X\) 가 누적분포함수 \(F(x)\) 를 가진다면:

\[ U = F(X) \sim \text{Uniform}(0, 1) \]

이를 생존 함수와 위험 함수로 바꾸면:

\[ S(X) = 1 - F(X) \sim \text{Uniform}(0, 1) \]

\[ H(X) = -\log S(X) \sim \text{Exp}(1) \quad \text{(단위 지수분포)} \]

마지막 등식이 핵심이다: 누적 위험에 자기 자신의 진짜 시간을 대입하면 단위 지수분포가 된다.

이제 Cox 모형이 정확하다고 가정한다. 그러면:

\[ H(T_j \mid Z_j) \sim \text{Exp}(1) \]

만약 \(\hat{H}(T_j \mid Z_j) \approx H(T_j \mid Z_j)\) 라면, \(r_j\) 는 단위 지수분포의 표본처럼 보여야 한다. 우중도절단이 있는 경우에는 단위 지수분포의 우중도절단 표본이 된다 (\(\delta_j = 1\) 이면 사건 시점, \(\delta_j = 0\) 이면 절단 시점).

3.4 점검 도표 — 45도 직선

\(r_j\) 가 단위 지수분포의 표본인지 확인하려면:

  1. \((r_j, \delta_j)\) 를 새로운 생존 데이터로 간주
  2. Nelson-Aalen 추정량 \(\hat{H}_r(r)\) 을 § 4.2 식 (4.2.3)으로 계산
  3. 단위 지수분포의 누적 위험은 \(H_E(r) = r\) 임을 이용
핵심 검증 식

단위 지수분포라면 \(H(r) = -\log e^{-r} = r\). 따라서:

\[ \hat{H}_r(r_j) \approx r_j \quad (\text{모든 } j) \]

도표 \(\hat{H}_r(r_j)\) vs \(r_j\) 를 그렸을 때 원점을 지나는 기울기 1 직선 (45도선)이 나오면 모형이 적합하다.

직선에서 벗어나는 패턴:

도표 형태 해석
45도선 적합
위로 휘어짐 모형이 위험을 과소예측 (실제 사망이 더 빠름)
아래로 휘어짐 모형이 위험을 과대예측
곡선 / 비선형 함수 형태 또는 PH 위반

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

시간의존 공변량 \(Z_j(t)\) 또는 층화 베이스라인 \(\hat{H}_{0j}(t)\) 가 있는 경우:

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

핵심 변경점:

  • 층화: 개체가 속한 층의 베이스라인 사용
  • 시간의존: 사건 시점에서의 공변량 값 사용

논리는 동일하다 — \(r_j\) 가 단위 지수처럼 보이면 모형이 적합하다.

3.6 Klein Example 11.1 — BMT 데이터 진단

데이터: 137명의 골수이식 환자, 공변량 \(Z_1, \ldots, Z_8\):

  • \(Z_1\): 환자 나이 (28세 중심화)
  • \(Z_2\): 기증자 나이 (28세 중심화)
  • \(Z_3\): \(Z_1 \times Z_2\) 교호작용
  • \(Z_4\): FAB 등급 (4 또는 5 + 급성 백혈병)
  • \(Z_5\): 급성 백혈병 (분류)
  • \(Z_6\): FAB 4-5 + 급성 백혈병
  • \(Z_7\): 진단-이식 대기 시간 (9개월 중심화)
  • \(Z_8\): MTX 사용 (GVHD 예방)

3.6.1 1단계 — 전체 도표 (Figure 11.1)

8개 공변량 모두를 포함한 Cox 모형 적합 후 식 (11.2.1)로 \(r_j\) 계산. \((r_j, \delta_j)\) 의 Nelson-Aalen 추정치를 그렸을 때 도표는 대체로 45도선에 가까웠다. 즉 모형이 “그리 나쁘지 않게” 적합한다.

3.6.2 2단계 — MTX 그룹별 분리 도표 (Figure 11.2)

그러나 환자를 MTX 사용 여부로 나누어 두 그룹의 Nelson-Aalen 추정치를 별도로 그리니, MTX 그룹의 도표가 45도선 위쪽으로 벗어났다. 이는:

  • 모형이 MTX 환자의 위험을 과소예측
  • \(Z_8\) 의 효과가 상수가 아닐 가능성 (PH 위반)

3.6.3 3단계 — 형식 검정으로 확인

PH 가정의 형식 검정을 위해 시간의존 공변량 \(Z_9(t) = Z_8 \log t\) 를 추가:

\[ h(t \mid Z) = h_0(t) \exp(\beta_1 Z_1 + \cdots + \beta_8 Z_8 + \beta_9 Z_8 \log t) \]

Wald 검정 결과 \(H_0: \beta_9 = 0\) 의 p-value = 0.0320. 따라서 MTX의 효과가 시간에 따라 변한다는 증거가 있다.

3.6.4 4단계 — 층화 모형으로 재적합 (Figure 11.3)

\(Z_8\) 을 모형에서 제외하고 MTX로 층화 한 모형을 적합. 식 (11.2.2)로 잔차 재계산. 이번에는 두 그룹 모두 45도선에 잘 맞는다.

Cox-Snell 진단의 실전 가치

이 예제가 보여주는 진단의 흐름:

  1. 전체 도표가 OK여도 그룹별 분리 도표를 봐야 한다 — 한 변수의 PH 위반이 전체 평균에서는 가려질 수 있다
  2. 잔차 도표는 가설을 시사, 형식 검정이 결론을 짓는다 — 도표만으로 PH 위반을 단정하지 않음
  3. 층화는 PH 위반의 첫 번째 처방 — 위반 변수를 모형에서 빼고 층화 변수로 변환

3단계 워크플로(전체 → 그룹별 → 형식 검정 → 층화)를 기억하면 실전에서 즉시 사용 가능하다.

3.7 한계와 주의 사항

Cox-Snell 잔차의 5가지 한계

Klein이 명시적으로 지적하는 사항들:

1. 위반 종류를 알려주지 않는다 (Crowley & Storer, 1983)

도표가 45도선에서 벗어나면 “어디선가 잘못됐다”는 신호이지만, 그것이: - 함수 형태 잘못인지 - PH 위반인지 - 다른 요인인지

알려주지 않는다. 따라서 다른 잔차로 후속 분석이 필요하다.

2. 단위 지수 근사가 작은 표본에서 부정확 (Lagakos, 1981)

이론은 \(\beta\) 와 \(H_0\) 의 참값을 사용할 때 성립한다. 추정값을 사용하면 추정 오차가 잔차에 침투하며, 작은 표본에서는 이 효과가 크다.

3. 꼬리 부분의 변동성

\(\hat{H}_r\) 의 분산은 시간에 따라 증가한다. 도표의 오른쪽 끝(큰 \(r_j\))에서는 변동이 커서 직선에서 벗어나도 문제로 단정하기 어렵다.

4. 누락 공변량 진단에는 마팅게일이 우월 (Kay, 1984)

누락 공변량 \(Z_*\) 를 빼고 적합한 모형의 Cox-Snell vs \(Z_*\) 도표보다, 마팅게일 잔차 vs \(Z_*\) 도표가 더 직관적이다 (마팅게일은 0 중심).

5. 우중도절단 처리의 모호성

원래 Cox & Snell의 정의는 절단된 잔차에 1을 더하는 보정을 제시했지만, 현대 구현은 보정 없이 사용한다. 이로 인해 절단 비율이 높을 때 도표 해석이 까다롭다.

4 마팅게일 잔차와의 관계 — 미리보기

§ 11.3에서 자세히 다룰 마팅게일 잔차는 다음과 같이 정의된다:

\[ \hat{M}_j = \delta_j - r_j \tag{식 11.3.2} \]

즉 마팅게일 = 사건 지시자 − Cox-Snell. 두 잔차는 동일한 정보의 다른 표현이지만 용도가 다르다:

잔차 분포 중심 주 용도
Cox-Snell \(r_j\) 단위 지수 (절단) 1 전반 적합도 (45도 도표)
Martingale \(\hat{M}_j\) 비대칭, 0 중심 0 함수 형태 (산점도 + LOWESS)

Cox-Snell은 “전체가 단위 지수에 맞는가”를 보고, 마팅게일은 “특정 공변량과 어떤 관계가 있는가”를 본다.

5 코드 예시

5.1 Step 1 — 순수 Python으로 Cox-Snell 잔차 구현

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def fit_cox_simple(T, delta, Z, max_iter=50, tol=1e-6):
    """단일 공변량 Cox 모형 — Newton-Raphson"""
    n = len(T)
    order = np.argsort(-T)  # 큰 시간부터
    T_s, d_s, Z_s = T[order], delta[order], Z[order]

    beta = 0.0
    for _ in range(max_iter):
        # 위험집합: 각 사건 시점에서 T_j ≥ T_(i) 인 개체들
        # 큰 시간부터 정렬했으므로, 누적합으로 위험집합 계산 (역순)
        eta = beta * Z_s
        r = np.exp(eta)

        # 사건 시점만 골라 score / information 계산
        score, info = 0.0, 0.0
        # 시간 오름차순 처리를 위해 다시 뒤집음
        T_a = T_s[::-1]; d_a = d_s[::-1]; Z_a = Z_s[::-1]; r_a = r[::-1]

        cum_r = 0.0; cum_rZ = 0.0; cum_rZ2 = 0.0
        # 시간 내림차순으로 누적
        for i in range(len(T_a) - 1, -1, -1):
            cum_r += r_a[i]
            cum_rZ += r_a[i] * Z_a[i]
            cum_rZ2 += r_a[i] * Z_a[i] ** 2
            if d_a[i] == 1:
                Z_bar = cum_rZ / cum_r
                score += Z_a[i] - Z_bar
                info += cum_rZ2 / cum_r - Z_bar ** 2

        if abs(score) < tol:
            break
        beta += score / info
    return beta

def breslow_baseline_H(T, delta, Z, beta):
    """Breslow 베이스라인 누적 위험 식 (8.8.2)"""
    n = len(T)
    order = np.argsort(T)
    T_s, d_s, Z_s = T[order], delta[order], Z[order]
    risk = np.exp(beta * Z_s)

    H0 = np.zeros(n)
    cum = 0.0
    # 위험집합 합: T_j ≥ T_(i) 의 risk 합
    for i in range(n):
        if d_s[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

def cox_snell_residuals(T, delta, Z, beta, H0_at_T):
    """Cox-Snell 잔차 식 (11.2.1)"""
    return H0_at_T * np.exp(beta * Z)

def nelson_aalen(times, events):
    """Nelson-Aalen 추정량 식 (4.2.3)"""
    order = np.argsort(times)
    t_s, e_s = times[order], events[order]
    n = len(t_s)
    H = np.zeros(n)
    cum = 0.0
    for i in range(n):
        if e_s[i] == 1:
            risk_set_size = n - i
            cum += 1.0 / risk_set_size
        H[i] = cum
    return t_s, H

# 합성 데이터: 진짜 모형이 Cox PH인 경우
n = 200
Z = np.random.binomial(1, 0.5, n).astype(float)
true_beta = 1.5
T_true = np.random.exponential(1.0 / np.exp(true_beta * Z))
C = np.random.exponential(2.0, n)
T = np.minimum(T_true, C)
delta = (T_true <= C).astype(int)

beta_hat = fit_cox_simple(T, delta, Z)
H0 = breslow_baseline_H(T, delta, Z, beta_hat)
r = cox_snell_residuals(T, delta, Z, beta_hat, H0)

print(f"추정 β = {beta_hat:.3f} (참값 {true_beta})")
print(f"검열 비율 = {1 - delta.mean():.2%}")

# Cox-Snell 도표
r_sorted, H_r = nelson_aalen(r, delta)
plt.figure(figsize=(6, 6))
plt.plot(r_sorted, H_r, "b-", label=r"$\hat{H}_r(r)$")
plt.plot([0, r.max()], [0, r.max()], "r--", label="45° 기준선")
plt.xlabel("Cox-Snell 잔차 $r_j$"); plt.ylabel("누적 위험 추정")
plt.title("Cox-Snell 적합도 도표 (45도선에 맞으면 적합)")
plt.legend(); plt.grid(alpha=0.3)

5.2 Step 2 — lifelines로 Cox-Snell 잔차 계산

lifelines는 직접 Cox-Snell 잔차 함수를 제공하지 않지만, 마팅게일 잔차로부터 변환할 수 있다 (\(r_j = \delta_j - \hat{M}_j\)):

from lifelines import CoxPHFitter, NelsonAalenFitter
from lifelines.datasets import load_rossi
import matplotlib.pyplot as plt

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

# 마팅게일 잔차 → Cox-Snell 변환
mart = cph.compute_residuals(df, kind="martingale")
cs = df["arrest"] - mart.iloc[:, 0]  # r_j = δ_j − M̂_j

# Nelson-Aalen으로 H_r 계산
naf = NelsonAalenFitter()
naf.fit(durations=cs.abs(), event_observed=df["arrest"])

plt.figure(figsize=(6, 6))
naf.cumulative_hazard_.plot(ax=plt.gca(), color="b")
m = cs.abs().max()
plt.plot([0, m], [0, m], "r--", label="45° 기준선")
plt.xlabel("Cox-Snell 잔차"); plt.ylabel(r"$\hat{H}_r(r)$")
plt.title("Rossi 데이터 — Cox-Snell 적합도 도표")
plt.legend()

5.3 Step 3 — R survival 패키지 (참고)

library(survival)

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

# Cox-Snell 잔차: martingale로부터 변환
mart <- residuals(fit, type = "martingale")
cs <- lung$status - mart  # r_j = δ_j − M̂_j

# Nelson-Aalen 추정
fit_cs <- survfit(Surv(cs, lung$status) ~ 1)
plot(fit_cs$time, -log(fit_cs$surv),
     type = "s", xlab = "Cox-Snell 잔차", ylab = "추정 누적 위험")
abline(0, 1, col = "red", lty = 2)

6 핵심 요약

§ 11.1–11.2 한 줄 요약

Cox 모형 진단은 함수 형태 / PH 가정 / 이상치 / 영향력 네 측면을 점검하며, Cox-Snell 잔차 \(r_j = \hat{H}_0(T_j) e^{b^\top Z_j}\) 는 모형이 정확하면 단위 지수분포의 (절단된) 표본처럼 보인다는 확률 적분 변환 원리에 기초해 전반적 적합도를 검증한다 — 도표가 45도선에 맞으면 OK.

항목 핵심 내용
진단 4측면 함수 형태 · PH · 이상치 · 영향력
Cox-Snell 식 \(r_j = \hat{H}_0(T_j) \exp(\sum b_k Z_{jk})\)
이론적 분포 단위 지수 (절단) — 확률 적분 변환
점검 도표 \(\hat{H}_r(r_j)\) vs \(r_j\) → 45도선 비교
시간의존 확장 식 (11.2.2): 층별 / 시간의존 베이스라인 사용
마팅게일 관계 \(r_j = \delta_j - \hat{M}_j\)
주 한계 위반 종류 미식별 · 작은 표본에서 부정확 · 꼬리 변동성

7 관련 주제

선행 지식

Ch.11 시리즈

관련 개념

Subscribe

Enjoy this blog? Get notified of new posts by email: