Determinan jam lembur karyawan dengan banyak nilai nol: model Tobit

# Contoh kode R: Tobit and Heckman Selection Models
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(AER)
  library(sampleSelection)
})
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 bernilai nol untuk yang tidak berpartisipasi ----
n <- max(300, 200)
d <- data.frame(work_target_load = round(rnorm(n, 50, 10), 3),
                hourly_wage = round(rnorm(n, 30, 6), 3),
                akses = round(rnorm(n, 50, 10), 1))      # variabel eksklusi: hanya memengaruhi partisipasi
z1 <- as.numeric(scale(d$work_target_load))
z2 <- as.numeric(scale(d$hourly_wage))
za <- as.numeric(scale(d$akses))
rho <- runif(1, 0.3, 0.7)                                # korelasi galat seleksi dan galat outcome
u <- rnorm(n)
e <- rho * u + sqrt(1 - rho^2) * rnorm(n)
# persamaan seleksi (partisipasi)
g0 <- qnorm(runif(1, 0.55, 0.72))
d$partisipasi <- as.integer(g0 + runif(1, 0.4, 0.8) * z1 + runif(1, 0.1, 0.3) * z2 + runif(1, 0.6, 1) * za + u > 0)
# persamaan outcome (hanya teramati bila partisipasi = 1)
b1 <- runif(1, 0.25, 0.45); b2 <- runif(1, 0.08, 0.25) * sample(c(-1, 1), 1)
y_star <- 100 * (1 + b1 * z1 + b2 * z2 + 0.3 * e)
d$monthly_overtime_hours <- round(ifelse(d$partisipasi == 1, pmax(y_star, 0), 0), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis 1: OLS naif dan Tobit (tersensor kiri di 0) ----
ols <- lm(monthly_overtime_hours ~ work_target_load + hourly_wage, data = d)
tob <- tobit(monthly_overtime_hours ~ work_target_load + hourly_wage, left = 0, data = d)
print(summary(ols)); print(summary(tob))
st <- coef(summary(tob))
sig_tob <- tob$scale
# ---- Analisis 2: Heckman dua tahap (heckit) dengan variabel eksklusi akses ----
hk <- selection(selection = partisipasi ~ work_target_load + hourly_wage + akses,
                outcome = monthly_overtime_hours ~ work_target_load + hourly_wage, data = d, method = "2step")
print(summary(hk))
sh <- coef(summary(hk))
io <- hk$param$index$betaO                   # baris persamaan outcome
b_hk <- sh[io, ]; rownames(b_hk) <- c("(Intercept)", "work_target_load", "hourly_wage")
imr <- sh[hk$param$index$Mills, ]           # koefisien inverse Mills ratio (lambda)
rho_hat <- unname(coef(hk)["rho"])
p_imr <- unname(imr[4])
nol <- 100 * mean(d$monthly_overtime_hours == 0)
cat(sprintf("Proporsi nol: %.1f%%, rho = %.3f, p IMR = %.4f\n", nol, rho_hat, p_imr))


# ---- Grafik ----
# 1) histogram Y dengan massa di nol
d$jenis <- factor(ifelse(d$monthly_overtime_hours == 0, "Zero (censored / not observed)", "Positive"),
                  levels = c("Zero (censored / not observed)", "Positive"))
p1 <- ggplot(d, aes(monthly_overtime_hours, fill = jenis)) +
  geom_histogram(bins = 30, colour = warna$panel, alpha = 0.9) +
  scale_fill_manual(values = c(pal[2], pal[1]), name = NULL) +
  labs(title = sprintf("Distribution of monthly_overtime_hours (%.0f%% zeros)", 100 * mean(d$monthly_overtime_hours == 0)), x = "monthly_overtime_hours", y = "Count")
# 2) koefisien work_target_load: OLS vs Tobit vs Heckman (IK 95%)
dk <- data.frame(model = factor(c("OLS", "Tobit", "Heckman"), levels = c("Heckman", "Tobit", "OLS")),
                 b = c(coef(ols)["work_target_load"], st["work_target_load", 1], b_hk["work_target_load", 1]),
                 se = c(coef(summary(ols))["work_target_load", 2], st["work_target_load", 2], b_hk["work_target_load", 2]))
dk$teks <- sprintf("%.3f", dk$b)
p2 <- ggplot(dk, aes(b, model, colour = model)) +
  geom_errorbarh(aes(xmin = b - 1.96 * se, xmax = b + 1.96 * se), height = 0.18, linewidth = 1) +
  geom_point(size = 3.2) +
  geom_text(aes(label = teks), vjust = -1.2, size = 3, colour = warna$teks) +
  scale_colour_manual(values = c(OLS = pal[2], Tobit = pal[1], Heckman = pal[3]), guide = "none") +
  scale_y_discrete(expand = expansion(add = 0.7)) +
  labs(title = "Coefficient of work_target_load (95% CI)", x = "Estimate", y = NULL)
# 3) prediksi E[Y] menurut work_target_load (X2 di rata-rata)
xs <- seq(min(d$work_target_load), max(d$work_target_load), length.out = 100)
nd <- data.frame(work_target_load = xs, hourly_wage = mean(d$hourly_wage))
xb <- as.numeric(model.matrix(~ work_target_load + hourly_wage, nd) %*% coef(tob))
dp <- rbind(data.frame(x = xs, y = predict(ols, nd), model = "OLS"),
            data.frame(x = xs, y = pnorm(xb / sig_tob) * xb + sig_tob * dnorm(xb / sig_tob), model = "Tobit E[Y]"),
            data.frame(x = xs, y = as.numeric(model.matrix(~ work_target_load + hourly_wage, nd) %*% b_hk[, 1]), model = "Heckman E[Y*]"))
p3 <- ggplot() +
  geom_point(data = d, aes(work_target_load, monthly_overtime_hours), colour = warna$redup, alpha = 0.3, size = 1.1) +
  geom_line(data = dp, aes(x, y, colour = model), linewidth = 1.1) +
  scale_colour_manual(values = c("OLS" = pal[2], "Tobit E[Y]" = pal[1], "Heckman E[Y*]" = pal[3]), name = NULL) +
  labs(title = "Predicted monthly_overtime_hours by work_target_load", x = "work_target_load", y = "monthly_overtime_hours")
p_hasil <- if (tata == 1) p1 / (p2 | p3) else p1 | (p2 / p3)

print(p_hasil)