Analisis Epidemiologi Insidensi Diare di Indonesia Tahun 2018 Menggunakan Pendekatan Regresi OLS dan Autokorelasi Spasial

1 LINK YOUTUBE

Link : https://youtu.be/jSW2G6Ywxe8

2 DASHBOARD

Link : https://ghozzypurwana.shinyapps.io/DashboardEpidemGhozzy/

3 ABSTRAK

Diare merupakan salah satu penyakit berbasis lingkungan yang masih menjadi permasalahan kesehatan masyarakat di Indonesia dengan tingkat kejadian yang bervariasi antarprovinsi. Penelitian ini bertujuan untuk menganalisis faktor-faktor yang berhubungan dengan insidensi diare serta mengidentifikasi pola spasial penyebarannya di Indonesia pada tahun 2018. Data yang digunakan meliputi jumlah kasus diare sebagai variabel respon, serta jumlah penduduk, jumlah tenaga gizi, dan persentase rumah tangga dengan akses air minum layak sebagai variabel penjelas, yang bersumber dari Badan Pusat Statistik. Metode analisis yang digunakan meliputi analisis deskriptif spasial, regresi Ordinary Least Squares (OLS), korelasi Spearman, dan uji autokorelasi spasial Global Moran’s I. Hasil analisis OLS menunjukkan bahwa secara simultan variabel penjelas belum berpengaruh signifikan terhadap insidensi diare, dengan nilai R-squared sebesar 0,17. Secara parsial, jumlah tenaga gizi berpengaruh positif dan signifikan, sementara kepadatan penduduk menunjukkan pengaruh negatif yang lemah, dan akses air minum layak tidak berpengaruh signifikan. Uji Moran’s I mengindikasikan tidak adanya autokorelasi spasial global yang signifikan pada insidensi diare. Temuan ini menunjukkan bahwa variasi insidensi diare antarprovinsi tidak membentuk klaster spasial yang kuat dan kemungkinan dipengaruhi oleh faktor lain di luar model, seperti perilaku higienis dan sanitasi lingkungan. Penelitian ini diharapkan dapat menjadi dasar awal bagi pengembangan analisis epidemiologi spasial yang lebih komprehensif dalam perencanaan kebijakan kesehatan masyarakat.

4 PENDAHULUAN

4.1 Latar Belakang

Penyakit diare merupakan salah satu masalah kesehatan masyarakat yang masih sering terjadi di Indonesia dan memiliki dampak yang signifikan terhadap kualitas hidup masyarakat. Diare tidak hanya menyebabkan peningkatan angka kesakitan, tetapi juga dapat menurunkan produktivitas penduduk serta membebani sistem pelayanan kesehatan. Berdasarkan data Badan Pusat Statistik (BPS), jumlah kasus diare pada tahun 2018 menunjukkan variasi yang cukup besar antarprovinsi di Indonesia, yang mengindikasikan adanya perbedaan kondisi sosial, demografi, dan lingkungan antarwilayah [1].

Dalam kajian epidemiologi, kejadian penyakit diare dipengaruhi oleh berbagai faktor risiko. Jumlah penduduk menjadi salah satu faktor penting karena wilayah dengan populasi yang besar cenderung memiliki tingkat kepadatan yang tinggi, sehingga meningkatkan risiko penularan penyakit berbasis lingkungan seperti diare. Data BPS tahun 2018 menunjukkan bahwa distribusi jumlah penduduk antarprovinsi di Indonesia tidak merata, sehingga berpotensi memengaruhi perbedaan jumlah kasus diare antarwilayah [2].

Selain faktor kependudukan, ketersediaan sumber daya kesehatan juga berperan dalam pengendalian penyakit diare. Tenaga gizi memiliki peran penting dalam upaya pencegahan dan penanganan diare melalui edukasi gizi, perbaikan pola makan, serta penanganan gizi pada kelompok rentan. Namun, jumlah tenaga gizi antarprovinsi pada tahun 2018 masih menunjukkan ketimpangan, yang dapat memengaruhi efektivitas pengendalian penyakit diare [3].

Faktor lingkungan, khususnya akses terhadap sumber air minum yang layak, merupakan determinan utama dalam kejadian diare. Rumah tangga yang tidak memiliki akses air minum layak lebih berisiko terpapar mikroorganisme penyebab diare. Data BPS tahun 2018 menunjukkan bahwa persentase rumah tangga dengan akses air minum layak juga bervariasi antarprovinsi, sehingga berpotensi memengaruhi pola kejadian diare secara regional [4]. Oleh karena itu, pendekatan epidemiologi diperlukan untuk memahami hubungan antara faktor-faktor tersebut dengan kejadian diare secara lebih komprehensif.

4.2 Rumusan Masalah

Berdasarkan latar belakang tersebut, permasalahan yang dapat diidentifikasi dalam penelitian ini adalah sebagai berikut:

1. Terdapat perbedaan jumlah kasus diare antarprovinsi di Indonesia pada tahun 2018.

2. Jumlah penduduk, jumlah tenaga gizi, dan akses air minum layak diduga memengaruhi variasi kasus diare antarprovinsi.

3. Hubungan antara faktor demografi, sumber daya kesehatan, dan lingkungan dengan kejadian diare belum dianalisis secara epidemiologis berdasarkan data antarprovinsi.

4.3 Tujuan

Penelitian ini bertujuan untuk menganalisis faktor-faktor yang memengaruhi jumlah kasus diare di Indonesia pada tahun 2018 menggunakan pendekatan epidemiologi. Secara khusus, penelitian ini bertujuan untuk mengetahui pengaruh jumlah penduduk, jumlah tenaga gizi, dan persentase rumah tangga dengan akses air minum layak terhadap jumlah kasus diare antarprovinsi.

4.4 Batasan Penelitian

• Penelitian menggunakan data sekunder yang bersumber dari Badan Pusat Statistik (BPS).
• Unit analisis adalah provinsi di Indonesia, Tahun pengamatan dibatasi pada tahun 2018.
• Variabel yang dianalisis terdiri dari jumlah kasus diare sebagai variabel respon, serta jumlah penduduk, jumlah tenaga gizi, dan persentase rumah tangga dengan akses air minum layak sebagai variabel penjelas.
• Pendekatan yang digunakan adalah analisis epidemiologi, tanpa membahas aspek klinis penyakit diare secara mendalam.

5 TINJAUAN PUSTAKA

5.1 Agent, Host, dan Environment

Dalam kajian epidemiologi, kejadian penyakit dianalisis melalui konsep segitiga epidemiologi yang terdiri atas agent, host, dan environment. Ketiga komponen ini saling berinteraksi dan menentukan pola kejadian penyakit di suatu wilayah.

Agent merupakan penyebab langsung terjadinya penyakit. Pada penelitian ini, agent yang dikaji adalah penyebab penyakit diare, yang umumnya berasal dari mikroorganisme patogen seperti bakteri, virus, dan parasit yang masuk ke tubuh manusia melalui air atau makanan yang terkontaminasi. Data jumlah kasus diare yang digunakan dalam penelitian ini diperoleh dari Badan Pusat Statistik (BPS) berdasarkan jumlah kasus penyakit menurut provinsi tahun 2018 [1]. Tingginya kasus diare mencerminkan masih adanya paparan agent penyakit yang belum sepenuhnya terkendali.

Host adalah individu atau populasi yang berpotensi terinfeksi penyakit. Dalam penelitian ini, karakteristik host direpresentasikan oleh jumlah penduduk dan jumlah tenaga gizi. Jumlah penduduk menggambarkan besarnya populasi yang berisiko terpapar diare, di mana wilayah dengan jumlah penduduk yang lebih besar cenderung memiliki risiko penyebaran penyakit yang lebih tinggi [2]. Sementara itu, jumlah tenaga gizi mencerminkan kapasitas sumber daya kesehatan dalam upaya pencegahan dan penanganan masalah gizi serta penyakit berbasis lingkungan, termasuk diare [3]. Kedua variabel tersebut diperoleh dari data resmi BPS tahun 2018.

Environment mencakup kondisi lingkungan fisik yang memengaruhi keberlangsungan agent dan kerentanan host. Dalam penelitian ini, faktor lingkungan direpresentasikan oleh persentase rumah tangga yang memiliki akses terhadap air minum layak. Akses air bersih merupakan faktor kunci dalam pencegahan penyakit diare karena berperan dalam menjaga kebersihan makanan, minuman, dan sanitasi rumah tangga. Data mengenai akses air minum layak diperoleh dari BPS dan digunakan untuk menggambarkan kualitas lingkungan antarprovinsi di Indonesia [4]. Lingkungan dengan akses air bersih yang rendah berpotensi meningkatkan risiko penularan diare secara signifikan.

5.2 Ukuran Asosiasi

Ukuran asosiasi digunakan untuk menilai hubungan antara jumlah kasus diare dengan jumlah penduduk, jumlah tenaga gizi, dan akses air minum layak. Analisis dilakukan menggunakan OLS, korelasi Spearman, serta asosiasi spasial.

Regresi OLS digunakan untuk melihat hubungan rata-rata antarvariabel dan memberikan gambaran awal pengaruh faktor demografi, kesehatan, dan lingkungan terhadap kejadian diare antarprovinsi. Korelasi Spearman digunakan untuk mengukur kekuatan dan arah hubungan monoton antarvariabel tanpa asumsi normalitas data.

Untuk menangkap ketergantungan spasial, digunakan Moran’s I Global guna mengidentifikasi adanya pola pengelompokan kasus diare antarwilayah. Selanjutnya, LISA (Local Moran’s I) digunakan untuk mengidentifikasi klaster lokal dan outlier spasial seperti wilayah hot spot dan cold spot kejadian diare.

5.3 Desain Studi

Penelitian ini menggunakan desain studi ekologi dengan unit analisis berupa provinsi di Indonesia tahun 2018. Seluruh variabel dianalisis secara cross-sectional menggunakan data sekunder berskala wilayah.

Analisis mengombinasikan metode statistik konvensional dan analisis spasial untuk menangkap hubungan faktor epidemiologi dan pola penyebaran diare antarprovinsi. Pendekatan ini sesuai untuk penyakit berbasis lingkungan yang dipengaruhi oleh kondisi wilayah dan keterkaitan spasial antar daerah.

6 METODOLOGI

6.1 Sumber Data

Penelitian ini menggunakan data sekunder yang bersumber dari Badan Pusat Statistik (BPS) tahun 2018. Unit analisis adalah provinsi di Indonesia, dan seluruh data diintegrasikan dengan peta batas administratif provinsi untuk keperluan analisis spasial epidemiologi.

Variabel yang digunakan terdiri atas satu variabel dependen dan tiga variabel independen, sebagaimana disajikan pada Tabel berikut.

6.2 Variabel Penelitian

Variabel	Keterangan
Y	Jumlah Kasus Diare (jumlah kasus per provinsi)
X1	Jumlah Penduduk (ribu jiwa)
X2	Jumlah Tenaga Gizi (orang)
X3	Persentase Rumah Tangga dengan Akses Air Minum Layak (%)

6.3 Metode Analisis

Analisis dilakukan dengan pendekatan epidemiologi kuantitatif dan spasial, yang meliputi perhitungan ukuran kejadian penyakit, pengujian pola spasial, uji asumsi statistik, dan pemodelan regresi.

6.3.1 Ukuran Kejadian Penyakit (Prevalensi)

Untuk menggambarkan tingkat kejadian diare antarprovinsi, digunakan ukuran prevalensi, yang dihitung dengan rumus:

\[ \text{Prevalensi Diare} = \frac{\text{Jumlah Kasus Diare}}{\text{Jumlah Penduduk}} \times 100{,}000 \]

Prevalensi menunjukkan beban penyakit diare pada populasi dan memungkinkan perbandingan antarwilayah.

6.3.2 Uji Autokorelasi Spasial Global (Moran’s I)

Uji Moran’s I digunakan untuk mengetahui apakah distribusi kasus diare memiliki pola spasial tertentu (mengelompok, menyebar, atau acak). Rumus Moran’s I adalah:

\[ I = \frac{n}{S_0} \, \frac{\displaystyle \sum_{i}\sum_{j} w_{ij}\,(x_i - \bar{x})(x_j - \bar{x})} {\displaystyle \sum_{i} (x_i - \bar{x})^2} \]

\[ \begin{aligned} n &:\ \text{jumlah provinsi} \\ w_{ij} &:\ \text{bobot spasial antara provinsi } i \text{ dan } j \\ x_i &:\ \text{nilai kasus diare di provinsi } i \\ \bar{x} &:\ \text{rata-rata kasus diare} \\ S_0 &:\ \text{jumlah total bobot spasial, } S_0 = \sum_i \sum_j w_{ij} \end{aligned} \]

6.3.3 Uji Asumsi Statistik

Sebelum dilakukan pemodelan regresi, dilakukan beberapa uji asumsi, yaitu:

1. Uji Multikolinearitas (Variance Inflation Factor / VIF)
   Digunakan untuk memastikan tidak terdapat hubungan linier yang kuat antarvariabel independen.

   • VIF < 10 menunjukkan tidak terjadi multikolinearitas.

2. Uji Heteroskedastisitas (Breusch–Pagan Test)
   Digunakan untuk menguji kesamaan varians residual.

   • p-value > 0,05 menunjukkan tidak terdapat heteroskedastisitas.

3. Uji Normalitas Residual (Shapiro–Wilk Test)
   Digunakan untuk menguji apakah residual berdistribusi normal.

   • p-value > 0,05 menunjukkan residual berdistribusi normal.

6.3.3.1 Model Regresi OLS

Model Ordinary Least Squares (OLS) digunakan untuk menganalisis pengaruh faktor demografi dan lingkungan terhadap jumlah kasus diare. Model yang digunakan adalah:

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{3i} + \varepsilon_i \]

\[ \begin{aligned} Y_i &:\ \text{jumlah kasus diare pada wilayah ke-} i \\ X_{1i} &:\ \text{jumlah penduduk pada wilayah ke-} i \\ X_{2i} &:\ \text{jumlah tenaga gizi pada wilayah ke-} i \\ X_{3i} &:\ \text{persentase rumah dengan akses air layak pada wilayah ke-} i \\ \varepsilon_i &:\ \text{error term acak dengan } E(\varepsilon_i)=0 \end{aligned} \]

6.4 Alur Kerja Penelitian

Penelitian dimulai dengan pengumpulan data kasus diare, jumlah penduduk, jumlah tenaga gizi, dan akses air minum layak dari BPS tahun 2018. Data kemudian digabungkan dengan peta batas administratif provinsi. Selanjutnya dilakukan perhitungan prevalensi dan analisis deskriptif. Uji Moran’s I digunakan untuk mendeteksi pola spasial kasus diare. Setelah itu dilakukan uji asumsi dan pemodelan regresi OLS untuk menganalisis pengaruh faktor demografi dan lingkungan terhadap kejadian diare. Hasil analisis kemudian diinterpretasikan dalam konteks epidemiologi dan spasial.

7 HASIL DAN PEMBAHASAN

7.1 Data

library(spdep)

## Warning: package 'spdep' was built under R version 4.4.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.4.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Warning: package 'sf' was built under R version 4.4.3

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

library(sp)

## Warning: package 'sp' was built under R version 4.4.3

library(sf)
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(spData)
library(geodata)

## Warning: package 'geodata' was built under R version 4.4.3

## Loading required package: terra

## Warning: package 'terra' was built under R version 4.4.3

## terra 1.8.60

library(tmap)

## Warning: package 'tmap' was built under R version 4.4.3

library(skimr)

## Warning: package 'skimr' was built under R version 4.4.3

library(psych)

## Warning: package 'psych' was built under R version 4.4.3

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:terra':
## 
##     describe, distance, rescale

## The following object is masked from 'package:car':
## 
##     logit

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.4.3

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(grid)

## 
## Attaching package: 'grid'

## The following object is masked from 'package:terra':
## 
##     depth

library(viridis)

## Warning: package 'viridis' was built under R version 4.4.3

## Loading required package: viridisLite

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.4.3

## Loading required package: Matrix

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.4.2

## Warning: package 'readr' was built under R version 4.4.2

## Warning: package 'stringr' was built under R version 4.4.3

## Warning: package 'forcats' was built under R version 4.4.2

## Warning: package 'lubridate' was built under R version 4.4.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.2
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%()         masks ggplot2::%+%()
## ✖ psych::alpha()       masks ggplot2::alpha()
## ✖ gridExtra::combine() masks dplyr::combine()
## ✖ tidyr::expand()      masks Matrix::expand()
## ✖ tidyr::extract()     masks terra::extract()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tidyr::pack()        masks Matrix::pack()
## ✖ car::recode()        masks dplyr::recode()
## ✖ purrr::some()        masks car::some()
## ✖ tidyr::unpack()      masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:terra':
## 
##     area
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library(broom)
library(stringr)

# SET WORKING DIRECTORY
# ===========================================
setwd("C:\\Users\\Ghozzy\\OneDrive\\Documents\\SMESTER 5\\Epidem UAS")

# Baca shapefile provinsi Indonesia (ADM1)
indo_adm1 <- st_read("IDN_AdminBoundaries_candidate.gdb",
                     layer = "idn_admbnda_adm1_bps_20200401")

## Reading layer `idn_admbnda_adm1_bps_20200401' from data source 
##   `C:\Users\Ghozzy\OneDrive\Documents\SMESTER 5\Epidem UAS\IDN_AdminBoundaries_candidate.gdb' 
##   using driver `OpenFileGDB'

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts.  The
## processing may be really slow.  You can skip the processing by setting
## METHOD=SKIP.

## Simple feature collection with 34 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84

# Cek nama kolom provinsi
names(indo_adm1)

##  [1] "admin1Name_en"     "admin1Pcode"       "admin1RefName"    
##  [4] "admin1AltName1_en" "admin1AltName2_en" "admin0Name_en"    
##  [7] "admin0Pcode"       "date"              "validOn"          
## [10] "validTo"           "Shape_Length"      "Shape_Area"       
## [13] "Shape"

data <- read.csv("DATA.csv", sep=";", dec=",")

# Samakan huruf besar kecil
data$Provinsi <- toupper(data$Provinsi)
indo_adm1$provinsi_sf <- toupper(indo_adm1$admin1Name_en)

data$Provinsi[data$Provinsi == "METRO JAYA"] <- "DKI JAKARTA"
data$Provinsi[data$Provinsi == "KEP. RIAU"] <- "KEPULAUAN RIAU"
data$Provinsi[data$Provinsi == "KEP. BANGKA BELITUNG"] <- "KEPULAUAN BANGKA BELITUNG"

indo_merged <- indo_adm1 %>%
  left_join(data, by = c("provinsi_sf" = "Provinsi"))

7.2 Insidensi

# ===========================
# INSIDENSI DIARE
# ===========================

indo_merged <- indo_merged %>%
  mutate(
    Diare = as.numeric(gsub(",", ".", Diare)),
    Kepadatan.Ribu. = as.numeric(gsub(",", ".", Kepadatan.Ribu.)),
    population = Kepadatan.Ribu. * 1000,
    incidence_per100k = (Diare / population) * 100000
  )

# Ringkasan insidensi
summary(indo_merged$incidence_per100k)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     231    1126    1529    1546    1956    2922

7.3 Peta Deskriptif

7.3.1 Diare

ggplot(data = indo_merged) +
  geom_sf(aes(fill = Diare),
          color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma", direction = -1,
                       name = "Jumlah Kasus Diare") +
  labs(
    title = "Peta Persebaran Kasus Diare di Indonesia",
    subtitle = "Jumlah Kasus Diare per Provinsi",
    caption = "Sumber: BPS"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

7.3.2 Kepadatan Penduduk

ggplot(data = indo_merged) +
  geom_sf(aes(fill = Kepadatan.Ribu.),
          color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma", direction = -1,
                       name = "Kepadatan Penduduk\n(per ribu)") +
  labs(
    title = "Peta Kepadatan Penduduk di Indonesia",
    subtitle = "Kepadatan Penduduk per Provinsi",
    caption = "Sumber: Data_UAS.csv"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

7.3.3 Jumlah Tenaga Gizi

ggplot(data = indo_merged) +
  geom_sf(aes(fill = Tenaga.Gizi.),
          color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma", direction = -1,
                       name = "Jumlah Tenaga Gizi") +
  labs(
    title = "Peta Persebaran Tenaga Gizi di Indonesia",
    subtitle = "Jumlah Tenaga Gizi per Provinsi",
    caption = "Sumber: Data_UAS.csv"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

7.3.4 Persentase Rumah tangga dengan air yang layak

ggplot(data = indo_merged) +
  geom_sf(
    aes(fill = PersentaseRumahAirLayak),
    color = "white",
    linewidth = 0.3
  ) +
  scale_fill_viridis_c(
    option = "plasma",
    direction = -1,
    name = "Persentase Rumah\nAir Layak (%)"
  ) +
  labs(
    title = "Peta Akses Air Layak di Indonesia",
    subtitle = "Persentase Rumah Tangga dengan Air Layak per Provinsi",
    caption = "Sumber: Data_UAS.csv"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

7.3.5 Peta Insidensi

indo_merged <- indo_merged %>%
  rename(Provinsi = provinsi_sf)

head(indo_merged[, c("Provinsi", "Diare", "population", "incidence_per100k")])

## Simple feature collection with 6 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -8.84919 xmax: 115.7125 ymax: 6.07693
## Geodetic CRS:  WGS 84
##                     Provinsi  Diare population incidence_per100k
## 1                       ACEH  76753    5281300          1453.297
## 2                       BALI  68142    4292200          1587.577
## 3                     BANTEN 248242   12689700          1956.248
## 4                   BENGKULU  21313    1963300          1085.570
## 5 DAERAH ISTIMEWA YOGYAKARTA  68043    3802900          1789.240
## 6                DKI JAKARTA 305841   10467600          2921.787
##                            Shape
## 1 MULTIPOLYGON (((97.39178 2....
## 2 MULTIPOLYGON (((115.1251 -8...
## 3 MULTIPOLYGON (((105.5498 -6...
## 4 MULTIPOLYGON (((102.3862 -5...
## 5 MULTIPOLYGON (((110.8198 -8...
## 6 MULTIPOLYGON (((106.8768 -6...

# Mode tampilan tmap
tmap_mode("plot")

## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`

# --- Persiapan & pembersihan data ---
# Pastikan kolom terkonversi ke numeric (tangani tanda koma desimal)
indo_merged$Diare <- as.numeric(gsub(",", ".", as.character(indo_merged$Diare)))
if ("Kepadatan.Ribu." %in% names(indo_merged)) {
  indo_merged$Kepadatan.Ribu. <- as.numeric(gsub(",", ".", as.character(indo_merged$Kepadatan.Ribu.)))
}

# Cari kolom populasi yang tersedia (prioritas: explicit pop, lalu Kepadatan.Ribu.)
# Asumsi: Kepadatan.Ribu. = jumlah penduduk dalam ribuan (contoh: 5281.3 -> 5,281,300)
pop_col <- NULL
if ("Population" %in% names(indo_merged)) {
  pop_col <- "Population"
} else if ("Pop" %in% names(indo_merged)) {
  pop_col <- "Pop"
} else if ("Kepadatan.Ribu." %in% names(indo_merged)) {
  # konversi ke total jiwa
  indo_merged$pop_total <- indo_merged$Kepadatan.Ribu. * 1000
  pop_col <- "pop_total"
} else {
  stop("Tidak menemukan kolom populasi. Tambahkan kolom populasi atau pastikan ada 'Kepadatan.Ribu.'")
}

# --- Hitung insidensi per 100.000 ---
# Hindari pembagian dengan nol / NA
indo_merged$incidence_per100k <- NA_real_
valid_idx <- !is.na(indo_merged$Diare) & !is.na(indo_merged[[pop_col]]) & indo_merged[[pop_col]] > 0
indo_merged$incidence_per100k[valid_idx] <-
  (indo_merged$Diare[valid_idx] / indo_merged[[pop_col]][valid_idx]) * 100000

# Opsional: pembulatan (untuk legend/label)
indo_merged$incidence_per100k_label <- round(indo_merged$incidence_per100k, 1)

# --- Peta tmap ---
tm_shape(indo_merged) +
  tm_polygons(
    col = "incidence_per100k",
    palette = "Reds",
    title = "Insidensi Diare\nper 100.000 penduduk",
    style = "quantile",    # bisa ganti ke "pretty" atau "jenks"
    n = 6,                 # jumlah kelas
    textNA = "Data tidak tersedia"
  ) +
  tm_borders(alpha = 0.4) +
  tm_layout(
    title = "Sebaran Insidensi Diare di Indonesia",
    legend.outside = TRUE,
    frame = FALSE,
    legend.outside.size = 0.25
  ) +
  tm_compass(type = "8star", position = c("left", "top")) +
  tm_scale_bar(position = c("left", "bottom"))

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values'),
##   'textNA' (rename to 'label.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

Interpretasi :

Peta sebaran menunjukkan bahwa insidensi diare per 100.000 penduduk di Indonesia pada tahun 2018 bervariasi antarprovinsi dan tidak tersebar secara merata. Beberapa wilayah tampak memiliki insidensi lebih tinggi yang mengindikasikan beban penyakit diare yang besar, kemungkinan berkaitan dengan faktor lingkungan, kepadatan penduduk, serta keterbatasan akses air minum layak dan sanitasi. Sementara itu, provinsi dengan insidensi lebih rendah menunjukkan kondisi kesehatan lingkungan yang relatif lebih baik. Pola perbedaan antarwilayah ini mengisyaratkan adanya keterkaitan spasial yang perlu dianalisis lebih lanjut untuk menentukan wilayah prioritas dalam upaya pencegahan dan pengendalian diare.

7.4 Global MORAN I

# ===============================
# GLOBAL MORAN'S I
# INSIDENSI DIARE INDONESIA
# ===============================

# Library
library(sf)
library(sp)
library(spdep)
library(dplyr)

# Konversi sf ke Spatial
map_sp <- as_Spatial(indo_merged)

# Membuat neighbour (Queen contiguity)
nb <- poly2nb(map_sp)

## Warning in poly2nb(map_sp): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(map_sp): neighbour object has 11 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

# Spatial weight matrix (atasi wilayah tanpa tetangga)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Uji Global Moran's I
moran.test(
  map_sp$incidence_per100k,
  lw,
  zero.policy = TRUE
)

## 
##  Moran I test under randomisation
## 
## data:  map_sp$incidence_per100k  
## weights: lw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.21465, p-value = 0.415
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.001409843      -0.037037037       0.032082859

Interpretasi :

Nilai Moran’s I = 0,0014 sangat mendekati nol dan p-value = 0,415 (> 0,05). Ini berarti tidak terdapat autokorelasi spasial yang signifikan pada insidensi diare antarprovinsi di Indonesia. Dengan kata lain, pola sebaran insidensi diare bersifat acak (random), sehingga provinsi dengan insidensi tinggi atau rendah tidak cenderung mengelompok secara spasial.

7.5 Model OLS

# ===============================
# REGRESI OLS INSIDENSI DIARE
# ===============================

library(dplyr)
library(car)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following object is masked from 'package:terra':
## 
##     time<-

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# Pastikan variabel numerik
indo_merged <- indo_merged %>%
  mutate(
    incidence_per100k = as.numeric(incidence_per100k),
    Kepadatan.Ribu. = as.numeric(gsub(",", ".", Kepadatan.Ribu.)),
    Tenaga.Gizi. = as.numeric(gsub(",", ".", Tenaga.Gizi.)),
    PersentaseRumahAirLayak = as.numeric(gsub(",", ".", PersentaseRumahAirLayak))
  )

# Model OLS
model_ols <- lm(
  incidence_per100k ~ Kepadatan.Ribu. + Tenaga.Gizi. + PersentaseRumahAirLayak,
  data = indo_merged
)

# Ringkasan model
summary(model_ols)

## 
## Call:
## lm(formula = incidence_per100k ~ Kepadatan.Ribu. + Tenaga.Gizi. + 
##     PersentaseRumahAirLayak, data = indo_merged)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1347.06  -345.44   -62.68   449.39  1316.54 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             1236.49724  823.17013   1.502   0.1435  
## Kepadatan.Ribu.           -0.02850    0.01554  -1.834   0.0766 .
## Tenaga.Gizi.               0.67170    0.27984   2.400   0.0228 *
## PersentaseRumahAirLayak    1.35037   11.37070   0.119   0.9063  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 591.9 on 30 degrees of freedom
## Multiple R-squared:  0.1701, Adjusted R-squared:  0.08713 
## F-statistic:  2.05 on 3 and 30 DF,  p-value: 0.128

Interpretasi :

Secara simultan, variabel kepadatan penduduk, tenaga gizi, dan persentase rumah dengan akses air layak tidak berpengaruh signifikan terhadap insidensi diare (p-value uji F = 0,128 > 0,05). Artinya, ketiga variabel tersebut belum mampu menjelaskan variasi insidensi diare secara bersama-sama pada tingkat provinsi.

• Kepadatan Penduduk memiliki koefisien negatif dan signifikan pada taraf 10% (p = 0,076). Ini menunjukkan bahwa peningkatan kepadatan penduduk cenderung berkaitan dengan penurunan insidensi diare, meskipun pengaruhnya lemah.
• Tenaga Gizi berpengaruh positif dan signifikan pada taraf 5% (p = 0,023). Artinya, provinsi dengan jumlah tenaga gizi lebih banyak cenderung memiliki insidensi diare yang lebih tinggi, yang kemungkinan mencerminkan respon terhadap tingginya kasus, bukan sebagai penyebab langsung.
• Persentase Rumah dengan Air Layak tidak berpengaruh signifikan terhadap insidensi diare (p = 0,906).

Nilai R-squared sebesar 0,17 menunjukkan bahwa model hanya mampu menjelaskan sekitar 17% variasi insidensi diare, sementara sisanya dipengaruhi oleh faktor lain di luar model, seperti perilaku higienis, sanitasi lingkungan, atau akses layanan kesehatan.

7.5.1 Uji Asumsi

# -------------------------------
# UJI ASUMSI DASAR
# -------------------------------

# Multikolinearitas
vif(model_ols)

##         Kepadatan.Ribu.            Tenaga.Gizi. PersentaseRumahAirLayak 
##                2.806380                2.870565                1.041940

# Heteroskedastisitas (Breusch-Pagan)
bptest(model_ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 1.3203, df = 3, p-value = 0.7243

# Normalitas residual (Shapiro-Wilk)
shapiro.test(residuals(model_ols))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_ols)
## W = 0.98433, p-value = 0.8955

Hasil uji Breusch–Pagan menghasilkan p-value = 0,7243 (> 0,05). Ini berarti tidak terdapat heteroskedastisitas, sehingga varians residual bersifat homoskedastis dan asumsi OLS terpenuhi.

Uji Shapiro–Wilk menunjukkan p-value = 0,8955 (> 0,05). Artinya, residual berdistribusi normal, sehingga asumsi normalitas residual dalam regresi OLS terpenuhi.

7.6 Uji Korelasi Spearman

# ===============================
# UJI KORELASI SPEARMAN
# ===============================

# Pastikan variabel numerik
indo_merged <- indo_merged %>%
  mutate(
    incidence_per100k = as.numeric(incidence_per100k),
    Kepadatan.Ribu. = as.numeric(gsub(",", ".", Kepadatan.Ribu.)),
    Tenaga.Gizi. = as.numeric(gsub(",", ".", Tenaga.Gizi.)),
    PersentaseRumahAirLayak = as.numeric(gsub(",", ".", PersentaseRumahAirLayak))
  )

# Korelasi Kepadatan Penduduk vs Insidensi Diare
cor.test(
  indo_merged$Kepadatan.Ribu.,
  indo_merged$incidence_per100k,
  method = "spearman",
  use = "complete.obs"
)

## 
##  Spearman's rank correlation rho
## 
## data:  indo_merged$Kepadatan.Ribu. and indo_merged$incidence_per100k
## S = 4640, p-value = 0.09509
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2910619

# Korelasi Tenaga Gizi vs Insidensi Diare
cor.test(
  indo_merged$Tenaga.Gizi.,
  indo_merged$incidence_per100k,
  method = "spearman",
  use = "complete.obs"
)

## 
##  Spearman's rank correlation rho
## 
## data:  indo_merged$Tenaga.Gizi. and indo_merged$incidence_per100k
## S = 5264, p-value = 0.2662
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1957219

# Korelasi Akses Air Layak vs Insidensi Diare
cor.test(
  indo_merged$PersentaseRumahAirLayak,
  indo_merged$incidence_per100k,
  method = "spearman",
  use = "complete.obs"
)

## 
##  Spearman's rank correlation rho
## 
## data:  indo_merged$PersentaseRumahAirLayak and indo_merged$incidence_per100k
## S = 6254, p-value = 0.8024
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04446142

Interpretasi :

• Kepadatan Penduduk memiliki korelasi positif lemah dengan insidensi diare (ρ = 0,29) dan tidak signifikan pada taraf 5% (p = 0,095), namun menunjukkan kecenderungan hubungan pada taraf 10%.

• Tenaga Gizi menunjukkan korelasi positif lemah dengan insidensi diare (ρ = 0,20) dan tidak signifikan (p = 0,266).

• Persentase Rumah dengan Air Layak memiliki korelasi sangat lemah dan tidak signifikan terhadap insidensi diare (ρ = 0,04; p = 0,802).

Kesimpulan: Tidak terdapat hubungan yang kuat dan signifikan antara ketiga variabel dengan insidensi diare pada tingkat provinsi; hubungan yang ada cenderung lemah dan tidak bermakna secara statistik.

8 KESIMPULAN

Berdasarkan hasil regresi OLS, dapat disimpulkan bahwa secara umum faktor kepadatan penduduk, jumlah tenaga gizi, dan persentase rumah tangga dengan akses air layak belum mampu menjelaskan variasi insidensi diare antarprovinsi secara kuat. Hal ini terlihat dari hasil uji simultan (uji F) yang tidak signifikan serta nilai R-squared yang relatif rendah, yaitu sekitar 17%. Artinya, sebagian besar variasi kasus diare masih dipengaruhi oleh faktor lain di luar model, seperti perilaku hidup bersih dan sehat, kondisi sanitasi, kualitas lingkungan, serta akses dan pemanfaatan layanan kesehatan. Dengan demikian, model OLS ini lebih bersifat eksploratif dan belum cukup kuat untuk digunakan sebagai dasar prediksi yang komprehensif.

Secara parsial, kepadatan penduduk menunjukkan hubungan negatif yang lemah terhadap insidensi diare, yang mengindikasikan bahwa wilayah dengan penduduk lebih padat tidak selalu memiliki kasus diare yang lebih tinggi, kemungkinan karena infrastruktur dan layanan kesehatan yang relatif lebih baik. Sementara itu, variabel tenaga gizi berpengaruh positif dan signifikan, yang lebih mencerminkan respons sistem kesehatan terhadap tingginya kasus diare daripada sebagai faktor penyebab langsung. Adapun persentase rumah dengan akses air layak tidak menunjukkan pengaruh yang berarti, yang mengindikasikan bahwa ketersediaan air layak saja belum cukup menekan insidensi diare tanpa didukung praktik sanitasi dan perilaku higienis yang baik. Temuan ini menegaskan bahwa permasalahan diare bersifat kompleks dan memerlukan pendekatan multidimensi dalam analisis maupun kebijakan penanggulangannya.

9 DAFTAR PUSTAKA

[1] Badan Pusat Statistik (BPS). Jumlah Kasus Penyakit Menurut Provinsi dan Jenis Penyakit, 2018.
https://www.bps.go.id/id/statistics-table/3/YTA1Q1ptRmhUMEpXWTBsQmQyZzBjVzgwUzB4aVp6MDkjMw==/jumlah-kasus-penyakit-menurut-provinsi-dan-jenis-penyakit--2018.html

[2] Badan Pusat Statistik (BPS). Penduduk, Laju Pertumbuhan Penduduk, Distribusi Persentase Penduduk, Kepadatan Penduduk, dan Rasio Jenis Kelamin Menurut Provinsi, 2018.
https://www.bps.go.id/id/statistics-table/3/V1ZSbFRUY3lTbFpEYTNsVWNGcDZjek53YkhsNFFUMDkjMw==/penduduk--laju-pertumbuhan-penduduk--distribusi-persentase-penduduk--kepadatan-penduduk--rasio-jenis-kelamin-penduduk-menurut-provinsi--2024.html?year=2018

[3] Badan Pusat Statistik (BPS). Jumlah Tenaga Kesehatan Menurut Provinsi, 2018.
https://www.bps.go.id/id/statistics-table/3/YVdwSFJHRjRVVkJqWlRWRU9EQkhNVFY0UjB4VVVUMDkjMyMwMDAw/jumlah-tenaga-kesehatan-menurut-provinsi.html?year=2018

[4] Badan Pusat Statistik (BPS). Persentase Rumah Tangga yang Memiliki Akses terhadap Sumber Air Minum Layak Menurut Provinsi, 2018.
https://www.bps.go.id/id/statistics-table/2/ODU0IzI=/persentase-rumah-tangga-yang-memiliki-akses-terhadap-sumber-air-minum-layak-menurut-provinsi-dan-klasifikasi-desa--persen-.html

[5] World Health Organization (WHO). Diarrhoeal Disease: Key Facts.
https://www.who.int/news-room/fact-sheets/detail/diarrhoeal-disease

[6] Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
https://doi.org/10.1111/j.1538-4632.1995.tb00338.x

[7] Cliff, A. D., & Ord, J. K. (1981). Spatial Processes: Models & Applications. London: Pion.