ANALISIS SPASIAL TINGKAT PENGANGGURAN TERBUKA PER KABUPATEN/KOTA DI PROVINSI JAWA BARAT TAHUN 2024

Disusun oleh: Muhammad Naufal Abdul Ghani (140610230063)

Dosen Pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.

Program Studi Statistika

Fakultas Matematika dan Ilmu Pengetahuan Alam

Universitas Padjadjaran

Jatinangor 2025

BAB 1. Pendahuluan

1.1. Latar Belakang

Pengangguran terbuka merupakan indikator kunci dalam analisis sosial-ekonomi suatu wilayah, yang mencerminkan ketidakmampuan tenaga kerja untuk memperoleh pekerjaan meskipun aktif mencari. Tingkat pengangguran terbuka yang tinggi sering kali terkait dengan masalah seperti kemiskinan, ketimpangan sosial, dan pertumbuhan ekonomi yang stagnan. Di Indonesia, data pengangguran terbuka bervariasi antarwilayah, dipengaruhi oleh faktor-faktor seperti tingkat pendidikan, akses infrastruktur, investasi, dan kebijakan pemerintah. Analisis spasial menjadi penting karena pengangguran tidak terdistribusi secara acak; pola spasialnya dapat menunjukkan kluster di daerah tertentu, seperti perkotaan versus pedesaan, yang memerlukan pendekatan pemodelan yang mempertimbangkan ketergantungan lokasi.

Untuk menganalisis dan memodelkan tingkat pengangguran terbuka serta faktor-faktornya, metode spasial seperti Spatial Autoregressive (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), dan Spatial Autocorrelation Model (SAC) dapat digunakan. Model-model ini memungkinkan pengujian hubungan spasial antara variabel dependen (pengangguran) dan independen (misalnya, tingkat pendidikan, PDB per kapita, atau jarak ke pusat ekonomi), dengan mempertimbangkan efek spillover antarwilayah. SAR fokus pada autokorelasi spasial dalam variabel dependen, SEM menangani kesalahan spasial, SDM menggabungkan efek langsung dan tidak langsung, SDEM mengintegrasikan elemen Durbin dengan kesalahan spasial, sedangkan SAC mengombinasikan autokorelasi dan kesalahan untuk model yang lebih komprehensif.

Selain pemodelan, interpolasi spasial seperti Inverse Distance Weighting (IDW) dan Kriging diperlukan untuk memperkirakan nilai pengangguran di lokasi yang tidak memiliki data langsung. IDW menggunakan bobot berbasis jarak untuk interpolasi sederhana, sementara Kriging memanfaatkan variogram untuk memperhitungkan variabilitas spasial dan memberikan estimasi yang lebih akurat dengan interval kepercayaan. Pendekatan ini akan membantu mengidentifikasi pola spasial pengangguran, memprediksi risiko di area baru, dan mendukung perumusan kebijakan seperti program pelatihan kerja atau investasi infrastruktur yang lebih tepat sasaran. Analisis ini diharapkan memberikan wawasan mendalam untuk mengurangi pengangguran melalui intervensi berbasis bukti spasial.

BAB 2. Tinjauan Pustaka

Penelitian mengenai analisis spasial fenomena sosial-ekonomi, seperti Tingkat Pengangguran Terbuka (TPT), telah banyak dilakukan untuk memahami pola persebaran dan faktor-faktor yang memengaruhi variabel tersebut secara geografis. Salah satu pendekatan yang sering digunakan dalam kajian spasial adalah permodelan spasial (spatial modeling) yang mampu menangkap ketergantungan antarwilayah serta efek lokasi terhadap variabel yang diteliti. Dalam studi yang diterbitkan pada Engineering, MAthematics and Computer Science Journal (EMACS), Mahmud (2021) mengembangkan permodelan spasial untuk menganalisis faktor-faktor yang memengaruhi suatu fenomena berdasarkan hubungan spasial antarwilayah. Studi tersebut menunjukkan pentingnya mempertimbangkan elemen spasial dalam analisis faktor untuk menangkap heterogenitas dan autokorelasi yang tidak dapat dijelaskan oleh model regresi klasik biasa.

Penelitian Mahmud (2021) dan kajian spasial lainnya memberikan landasan teoritis bahwa analisis fenomena sosial-ekonomi seperti TPT memerlukan model yang tidak hanya mempertimbangkan hubungan antarvariabel tetapi juga hubungan antarlokasi geografis untuk menghasilkan temuan yang komprehensif dan menggambarkan struktur spasial yang ada. Pendekatan semacam ini menjadi dasar dalam merancang metodologi penelitian spasial pada TPT di Jawa Barat tahun 2024, khususnya dalam penggunaan uji autokorelasi spasial serta interpolasi geostatistik sebagai bagian dari pemetaan dan analisis data spasial. [1]

BAB 3. Metode Penelitian

3.1. Data dan Sumber Data

Data diperoleh dari Badan Pusat Statistik Provinsi Jawa Barat dan Open Data Provinsi Jawa Barat, dengan variabel dependennya yaitu tingkat pengangguran terbuka per kabupaten/kota di Jawa Barat, dan untuk variabel independennya yaitu indeks pembangunan manusia, jumlah penduduk miskin, dan produk domestik regional bruto per kabupaten/kota di Jawa Barat. [2][3][4][5]

Variabel Keterangan
Y Tingkat Pengangguran Terbuka
X1 Indeks Pembangunan Manusia
X2 Jumlah Penduduk Miskin
X3 Produk Domestik Regional Bruto

3.2. Statistika Deskriptif

Statistika deskriptif merupakan metode analisis data yang digunakan untuk menggambarkan dan meringkas karakteristik utama dari data yang dikumpulkan tanpa melakukan penarikan kesimpulan secara inferensial. Tujuan utama statistika deskriptif adalah memberikan gambaran umum mengenai pola, kecenderungan, dan penyebaran data sehingga data lebih mudah dipahami.

Dalam penelitian ini, statistika deskriptif digunakan untuk menyajikan informasi mengenai data penelitian melalui ukuran-ukuran statistik, seperti:

  • Ukuran pemusatan, yang meliputi nilai rata-rata (mean), median, dan modus;
  • Ukuran penyebaran, seperti nilai minimum, maksimum, rentang (range), varians, dan simpangan baku (standar deviasi);
  • Penyajian data, yang dapat berupa tabel distribusi frekuensi, diagram batang, diagram lingkaran, atau grafik lainnya.

Melalui analisis statistika deskriptif, peneliti dapat memperoleh gambaran awal tentang kondisi data, mengidentifikasi pola tertentu, serta mendeteksi adanya data ekstrem (outlier) yang dapat memengaruhi analisis selanjutnya.

3.3. Peta Deskriptif

Peta deskriptif merupakan metode penyajian data dalam bentuk visual spasial yang bertujuan untuk menggambarkan persebaran fenomena atau variabel penelitian berdasarkan wilayah geografis tertentu. Metode ini digunakan untuk memudahkan interpretasi data yang memiliki keterkaitan dengan lokasi atau wilayah.

Dalam penelitian ini, peta deskriptif digunakan untuk menampilkan variasi nilai suatu variabel antar wilayah melalui simbol, warna, atau gradasi warna tertentu. Setiap wilayah pada peta diberi representasi visual sesuai dengan karakteristik data yang dimilikinya, sehingga perbedaan dan pola persebaran dapat terlihat secara jelas.

Peta deskriptif membantu peneliti dalam:

  • Mengidentifikasi pola spasial, seperti wilayah dengan nilai tinggi, sedang, atau rendah;
  • Menunjukkan ketimpangan atau kesenjangan antar wilayah;
  • Memberikan gambaran visual yang lebih informatif dan mudah dipahami dibandingkan penyajian data dalam bentuk tabel semata.

Dengan menggunakan peta deskriptif, hasil analisis menjadi lebih komunikatif dan dapat mendukung interpretasi temuan penelitian secara komprehensif.

3.4. Uji Autokorelasi Spasial

Uji autokorelasi spasial digunakan untuk mengetahui apakah terdapat hubungan atau keterkaitan antarwilayah berdasarkan lokasi geografisnya.

Autokorelasi spasial bisa bersifat:

  • Positif, jika wilayah yang berdekatan memiliki nilai yang mirip (tinggi–tinggi atau rendah–rendah).
  • Negatif, jika wilayah yang berdekatan memiliki nilai yang berlawanan (tinggi–rendah).
  • Netral (tidak ada autokorelasi), jika tidak ada pola spasial yang jelas. [4]

3.4.1. Moran’s I

Moran’s I merupakan ukuran global untuk mengetahui apakah terdapat pola spasial secara keseluruhan pada data. Nilai Moran’s I berkisar antara –1 hingga +1, di mana:

  • Nilai positif menunjukkan adanya autokorelasi spasial positif (wilayah yang berdekatan memiliki nilai serupa).
  • Nilai negatif menunjukkan autokorelasi spasial negatif (wilayah berdekatan memiliki nilai berbeda).
  • Nilai mendekati 0 menunjukkan tidak ada autokorelasi spasial. [5]

Rumus dari Moran’s I adalah,

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

Keterangan:

- \(I\) = indeks Moran’s I
- \(n\) = jumlah unit observasi (misalnya wilayah, kabupaten, atau kecamatan)
- \(x_i, x_j\) = nilai variabel pengamatan pada lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata dari seluruh nilai \(x\)
- \(w_{ij}\) = bobot spasial antara lokasi \(i\) dan \(j\)
- \(W = \sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}\) = jumlah total bobot spasial

3.4.2. Geary’s C

Geary’s C juga mengukur autokorelasi spasial global, tetapi lebih sensitif terhadap perbedaan lokal antarwilayah tetangga. [6]
Nilai Geary’s C berkisar antara 0 hingga 2, dengan interpretasi:

  • C<1: autokorelasi spasial positif,
  • C>1: autokorelasi spasial negatif,
  • C=1: tidak ada autokorelasi spasial.

Rumus dari Geary’s C adalah, 

\[ C = \frac{(n - 1)}{2W} \cdot \frac{\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}(x_i - x_j)^2} {\displaystyle\sum_{i=1}^{n}(x_i - \bar{x})^2} \]

Keterangan:

- \(C\) = indeks Geary’s C
- \(n\) = jumlah unit observasi
- \(x_i, x_j\) = nilai variabel pada lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata dari seluruh nilai \(x\)
- \(w_{ij}\) = bobot spasial antara lokasi \(i\) dan \(j\)
- \(W = \sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}\) = jumlah total bobot spasial

3.4.3. LISA (Local Indicator of Spatial Autocorelation)

LISA merupakan pengukuran autokorelasi spasial lokal, yang digunakan untuk mengetahui wilayah mana saja yang memiliki pola spasial signifikan. [7]

Berbeda dengan Moran’s I dan Geary’s C yang bersifat global, LISA fokus pada identifikasi klaster di tingkat wilayah.

LISA dapat mengidentifikasi empat jenis pola spasial:

  • High-High (HH): wilayah dengan nilai tinggi dikelilingi wilayah bernilai tinggi.
  • Low-Low (LL): wilayah dengan nilai rendah dikelilingi wilayah bernilai rendah.
  • High-Low (HL): wilayah bernilai tinggi dikelilingi wilayah bernilai rendah.
  • Low-High (LH): wilayah bernilai rendah dikelilingi wilayah bernilai tinggi.

Melalui analisis LISA, peneliti dapat memetakan wilayah yang menjadi hotspot atau coldspot, yang sangat berguna dalam kebijakan daerah berbasis spasial, seperti pemerataan gender, pendidikan, atau kesehatan.

3.5. Model Spasial

3.5.1. OLS

Model OLS (Ordinary Least Squares) digunakan ketika hubungan antara variabel dependen dan variabel independen diasumsikan tidak dipengaruhi oleh efek spasial atau ketergantungan antarwilayah. Artinya, nilai variabel dependen di suatu wilayah hanya dipengaruhi oleh variabel independen di wilayah tersebut, tanpa mempertimbangkan pengaruh wilayah sekitarnya.

Persamaan dari OLS adalah:

\[ \mathbf{y} =\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

- \(\mathbf{y}\): vektor variabel dependen
- \(\mathbf{X}\): matriks variabel independen
- \(\boldsymbol{\beta}\): vektor koefisien regresi
- \(\boldsymbol{\varepsilon}\): vektor error acak

3.5.2. SAR

Model SAR digunakan ketika nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Artinya, terdapat spatial lag dari variabel dependen.

Persamaan dari SAR adalah:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

Keterangan:
- \(\mathbf{y}\): vektor variabel dependen
- \(\mathbf{X}\): matriks variabel independen
- \(\boldsymbol{\beta}\): vektor koefisien regresi
- \(\rho\): parameter autokorelasi spasial (spatial lag coefficient)
- \(\mathbf{W}\): matriks bobot spasial (spatial weight matrix)
- \(\boldsymbol{\varepsilon}\): vektor error acak

Model SAR mengasumsikan bahwa nilai variabel dependen di suatu lokasi dipengaruhi oleh nilai-nilai variabel dependen di lokasi-lokasi tetangganya.

3.5.3. SEM

Model SEM digunakan ketika pengaruh spasial tidak muncul pada nilai variabel dependen secara langsung, tetapi melalui komponen error (gangguan) yang saling berkorelasi antarwilayah.

Persamaan dari SEM adalah:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} + \mathbf{u} \]

Keterangan:
- \(\lambda\): parameter autokorelasi spasial pada error
- \(\mathbf{u}\): komponen error acak yang bebas spasial

Model SEM digunakan ketika ketergantungan spasial tidak terdapat pada \(y\), melainkan pada komponen error.

3.5.4. SDM

Model SDM adalah model yang lebih umum yang memasukkan lag spasial dari variabel dependen WY dan juga lag spasial dari variabel-variabel independen WX. Model ini dapat menangkap efek limpahan yang berasal dari variabel dependen (interaksi endogen) dan variabel independen (interaksi eksogen). SDM mampu mengestimasi efek langsung (dampak perubahan variabel independen di suatu wilayah terhadap variabel dependen di wilayah itu sendiri) dan efek tidak langsung/limpahan (dampak perubahan di suatu wilayah terhadap wilayah lain). Persamaan model SDM adalah:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \]

Keterangan:
- \(\mathbf{W}\mathbf{X}\): lag spasial dari variabel independen
- \(\boldsymbol{\theta}\): koefisien efek tidak langsung (spillover effects)

Model SDM menggabungkan lag spasial pada variabel dependen (\(\mathbf{W}\mathbf{y}\)) dan lag pada variabel independen (\(\mathbf{W}\mathbf{X}\)),
sehingga menangkap pengaruh antarwilayah baik secara langsung maupun tidak langsung.

3.5.5. SDEM

Model SDEM adalah variasi dari model SEM yang juga memasukkan lag spasial dari variabel independen WX, namun tidak memasukkan lag spasial dari variabel dependen WY. Model ini cocok untuk situasi di mana autokorelasi spasial diyakini berada pada galat, namun tetap ada efek limpahan yang berasal dari variabel-variabel independen di wilayah tetangga. Persamaan model SDEM adalah:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} + \mathbf{u} \]

Keterangan:
- Model ini menggabungkan lag spasial pada variabel independen (\(\mathbf{W}\mathbf{X}\)) dengan autokorelasi error (\(\lambda \mathbf{W}\boldsymbol{\varepsilon}\)).
- Tidak ada lag pada variabel dependen (\(\mathbf{W}\mathbf{y}\)).

SDEM digunakan ketika efek spasial muncul pada variabel independen dan error, bukan pada \(y\).

3.5.6. SAC

Model SAC, juga dikenal sebagai General Spatial Model atau model Kelejian-Prucha, adalah model yang paling komprehensif karena memasukkan kedua jenis dependensi spasial secara bersamaan: lag spasial pada variabel dependen WY dan autokorelasi pada galat Wu. Model ini digunakan ketika ada alasan teoretis untuk meyakini bahwa kedua proses spasial (substantif dan residual) terjadi secara simultan. Persamaan model SAC adalah:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} + \mathbf{u} \]

Keterangan:
- Menggabungkan dua bentuk dependensi spasial:
- Lag spasial pada \(y\) seperti pada SAR
- Autokorelasi spasial pada error seperti pada SEM

Model SAC adalah bentuk paling umum dari model spasial, digunakan bila terdapat efek spasial baik pada variabel dependen maupun pada error.

3.6. Interpolasi Spasial

Interpolasi spasial merupakan metode analisis yang digunakan untuk memperkirakan nilai suatu variabel pada lokasi yang tidak memiliki data pengamatan secara langsung berdasarkan nilai variabel pada lokasi-lokasi terdekat. Dalam penelitian ini, interpolasi spasial digunakan untuk mengestimasi Tingkat Pengangguran Terbuka (TPT) di wilayah Jawa Barat secara kontinu pada seluruh area penelitian. Metode interpolasi yang digunakan meliputi Inverse Distance Weighting (IDW) dan Kriging, yang masing-masing memiliki pendekatan dan asumsi yang berbeda dalam memodelkan keterkaitan spasial antar lokasi.

3.6.1. IDW

Inverse Distance Weighting (IDW) merupakan metode interpolasi deterministik yang mengasumsikan bahwa pengaruh suatu titik pengamatan terhadap lokasi yang akan diprediksi berbanding terbalik dengan jaraknya. Semakin dekat jarak suatu titik pengamatan, semakin besar bobot yang diberikan terhadap nilai prediksi.

3.6.2. Kriging

Kriging merupakan metode interpolasi geostatistik yang mempertimbangkan struktur ketergantungan spasial melalui model variogram. Berbeda dengan IDW, Kriging tidak hanya memperhitungkan jarak antar titik, tetapi juga pola variabilitas spasial dari data.

BAB 4. Hasil dan Pembahasan

Load Library

library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
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(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(spatialreg)
## 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(tmap)
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.2
## corrplot 0.95 loaded
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.2
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
library(terra)
## terra 1.8.70

Persiapan Data

# === Membaca Data
data <- read.csv("D:/M Naufal/Spasial/Data UAS Spasial 1.csv", sep = ";", dec = ".")
jabar <- readRDS("D:/M Naufal/METOPEN/jabar_shp.rds")

# === Persiapan Data
dataku <- data.frame(
  P = toupper(data$Kabupaten.Kota),
  Y = data$TPT,
  X1 = data$IPM,
  X2 = data$JPM,
  X3 = data$PDRB
)

# Filter Jawa Barat dan samakan nama
if("admin1Name_en" %in% names(jabar)) {
  jabar <- jabar %>% filter(admin1Name_en == "Jawa Barat")
}

jabar$NAME_2 <- toupper(jabar$admin2Name_en)

# Tambahkan prefix KABUPATEN/KOTA
if(!any(grepl("KABUPATEN|KOTA", jabar$NAME_2))) {
  kota_list <- c("BANDUNG", "BOGOR", "SUKABUMI", "CIREBON", "BEKASI", 
                 "DEPOK", "CIMAHI", "TASIKMALAYA", "BANJAR")
  jabar$NAME_2 <- ifelse(jabar$NAME_2 %in% kota_list,
                         paste("KOTA", jabar$NAME_2),
                         paste("KABUPATEN", jabar$NAME_2))
}

# Merge data
jabar_merged <- left_join(jabar, dataku, by = c("NAME_2" = "P")) %>%
  filter(!is.na(Y))

4.1. Statistika Deskriptif

cat("\n")
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
cat("1. STATISTIKA DESKRIPTIF\n")
## 1. STATISTIKA DESKRIPTIF
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
cat("\nTingkat Pengangguran Terbuka (Y):\n")
## 
## Tingkat Pengangguran Terbuka (Y):
print(summary(jabar_merged$Y))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.580   6.205   6.730   6.525   7.590   8.970
cat("Standar Deviasi:", sd(jabar_merged$Y, na.rm = TRUE), "\n")
## Standar Deviasi: 1.706509
cat("\nIndeks Pembangunan Manusia (X1):\n")
## 
## Indeks Pembangunan Manusia (X1):
print(summary(jabar_merged$X1))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   68.89   71.20   73.82   74.68   77.25   83.75
cat("Standar Deviasi:", sd(jabar_merged$X1, na.rm = TRUE), "\n")
## Standar Deviasi: 4.333326
cat("\nJumlah Penduduk Miskin (X2):\n")
## 
## Jumlah Penduduk Miskin (X2):
print(summary(jabar_merged$X2))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    11.2    75.3   131.8   142.5   196.2   446.8
cat("Standar Deviasi:", sd(jabar_merged$X2, na.rm = TRUE), "\n")
## Standar Deviasi: 97.05356
cat("\nProduk Domestik Regional Bruto (X3):\n")
## 
## Produk Domestik Regional Bruto (X3):
print(summary(jabar_merged$X3))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5549   43330   61236  105053   99869  421324
cat("Standar Deviasi:", sd(jabar_merged$X3, na.rm = TRUE), "\n")
## Standar Deviasi: 112708.7

4.2. Peta Deskriptif

cat("\n--- Peta Sebaran TPT ---\n")
## 
## --- Peta Sebaran TPT ---
peta_sebaran <- tm_shape(jabar_merged) +
  tm_polygons(
    "Y",
    title = "Tingkat Pengangguran\nTerbuka (%)",
    palette = "YlGn",
    style = "quantile",
    n = 5,
    border.col = "black",
    border.lwd = 0.5
  ) +
  tm_text("NAME_2", size = 0.4, remove.overlap = TRUE) +
  tm_layout(
    main.title = "Sebaran TPT di Jawa Barat 2024",
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── 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()`: 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>)'
## [tm_polygons()] Argument `border.lwd` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_sebaran)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.

Peta sebaran Tingkat Pengangguran Terbuka (TPT) di Provinsi Jawa Barat tahun 2024 menunjukkan adanya variasi nilai TPT yang cukup jelas antar kabupaten/kota. Wilayah dengan kategori TPT tinggi hingga sangat tinggi umumnya terkonsentrasi di bagian utara dan barat Jawa Barat, khususnya pada wilayah perkotaan dan kawasan industri. Sebaliknya, wilayah dengan TPT rendah hingga sedang lebih banyak ditemukan di bagian selatan dan timur, yang didominasi oleh wilayah dengan karakteristik ekonomi berbasis pertanian dan jasa lokal. Pola ini mengindikasikan bahwa perbedaan struktur ekonomi dan tingkat urbanisasi berperan penting dalam membentuk variasi TPT antarwilayah di Jawa Barat.

4.3. Uji Autokorelasi Spasial

Membentuk Matriks Bobot

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

4.3.1. Moran’s I

# === Moran's I
cat("\n--- MORAN'S I GLOBAL ---\n")
## 
## --- MORAN'S I GLOBAL ---
moran_test <- moran.test(jabar_merged$Y, lw, zero.policy = TRUE)
print(moran_test)
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$Y  
## weights: lw    
## 
## Moran I statistic standard deviate = 4.0655, p-value = 2.397e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.51087888       -0.03846154        0.01825812
cat("\nInterpretasi:")
## 
## Interpretasi:
if(moran_test$estimate[1] > 0 & moran_test$p.value < 0.05) {
  cat("\n✓ Autokorelasi spasial POSITIF signifikan")
  cat("\n  (Wilayah dengan TPT tinggi bertetangga dengan TPT tinggi)\n")
} else if(moran_test$estimate[1] < 0 & moran_test$p.value < 0.05) {
  cat("\n✓ Autokorelasi spasial NEGATIF signifikan\n")
} else {
  cat("\n✗ Tidak ada autokorelasi spasial signifikan\n")
}
## 
## ✓ Autokorelasi spasial POSITIF signifikan
##   (Wilayah dengan TPT tinggi bertetangga dengan TPT tinggi)

Berdasarkan hasil pengujian Moran’s I Global, diperoleh nilai Moran’s I sebesar 0,5109 dengan standard deviate sebesar 4,0655 dan p-value sebesar 2,397 × 10⁻⁵. Nilai Moran’s I yang positif dan signifikan secara statistik (p-value < 0,05) menunjukkan adanya autokorelasi spasial positif yang signifikan pada Tingkat Pengangguran Terbuka (TPT) di Provinsi Jawa Barat. Hal ini mengindikasikan bahwa wilayah-wilayah dengan nilai TPT tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki TPT tinggi, demikian pula wilayah dengan TPT rendah cenderung berdekatan dengan wilayah bernilai rendah. Dengan demikian, pola persebaran TPT antarwilayah tidak bersifat acak, melainkan membentuk pola pengelompokan spasial.

4.3.2. Geary’s C

cat("\n--- GEARY'S C ---\n")
## 
## --- GEARY'S C ---
geary_test <- geary.test(jabar_merged$Y, lw, zero.policy = TRUE)
print(geary_test)
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$Y 
## weights: lw   
## 
## Geary C statistic standard deviate = 2.9199, p-value = 0.001751
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.46538345        1.00000000        0.03352364
cat("\nInterpretasi:")
## 
## Interpretasi:
if(geary_test$estimate[1] < 1 & geary_test$p.value < 0.05) {
  cat("\n✓ Autokorelasi spasial POSITIF signifikan (C < 1)")
  cat("\n  (Nilai mirip cenderung mengelompok)\n")
} else if(geary_test$estimate[1] > 1 & geary_test$p.value < 0.05) {
  cat("\n✓ Autokorelasi spasial NEGATIF signifikan (C > 1)\n")
} else {
  cat("\n✗ Tidak ada autokorelasi spasial signifikan\n")
}
## 
## ✓ Autokorelasi spasial POSITIF signifikan (C < 1)
##   (Nilai mirip cenderung mengelompok)

Hasil pengujian Geary’s C menunjukkan nilai statistik sebesar 0,4654, yang lebih kecil dari nilai ekspektasinya yaitu 1, dengan standard deviate sebesar 2,9199 dan p-value sebesar 0,001751. Nilai Geary’s C yang lebih kecil dari 1 dan signifikan pada taraf 5% menandakan adanya autokorelasi spasial positif yang signifikan. Temuan ini menunjukkan bahwa wilayah-wilayah yang saling bertetangga cenderung memiliki nilai TPT yang relatif mirip, sehingga mengindikasikan adanya kecenderungan pengelompokan nilai TPT pada wilayah-wilayah yang berdekatan.

4.3.3. LISA (Local Indicator of Spatial Autocorelation)

cat("\n--- LISA (Local Indicators of Spatial Association) ---\n")
## 
## --- LISA (Local Indicators of Spatial Association) ---
lisa <- localmoran(jabar_merged$Y, lw, zero.policy = TRUE)

jabar_merged$Ii <- lisa[, 1]
jabar_merged$Pvalue <- lisa[, 5]

# Klasifikasi cluster
mean_y <- mean(jabar_merged$Y)
lag_y <- lag.listw(lw, jabar_merged$Y, zero.policy = TRUE)

jabar_merged$cluster <- "Tidak Signifikan"
sig_level <- 0.05

jabar_merged$cluster[jabar_merged$Pvalue <= sig_level] <- "Signifikan"

# Untuk detail analisis tetap simpan tipe cluster
jabar_merged$cluster_type <- "Not Significant"
jabar_merged$cluster_type[jabar_merged$Y >= mean_y & lag_y >= mean_y & 
                            jabar_merged$Pvalue <= sig_level] <- "High-High"
jabar_merged$cluster_type[jabar_merged$Y < mean_y & lag_y < mean_y & 
                            jabar_merged$Pvalue <= sig_level] <- "Low-Low"
jabar_merged$cluster_type[jabar_merged$Y >= mean_y & lag_y < mean_y & 
                            jabar_merged$Pvalue <= sig_level] <- "High-Low"
jabar_merged$cluster_type[jabar_merged$Y < mean_y & lag_y >= mean_y & 
                            jabar_merged$Pvalue <= sig_level] <- "Low-High"

cat("\nDistribusi Signifikansi LISA:\n")
## 
## Distribusi Signifikansi LISA:
print(table(jabar_merged$cluster))
## 
##       Signifikan Tidak Signifikan 
##                4               23
cat("\n--- Peta Cluster LISA ---\n")
## 
## --- Peta Cluster LISA ---
peta_lisa <- tm_shape(jabar_merged) +
  tm_fill(
    "cluster",
    palette = c("Signifikan" = "red",
                "Tidak Signifikan" = "white"),
    title = "Indikator Local\nAutocorrelation"
  ) +
  tm_borders(col = "black", lwd = 0.5) +
  tm_text("NAME_2", size = 0.4, remove.overlap = TRUE) +
  tm_layout(
    main.title = "Pemetaan Hasil Pengujian\nSignifikansi Local Indicator of Spatial Autocorrelation (LISA)",
    main.title.size = 0.9,
    main.title.position = "center",
    legend.outside = TRUE,
    legend.outside.position = "right",
    legend.title.size = 0.9,
    legend.text.size = 0.8,
    frame = TRUE
  ) +
  tm_credits("Signifikansi: α = 0.05", position = c("left", "bottom"), size = 0.7)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_fill()`: 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_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_lisa)

Peta signifikansi Local Indicator of Spatial Autocorrelation (LISA) menunjukkan bahwa tidak seluruh wilayah di Jawa Barat memiliki autokorelasi spasial lokal yang signifikan. Wilayah yang ditandai dengan warna merah merupakan wilayah yang memiliki autokorelasi spasial lokal signifikan pada taraf signifikansi 5% (α = 0,05). Wilayah-wilayah signifikan tersebut terkonsentrasi di bagian selatan Jawa Barat, yang menunjukkan adanya keterkaitan kuat antara nilai TPT suatu wilayah dengan wilayah-wilayah di sekitarnya. Hal ini mengindikasikan adanya klaster spasial lokal, di mana wilayah dengan nilai TPT yang relatif serupa cenderung saling berdekatan. Sementara itu, wilayah yang tidak signifikan menunjukkan bahwa nilai TPT di wilayah tersebut tidak memiliki hubungan spasial lokal yang cukup kuat dengan wilayah sekitarnya, meskipun secara global masih terdapat autokorelasi spasial.

4.4. Model Spasial

4.4.1. OLS

cat("\n--- MODEL OLS (Baseline) ---\n")
## 
## --- MODEL OLS (Baseline) ---
model_ols <- lm(Y ~ X1 + X2 + X3, data = jabar_merged)
print(summary(model_ols))
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = jabar_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4777 -0.8545  0.1731  1.0099  2.0776 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -9.661e+00  7.490e+00  -1.290   0.2100  
## X1           2.042e-01  9.721e-02   2.100   0.0469 *
## X2           5.208e-03  4.739e-03   1.099   0.2831  
## X3           1.869e-06  3.808e-06   0.491   0.6282  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.512 on 23 degrees of freedom
## Multiple R-squared:  0.306,  Adjusted R-squared:  0.2155 
## F-statistic:  3.38 on 3 and 23 DF,  p-value: 0.0355

4.4.2. SAR

cat("\n--- MODEL SAR (Spatial Lag Model) ---\n")
## 
## --- MODEL SAR (Spatial Lag Model) ---
model_sar <- lagsarlm(Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, zero.policy = TRUE)
print(summary(model_sar))
## 
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.88777 -0.67077  0.23860  0.75826  2.45690 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -6.7168e+00  6.1706e+00 -1.0885   0.2764
## X1           1.3329e-01  8.2452e-02  1.6165   0.1060
## X2           2.5705e-03  3.9264e-03  0.6547   0.5127
## X3           8.0390e-07  3.1338e-06  0.2565   0.7975
## 
## Rho: 0.45353, LR test value: 4.6911, p-value: 0.030319
## Asymptotic standard error: 0.18156
##     z-value: 2.4979, p-value: 0.012492
## Wald statistic: 6.2397, p-value: 0.012492
## 
## Log likelihood: -44.9551 for lag model
## ML residual variance (sigma squared): 1.5448, (sigma: 1.2429)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 101.91, (AIC for lm: 104.6)
## LM test for residual autocorrelation
## test value: 0.017304, p-value: 0.89534

4.4.3. SEM

cat("\n--- MODEL SEM (Spatial Error Model) ---\n")
## 
## --- MODEL SEM (Spatial Error Model) ---
model_sem <- errorsarlm(Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, zero.policy = TRUE)
print(summary(model_sem))
## 
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93072 -0.76824  0.24650  0.86037  2.43362 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -5.5947e+00  7.2109e+00 -0.7759  0.43782
## X1           1.6065e-01  9.3273e-02  1.7224  0.08499
## X2           2.4899e-03  3.8346e-03  0.6493  0.51613
## X3          -1.0580e-07  3.2044e-06 -0.0330  0.97366
## 
## Lambda: 0.51102, LR test value: 4.0513, p-value: 0.044138
## Asymptotic standard error: 0.18006
##     z-value: 2.8381, p-value: 0.0045387
## Wald statistic: 8.0547, p-value: 0.0045387
## 
## Log likelihood: -45.27501 for error model
## ML residual variance (sigma squared): 1.5544, (sigma: 1.2468)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 102.55, (AIC for lm: 104.6)

4.4.4. SDM

cat("\n--- MODEL SDM (Spatial Durbin Model) ---\n")
## 
## --- MODEL SDM (Spatial Durbin Model) ---
model_sdm <- lagsarlm(Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
                      type = "mixed", zero.policy = TRUE)
print(summary(model_sdm))
## 
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
##     type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.8709027 -0.6280628 -0.0033014  0.7922163  2.7342228 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept)  1.9018e+01  1.3368e+01  1.4226  0.15485
## X1           1.7657e-01  9.4750e-02  1.8635  0.06239
## X2           4.4948e-03  3.8795e-03  1.1586  0.24663
## X3           3.2854e-07  2.9593e-06  0.1110  0.91160
## lag.X1      -3.9203e-01  2.0969e-01 -1.8696  0.06154
## lag.X2      -6.4356e-03  6.3862e-03 -1.0077  0.31358
## lag.X3       1.5604e-05  7.4756e-06  2.0873  0.03686
## 
## Rho: 0.29353, LR test value: 1.6656, p-value: 0.19685
## Asymptotic standard error: 0.21125
##     z-value: 1.3895, p-value: 0.16467
## Wald statistic: 1.9308, p-value: 0.16467
## 
## Log likelihood: -42.12734 for mixed model
## ML residual variance (sigma squared): 1.2971, (sigma: 1.1389)
## Number of observations: 27 
## Number of parameters estimated: 9 
## AIC: 102.25, (AIC for lm: 101.92)
## LM test for residual autocorrelation
## test value: 0.11183, p-value: 0.73807

4.4.5. SDEM

cat("\n--- MODEL SDEM (Spatial Durbin Error Model) ---\n")
## 
## --- MODEL SDEM (Spatial Durbin Error Model) ---
model_sdem <- errorsarlm(Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
                         etype = "emixed", zero.policy = TRUE)
print(summary(model_sdem))
## 
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
##     etype = "emixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.92036 -0.60092  0.06771  0.82892  2.67844 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept)  1.5928e+01  1.4016e+01  1.1364  0.25577
## X1           1.6202e-01  9.1577e-02  1.7692  0.07686
## X2           4.4338e-03  3.9647e-03  1.1183  0.26343
## X3           7.7457e-07  2.9353e-06  0.2639  0.79187
## lag.X1      -3.1480e-01  1.9661e-01 -1.6011  0.10935
## lag.X2      -4.9294e-03  6.0805e-03 -0.8107  0.41755
## lag.X3       1.6494e-05  7.4214e-06  2.2225  0.02625
## 
## Lambda: 0.31145, LR test value: 1.5111, p-value: 0.21898
## Asymptotic standard error: 0.21525
##     z-value: 1.4469, p-value: 0.14792
## Wald statistic: 2.0935, p-value: 0.14792
## 
## Log likelihood: -42.20459 for error model
## ML residual variance (sigma squared): 1.3006, (sigma: 1.1405)
## Number of observations: 27 
## Number of parameters estimated: 9 
## AIC: 102.41, (AIC for lm: 101.92)

4.4.6. SAC

cat("\n--- MODEL SAC (Spatial Autoregressive Combined) ---\n")
## 
## --- MODEL SAC (Spatial Autoregressive Combined) ---
model_sac <- sacsarlm(Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, zero.policy = TRUE)
print(summary(model_sac))
## 
## Call:sacsarlm(formula = Y ~ X1 + X2 + X3, data = jabar_merged, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.87025 -0.69967  0.25202  0.77599  2.46817 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -6.8207e+00  6.3498e+00 -1.0742   0.2827
## X1           1.3791e-01  9.7661e-02  1.4122   0.1579
## X2           2.5765e-03  4.0981e-03  0.6287   0.5295
## X3           6.8894e-07  3.1735e-06  0.2171   0.8281
## 
## Rho: 0.41988
## Asymptotic standard error: 0.51903
##     z-value: 0.80897, p-value: 0.41853
## Lambda: 0.062609
## Asymptotic standard error: 0.6654
##     z-value: 0.094092, p-value: 0.92504
## 
## LR test value: 4.7041, p-value: 0.095174
## 
## Log likelihood: -44.9486 for sac model
## ML residual variance (sigma squared): 1.5564, (sigma: 1.2475)
## Number of observations: 27 
## Number of parameters estimated: 7 
## AIC: 103.9, (AIC for lm: 104.6)

4.4.7. Pemilihan Model Terbaik

comparison <- data.frame(
  Model = c("OLS", "SAR", "SEM", "SDM", "SDEM", "SAC"),
  AIC = c(AIC(model_ols), AIC(model_sar), AIC(model_sem), 
          AIC(model_sdm), AIC(model_sdem), AIC(model_sac)),
  LogLik = c(logLik(model_ols), logLik(model_sar), logLik(model_sem),
             logLik(model_sdm), logLik(model_sdem), logLik(model_sac))
)

comparison <- comparison[order(comparison$AIC), ]
print(comparison)
##   Model      AIC    LogLik
## 2   SAR 101.9102 -44.95510
## 4   SDM 102.2547 -42.12734
## 5  SDEM 102.4092 -42.20459
## 3   SEM 102.5500 -45.27501
## 6   SAC 103.8972 -44.94860
## 1   OLS 104.6013 -47.30065

Model terbaik berdasarkan AIC

# Tentukan model terbaik berdasarkan AIC
best_model_name <- comparison$Model[1]
cat("\nModel terbaik:", best_model_name, "\n")
## 
## Model terbaik: SAR
# Pilih model sesuai nama
best_model <- switch(best_model_name,
                     "OLS" = model_ols,
                     "SAR" = model_sar,
                     "SEM" = model_sem,
                     "SDM" = model_sdm,
                     "SDEM" = model_sdem,
                     "SAC" = model_sac)

# Ambil nilai prediksi dan residual
if(best_model_name == "OLS") {
  jabar_merged$predicted <- fitted(best_model)
  jabar_merged$residual <- residuals(best_model)
} else {
  jabar_merged$predicted <- fitted(best_model)
  jabar_merged$residual <- residuals(best_model)
}
## This method assumes the response is known - see manual page
# Statistik prediksi dan residual
cat("\nStatistik Prediksi:\n")
## 
## Statistik Prediksi:
print(summary(jabar_merged$predicted))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.468   5.803   6.436   6.525   7.233   8.335
cat("Standar Deviasi:", sd(jabar_merged$predicted, na.rm = TRUE), "\n")
## Standar Deviasi: 1.051296
cat("\nStatistik Residual:\n")
## 
## Statistik Residual:
print(summary(jabar_merged$residual))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.8878 -0.6708  0.2386  0.0000  0.7583  2.4569
cat("Standar Deviasi:", sd(jabar_merged$residual, na.rm = TRUE), "\n")
## Standar Deviasi: 1.266575

4.4.8. Peta Prediksi Model Terbaik

# --- Peta Prediksi ---
cat("\n--- Peta Prediksi Model", best_model_name, "---\n")
## 
## --- Peta Prediksi Model SAR ---
peta_prediksi <- tm_shape(jabar_merged) +
  tm_polygons(
    "predicted",
    title = "TPT Prediksi (%)",
    palette = "YlGn",
    style = "quantile",
    n = 5,
    border.col = "black",
    border.lwd = 0.5
  ) +
  tm_text("NAME_2", size = 0.4, remove.overlap = TRUE) +
  tm_layout(
    main.title = paste("Peta Prediksi TPT - Model", best_model_name),
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── 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>)'
## [tm_polygons()] Argument `border.lwd` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_prediksi)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.

Peta prediksi Tingkat Pengangguran Terbuka (TPT) hasil Model Spatial Autoregressive (SAR) menunjukkan adanya pola spasial yang jelas antar kabupaten/kota di Jawa Barat. Wilayah dengan nilai TPT prediksi tinggi cenderung terkonsentrasi di kawasan perkotaan dan wilayah penyangga metropolitan seperti Bogor, Depok, Bekasi, dan Karawang. Kondisi ini mencerminkan tekanan pasar tenaga kerja yang lebih besar akibat urbanisasi, arus migrasi tenaga kerja, serta ketidakseimbangan antara penawaran dan permintaan kerja. Sebaliknya, wilayah bagian selatan dan timur Jawa Barat seperti Tasikmalaya, Ciamis, Pangandaran, dan Kuningan menunjukkan nilai TPT prediksi yang relatif lebih rendah, yang mengindikasikan peran sektor pertanian dan informal yang masih dominan serta tingkat urbanisasi yang lebih rendah.

4.4.9. Peta Residual Model Terbaik

# --- Peta Residual ---
cat("\n--- Peta Residual Model", best_model_name, "---\n")
## 
## --- Peta Residual Model SAR ---
peta_residual <- tm_shape(jabar_merged) +
  tm_polygons(
    "residual",
    title = "Residual",
    palette = "-RdBu",
    style = "jenks",
    n = 5,
    border.col = "black",
    border.lwd = 0.5,
    midpoint = 0
  ) +
  tm_text("NAME_2", size = 0.4, remove.overlap = TRUE) +
  tm_layout(
    main.title = paste("Peta Residual Model", best_model_name),
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'midpoint', 'palette' (rename to
##   'values') 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()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [tm_polygons()] Argument `border.lwd` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_residual)
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
## "rd_bu" (in long format "brewer.rd_bu")

Peta residual Model SAR memperlihatkan perbedaan antara nilai TPT aktual dan hasil prediksi model di masing-masing wilayah. Residual positif menunjukkan bahwa model cenderung meremehkan tingkat pengangguran yang sebenarnya, sedangkan residual negatif mengindikasikan bahwa model melebihkan prediksi TPT. Kabupaten Kuningan tampak memiliki residual positif yang cukup besar, yang menandakan adanya faktor-faktor lokal spesifik yang belum sepenuhnya terakomodasi dalam model, seperti struktur lapangan kerja atau karakteristik ekonomi daerah. Sementara itu, beberapa wilayah di bagian selatan Jawa Barat menunjukkan residual negatif, yang mengindikasikan bahwa kondisi pasar kerja aktual di wilayah tersebut relatif lebih baik dibandingkan dengan hasil prediksi model.

4.5. Interpolasi Spasial

4.5.1. Data Titik

jabar_merged <- st_transform(jabar_merged, crs = 32748)  # UTM Zona 48S (Jawa Barat)

# Hitung centroid tiap kab/kota
jabar_centroid <- st_centroid(jabar_merged)
## Warning: st_centroid assumes attributes are constant over geometries
# Ambil koordinat X dan Y
coords <- st_coordinates(jabar_centroid)
jabar_centroid$X <- coords[,1]
jabar_centroid$Y_coord <- coords[,2]

# Cek hasil
head(jabar_centroid[, c("NAME_2", "X", "Y_coord", "Y")])
## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 695438.7 ymin: 9192692 xmax: 878677.1 ymax: 9312518
## Projected CRS: WGS 84 / UTM zone 48S
##          NAME_2        X Y_coord    Y                    Shape
## 1       BANDUNG 788412.3 9214389 6.36 POINT (788412.3 9214389)
## 2 BANDUNG BARAT 766874.0 9236958 6.70   POINT (766874 9236958)
## 3        BEKASI 734650.5 9312518 8.82 POINT (734650.5 9312518)
## 4         BOGOR 695438.7 9274546 7.34 POINT (695438.7 9274546)
## 5        CIAMIS 878677.1 9192692 3.37 POINT (878677.1 9192692)
## 6       CIANJUR 738313.0 9210923 5.99   POINT (738313 9210923)

4.5.2. Matriks Bobot Spasial (Data Titik)

# Koordinat centroid
coords_centroid <- st_coordinates(jabar_centroid)

# K-nearest neighbors (misal k = 4)
nb_knn <- knn2nb(knearneigh(coords_centroid, k = 4))

# Listw
lw_knn <- nb2listw(nb_knn, style = "W")

# Moran's I berbasis jarak
moran_knn <- moran.test(jabar_centroid$Y, lw_knn)
print(moran_knn)
## 
##  Moran I test under randomisation
## 
## data:  jabar_centroid$Y  
## weights: lw_knn    
## 
## Moran I statistic standard deviate = 3.3661, p-value = 0.0003812
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.35483779       -0.03846154        0.01365192

4.5.3. IDW

cat("\n")
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
cat("9. INTERPOLASI INVERSE DISTANCE WEIGHTING (IDW)\n")
## 9. INTERPOLASI INVERSE DISTANCE WEIGHTING (IDW)
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
# === TAMBAHKAN INI: BUAT GRID UNTUK PREDIKSI ===
cat("\nMembuat grid untuk interpolasi...\n")
## 
## Membuat grid untuk interpolasi...
# Dapatkan bounding box dari data
bbox_jabar <- st_bbox(jabar_merged)

# Buat grid reguler (sesuaikan cellsize untuk resolusi)
# cellsize dalam meter (karena CRS UTM)
# Gunakan 5000m (5km) untuk resolusi sedang
grid_sf <- st_make_grid(
  jabar_merged,
  cellsize = 5000,  # 5 km
  what = "centers"
)

# Convert ke sf object
grid_nn_sf <- st_sf(geometry = grid_sf)

# Crop grid hanya dalam batas Jawa Barat
grid_nn_sf <- st_intersection(grid_nn_sf, st_union(jabar_merged))

cat("Grid berhasil dibuat dengan", nrow(grid_nn_sf), "titik\n")
## Grid berhasil dibuat dengan 1483 titik
# Siapkan data untuk gstat
jabar_sp <- as(jabar_centroid, "Spatial")

# Buat grid untuk prediksi
# Extract koordinat dari sf object
grid_coords_df <- as.data.frame(st_coordinates(grid_nn_sf))
colnames(grid_coords_df) <- c("X", "Y")
grid_idw <- grid_coords_df
coordinates(grid_idw) <- ~X+Y
gridded(grid_idw) <- TRUE
proj4string(grid_idw) <- proj4string(jabar_sp)

# IDW dengan power = 2 (default)
cat("\nMelakukan IDW dengan power = 2...\n")
## 
## Melakukan IDW dengan power = 2...
idw_result <- gstat::idw(
  formula = Y ~ 1,
  locations = jabar_sp,
  newdata = grid_idw,
  idp = 2.0
)
## [inverse distance weighted interpolation]
# Convert hasil ke raster
idw_rast <- rast(idw_result["var1.pred"])
idw_rast <- mask(idw_rast, vect(jabar_merged))

# Statistik hasil
cat("\nStatistik Hasil IDW:\n")
## 
## Statistik Hasil IDW:
cat("Min  :", round(min(values(idw_rast), na.rm=TRUE), 2), "\n")
## Min  : 1.65
cat("Max  :", round(max(values(idw_rast), na.rm=TRUE), 2), "\n")
## Max  : 8.8
cat("Mean :", round(mean(values(idw_rast), na.rm=TRUE), 2), "\n")
## Mean : 6.56
cat("\n--- Peta Interpolasi IDW ---\n")
## 
## --- Peta Interpolasi IDW ---
peta_idw <- tm_shape(idw_rast) +
  tm_raster(
    col = "var1.pred",
    title = "TPT (%) Predicted",
    palette = "YlGn",
    style = "cont",
    alpha = 0.7
  ) +
  tm_shape(jabar_merged) +
  tm_borders(col = "black", lwd = 1) +
  tm_shape(jabar_centroid) +
  tm_dots(col = "red", size = 0.3, alpha = 0.8) +
  tm_layout(
    main.title = "Interpolasi IDW (power=2) - TPT Jawa Barat",
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
## `tm_scale_continuous()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values') to
##   'tm_scale_continuous(<HERE>)'
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_idw)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"

Peta hasil interpolasi IDW dengan parameter power = 2 menunjukkan pola sebaran Tingkat Pengangguran Terbuka (TPT) di Provinsi Jawa Barat berdasarkan kedekatan spasial antarwilayah. Warna yang lebih gelap merepresentasikan nilai TPT yang lebih tinggi, sedangkan warna yang lebih terang menunjukkan nilai TPT yang lebih rendah. Terlihat bahwa wilayah dengan TPT tinggi cenderung membentuk area tertentu dan dipengaruhi kuat oleh titik pengamatan terdekat. Metode IDW menghasilkan permukaan prediksi yang relatif mengikuti nilai titik observasi, sehingga perubahan nilai TPT tampak cukup tajam di sekitar lokasi data. Hal ini menunjukkan bahwa variasi TPT antarwilayah bersifat lokal dan sangat dipengaruhi oleh wilayah sekitarnya.

4.5.4. Kriging

4.5.4.1. Peta Prediksi

cat("\n")
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
cat("10. INTERPOLASI KRIGING\n")
## 10. INTERPOLASI KRIGING
cat(rep("=", 70), "\n", sep = "")
## ======================================================================
# --- a. Variogram Empiris ---
cat("\n--- Menghitung Variogram Empiris ---\n")
## 
## --- Menghitung Variogram Empiris ---
vario_emp <- variogram(Y ~ 1, jabar_sp)
print(vario_emp)
##    np      dist    gamma dir.hor dir.ver   id
## 1   2  2784.943 0.206650       0       0 var1
## 2   1 10902.642 1.232450       0       0 var1
## 3   1 14266.798 2.576450       0       0 var1
## 4   7 19390.844 1.495429       0       0 var1
## 5   6 26155.570 1.686617       0       0 var1
## 6   9 32088.944 2.394339       0       0 var1
## 7  13 38162.428 1.853742       0       0 var1
## 8   9 42765.590 2.818411       0       0 var1
## 9  14 49172.720 1.044939       0       0 var1
## 10 10 54769.031 0.547760       0       0 var1
## 11 12 60725.888 1.733458       0       0 var1
## 12 10 66371.230 0.758220       0       0 var1
## 13 22 71782.092 2.579275       0       0 var1
## 14 15 77636.070 1.943613       0       0 var1
## 15 12 83509.756 2.533617       0       0 var1
# Plot variogram
plot(vario_emp, main = "Variogram Empiris TPT")

# --- b. Fitting Variogram Model ---
cat("\n--- Fitting Model Variogram ---\n")
## 
## --- Fitting Model Variogram ---
# Coba beberapa model
vario_sph <- fit.variogram(vario_emp, vgm("Sph"))  # Spherical
## Warning in fit.variogram(vario_emp, vgm("Sph")): No convergence after 200
## iterations: try different initial values?
vario_exp <- fit.variogram(vario_emp, vgm("Exp"))  # Exponential
vario_gau <- fit.variogram(vario_emp, vgm("Gau"))  # Gaussian
## Warning in fit.variogram(vario_emp, vgm("Gau")): No convergence after 200
## iterations: try different initial values?
cat("\nModel Spherical:\n")
## 
## Model Spherical:
print(vario_sph)
##   model    psill    range
## 1   Nug 0.000000     0.00
## 2   Sph 1.872871 28076.29
cat("\nModel Exponential:\n")
## 
## Model Exponential:
print(vario_exp)
##   model    psill    range
## 1   Nug 0.000000     0.00
## 2   Exp 2.033949 15873.67
cat("\nModel Gaussian:\n")
## 
## Model Gaussian:
print(vario_gau)
##   model      psill    range
## 1   Nug 0.02619629    0.000
## 2   Gau 1.78631688 8648.648
# Pilih model terbaik (misal Spherical)
vario_model <- vario_sph
cat("\n✓ Menggunakan model: Spherical\n")
## 
## ✓ Menggunakan model: Spherical
# Plot fitted variogram
plot(vario_emp, vario_model, main = "Fitted Variogram Model (Spherical)")

# --- c. Ordinary Kriging ---
cat("\n--- Melakukan Ordinary Kriging ---\n")
## 
## --- Melakukan Ordinary Kriging ---
kriging_result <- krige(
  formula = Y ~ 1,
  locations = jabar_sp,
  newdata = grid_idw,
  model = vario_model
)
## [using ordinary kriging]
# Convert hasil ke raster
kriging_pred <- rast(kriging_result["var1.pred"])
kriging_var <- rast(kriging_result["var1.var"])

kriging_pred <- mask(kriging_pred, vect(jabar_merged))
kriging_var <- mask(kriging_var, vect(jabar_merged))

# Statistik hasil
cat("\nStatistik Hasil Kriging:\n")
## 
## Statistik Hasil Kriging:
cat("Prediksi:\n")
## Prediksi:
cat("  Min  :", round(min(values(kriging_pred), na.rm=TRUE), 2), "\n")
##   Min  : 2.12
cat("  Max  :", round(max(values(kriging_pred), na.rm=TRUE), 2), "\n")
##   Max  : 8.69
cat("  Mean :", round(mean(values(kriging_pred), na.rm=TRUE), 2), "\n")
##   Mean : 6.44
cat("\nVariance:\n")
## 
## Variance:
cat("  Min  :", round(min(values(kriging_var), na.rm=TRUE), 2), "\n")
##   Min  : 0.11
cat("  Max  :", round(max(values(kriging_var), na.rm=TRUE), 2), "\n")
##   Max  : 1.96
cat("  Mean :", round(mean(values(kriging_var), na.rm=TRUE), 2), "\n")
##   Mean : 1.66
# --- d. Peta Hasil Kriging ---
cat("\n--- Peta Interpolasi Kriging ---\n")
## 
## --- Peta Interpolasi Kriging ---
# Peta Prediksi
peta_kriging_pred <- tm_shape(kriging_pred) +
  tm_raster(
    col = "var1.pred",
    title = "TPT (%) Predicted",
    palette = "YlGn",
    style = "cont",
    alpha = 0.7
  ) +
  tm_shape(jabar_merged) +
  tm_borders(col = "black", lwd = 1) +
  tm_shape(jabar_centroid) +
  tm_dots(col = "red", size = 0.3, alpha = 0.8) +
  tm_layout(
    main.title = "Ordinary Kriging - Prediksi TPT Jawa Barat",
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
## `tm_scale_continuous()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values') to
##   'tm_scale_continuous(<HERE>)'
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_kriging_pred)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"

Peta hasil Ordinary Kriging menggambarkan prediksi TPT di Jawa Barat dengan mempertimbangkan struktur ketergantungan spasial antarwilayah melalui model variogram. Pola sebaran TPT pada peta ini terlihat lebih halus dibandingkan hasil interpolasi IDW, yang menunjukkan adanya proses smoothing akibat pemanfaatan informasi spasial secara global. Wilayah dengan nilai TPT tinggi dan rendah masih menunjukkan pola pengelompokan yang konsisten, sejalan dengan hasil uji autokorelasi spasial global sebelumnya. Hal ini mengindikasikan bahwa TPT di Jawa Barat tidak tersebar secara acak, melainkan memiliki pola spasial yang terstruktur.

4.5.4.2. Peta Varians

# Peta Variance (Uncertainty)
peta_kriging_var <- tm_shape(kriging_var) +
  tm_raster(
    col = "var1.var",
    title = "Kriging\nVariance",
    palette = "Reds",
    style = "cont",
    alpha = 0.7
  ) +
  tm_shape(jabar_merged) +
  tm_borders(col = "black", lwd = 1) +
  tm_shape(jabar_centroid) +
  tm_dots(col = "blue", size = 0.3, alpha = 0.8) +
  tm_layout(
    main.title = "Ordinary Kriging - Uncertainty (Variance)",
    main.title.size = 1,
    legend.outside = TRUE,
    legend.outside.position = "right",
    frame = TRUE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
## `tm_scale_continuous()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values') to
##   'tm_scale_continuous(<HERE>)'
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
print(peta_kriging_var)
## [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"

Peta varians hasil Ordinary Kriging menunjukkan tingkat ketidakpastian dari nilai prediksi TPT pada setiap lokasi. Warna yang lebih terang menunjukkan varians prediksi yang lebih tinggi, sedangkan warna yang lebih gelap menandakan ketidakpastian yang lebih rendah. Terlihat bahwa varians cenderung lebih kecil di sekitar lokasi titik pengamatan, karena informasi data tersedia secara langsung, dan semakin meningkat pada wilayah yang jauh dari titik observasi. Hal ini menunjukkan bahwa hasil prediksi pada daerah yang minim data memiliki tingkat ketidakpastian yang lebih besar. Dengan demikian, peta varians ini penting sebagai pelengkap peta prediksi, karena memberikan informasi mengenai tingkat keandalan hasil estimasi TPT secara spasial.

BAB 5. Kesimpulan

Berdasarkan hasil analisis spasial Tingkat Pengangguran Terbuka (TPT) di Provinsi Jawa Barat tahun 2024, dapat disimpulkan bahwa persebaran TPT menunjukkan pola spasial yang tidak acak dan cenderung mengelompok secara signifikan. Hal ini dibuktikan oleh nilai Moran’s I dan Geary’s C yang menunjukkan autokorelasi spasial positif, serta analisis LISA yang mengidentifikasi adanya klaster lokal pada beberapa wilayah. Pemodelan spasial menunjukkan bahwa model Spatial Autoregressive (SAR) merupakan model terbaik berdasarkan kriteria AIC, yang mengindikasikan bahwa nilai TPT suatu wilayah dipengaruhi oleh TPT wilayah sekitarnya. Variabel Indeks Pembangunan Manusia, jumlah penduduk miskin, dan PDRB memiliki pengaruh yang berbeda-beda terhadap TPT, baik secara langsung maupun melalui efek spasial. Hasil pemetaan prediksi dan residual menunjukkan bahwa wilayah perkotaan dan kawasan penyangga metropolitan memiliki TPT yang relatif lebih tinggi, sedangkan wilayah selatan dan timur Jawa Barat cenderung lebih rendah, sehingga menegaskan pentingnya pendekatan kebijakan berbasis wilayah dalam upaya penanggulangan pengangguran.

Daftar Pustaka

[1] Mahmud, A., & Pasaribu, E. (2021). Permodelan spasial pada analisis faktor yang mempengaruhi tingkat pengangguran terbuka Provinsi Bangka Belitung tahun 2018. Engineering, Mathematics and Computer Science Journal (EMACS). https://journal.binus.ac.id/index.php/EMACS/article/view/7034

[2] Badan Pusat Statistik Kota Bogor. (n.d.). Tingkat pengangguran terbuka menurut kabupaten/kota di Provinsi Jawa Barat (persen). Badan Pusat Statistik. Diakses 21 Desember 2025, dari https://bogorkota.bps.go.id/id/statistics-table/2/MTA1MSMy/tingkat-pengangguran-terbuka-menurut-kabupaten-kota-di-provinsi-jawa-barat–persen-.html

[3] Badan Pusat Statistik Kabupaten Majalengka. (n.d.). Indeks Pembangunan Manusia (IPM) Provinsi Jawa Barat menurut kabupaten/kota, 2023–2024. Badan Pusat Statistik. Diakses 21 Desember 2025, dari https://majalengkakab.bps.go.id/id/statistics-table/1/MzM5MiMx/indeks-pembangunan-manusia–ipm–provinsi-jawa-barat-menurut-kabupaten-kota–2023-2024.html

[4] Badan Pusat Statistik. (2025). Jumlah penduduk miskin berdasarkan kabupaten/kota di Jawa Barat [Dataset]. Open Data Provinsi Jawa Barat. Diakses 21 Desember 2025, dari https://opendata.jabarprov.go.id/id/dataset/jumlah-penduduk-miskin-berdasarkan-kabupatenkota-di-jawa-barat

[5] Badan Pusat Statistik. (2025). Produk Domestik Regional Bruto atas dasar harga konstan berdasarkan kabupaten/kota di Jawa Barat [Dataset]. Open Data Provinsi Jawa Barat. Diakses 21 Desember 2025, dari https://opendata.jabarprov.go.id/id/dataset/produk-domestik-regional-bruto-atas-dasar-harga-konstan-berdasarkan-kabupatenkota-di-jawa-barat

[6] E-Book Analisis Spasial I Gede Nyoman Mindra Jaya, M.Si., Ph.D. https://rpubs.com/mindra/1365004