Processing math: 100%

Признаемся: часть сокровищ R мы утаили :)

Много пакетов R можно поставить с официального репозитория CRAN через Tools - Install packages. Однако несметные богатства пакетов лежат на github! Сайт rdrr.io позволяет искать сокровища :)

О новых прикольных пакетах обычно рассказывают на r-bloggers.com.

Чтобы поставить пакет нужно знать имя разработчика и название пакета. Мы поставим пакет patchwork от thomasp85. Томас придумал простой язык, чтобы рисовать несколько графиков рядом.

Ставим пакет Томаса:

devtools::install_github('thomasp85/patchwork')

Как правило, разработчики интригующе описывают пакет прямо на гитхабе. Например, описание пакета у Томаса.

Установка выполняется один раз :)

Теперь подключаем пакеты!

library(tidyverse) # обработка данных, графики...
## ── Attaching packages ──────────────────
## ✔ ggplot2 2.2.1.9000     ✔ purrr   0.2.5     
## ✔ tibble  1.4.2          ✔ dplyr   0.7.5     
## ✔ tidyr   0.8.1          ✔ stringr 1.3.1     
## ✔ readr   1.1.1          ✔ forcats 0.3.0
## Warning: package 'tibble' was built under R version 3.4.3
## Warning: package 'purrr' was built under R version 3.4.4
## Warning: package 'forcats' was built under R version 3.4.3
## ── Conflicts ── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(skimr)# описательные статистики
## Warning: package 'skimr' was built under R version 3.4.4
library(rio) # импорт фантастического количества форматов данных

library(cluster) # кластерный анализ
## Warning: package 'cluster' was built under R version 3.4.4
library(factoextra) # визуализации kmeans, pca,
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
library(dendextend) # визуализация дендрограмм
## 
## ---------------------
## Welcome to dendextend version 1.8.0
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## Or contact: <tal.galili@gmail.com>
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------
## 
## Attaching package: 'dendextend'
## The following object is masked from 'package:stats':
## 
##     cutree
library(corrplot) # визуализация корреляций
## corrplot 0.84 loaded
library(broom) # метла превращает результаты оценивания моделей в таблички
## Warning: package 'broom' was built under R version 3.4.4
library(naniar) # визуализация пропущенных значений
## Warning: package 'naniar' was built under R version 3.4.3
## 
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
## 
##     n_complete
library(visdat) # визуализация пропущенных значений

library(patchwork) # удобное расположение графиков рядом

library(nycflights13) # данные о полётах в Нью-Йорке

1 Основные классы данных

В R несколько разных классов данных:

1.0.1 Векторы

Создадим несколько векторов с пропущенными значениями!

a <- c(5, 6, 7, NA, 8)
b <- 11:15
b
## [1] 11 12 13 14 15
  • Что находится в векторе b?

С числовыми векторами работают школьные функции:

max(b)
## [1] 15
max(a)
## [1] NA
max(a, na.rm = TRUE)
## [1] 8
mean(a)
## [1] NA
mean(a, na.rm = TRUE)
## [1] 6.5
  • Чем отличаются функции max(a) и max(a, na.rm = TRUE)?

1.0.2 Таблички

Таблицы можно создавать самим c помощью функции tibble().

mx <- tibble(ID = c('A', 'B', 'C', 'D', 'E'),
           x = c(1, 2, 4, 1, 8))
mx
## # A tibble: 5 x 2
##   ID        x
##   <chr> <dbl>
## 1 A         1
## 2 B         2
## 3 C         4
## 4 D         1
## 5 E         8
  • Упражнение 1.

Создайте таблицу my, с двумя столбцами: ID и y, и пятью значениями в каждом. Столбец y можно заполнить числами произвольно.

my <- tibble(ID = c('A', 'C', 'D', 'F', 'E'),
             y = c(1, 2, -3, 4, -5))
my
## # A tibble: 5 x 2
##   ID        y
##   <chr> <dbl>
## 1 A         1
## 2 C         2
## 3 D        -3
## 4 F         4
## 5 E        -5

1.0.3 Матрицы

Матрицы очень похожи на таблицы. Их основное отличие в том, что в них могут быть только числа. С матрицами работают люди, которым нужна ну очень высокая скорость. Мы будем работать только с таблицами.

my_matrix <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2)
my_matrix
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
class(my_matrix)
## [1] "matrix"

А если нам вдруг встретится злобная, но ценная матрица, мы её тут же преобразуем в таблицу.

skim(my_matrix) # не работает :(
## No skim method exists for class matrix.
my_tibble <- as_tibble(my_matrix)
skim(my_tibble)
## Skim summary statistics
##  n obs: 2 
##  n variables: 3
## Warning: package 'bindrcpp' was built under R version 3.4.4
## 
## Variable type: numeric 
##  variable missing complete n mean   sd p0  p25 p50  p75 p100     hist
##        V1       0        2 2  1.5 0.71  1 1.25 1.5 1.75    2 ▇▁▁▁▁▁▁▇
##        V2       0        2 2  3.5 0.71  3 3.25 3.5 3.75    4 ▇▁▁▁▁▁▁▇
##        V3       0        2 2  5.5 0.71  5 5.25 5.5 5.75    6 ▇▁▁▁▁▁▁▇

1.0.4 Списки

Список - это чулан, в котором может быть полно ценных вещей. В списке может храниться всё что угодно. Как правило в списках хранятся результаты оценивания сложных моделей.

my_list <- list(kreks = 'Привет!', peks = c(5, 6, NA), feks = diamonds)
my_list$kreks
## [1] "Привет!"
my_list$peks
## [1]  5  6 NA
attributes(my_list)
## $names
## [1] "kreks" "peks"  "feks"

Узнать, что лежит в чулане-списке, можно либо с помощью функции attributes, либо нажав на значок лупы во вкладке Environment в Rstudio. Обычно вкладка Environment справа вверху.

2 Соединение таблиц

Таблицы можно соединять! Существует четыре способа сделать это: - left_join(x, y, by = 'colname') сохраняет все наблюдения в таблице x - right_join(x, y, by = 'colname') сохраняет все наблюдения в таблице y - full_join(x, y, by = 'colname') сохраняет наблюдения из обеих таблиц - inner_join(x, y, by = 'colname') объединяет только те строки обеих таблиц, в которых нет пропущенных наблюдений

Экспериментируем!

left_join(mx, my, by = 'ID')
## # A tibble: 5 x 3
##   ID        x     y
##   <chr> <dbl> <dbl>
## 1 A         1     1
## 2 B         2    NA
## 3 C         4     2
## 4 D         1    -3
## 5 E         8    -5
right_join(mx, my, by = 'ID')
## # A tibble: 5 x 3
##   ID        x     y
##   <chr> <dbl> <dbl>
## 1 A         1     1
## 2 C         4     2
## 3 D         1    -3
## 4 F        NA     4
## 5 E         8    -5
inner_join(mx, my, by = 'ID')
## # A tibble: 4 x 3
##   ID        x     y
##   <chr> <dbl> <dbl>
## 1 A         1     1
## 2 C         4     2
## 3 D         1    -3
## 4 E         8    -5
full_join(mx, my, by = 'ID')
## # A tibble: 6 x 3
##   ID        x     y
##   <chr> <dbl> <dbl>
## 1 A         1     1
## 2 B         2    NA
## 3 C         4     2
## 4 D         1    -3
## 5 E         8    -5
## 6 F        NA     4

Нужно объединить таблицы flights и weather о вылетах из Нью-Йорка сразу по нескольким столбцам!

glimpse(flights)
## Observations: 336,776
## Variables: 19
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2018c.1.0/
## zoneinfo/Europe/Moscow'
## $ year           <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,...
## $ month          <int> 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,...
## $ dep_time       <int> 517, 533, 542, 544, 554, 554, 555, 557, 557, 55...
## $ sched_dep_time <int> 515, 529, 540, 545, 600, 558, 600, 600, 600, 60...
## $ dep_delay      <dbl> 2, 4, 2, -1, -6, -4, -5, -3, -3, -2, -2, -2, -2...
## $ arr_time       <int> 830, 850, 923, 1004, 812, 740, 913, 709, 838, 7...
## $ sched_arr_time <int> 819, 830, 850, 1022, 837, 728, 854, 723, 846, 7...
## $ arr_delay      <dbl> 11, 20, 33, -18, -25, 12, 19, -14, -8, 8, -2, -...
## $ carrier        <chr> "UA", "UA", "AA", "B6", "DL", "UA", "B6", "EV",...
## $ flight         <int> 1545, 1714, 1141, 725, 461, 1696, 507, 5708, 79...
## $ tailnum        <chr> "N14228", "N24211", "N619AA", "N804JB", "N668DN...
## $ origin         <chr> "EWR", "LGA", "JFK", "JFK", "LGA", "EWR", "EWR"...
## $ dest           <chr> "IAH", "IAH", "MIA", "BQN", "ATL", "ORD", "FLL"...
## $ air_time       <dbl> 227, 227, 160, 183, 116, 150, 158, 53, 140, 138...
## $ distance       <dbl> 1400, 1416, 1089, 1576, 762, 719, 1065, 229, 94...
## $ hour           <dbl> 5, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5,...
## $ minute         <dbl> 15, 29, 40, 45, 0, 58, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ time_hour      <dttm> 2013-01-01 05:00:00, 2013-01-01 05:00:00, 2013...
glimpse(weather)
## Observations: 26,130
## Variables: 15
## $ origin     <chr> "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "E...
## $ year       <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 201...
## $ month      <dbl> 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, ...
## $ hour       <int> 0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...
## $ temp       <dbl> 37.04, 37.04, 37.94, 37.94, 37.94, 39.02, 39.02, 39...
## $ dewp       <dbl> 21.92, 21.92, 21.92, 23.00, 24.08, 26.06, 26.96, 28...
## $ humid      <dbl> 53.97, 53.97, 52.09, 54.51, 57.04, 59.37, 61.63, 64...
## $ wind_dir   <dbl> 230, 230, 230, 230, 240, 270, 250, 240, 250, 260, 2...
## $ wind_speed <dbl> 10.35702, 13.80936, 12.65858, 13.80936, 14.96014, 1...
## $ wind_gust  <dbl> 11.918651, 15.891535, 14.567241, 15.891535, 17.2158...
## $ precip     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ pressure   <dbl> 1013.9, 1013.0, 1012.6, 1012.7, 1012.8, 1012.0, 101...
## $ visib      <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,...
## $ time_hour  <dttm> 2013-01-01 00:00:00, 2013-01-01 01:00:00, 2013-01-...
left <- left_join(flights, weather, by = c('year', 'month', 'day', 'hour', 'origin'))
glimpse(left)
## Observations: 336,776
## Variables: 29
## $ year           <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,...
## $ month          <dbl> 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,...
## $ dep_time       <int> 517, 533, 542, 544, 554, 554, 555, 557, 557, 55...
## $ sched_dep_time <int> 515, 529, 540, 545, 600, 558, 600, 600, 600, 60...
## $ dep_delay      <dbl> 2, 4, 2, -1, -6, -4, -5, -3, -3, -2, -2, -2, -2...
## $ arr_time       <int> 830, 850, 923, 1004, 812, 740, 913, 709, 838, 7...
## $ sched_arr_time <int> 819, 830, 850, 1022, 837, 728, 854, 723, 846, 7...
## $ arr_delay      <dbl> 11, 20, 33, -18, -25, 12, 19, -14, -8, 8, -2, -...
## $ carrier        <chr> "UA", "UA", "AA", "B6", "DL", "UA", "B6", "EV",...
## $ flight         <int> 1545, 1714, 1141, 725, 461, 1696, 507, 5708, 79...
## $ tailnum        <chr> "N14228", "N24211", "N619AA", "N804JB", "N668DN...
## $ origin         <chr> "EWR", "LGA", "JFK", "JFK", "LGA", "EWR", "EWR"...
## $ dest           <chr> "IAH", "IAH", "MIA", "BQN", "ATL", "ORD", "FLL"...
## $ air_time       <dbl> 227, 227, 160, 183, 116, 150, 158, 53, 140, 138...
## $ distance       <dbl> 1400, 1416, 1089, 1576, 762, 719, 1065, 229, 94...
## $ hour           <dbl> 5, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5,...
## $ minute         <dbl> 15, 29, 40, 45, 0, 58, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ time_hour.x    <dttm> 2013-01-01 05:00:00, 2013-01-01 05:00:00, 2013...
## $ temp           <dbl> NA, NA, NA, NA, 39.92, NA, 39.02, 39.92, 39.02,...
## $ dewp           <dbl> NA, NA, NA, NA, 26.06, NA, 26.06, 26.06, 26.06,...
## $ humid          <dbl> NA, NA, NA, NA, 57.33, NA, 59.37, 57.33, 59.37,...
## $ wind_dir       <dbl> NA, NA, NA, NA, 260, NA, 270, 260, 260, 260, 26...
## $ wind_speed     <dbl> NA, NA, NA, NA, 13.80936, NA, 10.35702, 13.8093...
## $ wind_gust      <dbl> NA, NA, NA, NA, 15.89154, NA, 11.91865, 15.8915...
## $ precip         <dbl> NA, NA, NA, NA, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ pressure       <dbl> NA, NA, NA, NA, 1011.9, NA, 1012.0, 1011.9, 101...
## $ visib          <dbl> NA, NA, NA, NA, 10, NA, 10, 10, 10, 10, 10, 10,...
## $ time_hour.y    <dttm> NA, NA, NA, NA, 2013-01-01 06:00:00, NA, 2013-...
inner <- inner_join(flights, weather, by = c('year', 'month', 'day', 'hour', 'origin'))
glimpse(inner)
## Observations: 335,563
## Variables: 29
## $ year           <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,...
## $ month          <dbl> 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,...
## $ dep_time       <int> 554, 555, 557, 557, 558, 558, 558, 558, 558, 55...
## $ sched_dep_time <int> 600, 600, 600, 600, 600, 600, 600, 600, 600, 60...
## $ dep_delay      <dbl> -6, -5, -3, -3, -2, -2, -2, -2, -2, -1, -1, 0, ...
## $ arr_time       <int> 812, 913, 709, 838, 753, 849, 853, 924, 923, 94...
## $ sched_arr_time <int> 837, 854, 723, 846, 745, 851, 856, 917, 937, 91...
## $ arr_delay      <dbl> -25, 19, -14, -8, 8, -2, -3, 7, -14, 31, -8, -7...
## $ carrier        <chr> "DL", "B6", "EV", "B6", "AA", "B6", "B6", "UA",...
## $ flight         <int> 461, 507, 5708, 79, 301, 49, 71, 194, 1124, 707...
## $ tailnum        <chr> "N668DN", "N516JB", "N829AS", "N593JB", "N3ALAA...
## $ origin         <chr> "LGA", "EWR", "LGA", "JFK", "LGA", "JFK", "JFK"...
## $ dest           <chr> "ATL", "FLL", "IAD", "MCO", "ORD", "PBI", "TPA"...
## $ air_time       <dbl> 116, 158, 53, 140, 138, 149, 158, 345, 361, 257...
## $ distance       <dbl> 762, 1065, 229, 944, 733, 1028, 1005, 2475, 256...
## $ hour           <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,...
## $ minute         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 5...
## $ time_hour.x    <dttm> 2013-01-01 06:00:00, 2013-01-01 06:00:00, 2013...
## $ temp           <dbl> 39.92, 39.02, 39.92, 39.02, 39.92, 39.02, 39.02...
## $ dewp           <dbl> 26.06, 26.06, 26.06, 26.06, 26.06, 26.06, 26.06...
## $ humid          <dbl> 57.33, 59.37, 57.33, 59.37, 57.33, 59.37, 59.37...
## $ wind_dir       <dbl> 260, 270, 260, 260, 260, 260, 260, 260, 270, 26...
## $ wind_speed     <dbl> 13.80936, 10.35702, 13.80936, 12.65858, 13.8093...
## $ wind_gust      <dbl> 15.89154, 11.91865, 15.89154, 14.56724, 15.8915...
## $ precip         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ pressure       <dbl> 1011.9, 1012.0, 1011.9, 1012.6, 1011.9, 1012.6,...
## $ visib          <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,...
## $ time_hour.y    <dttm> 2013-01-01 06:00:00, 2013-01-01 06:00:00, 2013...

3 Длинные и широкие таблицы

Таблицы бывают длинными и широкими! При поступлении нового наблюдения длинная растёт в длину, а широкая может и в ширину :)

Возьмём квартальные данные Росстата о ВВП по источникам доходов.

gdp <- import('gdp.xls')
glimpse(gdp)
## Observations: 4
## Variables: 29
## $ X__1             <chr> "Валовой внутренний продукт", "оплата труда н...
## $ `1 квартал 2011` <dbl> 13028.844, 6055.775, 1807.421, 5165.648
## $ `2 квартал 2011` <dbl> 14481.142, 6469.798, 2205.521, 5805.824
## $ `3 квартал 2011` <dbl> 15805.582, 6678.376, 2280.732, 6846.474
## $ `4 квартал 2011` <dbl> 16966.973, 7182.727, 2453.331, 7330.951
## $ `1 квартал 2012` <dbl> 15182.987, 6991.065, 2194.436, 5997.486
## $ `2 квартал 2012` <dbl> 16472.215, 7448.678, 2582.311, 6441.226
## $ `3 квартал 2012` <dbl> 17733.545, 7615.371, 2423.306, 7694.848
## $ `4 квартал 2012` <dbl> 18775.136, 8146.048, 2630.635, 7998.463
## $ `1 квартал 2013` <dbl> 16375.261, 7828.690, 2335.080, 6211.491
## $ `2 квартал 2013` <dbl> 17538.839, 8350.152, 2572.417, 6616.269
## $ `3 квартал 2013` <dbl> 19058.114, 8569.038, 2522.510, 7966.566
## $ `4 квартал 2013` <dbl> 20161.681, 9044.402, 2632.313, 8484.967
## $ `1 квартал 2014` <dbl> 17390.202, 8637.304, 2497.615, 6255.283
## $ `2 квартал 2014` <dbl> 19127.958, 9261.309, 2879.702, 6986.947
## $ `3 квартал 2014` <dbl> 20758.591, 9463.617, 2712.466, 8582.508
## $ `4 квартал 2014` <dbl> 21922.908, 10024.947, 2914.458, 8983.503
## $ `1 квартал 2015` <dbl> 18579.1, 9102.5, 2416.8, 7059.8
## $ `2 квартал 2015` <dbl> 19865.7, 9689.5, 2163.7, 8012.5
## $ `3 квартал 2015` <dbl> 22017.6, 9877.0, 2397.9, 9742.7
## $ `4 квартал 2015` <dbl> 22924.8, 10440.8, 2292.1, 10191.9
## $ `1 квартал 2016` <dbl> 19110.1, 9543.9, 2243.0, 7323.2
## $ `2 квартал 2016` <dbl> 20624.9, 10224.9, 2208.8, 8191.2
## $ `3 квартал 2016` <dbl> 22515.8, 10422.7, 2403.1, 9690.0
## $ `4 квартал 2016` <dbl> 23897.8, 11085.0, 2566.0, 10246.8
## $ `1 квартал 2017` <dbl> 20549.8, 10288.2, 2254.7, 8006.9
## $ `2 квартал 2017` <dbl> 22035.1, 11012.3, 2393.0, 8629.8
## $ `3 квартал 2017` <dbl> 23948.9, 11147.0, 2519.2, 10282.7
## $ `4 квартал 2017` <dbl> 25503.4, 11815.7, 2787.9, 10899.8
gdp
##                                                                                                                         X__1
## 1                                                                                                 Валовой внутренний продукт
## 2 оплата труда наемных работников (включая оплату труда и смешанные доходы, не наблюдаемые прямыми статистическими методами)
## 3                                                                                     чистые налоги на производство и импорт
## 4                                                                       валовая прибыль экономики и валовые смешанные доходы
##   1 квартал 2011 2 квартал 2011 3 квартал 2011 4 квартал 2011
## 1      13028.844      14481.142      15805.582      16966.973
## 2       6055.775       6469.798       6678.376       7182.727
## 3       1807.421       2205.521       2280.732       2453.331
## 4       5165.648       5805.824       6846.474       7330.951
##   1 квартал 2012 2 квартал 2012 3 квартал 2012 4 квартал 2012
## 1      15182.987      16472.215      17733.545      18775.136
## 2       6991.065       7448.678       7615.371       8146.048
## 3       2194.436       2582.311       2423.306       2630.635
## 4       5997.486       6441.226       7694.848       7998.463
##   1 квартал 2013 2 квартал 2013 3 квартал 2013 4 квартал 2013
## 1      16375.261      17538.839      19058.114      20161.681
## 2       7828.690       8350.152       8569.038       9044.402
## 3       2335.080       2572.417       2522.510       2632.313
## 4       6211.491       6616.269       7966.566       8484.967
##   1 квартал 2014 2 квартал 2014 3 квартал 2014 4 квартал 2014
## 1      17390.202      19127.958      20758.591      21922.908
## 2       8637.304       9261.309       9463.617      10024.947
## 3       2497.615       2879.702       2712.466       2914.458
## 4       6255.283       6986.947       8582.508       8983.503
##   1 квартал 2015 2 квартал 2015 3 квартал 2015 4 квартал 2015
## 1        18579.1        19865.7        22017.6        22924.8
## 2         9102.5         9689.5         9877.0        10440.8
## 3         2416.8         2163.7         2397.9         2292.1
## 4         7059.8         8012.5         9742.7        10191.9
##   1 квартал 2016 2 квартал 2016 3 квартал 2016 4 квартал 2016
## 1        19110.1        20624.9        22515.8        23897.8
## 2         9543.9        10224.9        10422.7        11085.0
## 3         2243.0         2208.8         2403.1         2566.0
## 4         7323.2         8191.2         9690.0        10246.8
##   1 квартал 2017 2 квартал 2017 3 квартал 2017 4 квартал 2017
## 1        20549.8        22035.1        23948.9        25503.4
## 2        10288.2        11012.3        11147.0        11815.7
## 3         2254.7         2393.0         2519.2         2787.9
## 4         8006.9         8629.8        10282.7        10899.8

Обнаружив неудачное название переменной переименуем её сразу!

gdp <- rename(gdp, indicator = X__1)

Чтобы превратить широкую таблицу в длинную, нужна функция gather() из пакета tidyr.

В длинной таблице мы создадим два новых столбца: - столбец date с названиями старых столбцов - столбец quantity со значениями таблицы

Если столбец indicator нужно оставить как есть, упоминаем его с минусом!

long_gdp <- gather(gdp, key = 'date', value = 'quantity', -indicator)
glimpse(long_gdp)
## Observations: 112
## Variables: 3
## $ indicator <chr> "Валовой внутренний продукт", "оплата труда наемных ...
## $ date      <chr> "1 квартал 2011", "1 квартал 2011", "1 квартал 2011"...
## $ quantity  <dbl> 13028.844, 6055.775, 1807.421, 5165.648, 14481.142, ...

Создадим типичную длинную таблицу! Например, для каждого качества огранки cut и цвета color посчитаем медиану цены:

diamonds_med <- group_by(diamonds, cut, color) %>%
  summarise(med_price = median(price))
diamonds_med
## # A tibble: 35 x 3
## # Groups:   cut [?]
##    cut   color med_price
##    <ord> <ord>     <dbl>
##  1 Fair  D         3730 
##  2 Fair  E         2956 
##  3 Fair  F         3035 
##  4 Fair  G         3057 
##  5 Fair  H         3816 
##  6 Fair  I         3246 
##  7 Fair  J         3302 
##  8 Good  D         2728.
##  9 Good  E         2420 
## 10 Good  F         2647 
## # ... with 25 more rows

А теперь превратим её в широкую:

diamonds_wide <- spread(diamonds_med, key = 'color', value = 'med_price')
diamonds_wide
## # A tibble: 5 x 8
## # Groups:   cut [5]
##   cut           D     E     F     G     H     I     J
##   <ord>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Fair      3730  2956   3035 3057  3816  3246   3302
## 2 Good      2728. 2420   2647 3340  3468. 3640.  3733
## 3 Very Good 2310  1990.  2471 2437  3734  3888   4113
## 4 Premium   2009  1928   2841 2745  4511  4640   5063
## 5 Ideal     1576  1437   1775 1858. 2278  2659   4096

Здесь тупому компьютеру умный Homo Sapiens укажет: - какая переменная станет столбцами новой таблицы: color - какая переменная пойдёт внутрь таблицы: med_price

Не упомянутые переменные останутся строками новой широкой таблицы.

Из данных про воздух в Нью-Йорке airquality сделайте длинную таблицу. В качестве сохраняемых столбцов возьмите переменные Month и Day.

air_long <- gather(airquality, key = 'indicator', value = 'value', -Month, -Day)
air_long
##     Month Day indicator value
## 1       5   1     Ozone  41.0
## 2       5   2     Ozone  36.0
## 3       5   3     Ozone  12.0
## 4       5   4     Ozone  18.0
## 5       5   5     Ozone    NA
## 6       5   6     Ozone  28.0
## 7       5   7     Ozone  23.0
## 8       5   8     Ozone  19.0
## 9       5   9     Ozone   8.0
## 10      5  10     Ozone    NA
## 11      5  11     Ozone   7.0
## 12      5  12     Ozone  16.0
## 13      5  13     Ozone  11.0
## 14      5  14     Ozone  14.0
## 15      5  15     Ozone  18.0
## 16      5  16     Ozone  14.0
## 17      5  17     Ozone  34.0
## 18      5  18     Ozone   6.0
## 19      5  19     Ozone  30.0
## 20      5  20     Ozone  11.0
## 21      5  21     Ozone   1.0
## 22      5  22     Ozone  11.0
## 23      5  23     Ozone   4.0
## 24      5  24     Ozone  32.0
## 25      5  25     Ozone    NA
## 26      5  26     Ozone    NA
## 27      5  27     Ozone    NA
## 28      5  28     Ozone  23.0
## 29      5  29     Ozone  45.0
## 30      5  30     Ozone 115.0
## 31      5  31     Ozone  37.0
## 32      6   1     Ozone    NA
## 33      6   2     Ozone    NA
## 34      6   3     Ozone    NA
## 35      6   4     Ozone    NA
## 36      6   5     Ozone    NA
## 37      6   6     Ozone    NA
## 38      6   7     Ozone  29.0
## 39      6   8     Ozone    NA
## 40      6   9     Ozone  71.0
## 41      6  10     Ozone  39.0
## 42      6  11     Ozone    NA
## 43      6  12     Ozone    NA
## 44      6  13     Ozone  23.0
## 45      6  14     Ozone    NA
## 46      6  15     Ozone    NA
## 47      6  16     Ozone  21.0
## 48      6  17     Ozone  37.0
## 49      6  18     Ozone  20.0
## 50      6  19     Ozone  12.0
## 51      6  20     Ozone  13.0
## 52      6  21     Ozone    NA
## 53      6  22     Ozone    NA
## 54      6  23     Ozone    NA
## 55      6  24     Ozone    NA
## 56      6  25     Ozone    NA
## 57      6  26     Ozone    NA
## 58      6  27     Ozone    NA
## 59      6  28     Ozone    NA
## 60      6  29     Ozone    NA
## 61      6  30     Ozone    NA
## 62      7   1     Ozone 135.0
## 63      7   2     Ozone  49.0
## 64      7   3     Ozone  32.0
## 65      7   4     Ozone    NA
## 66      7   5     Ozone  64.0
## 67      7   6     Ozone  40.0
## 68      7   7     Ozone  77.0
## 69      7   8     Ozone  97.0
## 70      7   9     Ozone  97.0
## 71      7  10     Ozone  85.0
## 72      7  11     Ozone    NA
## 73      7  12     Ozone  10.0
## 74      7  13     Ozone  27.0
## 75      7  14     Ozone    NA
## 76      7  15     Ozone   7.0
## 77      7  16     Ozone  48.0
## 78      7  17     Ozone  35.0
## 79      7  18     Ozone  61.0
## 80      7  19     Ozone  79.0
## 81      7  20     Ozone  63.0
## 82      7  21     Ozone  16.0
## 83      7  22     Ozone    NA
## 84      7  23     Ozone    NA
## 85      7  24     Ozone  80.0
## 86      7  25     Ozone 108.0
## 87      7  26     Ozone  20.0
## 88      7  27     Ozone  52.0
## 89      7  28     Ozone  82.0
## 90      7  29     Ozone  50.0
## 91      7  30     Ozone  64.0
## 92      7  31     Ozone  59.0
## 93      8   1     Ozone  39.0
## 94      8   2     Ozone   9.0
## 95      8   3     Ozone  16.0
## 96      8   4     Ozone  78.0
## 97      8   5     Ozone  35.0
## 98      8   6     Ozone  66.0
## 99      8   7     Ozone 122.0
## 100     8   8     Ozone  89.0
## 101     8   9     Ozone 110.0
## 102     8  10     Ozone    NA
## 103     8  11     Ozone    NA
## 104     8  12     Ozone  44.0
## 105     8  13     Ozone  28.0
## 106     8  14     Ozone  65.0
## 107     8  15     Ozone    NA
## 108     8  16     Ozone  22.0
## 109     8  17     Ozone  59.0
## 110     8  18     Ozone  23.0
## 111     8  19     Ozone  31.0
## 112     8  20     Ozone  44.0
## 113     8  21     Ozone  21.0
## 114     8  22     Ozone   9.0
## 115     8  23     Ozone    NA
## 116     8  24     Ozone  45.0
## 117     8  25     Ozone 168.0
## 118     8  26     Ozone  73.0
## 119     8  27     Ozone    NA
## 120     8  28     Ozone  76.0
## 121     8  29     Ozone 118.0
## 122     8  30     Ozone  84.0
## 123     8  31     Ozone  85.0
## 124     9   1     Ozone  96.0
## 125     9   2     Ozone  78.0
## 126     9   3     Ozone  73.0
## 127     9   4     Ozone  91.0
## 128     9   5     Ozone  47.0
## 129     9   6     Ozone  32.0
## 130     9   7     Ozone  20.0
## 131     9   8     Ozone  23.0
## 132     9   9     Ozone  21.0
## 133     9  10     Ozone  24.0
## 134     9  11     Ozone  44.0
## 135     9  12     Ozone  21.0
## 136     9  13     Ozone  28.0
## 137     9  14     Ozone   9.0
## 138     9  15     Ozone  13.0
## 139     9  16     Ozone  46.0
## 140     9  17     Ozone  18.0
## 141     9  18     Ozone  13.0
## 142     9  19     Ozone  24.0
## 143     9  20     Ozone  16.0
## 144     9  21     Ozone  13.0
## 145     9  22     Ozone  23.0
## 146     9  23     Ozone  36.0
## 147     9  24     Ozone   7.0
## 148     9  25     Ozone  14.0
## 149     9  26     Ozone  30.0
## 150     9  27     Ozone    NA
## 151     9  28     Ozone  14.0
## 152     9  29     Ozone  18.0
## 153     9  30     Ozone  20.0
## 154     5   1   Solar.R 190.0
## 155     5   2   Solar.R 118.0
## 156     5   3   Solar.R 149.0
## 157     5   4   Solar.R 313.0
## 158     5   5   Solar.R    NA
## 159     5   6   Solar.R    NA
## 160     5   7   Solar.R 299.0
## 161     5   8   Solar.R  99.0
## 162     5   9   Solar.R  19.0
## 163     5  10   Solar.R 194.0
## 164     5  11   Solar.R    NA
## 165     5  12   Solar.R 256.0
## 166     5  13   Solar.R 290.0
## 167     5  14   Solar.R 274.0
## 168     5  15   Solar.R  65.0
## 169     5  16   Solar.R 334.0
## 170     5  17   Solar.R 307.0
## 171     5  18   Solar.R  78.0
## 172     5  19   Solar.R 322.0
## 173     5  20   Solar.R  44.0
## 174     5  21   Solar.R   8.0
## 175     5  22   Solar.R 320.0
## 176     5  23   Solar.R  25.0
## 177     5  24   Solar.R  92.0
## 178     5  25   Solar.R  66.0
## 179     5  26   Solar.R 266.0
## 180     5  27   Solar.R    NA
## 181     5  28   Solar.R  13.0
## 182     5  29   Solar.R 252.0
## 183     5  30   Solar.R 223.0
## 184     5  31   Solar.R 279.0
## 185     6   1   Solar.R 286.0
## 186     6   2   Solar.R 287.0
## 187     6   3   Solar.R 242.0
## 188     6   4   Solar.R 186.0
## 189     6   5   Solar.R 220.0
## 190     6   6   Solar.R 264.0
## 191     6   7   Solar.R 127.0
## 192     6   8   Solar.R 273.0
## 193     6   9   Solar.R 291.0
## 194     6  10   Solar.R 323.0
## 195     6  11   Solar.R 259.0
## 196     6  12   Solar.R 250.0
## 197     6  13   Solar.R 148.0
## 198     6  14   Solar.R 332.0
## 199     6  15   Solar.R 322.0
## 200     6  16   Solar.R 191.0
## 201     6  17   Solar.R 284.0
## 202     6  18   Solar.R  37.0
## 203     6  19   Solar.R 120.0
## 204     6  20   Solar.R 137.0
## 205     6  21   Solar.R 150.0
## 206     6  22   Solar.R  59.0
## 207     6  23   Solar.R  91.0
## 208     6  24   Solar.R 250.0
## 209     6  25   Solar.R 135.0
## 210     6  26   Solar.R 127.0
## 211     6  27   Solar.R  47.0
## 212     6  28   Solar.R  98.0
## 213     6  29   Solar.R  31.0
## 214     6  30   Solar.R 138.0
## 215     7   1   Solar.R 269.0
## 216     7   2   Solar.R 248.0
## 217     7   3   Solar.R 236.0
## 218     7   4   Solar.R 101.0
## 219     7   5   Solar.R 175.0
## 220     7   6   Solar.R 314.0
## 221     7   7   Solar.R 276.0
## 222     7   8   Solar.R 267.0
## 223     7   9   Solar.R 272.0
## 224     7  10   Solar.R 175.0
## 225     7  11   Solar.R 139.0
## 226     7  12   Solar.R 264.0
## 227     7  13   Solar.R 175.0
## 228     7  14   Solar.R 291.0
## 229     7  15   Solar.R  48.0
## 230     7  16   Solar.R 260.0
## 231     7  17   Solar.R 274.0
## 232     7  18   Solar.R 285.0
## 233     7  19   Solar.R 187.0
## 234     7  20   Solar.R 220.0
## 235     7  21   Solar.R   7.0
## 236     7  22   Solar.R 258.0
## 237     7  23   Solar.R 295.0
## 238     7  24   Solar.R 294.0
## 239     7  25   Solar.R 223.0
## 240     7  26   Solar.R  81.0
## 241     7  27   Solar.R  82.0
## 242     7  28   Solar.R 213.0
## 243     7  29   Solar.R 275.0
## 244     7  30   Solar.R 253.0
## 245     7  31   Solar.R 254.0
## 246     8   1   Solar.R  83.0
## 247     8   2   Solar.R  24.0
## 248     8   3   Solar.R  77.0
## 249     8   4   Solar.R    NA
## 250     8   5   Solar.R    NA
## 251     8   6   Solar.R    NA
## 252     8   7   Solar.R 255.0
## 253     8   8   Solar.R 229.0
## 254     8   9   Solar.R 207.0
## 255     8  10   Solar.R 222.0
## 256     8  11   Solar.R 137.0
## 257     8  12   Solar.R 192.0
## 258     8  13   Solar.R 273.0
## 259     8  14   Solar.R 157.0
## 260     8  15   Solar.R  64.0
## 261     8  16   Solar.R  71.0
## 262     8  17   Solar.R  51.0
## 263     8  18   Solar.R 115.0
## 264     8  19   Solar.R 244.0
## 265     8  20   Solar.R 190.0
## 266     8  21   Solar.R 259.0
## 267     8  22   Solar.R  36.0
## 268     8  23   Solar.R 255.0
## 269     8  24   Solar.R 212.0
## 270     8  25   Solar.R 238.0
## 271     8  26   Solar.R 215.0
## 272     8  27   Solar.R 153.0
## 273     8  28   Solar.R 203.0
## 274     8  29   Solar.R 225.0
## 275     8  30   Solar.R 237.0
## 276     8  31   Solar.R 188.0
## 277     9   1   Solar.R 167.0
## 278     9   2   Solar.R 197.0
## 279     9   3   Solar.R 183.0
## 280     9   4   Solar.R 189.0
## 281     9   5   Solar.R  95.0
## 282     9   6   Solar.R  92.0
## 283     9   7   Solar.R 252.0
## 284     9   8   Solar.R 220.0
## 285     9   9   Solar.R 230.0
## 286     9  10   Solar.R 259.0
## 287     9  11   Solar.R 236.0
## 288     9  12   Solar.R 259.0
## 289     9  13   Solar.R 238.0
## 290     9  14   Solar.R  24.0
## 291     9  15   Solar.R 112.0
## 292     9  16   Solar.R 237.0
## 293     9  17   Solar.R 224.0
## 294     9  18   Solar.R  27.0
## 295     9  19   Solar.R 238.0
## 296     9  20   Solar.R 201.0
## 297     9  21   Solar.R 238.0
## 298     9  22   Solar.R  14.0
## 299     9  23   Solar.R 139.0
## 300     9  24   Solar.R  49.0
## 301     9  25   Solar.R  20.0
## 302     9  26   Solar.R 193.0
## 303     9  27   Solar.R 145.0
## 304     9  28   Solar.R 191.0
## 305     9  29   Solar.R 131.0
## 306     9  30   Solar.R 223.0
## 307     5   1      Wind   7.4
## 308     5   2      Wind   8.0
## 309     5   3      Wind  12.6
## 310     5   4      Wind  11.5
## 311     5   5      Wind  14.3
## 312     5   6      Wind  14.9
## 313     5   7      Wind   8.6
## 314     5   8      Wind  13.8
## 315     5   9      Wind  20.1
## 316     5  10      Wind   8.6
## 317     5  11      Wind   6.9
## 318     5  12      Wind   9.7
## 319     5  13      Wind   9.2
## 320     5  14      Wind  10.9
## 321     5  15      Wind  13.2
## 322     5  16      Wind  11.5
## 323     5  17      Wind  12.0
## 324     5  18      Wind  18.4
## 325     5  19      Wind  11.5
## 326     5  20      Wind   9.7
## 327     5  21      Wind   9.7
## 328     5  22      Wind  16.6
## 329     5  23      Wind   9.7
## 330     5  24      Wind  12.0
## 331     5  25      Wind  16.6
## 332     5  26      Wind  14.9
## 333     5  27      Wind   8.0
## 334     5  28      Wind  12.0
## 335     5  29      Wind  14.9
## 336     5  30      Wind   5.7
## 337     5  31      Wind   7.4
## 338     6   1      Wind   8.6
## 339     6   2      Wind   9.7
## 340     6   3      Wind  16.1
## 341     6   4      Wind   9.2
## 342     6   5      Wind   8.6
## 343     6   6      Wind  14.3
## 344     6   7      Wind   9.7
## 345     6   8      Wind   6.9
## 346     6   9      Wind  13.8
## 347     6  10      Wind  11.5
## 348     6  11      Wind  10.9
## 349     6  12      Wind   9.2
## 350     6  13      Wind   8.0
## 351     6  14      Wind  13.8
## 352     6  15      Wind  11.5
## 353     6  16      Wind  14.9
## 354     6  17      Wind  20.7
## 355     6  18      Wind   9.2
## 356     6  19      Wind  11.5
## 357     6  20      Wind  10.3
## 358     6  21      Wind   6.3
## 359     6  22      Wind   1.7
## 360     6  23      Wind   4.6
## 361     6  24      Wind   6.3
## 362     6  25      Wind   8.0
## 363     6  26      Wind   8.0
## 364     6  27      Wind  10.3
## 365     6  28      Wind  11.5
## 366     6  29      Wind  14.9
## 367     6  30      Wind   8.0
## 368     7   1      Wind   4.1
## 369     7   2      Wind   9.2
## 370     7   3      Wind   9.2
## 371     7   4      Wind  10.9
## 372     7   5      Wind   4.6
## 373     7   6      Wind  10.9
## 374     7   7      Wind   5.1
## 375     7   8      Wind   6.3
## 376     7   9      Wind   5.7
## 377     7  10      Wind   7.4
## 378     7  11      Wind   8.6
## 379     7  12      Wind  14.3
## 380     7  13      Wind  14.9
## 381     7  14      Wind  14.9
## 382     7  15      Wind  14.3
## 383     7  16      Wind   6.9
## 384     7  17      Wind  10.3
## 385     7  18      Wind   6.3
## 386     7  19      Wind   5.1
## 387     7  20      Wind  11.5
## 388     7  21      Wind   6.9
## 389     7  22      Wind   9.7
## 390     7  23      Wind  11.5
## 391     7  24      Wind   8.6
## 392     7  25      Wind   8.0
## 393     7  26      Wind   8.6
## 394     7  27      Wind  12.0
## 395     7  28      Wind   7.4
## 396     7  29      Wind   7.4
## 397     7  30      Wind   7.4
## 398     7  31      Wind   9.2
## 399     8   1      Wind   6.9
## 400     8   2      Wind  13.8
## 401     8   3      Wind   7.4
## 402     8   4      Wind   6.9
## 403     8   5      Wind   7.4
## 404     8   6      Wind   4.6
## 405     8   7      Wind   4.0
## 406     8   8      Wind  10.3
## 407     8   9      Wind   8.0
## 408     8  10      Wind   8.6
## 409     8  11      Wind  11.5
## 410     8  12      Wind  11.5
## 411     8  13      Wind  11.5
## 412     8  14      Wind   9.7
## 413     8  15      Wind  11.5
## 414     8  16      Wind  10.3
## 415     8  17      Wind   6.3
## 416     8  18      Wind   7.4
## 417     8  19      Wind  10.9
## 418     8  20      Wind  10.3
## 419     8  21      Wind  15.5
## 420     8  22      Wind  14.3
## 421     8  23      Wind  12.6
## 422     8  24      Wind   9.7
## 423     8  25      Wind   3.4
## 424     8  26      Wind   8.0
## 425     8  27      Wind   5.7
## 426     8  28      Wind   9.7
## 427     8  29      Wind   2.3
## 428     8  30      Wind   6.3
## 429     8  31      Wind   6.3
## 430     9   1      Wind   6.9
## 431     9   2      Wind   5.1
## 432     9   3      Wind   2.8
## 433     9   4      Wind   4.6
## 434     9   5      Wind   7.4
## 435     9   6      Wind  15.5
## 436     9   7      Wind  10.9
## 437     9   8      Wind  10.3
## 438     9   9      Wind  10.9
## 439     9  10      Wind   9.7
## 440     9  11      Wind  14.9
## 441     9  12      Wind  15.5
## 442     9  13      Wind   6.3
## 443     9  14      Wind  10.9
## 444     9  15      Wind  11.5
## 445     9  16      Wind   6.9
## 446     9  17      Wind  13.8
## 447     9  18      Wind  10.3
## 448     9  19      Wind  10.3
## 449     9  20      Wind   8.0
## 450     9  21      Wind  12.6
## 451     9  22      Wind   9.2
## 452     9  23      Wind  10.3
## 453     9  24      Wind  10.3
## 454     9  25      Wind  16.6
## 455     9  26      Wind   6.9
## 456     9  27      Wind  13.2
## 457     9  28      Wind  14.3
## 458     9  29      Wind   8.0
## 459     9  30      Wind  11.5
## 460     5   1      Temp  67.0
## 461     5   2      Temp  72.0
## 462     5   3      Temp  74.0
## 463     5   4      Temp  62.0
## 464     5   5      Temp  56.0
## 465     5   6      Temp  66.0
## 466     5   7      Temp  65.0
## 467     5   8      Temp  59.0
## 468     5   9      Temp  61.0
## 469     5  10      Temp  69.0
## 470     5  11      Temp  74.0
## 471     5  12      Temp  69.0
## 472     5  13      Temp  66.0
## 473     5  14      Temp  68.0
## 474     5  15      Temp  58.0
## 475     5  16      Temp  64.0
## 476     5  17      Temp  66.0
## 477     5  18      Temp  57.0
## 478     5  19      Temp  68.0
## 479     5  20      Temp  62.0
## 480     5  21      Temp  59.0
## 481     5  22      Temp  73.0
## 482     5  23      Temp  61.0
## 483     5  24      Temp  61.0
## 484     5  25      Temp  57.0
## 485     5  26      Temp  58.0
## 486     5  27      Temp  57.0
## 487     5  28      Temp  67.0
## 488     5  29      Temp  81.0
## 489     5  30      Temp  79.0
## 490     5  31      Temp  76.0
## 491     6   1      Temp  78.0
## 492     6   2      Temp  74.0
## 493     6   3      Temp  67.0
## 494     6   4      Temp  84.0
## 495     6   5      Temp  85.0
## 496     6   6      Temp  79.0
## 497     6   7      Temp  82.0
## 498     6   8      Temp  87.0
## 499     6   9      Temp  90.0
## 500     6  10      Temp  87.0
## 501     6  11      Temp  93.0
## 502     6  12      Temp  92.0
## 503     6  13      Temp  82.0
## 504     6  14      Temp  80.0
## 505     6  15      Temp  79.0
## 506     6  16      Temp  77.0
## 507     6  17      Temp  72.0
## 508     6  18      Temp  65.0
## 509     6  19      Temp  73.0
## 510     6  20      Temp  76.0
## 511     6  21      Temp  77.0
## 512     6  22      Temp  76.0
## 513     6  23      Temp  76.0
## 514     6  24      Temp  76.0
## 515     6  25      Temp  75.0
## 516     6  26      Temp  78.0
## 517     6  27      Temp  73.0
## 518     6  28      Temp  80.0
## 519     6  29      Temp  77.0
## 520     6  30      Temp  83.0
## 521     7   1      Temp  84.0
## 522     7   2      Temp  85.0
## 523     7   3      Temp  81.0
## 524     7   4      Temp  84.0
## 525     7   5      Temp  83.0
## 526     7   6      Temp  83.0
## 527     7   7      Temp  88.0
## 528     7   8      Temp  92.0
## 529     7   9      Temp  92.0
## 530     7  10      Temp  89.0
## 531     7  11      Temp  82.0
## 532     7  12      Temp  73.0
## 533     7  13      Temp  81.0
## 534     7  14      Temp  91.0
## 535     7  15      Temp  80.0
## 536     7  16      Temp  81.0
## 537     7  17      Temp  82.0
## 538     7  18      Temp  84.0
## 539     7  19      Temp  87.0
## 540     7  20      Temp  85.0
## 541     7  21      Temp  74.0
## 542     7  22      Temp  81.0
## 543     7  23      Temp  82.0
## 544     7  24      Temp  86.0
## 545     7  25      Temp  85.0
## 546     7  26      Temp  82.0
## 547     7  27      Temp  86.0
## 548     7  28      Temp  88.0
## 549     7  29      Temp  86.0
## 550     7  30      Temp  83.0
## 551     7  31      Temp  81.0
## 552     8   1      Temp  81.0
## 553     8   2      Temp  81.0
## 554     8   3      Temp  82.0
## 555     8   4      Temp  86.0
## 556     8   5      Temp  85.0
## 557     8   6      Temp  87.0
## 558     8   7      Temp  89.0
## 559     8   8      Temp  90.0
## 560     8   9      Temp  90.0
## 561     8  10      Temp  92.0
## 562     8  11      Temp  86.0
## 563     8  12      Temp  86.0
## 564     8  13      Temp  82.0
## 565     8  14      Temp  80.0
## 566     8  15      Temp  79.0
## 567     8  16      Temp  77.0
## 568     8  17      Temp  79.0
## 569     8  18      Temp  76.0
## 570     8  19      Temp  78.0
## 571     8  20      Temp  78.0
## 572     8  21      Temp  77.0
## 573     8  22      Temp  72.0
## 574     8  23      Temp  75.0
## 575     8  24      Temp  79.0
## 576     8  25      Temp  81.0
## 577     8  26      Temp  86.0
## 578     8  27      Temp  88.0
## 579     8  28      Temp  97.0
## 580     8  29      Temp  94.0
## 581     8  30      Temp  96.0
## 582     8  31      Temp  94.0
## 583     9   1      Temp  91.0
## 584     9   2      Temp  92.0
## 585     9   3      Temp  93.0
## 586     9   4      Temp  93.0
## 587     9   5      Temp  87.0
## 588     9   6      Temp  84.0
## 589     9   7      Temp  80.0
## 590     9   8      Temp  78.0
## 591     9   9      Temp  75.0
## 592     9  10      Temp  73.0
## 593     9  11      Temp  81.0
## 594     9  12      Temp  76.0
## 595     9  13      Temp  77.0
## 596     9  14      Temp  71.0
## 597     9  15      Temp  71.0
## 598     9  16      Temp  78.0
## 599     9  17      Temp  67.0
## 600     9  18      Temp  76.0
## 601     9  19      Temp  68.0
## 602     9  20      Temp  82.0
## 603     9  21      Temp  64.0
## 604     9  22      Temp  71.0
## 605     9  23      Temp  81.0
## 606     9  24      Temp  69.0
## 607     9  25      Temp  63.0
## 608     9  26      Temp  70.0
## 609     9  27      Temp  77.0
## 610     9  28      Temp  75.0
## 611     9  29      Temp  76.0
## 612     9  30      Temp  68.0

А из получившейся длинной таблицы восстановите исходную!

air_wide <- spread(air_long, key = 'indicator', value = 'value')
air_wide
##     Month Day Ozone Solar.R Temp Wind
## 1       5   1    41     190   67  7.4
## 2       5   2    36     118   72  8.0
## 3       5   3    12     149   74 12.6
## 4       5   4    18     313   62 11.5
## 5       5   5    NA      NA   56 14.3
## 6       5   6    28      NA   66 14.9
## 7       5   7    23     299   65  8.6
## 8       5   8    19      99   59 13.8
## 9       5   9     8      19   61 20.1
## 10      5  10    NA     194   69  8.6
## 11      5  11     7      NA   74  6.9
## 12      5  12    16     256   69  9.7
## 13      5  13    11     290   66  9.2
## 14      5  14    14     274   68 10.9
## 15      5  15    18      65   58 13.2
## 16      5  16    14     334   64 11.5
## 17      5  17    34     307   66 12.0
## 18      5  18     6      78   57 18.4
## 19      5  19    30     322   68 11.5
## 20      5  20    11      44   62  9.7
## 21      5  21     1       8   59  9.7
## 22      5  22    11     320   73 16.6
## 23      5  23     4      25   61  9.7
## 24      5  24    32      92   61 12.0
## 25      5  25    NA      66   57 16.6
## 26      5  26    NA     266   58 14.9
## 27      5  27    NA      NA   57  8.0
## 28      5  28    23      13   67 12.0
## 29      5  29    45     252   81 14.9
## 30      5  30   115     223   79  5.7
## 31      5  31    37     279   76  7.4
## 32      6   1    NA     286   78  8.6
## 33      6   2    NA     287   74  9.7
## 34      6   3    NA     242   67 16.1
## 35      6   4    NA     186   84  9.2
## 36      6   5    NA     220   85  8.6
## 37      6   6    NA     264   79 14.3
## 38      6   7    29     127   82  9.7
## 39      6   8    NA     273   87  6.9
## 40      6   9    71     291   90 13.8
## 41      6  10    39     323   87 11.5
## 42      6  11    NA     259   93 10.9
## 43      6  12    NA     250   92  9.2
## 44      6  13    23     148   82  8.0
## 45      6  14    NA     332   80 13.8
## 46      6  15    NA     322   79 11.5
## 47      6  16    21     191   77 14.9
## 48      6  17    37     284   72 20.7
## 49      6  18    20      37   65  9.2
## 50      6  19    12     120   73 11.5
## 51      6  20    13     137   76 10.3
## 52      6  21    NA     150   77  6.3
## 53      6  22    NA      59   76  1.7
## 54      6  23    NA      91   76  4.6
## 55      6  24    NA     250   76  6.3
## 56      6  25    NA     135   75  8.0
## 57      6  26    NA     127   78  8.0
## 58      6  27    NA      47   73 10.3
## 59      6  28    NA      98   80 11.5
## 60      6  29    NA      31   77 14.9
## 61      6  30    NA     138   83  8.0
## 62      7   1   135     269   84  4.1
## 63      7   2    49     248   85  9.2
## 64      7   3    32     236   81  9.2
## 65      7   4    NA     101   84 10.9
## 66      7   5    64     175   83  4.6
## 67      7   6    40     314   83 10.9
## 68      7   7    77     276   88  5.1
## 69      7   8    97     267   92  6.3
## 70      7   9    97     272   92  5.7
## 71      7  10    85     175   89  7.4
## 72      7  11    NA     139   82  8.6
## 73      7  12    10     264   73 14.3
## 74      7  13    27     175   81 14.9
## 75      7  14    NA     291   91 14.9
## 76      7  15     7      48   80 14.3
## 77      7  16    48     260   81  6.9
## 78      7  17    35     274   82 10.3
## 79      7  18    61     285   84  6.3
## 80      7  19    79     187   87  5.1
## 81      7  20    63     220   85 11.5
## 82      7  21    16       7   74  6.9
## 83      7  22    NA     258   81  9.7
## 84      7  23    NA     295   82 11.5
## 85      7  24    80     294   86  8.6
## 86      7  25   108     223   85  8.0
## 87      7  26    20      81   82  8.6
## 88      7  27    52      82   86 12.0
## 89      7  28    82     213   88  7.4
## 90      7  29    50     275   86  7.4
## 91      7  30    64     253   83  7.4
## 92      7  31    59     254   81  9.2
## 93      8   1    39      83   81  6.9
## 94      8   2     9      24   81 13.8
## 95      8   3    16      77   82  7.4
## 96      8   4    78      NA   86  6.9
## 97      8   5    35      NA   85  7.4
## 98      8   6    66      NA   87  4.6
## 99      8   7   122     255   89  4.0
## 100     8   8    89     229   90 10.3
## 101     8   9   110     207   90  8.0
## 102     8  10    NA     222   92  8.6
## 103     8  11    NA     137   86 11.5
## 104     8  12    44     192   86 11.5
## 105     8  13    28     273   82 11.5
## 106     8  14    65     157   80  9.7
## 107     8  15    NA      64   79 11.5
## 108     8  16    22      71   77 10.3
## 109     8  17    59      51   79  6.3
## 110     8  18    23     115   76  7.4
## 111     8  19    31     244   78 10.9
## 112     8  20    44     190   78 10.3
## 113     8  21    21     259   77 15.5
## 114     8  22     9      36   72 14.3
## 115     8  23    NA     255   75 12.6
## 116     8  24    45     212   79  9.7
## 117     8  25   168     238   81  3.4
## 118     8  26    73     215   86  8.0
## 119     8  27    NA     153   88  5.7
## 120     8  28    76     203   97  9.7
## 121     8  29   118     225   94  2.3
## 122     8  30    84     237   96  6.3
## 123     8  31    85     188   94  6.3
## 124     9   1    96     167   91  6.9
## 125     9   2    78     197   92  5.1
## 126     9   3    73     183   93  2.8
## 127     9   4    91     189   93  4.6
## 128     9   5    47      95   87  7.4
## 129     9   6    32      92   84 15.5
## 130     9   7    20     252   80 10.9
## 131     9   8    23     220   78 10.3
## 132     9   9    21     230   75 10.9
## 133     9  10    24     259   73  9.7
## 134     9  11    44     236   81 14.9
## 135     9  12    21     259   76 15.5
## 136     9  13    28     238   77  6.3
## 137     9  14     9      24   71 10.9
## 138     9  15    13     112   71 11.5
## 139     9  16    46     237   78  6.9
## 140     9  17    18     224   67 13.8
## 141     9  18    13      27   76 10.3
## 142     9  19    24     238   68 10.3
## 143     9  20    16     201   82  8.0
## 144     9  21    13     238   64 12.6
## 145     9  22    23      14   71  9.2
## 146     9  23    36     139   81 10.3
## 147     9  24     7      49   69 10.3
## 148     9  25    14      20   63 16.6
## 149     9  26    30     193   70  6.9
## 150     9  27    NA     145   77 13.2
## 151     9  28    14     191   75 14.3
## 152     9  29    18     131   76  8.0
## 153     9  30    20     223   68 11.5

4 Работа с пропущенными данными :)

Команда skim() для описательных статистик показывает число пропущенных наблюдений в каждой перменной. Оказывается, их тоже можно визуализировать! Для этого будем использовать пакеты naniar и visdat. Посмотрим на пропущенные значения в сгенирированных данных typical_data командой vis_miss(). Сначала не будем добавлять аргументы, а затем сгруппируем пропуски с помощью иерархической кластеризации аргументом cluster = TRUE.

skim(typical_data)
## Skim summary statistics
##  n obs: 5000 
##  n variables: 9 
## 
## Variable type: character 
##  variable missing complete    n min max empty n_unique
##        ID    1000     4000 5000   4   4     0     4000
## 
## Variable type: factor 
##  variable missing complete    n n_unique
##      Race       0     5000 5000        8
##       Sex       0     5000 5000        2
##                               top_counts ordered
##  Whi: 3196, His: 811, Bla: 612, Asi: 244   FALSE
##              Mal: 2500, Fem: 2500, NA: 0   FALSE
## 
## Variable type: integer 
##  variable missing complete    n  mean   sd p0 p25 p50 p75 p100     hist
##       Age    1000     4000 5000 27.47 4.58 20  24  27  31   35 ▇▇▇▇▇▇▇▇
## 
## Variable type: logical 
##  variable missing complete    n mean                       count
##      Died       0     5000 5000 0.5  FAL: 2518, TRU: 2482, NA: 0
##    Smokes       0     5000 5000 0.18  FAL: 4088, TRU: 912, NA: 0
## 
## Variable type: numeric 
##    variable missing complete    n     mean       sd     p0      p25
##  Height(cm)       0     5000 5000   175.31     9.73 141.2    168.6 
##      Income    1000     4000 5000 40070.95 28533.8  403.77 19316.84
##          IQ    1000     4000 5000    99.76    10.09  63       93   
##       p50      p75      p100     hist
##    175.5    182.1     212    ▁▁▅▇▇▃▁▁
##  33128.29 53667.62 208210.62 ▇▇▃▂▁▁▁▁
##    100      107       135    ▁▁▃▇▇▅▁▁
vis_miss(typical_data)

vis_miss(typical_data, cluster = TRUE)

Визуализируйте пропуски во встроенных данных о воздухе в Нью-Йорке airquality, сгруппировав их по кластерам.

skim(airquality)
## Skim summary statistics
##  n obs: 153 
##  n variables: 6 
## 
## Variable type: integer 
##  variable missing complete   n   mean    sd p0    p25   p50    p75 p100
##       Day       0      153 153  15.8   8.86  1   8     16    23      31
##     Month       0      153 153   6.99  1.42  5   6      7     8       9
##     Ozone      37      116 153  42.13 32.99  1  18     31.5  63.25  168
##   Solar.R       7      146 153 185.93 90.06  7 115.75 205   258.75  334
##      Temp       0      153 153  77.88  9.47 56  72     79    85      97
##      hist
##  ▇▇▇▇▆▇▇▇
##  ▇▇▁▇▁▇▁▇
##  ▇▆▃▃▂▁▁▁
##  ▃▃▃▃▅▇▇▃
##  ▂▂▃▆▇▇▃▃
## 
## Variable type: numeric 
##  variable missing complete   n  mean   sd  p0 p25 p50  p75 p100     hist
##      Wind       0      153 153  9.96 3.52 1.7 7.4 9.7 11.5 20.7 ▁▃▇▇▅▅▁▁
vis_miss(airquality, cluster = TRUE)

Команда gg_miss_var() также позволяет посмотореть на количество пропущенных значений. По умолчанию пропуски считаются в абсолютных значениях, заказать проценты можно аргументом show_pct = TRUE.

gg_miss_var(typical_data)

gg_miss_var(typical_data, show_pct = TRUE) # число пропусков в процентах

Ещё больше красивых визуализаций — в виньетке пакета naniar.

Теперь, зная врага в лицо, приступим к борьбе с ним! Самый простой способ — удалить все строки данных, в которых есть пропуски. Сделаем это командой drop_na().

typical_removed <- drop_na(typical_data)
skim(typical_removed)
## Skim summary statistics
##  n obs: 2069 
##  n variables: 9 
## 
## Variable type: character 
##  variable missing complete    n min max empty n_unique
##        ID       0     2069 2069   4   4     0     2069
## 
## Variable type: factor 
##  variable missing complete    n n_unique
##      Race       0     2069 2069        8
##       Sex       0     2069 2069        2
##                              top_counts ordered
##  Whi: 1339, His: 336, Bla: 247, Asi: 94   FALSE
##             Fem: 1060, Mal: 1009, NA: 0   FALSE
## 
## Variable type: integer 
##  variable missing complete    n  mean   sd p0 p25 p50 p75 p100     hist
##       Age       0     2069 2069 27.46 4.62 20  23  27  32   35 ▇▇▇▇▇▇▇▇
## 
## Variable type: logical 
##  variable missing complete    n mean                       count
##      Died       0     2069 2069 0.49 FAL: 1060, TRU: 1009, NA: 0
##    Smokes       0     2069 2069 0.17  FAL: 1707, TRU: 362, NA: 0
## 
## Variable type: numeric 
##    variable missing complete    n     mean       sd     p0      p25
##  Height(cm)       0     2069 2069   175.07     9.61 141.5    168.5 
##      Income       0     2069 2069 39952.14 28482.09 403.77 18866.16
##          IQ       0     2069 2069    99.7     10.12  66       93   
##       p50      p75    p100     hist
##    175.2    181.9    208.1 ▁▁▃▇▇▅▁▁
##  33472.19 53274.31 2e+05   ▇▇▃▂▁▁▁▁
##    100      107      135   ▁▁▃▇▇▃▁▁

Команда replace_na из пакета tidyr позволяет заполнять пропуски любым числом, в том числе средним или медианой. Заполним пропущенные значения в переменной IQ минимальным значением, доход Income — медианным, ID — числом, а возраст Age — средним. Эти пожелания передадим аргументу replace списком.

typical_replaced <- replace_na(data = typical_data,
  replace = list(IQ =  min(typical_data$IQ, na.rm = TRUE),
  Income = median(typical_data$Income, na.rm = TRUE),
  ID = -1,
  Age = mean(typical_data$Age, na.rm = TRUE)))
skim(typical_replaced)
## Skim summary statistics
##  n obs: 5000 
##  n variables: 9 
## 
## Variable type: character 
##  variable missing complete    n min max empty n_unique
##        ID       0     5000 5000   2   4     0     4001
## 
## Variable type: factor 
##  variable missing complete    n n_unique
##      Race       0     5000 5000        8
##       Sex       0     5000 5000        2
##                               top_counts ordered
##  Whi: 3196, His: 811, Bla: 612, Asi: 244   FALSE
##              Mal: 2500, Fem: 2500, NA: 0   FALSE
## 
## Variable type: logical 
##  variable missing complete    n mean                       count
##      Died       0     5000 5000 0.5  FAL: 2518, TRU: 2482, NA: 0
##    Smokes       0     5000 5000 0.18  FAL: 4088, TRU: 912, NA: 0
## 
## Variable type: numeric 
##    variable missing complete    n     mean       sd     p0      p25
##         Age       0     5000 5000    27.47     4.09  20       24   
##  Height(cm)       0     5000 5000   175.31     9.73 141.2    168.6 
##      Income       0     5000 5000 38682.42 25671.45 403.77 22639.35
##          IQ       0     5000 5000    92.41    17.26  63       84   
##       p50      p75      p100     hist
##     27.47    30        35    ▂▃▃▇▂▂▃▂
##    175.5    182.1     212    ▁▁▅▇▇▃▁▁
##  33128.29 47762.38 208210.62 ▅▇▂▁▁▁▁▁
##     96      105       135    ▆▁▃▇▇▅▁▁
air_removed <- drop_na(airquality)
skim(air_removed)
## Skim summary statistics
##  n obs: 111 
##  n variables: 6 
## 
## Variable type: integer 
##  variable missing complete   n   mean    sd p0   p25 p50   p75 p100
##       Day       0      111 111  15.95  8.71  1   9    16  22.5   31
##     Month       0      111 111   7.22  1.47  5   6     7   9      9
##     Ozone       0      111 111  42.1  33.28  1  18    31  62    168
##   Solar.R       0      111 111 184.8  91.15  7 113.5 207 255.5  334
##      Temp       0      111 111  77.79  9.53 57  71    79  84.5   97
##      hist
##  ▇▆▆▇▇▇▆▇
##  ▆▂▁▇▁▆▁▇
##  ▇▆▃▂▂▁▁▁
##  ▅▃▃▂▆▇▇▃
##  ▂▃▃▅▇▅▃▂
## 
## Variable type: numeric 
##  variable missing complete   n mean   sd  p0 p25 p50  p75 p100     hist
##      Wind       0      111 111 9.94 3.56 2.3 7.4 9.7 11.5 20.7 ▂▂▇▇▂▂▁▁
air_replaced <- replace_na(data = airquality,
  replace = list(Ozone = mean(airquality$Ozone, na.rm = TRUE),
  Solar.R = median(airquality$Solar.R, na.rm = TRUE)))
skim(air_replaced)
## Skim summary statistics
##  n obs: 153 
##  n variables: 6 
## 
## Variable type: integer 
##  variable missing complete   n  mean   sd p0 p25 p50 p75 p100     hist
##       Day       0      153 153 15.8  8.86  1   8  16  23   31 ▇▇▇▇▆▇▇▇
##     Month       0      153 153  6.99 1.42  5   6   7   8    9 ▇▇▁▇▁▇▁▇
##      Temp       0      153 153 77.88 9.47 56  72  79  85   97 ▂▂▃▆▇▇▃▃
## 
## Variable type: numeric 
##  variable missing complete   n   mean    sd  p0   p25    p50   p75  p100
##     Ozone       0      153 153  42.13 28.69 1    21    42.13  46   168  
##   Solar.R       0      153 153 186.8  88.05 7   120   205    256   334  
##      Wind       0      153 153   9.96  3.52 1.7   7.4   9.7   11.5  20.7
##      hist
##  ▅▇▂▂▁▁▁▁
##  ▃▃▃▃▇▇▇▃
##  ▁▃▇▇▅▅▁▁

5 Кластеризация k-means

Для примера возьмём данные по потреблению протеинов Европе из книги Practial Machine Learning Cookbook. Сначала, как всегда, загрузим их и посмотрим описательные статистики.

protein <- import('Europenaprotein.csv')
skim(protein)
## Skim summary statistics
##  n obs: 25 
##  n variables: 10 
## 
## Variable type: character 
##  variable missing complete  n min max empty n_unique
##   Country       0       25 25   2  14     0       25
## 
## Variable type: numeric 
##   variable missing complete  n  mean    sd   p0  p25  p50  p75 p100
##    Cereals       0       25 25 32.25 10.97 18.6 24.3 28   40.1 56.7
##       Eggs       0       25 25  2.94  1.12  0.5  2.7  2.9  3.7  4.7
##       Fish       0       25 25  4.28  3.4   0.2  2.1  3.4  5.8 14.2
##     Fr&Veg       0       25 25  4.14  1.8   1.4  2.9  3.8  4.9  7.9
##       Milk       0       25 25 17.11  7.11  4.9 11.1 17.6 23.3 33.7
##       Nuts       0       25 25  3.07  1.99  0.7  1.5  2.4  4.7  7.8
##    RedMeat       0       25 25  9.83  3.35  4.4  7.8  9.5 10.6 18  
##     Starch       0       25 25  4.28  1.63  0.6  3.1  4.7  5.7  6.5
##  WhiteMeat       0       25 25  7.9   3.69  1.4  4.9  7.8 10.8 14  
##      hist
##  ▅▇▁▂▂▁▁▂
##  ▁▂▁▁▇▂▃▂
##  ▅▇▂▃▁▂▁▁
##  ▅▇▇▇▂▁▆▃
##  ▂▇▂▇▃▇▁▁
##  ▇▇▁▃▂▆▁▁
##  ▂▆▇▆▁▃▁▂
##  ▃▁▅▃▆▇▆▇
##  ▁▃▇▁▁▇▃▃

5.1 Масштабирование переменных

Для того, чтобы сравнивать переменные с разными единицами измерения, их масштабируют: вычитают среднее и делят на оценку стандартного отклонения.

x∗i=(xi−ˉx)/ˆσx, где ˉx=x1+x2+…+xnn, а ˆσx=(x1−ˉx)2+…+(xn−ˉx)2n−1.

Отмасштабируем данные с помощью встроенной функции scale(). Поскольку она может работать только с числами, первый столбец Country ей передавать не нужно. Результат сохраним в таблице protein_stand.

protein_stand <- mutate_if(protein, is.numeric, ~ as.vector(scale(.)))
skim(protein_stand)
## Skim summary statistics
##  n obs: 25 
##  n variables: 10 
## 
## Variable type: character 
##  variable missing complete  n min max empty n_unique
##   Country       0       25 25   2  14     0       25
## 
## Variable type: numeric 
##   variable missing complete  n     mean sd    p0   p25    p50  p75 p100
##    Cereals       0       25 25  2.6e-16  1 -1.24 -0.72 -0.39  0.72 2.23
##       Eggs       0       25 25  5.5e-17  1 -2.18 -0.21 -0.032 0.68 1.58
##       Fish       0       25 25  6.3e-17  1 -1.2  -0.64 -0.26  0.45 2.91
##     Fr&Veg       0       25 25 -6.6e-17  1 -1.52 -0.69 -0.19  0.42 2.09
##       Milk       0       25 25 -2.1e-16  1 -1.72 -0.85  0.069 0.87 2.33
##       Nuts       0       25 25 -1.1e-17  1 -1.19 -0.79 -0.34  0.82 2.38
##    RedMeat       0       25 25  1.8e-16  1 -1.62 -0.61 -0.098 0.23 2.44
##     Starch       0       25 25  1.7e-16  1 -2.25 -0.72  0.26  0.87 1.36
##  WhiteMeat       0       25 25  2.3e-17  1 -1.76 -0.81 -0.026 0.79 1.65
##      hist
##  ▅▇▁▂▂▁▁▂
##  ▁▂▁▁▇▂▃▂
##  ▅▇▂▃▁▂▁▁
##  ▅▇▇▇▂▁▆▃
##  ▂▇▂▇▃▇▁▁
##  ▇▇▁▃▂▆▁▁
##  ▂▆▇▆▁▃▁▂
##  ▃▁▅▃▆▇▆▇
##  ▁▃▇▁▁▇▃▃

Дополнение в виде функции as.vector нужно потому, что функция scale возвращает матрицу, а каждый столбец должен быть вектором :)

  • Посмотрите на описательные статистики. Что стало со средним значением и стандартным отклонением переменных?

  • Упражнение 8.

Хитрая исследовательница Даша хочет возвести в квадрат сразу все числовые столбцы в наборе данных diamonds. Помогите Даше возвести столбцы в квадрат!

diamonds2 <- mutate_if(diamonds, is.numeric, ~ .^2)

  • Упражнение 9.

  • Посмотрите на описательные статистики встроенного набора данных об арестах в США.
  • Отмасштабируйте данные и сохраните их в таблицу usa_stand.

  • Проверьте, что всё получилось :)

usa <- USArrests
skim(usa)
## Skim summary statistics
##  n obs: 50 
##  n variables: 4 
## 
## Variable type: integer 
##  variable missing complete  n   mean    sd p0   p25 p50    p75 p100
##   Assault       0       50 50 170.76 83.34 45 109   159 249     337
##  UrbanPop       0       50 50  65.54 14.47 32  54.5  66  77.75   91
##      hist
##  ▇▇▅▇▂▇▅▂
##  ▂▃▆▅▇▆▇▃
## 
## Variable type: numeric 
##  variable missing complete  n  mean   sd  p0   p25   p50   p75 p100
##    Murder       0       50 50  7.79 4.36 0.8  4.08  7.25 11.25 17.4
##      Rape       0       50 50 21.23 9.37 7.3 15.08 20.1  26.17 46  
##      hist
##  ▇▇▇▇▃▇▁▃
##  ▆▇▇▆▃▂▂▂
usa_stand <- mutate_if(usa , is.numeric, ~ as.vector(scale(.)))
glimpse(usa_stand)
## Observations: 50
## Variables: 4
## $ Murder   <dbl> 1.24256408, 0.50786248, 0.07163341, 0.23234938, 0.278...
## $ Assault  <dbl> 0.78283935, 1.10682252, 1.47880321, 0.23086801, 1.262...
## $ UrbanPop <dbl> -0.52090661, -1.21176419, 0.99898006, -1.07359268, 1....
## $ Rape     <dbl> -0.003416473, 2.484202941, 1.042878388, -0.184916602,...

Выполним кластеризацию методом k-средних с помощью функции kmeans. Название страны не используется для кластеризации, но нужно для меток на графиках. Поэтому мы уберем столбец Country из набора данных и превратим его в метки строк.

В качестве аргументов укажем отмасштабированные данные protein_no_country и количество кластеров centers. Пока мы не знаем, как выбирать оптимальное количество кластеров, поэтому предположим, что их три. Сохраним результат этого действия в список k_means_protein.

protein_no_country <- protein_stand %>% column_to_rownames(var = 'Country')
k_means_protein <- kmeans(protein_no_country, centers = 3)
k_means_protein
## K-means clustering with 3 clusters of sizes 4, 6, 15
## 
## Cluster means:
##      RedMeat  WhiteMeat       Eggs       Milk       Fish    Cereals
## 1 -0.5088020 -1.1088009 -0.4124850 -0.8320414  0.9819154  0.1300253
## 2 -0.7901419 -0.5267887 -1.1655757 -0.9047559 -0.9504683  1.4383272
## 3  0.4517373  0.5063957  0.5762263  0.5837801  0.1183432 -0.6100043
##       Starch       Nuts     Fr&Veg
## 1 -0.1842010  1.3108846  1.6292449
## 2 -0.7604664  0.8870168 -0.5373533
## 3  0.3533068 -0.7043759 -0.2195240
## 
## Clustering vector:
##        Albania        Austria        Belgium       Bulgaria Czechoslovakia 
##              2              3              3              2              3 
##        Denmark      E Germany        Finland         France         Greece 
##              3              3              3              3              1 
##        Hungary        Ireland          Italy    Netherlands         Norway 
##              2              3              1              3              3 
##         Poland       Portugal        Romania          Spain         Sweden 
##              3              1              2              1              3 
##    Switzerland             UK           USSR      W Germany     Yugoslavia 
##              3              3              2              3              2 
## 
## Within cluster sum of squares by cluster:
## [1] 18.92587 24.09113 62.96933
##  (between_SS / total_SS =  50.9 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

В списке k_means_protein лежит куча разной информации! Посмотрим на содержимое списка k_means_protein командой attributes().

attributes(k_means_protein)
## $names
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"      
## 
## $class
## [1] "kmeans"

Помните также, что в Rstudio можно кликнуть на лупу справа от k_means_protein во вкладке Environment.

Теперь, зная, что искать, мы можем посмотреть, например, на координаты центра кластеров или количество объектов в каждом из них.

k_means_protein$centers
##      RedMeat  WhiteMeat       Eggs       Milk       Fish    Cereals
## 1 -0.5088020 -1.1088009 -0.4124850 -0.8320414  0.9819154  0.1300253
## 2 -0.7901419 -0.5267887 -1.1655757 -0.9047559 -0.9504683  1.4383272
## 3  0.4517373  0.5063957  0.5762263  0.5837801  0.1183432 -0.6100043
##       Starch       Nuts     Fr&Veg
## 1 -0.1842010  1.3108846  1.6292449
## 2 -0.7604664  0.8870168 -0.5373533
## 3  0.3533068 -0.7043759 -0.2195240
k_means_protein$cluster
##        Albania        Austria        Belgium       Bulgaria Czechoslovakia 
##              2              3              3              2              3 
##        Denmark      E Germany        Finland         France         Greece 
##              3              3              3              3              1 
##        Hungary        Ireland          Italy    Netherlands         Norway 
##              2              3              1              3              3 
##         Poland       Portugal        Romania          Spain         Sweden 
##              3              1              2              1              3 
##    Switzerland             UK           USSR      W Germany     Yugoslavia 
##              3              3              2              3              2
k_means_protein$size
## [1]  4  6 15

Другой способ структурировать вывод kmeans — использовать команду tidy из пакета broom.

tidy(k_means_protein)
##           x1         x2         x3         x4         x5         x6
## 1 -0.5088020 -1.1088009 -0.4124850 -0.8320414  0.9819154  0.1300253
## 2 -0.7901419 -0.5267887 -1.1655757 -0.9047559 -0.9504683  1.4383272
## 3  0.4517373  0.5063957  0.5762263  0.5837801  0.1183432 -0.6100043
##           x7         x8         x9 size withinss cluster
## 1 -0.1842010  1.3108846  1.6292449    4 18.92587       1
## 2 -0.7604664  0.8870168 -0.5373533    6 24.09113       2
## 3  0.3533068 -0.7043759 -0.2195240   15 62.96933       3

Первые девять неназванных переменных — центры кластеров по каждой переменной.

  • Упражнение 10.

Кластеризуйте отмасштабированные данные по арестам в США. В выборе числа кластеров доверьтесь интуиции.

k_means_usa <- kmeans(usa_stand, centers = 3)

Осталось только визуализировать результаты! Для этого будем использовать команду fviz_cluster() из пакета factoextra. Её аргументы — результат кластеризации k_means_protein, исходные данные и ещё куча настроек вроде размера точек и цвета наблюдений-выбросов. Мы только попросим выделять цветом кластеры по их границам и укажем аргумент ellipse.type = 'convex'.

fviz_cluster(object = k_means_protein, data = protein_no_country,
             ellipse.type = 'convex')

  • Упражнение 11.

Визуализируйте результаты кластеризации данных по арестам k_means_usa.

fviz_cluster(object = k_means_usa, data = usa_stand, ellipse.type = 'convex')

  • Пересекаются ли кластеры?
  • Разумно ли изменить количество кластеров?

Как понять, сколько кластеров брать оптимально? Один из способов сделать это — воспользоваться командой fviz_nbclust из пакета factoextra.

g1 <- fviz_nbclust(protein_no_country, kmeans, method = 'wss') +
  labs(subtitle = 'Elbow method')
g1

g2 <- fviz_nbclust(protein_no_country, kmeans, method = 'silhouette') +
  labs(subtitle = 'Silhouette method')
g2

g3 <- fviz_nbclust(protein_no_country, kmeans, method = 'gap_stat') +
  labs(subtitle = 'Gap statistic method')
g3

С помощью хитрого пакета Томаса располагать графики легко! Попробуйте!

(g1 + g2) / g3

g1 + g2 + g3

g1 + (g2 / g3)

  • В чём разница?

  • Как расположит графики команда g1 / (g2 + g3)?

  • Упражнение 12.

Проверьте, совпадает ли выбранное вами число кластеров для usa_stand с оптимальным. Изобразите все три диаграммы вместе: первые две вместе наверху, вторую отеднльно внизу.

p1 <- fviz_nbclust(usa_stand, kmeans, method = 'wss') +
  labs(subtitle = 'Elbow method')
p2 <- fviz_nbclust(usa_stand, kmeans, method = 'silhouette') +
  labs(subtitle = 'Silhouette method')
p3 <- fviz_nbclust(usa_stand, kmeans, method = 'gap_stat') +
  labs(subtitle = 'Gap statistic method')
(p1 + p2) / p3

Метки кластерам легко добавить к исходным данным:

protein_plus <- mutate(protein, cluster = k_means_protein$cluster)
glimpse(protein_plus)
## Observations: 25
## Variables: 11
## $ Country   <chr> "Albania", "Austria", "Belgium", "Bulgaria", "Czecho...
## $ RedMeat   <dbl> 10.1, 8.9, 13.5, 7.8, 9.7, 10.6, 8.4, 9.5, 18.0, 10....
## $ WhiteMeat <dbl> 1.4, 14.0, 9.3, 6.0, 11.4, 10.8, 11.6, 4.9, 9.9, 3.0...
## $ Eggs      <dbl> 0.5, 4.3, 4.1, 1.6, 2.8, 3.7, 3.7, 2.7, 3.3, 2.8, 2....
## $ Milk      <dbl> 8.9, 19.9, 17.5, 8.3, 12.5, 25.0, 11.1, 33.7, 19.5, ...
## $ Fish      <dbl> 0.2, 2.1, 4.5, 1.2, 2.0, 9.9, 5.4, 5.8, 5.7, 5.9, 0....
## $ Cereals   <dbl> 42.3, 28.0, 26.6, 56.7, 34.3, 21.9, 24.6, 26.3, 28.1...
## $ Starch    <dbl> 0.6, 3.6, 5.7, 1.1, 5.0, 4.8, 6.5, 5.1, 4.8, 2.2, 4....
## $ Nuts      <dbl> 5.5, 1.3, 2.1, 3.7, 1.1, 0.7, 0.8, 1.0, 2.4, 7.8, 5....
## $ `Fr&Veg`  <dbl> 1.7, 4.3, 4.0, 4.2, 4.0, 2.4, 3.6, 1.4, 6.5, 6.5, 4....
## $ cluster   <int> 2, 3, 3, 2, 3, 3, 3, 3, 3, 1, 2, 3, 1, 3, 3, 3, 1, 2...

6 Иерархическая кластеризация

Другой способ разбить данные на группы — иерархическая кластеризация. Но, в отличие от метода k-средних, она работает с матрицей расстояний, поэтому первым делом посчитаем её! Для этого будем использовать функцию dist(). Передадим ей стандартизированные данные и укажем явно, как считать расстояния с помощью аргумента method. О всех остальных опциях можно узнать в справке.

protein_dist <- dist(protein_no_country, method = 'euclidian')

Расстояния тоже можно визуализировать! Сделаем это командой fviz_dist из пакета factoextra.

fviz_dist(protein_dist)

Посчитайте матрицу расстояний для таблицы usa_stand и визуализируйте её.

usa_dist <- dist(usa, method = 'euclidian')
fviz_dist(usa_dist)

Полученную матрицу расстояний можно передадать функции hclust(), которая кластеризует данные. Однако в пакете factoextra есть функция hcut(), которая работает с исходными данными. Будем использовать её и попросим выделить четыре кластера в аргументе k.

protein_hcl <- hcut(protein_no_country, k = 4)

С помощью функции fviz_dend визуализируем результат кластеризации. Укажем несколько аргументов, чтобы сделать дендрограмму красивее, а полный перечень найдётся в справке.

fviz_dend(protein_hcl,
          cex = 0.5, # размер подписи
          color_labels_by_k = TRUE) # цвет подписей по группам

Выявленные кластеры можно добавить к исходным данным!

protein_plus2 <- mutate(protein, cluster = protein_hcl$cluster)
glimpse(protein_plus2)
## Observations: 25
## Variables: 11
## $ Country   <chr> "Albania", "Austria", "Belgium", "Bulgaria", "Czecho...
## $ RedMeat   <dbl> 10.1, 8.9, 13.5, 7.8, 9.7, 10.6, 8.4, 9.5, 18.0, 10....
## $ WhiteMeat <dbl> 1.4, 14.0, 9.3, 6.0, 11.4, 10.8, 11.6, 4.9, 9.9, 3.0...
## $ Eggs      <dbl> 0.5, 4.3, 4.1, 1.6, 2.8, 3.7, 3.7, 2.7, 3.3, 2.8, 2....
## $ Milk      <dbl> 8.9, 19.9, 17.5, 8.3, 12.5, 25.0, 11.1, 33.7, 19.5, ...
## $ Fish      <dbl> 0.2, 2.1, 4.5, 1.2, 2.0, 9.9, 5.4, 5.8, 5.7, 5.9, 0....
## $ Cereals   <dbl> 42.3, 28.0, 26.6, 56.7, 34.3, 21.9, 24.6, 26.3, 28.1...
## $ Starch    <dbl> 0.6, 3.6, 5.7, 1.1, 5.0, 4.8, 6.5, 5.1, 4.8, 2.2, 4....
## $ Nuts      <dbl> 5.5, 1.3, 2.1, 3.7, 1.1, 0.7, 0.8, 1.0, 2.4, 7.8, 5....
## $ `Fr&Veg`  <dbl> 1.7, 4.3, 4.0, 4.2, 4.0, 2.4, 3.6, 1.4, 6.5, 6.5, 4....
## $ cluster   <int> 1, 2, 2, 1, 3, 2, 3, 2, 2, 4, 3, 2, 4, 2, 2, 3, 4, 1...
usa_hcl <- hcut(usa, k = 4)
fviz_dend(usa_hcl,
          cex = 0.5, # размер подписи
          color_labels_by_k = TRUE) # цвет подписей по группам

Иерархичская кластеризация полезна и для визуализаций корреляций. Если в функции corrplot() из одноимённого пакета указать аргумента order = hclust, то мы получим сгруппированные по кластерам переменные. Для красоты добавим ещё один аргумент — addrect = 3. Он обведёт прямоугольниками указанное число кластеров.

protein_cor <- cor(protein_no_country)
corrplot(protein_cor, order = 'hclust', addrect = 3)

Визуализируйте корреляции в данных об арестах usa и сгруппируйте их по двум кластерам. Замените кружочки на квадраты, передав аргументу method значение shade.

usa_cor <- cor(usa)
corrplot(usa_cor, order = 'hclust', addrect = 2, method = 'shade')

Добавьте к исходным данным usa кластеры, полученные с помощью иерархической кластеризации:

usa_plus2 <- mutate(usa, cluster = TRUE)
glimpse(usa_plus2)
## Observations: 50
## Variables: 5
## $ Murder   <dbl> 13.2, 10.0, 8.1, 8.8, 9.0, 7.9, 3.3, 5.9, 15.4, 17.4,...
## $ Assault  <int> 236, 263, 294, 190, 276, 204, 110, 238, 335, 211, 46,...
## $ UrbanPop <int> 58, 48, 80, 50, 91, 78, 77, 72, 80, 60, 83, 54, 83, 6...
## $ Rape     <dbl> 21.2, 44.5, 31.0, 19.5, 40.6, 38.7, 11.1, 15.8, 31.9,...
## $ cluster  <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,...

Визуализации [кластеров в известных наборах данных] (https://cran.r-project.org/web/packages/dendextend/vignettes/Cluster_Analysis.html)

Ура! :)