ANALISIS SPASIAL Kondisi Konektivitas Transportasi Indonesia tahun 2024
Disusun oleh: Raymond Emmanuel Krista (140610230067)
Dosen Pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.
Program Studi Statistika
Fakultas Matematika dan Ilmu Pengetahuan Alam
Universitas Padjadjaran
Jatinangor 2025
Bab 1. Pendahuluan
1.1 Latar Belakang
Konektivitas transportasi merupakan salah satu faktor penting dalam mendukung pertumbuhan ekonomi, pemerataan pembangunan, serta peningkatan kesejahteraan masyarakat. Infrastruktur jalan, khususnya jalan dengan kondisi baik, berperan penting dalam memperlancar mobilitas barang dan jasa, menurunkan biaya logistik, serta memperkuat integrasi antarwilayah. Di negara kepulauan seperti Indonesia, kualitas dan pemerataan infrastruktur jalan antar provinsi menjadi isu strategis dalam pembangunan nasional. Selain itu, konektivitas transportasi juga mendukung point 9 dan 11 SDGs yaitu Industry, Innovation, and Infrastructure serta Sustainable Cities and Communities.
Meskipun pembangunan infrastruktur jalan di Indonesia terus mengalami peningkatan, distribusi kondisi jalan yang baik antar provinsi belum tentu merata. Perbedaan karakteristik wilayah seperti tingkat aktivitas ekonomi, luas geografis, jumlah penduduk, serta tingkat kemiskinan berpotensi menyebabkan ketimpangan kualitas infrastruktur jalan. Selain itu, pembangunan jalan di suatu provinsi tidak berdiri sendiri, melainkan dapat dipengaruhi oleh kondisi wilayah sekitarnya, baik melalui keterkaitan jaringan transportasi, arus ekonomi antarwilayah, maupun kebijakan pembangunan regional.
Sejumlah penelitian terdahulu menunjukkan bahwa variabel infrastruktur, termasuk jalan, memiliki keterkaitan spasial yang signifikan. Penelitian terdahulu mengatakan bahwa fenomena pembangunan wilayah sering kali menunjukkan autokorelasi spasial, di mana kondisi suatu wilayah berkorelasi dengan kondisi wilayah tetangganya [1][2]. Dalam konteks infrastruktur transportasi, ditemukan bahwa kualitas jalan di suatu daerah dipengaruhi oleh pembangunan infrastruktur di daerah sekitarnya melalui efek limpahan (spillover effects) [3][4].
Selain itu, data kondisi jalan umumnya tersedia pada unit wilayah administratif, sementara fenomena pembangunan infrastruktur berlangsung secara kontinu di ruang geografis. Oleh karena itu, selain pemodelan berbasis wilayah, diperlukan pendekatan interpolasi spasial untuk menggambarkan variasi kondisi jalan secara lebih halus di seluruh wilayah Indonesia. Metode interpolasi spasial seperti Inverse Distance Weighting (IDW), Ordinary Kriging, dan Trend Surface Analysis dapat digunakan untuk mengestimasi nilai kondisi jalan pada lokasi yang tidak terobservasi secara langsung, sekaligus memberikan visualisasi pola spasial yang lebih komprehensif.
Dengan demikian, analisis kondisi panjang jalan dengan kondisi baik antarprovinsi di Indonesia tidak memadai jika hanya menggunakan pendekatan regresi klasik. Diperlukan pendekatan analisis spasial yang terintegrasi, meliputi pengujian autokorelasi spasial, pemodelan ekonometrika spasial, serta interpolasi spasial, agar mampu menangkap ketergantungan antarwilayah dan menggambarkan pola sebaran kondisi jalan secara lebih representatif.
1.2 Identifikasi Masalah
Berdasarkan latar belakang yang telah diuraikan, beberapa pertanyaan penelitian yang menjadi fokus utama adalah sebagai berikut:
Bagaimana pola sebaran spasial panjang jalan dengan kondisi baik antarprovinsi di Indonesia?
Apakah terdapat autokorelasi spasial dalam distribusi panjang jalan kondisi baik, baik secara global maupun lokal?
Model ekonometrika spasial apa yang paling sesuai untuk menjelaskan pengaruh faktor sosial-ekonomi dan demografis terhadap kondisi jalan antarprovinsi?
Bagaimana pola spasial kontinu dari kondisi jalan jika divisualisasikan menggunakan metode interpolasi spasial?
Metode interpolasi spasial manakah yang memberikan hasil estimasi terbaik dalam memodelkan variasi kondisi jalan ?
1.3 Tujuan Penelitian
Sejalan dengan identifikasi masalah, tujuan dari penelitian ini adalah:
Menganalisis pola sebaran spasial panjang jalan dengan kondisi baik di tingkat provinsi di Indonesia.
Mengidentifikasi keberadaan autokorelasi spasial menggunakan statistik Moran’s I, Geary’s C, Getis-Ord Gi*, dan analisis LISA.
Menentukan model ekonometrika spasial yang paling sesuai dalam menjelaskan hubungan antara kondisi jalan dan faktor-faktor penjelas seperti PDRB per kapita, luas wilayah, jumlah penduduk, dan tingkat kemiskinan.
Menerapkan metode interpolasi spasial, yaitu Inverse Distance Weighting (IDW), Ordinary Kriging, dan Trend Surface Analysis, untuk memvisualisasikan variasi kondisi jalan secara kontinu di wilayah Indonesia.
Membandingkan kinerja metode interpolasi spasial berdasarkan ukuran kesalahan prediksi untuk menentukan metode interpolasi yang paling representatif dalam menggambarkan pola spasial kondisi jalan.
1.4 Batasan Penelitian
Untuk menjaga agar penelitian ini tetap fokus dan mendalam, maka ditetapkan beberapa batasan sebagai berikut:
Unit analisis dalam penelitian ini adalah provinsi di Indonesia, sehingga variasi kondisi jalan pada level kabupaten/kota tidak dianalisis secara eksplisit. Setiap provinsi direpresentasikan oleh nilai agregat dan titik centroid wilayah.
Periode waktu penelitian terbatas pada tahun 2024, sehingga penelitian ini bersifat cross-sectional dan tidak menganalisis dinamika perubahan kondisi jalan dari waktu ke waktu.
Variabel dependen yang digunakan adalah panjang jalan dengan kondisi baik pada masing-masing provinsi. Variabel ini diperlakukan sebagai indikator utama konektivitas transportasi darat.
Variabel independen yang digunakan dalam pemodelan ekonometrika spasial dibatasi pada faktor-faktor makro, yaitu:
Produk Domestik Regional Bruto (PDRB) per kapita,
Luas wilayah,
Jumlah penduduk,
Tingkat kemiskinan.
Pendekatan spasial yang digunakan terbatas pada:
Pengujian autokorelasi spasial global (Moran’s I dan Geary’s C) dan lokal (LISA dan Getis-Ord Gi*),
Pemodelan ekonometrika spasial (SAR, SEM, SDEM, SDM, SAC),
Metode interpolasi spasial deterministik dan geostatistik, yaitu Inverse Distance Weighting (IDW), Ordinary Kriging, dan Trend Surface Analysis.
Interpolasi spasial dilakukan berdasarkan titik centroid provinsi dan bertujuan untuk menggambarkan pola sebaran spasial secara kontinu. Hasil interpolasi tidak dimaksudkan sebagai estimasi nilai aktual pada lokasi tertentu, melainkan sebagai pendekatan visual dan eksploratif terhadap variasi spasial kondisi jalan.
Pemilihan model interpolasi terbaik didasarkan pada kriteria statistik kesalahan prediksi seperti Root Mean Square Error (RMSE), Mean Absolute Error (MAE), dan Mean Error (ME), tanpa mempertimbangkan aspek biaya komputasi atau implementasi kebijakan secara operasional.
Bab 2. Tinjauan Pustaka
2.1 Dependensi Spasial
Dependensi spasial atau ketergantungan spasial, adalah kondisi dimana nilai suatu variabel pada satu lokasi memiliki korelasi dengan nilai variabel yang sama di lokasi-lokasi terdekatnya (tetangga). Dalam konteks sosial-ekonomi, fenomena ini dapat terjadi karena berbagai mekanisme, seperti efek limpahan (spillover), interaksi antarwilayah, persaingan, atau adanya faktor-faktor tak teramati yang memiliki struktur spasial.
Untuk merepresentasikan struktur ketetanggaan secara kuantitatif, digunakan Matriks Pembobot Spasial (W) yang mendefinisikan hubungan spasial antar unit observasi. Matriks ini dapat dibangun berdasarkan kriteria kontiguitas wilayah, jarak geografis, maupun k-nearest neighbors, dan berperan penting dalam hampir seluruh metode analisis spasial modern [5].
2.2 Autokorelasi Spasial
Autokorelasi spasial adalah ukuran formal untuk mengidentifikasi ada atau tidaknya pengelompokan (klaster) atau pola sistematis dalam distribusi data spasial.
Secara umum, autokorelasi spasial dibedakan menjadi autokorelasi global dan lokal. Autokorelasi global memberikan gambaran umum pola spasial pada seluruh wilayah studi, sedangkan autokorelasi lokal berfokus pada pola di sekitar masing-masing wilayah [6].
2.2.1 Matriks Bobot Spasial (W)
Matriks pembobot spasial (W) merupakan komponen kunci dalam analisis spasial yang merepresentasikan struktur hubungan atau ketetanggaan antar unit wilayah. Dalam bentuk matriks berukuran n x n, elemen wij menjelaskan derajat hubungan antara wilayah i dan j, di mana nilai yang lebih besar menunjukkan keterkaitan spasial yang lebih kuat [7]. Matriks ini digunakan untuk menentukan seberapa besar kontribusi nilai variabel dari lokasi tetangga terhadap lokasi pusat dalam berbagai teknik analisis spasial, termasuk pengujian autokorelasi dan model ekonometrika spasial. Pemilihan definisi ketetanggaan sangat memengaruhi hasil analisis, sehingga perlu disesuaikan dengan karakter data dan fenomena yang dianalisis [8][9].
Salah satu pendekatan populer untuk mendefinisikan ketetanggaan adalah melalui k-nearest neighbors (kNN), di mana setiap unit wilayah memiliki k tetangga terdekat berdasarkan ukuran jarak tertentu (misalnya Euclidean). Dalam pendekatan ini, bobot spasial bersifat simetris dan terfokus pada sejumlah terbatas tetangga terdekat, sehingga sangat cocok untuk dataset dengan distribusi spasial yang tidak homogen dan batas administratif yang komplek, seperti provinsi atau kabupaten/kota [10][11].
\[ \mathbf{W} = \begin{bmatrix} w_{11} & w_{12} & \cdots & w_{1n} \\ w_{21} & w_{22} & \cdots & w_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n1} & w_{n2} & \cdots & w_{nn} \end{bmatrix} \] \[ w_{ij} = \begin{cases} 1, & \text{jika } j \in k\text{-nearest neighbors dari } i \\ 0, & \text{lainnya} \end{cases} \] \[ w_{ij}^* = \frac{w_{ij}}{\sum_{j=1}^{n}w_{ij}} \] dengan,
\[ \begin{aligned} \mathbf{W} &: \text{matriks pembobot spasial } n\times n \\ w_{ij} &: \text{elemen pembobot antara wilayah } i \text{ dan } j \\ k &: \text{jumlah tetangga terdekat yang dipilih dalam kNN} \\ w_{ij}^* &: \text{bobot spasial yang telah dinormalisasi (row-standard)} \\ n &: \text{jumlah unit wilayah dalam analisis} \end{aligned} \]
2.2.2 Moran’s I
Indeks Moran (Moran’s I) adalah statistik yang paling umum digunakan untuk menguji autokorelasi spasial secara global di seluruh wilayah studi. Nilai Moran’s I berkisar antara -1 hingga +1 [12].
Nilai I > 0: Menandakan autokorelasi spasial positif, artinya wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi lainnya (pola mengelompok/klaster).
Nilai I < 0: Menandakan autokorelasi spasial negatif, artinya wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai rendah (pola menyebar/terdispersi).
Nilai I ≈ 0: Menandakan tidak ada autokorelasi spasial (pola acak).
2.2.2 Geary’s C
Geary’s C adalah ukuran autokorelasi spasial global lainnya yang sering digunakan sebagai pelengkap Moran’s I. Perbedaan utamanya adalah, jika Moran’s I menggunakan produk silang dari deviasi terhadap rata-rata, Geary’s C menggunakan kuadrat perbedaan nilai antara dua lokasi yang bertetangga. Dengan kata lain, Geary’s C lebih fokus pada disimilaritas atau variasi antar tetangga [13].
Nilai Geary’s C berkisar dari 0 hingga 2, dengan interpretasi yang berkebalikan dari Moran’s I:
0 < C < 1: Menandakan autokorelasi spasial positif (pola mengelompok). Semakin dekat nilai C ke 0, semakin kuat pengelompokannya.
C ≈ 1: Menandakan tidak ada autokorelasi spasial (pola acak).
1 < C < 2: Menandakan autokorelasi spasial negatif (pola menyebar). Semakin dekat nilai C ke 2, semakin kuat pola penyebarannya.
2.2.3 Getis-Ord Gi*
Statistik lokal Getis-Ord Gi* adalah ukuran autokorelasi spasial yang digunakan untuk mengidentifikasi cluster lokal dari nilai tinggi (hot spots) dan nilai rendah (cold spots) dalam data spasial [14][15]. Berbeda dengan statistik global seperti Moran’s I yang memberikan gambaran keseluruhan pola spasial, Gi* fokus pada intensitas lokal nilai variabel pada suatu wilayah serta tetangganya. Metode ini banyak dipakai dalam studi ekonomi regional, epidemiologi, geografi kesehatan, dan analisis infrastruktur untuk memetakan area prioritas intervensi [16][17]
Getis-Ord Gi* mendeteksi lokasi yang memiliki nilai variabel signifikan secara spasial lebih tinggi atau lebih rendah dari ekspektasi acak berdasarkan persentil statistik, sehingga sangat berguna dalam memetakan ketimpangan spasial maupun konsentrasi pembangunan [15].
\[ G_i^* = \frac{ \displaystyle\sum_{j=1}^{n} w_{ij} x_j - \bar{X} \displaystyle\sum_{j=1}^{n} w_{ij} }{S \sqrt{ \left( n \displaystyle\sum_{j=1}^{n} w_{ij}^2 - \left(\displaystyle\sum_{j=1}^{n} w_{ij}\right)^2 \right)/(n-1) }} \]
\[ \bar{X} = \frac{1}{n}\sum_{j=1}^{n} x_j \]
\[ S = \sqrt{ \frac{\displaystyle\sum_{j=1}^{n} x_j^2}{n} - \bar{X}^2 } \]
dengan,
\[ \begin{aligned} G_i^* &: \text{statistik Getis-Ord Gi* untuk lokasi } i \\ x_j &: \text{nilai variabel di lokasi } j \\ w_{ij} &: \text{elemen matriks bobot spasial antara } i \text{ dan } j \\ n &: \text{jumlah total wilayah observasi} \\ \bar{X} &: \text{rata-rata nilai variabel seluruh wilayah} \\ S &: \text{standar deviasi termodifikasi untuk Gi*} \end{aligned} \]
2.2.4 LISA
Sementara Moran’s I dan Geary’s C memberikan satu nilai ringkasan untuk seluruh area, Local Indicators of Spatial Association (LISA) memungkinkan identifikasi klaster lokal atau titik-titik anomali (outlier) spasial. LISA mengurai statistik global menjadi kontribusi spesifik dari setiap lokasi [18]. Hasil analisis LISA biasanya divisualisasikan dalam Peta Moran (Moran Map) yang mengklasifikasikan setiap wilayah ke dalam salah satu dari empat kuadran:
High-High (Hot Spot): Wilayah bernilai tinggi dikelilingi oleh wilayah bernilai tinggi.
Low-Low (Cold Spot): Wilayah bernilai rendah dikelilingi oleh wilayah bernilai rendah.
High-Low (Spatial Outlier): Wilayah bernilai tinggi dikelilingi oleh wilayah bernilai rendah.
Low-High (Spatial Outlier): Wilayah bernilai rendah dikelilingi oleh wilayah bernilai tinggi.
2.3 Model Spasial Ekonometrika
Model ekonometrika spasial digunakan ketika data menunjukkan adanya ketergantungan spasial antar unit observasi, yang menyebabkan pelanggaran asumsi independensi dalam regresi klasik. Dalam konteks pembangunan wilayah dan infrastruktur, fenomena spasial seperti efek limpahan dan keterkaitan antarwilayah sering kali tidak dapat diabaikan [19].
Pendekatan ekonometrika spasial memungkinkan pemodelan eksplisit interaksi antarwilayah melalui variabel dependen, variabel independen, maupun struktur galat, sehingga menghasilkan estimasi parameter yang lebih efisien dan inferensi statistik yang lebih valid [20].
2.3.1 Ordinary Least Squares (OLS)
Model Ordinary Least Squares (OLS) merupakan model dasar dalam analisis regresi yang mengasumsikan tidak adanya ketergantungan antar observasi. Dalam studi spasial, OLS tetap digunakan sebagai model awal (baseline) untuk mengidentifikasi adanya pelanggaran asumsi klasik sebelum menerapkan model spasial lanjutan [21].
\[ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \] dengan, \[ \begin{aligned} \mathbf{Y} &: \text{vektor variabel dependen berukuran } n \times 1 \\ \mathbf{X} &: \text{matriks variabel independen berukuran } n \times k \\ \boldsymbol{\beta} &: \text{vektor koefisien regresi} \\ \boldsymbol{\varepsilon} &: \text{vektor galat acak} \end{aligned} \]
Uji diagnostik terhadap residual OLS, khususnya uji autokorelasi spasial, menjadi langkah penting dalam menentukan apakah pendekatan ekonometrika spasial diperlukan [22].
2.3.1.1 Uji Normalitas
Uji normalitas residual digunakan untuk menilai apakah distribusi residual model regresi mengikuti distribusi normal. Normalitas residual merupakan salah satu asumsi klasik dalam OLS karena estimasi parameter dan uji signifikansi mengasumsikan galat berdistribusi normal [23][24]. Shapiro–Wilk adalah salah satu uji normalitas yang paling sensitif terhadap deviasi dari normalitas, terutama pada ukuran sampel kecil hingga menengah.
\[ W = \frac{\left(\sum_{i=1}^{n} a_i x_{(i)}\right)^2} {\sum_{i=1}^{n}(x_i - \bar{x})^2} \] dengan,
\[ \begin{aligned} W &: \text{statistik Shapiro–Wilk} \\ x_{(i)} &: \text{residual OLS yang diurutkan} \\ a_i &: \text{koefisien berdasarkan nilai harapan normal order statistics} \\ \bar{x} &: \text{rata-rata residual} \\ n &: \text{jumlah observasi} \end{aligned} \]
2.3.1.2 Uji Heteroskedastistias
Heteroskedastisitas terjadi ketika varians residual tidak konstan di seluruh rentang prediktor, melanggar asumsi OLS dan menghasilkan estimator yang tidak efisien. Uji Breusch–Pagan sering digunakan untuk mendeteksi heteroskedastisitas dengan menguji apakah varians residual dipengaruhi oleh variabel independen.
\[ BP = \frac{1}{2}\sum_{i=1}^{n} \left( \hat{u}_i^2 - \bar{\hat{u}}^2 \right)^2 \]\[ \hat{u}_i^2 = \alpha_0 + \alpha_1 X_{1i} + \cdots + \alpha_k X_{ki} + v_i \]
dengan,
\[ \begin{aligned} BP &: \text{statistik Breusch–Pagan} \\ \hat{u}_i^2 &: \text{kuadrat residual OLS} \\ \bar{\hat{u}}^2 &: \text{rata-rata kuadrat residual} \\ X_{ji} &: \text{variabel independen ke-}j \\ v_i &: \text{error regresi auxiliary} \end{aligned} \]
2.3.1.3 Uji Multikolinearitas
Multikolinearitas terjadi ketika dua atau lebih prediktor dalam regresi linier memiliki korelasi tinggi, yang menyebabkan varians koefisien meningkat dan interpretasi menjadi tidak stabil. Variance Inflation Factor (VIF) adalah ukuran yang umum digunakan untuk mendeteksi multikolinearitas. VIF tinggi menunjukkan prediktor memiliki redundansi informasi [25][26].
\[ VIF_j = \frac{1}{1 - R_j^2} \] dengan,
\[ \begin{aligned} VIF_j &: \text{variance inflation factor untuk variabel } j \\ R_j^2 &: \text{koefisien determinasi dari regresi } X_j \text{ pada regresor lain} \end{aligned} \]Kriteria nilai VIF yang sesuai dan tidak melanggar asumsi adalah < 10.
2.3.2 Spatial Autoregressive Model (SAR)
Model Spatial Autoregressive (SAR) digunakan untuk memodelkan dependensi spasial yang bersifat substantif, yaitu ketika nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah tetangganya. Model ini banyak digunakan untuk menganalisis interaksi antarwilayah dan efek limpahan dalam pembangunan ekonomi dan infrastruktur [27].
\[ \mathbf{Y} = \rho \mathbf{W}\mathbf{Y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \] dengan,
\[ \begin{aligned} \rho &: \text{koefisien lag spasial variabel dependen} \\ \mathbf{W} &: \text{matriks pembobot spasial} \\ \mathbf{W}\mathbf{Y} &: \text{lag spasial variabel dependen} \\ \mathbf{X} &: \text{matriks variabel independen} \\ \boldsymbol{\varepsilon} &: \text{galat acak} \end{aligned} \]
2.3.3 Spatial Error Model (SEM)
Model Spatial Error Model (SEM) digunakan ketika autokorelasi spasial berasal dari struktur galat, yang mencerminkan adanya variabel tidak teramati atau faktor eksternal yang memiliki pola spasial. Model ini umum digunakan dalam studi kebijakan publik dan regional ketika interaksi antarwilayah tidak terjadi secara langsung melalui variabel dependen [19].
\[ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{u} \]
\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]
.dengan,
\[ \begin{aligned} \lambda &: \text{koefisien autokorelasi spasial pada galat} \\ \mathbf{u} &: \text{galat yang mengandung struktur spasial} \\ \mathbf{W} &: \text{matriks pembobot spasial} \\ \boldsymbol{\varepsilon} &: \text{galat acak independen} \end{aligned} \]
2.3.4 Spatial Durbin Model (SDM)
Spatial Durbin Model (SDM) merupakan model yang fleksibel karena memasukkan lag spasial dari variabel dependen dan variabel independen. Model ini memungkinkan pengukuran efek langsung serta efek limpahan (spillover effects) dari variabel independen antarwilayah [2].
\[ \mathbf{Y} = \rho \mathbf{W}\mathbf{Y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \] dengan,
\[ \begin{aligned} \rho &: \text{koefisien lag spasial variabel dependen} \\ \mathbf{W}\mathbf{X} &: \text{lag spasial variabel independen} \\ \boldsymbol{\theta} &: \text{koefisien lag spasial variabel independen} \\ \mathbf{X} &: \text{matriks variabel independen} \\ \boldsymbol{\varepsilon} &: \text{galat acak} \end{aligned} \]
2.3.5 Spatial Durbin Error Model (SDEM)
Spatial Durbin Error Model (SDEM) merupakan pengembangan dari SEM dengan memasukkan lag spasial dari variabel independen. Model ini sesuai digunakan ketika autokorelasi spasial terjadi pada galat, namun terdapat efek limpahan dari karakteristik wilayah tetangga [19].
Model SDEM banyak diaplikasikan dalam analisis infrastruktur dan ekonomi regional karena mampu menangkap efek tidak langsung tanpa mengasumsikan adanya interaksi langsung pada variabel dependen [28].
\[ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \mathbf{u} \]
\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]
dengan,
\[ \begin{aligned} \mathbf{W}\mathbf{X} &: \text{lag spasial variabel independen} \\ \boldsymbol{\theta} &: \text{koefisien lag spasial variabel independen} \\ \lambda &: \text{koefisien autokorelasi spasial galat} \\ \mathbf{u} &: \text{komponen galat spasial} \end{aligned} \]
2.3.6 Spatial Autocorrelation Model (SAC)
Model Spatial Autocorrelation (SAC) merupakan model paling umum karena menggabungkan lag spasial pada variabel dependen dan autokorelasi spasial pada galat secara simultan. Model ini digunakan ketika terdapat bukti empiris bahwa kedua mekanisme spasial tersebut terjadi bersamaan [29].
Model SAC sering digunakan dalam studi lanjutan sebagai pembanding untuk menilai kompleksitas proses spasial yang terjadi dalam data wilayah [30].
\[ \mathbf{Y} = \rho \mathbf{W}\mathbf{Y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{u} \]
\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]
dengan,
\[ \begin{aligned} \rho &: \text{koefisien lag spasial variabel dependen} \\ \lambda &: \text{koefisien autokorelasi spasial galat} \\ \mathbf{W} &: \text{matriks pembobot spasial} \\ \mathbf{u} &: \text{galat dengan struktur spasial} \end{aligned} \]
2.4 Interpolasi Spasial
Interpolasi spasial merupakan teknik untuk memperkirakan nilai suatu variabel pada lokasi yang tidak teramati berdasarkan nilai pada lokasi yang diketahui, dengan memanfaatkan prinsip kedekatan geografis [31]. Metode ini banyak digunakan dalam analisis geografi dan infrastruktur untuk menggambarkan pola sebaran fenomena spasial secara kontinu, meskipun data awal tersedia dalam bentuk diskret atau agregat wilayah.
Dalam penelitian infrastruktur transportasi, interpolasi spasial digunakan sebagai alat eksploratif untuk memvisualisasikan variasi kondisi jalan antarwilayah dan mengidentifikasi pola global maupun lokal. Metode interpolasi dapat bersifat deterministik, seperti IDW dan Trend Surface, atau geostatistik seperti Kriging yang mempertimbangkan struktur autokorelasi spasial.
2.4.1 Inverse Distance Wighting (IDW)
Inverse Distance Weighting (IDW) merupakan metode interpolasi deterministik yang mengasumsikan bahwa pengaruh suatu titik pengamatan terhadap lokasi estimasi akan menurun seiring bertambahnya jarak. Nilai interpolasi pada lokasi S0 dihitung sebagai rata-rata berbobot dari nilai di lokasi sekitar, dengan bobot berbasis jarak [32].
Secara matematis, IDW dirumuskan sebagai:
\([\hat{Z}(s_0) = \frac{\displaystyle\sum_{i=1}^{n} w_i \, Z(s_i)}{\displaystyle\sum_{i=1}^{n} w_i},\quad \text{dengan} \ w_i = \frac{1}{d(s_0, s_i)^p}]\)
dengan,
\[ \begin{aligned} \hat{Z}(s_0) &: \text{nilai estimasi pada lokasi } s_0 \\ Z(s_i) &: \text{nilai variabel pada lokasi ke-}i \\ d(s_0, s_i) &: \text{jarak antara } s_0 \text{ dan } s_i \\ p &: \text{parameter pangkat (power parameter)} \\ n &: \text{jumlah titik tetangga} \end{aligned} \]
2.4.2 Ordinary Kriging
Ordinary Kriging merupakan metode interpolasi geostatistik yang memperkirakan nilai pada lokasi tak teramati dengan mempertimbangkan ketergantungan spasial antar pengamatan melalui variogram. Metode ini mengasumsikan bahwa nilai rata-rata variabel konstan tetapi tidak diketahui di seluruh wilayah studi [33].
Estimasi Ordinary Kriging dirumuskan sebagai:
\[ \hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i \, Z(s_i) \]
\[ \sum_{i=1}^{n} \lambda_i = 1 \]
\[ \sum_{j=1}^{n} \lambda_j \, \gamma(s_i - s_j) + \mu = \gamma(s_i - s_0), \quad i = 1, 2, \ldots, n \]
\[ \begin{aligned} \lambda_i &: \text{bobot kriging lokasi ke-}i \\ \gamma(h) &: \text{fungsi variogram pada jarak } h \\ \mu &: \text{pengali Lagrange} \\ s_0 &: \text{lokasi estimasi} \end{aligned} \]
2.4.3 Trend Surface Analysis
Trend Surface Analysis merupakan metode interpolasi deterministik yang memodelkan nilai variabel spasial sebagai fungsi polinomial dari koordinat geografis. Metode ini bertujuan untuk menangkap pola tren global atau kecenderungan umum dalam sebaran spasial suatu fenomena [34].
Secara umum, model Trend Surface dapat dituliskan sebagai:
\[ Z(s) = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 x^2 + \beta_4 y^2 + \beta_5 xy + \varepsilon \]
\[ \begin{aligned} x, y &: \text{koordinat geografis lokasi } s \\ \beta_0, \ldots, \beta_5 &: \text{parameter regresi} \\ \varepsilon &: \text{galat acak} \end{aligned} \]
Bab 3. Metodologi Penelitian
3.1 Jenis dan Sumber Data
Penelitian ini menggunakan data sekunder yang bersifat kuantitatif dengan dimensi cross-sectional pada tahun 2024. Sumber data utama berasal dari publikasi resmi Badan Pusat Statistik (BPS) Indonesia, serta lembaga terkait lainnya. Data yang dikumpulkan terdiri dari dua komponen utama:
Data Atribut/Statistik: Data tabular yang memuat nilai variabel penelitian untuk setiap provinsi di Indonesia, meliputi kondisi infrastruktur jalan dan indikator sosial-ekonomi.
Data Spasial: Data peta digital dalam format json (.kson) yang memuat batas administrasi provinsi di Indonesia. Data spasial ini digunakan untuk mendefinisikan hubungan ketetanggaan antar provinsi, membangun matriks pembobot spasial, serta melakukan visualisasi dan interpolasi spasial.
3.2 Variabel Penelitian
Variabel dalam penelitian ini dibedakan menjadi variabel dependen dan variabel independen.
3.2.1 Variabel Dependen (Y)
Variabel dependen dalam penelitian ini adalah Panjang Jalan Kondisi Baik.
Panjang jalan kondisi baik merepresentasikan total panjang jaringan jalan yang berada dalam kondisi mantap atau baik secara teknis di setiap provinsi. Variabel ini digunakan sebagai indikator utama konektivitas transportasi dan kualitas infrastruktur jalan, yang berperan penting dalam mendukung mobilitas, aktivitas ekonomi, serta integrasi antarwilayah.
Nilai variabel ini diukur dalam satuan kilometer (km) dan diperoleh dari publikasi resmi BPS [35].
3.2.2 Variabel Independen (X)
Variabel independen merupakan faktor-faktor sosial dan ekonomi yang secara teoretis diduga memengaruhi kondisi dan pembangunan infrastruktur jalan antar provinsi. Variabel yang digunakan dalam penelitian ini disajikan pada Tabel berikut.
| Dimensi | Nama Variabel | Simbol | Definisi Operasional | Sumber Data |
|---|---|---|---|---|
| Ekonomi | PDRB per Kapita ADHK | X1 | Nilai Produk Domestik Regional Bruto per kapita atas dasar harga konstan, sebagai indikator tingkat aktivitas dan kapasitas ekonomi wilayah | BPS [36] |
| Geografis | Luas Wilayah | X2 | Total luas administrasi provinsi dalam satuan km², mencerminkan tantangan geografis dalam penyediaan infrastruktur jalan | BPS [37] |
| Demografi | Jumlah Penduduk | X3 | Jumlah penduduk yang berdomisili di setiap provinsi, menggambarkan tingkat kebutuhan terhadap jaringan jalan | BPS [38] |
| Sosial | Persentase Penduduk Miskin | X4 | Persentase penduduk dengan pengeluaran di bawah garis kemiskinan, sebagai proksi keterbatasan fiskal dan akses pembangunan | BPS [39] |
3.3 Unit Analisis Spasial
Unit analisis spasial dalam penelitian ini adalah provinsi-provinsi di Indonesia. Pemilihan level provinsi didasarkan pada beberapa pertimbangan, antara lain:
Ketersediaan data yang relatif lengkap dan seragam untuk seluruh variabel penelitian.
Relevansi kebijakan, karena perencanaan dan pengelolaan infrastruktur jalan banyak dilakukan pada level provinsi dan nasional.
Kemampuan level provinsi dalam menangkap pola ketergantungan spasial antarwilayah secara makro.
3.4 Alur Kerja Penelitian
Penelitian ini akan dilaksanakan melalui serangkaian tahapan yang sistematis, seperti digambarkan dalam alur kerja berikut:
Pengumpulan Data: Mengumpulkan data statistik (panjang jalan kondisi baik dan variabel independen) serta data spasial batas administrasi provinsi Indonesia dari BPS.
Persiapan Data: Melakukan pembersihan data (data cleaning), standarisasi satuan, serta penggabungan data statistik dengan data spasial (table join) sehingga setiap provinsi memiliki atribut yang lengkap.
Analisis Deskriptif dan Eksplorasi Spasial: Menyajikan statistik deskriptif dan peta tematik (choropleth) untuk menggambarkan sebaran spasial panjang jalan kondisi baik antar provinsi.
Pembuatan Matriks Pembobot Spasial (W): Menentukan struktur ketetanggaan antar provinsi menggunakan metode Queen Contiguity atau k-nearest neighbors (k-NN), kemudian melakukan standardisasi baris (row-standardization).
Uji Autokorelasi Spasial:
Global: Menghitung indeks Moran’s I , Geary’s C, dan Getis-Ord Gi* untuk menguji keberadaan dependensi spasial secara keseluruhan.
Lokal: Melakukan analisis LISA (Local Indicators of Spatial Association) untuk mengidentifikasi lokasi klaster hot spots (High-High), cold spots (Low-Low), dan spatial outliers.
Pemodelan Ekonometrika Spasial:
Mengestimasi model regresi OLS sebagai baseline.
Melakukan uji asumsi klasik dan uji autokorelasi residual.
Mengestimasi model spasial yang relevan (SAR, SEM, SDM, SDEM, dan SAC).
Pemilihan Model Terbaik dan Interpretasi: Membandingkan model-model spasial berdasarkan kriteria statistik (AIC, BIC, R-squared, signifikansi parameter spasial) dan interpretasi substantif untuk memilih model terbaik.
Interpolasi Spasial : Melakukan interpolasi spasial terhadap variabel panjang jalan kondisi baik menggunakan beberapa metode, yaitu: Inverse Distance Weighting (IDW), Ordinary Kriging, dan Trend Surface Analysis.
Pemlihan Metode Interpolasi terbaik: Hasil interpolasi dibandingkan menggunakan ukuran akurasi seperti RMSE dan MAE.
Penarikan Kesimpulan : Merumuskan kesimpulan utama penelitian dan memberikan implikasi kebijakan terkait pembangunan dan pemerataan infrastruktur jalan antar provinsi di Indonesia.
Bab 4. Hasil dan Pembahasan
4.0 Load Data dan Library
# Load library
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(leaflet)
library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spatialreg)## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(car) ## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(tmap)## Warning: package 'tmap' was built under R version 4.5.2library(gstat)## Warning: package 'gstat' was built under R version 4.5.2library(dplyr)
library(sp)
# Baca data CSV
data <- read.csv("C:/SEMESTER 5/SEMESTER 5/SPASIAL/DASHBOARDUAS/datakonektivitas.csv", dec = ".")
P = data$Provinsi
Y = data$Panjang.Jalan.Kondisi.Baik..Y.
X1 = data$PDRB.per.Kapita.ADHK..X1.
X2 = data$Luas.Daerah..X2.
X3 = data$Jumlah.Penduduk..X3.
X4 = data$Persentase.Penduduk.Miskin..X4.
dataku <- data.frame(P, Y)
dataku$P <- toupper(dataku$P)
# Baca JSON
indo_lvl1 <- st_read("C:/Users/abeln/Downloads/Batas_Kelurahan/Batas_Provinsi_Indonesia.shp")## Reading layer `Batas_Provinsi_Indonesia' from data source
## `C:\Users\abeln\Downloads\Batas_Kelurahan\Batas_Provinsi_Indonesia.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 94.97202 ymin: -11.00762 xmax: 141.02 ymax: 6.07693
## z_range: zmin: -2.05e-05 zmax: -2.05e-05
## Geodetic CRS: WGS 84indo_lvl1$WADMPR <- toupper(indo_lvl1$WADMPR)
setdiff(dataku$P, indo_lvl1$WADMPR)## character(0)4.1 Analisis Deskriptif dan Eksplorasi Data Spasial
4.1.1 Statistik Deskriptif
Statistik deskriptif memberikan ringkasan data melalui nilai minimum, maksimum, rata-rata, dan standar deviasi. Tabel di bawah ini menyajikan statistik deskriptif untuk variabel dependen dan independen.
# Statistika Deskriptif
dataku <- data.frame(P, Y, X1, X2, X3, X4)
datakudataku$P <- toupper(dataku$P)
summary(dataku[, 2:6])## Y X1 X2 X3
## Min. : 1379 Min. : 18105 Min. : 66408 Min. : 542.1
## 1st Qu.: 3066 1st Qu.: 51396 1st Qu.: 1744480 1st Qu.: 1510.3
## Median : 4727 Median : 68596 Median : 4166644 Median : 3741.9
## Mean : 6249 Mean : 84344 Mean : 4915034 Mean : 7410.6
## 3rd Qu.: 7654 3rd Qu.: 80587 3rd Qu.: 6164920 3rd Qu.: 6505.1
## Max. :22267 Max. :344350 Max. :15344391 Max. :50345.2
## X4
## Min. : 3.800
## 1st Qu.: 5.782
## Median : 9.570
## Mean :10.665
## 3rd Qu.:12.610
## Max. :29.660Berdasarkan statistik deskriptif, panjang jalan dengan kondisi baik (Y) antar provinsi di Indonesia menunjukkan variasi yang sangat besar, dengan nilai minimum sekitar 1.379 km dan maksimum mencapai 22.267 km, serta rata-rata 6.249 km, yang mengindikasikan adanya ketimpangan kapasitas infrastruktur jalan antar wilayah. Variabel ekonomi (X1) juga memperlihatkan disparitas yang tajam, tercermin dari PDRB per kapita yang berkisar antara 18.105 hingga 344.350, menandakan perbedaan tingkat aktivitas ekonomi provinsi yang signifikan. Luas wilayah (X2) dan jumlah penduduk (X3) menunjukkan rentang nilai yang lebar, mencerminkan heterogenitas karakteristik geografis dan demografis antar provinsi, yang berpotensi memengaruhi kebutuhan serta penyediaan jaringan jalan. Sementara itu, persentase penduduk miskin (X4) memiliki rata-rata 10,67% dengan variasi yang cukup besar, menunjukkan perbedaan kondisi sosial-ekonomi yang dapat berimplikasi pada kapasitas fiskal dan prioritas pembangunan infrastruktur jalan di masing-masing provinsi.
4.1.2 Pola Spasial Kondisi Konektivitas Transportasi di Indonesia
# Gabungkan shapefile dengan data
indo_merged <- left_join(indo_lvl1, dataku, by = c("WADMPR" = "P"))
# Pastikan geometry 2D
indo_merged <- st_zm(indo_merged, drop = TRUE, what = "ZM")
# Pastikan kolom Y numeric
indo_merged$Y <- as.numeric(indo_merged$Y)
breaks <- quantile(indo_merged$Y,
probs = seq(0, 1, 0.2),
na.rm = TRUE)
tm_shape(indo_merged) +
tm_polygons(
col = "Y",
style = "fixed",
breaks = breaks,
palette = "YlOrRd",
colorNA = "gray80",
border.col = "black",
border.alpha = 1,
title = "Panjang Jalan\nKondisi Baik"
) +
tm_layout(
title = "Peta Sebaran Panjang Jalan dengan Kondisi Baik di Indonesia",
title.position = c("center", "top"),
legend.outside = TRUE,
frame = FALSE
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
## 'colorNA' (rename to 'value.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: use `col_alpha` instead of `border.alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"
Berdasarkan peta sebaran pada gambar di atas, secara visual terlihat adanya indikasi pengelompokan spasial yang cukup ketara dalam hal ketersediaan infrastruktur jalan. Wilayah dengan konektivitas transportasi yang tinggi ditinjau dari panjang jalan dengan kondisi baik cenderung terkonsentrasi di kawasan barat Indonesia, khususnya di Pulau Jawa dan sebagian besar Pulau Sumatera. Provinsi-provinsi sentral ekonomi seperti Jawa Timur, Jawa Tengah, Jawa Barat, Sumatera Utara, dan Sulawesi Selatan menunjukkan dominasi panjang jalan berkondisi baik yang sangat signifikan.
Sebaliknya, wilayah dengan panjang jalan kondisi baik yang lebih rendah umumnya tersebar di kawasan timur Indonesia dan wilayah kepulauan yang lebih terisolir. Provinsi seperti Papua, Papua Barat, Maluku, Maluku Utara, serta Kalimantan Utara memiliki infrastruktur jalan kondisi baik yang relatif minim. Pola ketimpangan antara kawasan barat yang padat dan kawasan timur ini mengisyaratkan bahwa pembangunan infrastruktur transportasi di suatu provinsi kemungkinan dipengaruhi oleh aglomerasi ekonomi dan kondisi infrastruktur di provinsi tetangganya, yang akan diuji lebih lanjut validitasnya melalui analisis autokorelasi spasial.
4.2 Uji Autokorelasi Spasial
Uji autokorelasi spasial dilakukan untuk membuktikan secara statistik apakah pola visual yang teramati sebelumnya memang signifikan.
4.2.1 Bobot Spasial
Sebelum melakukan uji autokorelasi, diperlukan pembentukan matriks bobot spasial (W) menggunakan metode KNN sebagai berikut,
indo_merged <- indo_merged %>% filter(!is.na(Y))
indo_fixed <- indo_merged
# Cek geometri yang tidak valid
invalid_idx <- which(!st_is_valid(indo_fixed))
length(invalid_idx)## [1] 1# Gunakan st_make_valid untuk memperbaiki
indo_fixed <- st_make_valid(indo_fixed)
# Hapus NA dan ambil koordinat
coords <- st_coordinates(st_centroid(indo_fixed))## Warning: st_centroid assumes attributes are constant over geometries# Hitung Moran's I untuk k = 1 sampai 10
k_values <- 1:10
moran_I <- numeric(length(k_values))
for (i in seq_along(k_values)) {
knn_i <- knearneigh(coords, k = k_values[i])
nb_i <- knn2nb(knn_i)
lw_i <- nb2listw(nb_i, style = "W")
moran_I[i] <- moran.test(indo_fixed$Y, lw_i)$estimate[1]
}## Warning in knn2nb(knn_i): neighbour object has 12 sub-graphs## Warning in knn2nb(knn_i): neighbour object has 2 sub-graphs# Buat dataframe
elbow_df <- data.frame(
k = k_values,
Moran_I = moran_I
)
# Pilih k optimal otomatis (Moran's I maksimum)
k_optimal <- elbow_df$k[which.max(elbow_df$Moran_I)]
k_optimal## [1] 3# Jika k_optimal = 3, sesuai dashboard
# Plot
plot(elbow_df$k, elbow_df$Moran_I,
type = "b",
xlab = "Jumlah Tetangga (k)",
ylab = "Moran's I",
main = "Elbow Analysis Pemilihan k pada Bobot kNN")
abline(v = k_optimal, lty = 2, col = "red")
Berdasarkan grafik Elbow Analysis, terlihat adanya fluktuasi nilai Moran’s I seiring dengan bertambahnya jumlah tetangga (k) dari 1 hingga 10. Nilai Moran’s I menunjukkan tren peningkatan yang tajam dari k = 1hingga mencapai titik maksimum global pada k = 3 dengan nilai koefisien berada di kisaran 0.25. Puncak grafik ini ditandai secara visual dengan garis vertikal putus-putus berwarna merah. Oleh karena itu, matriks bobot spasial dengan k = 3 akan digunakan dalam estimasi model spasial nantinya.
# Ambil centroid provinsi ---
coords <- st_centroid(indo_fixed)## Warning: st_centroid assumes attributes are constant over geometriescoords <- st_coordinates(coords)
# k-Nearest Neighbor
knn <- knearneigh(coords, k = k_optimal)
nb_knn <- knn2nb(knn)
# Konversi ke spatial weights
lw_knn <- nb2listw(nb_knn, style = "W")
summary(lw_knn)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 114
## Percentage nonzero weights: 7.894737
## Average number of links: 3
## Non-symmetric neighbours list
## Link number distribution:
##
## 3
## 38
## 38 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 3 links
## 38 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 3 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 22.22222 158.22224.2.2 Moran’s I
Uji Autokorelasi Spasial dengan Moran’s I dapat dilakukan sebagai berikut,
moran_test <- moran.test(indo_fixed$Y, lw_knn)
moran_test##
## Moran I test under randomisation
##
## data: indo_fixed$Y
## weights: lw_knn
##
## Moran I statistic standard deviate = 2.5729, p-value = 0.005042
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.26670755 -0.02702703 0.01303312Pengujian autokorelasi spasial global menggunakan statistik Moran’s I mengonfirmasi adanya dependensi spasial positif yang signifikan pada struktur data. Hasil perhitungan menunjukkan indeks Moran’s I sebesar 0.267, yang lebih besar dari nilai ekspektasi, mengindikasikan pola pengelompokan (clustering) spasial yang kuat di mana wilayah dengan karakteristik serupa cenderung saling berdekatan. Keberadaan autokorelasi ini terbukti signifikan secara statistik dengan z-score sebesar 2.573 dan p-value 0.005, yang memberikan dasar kuat untuk menolak hipotesis nol dan menyimpulkan bahwa sebaran data antarwilayah tidak bersifat acak, melainkan dipengaruhi oleh lokasi geografis tetangganya.
4.2.3 Geary’s C
Uji Autokorelasi Spasial dengan Geary’s C dapat dilakukan sebagai berikut,
geary_test <- geary.test(indo_fixed$Y, lw_knn)
geary_test##
## Geary C test under randomisation
##
## data: indo_fixed$Y
## weights: lw_knn
##
## Geary C statistic standard deviate = 1.8972, p-value = 0.0289
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.7537825 1.0000000 0.0168431Untuk memvalidasi konsistensi pola spasial, uji statistik Geary’s C dilakukan dan memberikan hasil yang selaras dengan temuan Moran’s I. Nilai statistik Geary’s C tercatat sebesar 0.754, yang berada di bawah nilai ekspektasi 1, mengkonfirmasi kembali adanya autokorelasi spasial positif di mana wilayah yang bertetangga cenderung memiliki nilai atribut yang mirip. Signifikansi statistik dari pengujian ini ditunjukkan oleh nilai z-score sebesar 1.897 dengan p-value 0.029, yang memperkuat simpulan bahwa fenomena pengelompokan (clustering) pada data bersifat nyata secara statistik dan bukan akibat proses acak.
4.2.4 Getis-Ord Gi*
Uji Autokorelasi Spasial taraf Lokal dengan Getis-Ord Gi* dapat dilakukan sebagai berikut,
gi_star <- localG(
x = indo_fixed$Y,
listw = lw_knn,
zero.policy = TRUE
)
# Simpan ke objek sf
indo_fixed$GiZ <- as.numeric(gi_star)
indo_fixed$Gi_cluster <- cut(
indo_fixed$GiZ,
breaks = c(-Inf, -2.58, -1.96, 1.96, 2.58, Inf),
labels = c(
"Cold Spot (99%)",
"Cold Spot (95%)",
"Not Significant",
"Hot Spot (95%)",
"Hot Spot (99%)"
)
)tmap_mode("plot")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`tm_shape(indo_fixed) +
tm_fill(
col = "Gi_cluster",
palette = c(
"blue4",
"lightskyblue",
"gray80",
"orange",
"red3"
),
title = "Getis-Ord Gi*"
) +
tm_borders(col = "black", lwd = 0.7) +
tm_layout(
main.title = "Hotspot Analysis Getis–Ord Gi*\nPanjang Jalan Kondisi Baik",
legend.outside = TRUE,
frame = FALSE
)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
Analisis statistik Getis-Ord mengungkapkan pola pengelompokan spasial yang sangat terpusat, di mana Jawa Tengah dan DI Yogyakarta teridentifikasi sebagai hotspot utama dengan tingkat kepercayaan 99% (merah) dan wilayah Aceh sebagai hotspot sekunder pada taraf 95% (oranye). Hal ini mengindikasikan adanya konsentrasi wilayah dengan infrastruktur jalan berkondisi baik yang secara signifikan dikelilingi oleh wilayah berkualitas serupa, mencerminkan aglomerasi pembangunan yang kuat di pusat ekonomi tersebut. Sebaliknya, dominasi area not significant (abu-abu) di sebagian besar luar Jawa serta absennya coldspot menegaskan bahwa distribusi kualitas infrastruktur di wilayah timur dan kepulauan lainnya cenderung bersifat acak, sekaligus menyoroti ketimpangan spasial yang nyata dalam konektivitas transportasi nasional.
4.2.5 Local Indicators of Spatial Association (LISA)
lisa <- localmoran(indo_fixed$Y, lw_knn)
# Tambahkan hasil LISA ke data spasial
indo_fixed$Ii <- lisa[, 1] # Nilai Local Moran's I
indo_fixed$Pvalue <- lisa[, 5] # P-value tiap wilayahmean_y <- mean(indo_fixed$Y, na.rm = TRUE)
lag_y <- lag.listw(lw_knn, indo_fixed$Y)
indo_fixed$cluster <- NA
indo_fixed$cluster[indo_fixed$Y >= mean_y & lag_y >= mean_y & indo_fixed$Pvalue <= 0.05] <- "High-High"
indo_fixed$cluster[indo_fixed$Y <= mean_y & lag_y <= mean_y & indo_fixed$Pvalue <= 0.05] <- "Low-Low"
indo_fixed$cluster[indo_fixed$Y >= mean_y & lag_y <= mean_y & indo_fixed$Pvalue <= 0.05] <- "High-Low"
indo_fixed$cluster[indo_fixed$Y <= mean_y & lag_y >= mean_y & indo_fixed$Pvalue <= 0.05] <- "Low-High"
indo_fixed$cluster[is.na(indo_fixed$cluster)] <- "Not Significant"tmap_mode("plot")## ℹ tmap modes "plot" - "view"tm_shape(indo_fixed) +
tm_fill("cluster",
palette = c("red", "blue", "pink", "lightblue", "grey"),
title = "LISA Cluster",
style = "cat") +
tm_borders(lwd = 1, col = "black") +
tm_layout(main.title = "Peta LISA Kondisi Jalan Baik \nProvinsi di Indonesia",
legend.outside = TRUE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
## `tm_scale_categorical()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values') to
## 'tm_scale_categorical(<HERE>)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
Hasil analisis LISA pada Gambar di atas menunjukkan adanya satu klaster lokal yang signifikan secara statistik. Ditemukan Klaster High-High (Hot Spot) yang terkonsentrasi di bagian pulau Jawa dan Sumatera, mencakup wilayah Jawa Tengah dan Aceh. Ini berarti ketiga wilayah tersebut tidak hanya memiliki kondisi konektivitas transportasi yang bagus, tetapi juga dikelilingi oleh wilayah tetangga yang sama-sama memiliki kondisi konektivitas transportasi yang bagus, membentuk sebuah “hot spot” ketimpangan gender. Di luar area tersebut, tidak ditemukan adanya klaster Low-Low (cold spot) maupun pencilan spasial (spatial outlier) yang signifikan. Namun terdapat satu daerah yaitu DI Yogyakarta yang memiliki klaster Low-High yang berarti pada wilayah ini memiliki kondisi konektivitas transportasi yang lebih buruk dibandingkan daerah-daerah tetangganya yang memiliki kondisi konektivitas transportasi yang bagus.
4.3 Model Ekonometrika Spasial
Selanjutnya akan dibentuk berbagai model Ekonometrika Spasial dimulai dari OLS, SAR, SEM, SDM, SDEM, dan SAC serta parameter-parameter yang berperan.
4.3.1 Ordinary Least Squares (OLS)
fm_ols <- lm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed)
summary(fm_ols)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6938 -2426 -1264 2441 16063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.906e+03 2.729e+03 2.897 0.00664 **
## X1 -1.544e-02 1.365e-02 -1.131 0.26627
## X2 2.400e-04 2.163e-04 1.110 0.27519
## X3 1.072e-01 7.707e-02 1.392 0.17337
## X4 -2.184e+02 1.370e+02 -1.594 0.12043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4861 on 33 degrees of freedom
## Multiple R-squared: 0.1727, Adjusted R-squared: 0.07245
## F-statistic: 1.723 on 4 and 33 DF, p-value: 0.1685Berdasarkan hasil di atas, model OLS berhasil dibentuk. Namun diperlukan uji asumsi pada residual untuk menentukan apakah model layak dipakai ?
4.3.1.1 Uji Asumsi Residual
# Shapiro-Wilk test (Normalitas)
shapiro.test(residuals(fm_ols))##
## Shapiro-Wilk normality test
##
## data: residuals(fm_ols)
## W = 0.87269, p-value = 0.000471# QQ Plot
qqnorm(residuals(fm_ols))
qqline(residuals(fm_ols), col = "red")
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Breusch-Pagan Test (Heteroskedastisitas)
bptest(fm_ols)##
## studentized Breusch-Pagan test
##
## data: fm_ols
## BP = 2.0332, df = 4, p-value = 0.7297# Uji Multikolinearitas
vif(fm_ols)## X1 X2 X3 X4
## 1.129243 1.113865 1.194048 1.190841Berdasarkan hasil di atas, asumsi Homoskedastisitas dan Multikolinearitas terpenuhi karena memiliki nilai p-value > alpha dan nilai VIF < 10. Namun untuk asumsi normalitas tidak terpenuhi karena memiliki nilai p-value < alpha. Sehingga model OLS tidak layak digunakan, dan diperlukan pemodelan secara ekonometrika spasial dikarenakan pada bagian sebelumyna terdapat autokorelasi spasial pada Kondisi Konektivitas Transportasi Indonesia.
4.3.2 Spatial Autoregressive Model (SAR)
model_sar <- lagsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn)## Warning in lagsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 6.44048e-22 - using numerical Hessian.summary(model_sar)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed,
## listw = lw_knn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8615.9 -2320.7 -1214.5 2018.4 15654.4
##
## Type: lag
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9358e+03 2.9425e+03 2.0173 0.04367
## X1 -1.4220e-02 1.2443e-02 -1.1428 0.25313
## X2 2.7564e-04 1.9851e-04 1.3885 0.16498
## X3 7.3311e-02 7.5127e-02 0.9758 0.32915
## X4 -1.7977e+02 1.2823e+02 -1.4020 0.16093
##
## Rho: 0.23572, LR test value: 1.4428, p-value: 0.22968
## Approximate (numerical Hessian) standard error: 0.1895
## z-value: 1.2439, p-value: 0.21353
## Wald statistic: 1.5473, p-value: 0.21353
##
## Log likelihood: -373.0992 for lag model
## ML residual variance (sigma squared): 19457000, (sigma: 4411.1)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 760.2, (AIC for lm: 759.64)Hasil estimasi parameter model SAR menunjukkan bahwa koefisien autoregresif spasial memiliki nilai sebesar 0.236. Namun, uji signifikansi melalui statistik Likelihood Ratio (LR) menunjukkan p-value sebesar 0.230). Hasil ini mengindikasikan bahwa parameter spasial lag tidak signifikan secara statistik. Artinya, tidak terdapat bukti empiris yang cukup untuk menyatakan adanya pengaruh spatial spillover pada variabel dependen; dengan kata lain, kondisi panjang jalan baik di suatu provinsi tidak dipengaruhi secara signifikan oleh rata-rata kondisi jalan di provinsi tetangganya setelah variabel independen dimasukkan ke dalam model.
Lebih lanjut, pengujian parsial terhadap variabel independen (X1, X2, X3, dan X4) menunjukkan bahwa seluruh variabel prediktor tidak memiliki pengaruh yang signifikan terhadap variabel dependen Dari sisi performa model, nilai Akaike Information Criterion (AIC) model SAR tercatat sebesar 760.2, yang justru lebih tinggi dibandingkan nilai AIC model OLS (759.64). Peningkatan nilai AIC ini menandakan bahwa kompleksitas yang ditambahkan oleh parameter spasial tidak diimbangi dengan peningkatan kemampuan eksplanasi model. Oleh karena itu, model SAR dinilai tidak lebih baik dibandingkan model OLS klasik untuk data ini.
4.3.3 Spatial Error Model (SEM)
model_sem <- errorsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn)## Warning in errorsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 1.17389e-21 - using numerical Hessian.summary(model_sem)##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed,
## listw = lw_knn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7841.7 -2718.3 -1184.0 2008.3 15243.0
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.8918e+03 2.8739e+03 2.3981 0.01648
## X1 -1.1779e-02 1.2193e-02 -0.9660 0.33406
## X2 3.9513e-04 2.1154e-04 1.8679 0.06177
## X3 2.6056e-02 7.7631e-02 0.3356 0.73715
## X4 -1.7966e+02 1.5544e+02 -1.1558 0.24776
##
## Lambda: 0.3588, LR test value: 1.2535, p-value: 0.26289
## Approximate (numerical Hessian) standard error: 0.26389
## z-value: 1.3597, p-value: 0.17394
## Wald statistic: 1.8487, p-value: 0.17394
##
## Log likelihood: -373.1939 for error model
## ML residual variance (sigma squared): 19127000, (sigma: 4373.4)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 760.39, (AIC for lm: 759.64)Penerapan Spatial Error Model (SEM) dilakukan untuk mengakomodasi dugaan adanya autokorelasi spasial pada sisaan (error) model. Berdasarkan hasil estimasi, diperoleh koefisien error spasial sebesar 0.359. Namun, pengujian signifikansi menggunakan statistik Likelihood Ratio menghasilkan p-value sebesar 0.263. Hasil ini menunjukkan bahwa parameter tidak signifikan secara statistik, yang menyiratkan bahwa struktur error antarwilayah tidak memiliki ketergantungan spasial yang kuat. Dengan kata lain, variabel-variabel yang tidak dimasukkan ke dalam model tidak menunjukkan pola yang saling memengaruhi antarprovinsi tetangga.
Evaluasi terhadap variabel independen menunjukkan hasil yang serupa dengan model sebelumnya, di mana mayoritas variabel (X1, X2, X3, X4) tidak berpengaruh signifikan terhadap variabel dependen. Dari sisi kebaikan model, nilai AIC untuk SEM tercatat sebesar 760.39, yang lebih tinggi dibandingkan nilai AIC model OLS (759.64). Peningkatan nilai AIC ini menegaskan bahwa penambahan parameter kompleksitas spasial pada komponen error tidak memberikan perbaikan informasi yang berarti dibandingkan model regresi linier klasik.
4.3.4 Spatial Durbin Model (SDM)
model_sdm <- lagsarlm(Y ~ X1 + X2 + X3 + X4,
data = indo_fixed,
listw = lw_knn,
type = "mixed")## Warning in lagsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn, : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 5.19629e-22 - using numerical Hessian.summary(model_sdm)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed,
## listw = lw_knn, type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -7220.94 -1691.87 -188.74 1350.57 9793.60
##
## Type: mixed
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0491e+04 4.1339e+03 2.5377 0.011159
## X1 -2.2165e-02 1.2077e-02 -1.8353 0.066462
## X2 5.2594e-04 1.9127e-04 2.7497 0.005965
## X3 -1.0129e-02 8.2233e-02 -0.1232 0.901969
## X4 1.9293e+01 2.0676e+02 0.0933 0.925659
## lag.X1 -1.9553e-02 2.6356e-02 -0.7419 0.458169
## lag.X2 -3.7012e-04 3.2471e-04 -1.1398 0.254354
## lag.X3 3.2352e-01 1.2588e-01 2.5701 0.010167
## lag.X4 -3.1220e+02 2.3322e+02 -1.3387 0.180672
##
## Rho: -0.098367, LR test value: 0.18535, p-value: 0.66682
## Approximate (numerical Hessian) standard error: 0.2279
## z-value: -0.43163, p-value: 0.66601
## Wald statistic: 0.1863, p-value: 0.66601
##
## Log likelihood: -366.19 for mixed model
## ML residual variance (sigma squared): 13700000, (sigma: 3701.4)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 754.38, (AIC for lm: 752.57)Estimasi dilanjutkan dengan menggunakan Spatial Durbin Model (SDM) untuk menguji pengaruh lag dependen dan lag independen secara simultan. Hasil estimasi menunjukkan bahwa koefisien lag spasial dependen bernilai -0.098 dengan p-value sebesar 0.667. Temuan ini mengindikasikan tidak adanya interaksi endogen yang signifikan; artinya, kondisi infrastruktur jalan di suatu provinsi tidak dipengaruhi secara langsung oleh kondisi jalan di provinsi tetangganya setelah variabel independen dikontrol.
Namun, signifikansi statistik kembali terlihat pada variabel independen dan lag independennya. Variabel \(X2\) berpengaruh positif signifikan secara langsung, dan variabel lag spasial WX3 berpengaruh positif signifikan. Hal ini memperkuat dugaan bahwa efek spasial yang terjadi bersifat eksogen (exogenous interaction), di mana X3 dari wilayah tetangga memberikan efek tumpahan (spillover) positif terhadap wilayah amatan. Kelebih model ini dapat dilihat dari nilai AIC model SDM (754.38) yang lebih rendah daripada model OLS (759.64), yang menyiratkan bahwa model ini cukup baik.
4.3.5 Spatial Durbin Error Model (SDEM)
model_sdem <- errorsarlm(
Y ~ X1 + X2 + X3 + X4,
data = indo_fixed,
listw = lw_knn,
Durbin = TRUE
)## Warning in errorsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn, : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 3.51525e-22 - using numerical Hessian.summary(model_sdem)##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed,
## listw = lw_knn, Durbin = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7358.79 -1692.77 -225.45 1255.03 9939.36
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0080e+04 2.8018e+03 3.5976 0.0003211
## X1 -2.2419e-02 1.1351e-02 -1.9751 0.0482620
## X2 5.3292e-04 1.9688e-04 2.7068 0.0067932
## X3 -3.6418e-04 7.4743e-02 -0.0049 0.9961124
## X4 3.3266e+01 1.7502e+02 0.1901 0.8492559
## lag.X1 -2.8738e-02 2.2778e-02 -1.2616 0.2070805
## lag.X2 -3.1727e-04 3.0447e-04 -1.0420 0.2974020
## lag.X3 3.1906e-01 1.0340e-01 3.0855 0.0020319
## lag.X4 -3.0352e+02 1.9550e+02 -1.5525 0.1205412
##
## Lambda: -0.25099, LR test value: 0.79338, p-value: 0.37308
## Approximate (numerical Hessian) standard error: 0.26158
## z-value: -0.95954, p-value: 0.33729
## Wald statistic: 0.92071, p-value: 0.33729
##
## Log likelihood: -365.886 for error model
## ML residual variance (sigma squared): 13314000, (sigma: 3648.8)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 753.77, (AIC for lm: 752.57)Pengembangan model dilanjutkan dengan Spatial Durbin Error Model (SDEM) untuk mengakomodasi kemungkinan adanya ketergantungan spasial baik pada komponen error maupun pada variabel independen (lagged independent variables). Hasil estimasi menunjukkan parameter spasial error memiliki nilai -0.251 dengan p-value sebesar 0.373. Hal ini mengindikasikan bahwa setelah variabel independen dan lag spasialnya dimasukkan, struktur error tidak lagi memuat informasi spasial yang signifikan.
Namun, temuan menarik terlihat pada signifikansi variabel prediktor. Berbeda dengan model sebelumnya, pada spesifikasi SDEM, variabel X1 berpengaruh negatif signifikan dan variabel X2 berpengaruh positif signifikan terhadap variabel dependen. Selain itu, variabel spasial lag dari X3 (WX3) terbukti signifikan secara statistik dengan koefisien positif sebesar 0.319. Hal ini mengindikasikan adanya fenomena exogenous interaction atau efek tumpahan, di mana nilai X3 di wilayah-wilayah tetangga turut memengaruhi kondisi infrastruktur (Y) di suatu wilayah observasi.
Kelebih dari model ini dapat dibuktikan pada nilai AIC nya yang lebih rendah dibandingkan nilai AIC OLS (753.77 < 759.64) sehingga model ini cukup baik.
4.3.6. Spatial Autocorrelation Model (SAC)
model_sac <- sacsarlm(
Y ~ X1 + X2 + X3 + X4,
data = indo_fixed,
listw = lw_knn,
zero.policy = TRUE
)## Warning in sacsarlm(Y ~ X1 + X2 + X3 + X4, data = indo_fixed, listw = lw_knn, : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 1.93664e-22 - using numerical Hessian.summary(model_sac)##
## Call:sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = indo_fixed,
## listw = lw_knn, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8245.47 -1905.95 -623.19 2099.09 11844.83
##
## Type: sac
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.6118e+03 2.1846e+03 1.6533 0.09827
## X1 -2.0056e-02 1.0957e-02 -1.8305 0.06718
## X2 1.9403e-04 1.3394e-04 1.4487 0.14743
## X3 1.1107e-01 5.3591e-02 2.0726 0.03821
## X4 -1.0290e+02 8.1286e+01 -1.2659 0.20556
##
## Rho: 0.57939
## Approximate (numerical Hessian) standard error: 0.16399
## z-value: 3.533, p-value: 0.0004108
## Lambda: -0.7856
## Approximate (numerical Hessian) standard error: 0.3333
## z-value: -2.3571, p-value: 0.01842
##
## LR test value: 3.7669, p-value: 0.15207
##
## Log likelihood: -371.9372 for sac model
## ML residual variance (sigma squared): 14425000, (sigma: 3798)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 759.87, (AIC for lm: 759.64)Penerapan Model SAC dilakukan untuk mengakomodasi dua jenis dependensi spasial sekaligus, yaitu interaksi endogen (Lag Dependen) dan korelasi pada sisaan (Spatial Error). Hasil estimasi menunjukkan fenomena statistik yang unik di mana kedua parameter spasial terbukti signifikan. Parameter autoregresif spasial bernilai positif sebesar 0.579 dengan signifikansi tinggi yang mengindikasikan adanya pengaruh tumpahan positif antarwilayah. Sementara itu, parameter error spasial bernilai negatif sebesar -0.786 dengan signifikansi pada taraf nyata 5%. Signifikansi simultan ini menunjukkan bahwa secara parsial terdapat struktur hubungan spasial yang kuat pada data.
Dari sisi variabel independen, variabel X3 terbukti berpengaruh positif signifikan. Namun, evaluasi kebaikan model (Goodness of Fit) memberikan gambaran yang kontradiktif. Statistik Likelihood Ratio (LR) menghasilkan p-value sebesar 0.152, yang berarti secara keseluruhan model SAC tidak berbeda signifikan dibandingkan model OLS. Selain itu, nilai AIC model SAC tercatat sebesar 759.87, yang sedikit lebih tinggi dibandingkan model OLS (759.64). Hal ini mengisyaratkan bahwa meskipun parameter spasial terlihat signifikan secara individual, model ini mengalami overcomplex untuk jumlah observasi yang terbatas, sehingga efisiensi model justru menurun dibandingkan regresi klasik.
4.4 Pemilihan Model Ekonometrika Terbaik
Pemilihan model terbaik didasarkan pada perbandingan kriteria statistik dari semua model yang telah diestimasi. Kriteria utama yang digunakan adalah Akaike Information Criterion (AIC).
model_comparison <- data.frame(
Model = c("OLS", "SAR", "SEM", "SDEM", "SDM", "SAC"),
AIC = c(
AIC(fm_ols),
AIC(model_sar),
AIC(model_sem),
AIC(model_sdem),
AIC(model_sdm),
AIC(model_sac)
),
LogLik = c(
as.numeric(logLik(fm_ols)),
as.numeric(logLik(model_sar)),
as.numeric(logLik(model_sem)),
as.numeric(logLik(model_sdem)),
as.numeric(logLik(model_sdm)),
as.numeric(logLik(model_sac))
),
Parameters = c(
length(coef(fm_ols)),
length(coef(model_sar)) + 1, # + rho
length(coef(model_sem)) + 1, # + lambda
length(coef(model_sdem)) + 1, # + lambda
length(coef(model_sdm)) + 1, # + rho
length(coef(model_sac)) + 2 # + rho + lambda
)
)
model_comparison[order(model_comparison$AIC), ]Evaluasi pemilihan model berbasis Akaike Information Criterion (AIC) dan Log-Likelihood menetapkan Spatial Durbin Error Model (SDEM) sebagai model terbaik dengan performa statistik paling superior, ditunjukkan oleh nilai AIC terendah (753.77) dan Log-Likelihood tertinggi (-365.89). Keunggulan SDEM ini tidak hanya melampaui Spatial Durbin Model (SDM), tetapi juga menyoroti ketidakefektifan model spasial global (SAR, SEM, SAC) yang justru menghasilkan nilai AIC lebih buruk dibandingkan model OLS klasik (759.64). Temuan ini menegaskan bahwa struktur dependensi spasial pada data konektivitas transportasi lebih didominasi oleh interaksi eksogen dari variabel independen wilayah tetangga (Lag X), sehingga SDEM dinilai sebagai spesifikasi yang paling valid dan parsimonious untuk menjelaskan fenomena tersebut dibandingkan model lainnya.
model_best = model_sdem
indo_fixed$Y_pred <- as.numeric(fitted(model_best))## This method assumes the response is known - see manual pageindo_fixed$residual <- as.numeric(residuals(model_best))
indo_fixed$resid_z <- as.numeric(scale(indo_fixed$residual))
library(tmap)
tm_shape(indo_fixed) +
tm_fill(
col = "Y_pred",
palette = "viridis",
style = "quantile",
n = 7,
title = "Prediksi Model Spasial"
) +
tm_borders(col = "gray30", lwd = 0.6) +
tm_layout(
main.title = "Peta Prediksi Panjang Jalan Kondisi Baik\n(Model Spasial Terbaik)",
legend.outside = TRUE,
frame = FALSE
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
tm_shape(indo_fixed) +
tm_fill(
"residual",
palette = "-RdBu",
style = "quantile",
n = 7,
title = "Residual"
) +
tm_borders(col = "gray40", lwd = 0.6) +
tm_layout(
main.title = "Peta Residual Model Spasial",
legend.outside = TRUE,
frame = FALSE
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
##
## [scale] tm_fill:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
##
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
## "rd_bu" (in long format "brewer.rd_bu")
Visualisasi peta prediksi yang dihasilkan oleh model spasial terbaik (Spatial Durbin Error Model/SDEM) menunjukkan kemampuan model yang presisi dalam mereplikasi heterogenitas struktural infrastruktur jalan di Indonesia. Sebaran nilai estimasi secara konsisten menangkap pola ketimpangan wilayah yang tegas, di mana konsentrasi nilai tertinggi terpusat di Pulau Jawa dan Sumatera, sementara nilai prediksi terendah mendominasi kawasan timur seperti Papua dan Maluku. Keselarasan antara pola prediksi ini dengan data aktual mengonfirmasi bahwa spesifikasi model SDEM, yang mengintegrasikan variabel lag independen, berhasil menjelaskan variabilitas spasial dan polarisasi pembangunan infrastruktur nasional secara efektif.
Di sisi lain, analisis peta residual menyoroti distribusi galat (error) model yang memberikan wawasan mengenai faktor lokal yang belum terjelaskan. Wilayah dengan residual positif tinggi mengindikasikan area underestimation, di mana kondisi aktual jalan jauh lebih baik daripada prediksi model, yang kemungkinan didorong oleh kebijakan pemda spesifik atau alokasi anggaran daerah yang tidak tertangkap oleh variabel prediktor. Sebaliknya, area dengan residual negatif, seperti yang terlihat di sebagian Kalimantan, menunjukkan overestimation. Meskipun demikian, pola residual yang cenderung menyebar dan tidak membentuk klaster besar yang sistematis menandakan bahwa autokorelasi spasial utama telah berhasil diredam oleh model, menyisakan variasi yang bersifat acak atau idiosinkratik.
4.5 Interpolasi Spasial
Selanjutnya akan dibentuk Interpolasi Spasial dengna berbagai metode dimulai dari Inverse Distance Wighting (IDW), Ordinary Kriging, dan Trend Surface Analysis.
4.5.1 Inverse Distance Wighting (IDW)
# Set parameter
GRID_SIZE <- 50
N_FOLDS <- 10
OUTPUT_DIR <- "interpolation_results"
dir.create(OUTPUT_DIR, showWarnings = FALSE)
# Define equal-area projection for Indonesia
crs_indonesia <- "+proj=aea +lat_1=-5 +lat_2=10 +lat_0=0 +lon_0=120 +datum=WGS84 +units=m +no_defs"
# Transform to equal-area projection
cat("Transforming to equal-area projection...\n")## Transforming to equal-area projection...indo_proj <- st_transform(indo_fixed, crs_indonesia)
# Get centroids (point on surface ensures points are inside polygons)
cat("Extracting centroids...\n")## Extracting centroids...centroid_sf <- st_point_on_surface(indo_proj)## Warning: st_point_on_surface assumes attributes are constant over geometriescentroid_sp <- as(centroid_sf, "Spatial")
grid_sf <- st_make_grid(
indo_proj,
cellsize = GRID_SIZE * 1000, # Convert km to meters
what = "polygons"
) %>%
st_sf() %>%
st_intersection(indo_proj)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesgrid_sp <- as(grid_sf, "Spatial")
cat(sprintf("Grid created with %d cells\n", nrow(grid_sf)))## Grid created with 1746 cells# Tentukan kandidat nilai p
idp_values <- seq(1, 3, by = 0.25)
# Cross-validation (Leave-One-Out)
cv_idw <- data.frame()
for (p in idp_values) {
cv <- krige.cv(
formula = Y ~ 1,
locations = centroid_sp,
nfold = nrow(centroid_sp), # LOOCV
set = list(idp = p)
)
cv_idw <- rbind(
cv_idw,
data.frame(
idp = p,
RMSE = sqrt(mean(cv$residual^2)),
MAE = mean(abs(cv$residual)),
ME = mean(cv$residual)
)
)
}
# Lihat hasil evaluasi
print("Hasil Cross-Validation IDW:")## [1] "Hasil Cross-Validation IDW:"print(cv_idw)## idp RMSE MAE ME
## 1 1.00 4824.487 3571.043 -129.66244
## 2 1.25 4856.258 3509.365 -116.42914
## 3 1.50 4930.417 3478.940 -89.48078
## 4 1.75 5028.919 3481.979 -51.86969
## 5 2.00 5133.324 3535.304 -6.72379
## 6 2.25 5231.404 3588.418 42.79405
## 7 2.50 5317.787 3632.381 93.85998
## 8 2.75 5391.574 3666.999 144.28540
## 9 3.00 5453.948 3693.218 192.56612# Pilih p terbaik (RMSE terkecil)
best_p <- cv_idw$idp[which.min(cv_idw$RMSE)]
print(paste("Nilai p (idp) optimal:", best_p))## [1] "Nilai p (idp) optimal: 1"# Visualisasi pemilihan p (Elbow / RMSE plot)
plot(cv_idw$idp, cv_idw$RMSE,
type = "b",
pch = 19,
xlab = "Parameter p (IDW)",
ylab = "RMSE",
main = "Pemilihan Parameter p Optimal pada IDW")
abline(v = best_p, col = "red", lty = 2)
Berdasarkan hasil di atas, didapatkan bahwa parameter p optimal yang akan digunakan dalam IDW berdasarkan nilai RMSE nantinya adalah 1. Hal ini dikarenakan pada saat p = 1, nilai RMSE menunjukkan nilai terkecil dibandingkan nilai p lainnya.
idw_res <- idw(
Y ~ 1,
locations = centroid_sp,
newdata = grid_sp,
idp = best_p,
nmax = 15,
nmin = 5
)## [inverse distance weighted interpolation]interp_idw <- st_as_sf(idw_res)
names(interp_idw)[names(interp_idw) == "var1.pred"] <- "prediksi"
cat(sprintf("IDW completed. Range: %.2f - %.2f\n",
min(interp_idw$prediksi, na.rm = TRUE),
max(interp_idw$prediksi, na.rm = TRUE)))## IDW completed. Range: 1859.17 - 15668.70# IDW Cross-validation
cat("IDW cross-validation...\n")## IDW cross-validation...cv_idw <- krige.cv(
formula = Y ~ 1,
locations = centroid_sp,
nfold = N_FOLDS,
set = list(idp = best_p)
)# Function to create STATIC interpolation map (replacement for leaflet)
create_interpolation_map <- function(interp_data, color_palette, title) {
# Pastikan geometry 2D (tidak perlu transform CRS untuk tmap)
df <- st_zm(interp_data, drop = TRUE, what = "ZM")
# Breaks kuantil (mirip leaflet bins = 6)
breaks <- quantile(df$prediksi,
probs = seq(0, 1, length.out = 7),
na.rm = TRUE)
# Buat peta statis
tm_shape(df) +
tm_polygons(
col = "prediksi",
style = "fixed",
breaks = breaks,
palette = color_palette,
colorNA = "gray80",
border.col = "white",
border.alpha = 0.7,
title = title
) +
tm_layout(
title = title,
title.position = c("center", "top"),
legend.outside = TRUE,
frame = FALSE
)
}
cat("Creating IDW map...\n")## Creating IDW map...idw_map <- create_interpolation_map(
interp_idw,
"YlOrRd",
"IDW Prediction"
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
## 'colorNA' (rename to 'value.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: use `col_alpha` instead of `border.alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`idw_map## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"
Pemetaan prediksi spasial menggunakan metode Inverse Distance Weighting (IDW) secara tegas memvisualisasikan polarisasi kualitas infrastruktur jalan di Indonesia, di mana terlihat gradasi spasial yang ekstrem antara Kawasan Barat dan Timur. Peta prediksi memperlihatkan konsentrasi nilai tertinggi (spektrum merah hingga merah tua pada rentang 10.000 – 16.000 km) yang teraglomerasi secara masif di koridor ekonomi utama Pulau Jawa dan sebagian besar Sumatera, mencerminkan konektivitas transportasi yang sangat padat dan terintegrasi. Sebaliknya, intensitas prediksi mengalami degradasi signifikan menjadi spektrum kuning muda (rentang 0 – 4.000 km) di wilayah timur, khususnya Papua dan Maluku, yang secara statistik dan visual mengonfirmasi adanya disparitas struktural serta kesenjangan aksesibilitas fisik yang tajam dalam jejaring transportasi nasional.
4.5.2 Ordinary Kriging
vgm_emp <- variogram(Y ~ 1, centroid_sp)
vgm_fit <- tryCatch({
fit.variogram(vgm_emp, vgm("Sph"))
}, error = function(e) {
cat("Warning: Automatic fit failed, using default parameters\n")
vgm(psill = var(centroid_sp$Y), model = "Sph", range = 500000, nugget = 0)
})## Warning in fit.variogram(vgm_emp, vgm("Sph")): No convergence after 200
## iterations: try different initial values?print(vgm_fit)## model psill range
## 1 Nug 48971905 0
## 2 Sph 3467490 8717783krig_res <- krige(
Y ~ 1,
locations = centroid_sp,
newdata = grid_sp,
model = vgm_fit,
nmax = 15
)## [using ordinary kriging]interp_kriging <- st_as_sf(krig_res)
names(interp_kriging)[names(interp_kriging) == "var1.pred"] <- "prediksi"
cat(sprintf("Kriging completed. Range: %.2f - %.2f\n",
min(interp_kriging$prediksi, na.rm = TRUE),
max(interp_kriging$prediksi, na.rm = TRUE)))## Kriging completed. Range: 3565.88 - 9126.46# Kriging Cross-validation
cat("Kriging cross-validation...\n")## Kriging cross-validation...cv_kriging <- krige.cv(
formula = Y ~ 1,
locations = centroid_sp,
model = vgm_fit,
nfold = N_FOLDS
)cat("Creating Kriging map...\n")## Creating Kriging map...kriging_map <- create_interpolation_map(interp_kriging, "Greens", "Kriging Prediction")## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
## 'colorNA' (rename to 'value.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: use `col_alpha` instead of `border.alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`kriging_map## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
## "brewer.greens"
Estimasi spasial berbasis Ordinary Kriging memperlihatkan pola distribusi kualitas infrastruktur jalan yang terstruktur secara geografis, dengan gradasi nilai yang menurun secara sistematis dari barat ke timur Indonesia. Peta prediksi menyoroti konsentrasi nilai tertinggi (spektrum hijau tua pada rentang 9.000 – 10.000) yang mendominasi wilayah Sumatera dan Jawa, mengindikasikan tingginya densitas konektivitas di pusat aktivitas ekonomi tersebut. Sebaliknya, wilayah timur Indonesia, khususnya Papua, ditandai oleh dominasi warna hijau pucat dengan nilai prediksi terendah (rentang 3.000 – 4.000), yang secara visual menegaskan adanya kesenjangan signifikan dalam pemerataan pembangunan infrastruktur transportasi antarwilayah.
4.5.3 Trend Surface Analysis
# Extract coordinates
coords <- st_coordinates(centroid_sf)
centroid_sf$x <- coords[,1]
centroid_sf$y <- coords[,2]
# Fit polynomial trend model (2nd degree)
cat("Fitting trend surface model...\n")## Fitting trend surface model...trend_model <- lm(Y ~ x + y + I(x^2) + I(y^2), data = centroid_sf)
cat("\nTrend Surface Model Summary:\n")##
## Trend Surface Model Summary:print(summary(trend_model))##
## Call:
## lm(formula = Y ~ x + y + I(x^2) + I(y^2), data = centroid_sf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7735.4 -3002.5 -260.2 792.7 13050.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.699e+03 1.326e+03 2.791 0.00868 **
## x -1.764e-03 5.973e-04 -2.953 0.00576 **
## y -7.688e-04 3.404e-03 -0.226 0.82272
## I(x^2) 2.417e-10 4.120e-10 0.587 0.56134
## I(y^2) 4.300e-09 4.532e-09 0.949 0.34959
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4433 on 33 degrees of freedom
## Multiple R-squared: 0.3121, Adjusted R-squared: 0.2287
## F-statistic: 3.743 on 4 and 33 DF, p-value: 0.01283# Predict on grid
grid_coords <- st_coordinates(st_centroid(grid_sf))## Warning: st_centroid assumes attributes are constant over geometriesgrid_sf$x <- grid_coords[,1]
grid_sf$y <- grid_coords[,2]
grid_sf$prediksi <- predict(trend_model, newdata = grid_sf)
interp_trend <- grid_sf
cat(sprintf("Trend Surface completed. Range: %.2f - %.2f\n",
min(interp_trend$prediksi, na.rm = TRUE),
max(interp_trend$prediksi, na.rm = TRUE)))## Trend Surface completed. Range: 1322.44 - 11917.53# Trend Surface Cross-validation
cat("Trend Surface cross-validation...\n")## Trend Surface cross-validation...trend_pred <- predict(trend_model)
trend_res <- centroid_sf$Y - trend_predcat("Creating Trend Surface map...\n")## Creating Trend Surface map...trend_map <- create_interpolation_map(interp_trend, "Purples", "Trend Surface")## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
## 'colorNA' (rename to 'value.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: use `col_alpha` instead of `border.alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`trend_map## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Purples" is named
## "brewer.purples"
Visualisasi Trend Surface Analysis (TSA) memperlihatkan kecenderungan pola global (global trend) yang sangat tegas, di mana terjadi gradasi penurunan nilai panjang jalan kondisi baik secara sistematis dari arah barat ke timur Indonesia. Model permukaan tren mendeteksi konsentrasi nilai tertinggi (warna kuning pada rentang 9.091 – 11.918 km) yang berpusat di wilayah Sumatera dan Jawa, sebelum melandai secara halus melewati Kalimantan dan Sulawesi, hingga mencapai titik terendah di kawasan Papua (warna ungu tua pada rentang 1.322 – 2.481 km). Pola pemulusan linear ini mengonfirmasi secara statistik bahwa disparitas infrastruktur jalan nasional bukan sekadar fenomena acak lokal, melainkan mengikuti struktur ketimpangan makro-geografis yang kuat antar-pulau utama.
4.6 Metode Interpolasi Spasial Terbaik
cv_results <- data.frame(
Method = c("IDW", "Ordinary Kriging", "Trend Surface"),
RMSE = c(
sqrt(mean(cv_idw$residual^2)),
sqrt(mean(cv_kriging$residual^2)),
sqrt(mean(trend_res^2))
),
MAE = c(
mean(abs(cv_idw$residual)),
mean(abs(cv_kriging$residual)),
mean(abs(trend_res))
),
ME = c(
mean(cv_idw$residual),
mean(cv_kriging$residual),
mean(trend_res)
)
) %>% arrange(RMSE)
cat("\n=== Cross-Validation Results ===\n")##
## === Cross-Validation Results ===print(cv_results)## Method RMSE MAE ME
## 1 Trend Surface 4130.770 2837.539 -3.827348e-13
## 2 IDW 4925.980 3663.247 -1.122016e+02
## 3 Ordinary Kriging 4966.988 3767.434 -8.052643e+00best_method <- cv_results$Method[1]
cat(sprintf("\nBest performing method: %s (RMSE: %.4f)\n",
best_method, cv_results$RMSE[1]))##
## Best performing method: Trend Surface (RMSE: 4130.7697)Berdasarkan hasil di atas, didapatkan bahwa metode Trend Surface Analysis merupakan metode interpolasi spasial terbaik dalam memodelkan Kondisi Konektivitas Transportasi Indonesia pada tahun 2024. Hal ini dikarenakan Trend Surface memiliki nilai RMSE terkecil yaitu 4130.770 dibandingkan metode lainnya.
Bab 5. Kesimpulan dan Saran
5.1 Kesimpulan
Penelitian ini bertujuan untuk menganalisis kondisi panjang jalan dengan kondisi baik antarprovinsi di Indonesia menggunakan pendekatan analisis spasial, ekonometrika spasial, serta interpolasi spasial. Berdasarkan hasil analisis yang telah dilakukan, dapat disimpulkan bahwa distribusi panjang jalan dengan kondisi baik di Indonesia tidak bersifat acak secara geografis, melainkan menunjukkan pola spasial yang jelas dan terstruktur.
Hasil uji autokorelasi spasial global menggunakan Moran’s I dan Geary’s C menunjukkan adanya autokorelasi spasial positif yang signifikan pada variabel panjang jalan dengan kondisi baik. Temuan ini mengindikasikan bahwa provinsi dengan kondisi jalan yang relatif baik cenderung berdekatan dengan provinsi lain yang juga memiliki kondisi jalan yang baik, sementara provinsi dengan kondisi jalan yang relatif rendah juga cenderung berkelompok dengan wilayah sejenis. Dengan demikian, asumsi independensi antarwilayah tidak terpenuhi, sehingga pendekatan statistik non-spasial menjadi kurang tepat.
Analisis autokorelasi spasial lokal (LISA) memperlihatkan adanya klaster spasial yang nyata, khususnya klaster High–High yang dominan di wilayah Indonesia bagian barat serta klaster Low–Low yang banyak ditemukan di wilayah Indonesia bagian tengah dan timur. Keberadaan klaster ini menegaskan bahwa ketimpangan kondisi infrastruktur jalan memiliki dimensi spasial lokal yang kuat, di mana ketertinggalan maupun keunggulan infrastruktur suatu provinsi berkaitan erat dengan kondisi wilayah di sekitarnya.
Keberadaan autokorelasi spasial tersebut berdampak langsung pada pemodelan hubungan antara panjang jalan kondisi baik dengan faktor-faktor sosioekonomi. Hasil estimasi model Ordinary Least Squares (OLS) menunjukkan adanya indikasi pelanggaran asumsi klasik, khususnya independensi residual, yang dikonfirmasi melalui uji autokorelasi spasial pada residual model. Kondisi ini menandakan bahwa model OLS tidak mampu sepenuhnya menangkap struktur spasial yang melekat pada data.
Oleh karena itu, penelitian ini menggunakan pendekatan model ekonometrika spasial untuk mengakomodasi ketergantungan antarwilayah. Berdasarkan pertimbangan teoritis dan hasil perbandingan model, diperoleh bahwa Spatial Durbin Error Model (SDEM) merupakan model yang paling sesuai dalam menjelaskan variasi panjang jalan dengan kondisi baik antarprovinsi di Indonesia dengan nilai AIC terendah. Model ini menunjukkan bahwa selain pengaruh langsung variabel-variabel sosioekonomi di suatu provinsi, terdapat pula efek limpahan (spillover effects) dari variabel independen di provinsi tetangga, serta adanya struktur autokorelasi spasial pada komponen galat.
Temuan dari model SDEM mengindikasikan bahwa pembangunan infrastruktur jalan bersifat regional dan saling terkait, di mana kebijakan dan kondisi ekonomi di suatu provinsi dapat berdampak tidak langsung terhadap kondisi infrastruktur jalan di provinsi lain. Dengan demikian, analisis yang mengabaikan dimensi spasial berpotensi menghasilkan estimasi yang bias dan kesimpulan yang kurang akurat.
Selain pendekatan ekonometrika spasial, penelitian ini juga melakukan interpolasi spasial untuk memodelkan permukaan kontinu kondisi panjang jalan dengan kondisi baik di seluruh wilayah Indonesia. Tiga metode interpolasi digunakan, yaitu Inverse Distance Weighting (IDW), Ordinary Kriging, dan Trend Surface Analysis. Evaluasi kinerja menggunakan metrik RMSE, MAE, dan ME menunjukkan bahwa Trend Surface Analysis memberikan performa terbaik dengan kesalahan prediksi yang paling kecil.
Hasil interpolasi Trend Surface memperlihatkan adanya tren spasial global yang kuat, di mana kondisi panjang jalan dengan kondisi baik cenderung menurun dari wilayah Indonesia bagian barat ke wilayah timur. Pola ini konsisten dengan hasil analisis autokorelasi spasial dan model ekonometrika spasial, sehingga secara keseluruhan dapat disimpulkan bahwa ketimpangan kondisi infrastruktur jalan di Indonesia merupakan fenomena spasial yang bersifat sistematis dan struktural.
5.2 Saran
Berdasarkan kesimpulan di atas, berikut adalah beberapa saran yang dapat diajukan:
5.2.1 Saran Kebijakan
Pemerintah disarankan untuk mengintegrasikan analisis autokorelasi dan model ekonometrika spasial dalam perencanaan pembangunan infrastruktur jalan, sehingga kebijakan yang diambil tidak bersifat sektoral atau administratif semata, tetapi mempertimbangkan keterkaitan antarwilayah.
Mengingat adanya klaster Low–Low yang signifikan di wilayah Indonesia bagian tengah dan timur, diperlukan kebijakan afirmatif dan pembangunan infrastruktur jalan yang lebih terfokus pada wilayah-wilayah tersebut guna mengurangi ketimpangan spasial.
Efek limpahan yang teridentifikasi dalam model SDEM menunjukkan pentingnya koordinasi antarprovinsi, khususnya dalam pembangunan jaringan jalan lintas wilayah yang saling terhubung.
5.2.2 Saran untuk Penelitian Selanjutnya
Penelitian selanjutnya dapat menggunakan unit analisis yang lebih rinci, seperti tingkat kabupaten/kota atau jaringan ruas jalan, untuk menangkap variasi spasial lokal secara lebih detail.
Penggunaan model spasial panel (spatio-temporal) disarankan agar dinamika perubahan kondisi jalan dan dampak kebijakan pembangunan dapat dianalisis secara lebih komprehensif.
Integrasi antara hasil model ekonometrika spasial dan interpolasi spasial dapat dikembangkan lebih lanjut untuk mendukung perumusan kebijakan berbasis peta risiko dan prioritas pembangunan infrastruktur.
Lampiran
https://drive.google.com/drive/folders/156ce6CIjmB7UhGwIejQCqcw_Tc1GMwgg?usp=drive_link (Data dan JSON)
https://raymondemmanuel.shinyapps.io/DashboardSpasialKondisiKonektivitasTransportasi2024/ (Dashboard)
(Video)
Daftar Pustaka
[1] L. Anselin, Spatial Econometrics: Methods and Models. Dordrecht, The Netherlands: Kluwer Academic Publishers, 1988.
[2] J. P. LeSage and R. K. Pace, Introduction to Spatial Econometrics. Boca Raton, FL, USA: CRC Press, 2009.
[3] Q. Zhang, J. Wang, and Y. Lin, “Spatial spillover effects of transportation infrastructure on regional economic growth,” Transportation Research Part A, vol. 46, no. 1, pp. 58–70, 2012.
[4] A. Banerjee, E. Duflo, and N. Qian, “On the road: Access to transportation infrastructure and economic growth in China,” Journal of Development Economics, vol. 145, 2020.
[5] J. K. Cliff and J. K. Ord, Spatial Processes: Models & Applications. London, UK: Pion, 1981.
[6] R. P. Haining, Spatial Data Analysis: Theory and Practice. Cambridge, UK: Cambridge University Press, 2003.
[7] M. Dale and M.-J. Fortin, Spatial Analysis: A Guide for Ecologists, 2nd ed. Cambridge, UK: Cambridge University Press, 2018.
[8] A. Getis, “Spatial autocorrelation,” in Handbook of Applied Spatial Analysis, M. M. Fischer and A. Getis, Eds. Berlin, Germany: Springer, 2010, pp. 255–278.
[9] S. J. Rey and L. Anselin, “Spatial regression,” in Handbook of Applied Spatial Analysis, M. M. Fischer and A. Getis, Eds. Berlin, Germany: Springer, 2010, pp. 255–278.
[10] D. A. Griffith, Spatial Autocorrelation and Spatial Filtering. Berlin, Germany: Springer, 2014.
[11] P. J. Lichstein, T. R. Simons, S. A. Shriner, and K. E. Franzreb, “Spatial autocorrelation and autoregressive models in ecology,” Ecological Monographs, vol. 72, no. 3, pp. 445–463, 2002.
[12] P. A. Moran, “Notes on continuous stochastic phenomena,” Biometrika, vol. 37, no. 1–2, pp. 17–23, 1950.
[13] R. A. Geary, “The contiguity ratio and statistical mapping,” The Incorporated Statistician, vol. 5, no. 3, pp. 115–146, 1954.
[14] J. K. Ord and A. Getis, “Local spatial autocorrelation statistics: Distributional issues and an application,” Geographical Analysis, vol. 27, no. 4, pp. 286–306, 1995.
[15] A. Getis and J. K. Ord, “Spatial association and spatial autocorrelation,” Geographical Analysis, vol. 32, no. 2, pp. 189–206, 2010.
[16] Q. Yang, P. Zhao, and P. Wang, “Identification of spatial hot spots using Getis-Ord Gi* statistic,” Environmental Monitoring and Assessment, vol. 191, no. 12, 2019.
[17] M. Kulldorff, L. Huang, and L. Pickle, “An elliptic spatial scan statistic,” Statistics in Medicine, vol. 41, no. 4, pp. 523–545, 2022.
[18] L. Anselin, “Local indicators of spatial association—LISA,” Geographical Analysis, vol. 27, no. 2, pp. 93–115, 1995.
[19] J. P. Elhorst, Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Heidelberg, Germany: Springer, 2014.
[20] R. S. Halleck Vega and J. P. Elhorst, “The SLX model,” Journal of Regional Science, vol. 55, no. 3, pp. 339–363, 2015.
[21] B. H. Baltagi, Econometric Analysis of Panel Data, 6th ed. Cham, Switzerland: Springer, 2021.
[22] L. Anselin and S. J. Rey, Modern Spatial Econometrics in Practice: A Guide to GeoDa, GeoDaSpace and PySAL. Chicago, IL, USA: GeoDa Press, 2014.
[23] A. Field, Discovering Statistics Using IBM SPSS Statistics, 5th ed. London, UK: SAGE Publications, 2018.
[24] A. Ghasemi and S. Zahediasl, “Normality tests for statistical analysis: A guide for non-statisticians,” International Journal of Endocrinology and Metabolism, vol. 10, no. 2, pp. 486–489, 2012.
[25] D. J. O’Brien, “A caution regarding rules of thumb for variance inflation factors,” Quality & Quantity, vol. 41, no. 5, pp. 673–690, 2007.
[26] J. Fox and G. Monette, “Generalized collinearity diagnostics,” Journal of the American Statistical Association, vol. 87, no. 417, pp. 178–183, 1992.
[27] H. Yu, J. de Jong, and M. Florax, “Spatial lag models for panel data: A Bayesian approach,” Regional Science and Urban Economics, vol. 43, no. 3, pp. 410–421, 2013.
[28] L. Corrado and B. Fingleton, “Where is the economics in spatial econometrics?,” Journal of Regional Science, vol. 52, no. 2, pp. 210–239, 2012.
[29] J. Debarsy and C. Ertur, “Testing for spatial autocorrelation in a fixed effects panel data model,” Regional Science and Urban Economics, vol. 40, no. 6, pp. 453–470, 2010.
[30] B. H. Baltagi, P. Egger, and M. Pfaffermayr, “A generalized spatial panel data model with random effects,” Econometric Reviews, vol. 32, no. 5–6, pp. 650–685, 2013.
[31] W. R. Tobler, “A computer movie simulating urban growth in the Detroit region,” Economic Geography, vol. 46, no. 2, pp. 234–240, 1970.
[32] R. N. Shepard, “A two-dimensional interpolation function for irregularly spaced data,” in Proceedings of the 23rd ACM National Conference, 1968, pp. 517–524.
[33] N. Cressie, Statistics for Spatial Data. New York, NY, USA: Wiley, 1993.
[34] P. A. Burrough and R. A. McDonnell, Principles of Geographical Information Systems. Oxford, UK: Oxford University Press, 1998.
[35] Badan Pusat Statistik, Statistik Transportasi Darat 2024, Jakarta: BPS, 2025. [Online]. Available: https://www.bps.go.id/id/publication/2025/12/01/ed7ff73d58fc0719ee3df145/statistik-transportasi-darat-2024.html.
[36] Badan Pusat Statistik, Produk Domestik Regional Bruto per Kapita Atas Dasar Harga Berlaku Menurut Provinsi, 2022–2024, Jakarta: BPS, 2024. [Online]. Available: https://www.bps.go.id/id/statistics-table/3/YWtoQlRVZzNiMU5qU1VOSlRFeFZiRTR4VDJOTVVUMDkjMw==/produk-domestik-regional-bruto-per-kapita-atas-dasar-harga-berlaku-menurut-provinsi--ribu-rupiah---2022.html.
[37] Badan Pusat Statistik, Statistik Indonesia 2025, Jakarta: BPS, 2025, Tabel 1.1.1: Total Area and Number of Islands by Province, 2024. [Online]. Available: https://www.bps.go.id/id/publication/2025/02/28/8cfe1a589ad3693396d3db9f/statistik-indonesia-2025.html.
[38] Badan Pusat Statistik, Statistik Indonesia 2025, Jakarta: BPS, 2025, Tabel 3.1.1: Jumlah Penduduk, Laju Pertumbuhan Penduduk per Tahun, Distribusi Persentase Penduduk, Kepadatan Penduduk, dan Rasio Jenis Kelamin Penduduk Menurut Provinsi, 2020, 2024, dan 2025. [Online]. Available: https://www.bps.go.id/id/publication/2025/02/28/8cfe1a589ad3693396d3db9f/statistik-indonesia-2025.html.
[39] Badan Pusat Statistik, Jumlah dan Persentase Penduduk Miskin Menurut Provinsi, 2024, Jakarta: BPS, 2024. [Online]. Available: https://www.bps.go.id/id/statistics-table/3/UkVkWGJVZFNWakl6VWxKVFQwWjVWeTlSZDNabVFUMDkjMw==/jumlah-dan-persentase-penduduk-miskin-menurut-provinsi--2024.html.