1 Hasil analisis

# ANALISIS MGWR FIXED BANDWIDTH - TB RATE 514 KABUPATEN/KOTA
# ---------------------- BAGIAN KODE 1 ----------------------
EXCEL_FILE <- "D:/STATISTIKA/SEM 5/Analisis Data Spasial/Materi Pekan 1-8/Tugas Project  GWR/Data Spasial_TBC/Data TBC Kabupaten Kota di Indonesia/Data TBC.xlsx"
SHP_FILE <- "D:/STATISTIKA/SEM 5/Analisis Data Spasial/Materi Pekan 1-8/Tugas Project  GWR/BATAS/Batas_Kabupaten_Kota.gpkg"

SHP_NAME_FIELD <- "WADMKK"

Y_TRANSFORM <- "raw"
STANDARDIZE <- TRUE

RUN_KERNEL_COMPARISON <- FALSE  
KERNEL_CANDIDATES <- c("exponential")
KERNEL_DEFAULT    <- "exponential"

ADAPTIVE <- FALSE 
LONG_LAT <- TRUE
K_MORAN  <- 8                

MGWR_CRITERION <- "dCVR"
MGWR_APPROACH  <- "AIC"
MGWR_MAX_ITER  <- 2000
MGWR_THRESHOLD <- 1e-5
MGWR_REOPT     <- 5
MGWR_NLOWER    <- 10

RUN_MONTE_CARLO <- TRUE    
MC_SIM <- 99
T_CRIT <- 1.96              
set.seed(123)

# ---------------------- BAGIAN KODE 2 ----------------------
pkgs <- c("readxl", "dplyr", "sf", "sp", "spdep", "car", "lmtest",
          "GWmodel", "ggplot2", "corrplot", "DT", "knitr", "rmarkdown")
new_pkgs <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if(length(new_pkgs) > 0) install.packages(new_pkgs, dependencies = TRUE)

library(readxl); library(dplyr); library(sf); library(sp); library(spdep)
library(car); library(lmtest); library(GWmodel); library(ggplot2); library(corrplot)

sf_use_s2(FALSE)

OUT <- "hasil_MGWR_TBC_fixed"
MAP <- file.path(OUT, "peta")
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)
dir.create(MAP, showWarnings = FALSE, recursive = TRUE)

# ---------------------- BAGIAN KODE 3 ----------------------
clean_text <- function(x){
  x <- toupper(trimws(as.character(x)))
  x <- gsub("[’']", "", x)
  x <- gsub("\\.", " ", x)
  x <- gsub("&", " DAN ", x, fixed = TRUE)
  x <- gsub("[^A-Z0-9 ]", " ", x)
  x <- gsub("\\s+", " ", x)
  trimws(x)
}
remove_prefix <- function(x){
  x <- clean_text(x)
  x <- gsub("^KABUPATEN\\s+", "", x)
  x <- gsub("^KAB\\s+", "", x)
  x <- gsub("^KOTA ADMINISTRATIF\\s+", "", x)
  x <- gsub("^KOTA\\s+", "", x)
  x <- gsub("^REGENCY\\s+", "", x)
  x <- gsub("^CITY\\s+", "", x)
  trimws(x)
}
get_type <- function(x){
  z <- clean_text(x)
  ifelse(grepl("KABUPATEN|KAB |REGENCY", z), "KABUPATEN",
         ifelse(grepl("KOTA|CITY", z), "KOTA", NA_character_))
}
make_key <- function(type, name) paste0(type, "|", remove_prefix(name))

save_txt <- function(x, file){
  capture.output(print(x), file = file.path(OUT, file))
}

get_stat <- function(obj, candidates){
  if(is.null(obj)) return(NA_real_)
  for(nm in candidates){
    if(nm %in% names(obj)){
      v <- suppressWarnings(as.numeric(obj[[nm]][1]))
      if(length(v) == 1 && is.finite(v)) return(v)
    }
  }
  NA_real_
}

find_col <- function(df, candidates){
  nm <- names(df)
  for(cand in candidates){
    hit <- nm[tolower(nm) == tolower(cand)]
    if(length(hit) > 0) return(hit[1])
  }
  NA_character_
}

# AICc OLS: AIC + koreksi sampel kecil.
# k menghitung parameter regresi dan parameter varians residual.
aicc_ols <- function(m){
  nn <- stats::nobs(m)
  kk <- length(stats::coef(m)) + 1L
  if(nn <= kk + 1L) return(NA_real_)
  stats::AIC(m) + (2 * kk * (kk + 1)) / (nn - kk - 1)
}

save_map <- function(sf_obj, col, title, file, diverging = FALSE, legend = "Nilai"){
  p <- ggplot(sf_obj) +
    geom_sf(aes(fill = .data[[col]]), color = "grey70", linewidth = 0.03) +
    theme_void() + labs(title = title, fill = legend)
  p <- if(diverging){
    p + scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                             midpoint = 0, na.value = "grey85")
  }else{
    p + scale_fill_viridis_c(option = "C", na.value = "grey85")
  }
  ggsave(file.path(MAP, file), p, width = 12, height = 6, dpi = 300)
}

# ---------------------- BAGIAN KODE 4 ----------------------
if(!file.exists(EXCEL_FILE)) stop("File Excel tidak ditemukan di working directory.")
raw <- read_excel(EXCEL_FILE, sheet = 1)
cat("Jumlah baris Excel:", nrow(raw), "\n"); print(names(raw))

required_cols <- c(
  "Kabupaten/Kota",
  "TB_Rate_per100.000_Penduduk",
  "Persentase_Penduduk_Miskin",
  "Persentase_Rumah_Tangga_dengan_Sanitasi_Layak",
  "Persentase_Penduduk15_yangMerokok",
  "Persentase_Rumah_Tangga_Memiliki_Akses_Hunian_Layak"
)
if(!all(required_cols %in% names(raw))){
  stop("Nama kolom Excel berbeda. Kolom yang hilang: ",
       paste(setdiff(required_cols, names(raw)), collapse = ", "))
}

data <- raw %>%
  transmute(
    id          = row_number(),
    wilayah     = as.character(`Kabupaten/Kota`),
    Y_TB        = as.numeric(`TB_Rate_per100.000_Penduduk`),
    X1_miskin   = as.numeric(`Persentase_Penduduk_Miskin`),
    X2_sanitasi = as.numeric(`Persentase_Rumah_Tangga_dengan_Sanitasi_Layak`),
    X3_merokok  = as.numeric(`Persentase_Penduduk15_yangMerokok`),
    X4_hunian   = as.numeric(`Persentase_Rumah_Tangga_Memiliki_Akses_Hunian_Layak`)
  ) %>%
  as.data.frame()

data$jenis      <- get_type(data$wilayah)
data$nama_dasar <- remove_prefix(data$wilayah)
data$key        <- make_key(data$jenis, data$nama_dasar)

xvars <- c("X1_miskin", "X2_sanitasi", "X3_merokok", "X4_hunian")
vars  <- c("Y_TB", xvars)

# Audit data
miss <- data.frame(variable = vars, n_missing = sapply(data[vars], function(x) sum(is.na(x))))
write.csv(miss, file.path(OUT, "01_missing.csv"), row.names = FALSE)
if(any(miss$n_missing > 0)) stop("Ada missing value. Lihat 01_missing.csv")
if(Y_TRANSFORM == "log" && any(data$Y_TB <= 0)) stop("TB rate harus > 0 untuk transformasi log.")
if(anyDuplicated(data$key) > 0) stop("Ada wilayah duplikat di Excel setelah normalisasi nama.")
if(nrow(data) != 514) warning("Jumlah baris bukan 514: ", nrow(data))
if(any(sapply(data[xvars], sd) == 0)) stop("Ada prediktor dengan variasi nol.")
write.csv(data, file.path(OUT, "02_data_clean.csv"), row.names = FALSE)

# ---------------------- BAGIAN KODE 5 ----------------------
# Memastikan file GeoPackage tersedia
if (!file.exists(SHP_FILE)) {
  stop("File GeoPackage tidak ditemukan: ", SHP_FILE)
}

# Baca layer kabupaten/kota
print(sf::st_layers(SHP_FILE))

admin <- sf::st_read(
  SHP_FILE,
  layer = "kabupaten_kota",
  quiet = TRUE
)

# Periksa kolom nama wilayah
print(names(admin))

# Gunakan kolom WADMKK sebagai nama wilayah
if (!SHP_NAME_FIELD %in% names(admin)) {
  stop("Kolom ", SHP_NAME_FIELD, " tidak ditemukan dalam GeoPackage.")
}

admin$name_shape <- as.character(admin[[SHP_NAME_FIELD]])

# Tentukan jenis wilayah dari nama pada peta
admin$type_shape <- ifelse(
  grepl("^KOTA\\s+", clean_text(admin$name_shape)),
  "KOTA",
  "KABUPATEN"
)

admin$name_only <- remove_prefix(admin$name_shape)

# Alias nama untuk menyamakan ejaan Excel dan peta
alias_nama <- c(
  "TOBA SAMOSIR" = "TOBA",
  "LABUHAN BATU" = "LABUHANBATU",
  "LABUHAN BATU SELATAN" = "LABUHANBATU SELATAN",
  "LABUHAN BATU UTARA" = "LABUHANBATU UTARA",
  "TANJUNG BALAI" = "TANJUNGBALAI",
  "PEMATANG SIANTAR" = "PEMATANGSIANTAR",
  "SAWAH LUNTO" = "SAWAHLUNTO",
  "BATANG HARI" = "BATANGHARI",
  "BANYU ASIN" = "BANYUASIN",
  "PENUKAL ABAB" = "PENUKAL ABAB LEMATANG ILIR",
  "LUBUKLINGGAU" = "LUBUK LINGGAU",
  "TULANGBAWANG" = "TULANG BAWANG",
  "KEPULAUAN SERIBU" = "ADMINISTRASI KEPULAUAN SERIBU",
  "JAKARTA SELATAN" = "ADMINISTRASI JAKARTA SELATAN",
  "JAKARTA TIMUR" = "ADMINISTRASI JAKARTA TIMUR",
  "JAKARTA PUSAT" = "ADMINISTRASI JAKARTA PUSAT",
  "JAKARTA BARAT" = "ADMINISTRASI JAKARTA BARAT",
  "JAKARTA UTARA" = "ADMINISTRASI JAKARTA UTARA",
  "GUNUNG KIDUL" = "GUNUNGKIDUL",
  "PALANGKA RAYA" = "PALANGKARAYA",
  "BANJAR BARU" = "BANJARBARU",
  "MAHAKAM HULU" = "MAHAKAM ULU",
  "BAUBAU" = "BAU BAU",
  "MALUKU TENGGARA BARAT" = "KEPULAUAN TANIMBAR",
  "FAKFAK" = "FAK FAK"
)

# Normalisasi nama dari Excel
data$name_match <- data$nama_dasar
hit_alias <- data$name_match %in% names(alias_nama)
data$name_match[hit_alias] <-
  unname(alias_nama[data$name_match[hit_alias]])

# Normalisasi nama pada peta
admin$name_match <- admin$name_only

# Buat kunci gabungan jenis dan nama
data$key_match <- paste(data$jenis, data$name_match, sep = "|")
admin$key_match <- paste(admin$type_shape, admin$name_match, sep = "|")

# Periksa jumlah wilayah
cat("Jumlah baris Excel:", nrow(data), "\n")
cat("Jumlah fitur GeoPackage:", nrow(admin), "\n")

if (nrow(admin) != 514) {
  stop("Jumlah fitur GeoPackage bukan 514. Periksa layer yang dibaca.")
}

# Audit geometri: simpan laporan dan hentikan proses bila geometri kosong/tidak valid.
geom_audit <- data.frame(
  n_fitur = nrow(admin),
  n_geometri_kosong = sum(sf::st_is_empty(admin)),
  n_geometri_tidak_valid = sum(!sf::st_is_valid(admin), na.rm = TRUE)
)
write.csv(geom_audit, file.path(OUT, "03_geometry_audit.csv"), row.names = FALSE)
if (geom_audit$n_geometri_kosong > 0 || geom_audit$n_geometri_tidak_valid > 0) {
  stop("Ada geometri kosong/tidak valid. Periksa 03_geometry_audit.csv sebelum analisis.")
}

# Periksa apakah kunci peta duplikat
if (anyDuplicated(admin$key_match) > 0) {
  print(admin[duplicated(admin$key_match),
              c("WADMKK", "WADMPR", "KDPKAB", "key_match")])
  stop("Kunci jenis dan nama wilayah pada peta duplikat.")
}

# Cocokkan urutan wilayah Excel dengan peta
idx <- match(data$key_match, admin$key_match)

cat("Wilayah tidak cocok:", sum(is.na(idx)), "\n")

if (any(is.na(idx))) {
  print(data[is.na(idx),
             c("id", "wilayah", "jenis", "name_match")])
  stop("Matching belum lengkap. Periksa nama wilayah yang tidak cocok.")
}

if (anyDuplicated(idx) > 0) {
  stop("Ada geometri yang digunakan lebih dari sekali.")
}

# Gabungkan data dan geometri
analysis_poly <- sf::st_sf(
  data,
  geometry = sf::st_geometry(admin)[idx]
)

stopifnot(
  nrow(analysis_poly) == 514,
  all(analysis_poly$id == data$id)
)

cat("\nMATCHING BERHASIL\n")
cat("Jumlah observasi:", nrow(analysis_poly), "\n")
cat("Jumlah geometri:", nrow(analysis_poly), "\n")
cat("Wilayah tidak cocok:", sum(is.na(idx)), "\n")

# ---------------------- BAGIAN KODE 6 ----------------------
desc <- do.call(rbind, lapply(vars, function(v){
  x <- data[[v]]
  data.frame(variable = v, n = length(x), mean = mean(x), sd = sd(x), min = min(x),
             Q1 = unname(quantile(x, .25)), median = median(x),
             Q3 = unname(quantile(x, .75)), max = max(x))
}))
write.csv(desc, file.path(OUT, "04_descriptive_statistics.csv"), row.names = FALSE)
print(desc, digits = 4)


# ---------------------- BAGIAN KODE 7 ----------------------
cor_all <- cor(data[vars])
write.csv(cor_all, file.path(OUT, "05_correlation_matrix.csv"))

cor_YX <- do.call(rbind, lapply(xvars, function(v){
  ct <- cor.test(data$Y_TB, data[[v]])
  data.frame(variable = v, r = unname(ct$estimate), p_value = ct$p.value)
}))
write.csv(cor_YX, file.path(OUT, "06_correlation_Y_X.csv"), row.names = FALSE)
print(cor_YX)

vif_awal <- car::vif(lm(Y_TB ~ X1_miskin + X2_sanitasi + X3_merokok + X4_hunian, data = data))
write.csv(data.frame(variable = names(vif_awal), VIF = as.numeric(vif_awal)),
          file.path(OUT, "07_VIF_awal.csv"), row.names = FALSE)
print(vif_awal) 

# ---------------------- BAGIAN KODE 8 ----------------------
data$Y_model <- if(Y_TRANSFORM == "log") log(data$Y_TB) else data$Y_TB
model_vars <- c("Y_model", xvars)
mdat <- data[, c("id", "wilayah", model_vars)]
if(STANDARDIZE) mdat[model_vars] <- lapply(mdat[model_vars], function(x) as.numeric(scale(x)))
formula_model <- Y_model ~ X1_miskin + X2_sanitasi + X3_merokok + X4_hunian

# ---------------------- BAGIAN KODE 9 ----------------------
OLS <- lm(formula_model, data = mdat)
save_txt(summary(OLS), "08_OLS_summary.txt")
write.csv(as.data.frame(coef(summary(OLS))), file.path(OUT, "09_OLS_coefficients.csv"))
n <- nrow(mdat)

save_txt(resettest(OLS, power = 2:3, type = "fitted"), "10_RESET.txt")
save_txt(shapiro.test(residuals(OLS)), "11_Shapiro_Wilk.txt")

# Uji Breusch-Pagan manual berdasarkan rumus matriks
bp_e <- as.numeric(residuals(OLS))
bp_n <- length(bp_e)
bp_sigma2 <- sum(bp_e^2) / bp_n
bp_Z <- model.matrix(OLS)[, -1, drop = FALSE]
bp_f <- (bp_e^2 / bp_sigma2) - 1

bp_ztz <- crossprod(bp_Z)
if (qr(bp_ztz)$rank < ncol(bp_ztz)) {
  stop("Matriks Z'Z singular pada uji BP manual.")
}
bp_stat <- as.numeric(0.5 * t(bp_f) %*% bp_Z %*%
                         solve(bp_ztz) %*% t(bp_Z) %*% bp_f)
bp_df <- ncol(bp_Z)
bp_p <- pchisq(bp_stat, df = bp_df, lower.tail = FALSE)
BP_MANUAL <- data.frame(
  statistic = bp_stat, df = bp_df,
  chi_square_critical = qchisq(0.95, df = bp_df),
  p_value = bp_p,
  decision = ifelse(bp_p < 0.05, "Tolak H0", "Gagal menolak H0"),
  method = "Breusch-Pagan manual berbasis rumus matriks"
)
write.csv(BP_MANUAL, file.path(OUT, "12_Breusch_Pagan_manual.csv"), row.names = FALSE)
save_txt(BP_MANUAL, "12_Breusch_Pagan_manual.txt")

vif_ols <- car::vif(OLS)
write.csv(data.frame(variable = names(vif_ols), VIF = as.numeric(vif_ols)),
          file.path(OUT, "13_VIF_OLS.csv"), row.names = FALSE)

cook <- cooks.distance(OLS)
write.csv(data.frame(id = mdat$id, wilayah = mdat$wilayah, Cook_D = cook) %>% arrange(desc(Cook_D)),
          file.path(OUT, "14_Cooks_distance.csv"), row.names = FALSE)

AICc_OLS <- aicc_ols(OLS)

# ---------------------- BAGIAN KODE 10 ----------------------
poly_proj <- st_transform(analysis_poly, 3857)
pts_geom  <- st_transform(
  st_point_on_surface(st_geometry(poly_proj)), 4326
)

xy <- st_coordinates(pts_geom)
colnames(xy) <- c("lon", "lat")

write.csv(
  data.frame(id = mdat$id, wilayah = mdat$wilayah, xy),
  file.path(OUT, "15_coordinates.csv"),
  row.names = FALSE
)

pts_sp <- SpatialPointsDataFrame(
  coords = xy,
  data = as.data.frame(mdat),
  proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs")
)

# Matriks bobot spasial KNN 8
knn_nb <- knn2nb(
  knearneigh(xy, k = K_MORAN, longlat = TRUE)
)

listw_knn <- nb2listw(
  knn_nb,
  style = "W",
  zero.policy = TRUE
)

# Ringkasan matriks bobot KNN 8
write.csv(
  data.frame(
    weights = "KNN-8",
    n = length(knn_nb),
    min_nb = min(card(knn_nb)),
    median_nb = median(card(knn_nb)),
    max_nb = max(card(knn_nb))
  ),
  file.path(OUT, "16_weights_summary.csv"),
  row.names = FALSE
)

# Analisis sensitivitas Moran's I untuk beberapa nilai KNN.
knn_sensitivity <- do.call(rbind, lapply(c(4, 8, 12), function(k) {
  nb_k <- knn2nb(knearneigh(xy, k = k, longlat = TRUE))
  lw_k <- nb2listw(nb_k, style = "W", zero.policy = TRUE)
  mt_k <- moran.test(mdat$Y_model, lw_k, randomisation = TRUE,
                     alternative = "two.sided", zero.policy = TRUE)
  data.frame(K = k, Moran_I = unname(mt_k$estimate[["Moran I statistic"]]),
             expectation = unname(mt_k$estimate[["Expectation"]]),
             variance = unname(mt_k$estimate[["Variance"]]),
             p_value = mt_k$p.value)
}))
write.csv(knn_sensitivity, file.path(OUT, "16a_Moran_sensitivity_KNN.csv"), row.names = FALSE)

# ---------------------- BAGIAN KODE 11 ----------------------
save_txt(
  moran.test(
    mdat$Y_model,
    listw_knn,
    randomisation = TRUE,
    alternative = "two.sided",
    zero.policy = TRUE
  ),
  "17_Moran_Y_KNN8.txt"
)

save_txt(
  lm.morantest(
    OLS,
    listw_knn,
    zero.policy = TRUE
  ),
  "19_Moran_OLS_residual_KNN8.txt"
)

# ---------------------- BAGIAN KODE 12 ----------------------
cat("\n===== GWR: PEMILIHAN KERNEL DAN BANDWIDTH =====\n")
kernels_to_run <- if(RUN_KERNEL_COMPARISON) KERNEL_CANDIDATES else KERNEL_DEFAULT

kernel_table <- do.call(rbind, lapply(kernels_to_run, function(k){
  cat("Kernel:", k, "\n")
  bw <- bw.gwr(formula_model, data = pts_sp, approach = "AICc", kernel = k,
               adaptive = ADAPTIVE, longlat = LONG_LAT)
  m  <- gwr.basic(formula_model, data = pts_sp, bw = bw, kernel = k,
                  adaptive = ADAPTIVE, longlat = LONG_LAT)
  data.frame(kernel = k, bandwidth = bw,
             AICc = get_stat(m$GW.diagnostic, "AICc"),
             R2   = get_stat(m$GW.diagnostic, c("gw.R2", "R2.val")),
             stringsAsFactors = FALSE)
}))
write.csv(kernel_table, file.path(OUT, "21_GWR_kernel_comparison.csv"), row.names = FALSE)
print(kernel_table)

best <- if(all(is.na(kernel_table$AICc))) 1 else which.min(kernel_table$AICc)
KERNEL  <- kernel_table$kernel[best]
BW_GWR  <- kernel_table$bandwidth[best]
cat("Kernel terpilih:", KERNEL, "| Bandwidth GWR:", BW_GWR, "| Mode:", if (ADAPTIVE) "adaptive" else "fixed", "\n")

GWR <- gwr.basic(formula_model, data = pts_sp, bw = BW_GWR, kernel = KERNEL,
                 adaptive = ADAPTIVE, longlat = LONG_LAT, F123.test = TRUE, cv = FALSE)
save_txt(GWR, "22_GWR_summary.txt")
# Simpan hasil uji F Leung (F1, F2, F3) jika tersedia dari GWmodel
save_txt(GWR$Ftest.res, "22a_GWR_Leung_F_tests.txt")
GWR_DF <- as.data.frame(GWR$SDF)
GWR_DF$id <- mdat$id; GWR_DF$wilayah <- mdat$wilayah
write.csv(GWR_DF, file.path(OUT, "23_GWR_local_results.csv"), row.names = FALSE)

# Diagnostik GWR: multikolinearitas lokal
GWR_COLLIN <- gwr.collin.diagno(formula_model, data = pts_sp, bw = BW_GWR, kernel = KERNEL,
                                adaptive = ADAPTIVE, longlat = LONG_LAT)
if(!is.null(GWR_COLLIN$SDF)){
  coll_df <- as.data.frame(GWR_COLLIN$SDF)
  write.csv(coll_df, file.path(OUT, "24_GWR_local_collinearity.csv"), row.names = FALSE)
  save_txt(summary(coll_df), "25_GWR_local_collinearity_summary.txt")   
  
  # Ringkasan condition number lokal, jika kolom tersedia
  cn_col <- find_col(coll_df, c("local_CN", "Local_CN"))
  if(!is.na(cn_col)){
    cn <- as.numeric(coll_df[[cn_col]])
    cn_valid <- is.finite(cn)
    if(any(cn_valid)){
      cn_summary <- data.frame(
        n_lokasi = sum(cn_valid),
        CN_min = min(cn[cn_valid]),
        CN_median = median(cn[cn_valid]),
        CN_max = max(cn[cn_valid]),
        n_CN_gt30 = sum(cn[cn_valid] > 30),
        persen_CN_gt30 = 100 * mean(cn[cn_valid] > 30)
      )
      write.csv(cn_summary, file.path(OUT, "25a_GWR_local_CN_summary.csv"), row.names = FALSE)
    }
  }
}

# Uji Monte Carlo (opsional) dengan validasi bandwidth dan penanganan error.
if (RUN_MONTE_CARLO) {
  dp_n <- nrow(pts_sp@data)
  bw_ok <- is.finite(BW_GWR) && length(BW_GWR) == 1 &&
    if (ADAPTIVE) BW_GWR < dp_n else BW_GWR > 0

  if (!bw_ok) {
    mc_note <- paste0(
      "Monte Carlo tidak dijalankan: bandwidth tidak memenuhi syarat. BW_GWR=",
      BW_GWR, "; jumlah observasi=", dp_n,
      ". Periksa bandwidth dan dokumentasi GWmodel."
    )
    writeLines(mc_note, file.path(OUT, "26_GWR_MonteCarlo.txt"))
    message(mc_note)
  } else {
    MC <- tryCatch(
      gwr.montecarlo(formula_model, data = pts_sp, nsim = MC_SIM,
                     kernel = KERNEL, adaptive = ADAPTIVE,
                     bw = BW_GWR, longlat = LONG_LAT),
      error = function(e) {
        msg <- paste("Uji Monte Carlo gagal:", conditionMessage(e))
        writeLines(msg, file.path(OUT, "26_GWR_MonteCarlo.txt"))
        message(msg)
        NULL
      }
    )
    if (!is.null(MC)) save_txt(MC, "26_GWR_MonteCarlo.txt")
  }
}

# Moran residual GWR
res_col_g <- find_col(GWR_DF, "residual")
if(!is.na(res_col_g)){
  save_txt(moran.test(GWR_DF[[res_col_g]], listw_knn, randomisation = TRUE,
                      alternative = "two.sided", zero.policy = TRUE), "27_Moran_GWR_residual.txt")
}

# ---------------------- BAGIAN KODE 13 ----------------------
cat("\n===== MGWR / BACKFITTING =====\n")
p_par <- length(xvars) + 1            

MGWR <- gwr.multiscale(
  formula            = formula_model,
  data               = pts_sp,
  kernel             = KERNEL,
  adaptive           = ADAPTIVE,
  criterion          = MGWR_CRITERION,
  approach           = MGWR_APPROACH,
  bws0               = rep(BW_GWR, p_par),     
  bw.seled           = rep(FALSE, p_par),      
  bws.reOpts         = MGWR_REOPT,
  max.iterations     = MGWR_MAX_ITER,
  threshold          = MGWR_THRESHOLD,
  predictor.centered = rep(TRUE, length(xvars)),
  hatmatrix          = TRUE,
  longlat            = LONG_LAT,
  verbose            = FALSE,
  nlower             = MGWR_NLOWER
)
save_txt(MGWR, "28_MGWR_summary.txt")
saveRDS(MGWR, file.path(OUT, "28_MGWR_object.rds"))   

# ---------------------- BAGIAN KODE 14 ----------------------
BW_FINAL <- MGWR$GW.arguments$bws
if(is.null(BW_FINAL)) BW_FINAL <- MGWR$bws.vars
if(is.null(BW_FINAL)) stop("Bandwidth tidak ditemukan. Jalankan str(MGWR$GW.arguments).")

BW_TABLE <- data.frame(
  parameter = c("Intercept", xvars),
  bandwidth = as.numeric(BW_FINAL)
)
BW_TABLE$persen_dari_n <- round(100 * BW_TABLE$bandwidth / n, 1)
write.csv(BW_TABLE, file.path(OUT, "29_MGWR_bandwidth.csv"), row.names = FALSE)
print(BW_TABLE)
if(length(unique(round(BW_TABLE$bandwidth, 3))) == 1){
  warning("Semua bandwidth MGWR sama. Periksa konvergensi/pengaturan bw.seled.")
}


# ---------------------- BAGIAN KODE 15 ----------------------
MGWR_DF <- as.data.frame(MGWR$SDF)
MGWR_DF$id <- mdat$id; MGWR_DF$wilayah <- mdat$wilayah
write.csv(MGWR_DF, file.path(OUT, "30_MGWR_all_local_results.csv"), row.names = FALSE)
writeLines(names(MGWR_DF), file.path(OUT, "31_MGWR_columns.txt"))

local_table <- data.frame(id = mdat$id, wilayah = mdat$wilayah)
sig_rows <- list()
for(v in xvars){
  c_col  <- find_col(MGWR_DF, v)
  se_col <- find_col(MGWR_DF, c(paste0(v, "_SE"), paste0(v, ".SE")))
  t_col  <- find_col(MGWR_DF, c(paste0(v, "_TV"), paste0(v, ".TV"), paste0(v, "_t")))
  if(is.na(c_col)) stop("Kolom koefisien ", v, " tidak ditemukan. Lihat 31_MGWR_columns.txt")
  
  b  <- MGWR_DF[[c_col]]
  tv <- if(!is.na(t_col)) MGWR_DF[[t_col]] else if(!is.na(se_col)) b / MGWR_DF[[se_col]] else rep(NA_real_, n)
  if(all(is.na(tv))) warning("t-value untuk ", v, " tidak tersedia; signifikansi tidak dapat dihitung.")
  
  # P-value dua sisi berbasis pendekatan normal untuk statistik t lokal
  p_raw <- 2 * pnorm(-abs(tv))
  p_raw[!is.finite(p_raw)] <- NA_real_
  
  sig <- abs(tv) > T_CRIT
  local_table[[paste0(v, "_coef")]] <- b
  local_table[[paste0(v, "_t")]] <- tv
  local_table[[paste0(v, "_p_raw")]] <- p_raw
  local_table[[paste0(v, "_sig")]] <- sig
  local_table[[paste0(v, "_coef_sig")]] <- ifelse(sig, b, NA_real_)
  
  sig_rows[[v]] <- data.frame(
    variable = v,
    bandwidth = BW_TABLE$bandwidth[BW_TABLE$parameter == v],
    coef_min = min(b, na.rm = TRUE), coef_Q1 = unname(quantile(b, .25, na.rm = TRUE)),
    coef_median = median(b, na.rm = TRUE), coef_mean = mean(b, na.rm = TRUE),
    coef_Q3 = unname(quantile(b, .75, na.rm = TRUE)), coef_max = max(b, na.rm = TRUE),
    n_sig_uncorrected = sum(sig, na.rm = TRUE),
    n_sig_positif_uncorrected = sum(sig & b > 0, na.rm = TRUE),
    n_sig_negatif_uncorrected = sum(sig & b < 0, na.rm = TRUE),
    persen_sig_uncorrected = round(100 * mean(sig, na.rm = TRUE), 1),
    n_sig_t1645 = sum(abs(tv) > 1.645, na.rm = TRUE),
    n_sig_t2576 = sum(abs(tv) > 2.576, na.rm = TRUE)
  )
}

# Koreksi Benjamini-Hochberg per variabel (FDR 5%)
# Setiap variabel dikoreksi atas 514 pengujian lokalnya masing-masing.
for(v in xvars){
  p_raw <- local_table[[paste0(v, "_p_raw")]]
  p_adj <- rep(NA_real_, length(p_raw))
  valid_p <- is.finite(p_raw)
  p_adj[valid_p] <- p.adjust(p_raw[valid_p], method = "BH")

  local_table[[paste0(v, "_p_BH")]] <- p_adj
  local_table[[paste0(v, "_sig_BH")]] <- !is.na(p_adj) & p_adj < 0.05
  local_table[[paste0(v, "_coef_sig_BH")]] <- ifelse(
    local_table[[paste0(v, "_sig_BH")]],
    local_table[[paste0(v, "_coef")]], NA_real_
  )
}

sig_table <- do.call(rbind, lapply(xvars, function(v){
  b <- local_table[[paste0(v, "_coef")]]
  tv <- local_table[[paste0(v, "_t")]]
  sig_bh <- local_table[[paste0(v, "_sig_BH")]]
  data.frame(
    variable = v,
    bandwidth = BW_TABLE$bandwidth[BW_TABLE$parameter == v],
    coef_min = min(b, na.rm = TRUE), coef_Q1 = unname(quantile(b, .25, na.rm = TRUE)),
    coef_median = median(b, na.rm = TRUE), coef_mean = mean(b, na.rm = TRUE),
    coef_Q3 = unname(quantile(b, .75, na.rm = TRUE)), coef_max = max(b, na.rm = TRUE),
    n_sig_BH = sum(sig_bh, na.rm = TRUE),
    n_sig_positif_BH = sum(sig_bh & b > 0, na.rm = TRUE),
    n_sig_negatif_BH = sum(sig_bh & b < 0, na.rm = TRUE),
    persen_sig_BH = round(100 * mean(sig_bh, na.rm = TRUE), 1),
    n_sig_t1645 = sum(abs(tv) > 1.645, na.rm = TRUE),
    n_sig_t2576 = sum(abs(tv) > 2.576, na.rm = TRUE)
  )
}))

write.csv(local_table, file.path(OUT, "32_MGWR_coef_t_signifikansi.csv"), row.names = FALSE)
write.csv(sig_table, file.path(OUT, "33_MGWR_ringkasan_koefisien.csv"), row.names = FALSE)
print(sig_table, digits = 3)

# ---------------------- BAGIAN KODE 16 ----------------------
res_col_m <- find_col(MGWR_DF, "residual")
fit_col_m <- find_col(MGWR_DF, c("yhat", "fitted"))
MGWR_RESID <- if(!is.na(res_col_m)) MGWR_DF[[res_col_m]] else
  if(!is.na(fit_col_m)) mdat$Y_model - MGWR_DF[[fit_col_m]] else rep(NA_real_, n)
GWR_RESID  <- if(!is.na(res_col_g)) GWR_DF[[res_col_g]] else rep(NA_real_, n)

if(all(is.finite(MGWR_RESID))){
  save_txt(moran.test(MGWR_RESID, listw_knn, randomisation = TRUE,
                      alternative = "two.sided", zero.policy = TRUE), "34_Moran_MGWR_residual.txt")
}

# ---------------------- BAGIAN KODE 17 ----------------------
GWR_DIAG  <- GWR$GW.diagnostic
MGWR_DIAG <- MGWR$GW.diagnostic
writeLines(capture.output(str(GWR_DIAG)),  file.path(OUT, "35_GWR_diag_structure.txt"))
writeLines(capture.output(str(MGWR_DIAG)), file.path(OUT, "36_MGWR_diag_structure.txt"))

tss <- sum((mdat$Y_model - mean(mdat$Y_model))^2)
rmse_fun <- function(r) {
  r <- as.numeric(r)
  if (!length(r) || any(!is.finite(r))) return(NA_real_)
  sqrt(mean(r^2))
}
rss_fun <- function(r) {
  r <- as.numeric(r)
  if (!length(r) || any(!is.finite(r))) return(NA_real_)
  sum(r^2)
}
res_ols <- as.numeric(residuals(OLS))
rss_ols <- rss_fun(res_ols)
rss_g <- rss_fun(GWR_RESID)
rss_m <- rss_fun(MGWR_RESID)

comparison <- data.frame(
  model = c("OLS", "GWR", "MGWR"),
  AICc  = c(AICc_OLS, get_stat(GWR_DIAG, "AICc"), get_stat(MGWR_DIAG, "AICc")),
  R2    = c(if(is.finite(rss_ols)) 1 - rss_ols/tss else NA_real_,
            if(is.finite(rss_g)) 1 - rss_g/tss else NA_real_,
            if(is.finite(rss_m)) 1 - rss_m/tss else NA_real_),
  Adj_R2 = c(summary(OLS)$adj.r.squared,
             get_stat(GWR_DIAG,  c("gwR2.adj", "R2adj")),
             get_stat(MGWR_DIAG, c("gwR2.adj", "R2adj"))),
  RSS   = c(rss_ols, rss_g, rss_m),
  RMSE_standardized_or_model_scale = c(rmse_fun(res_ols), rmse_fun(GWR_RESID), rmse_fun(MGWR_RESID))
)
writeLines(
  "Catatan: RMSE pada tabel berada pada skala respons model. Karena STANDARDIZE=TRUE, nilai ini berada pada skala Y terstandarisasi, bukan TB rate asli.",
  file.path(OUT, "37a_catatan_RMSE.txt")
)
write.csv(comparison, file.path(OUT, "37_perbandingan_model.csv"), row.names = FALSE)
print(comparison, digits = 4)

# RMSE pada skala TB rate asli dihitung dari fitted values yang tersedia.
fit_ols_original <- if (STANDARDIZE) {
  as.numeric(fitted(OLS)) * sd(data$Y_TB) + mean(data$Y_TB)
} else as.numeric(fitted(OLS))
fit_gwr_model <- mdat$Y_model - as.numeric(GWR_RESID)
fit_mgwr_model <- mdat$Y_model - as.numeric(MGWR_RESID)
if (Y_TRANSFORM == "raw") {
  fit_gwr_original <- if (STANDARDIZE) fit_gwr_model * sd(data$Y_TB) + mean(data$Y_TB) else fit_gwr_model
  fit_mgwr_original <- if (STANDARDIZE) fit_mgwr_model * sd(data$Y_TB) + mean(data$Y_TB) else fit_mgwr_model
  rmse_original <- data.frame(
    model = c("OLS", "GWR", "MGWR"),
    RMSE_TB_rate_original_scale = c(
      sqrt(mean((data$Y_TB - fit_ols_original)^2, na.rm = TRUE)),
      sqrt(mean((data$Y_TB - fit_gwr_original)^2, na.rm = TRUE)),
      sqrt(mean((data$Y_TB - fit_mgwr_original)^2, na.rm = TRUE))
    )
  )
  write.csv(rmse_original, file.path(OUT, "37b_RMSE_skala_TB_asli.csv"), row.names = FALSE)
} else {
  writeLines("RMSE skala TB asli tidak dihitung otomatis karena Y_TRANSFORM='log'; lakukan back-transform dengan mempertimbangkan bias retransformation.",
             file.path(OUT, "37b_RMSE_skala_TB_asli.txt"))
}

# ---------------------- BAGIAN KODE 18 ----------------------
poly_result <- left_join(analysis_poly, local_table, by = c("id", "wilayah"))
poly_result$resid_MGWR <- MGWR_RESID

st_write(poly_result, file.path(OUT, "38_MGWR_polygon_results.gpkg"),
         delete_dsn = TRUE, quiet = TRUE)

# ---------------------- BAGIAN KODE 19 ----------------------
cat("\n================ ANALISIS SELESAI ================\n")
cat("n =", n, "| Y:", Y_TRANSFORM, "| Standardisasi:", STANDARDIZE, "\n")
cat("Kernel:", KERNEL, "| Bandwidth GWR:", BW_GWR, "| Mode bandwidth:", if (ADAPTIVE) "adaptive" else "fixed", "\n")
print(BW_TABLE); print(comparison, digits = 4)
cat("Hasil di folder:", OUT, "\n")
writeLines(capture.output(sessionInfo()), file.path(OUT, "39_sessionInfo.txt"))



# ---------------------- BAGIAN KODE 20 ----------------------
# PETA KOMBINASI VARIABEL SIGNIFIKAN DAN PERSAMAAN LOKAL MGWR

if (length(BW_FINAL) != length(xvars) + 1L ||
    any(!is.finite(BW_FINAL)) || any(BW_FINAL <= 0)) {
  stop("Bandwidth MGWR tidak valid untuk menghitung R-square lokal evaluatif.")
}

BW_EVAL_R2 <- median(as.numeric(BW_FINAL[-1]))
if (!is.finite(BW_EVAL_R2) || BW_EVAL_R2 <= 0) {
  stop("Bandwidth evaluasi R-square lokal tidak valid.")
}

# Matriks jarak antartitik dalam kilometer.
dist_matrix_km <- sp::spDists(xy, longlat = TRUE)

# Bobot exponential fixed bandwidth: w_ij = exp(-d_ij / h_eval).
weight_matrix <- exp(-dist_matrix_km / BW_EVAL_R2)

y_eval <- as.numeric(mdat$Y_model)
yhat_eval <- as.numeric(MGWR_DF[[find_col(MGWR_DF, c("yhat", "fitted"))]])

if (length(yhat_eval) != n || any(!is.finite(y_eval)) ||
    any(!is.finite(yhat_eval))) {
  stop("Nilai aktual atau prediksi MGWR tidak lengkap untuk menghitung R-square lokal.")
}

local_r2_eval <- vapply(seq_len(n), function(i) {
  w <- weight_matrix[i, ]
  sw <- sum(w)
  if (!is.finite(sw) || sw <= 0) return(NA_real_)

  ybar_w <- sum(w * y_eval) / sw
  sse_w <- sum(w * (y_eval - yhat_eval)^2)
  sst_w <- sum(w * (y_eval - ybar_w)^2)

  if (!is.finite(sst_w) || sst_w <= 0) return(NA_real_)
  1 - sse_w / sst_w
}, numeric(1))

local_table$R2_lokal_evaluatif <- local_r2_eval
local_table$bandwidth_evaluasi_R2_km <- BW_EVAL_R2

if (all(is.na(local_table$R2_lokal_evaluatif))) {
  stop("Seluruh R-square lokal evaluatif bernilai NA. Periksa koordinat dan data.")
}

# Simpan dokumentasi definisi dan parameter R-square lokal evaluatif.
r2_note <- c(
  "R-square lokal evaluatif berbobot berbasis prediksi MGWR.",
  "Bukan Local R-square bawaan GWmodel::gwr.multiscale.",
  "Rumus: 1 - [sum_j w_ij (y_j - yhat_j)^2] / [sum_j w_ij (y_j - ybar_iw)^2].",
  paste0("Bobot: w_ij = exp(-d_ij / h); kernel exponential fixed."),
  paste0("Bandwidth evaluasi h = median bandwidth empat prediktor MGWR = ",
         signif(BW_EVAL_R2, 8), " km."),
  "Jarak dihitung menggunakan sp::spDists(..., longlat=TRUE), dalam kilometer.",
  "Y dan yhat berada pada skala respons model yang digunakan saat estimasi MGWR.",
  "Nilai dapat negatif jika kesalahan prediksi lokal berbobot melebihi variasi lokal berbobot.",
  "Interpretasikan sebagai ukuran kecocokan lokal evaluatif untuk pengurutan, bukan sebagai diagnostik Local R2 resmi GWmodel."
)
writeLines(r2_note, file.path(OUT, "40a_definisi_R2_lokal_evaluatif.txt"))
write.csv(
  data.frame(id = local_table$id, wilayah = local_table$wilayah,
             R2_lokal_evaluatif = local_table$R2_lokal_evaluatif,
             bandwidth_evaluasi_R2_km = BW_EVAL_R2),
  file.path(OUT, "40b_R2_lokal_evaluatif.csv"),
  row.names = FALSE
)

# Persamaan lokal: bintang diberikan jika koefisien signifikan setelah BH
intercept_col <- find_col(MGWR_DF, c("(Intercept)", "Intercept", "intercept"))
if (is.na(intercept_col)) stop("Kolom intersep MGWR tidak ditemukan.")
fmt_coef <- function(x) formatC(x, format = "f", digits = 4)
equations <- character(nrow(local_table))
for (i in seq_len(nrow(local_table))) {
  b0 <- as.numeric(MGWR_DF[[intercept_col]][i])
  terms <- vapply(seq_along(xvars), function(j) {
    v <- xvars[j]
    b <- as.numeric(local_table[[paste0(v, "_coef")]][i])
    sg <- isTRUE(local_table[[paste0(v, "_sig_BH")]][i])
    paste0(ifelse(b >= 0, "+", "-"),
           fmt_coef(abs(b)), ifelse(sg, "*", ""), " ", paste0("X", j))
  }, character(1))
  equations[i] <- paste0("Yhat = ", fmt_coef(b0), " ",
                         paste(terms, collapse = " "))
}
local_table$persamaan_MGWR <- equations

# Label kombinasi signifikan
sig_cols <- paste0(xvars, "_sig_BH")
local_table$kombinasi <- apply(local_table[, sig_cols, drop = FALSE], 1, function(z) {
  idx <- which(as.logical(z))
  if (!length(idx)) "Tidak ada variabel signifikan"
  else paste0("X", idx, collapse = ", ")
})

# Ringkasan jumlah wilayah untuk setiap kombinasi
kombinasi_table <- as.data.frame(table(local_table$kombinasi), stringsAsFactors = FALSE)
names(kombinasi_table) <- c("Kombinasi", "Jumlah_wilayah")
kombinasi_table$Persentase <- round(100 * kombinasi_table$Jumlah_wilayah / nrow(local_table), 2)
kombinasi_table <- kombinasi_table[order(-kombinasi_table$Jumlah_wilayah), ]
write.csv(kombinasi_table, file.path(OUT, "40_kombinasi_signifikansi_MGWR.csv"), row.names = FALSE)

# Tabel persamaan diurutkan berdasarkan R-square lokal evaluatif tertinggi.
persamaan_table <- local_table[, c(
  "id", "wilayah", "R2_lokal_evaluatif",
  "bandwidth_evaluasi_R2_km", "kombinasi", "persamaan_MGWR"
)]
persamaan_table <- persamaan_table[
  order(-persamaan_table$R2_lokal_evaluatif, na.last = TRUE),
]
write.csv(
  persamaan_table,
  file.path(OUT, "41_persamaan_lokal_MGWR_urut_R2.csv"),
  row.names = FALSE
)

# Peta kombinasi variabel signifikan
poly_result <- left_join(analysis_poly,
                         local_table[, c("id", "wilayah", "kombinasi")],
                         by = c("id", "wilayah"))
lev <- kombinasi_table$Kombinasi
poly_result$kombinasi <- factor(poly_result$kombinasi, levels = lev)

p_komb <- ggplot(poly_result) +
  geom_sf(aes(fill = kombinasi), color = "white", linewidth = 0.025) +
  scale_fill_viridis_d(
    option = "C",
    direction = 1,
    drop = FALSE
  ) +
  theme_void() +
  labs(title = "Peta Kombinasi Variabel Signifikan MGWR",
       subtitle = paste0("Kernel: ", KERNEL,
                         " | Koreksi BH per variabel (FDR 5%)"),
       fill = "Kombinasi variabel") +
  theme(plot.title = element_text(face = "bold"),
        plot.subtitle = element_text(size = 8),
        legend.position = "right")
ggsave(file.path(MAP, "Peta_Kombinasi_Variabel_Signifikan_MGWR.png"),
       p_komb, width = 12, height = 6, dpi = 300)

st_write(poly_result, file.path(OUT, "42_peta_kombinasi_MGWR.gpkg"),
         delete_dsn = TRUE, quiet = TRUE)

cat("\\nVisualisasi utama tersimpan di:", MAP, "\\n")
cat("Persamaan lokal terurut berdasarkan R-square MGWR tersimpan di:",
    file.path(OUT, "41_persamaan_lokal_MGWR_urut_R2.csv"), "\\n")
cat("Kombinasi signifikansi tersimpan di:",
    file.path(OUT, "40_kombinasi_signifikansi_MGWR.csv"), "\\n")

1 Hasil analisis

1.1 Statistik deskriptif

knitr::kable(desc, digits = 4, caption = "Statistik deskriptif variabel penelitian")

Statistik deskriptif variabel penelitian
variable	n	mean	sd	min	Q1	median	Q3	max
Y_TB	514	283.6167	185.4417	10.2809	168.1939	228.5771	349.6787	1478.435
X1_miskin	514	10.5813	6.9540	1.9000	5.9625	8.7450	12.7775	42.560
X2_sanitasi	514	82.7718	17.1565	0.0000	77.7625	88.1900	93.4650	100.000
X3_merokok	514	28.4714	5.5225	3.6400	25.2025	28.8350	32.3675	45.690
X4_hunian	514	66.3830	19.6847	0.0000	57.8975	70.3050	80.5850	98.410

1.2 Korelasi dan multikolinearitas

knitr::kable(cor_YX, digits = 4, caption = "Korelasi antara TB rate dan setiap prediktor")

Korelasi antara TB rate dan setiap prediktor
variable	r	p_value
X1_miskin	-0.0784	0.0758
X2_sanitasi	0.0294	0.5054
X3_merokok	-0.1014	0.0215
X4_hunian	0.0034	0.9385

knitr::kable(data.frame(Variabel = names(vif_awal), VIF = as.numeric(vif_awal)),
             digits = 4, caption = "VIF awal prediktor")

VIF awal prediktor
Variabel	VIF
X1_miskin	1.5466
X2_sanitasi	3.1967
X3_merokok	1.0191
X4_hunian	2.7384

knitr::kable(data.frame(Variabel = names(vif_ols), VIF = as.numeric(vif_ols)),
             digits = 4, caption = "VIF model OLS")

VIF model OLS
Variabel	VIF
X1_miskin	1.5466
X2_sanitasi	3.1967
X3_merokok	1.0191
X4_hunian	2.7384

1.3 Diagnostik model OLS

1.3.1 Uji RESET

cat(readLines(file.path(OUT, "10_RESET.txt")), sep = "\n")

#> 
#>  RESET test
#> 
#> data:  OLS
#> RESET = 2.7426, df1 = 2, df2 = 507, p-value = 0.06536

1.3.2 Uji normalitas Shapiro-Wilk

cat(readLines(file.path(OUT, "11_Shapiro_Wilk.txt")), sep = "\n")

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(OLS)
#> W = 0.81156, p-value < 0.00000000000000022

1.3.3 Uji heteroskedastisitas Breusch-Pagan manual

# Tabel khusus untuk tampilan: p-value ditampilkan dalam notasi ilmiah
BP_MANUAL_TABEL <- BP_MANUAL

BP_MANUAL_TABEL$p_value <- format.pval(
  BP_MANUAL$p_value,
  digits = 5,
  eps = 1e-16
)

knitr::kable(
  BP_MANUAL_TABEL,
  digits = 5,
  caption = "Hasil uji Breusch-Pagan manual",
  align = "c"
)

Hasil uji Breusch-Pagan manual
statistic	df	chi_square_critical	p_value	decision	method
56.98872	4	9.48773	0.000000000012439	Tolak H0	Breusch-Pagan manual berbasis rumus matriks

1.4 Hasil regresi OLS, GWR, dan MGWR

knitr::kable(comparison, digits = 4, caption = "Perbandingan kinerja model OLS, GWR, dan MGWR")

Perbandingan kinerja model OLS, GWR, dan MGWR
model	AICc	R2	Adj_R2	RSS	RMSE_standardized_or_model_scale
OLS	1459.176	0.0205	0.0128	502.4724	0.9887
GWR	1259.797	0.5186	0.3850	246.9413	0.6931
MGWR	1231.322	0.5458	0.4157	232.9917	0.6733

knitr::kable(BW_TABLE, digits = 3, caption = "Bandwidth setiap komponen MGWR")

Bandwidth setiap komponen MGWR
parameter	bandwidth	persen_dari_n
Intercept	121.365	23.6
X1_miskin	262.495	51.1
X2_sanitasi	5183.968	1008.6
X3_merokok	757.559	147.4
X4_hunian	338.314	65.8

1.5 Ringkasan koefisien lokal dan signifikansi MGWR

knitr::kable(sig_table, digits = 4, caption = "Ringkasan koefisien dan signifikansi lokal setelah koreksi BH")

Ringkasan koefisien dan signifikansi lokal setelah koreksi BH
variable	bandwidth	coef_min	coef_Q1	coef_median	coef_mean	coef_Q3	coef_max	n_sig_BH	n_sig_positif_BH	n_sig_negatif_BH	persen_sig_BH	n_sig_t1645	n_sig_t2576
X1_miskin	262.4950	-0.7943	-0.5371	-0.3816	-0.4227	-0.3143	-0.0849	471	0	471	91.6	482	448
X2_sanitasi	5183.9682	-0.1786	-0.1766	-0.1729	-0.1719	-0.1686	-0.1596	514	0	514	100.0	514	0
X3_merokok	757.5595	-0.3560	-0.3303	-0.2783	-0.2552	-0.2063	0.0352	471	0	471	91.6	474	466
X4_hunian	338.3136	-0.1379	-0.0331	0.0800	0.1240	0.1919	0.6796	92	92	0	17.9	153	93

1.6 Peta kombinasi variabel signifikan MGWR

Peta berikut menggunakan klasifikasi kombinasi variabel yang signifikan setelah koreksi Benjamini-Hochberg per variabel pada taraf FDR 5%.

print(p_komb)

1.7 Ringkasan kombinasi variabel signifikan

knitr::kable(kombinasi_table, caption = "Jumlah dan persentase wilayah menurut kombinasi variabel signifikan")

Jumlah dan persentase wilayah menurut kombinasi variabel signifikan
Kombinasi	Jumlah_wilayah	Persentase
X1, X2, X3	379	73.74
X1, X2, X3, X4	49	9.53
X1, X2, X4	43	8.37
X2, X3	43	8.37

1.8 Persamaan lokal MGWR: urutan R-square tertinggi hingga terendah

DT::datatable(
  persamaan_table,
  rownames = FALSE,
  filter = "top",
  extensions = "Buttons",
  options = list(
    pageLength = 10,
    lengthMenu = c(10, 25, 50, 100, 514),
    scrollX = TRUE,
    dom = "Bfrtip",
    buttons = c("copy", "csv", "excel")
  ),
  caption = "Seluruh persamaan lokal MGWR, diurutkan berdasarkan R-square lokal evaluatif dari tertinggi ke terendah"
)

1.8.1 Sepuluh persamaan dengan R-square tertinggi

knitr::kable(head(persamaan_table, 10), digits = 5,
             caption = "Sepuluh wilayah dengan R-square lokal evaluatif tertinggi")

Sepuluh wilayah dengan R-square lokal evaluatif tertinggi
	id	wilayah	R2_lokal_evaluatif	bandwidth_evaluasi_R2_km	kombinasi	persamaan_MGWR
486	486	Kab. Merauke	0.66180	547.9365	X1, X2, X4	Yhat = 3.4466 -0.2931* X1 -0.1624* X2 +0.0293 X3 +0.5289* X4
495	495	Kab. Boven Digoel	0.66021	547.9365	X1, X2, X4	Yhat = 3.7170 -0.3328* X1 -0.1616* X2 +0.0352 X3 +0.5553* X4
497	497	Kab. Asmat	0.65432	547.9365	X1, X2, X4	Yhat = 2.9547 -0.3833* X1 -0.1611* X2 +0.0188 X3 +0.6009* X4
506	506	Kab. Nduga	0.65327	547.9365	X1, X2, X4	Yhat = 2.1416 -0.4478* X1 -0.1598* X2 +0.0032 X3 +0.6644* X4
496	496	Kab. Mappi	0.65145	547.9365	X1, X2, X4	Yhat = 4.2418 -0.2921* X1 -0.1622* X2 +0.0289 X3 +0.5280* X4
507	507	Kab. Lanny Jaya	0.65051	547.9365	X1, X2, X4	Yhat = 2.0680 -0.4553* X1 -0.1596* X2 -0.0006 X3 +0.6690* X4
494	494	Kab. Mimika	0.65045	547.9365	X1, X2, X4	Yhat = 2.3991 -0.4238* X1 -0.1603* X2 -0.0120 X3 +0.6708* X4
498	498	Kab. Yahukimo	0.65019	547.9365	X1, X2, X4	Yhat = 2.1914 -0.4489* X1 -0.1600* X2 +0.0151 X3 +0.6422* X4
487	487	Kab. Jayawijaya	0.64943	547.9365	X1, X2, X4	Yhat = 2.0871 -0.4570* X1 -0.1597* X2 +0.0078 X3 +0.6556* X4
499	499	Kab. Pegunungan Bintang	0.64731	547.9365	X1, X2, X4	Yhat = 2.1595 -0.4483* X1 -0.1601* X2 +0.0239 X3 +0.6310* X4

Analisis MGWR Fixed Bandwidth pada TB Rate 514 Kabupaten/Kota di Indonesia

Analisis Data Spasial

04 November 2026

1 Hasil analisis

1.1 Statistik deskriptif

1.2 Korelasi dan multikolinearitas

1.3 Diagnostik model OLS

1.3.1 Uji RESET

1.3.2 Uji normalitas Shapiro-Wilk

1.3.3 Uji heteroskedastisitas Breusch-Pagan manual

1.4 Hasil regresi OLS, GWR, dan MGWR

1.5 Ringkasan koefisien lokal dan signifikansi MGWR

1.6 Peta kombinasi variabel signifikan MGWR

1.7 Ringkasan kombinasi variabel signifikan

1.8 Persamaan lokal MGWR: urutan R-square tertinggi hingga terendah

1.8.1 Sepuluh persamaan dengan R-square tertinggi