Klasifikasi karyawan berisiko burnout dengan Support Vector Machine

# Contoh kode R: Support Vector Machine
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(e1071)
  library(pROC)
})
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: kelas biner dengan batas keputusan nonlinear ----
n <- 200
positif <- function(x, m) if (m > 0) pmax(x, 0.05 * m) else x   # cegah nilai negatif yang tidak wajar
d <- data.frame(weekly_work_hours = positif(rnorm(n, 50, 10), 50),
                social_support = positif(rnorm(n, 30, 6), 30))
z1 <- as.numeric(scale(d$weekly_work_hours))
z2 <- as.numeric(scale(d$social_support))
arah <- sample(c(-1, 1), 2, replace = TRUE)
skor <- arah[1] * runif(1, 0.8, 1.3) * z1 + arah[2] * runif(1, 0.5, 1) * z2 +
  runif(1, 0.6, 1.1) * (z1^2 - 1) * sample(c(-1, 1), 1) + runif(1, 0.3, 0.7) * z1 * z2
skor <- skor - quantile(skor, runif(1, 0.4, 0.6))
d$at_risk_of_burnout <- rbinom(n, 1, plogis(runif(1, 2.5, 4) * skor))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Pembagian data latih/uji (70:30, bertingkat per kelas) ----
d$kelas <- factor(d$at_risk_of_burnout, levels = c(0, 1), labels = c("No", "Yes"))
idx <- unlist(lapply(split(seq_len(n), d$kelas), function(i) sample(i, round(0.7 * length(i)))))
latih <- d[idx, ]
uji <- d[-idx, ]

# ---- Analisis: SVM kernel radial, tuning cost dan gamma dengan CV 5 fold ----
grid_cost <- 2^(-2:5)
grid_gamma <- 2^(-4:1)
tun <- tune.svm(kelas ~ weekly_work_hours + social_support, data = latih, kernel = "radial", cost = grid_cost, gamma = grid_gamma,
                tunecontrol = tune.control(sampling = "cross", cross = 5))
print(summary(tun))
c_best <- tun$best.parameters$cost
g_best <- tun$best.parameters$gamma
svm_fit <- svm(kelas ~ weekly_work_hours + social_support, data = latih, kernel = "radial", cost = c_best, gamma = g_best,
               probability = TRUE)
print(summary(svm_fit))

# ---- Evaluasi pada data uji ----
pr <- predict(svm_fit, uji, probability = TRUE)
prob <- attr(pr, "probabilities")[, "Yes"]
konf <- table(Actual = uji$kelas, Predicted = pr)
print(konf)
tp <- konf["Yes", "Yes"]; fp <- konf["No", "Yes"]; fn <- konf["Yes", "No"]
akurasi <- sum(diag(konf)) / sum(konf)
presisi <- if (tp + fp > 0) tp / (tp + fp) else 0
recall <- if (tp + fn > 0) tp / (tp + fn) else 0
f1 <- if (presisi + recall > 0) 2 * presisi * recall / (presisi + recall) else 0
roc_svm <- pROC::roc(uji$kelas, prob, levels = c("No", "Yes"), direction = "<", quiet = TRUE)
auc_svm <- as.numeric(pROC::auc(roc_svm))
n_sv <- nrow(svm_fit$SV)
cat(sprintf("Akurasi = %.3f, presisi = %.3f, recall = %.3f, F1 = %.3f, AUC = %.3f, SV = %d\n",
            akurasi, presisi, recall, f1, auc_svm, n_sv))


# ---- Grafik ----
gx <- seq(min(d$weekly_work_hours), max(d$weekly_work_hours), length.out = 120)
gy <- seq(min(d$social_support), max(d$social_support), length.out = 120)
bidang <- expand.grid(weekly_work_hours = gx, social_support = gy)
bidang$prob <- attr(predict(svm_fit, bidang, probability = TRUE), "probabilities")[, "Yes"]
latih$sv <- seq_len(nrow(latih)) %in% svm_fit$index
p1 <- ggplot() +
  geom_raster(data = bidang, aes(weekly_work_hours, social_support, fill = prob), alpha = 0.55, interpolate = TRUE) +
  geom_contour(data = bidang, aes(weekly_work_hours, social_support, z = prob), breaks = 0.5, colour = warna$teks, linewidth = 0.8) +
  geom_point(data = latih[latih$sv, ], aes(weekly_work_hours, social_support), shape = 21, size = 3.4, stroke = 0.7,
             colour = warna$teks, fill = NA) +
  geom_point(data = latih, aes(weekly_work_hours, social_support, colour = kelas), size = 1.5, alpha = 0.85) +
  scale_fill_gradient2(low = pal[1], mid = warna$panel, high = pal[2], midpoint = 0.5, guide = "none") +
  scale_colour_manual(values = c(No = pal[1], Yes = pal[2]), name = "at_risk_of_burnout") +
  coord_cartesian(expand = FALSE) +
  labs(title = "SVM decision surface (training data)", subtitle = "Shading = P(Yes), line = 0.5 boundary, rings = support vectors",
       x = "weekly_work_hours", y = "social_support")
pf <- tun$performances
pf$akurasi <- 1 - pf$error
pf$terbaik <- pf$cost == c_best & pf$gamma == g_best
p2 <- ggplot(pf, aes(factor(log2(gamma)), factor(log2(cost)), fill = akurasi)) +
  geom_tile(colour = warna$panel, linewidth = 0.6, alpha = 0.8) +
  geom_tile(data = pf[pf$terbaik, ], colour = warna$teks, fill = NA, linewidth = 1.1) +
  geom_text(aes(label = sprintf("%.2f", akurasi)), size = 2.4, colour = warna$teks) +
  scale_fill_gradient(low = warna$panel, high = pal[3], name = "CV acc.") +
  labs(title = "CV accuracy (cost x gamma)", x = "log2(gamma)", y = "log2(cost)") +
  theme(legend.position = "none", panel.grid = element_blank())
rd <- data.frame(fpr = 1 - roc_svm$specificities, tpr = roc_svm$sensitivities)
rd <- rd[order(rd$fpr, rd$tpr), ]
p3 <- ggplot(rd, aes(fpr, tpr)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_area(fill = pal[4], alpha = 0.15) +
  geom_line(colour = pal[4], linewidth = 1.1) +
  annotate("text", x = 0.62, y = 0.12, label = sprintf("AUC = %.3f", auc_svm), colour = warna$teks, size = 3.4) +
  labs(title = "ROC curve (test data)", x = "False positive rate", y = "True positive rate")
p_hasil <- if (tata == 1) (p1 | (p2 / p3)) + plot_layout(widths = c(1.35, 1)) else (p1 / (p2 | p3)) + plot_layout(heights = c(1.3, 1))

print(p_hasil)