Survival Analysis-I
Sung Woo Ko
2022년 1월 14일
R을 이용한 생존곡선 그래프 그리기
1.’survival’, ‘survminer’ 패키지를 불러온다
library(survival)
library(survminer)## 필요한 패키지를 로딩중입니다: ggplot2## 필요한 패키지를 로딩중입니다: ggpubr2.“survminer” 패키지 내부에 있는 “lung” 데이터를 활용한다
“lung”데이터에서 성별에 따른 폐암 생존율 차이를 알아본다
head(lung) # time: 추적관찰 기간(days), status: 1(중단), 2(사망), sex : 1(남), 2(여)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 03.Log-rank test: ’survdiff’함수로 시행
Log-rank test란 생존확률이 관심있는 변수에 따라 차이가 있는지 검정하는 방법이다.
survdiff(Surv(time,status)~sex,data=lung) # 성별에 따른 생존곡선 차이 검정## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001
귀무가설: 두 군의 생존곡선은 차이가 나지 않는다
p-value <0.05로 귀무가설 기각된다
결론: 성별에 따라 생존곡선이 차이 난다
4. 생존 곡선 그래프 그리기: ’survfit’함수로 객체 생성
fit<-survfit(Surv(time,status)~sex,data=lung) # survfit함수로 생존분석에 필요한 객체를 'fit'이름으로 저장
fit ## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## n events median 0.95LCL 0.95UCL
## sex=1 138 112 270 212 310
## sex=2 90 53 426 348 550<해석>
sex=1의 median survival time은 270일
Sex=2의 median survival time은 426일
5. 생존곡선 그리기 1. - R default plot
코드 출처: ‘안과영역에서 생존분석 방법론 적용’ Journal of Retina 2018;3(1):1-12
plot(fit ,main = 'KM Curve by Sex',xlab = 'Days', ylab ='Survival rate', lty = 1:2 ,
col=1:2,bty = 'n', cex.lab = 1.4, cex.axis = 1.4, xlim = c(0, 1000),xaxt = 'n' ,mark.time = T)
# bty (박스타입, o,l,7,c,u,n 6가지 종류가 있음 )
# cex.lab ('Days', 'survival rate' 글자의 크기) <br>
# cex.axis (200,400,600 글자의 크기 ) <br>
# xaxt = 'n' (x축제거 ) <br>
# mark.time= 'T'(Censored data를 표시) <br>
axis(side = 1, at = seq(0, 1000, 100), labels = seq(0, 1000, 100), cex.axis = 1)
# axis (side=1 (아래쪽), 2(왼쪽) 3(위쪽) 4(오른쪽))
# at (x축의 눈금이 그려질 위치 지정)
text(10, 0.10, labels = 'p-value of Log-rank test = 0.001', cex = 0.8, adj = 0)
text(200, 0.85, labels = 'Female', cex = 0.8, adj = 0)
text(150, 0.55, labels = 'Male', cex = 0.8, adj = 0)
# labels (축 눈금에 라벨 데이터 입력 )<br>
# adj ( 0은 중앙, 양수면 글자가 왼쪽으로 밀림/ 음수면 오른쪽으로 밀림)6. 생존곡선 그리기2 : ggsurvplot 함수
’survminer’패키지의 ‘ggsurvplot’함수를 사용하면 멋진 그래프를 편리하게 그릴 수 있다
library(survminer)
library (survival)
ggsurvplot(fit,title="KM Curve by Sex", # 그래프의 제목
risk.table = TRUE, # 그래프 밑에 위험표 표시 여부
xlim= c(0,1000), # x축의 범위
break.time.by=200, # x축 시간을 200단위로 끊어서 표시
xlab="Time(Days)") # x축의 제목
그룹별 색상을 지정하지 않으면 r의 그룹별 default color대로 자동 배정됨
1) 생존곡선의 신뢰 구간 표시: ‘conf.int’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T) #'conf.int=T'면 95% 신뢰구간 표시, 'F'면 미표시
’conf.int’인자로 95% 신뢰구간 표시여부를 결정할 수 있다.
2) ‘Strata’ 제목 및 변수 이름 변경: ‘legend title’,‘legend.labs’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex", # strata 제목 변경
legend.labs =c("Male","Female")) # strata의 value 입력 
환자를 나눈 ‘strata’가 어떤 속성인지 표시 하기 위해’ legend.title’을 이용하여 ’sex’로 변경
Sex =1이 어떤 성별인지 알수 없으므로, 알아볼 수 있게 ‘legend.labs’ 으로 변경
‘Female’, ’male’은 누가봐도 ’sex’이므로, ’sex’를 지우고 싶다면 legend.title=""로 하면 된다.
3) 그룹별 색상 변경: ‘palette’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2")) #남자 색상을 'steelblue2', 여자를 'pink2'
그룹별 색상을 ‘palette’ argument로 지정할 수 있음
4) Risk table 높이 변경: ‘risk.table.height’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3)
‘risk.table.height’ argument를 이용하여 risk table 높이 조절 가능 (0-1사이, Defualt 는 0.25)
’risk.table.height’를 높이면, survival plot 그래프 높이는 그만큼 줄어든다
5) 그래프 바탕 디자인 변경: ‘ggtheme’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3,
ggtheme=theme_gray())
그래프 바탕디자인을 ‘ggtheme’ arguement로 이용하여 격자무늬 회색바탕으로 변경 가능
그 외에 ‘theme_bw()’, ’theme_light()’등으로 취향에 맞게 변경가능 하다
6) x,y축 제목의 크기, 색상 변경: ‘font.x’,‘font.’y’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3,
font.x=c(15,"bold.italic",'red'),
font.y=15)
x축 제목 글자크기를 15, 굵게, 빨간색으로 변경
y축 제목 글자크기를 15로 변경
7) x,y축 눈금 숫자의 크기, 색상 변경: ‘font.tickslab’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3,
font.x=c(15,"bold.italic",'red'),
font.y=15,
font.tickslab=c(15,'darkgreen'))
x 축 눈금의 숫자의 크기가 커지고, 색깔이 녹색으로 변경된다
8) 그래프에 medain survival time 표시: ‘surv.median.line=hv’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3, surv.median.line='hv')
50% 생존률에 해당하는 median survival time 이 눈금선으로 표시된다
9) Risk table x축 삭제: ‘tables.theme = theme_cleantable()’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3, surv.median.line='hv',
tables.theme = theme_cleantable())
’Time(days)’는 survival plot에서 risk table에도 이중으로 표시되어 있다
Risk.table X 축에 ‘Time’ (days)를 삭제하고 싶으면 tables.theme = theme_cleantable()로 설정한다
10) ‘log rank test의 p-value표시: ’pval=T’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.3,
surv.median.line='hv',
tables.theme = theme_cleantable(),
pval = T)
Log rank test의 p-value를 ’pval=T’로 추가할 수 있다
11) Y축 값을 0~100%로 변경: fun=‘pct’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.25,
surv.median.line='hv',
tables.theme = theme_cleantable(),
pval = T,
fun='pct')
Y축 생존율값이 기존 0-1에서 0-100%로 변경 된다
12) 그래프 안에 원하는 텍스트 표기: ‘annotate(text)’
survcurv<-ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T, #그래프를 'survcurv'로 저장
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.25,
surv.median.line='hv',
tables.theme = theme_cleantable(),
pval = F, #pval=F로 변경
fun='pct')
survcurv$plot<-survcurv$plot+annotate("text", # 그래프에 텍스트 추가
x=700, # 텍스트 x 축 위치
y=60, # 텍스트 y 축 위치
label="p-value for log rank test: 0.001")
survcurv
text표기는 ’ggsurvplot’함수 자체의 명령어로 작성할 수 없다.
그래프를 저장한다음, 거기에다가 ’annotate’함수를 덧붙여서 텍스트를 추가한다.
13) 생존곡선을 누적 사망률로 변경: ‘fun=’event’
ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)",
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.25,
surv.median.line='hv',
tables.theme = theme_cleantable(),
pval = F, #pval=F로 변경
fun='event')
fun=’event’로 하면 누적 사망률 곡선으로 변경 가능하다
14) 누적 사망률 곡선 y축을 %로 변경: “scale_y_continuous”
cumcurv<-ggsurvplot(fit,title="KM Curve by Sex",
risk.table = TRUE,xlim= c(0,1000),
break.time.by=200,
xlab="Time(Days)", #누적 사망 그래프를 'cumcurv'로 저장
conf.int = T,
legend.title="Sex",
legend.labs =c("Male","Female"),
palette = c("steelblue2", "pink2"),
risk.table.height=0.25,
surv.median.line='hv',
tables.theme = theme_cleantable(),
pval = F, #pval=F로 변경
fun='event')
cumcurv$plot = cumcurv$plot +
scale_y_continuous(
breaks = seq(0,1,by = 0.2), # 0-1까지를 0.2단위로 쪼갬
labels = paste0(seq(0,100,by = 20),"%"))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.cumcurv
fun=’event’로 하면, 누적 사망률로 변경되는 대신 y 값이 0-1로 자동으로 바뀐다.
따라서 누적 사망률 그래프에서 y 값0-100%로 변경하고 싶다면, 기존 그래프를 따로 저장하여 ‘scale_y_continujous’ 함수를 추가하여 변경한다.
7. R 그래프를 컴퓨터에 저장하기: ‘ggsave’
ggsave(print(survcurv),file="./Fig.1A.jpg",width=10,height=6,units="in",dpi = 300)
‘survcurv’.를 ’Fig.1A’로 저장한다 (파일 형식은 JPG).
’Fig.1A’의 너비/높이는 10/6인치 이며, 해상도는 300 dpi이다
8. 특정시간의 생존율 구하기: time= ()
summary(fit,time=c(1,3)*365.25) ## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## sex=1
## time n.risk n.event survival std.err lower 95% CI
## 365.2500 35.0000 85.0000 0.3361 0.0434 0.2609
## upper 95% CI
## 0.4329
##
## sex=2
## time n.risk n.event survival std.err lower 95% CI
## 365.2500 30.0000 36.0000 0.5265 0.0597 0.4215
## upper 95% CI
## 0.6576
summary 함수에 ‘time’ argument로 특정 시간의 생존율을 확인한다.
시간변수가 month, year이면 ’time’을 정수로 입력.
시간 변수가 day 이면 ‘365.25’ * n 으로 입력.
결과 해석
sex1 (남자)의 1년, 3년 생존율은 각각 33.6%, 3.5%.
sex2 (여자)의 1년, 3년 생존율은 각각 52.6%, 8.3%.