Семинар 3. Факторный анализ

Шаманское заклинание для настройки глобальных опций отчёта:

library(tidyverse) # обработка данных, графики...
library(skimr)# описательные статистики
library(rio) # импорт фантастического количества форматов данных

library(cluster) # кластерный анализ
library(factoextra) # визуализации kmeans, pca,
library(dendextend) # визуализация дендрограмм

library(corrplot) # визуализация корреляций

library(broom) # метла превращает результаты оценивания моделей в таблички

library(naniar) # визуализация пропущенных значений
library(visdat) # визуализация пропущенных значений

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

1 Немного формул для k-средних! :)

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

protein <- import('data/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
##  ▅▇▁▂▂▁▁▂
##  ▁▂▁▁▇▂▃▂
##  ▅▇▂▃▁▂▁▁
##  ▅▇▇▇▂▁▆▃
##  ▂▇▂▇▃▇▁▁
##  ▇▇▁▃▂▆▁▁
##  ▂▆▇▆▁▃▁▂
##  ▃▁▅▃▆▇▆▇
##  ▁▃▇▁▁▇▃▃

Отмасштабируем все числовые переменные с помощью функции scale(). Затем спрячем текстовую переменную Country в названия строк:

protein_no_country <- protein %>%
  mutate_if(is.numeric, ~ as.vector(scale(.))) %>%
  column_to_rownames(var = 'Country')

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

Выполним кластеризацию методом k-средних с помощью функции kmeans.

В качестве аргументов укажем отмасштабированные данные protein_no_country и количество кластеров centers. Берём три кластера! Сохраним результат кластеризации в список k_means_protein.

set.seed(13)
k_means_protein <- kmeans(protein_no_country, centers = 3)
k_means_protein

## K-means clustering with 3 clusters of sizes 10, 4, 11
## 
## Cluster means:
##        RedMeat  WhiteMeat       Eggs       Milk       Fish    Cereals
## 1 -0.677605893 -0.7595936 -0.8643394 -0.8756701 -0.1775148  0.9150064
## 2  0.006572897 -0.2290150  0.1914789  1.3458748  1.1582546 -0.8722721
## 3  0.613615213  0.7738178  0.7161344  0.3066547 -0.2598064 -0.5146342
##       Starch       Nuts     Fr&Veg
## 1 -0.5299602  1.0565639  0.3292860
## 2  0.1676780 -0.9553392 -1.1148048
## 3  0.4208082 -0.6131165  0.1060327
## 
## Clustering vector:
##        Albania        Austria        Belgium       Bulgaria Czechoslovakia 
##              1              3              3              1              3 
##        Denmark      E Germany        Finland         France         Greece 
##              2              3              2              3              1 
##        Hungary        Ireland          Italy    Netherlands         Norway 
##              1              3              1              3              2 
##         Poland       Portugal        Romania          Spain         Sweden 
##              3              1              1              1              2 
##    Switzerland             UK           USSR      W Germany     Yugoslavia 
##              3              3              1              3              1 
## 
## Within cluster sum of squares by cluster:
## [1] 70.957748  5.900318 37.885439
##  (between_SS / total_SS =  46.9 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

Внутригрупповая сумма квадратов расстояний, WSS (within): \[ WSS = ||x_1 - c_1||^2 + ||x_2 - c_2 ||^2 + \ldots + ||x_n - c_n||^2, \] где \(x_i\) — координаты \(i\)-го наблюдения, \(c_i\) — координаты центра кластера, к которому относится \(i\)-ое наблюдение. У наблюдений из одного кластера величины \(c_i\) равны.

Узнаем величину WSS с разбивкой по кластерам:

k_means_protein$withinss

## [1] 70.957748  5.900318 37.885439

Межгрупповая сумма квадратов расстояний, BSS (between): \[ BSS = ||c_1 - \bar x||^2 + ||c_2 - \bar x||^2 + \ldots + ||c_n - \bar x||^2 \]

k_means_protein$betweenss

## [1] 101.2565

Общая сумма квадратов расстояний, TSS (total): \[ TSS = ||x_1 - \bar x||^2 + ||x_2 - \bar x||^2 + \ldots + ||x_n - \bar x||^2 \]

k_means_protein$totss

## [1] 216

Суть алгоритма \(k\)-средних в терминах квадратов расстояний проста:

Кластеры определяются так, чтобы величина WSS была минимально возможной.

По теореме Пифагора: \[ TSS = WSS + BSS \]

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

g1 <- fviz_nbclust(protein_no_country, kmeans, method = 'wss') +
  labs(title = 'Зависимость WSS от числа кластеров',
       subtitle = 'Метод: на глаз!', x = 'Число кластеров',
       y = 'Внутригрупповая сумма квадратов расстояний')
g1

g2 <- fviz_nbclust(protein_no_country, kmeans, method = 'silhouette') +
  labs(subtitle = 'Метод силуэтов',
       title = 'Зависимость средней ширины силуэта от числа кластеров',
       y = 'Средняя ширина силуэта по всем точкам',
       x = 'Число кластеров')
g2

Обозначения:

• \(a_i\) — среднее расстояние от точки \(i\) до точек своего кластера;

• \(b_i\) — среднее расстояние от точки \(i\) до точек ближайшего кластера-конкурента;

• Ширина силуэта точки \(i\):

\[ s_i = \frac{b_i - a_i}{\max\{a_i, b_i\}} \]

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

Обозначения:

• \(k\) — число кластеров
• \(WSS_{random}\) — средняя величина внутригрупповой суммы квадратов расстояний, если генерировать точки равномерно в прямоугольнике имеющихся данных.

Статистика разрыва (gap-статистика): \[ Gap(k) = \ln WSS_{random}(k) - \ln WSS(k) \]

На графике показаны значения статистики разрыва плюс доверительный интервал шириной в одной стандартное отклонение.

Автоматически компьютер подбирает такое количество количество кластеров \(k\), начиная с которого статистика разрыва растёт слишком медленно или вовсе падает. «Слишком медленно» растёт означает, что следующее значение статистики попадает в доверительный интервал предыдущего.

Более формально, автоматически находится первое \(k\) при котором: \[ Gap(k) + 1 \cdot se(Gap(k)) \geq Gap(k+1) \]

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

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

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

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

fviz_dist(protein_dist)

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

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

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

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

protein_hcl <- hcut(protein_no_country, k = 4,
                    hc_metric = 'euclidean', hc_method = 'ward.D2')

Возможны бесчиленные вариации алгоритма!

По методу подсчёта расстояния между двумя точками, hc_metric =

• ‘euclidean’ — привычное Евклидово расстояние;
• ‘manhattan’ — расстояние по Майкопски! :)

Майкоп

По методу подсчёта расстояний между кластерами, hc_method =

• ‘single’ — расстояние между ближайшими точками этих кластеров;
• ‘complete’ — расстояние между самыми далёкими точками этих кластеров;
• ‘average’ — среднее расстояние между точками кластеров;
• ‘ward.D2’ — на сколько увеличится WSS, если слить два кластера в один;

С помощью функции 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...

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

• Сделайте иерархическую кластеризацию с четыремя группами на данных об арестах.

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

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)

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

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

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

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

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

usa_plus2 <- mutate(usa, cluster = usa_hcl$cluster)
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  <int> 1, 1, 1, 2, 1, 2, 3, 1, 1, 2, 4, 3, 1, 3, 4, 3, 3, 1,...

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

3 Метод главных компонент

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

Например, можно заменить 20 исходных переменных на две искусственные главные компоненты, чтобы изобразить многомерный набор данных на двумерном графике.

Как переменные превращают в главные компоненты?

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

Подход А. Максимизация разброса.

В первую компоненту берём исходные три переменные с такими весами, чтобы:

• максимизировать разброс первой главной компоненты;
• сумма квадратов весов равнялась бы единице.

Формально: \[ pc^{1}_i = \alpha_1 x_i + \alpha_2 y_i + \alpha_3 z_i \]

МГК максимизирует \(TSS(pc^1)=\sum (pc^1_i - 0)^2\) при ограничении \(\alpha_1^2 + \alpha_2^2 + \alpha_3^2=1\).

Подход Б. Минимизация расстояний от точек до новой системы координат.

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

\[ \sum (x_i - \hat x_i)^2 + \sum (y_i - \hat y_i)^2 + \sum (z_i - \hat z_i)^2 \to \min, \] при условии, что точки \(\hat x_i\), \(\hat y_i\), \(\hat z_i\) лежат на одной прямой.

Эти подходы эквивалентны.

После применения метода главных компонент оказывается, что:

\[ TSS(pc^1) + TSS(pc^2) + TSS(pc^3) = TSS(x) + TSS(y) + TSS(z), \] где \(TSS\) — сумма квадратов значений переменной. Напомним, что переменные у нас центрированы.

Реализуем метод главных компонент на данных о потреблении белка функцией prcomp(). Передадим ей отмасштабированные данные protein_no_country, хотя так поступать и необязательно. Другой путь — передать исходный набора данных и попросить функцию procmp() стандартизировать их аргументом scale = TRUE.

protein_pca <- prcomp(protein_no_country)
protein_pca

## Standard deviations (1, .., p=9):
## [1] 2.0016087 1.2786710 1.0620355 0.9770691 0.6810568 0.5702026 0.5211586
## [8] 0.3410160 0.3148204
## 
## Rotation (n x k) = (9 x 9):
##                  PC1         PC2         PC3          PC4         PC5
## RedMeat   -0.3026094 -0.05625165 -0.29757957 -0.646476536  0.32216008
## WhiteMeat -0.3105562 -0.23685334  0.62389724  0.036992271 -0.30016494
## Eggs      -0.4266785 -0.03533576  0.18152828 -0.313163873  0.07911048
## Milk      -0.3777273 -0.18458877 -0.38565773  0.003318279 -0.20041361
## Fish      -0.1356499  0.64681970 -0.32127431  0.215955001 -0.29003065
## Cereals    0.4377434 -0.23348508  0.09591750  0.006204117  0.23816783
## Starch    -0.2972477  0.35282564  0.24297503  0.336684733  0.73597332
## Nuts       0.4203344  0.14331056 -0.05438778 -0.330287545  0.15053689
## Fr&Veg     0.1104199  0.53619004  0.40755612 -0.462055746 -0.23351666
##                   PC6         PC7         PC8        PC9
## RedMeat   -0.45986989  0.15033385 -0.01985770  0.2459995
## WhiteMeat -0.12100707 -0.01966356 -0.02787648  0.5923966
## Eggs       0.36124872 -0.44327151 -0.49120023 -0.3333861
## Milk       0.61843780  0.46209500  0.08142193  0.1780841
## Fish      -0.13679059 -0.10639350 -0.44873197  0.3128262
## Cereals    0.08075842  0.40496408 -0.70299504  0.1522596
## Starch     0.14766670  0.15275311  0.11453956  0.1218582
## Nuts       0.44701001 -0.40726235  0.18379989  0.5182749
## Fr&Veg     0.11854972  0.44997782  0.09196337 -0.2029503

Посмотрим, что лежит в этом списке!

attributes(protein_pca)

## $names
## [1] "sdev"     "rotation" "center"   "scale"    "x"       
## 
## $class
## [1] "prcomp"

Сами новые искусственные главные компоненты лежат в матрице protein_pca$x. Например, первая главная компонента:

protein_pca$x[, 1]

##        Albania        Austria        Belgium       Bulgaria Czechoslovakia 
##      3.4853673     -1.4226694     -1.6220323      3.1340813     -0.3704646 
##        Denmark      E Germany        Finland         France         Greece 
##     -2.3652688     -1.4222108     -1.5638563     -1.4879824      2.2397000 
##        Hungary        Ireland          Italy    Netherlands         Norway 
##      1.4574398     -2.6634775      1.5345653     -1.6414454     -0.9747029 
##         Poland       Portugal        Romania          Spain         Sweden 
##     -0.1218695      1.7058540      2.7568124      1.3118074     -1.6337300 
##    Switzerland             UK           USSR      W Germany     Yugoslavia 
##     -0.9123182     -1.7353682      0.7825965     -2.0938353      3.6230077

Выборочные стандартные отклонения компонент (корни из \(\lambda_j\)) лежат в векторе protein_pca$sdev:

protein_pca$sdev

## [1] 2.0016087 1.2786710 1.0620355 0.9770691 0.6810568 0.5702026 0.5211586
## [8] 0.3410160 0.3148204

Матрица protein_pca$rotation содержит веса, с которыми исходные переменные входят в новые искуственные главные компоненты. Например, в первой главной компоненте лежат исходные переменные с весами:

protein_pca$rotation[, 1]

##    RedMeat  WhiteMeat       Eggs       Milk       Fish    Cereals 
## -0.3026094 -0.3105562 -0.4266785 -0.3777273 -0.1356499  0.4377434 
##     Starch       Nuts     Fr&Veg 
## -0.2972477  0.4203344  0.1104199

Визуализируем данные в осях первых двух главных компонент. Для этого воспользуемся функцией fviz_pca_ind() из пакета factoextra. Рисовать будем protein_pca, а аргумент repel = TRUE укажем для того, чтобы подписи на графике не перекрывали друг друга.

fviz_pca_ind(protein_pca, repel = TRUE)

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

• Выделите главные компоненты на данных об арестах в Америке. Будьте внимательны: эти данные мы не масштабировали!

• Визуализируйте данные в осях первых двух главных компонент. Проследите, чтобы подписи точек на графике были аккуратными :)

usa_pca <- prcomp(usa, scale. = TRUE)
usa_pca

## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432

fviz_pca_ind(usa_pca, repel = TRUE)

• Какой процент дисперсии объясняют в сумме первые две главные компоненты?
• Какой из штатов скорее можно назвать нетипичным, Калифорнию или Вирджинию?
• Какой из штатов наиболее похож на Техас?

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

fviz_pca_biplot(protein_pca, repel = TRUE)

• Какая исходная переменная входит в первую главную компоненту с большим положительным весом?
• Какая исходная переменная входит в первую главную компоненту с большим отрицательным весом?

• Как можно хотя бы примерно проинтепретировать первую главную компоненту?

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

Визуализируйте в осях первых двух главных компонент американские штаты и проекции исходных переменных.

fviz_pca_biplot(usa_pca, repel = TRUE)

• Какая исходная переменная входит в первую главную компоненту с большим отрицательным весом?
• Как можно хотя бы примерно проинтепретировать первую главную компоненту?
• Какие штаты можно считать самыми безопасными?

Теперь изобразим, какой процент разброса данных объясняет каждая главная компонента. Будем использовать команду fviz_eig() из пакета factoextra.

fviz_eig(protein_pca)

• Сколько нужно взять главных компонент, чтобы объяснить более 70% разброса исходных наблюдений?

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

Визиулизируйте процент разброса, который объясняет каждая главная компонента. Обратите внимание, что у функции появился новый аргумент :)

fviz_eig(usa_pca, addlabels = TRUE)

• Какой процент разброса в данных объясняют в сумме третья и четвёртая глаыне компоненты?
• Что меняет аргумент addlabels = TRUE?

Выясним, какой влкад вносит каждая переменная в конкретную главную компоненту. Для этого будем использовать команду fviz_contrib(). Передадим ей объект protein_pca и в качестве аргументов укажем choice = 'var', так как нам нужны именно переменные, а не наблюдения, и axes = 1, чтобы посмотреть на первую главную компоненту. Для визуализации вклада наблюдения значение аргумента choice поменяем на ind и для разнообразия посмотрим на третью главную компоненту :)

fviz_contrib(protein_pca, choice = 'var', axes = 1)

fviz_contrib(protein_pca, choice = 'ind', axes = 3)

• Какие переменные имеют наибольший вес в первой главной компоненте?

• Какая страна вносит наибольший вклад в сумму квадратов расстояний до центра для третьей главной компоненты?

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

Для данных по арестам визуализируйте вклад переменных во вторую главную компоненту и вклад наблюдений в первую.

fviz_contrib(usa_pca, choice = 'var', axes = 2)

fviz_contrib(usa_pca, choice = 'ind', axes = 1)

В осях первых двух главных компонент исходные данные можно раскрасить согласно кластерам:

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

Теперь мы понимаем, что на этом графике!

Ура! :)