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)

# 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
kmeans_data <- rfm[, c("recency_scaled", "frequency_scaled", "monetary_scaled")]

Error in exists(cacheKey, where = .rs.WorkingDataEnv, inherits = FALSE) : 
  invalid first argument
Error in assign(cacheKey, frame, .rs.CachedDataEnv) : 
  attempt to use zero-length variable name

# 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.9967277 

# 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: ")
8
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)

# 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 

# 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(kmeans_model$cluster, dist(kmeans_data))
avg_silhouette <- mean(sil[, "sil_width"])

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

Average Silhouette Width: 0.38318 

# 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] 8

$cluster.size
[1]   32  519    4 1066    3  352  635 1761

$min.cluster.size
[1] 3

$noisen
[1] 0

$diameter
[1] 10.586738  1.721251 14.771437  1.931527 12.673996 10.708891  3.480057  1.845847

$average.distance
[1]  3.3288854  0.5053873 10.7583198  0.4069616 10.5712937  1.1951336  0.5055392  0.4830964

$median.distance
[1]  2.7141943  0.4596258 12.4094830  0.3704346 11.8860221  0.9359738  0.4476375  0.4355743

$separation
[1]  0.21444031  0.00983200  9.25206163  0.01013293 10.39008523  0.07413154  0.00983200  0.01013293

$average.toother
[1]  6.687812  2.619226 25.705291  1.329533 28.374071  2.316040  1.685303  1.694311

$separation.matrix
           [,1]      [,2]      [,3]        [,4]     [,5]        [,6]        [,7]        [,8]
[1,]  0.0000000  4.089545 13.717897  3.38858823 10.66982  0.21444031  3.51212745  3.24143465
[2,]  4.0895453  0.000000 18.905945  1.14911109 22.85383  1.81823622  0.00983200  1.82964529
[3,] 13.7178974 18.905945  0.000000 17.98981473 10.39009  9.25206163 16.26086404 17.82189497
[4,]  3.3885882  1.149111 17.989815  0.00000000 22.88622  0.42066488  0.02173363  0.01013293
[5,] 10.6698195 22.853832 10.390085 22.88622046  0.00000 18.40516585 21.94562340 22.43432873
[6,]  0.2144403  1.818236  9.252062  0.42066488 18.40517  0.00000000  0.76451736  0.07413154
[7,]  3.5121275  0.009832 16.260864  0.02173363 21.94562  0.76451736  0.00000000  0.63126579
[8,]  3.2414346  1.829645 17.821895  0.01013293 22.43433  0.07413154  0.63126579  0.00000000

$ave.between.matrix
          [,1]      [,2]     [,3]       [,4]     [,5]      [,6]      [,7]       [,8]
[1,]  0.000000  7.565732 24.02495  6.7990596 22.58750  5.001686  7.045547  6.5033004
[2,]  7.565732  0.000000 25.99327  2.3146326 28.89569  3.648618  1.242046  2.9067023
[3,] 24.024954 25.993265  0.00000 25.8209380 22.54987 24.931342 25.848794 25.6892800
[4,]  6.799060  2.314633 25.82094  0.0000000 28.63445  2.113183  1.152452  0.7448827
[5,] 22.587505 28.895692 22.54987 28.6344502  0.00000 26.823880 28.698534 28.3739648
[6,]  5.001686  3.648618 24.93134  2.1131827 26.82388  0.000000  2.698112  1.7664069
[7,]  7.045547  1.242046 25.84879  1.1524516 28.69853  2.698112  0.000000  1.7377381
[8,]  6.503300  2.906702 25.68928  0.7448827 28.37396  1.766407  1.737738  0.0000000

$average.between
[1] 1.899229

$average.within
[1] 0.564919

$n.between
[1] 7039684

$n.within
[1] 2515322

$max.diameter
[1] 14.77144

$min.separation
[1] 0.009832

$within.cluster.ss
[1] 1639.373

$clus.avg.silwidths
        1         2         3         4         5         6         7         8 
0.2857795 0.5547946 0.4685329 0.3979470 0.4881997 0.2748271 0.4247718 0.3317208 

$avg.silwidth
[1] 0.38318

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.3247911

$dunn
[1] 0.0006656089

$dunn2
[1] 0.06923783

$entropy
[1] 1.493816

$wb.ratio
[1] 0.2974465

$ch
[1] 4363.246

$cwidegap
[1]  3.4331291  0.3625535 11.8906586  0.5992575 11.8860221  4.3314064  1.4109591  0.4994226

$widestgap
[1] 11.89066

$sindex
[1] 0.05512685

$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.0006656089 

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

$DB
[1] 0.8495511

$r
[1] 0.8518129 0.7093370 0.6272708 1.1689514 0.6272708 0.9334778 0.7093370 1.1689514

$R
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
[1,]       Inf 0.4550785 0.4217918 0.4895674 0.4136879 0.8518129 0.4883366 0.5217080
[2,] 0.4550785       Inf 0.2901002 0.3284309 0.2333918 0.4402072 0.7093370 0.2834922
[3,] 0.4217918 0.2901002       Inf 0.2882264 0.6272708 0.3309103 0.2918808 0.2921131
[4,] 0.4895674 0.3284309 0.2882264       Inf 0.2321807 0.7352239 0.6950362 1.1689514
[5,] 0.4136879 0.2333918 0.6272708 0.2321807       Inf 0.2782728 0.2351937 0.2365302
[6,] 0.8518129 0.4402072 0.3309103 0.7352239 0.2782728       Inf 0.6097411 0.9334778
[7,] 0.4883366 0.7093370 0.2918808 0.6950362 0.2351937 0.6097411       Inf 0.4901333
[8,] 0.5217080 0.2834922 0.2921131 1.1689514 0.2365302 0.9334778 0.4901333       Inf

$d
          1         2        3          4        5         6         7          8
1  0.000000  7.409534 23.66508  6.6947004 22.26535  4.784870  6.916591  6.4004317
2  7.409534  0.000000 25.71142  2.2977655 28.65593  3.527803  1.205033  2.8794675
3 23.665085 25.711420  0.00000 25.5510059 21.19958 24.667661 25.574112 25.4220845
4  6.694700  2.297766 25.55101  0.0000000 28.39877  1.983820  1.093989  0.6175568
5 22.265352 28.655926 21.19958 28.3987665  0.00000 26.563690 28.460641 28.1371885
6  4.784870  3.527803 24.66766  1.9838198 26.56369  0.000000  2.556281  1.6285349
7  6.916591  1.205033 25.57411  1.0939890 28.46064  2.556281  0.000000  1.6771179
8  6.400432  2.879468 25.42208  0.6175568 28.13719  1.628535  1.677118  0.0000000

$S
[1] 2.9473847 0.4245349 7.0343533 0.3301222 6.2635227 1.1284295 0.4302397 0.3917716

$centers
           [,1]        [,2]        [,3]
[1,] -0.8808462  6.08213030  1.78813906
[2,]  2.1450822 -0.38668175 -0.18655704
[3,] -0.2267805  2.80278496 25.21577909
[4,] -0.1500943 -0.28472428 -0.14788206
[5,] -0.9023141 22.33608281 17.00505334
[6,] -0.7642200  1.45089157  0.59105518
[7,]  0.9437050 -0.30419846 -0.14189884
[8,] -0.7108139 -0.04894461 -0.04121466

K-Means Redone

# 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)
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.9935658

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

Hopkins Statistic: 0.9935658 

# 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()

# 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

# 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)
library(factoextra)
library(cluster)
library(fpc)
library(flexclust)
library(mclust)
library(clusterSim)
library(hopkins)
install.packages("dbscan")

Error in install.packages : Updating loaded packages

library(dbscan)

# 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))
install.packages("dbscan")

Warning in install.packages :
  package ‘dbscan’ is in use and will not be installed

# 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 = 50)

# 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.7298603 

# 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]  131 4241

$min.cluster.size
[1] 131

$noisen
[1] 0

$diameter
[1] 35.126409  4.633622

$average.distance
[1] 5.390016 1.291642

$median.distance
[1] 2.805911 1.051016

$separation
[1] 0.1236574 0.1236574

$average.toother
[1] 5.284577 5.284577

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

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

$average.between
[1] 5.284577

$average.within
[1] 1.414443

$n.between
[1] 555571

$n.within
[1] 8999435

$max.diameter
[1] 35.12641

$min.separation
[1] 0.1236574

$within.cluster.ss
[1] 10586.09

$clus.avg.silwidths
         1          2 
-0.0759181  0.7547500 

$avg.silwidth
[1] 0.7298603

$g2
NULL

$g3
NULL

$pearsongamma
[1] 0.4874942

$dunn
[1] 0.003520354

$dunn2
[1] 0.9804381

$entropy
[1] 0.1346149

$wb.ratio
[1] 0.2676549

$ch
[1] 1043.123

$cwidegap
[1] 11.8906586  0.4802036

$widestgap
[1] 11.89066

$sindex
[1] 0.646103

$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.003520354 

# 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] 1.118255

$centers
           [,1]       [,2]        [,3]
[1,] 0.01878499 -0.1088257 -0.07521404