콕스 비례위험 실습¶

개요¶

콕스 비례위험 모형은 생존분석에서 가장 널리 쓰이는 회귀 틀이다. 기저위험의 모수적 형태를 지정하지 않고 공변량 효과를 위험함수에 연결하므로 준모수적 접근이다. 이 절에서는 모형을 정식화하고, 부분가능도를 유도하며, 위험비 해석을 논의하고, 비례위험 가정의 진단을 다룬다.

모형 정식화¶

콕스 모형은 공변량 벡터 \(\mathbf{x}_i = (x_{i1}, \ldots, x_{ip})^\top\)를 갖는 대상 \(i\)의 위험을 다음과 같이 설정한다.

여기서,

• \(h_0(t)\)는 기저위험으로, 전혀 지정되지 않은 임의의 음이 아닌 함수다.
• \(\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)^\top\)는 추정할 회귀계수다.
• \(\exp(\boldsymbol{\beta}^\top \mathbf{x}_i)\)는 대상 \(i\)의 상대위험 배수다.

비례위험 성질이 곧바로 따라온다. 임의의 두 대상 사이의 위험비가 시간에 걸쳐 일정하다.

비에서 기저위험 \(h_0(t)\)가 소거된다.

부분가능도¶

\(t_{(1)} < \cdots < t_{(K)}\)를 서로 다른 \(K\)개의 사건시간을 크기순으로 나열한 것이라 하고, \(t_{(j)}\)에서 사건을 겪는 대상을 \(i_j\)라 하자. \(t_{(j)}\)의 위험집합은

이다. 부분가능도는

이다. 기저위험 \(h_0(t_{(j)})\)이 분자와 분모에 모두 나타나 소거된다. 이것이 Cox(1972)의 핵심 착상이다. \(h_0(t)\)를 몰라도 공변량 효과를 추정할 수 있다.

부분 로그가능도는

이다.

추정¶

MLE \(\hat{\boldsymbol{\beta}}\)는 뉴턴-랩슨 반복으로 \(\ell_P\)를 최대화하여 얻는다.

여기서 \(U(\boldsymbol{\beta})\)는 점수벡터이고 \(\mathcal{I}(\boldsymbol{\beta})\)는 관측 정보행렬이다.

위험비¶

지수화한 계수 \(\exp(\hat{\beta}_j)\)는 다른 공변량을 고정한 채 공변량 \(x_j\)가 한 단위 증가할 때의 위험비다.

연속형 공변량이 \(c\) 단위 증가하면 위험비는 \(\exp(c \cdot \hat{\beta}_j)\)다.

신뢰구간¶

위험비의 \(100(1 - \alpha)\%\) 신뢰구간은

이다. 구간이 1을 포함하지 않으면 그 공변량 효과는 통계적으로 유의하다.

기저위험의 브레슬로 추정량¶

\(\hat{\boldsymbol{\beta}}\)를 추정한 뒤 기저 누적위험을 다음으로 추정한다.

대상별 생존함수는 그러면

이다.

각 집단 \(1{,}000\)명, 기저는 와이불 \(k = 1.5\)·척도 \(34\), 처리 효과는 \(\beta = -0.7\)(곧 \(\text{HR} = 0.497\))인 자료를 만들어 콕스 모형을 적합했다. 결과는 \(\hat\beta = -0.694\), \(\text{HR} = 0.500\)(95% 신뢰구간 \(0.440\)--\(0.568\))으로 참값을 거의 그대로 되찾는다. 기저위험의 모양을 한 번도 가정하지 않고 얻은 값이다.

왼쪽 그림이 이 절에서 가장 놀라운 부분이다. 보라색 계단이 브레슬로 추정치 \(\hat H_0(t)\)이고 회색 점선이 자료를 만든 참 \(H_0(t) = (t/34)^{1.5}\)인데, 둘이 거의 포개진다. \(t = 10\)에서 \(0.166\) 대 \(0.160\), \(t = 20\)에서 \(0.455\) 대 \(0.451\), \(t = 30\)에서 \(0.838\) 대 \(0.829\), \(t = 40\)에서 \(1.330\) 대 \(1.276\)이다. 부분가능도는 이 곡선을 계산에 한 번도 쓰지 않았다. \(\hat\beta\)를 먼저 구하고 나서, 각 사건시점의 사건 수를 그 시점 위험집합의 가중합으로 나누어 사후에 쌓아 올린 것이다.

이 순서가 콕스 모형의 설계다. 성가신 무한차원 모수 \(h_0(t)\)를 추정 단계에서는 소거해 버리고, 관심 모수 \(\boldsymbol\beta\)를 먼저 확정한 뒤, 필요하면 그때 가서 \(\boldsymbol\beta\)를 고정한 채 기저를 비모수적으로 채워 넣는다. 그래서 "기저를 지정하지 않는다"는 말이 "기저를 알 수 없다"는 뜻은 아니다.

오른쪽은 그렇게 얻은 \(\hat S_0(t)^{\exp(\hat\beta x)}\)를 집단별 카플란-마이어 위에 겹친 것이다. 굵은 선(콕스 예측)이 연한 선(카플란-마이어)을 거의 완벽히 따라간다. 모형이 옳게 지정되었을 때의 모습이며, 여기서 체계적으로 어긋난다면 비례위험이 깨졌거나 공변량의 함수 형태가 잘못된 것이다. 아래 일치도 지수는 \(0.585\)로, 두 집단만 구별하는 이항 공변량 하나로는 이 정도가 한계다. \(\text{HR} = 0.50\)처럼 강한 효과조차 개별 대상의 순서를 맞히는 능력으로는 \(0.5\)(동전 던지기)에서 크게 못 벗어난다는 사실은 21.5절에서 다시 다룬다.

비례위험 가정 점검하기¶

비례위험(PH) 가정은 콕스 모형의 타당성에 결정적이다. 위배되면 위험비가 시간에 걸쳐 일정하지 않고 표준적 해석이 무너진다.

그림 방법¶

• 로그-로그 그림: 집단마다 \(\ln(-\ln \hat{S}(t))\)를 \(\ln t\)에 대해 그린다. 비례위험 아래에서는 곡선들이 대략 평행해야 한다.
• 쇤펠트 잔차: 공변량마다 척도화된 쇤펠트 잔차를 시간에 대해 그린다. 기울기가 0이 아니면 시간 의존 효과를 시사한다.

형식적 검정¶

그램브시-서노 검정은 척도화된 쇤펠트 잔차를 시간의 함수에 회귀시킨다. 공변량 \(j\)의 기울기가 유의하면 \(\beta_j\)가 시간에 따라 변한다는 뜻이다.

구현 스케치¶

import numpy as np

def partial_log_likelihood(beta, X, times, events):
    """Cox 모형의 부분로그가능도를 계산한다.

    "부분"이라 부르는 까닭은 기저위험함수를 아예 셈에서 빼기 때문이다.
    사건이 일어난 시점마다 "그 순간 위험집합에 있던 사람들 중 하필 이
    사람에게 사건이 일어날 확률"만 곱해 나가면, 기저위험이 분자와 분모에서
    약분되어 사라진다. 그래서 위험함수의 모양을 가정하지 않고도 계수를
    추정할 수 있다.

    매개변수
    --------
    beta   : 길이 p 인 계수벡터
    X      : (n, p) 공변량 행렬
    times  : 관측된 시각
    events : 사건 지시자 (1 = 사건, 0 = 중도절단)

    돌려주는 값
    ----------
    ll : 부분로그가능도 값
    """
    # 선형예측자. exp 를 씌운 값이 그 사람의 상대적 위험이 된다.
    risk_scores = X @ beta
    exp_scores = np.exp(risk_scores)

    # 시각을 내림차순으로 정렬한다. 이러면 누적합이 곧 "그 시점 이후까지
    # 남아 있는 사람들", 곧 위험집합의 합이 된다.
    order = np.argsort(-times)
    sorted_events = events[order]
    sorted_exp = exp_scores[order]
    sorted_scores = risk_scores[order]

    # 위험집합의 합을 누적합 한 번으로 얻는다. 시점마다 집합을 다시
    # 만들면 O(n^2) 이 되는 계산이 O(n log n) 으로 끝난다.
    cumsum_exp = np.cumsum(sorted_exp)

    # 사건이 관측된 사람만 더한다. 중도절단된 사람은 위험집합에 기여할
    # 뿐 자기 항을 갖지 않는다.
    ll = np.sum(sorted_events * (sorted_scores - np.log(cumsum_exp)))
    return ll

# --- (1) 닫힌 꼴과 맞추어 본다 ---
times  = np.array([4., 7., 9., 12.])
events = np.array([1, 0, 1, 1])
X      = np.array([[1.], [0.], [0.], [1.]])

for b in (-1.0, -0.5, 0.0, 0.5, 1.0):
    closed = b - np.log(2) - 2 * np.log(1 + np.exp(b))
    print(f"beta = {b:+.1f}   코드 {partial_log_likelihood(np.array([b]), X, times, events):.6f}"
          f"   닫힌 꼴 {closed:.6f}")

grid = np.linspace(-3, 3, 60001)
vals = np.array([partial_log_likelihood(np.array([b]), X, times, events) for b in grid])
print(f"격자 최대점 beta = {grid[vals.argmax()]:.4f},  최댓값 {vals.max():.6f}")
print(f"-3 ln 2 = {-3 * np.log(2):.6f}")

# --- (2) 공변량을 통째로 옮겨도 값이 그대로여야 한다 ---
for c in (0.0, 5.0, -3.7):
    print(f"x 에 {c:+.1f} 를 더하면 ll(0.3) = "
          f"{partial_log_likelihood(np.array([0.3]), X + c, times, events):.9f}")

beta = -1.0   코드 -2.319671   닫힌 꼴 -2.319671
beta = -0.5   코드 -2.141301   닫힌 꼴 -2.141301
beta = +0.0   코드 -2.079442   닫힌 꼴 -2.079442
beta = +0.5   코드 -2.141301   닫힌 꼴 -2.141301
beta = +1.0   코드 -2.319671   닫힌 꼴 -2.319671
격자 최대점 beta = 0.0000,  최댓값 -2.079442
-3 ln 2 = -2.079442
x 에 +0.0 를 더하면 ll(0.3) = -2.101857669
x 에 +5.0 를 더하면 ll(0.3) = -2.101857669
x 에 -3.7 를 더하면 ll(0.3) = -2.101857669

해석¶

• 콕스 모형은 공변량이 위험에 미치는 상대적 효과를 추정하며 절대 위험 수준은 추정하지 않는다.
• 위험비는 순간 사건율의 곱셈적 변화를 정량화하며 누적확률의 변화가 아니다.
• 기저위험은 성가신 모수다. 절대 생존 예측이 필요하면 브레슬로 추정량으로 복원한다.
• 위험비를 시간에 걸쳐 일정한 효과로 해석하기 전에 언제나 비례위험 가정을 확인하라.

연습문제¶

정리하며¶

콕스 모형을 실제로 적합했다.

• lifelines 의 CoxPHFitter 가 표준 도구다. 시간·사건·공변량 열을 담은 데이터프레임을 넘기면 계수와 위험비, 신뢰구간이 나온다.
• check_assumptions() 를 반드시 호출한다. 쇤펠트 기반 진단을 자동으로 수행하며, 위반이 있으면 처방까지 제안한다.
• 범주형 변수는 더미로 바꾼다. 기준 범주가 무엇인지에 따라 위험비의 해석이 달라진다.
• 기저 생존함수를 따로 얻을 수 있다. baseline_survival_ 로 개별 예측 곡선을 그릴 수 있으며, 부분가능도가 \(h_0\) 를 쓰지 않았어도 사후에 추정된다.
• 일치도 지수(C-index)로 예측력을 본다. 순위 기반 지표이며 AUC 의 생존분석판이다.

다음 절 모형 비교로 21장을 마무리한다.