拓端資料|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. str(dataset)

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

2.2 資料集的單變量圖

生成一個資料集的所有單變量圖。

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

這是為了獲得對變量，對整個問題和資料集的了解，将通過多變量或至少雙變量的可視化來實作。

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

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

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

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

根據我們掌握的情況，男性與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. phi(matrix(table(dataset_cat_variables[,7],

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

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

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

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

1. # 教育與其他分類變量的Chi square獨立性測試
2. 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. dataset_cat_variables_2 <- subset(x = dataset_2, select = t )
4. plot(dataset_2)

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

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

在下一個資料集中，我把所有定量變量轉換成定性/分類變量。現在我們可以有一個全面的矩陣，盡管由于轉換，一些資訊會丢失。

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

這個模型是基于原始資料集的。有缺失值的記錄被從資料集中省略，模型顯示變量男性、年齡、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. 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. step_cv_error <- cv.glm(glmfit = step_aic_model, cost = cost, K = 5)$delta[1]
5. kable(bestglm_model_cv_error ,
6. step_model_cv_error )
7. )

交叉驗證誤分類誤差

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. #---- 每個標明的樹的CV錯誤分類率的最終結果被繪制出來
3. # 對于不同數量的樹，我們計算CV誤差。
4. for (n in seq(50,1000,50))
5. for (k in 1:5)
6. rf_dataset_train <- dataset_1[fold_seq != k ,]
7. rf_dataset_test <- dataset_1[fold_seq == k , ]
8. rf_model <- randomForest( formula,
9. kable(rf_df[sort(x = rf_df[,2])

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

在這裡，我同時使用了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. dataset_3 <- dataset_2[complete.cases(dataset_2),]
2. dataset_3_GK <-
3. 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)

結果大多符合預期。根據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機器學習分類預測銀行客戶流失

▍關注我們

【大資料部落】第三方資料服務提供商,提供全面的統計分析與資料挖掘咨詢服務,為客戶定制個性化的資料解決方案與行業報告等。

▍咨詢連結：http://y0.cn/teradat

▍聯系郵箱：[email protected]