探索性数据分析（EDA）
探索性数据分析exploratory data analysis
1 对分布进行可视化表示
分类变量在 R 中通常保存为因子或字符向量。要想检查分类变量的分布，可以使用条形图：
ggplot(data = diamonds) + geom_bar(mapping = aes(x = cut))
条形的高度表示每个 x 值中观测的数量，你可以使用 dplyr::count() 手动计算出这些值：
diamonds %>% count(cut)
要想检查连续变量的分布，可以使用直方图：
ggplot(data = diamonds) + geom_histogram(mapping = aes(x = carat), binwidth = 0.5)
最高的条形表示几乎有 30 000 个观测的 carat 值在 0.25 和 0.75 之间，这两个值分别是条形的左侧值和右侧值。可以通过 dplyr::count() 和 ggplot2::cut_width() 函数的组合来手动计算结果：
diamonds %>% count(cut_width(carat, 0.5))
如果只考虑重量小于 3 克拉的钻石，并选择一个更小的分箱宽度，那么直方图如下所示：
smaller <- diamonds %>% filter(carat < 3) ggplot(data = smaller, mapping = aes(x = carat)) + geom_histogram(binwidth = 0.1)
如果想要在同一张图上叠加多个直方图，那么我们建议你使用 geom_freqploy() 函数来代替
geom_histogram() 函数。 ggplot(data = smaller, mapping = aes(x = carat, color = cut)) + geom_freqpoly(binwidth = 0.1)
1.2 典型值
以下的直方图显示了美国黄石国家公园中的老忠实喷泉的 272 次喷发的时长（单位为分
钟）。喷发时间似乎聚集成了两组：短喷发（2 分钟左右）和长喷发（4~5 分钟），这两组
间几乎没有其他喷发时间：
ggplot(data = faithful, mapping = aes(x = eruptions)) + geom_histogram(binwidth = 0.25)
1.3 异常值
查看钻石数据集中 y 轴变量的分布，唯一能表示存在异常值的证据是，y 轴的取值范围出奇得宽：
ggplot(diamonds) + geom_histogram(mapping = aes(x = y), binwidth = 0.5)
正常值分箱中的观测太多了，以致于包括异常值的分箱高度太低，因此我们根本看不见
（如果仔细观察 x 轴 0 刻度附近，没准你能发现点什么）。为了更容易发现异常值，我们可以使用 coord_cartesian() 函数将 y 轴靠近 0 的部分放大：
ggplot(diamonds) + geom_histogram(mapping = aes(x = y), binwidth = 0.5) + coord_cartesian(ylim = c(0, 50))
（coord_cartesian() 函数中有一个用于放大 x 轴的 xlim() 参数。我们就可以看出有 3 个异常值，分别位于 0、30 左右和 60 左右。我们使用dplyr 将它们找出来：
unusual <- diamonds %>% filter(y < 3 | y > 20) %>% arrange(y) unusual
结果分析：y 变量测量钻石的三个维度之一，单位为毫米。我们知道钻石的宽度不可能是 0 毫米，因此这些值肯定是错误的。我们也完全可以认为 32 毫米和 59 毫米同样是令人难以置信的，这样的钻石长度超过 1 英寸（1 英寸 =2.54 厘米），简直就是无价之宝！
1.4 缺失值
如果在数据集中发现异常值，将带有可疑值的行全部丢弃：
diamonds2 <- diamonds %>% filter(between(y, 3, 20))
使用缺失值来代替异常值。最简单的做法就是使用 mutate() 函数创建一个新变量来代替原来的变量。你可以使用 ifelse() 函数将异常值替换为 NA：
diamonds2 <- diamonds %>% mutate(y = ifelse(y < 3 | y > 20, NA, y))
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)
在nycflights13::flights 中，dep_time 变量中的缺失值表示航班取消了。因此，你应该比较
一下已取消航班和未取消航班的计划出发时间。可以使用 is.na() 函数创建一个新变量来完成这个操作：
nycflights13::flights %>%
mutate(
cancelled = is.na(dep_time),
sched_hour = sched_dep_time %/% 100,
sched_min = sched_dep_time %% 100,
sched_dep_time = sched_hour + sched_min / 60
) %>%
ggplot(mapping = aes(sched_dep_time)) +
geom_freqpoly(
mapping = aes(color = cancelled),
binwidth = 1/4
)
2 相关变动
2.1 分组变量与连续变量
我们探索一下钻石价格是如何随着质量而变化的：
ggplot(data = diamonds, mapping = aes(x = price)) + geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
为了让比较变得更容易，需要改变 y 轴的显示内容，不再显示计数，而是显示密度。密度
是对计数的标准化，这样每个频率多边形下边的面积都是 1：
ggplot(
data = diamonds,
mapping = aes(x = price, y = ..density..)
) +
geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
结果分析：这张图的部分内容非常令人惊讶，其显示出一般钻石（质量最差）的平均价格是最高的！但这可能是因为频率多边形图很难解释，所以这张图还有很多可以改进的地方。
箱线图是对变量值分布的一种简单可视化表示，这种图在统计学家中非常流行。每张箱线图都包括以下内容。
使用 geom_boxplot() 函数查看按切割质量分类的价格分布：
ggplot(data = diamonds, mapping = aes(x = cut, y = price)) + geom_boxplot()
cut 是一个有序因子：“一般”不如“较好”、“较好”不如“很好”，以此类推。因为很多分类变量并没有这种内在的顺序，所以有时需要对其重新排序来绘制信息更丰富的图形。重新排序的其中一种方法是使用 reorder() 函数。
我们看一下 mpg 数据集中的 class 变量。你可能很想知道公路里程因汽车类别的不同会有怎样的变化：
ggplot(data = mpg, mapping = aes(x = class, y = hwy)) + geom_boxplot()
基于 hwy 值的中位数对 class 进行重新排序：
ggplot(data = mpg) +
geom_boxplot(
mapping = aes(
x = reorder(class, hwy, FUN = median),
y = hwy
)
)
如果变量名很长，那么将图形旋转 90 度效果会更好一些。你可以通过 coord_flip() 函数
完成这一操作：
ggplot(data = mpg) +
geom_boxplot(
mapping = aes(
x = reorder(class, hwy, FUN = median),
y = hwy
)
) +
coord_flip()
2.2 两个分类变量
要想对两个分类变量间的相关变动进行可视化表示，需要计算出每个变量组合中的观测数
量。完成这个任务的其中一种方法是使用内置的 geom_count() 函数：
ggplot(data = diamonds) + geom_count(mapping = aes(x = cut, y = color))
图中每个圆点的大小表示每个变量组合中的观测数量。
计算变量组合中的观测数量的另一种方法是使用 dplyr：
diamonds %>% count(color, cut)
接着使用 geom_tile() 函数和填充图形属性进行可视化表示：
diamonds %>%
count(color, cut) %>%
ggplot(mapping = aes(x = color, y = cut)) +
geom_tile(mapping = aes(fill = n))
2.3 两个连续变量
可以将相关变动看作点的模式。你可以看到钻石的克拉数和价值之间存在一种指数关系：
ggplot(data = diamonds) + geom_point(mapping = aes(x = carat, y = price), alpha = 1 / 100) #使用 alpha 图形属性添加透明度
geom_bin2d() 和 geom_hex() 函数将坐标平面分为二维分箱，并使用一种填充颜色表示落入
每个分箱的数据点。geom_bin2d() 创建长方形分箱。geom_hex() 创建六边形分箱。要想使
用 geom_hex()，需要安装 hexbin 包：
ggplot(data = smaller) + geom_bin2d(mapping = aes(x = carat, y = price))
library(hexbin) ggplot(data = smaller) + geom_hex(mapping = aes(x = carat, y = price))
以对 carat 进行分箱，然后为每个组生成一个箱线图：
ggplot(data = smaller, mapping = aes(x = carat, y = price)) + geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)))
3 模式和模型
的美国黄石国家公园中老忠实喷泉的喷发时长和两次喷发之间的等待时间做出一张散点图，该图会显示出一个模式：较长的等待时间与较长的喷发时间是相关的。图中还显示出两个簇：
ggplot(data = faithful) + geom_point(mapping = aes(x = eruptions, y = waiting))
以下代码拟合了一个模型，可以根据 carat 预测 price，并计算出残差（预测值和实际值之间的差别）。一旦去除克拉数对价格的影响，残差就能反映出钻石的价格：
library(modelr) mod <- lm(log(price) ~ log(carat), data = diamonds) diamonds2 <- diamonds %>% add_residuals(mod) %>% mutate(resid = exp(resid)) ggplot(data = diamonds2) + geom_point(mapping = aes(x = carat, y = resid))
