Pengaruh pengalaman dan sertifikasi terhadap sebaran gaji profesional teknologi: regresi kuantil

# Contoh kode R: Quantile Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(quantreg)
})
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: Y menceng ke kanan dan sebarannya melebar seiring X1 ----
n <- 200
d <- data.frame(work_experience = rnorm(n, 50, 10),
                certifications = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$work_experience))
z2 <- as.numeric(scale(d$certifications))
b1 <- runif(1, 0.3, 0.55)
b2 <- runif(1, 0.12, 0.3) * sample(c(-1, 1), 1)
skala_galat <- exp(runif(1, 0.3, 0.5) * z1)                 # heteroskedastis: galat membesar saat X1 naik
galat <- (rgamma(n, shape = 2, rate = 2) - 1) * runif(1, 0.6, 0.8)   # galat menceng ke kanan, rata-rata 0
y <- 100 + 15 * (b1 * z1 + b2 * z2 + skala_galat * galat)
d$monthly_salary <- round(pmax(y, 100 - 2.5 * 15, 0.02 * abs(100)), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi kuantil (quantreg::rq) dan OLS sebagai pembanding ----
taus <- c(0.1, 0.25, 0.5, 0.75, 0.9)
fit <- rq(monthly_salary ~ work_experience + certifications, tau = taus, data = d)
sq <- summary(fit, se = "nid")                              # SE metode Hendricks-Koenker (sandwich)
print(sq)
ols <- lm(monthly_salary ~ work_experience + certifications, data = d)
print(summary(ols))
koef_q <- sapply(sq, function(s) s$coefficients[2, 1])       # koefisien X1 tiap kuantil
p_q <- sapply(sq, function(s) s$coefficients[2, 4])

# Uji kesamaan slope antarkuantil (Wald, Koenker-Bassett) untuk X1
fq <- lapply(c(0.1, 0.5, 0.9), function(t) rq(monthly_salary ~ work_experience + certifications, tau = t, data = d))
uj <- anova(fq[[1]], fq[[2]], fq[[3]], joint = FALSE)
print(uj)
p_sama <- as.numeric(uj$table[1, "pvalue"])

# Pseudo R1 Koenker-Machado pada median
f50 <- fq[[2]]
f50_0 <- rq(monthly_salary ~ 1, tau = 0.5, data = d)
r1 <- 1 - sum(abs(residuals(f50))) / sum(abs(residuals(f50_0)))


# ---- Grafik ----
# Garis kuantil terhadap X1 (X2 ditahan pada rata-ratanya)
cf <- coef(fit)
gl <- data.frame(tau = factor(taus), a = cf[1, ] + cf[3, ] * mean(d$certifications), b = cf[2, ])
ab_ols <- c(coef(ols)[1] + coef(ols)[3] * mean(d$certifications), coef(ols)[2])
warna_tau <- setNames(c(pal[1], pal[3], pal[2], pal[3], pal[1]), levels(gl$tau))
p1 <- ggplot(d, aes(work_experience, monthly_salary)) +
  geom_point(colour = warna$redup, alpha = 0.35, size = 1.4) +
  geom_abline(data = gl, aes(intercept = a, slope = b, colour = tau, linetype = tau), linewidth = 1) +
  geom_abline(intercept = ab_ols[1], slope = ab_ols[2], colour = pal[4], linewidth = 0.9, linetype = "dotted") +
  scale_colour_manual(values = warna_tau, name = "Quantile") +
  scale_linetype_manual(values = c("dashed", "longdash", "solid", "longdash", "dashed"), name = "Quantile") +
  labs(title = "Quantile lines 0.1 to 0.9 (dotted = OLS)", x = "work_experience", y = "monthly_salary")
# Proses koefisien lintas tau dengan pita CI 95% dan garis OLS
tg <- seq(0.05, 0.95, by = 0.05)
proses <- do.call(rbind, lapply(tg, function(t) {
  s <- tryCatch(summary(rq(monthly_salary ~ work_experience + certifications, tau = t, data = d), se = "nid")$coefficients,
                error = function(e) summary(rq(monthly_salary ~ work_experience + certifications, tau = t, data = d), se = "iid")$coefficients)
  data.frame(tau = t, var = c("x1", "x2"), b = s[2:3, 1], lo = s[2:3, 1] - 1.96 * s[2:3, 2],
             hi = s[2:3, 1] + 1.96 * s[2:3, 2])
}))
ci_ols <- confint(ols)
garis_ols <- data.frame(var = c("x1", "x2"), b = coef(ols)[2:3], lo = ci_ols[2:3, 1], hi = ci_ols[2:3, 2])
plot_proses <- function(v, w, judul) {
  ggplot(proses[proses$var == v, ], aes(tau, b)) +
    geom_rect(data = garis_ols[garis_ols$var == v, ], aes(xmin = -Inf, xmax = Inf, ymin = lo, ymax = hi),
              inherit.aes = FALSE, fill = pal[4], alpha = 0.12) +
    geom_hline(data = garis_ols[garis_ols$var == v, ], aes(yintercept = b), colour = pal[4], linetype = "dotted", linewidth = 0.9) +
    geom_hline(yintercept = 0, colour = warna$redup, linewidth = 0.4) +
    geom_ribbon(aes(ymin = lo, ymax = hi), fill = w, alpha = 0.22) +
    geom_line(colour = w, linewidth = 1.1) + geom_point(colour = w, size = 1.6) +
    labs(title = judul, x = "Quantile (tau)", y = "Coefficient")
}
p2 <- plot_proses("x1", pal[2], paste("Effect of", substr("work_experience", 1, 30)))
p3 <- plot_proses("x2", pal[5], paste("Effect of", substr("certifications", 1, 30)))
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)