Часть 1

Реализуем в R кластерный анализ, который мы проделали на лекции вручную. У нас есть пять наблюдений и по ним известны значения \(x_1\) и \(x_2\).

x1 <- c(5, 5, 3, 2, 3.5)
x2 <- c(2, 1, 2, 1, 4)

“Склеим” их в маленькую базу:

dat <- as.data.frame(cbind(x1, x2))
View(dat)

Теперь посчитаем расстояния между всеми парами точек (по умолчанию функция dist() считает евклидово расстояние) – получим матрицу расстояний и назовем ее Mdist:

Mdist <- dist(dat)

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

Mdist
##          1        2        3        4
## 2 1.000000                           
## 3 2.000000 2.236068                  
## 4 3.162278 3.000000 1.414214         
## 5 2.500000 3.354102 2.061553 3.354102

На главной диагонали у этой матрицы ничего нет: зачем выводить значения, если и так ясно, что расстояние от точки до самой себя равно нулю? К тому же матрица симметрична (по свойству метрики, расстояние от точки A до точки B – это то же, что расстояние от точки B до точки A).

Теперь реализуем иерархический кластерный анализ (используем функцию hclust()) и подадим на вход этой функции нашу матрицу расстояний. В качестве способа агрегирования возьмем метод ближнего соседа, он же метод одиночной связи (method = "single").

hc <- hclust(Mdist, method = "single")

Построим дендрограмму:

plot(hc)

Дендрограмма получилась точно такой же, как и в нашей “ручной” реализации кластерного анализа на доске!

Маленький, но очень важный момент: вообще данные перед тем, как получить матрицу расстояний, нужно шкалировать, приводить к единой шкале. Для этого существует функция scale(). Она работает так: из каждого значения переменной вычитает среднее значение этой переменной и делит результат на ее стандартное отклонение. Проделаем ту же кластеризацию, но со шкалированными переменными:

Mdist1 <- dist(scale(dat))
hc1 <- hclust(Mdist1, method = "single")
plot(hc1)

Как можно заметить, результат глобально не изменился. Но так будет не всегда, просто у нас пример игрушечный. Шкалы у переменных могут быть разными и разброс значений тоже (например, явка на выборы, измеренная в процентах, может принимать значения от 0 до 100, а рейтинг губернатора в регионе от 1 до 10).

Часть 2

Поработаем с более серьезными данными. Загрузим базу данных, в которой хранятся следующие показатели за 2016 год:

Так:

df <- read.csv("wgi_fh.csv", dec = ",") # или вариант внизу

Или так (с file.choose):

df <- read.csv(file.choose(), dec = ",")

Не забудьте dec = "," – это разделитель в дробных числах: без него дробные числа с запятой (2,5) будут считываться как текст, и ничего хорошего мы не получим.

Посмотрим на базу:

View(df)

Удалим из нее пропущенные значения:

df <- na.omit(df)

Так как на первом занятии наша задача – познакомиться с кластерным анализом, а не заниматься настройкой параметров для красивых графиков, чтобы графики были читаемыми, выберем случайным образом 50 стран из базы (дендрограмма для всех 195 стран будет не очень наглядной).

set.seed(1234) # для воспроизводимости выделите эту строку и следующую и запустите
data <- df[sample(nrow(df), 50), ]

Посмотрим на нашу маленькую базу данных:

View(data)

Теперь отберем только те столбцы, которые нужны нам для кластерного анализа. Мы будем определять кластеры стран по значениям индексов-компонентов WGI и Freedom House. Обратимся к библиотеке dplyr и выберем соотвествующие столбцы. А потом сохраним их в базу d.

library(dplyr)
# от va до fh подряд - через :
d <- data %>% select(va:fh)

Назовем строки в базе данных по коду страны (аббревиатура). Зачем такие сложности? Почему нельзя просто добавить отдельный столбец с кодом страны? Этот столбцец будет текстовым, и он будет мешать, когда мы будем подавать базу d на вход функции dist().

rownames(d) <- data$cnt_code

Полюбуемся на базу еще раз:

View(d)

Теперь все готово к работе! В качестве метода агрегирования выберем метод Варда (Ward) – он эффективно работает на данных, измеренных в количественной шкале. Единственное, у него есть небольшой недостаток: он склонен создавать кластеры маленького размера (из малого числа наблюдений).

Метод Варда особенный: он требует, чтобы в качестве метрики был использован квадрат евклидова расстояния. Реализуем это! (Просто возведем все значения в обычной матрице расстояний в квадрат).

# матрица расстояний
M <- dist(scale(d))^2

А теперь реализуем сам иерархический кластерный анализ, используя метод Варда:

hc <- hclust(M, method = "ward.D") # пока ничего не произошло

Построим дендрограмму (а заодно поправим шрифт для меток – названий стран).

plot(hc, cex = 0.6) # cex = 0.6 - размер шрифта

Дендрограмма показывает нам все возможные кластеры, которые можно найти в наших данных. Сколько кластеров взять – на усмотрение исследователя и исходя их содержательных соображений. Включите в себе исследователя и подумайте: сколько кластеров можно выделить на основе этой дендрограммы?

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

Выделим их на дендрограмме:

# main - заголовок графика
plot(hc, cex = 0.6, main = "2 clusters") 
rect.hclust(hc, k = 2, border="red") # 2 кластера

plot(hc, cex = 0.6, main = "4 clusters")
rect.hclust(hc, k = 4, border="red") # 4 кластера

Случай с четырьмя кластерами интересный. Давайте его рассмотрим. В первом кластере (идем справа налево) находятся следующие страны: Эквадор, Коморы, Сьерра-Леоне, Украина, Кот д’Ивуар, Нигер, Филиппины, Кения, Мозамбик, Турция, Мексика. Во втором: Афганистан, Йемен, Саудовская Аравия, Таиланд, Иордания, Кувейт, Лаос, Куба, Мавритания, Республика Конго, Эритрея, Узбекистан. (Дальше можете посмотреть сами).

Кажется, что кластеры получились достаточно логичными. Посмотрим на них повнимательнее. Вытащим из полученного разбиения на кластеры метки для наблюдений, чтобы было ясно, какие страны в одном кластере, а какие – в разных. Воспользуемся функцией cutree().

groups4 <- cutree(hc, k = 4) 
groups4 # посмотрим на метки
## BRA MRT MNG MOZ SYC MUS AFG DEU NAM LAO NLD LIE ERI THA ESP SAU USA DOM 
##   1   2   1   3   1   4   2   4   1   2   4   4   2   2   4   2   4   1 
## CIV CUB FIN UKR BWA ATG COM PER KEN SGP PHL AUS HUN VUT ECU ISR CAN YEM 
##   3   2   4   3   1   1   3   1   3   4   3   4   1   1   3   4   4   2 
## CHL TUR SLE NER JOR LCA TWN KWT SWE SRB UZB HKG COG MEX 
##   4   3   3   3   2   1   4   2   4   1   2   4   2   3

Теперь добавим столбец с метками для кластеров в нашу базу d с помощью функции mutate() из библиотеки dplyr. Создадим в ней столбец groups4 типа “factor” (так как это условные метки, качественный показатель, не числа 1, 2, 3, 4). Заодно добавим столбец с кодами стран.

d <- d %>% mutate(groups4 = factor(groups4), country = data$cnt_code)
View(d)

А теперь будем строить диаграммы рассеяния для пар показателей и выделять точки, относящиеся к разным кластерам, разным цветом. Для начала возьмем Freedom House и Voice and Accountability.

library(ggplot2)

ggplot(data = d, aes(x = fh, y = va, color = groups4)) + geom_point()

Пока не будем наводить красоту на графике, только подпишем точки – добавим через + слой geom_text():

# vjust и hjust - чтобы подписи были чуть в стороне и не закрывали точки
ggplot(data = d, aes(x = fh, y = va, color = groups4)) + geom_point() +
  geom_text(aes(label = country, vjust = 0, hjust = 0))

Глядя на эту диаграмму рассеяния можно сказать, что деление на четыре кластера вполне удалось: точки разного цвета образуют сформированные, достаточно “плотные” группы. Но есть исключения: Сингапур и Гонконг (точки фиолетового цвета среди голубых точек), а также Санта-Лучия Сент-Лусия (розовая точка среди фиолетовых). Пока сложно судить, но, наверно, не стоит говорить о том, что мы поделили наблюдения на кластеры плохо: Сингапур и Гонконг во многом достаточно нетипичные страны, поэтому вряд ли они ровно лягут в любое деление на группы (если, конечно, не создать отдельную группу для них двоих).

Посмотрим на другие показатели: Freedom House и Control of Corruption.

ggplot(data = d, aes(x = fh, y = cc, color = groups4)) + geom_point() +
  geom_text(aes(label = country, vjust = 0, hjust = 0))

Не считая нескольких исключений (все те же Сингапур, Гонконг, Сент-Лусия), группы получились довольно однородные.

Возьмем Rule of Law and Government Effectiveness:

ggplot(data = d, aes(x = rl, y = ge, color = groups4)) + geom_point() +
  geom_text(aes(label = country, vjust = 0, hjust = 0))

Вот тут уже не так однозначно, нужно думать, почему страны так странно перемешались.

А давайте сгруппируем значения по каждому кластеру и выведем их характеристики (группировка - функция group_by()):

# summarise_at: вывести значения по переменным, указанным в vars
# mean - вывести средние
d %>% group_by(groups4) %>% summarise_at(vars(va:fh), mean)
## # A tibble: 4 x 8
##   groups4         va         ps         ge         rq         rl
##    <fctr>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
## 1       1  0.4658333  0.5158333  0.0225000  0.1066667  0.1408333
## 2       2 -1.3675000 -0.7550000 -0.6000000 -0.8641667 -0.6716667
## 3       3 -0.2072727 -0.9636364 -0.5445455 -0.4518182 -0.6636364
## 4       4  1.0326667  0.8080000  1.5780000  1.6213333  1.5400000
## # ... with 2 more variables: cc <dbl>, fh <dbl>
# summarise_at: вывести значения по переменным, указанным в vars
# median - вывести медиану
d %>% group_by(groups4) %>% summarise_at(vars(va:fh), median)
## # A tibble: 4 x 8
##   groups4     va     ps     ge     rq     rl     cc    fh
##    <fctr>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>
## 1       1  0.435  0.715  0.050 -0.005  0.225 -0.010  2.00
## 2       2 -1.395 -0.550 -0.495 -0.955 -0.775 -0.865  6.25
## 3       3 -0.180 -1.050 -0.580 -0.430 -0.670 -0.670  3.50
## 4       4  1.100  0.930  1.700  1.740  1.680  1.770  1.00
# summarise_at: вывести значения по переменным, указанным в vars
# min - вывести минимум
d %>% group_by(groups4) %>% summarise_at(vars(va:fh), min)
## # A tibble: 4 x 8
##   groups4    va    ps    ge    rq    rl    cc    fh
##    <fctr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1       1  0.16 -0.45 -0.88 -0.29 -0.49 -0.78     1
## 2       2 -2.10 -2.79 -1.82 -2.19 -1.62 -1.67     5
## 3       3 -0.63 -2.00 -1.54 -1.05 -1.13 -0.90     3
## 4       4 -0.28 -0.83  0.96  1.01  0.80  0.32     1
# summarise_at: вывести значения по переменным, указанным в vars
# max - вывести максимум
d %>% group_by(groups4) %>% summarise_at(vars(va:fh), max)
## # A tibble: 4 x 8
##   groups4    va    ps    ge    rq    rl    cc    fh
##    <fctr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1       1  1.09  1.09  0.51  0.60  0.52  0.93     3
## 2       2 -0.69  0.62  0.34  0.17  0.47  0.27     7
## 3       3  0.14 -0.02  0.14  0.29 -0.16 -0.20     4
## 4       4  1.50  1.53  2.21  2.18  2.04  2.28     4