# 第 39 章 探索性数据分析-anscombe数据集

?datasets::anscombe

Four x-y datasets which have the same traditional statistical properties (mean, variance, correlation, regression line, etc.), yet are quite different.

d <- datasets::anscombe
head(d)
##   x1 x2 x3 x4   y1   y2    y3   y4
## 1 10 10 10  8 8.04 9.14  7.46 6.58
## 2  8  8  8  8 6.95 8.14  6.77 5.76
## 3 13 13 13  8 7.58 8.74 12.74 7.71
## 4  9  9  9  8 8.81 8.77  7.11 8.84
## 5 11 11 11  8 8.33 9.26  7.81 8.47
## 6 14 14 14  8 9.96 8.10  8.84 7.04

## 39.1 探索anscombe

library(tidyverse)

• 规整数据
• 分组统计
• 建模
• 可视化

## 39.2 规整数据

head(d)
##   x1 x2 x3 x4   y1   y2    y3   y4
## 1 10 10 10  8 8.04 9.14  7.46 6.58
## 2  8  8  8  8 6.95 8.14  6.77 5.76
## 3 13 13 13  8 7.58 8.74 12.74 7.71
## 4  9  9  9  8 8.81 8.77  7.11 8.84
## 5 11 11 11  8 8.33 9.26  7.81 8.47
## 6 14 14 14  8 9.96 8.10  8.84 7.04

d %>%
ggplot(aes(x = x, y = y)) +
geom_point() +
facet_wrap(~set)

set x y
1 10 8.04
1 8 6.95
2 10 9.14
2 8 8.14

### 39.2.1 小小的回顾

dt <- tibble(id = c("a", "b"), x_1 = 1:2, x_2 = 3:4, y_1 = 5:6, y_2 = 8:9)
dt
## # A tibble: 2 x 5
##   id      x_1   x_2   y_1   y_2
##   <chr> <int> <int> <int> <int>
## 1 a         1     3     5     8
## 2 b         2     4     6     9
dt %>% pivot_longer(-id,
names_to = "name",
values_to = "vaules"
)
## # A tibble: 8 x 3
##   id    name  vaules
##   <chr> <chr>  <int>
## 1 a     x_1        1
## 2 a     x_2        3
## 3 a     y_1        5
## 4 a     y_2        8
## 5 b     x_1        2
## 6 b     x_2        4
## 7 b     y_1        6
## 8 b     y_2        9

dt %>% pivot_longer(
cols = -id,
names_to = "name",
names_pattern = "(.)_.",
values_to = "vaules"
)
## # A tibble: 8 x 3
##   id    name  vaules
##   <chr> <chr>  <int>
## 1 a     x          1
## 2 a     x          3
## 3 a     y          5
## 4 a     y          8
## 5 b     x          2
## 6 b     x          4
## 7 b     y          6
## 8 b     y          9

dt %>% pivot_longer(
cols = -id,
names_to = "name",
names_pattern = "._(.)",
values_to = "vaules"
)
## # A tibble: 8 x 3
##   id    name  vaules
##   <chr> <chr>  <int>
## 1 a     1          1
## 2 a     2          3
## 3 a     1          5
## 4 a     2          8
## 5 b     1          2
## 6 b     2          4
## 7 b     1          6
## 8 b     2          9

dt %>% pivot_longer(
cols = -id,
names_to = c("name", "group"),
names_pattern = "(.)_(.)",
values_to = "vaules"
)
## # A tibble: 8 x 4
##   id    name  group vaules
##   <chr> <chr> <chr>  <int>
## 1 a     x     1          1
## 2 a     x     2          3
## 3 a     y     1          5
## 4 a     y     2          8
## 5 b     x     1          2
## 6 b     x     2          4
## 7 b     y     1          6
## 8 b     y     2          9

dt %>% pivot_longer(
cols = -id,
names_to = c(".value", "group"),
names_pattern = "(.)_(.)",
values_to = "vaules"
)
## # A tibble: 4 x 4
##   id    group     x     y
##   <chr> <chr> <int> <int>
## 1 a     1         1     5
## 2 a     2         3     8
## 3 b     1         2     6
## 4 b     2         4     9

### 39.2.2 回到案例

knitr::include_graphics("images/pivot_longer_values.jpg")

tidy_d <- d %>%
pivot_longer(
cols = everything(),
names_to = c(".value", "set"),
names_pattern = "(.)(.)"
)
tidy_d
## # A tibble: 44 x 3
##    set       x     y
##    <chr> <dbl> <dbl>
##  1 1        10  8.04
##  2 2        10  9.14
##  3 3        10  7.46
##  4 4         8  6.58
##  5 1         8  6.95
##  6 2         8  8.14
##  7 3         8  6.77
##  8 4         8  5.76
##  9 1        13  7.58
## 10 2        13  8.74
## # ... with 34 more rows

• cols = everything() 表示选择所有列
• names_to = c(".value", "set") 希望变型后的列名是c(".value", "set"), 这里 ".value" 是个特殊的符号，代表着names_pattern匹配过来的值，一般情况下，是多个值，如果传给".value""x, y, z"，那么列名就会变成c("x", "y", "z", "set")
• names_pattern = "(.)(.)" 将变换前的列名按照指定的正则表达式匹配，并且传递给names_to的对应的参数，比如这里第一个(.)传递给.value；第二个(.)传递给set.

## 39.3 统计

tidy_d_summary <- tidy_d %>%
group_by(set) %>%
summarise(across(
.cols = everything(),
.fns = lst(mean, sd, var),
.names = "{col}_{fn}"
))
tidy_d_summary
## # A tibble: 4 x 7
##   set   x_mean  x_sd x_var y_mean  y_sd y_var
##   <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 1          9  3.32    11   7.50  2.03  4.13
## 2 2          9  3.32    11   7.50  2.03  4.13
## 3 3          9  3.32    11   7.5   2.03  4.12
## 4 4          9  3.32    11   7.50  2.03  4.12

## 39.4 建模

tidy_d %>%
group_nest(set) %>%
mutate(
fit = map(data, ~ lm(y ~ x, data = .x)),
tidy = map(fit, broom::tidy),
glance = map(fit, broom::glance)
) %>%
unnest(tidy)

tidy_d %>%
group_by(set) %>%
group_modify(
~ broom::tidy(lm(y ~ x, data = .))
)
## # A tibble: 8 x 6
## # Groups:   set [4]
##   set   term       estimate std.error statistic p.value
##   <chr> <chr>         <dbl>     <dbl>     <dbl>   <dbl>
## 1 1     (Intercep~    3.00      1.12       2.67 0.0257
## 2 1     x             0.500     0.118      4.24 0.00217
## 3 2     (Intercep~    3.00      1.13       2.67 0.0258
## 4 2     x             0.5       0.118      4.24 0.00218
## 5 3     (Intercep~    3.00      1.12       2.67 0.0256
## 6 3     x             0.500     0.118      4.24 0.00218
## 7 4     (Intercep~    3.00      1.12       2.67 0.0256
## 8 4     x             0.500     0.118      4.24 0.00216
tidy_d %>%
group_by(set) %>%
summarise(
broom::tidy(lm(y ~ x, data = cur_data()))
)
## # A tibble: 8 x 6
## # Groups:   set [4]
##   set   term       estimate std.error statistic p.value
##   <chr> <chr>         <dbl>     <dbl>     <dbl>   <dbl>
## 1 1     (Intercep~    3.00      1.12       2.67 0.0257
## 2 1     x             0.500     0.118      4.24 0.00217
## 3 2     (Intercep~    3.00      1.13       2.67 0.0258
## 4 2     x             0.5       0.118      4.24 0.00218
## 5 3     (Intercep~    3.00      1.12       2.67 0.0256
## 6 3     x             0.500     0.118      4.24 0.00218
## 7 4     (Intercep~    3.00      1.12       2.67 0.0256
## 8 4     x             0.500     0.118      4.24 0.00216

## 39.5 可视化看看

tidy_d %>%
ggplot(aes(x = x, y = y, colour = set)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme(legend.position = "none") +
facet_wrap(~set)