Prediksi nilai seumur hidup pelanggan ritel daring menggunakan XGBoost

# Contoh kode R: XGBoost
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(xgboost)
})
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: Y kontinu dengan hubungan nonlinear dan interaksi ----
n <- max(200, 300)
d <- data.frame(transaction_frequency = round(rnorm(n, 50, 10), 3),
                average_spend_per_order = round(rnorm(n, 30, 6), 3),
                fitur_3 = round(runif(n, 0, 10), 2),             # fitur pendukung (pengaruh ambang)
                fitur_4 = round(rnorm(n, 50, 12), 2),            # fitur pendukung (pengaruh lemah)
                fitur_5 = round(rnorm(n, 20, 5), 2))             # fitur pengganggu (tanpa pengaruh)
z1 <- as.numeric(scale(d$transaction_frequency)); z2 <- as.numeric(scale(d$average_spend_per_order))
z4 <- as.numeric(scale(d$fitur_4))
a1 <- runif(1, 0.7, 1.1); a2 <- runif(1, 0.4, 0.7) * sample(c(-1, 1), 1)
f <- a1 * tanh(1.3 * z1) + a2 * z2 + runif(1, 0.2, 0.4) * z2^2 * sign(a2) +   # efek jenuh dan kuadratik
  runif(1, 0.3, 0.5) * (d$fitur_3 > 6) + 0.25 * z1 * z2 + 0.15 * z4           # efek ambang dan interaksi
f <- (f - mean(f)) / sd(f)
d$customer_lifetime_value <- round(100 + 15 * (f * runif(1, 0.85, 0.95) + rnorm(n, 0, runif(1, 0.3, 0.45))), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Bagi data latih 75% dan uji 25% ----
fitur <- c("transaction_frequency", "average_spend_per_order", "fitur_3", "fitur_4", "fitur_5")
id_latih <- sample(n, round(0.75 * n))
latih <- d[id_latih, ]; uji <- d[-id_latih, ]
dlatih <- xgb.DMatrix(data = as.matrix(latih[, fitur]), label = latih$customer_lifetime_value)
duji <- xgb.DMatrix(data = as.matrix(uji[, fitur]), label = uji$customer_lifetime_value)

# ---- Validasi silang 5 lipatan untuk memilih jumlah iterasi (early stopping) ----
param <- list(objective = "reg:squarederror", eta = 0.05, max_depth = 4, subsample = 0.8,
              colsample_bytree = 0.8, min_child_weight = 3, nthread = 1,
              base_score = mean(latih$customer_lifetime_value))                           # titik awal = rata-rata Y latih
cv <- xgb.cv(params = param, data = dlatih, nrounds = 1000, nfold = 5, early_stopping_rounds = 30,
             metrics = "rmse", verbose = 0)
best <- cv$best_iteration
cat("Jumlah iterasi terbaik (CV):", best, "\n")

# ---- Model akhir dan evaluasi pada data uji ----
fit <- xgb.train(params = param, data = dlatih, nrounds = best, verbose = 0)
pred <- predict(fit, duji)
ukur <- function(y, p) c(RMSE = sqrt(mean((y - p)^2)), MAE = mean(abs(y - p)), R2 = 1 - sum((y - p)^2) / sum((y - mean(y))^2))
m_xgb <- ukur(uji$customer_lifetime_value, pred)
print(round(m_xgb, 4))

# ---- Pembanding: regresi linear berganda ----
fit_lm <- lm(customer_lifetime_value ~ transaction_frequency + average_spend_per_order + fitur_3 + fitur_4 + fitur_5, data = latih)
pred_lm <- predict(fit_lm, uji)
m_lm <- ukur(uji$customer_lifetime_value, pred_lm)
print(rbind(XGBoost = m_xgb, Linear = m_lm))

# ---- Kepentingan fitur (gain) ----
imp <- xgb.importance(feature_names = fitur, model = fit)
print(imp)


# ---- Grafik ----
dp <- rbind(data.frame(aktual = uji$customer_lifetime_value, prediksi = pred_lm, model = "Linear"),
            data.frame(aktual = uji$customer_lifetime_value, prediksi = pred, model = "XGBoost"))
dp$model <- factor(dp$model, levels = c("Linear", "XGBoost"))
p1 <- ggplot(dp, aes(aktual, prediksi, colour = model)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_point(alpha = 0.6, size = 1.6) +
  scale_colour_manual(values = c(Linear = warna$redup, XGBoost = pal[1]), name = NULL) +
  labs(title = "Predicted vs actual (test set)", x = "Actual customer_lifetime_value", y = "Predicted customer_lifetime_value") +
  theme(legend.position = "bottom")
nama_tampil <- c("transaction_frequency" = "transaction_frequency", "average_spend_per_order" = "average_spend_per_order", fitur_3 = "Feature 3", fitur_4 = "Feature 4", fitur_5 = "Feature 5")
nama_tampil[] <- make.unique(nama_tampil, sep = " ")              # cegah label kembar
di <- data.frame(fitur = nama_tampil[imp$Feature], gain = imp$Gain)
di$fitur <- factor(di$fitur, levels = di$fitur[order(di$gain)])
p2 <- ggplot(di, aes(gain, fitur)) +
  geom_col(fill = pal[2], width = 0.65, alpha = 0.9) +
  geom_text(aes(label = sprintf("%.1f%%", 100 * gain)), hjust = -0.15, colour = warna$teks, size = 3) +
  scale_x_continuous(labels = function(v) paste0(100 * v, "%"), expand = expansion(mult = c(0, 0.2))) +
  labs(title = "Feature importance (gain)", x = "Share of total gain", y = NULL)
ev <- as.data.frame(cv$evaluation_log)
dl <- rbind(data.frame(iter = ev$iter, rmse = ev$train_rmse_mean, sd = ev$train_rmse_std, set = "CV train"),
            data.frame(iter = ev$iter, rmse = ev$test_rmse_mean, sd = ev$test_rmse_std, set = "CV validation"))
p3 <- ggplot(dl, aes(iter, rmse, colour = set, fill = set)) +
  geom_ribbon(aes(ymin = rmse - sd, ymax = rmse + sd), alpha = 0.15, colour = NA) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = best, linetype = "dashed", colour = warna$redup) +
  scale_colour_manual(values = c("CV train" = pal[3], "CV validation" = pal[4]), name = NULL) +
  scale_fill_manual(values = c("CV train" = pal[3], "CV validation" = pal[4]), name = NULL) +
  labs(title = sprintf("Learning curve (best = %d rounds)", best), x = "Boosting round", y = "RMSE") +
  theme(legend.position = "bottom")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else (p2 / p3) | p1

print(p_hasil)