Determinan niat keluar karyawan pabrik tekstil: regresi logistik

# Contoh kode R: Binary Logistic Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  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 ----
n <- 200
d <- data.frame(pay_satisfaction = rnorm(n, 50, 10),
                tenure = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$pay_satisfaction))
z2 <- as.numeric(scale(d$tenure))
b0 <- runif(1, -0.8, 0.3)
b1 <- runif(1, 0.8, 1.5)
b2 <- runif(1, 0.25, 0.75) * sample(c(-1, 1), 1)
d$left_the_company <- rbinom(n, 1, plogis(b0 + b1 * z1 + b2 * z2))   # 1 = kejadian terjadi, 0 = tidak
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi logistik biner (glm, family binomial) ----
fit <- glm(left_the_company ~ pay_satisfaction + tenure, family = binomial, data = d)
s <- summary(fit)
print(s)
or_tab <- exp(cbind(OR = coef(fit), confint.default(fit)))   # odds ratio per satuan + CI Wald 95%
print(round(or_tab, 4))

# Odds ratio per 1 simpangan baku (lebih mudah dibandingkan antarprediktor)
sdx <- c(sd(d$pay_satisfaction), sd(d$tenure))
se_b <- s$coefficients[2:3, 2]
or_sd <- data.frame(var = c("pay_satisfaction", "tenure"),
                    or = exp(coef(fit)[2:3] * sdx),
                    lo = exp((coef(fit)[2:3] - 1.96 * se_b) * sdx),
                    hi = exp((coef(fit)[2:3] + 1.96 * se_b) * sdx),
                    p = s$coefficients[2:3, 4])
print(or_sd)

# ---- Uji kelayakan model ----
fit0 <- glm(left_the_company ~ 1, family = binomial, data = d)
lr <- anova(fit0, fit, test = "Chisq")                         # uji likelihood ratio (omnibus)
print(lr)
p_lr <- lr[2, "Pr(>Chi)"]
r2_mcf <- 1 - as.numeric(logLik(fit)) / as.numeric(logLik(fit0))   # pseudo R2 McFadden
d$peluang <- fitted(fit)
roc_obj <- pROC::roc(d$left_the_company, d$peluang, quiet = TRUE)
auc_v <- as.numeric(pROC::auc(roc_obj))
d$kelas <- as.integer(d$peluang >= 0.5)
print(table(Aktual = d$left_the_company, Prediksi = factor(d$kelas, levels = 0:1)))   # matriks klasifikasi
akurasi <- mean(d$kelas == d$left_the_company)
cat(sprintf("AUC = %.3f; akurasi (ambang 0.5) = %.3f; McFadden R2 = %.3f\n", auc_v, akurasi, r2_mcf))


# ---- Grafik ----
# Kurva peluang prediksi menurut X1 (X2 ditahan pada rata-ratanya) dengan pita 95%
g <- data.frame(pay_satisfaction = seq(min(d$pay_satisfaction), max(d$pay_satisfaction), length.out = 120), tenure = mean(d$tenure))
pr <- predict(fit, newdata = g, type = "link", se.fit = TRUE)
g$m <- plogis(pr$fit); g$lo <- plogis(pr$fit - 1.96 * pr$se.fit); g$hi <- plogis(pr$fit + 1.96 * pr$se.fit)
p1 <- ggplot() +
  geom_jitter(data = d, aes(pay_satisfaction, left_the_company), height = 0.04, width = 0, colour = pal[1], alpha = 0.35, size = 1.5) +
  geom_ribbon(data = g, aes(pay_satisfaction, ymin = lo, ymax = hi), fill = pal[2], alpha = 0.2) +
  geom_line(data = g, aes(pay_satisfaction, m), colour = pal[2], linewidth = 1.2) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1)) +
  labs(title = "Predicted probability (observed 0/1 jittered)", x = "pay_satisfaction", y = "P(event = 1)")
dr <- data.frame(fpr = 1 - roc_obj$specificities, tpr = roc_obj$sensitivities)
dr <- dr[order(dr$fpr, dr$tpr), ]
p2 <- ggplot(dr, aes(fpr, tpr)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_area(fill = pal[3], alpha = 0.15) +
  geom_path(colour = pal[3], linewidth = 1.1) +
  annotate("text", x = 0.62, y = 0.18, label = sprintf("AUC = %.3f", auc_v), colour = warna$teks, size = 3.6) +
  labs(title = "ROC curve", x = "1 - specificity", y = "Sensitivity")
p3 <- ggplot(or_sd, aes(or, var)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = warna$redup) +
  geom_segment(aes(x = lo, xend = hi, yend = var), colour = pal[1], linewidth = 1.1) +
  geom_point(colour = pal[2], size = 3.2) +
  geom_text(aes(label = sprintf("%.2f", or)), vjust = -1.1, colour = warna$teks, size = 3.2) +
  scale_x_log10() +
  labs(title = "Odds ratio per 1 SD (95% CI, log scale)", x = "Odds ratio", y = NULL)
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)