(R2)제18강(1.0) 생존분석(계속 Update...)

tongcon 생각 생존분석을 사용한 경험은 15년전 쯤인가? "교육 및 훈련의 경험"이 "직장의 근무기간"에 어떤 영향을 미치는가? 컨설팅한 것이었습니다. 그 당시에는 SAS의 PHREG 프로시져를 이용했습니다. 그 이후로 생존분석을 하지 않았습니다.  그런데 R을 접하면서 생존분석이 많이 활용되는 것을 보았습니다.  이런 면이 오픈소스인 R 의 특징인 것 같습니다. R이 무료이고 관련 자료도 많고... 생존분석 관련 서적을 보거나, R 패키지(예 survival) 내에 있는 데이터를 이용하면 생존분석에 대한 예제 데이터를 쉽게 구할 수 있고 R 패키지를 이용하여 생존분석을 쉽게 할 수 있습니다. 그러나 실제 상황에서는 데이터를 구하기가 어렵고, 자연히 해석하기에도 어려움이 있습니다. 패키지의 여러 옵션들을 주의해서 결과를 구해야 하고, 해석을 할 때에도 많은 주의가 필요합니다. 처음 "censored data"를 접했을 용어가 헷갈렸습니다. censored 와 uncensored 였습니다. 보통 uncensored 하면 뭔가 부족한 데이터인 줄 알았는데 uncensored 는 절단되지 않은 온전한 데이터이고, censored 가 중간에 중단된 온전하지 못한 데이터입니다... 참고문헌: R을 이용한 생존 분석기초(자유아카데미, 김재희) - 강추 ^^^    생존분석에 대해 상세히 설명하고 있습니다. 생존분석에 사용되는 이론은 만만하지 않습니다.    생존분석에 대해 보다 깊히 이해하고 싶은 분에게 권하고 싶습니다. 닥터 배의 술술 보건의학통계(한나래, 배정민)    보건통계 입문편으로 권장할 만 책입니다.    생존분석에 대한 부분도 설명하고 있습니다.

1. 생존분석이란?

생존할 때까지의 시간을 분석하는 기법

생존분석에서 가장 고려되는 데이터들은 주로 절단된 자료(censored data) 입니다.

사건발생여부와 사건 발생할 때 까지의 시간이 주된 관심

* 절단된 자료(censored data) :

어떤 사건이 발생하기 전에 자료 수집이 불가능해진 (여러가지) 경우,

사건 발생하기 전까지의 관찰된 자료

 

생존분석 데이터의 특징 : 중도절단(censored)된 데이터를 포함하고 있음

중도절단의 예:

연구기간이 끝날 때까지 생존

다른 병이나 사고로 사망

이사 등으로 연구 추적 불가

 

생존분석 용어

(1) 생존함수 : 

(사건이 발생할) 특정시점 t 까지 생존할 확률

 

처음 시작은 확률 1, 무한대로 가면 0

비증가함수.... 

연속 생존함수형태와 (데이터로부터 구하는) 비모수적 계단함수 형태가 있음

생존함수 그래프에서 ①각 시간에 따른 생존(확)율 구할 수 있고 ② 생존(확)율에 해당되는 생존시간을 구할 수 있음

 

확률밀도함수는 다음과 같습니다. (이를 생존함수에서 구해 보면)

생존함수 그리기 - 지수분포를 따르는 경우 exp(-λt) 

lmada가 클수록 생존함수가 가파른 것을 볼 수 있ㅅㅂ니다.

lamda1 <- 1  lamda2 <- 0.5  lamda3 <- 2  t <- seq(0,5,length=5)  plot(t,exp(-lamda1*t),ylab="f(t)",type="l",main="지수 생존함수")  lines(t,exp(-lamda2*t),lty=2)  lines(t,exp(-lamda3*t),lty=3)  legend(2,0.8,c("lamda=1","lamda=05","lamda=2"),lty=c(1,2,3))

(2) 위험함수

 

시점 t에서 (생존한 조건하에서) 사건이 발생한 확률 즉 조건부 확률

 

(3) 평균 잔여 수명( 참고문헌: R을 이용한 생존 분석기초(자유아카데미, 김재희))

특정시점(x) 이후의 평균잔여수명은

시점 x 까지 생존하는 조건하에서

x 시간 이후 생존 가능할 시간의 기대값으로 정의

 

(4) 생존분석에 사용되는 모수적 분포( 참고문헌: R을 이용한 생존 분석기초(자유아카데미, 김재희))

 

분포의 종류	확률밀도함수	위험함수 h(t)		생존함수 S(t)	평균
지수분포(exponential dist)	λexp(-λt)	λ (일정)		exp(-λt)	1/λ λ가 커지면 가파른 곡선 평균이 줄어들겠군 !!!
와이불분포(Weibull dist)
감마분포(Gamma dist)
곰페르츠분포(Gomperz dist)
로그-정규분포	생존시간에 로그를 취하고(ln T), 이 변환값이 정규분포를 따를 때
* 나머지는 참고문헌 참조:

지수분포가 가장 기본적인 분포

 

2. 생존분석 기법

생존율을 산출함

셍존율의 계산

(1) 사망환자의 경우, 치료 시작부터 사망할 때까지의 생존시간 고려

(2) 절단자료인 경우, 치료 시작부터 절단할 떄까지의 시간 고려

 

(1) 생명표법(Life Table method) :

        일정한 간격마다 (구간)생존율을 구함

        관측치 갯수가 50이 넘어야 가능

(2) Kaplan-Meier 방법(Kaplan Meier curve analysis):

        사건(또는 사망)이 발생하는 시점마다 구간생존율을 구함

        이를 누적하여 누적생존율을 구함

        작은 표본( < 50)에서도 적용 가능

 

3. Kaplan-Meier 방법(Kaplan Meier curve analysis):

 

 

생존분석 연습

setwd("d:/r_class") install.packages("survival") library(survival) surv_data <- read.csv("7_ovrian_cancer_survival_data.csv") str(surv_data)head(surv_data) head(surv_data)

# 예제 데이터는 "닥터배이 술술  보건 의학 통계" 에서 가져 왔습니다.

   일단 참고도서를 실행해 보고 하나씩 보강해 나가겠습니다.

bae_01

 

# Kaplan-Meier 생존분석 ====

?Surv attach(surv_data) Surv(month,death==1)

생존기간 변수(month)와 Censored 여부 변수(death) 만 있으면 생명표를 만들 수 있겠구나.

+ Censored 를 표시

bae_02

model_sur <- Surv(month,death==1) survfit(model_sur~treatment)

 

bae_03

model_result <- survfit(model_sur~treatment) plot(model_result,lty=1:2,col=c("blue","green"))

생존곡선: 

초록색: 새 치료법의 (누적)생존율

검은색: 표준치료법의 (누적)생존율

* 사건이 발생한 경우에만 단계가 낮아짐

 

* "척" 보면 새 치료법이 좋은 것인 것을 한 눈에 알 수 있습니다.

   이를 통계적을 설명하는 것이 "로그순위법"입니다.

bae_04

# 로그순위법 ====

 

귀무가설: 두 집단의 생존곡선이 같다.

대립가설: 두 집단의 생존곡선이 다르다.

 

이를 어떻게 검정하지?

만약 두 집단의 생존곡선이 같다면 어떤 형태가 될까?

각 구간마다 관찰대상수가 비슷할 것이다.

(물론 관찰대상수는 두 집단의 전체 대상수에 대한 상대 비율을 고려해야 할 것

여기에 비모수통계기법에서 사용된 순위합검정 방법을 응용해 보면 어떨까?

① 두 집단 합쳐서 하나의 전체 집단을 만든다

② 관찰기간 순으로 정렬한다

③ 이 때 중도절단된 항목은 모두 지운다 => 사건 발생한 구간만 남긴다.

④ 집단별로 각 구간의 사망에 대한 기대 빈도를 구한다,

⑤ 가만... 집단이 2개이고, 각 집단별로 실제 빈도와 기대빈도가 존재한다고?
    그러면 카이제곱분포 그리고 자유도 1인 분포를 따르겠구나

     Σ(관찰 사망자수-기대사망자)^2 / 기대사망자수

 

survdiff(model_sur~treatment)

bae_05

detach(surv_data)