Семинар 4. Кластерный анализ
Октября 1, 2018
Признаемся: часть сокровищ 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.4library(rio) # импорт фантастического количества форматов данных
library(cluster) # кластерный анализ## Warning: package 'cluster' was built under R version 3.4.4library(factoextra) # визуализации kmeans, pca,## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZlibrary(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':
##
## cutreelibrary(corrplot) # визуализация корреляций## corrplot 0.84 loadedlibrary(broom) # метла превращает результаты оценивания моделей в таблички## Warning: package 'broom' was built under R version 3.4.4library(naniar) # визуализация пропущенных значений## Warning: package 'naniar' was built under R version 3.4.3##
## Attaching package: 'naniar'## The following object is masked from 'package:skimr':
##
## n_completelibrary(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] 15max(a)## [1] NAmax(a, na.rm = TRUE)## [1] 8mean(a)## [1] NAmean(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 -51.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 6class(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 NAattributes(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') объединяет только те строки обеих таблиц, в которых нет пропущенных наблюдений
- Упражнение 2.
Экспериментируем!
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 -5right_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 -5inner_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 -5full_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Сколько наблюдений получается в объединённой таблице в каждом случае?
Упражнение 3.
Нужно объединить таблицы 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-...- Присоедините к таблице
flightsтаблицуweatherс помощью функцииleft_join(). В качестве аргументаbyиспользуйте вектор из названий общих столбцов:year,month,day,hourиorigin. - Сколько наблюдений в объединённой таблице?
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_join(). Сколько наблюдений в получившейся таблице? - В чём разница
left_joinиinner_join?
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.8gdp## 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
Не упомянутые переменные останутся строками новой широкой таблицы.
- Упражнение 4.
Из данных про воздух в Нью-Йорке 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.54 Работа с пропущенными данными :)
Команда 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)
- Упражнение 5.
Визуализируйте пропуски во встроенных данных о воздухе в Нью-Йорке 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)
- Какие признаки содержат пропущенные наблюдения?
- Какова доля пропусков в данных
airquality?
Команда 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 ▆▁▃▇▇▅▁▁Упражнение 7.
- Удалите все наблюдения с пропусками из данных
airquality. Сколько наблюдений осталось? - Заполните пропуски в данных. Для переменной
Ozone— средним значением, а дляSolar.R— медианным. Проверьте, что всё получилось :)
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^* = (x_i - \bar x) / \hat \sigma_x, \] где \(\bar x = \frac{x_1 + x_2 + \ldots + x_n}{n}\), а \[ \hat\sigma_x = \frac{(x_1 - \bar x)^2 + \ldots + (x_n - \bar x)^2}{n-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.2195240k_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 2k_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)
- Упражнение 13.
Посчитайте матрицу расстояний для таблицы 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...Упражнение 14.
Сделайте иерархическую кластеризацию с четыремя группами на данных об арестах.
Визуализируйте результат кластеризации и сделайте подписи цветными.
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)
- Упражнение 15.
Визуализируйте корреляции в данных об арестах usa и сгруппируйте их по двум кластерам. Замените кружочки на квадраты, передав аргументу method значение shade.
usa_cor <- cor(usa)
corrplot(usa_cor, order = 'hclust', addrect = 2, method = 'shade')
- Упражнение 16.
Добавьте к исходным данным 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)
Ура! :)