今天被粉絲發的文章給難住了，又偷偷去學習了一下競爭風險模型，想起之前寫的關於競爭風險模型的做法，真的都是皮毛喲，大家見笑了。想着就順便把所有的生存分析的知識和R語言的做法和論文報告方法都給大家梳理一遍。
什么時候用生存分析
當你關心結局和結局發生時間的時候，就要考慮生存分析了，這種既有結局又有時間的數據叫做生存數據，英文叫做Time-to-event data. 只不過因為這個方法醫學上用來分析存活情況用的多，所以得名生存分析，反正你就記住一個例子，我要研究汽車發生故障，我也應該用生存分析，因為我既關心是不是有故障，我還關心用了多久（跑了多遠）才出故障，就是既有time，又有event，Time-to-event data就用生存分析。
基本概念
首先是刪失，對象失訪了，脫落了，出現結局之前隨訪結束了，都叫做刪失：
刪失又分為左刪失，區間刪失和右刪失，圖示如下：
比如我想研究得了A病的人的生存情況，存在的所有可能情形為：
第一種，研究的開始的時候有人已經有A病，這個時候人家已經活了一段時間了，具體多久我不知道，叫做左刪失；
第二種，入組隨訪的時候沒病，中途得了A病死了，什么時候得的，沒記錄下來，叫區間刪失；
第三種，得了A病，一直活到了研究結束還沒死，叫做右刪失。
你看，所有的刪失情況造成的后果都是我們沒法准確估計發生結局的時間，這也是其名字刪失的由來，對於這類數據就需記錄為刪失數據。
生存分析的種類有哪些
具體的種類是為了回答具體的問題，我們做生存分析常常要回答的問題如下：
一個是描述生存情況，一個是比較，再一個就是探究影響因素。
比如我隨訪了很多病人，我就想知道隨着時間變化這群人的生存概率是如何變化的（描述）？我就可以用簡單粗暴的Kaplan-Meier method，又叫乘積極限法，基本思想就是此刻的生存概率等於上一時刻的生存概率乘以此刻的存活率。
比如我手上有如下數據：
time是隨訪時間，status是結局，我就可以寫出如下代碼擬合出總體人群的生存曲線：
fit1 <- survfit(Surv(time, status) ~ 1, data = mydata) summary（fit1）
並且得到每個時間點，病人生存概率的估計值和標准誤：
如果我還恰好收集了病人的性別，我又想看看男病人和女病人生存情況是不是有不同（比較），我可以對生存曲線分組比較，代碼如下：
fit2<- survfit(Surv(time, status) ~ sex, data = mydata) summary（fit2）
輸出兩組生存曲線比較的log-rank test的統計量：
survdiff(Surv(time, status) ~ sex, data = mydata)
並且我們還可以進行復雜分組的組間比較，大家可以翻看我之前的文章。
但是大家明白，KM法始終是在做單因素分析，而且都是做的分類變量的單因素分析，我們經常還會有的需求是考慮各種混雜的情況下去探討影響生存時間的因素，這個時候我們就要用到The Cox regression model了。
模型形式如下：
上面的式子把h0移項到左邊，等號兩邊同時取對數就成了一個線性模型，和廣義線性模型的理解一樣一樣的，b就是我們的影響因素的系數，解釋為lnHR的改變量，其為正就是危險因素，為負就是保護因素：
The quantities exp(bi) are called hazard ratios (HR). A value of bi greater than zero, or equivalently a hazard ratio greater than one, indicates that as the value of the ith covariate increases, the event hazard increases and thus the length of survival decreases.
再觀測一下上面的式子，我們擬合了預測因素對風險比例（t時刻風險比上基礎風險）的模型，這個時候暗含的假設就是就是對於每個人在任何時刻風險只有一個常數，就是在所有預測因素的作用下，各個時間的風險是不變的，這個就叫做Proportional Hazards Assumption, 比例風險假設。
做COX比例風險模型的時候都有必要驗證這個假設是不是滿足，具體方法如下：
我們需要先做一個cox模型，擬合cox模型需要用到coxph函數，代碼如下：
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung) res.cox模型輸出結果里面會有預測變量的β值（coef）和其標准誤，有HR（exp(coef)）值，還有預測變量的顯著性檢驗結果：
我們將這個模型對象直接喂給cox.zph，就可以得到風險比例假設的檢驗結果了：
test.ph <- cox.zph(res.cox) test.ph
可以看到，我們的p值都大於0.05，即都滿足ph假設。
我們還可以從圖上看，看scaled Schoenfeld residuals是不是和時間獨立，如果獨立則滿足ph假設：
ggcoxzph(test.ph)
有競爭事件時
上面寫了只有一個截距事件時生存分析的不同目的，描述，比較，和影響因素探討，接着再來看存在競爭事件的時候該如何描述，比較，和探討影響因素。
生存分析的另外情況就是競爭風險模型，就是time to event的event有多種，一種發生了另外一種就不可能發生了，這個就是競爭，比如好多文獻中都會列舉的經典例子：就是在造血干細胞移植人群中，我們關心其疾病復發風險，但是有些人因為移植死了（TRM），死了之后肯定也談不上復發，如果你把因為移植死了的人都作為刪失數據，肯定也是不對的，這個里面就會有兩個競爭結局相關性造成的問題，此時我們應該用競爭風險模型來分析。
For example, disease relapse is an event of interest in studies of allogeneic hematopoietic stem cell transplantation (HSCT), as is mortality related to complications of transplantation (transplant-related mortality or TRM). Relapse and TRM are not independent in this setting because these two events are likely related to immunologic effector mechanisms following HSCT, whereby efforts to reduce TRM may adversely affect the risk of relapse; moreover, patients who die from TRM cannot be at further risk of relapse.
在比例風險中我們的結局事件只有一個，我們關心的是哪些因素對事件發生的風險有影響。在競爭風險模型中我們關注的地方又變了，因為我們有好幾個結局事件，這個時候我們會關心在競爭事件存在的情況下各個結局事件的累計發病函數是如何隨着時間變化的，以及如何來比較不同的累計發病函數，以及如何探討影響因素（competing risks regression analysis）。
之前寫的在非競爭風險模型中累計發病率的組間比較可以用KM法，將縱坐標換一換，用log-rank檢驗，在競爭風險模型中我們需要用Fine and Gray提到的方法來做（Gray’s test），如果比如說我發現兩組（治療組和對照組）的累計發病風險不一樣，我肯定還想探討哪些因素影響累計發病風險，之前是用COX比例風險模型做的，在競爭結局存在的情況下我們依然是得用Fine and Gray提出的模型（Fine and Gray Model），叫做競爭風險回歸模型：
Fine and Gray (6) proposed a method for direct regression modeling of the effect of covariates on the cumulative incidence function for competing risks data. As in any other regression analysis, modeling cumulative incidence functions for competing risks can be used to identify potential prognostic factors for a particular failure in the presence of competing risks, or to assess a prognostic factor of interest after adjusting for other potential risk factors in the model.
首先看競爭風險時候累計發病曲線的比較方法。我手上有數據如下：
其中dis是疾病類型，2分類，ftime是時間，status是結局事件，結局事件有3個水平，多的1個水平是競爭事件。現在我關心不同疾病人群各個事件累積發病曲線有無不同，我可以用cuminc函數結合plot畫出各個組的累計發病曲線：
fit=cuminc (ftime, status, dis, cencode = 0) plot(fit)
fit對象的結果中還有每條曲線的比較，從比較結果可以看出兩組在事件2的累積發病曲線上是有顯著差異的：
上面介紹的方法相當於沒有競爭風險的時候的KM法，通過上面的方法我們知道有了不同風險結局累計發病曲線的差異，繼續我們會繼續看影響因素，要做的就是競爭風險回歸模型了，需要用到的函數就是crr。
比如我手上有如下數據，除了時間，結局還包括每個病人的像sex，age等等協變量，我想探討說這些因素是如何影響病人某個結局的，我就可以寫出一個競爭風險回歸模型：
mod1 <- crr(ftime,Status,x) summary(mod1)上面的代碼中x是自變量矩陣，運行代碼輸出競爭風險回歸模型的結果如下：
到這兒基本就寫完了所有的生存分析的情況，接着再結合一篇論文看看報告方法，論文就是下面這篇，是我隨意查的：
S. Chen, H. Sun, M. Heng, X. Tong, P. Geldsetzer, Z. Wang, P. Wu, J. Yang, Y. Hu, C.
Wang, T. Bärnighausen, Factors Predicting Progression to Severe COVID-19: A Competing Risk Survival Analysis of 1753 Patients in Community Isolation in Wuhan, China, Engineering (2021), doi: https://doi.org/10.1016/j.eng.2021.07.021
這個作者用競爭風險模型探討了重型新冠進程的影響因素，作者報告了每個影響因素的HR和置信區間，p值，這些都很容易做出來。還報告了固定時間點的累積發病的置信區間。在R語言中固定時間點的累計發病置信區間可以用CumIncidence這個函數計算出來：
The CumIncidence() function allows for the pointwise confidence intervals, by simply adding a further argument, level, where we specify the desired confidence level.
比如我想計算第10天各組的累計發病風險的置信區間，我就可以寫出如下代碼：
CumIncidence (ftime, status, dis, cencode = 0,t=10,level = 0.95) 輸出結果如下：
圖中就是各個組在時間為10的時候結局累計發病風險的置信區間。
小結
今天給大家寫了生存分析的做法，具體包括存在競爭事件和不存在競爭事件的情況下，生存曲線的描述，比較，影響因素探討，感謝大家耐心看完，自己的文章都寫的很細，重要代碼都在原文中，希望大家都可以自己做一做，請轉發本文到朋友圈后私信回復“數據鏈接”獲取所有數據和本人收集的學習資料。如果對您有用請先記得收藏，再點贊分享。
也歡迎大家的意見和建議，大家想了解什么統計方法都可以在文章下留言，說不定我看見了就會給你寫教程哦，有疑問歡迎私信。