Fenomena Kasus Diabetes Mellitus di Provinsi Jawa Barat Periode 2022-2024

Dosen Pembimbing:

Dr. I Gede Nyoman Mindra Jaya, M.Si.

Disusun Oleh:

Rasendriya Nandana Kurniawan 140610230009

PROGRAM STUDI S-1 STATISTIKA

DEPARTEMEN STATISTIKA

FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM

UNIVERSITAS PADJADJARAN

TAHUN AJARAN 2025/2026

Bab 1. Pendahuluan

1.1. Latar Belakang

Diabetes Mellitus (DM) merupakan salah satu penyakit tidak menular yang telah menjadi salah satu ancaman kesehatan global terbesar sebagai penyumbang angka morbiditas dan mortalitas yang signifikan setiap tahunnya. International Diabetes Federation (2021) melaporkan bahwa prevalensi diabetes terus mengalami peningkatan drastis, terutama di negara-negara berkembang, terutama Indonesia yang menduduki peringkat kelima dengan jumlah kasus terbanyak per tahun 2021. Salah satu provinsi dengan kasus Diabetes Mellitus terbanyak di Indonesia adalah Provinsi Jawa Barat. Per tahun 2024, tercatat ada total 668.220 kasus Diabetes Mellitus yang tercatat dari 27 Kabupaten dan Kota di Jawa Barat. Jumlah kasus tersebut lebih banyak dari tahun lalu, dimana pada tahun 2023 terdapat sebanyak 645.390 kasus yang tercatat.

Lonjakan kasus serta prevalensi kasus Diabetes Mellitus secara global di era sekarang sangat berkorelasi dengan perubahan gaya hidup yang dipicu oleh fenomena urbanisasi dan modernisasi secara cepat. Faktor-faktor seperti penurunan aktivitas fisik karena dominasi pekerjaan yang minim gerakan fisik, serta maraknya konsumsi makanan cepat saji yang rentan tinggi kalori dan rendah serat menjadi katalisator utama resistensi insulin. Fenomena ini diperburuk oleh pergeseran demografis penuaan penduduk (aging population), yang secara alami meningkatkan risiko degeneratif, namun kini diimbangi oleh tren mengkhawatirkan munculnya early-onset diabetes pada populasi usia produktif dan remaja.

Sebagai provinsi dengan populasi terbesar di Indonesia, Jawa Barat memiliki kerentanan yang tinggi terhadap peningkatan kasus penyakit tidak menular, salah satunya adalah Diabetes Mellitus. Kepadatan penduduk yang tinggi serta kesenjangan sosial-ekonomi antarwilayah menciptakan pola penyebaran penyakit yang tidak seragam. Faktor-faktor risiko seperti tingkat kemiskinan, rata-rata lama sekolah, rasio penduduk lansia diyakini, dan kepadatan penduduk memiliki kontribusi terhadap prevalensi kasus Diabetes Mellitus. Analisis dan Pemodelan Spasial Diabetes Mellitus di masa ini menjadi hal yang sangat penting bukan hanya untuk mengetahui besaran masalah, tetapi juga memahami distribusi dan korelasi antarwilayah bertetanggaan dalam konteks jumlah kasus Diabetes Mellitus

1.2. Identifikasi Masalah

Adapun rumusan masalah yang didapatkan dalam penyusunan penelitian ini adalah:

1. Apakah ada ketergantungan spasial pada kasus persebaran Diabetes Mellitus per kabupaten dan kota di Jawa Barat tahun 2022 sampai 2024?
2. Apakah persentase lansia, kepadatan penduduk, persentase penduduk miskin, dan rata-rata lama sekolah menjadi faktor yang mempengaruhi kasus Diabetes Mellitus secara signifikan?
3. Model spasial apakah yang paling cocok dan kuat untuk memodelkan distribusi kasus Diabetes Mellitus di Jawa Barat?

1.3. Tujuan

Adapun tujuan penyusunan penelitian ini adalah:

1. Mengetahui distribusi kasus Diabetes Mellitus di Jawa Barat lewat peta tematik
2. Mengidentifikasi ada tidaknya autokorelasi spasial secara global dan lokal pada kasus persebaran Diabetes Mellitus per kabupaten dan kota di Jawa Barat tahun 2022 sampai 2024
3. Mengidentifikasi apakah persentase lansia, kepadatan penduduk, persentase penduduk miskin, dan rata-rata lama sekolah menjadi faktor yang mempengaruhi kasus Diabetes Mellitus secara signifikan
4. Mengidentifikasikan model spasial apakah yang paling cocok dan kuat untuk memodelkan distribusi kasus Diabetes Mellitus di Jawa Barat

1.4. Batasan Penelitian

Adapun batasan penelitian ini adalah:

1. Penelitian ini berfokus pada kasus Diabetes Mellitus di Provinsi Jawa Barat per wilayah kabupaten dan kota pada tahun 2022 sampai 2024
2. Variabel dependen yang digunakan adalah kasus Diabetes Mellitus
3. Variabel independen yang digunakan dibatasi pada 4 faktor, yaitu Persentase Penduduk Miskin, Rata-Rata Lama Sekolah, Persentase Penduduk Lansia (Lanjut Usia), dan Kepadatan Penduduk

Bab 2. Tinjauan Pustaka

2.1. Teori Dasar Spatial Dependence

Spatial dependence atau ketergantungan spasial menggambarkan kondisi di mana suatu fenomena pada satu lokasi dipengaruhi oleh fenomena di lokasi sekitarnya. Dalam konteks kesehatan masyarakat, seperti kasus Diabetes Mellitus, wilayah dengan tingkat kejadian tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki tingkat kejadian tinggi. Prinsip ini sejalan dengan Tobler’s First Law of Geography yang menyatakan bahwa “everything is related to everything else, but near things are more related than distant things.” Oleh karena itu, analisis spasial diperlukan untuk mengidentifikasi adanya pola keterkaitan antarwilayah, yang tidak dapat dijelaskan hanya dengan analisis statistik konvensional.

2.2. Autokorelasi Spasial

Autokorelasi spasial digunakan untuk mengukur sejauh mana nilai suatu variabel di satu lokasi berkorelasi dengan nilai di lokasi lain yang berdekatan. Terdapat dua jenis analisis autokorelasi spasial, yaitu Global Moran’s I dan Local Moran’s I (LISA). Global Moran’s I digunakan untuk menilai apakah pola spasial secara keseluruhan bersifat mengelompok (clustered), acak (random), atau menyebar (dispersed), sedangkan Local Moran’s I digunakan untuk mengidentifikasi lokasi-lokasi spesifik yang membentuk klaster signifikan.

2.2.1 Global Moran’s I

Global Moran’s I merupakan sebuah statistik inferensial dengan tujuan untuk mendeteksi ada tidaknya autokorelasi spasial secara signifikan pada seluruh lokasi yang dijadikan bahan studi. Nilai Moran’s I memiliki rentang antara -1 hingga +1, dengan autokorelasi spasial positif apabila nilai statistiknya mendekati 1 dan autokorelasi spasial negatif apabila nilai statistiknya mendekati -1.

Adapun rumus untuk menghitung nilai Global Moran’s I adalah sebagai berikut:

\[ I = \frac{N}{S_0} \cdot \frac{\sum_{i=1}^{N} \sum_{j=1}^{N} w_{ij} (x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^{N} (x_i - \bar{x})^2} \]

2.2.2 Local Moran’s I

Local Moran’s I merupakan salah satu bagian dari Local Indicators of Spatial Association (LISA) yang berfungsi untuk mendeteksi adanya autokorelasi spasial secara lokal yang signifikan. Berbeda dengan Global Moran’s I yang mendeteksi secara global, Local Moran’s I mendeteksi autokorelasi spasial secara lokal, yaitu per wilayah uji untuk melihat autokorelasi spasial lebih dalam. Hasil pengujian Local Moran’s I umum ditampilkan dalam peta signifikansi LISA yang akan dibagi menjadi 4 kuadran utama: HH (High-High), LL (Low-Low), HL (High-Low), dan LH (Low-High)

Adapun rumus untuk menghitung nilai statistik Local Moran’s I adalah:

\[ I_i = \frac{x_i - \bar{x}}{m_2} \sum_{j=1}^{N} w_{ij}(x_j - \bar{x}), \quad \text{dengan: } m_2 = \frac{1}{N} \sum_{k=1}^{N} (x_k - \bar{x})^2 \]

2.3. Spatial Econometrics Model

Model ekonometrika spasial merupakan model yang digunakan untuk menganalisis ketergantungan dan heterogenitas spasial dalam data geografi. Tujuan utamanya adalah untuk mengatasi masalah-masalah interpretasi apabila asumsi model regresi standar (Ordinary Least Square) tidak cocok digunakan ketika asumsi non-autokorelasi spasial tidak terpenuhi.

2.3.1 Ordinary Least Square

Metode Kuadrat terkecil, atau biasa dikenal sebagai Ordinary Least Square (OLS) merupakan metode estimasi yang bertujuan untuk meminimalkan nilai SSR (Sum of Squared Residual) antara nilai variabel terikat yang diamati dan nilai yang diprediksi oleh model. Model OLS mengasumsikan bahwa objek pengamatan yang diuji bersifat independen dan iid (independent identically distributed)

2.3.2 Spatial Autoregressive

Spatial Autoregressive (SAR), atau juga dikenal sebagai Spatial Lag Model (SLM) merupakan model yang ditujukan untuk menangani fenomena spatial dependencies dalam variabel yang digunakan. Model ini menggunakkan variabel terikat ayng tertimbang secara spasial sebagai tambahan prediktor untuk menangkap efek spillover antara unit-unit spasial.

Adapun rumus untuk model Spatial Autoregressive adalah: \[ y = \rho W y + X \beta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I). \]

2.3.3 Spatial Error Model

Spatial Error Model (SEM) merupakan semua model dalam spasial yang ditujukan untuk menangani masalah autokorelasi spasial pada residual. Model ini berasumsi bahwa fenomena spatial dependencies muncul dari variabel yang relevan dan terstruktur secara spasial, namun dihilangkan dari model regresi

Adapun rumus untuk model Spatial Error Model adalah:

\[ y = X \beta + u, \quad \text{dengan:} \]

\[ \begin{aligned} u &= \lambda W u + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I). \end{aligned} \]

2.4. Interpolasi Spasial

Interpolasi spasial merupakan sebuah prosedur untuk memprediksi nilai pada lokasi yang tidak memiliki data dengan nilai-nilai dari lokasi yang ada disekitarnya dengan data yang sudah diketahui dengan asumsi lokasi yang berdekatan memiliki korelasi yang kuat. Tujuannya adalah untuk mentransformasikan wilayah dari yang memiliki batas administrasi menjadi sebuah permukaan kontinu yang halus, sehingga memungkinkan visualisasi gradasi sebaran kasus yang lebih realistis

Bab 3. Metodologi Penelitian

3.1. Data

Data yang digunakan pada penelitian ini merupakan data yang dapat diakses secara open source lewat website. Adapun variabel dependen dan independen beserta sumber yang digunakan pada penelitian ini sebagai berikut:

Nama Variabel	Jenis Variabel	Sumber Data
Jumlah Penderita Diabetes Melitus	Dependen	Open Data Jabar
Persentase Penduduk Miskin	Independen	Open Data Jabar
Rata-Rata Lama Sekolah	Independen	Open Data Jabar
Kepadatan Penduduk	Independen	Open Data Jabar
Persentase Penduduk Lansia (Lanjut Usia)	Independen	Badan Pusat Statistik Jawa Barat

3.2. Unit Spasial

Unit Spasial yang digunakan pada penelitian ini adalah 27 kabupaten dan kota yang ada dalam Provinsi Jawa Barat. Adapun nama-nama kabupaten dan kota yang menjadi unit spasial pada penelitian ini adalah Bandung, Bandung Barat, Bekasi, Bogor, Ciamis, Cianjur, Cirebon, Garut, Indramayu, Karawang, Kuningan, Majalengka, Pangandaran, Purwakarta, Subang, Sukabumi, Sumedang, Tasikmalaya, Kota Bandung, Kota Banjar, Kota Bekasi, Kota Bogor, Kota Cirebon, Kota Depok, Kota Sukabumi, Kota Tasikmalaya, dan Kota Cimahi. Setiap unit spasial direpresentasikan dalam bentuk poligon pada peta digital (shapefile), yang berfungsi untuk menganalisis hubungan spasial antar wilayah.

3.3. Metode Analisis

Metode analisis yang diterapkan dalam penelitian ini analisis autokorelasi spasial dan pemodelan Regresi Data Panel. Analisis diawali dengan statistika deskriptif spasial melalui peta choropleth serta penerapan metode interpolasi spasial untuk memvisualisasikan sebaran kasus Diabetes Mellitus dan mengestimasi kasus pada seluruh permukaan wilayah tanpa terpaku pada batas wilayah. Lalu, dilakukan deteksi autokorelasi spasial menggunakan uji Global Moran’s I dan Local Indicator of Spatial Association (LISA), yang menjadi dasar validasi untuk melakukan pemodelan Regresi Data Panel Spasial. Pada tahap pemodelan, spesifikasi model terbaik, antara Spatial Autoregressive Model (SAR) atau Spatial Error Model (SEM)

3.4. Alur Kerja Penelitian

Alur kerja penelitian dimulai dengan tahap pengumpulan dan pre-processing data yang wilayah antara data atribut statistik tahun 2022–2024 dengan data geometri peta (shapefile) Jawa Barat. Setelah data terintegrasi menjadi satu, dilakukan eksplorasi data panel yang tidak hanya melihat tren temporal, tetapi juga melibatkan pemetaan interpolasi untuk mengidentifikasi pola gradasi risiko penyebaran penyakit yang lebih halus antarwilayah. Hasil visualisasi tersebut kemudian diuji secara statistik melalui uji asumsi autokorelasi spasial, yang dilanjutkan dengan estimasi model regresi spasial sesuai spesifikasi terpilih (SAR/SEM) serta evaluasi model, sebelum akhirnya dilakukan interpretasi komprehensif terhadap faktor-faktor determinan penyebaran penyakit.

Bab 4. Hasil dan Pembahasan

Dashboard untuk hasil dan pembahasan dapat dilihat lewat tautan berikut: https://nanaurum.shinyapps.io/Spasial_Diabetes_Mellitus/

Persiapan untuk analisis spasial

# Packages
library(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)
library(spData)
library(geodata)
library(tmap)
library(skimr)
library(psych)
library(gridExtra)
library(grid)
library(viridis)
library(spatialreg)
library(stringr)
library(car)
library(plm)
library(tidyr)
library(splm)

#### ---------------------------------------------------------------------------
# Memanggil File .shp Untuk Peta Chloropeth (Geometri)
#### ---------------------------------------------------------------------------

# Baca dan persiapan dari file .shp
shp <- st_read("geoBoundaries-IDN-ADM2-all/geoBoundaries-IDN-ADM2.shp")
shp <- shp %>%
  mutate(
    NAME_join = shapeName 
  )

# Daftar kab/kota di Provinsi Jawa Barat
kabkota_jabar <- c(
  "Bandung", "Bandung Barat", "Bekasi", "Bogor", "Ciamis",
  "Cianjur", "Cirebon", "Garut", "Indramayu", "Karawang",
  "Kota Bandung", "Kota Banjar", "Kota Bekasi", "Kota Bogor",
  "Kota Cimahi", "Kota Cirebon", "Kota Depok", "Kota Sukabumi",
  "Kota Tasikmalaya", "Kuningan", "Majalengka", "Pangandaran",
  "Purwakarta", "Subang", "Sukabumi", "Sumedang", "Tasikmalaya"
)

# Filter Shapefile hanya untuk Kabupaten/Kota Jawa Barat
jabar_shp_filtered <- shp %>%
  filter(NAME_join %in% kabkota_jabar)

# Memanggil File .csv
data_csv <- read.csv("UAS - Spasial Statistics - Data Panel.csv", stringsAsFactors = FALSE)
col_join_csv_name <- "Kab_Kota" 

# Bersihkan nama kota di CSV agar cocok dengan Peta
data_csv <- data_csv %>%
  mutate(
    temp_name = str_to_title(str_trim(.data[[col_join_csv_name]])), 
    NAME_join = str_replace(temp_name, "Kabupaten ", "") 
  )

# Gabungkan Shapefile (Geometri) dan Data CSV (Atribut)
jabar_data <- left_join(jabar_shp_filtered, data_csv, by = "NAME_join") %>%
  rename(KabKota = NAME_join) 

# Variabel Y
vars_Y_panel <- c("DM_2022", "DM_2023", "DM_2024")

Persiapan Spasial Weight

nb_full <- poly2nb(jabar_data, queen = TRUE)
lw_full <- nb2listw(nb_full, style = "W", zero.policy = TRUE)

4.1. Peta Deskriptif

# Mode Plot
tmap_mode("plot")

# List Kosong
peta_list <- list()

# Buat Peta
for (var in vars_Y_panel) {
  tahun <- substr(var, nchar(var)-3, nchar(var))
  
  # Buat peta individual
  p <- tm_shape(jabar_data) +
    tm_fill(col = var, 
            style = "quantile", 
            palette = "Reds", 
            title = "Jumlah Kasus") + 
    tm_borders(col = "gray40", lwd = 0.5) +
    
    tm_layout(
      # Judul
      main.title = paste0("Persebaran Kasus Diabetes Mellitus\nJawa Barat Tahun ", tahun),
      main.title.position = "center",
      main.title.size = 0.9,
      main.title.fontface = "bold",
      
      # Legenda
      legend.outside = TRUE,             
      legend.outside.position = "right", 
      frame = FALSE, 
      legend.title.size = 0.8,
      legend.text.size = 0.6
    )
  
  # Simpan ke list
  peta_list[[var]] <- p
}


# Buat halaman baru
grid.newpage()

# Ukuran Grid
pushViewport(viewport(layout = grid.layout(2, 4)))

# Cetak Peta pada posisinya masing-masing
print(peta_list[[1]], vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
print(peta_list[[2]], vp = viewport(layout.pos.row = 1, layout.pos.col = 3:4))
print(peta_list[[3]], vp = viewport(layout.pos.row = 2, layout.pos.col = 2:3))

4.2. Uji Autokorelasi Spasial

4.2.1 Autokorelasi Spasial Global

# Hasil Rekap Uji
hasil_rekap <- data.frame(
  Tahun = character(),
  Moran_Index = numeric(),
  P_Value = numeric(),
  Kesimpulan = character(),
  stringsAsFactors = FALSE
)

# Uji Autokorelasi Spasial Global
for (var in vars_Y_panel) {
  
  uji <- moran.test(jabar_data[[var]], lw_full, zero.policy = TRUE)
  nilai_moran <- round(uji$estimate[1], 4)
  nilai_pval  <- round(uji$p.value, 4)
  
  # Tentukan Kesimpulan (Signifikan atau Tidak)
  kesimpulan <- ifelse(nilai_pval < 0.05, "Signifikan (Ada Pola)", "Tidak Signifikan (Acak)")
  
  # Masukkan ke Baris Baru
  baris_baru <- data.frame(
    Tahun = var,
    Moran_Index = nilai_moran,
    P_Value = nilai_pval,
    Kesimpulan = kesimpulan
  )
  
  hasil_rekap <- rbind(hasil_rekap, baris_baru)
}

print(hasil_rekap)

##                      Tahun Moran_Index P_Value            Kesimpulan
## Moran I statistic  DM_2022      0.2835  0.0099 Signifikan (Ada Pola)
## Moran I statistic1 DM_2023      0.2274  0.0272 Signifikan (Ada Pola)
## Moran I statistic2 DM_2024      0.2411  0.0217 Signifikan (Ada Pola)

Dengan taraf signifikansi sebesar 5% (0,05), terlihat bahwa terdapat indikasi autokorelasi spasial positif yang kuat pada ketiga tahun tersebut, dengan nilai p-value yang lebih kecil daripada taraf signifikansi dan nilai Moran’s I yang lebih besar dari 0.

4.2.2. Autokorelasi Spasial Lokal

# 1. Variabel 3 Tahun
plot_list <- list()

# 2. Looping
for (var in vars_Y_panel) {
  
  # --- A. PERHITUNGAN (Tetap sama) ---
  local_moran <- localmoran(jabar_data[[var]], lw_full, zero.policy = TRUE)
  
  vec_y    <- jabar_data[[var]]
  mean_y   <- mean(vec_y, na.rm = TRUE)
  lag_y    <- lag.listw(lw_full, vec_y, zero.policy = TRUE)
  mean_lag <- mean(lag_y, na.rm = TRUE)
  p_val    <- local_moran[, 5]
  
  quadrant <- rep("Not Sig.", length(vec_y))
  
  for (i in seq_along(vec_y)) {
    if (p_val[i] <= 0.05) { 
      if (vec_y[i] >= mean_y && lag_y[i] >= mean_lag) quadrant[i] <- "High-High"
      else if (vec_y[i] <= mean_y && lag_y[i] <= mean_lag) quadrant[i] <- "Low-Low"
      else if (vec_y[i] >= mean_y && lag_y[i] <= mean_lag) quadrant[i] <- "High-Low"
      else if (vec_y[i] <= mean_y && lag_y[i] >= mean_lag) quadrant[i] <- "Low-High"
    }
  }
  
  col_name_lisa <- paste0("LISA_", var)
  jabar_data[[col_name_lisa]] <- factor(quadrant, 
                                        levels = c("High-High", "Low-Low", "High-Low", "Low-High", "Not Sig."))
  
  # --- B. FILTER DATA (Hanya labeli yang Signifikan) ---
  data_signifikan <- jabar_data[jabar_data[[col_name_lisa]] != "Not Sig.", ]
  
  # --- C. VISUALISASI ---
  tahun <- substr(var, nchar(var)-3, nchar(var))
  
  p <- tm_shape(jabar_data) +
    tm_fill(col = col_name_lisa, 
            palette = c("red", "blue", "lightpink", "skyblue2", "white"),
            title = "Kategori LISA",
            drop.levels = FALSE) +
    tm_borders(col = "gray40", lwd = 0.5) +
    
    # --- PERBAIKAN DI SINI: Ganti "NAME_join" jadi "Kab_Kota" ---
    tm_shape(data_signifikan) +
    tm_text("Kab_Kota",             # Pastikan ini sesuai nama kolom di CSV Anda
            size = 0.5,             
            col = "black",          
            shadow = TRUE,          
            auto.placement = TRUE) +
    
    tm_layout(
      main.title = paste0("Peta Cluster LISA\nTahun ", tahun),
      main.title.position = "center",
      main.title.size = 0.9,
      main.title.fontface = "bold",
      legend.outside = TRUE,
      legend.outside.position = "right",
      frame = FALSE
    )
  
  plot_list[[var]] <- p
}

# 3. TAMPILKAN GRID
library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 4)))
print(plot_list[[1]], vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
print(plot_list[[2]], vp = viewport(layout.pos.row = 1, layout.pos.col = 3:4))
print(plot_list[[3]], vp = viewport(layout.pos.row = 2, layout.pos.col = 2:3))

Berdasarkan hasil analisis Local Moran’s I (LISA) pada peta kluster spasial, teridentifikasi beberapa pola keruangan kasus Diabetes Mellitus di Jawa Barat. Adapun hasil autokorelasi spasial lokal yang didapatkan adalah sebagai berikut:

• Hampir semua kota dan kabupaten yang terindikasi autokorelasi spasial merupakan wilayah yang sama pada 3 tahun tersebut, yaitu Kabupaten Bekasi, Kota Bekasi, Kota Depok, Kota Bogor, dan Kabupaten Ciamis. Hanya Kabupaten Ciamis yang tidak terindikasi autokorelasi spasial pada tahun 2023.
• Semua kota dan kabupaten yang terindikasi autokorelasi spasial memiliki tipe klaster yang konsisten pada rentang 3 tahun tersebut
• Klaster High-High teridentifikasi pada Kabupaten Bekasi, Kota Bekasi, dan Kota Depok periode 2022, 2023, dan 2024
• Klaster Low-High teridentifikasi pada Kota Bogor periode 2022, 2023, dan 2024
• Klaster Low-Low teridentifikasi pada Kabupaten Ciamis periode 2022 dan 2024

4.3. Pemodelan Regresi Spasial

4.3.1 Regresi Data Panel

Tahap 1: Penyusunan data untuk regresi data panel

df_raw <- st_drop_geometry(jabar_data) 

# Reshape ke Long Format
data_panel <- df_raw %>%
  pivot_longer(
    cols = matches("_202"),
    names_to = c(".value", "Tahun"),
    names_pattern = "(.*)_(\\d{4})"
  ) %>%
  mutate(Tahun = as.integer(Tahun)) %>%
  filter(Tahun >= 2022) # Ambil 2022-2024

# Deklarasikan sebagai Data Panel
pdata <- pdata.frame(data_panel, index = c("Kab_Kota", "Tahun"))

# --- B. ESTIMASI MODEL ---

# Rumus Regresi
rumus <- DM ~ PPM + RRLS + KpdtPddk + PerLansia

# 1. Common Effect Model (CEM)
cem <- plm(rumus, data = pdata, model = "pooling")
# 2. Fixed Effect Model (FEM)
fem <- plm(rumus, data = pdata, model = "within")
# 3. Random Effect Model (REM)
rem <- plm(rumus, data = pdata, model = "random")

Tahap 2: Pemilihan Model Data Panel

# UJI CHOW (CEM vs FEM)
uji_chow <- pFtest(fem, cem)
print(uji_chow)

## 
##  F test for individual effects
## 
## data:  rumus
## F = 19.123, df1 = 26, df2 = 50, p-value < 2.2e-16
## alternative hypothesis: significant effects

# UJI LAGRANGE MULTIPLIER (CEM vs REM)
uji_lm <- plmtest(cem, type = "bp") # Breusch-Pagan
print(uji_lm)

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  rumus
## chisq = 53.48, df = 1, p-value = 2.612e-13
## alternative hypothesis: significant effects

Dengan taraf signifikansi sebesar 5% (0.05), terlihat dari hasil kedua uji tersebut bahwa model REM dan FEM lebih baik daripada model CEM, dilihat dari hasil p-value yang yang dihasilkan dari kedua uji tersebut lebih kecil daripada taraf signifikansi, sehingga model FEM (dari uji Chow) dan model REM (dari uji LM-BP) lebih ideal dibandingkan dari model CEM, sehingga kedua model tersebut perlu dibandingkan kembali untuk menentukan model terbaik

# UJI HAUSMAN (FEM vs REM)
uji_hausman <- phtest(rem, fem)
print(uji_hausman)

## 
##  Hausman Test
## 
## data:  rumus
## chisq = 5.1814, df = 4, p-value = 0.2692
## alternative hypothesis: one model is inconsistent

Dengan taraf signifikansi sebesar 5% (0.05), terlihat dari hasil uji Hausman bahwa model REM lebih baik dibanding dengan model FEM dilihat dari hasil p-value yang yang dihasilkan lebih kecil daripada taraf signifikansi, sehingga model REM (dari uji Chow) lebih ideal untuk digunakan dibandingkan model FEM

4.3.2. Spatial Econometrics Model

Model Regresi Spasial

# Model Non-Spasial (REM)
model_rem <- plm(rumus, data = pdata, model = "random")

# Model SAR (REM)
model_sar <- spml(rumus, data = pdata, listw = lw_full, model = "random", lag = TRUE, spatial.error = "none")

# Model SEM (REM)
model_sem <- spml(rumus, data = pdata, listw = lw_full, model = "random", lag = FALSE, spatial.error = "b")

summary(model_rem);summary(model_sar);summary(model_sem)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = rumus, data = pdata, model = "random")
## 
## Balanced Panel: n = 27, T = 3, N = 81
## 
## Effects:
##                     var   std.dev share
## idiosyncratic  25175323      5018 0.127
## individual    172707817     13142 0.873
## theta: 0.7847
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -6301.2 -3293.2  -735.7  1785.3 28633.1 
## 
## Coefficients:
##                Estimate  Std. Error z-value Pr(>|z|)  
## (Intercept)  9.8338e+04  4.3692e+04  2.2507  0.02441 *
## PPM         -3.1931e+03  1.2588e+03 -2.5365  0.01120 *
## RRLS        -4.0197e+03  4.6071e+03 -0.8725  0.38294  
## KpdtPddk     8.0816e-02  1.3630e+00  0.0593  0.95272  
## PerLansia   -1.0372e+03  5.0253e+02 -2.0639  0.03903 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2172100000
## Residual Sum of Squares: 1.936e+09
## R-Squared:      0.10871
## Adj. R-Squared: 0.061803
## Chisq: 9.2699 on 4 DF, p-value: 0.054696

## ML panel with spatial lag, random effects 
## 
## Call:
## spreml(formula = formula, data = data, index = index, w = listw2mat(listw), 
##     w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl)
## 
## Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -29650  -11905   -3899   -3656    4124   45256 
## 
## Error variance parameters:
##     Estimate Std. Error t-value Pr(>|t|)   
## phi   6.6087     2.4433  2.7048 0.006834 **
## 
## Spatial autoregressive coefficient:
##        Estimate Std. Error t-value Pr(>|t|)
## lambda -0.13171    0.14999 -0.8781   0.3799
## 
## Coefficients:
##                Estimate  Std. Error t-value Pr(>|t|)   
## (Intercept) 97681.39419 41668.76633  2.3442 0.019066 * 
## PPM         -3213.70110  1200.43144 -2.6771 0.007426 **
## RRLS        -3386.80925  4393.40148 -0.7709 0.440775   
## KpdtPddk       -0.19971     1.29857 -0.1538 0.877775   
## PerLansia   -1040.15802   483.09118 -2.1531 0.031308 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## ML panel with , random effects, spatial error correlation 
## 
## Call:
## spreml(formula = formula, data = data, index = index, w = listw2mat(listw), 
##     w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl)
## 
## Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -25308   -7865     156       2    7869   48407 
## 
## Error variance parameters:
##     Estimate Std. Error t-value Pr(>|t|)   
## phi 6.968152   2.539192  2.7442 0.006065 **
## rho 0.012661   0.227801  0.0556 0.955678   
## 
## Coefficients:
##                Estimate  Std. Error t-value Pr(>|t|)   
## (Intercept)  9.8642e+04  4.2480e+04  2.3221 0.020227 * 
## PPM         -3.1921e+03  1.2246e+03 -2.6066 0.009145 **
## RRLS        -4.0608e+03  4.4794e+03 -0.9066 0.364641   
## KpdtPddk     9.2591e-02  1.3259e+00  0.0698 0.944327   
## PerLansia   -1.0368e+03  4.8863e+02 -2.1219 0.033846 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan hasil output dari pemodelan regresi data panel spasial REM, SAR, dan SEM, didapatkan model regresi sebagai berikut:

Model REM (Non-spatial Model): \(\hat Y = 9.8338 \times 10^4 - 3.1931 \times 10^3 \text{ (PPM)} - 4.0197 \times 10^3 \text{ (RRLS)} + 8.0816 \times 10^{-2} \text{ (KpdtPddk)} - 1.0372 \times 10^3 \text{ (PerLansia)} + \epsilon_i\)

Model SAR (Spatial Autoregressive): \(\hat Y = 9.7681 \times 10^4 - 3.214 \times 10^3 \text{ (PPM)}_i - 3.387 \times 10^3 \text{ (RRLS)}_i - 1.9971 \times 10^{-1} \text{ (KpdtPddk)}_i - 1.040 \times 10^3 \text{ (PerLansia)}_i -0.13171\text{ }WY_i + \epsilon_i\)

Model SEM (Spatial Error Model): \(\hat Y = 9.8642 \times 10^4 - 3.1921 \times 10^3 \text{ (PPM)}_i - 4.0608 \times 10^3 \text{ (RRLS)}_i + 9.2591 \times 10^{-2} \text{ (KpdtPddk)}_i - 1.0368 \times 10^3 \text{ (PerLansia)}_i + \hat u_i\) Dengan:

\(\hat u_i = 0.0127 \text{ }Wu_i + \epsilon_i\)

# Menghitung 4 varian LM Test

# Hitung satu per satu
lm_lag     <- slmtest(model_rem, listw = lw_full, test = "lml", type = "bp")
lm_err     <- slmtest(model_rem, listw = lw_full, test = "lme", type = "bp")
rlm_lag    <- slmtest(model_rem, listw = lw_full, test = "rlml", type = "bp")
rlm_err    <- slmtest(model_rem, listw = lw_full, test = "rlme", type = "bp")

# Buat Data Frame Gabungan
hasil_lm <- data.frame(
  Uji = c("LM Spatial Lag (RSlag)", 
          "LM Spatial Error (RSerr)", 
          "Robust LM Lag (RSadjlag)", 
          "Robust LM Error (RSadjerr)"),
  Statistic = c(lm_lag$statistic, 
                lm_err$statistic, 
                rlm_lag$statistic, 
                rlm_err$statistic),
  P_Value = c(lm_lag$p.value, 
              lm_err$p.value, 
              rlm_lag$p.value, 
              rlm_err$p.value)
)

# Tambahkan Kesimpulan Otomatis
hasil_lm$Kesimpulan <- ifelse(hasil_lm$P_Value < 0.05, "Signifikan", "Tidak Signifikan")

# Tampilkan Tabel
print(hasil_lm, row.names = FALSE)

##                         Uji Statistic   P_Value       Kesimpulan
##      LM Spatial Lag (RSlag) 0.6418867 0.4230284 Tidak Signifikan
##    LM Spatial Error (RSerr) 0.4425679 0.5058856 Tidak Signifikan
##    Robust LM Lag (RSadjlag) 0.4579695 0.4985745 Tidak Signifikan
##  Robust LM Error (RSadjerr) 0.2586508 0.6110486 Tidak Signifikan

# Perhitungan AIC dan BIC per model
model_rem_ml <- spml(rumus, 
                     data = pdata, 
                     listw = lw_full, 
                     model = "random", 
                     spatial.error = "none", 
                     lag = FALSE)

# 2. Fungsi Hitung Manual (Tetap diperlukan)
hitung_kriteria <- function(model, label, is_spatial = FALSE) {
  ll <- model$logLik
  k <- length(model$coefficients) + 1 + ifelse(is_spatial, 1, 0) # +1 utk sigma^2, +1 lagi jika spasial
  n <- nrow(model$model)
  
  return(c(Model = label, 
           AIC = round(-2 * ll + 2 * k, 4), 
           BIC = round(-2 * ll + k * log(n), 4)))
}

# 3. Langsung Gabung Jadi Tabel
tabel_aic_bic <- data.frame(rbind(
  hitung_kriteria(model_rem_ml, "REM (Non-Spasial)", is_spatial = FALSE),
  hitung_kriteria(model_sar,    "SAR - REM",         is_spatial = TRUE),
  hitung_kriteria(model_sem,    "SEM - REM",         is_spatial = TRUE)
))

# 4. Rapikan Format Angka
tabel_aic_bic$AIC <- as.numeric(tabel_aic_bic$AIC)
tabel_aic_bic$BIC <- as.numeric(tabel_aic_bic$BIC)

# 5. Tandai yang Terbaik (Minimum)
tabel_aic_bic$Ket_AIC <- ifelse(tabel_aic_bic$AIC == min(tabel_aic_bic$AIC), "Terbaik", "")
tabel_aic_bic$Ket_BIC <- ifelse(tabel_aic_bic$BIC == min(tabel_aic_bic$BIC), "Terbaik", "")

# Cetak Hasil
print(tabel_aic_bic, row.names = FALSE)

##              Model      AIC      BIC Ket_AIC Ket_BIC
##  REM (Non-Spasial) 1700.948 1715.315 Terbaik Terbaik
##          SAR - REM 1702.184 1718.945                
##          SEM - REM 1702.945 1719.706

Dengan taraf signifikansi sebesar 5% (0.05), dapat diketahui bahwa semua uji Lagrange Multiplier pada perhitungan Lag dan Error menghasilkan p-value lebih besar dari taraf signifikansi, baik nilai normal ataupun adjusted. Selain itu, didapatkan juga bahwa nilai AIC dan BIC yang dihasilkan juga menunjukkan bahwa nilai AIC dan BIC terkecil ada pada model REM. Dengan begitu, Random Effect Model menjadi model yang terbaik sebagai acuan untuk pemodelan spasial untuk kasus Diabetes Mellitus di Jawa Barat.

4.4. Interpolasi Spasial

Tahun 2022

# 1. DEFINISI FUNGSI
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]; z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

# 2. MEMBUAT 'df' (DARI DATA JABAR)
coords <- st_coordinates(st_centroid(jabar_data))
df <- data.frame(
  id = jabar_data$Kab_Kota,
  x  = coords[, "X"],
  y  = coords[, "Y"],
  z  = jabar_data$DM_2022
)

# 3. MEMBUAT GRID 'G'
xg <- seq(min(df$x), max(df$x), length.out = 100) 
yg <- seq(min(df$y), max(df$y), length.out = 100)
G <- expand.grid(x = xg, y = yg)

# 4. EKSEKUSI HITUNGAN (4 VARIASI POWER)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), 
                  G$x, G$y)
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), 
                 G$x, G$y)
G$z_p3 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 3), 
                 G$x, G$y)
G$z_p4 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 4), 
                 G$x, G$y)

# 5. VISUALISASI (GABUNGAN 4 GAMBAR)
# --- Langkah Tambahan: CLIPPING ---
G_sf <- st_as_sf(G, coords = c("x", "y"), crs = st_crs(jabar_data), remove = FALSE)
idx_inside <- st_intersects(G_sf, st_union(jabar_data), sparse = FALSE)
G_clipped <- G[idx_inside, ]

# --- Langkah A: Susun Data (Gabung 4 Model) ---
plot_data <- rbind(
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p1, Model = "IDW Power = 1"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p2,  Model = "IDW Power = 2"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p3,  Model = "IDW Power = 3"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p4,  Model = "IDW Power = 4")
)

# --- Langkah B: Plotting 2x2 ---
ggplot(plot_data, aes(x = x, y = y, fill = Prediksi)) +
  geom_tile() +
  scale_fill_gradientn(colours = c("black", "darkblue", "blue", "purple", "magenta", "orange", "yellow")) +

  # Facet wrap membagi jadi 2 kolom (hasilnya otomatis 2 baris karena ada 4 gambar)
  facet_wrap(~Model, ncol = 2) + 
  coord_equal() +
  labs(title = "Perbandingan Parameter Power (IDW) - Tahun 2022",
       x = NULL, y = NULL, fill = "Kasus DM") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "right"
  )

Tahun 2023

# 1. DEFINISI FUNGSI
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]; z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

# 2. MEMBUAT 'df' (DARI DATA JABAR)
coords <- st_coordinates(st_centroid(jabar_data))
df <- data.frame(
  id = jabar_data$Kab_Kota,
  x  = coords[, "X"],
  y  = coords[, "Y"],
  z  = jabar_data$DM_2023
)

# 3. MEMBUAT GRID 'G'
xg <- seq(min(df$x), max(df$x), length.out = 100) 
yg <- seq(min(df$y), max(df$y), length.out = 100)
G <- expand.grid(x = xg, y = yg)

# 4. EKSEKUSI HITUNGAN (4 VARIASI POWER)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), 
                  G$x, G$y)
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), 
                 G$x, G$y)
G$z_p3 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 3), 
                 G$x, G$y)
G$z_p4 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 4), 
                 G$x, G$y)

# 5. VISUALISASI (GABUNGAN 4 GAMBAR)
# --- Langkah Tambahan: CLIPPING ---
G_sf <- st_as_sf(G, coords = c("x", "y"), crs = st_crs(jabar_data), remove = FALSE)
idx_inside <- st_intersects(G_sf, st_union(jabar_data), sparse = FALSE)
G_clipped <- G[idx_inside, ]

# --- Langkah A: Susun Data (Gabung 4 Model) ---
plot_data <- rbind(
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p1, Model = "IDW Power = 1"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p2,  Model = "IDW Power = 2"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p3,  Model = "IDW Power = 3"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p4,  Model = "IDW Power = 4")
)

# --- Langkah B: Plotting 2x2 ---
ggplot(plot_data, aes(x = x, y = y, fill = Prediksi)) +
  geom_tile() +
  scale_fill_gradientn(colours = c("black", "darkblue", "blue", "purple", "magenta", "orange", "yellow")) +

  # Facet wrap membagi jadi 2 kolom (hasilnya otomatis 2 baris karena ada 4 gambar)
  facet_wrap(~Model, ncol = 2) + 
  coord_equal() +
  labs(title = "Perbandingan Parameter Power (IDW) - Tahun 2023",
       x = NULL, y = NULL, fill = "Kasus DM") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "right"
  )

Tahun 2024

# 1. DEFINISI FUNGSI
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]; z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

# 2. MEMBUAT 'df' (DARI DATA JABAR)
coords <- st_coordinates(st_centroid(jabar_data))
df <- data.frame(
  id = jabar_data$Kab_Kota,
  x  = coords[, "X"],
  y  = coords[, "Y"],
  z  = jabar_data$DM_2024
)

# 3. MEMBUAT GRID 'G'
xg <- seq(min(df$x), max(df$x), length.out = 100) 
yg <- seq(min(df$y), max(df$y), length.out = 100)
G <- expand.grid(x = xg, y = yg)

# 4. EKSEKUSI HITUNGAN (4 VARIASI POWER)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), 
                  G$x, G$y)
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), 
                 G$x, G$y)
G$z_p3 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 3), 
                 G$x, G$y)
G$z_p4 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 4), 
                 G$x, G$y)

# 5. VISUALISASI (GABUNGAN 4 GAMBAR)
# --- Langkah Tambahan: CLIPPING ---
G_sf <- st_as_sf(G, coords = c("x", "y"), crs = st_crs(jabar_data), remove = FALSE)
idx_inside <- st_intersects(G_sf, st_union(jabar_data), sparse = FALSE)
G_clipped <- G[idx_inside, ]

# --- Langkah A: Susun Data (Gabung 4 Model) ---
plot_data <- rbind(
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p1, Model = "IDW Power = 1"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p2,  Model = "IDW Power = 2"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p3,  Model = "IDW Power = 3"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p4,  Model = "IDW Power = 4")
)

# --- Langkah B: Plotting 2x2 ---
ggplot(plot_data, aes(x = x, y = y, fill = Prediksi)) +
  geom_tile() +
  scale_fill_gradientn(colours = c("black", "darkblue", "blue", "purple", "magenta", "orange", "yellow")) +

  # Facet wrap membagi jadi 2 kolom (hasilnya otomatis 2 baris karena ada 4 gambar)
  facet_wrap(~Model, ncol = 2) + 
  coord_equal() +
  labs(title = "Perbandingan Parameter Power (IDW) - Tahun 2024",
       x = NULL, y = NULL, fill = "Kasus DM") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "right"
  )

Rata-rata

# --- 1. PERSIAPAN DATA: HITUNG RATA-RATA (2022-2024) ---
# Pastikan nama kolom di data Anda sesuai (misal: DM_2022, DM_2023, DM_2024)
# Ambil data angkanya saja, buang geometry (peta)-nya
data_polos <- st_drop_geometry(jabar_data[, c("DM_2022", "DM_2023", "DM_2024")])

# Baru hitung rata-ratanya
jabar_data$DM_Rata <- rowMeans(data_polos, na.rm = TRUE)

# --- 2. DEFINISI FUNGSI IDW (Sama seperti sebelumnya) ---
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]; z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

# --- 3. MEMBUAT 'df' DENGAN DATA RATA-RATA ---
coords <- st_coordinates(st_centroid(jabar_data))

df <- data.frame(
  id = jabar_data$Kab_Kota,
  x  = coords[, "X"],
  y  = coords[, "Y"],
  z  = jabar_data$DM_Rata  # <--- BAGIAN INI DIGANTI JADI KOLOM RATA-RATA
)

# --- 4. MEMBUAT GRID 'G' ---
xg <- seq(min(df$x), max(df$x), length.out = 100) 
yg <- seq(min(df$y), max(df$y), length.out = 100)
G <- expand.grid(x = xg, y = yg)

# --- 5. EKSEKUSI HITUNGAN (4 VARIASI POWER) ---
# (Kode ini tetap sama, dia akan otomatis pakai 'z' baru yang isinya rata-rata)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), 
                 G$x, G$y)
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), 
                 G$x, G$y)
G$z_p3 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 3), 
                 G$x, G$y)
G$z_p4 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 4), 
                 G$x, G$y)

# --- 6. VISUALISASI ---
# Clipping agar bentuknya sesuai peta Jabar
G_sf <- st_as_sf(G, coords = c("x", "y"), crs = st_crs(jabar_data), remove = FALSE)
idx_inside <- st_intersects(G_sf, st_union(jabar_data), sparse = FALSE)
G_clipped <- G[idx_inside, ]

# Susun Data Plotting
plot_data <- rbind(
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p1, Model = "IDW Power = 1"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p2,  Model = "IDW Power = 2"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p3,  Model = "IDW Power = 3"),
  data.frame(x = G_clipped$x, y = G_clipped$y, Prediksi = G_clipped$z_p4,  Model = "IDW Power = 4")
)

# Plotting Final
ggplot(plot_data, aes(x = x, y = y, fill = Prediksi)) +
  geom_tile() +
  scale_fill_gradientn(colours = c("black", "darkblue", "blue", "purple", "magenta", "orange", "yellow")) +
  facet_wrap(~Model, ncol = 2) + 
  coord_equal() +
  labs(title = "Peta Interpolasi Risiko DM (Rata-rata 2022-2024)", # Judul disesuaikan
       subtitle = "Metode IDW dengan Variasi Parameter Power",
       x = NULL, y = NULL, fill = "Rata-rata Kasus") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "right"
  )

Berdasarkan output yang dihasilkan dengan metode Inverse Distannce Weighting (IWD), didapatkan empat hasil peta interpolasi dengan “power” yang berbeda, yaitu 1, 2, 3, dan 4. Peta tersebut menunjukkan adanya pola spasial yang cukup konsisten, terutama pada bagian daerah kiri atas dan kanan bawah, sesuai dengan hasil uji autokorelasi spasial lokal yang mendeteksi kota dan kabupaten dengan autokorelasi spasial pada daerah yang cenderung mirip, yaitu di sekitar Kabupaten Bekasi, Kota Bekasi, Kota Depok, Kota Bogor, dan Kabupaten Ciamis.

Selain itu, adanya perbedaan power yang digunakan sangat memengaruhi hasil peta interpolasi, dengan p = 1 yang menghasilkan peta cenderung halus dan p = 4 yang menunjukkan beberapa titik tertentu menjadi pusat gradasi warna.

4.5. Interpretasi Hasil

Dari hasil pengujian yang dilakukan, pemetaan kasus Diabetes Mellitus di Jawa Barat pada tahun 2022 sampai 2024 memperlihatkan pola pengklasteran yang kuat. Hal tersebut juga didukung dengan nilai Indeks Moran yang positif dan nilai p-value yang lebih kecil dari taraf signifikansi (5%), serta diperkuat dengan hasil Peta Cluster LISA yang konsisten mendeteksi zona hotspot (High-High) pada wilayah yang yang berdekatan dengan Kota Bekasi, Kota Depok, dan Kabupaten Bekasi, serta coldspot (Low-Low) pada wilayah bagian kiri bawah (seperti Kabupaten Ciamis). Pola ini didukung oleh peta hasil interpolasi dengan metode IDW yang memperlihatkan gradasi meluas dari wilayah yang memiliki kasus Diabetes Mellitus yang sangat tinggi, mengindikasikan bahwa secara visual, wilayah dengan kasus tinggi cenderung bertetangga dengan wilayah yang juga memiliki kasus tinggi.

Adapun model yang terbaik untuk pemodelan kasus Diabetes Mellitus dengan 4 variabel independen yang digunakan (Persentase Penduduk Miskin, Persentase Lansia, Kepadatan Penduduk, dan Rata-rata Lama sekolah) adalah dengan Random Effect Model (REM) Non-Spasial. Hal tersebut didapatkan dari hasil pada Uji Hausman yang mendukung penggunaan REM, serta Uji Lagrange Multiplier (LM) pada penentuan model spasial yang menunjukkan bahwa efek spasial (baik Lag maupun Error) ternyata tidak signifikan secara statistik setelah variabel independen dimasukkan. Kesimpulan ini diperkuat oleh nilai AIC/BIC model REM Non-Spasial yang lebih rendah (terbaik) dibandingkan model SAR atau SEM. Meskipun begitu, nilai R-squared yang dihasilkan dari REM cenderung kecil, menandakan model tidak cukup baik dalam menggambarkan kasus Diabetes Mellitus. Faktor seperti keterbatasan variabel yang bisa diakses menjadi salah satu hal yang membuat pemodelan kasus Diabetes Mellitus menjadi susah untuk dilakukan oleh orang umum.

Bab 5. Kesimpulan dan Saran

5.1. Kesimpulan

Penyakit Diabetes Mellitus merupakan salah satu hal yang sangat membebankan Provinsi Jawa Barat, meskipun bukan tergolong penyakit menular. Hasil uji menunjukkan bahwa terdapat autokorelasi spasial yang signifikan pada persebaran jumlah kasus Diabetes Mellitus cenderung menyebar, menandakan bahwa kasus penyakit Diabetes Mellitus cenderung mengumpul pada suatu bagian wilayah. Melakukan kajian lebih dalam spasial, model regresi, dan pemetaan interpolasi bisa menjadi bentuk penelitian yang baik untuk mengenali kasus Diabetes Mellitus lebih baik di Jawa Barat, meskipun pemodelan regresi yang dihasilkan dari penelitian ini masih belum bisa dianggap sebagai sebuah model yang menggambarkan kasus dengan baik. Hal tersebut dikarenakan nilai R-squared yang didapatkan dari Random Effect Model hanya sebesar 0.10871 menandakan bahwa variabel independen yang digunakan masih belum dapat mendeskripsikan atau memprediksi jumlah kasus Diabetes Mellitus pada suatu wilayah dengan baik/

5.2. Saran

Adapun saran untuk analisis kasus Diabetes Mellitus kedepannya adalah dengan menggunakan data dengan jumlah periode lebih banyak untuk mengenal lebih baik naik turunnya jumlah kasus dari waktu ke waktu, serta menggunakan variabel independen yang lebih memungkinkan untuk berkorelasi secara signifikan terhadap kasus Diabetes Mellitus.

Daftar Pustaka

Kementrian Kesehatan Republik Indonesia. (2025). Review kebijakan Diabetes Melitus berbasis transformasi sistem Kesehatan dan outlook 2025. https://lms.kemkes.go.id/courses/acdcbe95-9e14-4b6b-9ffd-305e13989c8e (Terakhir diakses 20 Desember 2025, 13:29 WIB)

Lovic, D., Piperidou, A., Zografou, I., Grassos, H., Pittaras, A., & Manolis, A. (2020). The Growing Epidemic Of Diabetes Mellitus.. Current vascular pharmacology. https://doi.org/10.2174/1570161117666190405165911.

Xu, L., Sun, X., Wang, N., Wang, Y., Li, Y., & Zhang, C. (2025). Risk factor assessment of prediabetes and diabetes based on epidemic characteristics in new urban areas: a retrospective and a machine learning study. Scientific Reports, 15. https://doi.org/10.1038/s41598-025-88073-6.

Braver, N., Lakerveld, J., Rutters, F., Schoonmade, L., Brug, J., Brug, J., Beulens, J., & Beulens, J. (2017). Built environmental characteristics and diabetes: a systematic review and meta-analysis. BMC Medicine, 16. https://doi.org/10.1186/s12916-017-0997-z.

Syavera, V., Syazali, M. (2024). Peta Risiko Diabetes Melitus di Jawa Barat Tahun 2019-2023 dengan Pemodelan Spatio-Temporal. SEHATRAKYAT, Vol. 3 No. 4. https://journal.yp3a.org/index.php/sehatrakyat/article/download/3296/1230

Open Data Jabar. (2025). Jumlah Penderita Diabetes Melitus Berdasarkan Kabupaten/Kota di Jawa Barat. [Dataset]. https://opendata.jabarprov.go.id/id/dataset/jumlah-penderita-diabetes-melitus-berdasarkan-kabupatenkota-di-jawa-barat

Open Data Jabar. (2025). Persentase Penduduk Miskin Berdasarkan Kabupaten/Kota di Jawa Barat. [Dataset]. https://opendata.jabarprov.go.id/id/dataset/persentase-penduduk-miskin-berdasarkan-kabupatenkota-di-jawa-barat

Open Data Jabar. (2025). Kepadatan Penduduk Berdasarkan Kabupaten/Kota di Jawa Barat. [Dataset]. https://opendata.jabarprov.go.id/id/dataset/kepadatan-penduduk-berdasarkan-kabupatenkota-di-jawa-barat

Open Data Jabar. (2025). Rata-Rata Lama Sekolah Berdasarkan Kabupaten/Kota Jawa Barat. [Dataset]. https://opendata.jabarprov.go.id/id/dataset/rata-rata-lama-sekolah-berdasarkan-kabupatenkota-di-jawa-barat

BPS Jawa Barat. (2023). Profil Penduduk Lanjut Usia Provinsi Jawa Barat 2022. https://jabar.bps.go.id/id/publication/2023/12/29/542e5fe5c1319fcefcf144ca/profil-penduduk-lanjut-usia-provinsi-jawa-barat-2022.html

BPS Jawa Barat. (2024). Profil Penduduk Lanjut Usia Provinsi Jawa Barat 2023. https://jabar.bps.go.id/id/publication/2024/12/23/18a003afe3f3acc802f7ee83/profil-penduduk-lanjut-usia-provinsi-jawa-barat-2023.html

BPS Jawa Barat. (2025). Profil Penduduk Lanjut Usia Provinsi Jawa Barat 2024. https://jabar.bps.go.id/id/publication/2025/05/23/db055b2f3428b1667f68ce8e/profil-penduduk-lanjut-usia-provinsi-jawa-barat-2024.html

Lampiran

Tautan Dashboard: https://nanaurum.shinyapps.io/Spasial_Diabetes_Mellitus/

Tautan YouTube: https://youtu.be/ipT4NVqePho