Efek pelatihan pemandu wisata bersertifikat terhadap pendapatan pemandu di Borobudur

# Contoh kode R: Propensity Score Matching
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(MatchIt)
  library(sandwich)
  library(lmtest)
})
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: partisipasi program tidak acak (bergantung kovariat) ----
n <- max(300, 200)
d <- data.frame(guiding_experience = round(rnorm(n, 30, 6), 3),
                skor_awal = round(rnorm(n, 50, 10), 1),        # kovariat: skor kondisi awal
                wilayah_kota = rbinom(n, 1, 0.45))             # kovariat: 1 = perkotaan
z2 <- as.numeric(scale(d$guiding_experience))
zs <- as.numeric(scale(d$skor_awal))
# seleksi: peluang ikut program bergantung kovariat
lp <- -0.6 + runif(1, 0.5, 0.9) * z2 + runif(1, 0.3, 0.6) * zs + runif(1, 0.3, 0.8) * d$wilayah_kota
d$attended_training <- rbinom(n, 1, plogis(lp))
att <- runif(1, 0.25, 0.6) * sample(c(-1, 1), 1, prob = c(0.3, 0.7))   # efek program (satuan SD)
d$guide_income <- round(100 + 15 * (runif(1, 0.4, 0.7) * z2 + runif(1, 0.3, 0.6) * zs + 0.25 * d$wilayah_kota +
                                      att * d$attended_training + rnorm(n, 0, 0.6)), 3)
d <- d[, c("guide_income", "attended_training", "guiding_experience", "skor_awal", "wilayah_kota")]
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: Propensity Score Matching (nearest neighbour 1:1, logit, caliper 0.1 SD) ----
naif <- t.test(guide_income ~ attended_training, data = d)                  # beda rata-rata tanpa penyesuaian
beda_naif <- unname(diff(naif$estimate))             # perlakuan - kontrol
ci_naif <- -rev(naif$conf.int)
m <- matchit(attended_training ~ guiding_experience + skor_awal + wilayah_kota, data = d, method = "nearest",
             distance = "glm", link = "logit", caliper = 0.1, replace = FALSE,
             exact = ~wilayah_kota)                         # wilayah dipasangkan persis
print(m)
sm <- summary(m, standardize = TRUE)
print(sm)
md <- match.data(m)
# ATT: regresi pada data hasil matching (bobot matching) + galat baku robust klaster per pasangan
fit_att <- lm(guide_income ~ attended_training + guiding_experience + skor_awal + wilayah_kota, data = md, weights = weights)
ct <- lmtest::coeftest(fit_att, vcov = sandwich::vcovCL(fit_att, cluster = ~subclass))
print(ct)
att_hat <- ct["attended_training", 1]; se_att <- ct["attended_training", 2]; p_att <- ct["attended_training", 4]
ci_att <- att_hat + c(-1, 1) * 1.96 * se_att
smd_sebelum <- sm$sum.all[, "Std. Mean Diff."]
smd_sesudah <- sm$sum.matched[, "Std. Mean Diff."]
n_pasang <- sum(md$attended_training == 1)
cat("Treated matched:", n_pasang, "dari", sum(d$attended_training == 1), "\n")


# ---- Grafik ----
# 1) love plot: |SMD| sebelum vs sesudah matching
lab_kov <- c("Propensity score", "guiding_experience", "Baseline score", "Urban area")
dl <- data.frame(kov = factor(rep(lab_kov, 2), levels = rev(lab_kov)),
                 smd = abs(c(smd_sebelum, smd_sesudah)),
                 tahap = factor(rep(c("Before", "After"), each = length(lab_kov)), levels = c("Before", "After")))
p1 <- ggplot(dl, aes(smd, kov, colour = tahap)) +
  geom_vline(xintercept = 0.1, linetype = "dashed", colour = warna$redup) +
  geom_line(aes(group = kov), colour = warna$redup, linewidth = 0.5) +
  geom_point(size = 3) +
  scale_colour_manual(values = c(Before = pal[2], After = pal[1]), name = NULL) +
  labs(title = "Covariate balance (love plot)", x = "Absolute standardized mean difference", y = NULL)
# 2) distribusi propensity score sebelum dan sesudah matching
d$ps <- m$distance
dps <- rbind(data.frame(ps = d$ps, kel = d$attended_training, tahap = "Before"),
             data.frame(ps = md$distance, kel = md$attended_training, tahap = "After"))
dps$tahap <- factor(dps$tahap, levels = c("Before", "After"))
dps$kel <- factor(ifelse(dps$kel == 1, "Treated", "Control"), levels = c("Treated", "Control"))
p2 <- ggplot(dps, aes(ps, fill = kel, colour = kel)) +
  geom_density(alpha = 0.3, linewidth = 0.7) +
  facet_wrap(~tahap, ncol = 2) +
  scale_fill_manual(values = c(Treated = pal[1], Control = pal[3]), name = NULL) +
  scale_colour_manual(values = c(Treated = pal[1], Control = pal[3]), name = NULL) +
  labs(title = "Propensity score overlap", x = "Propensity score", y = "Density")
# 3) estimasi efek: naif vs ATT hasil matching
de <- data.frame(est = factor(c("Naive difference", "Matched ATT"), levels = c("Matched ATT", "Naive difference")),
                 b = c(beda_naif, att_hat), lo = c(ci_naif[1], ci_att[1]), hi = c(ci_naif[2], ci_att[2]))
de$teks <- sprintf("%.2f", de$b)
p3 <- ggplot(de, aes(b, est, colour = est)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.15, linewidth = 1) +
  geom_point(size = 3.4) +
  geom_text(aes(label = teks), vjust = -1.2, size = 3.2, colour = warna$teks) +
  scale_colour_manual(values = c("Naive difference" = pal[2], "Matched ATT" = pal[1]), guide = "none") +
  scale_y_discrete(expand = expansion(add = 0.7)) +
  labs(title = "Effect on guide_income (95% CI)", x = "Difference in guide_income", y = NULL)
p_hasil <- if (tata == 1) (p1 | p3) / p2 else p1 | (p2 / p3)

print(p_hasil)