临床研究的新风口——利用机器学习方法建立和验证预测模型 | 疯狂统计学2.0
2018年盛夏， 26位科研大神作者以“局解”的方式回顾自身SCI论文发表经历，或介绍如何巧用公共数据库，或侧重某一种统计方法的应用。《疯狂统计学》一书由此横空出世，好评如潮。然而，高阶的统计学方法和数据库的利用需要因地制宜，广大科研初学者的迷思更多在于“科研思路从何而来”“如何推进一项SCI论文研究”。
为予广大读者指点迷津，制作能够“快乐做学术”的科研指导图书，AME出版社决定广纳各路SCI第一作者（欢迎广大读者参与作者阵营，投稿方式下拉至文末查阅），分享从开题到结题的SCI发表经验，汇编为《疯狂统计学（第二版）》。下文为新书《疯狂统计学（第二版）》中关于“临床研究新风口——利用机器学习方法建立和验证预测模型”的精彩篇章，请各位读者尽情享阅。
临床研究的新风口
利用机器学习方法建立和验证预测模型
吴进林，中国医学科学院阜外医院
我打算讲一讲自己在Journal of Thoracic Disease（JTD）上面发表的一篇利用机器学习（随机森林）建立和验证A型夹层术前破裂预测模型的文章，名为“Predicting in-hospital rupture of type A aortic dissection using Random Forest”。我主要从选题、研究背景、研究思路、数据分析、文章撰写、回复审稿人等几个方面来介绍。我的重点放在数据分析和数据可视化上面，并且把全部代码分享给大家，希望通过本文的介绍大家可以重复论文中的每一个图表，用在自己的论文当中。我并非专业统计人员，如果有错漏，还请大家不吝赐教。
选题
做临床研究，一个好的选题是成功（高分SCI）的保证。客观的讲，JTD的影响因子并不是特别高，但是近年来在业内的影响力持续上升。我习惯按照PICOS原则来分解或者设计自己的课题，即Patient population（人群），Intervention（干预措施）, Comparison（对照）, Outcomes（结局）, Study Design（研究设计方案）。 本文的“Patient population”为A型夹层患者。其发病率大概是3/10,000，算是一种发病率较低的病种。在世界上，年手术量能达到50例就算是大中心了，而阜外医院（心血管领域的超级航空母舰）年手术量能达到250~300左右，积累了大量的病例（本研究纳入1133例），这是研究的基础！我们做临床课题设计一定要根据自己医院的实际来出发，不要天马行空的胡思乱想，尤其是对于临床预测类型的文章，对于样本量要求大。本文作为观察性研究，没有“Intervention”。我们把“Outcomes”作为“Comparison”设置，术前A型夹层破裂组（test group）和不破裂组（control group)比较。“Study design”为观察性病例对照研究。对于大部分临床研究的初学者，很容易犯眼高手低的错误，一上手就要做双盲随机多中心临床试验，测试某个药物是否有效或者是否安全，或者某种术式等等。随机对照试验固然好，但是所需的资源（财力、物力、人力、时间成本）绝非普通研究者可以获得。病例对照研究依然是目前最经典、最实用、最容易开展、发表量最大的研究类型，也是基于“真实世界”的研究，是我等临床医生必须要掌握的看家本领。
研究背景
主动脉A型夹层异常凶险，急诊手术是唯一有效治疗方式。我国主动脉A型夹层手术的经典术式为“全弓置换+支架象鼻”。该术式难度大，培养周期长，目前只有少数心血管中心拥有成熟的技术和团队能够处理主动脉A型夹层。同时，主动脉A型夹层的发病和气候以及温度相关，以至于某一段时间内出现大量患者，让有限的医疗资源更加捉襟见肘。如果同时来了两名，甚至多名患者（经常发生这样的情况），但是我们无法为每一个患者都同时安排手术！如何做出最优选择/给谁先排手术？毫不夸张的讲，这是一个关乎生死抉择的问题。中国医学科学院阜外医院作为中国最大的心血管中心，积累了大量主动脉夹层诊治的经验。我们在临床实际中观察到，有些夹层患者入院后很快就发生了破裂（死亡），而有些则稳定得多。本研究试图通过机器学习算法（随机森林）对A型主动脉夹层患者进行特征提取分析，明确主动脉院内破裂的相关因素，并建立网页预测模型，为临床分诊提供依据。
研究思路
这是一篇典型的临床预测模型文章，该类研究的经典模型为Logistic回归（针对二分类结局变量binary outcomes）或者Cox回归（针对时间依赖性结局变量time-to-event）。临床预测模型文章思路为训练集建立模型，然后在验证集中验证模型。训练集和验证集的划分方式有两种：
（1） 根据年份来划分，比如2010-2015连续纳入的患者作为训练集，2016-2018的患者作为验证集；
（2）随机划分，即将所有患者随机分为训练集和验证集，一般比例为70%:30%。
本研究采用的是随机划分的方法。如果样本量太小，也可以不划分训练集和验证集，对同一个样本同时进行建模和验证（重抽样法）。这种验证叫做内部验证，相当于一个人既当运动员又当裁判，结论的可靠性和外推性较差，为迫不得已之法，决非上策。证据等级最高的验证，是多中心的前瞻性验证。
验证的方法最常见的有三种: 表明模型区分度（Discrimination）的一致性指数（C-index）; 表明模型拟合度（Calibration）的拟合曲线或者Hosmer-Lemeshow index（一种判断模型拟合优度的指数）; 还有比较新潮的决策曲线分析（Decision Curve Analysis）。另外，如果要探究的是引入新的变量对于预测模型的改善度还可以用净重新分类改善指数（NRI）和综合判别改善指数（IDI）。
本文的设计思路和传统Logistic模型是一致的，但是为什么要选择随机森林呢？原因有三个：
（1）机器学习是各个领域的研究热点，研究也讲究风潮，所以我们就“跟跟风”吧；
（2）Logistic回归作为老牌的预测模型工具，尽管没有过时，但是限制点太多（比如要求事件数10-15 events per variable; 对缺失值敏感，变量不能纳入太多等），对于目前的大数据处理多少有点捉襟见肘；
（3）有研究证实机器学习，尤其是随机森林的预测准确性较高（我们的研究结果也显示随机森林算法准确性很高，C-index可以在0.9以上！）。
随机森林顾名思义，是用随机的方式建立一个森林，森林里面有很多的决策树组成，随机森林的每一棵决策树之间是没有关联的。在得到森林之后，当有一个新的输入样本进入的时候，就让森林中的每一棵决策树分别进行一下判断，看看这个样本应该属于哪一类（对于分类算法），然后看看哪一类被选择最多，就预测这个样本为那一类。我相信如果之前对于随机森林没有了解的话读完上面这一段话一定懵圈。我也就更不打算深入的介绍随机森林的算法原理了，免得让大家产生畏惧心理。如果感兴趣的话网上有很多好的帖子，油管（Youtube）里也有很好的介绍。我们临床医生学习统计，重要的是“how”（怎么做）,而不是“why”（为什么）。对于学有余力的同学，当然，了解统计原理有助于自己对统计方法的全面理解和提升个人能力，但是大部分人对于数学公式和原理都是本能地抗拒。对于统计学，我们重要的是知道应用的条件，如何用软件实现，以及结果如何解读，至于理论原理，就交给统计学家、数学家去研究吧。我的文章统计都是我自己做的，从SPSS，到Stata，再到R，都是在实战中去学习，网上找帖子，找视频，完全无视统计学教科书的存在（当然，学习到一定阶段，发现统计学教科书也还是微言大义的。所谓看山还是山！）。在这种学习中，我很快乐。我知道一个个图表是怎么做出来的就够了——知足常足终身不辱，知止常止终身不殆，我们在知识的海洋里应该敢于承认自己所能够到的边际，集中精力，稳定心态，否则很容易迷失。我在大学期间对于统计学一直觉得无比高深，也因此对于临床研究畏首畏尾，那些统计学原理实在太可怕了，对于临床医生来说，如果深陷其中，绝对是灾难！
言归正传，说回预测模型工具。传统的回归模型有一个重大优势，由于它们历史悠久，所以有很多优秀的工具去呈现这些模型结果，比如网页计算器，甚至Excel计算器，当然还有新晋当红小生“Nomogram”。而随机森林的缺陷就在于此，由于其运算过程几乎像是一个黑匣子，没有一个好的呈现工具。相当于我告诉你，西瓜很好吃，但是没有告诉你如何吃，给你个瓜又有什么用呢？我于是联系了自己高中的好哥们徐达，他是计算机方面的专家，经过一番努力终于实现了把随机森林模型嵌入到网页当中（http://47.107.228.109/），输入患者的相关临床特征后，该网页预测器会自动进行判断，方便临床使用。如果大家有类似的需求我们可以一起合作。专业的人干专业的事儿，众人拾柴火焰高。
数据分析
本研究所采用的软件是R语言，该软件是近年来临床研究的“大杀器”。我结合文章的图片给大家讲解文章的实现过程。文章的表格就是破裂组和非破裂组的比较（分为全队列、训练集、测试集三个表），采用t-检验、卡方检验或者非参数检验，无特殊。
图1. 主动脉A型夹层患者院内破裂的放射学特征
我们的研究纳入的变量包括患者的人口学特征、临床表现、化验、影像学检查。对于主动脉夹层而言，影像学检查是核心变量。我们根据自己的临床经验，选取了一些可能和主动脉破裂相关的CT征象（见图1）。
图2. 随机森林特征选择和排序
从图2开始，我们就进入了机器学习的部分。A图展示了随机森林对于变量的重要性排序，绿色的就是和结局相关的重要变量，红色是被排除的变量，黄色代表状态两可（根据自己的需要可以归为重要或者不重要变量）。图B是变量重要性决定的过程，每一条线代表一个变量，我们的classifier跑了100次。
那么如何画出上面的图片呢？接下来的代码我将带着大家去实现和呈现随机森林。以下代码大家可以直接粘贴到R里面，可自行尝试。如果不懂代码的同学，建议直接跳过，不会影响阅读。
#加载R包
library(Boruta)
library(randomForest)
library(caret)
library(pROC)
library(ggplot2)
#导入数据(网络数据，可以直接加载），并随机拆分数据集mydata为train和test
library(mlbench); data(HouseVotes84)
mydata<-na.omit(HouseVotes84)
mydata<-rbind(mydata,mydata,mydata,mydata,mydata,mydata)
set.seed(123)
ind <- sample(2, nrow(mydata), replace = TRUE, prob = c(0.7, 0.3))
train <- mydata[ind==1,]
test <- mydata[ind==2,]
#利用Boruta函数进行特征选取，即图2
Boruta(Class~.,data=train,doTrace=2)->Bor.train
#random forest with selected variables
getConfirmedFormula(Bor.train) #确认为重要的变量（绿色）
getNonRejectedFormula(Bor.train) #没有被拒绝的变量（绿色+黄色）
print(Bor.train)
plot(Bor.train) #图2-A
plotImpHistory(Bor.train) #图2-B
## Prediction & Confusion Matrix,模型验证
rf <- randomForest(Class~., data=train,
ntree = 300,
mtry = 8,
importance = TRUE,
proximity = TRUE)
#注意，参数的设置需要自己花点功夫研究，此处不赘述
print(rf) #模型参数展示，比较重要的包括Number of trees, Number of variabels tried at each split, OOB estimate of the error rate
#训练集的C-index,敏感度，特异度，阴性预测值，阳性预测值，准确度等
p1 <- predict(rf, train)
confusionMatrix(p1, train$Class)
#测试集的C-index,敏感度，特异度，阴性预测值，阳性预测值，准确度等
p2 <- predict(rf, test)
confusionMatrix(p2, test$Class)
##训练集的ROC曲线（图4-A)
a<-train$Class
a<-factor(a,levels = c('democrat','republican'),labels = c('0','1'))
a=as.numeric(a)
b<-p1
b<-factor(b,levels = c('democrat','republican'),labels = c('0','1'))
b<-as.numeric(b)
r1=data.frame(a,b)
ROC1 <- roc(a~b,r1)
ROC1
plot(ROC1)
plot(ROC1,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c('green','red'),max.auc.polygon=TRUE,auc.polygon.col='gold',print.thres=TRUE)
##测试集的ROC（Figure 4-B)
c<-test$Class
c<-factor(c,levels = c('democrat','republican'),labels = c('0','1'))
c=as.numeric(c)
d<-p2
d<-factor(d,levels = c('democrat','republican'),labels = c('0','1'))
d<-as.numeric(d)
r2=data.frame(c,d)
ROC2 <- roc(c~d,r2)
ROC2
plot(ROC2)
plot(ROC2,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c('green','red'),max.auc.polygon=TRUE,auc.polygon.col='gold',print.thres=TRUE)
图3. ROC在训练和测试数据集中对于评估模型预测准确性评估
##尽管文章中没有展示，我们也可以探索一下calibration plot
cal1 <- calibration(as.factor(a) ~ b, data = r1)
cal2 <- calibration(as.factor(c) ~ d, data = r2)
xyplot(cal1)
xyplot(cal2)
##注意，由于随机森林给出来的预测是0或者1，而非连续性的概率，所以calibration，DCA等分析方法并不十分适合
##我们还可以查看一下变量的重要性排序（前十）,varImpPlot提供了两种方法：Accurary和Gini系数，都可以作为参考
varImpPlot(rf,
sort = T,
n.var = 10,
main = 'Top 10 - Variable Importance')
从以上代码可以看出来，随机森林还是比较容易实现的，关键是对于参数的设定还有模型网页化。
为了验证一下我们用随机森林选择的变量是否有效，我们同时用Lasso回归来进行了变量筛选，看两者是否有大部分交集。
#加载R包
library(glmnet)
library(rms)
library(foreign)
#数据转换，字符-数字
Class <-train$Class
Class<-factor(Class,levels = c('democrat','republican'),labels = c('0','1'))
Class=as.numeric(Class)
train$Class<-Class
#矩阵
x <- model.matrix(train$Class~ 0 + . , data=train)
#查看并去除x里面变量名后面有0的变量,比如x=x[,-c(5)]，这里没有要去除的变量
y <-train$Class
#varibale trace plot, number below the picture indicates lamda, while numbers above the picture indicates numbers of variable left
fit<-glmnet(x,y,alpha=1,family='binomial')
plot(fit, xvar = 'lambda', label = TRUE)
print(fit)
#cross validation of lasso
cv.fit <- cv.glmnet(x,y,alpha=1,family='binomial',type.measure = 'mse')
plot(cv.fit)
abline(v=log(c(cv.fit$lambda.min,cv.fit$lambda.1se)),lty=2)
#如果lamda取最小值时(preferred),选取到的变量如下
cv.fit$lambda.min
Coefficients <- coef(fit, s = cv.fit$lambda.min)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index
#coefficient of varibales being selected
Active.Coefficients
#varibales being selected
row.names(Coefficients)[Active.Index]
#如果lamda取min-1个标准差,选取到的变量如下
cv.fit$lambda.1se
Coefficients <- coef(fit, s = cv.fit$lambda.1se)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index
#coefficient of varibales being selected
Active.Coefficients
#varibales being selected
row.names(Coefficients)[Active.Index]
##注：本文中的测试数据为基于网络数据的多次重复生成，结果显得怪异，但是不影响代码运行和理解。如果代码有问题，我们可以一起探讨。
图4. 使用套索回归的特征选择
跑完Lasso回归就可以得到图4。图4-A, B两图就是Lasso回归的可视化。C图是通过Excel制作的每一个lasso回归选择出来的变量下术前破裂的比例直方图。
有了代码就有了一切，我希望我的分享能给大家科研助力。不懂代码的也没有关系（毕竟代码只是本文很小的一部分），我们可以汲取思路和方法，重要的是拓展眼界。具体操作可以寻求合作，不一定要事必躬亲。
文章撰写和回复审稿人
文章的撰写不想说太多，也说不清楚。大家学了那么多年的作文，也不见得每一个人都能成为作家。论文写作也是一个道理，需要自己多揣摩。我有几个心得跟大家分享下：
（1）多读文献，好词好句摘抄；
（2）文章写作过程中随时随地记录自己的想法，最后汇总起来；不要试图一蹴而就；
（3）分部分写：我一般先写结果，然后是方法、背景，最后写讨论，写完加参考文献。我个人觉得这个顺序蛮合理的，实战中也百试不爽。
文章撰写当中最好能有外国合作者，他们一般都比较认真，细节挑错很厉害，关键是语言润色。这篇文章我写完后拿给了一个外国的小伙伴Mohammad修改，他对于语言的润色提供了很大的帮助，所以我把他也列为了共同作者。给大家做个参考，我的英文托福110分，在国外留学过一年，英文其实还不错。不过从我和Mohammad打交道的情况来看，不管自己英文多么好，始终不是母语，他依旧给我挑出了一大堆语病。因此，润色是有必要的。
本篇文章审稿人几乎没有提出任何意见，很顺利就被接受了。但不是每一篇文章都这么幸运。我对于回复的经验是：
（1 ）简洁，而不是绕来绕去；
（2）回答过程最好能配以图、表，这是加分项；
（3）审稿人永远是对的；
（4）每一条回复都要说谢谢。
启发
（1）随机森林模型，听起来高大上，实现起来却很简单，值得一试；不过要提醒一下，好的统计学方法可以锦上添花，不可以雪中送炭；可以探求，不可以过度迷信。临床设计依然是一个研究命运的决定因素，所谓一开始就注定了“输赢”。
（2）善于合作共赢。和其他专业，尤其是计算机专业的小伙伴合作是随机森林模型网页化的前提；
（3） Lasso回归作为统计学新秀，可以作为各种回归模型的交叉验证手段；
（4） 英语母语者文章润色很有必要；
（5） 临床预测模型是临床研究的新风口，大数据在手，“paper”（文章）我有。
本文作者简介