Простой пример кластерного анализа алкогольных предпочтений по странам на 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. Проинтерпретируем результат.
Наиболее полный кластер, в смысле исходной классификации, образовали восточные страны. К ним попала Албания из центральной части. Стоит отметить, что Албания пространственно находится недалеко от восточных стран.
Второй кластер также очень похож на исходное разделение. В третьем же все плохо. В исходной выборке там меньше всего индивидов, и они очень близки по предпочтениям в алкоголе, согласно данным, к центральной части. Можно попробовать условно разделить на два кластера, так как видно, что для интерпретации лучше всего так и сделать и резюмировать биполярность Европы. Тогда кластера практически совпадут с Восточной и Западной Европой, где в Западную войдет центральная и Западная по исходным обозначениям.
Можно было бы провести кластеризацию на основе предположения о моделях кластеров, используя информационные критерии (тут описание), а также попробовать классический дискриминантный анализ для этого набора данных. Если эта статья была полезной, то планирую опубликовать продолжение.
