臨床研究的新風口

利用機器學習方法建立和驗證預測模型

我打算講一講自己在Journal of Thoracic Disease（JTD）上面發(fā)表的一篇利用機器學習（隨機森林）建立和驗證A型夾層術前破裂預測模型的文章，名為“Predicting in-hospital rupture of type A aortic dissection using Random Forest”。我主要從選題、研究背景、研究思路、數據分析、文章撰寫、回復審稿人等幾個方面來介紹。我的重點放在數據分析和數據可視化上面，并且把全部代碼分享給大家，希望通過本文的介紹大家可以重復論文中的每一個圖表，用在自己的論文當中。我并非專業(yè)統(tǒng)計人員，如果有錯漏，還請大家不吝賜教。

做臨床研究，一個好的選題是成功（高分SCI）的保證?？陀^的講，JTD的影響因子并不是特別高，但是近年來在業(yè)內的影響力持續(xù)上升。我習慣按照PICOS原則來分解或者設計自己的課題，即Patient population（人群），Intervention（干預措施）, Comparison（對照）, Outcomes（結局）, Study Design（研究設計方案）。 本文的“Patient population”為A型夾層患者。其發(fā)病率大概是3/10,000，算是一種發(fā)病率較低的病種。在世界上，年手術量能達到50例就算是大中心了，而阜外醫(yī)院（心血管領域的超級航空母艦）年手術量能達到250~300左右，積累了大量的病例（本研究納入1133例），這是研究的基礎！我們做臨床課題設計一定要根據自己醫(yī)院的實際來出發(fā)，不要天馬行空的胡思亂想，尤其是對于臨床預測類型的文章，對于樣本量要求大。本文作為觀察性研究，沒有“Intervention”。我們把“Outcomes”作為“Comparison”設置，術前A型夾層破裂組（test group）和不破裂組（control group)比較。“Study design”為觀察性病例對照研究。對于大部分臨床研究的初學者，很容易犯眼高手低的錯誤，一上手就要做雙盲隨機多中心臨床試驗，測試某個藥物是否有效或者是否安全，或者某種術式等等。隨機對照試驗固然好，但是所需的資源（財力、物力、人力、時間成本）絕非普通研究者可以獲得。病例對照研究依然是目前最經典、最實用、最容易開展、發(fā)表量最大的研究類型，也是基于“真實世界”的研究，是我等臨床醫(yī)生必須要掌握的看家本領。

主動脈A型夾層異常兇險，急診手術是唯一有效治療方式。我國主動脈A型夾層手術的經典術式為“全弓置換+支架象鼻”。該術式難度大，培養(yǎng)周期長，目前只有少數心血管中心擁有成熟的技術和團隊能夠處理主動脈A型夾層。同時，主動脈A型夾層的發(fā)病和氣候以及溫度相關，以至于某一段時間內出現大量患者，讓有限的醫(yī)療資源更加捉襟見肘。如果同時來了兩名，甚至多名患者（經常發(fā)生這樣的情況），但是我們無法為每一個患者都同時安排手術！如何做出最優(yōu)選擇/給誰先排手術？毫不夸張的講，這是一個關乎生死抉擇的問題。中國醫(yī)學科學院阜外醫(yī)院作為中國最大的心血管中心，積累了大量主動脈夾層診治的經驗。我們在臨床實際中觀察到，有些夾層患者入院后很快就發(fā)生了破裂（死亡），而有些則穩(wěn)定得多。本研究試圖通過機器學習算法（隨機森林）對A型主動脈夾層患者進行特征提取分析，明確主動脈院內破裂的相關因素，并建立網頁預測模型，為臨床分診提供依據。

這是一篇典型的臨床預測模型文章，該類研究的經典模型為Logistic回歸（針對二分類結局變量binary outcomes）或者Cox回歸（針對時間依賴性結局變量time-to-event）。臨床預測模型文章思路為訓練集建立模型，然后在驗證集中驗證模型。訓練集和驗證集的劃分方式有兩種：

（1） 根據年份來劃分，比如2010-2015連續(xù)納入的患者作為訓練集，2016-2018的患者作為驗證集；

（2）隨機劃分，即將所有患者隨機分為訓練集和驗證集，一般比例為70%:30%。

本研究采用的是隨機劃分的方法。如果樣本量太小，也可以不劃分訓練集和驗證集，對同一個樣本同時進行建模和驗證（重抽樣法）。這種驗證叫做內部驗證，相當于一個人既當運動員又當裁判，結論的可靠性和外推性較差，為迫不得已之法，決非上策。證據等級最高的驗證，是多中心的前瞻性驗證。

驗證的方法最常見的有三種: 表明模型區(qū)分度（Discrimination）的一致性指數（C-index）; 表明模型擬合度（Calibration）的擬合曲線或者Hosmer-Lemeshow index（一種判斷模型擬合優(yōu)度的指數）; 還有比較新潮的決策曲線分析（Decision Curve Analysis）。另外，如果要探究的是引入新的變量對于預測模型的改善度還可以用凈重新分類改善指數（NRI）和綜合判別改善指數（IDI）。

本文的設計思路和傳統(tǒng)Logistic模型是一致的，但是為什么要選擇隨機森林呢？原因有三個：

（1）機器學習是各個領域的研究熱點，研究也講究風潮，所以我們就“跟跟風”吧；

（2）Logistic回歸作為老牌的預測模型工具，盡管沒有過時，但是限制點太多（比如要求事件數10-15 events per variable; 對缺失值敏感，變量不能納入太多等），對于目前的大數據處理多少有點捉襟見肘；

（3）有研究證實機器學習，尤其是隨機森林的預測準確性較高（我們的研究結果也顯示隨機森林算法準確性很高，C-index可以在0.9以上?。?。

隨機森林顧名思義，是用隨機的方式建立一個森林，森林里面有很多的決策樹組成，隨機森林的每一棵決策樹之間是沒有關聯的。在得到森林之后，當有一個新的輸入樣本進入的時候，就讓森林中的每一棵決策樹分別進行一下判斷，看看這個樣本應該屬于哪一類（對于分類算法），然后看看哪一類被選擇最多，就預測這個樣本為那一類。我相信如果之前對于隨機森林沒有了解的話讀完上面這一段話一定懵圈。我也就更不打算深入的介紹隨機森林的算法原理了，免得讓大家產生畏懼心理。如果感興趣的話網上有很多好的帖子，油管（Youtube）里也有很好的介紹。我們臨床醫(yī)生學習統(tǒng)計，重要的是“how”（怎么做）,而不是“why”（為什么）。對于學有余力的同學，當然，了解統(tǒng)計原理有助于自己對統(tǒng)計方法的全面理解和提升個人能力，但是大部分人對于數學公式和原理都是本能地抗拒。對于統(tǒng)計學，我們重要的是知道應用的條件，如何用軟件實現，以及結果如何解讀，至于理論原理，就交給統(tǒng)計學家、數學家去研究吧。我的文章統(tǒng)計都是我自己做的，從SPSS，到Stata，再到R，都是在實戰(zhàn)中去學習，網上找帖子，找視頻，完全無視統(tǒng)計學教科書的存在（當然，學習到一定階段，發(fā)現統(tǒng)計學教科書也還是微言大義的。所謂看山還是山?。?。在這種學習中，我很快樂。我知道一個個圖表是怎么做出來的就夠了——知足常足終身不辱，知止常止終身不殆，我們在知識的海洋里應該敢于承認自己所能夠到的邊際，集中精力，穩(wěn)定心態(tài)，否則很容易迷失。我在大學期間對于統(tǒng)計學一直覺得無比高深，也因此對于臨床研究畏首畏尾，那些統(tǒng)計學原理實在太可怕了，對于臨床醫(yī)生來說，如果深陷其中，絕對是災難！

言歸正傳，說回預測模型工具。傳統(tǒng)的回歸模型有一個重大優(yōu)勢，由于它們歷史悠久，所以有很多優(yōu)秀的工具去呈現這些模型結果，比如網頁計算器，甚至Excel計算器，當然還有新晉當紅小生“Nomogram”。而隨機森林的缺陷就在于此，由于其運算過程幾乎像是一個黑匣子，沒有一個好的呈現工具。相當于我告訴你，西瓜很好吃，但是沒有告訴你如何吃，給你個瓜又有什么用呢？我于是聯系了自己高中的好哥們徐達，他是計算機方面的專家，經過一番努力終于實現了把隨機森林模型嵌入到網頁當中（http://47.107.228.109/），輸入患者的相關臨床特征后，該網頁預測器會自動進行判斷，方便臨床使用。如果大家有類似的需求我們可以一起合作。專業(yè)的人干專業(yè)的事兒，眾人拾柴火焰高。

本研究所采用的軟件是R語言，該軟件是近年來臨床研究的“大殺器”。我結合文章的圖片給大家講解文章的實現過程。文章的表格就是破裂組和非破裂組的比較（分為全隊列、訓練集、測試集三個表），采用t-檢驗、卡方檢驗或者非參數檢驗，無特殊。

圖1. 主動脈A型夾層患者院內破裂的放射學特征

我們的研究納入的變量包括患者的人口學特征、臨床表現、化驗、影像學檢查。對于主動脈夾層而言，影像學檢查是核心變量。我們根據自己的臨床經驗，選取了一些可能和主動脈破裂相關的CT征象（見圖1）。

圖2. 隨機森林特征選擇和排序

從圖2開始，我們就進入了機器學習的部分。A圖展示了隨機森林對于變量的重要性排序，綠色的就是和結局相關的重要變量，紅色是被排除的變量，黃色代表狀態(tài)兩可（根據自己的需要可以歸為重要或者不重要變量）。圖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，而非連續(xù)性的概率，所以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）分部分寫：我一般先寫結果，然后是方法、背景，最后寫討論，寫完加參考文獻。我個人覺得這個順序蠻合理的，實戰(zhàn)中也百試不爽。

文章撰寫當中最好能有外國合作者，他們一般都比較認真，細節(jié)挑錯很厲害，關鍵是語言潤色。這篇文章我寫完后拿給了一個外國的小伙伴Mohammad修改，他對于語言的潤色提供了很大的幫助，所以我把他也列為了共同作者。給大家做個參考，我的英文托福110分，在國外留學過一年，英文其實還不錯。不過從我和Mohammad打交道的情況來看，不管自己英文多么好，始終不是母語，他依舊給我挑出了一大堆語病。因此，潤色是有必要的。

本篇文章審稿人幾乎沒有提出任何意見，很順利就被接受了。但不是每一篇文章都這么幸運。我對于回復的經驗是：

（1 ）簡潔，而不是繞來繞去；

（2）回答過程最好能配以圖、表，這是加分項；

（3）審稿人永遠是對的；

（4）每一條回復都要說謝謝。

（1）隨機森林模型，聽起來高大上，實現起來卻很簡單，值得一試；不過要提醒一下，好的統(tǒng)計學方法可以錦上添花，不可以雪中送炭；可以探求，不可以過度迷信。臨床設計依然是一個研究命運的決定因素，所謂一開始就注定了“輸贏”。

（2）善于合作共贏。和其他專業(yè)，尤其是計算機專業(yè)的小伙伴合作是隨機森林模型網頁化的前提；

（3） Lasso回歸作為統(tǒng)計學新秀，可以作為各種回歸模型的交叉驗證手段；

（4） 英語母語者文章潤色很有必要；

（5） 臨床預測模型是臨床研究的新風口，大數據在手，“paper”（文章）我有。