생존분석 그래프(ggplot2, ggsurvminer) 그리기

Notice

Recent Posts

Recent Comments

Link

« 2024/04 »
일	월	화	수	목	금	토
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30

Tags more

Archives

Today

Total

관리 메뉴

Must Learning With Statistics

생존분석 그래프(ggplot2, ggsurvminer) 그리기 본문

R ggplot2

생존분석 그래프(ggplot2, ggsurvminer) 그리기

Doublek Park 2020. 1. 30. 01:52

1. 생존분석 그래프의 중요성

생존분석은 의학계열에서 가장 자주 쓰이는 분석방법 중 하나입니다. 논문에 많이 활용되는 만큼, 생존분석 그래프또한 논문에 실리게 됩니다. 관계자의 말에 따르면 생존분석 그래프를 잘 그려야 논문 Accept에 유리하다는 소식을 들었습니다. 이번 포스팅은 생존 그래프를 어떻게 하면 잘 그릴 수 있을지에 대해 다룹니다.

2. 생존분석 시각화(카플란 마이어)

이번 글은 생존분석의 시각화가 목적이기 때문에 생존분석의 이론은 다루지 않고 바로 진행합니다.

library(survival)
library(survminer)

survival 패키지 : 생존분석 모델링을 위한 패키지입니다.
survminer 패키지 : 생존분석 시각화를 위한 패키지입니다.

데이터 다운로드 링크: https://www.dropbox.com/sh/vtqlvrgdts2yfez/AAD_cd49dBcvgBNdz-C-A6TFa?dl=0

Survival = read.csv("D:\\Dropbox\\DATA SET(Dropbox)\\Survival_Sample.csv")

시각화에 사용할 데이터는 다음과 같습니다.

- TIME : 생존기간
- EVENT : 이벤트 발생여부(사망)
- GENDER : 성별

Kaplan Meier

카플란 마이어 생존분석은 관심 사건(환자의 사망, 기업의 파산, 유저의 복귀 등)이 발생한 시점마다 구간생존율을 계산하여 이의 누적생존율을 추정하는 비모수적 분석방법입니다.

\[ p[t]=\frac{n(Survived_{t})}{n(Observed_{t})} \\ S[t]=S(t-1)P[t],\ t는\ 이벤트\ 발생\ 시점. \]

KM = survfit(Surv(event = EVENT,time = TIME) ~ 1, data = Survival)
summary(KM)

Call: survfit(formula = Surv(event = EVENT, time = TIME) ~ 1, data = Survival)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     37       1   0.9730  0.0267      0.92210        1.000
   11     30       1   0.9405  0.0410      0.86352        1.000
   13     27       1   0.9057  0.0522      0.80892        1.000
   19     23       2   0.8269  0.0715      0.69812        0.980
   21     20       1   0.7856  0.0789      0.64516        0.957
   24     17       1   0.7394  0.0868      0.58746        0.931
   27     14       1   0.6866  0.0953      0.52304        0.901
   28     13       2   0.5809  0.1059      0.40636        0.831
   30     11       1   0.5281  0.1087      0.35285        0.791
   31     10       1   0.4753  0.1099      0.30212        0.748
   34      9       3   0.3169  0.1046      0.16590        0.605
   35      5       2   0.1901  0.0936      0.07245        0.499
   36      3       1   0.1268  0.0811      0.03619        0.444
   37      2       1   0.0634  0.0604      0.00978        0.411
 [ getOption("max.print") 에 도달했습니다 -- 1 행들을 생략합니다 ]

카플란 마이어 분석을 진행한 결과, 시점에 따른 생존율이 계산이 되는 것을 확인할 수가 있습니다. 분석결과에 대해 시각화를 진행해보도록 하겠습니다. ggsurvplot이라는 함수를 주로 사용합니다. 이 함수는 ggplot2를 기반으로 하지만, 조금은 다른 점이 있습니다.

ggsurvplot(KM,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           palette = c("#2E9FDF"), 
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           surv.median.line = "hv", # 50% 생존지점 표시
           break.time.by = 5,
           xlim = c(0,40)
           )

ggsurvplot()에 카플란 마이어 식을 넣어줌과 함께, 다양한 옵션을 통해 그래프를 그릴 수가 있습니다. 특히 50% 생존지점을 표현하는 surv.median.line 옵션의 경우에는 ’hv’를 주면 수직과 대각선으로 선이 그어지며, ’h’를 줄 경우에는 수평선으로, ’v’를 줄 경우에는 수직선으로 표현이 됩니다.

ggsurvplot(KM,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           palette = c("#2E9FDF"), 
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           surv.median.line = "h", # 50% 생존지점 표시
           break.time.by = 5,
           xlim = c(0,40)
           )

ggsurvplot(KM,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           palette = c("#2E9FDF"), 
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           surv.median.line = "v", # 50% 생존지점 표시
           break.time.by = 5,
           xlim = c(0,40)
           )

생존분석과 함께 표현할 수 있는 risk table의 경우에는 ggsurvplot 내부의 옵션에서 수정할 수는 없고, table 속성을 따로 불러야 합니다.

SURVPLOT = ggsurvplot(KM,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           palette = c("#2E9FDF"), 
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           surv.median.line = "hv", # 50% 생존지점 표시
           break.time.by = 5,
           xlim = c(0,40)
           )

SURVPLOT$table = SURVPLOT$table +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT

그래프를 SURVPLOT으로 저장하고 $를 활용해주면 table 속성을 따로 다룰 수가 있습니다. 이 방식을 통해 Table 값을 수정하게 됩니다.
또한 마찬가지로 생존곡선 속성또한 다룰 수가 있습니다. 사실은 이렇게 따로 속성을 빼내어 옵션을 바꾼다면 우리는 그래프를 더 익숙한 ggplot2 옵션 방식으로 바꿀 수 있습니다.

SURVPLOT$plot = SURVPLOT$plot +
  guides(col = FALSE) + # STRATA 범례 삭제
  scale_y_continuous(breaks = seq(0,1,by = 0.1),  
                     labels = paste0(seq(0,100,by = 10),"%")) +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT

반대로, fun = ‘event’ 옵션을 주면 누적 발생확률을 기준으로 그래프를 그릴 수가 있습니다.

SURVPLOT = ggsurvplot(KM,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           palette = c("#2E9FDF"), 
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           surv.median.line = "hv", # 50% 생존지점 표시
           break.time.by = 5,
           xlim = c(0,40),
           fun = 'event'
           )

SURVPLOT$plot = SURVPLOT$plot +
  guides(col = FALSE) + # STRATA 범례 삭제
  scale_y_continuous(breaks = seq(0,1,by = 0.1),  
                     labels = paste0(seq(0,100,by = 10),"%")) +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT$table = SURVPLOT$table +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT

표현방식은 매우 다양하기 때문에, 제가 다 표현할 수 없는 부분도 있습니다. 여러분들의 기호에 맞게 옵션을 추가하면 됩니다.

성별에 따라 다르게 줄 경우

이번에는 성별에 따라 다르게 그려보도록 하겠습니다.

KM2 = survfit(Surv(event = EVENT,time = TIME) ~ GENDER, data = Survival)
summary(KM2)

Call: survfit(formula = Surv(event = EVENT, time = TIME) ~ GENDER, 
    data = Survival)

                GENDER=F 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   11     11       1    0.909  0.0867       0.7541            1
   13      9       1    0.808  0.1225       0.6004            1
   21      6       1    0.673  0.1598       0.4229            1
   27      4       1    0.505  0.1887       0.2428            1
   35      2       1    0.253  0.2020       0.0527            1
   37      1       1    0.000     NaN           NA           NA

                GENDER=M 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     24       1   0.9583  0.0408       0.8816        1.000
   19     17       2   0.8456  0.0831       0.6975        1.000
   24     12       1   0.7751  0.1017       0.5993        1.000
   28     10       2   0.6201  0.1274       0.4145        0.928
   30      8       1   0.5426  0.1330       0.3356        0.877
   31      7       1   0.4651  0.1347       0.2636        0.820
   34      6       3   0.2325  0.1164       0.0872        0.620
   35      3       1   0.1550  0.1001       0.0437        0.550
   36      2       1   0.0775  0.0742       0.0119        0.507
   38      1       1   0.0000     NaN           NA           NA

SURVPLOT2 = ggsurvplot(KM2,
           conf.int = TRUE, # 신뢰구간 표현 여부
           risk.table = TRUE, # 테이블 표시 여부
           risk.table.height = 0.4, # 테이블 높이 설정
           ggtheme = theme_bw(), # 데이터 테마 설정
           font.x = c(10), # x축 제목 크기 설정
           font.y = c(10), # y축 제목 크기 설정
           font.tickslab = c(10), # 축 값 크기 설정
           break.time.by = 5,
           xlim = c(0,40),
           pval = TRUE,
           palette = c("#2E9FDF", "#FC4E07"), # 색상 조절
           legend  = c(0.9,0.8), # legend 위치 조정
           legend.title = "Sex", # legend 제목 조정
           legend.labs = c("Male", "Female"), # legend 내용 수정
           conf.int.style = "step" # 신뢰구간 표현 방식 변경
           )

SURVPLOT2$plot = SURVPLOT2$plot +
  scale_y_continuous(breaks = seq(0,1,by = 0.1),  
                     labels = paste0(seq(0,100,by = 10),"%")) +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT2$table = SURVPLOT2$table +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

SURVPLOT2

생존곡선이 여러개가 추가 될 경우, 신뢰구간이 겹치기 때문에 그래프가 자칫 더러워 보일 수가 있습니다. 그렇기에 conf.int.style = “step” 옵션을 줌으로써 신뢰구간을 점선(계단)형태로 표현하여 깔끔하게 나타낼 수가 있습니다.

3. 생존분석 시각화(COX비례위험모형)

다음으로는 COX비례위험모형을 활용하여 생존분석 시각화를 해보도록 하겠습니다. 회귀모형이라고 생각하면 편합니다. 성별 및 약물치료 여부를 변수로 주어 시각화를 진행하겠습니다.

# STARTA는 GENDER로 설정

COX = coxph(Surv(event = EVENT,time = TIME) ~ strata(GENDER) , data = Survival)
summary(COX)

Call:  coxph(formula = Surv(event = EVENT, time = TIME) ~ strata(GENDER), 
    data = Survival)

Null model
  log likelihood= -34.8395 
  n= 40

COX 비례모형을 시각화하는 방식은 구글에 치면 나오지만, 현재 그 방법은 오류가 뜨기 때문에 써먹을 수가 없습니다.

# 현재 오류가 나는 코드

data(lung)

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)

sex_df <- with(lung,
               data.frame(sex = c(1, 2), 
                          age = rep(mean(age, na.rm = TRUE), 2),
                          ph.ecog = c(1, 1)
               )
)
fit <- survfit(res.cox, newdata = sex_df)

ggsurvplot(fit, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
           risk.table = TRUE,
           palette = "Dark2",
           ggtheme = theme_minimal())

이런 상황이기 때문에, 우리는 직접 생존그래프를 작성해야됩니다. 괜찮습니다. 비례모형결과에 그래프 그리기에 필요한 모든 값이 저장되어 있습니다.

SURVFIT = survfit(COX,data = Survival)
SURVFIT$surv # 생존확률

 [1] 1.00000000 1.00000000 1.00000000 0.91310072 0.81707842 0.81707842
 [7] 0.69164195 0.69164195 0.53865129 0.53865129 0.32670852 0.12018935
[13] 1.00000000 0.95918946 0.95918946 0.95918946 0.95918946 0.95918946
[19] 0.95918946 0.84959952 0.84959952 0.84959952 0.84959952 0.78166929
[25] 0.78166929 0.63290519 0.55853687 0.48418327 0.26133336 0.18725353
[31] 0.11357501 0.04178191

SURVFIT$time # 생존기간

 [1]  2  5  7 11 13 17 21 25 27 34 35 37  1  5  6  8 10 11 14 19 20 22 23 24 26
[26] 28 30 31 34 35 36 38

SURVFIT$strata # STRATA(현재, 성별로 설정)

 F  M 
12 20

SURVFIT$upper # 신뢰구간(상)

 [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
 [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[15] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[22] 1.0000000 1.0000000 1.0000000 1.0000000 0.9292307 0.8808169 0.8264697
[29] 0.6343283 0.5633563 0.4960920 0.4854206

SURVFIT$lower # 신뢰구간(하)

 [1] 1.00000000 1.00000000 1.00000000 0.76407629 0.61668483 0.61668483
 [7] 0.44940755 0.44940755 0.28045228 0.28045228 0.10064915 0.01221459
[13] 1.00000000 0.88397028 0.88397028 0.88397028 0.88397028 0.88397028
[19] 0.88397028 0.70469850 0.70469850 0.70469850 0.70469850 0.60981143
[25] 0.60981143 0.43107593 0.35417515 0.28365642 0.10766527 0.06224105
[31] 0.02600180 0.00359632

DF = data.frame(
  GENDER = rep(c("F","M"),c(SURVFIT$strata[1],SURVFIT$strata[2])),
  TIME = SURVFIT$time,
  SURV = SURVFIT$surv,
  UPPER = SURVFIT$upper,
  LOWER = SURVFIT$lower
)

이렇게 생존분석 결과를 활용해 우리는 그래프를 그리기 위한 재료들을 다시 수집할 수가 있습니다.

ggplot(DF) +
  geom_step(aes(x = TIME, y = SURV, col = GENDER),size = 1.2) +
  geom_step(aes(x = TIME, y = LOWER, col = GENDER),linetype = 'dashed',size = 1) +
  geom_step(aes(x = TIME, y = UPPER, col = GENDER),linetype = 'dashed',size = 1) +
  geom_point(aes(x = TIME, y = SURV , col = GENDER), size = 3) +
  scale_x_continuous(breaks = seq(0,40,5)) +
  scale_y_continuous(breaks = seq(0,1,by = 0.1),  
                     labels = paste0(seq(0,100,by = 10),"%")) +
  scale_color_manual(values = c("#2E9FDF", "#FC4E07"),labels = c("Female","Male")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.position = c(0.2,0.2),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15))

ggplot2를 활용하면 굳이 ggsurvplot을 활용하지 않아도 생존그래프를 문제 없이 그릴 수가 있습니다.

여기까지 생존분석 시각화에 대해 다루었으며, 다음 생존분석 그래프에 대한 주제는 생존분석을 더 괴상하게 만들어보고자 합니다.

저작자표시 비영리 변경금지

'R ggplot2' 카테고리의 다른 글

생존분석 그래프 그리기, 생존율 지점 표시(R ggpot2) (0)	2020.02.03

'R ggplot2' Related Articles

생존분석 그래프 그리기, 생존율 지점 표시(R ggpot2) 2020.02.03

Comments

Must Learning With Statistics

생존분석 그래프(ggplot2, ggsurvminer) 그리기 본문

생존분석 그래프(ggplot2, ggsurvminer) 그리기

1. 생존분석 그래프의 중요성

2. 생존분석 시각화(카플란 마이어)

3. 생존분석 시각화(COX비례위험모형)

'R ggplot2' 카테고리의 다른 글

티스토리툴바