R語言ICU患者死亡率預測實戰

【圖書推薦】《R語言醫學數據分析實踐》-CSDN博客

《R語言醫學數據分析實踐李丹宋立桓蔡偉祺清華大學出版社9787302673484》【摘要書評試讀】- 京東圖書 (jd.com)

預測ICU患者死亡率對比較藥物的療效、比較護理的有效性、比較手術的有效性有重要意義，利用機器學習來構建預測模型，輔助臨床預測有著重要的意義。

6.10.1? 數據集背景

本數據集來自醫學公開數據庫MIMIC III ICU中的數據，有4000患者。納入的指標除了RecordID住院號，年齡age，性別Gender，身高baseWeight ，體重ICU類別baseICUType，SAPS評分，SOFA評分，住院時間 Length_of_stay , Surivial生存天數，院內死亡In.hospital_death。

其他指標處理，根據臨床的相關性選擇納入的指標，這些納入的指標在機器學習中的術語叫特征feature。比如根據臨床經驗，格拉斯哥昏迷評分“GCS”的最低值min及中位數mean可能預測死亡。其它指標也按照“GCS”同樣處理，最后納入的指標有94項，如圖6-19所示。當然某項臨床指標是否納入，如何納入，沒有明確答案，靠臨床經驗和業務知識。

圖6-19

6.10.2? 數據預處理

原始數據中的某些指標在某些患者中沒有觀測值，發現很多缺失值，數據表中以“NA”表示。讀入數據，查看每一列數據概況，R代碼如下所示：

#數據文件icumortality.csv的路徑根據實際情況修改sur<-read.csv("c:/icu_demo/icumortality.csv", header = T, sep = ",")summary(sur)

查看每一列數據概況，summary函數會輸出每一列數據平均值、最大值、最小值以及缺失情況，如圖6-20所示，由于結果輸出太長，這里截圖并不完整。

圖6-20

比如患者的性別Gender 值是0或1，有些患者值為-1，明顯是推斷錯誤，可以改為1。有些患者的height身高值小于100，weight體重小于10，也是推斷不可能，可以賦值NA。計算患者的體重指數BMI，計算完再刪除身高和體重兩列。患者minTemp，meanTemp體溫小于33°，PH值minPH 小于6或maxPH大于8，也是不可能，可以都賦值“NA”。

如圖6-21所示，medianTroponinI缺失值達3795個（90%以上），可以刪除，暫時不用的RecordID，SAPS.I、SOFA、length_of_stay、survival可以刪除。

圖6-21

MechVent指標是機械通氣，0表示無，1表示有，如圖6-22所示，發現部分患者MechVent值為“NA”，推斷沒有機械通氣，可以賦值為0。

圖6-22

血壓有有創血壓和無創血壓兩種，可以先找出有創血壓DiasABP，SysABP 、MAP的缺失值，用此患者的的無創血壓（帶NI的變量）替換，然后刪除無創血壓的數據。

以上是利用經驗判斷缺失值是什么，此外可以使用mice程序包進行插補。mice即是基于多重填補法構造的，基本思想是對于一個具有缺失值的變量，用其他變量的數據對這個變量進行擬合，再用擬合的預測值對這個變量的缺失值進行填補。

數據清洗R代碼如下：

#數據清洗
sur$baseGender[sur$baseGender<0] <- 1
sur$baseHeight[which(sur$baseHeight < 100)] <- NA
sur$baseHeight[which(sur$baseHeight > 210)] <- NA
sur$baseWeight[which(sur$baseWeight < 10)] <- NA
sur$BMI <- (sur$baseWeight/sur$baseHeight)/sur$baseHeight * 10000
sur[,"baseWeight"] <- NULL
sur[,"baseHeight"] <- NULL
sur[,"medianTroponinI"] <- NULLsub=which(is.na(sur$meanDiasABP))
sur$meanDiasABP[sub] <- sur$meanNIDiasABP[sub]
sur$minDiasABP[sub] <- sur$minNIDiasABP[sub]
sur$maxDiasABP[sub] <- sur$maxNIDiasABP[sub]
sur[,"meanNIDiasABP"] <- NULL
sur[,"minNIDiasABP"] <- NULL
sur[,"maxNIDiasABP"] <- NULL
sub=which(is.na(sur$meanMAP))
sur$meanMAP[sub] <- sur$meanNIMAP[sub]
sur$minMAP[sub] <- sur$minNIMAP[sub]
sur$maxMAP[sub] <- sur$maxNIMAP[sub]
sur[,"meanNIMAP"] <- NULL
sur[,"minNIMAP"] <- NULL
sur[,"maxNIMAP"] <- NULL
sub=which(is.na(sur$meanSysABP))
sur$meanSysABP[sub] <- sur$meanNISysABP[sub]
sur$minSysABP[sub] <- sur$minNISysABP[sub]
sur$maxSysABP[sub] <- sur$maxNISysABP[sub]
sur[,"meanNISysABP"] <- NULL
sur[,"minNISysABP"] <- NULL
sur[,"maxNISysABP"] <- NULLsur$minTemp[sur$minTemp < 33] <- NA
sur$meanTemp[sur$meanTemp < 33] <- NAsur$minpH[sur$minpH>8] <- NA
sur$minpH[sur$minpH<6] <- NA
sur$meanpH[sur$meanpH>8] <- NA
sur$maxpH[sur$maxpH>8] <- NAsur$meanMechVent[is.na(sur$meanMechVent)] <- 0sur[,"RecordID"] <- NULL
sur[,"SAPS.I"] <- NULL
sur[,"SOFA"] <- NULL
sur[,"Length_of_stay"] <- NULL
sur[,"Survival"] <- NULL

對于缺失值插補，可以利用mice程序包進行插補：

library(mice)
sur_clean <- mice(sur, m=5, seed=123456)
completedData <- complete(sur_clean,1)

6.10.3? 模型建立

首先我們先用大家熟悉的Logistic回歸分析模型，載入glmnet包，這個R包包括廣義線性回歸方程所有功能，R代碼如圖6-23所示：

圖6-23

R代碼運行，結果如圖6-24所示。

圖6-24

因為我們最后要比較各個機器學習模型的表現，需要新建一個表格用來存儲模型的AUC，新建一個表格final,新表格只有一列label，是真實的死亡情況，新表格增加一列lasso,存儲lasso模型的預測值，R代碼如下所示：

label<-completedData$In.hospital_death
final <- as.data.frame(label)
final$Lasso <- completedData$test_pred

我們接著看支持向量機SVM ，給定一組訓練樣本集，樣本數據集是二維的，分散在平面上，需要找到一條直線將數據集分割開。可以分開的直線有很多，我們要找到其中泛化能力最好，魯棒性最強的直線。這是在平面上的點，如果是在三維空間中，則需要找到一個平面；如果是超過三維以上的維數，則需要找到一個超平面。

SVM的三個重要參數是kernel參數、C參數cost和gamma參數來調節模型的表現：

（1）C參數：C參數（Cost）它是控制著模型的懲罰力度，即決定了模型在訓練時對誤差和復雜度之間的折衷關系。通常情況下，C值越大，模型的復雜度越高，預測準確度也相應提高。但是，如果C值過大，可能會導致模型出現過擬合的情況，因此需要在實際應用中進行試探和調整。

（2）gamma參數：gamma參數決定了樣本點對于模型的影響程度，即越接近樣本點的數據在模型中的權重也越大。通常情況下，gamma值越小，對于遠離樣本點的數據的影響也就越小，而對于靠近樣本點的數據的影響則更加強烈。但是，如果gamma值過小，可能會導致模型出現欠擬合的情況，因此需要進行合理的選擇。

（3）kernel參數：Kernel指的是支持向量機的類型，它可能是線性還是非線性，首先選擇線性的kernel作為基準。核函數是支持向量機回歸算法的核心，不同的核函數適用于不同的數據集。常見的核函數有線性核函數、多項式核函數和徑向基核函數等。在選擇核函數時需要考慮數據集的特點和所需的預測準確度。

C參數即Cost是違反約束時的成本函數，gamma是除線性SVM外其余所有SVM都使用的一個參數，它的選擇與直線或平面的非線性程度相關，gamma選擇越大，其非線性程度越大。

使用支持向量機模型R代碼如圖6-25所示。

代碼運行，結果如圖6-26所示。

圖6-26

在final數據表中存儲SVM模型的AUC

final$SVM <- completedData$test_pred

隨機森林Radom Forest算法介紹：簡單的說，隨機森林就是用隨機的方式建立一個森林，森林里面有很多的決策樹，并且每棵樹之間是沒有關聯的。得到一個森林后，當有一個新的樣本輸入，森林中的每一棵決策樹會分別進行一下判斷，進行類別歸類（針對分類算法），最后比較一下被判定哪一類最多，就預測該樣本為哪一類。在隨機森林方法中，創建了大量的決策樹。每個觀察結果都被送入每個決策樹。每個觀察結果最常用作最終輸出。對所有決策樹進行新的觀察，并對每個分類模型進行多數投票。

隨機森林Radom Forest的模型R代碼，如圖6-27所示：

運行結果，如圖6-27所示

圖6-27

同樣，在final數據表中存儲randomforest模型的AUC：

final$randomforest <- completedData$test_pred

6.10.4? 模型評估

利用基于plotROC包和ggplot2軟件包可以將模型得到的AUC都展示在同一張ROC曲線圖中，這樣比較直觀比較。如果提示出現“Error in library(plotROC) : 不存在叫‘plotROC’這個名字的程輯包”，可以通過 “install.packages("plotROC", dependencies = TRUE)”命令安裝plotROC包。

R代碼如下所示：

library(plotROC)
library(ggplot2)
longtest <- melt_roc(final, "label", c("Lasso", "randomforest"))
cbPalette <- c("#FF0000", "#FFCC00")
cbPalette <- c("#FF0000", "#FFCC00")
ggplot(longtest, aes(d = D, m = M, color = name)) +scale_colour_manual(values=cbPalette) + geom_roc(size = 1.75, n.cuts = 0, labels = FALSE) + style_roc() + theme(axis.text.x= element_text(size = 22), axis.title.x= element_text(size = 22), axis.text.y= element_text(size = 22), axis.title.y= element_text(size = 22), legend.text = element_text(size = 22))

運行結果，如圖6-28所示。可以看出隨機森林的AUC最大。