packages <- c(
  "sf", "spdep", "MASS", "car", "lmtest",
  "dplyr", "ggplot2", "readxl", "knitr",
  "kableExtra", "mgwnbr"
)

for (p in packages) {
  if (!requireNamespace(p, quietly = TRUE)) {
    install.packages(p)
  }
  library(p, character.only = TRUE)
}

sf_use_s2(FALSE)

library(readxl)

data_path <- "D:/UNY/SEMESTER 6/Statistik Spasial/Projek/Projek akhir/data1.xlsx"
shp_path  <- "D:/UNY/SEMESTER 6/Statistik Spasial/Projek/Projek akhir/jateng.shp"

# Kolom nama wilayah pada data dan shapefile
kolom_wilayah_data <- "Provinsi"
kolom_wilayah_shp  <- "WADMKK"

bersih_nama <- function(x) {
  x <- as.character(x)
  x <- iconv(x, to = "ASCII//TRANSLIT")
  x <- tolower(trimws(x))
  x <- gsub("\\s+", " ", x)
  return(x)
}

ubah_numeric <- function(x) {
  x <- as.character(x)
  x <- gsub("\\.", "", x)       # menghapus titik ribuan jika ada
  x <- gsub(",", ".", x)        # mengubah koma desimal menjadi titik
  x <- trimws(x)
  as.numeric(x)
}

ambil_aic_mgwnbr <- function(x) {
  if (is.null(x)) return(NA_real_)
  
  # kalau output berupa numeric vector
  if (is.numeric(x) && !is.null(names(x))) {
    idx <- grep("aic", names(x), ignore.case = TRUE)
    if (length(idx) > 0) return(as.numeric(x[idx[1]]))
  }
  
  # kalau output berupa data frame / matrix
  x_df <- tryCatch(as.data.frame(x), error = function(e) NULL)
  if (is.null(x_df)) return(NA_real_)
  
  kolom_aic <- grep("aic", names(x_df), ignore.case = TRUE, value = TRUE)
  if (length(kolom_aic) > 0) {
    return(as.numeric(x_df[1, kolom_aic[1]]))
  }
  
  baris_aic <- grep("aic", rownames(x_df), ignore.case = TRUE)
  if (length(baris_aic) > 0) {
    return(as.numeric(x_df[baris_aic[1], 1]))
  }
  
  return(NA_real_)
}

shp <- st_read(shp_path, quiet = TRUE)

if (grepl("\\.csv$", data_path, ignore.case = TRUE)) {
  data_raw <- read.csv2(data_path, stringsAsFactors = FALSE, check.names = FALSE)
} else if (grepl("\\.xlsx?$", data_path, ignore.case = TRUE)) {
  data_raw <- readxl::read_excel(data_path)
} else {
  stop("Format data harus .csv, .xls, atau .xlsx")
}

names(data_raw) <- trimws(names(data_raw))
names(shp)      <- trimws(names(shp))

cat("Jumlah wilayah shapefile :", nrow(shp), "\n")

cat("Jumlah baris data        :", nrow(data_raw), "\n")

cat("Kolom data               :", paste(names(data_raw), collapse = ", "), "\n")

if (!(kolom_wilayah_data %in% names(data_raw))) {
  stop("Kolom wilayah pada data tidak ditemukan.")
}

if (!(kolom_wilayah_shp %in% names(shp))) {
  stop("Kolom wilayah pada shapefile tidak ditemukan.")
}

data_raw$join_key <- bersih_nama(data_raw[[kolom_wilayah_data]])
shp$join_key      <- bersih_nama(shp[[kolom_wilayah_shp]])

cek_tidak_cocok <- setdiff(data_raw$join_key, shp$join_key)

if (length(cek_tidak_cocok) > 0) {
  cat("\nNama wilayah di data yang tidak cocok dengan shapefile:\n")
  print(cek_tidak_cocok)
}

data_sf <- shp %>%
  left_join(data_raw, by = "join_key")

cat("\nHasil join:", nrow(data_sf), "wilayah\n")

variabel_model <- c("Y", "X1", "X2", "X3", "X4")

for (v in variabel_model) {
  if (!(v %in% names(data_sf))) {
    stop(paste("Variabel", v, "tidak ditemukan pada data."))
  }
}

data_sf <- data_sf %>%
  mutate(
    Y  = ubah_numeric(Y),
    X1 = ubah_numeric(X1),
    X2 = ubah_numeric(X2),
    X3 = ubah_numeric(X3),
    X4 = ubah_numeric(X4),
  )

# Y harus berupa data cacah
data_sf$Y_count <- as.integer(round(data_sf$Y))

# Cek missing value
cat("\nMissing value sebelum dibersihkan:\n")

print(colSums(is.na(st_drop_geometry(data_sf)[, c("Y_count", "X1", "X2", "X3", "X4")])))

# Pakai hanya data lengkap
cc <- complete.cases(st_drop_geometry(data_sf)[, c("Y_count", "X1", "X2", "X3", "X4")])
data_sf <- data_sf[cc, ]

if (any(data_sf$Y_count < 0, na.rm = TRUE)) {
  stop("Y_count mengandung nilai negatif. Data count tidak boleh negatif.")
}

cat("\nJumlah data setelah cleaning:", nrow(data_sf), "\n")

if (is.na(st_crs(data_sf))) {
  stop("CRS shapefile tidak terbaca. Tetapkan CRS dulu sebelum analisis spasial.")
}

data_sf <- st_make_valid(data_sf)
data_sf <- st_transform(data_sf, 4326)

centroid_sf <- suppressWarnings(st_point_on_surface(data_sf))
coords <- st_coordinates(centroid_sf)

data_sf$long <- coords[, 1]
data_sf$lat  <- coords[, 2]

# Data tanpa geometri untuk model
df_model <- data_sf %>%
  st_drop_geometry() %>%
  mutate(
    wilayah = .data[[kolom_wilayah_shp]]
  )

x_vars <- c("X1", "X2", "X3", "X4")
form_model <- as.formula(
  paste("Y_count ~", paste(x_vars, collapse = " + "))
)

cat("\nFormula model yang digunakan:\n")

print(form_model)

df_desc <- df_model[, c("Y_count", x_vars)]

desc_tbl <- data.frame(
  Variabel = names(df_desc),
  N        = sapply(df_desc, function(x) sum(!is.na(x))),
  Min      = sapply(df_desc, min, na.rm = TRUE),
  Mean     = sapply(df_desc, mean, na.rm = TRUE),
  Median   = sapply(df_desc, median, na.rm = TRUE),
  Max      = sapply(df_desc, max, na.rm = TRUE),
  SD       = sapply(df_desc, sd, na.rm = TRUE),
  row.names = NULL
)

kable(desc_tbl, caption = "Statistik Deskriptif Variabel Penelitian", digits = 4) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

model_vif <- lm(form_model, data = df_model)
vif_val <- car::vif(model_vif)

vif_tbl <- data.frame(
  Variabel   = names(vif_val),
  VIF        = round(vif_val, 4),
  Keputusan  = ifelse(vif_val > 10, "Ada multikolinearitas", "Tidak ada multikolinearitas"),
  row.names  = NULL
)

kable(vif_tbl, caption = "Nilai VIF Prediktor") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

model_poisson <- glm(
  form_model,
  family = poisson(link = "log"),
  data = df_model
)

cat("\n==============================\n")

cat("MODEL POISSON\n")

cat("==============================\n")

print(summary(model_poisson))

dev_pois  <- deviance(model_poisson)
df_pois   <- df.residual(model_poisson)

rasio_dev <- dev_pois / df_pois
rasio_pearson <- sum(residuals(model_poisson, type = "pearson")^2) / df_pois

cat("\n==============================\n")

cat("UJI OVERDISPERSI\n")

cat("==============================\n")

cat("Deviance Poisson        :", round(dev_pois, 4), "\n")

cat("df residual             :", df_pois, "\n")

cat("Rasio Deviance/df       :", round(rasio_dev, 4), "\n")

cat("Rasio Pearson ChiSq/df  :", round(rasio_pearson, 4), "\n")

if (rasio_dev > 1 || rasio_pearson > 1) {
  cat("Keputusan: Terjadi overdispersi.\n")
  cat("Artinya model Negative Binomial lebih sesuai dibanding Poisson.\n")
} else {
  cat("Keputusan: Tidak terindikasi overdispersi.\n")
}

model_nb <- MASS::glm.nb(
  form_model,
  data = df_model
)

cat("\n==============================\n")

cat("MODEL NEGATIVE BINOMIAL GLOBAL\n")

cat("==============================\n")

print(summary(model_nb))

theta_nb <- model_nb$theta
alpha_nb <- 1 / theta_nb

cat("\nTheta NB dari glm.nb :", round(theta_nb, 6), "\n")

cat("Alpha NB = 1/theta   :", round(alpha_nb, 6), "\n")

cat("AIC Poisson          :", round(AIC(model_poisson), 4), "\n")

cat("AIC NB Global        :", round(AIC(model_nb), 4), "\n")

model_nb_null <- MASS::glm.nb(Y_count ~ 1, data = df_model)

D_nb <- as.numeric(2 * (logLik(model_nb) - logLik(model_nb_null)))
df_nb <- length(coef(model_nb)) - 1
chi_nb <- qchisq(0.95, df_nb)
p_nb_simultan <- pchisq(D_nb, df = df_nb, lower.tail = FALSE)

cat("\n==============================\n")

cat("UJI SIMULTAN MODEL NB\n")

cat("==============================\n")

cat("D hitung        :", round(D_nb, 4), "\n")

cat("df              :", df_nb, "\n")

cat("Chi-square 0.95 :", round(chi_nb, 4), "\n")

cat("p-value         :", round(p_nb_simultan, 6), "\n")

if (p_nb_simultan < 0.05) {
  cat("Keputusan: Tolak H0. Minimal ada satu prediktor berpengaruh.\n")
} else {
  cat("Keputusan: Gagal tolak H0.\n")
}

coef_nb <- summary(model_nb)$coefficients

parsial_nb <- data.frame(
  Variabel  = rownames(coef_nb),
  Koefisien = round(coef_nb[, 1], 6),
  SE        = round(coef_nb[, 2], 6),
  Z_hitung  = round(coef_nb[, 3], 6),
  P_value   = round(coef_nb[, 4], 6),
  Keputusan = ifelse(coef_nb[, 4] < 0.05, "Signifikan", "Tidak signifikan"),
  row.names = NULL
)

kable(parsial_nb, caption = "Uji Parsial Model Negative Binomial Global") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

nb_poly <- spdep::poly2nb(data_sf, queen = TRUE)

lw_poly <- spdep::nb2listw(nb_poly, style = "W", zero.policy = TRUE)

res_nb <- residuals(model_nb, type = "pearson")

moran_nb <- spdep::moran.test(
  res_nb,
  lw_poly,
  zero.policy = TRUE
)

cat("\n==============================\n")

cat("MORAN'S I RESIDUAL NB\n")

cat("==============================\n")

print(moran_nb)

p_moran <- moran_nb$p.value
I_moran <- moran_nb$estimate[["Moran I statistic"]]

cat("\nMoran's I :", round(I_moran, 6), "\n")

cat("p-value   :", round(p_moran, 6), "\n")

if (p_moran < 0.05) {
  cat("Keputusan: Ada autokorelasi spasial pada residual.\n")
  cat("Artinya analisis GWNBR layak dilanjutkan.\n")
} else {
  cat("Keputusan: Tidak terdapat autokorelasi spasial yang signifikan.\n")
  cat("Catatan: GWNBR masih bisa dicoba, tetapi alasan spasialnya lebih lemah.\n")
}

bp_test <- lmtest::bptest(model_vif)

cat("\n==============================\n")

cat("UJI HETEROGENITAS BREUSCH-PAGAN\n")

cat("==============================\n")

print(bp_test)

if (bp_test$p.value < 0.05) {
  cat("Keputusan: Tolak H0. Ada heterogenitas.\n")
  cat("Artinya parameter lokal GWNBR layak dipertimbangkan.\n")
} else {
  cat("Keputusan: Gagal tolak H0. Heterogenitas belum signifikan.\n")
}

options(warn = -1)
model_gwnbr <- mgwnbr::mgwnbr(
  data           = df_model,
  formula        = form_model,
  long           = "long",
  lat            = "lat",
  band_method    = "adaptive_bsq",
  band_criterion = "aic",
  distribution   = "negbin",
  globalmin      = TRUE,
  multiscale     = FALSE,
  distancekm     = TRUE,
  id             = "wilayah",
  max_int        = 50
)

cat("\n==============================\n")

cat("MODEL GWNBR\n")

cat("==============================\n")

print(model_gwnbr)

cat("\nRingkasan model GWNBR:\n")

print(summary(model_gwnbr))

cat("\nNama komponen output mgwnbr:\n")

print(names(model_gwnbr))

cat("\n==============================\n")

cat("GOODNESS OF FIT GWNBR\n")

cat("==============================\n")

if (!is.null(model_gwnbr$general_bandwidth)) {
  cat("\nBandwidth umum:\n")
  print(model_gwnbr$general_bandwidth)
}

if (!is.null(model_gwnbr$band)) {
  cat("\nBandwidth:\n")
  print(model_gwnbr$band)
}

if (!is.null(model_gwnbr$measures)) {
  cat("\nMeasures GWNBR:\n")
  print(model_gwnbr$measures)
}

if (!is.null(model_gwnbr$global_measures)) {
  cat("\nGlobal measures dari mgwnbr:\n")
  print(model_gwnbr$global_measures)
}

if (!is.null(model_gwnbr$mgwr_param_estimates)) {
  
  coef_lokal <- as.data.frame(model_gwnbr$mgwr_param_estimates)
  coef_lokal$Wilayah <- df_model$wilayah
  
  coef_lokal <- coef_lokal %>%
    dplyr::select(Wilayah, everything())
  
  kable(coef_lokal,
        caption = "Koefisien Lokal GWNBR per Kabupaten/Kota",
        digits = 6) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE) %>%
    scroll_box(height = "400px")
  
} else {
  cat("Komponen mgwr_param_estimates tidak ditemukan pada output mgwnbr.\n")
}

if (!is.null(model_gwnbr$mgwr_param_estimates) &&
    !is.null(model_gwnbr$mgwr_se)) {
  
  beta_mat <- as.matrix(model_gwnbr$mgwr_param_estimates)
  se_mat   <- as.matrix(model_gwnbr$mgwr_se)
  
  z_mat <- beta_mat / se_mat
  p_mat <- 2 * (1 - pnorm(abs(z_mat)))
  
  sig_mat <- p_mat < 0.05
  
  sig_tbl <- as.data.frame(sig_mat)
  sig_tbl$Wilayah <- df_model$wilayah
  
  sig_tbl <- sig_tbl %>%
    dplyr::select(Wilayah, everything())
  
  kable(sig_tbl,
        caption = "Signifikansi Lokal Parameter GWNBR, TRUE = Signifikan") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE) %>%
    scroll_box(height = "400px")
  
  # Pengelompokan wilayah berdasarkan variabel signifikan
  kelompok_sig <- apply(sig_mat, 1, function(x) {
    vars <- colnames(sig_mat)[x]
    vars <- vars[!vars %in% c("(Intercept)", "Intercept")]
    if (length(vars) == 0) {
      return("Tidak ada variabel signifikan")
    } else {
      return(paste(vars, collapse = ", "))
    }
  })
  
  data_sf$Kelompok_Signifikan <- kelompok_sig
  
  kelompok_tbl <- data.frame(
    Wilayah = df_model$wilayah,
    Kelompok_Signifikan = kelompok_sig
  )
  
  kable(kelompok_tbl,
        caption = "Kelompok Wilayah Berdasarkan Variabel Signifikan pada GWNBR") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
  
} else {
  cat("Komponen mgwr_param_estimates atau mgwr_se tidak ditemukan.\n")
}

AIC_pois <- AIC(model_poisson)
AIC_nb   <- AIC(model_nb)
AIC_gwnbr <- ambil_aic_mgwnbr(model_gwnbr$measures)

perbandingan_aic <- data.frame(
  Model = c("Regresi Poisson", "Regresi Negative Binomial Global", "GWNBR"),
  AIC   = c(round(AIC_pois, 4), round(AIC_nb, 4), round(AIC_gwnbr, 4))
)

kable(perbandingan_aic,
      caption = "Perbandingan AIC Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

cat("\n==============================\n")

cat("PERBANDINGAN AIC\n")

cat("==============================\n")

print(perbandingan_aic)

if (!is.na(AIC_gwnbr)) {
  if (AIC_gwnbr < AIC_nb) {
    cat("\nKesimpulan AIC: AIC GWNBR lebih kecil daripada AIC NB global.\n")
    cat("Artinya GWNBR memberikan model yang lebih baik berdasarkan AIC.\n")
  } else {
    cat("\nKesimpulan AIC: AIC GWNBR belum lebih kecil daripada AIC NB global.\n")
    cat("Coba cek kembali bandwidth, variabel, atau struktur spasial.\n")
  }
} else {
  cat("\nAIC GWNBR tidak berhasil diekstrak otomatis.\n")
  cat("Lihat langsung pada output model_gwnbr$measures.\n")
}

cat("1. Model yang digunakan:\n")

cat("   Y_count ~ X1 + X2 + X3 + X4\n\n")

cat("2. Overdispersi:\n")

cat("   Rasio Deviance/df      =", round(rasio_dev, 4), "\n")

cat("   Rasio Pearson ChiSq/df =", round(rasio_pearson, 4), "\n")

if (rasio_dev > 1 || rasio_pearson > 1) {
  cat("   Kesimpulan: terjadi overdispersi, sehingga NB lebih sesuai daripada Poisson.\n\n")
} else {
  cat("   Kesimpulan: overdispersi belum terlihat kuat.\n\n")
}

cat("3. Model Negative Binomial global:\n")

cat("   AIC NB =", round(AIC_nb, 4), "\n")

cat("   Theta =", round(theta_nb, 6), "\n")

cat("   Alpha =", round(alpha_nb, 6), "\n\n")

cat("4. Moran's I residual NB:\n")

cat("   Moran's I =", round(I_moran, 6), "\n")

cat("   p-value   =", round(p_moran, 6), "\n")

if (p_moran < 0.05) {
  cat("   Kesimpulan: terdapat autokorelasi spasial.\n\n")
} else {
  cat("   Kesimpulan: autokorelasi spasial belum signifikan.\n\n")
}

cat("5. Heterogenitas Breusch-Pagan:\n")

cat("   BP statistic =", round(bp_test$statistic, 6), "\n")

cat("   p-value      =", round(bp_test$p.value, 6), "\n")

if (bp_test$p.value < 0.05) {
  cat("   Kesimpulan: terdapat heterogenitas.\n\n")
} else {
  cat("   Kesimpulan: heterogenitas belum signifikan.\n\n")
}

cat("6. Perbandingan AIC:\n")

cat("   AIC Poisson =", round(AIC_pois, 4), "\n")

cat("   AIC NB      =", round(AIC_nb, 4), "\n")

cat("   AIC GWNBR   =", round(AIC_gwnbr, 4), "\n")

if (!is.na(AIC_gwnbr) && AIC_gwnbr < AIC_nb) {
  cat("   Kesimpulan: GWNBR lebih baik daripada NB global berdasarkan AIC.\n")
} else if (!is.na(AIC_gwnbr) && AIC_gwnbr >= AIC_nb) {
  cat("   Kesimpulan: GWNBR belum lebih baik daripada NB global berdasarkan AIC.\n")
} else {
  cat("   Kesimpulan: AIC GWNBR perlu dilihat manual pada model_gwnbr$measures.\n")
}

# ============================================================
# VISUALISASI PETA TEMATIK
# ============================================================

library(ggplot2)
library(sf)
library(dplyr)
library(ggrepel)
library(RColorBrewer)
library(scales)

# ── Label wilayah (centroid) ─────────────────────────────────
# Gunakan st_point_on_surface agar titik selalu jatuh di dalam polygon
centroids_label <- suppressWarnings(st_point_on_surface(st_geometry(data_sf)))
coords_label    <- st_coordinates(centroids_label)

df_label <- data.frame(
  wilayah = data_sf[[kolom_wilayah_shp]],
  lx      = coords_label[, 1],
  ly      = coords_label[, 2],
  stringsAsFactors = FALSE
)

# ── Tema dasar peta ──────────────────────────────────────────
tema_peta <- theme_void(base_family = "sans") +
  theme(
    plot.title        = element_text(face = "bold", size = 14, hjust = 0.5,
                                     margin = margin(b = 8)),
    plot.subtitle     = element_text(size = 10, hjust = 0.5, color = "grey40",
                                     margin = margin(b = 6)),
    legend.position   = "right",
    legend.title      = element_text(size = 9, face = "bold"),
    legend.text       = element_text(size = 8),
    plot.margin       = margin(10, 10, 10, 10),
    panel.background  = element_rect(fill = "#f0f4f8", color = NA),
    plot.background   = element_rect(fill = "white", color = NA)
  )
# ────────────────────────────────────────────────────────────
# FUNGSI UTAMA: peta khoropleth satu variabel
# ────────────────────────────────────────────────────────────
buat_peta_khoropleth <- function(sf_data, var_col, judul, subjudul = "",
                                  label_df, palet = "YlOrRd",
                                  arah_palet = 1) {
  
  sf_data$nilai_plot <- sf_data[[var_col]]
  
  ggplot(sf_data) +
    geom_sf(aes(fill = nilai_plot), color = "white", linewidth = 0.4) +
    scale_fill_distiller(
      palette   = palet,
      direction = arah_palet,
      name      = var_col,
      labels    = label_comma(big.mark = ".", decimal.mark = ","),
      guide     = guide_colorbar(barwidth = 0.8, barheight = 8,
                                  title.position = "top")
    ) +
    geom_text_repel(
      data          = label_df,
      aes(x = lx, y = ly, label = wilayah),
      size          = 2.2,
      fontface      = "bold",
      color         = "black",
      bg.color      = "white",
      bg.r          = 0.15,
      segment.color = "grey50",
      segment.size  = 0.3,
      max.overlaps  = 50,
      min.segment.length = 0.1,
      force         = 1.5,
      seed          = 42
    ) +
    labs(title = judul, subtitle = subjudul) +
    tema_peta
}

# ============================================================
# PETA 1 – SEBARAN VARIABEL Y (Jumlah Kasus / Count)
# ============================================================
p_Y <- buat_peta_khoropleth(
  sf_data   = data_sf,
  var_col   = "Y_count",
  judul     = "Peta Sebaran Variabel Y",
  subjudul  = "Jumlah Kasus per Kabupaten/Kota di Jawa Tengah",
  label_df  = df_label,
  palet     = "YlOrRd",
  arah_palet = 1
)

print(p_Y)

ggsave("peta_Y.png", p_Y, width = 12, height = 8, dpi = 200)


# ============================================================
# PETA 2–5 – SEBARAN VARIABEL X1, X2, X3, X4
# ============================================================

info_x <- list(
  X1 = list(judul    = "Peta Sebaran Variabel X1",
             subjudul = "Kabupaten/Kota di Jawa Tengah",
             palet    = "Blues",   arah = 1),
  X2 = list(judul    = "Peta Sebaran Variabel X2",
             subjudul = "Kabupaten/Kota di Jawa Tengah",
             palet    = "Greens",  arah = 1),
  X3 = list(judul    = "Peta Sebaran Variabel X3",
             subjudul = "Kabupaten/Kota di Jawa Tengah",
             palet    = "Purples", arah = 1),
  X4 = list(judul    = "Peta Sebaran Variabel X4",
             subjudul = "Kabupaten/Kota di Jawa Tengah",
             palet    = "Oranges", arah = 1)
)

peta_list <- list()

for (xv in names(info_x)) {
  info <- info_x[[xv]]
  p <- buat_peta_khoropleth(
    sf_data    = data_sf,
    var_col    = xv,
    judul      = info$judul,
    subjudul   = info$subjudul,
    label_df   = df_label,
    palet      = info$palet,
    arah_palet = info$arah
  )
  peta_list[[xv]] <- p
  print(p)
  ggsave(paste0("peta_", xv, ".png"), p, width = 12, height = 8, dpi = 200)
}

# ============================================================
# PETA 6 – PENGELOMPOKAN WILAYAH BERDASARKAN SIGNIFIKANSI GWNBR
# ============================================================

# Pastikan kolom Kelompok_Signifikan sudah ada (dibuat di chunk 21)
# Bila belum ada, jalankan ulang chunk 21 terlebih dahulu.

if (!("Kelompok_Signifikan" %in% names(data_sf))) {
  stop("Kolom Kelompok_Signifikan belum ada. Jalankan chunk 21 terlebih dahulu.")
}

# Buat palet warna kategorikal yang kontras
kelompok_unik <- sort(unique(data_sf$Kelompok_Signifikan))
n_kelompok    <- length(kelompok_unik)

# Palet kontras: gabungkan beberapa set warna yang tajam
palet_kontras <- c(
  "#E63946",  # merah terang
  "#2A9D8F",  # tosca
  "#F4A261",  # oranye hangat
  "#457B9D",  # biru tua
  "#8338EC",  # ungu cerah
  "#06D6A0",  # hijau mint
  "#FFB703",  # kuning amber
  "#FB5607",  # oranye tua
  "#3A86FF",  # biru cerah
  "#FFBE0B",  # kuning
  "#D62246",  # merah tua
  "#4CC9F0",  # biru langit
  "#B5179E",  # magenta
  "#4D908E",  # hijau kebiruan
  "#F72585",  # pink terang
  "#90E0EF",  # biru pastel kontras
  "#43AA8B",  # hijau sage
  "#577590"   # abu kebiruan
)

# Kalau jumlah kelompok > panjang vektor, pakai rainbow
if (n_kelompok > length(palet_kontras)) {
  palet_pakai <- hcl.colors(n_kelompok, palette = "Dynamic")
} else {
  palet_pakai <- palet_kontras[seq_len(n_kelompok)]
}
names(palet_pakai) <- kelompok_unik

p_sig <- ggplot(data_sf) +
  geom_sf(
    aes(fill = Kelompok_Signifikan),
    color    = "white",
    linewidth = 0.45
  ) +
  scale_fill_manual(
    values = palet_pakai,
    name   = "Variabel Signifikan",
    guide  = guide_legend(
      ncol         = 1,
      byrow        = TRUE,
      keywidth     = unit(0.6, "cm"),
      keyheight    = unit(0.5, "cm"),
      override.aes = list(color = "white", linewidth = 0.3)
    )
  ) +
  geom_text_repel(
    data          = df_label,
    aes(x = lx, y = ly, label = wilayah),
    size          = 2.2,
    fontface      = "bold",
    color         = "black",
    bg.color      = "white",
    bg.r          = 0.15,
    segment.color = "grey50",
    segment.size  = 0.3,
    max.overlaps  = 50,
    min.segment.length = 0.1,
    force         = 1.5,
    seed          = 42
  ) +
  labs(
    title    = "Peta Pengelompokan Wilayah Berdasarkan Variabel Signifikan GWNBR",
    subtitle = "Kabupaten/Kota di Jawa Tengah"
  ) +
  tema_peta +
  theme(
    legend.text  = element_text(size = 7.5),
    legend.title = element_text(size = 9, face = "bold")
  )

print(p_sig)

ggsave("peta_signifikansi_gwnbr.png", p_sig, width = 14, height = 9, dpi = 200)


# ============================================================
# BONUS: Panel 6 peta sekaligus (opsional, butuh patchwork)
# ============================================================

if (requireNamespace("patchwork", quietly = TRUE)) {
  library(patchwork)
  
  panel_semua <- (p_Y | peta_list$X1 | peta_list$X2) /
                 (peta_list$X3 | peta_list$X4 | p_sig) +
    plot_annotation(
      title   = "Analisis Spasial GWNBR – Jawa Tengah",
      caption = "Sumber: Data Penelitian | Disusun: Nisrina Aisyah & Nayli Kurnia Ilahi",
      theme   = theme(
        plot.title   = element_text(face = "bold", size = 16, hjust = 0.5),
        plot.caption = element_text(size = 8, color = "grey50")
      )
    )
  
  ggsave("panel_semua_peta.png", panel_semua,
         width = 26, height = 16, dpi = 180)
  cat("Panel lengkap disimpan: panel_semua_peta.png\n")
  
} else {
  cat("Install patchwork untuk panel gabungan: install.packages('patchwork')\n")
}

cat("\nSemua peta berhasil dibuat dan disimpan.\n")