分类: R语言

  • R语言做散点图

    plot(
    x = mtcars$wt,
    y = mtcars$mpg,
    main = "mpg vs wt",
    xlab = "wt (重量)",
    ylab = "mpg (每加仑英里数)",
    pch = 19, # 点的形状(19=实心圆)
    col = "steelblue"
    )
    abline(lm(mtcars$mpg ~ mtcars$wt), col = "red", lwd = 2) # 线性回归拟合线
    # install.packages("ggplot2")
    library(ggplot2)
    ggplot(mtcars, aes(x = wt, y = mpg)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE, color = "red")
  • R语言数据转换

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

    引言

    可视化是生成洞见的重要工具,但你很少能拿到正好符合绘图需求的数据形式。通常你需要创建一些新的变量或汇总,以便用数据回答你的问题;或者你只是想重命名变量、重新排列观测值,让数据更易于处理。接下来,你将学习如何完成这些操作(以及更多内容)。通过 dplyr 包介绍数据转换,并使用一个关于 2013 年从纽约市起飞航班的新数据集。

    本章的目标是为你概述用于转换数据框的所有关键工具。我们将从对数据框的行进行操作的函数开始,然后介绍对列进行操作的函数,接着回过头来进一步讨论管道(pipe),这是一个用于组合各种动词的重要工具。随后我们将介绍如何按分组进行操作。最后,本章将以一个案例研究收尾,展示这些函数的实际应用。

    先决条件

    在本章中,我们将重点介绍 dplyr 包。我们将使用 nycflights13 包中的数据来说明关键思想,并使用 ggplot2 帮助我们理解这些数据。

    # install.packages("nycflights13")
    library(nycflights13)
    # install.packages("dplyr")
    library(dplyr)

    请仔细留意在加载 包时打印的冲突信息,如果有的话。它会告诉你 dplyr 覆盖了 base R 中的一些函数。如果在加载 dplyr 之后你仍想使用这些函数的 base 版本,就需要使用它们的完整名称:stats::filter() 和 stats::lag()。到目前为止,我们大多忽略了函数来自哪个包,因为通常这并不重要。然而,了解函数所属的包可以帮助你查找帮助文档以及发现相关函数,因此当我们需要明确函数来源时,会使用与 R 相同的语法:packagename::functionname()。

    nycflights13

    为了探索 dplyr 的基本动词,我们将使用 nycflights13::flights。该数据集包含了 2013 年从纽约市出发的全部 336,776 趟航班。数据来源于美国运输统计局(US Bureau of Transportation Statistics),并在 ?flights 中有文档说明。

    > flights
    # A tibble: 336,776 × 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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 336,766 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    flights 是一个 tibble(tibble,一种数据框的特殊类型)。tidyverse 使用 tibble 来避免一些常见的“坑”。tibbles 和数据框之间最重要的区别在于它们的打印方式:tibble 是为大型数据集设计的,所以它们只显示前几行,以及能够在一屏中放下的列。要查看全部内容,有几种选项。如果你在使用 RStudio,那么最方便的可能是 View(flights),它会打开一个交互式、可滚动且可筛选的视图。否则你也可以使用 print(flights, width = Inf) 来显示所有列,或者使用 glimpse():

    > glimpse(flights)
    Rows: 336,776
    Columns: 19
    $ year <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2
    $ month <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
    $ day <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
    $ dep_time <int> 517, 533, 542, 544, 554, 554, 555, 557, 557, 558, 558, …
    $ sched_dep_time <int> 515, 529, 540, 545, 600, 558, 600, 600, 600, 600, 600, …
    $ dep_delay <dbl> 2, 4, 2, -1, -6, -4, -5, -3, -3, -2, -2, -2, -2, -2, -1
    $ arr_time <int> 830, 850, 923, 1004, 812, 740, 913, 709, 838, 753, 849,…
    $ sched_arr_time <int> 819, 830, 850, 1022, 837, 728, 854, 723, 846, 745, 851,…
    $ arr_delay <dbl> 11, 20, 33, -18, -25, 12, 19, -14, -8, 8, -2, -3, 7, -1
    $ carrier <chr> "UA", "UA", "AA", "B6", "DL", "UA", "B6", "EV", "B6", "…
    $ flight <int> 1545, 1714, 1141, 725, 461, 1696, 507, 5708, 79, 301, 4…
    $ tailnum <chr> "N14228", "N24211", "N619AA", "N804JB", "N668DN", "N394
    $ origin <chr> "EWR", "LGA", "JFK", "JFK", "LGA", "EWR", "EWR", "LGA",…
    $ dest <chr> "IAH", "IAH", "MIA", "BQN", "ATL", "ORD", "FLL", "IAD",…
    $ air_time <dbl> 227, 227, 160, 183, 116, 150, 158, 53, 140, 138, 149, 1
    $ distance <dbl> 1400, 1416, 1089, 1576, 762, 719, 1065, 229, 944, 733, …
    $ hour <dbl> 5, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 6, 6, 6
    $ minute <dbl> 15, 29, 40, 45, 0, 58, 0, 0, 0, 0, 0, 0, 0, 0, 0, 59, 0
    $ time_hour <dttm> 2013-01-01 05:00:00, 2013-01-01 05:00:00, 2013-01-01 0

    在这两种视图中,变量名后面会跟着一些缩写,用来告诉你每个变量的类型:<int> 表示 integer(整数),<dbl> 表示 double(也就是实数),<chr> 表示 character(也就是字符串),而 <dttm> 表示 date-time(日期时间)。这些信息很重要,因为你能对某一列执行的操作在很大程度上取决于它的“类型”。

    dplyr 基础知识

    你即将学习 dplyr 的主要动词(函数)。它们将使你能够解决绝大多数数据处理方面的挑战。但在讨论它们各自的差异之前,先说明它们的共同点:

    • 第一个参数始终是一个数据框(data frame)。
    • 后续参数通常会用变量名(不加引号)来说明要对哪些列进行操作。
    • 输出总是一个新的数据框。
    • 因为每个动词只擅长做一件事,所以要解决复杂问题通常需要把多个动词组合起来,而我们将用管道 |> 来完成。先简要说明:管道会把左边的内容传递给右边的函数,因此 x>f(y) 等价于 f(x,y),而 x>f(y)>g(z) 等价于 g(f(x,y),z)。最容易读出管道的方式是“then”(然后)。即使你还没有学到具体细节,你也依然能对下面的代码有一个大致的理解:

    dplyr 的动词根据它们作用的对象被分为四组:行(rows)、列(columns)、分组(groups)或表(tables)。接下来,你将学习针对行、列和分组最重要的动词。让我们开始吧!

    作用于数据集“行”的最重要动词是 filter() 和 arrange():filter() 会改变哪些行被保留(但不改变它们的顺序),而 arrange() 会改变行的顺序(但不改变哪些行被保留)。这两个函数都只会影响行,列保持不变。我们还会讨论 distinct(),它用于查找具有唯一值的行。与 arrange() 和 filter() 不同,distinct() 除了可以可选地修改列之外,还能实现其功能。

    filter()

    filter() 允许你根据各列的取值来保留行。第一个参数是数据框(data frame)。第二个及后续参数是必须为真(成立)才能保留该行的条件。比如,我们可以找到所有延误超过 120 分钟(两个小时)的航班:

    > flights |>
    + filter(dep_delay > 120)
    # A tibble: 9,723 × 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 848 1835 853 1001 1950
    2 2013 1 1 957 733 144 1056 853
    3 2013 1 1 1114 900 134 1447 1222
    4 2013 1 1 1540 1338 122 2020 1825
    5 2013 1 1 1815 1325 290 2120 1542
    6 2013 1 1 1842 1422 260 1958 1535
    7 2013 1 1 1856 1645 131 2212 2005
    8 2013 1 1 1934 1725 129 2126 1855
    9 2013 1 1 1938 1703 155 2109 1823
    10 2013 1 1 1942 1705 157 2124 1830
    # ℹ 9,713 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    除了 >(大于),你还可以使用 >=(大于等于)、<(小于)、<=(小于等于)、==(等于)和 !=(不等于)。你也可以用 & 或 , 来组合条件,表示“and”(同时满足两个条件);或者用 | 来表示“or”(满足任一条件):

    > flights |>
    + filter(month == 1 & day == 1)
    # A tibble: 842 × 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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 832 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    当你把 | 和 == 结合使用时,有一个很有用的快捷方式:%in%。它会保留那些变量等于右侧任意一个取值的行:

    > flights |>
    + filter(month %in% c(1, 2))
    # A tibble: 51,955 × 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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 51,945 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    当你运行 filter() 时,dplyr 会执行过滤操作,生成一个新的数据框,然后把它打印出来。它不会修改现有的航班数据集,因为 dplyr 的函数从不修改它们的输入。要保存结果,你需要使用赋值运算符 <-

    jan1 <- flights |>
    filter(month == 1 & day == 1)

    常见错误

    当你刚开始学习 R 时,最容易犯的错误是:在测试是否相等时用 =,而不是 ==。一旦发生这种情况,filter() 会提示你:

    > flights |>
    + filter(month = 1)
    Error in `filter()`:
    ! We detected a named input.
    This usually means that you've used `=` instead of `==`.
    ℹ Did you mean `month == 1`?
    Run `rlang::last_trace()` to see where the error occurred.

    另一个错误是你会像在英语里那样去写“or”(或)语句:

    flights |>
    filter(month == 1 | 2)

    这“能用”,就指它不会报错,但它并没有做到你想要的事情:| 会先检查条件 month == 1,然后再检查条件 2,而 2 并不是一个合理的条件来检查。

    arrange()

    arrange() 会根据列的取值来改变行的顺序。它接受一个数据框,以及一组列名(或更复杂的表达式)用于排序。若你提供不止一个列名,那么后续每增加一个列,都将用于在前一个列的取值相同的情况下打破“并列”。例如,下面的代码会按出发时间进行排序,而这个出发时间由四个列共同表示。我们会先得到最早的年份;然后在同一年内部,再得到最早的月份;依此类推。

    > flights |>
    + arrange(year, month, day, dep_time)
    # A tibble: 336,776 × 19
    year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay carrier flight tailnum origin dest air_time distance hour minute
    <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr> <int> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
    1 2013 1 1 517 515 2 830 819 11 UA 1545 N14228 EWR IAH 227 1400 5 15
    2 2013 1 1 533 529 4 850 830 20 UA 1714 N24211 LGA IAH 227 1416 5 29
    3 2013 1 1 542 540 2 923 850 33 AA 1141 N619AA JFK MIA 160 1089 5 40
    4 2013 1 1 544 545 -1 1004 1022 -18 B6 725 N804JB JFK BQN 183 1576 5 45
    5 2013 1 1 554 600 -6 812 837 -25 DL 461 N668DN LGA ATL 116 762 6 0
    6 2013 1 1 554 558 -4 740 728 12 UA 1696 N39463 EWR ORD 150 719 5 58
    7 2013 1 1 555 600 -5 913 854 19 B6 507 N516JB EWR FLL 158 1065 6 0
    8 2013 1 1 557 600 -3 709 723 -14 EV 5708 N829AS LGA IAD 53 229 6 0
    9 2013 1 1 557 600 -3 838 846 -8 B6 79 N593JB JFK MCO 140 944 6 0
    10 2013 1 1 558 600 -2 753 745 8 AA 301 N3ALAA LGA ORD 138 733 6 0
    # ℹ 336,766 more rows
    # ℹ 1 more variable: time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    你可以在 arrange() 里对某个列使用 desc(),根据该列的取值将数据框重新排序为降序(从大到小)。例如,下面这段代码会把航班按延误程度从最多到最少进行排序:

    > flights |>
    + arrange(desc(dep_delay))
    # A tibble: 336,776 × 19
    year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay carrier flight tailnum origin dest air_time distance hour minute
    <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr> <int> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
    1 2013 1 9 641 900 1301 1242 1530 1272 HA 51 N384HA JFK HNL 640 4983 9 0
    2 2013 6 15 1432 1935 1137 1607 2120 1127 MQ 3535 N504MQ JFK CMH 74 483 19 35
    3 2013 1 10 1121 1635 1126 1239 1810 1109 MQ 3695 N517MQ EWR ORD 111 719 16 35
    4 2013 9 20 1139 1845 1014 1457 2210 1007 AA 177 N338AA JFK SFO 354 2586 18 45
    5 2013 7 22 845 1600 1005 1044 1815 989 MQ 3075 N665MQ JFK CVG 96 589 16 0
    6 2013 4 10 1100 1900 960 1342 2211 931 DL 2391 N959DL JFK TPA 139 1005 19 0
    7 2013 3 17 2321 810 911 135 1020 915 DL 2119 N927DA LGA MSP 167 1020 8 10
    8 2013 6 27 959 1900 899 1236 2226 850 DL 2007 N3762Y JFK PDX 313 2454 19 0
    9 2013 7 22 2257 759 898 121 1026 895 DL 2047 N6716C LGA ATL 109 762 7 59
    10 2013 12 5 756 1700 896 1058 2020 878 AA 172 N5DMAA EWR MIA 149 1085 17 0
    # ℹ 336,766 more rows
    # ℹ 1 more variable: time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    请注意,行数并没有变化——我们只是对数据进行重新排序,并没有进行筛选。

    distinct()

    distinct() 用于查找数据集中所有不重复的行,因此从技术上讲,它主要是基于“行”来操作。大多数时候,你会希望得到某些变量的不同组合,因此你也可以(可选地)提供列名:

    > # 移除重复行
    > flights |>
    + distinct()
    # A tibble: 336,776 × 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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 336,766 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows
    > # 移除重复行,保留 origin 和 dest 列
    > flights |>
    + distinct(origin, dest)
    # A tibble: 224 × 2
    origin dest
    <chr> <chr>
    1 EWR IAH
    2 LGA IAH
    3 JFK MIA
    4 JFK BQN
    5 LGA ATL
    6 EWR ORD
    7 EWR FLL
    8 LGA IAD
    9 JFK MCO
    10 LGA ORD
    # ℹ 214 more rows
    # ℹ Use `print(n = ...)` to see more rows

    如果你在筛选不重复的行时还想保留其他列,也可以使用 .keep_all = TRUE 选项。

    > # 移除重复行,保留 origin 和 dest 列,并保留其他列
    > flights |>
    + distinct(origin, dest, .keep_all = TRUE)
    # A tibble: 224 × 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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 214 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    这些航班之所以都出现在 1 月 1 日并非巧合:distinct() 会在数据集中找到某个“唯一行”的第一次出现,并丢弃其余重复项。

    如果你想找的是出现次数,那么最好把 distinct() 换成 count()。通过设置 sort = TRUE,你还可以按“出现次数”从高到低对结果进行排序。

    > # 计算每个 origin 和 dest 组合的航班数量,并按数量降序排序
    > flights |>
    + count(origin, dest, sort = TRUE)
    # A tibble: 224 × 3
    origin dest n
    <chr> <chr> <int>
    1 JFK LAX 11262
    2 LGA ATL 10263
    3 LGA ORD 8857
    4 JFK SFO 8204
    5 LGA CLT 6168
    6 EWR ORD 6100
    7 JFK BOS 5898
    8 LGA MIA 5781
    9 JFK MCO 5464
    10 EWR BOS 5327
    # ℹ 214 more rows
    # ℹ Use `print(n = ...)` to see more rows

    练习

    • 在每个条件的单个流水线(pipeline)中,找到所有满足该条件的航班:
      • 到达延误达到两小时或以上
      • 飞往休斯敦(IAH 或 HOU)
      • 由联合航空(United)、美国航空(American)或达美航空(Delta)运营
      • 在夏季出发(7 月、8 月和 9 月)
      • 到达晚于两小时以上,但出发并没有晚点
      • 至少延误一小时,但在飞行过程中弥补了超过 30 分钟
    library(nycflights13)
    library(dplyr)
    # 1) 到达延误达到两小时或以上
    flights |>
    filter(arr_delay >= 120)
    # 2) 飞往休斯敦(IAH 或 HOU)
    flights |>
    filter(dest %in% c("IAH", "HOU"))
    # 3) 由联合航空(United)、美国航空(American)或达美航空(Delta)运营
    flights |>
    filter(carrier %in% c("UA", "AA", "DL"))
    # 4) 在夏季出发(7 月、8 月和 9 月)
    flights |>
    filter(month %in% c(7, 8, 9))
    # 5) 到达晚于两小时以上,但出发并没有晚点
    flights |>
    filter(arr_delay > 120, dep_delay <= 0)
    # 6) 至少延误一小时,但在飞行过程中弥补了超过 30 分钟
    # (即:dep_delay >= 60 且 arrived_delay 比 dep_delay 少超过 30)
    flights |>
    filter(
    dep_delay >= 60,
    dep_delay - arr_delay > 30
    )
    • 对航班进行排序,以找出出发延误最长的航班。再找出在早晨最早出发的航班。
    flights |>
    arrange(desc(dep_delay), dep_time) |>
    slice(1) # 选择第一行,即 dep_delay 最大的航班
    • 对航班进行排序以找出最快的航班。(提示:试着在函数里加入一个数学计算。)
    flights |>
    mutate(speed = distance / (air_time / 60)) |> # 计算速度,单位为英里/小时
    arrange(desc(speed)) |> # 按照速度降序排序
    slice(1) # 选择第一行,即速度最快的航班
    • 哪些航班飞行距离最远?哪些航班飞行距离最短?
    # 飞行距离最远的一班(含并列)
    flights |>
    filter(distance == max(distance, na.rm = TRUE))
    # 飞行距离最短的一班(含并列)
    flights |>
    filter(distance == min(distance, na.rm = TRUE))
    • 如果你同时使用 filter() 和 arrange(),它们使用的先后顺序会有影响吗?为什么?想一想结果会怎样变化,以及这两个函数需要做多少工作。

    一般建议把 filter() 尽量放前面(越早筛掉越好),让 arrange()/group_by() 等“昂贵”的操作尽可能在更少的行上进行。

    有四个重要的动词会影响列而不改变行:mutate() 创建从现有列派生出来的新列,select() 改变哪些列会出现,rename() 改变列名,relocate() 改变列的顺序位置。

    mutate()

    mutate() 的作用是添加新的列,这些列是根据现有的列计算得来的。你将学习一大组函数,用来处理不同类型的变量并对它们进行操作。就目前而言,我们先专注于基础代数:它让我们可以计算收益,也就是延误航班在空中补回了多少时间,以及速度(英里每小时):

    flights |>
    mutate(
    gain = dep_delay - arr_delay,
    speed = distance / air_time * 60
    )

    默认情况下,mutate() 会把新列添加到数据集的右侧,这使得你很难看清这里发生了什么。我们可以使用 .before 参数,而是把这些变量添加到左侧:

    > flights |>
    + mutate(
    + gain = dep_delay - arr_delay,
    + speed = distance / air_time * 60,
    + .before = 1
    + )
    # A tibble: 336,776 × 21
    gain speed year month day dep_time sched_dep_time dep_delay arr_time
    <dbl> <dbl> <int> <int> <int> <int> <int> <dbl> <int>
    1 -9 370. 2013 1 1 517 515 2 830
    2 -16 374. 2013 1 1 533 529 4 850
    3 -31 408. 2013 1 1 542 540 2 923
    4 17 517. 2013 1 1 544 545 -1 1004
    5 19 394. 2013 1 1 554 600 -6 812
    6 -16 288. 2013 1 1 554 558 -4 740
    7 -24 404. 2013 1 1 555 600 -5 913
    8 11 259. 2013 1 1 557 600 -3 709
    9 5 405. 2013 1 1 557 600 -3 838
    10 -10 319. 2013 1 1 558 600 -2 753
    # ℹ 336,766 more rows
    # ℹ 12 more variables: sched_arr_time <int>, arr_delay <dbl>, carrier <chr>,
    # flight <int>, tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>,
    # distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    . 表示 .before 是传给该函数的一个参数,而不是我们正在创建的第三个新变量的名字。你也可以使用 .after 来把内容加在某个变量之后;在 .before 和 .after 中,你都可以用变量名来代替位置参数。比如,我们可以把新变量加在 day 之后:

    flights |>
    mutate(
    gain = dep_delay - arr_delay,
    speed = distance / air_time * 60,
    .after = day
    )

    另外,你也可以使用 .keep 参数来控制保留哪些变量。尤其有用的参数是 "used":它指定我们只保留在 mutate() 步骤中所涉及或创建的列。例如,下面的输出将只包含变量 dep_delayarr_delayair_timegainhours 和 gain_per_hour

    flights |>
    mutate(
    gain = dep_delay - arr_delay,
    hours = air_time / 60,
    gain_per_hour = gain / hours,
    .keep = "used"
    )

    请注意,由于我们并没有把上面这次计算的结果重新赋值回 flights,因此新的变量 gain、hours 和 gain_per_hour 只会被打印出来,但不会被存储到数据框中。并且,如果我们希望让它们能够在未来使用时继续出现在数据框里,我们就需要认真考虑:是否要把结果赋值回 flights,用更多变量覆盖原始数据框;还是将结果保存到一个新的对象中。很多时候,正确做法是创建一个新的对象,并用具有信息量的名称来表明其内容,例如 delay_gain。当然,你也可能有充分理由选择覆盖掉 flights。

    select()

    拿到包含数百甚至上千个变量的数据集并不罕见。在这种情况下,最先面临的挑战往往只是:关注你真正感兴趣的变量。select() 允许你基于变量名所进行的操作,快速缩小到一个有用的子集:

    按名称选择列:

    flights |>
    select(year, month, day)

    选择从 year 到 day(包含 day)之间的所有列:

    flights |>
    select(year:day)

    选择所有列,但不包括从 year 到 day(含)之间的列:

    flights |>
    select(!year:day)

    从历史上看,这个操作使用的是 - 而不是 !,所以你很可能在实际中会看到这种写法。两个运算符的作用相同,但在行为上有一些微妙差异。我们建议使用 !,因为它读作 “not”,并且与 & 和 | 组合得很自然。

    选择所有字符列:

    flights |>
    select(where(is.character))

    在 select() 内部,你可以使用一些辅助函数:

    • starts_with("abc"):匹配以 “abc” 开头的名称。
    • ends_with("xyz"):匹配以 “xyz” 结尾的名称。
    • contains("ijk"):匹配包含 “ijk” 的名称。
    • num_range("x", 1:3):匹配 x1、x2 和 x3。

    更多细节请查看 ?select。一旦你了解正则表达式,你也就能使用 matches() 来选择与某种模式相匹配的变量。

    你可以在选择 select() 变量时使用 = 来重命名:新的名称出现在 = 的左侧,旧的变量出现在右侧:

    flights |>
    select(tail_num = tailnum)

    rename()

    如果你想保留所有现有变量,只是想重命名其中一些,那么就可以使用 rename() 而不是 select()

    flights |>
    rename(tail_num = tailnum)

    如果你有一堆命名不一致的列,而且手动逐一修正会很痛苦,可以看看 janitor::clean_names(),它提供了一些有用的自动清理功能。

    relocate()

    使用 relocate() 来移动变量。你可能想把相关的变量收集到一起,或者把重要的变量移到前面。默认情况下,relocate() 会把变量移动到最前面:

    flights |>
    relocate(time_hour, air_time)

    You can also specify where to put them using the .before and .after arguments, just like in mutate():

    flights |>
    relocate(year:dep_time, .after = time_hour)
    flights |>
    relocate(starts_with("arr"), .before = dep_time)

    练习

    • 从 flights 中选择 dep_time、dep_delay、arr_time 和 arr_delay。
    flights |>
    select(dep_time, dep_delay, arr_time, arr_delay)
    • 如果在一次 select() 调用中多次指定同一个变量名,会发生什么?

    在 R 的 dplyr::select() 里同一个变量名在同一次 select() 调用中出现多次,通常会导致以下效果之一(取决于版本与具体写法):

    最常见情况:只保留一次
    select() 会对结果进行去重(不会让同一列在输出里重复出现多份)。
    也可能:报错或忽略重复
    如果你的写法属于某些“选择语法”(比如和负选择、重命名、范围选择混在一起),有时会触发提示/报错,或以“后面覆盖/忽略”的方式处理。

    • 运行下面这段代码的结果是否让你感到意外?默认情况下,select 的这些辅助函数是如何处理大写和小写的?你又该如何改变这种默认行为?
    flights |> select(contains("TIME"))

    在 dplyr 的 select() 辅助函数(如 contains())里,默认是不区分大小写(也就是 TIME 会匹配到 time)。

    • 将 air_time 重命名为 air_time_min 以表明其计量单位,并把它移动到数据框的开头。
    flights |>
    rename(air_time_min = air_time) |> # 将 air_time 列重命名为 air_time_min
    relocate(air_time_min, .before = 1)
    • 为什么下面这段代码不起作用?这个报错信息是什么意思?
    > flights |>
    + select(tailnum) |>
    + arrange(arr_delay)
    Error in `arrange()`:
    In argument: `..1 = arr_delay`.
    Caused by error:
    ! object 'arr_delay' not found
    Run `rlang::last_trace()` to see where the error occurred.
    flights |>
    arrange(arr_delay) |> # 先排序再选择
    select(tailnum) # 选择后只有一列

    管道

    我们已经向你展示了上面管道的一些简单示例,但它真正的强大之处在于当你开始把多个动词组合在一起时就体现出来了。比如,假设你想找出飞往休斯顿 IAH 机场的最快航班:你需要把 filter()mutate()select() 和 arrange() 组合在一起:

    flights |>
    filter(dest == "IAH") |>
    mutate(speed = distance / air_time * 60) |>
    select(year:day, dep_time, carrier, flight, speed) |>
    arrange(desc(speed))

    尽管这个管道有四个步骤,但因为每行开头都是动词,所以很容易快速浏览:先从 flights 数据开始,然后 filter,接着 mutate,然后 select,最后 arrange。

    如果我们没有管道(pipe)会怎样呢?我们就需要把每次函数调用嵌套到前一次调用之中:

    arrange(
    select(
    mutate(
    filter(
    flights,
    dest == "IAH"
    ),
    speed = distance / air_time * 60
    ),
    year:day, dep_time, carrier, flight, speed
    ),
    desc(speed)
    )

    或者,我们也可以使用一大堆中间对象:

    flights1 <- filter(flights, dest == "IAH")
    flights2 <- mutate(flights1, speed = distance / air_time * 60)
    flights3 <- select(flights2, year:day, dep_time, carrier, flight, speed)
    arrange(flights3, desc(speed))

    尽管这两种写法各有其适用场景,但管道(pipe)通常会产出更容易编写、也更容易阅读的数据分析代码。

    要在你的代码中加入管道,我们建议使用内置的键盘快捷键 Ctrl/Cmd + Shift + M。你需要在 RStudio 的选项里做一个小调整,以便在代码中使用 |> 代替 %>%,如图 3.1 所示;关于 %>%,我们稍后再讲。

    到目前为止,你已经学习了那些在“行”和“列”上工作的函数。当你加入对“分组(groups)”的处理能力时,dplyr 会变得更加强大。在这一节中,我们将重点介绍最重要的函数:group_by()summarize() 以及 slice 系列的函数。

    group_by()

    使用 group_by() 将你的数据集划分为对分析而言有意义的组:

    > flights |>
    + group_by(month)
    # A tibble: 336,776 × 19
    # Groups: month [12]
    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 533 529 4 850 830
    3 2013 1 1 542 540 2 923 850
    4 2013 1 1 544 545 -1 1004 1022
    5 2013 1 1 554 600 -6 812 837
    6 2013 1 1 554 558 -4 740 728
    7 2013 1 1 555 600 -5 913 854
    8 2013 1 1 557 600 -3 709 723
    9 2013 1 1 557 600 -3 838 846
    10 2013 1 1 558 600 -2 753 745
    # ℹ 336,766 more rows
    # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
    # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    group_by() 不会改变数据,但如果你仔细看输出,会发现输出表明它是“按月分组”的(Groups: month [12])。这意味着后续的操作将改为“按月”进行。group_by() 会为数据框添加这种分组特性(称为 class),从而改变对数据应用的后续动词的行为。

    summarize()

    最重要的分组操作是汇总。如果只用来计算一个汇总统计量,它会把数据框减少为:每个组只有一行。在 dplyr 中,这个操作由 summarize()执行,正如下面的示例所示,它按月计算平均离港延误时间:

    > flights |>
    + group_by(month) |>
    + summarize(
    + avg_delay = mean(dep_delay)
    + )
    # A tibble: 12 × 2
    month avg_delay
    <int> <dbl>
    1 1 NA
    2 2 NA
    3 3 NA
    4 4 NA
    5 5 NA
    6 6 NA
    7 7 NA
    8 8 NA
    9 9 NA
    10 10 NA
    11 11 NA
    12 12 NA

    哎呀!出错了,我们所有的结果都是 NA(读作 “N-A”),这是 R 中表示缺失值的符号。这是因为有些观测到的航班在延误(delay)列中缺少数据,所以当我们在计算均值时把这些值也包含进去,就得到了 NA 的结果。我们会通过将参数 na.rm 设置为 TRUE,来告诉 mean() 函数忽略所有缺失值:

    > flights |>
    + group_by(month) |>
    + summarize(
    + avg_delay = mean(dep_delay, na.rm = TRUE)
    + )
    # A tibble: 12 × 2
    month avg_delay
    <int> <dbl>
    1 1 10.0
    2 2 10.8
    3 3 13.2
    4 4 13.9
    5 5 13.0
    6 6 20.8
    7 7 21.7
    8 8 12.6
    9 9 6.72
    10 10 6.24
    11 11 5.44
    12 12 16.6

    你可以在一次对 summarize() 的调用中创建任意数量的汇总结果。接下来的章节会学到各种有用的汇总方式,但其中一个非常有用的汇总是 n():它会返回每个组中的行数:

    > flights |>
    + group_by(month) |>
    + summarize(
    + avg_delay = mean(dep_delay, na.rm = TRUE),
    + n = n()
    + )
    # A tibble: 12 × 3
    month avg_delay n
    <int> <dbl> <int>
    1 1 10.0 27004
    2 2 10.8 24951
    3 3 13.2 28834
    4 4 13.9 28330
    5 5 13.0 28796
    6 6 20.8 28243
    7 7 21.7 29425
    8 8 12.6 29327
    9 9 6.72 27574
    10 10 6.24 28889
    11 11 5.44 27268
    12 12 16.6 28135

    在数据科学中,均值和计数往往能让你走得出奇地远!

    切片(slice)函数.

    有五个方便的函数,允许你在每个组内提取特定的行:

    • df |> slice_head(n = 1) 从每个组中取出第一行。
    • df |> slice_tail(n = 1) 从每个组中取出最后一行。
    • df |> slice_min(x, n = 1) 取出在列 x 上数值最小的那一行。
    • df |> slice_max(x, n = 1) 取出在列 x 上数值最大的那一行。
    • df |> slice_sample(n = 1) 随机取出一行。

    你可以改变 n 来选择多于一行;或者不使用 n = ...,而是用 prop = 0.1 来选择(例如)每个组中 10% 的行。比如,下面这段代码会找出:到达每个目的地时延误最严重的航班:

    flights |>
    group_by(dest) |> # 按目的地 dest 对数据分组
    slice_max(arr_delay, n = 1) |> # 在每个目的地组内,取 arr_delay 最大的那一行(延误时间最长)
    relocate(dest) # 将列 dest 移动到数据框的前面(放置位置更靠前,方便查看)
    # A tibble: 108 × 19
    # Groups: dest [105]
    dest year month day dep_time sched_dep_time dep_delay arr_time
    <chr> <int> <int> <int> <int> <int> <dbl> <int>
    1 ABQ 2013 7 22 2145 2007 98 132
    2 ACK 2013 7 23 1139 800 219 1250
    3 ALB 2013 1 25 123 2000 323 229
    4 ANC 2013 8 17 1740 1625 75 2042
    5 ATL 2013 7 22 2257 759 898 121
    6 AUS 2013 7 10 2056 1505 351 2347
    7 AVL 2013 8 13 1156 832 204 1417
    8 BDL 2013 2 21 1728 1316 252 1839
    9 BGR 2013 12 1 1504 1056 248 1628
    10 BHM 2013 4 10 25 1900 325 136
    # ℹ 98 more rows
    # ℹ 11 more variables: sched_arr_time <int>, arr_delay <dbl>, carrier <chr>,
    # flight <int>, tailnum <chr>, origin <chr>, air_time <dbl>, distance <dbl>,
    # hour <dbl>, minute <dbl>, time_hour <dttm>
    # ℹ Use `print(n = ...)` to see more rows

    注意,这里有 105 个目的地,但我们得到的是 108 行。怎么回事?slice_min() 和 slice_max() 会保留并列的值,因此当 n = 1 时,它会把所有具有最高值的行都返回。如果你想做到每个组恰好只有一行,可以把 with_ties = FALSE 设置上。

    这和用 summarize() 计算最大延误是类似的,但你得到的是完整的对应行(如果有并列,则是对应的多行),而不是单个汇总统计量。

    按多个变量进行分组

    你可以使用多个变量来创建分组。例如,我们可以为每个日期创建一个组。

    daily <- flights |>
    group_by(year, month, day)
    daily
    #> # A tibble: 336,776 × 19
    #> # 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 533 529 4 850 830
    #> 3 2013 1 1 542 540 2 923 850
    #> 4 2013 1 1 544 545 -1 1004 1022
    #> 5 2013 1 1 554 600 -6 812 837
    #> 6 2013 1 1 554 558 -4 740 728
    #> # ℹ 336,770 more rows
    #> # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>, …

    当你对一个按多个变量分组的 tibble 进行汇总时,每次汇总都会“剥离”掉最后一层分组。回过头看,这并不是让这个函数工作的理想方式,但在不破坏现有代码的前提下很难进行更改。为了让你清楚地知道发生了什么,dplyr 会显示一条消息,告诉你可以如何修改这种行为:

    daily_flights <- daily |>
    summarize(n = n()) # 每个分组里数一共有多少行
    #> `summarise()` has regrouped the output.
    #> ℹ Summaries were computed grouped by year, month, and day.
    #> ℹ Output is grouped by year and month.
    #> ℹ Use `summarise(.groups = "drop_last")` to silence this message.
    #> ℹ Use `summarise(.by = c(year, month, day))` for per-operation grouping
    #> (`?dplyr::dplyr_by`) instead.

    如果你对这种行为感到满意,你也可以明确地要求它,从而抑制该消息:

    daily_flights <- daily |>
    summarize(
    n = n(),
    .groups = "drop_last" # 默认情况,排除最后一层分组
    )

    或者,你可以通过设置不同的值来改变默认行为,例如使用“drop”来删除所有分组,或使用“keep”来保留相同的分组。

    取消分组

    你可能还想在不使用 summarize() 的情况下,从数据框中移除分组。你可以使用 ungroup() 来实现。

    daily |>
    ungroup()
    #> # A tibble: 336,776 × 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 533 529 4 850 830
    #> 3 2013 1 1 542 540 2 923 850
    #> 4 2013 1 1 544 545 -1 1004 1022
    #> 5 2013 1 1 554 600 -6 812 837
    #> 6 2013 1 1 554 558 -4 740 728
    #> # ℹ 336,770 more rows
    #> # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>, …

    现在我们来看看:当你对一个未分组的数据框进行汇总时,会发生什么。

    daily |>
    ungroup() |>
    summarize(
    avg_delay = mean(dep_delay, na.rm = TRUE),
    flights = n() # 分组里有多少行
    )
    #> # A tibble: 1 × 2
    #> avg_delay flights
    #> <dbl> <int>
    #> 1 12.6 336776

    你会得到一行结果,因为 dplyr 会把未分组的数据框中的所有行都视为属于同一个分组。

    .by

    dplyr 1.1.0 引入了一种新的、实验性的、用于“按操作(per-operation)分组”的语法,也就是 .by 参数。group_by() 和 ungroup() 不会被移除,但现在你也可以使用 .by 参数在单个操作内部进行分组:

    > flights |>
    + summarize(
    + delay = mean(dep_delay, na.rm = TRUE),
    + n = n(),
    + .by = month
    + )
    # A tibble: 12 × 3
    month delay n
    <int> <dbl> <int>
    1 1 10.0 27004
    2 10 6.24 28889
    3 11 5.44 27268
    4 12 16.6 28135
    5 2 10.8 24951
    6 3 13.2 28834
    7 4 13.9 28330
    8 5 13.0 28796
    9 6 20.8 28243
    10 7 21.7 29425
    11 8 12.6 29327
    12 9 6.72 27574

    或者,如果你想按多个变量进行分组:

    > flights |>
    + summarize(
    + delay = mean(dep_delay, na.rm = TRUE),
    + n = n(),
    + .by = c(origin, dest)
    + )
    # A tibble: 224 × 4
    origin dest delay n
    <chr> <chr> <dbl> <int>
    1 EWR IAH 11.8 3973
    2 LGA IAH 9.06 2951
    3 JFK MIA 9.34 3314
    4 JFK BQN 6.67 599
    5 LGA ATL 11.4 10263
    6 EWR ORD 14.6 6100
    7 EWR FLL 13.5 3793
    8 LGA IAD 16.7 1803
    9 JFK MCO 10.6 5464
    10 LGA ORD 10.7 8857
    # ℹ 214 more rows
    # ℹ Use `print(n = ...)` to see more rows

    .by 适用于所有动词(verbs),并且它的优点是:当你需要抑制分组提示消息时,不必使用 .groups 参数;等你做完之后也不必再调用 ungroup()。

    案例研究:汇总值与样本量

    每当你进行任何汇总时,最好始终包含计数(n())。这样,你就能确保自己不是基于极少量的数据得出结论。我们将用 Lahman 包中的一些棒球数据来演示这一点。具体来说,我们会比较球员击出安打(H)的比例与他们尝试把球打进场内的次数(AB):

    # install.packages("Lahman")
    library(Lahman)
    # 从 Lahman::Batting 数据集中,按球员 ID 分组
    batters <- Lahman::Batting |>
    group_by(playerID) |>
    summarize(
    # 计算每个球员的表现:安打总数 / 试图击球总数
    performance = sum(H, na.rm = TRUE) / sum(AB, na.rm = TRUE),
    # 统计每个球员的击球尝试次数(样本量)
    n = sum(AB, na.rm = TRUE)
    )
    # 查看结果
    batters
    #> # A tibble: 24,011 × 3
    #> playerID performance n
    #> <chr> <dbl> <int>
    #> 1 aardsda01 0 4
    #> 2 aaronha01 0.305 12364
    #> 3 aaronto01 0.229 944
    #> 4 aasedo01 0 5
    #> 5 abadan01 0.0952 21
    #> 6 abadfe01 0.111 9
    #> # ℹ 24,005 more rows

    当我们将击球手的水平(用打击率 performance 衡量)与击球机会数量(用打数 n 衡量)作图时,你会看到两种模式:

    打数较少的球员,其 performance 的波动更大。这种图的形状非常典型:只要你把均值(或其他汇总统计量)与组大小作图,就会看到随着样本量增加,波动会减小4。

    skill(performance)与击球机会(n)之间存在正相关,因为球队会把最多的击球机会留给最优秀的击球手。

    batters |>
    filter(n > 100) |> # 只保留打数(AB)超过 100 次的球员,排除样本过小、表现波动过大的噪声数据
    ggplot(aes(x = n, y = performance)) + # 创建散点图:横轴为打数(n),纵轴为打击率(performance)
    geom_point(alpha = 1 / 10) + # 绘制透明度为 0.1 的散点,避免大量重叠点导致图像过密("overplotting")
    geom_smooth(se = FALSE) # 添加平滑趋势线(默认用 loess 局部回归),不显示置信区间(se = FALSE)

    注意 ggplot2 与 dplyr 结合使用的便捷模式:只需记住,对数据集进行处理时用 |>,而向图形中添加图层时则切换为 +

    这也对排序(排名)有重要影响。如果你简单地按 desc(performance) 排序,那么打击率最高的人显然往往是那些打数极少却恰好击出安打的球员——他们未必是最优秀的击球手。

    batters |>
    arrange(desc(performance)) # 按打击率(performance)从高到低排序,最高者排在最前面
    #> # A tibble: 24,011 × 3
    #> playerID performance n
    #> <chr> <dbl> <int>
    #> 1 abramge01 1 1 # 注意:这些“100% 打击率”球员仅打过 1–2 次球(n 极小),结果高度不稳定,属统计噪声
    #> 2 alberan01 1 1
    #> 3 averysk01 1 1
    #> 4 banisje01 1 1
    #> 5 barbesa01 1 2
    #> 6 bartocl01 1 1
    #> # ℹ 24,005 more rows # 共 24,011 名球员(按 playerID 去重汇总后)

    总结

    在本章中,你学习了 dplyr 提供的用于处理数据框(data frames)的工具。按用途大致可以分为三类:操作行的工具(例如 filter() 和 arrange())、操作列的工具(例如 select() 和 mutate()),以及操作分组的工具(例如 group_by() 和 summarize())。

  • R语言工作流程:基础

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

    你现在已经有了一些运行 R 代码的经验。我们并没有给出太多细节,但你显然已经掌握了基础!当你开始使用 R 进行编程时,感到沮丧是很正常的,因为它对标点符号非常严格,哪怕一个字符放错位置都会引发报错。不过,虽然你可能会感到些许挫败,但可以放心,这种经历是典型且暂时的:每个人都会经历,而克服它的唯一办法就是不断尝试。

    在我们继续之前,先确保你已经打下了运行 R 代码的坚实基础,并且了解一些最有用的 RStudio 功能。

    代码基础

    让我们回顾一些此前为了尽快让你开始绘图而省略的基础内容。你可以使用 R 来进行基本的数学计算:

    1 / 200 * 30
    (59 + 73 + 2) / 3
    sin(pi / 2)

    你可以使用赋值运算符 <- 来创建新对象:

    x <- 3 * 4
    x

    请注意,x 的值不会被打印出来,它只是被存储了。如果你想查看该值,可以在控制台中输入 x。

    你可以使用 c() 将多个元素组合成一个向量

    primes <- c(2, 3, 5, 7, 11, 13)

    并且,对向量进行的基本算术运算会作用于该向量的每一个元素:

    primes <- c(2, 3, 5, 7, 11, 13)
    primes * 2
    primes - 1

    所有用于创建对象的 R 语句(即赋值语句)都具有相同的形式:

    对象名称 <-

    阅读这段代码时,可以在脑海中将其理解为“对象名获得值”。

    你会进行大量赋值操作,而输入 <- 会有点麻烦。你可以使用 RStudio 的键盘快捷键来节省时间:Alt + -(减号)。注意,RStudio 会自动在 <- 两侧添加空格,这是良好的代码格式习惯。代码即使在状态良好时也可能难以阅读,所以给你的眼睛一点缓解,记得使用空格。

    注释

    在一行中,# 之后的任何文本都会被 R 忽略。这使你可以编写注释,这些文本会被 R 忽略,但供其他人阅读。我们有时会在示例中加入注释来解释代码中发生了什么。

    注释可以帮助简要说明接下来的代码在做什么。

    # create vector of primes
    primes <- c(2, 3, 5, 7, 11, 13)
    # multiply primes by 2
    primes * 2

    对于像这样的简短代码片段,不一定需要为每一行代码都添加注释。但随着你编写的代码变得更加复杂,注释可以为你(以及你的协作者)节省大量理解代码所做工作的时间。

    使用注释来解释代码的“为什么”,而不是“做了什么”或“如何做”。通过仔细阅读代码,总是可以弄清楚代码在做什么以及如何实现,即使过程可能有些繁琐。如果你在注释中描述了每一个步骤,然后又修改了代码,就需要记得同时更新注释,否则将来回看代码时会造成困惑。

    理解为什么要这样做要困难得多,甚至可能无法从代码本身看出来。例如,geom_smooth() 有一个名为 span 的参数,用于控制曲线的平滑程度,值越大曲线越平滑。假设你决定将 span 的值从默认的 0.75 改为 0.9:未来的读者很容易理解发生了什么,但如果你没有在注释中说明你的思考过程,就没有人会知道你为什么要修改默认值。

    对于数据分析代码,使用注释来解释你的整体思路,并在遇到重要见解时将其记录下来。这些信息无法从代码本身重新获取。

    变量名中有什么含义?

    对象名称必须以字母开头,并且只能包含字母、数字、_ 和 .。你应该让对象名称具有描述性,因此需要为多个单词采用一种命名约定。我们推荐使用 snake_case,即用下划线 _ 分隔小写单词。

    i_use_snake_case
    otherPeopleUseCamelCase

    进行一次赋值:

    this_is_a_really_long_name <- 2.5

    要查看这个对象,可以试试 RStudio 的自动补全功能:输入“this”,按 TAB,继续输入字符直到形成唯一前缀,然后按回车。

    假设你犯了一个错误,this_is_a_really_long_name 的值应该是 3.5,而不是 2.5。你可以使用另一个键盘快捷键来帮助修正。例如,可以按 ↑ 调出你刚刚输入的命令并进行编辑。或者,输入“this”,然后按 Cmd/Ctrl + ↑ 来列出所有以这些字母开头的历史命令。使用方向键进行选择,然后按回车重新输入该命令。把 2.5 改为 3.5 并重新运行。

    再进行一次赋值:

    r_rocks <- 2^3

    让我们试着查看它:

    r_rock
    R_rocks

    这说明了你与 R 之间的一种隐含约定:R 会为你完成繁琐的计算,但作为交换,你必须在指令上做到完全精确。否则,你很可能会收到一条错误提示,说明找不到你要的对象。拼写错误很重要;R 无法读懂你的想法,也不会在你输入 r_rock 时自动理解为 r_rocks。大小写同样重要;类似地,当你输入 R_rocks 时,R 也不会自动理解为 r_rocks。

    调用函数

    R 拥有大量内置函数,调用方式如下:

    函数名(参数1 =1, 参数2 =2, …)

    让我们试试使用 seq(),它可以生成规则的数值序列,同时也顺便学习一些 RStudio 的实用功能。输入 se 并按 TAB,会弹出一个列表显示可能的补全项。通过继续输入(如输入 q)来明确选择 seq(),或者使用 ↑/↓ 方向键进行选择。注意弹出的浮动提示,它会提醒你该函数的参数和用途。如果你需要更多帮助,可以按 F1,在右下角的帮助面板中查看完整说明。

    当你选定所需函数后,再按一次 TAB。RStudio 会自动为你补全匹配的开括号 ( 和 闭括号 )。输入第一个参数的名称 from,并将其设为 1。然后输入第二个参数的名称 to,并将其设为 10。最后按回车。

    seq(from = 1, to = 10)

    在函数调用中,我们常常省略前几个参数的名称,因此可以将其改写如下:

    seq(1, 10)

    输入以下代码,并注意 RStudio 也会为成对的引号提供类似的辅助:

    x <- "hello world"

    引号和括号必须始终成对出现。RStudio 会尽力帮助你,但仍然有可能出错并导致不匹配。如果发生这种情况,R 会显示续行符“+”:

    符号 + 表示 R 正在等待更多输入;它认为你还没有完成。通常这意味着你遗漏了一个 ” 或一个 )。要么补上缺失的成对符号,要么按 ESCAPE 取消该表达式并重新尝试。

    请注意,右上角窗格中的环境(Environment)选项卡会显示你创建的所有对象:

    练习

    • 为什么这段代码无法运行?
    my_variable <- 10
    my_varıable

    仔细观察!(这看起来也许毫无意义,但训练你的大脑去注意哪怕最细微的差别,在编程时会带来回报。)

    • 对以下每条 R 命令进行修改,使其能够正确运行:
    libary(todyverse)
    ggplot(dTA = mpg) +
    geom_point(maping = aes(x = displ y = hwy)) +
    geom_smooth(method = "lm)

    正确的版本

    library(ggplot2)
    ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
    geom_point() +
    geom_smooth(method = "lm")
    • 按下 Option + Shift + K / Alt + Shift + K。会发生什么?如何通过菜单进入相同的位置?
    • 运行以下代码。两张图中哪一张被保存为 mpg-plot.png?为什么?
    my_bar_plot <- ggplot(mpg, aes(x = class)) +
    geom_bar()
    my_scatter_plot <- ggplot(mpg, aes(x = cty, y = hwy)) +
    geom_point()
    ggsave(filename = "mpg-plot.png", plot = my_bar_plot)

    保存的是 my_bar_plot(柱状图)

    原因是你在 ggsave() 里明确指定了 plot = my_bar_plot,所以保存的就是这个对象,而不是最近绘制的图(my_scatter_plot)。

    总结

    现在你已经对 R 代码的工作方式有了更多了解,也掌握了一些在将来回看代码时帮助理解它的技巧。接下来,我们将继续你的数据科学之旅,向你介绍 dplyr——tidyverse 中用于数据转换的包,无论是选择重要变量、筛选感兴趣的行,还是计算汇总统计量。

  • R语言数据可视化

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

    引言

    鳍更长的企鹅比鳍更短的企鹅更重还是更轻?你可能已经有了答案,但试着让你的回答更精确一些。鳍长与体重之间的关系是什么样的?是正相关还是负相关?是线性的还是非线性的?这种关系是否会因企鹅的物种不同而变化?又是否会因企鹅所生活的岛屿不同而变化?让我们创建一些可视化图表来回答这些问题。

    你可以使用palmerpenguins中的penguins数据框(也称为palmerpenguins::penguins)来检验你对这些问题的回答。数据框是一种矩形结构的数据集合,由变量(列)和观测(行)组成。penguins数据集包含了由Kristen Gorman博士和南极帕尔默站LTER项目收集并发布的344条观测数据。

    为了让讨论更清晰,我们先定义一些术语:

    变量是你可以测量的数量、特征或属性。

    数值是变量在被测量时所呈现的状态。变量的数值可能会随着每次测量而变化。

    观测是在相似条件下进行的一组测量(通常是在同一时间、针对同一对象完成所有测量)。一次观测会包含多个数值,每个数值对应一个不同的变量。我们有时也将观测称为数据点。

    表格数据是一组数值的集合,每个数值都与一个变量和一个观测相关联。当每个数值都位于独立的“单元格”中、每个变量位于各自的列中、每个观测位于各自的行中时,这种表格数据就是整洁的(tidy)。

    在此语境中,变量指的是所有企鹅的某一属性,而观测指的是单只企鹅的所有属性。

    在控制台中输入数据框的名称,R会打印其内容的预览。请注意,在该预览的顶部显示为tibble。在tidyverse中,我们使用一种称为tibble的特殊数据框,你很快就会进一步了解它。

    > penguins
    # A tibble: 344 × 8
    species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex year
    <fct> <fct> <dbl> <dbl> <int> <int> <fct> <int>
    1 Adelie Torgersen 39.1 18.7 181 3750 male 2007
    2 Adelie Torgersen 39.5 17.4 186 3800 female 2007
    3 Adelie Torgersen 40.3 18 195 3250 female 2007
    4 Adelie Torgersen NA NA NA NA NA 2007
    5 Adelie Torgersen 36.7 19.3 193 3450 female 2007
    6 Adelie Torgersen 39.3 20.6 190 3650 male 2007
    7 Adelie Torgersen 38.9 17.8 181 3625 female 2007
    8 Adelie Torgersen 39.2 19.6 195 4675 male 2007
    9 Adelie Torgersen 34.1 18.1 193 3475 NA 2007
    10 Adelie Torgersen 42 20.2 190 4250 NA 2007
    # ℹ 334 more rows
    # ℹ Use `print(n = ...)` to see more rows

    该数据框包含8列。若想以另一种方式查看数据,从而看到所有变量以及每个变量的前几条观测,可以使用glimpse()。或者,如果你使用的是RStudio,可以运行View(penguins)来打开一个交互式数据查看器。

    该数据框包含8列。若想以另一种方式查看数据,从而看到所有变量以及每个变量的前几条观测,可以运行View(penguins)来打开一个交互式数据查看器。

    penguins中的变量包括:

    species:企鹅的物种(Adelie、Chinstrap或Gentoo)。

    flipper_length_mm:企鹅鳍的长度(单位:毫米)。

    body_mass_g:企鹅的体重(单位:克)。

    要了解更多关于penguins的数据,可以运行 ?penguins 打开其帮助页面。

    最终目标是重现如下可视化图表,用于展示这些企鹅的鳍长与体重之间的关系,并将企鹅的物种纳入考虑。

    创建一个ggplot
    让我们一步一步地重现这个图。

    使用ggplot2时,你从函数ggplot()开始绘图,先定义一个绘图对象,然后在其上添加各个图层。ggplot()的第一个参数是用于绘图的数据集,因此ggplot(data = penguins)会创建一个空的图形,该图已准备好展示penguins数据,但由于我们尚未指定如何进行可视化,目前它仍然是空的。这个图并不十分有趣,但你可以把它看作是一块空白画布,接下来你将在其上绘制图形的各个图层。

    # install.packages("ggplot2")
    library(ggplot2)
    ggplot(data = penguins)

    接下来,我们需要告诉ggplot()如何将数据中的信息以可视化形式呈现。ggplot()函数的mapping参数用于定义数据集中变量如何映射到图形的视觉属性(美学属性)。mapping参数始终在aes()函数中进行定义,而aes()中的x和y参数用于指定映射到x轴和y轴的变量。现在,我们只将鳍长映射到x美学属性,将体重映射到y美学属性。ggplot2会在data参数中查找这些被映射的变量,在这里即为penguins数据。

    下面的图展示了添加这些映射后的结果。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    )

    我们的空白画布现在有了更多结构——可以清楚地看到鳍长将显示在何处(x轴),体重将显示在何处(y轴)。但企鹅本身尚未出现在图中。这是因为我们还没有在代码中说明如何将数据框中的观测值呈现在图上。

    为此,我们需要定义一个几何对象(geom):也就是图中用来表示数据的几何元素。这些几何对象在ggplot2中通过以geom_开头的函数提供。人们通常会根据图中使用的geom类型来描述图表。例如,柱状图使用柱形几何对象(geom_bar()),折线图使用线形几何对象(geom_line()),箱线图使用箱线几何对象(geom_boxplot()),散点图使用点几何对象(geom_point())等等。

    geom_point()函数会为图添加一个点图层,从而创建一个散点图。ggplot2提供了许多geom函数,每个函数都会向图中添加不同类型的图层。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point()
    #> Warning: Removed 2 rows containing missing values or values outside the scale range
    #> (`geom_point()`).

    现在我们得到了一个看起来像“散点图”的图形。它还没有完全达到我们的“最终目标”图,但通过这个图,我们已经可以开始回答推动本次探索的问题:“鳍长与体重之间的关系是什么样的?”这种关系似乎是正相关的(鳍越长,体重也越大),相当线性(点大致分布在一条直线附近而不是曲线附近),并且强度适中(围绕这条线的离散程度不大)。通常来说,鳍更长的企鹅体重也更大。

    在为该图添加更多图层之前,让我们先停下来回顾一下我们收到的警告信息:

    Removed 2 rows containing missing values (geom_point()).

    我们之所以会看到这条消息,是因为数据集中有两只企鹅缺少体重和/或鳍长的数值,而在缺少这两个值的情况下,ggplot2无法在图中表示它们。与R一样,ggplot2遵循这样一种理念:缺失值不应被悄然忽略。这类警告可能是你在处理真实数据时最常见的警告之一——缺失值是一个非常普遍的问题,在后续的图中,我们将抑制此警告,以避免它在我们生成的每一个图旁边都被打印出来。

    散点图对于展示两个数值变量之间的关系非常有用,但对于任何看起来存在的变量关系,都应保持一定的怀疑,并思考是否存在其他变量能够解释或改变这种表面关系的性质。例如,鳍长与体重之间的关系是否会因物种不同而有所差异?让我们将物种信息纳入图中,看看这是否能为这些变量之间的关系带来更多见解。我们将通过使用不同颜色的点来表示不同的物种。

    要实现这一点,我们需要修改美学映射还是几何对象呢?如果你猜的是“在美学映射中,也就是在aes()里面”,那你已经开始掌握使用ggplot2创建数据可视化的要领了!如果没有,也不用担心。你会绘制更多的ggplot图,并在实践中不断检验和提升你的直觉。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)
    ) +
    geom_point()

    当将一个分类变量映射到某个美学属性时,ggplot2会自动为该变量的每一个唯一水平(这里是三个物种)分配一个独特的美学取值(此处为不同颜色),这一过程称为缩放(scaling)。ggplot2还会添加一个图例,用于说明这些取值分别对应哪些水平。

    现在让我们再添加一个图层:一条平滑曲线,用于展示体重与鳍长之间的关系。在继续之前,请回顾上面的代码,并思考我们如何将其添加到现有的图中。

    由于这是一个用于表示数据的新几何对象,我们将添加一个新的geom作为点几何对象之上的一层:geom_smooth()。并且我们将指定使用线性模型绘制最佳拟合直线,即设置method = “lm”。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)
    ) +
    geom_point() +
    geom_smooth(method = "lm")

    我们已经成功添加了线条,但这个图看起来并不像前面的完成图,后者只有一条表示整个数据集的线,而不是为每个企鹅物种分别绘制多条线。

    当美学映射在ggplot()中(即全局层面)定义时,它们会被传递到图中的每一个后续geom图层中。不过,ggplot2中的每个geom函数也可以接受一个mapping参数,从而允许在局部层面定义美学映射,并将其与从全局继承的映射叠加。由于我们希望点根据物种进行着色,但不希望线条按物种分开绘制,因此应只在geom_point()中指定color = species。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point(mapping = aes(color = species)) +
    geom_smooth(method = "lm")

    好了!我们现在得到的图已经非常接近我们的最终目标了,虽然还不够完美。我们仍然需要为每种企鹅使用不同的形状,并改进标签。

    通常来说,仅使用颜色来表示图中的信息并不是一个好主意,因为由于色盲或其他色觉差异,人们对颜色的感知可能不同。因此,除了颜色之外,我们还可以将物种映射到形状这一美学属性上。

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point(mapping = aes(color = species, shape = species)) +
    geom_smooth(method = "lm")

    请注意,图例也会自动更新,以反映点形状的不同。

    最后,我们可以通过在新图层中使用labs()函数来改进图的标签。labs()中的一些参数可能不言自明:title用于添加标题,subtitle用于添加副标题。其他参数与美学映射相对应,x表示x轴标签,y表示y轴标签,而color和shape则定义图例的标签。此外,我们还可以使用ggthemes包中的scale_color_colorblind()函数来优化配色,使其对色盲更友好。

    library(ggthemes)
    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point(aes(color = species, shape = species)) +
    geom_smooth(method = "lm") +
    labs(
    title = "Body mass and flipper length",
    subtitle = "Dimensions for Adelie, Chinstrap, and Gentoo Penguins",
    x = "Flipper length (mm)", y = "Body mass (g)",
    color = "Species", shape = "Species"
    ) +
    scale_color_colorblind()

    我们终于有了一个与我们的“终极目标”完美匹配的图表!

    练习

    • penguins 有多少行?有多少列?
    > dim(penguins)
    [1] 344 8
    • penguins 数据框中的 bill_depth_mm 变量描述的是什么?请查看 ?penguins 的帮助文档来找出答案。

    bill_depth_mm 表示企鹅喙的深度,单位是毫米。

    • 绘制 bill_depth_mm 与 bill_length_mm 的散点图。也就是说,绘制一个以 bill_depth_mm 为纵轴、bill_length_mm 为横轴的散点图。描述这两个变量之间的关系。
    ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm)) +
    geom_point()
    • 如果你绘制 species 与 bill_depth_mm 的散点图,会发生什么?
    ggplot(penguins, aes(x = species, y = bill_depth_mm)) +
    geom_plot()
    • 为什么下面的代码会报错?你会如何修复它?
    ggplot(data = penguins) +
    geom_point()

    这段代码会报错是因为在 ggplot() 中没有指定美学映射(aes),而 geom_point() 需要明确的 x 和 y 变量才能绘制散点图。

    • geom_point() 中的 na.rm 参数有什么作用?该参数的默认值是什么?创建一个散点图,在其中成功地将该参数设置为 TRUE。

    na.rm 参数用于控制是否在绘图时自动移除包含缺失值(NA)的观测。如果设置为 TRUE,会忽略这些缺失值而不报错或警告;如果为 FALSE,则可能会出现警告信息。该参数的默认值是 FALSE。

    ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm)) +
    geom_point(na.rm = TRUE)
    • 将以下说明添加到你在上一个练习中制作的图中:“Data come from the palmerpenguins package.” 提示:查看 labs() 的文档。
    ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm)) +
    geom_point(na.rm = TRUE) +
    labs(caption = "Data come from the palmerpenguins package.")
    • 重现以下可视化。bill_depth_mm 应该映射到哪个美学属性?它应该在全局层级映射,还是在 geom 层级映射?
    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = bill_depth_mm)) +
    geom_point() +
    geom_smooth(se = TRUE)
    • 在脑海中运行这段代码并预测输出结果会是什么样子。然后在 R 中运行代码,检查你的预测。
    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g, color = island)
    ) +
    geom_point() +
    geom_smooth(se = FALSE)
    • 这两个图看起来会不同吗?为什么/为什么不?
    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point() +
    geom_smooth()
    ggplot() +
    geom_point(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_smooth(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    )

    这两段代码会生成相同的图:一个以 flipper_length_mm 为 x 轴、body_mass_g 为 y 轴的散点图,并叠加一条带置信区间的平滑曲线。

    区别在于映射和数据的指定位置:
    第一段在 ggplot() 中全局定义了 data 和 aes,因此后续所有几何对象都会继承。
    第二段在每个 geom 中分别定义 data 和 aes,但由于设置相同,最终效果一致。

    随着我们从这些入门部分继续前进,我们将过渡到更简洁地表达 ggplot2 代码的方式。到目前为止,我们一直非常明确,这在你学习时是有帮助的:

    ggplot(
    data = penguins,
    mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
    geom_point()

    通常,一个函数的前一两个参数非常重要,你应该把它们牢记于心。ggplot() 的前两个参数是 data 和 mapping,其余部分,我们将不再显式写出这些参数名。这样可以减少输入量,并通过减少多余文本,使你更容易看出不同图之间的差异。

    将前一个图用更简洁的方式重写如下:

    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point()

    可视化分布

    如何对一个变量的分布进行可视化取决于该变量的类型:分类变量或数值变量。

    分类变量

    如果一个变量只能取一小部分可能值中的一个,那么它就是分类变量。要查看分类变量的分布,可以使用条形图。条形的高度表示每个 x 值对应的观测数量。

    ggplot(penguins, aes(x = species)) +
    geom_bar()

    对于具有无序水平的分类变量(例如上面的企鹅物种)的条形图,通常更倾向于根据频数对条形进行重新排序。为此,需要将该变量转换为因子(R 处理分类数据的方式),然后对该因子的水平进行重新排序。

    library(forcats)
    ggplot(penguins, aes(x = fct_infreq(species))) +
    geom_bar()

    数值变量

    如果一个变量可以取广泛范围的数值,并且对这些数值进行加减或求平均是有意义的,那么它就是数值型变量(或定量变量)。数值型变量可以是连续的,也可以是离散的。

    连续变量分布的一种常用可视化方式是直方图。

    ggplot(penguins, aes(x = body_mass_g)) +
    geom_histogram(binwidth = 200)

    直方图将 x 轴划分为等宽的区间(bins),并用柱子的高度表示每个区间中的观测数量。在上面的图中,最高的柱子表示有 39 个观测的 body_mass_g 值位于 3500 到 3700 克之间,这两个数值分别是该柱子的左右边界。

    你可以使用 binwidth 参数来设置直方图中区间的宽度,该宽度以 x 变量的单位来衡量。在使用直方图时,你应该始终尝试多种不同的 binwidth,因为不同的区间宽度可能会揭示不同的模式。在下面的图中,binwidth 为 20 时过于狭窄,导致柱子过多,使得难以判断分布的形状。同样,binwidth 为 2,000 时又过大,导致所有数据只被分成三根柱子,也难以判断分布的形状。而 binwidth 为 200 则提供了一个较为合理的平衡。

    数值变量分布的另一种可视化方式是密度图。密度图是直方图的平滑版本,是一种实用的替代方案,特别适用于来自潜在平滑分布的连续数据。我们不会详细介绍 geom_density() 是如何估计密度的(你可以在函数文档中了解更多),但可以通过一个类比来解释密度曲线是如何绘制的。想象一个由木块组成的直方图,然后再想象把一根煮熟的意大利面条从上面放下。面条覆盖在这些木块上的形状,可以看作是密度曲线的形状。它比直方图展示的细节更少,但可以更容易快速看出分布的形状,特别是在峰值和偏度方面。

    ggplot(penguins, aes(x = body_mass_g)) +
    geom_density()

    练习

    • 绘制一个企鹅物种的条形图,并将 species 映射到 y 美学属性。这个图有什么不同?
    ggplot(penguins, aes(y = body_mass_g)) +
    geom_histogram(binwidth = 200)

    将 species 映射到 y 美学属性后,条形图会变为水平条形图,而不是默认的垂直条形图。本质上显示的信息没有变化,仍然是各个企鹅物种的计数,但方向发生了改变,这通常可以让类别名称更容易阅读。

    • 以下两个图有什么不同?在改变柱子的颜色时,哪种美学属性(color 还是 fill)更有用?
    ggplot(penguins, aes(x = species)) +
    geom_bar(color = "red")
    ggplot(penguins, aes(x = species)) +
    geom_bar(fill = "red")

    color = “red” 只会改变柱子的边框颜色为红色,柱子内部仍是默认填充(通常是灰色)。
    fill = “red” 会将柱子的内部填充为红色,而边框颜色保持默认。

    • geom_histogram() 中的 bins 参数有什么作用?

    bins 参数用于指定直方图中分箱(柱子)的数量。它会将数据范围划分为指定数量的区间,每个区间对应一根柱子。调整 bins 的值会影响图形的细致程度:值越大,柱子越多、越细;值越小,柱子越少、越宽。

    • 在加载 ggplot2包后,使用 diamonds 数据集中的 carat 变量绘制直方图。尝试不同的 binwidth。哪个 binwidth 能揭示最有趣的模式?
    ggplot(diamonds,aes(x=carat)) +
    geom_histogram(binwidth = 0.1)

    可视化关系

    为了可视化变量之间的关系,我们需要将至少两个变量映射到图形的美学属性上。接下来,你将学习用于可视化两个或多个变量关系的常用图形,以及用于创建这些图形的几何对象(geoms)。

    一个数值变量和一个分类变量

    为了可视化一个数值变量与一个分类变量之间的关系,我们可以使用并排的箱线图。箱线图是一种用于表示位置度量(百分位数)的可视化简写形式,用来描述数据分布。同时,它也有助于识别潜在的异常值。如图 1.1 所示,每个箱线图由以下部分组成:

    一个盒子表示数据中间一半范围,这个距离称为四分位距(IQR),从分布的第25百分位延伸到第75百分位。盒子中间有一条线表示中位数,即分布的第50百分位。这三条线可以让你了解分布的离散程度,以及分布是围绕中位数对称还是向一侧偏斜。

    用可视化的点表示那些距离盒子边缘超过1.5倍IQR的观测值。这些离群点较为异常,因此会被单独绘制。

    一条线(或称“须”)从盒子的两端延伸,直到分布中最远的非离群点。

    让我们使用 geom_boxplot() 来查看不同物种的体重分布:

    ggplot(penguins, aes(x = species, y = body_mass_g)) +
    geom_boxplot()

    或者,我们可以使用 geom_density() 绘制密度图。

    ggplot(penguins, aes(x = body_mass_g, color = species)) +
    geom_density(linewidth = 0.75)

    我们还通过使用 linewidth 参数自定义了线条的粗细,以便它们在背景中更加突出。

    此外,我们可以将物种映射到颜色和填充美学属性上,并使用 alpha 美学属性为填充的密度曲线添加透明度。该属性的取值范围在 0(完全透明)到 1(完全不透明)之间。在下面的图中,它被设置为 0.5。

    ggplot(penguins, aes(x = body_mass_g, color = species, fill = species)) +
    geom_density(alpha = 0.5)

    请注意我们在这里使用的术语:

    如果我们希望某个美学属性所表示的视觉特征根据变量的取值而变化,就将变量映射到该美学属性上。
    否则,我们是为该美学属性设置一个固定值。

    两个分类变量

    我们可以使用堆叠条形图来可视化两个分类变量之间的关系。例如,下面的两个堆叠条形图都展示了岛屿与物种之间的关系,具体来说,是可视化每个岛屿内物种的分布情况。

    第一个图显示了每个岛屿上各种企鹅物种的频数。该频数图表明每个岛屿上的阿德利企鹅数量是相同的,但我们无法很好地了解各个岛屿内部的百分比分布情况。

    第二个图是通过在几何对象中设置 position = “fill” 创建的相对频率图,对于比较不同岛屿上的物种分布更为有用,因为它不受各岛屿企鹅数量不均的影响。通过该图我们可以看到,Gentoo 企鹅全部生活在 Biscoe 岛,并且约占该岛企鹅数量的 75%;Chinstrap 企鹅全部生活在 Dream 岛,并且约占该岛企鹅数量的 50%;而 Adelie 企鹅生活在三个岛屿上,并且在 Torgersen 岛上占全部企鹅数量。

    在创建这些条形图时,我们将用于分隔成不同条形的变量映射到 x 美学属性,将用于改变条形内部颜色的变量映射到 fill 美学属性。不幸的是,ggplot2 默认将 y 轴标记为 “count”,但我们可以通过添加一个 labs() 图层并将 y 轴标签指定为 “proportion” 来覆盖这一设置。

    ggplot(penguins, aes(x = island, fill = species)) +
    geom_bar(position = "fill") +
    labs(y = "proportion")

    两个数值变量

    到目前为止,你已经学习了如何使用散点图(通过 geom_point() 创建)和光滑曲线(通过 geom_smooth() 创建)来可视化两个数值变量之间的关系。散点图可能是用于展示两个数值变量关系的最常用图形。

    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point()

    三个或更多变量

    正如我们所看到的,我们可以通过将变量映射到其他美学属性上,在图形中加入更多变量。例如,在下面的散点图中,点的颜色代表物种,点的形状代表岛屿。

    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point(aes(color = species, shape = island))

    然而,在图形中添加过多的美学映射会使其显得杂乱且难以理解。另一种对分类变量特别有用的方法是将图形拆分为分面,即每个子图分别展示数据的一个子集。

    要按单个变量对图形进行分面,可以使用 facet_wrap()。facet_wrap() 的第一个参数是一个公式,你可以通过在变量名前加上 ~ 来创建。传递给 facet_wrap() 的变量应当是分类变量。

    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point(aes(color = species, shape = species)) +
    facet_wrap(~island)

    练习

    • ggplot2 包中附带的 mpg 数据框包含了美国环境保护署针对 38 款汽车模型收集的 234 条观测数据。mpg 中哪些变量是分类变量?哪些变量是数值变量?(提示:输入 ?mpg 查看该数据集的文档。)当你运行 mpg 时,如何查看这些信息?
    str(mpg)
    • 使用 mpg 数据框绘制 hwy 与 displ 的散点图。接下来,将第三个数值变量分别映射到颜色、大小、同时映射到颜色和大小,以及形状。对于分类变量和数值变量,这些美学属性的表现有何不同?
    ggplot(mpg, aes(x = displ, y = hwy)) +
    geom_point(aes(color = cty, size = cty, shape = class))
    • 在 hwy 与 displ 的散点图中,如果将第三个变量映射到 linewidth,会发生什么?

    如果在 hwy 与 displ 的散点图中将第三个变量映射到 linewidth,基本不会产生预期效果。

    原因是 linewidth 主要用于线状几何对象(例如 geom_line()、geom_path()),而不是用于点(geom_point())。对于点图,控制大小的美学属性是 size,而不是 linewidth。

    因此:

    在 geom_point() 中使用 linewidth 通常会被忽略或不起作用
    有时可能会出现警告,提示该美学属性未被使用
    如果想根据变量改变点的大小,应使用 size 而不是 linewidth

    • 如果将同一个变量映射到多个美学属性,会发生什么?

    如果将同一个变量映射到多个美学属性(例如同时映射到颜色和大小),这些属性会同时根据该变量发生变化,从而增强数据模式的表达。

    具体来说:

    对于数值变量:颜色会呈现渐变,大小也会连续变化,两者共同反映数值大小规律
    对于分类变量:不同类别会对应不同的颜色和大小组合,更容易区分各组
    不过需要注意:

    映射过多美学属性可能会让图形显得杂乱,降低可读性
    有些组合(如形状 + 大量类别)可能不易区分

    • 绘制 bill_depth_mm 与 bill_length_mm 的散点图,并根据 species 为点着色。按物种着色揭示了这两个变量之间的关系的什么信息?那么按物种分面又会如何?
    # 按物种着色
    ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
    geom_point()
    # 按物种分面
    ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm)) +
    geom_point() +
    facet_wrap(~ species)

    按物种着色可以在同一坐标系中比较不同物种的分布情况,显示各物种之间的差异和可能的分组趋势,但可能会有重叠。

    按物种分面则将每个物种分开展示,使每一类内部的关系更清晰,更容易观察各自的模式,但不如着色方式直观对比不同物种之间的相对位置。

    • 为什么下面的代码会产生两个独立的图例?如何修复以将这两个图例合并?
    ggplot(
    data = penguins,
    mapping = aes(
    x = bill_length_mm, y = bill_depth_mm,
    color = species, shape = species
    )
    ) +
    geom_point() +
    labs(color = "Species") # 注意这里大写了,会被判别为另一个变量,删掉这句即可。
    • 创建以下两个堆叠柱状图。第一个可以回答什么问题?第二个可以回答什么问题?
    ggplot(penguins, aes(x = island, fill = species)) +
    geom_bar(position = "fill")
    ggplot(penguins, aes(x = species, fill = island)) +
    geom_bar(position = "fill")

    第一个图(x = island, fill = species)可以回答:在每个岛上,不同物种的相对比例是多少?

    第二个图(x = species, fill = island)可以回答:对于每个物种,它们分布在不同岛上的相对比例是多少?

    保存你的图形

    一旦你绘制好了图形,可能希望将其从 R 中导出并保存为图像,以便在其他地方使用。这正是 ggsave() 的作用,它会将最近创建的图形保存到磁盘:

    ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
    geom_point()
    ggsave(filename = "penguin-plot.png")

    这会将你的图形保存到你的工作目录中,这个概念你将在后面进一步学习。

    如果你没有指定宽度和高度,它们将取自当前绘图设备的尺寸。为了实现可复现的代码,建议你明确指定这些参数。你可以在文档中了解更多关于 ggsave() 的内容。

    不过,一般来说,我们建议你使用 Quarto 来整合最终报告。Quarto 是一个可复现的写作系统,允许你将代码与文字交织在一起,并自动将图形包含在报告中。

    练习

    • 运行以下代码行。两个图中哪个被保存为 mpg-plot.png?为什么?
    ggplot(mpg, aes(x = class)) +
    geom_bar()
    ggplot(mpg, aes(x = cty, y = hwy)) +
    geom_point()
    ggsave("mpg-plot.png") # 只保存当前图片
    • 在上面的代码中,你需要做哪些修改才能将图形保存为 PDF 而不是 PNG?你可以如何查找 ggsave() 支持哪些图像文件类型?

    常见问题

    当你开始运行 R 代码时,很可能会遇到问题。别担心——这对每个人来说都会发生。我们都已经写了多年的 R 代码,但每天仍然会写出第一次运行就出错的代码!

    首先,仔细对比你正在运行的代码和书中的代码。R 对细节非常敏感,一个放错位置的字符就可能带来完全不同的结果。确保每个 ( 都有对应的 ),每个 ” 都成对出现。有时你运行代码却没有任何反应,这时请查看控制台左侧:如果显示的是 +,说明 R 认为你还没有输入完整的表达式,它正在等待你继续输入。在这种情况下,通常可以按下 ESCAPE 键中止当前命令,然后从头重新开始。

    在创建 ggplot2 图形时,一个常见的问题是把 + 放在了错误的位置:它必须位于行尾,而不是行首。换句话说,确保你没有不小心写出像下面这样的代码:

    ggplot(data = mpg)
    + geom_point(mapping = aes(x = displ, y = hwy))

    如果你仍然卡住了,可以尝试查看帮助。你可以在控制台运行 ?function_name 来获取任何 R 函数的帮助,或者在 RStudio 中选中函数名并按下 F1。即使帮助看起来不太有用也没关系——可以直接跳到示例部分,寻找与你想要实现的内容相匹配的代码。

    如果这仍然没有帮助,请仔细阅读错误信息。有时候答案就隐藏在其中!不过当你刚接触 R 时,即使答案就在错误信息里,你也可能还不知道如何理解它。另一个很有用的工具是 Google:尝试搜索错误信息,很可能已经有人遇到过相同的问题,并在网上获得了解答。

    总结

    在本章中,你已经学习了使用 ggplot2 进行数据可视化的基础知识。我们从 ggplot2 的基本理念开始:可视化是将数据中的变量映射到位置、颜色、大小和形状等美学属性上。随后,你学习了如何通过逐层添加来增加图形的复杂性并提升其表现效果。你还了解了如何通过使用额外的美学映射以及使用分面(faceting)将图形拆分为多个小图,从而可视化单个变量的分布,以及两个或多个变量之间的关系等常见图形。

    在后续内容中,我们会反复使用可视化,并在需要时介绍新的技术。

    在掌握了可视化的基础之后,下面我们将稍微转换一下重点,为你提供一些实用的工作流程建议。

  • R语言的MNIST预测

    # Install and load the required package  
    # install.packages("keras")  
    library(keras)  
    
    # Function to get data  
    get_data <- function() {  
      mnist <- dataset_mnist()  
      x_test <- mnist$test$x  
      y_test <- mnist$test$y  
      
      # Reshape and normalize  
      x_test <- array_reshape(x_test, c(nrow(x_test), 28 * 28)) / 255  
      return(list(x_test = x_test, y_test = y_test))  
    }  
    
    # Sigmoid activation function  
    sigmoid <- function(x) {  
      return(1 / (1 + exp(-x)))  
    }  
    
    # Softmax activation function  
    softmax <- function(x) {  
      exp_x <- exp(x - max(x))  # Prevent overflow  
      return(exp_x / rowSums(exp_x))  
    }  
    
    # Predict function  
    predict <- function(network, x) {  
      W1 <- network$W1  
      W2 <- network$W2  
      W3 <- network$W3  
      b1 <- network$b1  
      b2 <- network$b2  
      b3 <- network$b3  
      
      a1 <- x %*% W1 + b1  
      z1 <- sigmoid(a1)  
      a2 <- z1 %*% W2 + b2  
      z2 <- sigmoid(a2)  
      a3 <- z2 %*% W3 + b3  
      y <- softmax(a3)  
      
      return(y)  
    }  
    
    # Initialize network  
    initialize_network <- function(input_size, hidden_size1, hidden_size2, output_size) {  
      network <- list()  
      network$W1 <- matrix(rnorm(input_size * hidden_size1), nrow = input_size) * 0.01  
      network$b1 <- matrix(0, nrow = 1, ncol = hidden_size1)  
      network$W2 <- matrix(rnorm(hidden_size1 * hidden_size2), nrow = hidden_size1) * 0.01  
      network$b2 <- matrix(0, nrow = 1, ncol = hidden_size2)  
      network$W3 <- matrix(rnorm(hidden_size2 * output_size), nrow = hidden_size2) * 0.01  
      network$b3 <- matrix(0, nrow = 1, ncol = output_size)  
      return(network)  
    }  
    
    # Main program  
    data <- get_data()  # Get test data and labels  
    x <- data$x_test  
    y_true <- data$y_test  
    network <- initialize_network(28 * 28, 128, 64, 10)  # Initialize network  
    
    accuracy_cnt <- 0  
    for (i in 1:nrow(x)) {  
      y <- predict(network, x[i, , drop = FALSE])  # Keep input shape  
      p <- which.max(y) - 1  # R indices start at 1  
      if (p == y_true[i]) {  
        accuracy_cnt <- accuracy_cnt + 1  
      }  
    }  
    
    accuracy <- accuracy_cnt / nrow(x)  
    cat(sprintf("Accuracy: %.2f%%\n", accuracy * 100))  
    
    展开/折叠结果
    Accuracy: 10.10%
    
  • R语言显示MNIST图像和标签

    # Install and load the required package  
    # install.packages("keras")  
    library(keras)  
    
    # Download the MNIST dataset  
    mnist <- dataset_mnist()  
    
    # Separate images and labels  
    train_images <- mnist$train$x  
    train_labels <- mnist$train$y  
    
    # Select an index of the image to view  
    i <- 1  # you can choose any index from 1:60000  
    
    # Display the image  
    image(  
      t(apply(train_images[i,,], 2, rev)),  
      col = gray(seq(0, 1, length = 256)),  
      main = paste("Label:", train_labels[i])  
    )
    
  • R语言加载MNIST

    MNIST是手写数字训练集。需要注意的是,主机要安装Python,如果最新版本的Python无法正确导入数据集,请给Python降级。

    # install.packages("tensorflow")  
    library(tensorflow)  
    install_tensorflow()  
    # install.packages("keras")
    library(keras) 
    mnist <- dataset_mnist()
    

  • R语言实现softmax

    softmax <- function(a){
      max_a <- max(a)
      a <- a-max_a
      exp_a <- exp(a)
      sum_exp_a <- sum(exp_a)
      y <- exp_a/sum_exp_a
      return(y)
    }
    a <- c(1010,1000,990)
    y <- softmax(a)
    y
    
    展开/折叠结果
    [1] 9.999546e-01 4.539787e-05 2.061060e-09
    
  • R语言简单的三层神经网络

    sigmoid <- function(x){
      return(1/(1+exp(-x)))
    }
    X <- c(1,0.5)
    W1 <- matrix(c(0.1,0.2,0.3,0.4,0.5,0.6),nrow = 2)
    B1 <- c(0.1,0.2,0.3)
    A1 <- X %*% W1 + B1
    A1
    Z1 <- sigmoid(A1)
    Z1
    W2 <-  matrix(c(0.1,0.2,0.3,0.4,0.5,0.6),nrow=3,byrow = F)
    B2 <- c(0.1,0.2)
    A2 <-  Z1 %*% W2 +B2
    Z2 <- sigmoid(A2)
    identity <- function(x){
      return(x)
    }
    W3 <- matrix(c(0.1,0.2,0.3,0.4),nrow=2,byrow=F)
    B3 <- c(0.1,0.2)
    A3 <- Z2 %*% W3 +B3
    Y <- identity(A3)
    
    展开/折叠结果
              [,1]      [,2]
    [1,] 0.3168271 0.6962791
    

    用函数表示:

    sigmoid <- function(x){
      return(1/(1+exp(-x)))
    
    identity <- function(x){
      return(x)
    }
    
    init_network_forward <- function(x){
      W1 <- matrix(c(0.1,0.2,0.3,0.4,0.5,0.6),nrow = 2)
      b1 <- c(0.1,0.2,0.3)
      W2 <- matrix(c(0.1,0.2,0.3,0.4,0.5,0.6),nrow=3,byrow = F)
      b2 <- c(0.1,0.2)
      W3 <- matrix(c(0.1,0.2,0.3,0.4),nrow=2,byrow=F)
      b3 <- c(0.1,0.2)
      a1 <- x %*% W1 + b1
      z1 <- sigmoid(a1)
      a2 <- z1 %*% W2 + b2
      z2 <- sigmoid(a2)
      a3 <- z2 %*% W3 +b3
      y <- identity(a3)
      return(y)
    }
    x <- c(1,0.5)
    y <- init_network_forward(x)
    
  • R语言实现最简单神经网络

    该神经网路应用矩阵运算,只有权重,没有激活函数。

    
    

    X <- c(1:2)
    W <- matrix(1:6,nrow = 2,byrow = F)
    Y <- X %*% W

    展开/折叠结果
         [,1] [,2] [,3]
    [1,]    5   11   17