RFM Customer Segmentation

Customer Segmentation K-Means Clustering

# Load required libraries
library(ggplot2)
install.packages("dplyr")

Error in install.packages : Updating loaded packages

library(dplyr)
library(factoextra)
library(cluster)
library(fpc)
install.packages("flexclust")

Error in install.packages : Updating loaded packages

library(flexclust)
install.packages("mclust")

Error in install.packages : Updating loaded packages

library(mclust)
install.packages("clusterSim")

Error in install.packages : Updating loaded packages

library(cluster)
install.packages("hopkins")

Error in install.packages : Updating loaded packages

library(hopkins)
library(tidytext)
library(scatterplot3d)
library(tidyverse)

── Attaching core tidyverse packages ─────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ lubridate 1.9.2     ✔ tibble    3.2.1
✔ purrr     1.0.1     ✔ tidyr     1.3.0
✔ readr     2.1.4     ── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ purrr::map()    masks mclust::map()
✖ MASS::select()  masks dplyr::select()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

# Load the Online Retail dataset
data <- Online_Retail

# Data preprocessing
data <- data %>%
  dplyr::filter(!is.na(CustomerID)) %>%
  dplyr::select(CustomerID, InvoiceNo, InvoiceDate, UnitPrice, Description)

# Concatenate the transaction descriptions for each customer
rfm_with_descriptions <- data %>%
  group_by(CustomerID) %>%
  summarise(Descriptions = paste(Description, collapse = ", "))

# Calculate total spending (monetary value) per customer
monetary <- data %>%
  group_by(CustomerID) %>%
  summarise(monetary = sum(UnitPrice))

# Calculate the recency and frequency variables
# Recency is calculated by the time elapsed since the last day of the dataset
# Frequency is calculated by summing up the distinct invoices of a customer
recency <- data %>%
  group_by(CustomerID) %>%
  summarise(recency = as.numeric(difftime(max(data$InvoiceDate), max(InvoiceDate), units = "days")),
            frequency = n_distinct(InvoiceNo))

# Merge RFM variables with monetary value
rfm <- left_join(recency, monetary, by = "CustomerID")

# Data normalization
rfm$recency_scaled <- scale(rfm$recency)
rfm$frequency_scaled <- scale(rfm$frequency)
rfm$monetary_scaled <- scale(rfm$monetary)

# Prepare data frame to use for clustering
kmeans_data <- rfm[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]

# Calculate the Hopkins statistic
hopkins_stat <- hopkins::hopkins(X = as.matrix(kmeans_data), m = nrow(kmeans_data) - 1, method = "simple")
cat("Hopkins Statistic:", hopkins_stat, "\n")

Hopkins Statistic: 0.9965223 

# Prompt the user to choose the number of clusters based on the elbow plot
chosen_k <- readline(prompt = "Enter the optimal number of clusters based on the elbow plot: ")
4
chosen_k <- as.integer(chosen_k)

# Perform K-means clustering with chosen number of clusters
set.seed(123)
kmeans_model <- kmeans(kmeans_data, centers = chosen_k, nstart = 25)

# Add cluster labels to the original dataset
rfm$cluster <- as.factor(kmeans_model$cluster)

# Visualize the clusters
fviz_cluster(kmeans_model, data = kmeans_data)

# Create a 3D scatter plot
scatterplot3d(kmeans_data$recency_scaled, kmeans_data$frequency_scaled, kmeans_data$monetary_scaled,
              color = kmeans_model$cluster,
              main = "3D Scatter Plot of RFM Clusters",
              xlab = "Recency",
              ylab = "Frequency",
              zlab = "Monetary")
legend("topright", legend = levels(as.factor(kmeans_model$cluster)), col = 1:max(kmeans_model$cluster), pch = 1, cex = 1.2)

# Determine and visualise the optimal number of clusters
fviz_nbclust(kmeans_data, kmeans, method = "wss")

fviz_nbclust(kmeans_data, kmeans, method = "silhouette")

fviz_nbclust(kmeans_data, kmeans, method = "gap_stat")

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.......

Warning: did not converge in 10 iterations

............

Warning: did not converge in 10 iterations

..........

Warning: did not converge in 10 iterations

.........

Warning: Quick-TRANSfer stage steps exceeded maximum (= 218600)

.....

Warning: did not converge in 10 iterations

....... 50 
..................

Warning: did not converge in 10 iterationsWarning: did not converge in 10 iterations

....

Warning: did not converge in 10 iterations

............

Warning: did not converge in 10 iterations

................ 100 

View(cluster_word_counts)

# Visualize data points plotted against recency, monetary, and frequency
ggplot(rfm, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Perform PCA and extract results
pca <- prcomp(kmeans_data, scale. = TRUE)
pca_results <- factoextra::get_pca(pca)
plot(pca)

summary(pca)

Importance of components:
                          PC1    PC2    PC3
Standard deviation     1.2827 0.9555 0.6646
Proportion of Variance 0.5484 0.3044 0.1472
Cumulative Proportion  0.5484 0.8528 1.0000

biplot(pca, scale. = TRUE)

# Cluster analysis
cluster_analysis <- rfm %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis)

# Silhouette analysis
sil <- silhouette(kmeans_model$cluster, dist(kmeans_data))
avg_silhouette <- mean(sil[, "sil_width"])

cat("Average Silhouette Width:", avg_silhouette, "\n")

Average Silhouette Width: 0.6064892 

# Calculate clustering indices using cluster.stats
clustering_indices <- cluster.stats(dist(kmeans_data), kmeans_model$cluster)
print(clustering_indices)

$n
[1] 4372

$cluster.number
[1] 4

$cluster.size
[1]  227 1066    7 3072

$min.cluster.size
[1] 7

$noisen
[1] 0

$diameter
[1] 13.829273  2.512074 30.606993  3.874338

$average.distance
[1]  2.1388173  0.8135454 17.4696282  0.6914000

$median.distance
[1]  1.3177458  0.7170029 17.3298115  0.6229301

$separation
[1] 0.08155842 0.02178571 9.25206163 0.02178571

$average.toother
[1]  3.277427  2.310772 26.852580  2.400118

$separation.matrix
           [,1]        [,2]      [,3]        [,4]
[1,] 0.00000000  0.95490609  9.252062  0.08155842
[2,] 0.95490609  0.00000000 18.739599  0.02178571
[3,] 9.25206163 18.73959866  0.000000 16.26086404
[4,] 0.08155842  0.02178571 16.260864  0.00000000

$ave.between.matrix
          [,1]      [,2]     [,3]      [,4]
[1,]  0.000000  4.085878 25.15694  2.947035
[2,]  4.085878  0.000000 27.16396  2.122972
[3,] 25.156939 27.163963  0.00000 26.869825
[4,]  2.947035  2.122972 26.86983  0.000000

$average.between
[1] 2.548275

$average.within
[1] 0.8231975

$n.between
[1] 4244633

$n.within
[1] 5310373

$max.diameter
[1] 30.60699

$min.separation
[1] 0.02178571

$within.cluster.ss
[1] 3669.908

$clus.avg.silwidths
        1         2         3         4 
0.2301737 0.5796220 0.2833442 0.6443557 

$avg.silwidth
[1] 0.6064892

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.4766032

$dunn
[1] 0.0007117885

$dunn2
[1] 0.1215236

$entropy
[1] 0.7559626

$wb.ratio
[1] 0.323041

$ch
[1] 3746.454

$cwidegap
[1]  4.3314064  0.6296994 11.8906586  1.4109591

$widestgap
[1] 11.89066

$sindex
[1] 0.1784922

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list
dunn_index <- clustering_indices$dunn
cat("Dunn Index:", dunn_index, "\n")

Dunn Index: 0.0007117885 

# Calculate the Davies-Bouldin Index
db_index <- clusterSim::index.DB(kmeans_data, kmeans_model$cluster)
print(db_index)

$DB
[1] 0.8493978

$r
[1] 0.9988550 0.7495313 0.6503498 0.9988550

$R
          [,1]      [,2]      [,3]      [,4]
[1,]       Inf 0.7495313 0.6503498 0.9988550
[2,] 0.7495313       Inf 0.5297251 0.6043077
[3,] 0.6503498 0.5297251       Inf 0.5315909
[4,] 0.9988550 0.6043077 0.5315909       Inf

$d
          1         2        3         4
1  0.000000  3.852261 22.53791  2.770151
2  3.852261  0.000000 24.80879  2.070622
3 22.537912 24.808791  0.00000 24.495207
4  2.770151  2.070622 24.49521  0.000000

$S
[1]  2.2015383  0.6858520 12.4559884  0.5654409

$centers
           [,1]        [,2]        [,3]
[1,] -0.8005187  2.49251048  0.90014256
[2,]  1.5671071 -0.35021827 -0.17370844
[3,] -0.5162949 11.17419832 21.69689662
[4,] -0.4834649 -0.08811413 -0.05567625

# Merge the cluster labels with the transaction descriptions
rfm_with_descriptions <- left_join(rfm_with_descriptions, rfm, by = "CustomerID")

# Text Mining on the concatenated transaction descriptions of each cluster
cluster_texts <- rfm_with_descriptions %>%
  group_by(cluster) %>%
  summarise(Descriptions = paste(Descriptions, collapse = ", ")) %>%
  unnest_tokens(output = "word", input = Descriptions) %>%
  count(cluster, word, sort = TRUE)

# Get the top 10 words for each cluster
top_words <- cluster_texts %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = n)

# Rename the occurrence count column
top_words <- top_words %>% rename(frequency = n)

# Print the data frames for each cluster
cluster_word_counts <- split(top_words, top_words$cluster)
for (i in seq_along(cluster_word_counts)) {
  cat("Cluster", i, "Top Words:\n")
  print(cluster_word_counts[[i]])
  cat("\n")
}

Cluster 1 Top Words:

Cluster 2 Top Words:

Cluster 3 Top Words:

Cluster 4 Top Words:

NA

K-Means Redone after filtering out outliers

# Remove observations in cluster 3
rfm_filtered <- rfm[!(rfm$cluster %in% c(3, 5)), ]

# Prepare filtered data frame for clustering
kmeans_data_filtered <- rfm_filtered[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]

# Calculate the Hopkins statistic
hopkins_stat_filtered <- hopkins::hopkins(X = as.matrix(kmeans_data_filtered), m = nrow(kmeans_data_filtered) - 1, method = "simple")
cat("Hopkins Statistic (Filtered):", hopkins_stat_filtered, "\n")

Hopkins Statistic (Filtered): 0.9991269 

# Perform K-means clustering on filtered data
set.seed(123)
kmeans_model_filtered <- kmeans(kmeans_data_filtered, centers = chosen_k, nstart = 25)

# Add cluster labels to the original dataset
rfm_filtered$cluster <- as.factor(kmeans_model_filtered$cluster)

# Visualize the clusters after filtering
fviz_cluster(kmeans_model_filtered, data = kmeans_data_filtered)

# Determine and visualise the optimal number of clusters
fviz_nbclust(kmeans_data_filtered, kmeans, method = "wss")

fviz_nbclust(kmeans_data_filtered, kmeans, method = "silhouette")

fviz_nbclust(kmeans_data_filtered, kmeans, method = "gap_stat")

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
....

Warning: Quick-TRANSfer stage steps exceeded maximum (= 218250)

.............................................. 50 
............................

Warning: did not converge in 10 iterations

...................... 100 

# Visualize data points plotted against recency, monetary, and frequency after filtering
ggplot(rfm_filtered, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Cluster analysis on filtered data
cluster_analysis_filtered <- rfm_filtered %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis_filtered)

# Silhouette analysis on filtered data
sil_filtered <- silhouette(kmeans_model_filtered$cluster, dist(kmeans_data_filtered))
avg_silhouette_filtered <- mean(sil_filtered[, "sil_width"])

cat("Average Silhouette Width (filtered data):", avg_silhouette_filtered, "\n")

Average Silhouette Width (filtered data): 0.4118717 

# Calculate clustering indices using cluster.stats on filtered data
clustering_indices_filtered <- cluster.stats(dist(kmeans_data_filtered), kmeans_model_filtered$cluster)
print(clustering_indices_filtered)

$n
[1] 4365

$cluster.number
[1] 8

$cluster.size
[1]    3   10  505  513  955 1671  618   90

$min.cluster.size
[1] 3

$noisen
[1] 0

$diameter
[1] 4.637608 7.651552 1.721251 3.216320 1.756549 1.178992 3.486484 4.484979

$average.distance
[1] 3.1414414 4.0383326 0.4982096 0.7722308 0.4122750 0.3888178 0.4997264 1.5229407

$median.distance
[1] 4.3314064 4.1604236 0.4549806 0.6923219 0.3824275 0.3639188 0.4411893 1.3517279

$separation
[1] 3.132319092 1.108150582 0.016228707 0.022005240 0.008015243 0.008015243 0.016228707 0.113300538

$average.toother
[1] 7.661483 9.985118 2.583039 1.763605 1.253097 1.585100 1.659767 3.741790

$separation.matrix
         [,1]     [,2]       [,3]       [,4]        [,5]        [,6]       [,7]      [,8]
[1,] 0.000000 4.669482 6.26605446 3.15982900 5.101067938 5.390615733 3.72321134 3.1323191
[2,] 4.669482 0.000000 6.87209694 4.85726632 6.747912838 6.784611372 6.38069184 1.1081506
[3,] 6.266054 6.872097 0.00000000 1.50301034 1.141197353 1.929722331 0.01622871 2.4331014
[4,] 3.159829 4.857266 1.50301034 0.00000000 0.110595653 0.022005240 0.37930242 0.1133005
[5,] 5.101068 6.747913 1.14119735 0.11059565 0.000000000 0.008015243 0.01816411 1.8173503
[6,] 5.390616 6.784611 1.92972233 0.02200524 0.008015243 0.000000000 0.80226028 1.7291999
[7,] 3.723211 6.380692 0.01622871 0.37930242 0.018164110 0.802260278 0.00000000 1.9174183
[8,] 3.132319 1.108151 2.43310143 0.11330054 1.817350319 1.729199930 1.91741826 0.0000000

$ave.between.matrix
         [,1]      [,2]      [,3]     [,4]       [,5]       [,6]      [,7]     [,8]
[1,] 0.000000  9.242087  8.300834 7.030313  7.7010301  7.5913299  7.886982 6.830516
[2,] 9.242087  0.000000 10.723893 8.945040 10.1731496 10.0156326 10.353049 6.704747
[3,] 8.300834 10.723893  0.000000 3.253159  2.2656684  2.8804640  1.218974 4.880249
[4,] 7.030313  8.945040  3.253159 0.000000  1.4925997  1.2052657  2.219068 2.546709
[5,] 7.701030 10.173150  2.265668 1.492600  0.0000000  0.7188965  1.128866 3.771610
[6,] 7.591330 10.015633  2.880464 1.205266  0.7188965  0.0000000  1.726342 3.566322
[7,] 7.886982 10.353049  1.218974 2.219068  1.1288662  1.7263417  0.000000 4.168953
[8,] 6.830516  6.704747  4.880249 2.546709  3.7716103  3.5663222  4.168953 0.000000

$average.between
[1] 1.769639

$average.within
[1] 0.5010059

$n.between
[1] 7220316

$n.within
[1] 2304114

$max.diameter
[1] 7.651552

$min.separation
[1] 0.008015243

$within.cluster.ss
[1] 910.0416

$clus.avg.silwidths
        1         2         3         4         5         6         7         8 
0.5392068 0.3456493 0.5525951 0.3012711 0.3600885 0.4324474 0.4211287 0.3596871 

$avg.silwidth
[1] 0.4118717

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.4666251

$dunn
[1] 0.001047532

$dunn2
[1] 0.1780181

$entropy
[1] 1.576961

$wb.ratio
[1] 0.2831119

$ch
[1] 4745.983

$cwidegap
[1] 4.3314064 3.4331291 0.3625535 0.7670103 0.5992575 0.2119316 1.4109591 2.8238771

$widestgap
[1] 4.331406

$sindex
[1] 0.05342293

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list for filtered data
dunn_index_filtered <- clustering_indices_filtered$dunn
cat("Dunn Index (filtered data):", dunn_index_filtered, "\n")

Dunn Index (filtered data): 0.001047532 

# Calculate the Davies-Bouldin Index for filtered data
db_index_filtered <- clusterSim::index.DB(kmeans_data_filtered, kmeans_model_filtered$cluster)
print(db_index_filtered)

$DB
[1] 0.7945017

$r
[1] 0.5743004 0.6650391 0.7158673 0.8434492 1.0271850 1.0271850 0.7158673 0.7871202

$R
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
[1,]       Inf 0.5743004 0.3098843 0.3947608 0.3222320 0.3232245 0.3270802 0.5012482
[2,] 0.5743004       Inf 0.3229575 0.4131000 0.3317137 0.3343349 0.3353985 0.6650391
[3,] 0.3098843 0.3229575       Inf 0.3267631 0.3340353 0.2529788 0.7158673 0.3504934
[4,] 0.3947608 0.4131000 0.3267631       Inf 0.6875080 0.8434492 0.4938074 0.7871202
[5,] 0.3222320 0.3317137 0.3340353 0.6875080       Inf 1.0271850 0.7121625 0.4312429
[6,] 0.3232245 0.3343349 0.2529788 0.8434492 1.0271850       Inf 0.4359715 0.4483821
[7,] 0.3270802 0.3353985 0.7158673 0.4938074 0.7121625 0.4359715       Inf 0.4140183
[8,] 0.5012482 0.6650391 0.3504934 0.7871202 0.4312429 0.4483821 0.4140183       Inf

$d
         1         2         3        4          5         6         7        8
1 0.000000  8.899824  8.193446 6.959811  7.6128924 7.5079290  7.790317 6.739153
2 8.899824  0.000000 10.555013 8.756380 10.0174267 9.8600126 10.190438 6.387276
3 8.193446 10.555013  0.000000 3.198366  2.2475634 2.8634578  1.181372 4.780792
4 6.959811  8.756380  3.198366 0.000000  1.3951913 1.1059753  2.134730 2.393632
5 7.612892 10.017427  2.247563 1.395191  0.0000000 0.6215934  1.066896 3.686398
6 7.507929  9.860013  2.863458 1.105975  0.6215934 0.0000000  1.682293 3.486674
7 7.790317 10.190438  1.181372 2.134730  1.0668955 1.6822927  0.000000 4.069080
8 6.739153  6.387276  4.780792 2.393632  3.6863982 3.4866737  4.069080 0.000000

$S
[1] 2.1206864 2.9904864 0.4183342 0.6267738 0.3324313 0.3060601 0.4273717 1.2573018

$centers
            [,1]       [,2]        [,3]
[1,] -0.78621837  1.1341113  7.30636307
[2,] -0.87273658  9.0937740  3.32601727
[3,]  2.16148705 -0.3861514 -0.18579211
[4,] -0.73501559  0.8589995  0.35218043
[5,] -0.08269356 -0.2706824 -0.14264619
[6,] -0.69091896 -0.1539326 -0.08963987
[7,]  0.98359034 -0.3067994 -0.14216607
[8,] -0.82387529  3.1900600  0.88858769

Hierarchical Clustering

# Load required libraries
library(ggplot2)
install.packages("dplyr")

Error in install.packages : Updating loaded packages

library(dplyr)
library(factoextra)
library(cluster)
library(fpc)
library(flexclust)
library(mclust)
library(clusterSim)

Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:dplyr’:

    select

library(hopkins)

# Load the Online Retail dataset
data <- Online_Retail

# Data preprocessing
data <- data %>%
  dplyr::filter(!is.na(CustomerID)) %>%
  dplyr::select(CustomerID, InvoiceNo, InvoiceDate, UnitPrice)

# Calculate total spending (monetary value) per customer
monetary <- data %>%
  group_by(CustomerID) %>%
  summarise(monetary = sum(UnitPrice))

# Calculate the recency and frequency variables
# Recency is calculated by the time elapsed since the last day of the dataset
# Frequency is calculated by summing up the distinct invoices of a customer
recency <- data %>%
  group_by(CustomerID) %>%
  summarise(recency = as.numeric(difftime(max(data$InvoiceDate), max(InvoiceDate), units = "days")),
            frequency = n_distinct(InvoiceNo))

# Merge RFM variables with monetary value
rfm <- left_join(recency, monetary, by = "CustomerID")

# Calculate the Hopkins statistic
hopkins_stat <- hopkins::hopkins(X = as.matrix(rfm), m = nrow(rfm) - 1, method = "simple")
print(hopkins_stat)

[1] 0.9938167

# Check the value of the Hopkins statistic
cat("Hopkins Statistic:", hopkins_stat, "\n")

Hopkins Statistic: 0.9938167 

# Data normalization
rfm$recency_scaled <- scale(rfm$recency)
rfm$frequency_scaled <- scale(rfm$frequency)
rfm$monetary_scaled <- scale(rfm$monetary)

# Perform hierarchical clustering
hclust_data <- rfm[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]
hclust_model <- hclust(dist(hclust_data))

# Specify the desired number of clusters
desired_clusters <- 7

# Determine the cutoff height for the desired number of clusters
cutoff_height <- hclust_model$height[length(hclust_model$height) - (desired_clusters - 1)]
cat("Cut-off height at ", desired_clusters, " clusters:", cutoff_height, "\n")

Cut-off height at  7  clusters: 10.17666 

# Plot the dendogram with an indication of the cut-off
plot(hclust_model, labels = FALSE)
abline(h = cutoff_height, col = "red", lty = 2)

# Determine and visualise the optimal number of clusters up to 10 with bootstrapping up to 10
fviz_nbclust(hclust_data, hcut, method = "wss")
fviz_nbclust(hclust_data, hcut, method = "silhouette")

gap_stat_hclust <- clusGap(hclust_data, hcut, K.max = 10, B = 10)

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 10)  [one "." per sample]:
.......... 10 

fviz_gap_stat(gap_stat_hclust)

# Determine the tree cut for a desired number of clusters
cut <- cutree(hclust_model, k = desired_clusters)

# Add cluster labels to the original dataset
rfm$cluster <- as.factor(cut)

# Visualize data points plotted against recency, monetary, and frequency
ggplot(rfm, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Perform PCA and extract results
pca <- prcomp(hclust_data, scale. = TRUE)
pca_results <- factoextra::get_pca(pca)
plot(pca)
summary(pca)

Importance of components:
                          PC1    PC2    PC3
Standard deviation     1.2827 0.9555 0.6646
Proportion of Variance 0.5484 0.3044 0.1472
Cumulative Proportion  0.5484 0.8528 1.0000

biplot(pca, scale. = TRUE)

# Cluster analysis
cluster_analysis <- rfm %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis)

# Silhouette analysis
sil <- silhouette(cut, dist(hclust_data))
avg_silhouette <- mean(sil[, "sil_width"])

cat("Average Silhouette Width:", avg_silhouette, "\n")

Average Silhouette Width: 0.8036549 

# Calculate clustering indices using cluster.stats
clustering_indices <- cluster.stats(dist(hclust_data), cut)
print(clustering_indices)

$n
[1] 4372

$cluster.number
[1] 7

$cluster.size
[1] 4337   25    3    2    2    2    1

$min.cluster.size
[1] 1

$noisen
[1] 0

$diameter
[1]  6.125128 10.176661  4.637608  7.111027  7.153864  3.673587        NA

$average.distance
[1] 1.370173 3.601903 3.141441 7.111027 7.153864 3.673587      NaN

$median.distance
[1] 1.108700 3.180654 4.331406 7.111027 7.153864 3.673587       NA

$separation
[1]  0.6352564  0.6352564  3.1323191  9.2520616 10.3900852 11.8906586 11.8860221

$average.toother
[1] 11.233361  7.307411  7.682402 19.735676 24.859943 31.663681 35.385990

$separation.matrix
           [,1]       [,2]      [,3]      [,4]     [,5]     [,6]     [,7]
[1,]  0.0000000  0.6352564  3.132319 16.260864 19.29452 28.01929 31.17255
[2,]  0.6352564  0.0000000  3.897518 13.717897 10.66982 27.11160 22.17650
[3,]  3.1323191  3.8975177  0.000000  9.252062 18.40517 21.10512 29.48117
[4,] 16.2608640 13.7178974  9.252062  0.000000 10.39009 11.89066 19.19929
[5,] 19.2945174 10.6698195 18.405166 10.390085  0.00000 21.84576 11.88602
[6,] 28.0192930 27.1116034 21.105120 11.890659 21.84576  0.00000 24.28120
[7,] 31.1725465 22.1764961 29.481169 19.199295 11.88602 24.28120  0.00000

$ave.between.matrix
          [,1]      [,2]      [,3]     [,4]     [,5]     [,6]     [,7]
[1,]  0.000000  7.281518  7.659358 19.75646 24.90601 31.69022 35.44623
[2,]  7.281518  0.000000  8.030149 17.57513 18.37130 30.07522 29.23605
[3,]  7.659358  8.030149  0.000000 12.54627 20.65035 24.31985 29.99149
[4,] 19.756461 17.575126 12.546266  0.00000 17.33963 13.44133 22.55345
[5,] 24.906008 18.371296 20.650353 17.33963  0.00000 26.20010 12.28001
[6,] 31.690222 30.075221 24.319853 13.44133 26.20010  0.00000 25.66630
[7,] 35.446227 29.236051 29.991487 22.55345 12.28001 25.66630  0.00000

$average.between
[1] 11.24764

$average.within
[1] 1.39048

$n.between
[1] 152084

$n.within
[1] 9402922

$max.diameter
[1] 10.17666

$min.separation
[1] 0.6352564

$within.cluster.ss
[1] 6490.281

$clus.avg.silwidths
        1         2         3         4         5         6         7 
0.8063232 0.4646934 0.5955956 0.4322680 0.4168380 0.7266762 0.0000000 

$avg.silwidth
[1] 0.8036549

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.6455611

$dunn
[1] 0.06242287

$dunn2
[1] 1.017844

$entropy
[1] 0.05497188

$wb.ratio
[1] 0.1236241

$ch
[1] 742.345

$cwidegap
[1] 2.433101 3.433129 4.331406 7.111027 7.153864 3.673587 0.000000

$widestgap
[1] 7.153864

$sindex
[1] 2.791934

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list
dunn_index <- clustering_indices$dunn
cat("Dunn Index:", dunn_index, "\n")

Dunn Index: 0.06242287 

# Calculate the Davies-Bouldin Index
db_index <- clusterSim::index.DB(hclust_data, cut)
print(db_index)

$DB
[1] 0.5097034

$r
[1] 0.5979024 0.6801343 0.6801343 0.4743650 0.4124213 0.4186534 0.3043130

$R
           [,1]      [,2]       [,3]      [,4]      [,5]       [,6]       [,7]
[1,]        Inf 0.5979024 0.44223797 0.2447884 0.1934706 0.09586504 0.03371190
[2,] 0.59790238       Inf 0.68013435 0.3876336 0.3681847 0.16324043 0.10425155
[3,] 0.44223797 0.6801343        Inf 0.4743650 0.2800601 0.16327935 0.07088147
[4,] 0.24478838 0.3876336 0.47436500       Inf 0.4124213 0.41865345 0.15786441
[5,] 0.19347064 0.3681847 0.28006012 0.4124213       Inf 0.20708594 0.30431295
[6,] 0.09586504 0.1632404 0.16327935 0.4186534 0.2070859        Inf 0.07164359
[7,] 0.03371190 0.1042516 0.07088147 0.1578644 0.3043130 0.07164359        NaN

$d
          1         2         3        4        5        6        7
1  0.000000  7.087909  7.496264 19.40435 24.66203 31.61987 35.43100
2  7.087909  0.000000  7.592790 17.02367 17.98109 29.89593 29.19315
3  7.496264  7.592790  0.000000 11.96589 20.34427 24.23748 29.91877
4 19.404353 17.023667 11.965892  0.00000 17.29408 12.88012 22.52258
5 24.662028 17.981090 20.344268 17.29408  0.00000 26.14241 11.75412
6 31.619867 29.895930 24.237480 12.88012 26.14241  0.00000 25.63793
7 35.431003 29.193149 29.918770 22.52258 11.75412 25.63793  0.00000

$S
[1] 1.194446 3.043431 2.120686 3.555514 3.576932 1.836794 0.000000

$centers
             [,1]        [,2]        [,3]
[1,]  0.006427956 -0.05711689 -0.05207893
[2,] -0.876213115  6.64376839  2.08275854
[3,] -0.786218366  1.13411127  7.30636307
[4,] -0.466269341  4.06098276 18.90436356
[5,] -0.903461530 20.49786474 13.54499588
[6,]  0.012708362  1.54458715 31.52719462
[7,] -0.900019316 26.01251895 23.92516824

Hierarchical Clustering Redone after filtering out outliers

# Filter out observations in clusters other than cluster 1
rfm_filtered <- rfm[rfm$cluster == 1, ]

# Perform hierarchical clustering on the filtered dataset
hclust_data_filtered <- rfm_filtered[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]
hclust_model_filtered <- hclust(dist(hclust_data_filtered))

# Specify the desired number of clusters
desired_clusters_filtered <- 7

# Determine the cutoff height for the desired number of clusters
cutoff_height_filtered <- hclust_model_filtered$height[length(hclust_model_filtered$height) - (desired_clusters_filtered - 1)]
cat("Cut-off height at", desired_clusters_filtered, "cluster(s):", cutoff_height_filtered, "\n")

Cut-off height at 7 cluster(s): 2.968331 

# Plot the dendrogram with an indication of the cut-off
plot(hclust_model_filtered, labels = FALSE)
abline(h = cutoff_height_filtered, col = "red", lty = 2)

# Determine and visualize the optimal number of clusters up to 10 with bootstrapping up to 10
fviz_nbclust(hclust_data_filtered, hcut, method = "wss")

fviz_nbclust(hclust_data_filtered, hcut, method = "silhouette")

gap_stat_hclust_filtered <- clusGap(hclust_data_filtered, hcut, K.max = 10, B = 10)

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 10)  [one "." per sample]:
.......... 10 

fviz_gap_stat(gap_stat_hclust_filtered)

# Determine the tree cut for a desired number of clusters
cut_filtered <- cutree(hclust_model_filtered, k = desired_clusters_filtered)

# Add cluster labels to the filtered dataset
rfm_filtered$cluster <- as.factor(cut_filtered)

# Visualize data points plotted against recency, monetary, and frequency
ggplot(rfm_filtered, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Cluster analysis
cluster_analysis_filtered <- rfm_filtered %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis_filtered)

# Silhouette analysis
sil_filtered <- silhouette(cut_filtered, dist(hclust_data_filtered))
avg_silhouette_filtered <- mean(sil_filtered[, "sil_width"])

cat("Average Silhouette Width:", avg_silhouette_filtered, "\n")

Average Silhouette Width: 0.5357976 

# Calculate clustering indices using cluster.stats
clustering_indices_filtered <- cluster.stats(dist(hclust_data_filtered), cut_filtered)
print(clustering_indices_filtered)

$n
[1] 4337

$cluster.number
[1] 7

$cluster.size
[1]  697 3437   16  167   15    4    1

$min.cluster.size
[1] 1

$noisen
[1] 0

$diameter
[1] 1.6371974 2.9683314 2.4724298 2.7461669 1.3157118 0.8605638        NA

$average.distance
[1] 0.5996526 0.8262783 1.1791495 0.9765076 0.5446674 0.6090913       NaN

$median.distance
[1] 0.5274782 0.7246795 1.2816614 0.9385098 0.5243770 0.5926533        NA

$separation
[1] 0.02074356 0.02074356 0.31379656 0.10850601 0.32610236 0.61479999 2.43310143

$average.toother
[1] 2.422730 2.357982 2.371085 2.480977 4.175850 3.999459 4.100689

$separation.matrix
           [,1]       [,2]      [,3]      [,4]      [,5]     [,6]     [,7]
[1,] 0.00000000 0.02074356 1.5151881 1.9544983 3.9808810 4.166508 2.937616
[2,] 0.02074356 0.00000000 0.3546532 0.1085060 2.0362918 1.943674 2.433101
[3,] 1.51518813 0.35465321 0.0000000 0.3137966 2.9035891 1.441290 3.429152
[4,] 1.95449833 0.10850601 0.3137966 0.0000000 0.3261024 0.614800 2.821362
[5,] 3.98088100 2.03629178 2.9035891 0.3261024 0.0000000 1.632430 2.840101
[6,] 4.16650807 1.94367383 1.4412905 0.6148000 1.6324298 0.000000 3.341107
[7,] 2.93761609 2.43310143 3.4291521 2.8213618 2.8401015 3.341107 0.000000

$ave.between.matrix
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
[1,] 0.000000 2.341808 3.284056 3.690433 5.158477 4.996068 3.732235
[2,] 2.341808 0.000000 2.188476 2.238880 4.081119 3.899362 4.220792
[3,] 3.284056 2.188476 0.000000 2.171796 3.837775 2.600005 4.019626
[4,] 3.690433 2.238880 2.171796 0.000000 2.108182 2.189105 3.282149
[5,] 5.158477 4.081119 3.837775 2.108182 0.000000 2.306935 3.060305
[6,] 4.996068 3.899362 2.600005 2.189105 2.306935 0.000000 3.504820
[7,] 3.732235 4.220792 4.019626 3.282149 3.060305 3.504820 0.000000

$average.between
[1] 2.420411

$average.within
[1] 0.7957624

$n.between
[1] 3241202

$n.within
[1] 6161414

$max.diameter
[1] 2.968331

$min.separation
[1] 0.02074356

$within.cluster.ss
[1] 1895.045

$clus.avg.silwidths
        1         2         3         4         5         6         7 
0.7351321 0.5047851 0.3743677 0.3391192 0.7295284 0.7138247 0.0000000 

$avg.silwidth
[1] 0.5357976

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.7707657

$dunn
[1] 0.006988289

$dunn2
[1] 1.787883

$entropy
[1] 0.6521774

$wb.ratio
[1] 0.3287716

$ch
[1] 1634.683

$cwidegap
[1] 0.3619369 0.4913610 1.4109591 0.4526189 0.3599992 0.5721821 0.0000000

$widestgap
[1] 1.410959

$sindex
[1] 0.1843331

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list
dunn_index_filtered <- clustering_indices_filtered$dunn
cat("Dunn Index:", dunn_index_filtered, "\n")

Dunn Index: 0.006988289 

# Calculate the Davies-Bouldin Index
db_index_filtered <- clusterSim::index.DB(hclust_data_filtered, cut_filtered)
print(db_index_filtered)

$DB
[1] 0.6298347

$r
[1] 0.5154591 0.8014672 0.8500708 0.8500708 0.5918395 0.5583761 0.2415594

$R
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,]       Inf 0.5154591 0.4419558 0.3528346 0.1795997 0.1804210 0.1365571
[2,] 0.5154591       Inf 0.8014672 0.6893572 0.2724894 0.2789222 0.1633706
[3,] 0.4419558 0.8014672       Inf 0.8500708 0.3516702 0.5206531 0.2293738
[4,] 0.3528346 0.6893572 0.8500708       Inf 0.5918395 0.5583761 0.2415594
[5,] 0.1795997 0.2724894 0.3516702 0.5918395       Inf 0.3575075 0.1373083
[6,] 0.1804210 0.2789222 0.5206531 0.5583761 0.3575075       Inf 0.1123374
[7,] 0.1365571 0.1633706 0.2293738 0.2415594 0.1373083 0.1123374       NaN

$d
         1        2        3        4        5        6        7
1 0.000000 2.304298 3.181986 3.621064 5.131752 4.970379 3.699563
2 2.304298 0.000000 1.975958 2.110669 4.033292 3.850998 4.178045
3 3.181986 1.975958 0.000000 1.968696 3.746566 2.482756 3.928504
4 3.621064 2.110669 1.968696 0.000000 2.008813 2.084606 3.197703
5 5.131752 4.033292 3.746566 2.008813 0.000000 2.260145 3.033027
6 4.970379 3.850998 2.482756 2.084606 2.260145 0.000000 3.485562
7 3.699563 4.178045 3.928504 3.197703 3.033027 3.485562 0.000000

$S
[1] 0.5052015 0.6825697 0.9010957 0.7724352 0.4164596 0.3915591 0.0000000

$centers
           [,1]       [,2]        [,3]
[1,]  1.9435668 -0.3838634 -0.18685373
[2,] -0.3410164 -0.1062839 -0.07101754
[3,] -0.3849143  0.3733924  1.84533135
[4,] -0.8030149  1.8405017  0.60092543
[5,] -0.8715853  3.8468214  0.52806302
[6,] -0.8665533  2.7224745  2.48869383
[7,]  2.0874779  3.2043374  0.70245458

DBSCAN

# Load required libraries
library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(factoextra)

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(cluster)
library(fpc)
library(flexclust)

Loading required package: grid
Loading required package: lattice
Loading required package: modeltools
Loading required package: stats4

library(mclust)

                   __           __ 
   ____ ___  _____/ /_  _______/ /_
  / __ `__ \/ ___/ / / / / ___/ __/
 / / / / / / /__/ / /_/ (__  ) /_  
/_/ /_/ /_/\___/_/\__,_/____/\__/   version 6.0.0
Type 'citation("mclust")' for citing this R package in publications.

library(clusterSim)

Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:dplyr’:

    select

library(hopkins)
install.packages("dbscan")

Installing package into ‘C:/Users/Frederick/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/dbscan_1.1-11.zip'
Content type 'application/zip' length 2837428 bytes (2.7 MB)
downloaded 2.7 MB

package ‘dbscan’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\Frederick\AppData\Local\Temp\Rtmpklsq3t\downloaded_packages

library(dbscan)

Attaching package: ‘dbscan’

The following object is masked from ‘package:fpc’:

    dbscan

The following object is masked from ‘package:stats’:

    as.dendrogram

# Load the Online Retail dataset
data <- Online_Retail

# Data preprocessing
data <- data %>%
  dplyr::filter(!is.na(CustomerID)) %>%
  dplyr::select(CustomerID, InvoiceNo, InvoiceDate, UnitPrice)

# Calculate total spending (monetary value) per customer
monetary <- data %>%
  group_by(CustomerID) %>%
  summarise(monetary = sum(UnitPrice))

# Calculate the recency and frequency variables
# Recency is calculated by the time elapsed since the last day of the dataset
# Frequency is calculated by summing up the distinct invoices of a customer
recency <- data %>%
  group_by(CustomerID) %>%
  summarise(recency = as.numeric(difftime(max(data$InvoiceDate), max(InvoiceDate), units = "days")),
            frequency = n_distinct(InvoiceNo))

# Merge RFM variables with monetary value
rfm <- left_join(recency, monetary, by = "CustomerID")

# Data normalization
rfm$recency_scaled <- scale(rfm$recency)
rfm$frequency_scaled <- scale(rfm$frequency)
rfm$monetary_scaled <- scale(rfm$monetary)

# Prepare data frame to use for clustering
dbscan_data <- rfm[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]

# Perform DBSCAN clustering
dbscan_model <- dbscan(dbscan_data, eps = 0.5, minPts = 700)

# Add cluster labels to the original dataset
rfm$cluster <- as.factor(dbscan_model$cluster)

# Visualize the clusters
fviz_cluster(dbscan_model, data = dbscan_data)

# Visualize data points plotted against recency, monetary, and frequency
ggplot(rfm, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Cluster analysis
cluster_analysis <- rfm %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis)

# Silhouette analysis
sil <- silhouette(dbscan_model$cluster, dist(dbscan_data))
avg_silhouette <- mean(sil[, "sil_width"])

cat("Average Silhouette Width:", avg_silhouette, "\n")

Average Silhouette Width: 0.5010976 

# Calculate clustering indices using cluster.stats
clustering_indices <- cluster.stats(dist(dbscan_data), dbscan_model$cluster)

Warning: clustering renumbered because maximum != number of clusters

print(clustering_indices)

$n
[1] 4372

$cluster.number
[1] 2

$cluster.size
[1] 1345 3027

$min.cluster.size
[1] 1345

$noisen
[1] 0

$diameter
[1] 36.021083  2.089346

$average.distance
[1] 2.2562992 0.6607006

$median.distance
[1] 1.4529015 0.6026901

$separation
[1] 0.01714471 0.01714471

$average.toother
[1] 2.340683 2.340683

$separation.matrix
           [,1]       [,2]
[1,] 0.00000000 0.01714471
[2,] 0.01714471 0.00000000

$ave.between.matrix
         [,1]     [,2]
[1,] 0.000000 2.340683
[2,] 2.340683 0.000000

$average.between
[1] 2.340683

$average.within
[1] 1.15157

$n.between
[1] 4071315

$n.within
[1] 5483691

$max.diameter
[1] 36.02108

$min.separation
[1] 0.01714471

$within.cluster.ss
[1] 10954.77

$clus.avg.silwidths
         1          2 
0.03345671 0.70888653 

$avg.silwidth
[1] 0.5010976

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.3659254

$dunn
[1] 0.0004759633

$dunn2
[1] 1.037399

$entropy
[1] 0.617199

$wb.ratio
[1] 0.4919802

$ch
[1] 860.9468

$cwidegap
[1] 11.8906586  0.2174164

$widestgap
[1] 11.89066

$sindex
[1] 0.1453776

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list
dunn_index <- clustering_indices$dunn
cat("Dunn Index:", dunn_index, "\n")

Dunn Index: 0.0004759633 

# Calculate the Davies-Bouldin Index
db_index <- clusterSim::index.DB(dbscan_data, dbscan_model$cluster)

Warning: no non-missing arguments to max; returning -Inf

print(db_index)

$DB
[1] NaN

$r
[1] -Inf

$R
     [,1]
[1,]  Inf

$d
  1
1 0

$S
[1] 0.534274

$centers
           [,1]       [,2]       [,3]
[1,] -0.4375446 -0.1420129 -0.0879324

DBSCAN Redone after filtering out outliers

# Refilter the dataset to remove specific entries that are outliers
rfm_filtered <- rfm[-c(329, 2723, 4043, 331, 1301, 2028, 1301, 1896), ]

# Prepare data frame to use for clustering
dbscan_data_filtered <- rfm_filtered[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]

# Perform DBSCAN clustering
dbscan_model_filtered <- dbscan(dbscan_data_filtered, eps = 0.5, minPts = 50)

# Add cluster labels to the original dataset
rfm_filtered$cluster <- as.factor(dbscan_model_filtered$cluster)

# Visualize the clusters
fviz_cluster(dbscan_model_filtered, data = dbscan_data_filtered)

# Visualize data points plotted against recency, monetary, and frequency
ggplot(rfm_filtered, aes(x = recency, y = monetary, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Monetary", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = recency, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Recency", y = "Frequency", color = "Cluster") +
  theme_minimal()

ggplot(rfm_filtered, aes(x = monetary, y = frequency, color = cluster)) +
  geom_point() +
  labs(x = "Monetary", y = "Frequency", color = "Cluster") +
  theme_minimal()

# Cluster analysis
cluster_analysis_filtered <- rfm_filtered %>%
  group_by(cluster) %>%
  summarise(average_recency = mean(recency),
            average_frequency = mean(frequency),
            average_monetary = mean(monetary),
            count_customers = n())

print(cluster_analysis_filtered)

# Silhouette analysis
sil_filtered <- silhouette(dbscan_model_filtered$cluster, dist(dbscan_data_filtered))
avg_silhouette_filtered <- mean(sil[, "sil_width"])
cat("Average Silhouette Width (Filtered):", avg_silhouette, "\n")

Average Silhouette Width (Filtered): 0.5010976 

# Calculate clustering indices using cluster.stats
clustering_indices_filtered <- cluster.stats(dist(dbscan_data_filtered), dbscan_model_filtered$cluster)

Warning: clustering renumbered because maximum != number of clusters

print(clustering_indices_filtered)

$n
[1] 4365

$cluster.number
[1] 2

$cluster.size
[1]  124 4241

$min.cluster.size
[1] 124

$noisen
[1] 0

$diameter
[1] 14.566632  4.633622

$average.distance
[1] 3.166538 1.291642

$median.distance
[1] 2.510941 1.051016

$separation
[1] 0.1236574 0.1236574

$average.toother
[1] 4.063365 4.063365

$separation.matrix
          [,1]      [,2]
[1,] 0.0000000 0.1236574
[2,] 0.1236574 0.0000000

$ave.between.matrix
         [,1]     [,2]
[1,] 0.000000 4.063365
[2,] 4.063365 0.000000

$average.between
[1] 4.063365

$average.within
[1] 1.344903

$n.between
[1] 525884

$n.within
[1] 8998546

$max.diameter
[1] 14.56663

$min.separation
[1] 0.1236574

$within.cluster.ss
[1] 6323.3

$clus.avg.silwidths
        1         2 
0.1813089 0.6783562 

$avg.silwidth
[1] 0.6642361

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.5156053

$dunn
[1] 0.008489087

$dunn2
[1] 1.28322

$entropy
[1] 0.1291632

$wb.ratio
[1] 0.3309827

$ch
[1] 1052.754

$cwidegap
[1] 4.3314064 0.4802036

$widestgap
[1] 4.331406

$sindex
[1] 0.6456025

$corrected.rand
NULL

$vi
NULL

# Extract the Dunn index from the clustering indices list
dunn_index_filtered <- clustering_indices_filtered$dunn
cat("Dunn Index (Filtered):", dunn_index, "\n")

Dunn Index (Filtered): 0.0004759633 

# Calculate the Davies-Bouldin Index
db_index_filtered <- clusterSim::index.DB(dbscan_data_filtered, dbscan_model_filtered$cluster)

Warning: no non-missing arguments to max; returning -Inf

print(db_index)

$DB
[1] NaN

$r
[1] -Inf

$R
     [,1]
[1,]  Inf

$d
  1
1 0

$S
[1] 0.534274

$centers
           [,1]       [,2]       [,3]
[1,] -0.4375446 -0.1420129 -0.0879324