R数据科学，第一部分：探索

R data science handbook
R for data science
发布日期:   2022-01-11
更新日期:   2022-10-10

Part1 探索

suppressMessages(library(tidyverse))
suppressMessages(library(ggplot2))
suppressMessages(library(dplyr))

使用ggplot2进行数据可视化

head(mpg,3)
A tibble: 3 × 11
manufacturermodeldisplyearcyltransdrvctyhwyflclass
<chr><chr><dbl><int><int><chr><chr><int><int><chr><chr>
audia41.819994auto(l5) f1829pcompact
audia41.819994manual(m5)f2129pcompact
audia42.020084manual(m6)f2031pcompact
 
options(repr.plot.width=8, repr.plot.height=6)
ggplot(data = mpg, aes(x=displ, y=hwy, color=class)) + geom_point()
ggplot(data = mpg, aes(x=displ, y=hwy, alpha=class)) + geom_point()
Warning message:
“Using alpha for a discrete variable is not advised.”

# 分面
ggplot(data = mpg, aes(x=displ, y=hwy)) + geom_point() + facet_wrap(~ class, nrow = 2)
ggplot(data = mpg, aes(x=displ, y=hwy)) + geom_smooth(method = 'loess') + geom_point()
`geom_smooth()` using formula 'y ~ x'

geom_smooth()函数可以按照不同的线型绘制出不同的曲线，每条曲线对应映射到线型的变量的一个唯一值

ggplot(data = mpg, aes(x=displ, y=hwy)) + geom_smooth(method = 'loess',aes(color=drv, linetype=drv)) + geom_point(aes(color=drv))
`geom_smooth()` using formula 'y ~ x'

也可以为不同的图层指定不同的数据。

ggplot(data = mpg, aes(x=displ, y=hwy)) + 
 geom_point(aes(color=class)) + 
 geom_smooth(data = filter(mpg, class=='subcompact'), se = FALSE) # se: whether show confidence level or not
`geom_smooth()` using method = 'loess' and formula 'y ~ x'

绘图时用来计算新数据的算法称为stat（statistical transformtaion，统计变换）。

通过查看stat参数的默认值，可以知道几何对象函数使用了哪些统计变换。例如，?geom_bar显示出stat的默认值是count，这说明geom_bar()使用stat_count()函数进行统计变换。

通常来说，几何对象函数和统计变换函数可以互换使用。

gplot(data = diamonds) + geom_bar(mapping = aes(x=cut))
ggplot(data = diamonds, aes(x=cut)) + stat_count()

可以这样做的原因是，每个几何对象函数都有一个默认统计变换，每个统计变换函数都有一个默认几何对象

想要显示使用某个统计变换的3个原因：

  • 覆盖默认的统计变换
demo <- tribble(
 ~a, ~b,
 "bar_1", 20,
 "bar_2", 30,
 "bar_3", 40)
ggplot(data = demo, aes(x=a, y=b)) + geom_bar(stat = "identity")
  • 覆盖从统计变换生成的变量到图形属性的默认映射。（例如，想显示一张表示比例的条性图，而不是计数）
ggplot(data = diamonds, aes(x=cut)) + geom_bar()
ggplot(data = diamonds, aes(x=cut, y=..prop.., group=1)) + geom_bar()

ggplot2提供了20多个统计变换以供你使用。每个统计变换都是一个函数，通过?stat_bin来获取帮助。如果想要查看全部的统计变换，可以使用ggplot2速查表。

位置探索

用fill和color为图形属性上色

ggplot(data = diamonds, aes(x=cut, fill=clarity)) + geom_bar()

这种堆叠是由position参数设定的位置调整功能自动完成的，如果不想生成堆叠条形图，有以下三种选项供选择：

  • identity
  • fill
  • dodge

position = “fill”的效果与堆叠相似，但每组堆叠条形具有同样的高度，可以轻松看出各组间的比例

ggplot(data = diamonds, aes(x=cut, fill=clarity)) + geom_bar(position = "fill")

position = “dodge”将每组条形依次并列放置

使用 position = position_dodge() 可以设置分离的宽度，默认 position_dodge(width=0.9)

ggplot(data = diamonds, aes(x=cut, fill=clarity)) + geom_bar(position = "dodge")
ggplot(data = mpg, aes(x = displ, y=hwy)) + 
 geom_point(position = 'jitter')

使用dplyr进行数据转换

dplyr基础

suppressMessages(library(nycflights13))
Warning message:
"package 'nycflights13' was built under R version 4.0.5"

head(flights,2)
A tibble: 2 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
201311517515283081911UA1545N14228EWRIAH22714005152013-01-01 05:00:00
201311533529485083020UA1714N24211LGAIAH22714165292013-01-01 05:00:00

5个dplyr核心函数：

这些函数都可以和group_by()函数联合起来使用，group_by函数可以改变以上每个函数的作用范围，让其从在整个数据集上操作变为在每个分组上分别操作。

工作方式：
1）第一个参数是数据框
2）随后的参数使用变量名称（不带引号）描述了在数据框上进行的操作
3）输出结果是一个新的数据框

使用filter()筛选行

filter()函数可以基于观测值筛选出一个观测子集。第一个参数是数据框名称，第二个及以后的参数是用来筛选数据框的表达式。

nrow(filter(flights, month == 1, day == 1))

842

dplyr不修改输入，所以想要得到新的数据框需要重新赋值。

jan1 <- filter(flights, month == 1, day == 1)

R要么输出结果，要么将结果保存在一个变量中。如果想同时完成这两项操作，可以用括号将赋值语句括起来。

head((jan1 <- filter(flights, month == 1, day == 1)))
A tibble: 6 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
201311517515 2 830 819 11UA1545N14228EWRIAH22714005152013-01-01 05:00:00
201311533529 4 850 830 20UA1714N24211LGAIAH22714165292013-01-01 05:00:00
201311542540 2 923 850 33AA1141N619AAJFKMIA16010895402013-01-01 05:00:00
201311544545-110041022-18B6 725N804JBJFKBQN18315765452013-01-01 05:00:00
201311554600-6 812 837-25DL 461N668DNLGAATL116 7626 02013-01-01 06:00:00
201311554558-4 740 728 12UA1696N39463EWRORD150 7195582013-01-01 05:00:00

比较运算符

= <= == !=

逻辑运算符

与（&）、或（|）、非（！）

缺失值

对于缺失值NA的操作，一般结果也是NA

如果想要确定一个值是否为缺失值，可以使用is.na()函数

使用arrange()排列行

arrange()函数与filter()函数工作方式相似，只是功能上前者是改变行的顺序

多个排列条件，靠前的条件，优先级越高

head(
    arrange(flights, year, month, day)
,1)
A tibble: 1 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
201311517515283081911UA1545N14228EWRIAH22714005152013-01-01 05:00:00

使用desc()可以进行降序排列

head(
    arrange(flights, desc(arr_delay))
    ,1)
A tibble: 1 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
2013196419001301124215301272HA51N384HAJFKHNL6404983902013-01-09 09:00:00

缺失值总是排在最后

使用select()选择列

通过基于变量名的操作，select()函数可诶让你快速生成一个有用的变量集

head(
    select(flights, year, month, day)
    ,2)
A tibble: 2 × 3
yearmonthday
<int><int><int>
201311
201311
 
head(
    select(flights, year:day)
    ,2)
A tibble: 2 × 3
yearmonthday
<int><int><int>
201311
201311

用-进行反向选择

head(
    select(flights, -(year:day))
,2)
A tibble: 2 × 16
dep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
517515283081911UA1545N14228EWRIAH22714005152013-01-01 05:00:00
533529485083020UA1714N24211LGAIAH22714165292013-01-01 05:00:00

还可以在select()函数中使用一些辅助函数

  • starts_with()
  • ends_with()
  • contains()
  • matches(“(.)\1”):选择匹配正则表达式的变量
head(
    rename(flights, tail_num = tailnum)
    ,2)
A tibble: 2 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttail_numorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
201311517515283081911UA1545N14228EWRIAH22714005152013-01-01 05:00:00
201311533529485083020UA1714N24211LGAIAH22714165292013-01-01 05:00:00

另一种用法是将select()函数和everything()辅助函数结合起来使用。当想要将几个变量移到开头时，这种用法非常奏效

head(
    select(flights, time_hour, air_time, everything())
    ,2)
A tibble: 2 × 19
time_hourair_timeyearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestdistancehourminute
<dttm><dbl><int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl>
2013-01-01 05:00:00227201311517515283081911UA1545N14228EWRIAH1400515
2013-01-01 05:00:00227201311533529485083020UA1714N24211LGAIAH1416529

select 使用的注意事项

MASS 也有 select 函数，并且优先级更高，哪怕再次载入 tidyverse也不行，载入 dplyr 也不行。
解决办法：

  • 绝对引用函数
    > dplyr::select
  • 放到环境变量中
    > select = dplyr::select

推荐载入包时，将下面代码放在开头，就可以肆无忌惮的应用select了，毕竟，环境变量的优先级是第一位的：

> library(tidyverse)
> select = dplyr::select

参考 R select

使用mutate()添加新变量

除了选择现有的列，我们还经常需要添加新列，新列是现有列的函数。这就是mutate()函数的作用

flights_sml <- select(flights,
                     year:day,
                     ends_with('delay'),
                     distance,
                     air_time)
head(flights_sml,2)
A tibble: 2 × 7
yearmonthdaydep_delayarr_delaydistanceair_time
<int><int><int><dbl><dbl><dbl><dbl>
2013112111400227
2013114201416227
 
mutate(flights_sml,
      gain = arr_delay - dep_delay,
      speed = distance / air_time * 60)
A tibble: 336776 × 9
yearmonthdaydep_delayarr_delaydistanceair_timegainspeed
<int><int><int><dbl><dbl><dbl><dbl><dbl><dbl>
201311 2 111400227 9370.0441
201311 4 201416227 16374.2731
201311 2 331089160 31408.3750
201311-1-181576183-17516.7213
201311-6-25 762116-19394.1379
201311-4 12 719150 16287.6000
201311-5 191065158 24404.4304
201311-3-14 229 53-11259.2453
201311-3 -8 944140 -5404.5714
201311-2 8 733138 10318.6957
201311-2 -21028149 0413.9597
201311-2 -31005158 -1381.6456
201311-2 72475345 9430.4348
201311-2-142565361-12426.3158
201311-1 311389257 32324.2802
201311 0 -4 187 44 -4255.0000
201311-1 -82227337 -7396.4985
201311 0 -71076152 -7424.7368
201311 0 12 762134 12341.1940
201311 1 -61023147 -7417.5510
201311-8 -81020170 0360.0000
201311-3 16 502105 19286.8571
201311-4-121085152 -8428.2895
201311-4 -8 760128 -4356.2500
201311 0-171085157-17414.6497
201311 8 32 719139 24310.3597
20131111 142586366 3423.9344
201311 3 41074175 1368.2286
201311 0-211598182-21526.8132
201311 0 -9 746120 -9373.0000
⋮⋮⋮⋮⋮⋮⋮⋮⋮
2013930 -2-24 305 45-22406.6667
2013930 -2 -9 529 72 -7440.8333
2013930 -2-311626213-29458.0282
2013930 30 -2 292 45-32389.3333
2013930 -9-30 213 36-21355.0000
2013930 0-302475298-30498.3221
2013930 13 11 284 47 -2362.5532
2013930 0-251598192-25499.3750
2013930 10 31076139 -7464.4604
2013930 -7-23 200 37-16324.3243
2013930 -9-16 209 39 -7321.5385
2013930194194 301 50 0361.2000
2013930 -2 8 378 61 10371.8033
2013930 27 7 764 97-20472.5773
2013930 72 57 872120-15436.0000
2013930-14-21 273 48 -7341.2500
2013930 80 422565318-38483.9623
2013930154130 944123-24460.4878
2013930 -8 -8 266 43 0371.1628
2013930 -5-17 209 41-12305.8537
2013930-10-20 301 52-10347.3077
2013930 -5-16 264 47-11337.0213
2013930 12 1 187 33-11340.0000
2013930-10-251617196-15495.0000
2013930 NA NA 764 NA NA NA
2013930 NA NA 213 NA NA NA
2013930 NA NA 198 NA NA NA
2013930 NA NA 764 NA NA NA
2013930 NA NA 419 NA NA NA
2013930 NA NA 431 NA NA NA

一旦创建，新列就可以立即使用

head(
    mutate(flights_sml,
          gain = arr_delay - dep_delay,
          hours = air_time / 60,
          gain_per_hour = gain / hours),
    2)
A tibble: 2 × 10
yearmonthdaydep_delayarr_delaydistanceair_timegainhoursgain_per_hour
<int><int><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
2013112111400227 93.7833332.378855
2013114201416227163.7833334.229075

如果只想保留新变量，可以使用transmute()函数

head(
    transmute(flights_sml,
          gain = arr_delay - dep_delay,
          hours = air_time / 60,
          gain_per_hour = gain / hours),
    2)
A tibble: 2 × 3
gainhoursgain_per_hour
<dbl><dbl><dbl>
93.7833332.378855
163.7833334.229075

使用summarize()进行分组摘要

summarize()，可以将数据框折叠成一行：

summarize(flights, delay = mean(dep_delay, na.rm = T))
A tibble: 1 × 1
delay
<dbl>
12.63907

如果不与 group_by() 一起使用，那么 summarize() 也就没什么大用，group_by() 可以将分析单位从整个数据集更改为单个分组。
在分组后的数据框上使用 dplyr 函数时，它们会自动地应用到每个分组。 

by_day <- group_by(flights, year, month, day)
delay_everyday <- summarize(by_day, delay = mean(dep_delay, na.rm = T))
head(delay_everyday)
`summarise()` has grouped output by 'year', 'month'. You can override using the `.groups` argument.
A grouped_df: 6 × 4
yearmonthdaydelay
<int><int><int><dbl>
20131111.548926
20131213.858824
20131310.987832
201314 8.951595
201315 5.732218
201316 7.148014

✏ group_by() 和 summarize() 的组合构成了使用 dplyr 包时最常用的操作之一：分组摘要。

在此之前，先介绍一个功能强大的新概念：管道

使用管道组合多种操作

假设我们想要研究每个目的地的距离和平均延误时间之间的关系。一般会写出以下代码：

by_dest <- group_by(flights, dest)
delay <- summarise(by_dest,
                   count = n(),
                   dist = mean(distance, na.rm = TRUE),
                   delay = mean(arr_delay, na.rm = TRUE)
                  )

delay <- filter(delay, count > 20, dest != 'NHL')

完成数据准备需要三步：

  1. 按照目的地对航安进行分组
  2. 进行摘要统计，计算距离、平均延误时间和航班数量
  3. 通过筛选除去噪声点

在这个过程中，我们不得不对不关心的中间变量进行命名，这会影响我们的分析速度

所以就引出了管道，%>%：

head(flights,2)
A tibble: 2 × 19
yearmonthdaydep_timesched_dep_timedep_delayarr_timesched_arr_timearr_delaycarrierflighttailnumorigindestair_timedistancehourminutetime_hour
<int><int><int><int><int><dbl><int><int><dbl><chr><int><chr><chr><chr><dbl><dbl><dbl><dbl><dttm>
201311517515283081911UA1545N14228EWRIAH22714005152013-01-01 05:00:00
201311533529485083020UA1714N24211LGAIAH22714165292013-01-01 05:00:00
 
delays <- flights %>%
    group_by(dest) %>% 
    summarize(
        count = n(),
        dist = mean(distance, na.rm = T),
        delay = mean(arr_delay, na.rm = T)
    ) %>%
    filter(count > 20, dest != 'HNL')

这种方法的重点在于转换的过程，而不是转换的对象，这使得代码具有更好的可读性。
将上述读作一串命令式语句：分组，然后摘要，然后进行筛选。 阅读时%>%可读作然后。

使用这种方法时，x %>% y 会转换为 f(x, y), x %>% f(y) %>% g(z) 会转换为 g(f(x, y), z),以此类推。

如果你熟悉 Linux 命令行的话，上述看起来就像：

cat flights | sort | awk ''

支持管道操作是 tidyverse 中的R包的核心原则之一。唯一的例外就是 ggplot2 ：ggplot2 的下一个版本 ggvis 支持管道操作

缺失值(NA)

聚合函数遵循缺失值的一般规则，如果输入中有缺失值，则输出也是缺失值。
好在，所有聚合函数都有一个 na.rm 参数，它可以在计算前，排除 NA。

is.na() 用来判断是否是缺失值

not_cancelled <- flights %>% 
    filter(!is.na(dep_delay), !is.na(arr_delay))

mean_delay <- not_cancelled %>% 
    group_by(year,month,day) %>%
    summarize(mean = mean(dep_delay))
`summarise()` has grouped output by 'year', 'month'. You can override using the `.groups` argument.

计数, n()

mean_delay个计数（n()）或非缺失值计数（sum(!is.na())）

例如，查看一下具有最长平均延误时间的飞机（通过机尾编号进行识别）：

delays <- not_cancelled %>%
    group_by(tailnum) %>%
    summarize(delay = mean(arr_delay, na.rm = T),
             n = n())

ggplot(data = delays, mapping = aes(x=n, y = delay)) + 
 geom_point(alpha = 1/10)

筛掉样本数很小的分组。

下面展示了将 ggplot2 集成到 dplyr 工作流的一种有效方式。

delays %>%
    filter(n>25) %>%
    ggplot(mapping = aes(x=n, y=delay)) + 
        geom_point(alpha = 1/10)

常用的摘要函数

除了 均值mean()、计数n()、求和sum()，R中还提供了很多其他常用的摘要函数

  • 位置度量
    • median()
    • 均方误差（标准误差, standard deviation，sd）是分散程度的标准度量方式
    • 四分位距IQR(),基本等价于中位差mad(x)
  • 秩的度量
    • min(x)
    • max(x)
    • quantile(x, 0.25)
  • 定位度量
    • first()
    • nth(x, 2)
    • last(x)
  • 计数
    • n()
    • count()
  • 逻辑值的计数和比例
    • TRUE 会 转换为1， FALSE会转换为 0.
    • sum(x>10)，找出 x 中TRUE的数量
    • mean(y==0)，找出 y 中为了 0 的比例

按照多个变量分组

取消分组

如果想要取消分组，并回到为分组的数据继续操作，那么可以使用ungroup()函数：

daily <- group_by(flights, year, month, day)
daily %>%
    ungroup() %>%   # 不再按日期分组
    summarize(flights = n()) # 所有航班
A tibble: 1 × 1
flights
<int>
336776

工作流：脚本

这里有几点写R脚本的建议：

  1. 建议将所用的包放在脚本开头，这样就很容易知道需要什么包了
  2. 不要在分享的脚本中包含 install.packages() 或 setwd()函数，这种改变别人计算机设置的行为会引起众怒。

探索性数据分析

探索性数据分析(exploratory data analysis, EDA)。

这一部分，我就不再复刻代码了。想表达的就是快速实现各种可能的想法，通过考察结果，不断迭代来进一步过滤、找到更进一步的问题，整理成文，与他人交流。

数据清洗只是EDA的一项具体应用，此时你提出的问题是，数据是否符合预期。想要使用所有的EDA工具：可视化、数据转换和建模。

经常数据的以下几方面特征：

  • 数据分布
  • 典型值
  • 异常值
  • 缺失值
  • 相关变动
    • 变量变量与连续变量
    • 两个分类变量
    • 两个连续变量
  • 模式和模型
    • 线性
    • 非线性
文章作者: 梁绍波
文章链接: https://liangshaw.github.io/2022/01/11/RdataScience-part1/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 梁绍波 !
R data science handbook
评论
  目录