其他
使用 R 语言进行数据处理（二）
为了让大家更好的理解本文的内容，欢迎各位培训班会员参加明天晚上 9 点的直播课：「使用 R 语言进行数据处理（二）」
该课程是「R 语言数据科学的第 5 课时」，课程主页(点击文末的阅读原文即可跳转)：
https://rstata.duanshu.com/#/brief/course/229b770183e44fbbb64133df929818ec
继续上次课的内容。上次课我们学习了filter、arrange、select 和 mutate 四个函数，本次课程中我们将学习 summarise 和 count 函数以及这些函数（包括前面介绍的）的搭配使用。
再进入进入的课程前我们先通过一个案例回顾下上节课学习的内容。有一个小伙伴提了个问题：他有 Q2 到 Q8 共 7 个变量。每个变量取值 1, 2, 3, 4。现在想把所有变量的 1 转成 100, 2 转成 80, 3 转成 60, 4 转成 0。显然这个问题可以使用 mutate_all 解决。
不过我们还需要引入两个函数，if_else 和 case_when。我举两个例子大家就明白它俩的功能了：
library(tidyverse)
tibble(x = 1:10) %>%
mutate(y = if_else(x >= 5, 1, 0))
tibble(x = 1:10) %>%
mutate(y = case_when(
x <= 3 ~ 1,
between(x, 4, 6) ~ 2,
T ~ 3
))
#> # A tibble: 10 × 2
#> x y
#> <int> <dbl>
#> 1 1 0
#> 2 2 0
#> 3 3 0
#> 4 4 0
#> 5 5 1
#> 6 6 1
#> 7 7 1
#> 8 8 1
#> 9 9 1
#> 10 10 1
为了应用 mutate_all 函数，我们先定义一个替换函数：
replace_fun <- function(x){
case_when(
x == 1 ~ 100,
x == 2 ~ 80,
x == 3 ~ 60,
T ~ 0
)
}
#> # A tibble: 10 × 2
#> x y
#> <int> <dbl>
#> 1 1 1
#> 2 2 1
#> 3 3 1
#> 4 4 2
#> 5 5 2
#> 6 6 2
#> 7 7 3
#> 8 8 3
#> 9 9 3
#> 10 10 3
然后我们造个类似的数据框然后应用 mutate_all 和上面定义的函数：
tibble(
Q2 = sample(1:4, 4),
Q3 = sample(1:4, 4),
Q4 = sample(1:4, 4),
Q5 = sample(1:4, 4),
Q6 = sample(1:4, 4),
Q7 = sample(1:4, 4),
Q8 = sample(1:4, 4)
) %>%
mutate_all(replace_fun)
或者使用 mutate 和 across 函数的组合：
tibble(
Q2 = sample(1:4, 4),
Q3 = sample(1:4, 4),
Q4 = sample(1:4, 4),
Q5 = sample(1:4, 4),
Q6 = sample(1:4, 4),
Q7 = sample(1:4, 4),
Q8 = sample(1:4, 4)
) %>%
mutate(across(everything(), ~replace_fun(.x)))
#> # A tibble: 4 × 7
#> Q2 Q3 Q4 Q5 Q6 Q7 Q8
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 60 80 80 80 60 0 0
#> 2 80 0 0 100 80 60 100
#> 3 100 60 60 0 100 80 60
#> 4 0 100 100 60 0 100 80
是不是很方便。下面我们开始今天的学习吧！
分组聚合
和上节课一样，我们继续使用 flights 数据集演示：
nycflights13::flights -> flights
flights
#> # A tibble: 336,776 × 19
#> year month day dep_time sched_de…¹ dep_d…² arr_t…³ sched…⁴ arr_d…⁵ carrier
#> <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr>
#> 1 2013 1 1 517 515 2 830 819 11 UA
#> 2 2013 1 1 533 529 4 850 830 20 UA
#> 3 2013 1 1 542 540 2 923 850 33 AA
#> 4 2013 1 1 544 545 -1 1004 1022 -18 B6
#> 5 2013 1 1 554 600 -6 812 837 -25 DL
#> 6 2013 1 1 554 558 -4 740 728 12 UA
#> 7 2013 1 1 555 600 -5 913 854 19 B6
#> 8 2013 1 1 557 600 -3 709 723 -14 EV
#> 9 2013 1 1 557 600 -3 838 846 -8 B6
#> 10 2013 1 1 558 600 -2 753 745 8 AA
#> # … with 336,766 more rows, 9 more variables: flight <int>, tailnum <chr>,
#> # origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, hour <dbl>,
#> # minute <dbl>, time_hour <dttm>, and abbreviated variable names
#> # ¹sched_dep_time, ²dep_delay, ³arr_time, ⁴sched_arr_time, ⁵arr_delay
我们可以把所有的观测值看做一组，然后计算 dep_delay 变量的均值：
flights %>%
summarise(delay = mean(dep_delay, na.rm = TRUE))
#> # A tibble: 1 × 1
#> delay
#> <dbl>
#> 1 12.6
我们也可以按照年月日进行分组，也就是每天一组，然后计算每天所有航班的 dep_delay 的均值：
flights %>%
group_by(year, month, day) %>%
summarise(delay = mean(dep_delay, na.rm = TRUE))
#> # A tibble: 365 × 4
#> # Groups: year, month [12]
#> year month day delay
#> <int> <int> <int> <dbl>
#> 1 2013 1 1 11.5
#> 2 2013 1 2 13.9
#> 3 2013 1 3 11.0
#> 4 2013 1 4 8.95
#> 5 2013 1 5 5.73
#> 6 2013 1 6 7.15
#> 7 2013 1 7 5.42
#> 8 2013 1 8 2.55
#> 9 2013 1 9 2.28
#> 10 2013 1 10 2.84
#> # … with 355 more rows
再例如按照 dest 变量进行分组：
flights %>%
group_by(dest) %>%
summarise(count = n(),
dist = mean(distance, na.rm = TRUE),
delay = mean(arr_delay, na.rm = TRUE)) %>%
dplyr::filter(count > 20 & dest != "HNL") %>%
ggplot(aes(dist, delay)) +
geom_point(aes(size = count), alpha = 1/3) +
geom_smooth()
计算均值的时候一定要注意缺失值的问题，假如我们不使用 na.rm = TRUE 参数：
flights %>%
group_by(year, month, day) %>%
summarise(mean = mean(dep_delay))
#> # A tibble: 365 × 4
#> # Groups: year, month [12]
#> year month day mean
#> <int> <int> <int> <dbl>
#> 1 2013 1 1 NA
#> 2 2013 1 2 NA
#> 3 2013 1 3 NA
#> 4 2013 1 4 NA
#> 5 2013 1 5 NA
#> 6 2013 1 6 NA
#> 7 2013 1 7 NA
#> 8 2013 1 8 NA
#> 9 2013 1 9 NA
#> 10 2013 1 10 NA
#> # … with 355 more rows
结果我们得到了一堆 NAs。因为缺失值在计算机里面是使用一个很大的数表示的，和任何数的和都是缺失值，均值自然也是缺失值。所以我们需要 na.rm = TRUE 参数来预先去除 NAs。
思考：能否先 filter 非缺失的值然后再求均值？代码怎么写？
观测值计数
首先我们从数据集中筛选出没有被取消的航班：
flights %>%
dplyr::filter(!is.na(dep_delay) & !is.na(arr_delay))
#> # A tibble: 327,346 × 19
#> year month day dep_time sched_de…¹ dep_d…² arr_t…³ sched…⁴ arr_d…⁵ carrier
#> <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr>
#> 1 2013 1 1 517 515 2 830 819 11 UA
#> 2 2013 1 1 533 529 4 850 830 20 UA
#> 3 2013 1 1 542 540 2 923 850 33 AA
#> 4 2013 1 1 544 545 -1 1004 1022 -18 B6
#> 5 2013 1 1 554 600 -6 812 837 -25 DL
#> 6 2013 1 1 554 558 -4 740 728 12 UA
#> 7 2013 1 1 555 600 -5 913 854 19 B6
#> 8 2013 1 1 557 600 -3 709 723 -14 EV
#> 9 2013 1 1 557 600 -3 838 846 -8 B6
#> 10 2013 1 1 558 600 -2 753 745 8 AA
#> # … with 327,336 more rows, 9 more variables: flight <int>, tailnum <chr>,
#> # origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, hour <dbl>,
#> # minute <dbl>, time_hour <dttm>, and abbreviated variable names
#> # ¹sched_dep_time, ²dep_delay, ³arr_time, ⁴sched_arr_time, ⁵arr_delay
然后使用 tailnum 变量分组计算每架航班的平均延误时间：
flights %>%
dplyr::filter(!is.na(dep_delay) & !is.na(arr_delay)) %>%
group_by(tailnum) %>%
summarise(delay = mean(arr_delay))
#> # A tibble: 4,037 × 2
#> tailnum delay
#> <chr> <dbl>
#> 1 D942DN 31.5
#> 2 N0EGMQ 9.98
#> 3 N10156 12.7
#> 4 N102UW 2.94
#> 5 N103US -6.93
#> 6 N104UW 1.80
#> 7 N10575 20.7
#> 8 N105UW -0.267
#> 9 N107US -5.73
#> 10 N108UW -1.25
#> # … with 4,027 more rows
但是我们知道，均值的最大缺点就是易受极端值的影响，因此我们先观察下 delay 的分布:
flights %>%
dplyr::filter(!is.na(dep_delay) & !is.na(arr_delay)) %>%
group_by(tailnum) %>%
summarise(delay = mean(arr_delay)) %>%
ggplot(aes(x = delay)) +
geom_freqpoly(binwidth = 10)
可以看到甚至有航班平均延误时间超过 300 分钟，因此这些航班会极大的影响我们的均值分析。比较下面两幅图：
flights %>%
dplyr::filter(!is.na(dep_delay) & !is.na(arr_delay)) %>%
group_by(tailnum) %>%
summarise(
delay = mean(arr_delay, na.rm = TRUE),
n = n()
) -> delays
ggplot(data = delays, mapping = aes(x = n, y = delay)) +
geom_point(alpha = 1/10)
delays %>%
dplyr::filter(n > 25) %>%
ggplot(mapping = aes(x = n, y = delay)) +
geom_point(alpha = 1/10)
有用的汇总函数
mean() & median()：计算均值和中位数
flights %>%
dplyr::filter(!is.na(dep_delay) & !is.na(arr_delay)) -> not_cancelled
not_cancelled %>%
group_by(year, month, day) %>%
summarise(
avg_delay1 = mean(arr_delay),
avg_delay2 = mean(arr_delay[arr_delay > 0]), # 计算所有正值的均值
median_delay = median(arr_delay)
) %>%
ungroup()
#> # A tibble: 365 × 6
#> year month day avg_delay1 avg_delay2 median_delay
#> <int> <int> <int> <dbl> <dbl> <dbl>
#> 1 2013 1 1 12.7 32.5 3
#> 2 2013 1 2 12.7 32.0 4
#> 3 2013 1 3 5.73 27.7 1
#> 4 2013 1 4 -1.93 28.3 -8
#> 5 2013 1 5 -1.53 22.6 -7
#> 6 2013 1 6 4.24 24.4 -1
#> 7 2013 1 7 -4.95 27.8 -10
#> 8 2013 1 8 -3.23 20.8 -7
#> 9 2013 1 9 -0.264 25.6 -6
#> 10 2013 1 10 -5.90 27.3 -11
#> # … with 355 more rows
sd(), IQR(), mad()：计算标准差、内距和中位数绝对偏差
例如看看哪个目的地的飞机飞行航程的波动最大：
not_cancelled %>%
group_by(dest) %>%
summarise(distance_sd = sd(distance)) %>%
arrange(desc(distance_sd))
#> # A tibble: 104 × 2
#> dest distance_sd
#> <chr> <dbl>
#> 1 EGE 10.5
#> 2 SAN 10.4
#> 3 SFO 10.2
#> 4 HNL 10.0
#> 5 SEA 9.98
#> 6 LAS 9.91
#> 7 PDX 9.87
#> 8 PHX 9.86
#> 9 LAX 9.66
#> 10 IND 9.46
#> # … with 94 more rows
min(), max(), quantile(x, 0.25)：计算最小值，最大值和分位数
not_cancelled %>%
group_by(year, month, day) %>%
summarise(
first = min(dep_time),
last = max(dep_time)
)
#> # A tibble: 365 × 5
#> # Groups: year, month [12]
#> year month day first last
#> <int> <int> <int> <int> <int>
#> 1 2013 1 1 517 2356
#> 2 2013 1 2 42 2354
#> 3 2013 1 3 32 2349
#> 4 2013 1 4 25 2358
#> 5 2013 1 5 14 2357
#> 6 2013 1 6 16 2355
#> 7 2013 1 7 49 2359
#> 8 2013 1 8 454 2351
#> 9 2013 1 9 2 2252
#> 10 2013 1 10 3 2320
#> # … with 355 more rows
first(), last(), nth(x, 2)：第一个、最后一个和第 n 个
not_cancelled %>%
group_by(year, month, day) %>%
summarise(
first_dep = first(dep_time),
last_dep = last(dep_time)
)
#> # A tibble: 365 × 5
#> # Groups: year, month [12]
#> year month day first_dep last_dep
#> <int> <int> <int> <int> <int>
#> 1 2013 1 1 517 2356
#> 2 2013 1 2 42 2354
#> 3 2013 1 3 32 2349
#> 4 2013 1 4 25 2358
#> 5 2013 1 5 14 2357
#> 6 2013 1 6 16 2355
#> 7 2013 1 7 49 2359
#> 8 2013 1 8 454 2351
#> 9 2013 1 9 2 2252
#> 10 2013 1 10 3 2320
#> # … with 355 more rows
n(), n_distinct() 和 count()：计数、计算互不相同的值的个数和变量计数
not_cancelled %>%
group_by(dest) %>%
summarise(carriers = n_distinct(carrier)) %>%
arrange(desc(carriers))
not_cancelled %>%
count(dest)
#> # A tibble: 104 × 2
#> dest carriers
#> <chr> <int>
#> 1 ATL 7
#> 2 BOS 7
#> 3 CLT 7
#> 4 ORD 7
#> 5 TPA 7
#> 6 AUS 6
#> 7 DCA 6
#> 8 DTW 6
#> 9 IAD 6
#> 10 MSP 6
#> # … with 94 more rows
count 函数还可以进行加权计数：
not_cancelled %>%
count(tailnum, wt = distance)
#> # A tibble: 104 × 2
#> dest n
#> <chr> <int>
#> 1 ABQ 254
#> 2 ACK 264
#> 3 ALB 418
#> 4 ANC 8
#> 5 ATL 16837
#> 6 AUS 2411
#> 7 AVL 261
#> 8 BDL 412
#> 9 BGR 358
#> 10 BHM 269
#> # … with 94 more rows
# 这等价于
not_cancelled %>%
group_by(tailnum) %>%
summarise(n = sum(distance)) %>%
ungroup()
#> # A tibble: 4,037 × 2
#> tailnum n
#> <chr> <dbl>
#> 1 D942DN 3418
#> 2 N0EGMQ 239143
#> 3 N10156 109664
#> 4 N102UW 25722
#> 5 N103US 24619
#> 6 N104UW 24616
#> 7 N10575 139903
#> 8 N105UW 23618
#> 9 N107US 21677
#> 10 N108UW 32070
#> # … with 4,027 more rows
对逻辑值的运算：sum(x > 10), mean(y == 0)
例如计数每天 dep_time < 500 的航班数量：
not_cancelled %>%
group_by(year, month, day) %>%
summarise(n_early = sum(dep_time < 500))
#> # A tibble: 365 × 4
#> # Groups: year, month [12]
#> year month day n_early
#> <int> <int> <int> <int>
#> 1 2013 1 1 0
#> 2 2013 1 2 3
#> 3 2013 1 3 4
#> 4 2013 1 4 3
#> 5 2013 1 5 3
#> 6 2013 1 6 2
#> 7 2013 1 7 2
#> 8 2013 1 8 1
#> 9 2013 1 9 3
#> 10 2013 1 10 3
#> # … with 355 more rows
再例如计算 arr_delay > 60 的概率(比例)：
not_cancelled %>%
group_by(year, month, day) %>%
summarise(hour_prop = mean(arr_delay > 60))
#> # A tibble: 365 × 4
#> # Groups: year, month [12]
#> year month day hour_prop
#> <int> <int> <int> <dbl>
#> 1 2013 1 1 0.0722
#> 2 2013 1 2 0.0851
#> 3 2013 1 3 0.0567
#> 4 2013 1 4 0.0396
#> 5 2013 1 5 0.0349
#> 6 2013 1 6 0.0470
#> 7 2013 1 7 0.0333
#> 8 2013 1 8 0.0213
#> 9 2013 1 9 0.0202
#> 10 2013 1 10 0.0183
#> # … with 355 more rows
mutate 和 filter 函数在分组中的应用
例如筛选出来每天延误最严重的 10 次航班：
flights %>%
group_by(year, month, day) %>%
dplyr::filter(rank(desc(arr_delay)) < 10)
#> # A tibble: 3,306 × 19
#> # Groups: year, month, day [365]
#> year month day dep_time sched_de…¹ dep_d…² arr_t…³ sched…⁴ arr_d…⁵ carrier
#> <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr>
#> 1 2013 1 1 848 1835 853 1001 1950 851 MQ
#> 2 2013 1 1 1815 1325 290 2120 1542 338 EV
#> 3 2013 1 1 1842 1422 260 1958 1535 263 EV
#> 4 2013 1 1 1942 1705 157 2124 1830 174 MQ
#> 5 2013 1 1 2006 1630 216 2230 1848 222 EV
#> 6 2013 1 1 2115 1700 255 2330 1920 250 9E
#> 7 2013 1 1 2205 1720 285 46 2040 246 AA
#> 8 2013 1 1 2312 2000 192 21 2110 191 EV
#> 9 2013 1 1 2343 1724 379 314 1938 456 EV
#> 10 2013 1 2 1244 900 224 1431 1104 207 EV
#> # … with 3,296 more rows, 9 more variables: flight <int>, tailnum <chr>,
#> # origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, hour <dbl>,
#> # minute <dbl>, time_hour <dttm>, and abbreviated variable names
#> # ¹sched_dep_time, ²dep_delay, ³arr_time, ⁴sched_arr_time, ⁵arr_delay
筛选出平均每天至少一趟航班的目的地，也就是一年中的航班数量大于 365 的：
flights %>%
group_by(dest) %>%
dplyr::filter(n() > 365)
#> # A tibble: 332,577 × 19
#> # Groups: dest [77]
#> year month day dep_time sched_de…¹ dep_d…² arr_t…³ sched…⁴ arr_d…⁵ carrier
#> <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr>
#> 1 2013 1 1 517 515 2 830 819 11 UA
#> 2 2013 1 1 533 529 4 850 830 20 UA
#> 3 2013 1 1 542 540 2 923 850 33 AA
#> 4 2013 1 1 544 545 -1 1004 1022 -18 B6
#> 5 2013 1 1 554 600 -6 812 837 -25 DL
#> 6 2013 1 1 554 558 -4 740 728 12 UA
#> 7 2013 1 1 555 600 -5 913 854 19 B6
#> 8 2013 1 1 557 600 -3 709 723 -14 EV
#> 9 2013 1 1 557 600 -3 838 846 -8 B6
#> 10 2013 1 1 558 600 -2 753 745 8 AA
#> # … with 332,567 more rows, 9 more variables: flight <int>, tailnum <chr>,
#> # origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>, hour <dbl>,
#> # minute <dbl>, time_hour <dttm>, and abbreviated variable names
#> # ¹sched_dep_time, ²dep_delay, ³arr_time, ⁴sched_arr_time, ⁵arr_delay
因此在 group_by() 之后记得及时 ungroup，否则后面的运行仍然是分组进行的。
案例时间：使用全球新冠疫情数据进行练习并解决下面的问题：
library(haven)
read_dta('world-covid19.dta') -> df
df
#> # A tibble: 29,673 × 7
#> iso country country_en date confirmed recovered deaths
#> <chr> <chr> <chr> <date> <dbl> <dbl> <dbl>
#> 1 BTN 不丹 Bhutan 2020-01-22 0 0 0
#> 2 BTN 不丹 Bhutan 2020-01-23 0 0 0
#> 3 BTN 不丹 Bhutan 2020-01-24 0 0 0
#> 4 BTN 不丹 Bhutan 2020-01-25 0 0 0
#> 5 BTN 不丹 Bhutan 2020-01-26 0 0 0
#> 6 BTN 不丹 Bhutan 2020-01-27 0 0 0
#> 7 BTN 不丹 Bhutan 2020-01-28 0 0 0
#> 8 BTN 不丹 Bhutan 2020-01-29 0 0 0
#> 9 BTN 不丹 Bhutan 2020-01-30 0 0 0
#> 10 BTN 不丹 Bhutan 2020-01-31 0 0 0
#> # … with 29,663 more rows
筛选出来美国的数据； 筛选出来数据集中最后一天的数据； 计算每天的全球总确诊人数、总治愈人数和总病死人数； 筛选出来截止最后一天确诊人数最多的 10 个国家的数据并绘图展示； 筛选出来每天确诊人数最多的 10 个国家； 将处理后的数据导出为 xlsx 格式，csv 格式 和 rds 格式； 计算每个月各个国家的确诊人数。 大家可以试试自己再提一些新的问题。
我提供了一份参考答案，见附件中的 answers.R 文件。
直播信息
为了让大家更好的理解上面的内容，欢迎各位培训班会员参加明天晚上 9 点的直播课：「使用 R 语言进行数据处理（二）」
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要报名 RStata 培训班，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
预定纸质版讲义材料
在之前的课程结尾我把讲义材料整理成了一本讲义材料，感兴趣的小伙伴可以联系微信号 r_stata 李老师预定（R 语言数据科学是 75 元/本，355 页彩色胶装）：
今日购买长期会员有赠送纸质书哦！