3 ABSTRAK

Penelitian ini menganalisis pengaruh faktor sosial ekonomi terhadap tingkat kriminalitas antarprovinsi di Indonesia dengan mempertimbangkan keterkaitan spasial antarwilayah. Variabel yang digunakan meliputi Upah Minimum Regional (UMR), PDRB per kapita, dan rata-rata lama sekolah. Hasil regresi OLS menunjukkan bahwa model awal belum mampu menjelaskan adanya ketergantungan spasial, yang kemudian dikonfirmasi melalui uji Lagrange Multiplier. Oleh karena itu, digunakan model regresi spasial, dan berdasarkan perbandingan kriteria AIC, Spatial Durbin Error Model (SDEM) dipilih sebagai model terbaik. Hasil SDEM menunjukkan bahwa UMR dan PDRB berpengaruh signifikan terhadap tingkat kriminalitas, baik secara langsung maupun melalui pengaruh wilayah tetangga, sedangkan pengaruh pendidikan relatif tidak signifikan. Selain itu, interpolasi spasial dengan metode kriging memberikan akurasi terbaik dibandingkan metode lain berdasarkan nilai RMSE terendah, sehingga mampu merepresentasikan pola sebaran kriminalitas secara lebih akurat. Temuan ini menegaskan bahwa analisis spasial penting untuk memahami pola dan determinan kriminalitas di Indonesia.

4 PENDAHULUAN

4.1 Latar Belakang

Tingkat kriminalitas merupakan salah satu indikator penting dalam menggambarkan kondisi keamanan dan stabilitas sosial suatu wilayah. Tingginya jumlah tindak pidana dapat menghambat aktivitas ekonomi, menurunkan kualitas hidup masyarakat, serta mengurangi efektivitas pembangunan daerah. Data Badan Pusat Statistik menunjukkan bahwa jumlah tindak pidana di Indonesia pada tahun 2023 masih menunjukkan variasi yang cukup besar antar provinsi, mencerminkan adanya perbedaan kondisi sosial dan ekonomi antarwilayah [1].

Berbagai faktor sosial ekonomi diduga berpengaruh terhadap variasi tingkat kriminalitas tersebut. Upah Minimum Regional (UMR) mencerminkan tingkat kesejahteraan tenaga kerja dan kemampuan masyarakat dalam memenuhi kebutuhan hidup dasar [2]. Produk Domestik Regional Bruto (PDRB) per kapita menggambarkan kapasitas ekonomi dan tingkat kemakmuran suatu daerah [3]. Selain itu, rata-rata lama sekolah merepresentasikan tingkat pendidikan penduduk yang berperan dalam membentuk kualitas sumber daya manusia dan perilaku sosial masyarakat [4]. Perbedaan kondisi ekonomi dan pendidikan antarwilayah tersebut berpotensi memengaruhi tingkat kriminalitas secara langsung maupun tidak langsung.

Selain dipengaruhi oleh karakteristik internal wilayah, tingkat kriminalitas juga dapat membentuk pola spasial, di mana wilayah dengan tingkat kriminal tinggi cenderung berdekatan dengan wilayah lain yang memiliki karakteristik serupa. Hal ini menunjukkan adanya keterkaitan spasial antarwilayah yang perlu diperhatikan dalam analisis. Oleh karena itu, pendekatan analisis spasial seperti Moran’s I, Local Indicators of Spatial Association (LISA), dan Getis-Ord Gi* digunakan untuk mengidentifikasi pola klaster, hot spot, dan outlier kriminalitas antarprovinsi, sehingga pemahaman terhadap distribusi spasial kriminalitas dapat diperoleh secara lebih komprehensif [5][6].

4.2 Identifikasi Masalah

Berdasarkan latar belakang tersebut, permasalahan yang dikaji dalam penelitian ini adalah bagaimana pola persebaran tingkat kriminal antar provinsi di Indonesia pada tahun 2023 serta apakah terdapat keterkaitan spasial antarwilayah. Selain itu, penelitian ini juga berupaya mengetahui sejauh mana faktor sosial ekonomi, yaitu UMR, PDRB per kapita, dan rata-rata lama sekolah, memengaruhi variasi tingkat kriminalitas antarprovinsi.

4.3 Tujuan Penelitian

Tujuan penelitian ini adalah untuk memetakan pola spasial tingkat kriminal antar provinsi di Indonesia tahun 2023, mengidentifikasi adanya klaster kriminalitas menggunakan uji Moran’s I, LISA, dan Getis-Ord Gi*, serta menganalisis pengaruh UMR, PDRB per kapita, dan rata-rata lama sekolah terhadap tingkat kriminalitas melalui perbandingan model regresi linier dan model spasial.

4.4 Batasan Penelitian

Penelitian ini dibatasi pada wilayah administratif provinsi di Indonesia dengan periode pengamatan tahun 2023. Data yang digunakan merupakan data sekunder yang bersumber dari Badan Pusat Statistik. Variabel yang dianalisis terdiri atas tingkat kriminal sebagai variabel dependen, serta UMR, PDRB per kapita, dan rata-rata lama sekolah sebagai variabel independen. Analisis spasial yang digunakan meliputi uji Moran’s I, LISA, dan Getis-Ord Gi*, tanpa membahas faktor lain seperti penegakan hukum, budaya lokal, atau karakteristik demografis secara mendalam.

5 TINJAUAN PUSTAKA

5.1 Teori Dasar Spatial Dependence

Spatial dependence menjelaskan bahwa nilai suatu variabel pada satu wilayah dapat dipengaruhi oleh nilai variabel yang sama pada wilayah tetangga, sehingga analisis yang mengabaikan hubungan spasial dapat menghasilkan estimasi yang bias. Pembobot spasial membentuk matriks hubungan antarwilayah, biasanya berdasarkan jarak, yang menjadi dasar perhitungan statistik dan model ekonometrika spasial [7].

5.2 Autokorelasi Spasial

Autokorelasi spasial mengukur apakah nilai variabel menunjukkan pola klaster positif, klaster negatif, atau acak. Pengujian dilakukan secara global untuk melihat pola keseluruhan dan secara lokal untuk mendeteksi klaster serta outlier, sehingga kombinasi ukuran global dan lokal memberikan gambaran lengkap pola spasial [8][9][10].

5.2.1 Moran’s I

Moran’s I digunakan untuk menilai ada tidaknya autokorelasi spasial secara global dengan membandingkan nilai suatu wilayah terhadap rata-rata keseluruhan dan kedekatan antarwilayah. Nilai positif menunjukkan klaster wilayah dengan nilai serupa, sedangkan nilai negatif menunjukkan pola berlawanan. Uji ini menentukan apakah pola spasial yang muncul bersifat signifikan atau acak [7][8].

5.2.2 Geary’s C

Geary’s C mengukur autokorelasi spasial dengan fokus pada perbedaan lokal antarwilayah. Nilai mendekati nol menunjukkan kemiripan tinggi, sedangkan nilai mendekati dua menunjukkan perbedaan besar. Uji ini melengkapi hasil Moran’s I untuk memperkuat kesimpulan tentang pola spasial yang ada [8][9].

5.2.3 Getis-Ord Gi

Getis-Ord Gi digunakan untuk mengidentifikasi area dengan konsentrasi nilai tinggi (hot spot) atau rendah (cold spot). Nilai positif menunjukkan wilayah dengan nilai tinggi yang berdekatan, sedangkan nilai negatif menunjukkan wilayah dengan nilai rendah yang saling berdekatan. Analisis ini membantu menemukan wilayah yang menonjol dalam konteks pembangunan [10].

5.3 Model Spatial Econometrics

Model ekonometrika spasial dikembangkan untuk menangani ketergantungan spasial pada data lintas wilayah. Ketika hubungan antarwilayah diabaikan, model konvensional seperti OLS bisa memberikan hasil estimasi yang bias. Oleh karena itu, berbagai model spasial digunakan untuk menggambarkan sumber pengaruh antarwilayah secara lebih akurat [7][8].

5.3.1 Ordinary Least Squares (OLS)

Model OLS digunakan sebagai titik awal dalam analisis spasial untuk menguji hubungan antara variabel dependen dan independen tanpa mempertimbangkan pengaruh spasial. Hasil model ini penting untuk melihat apakah terdapat autokorelasi pada residual yang menjadi dasar untuk memilih model spasial yang lebih sesuai [9].\[Y\ = X\beta\ + \ \varepsilon\ \]

dengan keterangan:

Y = vektor variabel dependen

X= matriks variabel independen

β = vektor koefisien regresi

ε = error

5.3.2 Spatial Lag of X (SLX)

Model SLX menambahkan lag spasial dari variabel independen ke dalam model regresi. Artinya, nilai variabel bebas pada wilayah tetangga ikut memengaruhi nilai variabel dependen pada wilayah tertentu. Model ini berguna untuk melihat seberapa jauh pengaruh eksternal suatu wilayah terhadap wilayah lain tanpa memasukkan efek spasial pada variabel dependen [8].

\(Y\ = \ X\beta\ + \ W\ X\theta\ + \ \varepsilon\)

Keterangan:

W = matriks pembobot spasial

W X = variabel independen yang dipengaruhi oleh tetangga (spatially lagged X)

5.3.3 Spatial Lag Model (SAR)

Model SAR atau spatial lag model memasukkan lag spasial dari variabel dependen, yang berarti nilai Y di suatu wilayah dipengaruhi oleh nilai Y di wilayah tetangga. Koefisien lag (ρ) menunjukkan besarnya penyebaran antarwilayah. Model ini cocok digunakan bila interaksi antarwilayah terjadi melalui fenomena yang menular atau saling memengaruhi secara langsung. [7].

\(Y\ = \ \rho WY\ + \ X\beta\ + \ \varepsilon\)

5.3.4 Spatial Error Model (SEM)

Model SEM digunakan ketika efek spasial muncul melalui komponen error. Dalam model ini, korelasi antarwilayah terjadi karena faktor yang tidak teramati, bukan dari hubungan langsung antarvariabel utama. SEM membantu memperbaiki kesalahan model OLS ketika autokorelasi spasial muncul pada error residual [8].

\(Y\ = \ X\beta\ + \ \mu\)

\(\mu\ = \ \lambda W\mu\ + \ \varepsilon\)

5.3.5 Spatial Durbin Model (SDM)

SDM merupakan pengembangan dari SAR dengan menambahkan lag spasial dari variabel independen. Model ini memungkinkan efek tidak hanya menyebar melalui variabel dependen, tetapi juga melalui faktor-faktor penjelas. SDM sering dianggap sebagai model yang paling fleksibel karena dapat menggambarkan berbagai bentuk keterkaitan spasial antarwilayah [7].

\(Y\ = \ \rho WY\ + \ X\beta\ + W\ X\theta\ + \ \varepsilon\)

5.3.6 Spatial Durbin Error Model (SDEM)

Model SDEM adalah pengembangan dari SDM yang menggabungkan efek lag pada variabel independen serta autokorelasi pada komponen error. Model ini digunakan jika keterkaitan spasial berasal dari gabungan pengaruh langsung antarwilayah dan kesalahan yang saling berkaitan [8].

\(Y\ = \ X\beta\ + W\ X\theta\ + \mu\)

5.4 Interpolasi Spasial

Interpolasi spasial merupakan metode analisis yang digunakan untuk memperkirakan nilai suatu variabel pada lokasi yang tidak teramati berdasarkan nilai pada lokasi terdekat yang telah diketahui. Dalam konteks analisis kriminalitas, interpolasi spasial membantu menggambarkan pola sebaran tingkat kriminal secara kontinu antarwilayah, sehingga memudahkan identifikasi area dengan potensi kerawanan tinggi maupun rendah. Metode interpolasi banyak digunakan pada analisis spasial sosial untuk mendukung visualisasi dan eksplorasi awal sebelum dilakukan pemodelan statistik lebih lanjut [11].

5.4.1 Kriging

Kriging merupakan metode interpolasi geostatistik yang mempertimbangkan struktur ketergantungan spasial antar lokasi melalui fungsi semivariogram. Metode ini tidak hanya memperhitungkan jarak antar titik, tetapi juga pola variabilitas data sehingga menghasilkan estimasi yang bersifat tidak bias dan memiliki varians minimum. Kriging banyak digunakan dalam studi spasial karena mampu menangkap pola spasial yang kompleks secara lebih akurat dibandingkan metode deterministik [12].

Rumus umum Kriging adalah:

\[ \hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i \, Z(s_i) \]

\[ \begin{aligned} \hat{Z}(s_0) &:\ \text{nilai hasil interpolasi (estimasi) pada lokasi } s_0 \\ Z(s_i) &:\ \text{nilai variabel teramati pada lokasi ke-} i \\ \lambda_i &:\ \text{bobot Kriging yang ditentukan berdasarkan model semivariogram} \\ n &:\ \text{jumlah titik pengamatan yang digunakan dalam interpolasi} \end{aligned} \]

5.4.2 Nearest Neighbor

Metode Nearest Neighbor merupakan teknik interpolasi paling sederhana, di mana nilai pada suatu lokasi ditentukan berdasarkan nilai titik pengamatan terdekat. Metode ini tidak melakukan perataan atau pembobotan jarak, sehingga perubahan nilai antarwilayah dapat terlihat secara tajam. Nearest Neighbor sering digunakan sebagai pembanding atau visualisasi awal karena sifatnya yang mudah diterapkan dan tidak memerlukan asumsi statistik [12].

Secara matematis, interpolasi Nearest Neighbor dapat dituliskan sebagai:

\[ \hat{Z}(s_0) = Z(s_k) \]

\[ \begin{aligned} Z(s_k) &:\ \text{nilai variabel teramati pada lokasi } s_k \\ s_k &:\ \text{lokasi pengamatan yang memiliki jarak terdekat terhadap } s_0 \end{aligned} \]

5.4.3 Natural Neighbor

Natural Neighbor merupakan metode interpolasi berbasis geometri yang menggunakan konsep Voronoi tessellation. Nilai interpolasi diperoleh dari rata-rata berbobot titik-titik tetangga alami yang berbagi area pengaruh terhadap lokasi estimasi. Metode ini menghasilkan permukaan interpolasi yang halus tanpa menghasilkan nilai di luar rentang data asli dan bersifat adaptif terhadap distribusi titik data [12].

Rumus umum Natural Neighbor dapat dituliskan sebagai:

\[ \hat{Z}(s_0) = \sum_{i=1}^{n} w_i \, Z(s_i) \]

\[ \begin{aligned} w_i &:\ \text{bobot proporsional yang ditentukan berdasarkan luas irisan poligon Voronoi} \\ &\quad \text{antara lokasi } s_0 \text{ dan titik pengamatan } s_i \\ \sum_{i=1}^{n} w_i &= 1 \end{aligned} \]

5.4.4 Inverse Distance Weighting (IDW)

Inverse Distance Weighting (IDW) adalah metode interpolasi deterministik yang mengasumsikan bahwa pengaruh suatu titik terhadap lokasi estimasi berkurang seiring bertambahnya jarak. Titik yang lebih dekat akan memiliki bobot lebih besar dibandingkan titik yang lebih jauh. IDW banyak digunakan karena sederhana, intuitif, dan efektif untuk menggambarkan pola spasial awal [11].

Rumus IDW adalah:

\[ \hat{Z}(s_0) = \frac{\displaystyle\sum_{i=1}^{n} \frac{Z(s_i)}{d_{i0}^{\,p}}} {\displaystyle\sum_{i=1}^{n} \frac{1}{d_{i0}^{\,p}}} \]

\[ \begin{aligned} d_{i0} &:\ \text{jarak antara lokasi pengamatan } s_i \text{ dan lokasi estimasi } s_0 \\ p &:\ \text{parameter pangkat (biasanya } p = 2\text{)} \\ n &:\ \text{jumlah titik pengamatan yang digunakan} \end{aligned} \]

6 METODOLOGI PENELITIAN

6.1 DATA

Penelitian ini menggunakan data sekunder yang diperoleh dari sumber resmi, yaitu Badan Pusat Statistik (BPS). Seluruh variabel disusun dalam satu tabel data dan diintegrasikan dengan data spasial tingkat provinsi di Indonesia. Data yang digunakan merepresentasikan kondisi sosial ekonomi dan kriminalitas pada tahun 2023, dengan satu variabel dependen dan tiga variabel independen yang relevan dalam menjelaskan variasi tingkat kriminalitas antarwilayah.

Sumber Data:

  • Y : Tingkat Kriminal

  • X1 : Upah Minimum Regional (UMR)

  • X2 : Rata-rata Lama Sekolah

  • X3 : Produk Domestik Regional Bruto (PDRB) per Kapita

Seluruh data numerik diproses dan diselaraskan agar dapat dianalisis secara statistik maupun spasial.

6.1.1 Variabel Dependen

Tingkat Kriminal merupakan variabel utama yang menjadi fokus penelitian. Variabel ini menggambarkan jumlah tindak pidana yang tercatat oleh kepolisian pada masing-masing wilayah. Tingkat kriminal digunakan sebagai indikator kondisi keamanan dan ketertiban masyarakat, serta mencerminkan tingkat kerawanan sosial suatu daerah. Variasi tingkat kriminal antarwilayah diduga tidak hanya dipengaruhi oleh faktor internal wilayah tersebut, tetapi juga oleh kondisi sosial ekonomi wilayah sekitarnya.

6.1.2 Variabel Independen

Penelitian ini menggunakan tiga variabel independen yang secara teoritis dan empiris berkaitan dengan tingkat kriminalitas, yaitu:

  1. Upah Minimum Regional (UMR)
    UMR mencerminkan tingkat kesejahteraan minimum tenaga kerja formal di suatu wilayah. UMR yang lebih tinggi diharapkan dapat meningkatkan kesejahteraan ekonomi masyarakat dan menurunkan dorongan melakukan tindak kriminal, meskipun dalam praktiknya efek ini dapat bervariasi antarwilayah.

  2. Produk Domestik Regional Bruto (PDRB) per Kapita
    PDRB per kapita menggambarkan tingkat kemakmuran dan kapasitas ekonomi suatu wilayah. Wilayah dengan PDRB per kapita yang lebih tinggi umumnya memiliki aktivitas ekonomi yang lebih baik, peluang kerja yang lebih luas, dan kondisi sosial yang relatif lebih stabil, sehingga berpotensi menekan tingkat kriminalitas.

  3. Rata-rata Lama Sekolah
    Rata-rata lama sekolah merepresentasikan tingkat pendidikan penduduk usia 15 tahun ke atas. Pendidikan berperan penting dalam membentuk perilaku sosial, meningkatkan kesempatan kerja, serta menurunkan kecenderungan individu untuk terlibat dalam aktivitas kriminal. Oleh karena itu, variabel ini digunakan untuk menangkap pengaruh kualitas sumber daya manusia terhadap tingkat kriminalitas.

6.2 UNIT SPASIAL

Unit analisis dalam penelitian ini adalah wilayah administratif provinsi di Indonesia. Setiap provinsi direpresentasikan sebagai satu unit spasial berbentuk poligon yang memiliki batas administratif resmi. Pendekatan ini memungkinkan analisis perbandingan tingkat kriminal antarwilayah serta pengujian adanya keterkaitan spasial antarprovinsi yang berdekatan.

Pembentukan Bobot Spasial

Bobot spasial (spatial weights) digunakan untuk mendefinisikan hubungan kedekatan antarwilayah. Bobot ini disusun menggunakan metode Queen Contiguity, yaitu dua wilayah dianggap bertetangga bila berbagi sisi maupun titik sudut peta.
Dalam penelitian ini dibuat dua tipe bobot:

  • lwW: bobot row-standardized, digunakan pada analisis autokorelasi dan sebagian besar model spasial.

  • lwB: bobot binary (0/1), digunakan pada model lain yang tidak memerlukan normalisasi bobot.

6.3 Interpolasi Spasial

Sebelum dilakukan analisis autokorelasi dan pemodelan, penelitian ini menerapkan interpolasi spasial untuk menggambarkan pola kontinu tingkat kriminalitas antarwilayah. Interpolasi spasial bertujuan untuk memperkirakan nilai variabel pada lokasi yang tidak teramati berdasarkan nilai pada lokasi terdekat yang telah diketahui.

Dua metode interpolasi yang digunakan adalah Kriging dan Nearest Neighbor. Kriging digunakan untuk menghasilkan estimasi yang mempertimbangkan struktur ketergantungan spasial antarwilayah, sedangkan Nearest Neighbor digunakan sebagai pendekatan sederhana dengan menetapkan nilai berdasarkan lokasi terdekat. Hasil interpolasi membantu visualisasi awal dan pemahaman pola spasial tingkat kriminalitas sebelum analisis statistik lanjutan dilakukan.

6.4 METODE ANALISIS

Analisis dilakukan melalui beberapa tahap mulai dari eksplorasi data, pengujian autokorelasi spasial, hingga pemodelan ekonometrika spasial.

6.4.1 Analisis Deskriptif

Langkah pertama adalah melihat gambaran umum variabel. Statistik deskriptif seperti nilai rata-rata, minimum, maksimum, dan standar deviasi dihitung menggunakan fungsi pada R.Selain itu, peta tematik dibuat dengan untuk menampilkan distribusi spasial masing-masing variabe

6.4.2 Uji Autokorelasi Spasial Global

6.4.2.1 Moran’s I

Uji Moran’s I digunakan untuk mendeteksi apakah nilai antarwilayah memiliki hubungan spasial. Rumus umum Moran’s I adalah:

\(I\ = \ \frac{n}{S_{0}}\ \cdot \ \frac{\sum_{i}^{}\sum_{j}^{}\ w_{ij}^{}\ (x_{i}\ - \ \underline{x}\ )\ (x_{j}\ - \ \ \underline{x})}{\sum_{i}^{}\ {(x_{i}\ - \ \underline{x}\ )}^{2}}\)

Keterangan :

n = Jumlah Wilayah

\(w_{ij}^{}\) = bobot spasial antara wilayah i dan j

\(x_{i}\) ​= nilai variabel pada wilayah iii

\(\ \underline{x}\) = rata-rata variabel

\(S_{0}\) = \(\sum_{i}^{}\sum_{j}^{}\ w_{ij}^{}\)

Interpretasinya:

  • I > 0 = terdapat klaster wilayah dengan nilai serupa (autokorelasi positif).

  • I < 0 = wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai rendah (autokorelasi negatif).

  • I sama dengan 0 = pola acak tanpa pengaruh spasial.

6.4.2.2 Geary’s C

Uji Geary’s C juga digunakan untuk mendeteksi autokorelasi spasial, namun lebih sensitif terhadap variasi lokal antarwilayah. Rumusnya adalah:

\(C\ = \ \frac{(n\ - \ 1)\ \sum_{i}^{}\sum_{j}^{}\ w_{ij}^{}\ {(x_{i}\ - \ \underline{x}\ )}^{2}}{2\ \ S_{0}\ \ \sum_{i}^{}\ {(x_{i}\ - \ \underline{x}\ )}^{2}}\)

Keterangan :

n = Jumlah Wilayah

\(w_{ij}^{}\) = bobot spasial antara wilayah i dan j

\(x_{i}\) ​= nilai variabel pada wilayah iii

\(\ \underline{x}\) = rata-rata variabel

\(S_{0}\) = \(\sum_{i}^{}\sum_{j}^{}\ w_{ij}^{}\)

Interpretasinya:

  • C < 1 = autokorelasi positif (nilai mirip antarwilayah).

  • C > 1 = autokorelasi negatif (perbedaan besar antarwilayah).

  • C sama dengan 1 = tidak ada pola spasial.

6.4.3 Analisis Autokorelasi Lokal

Untuk melihat klaster atau outlier lokal digunakan metode LISA dan Getis-Ord Gi.

LISA (Local Moran’s I) = mengidentifikasi wilayah dengan pola High-High, Low-Low, High-Low, atau Low-High.

Getis-Ord Gi = mendeteksi area dengan konsentrasi nilai tinggi (hot spot) dan rendah (cold spot).

6.4.4 Pemodelan Regresi Spasial

6.4.4.1 Model OLS

Model OLS digunakan sebagai langkah awal untuk mengetahui hubungan dasar antara IPM dan variabel-variabel independen. Setelah model OLS dijalankan, residualnya diuji menggunakan Moran’s I untuk memastikan apakah terdapat pengaruh spasial yang belum tertangkap model.

6.4.4.2 Model Spasial

Jika ditemukan autokorelasi pada residual OLS, digunakan model-model spasial (SAR, SEM, SLX, SDM, SDEM, SAC, GNC), Perbandingan antar model dilakukan dengan melihat nilai AIC (Akaike Information Criterion) dan Log Likelihood. Model dengan AIC terendah dianggap paling sesuai dengan data.

6.5 ALUR KERJA PENELITIAN

Alur kerja penelitian dimulai dari pengumpulan data tingkat kriminal, UMR, rata-rata lama sekolah, dan tingkat pengangguran terbuka dari BPS, kemudian data tersebut diintegrasikan dengan peta administratif wilayah. Selanjutnya dilakukan interpolasi spasial untuk visualisasi awal, diikuti eksplorasi data melalui statistik deskriptif dan peta tematik. Tahap berikutnya adalah uji autokorelasi spasial global dan lokal untuk mengidentifikasi pola dan klaster kriminalitas. Setelah itu dilakukan pemodelan regresi OLS dan model spasial, evaluasi model menggunakan AIC dan log likelihood, serta penyajian hasil analisis dalam bentuk visualisasi yang informatif.

7 HASIL DAN PEMBAHASAN

7.1 Statistika Deskriptif

Mengimport Library dan pengInputan Data

library(spdep)
## Warning: package 'spdep' was built under R version 4.4.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.4.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.4.3
library(sf)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spData)
library(geodata)
## Warning: package 'geodata' was built under R version 4.4.3
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.4.3
## terra 1.8.60
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
library(skimr)
## Warning: package 'skimr' was built under R version 4.4.3
library(psych)
## Warning: package 'psych' was built under R version 4.4.3
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:terra':
## 
##     describe, distance, rescale
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.4.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(grid)
## 
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
## 
##     depth
library(viridis)
## Warning: package 'viridis' was built under R version 4.4.3
## Loading required package: viridisLite
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.4.3
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library(stringr)
## Warning: package 'stringr' was built under R version 4.4.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.3
library(raster)
## Warning: package 'raster' was built under R version 4.4.3
## 
## Attaching package: 'raster'
## The following object is masked from 'package:skimr':
## 
##     bind
## The following object is masked from 'package:dplyr':
## 
##     select
library (deldir)
## deldir 2.0-4      Nickname: "Idol Comparison"
## 
##      The syntax of deldir() has changed since version 
##      0.0-10.  Read the help!!!.
library (psych)
# SET WORKING DIRECTORY
# ===========================================
setwd("C:\\Users\\Ghozzy\\OneDrive\\Documents\\SMESTER 5\\Spasial UAS")

# Baca shapefile provinsi Indonesia (ADM1)
indo_adm1 <- st_read("IDN_AdminBoundaries_candidate.gdb",
                     layer = "idn_admbnda_adm1_bps_20200401")
## Reading layer `idn_admbnda_adm1_bps_20200401' from data source 
##   `C:\Users\Ghozzy\OneDrive\Documents\SMESTER 5\Spasial UAS\IDN_AdminBoundaries_candidate.gdb' 
##   using driver `OpenFileGDB'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts.  The
## processing may be really slow.  You can skip the processing by setting
## METHOD=SKIP.
## Simple feature collection with 34 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84
# Cek nama kolom provinsi
names(indo_adm1)
##  [1] "admin1Name_en"     "admin1Pcode"       "admin1RefName"    
##  [4] "admin1AltName1_en" "admin1AltName2_en" "admin0Name_en"    
##  [7] "admin0Pcode"       "date"              "validOn"          
## [10] "validTo"           "Shape_Length"      "Shape_Area"       
## [13] "Shape"
data <- read.csv("Data_UAS.csv", sep=";", dec=".")

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

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

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

indo_adm1 <- indo_adm1 %>%
  mutate(provinsi_sf = toupper(str_trim(admin1Name_en)))

# Normalisasi istilah umum (tambah sesuai kebutuhan)
data$Provinsi[data$Provinsi == "METRO JAYA"] <- "DKI JAKARTA"
data$Provinsi[data$Provinsi == "KEP. RIAU"] <- "KEPULAUAN RIAU"
data$Provinsi[data$Provinsi == "KEP. BANGKA BELITUNG"] <- "KEPULAUAN BANGKA BELITUNG"

# Bersihkan angka yang memiliki koma sebagai desimal dan ubah jadi numeric
#    (sesuaikan nama kolom yang ada di CSV)
data <- data %>%
  mutate(
    Y_Tindak = as.numeric(gsub(",", ".", as.character(`Y..Tindak.Kejahatan.`))),
    X1_UMR   = as.numeric(gsub(",", ".", as.character(`X1.UMR.Jam.`))),
    X2_PDRB = as.numeric(gsub(",", ".", as.character(`X2.PDRB.`))),
    X3_RLS = as.numeric(gsub(",", ".", as.character(`X3.RLS.`)))
  )

# Rename kolom supaya lebih nyaman dipakai
data <- dplyr::select(
  data,
  Provinsi,
  Y_Tindak,
  X1_UMR,
  X2_PDRB,
  X3_RLS
)

# Tampilkan provinsi di data tapi TIDAK ADA di shapefile (mismatch dari sisi CSV)
not_in_shp <- data %>% filter(!Provinsi %in% indo_adm1$provinsi_sf)
if(nrow(not_in_shp)>0) {
  message("Provinsi di CSV yang TIDAK ditemukan di shapefile:")
  print(not_in_shp$Provinsi)
} else { message("Semua Provinsi di CSV ada di shapefile (nama cocok).") }
## Semua Provinsi di CSV ada di shapefile (nama cocok).
# Provinsi di shapefile tapi TIDAK ADA di CSV
not_in_data <- indo_adm1 %>% filter(!provinsi_sf %in% data$Provinsi) %>% pull(provinsi_sf)
if(length(not_in_data)>0) {
  message("Provinsi di shapefile yang TIDAK ada di CSV:")
  print(not_in_data)
} else { message("Semua provinsi shapefile ada data CSV.") }
## Semua provinsi shapefile ada data CSV.
# JOIN (left join shapefile dengan data CSV)
indo_merged <- indo_adm1 %>%
  left_join(data, by = c("provinsi_sf" = "Provinsi"))

# Cek baris hasil join yang menghasilkan NA pada nilai (tanda tidak ter-match atau data kosong)
missing_vals <- indo_merged %>%
  dplyr::filter(is.na(Y_Tindak)) %>%
  sf::st_drop_geometry() %>%
  dplyr::select(provinsi_sf)


# Tabel perbandingan nama + nilai (siap dicek)
compare_tbl <- indo_merged %>%
  sf::st_drop_geometry() %>%
  dplyr::select(
    provinsi_sf,
    Y_Tindak,
    X1_UMR,
    X2_PDRB,
    X3_RLS
  )

# Tampilkan 10 baris pertama
print(head(compare_tbl, 10))
##                   provinsi_sf Y_Tindak X1_UMR   X2_PDRB X3_RLS
## 1                        ACEH    12420  17585  41407.59   9.89
## 2                        BALI    11916  18521  62293.79   9.74
## 3                      BANTEN     7392  24685  66147.20   9.48
## 4                    BENGKULU     5579  17410  46300.49   9.35
## 5  DAERAH ISTIMEWA YOGYAKARTA    12061  16478  48359.85  10.16
## 6                 DKI JAKARTA    87426  42354 322619.37  11.42
## 7                   GORONTALO     3574  15138  42340.50   8.48
## 8                       JAMBI     7432  17596  79849.56   9.16
## 9                  JAWA BARAT    45694  21194  52651.59   9.16
## 10                JAWA TENGAH    42304  13381  45167.24   8.44
# =========================================================
# DESKRIPTIF DATA
# =========================================================
data_desc <- indo_merged %>%
  sf::st_drop_geometry() %>%
  dplyr::select(
    Y_Tindak,
    X1_UMR,
    X2_PDRB,
    X3_RLS
  )

# Statistik deskriptif
describe(data_desc)
##          vars  n     mean       sd   median  trimmed      mad      min
## Y_Tindak    1 34 17205.62 20753.99  8247.00 12939.29  6612.40  1701.00
## X1_UMR      2 34 19662.03  5664.10 18094.50 18916.89  2867.35 12933.00
## X2_PDRB     3 34 79973.62 62019.09 63921.37 67696.95 23196.64 23077.98
## X3_RLS      4 34     9.30     0.82     9.29     9.29     0.91     7.34
##                max     range skew kurtosis       se
## Y_Tindak  87426.00  85725.00 1.90     2.75  3559.28
## X1_UMR    42354.00  29421.00 2.01     5.21   971.39
## X2_PDRB  322619.37 299541.39 2.29     5.20 10636.19
## X3_RLS       11.42      4.08 0.10     0.03     0.14

Interpretasi :

Statistik deskriptif menunjukkan bahwa Tingkat Kriminal memiliki variasi yang sangat tinggi antarwilayah, ditandai oleh standar deviasi besar dan median yang jauh di bawah rata-rata. UMR memiliki variasi sedang dengan distribusi relatif seimbang. PDRB menunjukkan ketimpangan ekonomi yang kuat antarwilayah karena sebaran nilainya sangat lebar. Sementara itu, Rata-rata Lama Sekolah relatif homogen dengan variasi paling kecil dibandingkan variabel lainnya.

indo_spatial_clean <- indo_merged %>%
  dplyr::select(
    Y_Tindak,
    X1_UMR,
    X2_PDRB,
    X3_RLS
  ) %>%
  na.omit()
# =========================================================
# Spatial weights (queen) — DATA BERSIH
# =========================================================
nb <- poly2nb(as_Spatial(indo_spatial_clean), queen = TRUE)
## Warning in poly2nb(as_Spatial(indo_spatial_clean), queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(as_Spatial(indo_spatial_clean), queen = TRUE): neighbour object has 11 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
lwW <- nb2listw(nb, style = "W", zero.policy = TRUE)
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)

summary(nb)
## Neighbour list object:
## Number of regions: 34 
## Number of nonzero links: 66 
## Percentage nonzero weights: 5.709343 
## Average number of links: 1.941176 
## 6 regions with no links:
## 2, 17, 20, 21, 22, 23
## 11 disjoint connected subgraphs
## Link number distribution:
## 
## 0 1 2 3 4 5 
## 6 8 8 7 4 1 
## 8 least connected regions:
## 1 5 11 16 18 24 25 31 with 1 link
## 1 most connected region:
## 8 with 5 links
which(card(nb) == 0)
## [1]  2 17 20 21 22 23
# ================================
#  FUNCTION (tidak diubah)
# ================================
analyze_spatial <- function(data_sf, var_col, label) {
  cat("\n\n==============================\n", label, "\n==============================\n")
  
  x_raw <- data_sf[[var_col]]
  
  # ---------- Global tests ----------
  moran_res <- moran.test(x_raw, lwW, randomisation = TRUE, alternative = "two.sided")
  geary_res <- geary.test(x_raw, lwW, randomisation = TRUE, alternative = "two.sided")
  print(moran_res)
  print(geary_res)
  
  # ---------- Local Moran ----------
  x <- scale(x_raw)[, 1]
  lagx <- lag.listw(lwW, x)
  lisa <- localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
  lisa_df <- as.data.frame(lisa)
  names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
  
  alpha <- 0.05
  quad <- dplyr::case_when(
    x >= 0 & lagx >= 0 ~ "High-High",
    x < 0 & lagx < 0 ~ "Low-Low",
    x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
    x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
  )
  
  LISA_sf <- bind_cols(data_sf, lisa_df) |>
    mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
  
  # ---------- Local Geary’s C ----------
  localC_vals <- localC(x_raw, lwW)
  Geary_sf <- mutate(data_sf,
                     LocalC = as.numeric(localC_vals))
  
  # ---------- Getis–Ord Gi* ----------
  sum_x <- sum(x_raw)
  Wb <- listw2mat(lwB)
  Wb_star <- Wb; diag(Wb_star) <- 1
  G_star_raw <- as.numeric(Wb_star %*% x_raw) / sum_x
  Gz <- localG(x_raw, listw = lwW)
  G_sf <- mutate(data_sf,
                 z_Gistar = as.numeric(Gz),
                 hotcold = case_when(
                   z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
                   z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
                   TRUE ~ "Not significant"
                 ))
  
  # ---------- MAPS ----------
  p1 <- ggplot(LISA_sf) +
    geom_sf(aes(fill = quad), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "High-High" = "#d73027",
      "Low-Low" = "#4575b4",
      "High-Low (Outlier)" = "#fdae61",
      "Low-High (Outlier)" = "#74add1",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Local Moran's I (LISA) —", label), fill = "Kategori") +
    theme_minimal()
  
  p2 <- ggplot(Geary_sf) +
    geom_sf(aes(fill = LocalC), color = "white", size = 0.2) +
    scale_fill_viridis_c(option = "magma", direction = -1,
                         name = "Local Geary’s C") +
    labs(title = paste("Local Geary’s C —", label)) +
    theme_minimal()
  
  p3 <- ggplot(G_sf) +
    geom_sf(aes(fill = hotcold), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "Hot spot (p<0.05)" = "#b2182b",
      "Cold spot (p<0.05)" = "#2166ac",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Getis–Ord Gi* —", label), fill = NULL) +
    theme_minimal()
  
  print(p1); print(p2); print(p3)
}

7.2 Peta Deskriptif

7.2.1 Tingkat Kejahatan

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

Interpretasi:

peta menunjukkan bahwa tingkat kejahatan di Indonesia tidak merata antar wilayah. Wilayah dengan jumlah kejahatan relatif tinggi terkonsentrasi di provinsi-provinsi berpenduduk besar dan pusat aktivitas ekonomi, terutama di Pulau Jawa serta beberapa provinsi di Sumatra dan Sulawesi. Sementara itu, wilayah Indonesia bagian timur dan provinsi dengan kepadatan penduduk lebih rendah cenderung menunjukkan tingkat kejahatan yang lebih rendah. Pola ini mengindikasikan bahwa kepadatan penduduk, urbanisasi, dan intensitas aktivitas ekonomi berperan penting dalam perbedaan tingkat kejahatan antar daerah.

analyze_spatial(
  indo_spatial_clean,
  "Y_Tindak",
  "Tindak Kejahatan"
)
## 
## 
## ==============================
##  Tindak Kejahatan 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.4472, p-value = 0.1478
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.20369494       -0.03703704        0.02766961 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 1.6335, p-value = 0.1024
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.71176587        1.00000000        0.03113649

  • Local Moran’s I (LISA)
    Terlihat beberapa wilayah di Pulau Jawa tergolong High–High, yaitu daerah dengan tingkat kejahatan tinggi yang dikelilingi wilayah bertingkat kejahatan tinggi.

  • Local Geary’s C
    Warna gelap menunjukkan perbedaan tingkat kejahatan yang tinggi dengan wilayah tetangga, sedangkan warna terang menandakan kondisi yang lebih homogen.

  • Getis–Ord Gi*
    Beberapa wilayah di Pulau Jawa teridentifikasi sebagai hot spot, yaitu daerah dengan tingkat kejahatan tinggi yang dikelilingi wilayah dengan tingkat kejahatan tinggi.

7.2.2 Upah Minimum Rakyat (Jam)

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

Interpretasi :

Penyebaran UMR di Indonesia tidak merata antar wilayah. UMR relatif tinggi terkonsentrasi di wilayah barat Indonesia, terutama Pulau Jawa dan sebagian Sumatra, yang merupakan pusat kegiatan ekonomi dan industri. Sebaliknya, UMR relatif lebih rendah banyak ditemukan di wilayah Indonesia timur, seperti Maluku dan Papua, yang mencerminkan perbedaan tingkat pembangunan dan aktivitas ekonomi antar daerah.

analyze_spatial(
  indo_spatial_clean,
  "X1_UMR",
  "Upah Minimum Rata-rata per Jam"
)
## 
## 
## ==============================
##  Upah Minimum Rata-rata per Jam 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 2.6298, p-value = 0.008543
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.37130250       -0.03703704        0.02410947 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 3.0071, p-value = 0.002637
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.44043749        1.00000000        0.03462543

  • Local Moran’s I (LISA)
    Terlihat satu klaster High–High yang signifikan (sekitar wilayah Jawa bagian barat–tengah), menunjukkan daerah dengan upah minimum rata-rata per jam tinggi yang dikelilingi oleh daerah dengan nilai tinggi pula.

  • Local Geary’s C
    Wilayah dengan warna lebih gelap menunjukkan perbedaan upah yang tinggi dibandingkan tetangganya (heterogen), sedangkan warna terang menunjukkan kondisi yang relatif homogen.

  • Getis–Ord Gi*
    Terdapat satu wilayah hot spot yang signifikan, yaitu daerah dengan upah minimum rata-rata per jam tinggi di antara wilayah-wilayah sekitarnya.

7.2.3 Produk Domestik Regional Bruto (PDRB) per Kapita

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

Interpretasi :

Peta PDRB menunjukkan pola spasial yang tidak merata di Indonesia. Wilayah dengan PDRB tinggi terkonsentrasi di Pulau Jawa dan sebagian Kalimantan, sedangkan wilayah di Indonesia timur cenderung memiliki PDRB lebih rendah. Hal ini mengindikasikan adanya konsentrasi aktivitas ekonomi di wilayah barat Indonesia dibandingkan wilayah lainnya.

analyze_spatial(
  indo_spatial_clean,
  "X2_PDRB",
  "PDRB"
)
## 
## 
## ==============================
##  PDRB 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.73021, p-value = 0.4653
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07637069       -0.03703704        0.02412089 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 0.99145, p-value = 0.3215
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.81554211        1.00000000        0.03461424

  • Local Moran’s I (LISA)
    Terlihat satu wilayah di Kalimantan bagian utara tergolong High–High, yaitu wilayah dengan nilai PDRB tinggi yang dikelilingi oleh wilayah lain dengan PDRB tinggi. Selain itu, terdapat satu wilayah di Jawa bagian barat–selatan yang termasuk Low–High (outlier).

  • Local Geary’s C
    Wilayah dengan warna lebih gelap menunjukkan perbedaan PDRB yang tinggi dengan wilayah tetangganya (heterogen), sedangkan wilayah berwarna terang menunjukkan PDRB yang relatif homogen dengan daerah sekitar.

  • Getis–Ord Gi*
    Teridentifikasi hot spot PDRB di Kalimantan bagian utara dan satu wilayah di Jawa bagian barat–selatan, yang berarti wilayah tersebut memiliki PDRB tinggi dan dikelilingi oleh wilayah dengan PDRB tinggi. Wilayah lainnya tidak signifikan.

7.2.4 Rata Lama Sekolah (RLS)

ggplot(data = indo_merged) +
  geom_sf(aes(fill = X3_RLS),
          color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(
    option = "cividis",
    direction = -1,
    name = "RLS\n(Tahun)"
  ) +
  labs(
    title = "Peta Rata-rata Lama Sekolah di Indonesia",
    subtitle = "Rata-rata Lama Sekolah per Provinsi",
    caption = "Sumber: Data_UAS.csv"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

Interpretasi:

Peta menunjukkan rata-rata lama sekolah di Indonesia belum merata. Wilayah dengan warna lebih terang (kuning–oranye) menandakan lama sekolah lebih tinggi, terutama di sebagian Jawa dan beberapa provinsi di Indonesia timur. Sementara itu, wilayah berwarna lebih gelap (ungu) menunjukkan lama sekolah lebih rendah, yang masih banyak ditemukan di kawasan timur Indonesia dan beberapa wilayah luar Jawa.

analyze_spatial(
  indo_spatial_clean,
  "X3_RLS",
  "Rata Lama Sekolah"
)
## 
## 
## ==============================
##  Rata Lama Sekolah 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = -0.67845, p-value = 0.4975
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.15766245       -0.03703704        0.03161173 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 0.50419, p-value = 0.6141
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##         0.9167350         1.0000000         0.0272732

  • Local Moran’s I
    Terlihat beberapa wilayah di Papua tergolong High–Low (outlier), yaitu wilayah dengan nilai rata-rata lama sekolah relatif tinggi namun dikelilingi oleh wilayah lain yang memiliki nilai rata-rata lama sekolah lebih rendah, sementara sebagian besar wilayah lainnya tidak signifikan secara spasial.

  • Local Geary’s C
    Wilayah dengan warna gelap (ungu hingga hitam), terutama di Papua, menunjukkan perbedaan rata-rata lama sekolah yang tinggi dibandingkan wilayah sekitarnya (heterogen), sedangkan wilayah berwarna kuning menunjukkan kondisi yang lebih homogen dengan tetangganya.

  • Getis–Ord Gi*
    Terlihat beberapa wilayah di Papua berwarna biru sebagai cold spot (p < 0,05), yaitu wilayah dengan nilai rata-rata lama sekolah rendah yang dikelilingi oleh wilayah lain dengan nilai yang juga relatif rendah.

7.3 HASIL UJI AUTOKORELASI

7.3.1 UJi Moran I

Variabel Moran I P-Value Interpretasi
Tingkat Kriminal 0.2037 0.1478 Tidak terdapat autokorelasi spasial yang signifikan (p > 0,05), pola sebaran cenderung acak
Upah Minimum Rata-rata per Jam 0.3713 0.0085 Terdapat autokorelasi spasial positif dan signifikan, menunjukkan klaster wilayah dengan nilai UMR yang serupa
PDRB 0.0764 0.4653 Tidak terdapat autokorelasi spasial yang signifikan, distribusi PDRB relatif acak
Rata-rata Lama Sekolah -0.1577 0.4975 Tidak terdapat autokorelasi spasial yang signifikan, tidak terbentuk pola klaster

7.3.2 Uji Geary C

Variabel Geary C P-Value Interpretasi
Tingkat Kriminal 0.7118 0.1024 Tidak terdapat autokorelasi spasial yang signifikan meskipun nilai menunjukkan kecenderungan positif
Upah Minimum Rata-rata per Jam 0.4404 0.0026 Terdapat autokorelasi spasial positif dan signifikan, menunjukkan kemiripan nilai antarwilayah bertetangga
PDRB 0.8155 0.3215 Tidak terdapat autokorelasi spasial yang signifikan
Rata-rata Lama Sekolah 0.9167 0.6141 Tidak terdapat autokorelasi spasial yang signifikan, variasi antarwilayah relatif kecil

7.4 ESTIMASI MODEL

indo_merged_clean <- indo_merged %>%
  dplyr::select(
    Y_Tindak,
    X1_UMR,
    X2_PDRB,
    X3_RLS
  ) %>%
  na.omit()
indo_merged_scaled <- indo_merged_clean %>%
  mutate(
    Y_Tindak_scaled = as.numeric(scale(Y_Tindak)),
    X1_UMR_scaled   = as.numeric(scale(X1_UMR)),
    X2_PDRB_scaled  = as.numeric(scale(X2_PDRB)),
    X3_RLS_scaled  = as.numeric(scale(X3_RLS))
  )
# =========================================================
# Model Regresi Spasial – OLS (Indonesia)
# =========================================================
model_ols <- lm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled
)

summary(model_ols)
## 
## Call:
## lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.36149 -0.46505 -0.35813  0.07895  2.48261 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -7.420e-17  1.688e-01   0.000    1.000
## X1_UMR_scaled   6.180e-02  2.386e-01   0.259    0.797
## X2_PDRB_scaled  2.973e-01  2.554e-01   1.164    0.254
## X3_RLS_scaled   4.466e-03  2.107e-01   0.021    0.983
## 
## Residual standard error: 0.9843 on 30 degrees of freedom
## Multiple R-squared:  0.1193, Adjusted R-squared:  0.03121 
## F-statistic: 1.354 on 3 and 30 DF,  p-value: 0.2756
# =========================================================
# Lagrange Multiplier Test
# =========================================================
lm_tests <- lm.LMtests(
  model_ols,
  listw = lwW,
  test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA")
)
## Please update scripts to use lm.RStests in place of lm.LMtests
print(lm_tests)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled +
## X3_RLS_scaled, data = indo_merged_scaled)
## test weights: listw
## 
## RSerr = 4.7846, df = 1, p-value = 0.02871
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled +
## X3_RLS_scaled, data = indo_merged_scaled)
## test weights: listw
## 
## RSlag = 2.8671, df = 1, p-value = 0.09041
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled +
## X3_RLS_scaled, data = indo_merged_scaled)
## test weights: listw
## 
## adjRSerr = 5.2093, df = 1, p-value = 0.02247
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled +
## X3_RLS_scaled, data = indo_merged_scaled)
## test weights: listw
## 
## adjRSlag = 3.2917, df = 1, p-value = 0.06963
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled +
## X3_RLS_scaled, data = indo_merged_scaled)
## test weights: listw
## 
## SARMA = 8.0763, df = 2, p-value = 0.01763
#SAR
slag_model <- lagsarlm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  method = "Matrix",
  zero.policy = TRUE
)
summary(slag_model)
## 
## Call:lagsarlm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled, listw = lwB, method = "Matrix", 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.71540 -0.74410 -0.31316  0.16985  2.90606 
## 
## Type: lag 
## Regions with no neighbours included:
##  2 17 20 21 22 23 
## Coefficients: (asymptotic standard errors) 
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.060132   0.198659  0.3027   0.7621
## X1_UMR_scaled  0.160883   0.280936  0.5727   0.5669
## X2_PDRB_scaled 0.121363   0.300642  0.4037   0.6864
## X3_RLS_scaled  0.077217   0.247976  0.3114   0.7555
## 
## Rho: -0.29862, LR test value: 9.5451, p-value: 0.0020048
## Asymptotic standard error: 0.065086
##     z-value: -4.5881, p-value: 4.4722e-06
## Wald statistic: 21.051, p-value: 4.4722e-06
## 
## Log likelihood: -40.8046 for lag model
## ML residual variance (sigma squared): 1.341, (sigma: 1.158)
## Number of observations: 34 
## Number of parameters estimated: 6 
## AIC: 93.609, (AIC for lm: 101.15)
## LM test for residual autocorrelation
## test value: 90.314, p-value: < 2.22e-16
#SEM
serr_model <- errorsarlm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  method = "Matrix",
  zero.policy = TRUE
)
summary(serr_model)
## 
## Call:errorsarlm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled, listw = lwB, method = "Matrix", 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.52391 -0.81160 -0.30647  0.17383  3.68175 
## 
## Type: error 
## Regions with no neighbours included:
##  2 17 20 21 22 23 
## Coefficients: (asymptotic standard errors) 
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)     0.10296    0.11473  0.8974   0.3695
## X1_UMR_scaled   0.29448    0.21887  1.3454   0.1785
## X2_PDRB_scaled -0.29899    0.29224 -1.0231   0.3063
## X3_RLS_scaled   0.38984    0.25651  1.5198   0.1286
## 
## Lambda: -0.5, LR test value: 11.399, p-value: 0.0007347
## Asymptotic standard error: 0.012525
##     z-value: -39.919, p-value: < 2.22e-16
## Wald statistic: 1593.5, p-value: < 2.22e-16
## 
## Log likelihood: -39.87749 for error model
## ML residual variance (sigma squared): 1.7148, (sigma: 1.3095)
## Number of observations: 34 
## Number of parameters estimated: 6 
## AIC: 91.755, (AIC for lm: 101.15)
# SDM
sdm_model <- lagsarlm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  Durbin = TRUE,
  method = "Matrix",
  zero.policy = TRUE
)
summary(sdm_model)
## 
## Call:lagsarlm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled, listw = lwB, Durbin = TRUE, 
##     method = "Matrix", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.75547 -0.68796 -0.16844  0.63757  3.61064 
## 
## Type: mixed 
## Regions with no neighbours included:
##  2 17 20 21 22 23 
## Coefficients: (asymptotic standard errors) 
##                    Estimate Std. Error z value  Pr(>|z|)
## (Intercept)         0.11270    0.33804  0.3334 0.7388454
## X1_UMR_scaled      -0.54532    0.35388 -1.5410 0.1233202
## X2_PDRB_scaled      0.34034    0.31299  1.0874 0.2768688
## X3_RLS_scaled       0.21746    0.23076  0.9424 0.3460101
## lag.(Intercept)     0.22336    0.15322  1.4578 0.1448924
## lag.X1_UMR_scaled   0.99785    0.27968  3.5679 0.0003599
## lag.X2_PDRB_scaled -1.21344    0.30525 -3.9753  7.03e-05
## lag.X3_RLS_scaled   0.69411    0.21596  3.2141 0.0013084
## 
## Rho: -0.48166, LR test value: 16.635, p-value: 4.5296e-05
## Asymptotic standard error: 0.0084426
##     z-value: -57.051, p-value: < 2.22e-16
## Wald statistic: 3254.8, p-value: < 2.22e-16
## 
## Log likelihood: -33.16592 for mixed model
## ML residual variance (sigma squared): 1.1278, (sigma: 1.062)
## Number of observations: 34 
## Number of parameters estimated: 10 
## AIC: 86.332, (AIC for lm: 100.97)
## LM test for residual autocorrelation
## test value: 6.4521, p-value: 0.011082
# SLX
slx_model <- lmSLX(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  zero.policy = TRUE
)
summary(slx_model)
## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                     Estimate  Std. Error  t value   Pr(>|t|)
## (Intercept)         -0.08400   0.29835    -0.28156   0.78051
## X1_UMR_scaled       -0.14684   0.31223    -0.47031   0.64206
## X2_PDRB_scaled       0.36084   0.27580     1.30834   0.20221
## X3_RLS_scaled        0.02050   0.20364     0.10065   0.92060
## lag..Intercept.      0.15969   0.13518     1.18128   0.24818
## lag.X1_UMR_scaled    0.39816   0.24665     1.61429   0.11854
## lag.X2_PDRB_scaled  -0.67107   0.26859    -2.49849   0.01912
## lag.X3_RLS_scaled    0.31434   0.19059     1.64932   0.11111
# SDEM
sdem_model <- errorsarlm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  Durbin = TRUE,
  method = "Matrix",
  zero.policy = TRUE
)
summary(sdem_model)
## 
## Call:errorsarlm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled, listw = lwB, Durbin = TRUE, 
##     method = "Matrix", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.26998 -0.57894 -0.19136  0.45413  2.48493 
## 
## Type: error 
## Regions with no neighbours included:
##  2 17 20 21 22 23 
## Coefficients: (asymptotic standard errors) 
##                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)        -0.064979   0.260887 -0.2491   0.80331
## X1_UMR_scaled      -1.086681   0.259775 -4.1832 2.875e-05
## X2_PDRB_scaled      0.638103   0.257217  2.4808   0.01311
## X3_RLS_scaled       0.184108   0.213857  0.8609   0.38930
## lag.(Intercept)     0.268582   0.112755  2.3820   0.01722
## lag.X1_UMR_scaled   1.170248   0.187281  6.2486 4.141e-10
## lag.X2_PDRB_scaled -1.130134   0.187848 -6.0162 1.785e-09
## lag.X3_RLS_scaled   0.285559   0.137191  2.0815   0.03739
## 
## Lambda: -0.5, LR test value: 36.609, p-value: 1.4437e-09
## Asymptotic standard error: 0.012525
##     z-value: -39.919, p-value: < 2.22e-16
## Wald statistic: 1593.5, p-value: < 2.22e-16
## 
## Log likelihood: -23.17923 for error model
## ML residual variance (sigma squared): 0.64214, (sigma: 0.80134)
## Number of observations: 34 
## Number of parameters estimated: 10 
## AIC: 66.358, (AIC for lm: 100.97)
Model Hasil Uji Utama Signifikansi Interpretasi
OLS (Ordinary Least Squares) R² = 0.1193; Adj R² = 0.0312; F p-value = 0.2756 Tidak signifikan Tidak terdapat hubungan linear yang kuat antara UMR, PDRB, dan RLS terhadap tingkat kriminal. Model OLS belum mampu menjelaskan variasi data secara memadai.
LM Test (Spatial Dependence) RSerr p = 0.0287; adjRSerr p = 0.0225; SARMA p = 0.0176 Signifikan Terdapat indikasi dependensi spasial pada residual OLS, sehingga model OLS tidak memadai dan perlu digunakan model regresi spasial.
SAR (Spatial Autoregressive Model) ρ = −0.2986 (p < 0.01); AIC = 93.609 Signifikan (parameter spasial) Terdapat pengaruh spasial pada variabel dependen, namun variabel independen tidak signifikan secara langsung. Hubungan kriminal antarwilayah bersifat saling memengaruhi.
SEM (Spatial Error Model) λ = −0.5 (p < 0.001); AIC = 91.755 Signifikan (error spasial) Ketergantungan spasial terdapat pada komponen error, menunjukkan adanya faktor spasial tak teramati yang memengaruhi tingkat kriminal.
SDM (Spatial Durbin Model) ρ = −0.4817 (p < 0.001); AIC = 86.332 Signifikan Efek spasial baik langsung maupun tidak langsung berperan penting. Variabel di wilayah tetangga (spatial spillover) berpengaruh signifikan terhadap tingkat kriminal.
SLX (Spatial Lag of X) Lag PDRB signifikan (p = 0.0191) Parsial signifikan Pengaruh tidak langsung variabel ekonomi wilayah sekitar, khususnya PDRB, berperan terhadap tingkat kriminal suatu wilayah.
SDEM (Spatial Durbin Error Model) λ signifikan (p < 0.001); AIC = 66.358 Sangat signifikan Model terbaik. Menggabungkan efek error spasial dan spillover variabel independen. Menunjukkan pengaruh kuat faktor ekonomi lokal dan wilayah tetangga terhadap kriminalitas.

7.5 PERBANDINGAN MODEL

model_list <- list(
  OLS  = model_ols,
  SAR  = slag_model,
  SEM  = serr_model,
  SDM  = sdm_model,
  SLX  = slx_model,
  SDEM = sdem_model
)

model_comparison <- data.frame(
  Model  = names(model_list),
  AIC    = sapply(model_list, AIC),
  LogLik = sapply(model_list, logLik)
)

print(model_comparison)
##      Model       AIC    LogLik
## OLS    OLS 101.15431 -45.57716
## SAR    SAR  93.60919 -40.80460
## SEM    SEM  91.75497 -39.87749
## SDM    SDM  86.33185 -33.16592
## SLX    SLX 100.96732 -41.48366
## SDEM  SDEM  66.35847 -23.17923
Model AIC Log Likelihood Kinerja Model
OLS 101.15 −45.58 Terburuk, tidak mempertimbangkan dependensi spasial
SAR 93.61 −40.80 Lebih baik dari OLS, menangkap efek lag dependen
SEM 91.75 −39.88 Menangkap dependensi spasial pada error
SDM 86.33 −33.17 Baik, menangkap spillover variabel independen
SLX 100.97 −41.48 Spillover terbatas, kinerja relatif lemah
SDEM 66.36 −23.18 Terbaik (AIC terendah, LogLik tertinggi)

Interpretasi :

Hasil perbandingan menunjukkan bahwa model SDEM merupakan model paling optimal dalam menjelaskan variasi tingkat kriminal antarwilayah. Nilai AIC yang paling rendah (66.36) dan Log Likelihood tertinggi (−23.18) mengindikasikan bahwa SDEM memberikan keseimbangan terbaik antara ketepatan model dan kompleksitas parameter.

Secara substantif, model SDEM menunjukkan bahwa tingkat kriminal di suatu wilayah tidak hanya dipengaruhi oleh kondisi sosial-ekonomi wilayah itu sendiri, seperti UMR, PDRB, dan Rata Lama Sekolah, tetapi juga oleh pengaruh tidak langsung (spillover) dari wilayah tetangga. Selain itu, signifikansi komponen error spasial menandakan adanya faktor spasial lain yang tidak teramati seperti karakteristik sosial, mobilitas penduduk, atau efektivitas penegakan hukum yang menyebar lintas wilayah dan turut memengaruhi kriminalitas.

Dengan demikian, hasil ini menegaskan bahwa pendekatan spasial sangat krusial dalam analisis kriminalitas, dan kebijakan pengendalian kejahatan tidak dapat dirancang secara parsial per wilayah, melainkan perlu mempertimbangkan keterkaitan dan dinamika antarwilayah secara regional.

7.6 INTERPOLASI SPASIAL

library(sf)
library(sp)
library(gstat)
## Warning: package 'gstat' was built under R version 4.4.3
library(spdep)

# Centroid
indo_centroid <- st_centroid(st_transform(indo_merged_scaled, 4326))
## Warning: st_centroid assumes attributes are constant over geometries
sp_pts <- as(indo_centroid, "Spatial")

# Model regresi
reg_model <- lm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = sp_pts@data
)

# Residual
sp_pts@data$resid <- residuals(reg_model)
moran_resid <- moran.test(sp_pts@data$resid, lwB, zero.policy = TRUE)
moran_resid
## 
##  Moran I test under randomisation
## 
## data:  sp_pts@data$resid  
## weights: lwB  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.9402, p-value = 0.02618
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.27139647       -0.03703704        0.02527065
vgm_emp <- variogram(resid ~ 1, sp_pts)

vgm_fit <- fit.variogram(
  vgm_emp,
  vgm(model = "Sph")
)
## Warning in fit.variogram(vgm_emp, vgm(model = "Sph")): No convergence after 200
## iterations: try different initial values?
plot(vgm_emp, vgm_fit, main = "Variogram Residual Regression")

7.6.1 Kriging

library(sf)
library(gstat)
library(dplyr)
library(ggplot2)
library(viridis)

# =========================================================
# 1) Data titik: centroid provinsi + Y
# =========================================================
indo_proj <- st_transform(indo_merged_scaled, 3857)
pts_proj  <- st_transform(st_centroid(indo_merged_scaled), 3857)
## Warning: st_centroid assumes attributes are constant over geometries
# =========================================================
# 2) Variogram
# =========================================================
vgm_emp <- variogram(Y_Tindak_scaled ~ 1, pts_proj)

vgm_fit <- fit.variogram(
  vgm_emp,
  vgm(psill = 1, model = "Exp", range = 300000, nugget = 0.1)
)
## Warning in fit.variogram(vgm_emp, vgm(psill = 1, model = "Exp", range = 3e+05,
## : No convergence after 200 iterations: try different initial values?
# =========================================================
# 3) Ordinary Kriging langsung ke POLYGON provinsi
# =========================================================
ok <- krige(
  Y_Tindak_scaled ~ 1,
  locations = pts_proj,
  newdata   = indo_proj,
  model     = vgm_fit
)
## [using ordinary kriging]
# =========================================================
# 4) Visualisasi Prediksi per Provinsi
# =========================================================
ggplot(st_as_sf(ok)) +
  geom_sf(aes(fill = var1.pred), color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(title = "Ordinary Kriging — Prediksi Y_Tindak_scaled")

# =========================================================
# 5) Cross-Validation + RMSE
# =========================================================
cv_ok <- krige.cv(
  Y_Tindak_scaled ~ 1,
  locations = pts_proj,
  model = vgm_fit,
  nfold = nrow(pts_proj)
)

rmse_ok <- sqrt(mean(cv_ok$residual^2, na.rm = TRUE))
rmse_ok
## [1] 0.9936867

Hasil Ordinary Kriging untuk prediksi tingkat kriminal terstandarisasi (Y_Tindak_scaled) menunjukkan adanya pola spasial yang jelas antarwilayah di Indonesia. Nilai prediksi yang lebih tinggi (warna kuning–oranye) terkonsentrasi terutama di wilayah Jawa bagian tengah–timur dan sebagian Sumatra, yang mengindikasikan area dengan tingkat kriminal relatif lebih tinggi dibandingkan rata-rata nasional. Sebaliknya, wilayah Kalimantan, Sulawesi, Maluku, dan Papua didominasi nilai prediksi lebih rendah (warna ungu–biru), mencerminkan tingkat kriminal relatif lebih rendah atau mendekati rata-rata setelah distandarisasi.

Nilai RMSE sebesar 0,993 menunjukkan bahwa kesalahan prediksi model kriging berada sekitar satu satuan deviasi standar dari nilai aktual, sehingga akurasi prediksi tergolong sedang. Artinya, kriging cukup mampu menangkap pola spasial umum dan gradien regional tingkat kriminal, namun belum optimal untuk prediksi presisi di tingkat wilayah kecil. Oleh karena itu, hasil kriging ini lebih tepat digunakan untuk analisis eksploratif dan identifikasi pola spasial makro, bukan sebagai alat prediksi tunggal untuk penentuan kebijakan yang bersifat operasional tanpa dukungan analisis tambahan.

7.6.2 Inverse Distance Weighting (IDW)

library(sf)
library(sp)
library(gstat)
library(raster)
library(ggplot2)
library(viridis)

# =========================================================
# 1) Data provinsi & titik (centroid)
# =========================================================
# pastikan sf
indo_sf <- st_transform(indo_merged_scaled, 4326)

# centroid provinsi
pts_sf <- st_centroid(indo_sf)
## Warning: st_centroid assumes attributes are constant over geometries
# ubah ke SpatialPointsDataFrame
pts_sp <- as(pts_sf, "Spatial")

# sanity check
stopifnot("Y_Tindak_scaled" %in% names(pts_sp@data))

# =========================================================
# 2) Grid raster (long-lat)
# =========================================================
r <- raster(
  extent(pts_sp),
  res = 0.25,  # ~25 km
  crs = CRS("+proj=longlat +datum=WGS84")
)

# =========================================================
# 3) IDW interpolation
# =========================================================
idw_r <- interpolate(
  r,
  gstat(
    formula   = Y_Tindak_scaled ~ 1,
    locations = pts_sp,
    set = list(idp = 2)  # power IDW
  )
)
## [inverse distance weighted interpolation]
# mask ke batas Indonesia
idw_r <- mask(idw_r, as(indo_sf, "Spatial"))

# =========================================================
# 4) Data frame untuk ggplot
# =========================================================
idw_df <- as.data.frame(idw_r, xy = TRUE, na.rm = TRUE)
colnames(idw_df) <- c("x", "y", "Y_idw")

# =========================================================
# 5) Visualisasi (BERWARNA seperti choropleth)
# =========================================================
ggplot() +
  geom_tile(
    data = idw_df,
    aes(x = x, y = y, fill = Y_idw)
  ) +
  geom_sf(
    data = indo_sf,
    fill = NA,
    color = "black",
    linewidth = 0.4
  ) +
  scale_fill_viridis_c(
    option = "inferno",
    name = "Prediksi\n(Y scaled)"
  ) +
  coord_sf(expand = FALSE) +
  theme_minimal(base_size = 14) +
  labs(
    title = "Inverse Distance Weighting (IDW)",
    subtitle = "Interpolasi Tingkat Kriminalitas Antar Provinsi"
  )

coords <- st_coordinates(pts_sf)
z <- pts_sf$Y_Tindak_scaled
n <- length(z)
zhat <- numeric(n)
p <- 2

for (i in seq_len(n)) {
  d <- sqrt((coords[-i,1]-coords[i,1])^2 + (coords[-i,2]-coords[i,2])^2)
  w <- 1/(d^p)
  zhat[i] <- sum(z[-i]*w)/sum(w)
}

rmse_idw <- sqrt(mean((z - zhat)^2, na.rm = TRUE))
rmse_idw
## [1] 1.095725

7.6.3 Natural Neighbor Interpolation

# Titik observasi (sf, CRS meter)
pts_sf <- st_as_sf(sp_pts) |> st_transform(3857)

# Buat grid prediksi
res <- 50000  # 50 km
bb  <- st_bbox(pts_sf)
grid_df <- expand.grid(
  x = seq(bb["xmin"], bb["xmax"], by = res),
  y = seq(bb["ymin"], bb["ymax"], by = res)
)
grid_sf <- st_as_sf(grid_df, coords = c("x","y"), crs = st_crs(pts_sf))

# Koordinat
Xobs  <- st_coordinates(pts_sf)
Xgrid <- st_coordinates(grid_sf)

# Nearest Neighbor (1-NN)
library(FNN)
## Warning: package 'FNN' was built under R version 4.4.3
nn_ix <- get.knnx(data = Xobs, query = Xgrid, k = 1)$nn.index[,1]
grid_sf$Y_nn <- pts_sf$Y_Tindak_scaled[nn_ix]

# Visualisasi
ggplot() +
  geom_sf(data = grid_sf, aes(color = Y_nn), size = 1) +
  geom_sf(data = st_transform(indo_merged_scaled, 3857),
          fill = NA, color = "black", linewidth = 0.3) +
  scale_color_viridis_c(name = "Y (NN)") +
  theme_minimal() +
  labs(title = "Natural Neighbor (Voronoi / 1-NN) — Y_Tindak_scaled")

Interpretasi :

Peta ini menunjukkan hasil interpolasi Natural Neighbor, di mana nilai Y pada setiap lokasi grid ditentukan oleh provinsi terdekat secara geografis (berdasarkan pembagian Voronoi). Pola warna yang tajam dan tersegmentasi mencerminkan bahwa perubahan nilai terjadi secara diskrit mengikuti batas pengaruh provinsi terdekat, tanpa proses penghalusan. Metode ini cocok untuk verifikasi lokal dan eksplorasi pola spasial kasar, tetapi kurang merepresentasikan transisi gradual antarwilayah.

coords <- Xobs
z <- pts_sf$Y_Tindak_scaled
n <- length(z)
zhat <- numeric(n)

for (i in seq_len(n)) {
  kn <- get.knnx(coords[-i,], matrix(coords[i,], nrow=1), k=1)
  zhat[i] <- z[-i][kn$nn.index[1]]
}

rmse_nn <- sqrt(mean((z - zhat)^2, na.rm = TRUE))
rmse_nn
## [1] 1.434082

7.6.4 Nearest Neighbor

# =========================================
# Nearest Neighbor — Manual (Konseptual)
# =========================================

# Ambil centroid dan koordinat
pts_sf <- indo_merged_scaled |>
  st_centroid() |>
  st_transform(4326)
## Warning: st_centroid assumes attributes are constant over geometries
coords <- st_coordinates(pts_sf)

pts <- data.frame(
  id = seq_len(nrow(pts_sf)),
  x  = coords[,1],
  y  = coords[,2],
  Z  = pts_sf$Y_Tindak_scaled
)


# Titik target (contoh: centroid provinsi ke-i)
i <- 1
s0 <- data.frame(x = pts$x[i], y = pts$y[i])

# Fungsi jarak Euclidean
edist <- function(x1, y1, x2, y2)
  sqrt((x1 - x2)^2 + (y1 - y2)^2)

# Hitung jarak ke semua titik lain
pts$dist_s0 <- edist(pts$x, pts$y, s0$x, s0$y)
pts$dist_s0[i] <- Inf   # hindari dirinya sendiri

# Urutkan
pts_ord <- pts[order(pts$dist_s0), ]

# Tetangga terdekat
nearest <- pts_ord[1, ]

Zhat <- nearest$Z
Zhat
## [1] 2.171746
# =========================================
# Nearest Neighbor (1-NN) — Grid Nasional
# =========================================

library(FNN)
library(sp)

# 1) Titik observasi (centroid provinsi)
indo_centroid <- indo_merged_scaled |>
  st_centroid() |>
  st_transform(4326)
## Warning: st_centroid assumes attributes are constant over geometries
sp_obs <- as(indo_centroid, "Spatial")
Xobs <- coordinates(sp_obs)

# 2) Grid interpolasi
r <- raster::raster(
  extent(sp_obs),
  res = 0.25,
  crs = CRS("+proj=longlat +datum=WGS84")
)

grid_sp <- as(r, "SpatialPixels")

# 3) Cari nearest neighbor
knn <- get.knnx(
  data  = Xobs,
  query = coordinates(grid_sp),
  k = 1
)

nn_ix <- knn$nn.index[,1]

# 4) Prediksi NN
pred_nn <- sp_obs$Y_Tindak_scaled[nn_ix]

nn_grid <- SpatialPixelsDataFrame(
  points = grid_sp,
  data = data.frame(Y_nn = pred_nn),
  proj4string = CRS(proj4string(sp_obs))
)

# 5) Plot
spplot(
  nn_grid, "Y_nn",
  main = "Nearest Neighbor (1-NN) — Tindak Kejahatan (Scaled)"
)

# =========================================
# RMSE Nearest Neighbor (LOOCV)
# =========================================

rmse_nn <- function(sf_data, var) {

  coords <- st_coordinates(st_centroid(sf_data))
  z <- sf_data[[var]]

  valid <- !is.na(z)
  coords <- coords[valid, , drop = FALSE]
  z <- z[valid]

  n <- length(z)
  zhat <- numeric(n)

  for (i in seq_len(n)) {

    coords_train <- coords[-i, , drop = FALSE]
    z_train <- z[-i]

    nn_i <- get.knnx(
      data  = coords_train,
      query = matrix(coords[i, ], nrow = 1),
      k = 1
    )$nn.index[1]

    zhat[i] <- z_train[nn_i]
  }

  sqrt(mean((z - zhat)^2))
}

rmse_nn_val <- rmse_nn(indo_merged_scaled, "Y_Tindak_scaled")
## Warning: st_centroid assumes attributes are constant over geometries
rmse_nn_val
## [1] 1.434082
Metode Interpolasi RMSE
Ordinary Kriging 0.993
Inverse Distance Weighting (IDW) 1.095766
Natural Neighbor / Nearest Neighbor 1.434082

Interpretasi :

Berdasarkan nilai RMSE, Ordinary Kriging merupakan metode interpolasi terbaik pada kasus tingkat kriminal tahun 2023 karena menghasilkan kesalahan prediksi paling kecil dibandingkan IDW serta Natural Neighbor/Nearest Neighbor. Hal ini menunjukkan bahwa Kriging lebih mampu menangkap struktur spasial dan ketergantungan antarwilayah pada data kriminalitas. Dengan demikian, hasil interpolasi Kriging lebih andal untuk menggambarkan pola spasial tingkat kriminal secara regional, sedangkan metode lain cenderung kurang akurat karena tidak secara eksplisit memodelkan korelasi spasial dalam data.

7.7 INTERPRETASI HASIL

sdem_model <- errorsarlm(
  Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + X3_RLS_scaled,
  data = indo_merged_scaled,
  listw = lwB,
  Durbin = TRUE,
  method = "Matrix",
  zero.policy = TRUE
)
summary(sdem_model)
## 
## Call:errorsarlm(formula = Y_Tindak_scaled ~ X1_UMR_scaled + X2_PDRB_scaled + 
##     X3_RLS_scaled, data = indo_merged_scaled, listw = lwB, Durbin = TRUE, 
##     method = "Matrix", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.26998 -0.57894 -0.19136  0.45413  2.48493 
## 
## Type: error 
## Regions with no neighbours included:
##  2 17 20 21 22 23 
## Coefficients: (asymptotic standard errors) 
##                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)        -0.064979   0.260887 -0.2491   0.80331
## X1_UMR_scaled      -1.086681   0.259775 -4.1832 2.875e-05
## X2_PDRB_scaled      0.638103   0.257217  2.4808   0.01311
## X3_RLS_scaled       0.184108   0.213857  0.8609   0.38930
## lag.(Intercept)     0.268582   0.112755  2.3820   0.01722
## lag.X1_UMR_scaled   1.170248   0.187281  6.2486 4.141e-10
## lag.X2_PDRB_scaled -1.130134   0.187848 -6.0162 1.785e-09
## lag.X3_RLS_scaled   0.285559   0.137191  2.0815   0.03739
## 
## Lambda: -0.5, LR test value: 36.609, p-value: 1.4437e-09
## Asymptotic standard error: 0.012525
##     z-value: -39.919, p-value: < 2.22e-16
## Wald statistic: 1593.5, p-value: < 2.22e-16
## 
## Log likelihood: -23.17923 for error model
## ML residual variance (sigma squared): 0.64214, (sigma: 0.80134)
## Number of observations: 34 
## Number of parameters estimated: 10 
## AIC: 66.358, (AIC for lm: 100.97)
ggplot(st_as_sf(ok)) +
  geom_sf(aes(fill = var1.pred), color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(title = "Ordinary Kriging — Prediksi Y Tindak Kriminal")

Model akhir yang digunakan adalah Spatial Durbin Error Model (SDEM) :

\[ \text{Tindak Kriminal}_i = \beta_0 + \beta_1 \text{UMR}_i + \beta_2 \text{PDRB}_i + \beta_3 \text{RLS}_i + \theta_1 W(\text{UMR})_i + \theta_2 W(\text{PDRB})_i + \theta_3 W(\text{RLS})_i + u_i \]

\[ u_i = \lambda W u_i + \varepsilon_i \]

karena memiliki nilai AIC terendah (66,358) dan log-likelihood tertinggi, sehingga paling sesuai untuk menjelaskan variasi tingkat tindak kriminal antarwilayah. Hasil estimasi menunjukkan bahwa UMR berpengaruh negatif dan signifikan, yang berarti peningkatan upah minimum cenderung menurunkan kriminalitas. Sebaliknya, PDRB berpengaruh positif dan signifikan, mengindikasikan bahwa wilayah dengan aktivitas ekonomi tinggi berpotensi menarik lebih banyak tindak kriminal. Rata-rata Lama Sekolah tidak signifikan, sehingga belum menunjukkan pengaruh langsung terhadap variasi kriminalitas.

Efek spasial tidak langsung dalam SDEM juga signifikan, yang tercermin dari koefisien variabel lag dan nilai lambda = −0,5 yang signifikan. Hal ini menegaskan adanya ketergantungan spasial pada komponen error, sehingga kriminalitas di suatu wilayah dipengaruhi pula oleh kondisi ekonomi wilayah tetangga. Dengan demikian, faktor-faktor tak teramati seperti mobilitas penduduk dan karakteristik regional berperan dalam membentuk pola kriminal secara spasial, sehingga kebijakan pengendalian kriminal perlu mempertimbangkan keterkaitan antarwilayah.

Hasil Ordinary Kriging menunjukkan pola sebaran kriminalitas yang membentuk klaster spasial, dengan wilayah tertentu memiliki tingkat prediksi kriminal yang lebih tinggi dibandingkan wilayah lain. Nilai RMSE sebesar 0,993 menandakan bahwa Kriging merupakan metode interpolasi paling akurat dibandingkan metode lain. Dalam konteks penelitian ini, SDEM menjelaskan hubungan dan mekanisme pengaruh spasial kriminalitas, sedangkan Kriging memvisualisasikan pola sebarannya, sehingga keduanya saling melengkapi dalam analisis kriminalitas berbasis spasial.

8 KESIMPULAN DAN SARAN

Berdasarkan hasil analisis, model Spatial Durbin Error Model (SDEM) merupakan model terbaik dalam menjelaskan variasi tindak kriminal antarwilayah karena mampu menangkap pengaruh faktor lokal (UMR, PDRB, dan RLS), pengaruh wilayah tetangga, serta ketergantungan spasial pada error. Hasil menunjukkan bahwa kondisi ekonomi suatu wilayah dan wilayah sekitarnya memiliki peran penting terhadap tingkat kriminalitas, sementara metode Ordinary Kriging memberikan hasil interpolasi spasial paling akurat dibandingkan metode lain. Oleh karena itu, disarankan agar perumusan kebijakan penanggulangan kriminalitas tidak hanya berfokus pada satu wilayah secara terpisah, tetapi juga mempertimbangkan keterkaitan antarwilayah, khususnya dalam aspek ekonomi dan pembangunan sosial.

9 LAMPIRAN

Daftar Pustaka

[1] Badan Pusat Statistik Indonesia. (12 Desember 2025). Jumlah Tindak Pidana Menurut Kepolisian Daerah, 2023. Diakses pada 21 Desember 2025, dari https://www.bps.go.id/id/statistics-table/2/MTAxIzI=/jumlah-tindak-pidana-menurut-kepolisian-daerah.html

[2] Badan Pusat Statistik Indonesia. (2 Desember 2025). Upah Rata - Rata Per Jam Pekerja Menurut Provinsi, 2023. Diakses pada 21 Desember 2025, dari https://www.bps.go.id/id/statistics-table/2/MTE3MiMy/upah-rata---rata-per-jam-pekerja-menurut-provinsi--rupiah-jam-.html

[3] Badan Pusat Statistik Indonesia. (28 Februari 2025). [Seri 2010] Produk Domestik Regional Bruto Per Kapita, 2023. Diakses pada 21 Desember 2025, dari https://www.bps.go.id/id/statistics-table/2/Mjg4IzI=/-seri-2010--produk-domestik-regional-bruto-per-kapita--ribu-rupiah-.html

[4] Badan Pusat Statistik Indonesia. (12 Desember 2025). Rata-Rata Lama Sekolah Penduduk Umur 15 Tahun ke Atas Menurut Provinsi, 2023. Diakses pada 21 Desember 2025, dari https://www.bps.go.id/id/statistics-table/2/MTQyOSMy/rata-rata-lama-sekolah-penduduk-umur-15-tahun-ke-atas-menurut-provinsi.html

[5] Sari TNS dan Yulianto S. Analisis Metode Spasial untuk Menentukan Faktor-Faktor yang Mempengaruhi Kriminalitas pada Provinsi Jawa Timur Tahun 2022. Jurnal Ilmiah Dinamika Sosial. Artikel ini menggunakan pendekatan spasial dan regresi spasial untuk memodelkan tingkat kriminalitas dengan variabel sosial ekonomi serta hubungan antarwilayah menggunakan bobot Queen Contiguity. https://journal.undiknas.ac.id/index.php/fisip/article/view/5326

[6] Trisnawati & Khoirunurrofik. Inter-Provincial Spatial Linkages of Crime Pattern in Indonesia: Looking at Education and Economic Inequality Effects on Crime. Indonesian Journal of Geography. Penelitian ini menelaah keterkaitan spasial tingkat kriminal antarprovinsi di Indonesia, termasuk pengaruh pendidikan dan ketimpangan ekonomi terhadap kejadian kriminal, serta menunjukkan hubungan spasial antar wilayah. https://jurnal.ugm.ac.id/ijg/article/view/34026

[7] LeSage, J. dan Pace, R. K. Introduction to Spatial Econometrics. CRC Press, 2009.

https://www.crcpress.com/Introduction-to-Spatial-Econometrics/LeSage-Pace/p/book/9781584886105

[8] Elhorst, J. P. Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Springer, 2014.

https://link.springer.com/book/10.1007/978-3-319-23536-7

[9] Anselin, L. Local Indicators of Spatial Association LISA. Geographical Analysis, 1995.

https://dces.webhosting.cals.wisc.edu/wp-content/uploads/sites/128/2013/08/W4_Anselin1995.pdf

[10] Getis, A. dan Ord, J. K. The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis, 1992.

https://onlinelibrary.wiley.com/doi/10.1111/j.1538-4632.1992.tb00261.x

[11] Bivand, R., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R. Springer.

[12] Cressie, N. (1993). Statistics for Spatial Data. Wiley.