R数据科学（五）探索性数据分析

Posted 2023-05-01

参考技术A 明确概念：探索性数据分析（exploratory data analysis， EDA），一般过程为：
(1) 对数据提出问题。
(2) 对数据进行可视化、转换和建模，进而找出问题的答案。
(3) 使用上一个步骤的结果来精炼问题，并提出新问题。

确定变量是分类变量还是连续变量，要想检查分类变量的分布，可以使用条形图：

条形的高度表示每个 x 值中观测的数量，可以使用 dplyr::count() 手动计算出这些值：

要想检查连续变量的分布，可以使用直方图：

可以通过 dplyr::count() 和 ggplot2::cut_width() 函数的组合来手动计算结果. binwidth 参数来设定直方图中的间隔的宽度，该参数是用 x 轴变量的单位来度量的。

在同一张图上叠加多个直方图, 用geom_freqploy()代替geom_histogram()，用折线表示。

相似值聚集形成的簇表示数据中存在子组。

coord_cartesian() 函数中有一个用于放大 x 轴的 xlim() 参数。 ggplot2 中也有功能稍有区
别的 xlim() 和 ylim() 函数：它们会忽略溢出坐标轴范围的那些数据。
如果带有异常值和不带异常值的数据分别进行分析，结果差别较大的话要找出异常值的原因，如果差别不大，可以用NA代替。

5.3.4　练习
(1)研究 x、 y 和 z 变量在 diamonds 数据集中的分布。你能发现什么？思考一下，对于一条
钻石数据，如何确定表示长、宽和高的变量？

(2)研究 price 的分布，你能发现不寻常或令人惊奇的事情吗？（提示：仔细考虑一下
binwidth 参数，并确定试验了足够多的取值。）

(3) 0.99 克拉的钻石有多少？ 1 克拉的钻石有多少？造成这种区别的原因是什么？

(4)比较并对比 coord_cartesina() 和 xlim()/ylim() 在放大直方图时的功能。如果不设置
binwidth 参数，会发生什么情况？如果将直方图放大到只显示一半的条形，那么又会发
生什么情况？

数据中有异常值，可以将异常值去掉：

一般不建议去掉，建议使用缺失值来代替异常值。

ifelse函数参数1放入逻辑判断，如果为T，结果就是第二个参数的值，如果为F，就是第三个参数的值。
ggplot2会忽略缺失值：

练习
(1) 直方图如何处理缺失值？条形图如何处理缺失值？为什么会有这种区别？

(2) na.rm = TRUE 在 mean() 和 sum() 函数中的作用是什么？
移除缺失值再进行统计

按分类变量的分组显示连续变量分布的另一种方式是使用箱线图

练习
(1) 前面对比了已取消航班和未取消航班的出发时间，使用学习到的知识对这个对比的可视
化结果进行改善。

(2) 在钻石数据集中，哪个变量对于预测钻石的价格最重要？这个变量与切割质量的关系是
怎样的？为什么这两个变量的关系组合会导致质量更差的钻石价格更高呢？

(3) 安装 ggstance 包，并创建一个横向箱线图。这种方法与使用 coord_flip() 函数有何区别？

(4) 箱线图存在的问题是，在小数据集时代开发而成，对于现在的大数据集会显示出数量极
其庞大的异常值。解决这个问题的一种方法是使用字母价值图。安装 lvplot 包，并尝试
使用 geom_lv() 函数来显示价格基于切割质量的分布。你能发现什么问题？如何解释这
种图形？

(5) 比较并对比 geom_violin()、分面的 geom_histogram() 和着色的 geom_freqploy()。每种方法的优缺点是什么？

(6) 对于小数据集，如果要观察连续变量和分类变量间的关系，有时使用 geom_jitter() 函数是特别有用的。 ggbeeswarm 包提供了和 geom_jitter()相似的一些方法。列出这些方法
并简单描述每种方法的作用。

两个分类变量的关系肯定要先计数，可以用geom_count()函数
d3heatmap 或 heatmaply 包可以生成交互式图

练习
(1) 如何调整 count 数据，使其能更清楚地表示出切割质量在颜色间的分布，或者颜色在切
割质量间的分布？

(2) 使用 geom_tile() 函数结合 dplyr来探索平均航班延误数量是如何随着目的地和月份的
变化而变化的。为什么这张图难以阅读？如何改进？

(3) 为什么在以上示例中使用 aes(x = color, y = cut) 要比 aes(x = cut, y = color) 更好？

连续变量之间的关系一般用散点图来表示。geom_point()

对于大数据集，为了避免重合，可以用geom_bin2d() 和 geom_hex()函数将坐标平面分为二维分箱，并使用一种填充颜色表示落入
每个分箱的数据点。

另一种方式是对一个连续变量进行分箱，因此这个连续变量的作用就相当于分类变量。
cut_width(x, width) 函数将 x 变量分成宽度为 width 的分箱。参数 varwidth = TRUE 让箱线图的宽度与观测数量成正比。
cut_number() 函数近似地显示每个分箱中的数据点的数量

练习
(1) 除了使用箱线图对条件分布进行摘要统计，你还可以使用频率多边形图。使用 cut_
width() 函数或 cut_number() 函数时需要考虑什么问题？这对 carat 和 price 的二维分
布的可视化表示有什么影响？

(2) 按照 price 分类对 carat 的分布进行可视化表示。

(3) 比较特别大的钻石和比较小的钻石的价格分布。结果符合预期吗？还是出乎意料？

(4) 组合使用你学习到的两种技术，对 cut、 carat 和 price 的组合分布进行可视化表示。

(5) 二维图形可以显示一维图形中看不到的离群点。例如，以下图形中的有些点具有异常的
x 值和 y 值组合，这使得这些点成为了离群点，即使这些点的 x 值和 y 值在单独检验时
似乎是正常的。
ggplot(data = diamonds) +
geom_point(mapping = aes(x = x, y = y)) +
coord_cartesian(xlim = c(4, 11), ylim = c(4, 11))

数据中的模式提供了关系线索,用于探索两个变量的相关性。
模型是用于从数据中抽取模式的一种工具。
残差（预测值和实际值之间的差别）

阅读推荐：
生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！
B站链接： https://m.bilibili.com/space/338686099
YouTube链接： https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程师入门最佳指南： https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA

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

本报告是对心脏研究的机器学习/数据科学调查分析。更具体地说，我们的目标是在心脏研究的数据集上建立一些预测模型，并建立探索性和建模方法。但什么是心脏研究？

研究大纲

• 介绍数据集和研究的目标
• 探索数据集
• 可视化
• 使用Chi-Square独立检验、Cramers 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)

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

---

点击标题查阅往期内容

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

左右滑动查看更多

01

02

03

04

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测量，这是两个无序因子，即两个分类/名义变量之间的关联测量。在我们这个数据集中的因子变量中，只有教育是_序数变量_，即它的类别有意义。这种测量方法比Cramers V或chi-square测量方法更具信息量。

GKtauData(cat_variables)
plot(dataset)

可以看出，关于因变量的变异性，预测因素的解释力非常小。换句话说，根据Goodman和Kruskals tau度量，我们的预测因素和因变量之间几乎没有关联。这可以从TenYearCHD一栏的数值中看出。

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

为了检验这些发现，我们可以用Chi-square检验来检验分类变量与因变量的关联的显著性，然后用Phi相关系数来评估可能的关联的强度。Phi用于2x2等值表。对于更大的表格，即有更多层次的变量，可以利用Cramers 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和Kruskals tau图，我们不应该担心共线性。但是，有序变量的教育变量呢？Cramers V检验显示，其强度不大。

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

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

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

但是诸如currentSmoker和cigsPerDay这样的变量呢？很明显，其中一个是可以预测的。有一个数字变量和一个分类变量，我们可以把数字变量分成几个类别，然后使用Goodman和Kruskals 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)

左右滑动查看更多

01

02

03

04

结果大多符合预期。根据GKtau值，预测因子之间的关联最小。这正是我们想要的，以避免共线性现象。
然而，平行坐标仍然显示了一些有趣的点。例如，年龄组与 "十年健康发展 "结果之间的关联很有意思。较低的年龄组在TenYearCHD==TRUE中的参与度很低，这意味着年龄与该疾病有正相关。另一方面，与男性相比，女性（男性==FALSE）在0支烟和[1,20]支烟组的贡献更大。换句话说，男性倾向于抽更多的烟，而且是重度吸烟者。

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

5.结论

在这项研究中，为了建立预测模型，使用了包括4240个观测值和16个变量的心脏研究的数据集。这些模型旨在预测十年后的冠心病（CHD）。
在对数据集进行探索后，利用逻辑回归和随机森林模型来建立模型。使用K-Fold Cross-Validation对模型进行了评估。

为了扩展这项研究，可以使用进一步的分类方法，如支持向量机（SVM）、梯度提升（GB）、神经网络模型、K-近邻算法，甚至决策树。

---

---