COVER

Judul
Analisis Spasial Kasus Pneumonia Pada Balita di Jawa Timur Tahun 2022

Nama
Stephani Julieta

Program Studi
S-1 Statistika

NPM
140610230050

Mata Kuliah
Spasial

Dosen Pembimbing Mata Kuliah
I Gede Nyoman Mindra Jaya, P.hd

Abstrak

(Isi abstrak singkat ±150–200 kata: latar belakang, tujuan, metode, hasil utama, dan kesimpulan. Akan diisi setelah seluruh analisis selesai.)

1 Pendahuluan

1.1 Latar Belakang

Pneumonia merupakan salah satu penyebab utama kesakitan dan kematian pada balita di berbagai negara, termasuk Indonesia. Pada kelompok usia balita, pneumonia berkontribusi besar terhadap beban penyakit karena sistem imun yang belum matang, kerentanan terhadap infeksi saluran napas, serta faktor lingkungan dan sosial ekonomi yang dapat memperburuk risiko penularan dan keparahan. Selain berdampak pada kesehatan anak, pneumonia pada balita juga menimbulkan beban ekonomi bagi keluarga dan sistem layanan kesehatan akibat meningkatnya kebutuhan kunjungan pelayanan, pengobatan, dan perawatan.

Di tingkat daerah, beban pneumonia balita tidak selalu merata. Provinsi Jawa Timur sebagai salah satu provinsi dengan jumlah penduduk besar serta keragaman karakteristik wilayah—mulai dari kawasan perkotaan padat, wilayah industri, hingga kawasan perdesaan dan kepulauan—memiliki potensi variasi risiko pneumonia yang berbeda antar kabupaten/kota. Perbedaan tersebut dapat dipengaruhi oleh faktor kepadatan penduduk, kualitas hunian dan ventilasi rumah, akses air bersih dan sanitasi, cakupan imunisasi, status gizi, perilaku pencarian pengobatan, hingga aksesibilitas dan kualitas layanan kesehatan. Selain itu, faktor lingkungan seperti kualitas udara (misalnya polusi perkotaan, pembakaran sampah, atau asap rokok dalam rumah) dan kondisi iklim setempat juga dapat berperan dalam memengaruhi pola kejadian pneumonia.

Pada tahun 2022, kondisi pemulihan layanan kesehatan pascapandemi turut menjadi perhatian karena dapat memengaruhi cakupan layanan dasar (seperti imunisasi, pemantauan tumbuh kembang, dan penemuan kasus). Situasi ini berpotensi memunculkan perubahan pola pelaporan dan distribusi kasus antarwilayah. Oleh karena itu, pemetaan dan analisis distribusi kasus pneumonia balita tahun 2022 penting dilakukan untuk memahami apakah terdapat konsentrasi kasus pada wilayah tertentu (cluster) dan apakah wilayah dengan beban tinggi berdekatan secara geografis, yang dapat mengindikasikan adanya pengaruh faktor spasial.

Pendekatan analisis spasial menjadi relevan karena mampu menggambarkan pola persebaran kasus secara geografis serta mengidentifikasi wilayah prioritas secara lebih presisi dibandingkan analisis nonspasial. Melalui pemetaan tematik dan uji autokorelasi spasial (misalnya Moran’s I) maupun analisis klaster (misalnya LISA/hotspot), penelitian dapat menunjukkan apakah kasus pneumonia balita cenderung mengelompok pada lokasi tertentu atau menyebar secara acak. Hasil analisis spasial ini dapat menjadi dasar untuk intervensi kesehatan masyarakat yang lebih tepat sasaran—misalnya penguatan penemuan kasus dan tata laksana di kabupaten/kota tertentu, peningkatan cakupan imunisasi dan edukasi pencegahan di wilayah berisiko, serta perbaikan determinan lingkungan dan sosial yang berkontribusi terhadap tingginya kejadian pneumonia.

Berdasarkan uraian tersebut, analisis spasial kasus pneumonia pada balita di Jawa Timur tahun 2022 perlu dilakukan untuk (1) memetakan distribusi kasus menurut kabupaten/kota, (2) menilai ada tidaknya autokorelasi spasial, dan (3) mengidentifikasi klaster wilayah risiko tinggi/rendah sebagai dasar penyusunan kebijakan dan program pencegahan pneumonia balita yang lebih efektif dan berbasis wilayah.

1.2 Identifikasi Masalah

RUMUSAN MASALAH

1.3 Tujuan Penelitian

TUJUAN

1.4 Batasan Penelitian

Penelitian ini hanya mencakup 38 kabupaten/kota di Provinsi Jawa Timur menggunakan data sekunder tahun 2022, sehingga hasil analisis yang diperoleh tidak dapat digeneralisasikan untuk provinsi lain di Indonesia. Penelitian ini juga bersifat cross-sectional dan dilakukan tanpa mempertimbangkan dinamika waktu kasus Pneumonia dari waktu ke waktu.

2 Tinjauan Pustaka

2.1 Dependensi Spasial

Dependensi spasial menyatakan bahwa nilai suatu variabel pada satu wilayah tidak independen dari nilai variabel yang sama pada wilayah tetangga. Hal ini mencerminkan bahwa “barang-apa yang berada di dekat lebih berkaitan dibanding yang jauh” sesuai dengan Tobler’s First Law of Geography. Pada penelitian ini, dependensi spasial dapat muncul karena jumlah penduduk, kondisi lingkungan, akses layanan kesehatan antar daerah, dan kebiasaan hidup.

Beberapa implikasi dependensi spasial diantara lain :

Variabel dalam satu wilayah dapat memengaruhi variabel di wilayah tetangga (spillover).
Observasi pada satu unit wilayah tidak independen, maka terjadi pelanggaran asumsi klasik regresi.
Memerlukan pendekatan analisis spasial untuk menangkap efek ini.

2.2 Autokorelasi Spasial

Autokorelasi spasial mengukur sejauh mana suatu wilayah memiliki kesamaan nilai dengan wilayah di sekitarnya. Autokorelasi positif menunjukkan bahwa wilayah dengan nilai tinggi (atau rendah) cenderung berdekatan dengan wilayah lain yang juga memiliki nilai serupa (cluster), sedangkan autokorelasi negatif menunjukkan pola penyebaran yang saling berbeda (dispersion).

Ukuran utama yang digunakan untuk mengukur autokorelasi spasial adalah Moran’s I dan Local Indicators of Spatial Association (LISA) dengan formulasi sebagai berikut.

2.2.1 Rumus Moran’s I (Global Spasial Autocorrelation)

\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij} (y_i - \bar{y})(y_j - \bar{y})} {\sum_i (y_i - \bar{y})^2} \]

Keterangan:

\(n\) = jumlah observasi
\(y_i, y_j\) = nilai variabel pada lokasi i dan j
\(\bar{y}\) = rata-rata keseluruhan
\(w_{ij}\) = bobot spasial antara wilayah i dan j

Nilai \(I > 0\) menunjukkan autokorelasi spasial positif,
\(I < 0\) menunjukkan negatif, dan \(I = 0\) berarti tidak ada autokorelasi
(Anselin, 1995).

2.2.2 Rumus Geary’s C

\[ C = \frac{(n - 1)\sum_i \sum_j w_{ij}(y_i - y_j)^2} {2 \sum_i \sum_j w_{ij} \sum_i (y_i - \bar{y})^2} \]

Nilai \(C < 1\) menunjukkan bahwa terdapat autokorelasi positif,
sedangkan \(C > 1\) menunjukkan autokorelasi negatif.

2.2.3 Rumus LISA (Local Indicators Spatial Autocorrelation)

\[ I_i = (y_i - \bar{y}) \sum_j w_{ij} (y_j - \bar{y}) \]

LISA digunakan untuk mengidentifikasi klaster lokal, seperti: - High-High (daerah tinggi dikelilingi tinggi) - Low-Low - High-Low
- Low-High

2.2.4 Scatterplot Moran & Hotspot

\[ G_i^* = \frac {\sum_j w_{ij}x_j - \bar{X}\sum_j w_{ij}} {S \sqrt{ \frac{ \left[ n\sum_j w_{ij}^2 - (\sum_j w_{ij})^2 \right] }{n-1} }} \]

dimana: - \(x_j\) = nilai variabel tetangga - \(S\) = standar deviasi

Hotspot menandakan klaster tinggi yang signifikan secara spasial.

2.3 Model Ekonometrika Spasial

Model ekonometrik spasial digunakan untuk mengatasi pelanggaran asumsi independensi residual dalam model regresi konvensional. Menurut Anselin (1988),model-model ini memperhitungkan adanya efek spasial yang mungkin berasal dari pengaruh tetangga (spatial lag) atau dari error yang berkorelasi secara spasial (spatial error).

Beberapa model utama yang digunakan meliputi:

2.3.1 Spatial Autoregressive Model

\[ y = \rho W y + X \beta + \varepsilon \]

Keterangan:

\(\rho\) = koefisien autokorelasi spasial pada variabel dependen
\(W y\) = lag spasial dari \(y\)
\(X\beta\) = variabel independen dan koefisiennya
\(\varepsilon\) = error

Model ini menunjukkan bahwa nilai \(y\) pada suatu lokasi dipengaruhi oleh nilai \(y\) pada lokasi di sekitarnya.

2.3.2 Spatial Error Model (SEM)

\[ y = X\beta + u, \quad u = \lambda W u + \varepsilon \]

Keterangan:

\(y\) = vektor variabel dependen
\(X\) = matriks variabel independen
\(\beta\) = vektor koefisien regresi
\(u\) = komponen error yang mengandung dependensi spasial
\(\lambda\) = koefisien spatial error dependence
\(W\) = matriks bobot spasial
\(\varepsilon\) = error

Model ini digunakan ketika efek spasial muncul pada residual, bukan pada variabel dependen.
Nilai \(\lambda\) menunjukkan tingkat ketergantungan error spasial.

2.3.3 IDW (Inverse Distance Weighting)

\[ \hat{z}(s_0)= \frac{\sum_{i=1}^{n} z(s_i)\, d(s_0,s_i)^{-p}} {\sum_{i=1}^{n} d(s_0,s_i)^{-p}} \]

2.3.4 Spatial Durbin Model

\[ y = \rho W y + X \beta + W X \theta + \varepsilon \]

Keterangan:

\(\rho\) = koefisien autokorelasi spasial pada variabel dependen
\(W y\) = lag spasial dari \(y\)
\(\theta\) = vektor koefisien untuk lag spasial dari variabel independen
\(W X\) = lag spasial dari variabel independen
\(\varepsilon\) = error

Model ini menggabungkan pengaruh spasial baik pada variabel dependen maupun independen, sehingga memberikan fleksibilitas lebih dalam menangkap hubungan spasial antar wilayah.

2.3.5 Spatial Durbin Error Model

\[ y = X\beta + WX\theta + u, \quad u = \lambda W u + \varepsilon \]

Keterangan:

\(X\beta\) = pengaruh langsung dari variabel independen
\(WX\theta\) = pengaruh tidak langsung dari variabel independen tetangga
\(u\) = komponen error yang juga memiliki autokorelasi spasial
\(\lambda\) = koefisien dependensi spasial pada error

SDEM menggabungkan lag spasial dari variabel independen serta error yang bersifat spasial, menjadikannya cocok untuk fenomena lingkungan dan sosial seperti penyebaran penyakit.

3 Metodologi Penelitian

3.1 Data dan Unit Spasial

SUMBER DATA

3.2 Metode Analisis

METODE ANALISIS

3.3 Alur Kerja Penelitian

ALUR KERJA

4 Hasil dan Pembahasan

4.1 Import Data

## Package
library(readxl)

## Warning: package 'readxl' was built under R version 4.3.3

library(sf)

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

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## 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(stringr)
library(spdep)

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

## Loading required package: spData

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

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

library(tmap)
library(sf)
library(gstat)

## Warning: package 'gstat' was built under R version 4.3.3

library(stars)

## Warning: package 'stars' was built under R version 4.3.3

## Loading required package: abind

## Warning: package 'abind' was built under R version 4.3.3

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.3.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(lmtest)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

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

library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.3.3

## 
## Attaching package: 'kableExtra'

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

library(ggplot2)
library(car)

## Warning: package 'car' was built under R version 4.3.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.3

## 
## Attaching package: 'car'

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

## 1) Baca Excel
data <- read_excel("C:/Users/ASUS/Downloads/pneumonia jatim data.xlsx")

attach(data)

## 2) Baca shapefile (pastikan 4 file .shp .shx .dbf .prj satu folder)
jatim_shp <- st_read("C:/Users/ASUS/Downloads/jatim_shp")

## Reading layer `Kabupaten-Kota (Provinsi Jawa Timur)' from data source 
##   `C:\Users\ASUS\Downloads\jatim_shp' using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 110.8987 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84

## 3) Buat kunci join yang seragam (TIDAK mengubah nama asli, hanya tambah kolom bantu)
data2 <- data %>%
  mutate(
    join_kabkota = `Kabupaten/Kota` %>%
      str_to_upper() %>%
      str_trim() %>%
      str_replace_all("^KAB\\.?\\s+", "KABUPATEN ") %>%  # KAB. -> KABUPATEN
      str_replace_all("^KOTA\\s+", "KOTA ")
  )

jatim_shp2 <- jatim_shp %>%
  mutate(
    join_kabkota = `Kabupaten/Kota` %>%
      str_to_upper() %>%
      str_trim() %>%
      str_replace_all("^KAB\\.?\\s+", "KABUPATEN ") %>%
      str_replace_all("^KOTA\\s+", "KOTA ")
  )
names(jatim_shp)

##  [1] "ID_0"      "COUNTRY"   "NAME_1"    "NL_NAME_1" "ID_2"      "NAME_2"   
##  [7] "VARNAME_2" "NL_NAME_2" "TYPE_2"    "ENGTYPE_2" "CC_2"      "HASC_2"   
## [13] "geometry"

## 4) Join jadi sf
jatim_data <- jatim_shp2 %>%
  left_join(data2, by = "join_kabkota")

## 5) Cek join (harus 0 NA untuk variabel penting)
stopifnot(sum(is.na(jatim_data$`Jumlah Kasus Pneumonia`)) == 0)
stopifnot(sum(is.na(jatim_data$`Jumlah Penduduk Balita`)) == 0)
stopifnot(sum(is.na(jatim_data$`Kepadatan Penduduk`)) == 0)
stopifnot(sum(is.na(jatim_data$`Cakupan Imunisasi`)) == 0)
stopifnot(sum(is.na(jatim_data$`Prevalensi Balita Gizi Kurang`)) == 0)

## 6) Buat Y prevalensi pneumonia balita per 100.000 balita (masuk ke jatim_data)
jatim_data <- jatim_data %>%
  mutate(
    prevalensi_pneumonia_balita =
      (`Jumlah Kasus Pneumonia` / `Jumlah Penduduk Balita`) * 100000
  )

summary(jatim_data$prevalensi_pneumonia_balita)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   329.9  2154.7  3183.2  3321.4  4278.3 10033.2

jatim_data_utm <- st_transform(jatim_data, 32749)

4.2 Peta Deskriptif dan Statistik Deskriptif

4.3 3) Statistik Deskriptif (yang kamu bilang belum ada)

4.3.1 3.1 Ringkasan numerik variabel (tabel rapi)

desc_tbl <- jatim_data_utm %>%
  st_drop_geometry() %>%
  select(
    prevalensi_pneumonia_balita,
    `Kepadatan Penduduk`,
    `Prevalensi Balita Gizi Kurang`,
    `Cakupan Imunisasi`
  ) %>%
  summarise(
    n = n(),
    mean_Y = mean(prevalensi_pneumonia_balita, na.rm = TRUE),
    sd_Y   = sd(prevalensi_pneumonia_balita, na.rm = TRUE),
    min_Y  = min(prevalensi_pneumonia_balita, na.rm = TRUE),
    q1_Y   = quantile(prevalensi_pneumonia_balita, 0.25, na.rm = TRUE),
    med_Y  = median(prevalensi_pneumonia_balita, na.rm = TRUE),
    q3_Y   = quantile(prevalensi_pneumonia_balita, 0.75, na.rm = TRUE),
    max_Y  = max(prevalensi_pneumonia_balita, na.rm = TRUE)
  )

kable(desc_tbl, digits = 3, caption = "Statistik deskriptif prevalensi pneumonia balita (per 100.000).") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Statistik deskriptif prevalensi pneumonia balita (per 100.000).
n	mean_Y	sd_Y	min_Y	q1_Y	med_Y	q3_Y	max_Y
38	3321.397	1831.936	329.927	2154.728	3183.153	4278.291	10033.2

### 3.2 Histogram Y (ringan dan aman render)

ggplot(st_drop_geometry(jatim_data_utm), aes(x = prevalensi_pneumonia_balita)) +
  geom_histogram(bins = 12) +
  labs(
    title = "Distribusi Prevalensi Pneumonia Balita",
    x = "Prevalensi per 100.000 balita",
    y = "Frekuensi"
  ) +
  theme_minimal()

## 7) Perbaikan spasial: transform ke proyeksi (UTM) agar tidak warning planar
##    Untuk Jawa Timur biasanya aman pakai UTM 49S (EPSG:32749).
##    Kalau hasilnya aneh, ganti ke EPSG:32750.
jatim_data_utm <- st_transform(jatim_data, 32749)

## 8) Bobot spasial (queen contiguity)
##    Gunakan snap untuk membantu polygon yang "hampir bersentuhan"
nb <- poly2nb(jatim_data_utm, queen = TRUE, snap = 100)  # kalau masih sub-graphs, coba 500 atau 1000

## Warning in poly2nb(jatim_data_utm, queen = TRUE, snap = 100): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

## 9) Moran’s I GLOBAL untuk Y
moran_global <- moran.test(jatim_data_utm$prevalensi_pneumonia_balita, lw, zero.policy = TRUE)
moran_global

## 
##  Moran I test under randomisation
## 
## data:  jatim_data_utm$prevalensi_pneumonia_balita  
## weights: lw    
## 
## Moran I statistic standard deviate = -3.3072, p-value = 0.9995
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.43742195       -0.02702703        0.01539902

## 10) Geary’s C GLOBAL untuk Y
geary_global <- geary.test(jatim_data_utm$prevalensi_pneumonia_balita, lw, zero.policy = TRUE)
geary_global

## 
##  Geary C test under randomisation
## 
## data:  jatim_data_utm$prevalensi_pneumonia_balita 
## weights: lw   
## 
## Geary C statistic standard deviate = -2.6068, p-value = 0.9954
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.45861944        1.00000000        0.03095312

## 11) LISA (Local Moran)
lisa <- localmoran(jatim_data_utm$prevalensi_pneumonia_balita, lw, zero.policy = TRUE)

# Tambahkan hasil LISA ke data
jatim_data_utm <- jatim_data_utm %>%
  mutate(
    lisa_I    = lisa[, "Ii"],
    lisa_EI   = lisa[, "E.Ii"],
    lisa_VI   = lisa[, "Var.Ii"],
    lisa_Z    = lisa[, "Z.Ii"],
    lisa_p    = lisa[, "Pr(z != E(Ii))"]
  )

## 12) Buat peta klaster LISA (HH, LL, HL, LH) + Non-signifikan
##     Metode: gunakan mean sebagai pemisah High/Low (umum untuk tugas kuliah)
y <- jatim_data_utm$prevalensi_pneumonia_balita
y_mean <- mean(y, na.rm = TRUE)

# spatial lag dari y
lag_y <- lag.listw(lw, y, zero.policy = TRUE)

jatim_data_utm <- jatim_data_utm %>%
  mutate(
    y_c    = y - y_mean,
    lag_c  = lag_y - mean(lag_y, na.rm = TRUE),
    lisa_cluster = case_when(
      lisa_p <= 0.05 & y_c >= 0 & lag_c >= 0 ~ "High-High",
      lisa_p <= 0.05 & y_c <= 0 & lag_c <= 0 ~ "Low-Low",
      lisa_p <= 0.05 & y_c >= 0 & lag_c <= 0 ~ "High-Low",
      lisa_p <= 0.05 & y_c <= 0 & lag_c >= 0 ~ "Low-High",
      TRUE                                   ~ "Not significant"
    )
  )

table(jatim_data_utm$lisa_cluster)

## 
##        High-Low        Low-High Not significant 
##               1               2              35

## 13) PETA: Choropleth Y
tmap_mode("plot")

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

tm_shape(jatim_data_utm) +
  tm_polygons(
    "prevalensi_pneumonia_balita",
    style = "quantile",
    n = 5,
    palette = "Blues",
    title = "Prevalensi Pneumonia Balita\n(per 100.000 balita)"
  ) +
  tm_text(
    "Kabupaten/Kota",
    size = 0.5,
    auto.placement = TRUE
  ) +
  tm_layout(
    legend.outside = TRUE,
    main.title = "Peta Prevalensi Pneumonia Balita Jawa Timur",
    main.title.position = "center"
  )

## 
## ── 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') to
##   'tm_scale_intervals(<HERE>)'[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_text()`: migrate the layer options 'auto.placement' (rename to
## 'point.label') to 'options = opt_tm_text(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

## 14) PETA: LISA Cluster Map
tm_shape(jatim_data_utm) +
  tm_polygons(
    "lisa_cluster",
    palette = c(
      "High-High" = "red",
      "Low-Low" = "blue",
      "High-Low" = "orange",
      "Low-High" = "lightblue",
      "Not significant" = "grey80"
    ),
    title = "Klaster LISA (p ≤ 0.05)"
  ) +
  tm_text(
    "Kabupaten/Kota",
    size = 0.5,
    auto.placement = TRUE
  ) +
  tm_layout(
    legend.outside = TRUE,
    main.title = "Peta Klaster LISA Prevalensi Pneumonia Balita",
    main.title.position = "center"
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[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_text()`: migrate the layer options 'auto.placement' (rename to
## 'point.label') to 'options = opt_tm_text(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

## (Opsional) PETA: p-value LISA
tm_shape(jatim_data_utm) +
  tm_polygons(
    "lisa_p",
    style = "fixed",
    breaks = c(0, 0.01, 0.05, 0.1, 1),
    palette = "-RdYlGn",
    title = "p-value LISA"
  ) +
  tm_text(
    "Kabupaten/Kota",
    size = 0.5,
    auto.placement = TRUE
  ) +
  tm_layout(
    legend.outside = TRUE,
    main.title = "Peta Signifikansi LISA",
    main.title.position = "center"
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'[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_text()`: migrate the layer options 'auto.placement' (rename to
## 'point.label') to 'options = opt_tm_text(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`Multiple palettes called "rd_yl_gn" found: "brewer.rd_yl_gn", "matplotlib.rd_yl_gn". The first one, "brewer.rd_yl_gn", is returned.
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-RdYlGn" is named
## "rd_yl_gn" (in long format "brewer.rd_yl_gn")

4.4 Multikolinearitas

ols2 <- lm(
  prevalensi_pneumonia_balita ~
    `Kepadatan Penduduk` +
    `Prevalensi Balita Gizi Kurang` +
    `Cakupan Imunisasi`,
  data = jatim_data_utm
)

summary(ols2)

## 
## Call:
## lm(formula = prevalensi_pneumonia_balita ~ `Kepadatan Penduduk` + 
##     `Prevalensi Balita Gizi Kurang` + `Cakupan Imunisasi`, data = jatim_data_utm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2878.2 -1073.0  -219.6  1145.6  5594.8 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                     1604.5033  1610.3626   0.996   0.3261  
## `Kepadatan Penduduk`               0.2921     0.1249   2.340   0.0253 *
## `Prevalensi Balita Gizi Kurang`  -59.4644   101.3149  -0.587   0.5611  
## `Cakupan Imunisasi`               24.3892    23.5315   1.036   0.3073  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1716 on 34 degrees of freedom
## Multiple R-squared:  0.1941, Adjusted R-squared:  0.123 
## F-statistic: 2.729 on 3 and 34 DF,  p-value: 0.0591

# Pastikan ols2 sudah dibuat
stopifnot(exists("ols2"))

# VIF (untuk multikolinearitas antar X)
vif_vals <- car::vif(ols2)

vif_tbl <- data.frame(
  Variabel = names(vif_vals),
  VIF = as.numeric(vif_vals)
) %>%
  mutate(
    Kategori = case_when(
      VIF < 5 ~ "Rendah (<5)",
      VIF < 10 ~ "Sedang (5–10)",
      TRUE ~ "Tinggi (≥10)"
    )
  )

knitr::kable(vif_tbl, digits = 3,
             caption = "Uji multikolinearitas menggunakan Variance Inflation Factor (VIF).")

Uji multikolinearitas menggunakan Variance Inflation Factor (VIF).
Variabel	VIF	Kategori
`Kepadatan Penduduk`	1.034	Rendah (<5)
`Prevalensi Balita Gizi Kurang`	1.010	Rendah (<5)
`Cakupan Imunisasi`	1.024	Rendah (<5)

4.5 Uji Autokorelasi Spasial

4.6 B) Uji Autokorelasi Residual OLS

4.6.1 B1) Autokorelasi residual (Durbin–Watson, non-spasial)

Ini standar untuk cek autokorelasi residual “urutan” (time/sequence). Untuk data spasial, yang lebih relevan adalah Moran residual (B2).

library(lmtest)

# Durbin-Watson untuk residual OLS
dw <- lmtest::dwtest(ols2)
dw

## 
##  Durbin-Watson test
## 
## data:  ols2
## DW = 2.3271, p-value = 0.7879
## alternative hypothesis: true autocorrelation is greater than 0

library(spdep)

# Pastikan lw sudah dibuat
stopifnot(exists("lw"))

res_ols <- residuals(ols2)

moran_res <- moran.test(res_ols, lw, zero.policy = TRUE)
moran_res

## 
##  Moran I test under randomisation
## 
## data:  res_ols  
## weights: lw    
## 
## Moran I statistic standard deviate = -2.3656, p-value = 0.991
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.32606140       -0.02702703        0.01597882

# Tabel ringkas
res_tbl <- data.frame(
  Uji = "Moran's I Residual (OLS)",
  Statistik = unname(moran_res$estimate[["Moran I statistic"]]),
  P_value = moran_res$p.value
)

knitr::kable(res_tbl, digits = 6,
             caption = "Uji autokorelasi spasial pada residual model OLS menggunakan Moran's I.")

Uji autokorelasi spasial pada residual model OLS menggunakan Moran’s I.
Uji	Statistik	P_value
Moran’s I Residual (OLS)	-0.326061	0.991001

library(spdep)

lm_tests <- lm.LMtests(
  ols2, lw,
  test = c("LMlag","LMerr","RLMlag","RLMerr"),
  zero.policy = TRUE
)

## Please update scripts to use lm.RStests in place of lm.LMtests

lm_tests

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = prevalensi_pneumonia_balita ~ `Kepadatan Penduduk`
## + `Prevalensi Balita Gizi Kurang` + `Cakupan Imunisasi`, data =
## jatim_data_utm)
## test weights: listw
## 
## RSlag = 8.0246, df = 1, p-value = 0.004615
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = prevalensi_pneumonia_balita ~ `Kepadatan Penduduk`
## + `Prevalensi Balita Gizi Kurang` + `Cakupan Imunisasi`, data =
## jatim_data_utm)
## test weights: listw
## 
## RSerr = 5.6858, df = 1, p-value = 0.0171
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = prevalensi_pneumonia_balita ~ `Kepadatan Penduduk`
## + `Prevalensi Balita Gizi Kurang` + `Cakupan Imunisasi`, data =
## jatim_data_utm)
## test weights: listw
## 
## adjRSlag = 3.2699, df = 1, p-value = 0.07056
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = prevalensi_pneumonia_balita ~ `Kepadatan Penduduk`
## + `Prevalensi Balita Gizi Kurang` + `Cakupan Imunisasi`, data =
## jatim_data_utm)
## test weights: listw
## 
## adjRSerr = 0.93117, df = 1, p-value = 0.3346

ANALISIS SPASIAL KASUS PNEUMONIA PADA BALITA DI JAWA TIMUR TAHUN 2022

Stephani Julieta

2025-12-22

1 Pendahuluan

1.1 Latar Belakang

1.2 Identifikasi Masalah

1.3 Tujuan Penelitian

1.4 Batasan Penelitian

2 Tinjauan Pustaka

2.1 Dependensi Spasial

2.2 Autokorelasi Spasial

2.2.1 Rumus Moran’s I (Global Spasial Autocorrelation)

2.2.2 Rumus Geary’s C

2.2.3 Rumus LISA (Local Indicators Spatial Autocorrelation)

2.2.4 Scatterplot Moran & Hotspot

2.3 Model Ekonometrika Spasial

2.3.1 Spatial Autoregressive Model

2.3.2 Spatial Error Model (SEM)

2.3.3 IDW (Inverse Distance Weighting)

2.3.4 Spatial Durbin Model

2.3.5 Spatial Durbin Error Model

3 Metodologi Penelitian

3.1 Data dan Unit Spasial

3.2 Metode Analisis

3.3 Alur Kerja Penelitian

4 Hasil dan Pembahasan

4.1 Import Data

4.2 Peta Deskriptif dan Statistik Deskriptif

4.3 3) Statistik Deskriptif (yang kamu bilang belum ada)

4.3.1 3.1 Ringkasan numerik variabel (tabel rapi)

4.4 Multikolinearitas

4.5 Uji Autokorelasi Spasial

4.6 B) Uji Autokorelasi Residual OLS

4.6.1 B1) Autokorelasi residual (Durbin–Watson, non-spasial)

4.7 Estimasi Model

4.8 Pemilihan Model Terbaik

4.9 Uji Autokorelasi Spasial Model Terbaik

4.10 Peta Prediksi dan Residual

4.11 Ringkasan Output

5 Kesimpulan dan Saran

6 Daftar Pustaka

7 Lampiran