乐趣区

关于算法:R语言随机森林RandomForest逻辑回归Logisitc预测心脏病数据和可视化分析

jiezi

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

钻研纲要

  • 介绍数据集和钻研的指标
  • 摸索数据集
  • 可视化
  • 应用 Chi-Square 独立测验、Cramer’s V 测验和 GoodmanKruskal tau 值对数据集进行摸索
  • 预测模型，Logisitic 回归和 RandomForest
  • step()
  • bestglm()
  • 两个逻辑回归的实例
  • 应用 5 折穿插验证对模型实例进行评估
  • 变量抉择改良
  • 随机森林模型
  • 用 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. 理解数据的意义

在每一步之前，要加载所需的库。

require(knitr)
require(dplyr)
require(ggplot2)
require(readr)
require(gridExtra) #出现多幅图

而后，加载心脏钻研的数据集。

2.1 变量和数据集构造的查看

咱们对数据集进行一次查看。

dim(dataset)

kable(head(dataset))

str(dataset)

## 查看变量的摘要
summary(dataset)

2.2 数据集的单变量图

生成一个数据集的所有单变量图。

# 须要删除字符、工夫和日期等变量

                                                     geom_bar(data = dataset,
                                                     theme_linedraw()+

#colnames(dataset)        
marrangeGrob(grobs=all_plots, nrow=2, ncol=2)

这是为了取得对变量，对整个问题和数据集的了解，将通过多变量或至多双变量的可视化来实现。

2.3 数据集的双变量图：因变量和预测因素之间的关系

当初咱们能够进行一些双变量的可视化，特地是为了看到因变量（TenYearCHD）和预测因素之间的关系。因为图的数量太多，不是所有的一对变量都能被考察到！咱们能够在前面的步骤中持续考察。咱们能够稍后再回到这一步，深刻理解。

上面的代码能够生成因变量的所有双变量图。因为因变量是一个二元变量，所以当预测变量是定量的时候，咱们会有 boxplots，或者当预测变量是定性的时候，咱们会有分段的 bar 图。

for (var in colnames(dataset) ){if (class(dataset\[,var\]) %in% c("factor","logical") ) {ggplot(data = dataset)  + 
                        geom\_bar(aes\_string(x = var,} else if (class(dataset\[,var\]) %in% c("numeric","double","integer") ) {ggplot(data = dataset) + 
                        geom_boxplot()

依据咱们把握的状况，男性与 TenYearCHD 间接相干，因而男性这个变量仿佛是一个绝对较好的预测因素。同样，年龄仿佛也是一个很好的预测因素，因为 TenYearCHD == TRUE 的病人有较高的年龄中位数，其散布简直类似。相同，不同类别的教育和因变量之间仿佛没有关系。目前的吸烟者变量与因变量有轻微的关系，因为目前的吸烟者患 TenYearCHD 的危险略高。

2.4 应用 Goodman&Kruskal tau 测验定性变量之间的关系

然而，除了这些实质上是定性办法的图表外，人们可能心愿对这种关联有一个数字值。为了有这样的数字测量，我想应用 Goodman&Kruskal 的 tau 测量，这是两个无序因子，即两个分类 / 名义变量之间的关联测量。在咱们这个数据集中的因子变量中，只有教育是_序数变量_，即它的类别有意义。这种测量方法比 Cramer’s V 或 chi-square 测量方法更具信息量。

GKtauData(cat_variables)
plot(dataset)

能够看出，对于因变量的变异性，预测因素的解释力十分小。换句话说，依据 Goodman 和 Kruskal’s tau 度量，咱们的预测因素和因变量之间简直没有关联。这能够从 TenYearCHD 一栏的数值中看出。

假如我的 G &Ktau 测验正确的话，这对模型来说并不是一个好消息。

为了测验这些发现，咱们能够用 Chi-square 测验来测验分类变量与因变量的关联的显著性，而后用 Phi 相关系数来评估可能的关联的强度。Phi 用于 2 ×2 等值表。对于更大的表格，即有更多层次的变量，能够利用 Cramer’s V。

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

        phi(matrix(table(dataset\_cat\_variables\[,7\],

奇怪的是，当 Chi-square 的 P 值如此之低时，可能的关联的显著性为零。这两个测试（Chi-square 和 Phi 相干）在大量的察看中基本上得出雷同的后果，因为一个是基于正态分布的，另一个是基于 t 散布的。

2.5 多重共线性的双变量剖析

该模型的真正问题在于共线性现象。共线性关系产生在两个预测因子高度相干的状况下。咱们须要查看这种个性，而后持续建设对数回归模型。

依据 Goodman 和 Kruskal’s tau 图，咱们不应该放心共线性。然而，有序变量的教育变量呢？Cramer’s V 测验显示，其强度不大。

#  教育与其余分类变量的 Chi square 独立性测试  
chisq.test(table(education,variables\[,x\]))$p.value )

# 将教育变量从新定位到数据集的第一个变量上 
        assocstats(x = table(dataset\_cat\_variables\[,1\], dataset_$cramer ) )

没有一个变量显示与教育有很强的关联。Cramer’s V 的最高值是 0.145，这在教育和性别之间是相当弱的。

然而诸如 currentSmoker 和 cigsPerDay 这样的变量呢？很显著，其中一个是能够预测的。有一个数字变量和一个分类变量，咱们能够把数字变量分成几个类别，而后应用 Goodman 和 Kruskal’s tau。GroupNumeric()函数能够帮忙将定量变量转换成定性变量，然而，基于对数据的主观了解，以及之前看到的 cigsPerDay 的多模态散布，在这里应用 cut()函数很容易。
当初让咱们检查一下 GKtau 的数值

class\_list <- lapply(X = 1:ncol(dataset\_2), function(x) class(dataset_2\[,x\]))
t <- sapply(X = names(class\_list) , FUN = function(x) TRUE %in% (class\_list\[x\] %in% c("factor","logical")) )

dataset\_cat\_variables\_2 <- subset(x = dataset\_2, select = t)

plot(dataset_2)

从矩阵图上的 tau 值及其背景形态，咱们能够看到 cigsPerDay 能够齐全解释 currentSmoker 的变异性。这并不奇怪，因为如果咱们晓得一个人每天抽多少支烟就能够断言咱们晓得一个人是否是吸烟者!

第二个关联是 cigsPerDay 与男性的关系，但它并不强烈。因而，前者能够解释后者的较小的变动性。

在下一个数据集中，我把所有定量变量转换成定性 / 分类变量。当初咱们能够有一个全面的矩阵，只管因为转换，一些信息会失落。

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 回归模型实例

# 因为下一步的 cv.glm()不能解决缺失值。# 我只保留模型中的残缺案例。dataset_1 <- dataset\[complete.cases(dataset),\] 
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 模型实例的穿插验证评估

model1\_cv\_delta <- cv.glm(model1, cost = cost, K = 5)$delta\[1\]

kable(data.frame("model1" = model1\_cv\_delta ,

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

咱们能够看到，两个模型十分类似，然而，模型 2 显示出轻微的劣势。准确率的确相当高。然而，让咱们看看咱们是否能够通过删除一些变量来改良 model1。

3.3 通过变量抉择改良模型

咱们看一下 model1 的总结。

summary(model1)

到当初为止，咱们始终假如所有的变量都必须蕴含在模型中，除非是共线性的状况。当初，咱们被容许通过删除不重要的变量。这里有几种办法，如前向抉择和后向抉择。

例如，后向抉择法是基于不显著变量的 P 值。淘汰持续进行，直到 AIC 显示没有进一步改善。还有 stats::step()和 bestglm::bestglm()函数来主动进行变量抉择过程。后者的软件包及其次要函数有许多抉择信息规范的选项，如 AIC、BIC、LOOCV 和 CV，而前者的逐渐算法是基于 AIC 的。

bestglm(Xy = dataset_1 , family = binomial , IC = "BIC")
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()找到的最佳模型实例具备不显著的变量。

glm\_cv\_error <- cv.glm(
                                 glmfit = glm(formula 
                                              family = binomial, data = dataset_1),

step\_cv\_error <- cv.glm(glmfit = step\_aic\_model, cost = cost, K = 5)$delta\[1\]

kable(bestglm\_model\_cv_error ,
                step\_model\_cv_error )
        )

穿插验证误分类误差

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

穿插验证 - 准确度

AIC 办法和 BIC 办法都能产生雷同的准确性。该抉择哪种办法呢？我宁愿抉择 AIC，因为该模型实例有更多的预测因素，因而更有洞察力。然而，抉择 BIC 模型实例也是正当的，因为它更扼要。与 model1 的准确度相比，咱们通过变量抉择在准确度上有 0.8475-0.842=0.00550.8475-0.842=0.0055 的进步。然而，咱们失去了对于其余预测因子和因变量关系的信息。

3.4 RandomForest 模型

到目前为止，我只做了逻辑回归模型。有更多的模型能够用来为以后的问题建模，而 RandomForest 是一个受欢迎的模型。让咱们试一试，并将后果与之前的模型进行比拟。

#---- 差是每个 RF 模型实例的 CV 输入的谬误分类率

#---- 每个选定的树的 CV 谬误分类率的最终后果被绘制进去 

# 对于不同数量的树，咱们计算 CV 误差。for (n in seq(50,1000,50))     
        for (k in 1:5) 
                     rf\_dataset\_train <- dataset\_1\[fold\_seq != k ,\]
                     rf\_dataset\_test <- dataset\_1\[fold\_seq == k , \]
                     rf_model <- randomForest( formula,

kable(rf\_df\[sort(x = rf\_df\[,2\])

#----- 误差基于 RandomForest OOB，即 RandomForest 输入的混同矩阵

for (n in seq(50,1000,50)) {
        counter <- counter + 1 
                rf_model <- randomForest(formula ntree = n, x =}
ggplot() + 
        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。

> pred_data$ 年龄 <- 31
> pred_data$ 教育 <- factor(4, levels = c(1,2,3,4)) 
> pred_data$ 以后吸烟者 <- FALSE 
> pred_data$ 每日吸烟量 <- 0 
> pred_data$ 抗高血压药物 <- FALSE
> pred_data$ 流行性中风 <- FALSE 
> pred_data$ 流行性高血压 <- FALSE

逻辑回归模型的预测输入。

glm\_BIC\_opt <- glm(data = dataset_1 , formula ,family = binomial)
predict(glm\_BIC\_opt, newdata = pred_data)

随机森林预测。

rf_model <- randomForest( formula = . ,
predict(rf\_model, pred\_data)

因而，当初看来，我没有危险! 然而，正如我之前提到的，这些模型是为了教育和机器学习的实际，而不是为了医学预测！所以，我认为这些模型是有价值的。

4. 最终模型摸索

让咱们最初看一下这个模型

dataset\_3 <- dataset\_2\[complete.cases(dataset_2),\]
dataset\_3\_GK <- 
plot(dataset\_3\_GK)

ggpplot(data = dataset , text.angle = 0,label.size =2 , order = 0) +
        scale\_colour\_manual(values = color)+
        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 机器学习分类预测银行客户散失

退出移动版