Table of Contents
引言
本章将向你展示如何利用可视化和转换系统性地探索数据，统计学家称之为探索性数据分析（简称EDA）。EDA是一个迭代循环。你：
- 对你的数据提出疑问。
- 通过可视化、转换和建模你的数据来寻找答案。
- 用你学到的东西来完善你的问题和/或生成新的问题。
EDA不是一个有严格规则的正式流程。EDA最重要的是心态。在EDA的初期阶段，你应该自由地探索所有想到的想法。其中一些想法会实现，有些则会是死胡同。随着探索的深入，你会聚焦于一些特别有成效的见解，最终将它们写出来并传达给他人。
EDA是任何数据分析的重要组成部分，即使主要研究问题是直接递给你的，因为你总需要调查数据的质量。数据清理只是EDA的一个应用：你会询问你的数据是否符合你的期望。要进行数据清理，你需要部署EDA的所有工具：可视化、转换和建模。
问题
“没有常规的统计问题，只有值得质疑的统计程序。”——大卫·考克斯爵士
“对正确问题的近似答案（往往模糊）远比对错误问题的准确答案要好得多，错误问题总能被精确化。”——约翰·图基
你在EDA期间的目标是对你的数据有更深入的理解。最简单的方法是用问题作为引导调查的工具。当你提问时，问题会让你关注数据集的某个特定部分，帮助你决定要制作哪些图表、模型或转换。
EDA本质上是一个创造性的过程。像大多数创意过程一样，提出高质量问题的关键是产生大量问题。在分析开始时很难提出揭示性的问题，因为你不知道从数据集中能获得哪些洞见。另一方面，每问一个新问题，都会让你接触到数据的新方面，增加你发现新事物的机会。如果你根据发现为每个问题提出新问题，就能快速深入挖掘数据中最有趣的部分，并设计出一系列发人深省的问题。
没有规定你应该问哪些问题来指导你的研究。然而，有两种类型的问题总是有助于你在数据中发现真相。你可以大致地这样表述这些问题：
- 我的变量中会发生什么样的变化？
- 我的变量之间会发生什么样的协变？
本章的其余部分将探讨这两个问题。我们将解释什么是变异和协变，并展示几种回答每个问题的方法。
变异
变异是指变量值在不同测量过程中变化的趋势。现实生活中很容易看到差异;如果你测量任意一个连续变量两次，结果会不同。即使测量的是恒定的量，比如光速，这一点也成立。你的每次测量都会包含一个小幅误差，且误差会因测量而异。变量也会变化，比如测量不同受试者（例如不同人的眼睛颜色）或不同时间（例如电子在不同时刻的能级）。每个变量都有其独特的变异模式，这可以揭示它在同一观测值以及跨观测值间的变化有趣信息。理解这种模式的最佳方法是可视化变量值的分布。
我们将从可视化钻石数据集中~54,000颗钻石的重量（克拉）分布开始探索。由于克拉是一个数值变量，我们可以使用直方图：
library(ggplot2)# 使用 ggplot2 对 diamonds 数据集中的 carat（克拉数）绘制直方图ggplot(diamonds, aes(x = carat)) + geom_histogram(binwidth = 0.5) # 设置箱宽为 0.5，表示每个柱子的区间宽度为 0.5 克拉
既然你已经能直观地看到变化，那么你应该在图表中关注什么？你应该问哪些后续问题？我们整理了一份清单，列出你在图表中会发现的最有用信息类型，并为每种信息类型准备了一些后续问题。提出好的后续问题的关键是依靠你的好奇心（你想了解更多？）以及怀疑态度（这怎么可能误导人呢？）。
经典数值
在条形图和直方图中，高条表示变量的常见值，较短的条表示较少的数值。没有条形的地方会显示出你数据中未曾见过的数值。为了将这些信息转化为有用的问题，请留意任何意想不到的情况：
- 哪些数值最常见？为什么？
- 哪些数值较为稀有？为什么？这符合你的期望吗？
- 你能看到什么异常的规律吗？什么能解释它们？
让我们看看小钻石的克拉分布情况。
library(dplyr)# 从 diamonds 数据集中筛选出克拉数小于 3 的钻石，保存为 smallersmaller <- diamonds |> filter(carat < 3)# 使用 smaller 数据集，按 carat 绘制直方图ggplot(smaller, aes(x = carat)) + geom_histogram(binwidth = 0.01) # 设置箱宽为 0.01，更细致地展示 carat 的分布
这张直方图提出了几个有趣的问题：
- 为什么整克拉和普通克拉比例的钻石会更多？
- 为什么每个峰稍微偏右的钻石比稍微左一点的多？
可视化还可以揭示聚类，表明数据中存在子群。要理解这些子群体，可以问：
- 每个子群体内的观察结果彼此之间有何相似之处？
- 分离聚类中的观测值彼此有何不同？
- 你如何解释或描述这些簇？
- 为什么聚类的外观可能会产生误导？
其中一些问题可以通过数据回答，而另一些则需要对数据的领域知识。许多问题会促使你探索变量之间的关系，例如，看看一个变量的取值是否可以解释另一个变量的行为。我们很快就会讨论到这一点。
异常值
异常值是指不寻常的观测值；也就是那些看起来不符合模式的数据点。有时异常值是数据录入错误，有时它们只是这次数据收集中碰巧观测到的极端值，而另一些时候，它们则暗示着重要的新发现。当数据量很大时，异常值有时很难在直方图中看出来。例如，看看 diamonds 数据集中 y 变量的分布。异常值存在的唯一证据，就是 x 轴范围异常宽。
# 使用 diamonds 数据集中的 y 变量绘制直方图ggplot(diamonds, aes(x = y)) + geom_histogram(binwidth = 0.5) # 设置箱宽为 0.5
常见区间中的观测值太多，以至于罕见区间显得非常短，这使得它们很难被看见（不过也许如果你盯着 0 看得足够仔细，能发现一些东西）。为了更容易看出这些不寻常的值，我们需要使用 coord_cartesian() 将 y 轴缩放到较小的值范围：
# 使用 diamonds 数据集中的 y 变量绘制直方图ggplot(diamonds, aes(x = y)) + geom_histogram(binwidth = 0.5) + # 设置箱宽为 0.5 coord_cartesian(ylim = c(0, 50)) # 将 y 轴显示范围限制在 0 到 50 之间
coord_cartesian() 也有一个 xlim() 参数，当你需要放大 x 轴时可以使用。ggplot2 也有 xlim() 和 ylim() 函数，但它们的工作方式略有不同：它们会丢弃超出范围的数据。
这使我们能够看到有三个不寻常的值：0、约 30 和约 60。我们用 dplyr 将它们提取出来：
# 从 diamonds 数据集中筛选出 y 小于 3 或大于 20 的异常值unusual <- diamonds |> filter(y < 3 | y > 20) |> select(price, x, y, z) |> arrange(y) # 按 y 升序排列# 查看筛选出的异常值数据unusual
y 变量测量的是这些钻石的三个尺寸之一，单位是毫米。我们知道钻石不可能有 0 毫米的宽度，所以这些值一定是错误的。通过做 EDA，我们发现了被编码为 0 的缺失数据，而如果只是简单地搜索 NA，我们永远不会发现这些数据。接下来，我们也许会选择将这些值重新编码为 NA，以避免产生误导性的计算。我们还可能怀疑 32 毫米和 59 毫米的测量值不合理：这些钻石都超过一英寸长，但价格却并没有高到几十万美元！
最好分别在有异常值和没有异常值的情况下重复你的分析。如果异常值对结果影响很小，而且你也弄不清它们为什么会出现，那么可以合理地省略它们，然后继续进行。然而，如果它们对你的结果有显著影响，你就不应该在没有充分理由的情况下删掉它们。你需要弄清楚是什么导致了这些异常值（例如，数据录入错误），并在你的书面报告中说明你已经将它们移除。
异常值
如果你在数据集中遇到了不寻常的值，并且只是想继续进行其余的分析，你有两个选择。
- 删除包含这些异常值的整行：
# 筛选 diamonds 数据集中 y 在 3 到 20 之间的观测值diamonds2 <- diamonds |> filter(between(y, 3, 20))我们不推荐这种做法，因为一个无效值并不意味着该观测中的其他所有值也都是无效的。此外，如果你的数据质量较低，当你把这种方法应用到每个变量时，可能会发现最后根本没有数据剩下了！
2. 相反，我们建议将这些不寻常的值替换为缺失值。最简单的方法是使用 mutate() 来用修改后的副本替换该变量。你可以使用 if_else() 函数将不寻常的值替换为 NA：
diamonds2 <- diamonds |> mutate(y = if_else(y < 3 | y > 20, NA, y))不太容易确定应该把缺失值画在什么位置，所以 ggplot2 不会将它们包含在图中，但它会警告你它们已被移除：
ggplot(diamonds2, aes(x = x, y = y)) + geom_point()
为抑制该警告，将 na.rm 设为 TRUE：
ggplot(diamonds2, aes(x = x, y = y)) + geom_point(na.rm = TRUE)有时候，你想了解带有缺失值的观测与有记录值的观测有什么不同。例如，在 nycflights13::flights1 中，dep_time 变量里的缺失值表示航班被取消了。因此，你可能想比较已取消和未取消航班的计划起飞时间。你可以通过创建一个新变量，并使用 is.na() 来检查 dep_time 是否缺失来实现这一点。
# install.packages("nycflights13")library(nycflights13)# 从 nycflights13::flights 航班数据中出发，# 创建一个新的数据版本，用于分析“是否取消”与“计划起飞时间”的关系nycflights13::flights |> mutate( # 新建逻辑变量 cancelled： # 如果 dep_time 是缺失值（NA），说明实际起飞时间不存在， # 通常表示该航班被取消了；否则为 FALSE cancelled = is.na(dep_time), # 将计划起飞时间 sched_dep_time 拆成“小时”和“分钟”两部分 # 例如 1345 表示 13 点 45 分 sched_hour = sched_dep_time %/% 100, sched_min = sched_dep_time %% 100, # 把原来的“时分格式”转换为十进制小时，便于绘图和比较 # 例如 13:30 转成 13.5 sched_dep_time = sched_hour + (sched_min / 60) ) |> # 以转换后的计划起飞时间为 x 轴， # 绘制频率多边形（frequency polygon），用于查看不同时间段的航班分布 # 按 cancelled 着色，从而对比“取消”和“未取消”的航班 ggplot(aes(x = sched_dep_time)) + geom_freqpoly(aes(color = cancelled), binwidth = 1/4)
不过，这张图并不好，因为未取消的航班远远多于已取消的航班。在下一节中，我们将探讨一些改进这种比较的方法。
缺失值在直方图中会怎样？缺失值在条形图中会怎样？为什么直方图和条形图对缺失值的处理方式不同？
- 直方图（histogram）：缺失值
NA一般会被自动忽略，不会画在图里；很多函数还会提示有缺失值被移除。 - 条形图（bar chart）：缺失值通常也不会自动作为一类显示，除非你先把它明确转换成一个类别，比如用
forcats::fct_explicit_na()或自己mutate()出一个"Missing"水平。
为什么不同？
因为：
- 直方图处理的是连续变量，要把数值分到区间里；
NA没有具体数值，无法归入任何区间。 - 条形图处理的是离散/分类变量，本质上是在数每个类别有多少个观测值；所以如果你愿意，可以把缺失值显式设成一个类别来统计。
na.rm = TRUE 在 mean() 和 sum() 中的作用是什么？
na.rm = TRUE 的作用是在计算前先移除缺失值 NA。
mean(x, na.rm = TRUE)：忽略NA后再计算平均值sum(x, na.rm = TRUE)：忽略NA后再求和
如果不写 na.rm = TRUE，只要数据里有 NA，结果通常就会变成 NA。
协变异
如果变异描述的是单个变量内部的行为，那么协变异描述的是变量之间的行为。协变异是指两个或多个变量的取值以相关方式共同变化的趋势。发现协变异的最佳方法是将两个或多个变量之间的关系可视化。
一个分类变量和一个数值变量
例如，让我们使用 geom_freqpoly() 来探索钻石价格如何随其质量（以切工衡量）而变化：
ggplot(diamonds, aes(x = price)) + # 绘制价格的频率多边形，并按切工 cut 着色 geom_freqpoly(aes(color = cut), binwidth = 500, linewidth = 0.75)
注意，ggplot2 对 cut 使用有序颜色刻度，因为它在数据中被定义为有序因子变量。
geom_freqpoly() 的默认外观在这里并不是很有用，因为由总计数决定的高度在不同 cut 之间差异太大，导致很难看出它们分布形状的差异。
为了简化比较，我们需要交换y轴显示的内容。我们不会显示计数，而是显示密度，也就是标准化的计数，使得每个频率多边形下的面积为一。
# 使用 diamonds 数据集，映射价格到 x 轴，密度到 y 轴ggplot(diamonds, aes(x = price, y = after_stat(density))) + # 按 cut 分组着色，绘制频率多边形；设置箱宽为 500，线宽为 0.75 geom_freqpoly(aes(color = cut), binwidth = 500, linewidth = 0.75)
注意，我们将密度映射到 y 轴，但由于密度并不是 diamonds 数据集中的变量，因此我们需要先计算它。为此，我们使用 after_stat() 函数。
这张图有点令人惊讶——看起来普通钻石（质量最低）的平均价格最高！但也许这是因为频率多边形有点难以解读——这张图里有很多内容。
一个视觉上更简单的情节是使用并排的箱线情节来探讨这段关系。
ggplot(diamonds, aes(x = cut, y = price)) + geom_boxplot()
我们看到的关于分布的信息少得多，但箱形图更紧凑，因此我们更容易比较它们（并且能在一张图中容纳更多内容）。这也支持了一个反直觉的发现：质量更好的钻石通常更便宜！在练习中，你将被挑战去弄清楚原因。
cut是一个有序因素：公平比好差，好比非常好差，等等。许多类别变量没有这样的内在顺序，所以你可能需要重新排序它们，以使显示更具信息量。其中一种方法是用fct_reorder（）。你可以在第16.4节了解更多关于这个功能的内容，但我们这里想给你一个快速预览，因为它非常有用。例如，取 mpg 数据集中的 class 变量。你可能会感兴趣，了解不同级别的高速公路里程差异：
# 使用 mpg 数据集，x 轴为车型类别 class，y 轴为高速公路油耗 hwyggplot(mpg, aes(x = class, y = hwy)) + # 绘制箱线图，查看不同车型类别的油耗分布 geom_boxplot()
library(forcats)# 使用 mpg 数据集，将 class 按 hwy 的中位数重新排序后映射到 x 轴# 将高速公路油耗 hwy 映射到 y 轴# 绘制箱线图，比较不同车型类别的 hwy 分布ggplot(mpg, aes(x = fct_reorder(class, hwy, median), y = hwy)) + geom_boxplot()
如果你的变量名很长，geom_boxplot（）如果把它翻转90°会更好。你可以通过交换x和y的美学映射来实现这一点。
ggplot(mpg, aes(x = hwy, y = fct_reorder(class, hwy, median))) + geom_boxplot()
两个分类变量
要可视化分类变量之间的协变关系，你需要统计这些分类变量各个水平组合的观测值数量。实现这一点的一种方法是使用内置的 geom_count()：
# 使用 diamonds 数据集，将 cut 映射到 x 轴，color 映射到 y 轴# geom_count() 会统计每种组合出现的次数，并用点的大小表示频数ggplot(diamonds, aes(x = cut, y = color)) + geom_count()
图中每个圆的大小表示在每组取值组合下有多少观测值。协变关系会表现为特定 x 值与特定 y 值之间的强相关。
探索这些变量之间关系的另一种方法是使用 dplyr 计算计数：
library(dplyr)library(ggplot2)diamonds |> count(color, cut)
然后使用 geom_tile() 和 fill 美学映射进行可视化：
diamonds |> count(color, cut) |> ggplot(aes(x = color, y = cut)) + geom_tile(aes(fill = n))
如果分类变量是无序的，你可能会想使用 seriation 包，同时重新排列行和列，以更清楚地揭示有趣的模式。对于更大的图，你可能会想尝试 heatmaply 包，它可以创建交互式图。
两个数值变量
你已经见过一种很好的方法来可视化两个数值变量之间的协变关系：使用 geom_point() 绘制散点图。你可以把协变关系看作点中的一种模式。例如，你可以看到钻石克拉大小和价格之间存在正相关：克拉越大的钻石，价格越高。这种关系是指数型的。
smaller <- diamonds |> filter(carat < 3)ggplot(smaller, aes(x = carat, y = price)) + geom_point()
随着数据集规模的增大，散点图会变得不那么有用，因为点会开始重叠，并堆积成一片均匀的黑色区域，这使得很难判断数据在二维空间中各处密度的差异，也很难看出趋势。你已经见过一种解决这个问题的方法：使用 alpha 美学映射来增加透明度。
ggplot(smaller, aes(x = carat, y = price)) + geom_point(alpha = 1 / 100)
但对于非常大的数据集，使用透明度可能会很有挑战性。另一种解决方案是使用分箱。此前你已经使用过 geom_histogram() 和 geom_freqpoly() 在一个维度上进行分箱。现在你将学习如何使用 geom_bin2d() 和 geom_hex() 在两个维度上进行分箱。
geom_bin2d() 和 geom_hex() 会将坐标平面划分为二维分箱，然后使用填充颜色来显示每个分箱中落入了多少点。geom_bin2d() 创建矩形分箱。geom_hex() 创建六边形分箱。你需要安装 hexbin 包才能使用 geom_hex()。
ggplot(smaller, aes(x = carat, y = price)) + geom_bin2d()ggplot(smaller, aes(x = carat, y = price)) + geom_hex()
另一种选择是对一个连续变量进行分箱，使其表现得像一个分类变量。然后，你就可以使用你已经学过的、用于可视化分类变量和连续变量组合的技术之一。例如，你可以对克拉进行分箱，然后为每个组绘制箱线图：
# 使用较小的数据集 smaller，绘制克拉（carat）与价格（price）的关系ggplot(smaller, aes(x = carat, y = price)) + # 按 carat 每 0.1 为一组进行分箱分组，绘制箱线图 geom_boxplot(aes(group = cut_width(carat, 0.1)))
cut_width(x, width)，如上所示，会将 x 按宽度为 width 的区间进行分箱。默认情况下，无论观测值有多少，箱线图看起来都大致相同（除了离群值的数量），因此很难看出每个箱线图汇总了不同数量的点。一种显示这一点的方法是将箱线图的宽度设为与点的数量成正比，使用 varwidth = TRUE。
发表评论