Intro
Survival analysis, 한국어로 생존분석은 관심있는 사건이 발생하는 시간을 통계적으로 분석하는 방법이다. 생존이라 이름이 붙은 이유는, 일반적으로 의료 분야에서 환자의 사망까지 걸리는 시간에 대한 분석을 하기때문으로 생각된다. (뇌피셜임!) 때문에 꼭 사망을 관심있는 사건으로 보지 않아도 된다. 가령, 환자의 예후예측뿐 아니라, 고객의 상품에 대한 이탈 분석에도 사용이 가능하다. 이 포스터에서는 관심있는 사건을 사망으로 보며, 간단히 생존분석의 개념과 비모수 방법중 Kaplan-Meier 추정에 대해 알아볼 예정이다. 이 과목이 디테일하게는 대학원 이상의 수준을 다루고있어서 정보 손실의 최소화라거나, 다른 추정 방식 및 수학적 엄밀함을 포함하는 전반적인 내용은 따로 정리할 예정이다.
두 그룹 간의 생존시간 비교를 위하여 어떻게 분석을 시작해야할까? 우선 가설검정의 관점에서 봤을때 다음과 같은 귀무가설을 갖는다, "$H_0:$ There is no difference in survival between the two groups". 이에따라 생존분석은 생존시간(Time), 중도절단을 포함하는 사건발생 여부(event), 그리고 그룹정보(group) 총 세 가지 정보를 필수로 사용한다.
Concept
Time
Time($t$)은 관찰 시작으로 부터의 경과시간으로, 절대적인 시점에서의 값이 아니다. 예를들어, 암 진단받고 6개월간 살아남았다면, $t=6$ 이 된다. 약간 상대적인 느낌이 든다.
Event
Event는 사건의 발생, 대상의 입장에서는 일어나지 않았으면 하는 사건, 죽음이다. 생존분석의 뉘앙스는 그냥 별 말없으면 계속 살아있는 것이고, 사건이 발생하면 죽은 것이다. 중도절단(Sensored)이란, 사건이 일어나지 않은 상태에서 더 이상 관측이 불가능한 상황을 의미한다. 주로 관심 있는 사건이 연구 종료 시점까지 발생하지 않은 상태에서, 참여자에 대한 추적을 중단한 경우에 발생한다. 비록 완전한 정보는 아니지만, 관측된 시점까지는 살아남았다는 자료를 반영하여 정보의 손실을 최소화 한다. Event는 이러한 sensored 자료를 포함한다. 때문에 event 발생은 1로, 중도절단이나 event 발생하지 않는다면 0으로 묶어서 처리한다.
참고로 이 중도절단된 자료 자체는 무작위적이며, 독립적이며, 정보성을 보이면 안된다. 자료가 이 세가지 가정만 만족된다면, 생존분석에 사용 가능하다.
Data
아래의 왼쪽 그래프에서 각각의 화살표는 sample이다. 오른쪽은 time과 event 정보를 포함하는 tabular 자료로 왼쪽의 sample 순서와 일치한다. Event가 일어나는 경우만 1이고, 아닌경우는 그냥 0이다.
Survival Function
생존 함수(Survival function), $S(t)$는 주어진 시간 $t$ 동안 사건이 발생하지 않을 확률, 즉 개체의 생존할 확률을 나타낸다. 시점 $t$ 직전까지 사건이 한번이라도 일어날 확률을 $F(t)=P(0\leq T < t)$라고 할때, 아래와 같이 나타낼 수 있다:
$$ S(t)=P(T \geq t)=1-F(t) $$
생존 함수는, event가 적어도 $t$ 이전에는 일어나지 않을 확률이기 때문이다. 그렇다면 $S(0)=1$, $\lim_{t\to\infty}S(t)=0$ 라는 기본 사실을 생각해볼 수 있을 것이다. 간단하게, $t$ 시간이 지난 시점에서의 살아있는 사람의 비율을 생각하면 된다.
Hazard Function
위험 함수는 $h(t)$로 나타내며 $t$시점에서의 "Intensity of death"를 의미한다. 해당 시간까지의 생존을 고려할 때, 사건이 발생할 순간적인 위험을 설명한다. 조금 더 풀어 이야기하자면 위험 함수는, $t$ 시점까지 살아남았을 때, 그 때의 사망확률을 이야기한다. 결코 가볍게 다룰 내용은 아니지만, 설명을 위해 아래와 같은 상황을 예로들수있다.
"몇시지?!"
"진단 받고나서 $t$만큼 시간이 흘렀군... 난 살아남았어..!!"
($\epsilon$만큼의 시간이 흐른 뒤, 거의 바로 그 순간...!)
"킄...!!!(사망)"
이런 일이 일어날만한 확률인 셈이다.
수학적으로 식을 세우자면, 조건부 확률과 미분의 개념이 들어가면 되어 아래와같이 나타낼 수 있다.
$$h(t)=\lim_{\Delta t\to 0}{\frac{P(t \leq T \leq t + \Delta t | T \geq t)} {\Delta t}}$$
Relationships
이 두 함수는 관계성이 있다.
$$h(t)=\lim_{\Delta t\to 0}{\frac{P(t \leq T \leq t + \Delta t | T \geq t)} {\Delta t}}$$
조건부 확률,
$$=\lim_{\Delta t\to 0}{\frac{P(t \leq T \leq t + \Delta t, T \geq t)} {\Delta t \times P(T \geq t)}}$$
생존함수의 정의,
$$=\lim_{\Delta t\to 0}{\frac{P(t \leq T \leq t + \Delta t, T \geq t)} {\Delta t \times S(t)}}$$
$S(t)$제외한 부분은 결국 $t$ 시점에 사건이 일어날 확률이다보니,
$$=\frac{p(t)} {S(t)}$$
where $p(t)=\frac{d}{dt}F(t)$.
그리고 합성함수 미분법을 생각하면,
$$=-\frac{\partial \log S(t)}{\partial t}$$
Cumulative Hazard Function
위의 두 함수 말고도, 위험도를 재는 다른 함수가 또 있다. 바로 누적 위험함수($H(t)$)이다. 말그대로 아래와 같이 정의된다.
$$H(t)=\int_0^t h(s) ds$$
$h(t)$가 $t$시점에서의 "Intensity of death"를 의미했다면, $H(t)$는 $t$시점 까지의 "Amount of hazard"를 나타낸다. 이 역시 다른 두 함수와 다음과 같은 연관성을 갖는다.
$$H(t)=-\log S(t)$$
$$S(t)=e^{-H(t)}=e^{-\int_0^t h(s)ds}$$
Kaplan-Meier Estimator
Kaplan-Meier (KM) estimator는 생존함수를 추정하는 가장 간단한 방법중 하나이다. 아래와 같이 추정한다:
$$\hat{S(t)_{KM}}=\prod_{i=1}^t (1-d_i/n_i)$$
where $d_i:$ # of events at time $i$ and $n_i:$ # of records at time $i$
$d_i/n_i$ 가 결국 시점 $i$에서 죽은 확률이므로, $ \hat{S(t)_{KM}} $가 관찰 시간에 따른 사건발생 시점마다의 사건 발생률을 계산한다는 의미가 된다. 아래 그림의 x축은 time, y축은 추정된 survival function value이다. 그래프 선의 색은 두 그룹을 분리하여 나타낸 것이다. 일종의 누적확률분포 그림같다.
이 외에도 Nelson-Aalen estimator와 같은 비모수 방법이 더 있지만, 이 포스터에서는 소개하지 않을 예정이다.
Log-rank test
Intro에서 이야기한 가설의 검증 차례이다. 두 그룹의 KM function에 대한 차이를 검정하기 위해, log-rank test를 진행한다. 검정 대상이되는 두 그룹은 각각 Cumulative distribution의 한 형태를 띄고있어, Kolmogorov-Smirnov test를 하여도 되지만, 대중적으로 log-rank test가 가장 많이 사용된다. 조금 더 자세한 내용은 첨부된 Detail.html/Log Rank Test 를 참고 바란다. 생각보다 재미있다!! ㅎㅎ
Implementations
기존에 구현되어있는 survival analysis 라이브러리를 사용하지 않고, scratch부터 구현해보았다. 아래의 첨부된 문서에 구현해 놓았다. 데이터는 설명 파일에 어떤 형태인지 확인 가능할 것이다.
Data loading
주어진 자료의 형태를 고려하여 읽어들이는 함수가 들어있는 스크립트.
Survival Analysis
중도절단 열을 추가하는 자잘한 함수부터 해서, KM 추정과 log-rank test하는 일련의 모든 계산을 함수로 구현해놓은 스크립트.
Driver
구현된 모듈을 돌리고 검증하는 driver 스크립트와 마크다운.
Details
이론과 함께 전반적인 설명이 들어있는 마크다운.
파일을 정리하자면,
Detail.nb.html: 구체적인 작업 흐름을 보여주기 위한 파일
Driver.nb.html: 기존의 라이브러리와의 비교를 위한 파일
mySurv.R: 생존분석을 위해 구현한 함수가 있는 소스코드
dataLoading.R: 주어진 자료파일을 전처리하기 위한 함수가 있는 소스코드
Outro
생존분석은 학부때부터 사용은 해왔지만, 이론을 이해하기가 쉽지 않았다. 그래도 시간이 지나고 다시 보니 은근 재미있고, 조금은 더 수식이 이해가 되는 부분이 많았다. 개념과는 별개로, KM을 추정하여 계산하는 일련의 과정을 코드로 구현하였다. 어쩌다 구현할 일이 있어 하게되었는데, 당시에는 이틀만에 급하게 구현하여 허접해보이는 코드가 나왔다. 결과가 기존의 라이브러리와 같아 나름 자부심있었지만, 다른 자료에 적용했을때에는 다른 결과가 나올때가 있다고 전달받아, 여러모로 아쉬움이 많이 남는 스크립트다. 나는 아직 그 경우를 못 찾았다..! ㅠㅜ 그래도 역시 scratch 부터의 구현은 대부분 재미있는것 같다! ㅎㅎ
References
Wiki
https://faculty.washington.edu/yenchic/18W_425/Lec5_survival.pdf
https://webfocusinfocenter.informationbuilders.com/wfappent/TLs/TL_rstat/source/Survival45.htm
Explanation of Survival Analysis
In this example, we are going to create a model for matching children with foster parents. We will use historical data that contains observations on the length of stay of children with foster parents. The predictor variables for the length of stay include
webfocusinfocenter.informationbuilders.com
'정리 조금 > Bioinformatics & Biostatistics' 카테고리의 다른 글
Weighted Gene Co-expression Network Analysis (0) | 2023.10.26 |
---|---|
Gene Set Enrichment Analysis (0) | 2023.08.29 |
The Characteristic Direction (0) | 2023.08.28 |
ZIBseq (0) | 2023.01.12 |
edgeR & DESeq2 (0) | 2023.01.11 |