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

1. 생존 그래프에 구간 표시하기

 이번 포스팅에서는 지난 포스팅에서 그렸던 생존분석 그래프를 더 다양하게 그려보도록 하겠습니다. 생존분석을 진행하는 의사선생님들은 생존율이 90%, 80%, 70% 등 떨어지는 구간을 생존그래프에 표현하고 싶어합니다. 데이터는 저번 포스팅에서 썼던 데이터를 그대로 사용하도록 하겠습니다. 이번 포스팅부터는 survminer패키지는 사용하지 않도록 하겠습니다. ggplot2만으로도 모든 것을 다 그릴 수 있기 때문입니다.

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

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

TIME	EVENT	GENDER	TREATMENT
5	1	F	T
9	1	M	T
2	0	F	N
5	0	F	T
9	0	F	N
9	0	M	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     36       1    0.972  0.0274       0.9200        1.000
    9     33       1    0.943  0.0393       0.8687        1.000
   11     30       1    0.911  0.0490       0.8202        1.000
   12     27       1    0.878  0.0576       0.7716        0.998
   15     25       1    0.842  0.0652       0.7240        0.980
   19     22       1    0.804  0.0726       0.6738        0.960
   20     21       1    0.766  0.0786       0.6264        0.936
   22     20       1    0.728  0.0835       0.5811        0.911
   23     19       1    0.689  0.0874       0.5376        0.884
   25     15       1    0.643  0.0929       0.4848        0.854
   26     14       1    0.597  0.0970       0.4346        0.821
   30     12       1    0.548  0.1008       0.3817        0.786
   32      9       1    0.487  0.1064       0.3171        0.747
   34      7       2    0.348  0.1126       0.1843        0.656
 [ getOption("max.print") 에 도달했습니다 -- 3 행들을 생략합니다 ]

1. 생존분석 그래프 그리기

TIME = KM$time
SURV = KM$surv
UPPER = KM$upper
LOWER = KM$lower

SURV_PLOT = ggplot(NULL) +
  geom_step(aes(x = TIME, y = SURV), col = 'royalblue', size = 1.2) +
  geom_point(aes(x = TIME, y = SURV),col = 'red', size = 3, alpha = 0.5) + 
  geom_step(aes(x = TIME, y = UPPER), linetype = 'dashed',size = 1,alpha = 0.5) +
  geom_step(aes(x = TIME, y = LOWER), linetype = 'dashed',size = 1,alpha = 0.5) +
  theme_classic() +
  scale_x_continuous(expand = c(0,0),limits = c(0,45),
                     breaks = seq(0,50,by = 10)) +
  scale_y_continuous(expand = c(0,0),limits = c(0,1.1),
                     breaks = seq(0,1,by = 0.1)) +
  xlab("Time") + ylab("Surv") + ggtitle("Survival Plot") +
  theme(text = element_text(size = 15, face = "bold")) 

SURV_PLOT

이런 방식으로 ggplot2를 활용해서도 생존그래프를 깔끔하게 그릴 수가 있습니다. 그럼 여기서 이것저것 그래프들을 추가를 해보도록 하겠습니다.

2. 90%, 80%,….,10% 생존율 구간 표시

 생존율이 10% ~ 90%에 근접하였을 때를 그래프에 표현하도록 해보겠습니다. 먼저 각 생존율에 가장 근접한 생존기간을 구해야 합니다. 간단하게 편차의 개념을 이용하면 우리는 이 값들을 계산할 수가 있습니다.

\[ Time=which[min(abs(Surv-0.9))] \]

간단하게 이런 형태의 식을 생각하시면 가장 근접한 생존 기간을 계산할 수가 있습니다. 보고싶은 생존 구간이 0.1 ~ 0.9이니깐 우리는 이걸 for(loop)문을 활용해서 한번에 구하면 될 것 같습니다. 가장 먼저 해주어야 하는 일들은 SURV(생존확률)벡터에 name들을 붙여줘야 합니다.

# name 적용 전
SURV

 [1] 1.0000000 1.0000000 1.0000000 0.9722222 0.9722222 0.9427609 0.9113356
 [8] 0.8775824 0.8424791 0.8424791 0.8041846 0.7658901 0.7275956 0.6893011
[15] 0.6893011 0.6433477 0.5973943 0.5973943 0.5476114 0.5476114 0.4867657
[22] 0.3476898 0.2086139 0.1043069 0.0000000

# name 적용 후
names(SURV) = TIME
SURV

        1         2         4         5         6         9        11 
1.0000000 1.0000000 1.0000000 0.9722222 0.9722222 0.9427609 0.9113356 
       12        15        17        19        20        22        23 
0.8775824 0.8424791 0.8424791 0.8041846 0.7658901 0.7275956 0.6893011 
       24        25        26        27        30        31        32 
0.6893011 0.6433477 0.5973943 0.5973943 0.5476114 0.5476114 0.4867657 
       34        38        39        40 
0.3476898 0.2086139 0.1043069 0.0000000 

이렇게 되면 우리는 각 생존율이 근접하였을때 SURV벡터의 각 원소에 붙은 name(생존기간)을 호출할 수가 있게 됩니다.

R = seq(0.1,0.9, by = 0.1)
R = rev(R)

TIME2 = c()
SURV2 = c()

for(k in 1:length(R)){
  
  TIME2[k] = as.numeric(names(sort(abs(SURV - R[k]))[1]))
  SURV2[k] = SURV[names(SURV) == as.numeric(names(sort(abs(SURV - R[k]))[1]))]
}

TIME2

[1] 11 19 23 26 32 34 34 38 39

SURV2

[1] 0.9113356 0.8041846 0.6893011 0.5973943 0.4867657 0.3476898 0.3476898
[8] 0.2086139 0.1043069

이러면 각 생존율 기준에서 가장 근접한 생존기간을 계산할 수가 있습니다. 여기서 우리가 만들어둔 생존그래프에 간단하게 geom_segment를 추가시키면 구간 별 표시를 할 수 있게 됩니다.

SG_DF = data.frame(
  TIME2 = TIME2,
  SURV2 = SURV2
)

SG_DF$y = 0

SURV_PLOT +
  geom_segment(data = SG_DF,aes(x = TIME2, xend = TIME2, y= y, yend = SURV2),
               linetype = 'dashed', col = 'tan3', size = 1) 

정말 간단하게 생존그래프에 표시를 할 수가 있습니다. 그럼 여기서 더 궁금증이 생길 것이 생존율이 10% ~ 90%에 해당되는 기간을 x축에 표시하고 싶어질 것입니다. 그 부분은 scale_x_continuous()를 다시 활용해주면 말끔하게 수정을 할 수 있습니다.(SURVPLOT에 이미 scale 명령어가 적용되어 있는 상황이지만, 다시 scale을 정해주면 가장 마지막 명령어로 덮어씌이기 때문에 큰 상관이 없습니다.) 여기에 생존율도 그래프 상에 표현을 해주면 좋을 것 같군요.

SURV_PLOT +
  geom_segment(data = SG_DF,aes(x = TIME2, xend = TIME2, y= y, yend = SURV2),
               linetype = 'dashed', col = 'tan3', size = 1)  +
  geom_segment(data = SG_DF,aes(x = 0, xend = TIME2, 
                                y = SURV2, yend = SURV2), 
               col = 'tan3', size = 1, linetype = 'dashed') +
  scale_x_continuous(breaks = SG_DF$TIME2, limits = c(0,45),
                     expand = c(0,0))

이렇게 보니깐, 신뢰구간이랑 너무 겹쳐서 볼품이 없습니다. 다시 깔끔하게 조정하면 다음과 같은 그래프가 생깁니다.

ggplot(NULL) +
  geom_step(aes(x = TIME, y = SURV), col = 'royalblue', size = 1.2) +
  geom_point(aes(x = TIME, y = SURV),col = 'red', size = 3, alpha = 0.5) + 
  geom_segment(data = SG_DF,aes(x = TIME2, xend = TIME2, y= y, yend = SURV2),
               linetype = 'dashed', col = 'tan3', size = 1)  +
  geom_segment(data = SG_DF,aes(x = 0, xend = TIME2, 
                                y = SURV2, yend = SURV2), 
               col = 'tan3', size = 1, linetype = 'dashed') +
  geom_text(data = SG_DF,aes(x = TIME2 + 1, y = SURV2 + 0.05, 
                             label = round(SURV2,2)),
            fontface = 'bold') +
  scale_y_continuous(expand = c(0,0),limits = c(0,1.1),
                     breaks = seq(0,1,by = 0.1)) +
  scale_x_continuous(breaks = SG_DF$TIME2, limits = c(0,45),
                     expand = c(0,0)) +
  xlab("Time") + ylab("Surv") + ggtitle("Survival Plot") +
  theme_classic() +
  theme(text = element_text(size = 15, face = "bold")) 

이렇게 생존그래프를 괴상하게 커스텀마이징을 할 수 있습니다. 다음 포스팅에서는 생존그래프를 여러개 짬뽕시키는 방법과 테이블을 직접 그리는 방법에 대해 다루도록 하겠습니다.