Survival analysis_1
Cato LEE
2021 1 11
원출처 : https://radhakrishna.typepad.com/summary---survivalanalysis-a-self-learning-text.pdf
Survival Analysis (SA) 은 cacner related study를 하기 위한 전체적인 개론과 코딩을 설명한다.
#INTRO
SA는 time에 따라 event가 발생했는지 여부를 확인하는 분석이다. 암환자를 예시로 들자면 약물이 투여된 후 시간에 지남에 따라 생존할 확률이 어떻게 변화되는지, 위험에 처한 정도는 얼마나인지 분석을 한다. 이 경우 event는 사망이 된다.
시간에 따라 추적하는 연구이기 때문에 필연적으로 추적이 되지 않는 경우가 발생할 수 있다 (censoring).
censoring이 발생하는 경우를 예시로 들어보면,
- 연구가 종료될 때까지 죽음이 발생하지 않는 경우
- 연구가 진행되는 동안 환자에 대한 추적이 중단된 경우
크게 두 가지를 들 수 있고 이를 Right-censoring 이라고 부른다.
(Left-censoring은 실험에서는 거의 발생하지 않으므로 논외로 한다.)
Survival function은 SA에 있어 기초가 되는 두 개의 기둥 중 하나이다. 이를 통해서 시간의 변화에 따른 생존확률을 추정할 수 있기 때문이다. 확률은 연속적이기 때문에 Survival function의 그래프 또한 곡선이 나와야 하지만 연구 데이터는 countable data이다. 즉 10명이 참여한다면 시간에 따라 1명이 사망한 시점에서 전체 생존은 9명이 되는 식이다. 따라서 Survival function의 그래프는 일반적으로 계단식으로 표현이 된다.
Hazard functon은 나머지 기둥 한개이다. 말 그대로 개체가 특정 시간 t 까지 생존한 찰나의 시점에서 이 개체가 처한 위험을 나타내는 함수이다. 이는 Conditional failure rate라고 불리기도 한다. 즉, 비율이라는 것이며 따라서 Hazard function의 값은 항상 양수이다.
\(S(t) = exp(-\int_{0}^{t} h(u)du)\)
SA는 주어진 결과를 토대로 상기 식을 도출 및 해석함으로써 설명변수들이 S(t)에 미치는 영향을 분석하는 것이 목표이다.
library(tidyverse)
library(survival)생존분석의 결과를 표시하는 몇 가지 방식을 살펴보자.
data_leukemia <- read_csv("leukemia.csv", col_names = T)## Parsed with column specification:
## cols(
## id = col_double(),
## group = col_character(),
## time = col_double(),
## failure = col_double(),
## logWBC = col_double()
## )head(data_leukemia)## # A tibble: 6 x 5
## id group time failure logWBC
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 treat 6 1 2.31
## 2 2 treat 6 1 4.06
## 3 3 treat 6 1 3.28
## 4 4 treat 7 1 4.43
## 5 5 treat 10 1 2.96
## 6 6 treat 13 1 2.88# id, 처리군, 생존시간, 사망발생여부, 혈액 내 백혈구 수치로 테이블이 구성되어 있다.
# time과 failure에 대해서 기록되어 있는 것에 주의하자
# 일반적으로 사건의 발생 시 1, 없을 시 0으로 표기한다.
table(data_leukemia$group, data_leukemia$failure)##
## 0 1
## placebo 0 21
## treat 12 9# pacebo군은 모두 사건이 발생했고 treat 군은 12명이 발생하지 않고 9명이 발생했다
data_bladder <- read.table("bladder.txt", header = T)
head(data_bladder)## ID EVENT INTERVAL INTTIME START STOP TX NUM SIZE
## 1 1 0 1 0 0 0 0 1 1
## 2 2 0 1 1 0 1 0 1 3
## 3 3 0 1 4 0 4 0 2 1
## 4 4 0 1 7 0 7 0 1 1
## 5 5 0 1 10 0 10 0 5 1
## 6 6 1 1 6 0 6 0 4 1# event가 사건의 발생이며 start time과 stop time이 따로 표기되어 있다. 생존시간은 INTTIME 변수이다.
subset(data_bladder, ID == 14)## ID EVENT INTERVAL INTTIME START STOP TX NUM SIZE
## 22 14 1 1 3 0 3 0 3 1
## 23 14 1 2 6 3 9 0 3 1
## 24 14 1 3 12 9 21 0 3 1
## 25 14 0 4 2 21 23 0 3 1# 이와 같이 표기하는 이유는 14번과 같이 사건이 죽음이 아닌 경우, 시험 기간 내 여러차례 발생하는 경우도 있기 때문이다.censoring data를 표현하기 위해서는 survival packages 내 Surv 함수를 이용한다.
treat.group <- subset(data_leukemia, group == "treat" )
# 데이터 내 약물처리군만 뽑아낸다.
Surv(time = treat.group$time, event = treat.group$failure)## [1] 6 6 6 7 10 13 16 22 23 6+ 9+ 10+ 11+ 17+ 19+ 20+ 25+ 32+ 32+
## [20] 34+ 35+# censoring된 data는 +가 숫자 옆에 붙는다.survival rate를 구하기 위해서 survfit 함수를 이용하여 모델을 적합한다.
fit_1 <- survfit(Surv(time = time, event = failure, type = "right") ~ 1,
data = treat.group,
conf.type = "none")
summary(fit_1)## Call: survfit(formula = Surv(time = time, event = failure, type = "right") ~
## 1, data = treat.group, conf.type = "none")
##
## time n.risk n.event survival std.err
## 6 21 3 0.857 0.0764
## 7 17 1 0.807 0.0869
## 10 15 1 0.753 0.0963
## 13 12 1 0.690 0.1068
## 16 11 1 0.627 0.1141
## 22 7 1 0.538 0.1282
## 23 6 1 0.448 0.1346survival curve에 대한 그래프는 다음과 같이 그릴 수 있다.
fit_2 <- survfit(Surv(time = time, event = failure, type = "right") ~ group,
data = data_leukemia)
summary(fit_2)## Call: survfit(formula = Surv(time = time, event = failure, type = "right") ~
## group, data = data_leukemia)
##
## group=placebo
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 21 2 0.9048 0.0641 0.78754 1.000
## 2 19 2 0.8095 0.0857 0.65785 0.996
## 3 17 1 0.7619 0.0929 0.59988 0.968
## 4 16 2 0.6667 0.1029 0.49268 0.902
## 5 14 2 0.5714 0.1080 0.39455 0.828
## 8 12 4 0.3810 0.1060 0.22085 0.657
## 11 8 2 0.2857 0.0986 0.14529 0.562
## 12 6 2 0.1905 0.0857 0.07887 0.460
## 15 4 1 0.1429 0.0764 0.05011 0.407
## 17 3 1 0.0952 0.0641 0.02549 0.356
## 22 2 1 0.0476 0.0465 0.00703 0.322
## 23 1 1 0.0000 NaN NA NA
##
## group=treat
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 6 21 3 0.857 0.0764 0.720 1.000
## 7 17 1 0.807 0.0869 0.653 0.996
## 10 15 1 0.753 0.0963 0.586 0.968
## 13 12 1 0.690 0.1068 0.510 0.935
## 16 11 1 0.627 0.1141 0.439 0.896
## 22 7 1 0.538 0.1282 0.337 0.858
## 23 6 1 0.448 0.1346 0.249 0.807# 그룹에 따라 생존함수를 구한다.
plot(fit_2, lty = 1:2)
legend("topright", lty = 1:2, legend = c("placebo", "treat"), bty = "n")
abline(h = 0.5, col = "blue", lty = 4)
treat와 control과 같이 두 개의 그룹으로 나눈 후 logWBC의 영향과 관련되어 두 그룹의 차이가 실제로 존재하는지도 확인해볼 필요가 있다.
data_er_1 <- subset(data_leukemia, group == "treat")
data_er_2 <- subset(data_leukemia, group == "placebo")
# placebo와 treat 그룹으로 분리한다.
plot(density(data_er_1$logWBC), xlim = c(0, 8), ylim = c(0, 0.6), xlab = "", main = "", lty = 2)
par(new = T)
plot(density(data_er_2$logWBC), xlim = c(0, 8), ylim = c(0, 0.6), xlab = " ", main = " ")
legend("topright", legend = c("treated", "control"), lty = c(2, 1), bty = "n" )
par(mfrow = c(1, 1))
# 밀도함수그래프를 그려본다.
# logWBC가 낮은경우 treat와 control이 겹치고 높은 경우 겹치지 않는 것을 볼 수 있다.
# 따라서 logWBC의 값의 높낮이에 따라 분석을 해 볼 필요가 존재한다.
# 그래프 상 중앙 값이 2.5 정도 되므로 이를 기준으로 군을 나누어 분석해본다.
data_er_1 <- data_leukemia[data_leukemia$logWBC <= 2.5, ]
head(data_er_1)## # A tibble: 6 x 5
## id group time failure logWBC
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 treat 6 1 2.31
## 2 8 treat 22 1 2.32
## 3 14 treat 17 0 2.16
## 4 15 treat 19 0 2.05
## 5 16 treat 20 0 2.01
## 6 17 treat 25 0 1.78fit_1 <- survfit(Surv(time = time, event = failure) ~ group, data = data_er_1)
data_er_2 <- data_leukemia[data_leukemia$logWBC > 2.5, ]
head(data_er_2)## # A tibble: 6 x 5
## id group time failure logWBC
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 2 treat 6 1 4.06
## 2 3 treat 6 1 3.28
## 3 4 treat 7 1 4.43
## 4 5 treat 10 1 2.96
## 5 6 treat 13 1 2.88
## 6 7 treat 16 1 3.6fit_2 <- survfit(Surv(time = time, event = failure) ~ group, data = data_er_2)
par(mfrow = c(1, 2))
plot(fit_1, lty = 1:2, main = "low logWBC")
legend("topright", legend = c("treated", "control"), lty = 1:2, bty = "n")
plot(fit_2, lty = 1:2, main = "high logWBC")
legend("topright", legend = c("treated", "control"), lty = 1:2, bty = "n")
# 낮은 logWBC 값을 갖는 실험자들의 경우 속칭 약빨을 더 잘 받은 것처럼 보인다.이처럼 SA는 생존율의 분석을 설명변수에 따라 추정해 볼 수 있다. 앞으로 SA의 기법과 원리 그리고 R로 구현하는 방법에 대해 공유하려고 한다.