# ══════════════════════════════════════════════════════════════════════
# ⚠  SESUAIKAN KEDUA PATH INI
# ══════════════════════════════════════════════════════════════════════
PATH_EXCEL <- "C:/Users/User/Documents/AKA/Projek Paper/Pustaka 2026/Data Olah Pusaka - 2025.xlsx"
PATH_SHP   <- "C:/Users/User/Documents/AKA/Projek Paper/Pustaka 2026/SHP Jateng/Kab_Jateng_proj.shp"
# ══════════════════════════════════════════════════════════════════════

data_raw <- read_excel(PATH_EXCEL) |>
  rename(jml_penduduk = `jumlah/penduduk`)

data_pusaka <- data_raw |>
  filter(kodekab >= 3301 & kodekab <= 3376) |>
  mutate(
    # kdkab: 2 digit terakhir kode BPS; diawali 7x = kota
    kdkab2  = kodekab %% 100,
    is_kota = (floor(kodekab %% 100 / 10) == 7) |
              (kodekab >= 3371 & kodekab <= 3376),
    # Label peta: "Kab. Cilacap" atau "Kota Semarang"
    label_peta = paste0(
      ifelse(is_kota, "Kota ", "Kab. "),
      str_to_title(nmkab)
    ),
    dmu_id = label_peta   # pakai label_peta sebagai DMU identifier
  )

cat("Observasi :", nrow(data_pusaka), "\n")

cat("Kota      :", sum(data_pusaka$is_kota), "\n")

cat("Kabupaten :", sum(!data_pusaka$is_kota), "\n")

cat("Cek NA DEA:\n")

data_pusaka |> select(x1_TFI, x2_TII, x3_TRI_sum, y1_TNLI_sum, y2_TSI) |>
  summarise(across(everything(), ~sum(is.na(.)))) |> print()

shp_jateng <- st_read(PATH_SHP, quiet = TRUE) |>
  filter(kdkab != "00") |>
  mutate(kodekab = as.numeric(kodekab)) |>
  st_transform(crs = 4326)

shp_merged <- shp_jateng |>
  left_join(data_pusaka, by = "kodekab", suffix = c("_shp", "_data")) |>
  rename(nmkab = nmkab_data) |>
  select(-nmkab_shp)

cat("Baris shapefile :", nrow(shp_merged), "\n")

cat("Cek NA x1_TFI   :", sum(is.na(shp_merged$x1_TFI)), "\n")

# ── Helper: centroid + label_peta untuk semua peta ──────────────────
get_centroid_labels <- function(shp) {
  shp |>
    st_centroid() |>
    mutate(
      lon = st_coordinates(geometry)[, 1],
      lat = st_coordinates(geometry)[, 2]
    ) |>
    as.data.frame() |>
    select(label_peta, lon, lat)
}

centroid_labels <- get_centroid_labels(shp_merged)

data_pusaka |>
  select(
    `TFI (Kunjungan)`       = x1_TFI,
    `TII (Infrastruktur)`   = x2_TII,
    `TRI (POI Wisata)`      = x3_TRI_sum,
    `NTL (Nighttime Light)` = y1_TNLI_sum,
    `TSI (Rating)`          = y2_TSI,
    `Kepadatan Pddk`        = z_density,
    `PDRB Non-Wisata`       = z_pdrb_nonwisata,
    `Digitalisasi`          = reg_dig,
    `Aksesibilitas`         = reg_akses,
    `Elevasi (mdpl)`        = reg_elevasi,
    `Curah Hujan (mm)`      = x_hujan,
    `PDRB Tersier`          = X_pdrb3
  ) |>
  pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai") |>
  group_by(Variabel) |>
  summarise(
    N    = n(),
    Min  = min(Nilai, na.rm = TRUE),
    Mean = mean(Nilai, na.rm = TRUE),
    Max  = max(Nilai, na.rm = TRUE),
    SD   = sd(Nilai, na.rm = TRUE),
    .groups = "drop"
  ) |>
  kable(digits = 3,
        caption = "Tabel 1. Statistik Deskriptif Variabel Penelitian (n = 35 kab/kota)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

df_dea <- data_pusaka |>
  select(dmu_id, x1_TFI, x2_TII, x3_TRI_sum, y1_TNLI_sum, y2_TSI) |>
  as.data.frame()

deadata_s1 <- make_deadata(
  datadea = df_dea, dmus = "dmu_id",
  inputs  = c("x1_TFI", "x2_TII", "x3_TRI_sum"),
  outputs = c("y1_TNLI_sum", "y2_TSI")
)

model_s1 <- model_sbmeff(deadata_s1, orientation = "oo", rts = "vrs")

# deaR mengembalikan θ_raw ∈ (0,1]. Konversi ke FEO = 1/θ_raw (≥1).
# FEO = 1 → efisien sempurna | FEO > 1 → butuh ekspansi output
eff_s1_raw  <- sapply(model_s1$DMU, function(x) x$efficiency)
eff_s1_norm <- 1 / eff_s1_raw   # FEO Stage 1

slack_input_s1  <- do.call(rbind, lapply(model_s1$DMU, function(x) x$slack_input))
slack_output_s1 <- do.call(rbind, lapply(model_s1$DMU, function(x) x$slack_output))
rownames(slack_input_s1)  <- data_pusaka$dmu_id
rownames(slack_output_s1) <- data_pusaka$dmu_id
colnames(slack_input_s1)  <- c("slack_TFI", "slack_TII", "slack_TRI")
colnames(slack_output_s1) <- c("slack_NTL", "slack_TSI")

cat("Stage 1 — FEO rata-rata (1/θ_raw):", round(mean(eff_s1_norm, na.rm = TRUE), 4), "\n")

cat("  → Artinya rata-rata output perlu ×",
    round(mean(eff_s1_norm, na.rm = TRUE), 2), "untuk mencapai frontier\n")

cat("DMU efisien (FEO = 1)  :", sum(eff_s1_norm <= 1.0001, na.rm = TRUE), "\n")

cat("Rentang FEO            :", round(min(eff_s1_norm, na.rm = TRUE), 4),
    "–", round(max(eff_s1_norm, na.rm = TRUE), 4), "\n")

df_stage2 <- data.frame(
  dmu_id           = data_pusaka$dmu_id,
  z_density        = data_pusaka$z_density,
  z_pdrb_nonwisata = data_pusaka$z_pdrb_nonwisata,
  x1_TFI           = data_pusaka$x1_TFI,
  x2_TII           = data_pusaka$x2_TII,
  x3_TRI_sum       = data_pusaka$x3_TRI_sum,
  slack_TFI        = slack_input_s1[, "slack_TFI"],
  slack_TII        = slack_input_s1[, "slack_TII"],
  slack_TRI        = slack_input_s1[, "slack_TRI"]
)

adjust_input_sfa <- function(slack_vec, z1, z2, x_vec, var_name) {
  df_sfa <- data.frame(slack = slack_vec, z1 = z1, z2 = z2)
  result <- tryCatch({
    fit_sfa <- sfa(slack ~ z1 + z2, data = df_sfa, ineffDecrease = FALSE)
    b       <- coef(fit_sfa)[c("(Intercept)", "z1", "z2")]
    env_hat <- as.numeric(b[1] + b[2] * z1 + b[3] * z2)
    v_hat   <- as.numeric(residuals(lm(slack ~ z1 + z2, data = df_sfa)))
    cat("  [SFA OK]:", var_name, "\n")
    list(env_hat = env_hat, v_hat = v_hat, metode = "SFA")
  }, error = function(e) {
    cat("  [OLS fallback]:", var_name, "\n")
    fit_ols <- lm(slack ~ z1 + z2, data = df_sfa)
    list(env_hat = as.numeric(fitted(fit_ols)),
         v_hat   = as.numeric(residuals(fit_ols)),
         metode  = "OLS fallback")
  })
  x_adj <- pmax(x_vec +
    (max(result$env_hat) - result$env_hat) +
    (max(result$v_hat)   - result$v_hat), 1e-6)
  list(x_adj = x_adj, metode = result$metode)
}

cat("Purifikasi Stage 2:\n")

adj_TFI <- adjust_input_sfa(df_stage2$slack_TFI, df_stage2$z_density,
                             df_stage2$z_pdrb_nonwisata, df_stage2$x1_TFI, "TFI")

adj_TII <- adjust_input_sfa(df_stage2$slack_TII, df_stage2$z_density,
                             df_stage2$z_pdrb_nonwisata, df_stage2$x2_TII, "TII")

adj_TRI <- adjust_input_sfa(df_stage2$slack_TRI, df_stage2$z_density,
                             df_stage2$z_pdrb_nonwisata, df_stage2$x3_TRI_sum, "TRI")

df_adj <- data.frame(
  dmu_id     = df_stage2$dmu_id,
  x1_TFI_adj = adj_TFI$x_adj,
  x2_TII_adj = adj_TII$x_adj,
  x3_TRI_adj = adj_TRI$x_adj
)
cat("Metode — TFI:", adj_TFI$metode, "| TII:", adj_TII$metode,
    "| TRI:", adj_TRI$metode, "\n")

df_dea_s3 <- data.frame(
  dmu_id      = df_adj$dmu_id,
  x1_TFI_adj  = df_adj$x1_TFI_adj,
  x2_TII_adj  = df_adj$x2_TII_adj,
  x3_TRI_adj  = df_adj$x3_TRI_adj,
  y1_TNLI_sum = data_pusaka$y1_TNLI_sum,
  y2_TSI      = data_pusaka$y2_TSI
)

deadata_s3 <- make_deadata(
  datadea = df_dea_s3, dmus = "dmu_id",
  inputs  = c("x1_TFI_adj", "x2_TII_adj", "x3_TRI_adj"),
  outputs = c("y1_TNLI_sum", "y2_TSI")
)

model_s3 <- model_sbmeff(deadata_s3, orientation = "oo", rts = "vrs")

# FEO Stage 3 = efisiensi murni, bebas bias urban (setelah purifikasi SFA)
eff_s3_raw  <- sapply(model_s3$DMU, function(x) x$efficiency)
eff_s3_norm <- 1 / eff_s3_raw   # FEO Stage 3 (≥1; semakin kecil = semakin efisien)

slack_input_s3  <- do.call(rbind, lapply(model_s3$DMU, function(x) x$slack_input))
slack_output_s3 <- do.call(rbind, lapply(model_s3$DMU, function(x) x$slack_output))
rownames(slack_input_s3)  <- data_pusaka$dmu_id
rownames(slack_output_s3) <- data_pusaka$dmu_id
colnames(slack_input_s3)  <- c("slack_TFI", "slack_TII", "slack_TRI")
colnames(slack_output_s3) <- c("slack_NTL", "slack_TSI")

cat("=== Stage 3 (Efisiensi Murni) ===\n")

cat("FEO₃ rata-rata (DMU inefisien) :", round(mean(eff_s3_norm[eff_s3_norm > 1.0001],
    na.rm = TRUE), 4), "\n")

cat("  → Rata-rata butuh ekspansi output ×",
    round(mean(eff_s3_norm[eff_s3_norm > 1.0001], na.rm = TRUE), 2),
    "(", round((mean(eff_s3_norm[eff_s3_norm > 1.0001], na.rm = TRUE)-1)*100, 1), "% lebih tinggi)\n")

cat("DMU efisien (FEO = 1)   :", sum(eff_s3_norm <= 1.0001, na.rm = TRUE), "\n")

cat("Rentang FEO₃            :", round(min(eff_s3_norm, na.rm = TRUE), 4),
    "–", round(max(eff_s3_norm, na.rm = TRUE), 4), "\n")

normalize <- function(x) {
  r <- range(x, na.rm = TRUE)
  if (diff(r) == 0) return(rep(0, length(x)))
  (x - r[1]) / diff(r)
}

df_slack <- data.frame(
  kodekab   = data_pusaka$kodekab,
  label_peta = data_pusaka$label_peta,
  nmkab     = data_pusaka$nmkab,
  s_TFI = as.numeric(slack_input_s3[, "slack_TFI"]),
  s_TII = as.numeric(slack_input_s3[, "slack_TII"]),
  s_TRI = as.numeric(slack_input_s3[, "slack_TRI"]),
  s_NTL = as.numeric(slack_output_s3[, "slack_NTL"]),
  s_TSI = as.numeric(slack_output_s3[, "slack_TSI"])
) |>
  mutate(
    Fenomena = case_when(
      s_TII > quantile(s_TII, 0.75, na.rm = TRUE) &
        s_NTL > quantile(s_NTL, 0.5,  na.rm = TRUE) ~ "Ghost Hotel Syndrome",
      s_TFI > quantile(s_TFI, 0.75, na.rm = TRUE) &
        s_TSI > quantile(s_TSI, 0.5,  na.rm = TRUE) ~ "Mass Tourism Tanpa Kualitas",
      s_TRI > quantile(s_TRI, 0.75, na.rm = TRUE) &
        s_NTL > quantile(s_NTL, 0.5,  na.rm = TRUE) ~ "Connectivity Gap",
      s_TSI > quantile(s_TSI, 0.75, na.rm = TRUE)   ~ "Kualitas Layanan Rendah",
      s_NTL > quantile(s_NTL, 0.75, na.rm = TRUE)   ~ "Aktivitas Malam Rendah",
      TRUE                                            ~ "Relatif Efisien"
    ),
    Prioritas = case_when(
      Fenomena == "Ghost Hotel Syndrome"        ~ "Aktivasi night economy & festival malam",
      Fenomena == "Mass Tourism Tanpa Kualitas" ~ "Pengendalian kuantitas, paket multi-hari",
      Fenomena == "Connectivity Gap"            ~ "Evaluasi konektivitas antar POI",
      Fenomena == "Kualitas Layanan Rendah"     ~ "Peningkatan kualitas hospitality",
      Fenomena == "Aktivitas Malam Rendah"      ~ "Evaluasi produk malam & keamanan",
      TRUE                                      ~ "Pertahankan & kembangkan sirkuit wisata"
    )
  )

df_slack |>
  arrange(desc(s_NTL + s_TSI)) |>
  select(label_peta, s_TFI, s_TII, s_TRI, s_NTL, s_TSI, Fenomena) |>
  mutate(across(s_TFI:s_TSI, ~round(., 3))) |>
  kable(
    caption   = "Tabel 4. Matriks Diagnostik Slack per Kab/Kota (Stage 3)",
    col.names = c("Kab/Kota", "s_TFI", "s_TII", "s_TRI", "s_NTL", "s_TSI", "Diagnosis")
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE) |>
  add_header_above(c(" " = 1,
                     "Input Excess" = 3, "Output Shortfall" = 2, " " = 1))

df_slack |>
  select(label_peta, Fenomena, Prioritas) |>
  arrange(Fenomena, label_peta) |>
  kable(
    caption   = "Tabel 5. Rekomendasi Kebijakan Berdasarkan Pola Slack",
    col.names = c("Kab/Kota", "Fenomena Terdeteksi", "Prioritas Kebijakan")
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

nb_queen <- poly2nb(shp_merged, queen = TRUE)
listw    <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
cat("Ringkasan matriks bobot:\n"); summary(nb_queen)

# PENTING: Negasi FEO agar interpretasi LISA konvensional berlaku.
# FEO tinggi = inefisiensi tinggi. Dengan negasi (-FEO):
# nilai tinggi = efisiensi tinggi → High-High = Hot-Spot efisiensi (konvensional)
eff_s3_spasial <- -as.numeric(shp_merged$Theta_S3)   # negasi FEO
if (any(is.na(eff_s3_spasial))) {
  eff_s3_spasial[is.na(eff_s3_spasial)] <- median(eff_s3_spasial, na.rm = TRUE)
}

moran_global <- moran.test(eff_s3_spasial, listw, zero.policy = TRUE)
print(moran_global)

data.frame(
  Statistik = c("Moran's I", "E(I)", "Varians", "p-value"),
  Nilai     = c(
    round(moran_global$estimate["Moran I statistic"], 4),
    round(moran_global$estimate["Expectation"],       4),
    round(moran_global$estimate["Variance"],          6),
    format.pval(moran_global$p.value, digits = 4, eps = 0.001)
  )
) |>
  kable(caption = "Tabel 6. Uji Moran's I Global — Efisiensi Murni (berbasis FEO₃, dinegasi)") |>
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) |>
  footnote(general = "Variabel yang diuji = -FEO₃ sehingga nilai tinggi = daerah efisien. High-High = klaster daerah efisien (Hot-Spot efisiensi).")

lisa_result <- localmoran(eff_s3_spasial, listw, zero.policy = TRUE)

z_score <- scale(eff_s3_spasial)[, 1]
lag_z   <- scale(lag.listw(listw, eff_s3_spasial, zero.policy = TRUE))[, 1]
p_val   <- lisa_result[, "Pr(z != E(Ii))"]

klaster_lisa <- case_when(
  p_val > 0.05                   ~ "Not Significant",
  z_score >= 0 & lag_z >= 0     ~ "High-High (Hot-Spot)",
  z_score <  0 & lag_z <  0     ~ "Low-Low (Cold-Spot)",
  z_score >= 0 & lag_z <  0     ~ "High-Low (Outlier+)",
  z_score <  0 & lag_z >= 0     ~ "Low-High (Outlier−)",
  TRUE                           ~ "Not Significant"
)

shp_merged$Ii_local <- lisa_result[, "Ii"]
shp_merged$p_LISA   <- p_val
shp_merged$klaster  <- klaster_lisa

shp_merged |>
  as.data.frame() |>
  select(label_peta, Theta_S3, Ii_local, p_LISA, klaster) |>
  arrange(klaster, Theta_S3) |>   # ascending: daerah efisien (FEO kecil) duluan
  mutate(across(c(Theta_S3, Ii_local, p_LISA), ~round(., 4))) |>
  kable(
    caption   = "Tabel 7. Hasil LISA per Kab/Kota Jawa Tengah",
    col.names = c("Kab/Kota", "FEO₃", "Ii Lokal", "p-value", "Klaster LISA")
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE) |>
  footnote(general = "FEO₃ = 1/θ_raw; nilai kecil = efisien. Klaster LISA dihitung dari -FEO₃ sehingga High-High = Hot-Spot daerah efisien.")

df_reg <- shp_merged |>
  as.data.frame() |>
  dplyr::select(
    kodekab,
    eff_s3      = Theta_S3,
    reg_dig, reg_akses, reg_elevasi, x_hujan
  ) |>
  mutate(eff_s3 = ifelse(is.na(eff_s3), median(eff_s3, na.rm = TRUE), eff_s3))

formula_ols <- eff_s3 ~ reg_dig + reg_akses + reg_elevasi + x_hujan
ols_model   <- lm(formula_ols, data = df_reg)

summary(ols_model)

bp_test <- bptest(ols_model)
kb_test <- bptest(ols_model, studentize = TRUE)
gq_test <- gqtest(ols_model, order.by = ~ reg_elevasi, data = df_reg)

data.frame(
  Uji       = c("Breusch-Pagan", "Koenker-Bassett", "Goldfeld-Quandt"),
  Statistik = round(c(bp_test$statistic, kb_test$statistic, gq_test$statistic), 4),
  df = c(as.character(bp_test$parameter), as.character(kb_test$parameter),
         paste(gq_test$parameter, collapse = ", ")),
  p_value = c(format.pval(bp_test$p.value, digits = 4, eps = 0.001),
              format.pval(kb_test$p.value, digits = 4, eps = 0.001),
              format.pval(gq_test$p.value, digits = 4, eps = 0.001)),
  Kesimpulan = c(
    ifelse(bp_test$p.value < 0.05, "Heteroskedastik ✓", "Homoskedastik"),
    ifelse(kb_test$p.value < 0.05, "Heteroskedastik ✓", "Homoskedastik"),
    ifelse(gq_test$p.value < 0.05, "Heteroskedastik ✓", "Homoskedastik")
  )
) |>
  kable(
    caption   = "Tabel 8a. Uji Homoskedastisitas Residual OLS",
    col.names = c("Uji", "Statistik", "df", "p-value", "Kesimpulan")
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) |>
  footnote(general = "H₀: homoskedastik. Tolak H₀ (p<0.05) → heteroskedastisitas → GWR lebih tepat dari OLS.")

# Pastikan eff_s3 tersedia di shp_merged
shp_merged <- shp_merged |> mutate(eff_s3 = Theta_S3)

gwr_sf_ll <- shp_merged |>
  filter(!is.na(eff_s3), !is.na(reg_dig), !is.na(reg_akses),
         !is.na(reg_elevasi), !is.na(x_hujan)) |>
  st_make_valid() |>
  st_transform(crs = 4326)

centroid_gwr      <- st_centroid(gwr_sf_ll)
coords_gwr        <- st_coordinates(centroid_gwr)
gwr_sf_ll$lon     <- coords_gwr[, 1]
gwr_sf_ll$lat     <- coords_gwr[, 2]

gwr_df            <- st_drop_geometry(gwr_sf_ll)
coordinates(gwr_df) <- ~ lon + lat
proj4string(gwr_df)  <- CRS("+proj=longlat +datum=WGS84")

cat("Observasi GWR:", nrow(gwr_df), "\n")

bw_gauss <- gwr.sel(formula_ols, data = gwr_df, adapt = FALSE,
                    gweight = gwr.Gauss, verbose = FALSE)
bw_bisq  <- gwr.sel(formula_ols, data = gwr_df, adapt = FALSE,
                    gweight = gwr.bisquare, verbose = FALSE)
bw_adapt <- gwr.sel(formula_ols, data = gwr_df, adapt = TRUE,
                    gweight = gwr.Gauss, verbose = FALSE)

cat("Bandwidth Gaussian fixed   :", round(bw_gauss, 6), "\n")

cat("Bandwidth Bisquare fixed   :", round(bw_bisq,  6), "\n")

cat("Bandwidth Gaussian adaptive:", round(bw_adapt, 6), "\n")

gwr_gauss <- gwr(formula_ols, data = gwr_df,
                 bandwidth = bw_gauss, gweight = gwr.Gauss,
                 hatmatrix = TRUE, se.fit = TRUE)
gwr_bisq  <- gwr(formula_ols, data = gwr_df,
                 bandwidth = bw_bisq, gweight = gwr.bisquare,
                 hatmatrix = TRUE, se.fit = TRUE)
gwr_adapt <- gwr(formula_ols, data = gwr_df,
                 adapt = bw_adapt, gweight = gwr.Gauss,
                 hatmatrix = TRUE, se.fit = TRUE)

ambil_metrik <- function(model, nama) {
  r    <- capture.output(print(model))
  get  <- function(k) {
    b <- grep(k, r, value = TRUE)
    if (!length(b)) return(NA_real_)
    as.numeric(sub(".*:\\s*", "", b[1]))
  }
  n  <- length(model$SDF$pred)
  k  <- length(all.vars(formula_ols)) - 1
  r2 <- get("Quasi-global R2")
  data.frame(
    Model  = nama,
    AIC    = round(get("^AIC"),           4),
    AICc   = round(get("^AICc"),          4),
    RSS    = round(get("Residual sum"),    4),
    R2     = round(r2,                    4),
    Adj_R2 = round(1 - (1 - r2) * (n - 1) / (n - k - 1), 4)
  )
}

tbl_bw <- rbind(
  ambil_metrik(gwr_gauss, "Gaussian Fixed"),
  ambil_metrik(gwr_bisq,  "Bisquare Fixed"),
  ambil_metrik(gwr_adapt, "Gaussian Adaptive")
)

tbl_bw |>
  kable(caption = "Tabel 9. Perbandingan Model GWR") |>
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) |>
  footnote(general = "Model terbaik: AICc terkecil dan Adj. R² tertinggi.")

idx_best  <- which.min(tbl_bw$AICc)
nama_best <- tbl_bw$Model[idx_best]
gwr_best  <- list(gwr_gauss, gwr_bisq, gwr_adapt)[[idx_best]]

cat("\nModel GWR terpilih:", nama_best,
    "| AICc =", tbl_bw$AICc[idx_best],
    "| R² =",   tbl_bw$R2[idx_best], "\n")

vars_gwr <- c("reg_dig", "reg_akses", "reg_elevasi", "x_hujan")
label_var <- c(
  reg_dig     = "Digitalisasi",
  reg_akses   = "Aksesibilitas",
  reg_elevasi = "Elevasi",
  x_hujan     = "Curah Hujan"
)

gwr_hasil <- as.data.frame(gwr_best$SDF)

for (v in vars_gwr) {
  se_col  <- paste0(v, "_se")
  gwr_hasil[[paste0("t_",   v)]]     <- gwr_hasil[[v]] / gwr_hasil[[se_col]]
  gwr_hasil[[paste0("sig_", v)]]     <- ifelse(abs(gwr_hasil[[paste0("t_", v)]]) > 1.96, "Ya", "Tidak")
  gwr_hasil[[paste0("b_",   v, "_sig")]] <- ifelse(
    abs(gwr_hasil[[paste0("t_", v)]]) > 1.96, gwr_hasil[[v]], NA_real_)
}

gwr_hasil$kodekab <- gwr_sf_ll$kodekab

gwr_sf_poly <- gwr_sf_ll |>
  left_join(
    gwr_hasil |> select(kodekab, starts_with("t_"), starts_with("sig_"),
                         starts_with("b_"), pred, localR2),
    by = "kodekab"
  )

gwr_sp <- as(gwr_sf_poly, "Spatial")
cat("Join selesai — kolom sig:", grep("^sig_", names(gwr_sp@data), value = TRUE), "\n")

sig_counts <- sapply(vars_gwr, function(v)
  sum(gwr_sp@data[[paste0("sig_", v)]] == "Ya", na.rm = TRUE))

data.frame(
  Variabel         = label_var[vars_gwr],
  Signifikan       = sig_counts,
  Tidak_Signifikan = nrow(gwr_sp@data) - sig_counts,
  Persen           = paste0(round(sig_counts / nrow(gwr_sp@data) * 100, 1), "%")
) |>
  kable(
    row.names = FALSE,
    col.names = c("Variabel", "Kab/Kota Signifikan", "Tidak Signifikan", "% Sig."),
    caption   = "Tabel 10. Jumlah Kab/Kota Berdasarkan Signifikansi Lokal GWR"
  ) |>
  kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE)

gwr_sf_poly |>
  st_drop_geometry() |>
  group_by(`Kombinasi Variabel Signifikan` = var_sig) |>
  summarise(
    `Jumlah Kab/Kota` = n(),
    `Daftar Kabupaten` = paste(sort(label_peta), collapse = ", "),
    .groups = "drop"
  ) |>
  arrange(desc(`Jumlah Kab/Kota`)) |>
  kable(caption = "Tabel 11. Rekap Kab/Kota per Kombinasi Variabel Signifikan GWR") |>
  kable_styling(bootstrap_options = c("striped", "condensed")) |>
  column_spec(3, width = "45%")