Prediksi penyelesaian kursus daring terbuka (MOOC) dengan Naive Bayes

# Contoh kode R: Naive Bayes
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: fitur kontinu berdistribusi normal bersyarat kelas ----
n <- 200
positif <- function(x, m) if (m > 0) pmax(x, 0.05 * m) else x   # cegah nilai negatif yang tidak wajar
prevalensi <- runif(1, 0.3, 0.5)
kls <- rbinom(n, 1, prevalensi)
geser1 <- sample(c(-1, 1), 1) * runif(1, 0.9, 1.5)    # beda rerata antar kelas (dalam satuan SB)
geser2 <- sample(c(-1, 1), 1) * runif(1, 0.6, 1.2)
rasio_sb <- runif(1, 1.2, 1.7)                        # kelas 1 lebih bervariasi pada X2
d <- data.frame(video_watching_time = positif(50 + 10 * (geser1 * kls + rnorm(n)), 50),
                quizzes_attempted = positif(30 + 6 * (geser2 * kls + rnorm(n) * ifelse(kls == 1, rasio_sb, 1)), 30))
d$mooc_completion <- kls
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Pembagian data latih/uji (70:30, bertingkat per kelas) ----
d$kelas <- factor(d$mooc_completion, 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: Naive Bayes Gaussian (e1071::naiveBayes) ----
nb <- naiveBayes(kelas ~ video_watching_time + quizzes_attempted, data = latih)
print(nb)                                   # prior dan rerata/SB bersyarat kelas
prob_nb <- predict(nb, uji, type = "raw")[, "Yes"]
pred_nb <- factor(ifelse(prob_nb >= 0.5, "Yes", "No"), levels = c("No", "Yes"))
konf <- table(Actual = uji$kelas, Predicted = pred_nb)
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

# ---- Pembanding: regresi logistik dengan prediktor yang sama ----
lr <- glm(kelas ~ video_watching_time + quizzes_attempted, data = latih, family = binomial)
print(summary(lr))
prob_lr <- predict(lr, uji, type = "response")
akurasi_lr <- mean(ifelse(prob_lr >= 0.5, "Yes", "No") == uji$kelas)
roc_nb <- pROC::roc(uji$kelas, prob_nb, levels = c("No", "Yes"), direction = "<", quiet = TRUE)
roc_lr <- pROC::roc(uji$kelas, prob_lr, levels = c("No", "Yes"), direction = "<", quiet = TRUE)
auc_nb <- as.numeric(pROC::auc(roc_nb))
auc_lr <- as.numeric(pROC::auc(roc_lr))
uji_delong <- pROC::roc.test(roc_nb, roc_lr, method = "delong")   # beda AUC (sampel berpasangan)
print(uji_delong)
p_delong <- if (is.finite(uji_delong$p.value)) uji_delong$p.value else 1


# ---- Grafik ----
kurva_nb <- function(kolom) {
  # kurva normal yang diestimasi Naive Bayes untuk tiap kelas
  tb <- nb$tables[[kolom]]
  xs <- seq(min(d[[kolom]]), max(d[[kolom]]), length.out = 200)
  do.call(rbind, lapply(rownames(tb), function(k) data.frame(x = xs, dens = dnorm(xs, tb[k, 1], tb[k, 2]), kelas = k)))
}
k1 <- kurva_nb("video_watching_time")
k2 <- kurva_nb("quizzes_attempted")
warna_kelas <- c(No = pal[1], Yes = pal[2])
p1 <- ggplot(latih, aes(video_watching_time, fill = kelas)) +
  geom_density(alpha = 0.35, colour = NA) +
  geom_line(data = k1, aes(x, dens, colour = kelas), linewidth = 0.9, linetype = "dashed", inherit.aes = FALSE) +
  scale_fill_manual(values = warna_kelas, name = "mooc_completion") +
  scale_colour_manual(values = warna_kelas, guide = "none") +
  labs(title = "Class-conditional density", subtitle = "Dashed = Gaussian fitted by Naive Bayes", x = "video_watching_time", y = "Density") +
  theme(legend.position = "top")
p2 <- ggplot(latih, aes(quizzes_attempted, fill = kelas)) +
  geom_density(alpha = 0.35, colour = NA) +
  geom_line(data = k2, aes(x, dens, colour = kelas), linewidth = 0.9, linetype = "dashed", inherit.aes = FALSE) +
  scale_fill_manual(values = warna_kelas, guide = "none") +
  scale_colour_manual(values = warna_kelas, guide = "none") +
  labs(title = "Class-conditional density", x = "quizzes_attempted", y = "Density")
rd <- rbind(data.frame(fpr = 1 - roc_nb$specificities, tpr = roc_nb$sensitivities, model = sprintf("Naive Bayes (AUC %.3f)", auc_nb)),
            data.frame(fpr = 1 - roc_lr$specificities, tpr = roc_lr$sensitivities, model = sprintf("Logistic (AUC %.3f)", auc_lr)))
rd <- rd[order(rd$model, rd$fpr, rd$tpr), ]
p3 <- ggplot(rd, aes(fpr, tpr, colour = model)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_line(linewidth = 1) +
  scale_colour_manual(values = c(pal[4], pal[3]), name = NULL) +
  labs(title = "ROC curve (test data)", x = "False positive rate", y = "True positive rate") +
  theme(legend.position = "inside", legend.position.inside = c(0.97, 0.03), legend.justification = c(1, 0),
        legend.background = element_rect(fill = warna$panel, colour = warna$grid))
kal <- rbind(data.frame(prob = prob_nb, model = "Naive Bayes"), data.frame(prob = prob_lr, model = "Logistic"))
kal$aktual <- rep(as.numeric(uji$kelas == "Yes"), 2)
kal$bin <- cut(kal$prob, breaks = seq(0, 1, by = 0.2), include.lowest = TRUE)
kb <- aggregate(cbind(prob, aktual) ~ bin + model, data = kal, FUN = mean)
kb$jumlah <- aggregate(prob ~ bin + model, data = kal, FUN = length)$prob
p4 <- ggplot(kb, aes(prob, aktual, colour = model)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_line(linewidth = 0.9) +
  geom_point(aes(size = jumlah), alpha = 0.85) +
  scale_colour_manual(values = c(Logistic = pal[3], `Naive Bayes` = pal[4]), guide = "none") +
  scale_size_area(max_size = 4.5, guide = "none") +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(title = "Calibration (test data)", subtitle = "Point size = cases per probability bin",
       x = "Mean predicted P(Yes)", y = "Observed share Yes")
p_hasil <- if (tata == 1) (p1 / p2) | (p3 / p4) else (p1 | p2) / (p3 | p4)

print(p_hasil)