大家好,，我是阿琛。在之前的推文中,，對于一個分子或者新的風險模型,，我給大家介紹了如何將新指標與臨床病理特征相結合,，通過建立Nomogram圖來可視化地展示患者的預后情況,。但是，當臨床預測模型成功建立,，就會面臨一個問題,，如何去評價模型的質(zhì)量。模型的評價一般是通過比較人群的模型預測結果與實際觀測結果,，來評價預測模型的效果,。
1. 區(qū)分度評價指標：C指數(shù)（C-Index）,，重新分類指數(shù)（Net reclassification index,，NRI）；2. 一致性評價指標：校正曲線（Calibration plot）,；3. 臨床有效性評價指標：決策分析曲線（Decision Curve Analysis,， DCA）；今天就讓阿琛來給大家分別講講如何評價前期建立的臨床預測模型,。首先,，我們一起快速的回顧一下如何構建Cox回歸模型。#install.packages("rms")#install.packages("foreign")#install.packages("survival")library(rms)library(foreign)library(survival)
setwd("C:\Users\000\Desktop\13_clinical") #設置工作目錄rt <- read.table("clinical.txt",header=T,sep=" ") #讀取數(shù)據(jù)head(rt) #查看前5行的數(shù)據(jù)
在該數(shù)據(jù)集中,，共包含6個不同的臨床病理特征以及患者的生存狀態(tài)和生存時間,。首先，對數(shù)據(jù)的相關參數(shù)進行整理與設置：rt$gender <- factor(rt$gender,labels=c("F", "M"))rt$stage <- factor(rt$stage,labels=c("Stage1", "Stage2", "Stage3", "Stage4"))rt$T <- factor(rt$T,labels=c("T1", "T2", "T3", "T4"))rt$M <- factor(rt$M,labels=c("M0", "M1"))rt$N <- factor(rt$N,labels=c("N0", "N1", "N2", "N3"))rt$risk <- factor(rt$risk,labels=c("low", "high"))str(rt)ddist <- datadist(rt) #將數(shù)據(jù)打包options(datadist='ddist')units(rt$futime) <- "year" #定義時間的單位
接著,，將所有變量納入模型中,，構建Cox回歸模型：#構建回歸模型f <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk, x=T, y=T, surv=T, data=rt)f #查看模型的相關參數(shù)
自此，我們的Cox回歸模型就構建完成了,。接下來就是從不同的角度來評價模型的質(zhì)量,。validate(f, method="boot", B=1000, dxy=T) #重復模擬1000次rcorrcens(Surv(futime, fustat) ~ predict(f), data = rt)
在此模型中，C-index為0.664,，即1-0.336=0.664,。當然，根據(jù)C指數(shù)而言,，該模型的預測區(qū)分度并不是太高,。一般而言，我們將0.51-0.7認為是低可信度,，0.71-0.9為中等可信度,，> 0.9為高可信度。#計算三年生存率的校準曲線f3 <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk, x=T, y=T, surv=T, time.inc = 3, data=rt) cal3 <- calibrate(f3, cmethod="KM", method="boot", u= 3, m= 90, B= 1000)#一般將患者分成3組,，即m約等于患者總數(shù)/3,，且重復模擬1000次plot(cal3)
#計算五年生存率的校準曲線f5 <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk, x=T, y=T, surv=T, time.inc = 5, data=rt) cal5 <- calibrate(f5, cmethod="KM", method="boot", u= 5, m= 90, B= 1000)plot(cal5)
對于Calibration plot，橫坐標為模型預測得到的患者生存率,，縱坐標為患者實際觀察得到的生存率,，三點之間的聯(lián)系與虛線之間越相近，越能說明模型良好的一致性,。source("stdca.R") #把stdca.R放在之前設定的工作目錄中,，并進行引用rt$three.years.Survival.Probabilitynew = c(1- (summary(survfit(f, newdata = rt),times = 3)$surv)) #計算3年生存率rt$five.years.Survival.Probabilitynew = c(1- (summary(survfit(f, newdata = rt),times = 5)$surv)) #計算5年生存率write.csv(rt, "DCA.Survival.csv")#將患者的3年和5年生存率結果保存
隨后，根據(jù)患者的生存情況,，分別繪制3年和5年的DCA曲線,；#繪制3年生存率的DCA曲線stdca(data=rt, outcome="fustat", ttoutcome="futime", timepoint=3, predictors="three.years.Survival.Probabilitynew", xstop=0.9, smooth=TRUE)
#繪制5年生存率的DCA曲線stdca(data=rt, outcome="fustat", ttoutcome="futime", timepoint=5, predictors="five.years.Survival.Probabilitynew", xstop=0.9,smooth=TRUE)
對于NRI指數(shù)的計算，有2個專門的R包,，分別為nricens和PredictABEL,；其中，nricens計算出來的是絕對的NRI,，PredictABEL則為相對的NRI,，一般在使用過程中選擇其中之一即可,。#install.packages("nricens")library(nricens)
mstd <- cph(Surv(futime, fustat) ~gender + stage +T + M + N, x=T, y=T, surv=T, data=rt) #使用除risk外的變量進行建模mnew <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk, x=T, y=T, surv=T, data=rt) #該模型包含所有的變量
#計算連續(xù)的NRI,，取值為5%時就定義為重分類#3年生存率NRInricens(mdl.std = mstd, mdl.new = mnew, t0 = 3, updown = 'diff',cut = 0.05, niter = 200)#5年生存率NRInricens(mdl.std = mstd, mdl.new = mnew, t0 = 5, updown = 'diff',cut = 0.05, niter = 200)#計算分類變量計算的NRInricens(mdl.std = mstd, mdl.new = mnew, t0 = 5, updown = 'category',cut = c(0.3,0.6), niter = 200)
對于NRI的計算,，其主要是比較新舊兩個模型之間存在區(qū)分度。其中,，updown主要定義了一個樣本的風險是否變動的方式,，category是指分類值，即熟悉的低,、中,、高風險，另有一種diff,，為連續(xù)值,。Cut是判斷風險高低的臨界值。當updown為diff時,，cut只需設置1個值,，比如0.05，即認為當預測的風險在新舊模型中相差5%時,，即被認為是重新分類了,；而當updown為category時，可以對cut設置兩個不同的值,，即0~29%為低風險,，30%~59%為中風險，60%~100%為高風險,。當updown為diff時,，3年和5年的NRI值分別為0.26和0.53，均大于5%,，可以認為預測的風險在新舊模型中發(fā)生了重新分類,；同時，由于做了200次重復取樣,，相當于有200個NRI,，于是NRI就有了標準誤和95%的置信區(qū)間。同理,，可以得到當updown為category時的NRI值,。到此，對于臨床預測模型評價指標的介紹就到此結束了,。大家可以對照著整個分析過程,，使用示例數(shù)據(jù)或者自己的模型數(shù)據(jù),，進行相應的學習~
|
