今天被粉絲發(fā)的文章給難住了，又偷偷去學(xué)習(xí)了一下競爭風(fēng)險模型，想起之前寫的關(guān)于競爭風(fēng)險模型的做法，真的都是皮毛喲，大家見笑了。想著就順便把所有的生存分析的知識和R語言的做法和論文報告方法都給大家梳理一遍。

什么時候用生存分析

當(dāng)你關(guān)心結(jié)局和結(jié)局發(fā)生時間的時候，就要考慮生存分析了，這種既有結(jié)局又有時間的數(shù)據(jù)叫做生存數(shù)據(jù)，英文叫做Time-to-event data. 只不過因?yàn)檫@個方法醫(yī)學(xué)上用來分析存活情況用的多，所以得名生存分析，反正你就記住一個例子，我要研究汽車發(fā)生故障，我也應(yīng)該用生存分析，因?yàn)槲壹汝P(guān)心是不是有故障，我還關(guān)心用了多久（跑了多遠(yuǎn)）才出故障，就是既有time，又有event，Time-to-event data就用生存分析。

基本概念

首先是刪失，對象失訪了，脫落了，出現(xiàn)結(jié)局之前隨訪結(jié)束了，都叫做刪失：

刪失又分為左刪失，區(qū)間刪失和右刪失，圖示如下：

比如我想研究得了A病的人的生存情況，存在的所有可能情形為：

第一種，研究的開始的時候有人已經(jīng)有A病，這個時候人家已經(jīng)活了一段時間了，具體多久我不知道，叫做左刪失；

第二種，入組隨訪的時候沒病，中途得了A病死了，什么時候得的，沒記錄下來，叫區(qū)間刪失；

第三種，得了A病，一直活到了研究結(jié)束還沒死，叫做右刪失。

你看，所有的刪失情況造成的后果都是我們沒法準(zhǔn)確估計發(fā)生結(jié)局的時間，這也是其名字刪失的由來，對于這類數(shù)據(jù)就需記錄為刪失數(shù)據(jù)。

生存分析的種類有哪些

具體的種類是為了回答具體的問題，我們做生存分析常常要回答的問題如下：

一個是描述生存情況，一個是比較，再一個就是探究影響因素。

比如我隨訪了很多病人，我就想知道隨著時間變化這群人的生存概率是如何變化的（描述）？我就可以用簡單粗暴的Kaplan-Meier method，又叫乘積極限法，基本思想就是此刻的生存概率等于上一時刻的生存概率乘以此刻的存活率。

比如我手上有如下數(shù)據(jù)：

time是隨訪時間，status是結(jié)局，我就可以寫出如下代碼擬合出總體人群的生存曲線：

fit1 <- survfit(Surv(time,status) ~1, data = mydata)summary（fit1）

并且得到每個時間點(diǎn)，病人生存概率的估計值和標(biāo)準(zhǔn)誤：

如果我還恰好收集了病人的性別，我又想看看男病人和女病人生存情況是不是有不同（比較），我可以對生存曲線分組比較，代碼如下：

fit2<- survfit(Surv(time,status) ~ sex, data = mydata)summary（fit2）

輸出兩組生存曲線比較的log-rank test的統(tǒng)計量：

survdiff(Surv(time,status) ~ sex, data = mydata)

并且我們還可以進(jìn)行復(fù)雜分組的組間比較，大家可以翻看我之前的文章。

但是大家明白，KM法始終是在做單因素分析，而且都是做的分類變量的單因素分析，我們經(jīng)常還會有的需求是考慮各種混雜的情況下去探討影響生存時間的因素，這個時候我們就要用到The Cox regression model了。

模型形式如下：

上面的式子把h0移項(xiàng)到左邊，等號兩邊同時取對數(shù)就成了一個線性模型，和廣義線性模型的理解一樣一樣的，b就是我們的影響因素的系數(shù)，解釋為lnHR的改變量，其為正就是危險因素，為負(fù)就是保護(hù)因素：

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.

再觀測一下上面的式子，我們擬合了預(yù)測因素對風(fēng)險比例（t時刻風(fēng)險比上基礎(chǔ)風(fēng)險）的模型，這個時候暗含的假設(shè)就是就是對于每個人在任何時刻風(fēng)險只有一個常數(shù)，就是在所有預(yù)測因素的作用下，各個時間的風(fēng)險是不變的，這個就叫做Proportional Hazards Assumption, 比例風(fēng)險假設(shè)。

做COX比例風(fēng)險模型的時候都有必要驗(yàn)證這個假設(shè)是不是滿足，具體方法如下：

我們需要先做一個cox模型，擬合cox模型需要用到coxph函數(shù)，代碼如下：

res.cox <- coxph(Surv(time,status) ~ age + sex + wt.loss, data =? lung)res.cox

模型輸出結(jié)果里面會有預(yù)測變量的β值（coef）和其標(biāo)準(zhǔn)誤，有HR（exp(coef)）值，還有預(yù)測變量的顯著性檢驗(yàn)結(jié)果：

我們將這個模型對象直接喂給cox.zph，就可以得到風(fēng)險比例假設(shè)的檢驗(yàn)結(jié)果了：

test.ph<-cox.zph(res.cox)test.ph

可以看到，我們的p值都大于0.05，即都滿足ph假設(shè)。

我們還可以從圖上看，看scaled Schoenfeld residuals是不是和時間獨(dú)立，如果獨(dú)立則滿足ph假設(shè)：

ggcoxzph(test.ph)

有競爭事件時

上面寫了只有一個截距事件時生存分析的不同目的，描述，比較，和影響因素探討，接著再來看存在競爭事件的時候該如何描述，比較，和探討影響因素。

生存分析的另外情況就是競爭風(fēng)險模型，就是time to event的event有多種，一種發(fā)生了另外一種就不可能發(fā)生了，這個就是競爭，比如好多文獻(xiàn)中都會列舉的經(jīng)典例子：就是在造血干細(xì)胞移植人群中，我們關(guān)心其疾病復(fù)發(fā)風(fēng)險，但是有些人因?yàn)橐浦菜懒耍═RM），死了之后肯定也談不上復(fù)發(fā)，如果你把因?yàn)橐浦菜懒说娜硕甲鳛閯h失數(shù)據(jù)，肯定也是不對的，這個里面就會有兩個競爭結(jié)局相關(guān)性造成的問題，此時我們應(yīng)該用競爭風(fēng)險模型來分析。

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.

在比例風(fēng)險中我們的結(jié)局事件只有一個，我們關(guān)心的是哪些因素對事件發(fā)生的風(fēng)險有影響。在競爭風(fēng)險模型中我們關(guān)注的地方又變了，因?yàn)槲覀冇泻脦讉€結(jié)局事件，這個時候我們會關(guān)心在競爭事件存在的情況下各個結(jié)局事件的累計發(fā)病函數(shù)是如何隨著時間變化的，以及如何來比較不同的累計發(fā)病函數(shù)，以及如何探討影響因素（competing risks regression analysis）。

之前寫的在非競爭風(fēng)險模型中累計發(fā)病率的組間比較可以用KM法，將縱坐標(biāo)換一換，用log-rank檢驗(yàn)，在競爭風(fēng)險模型中我們需要用Fine and Gray提到的方法來做（Gray’s test），如果比如說我發(fā)現(xiàn)兩組（治療組和對照組）的累計發(fā)病風(fēng)險不一樣，我肯定還想探討哪些因素影響累計發(fā)病風(fēng)險，之前是用COX比例風(fēng)險模型做的，在競爭結(jié)局存在的情況下我們依然是得用Fine and Gray提出的模型（Fine and Gray Model），叫做競爭風(fēng)險回歸模型：

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.

首先看競爭風(fēng)險時候累計發(fā)病曲線的比較方法。我手上有數(shù)據(jù)如下：

其中dis是疾病類型，2分類，ftime是時間，status是結(jié)局事件，結(jié)局事件有3個水平，多的1個水平是競爭事件?，F(xiàn)在我關(guān)心不同疾病人群各個事件累積發(fā)病曲線有無不同，我可以用cuminc函數(shù)結(jié)合plot畫出各個組的累計發(fā)病曲線：

fit=cuminc (ftime,status, dis, cencode =0) plot(fit)

fit對象的結(jié)果中還有每條曲線的比較，從比較結(jié)果可以看出兩組在事件2的累積發(fā)病曲線上是有顯著差異的：

上面介紹的方法相當(dāng)于沒有競爭風(fēng)險的時候的KM法，通過上面的方法我們知道有了不同風(fēng)險結(jié)局累計發(fā)病曲線的差異，繼續(xù)我們會繼續(xù)看影響因素，要做的就是競爭風(fēng)險回歸模型了，需要用到的函數(shù)就是crr。

比如我手上有如下數(shù)據(jù)，除了時間，結(jié)局還包括每個病人的像sex，age等等協(xié)變量，我想探討說這些因素是如何影響病人某個結(jié)局的，我就可以寫出一個競爭風(fēng)險回歸模型：

mod1<- crr(ftime,Status,x)summary(mod1)

上面的代碼中x是自變量矩陣，運(yùn)行代碼輸出競爭風(fēng)險回歸模型的結(jié)果如下：

到這兒基本就寫完了所有的生存分析的情況，接著再結(jié)合一篇論文看看報告方法，論文就是下面這篇，是我隨意查的：

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

這個作者用競爭風(fēng)險模型探討了重型新冠進(jìn)程的影響因素，作者報告了每個影響因素的HR和置信區(qū)間，p值，這些都很容易做出來。還報告了固定時間點(diǎn)的累積發(fā)病的置信區(qū)間。在R語言中固定時間點(diǎn)的累計發(fā)病置信區(qū)間可以用CumIncidence這個函數(shù)計算出來：

The CumIncidence() function allows for the pointwise confidence intervals, by simply adding a further argument, level, where we specify the desired confidence level.

比如我想計算第10天各組的累計發(fā)病風(fēng)險的置信區(qū)間，我就可以寫出如下代碼：

CumIncidence(ftime,status,dis,cencode=0,t=10,level=0.95)

輸出結(jié)果如下：

圖中就是各個組在時間為10的時候結(jié)局累計發(fā)病風(fēng)險的置信區(qū)間。

小結(jié)

今天給大家寫了生存分析的做法，具體包括存在競爭事件和不存在競爭事件的情況下，生存曲線的描述，比較，影響因素探討，感謝大家耐心看完，自己的文章都寫的很細(xì)，重要代碼都在原文中，希望大家都可以自己做一做，請轉(zhuǎn)發(fā)本文到朋友圈后私信回復(fù)“數(shù)據(jù)鏈接”獲取所有數(shù)據(jù)和本人收集的學(xué)習(xí)資料。如果對您有用請先記得收藏，再點(diǎn)贊分享。

也歡迎大家的意見和建議，大家想了解什么統(tǒng)計方法都可以在文章下留言，說不定我看見了就會給你寫教程哦，有疑問歡迎私信。