Сегодня мы поработаем с базой данных, которая не очень новая (1977 года), но достаточно интересная. В файле USJR.csv сохранены данные по оценкам судей Верховного суда США. Оценки выставлялись адвокатами. Оценивались следующие аспекты поведения в ходе судебного процесса:

Оригинал описания данных см. здесь.

Загрузим базу данных:

d <- read.csv("USJR.csv")

Посмотрим на нее:

head(d)
##                X CONT INTG DMNR DILG CFMG DECI PREP FAMI ORAL WRIT PHYS
## 1  AARONSON,L.H.  5.7  7.9  7.7  7.3  7.1  7.4  7.1  7.1  7.1  7.0  8.3
## 2 ALEXANDER,J.M.  6.8  8.9  8.8  8.5  7.8  8.1  8.0  8.0  7.8  7.9  8.5
## 3 ARMENTANO,A.J.  7.2  8.1  7.8  7.8  7.5  7.6  7.5  7.5  7.3  7.4  7.9
## 4    BERDON,R.I.  6.8  8.8  8.5  8.8  8.3  8.5  8.7  8.7  8.4  8.5  8.8
## 5   BRACKEN,J.J.  7.3  6.4  4.3  6.5  6.0  6.2  5.7  5.7  5.1  5.3  5.5
## 6     BURNS,E.B.  6.2  8.8  8.7  8.5  7.9  8.0  8.1  8.0  8.0  8.0  8.6
##   RTEN
## 1  7.8
## 2  8.7
## 3  7.8
## 4  8.7
## 5  4.8
## 6  8.6

Используя библиотекуdplyr, выберем из базы d все столбцы, кроме X, X - это имя судьи, текстовый показатель, который нам не понадобится в реализации метода главных компонент. Сохраним нужные переменные в базу m.

library(dplyr)
m <- d %>% select(-X)

Пояснения к коду: возьмем базу d, подадим ее на вход функции select(). Аselect(-X) - выбрать все столбцы, кроме X (поэтому перед X стоит минус).

Установим (если еще не устанавливали) библиотеку GGally для графиков:

install.packages("GGally")

Обратимся к библиотеке и построим матрицу диаграмм рассеяния:

library(GGally)
ggpairs(m)

Теперь реализуем метод главных компонент и посмотрим на результаты.

pca <- prcomp(m, center = TRUE, scale = TRUE)

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

pca
## Standard deviations:
##  [1] 3.18331647 1.05078398 0.57697626 0.50383231 0.29060762 0.19309598
##  [7] 0.14029545 0.12415832 0.08850690 0.07491146 0.05708042 0.04539134
## 
## Rotation:
##               PC1          PC2          PC3          PC4          PC5
## CONT  0.003075143 -0.932890644  0.334756548  0.058576867  0.093438368
## INTG -0.288550775  0.182040993  0.549360126  0.173977074 -0.014543880
## DMNR -0.286884206  0.197565743  0.556490386 -0.124412022 -0.228832817
## DILG -0.304354091 -0.036304667 -0.163629910  0.321395544 -0.301936920
## CFMG -0.302572733 -0.168393523 -0.207341904  0.012949223 -0.448430522
## DECI -0.301891969 -0.127877299 -0.297902771  0.030491779 -0.424003128
## PREP -0.309406446 -0.032230248 -0.151869345  0.213656069  0.202853400
## FAMI -0.306679527  0.001315183 -0.195290454  0.200651140  0.507470003
## ORAL -0.312708348  0.003625720 -0.002150634 -0.007441042  0.246059421
## WRIT -0.311061231  0.031378756 -0.056045596  0.137104995  0.305562842
## PHYS -0.280723624 -0.089037698 -0.154000444 -0.841266046  0.118424976
## RTEN -0.309790218  0.039381306  0.172869757 -0.184223629  0.006717911
##               PC6          PC7           PC8         PC9        PC10
## CONT -0.004064432  0.005214784 -6.006597e-02  0.02514533 -0.03038881
## INTG  0.369937339 -0.449810741  3.341645e-01  0.27537794  0.10897641
## DMNR -0.394724667  0.466747889 -2.470974e-01  0.19910004 -0.07241282
## DILG  0.598676072  0.209710731 -3.548587e-01 -0.03977180 -0.38339165
## CFMG -0.085728870  0.246903359  7.135261e-01 -0.14342471  0.09850310
## DECI -0.392609484 -0.536429933 -3.024227e-01  0.25823773  0.06743847
## PREP  0.083216652  0.335390036 -1.536754e-01  0.10876864  0.67986284
## FAMI -0.101538704 -0.036378004  2.038889e-02  0.22306628  0.04004599
## ORAL -0.150272440  0.057580177  9.062990e-02 -0.29951714 -0.25599455
## WRIT -0.238172386 -0.060899994  1.261203e-01 -0.02497324 -0.47478254
## PHYS  0.299281534  0.024959951 -1.364511e-05  0.26627286 -0.05900837
## RTEN  0.036126847 -0.256194180 -2.213898e-01 -0.75645893  0.24993250
##               PC11         PC12
## CONT  0.0145329260 -0.007940919
## INTG -0.1125535650  0.009848658
## DMNR  0.1343234234  0.059121657
## DILG  0.0709517642  0.053790339
## CFMG  0.1658680927  0.025082947
## DECI -0.1284999526  0.044141604
## PREP -0.3187612119 -0.273286884
## FAMI  0.5733628652  0.421739844
## ORAL -0.6386061655  0.494391025
## WRIT  0.0004056397 -0.696107204
## PHYS -0.0181381019 -0.053783960
## RTEN  0.2855143026 -0.080267574

Видно, что в первую главную компоненту все переменные, кроме CONT, входят с отриательным весом. Вес у CONT, хоть и маленький, но относительно других переменных более значительный. Поэтому можем считать, что первая главная компонента отвечает за “контактность” судьи, его общение с адвокатом. Во второй главной компоненте наибольший вес имеют открытость и поведение судьи, давайте считать эту компонентом индикатором “хороших манер” судьи. Третья главная компонента несильно отличается от второй, а вот четвертая интереснее. Можно сказать, что эта главная компонента отражает компетентность судьи: его добросовестность, подготовленность к судебному процессу, знание законов. Пятая главная компонента похожа на показатель “общего впечатления” (performance): умение делать внятные устные и письменные постановления, хорошие физические способности (способность долгое время находиться в тонусе во время заседания). На этом можно пока остановиться.

Главные коппоненты можно интерпретировать и дальше, но тут встает вопрос о целесообразности такого занятия. Ведь наша задача в методе главных компонент - снизить размерность, то есть уменьшить число переменных, с которыми мы будем работать. Поэтому, будет странно, если мы решим выбрать, например, 10 компонент вместо 12 переменных. Хочется же значительно сократить число исследуемых показателей и выбрать наиболее информативные…

Посмотрим на информативность (важность) каждой главной компоненты.

summary(pca)
## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6
## Standard deviation     3.1833 1.05078 0.57698 0.50383 0.29061 0.19310
## Proportion of Variance 0.8445 0.09201 0.02774 0.02115 0.00704 0.00311
## Cumulative Proportion  0.8445 0.93647 0.96421 0.98537 0.99240 0.99551
##                            PC7     PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.14030 0.12416 0.08851 0.07491 0.05708 0.04539
## Proportion of Variance 0.00164 0.00128 0.00065 0.00047 0.00027 0.00017
## Cumulative Proportion  0.99715 0.99844 0.99909 0.99956 0.99983 1.00000

Результаты представлены в виде сводных показателей. Первый - стандартное отклонение главной компоненты. То есть, возведя это значение в квадрат, мы узнаем, дисперсию показателя в новых координатах. Еще, стандартное отклонение здесь - это корень из собственного числа ковариационной матрицы, которое показывает, во сколько раз нужно растянуть векторы, соответствующие центрированно-нормированным данным, чтобы воспроизвести имеющуюся у нас структуру данных (вариация и связи переменных в базе m).

Второй показатель - доля объясненной дисперсии, которая позволяет получить ответ на следующий вопрос: какая доля дисперсии исходных данных объясняется конкретной главной компонентой. Третий показатель - кумулятивная (совокупная) доля объясненной дисперсии, то есть, какую долю дисперсии объясняет главная компонента в совокупности с предыдущими ГК.

Как определить, сколько главных компонент выбрать, если мы хотим не только привести содержательное обоснование, но и статистическое подтверждение? Во-первых, есть правило Кайзера: нужно выбирать столько главных компонент, сколько получается собственных значений больше 1. В нашем случае видно, что таких главных компонент две: PC1 (\(3.1833^2 \simeq 9\) и \(1.05078^2 \simeq 1\)).

Во-вторых, есть эмпирическое правило: нужно брать столько главных компонент, чтобы кумулятивный процент объясненной дисперсии был не менее 70%-80%. По этому правилу тоже получаются две первые главные компоненты.

Наконец, существует правило Кеттела, которое имеет логику, схожую с методом согнутого колена (локтя) в кластерном анализе. Нужно построить график, по горизонтальной оси которого отложены номера главных компонент, а по вертикальной – их дисперсии, и выбрать столько главных компонент, сколько находятся “до колена”. Такой график еще назвается scree plot:

plot(pca, type = "l", main = "Scree plot")

И вновь две главных компоненты. Что нам это дает? Знание о том, что для достаточно точного описания наших 12 показателей достаточно иметь две определенные комбинации этих показателей:

\[PC1 = 0.003 \cdot CONT -0.289 \cdot INTG - 0.286 \cdot DMNR - \dots - 0.310 \cdot RTEN \]

\[ PC2 = -0.933 \cdot CONT + 0.182 \cdot INTG + 0.198 \cdot DMNR + \dots + 0.039 \cdot RTEN\]

Эти две компоненты можно считать интегральными индексами. Первая - индекс общительности и вторая –индекс поведения (хороших манер) судьи. Можно рассчитать значения этих индексов для каждого судьи, чтобы потом ранжировать их (упорядочивать по возврастанию или убыванию), сравнивать средние значения по индексов по группам (если найдем демографическую информаию по этим людям), и так далее.

Определим значения каждой главной компоненты для наблюдений в базе (в терминах перемножения матриц мы просто умножим каждую строку в m на столбец из матрицы Rotation):

m2 <- predict(pca, newdata = m)

Посмотрим на объект m2:

head(m2)
##             PC1        PC2         PC3         PC4         PC5
## [1,]  0.5857508  1.8130995 -0.25849345 -0.88541465 -0.19810313
## [2,] -2.3800115  0.8678932  0.46099551 -0.05991797 -0.30448920
## [3,] -0.2251390  0.2987803  0.10084841  0.02221771 -0.13017390
## [4,] -3.6178993  0.5893672 -0.42073082  0.19730787  0.08296828
## [5,]  6.8711315 -0.1357495 -0.95544940  1.06915077 -0.30075003
## [6,] -2.4454791  1.4057858  0.08624784 -0.14887869 -0.22772257
##                PC6         PC7         PC8         PC9         PC10
## [1,] -0.0005757578 -0.09137858 -0.13526907 -0.12857586  0.119341425
## [2,]  0.2291020429 -0.14693431 -0.29068090  0.02688336  0.033482298
## [3,] -0.0150672307  0.03934473 -0.09174440 -0.06741825  0.021651215
## [4,]  0.0823548542 -0.04125003 -0.08852460  0.10683667  0.062934470
## [5,]  0.3662956835 -0.11568392  0.09722221  0.03493898 -0.031586001
## [6,]  0.2398656404  0.02902162 -0.11848363 -0.04460288 -0.009661878
##             PC11          PC12
## [1,] -0.01091350  0.0598665378
## [2,]  0.13678044  0.0064969215
## [3,]  0.06865240  0.0016395476
## [4,] -0.02381327 -0.0034183317
## [5,] -0.08729293  0.1160346423
## [6,] -0.02203656  0.0009017388

В m2 сохранены значения главных компонент для каждого наблюдения, то есть значения новых интегральных “индексов компетентности” для судей. Теперь посмотрим на корреляции между главными компонентами (“индексами компетентности”):

cor(m2)
##                PC1           PC2           PC3           PC4           PC5
## PC1   1.000000e+00 -3.144469e-16 -1.533744e-16 -5.415659e-16 -1.894434e-15
## PC2  -3.144469e-16  1.000000e+00  3.656561e-16 -3.736549e-17  1.004348e-16
## PC3  -1.533744e-16  3.656561e-16  1.000000e+00  1.235125e-16  9.875528e-16
## PC4  -5.415659e-16 -3.736549e-17  1.235125e-16  1.000000e+00  4.966328e-16
## PC5  -1.894434e-15  1.004348e-16  9.875528e-16  4.966328e-16  1.000000e+00
## PC6   3.269063e-16 -1.962836e-16  4.255859e-17  9.844700e-16  2.488013e-17
## PC7   1.668104e-15  3.115529e-16  2.824472e-15 -7.492448e-16  4.299730e-16
## PC8  -1.990301e-15  4.693585e-16 -1.846566e-15  3.567437e-16 -1.434178e-15
## PC9   4.392010e-15 -7.434676e-16 -1.761981e-15 -6.846485e-16  1.816207e-16
## PC10 -1.299895e-15 -6.423597e-17  1.134284e-15  1.060094e-15  9.588901e-16
## PC11  6.089771e-15 -6.559170e-16 -3.045706e-15 -1.473830e-15 -1.728038e-15
## PC12  1.186571e-15 -4.869973e-17  2.333852e-15  1.116960e-16  1.199785e-15
##                PC6           PC7           PC8           PC9          PC10
## PC1   3.269063e-16  1.668104e-15 -1.990301e-15  4.392010e-15 -1.299895e-15
## PC2  -1.962836e-16  3.115529e-16  4.693585e-16 -7.434676e-16 -6.423597e-17
## PC3   4.255859e-17  2.824472e-15 -1.846566e-15 -1.761981e-15  1.134284e-15
## PC4   9.844700e-16 -7.492448e-16  3.567437e-16 -6.846485e-16  1.060094e-15
## PC5   2.488013e-17  4.299730e-16 -1.434178e-15  1.816207e-16  9.588901e-16
## PC6   1.000000e+00 -1.902228e-16 -1.044185e-15 -1.812230e-15  8.156293e-16
## PC7  -1.902228e-16  1.000000e+00  9.807703e-16  1.354460e-15  3.514828e-16
## PC8  -1.044185e-15  9.807703e-16  1.000000e+00  1.209401e-15  1.845993e-16
## PC9  -1.812230e-15  1.354460e-15  1.209401e-15  1.000000e+00 -4.674462e-16
## PC10  8.156293e-16  3.514828e-16  1.845993e-16 -4.674462e-16  1.000000e+00
## PC11 -3.175632e-16 -6.032765e-18 -8.961713e-16  2.167627e-16 -1.171468e-15
## PC12 -4.097175e-16  8.988985e-16  1.548460e-15 -3.043606e-16 -1.003605e-15
##               PC11          PC12
## PC1   6.089771e-15  1.186571e-15
## PC2  -6.559170e-16 -4.869973e-17
## PC3  -3.045706e-15  2.333852e-15
## PC4  -1.473830e-15  1.116960e-16
## PC5  -1.728038e-15  1.199785e-15
## PC6  -3.175632e-16 -4.097175e-16
## PC7  -6.032765e-18  8.988985e-16
## PC8  -8.961713e-16  1.548460e-15
## PC9   2.167627e-16 -3.043606e-16
## PC10 -1.171468e-15 -1.003605e-15
## PC11  1.000000e+00 -7.467347e-16
## PC12 -7.467347e-16  1.000000e+00

Из теории мы знаем, что главные компоненты никак не связаны друг с другом - корреляция между ними равна 0. То же наблюдается и в R, просто из-за погрешностей мы получаем не “чистый” нуль, а значения, очень близкие к нулю.

Объект m2 - это матрица. Для удобства превратим матрицу в датафрейм:

m2 <- as.data.frame(m2)

Теперь сделаем следующее: в исходную базу данных m добавим столбец, который отвечает за значения первой главной компоненты. Так будет удобнее оенить корреляцию между исходными показателями и первой главной компонентой (корреляция между главными компонентами равна 0, но оценивать корреляцию между исходными переменными и их линейной комбинацией никто не запрещал).

with1 <- cbind(m, m2$PC1)
colnames(with1)[ncol(with1)] <-  "PC1"

Пояснения к коду: функция cbind() позволяет объединить два объекта по столбцам (c - от columns). Функция ncol() возвращает число столбцов в базе данных или матрице, а функция colnames() - названия столбцов. Здесь мы переименовываем последний столбец: берем последнее название - название с номером, совпадающим с числом столбцов в базе, и присваиваем ему значение "PC1".

Построим матрицу диаграмм рассеяния для новой базы:

ggpairs(with1)

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

Кроме того, можем явно получить корреляионную матрицу, без всяких графиков:

cor(with1)
##              CONT       INTG       DMNR       DILG       CFMG        DECI
## CONT  1.000000000 -0.1331909 -0.1536885  0.0123920  0.1369123  0.08653823
## INTG -0.133190890  1.0000000  0.9646153  0.8715111  0.8140858  0.80284636
## DMNR -0.153688533  0.9646153  1.0000000  0.8368510  0.8133582  0.80411683
## DILG  0.012391997  0.8715111  0.8368510  1.0000000  0.9587988  0.95616608
## CFMG  0.136912301  0.8140858  0.8133582  0.9587988  1.0000000  0.98113590
## DECI  0.086538232  0.8028464  0.8041168  0.9561661  0.9811359  1.00000000
## PREP  0.011469213  0.8777965  0.8558175  0.9785684  0.9579140  0.95708831
## FAMI -0.025636564  0.8688580  0.8412415  0.9573634  0.9354684  0.94280452
## ORAL -0.011996807  0.9113992  0.9067729  0.9544758  0.9505657  0.94825640
## WRIT -0.043810249  0.9088347  0.8930611  0.9592503  0.9422470  0.94610093
## PHYS  0.054248275  0.7419360  0.7886804  0.8129211  0.8794874  0.87176277
## RTEN -0.033643434  0.9372632  0.9437002  0.9299652  0.9270827  0.92499241
## PC1   0.009789154 -0.9185484 -0.9132432 -0.9688554 -0.9631848 -0.96101767
##             PREP        FAMI        ORAL        WRIT        PHYS
## CONT  0.01146921 -0.02563656 -0.01199681 -0.04381025  0.05424827
## INTG  0.87779650  0.86885798  0.91139915  0.90883469  0.74193597
## DMNR  0.85581749  0.84124150  0.90677295  0.89306109  0.78868038
## DILG  0.97856839  0.95736345  0.95447583  0.95925032  0.81292115
## CFMG  0.95791402  0.93546838  0.95056567  0.94224697  0.87948744
## DECI  0.95708831  0.94280452  0.94825640  0.94610093  0.87176277
## PREP  1.00000000  0.98986345  0.98310045  0.98679918  0.84867350
## FAMI  0.98986345  1.00000000  0.98133905  0.99069557  0.84374436
## ORAL  0.98310045  0.98133905  1.00000000  0.99342943  0.89116392
## WRIT  0.98679918  0.99069557  0.99342943  1.00000000  0.85594002
## PHYS  0.84867350  0.84374436  0.89116392  0.85594002  1.00000000
## RTEN  0.95029259  0.94164495  0.98213227  0.96755639  0.90654782
## PC1  -0.98493864 -0.97625799 -0.99544963 -0.99020634 -0.89363214
##             RTEN          PC1
## CONT -0.03364343  0.009789154
## INTG  0.93726315 -0.918548434
## DMNR  0.94370017 -0.913243218
## DILG  0.92996523 -0.968855390
## CFMG  0.92708271 -0.963184763
## DECI  0.92499241 -0.961017675
## PREP  0.95029259 -0.984938636
## FAMI  0.94164495 -0.976257987
## ORAL  0.98213227 -0.995449633
## WRIT  0.96755639 -0.990206339
## PHYS  0.90654782 -0.893632136
## RTEN  1.00000000 -0.986160301
## PC1  -0.98616030  1.000000000

Напоследок посмотрим на еще одну возможную визуализацию для МГК. Для нее нам потребуется библиотека ggbiplot. Ее лучше установить с Github таким образом:

install.packages("devtools")
library(devtools)
install_github("vqv/ggbiplot")

Обратимся к библиотеке и построим график:

library(ggbiplot)
ggbiplot(pca, scale = 0)

Пояснения к коду: pca - объект, где сохранены результаты МГК, а scale=0 означает, что

Что это за график? На графике изображены векторы, соответствующие переменным в исходной базе данных (их названия подписаны). По осям \(x\) и \(y\) идут две первые (наиболее информативные) главные компоненты. По этому графику удобно смотреть, как связаны исходные показатели и главные компоненты. Если вектор, соответствующий переменной, сонаправлен с осью, то между ними есть положительная связь, если противоположно направлены - отрицательная. Кроме того, угол между векторами тоже может сказать о многом: косинус угла между векторами - не что иное как коэффициент корреляции между показателями, соответствующими этим векторам. Поэтому, чем меньше угол, тем сильнее связь. Если угол равен 90 градусам (векторы перпендикулярны друг другу), то показатели совсем не связаны, коэффициент корреляции равен 0 (косинус 90 градусов равен 0). Если угол равен 45 градусам, то между показателями наблюдается сильная положительная корреляция (косинус 45 градусов примерно равен 0.7). Так что, можно смотреть на картинку и размышлять.

На этом всё.