В третьей домашке будем разбираться с млекопитающими!
Данные возьмём из пакета VIM
. Набор называется sleep
, и он содержит информацию о сне млекопитающих.
Описание данных:
BodyWgt — вес, кг
BrainWgt — вес мозга, г
NonD — длительность сна без сноведений, часов в день
Dream — длительность сна со сновидениями, часов в день
Sleep — общая длительность сна, часов в день
Span — максимальная продолжительность жизни в годах
Gest — длительность беременности в днях
Pred — степень хищничества, от 1 до 5, где 5 соответсвует наибольшей вероятности быть убитым другим животным
Exp — степень склонности ко сну, от 1 до 5, где 5 соответствует наибольшей склонности
Danger — общий индекс опасности, от 1 до 5, где 5 соответствует наибольшей опасности
Упражнение 0.
Загрузите необходимые пакеты и скройте все сообщения и предупреждения. За “все” отвечает опция global_options
, а остальное нужно вспомнить :)
library(tidyverse) # обработка данных, графики...
library(skimr)# описательные статистики
library(rio) # импорт фантастического количества форматов данных
library(cluster) # кластерный анализ
library(factoextra) # визуализации kmeans, pca,
library(dendextend) # визуализация дендрограмм
library(corrplot) # визуализация корреляций
library(naniar) # визуализация пропущенных значений
library(visdat) # визуализация пропущенных значений
library(patchwork) # удобное расположение графиков рядом
library(VIM) # данные о сне млекопитающих
Посмотрите на описательные статистики данных sleep
.
Сколько наблюдений в данных?
Сколько переменных?
Есть ли пропущенные значения?
skim(sleep)
## Skim summary statistics
## n obs: 62
## n variables: 10
## Warning: package 'bindrcpp' was built under R version 3.4.4
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100 hist
## Danger 0 62 62 2.61 1.44 1 1 2 4 5 ▇▆▁▅▁▅▁▃
## Exp 0 62 62 2.42 1.6 1 1 2 4 5 ▇▃▁▁▁▂▁▃
## Pred 0 62 62 2.87 1.48 1 2 3 4 5 ▇▇▁▆▁▃▁▇
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75
## BodyWgt 0 62 62 198.79 899.16 0.005 0.6 3.34 48.2
## BrainWgt 0 62 62 283.13 930.28 0.14 4.25 17.25 166
## Dream 12 50 62 1.97 1.44 0 0.9 1.8 2.55
## Gest 4 58 62 142.35 146.81 12 35.75 79 207.5
## NonD 14 48 62 8.67 3.67 2.1 6.25 8.35 11
## Sleep 4 58 62 10.53 4.61 2.6 8.05 10.45 13.2
## Span 4 58 62 19.88 18.21 2 6.62 15.1 27.75
## p100 hist
## 6654 ▇▁▁▁▁▁▁▁
## 5712 ▇▁▁▁▁▁▁▁
## 6.6 ▇▇▇▃▃▁▁▁
## 645 ▇▃▂▁▁▁▁▁
## 17.9 ▃▅▇▇▇▃▂▁
## 19.9 ▅▃▃▇▅▂▂▂
## 100 ▇▃▂▂▁▁▁▁
Постройте диаграмму рассеяния для веса млекопитающего BodyWgt
и веса его мозга BrainWgt
. Раскрасьте точки в соответвии со степенью склонности ко сну Exp
, но помните, что на самом деле эта переменная факторная! Также не забудьте подписать оси и придумать клёвое название. Для удобства оси графика будем строить в логарифмах — за это отвечают слои scale_x_log10()
и scale_y_log10()
.
Для построения графика используйте данные sleep2
.
sleep2 <- mutate(sleep, Exp = factor(Exp))
ggplot(data = sleep2) +
geom_point(aes(x = BodyWgt, y = BrainWgt, color = Exp)) +
scale_x_log10() + # для удобства отмасштабируем оси
scale_y_log10() +
labs(x = 'Логарифм массы тела', y = 'Логарифм массы мозга', title = 'Клёвое название!')
Нарисуйте гистограмму логарифма веса мозга млекопитающих для каждой категории склонности ко сну. Как всегда, подпишите оси и добавьте название.
Для построения графика используйте данные sleep2
.
ggplot(data = sleep2) +
geom_histogram(aes(x = log(BrainWgt))) +
facet_grid(Exp ~ .) +
labs(x = 'Масса мозга', title = 'Вес мозга млекопитащих')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Визуализируйте пропуски в данных sleep
:
С помощью функции vis_miss
, и сгруппируйте их по кластерам;
С помощью функции gg_miss_var
, и количество пропусков считайте в процентах.
vis_miss(sleep, cluster = TRUE)
gg_miss_var(sleep, show_pct = TRUE)
Удалите все строки, содержащие пропуски, и сохраните результат в переменную sleep_remove
.
sleep_remove <- drop_na(sleep)
skim(sleep_remove)
## Skim summary statistics
## n obs: 42
## n variables: 10
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100 hist
## Danger 0 42 42 2.69 1.39 1 1.25 2.5 4 5 ▇▇▁▅▁▆▁▃
## Exp 0 42 42 2.36 1.53 1 1 2 3.75 5 ▇▃▁▂▁▂▁▃
## Pred 0 42 42 2.95 1.45 1 2 3 4 5 ▆▇▁▅▁▅▁▆
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75
## BodyWgt 0 42 42 100.81 402.08 0.005 0.32 2.25 10.41
## BrainWgt 0 42 42 218.68 732.85 0.14 3.6 12.2 155.5
## Dream 0 42 42 1.9 1.39 0 0.9 1.65 2.38
## Gest 0 42 42 129.94 127.84 12 32 90 195
## NonD 0 42 42 8.74 3.84 2.1 6.15 8.5 11
## Sleep 0 42 42 10.64 4.71 2.9 8.05 9.8 13.6
## Span 0 42 42 19.37 20.26 2 5.25 11.2 27
## p100 hist
## 2547 ▇▁▁▁▁▁▁▁
## 4603 ▇▁▁▁▁▁▁▁
## 6.6 ▇▇▇▂▃▁▁▁
## 624 ▇▃▂▁▁▁▁▁
## 17.9 ▅▅▇▆▇▂▂▁
## 19.9 ▅▅▇▇▅▃▃▃
## 100 ▇▂▂▂▁▁▁▁
Заполните пропуски в данных с помощью функции replace_na
следующим образом:
длительность сна без сновидений NonD
средним значением;
длительность сна со сновидениями Dream
медианным значением;
продолжительность жизни Span
— минимальным;
общую длительность сна Sleep
— числом 11;
длительность беременности — числом 90;
и проверьте, что всё получилось :)
sleep_replace <- replace_na(data = sleep,
replace = list(NonD = mean(sleep$NonD, na.rm = TRUE),
Dream = median(sleep$Dream, na.rm = TRUE),
Span = min(sleep$Span, na.rm = TRUE),
Sleep = 11,
Gest = 90))
skim(sleep_replace)
## Skim summary statistics
## n obs: 62
## n variables: 10
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100 hist
## Danger 0 62 62 2.61 1.44 1 1 2 4 5 ▇▆▁▅▁▅▁▃
## Exp 0 62 62 2.42 1.6 1 1 2 4 5 ▇▃▁▁▁▂▁▃
## Pred 0 62 62 2.87 1.48 1 2 3 4 5 ▇▇▁▆▁▃▁▇
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75
## BodyWgt 0 62 62 198.79 899.16 0.005 0.6 3.34 48.2
## BrainWgt 0 62 62 283.13 930.28 0.14 4.25 17.25 166
## Dream 0 62 62 1.94 1.29 0 1.05 1.8 2.27
## Gest 0 62 62 138.98 142.5 12 39 90 195
## NonD 0 62 62 8.67 3.22 2.1 6.8 8.67 10.55
## Sleep 0 62 62 10.56 4.45 2.6 8.22 10.7 13.15
## Span 0 62 62 18.72 18.15 2 5.25 13.35 27
## p100 hist
## 6654 ▇▁▁▁▁▁▁▁
## 5712 ▇▁▁▁▁▁▁▁
## 6.6 ▃▅▇▂▂▁▁▁
## 645 ▇▂▂▁▁▁▁▁
## 17.9 ▂▂▃▇▃▁▁▁
## 19.9 ▃▂▃▇▃▂▂▂
## 100 ▇▃▂▂▁▁▁▁
Приступим к кластерному анализу!
Первым шагом отмасштабируйте все столбцы sleep_replace
с числовыми значениями так, чтобы получившиеся данные остались в формате data.frame
. Сохраните их в переменную sleep_stand
. Затем сделайте кластерный анализ методом k-средних, выбрав произвольное количество кластеров. И наконец, визуализируйте результат!
sleep_stand <- mutate_if(sleep_replace, is.numeric, ~ as.vector(scale(.)))
skim(sleep_stand)
## Skim summary statistics
## n obs: 62
## n variables: 10
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75 p100
## BodyWgt 0 62 62 -2.2e-18 1 -0.22 -0.22 -0.22 -0.17 7.18
## BrainWgt 0 62 62 6e-18 1 -0.3 -0.3 -0.29 -0.13 5.84
## Danger 0 62 62 7.9e-17 1 -1.12 -1.12 -0.43 0.96 1.66
## Dream 0 62 62 -1.2e-18 1 -1.5 -0.69 -0.11 0.26 3.6
## Exp 0 62 62 -3.2e-17 1 -0.88 -0.88 -0.26 0.98 1.61
## Gest 0 62 62 3e-17 1 -0.89 -0.7 -0.34 0.39 3.55
## NonD 0 62 62 -1.9e-16 1 -2.04 -0.58 0 0.58 2.87
## Pred 0 62 62 -6.7e-17 1 -1.27 -0.59 0.087 0.76 1.44
## Sleep 0 62 62 7.8e-17 1 -1.79 -0.52 0.031 0.58 2.1
## Span 0 62 62 7.8e-17 1 -0.92 -0.74 -0.3 0.46 4.48
## hist
## ▇▁▁▁▁▁▁▁
## ▇▁▁▁▁▁▁▁
## ▇▆▁▅▁▅▁▃
## ▃▅▇▂▂▁▁▁
## ▇▃▁▁▁▂▁▃
## ▇▂▂▁▁▁▁▁
## ▂▂▃▇▃▂▁▁
## ▇▇▁▆▁▁▃▇
## ▃▂▃▇▃▂▂▂
## ▇▃▂▂▁▁▁▁
sleep_kmeans <- kmeans(sleep_stand, centers = 2)
fviz_cluster(object = sleep_kmeans, data = sleep_stand, ellipse.type = 'convex')
Найдите оптимальное число кластеров. Если оно не совпало с выбранным вами, перезапустите предыдущий кусок кода, изменив нужный параметр функции kmeans
.
fviz_nbclust(sleep_stand, kmeans, method = 'gap_stat') +
labs(subtitle = 'Gap statistic method')
Добавьте центры кластеров к данным sleep_replace
.
sleep_plus <- mutate(sleep_replace, cluster = sleep_kmeans$cluster)
glimpse(sleep_plus)
## Observations: 62
## Variables: 11
## $ BodyWgt <dbl> 6654.000, 1.000, 3.385, 0.920, 2547.000, 10.550, 0.02...
## $ BrainWgt <dbl> 5712.0, 6.6, 44.5, 5.7, 4603.0, 179.5, 0.3, 169.0, 25...
## $ NonD <dbl> 8.672917, 6.300000, 8.672917, 8.672917, 2.100000, 9.1...
## $ Dream <dbl> 1.8, 2.0, 1.8, 1.8, 1.8, 0.7, 3.9, 1.0, 3.6, 1.4, 1.5...
## $ Sleep <dbl> 3.3, 8.3, 12.5, 16.5, 3.9, 9.8, 19.7, 6.2, 14.5, 9.7,...
## $ Span <dbl> 38.6, 4.5, 14.0, 2.0, 69.0, 27.0, 19.0, 30.4, 28.0, 5...
## $ Gest <dbl> 645, 42, 60, 25, 624, 180, 35, 392, 63, 230, 112, 281...
## $ Pred <int> 3, 3, 1, 5, 3, 4, 1, 4, 1, 1, 5, 5, 2, 5, 1, 2, 2, 2,...
## $ Exp <int> 5, 1, 1, 2, 5, 4, 1, 5, 2, 1, 4, 5, 1, 5, 1, 2, 2, 2,...
## $ Danger <int> 3, 3, 1, 3, 4, 4, 1, 4, 1, 1, 4, 5, 2, 5, 1, 2, 2, 2,...
## $ cluster <int> 2, 1, 1, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 1, 1, 1,...
Сделайте иерархическую кластеризацию c пятью кластерами и визуализируйте результат.
sleep_hcl <- hcut(sleep_stand, k = 5)
fviz_dend(sleep_hcl,
cex = 0.5, # размер подписи
color_labels_by_k = TRUE) # цвет подписей по группам
Посчитайте корреляции по матрице без пропусков sleep_replace
и визулизируйте их. Выделите квадратами две близкие по корреляциям группы.
sleep_cor <- cor(sleep_replace)
corrplot(sleep_cor, order = 'hclust', addrect = 2)
Ура, домашка сделана! :)