R에서 생존분석: Kaplan-Meier 방법

R에서 생존분석: Kaplan-Meier 방법

준비

생존 분석 패키지를 다음 명령으로 읽어들인다.

> require(survival)

급성 골수성 백혈병(Acute Myelogenous Leukemia) 데이터가 예제로 들어있다. 작은 데이터이다. 내용을 보면:

> aml
   time status             x
1     9      1    Maintained
2    13      1    Maintained
3    13      0    Maintained
4    18      1    Maintained
5    23      1    Maintained
6    28      0    Maintained
7    31      1    Maintained
8    34      1    Maintained
9    45      0    Maintained
10   48      1    Maintained
11  161      0    Maintained
12    5      1 Nonmaintained
13    5      1 Nonmaintained
14    8      1 Nonmaintained
15    8      1 Nonmaintained
16   12      1 Nonmaintained
17   16      0 Nonmaintained
18   23      1 Nonmaintained
19   27      1 Nonmaintained
20   30      1 Nonmaintained
21   33      1 Nonmaintained
22   43      1 Nonmaintained
23   45      1 Nonmaintained

이 데이터는 단순한 데이터프레임이다. 확인해 보자.

> class(aml)
[1] "data.frame"
> str(aml)
'data.frame': 23 obs. of  3 variables:
 $ time  : num  9 13 13 18 23 28 31 34 45 48 ...
 $ status: num  1 1 0 1 1 0 1 1 0 1 ...
 $ x     : Factor w/ 2 levels "Maintained","Nonmaintained": 1 1 1 1 1 1 1 1 1 1 ...

즉, 단순한 테이블 데이터라는 뜻의다. 아직 특별히 생존 분석을 위한 가공이 이루어지 않은 것이다. 위와 같은 형태로만 데이터가 입력되어 있으면 된다.

Survival Object

생존 분석을 하기 위해 데이터를 Survival Object로 만들어 보자. 이 객체는 survival 패키지에서 제공하는 것이다. Survival Object는 Surv()를 이용하여 만든다. 이 함수는 기본적으로 Surv(time, event)로 사용할 수 있다. 두 변수가 필요하다. time은 생존시간이고 event는 상태(status)를 나타내는 변수이다. 상태는 기본적으로

• 0=right censored
• 1=event at time
• 2=left censored
• 3=interval censored

로 코딩한다. 지금 주어진 aml데이터는 status 변수에 0/1밖에 없다. 다음과 같이 Survival Object를 만들면 status 값이 0인 데이터 즉 right censored 데이터인 경우에 시간 뒤에 +가 붙어 표시되는 것을 확인할 수 있다.

> (aml.Surv <- with(aml, Surv(time, status)))
 [1]   9   13   13+  18   23   28+  31   34   45+  48  161+   5    5    8    8   12   16+  23   27   30   33   43 
[23]  45 

Kaplan-Meier method

Kaplan-Meier 방법으로 중도 절단이 있는 자료(censored data)의 생존 함수를 추정하는 것은 survfit()로 할 수 있다. 앞서 만든 Survival Object를 이용해서

> aml.survfit <- survfit(aml.Surv ~ aml$x)

라고 할 수도 있지만,

> aml.survfit <- survfit(Surv(time, status) ~ x, aml)

라고 하는 것이 보기 좋다. 분석 자료 개요는 다음과 같다.

> aml.survfit
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                records n.max n.start events median 0.95LCL 0.95UCL
x=Maintained         11    11      11      7     31      18      NA
x=Nonmaintained      12    12      12     11     23       8      NA

분석 결과를 보려면 summary() 함수를 이용하면 된다. 결과는 다음과 같다.

> summary(aml.survfit)
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                x=Maintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    9     11       1    0.909  0.0867       0.7541        1.000
   13     10       1    0.818  0.1163       0.6192        1.000
   18      8       1    0.716  0.1397       0.4884        1.000
   23      7       1    0.614  0.1526       0.3769        0.999
   31      5       1    0.491  0.1642       0.2549        0.946
   34      4       1    0.368  0.1627       0.1549        0.875
   48      2       1    0.184  0.1535       0.0359        0.944

                x=Nonmaintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     12       2   0.8333  0.1076       0.6470        1.000
    8     10       2   0.6667  0.1361       0.4468        0.995
   12      8       1   0.5833  0.1423       0.3616        0.941
   23      6       1   0.4861  0.1481       0.2675        0.883
   27      5       1   0.3889  0.1470       0.1854        0.816
   30      4       1   0.2917  0.1387       0.1148        0.741
   33      3       1   0.1944  0.1219       0.0569        0.664
   43      2       1   0.0972  0.0919       0.0153        0.620
   45      1       1   0.0000     NaN           NA           NA

생존 곡선 그리기

생존 곡선은 plot() 명령으로 간단히 그릴 수 있다.

> plot(aml.survfit)

위 데이터에서 이미 보았듯이 aml$x를 보면

> levels(aml$x)
[1] "Maintained"    "Nonmaintained"

두 개의 집단으로 이루어져 있다. 두 생존곡선을 다른 형태의 선으로 그려서 구별하고 싶다면

> plot(aml.survfit, lty=1:2)

라고 lty (line type) 옵션을 주는 것만으로도 충분하다. 범례도 달고 싶다면 일반적으로는 다음과 같이 하면 된다.

K <- length(aml.survfit$strata)
LEGEND <- attr(aml.survfit$strata, "names")
plot(aml.survfit, lty=1:K)
legend("topright", LEGEND, lty=1:K)

물론, 여기서 K나 LEGEND를 직접 타이핑해도 된다. 직접 써줄 때는 범례 이름을 바꾸어 쓰는 일이 없도록 주의해야 한다.

plot(aml.survfit, lty=1:2)
legend("topright", c("Maintained", "Nonmaintained"), lty=1:2)

생존 분석 기본 개념

생존 분석 기본 개념
생존 분석(survival analysis)이라고 하면 주로 생명체의 죽음 또는 기계 장치의 고장에 관련된 통계 분석이지만 이 분야에 한정된 것은 아니다. 공학 분야에서 reliability analysis, failure time analysis, 경제학 분야에서 duration analysis, transition analysis, 사회학 분야에서 event history analysis 등의 이름으로 불리며 이들 사이에 실질적인 분석 방법의 차이는 없다. 질병의 발생이나 기계의 오작동 뿐만 아니라 지진의 발생, 자동차 사고 발생, 주가의 추락, 등을 분석하는 데에 사용된다[1].
생존 분석을 더 일반적인 개념을 사용하여 'time to event' 자료의 분석이라고도 하는데 사망, 고장, 사고, 폭락 등을 사건(event)이라고 보고 해당 사건이 일어날 때까지 걸린 시간을 데이터로 다룬다는 의미이다. 이제 이런 사건을 생존 분석에서는 사망(death)라고 지칭하는 경우가 많고 좀더 일반적으로는 고장/실패(failure)로 지칭하는 경우가 많다. 하지만 반드시 '사망'이거나 '실패'일 필요는 없다. 무엇이든 우리가 관심을 가지는 사건이면 된다.

생존 시간

우리는 사건이 일어나는 시간에 관심이 있다. 이 사건이 발생하는 시간을 \(T\)라고 하자. 이것은 확률 변수이다. 생존 분석에서는 이것을 생존 시간(survival time)이라고 한다. 시간이 항상 0에서부터 시작한다고 보면 사건이 일어나는 시점 \(T\)는 사건이 일어나기까지 걸린 시간이기도 하다.

생존 함수

생존 함수(survival function 또는 suvivor function) \(S(t)\)는 생존 시간 \(T\)가 특정 시간 \(t\) 보다 클 확률을 가리킨다. \(T\)의 cdf를 \(F(t)\)라고 하면 다음 관계가 성립한다. \[ S(t) = P(T > t) = 1 - F(t) \] 생존 함수 \(S(t)\)는 \(t\) 시점에서 살아있을 또는 \(t\) 시점까지 살아남을 확률이다. 이때 \(t\) 시점까지 살아남는다는 것은 \(t\) 시점까지만 살고 바로 그때 죽는다는 뜻이 아니다. \(t\) 시점 이후에 죽는다는 뜻이다. 더 일반적으로 말로 표현하면 \(t\) 시점까지 관심 사건이 일어나지 않을 확률이다.

위험 함수

위험 함수(hazard function) \(h(t)\)는 주어진 시간 \(t\)에서 단위 시간 당 사건 발생율이다. \[ h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t+ \Delta t |T \geq t)}{\Delta t} \] 사건이 일어날 순간 가능성(instantaneous potential)이라고도 할 수 있다. 위험 함수를 조건부 실패율(conditional failure rate)라고도 한다.

함수들 사이의 관계

생존 함수 \(S(t) = P(T>t)\)와 생존 시간 \(T\)의 cdf인 \(F(t)\)사이에는 \[ S(t) = 1 - F(t) \] 관계가 성립하므로 양변을 \(t\)에 관해 미분하면 \[ \frac{d}{dt}S(t) = - f(t) \] 와 같은 관계를 얻을 수 있다. 여기에서 \(f(t)\)는 \(T\)의 pdf이다.

위험 함수 \(h(t)\)의 정의를 따라 식을 바꾸어 보면 \begin{align*} h(t) &= \lim_{\Delta t \to 0} \frac{P(t \leq T < t+ \Delta t |T \geq t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{P(t \leq T < t+\Delta t)}{P(T \geq t)}\frac{1}{\Delta t}\\ & = \frac{1}{P(T \geq t)}\lim_{\Delta t \to 0} \frac{P(t\leq T < t+\Delta t)}{\Delta t}\\ & = \frac{1}{P(T \geq t)}f(t)\\ & = \frac{f(t)}{S(t)}\\ & = - \frac{d}{dt}\ln S(t) \end{align*}

누적 위험 함수

누적 위험 함수 \(H(t)\)는 다음과 같다. \[ H(t) = \int_0^t h(u)du = - \ln S(t) \] 따라서 \[ S(t) = e^{H(t)} \] 이다.

참고

1. Survival Analysis Using SAS: A Practical Guide

생존 분석 좋은 책

생존 분석 좋은 책

• Klein and Moeschberger (2010), Survival Analysis: Techniques for Censored and Truncated Data, Springer. 대학원 석사 수준 생존 분석. 기본적인 통계 지식만 있어도 이해할수 있음.
• Kleinbaum and Klein (2005), Survival analysis: a self-learning text, Springer.
• Hosmer, Stanley, and May (2008), Applied survival analysis: regression modelling of time-to-event data, Wiley-Interscience.
• Themeau and Grambsch (2000), Modelling survival data: extending the Cox model, Springer.