Household tutoring expenditure: Tobit and Heckman selection models

# 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(parental_income = round(rnorm(n, 50, 10), 3),
                child_grade_level = round(rnorm(n, 30, 6), 3),
                akses = round(rnorm(n, 50, 10), 1))      # variabel eksklusi: hanya memengaruhi partisipasi
z1 <- as.numeric(scale(d$parental_income))
z2 <- as.numeric(scale(d$child_grade_level))
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$tutoring_expenditure <- 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(tutoring_expenditure ~ parental_income + child_grade_level, data = d)
tob <- tobit(tutoring_expenditure ~ parental_income + child_grade_level, 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 ~ parental_income + child_grade_level + akses,
                outcome = tutoring_expenditure ~ parental_income + child_grade_level, 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)", "parental_income", "child_grade_level")
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$tutoring_expenditure == 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$tutoring_expenditure == 0, "Zero (censored / not observed)", "Positive"),
                  levels = c("Zero (censored / not observed)", "Positive"))
p1 <- ggplot(d, aes(tutoring_expenditure, 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 tutoring_expenditure (%.0f%% zeros)", 100 * mean(d$tutoring_expenditure == 0)), x = "tutoring_expenditure", y = "Count")
# 2) koefisien parental_income: 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)["parental_income"], st["parental_income", 1], b_hk["parental_income", 1]),
                 se = c(coef(summary(ols))["parental_income", 2], st["parental_income", 2], b_hk["parental_income", 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 parental_income (95% CI)", x = "Estimate", y = NULL)
# 3) prediksi E[Y] menurut parental_income (X2 di rata-rata)
xs <- seq(min(d$parental_income), max(d$parental_income), length.out = 100)
nd <- data.frame(parental_income = xs, child_grade_level = mean(d$child_grade_level))
xb <- as.numeric(model.matrix(~ parental_income + child_grade_level, 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(~ parental_income + child_grade_level, nd) %*% b_hk[, 1]), model = "Heckman E[Y*]"))
p3 <- ggplot() +
  geom_point(data = d, aes(parental_income, tutoring_expenditure), 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 tutoring_expenditure by parental_income", x = "parental_income", y = "tutoring_expenditure")
p_hasil <- if (tata == 1) p1 / (p2 | p3) else p1 | (p2 / p3)

print(p_hasil)