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

阅读量：

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

本报告旨在对心脏病进行系统性分析

研究大纲

介绍数据集和研究的目标
探索数据集
可视化

通过卡方独立性检验、Cramer's V统计量以及GoodmanKruskal tau相关系数对数据集进行深入分析

预测模型，Logisitic回归和RandomForest
step()
bestglm()
两个逻辑回归的实例
使用5折交叉验证对模型实例进行评估
变量选择改进
随机森林模型
用RandomForest和Logisitc回归进行预测
使用可视化进行最终的模型探索
结论和下一步改进

相关视频

1.简介

我们阅读了关于FHS的资料：

心肌缺血性疾病的流行病学研究是对具有自由生活方式的社区人群的心血管疾病病因进行长期前瞻性调查的研究项目。
作为流行病学领域的重要研究成果之一,心脏研究因其开创性地开创了心血管疾病前瞻性研究的新纪元而闻名。
因为它开创性地开创了心血管疾病前瞻性研究的新纪元,心脏研究不仅在方法论上实现了革命性的突破,而且首次系统地确立了心血管疾病的风险因素概念。

该数据集是FHS数据集中的一部分，在规模上相对较小，在本研究中使用了4240个观测样本以及16个相关变量作为分析依据

观测值的性别。该变量在数据集中是一个名为 "男性 "的二值。
年龄：体检时的年龄，单位为岁。

教育 : 参与者的教育程度是一个分类变量,分为不同的等级. 其中,部分高中(1),高中/GED(2),部分大学/职业学校(3),以及大学(4)

目前吸烟者。
每天抽的烟的数量
检查时使用抗高血压药物的情况
流行性中风。流行性中风（0 = 无病）。
common hypertension (common Hyp). If enrolled in treatment, participants will be classified as hypertensive.
糖尿病。根据第一次检查的标准治疗的糖尿病患者
总胆固醇(mg/dL)
收缩压(mmHg)
舒张压(mmHg)
BMI: 身体质量指数，体重（公斤）/身高（米）^2
心率（次/分钟）
葡萄糖。血糖水平(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；而在預測變數為定性時則會呈現分組條形圖。

复制代码

 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 病人中观察到较高年龄值的群体其病灶分布特征基本一致；相比之下, 不同教育水平与其他因变量之间的关联程度较低; 通过分析当前调查数据发现, 吸烟者的比例与其对应的疾病风险呈现弱相关性

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

然而，在这些定性分析的基础上，研究者可能更倾向于获得一个具体的数值指标来量化这种关联关系。基于其在区分不同类别间的差异方面的能力，Goodman&Kruskal的tau系数被选作衡量两个无序分类变量之间相关程度的理想工具。在我们所分析的数据集中，仅有教育水平作为唯一的一个_顺序变量_存在这一特点。相较于Cramer's V和卡方检验等方法，在本研究中采用tau系数能够提供更为丰富的信息量。

复制代码

 GKtauData(cat_variables)

    
 plot(dataset)

可见地看出，在分析因变量Y的变化时（或：在考察因变量Y的变化情况时），预测因子X的作用范围极为有限（或：预测因子X对因变量Y的影响程度极微）。换言之，在应用Goodman and Kruskal's τ系数计算方法时（或：在采用Goodman and Kruskal's tau统计指标时），我们发现预测因子X与因变量Y之间的相关关系几乎不存在显著关联（或：两者之间的关联强度极其微弱）。通过查看TenYearCHD一栏的数据（或：通过参考TenYearCHD这一栏的数据），这一现象得到了进一步印证。

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

为了评估分类分组变量与响应变量之间的统计学意义, 该研究采用Pearson’s Chi-square test to evaluate the statistical significance between分组变量 and响应变量. 研究者随后应用Phi相关系数量化分组间差异的大小. 对于仅有两个类别的情况, Phi指数特别适合. 当研究涉及多于两个类别时, 则可考虑使用Cramer's V作为替代指标.

复制代码

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

复制代码

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

奇怪的是，在 Chi-square 检验中的 P 值极其低时却提示可能存在关联实际上并不存在的可能性。从统计学上讲，在处理大样本数据时这两种方法得出的结果基本一致。
这是因为它们分别采用了不同的统计理论基础：一种基于正态分布理论（Chi-square），另一种基于 t 分布理论（Phi 相关）。因此，在应用这些方法分析大规模观测数据时所得出的结果几乎相同。

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

该模型的本质问题是共线性现象的表现。共线性关系表现在两个预测因子之间存在高度相关性的关系。为了更好地理解和应用这一特性，建议首先进行相关的属性检查和评估。接着可以基于这些分析结果来构建相应的对数回归模型

基于Goodman与Kruskal提出的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 and Kruskal's tau指标能够有效衡量两者之间的关联程度此处分割的方式应基于对数据特性的深入理解以及以往经验在此情境下应用cut()函数是否可行

复制代码

 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)

从矩阵图中的相关系数及其分布特征可以看出cigsPerDay能够充分说明currentSmoker的变化程度。这其实并不令人意外,因为如果一个人每天抽的香烟数量已知,则可判断其是否属于吸烟者

除了与男性相关的关联外，cigsPerDay与男性的关系强度较弱。然而，前者能够部分解释后者微小的变化幅度。

对于下一个数据集而言，在将所有数值型变量转化为分类型变量之后，我们能够构建一个涵盖所有属性的矩阵。然而，在这种转化过程中可能会导致部分信息被丢失。

复制代码

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

观察到sysBP与diaBP具有对prevousHyp预测的能力，并且其预测能力相对较弱（约为0.5）。由此可见，在构建模型时应考虑prevousHyp这一因素。第二个重点涉及GK tau指标的表现。

3.预测模型:Logistic回归和RandomForest

现在是评估模型实例的时候了。在这里，我们把逻辑回归称为模型。

我们有两个实例。

该模型实例整合了全部基础变量

基于所有原始变量构建的一个模型实例中,其中currentSmoker这一变量未进行处理而被保留为原始数据,在这种情况下,我们将其转换为分类变量以符合分析需求

3.1 两个Logistic回归模型实例

复制代码

 # 因为下一步的cv.glm()不能处理缺失值。

    
 # 我只保留模型中的完整案例。
    
 dataset_1 <- dataset\[complete.cases(dataset),\] 
    
 glm(TenYearCHD ~ . , family = "binomial")

该模型建立在原始数据集的基础之上。未包含缺失值的数据将从该数据集中排除。研究表明，在该模型中显著影响结果的因素包括男性性别、年龄以及吸烟频率（cigsPerDay）、总胆固醇水平（totChol）、收缩压（sysBP）等参数；另外，在一定程度上也显示出高血压（prevalentHyp）的相关性

复制代码

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

在第二个模型实例中，重要变量与前一个模型实例相同。

一个关键问题是：如何评估这两个模型实例的性能并鉴别它们之间的差异？特别重要的是：有多种方法可用来进行性能评估；但在这里采用了五折交叉验证法

在进行交叉验证以评估模型实例时，成本函数是一个必需的工具。在统计学领域中，默认情况下我们使用R语言进行数据分析，并且对于缺失数据通常采用排除法处理。在回归分析中，默认情况下我们假设变量之间存在线性关系，并且这种假设可能并不适用于所有实际情况。为了提高模型预测精度，在构建回归模型时需要考虑多重共线性问题，并且我们需要对数据分布情况进行详细描述。此外，在选择合适的统计方法时还需要考虑样本量大小以及数据类型特征等因素的影响

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)

到目前为止，我们一直以来都假设所有的变量都必须包含在模型中，并且只有当遇到共线性问题时才会有例外情况存在。然而现在我们是可以删除那些无关紧要的变量来优化模型构建过程的。这里介绍几种常用的方法包括前向筛选法和后向筛选法等

比如，在使用统计模型时

复制代码

 bestglm(Xy = dataset_1 , family = binomial , IC = "BIC")

    
 step(object = model1 )

现在让我们来看看这两个模型和它们的交叉验证误差。

复制代码

    bestglm\_bic\_model

采用BIC准则下的bestglm包函数将模型自变量精简至5项：包括男性参与者、不同年龄的人群、每日吸烟支数以及高血压患者和血糖水平较高的个体。所有选定的自变量均为高度显著的结果

复制代码

    summary(step\_aic\_model)

按照AIC标准运用step()函数后，在线性回归模型中使模型中的自变量数量缩减为8个：其中包括男性这一人口统计特征以及涉及年龄因素的相关指标。值得注意的是，在应用逐步筛选过程中所得到的最佳回归模型包含了若干不具统计学意义的自变量

复制代码

 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）均能提供类似的准确性水平。那么该如何选择两者中的一个呢？相较于BIC方法而言, AIC表现出更强的优势。其优势在于能够捕捉到更多的影响因素,从而提供更为深入的数据分析能力。相比之下, BIC标准虽然较为保守，在计算过程中相对更为简洁。然而, 通过引入变量筛选机制,我们可以将分类精度提升约...

3.4 RandomForest模型

截至目前为止，我的工作仅限于构建逻辑回归模型。
存在多种替代方案来构建适用于当前问题的模型。
其中Random Forest作为一种备受推崇的方法。
让我们尝试一下，并对比之前的结果以获取新的见解。

复制代码

 #---- 差是每个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)来评估随机森林性能。

观察发现，在 $5\times 1\text{E}3$ 棵决策树范围内，默认交叉验证（CV）设置为4个树节点时可实现随机森林模型的最大精度水平。图中红色线条标记了基于逻辑回归实例所得的最佳交叉验证（CV）精度水平值。值得注意的是，在此情境下尽管交叉验证（CV）方法达到的高度准确率为 $OOB$ 方法提供了更高的准确度值（即更高），但出于更为谨慎的角度考量我还是选择了基于交叉验证（CV）计算出的结果数值作为最终确定的标准参数值设定为 $ntree=4\times 1\text{E}3$ 时交叉验证（cv_accuracy）值为 $O(8\times 1\text{E-3})$ 级别并且与最好逻辑回归模型相比略低了一个微小数值单位量级约为 $O(7\times 1\text{E-6})$ 差异程度而言相对较小的影响因素应当被合理考虑进去这样的结果表现仍然能够满足实际应用需求以及项目目标要求

在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值计算得出的结果表明预测因子间的相关性最低；这正是我们所期望的结果；以防止多重共线性现象的发生

桑吉图可以产生更好的洞察力，因为我们可以沿着坐标轴观察样本。

5.结论

在该研究项目中，在深入分析心脏相关数据集的基础上（包含4240个观测样本及16个变量），我们构建了一个预测体系以期在未来十年内准确预判心肌梗死的发生情况（CHD）。该研究采用了多组验证方法（K-Fold Cross-Validation）以确保结果的有效性与可靠性。

作为扩展这项研究的延续工作之一，在现有基础上可采用一系列分类技术来进一步探索数据特征，并结合现有的机器学习模型框架进行优化与改进。具体而言，支持向量机模型（SVM）是一个强大的二类分类器，在小样本数据集上表现突出；梯度提升方法（GB）则通过迭代优化弱学习器的方式显著提升了预测准确性；神经网络架构提供了高度灵活的数据建模能力；此外还涉及到了K最近邻分类法以及决策树模型等基础算法的应用。

点击文末**“阅读原文”**

获取全文完整代码 资料。

本文源自《R语言中的随机森林模型与逻辑回归分析及其对心脏病数据的预测与可视化研究》

本文中的数据分享到会员群，扫描下面二维码即可加群！

点击标题查阅往期内容

数据呈现|采用R语言中的逻辑回归分析方法进行葡萄酒分类研究；运用线性判别式分析LDA技术实现特征降维；基于广义可加模型（Generalized Additive Model）构建非线性关系预测模型；采用多变量自适应回归样条法（MARS）处理复杂数据模式；结合k近邻分类法（K-Nearest Neighbors）进行非参数化预测；使用二次判别分析法（Quadratic Discriminant Analysis）区分不同类别样本；基于决策树模型构建分类树结构用于模式识别；通过随机森林算法（Random Forest Algorithm）提升分类精度并减少过拟合风险；在葡萄酒交叉验证过程中绘制并评估ROC曲线表现

MATLAB随机森林优化贝叶斯预测分析汽车燃油经济性

R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病

R语言中贝叶斯网络（BN）、动态贝叶斯网络、线性模型分析错颌畸形数据

R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

Python贝叶斯回归分析住房负担能力数据集

利用R语言进行贝叶斯分位数回归的实现，并对least absolute shrinkage and selection operator以及其自适应版本的贝叶斯分位数回归分析

Python用PyMC3实现贝叶斯线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

R语言和STAN,JAGS：用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究