Pengaruh kepemimpinan melayani terhadap kinerja: mediasi kepuasan kerja dan moderasi iklim kerja

# Contoh kode R: Mediation and Moderation Analysis
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(lavaan)
})
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: X (variabel bebas), M (mediator), W (moderator), Y (hasil) ----
n <- 200
zx <- rnorm(n); zw <- rnorm(n)
a <- runif(1, 0.35, 0.6)                                        # jalur X -> M
b <- runif(1, 0.25, 0.45)                                       # jalur M -> Y
cp <- runif(1, 0.1, 0.3)                                        # efek langsung X -> Y
h <- runif(1, 0.12, 0.25) * sample(c(-1, 1), 1)                 # interaksi X x W
zm <- a * zx + rnorm(n, 0, sqrt(1 - a^2))
zy <- cp * zx + b * zm + 0.2 * zw + h * zx * zw + rnorm(n, 0, runif(1, 0.6, 0.8))
zy <- as.numeric(scale(zy))
d <- data.frame(servant_leadership = round(50 + 10 * zx, 2), job_satisfaction = round(30 + 6 * zm, 2),
                x3 = round(20 + 5 * zw, 2), employee_performance = round(100 + 15 * zy, 2))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Mediasi: X -> M -> Y dengan SE bootstrap (variabel dibakukan agar koefisien sebanding) ----
ds <- as.data.frame(scale(d))
model <- "
  job_satisfaction ~ a*servant_leadership
  employee_performance ~ cp*servant_leadership + b*job_satisfaction
  tidak_langsung := a*b
  total := cp + a*b
"
fit <- sem(model, data = ds, se = "bootstrap", bootstrap = 500)
pe <- parameterEstimates(fit, boot.ci.type = "perc")
print(pe[pe$op %in% c("~", ":="), c("lhs", "op", "rhs", "label", "est", "se", "pvalue", "ci.lower", "ci.upper")])
ambil <- function(lb) pe[pe$label == lb | pe$lhs == lb, ][1, ]
e_a <- ambil("a"); e_b <- ambil("b"); e_c <- ambil("cp"); e_i <- ambil("tidak_langsung"); e_t <- ambil("total")
bb <- lavInspect(fit, "boot")                                   # sampel bootstrap koefisien
ab_boot <- bb[, "a"] * bb[, "b"]

# ---- Moderasi: regresi dengan interaksi (prediktor dipusatkan) ----
dm <- data.frame(yy = d$employee_performance, xp = d$servant_leadership - mean(d$servant_leadership), wp = d$x3 - mean(d$x3))
mod <- lm(yy ~ xp * wp, data = dm)
print(summary(mod))
sw <- sd(d$x3)
V <- vcov(mod)
lereng <- function(w) {                                         # simple slope X pada nilai W tertentu (terpusat)
  est <- coef(mod)[["xp"]] + coef(mod)[["xp:wp"]] * w
  se <- sqrt(V["xp", "xp"] + w^2 * V["xp:wp", "xp:wp"] + 2 * w * V["xp", "xp:wp"])
  c(est = est, se = se, p = 2 * pt(-abs(est / se), df.residual(mod)))
}
ss <- t(sapply(c(-sw, 0, sw), lereng)); rownames(ss) <- c("W -1 SD", "W mean", "W +1 SD")
print(round(ss, 4))
# Johnson-Neyman: nilai W saat simple slope X berubah signifikan
wg <- seq(min(dm$wp), max(dm$wp), length.out = 200)
jn <- as.data.frame(t(sapply(wg, lereng))); jn$w <- wg + mean(d$x3)
p_int <- summary(mod)$coefficients["xp:wp", 4]
b_int <- coef(lm(scale(yy) ~ scale(xp) * scale(wp), data = dm))[4]   # koefisien interaksi terstandar
sig_pct <- 100 * mean(jn$p < 0.05)


# ---- Grafik ----
bintang <- function(p) ifelse(p < 0.001, "***", ifelse(p < 0.01, "**", ifelse(p < 0.05, "*", " ns")))
bungkus <- function(s) paste(head(strwrap(s, 14), 2), collapse = "\n")   # nama maks. 2 baris
singkat <- function(s, w = 18) {                                # nama pendek untuk legenda
  k <- strsplit(s, " ")[[1]]
  if (nchar(s) > w) k <- k[!tolower(k) %in% c("of", "to", "the", "and", "for", "in", "on", "with")]
  i <- length(k)
  while (nchar(paste(k, collapse = " ")) > w && i >= 1) {
    if (nchar(k[i]) > 5) k[i] <- paste0(substr(k[i], 1, 4), ".")
    i <- i - 1
  }
  substr(paste(k, collapse = " "), 1, w)
}
nd <- data.frame(x = c(0, 1, 2), y = c(0, 1, 0), lab = c(bungkus("servant_leadership"), bungkus("job_satisfaction"), bungkus("employee_performance")))
pan <- data.frame(x = c(0.2, 1.2, 0.36), y = c(0.2, 0.8, 0), xend = c(0.8, 1.8, 1.64), yend = c(0.8, 0.2, 0),
                  lab = c(sprintf("a = %.2f%s", e_a$est, bintang(e_a$pvalue)), sprintf("b = %.2f%s", e_b$est, bintang(e_b$pvalue)),
                          sprintf("c' = %.2f%s", e_c$est, bintang(e_c$pvalue))),
                  lx = c(0.28, 1.72, 1), ly = c(0.62, 0.62, -0.1))
p1 <- ggplot() +
  geom_segment(data = pan, aes(x, y, xend = xend, yend = yend), colour = warna$redup, linewidth = 0.9,
               arrow = arrow(length = unit(0.13, "inches"), type = "closed")) +
  geom_text(data = pan, aes(lx, ly, label = lab), colour = warna$teks, size = 3.2, fontface = "bold") +
  geom_tile(data = nd, aes(x, y), width = 0.8, height = 0.38, fill = pal[1:3]) +
  geom_text(data = nd, aes(x, y, label = lab), colour = warna$panel, fontface = "bold", size = 2.8, lineheight = 0.9) +
  annotate("text", x = 1, y = -0.34, colour = warna$teks, size = 3,
           label = sprintf("Indirect a*b = %.3f [%.3f, %.3f]", e_i$est, e_i$ci.lower, e_i$ci.upper)) +
  coord_cartesian(xlim = c(-0.45, 2.45), ylim = c(-0.42, 1.22)) +
  theme_void() + theme(plot.background = element_rect(fill = warna$panel, colour = NA),
                       plot.title = element_text(face = "bold", colour = warna$teks, margin = margin(4, 0, 4, 6))) +
  labs(title = "Mediation model (standardised)")
xs <- seq(min(dm$xp), max(dm$xp), length.out = 50)
dsl <- do.call(rbind, lapply(1:3, function(i) data.frame(
  x = xs + mean(d$servant_leadership), y = predict(mod, data.frame(xp = xs, wp = c(-sw, 0, sw)[i])),
  tingkat = c("Low W (-1 SD)", "Mean W", "High W (+1 SD)")[i])))
dsl$tingkat <- factor(dsl$tingkat, levels = c("Low W (-1 SD)", "Mean W", "High W (+1 SD)"))
p2 <- ggplot() +
  geom_point(data = d, aes(servant_leadership, employee_performance), colour = warna$redup, alpha = 0.25, size = 1.1) +
  geom_line(data = dsl, aes(x, y, colour = tingkat), linewidth = 1.1) +
  scale_colour_manual(values = pal[c(4, 3, 5)], name = singkat("x3")) +
  labs(title = "Simple slopes at low, mean and high W", x = "servant_leadership", y = "employee_performance")
p3 <- ggplot(data.frame(ab = ab_boot), aes(ab)) +
  geom_histogram(bins = 30, fill = pal[2], colour = warna$panel, alpha = 0.85) +
  geom_vline(xintercept = c(e_i$ci.lower, e_i$ci.upper), linetype = "dashed", colour = warna$teks) +
  geom_vline(xintercept = 0, colour = warna$redup) +
  labs(title = "Bootstrap distribution of a*b (500)", x = "Indirect effect (std.)", y = "Count")
p4 <- ggplot(jn, aes(w, est)) +
  geom_ribbon(aes(ymin = est - 1.96 * se, ymax = est + 1.96 * se), fill = pal[1], alpha = 0.2) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_line(aes(colour = p < 0.05, group = 1), linewidth = 1.1) +
  scale_colour_manual(values = c(`TRUE` = pal[1], `FALSE` = pal[2]), labels = c(`TRUE` = "p < .05", `FALSE` = "n.s."),
                      name = NULL, drop = FALSE) +
  labs(title = "Johnson-Neyman: slope of X across W", x = "x3", y = "Simple slope of X")
p_hasil <- if (tata == 1) (p1 | p2) / (p3 | p4) else (p1 | p3) / (p2 | p4)

print(p_hasil)