本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程

R小盐准备介绍R语言机器学习与预测模型的学习笔记

你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】

01 C-statistics&C-index计算

C-statistics即为C统计量是ROC曲线下面积。如果这个logistic回归用于预测或诊断，roc曲线下面积可以估计用logistic回归模型进行诊断或预测的能力。
Logistic回归模型根据预测概率绘制ROC曲线，其曲线下面积AUC = C-Statistics
C-index，即一致性指数（index of concordance），通过评估模型预测结果与实际观察结果的符合程度，以评价模型的预测准确性。
所谓一致性，就是把研究对象随机地两两组对，对于一对病人，如果A实际生存时间长于B，且模型预测的A的生存时间也长于B，则称之为预测结果与实际结果一致。
C-index本质是计算预测结果与实际结果一致的情况所占的比例，类似于ROC曲线下面积AUC。理论上，C-index取值范围是[0.5-1]，若预测与实际完全不一致，则C-index = 0.5，也即模型并无比随机好，没啥实际效用。若模型预测结果与实际结果完全一致，则C-index = 1。
Cox回归中的C-Statistics 一般称为C-index。

#C-statistics计算
library(foreign) 
library(rms)
mydata<-read.spss("lweight.sav")
mydata<-as.data.frame(mydata)
head(mydata)
mydata$low <- ifelse(mydata$low =="低出生体重",1,0)
mydata$race1 <- ifelse(mydata$race =="白种人",1,0)
mydata$race2 <- ifelse(mydata$race =="黑种人",1,0)
mydata$race3 <- ifelse(mydata$race =="其他种族",1,0)
attach(mydata)
dd<-datadist(mydata)
options(datadist='dd')
fit1<-lrm(low~age+ftv+ht+lwt+ptl+smoke+ui+race1+race2,data=mydata,x=T,y=T)
fit1 #直接读取模型中Rank Discrim.参数 C
mydata$predvalue<-predict(fit1)
library(ROCR)
pred <- prediction(mydata$predvalue, mydata$low)
perf<- performance(pred,"tpr","fpr")
plot(perf)
abline(0,1)
auc <- performance(pred,"auc")
auc #auc即是C-statistics
somers2(mydata$predvalue, mydata$low) #somers2 {Hmisc}
  se <- x["S.D."]/sqrt(x["n"])
  Low95 <- x["C Index"] - 1.96*se 
  Upper95 <- x["C Index"] + 1.96*se 
  cbind(x["C Index"], Low95, Upper95) 
cindex <- rcorr.cens(data$prediction,data$event)
cindex
print(CstatisticCI(cindex))
# 该方法可以获得Dxy，但是这个跟cox直接出来的Cindex相比要少，是因为Dxy的原因吗？
# COX的计算
# 根据survival包，可以计算Dxy
# 把-值去掉， |Dxy|/2+0.5
02 模型比较方法
似然比检验
  # 方法1
  # model1 vs model2 (直接使用model名称即可)
  anova(model1,model2) 
  # 方法2
  library(rms)
all.X <-data.frame(x.T=data.T, x.N=data.N, x.S=data.S, x.G=data.G, x.V=data.V, x.P=data.P, x.CEA2=data.CEA2, x.CA1992=data.CA1992)
TN.model <- cph(Surv(survival.time,survival.status)~ x.T+x.N, 
                data=all.X, na.action=na.omit )
TNC.model <- cph(Surv(survival.time,survival.status)~ x.T+x.N+x.CEA2, 
                 data=all.X, na.action=na.omit )
TN2TNC <- lrtest(TN.model, TNC.model)
  # 本质上就是对LIKEHOOD进行卡方检验
  # 经过检验发现，这两者的结果是一致的，均可采用。
基于survcomp包计算c.index
library("survival")
library("prodlim")
library("survcomp")
C_index1 <- concordance.index(x=data$model1_prediction, surv.time=data$time, surv.event=data$event,method="noether")
# 往上翻可以直接查看到C_index以及置信区间
C_index1
C_index2 <- concordance.index(x=data$model2_prediction, surv.time=data$time, surv.event=data$event, method="noether")
C_index2
cindex.comp(C_index1, C_index2)
compareC包
nricens包