拓端數據|R語言隨機森林RandomForest、邏輯回歸Logisitc預測心臟病數據和可視化分析

原文鏈接：http://tecdat.cn/?p=22596 

原文出處：拓端數據部落公眾號

 

研究大綱

• 介紹數據集和研究的目標
• 探索數據集
  • 可視化
  • 使用Chi-Square獨立檢驗、Cramer's V檢驗和GoodmanKruskal tau值對數據集進行探索
• 預測模型，Logisitic回歸和RandomForest
  • 兩個邏輯回歸的實例
  • 使用5折交叉驗證對模型實例進行評估
  • 變量選擇改進
    • step()
    • bestglm()
  • 隨機森林模型
  • 用RandomForest和Logisitc回歸進行預測
• 使用可視化進行最終的模型探索
• 結論和下一步措施

 

1.簡介

本報告是對心臟研究的機器學習/數據科學調查分析。更具體地說，我們的目標是在心臟研究的數據集上建立一些預測模型，並建立探索性和建模方法。但什么是心臟研究？

我們閱讀了關於FHS的資料：

心臟研究是對社區自由生活的人群中心血管疾病病因的長期前瞻性研究。心臟研究是流行病學的一個里程碑式的研究，因為它是第一個關於心血管疾病的前瞻性研究，並確定了風險因素的概念。

該數據集是FHS數據集的一個相當小的子集，有4240個觀測值和16個變量。這些變量如下：

1. 觀測值的性別。該變量在數據集中是一個名為 "男性 "的二值。
2. 年齡：體檢時的年齡，單位為歲。
3. 教育 : 參與者教育程度的分類變量，有不同的級別。一些高中（1），高中/GED（2），一些大學/職業學校（3），大學（4)
4. 目前吸煙者。
5. 每天抽的煙的數量
6. 檢查時使用抗高血壓葯物的情況
7. 流行性中風。流行性中風（0 = 無病）。
8. 流行性高血壓（prevalentHyp）。流行性高血壓。如果接受治療，受試者被定義為高血壓
9. 糖尿病。根據第一次檢查的標准治療的糖尿病患者
10. 總膽固醇(mg/dL)
11. 收縮壓(mmHg)
12.  舒張壓(mmHg)
13. BMI: 身體質量指數，體重（公斤）/身高（米）^2
14. 心率（次/分鍾）
15. 葡萄糖。血糖水平(mg/dL)

最后是因變量：冠心病（CHD）的10年風險。

 

這4240條記錄中有3658條是完整的病例，其余的有一些缺失值。

---

2.了解數據的意義

在每一步之前，要加載所需的庫。

1.  
   require(knitr)
2.  
   require(dplyr)
3.  
   require(ggplot2)
4.  
   require(readr)
5.  
   require(gridExtra) #呈現多幅圖

 

然后，加載心臟研究的數據集。

2.1 變量和數據集結構的檢查

我們對數據集進行一次檢查。

dim(dataset)

kable(head(dataset))

1.  
    
2.  
   str(dataset)

1.  
   ##檢查變量的摘要
2.  
   summary(dataset)

2.2 數據集的單變量圖

生成一個數據集的所有單變量圖。

1.  
   # 需要刪除字符、時間和日期等變量
2.  
    
3.  
   geom_bar(data = dataset,
4.  
   theme_linedraw()+
5.  
    
6.  
   #colnames(dataset)
7.  
   marrangeGrob(grobs=all_plots, nrow=2, ncol=2)

這是為了獲得對變量，對整個問題和數據集的理解，將通過多變量或至少雙變量的可視化來實現。
 

2.3 數據集的雙變量圖：因變量和預測因素之間的關系

 

現在我們可以進行一些雙變量的可視化，特別是為了看到因變量（TenYearCHD）和預測因素之間的關系。由於圖的數量太多，不是所有的一對變量都能被調查到！我們可以在后面的步驟中繼續調查。我們可以稍后再回到這一步，深入了解。

下面的代碼可以生成因變量的所有雙變量圖。由於因變量是一個二元變量，所以當預測變量是定量的時候，我們會有boxplots，或者當預測變量是定性的時候，我們會有分段的bar圖。

1.  
    
2.  
   for (var in colnames(dataset) ){
3.  
   if (class(dataset[,var]) %in% c("factor","logical") ) {
4.  
   ggplot(data = dataset) +
5.  
   geom_bar( aes_string(x = var,
6.  
    
7.  
    
8.  
   } else if (class(dataset[,var]) %in% c("numeric","double","integer") ) {
9.  
   ggplot(data = dataset) +
10.  
    geom_boxplot()
11.  
     

根據我們掌握的情況，男性與TenYearCHD直接相關，因此男性這個變量似乎是一個相對較好的預測因素。同樣，年齡似乎也是一個很好的預測因素，因為TenYearCHD == TRUE的病人有較高的年齡中位數，其分布幾乎相似。相反，不同類別的教育和因變量之間似乎沒有關系。目前的吸煙者變量與因變量有輕微的關系，因為目前的吸煙者患TenYearCHD的風險略高。

2.4 使用Goodman&Kruskal tau檢驗定性變量之間的關系

然而，除了這些本質上是定性方法的圖表外，人們可能希望對這種關聯有一個數字值。為了有這樣的數字測量，我想使用Goodman&Kruskal的tau測量，這是兩個無序因子，即兩個分類/名義變量之間的關聯測量。在我們這個數據集中的因子變量中，只有教育是序數變量，即它的類別有意義。這種測量方法比Cramer's V或chi-square測量方法更具信息量。

1.  
   GKtauData(cat_variables)
2.  
   plot(dataset)

可以看出，關於因變量的變異性，預測因素的解釋力非常小。換句話說，根據Goodman和Kruskal's tau度量，我們的預測因素和因變量之間幾乎沒有關聯。這可以從TenYearCHD一欄的數值中看出。

假設我的G&Ktau檢驗正確的話，這對模型來說並不是一個好消息。

為了檢驗這些發現，我們可以用Chi-square檢驗來檢驗分類變量與因變量的關聯的顯著性，然后用Phi相關系數來評估可能的關聯的強度。Phi用於2x2等值表。對於更大的表格，即有更多層次的變量，可以利用Cramer's V。

chisq.test(table(dataset_cat$p.value ))

1.  
    
2.  
   phi(matrix(table(dataset_cat_variables[,7],

奇怪的是，當Chi-square的P值如此之低時，可能的關聯的顯著性為零。這兩個測試（Chi-square和Phi相關）在大量的觀察中基本上得出相同的結果，因為一個是基於正態分布的，另一個是基於t分布的。
 

2.5 多重共線性的雙變量分析

該模型的真正問題在於共線性現象。共線性關系發生在兩個預測因子高度相關的情況下。我們需要檢查這種特性，然后繼續建立對數回歸模型。

根據Goodman和Kruskal's tau圖，我們不應該擔心共線性。但是，有序變量的教育變量呢？Cramer's V檢驗顯示，其強度不大。

1.  
    
2.  
   # 教育與其他分類變量的Chi square獨立性測試
3.  
   chisq.test(table(education,variables[,x]))$p.value )

1.  
   #將教育變量重新定位到數據集的第一個變量上
2.  
   assocstats(x = table(dataset_cat_variables[,1], dataset_$cramer ) )

沒有一個變量顯示與教育有很強的關聯。Cramer's V的最高值是0.145，這在教育和性別之間是相當弱的。
 

但是諸如currentSmoker和cigsPerDay這樣的變量呢？很明顯，其中一個是可以預測的。有一個數字變量和一個分類變量，我們可以把數字變量分成幾個類別，然后使用Goodman和Kruskal's tau。GroupNumeric()函數可以幫助將定量變量轉換成定性變量，然而，基於對數據的主觀理解，以及之前看到的cigsPerDay的多模態分布，在這里使用cut()函數很容易。
現在讓我們檢查一下GKtau的數值

 

1.  
   class_list <- lapply(X = 1:ncol(dataset_2), function(x) class(dataset_2[,x]))
2.  
   t <- sapply(X = names(class_list) , FUN = function(x) TRUE %in% ( class_list[x] %in% c("factor","logical")) )
3.  
    
4.  
   dataset_cat_variables_2 <- subset(x = dataset_2, select = t )
5.  
    
6.  
   plot(dataset_2)

從矩陣圖上的tau值及其背景形狀，我們可以看到cigsPerDay可以完全解釋currentSmoker的變異性。這並不奇怪，因為如果我們知道一個人每天抽多少支煙就可以斷言我們知道一個人是否是吸煙者!
 

第二個關聯是cigsPerDay與男性的關系，但它並不強烈。因此，前者可以解釋后者的較小的變化性。
 

在下一個數據集中，我把所有定量變量轉換成定性/分類變量。現在我們可以有一個全面的矩陣，盡管由於轉換，一些信息會丟失。

1.  
    
2.  
   dataset_3$totChol <- GroupNumeric(x = dataset$totChol , n = 5 )

我們可以看到，sysBP和diaBP可以預測prevalentHyp，但不是很強。(0.5左右)。因此我們可以在模型中保留prevalentHyp。第二點是關於GK tau的輸出。

 

3.預測模型:Logistic回歸和RandomForest

現在是評估模型實例的時候了。在這里，我們把邏輯回歸稱為模型。

我們有兩個實例。

1. 一個包括所有原始變量的模型實例，特別是cigsPerday和currentSmoker變量

2. 一個包括所有原始變量的模型實例，除了currentSmoker，cigsPerday被轉換為一個因子變量
為了評估模型實例，我們可以使用數學調整訓練誤差率的方法，如AIC。另一種方法是使用驗證數據集，根據模型在這個數據集上的表現來評估模型。在后一種方法中，我選擇使用K-fold Cross-Validation(CV)技術，更具體地說是5-fold CV。在這里，還有其他一些技術，如留一法交叉驗證。

3.1 兩個Logistic回歸模型實例

1.  
    
2.  
    
3.  
   # 因為下一步的cv.glm()不能處理缺失值。
4.  
   # 我只保留模型中的完整案例。
5.  
   dataset_1 <- dataset[complete.cases(dataset),]
6.  
   glm(TenYearCHD ~ . , family = "binomial")
7.  
    

 

這個模型是基於原始數據集的。有缺失值的記錄被從數據集中省略，模型顯示變量男性、年齡、cigsPerDay、totChol、sysBP和葡萄糖是顯著的，而prevalentHyp在某種程度上是顯著的。

 glm(formula = TenYearCHD ~ . , family = "binomial")

 

在第二個模型實例中，重要變量與前一個模型實例相同。

一個非常重要的問題是，如何衡量這兩個模型實例的性能以及如何比較它們？有各種方法來衡量性能，但我在這里選擇了5折交叉驗證法。

為了進行交叉驗證和評估模型實例，我們需要一個成本函數。boot軟件包推薦的一個函數，是一個簡單的函數，它可以根據一個閾值返回錯誤分類的平均數。閾值默認設置為0.5，這意味着任何觀察到的超過50%的CHD機會都被標記為有持續疾病的TRUE病例。從醫學的角度來看，我把閾值降低到0.4，這樣即使是有40%機會得心臟病的病例，也會被標記為接受進一步的醫療關注。降低閾值，增加了假陽性率，從而增加了醫療費用，但減少了假陰性率，挽救了生命。我們可以使用敏感度或特異性作為成本函數。此外，也可以使用cvAUC軟件包將曲線下面積（AUC）與CV結合起來。

 

3.2 模型實例的交叉驗證評估

1.  
   model1_cv_delta <- cv.glm( model1, cost = cost, K = 5)$delta[1]
2.  
    
3.  
   kable(data.frame("model1" = model1_cv_delta ,

1.  
   kable(
2.  
   caption = "CV-Accuracy", digits = 4)

我們可以看到，兩個模型非常相似，然而，模型2顯示出輕微的優勢。准確率確實相當高。但是，讓我們看看我們是否可以通過刪除一些變量來改進model1。
 

3.3 通過變量選擇改進模型

我們看一下model1的總結。

summary(model1)

到現在為止，我們一直假設所有的變量都必須包含在模型中，除非是共線性的情況。現在，我們被允許通過刪除不重要的變量。這里有幾種方法，如前向選擇和后向選擇。

例如，后向選擇法是基於不顯著變量的P值。淘汰繼續進行，直到AIC顯示沒有進一步改善。還有stats::step()和bestglm::bestglm()函數來自動進行變量選擇過程。后者的軟件包及其主要函數有許多選擇信息標准的選項，如AIC、BIC、LOOCV和CV，而前者的逐步算法是基於AIC的。

 

1.  
   bestglm(Xy = dataset_1 , family = binomial , IC = "BIC")
2.  
   step(object = model1 )

現在讓我們來看看這兩個模型和它們的交叉驗證誤差。

bestglm_bic_model

基於BIC的bestglm::bestglm()將模型變量減少到5個：男性、年齡、cigsPerDay、sysBP和葡萄糖。所有的變量都是非常顯著的，正如預期的那樣。

summary(step_aic_model)

基於AIC的step()函數將模型變量減少到8個：男性、年齡、cigsPerDay,prevalentStroke、prevalentHyp、totChol、sysBP和glucose。值得注意的是，通過step()找到的最佳模型實例具有不顯著的變量。
 

1.  
   glm_cv_error <- cv.glm(
2.  
   glmfit = glm(formula
3.  
   family = binomial, data = dataset_1),
4.  
    
5.  
   step_cv_error <- cv.glm(glmfit = step_aic_model, cost = cost, K = 5)$delta[1]
6.  
    
7.  
   kable(bestglm_model_cv_error ,
8.  
   step_model_cv_error )
9.  
   )

交叉驗證誤分類誤差

 

1.  
   kable(data.frame("bestglm() bic model"
2.  
   "step() aic model"

交叉驗證-准確度

AIC方法和BIC方法都能產生相同的准確性。該選擇哪種方法呢？我寧願選擇AIC，因為該模型實例有更多的預測因素，因此更有洞察力。然而，選擇BIC模型實例也是合理的，因為它更簡明。與model1的准確度相比，我們通過變量選擇在准確度上有0.8475-0.842=0.00550.8475-0.842=0.0055的提高。然而，我們失去了關於其他預測因子和因變量關系的信息。
 

3.4 RandomForest模型

到目前為止，我只做了邏輯回歸模型。有更多的模型可以用來為當前的問題建模，而RandomForest是一個受歡迎的模型。讓我們試一試，並將結果與之前的模型進行比較。

1.  
   #---- 差是每個RF模型實例的CV輸出的錯誤分類率
2.  
    
3.  
   #---- 每個選定的樹的CV錯誤分類率的最終結果被繪制出來
4.  
    
5.  
   # 對於不同數量的樹，我們計算CV誤差。
6.  
   for (n in seq(50,1000,50))
7.  
   for (k in 1:5)
8.  
   rf_dataset_train <- dataset_1[fold_seq != k ,]
9.  
   rf_dataset_test <- dataset_1[fold_seq == k , ]
10.  
    rf_model <- randomForest( formula,
11.  
     
12.  
    kable(rf_df[sort(x = rf_df[,2])

1.  
    
2.  
    
3.  
   #----- 誤差基於RandomForest OOB，即RandomForest輸出的混淆矩陣
4.  
    
5.  
   for (n in seq(50,1000,50)) {
6.  
   counter <- counter + 1
7.  
   rf_model <- randomForest( formula ntree = n, x =
8.  
    
9.  
   }
10.  
    ggplot() +
11.  
    geom_point(data = rf_df , aes(x = ntree , y = accuracy)
12.  
     

在這里，我同時使用了CV和out-of-bag(OOB)來評估隨機森林性能。

我們可以看到，在50到1000棵樹的范圍內，RandomForest模型的最高精度可以通過設置CV方法的樹數等於400來獲得。圖中的紅線顯示了我們從邏輯回歸模型實例中得到的最佳CV精度。由於OOB的最高准確率高於CV的最高准確率，所以我選擇了CV的准確率，使其更加謹慎。ntree=400的CVaccuracy=0.8486CVaccuracy=0.8486，比最好的邏輯回歸模型差0.00020.0002! 然而，如果我們考慮OOB的准確性，那么RandomForest模型比最佳邏輯回歸模型好0.00120.0012。
 

在RF中，模型的准確性有所提高，但代價是失去了可解釋性。RF是一個黑箱，我們無法解釋預測因子和因變量之間的關系。

3.5 模型對個人數據如何預測？

這里為了完成這個報告，我想在一個新的數據集上增加一個預測部分。該數據集只有一條記錄，其中包括我自己的個人數據。換句話說，我已經創建了一個模型，我想知道它是否預測了我的CHD。

1.  
   > pred_data$年齡 <- 31
2.  
   > pred_data$教育 <- factor(4, levels = c(1,2,3,4))
3.  
   > pred_data$當前吸煙者 <- FALSE
4.  
   > pred_data$每日吸煙量 <- 0
5.  
   > pred_data$抗高血壓葯物 <- FALSE
6.  
   > pred_data$流行性中風 <- FALSE
7.  
   > pred_data$流行性高血壓 <- FALSE

邏輯回歸模型的預測輸出。

1.  
   glm_BIC_opt <- glm(data = dataset_1 , formula ,family = binomial )
2.  
   predict(glm_BIC_opt, newdata = pred_data)

隨機森林預測。

1.  
   rf_model <- randomForest( formula = . ,
2.  
   predict(rf_model, pred_data)

因此，現在看來，我沒有風險! 然而，正如我之前提到的，這些模型是為了教育和機器學習的實踐，而不是為了醫學預測！所以，我認為這些模型是有價值的。

4.最終模型探索

讓我們最后看一下這個模型

1.  
    
2.  
   dataset_3 <- dataset_2[complete.cases(dataset_2),]
3.  
   dataset_3_GK <-
4.  
   plot(dataset_3_GK)

1.  
   ggpplot(data = dataset , text.angle = 0,label.size =2 , order = 0 ) +
2.  
   scale_colour_manual(values = color)+
3.  
   scale_fill_manual(values = color)
4.  
    

結果大多符合預期。根據GKtau值，預測因子之間的關聯最小。這正是我們想要的，以避免共線性現象。
然而，平行坐標仍然顯示了一些有趣的點。例如，年齡組與 "十年健康發展 "結果之間的關聯很有意思。較低的年齡組在TenYearCHD==TRUE中的參與度很低，這意味着年齡與該疾病有正相關。另一方面，與男性相比，女性（男性==FALSE）在0支煙和[1,20]支煙組的貢獻更大。換句話說，男性傾向於抽更多的煙，而且是重度吸煙者。

桑吉圖可以產生更好的洞察力，因為我們可以沿着坐標軸觀察樣本。

5.結論

在這項研究中，為了建立預測模型，使用了包括4240個觀測值和16個變量的心臟研究的數據集。這些模型旨在預測十年后的冠心病（CHD）。
在對數據集進行探索后，利用邏輯回歸和隨機森林模型來建立模型。使用K-Fold Cross-Validation對模型進行了評估。

為了擴展這項研究，可以使用進一步的分類方法，如支持向量機（SVM）、梯度提升（GB）、神經網絡模型、K-近鄰算法，甚至決策樹。

---

最受歡迎的見解

1.從決策樹模型看員工為什么離職

2.R語言基於樹的方法：決策樹，隨機森林

3.python中使用scikit-learn和pandas決策樹

4.機器學習：在SAS中運行隨機森林數據分析報告

5.R語言用隨機森林和文本挖掘提高航空公司客戶滿意度

6.機器學習助推快時尚精准銷售時間序列

7.用機器學習識別不斷變化的股市狀況——隱馬爾可夫模型的應用

8.python機器學習：推薦系統實現（以矩陣分解來協同過濾）

9.python中用pytorch機器學習分類預測銀行客戶流失