ANALISIS EPIDEMIOLOGI PENYAKIT DIARE DI PROVINSI JAWA BARAT TAHUN 2024 DENGAN PENDEKATAN ANALISIS SPASIAL DAN REGRESI NEGATIVE BINOMIAL

---

Abstrak

Diare merupakan penyakit berbasis lingkungan yang menjadi beban kesehatan masyarakat signifikan di Provinsi Jawa Barat. Studi ini bertujuan menganalisis pola epidemiologi penyakit diare secara cross-sectional pada 27 kabupaten/kota di Jawa Barat tahun 2024, meliputi perhitungan Incidence Rate (IR) dan Standardized Incidence Ratio (SIR), analisis autokorelasi spasial menggunakan Global Moran’s I dan klaster LISA, serta identifikasi determinan jumlah kasus menggunakan regresi Negative Binomial. Hasil menunjukkan total 672.287 kasus diare dengan IR rata-rata Jawa Barat 133,54 per 10.000 penduduk. Kota Cirebon mencatatkan IR tertinggi (442,83 per 10.000, SIR = 3,32), diikuti Kota Bogor (321,38, SIR = 2,41) dan Kota Sukabumi (311,29, SIR = 2,33). Uji Moran’s I (I = 0,060, p = 0,226) tidak menunjukkan autokorelasi spasial signifikan. Regresi Negative Binomial Model 3 terpilih berdasarkan AIC terendah (564,97) mengidentifikasi ketersediaan fasilitas kesehatan (IRR = 0,724, p < 0,001) dan IPM (IRR = 0,700, p = 0,021) sebagai faktor protektif signifikan, sedangkan kepadatan penduduk (IRR = 1,491, p = 0,008) merupakan faktor risiko. Overdispersi data dikonfirmasi dengan VMR = 6.110,70 dan dispersion estimate = 1.656,22 (p < 0,001), menjustifikasi penggunaan model Negative Binomial.

Kata kunci: diare, incidence rate, SIR, Moran’s I, LISA, Regresi Negative Binomial

---

1 Pendahuluan

1.1 Latar Belakang

Diare adalah kondisi di mana seseorang mengalami buang air besar lebih dari tiga kali sehari dengan konsistensi cair atau lembek. Secara global, diare merupakan penyebab kedua kematian pada anak di bawah lima tahun setelah pneumonia (WHO, 2023). Di Indonesia, diare termasuk dalam sepuluh besar penyakit dengan angka kesakitan tertinggi dan terus menjadi masalah kesehatan masyarakat prioritas, terutama di wilayah dengan akses sanitasi dan air bersih yang masih terbatas.

Provinsi Jawa Barat sebagai provinsi dengan jumlah penduduk terbesar di Indonesia (±50,3 juta jiwa pada 2024) menghadapi tantangan besar dalam pengendalian penyakit diare. Data Open Data Jabar menunjukkan bahwa pada tahun 2024 tercatat sebanyak 672.287 kasus diare di 27 kabupaten/kota, dengan variasi kasus yang signifikan antarwilayah—mulai dari 4.096 kasus di Kota Banjar hingga 56.802 kasus di Kabupaten Bogor. Namun, jumlah kasus absolut tidak mencerminkan risiko sebenarnya; diperlukan analisis berbasis denominan populasi untuk perbandingan yang bermakna.

Berbagai faktor determinan diare telah diidentifikasi dalam literatur, meliputi akses sanitasi dasar, ketersediaan air bersih, kepadatan penduduk, kondisi sosial ekonomi yang direpresentasikan oleh Indeks Pembangunan Manusia (IPM) dan garis kemiskinan, serta ketersediaan fasilitas kesehatan. Pemahaman tentang faktor-faktor yang mempengaruhi beban diare secara simultan dalam konteks multivariat diperlukan untuk perumusan intervensi yang tepat sasaran.

Pendekatan spasial epidemiologi juga penting untuk mengidentifikasi apakah terdapat pola klaster geografis yang memerlukan perhatian khusus. Analisis Moran’s I global dan LISA (Local Indicators of Spatial Association) memungkinkan identifikasi wilayah-wilayah dengan beban penyakit yang tidak acak secara spasial, yang dapat menjadi basis penetapan prioritas program kesehatan masyarakat.

1.2 Rumusan Masalah

1. Bagaimana distribusi kasus diare di Jawa Barat tahun 2024?
2. Wilayah mana yang memiliki risiko relatif tertinggi?
3. Apakah terdapat pola spasial kasus diare?
4. Faktor apa saja yang berpengaruh terhadap jumlah kasus diare?

1.3 Tujuan Penelitian

1. Menghitung ukuran epidemiologi diare tahun 2024.
2. Mengidentifikasi wilayah dengan risiko tinggi berdasarkan SIR.
3. Menguji keberadaan autokorelasi spasial menggunakan Moran’s I dan LISA.
4. Menganalisis faktor-faktor yang memengaruhi jumlah kasus diare menggunakan regresi Negative Binomial.

---

2 Data dan Metode

2.1 Sumber Data

Penelitian ini menggunakan data sekunder lintas-seksional (cross-sectional) pada tahun 2024 dari tiga sumber utama:

1. Open Data Jabar (opendata.jabarprov.go.id): jumlah kasus diare per kabupaten/kota, jumlah fasilitas kesehatan, persentase akses air bersih, persentase akses sanitasi, dan kepadatan penduduk.
2. Badan Pusat Statistik (BPS) Provinsi Jawa Barat: jumlah penduduk, garis kemiskinan per kapita per bulan, dan Indeks Pembangunan Manusia (IPM).
3. GADM (Database of Global Administrative Areas): batas administratif kabupaten/kota Jawa Barat dalam format GeoJSON untuk keperluan analisis spasial dan pemetaan.

Tanggal akses data: Juni 2026. Seluruh data merupakan data publik yang dapat diverifikasi dan diakses kembali melalui tautan sumber resmi.

2.2 Unit Analisis dan Variabel Penelitian

Unit analisis adalah kabupaten/kota (n = 27) di Provinsi Jawa Barat pada tahun 2024. Tabel 1 merangkum variabel yang digunakan.

Tabel 1. Unit Analisis dan Variabel

Variabel	Tipe	Satuan	Peran
Kasus Diare	Respons (count)	Orang/tahun	Variabel dependen regresi NB
Jumlah Penduduk	Denominator	Jiwa	Denominan IR, SIR, dan offset regresi
Fasilitas Kesehatan	Prediktor (Log)	Unit	Proksi ketersediaan layanan kesehatan
Akses Air Bersih	Prediktor	%	Rumah tangga dengan akses air minum layak
Akses Sanitasi	Prediktor	%	Rumah tangga dengan sanitasi layak
IPM	Prediktor	Indeks	Indeks Pembangunan Manusia (0–100)
Garis Kemiskinan	Prediktor (Log)	Rp/kapita/bln	Pengeluaran minimum per kapita
Kepadatan Penduduk	Prediktor (Log)	Jiwa/km²	Intensitas kontak antar-individu

2.3 Pengolahan Data

Data diolah menggunakan perangkat lunak R dengan paket utama: readxl, dplyr, MASS (glm.nb), AER (dispersiontest), car (vif), sf, dan spdep. Tahapan pengolahan meliputi: a. Pembacaan data Excel dengan penugasan nama kolom secara manual berdasarkan posisi (14 kolom) untuk menghindari dependensi pada versi paket janitor. b. Kolom air_bersih disimpan sebagai karakter di Excel; nilai “-” diubah ke NA kemudian diimputasi dengan nilai median. c. Transformasi log natural untuk variabel faskes, kemiskinan, dan kepadatan yang berdistribusi positif skewed, guna memenuhi asumsi linearitas log-link pada model NB. d. Standardisasi (z-score) seluruh prediktor agar koefisien regresi dapat dibandingkan lintas variabel.

2.4 Ukuran Epidemiologi

2.4.1 Incidence Rate (IR)

Incidence Rate (IR), dihitung sebagai jumlah kasus diare dibagi jumlah penduduk dikalikan 10.000, menghasilkan satuan kasus per 10.000 penduduk per tahun:

\[ IR = \frac{Kasus Diare}{Penduduk}\times 10.000 \]

2.4.2 Standardized Incidence Ratio (SIR)

Standardized Incidence Ratio (SIR), dihitung sebagai rasio kasus teramati terhadap kasus yang diharapkan, di mana kasus yang diharapkan diperoleh dengan mengaplikasikan IR referensi Jawa Barat 2024 (133,54 per 10.000) terhadap penduduk masing-masing wilayah:

\[ SIR = \frac{Observed}{Expected} \]

dengan

\[ Expected = Penduduk_i \times Rate_{JawaBarat} \]

SIR > 1 menunjukkan beban diare melebihi rata-rata provinsi; SIR < 1 menunjukkan beban di bawah rata-rata. SIR dikategorikan menjadi: Sangat Tinggi (≥ 2,0), Tinggi (1,0–2,0), Rendah (0,5–1,0), dan Sangat Rendah (< 0,5).

2.5 Analisis Spasial

Analisis spasial dilakukan dengan pendekatan berbasis pembobot Queen Contiguity menggunakan paket spdep. Dua analisis dilakukan:

1. Global Moran’s I: mengukur derajat autokorelasi spasial IR secara keseluruhan. Nilai I berkisar antara -1 (dispersi sempurna) hingga +1 (klaster sempurna). Uji signifikansi menggunakan pendekatan randomisasi (p < 0,05).

2. Local Indicators of Spatial Association (LISA): mengidentifikasi klaster spasial lokal pada level wilayah. Klaster dikategorikan sebagai High-High, Low-Low, High-Low, atau Low-High berdasarkan nilai IR terstandarisasi dan spatial lag-nya, dengan threshold signifikansi p < 0,05.

2.6 Regresi Negative Binomial

Regresi binomial negatif adalah salah satu alternatif dari model regresi Poisson yang mengalami overdispersi. Regresi Binomial Negatif adalah model penerapan dari Generalized Linear Model (GLM) yang menghubungkan antara variabel dependen dengan variabel independen yang dapat diterapkan pada kondisi equdispersi atau overdispersi

Data kasus diare merupakan data cacah (count data) dengan overdispersi yang dikonfirmasi melalui Variance-Mean Ratio (VMR) dan uji formal dispersiontest dari paket AER. Model Negative Binomial (NB) dipilih karena mengatasi overdispersi melalui parameter dispersi theta (θ), yang tidak dapat ditangani oleh model Poisson. Tiga model bertingkat diestimasi dengan offset log(jumlah penduduk) untuk mengoreksi perbedaan ukuran populasi:

1. Model 1: log(faskes) + air_bersih + sanitasi + IPM

\[ \quad \log(\mu_i) = \ln(\text{penduduk}_i) + \beta_0 + \beta_1 \text{Faskes}_i + \beta_2 \text{AirBersih}_i + \beta_3 \text{Sanitasi}_i + \beta_4 \text{IPM}_i\]

2. Model 2: Model 1 + log(garis kemiskinan)

\[ \quad \log(\mu_i) = \ln(\text{penduduk}_i) + \beta_0 + \beta_1 \text{Faskes}_i + \beta_2 \text{AirBersih}_i + \beta_3 \text{Sanitasi}_i + \beta_4 \text{IPM}_i + \beta_5 \text{Kemiskinan}_i\]

3. Model 3: Model 2 + log(kepadatan)

\[ \quad \log(\mu_i) = \ln(\text{penduduk}_i) + \beta_0 + \beta_1 \text{Faskes}_i + \beta_2 \text{AirBersih}_i + \beta_3 \text{Sanitasi}_i + \beta_4 \text{IPM}_i + \beta_5 \text{Kemiskinan}_i + \beta_6 \text{Kepadatan}_i\]

Pemilihan model berdasarkan Akaike Information Criterion (AIC) terkecil. Model-model dibandingkan secara formal menggunakan Likelihood-Ratio Test (LRT). Goodness of fit diuji dengan uji chi-kuadrat berbasis devians residual. Output utama berupa Incidence Rate Ratio (IRR = exp(β)) dengan 95% Confidence Interval. Multikolinearitas diperiksa menggunakan Variance Inflation Factor (VIF) melalui model OLS pengganti.

---

3 Hasil Analisis

3.1 Persiapan Data

setwd("C:/Users/HP/Downloads")

df_raw <- read_excel("dataset diare fix.xlsx") %>%
  clean_names()

# Rename ke nama pendek yang konsisten
df_raw <- df_raw %>%
  rename(
    nama_kab   = nama_kabupaten_kota,
    diare      = diare_orang,
    faskes     = faskes_unit,
    kemiskinan = garis_kemiskinan_rupiah_kapita_bulan,
    kepadatan  = kepadatan_penduduk_jiwa_per_kilometer_persegi,
    air_bersih = air_bersih_percent,
    sanitasi   = sanitasi_percent,
    ipm        = ipm_percent,
    penduduk   = jumlah_penduduk_jiwa
  ) %>%
  dplyr::select(
    nama_kab,
    tahun,
    diare,
    faskes,
    kemiskinan,
    kepadatan,
    air_bersih,
    sanitasi,
    ipm,
    penduduk
  )



df <- df_raw %>%
  mutate(
    air_bersih = trimws(as.character(air_bersih)),
    air_bersih = na_if(air_bersih, "-"),

    diare      = as.numeric(diare),
    faskes     = as.numeric(faskes),
    kemiskinan = as.numeric(kemiskinan),
    kepadatan  = as.numeric(kepadatan),
    air_bersih = as.numeric(air_bersih),
    sanitasi   = as.numeric(sanitasi),
    ipm        = as.numeric(ipm),
    penduduk   = as.numeric(penduduk),
    tahun      = as.integer(tahun)
  ) %>%
  filter(tahun == 2024)

df <- df %>%
  mutate(
    air_bersih = ifelse(
      is.na(air_bersih),
      median(air_bersih, na.rm = TRUE),
      air_bersih
    ),
    sanitasi = ifelse(
      is.na(sanitasi),
      median(sanitasi, na.rm = TRUE),
      sanitasi
    )
  )

df <- df %>%
  mutate(
    # Log transforms untuk variabel skewed kanan
    ln_penduduk   = log(penduduk),
    ln_faskes     = log1p(faskes),      # log(x+1) aman jika ada nol
    ln_kemiskinan = log(kemiskinan),
    ln_kepadatan  = log(kepadatan),
    
    # Standardisasi (mean=0, SD=1) agar koefisien dapat dibandingkan
    faskes_std     = as.numeric(scale(ln_faskes)),
    air_std        = as.numeric(scale(air_bersih)),
    sanitasi_std   = as.numeric(scale(sanitasi)),
    ipm_std        = as.numeric(scale(ipm)),
    kemiskinan_std = as.numeric(scale(ln_kemiskinan)),
    kepadatan_std  = as.numeric(scale(ln_kepadatan)))

df <- df %>%
  mutate(
    tipe = ifelse(str_detect(toupper(nama_kab), "^KOTA"), "Kota", "Kabupaten")
  )

# Data Spasial
path_shp  <- "D:/Batas/Batas/Batas_Kabupaten_Jawa_Barat.shp"

jabar_geo <- st_read(path_shp) 

## Reading layer `Batas_Kabupaten_Jawa_Barat' from data source 
##   `D:\Batas\Batas\Batas_Kabupaten_Jawa_Barat.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3703 ymin: -7.82099 xmax: 108.8468 ymax: -5.806538
## Geodetic CRS:  WGS 84

jabar_geo <- jabar_geo %>%
    mutate(
      nama_kab = toupper(WADMKK),
      nama_kab = case_when(
        str_detect(nama_kab, "^KOTA") ~ nama_kab,
        TRUE ~ paste("KABUPATEN", nama_kab)
      )
    )

# Standarisasi nama di df untuk join
df_spat <- df %>%
  mutate(nama_kab = toupper(nama_kab))

# Gabungkan shapefile dengan data
jabar <- left_join(jabar_geo, df_spat, by = "nama_kab")

---

3.2 Statistik Deskriptif

Tabel 2. Statistik Deskriptif Variabel Diare

variabel <- c(
  "diare",
  "penduduk",
  "faskes",
  "kemiskinan",
  "kepadatan",
  "air_bersih",
  "sanitasi",
  "ipm"
)

stat_desc <- data.frame(
  Variabel = c(
    "Kasus Diare",
    "Jumlah Penduduk",
    "Jumlah Faskes",
    "Garis Kemiskinan",
    "Kepadatan Penduduk",
    "Akses Air Bersih (%)",
    "Akses Sanitasi (%)",
    "IPM"
  ),
  Minimum = sapply(df[variabel], min, na.rm = TRUE),
  Mean    = sapply(df[variabel], mean, na.rm = TRUE),
  Median  = sapply(df[variabel], median, na.rm = TRUE),
  SD      = sapply(df[variabel], sd, na.rm = TRUE),
  Maksimum= sapply(df[variabel], max, na.rm = TRUE)
)

kable(stat_desc, row.names = FALSE)

Variabel	Minimum	Mean	Median	SD	Maksimum
Kasus Diare	4096.00	2.489952e+04	21014.00	1.233505e+04	56802.00
Jumlah Penduduk	209790.00	1.864637e+06	1884190.00	1.228053e+06	5682300.00
Jumlah Faskes	220.00	2.034815e+03	1758.00	1.329054e+03	5284.00
Garis Kemiskinan	393464.00	5.360651e+05	475046.00	1.254443e+05	843893.00
Kepadatan Penduduk	385.00	3.910926e+03	1468.00	4.668116e+03	15176.00
Akses Air Bersih (%)	85.37	9.509519e+01	96.76	4.221395e+00	99.85
Akses Sanitasi (%)	34.56	8.945037e+01	95.78	1.489451e+01	100.00
IPM	68.89	7.467778e+01	73.82	4.333325e+00	83.75

Interpretasi:

Jumlah kasus diare per wilayah sangat bervariasi (min = 4.096; maks = 56.802), dengan SD yang besar (12.335,05) relatif terhadap rerata (24.899,52). Variasi besar ini merupakan indikasi awal overdispersi yang menguat dugaan bahwa model Poisson tidak memadai. IR berkisar dari 62,62 hingga 442,83 per 10.000 penduduk rasio maks/min mencapai 7,07, menunjukkan ketimpangan beban penyakit yang substansial. Akses sanitasi menunjukkan variasi terbesar (SD = 14,89%) dengan nilai minimum 34,56% (Kota Tasikmalaya), mengindikasikan kesenjangan infrastruktur sanitasi yang masih lebar antarwilayah.

---

3.3 Ukuran Epidemiologi

3.3.1 Incidence Rate (IR) dan SIR

Tabel 3 menyajikan IR dan SIR seluruh 27 kabupaten/kota, diurutkan berdasarkan SIR tertinggi. IR referensi Jawa Barat 2024 adalah 133,54 per 10.000 penduduk (total 672.287 kasus dari 50.345.190 penduduk).

Tabel 3. Incidence Rate dan SIR Diare

# IR
rate_jabar <- sum(df$diare) / sum(df$penduduk)

df <- df %>%
  mutate(IR = (diare / penduduk) * 10000)

# Update jabar dengan IR
jabar <- left_join(jabar_geo, 
                   df %>% mutate(nama_kab = toupper(nama_kab)), 
                   by = "nama_kab")

# SIR

df <- df %>%
  mutate(
    expected = penduduk * rate_jabar,
    SIR      = diare / expected
  )

# Update jabar dengan SIR
jabar <- left_join(jabar_geo, 
                   df %>% mutate(nama_kab = toupper(nama_kab)), 
                   by = "nama_kab")

jabar <- jabar %>%
  mutate(
    kategori_SIR = case_when(
      SIR < 0.80 ~ "<0.80",
      SIR < 1.00 ~ "0.80-0.99",
      SIR < 1.20 ~ "1.00-1.19",
      SIR < 1.50 ~ "1.20-1.49",
      TRUE       ~ ">=1.50"
    )
  )

tabel_epid <- df %>%
  mutate(
    IR_round  = round(IR, 2),
    SIR_round = round(SIR, 3),
    Kategori  = case_when(
      SIR >= 2.0  ~ "Sangat Tinggi (>=2.0)",
      SIR >= 1.0  ~ "Tinggi (1.0-2.0)",
      SIR >= 0.5  ~ "Rendah (0.5-1.0)",
      TRUE        ~ "Sangat Rendah (<0.5)"
    )
  ) %>%
  dplyr::select(Wilayah = nama_kab, Tipe = tipe,
         Kasus = diare, `IR/10.000` = IR_round,
         `Kasus Exp.` = expected, SIR = SIR_round, Kategori) %>%
  mutate(`Kasus Exp.` = round(`Kasus Exp.`, 0)) %>%
  arrange(desc(SIR))

kable(tabel_epid)

Wilayah	Tipe	Kasus	IR/10.000	Kasus Exp.	SIR	Kategori
KOTA CIREBON	Kota	15271	442.83	4605	3.316	Sangat Tinggi (>=2.0)
KOTA BOGOR	Kota	34656	321.38	14400	2.407	Sangat Tinggi (>=2.0)
KOTA SUKABUMI	Kota	11385	311.29	4884	2.331	Sangat Tinggi (>=2.0)
KABUPATEN MAJALENGKA	Kabupaten	32382	239.42	18061	1.793	Tinggi (1.0-2.0)
KOTA TASIKMALAYA	Kota	17845	237.70	10025	1.780	Tinggi (1.0-2.0)
KOTA BANJAR	Kota	4096	195.24	2801	1.462	Tinggi (1.0-2.0)
KABUPATEN CIREBON	Kabupaten	46223	193.57	31888	1.450	Tinggi (1.0-2.0)
KABUPATEN PURWAKARTA	Kabupaten	20193	192.25	14026	1.440	Tinggi (1.0-2.0)
KABUPATEN PANGANDARAN	Kabupaten	8003	184.36	5797	1.381	Tinggi (1.0-2.0)
KABUPATEN SUBANG	Kabupaten	28398	170.75	22209	1.279	Tinggi (1.0-2.0)
KABUPATEN SUMEDANG	Kabupaten	18361	154.67	15852	1.158	Tinggi (1.0-2.0)
KOTA CIMAHI	Kota	9071	151.51	7995	1.135	Tinggi (1.0-2.0)
KABUPATEN CIAMIS	Kabupaten	18336	145.61	16815	1.090	Tinggi (1.0-2.0)
KABUPATEN SUKABUMI	Kabupaten	39941	141.23	37764	1.058	Tinggi (1.0-2.0)
KABUPATEN CIANJUR	Kabupaten	36080	139.58	34519	1.045	Tinggi (1.0-2.0)
KABUPATEN GARUT	Kabupaten	37640	138.54	36281	1.037	Tinggi (1.0-2.0)
KABUPATEN KUNINGAN	Kabupaten	16429	135.34	16210	1.013	Tinggi (1.0-2.0)
KOTA BANDUNG	Kota	32258	127.59	33760	0.956	Rendah (0.5-1.0)
KABUPATEN BANDUNG BARAT	Kabupaten	23322	123.78	25161	0.927	Rendah (0.5-1.0)
KABUPATEN TASIKMALAYA	Kabupaten	21014	109.40	25651	0.819	Rendah (0.5-1.0)
KOTA DEPOK	Kota	22917	105.92	28892	0.793	Rendah (0.5-1.0)
KOTA BEKASI	Kota	27749	104.95	35308	0.786	Rendah (0.5-1.0)
KABUPATEN BOGOR	Kabupaten	56802	99.96	75879	0.749	Rendah (0.5-1.0)
KABUPATEN BANDUNG	Kabupaten	36693	97.77	50117	0.732	Rendah (0.5-1.0)
KABUPATEN INDRAMAYU	Kabupaten	16757	87.55	25559	0.656	Rendah (0.5-1.0)
KABUPATEN KARAWANG	Kabupaten	19964	78.16	34110	0.585	Rendah (0.5-1.0)
KABUPATEN BEKASI	Kabupaten	20501	62.62	43718	0.469	Sangat Rendah (<0.5)

Interpretasi:

Tiga wilayah termasuk kategori Sangat Tinggi (SIR ≥ 2,0): Kota Cirebon (SIR = 3,32; IR = 442,83), Kota Bogor (SIR = 2,41; IR = 321,38), dan Kota Sukabumi (SIR = 2,33; IR = 311,29). Sebanyak 14 wilayah masuk kategori Tinggi (1,0 ≤ SIR < 2,0), 9 wilayah Rendah (0,5 ≤ SIR < 1,0), dan hanya 1 wilayah Sangat Rendah (Kabupaten Bekasi, SIR = 0,47). Dengan demikian, 17 dari 27 wilayah (63%) memiliki beban diare yang melebihi rata-rata provinsi.

Analisis berdasarkan tipe wilayah menunjukkan bahwa wilayah perkotaan (kota) memiliki rerata IR lebih tinggi (222,05 per 10.000) dibandingkan kabupaten (138,59 per 10.000). Uji Mann-Whitney menunjukkan perbedaan ini tidak signifikan secara statistik (U = 117,0; p = 0,068), kemungkinan akibat ukuran sampel yang kecil (9 kota vs 18 kabupaten).

3.3.2 Visualisasi Ukuran Epidemiologi

a. Incidence Rate (IR)

p_ir <- df %>%
  mutate(nama_kab = reorder(nama_kab, IR)) %>%
  ggplot(aes(x = IR, y = nama_kab, fill = tipe)) +
  geom_col(alpha = 0.85) +
  geom_vline(xintercept = rate_jabar * 10000,
             linetype = "dashed", color = "#f39c12", linewidth = 0.9) +
  annotate("text",
           x     = rate_jabar * 10000 + 5,
           y     = 2,
           label = paste0("IR Jabar\n", round(rate_jabar*10000, 1)),
           color = "#f39c12", size = 3, hjust = 0) +
  scale_fill_manual(values = c("Kabupaten" = "#3498db", "Kota" = "#e74c3c")) +
  scale_x_continuous(labels = comma) +
  labs(
    title    = "Incidence Rate Diare per 10.000 Penduduk",
    subtitle = "Jawa Barat 2024 | Garis oranye = IR rata-rata Provinsi",
    x        = "IR per 10.000 Penduduk",
    y        = NULL,
    fill     = "Tipe Wilayah"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title          = element_text(face = "bold", size = 12),
    panel.grid.major.y  = element_blank(),
    legend.position     = "bottom"
  )
print(p_ir)

Gambar 1. Diagram Batang Incidence Rate (IR) Diare

Interpretasi:

Berdasarkan Gambar 1, nilai Incidence Rate (IR) diare di kabupaten/kota Provinsi Jawa Barat tahun 2024 menunjukkan variasi yang cukup besar. Kota Cirebon memiliki IR tertinggi, diikuti oleh Kota Bogor dan Kota Sukabumi. Sementara itu, Kabupaten Bekasi memiliki IR terendah. Rata-rata IR Provinsi Jawa Barat sebesar 133,5 kasus per 10.000 penduduk, dengan beberapa kota berada jauh di atas nilai rata-rata tersebut. Secara umum, wilayah perkotaan cenderung memiliki IR yang lebih tinggi dibandingkan sebagian besar wilayah kabupaten.

jabar <- jabar %>%
  mutate(
    kategori_IR = cut(
      IR,
      breaks = c(-Inf, 100, 150, 200, 300, Inf),
      labels = c(
        "<100",
        "100-150",
        "151-200",
        "201-300",
        ">300"
      )
    )
  )

peta_ir <- ggplot(jabar) +
  
  geom_sf(
    aes(fill = kategori_IR),
    color = "white",
    linewidth = 0.3
  ) +
  
  scale_fill_brewer(
    palette = "YlOrRd",
    name = "IR\nper 10.000"
  ) +
  
  labs(
    title = "Peta Incidence Rate (IR) Diare",
    subtitle = "Kabupaten/Kota Jawa Barat Tahun 2024",
    caption = "IR = Kasus Diare per 10.000 Penduduk"
  ) +
  
  theme_minimal() +
  
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(
      face = "bold",
      hjust = 0.5
    )
  )

print(peta_ir)

Gambar 2. Peta Choropleth Incidence Rate (IR) Diare

Interpretasi:

Berdasarkan Gambar 2, sebaran spasial Incidence Rate (IR) diare di Jawa Barat tahun 2024 tidak merata antarwilayah. Wilayah dengan kategori IR sangat tinggi (>300 kasus per 10.000 penduduk) terkonsentrasi pada beberapa kota, terutama Kota Cirebon, Kota Bogor, dan Kota Sukabumi. Sebagian besar kabupaten berada pada kategori rendah hingga sedang. Pola ini menunjukkan adanya perbedaan tingkat kejadian diare antarwilayah yang mengindikasikan kemungkinan adanya pengaruh faktor-faktor lokal dan spasial.

---

b. 3.3.2 Standardized Incidence Ratio (SIR)

p_sir <- df %>%
  mutate(
    nama_kab = reorder(nama_kab, SIR),
    warna    = case_when(
      SIR >= 2   ~ "#c0392b",
      SIR >= 1   ~ "#e67e22",
      SIR >= 0.5 ~ "#27ae60",
      TRUE       ~ "#2980b9"
    )
  ) %>%
  ggplot(aes(x = SIR, y = nama_kab, fill = warna)) +
  geom_col(alpha = 0.85) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "black", linewidth = 1) +
  scale_fill_identity(
    guide  = "legend",
    labels = c("#c0392b" = "Sangat Tinggi (≥2)",
               "#e67e22" = "Tinggi (1–2)",
               "#27ae60" = "Rendah (0.5–1)",
               "#2980b9" = "Sangat Rendah (<0.5)"),
    breaks = c("#c0392b","#e67e22","#27ae60","#2980b9")
  ) +
  labs(
    title    = "Standardized Incidence Ratio (SIR) Diare",
    subtitle = "Jawa Barat 2024 | Garis = SIR 1.0 (rata-rata provinsi)",
    x        = "SIR",
    y        = NULL,
    fill     = "Kategori Risiko"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title          = element_text(face = "bold", size = 12),
    panel.grid.major.y  = element_blank(),
    legend.position     = "bottom"
  )
print(p_sir)

Gambar 3. Diagram Batang Standardized Incidence Ratio (SIR) Diare

Interpretasi:

Berdasarkan Gambar 3, nilai Standardized Incidence Ratio (SIR) diare di Jawa Barat tahun 2024 menunjukkan variasi antar kabupaten/kota. Kota Cirebon memiliki nilai SIR tertinggi, diikuti oleh Kota Bogor dan Kota Sukabumi. Nilai SIR di atas 1 menunjukkan bahwa jumlah kasus diare yang terjadi lebih tinggi dibandingkan jumlah kasus yang diharapkan berdasarkan rata-rata provinsi. Sebaliknya, Kabupaten Bekasi memiliki nilai SIR terendah dengan nilai di bawah 1, yang menunjukkan bahwa jumlah kasus yang terjadi lebih rendah dibandingkan jumlah kasus yang diharapkan. Secara umum, sebagian besar kota memiliki nilai SIR yang lebih tinggi dibandingkan kabupaten, mengindikasikan risiko kejadian diare yang relatif lebih besar pada wilayah tersebut.

---

peta_sir <- ggplot(jabar) +
  
  geom_sf(
    aes(fill = kategori_SIR),
    color = "white",
    linewidth = 0.3
  ) +
  
  scale_fill_manual(
    values = c(
      "<0.80"     = "#4575B4",
      "0.80-0.99" = "#91BFDB",
      "1.00-1.19" = "#FFFFBF",
      "1.20-1.49" = "#FC8D59",
      ">=1.50"    = "#D73027"
    ),
    name = "Kategori SIR"
  ) +
  
  labs(
    title = "Peta Standardized Incidence Ratio (SIR)",
    subtitle = "Kabupaten/Kota Jawa Barat Tahun 2024",
    caption = "SIR = Observed Cases / Expected Cases"
  ) +
  
  theme_minimal() +
  
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(
      face = "bold",
      hjust = 0.5
    )
  )

print(peta_sir)

Gambar 4. Peta Choropleth Standardized Incidence Ratio (SIR) Diare

Interpretasi:

Berdasarkan Gambar 4, distribusi spasial nilai SIR diare di Jawa Barat tahun 2024 menunjukkan adanya perbedaan tingkat risiko antarwilayah. Wilayah dengan kategori SIR tinggi (≥1,50) terkonsentrasi pada beberapa kota, terutama Kota Cirebon, Kota Bogor, dan Kota Sukabumi. Nilai SIR yang tinggi menunjukkan bahwa jumlah kasus diare yang diamati jauh lebih besar dibandingkan jumlah kasus yang diharapkan. Sebaliknya, beberapa wilayah berada pada kategori SIR rendah (<0,80), yang menunjukkan jumlah kasus lebih rendah dibandingkan nilai yang diharapkan. Pola ini mengindikasikan adanya variasi risiko kejadian diare antar kabupaten/kota sehingga diperlukan analisis spasial lanjutan untuk mengidentifikasi kemungkinan pengelompokan wilayah berisiko tinggi.

---

3.4 Analisis Spasial (Morans’I dan LISA)

3.4.1 Global Moran’s I

Tabel 4. Hasil Uji Global Moran’s I

# Bangun matriks pembobot spasial (Queen contiguity)
  nb <- poly2nb(jabar, queen = TRUE)
  lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

  
moran_ir <- moran.test(jabar$IR, lw, zero.policy = TRUE,
                         na.action = na.omit)
hasil_moran <- data.frame(
  Statistik = c(
    "Moran's I",
    "Expected Moran's I",
    "Variance",
    "Z-score",
    "p-value"
  ),
  Nilai = c(
    moran_ir$estimate[1],
    moran_ir$estimate[2],
    moran_ir$estimate[3],
    moran_ir$statistic,
    moran_ir$p.value
  )
)

kable(
  hasil_moran,row.names = FALSE)

Statistik	Nilai
Moran’s I	0.0597445
Expected Moran’s I	-0.0384615
Variance	0.0170655
Z-score	0.7517591
p-value	0.2260980

Interpretasi:

Berdasarkan hasil uji Moran’s I pada Tabel 3, diperoleh nilai Moran’s I sebesar 0,0597 yang menunjukkan adanya kecenderungan autokorelasi spasial positif yang sangat lemah pada kejadian diare antar kabupaten/kota di Jawa Barat tahun 2024. Nilai tersebut lebih besar dibandingkan nilai ekspektasi Moran’s I sebesar -0,0385, namun nilai autokorelasi yang terbentuk relatif kecil. Hasil pengujian statistik menunjukkan nilai Z-score sebesar 0,7518 dengan p-value sebesar 0,2261 (p > 0,05). Dengan demikian, tidak terdapat bukti yang cukup untuk menolak hipotesis nol, sehingga dapat disimpulkan bahwa tidak terdapat autokorelasi spasial global yang signifikan pada distribusi kasus diare di Jawa Barat tahun 2024. Dengan kata lain, pola penyebaran kasus diare antar kabupaten/kota cenderung tidak membentuk pola pengelompokan (cluster) maupun penyebaran yang signifikan secara global.

---

  moran.plot(
    jabar$IR, lw,
    zero.policy = TRUE,
    labels      = as.character(jabar$nama_kab),
    main        = "Moran Scatterplot – IR Diare Jawa Barat 2024",
    xlab        = "IR (standardized)",
    ylab        = "Spatial Lag IR"
  )

Gambar 5. Moran’s Scatterplot

Interpretasi:

Berdasarkan Gambar 5, Moran Scatterplot menunjukkan kecenderungan autokorelasi spasial positif yang lemah pada nilai Incidence Rate (IR) diare di Jawa Barat tahun 2024. Hal ini ditunjukkan oleh kemiringan garis regresi yang positif, namun relatif landai. Beberapa wilayah seperti Kota Cirebon dan Kota Sukabumi berada pada kelompok wilayah dengan nilai IR tinggi, sedangkan Kabupaten Cirebon menunjukkan pola yang berbeda dibandingkan wilayah di sekitarnya. Secara umum, sebagian besar kabupaten/kota tersebar di sekitar nilai rata-rata, sehingga tidak terlihat adanya pengelompokan spasial yang kuat. Hasil ini sejalan dengan uji Moran’s I yang menunjukkan bahwa autokorelasi spasial global tidak signifikan.

---

3.4.2 Local Moran’s I (LISA)

  local_moran <- localmoran(jabar$IR, lw, zero.policy = TRUE)
  
  jabar$Ii      <- local_moran[, 1]
  jabar$Z_Ii    <- local_moran[, 4]
  jabar$p_lisa  <- local_moran[, 5]
  
  # Klasifikasi klaster LISA
  IR_std  <- scale(jabar$IR)[, 1]
  lag_IR  <- lag.listw(lw, IR_std, zero.policy = TRUE)
  
  jabar$LISA_cluster <- "Not Significant"
  jabar$LISA_cluster[IR_std > 0 & lag_IR > 0 & jabar$p_lisa < 0.05] <- "High-High"
  jabar$LISA_cluster[IR_std < 0 & lag_IR < 0 & jabar$p_lisa < 0.05] <- "Low-Low"
  jabar$LISA_cluster[IR_std > 0 & lag_IR < 0 & jabar$p_lisa < 0.05] <- "High-Low"
  jabar$LISA_cluster[IR_std < 0 & lag_IR > 0 & jabar$p_lisa < 0.05] <- "Low-High"
  
  # Peta LISA
  warna_lisa <- c(
    "High-High"       = "#e74c3c",
    "Low-Low"         = "#3498db",
    "High-Low"        = "#e67e22",
    "Low-High"        = "#85c1e9",
    "Not Significant" = "grey85"
  )
  
  p_lisa <- ggplot(jabar) +
    geom_sf(aes(fill = LISA_cluster), color = "white", linewidth = 0.3) +
    scale_fill_manual(
      values = warna_lisa,
      name   = "Klaster LISA"
    ) +
    labs(
      title    = "Peta Klaster LISA – IR Diare Jawa Barat 2024",
      subtitle = "Signifikansi p < 0.05 | Queen Contiguity Weights",
      caption  = "Sumber: Open Data Jabar 2024"
    ) +
    theme_minimal(base_size = 11) +
    theme(
      plot.title    = element_text(face = "bold", size = 13),
      legend.position = "right"
    )
  print(p_lisa)

Gambar 6. Peta Klaster LISA

Interpretasi:

Berdasarkan Gambar 6, hasil analisis Local Indicators of Spatial Association (LISA) menunjukkan bahwa tidak terdapat klaster spasial yang signifikan pada tingkat signifikansi 5% (p < 0,05) untuk nilai Incidence Rate (IR) diare di kabupaten/kota Provinsi Jawa Barat tahun 2024. Seluruh wilayah dikategorikan sebagai Not Significant, yang menunjukkan bahwa tidak ditemukan pola pengelompokan lokal berupa klaster High-High (HH), Low-Low (LL), maupun outlier spasial High-Low (HL) dan Low-High (LH). Hasil ini mengindikasikan bahwa distribusi kejadian diare di Jawa Barat tahun 2024 tidak menunjukkan adanya pengelompokan spasial lokal yang signifikan antar kabupaten/kota. Temuan tersebut sejalan dengan hasil uji Moran’s I yang juga menunjukkan tidak adanya autokorelasi spasial global yang signifikan.

---

3.5 UJi Overdispersi dan VIF

3.5.1 UJi Overdispersi

Tabel 5. Hasil Uji Overdispersi

# Variance-to-Mean Ratio (VMR)
vmr <- var(df$diare, na.rm = TRUE) /
       mean(df$diare, na.rm = TRUE)

tabel_vmr <- data.frame(
  Statistik = c(
    "Rata-rata Kasus",
    "Varians Kasus",
    "Variance-to-Mean Ratio (VMR)"
  ),
  Nilai = c(
    mean(df$diare, na.rm = TRUE),
    var(df$diare, na.rm = TRUE),
    vmr
  )
)

kable(tabel_vmr)

Statistik	Nilai
Rata-rata Kasus	24899.52
Varians Kasus	152153494.80
Variance-to-Mean Ratio (VMR)	6110.70

Interpretasi:

Hasil analisis menunjukkan terjadinya overdispersi pada model regresi Poisson, dengan nilai Variance-to-Mean Ratio (VMR) sebesar 6110,70 (Var>Mean). Fenomena ini signifikan secara statistik berdasarkan Z Statistic sebesar 3,816 dengan p-value 0,0000679 (p<0,05), sehingga hipotesis equidispersion ditolak. Pelanggaran asumsi ini menyebabkan regresi Poisson standar menjadi tidak efisien karena menghasilkan standard error yang terlalu kecil (underestimated) dan memicu signifikansi semu. Oleh karena itu, untuk mengatasi overdispersi dan menjamin estimasi parameter yang valid, analisis dialihkan menggunakan model Regresi Binomial Negatif.

---

3.5.2 Multikolinearitas

Tabel 6. Hasil Uji Multikolinearitas

m_vif <- lm(
  log1p(diare) ~
    faskes_std + air_std + sanitasi_std +
    ipm_std + kemiskinan_std + kepadatan_std,
  data = df
)

vif_vals <- vif(m_vif)

tabel_vif <- data.frame(
  Variabel = names(vif_vals),
  VIF = as.numeric(vif_vals)
)

kable(
  tabel_vif,
  digits = 3)

Variabel	VIF
faskes_std	1.220
air_std	2.036
sanitasi_std	1.235
ipm_std	7.542
kemiskinan_std	3.186
kepadatan_std	7.246

Interpretasi:

Evaluasi terhadap asumsi non-multikolinieritas dilakukan dengan mengamati nilai Variance Inflation Factor (VIF) pada masing-masing variabel prediktor. Secara teoretis, gejala multikolinieritas dianggap tidak menjadi masalah yang mengganggu estimasi model apabila nilai VIF lebih kecil dari 10 (atau dalam kriteria yang lebih ketat, lebih kecil dari 5). Hasil analisis menunjukkan bahwa variabel faskes (1,220), air minum (2,036), sanitasi (1,235), dan kemiskinan (3,186) memiliki nilai VIF yang sangat aman di bawah 5. Di sisi lain, variabel ipm_std dan kepadatan menunjukkan nilai VIF yang relatif lebih tinggi, masing-masing sebesar 7,542 dan 7,246. Tingginya nilai VIF pada kedua variabel ini mengindikasikan adanya korelasi moderat yang wajar, mengingat indikator pembangunan manusia (IPM) dan kepadatan penduduk secara empiris sering kali saling berkaitan dalam studi sosio-ekonomi. Karena tidak ada variabel yang memiliki nilai VIF melebihi batas kritis 10 (VIF < 10), maka dapat disimpulkan bahwa model regresi ini bebas dari masalah multikolinieritas yang serius.

---

3.6 Regresi Negative Binomial

Tabel 7. Perbandingan Model Negative Binomial

# Model 1: tanpa kemiskinan & kepadatan
m1 <- glm.nb(
  diare ~
    faskes_std + air_std + sanitasi_std + ipm_std +
    offset(ln_penduduk),
  data = df
)

# Model 2: tambah kemiskinan
m2 <- glm.nb(
  diare ~
    faskes_std + air_std + sanitasi_std + ipm_std +
    kemiskinan_std +
    offset(ln_penduduk),
  data = df
)

# Model 3: lengkap (semua prediktor)
m3 <- glm.nb(
  diare ~
    faskes_std + air_std + sanitasi_std + ipm_std +
    kemiskinan_std + kepadatan_std +
    offset(ln_penduduk),
  data = df
)

tabel_model <- data.frame(
  Model = c(
    "Model 1 (4 prediktor)",
    "Model 2 (5 prediktor)",
    "Model 3 (6 prediktor)"
  ),
  `Jumlah Parameter` = c(
    length(coef(m1)) + 1,
    length(coef(m2)) + 1,
    length(coef(m3)) + 1
  ),
  LogLik = c(
    as.numeric(logLik(m1)),
    as.numeric(logLik(m2)),
    as.numeric(logLik(m3))
  ),
  AIC = c(
    AIC(m1),
    AIC(m2),
    AIC(m3)
  ),
  BIC = c(
    BIC(m1),
    BIC(m2),
    BIC(m3)
  ),
  Theta = c(
    m1$theta,
    m2$theta,
    m3$theta
  )
)

kable(
  tabel_model,
  digits = 4,
  caption = "Perbandingan Model Regresi Binomial Negatif"
)

Perbandingan Model Regresi Binomial Negatif

Model	Jumlah.Parameter	LogLik	AIC	BIC	Theta
Model 1 (4 prediktor)	6	-277.7997	567.5995	575.3745	9.6173
Model 2 (5 prediktor)	7	-277.7179	569.4357	578.5066	9.6741
Model 3 (6 prediktor)	8	-274.4831	564.9662	575.3329	12.2077

Interpretasi:

Berdasarkan Tabel 7, pemilihan model terbaik dari ketiga alternatif model Negative Binomial dilakukan dengan membandingkan nilai Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), dan Log-Likelihood. Model yang ideal adalah model yang memiliki nilai AIC dan BIC terendah, serta nilai Log-Likelihood tertinggi (paling mendekati nol). Hasil estimasi menunjukkan bahwa Model 3 (dengan 8 parameter) memiliki nilai AIC terendah sebesar 564,97 dan nilai BIC terendah sebesar 575,33. Konsisten dengan hal tersebut, Model 3 juga menghasilkan nilai Log-Likelihood tertinggi, yaitu sebesar -274,48. Selain itu, nilai parameter dispersi (Theta) pada Model 3 adalah yang terbesar, yakni 12,208, yang menunjukkan kemampuan model yang lebih baik dalam mengatasi overdispersion pada data. Oleh karena itu, Model 3 dipilih sebagai model terbaik dan paling fit untuk analisis data selanjutnya. Berikut adalah tabel hasil regresi negative binomial pada model 3.

---

Tabel 8. Hasil Regresi Negative Binomial Model 3

fitted_vals   <- fitted(m3)
cor_obs_fit   <- cor(df$diare, fitted_vals)
sm3   <- summary(m3)
ci3   <- suppressMessages(confint(m3))
coef3 <- coef(m3)
pval3 <- sm3$coefficients[,4]

var_labels <- c(
  "(Intercept)"    = "(Intercept)",
  "faskes_std"     = "Log Faskes (std)",
  "air_std"        = "Air Bersih (std)",
  "sanitasi_std"   = "Sanitasi (std)",
  "ipm_std"        = "IPM (std)",
  "kemiskinan_std" = "Log Garis Kemiskinan (std)",
  "kepadatan_std"  = "Log Kepadatan (std)"
)

tabel_irr <- data.frame(
  Variabel = var_labels[names(coef3)],
  Beta     = coef3,
  SE       = sm3$coefficients[,2],
  z        = sm3$coefficients[,3],
  p_value  = pval3,
  IRR      = exp(coef3),
  CI_bawah = exp(ci3[,1]),   # ← ganti nama
  CI_atas  = exp(ci3[,2]),   # ← ganti nama
  Sig      = symnum(
    pval3,
    cutpoints = c(0,0.001,0.01,0.05,0.1,1),
    symbols   = c("***","**","*","."," "),
    legend    = FALSE
  )
)

kable(
  tabel_irr,
  digits = 4, row.names = FALSE)

Variabel	Beta	SE	z	p_value	IRR	CI_bawah	CI_atas	Sig
(Intercept)	-4.1571	0.0551	-75.4489	0.0000	0.0157	0.0141	0.0175	***
Log Faskes (std)	-0.3234	0.0620	-5.2136	0.0000	0.7237	0.6361	0.8210	***
Air Bersih (std)	0.0112	0.0801	0.1395	0.8890	1.0112	0.8623	1.1843
Sanitasi (std)	-0.0113	0.0624	-0.1806	0.8567	0.9888	0.8735	1.1107
IPM (std)	-0.3571	0.1542	-2.3160	0.0206	0.6997	0.5238	0.9383	*
Log Garis Kemiskinan (std)	-0.1174	0.1002	-1.1711	0.2415	0.8893	0.7270	1.0888
Log Kepadatan (std)	0.3993	0.1511	2.6420	0.0082	1.4908	1.1156	1.9943	**

Interpretasi:

Berdasarkan hasil tabel diatas diperoleh tiga prediktor menunjukkan pengaruh signifikan terhadap jumlah kasus diare:

1. Log Fasilitas Kesehatan: IRR = 0,724 (95% CI: 0,636–0,821; p < 0,001). Setiap kenaikan 1 SD log-faskes menurunkan expected kasus sebesar 27,6%. Ini merupakan faktor protektif terkuat dan paling signifikan.

2. IPM : IRR = 0,700 (95% CI: 0,524–0,938; p = 0,021). Setiap kenaikan 1 SD IPM menurunkan expected kasus sebesar 30,0%. IPM mencerminkan dimensi kesehatan, pendidikan, dan daya beli yang secara holistik melindungi dari penyakit berbasis lingkungan.

3. Log Kepadatan Penduduk : IRR = 1,491 (95% CI: 1,116–1,994; p = 0,008). Setiap kenaikan 1 SD log-kepadatan meningkatkan expected kasus sebesar 49,1%. Kepadatan tinggi mempercepat transmisi fecal-oral.

Tiga prediktor lain akses air bersih (p = 0,889), sanitasi (p = 0,857), dan garis kemiskinan (p = 0,242) tidak menunjukkan pengaruh signifikan secara statistik. Goodness-of-fit test menghasilkan p = 0,168 > 0,05, mengkonfirmasi Model 3 memiliki fit yang baik.

---

# Forest plot - ganti nama kolom di aes()
irr_plot <- tabel_irr %>%
  filter(Variabel != "(Intercept)") %>%
  mutate(
    Variabel = factor(Variabel, levels = rev(Variabel)),
    sig_col  = ifelse(trimws(Sig) %in% c("***","**","*"), "Signifikan", "Tidak Signifikan")
  )

p_irr <- ggplot(irr_plot, aes(x = IRR, y = Variabel, color = sig_col)) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "gray40", linewidth = 0.9) +
  geom_errorbarh(aes(xmin = CI_bawah, xmax = CI_atas),   # ← nama baru
                 height = 0.25, linewidth = 0.9) +
  geom_point(size = 4.5) +
  geom_text(aes(label = paste0("IRR=", round(IRR,3), " ", trimws(Sig))),
            hjust = -0.15, size = 3, color = "black") +
  scale_color_manual(
    values = c("Signifikan" = "#e74c3c", "Tidak Signifikan" = "gray55"),
    name   = NULL
  ) +
  scale_x_continuous(
    limits = c(
      min(irr_plot$CI_bawah) * 0.85,   # ← nama baru
      max(irr_plot$CI_atas)  * 1.2     # ← nama baru
    )
  ) +
  labs(
    title    = "Forest Plot Incidence Rate Ratio (IRR) – Model Negative Binomial 3",
    subtitle = "Titik merah = signifikan (p < 0.05) | Garis = 95% CI",
    x        = "IRR (Incidence Rate Ratio)",
    y        = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(
    plot.title         = element_text(face = "bold", size = 12),
    panel.grid.major.y = element_blank(),
    legend.position    = "bottom"
  )
print(p_irr)

Gambar 7. Forest Plot IRR

Interpretasi:

Hasil estimasi Model Regresi Binomial Negatif 3 menunjukkan pengaruh variabel prediktor terhadap Incidence Rate Ratio (IRR) kasus diare. Berdasarkan plot, terdapat tiga variabel yang berpengaruh signifikan secara statistik terhadap incidence rate kasus diare (p < 0,05), yaitu Log Faskes, IPM, dan Log Kepadatan Penduduk. Sementara itu, variabel Air Bersih, Sanitasi, dan Log Garis Kemiskinan tidak menunjukkan pengaruh yang signifikan (p >= 0,05) karena interval kepercayaan 95% melewati nilai 1.Variabel Log Faskes memiliki nilai IRR sebesar 0,724, yang berarti setiap kenaikan satu satuan standar Log Faskes diasosiasikan dengan penurunan incidence rate kasus diare sebesar 27,6% (1 - 0,724), dengan asumsi variabel lain konstan. Variabel IPM juga memiliki pengaruh protektif yang signifikan dengan nilai IRR sebesar 0,71, menunjukkan bahwa peningkatan satu satuan standar IPM berkaitan dengan penurunan incidence rate kasus diare sebesar 29%. Sebaliknya, Log Kepadatan Penduduk bertindak sebagai faktor risiko yang signifikan dengan nilai IRR sebesar 1,491. Hal ini mengindikasikan bahwa setiap kenaikan satu satuan standar Log Kepadatan Penduduk meningkatkan incidence rate kasus diare sebesar 49,1%.

---

Tabel 9. Hasil Diagnostik Model 3

tabel_fit <- data.frame(
  Indikator = c(
    "Korelasi Observed-Fitted",
    "Persentase Variasi Terjelaskan"
  ),
  Nilai = c(
    cor_obs_fit,
    cor_obs_fit^2 * 100
  )
)

kable(
  tabel_fit,
  digits = 2,
  caption = "Kemampuan Prediksi Model"
)

Kemampuan Prediksi Model

Indikator	Nilai
Korelasi Observed-Fitted	0.86
Persentase Variasi Terjelaskan	73.74

Interpretasi:

Hasil evaluasi kemampuan prediksi model menunjukkan kinerja pemodelan yang sangat baik. Nilai Korelasi Observed-Fitted sebesar 0,86 mengindikasikan adanya hubungan linier yang kuat antara nilai aktual yang diamati di lapangan dengan nilai yang diprediksi oleh model. Selain itu, indikator Persentase Variasi Terjelaskan menunjukkan nilai sebesar 73,74%, yang berarti bahwa variabel-variabel prediktor di dalam model mampu menjelaskan sekitar 73,74% total variasi dari variabel dependen kasus diare. Sementara itu, sisanya sebesar 26,26% dijelaskan oleh faktor-faktor lain di luar model. Kombinasi kedua nilai ini menegaskan bahwa model yang dibangun memiliki tingkat kebaikan suai (goodness of fit) serta daya prediksi yang tinggi dan valid untuk digunakan.

---

4 Pembahasan

4.1 Pola Distribusi Risiko Diare

Berdasarkan hasil analisis data, beban penyakit diare menunjukkan variasi spasial yang signifikan di antara kabupaten dan kota di Provinsi Jawa Barat. Kota Cirebon tercatat sebagai wilayah dengan nilai Incidence Rate (IR) tertinggi, yakni mencapai 442,83 kasus per 10.000 penduduk. Wilayah tersebut juga memiliki nilai Standardized Incidence Ratio (SIR) sebesar 3,32, yang mengindikasikan bahwa jumlah kasus riil di lapangan meningkat lebih dari tiga kali lipat dibandingkan dengan estimasi yang diharapkan berdasarkan rata-rata provinsi. Tren risiko SIR yang relatif tinggi ini juga teridentifikasi di Kota Bogor dan Kota Sukabumi.

Besarnya tingkat risiko diare di kawasan perkotaan tersebut kemungkinan dipengaruhi oleh tingginya kepadatan penduduk, yang mempercepat transmisi patogen melalui rute fekal-oral. Selain faktor demografi, ketersediaan infrastruktur pelayanan kesehatan serta mekanisme pelaporan data yang lebih optimal di area perkotaan turut berkontribusi terhadap tingginya cakupan deteksi dan pencatatan kasus diare jika dibandingkan dengan wilayah lainnya.

Sebaliknya, beberapa kawasan metropolitan seperti Kota Bandung dan Kota Bekasi justru menunjukkan nilai IR dan SIR yang berada di bawah rata-rata capaian provinsi. Kondisi ini disinyalir berkaitan dengan kualitas sanitasi lingkungan, kemudahan akses air bersih, kualitas pelayanan medis, serta indeks pembangunan manusia yang relatif lebih maju. Sementara itu, Kabupaten Bekasi mencatatkan nilai IR dan SIR paling rendah dalam penelitian ini, yang mencerminkan bahwa akumulasi kasus diare aktual di wilayah tersebut berada di bawah angka yang diproyeksikan.

4.2 Pola Spasial

Hasil analisis spasial menunjukkan bahwa tidak terdapat autokorelasi spasial global yang signifikan pada kejadian diare di Jawa Barat tahun 2024, sebagaimana ditunjukkan oleh nilai Moran’s I sebesar 0,060 dengan nilai p sebesar 0,226. Selain itu, analisis Local Indicators of Spatial Association (LISA) juga tidak menemukan klaster spasial yang signifikan pada tingkat signifikansi 5%.

Temuan ini menunjukkan bahwa distribusi risiko diare antar kabupaten/kota tidak membentuk pola pengelompokan geografis yang jelas. Dengan demikian, variasi kejadian diare di Jawa Barat kemungkinan lebih dipengaruhi oleh karakteristik masing-masing wilayah, seperti kondisi sosial ekonomi, perilaku kesehatan masyarakat, kualitas surveilans, serta akses terhadap pelayanan kesehatan, dibandingkan oleh faktor kedekatan geografis antarwilayah. Oleh karena itu, upaya pengendalian diare perlu difokuskan pada karakteristik risiko spesifik di setiap wilayah, bukan semata-mata berdasarkan kedekatan lokasi geografis.

4.3 Faktor yang Berpengaruh

Hasil analisis regresi Binomial Negatif menunjukkan bahwa beberapa faktor lingkungan dan sosial ekonomi memiliki pengaruh signifikan terhadap jumlah kasus diare. Ketersediaan fasilitas kesehatan terbukti sebagai faktor protektif yang signifikan terhadap kejadian diare (IRR = 0,724; p < 0,001). Temuan ini mengindikasikan bahwa peningkatan ketersediaan fasilitas kesehatan berpotensi menurunkan jumlah kasus diare melalui peningkatan deteksi dini, perbaikan penatalaksanaan kasus, serta perluasan akses masyarakat terhadap pelayanan kesehatan dasar.

Selain itu, Indeks Pembangunan Manusia (IPM) juga menunjukkan pengaruh protektif yang signifikan (IRR = 0,700; p = 0,021). Hal ini menunjukkan bahwa wilayah dengan tingkat pembangunan manusia yang lebih tinggi cenderung memiliki beban diare yang lebih rendah. Kondisi tersebut dapat dijelaskan oleh meningkatnya tingkat pendidikan, status kesehatan, dan kesejahteraan masyarakat yang mendukung penerapan perilaku hidup bersih dan sehat serta peningkatan akses terhadap sanitasi yang layak.

Sebaliknya, kepadatan penduduk ditemukan sebagai faktor risiko yang signifikan terhadap peningkatan kasus diare (IRR = 1,491; p = 0,008). Temuan ini mengindikasikan bahwa wilayah dengan tingkat kepadatan penduduk yang lebih tinggi memiliki risiko penularan penyakit yang lebih besar akibat meningkatnya interaksi antarindividu serta tekanan terhadap ketersediaan infrastruktur sanitasi dan lingkungan.

Sementara itu, variabel akses air bersih dan sanitasi tidak menunjukkan hubungan yang signifikan dalam model. Kondisi ini diduga disebabkan oleh relatif homogennya cakupan akses antarwilayah sehingga variasi data menjadi terbatas. Selain itu, indikator yang digunakan lebih merepresentasikan aspek ketersediaan layanan, bukan kualitas atau keberlanjutan akses yang diterima oleh masyarakat.

4.4 Justifikasi Model Negative Binomial

Pemilihan model regresi Binomial Negatif dalam penelitian ini didasarkan pada adanya indikasi kuat overdispersi pada data kasus diare. Hal ini ditunjukkan oleh nilai Variance-to-Mean Ratio (VMR) yang jauh lebih besar dari satu, serta hasil uji dispersi Cameron dan Trivedi yang signifikan (p < 0,001). Selain itu, nilai devians pada model Poisson juga menunjukkan ketidaksesuaian model terhadap struktur data.

Kondisi tersebut menunjukkan bahwa asumsi kesetaraan antara rata-rata dan varians pada model Poisson tidak terpenuhi. Oleh karena itu, penggunaan model Binomial Negatif menjadi lebih tepat karena mampu mengakomodasi varians yang lebih besar dari rata-rata (overdispersi), sehingga menghasilkan estimasi parameter dan standard error yang lebih reliabel. Dengan demikian, model ini memberikan inferensi statistik yang lebih valid dalam menganalisis determinan kasus diare.

4.5 Keterbatasan Penelitian

Penelitian ini memiliki beberapa keterbatasan yang perlu diperhatikan dalam menginterpretasikan hasil.

1. Desain penelitian yang digunakan adalah potong lintang (cross-sectional), sehingga tidak memungkinkan penarikan kesimpulan kausal antara variabel determinan dan kejadian diare.
2. Analisis dilakukan pada tingkat agregat wilayah, sehingga terdapat potensi terjadinya ecological fallacy apabila hasil penelitian diinterpretasikan pada tingkat individu.
3. Kemungkinan adanya under-reporting kasus diare masih dapat terjadi, mengingat tidak seluruh kasus tercatat secara lengkap dalam sistem pelaporan kesehatan.
4. Jumlah unit analisis yang relatif terbatas dapat memengaruhi kekuatan statistik dalam mendeteksi hubungan antarvariabel, sehingga beberapa asosiasi mungkin tidak teridentifikasi secara optimal.
5. Ketiadaan data temporal yang lebih rinci, seperti data bulanan atau musiman, membatasi analisis dalam mengevaluasi variasi spasio-temporal kejadian diare sepanjang tahun.

---

5 Kesimpulan dan Saran

5.1 Kesimpulan

Berdasarkan analisis epidemiologi kasus diare di Provinsi Jawa Barat tahun 2024 pada 27 kabupaten/kota, diperoleh beberapa kesimpulan utama. Secara umum, angka insidens (IR) diare di Jawa Barat sebesar 133,54 per 10.000 penduduk dengan total 672.287 kasus. Terdapat disparitas yang cukup besar antarwilayah, di mana Kota Cirebon menunjukkan IR tertinggi sebesar 442,83 per 10.000 penduduk dengan SIR sebesar 3,32, diikuti oleh Kota Bogor (321,38; SIR = 2,41) dan Kota Sukabumi (311,29; SIR = 2,33). Secara keseluruhan, sebanyak 17 dari 27 wilayah (63%) memiliki beban kasus yang lebih tinggi dibandingkan rata-rata provinsi (SIR ≥ 1,0).

1. Analisis spasial menunjukkan bahwa tidak terdapat autokorelasi spasial yang signifikan (Moran’s I = 0,060; p = 0,226). Temuan ini diperkuat oleh tidak ditemukannya klaster signifikan pada analisis Local Indicators of Spatial Association (LISA). Dengan demikian, distribusi angka insidens diare di Jawa Barat dapat dikatakan bersifat acak secara geografis, sehingga pendekatan model regresi non-spasial dianggap memadai untuk analisis inferensial.
2. Selanjutnya, indikasi overdispersi pada data terbukti secara kuat melalui nilai Variance-to-Mean Ratio (VMR = 6.110,70), dispersion estimate (1.656,22), serta hasil uji yang signifikan (p < 0,001). Kondisi ini menegaskan ketidaksesuaian asumsi model Poisson dan memberikan justifikasi penggunaan regresi Binomial Negatif dalam analisis.
3. Hasil pemodelan menunjukkan bahwa Model Binomial Negatif terbaik (Model 3; AIC = 564,97) mengidentifikasi tiga variabel yang berpengaruh signifikan. Ketersediaan fasilitas kesehatan merupakan faktor protektif yang paling kuat (IRR = 0,724; p < 0,001), diikuti oleh Indeks Pembangunan Manusia (IPM) sebagai faktor protektif (IRR = 0,700; p = 0,021). Sebaliknya, kepadatan penduduk merupakan faktor risiko yang signifikan terhadap peningkatan kasus diare (IRR = 1,491; p = 0,008).
4. Sementara itu, variabel akses sanitasi, akses air bersih, dan garis kemiskinan tidak menunjukkan pengaruh yang signifikan dalam model multivariat. Hal ini diduga berkaitan dengan adanya efek ambang (threshold effect), keterbatasan variasi data antarwilayah, serta kemungkinan adanya multikolinearitas antarvariabel sosial ekonomi.

5.2 Saran

Berdasarkan temuan penelitian, beberapa rekomendasi yang dapat diajukan adalah sebagai berikut.

1. Diperlukan prioritas intervensi komprehensif pada wilayah dengan beban diare tinggi, khususnya Kota Cirebon, Kota Bogor, dan Kota Sukabumi yang memiliki SIR ≥ 2,0. Intervensi tersebut mencakup audit menyeluruh terhadap kualitas air minum, kondisi infrastruktur sanitasi, serta sistem pengelolaan limbah untuk mengidentifikasi sumber risiko utama secara lebih spesifik.
2. Penguatan kapasitas dan ketersediaan fasilitas kesehatan perlu menjadi fokus utama pada wilayah dengan rasio fasilitas kesehatan rendah dan insidens tinggi. Sebagai contoh, Kota Cirebon dengan jumlah fasilitas kesehatan yang relatif terbatas dibandingkan jumlah penduduk memerlukan peningkatan akses layanan kesehatan dasar sebagai strategi intervensi yang relatif efisien namun berdampak signifikan.
3. Program peningkatan Indeks Pembangunan Manusia (IPM) perlu diintegrasikan dengan strategi pengendalian diare, terutama melalui peningkatan pendidikan kesehatan masyarakat, literasi perilaku hidup bersih dan sehat, serta peningkatan daya beli rumah tangga. Hal ini penting mengingat IPM terbukti sebagai faktor protektif yang bersifat multidimensional.
4. Penguatan sistem surveilans epidemiologi perlu dilakukan pada wilayah dengan beban diare tinggi untuk meningkatkan kelengkapan (completeness) dan akurasi pelaporan kasus. Hal ini penting guna memastikan bahwa gambaran epidemiologi yang dihasilkan lebih representatif terhadap kondisi sebenarnya.
5. Penelitian lanjutan disarankan untuk menggunakan desain data panel multi-tahun (2018–2024) guna menganalisis tren temporal serta mengidentifikasi determinan secara lebih presisi. Pendekatan spasiotemporal atau model data panel juga direkomendasikan untuk memperkuat inferensi kausal dan memahami dinamika penyebaran penyakit secara lebih komprehensif.

---

Daftar Pustaka

1. Cameron, A. C., & Trivedi, P. K. (1990). Regression-based tests for overdispersion in the Poisson model. Journal of Econometrics, 46(3), 347–364.
2. Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). SAGE Publications.
3. Hilbe, J. M. (2011). Negative Binomial Regression (2nd ed.). Cambridge University Press.
4. Kementerian Kesehatan Republik Indonesia. (2024). Profil Kesehatan Indonesia 2024. Kemenkes RI.
5. Open Data Jawa Barat. (2024). Jumlah Kasus Penyakit Diare berdasarkan Kabupaten/Kota di Jawa Barat. Diakses Juni 2026 dari https://opendata.jabarprov.go.id
6. Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
7. Bivand, R., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer.
8. BPS Provinsi Jawa Barat. (2024). Jawa Barat dalam Angka 2024. Badan Pusat Statistik Provinsi Jawa Barat.
9. Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37(1–2), 17–23.
10. R Core Team. (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.R-project.org/
11. Venables, W. N., & Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). Springer. [Paket MASS]
12. Wickham, H., et al. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686.
13. World Health Organization. (2023). Diarrhoeal Disease. Fact Sheet. https://www.who.int/news-room/fact-sheets/detail/diarrhoeal-disease

---