연금술사

Cox proportional hazard model

각 i번째 환자(데이터포인트) X가 있다고 하자. $X_{i}=(X_{i1}, ...X_{ip})$ 와 같이 표기할 수 있으며 i번쨰 환자의 각 p라는 변수(공변량, covariates)를 의미한다. Cox비례위험모형은 아래와 같이 작성할 수 있다[식(2)]. 

$\lambda(t|X_{i}) = \lambda_{0}(t)exp(\beta_{1}X_{i1}+...+\beta_{p}X_{ip}$ [식2]

이를 Matrix 표기($X_{i}$)로 변경하면 흔히 보는 식이다[식3]. 보통 이를 hazard function이라고 한다.

$\lambda(t|X_{i}) = \lambda_{0}(t)exp(X_{i}\beta)$ [식3]

위의 식을 들여다보면, X가 커지든 작아지든 t보다는 $\beta$와 연관이 있다.

1) 즉, 시점이 고정이라면(환자들끼리 같은 조건이라면), X가 커질수록 $\beta$가 커져서, 위험이 올라간다고 할 수 있다. 즉 시간에 따라서는 위험에 대한 비율이 일정하게 유지가 된다(비례위험ㄱ자ㅓㅇ)

2) 또한, t에 영향을 받는것은 $\lambda_{0}$만 영향을 받기때문에, 시간이 증가할수록 커진다.

그리고, 흔히 위험비(harzard ratio, HR)라고 불리는 위험에 대한 비율(집단A에 대한 집단 B에 대한 이벤트 발생에 대한 비율)을 다음과 같이 구할 수 있다. 위의 hazard function을 A집단에 있는 환자a와 B집단에 있는 환자b에 대해서 적용하고자면 아래와 같이 각각 구할 수 있다.

$\lambda(t|X_{a}) = \lambda_{0}(t)exp(\beta_{1}X_{i1a}+...+\beta_{p}X_{ia})$

$\lambda(t|X_{b}) = \lambda_{0}(t)exp(\beta_{1}X_{i1a}+...+\beta_{p}X_{ib})$

위의 두 식을 나누면 아래와 같은 식을 구할 수 있는데, 보면 $\lambda(t)$은 없어지고, $X, \beta$로만 이뤄져있다. 즉, 시점 T에는 무관하게 위험이 동일하다고 할 수 있다. 

$\frac{\lambda(t|X_{a})}{\lambda(t|X_{b})} 
= \frac{exp(\beta_{1}X_{i1a}+...+\beta_{p}X_{ia})}{exp(\beta_{1}X_{i1a}+...+\beta_{p}X_{ib})} 
= exp(\beta(X_{a}-X_{b})) $

$HR=exp(\beta(X_{a}-X_{b}) = exp(\beta_{1}(X_{a1}-X_{b1})+...+beta_{p}(X_{ap}-X_{bp}))$을 의미하고, 각 $\beta_{p}$은 $X_{p}$에 따른 위험을 의미한다. 만약 다른 환자 a와 환자b가 모든 변수가 동일하고 $x_{1}$만 다르면 어떻게될까?

$\frac{\lambda(t|X_{a})}{\lambda(t|X_{b})} 

= \frac{h(t|x_{1}=1)}{h(t|x_{1}=0)}
= exp(\beta_{1}(1-0))
= exp(\beta_{1})
$

위처럼 모든 $x_{1}$ 변수만 제외하고 모든 변수가 동일하였을때 계수 $\beta_{1}$이 $x_{1}$가 1단위 변화할때마다 $exp(\beta)$한다고 해석하고, 이를 harzard ratio 라고 한다.

Python3에서는 lifelines패키지가 생존분석을 지원하는데, 시각화가 R보다는 약하다. 아래의 예시는 Cox을 이용한 시각화 및 통계플롯인데, CI값이 안나와서 다소 아쉽지만. python3에서 지원하는게 어딘가 싶다. (https://buildmedia.readthedocs.org/media/pdf/lifelines/latest/lifelines.pdf)

from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter

rossi = load_rossi()
cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')
cph.print_summary()


cph.plot()

cph.plot_partial_effects_on_outcome(covariates='race', values=list(set(rossi.race)), cmap='coolwarm')