13  Numbers

13.1 Introduction

数值向量(Numeric vectors)是数据科学的支柱,您已在本节前文中多次使用过它们。 现在需要系统性地考察如何在 R 语言中操作数值向量,确保您为今后处理涉及数值向量的任何问题做好充分准备。

我们将首先介绍几种将字符串转换为数字的工具,然后更详细地讲解 count() 函数。 接下来将深入探讨与 mutate() 配合使用的多种数值变换方法,包括那些可应用于其他类型向量但常用于数值向量的通用转换。 最后将介绍与 summarize() 搭配使用的汇总函数,并演示如何将其与 mutate() 结合使用。

13.1.1 Prerequisites

本章主要使用 base R 中的函数,这些函数无需加载任何包即可使用。 但我们仍然需要 tidyverse,因为我们会在 mutate()filter() 等 tidyverse 函数内部使用这些基础 R 函数。 与上一章类似,我们将使用 nycflights13 中的真实示例,以及通过 c()tribble() 创建的模拟示例。

13.2 Making numbers

在大多数情况下,您获得的数字已以 R 的数值类型(整数或双精度)记录。 然而在某些情况下,您会遇到字符串形式的数字,这可能是因为您通过从列标题转换而创建了它们,或者是因为数据导入过程中出现了问题。

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

x <- c("1.2", "5.6", "1e3")
parse_double(x)
#> [1]    1.2    5.6 1000.0

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

x <- c("$1,234", "USD 3,513", "59%")
parse_number(x)
#> [1] 1234 3513   59

13.3 Counts

令人惊讶的是,仅通过计数和少量基础运算就能完成大量数据科学工作,因此 dplyr 力求通过 count() 使计数尽可能简单。 该函数非常适用于快速探索和分析过程中的检查:

flights |> count(dest)
#> # A tibble: 105 × 2
#>   dest      n
#>   <chr> <int>
#> 1 ABQ     254
#> 2 ACK     265
#> 3 ALB     439
#> 4 ANC       8
#> 5 ATL   17215
#> 6 AUS    2439
#> # ℹ 99 more rows

(尽管 Chapter 4 提出了相关建议,但我们通常将 count() 单独放在一行,因为它通常在控制台中使用,用于快速检查计算是否符合预期。)

如果想查看最常见的值,可添加 sort = TRUE 参数:

flights |> count(dest, sort = TRUE)
#> # A tibble: 105 × 2
#>   dest      n
#>   <chr> <int>
#> 1 ORD   17283
#> 2 ATL   17215
#> 3 LAX   16174
#> 4 BOS   15508
#> 5 MCO   14082
#> 6 CLT   14064
#> # ℹ 99 more rows

请记住,若要查看所有值,可以使用 |> View()|> print(n = Inf)

您可以通过 group_by()summarize()n() 手动执行相同的计算。这 种方式很有用,因为它允许您同时计算其他摘要信息:

flights |> 
  group_by(dest) |> 
  summarize(
    n = n(),
    delay = mean(arr_delay, na.rm = TRUE)
  )
#> # A tibble: 105 × 3
#>   dest      n delay
#>   <chr> <int> <dbl>
#> 1 ABQ     254  4.38
#> 2 ACK     265  4.85
#> 3 ALB     439 14.4 
#> 4 ANC       8 -2.5 
#> 5 ATL   17215 11.3 
#> 6 AUS    2439  6.02
#> # ℹ 99 more rows

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

n()
#> Error in `n()`:
#> ! Must only be used inside data-masking verbs like `mutate()`,
#>   `filter()`, and `group_by()`.

n()count() 有几个变体可能会对您有所帮助:

  • n_distinct(x) 用于计算一个或多个变量的不同(唯一)值的数量。 例如,我们可以找出由最多航空公司(carriers)服务的目的地(destinations):

    flights |> 
      group_by(dest) |> 
      summarize(carriers = n_distinct(carrier)) |> 
      arrange(desc(carriers))
    #> # A tibble: 105 × 2
    #>   dest  carriers
    #>   <chr>    <int>
    #> 1 ATL          7
    #> 2 BOS          7
    #> 3 CLT          7
    #> 4 ORD          7
    #> 5 TPA          7
    #> 6 AUS          6
    #> # ℹ 99 more rows
  • 加权计数是一种求和运算。 例如,您可以”统计”每架飞机飞行的英里数:

    flights |> 
      group_by(tailnum) |> 
      summarize(miles = sum(distance))
    #> # A tibble: 4,044 × 2
    #>   tailnum  miles
    #>   <chr>    <dbl>
    #> 1 D942DN    3418
    #> 2 N0EGMQ  250866
    #> 3 N10156  115966
    #> 4 N102UW   25722
    #> 5 N103US   24619
    #> 6 N104UW   25157
    #> # ℹ 4,038 more rows

    加权计数是常见问题,因此 count() 设有 wt 参数来实现相同功能:

    flights |> count(tailnum, wt = distance)
  • 您可以通过结合 sum()is.na() 来统计缺失值。 在 flights 数据集中,这代表被取消的航班:

    flights |> 
      group_by(dest) |> 
      summarize(n_cancelled = sum(is.na(dep_time))) 
    #> # A tibble: 105 × 2
    #>   dest  n_cancelled
    #>   <chr>       <int>
    #> 1 ABQ             0
    #> 2 ACK             0
    #> 3 ALB            20
    #> 4 ANC             0
    #> 5 ATL           317
    #> 6 AUS            21
    #> # ℹ 99 more rows

13.3.1 Exercises

  1. 你如何使用 count() 来统计给定变量存在缺失值的行数?
  2. 将以下 count() 调用展开为使用 group_by()summarize()arrange()
    1. flights |> count(dest, sort = TRUE)

    2. flights |> count(tailnum, wt = distance)

13.4 Numeric transformations

转换函数与 mutate() 配合良好,因为它们的输出长度与输入相同。 绝大多数转换函数已内置在 base R 中。 逐一列举它们并不现实,因此本节将展示最实用的函数。 例如,虽然 R 提供了您可能需要的所有三角函数,但我们在此不列举它们,因为数据科学中很少需要这些函数。

13.4.1 Arithmetic and recycling rules

我们在 Chapter 2 介绍了算术运算的基础(+, -, *, /, ^)并多次使用过它们。 这些函数无需过多解释,因为它们的功能与你在小学所学的一致。 但我们需要简要讨论回收规则(recycling rules),该规则决定了当左右两侧长度不同时的处理方式。 这对于像 flights |> mutate(air_time = air_time / 60) 这样的操作很重要,因为 / 左侧有 336,776 个数字而右侧只有一个。

R 通过回收(recycling)或重复短向量来处理长度不匹配的情况。 如果在数据框外创建一些向量,我们可以更直观地看到这个机制:

x <- c(1, 2, 10, 20)
x / 5
#> [1] 0.2 0.4 2.0 4.0
# is shorthand for
x / c(5, 5, 5, 5)
#> [1] 0.2 0.4 2.0 4.0

通常,您只会想要回收单个数字(即长度为 1 的向量),但 R 会回收任何较短的向量。 如果较长向量的长度不是较短向量的整数倍,R 通常会(但不总是)给出警告:

x * c(1, 2)
#> [1]  1  4 10 40
x * c(1, 2, 3)
#> Warning in x * c(1, 2, 3): longer object length is not a multiple of shorter
#> object length
#> [1]  1  4 30 20

这些回收规则同样适用于逻辑比较(==, <, <=, >, >=, !=),如果您意外使用 == 而不是 %in%,且数据框的行数恰好不符合倍数关系,可能会导致出人意料的结果。 例如以下试图查找一月和二月所有航班的代码:

flights |> 
  filter(month == c(1, 2))
#> # A tibble: 25,977 × 19
#>    year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
#>   <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
#> 1  2013     1     1      517            515         2      830            819
#> 2  2013     1     1      542            540         2      923            850
#> 3  2013     1     1      554            600        -6      812            837
#> 4  2013     1     1      555            600        -5      913            854
#> 5  2013     1     1      557            600        -3      838            846
#> 6  2013     1     1      558            600        -2      849            851
#> # ℹ 25,971 more rows
#> # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>, …

这段代码运行时不会报错,但返回的结果并非您所期望的。 由于回收规则的存在,系统会查找奇数行中一月起飞的航班和偶数行中二月起飞的航班。 不幸的是,由于 flights 数据集的行数恰好是偶数,这个操作不会触发任何警告。

为防止这类静默失败,大多数 tidyverse 函数采用更严格的回收机制——仅允许回收单个值。 但遗憾的是,这种方法在此处或许多其他场景中并不奏效,因为关键计算是由 base R 函数 == 而非 filter() 执行的。

13.4.2 Minimum and maximum

算术函数可处理成对变量。 两个密切相关的函数是 pmin()pmax(),当输入两个或多个变量时,它们会返回每行中的最小值或最大值:

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)
  )
#> # A tibble: 3 × 4
#>       x     y   min   max
#>   <dbl> <dbl> <dbl> <dbl>
#> 1     1     3     1     3
#> 2     5     2     2     5
#> 3     7    NA     7     7

需要注意的是,这些函数不同于汇总函数 min()max() 接收多个观测值但返回单个值。 当所有最小值或所有最大值都相同时,就能发现使用了错误的形式:

df |> 
  mutate(
    min = min(x, y, na.rm = TRUE),
    max = max(x, y, na.rm = TRUE)
  )
#> # A tibble: 3 × 4
#>       x     y   min   max
#>   <dbl> <dbl> <dbl> <dbl>
#> 1     1     3     1     7
#> 2     5     2     1     7
#> 3     7    NA     1     7

13.4.3 Modular arithmetic

模数运算是学习小数点之前所接触的数学运算的技术名称,即产生整数和余数的除法。 在 R 中,%/% 执行整数除法,%% 计算余数:

1:10 %/% 3
#>  [1] 0 0 1 1 1 2 2 2 3 3
1:10 %% 3
#>  [1] 1 2 0 1 2 0 1 2 0 1

模运算在处理 flights 数据集时非常方便,因为我们可以用它来将 sched_dep_time 变量分解为 hourminute

flights |> 
  mutate(
    hour = sched_dep_time %/% 100,
    minute = sched_dep_time %% 100,
    .keep = "used"
  )
#> # A tibble: 336,776 × 3
#>   sched_dep_time  hour minute
#>            <int> <dbl>  <dbl>
#> 1            515     5     15
#> 2            529     5     29
#> 3            540     5     40
#> 4            545     5     45
#> 5            600     6      0
#> 6            558     5     58
#> # ℹ 336,770 more rows

我们可以将此与 Section 12.4 中的 mean(is.na(x)) 技巧结合,来观察航班取消比例在一天内的变化情况。 结果如 Figure 13.1 所示。

flights |> 
  group_by(hour = sched_dep_time %/% 100) |> 
  summarize(prop_cancelled = mean(is.na(dep_time)), n = n()) |> 
  filter(hour > 1) |> 
  ggplot(aes(x = hour, y = prop_cancelled)) +
  geom_line(color = "grey50") + 
  geom_point(aes(size = n))
A line plot showing how proportion of cancelled flights changes over the course of the day. The proportion starts low at around 0.5% at 6am, then steadily increases over the course of the day until peaking at 4% at 7pm. The proportion of cancelled flights then drops rapidly getting down to around 1% by midnight.
Figure 13.1: x 轴显示计划出发小时,y 轴显示航班取消比例的折线图。 航班取消数量似乎在白天逐渐累积直至晚上 8 点,而深夜航班的取消概率显著降低。

13.4.4 Logarithms

对数是一种极其有用的转换方法,适用于处理跨越多个数量级的数据以及将指数增长转换为线性增长。 在 R 中,您可以选择三种对数函数:log()(自然对数,底数为 e)、log2()(底数为 2)和 log10()(底数为 10)。 我们推荐使用 log2()log10()log2() 易于解释,因为对数尺度上 1 的差异对应原始尺度的倍增,-1 的差异对应减半;而 log10() 则易于反向转换,例如 3 对应 10^3 = 1000。 log() 的反函数是 exp();要计算 log2()log10() 的反函数,您需要使用 2^10^

13.4.5 Rounding

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

round(123.456)
#> [1] 123

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

round(123.456, 2)  # two digits
#> [1] 123.46
round(123.456, 1)  # one digit
#> [1] 123.5
round(123.456, -1) # round to nearest ten
#> [1] 120
round(123.456, -2) # round to nearest hundred
#> [1] 100

关于 round() 有一个奇怪的地方,乍一看可能会让人感到意外:

round(c(1.5, 2.5))
#> [1] 2 2

round() 采用所谓的”四舍五入”或银行家舍入法:若数字恰好处于两个整数中间,则舍入到最近的偶数。 这是一种良好的策略,因为它能保持舍入的无偏性:所有 0.5 中有一半向上舍入,另一半向下舍入。

round() 与始终向下舍入的 floor() 以及始终向上舍入的 ceiling() 配对使用:

x <- 123.456

floor(x)
#> [1] 123
ceiling(x)
#> [1] 124

这些函数没有 digits 参数,因此您可以先缩小数值,进行舍入,然后再放大还原:

# Round down to nearest two digits
floor(x / 0.01) * 0.01
#> [1] 123.45
# Round up to nearest two digits
ceiling(x / 0.01) * 0.01
#> [1] 123.46

如果您想将数值 round() 到其他数字的倍数,也可以采用相同的技术:

# Round to nearest multiple of 4
round(x / 4) * 4
#> [1] 124

# Round to nearest 0.25
round(x / 0.25) * 0.25
#> [1] 123.5

13.4.6 Cutting numbers into ranges

使用 cut()1 将数值向量分割(也称为分 bin)成离散的区间:

x <- c(1, 2, 5, 10, 15, 20)
cut(x, breaks = c(0, 5, 10, 15, 20))
#> [1] (0,5]   (0,5]   (0,5]   (5,10]  (10,15] (15,20]
#> Levels: (0,5] (5,10] (10,15] (15,20]

breaks 不需要等距分布:

cut(x, breaks = c(0, 5, 10, 100))
#> [1] (0,5]    (0,5]    (0,5]    (5,10]   (10,100] (10,100]
#> Levels: (0,5] (5,10] (10,100]

您可以选择自定义 labels。 注意 labels 数量应比 breaks 少一个。

cut(x, 
  breaks = c(0, 5, 10, 15, 20), 
  labels = c("sm", "md", "lg", "xl")
)
#> [1] sm sm sm md lg xl
#> Levels: sm md lg xl

任何超出 breaks 范围的值将变为 NA

y <- c(NA, -10, 5, 10, 30)
cut(y, breaks = c(0, 5, 10, 15, 20))
#> [1] <NA>   <NA>   (0,5]  (5,10] <NA>  
#> Levels: (0,5] (5,10] (10,15] (15,20]

请参阅文档了解其他有用参数,如 rightinclude.lowest,这些参数控制区间是 [a, b) 还是 (a, b],以及是否应包含最低区间 [a, b]

13.4.7 Cumulative and rolling aggregates

Base R 提供了 cumsum()cumprod()cummin()cummax() 函数,用于计算累加和、累积乘积、累计最小值和累计最大值。 dplyr 则提供了 cummean() 函数用于计算累积平均值。 累计求和在实际应用中出现频率最高:

x <- 1:10
cumsum(x)
#>  [1]  1  3  6 10 15 21 28 36 45 55

如果您需要更复杂的滚动或滑动聚合计算,可以尝试使用 slider 包。

13.4.8 Exercises

  1. 用文字解释生成 Figure 13.1 的代码每行分别实现什么功能。

  2. R 提供哪些三角函数? 猜测几个函数名并查阅文档。 它们使用角度还是弧度?

  3. 目前 dep_timesched_dep_time 便于查看但难以计算,因为它们不是真正的连续数值。 运行以下代码可以看到核心问题:每个小时之间存在间隔。

    flights |> 
      filter(month == 1, day == 1) |> 
      ggplot(aes(x = sched_dep_time, y = dep_delay)) +
      geom_point()

    将它们转换为更真实的时间表示形式(自午夜起的分钟数或小数小时数)。

  4. dep_timearr_time 四舍五入到最接近的五分钟。

13.5 General transformations

以下部分介绍了一些常用于数值向量的通用转换方法,但这些方法也可应用于所有其他列类型。

13.5.1 Ranks

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

x <- c(1, 2, 2, 3, 4, NA)
min_rank(x)
#> [1]  1  2  2  4  5 NA

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

min_rank(desc(x))
#> [1]  5  3  3  2  1 NA

如果 min_rank() 无法满足需求,可以查看其变体 dplyr::row_number()dplyr::dense_rank()dplyr::percent_rank()dplyr::cume_dist()。 详情请参阅文档。

df <- tibble(x = x)
df |> 
  mutate(
    row_number = row_number(x),
    dense_rank = dense_rank(x),
    percent_rank = percent_rank(x),
    cume_dist = cume_dist(x)
  )
#> # A tibble: 6 × 5
#>       x row_number dense_rank percent_rank cume_dist
#>   <dbl>      <int>      <int>        <dbl>     <dbl>
#> 1     1          1          1         0          0.2
#> 2     2          2          2         0.25       0.6
#> 3     2          3          2         0.25       0.6
#> 4     3          4          3         0.75       0.8
#> 5     4          5          4         1          1  
#> 6    NA         NA         NA        NA         NA

通过为 base R 的 rank() 函数选择适当的 ties.method 参数,您可以实现许多相同的结果;您可能还需要设置 na.last = "keep" 以将 NA 值保持为 NA

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

df <- tibble(id = 1:10)

df |> 
  mutate(
    row0 = row_number() - 1,
    three_groups = row0 %% 3,
    three_in_each_group = row0 %/% 3
  )
#> # A tibble: 10 × 4
#>      id  row0 three_groups three_in_each_group
#>   <int> <dbl>        <dbl>               <dbl>
#> 1     1     0            0                   0
#> 2     2     1            1                   0
#> 3     3     2            2                   0
#> 4     4     3            0                   1
#> 5     5     4            1                   1
#> 6     6     5            2                   1
#> # ℹ 4 more rows

13.5.2 Offsets

dplyr::lead()dplyr::lag() 允许您引用”当前”值之前或之后的值。 它们会返回与输入长度相同的向量,并在首尾用 NA 填充:

x <- c(2, 5, 11, 11, 19, 35)
lag(x)
#> [1] NA  2  5 11 11 19
lead(x)
#> [1]  5 11 11 19 35 NA
  • x - lag(x) 可计算当前值与先前值的差值。

    x - lag(x)
    #> [1] NA  3  6  0  8 16
  • x == lag(x) 可标识当前值何时发生变化。

    x == lag(x)
    #> [1]    NA FALSE FALSE  TRUE FALSE FALSE

您可以通过第二个参数 n 来指定向前或向后移动多个位置。

13.5.3 Consecutive identifiers

有时您希望在某个事件每次发生时创建新分组。 例如,当查看网站数据时,通常需要将事件划分为多个会话(sessions),若距离上次活动超过 x 分钟则开始新会话。 例如,假设您记录了某人访问网站的时间点:

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

并且您已经计算了每个事件之间的时间间隔,并判断是否存在足够大的间隔:

events <- events |> 
  mutate(
    diff = time - lag(time, default = first(time)),
    has_gap = diff >= 5
  )
events
#> # A tibble: 14 × 3
#>    time  diff has_gap
#>   <dbl> <dbl> <lgl>  
#> 1     0     0 FALSE  
#> 2     1     1 FALSE  
#> 3     2     1 FALSE  
#> 4     3     1 FALSE  
#> 5     5     2 FALSE  
#> 6    10     5 TRUE   
#> # ℹ 8 more rows

但我们如何从这个逻辑向量转换为可以用于 group_by() 的分组变量? 来自 Section 13.4.7cumsum() 可以解决这个问题,因为间隔(即 has_gapTRUE)会使 group 编号增加 1(Section 12.4.2):

events |> mutate(
  group = cumsum(has_gap)
)
#> # A tibble: 14 × 4
#>    time  diff has_gap group
#>   <dbl> <dbl> <lgl>   <int>
#> 1     0     0 FALSE       0
#> 2     1     1 FALSE       0
#> 3     2     1 FALSE       0
#> 4     3     1 FALSE       0
#> 5     5     2 FALSE       0
#> 6    10     5 TRUE        1
#> # ℹ 8 more rows

另一种创建分组变量的方法是使用 consecutive_id(),每当其参数发生变化时就会启动新分组。 例如,受 this stackoverflow question 启发,假设您有一个包含大量重复值的数据框:

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()

df |> 
  group_by(id = consecutive_id(x)) |> 
  slice_head(n = 1)
#> # A tibble: 7 × 3
#> # Groups:   id [7]
#>   x         y    id
#>   <chr> <dbl> <int>
#> 1 a         1     1
#> 2 b         2     2
#> 3 c         4     3
#> 4 d         3     4
#> 5 e         9     5
#> 6 a         4     6
#> # ℹ 1 more row

13.5.4 Exercises

  1. 使用排名函数找出延误最严重的 10 个航班。 您希望如何处理并列值? 请仔细阅读 min_rank() 的文档。

  2. 哪架飞机(tailnum)的准点记录最差?

  3. 如果想尽量避免延误,您应该在一天中的什么时间乘坐飞机?

  4. flights |> group_by(dest) |> filter(row_number() < 4) 会实现什么功能? flights |> group_by(dest) |> filter(row_number(dep_delay) < 4) 会实现什么功能?

  5. 针对每个目的地,计算总延误分钟数。 针对每个航班,计算其航班延误时间占目的地总延误时间的比例。

  6. 延误通常具有时间相关性:即使导致初始延误的问题已经解决,后续航班仍会因等待前序航班离港而延误。 使用 lag() 函数探索某小时的平均航班延误与前一小时平均延误之间的关系。

    flights |> 
      mutate(hour = dep_time %/% 100) |> 
      group_by(year, month, day, hour) |> 
      summarize(
        dep_delay = mean(dep_delay, na.rm = TRUE),
        n = n(),
        .groups = "drop"
      ) |> 
      filter(n > 5)
  7. 检查每个目的地。 您能否找到快得可疑的航班(即可能存在数据录入错误的航班)? 计算航班相对于该目的地最短航班的空中时间。 哪些航班在空中的延误最严重?

  8. 找出所有至少由两家航空公司飞行的目的地。 利用这些目的地,根据各航空公司在同一目的地的表现进行相对排名。

13.6 Numeric summaries

仅使用我们已经介绍的计数、均值和求和等函数就能解决大部分问题,但 R 还提供了许多其他有用的汇总函数。 以下是一些您可能会觉得有用的精选函数。

13.6.1 Center

到目前为止,我们主要使用 mean() 来概括一组值的中心趋势。 正如我们在 Section 3.6 中所见,由于均值是总和除以计数,即使只有几个异常高或低的值也会对其产生影响。 另一种方法是使用 median(),它可以找到位于向量”中间”的值,即 50% 的值高于它,50% 的值低于它。 根据您感兴趣的变量分布形态,均值或中位数可能是更好的中心度量指标。 例如,对于对称分布我们通常报告均值,而对于偏态分布我们通常报告中位数。

Figure 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)) + 
  geom_abline(slope = 1, intercept = 0, color = "white", linewidth = 2) +
  geom_point()
All points fall below a 45° line, meaning that the median delay is always less than the mean delay. Most points are clustered in a  dense region of mean [0, 20] and median [0, 5]. As the mean delay increases, the spread of the median also increases. There are two outlying points with mean ~60, median ~50, and mean ~85, median ~55.
Figure 13.2: 展示使用中位数而非均值汇总每日出发延误差异的散点图。

您可能还会想到众数(mode),即最常见的值。 这种汇总方法仅适用于非常简单的场景(这或许是您在高中学习过它的原因),但对许多真实数据集效果不佳。 如果数据是离散的,可能会出现多个最常见值;如果数据是连续的,则可能不存在最常见值,因为每个值都存在微小差异。 由于这些原因,统计学家往往不使用众数,base R2 中也没有内置众数函数。

13.6.2 Minimum, maximum, and quantiles

如果您关注的是中心点以外的其他位置呢? min()max() 可以给出最大值和最小值。 另一个强大的工具是 quantile(),它是中位数的一种泛化:quantile(x, 0.25) 会找出大于 25% 数值的 x 值,quantile(x, 0.5) 等同于中位数,而 quantile(x, 0.95) 则会找出大于 95% 数值的值。

对于 flights 数据,您可能需要查看延误的 95% 分位数而非最大值,因为这样可以忽略最严重的 5% 极端延误航班。

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"
  )
#> # A tibble: 365 × 5
#>    year month   day   max   q95
#>   <int> <int> <int> <dbl> <dbl>
#> 1  2013     1     1   853  70.1
#> 2  2013     1     2   379  85  
#> 3  2013     1     3   291  68  
#> 4  2013     1     4   288  60  
#> 5  2013     1     5   327  41  
#> 6  2013     1     6   202  51  
#> # ℹ 359 more rows

13.6.3 Spread

有时您不太关心数据主体的分布位置,而更关注其离散程度。 两个常用的汇总指标是标准差 sd(x) 和四分位距 IQR()。 由于您可能已熟悉sd(),这里不再解释,但 IQR() 可能是新概念——它计算 quantile(x, 0.75) - quantile(x, 0.25),给出包含中间 50% 数据的范围。

我们可以借此揭示 flights 数据中的一个小异常。 您可能认为起飞机场和目的地机场之间的距离离散度应该为零,因为机场位置是固定不变的。 但下面的代码揭示了 EGE 机场的数据异常:

flights |> 
  group_by(origin, dest) |> 
  summarize(
    distance_sd = IQR(distance), 
    n = n(),
    .groups = "drop"
  ) |> 
  filter(distance_sd > 0)
#> # A tibble: 2 × 4
#>   origin dest  distance_sd     n
#>   <chr>  <chr>       <dbl> <int>
#> 1 EWR    EGE             1   110
#> 2 JFK    EGE             1   103

13.6.4 Distributions

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

Figure 13.3 显示了出发延误的整体分布情况。 该分布极度偏斜,我们必须放大图表才能看到主要数据。 这表明均值不太可能是一个好的汇总指标,我们可能更倾向于使用中位数。

Two histograms of `dep_delay`. On the left, it's very hard to see any pattern except that there's a very large spike around zero, the bars rapidly decay in height, and for most of the plot, you can't see any bars because they are too short to see. On the right, where we've discarded delays of greater than two hours, we can see that the spike occurs slightly below zero (i.e. most flights leave a couple of minutes early), but there's still a very steep decay after that.
Figure 13.3: (左)完整数据的直方图极度偏斜,难以观察任何细节。(右)将延迟时间限制在两小时内后,可以观察到大多数观测值的分布情况。

检查子组分布是否与整体分布相似也是个好方法。 在下图中,我们叠加了 365 条 dep_delay 的频率多边形(每天一条)。 这些分布似乎遵循共同模式,表明每天使用相同的汇总方法是可行的。

flights |>
  filter(dep_delay < 120) |> 
  ggplot(aes(x = dep_delay, group = interaction(day, month))) + 
  geom_freqpoly(binwidth = 5, alpha = 1/5)

The distribution of `dep_delay` is highly right skewed with a strong peak slightly less than 0. The 365 frequency polygons are mostly  overlapping forming a thick black bland.

不要害怕探索专门为您正在处理的数据量身定制的自定义汇总方法。 在这种情况下,可能意味着分别汇总提前起飞和延误起飞的航班,或者鉴于数值严重偏斜,您可以尝试对数转换。 最后,不要忘记 Section 3.6 中学到的内容:创建数值汇总时,最好包含每个组中的观测数量。

13.6.5 Positions

最后还有一种对数值向量特别有用、但也适用于其他所有值类型的汇总方法:提取特定位置的值,first(x)last(x)nth(x, n)

例如,我们可以查找每天最早和最晚的离港航班:

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)
  )
#> `summarise()` has grouped output by 'year', 'month'. You can override using
#> the `.groups` argument.
#> # A tibble: 365 × 6
#> # Groups:   year, month [12]
#>    year month   day first_dep fifth_dep last_dep
#>   <int> <int> <int>     <int>     <int>    <int>
#> 1  2013     1     1       517       554     2356
#> 2  2013     1     2        42       535     2354
#> 3  2013     1     3        32       520     2349
#> 4  2013     1     4        25       531     2358
#> 5  2013     1     5        14       534     2357
#> 6  2013     1     6        16       555     2355
#> # ℹ 359 more rows

(注:由于 dplyr 函数使用 _ 分隔函数名和参数名,这些函数使用 na_rm 而非 na.rm。)

如果您熟悉 [ 符号(我们将在 Section 27.2 中讨论),可能会疑惑是否还需要这些函数。 原因有三:default 参数允许在指定位置不存在时提供默认值;order_by 参数允许局部覆盖行排序;na_rm 参数允许删除缺失值。

按位置提取值与按排名过滤是互补操作。 过滤操作会提供所有变量,每个观测值单独成行:

flights |> 
  group_by(year, month, day) |> 
  mutate(r = min_rank(sched_dep_time)) |> 
  filter(r %in% c(1, max(r)))
#> # A tibble: 1,195 × 20
#> # Groups:   year, month, day [365]
#>    year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
#>   <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
#> 1  2013     1     1      517            515         2      830            819
#> 2  2013     1     1     2353           2359        -6      425            445
#> 3  2013     1     1     2353           2359        -6      418            442
#> 4  2013     1     1     2356           2359        -3      425            437
#> 5  2013     1     2       42           2359        43      518            442
#> 6  2013     1     2      458            500        -2      703            650
#> # ℹ 1,189 more rows
#> # ℹ 12 more variables: arr_delay <dbl>, carrier <chr>, flight <int>, …

13.6.6 With mutate()

正如其名称所示,汇总函数通常与 summarize() 配对使用。 但由于我们在 Section 13.4.1 讨论过的回收规则,它们也可以有效地与 mutate() 配合使用,特别是在需要进行某种分组标准化时。 例如:

  • x / sum(x) 计算占总量的比例。
  • (x - mean(x)) / sd(x) 计算 Z-score(标准化为均值为 0、标准差为 1)。
  • (x - min(x)) / (max(x) - min(x)) 标准化到 [0, 1] 范围。
  • x / first(x) 基于第一个观测值计算指数。

13.6.7 Exercises

  1. 构思至少五种不同的方法来评估一组航班的典型延误特征。 mean() 在何时有用? median() 在何时有用? 何时可能需要使用其他指标? 应该使用到达延误还是出发延误? 为什么可能需要使用 planes 数据?

  2. 哪些目的地在飞行速度上表现出最大差异?

  3. 创建图表以进一步探索 EGE 机场的异常情况。 能否找到机场搬迁位置的证据? 能否找到其他可能解释这种差异的变量?

13.7 Summary

您已经熟悉许多处理数字的工具,在阅读本章后也了解了如何在 R 中使用它们。 您还学习了一些常用的通用转换方法(虽然不仅限于数值向量),例如排名和偏移量计算。 最后,您实践了多种数值汇总方法,并讨论了若干需要考虑的统计挑战。

接下来两章我们将深入探讨 stringr 包中的字符串处理。 字符串是个庞大的主题,因此将占用两个章节:一章讲解字符串基础,另一章专攻正则表达式。


  1. ggplot2 在 cut_interval()cut_number()cut_width() 中为常见情况提供了一些辅助函数。 必须承认这些函数存在于 ggplot2 中确实有些奇怪,但它们作为直方图计算的一部分非常实用,且编写时间早于 tidyverse 的其他任何组件。↩︎

  2. mode() 函数的功能完全不同!↩︎