Hubungan nonlinear usia dan produktivitas pekerja pabrik: pendekatan GAM

# Contoh kode R: Generalized Additive Model (GAM)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(mgcv)
})
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: hubungan nonlinear ----
n <- 200
d <- data.frame(worker_age = runif(n, 50 - 1.9 * 10, 50 + 1.9 * 10),
                tenure = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$worker_age))
z2 <- as.numeric(scale(d$tenure))
# Bentuk kurva dipilih acak: U terbalik, gelombang, atau jenuh (saturasi)
bentuk1 <- sample(1:3, 1)
f1 <- switch(bentuk1, -(z1 - runif(1, -0.4, 0.4))^2, sin(runif(1, 1.2, 1.7) * z1), log1p(exp(2.2 * z1)) - 0.25 * z1^2)
f2 <- if (runif(1) < 0.5) tanh(runif(1, 1, 1.8) * z2) else exp(0.6 * z2)
f1 <- as.numeric(scale(f1)); f2 <- as.numeric(scale(f2))
y <- runif(1, 0.6, 0.85) * f1 + runif(1, 0.3, 0.5) * sample(c(-1, 1), 1) * f2 + rnorm(n, 0, runif(1, 0.45, 0.65))
d$hourly_productivity <- round(100 + 15 * y, 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: Generalized Additive Model (mgcv::gam, REML) ----
fit <- gam(hourly_productivity ~ s(worker_age) + s(tenure), data = d, method = "REML")
s <- summary(fit)
print(s)
print(k.check(fit))                                           # cek kecukupan dimensi basis k (k-index)
lin <- gam(hourly_productivity ~ worker_age + tenure, data = d, method = "REML")  # pembanding: model linear
aic <- c(Linear = AIC(lin), GAM = AIC(fit))
print(aic)
uji <- anova(lin, fit, test = "F")                            # uji F aproksimasi: linear vs GAM
print(uji)
edf <- unname(s$s.table[, "edf"])
p_s <- unname(s$s.table[, "p-value"])


# ---- Grafik ----
# Efek parsial tiap smooth (skala Y, variabel lain pada rata-rata) dengan pita 95% dan residual parsial
efek <- function(v, w, lab) {
  g <- data.frame(worker_age = rep(mean(d$worker_age), 150), tenure = rep(mean(d$tenure), 150))
  g[[v]] <- seq(min(d[[v]]), max(d[[v]]), length.out = 150)
  pr <- predict(fit, newdata = g, type = "terms", se.fit = TRUE)
  j <- match(paste0("s(", v, ")"), colnames(pr$fit))
  gg <- data.frame(x = g[[v]], m = pr$fit[, j], lo = pr$fit[, j] - 1.96 * pr$se.fit[, j], hi = pr$fit[, j] + 1.96 * pr$se.fit[, j])
  tr <- predict(fit, type = "terms")
  dp <- data.frame(x = d[[v]], r = tr[, j] + residuals(fit))   # residual parsial
  ggplot() +
    geom_hline(yintercept = 0, colour = warna$redup, linewidth = 0.4) +
    geom_point(data = dp, aes(x, r), colour = warna$redup, alpha = 0.3, size = 1.1) +
    geom_ribbon(data = gg, aes(x, ymin = lo, ymax = hi), fill = w, alpha = 0.25) +
    geom_line(data = gg, aes(x, m), colour = w, linewidth = 1.3) +
    geom_rug(data = dp, aes(x), sides = "b", colour = w, alpha = 0.3, length = unit(0.02, "npc")) +
    labs(title = sprintf("Partial effect s(%s), edf = %.2f", substr(lab, 1, 14), s$s.table[j, "edf"]),
         x = if (v == "worker_age") "worker_age" else "tenure", y = "Partial effect")
}
p1 <- efek("worker_age", pal[1], "worker_age")
p2 <- efek("tenure", pal[2], "tenure")
d$prediksi <- fitted(fit)
p3 <- ggplot(d, aes(prediksi, hourly_productivity)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_point(colour = pal[3], alpha = 0.55, size = 1.5) +
  annotate("text", x = min(d$prediksi), y = max(d$hourly_productivity), hjust = 0, vjust = 1, colour = warna$teks, size = 3.3,
           label = sprintf("Deviance explained = %.1f%%", 100 * summary(fit)$dev.expl)) +
  labs(title = "Predicted vs actual", x = "Predicted hourly_productivity", y = "hourly_productivity")
p_hasil <- if (tata == 1) (p1 | p2) / p3 else p1 | (p2 / p3)

print(p_hasil)