Tingkat kepuasan kerja guru honorer di sekolah dasar: regresi logistik ordinal

# Contoh kode R: Ordinal Logistic Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(MASS)
  library(ordinal)
})
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(monthly_salary = rnorm(n, 50, 10),
                supervisor_support = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$monthly_salary))
z2 <- as.numeric(scale(d$supervisor_support))
b1 <- runif(1, 0.9, 1.6)
b2 <- runif(1, 0.3, 0.8) * sample(c(-1, 1), 1)
laten <- b1 * z1 + b2 * z2 + rlogis(n)                    # variabel laten kontinu di balik kategori
batas <- quantile(laten, c(runif(1, 0.25, 0.4), runif(1, 0.6, 0.75)))
kat <- c("A", "B", "C")                     # urutan rendah ke tinggi
d$job_satisfaction_level <- factor(kat[findInterval(laten, batas) + 1], levels = kat, ordered = TRUE)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi logistik ordinal (proportional odds, MASS::polr) ----
d$job_satisfaction_level <- factor(d$job_satisfaction_level, levels = kat, ordered = TRUE)  # pastikan urutan kategori benar (mis. setelah read.csv)
fit <- MASS::polr(job_satisfaction_level ~ monthly_salary + supervisor_support, data = d, Hess = TRUE)
print(summary(fit))
ct <- coef(summary(fit))
pval <- 2 * pnorm(abs(ct[, "t value"]), lower.tail = FALSE)   # p-value aproksimasi normal
print(cbind(ct, "p value" = round(pval, 4)))
or_tab <- exp(cbind(OR = coef(fit), confint.default(fit)))    # odds ratio per satuan + CI Wald 95%
print(round(or_tab, 4))

# OR per 1 simpangan baku (sebanding antarprediktor)
sdx <- c(sd(d$monthly_salary), sd(d$supervisor_support))
se_b <- ct[1:2, "Std. Error"]
or_sd <- data.frame(var = c("monthly_salary", "supervisor_support"),
                    or = exp(coef(fit) * sdx), lo = exp((coef(fit) - 1.96 * se_b) * sdx),
                    hi = exp((coef(fit) + 1.96 * se_b) * sdx), p = pval[1:2])

# ---- Uji model: LR terhadap model null dan uji asumsi odds proporsional (ordinal::clm) ----
f_clm <- ordinal::clm(job_satisfaction_level ~ monthly_salary + supervisor_support, data = d)
f_null <- ordinal::clm(job_satisfaction_level ~ 1, data = d)
lr <- anova(f_null, f_clm)
print(lr)
p_lr <- lr[2, "Pr(>Chisq)"]
r2_mcf <- 1 - as.numeric(logLik(f_clm)) / as.numeric(logLik(f_null))
p_po <- tryCatch({
  f_nom <- ordinal::clm(job_satisfaction_level ~ 1, nominal = ~ monthly_salary + supervisor_support, data = d)   # ambang berbeda per prediktor
  a <- anova(f_clm, f_nom); print(a); a[2, "Pr(>Chisq)"]
}, error = function(e) {
  nt <- ordinal::nominal_test(f_clm); print(nt); min(nt[-1, "Pr(>Chi)"], na.rm = TRUE)
})
if (!is.finite(p_po)) p_po <- 1
prob <- predict(fit, type = "probs")
d$prediksi <- factor(kat[max.col(prob)], levels = kat)
tab <- table(Aktual = d$job_satisfaction_level, Prediksi = d$prediksi)
print(tab)
akurasi <- mean(as.character(d$prediksi) == as.character(d$job_satisfaction_level))


# ---- Grafik ----
# Peluang prediksi tiap kategori menurut X1 (X2 pada rata-rata), area bertumpuk
g <- data.frame(monthly_salary = seq(min(d$monthly_salary), max(d$monthly_salary), length.out = 120), supervisor_support = mean(d$supervisor_support))
pg <- predict(fit, newdata = g, type = "probs")
dg <- data.frame(x = rep(g$monthly_salary, 3), kategori = factor(rep(kat, each = nrow(g)), levels = rev(kat)),
                 p = as.numeric(pg))
p1 <- ggplot(dg, aes(x, p, fill = kategori)) +
  geom_area(alpha = 0.8, colour = warna$panel, linewidth = 0.3) +
  scale_fill_manual(values = setNames(pal[1:3], kat), breaks = kat, name = NULL,
                    labels = function(s) substr(s, 1, 18)) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Predicted category probabilities", x = "monthly_salary", y = "Probability")
# Distribusi kategori: teramati vs prediksi model
dk <- rbind(data.frame(kategori = kat, jenis = "Observed", jml = as.numeric(table(d$job_satisfaction_level))),
            data.frame(kategori = kat, jenis = "Predicted", jml = as.numeric(colSums(prob))))
dk$kategori <- factor(dk$kategori, levels = kat)
p2 <- ggplot(dk, aes(kategori, jml, fill = jenis)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.7, alpha = 0.85) +
  scale_fill_manual(values = c(Observed = pal[4], Predicted = pal[5]), name = NULL) +
  scale_x_discrete(labels = function(s) substr(s, 1, 12)) +
  labs(title = "Category distribution", x = NULL, y = "Count")
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 = "Proportional odds ratio per 1 SD (95% CI)", x = "Odds ratio (log scale)", y = NULL)
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)