Анализ личности клиентов II: Кластерный анализ и ранжирование клиентов

Машинное обучение — это широкая область исследований, которую можно условно разделить на несколько подобластей, таких как обучение с учителем, обучение без учителя, обучение с подкреплением и так далее. Большинство связанных книг посвятили большую часть своего содержания обучению с учителем, потому что это хорошо изученная область. Методы контролируемого обучения часто имеют четкую цель (целевая функция или функция стоимости), а также систематические способы оценки производительности обученных моделей в реальном времени (например, точность, полнота, воспроизводимость и т. д.). С другой стороны, неконтролируемое обучение относится к методам, которые выявляют шаблон данных без какой-либо цели или переменной отклика. Неконтролируемое обучение, как правило, более субъективно и в большинстве случаев является частью исследовательского анализа данных¹.

Как одна из ветвей неконтролируемого обучения, кластерный анализ (CA) (также известный как кластеризация или сегментация данных) направлен на обнаружение скрытых подгрупп или сегментов в данных, чтобы точки данных в одной и той же подгруппе (кластере) были более похожи и тесно связаны с друг друга по сравнению с точками данных, назначенными другим кластерам. Понятие сходства/несходства/близости между парами данных, являющееся неотъемлемой частью алгоритма CA, обычно измеряется расстоянием. Чем больше расстояние между двумя экземплярами, тем ниже мера сходства. В нашем случае CA используется для разделения клиентов на несколько отдельных однородных групп с целью, чтобы различные маркетинговые кампании и рекламные объявления могли быть адаптированы для каждой группы клиентов.

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

Загрузить пакеты библиотек

library(tidyverse)
library(e1071)   # for skewness calculation
require(cluster)  # gower distance matrix and pam algorithm
library(fpc)      # cluster evaluation metrics
library(mdw)      # entropy-based feature weights
library(MCDA)      # TOPSIS
library(FactoMineR)  # Factor analysis for mixed data (FAMD)
library(factoextra)  # Visualization of factor analysis
library(knitr)       # For beatiful table display
library(car)         # For interactive 3D scatter plot
library(kableExtra)  # Custom html table styles

Импортировать предыдущие сохраненные данные

# set the working directory
setwd("~/ml projects/customer segmentation")
cust_data = readRDS("preprocessed_cust_data.rds")

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

Вычисление матрицы различий с использованием меры расстояния Гауэра

Расстояние Гауэра измеряет несходство двух точек данных со смешанными числовыми и категориальными данными. Я не буду здесь говорить о формуле расстояния Гауэра, но простой пример расчета расстояния Гауэра можно найти здесь.

Для смешанных типов данных можно использовать функцию daisy() из пакета cluster для вычисления попарного различия между наблюдениями. Аргумент type в функции daisy используется для указания типов переменных во входных данных cust_data_without_ID. Дополнительную информацию о функции daisy можно получить из функции help() или ввести команду ?daisy в консоли.

idx_col_retain = !(colnames(cust_data) %in% c("ID", "min_num_household", "tot_AcceptedCmp"))
cust_data_without_ID = cust_data[,idx_col_retain]

# Boolean attributes
seq_binary = c("Complain", "Response", "Accepted")
idx_binary = sapply(seq_binary, 
                    function(x) grep(x, colnames(cust_data_without_ID)))
idx_binary = unlist(idx_binary)

# continuos attributes
cont_features_patterns = c("Mnt", "Num", "Income","Recency", "age", 
                           "days_enroll") 
idx_cont = sapply(cont_features_patterns, 
                  function(x) grep(x, colnames(cust_data_without_ID)))
idx_cont = unlist(idx_cont)

skewness_col = apply(cust_data_without_ID[, idx_cont], 2, skewness)
idx_logtrans = idx_cont[which(abs(skewness_col)>1)]

# Ordinal attributes

dissimilarity_matrix = daisy(cust_data_without_ID, metric = "gower", 
                             type = list(ordratio=grep("home", colnames(cust_data)),
                                         asymm = idx_binary,
                                         logratio = idx_logtrans))

Кластеризация K-medoid (алгоритм PAM)

Существуют различные подходы к кластеризации, например, на основе разделения, на основе иерархии, на основе плотности и на основе моделей.

Подобно популярной кластеризации k-средних, кластеризация K-medoid представляет собой метод групповой кластеризации, в котором точки данных разделяются на группы или так называемые кластеры. Гиперпараметр, количество кластеров выбирается пользователем. Однако k-medoid, как правило, более надежен и менее чувствителен к выбросам, чем кластеризация k-средних. Почему? Ответы двоякие:

k-medoid выбирает медоиды (фактические точки данных) в качестве центроида (т. е. центра кластера, тогда как k-means выбирает среднее значение в качестве центроида. Это похоже на аналогию между медианой и средним значением, где медиана более надежна, чем среднее, как значение разбивки Медиана равна 0,5, тогда как точка разбивки среднего равна 1 / n. Вы можете рассматривать значение разбивки как сопротивление выбросам. Чем оно выше, тем надежнее статистическая мера.
Если мы посмотрим на функции стоимости обоих методов кластеризации, как показано ниже, мы заметим, что функция стоимости k-средних фактически находится в пределах суммы квадратов кластера, в то время как функция расстояния d(x, z) для k-medoid является произвольной. На функцию стоимости L2-нормы могут влиять выбросы.

В этом исследовании для кластеризации k-medoid выбран алгоритм разбиения на основе алгоритма medoids (PAM)². Алгоритм оказывается очень похожим на алгоритм k-средних: 1) Инициализация: случайным образом выбирает k точек данных в качестве медоидов, m 2) Назначает все немедоидные точки, o к ближайшим медоидам, 3) Для каждого m и o поменять местами m и o. Повторяйте шаги 2–3, пока функция стоимости не перестанет уменьшаться.

Выбор оптимального количества кластеров, k

После принятия решения о том, какой метод кластеризации использовать, теперь возникает вопрос, как мы можем измерить «качество» кластеров, сгенерированных для разных k? Ответ заключается в индексах проверки кластера (внутренних или внешних).

Внутренние: типичные целевые функции при кластеризации формализуют цель достижения высокого сходства внутри кластера (компактность) и низкого сходства между кластерами (разделение)³. Сходство — это числовая мера того, насколько похожи два экземпляра данных, как правило, измеряются мерами расстояния (евклидово расстояние). Чем больше расстояние между двумя точками данных, тем меньше сходство.
Внешний: должна быть доступна основная истина или золотой стандарт для оценки того, насколько хорошо результаты кластеризации совпадают с метками основной истины.

Поскольку данные не имеют аннотированных меток, у нас остаются внутренние метрики оценки. Чтобы выбрать лучший k, я буду использовать внутренние метрики оценки кластера, найденные в пакете fpc, такие как сумма квадратов внутри кластера (метод локтя), средняя ширина силуэта, Индекс Калински-Харабаша (CH) и Индекс Данна.

Наилучшее k можно выбрать из изгиба графика зависимости суммы квадратов от k (точка, в которой сумма квадратов начинает неуклонно уменьшаться); для случая средней ширины силуэта, индекса CH и индекса Данна мы пытаемся максимизировать индексы проверки и, таким образом, в диапазоне экспериментальных значений k наилучшее значение k достигается при наибольшем значении этих индексов. Я не показываю здесь формулу, но для ясности ключевой вывод заключается в том, что средняя ширина силуэта, индекс CH и индекс Данна являются скалярными значениями, полученными с учетом попарного расстояния между точками данных. Для заинтересованных читателей, желающих копнуть глубже, формулы этих индексов приведены в этой бумаге.

# Possible number of clusters are assumed to be in the range of 2 to 8. 
k_array = seq(from = 2, to = 8, by=1)
cluster_eval_df = data.frame(matrix(, nrow = length(k_array), ncol = 4))
colnames(cluster_eval_df) = c("silhouette", "CH_index", "Dunn_index", "wc_SOS")
cluster_eval_df$k = k_array

for (i in (1:length(k_array))){
  set.seed(i+100)
  kmedoid = pam(dissimilarity_matrix, k = k_array[i], diss = TRUE,
                nstart = 10)
  # set diss to TRUE, and set the number of random start as 10.
  clust_stat = cluster.stats(dissimilarity_matrix, 
                             clustering = kmedoid$clustering)
  cluster_eval_df[i,"silhouette"] = clust_stat$avg.silwidth
  cluster_eval_df[i, "CH_index"] = clust_stat$ch
  cluster_eval_df[i, "Dunn_index"] = clust_stat$dunn
  # Add in the within cluster sum of squares
  cluster_eval_df[i, "wc_SOS"] = clust_stat$within.cluster.ss
}

# Line plot for all evaluation metrics (internal cluster evaluation)
cluster_eval_df %>% 
  pivot_longer(!k, names_to = "cluster_eval", values_to = "value") %>% 
  ggplot(aes(x = k, y = value)) +
  geom_line() + geom_point() + facet_grid(rows = vars(cluster_eval), 
                                          scales = "free_y")

Мы продолжаем, запустив алгоритм PAM, передав аргумент k равным 3.

# find the best_k from stackoverflow using the mode function found on
# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode?page=1&tab=scoredesc#tab-top
mode = function(x){
  ux = unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}
best_k = mode(c(which.max(cluster_eval_df$silhouette), 
                which.max(cluster_eval_df$CH_index),
                which.max(cluster_eval_df$Dunn_index))) + 1
set.seed(100 + best_k - 1)  # for reproducibility
kmedoid = pam(dissimilarity_matrix, k = best_k, diss = TRUE, nstart = 10)   # set nstart as 10 for cluster stability
# Save the cluster label in the dataframe. Change it to factor to facilitate data wrangling.
cust_data$cluster_label = as.factor(kmedoid$clustering)

Статистические сводки функций для различных сегментов клиентов

Количество клиентов в каждом кластере.

#1: Number of observations (customers)
cust_data %>% group_by(cluster_label) %>% 
  summarise(n = n()) %>% kable(caption = "Number of observations in each cluster")

2. Среднее значение числовых признаков в каждом кластере.

#2: Average of numerical features for each cluster
summary_cont_features_per_cluster = 
  cust_data %>% group_by(cluster_label) %>% 
  summarise_if(is.numeric, mean) %>% select(-ID)
summary_cont_features_per_cluster

3. Распределение (доля) категориальных признаков для каждого кластера.

#3: Distribution of categorical features for each cluster
cate_features_names = names(Filter(is.factor, cust_data))
# Change the second argument of group_by() and third argument aes() to desired categorical feature names: Education, Marital_status, Kidhome, Teenhome.
cust_data %>% select(one_of(cate_features_names)) %>% 
  group_by(cluster_label, Teenhome) %>% summarise(n = n()) %>% 
  mutate(prop = (n/sum(n))) %>% 
  ggplot(aes(x = cluster_label, y = prop, fill = Teenhome)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_text(aes(label = round(prop, 2)), vjust = 1.5, size = 2,
            position = position_dodge(0.9)) +
  theme_minimal()

4. Распределение производных признаков (например, минимальное количество членов домохозяйства) по кластерам.

cust_data$min_num_household = factor(cust_data$min_num_household, ordered = TRUE)
cust_data %>% group_by(cluster_label, min_num_household) %>% summarise(n = n()) %>% 
  mutate(prop = (n/sum(n))) %>% 
  ggplot(aes(x = cluster_label, y = prop, fill = min_num_household)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_text(aes(label = round(prop, 2)), vjust = 1.5, size = 2,
            position = position_dodge(0.9)) +
  theme_minimal()

5. Занесите в таблицу бинарные (логические) переменные для каждого кластера.

table1 = 
  cust_data %>% 
  group_by(cluster_label, Complain) %>% summarise(n = n()) %>% 
  mutate(prop = (n/sum(n))) %>% select(-n)
colnames(table1)[2] = "Binary_outcomes"

#pattern = c("Complain", "Response", "Accept", "cluster")
#idx_select = lapply(pattern, function (x) grep(x, colnames(cust_data)))
#idx = unlist(idx_select)
cust_data_1 = cust_data %>% 
  select(contains("Response") | contains("Accept") | contains("cluster"))

n_col = ncol(cust_data_1)
for (i in (1:(n_col-2))){
  table2 = cust_data_1 %>%  select(c(n_col, all_of(i))) %>% 
    group_by_all() %>% 
    summarise(n = n()) %>% 
    mutate(prop = (n/sum(n))) %>% 
    select(-n)
  colnames(table2)[2] = "Binary_outcomes"
  table1 = table1 %>% 
    left_join(table2, by = c("cluster_label" = "cluster_label", 
                             "Binary_outcomes" = "Binary_outcomes"))
}
oldname = colnames(table1)[3:length(colnames(table1))]
newname = c("Complain", 
            colnames(cust_data_1)[-length(colnames(cust_data_1))])
# Rename
for (i in (1:length(newname))){
  names(table1)[names(table1) == oldname[i]] = newname[i]
}
table1 %>% as_tibble()

Характеристики каждого кластера, полученные из приведенных выше статистических сводок, сведены в следующую таблицу⁴:

Анализ главных компонентов (PCA) для смешанных данных

PCA реализован в первую очередь для уменьшения размерности, чтобы мы могли визуализировать данные в пространстве с меньшими размерностями.

res.famd = FAMD(cust_data_without_ID, graph = F)

eig.value = get_eigenvalue(res.famd)
head(eig.value)
fviz_screeplot(res.famd)

# predict to get the transformed feature space and plot on 3 dimensional scatter plot
transformed_data = predict(res.famd, cust_data_without_ID)
library(rgl)
scatter3d(x = transformed_data$coord[,1], y = transformed_data$coord[,2],
          z = transformed_data$coord[,3], 
          groups = cust_data$cluster_label, grid = FALSE, 
          surface = FALSE, ellipsoid = TRUE, 
          surface.col = c("#80FF00", "#009999", "#FF007F"))

В Rstudio появится всплывающее окно, если вы запустите приведенный выше фрагмент кода, и вы сможете свободно вращать интерактивную трехмерную точечную диаграмму и просматривать ее под разными углами. Если вы открываете коды уценки R в формате html, убедитесь, что в вашем браузере включен WebGL. Вы можете обратиться к этой статье о том, как включить WebGL в разных браузерах.

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

Рейтинг клиентов для определенного кластера (сегмента клиентов)

Согласно принципу Парето (также известному как правило 80–20), одна пятая часть (20%) клиентов приносит около 80% дохода компании. Таким образом, выявление наиболее ценных клиентов имеет решающее значение для компании, чтобы выделить больше ресурсов для развития тесных отношений с этими клиентами.

Техника предпочтения порядка по сходству (TOPSIS)

TOPSIS — это один из инструментов многокритериального анализа решений (MCDA), используемых для ранжирования или выбора альтернатив (клиентов) на основе нескольких критериев (выбранного набора характеристик). Интуитивно говоря, альтернативы, которые находятся на небольшом расстоянии от положительного идеального решения и на большом расстоянии от отрицательного идеального решения, должны получать высокие баллы. Вы можете найти подробное описание алгоритма TOPSIS на этой Вики-странице, так что здесь я бы не стал останавливаться на деталях. Суть внедрения TOPSIS заключается в следующем:

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

Предположим, что мы заинтересованы в ранжировании клиента в первом кластере,

# Calculate feature weights for continuous and categorical features for first cluster
i = 1
data_analysis_cont = cust_data %>% filter(cluster_label==i) %>% dplyr::select(starts_with(c("Mnt", "Num")))
data_analysis_cate = cust_data %>% filter(cluster_label==i) %>% 
  dplyr::select(starts_with(c("Accepted", "Response")))

h_cont = get.bw(scale(data_analysis_cont), bw = "nrd", nb = 'na')
w_cont = entropy.weight(scale(data_analysis_cont), h = h_cont)
w_cate = entropy.weight(data_analysis_cate, h='na')

Выполнение приведенного выше фрагмента кода может занять некоторое время.

data_analysis_cate <- lapply(data_analysis_cate, 
                             function(x) as.numeric(as.character(x)))
data_analysis = cbind(data_analysis_cont, data_analysis_cate)

feat_imp = c(w_cont, w_cate)
overall = TOPSIS(data_analysis, 
                 feat_imp, 
                 criteriaMinMax = rep("max", length(feat_imp)))
data_analysis$topsis_score = overall
data_analysis$ID = cust_data %>% filter(cluster_label==i) %>% select(ID)

# Top 10 customers for cluster 1
data_analysis %>% arrange(desc(topsis_score)) %>% as_tibble() %>% head(10) %>%  relocate(ID, topsis_score)

Заключение

Давайте кратко рассмотрим рабочий процесс этого эксперимента:

Получение данных. Загрузите данные и импортируйте их в рабочую область R.
Предварительная обработка данных.
- Удаление переменных с нулевой дисперсией.
- Отбрасывание выбросов / недействительных наблюдений.
- Извлечение признаков.
- Вменение пропущенных значений.
Визуализация данных.
Кластерный анализ методом k-medoid кластеризации.
Рейтинг клиентов TOPSIS.

У этого подхода к сегментации клиентов есть некоторые недостатки, которые я хотел бы выделить:

Это не инкрементное обучение. Если мы хотим сегментировать или сгруппировать нового клиента, нам нужно будет запустить модель кластеризации для всего фрагмента набора данных (включая новые данные).
Результаты могут быть не воспроизведены другими методами кластеризации. Стабильность и непротиворечивость кластера всегда были неотъемлемыми проблемами CA. Для этого эксперимента я использовал метод кластеризации на основе разделов, но нет гарантии, что это лучший метод. Лучшим подходом было бы попробовать другие методы кластеризации, которые исходят из других предположений, и сравнить результаты кластеризации.

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

Наконец, все коды в этой статье можно найти на Github и RPubs. Спасибо за прочтение!

[1]: Джеймс Г., Виттен Д., Хасти Т., Тибширани Р. (2021). Введение. В: Введение в статистическое обучение. Тексты Спрингера в статистике. Спрингер, Нью-Йорк, штат Нью-Йорк.

[2]: задача K-medoid является NP-сложной, а алгоритм PAM является эвристическим решением.

[3]: Помимо обычных метрик проверки, которые учитывают внутрикластерное и межкластерное сходство, стабильность кластера и сила прогнозирования являются некоторыми другими жизнеспособными вариантами.

[4]: Характеристики каждого кластера относятся к большинству членов кластера, а не ко всем.

Предложения по подаче заявок на Mlearning.ai
Как стать писателем на Mlearning.aimedium.com