---
title: "Artikel MKMI"
author: "SILFIANA LIS SETYOWATI (G1501231093)"
date: "2025-03-13"
output: html_document
---

```{r}
# Baca file Excel dan ambil data dari sheet 1
library(readxl)
library(dplyr)
library(tidyr)
data <- read_excel("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/Bimbingan/Ringkasan Analisis Empiris.xlsx", sheet = "Data Fix")

head(data)
```

```{r}
summary(data)
```

```{r}
# Load library
library(tidyr)
library(ggplot2)

# Baca data dari file Excel
data <- read_excel("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/Bimbingan/Ringkasan Analisis Empiris.xlsx", 
                   sheet = "Data Fix")

# Seleksi hanya kolom numerik
numeric_data <- data %>% select(where(is.numeric))

# Ubah data ke format long agar bisa dibuat boxplot dengan ggplot
data_long <- numeric_data %>% 
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")

# Buat boxplot
ggplot(data_long, aes(x = Variable, y = Value)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Boxplot of All Numeric Variables", x = "Variables", y = "Values")

```

```{r}
library(ggplot2)
# Pisahkan data berdasarkan wilayah
data_kabupaten <- data %>% filter(Lokasi == "Kab")
data_kota <- data %>% filter(Lokasi == "Kota")
summary(data_kabupaten)

# Pilih hanya kolom numerik
numeric_kabupaten <- data_kabupaten %>% select(where(is.numeric))
numeric_kota <- data_kota %>% select(where(is.numeric))

# Ubah data ke format long agar bisa dibuat boxplot
data_long_kabupaten <- numeric_kabupaten %>% 
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")

data_long_kota <- numeric_kota %>% 
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")

# Buat boxplot untuk Kabupaten
ggplot(data_long_kabupaten, aes(x = Variable, y = Value)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Boxplot Kabupaten", x = "Variables", y = "Values")

# Buat boxplot untuk Kota
ggplot(data_long_kota, aes(x = Variable, y = Value)) +
  geom_boxplot(fill = "tomato", color = "black") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Boxplot Kota", x = "Variables", y = "Values")

```

Klasifikasi WHO

```{r}
library(ggplot2)
library(dplyr)

# Buat data ringkasan berdasarkan jumlah kategori WHO Stunting
data_summary <- data.frame(
  WHO_Stunting = c("High", "Low", "Moderate", "Very High"),
  Count = c(48, 1, 29, 7)
) %>%
  mutate(percentage = Count / sum(Count) * 100) %>%
  mutate(WHO_Stunting = factor(WHO_Stunting, levels = c("Low", "Moderate", "High", "Very High")))

# Membuat pie chart dengan angka dan persentase
ggplot(data_summary, aes(x = "", y = percentage, fill = WHO_Stunting)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar(theta = "y") +
  theme_void() +
  scale_fill_manual(values = c("Low" = "green", "Moderate" = "yellow", "High" = "darkorange", "Very High" = "red")) +
  labs(title = "Pie Chart WHO Stunting di Kabupaten", fill = "WHO Stunting") +
  theme(axis.text.x = element_blank()) +
  geom_text(aes(label = paste(Count, "\n", round(percentage, 1), "%")), 
            position = position_stack(vjust = 0.5), color = "black", size = 2, fontface="bold") +  # Kecilkan ukuran teks
  theme(legend.position = "right")  # Urutkan legenda

```

```{r}
summary(data_kabupaten)
summary(data_kota)
```

```{r}
library(sf)
shp_banten <- st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]Banten_Kab/Banten_ADMIN_BPS.shp")
shp_diy <- st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]DIY_Kab/Daerah_Istimewa_Yogyakarta_ADMIN_BPS.shp")
shp_dki <- st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]DKI_Kab/Dki_Jakarta_ADMIN_BPS.shp")
shp_jabar <- st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]Jawa_Barat_Kab/Jawa_Barat_ADMIN_BPS.shp")
shp_jateng <- st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]Jawa_Tengah_Kab/Jawa_Tengah_ADMIN_BPS.shp")
shp_jatim=st_read("/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/[geosai.my.id]Jawa_Timur_Kab/Jawa_Timur_ADMIN_BPS.shp")

shp_pulau_jawa <- rbind(shp_banten,shp_dki,shp_diy,shp_jabar,shp_jateng,shp_jatim)
# Filter data untuk menghapus nama tertentu
shp_pulau_jawa<- shp_pulau_jawa %>%
  filter(!Kabupaten %in% c("Hutan", "Waduk Cirata", "Wadung Kedungombo"))

# Simpan kembali file tanpa data yang dihapus
st_write(shp_pulau_jawa, "/Users/bmn/Documents/STATISTIKA DAN SAINS DATA/THESIS/SHP Pulau Jawa/pulau jawa gab.shp", delete_dsn = TRUE)
```

```{r}
# Gabungkan data dengan shapefile
shp_pulau_jawa_merged <- shp_pulau_jawa %>%
  left_join(data, by = c("Kabupaten" = "NamaKab/Kota"))
shp_pulau_jawa_merged <- shp_pulau_jawa_merged %>%
  select(-Provinsi)
shp_pulau_jawa_merged <- na.omit(shp_pulau_jawa_merged)
shp_pulau_jawa_merged
summary(shp_pulau_jawa_merged)
```

```{r}
ggplot(data = shp_pulau_jawa_merged) +
  geom_sf(aes(fill = Y), color = "black", size = 0.001) + 
  scale_fill_steps(
    low = "yellow", 
    high = "red", 
    n.breaks = 5,  # Jumlah rentang dalam legenda
    na.value = "grey50"
  ) +
  labs(
    title = "Kondisi Stunting di Pulau Jawa",
    fill = "Prevalence (%)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    legend.position = "right"
  )
```
```{r}
library(ggplot2)
library(sf)
library(dplyr)

# Ubah WHO Stunting menjadi angka (1 = Low, 2 = Moderate, 3 = High, 4 = Very High)
shp_pulau_jawa_merged <- shp_pulau_jawa_merged %>%
  mutate(WHO_Stunting_numeric = case_when(
    `WHO Stunting` == "Low" ~ 1,
    `WHO Stunting` == "Moderate" ~ 2,
    `WHO Stunting` == "High" ~ 3,
    `WHO Stunting` == "Very High" ~ 4,
    TRUE ~ NA_real_  # Untuk mengatasi jika ada nilai yang tidak terdeteksi
  ))

# Periksa hasil untuk memastikan bahwa nilai sudah diubah dengan benar
table(shp_pulau_jawa_merged$WHO_Stunting_numeric)

# Membuat peta berdasarkan kategori WHO Stunting
ggplot(shp_pulau_jawa_merged) +
  geom_sf(aes(fill = factor(WHO_Stunting_numeric))) +  # Memetakan berdasarkan angka WHO Stunting
  scale_fill_manual(values = c("1" = "lightgreen", "2" = "yellow", "3" = "orange", "4" = "red"),
                    labels = c("1" = "Low", "2" = "Moderate", "3" = "High", "4" = "Very High")) +  # Mengatur warna dan label
  labs(title = "Pemetaan WHO Stunting di Kabupaten Pulau Jawa", fill = "WHO Stunting") +
  theme_minimal() +
  theme(legend.position = "right")  # Menambahkan posisi legenda
```

Untuk mendeteksi perbedaan rata-rata dua sampel (Kabupaten vs Kota) dalam dataset, kita bisa menggunakan uji t (t-test) jika asumsi distribusi normal terpenuhi atau uji Mann-Whitney jika data tidak normal.

```{r}
# Pisahkan data Kabupaten dan Kota untuk variabel Y
data_kabupaten_Y <- data %>% filter(Lokasi == "Kab") %>% pull(Y)
data_kota_Y <- data %>% filter(Lokasi == "Kota") %>% pull(Y)

# Uji normalitas dengan Shapiro-Wilk
p_kab <- shapiro.test(data_kabupaten_Y)$p.value
p_kota <- shapiro.test(data_kota_Y)$p.value

cat("\n===== UJI NORMALITAS =====\n")
cat("Shapiro-Wilk p-value (Kabupaten):", p_kab, "\n")
cat("Shapiro-Wilk p-value (Kota):", p_kota, "\n")

# Jika data normal, gunakan uji t, jika tidak normal gunakan Mann-Whitney
if (p_kab > 0.05 & p_kota > 0.05) {
  # Uji t
  test_result <- t.test(data_kabupaten_Y, data_kota_Y, var.equal = FALSE)
  cat("\n===== UJI T (t-test) =====\n")
} else {
  # Uji Mann-Whitney (Wilcoxon Rank Sum Test)
  test_result <- wilcox.test(data_kabupaten_Y, data_kota_Y)
  cat("\n===== UJI MANN-WHITNEY =====\n")
}

# Tampilkan hasil uji
print(test_result)
```

```{r}
library(car)
# Pilih hanya kolom numerik yang digunakan dalam regresi
predictors <- c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")

# Cek VIF untuk Kabupaten
model_kab <- lm(Y ~ ., data = data_kabupaten %>% select(Y, all_of(predictors)))
vif_kab <- vif(model_kab)
print(vif_kab)
library(lmtest)
bptest(model_kab, studentize = FALSE)

# Cek VIF untuk Kota
model_kota <- lm(Y ~ ., data = data_kota %>% select(Y, all_of(predictors)))
vif_kota <- vif(model_kota)
print(vif_kota)
bptest(model_kota)

summary(model_kab)
summary(model_kota, studentize= FALSE)

aic_kab <- AIC(model_kab)
aic_kota <- AIC(model_kota)
aic_kab
aic_kota
```

MODEL GWR KAB

```{r}
library(GWmodel)
shp_kabupaten <- shp_pulau_jawa_merged %>%
filter(Lokasi == "Kab") %>%
select(Kabupaten, everything())

# Tentukan formula regresi
formula_gwr <- Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 + X12

# Konversi sf ke Spatial untuk GWModel
shp_kabupaten_sp <- as(shp_kabupaten, "Spatial")

# Tentukan bandwidth optimal
set.seed(123)
bw_gwr <- bw.gwr(formula_gwr, 
                 data = shp_kabupaten_sp, 
                 adaptive = TRUE,  # Adaptive kernel
                 kernel = "gaussian", 
                 approach = "CV")
print(bw_gwr)

# Jalankan model GWR
gwr_model <- gwr.basic(formula_gwr, 
                        data = shp_kabupaten_sp, 
                        bw = bw_gwr, 
                        adaptive = TRUE, 
                        kernel = "gaussian")

# Lihat hasilnya
print(gwr_model)

```

```{r}
# Konversi gwr_model$SDF ke data frame
coef_gwr <- as.data.frame(gwr_model$SDF)

# Buat list untuk menyimpan p-value
p_values <- list()

# Dapatkan daftar variabel berdasarkan kolom yang memiliki "_TV"
tv_vars <- grep("_TV$", colnames(coef_gwr), value = TRUE)

# Hitung derajat kebebasan (df)
n <- nrow(coef_gwr)  # Jumlah wilayah (Kabupaten)
k <- length(tv_vars)  # Jumlah variabel prediktor dalam model
df <- n - k  # Derajat kebebasan

# Loop untuk menghitung p-value dari t-statistik (TV)
for (var in tv_vars) {
  var_name <- gsub("_TV", "", var)  # Ambil nama variabel asli tanpa _TV
  t_stat <- coef_gwr[[var]]  # Ambil t-statistik dari kolom _TV
  p_values[[var_name]] <- 2 * (1 - pt(abs(t_stat), df))  # Hitung p-value
}

# Konversi hasil ke dataframe
p_values_df <- as.data.frame(p_values)

# Tampilkan ringkasan p-value
summary(p_values_df)
```

```{r}
library(sf)
library(ggplot2)
library(dplyr)
library(gridExtra)

# Konversi hasil GWR ke sf dengan geometri dari shp_kabupaten
shp_kabupaten <- st_as_sf(as.data.frame(gwr_model$SDF), geometry = st_geometry(shp_kabupaten))

# Daftar variabel koefisien yang akan dipetakan
vars_to_plot <- c("X1", "X2", "X4", "X6", "X11")

# Daftar variabel asli dan nama yang akan digunakan dalam visualisasi
var_labels <- c(
  "X1"  = "Complete Immunization",
  "X2"  = "Growth Monitoring",
  "X4"  = "ARI",
  "X6"  = "Diarrhoea",
  "X11" = "Poverty Prevalence"
)

# Loop untuk membuat peta tiap variabel dengan nama yang diperbarui
plots <- list()
for (var in names(var_labels)) {
  p <- ggplot(data = shp_kabupaten) +
    geom_sf(aes_string(fill = var), color = "black", size = 10) +  # Ukuran garis lebih tipis
    scale_fill_steps(
      low = "red",   # Warna untuk koefisien negatif
      high = "blue",  # Warna untuk koefisien positif
      n.breaks = 5,
      na.value = "white"
    ) +
    labs(
      title = paste( var_labels[var]),
      fill = "coefficients"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),  # Judul lebih besar
      legend.position = "right",
      legend.title = element_text(size = 8, face = "bold"),  # Ukuran judul legend lebih kecil
      legend.text = element_text(size = 7),  # Ukuran teks legend lebih kecil
      legend.key.size = unit(0.4, "cm"),
      plot.margin = margin(5, 5, 5, 5)  # Memberikan ruang ekstra antar peta
    ) +
    coord_sf(expand = TRUE)  # Memberikan ruang ekstra di sekitar peta
  
  plots[[var]] <- p
}

# Tampilkan semua peta dalam grid (3 kolom agar lebih proporsional)
do.call(grid.arrange, c(plots, ncol = 2))

```

```{r}
# Konversi hasil GWR ke data frame
coef_gwr <- as.data.frame(gwr_model$SDF)
coef_gwr <- coef_gwr %>%
  mutate(Kabupaten = shp_kabupaten_sp@data[["Kabupaten"]])

# Buat list untuk menyimpan p-value
p_values <- list()

# Dapatkan daftar variabel berdasarkan kolom yang memiliki "_TV" (T-statistic Value)
tv_vars <- grep("_TV$", colnames(coef_gwr), value = TRUE)

# Hitung derajat kebebasan (df)
n <- nrow(coef_gwr)  # Jumlah wilayah (Kabupaten)
k <- length(tv_vars)  # Jumlah variabel prediktor dalam model
df <- n - k  # Derajat kebebasan

for (var in tv_vars) {
  var_name <- gsub("_TV", "", var)  # Ambil nama variabel asli tanpa _TV
  t_stat <- coef_gwr[[var]]  # Ambil nilai t-statistik dari kolom _TV
  p_values[[var_name]] <- 2 * (1 - pt(abs(t_stat), df))  # Hitung p-value
}
p_values_df <- as.data.frame(p_values)
p_values_df$Kabupaten <- coef_gwr$Kabupaten  # Tambahkan kolom nama Kabupaten
significant_vars <- p_values_df %>%
  mutate(across(-Kabupaten, ~ ifelse(. < 0.05, "Significant", "Not Significant")))
significant_only <- significant_vars %>%
  mutate(Significant_Variables = apply(select(., -Kabupaten), 1, function(row) {
    paste(names(row)[row == "Significant"], collapse = ", ")
  })) %>%
  select(Kabupaten, Significant_Variables)
print(significant_only)
```

```{r}
# Hitung jumlah kelompok unik berdasarkan pola variabel signifikan
unique_groups <- significant_only %>%
  distinct(Significant_Variables)  # Ambil kombinasi unik

# Tampilkan jumlah kelompok
num_groups <- nrow(unique_groups)
print(paste("Jumlah kelompok unik berdasarkan pola signifikansi:", num_groups))

# Tampilkan daftar kelompok unik
print(unique_groups)
```

```{r}

shp_kabupaten$Nama_Kabupaten <- coef_gwr$Kabupaten

library(stringr)

significant_only <- significant_only %>%
  mutate(Significant_Variables = str_remove_all(Significant_Variables, "Intercept,?\\s*"))  # Hapus "Intercept"

significant_only <- significant_only %>%
  mutate(Group = case_when(
    Significant_Variables == "X1, X2, X4, X11" ~ 1,
    Significant_Variables == "X1, X2, X4, X6, X11" ~ 2,
    Significant_Variables == "X2, X4, X11" ~ 3,
    Significant_Variables == "X2, X4, X6, X11" ~ 4,
    TRUE ~ NA_real_  # Jika pola tidak cocok, dibiarkan NA
  ))

library(sf)

# Gabungkan data kelompok dengan shapefile berdasarkan nama Kabupaten
shp_kabupaten_grouped <- shp_kabupaten %>%
  left_join(significant_only, by = c("Nama_Kabupaten" = "Kabupaten"))

library(ggplot2)

ggplot(data = shp_kabupaten_grouped) +
  geom_sf(aes(fill = factor(Group)), color = "black", size = 0.3) +
  scale_fill_manual(
    values = c("1" = "red", "2" = "blue", "3" = "green", "4" = "purple"),  # Warna untuk setiap kelompok
    na.value = "grey80"  # Warna untuk wilayah tanpa data
  ) +
  labs(
    title = "Peta Pengelompokan Kabupaten Berdasarkan Variabel Signifikan",
    fill = "Group"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
    legend.position = "right",
    legend.title = element_text(size = 8, face = "bold"),
    legend.text = element_text(size = 7),
    legend.key.size = unit(0.5, "cm"),
    plot.margin = margin(10, 10, 10, 10)
  ) +
  coord_sf(expand = TRUE)
```

```{r}
print(significant_only)
```

Cek Residual GWR

```{r}
library(sf)

# Hilangkan kolom geometri dari shp_kabupaten
data_no_geom <- st_drop_geometry(shp_kabupaten)

# Konversi menjadi matriks (termasuk Intercept)
X <- as.matrix(cbind(1, data_no_geom[, c("X1","X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")]))

# Konversi sf ke Spatial
shp_kabupaten_sp <- as(shp_kabupaten, "Spatial")

# Koordinat dari shp_kabupaten
coords <- as.matrix(coordinates(shp_kabupaten_sp))

# Ambil bandwidth dari model GWR
bandwidth <- gwr_model$GW.arguments$bw



#KERNEL dan GAUSIIAN ADAPTIVE
# Fungsi Kernel Gaussian
gaussian_kernel <- function(d, bw) {
  exp(- (d^2) / (bw^2))
}

# Fungsi untuk menghitung bobot
compute_weights <- function(coords, target_idx, bw) {
  dists <- sqrt(rowSums((coords - coords[target_idx, ])^2))
  weights <- gaussian_kernel(dists, bw)
  if (any(is.na(weights))) stop("Weights contain NA values.")
  return(weights)
}

#HAT MATRIX
# Inisialisasi Hat Matrix
hat_matrix <- matrix(0, nrow = nrow(X), ncol = nrow(X))

for (i in 1:nrow(X)) {
  print(paste("Processing location:", i))
  
  # Matriks bobot
  W <- diag(compute_weights(coords, i, bandwidth))
  
  # Perhitungan XtWX
  XtWX <- t(X) %*% W %*% X
  if (det(XtWX) == 0) stop("XtWX is singular. Check your X or W matrices.")
  XtWX_inv <- solve(XtWX)
  
  # Simpan leverage
  hat_matrix[i, ] <- diag(X %*% XtWX_inv %*% t(X) %*% W)
}

# Simpan leverage dalam vektor
leverage <- diag(hat_matrix)
leverage

#COOK DISTANCE
# Ambil residual dari model GWR
residuals <- gwr_model$SDF$residual

# Mean Squared Error (MSE)
MSE <- mean(residuals^2, na.rm = TRUE)

# Jumlah parameter dalam model
p <- ncol(X)  # Termasuk kolom intercept

# Hitung Cook's Distance
cook_distance <- (residuals^2 / (p * MSE)) * (leverage / (1 - leverage)^2)

# Hitung jumlah pengamatan (n)
n <- nrow(data_no_geom)

# Cut-off untuk Cook's Distance
cutoff <- 4 / n

# Identifikasi pencilan (outliers)
outliers <- which(cook_distance > cutoff)

# Simpan hasil pencilan dalam dataframe
outlier_data <- data_no_geom[outliers, ]
outlier_data$cook_distance <- cook_distance[outliers]

# Cetak hasil outlier
print(outlier_data)


#VISUALISASI
plot(cook_distance, 
     main = "Cook's Distance with Outliers", 
     xlab = "Observasi", 
     ylab = "Cook's Distance", 
     type = "h")

abline(h = cutoff, col = "red", lty = 2)  # Garis cut-off
points(outliers, cook_distance[outliers], col = "blue", pch = 19)  # Tandai pencilan
```

```{r}
library(sf)
library(ggplot2)

# Pastikan cook_distance telah dihitung dengan benar
summary(cook_distance)

# Gabungkan Cook's Distance dengan shp_kabupaten
shp_kabupaten_CookDistance <- cbind(shp_kabupaten, cook_distance)

# Tambahkan flag outlier berdasarkan cutoff
cutoff <- 4 / nrow(shp_kabupaten)
shp_kabupaten_CookDistance$outlier_flag <- ifelse(shp_kabupaten_CookDistance$cook_distance > cutoff, "Outlier", "Not Outlier")

# Konversi kembali ke sf (jika perlu)
shp_kabupaten_CookDistance <- st_as_sf(shp_kabupaten_CookDistance)

ggplot(data = shp_kabupaten_CookDistance) +
  geom_sf(aes(fill = outlier_flag), color = "black", size = 0.2) +  # Garis batas wilayah
  scale_fill_manual(values = c("Outlier" = "red", "Not Outlier" = "yellow")) +  # Warna: merah untuk outlier, kuning untuk tidak
  labs(title = "Outlier Map Based on Cook's Distance",
       fill = "Category") +
  theme_minimal()
```

Penduga dengan model M Robust GWR Using M Estimator

```{r}
shp_kabupaten_RGWRM <- shp_pulau_jawa_merged %>%
filter(Lokasi == "Kab") %>%
select(Kabupaten, everything())
shp_kabupaten_RGWRM <- st_make_valid(shp_kabupaten_RGWRM)

# Menjalankan Robust GWR
RGWRM <- gwr.robust(
  formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 + X12,
  data = shp_kabupaten_RGWRM,  # Data dalam format Spatial
  bw = bandwidth,          # Bandwidth yang telah ditentukan
  kernel = "bisquare",     # Kernel Bisquare untuk pembobotan spasial
  adaptive = TRUE,         # Menggunakan adaptive bandwidth
  maxiter = 20,            # Iterasi maksimum untuk konvergensi
  cut1 = 2,                # Ambang batas residual awal
  cut2 = 3                 # Ambang batas untuk filtering iteratif
)

# Tampilkan ringkasan hasil model
summary(RGWRM)

# Ekstrak koefisien dari model
RGWRM_parameters <- as.data.frame(RGWRM$SDF[, c(1:13)])

# Ekstrak Local R²
RGWRM_local_r2 <- RGWRM$SDF$Local_R2

# Tampilkan ringkasan
summary(RGWRM_parameters)
summary(RGWRM_local_r2)
```

P Value RGWRM

```{r}
# Konversi hasil Robust GWR ke data frame
coef_rgwrm <- as.data.frame(RGWRM[["SDF"]])

# Dapatkan daftar variabel yang memiliki "_TV" (t-statistic value)
tv_vars <- grep("_TV$", colnames(coef_rgwrm), value = TRUE)

# Hitung derajat kebebasan (df)
n <- nrow(coef_rgwrm)  # Jumlah wilayah (Kabupaten)
k <- length(tv_vars)  # Jumlah variabel prediktor dalam model
df <- n - k  # df = jumlah observasi - jumlah variabel

# Buat list untuk menyimpan p-values
p_values <- list()

# Loop untuk menghitung p-value dari setiap t-statistik
for (var in tv_vars) {
  var_name <- gsub("_TV", "", var)  # Ambil nama variabel tanpa "_TV"
  t_stat <- coef_rgwrm[[var]]  # Ambil nilai t-statistik
  p_values[[var_name]] <- 2 * (1 - pt(abs(t_stat), df))  # Hitung p-value (uji dua sisi)
}

# Konversi hasil ke data frame
p_values_df_RGWRM <- as.data.frame(p_values)

# Tambahkan nama Kabupaten
p_values_df_RGWRM$Kabupaten <- coef_rgwrm$Kabupaten

# Tampilkan ringkasan p-value
summary(p_values_df_RGWRM)
```

RGWRMM : Evaluasi Model

```{r}
AIC_RGWRM <- RGWRM[["GW.diagnostic"]][["AIC"]]
print(AIC_RGWRM)
# Ambil nilai aktual dan prediksi
Y_actual <- RGWRM$SDF$y  # Nilai Y yang sebenarnya
Y_pred <- RGWRM$SDF$yhat  # Nilai Y hasil prediksi dari model

# Hitung RMSE
RMSE_RGWRM <- sqrt(mean((Y_actual - Y_pred)^2))
print(RMSE_RGWRM)
```

RGWRS : Persiapan

```{r}
library(GWmodel)
library(sf)
library(dplyr)

# Buat shapefile khusus untuk model RGWRS
shp_kabupaten_RGWRS <- shp_pulau_jawa_merged %>%
  filter(Lokasi == "Kab") %>%
  select(Kabupaten, everything())

# Perbaiki geometri jika ada masalah
shp_kabupaten_RGWRS <- st_make_valid(shp_kabupaten_RGWRS)

# Hilangkan geometri untuk perhitungan non-spasial
data_no_geom_RGWRS <- st_drop_geometry(shp_kabupaten_RGWRS)

# Konversi koordinat menjadi matrix
coords_RGWRS <- as.matrix(st_coordinates(shp_kabupaten_RGWRS))

# Variabel respons dan prediktor (X1 hingga X12)
Y_RGWRS <- data_no_geom_RGWRS$Y
X_RGWRS <- as.matrix(cbind(1, data_no_geom_RGWRS[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")]))

# Tentukan bandwidth optimal
set.seed(123)
bandwidth_RGWRS <- bw.gwr(
  formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 + X12,
  data = as(shp_kabupaten_RGWRS, "Spatial"),
  approach = "AICc",
  kernel = "gaussian",
  adaptive = TRUE
)
```

RGWRS : MODEL

```{r}
# Maksimum iterasi & toleransi konvergensi
max_iter_RGWRS <- 20
tolerance_RGWRS <- 1e-5  
n_RGWRS <- nrow(coords_RGWRS)
p_RGWRS <- ncol(X_RGWRS)

# Inisialisasi parameter
beta_hat_RGWRS <- matrix(0, n_RGWRS, p_RGWRS)
residuals_RGWRS <- rep(0, n_RGWRS)

# Fungsi Kernel Gaussian
gaussian_kernel_RGWRS <- function(d, bw) {
  exp(- (d^2) / (2 * bw^2))
}

# Pastikan koordinat diambil dalam format matrix
coords_RGWRS <- coords_RGWRS[, c("X", "Y")]

for (iter in 1:max_iter_RGWRS) {
  cat("Iterasi ke-", iter, "\n")
  
  if (iter == 1) {
    for (i in 1:n_RGWRS) {
      dists <- sqrt(rowSums((coords_RGWRS - matrix(coords_RGWRS[i, ], nrow = n_RGWRS, ncol = 2, byrow = TRUE))^2))
      spatial_weights <- gaussian_kernel_RGWRS(dists, bandwidth_RGWRS)
      W <- diag(spatial_weights)
      
      XtWX <- t(X_RGWRS) %*% W %*% X_RGWRS
      XtWY <- t(X_RGWRS) %*% W %*% Y_RGWRS
      beta_hat_RGWRS[i, ] <- solve(XtWX, XtWY)
      residuals_RGWRS[i] <- Y_RGWRS[i] - X_RGWRS[i, ] %*% beta_hat_RGWRS[i, ]
    }
  }
  
  mad_residuals <- median(abs(residuals_RGWRS - median(residuals_RGWRS)))
  residuals_std <- residuals_RGWRS / (mad_residuals / 0.675)

  robust_weights <- ifelse(
    abs(residuals_std) <= 1.547,
    (1 - (residuals_std / 1.547)^2)^2,
    0
  )

  for (i in 1:n_RGWRS) {
    dists <- sqrt(rowSums((coords_RGWRS - matrix(coords_RGWRS[i, ], nrow = n_RGWRS, ncol = 2, byrow = TRUE))^2))
    spatial_weights <- gaussian_kernel_RGWRS(dists, bandwidth_RGWRS)
    
    combined_weights <- spatial_weights * robust_weights
    W <- diag(combined_weights)
    
    XtWX <- t(X_RGWRS) %*% W %*% X_RGWRS
    XtWY <- t(X_RGWRS) %*% W %*% Y_RGWRS
    beta_hat_RGWRS[i, ] <- solve(XtWX, XtWY)
    residuals_RGWRS[i] <- Y_RGWRS[i] - X_RGWRS[i, ] %*% beta_hat_RGWRS[i, ]
  }
  
  if (iter > 1 && max(abs(beta_hat_RGWRS - beta_hat_old_RGWRS)) < tolerance_RGWRS) {
    cat("Model RGWRS konvergen pada iterasi:", iter, "\n")
    break 
  }
  
  beta_hat_old_RGWRS <- beta_hat_RGWRS
}
```

```{r}
# Konversi hasil pendugaan koefisien menjadi data frame
coef_RGWRS <- as.data.frame(beta_hat_RGWRS)

# Tambahkan nama kabupaten dari shapefile
coef_RGWRS$Kabupaten <- shp_kabupaten_RGWRS$Kabupaten

# Tambahkan nama kolom untuk setiap koefisien
colnames(coef_RGWRS) <- c("Intercept", "X1", "X2", "X3", "X4", "X5", 
                          "X6", "X7", "X8", "X9", "X10", "X11", "X12", "Kabupaten")

# Ringkasan statistik koefisien
summary(coef_RGWRS)
```

RGWRS : SE dan Pval

```{r}
standard_errors_RGWRS <- matrix(0, n_RGWRS, p_RGWRS)

for (i in 1:n_RGWRS) {
  dists <- sqrt(rowSums((coords_RGWRS - matrix(coords_RGWRS[i, ], nrow = n_RGWRS, ncol = 2, byrow = TRUE))^2))
  spatial_weights <- gaussian_kernel_RGWRS(dists, bandwidth_RGWRS)
  combined_weights <- spatial_weights * robust_weights
  W <- diag(combined_weights)
  
  XtWX_inv <- solve(t(X_RGWRS) %*% W %*% X_RGWRS)  

  for (j in 1:p_RGWRS) {
    standard_errors_RGWRS[i, j] <- sqrt(XtWX_inv[j, j])  
  }
}
t_values_RGWRS <- beta_hat_RGWRS / standard_errors_RGWRS  
pvalue_RGWRS <- 2 * pt(-abs(t_values_RGWRS), df = n_RGWRS - p_RGWRS)  

# Konversi ke DataFrame
pvalue_RGWRS <- as.data.frame(pvalue_RGWRS)
colnames(pvalue_RGWRS) <- colnames(X_RGWRS)

# Tambahkan Nama Kabupaten
pvalue_RGWRS$Kabupaten <- shp_kabupaten_RGWRS$Kabupaten

# Ringkasan p-value
summary(pvalue_RGWRS)
```

RGWRS: Evaluasi Model

```{r}
# Hitung AIC RGWRS secara manual
AIC_RGWRS <- n_RGWRS * log(mean(residuals_RGWRS^2)) + 2 * p_RGWRS
print(AIC_RGWRS)

# Hitung RMSE
RMSE_RGWRS <- sqrt(mean(residuals_RGWRS^2))
print(RMSE_RGWRS)
```

RGWRMM : Persiapan

```{r}
library(GWmodel)
library(sf)
library(dplyr)

# Buat shapefile khusus untuk model RGWRMM
shp_kabupaten_RGWRMM <- shp_pulau_jawa_merged %>%
  filter(Lokasi == "Kab") %>%
  select(Kabupaten, everything())

# Perbaiki geometri jika ada masalah
shp_kabupaten_RGWRMM <- st_make_valid(shp_kabupaten_RGWRMM)

# Hilangkan geometri untuk perhitungan non-spasial
data_no_geom_RGWRMM <- st_drop_geometry(shp_kabupaten_RGWRMM)

# Konversi koordinat menjadi matrix
coords_RGWRMM <- as.matrix(st_coordinates(st_centroid(shp_kabupaten_RGWRMM)))

# Variabel respons dan prediktor (X1 hingga X12)
Y_RGWRMM <- data_no_geom_RGWRMM$Y
X_RGWRMM <- as.matrix(cbind(1, data_no_geom_RGWRMM[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")]))

# Tentukan bandwidth optimal menggunakan AICc
set.seed(123)
bandwidth_RGWRMM <- bw.gwr(
  formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 + X12,
  data = as(shp_kabupaten_RGWRMM, "Spatial"),
  approach = "AICc",
  kernel = "gaussian",
  adaptive = TRUE
)

print(bandwidth_RGWRMM)
```

RGWRMM: Model

```{r}
# Maksimum iterasi & toleransi konvergensi
max_iter_RGWRMM <- 20
tolerance_RGWRMM <- 1e-5  
n_RGWRMM <- nrow(coords_RGWRMM)
p_RGWRMM <- ncol(X_RGWRMM)

# Inisialisasi parameter
beta_hat_RGWRMM <- matrix(0, n_RGWRMM, p_RGWRMM)
residuals_RGWRMM <- rep(0, n_RGWRMM)

# Fungsi Kernel Gaussian
gaussian_kernel_RGWRMM <- function(d, bw) {
  exp(- (d^2) / (2 * bw^2))
}

# Iterasi RGWRMM
for (iter in 1:max_iter_RGWRMM) {
  cat("Iterasi ke-", iter, "\n")
  
  if (iter == 1) {
    for (i in 1:n_RGWRMM) {
      dists <- sqrt(rowSums((coords_RGWRMM - coords_RGWRMM[rep(i, n_RGWRMM), , drop = FALSE])^2))
      spatial_weights <- gaussian_kernel_RGWRMM(dists, bandwidth_RGWRMM)
      W <- diag(spatial_weights)
      
      XtWX <- t(X_RGWRMM) %*% W %*% X_RGWRMM
      XtWY <- t(X_RGWRMM) %*% W %*% Y_RGWRMM
      beta_hat_RGWRMM[i, ] <- solve(XtWX, XtWY)
      residuals_RGWRMM[i] <- Y_RGWRMM[i] - X_RGWRMM[i, ] %*% beta_hat_RGWRMM[i, ]
    }
  }
  
  # Hitung simpangan baku robust (MAD)
  mad_residuals <- median(abs(residuals_RGWRMM - median(residuals_RGWRMM)))
  residuals_std <- residuals_RGWRMM / (mad_residuals / 0.675)

  # Bobot robust berdasarkan M-Estimator
  robust_weights <- ifelse(
    abs(residuals_std) <= 4.685,
    (1 - (residuals_std / 4.685)^2)^2,
    0
  )

  for (i in 1:n_RGWRMM) {
    dists <- sqrt(rowSums((coords_RGWRMM - coords_RGWRMM[rep(i, n_RGWRMM), , drop = FALSE])^2))
    spatial_weights <- gaussian_kernel_RGWRMM(dists, bandwidth_RGWRMM)
    
    combined_weights <- spatial_weights * robust_weights
    W <- diag(combined_weights)
    
    XtWX <- t(X_RGWRMM) %*% W %*% X_RGWRMM
    XtWY <- t(X_RGWRMM) %*% W %*% Y_RGWRMM
    beta_hat_RGWRMM[i, ] <- solve(XtWX, XtWY)
    residuals_RGWRMM[i] <- Y_RGWRMM[i] - X_RGWRMM[i, ] %*% beta_hat_RGWRMM[i, ]
  }
  
  if (iter > 1 && max(abs(beta_hat_RGWRMM - beta_hat_old_RGWRMM)) < tolerance_RGWRMM) {
    cat("Model RGWRMM konvergen pada iterasi:", iter, "\n")
    break 
  }
  
  beta_hat_old_RGWRMM <- beta_hat_RGWRMM
}
```

RGWRMM: Koefisien

```{r}
# Konversi hasil pendugaan koefisien menjadi data frame
coef_RGWRMM <- as.data.frame(beta_hat_RGWRMM)

# Tambahkan nama kabupaten dari shapefile
coef_RGWRMM$Kabupaten <- shp_kabupaten_RGWRMM$Kabupaten

# Tambahkan nama kolom untuk setiap koefisien
colnames(coef_RGWRMM) <- c("Intercept", "X1", "X2", "X3", "X4", "X5", 
                           "X6", "X7", "X8", "X9", "X10", "X11", "X12", "Kabupaten")

# Tampilkan ringkasan parameter
summary(coef_RGWRMM)
```

RGWRMM : P value

```{r}
standard_errors_RGWRMM <- matrix(0, n_RGWRMM, p_RGWRMM)

for (i in 1:n_RGWRMM) {
  dists <- sqrt(rowSums((coords_RGWRMM - coords_RGWRMM[rep(i, n_RGWRMM), , drop = FALSE])^2))
  spatial_weights <- gaussian_kernel_RGWRMM(dists, bandwidth_RGWRMM)
  combined_weights <- spatial_weights * robust_weights
  W <- diag(combined_weights)
  
  XtWX_inv <- solve(t(X_RGWRMM) %*% W %*% X_RGWRMM)  

  for (j in 1:p_RGWRMM) {
    standard_errors_RGWRMM[i, j] <- sqrt(XtWX_inv[j, j])  
  }
}

t_values_RGWRMM <- beta_hat_RGWRMM / standard_errors_RGWRMM  
pvalue_RGWRMM <- 2 * pt(-abs(t_values_RGWRMM), df = n_RGWRMM - p_RGWRMM)  

# Konversi ke DataFrame
pvalue_RGWRMM <- as.data.frame(pvalue_RGWRMM)
colnames(pvalue_RGWRMM) <- colnames(X_RGWRMM)

# Tambahkan Nama Kabupaten
pvalue_RGWRMM$Kabupaten <- shp_kabupaten_RGWRMM$Kabupaten

# Ringkasan p-value
summary(pvalue_RGWRMM)
```

RGWRMM : Evaluasi Model

```{r}
# Hitung AIC RGWRMM
AIC_RGWRMM <- n_RGWRMM * log(mean(residuals_RGWRMM^2)) + 2 * p_RGWRMM
print(AIC_RGWRMM)

# Hitung RMSE
RMSE_RGWRMM <- sqrt(mean(residuals_RGWRMM^2))
print(RMSE_RGWRMM)
```

Pemataan Pendugaan Koefisien

```{r}
library(sf)
library(ggplot2)
library(dplyr)
library(gridExtra)

# Konversi hasil RGWRMM ke sf dengan geometri dari shp_kabupaten_RGWRMM
shp_kabupaten_RGWRMM_mapped <- st_as_sf(
  cbind(shp_kabupaten_RGWRMM, coef_RGWRMM),
  geometry = st_geometry(shp_kabupaten_RGWRMM)
)

# Daftar variabel koefisien yang akan dipetakan (tanpa Intercept)
vars_to_plot <- c("X1.1", "X2.1", "X3.1", "X4.1", 
                  "X5.1", "X6.1", "X7.1", "X8.1", 
                  "X9.1", "X10.1", "X11.1", "X12.1")

# Nama variabel dalam visualisasi
var_labels <- c(
  "X1.1"  = "Complete Immunization",
  "X2.1"  = "Growth Monitoring",
  "X3.1"  = "Maternal Education",
  "X4.1"  = "ARI Prevalence",
  "X5.1"  = "Birth Weight",
  "X6.1"  = "Diarrhoea Incidence",
  "X7.1"  = "Exclusive Breastfeeding",
  "X8.1"  = "Sanitation Access",
  "X9.1"  = "Health Workers Ratio",
  "X10.1" = "Food Security Index",
  "X11.1" = "Poverty Rate",
  "X12.1" = "Child Mortality"
)

# **Pemisahan ke dalam 3 bagian (setiap frame berisi 4 variabel)**
vars_part1 <- vars_to_plot[1:4]   # Frame 1
vars_part2 <- vars_to_plot[5:8]   # Frame 2
vars_part3 <- vars_to_plot[9:12]  # Frame 3

# **Fungsi untuk membuat plot berdasarkan kelompok variabel**
plot_rgwrmm_maps <- function(var_subset) {
  plots <- list()
  for (var in var_subset) {
    p <- ggplot(data = shp_kabupaten_RGWRMM_mapped) +
      geom_sf(aes_string(fill = var), color = "black", size = 0.2) +  
      scale_fill_gradient2(
        low = "red", mid = "white", high = "blue", midpoint = 0, 
        name = "Coefficient"
      ) +
      labs(
        title = var_labels[var]
      ) +
      theme_minimal() +
      theme(
        plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
        legend.position = "right",
        legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(0.4, "cm"),
        plot.margin = margin(5, 5, 5, 5)
      ) +
      coord_sf(expand = TRUE)  
  
    plots[[var]] <- p
  }
  do.call(grid.arrange, c(plots, ncol = 2))  # Menampilkan 4 variabel dalam satu frame (2x2)
}

# **Menampilkan masing-masing bagian secara terpisah**
cat("\n===== PETA PENDUGAAN KOEFISIEN RGWRMM (Frame 1: X1 - X4) =====\n")
plot_rgwrmm_maps(vars_part1)

cat("\n===== PETA PENDUGAAN KOEFISIEN RGWRMM (Frame 2: X5 - X8) =====\n")
plot_rgwrmm_maps(vars_part2)

cat("\n===== PETA PENDUGAAN KOEFISIEN RGWRMM (Frame 3: X9 - X12) =====\n")
plot_rgwrmm_maps(vars_part3)
```

Pengelompokan Peubah

```{r}
# Bagi tiap variabel ke dalam 3 kuartil: Q1, Q2, Q3 (tertile)
for (var in vars_to_plot) {
  q <- quantile(shp_kabupaten_RGWRMM_mapped[[var]], probs = c(1/3, 2/3), na.rm = TRUE)
  shp_kabupaten_RGWRMM_mapped[[paste0(var, "_q3")]] <- cut(
    shp_kabupaten_RGWRMM_mapped[[var]],
    breaks = c(-Inf, q[1], q[2], Inf),
    labels = c("Q1", "Q2", "Q3"),
    include.lowest = TRUE
  )
}

# Fungsi plotting baru berdasarkan Q1-Q3
plot_rgwrmm_q3_maps <- function(var_subset) {
  plots <- list()
  for (var in var_subset) {
    var_q3 <- paste0(var, "_q3")
    
    p <- ggplot(data = shp_kabupaten_RGWRMM_mapped) +
      geom_sf(aes_string(fill = var_q3), color = "black", size = 0.2) +
      scale_fill_manual(
        values = c("Q1" = "#d73027",  # merah
                   "Q2" = "#fc8d59",  # oranye
                   "Q3" = "#1a9850"), # hijau
        name = "Tertile"
      ) +
      labs(title = var_labels[var]) +
      theme_minimal() +
      theme(
        plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
        legend.position = "right",
        legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(0.4, "cm"),
        plot.margin = margin(5, 5, 5, 5)
      ) +
      coord_sf(expand = TRUE)
    
    plots[[var]] <- p
  }
  do.call(grid.arrange, c(plots, ncol = 2))
}

# Tampilkan hasil hanya Q1-Q3
cat("\n===== PETA KOEFISIEN RGWRMM BERDASARKAN Q1-Q3 (Frame 1: X1 - X4) =====\n")
plot_rgwrmm_q3_maps(vars_part1)

cat("\n===== PETA KOEFISIEN RGWRMM BERDASARKAN Q1-Q3 (Frame 2: X5 - X8) =====\n")
plot_rgwrmm_q3_maps(vars_part2)

cat("\n===== PETA KOEFISIEN RGWRMM BERDASARKAN Q1-Q3 (Frame 3: X9 - X12) =====\n")
plot_rgwrmm_q3_maps(vars_part3)
```
```{r}
# Tambahkan kolom kombinasi tanda dan tertile (Q1-Q3)
for (var in vars_to_plot) {
  value <- shp_kabupaten_RGWRMM_mapped[[var]]
  
  # Hitung batas tertile (1/3 dan 2/3 quantile)
  tertiles <- quantile(value, probs = c(1/3, 2/3), na.rm = TRUE)
  tertile_class <- cut(
    value,
    breaks = c(-Inf, tertiles[1], tertiles[2], Inf),
    labels = c("Q1", "Q2", "Q3"),
    include.lowest = TRUE
  )
  
  # Tentukan tanda: Positif / Negatif
  sign_class <- ifelse(value < 0, "Neg", "Pos")
  
  # Gabungkan menjadi satu label
  combined_class <- paste(sign_class, tertile_class, sep = "_")
  
  # Simpan dalam kolom baru
  shp_kabupaten_RGWRMM_mapped[[paste0(var, "_sign_q3")]] <- combined_class
}

plot_rgwrmm_sign_q3_maps <- function(var_subset) {
  plots <- list()
  
  for (var in var_subset) {
    var_sign_q3 <- paste0(var, "_sign_q3")
    
    p <- ggplot(data = shp_kabupaten_RGWRMM_mapped) +
      geom_sf(aes_string(fill = var_sign_q3), color = "black", size = 0.2) +
      scale_fill_manual(
        values = c(
          "Neg_Q1" = "#67001f",  # merah tua (negatif kecil)
          "Neg_Q2" = "#b2182b",  # merah (negatif sedang)
          "Neg_Q3" = "#ef8a62",  # merah muda (negatif besar)
          "Pos_Q1" = "#92c5de",  # biru muda (positif kecil)
          "Pos_Q2" = "#4393c3",  # biru (positif sedang)
          "Pos_Q3" = "#2166ac"   # biru tua (positif besar)
        ),
        name = "Tanda & Tertile"
      ) +
      labs(title = var_labels[var]) +
      theme_minimal() +
      theme(
        plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
        legend.position = "right",
        legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(0.4, "cm"),
        plot.margin = margin(5, 5, 5, 5)
      ) +
      coord_sf(expand = TRUE)
    
    plots[[var]] <- p
  }
  
  do.call(grid.arrange, c(plots, ncol = 2))
}

cat("\n===== PETA RGWRMM TERTILE + TANDA (Frame 1: X1 - X4) =====\n")
plot_rgwrmm_sign_q3_maps(vars_part1)

cat("\n===== PETA RGWRMM TERTILE + TANDA (Frame 2: X5 - X8) =====\n")
plot_rgwrmm_sign_q3_maps(vars_part2)

cat("\n===== PETA RGWRMM TERTILE + TANDA (Frame 3: X9 - X12) =====\n")
plot_rgwrmm_sign_q3_maps(vars_part3)
```