Простой пример кластерного анализа алкогольных предпочтений по странам на R

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

Сразу оговорюсь, эта статья ни в коем случае не претендует на академическую полноту, уникальность полученных результатов или полноту освещения вопроса. Статья призвана продемонстрировать основные шаги классического кластерного анализа, которые могут быть использованы для простого и осмысленного (возможно, предваряющего более детальное) исследования. Любые исправления, замечания и дополнения по существу приветствуются.

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

В то же время каждая страна условно относится к одной из географических групп: восток, центр и запад. Разделение весьма условное и очень спорное по разным причинам, но будем отталкиваться от того, что есть. Источник данных — Global status report on alcohol and health 2014, С. 289-364

(Раскрашивал вручную, могут быть погрешности, но общая мысль, думаю, понятна)

Предварительный анализ

Подключим используемые библиотеки.

library(rgl)
library(heplots)
library(MVN)
library(klaR)
library('Morpho')
library(caret)
library(mclust)
library(ggplot2)
library(GGally)
library(plyr)
library(psych)
library(GPArotation)
library(ggpubr)

Загрузим данные и посмотрим, как они выглядят.

# Загружаем данные из файла
data <- read.table("alcohol_data.csv", header=TRUE,  sep=",")
# Делаем из первого столбца названия строк
rownames(data) <- make.names(data[,1], unique = TRUE)
# Записываем данные и удаляем пропуска, если они есть
data <- data[,-1]
data <- na.omit(data)
# Выводим несколько первых строчек
head(data)

Посчитаем стандартные статистики

summary(data)

Глядя на эти статистики, можно многое сказать о характере данных. Например, что у Other очень большой разброс, и максимум сильно отстоит от третьего квартиля, значит, там есть как минимум одно сильно отличающееся наблюдение, так называемый выброс. Также видно по первому квартилю и среднему, что в этом столбце существенное количество нулевых значений, что, скорее всего, связано с недостатком данных. Также сразу видим сколько стран в каждой из групп, группы не сбалансированны по количеству. Аналогично можно рассматривать и статистики по другим признакам и делать какие-то полезные выводы и предположения.

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

options(rgl.useNULL=TRUE)
open3d()
mfrow3d(2,2)
levelColors <- c('west'='blue', 'east'='red', 'center'='yellow')
plot3d(data$Beer, data$Wine, data$Spirit, xlab="Beer", ylab="Wine", zlab="Spirit", col = levelColors[data$Group], size=3)

widget <- rglwidget()
widget

Получаются довольно наглядные графики, которые можно покрутить и помасштабировать интерактивно. Судя по этому графику понятно, что по этим трем признакам группы визуально различаются.

Посмотрим эмпирические плотности по группам

ggpairs(
  data,
  mapping = ggplot2::aes(color = data$Group),
  upper = list(continuous = wrap("cor", alpha = 0.5), combo = "box"),
  lower = list(continuous = wrap("points", alpha = 0.3), combo = wrap("dot", alpha = 0.4)),
  diag = list(continuous = wrap("densityDiag",alpha = 0.5)),
  title = "Alcohol"
)

Так как Average и Total сильно коррелируют, исключим из рассмотрения Average.

data <- data[, -6]

У нас несколько групп, и не только предполагается, а даже видно, что они разные, поэтому нужно рассматривать распределение отдельно по группам. Посмотрим отличающиеся данные по группам.

data[data$Wine>60,]

В том что итальянцы пьют вина больше всех, даже без учета разделения на группы, думаю, нет ничего удивительного, поэтому из-за того, что и так мало данных, оставим это наблюдение.

data[data$Spirit>70,]
data[data$Spirit<10,]

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

Данных не так много, поэтому посмотрим на данные по группам

split(data[,1:5],data$Group)

$center

$east

$west

ggpairs(
  data,
  mapping = ggplot2::aes(color = data$Group),
  diag=list(continuous="bar", alpha=0.4)
)

Из гистограмм видно, что группы лучше всего различаются по пиву, вину и спиртным напиткам. Если брать признак Other, заметим: мало того, что у него весьма странное распределение, так еще и в большинстве стран маленькие и даже нулевые значения, кроме нескольких (примерно 10-12 стран, а так как индивидов всего 45, исключив их, мы существенно сократим объем данных). Больше всего отличается Республика Беларусь. Немного подумав, можно сделать предположение, что там предпочитают крамбамбуля вместо пива, вина и спиртных напитков (шутка). Ну, а если более формально, то скорее всего это связано с неполнотой исходных данных или несовершенством регистрации данных. Уберем признак Other из рассмотрения.

Также на гистограммах можно заметить, что для центра превалирует пиво, для запада — вино, а для востока — спиртные напитки. Это вполне укладывается в общеизвестные представления, можно даже сказать — стереотипы, о культуре потребления спиртных напитков в этих регионах.
На диаграммах рассеяния признака Total и Other, визуально группы не выделяются. Будем иметь этот факт в виду в дальнейшем.

Примечательно, что между признаками Beer, Spirit и Wine отрицательные корреляции. Возможно, это также относится к тому, что по этим переменным можно выделять группы предпочтения в алкоголе, и они будут близки к географическим. После того как изучили данные, получили некие априорные представления, убрали лишние, на наш взгляд, признаки, перейдем к кластерному анализу.

Кластерный анализ

Уберем разметку данных на группы и уберем признак Total. Это позволит не нормировать данные, так как остальные признаки — в одной шкале.

data.group = data[,5]
data <- data[,-5]
data<- data[,-4]

Определим число кластеров Elbow method (“метод согнутого колена”, он же “метод каменистой осыпи”). Построим график, где по оси абсцисс отмечено число кластеров k, а по оси ординат – значения функции W(K), которая определяет внутригрупповой разброс в зависимости от числа кластеров.

library(factoextra)
fviz_nbclust(data, kmeans, method = "wss") +
  labs(subtitle = "Elbow method") +
  geom_vline(xintercept = 4, linetype = 2)

data.dist <- dist((data))
hc <- hclust(data.dist, method = "ward.D2")
plot(hc, cex = 0.7)

Нарисуем график поинтереснее. Цвета будут означать исходную географическую маркировку.

colors=c('green', 'red', 'blue')
hcd = as.dendrogram(hc)
clusMember = cutree(hc, 4)
colLab <- function(n) {
    if (is.leaf(n)) {
        a <- attributes(n)
        labCol <- colors[data.group[n]]
        attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
    }
    n
}
clusDendro = dendrapply(hcd, colLab)
plot(clusDendro, main = "Cool Dendrogram", type = "triangle")

rect.hclust(hc, k = 4)

Грузия из восточного кластера единственная не попала в свой географический кластер. Пока отложим интерпретацию, посмотрим на другие методы.
Причем здесь, наверное, лучше использовать три кластера, так как в четвертый кластер странно выделились всего 4 страны.

plot(clusDendro, main = "Cool Dendrogram", type = "triangle")
data.hclas_group <- factor(cutree(hc, k = 3))

rect.hclust(hc, k = 3)

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

library(FactoMineR)
res.pca <- PCA(data,scale.unit=T, graph = F)
fviz_pca_biplot(res.pca, 
                col = colors[data.hclas_group], palette = "jco", 
                label = "var",
                ellipse.level = 0.8,
                 addEllipses = T,
                col.var = "black",
                legend.title = "groups4")

Третья группа, которая по центру, портит всю картину. Если бы она была покучнее, можно было говорить о кластеризации, а так это, скорее, сегментация. В целом, между группами заметно различие, посмотрим, как справится метод k-средних++.

library(flexclust)
data.kk <- kcca(data, k=3, family=kccaFamily("kmeans"),
control=list(initcent="kmeanspp"))

fviz_pca_biplot(res.pca, 
                col.ind =as.factor(data.kk@cluster), palette = "jco", 
                label = "var",
                ellipse.level = 0.8,
                 addEllipses = T,
                col.var = "black", repel = TRUE,
                legend.title = "clusters")

Эллипсоид третьей группы получился слишком широким, так как k-средних пытается разбить на равные кластеры. Из графика видно, что метод включил в третий кластер лишние точки, которые скорее должны относиться к другим кластерам.

Больше всего, исходя из исходных классов, на правду похож hclust. Проинтерпретируем результат.

Наиболее полный кластер, в смысле исходной классификации, образовали восточные страны. К ним попала Албания из центральной части. Стоит отметить, что Албания пространственно находится недалеко от восточных стран.

Второй кластер также очень похож на исходное разделение. В третьем же все плохо. В исходной выборке там меньше всего индивидов, и они очень близки по предпочтениям в алкоголе, согласно данным, к центральной части. Можно попробовать условно разделить на два кластера, так как видно, что для интерпретации лучше всего так и сделать и резюмировать биполярность Европы. Тогда кластера практически совпадут с Восточной и Западной Европой, где в Западную войдет центральная и Западная по исходным обозначениям.

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