Segmentasi pelanggan toko daring berdasarkan recency, frequency, dan monetary

# Contoh kode R: K-Means and Hierarchical Cluster Analysis
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(cluster)
  library(ggdendro)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: objek dari 3 atau 4 kelompok laten, 4 variabel karakteristik ----
n <- 200
positif <- function(x, m) if (m > 0) pmax(x, 0.05 * m) else x   # cegah nilai negatif yang tidak wajar
k_benar <- sample(3:4, 1)
repeat {                                       # pusat kelompok yang cukup terpisah (skala z)
  pusat <- matrix(rnorm(k_benar * 4, 0, 1.6), k_benar, 4)
  if (min(dist(pusat)) > 3.2) break
}
anggota <- sample(seq_len(k_benar), n, replace = TRUE, prob = runif(k_benar, 0.6, 1.4))
sebaran <- runif(1, 0.5, 0.72)
z <- pusat[anggota, ] + matrix(rnorm(n * 4, 0, sebaran), n, 4)
d <- data.frame(purchase_recency = positif(100 + 0.6 * 15 * z[, 1], 100),
                purchase_frequency = positif(50 + 0.6 * 10 * z[, 2], 50),
                monetary_value = positif(30 + 0.6 * 6 * z[, 3], 30),
                x3 = round(50 + 6 * z[, 4], 2))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Standarisasi variabel ----
vars <- c("purchase_recency", "purchase_frequency", "monetary_value", "x3")
xs <- scale(d[vars])
jarak <- dist(xs)

# ---- Penentuan jumlah klaster: elbow (WSS) dan silhouette rata-rata, k = 2..8 ----
daftar_k <- 2:8
uji_k <- lapply(daftar_k, function(k) {
  km <- kmeans(xs, centers = k, nstart = 25, iter.max = 50)
  sil <- cluster::silhouette(km$cluster, jarak)
  c(k = k, wss = km$tot.withinss, sil = mean(sil[, "sil_width"]))
})
uji_k <- as.data.frame(do.call(rbind, uji_k))
print(round(uji_k, 3))
k_best <- uji_k$k[which.max(uji_k$sil)]
cat("Jumlah klaster terpilih (silhouette maksimum):", k_best, "\n")

# ---- K-Means akhir (nstart 25); label klaster diurutkan menurut rerata Y ----
km <- kmeans(xs, centers = k_best, nstart = 25, iter.max = 50)
urut <- order(tapply(d$purchase_recency, km$cluster, mean))
d$klaster <- factor(paste0("C", match(km$cluster, urut)), levels = paste0("C", seq_len(k_best)))
sil_akhir <- mean(cluster::silhouette(km$cluster, jarak)[, "sil_width"])
rasio_bss <- km$betweenss / km$totss

# ---- Klaster hierarki Ward (ward.D2) dan kesesuaiannya dengan K-Means ----
hc <- hclust(jarak, method = "ward.D2")
ward <- cutree(hc, k = k_best)
print(table(KMeans = d$klaster, Ward = ward))
ari <- function(a, b) {
  # adjusted Rand index (Hubert dan Arabie)
  tb <- table(a, b); c2 <- function(x) sum(x * (x - 1) / 2)
  s_ij <- c2(tb); s_a <- c2(rowSums(tb)); s_b <- c2(colSums(tb)); tot <- c2(sum(tb))
  harap <- s_a * s_b / tot
  (s_ij - harap) / (0.5 * (s_a + s_b) - harap)
}
ari_km_ward <- ari(d$klaster, ward)
cat(sprintf("Adjusted Rand index K-Means vs Ward = %.3f\n", ari_km_ward))

# ---- Profil klaster (rerata satuan asli) ----
profil <- aggregate(d[vars], by = list(klaster = d$klaster), FUN = mean)
profil$n <- as.numeric(table(d$klaster))
print(profil, digits = 4)
k_atas <- as.character(profil$klaster[which.max(profil$purchase_recency)])


# ---- Grafik ----
dk <- rbind(data.frame(k = uji_k$k, nilai = uji_k$sil, ukuran = "Mean silhouette"),
            data.frame(k = uji_k$k, nilai = uji_k$wss, ukuran = "Within SS (elbow)"))
p1 <- ggplot(dk, aes(k, nilai)) +
  geom_vline(xintercept = k_best, linetype = "dashed", colour = warna$redup) +
  geom_line(colour = pal[1], linewidth = 0.9) +
  geom_point(aes(colour = k == k_best), size = 2.2) +
  scale_colour_manual(values = c(`FALSE` = pal[1], `TRUE` = pal[2]), guide = "none") +
  facet_wrap(~ukuran, ncol = 1, scales = "free_y") +
  scale_x_continuous(breaks = daftar_k) +
  labs(title = "Choosing k", x = "Number of clusters k", y = NULL)
pc <- prcomp(xs)
ve <- 100 * pc$sdev^2 / sum(pc$sdev^2)
sk <- data.frame(PC1 = pc$x[, 1], PC2 = pc$x[, 2], klaster = d$klaster)
pusat_pc <- aggregate(cbind(PC1, PC2) ~ klaster, data = sk, FUN = mean)
p2 <- ggplot(sk, aes(PC1, PC2, colour = klaster)) +
  stat_ellipse(aes(fill = klaster), geom = "polygon", alpha = 0.12, level = 0.9, show.legend = FALSE) +
  geom_point(size = 1.5, alpha = 0.75) +
  geom_point(data = pusat_pc, shape = 4, size = 4, stroke = 1.4, colour = warna$teks) +
  geom_text(data = pusat_pc, aes(label = klaster), colour = warna$teks, fontface = "bold", size = 3.6, vjust = -1.1) +
  scale_colour_manual(values = pal, name = "Cluster") +
  scale_fill_manual(values = pal) +
  labs(title = "K-Means clusters on the first two PCs", x = sprintf("PC1 (%.1f%%)", ve[1]), y = sprintf("PC2 (%.1f%%)", ve[2]))
zp <- aggregate(as.data.frame(xs), by = list(klaster = d$klaster), FUN = mean)
nama_var <- sapply(c("purchase_recency", "purchase_frequency", "monetary_value", "x3 (extra)"), function(s) paste(strwrap(s, 14), collapse = "
"))
pl <- data.frame(klaster = rep(zp$klaster, length(vars)), variabel = rep(nama_var, each = nrow(zp)), z = unlist(zp[vars]))
pl$variabel <- factor(pl$variabel, levels = unique(nama_var))
p3 <- ggplot(pl, aes(variabel, klaster, fill = z)) +
  geom_tile(colour = warna$panel, linewidth = 0.8) +
  geom_text(aes(label = sprintf("%.2f", round(z, 2) + 0)), colour = warna$teks, size = 3.1) +
  scale_fill_gradient2(low = pal[1], mid = warna$panel, high = pal[2], midpoint = 0, name = "z") +
  labs(title = "Cluster profile (mean z-score)", x = NULL, y = NULL) +
  theme(panel.grid = element_blank(), legend.position = "none", axis.text.x = element_text(size = rel(0.8)))
dd <- ggdendro::dendro_data(as.dendrogram(hc), type = "rectangle")
h_tinggi <- sort(hc$height, decreasing = TRUE)
h_potong <- mean(h_tinggi[(k_best - 1):k_best])
p4 <- ggplot(ggdendro::segment(dd)) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend), colour = pal[3], linewidth = 0.35) +
  geom_hline(yintercept = h_potong, linetype = "dashed", colour = pal[2], linewidth = 0.8) +
  annotate("text", x = nrow(d) * 0.98, y = h_potong, label = sprintf("cut: k = %d", k_best), hjust = 1, vjust = -0.5,
           colour = warna$teks, size = 3) +
  labs(title = "Ward (ward.D2) dendrogram", x = "Objects", y = "Height") +
  theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank())
p_hasil <- if (tata == 1) ((p2 | p1) / (p4 | p3)) + plot_layout(widths = c(1.5, 1), heights = c(1.3, 1)) else
  ((p2 / p4) | (p1 / p3)) + plot_layout(widths = c(1.5, 1))

print(p_hasil)