생존분석: R을 이용한 상대생존률 구하기
Keon-Woong Moon
October 20, 2016
생존분석 : 관찰생존률
의학연구에 있어 생존분석은 상당히 중요하다. 일반적으로 의학자료의 생존기간은 진단일로부터 사망일까지를 생존기간으로 한다. 생존률은 전체환자 중 사망한 환자의 비율로 예를 들어 2000년에 진단된 10명의 폐암 환자를 1년간 추적관찰하였더니 4명이 사망하였다면 1년 사망률은 40%(= \({4\div10}\times100\))이 되고 생존률은 60%(= \({6\div10}\times100\))이 된다. 이 때의 생존률은 관찰생존률이다.
원인별 생존확률, 순생존률, 교정생존률
하지만 관찰기간 중 사망한 4명 가운데 3명은 질병으로 인해 사망하였지만 1명은 교통사고로 인해 사망하였다면 암 환자의 1년 생존률이 60%라고 할 수는 없다. 교정생존률(correctred survival rate)은 순생존률(net survival rate) 또는 원인별생존확률(cause-specific survival rate)라고도 하는데 사망원인이 교통사고와 같이 해당 질병과 관련 없는 경우 이를 불완전관찰치로 처리하여 생존률을 계산하는 것이다. 예를 들어 악성 흑색종 환자 50명을 2년간 추적 관찰했을 때 첫해에 질병으로 인한 사망이 8명, 다른 원인으로 인한 사망이 1명이고 누락 환자는 없었고 1년-2년 사이에 질병으로 인한 사망이 4명 다른원인으로 인한 사망이 2명, 두명의 누락(불완전 관찰)이 있었다고 할 때 사망률은 다음과 같이 계산할 수 있다.
| 년도(i) | 년도시작시 생존자(li) | 1년간 질환사망(ddi) | 1년간 다른 원인사망(doi) | 누락(wi) | effective n at risk(ri) | 년간사망률(qi) | 년간 생존률(pi) | 누적 사망률 (\(\Pi\)i) |
|---|---|---|---|---|---|---|---|---|
| 0 | 50 | 8 | 1 | 0 | 49.5 | 0.162 | 0.838 | 0.838 |
| 1 | 41 | 4 | 2 | 1 | 39.5 | 0.101 | 0.899 | 0.754 |
이 때 \(ri=li-{wi+doi}\div2, qi=ddi\div ri\) 로 계산된다.
상대생존률(relative survival rate)
임상시험인 경우 대부분 연구대상자 모두에 대한 완전한 추적이 가능하며 위와 같이 순사망률을 계산할 수 있다. 하지만 지역주민을 대상으로 통계청의 사망통계를 근거로 사망 여부를 추적하는 역학조사라면 사망자 개개인별로 사망원인을 정확하게 알 수 없는 한계에 접하게 된다. 또한 사망원인이 질병으로 인한 것인지 명확하지 않을 경우도 있다. 이러한 경우 상대생존률을 사용할 수가 있는데 상대생존률은 연구대상자로부터 얻어진 관찰생존률을 기대생존률로 나누어 얻을 수 있다.
\(Relative\ Survival\ Rate = Observed\ Survival\ Rate \div Expected\ Survival\ Rate\)
여기서 기대생존률은 연구대상자가 속해있는 집단에서 해당 질병이 없는 가운데 생존할 확률을 의미한다. 우리나라의 경우 통계청에서 제공하는 완전생명표를 이용하여 기대생존률을 구할 수 있다.
기대생존률 구하기
완전생명표
국가통계포털(kosis.kr)에서 국내통계-> 주제별통계->인구가구->생명표->완전생명표에서 남여 성별 1세간격으로 조사된 완전생명표를 다운로드 받을 수 있다. 1997년, 1999년자료가 있고 2001년부터는 매년 자료가 있다. 이 생명표를 이용하면 우리나라의
R로 관찰생존률 구하기
먼저 생존분석에 필요한 survival 패키지와 relsurv 패키지를 설치한다. 이미 설치한 경우는 다시 설치할 필요가 없다.
install.packages(c("survival","relsurv"))패키지를 불러온다.
library(survival) # surexp() 함수의 사용을 위해
library(relsurv) # rdata, slopop의 사용을 위해 Loading required package: splinesLoading required package: daterelsurv패키지에 있는 rdata를 사용할 예정인데 이 데이터는 1040명의 데이터로 다음과 같은 구조를 갖는다.
data(rdata)
str(rdata)'data.frame': 1040 obs. of 6 variables:
$ time : int 2657 1097 3764 3724 5076 139 4940 5078 596 4635 ...
$ cens : int 1 1 1 1 0 1 1 0 1 1 ...
$ age : int 68 63 60 66 57 57 67 59 43 50 ...
$ sex : num 2 2 1 2 2 2 1 1 1 2 ...
$ year : int 8210 8278 8254 8054 8224 8233 8335 8177 8288 8281 ...
$ agegr: Factor w/ 4 levels "<54","54-61",..: 3 3 2 3 2 2 3 2 1 1 ...head(rdata) time cens age sex year agegr
1 2657 1 68 2 8210 62-70
2 1097 1 63 2 8278 62-70
3 3764 1 60 1 8254 54-61
4 3724 1 66 2 8054 62-70
5 5076 0 57 2 8224 54-61
6 139 1 57 2 8233 54-61여기서 눈여겨봐야 할 것은 time열이 생존시간을 나타내는데 일수(days)로 기록되어 있다는 점이다. 이 데이터를 이용하여 성별로 생존곡선을 만들면 다음과 같다.
fit=survfit(Surv(time,cens)~sex,data=rdata)
fitCall: survfit(formula = Surv(time, cens) ~ sex, data = rdata)
n events median 0.95LCL 0.95UCL
sex=1 751 360 4204 3810 4596
sex=2 289 187 2203 1720 2823이 데이타에서 1년, 3년, 5년 관찰생존률을 구하려면 다음과 같이 볼 수 있다.
summary(fit,time=c(1,3,5)*365.24)Call: survfit(formula = Surv(time, cens) ~ sex, data = rdata)
sex=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
365 683 66 0.912 0.0103 0.892 0.933
1096 609 67 0.822 0.0140 0.795 0.850
1826 512 61 0.738 0.0162 0.707 0.771
sex=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
365 244 45 0.844 0.0213 0.804 0.887
1096 189 50 0.669 0.0278 0.617 0.726
1826 148 34 0.548 0.0296 0.493 0.609여기서 남성의 1년 관찰생존률은 0.912(95%신뢰구간 0.892-0.933)이고 여성의 5년 관찰생존률은 0.548(95% 신뢰구간 0.493-0.609)임을 알수 있다.
library(survminer)Loading required package: ggplot2ggsurvplot(fit)
censor 마크(+)를 없애고 95%신뢰구간을 표시하려면 다음과 같이 한다.
ggsurvplot(fit,censor=FALSE,conf.int=TRUE)
risk.table을 표시하려면 다음과 같이 한다.
ggsurvplot(fit,censor=FALSE,conf.int=TRUE,risk.table=TRUE)
R로 기대생존률 구하기
R로 기대생존률을 구하는 함수는 survexp()함수로 survival 패키지에 포함되어 있다. 기대생존률을 계산하기 위해서는 성별,나이별,연도별로 생존률을 정리한 ratetable이 필요한데 미국의 생존률표는 survexp.us로 역시 survival패키지에 포함되어 있다.
위의 데이타에서 미국인구의 생존률을 기준으로 상대생존률을 구하려면 다음과 같이 할 수 있다.
fit1=survexp(Surv(time,cens)~sex,data=rdata,ratetable=survexp.us)
summary(fit1,time=c(1,3,5)*365.24)Call: survexp(formula = Surv(time, cens) ~ sex, data = rdata, ratetable = survexp.us)
time nrisk1 nrisk2 survival1 survival2
365 683 244 0.996 0.997
1096 609 189 0.994 0.996
1826 512 148 0.994 0.995survexp.us데이터는 미국인구의 생존데이터로 1960부터 1980년 까지의 사망률 데이터이다. relsurv패키지에 포함되어 있는 slopop데이타는 슬로베니아 인구의 1930-2014의 인구통계다.
data(slopop)
str(slopop) ratetable [1:104, 1:39, 1:2] 4.43e-04 7.15e-05 3.34e-05 2.16e-05 1.76e-05 ...
- attr(*, "dimnames")=List of 3
..$ : chr [1:104] "0" "1" "2" "3" ...
..$ : chr [1:39] "1930" "1948" "1952" "1960" ...
..$ : chr [1:2] "male" "female"
- attr(*, "dimid")= chr [1:3] "age" "year" "sex"
- attr(*, "factor")= num [1:3] 0 0 1
- attr(*, "type")= num [1:3] 2 3 1
- attr(*, "cutpoints")=List of 3
..$ : num [1:104] 0 365 730 1096 1461 ...
..$ :Class 'date' int [1:39] -10957 -4383 -2922 0 3653 7305 8036 8401 8766 9132 ...
..$ : NULLattr(slopop,"dimnames")[[1]]
[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
[12] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21"
[23] "22" "23" "24" "25" "26" "27" "28" "29" "30" "31" "32"
[34] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43"
[45] "44" "45" "46" "47" "48" "49" "50" "51" "52" "53" "54"
[56] "55" "56" "57" "58" "59" "60" "61" "62" "63" "64" "65"
[67] "66" "67" "68" "69" "70" "71" "72" "73" "74" "75" "76"
[78] "77" "78" "79" "80" "81" "82" "83" "84" "85" "86" "87"
[89] "88" "89" "90" "91" "92" "93" "94" "95" "96" "97" "98"
[100] "99" "100" "101" "102" "103"
[[2]]
[1] "1930" "1948" "1952" "1960" "1970" "1980" "1982" "1983" "1984" "1985"
[11] "1986" "1987" "1988" "1989" "1990" "1991" "1992" "1993" "1994" "1995"
[21] "1996" "1997" "1998" "1999" "2000" "2001" "2002" "2003" "2004" "2005"
[31] "2006" "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014"
[[3]]
[1] "male" "female"같은 rdata를 이용하더라도 ratetable이 달라지면 기대생존률 또한 달라진다.
fit2=survexp(Surv(time,cens)~sex,data=rdata,ratetable=slopop)
summary(fit2,time=c(1,3,5)*365.24)Call: survexp(formula = Surv(time, cens) ~ sex, data = rdata, ratetable = slopop)
time nrisk1 nrisk2 survival1 survival2
365 683 244 0.987 0.990
1096 609 189 0.986 0.990
1826 512 148 0.986 0.989우리나라의 인구통계
우리나라의 완전생명표는 국가통계포털(kosis.kr)에서 구할 수 있다. 국내통계-> 주제별통계->인구가구->생명표->완전생명표(각세별)에서 1997년부터 2014년까지의 완전생명표를 다운로드 받을 수 있는데 이를 이용하여 우리나라의 ratetable을 만들수 있다. relsurv패키지의 transrate()함수를 이용하면 우리나라의 ratetable을 만들 수 있는데 1997년자료는 85세 이상자료가 없으며 1999년에는 95세 이상의 자료가 없으므로 약간의 data handling이 필요하다 또한 1998년과 2000년자료는 누락이 되어 있어 1998년은 1997년 자료로, 2000년 자료는 1999년 자료로 대신 입력하여 우리나라의 ratetable을 만들 수 있었다. 완성된 ratetable은 github에 있는 moonBook2패키지에 ratetable.kor로 포함되어 있다. R에서는 다음과 같이 읽어서 사용하면 된다.
install.packages("devtools") ## 설치되어 있지 않은 경우 한번만 실행
devtools::install_github("cardiomoon/moonBook2") ## 설치되어 있지 않은 경우 한번만 실행library(moonBook2)
attr(ratetable.kor,"dimnames")[[1]]
[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
[12] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21"
[23] "22" "23" "24" "25" "26" "27" "28" "29" "30" "31" "32"
[34] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43"
[45] "44" "45" "46" "47" "48" "49" "50" "51" "52" "53" "54"
[56] "55" "56" "57" "58" "59" "60" "61" "62" "63" "64" "65"
[67] "66" "67" "68" "69" "70" "71" "72" "73" "74" "75" "76"
[78] "77" "78" "79" "80" "81" "82" "83" "84" "85" "86" "87"
[89] "88" "89" "90" "91" "92" "93" "94" "95" "96" "97" "98"
[100] "99" "100"
[[2]]
[1] "male" "female"
[[3]]
[1] "1997" "1998" "1999" "2000" "2001" "2002" "2003" "2004" "2005" "2006"
[11] "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014"rdata의 생존을 우리나라 인구통계를 기준을 본 1년,3년,5년의 기대 생존률은 다음과 같다.
fit3=survexp(Surv(time,cens)~sex,data=rdata,ratetable=ratetable.kor)
summary(fit3,time=c(1,3,5)*365.24)Call: survexp(formula = Surv(time, cens) ~ sex, data = rdata, ratetable = ratetable.kor)
time nrisk1 nrisk2 survival1 survival2
365 683 244 0.994 0.994
1096 609 189 0.992 0.993
1826 512 148 0.991 0.992상대생존률 구하기
이 데이타의 성별 관찰생존률은 다음과 같았다.
fit=survfit(Surv(time,cens)~sex,data=rdata)
summary(fit,time=c(1,3,5)*365.24)Call: survfit(formula = Surv(time, cens) ~ sex, data = rdata)
sex=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
365 683 66 0.912 0.0103 0.892 0.933
1096 609 67 0.822 0.0140 0.795 0.850
1826 512 61 0.738 0.0162 0.707 0.771
sex=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
365 244 45 0.844 0.0213 0.804 0.887
1096 189 50 0.669 0.0278 0.617 0.726
1826 148 34 0.548 0.0296 0.493 0.609따라서 우리나라 인구통계를 기준으로 본 1,3,5년 상대생존률은 다음과 같다.
남자1년 OSR : 0.912 남자1년: 기대생존률 0.994 상대생존률: 0.918(=0.912/0.994)