최근 병원 타과 선생님께서 요청하셔서 정리했던 내용을 업로드합니다.
당연히 예제로 쓰인 엑셀파일의 값은 가상의 값입니다.
ipynb 파일 [링크]
엑셀 파일 [링크]
- 여기의 예제를 이용해서 작성해보았으며, 바꾼 부분은 아래 소스 코드에 주석으로 표시해보았습니다.
- 원문 예제의 경우 linear prediction을 시행하여 얻은 값으로 다시 생존분석에 들어갔으나, 일반적인 경우에는 이미 존재하는 변수로 바로 분석하는 경우가 대부분인 듯 하여, 이 부분을 바꾸었습니다.
- 또한 통계 분석 방법으로 NNE 를 썼으나, Kaplan Meier 용어가 익숙하기 때문에 KM으로 바꿔서 진행하였습니다.
- 참고로 사용한 라이브러리는 pROC가 아니라 survivalROC 입니다.
예제 작성 by Junn (https://junn.net)¶
작성일: 2021.01.02¶
Original example: https://datascienceplus.com/time-dependent-roc-for-survival-prediction-models-in-r/¶
In [2]:
library(survival)
require(dplyr)
library(psych)
require(ggplot2)
require(survminer)
library(data.table)
require(EnvStats)
library(openxlsx)
library(pROC)
In [4]:
library(repr)
In [5]:
df.basic <- read.xlsx("./01_TimeDependentROC.xlsx", sheet = 1, startRow = 1, colNames = TRUE)
In [6]:
head(df.basic)
Step 1:전체 대상자의 event 커브를 그려본다¶
In [8]:
options(repr.plot.width=6, repr.plot.height=5)
surv_obj <- with(df.basic, Surv(day,event==1))
ggsurvplot(survfit(surv_obj ~ 1,
data = df.basic),
risk.table = TRUE,
break.time.by = 365.25) # 365.25 가 1년이다
In [ ]:
In [9]:
surv_obj <- with(df.basic, Surv(day,event==1))
surv_cox <- coxph(surv_obj~ var2,data=df.basic)
summary(surv_cox)
In [ ]:
In [10]:
library(survivalROC)
## Define a helper functio nto evaluate at various t
survivalROC_helper <- function(t) {
survivalROC(Stime = df.basic$day,
status = (df.basic$event == 1),
marker = df.basic$var2,
predict.time = t,
method = "KM") #NNE
#span = 0.25 * nrow(df.basic)^(-0.20))
}
아래의 예제는 기간(t)만 변경했다¶
In [12]:
library(tidyverse)
## Evaluate every 180 days
survivalROC_data <- data_frame(t = 365.25 * c(1,2,3,4,5,6,7,8,9,10,11)) %>%
mutate(survivalROC = map(t, survivalROC_helper),
## Extract scalar AUC
auc = map_dbl(survivalROC, magrittr::extract2, "AUC"),
## Put cut off dependent values in a data_frame
df_survivalROC = map(survivalROC, function(obj) {
as_data_frame(obj[c("cut.values","TP","FP")])
})) %>%
dplyr::select(-survivalROC) %>%
unnest() %>%
arrange(t, FP, TP)
## Plot
survivalROC_data %>%
ggplot(mapping = aes(x = FP, y = TP)) +
geom_point() +
geom_line() +
geom_label(data = survivalROC_data %>% dplyr::select(t,auc) %>% unique,
mapping = aes(label = sprintf("%.3f", auc)), x = 0.5, y = 0.5) +
facet_wrap( ~ t) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5),
strip.background = element_blank())
Step 4: (테스트) 1년째 (t=365.25)에서 1-TP+FP를 최소화시키는 위치를 찾음 (TP + (1-FP)의 최대화와 동일)¶
In [13]:
time = 365.25
index = with(survivalROC_data %>% filter(t == time), which.min(1-TP+ FP))
# 동일한 결과 with(survivalROC_data %>% filter(t == time), which.max(TP+1-FP))
index
결과를 보니, cut off가 1.0373¶
In [14]:
(survivalROC_data %>% filter(t == 365.25))[index,]
In [ ]:
Step 5: 5년째 부터는 모두 AUC가 0.690, 그래서 일단 5년 event rate을 보자¶
In [15]:
time = 365.25 * 5
index = with(survivalROC_data %>% filter(t == time), which.min(1-TP+ FP))
index
In [16]:
(survivalROC_data %>% filter(t == time))[index,]
Cut-off가 0.7290579로 나온다¶
Step 6: 이 기준으로 그룹을 지어보자¶
In [17]:
df.basic$group = ifelse(df.basic$var2 >= 0.7290579, "Over","Under")
In [18]:
ggsurvplot(survfit(surv_obj ~ group,
data = df.basic),
risk.table = TRUE,
break.time.by = 365.25)
In [ ]:
