BAB 1. Pendahuluan
BAB 2. Tinjauan Pustaka
BAB 3. Metode Penelitian
- 3.1. Data dan Sumber Data
- 3.2. Tahapan Analisis
BAB 4. Hasil dan Pembahasan
BAB 5. Kesimpulan
Daftar Pustaka

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 lokasi pengamatan membentuk pola acak, mengelompok, atau seragam. Salah satu ukuran yang dapat digunakan adalah Average Nearest Neighbor Index (ANNI), yaitu perbandingan antara rata-rata jarak tetangga terdekat yang diamati dan rata-rata jarak harapan pada pola acak. Nilai rasio kurang dari 1 menunjukkan pola mengelompok, nilai mendekati 1 menunjukkan pola acak, dan nilai lebih dari 1 menunjukkan pola seragam atau teratur.

2.2. Matriks Pembobot Spasial

Matriks pembobot spasial merupakan matriks yang menggambarkan hubungan ketetanggaan antarunit wilayah. Pembobot dapat dibentuk berdasarkan kontiguitas, jarak, atau jumlah tetangga terdekat. Pada data kabupaten/kota, pembobot kontiguitas sering digunakan karena wilayah administratif dapat dikatakan bertetangga jika memiliki batas wilayah yang bersinggungan. Dalam penelitian ini, pembobot utama yang digunakan adalah Queen Contiguity, yaitu dua wilayah dianggap bertetangga jika berbagi sisi atau titik sudut.

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 yang berdekatan memiliki nilai yang berbeda.

2.3.1. Moran’s I

Moran’s I merupakan ukuran autokorelasi spasial global yang digunakan untuk mengetahui apakah suatu variabel memiliki pola spasial secara keseluruhan. Nilai Moran’s I positif menunjukkan kecenderungan pengelompokan nilai yang mirip, nilai negatif menunjukkan pola menyebar atau berbeda antarwilayah bertetangga, dan nilai mendekati nol menunjukkan pola acak.

2.3.2. Geary’s C

Geary’s C juga digunakan untuk mengukur autokorelasi spasial global. Nilai Geary’s C kurang dari 1 menunjukkan autokorelasi positif, nilai lebih dari 1 menunjukkan autokorelasi negatif, dan nilai mendekati 1 menunjukkan tidak adanya autokorelasi spasial.

2.3.3. Getis-Ord G

Getis-Ord G digunakan untuk mengidentifikasi konsentrasi nilai tinggi atau rendah secara spasial. Secara lokal, Getis-Ord dapat menunjukkan wilayah Hot Spot jika nilai tinggi terkonsentrasi di sekitar wilayah tersebut, atau Cold Spot jika nilai rendah terkonsentrasi.

2.4. Kriging

Kriging merupakan metode interpolasi geostatistik yang memperkirakan nilai pada lokasi yang tidak teramati berdasarkan nilai pada lokasi pengamatan dan struktur variogram. Ordinary Kriging global menggunakan seluruh titik pengamatan dalam proses prediksi, sedangkan Kriging lokal membatasi jumlah tetangga terdekat yang digunakan dalam prediksi, sehingga lebih menekankan pengaruh lokal.

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. Persamaan umum model SAR adalah:

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

Keterangan:

\(y\) adalah vektor variabel dependen,
\(W\) adalah matriks pembobot spasial,
\(\rho\) adalah koefisien lag spasial,
\(X\) adalah matriks variabel independen,
\(\beta\) adalah vektor parameter regresi,
\(\varepsilon\) adalah error acak.

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

base_dir <- "C:/Users/ACER SWIFT/Documents/Vindy Andrea/Kuliah/Semester 6/STATISTIKA SPASIAL/Tugas Poyek"
data_path <- file.path(base_dir, "Data Spasial Vindy.xlsx")
shp_path  <- file.path(base_dir, "gadm41_IDN_2.shp")

cat("Working directory saat Knit:\n", getwd(), "\n\n")

## Working directory saat Knit:
##  C:/Users/ACER SWIFT/Downloads

cat("Folder data yang dipakai:\n", base_dir, "\n\n")

## Folder data yang dipakai:
##  C:/Users/ACER SWIFT/Documents/Vindy Andrea/Kuliah/Semester 6/STATISTIKA SPASIAL/Tugas Poyek

cat("File Excel dipakai:\n", data_path, "\n")

## File Excel dipakai:
##  C:/Users/ACER SWIFT/Documents/Vindy Andrea/Kuliah/Semester 6/STATISTIKA SPASIAL/Tugas Poyek/Data Spasial Vindy.xlsx

cat("File shapefile dipakai:\n", shp_path, "\n\n")

## File shapefile dipakai:
##  C:/Users/ACER SWIFT/Documents/Vindy Andrea/Kuliah/Semester 6/STATISTIKA SPASIAL/Tugas Poyek/gadm41_IDN_2.shp

if (!file.exists(data_path)) {
  stop(paste0(
    "File Excel tidak ditemukan di path berikut:\n", data_path,
    "\nPastikan nama file dan folder sudah sama persis."
  ))
}

if (!file.exists(shp_path)) {
  stop(paste0(
    "File .shp tidak ditemukan di path berikut:\n", shp_path,
    "\nPastikan nama file shapefile sudah sama persis."
  ))
}

# 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")
}

# 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)
}

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)

# 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)

# Membuat titik label agar seluruh peta menampilkan nama kabupaten/kota.
# LABEL dibuat lebih ringkas dengan menghapus awalan KABUPATEN/KOTA.
jateng <- jateng %>%
  dplyr::mutate(
    LABEL = stringr::str_replace_all(NAMA, "KABUPATEN\\s+|KOTA\\s+", ""),
    LABEL = stringr::str_to_title(LABEL)
  )

label_points <- sf::st_point_on_surface(jateng)

cat("Jumlah kabupaten/kota yang berhasil digabung:", nrow(jateng), "
")

## Jumlah kabupaten/kota yang berhasil digabung: 35

print(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 dicocokkan dengan batas administrasi kabupaten/kota Provinsi Jawa Tengah dari shapefile. Tahap ini penting karena analisis spasial hanya dapat dilakukan jika setiap baris data sosial-ekonomi memiliki pasangan geometri yang tepat. Nama wilayah dibersihkan terlebih dahulu agar perbedaan penulisan, spasi, atau awalan administratif tidak menghambat proses penggabungan. Nilai kosong pada variabel numerik diisi dengan rata-rata masing-masing variabel agar jumlah unit analisis tetap lengkap dan tidak ada kabupaten/kota yang terhapus dari peta. Setelah data digabung, sistem koordinat diubah ke UTM Zona 49S supaya perhitungan jarak, tetangga terdekat, dan kriging menggunakan satuan meter yang lebih sesuai untuk analisis wilayah.

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")

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")

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 memahami karakter awal setiap variabel sebelum masuk ke analisis spasial dan pemodelan. Nilai minimum dan maksimum TPT menunjukkan rentang perbedaan tingkat pengangguran antar kabupaten/kota. Jika rentangnya lebar, berarti kondisi ketenagakerjaan antarwilayah tidak homogen dan ada wilayah yang jauh lebih bermasalah dibanding wilayah lain. Nilai mean menggambarkan tingkat TPT rata-rata Jawa Tengah, sedangkan median membantu melihat nilai tengah yang lebih tahan terhadap pengaruh wilayah ekstrem. Standar deviasi menunjukkan besar kecilnya variasi; standar deviasi yang tinggi menandakan perbedaan TPT antarwilayah cukup kuat. Untuk variabel penjelas, RLS dan IPM merepresentasikan kualitas pendidikan dan pembangunan manusia, Miskin menggambarkan tekanan sosial-ekonomi, sedangkan PDRB menunjukkan kapasitas ekonomi wilayah. Matriks korelasi dibaca sebagai indikasi awal hubungan linear antarvariabel. Korelasi yang kuat antara variabel independen perlu diperhatikan karena dapat mengarah pada multikolinearitas saat model OLS dan SAR dijalankan.

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 = 1.8,
    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 TPT memperlihatkan sebaran tingkat pengangguran terbuka pada setiap kabupaten/kota di Jawa Tengah. Setiap wilayah diberi label nama kabupaten/kota agar posisi wilayah dapat dikenali langsung tanpa harus melihat tabel atribut. Warna yang semakin pekat menunjukkan TPT yang semakin tinggi, sedangkan warna yang lebih muda menunjukkan TPT yang lebih rendah. Jika beberapa wilayah berwarna pekat terlihat saling berdekatan, kondisi tersebut memberi indikasi awal adanya klaster TPT tinggi. Sebaliknya, jika wilayah berwarna muda tampak terkumpul, terdapat indikasi klaster TPT rendah. Peta ini belum membuktikan autokorelasi secara statistik, tetapi sangat membantu membaca pola visual awal. Wilayah yang tampak memiliki TPT tinggi dapat dipandang sebagai area prioritas karena pengangguran tidak hanya muncul sebagai angka tunggal, tetapi juga memiliki konteks lokasi dan kedekatan dengan wilayah sekitar.

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")

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 ANN Index kurang dari 1, jarak antartitik representatif kabupaten/kota lebih dekat dibandingkan pola acak, sehingga posisi wilayah cenderung mengelompok. Jika nilainya lebih dari 1, jarak antartitik relatif lebih jauh dan pola cenderung seragam. Jika p-value tidak signifikan, pola lokasi belum cukup kuat untuk dinyatakan berbeda dari pola acak. Interpretasi ini hanya membaca susunan lokasi kabupaten/kota, bukan nilai TPT-nya. Karena itu, hasil ANN dipakai sebagai informasi awal tentang struktur posisi wilayah, sedangkan hubungan nilai TPT antarwilayah tetap harus diuji melalui Moran’s I, Geary’s C, Getis-Ord, dan LISA.

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 = 1.8,
    color = "black"
  ) +
  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")

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.228571

cat("Wilayah tanpa tetangga Queen:", sum(spdep::card(nb_queen) == 0), "\n")

## Wilayah tanpa tetangga Queen: 0

plot(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 queen dipilih karena data berbentuk wilayah administrasi, bukan titik murni. Dua kabupaten/kota dianggap bertetangga jika berbagi sisi batas ataupun titik sudut. Pendekatan ini lebih inklusif dibanding rook contiguity karena tetap menangkap hubungan wilayah yang bersentuhan pada sudut. Jumlah tetangga pada tabel menunjukkan seberapa banyak wilayah lain yang berhubungan langsung dengan suatu kabupaten/kota. Wilayah dengan banyak tetangga memiliki potensi interaksi spasial yang lebih luas, sedangkan wilayah di tepi provinsi biasanya memiliki tetangga lebih sedikit. Standardisasi baris membuat total bobot setiap wilayah bernilai satu, sehingga lag spasial dapat dibaca sebagai rata-rata kondisi wilayah tetangga. Pada peta struktur tetangga, garis merah memperlihatkan hubungan antarwilayah yang menjadi dasar perhitungan autokorelasi dan model SAR.

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.01314389

print(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: greater

Nilai Moran’s I sebesar 0.6318 dengan p-value sebesar 0. Moran’s I membaca apakah TPT antarwilayah yang bertetangga cenderung mirip secara keseluruhan. Jika nilai Moran’s I positif dan p-value kurang dari 0,05, maka terdapat autokorelasi spasial positif; kabupaten/kota dengan TPT tinggi cenderung berada dekat dengan kabupaten/kota yang juga tinggi, sedangkan wilayah dengan TPT rendah cenderung berada dekat dengan wilayah yang juga rendah. Jika nilai Moran’s I negatif dan signifikan, pola yang muncul adalah ketidaksamaan, misalnya wilayah TPT tinggi dikelilingi wilayah rendah atau sebaliknya. Jika tidak signifikan, maka persebaran TPT belum menunjukkan pola global yang kuat. Hasil Moran global penting karena menjadi dasar apakah analisis lokal dan model spasial perlu diperhatikan lebih lanjut.

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.01773194

Nilai Geary’s C sebesar 0.3476 dengan p-value sebesar 0. Geary’s C menekankan perbedaan nilai TPT antara suatu wilayah dan tetangganya. Nilai kurang dari 1 menunjukkan autokorelasi spasial positif, artinya wilayah bertetangga cenderung memiliki TPT yang tidak terlalu berbeda. Nilai lebih dari 1 menunjukkan autokorelasi spasial negatif, artinya wilayah bertetangga cenderung memiliki perbedaan TPT yang besar. Jika p-value signifikan, pola tersebut dapat dinyatakan kuat secara statistik. Geary’s C melengkapi Moran’s I karena lebih sensitif pada ketidaksamaan lokal, sehingga berguna untuk mendeteksi apakah ada wilayah tertentu yang berbeda tajam dibanding daerah sekitarnya.

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.00003863361

Getis-Ord G global digunakan untuk melihat apakah terdapat konsentrasi nilai TPT tinggi atau rendah secara spasial. Jika p-value signifikan dan arah statistik menunjukkan kecenderungan nilai tinggi, maka TPT tinggi tidak tersebar acak, melainkan terkonsentrasi pada wilayah yang saling berdekatan. Jika kecenderungan yang muncul adalah nilai rendah, maka terdapat konsentrasi wilayah dengan TPT rendah. Perbedaan utama Getis-Ord dengan Moran’s I adalah fokusnya pada keberadaan hot spot atau cold spot. Dengan demikian, hasil Getis-Ord membantu memperjelas apakah masalah pengangguran tinggi terkumpul pada zona tertentu yang perlu menjadi prioritas kebijakan.

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")

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

Ringkasan autokorelasi global memperlihatkan tiga ukuran yang saling melengkapi. Moran’s I digunakan untuk membaca kemiripan nilai TPT antarwilayah bertetangga secara umum. Geary’s C lebih menekankan perbedaan nilai antarwilayah tetangga, sehingga sensitif terhadap wilayah yang berbeda dari lingkungan sekitarnya. Getis-Ord G digunakan untuk melihat apakah nilai tinggi atau rendah terkonsentrasi dalam ruang. Jika ketiga statistik mengarah pada kesimpulan yang sama, maka bukti pola spasial TPT menjadi lebih kuat. Jika hasilnya berbeda, interpretasi perlu dilakukan hati-hati karena setiap statistik menekankan aspek spasial yang tidak sama.

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")

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               29

peta_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 = 1.8,
    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. Peta LISA sudah menampilkan nama kabupaten/kota sehingga wilayah yang masuk kategori tertentu dapat langsung dikenali. Kategori High-High menunjukkan kabupaten/kota dengan TPT tinggi yang dikelilingi tetangga bertingkat TPT tinggi; kategori ini dapat dipahami sebagai klaster pengangguran tinggi. Kategori Low-Low menunjukkan wilayah TPT rendah yang dikelilingi tetangga TPT rendah; kategori ini dapat disebut klaster pengangguran rendah. Kategori High-Low menunjukkan wilayah dengan TPT tinggi tetapi berada di sekitar wilayah TPT rendah, sehingga dapat dipandang sebagai pencilan lokal. Kategori Low-High menunjukkan wilayah TPT rendah di tengah lingkungan tetangga yang relatif tinggi. Wilayah tidak signifikan berarti hubungan lokal antara nilai TPT wilayah tersebut dan tetangganya belum cukup kuat pada taraf signifikansi 5%. Hasil ini lebih operasional dibanding Moran global karena dapat menunjukkan kabupaten/kota mana yang perlu diperhatikan secara spesifik.

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")

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 
##               35

ggplot2::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 = 1.8,
    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. Pada peta, wilayah signifikan menunjukkan kabupaten/kota yang nilai TPT-nya relatif tidak mirip dengan tetangganya. Artinya, wilayah tersebut berpotensi menjadi pencilan lokal atau daerah transisi antara kelompok TPT tinggi dan rendah. Jika suatu kabupaten/kota signifikan pada Local Geary, wilayah tersebut perlu dibaca lebih hati-hati karena masalah penganggurannya mungkin tidak mengikuti pola lingkungan sekitar. Local Geary penting sebagai pelengkap LISA: LISA lebih mudah membaca jenis klaster tinggi-rendah, sedangkan Local Geary lebih menekankan besar kecilnya perbedaan nilai antarwilayah bertetangga.

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")

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               29

ggplot2::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 = 1.8,
    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 secara lebih langsung. Pada peta, Hot Spot menunjukkan kabupaten/kota yang berada dalam konsentrasi TPT tinggi, sedangkan Cold Spot menunjukkan kabupaten/kota yang berada dalam konsentrasi TPT rendah. Jika suatu wilayah menjadi hot spot, maka persoalan pengangguran di wilayah tersebut tidak berdiri sendiri, melainkan berkaitan dengan kondisi wilayah di sekitarnya yang juga tinggi. Wilayah hot spot dapat menjadi prioritas intervensi ketenagakerjaan regional, misalnya perluasan kesempatan kerja lintas kabupaten/kota, penguatan pelatihan kerja, atau dukungan aktivitas ekonomi. Cold spot dapat dijadikan pembanding karena wilayah tersebut menunjukkan kondisi pengangguran yang relatif lebih rendah dan dapat ditelusuri faktor pendukungnya.

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: 1383

4.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 var1

initial_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: Exponential

print(vario_model)

##   model     psill    range
## 1   Nug 0.1798728     0.00
## 2   Exp 0.9790047 24044.99

plot(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")

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 = 1.8,
    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 = 1.8,
    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")

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 = 1.8,
    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 = 1.8,
    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")

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 dalam proses estimasi. Hasilnya cocok untuk membaca kecenderungan umum TPT di seluruh Jawa Tengah. Warna pada peta prediksi menunjukkan area yang diperkirakan memiliki TPT lebih tinggi atau lebih rendah, sedangkan nama kabupaten/kota membantu menghubungkan permukaan prediksi dengan wilayah administrasi. Peta varians global menunjukkan tingkat ketidakpastian prediksi. Area dengan varians lebih tinggi berarti prediksinya kurang pasti, biasanya karena posisi area tersebut relatif jauh dari titik pengamatan atau karena variasi lokal tidak sepenuhnya dijelaskan oleh model variogram.

Kriging lokal menggunakan maksimum 8 tetangga terdekat sehingga lebih peka terhadap perubahan setempat. Jika peta lokal terlihat lebih bervariasi dibanding peta global, hal itu wajar karena pengaruh titik yang jauh dikurangi. Kriging lokal berguna untuk melihat detail spasial yang mungkin tersamarkan pada kriging global. Perbandingan ringkasan kriging global dan lokal menunjukkan apakah metode lokal menghasilkan rentang prediksi dan rata-rata varians yang berbeda. Jika mean variance kriging lokal lebih rendah, prediksi lokal cenderung lebih pasti pada lingkungan terdekat. Jika lebih tinggi, berarti pembatasan tetangga membuat ketidakpastian meningkat pada beberapa area.

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.00001002

Model OLS digunakan sebagai model dasar untuk melihat hubungan TPT dengan RLS, IPM, Miskin, dan PDRB tanpa mempertimbangkan efek spasial. Koefisien pada OLS menunjukkan arah hubungan masing-masing variabel penjelas terhadap TPT ketika variabel lain dianggap tetap. Nilai p-value pada ringkasan OLS digunakan untuk melihat apakah pengaruh variabel tersebut signifikan secara statistik. VIF digunakan untuk memeriksa multikolinearitas; nilai VIF yang tinggi menunjukkan adanya hubungan kuat antarvariabel independen sehingga interpretasi koefisien perlu hati-hati. Uji LM atau Rao’s Score digunakan untuk mengecek apakah residual model OLS masih memiliki indikasi ketergantungan spasial. Jika uji lag spasial signifikan, maka model SAR lebih tepat dipertimbangkan karena TPT suatu wilayah mungkin dipengaruhi oleh 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.35679

Nilai koefisien rho pada model SAR menunjukkan kekuatan pengaruh spasial dari TPT wilayah tetangga terhadap TPT suatu wilayah. Jika rho bernilai positif dan signifikan, maka kenaikan TPT pada wilayah tetangga cenderung diikuti oleh kenaikan TPT pada wilayah yang diamati. Ini berarti masalah pengangguran bersifat saling berkaitan antarwilayah dan tidak cukup dijelaskan hanya oleh kondisi internal kabupaten/kota. Jika rho tidak signifikan, maka setelah RLS, IPM, Miskin, dan PDRB dimasukkan ke model, pengaruh TPT tetangga belum cukup kuat secara statistik. Koefisien variabel independen pada SAR tidak hanya dibaca sebagai pengaruh langsung, karena adanya komponen spasial dapat menghasilkan dampak tidak langsung. Oleh karena itu, bagian impacts SAR digunakan untuk membedakan dampak langsung pada wilayah itu sendiri, dampak tidak langsung ke wilayah lain, dan dampak total.

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")

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: SAR

Model dengan nilai AIC paling kecil dipilih sebagai model yang lebih baik karena AIC mempertimbangkan kecocokan model sekaligus jumlah parameter. Jika SAR memiliki AIC lebih kecil daripada OLS, maka penambahan efek spasial melalui lag variabel dependen mampu meningkatkan kesesuaian model. Hal tersebut menunjukkan bahwa TPT wilayah tetangga membawa informasi penting untuk menjelaskan TPT suatu kabupaten/kota. Jika OLS memiliki AIC lebih kecil, model nonspasial sudah cukup kompetitif dan efek spasial tidak meningkatkan model secara berarti. Perbandingan ini membuat pemilihan model lebih objektif, tidak hanya berdasarkan dugaan teoritis, tetapi juga berdasarkan ukuran performa 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 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.01345142

ggplot2::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 = 1.8,
    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 = 1.8,
    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 estimasi model pada setiap kabupaten/kota dengan mempertimbangkan variabel RLS, IPM, Miskin, PDRB, serta efek spasial dari wilayah tetangga. Wilayah dengan warna prediksi lebih tinggi menunjukkan kabupaten/kota yang menurut model memiliki risiko TPT lebih tinggi. Peta residual memperlihatkan selisih antara TPT aktual dan TPT prediksi. Residual positif berarti TPT aktual lebih tinggi daripada prediksi model, sehingga model masih meremehkan tingkat pengangguran pada wilayah tersebut. Residual negatif berarti TPT aktual lebih rendah daripada prediksi model, sehingga model memperkirakan nilai yang terlalu tinggi. Jika peta residual masih memperlihatkan pola mengelompok, berarti masih ada faktor lain yang belum tertangkap oleh model, misalnya struktur industri, migrasi tenaga kerja, urbanisasi, akses transportasi, atau kebijakan lokal. Uji Moran residual digunakan untuk memeriksa apakah sisa kesalahan model masih memiliki autokorelasi spasial. Model yang baik diharapkan menghasilkan residual yang tidak lagi berpola spasial.

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.

Analisis Spasial Tingkat Pengangguran Terbuka (TPT) per Kabupaten/Kota di Provinsi Jawa Tengah Tahun 2024

Disusun oleh: Vindy Andrea Pratiwi (24050123120014) Dosen Pengampu: Iut Tri Utami, S.Si., M.Si. Program Studi Statistika Fakultas Sains dan Matematika Universitas Diponegoro Semarang 2025

2025-12-20