本章将向大家展示如何使用可视化和转换,以系统的方式探索数据,统计学家将这项任务称为探索性数据分析(Exploratory
Data Analysis),简称 EDA。
EDA
是一个迭代循环:1)对数据提出相关问题;2)通过可视化、转换和建模数据来搜索答案;3)使用您学到的知识来完善您的问题和/或产生新问题。依次循环达到对数据充分理解的目的。
EDA不是一个具有严格规则的正式流程。最重要的是,EDA是一种心态。在EDA的初始阶段,你应该随意测试您的每个想法。其中一些想法会成功,而另一些则会是死胡同。随着探索的继续,你将专注于一些特别富有成效的点,并最终把它们写下来并与其他人交流。
EDA期间的目标: 加深对数据的理解
最简单的方法是使用问题作为指导调查的工具。当您提出问题时,该问题会将您的注意力集中在数据集的特定部分上,并帮助您决定要进行哪些图形、模型或转换。
EDA从根本上来说是一个创造性的过程。与大多数创意过程一样,提出高质量问题的关键是产生大量问题。
在分析开始时很难提出有启发性的问题,因为您不知道数据集中包含哪些见解。
另一方面,你提出的每个新问题都会让你接触到数据的新方面,并增加你新发现的机会。
如果你根据新发现提出一个新问题来跟进每个问题,则可以快速深入了解数据中最有趣的部分,并提出一组发人深省的问题。
对于应该提出哪些问题来指导你的研究,没有任何规定。然而,有两种类型的问题一般是通用的,可大致表述为:
1)我的变量中发生什么类型的变化?
2)我的变量之间发生什么类型的协变?
下面我们将探讨这两个问题,并以此延申出许多小问题。我将解释什么是变异和协变,并向你展示几种回答每个问题的方法。
1)变量是指研究对象的某个特征或属性,在不同个体或观察单位之间可以有差异或变化的量。它是我们感兴趣的、需要收集和测量的事物或现象的表示。
2)变异是变量值在不同测量之间发生的变化。每个变量都有其自身的变异模式,可以揭示有趣的信息。理解该模式的最佳方法是通过可视化变量值的分布来实现。
此处仅代表我的环境中需要解除的冲突,可根据自己的情况选用。
library(tidyverse, warn.conflicts = FALSE)
data("diamonds")
1)数据中包含哪些变量,以及变量的类型?
2)每个变量的值大概是个什么样子的?
3)数值型变量的最大值、最小值、均值、中值、四分位值是多少?变异范围是指?
4)定性变量的各描述性状的分布情况?(定性变量在R中的数据框中标注为factor)
head(diamonds,5)
## # A tibble: 5 × 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
summary(diamonds)
## carat cut color clarity depth
## Min. :0.2000 Fair : 1610 D: 6775 SI1 :13065 Min. :43.00
## 1st Qu.:0.4000 Good : 4906 E: 9797 VS2 :12258 1st Qu.:61.00
## Median :0.7000 Very Good:12082 F: 9542 SI2 : 9194 Median :61.80
## Mean :0.7979 Premium :13791 G:11292 VS1 : 8171 Mean :61.75
## 3rd Qu.:1.0400 Ideal :21551 H: 8304 VVS2 : 5066 3rd Qu.:62.50
## Max. :5.0100 I: 5422 VVS1 : 3655 Max. :79.00
## J: 2808 (Other): 2531
## table price x y
## Min. :43.00 Min. : 326 Min. : 0.000 Min. : 0.000
## 1st Qu.:56.00 1st Qu.: 950 1st Qu.: 4.710 1st Qu.: 4.720
## Median :57.00 Median : 2401 Median : 5.700 Median : 5.710
## Mean :57.46 Mean : 3933 Mean : 5.731 Mean : 5.735
## 3rd Qu.:59.00 3rd Qu.: 5324 3rd Qu.: 6.540 3rd Qu.: 6.540
## Max. :95.00 Max. :18823 Max. :10.740 Max. :58.900
##
## z
## Min. : 0.000
## 1st Qu.: 2.910
## Median : 3.530
## Mean : 3.539
## 3rd Qu.: 4.040
## Max. :31.800
##
或者查看某一变量,如定性变量(Qualitative Variable) cut.
summary(diamonds$cut)
## Fair Good Very Good Premium Ideal
## 1610 4906 12082 13791 21551
定性变量(Qualitative Variables)是用于描述特征、属性或类别的变量,不能用数值进行度量,可分为:<br>
1)名义变量(Nominal Variables):名义变量是没有顺序或等级关系的变量,只能用于标识不同的类别或群组。例如,性别(男、女)、民族(汉族、满族、回族等)、颜色(红、绿、蓝等)都是名义变量。
2)有序变量(Ordinal Variables):有序变量是具有顺序或等级关系的变量,可以用于表示不同类别之间的相对大小或顺序。例如,教育水平(小学、初中、高中、大学等)、满意度评分(非常满意、满意、一般、不满意等)都属于有序变量。
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut))
#这段代码使用了tidyverse包中的库,主要目的是创建一个基于diamonds数据集的图形。图形使用ggplot函数创建,并通过geom_bar函数添加柱状图层。在aes函数内部,使用x参数指定了x轴的变量为cut,即钻石的切割质量。最终,你将得到一个显示不同钻石切割质量频率的柱状图。
定量变量(Quantitative Variable):也称为连续变量或数值变量,表示可量化或测量的数值属性的变量。进一步分为两类:<br>
1)离散变量(Discrete Variable):取有限个或可数个数值的变量,只能为整数。<br>
2)连续变量(Continuous Variable):可以取任意数值的变量,可以有小数。<br>
例如:该diamonds数据集中,定量变量carat(克拉数)的分布情况,也用直方图表示:
ggplot(data = diamonds) +
geom_histogram(mapping = aes(x = carat), binwidth = 0.5)
#在aes函数内部,使用x参数指定了x轴的变量为carat,即钻石的克拉数。binwidth参数设置了直方图的箱宽为0.5,即将数据分成了大小为0.5的小组。直方图将x轴分为等距的区间,并使用每个区间的柱高来显示落在该区间内的观测数量。
在使用直方图时,可应尝试不同的 binwidth 来探索不同的模式,找到适合的箱宽,例如:
ggplot(data = diamonds) +
geom_histogram(mapping = aes(x = carat), binwidth = 0.1)
例如:在指定克拉数的数值区间内的各有多少个钻石?
diamonds %>%
count(cut_width(carat, width = 0.5))# width=0.5的数值区间根据自己的兴趣和好奇心置换
## # A tibble: 11 × 2
## `cut_width(carat, width = 0.5)` n
## <fct> <int>
## 1 [-0.25,0.25] 785
## 2 (0.25,0.75] 29498
## 3 (0.75,1.25] 15977
## 4 (1.25,1.75] 5313
## 5 (1.75,2.25] 2002
## 6 (2.25,2.75] 322
## 7 (2.75,3.25] 32
## 8 (3.25,3.75] 5
## 9 (3.75,4.25] 4
## 10 (4.25,4.75] 1
## 11 (4.75,5.25] 1
#在count函数的参数中,cut_width函数将连续变量carat按照指定的宽度(0.5)划分成不同的分组。
diamonds %>%
count(cut_width(carat, center = 0.25, width = 0.5))
## # A tibble: 10 × 2
## `cut_width(carat, center = 0.25, width = 0.5)` n
## <fct> <int>
## 1 [0,0.5] 18932
## 2 (0.5,1] 17506
## 3 (1,1.5] 12060
## 4 (1.5,2] 3553
## 5 (2,2.5] 1763
## 6 (2.5,3] 94
## 7 (3,3.5] 23
## 8 (3.5,4] 4
## 9 (4,4.5] 4
## 10 (5,5.5] 1
%>%:允许将左侧表达式的结果传递给右侧表达式作为第一个参数。这种符号常用于dplyr包和magrittr包中,用于链式操作数据集或函数。
%in%:用于判断一个元素是否出现在向量或列表中;
%*%:用于矩阵的乘法操作;
%/%:用于整数除法,返回商的整数部分。
通过不断尝试替换width的区间值,可以找到你想要的答案。
例如:我们项了解克拉数小于等于3的钻石分布情况以及占比
smaller <- diamonds %>%
filter(carat <= 3) # 符号<= 代表小于等于
ggplot(data = smaller, mapping = aes(x = carat)) +
geom_histogram(binwidth = 0.1)
P <- (nrow(smaller) / nrow(diamonds)) * 100
P
## [1] 99.94067
#在filter函数中,通过carat < 3的条件筛选出diamonds数据集中克拉数小于3的钻石,得到一个新的数据集smaller。
#接下来,使用ggplot函数创建一个图形,并设置数据为smaller。
#然后,使用geom_histogram函数添加直方图层,binwidth参数设置直方图的箱宽为0.1
例如:在小于3克拉的钻石中,根据cut变量对钻石的质量进行分组,每组中的钻石克拉数的分布情况是怎样的?
用折线图比直方图的可视化效果更好,因此建议使用 geom_freqpoly() 而不是
geom_histogram()。
ggplot(data = smaller, mapping = aes(x = carat, colour = cut)) +
geom_freqpoly(binwidth = 0.1)
# 设置数据为smaller。在aes函数中,指定x轴的变量为carat,即钻石的克拉数。另外,还指定了颜色映射变量为cut,即钻石的切工。
# 注意:函数中的colour是指aes()函数中的参数设置,而不是diamonds数据集中的color,不要混淆。
# 使用geom_freqpoly函数添加频率多边形层。频率多边形是一种用直线段连接各个频数点的图形,它反映了变量的分布情况。
我们上面的一些实例问题,目的是引发你思考,引导你提出自己感兴趣的其他的问题。<br>
通过上面的例子,现在你已经能够可视化变量的分布图了,那么在你的图表中应该寻找哪些信息?以及你应该提出什么类型的后续问题?<br>
提出良好的后续问题的关键将是依靠你的好奇心(你想进一步了解什么?)以及你的怀疑心态(这可能会产生误导吗?)。<br>
我下面列出了的一些有用的信息类型,以及每种信息类型的一些后续问题。<br>
Typical values是指描述数据集中观测值的一些典型特征或常见值。这些值可用于了解数据的中心趋势、离散程度和分布情况。
在条形图和直方图中,高柱表示变量的常见值,较短的柱表示不太常见的值。没有柱的地方表示数据中未出现的值。要将这些信息转化为有用的问题,请寻找任何意外情况:<br>
1)哪些值最常见?为什么?<br>
2)哪些值很少见?为什么?这是否符合你的预期?<br>
3)你能看到任何异常模式吗?可能的解释是什么?<br>
ggplot(data = smaller, mapping = aes(x = carat)) +
geom_histogram(binwidth = 0.01)
1)为什么整克拉的钻石和常见分数克拉的钻石更多?整数克拉如1、2、3、4,常见分数如0.5、1.5、2.5,不常见分数如0.48、0.88、1.47等。
这个问题涉及到钻石的市场需求和生产的一些因素,以下是一个可能的解释:
完整的克拉数和常见分数的钻石更受欢迎和更常见的原因是市场需求的影响。大多数消费者在购买钻石时更倾向于选择常见的克拉数或分数,因为它们比较容易理解和比较。例如,人们更容易理解0.5克拉的钻石相对于0.48克拉的钻石,即使两者几乎没有实质上的区别。这种心理因素导致了对完整克拉数和常见分数的钻石的更高需求,这也使得它们更广泛地出现在市场上。
2)为什么每个峰值右侧的钻石稍多一些,而左侧的稍少一些?
这个问题涉及到统计学中的偏斜(Skewness)概念。如果每个峰值的右侧比左侧有更多的钻石,那意味着分布是向右倾斜的,并具有正偏斜。
正偏斜表示数据分布的尾部相对较长,也就是说,在峰值的右侧,有更多的极端观测值。对于钻石的情况来说,即表示有更多的钻石略微大于每个峰值的克拉数。
这种正偏斜现象主要受到市场需求和价格的影响。通常情况下,较大克拉数的钻石更为稀有和昂贵。由于供应受限,钻石价格会随着克拉数的增加而上升。因此,相对而言,在峰值右侧的略微大于克拉数的钻石会比峰值左侧的略微小于克拉数的钻石更多。
此外,人们在购买钻石时也有一种心理倾向,即倾向于选择略微大于常见克拉数或分数的钻石。这种心理因素进一步增加了右侧略微大于峰值的钻石需求。
3)为什么没有大于 3 克拉的钻石?
相似值的聚集表明在您的数据中存在子群体。为了理解这些子群体,可以提出以下问题:
1)每个聚类中的观察结果相互之间有何相似之处?
2)不同聚类中的观察结果有何不同之处?
3)如何解释或描述这些聚类?
4)为什么聚类的出现可能具有误导性?
Unusual
values是不寻常的观测值,即不符合整体规律的数据点,包括离群点、高杠杆点和强影响点。
离群点是指那些模型预测效果不佳的观测点。
# 创建一个包含一些离群点的示例数据集
data <- c(1, 2, 3, 4, 5, 7, 10, 15, 20, 22, 25, 60, 70)
# 绘制箱线图以检测离群点
boxplot(data)
# 计算上下边界
Q1 <- quantile(data, 0.25)
Q3 <- quantile(data, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
# 打印出超过上下边界的离群点
outliers <- data[data < lower_bound | data > upper_bound]
print(outliers)
## [1] 60 70
在回归分析中,对模型拟合程度产生较大影响的观测点。具体而言,高杠杆值点在自变量取值上远离其他观测点,因此在对回归模型进行参数估计和预测时具有较大的权重。
e.g.
set.seed(123)
n <- 50 # 观测值数量
x <- seq(1, n)
y <- 2 * x + rnorm(n) * 10 # 添加随机干扰项
# 将第50个观测值y设为一个极端值以创建高杠杆值点
y[47:50] <- c(165, 175, 160, 200)
# 执行线性回归
model <- lm(y ~ x)
# 计算每个观测点的杠杆值
leverage <- hatvalues(model)
# 定义杠杆值阈值为2 * (p + 1) / n,其中p是自变量的数量
leverage_threshold <- 2 * (length(coef(model)) + 1) / length(residuals(model))
# 打印出高杠杆值点的索引
high_leverage_points <- which(leverage > leverage_threshold)
print(high_leverage_points)
## named integer(0)
# 绘制原始数据和回归线
plot(x, y, main = "高杠杆值点示例", xlab = "x", ylab = "y")
abline(model, col = "red")
指在统计模型中对参数估计和模型拟合具有显著影响的观测点。若移除该点时,模型会发生巨大的变化。
# 创建包含强影响点的示例数据集
set.seed(123)
n <- 10 # 观测值数量
x <- seq(1, n)
y <- 2 * x + rnorm(n) * 10 # 添加随机干扰项
# 将第50个观测值y设为一个极端值以创建强影响点
y[10] <- 200
# 执行线性回归
model <- lm(y ~ x)
# 计算每个观测点的影响度指标(Cook's距离)
influence <- influence.measures(model)
# 定义影响度阈值为4 / (n - p - 1),其中p是自变量的数量
influence_threshold <- 4 / (length(y) - length(coef(model)) - 1)
# 打印出强影响点的索引
high_influence_points <- which(influence$cook.d > influence_threshold)
print(high_influence_points)
## integer(0)
# 绘制原始数据和回归线
plot(x, y, main = "强影响点例", xlab = "x", ylab = "y")
abline(model, col = "red")
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5)
有时异常值是数据录入错误,<br>
而其他情况下则可能暗示着重要的新科学发现。<br>
当您拥有大量数据时,在直方图中很难看到异常值。例如,看看钻石数据集中的 y变量的分布。唯一能看出异常值的证据是 x 轴上异常宽广的极限范围。<br>
常见值的直方柱中有如此多的观测值,以至于罕见值的直方柱非常短,你看不到它们(尽管也许如果你专注地盯着0,可能会发现一些)。<br>
为了方便查看异常值,我们可以使用coord_cartesian()将y轴缩小到较小的值范围内:
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5) +
coord_cartesian(ylim = c(0, 50))#定义y轴的显示范围
#(coord_cartesian() 还有一个 xlim() 参数,用于缩小到 x 轴的范围。ggplot2 还有 xlim() 和 ylim() 函数,它们的工作方式略有不同:它们会丢弃超出限制范围之外的数据。)
这样我们就可以看到有三个异常值:0、约30和约60。我们可以使用dplyr将它们提取出来:
unusual <- diamonds %>%
filter(y < 3 | y > 20) %>%
select(carat, price, x, y, z)
#通过filter函数对diamonds数据集进行筛选。使用逻辑运算符|(或)筛选出y轴的值小于3或大于20的数据,得到一个新的数据集unusual。<br>
#然后,使用select函数选择了unusual数据集中的carat、price、x、y和z这四个变量,即钻石的克拉、价格、长、宽和高。
显示稀有值(也叫离散数值)
unusual
## # A tibble: 9 × 5
## carat price x y z
## <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 5139 0 0 0
## 2 1.14 6381 0 0 0
## 3 2 12210 8.09 58.9 8.06
## 4 1.56 12800 0 0 0
## 5 1.2 15686 0 0 0
## 6 2.25 18034 0 0 0
## 7 0.51 2075 5.15 31.8 5.12
## 8 0.71 2130 0 0 0
## 9 0.71 2130 0 0 0
从上面的数据unusual中,请分析是什么原因导致的这些异常值?
Notice:
重复检测数据中是否包括异常值,是一个良好的数据分析习惯。
1)如果异常值对结果的影响很小,而且您无法弄清楚它们存在的原因,将其替换为缺失值并继续分析是合理的。
2)如果异常值对结果产生重大影响,您不能无故将其删除。您需要弄清楚是什么导致了这些异常值(例如数据输入错误或则其他什么原因),并在报告中说明您已经移除了它们。
1)探索钻石中每个变量(x、y和z)的分布。你学到了什么?想想一颗钻石,你如何确定它的长度、宽度和深度。<br>
2) 探索价格的分布。你发现了什么异常或令人惊讶的情况吗?(提示:仔细考虑一下直方图的柱宽(binwidth),确保尝试一系列不同的值。)<br>
3) 有多少颗钻石重量为0.99克拉?有多少颗钻石重量为1克拉?你认为造成这种差异的原因是什么?<br>
4)比较和对比在直方图上使用 coord_cartesian() 和 xlim()、ylim()<br>
5) 进行缩放时的效果。如果不设置柱宽(binwidth),会发生什么?如果尝试只显示一半柱子的区域,会发生什么?<br>
数据中的缺失值可能由两类情况:<br>
1)因为本身在收集数据时,无法获取某变量的某个值,造成缺失,这时我们需要在相应的单元格中用NA代表缺失值。<br>
2)因为我们的数据中有异常值,需要移除异常值造成的缺失。下面我们针对异常值的处理进行介绍。
如果在数据集中遇到缺失值或异常值,而且想继续进行其他分析,主要有两种选择。<br>
diamonds2 <- diamonds %>%
filter(between(y, 3, 20))
#通过filter函数对diamonds数据集进行筛选。使用between函数,筛选出y轴的值在3和20之间(包括3和20)的数据。
#将筛选结果赋值给diamonds2,即得到一个新的数据集diamonds2。
这种方法一般不推荐,因为仅仅因为一个测量值无效,并不意味着所有测量值都无效。此外,如果您的数据质量较低,在对每个变量应用此方法之后,可能会发现几乎没有任何数据剩余下!因此,这种方法仅用于数据集很大,删除的少量行数据对总体样本的分析构成的影响极小。就如以上示例,总共有5万多行,我们删除了存在异常值的8行,对总体数据的影响很小。<br>
通常情况下,我们建议用缺失值替换异常值。最简单的方法是使用 mutate() 函数将变量替换为修改后的副本。可以使用 ifelse() <br>
diamonds3 <- diamonds %>%
mutate(y = ifelse(y < 3 | y > 20, NA, y))
对比异常值替换前后的数据
#输出原始数据的异常值
unusual <- diamonds %>%
filter(y < 3 | y > 20) %>%
select(price, x, y, z)
unusual
## # A tibble: 9 × 4
## price x y z
## <int> <dbl> <dbl> <dbl>
## 1 5139 0 0 0
## 2 6381 0 0 0
## 3 12210 8.09 58.9 8.06
## 4 12800 0 0 0
## 5 15686 0 0 0
## 6 18034 0 0 0
## 7 2075 5.15 31.8 5.12
## 8 2130 0 0 0
## 9 2130 0 0 0
#异常值替换为NA后,检测是否成功
unusual2 <- diamonds2 %>%
filter(y < 3 | y > 20) %>%
select(price, x, y, z)
unusual2
## # A tibble: 0 × 4
## # ℹ 4 variables: price <int>, x <dbl>, y <dbl>, z <dbl>
#查看原始数据中是否带有NA
na_rows <- diamonds[!complete.cases(diamonds), ]
print(na_rows)
## # A tibble: 0 × 10
## # ℹ 10 variables: carat <dbl>, cut <ord>, color <ord>, clarity <ord>,
## # depth <dbl>, table <dbl>, price <int>, x <dbl>, y <dbl>, z <dbl>
#检查异常值替换后具有NA的行
na_rows2 <- diamonds2[!complete.cases(diamonds2), ]
print(na_rows2)
## # A tibble: 0 × 10
## # ℹ 10 variables: carat <dbl>, cut <ord>, color <ord>, clarity <ord>,
## # depth <dbl>, table <dbl>, price <int>, x <dbl>, y <dbl>, z <dbl>
对于存在缺失值数据的分析或可视化:ggplot2认同缺失值不应该被悄无声息地丢弃的理念。尽管ggplot2不在图中包含它们,但在默认参数下会发出警告说明它们已被移除,例如:
ggplot(data = diamonds2, mapping = aes(x = x, y = y)) +
geom_point()
要抑制此警告,请将na.rm = TRUE。
ggplot(data = diamonds2, mapping = aes(x = x, y = y)) +
geom_point(na.rm = TRUE) #如果不设置,默认情况是na.rm = FALSE
1)在直方图中,缺失值会如何处理?在条形图中,缺失值会如何处理?
备注:条形图是用条形的“长度” 表示各类别频数的多少,其宽度(表示类别)则是固定的;直方图是用“面积” 表示各组频数的多少,矩形的高度表示每一组的频数或频率,宽度则表示各组的组距,因此其高度与宽度均有意义;由于分组数据具有连续性,直方图的各矩形通常是连续排列,而条形图则是分开排列。 条形图可用barplot(data),直方图可用hist(data, na.rm = TRUE)
2)在mean()和sum()函数中,na.rm = TRUE的作用是什么?
推荐:Robert I. Kabacoff. R语言实战(R in Action)第二版本, 王小宁等译,2016,P:182-197