18.2 생존분석 예제 - 패키지 {survival}, 데이터 retinopath, 함수 Surv

생존분석에 많이 사용되는 패키지로는 survival, KMsurv, rms 가 있습니다.

이중에서 패키지 {survival}에 들어있는 데이터 retinopath 를 사용한 예입니다.

당뇨병 환자의 시력 손실을 지연시키기 위한 레이저치료법의 효과를 확인하기 위해 수집된 자료로

총 197명의 당뇨병 환자 한쪽 안에는 레이저치료를 진행하고 다른 안에는 치료를 진행하지 않은 상태에서

각 안의 시력상실까지의 시간을 관찰한 자료입니다.

# 참고문헌: https://doi.org/10.21561/jor.2018.3.1.1.

 

사용되는 Surv( ) 함수에는 "생존시간"을 나타내는 변수와 "사건 발생 여부"를 나타내는 변수를 사용합니다.

 

install.packages("survival")

library(survival)

data(retinopathy,package='survival')

head(retinopathy)

# id   laser    eye age type    trt futime status risk

#1 5 argon    left 28  adult   1    46.23   0 9

#2 5 argon    left 28  adult   0   46.23    0 9

#3 14 argon right 12 juvenile 1   42.50   0 8

#4 14 argon right 12 juvenile 0   31.30   1 6

#5 16 xenon right 9  juvenile 1   42.27   0 11

#6 16 xenon right 9  juvenile 0   42.27   0 11

 

# 사용되는 Surv( ) 함수에는 "생존시간"을 나타내는 변수와 "사건 발생 여부"를 나타내는 변수를 사용합니다.

# 치료 여부에 따라 이식편의 생존율이 차이가 있는지 ‘survdiff’ 함수를 이용하여 로그순위 검정법을 실행합니다

# 로그순위 검정법

log <- survdiff(Surv(futime, status) ~ trt, data=retinopathy)

log

 

# 생존곡선 개체를 ‘survfit’ 함수로 생성합니다.

km <- survfit (Surv(futime, status) ~ trt, data = retinopathy)

 

#  ‘plot’ 함수를 이용하여 생존곡선을 그립니다.

plot(km,xlab = "Months", ylab ="Survival rate for visual loss", lty = 1:2, bty = "n", cex.lab = 1.4, cex.axis = 1.4,

xlim = c(0, 80), xaxt = "n")

 

axis(side = 1, at = seq(0, 80, 8), labels = seq(0, 80, 8), cex.axis = 1.4)

text(10, 0.2, labels = "p-value of Log-rank test ＜ 0.001", cex = 1.4, adj = 0)

text(42, 0.76, labels = "Laser treatment (trt = 1)", cex = 0.8, adj = 0)

text(42, 0.36, labels = "No treatment eyes (trt = 0)", cex = 0.8, adj = 0)