Analisis Spasial Tingkat Pengangguran Terbuka (TPT) per Kabupaten/Kota di Provinsi Jawa Tengah Tahun 2024
Nama: Vindy Andrea Pratiwi
NIM: 24050123120014
Mata Kuliah: Statistika Spasial A
Program Studi: Statistika
Universitas Diponegoro
15-05-2026
- BAB 1. Pendahuluan
- BAB 2. Tinjauan Pustaka
- BAB 3. Metode Penelitian
- BAB 4. Hasil dan Pembahasan
- 4.1. Load Library
- 4.2. Persiapan Data
- 4.3. Statistika Deskriptif
- 4.4. Peta Deskriptif TPT
- 4.5. Identifikasi Pola Spasial: Acak, Cluster, atau Seragam
- 4.6. Matriks Pembobot Spasial
- 4.7. Autokorelasi Spasial Global
- 4.8. Autokorelasi Spasial Lokal
- 4.9. Kriging Global dan Kriging Lokal
- 4.10. Pemodelan Spasial Menggunakan SAR
- BAB 5. Kesimpulan
- Daftar Pustaka
- Lampiran: Session Info
BAB 1. Pendahuluan
1.1. Latar Belakang
Tingkat Pengangguran Terbuka (TPT) merupakan salah satu indikator utama dalam menggambarkan kondisi ketenagakerjaan suatu wilayah. TPT menunjukkan persentase angkatan kerja yang belum memperoleh pekerjaan, sedang mencari pekerjaan, atau sedang mempersiapkan usaha. Nilai TPT yang tinggi dapat mengindikasikan adanya ketidakseimbangan antara penawaran dan permintaan tenaga kerja, keterbatasan kesempatan kerja, serta adanya perbedaan karakteristik sosial-ekonomi antarwilayah.
Provinsi Jawa Tengah terdiri atas kabupaten/kota dengan karakteristik ekonomi, sosial, pendidikan, dan kependudukan yang beragam. Perbedaan karakteristik tersebut menyebabkan TPT tidak selalu tersebar secara merata. Wilayah yang berdekatan secara geografis dapat memiliki kondisi ketenagakerjaan yang mirip karena adanya keterkaitan aktivitas ekonomi, mobilitas penduduk, konektivitas transportasi, dan pengaruh pusat-pusat pertumbuhan. Oleh sebab itu, analisis yang hanya menggunakan pendekatan statistik biasa belum cukup untuk menjelaskan kemungkinan adanya ketergantungan antarwilayah.
Analisis spasial digunakan untuk mengidentifikasi apakah persebaran TPT membentuk pola acak, mengelompok (cluster), atau seragam. Selain itu, analisis autokorelasi spasial seperti Moran’s I, Geary’s C, dan Getis-Ord dapat digunakan untuk mengukur hubungan spasial secara global maupun lokal. Melalui pendekatan ini, wilayah dengan karakteristik TPT tinggi atau rendah yang berdekatan dapat diidentifikasi secara lebih jelas.
Selain pengujian autokorelasi, interpolasi spasial seperti Kriging digunakan untuk membentuk permukaan prediksi TPT berdasarkan lokasi kabupaten/kota. Kriging mempertimbangkan struktur ketergantungan spasial melalui variogram sehingga dapat memberikan prediksi dan ukuran ketidakpastian. Selanjutnya, pemodelan spasial seperti Spatial Autoregressive Model (SAR) digunakan untuk menjelaskan hubungan TPT dengan variabel penjelas seperti Rata-rata Lama Sekolah (RLS), Indeks Pembangunan Manusia (IPM), jumlah penduduk miskin, dan Produk Domestik Regional Bruto (PDRB), dengan tetap mempertimbangkan pengaruh wilayah tetangga.
1.2. Rumusan Masalah
Berdasarkan latar belakang tersebut, rumusan masalah dalam penelitian ini adalah sebagai berikut.
- Bagaimana pola spasial TPT kabupaten/kota di Provinsi Jawa Tengah, apakah acak, mengelompok, atau seragam?
- Matriks pembobot spasial apa yang sesuai untuk data kabupaten/kota di Provinsi Jawa Tengah?
- Apakah terdapat autokorelasi spasial global dan lokal berdasarkan Moran’s I, Geary’s C, dan Getis-Ord?
- Bagaimana hasil interpolasi TPT menggunakan Kriging global dan Kriging lokal?
- Bagaimana penerapan model SAR dalam menjelaskan TPT berdasarkan variabel sosial-ekonomi dan efek spasial?
1.3. Tujuan Penelitian
Tujuan penelitian ini adalah sebagai berikut.
- Mengidentifikasi pola spasial TPT kabupaten/kota di Provinsi Jawa Tengah.
- Membentuk matriks pembobot spasial berbasis kontiguitas dan jarak.
- Menguji autokorelasi spasial global dan lokal menggunakan Moran’s I, Geary’s C, dan Getis-Ord.
- Melakukan interpolasi spasial menggunakan Kriging global dan Kriging lokal.
- Menerapkan model SAR untuk menganalisis TPT dengan variabel RLS, IPM, Jumlah Penduduk Miskin, dan PDRB.
BAB 2. Tinjauan Pustaka
2.1. Analisis Pola Spasial
Analisis pola spasial digunakan untuk mengetahui apakah susunan
lokasi pengamatan membentuk pola acak (random), mengelompok
(clustered), atau seragam/teratur (dispersed). Pada
data wilayah kabupaten/kota, pola spasial dapat dilihat melalui titik
representatif dari setiap wilayah, misalnya titik
st_point_on_surface, kemudian dihitung jarak tetangga
terdekatnya. Salah satu ukuran yang digunakan adalah Average Nearest
Neighbor Index (ANNI), yaitu perbandingan rata-rata jarak tetangga
terdekat yang diamati dengan rata-rata jarak harapan pada pola acak.
Rumus Average Nearest Neighbor Index adalah:
\[ ANNI = \frac{\bar{d}_{obs}}{\bar{d}_{exp}} \]
dengan:
\[ \bar{d}_{obs}=\frac{1}{n}\sum_{i=1}^{n}d_i \]
\[ \bar{d}_{exp}=\frac{0.5}{\sqrt{n/A}} \]
\[ SE=\frac{0.26136}{\sqrt{n^2/A}} \]
\[ Z=\frac{\bar{d}_{obs}-\bar{d}_{exp}}{SE} \]
Keterangan:
- \(ANNI\) adalah indeks tetangga terdekat rata-rata.
- \(\bar{d}_{obs}\) adalah rata-rata jarak tetangga terdekat hasil pengamatan.
- \(\bar{d}_{exp}\) adalah rata-rata jarak harapan jika titik tersebar acak.
- \(d_i\) adalah jarak dari lokasi ke-\(i\) menuju tetangga terdekatnya.
- \(n\) adalah jumlah wilayah pengamatan.
- \(A\) adalah luas wilayah studi.
- \(Z\) digunakan untuk menguji apakah pola berbeda signifikan dari pola acak.
Interpretasi ANNI yaitu \(ANNI < 1\) menunjukkan pola mengelompok, \(ANNI \approx 1\) menunjukkan pola acak, dan \(ANNI > 1\) menunjukkan pola seragam/teratur. Jika p-value < 0,05, maka pola tersebut dianggap signifikan secara statistik.
2.2. Matriks Pembobot Spasial
Matriks pembobot spasial (\(W\)) merupakan matriks yang menggambarkan hubungan kedekatan atau ketetanggaan antarunit wilayah. Elemen \(w_{ij}\) menunjukkan kekuatan hubungan spasial antara wilayah \(i\) dan wilayah \(j\). Pada data kabupaten/kota, pembobot kontiguitas sering digunakan karena wilayah administratif dapat dikatakan bertetangga jika memiliki batas wilayah yang bersinggungan.
Secara umum, matriks pembobot spasial ditulis sebagai:
\[ W = [w_{ij}], \quad i,j=1,2,\ldots,n \]
Untuk pembobot kontiguitas biner:
\[ w_{ij}=\begin{cases} 1, & \text{jika wilayah } i \text{ bertetangga dengan wilayah } j \\ 0, & \text{jika wilayah } i \text{ tidak bertetangga dengan wilayah } j \end{cases} \]
Pada Queen Contiguity, dua wilayah dianggap bertetangga jika berbagi sisi atau titik sudut. Pada Rook Contiguity, dua wilayah dianggap bertetangga hanya jika berbagi sisi. Dalam penelitian ini, pembobot utama yang digunakan adalah Queen Contiguity karena lebih inklusif untuk menggambarkan hubungan antarwilayah kabupaten/kota.
Agar setiap baris memiliki total bobot yang sama, dilakukan standardisasi baris:
\[ w_{ij}^{*}=\frac{w_{ij}}{\sum_{j=1}^{n}w_{ij}} \]
Keterangan:
- \(w_{ij}\) adalah bobot awal antara wilayah \(i\) dan \(j\).
- \(w_{ij}^{*}\) adalah bobot setelah standardisasi baris.
- \(\sum_{j=1}^{n}w_{ij}\) adalah jumlah tetangga atau total bobot pada wilayah \(i\).
Standardisasi baris membuat nilai spasial lag, misalnya \(Wy\), dapat diinterpretasikan sebagai rata-rata nilai variabel pada wilayah tetangga.
2.3. Autokorelasi Spasial
Autokorelasi spasial menunjukkan hubungan antara nilai suatu variabel pada satu wilayah dengan nilai variabel pada wilayah lain yang berdekatan. Autokorelasi positif menunjukkan bahwa wilayah yang berdekatan memiliki nilai yang mirip, sedangkan autokorelasi negatif menunjukkan bahwa wilayah berdekatan memiliki nilai yang berbeda. Pada penelitian ini, autokorelasi spasial dianalisis secara global dan lokal menggunakan Moran’s I, Geary’s C, dan Getis-Ord.
2.3.1. Moran’s I Global
Moran’s I merupakan ukuran autokorelasi spasial global yang digunakan untuk mengetahui apakah suatu variabel memiliki pola spasial secara keseluruhan. Rumus Moran’s I adalah:
\[ I = \frac{n}{S_0}\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2} \]
dengan:
\[ S_0 = \sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij} \]
Keterangan:
- \(I\) adalah indeks Moran.
- \(n\) adalah jumlah wilayah.
- \(x_i\) dan \(x_j\) adalah nilai variabel pada wilayah \(i\) dan \(j\).
- \(\bar{x}\) adalah rata-rata variabel.
- \(w_{ij}\) adalah bobot spasial antara wilayah \(i\) dan \(j\).
- \(S_0\) adalah jumlah seluruh bobot spasial.
Nilai \(I>0\) menunjukkan autokorelasi spasial positif, \(I<0\) menunjukkan autokorelasi spasial negatif, dan \(I\approx 0\) menunjukkan tidak terdapat pola spasial yang jelas.
2.3.2. Moran’s I Lokal / LISA
Moran’s I lokal atau Local Indicator of Spatial Association (LISA) digunakan untuk mengidentifikasi wilayah mana saja yang membentuk klaster lokal. Rumus Moran lokal untuk wilayah ke-\(i\) adalah:
\[ I_i = z_i \sum_{j=1}^{n}w_{ij}z_j \]
atau dalam bentuk standar:
\[ z_i = \frac{x_i-\bar{x}}{s} \]
Keterangan:
- \(I_i\) adalah Moran lokal wilayah ke-\(i\).
- \(z_i\) adalah nilai baku wilayah ke-\(i\).
- \(z_j\) adalah nilai baku wilayah tetangga ke-\(j\).
- \(w_{ij}\) adalah bobot spasial antara wilayah \(i\) dan \(j\).
- \(s\) adalah simpangan baku variabel.
LISA menghasilkan empat kategori utama, yaitu High-High (wilayah bernilai tinggi dikelilingi wilayah tinggi), Low-Low (wilayah rendah dikelilingi wilayah rendah), High-Low (wilayah tinggi dikelilingi wilayah rendah), dan Low-High (wilayah rendah dikelilingi wilayah tinggi).
2.3.3. Geary’s C Global
Geary’s C digunakan untuk mengukur autokorelasi spasial global, tetapi lebih sensitif terhadap perbedaan lokal antarwilayah bertetangga. Rumus Geary’s C adalah:
\[ C = \frac{(n-1)}{2S_0}\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-x_j)^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2} \]
Keterangan:
- \(C\) adalah indeks Geary.
- \(n\) adalah jumlah wilayah.
- \(x_i\) dan \(x_j\) adalah nilai variabel pada wilayah \(i\) dan \(j\).
- \(\bar{x}\) adalah rata-rata variabel.
- \(w_{ij}\) adalah bobot spasial.
- \(S_0\) adalah jumlah seluruh bobot spasial.
Interpretasi Geary’s C yaitu \(C<1\) menunjukkan autokorelasi spasial positif, \(C>1\) menunjukkan autokorelasi spasial negatif, dan \(C\approx 1\) menunjukkan tidak terdapat autokorelasi spasial.
2.3.4. Geary’s C Lokal
Geary’s C lokal digunakan untuk melihat wilayah yang memiliki perbedaan lokal kuat dibandingkan wilayah tetangganya. Salah satu bentuk rumus local Geary untuk wilayah ke-\(i\) adalah:
\[ C_i = \sum_{j=1}^{n}w_{ij}(x_i-x_j)^2 \]
Keterangan:
- \(C_i\) adalah nilai Geary lokal wilayah ke-\(i\).
- \(x_i\) adalah nilai variabel pada wilayah \(i\).
- \(x_j\) adalah nilai variabel pada wilayah tetangga \(j\).
- \(w_{ij}\) adalah bobot spasial antara wilayah \(i\) dan \(j\).
Nilai \(C_i\) yang besar menunjukkan perbedaan nilai antara suatu wilayah dengan tetangganya relatif kuat. Jika signifikan, wilayah tersebut dapat diinterpretasikan sebagai lokasi yang memiliki ketidakmiripan spasial lokal.
2.3.5. Getis-Ord G Global
Getis-Ord G global digunakan untuk menguji apakah nilai tinggi atau nilai rendah suatu variabel terkonsentrasi secara spasial. Rumus Getis-Ord G global adalah:
\[ G = \frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}x_ix_j}{\sum_{i=1}^{n}\sum_{j=1, j\neq i}^{n}x_ix_j} \]
Keterangan:
- \(G\) adalah statistik Getis-Ord global.
- \(x_i\) dan \(x_j\) adalah nilai variabel pada wilayah \(i\) dan \(j\).
- \(w_{ij}\) adalah bobot spasial antara wilayah \(i\) dan \(j\).
Nilai Getis-Ord G yang signifikan menunjukkan adanya konsentrasi spasial nilai tinggi atau rendah. Dalam penelitian TPT, statistik ini membantu melihat apakah wilayah dengan TPT tinggi cenderung berkumpul di area tertentu.
2.3.6. Getis-Ord Lokal
Getis-Ord lokal digunakan untuk mengidentifikasi hot spot dan cold spot. Rumus statistik Getis-Ord lokal adalah:
\[ G_i = \frac{\sum_{j=1}^{n}w_{ij}x_j}{\sum_{j=1}^{n}x_j} \]
Bentuk baku yang sering digunakan adalah \(Z(G_i)\), yaitu skor-Z dari statistik lokal. Interpretasinya adalah:
- \(Z(G_i)>0\) dan signifikan menunjukkan hot spot, yaitu wilayah yang dikelilingi nilai tinggi.
- \(Z(G_i)<0\) dan signifikan menunjukkan cold spot, yaitu wilayah yang dikelilingi nilai rendah.
- Tidak signifikan menunjukkan tidak ada konsentrasi lokal yang cukup kuat secara statistik.
2.4. Kriging
Kriging merupakan metode interpolasi geostatistik yang memperkirakan nilai pada lokasi tidak teramati berdasarkan nilai pada lokasi pengamatan dan struktur ketergantungan spasial melalui variogram. Dalam penelitian ini digunakan Ordinary Kriging, baik global maupun lokal. Ordinary Kriging mengasumsikan rata-rata lokal tidak diketahui tetapi konstan pada lingkungan prediksi.
Prediksi Ordinary Kriging ditulis sebagai:
\[ \hat{Z}(s_0)=\sum_{i=1}^{n}\lambda_i Z(s_i) \]
dengan syarat:
\[ \sum_{i=1}^{n}\lambda_i=1 \]
Keterangan:
- \(\hat{Z}(s_0)\) adalah nilai prediksi pada lokasi \(s_0\).
- \(Z(s_i)\) adalah nilai pengamatan pada lokasi \(s_i\).
- \(\lambda_i\) adalah bobot kriging untuk lokasi pengamatan ke-\(i\).
- \(n\) adalah jumlah titik pengamatan yang digunakan.
Kriging menggunakan variogram untuk menggambarkan hubungan jarak dengan variasi nilai. Variogram empiris dirumuskan sebagai:
\[ \hat{\gamma}(h)=\frac{1}{2N(h)}\sum_{N(h)}[Z(s_i)-Z(s_i+h)]^2 \]
Keterangan:
- \(\hat{\gamma}(h)\) adalah semivariogram pada jarak \(h\).
- \(N(h)\) adalah jumlah pasangan lokasi yang memiliki jarak sekitar \(h\).
- \(Z(s_i)\) dan \(Z(s_i+h)\) adalah nilai pada dua lokasi yang dipisahkan jarak \(h\).
Ordinary Kriging global menggunakan seluruh titik pengamatan dalam proses prediksi. Ordinary Kriging lokal membatasi titik yang digunakan, misalnya hanya maksimum 8 tetangga terdekat, sehingga hasilnya lebih sensitif terhadap variasi setempat.
2.5. Spatial Autoregressive Model (SAR)
Spatial Autoregressive Model (SAR) digunakan ketika variabel dependen pada suatu wilayah dipengaruhi oleh nilai variabel dependen pada wilayah tetangga. Model ini memasukkan komponen lag spasial \(Wy\) sehingga mampu menangkap efek ketergantungan antarwilayah. Persamaan umum model SAR adalah:
\[ y = \rho Wy + X\beta + \varepsilon \]
atau dapat ditulis sebagai:
\[ (I-\rho W)y = X\beta + \varepsilon \]
Keterangan:
- \(y\) adalah vektor variabel dependen, yaitu TPT.
- \(W\) adalah matriks pembobot spasial.
- \(Wy\) adalah lag spasial variabel dependen atau rata-rata TPT wilayah tetangga.
- \(\rho\) adalah koefisien autoregresif spasial.
- \(X\) adalah matriks variabel independen, yaitu RLS, IPM, jumlah penduduk miskin, dan PDRB.
- \(\beta\) adalah vektor parameter regresi.
- \(\varepsilon\) adalah galat acak.
Jika \(\rho\) signifikan dan bernilai positif, maka TPT suatu wilayah cenderung meningkat ketika TPT wilayah tetangganya juga tinggi. Model SAR dipilih dalam penelitian ini karena tujuan analisis adalah melihat pengaruh karakteristik sosial-ekonomi dan pengaruh TPT wilayah sekitar terhadap TPT kabupaten/kota di Jawa Tengah.
BAB 3. Metode Penelitian
3.1. Data dan Sumber Data
Data yang digunakan dalam penelitian ini adalah data kabupaten/kota
di Provinsi Jawa Tengah yang diperoleh dari file Excel yang dilampirkan,
yaitu Data Spasial Vindy.xlsx. Data batas wilayah
kabupaten/kota diperoleh dari shapefile
gadm41_IDN_2.shp.
Variabel yang digunakan adalah sebagai berikut.
| Variabel | Keterangan |
|---|---|
| Y | Tingkat Pengangguran Terbuka (TPT) |
| X1 | Rata-rata Lama Sekolah (RLS) |
| X2 | Indeks Pembangunan Manusia (IPM) |
| X3 | Jumlah Penduduk Miskin |
| X4 | Produk Domestik Regional Bruto (PDRB) |
3.2. Tahapan Analisis
Tahapan analisis dalam penelitian ini adalah sebagai berikut.
- Membaca dan membersihkan data Excel serta shapefile.
- Menggabungkan data atribut dengan data spasial kabupaten/kota Jawa Tengah.
- Melakukan statistika deskriptif dan pemetaan sebaran TPT.
- Mengidentifikasi pola spasial menggunakan Average Nearest Neighbor Index.
- Membentuk matriks pembobot spasial berbasis Queen Contiguity dan K-nearest neighbors.
- Menguji autokorelasi spasial global dan lokal menggunakan Moran’s I, Geary’s C, dan Getis-Ord.
- Melakukan interpolasi Kriging global dan Kriging lokal.
- Menerapkan model SAR sebagai model regresi spasial.
- Menyusun interpretasi dan kesimpulan.
BAB 4. Hasil dan Pembahasan
4.1. Load Library
required_packages <- c(
"sf", "readxl", "dplyr", "tidyr", "stringr", "spdep", "spatialreg",
"ggplot2", "gstat", "sp", "knitr", "car", "broom"
)
new_packages <- required_packages[!(required_packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) {
install.packages(new_packages, dependencies = TRUE)
}
invisible(lapply(required_packages, library, character.only = TRUE))
sf::sf_use_s2(FALSE)
# Fungsi pembulatan aman untuk tabel: hanya kolom numerik yang dibulatkan.
round_df <- function(x, digits = 3) {
if (is.data.frame(x)) {
return(dplyr::mutate(x, dplyr::across(where(is.numeric), ~ round(.x, digits))))
}
if (is.matrix(x)) return(round(x, digits))
x
}4.2. Persiapan Data
# 1. Mencari path file secara otomatis
# Letakkan file Rmd, Excel, dan komponen shapefile dalam satu folder agar bagian ini langsung berjalan.
# Jika file berada di folder lain, cukup ubah base_dir secara manual.
input_dir <- tryCatch(dirname(knitr::current_input()), error = function(e) NA_character_)
base_dir_candidates <- unique(c(
getwd(),
input_dir,
"C:/Users/ACER SWIFT/Documents/Vindy Andrea/Kuliah/Semester 6/STATISTIKA SPASIAL/Tugas Poyek",
"/mnt/data"
))
base_dir_candidates <- base_dir_candidates[!is.na(base_dir_candidates) & nzchar(base_dir_candidates)]
excel_candidates <- c("Data Spasial Vindy.xlsx", "Data Spasial Vindy(4).xlsx")
shp_candidates <- c("gadm41_IDN_2.shp", "gadm41_IDN_2(3).shp")
find_file <- function(dirs, files) {
all_paths <- as.vector(outer(dirs, files, file.path))
found <- all_paths[file.exists(all_paths)]
if (length(found) == 0) return(NA_character_)
found[1]
}
data_path <- find_file(base_dir_candidates, excel_candidates)
shp_path <- find_file(base_dir_candidates, shp_candidates)
base_dir <- ifelse(!is.na(data_path), dirname(data_path), getwd())
cat("Working directory saat Knit:\n", getwd(), "\n\n")## Working directory saat Knit:
## C:/Users/ACER SWIFT/Downloadscat("Folder data yang dipakai:\n", base_dir, "\n\n")## Folder data yang dipakai:
## C:/Users/ACER SWIFT/Downloadscat("File Excel dipakai:\n", data_path, "\n")## File Excel dipakai:
## C:/Users/ACER SWIFT/Downloads/Data Spasial Vindy.xlsxcat("File shapefile dipakai:\n", shp_path, "\n\n")## File shapefile dipakai:
## C:/Users/ACER SWIFT/Downloads/gadm41_IDN_2.shpif (is.na(data_path) || !file.exists(data_path)) {
stop(paste0(
"File Excel tidak ditemukan. Pastikan salah satu nama file berikut ada di folder Rmd: ",
paste(excel_candidates, collapse = ", ")
))
}
if (is.na(shp_path) || !file.exists(shp_path)) {
stop(paste0(
"File .shp tidak ditemukan. Pastikan salah satu nama file berikut ada di folder Rmd: ",
paste(shp_candidates, collapse = ", "),
" dan komponen .dbf, .shx, .prj tersedia."
))
}
# Cek komponen shapefile wajib agar st_read tidak gagal.
shp_base <- tools::file_path_sans_ext(shp_path)
required_shp_parts <- paste0(shp_base, c(".shp", ".shx", ".dbf"))
missing_shp_parts <- required_shp_parts[!file.exists(required_shp_parts)]
if (length(missing_shp_parts) > 0) {
stop(paste("Komponen shapefile wajib belum lengkap:", paste(basename(missing_shp_parts), collapse = ", ")))
}
# File .prj dipakai untuk CRS. Jika tidak terbaca, nanti CRS tetap diset manual ke EPSG:4326.
optional_prj <- paste0(shp_base, ".prj")
if (!file.exists(optional_prj)) {
cat("Catatan: file .prj tidak ditemukan. CRS akan diset manual ke EPSG:4326 setelah st_read().\n")
}## Catatan: file .prj tidak ditemukan. CRS akan diset manual ke EPSG:4326 setelah st_read().# 2. Fungsi membersihkan nama wilayah
clean_name <- function(x) {
x %>%
as.character() %>%
stringr::str_to_upper() %>%
stringr::str_replace_all("^\\s*[0-9]{4}\\s*", "") %>%
stringr::str_replace_all("\\s+", " ") %>%
stringr::str_trim()
}
make_admin_name <- function(type, name) {
type2 <- stringr::str_to_upper(type)
name2 <- clean_name(name)
dplyr::case_when(
stringr::str_detect(name2, "^KOTA\\s+") ~ name2,
stringr::str_detect(name2, "^KABUPATEN\\s+") ~ name2,
stringr::str_detect(type2, "KOTA|MUNICIPALITY") ~ clean_name(paste("Kota", name2)),
stringr::str_detect(type2, "KABUPATEN|REGENCY") ~ clean_name(paste("Kabupaten", name2)),
TRUE ~ name2
)
}
# 3. Membaca data Excel
data_raw <- readxl::read_excel(data_path)
# Samakan nama kolom agar aman jika ada spasi di header Excel.
names(data_raw) <- names(data_raw) %>%
stringr::str_replace_all("\\s+", "") %>%
stringr::str_trim()
needed_cols <- c("KabKota", "TPT", "RLS", "IPM", "Miskin", "PDRB")
missing_cols <- setdiff(needed_cols, names(data_raw))
if (length(missing_cols) > 0) {
stop(paste("Kolom berikut tidak ditemukan pada Excel:", paste(missing_cols, collapse = ", ")))
}
# Format data sesuai variabel pada file Vindy.
data_tbl <- data_raw %>%
dplyr::mutate(
NAMA = clean_name(KabKota),
TPT = as.numeric(TPT),
RLS = as.numeric(RLS),
IPM = as.numeric(IPM),
Miskin = as.numeric(Miskin),
PDRB = as.numeric(PDRB)
) %>%
dplyr::select(NAMA, KabKota, TPT, RLS, IPM, Miskin, PDRB)
# Cek NA sebelum imputasi.
na_awal <- sapply(data_tbl[, c("TPT", "RLS", "IPM", "Miskin", "PDRB")], function(x) sum(is.na(x)))
cat("Jumlah NA sebelum imputasi:
")## Jumlah NA sebelum imputasi:print(na_awal)## TPT RLS IPM Miskin PDRB
## 0 1 0 0 0# Agar analisis regresi dan spasial tidak error, nilai kosong pada variabel numerik
# diisi dengan rata-rata variabel tersebut.
for (v in c("TPT", "RLS", "IPM", "Miskin", "PDRB")) {
data_tbl[[v]][is.na(data_tbl[[v]])] <- mean(data_tbl[[v]], na.rm = TRUE)
}
# 4. Membaca shapefile GADM Level 2 dan memilih Jawa Tengah
shp_raw <- sf::st_read(shp_path, quiet = TRUE)
# GADM Indonesia level 2 menggunakan sistem koordinat geografis WGS84 (EPSG:4326).
# Pada beberapa komputer, CRS shapefile bisa tidak terbaca walaupun file .prj ada.
# Karena itu, jika CRS kosong, CRS diset manual terlebih dahulu agar st_transform() tidak error.
if (is.na(sf::st_crs(shp_raw))) {
cat("CRS shapefile tidak terbaca. CRS diset manual ke EPSG:4326 (WGS84) sebelum transformasi.\n")
shp_raw <- sf::st_set_crs(shp_raw, 4326)
}## CRS shapefile tidak terbaca. CRS diset manual ke EPSG:4326 (WGS84) sebelum transformasi.if (!all(c("NAME_1", "NAME_2", "TYPE_2") %in% names(shp_raw))) {
stop("Kolom NAME_1, NAME_2, atau TYPE_2 tidak ditemukan pada shapefile GADM level 2.")
}
jateng_shp <- shp_raw %>%
dplyr::filter(NAME_1 == "Jawa Tengah", TYPE_2 != "Water Body") %>%
sf::st_make_valid() %>%
dplyr::mutate(NAMA = make_admin_name(TYPE_2, NAME_2))
if (nrow(jateng_shp) == 0) {
stop("Tidak ada wilayah dengan NAME_1 == 'Jawa Tengah' pada shapefile. Pastikan shapefile yang digunakan adalah GADM Indonesia level 2.")
}
# 5. Cek kecocokan nama wilayah
cek_data_tidak_ada_peta <- dplyr::anti_join(data_tbl, sf::st_drop_geometry(jateng_shp), by = "NAMA")
cek_peta_tidak_ada_data <- dplyr::anti_join(sf::st_drop_geometry(jateng_shp), data_tbl, by = "NAMA")
cat("Wilayah pada data tetapi tidak ada di peta:
")## Wilayah pada data tetapi tidak ada di peta:print(cek_data_tidak_ada_peta[, c("NAMA", "KabKota")])## # A tibble: 0 × 2
## # ℹ 2 variables: NAMA <chr>, KabKota <chr>cat("Wilayah pada peta tetapi tidak ada di data:
")## Wilayah pada peta tetapi tidak ada di data:print(cek_peta_tidak_ada_data[, c("NAMA", "NAME_2", "TYPE_2")])## [1] NAMA NAME_2 TYPE_2
## <0 rows> (or 0-length row.names)if (nrow(cek_data_tidak_ada_peta) > 0 | nrow(cek_peta_tidak_ada_data) > 0) {
stop("Masih ada nama wilayah yang belum cocok. Periksa tabel hasil anti_join di atas.")
}
# 6. Merge data atribut (kalau yang saya pake merupakan data yang terdapat di Excel) dan spasial
jateng <- jateng_shp %>%
dplyr::left_join(data_tbl, by = "NAMA") %>%
dplyr::arrange(NAMA) %>%
dplyr::mutate(
LABEL = stringr::str_replace_all(NAMA, "KABUPATEN\\s+|KOTA\\s+", ""),
LABEL = stringr::str_to_title(LABEL)
)
# Gunakan proyeksi UTM Zona 49S agar analisis jarak dan kriging menggunakan satuan meter.
# Jika CRS masih kosong setelah merge, set kembali ke EPSG:4326 sebelum transformasi.
if (is.na(sf::st_crs(jateng))) {
sf::st_crs(jateng) <- 4326
}
jateng <- sf::st_transform(jateng, 32749)
# Titik label dipakai pada seluruh peta agar setiap kabupaten/kota memiliki tulisan nama.
label_points <- sf::st_point_on_surface(jateng)
cat("Jumlah kabupaten/kota yang berhasil digabung:", nrow(jateng), "
")## Jumlah kabupaten/kota yang berhasil digabung: 35print(jateng[, c("NAMA", "TPT", "RLS", "IPM", "Miskin", "PDRB")])## Simple feature collection with 35 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 230132 ymin: 9092271 xmax: 576407.7 ymax: 9367063
## Projected CRS: WGS 84 / UTM zone 49S
## First 10 features:
## NAMA TPT RLS IPM Miskin PDRB
## 1 KABUPATEN BANJARNEGARA 5.57 6.870000 69.60 137.68 27685
## 2 KABUPATEN BANYUMAS 6.18 8.332941 74.52 207.78 40236
## 3 KABUPATEN BATANG 5.67 7.080000 70.74 68.85 36295
## 4 KABUPATEN BLORA 3.67 7.260000 71.39 99.14 36598
## 5 KABUPATEN BOYOLALI 3.16 8.170000 75.96 95.96 41699
## 6 KABUPATEN BREBES 8.35 6.410000 68.46 283.28 30040
## 7 KABUPATEN CILACAP 7.83 7.400000 72.38 186.08 67203
## 8 KABUPATEN DEMAK 4.75 8.280000 74.57 142.91 28313
## 9 KABUPATEN GROBOGAN 3.23 7.290000 72.02 159.00 24908
## 10 KABUPATEN JEPARA 3.34 8.270000 74.32 80.84 32894
## geometry
## 1 POLYGON ((325676.6 9166909,...
## 2 MULTIPOLYGON (((328016 9158...
## 3 POLYGON ((380429.7 9204672,...
## 4 POLYGON ((548863.8 9186970,...
## 5 POLYGON ((448766 9156579, 4...
## 6 POLYGON ((282776.3 9187000,...
## 7 MULTIPOLYGON (((266525.8 91...
## 8 POLYGON ((451609.5 9215806,...
## 9 MULTIPOLYGON (((455031.8 92...
## 10 MULTIPOLYGON (((427149.5 93...Berdasarkan hasil persiapan data, seluruh kabupaten/kota pada data atribut berhasil dicocokkan dengan shapefile Provinsi Jawa Tengah. Nilai kosong pada variabel numerik ditangani dengan imputasi rata-rata agar seluruh prosedur analisis dapat dijalankan tanpa menghapus wilayah pengamatan.
4.3. Statistika Deskriptif
vars <- c("TPT", "RLS", "IPM", "Miskin", "PDRB")
stat_deskriptif <- data.frame(
Variabel = vars,
Minimum = sapply(vars, function(v) min(jateng[[v]], na.rm = TRUE)),
Q1 = sapply(vars, function(v) quantile(jateng[[v]], 0.25, na.rm = TRUE)),
Median = sapply(vars, function(v) median(jateng[[v]], na.rm = TRUE)),
Mean = sapply(vars, function(v) mean(jateng[[v]], na.rm = TRUE)),
Q3 = sapply(vars, function(v) quantile(jateng[[v]], 0.75, na.rm = TRUE)),
Maximum = sapply(vars, function(v) max(jateng[[v]], na.rm = TRUE)),
SD = sapply(vars, function(v) sd(jateng[[v]], na.rm = TRUE)),
row.names = NULL
)
knitr::kable(round_df(stat_deskriptif, 3), caption = "Statistika Deskriptif Variabel Penelitian")| Variabel | Minimum | Q1 | Median | Mean | Q3 | Maximum | SD |
|---|---|---|---|---|---|---|---|
| TPT | 2.35 | 3.500 | 3.97 | 4.520 | 5.32 | 8.35 | 1.495 |
| RLS | 6.41 | 7.380 | 7.87 | 8.333 | 9.27 | 11.48 | 1.373 |
| IPM | 68.46 | 71.855 | 74.32 | 74.749 | 77.22 | 85.72 | 4.502 |
| Miskin | 7.25 | 72.860 | 95.96 | 105.838 | 139.76 | 283.28 | 59.435 |
| PDRB | 22652.00 | 30632.000 | 38792.00 | 51391.371 | 54210.50 | 156607.00 | 33188.236 |
# Korelasi antarvariabel numerik
cor_matrix <- cor(sf::st_drop_geometry(jateng[, vars]), use = "complete.obs")
knitr::kable(round_df(cor_matrix, 3), caption = "Matriks Korelasi Variabel Penelitian")| TPT | RLS | IPM | Miskin | PDRB | |
|---|---|---|---|---|---|
| TPT | 1.000 | -0.182 | -0.198 | 0.464 | 0.002 |
| RLS | -0.182 | 1.000 | 0.976 | -0.644 | 0.759 |
| IPM | -0.198 | 0.976 | 1.000 | -0.625 | 0.781 |
| Miskin | 0.464 | -0.644 | -0.625 | 1.000 | -0.464 |
| PDRB | 0.002 | 0.759 | 0.781 | -0.464 | 1.000 |
Statistika deskriptif digunakan untuk melihat karakteristik awal dari variabel penelitian. Variabel TPT menjadi fokus utama sebagai variabel dependen, sedangkan RLS, IPM, jumlah penduduk miskin, dan PDRB digunakan sebagai variabel penjelas dalam pemodelan spasial. Nilai minimum dan maksimum menunjukkan rentang ketimpangan antarwilayah, sedangkan standar deviasi menunjukkan seberapa besar variasi data. Jika standar deviasi TPT relatif besar, maka perbedaan tingkat pengangguran antar kabupaten/kota cukup menonjol. Matriks korelasi digunakan sebagai gambaran awal hubungan linear antarvariabel; korelasi yang tinggi antarvariabel independen perlu diperhatikan karena dapat memunculkan indikasi multikolinearitas pada model regresi.
4.4. Peta Deskriptif TPT
peta_tpt <- ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = TPT), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "C", name = "TPT (%)") +
ggplot2::labs(
title = "Sebaran Tingkat Pengangguran Terbuka (TPT)",
subtitle = "Kabupaten/Kota di Provinsi Jawa Tengah",
caption = "Sumber: Data Spasial Vindy.xlsx dan GADM Level 2"
) +
ggplot2::theme_minimal()
peta_tpt
Peta deskriptif menunjukkan variasi nilai TPT antar kabupaten/kota di Provinsi Jawa Tengah. Warna yang semakin pekat menunjukkan wilayah dengan nilai TPT yang semakin tinggi. Pada peta ini, setiap kabupaten/kota sudah diberi label nama wilayah sehingga pembacaan lokasi menjadi lebih jelas. Secara substantif, wilayah dengan warna lebih pekat dapat ditafsirkan sebagai wilayah yang membutuhkan perhatian lebih besar dalam isu ketenagakerjaan, sedangkan wilayah dengan warna lebih muda menunjukkan TPT yang relatif lebih rendah. Peta ini menjadi dasar awal untuk melihat apakah wilayah bernilai tinggi atau rendah tampak berdekatan, sebelum pola tersebut diuji secara statistik melalui ANN, Moran’s I, Geary’s C, dan Getis-Ord.
4.5. Identifikasi Pola Spasial: Acak, Cluster, atau Seragam
Analisis pola spasial dilakukan menggunakan Average Nearest
Neighbor Index (ANNI) berdasarkan titik representatif
kabupaten/kota. Karena data berupa wilayah administratif, titik
representatif yang digunakan adalah st_point_on_surface,
yaitu titik yang berada di dalam setiap wilayah.
# Titik representatif wilayah
jateng_point <- sf::st_point_on_surface(jateng)
coords <- sf::st_coordinates(jateng_point)
# Nearest neighbor k = 1
knn_1 <- spdep::knearneigh(coords, k = 1)
nb_1 <- spdep::knn2nb(knn_1)
dist_nn <- unlist(spdep::nbdists(nb_1, coords))
n <- nrow(jateng)
area_total <- as.numeric(sf::st_area(sf::st_union(jateng)))
mean_obs <- mean(dist_nn)
mean_exp <- 0.5 / sqrt(n / area_total)
se_exp <- 0.26136 / sqrt((n^2) / area_total)
ann_index <- mean_obs / mean_exp
z_ann <- (mean_obs - mean_exp) / se_exp
p_ann <- 2 * pnorm(-abs(z_ann))
pola_spasial <- ifelse(
p_ann < 0.05 & ann_index < 1, "Mengelompok / Cluster",
ifelse(p_ann < 0.05 & ann_index > 1, "Seragam / Teratur", "Acak / Random")
)
ann_table <- data.frame(
Jumlah_Wilayah = n,
Rata_rata_Jarak_Terdekat_Observasi = mean_obs,
Rata_rata_Jarak_Terdekat_Harapan = mean_exp,
ANN_Index = ann_index,
Z_score = z_ann,
P_value = p_ann,
Keputusan_Pola = pola_spasial
)
knitr::kable(round_df(ann_table, 4), caption = "Average Nearest Neighbor Index")| Jumlah_Wilayah | Rata_rata_Jarak_Terdekat_Observasi | Rata_rata_Jarak_Terdekat_Harapan | ANN_Index | Z_score | P_value | Keputusan_Pola |
|---|---|---|---|---|---|---|
| 35 | 24395.33 | 15673.7 | 1.5564 | 6.2978 | 0 | Seragam / Teratur |
Berdasarkan nilai ANN Index sebesar 1.5564 dengan p-value sebesar 0, pola spasial lokasi kabupaten/kota dikategorikan sebagai Seragam / Teratur. Jika nilai indeks kurang dari 1, maka jarak antartitik kabupaten/kota secara rata-rata lebih dekat dibandingkan kondisi acak sehingga menunjukkan kecenderungan mengelompok. Jika nilai indeks lebih dari 1, maka jarak antartitik relatif lebih berjauhan dan cenderung seragam. Jika p-value tidak signifikan, maka pola titik belum cukup kuat untuk dinyatakan berbeda dari pola acak. Hasil ini penting sebagai gambaran awal karena analisis ANN hanya menggunakan posisi wilayah, sedangkan analisis autokorelasi berikutnya mempertimbangkan nilai TPT di setiap wilayah.
ggplot2::ggplot() +
ggplot2::geom_sf(data = jateng, fill = "white", color = "grey40") +
ggplot2::geom_sf(data = jateng_point, color = "red", size = 1.8) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE
) +
ggplot2::labs(
title = "Titik Representatif Kabupaten/Kota Jawa Tengah",
subtitle = "Digunakan untuk analisis pola spasial nearest neighbor"
) +
ggplot2::theme_minimal()
4.6. Matriks Pembobot Spasial
Matriks pembobot utama yang digunakan adalah queen contiguity. Pembobot ini dipilih karena unit analisis berupa wilayah administratif kabupaten/kota. Dua wilayah dianggap bertetangga jika memiliki sisi atau titik sudut yang bersinggungan. Selain itu, pembobot K-nearest neighbors dengan k = 4 juga dibentuk sebagai pembanding berbasis jarak.
# Pembobot queen contiguity
nb_queen <- spdep::poly2nb(jateng, queen = TRUE, row.names = jateng$NAMA)
lw_queen <- spdep::nb2listw(nb_queen, style = "W", zero.policy = TRUE)
lw_queen_b <- spdep::nb2listw(nb_queen, style = "B", zero.policy = TRUE)
# Pembobot rook contiguity sebagai pembanding
nb_rook <- spdep::poly2nb(jateng, queen = FALSE, row.names = jateng$NAMA)
lw_rook <- spdep::nb2listw(nb_rook, style = "W", zero.policy = TRUE)
# Pembobot K-nearest neighbors berbasis jarak centroid/point-on-surface
nb_knn4 <- spdep::knn2nb(spdep::knearneigh(coords, k = 4), row.names = jateng$NAMA)
lw_knn4 <- spdep::nb2listw(nb_knn4, style = "W", zero.policy = TRUE)
ringkasan_tetangga <- data.frame(
Wilayah = jateng$NAMA,
Jumlah_Tetangga_Queen = spdep::card(nb_queen),
Jumlah_Tetangga_Rook = spdep::card(nb_rook),
Jumlah_Tetangga_KNN4 = spdep::card(nb_knn4)
)
knitr::kable(ringkasan_tetangga, caption = "Jumlah Tetangga Berdasarkan Beberapa Matriks Pembobot")| Wilayah | Jumlah_Tetangga_Queen | Jumlah_Tetangga_Rook | Jumlah_Tetangga_KNN4 |
|---|---|---|---|
| KABUPATEN BANJARNEGARA | 6 | 6 | 4 |
| KABUPATEN BANYUMAS | 7 | 7 | 4 |
| KABUPATEN BATANG | 5 | 5 | 4 |
| KABUPATEN BLORA | 3 | 3 | 4 |
| KABUPATEN BOYOLALI | 8 | 8 | 4 |
| KABUPATEN BREBES | 4 | 4 | 4 |
| KABUPATEN CILACAP | 3 | 3 | 4 |
| KABUPATEN DEMAK | 5 | 5 | 4 |
| KABUPATEN GROBOGAN | 7 | 7 | 4 |
| KABUPATEN JEPARA | 3 | 3 | 4 |
| KABUPATEN KARANGANYAR | 5 | 5 | 4 |
| KABUPATEN KEBUMEN | 5 | 5 | 4 |
| KABUPATEN KENDAL | 5 | 5 | 4 |
| KABUPATEN KLATEN | 3 | 2 | 4 |
| KABUPATEN KUDUS | 4 | 4 | 4 |
| KABUPATEN MAGELANG | 7 | 6 | 4 |
| KABUPATEN PATI | 5 | 5 | 4 |
| KABUPATEN PEKALONGAN | 5 | 5 | 4 |
| KABUPATEN PEMALANG | 4 | 4 | 4 |
| KABUPATEN PURBALINGGA | 4 | 4 | 4 |
| KABUPATEN PURWOREJO | 3 | 3 | 4 |
| KABUPATEN REMBANG | 2 | 2 | 4 |
| KABUPATEN SEMARANG | 8 | 8 | 4 |
| KABUPATEN SRAGEN | 3 | 3 | 4 |
| KABUPATEN SUKOHARJO | 5 | 5 | 4 |
| KABUPATEN TEGAL | 4 | 4 | 4 |
| KABUPATEN TEMANGGUNG | 4 | 4 | 4 |
| KABUPATEN WONOGIRI | 2 | 2 | 4 |
| KABUPATEN WONOSOBO | 7 | 7 | 4 |
| KOTA MAGELANG | 1 | 1 | 4 |
| KOTA PEKALONGAN | 2 | 2 | 4 |
| KOTA SALATIGA | 1 | 1 | 4 |
| KOTA SEMARANG | 3 | 3 | 4 |
| KOTA SURAKARTA | 3 | 3 | 4 |
| KOTA TEGAL | 2 | 2 | 4 |
cat("Rata-rata jumlah tetangga Queen:", mean(spdep::card(nb_queen)), "\n")## Rata-rata jumlah tetangga Queen: 4.228571cat("Wilayah tanpa tetangga Queen:", sum(spdep::card(nb_queen) == 0), "\n")## Wilayah tanpa tetangga Queen: 0plot(sf::st_geometry(jateng), border = "grey50", main = "Struktur Tetangga Queen Contiguity")
plot(nb_queen, coords, add = TRUE, col = "red", lwd = 1)
points(coords, pch = 19, col = "blue", cex = 0.7)
text(coords[, 1], coords[, 2], labels = jateng$LABEL, cex = 0.55, pos = 3)
Dalam analisis berikutnya, matriks pembobot yang digunakan adalah queen contiguity dengan standardisasi baris (style = W). Pembobot ini sesuai untuk data areal karena mempertimbangkan hubungan batas antar kabupaten/kota, baik yang bersinggungan melalui sisi maupun titik sudut. Standardisasi baris membuat kontribusi seluruh tetangga pada suatu wilayah berjumlah satu, sehingga lag spasial dapat dibaca sebagai rata-rata kondisi TPT pada wilayah tetangga. Jika suatu wilayah memiliki banyak tetangga, setiap tetangga memperoleh bobot yang lebih kecil; sebaliknya, jika tetangganya sedikit, bobot setiap tetangga menjadi lebih besar.
4.7. Autokorelasi Spasial Global
4.7.1. Moran’s I Global
moran_global <- spdep::moran.test(jateng$TPT, lw_queen, zero.policy = TRUE, alternative = "greater")
moran_global_perm <- spdep::moran.mc(jateng$TPT, lw_queen, nsim = 999, zero.policy = TRUE)
print(moran_global)##
## Moran I test under randomisation
##
## data: jateng$TPT
## weights: lw_queen
##
## Moran I statistic standard deviate = 5.767, p-value = 0.000000004036
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.63175126 -0.02941176 0.01314389print(moran_global_perm)##
## Monte-Carlo simulation of Moran I
##
## data: jateng$TPT
## weights: lw_queen
## number of simulations + 1: 1000
##
## statistic = 0.63175, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterNilai Moran’s I sebesar 0.6318 dengan p-value sebesar 0. Apabila p-value kurang dari 0,05 dan nilai Moran’s I positif, maka TPT memiliki autokorelasi spasial positif yang signifikan. Artinya, kabupaten/kota dengan TPT tinggi cenderung berdekatan dengan kabupaten/kota lain yang juga memiliki TPT tinggi, dan kabupaten/kota dengan TPT rendah cenderung berdekatan dengan wilayah bernilai rendah. Sebaliknya, jika Moran’s I negatif dan signifikan, pola yang muncul adalah ketidaksamaan antarwilayah bertetangga. Jika tidak signifikan, maka variasi TPT belum menunjukkan pola spasial global yang kuat.
4.7.2. Geary’s C Global
geary_global <- spdep::geary.test(jateng$TPT, lw_queen, zero.policy = TRUE, alternative = "greater")
print(geary_global)##
## Geary C test under randomisation
##
## data: jateng$TPT
## weights: lw_queen
##
## Geary C statistic standard deviate = 4.8995, p-value = 0.0000004804
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.34757812 1.00000000 0.01773194Nilai Geary’s C sebesar 0.3476 dengan p-value sebesar 0. Nilai Geary’s C kurang dari 1 menunjukkan autokorelasi spasial positif, sedangkan nilai lebih dari 1 menunjukkan autokorelasi spasial negatif. Karena Geary’s C menekankan selisih nilai antarwilayah bertetangga, hasil ini melengkapi Moran’s I. Jika hasilnya signifikan dan kurang dari 1, maka wilayah bertetangga cenderung memiliki nilai TPT yang mirip. Jika lebih dari 1 dan signifikan, maka wilayah bertetangga cenderung memiliki perbedaan TPT yang besar.
4.7.3. Getis-Ord G Global
# Getis-Ord G lebih sesuai menggunakan pembobot biner.
getis_global <- spdep::globalG.test(jateng$TPT, lw_queen_b, zero.policy = TRUE, alternative = "greater")
print(getis_global)##
## Getis-Ord global G statistic
##
## data: jateng$TPT
## weights: lw_queen_b
##
## standard deviate = 1.0011, p-value = 0.1584
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 0.13059241499 0.12436974790 0.00003863361Getis-Ord G global digunakan untuk melihat apakah terdapat konsentrasi nilai TPT tinggi secara spasial. Jika p-value signifikan dan statistik mengarah positif, maka terdapat indikasi bahwa nilai TPT tinggi terkonsentrasi pada wilayah yang saling berdekatan. Jika statistik mengarah negatif, maka konsentrasi yang lebih menonjol adalah nilai rendah. Statistik ini berbeda dari Moran’s I karena lebih fokus pada keberadaan hot spot atau konsentrasi nilai tinggi/rendah, bukan sekadar kemiripan nilai antarwilayah.
4.7.4. Ringkasan Autokorelasi Global
ringkasan_global <- data.frame(
Statistik = c("Moran's I", "Geary's C", "Getis-Ord G"),
Nilai = c(
as.numeric(moran_global$estimate[1]),
as.numeric(geary_global$estimate[1]),
as.numeric(getis_global$estimate[1])
),
Z_score = c(
as.numeric(moran_global$statistic),
as.numeric(geary_global$statistic),
as.numeric(getis_global$statistic)
),
P_value = c(
moran_global$p.value,
geary_global$p.value,
getis_global$p.value
)
)
knitr::kable(round_df(ringkasan_global, 5), caption = "Ringkasan Autokorelasi Spasial Global")| Statistik | Nilai | Z_score | P_value |
|---|---|---|---|
| Moran’s I | 0.63175 | 5.76696 | 0.00000 |
| Geary’s C | 0.34758 | 4.89948 | 0.00000 |
| Getis-Ord G | 0.13059 | 1.00114 | 0.15838 |
4.8. Autokorelasi Spasial Lokal
4.8.1. Moran’s I Lokal / LISA
lisa <- spdep::localmoran(jateng$TPT, lw_queen, zero.policy = TRUE)
jateng$Local_Moran_I <- lisa[, 1]
jateng$Local_Moran_Z <- lisa[, 4]
jateng$P_Moran <- lisa[, 5]
mean_tpt <- mean(jateng$TPT, na.rm = TRUE)
lag_tpt <- spdep::lag.listw(lw_queen, jateng$TPT, zero.policy = TRUE)
jateng$Cluster_Moran <- "Tidak Signifikan"
jateng$Cluster_Moran[jateng$TPT >= mean_tpt & lag_tpt >= mean_tpt & jateng$P_Moran <= 0.05] <- "High-High"
jateng$Cluster_Moran[jateng$TPT < mean_tpt & lag_tpt < mean_tpt & jateng$P_Moran <= 0.05] <- "Low-Low"
jateng$Cluster_Moran[jateng$TPT >= mean_tpt & lag_tpt < mean_tpt & jateng$P_Moran <= 0.05] <- "High-Low"
jateng$Cluster_Moran[jateng$TPT < mean_tpt & lag_tpt >= mean_tpt & jateng$P_Moran <= 0.05] <- "Low-High"
lisa_table <- sf::st_drop_geometry(jateng) %>%
dplyr::select(NAMA, TPT, Local_Moran_I, Local_Moran_Z, P_Moran, Cluster_Moran) %>%
dplyr::arrange(P_Moran)
knitr::kable(round_df(lisa_table, 4), caption = "Hasil Moran's I Lokal / LISA")| NAMA | TPT | Local_Moran_I | Local_Moran_Z | P_Moran | Cluster_Moran |
|---|---|---|---|---|---|
| KABUPATEN BANYUMAS | 6.18 | 1.5633 | 4.1729 | 0.0000 | High-High |
| KABUPATEN BREBES | 8.35 | 4.1223 | 3.8385 | 0.0001 | High-High |
| KABUPATEN TEGAL | 7.53 | 3.1080 | 3.4899 | 0.0005 | High-High |
| KOTA TEGAL | 5.88 | 2.1442 | 3.3676 | 0.0008 | High-High |
| KABUPATEN CILACAP | 7.83 | 3.0719 | 2.7354 | 0.0062 | High-High |
| KABUPATEN PATI | 3.87 | 0.3790 | 2.0565 | 0.0397 | Low-Low |
| KABUPATEN BOYOLALI | 3.16 | 0.5027 | 1.8186 | 0.0690 | Tidak Signifikan |
| KABUPATEN MAGELANG | 3.55 | 0.3906 | 1.7777 | 0.0755 | Tidak Signifikan |
| KABUPATEN KARANGANYAR | 3.47 | 0.5078 | 1.7376 | 0.0823 | Tidak Signifikan |
| KABUPATEN GROBOGAN | 3.23 | 0.4872 | 1.6978 | 0.0896 | Tidak Signifikan |
| KABUPATEN SUKOHARJO | 3.65 | 0.3999 | 1.6416 | 0.1007 | Tidak Signifikan |
| KABUPATEN PEKALONGAN | 3.30 | -0.5780 | -1.6001 | 0.1096 | Tidak Signifikan |
| KABUPATEN KEBUMEN | 5.07 | 0.2481 | 1.5904 | 0.1117 | Tidak Signifikan |
| KABUPATEN SRAGEN | 3.53 | 0.5624 | 1.5193 | 0.1287 | Tidak Signifikan |
| KABUPATEN PEMALANG | 6.63 | 0.9461 | 1.4984 | 0.1340 | Tidak Signifikan |
| KABUPATEN BLORA | 3.67 | 0.4724 | 1.4796 | 0.1390 | Tidak Signifikan |
| KOTA SURAKARTA | 4.61 | -0.0455 | -1.3039 | 0.1923 | Tidak Signifikan |
| KABUPATEN KLATEN | 3.97 | 0.2702 | 1.2971 | 0.1946 | Tidak Signifikan |
| KABUPATEN PURBALINGGA | 4.96 | 0.1827 | 1.2834 | 0.1993 | Tidak Signifikan |
| KABUPATEN SEMARANG | 3.73 | 0.2014 | 1.2348 | 0.2169 | Tidak Signifikan |
| KABUPATEN KUDUS | 3.19 | 0.4426 | 1.0816 | 0.2794 | Tidak Signifikan |
| KABUPATEN DEMAK | 4.75 | -0.0698 | -1.0393 | 0.2987 | Tidak Signifikan |
| KABUPATEN WONOGIRI | 2.40 | 0.9375 | 1.0135 | 0.3108 | Tidak Signifikan |
| KABUPATEN REMBANG | 2.84 | 0.5804 | 0.7830 | 0.4336 | Tidak Signifikan |
| KABUPATEN JEPARA | 3.34 | 0.3170 | 0.7457 | 0.4558 | Tidak Signifikan |
| KABUPATEN TEMANGGUNG | 2.35 | 0.4422 | 0.7340 | 0.4629 | Tidak Signifikan |
| KABUPATEN BANJARNEGARA | 5.57 | 0.1679 | 0.6774 | 0.4982 | Tidak Signifikan |
| KOTA MAGELANG | 4.40 | 0.0535 | 0.6513 | 0.5149 | Tidak Signifikan |
| KOTA SALATIGA | 3.86 | 0.2401 | 0.5430 | 0.5871 | Tidak Signifikan |
| KABUPATEN PURWOREJO | 3.89 | 0.0889 | 0.3895 | 0.6969 | Tidak Signifikan |
| KABUPATEN KENDAL | 5.01 | -0.0456 | -0.2994 | 0.7646 | Tidak Signifikan |
| KABUPATEN WONOSOBO | 4.02 | 0.0174 | 0.1767 | 0.8598 | Tidak Signifikan |
| KABUPATEN BATANG | 5.67 | 0.0224 | 0.1226 | 0.9024 | Tidak Signifikan |
| KOTA PEKALONGAN | 4.91 | -0.0062 | -0.0224 | 0.9822 | Tidak Signifikan |
| KOTA SEMARANG | 5.82 | -0.0138 | 0.0184 | 0.9853 | Tidak Signifikan |
table(jateng$Cluster_Moran)##
## High-High Low-Low Tidak Signifikan
## 5 1 29peta_lisa <- ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = Cluster_Moran), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_manual(
values = c(
"High-High" = "red",
"Low-Low" = "blue",
"High-Low" = "orange",
"Low-High" = "skyblue",
"Tidak Signifikan" = "grey85"
),
name = "Cluster LISA"
) +
ggplot2::labs(
title = "Peta Local Moran's I / LISA",
subtitle = "Klasifikasi cluster lokal TPT, alpha = 0,05"
) +
ggplot2::theme_minimal()
peta_lisa
LISA digunakan untuk mengidentifikasi wilayah yang memiliki pola lokal signifikan. Kategori High-High menunjukkan wilayah dengan TPT tinggi yang dikelilingi wilayah bertetangga dengan TPT tinggi, sehingga dapat disebut sebagai klaster TPT tinggi. Kategori Low-Low menunjukkan wilayah dengan TPT rendah yang dikelilingi wilayah bertetangga dengan TPT rendah, sehingga dapat disebut sebagai klaster TPT rendah. Kategori High-Low dan Low-High menunjukkan wilayah yang berbeda dari lingkungan sekitarnya, atau dapat dipandang sebagai indikasi spatial outlier. Wilayah yang tidak signifikan berarti nilai TPT wilayah tersebut belum menunjukkan hubungan lokal yang cukup kuat dengan tetangganya pada taraf signifikansi 5%.
4.8.2. Geary’s C Lokal
# Local Geary's C
local_c_obs <- as.numeric(spdep::localC(jateng$TPT, lw_queen, zero.policy = TRUE))
# P-value pseudo-permutation sederhana untuk local Geary.
# Nilai local C yang besar menunjukkan perbedaan lokal yang kuat.
set.seed(123)
nsim_c <- 999
perm_local_c <- replicate(
nsim_c,
as.numeric(spdep::localC(sample(jateng$TPT), lw_queen, zero.policy = TRUE))
)
p_local_c <- sapply(seq_along(local_c_obs), function(i) {
valid_perm <- perm_local_c[i, !is.na(perm_local_c[i, ])]
(sum(valid_perm >= local_c_obs[i], na.rm = TRUE) + 1) / (length(valid_perm) + 1)
})
jateng$Local_Geary_C <- local_c_obs
jateng$P_Geary_Local <- p_local_c
jateng$Cluster_Geary <- ifelse(jateng$P_Geary_Local <= 0.05, "Signifikan", "Tidak Signifikan")
geary_local_table <- sf::st_drop_geometry(jateng) %>%
dplyr::select(NAMA, TPT, Local_Geary_C, P_Geary_Local, Cluster_Geary) %>%
dplyr::arrange(P_Geary_Local)
knitr::kable(round_df(geary_local_table, 4), caption = "Hasil Geary's C Lokal")| NAMA | TPT | Local_Geary_C | P_Geary_Local | Cluster_Geary |
|---|---|---|---|---|
| KABUPATEN PEKALONGAN | 3.30 | 2.4360 | 0.264 | Tidak Signifikan |
| KOTA TEGAL | 5.88 | 1.9751 | 0.357 | Tidak Signifikan |
| KABUPATEN PEMALANG | 6.63 | 1.6666 | 0.423 | Tidak Signifikan |
| KABUPATEN CILACAP | 7.83 | 1.5834 | 0.443 | Tidak Signifikan |
| KABUPATEN TEMANGGUNG | 2.35 | 1.4784 | 0.481 | Tidak Signifikan |
| KABUPATEN BREBES | 8.35 | 1.3154 | 0.559 | Tidak Signifikan |
| KOTA SEMARANG | 5.82 | 0.9206 | 0.606 | Tidak Signifikan |
| KABUPATEN KEBUMEN | 5.07 | 1.0382 | 0.649 | Tidak Signifikan |
| KOTA MAGELANG | 4.40 | 0.3235 | 0.666 | Tidak Signifikan |
| KABUPATEN KENDAL | 5.01 | 0.9657 | 0.673 | Tidak Signifikan |
| KOTA PEKALONGAN | 4.91 | 0.7095 | 0.677 | Tidak Signifikan |
| KABUPATEN WONOGIRI | 2.40 | 0.6060 | 0.693 | Tidak Signifikan |
| KABUPATEN PURBALINGGA | 4.96 | 0.8288 | 0.696 | Tidak Signifikan |
| KOTA SURAKARTA | 4.61 | 0.6452 | 0.709 | Tidak Signifikan |
| KABUPATEN BATANG | 5.67 | 0.8383 | 0.725 | Tidak Signifikan |
| KABUPATEN DEMAK | 4.75 | 0.7985 | 0.746 | Tidak Signifikan |
| KABUPATEN TEGAL | 7.53 | 0.6746 | 0.761 | Tidak Signifikan |
| KABUPATEN REMBANG | 2.84 | 0.3917 | 0.767 | Tidak Signifikan |
| KABUPATEN BANYUMAS | 6.18 | 0.8026 | 0.784 | Tidak Signifikan |
| KABUPATEN WONOSOBO | 4.02 | 0.6545 | 0.846 | Tidak Signifikan |
| KABUPATEN BANJARNEGARA | 5.57 | 0.6387 | 0.848 | Tidak Signifikan |
| KABUPATEN JEPARA | 3.34 | 0.3420 | 0.872 | Tidak Signifikan |
| KABUPATEN KUDUS | 3.19 | 0.3268 | 0.910 | Tidak Signifikan |
| KABUPATEN SEMARANG | 3.73 | 0.5359 | 0.911 | Tidak Signifikan |
| KABUPATEN PURWOREJO | 3.89 | 0.2276 | 0.917 | Tidak Signifikan |
| KOTA SALATIGA | 3.86 | 0.0076 | 0.947 | Tidak Signifikan |
| KABUPATEN BLORA | 3.67 | 0.1377 | 0.948 | Tidak Signifikan |
| KABUPATEN KLATEN | 3.97 | 0.1395 | 0.959 | Tidak Signifikan |
| KABUPATEN SUKOHARJO | 3.65 | 0.2560 | 0.961 | Tidak Signifikan |
| KABUPATEN KARANGANYAR | 3.47 | 0.2307 | 0.967 | Tidak Signifikan |
| KABUPATEN PATI | 3.87 | 0.2018 | 0.973 | Tidak Signifikan |
| KABUPATEN SRAGEN | 3.53 | 0.0344 | 0.990 | Tidak Signifikan |
| KABUPATEN GROBOGAN | 3.23 | 0.2085 | 0.993 | Tidak Signifikan |
| KABUPATEN MAGELANG | 3.55 | 0.1829 | 0.994 | Tidak Signifikan |
| KABUPATEN BOYOLALI | 3.16 | 0.2078 | 0.999 | Tidak Signifikan |
table(jateng$Cluster_Geary)##
## Tidak Signifikan
## 35ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = Cluster_Geary), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_manual(
values = c("Signifikan" = "purple", "Tidak Signifikan" = "grey85"),
name = "Local Geary"
) +
ggplot2::labs(
title = "Peta Signifikansi Local Geary's C",
subtitle = "Wilayah signifikan menunjukkan ketidakmiripan lokal yang relatif kuat"
) +
ggplot2::theme_minimal()
Local Geary’s C digunakan untuk melihat wilayah yang memiliki perbedaan lokal kuat dibandingkan wilayah sekitarnya. Wilayah signifikan menunjukkan adanya indikasi ketidakmiripan spasial lokal, yaitu nilai TPT wilayah tersebut relatif berbeda dibandingkan nilai TPT tetangganya. Dengan demikian, peta Local Geary berguna untuk menandai wilayah yang berpotensi menjadi pencilan lokal. Jika suatu wilayah signifikan pada Local Geary tetapi tidak termasuk klaster LISA High-High atau Low-Low, maka wilayah tersebut perlu ditinjau lebih lanjut karena karakteristiknya dapat berbeda dari pola umum di sekitarnya.
4.8.3. Getis-Ord Lokal
local_getis_z <- as.numeric(spdep::localG(jateng$TPT, lw_queen_b, zero.policy = TRUE))
local_getis_p <- 2 * pnorm(-abs(local_getis_z))
jateng$Getis_Z <- local_getis_z
jateng$P_Getis <- local_getis_p
jateng$Cluster_Getis <- "Tidak Signifikan"
jateng$Cluster_Getis[jateng$Getis_Z > 0 & jateng$P_Getis <= 0.05] <- "Hot Spot"
jateng$Cluster_Getis[jateng$Getis_Z < 0 & jateng$P_Getis <= 0.05] <- "Cold Spot"
getis_local_table <- sf::st_drop_geometry(jateng) %>%
dplyr::select(NAMA, TPT, Getis_Z, P_Getis, Cluster_Getis) %>%
dplyr::arrange(P_Getis)
knitr::kable(round_df(getis_local_table, 4), caption = "Hasil Getis-Ord Lokal")| NAMA | TPT | Getis_Z | P_Getis | Cluster_Getis |
|---|---|---|---|---|
| KABUPATEN BANYUMAS | 6.18 | 4.1729 | 0.0000 | Hot Spot |
| KABUPATEN BREBES | 8.35 | 3.8385 | 0.0001 | Hot Spot |
| KABUPATEN TEGAL | 7.53 | 3.4899 | 0.0005 | Hot Spot |
| KOTA TEGAL | 5.88 | 3.3676 | 0.0008 | Hot Spot |
| KABUPATEN CILACAP | 7.83 | 2.7354 | 0.0062 | Hot Spot |
| KABUPATEN PATI | 3.87 | -2.0565 | 0.0397 | Cold Spot |
| KABUPATEN BOYOLALI | 3.16 | -1.8186 | 0.0690 | Tidak Signifikan |
| KABUPATEN MAGELANG | 3.55 | -1.7777 | 0.0755 | Tidak Signifikan |
| KABUPATEN KARANGANYAR | 3.47 | -1.7376 | 0.0823 | Tidak Signifikan |
| KABUPATEN GROBOGAN | 3.23 | -1.6978 | 0.0896 | Tidak Signifikan |
| KABUPATEN SUKOHARJO | 3.65 | -1.6416 | 0.1007 | Tidak Signifikan |
| KABUPATEN PEKALONGAN | 3.30 | 1.6001 | 0.1096 | Tidak Signifikan |
| KABUPATEN KEBUMEN | 5.07 | 1.5904 | 0.1117 | Tidak Signifikan |
| KABUPATEN SRAGEN | 3.53 | -1.5193 | 0.1287 | Tidak Signifikan |
| KABUPATEN PEMALANG | 6.63 | 1.4984 | 0.1340 | Tidak Signifikan |
| KABUPATEN BLORA | 3.67 | -1.4796 | 0.1390 | Tidak Signifikan |
| KOTA SURAKARTA | 4.61 | -1.3039 | 0.1923 | Tidak Signifikan |
| KABUPATEN KLATEN | 3.97 | -1.2971 | 0.1946 | Tidak Signifikan |
| KABUPATEN PURBALINGGA | 4.96 | 1.2834 | 0.1993 | Tidak Signifikan |
| KABUPATEN SEMARANG | 3.73 | -1.2348 | 0.2169 | Tidak Signifikan |
| KABUPATEN KUDUS | 3.19 | -1.0816 | 0.2794 | Tidak Signifikan |
| KABUPATEN DEMAK | 4.75 | -1.0393 | 0.2987 | Tidak Signifikan |
| KABUPATEN WONOGIRI | 2.40 | -1.0135 | 0.3108 | Tidak Signifikan |
| KABUPATEN REMBANG | 2.84 | -0.7830 | 0.4336 | Tidak Signifikan |
| KABUPATEN JEPARA | 3.34 | -0.7457 | 0.4558 | Tidak Signifikan |
| KABUPATEN TEMANGGUNG | 2.35 | -0.7340 | 0.4629 | Tidak Signifikan |
| KABUPATEN BANJARNEGARA | 5.57 | 0.6774 | 0.4982 | Tidak Signifikan |
| KOTA MAGELANG | 4.40 | -0.6513 | 0.5149 | Tidak Signifikan |
| KOTA SALATIGA | 3.86 | -0.5430 | 0.5871 | Tidak Signifikan |
| KABUPATEN PURWOREJO | 3.89 | -0.3895 | 0.6969 | Tidak Signifikan |
| KABUPATEN KENDAL | 5.01 | -0.2994 | 0.7646 | Tidak Signifikan |
| KABUPATEN WONOSOBO | 4.02 | -0.1767 | 0.8598 | Tidak Signifikan |
| KABUPATEN BATANG | 5.67 | 0.1226 | 0.9024 | Tidak Signifikan |
| KOTA PEKALONGAN | 4.91 | -0.0224 | 0.9822 | Tidak Signifikan |
| KOTA SEMARANG | 5.82 | 0.0184 | 0.9853 | Tidak Signifikan |
table(jateng$Cluster_Getis)##
## Cold Spot Hot Spot Tidak Signifikan
## 1 5 29ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = Cluster_Getis), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_manual(
values = c(
"Hot Spot" = "red",
"Cold Spot" = "blue",
"Tidak Signifikan" = "grey85"
),
name = "Getis-Ord"
) +
ggplot2::labs(
title = "Peta Getis-Ord Lokal",
subtitle = "Hot spot dan cold spot TPT, alpha = 0,05"
) +
ggplot2::theme_minimal()
Getis-Ord lokal mengidentifikasi konsentrasi nilai tinggi dan rendah. Hot spot menunjukkan wilayah dengan konsentrasi TPT tinggi di sekitar wilayah tersebut, sedangkan cold spot menunjukkan konsentrasi TPT rendah. Dalam konteks kebijakan ketenagakerjaan, wilayah hot spot dapat diprioritaskan untuk intervensi pengurangan pengangguran karena wilayah tersebut tidak hanya memiliki nilai tinggi secara lokal, tetapi juga berada dalam lingkungan tetangga yang relatif tinggi. Cold spot menunjukkan kelompok wilayah dengan kondisi TPT relatif rendah dan dapat dikaji sebagai pembanding untuk memahami faktor yang mendukung rendahnya pengangguran.
4.9. Kriging Global dan Kriging Lokal
Kriging dilakukan menggunakan titik representatif kabupaten/kota.
Ordinary Kriging global menggunakan seluruh titik pengamatan, sedangkan
Kriging lokal menggunakan maksimal 8 titik tetangga terdekat
(nmax = 8).
4.9.1. Persiapan Titik dan Grid Prediksi
# Titik pengamatan untuk kriging
points_sf <- sf::st_point_on_surface(jateng)
points_sp <- as(points_sf, "Spatial")
# Grid prediksi dalam wilayah Jawa Tengah
cellsize <- 5000 # meter; dapat diperkecil untuk peta lebih halus, tetapi komputasi lebih lama
boundary_jateng <- sf::st_union(jateng)
grid_sfc <- sf::st_make_grid(boundary_jateng, cellsize = cellsize, what = "centers")
grid_sf <- sf::st_sf(id = seq_along(grid_sfc), geometry = grid_sfc, crs = sf::st_crs(jateng))
grid_sf <- sf::st_filter(grid_sf, boundary_jateng, .predicate = sf::st_intersects)
grid_sp <- as(grid_sf, "Spatial")
cat("Jumlah titik grid prediksi:", nrow(grid_sf), "\n")## Jumlah titik grid prediksi: 13834.9.2. Variogram Empiris dan Fitting Variogram
vario_emp <- gstat::variogram(TPT ~ 1, points_sp)
print(vario_emp)## np dist gamma dir.hor dir.ver id
## 1 1 4544.197 0.3612500 0 0 var1
## 2 2 13277.745 0.2346250 0 0 var1
## 3 2 19128.100 1.3286500 0 0 var1
## 4 24 26560.347 0.9184500 0 0 var1
## 5 28 33422.476 0.8873429 0 0 var1
## 6 21 40133.908 0.9633238 0 0 var1
## 7 19 48519.816 0.5361000 0 0 var1
## 8 35 55142.637 1.0912500 0 0 var1
## 9 27 62534.752 1.2578278 0 0 var1
## 10 26 68974.863 1.0308442 0 0 var1
## 11 22 75938.070 0.9561250 0 0 var1
## 12 38 84254.783 1.3208882 0 0 var1
## 13 31 90685.718 1.1827710 0 0 var1
## 14 29 97976.854 1.2206500 0 0 var1
## 15 19 105703.743 1.2157079 0 0 var1initial_psill <- stats::var(points_sf$TPT, na.rm = TRUE)
initial_range <- max(sp::spDists(sp::coordinates(points_sp))) / 3
fit_variogram_safely <- function(model_name) {
initial_model <- gstat::vgm(
psill = initial_psill,
model = model_name,
range = initial_range,
nugget = 0
)
suppressWarnings(
tryCatch(
gstat::fit.variogram(vario_emp, initial_model),
error = function(e) NULL
)
)
}
vario_models <- list(
Spherical = fit_variogram_safely("Sph"),
Exponential = fit_variogram_safely("Exp"),
Gaussian = fit_variogram_safely("Gau")
)
sse_models <- sapply(vario_models, function(m) {
if (is.null(m)) Inf else attr(m, "SSErr")
})
best_vario_name <- names(which.min(sse_models))
vario_model <- vario_models[[best_vario_name]]
if (is.null(vario_model) || !is.finite(min(sse_models))) {
best_vario_name <- "Spherical_Default"
vario_model <- gstat::vgm(psill = initial_psill, model = "Sph", range = initial_range, nugget = 0)
}
cat("Model variogram terbaik:", best_vario_name, "\n")## Model variogram terbaik: Exponentialprint(vario_model)## model psill range
## 1 Nug 0.1798728 0.00
## 2 Exp 0.9790047 24044.99plot(vario_emp, vario_model, main = paste("Variogram Empiris dan Model", best_vario_name))
4.9.3. Ordinary Kriging Global
kriging_global <- gstat::krige(
formula = TPT ~ 1,
locations = points_sp,
newdata = grid_sp,
model = vario_model
)## [using ordinary kriging]kriging_global_sf <- sf::st_as_sf(kriging_global)
kriging_global_df <- data.frame(
sf::st_coordinates(kriging_global_sf),
sf::st_drop_geometry(kriging_global_sf)
)
summary_kriging_global <- data.frame(
Statistik = c("Minimum", "Mean", "Maximum", "Mean Variance"),
Nilai = c(
min(kriging_global_df$var1.pred, na.rm = TRUE),
mean(kriging_global_df$var1.pred, na.rm = TRUE),
max(kriging_global_df$var1.pred, na.rm = TRUE),
mean(kriging_global_df$var1.var, na.rm = TRUE)
)
)
knitr::kable(round_df(summary_kriging_global, 4), caption = "Ringkasan Ordinary Kriging Global")| Statistik | Nilai |
|---|---|
| Minimum | 2.8320 |
| Mean | 4.5772 |
| Maximum | 7.8274 |
| Mean Variance | 0.7909 |
ggplot2::ggplot() +
ggplot2::geom_raster(
data = kriging_global_df,
ggplot2::aes(x = X, y = Y, fill = var1.pred)
) +
ggplot2::geom_sf(data = jateng, fill = NA, color = "black", linewidth = 0.25) +
ggplot2::geom_sf(data = points_sf, color = "red", size = 1) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "C", name = "Prediksi TPT") +
ggplot2::coord_sf(crs = sf::st_crs(jateng)) +
ggplot2::labs(
title = "Ordinary Kriging Global - Prediksi TPT",
subtitle = "Menggunakan seluruh titik pengamatan"
) +
ggplot2::theme_minimal()
ggplot2::ggplot() +
ggplot2::geom_raster(
data = kriging_global_df,
ggplot2::aes(x = X, y = Y, fill = var1.var)
) +
ggplot2::geom_sf(data = jateng, fill = NA, color = "black", linewidth = 0.25) +
ggplot2::geom_sf(data = points_sf, color = "blue", size = 1) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "A", name = "Variance") +
ggplot2::coord_sf(crs = sf::st_crs(jateng)) +
ggplot2::labs(
title = "Ordinary Kriging Global - Variance",
subtitle = "Varians menunjukkan ketidakpastian prediksi"
) +
ggplot2::theme_minimal()
4.9.4. Ordinary Kriging Lokal
kriging_local <- gstat::krige(
formula = TPT ~ 1,
locations = points_sp,
newdata = grid_sp,
model = vario_model,
nmax = 8
)## [using ordinary kriging]kriging_local_sf <- sf::st_as_sf(kriging_local)
kriging_local_df <- data.frame(
sf::st_coordinates(kriging_local_sf),
sf::st_drop_geometry(kriging_local_sf)
)
summary_kriging_local <- data.frame(
Statistik = c("Minimum", "Mean", "Maximum", "Mean Variance"),
Nilai = c(
min(kriging_local_df$var1.pred, na.rm = TRUE),
mean(kriging_local_df$var1.pred, na.rm = TRUE),
max(kriging_local_df$var1.pred, na.rm = TRUE),
mean(kriging_local_df$var1.var, na.rm = TRUE)
)
)
knitr::kable(round_df(summary_kriging_local, 4), caption = "Ringkasan Ordinary Kriging Lokal")| Statistik | Nilai |
|---|---|
| Minimum | 2.6642 |
| Mean | 4.5769 |
| Maximum | 7.9856 |
| Mean Variance | 0.8043 |
ggplot2::ggplot() +
ggplot2::geom_raster(
data = kriging_local_df,
ggplot2::aes(x = X, y = Y, fill = var1.pred)
) +
ggplot2::geom_sf(data = jateng, fill = NA, color = "black", linewidth = 0.25) +
ggplot2::geom_sf(data = points_sf, color = "red", size = 1) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "C", name = "Prediksi TPT") +
ggplot2::coord_sf(crs = sf::st_crs(jateng)) +
ggplot2::labs(
title = "Ordinary Kriging Lokal - Prediksi TPT",
subtitle = "Menggunakan maksimum 8 tetangga terdekat"
) +
ggplot2::theme_minimal()
ggplot2::ggplot() +
ggplot2::geom_raster(
data = kriging_local_df,
ggplot2::aes(x = X, y = Y, fill = var1.var)
) +
ggplot2::geom_sf(data = jateng, fill = NA, color = "black", linewidth = 0.25) +
ggplot2::geom_sf(data = points_sf, color = "blue", size = 1) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "A", name = "Variance") +
ggplot2::coord_sf(crs = sf::st_crs(jateng)) +
ggplot2::labs(
title = "Ordinary Kriging Lokal - Variance",
subtitle = "Varians menunjukkan ketidakpastian prediksi lokal"
) +
ggplot2::theme_minimal()
perbandingan_kriging <- data.frame(
Metode = c("Kriging Global", "Kriging Lokal nmax=8"),
Minimum = c(min(kriging_global_df$var1.pred, na.rm = TRUE), min(kriging_local_df$var1.pred, na.rm = TRUE)),
Mean = c(mean(kriging_global_df$var1.pred, na.rm = TRUE), mean(kriging_local_df$var1.pred, na.rm = TRUE)),
Maximum = c(max(kriging_global_df$var1.pred, na.rm = TRUE), max(kriging_local_df$var1.pred, na.rm = TRUE)),
Mean_Variance = c(mean(kriging_global_df$var1.var, na.rm = TRUE), mean(kriging_local_df$var1.var, na.rm = TRUE))
)
knitr::kable(round_df(perbandingan_kriging, 4), caption = "Perbandingan Kriging Global dan Lokal")| Metode | Minimum | Mean | Maximum | Mean_Variance |
|---|---|---|---|---|
| Kriging Global | 2.8320 | 4.5772 | 7.8274 | 0.7909 |
| Kriging Lokal nmax=8 | 2.6642 | 4.5769 | 7.9856 | 0.8043 |
Kriging global menghasilkan permukaan prediksi yang lebih halus karena menggunakan seluruh titik pengamatan pada proses estimasi. Oleh karena itu, hasil global lebih stabil dan cocok untuk melihat kecenderungan umum TPT di Jawa Tengah. Kriging lokal lebih sensitif terhadap variasi setempat karena hanya menggunakan sejumlah tetangga terdekat, sehingga perubahan antarwilayah dapat terlihat lebih tajam. Peta varians digunakan untuk melihat wilayah dengan ketidakpastian prediksi yang relatif lebih tinggi. Varians yang besar menunjukkan bahwa prediksi di area tersebut lebih tidak pasti, biasanya karena jaraknya relatif jauh dari titik pengamatan atau pola variogram kurang mampu menjelaskan variasi lokal pada area tersebut.
4.10. Pemodelan Spasial Menggunakan SAR
Pada bagian ini digunakan model SAR karena model tersebut memasukkan
pengaruh TPT wilayah tetangga melalui komponen lag spasial
WY. Sebelum membentuk model SAR, model OLS digunakan
sebagai pembanding awal.
4.10.1. Model OLS
model_ols <- lm(TPT ~ RLS + IPM + Miskin + PDRB, data = sf::st_drop_geometry(jateng))
summary(model_ols)##
## Call:
## lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data = sf::st_drop_geometry(jateng))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4211 -0.6294 -0.1513 0.7012 3.1865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.36682855 12.18160827 1.261 0.21686
## RLS 0.72849132 0.77370633 0.942 0.35394
## IPM -0.25848784 0.24144307 -1.071 0.29289
## Miskin 0.01459710 0.00501563 2.910 0.00674 **
## PDRB 0.00001672 0.00001101 1.518 0.13950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.329 on 30 degrees of freedom
## Multiple R-squared: 0.3027, Adjusted R-squared: 0.2097
## F-statistic: 3.256 on 4 and 30 DF, p-value: 0.02483# Cek multikolinearitas
car::vif(model_ols)## RLS IPM Miskin PDRB
## 21.751065 22.755978 1.711626 2.573365# Uji Lagrange Multiplier / Rao's Score untuk indikasi model spasial.
# Perbaikan: beberapa versi spdep tidak menerima daftar test manual.
# Gunakan test = "all" agar kompatibel lintas versi spdep.
lm_spatial_test <- tryCatch({
if ("lm.RStests" %in% getNamespaceExports("spdep")) {
spdep::lm.RStests(
model_ols,
listw = lw_queen,
test = "all",
zero.policy = TRUE
)
} else {
spdep::lm.LMtests(
model_ols,
listw = lw_queen,
test = "all",
zero.policy = TRUE
)
}
}, error = function(e) {
cat("Uji LM/Rao's Score tidak dapat dijalankan pada versi spdep ini.
")
cat("Pesan error:", e$message, "
")
NULL
})
if (!is.null(lm_spatial_test)) {
print(lm_spatial_test)
}##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data =
## sf::st_drop_geometry(jateng))
## test weights: lw_queen
##
## RSerr = 13.108, df = 1, p-value = 0.000294
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data =
## sf::st_drop_geometry(jateng))
## test weights: lw_queen
##
## RSlag = 21.678, df = 1, p-value = 0.000003224
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data =
## sf::st_drop_geometry(jateng))
## test weights: lw_queen
##
## adjRSerr = 1.343, df = 1, p-value = 0.2465
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data =
## sf::st_drop_geometry(jateng))
## test weights: lw_queen
##
## adjRSlag = 9.9134, df = 1, p-value = 0.001641
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = TPT ~ RLS + IPM + Miskin + PDRB, data =
## sf::st_drop_geometry(jateng))
## test weights: lw_queen
##
## SARMA = 23.021, df = 2, p-value = 0.00001002Model OLS digunakan sebagai model dasar tanpa mempertimbangkan efek spasial. Jika uji lag spasial signifikan, maka model SAR dapat digunakan untuk menangkap pengaruh TPT wilayah tetangga.
4.10.2. Model SAR
model_sar <- spatialreg::lagsarlm(
TPT ~ RLS + IPM + Miskin + PDRB,
data = sf::st_drop_geometry(jateng),
listw = lw_queen,
zero.policy = TRUE,
method = "eigen"
)
summary(model_sar)##
## Call:spatialreg::lagsarlm(formula = TPT ~ RLS + IPM + Miskin + PDRB,
## data = sf::st_drop_geometry(jateng), listw = lw_queen, method = "eigen",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.42091 -0.57402 -0.11877 0.61847 1.80233
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4868763672 7.6962406342 -0.0633 0.949558
## RLS 0.1630976682 0.4852828730 0.3361 0.736805
## IPM -0.0143924195 0.1515216617 -0.0950 0.924326
## Miskin 0.0104549395 0.0031948248 3.2725 0.001066
## PDRB 0.0000090320 0.0000068993 1.3091 0.190498
##
## Rho: 0.70266, LR test value: 22.189, p-value: 0.0000024711
## Asymptotic standard error: 0.11157
## z-value: 6.2978, p-value: 0.00000000030196
## Wald statistic: 39.662, p-value: 0.00000000030196
##
## Log likelihood: -45.8159 for lag model
## ML residual variance (sigma squared): 0.68907, (sigma: 0.8301)
## Number of observations: 35
## Number of parameters estimated: 7
## AIC: 105.63, (AIC for lm: 125.82)
## LM test for residual autocorrelation
## test value: 1.6155, p-value: 0.20372# Dampak langsung, tidak langsung, dan total pada model SAR.
impacts_sar <- tryCatch(
spatialreg::impacts(model_sar, listw = lw_queen, R = 199),
error = function(e) e
)
if (inherits(impacts_sar, "error")) {
cat("Perhitungan impacts SAR tidak berhasil:", impacts_sar$message, "\n")
} else {
print(summary(impacts_sar, zstats = TRUE))
}## Impact measures (lag, exact):
## Direct Indirect Total
## RLS 0.19586250692 0.35266077714 0.54852328406
## IPM -0.01728372575 -0.03112026008 -0.04840398583
## Miskin 0.01255524181 0.02260637527 0.03516161708
## PDRB 0.00001084641 0.00001952953 0.00003037594
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:198
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 198
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## RLS 0.26299646 0.573048018 0.0407247493 0.0407247493
## IPM -0.03048604 0.188581795 0.0134019246 0.0134019246
## Miskin 0.01363268 0.003865515 0.0002747102 0.0003181798
## PDRB 0.00001075 0.000007827 0.0000005562 0.0000005562
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## RLS -0.879528919 -0.08061172 0.26947212 0.65568110 1.50666279
## IPM -0.377760712 -0.16179809 -0.03822411 0.10879834 0.27894662
## Miskin 0.006802802 0.01122439 0.01354080 0.01588818 0.02112473
## PDRB -0.000003687 0.00000542 0.00001128 0.00001592 0.00002573
##
## ========================================================
## Indirect:
##
## Iterations = 1:198
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 198
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## RLS 0.40953554 1.70036613 0.120839759 0.120839759
## IPM -0.00262626 0.75693642 0.053793129 0.053793129
## Miskin 0.02850611 0.02652551 0.001885086 0.003038327
## PDRB 0.00002024 0.00002751 0.000001955 0.000001561
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## RLS -2.281378842 -0.142707100 0.42132563 1.01727961 3.33806260
## IPM -0.932819455 -0.246635547 -0.05672010 0.15190222 0.89065483
## Miskin 0.007922220 0.016011020 0.02252719 0.03168790 0.08954180
## PDRB -0.000006435 0.000008243 0.00001836 0.00002867 0.00007673
##
## ========================================================
## Total:
##
## Iterations = 1:198
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 198
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## RLS 0.672532 2.19144660 0.155739328 0.15573933
## IPM -0.033112 0.90575710 0.064369354 0.06436935
## Miskin 0.042139 0.02884402 0.002049856 0.00327729
## PDRB 0.000031 0.00003364 0.000002391 0.00000195
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## RLS -3.00472343 -0.23053190 0.74205022 1.72718529 4.40406839
## IPM -1.25610007 -0.41089727 -0.11372097 0.26259284 1.14585505
## Miskin 0.01597635 0.02752009 0.03640328 0.04568676 0.10918401
## PDRB -0.00001034 0.00001394 0.00003055 0.00004435 0.00009896
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## RLS 0.573048017690 1.70036613460 2.19144659881
## IPM 0.188581795431 0.75693642308 0.90575709554
## Miskin 0.003865515327 0.02652550886 0.02884402429
## PDRB 0.000007826774 0.00002751145 0.00003363807
##
## Simulated z-values:
## Direct Indirect Total
## RLS 0.4589431 0.240851388 0.3068895
## IPM -0.1616595 -0.003469591 -0.0365576
## Miskin 3.5267444 1.074667790 1.4609194
## PDRB 1.3740339 0.735819406 0.9215067
##
## Simulated p-values:
## Direct Indirect Total
## RLS 0.6462750 0.80967 0.75893
## IPM 0.8715740 0.99723 0.97084
## Miskin 0.0004207 0.28252 0.14404
## PDRB 0.1694311 0.46184 0.35679Nilai koefisien rho pada model SAR menunjukkan kekuatan
pengaruh spasial dari TPT wilayah tetangga terhadap TPT suatu wilayah.
Jika rho bernilai positif dan signifikan, maka peningkatan
TPT pada wilayah tetangga cenderung diikuti oleh peningkatan TPT pada
wilayah yang diamati. Artinya, masalah pengangguran tidak berdiri
sendiri pada satu kabupaten/kota, tetapi berkaitan dengan kondisi
wilayah sekitarnya. Jika rho tidak signifikan, maka setelah
variabel RLS, IPM, Miskin, dan PDRB dimasukkan ke model, pengaruh TPT
tetangga belum cukup kuat secara statistik.
4.10.3. Perbandingan OLS dan SAR
perbandingan_model <- data.frame(
Model = c("OLS", "SAR"),
AIC = c(AIC(model_ols), AIC(model_sar)),
LogLik = c(as.numeric(logLik(model_ols)), as.numeric(logLik(model_sar)))
) %>%
dplyr::arrange(AIC)
knitr::kable(round_df(perbandingan_model, 4), caption = "Perbandingan Model OLS dan SAR")| Model | AIC | LogLik |
|---|---|---|
| SAR | 105.6318 | -45.8159 |
| OLS | 125.8206 | -56.9103 |
best_model <- perbandingan_model$Model[1]
cat("Model terbaik berdasarkan AIC adalah:", best_model, "\n")## Model terbaik berdasarkan AIC adalah: SARModel dengan nilai AIC paling kecil dipilih sebagai model yang lebih baik karena menunjukkan keseimbangan terbaik antara kecocokan model dan kompleksitas parameter. Jika SAR memiliki AIC lebih kecil daripada OLS, maka penambahan efek spasial melalui lag variabel dependen mampu meningkatkan kesesuaian model. Sebaliknya, jika OLS memiliki AIC lebih kecil, maka model tanpa efek spasial sudah relatif cukup menjelaskan variasi TPT berdasarkan variabel penjelas yang digunakan. Perbandingan ini penting agar pemilihan model tidak hanya berdasarkan teori, tetapi juga berdasarkan ukuran kinerja model.
4.10.4. Prediksi dan Residual Model SAR
jateng$Prediksi_SAR <- as.numeric(fitted(model_sar))
jateng$Residual_SAR <- as.numeric(residuals(model_sar))
summary_pred_res <- data.frame(
Statistik = c("Min", "Mean", "Median", "Max", "SD"),
Prediksi_SAR = c(
min(jateng$Prediksi_SAR, na.rm = TRUE),
mean(jateng$Prediksi_SAR, na.rm = TRUE),
median(jateng$Prediksi_SAR, na.rm = TRUE),
max(jateng$Prediksi_SAR, na.rm = TRUE),
sd(jateng$Prediksi_SAR, na.rm = TRUE)
),
Residual_SAR = c(
min(jateng$Residual_SAR, na.rm = TRUE),
mean(jateng$Residual_SAR, na.rm = TRUE),
median(jateng$Residual_SAR, na.rm = TRUE),
max(jateng$Residual_SAR, na.rm = TRUE),
sd(jateng$Residual_SAR, na.rm = TRUE)
)
)
knitr::kable(round_df(summary_pred_res, 4), caption = "Statistik Prediksi dan Residual SAR")| Statistik | Prediksi_SAR | Residual_SAR |
|---|---|---|
| Min | 3.3649 | -1.4209 |
| Mean | 4.5197 | 0.0000 |
| Median | 4.1136 | -0.1188 |
| Max | 7.6230 | 1.8023 |
| SD | 1.1133 | 0.8422 |
# Cek autokorelasi residual model SAR
moran_residual_sar <- spdep::moran.test(jateng$Residual_SAR, lw_queen, zero.policy = TRUE)
print(moran_residual_sar)##
## Moran I test under randomisation
##
## data: jateng$Residual_SAR
## weights: lw_queen
##
## Moran I statistic standard deviate = -0.4654, p-value = 0.6792
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.08338918 -0.02941176 0.01345142ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = Prediksi_SAR), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_viridis_c(option = "C", name = "Prediksi TPT") +
ggplot2::labs(
title = "Peta Prediksi TPT Model SAR",
subtitle = "Prediksi berdasarkan RLS, IPM, Miskin, PDRB, dan efek spasial"
) +
ggplot2::theme_minimal()
ggplot2::ggplot(jateng) +
ggplot2::geom_sf(ggplot2::aes(fill = Residual_SAR), color = "grey20", linewidth = 0.25) +
ggplot2::geom_sf_text(
data = label_points,
ggplot2::aes(label = LABEL),
size = 2,
check_overlap = TRUE,
color = "black"
) +
ggplot2::scale_fill_gradient2(name = "Residual", midpoint = 0) +
ggplot2::labs(
title = "Peta Residual Model SAR",
subtitle = "Residual positif berarti TPT aktual lebih tinggi daripada prediksi model"
) +
ggplot2::theme_minimal()
Peta prediksi SAR menunjukkan nilai TPT hasil model pada setiap kabupaten/kota dengan mempertimbangkan RLS, IPM, jumlah penduduk miskin, PDRB, dan pengaruh TPT wilayah tetangga. Wilayah dengan prediksi tinggi dapat diinterpretasikan sebagai wilayah yang menurut struktur model memiliki risiko pengangguran relatif lebih tinggi. Peta residual menunjukkan selisih antara nilai aktual dan nilai prediksi. Residual positif berarti TPT aktual lebih tinggi daripada prediksi, sehingga model masih meremehkan nilai TPT di wilayah tersebut. Residual negatif berarti TPT aktual lebih rendah daripada prediksi, sehingga model melebihkan nilai TPT. Jika residual masih membentuk pola spasial, maka masih ada faktor lain yang belum dimasukkan dalam model, misalnya struktur industri, mobilitas tenaga kerja, urbanisasi, atau akses transportasi.
BAB 5. Kesimpulan
Berdasarkan hasil analisis spasial TPT kabupaten/kota di Provinsi Jawa Tengah, diperoleh beberapa kesimpulan sebagai berikut.
- Pola spasial berdasarkan Average Nearest Neighbor Index menunjukkan bahwa lokasi kabupaten/kota memiliki pola Seragam / Teratur dengan nilai ANN Index sebesar 1.5564.
- Matriks pembobot utama yang digunakan adalah queen contiguity karena data berbentuk wilayah administratif kabupaten/kota. Pembobot ini mampu menggambarkan hubungan ketetanggaan berdasarkan batas wilayah.
- Uji autokorelasi spasial global menggunakan Moran’s I, Geary’s C, dan Getis-Ord menunjukkan apakah TPT memiliki keterkaitan spasial secara keseluruhan. Moran’s I sebesar 0.6318 dengan p-value 0 menjadi indikator utama adanya pola pengelompokan spasial.
- Analisis lokal menggunakan LISA, Local Geary’s C, dan Getis-Ord lokal mampu mengidentifikasi wilayah yang menjadi klaster lokal, wilayah berbeda dari tetangganya, serta hot spot atau cold spot TPT.
- Kriging global dan lokal menghasilkan permukaan prediksi TPT di wilayah Jawa Tengah. Kriging global cenderung lebih halus, sedangkan Kriging lokal lebih menekankan variasi setempat.
- Model SAR digunakan untuk menjelaskan TPT berdasarkan RLS, IPM, jumlah penduduk miskin, PDRB, serta efek spasial dari wilayah tetangga. Perbandingan AIC menunjukkan model terbaik antara OLS dan SAR adalah SAR.
Secara umum, pendekatan spasial penting digunakan karena TPT antarwilayah tidak hanya dipengaruhi oleh karakteristik internal wilayah, tetapi juga dapat berkaitan dengan kondisi wilayah sekitar. Hasil analisis ini dapat menjadi dasar dalam perumusan kebijakan ketenagakerjaan berbasis wilayah, terutama untuk mengidentifikasi daerah prioritas penanganan pengangguran.
Daftar Pustaka
Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115. Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R. Springer. Cressie, N. (1993). Statistics for Spatial Data. Wiley. Getis, A., & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189–206. GADM. (2024). Database of Global Administrative Areas, Version 4.1.
Lampiran: Session Info
sessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8 LC_CTYPE=English_Indonesia.utf8
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Indonesia.utf8
##
## time zone: Asia/Jakarta
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] broom_1.0.10 car_3.1-3 carData_3.0-5 knitr_1.51
## [5] sp_2.2-0 gstat_2.1-3 ggplot2_4.0.2 spatialreg_1.3-6
## [9] Matrix_1.5-4.1 spdep_1.3-10 spData_2.3.4 stringr_1.5.2
## [13] tidyr_1.3.1 dplyr_1.1.4 readxl_1.4.5 sf_1.0-20
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 farver_2.1.2 S7_0.2.0
## [5] fastmap_1.2.0 TH.data_1.1-5 digest_0.6.37 lifecycle_1.0.4
## [9] LearnBayes_2.15.2 survival_3.5-5 magrittr_2.0.3 compiler_4.3.1
## [13] rlang_1.1.5 sass_0.4.9 tools_4.3.1 igraph_2.1.4
## [17] yaml_2.3.10 FNN_1.1.4.1 labeling_0.4.3 classInt_0.4-11
## [21] RColorBrewer_1.1-3 multcomp_1.4-29 abind_1.4-8 KernSmooth_2.23-21
## [25] withr_3.0.2 purrr_1.0.4 grid_4.3.1 xts_0.14.1
## [29] e1071_1.7-16 scales_1.4.0 MASS_7.3-60 cli_3.6.4
## [33] mvtnorm_1.3-3 rmarkdown_2.30 intervals_0.15.5 generics_0.1.4
## [37] otel_0.2.0 rstudioapi_0.17.1 DBI_1.2.3 cachem_1.1.0
## [41] proxy_0.4-27 splines_4.3.1 s2_1.1.7 cellranger_1.1.0
## [45] vctrs_0.6.5 boot_1.3-28.1 sandwich_3.1-1 jsonlite_2.0.0
## [49] Formula_1.2-5 jquerylib_0.1.4 units_0.8-7 glue_1.8.0
## [53] codetools_0.2-19 stringi_1.8.7 gtable_0.3.6 deldir_2.0-4
## [57] tibble_3.2.1 pillar_1.11.0 htmltools_0.5.8.1 R6_2.6.1
## [61] wk_0.9.4 evaluate_1.0.5 lattice_0.21-8 backports_1.5.0
## [65] bslib_0.9.0 class_7.3-22 Rcpp_1.0.14 coda_0.19-4.1
## [69] nlme_3.1-162 spacetime_1.3-3 xfun_0.52 zoo_1.8-13
## [73] pkgconfig_2.0.3