Perbedaan kepuasan kerja, komitmen, dan niat bertahan antar tiga generasi karyawan: MANOVA

# Contoh kode R: MANOVA
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(MASS)
  library(scales)
})
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: tiga variabel terikat yang diukur bersama pada tiga kelompok ----
n <- 200
lv <- c("A", "B", "C")
kelompok <- factor(sample(rep(lv, length.out = n)), levels = lv)
r12 <- runif(1, 0.2, 0.55); r13 <- runif(1, -0.3, 0.4); r23 <- runif(1, 0.1, 0.45)
S <- matrix(c(1, r12, r13, r12, 1, r23, r13, r23, 1), 3)              # korelasi antar-variabel terikat
galat <- MASS::mvrnorm(n, mu = rep(0, 3), Sigma = S)
geser <- matrix(runif(9, -0.55, 0.55), 3)                              # profil rata-rata tiap kelompok (dalam SD)
geser[, 1] <- geser[, 1] + c(-0.35, 0, 0.35)[sample(3)]                # minimal satu DV jelas berbeda
geser <- sweep(geser, 2, colMeans(geser))
z <- geser[as.integer(kelompok), ] + galat * runif(1, 0.75, 0.95)
d <- data.frame(kelompok = kelompok,
                job_satisfaction = round(100 + 15 * z[, 1], 2),
                organisational_commitment = round(50 + 10 * z[, 2], 2),
                intention_to_stay = round(30 + 6 * z[, 3], 2))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: MANOVA satu arah ----
fit <- manova(cbind(job_satisfaction, organisational_commitment, intention_to_stay) ~ kelompok, data = d)
s_pillai <- summary(fit, test = "Pillai")
s_wilks <- summary(fit, test = "Wilks")
print(s_pillai); print(s_wilks)

# ---- ANOVA lanjutan per variabel terikat + partial eta kuadrat ----
s_aov <- summary.aov(fit)
print(s_aov)
lanjut <- do.call(rbind, lapply(s_aov, function(a) data.frame(F = a[1, "F value"], p = a[1, "Pr(>F)"],
                                                              peta2 = a[1, "Sum Sq"] / sum(a[, "Sum Sq"]))))
lanjut$dv <- make.unique(c("job_satisfaction", "organisational_commitment", "intention_to_stay"), sep = " ")
print(lanjut)

# ---- Analisis diskriminan linear (visualisasi pemisahan kelompok) ----
ld <- MASS::lda(kelompok ~ job_satisfaction + organisational_commitment + intention_to_stay, data = d)
print(ld)
skor_ld <- predict(ld)
akurasi <- mean(skor_ld$class == d$kelompok)
cat("Akurasi klasifikasi LDA (resubstitusi):", round(akurasi, 3), "\n")


# ---- Grafik ----
pendek <- function(s, w = 14) ifelse(nchar(s) > w, paste0(trimws(substr(s, 1, w - 1)), "."), s)
dv <- c("job_satisfaction", "organisational_commitment", "intention_to_stay")
lab_dv <- make.unique(c("job_satisfaction", "organisational_commitment", "intention_to_stay"), sep = " ")
panjang <- do.call(rbind, lapply(seq_along(dv), function(j) {
  zz <- as.numeric(scale(d[[dv[j]]]))                                  # skala terstandar (z) karena satuan berbeda
  a <- aggregate(zz, list(kelompok = d$kelompok), function(v) c(m = mean(v), se = sd(v) / sqrt(length(v))))
  data.frame(dv = lab_dv[j], kelompok = a$kelompok, m = a$x[, "m"], se = a$x[, "se"])
}))
panjang$dv <- factor(panjang$dv, levels = unique(lab_dv))
p1 <- ggplot(panjang, aes(dv, m, colour = kelompok)) +
  geom_hline(yintercept = 0, colour = warna$redup, linetype = "dashed") +
  geom_pointrange(aes(ymin = m - 1.96 * se, ymax = m + 1.96 * se), position = position_dodge(width = 0.55),
                  size = 0.45, linewidth = 0.9) +
  scale_colour_manual(values = pal[1:3], name = NULL, labels = pendek) +
  scale_x_discrete(labels = scales::label_wrap(12)) +
  labs(title = "Group means of each outcome (z, 95% CI)", x = NULL, y = "Standardised mean (z)") +
  theme(legend.position = "bottom")
dl <- data.frame(LD1 = skor_ld$x[, 1], LD2 = skor_ld$x[, 2], kelompok = d$kelompok)
prop <- 100 * ld$svd^2 / sum(ld$svd^2)
p2 <- ggplot(dl, aes(LD1, LD2, colour = kelompok)) +
  geom_point(alpha = 0.55, size = 1.4) +
  stat_ellipse(level = 0.68, linewidth = 1) +
  geom_point(data = aggregate(cbind(LD1, LD2) ~ kelompok, dl, mean), shape = 23, size = 3.2,
             aes(fill = kelompok), colour = warna$teks, stroke = 1) +
  scale_colour_manual(values = pal[1:3], guide = "none") +
  scale_fill_manual(values = pal[1:3], guide = "none") +
  labs(title = "Discriminant function scores (68% ellipses)",
       x = sprintf("LD1 (%.1f%%)", prop[1]), y = sprintf("LD2 (%.1f%%)", prop[2]))
kr <- cor(d[, dv])
dk <- expand.grid(a = seq_along(dv), b = seq_along(dv))
dk$r <- kr[cbind(dk$a, dk$b)]
dk$a <- factor(lab_dv[dk$a], levels = unique(lab_dv)); dk$b <- factor(lab_dv[dk$b], levels = rev(unique(lab_dv)))
p3 <- ggplot(dk, aes(a, b, fill = r)) +
  geom_tile(colour = warna$panel, linewidth = 1.2) +
  geom_label(aes(label = sprintf("%.2f", r)), fill = warna$panel, colour = warna$teks, size = 3.3,
             label.size = 0, alpha = 0.85) +
  scale_fill_gradient2(low = pal[1], mid = warna$grid, high = pal[2], limits = c(-1, 1), guide = "none") +
  scale_x_discrete(labels = scales::label_wrap(12)) + scale_y_discrete(labels = scales::label_wrap(12)) +
  labs(title = "Correlation among outcomes", x = NULL, y = NULL) +
  theme(panel.grid = element_blank())
dd <- data.frame(dv = factor(lanjut$dv, levels = unique(lab_dv)), peta2 = lanjut$peta2, p = lanjut$p)
p4 <- ggplot(dd, aes(dv, peta2)) +
  geom_col(fill = pal[4], width = 0.6, alpha = 0.85) +
  geom_text(aes(label = ifelse(p < 0.001, "p<0.001", sprintf("p=%.3f", p))), vjust = -0.4, colour = warna$teks, size = 3) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  scale_x_discrete(labels = scales::label_wrap(12)) +
  labs(title = "Follow-up ANOVA: partial eta squared", x = NULL, y = "Partial eta squared")
p_hasil <- if (tata == 1) (p1 | p2) / (p3 | p4) else (p2 | p1) / (p4 | p3)

print(p_hasil)