原創博客，未經允許，不得轉載。
生存分析，survival analysis，顧名思義是用來研究個體的存活概率與時間的關系。例如研究病人感染了病毒后，多長時間會死亡；工作的機器多長時間會發生崩潰等。 這里“個體的存活”可以推廣抽象成某些關注的事件。 所以SA就成了研究某一事件與它的發生時間的聯系的方法。這個方法廣泛的用在醫學、生物學等學科上，近年來也越來越多人用在互聯網數據挖掘中，例如用survival analysis去預測信息在社交網絡的傳播程度，或者去預測用戶流失的概率。 R里面有很成熟的SA工具。 本文介紹生存分析的基本概念和一些公式，以及R語言應用示例。
相比於logistics regression， survival analysis有個很大的優點是可以處理缺失數據（censored data， 只有個體的部分數據）。比如說醫生想研究病人在服葯后的健康狀況，跟蹤期為2015一整年。有些人接近2015年底才開始服葯，那么他就只有很短的數據。如果在一年內知道病人死亡了，那么這個病人的數據是完全的；而健康的病人的數據到2015年末就是缺失的了（right censored），因為假如病人在一年后死亡了，我們並不知道這個事。或者說病人在一年內突然失聯了，那么他的數據也是缺失。
用T表示事件的發生時刻
存活函數 ：
S(t)= Pr(T>t) = 1-F(t)
表示個體在t時刻還活着的概率，即事件在t時刻還未發生的概率。S(0)=1, S(無窮大) = 0， 並且S(t)是單調非增函數。實際數據中，t通常取離散值，例如天，周等。
風險函數：
（2） 模型
2.1 如果每個個體都遵循相同的規律，即個體間沒有差異，那么問題比較簡單。Kaplan-Meier是一種無參數的模型，它在每個興趣時間點做一次存活統計，估計存活函數. 令 表示第i時刻還存活的個體數， 表示第i時刻發生事件的個體數量，那么：
這種方法學出來的S(t)是不平滑的。 帶參數的模型會假設模型服從某個分布，使學得的函數平滑。常用的有：
install.packages("OIsurv") library(OIsurv) data(tongue) attach(tongue) my.surv<-Surv(time[type==1], delta[type==1]) my.fit<-survfit(my.surv~1) #Kaplan-Meier summary(my.fit) plot(my.fit) #比較type=1和type=2這兩個組的存活函數 my.fit1<-survfit(Surv(time,delta)~type) plot(my.fit1) #計算風險函數 H.hat<--log(my.fit$surv) H.hat<-c(H.hat, tail(H.hat,1)) print(my.fit, print.rmean=TRUE) #檢驗兩個存活函數是否有區別 survdiff(Surv(time, delta) ~ type) # output omitted detach(tongue)
2.2 Cox proportional hazards model
# cox PH model data(burn) attach(burn) my.surv <- Surv(T1, D1) coxph.fit <- coxph(my.surv ~ Z1 + as.factor(Z11), method="breslow") detach(burn)
預測新數據的值感覺比較不正規，因為survival analysis本身不是針對預測的：
risk = function(model, newdata, time) {
as.numeric(1-summary(survfit(model, newdata = newdata, se.fit = F, conf.int = F), times = time)$surv)
}2.3 extended Cox model
如果Cox PH Model中的變量會隨時間變化，那么就成了extended Cox model，此時HR不再是一個常量。很簡單的例子，如果病人的居住地也是一個變量，病人有可能會搬家，例如在北京吸霾了5年，再跑去廈門生活，那么他舊病復發的概率肯定會降低。所以住所這個變量是和時間相關的。一種簡單的做法是，按照變量改變的時刻，把時間切割成區間，使得每個區間內的變量沒有變化。然后再套用Cox PH模型。
# extended cox data(relapse) N <- dim(relapse)[1] t1 <- rep(0, N+sum(!is.na(relapse$int))) # initialize start time at 0 t2 <- rep(-1, length(t1)) # build vector for end times d <- rep(-1, length(t1)) # whether event was censored g <- rep(-1, length(t1)) # gender covariate i <- rep(FALSE, length(t1)) # initialize intervention at FALSE j <- 1 for(ii in 1:dim(relapse)[1]){ if(is.na(relapse$int[ii])){ # no intervention, copy survival record t2[j] <- relapse$event[ii] d[j] <- relapse$delta[ii] g[j] <- relapse$gender[ii] j <- j+1 } else { # intervention, split records g[j+0:1] <- relapse$gender[ii] # gender is same for each time d[j] <- 0 # no relapse observed pre-intervention d[j+1] <- relapse$delta[ii] # relapse occur post-intervention? i[j+1] <- TRUE # intervention covariate, post-intervention t2[j] <- relapse$int[ii]-1 # end of pre-intervention t1[j+1] <- relapse$int[ii]-1 # start of post-intervention t2[j+1] <- relapse$event[ii] # end of post-intervention j <- j+2 # two records added } } mySurv <- Surv(t1, t2, d) # pg 3 discusses left-trunc. right-cens. data myCPH <- coxph(mySurv ~ g + i)
參考文獻：
http://www.stat.columbia.edu/~madigan/W2025/notes/survival.pdf
http://www.openintro.org/stat/down/Survival-Analysis-in-R.pdf
http://www.ispor.org/congresses/Spain1111/presentations/W29_Ishak-Jack.pdf
http://data.princeton.edu/pop509/ParametricSurvival.pdf