LAPORAN UJIAN AKHIR SEMESTER
Determinasi Indeks Pembangunan Manusia Kabupaten/Kota Jawa Timur Melalui Pendekatan Ekonometrika Spasial

Disusun oleh :
Sefia Herlina (140610230003)

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

PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

Abstrak

Pembangunan manusia merupakan indikator kunci keberhasilan pembangunan daerah karena mencerminkan kualitas hidup masyarakat secara menyeluruh, tidak hanya dari sisi ekonomi tetapi juga pendidikan dan kesejahteraan sosial. Penelitian bertujuan untuk menganalisis pola persebaran spasial Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur, menguji keberadaan ketergantungan spasial (spatial dependence) pada IPM, mengidentifikasi pengaruh Produk Domestik Regional Bruto (PDRB), Rata-Rata Lama Sekolah (RLS), Harapan Lama Sekolah (HLS), dan jumlah penduduk miskin terhadap IPM, serta menentukan model ekonometrika spasial yang paling tepat dalam menjelaskan determinasi IPM. Metode yang digunakan adalah pendekatan kuantitatif dengan data cross-section 38 kabupaten/kota tahun 2024 yang bersumber dari Badan Pusat Statistik dan data spasial administrasi. Analisis dilakukan melalui integrasi data spasial dan statistik, analisis deskriptif dan pemetaan tematik, pengujian autokorelasi spasial menggunakan Moran’s I, Geary’s C, Getis-Ord Gi*, dan Local Moran’s I (LISA), serta estimasi dan perbandingan berbagai model ekonometrika spasial, yaitu OLS, SAR, SEM, SDM, SDEM, SAC, dan GNS, dengan pemilihan model terbaik berdasarkan uji Lagrange Multiplier (LM) serta kriteria informasi AIC dan BIC. Sebagai analisis tambahan digunakan Geographically Weighted Regression (GWR) dan metode interpolasi spasial. Hasil penelitian menunjukkan bahwa IPM di Jawa Timur tidak tersebar secara acak dan membentuk pengelompokan spasial yang signifikan, meskipun secara lokal klaster kuat hanya muncul pada wilayah tertentu. Estimasi model menunjukkan bahwa RLS dan HLS berpengaruh positif dan signifikan terhadap IPM, sedangkan PDRB dan jumlah penduduk miskin tidak signifikan setelah mempertimbangkan efek spasial. Berdasarkan perbandingan model, Spatial Autoregressive Combined (SAC) merupakan model terbaik dengan nilai AIC dan BIC terendah serta koefisien ketergantungan spasial yang signifikan, yang mengindikasikan bahwa IPM dipengaruhi oleh interaksi antarwilayah dan faktor tak teramati yang berpola spasial. Temuan ini mengimplikasikan bahwa kebijakan pembangunan manusia di Provinsi Jawa Timur perlu dirancang secara terintegrasi berbasis wilayah, dengan penekanan pada peningkatan kualitas pendidikan dan penguatan sinergi antar daerah guna mengurangi ketimpangan pembangunan manusia secara berkelanjutan.

Kata kunci: Indeks Pembangunan Manusia; Ekonometrika Spasial; Autokorelasi Spasial; Spatial Autoregressive Combined (SAC); Provinsi Jawa Timur

Daftar Isi

BAB I (PENDAHULUAN)

1.1 Latar belakang

Pembangunan manusia merupakan inti dari pembangunan berkelanjutan karena keberhasilannya tidak hanya diukur melalui pertumbuhan ekonomi, tetapi juga melalui peningkatan kualitas hidup masyarakat secara menyeluruh. Paradigma pembangunan modern menempatkan manusia sebagai tujuan akhir pembangunan, bukan sekadar sebagai alat untuk mencapai pertumbuhan ekonomi. Oleh karena itu, pembangunan yang berorientasi pada manusia menekankan pentingnya peningkatan kapasitas individu agar mampu menjalani kehidupan yang sehat, berpendidikan, dan sejahtera, serta memiliki kesempatan untuk berpartisipasi secara produktif dalam aktivitas sosial dan ekonomi. United Nations Development Programme (UNDP) memperkenalkan Indeks Pembangunan Manusia (IPM) atau Human Development Index (HDI) sebagai indikator komposit untuk mengukur tingkat keberhasilan pembangunan manusia. IPM adalah ukuran pencapaian rata-rata suatu wilayah dalam dimensi kesehatan, pendidikan, dan ekonomi yang merefleksikan kualitas hidup penduduknya.

Provinsi Jawa Timur merupakan salah satu provinsi dengan kontribusi ekonomi terbesar di Pulau Jawa, dengan konstribusi sekitar ±25,55% terhadap perekonomian Pulau Jawa pada triwulan III tahun 2024. Besarnya kontribusi tersebut mencerminkan peran strategis Jawa Timur sebagai motor pertumbuhan ekonomi regional. Namun demikian, capaian pembangunan manusia di provinsi ini belum sepenuhnya menunjukkan kondisi yang merata antar wilayah. Pada tahun 2024, Indeks Pembangunan Manusia (IPM) Provinsi Jawa Timur tercatat sebesar 75,31, yang berada dalam kategori tinggi berdasarkan klasifikasi BPS. Meskipun demikian, nilai IPM tersebut merupakan capaian rata-rata provinsi yang menyembunyikan adanya disparitas pembangunan manusia yang cukup signifikan antar kabupaten/kota. Ketimpangan tersebut tercermin dari perbedaan nilai IPM yang mencolok antar wilayah, seperti Kota Surabaya dengan nilai IPM tercatat sekitar 84,69, yang menunjukkan kualitas hidup masyarakat yang relatif lebih baik. Namun, beberapa wilayah kabupaten masih berada jauh di bawah rata-rata provinsi seperti Kabupaten Sampang yang tercatat nilai IPM sekitar 66,72, yang berada pada kategori menengah dan tertinggal cukup jauh dibandingkan daerah lainnya. Perbedaan ini mengindikasikan bahwa tingginya kinerja ekonomi di tingkat provinsi belum sepenuhnya diikuti oleh pemerataan kualitas pembangunan manusia di seluruh wilayah.

Gambar 1. Peta Persebaran Indeks Pembangunan Manusia (IPM) per Kabupaten/Kota di Provinsi Jawa Timur

Berdasarkan Gambar 1, terlihat bahwa nilai Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur menunjukkan pola pengelompokan wilayah dengan capaian IPM yang relatif tinggi dan rendah. Distribusi IPM tersebut tidak tersebar secara acak, melainkan membentuk klaster spasial tertentu. Wilayah dengan IPM tinggi hingga sangat tinggi cenderung terkonsentrasi di kawasan perkotaan dan pusat pertumbuhan ekonomi, sementara wilayah dengan IPM sedang hingga rendah lebih banyak ditemukan di daerah perdesaan. Pola pengelompokan ini mengindikasikan adanya ketergantungan spasial dalam pembangunan manusia, di mana capaian IPM suatu wilayah tidak hanya ditentukan oleh karakteristik internal wilayah tersebut, tetapi juga oleh kondisi sosial ekonomi wilayah sekitarnya. Kedekatan geografis antar kabupaten/kota berpotensi mendorong kesamaan capaian pembangunan manusia melalui berbagai mekanisme interaksi. Dengan demikian, ketimpangan pembangunan manusia di Jawa Timur menjadi permasalahan serius yang harus segera ditangani.

Ketimpangan pembangunan manusia ini menjadi isu strategis dalam pencapaian Sustainable Development Goals (SDGs), khususnya Tujuan 1 (Tanpa Kemiskinan), Tujuan 4 (Pendidikan Berkualitas), Tujuan 8 (Pekerjaan Layak dan Pertumbuhan Ekonomi), serta Tujuan 10 (Mengurangi Ketimpangan). Ketidakmerataan IPM menunjukkan bahwa manfaat pembangunan belum dirasakan secara seimbang oleh seluruh wilayah, sehingga berpotensi menghambat pencapaian target pembangunan yang bersifat inklusif dan berkelanjutan. Jika ketimpangan ini tidak ditangani secara tepat, maka kemajuan pembangunan di tingkat provinsi maupun nasional berisiko bersifat parsial dan tidak merata. Oleh karena itu, pencapaian SDGs menuntut kebijakan pembangunan yang tidak hanya berorientasi pada peningkatan capaian agregat, tetapi juga pada pemerataan kualitas pembangunan manusia antarwilayah. Pendekatan pembangunan berbasis wilayah menjadi penting untuk mengidentifikasi daerah-daerah yang tertinggal serta memahami keterkaitan antarwilayah dalam proses pembangunan. Dengan mempertimbangkan karakteristik dan interaksi spasial, kebijakan yang dirumuskan diharapkan mampu mengurangi ketimpangan pembangunan manusia dan mendorong pencapaian SDGs secara lebih efektif, adil, dan berkelanjutan.

Sejumlah penelitian terdahulu telah mengkaji faktor-faktor yang memengaruhi Indeks Pembangunan Manusia (IPM) di Provinsi Jawa Timur dengan menggunakan berbagai pendekatan analisis. Penelitian oleh Santoso dkk. (2022) menunjukkan adanya dependensi spasial dalam pembentukan IPM, di mana capaian IPM suatu kabupaten/kota tidak hanya dipengaruhi oleh karakteristik internal wilayah tersebut, tetapi juga oleh kondisi wilayah sekitarnya melalui komponen error spasial. Temuan tersebut menegaskan bahwa pembangunan manusia memiliki dimensi kewilayahan yang kuat dan pendekatan ekonometrika konvensional yang mengabaikan interaksi spasial berpotensi menghasilkan estimasi yang bias dan kurang akurat. Komponen pengaruh IPM, seperti Produk Domestik Regional Bruto (PDRB) yang mencerminkan kapasitas dan kinerja ekonomi suatu wilayah melalui nilai tambah barang dan jasa yang dihasilkan. Rata-Rata Lama Sekolah (RLS) yang menggambarkan tingkat pendidikan aktual penduduk berusia 15 tahun ke atas. Harapan Lama Sekolah (HLS) yang mencerminkan peluang pendidikan yang diharapkan dapat ditempuh oleh penduduk usia sekolah di masa mendatang. Jumlah penduduk miskin yang menunjukkan tingkat keterbatasan ekonomi masyarakat berdasarkan pengeluaran per kapita di bawah garis kemiskinan. Meskipun variabel-variabel tersebut telah banyak digunakan dalam penelitian IPM, sebagian besar studi masih menerapkan pendekatan ekonometrika konvensional atau hanya membandingkan model spasial tertentu, seperti Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM). Padahal, setiap spesifikasi model spasial merepresentasikan mekanisme ketergantungan antarwilayah yang berbeda, baik melalui variabel dependen, variabel independen, maupun komponen error. Oleh karena itu, diperlukan analisis yang lebih komprehensif dengan membandingkan berbagai model ekonometrika spasial untuk memperoleh pemahaman yang lebih akurat mengenai determinasi IPM.

Untuk itu pada penelitian ini menggunakan analisis dengan membandingkan beberapa model ekonometrika spasial, yaitu Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS). Pemilihan model terbaik dilakukan secara sistematis menggunakan Lagrange Multiplier (LM) test untuk mendeteksi ketergantungan spasial, serta Akaike Information Criterion (AIC) dan Bayesian Information Criterion (BIC) untuk mengevaluasi kinerja dan kesesuaian model. Dengan demikian, penelitian ini diharapkan mampu memberikan gambaran yang lebih komprehensif mengenai faktor-faktor penentu pembangunan manusia serta dinamika spasial antar kabupaten/kota di Provinsi Jawa Timur, yang dapat menjadi dasar bagi perumusan kebijakan pembangunan manusia berbasis wilayah.

1.2 Identifikasi Masalah

Berdasarkan hal tersebut, penelitian ini mengidentifikasi permasalahan sebagai berikut:

Bagaimana pola persebaran spasial Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur?
Apakah terdapat ketergantungan spasial (spatial dependence) IPM antarwilayah di Provinsi Jawa Timur?
Bagaimana pengaruh Produk Domestik Regional Bruto (PDRB), Rata-Rata Lama Sekolah (RLS), Harapan Lama Sekolah (HLS), dan jumlah penduduk miskin terhadap IPM di Provinsi Jawa Timur?
Model ekonometrika spasial manakah yang paling tepat dalam menjelaskan determinasi IPM di Provinsi Jawa Timur?

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk:

Menganalisis pola persebaran spasial IPM antar kabupaten/kota di Provinsi Jawa Timur.
Menguji adanya ketergantungan spasial IPM antarwilayah di Provinsi Jawa Timur.
Menganalisis pengaruh PDRB, Rata-Rata Lama Sekolah, Harapan Lama Sekolah, dan jumlah penduduk miskin terhadap IPM di Provinsi Jawa Timur.
Menentukan model ekonometrika spasial terbaik dalam menjelaskan determinasi IPM di Provinsi Jawa Timur.

1.4 Batasan Penelitian

Untuk memperjelas ruang lingkup penelitian, maka batasan penelitian ini adalah sebagai berikut:

Penelitian dilakukan pada 38 kabupaten/kota di Provinsi Jawa Timur.
Variabel yang digunakan terdiri atas variabel dependen, yaitu Indeks Pembangunan Manusia (IPM) dan variabel independen, yaitu PDRB, Rata-Rata Lama Sekolah, Harapan Lama Sekolah, dan jumlah penduduk miskin..
Data yang digunakan merupakan data sekunder yang diperoleh dari Badan Pusat Statistik (BPS).
Penelitian menggunakan pendekatan cross-section dengan analisis ekonometrika spasial.
Model spasial yang dibandingkan meliputi SAR, SEM, SDM, SDEM, SAC, dan GNS, dengan pemilihan model terbaik berdasarkan uji LM, AIC, dan BIC.

BAB II (TINJAUAN PUSTAKA)

2.1 Index Pembangunan Manusia (IPM)

Indeks Pembangunan Manusia (IPM) merupakan indikator komposit yang dikembangkan oleh United Nations Development Programme (UNDP) untuk mengukur tingkat pencapaian pembangunan manusia secara multidimensi. Konsep IPM berangkat dari paradigma pembangunan manusia yang menempatkan manusia sebagai tujuan utama pembangunan, bukan sekadar sebagai instrumen pertumbuhan ekonomi. IPM mengintegrasikan tiga dimensi dasar, yaitu umur panjang dan hidup sehat yang direpresentasikan oleh angka harapan hidup, pengetahuan yang diukur melalui rata-rata lama sekolah dan harapan lama sekolah, serta standar hidup layak yang dicerminkan oleh kemampuan daya beli masyarakat. Dengan pendekatan tersebut, IPM memberikan kerangka analitis untuk menilai kualitas hidup penduduk serta membandingkan tingkat pembangunan manusia antarwilayah dan antarwaktu.

2.2 Teori Dasar Spatial Dependence

Spatial dependence atau ketergantungan spasial merupakan konsep utama dalam analisis spasial yang menjelaskan bahwa nilai suatu variabel pada suatu lokasi tidak bersifat independen, melainkan dipengaruhi oleh nilai variabel pada lokasi-lokasi lain yang berdekatan secara geografis. Ketergantungan spasial dapat disebabkan oleh interaksi ekonomi regional, seperti perdagangan antar daerah, aglomerasi aktivitas ekonomi, serta penyebaran manfaat pembangunan dari suatu wilayah ke wilayah sekitarnya (spillover effects). Selain itu, faktor sosial seperti migrasi penduduk, kesamaan budaya, dan akses bersama terhadap fasilitas publik juga berkontribusi terhadap terbentuknya ketergantungan spasial. Dalam konteks pembangunan manusia, capaian Indeks Pembangunan Manusia (IPM) suatu wilayah dapat dipengaruhi oleh ketersediaan fasilitas pendidikan dan kesehatan di wilayah sekitar, sehingga pembangunan manusia bersifat saling terkait secara spasial.

Secara metodologis, keberadaan spatial dependence memiliki implikasi penting terhadap pemodelan statistik dan ekonometrika. Model regresi klasik mengasumsikan bahwa observasi bersifat independen satu sama lain. Namun, ketika spatial dependence hadir, asumsi ini menjadi tidak terpenuhi, sehingga estimasi parameter dengan metode Ordinary Least Squares (OLS) berpotensi menjadi bias dan tidak efisien. Namun, pengabaian ketergantungan spasial dapat menghasilkan kesimpulan statistik yang menyesatkan dan mengurangi validitas inferensi. Oleh karena itu, diperlukan pendekatan ekonometrika spasial yang secara eksplisit memasukkan struktur ketergantungan spasial ke dalam model. Dalam kerangka ekonometrika spasial, spatial dependence umumnya dimodelkan melalui dua bentuk utama, yaitu spatial lag dependence dan spatial error dependence. Spatial lag dependence terjadi ketika nilai variabel dependen pada suatu wilayah dipengaruhi secara langsung oleh nilai variabel dependen di wilayah sekitarnya, sebagaimana direpresentasikan dalam model Spatial Autoregressive (SAR). Sementara itu, spatial error dependence terjadi ketika ketergantungan spasial muncul melalui komponen error, yang mencerminkan adanya faktor-faktor tidak terobservasi yang bersifat spasial, sebagaimana dimodelkan dalam Spatial Error Model (SEM). Perkembangan lebih lanjut dari teori ini melahirkan berbagai model hibrida, seperti Spatial Durbin Model (SDM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS), yang memungkinkan pemodelan ketergantungan spasial secara lebih komprehensif.

2.3 Autokorelasi Spasial

Autokorelasi Spasial merupakan konsep yang menggambarkan sejauh mana nilai suatu variabel di suatu lokasi berhubungan dengan nilai variabel yang sama di lokasi lain yang berdekatan secara geografis. Secara sederhana, fenomena ini menunjukkan adanya pola keterkaitan antarwilayah, di mana wilayah yang berdekatan cenderung memiliki karakteristik yang mirip (autokorelasi positif), atau justru sangat berbeda (autokorelasi negatif). Ketika nilai suatu variabel di satu daerah memengaruhi nilai di daerah lain, maka data tersebut tidak lagi bersifat independen secara spasial.

Secara matematis, spatial autocorrelation dapat dijelaskan melalui kovarians antara dua observasi yi dan yj yang berjarak dij, dengan Cov(yi, yj ) ≠ 0. Kondisi ini melanggar asumsi dasar regresi OLS yang mengharuskan komponen error antar observasi bersifat bebas (independent and identically distributed). Jika ketergantungan spasial diabaikan, maka hasil estimasi akan menjadi bias dan tidak efisien.

Untuk mendeteksi keberadaan autokorelasi spasial, digunakan berbagai ukuran statistik yang dapat diklasifikasikan ke dalam ukuran global dan ukuran lokal. Ukuran autokorelasi spasial global digunakan untuk mengidentifikasi apakah secara keseluruhan terdapat pola keterkaitan spasial dalam data, tanpa memperhatikan lokasi spesifik dari pengelompokan tersebut. Beberapa ukuran global yang umum digunakan antara lain Moran’s I, Geary’s C, dan Getis-Ord G. Nilai Moran’s I yang positif menunjukkan adanya kecenderungan pengelompokan wilayah dengan nilai variabel yang serupa (positive spatial autocorrelation), sedangkan nilai negatif mengindikasikan pola penyebaran atau ketidaksamaan antarwilayah. Sementara itu, Geary’s C lebih sensitif terhadap perbedaan nilai antarwilayah yang berdekatan, dan Getis-Ord G digunakan untuk mengidentifikasi kecenderungan konsentrasi nilai tinggi atau rendah secara global.Namun, ukuran autokorelasi spasial global hanya mampu menunjukkan keberadaan autokorelasi secara keseluruhan, tanpa mengidentifikasi lokasi spesifik terjadinya pengelompokan spasial. Oleh karena itu, analisis perlu dilengkapi dengan ukuran autokorelasi spasial lokal, yang dikenal sebagai Local Indicators of Spatial Association (LISA). LISA digunakan untuk mengidentifikasi klaster spasial pada tingkat lokal serta mendeteksi keberadaan hot spots dan cold spots, seperti klaster High–High, Low–Low, serta pola pencilan spasial High–Low dan Low–High. Dengan demikian, LISA memungkinkan pemahaman yang lebih mendalam mengenai variasi spasial dan lokasi spesifik terjadinya ketergantungan antarwilayah.

2.3.1 Moran’s I

Uji Moran’s I Global digunakan untuk mendeteksi keberadaan autokorelasi spasial dalam distribusi IPM antarkabupaten/kota. Rumusnya sebagai berikut:

\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij} (x_i - \bar{x})(x_j - \bar{x})} {\sum_i (x_i - \bar{x})^2} \]

Keterangan:
- \(n\) = jumlah unit observasi
- \(x_i, x_j\) = nilai variabel di lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata nilai variabel
- \(w_{ij}\) = elemen matriks pembobot spasial

Kriteria Pengujian:
- \(I > 0\): autokorelasi spasial positif (wilayah dengan nilai serupa saling berdekatan)
- \(I < 0\): autokorelasi spasial negatif (wilayah dengan nilai berbeda berdekatan)
- \(I \approx 0\): tidak terdapat autokorelasi spasial (pola acak)

2.3.2 Geary’s C

Statistik Geary’s C diperkenalkan oleh R. C. Geary pada tahun 1954 untuk mengukur autokorelasi spasial, namun lebih sensitif terhadap perbedaan lokal antar pasangan lokasi. Rumusnya sebagai berikut:

\[ C = \frac{(n - 1)}{2 \sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij} (x_i - x_j)^2} {\sum_i (x_i - \bar{x})^2} \]

Keterangan:
- \(n\) = jumlah unit observasi
- \(x_i, x_j\) = nilai variabel di lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata nilai variabel
- \(w_{ij}\) = elemen matriks pembobot spasial

Kriteria Pengujian:
- \(C = 1\): tidak ada autokorelasi spasial
- \(C < 1\): terdapat autokorelasi spasial positif
- \(C > 1\): terdapat autokorelasi spasial negatif

2.3.3 Getis-Ord G

Statistik Getis-Ord G (atau sering disebut Getis-Ord General G) dikembangkan oleh Getis dan Ord (1992) untuk mengukur autokorelasi spasial global maupun lokal dengan fokus pada pengelompokan nilai tinggi (hot spots) dan nilai rendah (cold spots). Rumusnya sebagai berikut:

\[ G = \frac{\sum_i \sum_j w_{ij} x_i x_j} {\sum_i \sum_j x_i x_j} \]

Keterangan:
- \(x_i, x_j\) = nilai variabel pada lokasi \(i\) dan \(j\)
- \(w_{ij}\) = bobot spasial antara lokasi \(i\) dan \(j\)

Kriteria Pengujian:
- Nilai \(G\) tinggi menunjukkan hot spot (pengelompokan nilai besar).
- Nilai \(G\) rendah menunjukkan cold spot (pengelompokan nilai kecil).
- Nilai \(G\) mendekati rata-rata global menunjukkan pola acak.

Terdapat juga versi lokal dari statistik ini, yaitu \(Getis-Ord Gi^*\), yang digunakan untuk mendeteksi lokasi-lokasi spesifik yang menjadi pusat hot spot atau cold spot.

2.3.4 Local Moran’s I

Statistik Local Moran’s I dikembangkan oleh Luc Anselin (1995) untuk mengidentifikasi klaster spasial pada tingkat lokal serta mendeteksi keberadaan hot spots dan cold spots, seperti klaster High–High, Low–Low, serta pola pencilan spasial High–Low dan Low–High.. Rumusnya sebagai berikut:

\[ I_i = \frac{x_i - \bar{x}}{m_2} \sum_{j=1}^{N} w_{ij}(x_j - \bar{x}), \]

dengan

\[ m_2 = \frac{1}{N} \sum_{k=1}^{N} (x_k - \bar{x})^2. \]

Kriteria Pengujian:
- Jika \(I_i > 0\), maka unit spasial ke-\(i\) memiliki nilai yang mirip dengan wilayah tetangganya, sehingga membentuk klaster High–High atau Low–Low. - Jika \(I_i < 0\), maka unit spasial ke-\(i\) memiliki nilai yang berbeda dengan wilayah tetangganya, yang menunjukkan pola High–Low atau Low–High. - Jika \(I_i \approx 0\), maka tidak terdapat pola spasial lokal yang signifikan atau pola bersifat acak.

Temuan adanya autokorelasi spasial yang signifikan, baik secara global maupun lokal, mengindikasikan bahwa asumsi independensi antarobservasi tidak terpenuhi. Kondisi ini menegaskan perlunya penerapan model ekonometrika spasial yang secara eksplisit mampu mengakomodasi hubungan dan interaksi antarwilayah, sehingga menghasilkan estimasi parameter yang lebih akurat dan inferensi statistik yang lebih valid.

2.4 Estimasi Model Spasial Ekonometrik

Model ekonometrika spasial dikembangkan untuk mengatasi permasalahan spatial dependence dan spatial spillover dalam data yang memiliki dimensi geografis. Anselin (1988) dan LeSage dan Pace (2009) menegaskan bahwa pengabaian ketergantungan spasial dalam model regresi konvensional dapat menghasilkan estimasi parameter yang bias, tidak efisien, serta kesimpulan statistik yang menyesatkan. Oleh karena itu, model ekonometrika spasial digunakan untuk secara eksplisit memasukkan struktur ketergantungan spasial ke dalam proses estimasi. Bentuk paling dasar dari model ekonometrika spasial adalah Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM). Model SAR mengakomodasi ketergantungan spasial melalui variabel dependen yang terlambat secara spasial (spatial lag), di mana nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Model ini sesuai digunakan ketika interaksi antarwilayah secara langsung memengaruhi variabel yang diteliti. Sementara itu, model SEM menangkap ketergantungan spasial melalui komponen error, yang mencerminkan adanya faktor-faktor tidak teramati yang bersifat spasial dan memengaruhi variabel dependen secara tidak langsung.

Perkembangan lebih lanjut dari model SAR dan SEM melahirkan model yang lebih fleksibel, salah satunya adalah Spatial Durbin Model (SDM). Model SDM memperluas model SAR dengan memasukkan spatial lag dari variabel independen, sehingga memungkinkan identifikasi pengaruh langsung dan tidak langsung (spillover effects) dari variabel penjelas. Selain itu, terdapat pula Spatial Durbin Error Model (SDEM) yang mengombinasikan spatial lag dari variabel independen dengan komponen error spasial. Model ini digunakan ketika pengaruh spasial terutama berasal dari variabel penjelas dan faktor tidak terobservasi, tetapi tidak melalui variabel dependen secara langsung. Selanjutnya, Spatial Autoregressive Combined Model (SAC) merupakan gabungan antara model SAR dan SEM, di mana ketergantungan spasial dimodelkan baik melalui variabel dependen maupun komponen error secara simultan. Model ekonometrika spasial yang paling umum dan fleksibel adalah General Nested Spatial Model (GNS). Model ini mencakup seluruh bentuk ketergantungan spasial, baik melalui variabel dependen, variabel independen, maupun komponen error. GNS dianggap sebagai model paling lengkap karena SAR, SEM, SDM, SDEM, dan SAC dapat dipandang sebagai bentuk khusus (restricted forms) dari model GNS.

Dengan menggunakan berbagai spesifikasi model ekonometrika spasial tersebut, analisis dapat dilakukan secara lebih komprehensif untuk memahami pola dan mekanisme ketergantungan spasial dalam data. Pemilihan model yang paling tepat menjadi aspek krusial, karena setiap model merepresentasikan asumsi dan mekanisme spasial yang berbeda. Oleh karena itu, perbandingan antar model melalui uji statistik dan kriteria informasi menjadi langkah penting dalam memperoleh estimasi yang akurat dan interpretasi yang bermakna terhadap fenomena spasial yang dianalisis.

2.4.1 Spatial Autoregressive Model (SAR)

Spatial Autoregressive Model (SAR) digunakan ketika nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah lain yang berdekatan. Artinya, terdapat efek limpahan (spillover effect) antarwilayah yang menyebabkan suatu daerah tidak dapat dipelajari secara terpisah. Model SAR dapat ditulis sebagai berikut:

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

Keterangan:

\(y\) = vektor variabel dependen (misalnya IPM)
\(X\) = matriks variabel independen (misalnya PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(\beta\) = vektor koefisien regresi
\(W\) = matriks pembobot spasial (spatial weights matrix)
\(\rho\) = koefisien autokorelasi spasial (spatial autoregressive coefficient)
\(\varepsilon\) = vektor error, dengan asumsi \(\varepsilon \sim N(0, \sigma^2 I)\)

Koefisien \(\rho\) menunjukkan sejauh mana pengaruh spasial terjadi antarwilayah. Jika \(\rho\) signifikan, maka variabel dependen di suatu wilayah dipengaruhi oleh nilai di wilayah tetangga. Dalam konteks penelitian IPM kabupaten/kota Jawa Timur, model SAR relevan apabila IPM di satu daerah dipengaruhi oleh kondisi IPM di daerah sekitar.

2.4.2 Spatial Error Model (SEM)

Berbeda dengan SAR, Spatial Error Model (SEM) digunakan ketika ketergantungan spasial tidak terjadi pada variabel dependen, tetapi pada komponen error. Artinya, pengaruh antarwilayah muncul akibat faktor yang tidak terobservasi namun saling berkaitan secara spasial, seperti kebijakan daerah, infrastruktur, atau kondisi sosial ekonomi yang mirip. Model SEM dinyatakan sebagai berikut:

\[ y = X\beta + (I - \lambda W)^{-1} \varepsilon \]

Keterangan:

\(y\) = vektor variabel dependen (IPM)
\(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(\beta\) = vektor koefisien regresi
\(W\) = matriks pembobot spasial
\(\lambda\) = koefisien autokorelasi spasial pada komponen error
\(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika \(\lambda\) signifikan, maka terdapat hubungan spasial dalam error term, yang menandakan bahwa variabel-variabel yang tidak dimasukkan ke dalam model berhubungan antarwilayah. Dalam konteks IPM, SEM cocok digunakan ketika faktor-faktor eksternal memiliki keterkaitan yang menyebabkan error model saling berhubungan.

2.4.3 Spatial Durbin Model (SDM)

Spatial Durbin Model digunakan untuk menganalisis pengaruh langsung dan tidak langsung (spillover effect) dari variabel independen antarwilayah dalam data spasial. Model SDM dinyatakan sebagai berikut:

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

Keterangan:

\(y\) = vektor variabel dependen (IPM)
\(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(W\) = matriks pembobot spasial
\(\rho\) = koefisien autoregresif spasial yang menunjukkan pengaruh variabel dependen di wilayah sekitar terhadap wilayah ke-\(i\).
\(\beta\) = vektor koefisien regresi
\(WX\) = spatial lag dari variabel independen
\(\theta\) = vektor koefisien spatial lag variabel independen
\(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika koefisien spasial \(\rho\) dan koefisien spatial lag variabel independen (\(\theta\)) signifikan, maka terdapat pengaruh spasial yang bersumber dari interaksi antarwilayah, baik melalui variabel dependen maupun variabel independen. Hal ini menandakan bahwa perubahan nilai variabel penjelas di suatu wilayah tidak hanya memengaruhi variabel dependen di wilayah tersebut, tetapi juga berdampak pada wilayah-wilayah sekitarnya (spillover effect). Dalam konteks IPM, SDM cocok digunakan ketika faktor-faktor seperti PDRB, tingkat pendidikan, dan kemiskinan di suatu kabupaten/kota berpotensi memengaruhi capaian IPM di kabupaten/kota lain melalui keterkaitan ekonomi, mobilitas penduduk, serta akses lintas wilayah terhadap layanan pendidikan dan kesehatan.

2.4.4 Spatial Durbin Error Model (SDEM)

Spatial Durbin Error Model digunakan untuk menganalisis pengaruh spasial yang bersumber dari variabel independen dan faktor-faktor tidak terobservasi yang saling berkaitan antarwilayah, tanpa melibatkan ketergantungan langsung pada variabel dependen. Model SDM dinyatakan sebagai berikut:

\[ y = X{\beta} + WX{\theta} + u, \]

dengan komponen error spasial:

\[ u = \lambda Wu + {\varepsilon}. \]

Keterangan:

\(y\) = vektor variabel dependen (IPM)
\(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(W\) = matriks pembobot spasial
\(\lambda\) = koefisien error spasial yang menunjukkan tingkat ketergantungan antar error wilayah
\(\beta\) = vektor koefisien regresi
\(WX\) = spatial lag dari variabel independen
\(\theta\) = vektor koefisien spatial lag variabel independen
\(u\) = komponen error yang mengandung ketergantungan spasial
\(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika koefisien spatial lag variabel independen (\(\theta\)) dan koefisien error spasial (\(\lambda\)) signifikan, maka terdapat pengaruh spasial yang bersumber dari variabel penjelas dan faktor-faktor tidak terobservasi yang bersifat spasial. Hal ini menandakan bahwa perubahan nilai variabel independen di suatu wilayah tidak hanya memengaruhi variabel dependen di wilayah tersebut, tetapi juga berdampak pada wilayah-wilayah sekitarnya melalui mekanisme spillover effect. Selain itu, signifikansi koefisien \(\lambda\) mengindikasikan bahwa terdapat keterkaitan antarwilayah yang berasal dari faktor-faktor eksternal yang tidak dimasukkan ke dalam model, sehingga komponen error antar wilayah saling berhubungan.

2.4.5 Spatial Autoregressive Combined (SAC)

Spatial Autoregressive Combined digunakan untuk menganalisis ketergantungan spasial yang terjadi secara simultan melalui variabel dependen dan komponen error. Model SDM dinyatakan sebagai berikut:

\[ y = \rho Wy + X{\beta} + u, \]

dengan komponen error spasial:

\[ u = \lambda Wu + {\varepsilon}. \]

Keterangan:

\(y\) = vektor variabel dependen (IPM)
\(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(W\) = matriks pembobot spasial
\(\rho\) = koefisien autoregresif spasial yang menunjukkan pengaruh variabel dependen di wilayah sekitar terhadap wilayah ke-\(i\).
\(\beta\) = vektor koefisien regresi
\(WX\) = spatial lag dari variabel independen
\(u\) = komponen error yang mengandung ketergantungan spasial
\(\lambda\) = koefisien error spasial yang menunjukkan tingkat ketergantungan antar error wilayah
\(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika koefisien spasial \(\rho\) dan koefisien error spasial \(\lambda\) signifikan, maka terdapat ketergantungan spasial yang terjadi secara simultan melalui variabel dependen dan komponen error. Signifikansi \(\rho\) menunjukkan adanya interaksi langsung antarwilayah, yaitu nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Sementara itu, signifikansi \(\lambda\) mengindikasikan bahwa komponen error antarwilayah saling berkaitan, yang menandakan adanya faktor-faktor eksternal atau variabel yang tidak dimasukkan ke dalam model tetapi memiliki pola spasial, sehingga memengaruhi variabel dependen secara tidak langsung.

2.4.6 General Nested Spatial Model (GNS)

General Nested Spatial Model digunakan untuk menganalisis seluruh bentuk ketergantungan spasial secara simultan, baik yang terjadi melalui variabel dependen, variabel independen, maupun komponen error. Model SDM dinyatakan sebagai berikut:

\[ y = \rho {W}{y} + {X}{\beta} + WX {\theta} + {\varepsilon}. \]

dengan komponen error spasial:

\[ u = \lambda Wu + {\varepsilon}. \] Keterangan:

\(y\) = vektor variabel dependen (IPM)
\(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin)
\(W\) = matriks pembobot spasial
\(\rho\) = koefisien autoregresif spasial yang menunjukkan pengaruh variabel dependen di wilayah sekitar terhadap wilayah ke-\(i\).
\(\beta\) = vektor koefisien regresi
\(WX\) = spatial lag dari variabel independen
\(\theta\) = vektor koefisien spatial lag variabel independen
\(u\) = komponen error yang mengandung ketergantungan spasial
\(\lambda\) = koefisien error spasial yang menunjukkan tingkat ketergantungan antar error wilayah
\(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika koefisien spasial \(\rho\), koefisien spatial lag variabel independen (\(\theta\)), dan koefisien error spasial \(\lambda\) signifikan, maka terdapat ketergantungan spasial yang bersifat komprehensif, yaitu melalui variabel dependen, variabel independen, dan komponen error secara simultan. Signifikansi \(\rho\) menunjukkan adanya interaksi langsung antarwilayah, di mana nilai variabel dependen pada suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Signifikansi \(\theta\) mengindikasikan adanya spillover effect dari variabel independen, yaitu perubahan variabel penjelas di suatu wilayah dapat memengaruhi variabel dependen di wilayah-wilayah tetangga. Sementara itu, signifikansi \(\lambda\) menandakan bahwa terdapat faktor-faktor eksternal yang tidak dimasukkan dalam model tetapi memiliki pola spasial, sehingga menyebabkan error antarwilayah saling berkaitan.

BAB III (METODOLOGI PENELITIAN)

3.1 Data dan Unit Spasial

Jenis data yang digunakan dalam penelitian ini merupakan data sekunder yang terdiri atas Indeks Pembangunan Manusia (IPM) sebagai variabel dependen serta Produk Domestik Regional Bruto (PDRB), Rata-Rata Lama Sekolah (RLS), Harapan Lama Sekolah (HLS), dan jumlah penduduk miskin sebagai variabel independen. Data tersebut dikumpulkan pada tingkat kabupaten/kota, sehingga memungkinkan analisis pembangunan manusia pada skala wilayah yang lebih rinci. Kelima variabel tersebut disusun dalam bentuk data cross section pada satu periode pengamatan yang sama, sehingga setiap kabupaten/kota diperlakukan sebagai satu unit observasi. Struktur data ini memungkinkan dilakukannya analisis hubungan spasial antarwilayah dengan mempertimbangkan kedekatan geografis antar kabupaten/kota.

Unit spasial dalam penelitian ini adalah wilayah administrasi kabupaten/kota. Setiap kabupaten/kota diperlakukan sebagai satu unit analisis yang merepresentasikan kondisi pembangunan manusia dan karakteristik sosial ekonomi wilayah tersebut. Pemilihan unit spasial pada tingkat kabupaten/kota didasarkan pada ketersediaan dan konsistensi data statistik yang bersumber dari publikasi resmi, serta relevansinya untuk menggambarkan variasi pembangunan manusia antarwilayah secara lebih rinci. Dengan menggunakan unit spasial kabupaten/kota, penelitian ini memungkinkan analisis pola persebaran dan ketergantungan spasial IPM serta faktor-faktor sosial ekonomi antarwilayah yang memiliki kedekatan geografis.

Pengolahan dan analisis data dilakukan menggunakan perangkat lunak R Studio, dengan memanfaatkan paket-paket pendukung analisis spasial dan ekonometrika spasial. Penggunaan R Studio memungkinkan integrasi antara pengolahan data, analisis statistik, pemodelan ekonometrika spasial, serta visualisasi peta tematik secara sistematis dan reprodusibel. Melalui pendekatan ini, diharapkan hasil analisis mampu memberikan gambaran yang komprehensif mengenai dinamika spasial pembangunan manusia antar kabupaten/kota.

3.2 Variabel Penelitian

Penelitian ini melibatkan satu variabel dependen dan empat variabel independen. Definisi operasional variabel dijelaskan sebagai berikut:

Jenis Variabel	Nama Variabel	Simbol	Satuan	Keterangan
Dependen	Indeks Pembangunan Manusia	\(Y\)	Indeks (0-100)	Indikator komposit yang mencerminkan capaian pembangunan manusia melalui dimensi kesehatan, pendidikan, dan standar hidup layak pada tingkat kabupaten/kota.
Independen	Produk Domestik Regional Bruto	\(X_1\)	Rupiah (Ribu)	Nilai tambah bruto seluruh barang dan jasa yang dihasilkan oleh berbagai sektor ekonomi dalam suatu wilayah pada periode tertentu.
Independen	Rata-Rata Lama Sekolah	\(X_2\)	Tahun	Rata-rata jumlah tahun pendidikan formal yang telah ditempuh oleh penduduk berusia 15 tahun ke atas.
Independen	Harapan Lama Sekolah	\(X_3\)	Tahun	Perkiraan jumlah tahun pendidikan formal yang diharapkan dapat ditempuh oleh anak usia sekolah di masa mendatang.
Independen	Jumlah Penduduk Miskin	\(X_4\)	Jiwa (orang)	Jumlah penduduk yang memiliki pengeluaran per kapita per bulan di bawah garis kemiskinan.

Fokus utama penelitian ini adalah untuk melihat pola keterkaitan spasial (spatial dependence) dan autokorelasi spasial (spatial autocorrelation) antarkabupaten/kota dalam hal Indeks Pembangunan Manusia, baik secara global maupun lokal. Serta model ekonometrika spasial yang paling tepat menggambarkan hubungan tersebut.

3.3 Metode Analisis Data

Analisis dalam penelitian ini dilakukan menggunakan pendekatan spasial dengan data cross-section spasial, yaitu data yang dikumpulkan pada satu periode waktu tertentu dengan unit observasi berupa kabupaten/kota. Pendekatan ini digunakan untuk mengkaji Indeks Pembangunan Manusia (IPM) antar kabupaten/kota dengan mempertimbangkan adanya keterkaitan dan interaksi spasial antarwilayah. Analisis dilakukan melalui beberapa tahapan, yaitu persiapan data (data preparation), analisis deskriptif, pengujian autokorelasi spasial, estimasi model ekonometrika spasial dan pemilihan model terbaik. Selain itu, sebagai analisis tambahan, penelitian ini juga menerapkan Geographically Weighted Regression (GWR) dan pemodelan interpolasi spasial.

3.3.1 Persiapan Data (Data Preparation)

Tahap persiapan data (data preparation) dalam penelitian ini dilakukan dengan mengintegrasikan data statistik dan data spasial kabupaten/kota di Provinsi Jawa Timur secara sistematis. Proses diawali dengan pemilihan wilayah studi, yaitu memfilter data spasial administrasi kabupaten/kota Indonesia sehingga hanya mencakup wilayah Provinsi Jawa Timur. Setiap poligon kabupaten/kota kemudian diberi identitas unik dan dikonversi ke dalam format simple features (sf), serta dilakukan pemeriksaan dan perbaikan validitas geometri untuk memastikan bahwa seluruh objek spasial layak digunakan dalam analisis spasial.

Selanjutnya dilakukan proses standarisasi nama wilayah pada data statistik dan data spasial menggunakan fungsi pembersihan teks, yang mencakup konversi huruf ke format huruf kecil, penghapusan tanda baca, serta penyelarasan istilah administratif seperti “kota” dan “kabupaten”. Langkah ini bertujuan untuk menghindari ketidaksesuaian penamaan wilayah yang dapat menghambat proses penggabungan data. Setelah itu, dilakukan pemeriksaan kesesuaian nama kabupaten/kota melalui proses pencocokan (matching) untuk mengidentifikasi wilayah yang tidak memiliki pasangan antara data statistik dan data spasial.

Dilanjutkan dengan penggabungan data statistik dengan data spasial menggunakan kunci nama wilayah yang telah distandarisasi. Proses ini menghasilkan satu set data spasial terintegrasi yang memuat informasi geografis serta variabel penelitian pada tingkat kabupaten/kota. Untuk memastikan kelengkapan data, dilakukan pengecekan terhadap keberadaan nilai kosong (missing values) pada variabel dependen, khususnya Indeks Pembangunan Manusia (IPM). Dengan demikian, tahap persiapan data ini memastikan bahwa data yang digunakan dalam analisis selanjutnya telah bersih, konsisten, dan siap digunakan untuk analisis deskriptif, pengujian autokorelasi spasial, serta estimasi model ekonometrika spasial.

3.3.2 Eksplorasi Data dan Visualisasi

Pada tahap ini, dilakukan perhitungan statistik deskriptif terhadap seluruh variabel penelitian, yaitu Indeks Pembangunan Manusia (IPM) sebagai variabel dependen serta Produk Domestik Regional Bruto (PDRB), Rata-Rata Lama Sekolah (RLS), Harapan Lama Sekolah (HLS), dan jumlah penduduk miskin sebagai variabel independen. Statistik deskriptif yang dihasilkan meliputi nilai minimum, maksimum, median, rata-rata, dan kuartil, yang bertujuan untuk menggambarkan karakteristik data, tingkat variasi antar kabupaten/kota, serta potensi ketimpangan nilai antarwilayah.

Selain analisis numerik, dilakukan pula visualisasi spasial terhadap variabel dependen, yaitu IPM, dalam bentuk peta tematik berbasis poligon administrasi kabupaten/kota di Provinsi Jawa Timur. Pemetaan ini bertujuan untuk memberikan gambaran awal mengenai pola persebaran IPM secara spasial serta mengidentifikasi indikasi awal adanya pengelompokan wilayah dengan capaian IPM yang relatif tinggi maupun rendah. Visualisasi spasial ini menjadi langkah awal untuk memahami dinamika pembangunan manusia antarwilayah sebelum dilakukan pengujian autokorelasi spasial dan estimasi model ekonometrika spasial. Dengan demikian, tahap analisis deskriptif ini memberikan dasar pemahaman yang kuat mengenai karakteristik dan pola awal data, yang selanjutnya menjadi pijakan bagi analisis spasial lanjutan dalam penelitian ini.

3.3.2 Uji Autokorelasi Spasial

Setelah dilakukan pemetaan untuk melihat pola sebaran IPM antarkabupaten/kota, langkah berikutnya adalah melakukan uji autokorelasi spasial. Tujuan pengujian ini adalah untuk memastikan apakah IPM di suatu kabupaten/kota memiliki hubungan atau ketergantungan dengan nilai IPM di kabupaten/kota yang berdekatan. Dengan kata lain, apakah Indeks Pembangunan Manusia di suatu kabupaten/kota dipengaruhi oleh kondisi kabupaten/kota sekitarnya. Dalam penelitian ini digunakan empat jenis uji autokorelasi spasial, yaitu Moran’s I, Geary’s C, Getis-Ord G, dan LISA.

Hasil keemapat uji autokorelasi spasial ini menjadi dasar penting dalam menentukan pendekatan analisis selanjutnya. Apabila ditemukan autokorelasi spasial yang signifikan, maka asumsi independensi antarobservasi tidak terpenuhi, sehingga diperlukan penggunaan model ekonometrika spasial untuk memperoleh estimasi yang lebih akurat dan inferensi statistik yang valid dalam menganalisis determinasi IPM antarkabupaten/kota.

3.3.3 Estimasi Model Ekonometrika Spatial

Setelah hasil uji Moran’s I menunjukkan adanya autokorelasi spasial pada Indeks Pembangunan Manusia (IPM) antar kabupaten/kota, langkah selanjutnya adalah menentukan bentuk ketergantungan spasial yang paling sesuai untuk memodelkan fenomena tersebut. Dalam analisis ekonometrika spasial, ketergantungan antarwilayah dapat muncul melalui berbagai mekanisme, baik secara langsung melalui variabel dependen, melalui variabel independen, maupun melalui faktor-faktor tidak teramati yang tercermin dalam komponen error. Oleh karena itu, penelitian ini tidak hanya membatasi analisis pada model spasial dasar, tetapi juga mempertimbangkan berbagai spesifikasi model yang mampu menangkap bentuk dependensi spasial secara lebih komprehensif.

Model Spatial Autoregressive Model (SAR) digunakan ketika ketergantungan spasial muncul secara langsung melalui variabel dependen, di mana nilai IPM suatu kabupaten/kota dipengaruhi oleh nilai IPM kabupaten/kota yang berdekatan. Sementara itu, Spatial Error Model (SEM) digunakan apabila ketergantungan spasial tidak terjadi secara langsung pada variabel dependen, melainkan melalui komponen error, yang mencerminkan adanya faktor-faktor tidak terobservasi yang bersifat spasial dan memengaruhi IPM secara tidak langsung antarwilayah. Selanjutnya, analisis diperluas dengan menggunakan Spatial Durbin Model (SDM) dan Spatial Durbin Error Model (SDEM). Model SDM memungkinkan ketergantungan spasial terjadi baik melalui variabel dependen maupun variabel independen, sehingga mampu menangkap pengaruh langsung dan tidak langsung (spillover effects) dari faktor-faktor penjelas antarwilayah. Sementara itu, model SDEM digunakan ketika pengaruh spasial terutama berasal dari variabel independen dan faktor tidak teramati yang tercermin dalam error spasial, tanpa adanya ketergantungan langsung pada variabel dependen. Untuk mengakomodasi kemungkinan adanya ketergantungan spasial yang lebih kompleks, penelitian ini juga mempertimbangkan Spatial Autoregressive Combined Model (SAC) dan General Nested Spatial Model (GNS). Model SAC digunakan ketika ketergantungan spasial terjadi secara simultan melalui variabel dependen dan komponen error. Adapun model GNS merupakan spesifikasi paling umum yang mencakup seluruh bentuk ketergantungan spasial, baik melalui variabel dependen, variabel independen, maupun komponen error, sehingga model-model lainnya dapat dipandang sebagai bentuk khusus dari model ini.

Perbandingan antar model tersebut dilakukan untuk mengidentifikasi spesifikasi model yang paling sesuai dengan kondisi empiris pembangunan manusia di Provinsi Jawa Timur.

3.3.4 Pemilihan Model Terbaik

Pemilihan model terbaik dilakukan untuk menentukan spesifikasi model ekonometrika spasial yang paling sesuai dalam menggambarkan hubungan antarkabupaten/kota pada data Indeks Pembangunan Manusia (IPM) di Provinsi Jawa Timur tahun 2024. Perbandingan antar model dilakukan dengan mempertimbangkan berbagai spesifikasi model spasial, yaitu Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS).

Evaluasi kinerja model dilakukan dengan menggunakan Akaike Information Criterion (AIC) dan Bayesian Information Criterion (BIC), di mana model dengan nilai AIC dan BIC yang lebih rendah menunjukkan tingkat kecocokan model yang lebih baik dengan tetap mempertimbangkan kompleksitas model. Selain itu, hasil Lagrange Multiplier (LM) test digunakan sebagai dasar untuk mengidentifikasi bentuk ketergantungan spasial yang dominan, baik yang muncul melalui variabel dependen, variabel independen, maupun komponen error.

Dengan mengombinasikan hasil LM test serta perbandingan nilai AIC dan BIC, model terbaik dipilih secara sistematis untuk memastikan spesifikasi model yang paling mampu merepresentasikan mekanisme ketergantungan spasial dalam pembentukan IPM antar kabupaten/kota. Pendekatan ini memungkinkan identifikasi sumber utama pengaruh spasial, apakah terjadi melalui interaksi langsung antarwilayah, spillover dari variabel penjelas, atau faktor-faktor eksternal yang tidak terobservasi, sehingga menghasilkan estimasi yang lebih akurat dan interpretasi yang lebih bermakna.

3.3.5 Model Geographically Weighted Regression (GWR)

Sebagai analisis tambahan, penelitian ini juga menerapkan Geographically Weighted Regression (GWR) untuk mengeksplorasi heterogenitas spasial, yaitu variasi pengaruh variabel-variabel penjelas terhadap Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur. Geographically Weighted Regression (GWR) dirumuskan sebagai berikut:

\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) x_{ik} + \varepsilon_i \] Keterangan:

\(y_i\) adalah nilai variabel dependen pada lokasi ke-\(i\).
\(\beta_0(u_i, v_i)\) adalah intersep regresi pada lokasi geografis dengan koordinat \((u_i, v_i)\).
\(\beta_k(u_i, v_i)\) adalah parameter regresi untuk variabel ke-\(k\) pada lokasi geografis \((u_i, v_i)\).
\(x_{ik}\) adalah variabel independen ke-\(k\) pada lokasi ke-\(i\).
\(\varepsilon_i\) adalah error residual pada lokasi ke-\(i\).

Analisis GWR ini dilakukan dengan penyiapan titik lokasi spasial yang diperoleh dari titik pusat (centroid) masing-masing kabupaten/kota. Data spasial tersebut kemudian dikonversi ke dalam format titik dan diproyeksikan ke sistem koordinat UTM guna memastikan perhitungan jarak antarwilayah dilakukan secara akurat. Selanjutnya, dilakukan penentuan bandwidth optimal menggunakan metode adaptive bandwidth, yang memungkinkan jumlah tetangga yang dipertimbangkan dalam estimasi bervariasi antar lokasi. Pemilihan bandwidth ini bertujuan untuk memperoleh keseimbangan antara tingkat kedekatan spasial dan kestabilan estimasi parameter. Setelah bandwidth optimal diperoleh, model GWR diestimasi dengan memasukkan seluruh variabel penjelas, sehingga dihasilkan koefisien regresi yang bersifat lokal untuk setiap kabupaten/kota.

Kinerja model GWR kemudian dibandingkan dengan model regresi global melalui perbandingan nilai Akaike Information Criterion corrected (AICc) pada GWR dan nilai AIC pada model OLS. Perbandingan ini digunakan untuk menilai apakah pendekatan lokal GWR mampu memberikan peningkatan kecocokan model dibandingkan pendekatan global. Selanjutnya, hasil estimasi GWR diekstraksi untuk memperoleh nilai koefisien lokal masing-masing variabel serta residual model pada setiap kabupaten/kota. Selanjutnya, dilakukan visualisasi spasial terhadap koefisien regresi lokal dan residual GWR dalam bentuk peta tematik. Pemetaan ini bertujuan untuk mengidentifikasi variasi geografis pengaruh masing-masing variabel terhadap IPM serta untuk mengevaluasi pola sebaran residual model. Dengan demikian, analisis GWR memberikan pemahaman tambahan mengenai perbedaan pengaruh faktor-faktor penentu IPM antarwilayah, yang tidak sepenuhnya dapat ditangkap oleh model ekonometrika spasial global.

3.3.5 Interpolasi Spasial

Sebagai analisis tambahan, penelitian ini juga melakukan pemodelan interpolasi spasial untuk memperoleh gambaran permukaan spasial Indeks Pembangunan Manusia (IPM) yang bersifat kontinu di Provinsi Jawa Timur. Interpolasi dilakukan dengan menggunakan titik observasi yang direpresentasikan oleh titik pusat (centroid) masing-masing kabupaten/kota. Data titik tersebut kemudian diproyeksikan ke dalam sistem koordinat UTM guna memastikan perhitungan jarak antar lokasi dilakukan secara akurat. Selanjutnya, dibentuk grid prediksi yang mencakup seluruh wilayah Provinsi Jawa Timur sebagai dasar untuk estimasi nilai IPM pada lokasi-lokasi yang tidak memiliki pengamatan langsung.

Metode interpolasi pertama yang digunakan adalah Inverse Distance Weighting (IDW), di mana estimasi nilai IPM pada suatu lokasi ditentukan berdasarkan rata-rata berbobot dari nilai IPM di lokasi sekitarnya, dengan bobot yang berbanding terbalik terhadap jarak. Parameter pangkat jarak (power parameter) dan jumlah tetangga terdekat digunakan untuk mengontrol pengaruh kedekatan spasial dalam proses interpolasi. Metode ini memberikan pendekatan deterministik yang sederhana untuk menggambarkan variasi spasial IPM.

Metode interpolasi kedua adalah Thin Plate Splines (TPS), yang merupakan pendekatan berbasis fungsi halus (smoothing function) untuk memodelkan permukaan spasial secara kontinu. TPS memungkinkan pembentukan permukaan IPM yang lebih mulus dengan meminimalkan kelengkungan permukaan, sehingga cocok untuk menggambarkan tren spasial global tanpa terlalu dipengaruhi oleh fluktuasi lokal yang ekstrem.

Selanjutnya, dilakukan interpolasi menggunakan Ordinary Kriging, yaitu metode geostatistik yang mempertimbangkan struktur ketergantungan spasial melalui pemodelan variogram. Pada tahap ini, dilakukan estimasi variogram empiris dan penyesuaian model variogram teoretis untuk menggambarkan pola variabilitas spasial IPM. Model variogram yang telah diperoleh digunakan dalam proses kriging untuk menghasilkan estimasi IPM pada grid prediksi. Selain itu, dilakukan validasi silang (cross-validation) dengan pendekatan leave-one-out untuk mengevaluasi kinerja model kriging, yang ditunjukkan melalui nilai Root Mean Square Error (RMSE).

Melalui penerapan berbagai metode interpolasi spasial tersebut, penelitian ini bertujuan untuk melengkapi hasil analisis utama dengan memberikan visualisasi dan pemahaman tambahan mengenai variasi spasial IPM di Provinsi Jawa Timur. Analisis interpolasi ini tidak digunakan sebagai dasar penarikan kesimpulan utama, tetapi berfungsi sebagai alat eksploratif untuk memperkaya interpretasi pola spasial pembangunan manusia antarwilayah.

BAB IV (HASIL DAN PEMBAHASAN)

4.1 Persiapan Data (Data Preparation)

# Filter for East Java province
JaTim_sp <- IDN_kab[IDN_kab$NAME_1 == "Jawa Timur", ]
cat("GADM East Java polygons:", nrow(JaTim_sp), "\n")

## GADM East Java polygons: 38

JaTim_sp$id <- seq_len(nrow(JaTim_sp))
JaTim_sf <- st_as_sf(JaTim_sp) |> st_make_valid()

# Enhanced cleaning function that preserves more distinctions
clean_nm <- function(x){
  x |>
    str_to_lower() |>
    str_trim() |>
    # Remove "Kabupaten" or "Kab." prefix but keep "Kota"
    str_replace_all("^kabupaten\\s+", "") |>
    str_replace_all("^kab\\.\\s+", "") |>
    str_replace_all("^kab\\s+", "") |>
    # Standardize "Kota" prefix
    str_replace_all("^kota\\s+", "kota ") |>
    # Clean up punctuation but preserve spaces
    str_replace_all("\\.", "") |>
    str_replace_all(",", "") |>
    str_squish()
}

# Apply cleaning
jatim_dat <- jatim_dat |>
  mutate(
    key_nm = clean_nm(`Kab/Kota`),
    original_name = `Kab/Kota`
  )

JaTim_sf <- JaTim_sf |>
  mutate(
    key_nm = clean_nm(NAME_2),
    original_name = NAME_2
  )

# Check matching status
cat("\n=== MATCHING CHECK ===\n")

## 
## === MATCHING CHECK ===

cat("Unique regions in CSV:", length(unique(jatim_dat$key_nm)), "\n")

## Unique regions in CSV: 38

cat("Unique regions in shapefile:", length(unique(JaTim_sf$key_nm)), "\n")

## Unique regions in shapefile: 38

# Show all names for comparison
cat("\n--- CSV Names (cleaned) ---\n")

## 
## --- CSV Names (cleaned) ---

csv_names <- jatim_dat |>
  select(original_name, key_nm) |>
  arrange(key_nm) |>
  tibble::as_tibble()
print(csv_names, n = Inf)

## # A tibble: 38 × 2
##    original_name         key_nm          
##    <chr>                 <chr>           
##  1 Kabupaten Bangkalan   bangkalan       
##  2 Kabupaten Banyuwangi  banyuwangi      
##  3 Batu                  batu            
##  4 Kabupaten Blitar      blitar          
##  5 Kabupaten Bojonegoro  bojonegoro      
##  6 Kabupaten Bondowoso   bondowoso       
##  7 Kabupaten Gresik      gresik          
##  8 Kabupaten Jember      jember          
##  9 Kabupaten Jombang     jombang         
## 10 Kabupaten Kediri      kediri          
## 11 Kota Blitar           kota blitar     
## 12 Kota Kediri           kota kediri     
## 13 Kota Madiun           kota madiun     
## 14 Kota Malang           kota malang     
## 15 Kota Mojokerto        kota mojokerto  
## 16 Kota Pasuruan         kota pasuruan   
## 17 Kota Probolinggo      kota probolinggo
## 18 Kabupaten Lamongan    lamongan        
## 19 Kabupaten Lumajang    lumajang        
## 20 Kabupaten Madiun      madiun          
## 21 Kabupaten Magetan     magetan         
## 22 Kabupaten Malang      malang          
## 23 Kabupaten Mojokerto   mojokerto       
## 24 Kabupaten Nganjuk     nganjuk         
## 25 Kabupaten Ngawi       ngawi           
## 26 Kabupaten Pacitan     pacitan         
## 27 Kabupaten Pamekasan   pamekasan       
## 28 Kabupaten Pasuruan    pasuruan        
## 29 Kabupaten Ponorogo    ponorogo        
## 30 Kabupaten Probolinggo probolinggo     
## 31 Kabupaten Sampang     sampang         
## 32 Kabupaten Sidoarjo    sidoarjo        
## 33 Kabupaten Situbondo   situbondo       
## 34 Kabupaten Sumenep     sumenep         
## 35 Surabaya              surabaya        
## 36 Kabupaten Trenggalek  trenggalek      
## 37 Kabupaten Tuban       tuban           
## 38 Kabupaten Tulungagung tulungagung

cat("\n--- Shapefile Names (cleaned) ---\n")

## 
## --- Shapefile Names (cleaned) ---

shp_names <- JaTim_sf |> 
  st_drop_geometry() |>
  select(original_name, key_nm) |> 
  arrange(key_nm) |>
  tibble::as_tibble()
print(shp_names, n = Inf)

## # A tibble: 38 × 2
##    original_name    key_nm          
##    <chr>            <chr>           
##  1 Bangkalan        bangkalan       
##  2 Banyuwangi       banyuwangi      
##  3 Batu             batu            
##  4 Blitar           blitar          
##  5 Bojonegoro       bojonegoro      
##  6 Bondowoso        bondowoso       
##  7 Gresik           gresik          
##  8 Jember           jember          
##  9 Jombang          jombang         
## 10 Kediri           kediri          
## 11 Kota Blitar      kota blitar     
## 12 Kota Kediri      kota kediri     
## 13 Kota Madiun      kota madiun     
## 14 Kota Malang      kota malang     
## 15 Kota Mojokerto   kota mojokerto  
## 16 Kota Pasuruan    kota pasuruan   
## 17 Kota Probolinggo kota probolinggo
## 18 Lamongan         lamongan        
## 19 Lumajang         lumajang        
## 20 Madiun           madiun          
## 21 Magetan          magetan         
## 22 Malang           malang          
## 23 Mojokerto        mojokerto       
## 24 Nganjuk          nganjuk         
## 25 Ngawi            ngawi           
## 26 Pacitan          pacitan         
## 27 Pamekasan        pamekasan       
## 28 Pasuruan         pasuruan        
## 29 Ponorogo         ponorogo        
## 30 Probolinggo      probolinggo     
## 31 Sampang          sampang         
## 32 Sidoarjo         sidoarjo        
## 33 Situbondo        situbondo       
## 34 Sumenep          sumenep         
## 35 Surabaya         surabaya        
## 36 Trenggalek       trenggalek      
## 37 Tuban            tuban           
## 38 Tulungagung      tulungagung

# Show unmatched regions from CSV
cek_unmatch_csv <- anti_join(jatim_dat, st_drop_geometry(JaTim_sf), by = "key_nm")
if(nrow(cek_unmatch_csv) > 0) {
  cat("\n⚠️ UNMATCHED in CSV (", nrow(cek_unmatch_csv), "regions):\n")
  print(cek_unmatch_csv[, c("original_name", "key_nm")])
}

# Show unmatched regions from shapefile
cek_unmatch_shp <- anti_join(st_drop_geometry(JaTim_sf), jatim_dat, by = "key_nm")
if(nrow(cek_unmatch_shp) > 0) {
  cat("\n⚠️ UNMATCHED in Shapefile (", nrow(cek_unmatch_shp), "regions):\n")
  print(cek_unmatch_shp[, c("original_name", "key_nm")])
}

# Manual matching dictionary for problematic cases
manual_fixes <- data.frame(
  csv_name = c(
    # Add manual mappings here if needed after seeing the output
    # Example: "special case in CSV"
  ),
  shp_name = c(
    # Corresponding shapefile names
    # Example: "special case in shapefile"
  ),
  stringsAsFactors = FALSE
)

# Apply manual fixes if any exist
if(nrow(manual_fixes) > 0) {
  # Create lookup for CSV
  csv_lookup <- setNames(manual_fixes$shp_name, manual_fixes$csv_name)
  
  jatim_dat <- jatim_dat |>
    mutate(key_nm_fixed = ifelse(
      key_nm %in% names(csv_lookup),
      csv_lookup[key_nm],
      key_nm
    ))
  
  # Use fixed names for merging
  JaTim_merged <- JaTim_sf |>
    inner_join(jatim_dat, by = c("key_nm" = "key_nm_fixed"))
} else {
  # Standard merge without manual fixes
  JaTim_merged <- JaTim_sf |>
    inner_join(jatim_dat, by = "key_nm")
}

cat("\n=== AFTER MERGE ===\n")

## 
## === AFTER MERGE ===

cat("Final observations:", nrow(JaTim_merged), "\n")

## Final observations: 38

cat("Missing Y values:", sum(is.na(JaTim_merged$Y)), "\n")

## Missing Y values: 0

# Verify we have 38 observations
if(nrow(JaTim_merged) != 38) {
  cat("\n!!! WARNING: Expected 38 observations but got", nrow(JaTim_merged), "!!!\n")
  cat("\nDIAGNOSTIC INFO:\n")
  cat("- CSV has", nrow(jatim_dat), "rows\n")
  cat("- Shapefile has", nrow(JaTim_sf), "polygons\n")
  cat("- Matched regions:", nrow(JaTim_merged), "\n")
  cat("- Unmatched CSV:", nrow(cek_unmatch_csv), "\n")
  cat("- Unmatched Shapefile:", nrow(cek_unmatch_shp), "\n")
  cat("\nPlease check the name lists above and add manual fixes if needed.\n\n")
  
  # Stop execution if count is wrong
  stop("Observation count mismatch! Expected 38, got ", nrow(JaTim_merged))
} else {
  cat("\n✅ SUCCESS: Got exactly 38 observations as expected!\n\n")
}

## 
## ✅ SUCCESS: Got exactly 38 observations as expected!

Interpretasi:

Output tersebut menunjukkan bahwa proses penggabungan data berjalan dengan sangat baik dan konsisten. Shapefile GADM Provinsi Jawa Timur memiliki 38 polygon wilayah, yang merepresentasikan seluruh kabupaten dan kota di Jawa Timur, dan jumlah ini tepat sama dengan jumlah wilayah unik yang terdapat dalam data CSV. Setelah dilakukan proses pembersihan nama wilayah (cleaning), tidak ditemukan perbedaan antara nama wilayah di CSV dan di shapefile, sehingga seluruh wilayah dapat dipasangkan secara satu-ke-satu tanpa konflik atau ketidaksesuaian. Hasil penggabungan (merge) menghasilkan 38 observasi akhir, sesuai dengan jumlah wilayah administratif, dan tidak terdapat nilai variabel Y yang hilang. Dengan demikian, integritas data spasial dan nonspasial terjaga dengan baik, proses matching berhasil 100%, dan dataset tersebut siap digunakan untuk analisis lanjutan.

4.2 Eksplorasi Data dan Visualisasi

# Basic Statistics
summary(
  JaTim_merged |>
    dplyr::select(Y, X1, X2, X3, X4) |>
    sf::st_drop_geometry()
)

##        Y               X1               X2               X3       
##  Min.   :66.72   Min.   : 24484   Min.   : 5.080   Min.   :11.98  
##  1st Qu.:71.69   1st Qu.: 34058   1st Qu.: 7.478   1st Qu.:12.91  
##  Median :74.56   Median : 39509   Median : 8.060   Median :13.49  
##  Mean   :75.31   Mean   : 75005   Mean   : 8.464   Mean   :13.58  
##  3rd Qu.:78.62   3rd Qu.: 75381   3rd Qu.: 9.832   3rd Qu.:14.02  
##  Max.   :84.69   Max.   :565841   Max.   :12.110   Max.   :15.79  
##        X4        
##  Min.   :  6.59  
##  1st Qu.: 68.07  
##  Median :107.49  
##  Mean   :104.81  
##  3rd Qu.:146.44  
##  Max.   :240.14

Interpretasi:

Hasil statistik deskriptif menunjukkan bahwa rata-rata Y yang menunjukan Indeks Pembangunan Manusia (IPM) kabupaten/kota di Provinsi Jawa Timur sebesar 75,31, dengan nilai minimum 66,72 dan maksimum 84,69. Rentang nilai ini menunjukkan adanya kesenjangan pembangunan manusia antarwilayah. Dari sisi ekonomi, pada variabel X1 yang menunjukan PDRB memiliki variasi yang sangat besar antar kabupaten/kota, mencerminkan ketimpangan kapasitas ekonomi daerah. Pada dimensi pendidikan yang ditunjukan oleh X2 yaitu Rata-Rata Lama Sekolah (RLS) dan X3 yaitu Harapan Lama Sekolah (HLS) juga menunjukkan perbedaan yang cukup nyata antarwilayah, yang mengindikasikan ketimpangan akses dan capaian pendidikan. Sementara itu, X4 yang menunjukan jumlah penduduk miskin memperlihatkan variasi yang signifikan antar kabupaten/kota, menggambarkan ketidakmerataan tingkat kesejahteraan masyarakat. Secara keseluruhan, hasil ini mengindikasikan bahwa pembangunan manusia di Jawa Timur memiliki pola yang tidak merata dan berpotensi dipengaruhi oleh faktor-faktor yang bersifat spasial.

# Map of Variable Y
ggplot(JaTim_merged) +
  geom_sf(aes(fill = Y), color="white", size=0.2) +
  scale_fill_viridis_c(option="magma", na.value="grey85") +
  labs(title="East Java — Distribution of Y Variable (n=38)", fill="Y") +
  theme_minimal()

Interpretasi:

Peta persebaran Indeks Pembangunan Manusia (IPM) di Provinsi Jawa Timur memperlihatkan adanya indikasi awal pengelompokan spasial yang cukup jelas antar kabupaten/kota. Wilayah dengan nilai IPM tinggi hingga sangat tinggi (ditunjukkan oleh warna kuning muda hingga merah terang) cenderung terkonsentrasi di kawasan perkotaan dan pusat pertumbuhan, seperti Kota Surabaya dan beberapa kabupaten/kota di sekitarnya. Sebaliknya, wilayah dengan nilai IPM relatif lebih rendah (ditunjukkan oleh warna ungu tua hingga gelap) lebih banyak tersebar di wilayah perdesaan dan kawasan pinggiran, termasuk beberapa kabupaten di bagian timur dan selatan Jawa Timur. Pola persebaran ini menunjukkan bahwa capaian pembangunan manusia di Jawa Timur tidak tersebar secara merata, melainkan membentuk kecenderungan pengelompokan wilayah dengan karakteristik yang serupa.

4.2 Uji Autokorelasi Spasial

4.2.1 Uji Moran’s I

# Queen contiguity
nb_queen <- poly2nb(JaTim_merged, queen = TRUE)

## Warning in poly2nb(JaTim_merged, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

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

coords <- st_coordinates(st_centroid(JaTim_merged))

## Warning: st_centroid assumes attributes are constant over geometries

# Fix isolates with 1-NN
iso <- which(card(nb_queen) == 0)
if(length(iso) > 0){
  knn_nb <- knn2nb(knearneigh(coords, k = 1))
  for(i in iso) nb_queen[[i]] <- knn_nb[[i]]
}

## Warning in knn2nb(knearneigh(coords, k = 1)): neighbour object has 13
## sub-graphs

# Final weights (NO zero.policy needed anymore)
lwW <- nb2listw(nb_queen, style = "W", zero.policy = FALSE)
lwB <- nb2listw(nb_queen, style = "B", zero.policy = FALSE)

cat("Isolated:", sum(card(nb_queen) == 0), "\n")

## Isolated: 0

cat("Avg neighbors:", mean(card(nb_queen)), "\n")

## Avg neighbors: 2.447368

moran_res <- moran.test(JaTim_merged$Y, lwW,
                        randomisation = TRUE,
                        alternative = "two.sided")
print(moran_res)

## 
##  Moran I test under randomisation
## 
## data:  JaTim_merged$Y  
## weights: lwW    
## 
## Moran I statistic standard deviate = 2.663, p-value = 0.007745
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.36946012       -0.02702703        0.02216753

Interpretasi:

Hasil uji Moran’s I menunjukkan bahwa nilai statistik Moran’s I sebesar 0,3695 dengan p-value sebesar 0,0077, yang signifikan pada tingkat signifikansi 5 persen. Nilai Moran’s I bernilai positif dan signifikan, mengindikasikan adanya autokorelasi spasial positif pada Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur. Artinya, wilayah dengan nilai IPM tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki IPM tinggi, begitupun sebaliknya. Dengan demikian, persebaran IPM di Jawa Timur tidak terjadi secara acak, melainkan membentuk pola pengelompokan spasial (spatial clustering). Hal ini menunjukkan adanya ketergantungan spasial dalam capaian pembangunan manusia, di mana kondisi IPM suatu kabupaten/kota berpotensi dipengaruhi oleh kondisi wilayah sekitarnya. Oleh karena itu, hasil ini memperkuat pentingnya penggunaan pendekatan dan model ekonometrika spasial dalam analisis lanjutan agar estimasi yang diperoleh lebih tepat dan tidak bias akibat pengabaian efek spasial.

4.2.2 Uji Geary’s C

geary_res <- geary.test(JaTim_merged$Y, lwW,
                        randomisation = TRUE,
                        alternative = "two.sided")
print(geary_res)

## 
##  Geary C test under randomisation
## 
## data:  JaTim_merged$Y 
## weights: lwW   
## 
## Geary C statistic standard deviate = 3.1068, p-value = 0.001891
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.49857222        1.00000000        0.02604891

Interpretasi:

Hasil uji Geary’s C menunjukkan nilai sebesar 0,4986 dengan p-value sebesar 0,0019, yang signifikan pada tingkat signifikansi 5 persen. Nilai Geary’s C lebih kecil dari 1 dan signifikan, mengindikasikan adanya autokorelasi spasial positif pada Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur. Hal ini berarti bahwa wilayah-wilayah yang berdekatan secara geografis cenderung memiliki nilai IPM yang relatif mirip. Dengan demikian, hasil uji Geary’s C semakin menegaskan adanya ketergantungan spasial dalam capaian pembangunan manusia antar wilayah di Jawa Timur.

4.2.3 Getis Ord

Gi_z <- localG(JaTim_merged$Y, lwW)

JaTim_G <- JaTim_merged |>
  dplyr::mutate(
    z_Gistar = as.numeric(Gi_z),
    hotcold = dplyr::case_when(
      z_Gistar >=  1.96 ~ "Hot spot (p≈0.05)",
      z_Gistar <= -1.96 ~ "Cold spot (p≈0.05)",
      TRUE ~ "Not significant"
    )
  )
ggplot(JaTim_G) +
  geom_sf(aes(fill = hotcold), color = "white", linewidth = 0.2) +
  scale_fill_manual(values = c(
    "Hot spot (p≈0.05)" = "#b2182b",
    "Cold spot (p≈0.05)" = "#2166ac",
    "Not significant" = "grey85"
  )) +
  labs(title = "East Java — Getis-Ord Gi* Hot/Cold Spots", fill = NULL) +
  theme_minimal()

Interpretasi:

Peta Getis-Ord Gi* memperlihatkan bahwa hanya sebagian kecil wilayah di Provinsi Jawa Timur yang menunjukkan pola hot spot atau cold spot IPM yang signifikan secara statistik pada tingkat signifikansi 5 persen. Cold spot IPM (warna biru) terlihat terkonsentrasi pada beberapa kabupaten yang saling berdekatan di wilayah timur laut Jawa Timur, yang menunjukkan adanya klaster wilayah dengan nilai IPM relatif rendah yang dikelilingi oleh wilayah dengan karakteristik serupa. Kondisi ini mengindikasikan bahwa rendahnya capaian pembangunan manusia di wilayah tersebut bersifat terlokalisasi dan membentuk pengelompokan spasial yang signifikan. Sementara itu, hot spot IPM (warna merah) hanya muncul secara terbatas, yaitu pada satu wilayah yang teridentifikasi sebagai pusat dengan konsentrasi IPM tinggi, yang dikelilingi oleh wilayah dengan nilai IPM relatif tinggi pula. Di luar area-area tersebut, mayoritas kabupaten/kota di Jawa Timur ditampilkan dengan warna abu-abu, yang menunjukkan bahwa tidak terdapat pengelompokan IPM yang signifikan secara statistik pada tingkat lokal. Hal ini mengindikasikan bahwa meskipun secara global IPM di Jawa Timur menunjukkan adanya autokorelasi spasial, pola pengelompokan yang kuat secara lokal hanya terjadi di wilayah-wilayah tertentu saja.

4.2.4 Local Moran’s I (LISA)

x <- as.numeric(scale(JaTim_merged$Y))
lagx <- lag.listw(lwW, x)

lisa <- localmoran(x, lwW)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")

alpha <- 0.05
quad_raw <- 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)"
)

JaTim_LISA <- dplyr::bind_cols(JaTim_merged, lisa_df) |>
  dplyr::mutate(
    quad = ifelse(Pi.two.sided <= alpha, quad_raw, "Not significant")
  )
ggplot(JaTim_LISA) +
  geom_sf(aes(fill = quad), color = "white", linewidth = 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 = "East Java — Local Moran's I (LISA)", fill = "Category") +
  theme_minimal()

Interpretasi:

Peta Local Moran’s I (LISA) menunjukkan bahwa hanya sebagian kecil kabupaten/kota di Provinsi Jawa Timur yang memiliki autokorelasi spasial lokal yang signifikan pada tingkat signifikansi 5 persen. Klaster High–High, yang ditandai dengan warna merah, mengindikasikan wilayah dengan nilai IPM tinggi yang dikelilingi oleh wilayah lain dengan IPM relatif tinggi, sehingga mencerminkan adanya pengelompokan capaian pembangunan manusia yang tinggi pada area tersebut. Sebaliknya, klaster Low–Low, yang ditandai dengan warna biru, menunjukkan wilayah dengan nilai IPM rendah yang dikelilingi oleh wilayah dengan IPM rendah pula, sehingga membentuk klaster wilayah dengan capaian pembangunan manusia yang relatif tertinggal. Sementara itu, mayoritas wilayah lainnya ditunjukkan dengan warna abu-abu, yang berarti tidak memiliki autokorelasi spasial lokal yang signifikan. Hal ini menunjukkan bahwa meskipun secara global IPM di Jawa Timur memperlihatkan ketergantungan spasial, pengelompokan IPM yang signifikan secara lokal hanya terjadi pada wilayah-wilayah tertentu dan tidak merata di seluruh provinsi.

4.3 Estimasi Model Spasial

4.3.1 Ordinary Least Squares (OLS)

vars_use <- c("Y","X1","X2","X3","X4")

dat_jatim <- JaTim_merged |>
  dplyr::select(dplyr::all_of(vars_use)) |>
  sf::st_drop_geometry()

cc <- complete.cases(dat_jatim)
dat_jatim <- dat_jatim[cc, , drop = FALSE]

# SUBSET WEIGHTS YANG SAMA
lw_jatim <- spdep::subset.listw(lwW, subset = cc)

ols_jatim <- lm(Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
summary(ols_jatim)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5936 -0.4965 -0.1574  0.5287  2.1325 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.877e+01  2.876e+00  13.482 5.64e-15 ***
## X1          -9.190e-07  1.745e-06  -0.527 0.601909    
## X2           2.426e+00  1.597e-01  15.196  < 2e-16 ***
## X3           1.140e+00  2.589e-01   4.404 0.000106 ***
## X4           5.565e-03  2.936e-03   1.896 0.066813 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8039 on 33 degrees of freedom
## Multiple R-squared:  0.9745, Adjusted R-squared:  0.9715 
## F-statistic: 315.8 on 4 and 33 DF,  p-value: < 2.2e-16

lm.morantest(ols_jatim, lw_jatim)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## weights: lw_jatim
## 
## Moran I statistic standard deviate = -1.066, p-value = 0.8568
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.20384553      -0.04703029       0.02163859

lm.LMtests(ols_jatim, lw_jatim,
           test = c("LMlag","LMerr","RLMlag","RLMerr"))

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

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## test weights: listw
## 
## RSlag = 3.4621, df = 1, p-value = 0.06279
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## test weights: listw
## 
## RSerr = 1.7577, df = 1, p-value = 0.1849
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## test weights: listw
## 
## adjRSlag = 5.0003, df = 1, p-value = 0.02534
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
## test weights: listw
## 
## adjRSerr = 3.2959, df = 1, p-value = 0.06945

Interpretasi:

Hasil estimasi model OLS menunjukkan bahwa model memiliki daya jelas yang sangat tinggi, dengan nilai Adjusted R-squared sebesar 0,9715 atau sekitar 97 persen variasi IPM antar kabupaten/kota di Jawa Timur dapat dijelaskan oleh variabel independen dalam model. Secara parsial, X2 dan X3 berpengaruh positif dan signifikan terhadap IPM pada tingkat signifikansi 5 persen, sementara X4 signifikan pada tingkat 10 persen, dan X1 tidak berpengaruh signifikan. Uji F juga menunjukkan bahwa model secara keseluruhan signifikan. Uji Moran’s I terhadap residual model OLS menghasilkan p-value sebesar 0,8568, sehingga tidak terdapat bukti adanya autokorelasi spasial pada residual. Hal ini mengindikasikan bahwa secara umum model OLS telah mampu menangkap variasi spasial dalam data.

Namun demikian, hasil uji Lagrange Multiplier (LM) menunjukkan indikasi yang lebih spesifik. Uji LM-lag dan LM-error tidak signifikan pada tingkat 5 persen, tetapi uji Robust LM-lag (adjRSlag) signifikan pada tingkat 5 persen, sementara Robust LM-error tidak signifikan. Temuan ini mengindikasikan bahwa jika diperlukan pemodelan spasial lanjutan, model Spatial Lag (SAR) lebih tepat dibandingkan model Spatial Error. Dengan demikian, meskipun residual OLS tidak menunjukkan autokorelasi spasial yang kuat, terdapat indikasi ketergantungan spasial dalam bentuk interaksi antarwilayah pada variabel dependen, sehingga estimasi menggunakan model SAR layak dipertimbangkan sebagai pembanding atau model utama.

4.3.2 Spatial Autoregressive Model (SAR)

sar_jatim <- lagsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sar_jatim)

## 
## Call:
## lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, listw = lw_jatim, 
##     method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.30234 -0.31942 -0.23323  0.55119  1.95838 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  3.3113e+01  3.7202e+00  8.9009 < 2.2e-16
## X1          -1.2377e-06  1.5532e-06 -0.7968    0.4255
## X2           2.2687e+00  1.5956e-01 14.2183 < 2.2e-16
## X3           1.2321e+00  2.3305e-01  5.2867 1.245e-07
## X4           3.7309e-03  2.7590e-03  1.3523    0.1763
## 
## Rho: 0.080241, LR test value: 3.8945, p-value: 0.048445
## Asymptotic standard error: 0.038543
##     z-value: 2.0819, p-value: 0.037353
## Wald statistic: 4.3342, p-value: 0.037353
## 
## Log likelihood: -40.99531 for lag model
## ML residual variance (sigma squared): 0.50543, (sigma: 0.71094)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 95.991, (AIC for lm: 97.885)
## LM test for residual autocorrelation
## test value: 3.6403, p-value: 0.056395

Interpretasi:

Hasil estimasi Spatial Autoregressive Model (SAR) menunjukkan bahwa terdapat ketergantungan spasial yang signifikan pada IPM kabupaten/kota di Provinsi Jawa Timur. Hal ini ditunjukkan oleh nilai koefisien spasial ρ (rho) sebesar 0,0802 yang signifikan pada tingkat 5 persen, mengindikasikan bahwa IPM suatu wilayah dipengaruhi secara positif oleh IPM wilayah-wilayah tetangganya. Secara parsial, variabel X2 dan X3 berpengaruh positif dan signifikan terhadap IPM, sementara X1 dan X4 tidak signifikan. Model SAR menunjukkan kinerja yang lebih baik, ditunjukkan oleh nilai AIC yang lebih rendah yaitu 95,991 serta varians residual yang lebih kecil. Selain itu, uji autokorelasi residual menunjukkan tidak adanya autokorelasi spasial yang signifikan pada residual model, sehingga model SAR dapat dianggap telah menangkap ketergantungan spasial dalam data secara memadai.

4.3.3 Spatial Error Model (SEM)

sem_jatim <- errorsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sem_jatim)

## 
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, 
##     listw = lw_jatim, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.51516 -0.48260 -0.10926  0.49160  1.86113 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  3.8398e+01  2.5741e+00 14.9169 < 2.2e-16
## X1          -1.2935e-06  1.6086e-06 -0.8041   0.42133
## X2           2.4892e+00  1.2724e-01 19.5629 < 2.2e-16
## X3           1.1307e+00  2.3091e-01  4.8966  9.75e-07
## X4           5.5814e-03  2.3735e-03  2.3516   0.01869
## 
## Lambda: -0.2903, LR test value: 2.4034, p-value: 0.12107
## Asymptotic standard error: 0.1648
##     z-value: -1.7616, p-value: 0.078139
## Wald statistic: 3.1032, p-value: 0.078139
## 
## Log likelihood: -41.74083 for error model
## ML residual variance (sigma squared): 0.51247, (sigma: 0.71587)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 97.482, (AIC for lm: 97.885)

Interpretasi:

Hasil estimasi Spatial Error Model (SEM) menunjukkan bahwa komponen ketergantungan spasial pada error tidak signifikan. Hal ini ditunjukkan oleh nilai parameter λ (lambda) sebesar −0,2903 yang tidak signifikan pada tingkat 5 persen. Dengan demikian, tidak terdapat bukti kuat bahwa korelasi spasial pada IPM Jawa Timur berasal dari faktor tak teramati yang bersifat spasial. Secara parsial, variabel X2, X3, dan X4 berpengaruh positif dan signifikan terhadap IPM, sementara X1 tidak signifikan. Model SEM hanya memberikan perbaikan yang sangat terbatas, yang tercermin dari nilai AIC 97,482 yang hanya sedikit lebih rendah daripada OLS dan masih lebih tinggi dibandingkan model SAR. Oleh karena itu, hasil ini mengindikasikan bahwa model SEM kurang tepat untuk merepresentasikan struktur ketergantungan spasial dalam data.

4.3.4 Spatial Durbin Model (SDM)

sdm_jatim <- lagsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(sdm_jatim)

## 
## Call:
## lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, listw = lw_jatim, 
##     Durbin = TRUE, method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.437039 -0.421597 -0.017305  0.326250  1.650223 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  4.0274e+01  8.3837e+00  4.8038 1.557e-06
## X1          -2.0002e-06  1.5511e-06 -1.2895   0.19722
## X2           2.2056e+00  1.5879e-01 13.8896 < 2.2e-16
## X3           1.2388e+00  2.1705e-01  5.7075 1.146e-08
## X4           1.2611e-03  2.8269e-03  0.4461   0.65552
## lag.X1       2.1497e-06  4.0740e-06  0.5277   0.59773
## lag.X2       6.5633e-01  4.2880e-01  1.5306   0.12587
## lag.X3       4.4442e-01  4.1494e-01  1.0710   0.28415
## lag.X4       5.3121e-03  3.1428e-03  1.6902   0.09098
## 
## Rho: -0.16852, LR test value: 0.84603, p-value: 0.35768
## Asymptotic standard error: 0.16696
##     z-value: -1.0094, p-value: 0.31279
## Wald statistic: 1.0189, p-value: 0.31279
## 
## Log likelihood: -38.16477 for mixed model
## ML residual variance (sigma squared): 0.4324, (sigma: 0.65757)
## Number of observations: 38 
## Number of parameters estimated: 11 
## AIC: 98.33, (AIC for lm: 97.176)
## LM test for residual autocorrelation
## test value: 1.0534, p-value: 0.30472

Interpretasi:

Hasil estimasi Spatial Durbin Model (SDM) menunjukkan bahwa tidak terdapat ketergantungan spasial yang signifikan, baik melalui variabel dependen maupun melalui spillover variabel independen. Hal ini tercermin dari nilai koefisien spasial ρ (rho) sebesar −0,1685 yang tidak signifikan, serta seluruh koefisien lag variabel independen yang juga tidak signifikan pada tingkat 5 persen. Secara parsial, hanya variabel X2 dan X3 yang berpengaruh positif dan signifikan terhadap IPM, sementara X1 dan X4 tidak signifikan, dan pola ini konsisten dengan hasil model sebelumnya. Meskipun model SDM menghasilkan varians residual yang relatif lebih kecil, namun tidak diimbangi dengan peningkatan kinerja yang lebih baik, yang ditunjukkan oleh nilai AIC 98,33 yang lebih tinggi. Oleh karena itu, model SDM tidak memberikan nilai tambah dalam menjelaskan struktur spasial IPM di Jawa Timur.

4.3.5 Spatial Durbin Error Model (SDEM)

sdem_jatim <- errorsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(sdem_jatim)

## 
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, 
##     listw = lw_jatim, Durbin = TRUE, method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.421883 -0.409431  0.031956  0.340269  1.552964 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  3.3538e+01  4.5922e+00  7.3031 2.811e-13
## X1          -1.9891e-06  1.5160e-06 -1.3120    0.1895
## X2           2.1735e+00  1.6576e-01 13.1123 < 2.2e-16
## X3           1.2391e+00  2.1721e-01  5.7045 1.167e-08
## X4           7.5016e-04  2.9237e-03  0.2566    0.7975
## lag.X1       2.3505e-06  3.5196e-06  0.6678    0.5042
## lag.X2       2.5286e-01  2.0253e-01  1.2485    0.2118
## lag.X3       2.8709e-01  3.7181e-01  0.7721    0.4400
## lag.X4       4.4559e-03  3.2259e-03  1.3813    0.1672
## 
## Lambda: -0.23553, LR test value: 1.5239, p-value: 0.21703
## Asymptotic standard error: 0.16779
##     z-value: -1.4037, p-value: 0.1604
## Wald statistic: 1.9704, p-value: 0.1604
## 
## Log likelihood: -37.82583 for error model
## ML residual variance (sigma squared): 0.42102, (sigma: 0.64886)
## Number of observations: 38 
## Number of parameters estimated: 11 
## AIC: 97.652, (AIC for lm: 97.176)

Interpretasi:

Hasil estimasi Spatial Durbin Error Model (SDEM) menunjukkan bahwa tidak terdapat ketergantungan spasial yang signifikan pada komponen error, maupun efek spillover dari variabel independen. Hal ini ditunjukkan oleh parameter λ (lambda) sebesar −0,2355 yang tidak signifikan, serta seluruh koefisien lag variabel independen yang juga tidak signifikan pada tingkat 5 persen. Secara parsial, variabel X2 dan X3 berpengaruh positif dan signifikan terhadap IPM, sedangkan X1 dan X4 tidak signifikan. Dari sisi kinerja model, nilai AIC 97,652 lebih tinggi dibandingkan model OLS, sehingga penambahan komponen Durbin pada model error tidak memberikan peningkatan kualitas model yang berarti. Dengan demikian, model SDEM kurang tepat untuk merepresentasikan struktur spasial IPM di Jawa Timur.

4.3.6 Spatial Autoregressive Combined Model (SAC)

sac_jatim <- sacsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sac_jatim)

## 
## Call:
## sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, listw = lw_jatim, 
##     method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.345713 -0.465782 -0.047563  0.437583  1.631009 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  3.1839e+01  3.5215e+00  9.0413 < 2.2e-16
## X1          -1.4537e-06  1.4912e-06 -0.9748    0.3297
## X2           2.2568e+00  1.5037e-01 15.0085 < 2.2e-16
## X3           1.2503e+00  2.1773e-01  5.7424 9.333e-09
## X4           3.7138e-03  2.2805e-03  1.6285    0.1034
## 
## Rho: 0.095475
## Asymptotic standard error: 0.038372
##     z-value: 2.4882, p-value: 0.012841
## Lambda: -0.35314
## Asymptotic standard error: 0.1641
##     z-value: -2.152, p-value: 0.031399
## 
## LR test value: 7.993, p-value: 0.01838
## 
## Log likelihood: -38.94607 for sac model
## ML residual variance (sigma squared): 0.4351, (sigma: 0.65962)
## Number of observations: 38 
## Number of parameters estimated: 8 
## AIC: 93.892, (AIC for lm: 97.885)

Interpretasi:

Hasil estimasi Spatial Autoregressive Combined Model (SAC) menunjukkan bahwa terdapat dua bentuk ketergantungan spasial yang signifikan dalam IPM kabupaten/kota di Provinsi Jawa Timur, yaitu melalui interaksi spasial pada variabel dependen dan korelasi spasial pada komponen error. Hal ini ditunjukkan oleh nilai ρ (rho) sebesar 0,0955 yang signifikan pada tingkat 5 persen, serta λ (lambda) sebesar −0,3531 yang juga signifikan. Uji LR yang signifikan semakin menegaskan bahwa memasukkan kedua komponen spasial tersebut secara simultan memberikan perbaikan model yang bermakna. Secara parsial, variabel X2 dan X3 berpengaruh positif dan signifikan terhadap IPM, sementara X1 dan X4 tidak signifikan, dengan hasil yang konsisten dibandingkan model-model sebelumnya. Dari sisi kinerja model, model SAC merupakan model spesifikasi terbaik, yang ditunjukkan oleh nilai AIC terendah 93,892 dibandingkan OLS. Hal ini mengindikasikan bahwa variasi IPM di Jawa Timur paling baik dijelaskan oleh model yang mengakomodasi sekaligus ketergantungan spasial langsung antarwilayah dan efek spasial dari faktor tak teramati.

4.3.7 General Nested Spatial Model (GNS)

gns_jatim <- sacsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(gns_jatim)

## 
## Call:
## sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_jatim, listw = lw_jatim, 
##     Durbin = TRUE, method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.083842 -0.356290  0.057358  0.263559  1.469029 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  1.6582e+01  8.8340e+00  1.8771   0.06051
## X1          -1.7908e-06  1.3460e-06 -1.3304   0.18337
## X2           2.1593e+00  1.7439e-01 12.3824 < 2.2e-16
## X3           1.2206e+00  2.1127e-01  5.7775 7.583e-09
## X4           1.1699e-03  2.9354e-03  0.3986   0.69022
## lag.X1       1.6503e-06  2.8688e-06  0.5752   0.56513
## lag.X2      -8.2466e-01  5.2337e-01 -1.5757   0.11510
## lag.X3      -2.5281e-01  4.1170e-01 -0.6141   0.53917
## lag.X4       2.1593e-03  3.4118e-03  0.6329   0.52681
## 
## Rho: 0.45137
## Asymptotic standard error: 0.19909
##     z-value: 2.2672, p-value: 0.023378
## Lambda: -0.64369
## Asymptotic standard error: 0.18537
##     z-value: -3.4725, p-value: 0.00051572
## 
## LR test value: 11.898, p-value: 0.064272
## 
## Log likelihood: -36.9933 for sacmixed model
## ML residual variance (sigma squared): 0.32804, (sigma: 0.57275)
## Number of observations: 38 
## Number of parameters estimated: 12 
## AIC: 97.987, (AIC for lm: 97.885)

Interpretasi:

Hasil estimasi General Nested Spatial Model (GNS) menunjukkan adanya ketergantungan spasial yang signifikan baik melalui komponen spatial lag maupun spatial error, yang ditunjukkan oleh nilai ρ (rho) sebesar 0,4514 dan λ (lambda) −0,6437, keduanya signifikan pada tingkat 5 persen. Hal ini mengindikasikan bahwa IPM suatu kabupaten/kota di Jawa Timur dipengaruhi oleh IPM wilayah sekitarnya sekaligus oleh faktor tak teramati yang bersifat spasial. Namun, seluruh koefisien lag variabel independen (spillover effects) tidak signifikan, sehingga tidak terdapat bukti kuat adanya efek tidak langsung dari variabel penjelas antarwilayah. Secara parsial, hanya X2 dan X3 yang berpengaruh positif dan signifikan terhadap IPM, sedangkan X1 dan X4 tidak signifikan. Dari sisi kinerja model, nilai AIC 97,987 lebih tinggi, sehingga peningkatan kompleksitas model tidak diikuti oleh perbaikan kualitas model. Oleh karena itu, model GNS kurang efisien untuk analisis ini.

4.4 Perbandingan Model

# Comparison Table (AIC & BIC)
model_tbl <- data.frame(
  Model = c("OLS", "SAR", "SEM", "SDM", "SDEM", "SAC", "GNS"),
  AIC = c(
    AIC(ols_jatim),
    AIC(sar_jatim),
    AIC(sem_jatim),
    AIC(sdm_jatim),
    AIC(sdem_jatim),
    AIC(sac_jatim),
    AIC(gns_jatim)
  ),
  BIC = c(
    BIC(ols_jatim),
    BIC(sar_jatim),
    BIC(sem_jatim),
    BIC(sdm_jatim),
    BIC(sdem_jatim),
    BIC(sac_jatim),
    BIC(gns_jatim)
  )
)
print(model_tbl)

##   Model      AIC      BIC
## 1   OLS 97.88509 107.7106
## 2   SAR 95.99062 107.4537
## 3   SEM 97.48166 108.9448
## 4   SDM 98.32953 116.3430
## 5  SDEM 97.65165 115.6651
## 6   SAC 93.89214 106.9928
## 7   GNS 97.98660 117.6376

# Find best model
best_aic <- model_tbl$Model[which.min(model_tbl$AIC)]
best_bic <- model_tbl$Model[which.min(model_tbl$BIC)]
cat("\nBest model by AIC:", best_aic, "\n")

## 
## Best model by AIC: SAC

cat("Best model by BIC:", best_bic, "\n")

## Best model by BIC: SAC

Interpretasi:

Tabel perbandingan AIC dan BIC menunjukkan bahwa model SAC (Spatial Autoregressive Combined) memiliki nilai AIC dan BIC terendah dibandingkan seluruh model yang diuji. Hal ini mengindikasikan bahwa model SAC merupakan model terbaik dan paling efisien dalam menjelaskan variasi IPM kabupaten/kota di Provinsi Jawa Timur, karena mampu mencapai keseimbangan terbaik antara goodness of fit dan kompleksitas model. Dengan demikian, baik berdasarkan kriteria AIC maupun BIC, model SAC paling tepat dipilih sebagai model utama untuk analisis IPM di Jawa Timur karena mampu menangkap ketergantungan spasial secara lebih komprehensif.

4.5 Model Terbaik

4.5.1 Spatial Autoregressive Combined Model (SAC)

Didapatkan model terbaik, yaitu Spatial Autoregressive Combined Model (SAC). Dengan model sebagai berikut:

\[ \begin{aligned} y =\;& 0.095475\, Wy + 31.839 - 1.4537 \times 10^{-6}\, X_1 + 2.2568\, X_2 + 1.2503\, X_3 + 0.0037138\, X_4 + u. \end{aligned} \]

dengan komponen error spasialnya: \[ u = -0.35314\, Wu + \boldsymbol{\varepsilon}, \]

dengan ragam error sebesar: \[ \sigma^2 \approx 0.4351 \quad (\sigma \approx 0.6596). \]

Keterangan: - \(y\) = vektor variabel dependen (IPM) - \(X\) = matriks variabel independen (PDRB, RLS, HLS, Jumlah Penduduk Miskin) - \(W\) = matriks pembobot spasial - \(\rho\) = koefisien autoregresif spasial yang menunjukkan pengaruh variabel dependen di wilayah sekitar terhadap wilayah ke-\(i\). - \(\beta\) = vektor koefisien regresi - \(WX\) = spatial lag dari variabel independen - \(u\) = komponen error yang mengandung ketergantungan spasial - \(\lambda\) = koefisien error spasial yang menunjukkan tingkat ketergantungan antar error wilayah - \(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Model Spatial Autoregressive Combined (SAC) menunjukkan bahwa Indeks Pembangunan Manusia (IPM) kabupaten/kota di Provinsi Jawa Timur dipengaruhi tidak hanya oleh faktor-faktor internal wilayah, tetapi juga oleh interaksi spasial antarwilayah serta faktor tak teramati yang bersifat spasial. Koefisien spasial ρ (rho) sebesar 0,095 yang signifikan mengindikasikan adanya ketergantungan spasial positif pada variabel dependen, di mana peningkatan IPM di wilayah sekitar cenderung mendorong peningkatan IPM di suatu kabupaten/kota. Hal ini menunjukkan adanya efek limpahan (spillover effect) pembangunan manusia antarwilayah di Jawa Timur.

Selain itu, nilai koefisien λ (lambda) sebesar −0,353 yang signifikan menunjukkan adanya ketergantungan spasial pada komponen error, yang mengindikasikan bahwa terdapat faktor-faktor lain di luar model—seperti karakteristik sosial, budaya, atau kebijakan lokal—yang memiliki pola spasial dan memengaruhi IPM secara bersamaan di wilayah-wilayah yang berdekatan.

Secara parsial, variabel Rata-Rata Lama Sekolah (RLS) dan Harapan Lama Sekolah (HLS) berpengaruh positif dan signifikan terhadap IPM. Koefisien RLS sebesar 2,257 menunjukkan bahwa setiap peningkatan satu tahun rata-rata lama sekolah akan meningkatkan IPM sebesar sekitar 2,26 poin, dengan asumsi variabel lain konstan. Sementara itu, koefisien HLS sebesar 1,250 mengindikasikan bahwa peningkatan satu tahun harapan lama sekolah akan meningkatkan IPM sebesar sekitar 1,25 poin. Sebaliknya, PDRB dan jumlah penduduk miskin tidak menunjukkan pengaruh yang signifikan secara statistik setelah memperhitungkan ketergantungan spasial. Secara keseluruhan, hasil ini menunjukkan bahwa variasi IPM di Jawa Timur lebih kuat dipengaruhi oleh dimensi pendidikan serta keterkaitan spasial antarwilayah, baik secara langsung maupun melalui faktor laten. Hal ini menegaskan pentingnya kebijakan pembangunan manusia yang terintegrasi secara spasial, khususnya di bidang pendidikan, agar peningkatan kualitas manusia di suatu wilayah dapat memberikan dampak yang lebih luas ke wilayah sekitarnya.

4.5.2 Peta Prediksi IPM Kabupaten/Kota Jawa Timur Dengan SAC

# Map subset (sesuai complete cases)
JaTim_map <- JaTim_merged[cc, ]

# SAC predictions & residuals
JaTim_map$Yhat_SAC  <- as.numeric(sac_jatim$fitted.values)
JaTim_map$resid_SAC <- as.numeric(residuals(sac_jatim))

# Spatial lag of Y (Wy)
JaTim_map$Wy <- as.numeric(spdep::lag.listw(lw_jatim, dat_jatim$Y, zero.policy = TRUE))

# Map of SAC Predictions
ggplot(JaTim_map) +
  geom_sf(aes(fill = Yhat_SAC), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "East Java — SAC Model Predictions", fill = "Y-hat") +
  theme_minimal()

Interpretasi

Peta prediksi IPM berdasarkan model Spatial Autoregressive Combined (SAC) menunjukkan pola persebaran yang sangat mirip dengan IPM aktual, di mana nilai IPM yang lebih tinggi diperkirakan terkonsentrasi di kawasan perkotaan dan pusat pertumbuhan di bagian tengah Jawa Timur, sementara wilayah dengan IPM yang lebih rendah cenderung berada di daerah perdesaan dan wilayah pinggiran. Kesamaan pola ini mengindikasikan bahwa model SAC mampu menangkap variasi spasial IPM secara akurat. Selain itu, hasil prediksi yang relatif halus antarwilayah mencerminkan peran ketergantungan spasial ganda, baik melalui interaksi langsung antarwilayah maupun melalui faktor-faktor tak teramati yang bersifat spasial. Hal ini menunjukkan bahwa nilai IPM suatu kabupaten/kota tidak hanya dipengaruhi oleh karakteristik internal wilayah tersebut, tetapi juga oleh kondisi pembangunan manusia di wilayah sekitarnya. Dengan demikian, model SAC efektif dalam merepresentasikan kompleksitas interaksi spasial dalam pembangunan manusia di Provinsi Jawa Timur dan memberikan gambaran prediksi yang lebih realistis dibandingkan model spasial yang lebih sederhana.

4.5.3 Peta Residual SAC IPM Kabupaten/Kota Jawa Timur

ggplot(JaTim_map) +
  geom_sf(aes(fill = resid_SAC), color = "white", linewidth = 0.2) +
  scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b") +
  labs(title = "East Java — SAC Model Residuals", fill = "Residual") +
  theme_minimal()

Interpretasi

Peta residual hasil model SAC menunjukkan bahwa nilai residual IPM kabupaten/kota di Provinsi Jawa Timur tersebar relatif acak dan tidak membentuk pola spasial yang jelas. Sebagian besar wilayah memiliki residual yang mendekati nol, yang mengindikasikan bahwa perbedaan antara nilai IPM aktual dan nilai prediksi model relatif kecil. Beberapa wilayah memang menunjukkan residual positif atau negatif yang cukup besar, namun pola tersebut bersifat lokal dan tidak saling berkelompok. Kondisi ini menunjukkan bahwa model SAC telah mampu menangkap struktur spasial utama dalam data, sehingga sisa kesalahan prediksi tidak lagi mengandung ketergantungan spasial yang sistematis. Dengan demikian, peta residual ini mengonfirmasi bahwa model SAC memberikan hasil estimasi yang memadai dan tidak menyisakan pola spasial yang signifikan pada residual.

4.6 Model Geographically Weighted Regression (GWR)

# Prepare spatial points
cent_sf <- sf::st_centroid(JaTim_map)

## Warning: st_centroid assumes attributes are constant over geometries

cent_sp <- as(cent_sf, "Spatial")

gwr_sp <- SpatialPointsDataFrame(
  coords = sp::coordinates(cent_sp),
  data   = dat_jatim,
  proj4string = sp::CRS(sf::st_crs(JaTim_map)$wkt)
)

# Project to UTM for proper distance calculation
gwr_utm <- spTransform(gwr_sp, 
                       CRS("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs"))

# Find optimal bandwidth
bw_opt <- gwr.sel(
  Y ~ X1 + X2 + X3 + X4,
  data  = gwr_utm,
  adapt = TRUE
)

## Adaptive q: 0.381966 CV score: 58.5006 
## Adaptive q: 0.618034 CV score: 59.39088 
## Adaptive q: 0.236068 CV score: 58.04362 
## Adaptive q: 0.145898 CV score: 55.23632 
## Adaptive q: 0.09016994 CV score: 55.79186 
## Adaptive q: 0.1357269 CV score: 55.53533 
## Adaptive q: 0.1803399 CV score: 55.56327 
## Adaptive q: 0.1576742 CV score: 54.89829 
## Adaptive q: 0.1663317 CV score: 55.06881 
## Adaptive q: 0.1578453 CV score: 54.89362 
## Adaptive q: 0.1602256 CV score: 54.92379 
## Adaptive q: 0.1587545 CV score: 54.9019 
## Adaptive q: 0.1581648 CV score: 54.89502 
## Adaptive q: 0.157971 CV score: 54.89302 
## Adaptive q: 0.1579304 CV score: 54.89261 
## Adaptive q: 0.1578897 CV score: 54.8924 
## Adaptive q: 0.1578897 CV score: 54.8924

# Fit GWR model
gwr_fit <- gwr(
  Y ~ X1 + X2 + X3 + X4,
  data  = gwr_utm,
  adapt = bw_opt,
  hatmatrix = TRUE,
  se.fit = TRUE
)

summary(gwr_fit)

##           Length Class                  Mode     
## SDF         38   SpatialPointsDataFrame S4       
## lhat      1444   -none-                 numeric  
## lm          11   -none-                 list     
## results     14   -none-                 list     
## bandwidth   38   -none-                 numeric  
## adapt        1   -none-                 numeric  
## hatmatrix    1   -none-                 logical  
## gweight      1   -none-                 character
## gTSS         1   -none-                 numeric  
## this.call    6   -none-                 call     
## fp.given     1   -none-                 logical  
## timings     12   -none-                 numeric

# Model comparison
cat("\nOLS AIC:", AIC(ols_jatim), "\n")

## 
## OLS AIC: 97.88509

cat("GWR AICc:", gwr_fit$results$AICc, "\n")

## GWR AICc: 125.1082

Interpretasi:

Hasil estimasi Geographically Weighted Regression (GWR) menunjukkan bahwa bandwidth adaptif optimal diperoleh pada nilai adaptive q sekitar 0,158, yang mengindikasikan bahwa dalam setiap estimasi lokal digunakan sekitar 16% observasi terdekat. Nilai bandwidth adaptif yang relatif kecil ini menunjukkan bahwa hubungan antara Indeks Pembangunan Manusia (IPM) dan variabel penjelas bersifat lokal dan bervariasi antar kabupaten/kota di Provinsi Jawa Timur. Dengan kata lain, pengaruh variabel seperti PDRB, pendidikan, dan kemiskinan terhadap IPM tidak bersifat homogen secara spasial, melainkan berbeda-beda antarwilayah. Namun, nilai AICc GWR sebesar 125,11 jauh lebih tinggi dibandingkan nilai AIC model OLS yang menunjukkan bahwa secara keseluruhan model GWR tidak memberikan peningkatan kinerja model dibandingkan model regresi global. Perbedaan nilai AIC/AICc yang cukup besar ini mengindikasikan bahwa kompleksitas model GWR tidak diimbangi dengan peningkatan goodness of fit yang memadai. Dengan demikian, meskipun GWR mampu mengungkap adanya potensi variasi pengaruh lokal antarwilayah, model ini kurang efisien untuk menjelaskan pola IPM Jawa Timur secara agregat.

# Extract local coefficients
gwr_coef <- gwr_fit$SDF@data

JaTim_map$gwr_b0 <- gwr_coef[["(Intercept)"]]
JaTim_map$gwr_x1 <- gwr_coef[["X1"]]
JaTim_map$gwr_x2 <- gwr_coef[["X2"]]
JaTim_map$gwr_x3 <- gwr_coef[["X3"]]
JaTim_map$gwr_x4 <- gwr_coef[["X4"]]
JaTim_map$resid_gwr <- gwr_fit$SDF$gwr.e

# Maps of local coefficients
ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x1), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X1)", fill="β_X1") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x2), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X2)", fill="β_X2") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x3), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X3)", fill="β_X3") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x4), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X4)", fill="β_X4") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = resid_gwr), color="white", size=0.2) +
  scale_fill_gradient2(low="#2166ac", mid="white", high="#b2182b") +
  labs(title="East Java — GWR Residuals", fill="Residual") +
  theme_minimal()

Interpretasi:

Hasil Geographically Weighted Regression (GWR) menunjukkan bahwa pengaruh variabel penjelas terhadap Indeks Pembangunan Manusia (IPM) bervariasi secara spasial antar kabupaten/kota di Provinsi Jawa Timur. Koefisien lokal PDRB (X1) terlihat relatif kecil dan cenderung homogen di sebagian besar wilayah, yang mengindikasikan bahwa pengaruh faktor ekonomi regional terhadap IPM tidak dominan dan relatif lemah pada hampir seluruh kabupaten/kota. Hal ini sejalan dengan hasil model global yang juga menunjukkan bahwa PDRB tidak berpengaruh signifikan secara statistik. Sebaliknya, Rata-Rata Lama Sekolah (X2) menunjukkan variasi koefisien yang cukup jelas antarwilayah, dengan pengaruh yang relatif lebih kuat di wilayah bagian tengah dan utara Jawa Timur dibandingkan wilayah selatan dan timur. Hal ini mengindikasikan bahwa kontribusi pendidikan formal terhadap peningkatan IPM berbeda-beda antar kabupaten/kota. Pola serupa juga terlihat pada Harapan Lama Sekolah (X3), di mana koefisien lokal cenderung lebih besar di beberapa wilayah tertentu, khususnya di bagian barat dan selatan Jawa Timur, menunjukkan bahwa ekspektasi pendidikan jangka panjang memiliki peran yang lebih kuat pada wilayah-wilayah tersebut.Sementara itu, Jumlah Penduduk Miskin (X4) memperlihatkan variasi koefisien lokal yang cukup beragam, baik dari sisi besaran maupun intensitas pengaruh, meskipun secara umum pengaruhnya relatif kecil. Hal ini menunjukkan bahwa dampak kemiskinan terhadap IPM tidak seragam antarwilayah dan sangat bergantung pada karakteristik lokal masing-masing kabupaten/kota. Peta residual GWR menunjukkan bahwa sebagian besar wilayah memiliki residual yang relatif kecil dan tersebar tanpa membentuk pola spasial yang kuat, meskipun masih terdapat beberapa lokasi dengan residual yang cukup besar. Secara keseluruhan, hal ini menegaskan adanya heterogenitas spasial dalam faktor-faktor penentu IPM di Jawa Timur.

4.7 Interpolasi Spasial

4.7.1 Inverse Distance Weighting (IDW)

# 8.1 Prepare observation points
obs_sp <- SpatialPointsDataFrame(
  coords = coordinates(cent_sp),
  data   = dat_jatim,
  proj4string = CRS(st_crs(JaTim_map)$wkt)
)

obs_utm <- spTransform(obs_sp, 
                       CRS("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs"))

# Create prediction grid
bb <- bbox(obs_utm)
grid_df <- expand.grid(
  x = seq(bb[1,1], bb[1,2], length.out = 200),
  y = seq(bb[2,1], bb[2,2], length.out = 200)
)

# 8.2 Inverse Distance Weighting (IDW)
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]
    z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

X <- coordinates(obs_utm)
z <- obs_utm@data$Y

grid_df$Y_idw <- mapply(
  function(x0, y0) idw_predict(X[,1], X[,2], z, x0, y0, p = 2, k = 8),
  grid_df$x, grid_df$y
)

ggplot(grid_df, aes(x, y, fill = Y_idw)) +
  geom_raster() +
  coord_equal() +
  scale_fill_viridis_c() +
  labs(title = "IDW Interpolation (p=2) — East Java", 
       x = "UTM X", y = "UTM Y", fill = "Y-hat") +
  theme_minimal()

Interpretasi:

Hasil interpolasi Inverse Distance Weighting (IDW) dengan parameter p=2 menunjukkan pola permukaan IPM yang halus dan kontinu di Provinsi Jawa Timur. Nilai IPM yang lebih tinggi (ditunjukkan oleh warna lebih terang) terkonsentrasi pada area dengan kepadatan titik observasi yang tinggi, khususnya di kawasan tengah–barat Jawa Timur, sedangkan nilai IPM yang lebih rendah (warna lebih gelap) cenderung muncul di wilayah pinggiran dan area yang relatif jauh dari titik pengamatan. Pola ini mencerminkan asumsi dasar metode IDW, yaitu bahwa wilayah yang berdekatan secara geografis memiliki karakteristik IPM yang lebih mirip dibandingkan wilayah yang lebih jauh. Dengan demikian, peta IDW memberikan gambaran visual awal mengenai gradien spasial IPM.

4.7.2 Thin Plate Splines (TPS)

tps_fit <- fields::Tps(X, z)

grid_sp <- SpatialPoints(grid_df[,c("x","y")], 
                         proj4string = CRS(proj4string(obs_utm)))
pred_tps <- predict(tps_fit, coordinates(grid_sp))

grid_df$Y_tps <- as.numeric(pred_tps)

ggplot(grid_df, aes(x, y, fill = Y_tps)) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_equal() +
  labs(title = "Thin Plate Splines — East Java", 
       x = "UTM X", y = "UTM Y", fill = "Y-hat") +
  theme_minimal()

Interpretasi:

Hasil interpolasi Thin Plate Splines (TPS) menunjukkan permukaan IPM yang sangat halus dan kontinu di Provinsi Jawa Timur. Nilai IPM yang lebih tinggi diperkirakan terkonsentrasi di wilayah bagian tengah, kemudian secara bertahap menurun ke arah wilayah pinggiran. Pola gradien yang mulus ini mencerminkan kemampuan metode TPS dalam menangkap tren spasial global tanpa menghasilkan fluktuasi lokal yang tajam. Dengan demikian, peta TPS memberikan gambaran umum mengenai kecenderungan spasial IPM di Jawa Timur, namun cenderung menghaluskan variasi lokal sehingga lebih tepat digunakan untuk visualisasi tren regional dibandingkan untuk analisis detail pada tingkat kabupaten/kota.

4.7.3 Ordinary Kriging

# Handle duplicate coordinates
xy <- coordinates(obs_utm)
key <- paste(round(xy[,1], 3), round(xy[,2], 3), sep = "_")

df_obs <- data.frame(
  key = key,
  x   = xy[,1],
  y   = xy[,2],
  Y   = obs_utm$Y
)

agg_Y <- aggregate(Y ~ key, data = df_obs, FUN = mean)
xy_unique <- df_obs[!duplicated(df_obs$key), c("key","x","y")]
agg2 <- merge(agg_Y, xy_unique, by = "key", all.x = TRUE)

obs_utm2 <- SpatialPointsDataFrame(
  coords = agg2[, c("x","y")],
  data   = data.frame(Y = agg2$Y),
  proj4string = CRS(proj4string(obs_utm))
)

# Fit variogram
vg_emp <- variogram(Y ~ 1, data = obs_utm2, cutoff = 250000, width = 20000)

vg_init <- vgm(
  psill  = var(obs_utm2$Y, na.rm = TRUE) * 0.5,
  model  = "Sph",
  range  = 80000,
  nugget = var(obs_utm2$Y, na.rm = TRUE) * 0.5
)

vg_fit <- fit.variogram(vg_emp, vg_init, fit.method = 2)

## Warning in fit.variogram(vg_emp, vg_init, fit.method = 2): No convergence after
## 200 iterations: try different initial values?

plot(vg_emp, vg_fit, main = "Empirical Variogram + Spherical Model")

# Prepare grid for kriging
coordinates(grid_df) <- ~ x + y
proj4string(grid_df) <- CRS(proj4string(obs_utm2))

grid_pixdf <- SpatialPixelsDataFrame(
  points = SpatialPixels(grid_df),
  data   = data.frame(dummy = rep(1, nrow(coordinates(grid_df)))),
  proj4string = CRS(proj4string(obs_utm2))
)

# Ordinary kriging
ok_jatim <- krige(
  Y ~ 1,
  locations = obs_utm2,
  newdata   = grid_pixdf,
  model     = vg_fit,
  nmax      = 20,
  maxdist   = 150000
)

## [using ordinary kriging]

spplot(ok_jatim["var1.pred"], main = "Ordinary Kriging — East Java")

# Cross-validation
cv_ok <- krige.cv(
  formula   = Y ~ 1,
  locations = obs_utm2,
  model     = vg_fit,
  nfold     = nrow(obs_utm2@data),
  nmax      = 20,
  maxdist   = 150000
)

cv_df <- as.data.frame(cv_ok)
rmse_ok <- sqrt(mean((cv_df$observed - cv_df$var1.pred)^2, na.rm = TRUE))
cat("\nOrdinary Kriging RMSE (LOO):", rmse_ok, "\n")

## 
## Ordinary Kriging RMSE (LOO): 4.325255

ggplot(cv_df, aes(x = observed, y = var1.pred)) +
  geom_point(alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
  theme_minimal() +
  labs(
    title = paste0("Kriging Cross-validation (RMSE = ", round(rmse_ok, 2), ")"),
    x = "Observed Y",
    y = "Predicted Y (LOO)"
  )

Interpretasi:

Hasil analisis variogram empiris menunjukkan bahwa semivarian IPM meningkat seiring bertambahnya jarak antarwilayah hingga mencapai kondisi mendatar, yang mengindikasikan adanya ketergantungan spasial pada IPM kabupaten/kota di Provinsi Jawa Timur. Pemodelan dengan variogram spherical mampu menangkap pola tersebut secara umum, meskipun muncul peringatan tidak tercapainya konvergensi penuh, yang menunjukkan bahwa struktur spasial IPM relatif kompleks dan sensitif terhadap pemilihan parameter awal.

Berdasarkan variogram tersebut, estimasi Ordinary Kriging menghasilkan permukaan prediksi IPM yang halus, dengan nilai IPM lebih tinggi terkonsentrasi di wilayah bagian tengah Jawa Timur dan menurun secara bertahap ke arah wilayah pinggiran. Pola ini konsisten dengan hasil interpolasi sebelumnya dan mencerminkan adanya gradien spasial IPM yang kuat. Hasil validasi silang menunjukkan nilai RMSE sebesar 4,33, yang menandakan tingkat kesalahan prediksi sedang. Plot cross-validation memperlihatkan bahwa sebagian besar titik prediksi berada di sekitar garis diagonal, meskipun masih terdapat penyimpangan pada beberapa observasi. Secara keseluruhan, hasil ini menunjukkan bahwa metode kriging mampu menangkap struktur spasial IPM dengan cukup baik, namun ketelitiannya masih dipengaruhi oleh keterbatasan jumlah dan sebaran titik observasi.

4.8 Spatial Analysis Dashboard

Dashboard ini menyajikan analisis ekonometrika spasial untuk memahami pola dan determinasi Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur. Analisis dilakukan dengan membandingkan berbagai model regresi, yaitu OLS (Ordinary Least Squares) sebagai model regresi konvensional, Spatial Autoregressive Model (SAR) yang menangkap pengaruh limpahan spasial melalui variabel dependen, Spatial Error Model (SEM) yang memperhitungkan ketergantungan spasial pada komponen galat, serta pengembangan model spasial lanjutan seperti Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS).

Variabel dependen yang dianalisis adalah Indeks Pembangunan Manusia (IPM), sedangkan variabel independen meliputi Produk Domestik Regional Bruto (PDRB), Rata-Rata Lama Sekolah (RLS), Harapan Lama Sekolah (HLS), dan jumlah penduduk miskin. Dashboard ini dilengkapi dengan berbagai fitur visualisasi dan analisis, seperti peta persebaran IPM aktual dan hasil prediksi model spasial, perbandingan kinerja model berdasarkan indikator statistik (AIC dan BIC), analisis autokorelasi spasial global dan lokal (Moran’s I, Geary’s C, LISA, dan Getis-Ord Gi*), serta visualisasi residual model untuk mengevaluasi kualitas estimasi.

Selain itu, dashboard ini juga menyediakan analisis tambahan berupa Geographically Weighted Regression (GWR) untuk mengeksplorasi variasi pengaruh lokal antarwilayah, serta interpolasi spasial menggunakan metode Inverse Distance Weighting (IDW), Thin Plate Splines (TPS), dan Ordinary Kriging untuk menggambarkan pola IPM secara kontinu. Untuk menjalankan dashboard ini, pengguna perlu menyiapkan data IPM dan variabel penjelas dalam format CSV yang memuat kolom kabupaten/kota dan variabel penelitian, serta file batas wilayah kabupaten/kota Provinsi Jawa Timur dalam format spasial (.rds atau shapefile).

Link Dashboard : https://sefia-herlina724.shinyapps.io/dashboard/ Link Video: https://youtu.be/jJ_Y5WhiM5k?si=pv1lOvt9UddasP7s

BAB V (KESIMPULAN DAN SARAN)

5.1 Kesimpulan

Berdasarkan hasil analisis, dapat disimpulkan bahwa persebaran Indeks Pembangunan Manusia (IPM) antar kabupaten/kota di Provinsi Jawa Timur membentuk pola spasial yang tidak merata dan cenderung mengelompok. Peta tematik menunjukkan bahwa wilayah dengan IPM tinggi terkonsentrasi di kawasan perkotaan dan pusat pertumbuhan (terutama wilayah tengah–utara), sedangkan wilayah dengan IPM lebih rendah cenderung berada di daerah perdesaan dan wilayah pinggiran. Hal ini diperkuat oleh hasil uji autokorelasi spasial, di mana Moran’s I bernilai positif dan signifikan serta Geary’s C lebih kecil dari 1 dan signifikan, yang menegaskan adanya ketergantungan spasial (spatial dependence) pada IPM di Jawa Timur. Analisis lokal menggunakan Getis-Ord Gi* dan Local Moran’s I (LISA) juga menunjukkan bahwa klaster spasial yang signifikan hanya muncul pada beberapa wilayah tertentu, baik dalam bentuk klaster nilai tinggi (High–High) maupun nilai rendah (Low–Low), sementara sebagian besar wilayah tidak menunjukkan signifikansi lokal yang kuat.

Hasil estimasi model ekonometrika menunjukkan bahwa secara parsial variabel pendidikan, yaitu Rata-Rata Lama Sekolah (RLS) dan Harapan Lama Sekolah (HLS), berpengaruh positif dan signifikan terhadap IPM, sedangkan PDRB dan jumlah penduduk miskin tidak signifikan setelah efek spasial diperhitungkan. Dari berbagai spesifikasi model yang dibandingkan (SAR, SEM, SDM, SDEM, SAC, dan GNS), model terbaik yang paling tepat menjelaskan determinasi IPM di Jawa Timur adalah Spatial Autoregressive Combined (SAC), yang ditunjukkan oleh nilai AIC dan BIC terendah serta signifikansi parameter ρ dan λ. Hal ini mengindikasikan bahwa variasi IPM di Jawa Timur dipengaruhi oleh interaksi antarwilayah secara langsung melalui IPM tetangga sekaligus oleh faktor-faktor tak teramati yang memiliki pola spasial. Dengan demikian, pembangunan manusia di Jawa Timur tidak hanya ditentukan oleh karakteristik internal daerah, tetapi juga oleh keterkaitan spasial antar kabupaten/kota, sehingga kebijakan peningkatan pembangunan manusia perlu dirancang secara terintegrasi berbasis wilayah, terutama melalui penguatan aspek pendidikan pada daerah-daerah yang tertinggal.

5.2 Saran

Berdasarkan hasil penelitian ini, terdapat beberapa saran yang dapat diberikan:

Pemerintah daerah perlu menerapkan kebijakan pembangunan manusia berbasis wilayah
Peningkatan kualitas dan pemerataan pendidikan harus menjadi prioritas utama pembangunan manusia di Kabupaten/Kota Jawa Timur
Penelitian selanjutnya disarankan untuk menggunakan data panel dan variabel tambahan, seperti kualitas layanan kesehatan, infrastruktur, dan belanja pemerintah, sehingga mampu menangkap dinamika pembangunan manusia secara lebih komprehensif.

LAMPIRAN

Syntax R

# Load Required Packages
library(sf)
library(sp)
library(spdep)
library(spatialreg)
library(dplyr)
library(stringr)
library(ggplot2)
library(spgwr)
library(gstat)
library(car)
library(lmtest)
library(FNN)
library(deldir)
library(fields)

# =========================================================
# 1. DATA PREPARATION
# =========================================================

# Load CSV data
jatim_dat <- read.csv("C:/Users/hp/Downloads/Spasial UAS - JaTim.csv", 
                      check.names = FALSE)

cat("Original CSV rows:", nrow(jatim_dat), "\n")

# Load GADM polygon data (level 2: Kabupaten/Kota)
IDN_kab <- readRDS("C:/Users/hp/Downloads/Dashboard/gadm36_IDN_2_small.rds")

# Filter for East Java province
JaTim_sp <- IDN_kab[IDN_kab$NAME_1 == "Jawa Timur", ]
cat("GADM East Java polygons:", nrow(JaTim_sp), "\n")

JaTim_sp$id <- seq_len(nrow(JaTim_sp))
JaTim_sf <- st_as_sf(JaTim_sp) |> st_make_valid()

# Enhanced cleaning function that preserves more distinctions
clean_nm <- function(x){
  x |>
    str_to_lower() |>
    str_trim() |>
    # Remove "Kabupaten" or "Kab." prefix but keep "Kota"
    str_replace_all("^kabupaten\\s+", "") |>
    str_replace_all("^kab\\.\\s+", "") |>
    str_replace_all("^kab\\s+", "") |>
    # Standardize "Kota" prefix
    str_replace_all("^kota\\s+", "kota ") |>
    # Clean up punctuation but preserve spaces
    str_replace_all("\\.", "") |>
    str_replace_all(",", "") |>
    str_squish()
}

# Apply cleaning
jatim_dat <- jatim_dat |>
  mutate(
    key_nm = clean_nm(`Kab/Kota`),
    original_name = `Kab/Kota`
  )

JaTim_sf <- JaTim_sf |>
  mutate(
    key_nm = clean_nm(NAME_2),
    original_name = NAME_2
  )

# Check matching status
cat("\n=== MATCHING CHECK ===\n")
cat("Unique regions in CSV:", length(unique(jatim_dat$key_nm)), "\n")
cat("Unique regions in shapefile:", length(unique(JaTim_sf$key_nm)), "\n")

# Show all names for comparison
cat("\n--- CSV Names (cleaned) ---\n")
csv_names <- jatim_dat |>
  select(original_name, key_nm) |>
  arrange(key_nm) |>
  tibble::as_tibble()
print(csv_names, n = Inf)

cat("\n--- Shapefile Names (cleaned) ---\n")
shp_names <- JaTim_sf |> 
  st_drop_geometry() |>
  select(original_name, key_nm) |> 
  arrange(key_nm) |>
  tibble::as_tibble()
print(shp_names, n = Inf)

# Show unmatched regions from CSV
cek_unmatch_csv <- anti_join(jatim_dat, st_drop_geometry(JaTim_sf), by = "key_nm")
if(nrow(cek_unmatch_csv) > 0) {
  cat("\n⚠️ UNMATCHED in CSV (", nrow(cek_unmatch_csv), "regions):\n")
  print(cek_unmatch_csv[, c("original_name", "key_nm")])
}

# Show unmatched regions from shapefile
cek_unmatch_shp <- anti_join(st_drop_geometry(JaTim_sf), jatim_dat, by = "key_nm")
if(nrow(cek_unmatch_shp) > 0) {
  cat("\n⚠️ UNMATCHED in Shapefile (", nrow(cek_unmatch_shp), "regions):\n")
  print(cek_unmatch_shp[, c("original_name", "key_nm")])
}

# Manual matching dictionary for problematic cases
manual_fixes <- data.frame(
  csv_name = c(
    # Add manual mappings here if needed after seeing the output
    # Example: "special case in CSV"
  ),
  shp_name = c(
    # Corresponding shapefile names
    # Example: "special case in shapefile"
  ),
  stringsAsFactors = FALSE
)

# Apply manual fixes if any exist
if(nrow(manual_fixes) > 0) {
  # Create lookup for CSV
  csv_lookup <- setNames(manual_fixes$shp_name, manual_fixes$csv_name)
  
  jatim_dat <- jatim_dat |>
    mutate(key_nm_fixed = ifelse(
      key_nm %in% names(csv_lookup),
      csv_lookup[key_nm],
      key_nm
    ))
  
  # Use fixed names for merging
  JaTim_merged <- JaTim_sf |>
    inner_join(jatim_dat, by = c("key_nm" = "key_nm_fixed"))
} else {
  # Standard merge without manual fixes
  JaTim_merged <- JaTim_sf |>
    inner_join(jatim_dat, by = "key_nm")
}

cat("\n=== AFTER MERGE ===\n")
cat("Final observations:", nrow(JaTim_merged), "\n")
cat("Missing Y values:", sum(is.na(JaTim_merged$Y)), "\n")

# Verify we have 38 observations
if(nrow(JaTim_merged) != 38) {
  cat("\n!!! WARNING: Expected 38 observations but got", nrow(JaTim_merged), "!!!\n")
  cat("\nDIAGNOSTIC INFO:\n")
  cat("- CSV has", nrow(jatim_dat), "rows\n")
  cat("- Shapefile has", nrow(JaTim_sf), "polygons\n")
  cat("- Matched regions:", nrow(JaTim_merged), "\n")
  cat("- Unmatched CSV:", nrow(cek_unmatch_csv), "\n")
  cat("- Unmatched Shapefile:", nrow(cek_unmatch_shp), "\n")
  cat("\nPlease check the name lists above and add manual fixes if needed.\n\n")
  
  # Stop execution if count is wrong
  stop("Observation count mismatch! Expected 38, got ", nrow(JaTim_merged))
} else {
  cat("\n✅ SUCCESS: Got exactly 38 observations as expected!\n\n")
}

# =========================================================
# 2. EXPLORATORY SPATIAL DATA ANALYSIS
# =========================================================

# 2.1 Basic Statistics
summary(
  JaTim_merged |>
    dplyr::select(Y, X1, X2, X3, X4) |>
    sf::st_drop_geometry()
)

# 2.2 Map of Variable Y
ggplot(JaTim_merged) +
  geom_sf(aes(fill = Y), color="white", size=0.2) +
  scale_fill_viridis_c(option="magma", na.value="grey85") +
  labs(title="East Java — Distribution of Y Variable (n=38)", fill="Y") +
  theme_minimal()

# =========================================================
# 3. SPATIAL AUTOCORRELATION ANALYSIS
# =========================================================
# Queen contiguity
nb_queen <- poly2nb(JaTim_merged, queen = TRUE)

coords <- st_coordinates(st_centroid(JaTim_merged))

# Fix isolates with 1-NN
iso <- which(card(nb_queen) == 0)
if(length(iso) > 0){
  knn_nb <- knn2nb(knearneigh(coords, k = 1))
  for(i in iso) nb_queen[[i]] <- knn_nb[[i]]
}

# Final weights (NO zero.policy needed anymore)
lwW <- nb2listw(nb_queen, style = "W", zero.policy = FALSE)
lwB <- nb2listw(nb_queen, style = "B", zero.policy = FALSE)

cat("Isolated:", sum(card(nb_queen) == 0), "\n")
cat("Avg neighbors:", mean(card(nb_queen)), "\n")

moran_res <- moran.test(JaTim_merged$Y, lwW,
                        randomisation = TRUE,
                        alternative = "two.sided")
print(moran_res)

geary_res <- geary.test(JaTim_merged$Y, lwW,
                        randomisation = TRUE,
                        alternative = "two.sided")
print(geary_res)

x <- as.numeric(scale(JaTim_merged$Y))
lagx <- lag.listw(lwW, x)

lisa <- localmoran(x, lwW)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")

alpha <- 0.05
quad_raw <- 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)"
)

JaTim_LISA <- dplyr::bind_cols(JaTim_merged, lisa_df) |>
  dplyr::mutate(
    quad = ifelse(Pi.two.sided <= alpha, quad_raw, "Not significant")
  )
ggplot(JaTim_LISA) +
  geom_sf(aes(fill = quad), color = "white", linewidth = 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 = "East Java — Local Moran's I (LISA)", fill = "Category") +
  theme_minimal()

Gi_z <- localG(JaTim_merged$Y, lwW)

JaTim_G <- JaTim_merged |>
  dplyr::mutate(
    z_Gistar = as.numeric(Gi_z),
    hotcold = dplyr::case_when(
      z_Gistar >=  1.96 ~ "Hot spot (p≈0.05)",
      z_Gistar <= -1.96 ~ "Cold spot (p≈0.05)",
      TRUE ~ "Not significant"
    )
  )
ggplot(JaTim_G) +
  geom_sf(aes(fill = hotcold), color = "white", linewidth = 0.2) +
  scale_fill_manual(values = c(
    "Hot spot (p≈0.05)" = "#b2182b",
    "Cold spot (p≈0.05)" = "#2166ac",
    "Not significant" = "grey85"
  )) +
  labs(title = "East Java — Getis-Ord Gi* Hot/Cold Spots", fill = NULL) +
  theme_minimal()

# =========================================================
# 4. SPATIAL ECONOMETRIC MODELS
# =========================================================

vars_use <- c("Y","X1","X2","X3","X4")

dat_jatim <- JaTim_merged |>
  dplyr::select(dplyr::all_of(vars_use)) |>
  sf::st_drop_geometry()

cc <- complete.cases(dat_jatim)
dat_jatim <- dat_jatim[cc, , drop = FALSE]

# SUBSET WEIGHTS YANG SAMA
lw_jatim <- spdep::subset.listw(lwW, subset = cc)

ols_jatim <- lm(Y ~ X1 + X2 + X3 + X4, data = dat_jatim)
summary(ols_jatim)

lm.morantest(ols_jatim, lw_jatim)
lm.LMtests(ols_jatim, lw_jatim,
           test = c("LMlag","LMerr","RLMlag","RLMerr"))

sar_jatim <- lagsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sar_jatim)

sem_jatim <- errorsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sem_jatim)

sdm_jatim <- lagsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(sdm_jatim)

sdem_jatim <- errorsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(sdem_jatim)

sac_jatim <- sacsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  method = "eigen"
)
summary(sac_jatim)

gns_jatim <- sacsarlm(
  Y ~ X1 + X2 + X3 + X4,
  data = dat_jatim,
  listw = lw_jatim,
  Durbin = TRUE,
  method = "eigen"
)
summary(gns_jatim)

# =========================================================
# 5. MODEL COMPARISON AND SELECTION
# =========================================================

# 5.1 Comparison Table (AIC & BIC)
model_tbl <- data.frame(
  Model = c("OLS", "SAR", "SEM", "SDM", "SDEM", "SAC", "GNS"),
  AIC = c(
    AIC(ols_jatim),
    AIC(sar_jatim),
    AIC(sem_jatim),
    AIC(sdm_jatim),
    AIC(sdem_jatim),
    AIC(sac_jatim),
    AIC(gns_jatim)
  ),
  BIC = c(
    BIC(ols_jatim),
    BIC(sar_jatim),
    BIC(sem_jatim),
    BIC(sdm_jatim),
    BIC(sdem_jatim),
    BIC(sac_jatim),
    BIC(gns_jatim)
  )
)
print(model_tbl)

# Find best model
best_aic <- model_tbl$Model[which.min(model_tbl$AIC)]
best_bic <- model_tbl$Model[which.min(model_tbl$BIC)]
cat("\nBest model by AIC:", best_aic, "\n")
cat("Best model by BIC:", best_bic, "\n")

# =========================================================
# 6. VISUALIZATION OF BEST MODEL (SAC)
# =========================================================

# Map subset (sesuai complete cases)
JaTim_map <- JaTim_merged[cc, ]

# Actual Y
JaTim_map$Y_actual <- dat_jatim$Y

# SAC predictions & residuals
JaTim_map$Yhat_SAC  <- as.numeric(sac_jatim$fitted.values)
JaTim_map$resid_SAC <- as.numeric(residuals(sac_jatim))

# Spatial lag of Y (Wy)
JaTim_map$Wy <- as.numeric(spdep::lag.listw(lw_jatim, dat_jatim$Y, zero.policy = TRUE))

# ---------------------------------------------------------
# 6.1 Map of Actual Y
ggplot(JaTim_map) +
  geom_sf(aes(fill = Y_actual), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "East Java — Actual Y Values", fill = "Y") +
  theme_minimal()

# 6.2 Map of SAC Predictions
ggplot(JaTim_map) +
  geom_sf(aes(fill = Yhat_SAC), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "East Java — SAC Model Predictions", fill = "Y-hat") +
  theme_minimal()

# 6.3 Map of SAC Residuals
ggplot(JaTim_map) +
  geom_sf(aes(fill = resid_SAC), color = "white", linewidth = 0.2) +
  scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b") +
  labs(title = "East Java — SAC Model Residuals", fill = "Residual") +
  theme_minimal()

# 6.4 Map of Spatial Spillover (Wy)
ggplot(JaTim_map) +
  geom_sf(aes(fill = Wy), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "plasma") +
  labs(title = "East Java — Spatial Lag (Wy)", fill = "Wy") +
  theme_minimal()

# ---------------------------------------------------------
# Check residual spatial autocorrelation (SAC residuals)
moran_resid_sac <- spdep::moran.test(
  JaTim_map$resid_SAC,
  lw_jatim,
  randomisation = TRUE,
  alternative = "two.sided",
  zero.policy = TRUE
)
print(moran_resid_sac)


# =========================================================
# 7. GEOGRAPHICALLY WEIGHTED REGRESSION (GWR)
# =========================================================

# 7.1 Prepare spatial points
cent_sf <- sf::st_centroid(JaTim_map)
cent_sp <- as(cent_sf, "Spatial")

gwr_sp <- SpatialPointsDataFrame(
  coords = sp::coordinates(cent_sp),
  data   = dat_jatim,
  proj4string = sp::CRS(sf::st_crs(JaTim_map)$wkt)
)

# Project to UTM for proper distance calculation
gwr_utm <- spTransform(gwr_sp, 
                       CRS("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs"))

# 7.2 Find optimal bandwidth
bw_opt <- gwr.sel(
  Y ~ X1 + X2 + X3 + X4,
  data  = gwr_utm,
  adapt = TRUE
)

# 7.3 Fit GWR model
gwr_fit <- gwr(
  Y ~ X1 + X2 + X3 + X4,
  data  = gwr_utm,
  adapt = bw_opt,
  hatmatrix = TRUE,
  se.fit = TRUE
)

summary(gwr_fit)

# 7.4 Model comparison
cat("\nOLS AIC:", AIC(ols_jatim), "\n")
cat("GWR AICc:", gwr_fit$results$AICc, "\n")

# 7.5 Extract local coefficients
gwr_coef <- gwr_fit$SDF@data

JaTim_map$gwr_b0 <- gwr_coef[["(Intercept)"]]
JaTim_map$gwr_x1 <- gwr_coef[["X1"]]
JaTim_map$gwr_x2 <- gwr_coef[["X2"]]
JaTim_map$gwr_x3 <- gwr_coef[["X3"]]
JaTim_map$gwr_x4 <- gwr_coef[["X4"]]
JaTim_map$resid_gwr <- gwr_fit$SDF$gwr.e

# 7.6 Maps of local coefficients
ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x1), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X1)", fill="β_X1") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x2), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X2)", fill="β_X2") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x3), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X3)", fill="β_X3") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = gwr_x4), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="East Java — GWR Local Coefficient (X4)", fill="β_X4") +
  theme_minimal()

ggplot(JaTim_map) +
  geom_sf(aes(fill = resid_gwr), color="white", size=0.2) +
  scale_fill_gradient2(low="#2166ac", mid="white", high="#b2182b") +
  labs(title="East Java — GWR Residuals", fill="Residual") +
  theme_minimal()

# =========================================================
# 8. SPATIAL INTERPOLATION METHODS
# =========================================================

# 8.1 Prepare observation points
obs_sp <- SpatialPointsDataFrame(
  coords = coordinates(cent_sp),
  data   = dat_jatim,
  proj4string = CRS(st_crs(JaTim_map)$wkt)
)

obs_utm <- spTransform(obs_sp, 
                       CRS("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs"))

# Create prediction grid
bb <- bbox(obs_utm)
grid_df <- expand.grid(
  x = seq(bb[1,1], bb[1,2], length.out = 200),
  y = seq(bb[2,1], bb[2,2], length.out = 200)
)

# 8.2 Inverse Distance Weighting (IDW)
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  if (!is.null(k) && k < length(d)) {
    idx <- order(d)[1:k]
    d <- d[idx]
    z <- z[idx]
  }
  w <- 1 / (d^p)
  sum(w * z) / sum(w)
}

X <- coordinates(obs_utm)
z <- obs_utm@data$Y

grid_df$Y_idw <- mapply(
  function(x0, y0) idw_predict(X[,1], X[,2], z, x0, y0, p = 2, k = 8),
  grid_df$x, grid_df$y
)

ggplot(grid_df, aes(x, y, fill = Y_idw)) +
  geom_raster() +
  coord_equal() +
  scale_fill_viridis_c() +
  labs(title = "IDW Interpolation (p=2) — East Java", 
       x = "UTM X", y = "UTM Y", fill = "Y-hat") +
  theme_minimal()

# 8.3 Thin Plate Splines (TPS)
tps_fit <- fields::Tps(X, z)

grid_sp <- SpatialPoints(grid_df[,c("x","y")], 
                         proj4string = CRS(proj4string(obs_utm)))
pred_tps <- predict(tps_fit, coordinates(grid_sp))

grid_df$Y_tps <- as.numeric(pred_tps)

ggplot(grid_df, aes(x, y, fill = Y_tps)) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_equal() +
  labs(title = "Thin Plate Splines — East Java", 
       x = "UTM X", y = "UTM Y", fill = "Y-hat") +
  theme_minimal()

# 8.4 Ordinary Kriging
# Handle duplicate coordinates
xy <- coordinates(obs_utm)
key <- paste(round(xy[,1], 3), round(xy[,2], 3), sep = "_")

df_obs <- data.frame(
  key = key,
  x   = xy[,1],
  y   = xy[,2],
  Y   = obs_utm$Y
)

agg_Y <- aggregate(Y ~ key, data = df_obs, FUN = mean)
xy_unique <- df_obs[!duplicated(df_obs$key), c("key","x","y")]
agg2 <- merge(agg_Y, xy_unique, by = "key", all.x = TRUE)

obs_utm2 <- SpatialPointsDataFrame(
  coords = agg2[, c("x","y")],
  data   = data.frame(Y = agg2$Y),
  proj4string = CRS(proj4string(obs_utm))
)

# Fit variogram
vg_emp <- variogram(Y ~ 1, data = obs_utm2, cutoff = 250000, width = 20000)

vg_init <- vgm(
  psill  = var(obs_utm2$Y, na.rm = TRUE) * 0.5,
  model  = "Sph",
  range  = 80000,
  nugget = var(obs_utm2$Y, na.rm = TRUE) * 0.5
)

vg_fit <- fit.variogram(vg_emp, vg_init, fit.method = 2)

plot(vg_emp, vg_fit, main = "Empirical Variogram + Spherical Model")

# Prepare grid for kriging
coordinates(grid_df) <- ~ x + y
proj4string(grid_df) <- CRS(proj4string(obs_utm2))

grid_pixdf <- SpatialPixelsDataFrame(
  points = SpatialPixels(grid_df),
  data   = data.frame(dummy = rep(1, nrow(coordinates(grid_df)))),
  proj4string = CRS(proj4string(obs_utm2))
)

# Ordinary kriging
ok_jatim <- krige(
  Y ~ 1,
  locations = obs_utm2,
  newdata   = grid_pixdf,
  model     = vg_fit,
  nmax      = 20,
  maxdist   = 150000
)

spplot(ok_jatim["var1.pred"], main = "Ordinary Kriging — East Java")

# Cross-validation
cv_ok <- krige.cv(
  formula   = Y ~ 1,
  locations = obs_utm2,
  model     = vg_fit,
  nfold     = nrow(obs_utm2@data),
  nmax      = 20,
  maxdist   = 150000
)

cv_df <- as.data.frame(cv_ok)
rmse_ok <- sqrt(mean((cv_df$observed - cv_df$var1.pred)^2, na.rm = TRUE))
cat("\nOrdinary Kriging RMSE (LOO):", rmse_ok, "\n")

ggplot(cv_df, aes(x = observed, y = var1.pred)) +
  geom_point(alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
  theme_minimal() +
  labs(
    title = paste0("Kriging Cross-validation (RMSE = ", round(rmse_ok, 2), ")"),
    x = "Observed Y",
    y = "Predicted Y (LOO)"
  )