В третьей домашке будем разбираться с млекопитающими!

Данные возьмём из пакета VIM. Набор называется sleep, и он содержит информацию о сне млекопитающих.

Описание данных:

Загрузите необходимые пакеты и скройте все сообщения и предупреждения. За “все” отвечает опция 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:

  1. С помощью функции vis_miss, и сгруппируйте их по кластерам;

  2. С помощью функции 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 следующим образом:

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

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)

Ура, домашка сделана! :)