생존 모형 비교¶

(1) 해석적으로. \(T\)의 참 누적위험을 \(H\)라 하면 \(S(t) = e^{-H(t)}\)이고 \(S(T) \sim \text{Unif}(0,1)\)이다(연속분포의 확률적분변환). 따라서

이다. \(U \sim \text{Unif}(0,1)\)일 때 \(-\ln U\)가 지수\((1)\)이기 때문이다. 모형이 옳아 \(\hat H \approx H\)이면 잔차 \(r_i = \hat H(t_i)\)가 지수\((1)\) 표본처럼 보여야 한다.

지수\((1)\)의 누적위험은 \(H_r(r) = -\ln e^{-r} = r\)로 항등함수다. 그러므로 잔차를 하나의 자료로 보고 그 누적위험을 추정해 잔차 자체에 대해 그리면 45도선이 나와야 한다.

코드의 세로축을 보자. 지수\((1)\)의 분포함수는 \(F(x) = 1-e^{-x}\)이므로 \(p\)분위수는

다. 코드가 \(i\)번째 점에 주는 값 \(-\ln(1 - i/(n+1))\)은 바로 \(p = i/(n+1)\)에서의 지수\((1)\) 분위수다. 그러므로 이 그림은 잔차의 순서통계량을 지수\((1)\) 분위수에 대해 그린 Q-Q 그림이고, 두 분포가 같으면 점들이 \(y = x\) 위에 놓인다. 45도선이 기준인 까닭이 이것이다.

분모를 \(n\)이 아니라 \(n+1\)로 둔 것도 여기서 설명된다. \(i = n\)에서 \(-\ln(1 - n/n) = \infty\)로 발산하기 때문이며, \(i/(n+1)\)은 널리 쓰이는 작도위치 가운데 하나다.

(2) 수치적으로. 잘못 지정된 경우로 잔차가 \(\text{Exp}(1)^{1.5}\)를 따르는 상황을 쓴다. 이는 참 분포가 와이불인데 지수를 적합했을 때 생기는 꼴이다.

import numpy as np
import matplotlib.pyplot as plt

def cox_snell_diagnostic(residuals):
    """Cox-Snell 잔차 진단 그림.

    모형이 옳다면 잔차의 누적위험함수가 단위 지수분포를 따른다. 곧
    H(r) = r 이므로 점들이 45도선 위에 놓인다. 생존분석에서 모형의
    적합도를 보는 표준적인 그림이다.
    """
    sorted_r = np.sort(residuals)
    n = len(sorted_r)

    # 잔차의 누적위험함수를 Nelson-Aalen 방식으로 추정한다.
    # 분모에 n+1 을 쓰는 것은 마지막 점에서 로그가 발산하는 것을 피하기 위함이다.
    H_na = -np.log(1 - np.arange(1, n + 1) / (n + 1))

    fig, ax = plt.subplots(figsize=(6, 6))
    ax.plot(sorted_r, H_na, "o", markersize=3, label="Cox-Snell residuals")
    ax.plot([0, sorted_r.max()], [0, sorted_r.max()], "r--", label="45-degree line")
    ax.set_xlabel("Cox-Snell Residuals")
    ax.set_ylabel("Cumulative Hazard")
    ax.legend()
    plt.tight_layout()
    plt.show()

    return sorted_r, H_na


# 모형이 옳을 때: Cox-Snell 잔차는 단위 지수분포를 따르므로 45도선 위에 놓인다
rng = np.random.default_rng(0)
r, H = cox_snell_diagnostic(rng.exponential(1.0, 200))

# --- 45도선에서 얼마나 벗어나는지 수로 잰다 ---
def report(name, r, H):
    slope, intercept = np.polyfit(r, H, 1)
    print(f"{name}: 기울기 {slope:.4f}, 절편 {intercept:+.4f}, "
          f"최대 이탈 {np.abs(H - r).max():.4f}, 평균잔차 {r.mean():.4f}")

report("옳은 모형  ", r, H)

# 잘못 지정한 경우: 잔차가 Exp(1)^1.5 를 따른다(참이 와이불인데 지수를 적합한 꼴).
bad = rng.exponential(1.0, 200) ** 1.5
rb = np.sort(bad)
Hb = -np.log(1 - np.arange(1, 201) / 201)
report("잘못된 모형", rb, Hb)

출력:

옳은 모형  : 기울기 0.9232, 절편 -0.0553, 최대 이탈 0.7544, 평균잔차 1.1292
잘못된 모형: 기울기 0.3323, 절편 +0.4599, 최대 이탈 17.8726, 평균잔차 1.5866

그림은 모형이 옳을 때의 모습이다. 점들이 45도선을 따라 놓이고, 맞춘 직선의 기울기가 \(0.9232\), 절편이 \(-0.0553\)으로 이상적인 \(1\)과 \(0\)에 가깝다. 지수\((1)\)의 평균이 1인데 잔차의 평균이 \(1.1292\)인 것도 \(n = 200\)에서 표준오차가 \(1/\sqrt{200} = 0.071\)임을 생각하면 \(1.8\) 표준오차 안이다.

잘못 지정된 쪽은 기울기가 \(0.3323\)으로 떨어지고 최대 이탈이 \(17.87\)까지 벌어진다. 꼬리에서 잔차가 너무 크게 나오기 때문이며, 그림에서는 점들이 선 아래로 휘어 내려가는 모양으로 나타난다.

다만 옳은 모형에서도 최대 이탈이 \(0.7544\)다. 꼬리의 몇 점은 늘 크게 벌어지며, 이는 가장 큰 순서통계량의 분산이 크기 때문이지 모형의 결함이 아니다. 콕스-스넬 그림은 몸통의 추세로 읽고 꼬리의 한두 점으로 판정하지 말아야 한다.

(a) 비모수(카플란-마이어). 표본이 작고 공변량이 없으므로 분포 가정이 위험하다. 카플란-마이어 추정량이 모형에 의존하지 않는 생존곡선을 주고, 필요하면 로그순위 검정으로 집단을 비교할 수 있다.

(b) 준모수(콕스 모형). 목표가 공변량(처리, 나이, 병기)이 위험에 미치는 효과를 추정하는 것이다. 콕스 모형은 기저위험 분포를 지정하지 않고 여러 공변량을 다룬다.

(c) 모수 모형(예: 와이불). 관측 시간 범위를 넘어선 외삽에는 완전히 지정된 분포가 필요하다. 비모수 모형과 콕스 모형은 외삽할 수 없다.

(a) AIC \(= -2\hat{\ell} + 2p\)이고 BIC \(= -2\hat{\ell} + p\ln n\)이다 (\(\ln 150 = 5.011\)).

(b) AIC(815.6)와 BIC(821.6) 모두 와이불 모형을 선호한다. 와이불과 로그로지스틱의 차이는 AIC와 BIC 모두에서 \(2.8\)로 같다. 두 모형의 모수 개수가 같아 벌점이 소거되기 때문이며, 이 경우 AIC와 BIC의 순위는 반드시 일치한다.

(c) 가능도비 검정(지수 대 와이불)은

이다. \(H_0: k = 1\) 아래에서 \(\Lambda \sim \chi^2_1\)이고 \(\alpha = 0.05\)의 임계값은 3.84다. \(29.0 \gg 3.84\)이므로 지수 모형을 기각한다. 상수 위험 가정은 지지되지 않는다.

추정된 생존함수가 \(\hat{S}(t_i)\)인 적합된 모수 모형에서 대상 \(i\)의 콕스-스넬 잔차는

로, 관측 시점에서의 추정 누적위험이다. 절단된 관측치에서는 \(r_i\)도 절단된다.

핵심 성질: 모형이 옳게 지정되어 있으면 바탕 생존분포와 무관하게 콕스-스넬 잔차가 \(\text{Exp}(1)\) 분포를 따른다. \(T \sim F\)이면 \(-\ln S(T) \sim \text{Exp}(1)\)이기 때문이다.

이는 확률적분변환의 한 형태다. \(U = F(T) \sim \text{Uniform}(0,1)\)이므로 \(S(T) = 1 - U\) 역시 균등분포이고, \(-\ln(\text{Uniform}(0,1)) \sim \text{Exp}(1)\)이다.

진단: 잔차의 누적위험에 대한 넬슨-알렌 추정치를 계산하고 잔차 자체에 대해 그린다. 모형이 잘 맞으면 점들이 45도 직선 \(H(r) = r\)을 따라야 한다. 벗어나면 오지정을 가리킨다.

다만 이 진단에는 알려진 한계가 있다. 모수를 자료에서 추정했으므로 잔차가 정확히 \(\text{Exp}(1)\)을 따르지 않으며, 특히 꼬리에서 어긋난다. 큰 잔차 몇 개가 45도 직선에서 벗어나는 것은 흔한 일이므로 그것만으로 모형을 기각하지 말라.

(a) 통계적 효율은 추정의 정밀도를 뜻한다. 같은 표본 크기에서 분산이 작을수록(신뢰구간이 좁을수록) 효율적인 추정량이다.

(b) 옳게 지정된 모수 모형은 알려진 분포 형태를 이용해 모든 관측치에 걸쳐 "힘을 빌린다." 카플란-마이어 추정량은 분포 구조를 활용하지 않고 각 사건시간을 독립적으로 다룬다. 모수 모형의 작은 분산은 옳은 구조 가정에서 나오며, 이 가정이 추정해야 할 양의 실질적인 개수를 줄인다.

구체적으로, 사건이 100건인 자료에서 카플란-마이어는 100개의 조건부확률을 추정하지만 지수 모형은 모수 하나만 추정한다. 100개의 정보가 하나로 모이므로 표준오차가 \(1/\sqrt{100}\) 수준으로 작아진다.

(c) 모수 가정이 틀리면 효율의 이점이 단점이 된다. 잘못 지정된 모형은 비효율적일 뿐 아니라 편향된다. 틀린 생존함수로 수렴한다. 카플란-마이어 추정량은 참 분포와 무관하게 언제나 일치성을 갖는다. 즉 모수 모형의 효율 이득은 강건성을 대가로 얻는 것이다.

더 위험한 점은 표준오차가 이 편향을 전혀 알려 주지 않는다는 것이다. 잘못 지정된 모형도 좁은 신뢰구간을 보고하며, 그 구간이 참값을 포함하지 않을 뿐이다. 그림 진단이 필수적인 이유가 여기 있다.

1단계 --- 비모수적 탐색. 처리군과 대조군의 카플란-마이어 생존곡선을 계산하고 95% 신뢰띠와 함께 시각화한다. 로그순위 검정으로 예비 비교를 한다. 넬슨-알렌 누적위험 그림을 살펴 위험의 모양(일정, 단조, 비단조)을 평가한다.

2단계 --- 공변량 보정을 위한 콕스 모형. 처리 지시자, 나이, 성별, 병기를 공변량으로 갖는 콕스 비례위험 모형을 적합한다. 위험비와 신뢰구간을 추정한다. 쇤펠트 잔차로 비례위험 가정을 점검한다.

3단계 --- 예측을 위한 모수 모형. 1단계에서 파악한 위험 모양에 근거해 후보 모수 모형을 적합한다(예: 위험이 단조면 와이불). AIC로 비교한다. 가장 좋은 모수 모형으로 중앙 생존시간을 외삽하고 매끄러운 생존곡선 추정치를 만든다.

4단계 --- 통합과 보고. 모수적 생존곡선과 콕스 기반 브레슬로 생존 추정치를 카플란-마이어 곡선 위에 겹쳐 그린다. 세 방법이 일치하면 결론에 대한 신뢰가 커진다. 어긋나면 오지정 가능성을 알리는 신호다. 콕스 모형의 처리 위험비를 주된 결과로, 카플란-마이어 곡선을 모형에 의존하지 않는 기준선으로, 모수 모형을 매끄러운 추정과 예측용으로 보고한다.