#title 생존분석 [[TableOfContents]] [http://www.kyobobook.co.kr/product/detailViewKor.laf?ejkGb=KOR&mallGb=KOR&barcode=9788973381944&orderClick=LAH&Kc= 생존분석, 김양진 지음, 자유아카데미] 책을 정리함. ==== 개요 ==== 생존 분석의 주 관심 * 생존 함수(survival function)의 추정 * 생존 함수 또는 위험 함수(hazard function)에 영향을 주는 공변량(corvarate) 또는 예측 변수를 찾아내어 그 연관의 정도를 구하고자 함. 중도 절단 자료 * 우중도절단: t시점까지만 자료의 상태를 알 수 있다. * 연구 종료등의 이유로 t시점까지만 자료를 관측함. 모든 자료가 우중도절단자료 * 미리 정해놓은 사건 발생률에 도달함. 모든 자료가 우중도절단자료이지만, 언제 관측이 종료될지 모름. * 임의 우중도절단자료 * 행정상 중도절단(administrative censoring): 연구가 종료될 때까지 사간이 일어나지 않음. * 추적 실패(loss to follow up) * 경쟁 위험 모형(competing risk): 다른 종류의 사건 발생으로 인해 관심있는 사건이 중도절단된 경우 === 옛날꺼 === ==== 절단된 자료 ==== * 살았는지 죽었는지 모르는 자료. * 하지만 살아있을 당시까지의 데이터는 유효함. * 5년 동안 살아있음을 확인했으나 전화번호가 바뀌어 생존확인이 어려운 경우 ==== 생존함수S(t), 위험함수h(t) ==== * 생존함수 - t시점에 살아있을 확률 * 위험함수 - t시점에 살아있지만 곧 죽을 확률 ==== 예제:생존함수 ==== normally * 0=alive * 1=dead * Other choices TRUE/FALSE (TRUE = death) or 1/2 (2=death) For interval censored data * 0=right censored * 1=event at time * 2=left censored * 3=interval censored time은 일자다. status=1이면 죽음, status=0이면 생존이다. gender=1이면 남자, gender=0이면 여자다. {{{ library("RODBC") conn <- odbcConnect("192.168.201.36",uid="id", pwd="pw") x <- sqlQuery(conn, "select * from ods.dbo.v_h") #성별별로 차이가 있는가? out <- survfit(Surv(time, status==1) ~ gender, data=x) plot(out, lty=1:2, col=c("red", "blue")) survdiff(Surv(time, status==1) ~ gender, data=x) }}} attachment:생존분석/s01.png survdiff(Surv(time, status==1) ~ gender, data=x)의 결과는 다음과 같다. survdiff는 생존함수를 얻는다. {{{ > survdiff(Surv(time, status==1) ~ gender, data=x) Call: survdiff(formula = Surv(time, status == 1) ~ gender, data = x) N Observed Expected (O-E)^2/E (O-E)^2/V gender=0 881 804 852 2.70 6.14 gender=1 1142 1054 1006 2.28 6.14 Chisq= 6.1 on 1 degrees of freedom, p= 0.0132 }}} p-value가 0.0132로 유의수준 0.05에서 귀무가설 기각, 대립가설 채택. 즉, 성별별로 차이가 있음을 확인할 수 있다. summary(out)의 결과를 보자. {{{ > summary(out) Call: survfit(formula = Surv(time, status == 1) ~ gender, data = x) gender=0 time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 881 281 0.6810 0.01570 0.6510 0.7125 2 525 34 0.6369 0.01641 0.6056 0.6699 3 491 25 0.6045 0.01680 0.5725 0.6384 4 466 9 0.5928 0.01692 0.5606 0.6269 5 457 13 0.5760 0.01708 0.5435 0.6104 6 444 10 0.5630 0.01718 0.5303 0.5977 7 434 9 0.5513 0.01726 0.5185 0.5862 8 425 8 0.5409 0.01732 0.5080 0.5760 9 417 5 0.5345 0.01735 0.5015 0.5696 10 412 2 0.5319 0.01736 0.4989 0.5670 11 410 8 0.5215 0.01741 0.4885 0.5567 12 402 5 0.5150 0.01743 0.4819 0.5503 13 397 4 0.5098 0.01745 0.4767 0.5452 14 393 5 0.5033 0.01747 0.4702 0.5387 15 388 8 0.4929 0.01749 0.4598 0.5284 16 380 3 0.4891 0.01749 0.4559 0.5246 17 377 4 0.4839 0.01750 0.4508 0.5194 18 373 4 0.4787 0.01750 0.4456 0.5142 19 369 1 0.4774 0.01750 0.4443 0.5129 20 368 2 0.4748 0.01750 0.4417 0.5104 21 366 2 0.4722 0.01750 0.4391 0.5078 22 364 2 0.4696 0.01750 0.4365 0.5052 23 362 3 0.4657 0.01750 0.4326 0.5013 24 359 3 0.4618 0.01750 0.4288 0.4974 . . . gender=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 1142 437 0.6173 0.01438 0.5898 0.6462 2 618 57 0.5604 0.01490 0.5319 0.5904 3 561 24 0.5364 0.01505 0.5077 0.5667 4 537 20 0.5164 0.01514 0.4876 0.5470 5 517 7 0.5095 0.01516 0.4806 0.5400 6 510 10 0.4995 0.01519 0.4706 0.5301 7 500 11 0.4885 0.01521 0.4596 0.5192 8 489 6 0.4825 0.01522 0.4536 0.5133 9 483 11 0.4715 0.01523 0.4426 0.5023 10 472 9 0.4625 0.01523 0.4336 0.4933 11 463 4 0.4585 0.01523 0.4296 0.4894 12 459 11 0.4475 0.01522 0.4187 0.4784 13 448 4 0.4435 0.01522 0.4147 0.4744 14 444 7 0.4365 0.01520 0.4077 0.4674 15 437 6 0.4305 0.01519 0.4018 0.4614 16 431 7 0.4235 0.01517 0.3948 0.4544 17 424 5 0.4186 0.01516 0.3899 0.4493 18 419 9 0.4096 0.01512 0.3810 0.4403 19 410 2 0.4076 0.01512 0.3790 0.4383 20 408 5 0.4026 0.01510 0.3740 0.4333 21 403 5 0.3976 0.01507 0.3691 0.4282 22 398 4 0.3936 0.01505 0.3652 0.4242 . . . }}} 각 시점별로 죽을 확률이 나온다. 여자는 5일차까지 살아있을 확률이 0.5760이고, 남자는 0.5095다. ==== 예제:위험함수 ==== cox regression을 쓴다. {{{ > summary(coxph(Surv(time, status==1) ~ gender, data=x)) Call: coxph(formula = Surv(time, status == 1) ~ gender, data = x) n= 607, number of events= 605 coef exp(coef) se(coef) z Pr(>|z|) gender 0.10346 1.10900 0.08214 1.26 0.208 exp(coef) exp(-coef) lower .95 upper .95 gender 1.109 0.9017 0.9441 1.303 Concordance= 0.527 (se = 0.033 ) Rsquare= 0.003 (max possible= 1 ) Likelihood ratio test= 1.59 on 1 df, p=0.2069 Wald test = 1.59 on 1 df, p=0.2078 Score (logrank) test = 1.59 on 1 df, p=0.2076 }}} 여자에 비해 남자의 순간사망위험(hazard)가 1.10900배 높다. p-value가 0.208로 유의하지는 않다. ==== 예제 ==== {{{ out2 <- survfit(Surv(suvival_time, status==1) ~ mm, data=m) summary(out2) col <- colors()[c(26,31,33,51,76,153)] plot(out2, lty=1:2, col=col) legend("topright", levels(m$mm), col=col, text.col=col) }}} {{{ col <- colors()[c(26,31,33,36,50,62,70,128,142,258,367,477,275)] plot(fit2, col=col, lwd=1) legend("topright", levels(sur$가입월), col=col, text.col=col) }}} ==== interval ==== {{{ library(survival) x <- df[df$genre == "퍼즐",] out <- survfit(Surv(begin_dt, end_dt, event, type="interval") ~ nation, data=df) col <- colors()[c(26,31,33,36,50,62,70,128,142,258,367,477,275)] plot(out, col=col) legend("topright", levels(df$nation), col=col, text.col=col) s <- summary(out) data.frame(nation=gsub("nation=", "", s$strata), time=s$time, n.risk=s$n.risk, n.event=s$n.event, survival=s$surv) }}} ==== AUC(area under the curve) ==== {{{ #AUC install.packages("pracma") require(pracma) out <- survfit(Surv(time=diff, event=event, type="right") ~ nation, data=x) rs <- data.frame(time=summary(out)$time, surv=summary(out)$surv, nation=gsub("nation=", "",summary(out)$strata)) trapz(rs[rs$nation=="중국",]$time, rs[rs$nation=="중국",]$surv) }}} ==== ctree ==== {{{ library(party) library(survival) data("GBSG2", package = "TH.data") fit <- ctree(Surv(time, cens) ~ ., data = GBSG2) plot(fit) }}} ==== survfit 결과를 data.frame 으로 ==== https://github.com/kmiddleton/rexamples/blob/master/qplot_survival.R {{{ createSurvivalFrame <- function(f.survfit){ # initialise frame variable f.frame <- NULL # check if more then one strata if(length(names(f.survfit$strata)) == 0){ # create data.frame with data from survfit f.frame <- data.frame(time=f.survfit$time, n.risk=f.survfit$n.risk, n.event=f.survfit$n.event, n.censor = f.survfit$n.censor, surv=f.survfit$surv, upper=f.survfit$upper, lower=f.survfit$lower) # create first two rows (start at 1) f.start <- data.frame(time=c(0, f.frame$time[1]), n.risk=c(f.survfit$n, f.survfit$n), n.event=c(0,0), n.censor=c(0,0), surv=c(1,1), upper=c(1,1), lower=c(1,1)) # add first row to dataset f.frame <- rbind(f.start, f.frame) # remove temporary data rm(f.start) } else { # create vector for strata identification f.strata <- NULL for(f.i in 1:length(f.survfit$strata)){ # add vector for one strata according to number of rows of strata f.strata <- c(f.strata, rep(names(f.survfit$strata)[f.i], f.survfit$strata[f.i])) } # create data.frame with data from survfit (create column for strata) f.frame <- data.frame(time=f.survfit$time, n.risk=f.survfit$n.risk, n.event=f.survfit$n.event, n.censor = f.survfit$n.censor, surv=f.survfit$surv, upper=f.survfit$upper, lower=f.survfit$lower, strata=factor(f.strata)) # remove temporary data rm(f.strata) # create first two rows (start at 1) for each strata for(f.i in 1:length(f.survfit$strata)){ # take only subset for this strata from data f.subset <- subset(f.frame, strata==names(f.survfit$strata)[f.i]) # create first two rows (time: 0, time of first event) f.start <- data.frame(time=c(0, f.subset$time[1]), n.risk=rep(f.survfit[f.i]$n, 2), n.event=c(0,0), n.censor=c(0,0), surv=c(1,1), upper=c(1,1), lower=c(1,1), strata=rep(names(f.survfit$strata)[f.i],2)) # add first two rows to dataset f.frame <- rbind(f.start, f.frame) # remove temporary data rm(f.start, f.subset) } # reorder data f.frame <- f.frame[order(f.frame$strata, f.frame$time), ] # rename row.names rownames(f.frame) <- NULL } # return frame return(f.frame) } }}} ==== 카플란-마이어 방법 ==== attachment:생존분석/kaplan_meier.xlsx 위 엑셀에서 누적생존율의 합은 restricted mean survival time (RMST) 다. R코드로 하면.. {{{ library("survival") time <- c(3,2,2,1,0,0,0) status <- c(1,1,1,1,1,1,1) fit <- survfit(Surv(time,status) ~ 1) print(fit, print.rmean=TRUE) plot(fit) rmst <- sum(fit$surv) }}} ==== 참고자료 ==== * https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/ * http://www.r-statistics.com/2013/07/creating-good-looking-survival-curves-the-ggsurv-function/ --> ggsurv 이쁘게 그려줌. * http://www.ats.ucla.edu/stat/r/examples/asa/default.htm * [attachment:생존분석/suvival_lifetime.PDF 고객세분화에_기반한_생존분석을_활용한_고객수명예측모델] ---- 뭘 찾다가 보면 다시 이 사이트로 오는 경우가 많습니다. 좋은 자료 공유해 주셔셔 감사합니다. -- 장운호 2014-11-11 10:16:39 ---- ㅋ 그런가요~ -- 이재학 2014-11-12 16:14:52