Подготовка данных

Выберем индексы WGI и Freedom House – и на их основе будем кластеризовать данные):

df <- read.csv("wgi_fh.csv", dec = ",")
df <- na.omit(df)

library(dplyr)
dat <- df %>% select(va:fh)
rownames(dat) <- df$cnt_code

Во всех примерах будем использовать кластеризацию со следующими параметрами:

M <- dist(scale(dat))^2
hc <- hclust(M, method = "ward.D")

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

Как можно догадаться, дендрограмма для массива, состоящего из большого числа наблюдений, будет выглядеть не очень симпатично: метки для наблюдений (даже если это будут сокращенные названия или просто номера наблюдений) будут накладываться друг на друга, да и расстояние между “ветками” дендрограммы будет очень маленьким. Как быть? Есть несколько решений.

1. Построить дендрограмму сразу в pdf-файле, минуя окно для графиков в RStudio.

pdf("my_dend.pdf", width = 45, height = 20)
plot(hc)
dev.off()
## png 
##   2

Выглядеть это будет вполне прилично: см. полученный pdf-файл здесь.

Если вы не знаете, куда R сохранил pdf-файл. R сохраняет файлы в рабочую папку. Узнать, какая папка является рабочей (и где искать файл) можно так:

getwd() # wd - от working directory
## [1] "/home/oem/Рабочий стол"

2. Обрезать дендрограмму на определенном уровне и исследовать ее по частям.

(оригинальный ответ на stackoverflow – на нем основан код ниже)

Превратим объект типа hclust в дендрограмму:

hcd <- as.dendrogram(hc)
hcd # дендрограмма, только пока не построенная
## 'dendrogram' with 2 branches and 195 members total, at height 1426.79
plot(hcd) # не надо так

“Обрежем” дендрограмму на уровне 100 – проведем горизонтальную линию при y = 100 (по оси y отмечены, как мы помним, расстояние между кластерами).

# далее будут графики подряд, в 3 строки, 1 столбец
par(mfrow = c(3, 1))

# верхняя часть при обрезке
plot(cut(hcd, h = 100)$upper, main = "Верхняя часть дендрограммы")

# первая ветка нижней части (Branch 1)
plot(cut(hcd, h = 100)$lower[[1]], 
     main = "Branch 1")

# третья ветка нижней части (Branch 3)
plot(cut(hcd, h = 100)$lower[[3]], 
     main = "Branch 5")

Таким образом, мы можем обрезать дендрограмму на нужной нам высоте и “приближать” нужную нам ветку.

3. Приближать определенные части дендрограммы

Способ, аналогичный предыдущему, но здесь мы просто будем ограничивать значения по оси x и по оси y, которые должны быть на графике (xlim – границы окошка по горизонтали, ylim – границы окошка по вертикали).

hcd <- as.dendrogram(hc)
plot(hcd, xlim = c(1, 20), ylim = c(1, 10))

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

Например, так. Выделим 9 кластеров.

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

Создадим метки для наблюдений – номер кластера, в котором они находятся. Добавим столбец с метками в базу:

groups9 <- cutree(hc, k = 9) 
dat <- dat %>% mutate(groups9 = factor(groups9), 
                      country = df$country) # и название страны

Выберем из базы все страны, которые находятся в кластере номер 1:

dat %>% filter(groups9 == 1)
##       va   ps   ge   rq   rl   cc  fh groups9              country
## 1   1.20 1.40 1.86 0.87 1.56 1.23 1.0       1              Andorra
## 2   1.30 0.96 1.58 1.90 1.75 1.77 1.0       1            Australia
## 3   1.38 1.24 1.80 1.74 1.84 1.98 1.0       1               Canada
## 4   1.46 1.32 2.03 1.91 1.94 2.05 1.0       1          Switzerland
## 5   1.33 0.76 1.74 1.82 1.61 1.83 1.0       1              Germany
## 6   1.47 0.85 1.89 1.58 1.90 2.24 1.0       1              Denmark
## 7   1.49 0.96 1.85 1.82 2.02 2.28 1.0       1              Finland
## 8   0.27 0.84 1.86 2.15 1.70 1.58 3.5       1 Hong Kong SAR, China
## 9   1.34 1.33 1.41 1.28 1.51 1.99 1.0       1              Iceland
## 10  1.19 1.46 1.70 1.37 1.68 2.05 1.0       1        Liechtenstein
## 11  1.44 1.41 1.69 1.72 1.71 2.08 1.0       1           Luxembourg
## 12  1.48 0.89 1.84 1.98 1.89 1.95 1.0       1          Netherlands
## 13  1.58 1.17 1.88 1.70 2.02 2.20 1.0       1               Norway
## 14  1.44 1.49 1.86 2.04 1.93 2.30 1.0       1          New Zealand
## 15 -0.28 1.53 2.21 2.18 1.83 2.07 4.0       1            Singapore
## 16  1.50 0.98 1.79 1.85 2.04 2.22 1.0       1               Sweden

А теперь в кластере 5:

dat %>% filter(groups9 == 5)
##       va    ps    ge    rq    rl    cc  fh groups9           country
## 1  -0.62 -0.60 -0.15  0.25 -0.11 -0.57 4.5       5           Armenia
## 2  -1.45 -0.86  0.32  0.61  0.46 -0.06 6.5       5           Bahrain
## 3  -0.95  1.26  1.07  0.59  0.62  0.66 5.5       5 Brunei Darussalam
## 4  -0.05  0.99  0.49 -0.67  0.47  1.14 3.5       5            Bhutan
## 5  -1.62 -0.52  0.36 -0.26 -0.22 -0.25 6.5       5             China
## 6   0.22 -0.29  0.51  1.01  0.37  0.67 3.0       5           Georgia
## 7  -0.76 -0.53  0.14  0.05  0.31  0.27 5.5       5            Jordan
## 8  -1.29  0.04 -0.06 -0.10 -0.42 -0.80 5.5       5        Kazakhstan
## 9  -0.69 -0.15 -0.18 -0.07  0.03 -0.20 5.0       5            Kuwait
## 10 -0.65 -0.29 -0.10 -0.23 -0.14 -0.15 4.5       5           Morocco
## 11 -0.74  0.41 -0.33 -0.46 -0.41 -0.67 4.5       5          Maldives
## 12 -0.47  0.10  0.88  0.71  0.54  0.11 4.0       5          Malaysia
## 13 -1.11  0.80  0.19  0.61  0.43  0.37 5.5       5              Oman
## 14 -1.20  0.87  0.75  0.70  0.86  0.92 5.5       5             Qatar
## 15 -1.21 -0.05  0.11  0.11  0.07  0.69 6.0       5            Rwanda
## 16 -1.78 -0.50  0.24  0.08  0.47  0.23 7.0       5      Saudi Arabia
## 17  0.16  0.72  0.36 -0.26  0.13  0.79 3.0       5        Seychelles
## 18 -1.10 -0.93  0.34  0.17  0.01 -0.40 5.5       5          Thailand
## 19 -1.41  0.17  0.01 -0.45  0.05 -0.40 6.0       5           Vietnam

Подготовка данных 2

Чтобы не пришлось выполнять манипуляции, описанные в предыдущем разделе, выберем случайным образом 60 стран из базы данных для получения более читаемых дендрограмм:

set.seed(1234) # для воспроизводимости выделите эту строку и следующую и запустите
dat2 <- df[sample(nrow(df), 40), ]
d <- dat2 %>% select(va:fh)
rownames(d) <- dat2$cnt_code

Далее будем работать с урезанной версией базы.

m <- dist(scale(d))^2
hc2 <- hclust(m, method = "ward.D")
plot(hc2)

Как построить дендрограмму другой формы

Очень просто! Взять объект типа dendrogram и при построении добавить в plot() аргумент type.

hcd2 <- as.dendrogram(hc2)
plot(hcd2, type = "triangle")

Более радикальные варианты – с помощью библиотеки ape (Analyses of Phylogenetics and Evolution).

# install.packages("ape")
library(ape)

Например, круглая:

ph <- as.phylo(hc2)
plot(ph, type = "fan")

Или радиальная:

ph <- as.phylo(hc2)
plot(ph, type = "radial")

Как выделить кластеры на дендрограмме цветом

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

groups4 <- cutree(hc2, k = 4)
colors = c("red", "blue", "green", "magenta")

# строим график, выставляя tip.color и размер шрифта cex
plot(ph, tip.color = colors[groups4], cex = 0.6)

Как сравнить на дендрограмме ожидаемые кластеры с получившимися

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

В имеющейся базе у нас нет столбца с такими категориями стран, создадим его:

d <- d %>% mutate(status = cut(fh, breaks = c(-Inf, 3, 5.5, Inf), 
             labels = c('free', 'partly free', 'not free'), right = F))
head(d)
##      va    ps    ge    rq    rl    cc  fh      status
## 1  0.47 -0.45 -0.18 -0.21 -0.08 -0.44 2.0        free
## 2 -0.83 -0.74 -0.79 -0.74 -0.78 -0.80 5.5    not free
## 3  0.45  0.82 -0.11 -0.08 -0.22 -0.50 1.5        free
## 4 -0.39 -1.05 -0.85 -0.70 -1.02 -0.87 4.0 partly free
## 5  0.16  0.72  0.36 -0.26  0.13  0.79 3.0 partly free
## 6  0.86  1.05  0.96  1.03  0.80  0.32 1.5        free
to_clust <- d %>% select(va:fh)
rownames(to_clust) <- dat2$cnt_code
m <- dist(scale(to_clust))^2
hc3 <- hclust(m, method = "ward.D")
hcd3 <- as.dendrogram(hc3)

Выберем три цвета:

cols <- c("green", "orange", "red")
cols_vect <- cols[d$status]

Установим и загрузим библиотеку dendextend. Работать с ней можно так же, как и с dplyr, используя оператор %>%.

install.packages("dendextend")
library(dendextend)

Далее hcd3 – это наш объект типа dendrogram (результат функции as.dendrogram()). Выставляем цвета меток наблюдений labels_col, цвета веток дендрограммы branches_k_color, тип точек-листьев дендрограммы leaves_pch, размер точек-листьев дендрограммы nodes_cexи их цвета nodes_col. Цвета в nodes_col – те цвета, которые мы установили в соответствии с ожидаемыми тремя группами стран: free, partly free, not free.

hcd3 %>%
  set("labels_col", value = c("skyblue", "plum", "grey"), k = 3) %>%
  set("branches_k_color", value = c("skyblue", "plum", "grey"), k = 3) %>%
  set("leaves_pch", 19)  %>% set("nodes_cex", 0.7) %>% set("nodes_col",  cols[d$status][order.dendrogram(hcd3)]) %>%
  plot()

Три кластера стран выделить можно, но сказать, что эти три кластера соответствуют делению по признаку свободы (free, partly free, not free), нельзя: по дендрограмме видно, что точки разных цветов расположены в разных кластерах. Например, зеленые точки (свободные страны) разбросаны по всем трем кластерам.

Что еще почитать