Pengaruh PDRB dan jumlah pelanggan industri terhadap konsumsi listrik antar provinsi dengan regresi data panel

# Contoh kode R: Panel Data Regression (Pooled, Fixed and Random Effects)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(plm)
  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 entitas x tahun ----
n_ent <- sample(15:30, 1)                        # banyaknya entitas (Group)
n_thn <- sample(8:12, 1)                         # banyaknya tahun
d <- expand.grid(tahun = 2023 - n_thn + seq_len(n_thn), entitas = sprintf("Group %02d", seq_len(n_ent)),
                 stringsAsFactors = FALSE)[, c("entitas", "tahun")]
ai <- rnorm(n_ent)                                # efek individu (tak teramati)
a <- ai[match(d$entitas, sort(unique(d$entitas)))]
r_a <- runif(1, 0.35, 0.8)                       # korelasi efek individu dengan X1 (memicu FEM)
z1 <- r_a * a + sqrt(1 - r_a^2) * rnorm(nrow(d))
z2 <- 0.25 * a + rnorm(nrow(d))
tt <- (d$tahun - min(d$tahun)) / n_thn
b1 <- runif(1, 0.35, 0.7)
b2 <- runif(1, 0.15, 0.45) * sample(c(-1, 1), 1)
d$provincial_grdp <- round(50 + 10 * z1, 3)
d$industrial_customers <- round(30 + 6 * z2, 3)
d$electricity_consumption <- round(100 + 15 * (b1 * z1 + b2 * z2 + runif(1, 0.6, 1.1) * a + 0.3 * tt +
                                      rnorm(nrow(d), 0, runif(1, 0.3, 0.5))), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: CEM (pooled), FEM (within), REM (random) ----
pd <- pdata.frame(d, index = c("entitas", "tahun"))
f <- electricity_consumption ~ provincial_grdp + industrial_customers
cem <- plm(f, data = pd, model = "pooling")
fem <- plm(f, data = pd, model = "within")
rem <- tryCatch(plm(f, data = pd, model = "random"),
                error = function(e) plm(f, data = pd, model = "random", random.method = "walhus"))
print(summary(cem)); print(summary(fem)); print(summary(rem))
# ---- Pemilihan model ----
chow <- pFtest(fem, cem)                          # Chow: CEM vs FEM
haus <- phtest(fem, rem)                          # Hausman: FEM vs REM
lmbp <- plmtest(cem, effect = "individual", type = "bp")   # LM Breusch-Pagan: CEM vs REM
print(chow); print(haus); print(lmbp)
p_chow <- unname(chow$p.value); p_haus <- unname(haus$p.value); p_lm <- unname(lmbp$p.value)
terpilih <- if (p_chow < 0.05) { if (p_haus < 0.05) "FEM" else "REM" } else { if (p_lm < 0.05) "REM" else "CEM" }
mod <- list(CEM = cem, FEM = fem, REM = rem)[[terpilih]]
# galat baku robust (Arellano, klaster per entitas) untuk model terpilih
ct <- lmtest::coeftest(mod, vcov = vcovHC(mod, method = "arellano", cluster = "group"))
print(ct)
r2_sel <- unname(summary(mod)$r.squared["rsq"])


# ---- Grafik ----
# 1) perbandingan koefisien di tiga model (IK 95%)
lab_v <- make.unique(c("provincial_grdp", "industrial_customers"), sep = " ")
dk <- do.call(rbind, lapply(c("CEM", "FEM", "REM"), function(m) {
  s <- coef(summary(list(CEM = cem, FEM = fem, REM = rem)[[m]]))
  data.frame(model = m, var = lab_v, b = s[c("provincial_grdp", "industrial_customers"), 1], se = s[c("provincial_grdp", "industrial_customers"), 2])
}))
dk$var <- factor(dk$var, levels = lab_v)
dk$model <- factor(dk$model, levels = c("REM", "FEM", "CEM"))
p1 <- ggplot(dk, aes(b, model, colour = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_errorbarh(aes(xmin = b - 1.96 * se, xmax = b + 1.96 * se), height = 0.2, linewidth = 0.9) +
  geom_point(size = 3) +
  facet_wrap(~var, scales = "free_x", ncol = 1) +
  scale_colour_manual(values = c(CEM = pal[1], FEM = pal[2], REM = pal[3]), guide = "none") +
  labs(title = "Coefficients by panel model (95% CI)", x = "Estimate", y = NULL)
# 2) intersep efek tetap per entitas
fe <- fixef(fem)
dfe <- data.frame(ent = sub("^.* ", "", names(fe)), a = as.numeric(fe) - mean(fe))
dfe$ent <- factor(dfe$ent, levels = dfe$ent[order(dfe$a)])
p2 <- ggplot(dfe, aes(a, ent)) +
  geom_vline(xintercept = 0, colour = warna$redup) +
  geom_segment(aes(x = 0, xend = a, yend = ent), colour = pal[4], linewidth = 0.7) +
  geom_point(colour = pal[4], size = 2) +
  theme(axis.text.y = element_text(size = rel(0.6))) +
  labs(title = "Fixed-effect intercepts (vs mean)", x = "Deviation in electricity_consumption", y = "Group")
# 3) hubungan within: data dikurangi rata-rata entitas
dw <- data.frame(x = as.numeric(Within(pd$provincial_grdp)), y = as.numeric(Within(pd$electricity_consumption)))
p3 <- ggplot(dw, aes(x, y)) +
  geom_point(colour = pal[1], alpha = 0.45, size = 1.4) +
  geom_smooth(method = "lm", formula = y ~ x, colour = pal[2], fill = pal[2], alpha = 0.15, linewidth = 1) +
  labs(title = "Within-transformed relationship", x = "provincial_grdp (demeaned)", y = "electricity_consumption (demeaned)")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else (p1 | p2) / p3

print(p_hasil)