天天看點

R資料分析:生存資料預測模型的建立和評價(二)

作者:Codewar

上篇文章依照jama surgery的一篇文章給大家寫了生存資料預測模型評價的C指數、校準曲線和模型驗證結果的做法,其實生存資料預測模型的評價方法還有很多,本期接着往下看。

Time-dependent ROC

當結局是一個二分類變量的時候,考慮模型性能的兩個名額一個叫靈敏度和特異度,我們希望兩個都大,模型在不同的cutoff點時這兩個名額如何變化的圖示就是ROC曲線。詳細參考之前的邏輯回歸分類器的文章。

上面的說法是沒有考慮時間因素的,認為結局是固定的,這個在大多數情況下都是适用的,就是說一般情況下我們的結局都是測一次或者說結局是不随着時間變化的,比如結局隻有一個點,可以直接用ROC曲線評估模型。

The classical (standard) approach of ROC curve analysis considers event (disease) status and marker value for an individual as fixed over time, however in practice, both the disease status and marker value change over time

但是遇到生存資料,比如随訪的病人這個時間點沒死,随訪時間延長可能就死了,這個時候模型的靈敏度和特異度在各個時間都不一樣,光說ROC就沒有意義,需要加上時間,就叫做time-dependent ROC。

捋一捋生存分析模型判斷生存情況的原理:

模型會給每一個個案預測一個風險x(文獻中叫做個體的marker),t時刻歸類到死亡或者存活需要一個标準c,模型根據這個風險x和标準c的大小比較結果判定個案是否死亡(陽性陰性)。

模型判斷的是否死亡與個案的實際狀況進行對比就有了混淆矩陣和靈敏度特異度。

D i (t)是t時間的實際狀态,t時刻的敏感度和特異度的計算公式就是下圖:

R資料分析:生存資料預測模型的建立和評價(二)

就是說t時刻個案确實是陽性D i (t)=1的同時模型也說我是“陽性”、“死亡”(marker比标準c大)叫做模型的真陽性率也就是敏感度,模型能将陽性個案識别為陽性的能力。

同理sp(c,t)就是特異度,是真陰率,模型能将無病識别為無病的能力。

R資料分析:生存資料預測模型的建立和評價(二)

本質上ROC曲線可以根據靈敏度和特異度兩個名額來繪制的,我們知道了時間依賴的靈敏度和特異度,曲線也就可以做出來了。

上面的過程簡化下就是:時間依賴的ROC就是不同時刻的結局計算出來的靈敏度和1-特異度構成的ROC。理論上生存結局可以有無數個ROC。

研究目的的不同,Time-dependent ROC又可以進行細分:

(1) cumulative/dynamic (C/D),

(2) incident/dynamic (I/D)

(3) incident/static (I/S)

看下圖,ABCDEF6個個案進行随訪,模型的判定是否死亡的标準是c:

R資料分析:生存資料預測模型的建立和評價(二)

這三種ROC的計算就在于如何界定某時刻t的陽性個案,比如(C/D)就是t時刻所有累計的陽性都算數,是以此時算靈敏度的分母按照下圖就是ADE,随着時間的增大,陽性數量會增大,在計算特異度的時候隻考慮DF的機率,是以叫做累計動态cumulative/dynamic (C/D),這種方法也是臨床使用最多的,它可以評價到t時刻為止模型的表現。

第二種叫做發病動态incident/dynamic (I/D),在這種情況下隻有A作為t時刻的實際陽性,BE都不是(在計算的時候BE直接扔掉了),然後CDF當作陰性,在算特異度的時候隻考慮DF的機率;第三種叫做發病靜态incident/static (I/S),這種情況下依然是看t時刻發病,依然隻有A作為陽性,超過某個區間(0,t*)都沒有發病的作為陰性,即圖中的DF作為陰性。這兩種情況文章應用很少,不做重點記憶。

實操

畫時間依賴ROC用到的包叫做timeROC,主要參數如下所示:

R資料分析:生存資料預測模型的建立和評價(二)

timeROC函數接受的第一個參數T是時間,delta參數是生存狀态,删失一定要設定為0。還有marker參數,這個是模型的預測的自變量,我現在有資料如下,想做一個以

R資料分析:生存資料預測模型的建立和評價(二)

bili,chol,albumin為自變量形成的COX模型的時間依賴ROC,并且将時間設定在時間的每個5分位數上,我就可以寫出如下代碼:

timeROC(T=data$time,
                      delta=data$status,marker=data$bili,
                      other_markers=as.matrix(data[,c("chol","albumin")]),
                      cause=1,weighting="marginal",
                      times=quantile(data$time,probs=seq(0.2,0.8,0.1)),
                      iid=TRUE)           

輸出結果如下,并且可以用confint得到AUC的置信區間:

R資料分析:生存資料預測模型的建立和評價(二)

然後就可以出圖了,我們用plot函數設定好time參數(time一定要和timeROC的參數對應起來)就可以直接出圖了如下:

plot(ROC.bili.cox,time=999.2)           

輸出如下,time改動的話圖也會相應改動:

R資料分析:生存資料預測模型的建立和評價(二)

這樣出來的畫其實是很占空間的,因為你每個時間點都得畫個圖。常見的在文獻中都會将幾個時間點的ROC畫在一個圖上,然後加上圖例進行區分,比如下面這篇來自Diabetes Care的時間依賴ROC圖就是這樣表示的:

R資料分析:生存資料預測模型的建立和評價(二)

那麼對于我們的圖,我們也可以進行修改成像文獻中一樣,直接用plot函數将add參數改為TRUE,然後依次畫出各個時間點的ROC就可以,參考代碼及輸出如下:

plot(ROC.bili.cox,time=999.2)
plot(ROC.bili.cox,time=1307.4,add = T,col="blue")
plot(ROC.bili.cox,time=2555.7,add = T,col="green")
legend(.7, .2, legend=c("關注", "公衆号","Codewar"),
       col=c("red", "blue","green"), lty=1, cex=0.7,bty = "n")           
R資料分析:生存資料預測模型的建立和評價(二)

以上就是Time-dependent ROC的所有内容,我們接着往下看決策曲線的做法。

決策曲線

關于預測模型的決策曲線和校準曲線之前我有專門一篇文章做詳細的解釋,那麼具體到生存分析中,我們依然來回憶一下:

Diagnostic and prognostic models are typically evaluated with measures of accuracy that do not address clinical consequences. Decision-analytic techniques allow assessment of clinical outcomes but often require collection of additional information may be cumbersome to apply to models that yield a continuous result.

就是我們知道模型的評估名額,比如對于生存分析來講我們的模型Brier score和C-index表現的再好,其實也隻能告訴你,模型本身不錯。

但是指導臨床決策我們還要考慮受益的問題,比如對于癌症來講不必要的活檢的危害其實是比沒發現癌症的危害小很多的,我們甯願病人多做一個檢查也絕不願意漏診一個癌症。那麼多少個額外活檢的危害抵得上一次漏診的危害,這個值是可變的,也是我們不知道的,但是我希望的的這個值無論如何變,通過我們的模型做出臨床決策都是受益比瞎猜要好的,這個是我們的終極目的。而決策曲線就是看實際情況是不是符合剛剛的終極目的。

内在的邏輯請參看我之前文章。還有JAMA Guide to Statistics and Methods中的相應介紹,貼在下面就不給大家翻譯了:

R資料分析:生存資料預測模型的建立和評價(二)

文獻中常見的決策曲線圖長這樣

R資料分析:生存資料預測模型的建立和評價(二)

上圖是來自jama子刊的一篇文章,它是将兩個模型的dca畫在了一個圖上,我們依然是在R中進行對照複現。

做決策曲線用到的函數是dca函數,具體到cox模型我們需要指定時間點time:

R資料分析:生存資料預測模型的建立和評價(二)

同時經常出錯的地方還在formula參數的寫法,形成模型整體的決策曲線需要寫模型的marker就是模型的預測風險,範圍必須在0-1。這個是很多同學做錯的地方。

比如我對做好的cox模型将time設定為15然後做決策曲線的示例代碼如下:

dca(Surv(status, time) ~ cancerpredmarker, data = mydata, time = 15)           
R資料分析:生存資料預測模型的建立和評價(二)

其中cancerpredmarker是我拟合的cox模型對個案的風險預測值,上圖就顯示了我的預測模型的決策曲線,多數文章也會對多個預測模型的DCA放一起進行比較,比如我做了模型1得到marker命名為向量“關注”,模型2命名為“公衆号Codewar”,下面代碼即刻讓兩個模型的DCA同時顯示:

dca(Surv(status, time) ~ 關注+公衆号Codewar, 
    data = mydata,
    time = 15,
    thresholds = 1:50 / 100)            

形成的圖示如下:

R資料分析:生存資料預測模型的建立和評價(二)

具體圖的各種風格大家還可根據ggplot的文法進行圖裡位置和圖檔主題風格的修飾。此文略過。還要給大家提醒的是如果模型是競争風險做法依然是一樣的,但是一定要将結局編碼為因子同時删失水準編碼為“censor”,其餘做法都一樣。

Competing risks endpoints are handled similarly to survival endpoints. The outcome must be defined as a factor, with the lowest level called "censor", and the other levels defining the events of interest. The dca() function will treat the first outcome listed as the outcome of interest

好啦,生存資料的預測模型的做法和評價系列到這兒就告一段落了,整體還是比較全面的,全部吃透發個jama應該在做法上不成問題了,有idea的話盡快聯系我實施吧。

小結

今天接着前兩期給大家補上了生存資料預測模型的時間依賴ROC和DCA的原理和做法。感謝大家耐心看完,自己的文章都寫的很細,重要代碼都在原文中,希望大家都可以自己做一做,請轉發本文到朋友圈後私信回複“資料連結”擷取所有資料和本人收集的學習資料。如果對您有用請先記得收藏,再點贊分享。

也歡迎大家的意見和建議,大家想了解什麼統計方法都可以在文章下留言,說不定我看見了就會給你寫教程哦,有疑問歡迎私信,有合作意向請直接滴滴我。

如果你是一個大學大學生或研究所學生,如果你正在因為你的統計作業、資料分析、模型建構,科研統計設計等發愁,如果你在使用SPSS, R,Mplus中遇到任何問題,都可以聯系我。因為我可以給您提供最好的,最詳細和耐心的資料分析服務。

如果你對Z檢驗,t檢驗,方差分析,多元方差分析,回歸,卡方檢驗,相關,多水準模型,結構方程模型,中介調節,量表信效度等等統計技巧有任何問題,請私信我,擷取詳細和耐心的指導。

如果你或你的團隊需要專業的科研資料清洗,模組化服務,教學教育訓練需求等等。請聯系我。

If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #Reports, #Composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.

Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??

Then Contact Me. I will solve your Problem...

If You or Your Research Team Need Professional Scientific Data Cleaning, Model Building Services or Statistical Consulting... Please Contact Me.

往期精彩

R資料分析:生存資料的預測模型建立方法與評價

R資料分析:生存分析的列線圖的了解與繪制詳細教程

R資料分析:結合APA格式作圖大法講講ggplot2和ggsci,請收藏

R資料分析:變量間的非線性關系,多項式,樣條回歸和可加模型

Mplus資料分析:性别差異gendergap的相關研究如何做?

R機器學習:分類算法之logistics回歸分類器的原理和實作

R資料分析:PLS結構方程模型介紹,論文報告方法和實際操作

R資料分析:跟随top期刊手把手教你做一個臨床預測模型

R資料分析:Lasso回歸篩選變量建構Cox模型并繪制列線圖

R資料分析:如何用層次聚類分析做“症狀群”,執行個體操練

R資料分析:工具變量回歸與孟德爾随機化,執行個體解析

R資料分析:潛類别軌迹模型LCTM的做法,執行個體解析

R文本挖掘:中文詞雲生成,以2021新年賀詞為例

R機器學習:分類算法之判别分析LDA,QDA的原理與實作

R可視化:plot函數基礎操作,小白教程

R機器學習:重複抽樣在機器學習模型建立過程中的地位了解

R資料分析:用lme4包拟合線性和非線性混合效應模型

R資料分析:如何用mice做多重插補,執行個體解析

R資料分析:孟德爾随機化中介的原理和實操

R資料分析:生存分析的列線圖的了解與繪制詳細教程

R資料分析:cox模型如何做預測,高分文章複現

R資料分析:廣義估計方程式GEE的做法和解釋

R資料分析:潛類别軌迹模型LCTM的做法,執行個體解析

R資料分析:潛變量與降維方法(主成分分析與因子分析)

R資料分析:如何給結構方程畫路徑圖,tidySEM包詳解

R資料分析:自我報告的身高資料的離群值探索

R資料分析:生存分析與有競争事件的生存分析的做法和解釋

R機器學習:樸素貝葉斯與支援向量機的原理與實作

R資料分析:混合效應模型的可視化解釋,再不懂就真沒辦法

R資料分析:如何了解模型中的“控制”,圖例展示

R資料分析:tableone包的詳細使用介紹

R資料分析:如何用lavaan包做結構方程模型,執行個體解析

R機器學習:分類算法之K最鄰進算法(KNN)的原理與實作

R資料分析:潛增長模型LGM的做法和解釋,及其與混合模型對比

R資料分析:論文中的軌迹的做法,潛增長模型和增長混合模型

R資料分析:縱向分類結局的分析-馬爾可夫多态模型的了解與實操

R資料分析:臨床預測模型實操,校準曲線和DCA曲線做法示例

R資料分析:國産新冠口服藥比輝瑞好的文章的統計做法分享

R資料分析:再寫潛在類别分析LCA的做法與解釋

R資料分析:潛在轉化分析LTA的做法和解釋(一)

R機器學習:分類算法之K最鄰進算法(KNN)的原理與實作

R資料分析:互動作用的簡單斜率圖做法及解釋

R資料分析:雙連續變量互動作用的簡單斜率圖作圖及解釋

繼續閱讀