ggscidca包2.5版本上线，优雅绘制各种决策曲线，添加了xgboost模块

DCA(Decision Curve Analysis)临床决策曲线是一种用于评价诊断模型诊断准确性的方法，在2006年由AndrewVickers博士创建，我们通常判断一个疾病喜欢使用ROC曲线的AUC值来判定模型的准确性，但ROC曲线通常是通过特异度和敏感度来评价，实际临床中我们还应该考虑，假阳性和假阴性对病人带来的影响，因此在DCA曲线中引入了阈概率和净获益的概念。

ggscidca包是我既往编写的专门用于决策曲线的R包，能绘制多种类型的决策曲线，既往2.3版本已经发布很久了，好评不断，根据粉丝建议，目前上线了2.5版本，这个版本修正了一些前面版本的错误，优化了线条颜色，更加美观，增加了自定线条颜色功能，并且增加了xgboost机器学习模块。

下面我来讲介绍演示一下，安装方法

install.packages("ggscidca")

安装过旧版本的，重新安装一次就可以了

先导入数据和R包

library(sciml)
library(scitable)
library(ggscidca)
bc<-read.csv("E:/r/test/demo.csv",sep=',',header=TRUE)
bc <- na.omit(bc)

做机器学习和生存分析的决策曲线，最好不要有缺失值，不然容易出现各种各样的错误。

数据变量很多，我解释几个我等下要用的，HBP：是否发生高血压，结局指标，AGE：年龄，是我们的协变量，BMI肥胖指数，FEV1肺活量指标，WEIGHT体重，“SBP”,“DBP”：收缩压和舒张压。公众号回复：体检数据，可以获得数据。
有些变量用不到，我先精简一下

bc<-bc[,c("HBP","SEX","AGE","FEV1","OCCU","COUGH","EDU")]

使用organizedata2函数把分类变量转成因子，你自己手动转也是可以的

out2<-organizedata2(data = bc,username=username,token=token,explore = T)
data<-out2[["data"]]
allVars<-out2[["allVars"]]
fvars<-out2[["factorvarout"]]

把数据进行7：3拆分

set.seed(123)
tr1<- sample(nrow(data),0.7*nrow(data))##随机无放抽取
data_train <- data[tr1,]#70%数据集
data_test<- data[-tr1,]#30%数据集

建立xgboot模型，使用sciml包建立xgboot模型非常简单，一句话代码

out<-scixgboost(data = data_train,y="HBP",var = allVars,username=username,token=token)
fitxgboot<-out[["fitxgboot"]]

绘制xgboot模型决策曲线，非常简单，就是一句话代码

scidca(fitxgboot,data_train)

大家可以看到，颜色已经被我优化，比之前更加漂亮了，如果你想自定义颜色也是可以的

scidca(fit2,data_train,lincol = c("#E41A1C", "#377EB8", "#4DAF4A"))

如果你想查看验证集的决策曲线

scidca(fitxgboot,data_test)

此外，ggscidca还能绘制很多模型的决策曲线，

例如逻辑回归模型，在旧版本ggscidca包有个小bug，结局变量类型不能是因子，必须是数字类型，目前这个bug也修正了，函数会自动识别类型并转换。

###逻辑回归
fit2<- glm(HBP ~ SEX + AGE + FEV1+ OCCU+COUGH+EDU,
           family = binomial(link = logit), data = data_train)
scidca(fit2,data_train)

支持向量机模型，我们先用scisvm生成一个模型，你用常规方法生成也是一样的

out2<-scisvm(data = data_train,y="HBP",var = allVars,username=username,token=token)
fitsvm<-out2[["fit"]]

绘制决策曲线

scidca(fitsvm,data_train)

随机森林模型，我这里用scirandomForest函数生成，你用常规方法也是一样的

out<-scirandomForest(data=data_train,y="HBP",username=username,token=token,onlygetfit = T)
fitForest<-out[["fit"]]

绘制决策曲线

scidca(fitForest,data_train)

在SCI文章中，咱们经常可以看到多个决策曲线在一起的图形，如下面文章

咱们使用ggscidca包也是能轻易做到的，先把数据组成一个列表

newdata<-list(data_train,data_train,data_train,data_train)

绘图

tcdca(fitxgboot,fitForest,fitsvm,fit2,newdata = newdata,legend.position = c(0.9,0.8))

既往有粉丝问，这个图能不能自定义颜色，在新版本也支持这个功能了，这里要注意一下，要加上决策曲线的all和none，一共设置5条线

tcdca(fitxgboot,fitForest,fitsvm,fit2,newdata = newdata,legend.position = c(0.9,0.8),
      lincol = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33"))

联合scitable的m.sciroc函数可以事半功倍

out<-m.sciroc(fitxgboot,fitForest,fitsvm,fit2,newdata = newdata,legend.position = c(0.7,0.6))
out[["p"]]

另外生存分析和竞争风险模型也是支持，我简单演示一下，就不解释这么多了，基本一样

library(survival)
library(reshape2)
library(ggplot2)
library(foreign)
library(cmprsk)
library(ggscidca)
bc <- read.spss("E:/r/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc$histgrad<-as.factor(bc$histgrad)
bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
bc <- na.omit(bc)

这里我要说一下，生存分析绘制决策曲线不能有缺失值，不然会报错，这个和我ggscidca包无关，主要是survival包生成的概率值和原数据不一样，导致报错。

生成模型

f1<-coxph(Surv(time,status)~er+histgrad+pr+age+ln_yesno,bc)
f2<-coxph(Surv(time,status)~er+histgrad+ln_yesno,bc)
f3<-coxph(Surv(time,status)~ln_yesno,bc)

绘制单个模型决策曲线

##单个
scidca(f2)

多个模型

newdat<-list(bc,bc,bc)
cox.tcdca(f1,f2,f3,newdata=newdat)

更改颜色，这里要注意一下，要加上决策曲线的all和none，一共设置5条线

cox.tcdca(f1,f2,f3,newdata=newdat,lincol =c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00") )

最后再简单介绍一下竞争风险模型，具体可以看相关文章和视频，这里只上代码

#########竞争风险模型
df_surv<-ggscidca::df_surv
df_surv$cancer_cr<-ifelse(df_surv$cancer_cr=="diagnosed with cancer",1,
                          ifelse(df_surv$cancer_cr=="dead other causes",2,0))
cox_model <- coxph(Surv(ttcancer, cancer_cr==1) ~ age + famhistory + marker, data = df_surv)
cox_model2 <- coxph(Surv(ttcancer, cancer_cr==1) ~ age + famhistory , data = df_surv)
cox_model3 <- coxph(Surv(ttcancer, cancer_cr==1) ~ age  , data = df_surv)
f1<-newcrr(cox_model)
f2<-newcrr(cox_model2)
f3<-newcrr(cox_model3)
newdat<-list(df_surv,df_surv,df_surv)  
cox.tcdca(f1,f2,f3,newdata=newdat)

更改颜色

cox.tcdca(f1,f2,f3,newdata=newdat,lincol =c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00") )

总之就是更好用更方便了，更加详细的操作看相关文章或者下面视频，有问题和建议可以私信给我，本期结束。下面视频更加详细