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

Существует два эквивалентных взгляда на метод главных компонент.

  1. Последовательная максимизация выборочного разброса новых переменных:
  1. Последовательная минимизация суммы квадратов расстояний от старых точек до новых точек:

Формально:

Мы хотим вместо \(k\) переменных оставить \(p<k\). Другими словами, хотим заменить матрицу \(X_{n\times k}\) на матрицу \(\hat{X}_{n\times k}\) ранга \(p\), так чтобы минимизировать сумму квадратов ошибок: \[ \min \sum_{j=1}^k \sum_{i=1}^n (x_{ij}-\hat{x}_{ij})^2 \]

Частный случай: от двух переменных к одной

Что делает регрессия \(y\) на \(x\)? Строит прямую так, чтобы сумма квадратов вертикальных отклонений точек от прямой была бы минимальна.

Что делает регрессия \(x\) на \(y\)? Строит прямую так, чтобы сумма квадратов горизонтальных отклонений точек от прямой была бы минимальна.

Что делает выделение главной компоненты пары векторов \(x\) и \(y\)? Строит прямую так, чтобы сумма квадратов отклонений точек от прямой по перпендикуляру была бы минимальна. И от двух координат на плоскости переходим к одной координате на прямой.

Шампур — это первая главная компонента шашлыка!

Метод главных компонент на языке линейной алгебры

про квадратичные формы

про собственные числа и собственные векторы

про SVD

В R всё просто и быстро:

Загружаем нужные пакеты:

library("knitr")
# opts_chunk$set(cache=TRUE) # кэшируем фрагменты кода

library("ggplot2") # для построения графиков
library("HSAUR") # нам нужен только набор данных по многоборью
# install.packages("HSAUR") # может быть нужно установить пакеты...

opts_chunk$set(fig.align = 'center') 
options(width = 110)
theme_set(theme_bw())

Загружаем интересующий нас набор данных по многоборью:

data(heptathlon)
head(heptathlon)
##                     hurdles highjump  shot run200m longjump javelin run800m score
## Joyner-Kersee (USA)   12.69     1.86 15.80   22.56     7.27   45.66  128.51  7291
## John (GDR)            12.85     1.80 16.23   23.65     6.71   42.56  126.12  6897
## Behmer (GDR)          13.20     1.83 14.20   23.10     6.68   44.54  124.20  6858
## Sablovskaite (URS)    13.61     1.80 15.23   23.92     6.25   42.78  132.24  6540
## Choubenkova (URS)     13.51     1.74 14.76   23.93     6.32   47.46  127.90  6540
## Schulz (GDR)          13.75     1.83 13.50   24.65     6.33   42.82  125.79  6411

Корреляции исходных переменных:

print(cor(heptathlon), digits = 2)
##          hurdles highjump  shot run200m longjump javelin run800m score
## hurdles   1.0000  -0.8114 -0.65    0.77   -0.912 -0.0078    0.78 -0.92
## highjump -0.8114   1.0000  0.44   -0.49    0.782  0.0022   -0.59  0.77
## shot     -0.6513   0.4408  1.00   -0.68    0.743  0.2690   -0.42  0.80
## run200m   0.7737  -0.4877 -0.68    1.00   -0.817 -0.3330    0.62 -0.86
## longjump -0.9121   0.7824  0.74   -0.82    1.000  0.0671   -0.70  0.95
## javelin  -0.0078   0.0022  0.27   -0.33    0.067  1.0000    0.02  0.25
## run800m   0.7793  -0.5912 -0.42    0.62   -0.700  0.0200    1.00 -0.77
## score    -0.9232   0.7674  0.80   -0.86    0.950  0.2531   -0.77  1.00

Восьмой столбец, score, нам не нужен, т.к. это не результат соревнования, а итоговый балл по многоборью. Применяем метод главных компонент к стандартизированным переменным. Стандартизация здесь нужна, так как разные переменные измерены в несопоставимых единицах измерения:

h <- heptathlon[, -8]
hept_pca <- prcomp(h, scale = TRUE)

Теперь можно раздобыть сами главные компоненты:

pc <- hept_pca$x
head(pc)
##                           PC1         PC2        PC3         PC4        PC5         PC6        PC7
## Joyner-Kersee (USA) -4.121448 -1.24240435 -0.3699131 -0.02300174  0.4260062 -0.33932922  0.3479213
## John (GDR)          -2.882186 -0.52372600 -0.8974147  0.47545176 -0.7030659  0.23808730  0.1440158
## Behmer (GDR)        -2.649634 -0.67876243  0.4591767  0.67962860  0.1055252 -0.23919071 -0.1296478
## Sablovskaite (URS)  -1.343351 -0.69228324 -0.5952704  0.14067052 -0.4539282  0.09180564 -0.4865780
## Choubenkova (URS)   -1.359026 -1.75316563  0.1507013  0.83595001 -0.6871948  0.12630397  0.2394820
## Schulz (GDR)        -1.043847  0.07940725  0.6745305  0.20557253 -0.7379335 -0.35578939 -0.1034143

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

loadings <- hept_pca$rotation
head(loadings)
##                 PC1         PC2         PC3         PC4         PC5         PC6        PC7
## hurdles   0.4528710 -0.15792058 -0.04514996  0.02653873 -0.09494792 -0.78334101 -0.3802471
## highjump -0.3771992  0.24807386  0.36777902 -0.67999172 -0.01879888 -0.09939981 -0.4339311
## shot     -0.3630725 -0.28940743 -0.67618919 -0.12431725 -0.51165201  0.05085983 -0.2176249
## run200m   0.4078950  0.26038545  0.08359211 -0.36106580 -0.64983404  0.02495639  0.4533848
## longjump -0.4562318  0.05587394 -0.13931653 -0.11129249  0.18429810 -0.59020972  0.6120639
## javelin  -0.0754090 -0.84169212  0.47156016 -0.12079924 -0.13510669  0.02724076  0.1729467

Посмотрим как первая компонента связана с итоговым результатом спортсмена:

qplot(x = heptathlon$score, y = pc[, 1]) + 
  labs(x = "Сумма баллов по многоборью", 
       y = "Первая главная компонента",
       title = "Связь первой компоненты с итоговым результатом")

cor(heptathlon$score, pc[, 1])
## [1] -0.9910978

Мораль такая: не зная способ подсчета очков за семиборье мы его восстановили! Первая главная компонента вбирает в себя максимум отличий между спортсменами :)

Смотрим какую долю дисперсии объясняют главные компоненты:

summary(hept_pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.1119 1.0928 0.72181 0.67614 0.49524 0.27010 0.2214
## Proportion of Variance 0.6372 0.1706 0.07443 0.06531 0.03504 0.01042 0.0070
## Cumulative Proportion  0.6372 0.8078 0.88223 0.94754 0.98258 0.99300 1.0000

Для построения графика сначала получим вектор долей дисперсии, объясняемых переменными:

exp_percent <- summary(hept_pca)$importance[2, ]
qplot(y = exp_percent, x = names(exp_percent)) + 
  geom_bar(stat = 'identity') +
  labs(x = "Главные компоненты",
       y = "Доли диспесии",
       title = "Доли дисперсии, объясняемые главными компонентами \n (без предварительной стандартизации)")

График первых двух компонент с вкладом каждой переменной:

biplot(hept_pca, 
       xlim = c(-0.8, 0.8), ylim = c(-0.6, 0.6))

Красные векторы — это проекции единичных (?) векторов исходных координат на плоскость первых двух главных компонент.

На этом графике видно, например, что:

Вращение

После PCA иногда дополнительно вращают главные компоненты. Это уже выходит не PCA, а PCA плюс вращение! Кажется, SPSS делает по дефолту PCA + вращение, поэтому результаты не совпадают с R.

Здесь про varimax

Почиташки: