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)