• R语言中的数量

    引言

    数值向量是数据科学的基石,而你在本书前面的部分已经多次使用过它们。现在,是时候系统地了解在 R 中可以对它们做些什么了,以确保你具备良好的基础来应对未来任何涉及数值向量的问题。

    我们将先向你介绍几个工具,帮助你在只有字符串时生成数字,然后再更详细地讲解 count()。接着,我们会深入介绍与 mutate() 很搭配的各种数值转换,包括一些更通用、可应用于其他类型向量的转换,不过它们通常还是与数值向量一起使用。最后,我们会介绍与 summarize() 很搭配的汇总函数,并展示它们也可以如何与 mutate() 一起使用。

    生成数字

    在大多数情况下,你拿到的数字已经以 R 的数值类型之一记录好了:整数或双精度数。不过,在某些情况下,你会把它们当作字符串遇到,可能是因为你通过将列标题转换生成了它们,或者是因为在数据导入过程中出了问题。

    readr 提供了两个有用的函数,用于将字符串解析为数字:parse_double() 和 parse_number()。当你有以字符串形式写出的数字时,请使用 parse_double()

    library(readr)
    x <- c("1.2", "5.6", "1e3")
    parse_double(x) # 需要readr包

    当字符串包含你想忽略的非数字文本时,请使用 parse_number()。这对于货币数据和百分比尤其有用:

    x <- c("$1,234", "USD 3,513", "59%")
    parse_number(x)

    count

    令人惊讶的是,仅靠计数和一点基础算术,你就能完成很多数据科学工作,因此 dplyr 努力让使用 count() 进行计数尽可能简单。这个函数非常适合在分析过程中进行快速探索和检查:

    library(nycflights13)
    flights |> count(dest)

    (我们通常把 count() 写在一行,因为它一般是在控制台里用来快速检查某个计算是否按预期工作的。)

    如果你想查看最常见的值,请添加 sort = TRUE:

    library(dplyr)
    flights |> count(dest, sort = TRUE)

    并且请记住,如果你想查看所有值,可以使用 |> View() 或 |> print(n = Inf)

    你可以使用 group_by()、summarize() 和 n() 手动执行相同的计算。这很有用,因为它允许你同时计算其他摘要:

    # 按目的地分组,统计每个目的地的航班数量,并计算平均到达延误时间
    flights |>
    group_by(dest) |>
    summarize(
    n = n(), # 计算每个目的地的航班总数
    delay = mean(arr_delay, na.rm = TRUE) # 计算平均到达延误,忽略缺失值
    )

    n() 是一个特殊的摘要函数,它不接受任何参数,而是访问关于“当前”组的信息。这意味着它只能在 dplyr 动词内部使用:

    n()

    n() 和 count() 有几个变体,你可能会觉得很有用:

    n_distinct(x) 计算一个或多个变量中不同(唯一)值的数量。例如,我们可以找出哪些目的地由最多的航空公司提供服务:

    # 按目的地分组
    flights |>
    group_by(dest) |>
    # 统计每个目的地对应的不同航空公司数量
    summarize(carriers = n_distinct(carrier)) |>
    # 按航空公司数量降序排列
    arrange(desc(carriers))

    加权计数就是求和。例如,你可以“统计”每架飞机飞行的英里数:

    # 按飞机编号分组
    flights |>
    group_by(tailnum) |>
    # 统计每架飞机的总飞行里程
    summarize(miles = sum(distance))

    加权计数是一个常见问题,因此 count() 有一个 wt 参数,可以完成同样的事情:

    # 按飞机编号统计,并对飞行里程 distance 求和
    flights |> count(tailnum, wt = distance)

    你可以将 sum() 和 is.na() 结合起来统计缺失值。在 flights 数据集中,这表示被取消的航班:

    # 按目的地分组
    flights |>
    group_by(dest) |>
    # 统计每个目的地中出发时间缺失的航班数,即取消的航班数
    summarize(n_cancelled = sum(is.na(dep_time)))

    数值变换

    变换函数与 mutate() 很适配,因为它们的输出长度与输入相同。绝大多数变换函数已经内置在基础 R 中。把它们全部列出来是不现实的,因此这一节只展示最有用的那些。比如,虽然 R 提供了你能想到的所有三角函数,但我们这里不列出它们,因为数据科学很少需要用到。

    算术与回收规则

    我们在第 2 章介绍了算术的基础知识(+-*/^),并且之后也多次使用它们。这些函数不需要太多解释,因为它们的作用就是你在小学学过的那些。不过,我们需要简单谈一下回收规则,它决定了当左侧和右侧的长度不同时时会发生什么。这一点对像 flights |> mutate(air_time = air_time / 60) 这样的操作很重要,因为 / 左边有 336,776 个数,而右边只有一个。

    R 通过回收,或者重复,较短的向量来处理长度不匹配的问题。我们可以通过在数据框外创建一些向量,更容易看到这一机制的运行:

    x <- c(1, 2, 10, 20)
    x / 5
    # is shorthand for
    x / c(5, 5, 5, 5)

    通常,你只想回收单个数字(即长度为 1 的向量),但 R 会回收任何较短长度的向量。一般来说(但并非总是),如果较长的向量不是较短向量长度的倍数,R 会给你一个警告:

    x * c(1, 2)
    x * c(1, 2, 3)

    这些回收规则也适用于逻辑比较(==<<=>>=!=),如果你不小心把 %in% 写成了 ==,而数据框的行数又碰巧不合适,就可能得到一个令人意外的结果。例如,看看这段试图找出 1 月和 2 月所有航班的代码:

    # 过滤出 1 月和 2 月的航班
    flights |>
    filter(month == c(1, 2))

    这段代码运行时不会报错,但它返回的结果并不是你想要的。由于回收规则的存在,它找到了奇数行中 1 月出发的航班,以及偶数行中 2 月出发的航班。遗憾的是,因为 flights 有偶数行,所以不会出现警告。

    为防止这类静默失败,大多数 tidyverse 函数使用更严格的回收方式,只回收单个值。不幸的是,这在这里并没有帮助,在许多其他情况下也一样,因为关键的计算是由基础 R 函数 == 完成的,而不是 filter()

    最小值和最大值

    df <- tribble(
    ~x, ~y,
    1, 3,
    5, 2,
    7, NA,
    )
    df |>
    mutate(
    min = pmin(x, y, na.rm = TRUE),
    max = pmax(x, y, na.rm = TRUE)
    )

    请注意,这些与摘要函数 min() 和 max() 不同,后者会接收多个观测值并返回单个值。若所有最小值和所有最大值都相同,你就能判断自己用了错误的形式:

    df |>
    mutate(
    min = min(x, y, na.rm = TRUE),
    max = max(x, y, na.rm = TRUE)
    )

    模运算

    模算术是这种运算的技术名称:你在学习小数点之前就已经做过了,也就是除法会得到一个整数和一个余数。在 R 中,%/% 表示整数除法,%% 用来计算余数:

    1:10 %/% 3 # 整除
    1:10 %% 3 # 取余

    模运算在 flights 数据集中很有用,因为我们可以用它把 sched_dep_time 变量拆解为小时和分钟:

    flights |>
    mutate(
    hour = sched_dep_time %/% 100, # 提取计划起飞时间中的“小时”部分
    minute = sched_dep_time %% 100, # 提取计划起飞时间中的“分钟”部分
    .keep = "used" # 只保留本次计算中用到的原始列
    )

    我们可以将其与第 12.4 节中的 mean(is.na(x)) 技巧结合起来,看看取消航班的比例如何随一天中的时间变化。结果如图 13.1 所示。

    library(ggplot2)
    flights |>
    group_by(hour = sched_dep_time %/% 100) |> # 按计划起飞时间的“小时”分组
    summarize(prop_cancelled = mean(is.na(dep_time)), # 计算每小时取消航班的比例
    n = n()) |> # 统计每小时的航班数量
    filter(hour > 1) |> # 过滤掉 1 点及以前的数据
    ggplot(aes(x = hour, y = prop_cancelled)) + # 以小时为 x 轴,取消比例为 y 轴作图
    geom_line(color = "grey50") + # 绘制灰色折线
    geom_point(aes(size = n)) # 绘制点,并用点大小表示航班数量

    对数

    对数是一种非常有用的变换,适用于处理跨越多个数量级的数据,并将指数增长转换为线性增长。在 R 中,你可以选择三种对数函数:log()(自然对数,以 e 为底)、log2()(以 2 为底)和 log10()(以 10 为底)。我们建议使用 log2() 或 log10()log2() 很容易解释,因为在对数尺度上相差 1,对应于原始尺度上翻倍;相差 -1,则对应于减半;而 log10() 则便于反向转换,因为例如 3 表示 103=1000log() 的逆函数是 exp();要计算 log2() 或 log10() 的逆函数,则需要使用 2^ 或 10^

    四舍五入

    使用 round(x) 将数字四舍五入到最接近的整数:

    round(123.456)

    你可以使用第二个参数 digits 来控制四舍五入的精度。round(x, digits) 会四舍五入到最接近的 10^-n,因此 digits = 2 时会四舍五入到最接近的 0.01。这个定义很有用,因为它意味着 round(x, -3) 会四舍五入到最接近的千位,而它确实如此:

    round(123.456, 2) # two digits
    round(123.456, 1) # one digit
    round(123.456, -1) # round to nearest ten
    round(123.456, -2) # round to nearest hundred

    round() 有一个看起来一开始有些出乎意料的怪现象:

    round(c(1.5, 2.5))

    round() 使用的是所谓的“向偶数舍入”或“银行家舍入”:如果一个数字正好处于两个整数的中间,它会被舍入到偶数的那个整数。这是一个很好的策略,因为它能保持舍入结果不偏不倚:所有 0.5 中有一半会向上舍入,另一半会向下舍入。

    round() 与 floor() 和 ceiling() 成对使用:floor() 总是向下舍入,而 ceiling() 总是向上舍入:

    x <- 123.456
    floor(x)
    ceiling(x)

    将数字按范围分组

    使用 cut() 将数值向量划分(即分箱)为离散的区间:

    x <- c(1, 2, 5, 10, 15, 20)
    cut(x, breaks = c(0, 5, 10, 15, 20))

    断点不需要均匀分布:

    cut(x, breaks = c(0, 5, 10, 100))
    # 将 x 按指定断点分组,并为各区间指定标签
    cut(x,
    breaks = c(0, 5, 10, 15, 20),
    labels = c("sm", "md", "lg", "xl")
    )

    范围之外的任何值都会变成 NA:

    y <- c(NA, -10, 5, 10, 30)
    cut(y, breaks = c(0, 5, 10, 15, 20))

    查看文档了解其他有用的参数,例如 right 和 include.lowest,它们用于控制区间是 [a, b) 还是 (a, b],以及是否应将最低区间设为 [a, b]

    累积和滚动聚合

    Base R 提供了 cumsum()cumprod()cummin()cummax(),用于计算累计和、累计积、累计最小值和累计最大值。dplyr 提供了 cummean() 用于计算累计平均值。累计和在实际中最常见:

    x <- 1:10
    cumsum(x)

    通用变换

    以下各节描述了一些常用的通用变换,这些变换通常用于数值向量,但也可应用于所有其他列类型。

    排名

    dplyr 提供了许多受 SQL 启发的排名函数,但你应该始终从 dplyr::min_rank() 开始。它使用处理并列项的典型方法,例如:第 1、2、2、4。

    x <- c(1, 5, 5, 17, 22, NA)
    min_rank(x)

    请注意,最小的值会获得最低的排名;使用 desc(x) 可让最大的值获得最小的排名:

    min_rank(desc(x))

    如果 min_rank() 不能满足你的需求,可以看看 dplyr::row_number()dplyr::dense_rank()dplyr::percent_rank() 和 dplyr::cume_dist() 这些变体。详情请参阅文档。

    # 创建包含向量 x 的 tibble
    df <- tibble(x = x)
    # 计算不同的排名指标
    df |>
    mutate(
    # 按 x 的顺序生成行号排名
    row_number = row_number(x),
    # 密集排名:相同值共享同一名次,且名次不跳号
    dense_rank = dense_rank(x),
    # 百分位排名:返回每个值在分布中的相对位置
    percent_rank = percent_rank(x),
    # 累积分布:返回小于等于当前值的比例
    cume_dist = cume_dist(x)
    )

    你可以通过为 base R 的 rank() 选择合适的 ties.method 参数来实现许多相同的结果;你可能还希望将 na.last = "keep" 设为保留 NA 为 NA

    row_number() 也可以在 dplyr 动词内部不带任何参数使用。在这种情况下,它会给出“当前”行的编号。当它与 %% 或 %/% 结合使用时,可作为将数据划分为大小相近组别的有用工具:

    # 创建一个包含 1 到 10 的 id 列的 tibble
    df <- tibble(id = 1:10)
    # 基于行号进行分组计算
    df |>
    mutate(
    # 从 0 开始的行号
    row0 = row_number() - 1,
    # 将行号按 3 取余,得到 3 个循环分组
    three_groups = row0 %% 3,
    # 每 3 行划分为一组
    three_in_each_group = row0 %/% 3
    )

    偏移量

    dplyr::lead() 和 dplyr::lag() 允许你引用紧邻“当前”值之前或之后的值。它们返回与输入长度相同的向量,并在开头或结尾用 NA 填充:

    # 创建一个数值向量 x
    x <- c(2, 5, 11, 11, 19, 35)
    # lag(x) 返回前一个值,首位用 NA 填充
    lag(x)
    # lead(x) 返回后一个值,末位用 NA 填充
    lead(x)

    x - lag(x) 给出了当前值与前一个值之间的差值。

    x - lag(x)

    x == lag(x) 告诉你当前值何时发生变化。

    # 判断当前值是否与前一个值相同
    x == lag(x)

    你可以通过使用第二个参数 n 来向前或向后移动多个位置。

    连续标识符

    events <- tibble(
    time = c(0, 1, 2, 3, 5, 10, 12, 15, 17, 19, 20, 27, 28, 30)
    )

    并且你已经计算出每个事件之间的时间间隔,并判断是否存在足够大的间隔来满足条件:

    # 计算每个事件与前一个事件之间的时间差,并判断是否存在至少 5 分钟的间隔
    events <- events |>
    mutate(
    diff = time - lag(time, default = first(time)),
    has_gap = diff >= 5
    )
    events

    但我们如何从这个逻辑向量得到可以用于 group_by() 的结果呢?第 13.4.7 节介绍的 cumsum() 就派上用场了:当出现间隔,也就是 has_gap 为 TRUE 时,组号就会加一(第 12.4.2 节):

    # 根据 has_gap 的 TRUE 出现次数累计分组编号
    events |> mutate(
    group = cumsum(has_gap)
    )

    创建分组变量的另一种方法是使用 consecutive_id(),它会在其参数中的某个值发生变化时开始一个新的组。例如,受这个 Stack Overflow 问题的启发,假设你有一个包含一堆重复值的数据框:

    df <- tibble(
    x = c("a", "a", "a", "b", "c", "c", "d", "e", "a", "a", "b", "b"),
    y = c(1, 2, 3, 2, 4, 1, 3, 9, 4, 8, 10, 199)
    )

    如果你想保留每个重复的 x 的第一行,可以使用 group_by()consecutive_id() 和 slice_head()

    # 按连续相同的 x 值分组,并保留每组的第一行
    df |>
    group_by(id = consecutive_id(x)) |>
    slice_head(n = 1)

    数值摘要

    仅使用我们已经介绍过的计数、均值和求和,就能让你走很远,但 R 还提供了许多其他有用的摘要函数。下面列出一些你可能会觉得有用的函数。

    居中

    到目前为止,我们主要使用 mean() 来概括一组数值的中心。正如我们在第 3.6 节中看到的那样,由于均值是总和除以数量,因此即使只有几个异常高或异常低的值,它也会受到影响。另一种方法是使用 median(),它会找到位于向量“中间”的一个值,也就是说,50% 的值高于它,50% 的值低于它。根据你所关注变量的分布形状,均值或中位数可能是更好的中心度量。例如,对于对称分布,我们通常报告均值;而对于偏态分布,我们通常报告中位数。

    图 13.2 比较了每个目的地的平均起飞延误时间与中位起飞延误时间(单位:分钟)。中位延误总是小于平均延误,因为航班有时会晚点几个小时起飞,但从不会提前几个小时起飞。

    # 按年、月、日分组,计算每天的出发延误平均值、中位数和航班数量
    flights |>
    group_by(year, month, day) |>
    summarize(
    mean = mean(dep_delay, na.rm = TRUE),
    median = median(dep_delay, na.rm = TRUE),
    n = n(),
    .groups = "drop"
    ) |>
    # 绘制平均延误与中位延误的散点图
    ggplot(aes(x = mean, y = median)) +
    # 添加斜率为 1 的参考线
    geom_abline(slope = 1, intercept = 0, color = "white", linewidth = 2) +
    # 绘制每一天的点
    geom_point()

    你可能还会想到众数,也就是最常见的值。这种摘要方法只适用于非常简单的情况(这也是你可能在中学里学过它的原因),但对许多真实数据集来说并不好用。如果数据是离散的,可能会有多个最常见值;如果数据是连续的,则可能根本不存在最常见值,因为每个值都略有不同。基于这些原因,统计学家通常不太使用众数,而且 base R 中也没有包含 mode 函数。

    最小值、最大值和分位数

    如果你感兴趣的是中心以外的位置呢?min() 和 max() 会分别给出最大值和最小值。另一个强大的工具是 quantile(),它是中位数的推广:quantile(x, 0.25) 会找出 x 中大于 25% 数值的那个值,quantile(x, 0.5) 等同于中位数,而 quantile(x, 0.95) 会找出大于 95% 数值的那个值。

    对于 flights 数据,你可能更想查看延误的 95% 分位数,而不是最大值,因为它会忽略最严重延误的 5% 航班,而这些航班的延误可能非常极端。

    # 按年、月、日分组,计算每天出发延误的最大值和 95% 分位数
    flights |>
    group_by(year, month, day) |>
    summarize(
    max = max(dep_delay, na.rm = TRUE),
    q95 = quantile(dep_delay, 0.95, na.rm = TRUE),
    .groups = "drop"
    )

    离散程度

    我们可以用它来揭示 flights 数据中的一个小异常。你可能会以为起点和终点之间的距离分散程度应该为 0,因为机场总是在同一个位置。但下面的代码揭示了机场 EGE 的一个数据异常:

    # 按起点和终点分组,计算每条航线的距离四分位距和航班数量
    # 再筛选出距离四分位距大于 0 的航线,找出距离不一致的数据
    flights |>
    group_by(origin, dest) |>
    summarize(
    distance_iqr = IQR(distance),
    n = n(),
    .groups = "drop"
    ) |>
    filter(distance_iqr > 0)

    分布

    值得记住的是,上面描述的所有汇总统计量,都是把分布压缩成一个单一数字的方式。这意味着它们本质上是有简化性的,如果你选错了汇总方式,就很容易忽略组间的重要差异。这就是为什么,在确定汇总统计量之前,先把分布可视化总是一个好主意。

    图 13.3 展示了出发延误的整体分布。该分布严重偏斜,以至于我们不得不放大才能看清数据的大部分。这表明均值不太可能是一个好的摘要统计量,我们或许更倾向于使用中位数。

    另外,检查子组的分布是否与整体相似也是个好主意。在下面的图中,叠加了 365 条 dep_delay 的频率多边形,每条对应一天。各个分布看起来遵循相同的模式,这说明对每一天使用相同的摘要统计量是可行的。

    # 过滤出起飞延误小于 120 分钟的航班
    # 按“月 + 日”组合分组,绘制每一天的起飞延误频率多边形
    # 用较小的 binwidth 让分布更平滑,alpha 用于叠加时的透明度控制
    flights |>
    filter(dep_delay < 120) |>
    ggplot(aes(x = dep_delay, group = interaction(day, month))) +
    geom_freqpoly(binwidth = 5, alpha = 1/5)

    不要害怕为你正在处理的数据探索一些专门定制的摘要统计量。在这种情况下,这可能意味着分别汇总提前起飞的航班和延误起飞的航班;或者由于这些值严重偏斜,你也可以尝试进行对数变换。最后,不要忘记你在第 3.6 节学到的内容:无论何时创建数值摘要,最好都把每个组中的观测数量一起包含进去。

    位置

    还有一种最后的摘要类型对数值向量很有用,但它也适用于其他任何类型的值:提取特定位置上的值:first(x)last(x) 和 nth(x, n)

    例如,我们可以找出每天的第一班、第五班和最后一班出发航班:

    # 按年、月、日分组
    # 汇总每一天的起飞时间:
    # first_dep:当天第一个起飞时间(忽略缺失值)
    # fifth_dep:当天第五个起飞时间(忽略缺失值)
    # last_dep:当天最后一个起飞时间(忽略缺失值)
    flights |>
    group_by(year, month, day) |>
    summarize(
    first_dep = first(dep_time, na_rm = TRUE),
    fifth_dep = nth(dep_time, 5, na_rm = TRUE),
    last_dep = last(dep_time, na_rm = TRUE)
    )

    (注意:由于 dplyr 函数使用下划线 _ 来分隔函数名和参数名,因此这些函数使用 na_rm,而不是 na.rm。)

    如果你熟悉 [,我们将在第 27.2 节回到它,那么你可能会想,自己是否还需要这些函数。原因有三点:default 参数允许你在指定位置不存在时提供默认值,order_by 参数允许你在局部范围内覆盖行的顺序,而 na_rm 参数允许你删除缺失值。

    按位置提取值与按排名筛选是互补的。筛选会给出所有变量,并将每个观测值放在单独的一行中:

    # 按年、月、日分组
    # 计算每组中计划起飞时间的排名 r
    # 保留每组中排名最小(最早)和排名最大(最晚)的航班
    flights |>
    group_by(year, month, day) |>
    mutate(r = min_rank(sched_dep_time)) |>
    filter(r %in% c(1, max(r)))

    使用mutate()

    顾名思义,汇总函数通常与 summarize() 搭配使用。不过,由于我们在第 13.4.1 节讨论过的回收规则,它们也可以很有用地与 mutate() 搭配,尤其是在你想做某种分组标准化时。例如:

    • x / sum(x) 计算某个总量的占比。
    • (x – mean(x)) / sd(x) 计算 Z 分数(标准化为均值 0、标准差 1)。
    • (x – min(x)) / (max(x) – min(x)) 将数值标准化到 [0, 1] 范围。
    • x / first(x) 根据第一个观测值计算指数。
  • R语言的逻辑向量

    https://wp.me/p80aHo-2Ef

    引言

    在这一章中,你将学习处理逻辑向量的工具。逻辑向量是最简单的向量类型,因为每个元素只能取三种可能的值之一:TRUE、FALSE 和 NA。在你的原始数据中,逻辑向量相对少见,但在几乎每一次分析过程中,你都会创建和操作它们。

    我们将从讨论创建逻辑向量最常见的方法开始:使用数值比较。然后你将学习如何使用布尔代数来组合不同的逻辑向量,以及一些有用的汇总方法。最后,我们会介绍 if_else() 和 case_when(),这两个有用的函数可以借助逻辑向量进行条件变更。

    先修知识

    随着我们开始介绍更多工具,并不总能找到一个完美的真实示例。因此,我们将开始使用 c() 来虚构一些示例数据:

    x <- c(1, 2, 3, 5, 7, 11, 13)
    x * 2
    #> [1] 2 4 6 10 14 22 26

    这会更容易解释各个函数,但代价是更难看出它如何应用到你的数据问题中。只需记住,我们对一个独立向量所做的任何操作,你都可以借助 mutate() 及其相关函数对数据框中的变量进行相同的处理。

    library(dplyr)
    df <- tibble(x)
    df |>
    mutate(y = x * 2)

    比较

    创建逻辑向量的一种非常常见的方法,是使用 <、<=、>、>=、!= 和 == 进行数值比较。到目前为止,我们主要是在 filter() 中临时创建逻辑变量——它们被计算、使用,然后就被丢弃了。例如,下面的 filter() 会找出所有大致准时到达的白天出发航班:

    flights |>
    filter(dep_time > 600 & dep_time < 2000 & abs(arr_delay) < 20)

    了解这一点很有用:这是一种快捷方式,而你也可以使用 mutate() 显式创建底层的逻辑变量:

    flights |>
    mutate(
    daytime = dep_time > 600 & dep_time < 2000, # 判断是否为白天出发:起飞时间在 06:00 到 20:00 之间
    approx_ontime = abs(arr_delay) < 20, # 判断是否大致准点:到达延误绝对值小于 20 分钟
    .keep = "used" # 只保留本次计算中用到的列
    )

    这对于更复杂的逻辑尤其有用,因为为中间步骤命名,不仅能让代码更易读,也更容易检查每一步是否都已正确计算。

    总的来说,最初的 filter() 等同于:

    flights |>
    mutate(
    daytime = dep_time > 600 & dep_time < 2000, # 判断是否为白天出发:起飞时间在 06:00 到 20:00 之间
    approx_ontime = abs(arr_delay) < 20, # 判断是否大致准点:到达延误绝对值小于 20 分钟
    ) |>
    filter(daytime & approx_ontime) # 筛选出同时满足“白天出发”和“大致准点”的航班

    浮点数比较

    注意不要对数字使用 ==。例如,看起来这个向量包含数字 1 和 2:

    x <- c(1 / 49 * 49, sqrt(2) ^ 2) # 计算两个看似应为 1 和 2 的浮点数表达式,可能因精度问题得到近似值
    x

    但如果你测试它们是否相等,结果会是 FALSE:

    x == c(1, 2)

    怎么回事?计算机用固定数量的小数位来存储数字,因此无法精确表示 1/49 或 sqrt(2),后续计算会有非常轻微的偏差。我们可以通过将 digits 参数传给 print() 来查看精确值:

    print(x, digits = 16)

    既然你已经看到 == 为什么会失效,那该怎么办呢?一个选择是使用 dplyr::near(),它会忽略微小差异:

    library(dplyr)
    near(x, c(1, 2))

    缺失值

    缺失值表示未知,因此它们具有“传染性”:几乎任何涉及未知值的操作结果也会是未知:

    NA > 5
    10 == NA

    最令人困惑的结果是这一条:

    NA == NA

    如果我们人为提供一点更多上下文,就最容易理解为什么这是真的:

    # We don't know how old Mary is
    age_mary <- NA
    # We don't know how old John is
    age_john <- NA
    # Are Mary and John the same age?
    age_mary == age_john

    所以,如果你想找出所有 dep_time 缺失的航班,下面的代码不会起作用,因为 dep_time == NA 会对每一行都返回 NA,而 filter() 会自动丢弃缺失值:

    library(nycflights13)
    flights |>
    filter(dep_time == NA)

    我们需要一个新工具:is.na()

    is.na()

    is.na(x) 适用于任何类型的向量,并且会对缺失值返回 TRUE,对其他所有值返回 FALSE

    is.na(c(TRUE, NA, FALSE))
    is.na(c(1, NA, 3))
    is.na(c("a", NA, "b"))

    我们可以使用 is.na() 找出所有 dep_time 缺失的行:

    flights |>
    filter(is.na(dep_time))

    is.na() 也可以在 arrange() 中派上用场。arrange() 通常会把所有缺失值放到最后,但你可以先按 is.na() 排序,从而覆盖这一默认行为:

    flights |>
    filter(month == 1, day == 1) |>
    arrange(dep_time)
    flights |>
    filter(month == 1, day == 1) |>
    arrange(desc(is.na(dep_time)), dep_time)

    布尔代数

    一旦你有了多个逻辑向量,就可以使用布尔代数将它们组合起来。在 R 中,& 表示“与”,| 表示“或”,! 表示“非”,而 xor() 表示异或。例如,df |> filter(!is.na(x)) 会找到所有 x 不缺失的行,而 df |> filter(x < -10 | x > 0) 会找到所有 x 小于 -10 或大于 0 的行。图 12.1 展示了常用布尔运算的示例及其工作方式。

    除了 & 和 | 之外,R 还有 && 和 ||。不要在 dplyr 函数中使用它们!它们被称为短路运算符,并且只会返回单个 TRUE 或 FALSE。它们对编程很重要,但不适用于数据科学。

    缺失值

    布尔代数中缺失值的规则有点难以解释,因为乍一看它们似乎是不一致的:

    df <- tibble(x = c(TRUE, FALSE, NA))
    df |>
    mutate(
    and = x & NA,
    or = x | NA
    )

    要理解这是怎么回事,可以想一想 NA | TRUE(NA 或 TRUE)。逻辑向量中的缺失值表示这个值可能是 TRUE 也可能是 FALSE。TRUE | TRUE 和 FALSE | TRUE 都是 TRUE,因为它们中至少有一个为 TRUE。NA | TRUE 也必须是 TRUE,因为 NA 既可能是 TRUE,也可能是 FALSE。然而,NA | FALSE 是 NA,因为我们不知道 NA 是 TRUE 还是 FALSE。对于 &,类似的推理也适用,只不过这时两个条件都必须满足。因此,NA & TRUE 是 NA,因为 NA 既可能是 TRUE 也可能是 FALSE,而 NA & FALSE 是 FALSE,因为至少有一个条件为 FALSE。

    运算顺序

    请注意,运算顺序并不像英语那样起作用。看下面这段代码,它用于找出所有在 11 月或 12 月起飞的航班:

    flights |>
    filter(month == 11 | month == 12)

    你可能会想像用英语那样来写:“找出所有在 11 月或 12 月起飞的航班。”:

    flights |>
    filter(month == 11 | 12)

    这段代码不会报错,但看起来也没有达到预期效果。这是怎么回事?

    在这里,R 会先计算 month == 11,生成一个逻辑向量,我们把它叫做 nov。然后它再计算 nov | 12

    当你把数字和逻辑运算符一起使用时,R 会把除了 0 以外的所有数都转换成 TRUE。所以 12 会被当作 TRUE,因此这就等价于 nov | TRUE

    而 nov | TRUE 的结果永远是 TRUE,所以最终会选中每一行

    flights |>
    mutate(
    nov = month == 11,
    final = nov | 12,
    .keep = "used"
    )

    %in%

    避免把 == 和 | 的顺序弄错的一个简单方法是使用 %in%x %in% y 会返回一个与 x 长度相同的逻辑向量,只要 x 中的某个值出现在 y 中,就会是 TRUE

    1:12 %in% c(1, 5, 11)
    letters[1:10] %in% c("a", "e", "i", "o", "u")

    所以,要找出所有在 11 月和 12 月起飞的航班,我们可以这样写:

    flights |>
    filter(month %in% c(11, 12))

    请注意,%in% 对 NA 的处理规则与 == 不同,因为 NA %in% NA 的结果是 TRUE

    c(1, 2, NA) == NA
    c(1, 2, NA) %in% NA

    这可以成为一个有用的快捷写法:

    flights |>
    filter(dep_time %in% c(NA, 0800))

    汇总

    以下各节将介绍一些用于汇总逻辑向量的有用技巧。除了只适用于逻辑向量的函数外,你还可以使用适用于数值向量的函数。

    逻辑汇总

    有两个主要的逻辑汇总函数:any() 和 all()any(x) 相当于 |;如果 x 中存在任何 TRUE,它就会返回 TRUEall(x) 相当于 &;只有当 x 的所有值都是 TRUE 时,它才会返回 TRUE。和大多数汇总函数一样,你可以通过 na.rm = TRUE 来忽略缺失值。

    例如,我们可以使用 all() 和 any() 来判断是否所有航班的起飞延误都不超过一小时,或者是否有航班的到达延误达到了五小时或以上。而使用 group_by() 则可以让我们按天进行这样的分析:

    # 按年、月、日分组,汇总每天的航班延误情况
    flights |>
    group_by(year, month, day) |>
    summarize(
    # 如果当天所有航班的起飞延误都不超过 60 分钟,则为 TRUE
    all_delayed = all(dep_delay <= 60, na.rm = TRUE),
    # 如果当天有任一航班的到达延误达到或超过 300 分钟,则为 TRUE
    any_long_delay = any(arr_delay >= 300, na.rm = TRUE),
    # 汇总后取消分组
    .groups = "drop"
    )

    不过在大多数情况下,any() 和 all() 都有些过于粗略;如果能更详细地了解有多少值为 TRUE 或 FALSE 就好了。这就引出了数值汇总。

    逻辑向量的数值汇总

    当你在数值上下文中使用逻辑向量时,TRUE 会变成 1FALSE 会变成 0。这使得 sum() 和 mean() 在处理逻辑向量时非常有用,因为 sum(x) 会给出 TRUE 的个数,而 mean(x) 会给出 TRUE 的比例(因为 mean() 其实就是 sum() 除以 length())。

    这就让我们能够,例如,看到出发延误最多不超过一小时的航班所占比例,以及到达延误达到五小时或以上的航班数量:

    flights |>
    group_by(year, month, day) |>
    summarize(
    proportion_delayed = mean(dep_delay <= 60, na.rm = TRUE),
    count_long_delay = sum(arr_delay >= 300, na.rm = TRUE),
    .groups = "drop"
    )

    逻辑子集选择

    逻辑向量在汇总中的最后一种用途是:你可以用逻辑向量将单个变量筛选为感兴趣的子集。

    假设我们想只查看那些确实发生了延误的航班的平均延误时间。一种做法是先筛选出这些航班,然后再计算平均延误时间:

    # 按日期分组,筛选出到达延误大于 0 的航班
    flights |>
    filter(arr_delay > 0) |>
    # 按年、月、日分组
    group_by(year, month, day) |>
    # 计算每一天的平均到达延误和航班数量
    summarize(
    behind = mean(arr_delay),
    n = n(),
    .groups = "drop"
    )

    这就导致:

    # 按年、月、日分组,分别计算:
    # 1. 仅对到达延误大于 0 的航班求平均延误时间
    # 2. 仅对到达延误小于 0 的航班求平均提前时间
    # 3. 每组航班总数
    flights |>
    group_by(year, month, day) |>
    summarize(
    behind = mean(arr_delay[arr_delay > 0], na.rm = TRUE),
    ahead = mean(arr_delay[arr_delay < 0], na.rm = TRUE),
    n = n(),
    .groups = "drop"
    )

    还要注意组大小的差异:在第一个分组中,n() 表示每天延误的航班数量;而在第二个分组中,n() 表示航班总数。

    条件变换

    逻辑向量最强大的功能之一是它们可用于条件变换,也就是在条件 x 下做一件事,而在条件 y 下做另一件事。实现这一点有两个重要工具:if_else() 和 case_when()

    if_else()

    如果你想在条件为 TRUE 时使用一个值,而在条件为 FALSE 时使用另一个值,可以使用 dplyr::if_else()。你总是会用到 if_else() 的前三个参数。第一个参数 condition 是一个逻辑向量,第二个参数 true 表示条件为真时的输出,第三个参数 false 表示条件为假时的输出。

    先从一个简单的例子开始,把一个数值向量标记为“+ve”(正)或“-ve”(负):

    x <- c(-3:3, NA)
    if_else(x > 0, "+ve", "-ve")

    还有一个可选的第四个参数 missing,当输入为 NA 时会使用它:

    if_else(x > 0, "+ve", "-ve", "???")

    你也可以将向量用于 true 和 false 参数。例如,这使我们能够创建一个 abs() 的最简实现:

    if_else(x < 0, -x, x)

    到目前为止,所有参数都使用了相同的向量,但当然你也可以混合搭配。例如,你可以像这样实现一个简单版本的 coalesce()

    你可能已经注意到上面的标记示例中有一个小问题:零既不是正数,也不是负数。我们可以通过再添加一个 if_else() 来解决这个问题:

    if_else(x == 0, "0", if_else(x < 0, "-ve", "+ve"), "???")

    这已经有点难读了,而且你可以想象,如果条件更多,情况只会更难读。相反,你可以改用 dplyr::case_when()

    case_when()

    dplyr 的 case_when() 受 SQL 的 CASE 语句启发,提供了一种根据不同条件执行不同计算的灵活方式。它有一种特殊的语法,不幸的是,这种语法在 tidyverse 中你几乎不会在别处见到。它接受看起来像 condition ~ output 的成对表达式。condition 必须是一个逻辑向量;当它为 TRUE 时,就会使用 output

    这意味着我们可以按如下方式重现之前嵌套的 if_else()

    x <- c(-3:3, NA)
    case_when(
    x == 0 ~ "0",
    x < 0 ~ "-ve",
    x > 0 ~ "+ve",
    is.na(x) ~ "???"
    )

    这段代码更多,但也更明确。

    为了解释 case_when() 是如何工作的,我们先来看看一些更简单的情况。如果没有任何条件匹配,输出就会得到一个 NA

    case_when(
    x < 0 ~ "-ve",
    x > 0 ~ "+ve"
    )

    如果你想创建一个“默认”/兜底值,可以使用 .default

    case_when(
    x < 0 ~ "-ve",
    x > 0 ~ "+ve",
    .default = "???"
    )

    并且请注意,如果多个条件都匹配,只会使用第一个:

    case_when(
    x > 0 ~ "+ve",
    x > 2 ~ "big"
    )

    就像 if_else() 一样,你可以在 ~ 的两边都使用变量,并且可以根据需要灵活组合变量来解决你的问题。例如,我们可以使用 case_when() 为到达延误提供一些便于人类阅读的标签:

    flights |>
    mutate(
    status = case_when(
    is.na(arr_delay) ~ "cancelled",
    arr_delay < -30 ~ "very early",
    arr_delay < -15 ~ "early",
    abs(arr_delay) <= 15 ~ "on time",
    arr_delay < 60 ~ "late",
    arr_delay < Inf ~ "very late",
    ),
    .keep = "used"
    )

    在编写这种复杂的 case_when() 语句时要小心;我最初的两次尝试混用了 < 和 >,结果总是不小心创建出重叠的条件。

    兼容类型

    请注意,if_else() 和 case_when() 都要求输出中的类型彼此兼容。如果它们不兼容,你会看到类似这样的错误:

    if_else(TRUE, "a", 1)
    case_when(
    x < -1 ~ TRUE,
    x > 0 ~ now()
    )

    总体来说,兼容的类型相对较少,因为自动将一种向量类型转换为另一种向量类型是错误的常见来源。以下是最重要的兼容情况:

    • 数值向量和逻辑向量是兼容的。
    • 字符串和因子是兼容的,因为你可以把因子看作是取值范围受限的字符串。
    • 日期和日期时间是兼容的,因为你可以把日期看作日期时间的一种特殊情况。
    • NA 在技术上属于逻辑向量,它与任何类型都兼容,因为每种向量都有某种表示缺失值的方式。

    我们并不指望你把这些规则都背下来,但随着时间推移,它们会变得很自然。

  • R语言的交流

    Table of Contents

    https://wp.me/p80aHo-2BN

    引言

    当你制作探索性图表时,甚至在查看之前,你就知道该图表会展示哪些变量。你制作每一张图表都有明确目的,能够快速看一眼,然后继续看下一张图表。在大多数分析过程中,你会生成几十甚至几百张图表,其中大多数会被立刻丢弃。

    现在你已经理解了数据,接下来需要把你的理解传达给他人。你的受众很可能不会和你拥有相同的背景知识,也不会对这些数据投入很深的兴趣。为了帮助他人快速建立起对数据的良好认知,你需要花费相当多的精力,让你的图表尽可能做到自解释。在本章中,你将学习 ggplot2 提供的一些工具来实现这一点。

    本章聚焦于创建优秀图形所需的工具。我们假设你已经知道自己想要什么,只是需要知道如何去实现。因此,我们非常建议将本章与一本优秀的通用可视化书籍搭配阅读。我们尤其喜欢 Albert Cairo 的《The Truthful Art》。这本书并不教授创建可视化的具体操作方法,而是侧重于如何思考,才能制作出有效的图形。

    标签

    将探索性图形转换为说明性图形时,最容易入手的地方就是添加良好的标签。你可以使用 labs() 函数来添加标签。

    library(ggplot2)
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class)) +
    geom_smooth(se = FALSE) +
    labs(
    x = "Engine displacement (L)",
    y = "Highway fuel economy (mpg)",
    color = "Car type",
    title = "Fuel efficiency generally decreases with engine size",
    subtitle = "Two seaters (sports cars) are an exception because of their light weight",
    caption = "Data from fueleconomy.gov"
    )

    图表标题的目的是概括主要发现。应避免那种只描述图表内容的标题,例如“发动机排量与燃油经济性的散点图”。

    如果你需要添加更多文本,还有另外两个有用的标签:subtitle 会在标题下方以较小的字体添加更多细节,而 caption 会在图表右下角添加文本,通常用于说明数据来源。你也可以使用 labs() 来替换坐标轴和图例标题。通常,最好把简短的变量名替换为更详细的描述,并注明单位。

    也可以使用数学公式来代替文本字符串。只需把 "" 换成 quote(),并查看 ?plotmath 中可用的选项。

    library(dplyr)
    df <- tibble(
    x = 1:10,
    y = cumsum(x^2)
    )
    ggplot(df, aes(x, y)) +
    geom_point() +
    labs(
    x = quote(x[i]),
    y = quote(sum(x[i] ^ 2, i == 1, n))
    )

    注释

    除了为图表的主要组成部分添加标签之外,为单个观测值或观测组添加标签也常常很有用。你手头的第一个工具是 geom_text()geom_text() 与 geom_point() 类似,但它多了一个额外的美学映射:label。这使得你可以为图表添加文本标签。

    标签有两个可能的来源。首先,你可能有一个提供标签的 tibble。在下面的图中,我们筛选出每种驱动类型中发动机排量最大的汽车,并将它们的信息保存为一个名为 label_info 的新数据框。

    # 从 mpg 数据集中筛选出每种驱动类型(drv)里排量(displ)最大的汽车
    # 并为这些记录创建更容易理解的驱动类型标签
    label_info <- mpg |>
    # 按驱动类型分组,确保后续操作在每个驱动类型内部独立进行
    group_by(drv) |>
    # 在每个组内按排量从大到小排序,让排量最大的记录排在最前面
    arrange(desc(displ)) |>
    # 取每个组的第一条记录,也就是该驱动类型中排量最大的那辆车
    slice_head(n = 1) |>
    # 新建一个更可读的变量 drive_type,把简写的 drv 转成完整英文描述
    mutate(
    drive_type = case_when(
    drv == "f" ~ "front-wheel drive", # 前轮驱动
    drv == "r" ~ "rear-wheel drive", # 后轮驱动
    drv == "4" ~ "4-wheel drive" # 四轮驱动
    )
    ) |>
    # 只保留后续绘图或标注真正需要的变量
    select(displ, hwy, drv, drive_type)
    # 打印查看筛选和整理后的结果
    label_info
    # 使用 mpg 数据集创建一个散点图,查看发动机排量(displ)与高速油耗(hwy)的关系,
    # 并根据驱动类型(drv)着色
    ggplot(mpg, aes(x = displ, y = hwy, color = drv)) +
    # 绘制散点,alpha = 0.3 让点更透明,避免重叠过于拥挤
    geom_point(alpha = 0.3) +
    # 添加平滑趋势线,不显示标准误差带
    geom_smooth(se = FALSE) +
    # 添加文本标签,使用 label_info 数据框中的内容
    # 这里单独指定数据源,是为了只给少数关键点加注释
    geom_text(
    data = label_info,
    aes(x = displ, y = hwy, label = drive_type),
    # fontface = "bold" 让标签加粗,size 控制字号
    # hjust 和 vjust 用来微调标签相对点的位置
    fontface = "bold", size = 5, hjust = "right", vjust = "bottom"
    ) +
    # 隐藏图例,因为标签已经直接说明了驱动类型
    theme(legend.position = "none")

    你也可以使用相同的思路,通过 ggrepel 包中的 geom_text_repel() 来突出图中的某些点。再注意这里用到的另一个实用技巧:我们额外添加了第二层较大的空心点,以进一步突出这些被标注的点。

    # install.packages("ggrepel")
    # 加载 ggrepel 包,用于绘制自动避让的文本标签,避免标签重叠
    library(ggrepel)
    # 筛选出潜在离群点:
    # 1) 高速油耗 hwy 大于 40 的车辆;
    # 2) 或者同时满足高速油耗大于 20 且发动机排量 displ 大于 5 的车辆
    potential_outliers <- mpg |>
    filter(hwy > 40 | (hwy > 20 & displ > 5))
    # 绘制散点图:横轴为发动机排量(displ),纵轴为高速油耗(hwy)
    ggplot(mpg, aes(x = displ, y = hwy)) +
    # 绘制所有车辆的散点
    geom_point() +
    # 为潜在离群点添加文本标签,标签内容使用车型名称 model
    # geom_text_repel() 会自动调整标签位置,尽量避免与其他点和标签重叠
    geom_text_repel(data = potential_outliers, aes(label = model)) +
    # 用红色点突出显示这些潜在离群点,便于观察
    geom_point(data = potential_outliers, color = "red") +
    # 再次用更醒目的红色空心圆圈标出这些点
    # size 控制圈的大小,shape = "circle open" 表示空心圆
    geom_point(
    data = potential_outliers,
    color = "red", size = 3, shape = "circle open"
    )

    请记住,除了 geom_text() 和 geom_label() 之外,ggplot2 中还有许多其他几何对象可用于帮助你为图表添加注释。下面是一些想法:

    • 使用 geom_hline() 和 geom_vline() 来添加参考线。我们通常会把它们画得较粗(linewidth = 2)并设为白色(color = "white"),同时将它们放在主要数据层的下方。这样可以让它们清晰可见,同时又不会把注意力从数据本身上移开。
    • 使用 geom_rect() 可以在感兴趣的点周围绘制一个矩形。矩形的边界由美学映射 xminxmaxyminymax 来定义。或者,你也可以查看 ggforce 包,特别是 geom_mark_hull(),它可以让你用凸包来标注点的子集。
    • 使用带有 arrow 参数的 geom_segment() 来通过箭头突出某个点。使用美学映射 x 和 y 来定义起始位置,使用 xend 和 yend 来定义结束位置。

    另一个用于向图表添加注释的实用函数是 annotate()。一般来说,geom 通常适合突出显示数据的一个子集,而 annotate() 则适合向图表中添加一个或少数几个注释元素。

    为了演示如何使用 annotate(),我们先创建一些要添加到图中的文本。由于这段文本有点长,我们会使用 stringr::str_wrap(),根据你希望每行显示的字符数自动给它添加换行:

    library(stringr)
    # 将一段说明文字按指定宽度自动换行,便于在图表中显示
    trend_text <- "Larger engine sizes tend to have lower fuel economy." |>
    str_wrap(width = 30)
    # 输出处理后的文本
    trend_text

    然后,我们添加两层注释:一层使用标签几何对象,另一层使用线段几何对象。两者中的 x 和 y 美学映射定义注释应从何处开始,而线段注释中的 xend 和 yend 美学映射定义线段的结束位置。另请注意,这条线段被设置成了箭头样式。

    # 绘制 mpg 数据集中发动机排量(displ)与高速油耗(hwy)的散点图
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point() +
    # 添加标签注释:在图中指定位置放置一段说明文字,用来解释总体趋势
    annotate(
    geom = "label", # 使用标签几何对象,显示带边框的文本框
    x = 3.5, # 标签的横坐标位置
    y = 38, # 标签的纵坐标位置
    label = trend_text, # 标签内容,使用前面整理好的趋势说明文字
    hjust = "left", # 文本水平左对齐,使排版更自然
    color = "red" # 文本颜色设置为红色,便于突出显示
    ) +
    # 添加线段注释:用一条带箭头的线段,把文字和图中的目标位置联系起来
    annotate(
    geom = "segment", # 使用线段几何对象,绘制一条连接线
    x = 3, # 线段起点的横坐标
    y = 35, # 线段起点的纵坐标
    xend = 5, # 线段终点的横坐标
    yend = 25, # 线段终点的纵坐标
    color = "red" # 线段颜色设置为红色,与标签保持一致
    arrow = arrow(type = "closed") # 在线段末端添加闭合箭头,增强指向性
    )

    注释是传达可视化主要结论和有趣特征的强大工具。唯一的限制只是你的想象力(以及你为了让注释排版更美观而调整位置时的耐心)!

    比例尺

    改进图表以便更好地传达信息的第三种方法是调整比例尺。比例尺控制美学映射在视觉上的呈现方式。

    默认比例尺

    通常情况下,ggplot2 会自动为你添加比例尺。例如,当你输入:

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class))

    ggplot2 会在后台自动添加默认比例尺:

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class)) +
    scale_x_continuous() +
    scale_y_continuous() +
    scale_color_discrete()

    请注意尺度(比例尺)的命名规则:scale_ 后接美学属性名称,然后是 _,最后是尺度名称。默认尺度的命名取决于它们所对应的变量类型:continuous、discrete、datetime 或 date。scale_x_continuous() 会把 displ 中的数值映射到 x 轴上的连续数轴上,scale_color_discrete() 会为每一类汽车选择颜色,等等。下面你将学习很多非默认尺度。

    默认比例尺经过精心选择,能够很好地适用于各种输入。不过,出于两个原因,你可能希望覆盖这些默认设置:

    • 你可能会想调整默认尺度的一些参数。这样你就可以更改坐标轴上的刻度,或者图例中的键标签。
    • 你可能想完全替换该尺度,并使用一种完全不同的算法。通常,由于你对数据了解更多,因此可以做得比默认设置更好。

    坐标轴刻度和图例键

    统称起来,坐标轴和图例被称为导引。坐标轴用于 x 和 y 美学属性;图例用于其他所有美学属性。

    有两个主要参数会影响坐标轴上刻度以及图例中的键的外观:breaks 和 labelsbreaks 控制刻度的位置,或与键相关联的值。labels 控制与每个刻度/键对应的文本标签。breaks 最常见的用途是覆盖默认选择:

    ggplot(mpg, aes(x = displ, y = hwy, color = drv)) +
    geom_point() +
    scale_y_continuous(breaks = seq(15, 40, by = 5))

    你可以以同样的方式使用标签(一个长度与刻度相同的字符向量),但也可以将其设置为 NULL 以完全隐藏标签。这对于地图,或者用于发布那些你不能分享绝对数值的图表时会很有用。你也可以使用刻度和标签来控制图例的外观。对于分类变量的离散尺度,标签可以是一个命名列表,其中包含现有的级别名称及其对应的期望标签。

    # 绘制 mpg 数据集中排量与高速油耗的散点图,并按驱动方式着色
    ggplot(mpg, aes(x = displ, y = hwy, color = drv)) +
    # 添加散点层
    geom_point() +
    # 隐藏 x 轴刻度标签
    scale_x_continuous(labels = NULL) +
    # 隐藏 y 轴刻度标签
    scale_y_continuous(labels = NULL) +
    # 自定义图例标签:4=四驱,f=前驱,r=后驱
    scale_color_discrete(labels = c("4" = "4-wheel", "f" = "front", "r" = "rear"))

    labels 参数配合 scales 包中的标注函数,也可用于将数字格式化为货币、百分比等。左侧图展示了使用 label_dollar() 的默认标注方式,它会添加美元符号以及千位分隔符逗号。右侧图则通过将美元值除以 1,000,并添加后缀 “K”(表示“千”),再加上自定义的刻度,进一步实现了个性化设置。请注意,breaks 仍然采用数据的原始尺度。

    library(scales)
    # Left
    ggplot(diamonds, aes(x = price, y = cut)) +
    geom_boxplot(alpha = 0.05) +
    scale_x_continuous(labels = label_dollar())
    # Right
    ggplot(diamonds, aes(x = price, y = cut)) +
    geom_boxplot(alpha = 0.05) +
    scale_x_continuous(
    labels = label_dollar(scale = 1/1000, suffix = "K"),
    breaks = seq(1000, 19000, by = 6000)
    )

    另一个很实用的标签函数是 label_percent()

    # 使用 diamonds 数据集,按 cut 作为 x 轴,并按 clarity 填充颜色
    ggplot(diamonds, aes(x = cut, fill = clarity)) +
    # 绘制堆叠比例柱状图,按总量归一化为 100%
    geom_bar(position = "fill") +
    # 设置 y 轴名称为 Percentage,并将刻度标签显示为百分比格式
    scale_y_continuous(name = "Percentage", labels = label_percent())

    breaks 的另一个用途是:当你只有相对较少的数据点,并希望准确突出观测值出现的位置时。比如,下面这张图展示了每位美国总统开始和结束任期的时间。

    # 对 presidential 数据进行处理,生成按顺序递增的 id
    presidential |>
    mutate(id = 33 + row_number()) |>
    # 绘制开始时间与 id 的点线图
    ggplot(aes(x = start, y = id)) +
    # 添加每位总统任期开始时间的点
    geom_point() +
    # 添加从 start 到 end 的线段,表示任期区间
    geom_segment(aes(xend = end, yend = id)) +
    # 设置 x 轴为日期轴,去掉名称,并以总统任期开始日期作为刻度,标签显示为两位年份
    scale_x_date(name = NULL, breaks = presidential$start, date_labels = "'%y")

    请注意,对于 breaks 参数,我们使用 presidential$start 将 start 变量提取为一个向量,因为这个参数不能进行美学映射。另请注意,日期和日期时间尺度中 breaks 与 labels 的指定方式略有不同:

    • breaks 不能像 aes() 里的变量那样“按数据行自动映射”,所以这里要直接传入一个向量presidential$start
    • presidential$start 的意思是把数据框里 start 这一列取出来,作为 breaks 的具体位置。
    • date_labels 用来控制刻度标签怎么显示,它后面接的是一种格式字符串,比如年份、月份、日期的显示格式。
    • date_breaks 用来控制刻度间隔,比如每 2 天、每 1 个月一个刻度,所以它接收像 "2 days""1 month" 这样的字符串。

    图例布局

    你最常会使用 breaks 和 labels 来调整坐标轴。虽然它们也可以用于图例,但还有一些其他技巧你更可能会用到。

    要控制图例的整体位置,你需要使用 theme() 设置。我们会在本章结尾回到主题(theme)部分,但简单来说,它们用于控制图形中非数据的部分。theme 设置中的 legend.position 控制图例的绘制位置:

    
    # 创建一个基础散点图:x 轴为发动机排量 displ,y 轴为高速油耗 hwy
    # 并根据 class 给点上色
    base <- ggplot(mpg, aes(x = displ, y = hwy)) +
      geom_point(aes(color = class))
    
    # 将图例放在右侧(默认位置)
    base + theme(legend.position = "right") # the default
    
    # 将图例放在左侧
    base + theme(legend.position = "left")
    
    # 将图例放在上方,并将颜色图例分成 3 行显示
    base + 
      theme(legend.position = "top") +
      guides(color = guide_legend(nrow = 3))
    
    # 将图例放在下方,并将颜色图例分成 3 行显示
    base + 
      theme(legend.position = "bottom") +
      guides(color = guide_legend(nrow = 3))

    如果你的绘图区域偏矮宽,就将图例放置在顶部或底部;若绘图区域偏高窄,则把图例放在左侧或右侧。你也可以通过设置 legend.position = “none” 来彻底隐藏图例。

    若要单独控制各个图例的展示效果,可搭配guides()函数与guide_legend()函数或guide_colorbar()函数使用。下方示例展示了两项重要设置:通过nrow参数设置图例所占行数,以及修改某项美学参数来放大数据点。倘若你在绘图时设置了较低透明度以呈现大量数据点,该用法会格外实用。

    # 创建一个散点图:x 轴为发动机排量 displ,y 轴为高速油耗 hwy
    # 点的颜色按 class 分类
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class)) +
    # 添加平滑曲线,不显示置信区间
    geom_smooth(se = FALSE) +
    # 将图例放在底部
    theme(legend.position = "bottom") +
    # 调整颜色图例:
    # - 分成 2 行显示
    # - override.aes 用于修改图例中的点样式
    # 这里把图例点的大小设为 4
    guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))

    注意,guides() 中参数的名称与该美学的名称一致,就像在 labs() 中一样。

    替换比例尺

    与其只是稍微调整细节,不如直接完全替换比例尺。你最可能想要替换的比例尺主要有两类:连续位置比例尺和颜色比例尺。幸运的是,同样的原则也适用于其他所有美学,因此一旦你掌握了位置和颜色,就能很快学会其他比例尺的替换。

    绘制变量的变换通常非常有用。例如,如果对克拉和价格进行对数变换,就更容易看清它们之间的精确关系:

    # Left
    ggplot(diamonds, aes(x = carat, y = price)) +
    geom_bin2d()
    #> `stat_bin2d()` using `bins = 30`. Pick better value `binwidth`.
    # Right
    ggplot(diamonds, aes(x = log10(carat), y = log10(price))) +
    geom_bin2d()
    #> `stat_bin2d()` using `bins = 30`. Pick better value `binwidth`.

    然而,这种转换的缺点在于,轴线现在是用转换后的值来标记的,这使得解读图表变得困难。我们可以在美学映射中不进行转换,而是用刻度来进行转换。这在视觉上是一样的,只是轴线是按照原始数据刻度来标记的。

    ggplot(diamonds, aes(x = carat, y = price)) +
    geom_bin2d() +
    scale_x_log10() +
    scale_y_log10()

    另一个经常需要自定义的比例尺是颜色。默认的分类比例尺会选择在色轮上均匀分布的颜色。一个有用的替代方案是 ColorBrewer 比例尺,它们经过人工调校,对常见色盲人群更友好。下面这两幅图看起来相似,但红色和绿色的色调差异足够明显,因此右图中的点即使对红绿色盲的人也能区分出来。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv))
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv)) +
    scale_color_brewer(palette = "Set1")

    不要忘记一些更简单的可访问性改进技巧。如果颜色数量不多,你可以添加一个冗余的形状映射。这也有助于确保你的图在黑白情况下仍然易于解读。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv, shape = drv)) +
    scale_color_brewer(palette = "Set1")

    ColorBrewer 比例尺的文档可在 https://colorbrewer2.org/ 在线查看,并可通过 Erich Neuwirth 开发的 RColorBrewer 包在 R 中使用。图 11.1 展示了所有调色板的完整列表。如果你的分类值是有序的,或者有一个“中间值”,那么顺序型(上方)和发散型(下方)调色板尤其有用。这种情况通常会在你使用 cut() 将连续变量转换为分类变量时出现。

    当你已经预先定义了数值与颜色之间的映射时,可以使用 scale_color_manual()。例如,如果我们把总统所属政党映射为颜色,就希望采用标准对应关系:共和党用红色,民主党用蓝色。分配这些颜色的一种方法是使用十六进制颜色代码:

    presidential |>
    mutate(id = 33 + row_number()) |>
    ggplot(aes(x = start, y = id, color = party)) +
    geom_point() +
    geom_segment(aes(xend = end, yend = id)) +
    scale_color_manual(values = c(Republican = "#E81B23", Democratic = "#00AEF3"))

    对于连续颜色,可以使用内置的 scale_color_gradient() 或 scale_fill_gradient()。如果你使用的是发散型比例尺,可以使用 scale_color_gradient2()。这样可以让你为正值和负值指定不同的颜色。例如,当你想区分均值以上和均值以下的点时,这也很有用。

    另一种选择是使用 viridis 颜色比例尺。设计者 Nathaniel Smith 和 Stéfan van der Walt 精心设计了连续色彩方案,使其不仅能被各种色盲人群感知,而且在颜色和黑白显示中都具有知觉上的一致性。这些比例尺在 ggplot2 中提供连续型(c)、离散型(d)和分箱型(b)调色板。

    # 创建包含 10000 个随机正态分布样本的二维数据框
    df <- tibble(
    x = rnorm(10000),
    y = rnorm(10000)
    )
    # 默认连续色标的六边形分箱图
    ggplot(df, aes(x, y)) +
    geom_hex() +
    coord_fixed() +
    labs(title = "Default, continuous", x = NULL, y = NULL)
    # 使用 Viridis 连续色标的六边形分箱图
    ggplot(df, aes(x, y)) +
    geom_hex() +
    coord_fixed() +
    scale_fill_viridis_c() +
    labs(title = "Viridis, continuous", x = NULL, y = NULL)
    # 使用 Viridis 分箱色标的六边形分箱图
    ggplot(df, aes(x, y)) +
    geom_hex() +
    coord_fixed() +
    scale_fill_viridis_b() +
    labs(title = "Viridis, binned", x = NULL, y = NULL)

    请注意,所有颜色比例尺都有两种形式:分别对应 color 和 fill 美学映射的 scale_color_*() 和 scale_fill_*()(其中 color 比例尺同时支持英式和美式拼写)。

    缩放

    控制图形范围有三种方式:

    1. 调整要绘制的数据。
    2. 在每个比例尺中设置范围。
    3. 在 coord_cartesian() 中设置 xlim 和 ylim

    我们将通过一系列图形来演示这些选项。左侧的图展示了发动机排量与燃油效率之间的关系,并按驱动类型着色。右侧的图展示了相同的变量,但只绘制了数据的一个子集。对数据进行子集化不仅影响了 x 和 y 比例尺,也影响了平滑曲线。

    # Left
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv)) +
    geom_smooth()
    # Right
    mpg |>
    filter(displ >= 5 & displ <= 6 & hwy >= 10 & hwy <= 25) |>
    ggplot(aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv)) +
    geom_smooth()

    让我们把这些与下面的两个图进行比较:左图是在各个比例尺中设置范围,右图是在 coord_cartesian() 中设置范围。我们可以看到,缩小范围等同于对数据进行子集筛选。因此,如果要放大图中的某个区域,通常最好使用 coord_cartesian()

    # Left
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv)) +
    geom_smooth() +
    scale_x_continuous(limits = c(5, 6)) +
    scale_y_continuous(limits = c(10, 25))
    # Right
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = drv)) +
    geom_smooth() +
    coord_cartesian(xlim = c(5, 6), ylim = c(10, 25))

    另一方面,如果你想扩展范围,例如使不同图形之间的比例尺一致,那么通常在单个比例尺上设置范围会更有用。比如,如果我们提取两类汽车并分别绘图,由于三个比例尺(x 轴、y 轴和颜色美学映射)都具有不同的范围,因此很难比较这些图形。

    suv <- mpg |> filter(class == "suv")
    compact <- mpg |> filter(class == "compact")
    # Left
    ggplot(suv, aes(x = displ, y = hwy, color = drv)) +
    geom_point()
    # Right
    ggplot(compact, aes(x = displ, y = hwy, color = drv)) +
    geom_point()

    克服这一问题的一种方法是在多个图形之间共享比例尺,并使用完整数据的范围来训练这些比例尺。

    x_scale <- scale_x_continuous(limits = range(mpg$displ))
    y_scale <- scale_y_continuous(limits = range(mpg$hwy))
    col_scale <- scale_color_discrete(limits = unique(mpg$drv))
    # Left
    ggplot(suv, aes(x = displ, y = hwy, color = drv)) +
    geom_point() +
    x_scale +
    y_scale +
    col_scale
    # Right
    ggplot(compact, aes(x = displ, y = hwy, color = drv)) +
    geom_point() +
    x_scale +
    y_scale +
    col_scale

    在这个特定情况下,你本可以直接使用分面,但这种技术在更一般的场景下也很有用,比如当你想把图形分布到报告的多页中时。

    主题

    最后,你可以使用主题来自定义图形中的非数据元素:

    # 使用 mpg 数据集创建散点图,x 轴为发动机排量 displ,y 轴为高速油耗 hwy
    ggplot(mpg, aes(x = displ, y = hwy)) +
    # 绘制散点图,并按车辆类别 class 着色
    geom_point(aes(color = class)) +
    # 添加平滑曲线,不显示标准误差带
    geom_smooth(se = FALSE) +
    # 使用黑白主题
    theme_bw()

    ggplot2 包含图 11.2 所示的八种主题,默认主题是 theme_gray()

    也可以控制每个主题的单独组成部分,比如 y 轴所用字体的大小和颜色。我们已经看到,legend.position 控制图例的绘制位置。使用 theme() 还可以自定义图例的许多其他方面。例如,在下面的图中,我们改变了图例的方向,并为其添加了黑色边框。注意,图例框和图形标题元素的自定义是通过 element_*() 函数完成的。这些函数用于指定非数据组件的样式,例如,标题文本的粗体是在 element_text() 的 face 参数中设置的,而图例边框的颜色则在 element_rect() 的 color 参数中定义。控制标题和副标题位置的主题元素分别是 plot.title.position 和 plot.caption.position。在下面的图中,这些值被设置为 "plot",表示这些元素与整个绘图区对齐,而不是与图面板对齐(默认设置)。另外还使用了几个有用的 theme() 组件来调整标题和副标题文本的排列和格式。

    # 使用 mpg 数据集绘制散点图
    # x 轴表示发动机排量 displ,y 轴表示高速油耗 hwy
    # 并根据驱动方式 drv 对点进行着色
    ggplot(mpg, aes(x = displ, y = hwy, color = drv)) +
    # 绘制每辆车对应的散点
    geom_point() +
    # 添加标题和数据来源说明
    labs(
    title = "发动机越大,燃油经济性通常越低",
    caption = "数据来源:https://fueleconomy.gov。"
    ) +
    # 使用 theme() 自定义图形外观:
    # - 调整图例位置
    # - 设置图例横向排列
    # - 给图例外框加黑色边框
    # - 将标题设为加粗
    # - 让标题和注释相对于整个图形区域对齐
    # - 将注释左对齐
    theme(
    legend.position = c(0.6, 0.7),
    legend.direction = "horizontal",
    legend.box.background = element_rect(color = "black"),
    plot.title = element_text(face = "bold"),
    plot.title.position = "plot",
    plot.caption.position = "plot",
    plot.caption = element_text(hjust = 0)
    )

    要查看所有 theme() 组件的概览,请参阅 ?theme 的帮助文档。ggplot2 书籍也是了解主题设置完整细节的绝佳去处。

    到目前为止,我们讨论的是如何创建和修改单个图形。那么,如果你有多个图形,并且想以某种方式将它们排列在一起,该怎么办呢?patchwork 包允许你把多个独立的图形组合到同一张图中。我们在本章前面已经加载过这个包。

    要将两个图并排放置,你只需把它们相加即可。注意,你首先需要创建这些图并将它们保存为对象(在下面的示例中,它们被命名为 p1 和 p2)。然后,使用 + 将它们并排放置。

    install.packages("patchwork")
    library(patchwork)
    p1 <- ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point() +
    labs(title = "Plot 1")
    p2 <- ggplot(mpg, aes(x = drv, y = hwy)) +
    geom_boxplot() +
    labs(title = "Plot 2")
    p1 + p2

    需要注意的是,在上面的代码块中,我们并没有使用 patchwork 包中的新函数。相反,这个包为 + 运算符添加了新的功能。

    你也可以使用 patchwork 创建复杂的图布局。下面,| 将 p1 和 p3 并排放置,/ 则将 p2 移到下一行。

    p3 <- ggplot(mpg, aes(x = cty, y = hwy)) +
    geom_point() +
    labs(title = "Plot 3")
    (p1 | p3) / p2

    此外,patchwork 还允许你将多个图中的图例收集到一个公共图例中,自定义图例的位置以及各个图的尺寸,并为你的图添加共同的标题、副标题、说明等。下面我们创建 5 个图。我们已经关闭了箱线图和散点图的图例,并使用 & theme(legend.position = "top") 将密度图的图例收集到图的顶部。请注意这里使用的是 & 运算符,而不是通常的 +。这是因为我们修改的是 patchwork 图的主题,而不是单个 ggplot。图例被放置在顶部的 guide_area() 内。最后,我们还自定义了 patchwork 各个组成部分的高度——guide 的高度为 1,箱线图为 3,密度图为 2,分面散点图为 4。Patchwork 会根据这个比例分配你为图预留的区域,并相应地放置各个组件。

    p1 <- ggplot(mpg, aes(x = drv, y = cty, color = drv)) +
    geom_boxplot(show.legend = FALSE) +
    labs(title = "Plot 1")
    p2 <- ggplot(mpg, aes(x = drv, y = hwy, color = drv)) +
    geom_boxplot(show.legend = FALSE) +
    labs(title = "Plot 2")
    p3 <- ggplot(mpg, aes(x = cty, color = drv, fill = drv)) +
    geom_density(alpha = 0.5) +
    labs(title = "Plot 3")
    p4 <- ggplot(mpg, aes(x = hwy, color = drv, fill = drv)) +
    geom_density(alpha = 0.5) +
    labs(title = "Plot 4")
    p5 <- ggplot(mpg, aes(x = cty, y = hwy, color = drv)) +
    geom_point(show.legend = FALSE) +
    facet_wrap(~drv) +
    labs(title = "Plot 5")
    (guide_area() / (p1 + p2) / (p3 + p4) / p5) +
    plot_annotation(
    title = "City and highway mileage for cars with different drive trains",
    caption = "Source: https://fueleconomy.gov."
    ) +
    plot_layout(
    guides = "collect",
    heights = c(1, 3, 2, 4)
    ) &
    theme(legend.position = "top")

    如果你想进一步了解如何使用 patchwork 组合和排布多个图,我们建议查看该包官网上的指南:https://patchwork.data-imaginist.com。

  • R语言的探索性数据分析

    https://wp.me/p80aHo-2zN

    Table of Contents

    引言

    本章将向你展示如何利用可视化和转换系统性地探索数据,统计学家称之为探索性数据分析(简称EDA)。EDA是一个迭代循环。你:

    1. 对你的数据提出疑问。
    2. 通过可视化、转换和建模你的数据来寻找答案。
    3. 用你学到的东西来完善你的问题和/或生成新的问题。

    EDA不是一个有严格规则的正式流程。EDA最重要的是心态。在EDA的初期阶段,你应该自由地探索所有想到的想法。其中一些想法会实现,有些则会是死胡同。随着探索的深入,你会聚焦于一些特别有成效的见解,最终将它们写出来并传达给他人。

    EDA是任何数据分析的重要组成部分,即使主要研究问题是直接递给你的,因为你总需要调查数据的质量。数据清理只是EDA的一个应用:你会询问你的数据是否符合你的期望。要进行数据清理,你需要部署EDA的所有工具:可视化、转换和建模。

    问题

    “没有常规的统计问题,只有值得质疑的统计程序。”——大卫·考克斯爵士

    “对正确问题的近似答案(往往模糊)远比对错误问题的准确答案要好得多,错误问题总能被精确化。”——约翰·图基

    你在EDA期间的目标是对你的数据有更深入的理解。最简单的方法是用问题作为引导调查的工具。当你提问时,问题会让你关注数据集的某个特定部分,帮助你决定要制作哪些图表、模型或转换。

    EDA本质上是一个创造性的过程。像大多数创意过程一样,提出高质量问题的关键是产生大量问题。在分析开始时很难提出揭示性的问题,因为你不知道从数据集中能获得哪些洞见。另一方面,每问一个新问题,都会让你接触到数据的新方面,增加你发现新事物的机会。如果你根据发现为每个问题提出新问题,就能快速深入挖掘数据中最有趣的部分,并设计出一系列发人深省的问题。

    没有规定你应该问哪些问题来指导你的研究。然而,有两种类型的问题总是有助于你在数据中发现真相。你可以大致地这样表述这些问题:

    1. 我的变量中会发生什么样的变化?
    2. 我的变量之间会发生什么样的协变?

    本章的其余部分将探讨这两个问题。我们将解释什么是变异和协变,并展示几种回答每个问题的方法。

    变异

    变异是指变量值在不同测量过程中变化的趋势。现实生活中很容易看到差异;如果你测量任意一个连续变量两次,结果会不同。即使测量的是恒定的量,比如光速,这一点也成立。你的每次测量都会包含一个小幅误差,且误差会因测量而异。变量也会变化,比如测量不同受试者(例如不同人的眼睛颜色)或不同时间(例如电子在不同时刻的能级)。每个变量都有其独特的变异模式,这可以揭示它在同一观测值以及跨观测值间的变化有趣信息。理解这种模式的最佳方法是可视化变量值的分布。

    我们将从可视化钻石数据集中~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 的钻石,保存为 smaller
    smaller <- 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 毫米的测量值不合理:这些钻石都超过一英寸长,但价格却并没有高到几十万美元!

    最好分别在有异常值和没有异常值的情况下重复你的分析。如果异常值对结果影响很小,而且你也弄不清它们为什么会出现,那么可以合理地省略它们,然后继续进行。然而,如果它们对你的结果有显著影响,你就不应该在没有充分理由的情况下删掉它们。你需要弄清楚是什么导致了这些异常值(例如,数据录入错误),并在你的书面报告中说明你已经将它们移除。

    异常值

    如果你在数据集中遇到了不寻常的值,并且只是想继续进行其余的分析,你有两个选择。

    1. 删除包含这些异常值的整行:
    # 筛选 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 轴为高速公路油耗 hwy
    ggplot(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

  • R语言的层

    https://wp.me/p80aHo-2wM

    前言

    在前面,你学到的不仅仅是如何制作散点图、柱状图和箱线图。你学到的是一个基础,这个基础可以让你用 ggplot2 制作任何类型的图形。

    在本章中,你将继续扩展这一基础,学习分层的图形语法。我们将先深入了解美学映射、几何对象和分面。然后,你将学习 ggplot2 在绘图时在底层自动完成的统计变换。这些变换用于计算要绘制的新数值,例如柱状图中柱子的高度或箱线图中的中位数。你还将学习位置调整,它会改变几何对象在图中的显示方式。最后,我们将简要介绍坐标系。

    我们不会逐一介绍这些层中的每一个函数和选项,但会带你了解 ggplot2 提供的最重要、最常用的功能,并向你介绍一些可以扩展 ggplot2 的包。

    美学映射

    请记住,ggplot2 包自带的 mpg 数据框包含 38 款汽车模型的 234 个观测值。

    mpg 中的变量包括:

    displ:汽车发动机排量,单位为升。数值型变量。

    hwy:汽车在高速公路上的燃油效率,单位为每加仑英里数(mpg)。在相同距离下,燃油效率低的汽车比燃油效率高的汽车消耗更多燃料。数值型变量。

    class:汽车类型。分类变量。

    让我们先可视化不同汽车类型中 displ 和 hwy 之间的关系。我们可以用散点图来实现:将数值型变量映射到 x 和 y 美学,将分类变量映射到颜色或形状等美学属性。

    ggplot(mpg, aes(x = displ, y = hwy, color = class)) +
    geom_point()
    ggplot(mpg, aes(x = displ, y = hwy, shape = class)) +
    geom_point()

    当将 class 映射到形状时,我们会得到两个警告:

    1:形状调色板最多只能处理 6 个离散值,因为超过 6 个就会变得难以区分;而你有 7 个。若必须使用这些形状,请考虑手动指定。

    2:移除了 62 行包含缺失值的记录(geom_point())。

    由于 ggplot2 一次只会使用 6 种形状,默认情况下,当你使用形状美学时,额外的组将不会被绘制。第二个警告与此有关——数据集中有 62 辆 SUV,但它们没有被绘制出来。

    同样,我们也可以把 class 映射到 size 或 alpha 美学上,分别控制点的大小和透明度。

    ggplot(mpg, aes(x = displ, y = hwy, size = class)) +
    geom_point()
    #> Warning: Using size for a discrete variable is not advised.
    ggplot(mpg, aes(x = displ, y = hwy, alpha = class)) +
    geom_point()
    #> Warning: Using alpha for a discrete variable is not advised.

    这两种做法也都会产生警告:

    不建议将 alpha 用于离散变量。

    将无序的离散(分类)变量(class)映射到有序的美学属性(size 或 alpha)通常不是一个好主意,因为这会暗示一种实际上并不存在的排序关系。

    一旦你完成了美学映射,ggplot2 会自动处理其余部分。它会为该美学选择一个合适的尺度,并构建一个图例来解释水平与数值之间的映射。对于 x 和 y 美学,ggplot2 不会创建图例,但会创建带有刻度和标签的坐标轴。坐标轴提供的信息与图例相同;它解释了位置与数值之间的映射关系。

    你也可以直接在 geom 函数的参数中(也就是 aes() 外部)手动设置几何对象的视觉属性,而不是依赖变量映射来决定外观。例如,我们可以把图中的所有点都设为蓝色:

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(color = "blue")

    这里,颜色并不传达关于某个变量的信息,而只是改变图形的外观。你需要为该美学属性选择一个合适的值:

    • 颜色的名称,写成字符串,例如 color = "blue"
    • 点的大小,单位为毫米,例如 size = 1
    • 点的形状,写成数字,例如 shape = 1

    到目前为止,我们已经讨论了在使用点几何对象绘制散点图时,可以映射或设置的各种美学属性。你可以在美学规范简介(vignette)中了解更多关于所有可能的美学映射:https://ggplot2.tidyverse.org/articles/ggplot2-specs.html。

    你在图中可以使用的具体美学属性取决于你用来表示数据的几何对象。在下一节中,我们将更深入地介绍几何对象。

    练习

    • 创建一个 hwy 对 displ 的散点图,点为粉色填充三角形。
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(color = "pink", shape = 17)
    • 为什么下面的代码没有生成一个带有蓝色点的图?
    ggplot(mpg) +
    geom_point(aes(x = displ, y = hwy, color = "blue"))
    # 正确
    ggplot(mpg) +
    geom_point(aes(x = displ, y = hwy), color = "blue")
    • stroke 美学属性有什么作用?它适用于哪些形状?(提示:使用 ?geom_point

    stroke 用来控制点的边框线宽

    它主要适用于 21–25 号形状,也就是既有填充色又有边框的点形状。
    例如:

    • shape = 21 圆
    • shape = 22 方形
    • shape = 23 菱形
    • shape = 24 上三角
    • shape = 25 下三角

    对这些形状,stroke 会影响边框粗细;对其他大多数形状,作用不明显或不起作用。

    几何对象

    这两幅图有什么相似之处?

    这两个图都包含相同的 x 变量、相同的 y 变量,而且都描述了同一组数据。但这些图并不完全相同。每个图都使用不同的几何对象 geom 来表示数据。左边的图使用 point geom,右边的图使用 smooth geom,即一条拟合数据的平滑线。

    要更改图中的 geom,只需修改你添加到 ggplot() 中的 geom 函数。例如,要生成上面的图,你可以使用以下代码:

    # Left
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point()
    # Right
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_smooth()
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

    ggplot2 中的每个 geom 函数都接受一个映射参数,这个参数可以在 geom 层中局部定义,也可以在 ggplot() 层中全局定义。不过,并不是每种美学都能适用于每个 geom。你可以设置点的形状,但不能设置线的“形状”。如果你尝试这样做,ggplot2 会静默忽略该美学映射。另一方面,你可以设置线的线型。geom_smooth() 会根据你映射到 linetype 的变量的每个唯一值,绘制一条具有不同线型的不同线。

    ggplot(mpg, aes(x = displ, y = hwy, shape = drv)) +
    geom_smooth()
    ggplot(mpg, aes(x = displ, y = hwy, linetype = drv)) +
    geom_smooth()

    这里,geom_smooth() 根据汽车的 drv 值将汽车分成三条线,drv 值描述的是汽车的传动系统。一条线描述所有 4 值的点,一条线描述所有 f 值的点,还有一条线描述所有 r 值的点。这里,4 代表四轮驱动,f 代表前轮驱动,r 代表后轮驱动。

    如果这听起来有些奇怪,我们可以通过将这些线叠加在原始数据之上,并按 drv 对所有内容着色,来让它更清楚一些。

    ggplot(mpg, aes(x = displ, y = hwy, color = drv)) +
    geom_point() +
    geom_smooth(aes(linetype = drv))

    注意,这张图在同一个图形中包含了两个 geom。

    许多 geom,比如 geom_smooth(),使用单个几何对象来显示多行数据。对于这些 geom,你可以将 group 美学映射设置为一个分类变量,以绘制多个对象。ggplot2 会为分组变量的每个唯一值绘制一个单独的对象。实际上,只要你将某个美学映射到离散变量,ggplot2 就会自动对这些 geom 的数据进行分组(如线型示例所示)。依赖这一特性很方便,因为单独使用 group 美学映射并不会为 geom 添加图例或区分特征。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_smooth()
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_smooth(aes(group = drv))
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_smooth(aes(color = drv), show.legend = FALSE)

    如果你在某个 geom 函数中放置映射,ggplot2 会将它们视为该图层的局部映射。它会使用这些映射来扩展或覆盖全局映射,但仅适用于该图层。这使得在不同图层中显示不同的美学属性成为可能。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point()
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class))
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = class)) +
    geom_smooth()

    你可以使用相同的思路为每一层指定不同的数据。这里,我们使用红色点和空心圆圈来突出显示双座汽车。geom_point() 中的局部 data 参数只会覆盖该图层中的全局 data 参数。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point() +
    geom_point(
    data = dplyr::filter(mpg, class == "2seater"),
    color = "red"
    ) +
    geom_point(
    data = dplyr::filter(mpg, class == "2seater"),
    shape = 1, size = 3, color = "red"
    )

    Geom 是 ggplot2 的基本构建块。通过更改图中的 geom,你可以完全改变图形的外观,而不同的 geom 还能揭示数据的不同特征。例如,下方的直方图和密度图显示,高速公路里程的分布是双峰且右偏的,而箱线图则揭示了两个潜在的离群值。

    ggplot(mpg, aes(x = hwy)) +
    geom_histogram(binwidth = 2)
    ggplot(mpg, aes(x = hwy)) +
    geom_density()
    ggplot(mpg, aes(x = hwy)) +
    geom_boxplot()

    ggplot2 提供了 40 多种 geom,但这些并不能涵盖所有可能绘制的图形。如果你需要不同的 geom,我们建议先查看扩展包,看看是否有人已经实现了它(可参见 https://exts.ggplot2.tidyverse.org/gallery/ 获取一些示例)。例如,ggridges 包(https://wilkelab.org/ggridges)可用于制作山脊图,这类图在可视化数值变量在分类变量不同水平下的分布密度时很有用。在下面的图中,我们不仅使用了一个新的 geom(geom_density_ridges()),还将同一个变量映射到了多个美学属性(将 drv 映射到 yfill 和 color),并且设置了一个美学属性(alpha = 0.5)来使密度曲线透明。

    library(ggridges)
    ggplot(mpg, aes(x = hwy, y = drv, fill = drv, color = drv)) +
    geom_density_ridges(alpha = 0.5, show.legend = FALSE)

    了解 ggplot2 提供的所有 geom 以及包中所有函数的全面概览的最佳地点,是参考页面:https://ggplot2.tidyverse.org/reference。若想了解某个具体的 geom,可以使用帮助文档(例如,?geom_smooth)。

    练习

    • geom_smooth() 中的 se 参数有什么作用?

    se 用来控制是否显示平滑曲线的标准误差带

    1. se = TRUE:显示置信带(默认)
    2. se = FALSE:不显示置信带,只画平滑线

    分面

    在第 1 章中,你学习了使用 facet_wrap() 进行分面,它会根据一个分类变量把图形拆分成多个子图,每个子图显示数据的一个子集。

    # 加载 ggplot2
    ggplot(mpg, aes(x = displ, y = hwy)) +
    # 绘制散点图
    geom_point() +
    # 按 cyl 变量分面显示
    facet_wrap(~cyl)

    要根据两个变量的组合对图形进行分面,可以将 facet_wrap() 改为 facet_grid()facet_grid() 的第一个参数同样是一个公式,但这次是双边公式:rows ~ cols

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point() +
    facet_grid(drv ~ cyl)

    默认情况下,每个分面共享相同的 x 轴和 y 轴刻度与范围。当你想在不同分面之间比较数据时,这很有用,但如果你想更好地观察每个分面内部的关系,它可能会受到限制。在分面函数中将 scales 参数设置为 "free_x",可以让各列的 x 轴使用不同的刻度;设置为 "free_y",可以让各行的 y 轴使用不同的刻度;设置为 "free",则两者都可以变化。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point() +
    facet_grid(drv ~ cyl, scales = "free")

    练习

    • 如果对连续变量进行分面,会发生什么?

    如果对连续变量进行分面,通常会先把它离散化(分成若干组),然后再分别绘图。

    不过要注意:

    1. 分面更适合分类变量
    2. 连续变量如果直接用于分面,通常不太合适
    3. 常见做法是先把连续变量分箱,例如按区间分组,再用 facet_wrap() 或 facet_grid() 分面
    ggplot(mpg, aes(x = displ, y = cty)) +
    geom_point() +
    facet_wrap(~ hwy)
    • 上面使用 facet_grid(drv ~ cyl) 绘制的图中,空白单元格是什么意思?运行下面的代码。它们与生成的图有什么关系?

    在 facet_grid() 里,. 表示这一边不使用变量分面

    1. drv ~ .:只按行分面,不按列分面
    2. . ~ cyl:只按列分面,不按行分面

    统计变换

    考虑一个使用 geom_bar() 或 geom_col() 绘制的基本条形图。下图显示了 diamonds 数据集中按切工分组的钻石总数。diamonds 数据集位于 ggplot2 包中,包含约 54,000 颗钻石的信息,包括每颗钻石的价格、克拉、颜色、净度和切工。该图表显示,高质量切工的钻石数量比低质量切工的钻石数量更多。

    ggplot(diamonds, aes(x = cut)) +
    geom_bar()

    在 x 轴上,图表显示的是 cut,这是 diamonds 中的一个变量。在 y 轴上,显示的是 count,但 count 并不是 diamonds 中的变量!count 是从哪里来的?许多图形,比如散点图,绘制的是数据集的原始值。其他图形,比如条形图,会先计算新的值再进行绘图:

    条形图、直方图和频率多边形会先对数据进行分箱,然后绘制每个箱中的计数,也就是落入每个箱内的点的数量。

    平滑曲线会先对数据拟合一个模型,然后绘制模型的预测值。

    箱线图会计算分布的五数概括,然后把这个概括以一种特殊格式的箱体展示出来。

    用于为图形计算新值的算法称为 stat,即 statistical transformation(统计变换)的缩写。图中展示了这一过程在 geom_bar() 中是如何工作的。

    你可以通过查看 stat 参数的默认值来了解某个 geom 使用了哪种统计变换。例如,?geom_bar 显示 stat 的默认值是 “count”,这意味着 geom_bar() 使用的是 stat_count()stat_count() 的文档和 geom_bar() 在同一页面中。向下滚动,你会看到一个名为 “Computed variables” 的部分,它说明它会计算两个新变量:count 和 prop

    每个 geom 都有一个默认的 stat;每个 stat 也有一个默认的 geom。这意味着你通常可以直接使用 geom,而不必担心底层的统计变换。不过,你可能需要显式使用 stat,主要有三个原因:

    你可能想覆盖默认的 stat。在下面的代码中,我们把 geom_bar() 的 stat 从 count(默认值)改为 identity。这样就可以把柱子的高度映射到 y 变量的原始值。

    diamonds |>
    count(cut) |>
    ggplot(aes(x = cut, y = n)) +
    geom_bar(stat = "identity")

    你可能想要覆盖从变换后的变量到美学属性的默认映射。例如,你可能想显示一个比例条形图,而不是计数条形图:

    ggplot(diamonds, aes(x = cut, y = after_stat(prop), group = 1)) +
    geom_bar()

    要找出 stat 可以计算的可能变量,请查看 geom_bar() 帮助文档中标题为 “computed variables” 的部分。

    你可能希望在代码中更突出地显示统计变换。例如,你可以使用 stat_summary(),它会对每个唯一的 x 值汇总 y 值,从而突出你正在计算的汇总结果:

    # 使用 ggplot2 对 diamonds 数据集进行绘图
    ggplot(diamonds) +
    # 添加统计汇总层
    stat_summary(
    # 指定映射:x 轴为 cut,y 轴为 depth
    aes(x = cut, y = depth),
    # 计算每组的最小值
    fun.min = min,
    # 计算每组的最大值
    fun.max = max,
    # 计算每组的中位数
    fun = median
    )

    位置调整

    柱状图还有一个魔法。你可以使用 color 美学映射为柱状图着色,不过更有用的是使用 fill 美学映射:

    ggplot(mpg, aes(x = drv, color = drv)) +
    geom_bar()
    ggplot(mpg, aes(x = drv, fill = drv)) +
    geom_bar()

    注意,如果你将 fill 美学映射到另一个变量,比如 class,会发生什么:柱子会自动堆叠。每个有颜色的矩形代表 drv 和 class 的一种组合。

    ggplot(mpg, aes(x = drv, fill = class)) +
    geom_bar()

    堆叠是通过 position 参数指定的位置调整自动完成的。如果你不想要堆叠柱状图,可以使用另外三个选项之一:"identity""dodge" 或 "fill"

    ggplot(mpg, aes(x = drv, fill = class)) +
    geom_bar(alpha = 1/5, position = "identity")
    ggplot(mpg, aes(x = drv, color = class)) +
    geom_bar(fill = NA, position = "identity")

    identity 位置调整对二维几何对象(如点)更有用,因为它是默认设置。

    position = “fill” 的作用类似于堆叠,但会让每组堆叠柱子的高度相同。这样更容易比较不同组之间的比例。

    position = “dodge” 会将重叠的对象并排放置。这样更容易比较各个单独的数值。

    ggplot(mpg, aes(x = drv, fill = class)) +
    geom_bar(position = "fill")
    ggplot(mpg, aes(x = drv, fill = class)) +
    geom_bar(position = "dodge")

    还有一种调整对柱状图没什么用,但对散点图非常有用。回想一下我们的第一个散点图。你注意到了吗?尽管数据集中有 234 个观测值,但图中只显示了 126 个点。

    hwy 和 displ 的底层数值经过了四舍五入,因此这些点会显示在一个网格上,许多点彼此重叠。这个问题称为过度绘制(overplotting)。这种排列方式会使数据分布变得难以观察。数据点是均匀分布在整张图中,还是存在某个特殊的 hwy 和 displ 组合包含 109 个值?

    你可以通过将位置调整设置为 “jitter” 来避免这种网格状排列。position = "jitter" 会为每个点添加少量随机噪声。这样可以把点分散开,因为不太可能有两个点获得完全相同数量的随机噪声。

    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(position = "jitter")

    添加随机性似乎是一种奇怪的改进图形的方法,但虽然它会让你的图在小尺度上不那么精确,却会让你的图在大尺度上更能揭示信息。由于这是一个非常有用的操作,ggplot2 为 geom_point(position = "jitter") 提供了一个简写:geom_jitter()

    要了解更多关于位置调整的信息,请查看每种调整对应的帮助页面:?position_dodge?position_fill?position_identity?position_jitter 和 ?position_stack

    坐标系

    坐标系统大概是 ggplot2 中最复杂的部分。默认的坐标系统是笛卡尔坐标系,其中 x 和 y 的位置相互独立,用来确定每个点的位置。还有另外两种坐标系统,有时也会很有用。

    coord_quickmap() 会为地理地图正确设置纵横比。如果你使用 ggplot2 绘制空间数据,这一点非常重要。

    # 加载新西兰地图数据
    nz <- map_data("nz")
    # 绘制普通笛卡尔坐标下的新西兰地图
    ggplot(nz, aes(x = long, y = lat, group = group)) +
    geom_polygon(fill = "white", color = "black")
    # 绘制使用快速地图坐标系的新西兰地图
    # coord_quickmap() 会尽量保持地图的纵横比例,避免地理形状失真
    ggplot(nz, aes(x = long, y = lat, group = group)) +
    geom_polygon(fill = "white", color = "black") +
    coord_quickmap()

    coord_polar() 使用极坐标。极坐标揭示了条形图和 Coxcomb 图之间的一个有趣联系。

    # 创建一个基础柱状图对象 bar
    # 使用 diamonds 数据集,
    # x 轴为 clarity,
    # 并用 clarity 填充颜色
    # 不显示图例,柱宽设为 1
    # 同时设置图形纵横比为 1:1
    bar <- ggplot(data = diamonds) +
    geom_bar(
    mapping = aes(x = clarity, fill = clarity),
    show.legend = FALSE,
    width = 1
    ) +
    theme(aspect.ratio = 1)
    # 显示柱状图
    bar
    # 将柱状图翻转为水平图
    bar + coord_flip()
    # 将柱状图转换为极坐标图
    bar + coord_polar()

    图形的分层语法

    我们可以通过添加位置调整、统计变换、坐标系统和分面来扩展你在前面学到的绘图模板:

    ggplot(data = <DATA>) +
    <GEOM_FUNCTION>(
    mapping = aes(<MAPPINGS>),
    stat = <STAT>,
    position = <POSITION>
    ) +
    <COORDINATE_FUNCTION> +
    <FACET_FUNCTION>

    我们的新模板包含七个参数,也就是模板中括号里的那些词。实际上,你很少需要同时提供全部七个参数来绘制图形,因为 ggplot2 会为除数据、映射和几何对象函数之外的所有内容提供有用的默认值。

    模板中的这七个参数构成了图形语法(grammar of graphics),这是一种用于构建图形的形式化系统。图形语法基于这样一个观点:你可以将任何图形唯一地描述为数据集、几何对象、映射、统计变换、位置调整、坐标系统、分面方案和主题的组合。

    为了理解这一点,不妨考虑如何从头构建一个基础图形:你可以先从一个数据集开始,然后将其转换为你想要展示的信息(通过统计变换)。接着,你可以选择一种几何对象来表示转换后数据中的每个观测值。然后,你可以使用几何对象的美学属性来表示数据中的变量,并将每个变量的取值映射到某种美学属性的不同水平。这些步骤如图所示。之后,你还可以选择一个坐标系统来放置这些几何对象,并利用对象的位置(它本身也是一种美学属性)来显示 x 和 y 变量的值。

    此时,你已经拥有一个完整的图形,但你还可以进一步调整坐标系统中各个几何对象的位置(位置调整),或者将图形拆分成多个子图(分面)。你也可以通过添加一个或多个额外图层来扩展图形,其中每个额外图层都使用一个数据集、一个几何对象、一组映射、一个统计变换和一个位置调整。

    你可以用这种方法构建你能想象到的任何图形。换句话说,你可以使用本章学到的代码模板来构建成千上万种独特的图形。

  • R语言工作流程:获取帮助

    https://wp.me/p80aHo-2vT

    Table of Contents

    随着你开始将本书中介绍的技术应用到自己的数据上,你很快就会遇到一些我们没有回答的问题。本节介绍一些如何获取帮助以及持续学习的技巧。

    Google 是你的朋友

    如果你卡住了,先从 Google 开始。通常在查询中加上 “R” 就足以将结果限制在相关内容:如果搜索没有帮助,这通常意味着没有可用的 R 特定结果。此外,加入 “ggplot2”之类的包名,也能进一步缩小范围,让你找到更熟悉的代码,例如,“如何在 R 中制作箱线图” 与 “如何在带有 ggplot2 的 R 中制作箱线图”。Google 对错误信息尤其有用。如果你看到一条错误信息却完全不知道它是什么意思,试着去 Google 一下!很可能过去已经有人被它困扰过,网上某处会有帮助。(如果错误信息不是英文,运行 Sys.setenv(LANGUAGE = "en") 后重新执行代码;你更有可能找到针对英文错误信息的帮助。)

    如果 Google 没能帮上忙,试试 Stack Overflow。先花一点时间搜索已有答案,并加入 [R] 来限制你的搜索范围,使其只包含使用 R 的问题和答案。

    制作 reprex

    如果你的 Google 搜索没有找到任何有用的信息,那么准备一个 reprex(即 minimal reproducible example,最小可重现示例)是个非常好的主意。一个好的 reprex 能让别人更容易帮助你,而且你常常会在制作它的过程中自己就把问题弄清楚。创建 reprex 主要分两部分:

    第一,你需要让代码可重现。这意味着你必须捕捉所有内容,也就是包含任何 library() 调用,并创建所有必要的对象。确保这一点最简单的方法是使用 reprex 包。

    第二,你需要让它尽可能精简。删掉所有与问题没有直接关系的内容。这通常意味着要创建一个比你在实际中面对的对象小得多、简单得多的 R 对象,甚至直接使用内置数据。

    听起来工作量很大!确实可能如此,但回报也很大:

    80% 的情况下,制作一个优秀的 reprex 就能揭示问题的根源。把一个自包含且尽可能精简的示例写出来,常常神奇地就能让你自己回答自己的问题。

    另外 20% 的情况下,你也会把问题的核心以一种别人很容易上手的方式呈现出来。这会大大提高你获得帮助的机会!

    手工创建 reprex 时,很容易不小心遗漏某些东西,导致你的代码无法在别人的电脑上运行。要避免这个问题,可以使用 reprex 包。假设你把下面这段代码复制到剪贴板上(或者,在 RStudio Server 或 Cloud 中,选中它):

    y <- 1:4
    mean(y)
    library(reprex)
    reprex()

    一个渲染精美的 HTML 预览会显示在 RStudio 的 Viewer 中(如果你正在使用 RStudio),否则会在你的默认浏览器中打开。reprex 会自动复制到你的剪贴板中(在 RStudio Server 或 Cloud 中,你需要自己手动复制):

    ``` r
    y <- 1:4
    mean(y)
    #> [1] 2.5
    ```

    这段文本采用一种特殊的格式,称为 Markdown,可以粘贴到 Stack Overflow 或 GitHub 等网站上,它们会自动将其渲染成代码样式。下面是这段 Markdown 在 GitHub 上渲染后的样子:

    任何人都可以立即复制、粘贴并运行它。

    要让你的示例可重现,你需要包含三样东西:所需的包、数据和代码。

    包应该在脚本顶部加载,这样就很容易看出示例需要哪些包。这也是检查你是否使用了每个包的最新版本的好时机;你可能已经发现了一个自从你安装或上次更新该包以来就已修复的 bug。

    包含数据最简单的方法是使用 dput() 生成重建它所需的 R 代码。例如,要在 R 中重建 mtcars 数据集,请执行以下步骤:

    在 R 中运行 dput(mtcars)
    复制输出
    在 reprex 中输入 mtcars <-,然后粘贴。
    尽量使用能仍然揭示问题的最小数据子集。

    花一点时间确保你的代码易于他人阅读:

    确保你使用了空格,并且变量名简洁但有信息量。

    使用注释来指出问题所在。

    尽最大努力删除所有与问题无关的内容。

    你的代码越短,就越容易理解,也越容易修复。

    最后,通过启动一个全新的 R 会话并复制粘贴你的脚本,检查你是否真的制作出了一个可重现的示例。

    创建 reprex 并不简单,学会制作好的、真正最小的 reprex 需要一些练习。不过,学会在提问时包含代码,并投入时间让它可重现,会在你学习和掌握 R 的过程中持续带来回报。

    投资于自己

    你还应该花一些时间提前为解决问题做准备,而不是等问题出现后再去应对。每天花一点时间学习 R,长期来看会带来丰厚回报。一个方法是关注 tidyverse 团队在 tidyverse 博客上的动态。若想更广泛地了解 R 社区的情况,我们建议阅读 R Weekly:这是一个社区协作项目,每周汇总 R 社区中最有趣的新闻。

    总结

    本章结束了本书的 Whole Game 部分。到现在为止,你已经看到了数据科学流程中最重要的部分:可视化、转换、整理和导入。现在你已经对整个流程有了整体性的理解,接下来我们将开始深入讲解各个小部分的细节。

    本书的下一部分,Visualize,将更深入地探讨图形语法以及使用 ggplot2 创建数据可视化,展示如何运用你目前学到的工具进行探索性数据分析,并介绍用于制作便于交流的图表的良好实践。

  • R语言数据导入

    https://wp.me/p80aHo-2uJ

    Table of Contents

    引言

    使用 R 包提供的数据进行练习,是学习数据科学工具的好方法,但你总会在某个时候想把所学应用到自己的数据上。在本章中,你将学习将数据文件读取到 R 的基础知识。

    具体来说,重点介绍读取纯文本矩形文件。我们将从处理列名、类型和缺失数据等特性的实用建议开始。然后,你将学习如何一次从多个文件读取数据,以及如何将 R 中的数据写入文件。最后,你将学习如何在 R 中手工创建数据框。

    从文件读取数据

    首先,我们将重点介绍最常见的矩形数据文件类型:CSV,即逗号分隔值(comma-separated values)的缩写。下面是一个简单的 CSV 文件示例。第一行通常称为表头行,提供列名,接下来的六行提供数据。各列之间用逗号分隔,也就是用逗号作为分隔符。

    Student ID,Full Name,favourite.food,mealPlan,AGE
    1,Sunil Huffmann,Strawberry yoghurt,Lunch only,4
    2,Barclay Lynn,French fries,Lunch only,5
    3,Jayendra Lyne,N/A,Breakfast and lunch,7
    4,Leon Rossini,Anchovies,Lunch only,
    5,Chidiegwu Dunkel,Pizza,Breakfast and lunch,five
    6,Güvenç Attila,Ice cream,Lunch only,6
    Student IDFull Namefavourite.foodmealPlanAGE
    1Sunil HuffmannStrawberry yoghurtLunch only4
    2Barclay LynnFrench friesLunch only5
    3Jayendra LyneN/ABreakfast and lunch7
    4Leon RossiniAnchoviesLunch onlyNA
    5Chidiegwu DunkelPizzaBreakfast and lunchfive
    6Güvenç AttilaIce creamLunch only6

    我们可以使用 read_csv() 将这个文件读入 R。第一个参数最重要:它是文件路径。你可以把路径理解为文件的地址:这个文件名为 students.csv,位于 data 文件夹中。

    students <- read_csv("data/students.csv")
    #> Rows: 6 Columns: 5
    #> ── Column specification ─────────────────────────────────────────────────────
    #> Delimiter: ","
    #> chr (4): Full Name, favourite.food, mealPlan, AGE
    #> dbl (1): Student ID
    #>
    #> ℹ Use `spec()` to retrieve the full column specification for this data.
    #> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

    如果你的项目中有一个 data 文件夹,并且其中包含 students.csv 文件,那么上面的代码就可以正常运行。你可以从 https://pos.it/r4ds-students-csv 下载 students.csv 文件,或者直接使用以下方式从该 URL 读取它:

    library(readr)
    students <- read_csv("https://pos.it/r4ds-students-csv")

    当你运行 read_csv() 时,它会打印一条消息,告诉你数据的行数和列数、所使用的分隔符,以及列规范(按列中所包含的数据类型组织的列名)。它还会打印一些关于如何获取完整列规范的信息,以及如何静默这条消息的信息。

    实用建议

    当你读入数据后,第一步通常是以某种方式对其进行转换,使其在后续分析中更容易处理。让我们带着这一点,再来看看学生数据。

    students
    #> # A tibble: 6 × 5
    #> `Student ID` `Full Name` favourite.food mealPlan AGE
    #> <dbl> <chr> <chr> <chr> <chr>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne N/A Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only <NA>
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch five
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    在 favourite.food 这一列中,有一堆食物项,之后还有字符字符串 N/A,它本应是一个真正的 NA,R 会将其识别为“不可用”。我们可以使用 na 参数来解决这个问题。默认情况下,read_csv() 只会将该数据集中的空字符串(””)识别为 NA,而我们希望它也能识别字符字符串 “N/A”。

    students <- read_csv("https://pos.it/r4ds-students-csv", na = c("N/A", ""))

    你可能还会注意到,Student ID 和 Full Name 两列被反引号包围着。这是因为它们包含空格,打破了 R 对变量名的通常规则;它们属于非语法名称。要引用这些变量,你需要用反引号将它们括起来,`:

    students |>
    rename(
    student_id = `Student ID`,
    full_name = `Full Name`
    )
    #> # A tibble: 6 × 5
    #> student_id full_name favourite.food mealPlan AGE
    #> <dbl> <chr> <chr> <chr> <chr>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only <NA>
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch five
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    另一种方法是使用 janitor::clean_names(),借助一些启发式规则一次性把它们全部转换为 snake case。

    library(janitor)
    students |> janitor::clean_names()
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <chr> <chr>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only <NA>
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch five
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    在读入数据之后,另一个常见任务是考虑变量类型。例如,meal_plan 是一个分类变量,它有一组已知的可能取值,在 R 中应将其表示为因子:

    library(dplyr)
    students |>
    janitor::clean_names() |>
    mutate(meal_plan = factor(meal_plan))
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <fct> <chr>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only <NA>
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch five
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    请注意,meal_plan 变量中的值保持不变,但变量名下方标示的变量类型已经从字符型(<chr>)变为因子型(<fct>)。

    在分析这些数据之前,你可能需要先修正 age 列。目前,age 是一个字符型变量,因为其中一个观测值被写成了 five,而不是数字 5。

    students <- students |>
    janitor::clean_names() |>
    mutate(
    meal_plan = factor(meal_plan),
    age = parse_number(if_else(age == "five", "5", age))
    )
    students
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <fct> <dbl>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only NA
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch 5
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    这里一个新的函数是 if_else(),它有三个参数。第一个参数 test 应该是一个逻辑向量。当 test 为 TRUE 时,结果将包含第二个参数 yes 的值;当它为 FALSE 时,结果将包含第三个参数 no 的值。这里我们是在说,如果 age 是字符字符串 "five",就把它改成 "5";如果不是,就保持为 age

    其他参数

    还有另外几个重要参数需要提及,而如果我们先向你展示一个实用技巧,它们会更容易说明:read_csv() 可以读取你创建并格式化成 CSV 文件样式的文本字符串:

    read_csv(
    "a,b,c
    1,2,3
    4,5,6"
    )

    通常,read_csv() 会使用数据的第一行作为列名,这是一种非常常见的约定。但文件顶部包含几行元数据也并不少见。你可以使用 skip = n 来跳过前 n 行,或者使用 comment = "#" 来删除所有以(例如)# 开头的行:

    read_csv(
    "The first line of metadata
    The second line of metadata
    x,y,z
    1,2,3",
    skip = 2
    )
    read_csv(
    "# A comment I want to skip
    x,y,z
    1,2,3",
    comment = "#"
    )

    在其他情况下,数据可能没有列名。你可以使用 col_names = FALSE 告诉 read_csv() 不要将第一行当作标题,而是按顺序将它们标记为 X1 到 Xn

    read_csv(
    "1,2,3
    4,5,6",
    col_names = FALSE
    )

    或者,你也可以向 col_names 传入一个字符向量,它将被用作列名:

    read_csv(
    "1,2,3
    4,5,6",
    col_names = c("x", "y", "z")
    )

    这些参数已经足够让你读取在实践中遇到的大多数 CSV 文件。(至于其余的情况,你需要仔细检查你的 .csv 文件,并阅读 read_csv() 其他许多参数的文档。)

    其他文件类型

    一旦你掌握了 read_csv(),使用 readr 的其他函数就很简单了;关键只是知道该使用哪个函数:

    read_csv2() 读取以分号分隔的文件。这类文件使用 ; 而不是 , 来分隔字段,并且在以 , 作为小数点标记的国家很常见。

    read_tsv() 读取以制表符分隔的文件。

    read_delim() 读取使用任意分隔符的文件;如果你没有指定分隔符,它会尝试自动猜测。

    read_fwf() 读取固定宽度文件。你可以用 fwf_widths() 按字段宽度指定,或用 fwf_positions() 按字段位置指定。

    read_table() 读取固定宽度文件的一种常见变体,其中各列由空白分隔。

    read_log() 读取 Apache 风格的日志文件。

    练习

    • 你会使用什么函数来读取一个字段由“|”分隔的文件?
    read_delim()
    • 有时 CSV 文件中的字符串包含逗号。为防止它们引发问题,需要用引号字符将其括起来,例如 " 或 '。默认情况下,read_csv() 假定引号字符是 "。要将下面的文本读入数据框,你需要为 read_csv() 指定什么参数?”x,y\n1,’a,b’”
    read_csv(
    "x,y\n1,'a,b'",
    quote = "'"
    )

    控制列类型

    CSV 文件不包含任何关于每个变量类型的信息(即它是逻辑型、数值型、字符串型等),因此 readr 会尝试猜测其类型。本节将描述这种猜测过程是如何工作的,如何解决一些导致其失败的常见问题,以及在需要时如何自行提供列类型。最后,我们还会提到一些通用策略;当 readr 失败得很严重、你需要更深入了解文件结构时,这些策略会很有帮助。

    猜测类型

    readr 使用一种启发式方法来判断列类型。对于每一列,它会从第一行到最后一行等间隔抽取 1,000 行的值,并忽略缺失值。然后它会依次回答以下问题:

    • 它是否只包含 F、T、FALSE 或 TRUE(忽略大小写)?如果是,那么它是逻辑型。
    • 它是否只包含数字(例如 1、-4.5、5e6、Inf)?如果是,那么它是数值型。
    • 它是否符合 ISO8601 标准?如果是,那么它是日期或日期时间。

    否则,它一定是字符串。
    你可以在下面这个简单示例中看到这种行为:

    read_csv("
    logical,numeric,date,string
    TRUE,1,2021-01-15,abc
    false,4.5,2021-02-15,def
    T,Inf,2021-02-16,ghi
    ")

    如果你有一个干净的数据集,这种启发式方法效果很好;但在现实中,你会遇到各种奇怪而又“美丽”的失败情况。

    缺失值、列类型和问题

    列检测最常见的失败方式是:某一列包含了意料之外的值,结果你得到的是字符型列,而不是更具体的类型。造成这种情况最常见的原因之一,是缺失值使用了 readr 预期之外的其他标记方式来记录。

    以这个简单的单列表 CSV 文件为例:

    simple_csv <- "
    x
    10
    .
    20
    30"

    如果我们不添加任何额外参数来读取它,x 会变成字符列:

    read_csv(simple_csv)

    在这个非常小的例子里,你很容易就能看出缺失值“.”。但如果你有成千上万行数据,其中只有少数缺失值以散布其中的“.”来表示,又会怎样呢?一种做法是告诉 readr,x 是数值型列,然后看看它在哪里失败。你可以使用 col_types 参数来做到这一点,它接受一个命名列表,其中名称与 CSV 文件中的列名相匹配:

    df <- read_csv(
    simple_csv,
    col_types = list(x = col_double())
    )
    df

    现在 read_csv() 会报告存在一个问题,并告诉我们可以通过 problems() 了解更多信息:

    这告诉我们,在第 3 行、第 1 列出现了一个问题:readr 期望得到一个双精度数,但实际得到了一个“.”。这说明该数据集使用“.”表示缺失值。于是我们设置 na = ".",自动推断就成功了,得到了我们想要的数值列:

    read_csv(simple_csv, na = ".")

    列类型

    readr 为你提供了总共九种可用的列类型:

    • col_logical() 和 col_double() 用于读取逻辑值和实数。它们相对很少需要手动指定(如上所示的情况除外),因为 readr 通常会自动为你推断出来。
    • col_integer() 用于读取整数。我们很少区分整数和双精度数,因为它们在功能上等价,但显式读取整数有时会很有用,因为它们占用的内存只有双精度数的一半。
    • col_character() 用于读取字符串。当某一列是数值型标识符时,显式指定它会很有用,也就是那种由一长串数字组成、用于标识某个对象,但对其进行数学运算并没有意义的数据。例如电话号码、社会保险号、信用卡号码等。
    • col_factor()col_date() 和 col_datetime() 分别用于创建因子、日期和日期时间。
    • col_number() 是一种更宽松的数值解析器,它会忽略非数值部分,尤其适合处理货币。你会在第 13 章进一步了解它。
    • col_skip() 会跳过某一列,因此不会将其包含在结果中;如果你有一个很大的 CSV 文件,但只想使用其中部分列,这样可以加快读取速度。

    你也可以通过将 list() 改为 cols() 并指定 .default,来覆盖默认的启发式类型推断:

    another_csv <- "
    x,y,z
    1,2,3"
    read_csv(
    another_csv,
    col_types = cols(.default = col_character())
    )

    另一个有用的辅助函数是 cols_only(),它只会读取你指定的列:

    read_csv(
    another_csv,
    col_types = cols_only(x = col_character())
    )

    从多个文件读取数据

    有时你的数据会分散在多个文件中,而不是包含在一个文件里。例如,你可能有按月划分的销售数据,每个月的数据都放在单独的文件中:01-sales.csv 表示一月,02-sales.csv 表示二月,03-sales.csv 表示三月。使用 read_csv(),你可以一次性读取这些数据,并将它们按顺序叠加到同一个数据框中。

    sales_files <- c("data/01-sales.csv", "data/02-sales.csv", "data/03-sales.csv")
    read_csv(sales_files, id = "file")

    如果你的 CSV 文件位于项目中的 data 文件夹里,上面的代码同样可以正常工作。你可以从 https://pos.it/r4ds-01-sales、https://pos.it/r4ds-02-sales 和 https://pos.it/r4ds-03-sales 下载这些文件,或者直接用下面的方法读取它们:

    sales_files <- c(
    "https://pos.it/r4ds-01-sales",
    "https://pos.it/r4ds-02-sales",
    "https://pos.it/r4ds-03-sales"
    )
    read_csv(sales_files, id = "file")

    id 参数会在结果数据框中新增一列名为 file,用于标识数据来自哪个文件。当你读取的文件没有可用于追踪观测值原始来源的标识列时,这一点尤其有用。

    如果你有很多文件需要读取,把它们的文件名逐个写成列表会很麻烦。相反,你可以使用基础函数 list.files(),通过匹配文件名中的模式来帮你找到这些文件。

    sales_files <- list.files("data", pattern = "sales\\.csv$", full.names = TRUE)
    sales_files
    #> [1] "data/01-sales.csv" "data/02-sales.csv" "data/03-sales.csv"

    Writing to a file

    readr 还提供了两个用于将数据写回磁盘的实用函数:write_csv() 和 write_tsv()。这些函数最重要的参数是 x(要保存的数据框)和 file(保存位置)。你还可以通过 na 指定缺失值的写入方式,并且可以选择是否追加到已有文件中。

    write_csv(students, "students.csv")

    现在让我们把那个 CSV 文件再读回来。请注意,你刚刚设置的变量类型信息在保存为 CSV 时会丢失,因为你又是从一个纯文本文件重新开始读取的:

    students
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <fct> <dbl>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only NA
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch 5
    #> 6 6 Güvenç Attila Ice cream Lunch only 6
    write_csv(students, "students-2.csv")
    read_csv("students-2.csv")
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <chr> <dbl>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only NA
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch 5
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    这会使 CSV 在缓存中间结果时有些不太可靠——你每次加载时都需要重新创建列规格。有两个主要的替代方案:

    write_rds() 和 read_rds() 是对基础函数 readRDS() 和 saveRDS() 的统一封装。它们将数据存储为 R 的自定义二进制格式,称为 RDS。这意味着当你重新加载对象时,加载的是你存储时那个完全相同的 R 对象。

    write_rds(students, "students.rds")
    read_rds("students.rds")
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <fct> <dbl>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne <NA> Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only NA
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch 5
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    arrow 包允许你读写 Parquet 文件,这是一种可以跨编程语言共享的快速二进制文件格式。

    library(arrow)
    write_parquet(students, "students.parquet")
    read_parquet("students.parquet")
    #> # A tibble: 6 × 5
    #> student_id full_name favourite_food meal_plan age
    #> <dbl> <chr> <chr> <fct> <dbl>
    #> 1 1 Sunil Huffmann Strawberry yoghurt Lunch only 4
    #> 2 2 Barclay Lynn French fries Lunch only 5
    #> 3 3 Jayendra Lyne NA Breakfast and lunch 7
    #> 4 4 Leon Rossini Anchovies Lunch only NA
    #> 5 5 Chidiegwu Dunkel Pizza Breakfast and lunch 5
    #> 6 6 Güvenç Attila Ice cream Lunch only 6

    Parquet 通常比 RDS 快得多,而且可以在 R 之外使用,但它需要 arrow 包。

    数据输入

    有时候,你需要在 R 脚本里通过少量手工录入来“手动”组装一个 tibble。这里有两个很有用的函数可以帮助你完成这件事,它们的区别在于你是按列还是按行来构建 tibble。tibble() 是按列工作的:

    tibble(
    x = c(1, 2, 5),
    y = c("h", "m", "g"),
    z = c(0.08, 0.83, 0.60)
    )
    #> # A tibble: 3 × 3
    #> x y z
    #> <dbl> <chr> <dbl>
    #> 1 1 h 0.08
    #> 2 2 m 0.83
    #> 3 5 g 0.6

    按列排列数据会让人难以看出各行之间的关系,因此可以使用 tribble(),即转置 tibble 的简称,它允许你按行排列数据。tribble() 是专为在代码中输入数据而设计的:列名以 ~ 开头,条目之间用逗号分隔。这使得少量数据可以以一种易读的形式排布出来:

    tribble(
    ~x, ~y, ~z,
    1, "h", 0.08,
    2, "m", 0.83,
    5, "g", 0.60
    )
    #> # A tibble: 3 × 3
    #> x y z
    #> <dbl> <chr> <dbl>
    #> 1 1 h 0.08
    #> 2 2 m 0.83
    #> 3 5 g 0.6

    总结

    你已经学习了如何使用 read_csv() 加载 CSV 文件,以及如何使用 tibble() 和 tribble() 自己输入数据。你已经了解了 CSV 文件的工作方式、可能遇到的一些问题,以及如何克服这些问题。

  • R语言工作流程:脚本和项目

    https://wp.me/p80aHo-2u8

    Table of Contents

    脚本

    到目前为止,你一直使用控制台来运行代码。这是一个很好的起点,但当你创建更复杂的 ggplot2 图形以及更长的 dplyr 管道时,你会发现这里很快就会变得非常拥挤。为了给自己更多的工作空间,请使用脚本编辑器。你可以通过点击“File(文件)”菜单、选择“New File(新建文件)”,然后选择“R script(R 脚本)”,或者使用键盘快捷键 Cmd/Ctrl + Shift + N 来打开它。现在你会看到四个面板,就像图 6.1 所示一样。脚本编辑器非常适合用来试验你的代码。当你想要修改某些内容时,你不必重新输入整段内容;你只需要编辑脚本并重新运行即可。而当你已经编写出能够正常运行、并实现你想要效果的代码之后,你就可以将其保存为脚本文件,之后随时方便地返回继续使用。

    运行代码

    脚本编辑器是构建复杂的 ggplot2 图形或进行长串 dplyr 处理的绝佳场所。要高效地使用脚本编辑器,关键在于记住最重要的一个键盘快捷键:Cmd/Ctrl + Enter。这个快捷键会在控制台中执行当前的 R 表达式。例如,看看下面这段代码。

    # 1) 加载需要的 R 包
    library(dplyr) # 提供数据处理/管道(filter、group_by、summarize 等)函数
    library(nycflights13) # 提供 nycflights13 数据集(航班数据)
    # 2) 筛选:只保留“未被取消”的航班记录
    # - dep_delay:起飞延误(如果航班被取消,通常会对应为 NA)
    # - arr_delay:到达延误(如果航班被取消,通常也会对应为 NA)
    # - 这里的作用是:排除 dep_delay 或 arr_delay 为 NA 的行
    not_cancelled <- flights |>
    filter(!is.na(dep_delay)█, !is.na(arr_delay))
    # 3) 分组:按航班日期(年、月、日)分成多个组
    # group_by(year, month, day) 让后续的 summarize 在每个日期组内分别计算
    not_cancelled |>
    group_by(year, month, day) |>
    summarize(
    # 4) 汇总:对每个日期组计算平均起飞延误
    # - mean(dep_delay):该日期组内所有 dep_delay 的平均值
    # - summarize 会返回一个“每组一行”的结果数据框
    mean = mean(dep_delay)
    )

    如果你的光标位于 █ 处,按下 Cmd/Ctrl + Enter 将运行生成 not_cancelled 的完整命令。它还会把光标移动到下一条语句(以 not_cancelled |> 开头)。因此,你可以通过反复按 Cmd/Ctrl + Enter,轻松地一步步浏览/执行你的完整脚本。

    我们建议你始终在脚本的开头先写好你需要的包。这样,当你把代码分享给别人时,他们就能很容易地看出需要安装哪些包。不过请注意:你在分享的脚本中绝对不要包含 install.packages()。因为如果别人没有足够小心,你把脚本交给他们后,它可能会在他们的电脑上做出一些改变——这是不太体贴的做法。

    在后续章节中学习时,我们强烈建议你从脚本编辑器开始,并练习你的键盘快捷键。久而久之,你用这种方式把代码发送到控制台会变得非常自然,甚至都不用再去刻意思考。

    保存与命名

    RStudio 会在你退出时自动保存脚本编辑器中的内容,并在你重新打开时自动重新加载。然而,还是建议你避免使用 Untitled1、Untitled2、Untitled3 之类的文件名,而是要把你的脚本保存起来,并给它们起有信息量的名字。

    给文件起名为 code.R 或 myscript.R 可能会很诱人,但在选择文件名之前,你应该多想一想。文件命名有三个重要原则如下:

    • 文件名应当可被机器读取:避免使用空格、符号和特殊字符。不要依赖大小写敏感来区分不同文件。
    • 文件名应当便于人类阅读:用文件名来描述文件里包含的内容。
    • 文件名应当与默认排序方式友好兼容:让文件名以数字开头,这样按字母顺序排序时,它们就会按照实际使用的顺序排列。

    例如,假设你的项目文件夹中有如下这些文件。

    alternative model.R
    code for exploratory analysis.r
    finalreport.qmd
    FinalReport.qmd
    fig 1.png
    Figure_02.png
    model_first_try.R
    run-first.r
    temp.txt

    这里有多种问题:很难弄清应该先运行哪个文件,文件名里包含空格,有两个文件虽然名字相同但大小写不同(finalreport vs. FinalReport1),另外还有一些名字不能很好地描述它们的内容(run-first 和 temp)。

    下面是为同一组文件起名和整理的更好方法:

    01-load-data.R
    02-exploratory-analysis.R
    03-model-approach-1.R
    04-model-approach-2.R
    fig-01.png
    fig-02.png
    report-2022-03-20.qmd
    report-2022-04-02.qmd
    report-draft-notes.txt

    对关键脚本进行编号可以清楚地表明运行它们的顺序,而一致的命名规则也让你更容易看出哪些地方是变化的。此外,图形的标注也采用了类似的方式,报告通过文件名中包含的日期来区分,temp 也被重命名为 report-draft-notes,以更准确地描述其内容。如果某个目录里有很多文件,建议在组织方面再更进一步:把不同类型的文件(脚本、图形等)放到不同的目录中。

    项目

    总有一天,你需要退出 R,去做别的事情,然后稍后再回来继续分析。总有一天,你会同时进行多个分析,并希望把它们彼此分开。总有一天,你需要把外部世界的数据导入 R,并把数值结果和图形从 R 再输出到外部世界。

    为了应对这些现实生活中的情况,你需要做出两个决定:

    1. 什么是事实来源?你会把什么保存为对所发生事情的持久记录?
    2. 你的分析存放在哪里?

    什么是权威来源?

    作为初学者,你可以依赖当前的 Environment 来保存你在分析过程中创建的所有对象。不过,为了更方便地处理更大的项目或与他人协作,你的事实来源应该是 R 脚本。有了 R 脚本(以及你的数据文件),你就可以重建 Environment。只有 Environment 的话,要重建你的 R 脚本就困难得多:你要么得凭记忆重写大量代码(这一路上难免会出错),要么就得仔细翻查你的 R 历史记录。

    为了帮助你把 R 脚本作为分析的事实来源,我们强烈建议你指示 RStudio 不要在会话之间保留工作区。你可以通过运行 usethis::use_blank_slate(),或者按照图 6.2 中显示的选项进行设置。这样做在短期内会让你有些不便,因为当你重新启动 RStudio 时,它将不再记得你上次运行的代码,你创建的对象或读取的数据集也不会再可用。但这种短期的不便会为你节省长期的痛苦,因为它迫使你把所有重要的步骤都记录在代码中。没有什么比在三个月后发现你只把一个重要计算的结果保存在 Environment 里,而没有把计算过程本身写进代码中,更糟糕的了。

    有一组非常实用的键盘快捷键,可以配合使用,确保你已经把代码中的重要部分保存在编辑器里:

    1. 按 Cmd/Ctrl + Shift + 0/F10 重新启动 R。
    2. 按 Cmd/Ctrl + Shift + S 重新运行当前脚本。

    我们每周都会把这个操作模式一起使用上百次。

      或者,如果你不使用键盘快捷键,也可以依次选择 Session > Restart R,然后高亮并重新运行当前脚本。

      如果你使用的是 RStudio Server,那么默认情况下你的 R 会话永远不会被重启。当你关闭 RStudio Server 的标签页时,感觉上像是在关闭 R,但实际上服务器仍会在后台继续运行。等你下次回来时,你会回到和离开时完全一样的位置。这也就使得定期重启 R 变得更加重要,这样你才能从一个干净的状态开始。

      你的分析存放在哪里?

      你可以通过运行 getwd() 在 R 代码中将其打印出来:

      getwd()
      #> [1] "/Users/hadley/Documents/r4ds"

      在这个 R 会话中,当前工作目录(可以把它看作“家”)位于 hadley 的 Documents 文件夹中,名为 r4ds 的子文件夹里。当你运行这段代码时,它会返回不同的结果,因为你的电脑的目录结构和 Hadley 的不一样!

      作为初学者,把你的工作目录设为主目录、文档目录,或者电脑上的任何其他奇怪目录都没问题。很快你就应该逐渐把项目整理到各自的目录中,并且在处理某个项目时,把 R 的工作目录设置为对应的目录。

      你可以在 R 中设置工作目录,但我们不建议这样做:

      setwd("/path/to/my/CoolProject")

      有一种更好的方法;一种还能让你走上像专家一样管理 R 工作之路的方法。那就是 RStudio 项目。

      RStudio 项目

      将与某个项目相关的所有文件(输入数据、R 脚本、分析结果和图表)都放在同一个目录中,是一种非常明智且常见的做法,因此 RStudio 通过 project 提供了内置支持。让我们为你创建一个 project,供你在阅读本书其余部分时使用。点击 File > New Project,然后按照图 6.3 中所示的步骤进行操作。

      将你的项目命名为 r4ds,并仔细考虑把项目放在哪个子目录中。如果你不把它存放在一个合适的位置,将来就很难找到它!

      一旦这个过程完成,你就会得到一个新 RStudio 项目。请检查你的项目“home”是否就是当前工作目录:

      getwd()
      #> [1] /Users/hadley/Documents/r4ds

      现在在脚本编辑器中输入以下命令,并将文件保存为“diamonds.R”。然后,新建一个名为“data”的文件夹。你可以通过点击 RStudio 中 Files 面板里的“New Folder”按钮来完成。最后,运行整个脚本,这会将一个 PNG 文件和一个 CSV 文件保存到你的项目目录中。

      # 加载 tidyverse 套件
      library(tidyverse)
      # 使用钻石数据集创建 carat 与 price 的六边形分箱图
      ggplot(diamonds, aes(x = carat, y = price)) +
      geom_hex()
      # 将图形保存为 PNG 文件
      ggsave("diamonds.png")
      # 将 diamonds 数据集写入 CSV 文件
      write_csv(diamonds, "data/diamonds.csv")

      退出 RStudio。查看与你的项目相关联的文件夹——注意 .Rproj 文件。双击该文件以重新打开项目。注意你会回到离开时的状态:相同的工作目录和命令历史记录,所有正在处理的文件也仍然打开着。不过,由于你遵循了我们上面的说明,你将拥有一个完全全新的环境,这保证你是从一个干净的起点开始的。

      用你最喜欢的、针对特定操作系统的方式,在电脑中搜索 diamonds.png,你不仅会找到这个 PNG 文件(不奇怪),还会找到创建它的脚本(diamonds.R)。这可是一个巨大的收获!总有一天,你会想重新制作某个图形,或者只是想弄清楚它是从哪里来的。如果你严格地用 R 代码把图形保存到文件中,而不是使用鼠标或剪贴板,那么你就能轻松地重现旧的工作成果!

      相对路径和绝对路径

      一旦进入项目中,你应该始终只使用相对路径,而不要使用绝对路径。它们有什么区别?相对路径是相对于工作目录而言的,也就是项目的 home 目录。比如,Hadley 上面写的 data/diamonds.csv,其实是 /Users/hadley/Documents/r4ds/data/diamonds.csv 的简写。但更重要的是,如果 Mine 在她自己的电脑上运行这段代码,它指向的就会是 /Users/Mine/Documents/r4ds/data/diamonds.csv。这就是相对路径很重要的原因:无论 R 项目文件夹最终放在哪里,它们都能正常工作。

      绝对路径则无论你的工作目录是什么,都会指向同一个位置。它们的形式会因操作系统不同而略有差异。在 Windows 上,它们以驱动器字母开头(例如 C:),或者以两个反斜杠开头(例如 \\servername);而在 Mac/Linux 上,它们以斜杠 / 开头(例如 /users/hadley)。你不应该在脚本中使用绝对路径,因为这会妨碍共享:没有人的目录配置会和你完全一样。

      不同操作系统之间还有一个重要区别:如何分隔路径中的各个部分。Mac 和 Linux 使用正斜杠(例如 data/diamonds.csv),而 Windows 使用反斜杠(例如 data\diamonds.csv)。R 可以处理这两种写法(无论你当前使用的是哪个平台),但不幸的是,反斜杠在 R 中有特殊含义,因此如果你想在路径中输入一个反斜杠,就必须输入两个反斜杠!这会让人很烦,所以我们建议始终使用 Linux/Mac 风格的正斜杠。

      总结

      在这一章中,你已经学会了如何在脚本(文件)和项目(目录)中组织你的 R 代码。和代码风格一样,这一开始可能会让你觉得像是在做琐事。但随着你在多个项目中积累越来越多的代码,你会逐渐体会到,前期做一点组织工作能在后面节省大量时间。

      总之,脚本和项目能为你提供一套可靠的工作流程,并且在未来也会非常有用:

      • 为每个数据分析项目创建一个 RStudio 项目。
      • 将你的脚本(用有意义的名称命名)保存在项目中,编辑它们,并分段或整体运行它们。要经常重启 R,以确保你已经把脚本中的内容都保存好了。
      • 只使用相对路径,不要使用绝对路径。
      • 这样,你所需的一切都会集中在一个地方,并且与其他你正在处理的项目清晰分开。
    1. R语言数据清理

      https://wp.me/p80aHo-2sZ

      Table of Contents

      引言

      在本章中,你将学习一种在 R 里以一致方式组织数据的方法,也就是所谓的“整洁数据(tidy data)”体系。要把数据整理成这种格式,前期需要付出一些工作,但从长远来看这些付出会带来回报。等你拥有了整洁数据,你将花更少的时间把数据从一种表示形式“搓”到另一种表示形式,从而让你能够把更多时间投入到你真正关心的数据问题上。

      在本章中,你将首先学习整洁数据的定义,并用一个简单的玩具数据集来加以说明。然后,我们会深入讲解你在整理数据时将使用的主要工具:pivot(透视/旋转)。pivot 允许你在不改变任何数值的情况下,改变数据的形态。

      整洁数据

      你可以用多种方式来表示同一份底层数据。下面的例子展示了同一份数据以三种不同的方式组织起来。每个数据集都包含四个变量的相同取值:country(国家)、year(年份)、population(人口)以及 TB(结核病)的已记录病例数(number of documented cases of TB),但每个数据集对这些取值的组织方式都不同。

      library(tidyr)
      table1
      #> # A tibble: 6 × 4
      #> country year cases population
      #> <chr> <dbl> <dbl> <dbl>
      #> 1 Afghanistan 1999 745 19987071
      #> 2 Afghanistan 2000 2666 20595360
      #> 3 Brazil 1999 37737 172006362
      #> 4 Brazil 2000 80488 174504898
      #> 5 China 1999 212258 1272915272
      #> 6 China 2000 213766 1280428583
      table2
      #> # A tibble: 12 × 4
      #> country year type count
      #> <chr> <dbl> <chr> <dbl>
      #> 1 Afghanistan 1999 cases 745
      #> 2 Afghanistan 1999 population 19987071
      #> 3 Afghanistan 2000 cases 2666
      #> 4 Afghanistan 2000 population 20595360
      #> 5 Brazil 1999 cases 37737
      #> 6 Brazil 1999 population 172006362
      #> # ℹ 6 more rows
      table3
      #> # A tibble: 6 × 3
      #> country year rate
      #> <chr> <dbl> <chr>
      #> 1 Afghanistan 1999 745/19987071
      #> 2 Afghanistan 2000 2666/20595360
      #> 3 Brazil 1999 37737/172006362
      #> 4 Brazil 2000 80488/174504898
      #> 5 China 1999 212258/1272915272
      #> 6 China 2000 213766/1280428583

      这些都是对同一份底层数据的不同表示方式,但它们并不一样容易使用。其中一种表示(table1)会更容易进行处理,因为它是“整洁的(tidy)”。

      有三条相互关联的规则决定了数据集是否整洁:

      1. 每个变量都是一列;每一列都是一个变量。
      2. 每个观测值都是一行;每一行都是一个观测值。
      3. 每个值都是一个单元格;每个单元格都是一个单一的值。

      为什么要确保你的数据是整洁的(tidy)?主要有两个好处:

      1. 选择一种一致的数据存储方式带来的总体优势。如果你的数据结构是统一的,那么学习那些与之配套的工具会更容易,因为它们背后有一种共同的、相对统一的规律。
      2. 把变量放在列中带来的特定优势:这能让 R 的向量化特性充分发挥。大多数内置的 R 函数都能作用在一组值的向量上。因此,对整洁数据进行转换会感觉特别自然。

      dplyr、ggplot2包,都被设计成可以与整洁数据(tidy data)协同工作。下面是几个小例子,展示你可能会如何使用 table1。

      # Compute rate per 10,000
      table1 |>
      mutate(rate = cases / population * 10000)
      #> # A tibble: 6 × 5
      #> country year cases population rate
      #> <chr> <dbl> <dbl> <dbl> <dbl>
      #> 1 Afghanistan 1999 745 19987071 0.373
      #> 2 Afghanistan 2000 2666 20595360 1.29
      #> 3 Brazil 1999 37737 172006362 2.19
      #> 4 Brazil 2000 80488 174504898 4.61
      #> 5 China 1999 212258 1272915272 1.67
      #> 6 China 2000 213766 1280428583 1.67
      # Compute total cases per year
      table1 |>
      group_by(year) |>
      summarize(total_cases = sum(cases))
      #> # A tibble: 2 × 2
      #> year total_cases
      #> <dbl> <dbl>
      #> 1 1999 250740
      #> 2 2000 296920
      # Visualize changes over time
      ggplot(table1, aes(x = year, y = cases)) +
      geom_line(aes(group = country), color = "grey50") +
      geom_point(aes(color = country, shape = country)) +
      scale_x_continuous(breaks = c(1999, 2000)) # x-axis breaks at 1999 and 2000

      拉长数据

      整洁数据(tidy data)的原则看起来似乎如此显而易见,以至于你可能会想:你是不是再也不会遇到不整洁的数据集了。然而,不幸的是,大多数真实数据都是不整洁的。有两个主要原因:

      数据往往是为实现某种目标(而不是为了分析)而进行组织的。例如,数据被设计成便于数据录入而不是便于分析,这是很常见的。

      大多数人并不熟悉整洁数据的原则,而且除非你花很多时间亲自处理数据,否则很难自己推导出这些原则。

      这意味着,大多数真实分析都至少需要做一点整理(tidying)。你将从弄清楚底层的变量和观测值是什么开始。有时这很容易;但有时你需要去查阅最初生成这些数据的人。接下来,你会把数据转换成整洁的形式:把变量放在列中,把观测值放在行中。

      tidyr 提供了两个用于数据透视(pivot)的函数:pivot_longer() 和 pivot_wider()。我们会先从最常见的情况 pivot_longer() 开始。下面我们进入一些例子。

      列名中的数据

      广告牌数据集记录了2000年歌曲的广告牌排名:

      billboard
      #> # A tibble: 317 × 79
      #> artist track date.entered wk1 wk2 wk3 wk4 wk5
      #> <chr> <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
      #> 1 2 Pac Baby Don't Cry (Ke… 2000-02-26 87 82 72 77 87
      #> 2 2Ge+her The Hardest Part O… 2000-09-02 91 87 92 NA NA
      #> 3 3 Doors Down Kryptonite 2000-04-08 81 70 68 67 66
      #> 4 3 Doors Down Loser 2000-10-21 76 76 72 69 67
      #> 5 504 Boyz Wobble Wobble 2000-04-15 57 34 25 17 17
      #> 6 98^0 Give Me Just One N… 2000-08-19 51 39 34 26 26
      #> # ℹ 311 more rows
      #> # ℹ 71 more variables: wk6 <dbl>, wk7 <dbl>, wk8 <dbl>, wk9 <dbl>, …

      在该数据集中,每个观测值都是一首歌曲。前三列(artist,track 和 date.entered)是描述这首歌曲的变量。接着我们有 76 列(wk1-wk76),用于描述该歌曲在每周的排名。这里,列名是一个变量(周次),单元格的值是另一个变量(排名)。

      为了整理这份数据,我们将使用 pivot_longer():

      billboard |>
      pivot_longer(
      cols = starts_with("wk"), # 选择所有列名以“wk”开头的列(例如 wk1, wk2, ... wk76)
      names_to = "week", # 把这些列的列名(wk1, wk2, ...)整理到新变量 week 中
      values_to = "rank" # 把这些列对应的单元格值(排名数字)整理到新变量 rank 中
      )
      #> # A tibble: 24,092 × 5
      #> artist track date.entered week rank
      #> <chr> <chr> <date> <chr> <dbl>
      #> 1 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk1 87
      #> 2 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk2 82
      #> 3 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk3 72
      #> 4 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk4 77
      #> 5 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk5 87
      #> 6 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk6 94
      #> 7 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk7 99
      #> 8 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk8 NA
      #> 9 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk9 NA
      #> 10 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk10 NA
      #> # ℹ 24,082 more rows

      在数据之后,有三个关键参数:

      • cols 指定需要进行 pivot 的列,也就是哪些列不是变量。该参数使用与 select() 相同的语法,因此这里我们可以使用 !c(artist, track, date.entered) 或 starts_with("wk")
      • names_to 指定将存储在列名中的变量,我们把该变量命名为 week
      • values_to 指定将存储在单元格值中的变量,我们把该变量命名为 rank

      注意:在代码中之所以给 "week" 和 "rank" 加上引号,是因为它们是我们要新创建的变量;在运行 pivot_longer() 调用时,这些变量在原始数据中还不存在。

      现在让我们把注意力转向整理后的、更长的数据框。如果一首歌进入了前 100 名的时间少于 76 周,会发生什么?例如,2 Pac 的 “Baby Don’t Cry”。上面的输出表明它只在前 100 名中出现了 7 周,而其余周都会填入缺失值。这些 NA 并不真正表示未知的观测值;它们是由数据集的结构“被迫”存在的,所以我们可以通过设置 values_drop_na = TRUE 让 pivot_longer() 把它们去掉:

      billboard |>
      pivot_longer(
      cols = starts_with("wk"), # 选择所有列名以“wk”开头的列(即 wk1-wk76),把它们从宽表转换为长表
      names_to = "week", # 将这些列名(wk1、wk2、...)存到新变量 week 中
      values_to = "rank", # 将每个单元格的值(排名数字)存到新变量 rank 中
      values_drop_na = TRUE # 删除生成结果中的缺失值(NA);避免为不足 76 周的歌曲补出的无效观测
      )
      #> # A tibble: 5,307 × 5
      #> artist track date.entered week rank
      #> <chr> <chr> <date> <chr> <dbl>
      #> 1 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk1 87
      #> 2 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk2 82
      #> 3 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk3 72
      #> 4 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk4 77
      #> 5 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk5 87
      #> 6 2 Pac Baby Don't Cry (Keep... 2000-02-26 wk6 94
      #> # ℹ 5,301 more rows

      行数现在大幅减少,说明丢掉了很多包含 NA 的行。

      你可能还会好奇:如果一首歌在前 100 名中超过了 76 周会怎样?仅凭这份数据我们无法得知,但你可能会猜想,数据集中会新增更多列,如 wk77、wk78、……。

      这份数据现在已经整理得很整齐了,不过我们可以通过使用 mutate() 和 readr::parse_number(),把将 week 的值从字符字符串转换为数字,从而让后续计算更容易一些。parse_number() 是一个很方便的函数:它会从字符串中提取第一个数字,并忽略字符串中的其他所有文本。

      library(readr)
      billboard_longer <- billboard |>
      pivot_longer(
      cols = starts_with("wk"), # 选择所有以“wk”开头的列(wk1-wk76)进行透视展开
      names_to = "week", # 把原列名(wk1、wk2、...)存入新变量 week
      values_to = "rank", # 把原单元格值(排名)存入新变量 rank
      values_drop_na = TRUE # 删除由“未覆盖到所有周”产生的 NA,避免无效观测
      ) |>
      mutate(
      week = parse_number(week) # 将 week 由字符串(例如 "wk1")转换为数字(例如 1)
      )
      billboard_longer # 输出整理后的长表数据框
      #> # A tibble: 5,307 × 5
      #> artist track date.entered week rank
      #> <chr> <chr> <date> <dbl> <dbl>
      #> 1 2 Pac Baby Don't Cry (Keep... 2000-02-26 1 87
      #> 2 2 Pac Baby Don't Cry (Keep... 2000-02-26 2 82
      #> 3 2 Pac Baby Don't Cry (Keep... 2000-02-26 3 72
      #> 4 2 Pac Baby Don't Cry (Keep... 2000-02-26 4 77
      #> 5 2 Pac Baby Don't Cry (Keep... 2000-02-26 5 87
      #> 6 2 Pac Baby Don't Cry (Keep... 2000-02-26 6 94
      #> # ℹ 5,301 more rows

      现在我们已经把所有的周次编号放在同一个变量里、把所有的排名数值放在另一个变量里了,因此我们很适合用可视化来展示歌曲排名如何随时间变化。下面给出了代码,运行结果如图 5.2 所示。我们可以看到,很少有歌曲会在前 100 名中停留超过 20 周。

      billboard_longer |>
      ggplot(aes(x = week, y = rank, group = track)) + # 指定绘图数据与美学映射:横轴 week,纵轴 rank;按 track 分组以便连线
      geom_line(alpha = 0.25) + # 用折线图连接同一首歌(同一 track)在不同周的排名;alpha=0.25 设置透明度,便于观察重叠
      scale_y_reverse() # 反转 y 轴方向:排名 1(最好)显示在最上方

      透视(pivot)是如何工作的?

      现在你已经看到了我们如何使用透视(pivoting)来重塑数据,让我们花一点时间来建立直观理解:pivoting 到底会如何改变数据。我们先从一个非常简单的数据集开始,这样更容易看清发生了什么。假设我们有三位患者,id 分别是 A、B 和 C,并且每位患者都有两次血压测量。我们将使用 tribble() 来创建这些数据——tribble() 是一个很方便的函数,可以手动构造小的 tibble:

      df <- tribble(
      ~id, ~bp1, ~bp2,
      "A", 100, 120,
      "B", 140, 115,
      "C", 120, 125
      )

      我们希望新的数据集包含三个变量:id(已经存在)、measurement(列名)以及 value(单元格里的值)。为实现这一点,我们需要把 df 透视为更长的格式(pivot df longer):

      df |>
      pivot_longer(
      cols = bp1:bp2, # 选择要透视展开的列范围(bp1 到 bp2),这些列名会被转成新的变量 measurement
      names_to = "measurement",# 将被透视的列名(如 bp1、bp2)放入新列 measurement
      values_to = "value" # 将对应单元格的数值放入新列 value
      )
      #> # A tibble: 6 × 3
      #> id measurement value
      #> <chr> <chr> <dbl>
      #> 1 A bp1 100
      #> 2 A bp2 120
      #> 3 B bp1 140
      #> 4 B bp2 115
      #> 5 C bp1 120
      #> 6 C bp2 125

      重塑(reshaping)是如何工作的?如果我们按列来想,会更容易理解。正如图 5.3 所示:在原始数据集中已经是变量的某一列(id)里的值,需要为每一个被透视展开(pivoted)的列各重复一次。

      列名会变成一个新变量中的取值,而该新变量的名称由 names_to 定义,如图 5.4 所示。它们需要针对原始数据集中的每一行重复一次。

      单元格里的值也会变成新变量中的值,其名称由 values_to 指定。它们会按行逐行“解开”(unwound)。图 5.5 展示了这个过程。

      列名中有很多变量

      当列名中塞进了多条信息,而且你希望把这些信息分别存到若干个新的变量中时,就会遇到更具挑战性的情况。比如,来看一下 who2 数据集——也就是你在上面看到的 table1 及其相关内容的来源:

      who2
      #> # A tibble: 7,240 × 58
      #> country year sp_m_014 sp_m_1524 sp_m_2534 sp_m_3544 sp_m_4554
      #> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
      #> 1 Afghanistan 1980 NA NA NA NA NA
      #> 2 Afghanistan 1981 NA NA NA NA NA
      #> 3 Afghanistan 1982 NA NA NA NA NA
      #> 4 Afghanistan 1983 NA NA NA NA NA
      #> 5 Afghanistan 1984 NA NA NA NA NA
      #> 6 Afghanistan 1985 NA NA NA NA NA
      #> # ℹ 7,234 more rows
      #> # ℹ 51 more variables: sp_m_5564 <dbl>, sp_m_65 <dbl>, sp_f_014 <dbl>, …

      该数据集由世界卫生组织收集,记录了关于结核病(tuberculosis)诊断的信息。已有两列变量,且很容易理解:country 和 year。它们之后还有 56 列,例如 sp_m_014、ep_m_4554 和 rel_m_3544。只要盯着这些列名看够久,你就会发现其中存在规律。每个列名由三部分组成,它们之间用 _ 分隔:第一部分 sp/rel/ep 描述用于诊断的方法;第二部分 m/f 是性别(在该数据集中以二元变量编码);第三部分 014/1524/2534/3544/4554/5564/65 表示年龄范围(例如 014 代表 0-14)。

      所以在这个例子中,who2 里记录了六条信息:国家和年份(已经是列);诊断方法、性别类别和年龄段类别(包含在其他列名中);以及该类别下患者的数量(单元格的值)。为了把这六条信息组织到六个独立的列中,我们使用 pivot_longer():对 names_to 使用一个列名向量,并对 instructors 使用用于拆分原始变量名为若干片段的设置,同时还需要一个用于指定 values_to 的列名:

      who2 |>
      pivot_longer(
      cols = !(country:year), # 选择除 country 和 year 之外的所有列作为要被“转成长表”的列
      names_to = c("diagnosis", "gender", "age"), # 将列名按分隔符拆成三部分,分别存入这三个新变量
      names_sep = "_", # 三部分之间用 "_" 分隔(对应列名格式:diagnosis_gender_age)
      values_to = "count" # 被透视后的数值放到新列 count 中
      )
      #> # A tibble: 405,440 × 6
      #> country year diagnosis gender age count
      #> <chr> <dbl> <chr> <chr> <chr> <dbl>
      #> 1 Afghanistan 1980 sp m 014 NA
      #> 2 Afghanistan 1980 sp m 1524 NA
      #> 3 Afghanistan 1980 sp m 2534 NA
      #> 4 Afghanistan 1980 sp m 3544 NA
      #> 5 Afghanistan 1980 sp m 4554 NA
      #> 6 Afghanistan 1980 sp m 5564 NA
      #> # ℹ 405,434 more rows

      从概念上来说,这只是对你之前已经见过的更简单情况的一种小幅变化。图 5.6 展示了这个基本思路:现在,不再是把列名透视(pivot)到一个单独的列中,而是透视到多个列中。你可以把这个过程想象为两步完成(先透视,再分离),但在底层它会在一步内完成,因为这样更快。

      数据与变量名出现在列标题中

      复杂度更进一步时,列名中会同时包含变量值和变量名。例如,来看家庭数据集:

      library(tidyr)
      household
      #> # A tibble: 5 × 5
      #> family dob_child1 dob_child2 name_child1 name_child2
      #> <int> <date> <date> <chr> <chr>
      #> 1 1 1998-11-26 2000-01-29 Susan Jose
      #> 2 2 1996-06-22 NA Mark <NA>
      #> 3 3 2002-07-11 2004-04-05 Sam Seth
      #> 4 4 2004-10-10 2009-08-27 Craig Khai
      #> 5 5 2000-12-05 2005-02-28 Parker Gracie

      这个数据集包含关于五个家庭的数据,以及最多两个孩子的姓名和出生日期。该数据集中的新挑战在于,列名同时包含两个变量(dob、name)的名称,以及另一个变量(child)的取值(其值为 1 或 2)。为了解决这个问题,我们同样需要向 names_to 提供一个向量,但这次使用特殊的 “.value” 哨兵;这并不是某个变量的名称,而是一个独特的值,用来告诉 pivot_longer() 采取不同的处理方式。它会覆盖通常的 values_to 参数,把透视后列名的第一个组成部分作为输出中的变量名。

      household |> # 以 household 数据框为起点,接着把结果传给 pivot_longer()
      pivot_longer( # 将数据从宽格式转换为长格式
      cols = !family, # 除了 family 列以外,其余列都作为要整理的列
      names_to = c(".value", "child"), # 将列名拆成两部分:前半部分作为输出列名,后半部分存入 child
      names_sep = "_", # 按下划线 "_" 拆分列名
      values_drop_na = TRUE # 删除转换后产生的缺失值行
      ) # 结束 pivot_longer() 调用
      #> # A tibble: 9 × 4
      #> family child dob name
      #> <int> <chr> <date> <chr>
      #> 1 1 child1 1998-11-26 Susan
      #> 2 1 child2 2000-01-29 Jose
      #> 3 2 child1 1996-06-22 Mark
      #> 4 3 child1 2002-07-11 Sam
      #> 5 3 child2 2004-04-05 Seth
      #> 6 4 child1 2004-10-10 Craig
      #> # ℹ 3 more rows

      我们再次使用 values_drop_na = TRUE,因为输入数据的形状会强制创建显式的缺失变量(例如,对于只有一个孩子的家庭)。

      图 5.7 用一个更简单的例子说明了这一基本思路。当你在 names_to 中使用 “.value” 时,输入中的列名会同时参与输出中的值和变量名。

      数据变宽

      到目前为止,我们一直使用 pivot_longer() 来解决一种常见问题:值被放到了列名中。接下来我们转向 pivot_wider(),它通过增加列数、减少行数来让数据集变宽,并且在一条观察被分散在多行中时很有帮助。这种情况在实际中似乎没那么常见,但在处理政府数据时似乎经常出现。

      我们先来看 cms_patient_experience,这是来自医疗保险和医疗补助服务中心(Centers of Medicare and Medicaid Services)的一个数据集,用于收集患者体验方面的数据:

      cms_patient_experience
      #> # A tibble: 500 × 5
      #> org_pac_id org_nm measure_cd measure_title prf_rate
      #> <chr> <chr> <chr> <chr> <dbl>
      #> 1 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_1 CAHPS for MIPS… 63
      #> 2 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_2 CAHPS for MIPS… 87
      #> 3 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_3 CAHPS for MIPS… 86
      #> 4 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_5 CAHPS for MIPS… 57
      #> 5 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_8 CAHPS for MIPS… 85
      #> 6 0446157747 USC CARE MEDICAL GROUP INC CAHPS_GRP_12 CAHPS for MIPS… 24
      #> # ℹ 494 more rows

      研究的核心单位是一个组织,但每个组织分布在六行中,调查中对组织所做的每一次测量各占一行。我们可以通过使用 distinct() 来查看 measure_cd 和 measure_title 的完整取值集合:

      cms_patient_experience |> # 以 cms_patient_experience 数据框为起点
      distinct(measure_cd, measure_title) # 提取 measure_cd 和 measure_title 的唯一组合
      #> # A tibble: 6 × 2
      #> measure_cd measure_title
      #> <chr> <chr>
      #> 1 CAHPS_GRP_1 CAHPS for MIPS SSM: Getting Timely Care, Appointments, and In…
      #> 2 CAHPS_GRP_2 CAHPS for MIPS SSM: How Well Providers Communicate
      #> 3 CAHPS_GRP_3 CAHPS for MIPS SSM: Patient's Rating of Provider
      #> 4 CAHPS_GRP_5 CAHPS for MIPS SSM: Health Promotion and Education
      #> 5 CAHPS_GRP_8 CAHPS for MIPS SSM: Courteous and Helpful Office Staff
      #> 6 CAHPS_GRP_12 CAHPS for MIPS SSM: Stewardship of Patient Resources

      这两个列名都不太适合作为变量名:measure_cd 不能说明变量含义,而 measure_title 是一个包含空格的长句子。现在我们先用 measure_cd 作为新列名的来源,但在真正的分析中,你可能会想自己创建既简短又有意义的变量名。

      pivot_wider() 的接口与 pivot_longer() 正好相反:我们不是选择新的列名,而是需要提供用于定义值的现有列(values_from)以及用于定义列名的列(names_from):

      cms_patient_experience |> # 以 cms_patient_experience 数据框为起点
      pivot_wider( # 将数据从长格式转换为宽格式
      names_from = measure_cd, # 使用 measure_cd 的值作为新列名
      values_from = prf_rate # 使用 prf_rate 的值填充新生成的列
      ) # 结束 pivot_wider() 调用
      #> # A tibble: 500 × 9
      #> org_pac_id org_nm measure_title CAHPS_GRP_1 CAHPS_GRP_2
      #> <chr> <chr> <chr> <dbl> <dbl>
      #> 1 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… 63 NA
      #> 2 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… NA 87
      #> 3 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… NA NA
      #> 4 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… NA NA
      #> 5 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… NA NA
      #> 6 0446157747 USC CARE MEDICAL GROUP … CAHPS for MIPS… NA NA
      #> # ℹ 494 more rows
      #> # ℹ 4 more variables: CAHPS_GRP_3 <dbl>, CAHPS_GRP_5 <dbl>, …

      输出结果看起来不太对;我们似乎仍然为每个组织保留了多行。这是因为,我们还需要告诉 pivot_wider() 哪一列或哪些列的值能够唯一标识每一行;在这个例子中,这些列是以 “org” 开头的变量:

      cms_patient_experience |> # 以 cms_patient_experience 数据框为起点
      pivot_wider( # 将数据从长格式转换为宽格式
      id_cols = starts_with("org"), # 使用以 "org" 开头的列作为标识列
      names_from = measure_cd, # 使用 measure_cd 的值作为新列名
      values_from = prf_rate # 使用 prf_rate 的值填充新生成的列
      ) # 结束 pivot_wider() 调用
      #> # A tibble: 95 × 8
      #> org_pac_id org_nm CAHPS_GRP_1 CAHPS_GRP_2 CAHPS_GRP_3 CAHPS_GRP_5
      #> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
      #> 1 0446157747 USC CARE MEDICA… 63 87 86 57
      #> 2 0446162697 ASSOCIATION OF … 59 85 83 63
      #> 3 0547164295 BEAVER MEDICAL … 49 NA 75 44
      #> 4 0749333730 CAPE PHYSICIANS… 67 84 85 65
      #> 5 0840104360 ALLIANCE PHYSIC… 66 87 87 64
      #> 6 0840109864 REX HOSPITAL INC 73 87 84 67
      #> # ℹ 89 more rows
      #> # ℹ 2 more variables: CAHPS_GRP_8 <dbl>, CAHPS_GRP_12 <dbl>

      这就得到了我们想要的输出。

      pivot_wider() 是如何工作的?

      为了理解 pivot_wider() 的工作方式,我们再次从一个非常简单的数据集开始。这次我们有两位患者,ID 分别为 A 和 B,患者 A 有三次血压测量,患者 B 有两次:

      df <- tribble(
      ~id, ~measurement, ~value,
      "A", "bp1", 100,
      "B", "bp1", 140,
      "B", "bp2", 115,
      "A", "bp2", 120,
      "A", "bp3", 105
      )

      我们将从 value 列中取值,并从 measurement 列中取列名:

      df |>
      pivot_wider(
      names_from = measurement,
      values_from = value
      )
      #> # A tibble: 2 × 4
      #> id bp1 bp2 bp3
      #> <chr> <dbl> <dbl> <dbl>
      #> 1 A 100 120 105
      #> 2 B 140 115 NA

      要开始这个过程,pivot_wider() 需要先确定哪些内容放在行里、哪些内容放在列里。新的列名将取自 measurement 的唯一值。

      df |>
      distinct(measurement) |>
      pull()
      #> [1] "bp1" "bp2" "bp3"

      默认情况下,输出中的行由所有不用于生成新名称或新值的变量决定。这些变量称为 id_cols。这里只有一列,但一般来说可以有任意数量。

      df |>
      select(!measurement & !value) |>
      distinct()
      #> # A tibble: 2 × 1
      #> id
      #> <chr>
      #> 1 A
      #> 2 B

      pivot_wider() 然后将这些结果组合起来,生成一个空的数据框:

      df |>
      select(!measurement & !value) |>
      distinct() |>
      mutate(x = NA, y = NA, z = NA)
      #> # A tibble: 2 × 4
      #> id x y z
      #> <chr> <lgl> <lgl> <lgl>
      #> 1 A NA NA NA
      #> 2 B NA NA NA

      然后,它会使用输入中的数据填充所有缺失的值。在这种情况下,输出中的每个单元格都不一定在输入中有对应值,因为患者 B 没有第三次血压测量,所以那个单元格仍然是缺失的。

      你可能还会想,如果输入中有多行对应输出中的一个单元格,会发生什么。下面的例子中有两行对应于 id “A” 和 measurement “bp1”:

      df <- tribble(
      ~id, ~measurement, ~value,
      "A", "bp1", 100,
      "A", "bp1", 102, # 这行重复了
      "A", "bp2", 120,
      "B", "bp1", 140,
      "B", "bp2", 115
      )

      如果我们尝试对此进行透视,就会得到一个包含列表列的输出

      df2 <- df |>
      pivot_wider(
      names_from = measurement,
      values_from = value
      )

      由于你还不知道如何处理这类数据,你需要按照警告中的提示来找出问题所在:

      # 按 id 和 measurement 分组,统计每组的行数
      df |>
      group_by(id, measurement) |>
      summarize(n = n(), .groups = "drop") |>
      # 只保留重复出现的组合(n > 1)
      filter(n > 1)
      #> # A tibble: 1 × 3
      #> id measurement n
      #> <chr> <chr> <int>
      #> 1 A bp1 2

      接下来就要由你来弄清楚数据出了什么问题,并修复底层损坏,或者运用分组和汇总技巧,确保每一种行值和列值的组合都只有一行。

      总结

      在这一章中,你学习了整洁数据:即变量放在列中、观测放在行中的数据。整洁数据能让你更轻松地工作,因为它是一种大多数函数都能理解的一致结构;而主要的挑战,是把你接收到的任何结构的数据转换成整洁格式。为此,你学习了 pivot_longer() 和 pivot_wider(),它们可以帮助你整理许多不整洁的数据集。我们这里展示的示例,是 vignette("pivot", package = "tidyr") 中示例的一部分,所以如果你遇到本章无法帮助你解决的问题,那篇 vignette 是下一个很好的尝试地点。

      另一个挑战是,对于某个数据集来说,长格式或宽格式版本有时都无法被明确称为“整洁”的那个。这在一定程度上反映了我们对整洁数据的定义:我们说整洁数据是每一列对应一个变量,但我们实际上并没有定义什么是变量(而且这件事出乎意料地难以界定)。实用一点地说,把变量定义为“最有利于你的分析的东西”完全没问题。所以如果你在弄清楚如何进行某个计算时遇到困难,不妨尝试重新组织数据;不要害怕在需要时先弄得不整洁、再转换、然后再整理回整洁!

      既然你已经开始编写大量的 R 代码,现在是时候进一步学习如何把代码组织到文件和目录中了。

    2. R语言堆叠面积图

      library(gcookbook) # Load gcookbook for the uspopage data set
      library(ggplot2)
      # 使用 ggplot2:以 uspopage 为数据源,设置美学映射
      ggplot(
      uspopage,
      aes(
      x = Year, # x轴:年份
      y = Thousands, # y轴:Thousands
      fill = AgeGroup, # 用 AgeGroup 决定填充颜色
      order = dplyr::desc(AgeGroup) # 设置堆叠/绘制顺序:AgeGroup 倒序
      )
      ) +
      # 画堆叠面积图
      # colour = NA:不画面积边界线(避免边界干扰)
      # alpha = .4:填充透明度 0.4
      geom_area(colour = NA, alpha = .4) +
      # 使用 RColorBrewer 的 Blues 配色方案为 fill 上色
      scale_fill_brewer(palette = "Blues") +
      # 在“堆叠结构一致”的基础上叠加折线
      # position = "stack":线的位置按堆叠方式计算
      # size = .2:线宽很细
      geom_line(position = "stack", size = .2)