Evaluasi dampak program digitalisasi UMKM terhadap omzet usaha dengan Difference-in-Differences

# Contoh kode R: Difference-in-Differences
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  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: panel unit x 8 periode, kebijakan mulai periode 5 ----
n_unit <- max(60, round(200 / 3))
n_per <- 8; mulai <- 5
d <- expand.grid(periode = seq_len(n_per), entitas = sprintf("U%03d", seq_len(n_unit)))[, c("entitas", "periode")]
d$entitas <- as.character(d$entitas)
idx <- match(d$entitas, sprintf("U%03d", seq_len(n_unit)))
grp <- rbinom(n_unit, 1, runif(1, 0.4, 0.55))   # 1 = kelompok perlakuan
d$programme_participant <- grp[idx]
d$pasca <- as.integer(d$periode >= mulai)       # 1 = setelah kebijakan berlaku
ui <- rnorm(n_unit)[idx]                         # efek unit
tren <- cumsum(c(0, rnorm(n_per - 1, runif(1, 0.05, 0.15), 0.08)))[d$periode]   # tren waktu bersama
z2 <- 0.3 * ui + rnorm(nrow(d))
d$business_age <- round(30 + 6 * z2, 3)
att <- runif(1, 0.35, 0.8) * sample(c(-1, 1), 1, prob = c(0.35, 0.65))   # efek kebijakan (satuan SD)
rampa <- pmin(1, 0.6 + 0.15 * (d$periode - mulai))                          # efek menguat bertahap
d$monthly_revenue <- round(100 + 15 * (0.5 * ui + runif(1, 0.2, 0.5) * d$programme_participant + tren + 0.2 * z2 +
                                      att * d$programme_participant * d$pasca * rampa + rnorm(nrow(d), 0, 0.25)), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: DiD dua arah dengan galat baku klaster per unit ----
fit <- lm(monthly_revenue ~ programme_participant * pasca + business_age, data = d)
vc <- sandwich::vcovCL(fit, cluster = ~entitas)
ct <- lmtest::coeftest(fit, vcov = vc)
print(ct)
did <- ct["programme_participant:pasca", ]
ci_did <- did[1] + c(-1, 1) * qt(0.975, fit$df.residual) * did[2]
# ---- Uji tren paralel: beda kemiringan tren pra-kebijakan antar kelompok ----
pra <- subset(d, pasca == 0)
fit_pra <- lm(monthly_revenue ~ programme_participant * periode + business_age, data = pra)
ct_pra <- lmtest::coeftest(fit_pra, vcov = sandwich::vcovCL(fit_pra, cluster = ~entitas))
print(ct_pra)
p_tren <- ct_pra["programme_participant:periode", 4]
# ---- Event study: efek per periode relatif terhadap periode sebelum kebijakan (periode 4) ----
de <- d
per_lain <- setdiff(seq_len(n_per), mulai - 1)
for (k in per_lain) de[[paste0("ev", k)]] <- de$programme_participant * (de$periode == k)
f_es <- as.formula(paste("monthly_revenue ~ factor(periode) + programme_participant + business_age +", paste0("ev", per_lain, collapse = " + ")))
fit_es <- lm(f_es, data = de)
ct_es <- lmtest::coeftest(fit_es, vcov = sandwich::vcovCL(fit_es, cluster = ~entitas))
es <- data.frame(periode = c(per_lain, mulai - 1),
                 b = c(ct_es[paste0("ev", per_lain), 1], 0), se = c(ct_es[paste0("ev", per_lain), 2], 0))
print(es[order(es$periode), ])
r2 <- summary(fit)$r.squared


# ---- Grafik ----
d$kel <- factor(ifelse(d$programme_participant == 1, "Treated", "Control"), levels = c("Treated", "Control"))
# 1) rata-rata Y per periode per kelompok + kontrafaktual kelompok perlakuan
rt <- aggregate(monthly_revenue ~ periode + kel, data = d, FUN = mean)
mt <- rt[rt$kel == "Treated", ]; mc <- rt[rt$kel == "Control", ]
celah <- mean(mt$monthly_revenue[mt$periode < mulai] - mc$monthly_revenue[mc$periode < mulai])
kf <- data.frame(periode = mc$periode[mc$periode >= mulai - 1], y = mc$monthly_revenue[mc$periode >= mulai - 1] + celah)
p1 <- ggplot(rt, aes(periode, monthly_revenue, colour = kel)) +
  annotate("rect", xmin = mulai - 0.5, xmax = n_per + 0.3, ymin = -Inf, ymax = Inf, fill = pal[5], alpha = 0.08) +
  geom_vline(xintercept = mulai - 0.5, linetype = "dashed", colour = warna$redup) +
  geom_line(data = kf, aes(periode, y), colour = pal[1], linetype = "dotted", linewidth = 0.9, inherit.aes = FALSE) +
  geom_line(linewidth = 1.1) + geom_point(size = 2.2) +
  annotate("text", x = mulai - 0.4, y = Inf, label = "Policy starts", hjust = 0, vjust = 1.5, size = 3,
           colour = warna$teks) +
  scale_colour_manual(values = c(Treated = pal[1], Control = pal[2]), name = NULL) +
  scale_x_continuous(breaks = seq_len(n_per)) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
  labs(title = "Group means by period (dotted = counterfactual)", x = "Period", y = "monthly_revenue")
# 2) event study dengan IK 95%
p2 <- ggplot(es, aes(periode, b)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_vline(xintercept = mulai - 0.5, linetype = "dashed", colour = warna$redup) +
  geom_ribbon(aes(ymin = b - 1.96 * se, ymax = b + 1.96 * se), fill = pal[3], alpha = 0.15) +
  geom_errorbar(aes(ymin = b - 1.96 * se, ymax = b + 1.96 * se), width = 0.15, colour = pal[3], linewidth = 0.8) +
  geom_point(size = 2.6, colour = pal[3]) +
  scale_x_continuous(breaks = seq_len(n_per)) +
  labs(title = sprintf("Event study (ref. period %d, 95%% CI)", mulai - 1), x = "Period", y = "Effect on monthly_revenue")
# 3) tabel 2x2 rata-rata sebelum/sesudah
r22 <- aggregate(monthly_revenue ~ kel + pasca, data = d, FUN = mean)
r22$waktu <- factor(ifelse(r22$pasca == 1, "After", "Before"), levels = c("Before", "After"))
p3 <- ggplot(r22, aes(waktu, monthly_revenue, colour = kel, group = kel)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_text(aes(label = sprintf("%.1f", monthly_revenue)), vjust = -1, size = 3, colour = warna$teks, show.legend = FALSE) +
  scale_colour_manual(values = c(Treated = pal[1], Control = pal[2]), name = NULL) +
  scale_y_continuous(expand = expansion(mult = c(0.12, 0.22))) +
  labs(title = sprintf("2x2 means (DiD = %.2f)", did[1]), x = NULL, y = "Mean monthly_revenue")
p_hasil <- if (tata == 1) p1 / (p2 | p3) + plot_layout(heights = c(1.2, 1)) else p1 | (p2 / p3)

print(p_hasil)