Analisis Spasial Determinan PDRB Pertanian, Perburuan, dan Peternakan di Indonesia Tahun 2023: Peran Komoditas Unggul dan Tenaga Kerja
Disusun oleh: Ariq Abdillah Syahlan (140610230065)
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
Pembangunan ekonomi di tingkat wilayah merupakan salah satu indikator penting dalam memahami dinamika pertumbuhan dan ketimpangan antar daerah[1]. Produk Domestik Regional Bruto (PDRB) sektor pertanian, perburuan, dan peternakan tidak hanya mencerminkan kapasitas produksi suatu provinsi, tetapi juga hubungan ekonomi antar wilayah yang mungkin memiliki keterikatan secara spasial[2]. Hal ini mungkin terjadi di Indonesia karena banyaknya keragaman geografis di tiap pulau maupun provinsi, sehingga nilai PDRB terutama sektor pertanian, perburuan, dan peternakan dapat dipengaruhi kondisi geografis dan juga kondisi provinsi tetangga seperti faktor-faktor integrasi pasar, mobilitas tenaga kerja, dan penyebaran teknologi[3].
Penelitian ini menggunakan variabel-variabel yang mencerminkan karakteristik ekonomi spasial sektor pertanian dan peternakan untuk menjelaskan variasi Produk Domestik Regional Bruto (PDRB) subsektor pertanian, perburuan, dan peternakan di tingkat provinsi tahun 2023. Variabel dependen dalam penelitian adalah PDRB subsektor pertanian, perburuan, dan peternakan. Variabel independen yang dipilih mencakup jumlah komoditas unggul tanaman pangan, jumlah komoditas unggul peternakan, dan jumlah tenaga kerja di bidang pertanian berumur di atas 15 tahun. Pemilihan komoditas unggul tanaman pangan dan peternakan dilakukan guna dapat menunjukkan pengaruh subsektor pertanian dalam hal ini tanaman pangan dan peternakan terhadap PDRB[4]. Lalu untuk ketenagakerjaan mencerminkan alokasi sumber daya manusia dalam sektor pertanian yang memiliki hubungan kuat dengan nilai tambah ekonomi regional[5].
Dalam penelitian yang menggunakan data wilayah perti PDRB per provinsi, penting untuk memahami sejauh mana suatu fenomena ekonomi menunjukkan keterkaitan antar wilayah secara spasial. Konsep autokorelasi spasial merujuk pada kecenderungan nilai suatu variabel untuk menunjukkan pola yang serupa atau berbeda di wilayah yang berdekatan, pola ini tidak dapat ditangkap oleh model statistik konvensional yang mengasumsikan independensi observasi[6]. Metode yang umum digunakan untuk mengukur autokorelasi spasial adalah Moran’s I, Geary’s C, dan juga LISA. Penggunaan metode tersebut bertujuan untuk mengidentifikasi pola pengelompokan wilayah (spatial clustering) serta mengetahui apakah nilai PDRB suatu provinsi dipengaruhi oleh nilai PDRB provinsi di sekitarnya[7].
Jika terdapat autokorelasi spasial maka perlu dilakukan analisis menggunakan model ekonometrik spasial. Model Ordinary Least Squares (OLS) digunakan sebagai model awal (baseline) untuk mengevaluasi asumsi independensi residual, namun sering kali tidak memadai ketika terdapat pengaruh spasial antar wilayah[8]. Oleh karena itu, digunakan model Spatial Autoregressive (SAR) yang menangkap pengaruh langsung nilai variabel dependen dari wilayah tetangga, serta Spatial Error Model (SEM) yang mengakomodasi ketergantungan spasial pada komponen error akibat faktor tak teramati[9]. Selain itu, model yang lebih fleksibel seperti Spatial Durbin Model (SDM) dan Spatial Durbin Error Model (SDEM) dipertimbangkan karena mampu menangkap efek limpahan spasial (spatial spillover) baik dari variabel dependen maupun variabel independen[10][11]. Penggunaan beberapa spesifikasi model ini memungkinkan pemilihan model terbaik yang sesuai dengan karakteristik data PDRB antar provinsi serta struktur ketergantungan spasial yang mendasarinya.
Berdasarkan uraian sebelumnya, permasalahan utama dalam penelitian ini adalah adanya potensi ketergantungan spasial pada PDRB subsektor pertanian, perburuan, dan peternakan antar provinsi yang belum sepenuhnya dapat dijelaskan oleh faktor internal wilayah saja. Selain itu, perlu diketahui juga sejauh mana komoditas unggul tanaman pangan, komoditas unggul peternakan, dan tenaga kerja pertanian usia di atas 15 tahun berkontribusi secara langsung maupun melalui efek limpahan spasial terhadap kinerja ekonomi subsektor pertanian, perburuan, dan peternakan. Oleh karena itu, tujuan penelitian ini adalah untuk mengidentifikasi pola autokorelasi spasial PDRB subsektor pertanian serta menganalisis pengaruh variabel-variabel penjelas menggunakan pendekatan model spasial ekonometrik yang sesuai. Adapun batasan penelitian ini meliputi penggunaan data penampang lintang (cross-section) pada tingkat provinsi di Indonesia tahun 2023, pembatasan variabel independen pada aspek komoditas unggul dan tenaga kerja pertanian, serta penggunaan matriks pembobot spasial berbasis kedekatan wilayah administratif, sehingga hasil penelitian diharapkan relevan dalam konteks analisis ekonomi regional berbasis spasial namun tidak dimaksudkan untuk inferensi kausal jangka panjang.
Bab 2. Tinjauan Pustaka
Penelitian ini difokuskan pada kajian empiris dan teoritis yang membahas determinan Produk Domestik Regional Bruto (PDRB) sektor pertanian, perburuan, dan peternakan serta pendekatan spasial yang digunakan dalam menganalisis fenomena ekonomi regional. Penelitian terdahulu menemukan bahwa peningkatan PDRB sektor pertanian di Provinsi Aceh berkontribusi langsung terhadap peningkatan PDRB provinsi secara keseluruhan, di mana subsektor kehutanan memberikan pengaruh terbesar[12]. Namun, penelitian tersebut masih menggunakan pendekatan konvensional tanpa mempertimbangkan keterkaitan antarwilayah.
Terdapat penelitian lainnya yang menganalisis pertumbuhan ekonomi di 119 kabupaten/kota di Pulau Jawa menggunakan Spatial Durbin Model (SDM)[13]. Hasil penelitian menunjukkan adanya spatial clustering dengan nilai Moran’s I positif serta adanya efek limpahan (spillover effect) dari faktor pendapatan awal dan pendidikan terhadap wilayah sekitarnya. Lalu ada penelitian yang melakukan penggabungan Stochastic Frontier Analysis (SFA) dengan model spasial (SAR, SEM, dan SDM) untuk menganalisis efisiensi teknis pertanian di Provinsi Hebei, Tiongkok [14]. Hasil penelitian mereka menyatakan adanya dependensi spasial signifikan dalam efisiensi pertanian, serta pengaruh penting dari faktor-faktor seperti teknologi, irigasi, dan struktur ekonomi terhadap produktivitas wilayah.
Dari hasil penelitian terdahulu menunjukkan analisis spasial mampu memberikan gambaran secara utuh terhadap dinamika kondisi ekonomi regional. Namun, kajian yang secara spesifik memodelkan faktor-faktor yang memengaruhi PDRB sektor pertanian dan peternakan antar provinsi di Indonesia dengan memperhatikan aspek keterkaitan spasial masih terbatas.
Bab 3. Metodologi
Data
Data yang digunakan pada penelitian ini merupakan data sekunder dan kuantitatif dengan dimensi cross-sectional pada tahun 2023. Sumber data berasal dari publikasi portal satu data pertanian. data terdiri dari 2 komponen utama :
Data Statistik
Data statistik yang digunakan dalam penelitian ini merupakan data sekunder yang berisi nilai-nilai variabel penelitian untuk setiap provinsi di Indonesia.Data Spasial
Peta digital dalam format file JSON yang memuat batas-batas administratif 38 Provinsi di Indonesia. Data ini diperlukan untuk mendefinisikan hubungan ketetanggaan dan melakukan visualisasi spasial.
Penelitian ini menggunakan satu variabel dependen (Y) dan tiga variabel independen (X) yang disajikan pada Tabel berikut.
| Kode | Variabel | Keterangan |
|---|---|---|
| X₁ | Komoditas Unggul Tanaman Pangan | Data Produksi Tanaman Pangan per provinsi tahun 2023 (Ton) |
| X₂ | Komoditas Unggul Peternakan | Data Komoditas Peternakan per provinsi tahun 2023 (ekor) |
| X₃ | Penduduk Berumur ≥15 Tahun yang Bekerja di Sektor Pertanian | Data ketenagakerjaan sektor pertanian per provinsi tahun 2023 (Jiwa) |
| Y | PDRB Pertanian, Perburuan, dan Peternakan | Data nilai PDRB sektor pertanian, perburuan, dan peternakan per provinsi tahun 2023 (Rupiah) |
Statistika Deskriptif
Statistika deskriptif merupakan cabang statistika yang bertujuan untuk menyajikan dan merangkum data sehingga informasi yang terkandung di dalamnya dapat dipahami dengan lebih mudah. Analisis ini mencakup perhitungan berbagai ukuran ringkasan data, seperti ukuran pemusatan dan ukuran penyebaran, antara lain mean, median, modus, kuartil atau kuantil, nilai maksimum, serta nilai minimum. Melalui statistika deskriptif, karakteristik umum dari setiap variabel penelitian dapat digambarkan secara ringkas dan sistematis.
Peta Deskriptif
Dalam analisis spasial, peta deskriptif digunakan untuk menampilkan variasi nilai suatu variabel antarwilayah. Sebagai contoh, peta ini dapat menunjukkan perbedaan PDRB di setiap provinsi
Peta deskriptif umumnya disajikan dalam bentuk peta tematik (choropleth map), di mana masing-masing wilayah diberi warna atau gradasi warna tertentu yang merepresentasikan besarnya nilai variabel yang diamati. Dengan demikian, pola keruangan dapat diamati secara visual.
Uji Autokorelasi Spasial
Uji autokorelasi spasial bertujuan untuk mengidentifikasi ada atau tidaknya keterkaitan antarwilayah yang dipengaruhi oleh kedekatan lokasi geografis. ## Autokorelasi Spasial dan Matriks Pembobot Spasial
Autokorelasi spasial atau spatial autocorrelation adalah kondisi di mana nilai suatu variabel pada satu lokasi berkorelasi dengan nilai variabel yang sama di lokasi-lokasi berdekatan. Fenomena ini sering muncul dalam data sosial-ekonomi karena efek limpahan (spillover), interaksi antarwilayah, persaingan, atau faktor tak teramati yang memiliki struktur spasial[15].
Untuk merepresentasikan struktur keterkaitan antar unit observasi, digunakan Matriks Pembobot Spasial (W), yaitu matriks persegi yang mendefinisikan hubungan spasial antar unit Matriks \(W\) secara umum dapat dituliskan sebagai: \[ W = \begin{pmatrix} 0 & w_{12} & w_{13} & \cdots & w_{1N} \\ w_{21} & 0 & w_{23} & \cdots & w_{2N} \\ w_{31} & w_{32} & 0 & \cdots & w_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ w_{N1} & w_{N2} & w_{N3} & \cdots & 0 \end{pmatrix} \] . Elemen $ w_{ij} $ mencerminkan derajat hubungan antara lokasi $ i $ dan $ j $, sedangkan diagonal biasanya nol \(( w_{ii}=0 )\). Matriks ini dapat dibangun berdasarkan kriteria kontiguitas wilayah, jarak geografis, maupun k-nearest neighbors, dan menjadi dasar perhitungan indeks autokorelasi spasial seperti Moran’s I atau Geary’s C[16].
Secara umum, autokorelasi spasial dapat diklasifikasikan menjadi:
- Positif, apabila wilayah-wilayah yang berdekatan memiliki nilai yang relatif serupa (tinggi–tinggi atau rendah–rendah).
- Negatif, apabila wilayah yang berdekatan justru menunjukkan nilai yang kontras (tinggi–rendah).
- Netral, apabila tidak ditemukan pola spasial yang jelas antarwilayah.
Moran’s I
Moran’s I merupakan instrumen statistik global yang berfungsi untuk mengidentifikasi kecenderungan pola spasial dari suatu variabel di seluruh wilayah amatan. Koefisien ini memiliki rentang nilai antara \(-1\) hingga \(+1\), di mana nilai tersebut mencerminkan derajat ketergantungan antar lokasi[17]. Secara spesifik, interpretasi indeks ini dibagi menjadi tiga kategori utama:
- Nilai Moran’s I > 1 : Mengindikasikan adanya pengelompokan (clustering), di mana lokasi-lokasi yang berdekatan cenderung memiliki nilai atau karakteristik yang serupa.
- Nilai Moran’s I < 1 :Menunjukkan adanya pola menyebar (dispersed), di mana lokasi-lokasi yang berdekatan cenderung memiliki nilai yang berbeda secara signifikan.
- Nilai Moran’s I \(\approx\) 0 : Menandakan tidak adanya keterkaitan spasial, yang berarti nilai suatu variabel terdistribusi secara independen dan tidak dipengaruhi oleh lokasi tetangganya.
Geary’s C
Geary’s C merupakan ukuran autokorelasi spasial global yang lebih sensitif terhadap perbedaan nilai antarwilayah yang bertetangga[18]. Indeks ini memiliki rentang nilai antara 0 hingga 2 dengan interpretasi sebagai berikut:
- Autokorelasi Spasial Positif (\(0 \le C < 1\)): Menggambarkan adanya kemiripan nilai antar wilayah yang bertetangga (homogenitas spasial). Semakin mendekati nol, semakin kuat pola pengelompokan (clustering) yang terbentuk.
- Autokorelasi Spasial Negatif (\(1 < C \le 2\)): Menunjukkan adanya perbedaan nilai yang signifikan antar wilayah yang berdekatan, yang mencerminkan pola persebaran atau ketimpangan.
- Independensi Spasial (\(C \approx 1\)): Menandakan tidak adanya keterkaitan spasial, di mana distribusi nilai pengamatan bersifat acak (random).
Getis-Ord Gi*
Getis-Ord \(G_i^*\) merupakan metode statistik spasial lokal yang dirancang untuk mendeteksi keberadaan klaster atau pengelompokan nilai ekstrem secara spasial, baik berupa titik panas (hot spots) untuk nilai tinggi maupun titik dingin (cold spots) untuk nilai rendah. Berbeda dengan indikator global seperti Moran’s I yang menyimpulkan pola spasial secara general di seluruh wilayah amatan, statistik \(G_i^*\) memberikan wawasan yang lebih mendalam mengenai intensitas nilai suatu variabel pada area tertentu beserta hubungannya dengan wilayah sekitar (neighboring areas)[19].
Dalam praktiknya, teknik ini sering diimplementasikan dalam berbagai disiplin ilmu, seperti ekonomi regional, epidemiologi, geografi kesehatan, hingga analisis infrastruktur, guna memetakan zona prioritas untuk intervensi kebijakan. Secara teknis, Getis-Ord \(G_i^*\) bekerja dengan mengidentifikasi lokasi yang memiliki signifikansi statistik secara spasial—di mana nilai-nilainya jauh lebih tinggi atau rendah dibandingkan ekspektasi distribusi acak—sehingga sangat andal dalam membedakan konsentrasi pembangunan atau mendeteksi ketimpangan antarwilayah.
LISA
LISA (Local Indicators of Spatial Association) merupakan instrumen statistik spasial yang digunakan untuk memetakan signifikansi hubungan antarwilayah pada tingkat lokal. Berbeda dengan Moran’s I atau Geary’s C yang memberikan kesimpulan mengenai pola spasial secara menyeluruh (global), LISA mampu mengidentifikasi heterogenitas spasial dengan mendeteksi konsentrasi nilai pada unit wilayah individual[20]. Analisis ini mengklasifikasikan wilayah ke dalam empat tipologi kuadran berdasarkan keterkaitannya dengan wilayah tetangga:
- High–High (HH): wilayah bernilai tinggi dikelilingi
wilayah bernilai tinggi
- Low–Low (LL): wilayah bernilai rendah dikelilingi
wilayah bernilai rendah
- High–Low (HL): wilayah bernilai tinggi dikelilingi
wilayah bernilai rendah
- Low–High (LH): wilayah bernilai rendah dikelilingi wilayah bernilai tinggi
Model Spasial
OLS (Ordinary Least Square)
Model Ordinary Least Squares (OLS) merupakan kerangka estimasi fundamental dalam analisis regresi yang beroperasi di bawah asumsi independensi antar observasi. Dalam konteks analisis data bereferensi geografis, OLS berfungsi sebagai model acuan (baseline) untuk mendeteksi potensi penyimpangan asumsi klasik sebelum beralih ke estimasi spasial yang lebih kompleks [21]. Struktur model ini secara matematis dinyatakan sebagai berikut:
\[\mathbf{Y} = \mathbf{X}\beta + \varepsilon\]
Di mana: - \(\mathbf{Y}\) : merupakan vektor untuk variabel dependen dengan dimensi \(n \times 1\). - \(\mathbf{X}\) : merupakan matriks variabel independen berukuran \(n \times k\). - \(\beta\) : merupakan vektor koefisien regresi yang diestimasi. - \(\varepsilon\) : merupakan vektor sisaan atau galat acak (random error).
Pelaksanaan uji diagnostik terhadap residual model OLS, terutama pengujian aspek autokorelasi spasial, merupakan tahapan krusial. Hasil pengujian ini menjadi dasar pertimbangan ilmiah dalam menentukan apakah model perlu dikembangkan menggunakan pendekatan ekonometrika spasial guna mengakomodasi ketergantungan antarwilayah.
SAR (Spatial Autoregressive Model)
SAR merupakan bentuk pengembangan dari regresi linier yang digunakan ketika variabel dependen menunjukkan adanya ketergantungan spasial. Dalam model ini, nilai variabel dependen pada suatu lokasi diasumsikan memiliki keterkaitan langsung dengan nilai variabel dependen di lokasi-lokasi tetangganya. Hubungan ini menciptakan efek lag spasial (\(\rho Wy\)), di mana fenomena di suatu wilayah tidak hanya dipengaruhi oleh faktor internal (variabel independen), tetapi juga dipengaruhi oleh interaksi atau difusi dari wilayah sekitarnya[22].
Persamaan model SAR adalah:
\[
\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} +
\boldsymbol{\varepsilon}
\] Keterangan:
- \(\mathbf{y}\): vektor variabel
dependen
- \(\mathbf{X}\): matriks variabel
independen
- \(\boldsymbol{\beta}\): vektor
koefisien regresi
- \(\rho\): parameter autokorelasi
spasial (spatial lag coefficient)
- \(\mathbf{W}\): matriks bobot spasial
(spatial weight matrix)
- \(\boldsymbol{\varepsilon}\): vektor
error acak
Model SAR mengasumsikan bahwa nilai variabel dependen di suatu lokasi dipengaruhi oleh nilai-nilai variabel dependen di lokasi-lokasi tetangganya.
SEM (Spatial Error Model)
SEM merupakan pendekatan ekonometrika spasial yang digunakan apabila ketergantungan antarwilayah tidak ditemukan pada variabel dependen secara langsung, melainkan terakumulasi dalam komponen sisaan atau galat (error term). Fenomena ini sering kali mencerminkan adanya variabel-variabel relevan yang tidak teramati (omitted variables) atau faktor eksternal yang polanya melintasi batas-batas wilayah amatan[23]. Dalam model ini, korelasi spasial ditangani dengan menyertakan koefisien autoregresif pada struktur galatnya.
Persamaan model SEM adalah:
\[
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},
\quad
\boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} +
\mathbf{u}
\] Keterangan:
- \(\lambda\): parameter autokorelasi
spasial pada error
- \(\mathbf{u}\): komponen error acak
yang bebas spasial
SDM (Spatial Durbin Model)
Spatial Durbin Model (SDM) merupakan perluasan dari model SAR yang menawarkan kerangka analisis lebih komprehensif dengan mengintegrasikan ketergantungan spasial pada variabel dependen maupun variabel independen. Dalam model ini, lag spasial tidak hanya diterapkan pada variabel respon (\(Wy\)), tetapi juga pada variabel penjelas (\(WX\)). Karakteristik utama yang membedakan SDM adalah kemampuannya dalam menangkap efek langsung (direct effects)—yaitu pengaruh variabel independen di suatu wilayah terhadap wilayah itu sendiri—serta efek limpahan (spillover effects)—yakni pengaruh variabel independen dari wilayah tetangga terhadap wilayah amatan[24].
Persamaan model SDM adalah:
\[
\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} +
\mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon}
\] Keterangan:
- \(\mathbf{W}\mathbf{X}\): lag spasial
dari variabel independen
- \(\boldsymbol{\theta}\): koefisien
efek tidak langsung (spillover effects)
SDEM (Spatial Durbin Error Model)
Spatial Durbin Error Model (SDEM) merupakan varian model ekonometrika spasial yang dikembangkan dari struktur Spatial Error Model (SEM) dengan mengintegrasikan efek lag spasial pada variabel independen (\(WX\)). Berbeda dengan model SDM, SDEM tidak menyertakan interaksi spasial pada variabel dependen (\(Wy\)). Model ini sangat efektif untuk menganalisis fenomena di mana pengaruh eksternal dari wilayah tetangga berasal dari variabel-variabel penjelas tertentu (efek spillover lokal) serta adanya keterkaitan pada komponen galat yang tidak terobservasi[25].
Persamaan model SDEM adalah:
\[
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} +
\mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon},
\quad
\boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} +
\mathbf{u}
\] Keterangan:
- Model ini menggabungkan lag spasial pada variabel independen (\(\mathbf{W}\mathbf{X}\)) dengan autokorelasi
error (\(\lambda
\mathbf{W}\boldsymbol{\varepsilon}\)).
- Tidak ada lag pada variabel dependen (\(\mathbf{W}\mathbf{y}\)).
SAC (Spatial Autocorrelation Model)
Model Spatial Autocorrelation (SAC) merupakan model ekonometrika spasial yang bersifat komprehensif karena mengintegrasikan efek lag spasial pada variabel dependen sekaligus autokorelasi pada komponen galat secara simultan. Instrumen ini diimplementasikan ketika terdapat indikasi empiris bahwa fenomena ketergantungan antarwilayah terjadi melalui dua mekanisme tersebut secara bersamaan. Dalam berbagai riset tingkat lanjut, SAC sering kali diposisikan sebagai model pembanding utama guna membedah kompleksitas proses serta interaksi spasial yang terkandung dalam data kewilayahan[26].
Persamaan model SAC adalah : \[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} = \lambda \mathbf{W}\boldsymbol{\varepsilon} + \mathbf{u} \]
Keterangan:
- Menggabungkan dua bentuk dependensi spasial:
- Lag spasial pada \(y\) seperti pada SAR
- Autokorelasi spasial pada error seperti pada SEM
Interpolasi
Interpolasi spasial merupakan prosedur statistika yang digunakan untuk mengestimasi nilai variabel pada lokasi-lokasi yang tidak tercakup dalam sampel berdasarkan data pada lokasi yang telah diketahui, dengan berlandaskan pada prinsip dependensi atau kedekatan geografis. Teknik ini secara luas diaplikasikan dalam analisis geografi serta perencanaan infrastruktur untuk mentransformasi data yang bersifat diskret atau agregat wilayah menjadi representasi spasial yang bersifat kontinu. Hal ini memungkinkan peneliti untuk memvisualisasikan pola sebaran suatu fenomena secara lebih menyeluruh dan mendalam, meskipun basis data awal yang tersedia terbatas pada titik-titik tertentu saja[27].
IDW (Inverse Distance Wighting)
Inverse Distance Weighting (IDW) adalah metode interpolasi deterministik yang mengestimasi nilai pada lokasi yang tidak diamati dengan asumsi bahwa pengaruh titik pengamatan akan menurun seiring bertambahnya jarak. Nilai di lokasi estimasi \(s_0\) dihitung sebagai rata-rata berbobot dari nilai variabel di titik-titik sekitar, di mana bobotnya ditentukan oleh jarak ke titik estimasi[28].
Secara matematis, estimasi IDW dituliskan sebagai:
\[ \hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i Z(s_i)}{\sum_{i=1}^{n} w_i}, \quad w_i = \frac{1}{d(s_0, s_i)^p} \]
Keterangan: - \(\hat{Z}(s_0)\) :
nilai estimasi pada lokasi \(s_0\)
- \(Z(s_i)\) : nilai variabel pada
lokasi pengamatan ke-\(i\)
- \(d(s_0, s_i)\) : jarak antara lokasi
estimasi \(s_0\) dan titik pengamatan
\(s_i\)
- \(p\) : parameter pangkat (power
parameter) yang mengontrol penurunan bobot terhadap jarak
- \(n\) : jumlah titik tetangga yang
digunakan dalam perhitungan
Ordinary Kriging
Ordinary Kriging adalah metode interpolasi geostatistik yang digunakan untuk memperkirakan nilai variabel pada lokasi yang tidak diamati dengan mempertimbangkan ketergantungan spasial antar titik pengamatan melalui variogram[29]. Metode ini berasumsi bahwa nilai rata-rata variabel bersifat konstan tetapi tidak diketahui di seluruh wilayah studi.
Estimasi nilai pada lokasi yang tidak diamati melalui Ordinary Kriging dirumuskan sebagai berikut:
\[ \hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i) \]
dengan kondisi pembobot:
\[ \sum_{i=1}^{n} \lambda_i = 1 \]
dan persamaan sistem untuk menentukan bobot \(\lambda_i\):
\[ \sum_{j=1}^{n} \lambda_j \gamma(s_i - s_j) + \mu = \gamma(s_i - s_0), \quad i = 1,2,\dots,n \]
Keterangan: - \(\lambda_i\) : bobot
Kriging pada lokasi ke-\(i\)
- \(\gamma(h)\) : fungsi variogram pada
jarak \(h\)
- \(\mu\) : pengali Lagrange
- \(s_0\) : lokasi estimasi
- \(Z(s_i)\) : nilai variabel pada
lokasi pengamatan \(i\)
Bab 4. Hasil dan Pembahasan
Load Data dan Library serta pengolahan data
# Library
library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## 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(sf)## Warning: package 'sf' was built under R version 4.4.3## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(sp)## Warning: package 'sp' was built under R version 4.4.3library(spdep)## Warning: package 'spdep' was built under R version 4.4.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.4.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.4.3## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(gstat)## Warning: package 'gstat' was built under R version 4.4.3library(tmap)## Warning: package 'tmap' was built under R version 4.4.3library(leaflet)## Warning: package 'leaflet' was built under R version 4.4.3library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode#Data
data <- read.csv("C:/my data/Kuliah/Sem 5/Spatial Statistics/uas/data uas spasial.csv",sep = ",",dec = ".")
data## Provinsi PDRB_Y Total_Komoditas_Unggulan_Pangan_X1
## 1 Aceh 54833655 5.708637e+04
## 2 Sumatera Utara 220428802 1.217843e+06
## 3 Sumatera Barat 48591004 2.654659e+05
## 4 Riau 196956661 6.207234e+04
## 5 Jambi 83422451 8.659505e+04
## 6 Sumatera Selatan 56520214 2.379734e+05
## 7 Bengkulu 20155726 5.280604e+04
## 8 Lampung 98698490 7.266867e+06
## 9 Kepulauan Bangka Belitung 12050335 5.059921e+04
## 10 Kepulauan Riau 3971384 1.829498e+04
## 11 DKI Jakarta 1341915 0.000000e+00
## 12 Jawa Barat 195134345 1.650598e+06
## 13 Jawa Tengah 198539321 2.990917e+06
## 14 DAERAH ISTIMEWA Yogyakarta 16887617 1.155494e+06
## 15 Jawa Timur 244174125 2.047216e+06
## 16 Banten 42073957 7.681086e+04
## 17 Bali 25582226 1.047284e+05
## 18 Nusa Tenggara Barat 29082888 1.109844e+05
## 19 Nusa Tenggara Timur 29558385 4.763711e+05
## 20 Kalimantan Barat 50267248 1.692355e+05
## 21 Kalimantan Tengah 38223491 5.380800e+04
## 22 Kalimantan Selatan 21020330 4.680864e+04
## 23 Kalimantan Timur 47126293 5.749416e+04
## 24 Kalimantan Utara 7106315 3.379932e+04
## 25 Sulawesi Utara 21231518 5.278240e+04
## 26 Sulawesi Tengah 41426258 5.420298e+04
## 27 Sulawesi Selatan 84955141 4.087608e+05
## 28 Sulawesi Tenggara 19447995 1.218174e+05
## 29 Gorontalo 13889772 5.951854e+03
## 30 Sulawesi Barat 19071683 2.926421e+04
## 31 Maluku 5614880 7.211234e+04
## 32 Maluku Utara 6719216 1.352929e+04
## 33 Papua Barat 1659179 1.141092e+04
## 34 Papua Barat Daya 968190 0.000000e+00
## 35 Papua 3831579 7.492688e+01
## 36 Papua Selatan 3190763 0.000000e+00
## 37 Papua Tengah 3281047 0.000000e+00
## 38 Papua Pegunungan 3840581 0.000000e+00
## Total_Komoditas_Unggulan_Peternakan_X2 Tenaga_Kerja_Petani_X3
## 1 36638331 946649
## 2 202408013 2559037
## 3 78334563 964493
## 4 90842075 1146960
## 5 52763949 814401
## 6 126883418 1967142
## 7 11083474 484033
## 8 104523223 1987969
## 9 22775942 188091
## 10 18784792 85171
## 11 3929 23705
## 12 770364401 3603172
## 13 679126990 488360
## 14 61650069 420457
## 15 559399243 7149540
## 16 208245113 662474
## 17 91536757 495732
## 18 42599944 970339
## 19 17176634 1421199
## 20 63422822 1318010
## 21 32907080 485968
## 22 102439156 612834
## 23 61727241 354868
## 24 4909397 123199
## 25 14012582 285522
## 26 11007986 627840
## 27 114962752 1609771
## 28 9390764 416507
## 29 6977372 188127
## 30 4216932 358508
## 31 799416 269604
## 32 150749 160825
## 33 381916 199599
## 34 1029279 0
## 35 4719348 1666431
## 36 721737 0
## 37 1869139 0
## 38 56779 0# Transformasi variabel
Provinsi <- toupper(data$Provinsi)
Y <- log(as.numeric(data$PDRB_Y) + 1)
X1 <- log(as.numeric(data$Total_Komoditas_Unggulan_Pangan_X1) + 1)
X2 <- log(as.numeric(data$Total_Komoditas_Unggulan_Peternakan_X2) + 1)
X3 <- log(as.numeric(data$Tenaga_Kerja_Petani_X3) + 1)
# Data frame analisis
data_tr <- data.frame(Provinsi, Y, X1, X2, X3)
# Baca shapefile provinsi
indo_lvl1 <- st_read("C:/my data/Kuliah/Sem 5/Spatial Statistics/uas/Batas_Provinsi_Indonesia (4)/Batas_Provinsi_Indonesia.shp")## Reading layer `Batas_Provinsi_Indonesia' from data source
## `C:\my data\Kuliah\Sem 5\Spatial Statistics\uas\Batas_Provinsi_Indonesia (4)\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(data_tr$P, indo_lvl1$WADMPR)## character(0)Analisis Deskriptif
Statistik Deskriptif Data Aktual
summary(data[,2:5])## PDRB_Y Total_Komoditas_Unggulan_Pangan_X1
## Min. : 968190 Min. : 0
## 1st Qu.: 6815991 1st Qu.: 21037
## Median : 23406872 Median : 57290
## Mean : 51865131 Mean : 501573
## 3rd Qu.: 53692053 3rd Qu.: 220789
## Max. :244174125 Max. :7266867
## Total_Komoditas_Unggulan_Peternakan_X2 Tenaga_Kerja_Petani_X3
## Min. : 3929 Min. : 0
## 1st Qu.: 4766860 1st Qu.: 190995
## Median : 27841511 Median : 487164
## Mean : 95022192 Mean : 922540
## 3rd Qu.: 91363086 3rd Qu.:1102805
## Max. :770364401 Max. :7149540Berdasarkan hasil variasi yang cukup besar pada tiap variabel, pada PDRB (Y) nilai minimum berada pada Rp968,190 tidak sampai menyentuh angka satu juta rupiah sementara nilai maksimum PDRB sebesar Rp244,174,125. Variabel komoditas unggulan tanaman pangan (X1) juga menunjukkan variasi yang sangat besar. Oleh karena terdapat nilai-nilai yang cukup ekstrem dilakukan transformasi data dengan menggunakan log.
Statistik Deskriptif Data Transformasi
summary(data_tr[, 2:5])## Y X1 X2 X3
## Min. :13.78 Min. : 0.000 Min. : 8.276 Min. : 0.00
## 1st Qu.:15.73 1st Qu.: 9.932 1st Qu.:15.377 1st Qu.:12.16
## Median :16.96 Median :10.956 Median :17.125 Median :13.10
## Mean :16.89 Mean : 9.990 Mean :16.541 Mean :11.87
## 3rd Qu.:17.80 3rd Qu.:12.295 3rd Qu.:18.330 3rd Qu.:13.91
## Max. :19.31 Max. :15.799 Max. :20.462 Max. :15.78Setelah dilakukan Transformasi terlihat variasi data tidak terlalu besar seperti pada data aktual sebelumnya.
Peta Spasial PDRB
# Merge Data
indo_merged <- indo_lvl1 %>%
left_join(data_tr, by = c("WADMPR" = "Provinsi")) %>%
st_zm(drop = TRUE, what = "ZM")
indo_merged$Y <- as.numeric(indo_merged$Y)# Peta Spasial
q_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 = q_breaks,
palette = "YlOrRd",
colorNA = "gray80",
border.col = "black",
border.alpha = 1,
title = "Nilai PDRB"
) +
tm_layout(
title = "Peta Sebaran PDRB Pertanian, Perburuan, dan Pendidikan 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"
Dari hasil visualisasi terlihat konsentrasi PDRB Pertanian, Perburuan,
dan Peternakan condong ada di Indonesia bagian barat yaitu Sumatera dan
Jawa, Sementara di Indonesia bagian tengah terlihat hasil yang juga
cukup menengah dengan ditandai oleh warna oranye. Sementara untuk di
daerah timur terlihat jauh lebih sedikit, bisa jadi ini terjadi karena
kondisi geografis yang berbeda antara Indonesia barat dan timur.
Autokorelasi Spasial
Membentuk Matriks Spasial
indo_fixed <- indo_merged %>% filter(!is.na(Y))
indo_fixed <- st_make_valid(indo_fixed)
coords <- st_coordinates(st_centroid(indo_fixed))## Warning: st_centroid assumes attributes are constant over geometries# Pemilihan k optimal
k_values <- 1:10
moran_I <- sapply(k_values, function(k) {
nb_k <- knn2nb(knearneigh(coords, k = k))
lw_k <- nb2listw(nb_k, style = "W")
moran.test(indo_fixed$Y, lw_k)$estimate[1]
})## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 12
## sub-graphs## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 2 sub-graphsk_optimal <- k_values[which.max(moran_I)]
k_optimal## [1] 3# Bobot spasial kNN
nb_knn <- knn2nb(knearneigh(coords, k = k_optimal))
lw_knn <- nb2listw(nb_knn, style = "W")Morans’ I
moran.test(indo_fixed$Y, lw_knn)##
## Moran I test under randomisation
##
## data: indo_fixed$Y
## weights: lw_knn
##
## Moran I statistic standard deviate = 4.1602, p-value = 1.59e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.46690667 -0.02702703 0.01409621Hasil uji autokorelasi global dengan menggunakan statistik Moran’s I menunjukkan adanya dependensi secara spasial yang positif pada data, dengan nilai statistik Morans’s I sebesar 0.4669. Ini mengindikasikan adanya pengelompokkan spasial yang kuat antar wilayah dengan karakteristik yang mirip. Hal ini juga terbukti signifikan secara statistik dengan p-value sebesar 1.59e-05, yang menjadi alasan untuk tolak h0 dan menyimpulkan bahwa sebaran data antarwilayah tidak bersifat acak, melainkan dipengaruhi oleh lokasi geografis tetangganya. #Geary’s C
geary.test(indo_fixed$Y, lw_knn)##
## Geary C test under randomisation
##
## data: indo_fixed$Y
## weights: lw_knn
##
## Geary C statistic standard deviate = 3.6372, p-value = 0.0001378
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.55662473 1.00000000 0.01485971Agar bukti autokorelasi spasial lebih konsisten dilakukan juga uji statistik Geary’s C. Hasil dari Uji Geary’s C memberikan nilai sebesar 0.5566, ini mengonfirmasi benar adanya autokorelasi spasial positif di mana wilayah yang bertetangga cenderung memiliki nilai atribut yang mirip. Signifkansi dari uji ini juga menunjukkan hasil yang sesuai yaitu dengan p-value sebesar 0.0001, yang menjadi alasan untuk tolak h0 dan memperkuat bahwa terdapat pengelompokkan pada data bersifat nyata secara statistik dan bukan akibat proses acak. ### Getis-Ord Gi*
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*\nNilai PDRB Pertanian, Perburuan, dan Peternakan",
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 DI Yogyakarta teridentifikasi sebagai
hotspot utama dengan tingkat kepercayaan 99% (merah) dan wilayah Aceh
sebagai hotspot sekunder pada taraf 95% (oranye), untuk wilayah timur
terutama Papua terdapat pengelompokkan daerah cold spot. Hal ini
mengindikasikan terdapat konsentrasi wilayah dengan PDRB Pertanian,
Perburuan, dan Peternakanyang secara signifikan dikelilingi oleh wilayah
berkualitas serupa, memperlihatkan pengaruh kondisi geografis yang
berbeda. ### 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 wilayah
mean_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", "darkgreen", "lightblue", "white", "grey"),
title = "LISA Cluster",
style = "cat") +
tm_borders(lwd = 1, col = "black") +
tm_layout(main.title = "Peta LISA PDRB \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 menunjukkan ada klaster High-High yang terdapat di
Pulau Sumatera yaitu Aceh, yang berarti daerah di sekitar aceh
menunjukkan nilai PDRB yang tinggi juga, kemudian ada klaster Low-High
di DI Yogtakarta, ini menunjukkan bahwa DI Yogyakarta memiliki nilai
PDRB yang rendah akan tetapi dikelilingi oleh daerah yang memiliki PDRB
tinggi. Lalu juga ada Papua dan sekitarnya yang merupakan klaster
Low-Low, ini merupakan kumpulan daerah yang memiliki PDRB yang
rendah.
Model Spasial
Model OLS
ols_model <- lm(Y ~ X1 + X2 + X3, data = indo_fixed)
summary(ols_model)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = indo_fixed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.73503 -0.42900 -0.00399 0.44646 1.44785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.01479 1.01814 9.836 1.78e-11 ***
## X1 0.08061 0.07276 1.108 0.276
## X2 0.34516 0.07594 4.545 6.62e-05 ***
## X3 0.03020 0.05859 0.515 0.610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.747 on 34 degrees of freedom
## Multiple R-squared: 0.7691, Adjusted R-squared: 0.7487
## F-statistic: 37.75 on 3 and 34 DF, p-value: 6.344e-11Model OLS berhasil dibuat, akan tetapi perlu dilakukan uji asumsi apakah OLS layak dipakai atau tidak. #### Asumsi Residual ##### Normalitas
# Shapiro-Wilk test (Normalitas)
shapiro.test(residuals(ols_model))##
## Shapiro-Wilk normality test
##
## data: residuals(ols_model)
## W = 0.98412, p-value = 0.8552# QQ Plot
qqnorm(residuals(ols_model))
qqline(residuals(ols_model), col = "red")
Heteroskedastisitas
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(ols_model)##
## studentized Breusch-Pagan test
##
## data: ols_model
## BP = 2.0671, df = 3, p-value = 0.5586Multikolinearitas
# Uji Multikolinearitas
vif(ols_model)## X1 X2 X3
## 6.808364 2.696101 4.141048Berdasarkan hasil uji di atas, asumsi Normalitas dan Heteroskedastisitas memenuhi asumsi dengan nilai p-value > alpha dan untuk Multikolinearitas juga terpenuhi karena semua nilai VIF < 10. Namun karena pada uji autokorelasi spasial terbukti ada pengaruh spasial, pemodelan ekonometrika spasial tetap dilakukan. ### Model SAR
sar_model <- lagsarlm(Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn)
summary(sar_model)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.856821 -0.325128 0.039842 0.451958 1.357313
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.100029 1.843370 4.3941 1.112e-05
## X1 0.082901 0.067428 1.2295 0.2189
## X2 0.326792 0.070543 4.6325 3.613e-06
## X3 0.012107 0.054816 0.2209 0.8252
##
## Rho: 0.14133, LR test value: 1.5979, p-value: 0.20619
## Asymptotic standard error: 0.11043
## z-value: 1.2798, p-value: 0.20061
## Wald statistic: 1.6379, p-value: 0.20061
##
## Log likelihood: -39.9234 for lag model
## ML residual variance (sigma squared): 0.47621, (sigma: 0.69008)
## Number of observations: 38
## Number of parameters estimated: 6
## AIC: 91.847, (AIC for lm: 91.445)
## LM test for residual autocorrelation
## test value: 0.18502, p-value: 0.66709Hasil estimasi model Spatial Autoregressive (SAR) menunjukkan bahwa koefisien autoregresif spasial sebesar 0.141, ini tidak signifikan secara statistik, yang dibuktikan oleh nilai p-value uji Likelihood Ratio sebesar 0.206, sehingga tidak terdapat bukti empiris yang kuat mengenai adanya efek limpahan spasial pada variabel dependen. Artinya, nilai variabel dependen di suatu provinsi tidak dipengaruhi secara signifikan oleh nilai rata-rata variabel dependen di provinsi tetangganya setelah variabel independen dimasukkan ke dalam model. Dari sisi pengaruh variabel independen, hanya X2 yang berpengaruh positif dan signifikan terhadap variabel dependen, sedangkan X1 dan X3 tidak signifikan. Selain itu, nilai AIC model SAR (91.847) yang sedikit lebih tinggi dibandingkan AIC model OLS (91.445) serta hasil LM test residual yang tidak signifikan (p-value 0.667) menunjukkan bahwa penambahan komponen spasial tidak meningkatkan kinerja model, sehingga model OLS dinilai lebih memadai untuk data yang dianalisis.
Model SEM
sem_model <- errorsarlm(Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn)
summary(sem_model)##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7626879 -0.4164904 -0.0061262 0.4409365 1.4408206
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.221082 0.969146 10.5465 < 2.2e-16
## X1 0.091358 0.068438 1.3349 0.1819
## X2 0.331462 0.071653 4.6259 3.729e-06
## X3 0.021800 0.054475 0.4002 0.6890
##
## Lambda: 0.081546, LR test value: 0.08304, p-value: 0.77322
## Asymptotic standard error: 0.20454
## z-value: 0.39868, p-value: 0.69013
## Wald statistic: 0.15894, p-value: 0.69013
##
## Log likelihood: -40.68085 for error model
## ML residual variance (sigma squared): 0.49733, (sigma: 0.70522)
## Number of observations: 38
## Number of parameters estimated: 6
## AIC: 93.362, (AIC for lm: 91.445)Hasil estimasi model Spatial Error Model (SEM) menunjukkan bahwa parameter dependensi spasial pada error sebesar 0.082, ini tidak signifikan secara statistik, yang ditunjukkan oleh nilai p-value uji Likelihood Ratio sebesar 0.773, sehingga tidak terdapat bukti adanya ketergantungan spasial yang tersisa pada komponen galat setelah variabel independen dimasukkan ke dalam model. Hal ini mengindikasikan bahwa variasi variabel dependen tidak dipengaruhi oleh faktor-faktor tak teramati yang memiliki pola spasial. Secara parsial, hanya X2 yang berpengaruh positif dan signifikan terhadap variabel dependen, sedangkan X1 dan X3 tidak signifikan secara statistik. Dari sisi kinerja model, nilai AIC model SEM (93.362) lebih tinggi dibandingkan AIC model OLS (91.445), yang menunjukkan bahwa penambahan komponen spasial pada error tidak meningkatkan kualitas model, sehingga model regresi OLS dinilai lebih sesuai untuk data yang dianalisis. ### Model SDM
sdm_model <- lagsarlm(
Y ~ X1 + X2 + X3,
data = indo_fixed,
listw = lw_knn,
type = "mixed"
)
summary(sdm_model)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.954402 -0.324624 0.017642 0.349118 1.338053
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.218916 2.285911 3.1580 0.001589
## X1 0.038607 0.066443 0.5811 0.561206
## X2 0.331590 0.066629 4.9767 6.468e-07
## X3 0.064094 0.057023 1.1240 0.261017
## lag.X1 -0.319084 0.129285 -2.4681 0.013584
## lag.X2 0.327411 0.139761 2.3427 0.019147
## lag.X3 0.265756 0.118759 2.2378 0.025236
##
## Rho: -0.13992, LR test value: 0.3434, p-value: 0.55787
## Asymptotic standard error: 0.21585
## z-value: -0.64822, p-value: 0.51684
## Wald statistic: 0.42019, p-value: 0.51684
##
## Log likelihood: -36.56514 for mixed model
## ML residual variance (sigma squared): 0.39926, (sigma: 0.63187)
## Number of observations: 38
## Number of parameters estimated: 9
## AIC: 91.13, (AIC for lm: 89.474)
## LM test for residual autocorrelation
## test value: 0.14528, p-value: 0.70309Hasil estimasi Spatial Durbin Model (SDM/mixed) menunjukkan bahwa koefisien autoregresif spasial sebesesar 0.-140, ini tidak signifikan secara statistik, yang ditunjukkan oleh p-value uji Likelihood Ratio sebesar 0.558, sehingga tidak terdapat bukti kuat adanya efek ketergantungan spasial langsung pada variabel dependen. Secara parsial, variabel X2 berpengaruh positif dan signifikan terhadap variabel dependen, sedangkan X1 dan X3 tidak signifikan. Namun demikian, koefisien spasial dari variabel independen (lag.X1, lag.X2, dan lag.X3) menunjukkan pengaruh yang signifikan, di mana lag.X1 berpengaruh negatif, sedangkan lag.X2 dan lag.X3 berpengaruh positif, yang mengindikasikan adanya spillover spasial dari variabel penjelas antar wilayah. Dari sisi kinerja model, nilai AIC model SDM sebesar 91.13 lebih tinggi dibandingkan AIC model OLS (89.47) dan uji LM residual yang tidak signifikan (p-value 0.703) menunjukkan bahwa meskipun terdapat indikasi spillover pada variabel independen, penambahan kompleksitas spasial belum memberikan peningkatan kualitas model yang lebih baik dibandingkan model OLS. ### SDEM
sdem_model <- errorsarlm(
Y ~ X1 + X2 + X3,
data = indo_fixed,
listw = lw_knn,
Durbin = TRUE
)
summary(sdem_model)##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn,
## Durbin = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.969017 -0.291516 0.022835 0.345245 1.324386
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.195083 1.500211 4.1295 3.636e-05
## X1 0.049922 0.064636 0.7724 0.439904
## X2 0.323020 0.065643 4.9209 8.615e-07
## X3 0.054493 0.055067 0.9896 0.322384
## lag.X1 -0.307599 0.125696 -2.4472 0.014399
## lag.X2 0.265431 0.102111 2.5994 0.009337
## lag.X3 0.241177 0.111253 2.1678 0.030172
##
## Lambda: -0.090852, LR test value: 0.14095, p-value: 0.70734
## Asymptotic standard error: 0.21887
## z-value: -0.41509, p-value: 0.67808
## Wald statistic: 0.1723, p-value: 0.67808
##
## Log likelihood: -36.66637 for error model
## ML residual variance (sigma squared): 0.40249, (sigma: 0.63442)
## Number of observations: 38
## Number of parameters estimated: 9
## AIC: 91.333, (AIC for lm: 89.474)Hasil estimasi Spatial Durbin Error Model (SDEM) menunjukkan bahwa parameter ketergantungan spasial pada komponen galat −0.091, ini tidak signifikan secara statistik, yang ditunjukkan oleh p-value uji Likelihood Ratio sebesar 0.707, sehingga tidak terdapat bukti adanya dependensi spasial yang tersisa pada error model. Secara parsial, hanya X2 yang berpengaruh positif dan signifikan terhadap variabel dependen, sedangkan X1 dan X3 tidak signifikan. Namun demikian, koefisien lag spasial dari variabel independen (lag.X1, lag.X2, dan lag.X3) signifikan secara statistik, dengan lag.X1 berpengaruh negatif serta lag.X2 dan lag.X3 berpengaruh positif, yang mengindikasikan adanya efek limpahan spasial dari variabel penjelas antar wilayah. Dari sisi kinerja model, nilai AIC model SDEM sebesar 91.33 masih lebih tinggi dibandingkan AIC model OLS (89.47), sehingga meskipun terdapat indikasi spillover pada variabel independen, model ini belum menunjukkan perbaikan kinerja dibandingkan model OLS. ### SAC
sac_model <- sacsarlm(
Y ~ X1 + X2 + X3,
data = indo_fixed,
listw = lw_knn,
zero.policy = TRUE
)
summary(sac_model)##
## Call:sacsarlm(formula = Y ~ X1 + X2 + X3, data = indo_fixed, listw = lw_knn,
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.845578 -0.346317 0.066614 0.471372 1.370439
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.284679 1.950298 3.7352 0.0001876
## X1 0.057624 0.067143 0.8582 0.3907665
## X2 0.344042 0.070578 4.8746 1.09e-06
## X3 0.023282 0.056651 0.4110 0.6810959
##
## Rho: 0.18072
## Asymptotic standard error: 0.12945
## z-value: 1.3961, p-value: 0.1627
## Lambda: -0.18512
## Asymptotic standard error: 0.26836
## z-value: -0.68981, p-value: 0.49031
##
## LR test value: 1.9041, p-value: 0.38594
##
## Log likelihood: -39.7703 for sac model
## ML residual variance (sigma squared): 0.46689, (sigma: 0.68329)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 93.541, (AIC for lm: 91.445)Hasil estimasi Spatial Autoregressive Combined (SAC) menunjukkan bahwa baik ketergantungan spasial pada variabel dependen 0.181 maupun pada komponen galat 0.185) tidak signifikan secara statistik, yang diperkuat oleh p-value uji LR sebesar 0.386, sehingga tidak terdapat bukti kuat adanya dependensi spasial simultan dalam model. Secara parsial, hanya X2 yang berpengaruh positif dan signifikan terhadap variabel dependen, sedangkan X1 dan X3 tidak signifikan. Dari sisi kecocokan model, nilai AIC sebesar 93.54 lebih besar dibandingkan AIC model OLS (91.45), yang mengindikasikan bahwa model SAC belum memberikan peningkatan kinerja dibandingkan model regresi linier tanpa mempertimbangkan efek spasial.
Pemilihan Model Terbaik
model_comparison <- data.frame(
Model = c("OLS", "SAR", "SEM", "SDEM", "SDM", "SAC"),
AIC = c(
AIC(ols_model),
AIC(sar_model),
AIC(sem_model),
AIC(sdem_model),
AIC(sdm_model),
AIC(sac_model)
),
LogLik = c(
as.numeric(logLik(ols_model)),
as.numeric(logLik(sar_model)),
as.numeric(logLik(sem_model)),
as.numeric(logLik(sdem_model)),
as.numeric(logLik(sdm_model)),
as.numeric(logLik(sac_model))
)
)
model_comparison[order(model_comparison$AIC), ]## Model AIC LogLik
## 5 SDM 91.13028 -36.56514
## 4 SDEM 91.33273 -36.66637
## 1 OLS 91.44475 -40.72237
## 2 SAR 91.84681 -39.92340
## 3 SEM 93.36171 -40.68085
## 6 SAC 93.54060 -39.77030Berdasarkan perbandingan kriteria AIC dan nilai log-likelihood, model Spatial Durbin Model (SDM) merupakan model terbaik karena memiliki nilai AIC paling kecil (91.13) dan log-likelihood tertinggi, menunjukkan keseimbangan terbaik antara kecocokan model dan kompleksitas.Hasil ini menegaskan bahwa memasukkan efek lag spasial pada variabel independen sebagaimana pada SDM memberikan peningkatan kinerja model yang paling signifikan dalam menjelaskan variasi antarwilayah.
# Model terbaik
best_model <- sdm_modelindo_fixed$Y_pred_log <- as.numeric(predict(best_model, type = "response"))## This method assumes the response is known - see manual page# BACK-TRANSFORM KE SKALA ASLI
# Karena Y = log(Y + 1)
indo_fixed$Y_pred <- exp(indo_fixed$Y_pred_log) - 1
# RESIDUAL (TETAP DI SKALA LOG)
indo_fixed$residual_log <- as.numeric(residuals(best_model))
# Residual distandarisasi (untuk diagnostik spasial)
indo_fixed$resid_z <- as.numeric(scale(indo_fixed$residual_log))Visualisasi Model Terbaik
library(tmap)
tm_shape(indo_fixed) +
tm_fill(
col = "Y_pred",
palette = "viridis",
style = "quantile",
n = 7,
title = "Hasil Prediksi Model Spasial"
) +
tm_borders(col = "gray30", lwd = 0.6) +
tm_layout(
main.title = "Peta Prediksi PDRB Pertanian, Perburuan, dan Peternakan\n(Model SDM 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 = )`
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
tm_shape(indo_fixed) +
tm_fill(
"residual_log",
palette = "YlOrRd",
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 = )`
## [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"
## Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.
Visualisasi peta prediksi menunjukkan bahwa Spatial Durbin Model (SDM)
sebagai model terbaik mengidentifikasi Pulau Jawa dan sebagian besar
wilayah Sumatera (khususnya Sumatera Utara dan Sumatera Selatan) sebagai
daerah dengan estimasi PDRB sektor pertanian tertinggi. Hal ini ditandai
dengan dominasi warna kuning hingga hijau terang yang menunjukkan nilai
prediksi berada pada kategori tertinggi (di atas 50 juta/miliar sesuai
satuan data). Konsentrasi nilai tinggi di wilayah barat Indonesia ini
mencerminkan adanya ketergantungan spasial (spatial dependence), di mana
wilayah dengan produktivitas pertanian tinggi cenderung mengelompok
(cluster), kemungkinan besar didukung oleh faktor kesuburan lahan,
aksesibilitas pasar, serta kepadatan tenaga kerja sektor pertanian yang
lebih masif dibandingkan wilayah timur Indonesia yang didominasi warna
ungu gelap (nilai prediksi rendah). Sementara itu, peta residual
memberikan gambaran mengenai tingkat akurasi model dengan memperlihatkan
selisih antara nilai aktual dan hasil prediksi. Wilayah dengan warna
merah gelap (residual positif tinggi), seperti di Aceh, Riau, Jawa
Tengah, dan Sulawesi Selatan, menunjukkan kondisi underprediction,
artinya nilai PDRB pertanian aktual di wilayah tersebut lebih besar
daripada yang mampu diprediksi oleh variabel-variabel dalam model.
Sebaliknya, wilayah berwarna kuning pucat (residual negatif) menunjukkan
overprediction. Secara keseluruhan, persebaran residual yang cukup
variatif dan tidak menunjukkan pengelompokan yang sangat ekstrem di satu
pulau saja mengindikasikan bahwa model SDM telah berhasil menyerap
sebagian besar efek autokorelasi spasial, meskipun masih terdapat
faktor-faktor lokal spesifik di beberapa provinsi yang belum sepenuhnya
tertangkap oleh model global ini.
Interpolasi Spasial
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 1.279562 1.0123780 -0.07991148
## 2 1.25 1.260363 0.9819024 -0.07535527
## 3 1.50 1.261452 0.9696235 -0.06658047
## 4 1.75 1.278032 0.9729874 -0.05625022
## 5 2.00 1.303778 0.9822227 -0.04586438
## 6 2.25 1.333303 0.9982589 -0.03605827
## 7 2.50 1.363033 1.0144459 -0.02701623
## 8 2.75 1.391006 1.0287138 -0.01874601
## 9 3.00 1.416343 1.0403034 -0.01120716# 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.25"# 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)
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: 13.99 - 18.98# 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 PDRB sektor pertanian menggunakan metode Inverse
Distance Weighting (IDW) secara visual mempertegas adanya pola
konsentrasi produktivitas yang dominan di wilayah Barat dan Tengah
Indonesia, menciptakan gradasi spasial yang kontras terhadap wilayah
Timur. Peta tersebut memperlihatkan aglomerasi nilai prediksi tertinggi
(spektrum merah hingga merah tua pada rentang 17.79 – 18.98) yang
tersebar luas di sepanjang koridor Pulau Sumatera, seluruh Pulau Jawa,
serta kantong-kantong produktivitas di Kalimantan Selatan dan Sulawesi
Selatan, yang mencerminkan intensitas aktivitas ekonomi agraris yang
mapan. Sebaliknya, intensitas prediksi mengalami degradasi signifikan
menuju spektrum kuning muda (rentang 13.99 – 15.23) di sebagian besar
wilayah Kalimantan Utara, Maluku, hingga Papua, yang secara visual
mengonfirmasi adanya disparitas output sektoral serta kesenjangan
pengembangan potensi pertanian yang lebar dalam struktur ekonomi
antarwilayah di Indonesia.
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 3.048379 0
## 2 Sph 1.017862 3733300krig_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: 15.31 - 17.86# 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, "Blues", "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 "Blues" is named
## "brewer.blues"
Pemetaan prediksi PDRB sektor pertanian melalui metode Inverse Distance
Weighting (IDW) secara visual mempertegas adanya polarisasi
produktivitas agraris di Indonesia, yang menunjukkan gradasi spasial
mencolok antara wilayah Barat dan Timur. Peta prediksi menampilkan
konsentrasi nilai tertinggi (spektrum merah hingga merah tua pada
rentang 17.79 – 18.98) yang teraglomerasi secara intensif di sepanjang
Pulau Jawa, sebagian besar Sumatera, dan Sulawesi Selatan, mencerminkan
basis ekonomi pertanian yang sangat mapan dengan tingkat ketergantungan
spasial yang tinggi. Sebaliknya, intensitas prediksi mengalami degradasi
tajam menjadi spektrum kuning muda (rentang 13.99 – 15.23) di wilayah
Timur, terutama Papua dan Maluku, yang secara statistik dan visual
mengonfirmasi adanya disparitas output sektoral serta kesenjangan
pengembangan potensi pertanian yang lebar dalam struktur ekonomi
antarwilayah.
cv_results <- data.frame(
Method = c("IDW", "Ordinary Kriging"),
RMSE = c(
sqrt(mean(cv_idw$residual^2)),
sqrt(mean(cv_kriging$residual^2))
),
MAE = c(
mean(abs(cv_idw$residual)),
mean(abs(cv_kriging$residual))
),
ME = c(
mean(cv_idw$residual),
mean(cv_kriging$residual)
)
) %>% arrange(RMSE)
cat("\n=== Cross-Validation Results ===\n")##
## === Cross-Validation Results ===print(cv_results)## Method RMSE MAE ME
## 1 Ordinary Kriging 1.189184 0.9138499 -0.03849947
## 2 IDW 1.257816 0.9908859 -0.06177685Berdasarkan hasil di atas model Oedinary Kriging merupakan metode interpolasi spasial terbaik dalam memodelkan PDRB Pertanian, Perburuan, dan Peternakan tahun 2023.
Bab 5. Kesimpulan
Kesimpulan
Penelitian ini bertujuan untuk menganalisis PDRB sektor pertanian, perburuan, dan peternakan. Hasil penelitian menunjukkan bahwa PDRB sektor pertanian, perburuan, dan peternakan memiliki pola sebaran spasial.
Hasil uji autokorelasi spasial global menggunakan Moran’s I dan Geary’s C menunjukkan adanya autokorelasi spasial positif yang signifikan pada variabel PDRB sektor pertanian, perburuan, dan peternakan. Sementara itu uji autokorelasi spasial loka menggunakan 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 timur. Hal ini menunjukkan pengaruh perbedaan geografis dan iklim yang ada di daerah Indonesia.
Oleh karena itu, dilakukan pendekatan dengan model ekonometrika spasial untuk mengakomodasi ketergantungan antarwilayah. Berdasarkan pertimbangan teoritis dan hasil perbandingan model, diperoleh bahwa Spatial Durbin Model (SDM) merupakan model yang paling sesuai dalam menjelaskan variasi PDRB antarprovinsi di Indonesia dengan nilai AIC terendah.
Model ini menunjukkan bahwa selain pengaruh langsung variabel-variabel independen di suatu provinsi, terdapat pula efek limpahan (spillover effects) yang masif, di mana produktivitas suatu wilayah dipengaruhi secara simultan oleh variabel dependen (spatial lag Y) dan variabel independen (spatial lag X) dari provinsi tetangga. Penggunaan SDM memungkinkan penangkapan dinamika interaksi antarwilayah yang lebih komprehensif, mengasumsikan bahwa perubahan faktor produksi di satu lokasi tidak hanya berdampak secara lokal, tetapi juga memberikan dampak berantai terhadap kinerja ekonomi sektor pertanian di wilayah sekitarnya melalui keterkaitan spasial yang kuat.
Metode interpolasi juga digunakan pada peneltian ini, dengan menggunakan metode yaitu Inverse Distance Weighting (IDW) dan Ordinary Kriging. Hasil menunjukkan bahwa Ordinary Kriging memberikan performa terbaik dengan galat prediksi terkecil. Hasil interpolasi Kriging memperlihatkan adanya tren spasial global yang kuat, di mana PDRB 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 PDRB sektor pertanian, perburuan, dan peternakan merupakan sebuah fenomena spasial.
Saran
saran kebijakan - Mengingat adanya efek limpahan spasial yang signifikan, perencanaan pembangunan sektor pertanian sebaiknya dilakukan secara terkoordinasi antarprovinsi, sehingga peningkatan produktivitas di suatu wilayah dapat memberikan dampak positif bagi wilayah sekitarnya. - Pendekatan berbasis wilayah (spatial-based policy) perlu dikedepankan, terutama dalam pengembangan kawasan pertanian unggulan, agar kebijakan yang diterapkan mampu memperkuat keterkaitan ekonomi antarwilayah.
saran untuk penelitian - Penelitian selanjutnya disarankan menggunakan unit analisis yang lebih detail, seperti tingkat kabupaten/kota, untuk menangkap variasi spasial lokal yang tidak teridentifikasi pada tingkat provinsi. - Penggunaan data panel dan model spasial dinamis dapat dipertimbangkan untuk menganalisis perubahan pola spasial PDRB sektor pertanian dari waktu ke waktu serta mengevaluasi dampak kebijakan secara lebih komprehensif. - Disarankan untuk menambahkan variabel-variabel penjelas yang lebih komprehensif, seperti luas lahan pertanian, produktivitas komoditas utama, tingkat mekanisasi pertanian, akses irigasi, infrastruktur logistik, serta variabel iklim (curah hujan dan suhu), guna meningkatkan kemampuan model dalam menjelaskan variasi PDRB antarwilayah.
Daftar Pustaka
[1] Y. Ru, B. Blankespoor, U. Wood-Sichra, T. S. Thomas, L. You, and E. Kalvelagen, “Estimating local agricultural gross domestic product (AgGDP) across the world,” Earth Syst. Sci. Data, vol. 15, no. 3, pp. 1357–1387, Mar. 2023, doi: 10.5194/essd-15-1357-2023.
[2] D. Xu, “Spatial spillover and threshold effects of digital rural development on agricultural circular economy growth,” Front. Sustain. Food Syst., vol. 8, Art. no. 1337637, Mar. 2024, doi: 10.3389/fsufs.2024.1337637.
[3] G. Yang, Y. Zhao, and L. Wang, “The Regional Heterogeneity of the Impact of Agricultural Market Integration on Regional Economic Development: An Analysis of Pre-COVID-19 Data in China,” Sustainability, vol. 16, no. 5, p. 1734, Feb. 2024, doi: 10.3390/su16051734.
[4] R. F. Sari, “Analisis Pengaruh PDRB Sektor Pertanian, Kehutanan dan Perikanan Terhadap Penyerapan Tenaga Kerja di Provinsi Kalimantan Selatan,” ECOPLAN, vol. 3, no. 2, pp. 69–75, Oct. 2020, doi: 10.20527/ecoplan.v3i2.69.
[5] M. Y. Ariefqi and D. Ariyani, “The Effect Of Agricultural Sector Labor, Agricultural Sector Exports, And Agricultural Sector Investment On Economic Growth With Technology As A Moderating Variable,” J. Econ. Bus. Aseanomics (JEBA), vol. 10, no. 1, pp. 30–41, 2022, doi: 10.33476/jeba.v10i1.5478.
[6] D. A. Griffith, “Understanding Spatial Autocorrelation: An Everyday Metaphor and Additional New Interpretations,” Geographies, vol. 3, no. 3, pp. 543–562, Aug. 2023, doi: 10.3390/geographies3030028.
[7] L. Pregi and L. Novotný, “Spatial Autocorrelation Methods in Identifying Migration Patterns: Case Study of Slovakia,” Appl. Spat. Anal. Policy, vol. 17, no. 1, pp. 249–272, Mar. 2024, doi: 10.1007/s12061-024-09615-5.
[8] P. M. Robinson and F. Rossi, “Refined tests for spatial correlation,” Econometric Theory, vol. 31, no. 6, pp. 1249–1280, Dec. 2015, doi: 10.1017/S0266466614000498.
[9] M. M. Fischer and J. P. LeSage, “Spatial Econometrics,” Working Papers in Regional Science, no. 2023/01, 2023, doi: 10.57938/1019c1c5-41ba-4ffc-a61f-0c645889aa3b.
[10] W. Yamaka and P. Shi, “Spatial Spillover Effects of Internet Development on Foreign Trade in China,” Sustainability, vol. 15, no. 5, p. 4213, Feb. 2023, doi: 10.3390/su15054213.
[11] J. P. LeSage and R. K. Pace, “Bayesian estimation and model selection for spatial Durbin error models,” Reg. Sci. Urban Econ., vol. 43, no. 4, pp. 606–613, Jul. 2013, doi: 10.1016/j.regsciurbeco.2013.04.006.
[12] T. Juniarsih, “Pengaruh Produk Domestik Regional Bruto (PDRB) Sektor Pertanian di Aceh terhadap Produk Domestik Regional Bruto (PDRB) Provinsi Aceh,” VALUE, vol. 2, no. 1, pp. 29–44, 2021, doi: 10.36490/value.v2i1.119.
[13] M. Khotiawan, R. K. Sakti, and S. T. Wahyudi, “An Analysis of the Effects of Spatial Dependence on Economic Growth Among Regencies and Cities in Java,” Jurnal Ekonomi Pembangunan, vol. 24, no. 2, pp. 90–100, 2023, doi: 10.23917/jep.v24i2.22109.
[14] Z. Yin and J. Wu, “Spatial dependence evaluation of agricultural technical efficiency—based on the stochastic frontier and spatial econometric model,” Sustainability, vol. 13, no. 5, p. 2708, 2021, doi: 10.3390/su13052708.
[15] Y. Chen, “An analytical process of spatial autocorrelation functions based on Moran’s index,” PLOS ONE, vol. 16, no. 4, p. e0249589, Apr. 2021, doi: 10.1371/journal.pone.0249589.
[16] Rahmawati, D. Safitri, and O. U. Fairuzdhiya, “Analisis Spasial Pengaruh Tingkat Pengangguran terhadap Kemiskinan di Indonesia (Studi Kasus Provinsi Jawa Tengah),” Media Statistika, vol. 8, no. 1, pp. 23–30, Jun. 2015, doi: 10.14710/medstat.8.1.23-30.
[17] Y. Chen, “Spatial autocorrelation equation based on Moran’s index,” Sci. Rep., vol. 13, Art. no. 19296, Nov. 2023, doi: 10.1038/s41598-023-45947-x.
[18] H. Yamada, “Geary’s c for Multivariate Spatial Data,” Mathematics, vol. 12, no. 12, p. 1820, Jun. 2024, doi: 10.3390/math12121820.
[19] J. A. Taylor, B. Tisseyre, C. Jain, B. J. Lawrence, and A. Getis, “Getis–Ord’s hot- and cold-spot statistics as a basis for multivariate spatial clustering of orchard tree data,” Comput. Electron. Agric., vol. 111, pp. 20–27, Feb. 2015, doi: 10.1016/j.compag.2014.12.011.
[20] L. Anselin, “Local Indicators of Spatial Association—LISA,” Geogr. Anal., vol. 27, no. 2, pp. 93–115, Apr. 1995, doi: 10.1111/j.1538-4632.1995.tb00338.x.
[21] J. P. LeSage, “Spatial econometrics,” in Handbook of Research Methods and Applications in Economic Geography, C. Tammaru et al., Eds. Cheltenham, UK: Edward Elgar Publishing, 2021, pp. 276–291, doi: 10.4337/9781789903942.00014.
[22] J. S. Clark, “Why environmental data are Poisson: some diagnostics and strategies,” Ecol. Monogr., vol. 88, no. 2, pp. 1–15, May 2018, doi: 10.1002/ecm.1283.
[23] V. Yildirim and Y. M. Kantar, “Robust estimation approach for spatial error model,” J. Stat. Comput. Simul., vol. 90, no. 9, pp. 1602–1619, Jun. 2020, doi: 10.1080/00949655.2020.1740223.
[24] L. F. Lee and J. Yu, “Identification of Spatial Durbin Panel Models,” J. Appl. Econ., vol. 31, no. 1, pp. 133–162, Jan. 2016, doi: 10.1002/jae.2450.
[25] J. P. LeSage, N. J. Yang, and S. Fang, “Energy efficiency in the Chinese provinces: a fixed effects stochastic frontier approach,” Ann. Reg. Sci., vol. 59, no. 2, pp. 437–454, Sep. 2017, doi: 10.1007/s00168-016-0782-5.
[26] G. Gaspard, D. Kim, and Y. Chun, “Residual spatial autocorrelation in macroecological and biogeographical modeling: a review,” J. Ecol. Environ., vol. 43, Art. no. 19, Jun. 2019, doi: 10.1186/s41610-019-0118-3.
[27] A. Getis, “The Moran’s I Index,” in Handbook of Applied System Analysis, M. M. Fischer and P. Nijkamp, Eds. Cham, Switzerland: Springer, 2019, pp. 1–18, doi: 10.1007/978-981-13-8664-0_2.
[28] Z. Li, X. Zhang, R. Zhu, Z. Zhang, and Z. Weng, “Integrating data-to-data correlation into inverse distance weighting,” Comput. Geosci., vol. 24, no. 2, pp. 203–216, Apr. 2020, doi: 10.1007/s10596-019-09913-9.
[29] I. Mesić Kiš, “Comparison of Ordinary and Universal Kriging interpolation techniques on a depth variable (a case of linear spatial trend), case study of the Šandrovac Field,” Rud.-geol.-naft. zb., vol. 31, no. 2, pp. 41–58, May 2016, doi: 10.17794/rgn.2016.2.4.