Link Dashboard: https://chrisrafael.shinyapps.io/Dashboard_Spasial_IKLH/
Link Video: https://youtu.be/Ld4Qj07T6o8
Lingkungan hidup memiliki banyak kontribusi terhadap pembangunan berkelanjutan dan kualitas hidup masyarakat. Di Indonesia, tekanan terhadap lingkungan semakin dirasakan dengan peningkatan ekonomi, pertambahan populasi, dan pemanfaatan sumber daya alam yang lebih intensif. Masalah seperti deforestasi, pencemaran air dan udara, degradasi lahan, dan dampak perubahan iklim masih menjadi tantangan utama dalam pengelolaan lingkungan (Setyawan & Primanandha, 2022). Indeks Kualitas Lingkungan Hidup (IKLH/EQI) dikembangkan sebagai indikator komposit yang menggabungkan kualitas udara, udara, dan lahan, dan kini menjadi alat utama evaluasi kinerja pembangunan lingkungan lintas provinsi, termasuk dasar bagi penyusunan kebijakan pengelolaan dan komunikasi kondisi lingkungan kepada masyarakat (Wafa, 2025). IKLH telah menjadi alat penting dalam menilai sejauh mana pembangunan telah dicapai sejalan dengan perlindungan lingkungan. Hal tersebut didukung oleh penelitian oleh Perwithosuci dkk. (2025) yang menggunakan IKLH untuk menilai keterkaitan pertumbuhan ekonomi, indeks pembangunan manusia, kepadatan penduduk, dan tutupan hutan sebagai variabel independen. Dengan demikian, IKLH berperan sebagai indikator komprehensif untuk menilai keseimbangan antara pembangunan ekonomi, dinamika sosial, dan keberlanjutan lingkungan, sekaligus menjadi dasar penting dalam penyusunan kebijakan.
Struktur PDRB sektoral berperan penting dalam membentuk kualitas lingkungan. Sektor-sektor agraris seperti pertanian, kehutanan, dan perikanan dapat memberikan kontribusi positif apabila dikelola secara berkelanjutan, namun ekspansi yang tidak terkendali, termasuk konversi hutan menjadi lahan pertambangan atau perkebunan, terbukti menurunkan indeks kualitas tutupan lahan dan memicu deforestasi (Setyawan & Primanandha, 2022). Pengelolaan lingkungan yang tidak efisien akan mengancam kualitas lingkungan. Namun, penerapan praktik pembangunan berkelanjutan memiliki potensi untuk meningkatkan dan memperbaiki IKLH. Sebaliknya, sektor pertambangan dan industri pengolahan telah terkait dengan penurunan kualitas lingkungan akibat aktivitas mereka yang menghasilkan emisi dan limbah yang tinggi (Listiyani, 2017) Perbedaan dalam sifat dan intensitas dampak dari sektor-sektor ini menunjukkan bahwa kontribusi ekonomi dari sektor-sektor ini, khususnya sektor pertambangan, berpotensi memengaruhi kualitas lingkungan dengan cara yang berbeda.
Kondisi sosial-ekonomi masyarakat juga menentukan kualitas lingkungan. Berbagai penelitian menemukan bahwa kemiskinan memiliki efek negatif terhadap IKLH/EQI karena rumah tangga miskin lebih bergantung pada eksploitasi langsung sumber daya alam dan cenderung kurang memiliki pilihan mata pencaharian alternatif (Sumargo dkk., 2020). Sebaliknya, peningkatan kualitas sumber daya manusia yang tercermin dalam Indeks Pembangunan Manusia maupun indikator pendidikan berpengaruh positif dan signifikan terhadap kualitas lingkungan di banyak provinsi, antara lain melalui peningkatan kesadaran, perilaku pro-lingkungan, dan kemampuan mengadopsi teknologi yang lebih efisien (Rahayu & Handri, 2023). Hubungan antara pertumbuhan ekonomi dan kualitas lingkungan di Indonesia juga sering dibahas dalam kerangka Environmental Kuznets Curve (EKC), di mana pada tahap awal industrialisasi, pertumbuhan cenderung meningkatkan degradasi lingkungan, namun pada tingkat pendapatan tertentu, penguatan kebijakan lingkungan dan penerapan teknologi bersih dapat memperbaiki kualitas lingkungan (Prasetyanto & Sari, 2021). Dengan demikian, analisis determinan IKLH di Indonesia perlu memasukkan dinamika pertumbuhan, kemiskinan, pendidikan, dan struktur ekonomi secara simultan.
Penelitian ini dilakukan untuk mengetahui bagaimana faktor sosio-ekonomi yaitu PDRB Sektor Perhutanan, Pertanian, dan Perikanan (X1), PDRB Sektor Pertambangan (X2), PDRB Industri Pengolahan (X3), Persentase Kemiskinan (X4), Rata Lama Sekolah (X5), dan Pertumbuhan Ekonomi (X6) pada lingkungan hidup Indonesia yang direpresentasikan Indeks Kualitas Lingkungan Hidup (Y). Dengan menangani adanya faktor spasial yang terjadi pada variabel dependen, independen, ataupun keduanya, penelitian ini diharapkan bisa menjadi fondasi statistik untuk penanganan lingkungan hidup di negara Indonesia serta pengelolaan perekonomian dan sosial Indonesia untuk bisa berkelanjutan dan ramah lingkungan.
Riset terfokus pada wilayah Indonesia yang mencangkupi 34 Provinsi pada tahun 2022, dimana data juga terambil pada tahun yang sama. Provinsi tersebut adalah Aceh, Sumatera Utara, Sumatera Barat, Riau, Kepulauan Riau, Jambi, Sumatera Selatan, Kepulauan Bangka Belitung, Bengkulu, Lampung, DKI Jakarta, Jawa Barat, Jawa Tengah, Daerah Istimewa Yogyakarta, Jawa Timur, Banten, Bali, Nusa Tenggara Barat, Nusa Tenggara Timur, Kalimantan Barat, Kalimantan Tengah, Kalimantan Selatan, Kalimantan Timur, Kalimantan Utara, Sulawesi Utara, Sulawesi Tengah, Sulawesi Selatan, Sulawesi Tenggara, Gorontalo, Sulawesi Barat, Maluku, Maluku Utara, Papua Barat, Papua. Indikator sosio-ekonomi yang digunakan adalah PDRB Sektor Perhutanan, Pertanian, dan Perikanan (X1), PDRB Sektor Pertambangan (X2), PDRB Industri Pengolahan (X3), Persentase Kemiskinan (X4), Rata Lama Sekolah (X5), dan Pertumbuhan Ekonomi (X6) yang menggunakan Indeks Kualitas Lingkungan Hidup (Y) sebagai variabel dependen. Semua data diambil pada tahun 2022. Kapita.
Dependensi spasial menyatakan bahwa adanya kebergantungan pada data, terutama yang mempunyai status geografi, saling bergantung dengan variabel yang sama pada lokasi atau area yang berbeda. Dibandingkan dengan sebuah data cross-sectional dimana asumsi bahwa semua data independen (Anselin, 1988). Data spasial yang merupakan data yang diambil atau diobservasi berdasarkan lokasi pada sebuah space atau ruang. Dan untuk ruang tersebut bisa berupa sebuah titik, garis, atau berbentuk polygon seperti daerah administratif negara (Anselin, 1988).
Bisa disimpulkan bahwa sebuah data spasial yang memiliki kejauhan yang cukup dekat, tetangga, sering kali memiliki karakteristik yang sama karena kedekatan tersebut (ceteris paribus). Bisa diambil sebuah pola berdasarkan kedekatan yang bisa menjelaskan hubungan antara dua unit spasial, seperti hot-hot atau hot-cold (Bivand et al., 2008).
Untuk bisa melakukan analisis berdasarkan autokorelasi spasial, diperlukan sebuah spatial weight. Spatial weight merupakan sebuah matrix untuk menjelaskan hubungan antara unit spasial, serta seberapa besar pengaruh satu sama lain. Pada bentuk paling dasar, spatial weight bisa berbentuk matrix bobot, dimana 1 menjelaskan hubungan tetangga dan 0 menunjukkan tidak. Dan dari hal tersebut didapatkan banyak model weight seperti biner ataupun uniform.
Dan untuk melihat apakah adanya autokorelasi spasial maka diperlukan beberapa uji. Pada dasarnya ada tipe uji untuk melihat autokorelasi spasial pada sebuah data, sebuah uji untuk melihat semua hubungan spasial dan melihat apakah ada signifikansinya, serta uji lokal untuk melihat hubungan secara lokal dan membuat kluster berdasarkan pengujian tersebut (Bivand et al., 2008).
Untuk pemodelan Ekonometrika pada spasial, memang sebenarnya merupakan ilmu cabang dari ekonomterika, ilmu statistika untuk melakukan estimasi model serta menguji hubungan kausal antara fenomena terutama dalam bidang ekonomi. Dengan adanya spasialitas, estimasi model tersebut juga akan ditambahkan dengan penjelasan interdependensi spasial di sebuah model spasial serta juga estimasi hubungan spasial yang bersifat interregional. Pada model-model spasial ekonometrika juga terjadi spesifikasi modelm estimasi, serta juga testing hipotesa namun dengan tambahan adanya representasi dependensi spasial serta juga heterogenitas spasial. Model-model yang akan digunakan pada spasial ekonometrika ada bermacam-macam (Bivand et al., 2008). Digunakannya pada penelitian saat ini untuk menghandalkan efek dependensi spasial pada variabel model.
Interpolasi spasial merupakan sebuah metode pada analisis spasial yang digunakan untuk melakukan estimasi pada sebuah variabel yang tidak diobserved (Comber, 2019) (Ikechukwu, 2017). Pada dasarnya interpolasi bias dibagi menjadi dua jenis, yaitu interpolasi titik dan interpolasi area. Dimana interpolasi titik melakukan estimasi pada lokasi berdasarkan data titk yang sudah diketahui, dimana metode yang umum digunakan adalah Nearest Neighbour, Inverse Distance Weight, Thin-Spliced, dan Kriging (Comber, 2019). Walaupun, areal merupakan estimasi berdasarkan data area. Interpolasi spasial sering kali digunakan untuk bisa melakukan estimasi, terutama pada sistem GIS.
Data yang digunakan pada penelitian, dimana dibagikan menjadi variabel respon (Y) dan variabel prediktor (X) dalam pembuatan model, adalah sebagai berikut:
| Variabel | Nama Variabel | Indikator | Sumber |
|---|---|---|---|
| Y | Indeks Kualitas Lingkungan Hidup (2022) | Kualitas lingkungan hidup di sebuah wilayah | BPS (Badan Pusat Statistika) |
| X1 | PDRB Sektor Perhutanan, Pertanian, dan Perikanan (2022) | Pembangunan ekonomi wilayah | BPS (Badan Pusat Statistika) |
| X2 | PDRB Sektor Pertambangan (2022) | Pembangunan ekonomi wilayah | BPS (Badan Pusat Statistika) |
| X3 | PDRB Industri Pengolahan (2022) | Pembangunan ekonomi wilayah | BPS (Badan Pusat Statistika) |
| X4 | Persentase Kemiskinan (2022) | Indikator sosial provinsi | BPS (Badan Pusat Statistika) |
| X5 | Rata Lama Sekolah (2022) | Indikator sosial provinsi | BPS (Badan Pusat Statistika) |
| X6 | Laju Pertumbuhan Ekonomi (2022) | Pembangunan ekonomi wilayah | BPS (Badan Pusat Statistika) |
Unit spasial yang akan digunakan berupa semua provinsi Indonesia pada tahun 2022 sebelum pengesahan provinsi baru pada Papua. Provinsi tersebut adalah Aceh, Sumatera Utara, Sumatera Barat, Riau, Kepulauan Riau, Jambi, Sumatera Selatan, Kepulauan Bangka Belitung, Bengkulu, Lampung, DKI Jakarta, Jawa Barat, Jawa Tengah, Daerah Istimewa Yogyakarta, Jawa Timur, Banten, Bali, Nusa Tenggara Barat, Nusa Tenggara Timur, Kalimantan Barat, Kalimantan Tengah, Kalimantan Selatan, Kalimantan Timur, Kalimantan Utara, Sulawesi Utara, Sulawesi Tengah, Sulawesi Selatan, Sulawesi Tenggara, Gorontalo, Sulawesi Barat, Maluku, Maluku Utara, Papua Barat, Papua.
Metode analisis yang akan dilakukan akan berupa dari metode untuk mendeteksi autokorelasi spasial pada variabel respon serta prediktor sebagai referensi yang digunakan untuk menentukan model spasial ekonometrika yang tepat dengan metode Global Moran’s I, Local Moran’s I, Getis G, dan Geary’s C. Kemudian akan dilakukan pemodelan spasial ekonometrika berdasarkan sembilan model utama kemudian dilakukan penentuan model terbaik dengan LM test, perbandingan AIC, dan teoritis. Dimana digunakan untuk menghandalkan eksistensinya dependensi spasial pada variabel model. Kemudian dilakukan analisis interpolasi spasial pada masing-masing variabel menggunakan metode interpolasi Nearest Neighbor, Inverse Distance Weighting, dan Kriging.
Uji global Moran’s I merupakan uji global yang lebih banyak dipakai. Dimana Moran’s I digunakan untuk melihat seluruh hubungan spasial autokorelasi pada keseluruhan area. Berdasarkan tersebut didapatkan rumus Moran’s I sebagai berikut (Bivand et al., 2008):
\[I = \ \frac{N}{S_{0}}\ .\ \frac{\sum_{i = 1}^{N}{\sum_{j = 1}^{N}\omega_{ij}}(x_{i} - \overline{x})(x_{j} - \overline{x})}{\sum_{i = 1}^{N}{(x_{i} - \overline{x})}^{2}}\]
Dimana, N merupakan jumlah data yang ada, \(\omega_{ij}\) sebagai weight dan x merupakan datanya tersendiri. Dan untuk interpretasi adalah sebagai berikut:
I > 0: Autokorelasi spasial positif, data berkarakteristik sama cenderung saling berdekatan
I < 0: Autokorelasi spasial negatif, data berkarakteristik berbeda cenderung saling berdekatan
I ≈ 0: Pola acak
Untuk uji global kedua adalah uji Geary’s C, sebuah uji global untuk melihat autokorelasi spasial, namun dengan tambahan bahwa perhitungan menggunakan weighted sum. Dibandingkan dengan Moran, Geary lebih sensitif terhadap autokorelasi bersifat lokal. Didapatkan perhitungan sebagai berikut (Yamada, 2024):
\[C = \ \frac{N - 1\ \sum_{i = 1}^{N}{\sum_{j = 1}^{N}\omega_{ij}}{(x_{i} - x_{j})}^{2}}{2S_{0}\sum_{i = 1}^{N}{(x_{i} - \overline{x})}^{2}}\ \]
Dimana, N merupakan jumlah data yang ada, \(\omega_{ij}\) sebagai weight dan x merupakan datanya tersendiri. Dan untuk interpretasi adalah sebagai berikut:
C < 1: Autokorelasi spasial positif, data berkarakteristik sama cenderung saling berdekatan
C > 1: Autokorelasi spasial negatif, data berkarakteristik berbeda cenderung saling berdekatan
C = 1: Pola acak
Untuk uji lokal yang akan digunakan merupakan varian dari Moran’s I namun diberlakukan untuk uji lokal, dimana akan diperhitungkan berdasarkan dari hubungan lokal antar tetangga. Didapatkan perhitungan sebagai berikut (Bivand et al., 2008):
\[I_{i} = \frac{(x_{i} - \overline{x})\sum_{j = 1}^{N}\omega_{ij}(x_{j} - \overline{x})}{\frac{\sum_{i = 1}^{N}{(x_{i} - \overline{x})}^{2}}{n}}\]
Dimana, N merupakan jumlah data yang ada, \(\omega_{ij}\) sebagai weight dan x merupakan datanya tersendiri dan n sebagai unit spasial. Dan untuk interpretasi adalah sebagai berikut:
I > 0: Autokorelasi spasial positif, data berkarakteristik sama cenderung saling berdekatan
I < 0: Autokorelasi spasial negatif, data berkarakteristik berbeda cenderung saling berdekatan
I ≈ 0: Pola acak
Uji lokal selanjutnya merupakan Getis-Ord. Sebuah uji lokal yang dibandingkan dengan moran dimana Getis-Ord lebih sensitif terhadap perbedaan nilai. Didapatkan hasil (Chen, 2020):
\[G = \ \frac{\sum_{i = i}^{n}{\sum_{j = i}^{n}\omega_{ij}x_{i}x_{j}}}{\sum_{i = i}^{n}{x_{i}x_{j}}}\]
Dimana, N merupakan jumlah data yang ada, \(\omega_{ij}\) sebagai weight dan x merupakan datanya tersendiri dan n sebagai unit spasial.
OLS (Ordinary Least Squares) digunakan untuk membuat inferensi regresi dari data dan menghasilkan sebuah model yang bersifat BLUE (Best Linear Unbiased Estimator). Dengan OLS, tujuan utamanya adalah meminimalkan jumlah kuadrat error untuk mendapatkan estimasi parameter yang paling optimal. Untuk menghitung serta mendapatkan koefisien regresi melalui OLS, digunakan rumus sebagai berikut (Verbeek, 2004):
\[\mathbf{\beta =}\left( \mathbf{X`X} \right)^{\mathbf{- 1}}\mathbf{X`y,}\]
Rumus untuk estimasi koefisien menggunakan matrix, didapatkan model:
\[\mathbf{y = X\beta + e}\]
OLS merupakan bentuk terdasar dari model ekonometrika dan juga yang sangat sederhana. Pada riset ini jika memang tidak terbuktinya adanya spasial interdependensi pada modelm maka akan digunakan OLS serta menjadi baseline.
Model spasial dimana diasumsikan bahwa ada spasial interdependensi yang terdapat pada variabel prediktornya (X) dan dimana X dan W adalah non-stokastik. Model merupakan model tersederhana pada spasial ekonometrika dalam memodelkan spillover pada x, serta bisa digunakan estimasi OLS. didapatkan model seperti berikut (Arbia, 2014):
\[y = WX\gamma + X\beta + u\]
Model spasial dimana diasumsikan bahwa ada spasial interdependensi yang terdapat pada variabel respon (Y). Model terdapat endogeneitas antara y dan gangguan stokastik, sehingga dibutuhkan estimasi MLE untuk mendapatkan hasil serta remodeling. Didapatkan model (Arbia, 2014):
\[y = \ \lambda Wy + Z\beta + u\]
Untuk mencegah korelasi y:
\[y = {(I - \lambda W)}^{- 1}Z\beta + (I - \lambda W)\]
Model spasial dimana diasumsikan bahwa ada spasial heteroskedastisitas dan adanya korelasi dengan error karena adanya omitted variabel. Didapatkan model (Arbia, 2014):
\[y = Z\beta + u\]
\[u = \rho Wu + \varepsilon\ |\rho| < 1\]
Varian dari Spatial Lag Model, dimana selain adanya interdependensi spasial pada variabel respon model, juga terdapat interdepensi yang muncul pada variabel ekpslanatory, namun model sangat rawan terhadap adanya omitted variabel sehingga sering kali diturunkan menjadu model SEM. Didapatkan model sebagai berikut (Arbia, 2014):
\[y = \lambda Wy + WX\gamma + X\beta + u\]
Varian dari Spatial Durbin Model, dimana tidak adanya interdependensi spasial pada variabel respon model, namun tetap terdapat interdepensi yang muncul pada variabel ekpslanatory serta adanya omitted variabel. Didapatkan model sebagai berikut (Kholifia et al., 2021):
\[y = X\beta + WX\theta + u\]
\[u = \ \lambda Wu + \ \varepsilon,\ \varepsilon\sim N(0,\sigma^{2})\]
Model spasial dimana diasumsikan bahwa terdapat interdependensi spasial pada variabel respon serta adanya omitted variabel seperti pada model SEM. Didapatkan model (Anselin, 1988):
\[y = \ \rho Wy + X\beta + \ \epsilon,\]
\[\epsilon = \lambda W\epsilon + \ \mu\]
Model dimana digunakan saat terdapat interdependensi spasial pada variabel respon (Y) dan variabel prediktor (X) serta juga terdapat omitted variabel yang menyebabkan korelasi antar error di model. Model mencangkup semua permasalahan yang muncul saat melakukan pemodelan spasial. Model adalah sebagai berikut (Elhorst, 2022):
\[Y_{t} = \ \tau Y_{t - 1} + \rho WY_{t} + \eta WY_{t - 1} + X_{t}\beta + WX_{t}\theta + \sum_{r}^{}\Gamma_{r}^{T}f_{rt} + u_{t}\]
\[u_{t} = \lambda Wu_{t} + \varepsilon_{t}\]
Metode yang digunakan untuk membandingkan dua model, terutama sangat dipakai pada Ekonomentrika untuk penentuan model yang sangat efektif. Didapatkan statistik sebagai berikut (Astaiza-Gomez, 2020):
\[\xi^{LM} = \frac{1}{n}{\widehat{\lambda}}^{`}g\left( \widehat{\theta} \right)`I\left( \widehat{\theta} \right)^{- 1}{\widehat{\lambda}}^{`}g\left( \widehat{\theta} \right)\ \sim\ \chi_{r}^{2}\]
Metode yang digunakan untuk membandingkan dua model seperti Lagrange Multiplier Test, digunakan oleh peneliti dalam bidang Ekonometrika. Didapatkan statistik sebagai berikut (Bai et al., 2024):
\[LR = \ - k\log\left( \left| {\widehat{\Sigma}}_{1}(k) \right| \right) - (T - k)log(\left| {\widehat{\Sigma}}_{2}(k) \right|)\] ### NEAREST NEIGHBOR
Metode interpolasi yang mengukur estimasi berdasarkan kedekatan dengan tetangga lain menggunakan perhitungan jarak.
Metode interpolasi dengan melakukan estimasi sebagai rata-rata berbobot dari titik yag diketahui berdasarkan sebuah jarak tertentu. Didefinisikan jarakj sebagai proporsi terhadap jarak. Didapatkan estimator sebagai berikut (Ikechukwu, 2017):
\[ F(s) = \sum_{i=1}^{n} w_i \, z(s_i)\; \frac{\displaystyle \sum_{k=1}^{m} \frac{z(s_k)}{\lvert s - s_k \rvert^{p}}} {\displaystyle \sum_{j=1}^{m} \frac{1}{\lvert s - s_j \rvert^{p}}} \]
\[ \lambda_i = \frac{\dfrac{1}{d_i^{p}}}{\displaystyle \sum_{i=1}^{n} \frac{1}{d_i^{p}}} \]
Metode interpolasi yang menggunakan fungsi spline, untuk melakukan fit yang meminimalkan kelengkungan total permukaan berdasarkan beberapa poin. Didapatkan estimator sebagai berikut (Ikechukwu, 2017):
\[ f(s) = a_1 + a_2 x + a_3 y + \sum_{i=1}^{n} w_i \, z\bigl(\lvert s_i - s_0 \rvert\bigr) \]
Metode interpolasi spasial berdasarkan variabel yang terobservasi. Dimana Kriging diperlukannya melakukan pemodelan semivariogram. Dimana variabel akan diberlakukan sebagai probabilitas. Estimasi Kriging merupakan estimasi BLUE (Best Linear Unbiased Estimation). Dari Kriging juga terdapat tiga jenis metode. Yaitu Simple Kriging, Ordinary Kriging, dan Universal Kriging. Dimana akan digunakan Ordinary dan Universal (Zaresefat, 2019).
Ordinary Kriging mengasumsi bahwa rata-rata kontan dari unit spasial ke unit spasial. Dan metode tersebut bergantung pada adanya dependensi spasial untuk melakukan estimasi. Dengan estimasi sebagai berikut (Zaresefat, 2019):
\[ z(x_0) = m\left(1 - \sum_{i=1}^{n} \lambda_i \right) + \sum_{i=1}^{n} \lambda_i \, z(x_i) \]
Universal Kriging digunakan untuk data spasial yang terdapat sebuah global tren dengan melakukan detrend pada data. Serta digunakan juga model derajat pertama atau kedua polinomial. Didapatkan estimasi sebagai berikut, untuk derajat pertama diatas dan derajat kedua dibawah (Zaresefat, 2019):
\[ m(x_0) = a_0 + \sum_{i=1}^{k} \left( a_1 z_{1,i} + a_2 z_{2,i} \right) \]
\[ m(x_0) = a_0 + \sum_{i=1}^{k} \left( a_1 z_{1,i} + a_2 z_{2,i} + a_3 z_{1,i}^{2} + a_4 z_{2,i}^{2} + a_5 z_{1,i} z_{2,i} \right) \]
Penelitian dimulai dengan pengumpulan data-data sekunder yang akan diestimasi menjadi model. Analisis pada program lunak Rstudio. Pertama-tama dilakukan analisis deskriptif dalam bentuk perhitungan serta juga visualisasi untuk menggali identitas dari masing-masing variabel pada model. Kemudian bisa diikuti dengan penentuan bobot yang akan digunakan pada model kedepannya. Dimana untuk bobot yang digunakan merupakan Row Based Proportion. Kemudian dilakukan pengujian autokorelasi menggunakan metode yang dibahas sebelumnya. Kemudian dilakukan spesifikasi model menggunakan data yang telah distandardized pada kesembilan model ekonometrika. Pada tahap terakhir berupa penentuan model menggunakan pembandingan AIC kemudian LM Test dan LR Test pada OLS dan Nested. Kemudian dilakukan analisis interpolasi spasial dengan mengubah data menjadi data titik menurut sentroid pada masing-masing provinsi, dimana akan dilakukan spesifikasi berdasarkan empat metode interpolasi kemudian dianalisis yang paling optimal.
Analisis deskriptif pada semua variabel pada model mendapatkan hasil sebagai berikut:
# Package
library(raster)
## Loading required package: sp
library(gstat)
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')`
## Loading required package: sf
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(sp)
library(sf)
library(FNN)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(spData)
library(geodata)
## Loading required package: terra
## terra 1.8.86
library(skimr)
##
## Attaching package: 'skimr'
## The following object is masked from 'package:raster':
##
## bind
library(psych)
##
## Attaching package: 'psych'
## The following objects are masked from 'package:terra':
##
## describe, distance, rescale
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## The following object is masked from 'package:raster':
##
## distance
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(grid)
##
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
##
## depth
library(viridis)
## Loading required package: viridisLite
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.ZeroPolicyOption
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
##
## time<-
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
## The following object is masked from 'package:terra':
##
## extract
## The following object is masked from 'package:raster':
##
## extract
library(lwgeom)
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.13.1, PROJ 9.7.0
##
## Attaching package: 'lwgeom'
## The following objects are masked from 'package:sf':
##
## st_minimum_bounding_circle, st_perimeter
library(fields)
## Loading required package: spam
## Spam version 2.11-1 (2025-01-20) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following object is masked from 'package:Matrix':
##
## det
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: RColorBrewer
##
## Try help(fields) to get started.
##
## Attaching package: 'fields'
## The following object is masked from 'package:psych':
##
## describe
## The following object is masked from 'package:geodata':
##
## world
## The following object is masked from 'package:terra':
##
## describe
# Set up Data
setwd("C:\\Users\\HP\\Documents\\KULIAH\\SEMESTER 5\\Spatial\\Dashboard_UAS")
indo_adm1 <- st_read("IDN_AdminBoundaries_candidate.gdb",
layer = "idn_admbnda_adm1_bps_20200401")
## Reading layer `idn_admbnda_adm1_bps_20200401' from data source
## `C:\Users\HP\Documents\KULIAH\SEMESTER 5\Spatial\Dashboard_UAS\IDN_AdminBoundaries_candidate.gdb'
## using driver `OpenFileGDB'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts. The
## processing may be really slow. You can skip the processing by setting
## METHOD=SKIP.
## Simple feature collection with 34 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
data <- (read.csv("Data_UAS.csv", sep=";", dec="."))
indo_adm1$id <- c(1:34)
indo_sf <- st_as_sf(indo_adm1)
indo_merged <- indo_sf %>%
left_join(data, by = "id")
# Explorasi Data
## Peta Indeks Kualitas Lingkungan Hidup
ggplot(data = indo_merged) +
geom_sf(aes(fill = IKLH..Y.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(
option = "turbo",
direction = -1,
name = "Indeks Kualitas Lingkungan Hidup (IKLH)"
) +
labs(
title = "Indeks Kualitas Lingkungan Hidup (IKLH) 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta PDRB Sektor Hutan
ggplot(data = indo_merged) +
geom_sf(aes(fill = PDRB.Perhutanan.Pertanian.Perikanan..X1.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "PDRB Sektor Perhutanan/Pertanian/Perikanan") +
labs(
title = "PDRB Sektor Perhutanan/Pertanian/Perikanan 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta PDRB Sektor Pertambangan
ggplot(data = indo_merged) +
geom_sf(aes(fill = PDRB.Pertambangan..X2.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "inferno", direction = -1,
name = "PDRB Sektor Pertambangan") +
labs(
title = "PDRB Sektor Pertambangan 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta PDRB Industri Pengolahan
ggplot(data = indo_merged) +
geom_sf(aes(fill = PDRB.Industri.Pengolahan..X3.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "cividis", direction = -1,
name = "PDRB Sektor Industri") +
labs(
title = "PDRB Sektor Industri Pengolahan 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta Persentase Kemiskinan
ggplot(data = indo_merged) +
geom_sf(aes(fill = Persentase.Kemiskinan..X4.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "rocket", direction = -1,
name = "Persentase Populasi Miskin") +
labs(
title = "Persentase Populasi Miskin 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta RLS
ggplot(data = indo_merged) +
geom_sf(aes(fill = RLS..X5.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Rata Lama Sekolah") +
labs(
title = "Rata Lama Sekolah 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
## Peta Pertumbuhan Ekonomi
ggplot(data = indo_merged) +
geom_sf(aes(fill = Pertumbuhan.Ekonomi..X6.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "plasma", direction = -1,
name = "Laju Tumbuh PDRB") +
labs(
title = "Laju Tumbuh PDRB 2022",
caption = "Sumber: BPS (2022)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Bisa dilihat dari hasil, bahwa kualitas lingkungan hidup pada Indonesia memang sudah cukup setara dengan pengecualian pada Papua serta Kalimantan Utara yang masih memiliki lingkungan yang cukup sehat, namun kebalikannya terjadi pada Jawa dimana mereka provinsi memiliki IKLH terendah, terutama DKI Jakarta. Dilihat juga bahwa industri perhutanan-pertanian-perikanan memang terpusat pada Pulau Sumatera dan Jawa yang memiliki industri terbesar. Hal tersebut tetap benar untuk industri pengolahan yang secara utama terpusat pada pulau Jawa. Industri pertambangan bisa dilihat tersebar di seluruh Indonesia dengan perkembangan besar di Papua Barat, Sumatera Utara, Sumatera Selatan, dan Kalimantan Utara (dimana PDRB pertambangan merupakan yang terbesar). Untuk persentase kemiskinan bisa dilihat jelas sebuah perbedaan besar dari provinsi Indonesia lainnya yang memiliki persentase cukup stabil dengan Papua, Papua Barat, dan NTT yang memiliki persentase kemiskinan sangatlah tinggi. Juga bisa dilihat pada Rata Lama Sekolah dimana Papua Barat dan Papua memiliki tingkat yang cukup rendah dibandingkan provinsi lain, kecuali Kalimantan Barat. Untuk laju pertumbuhan PDRB bisa dilihat jelas bahwa semua perekonomian daerah Indonesia berkembang secara stabil, dengan provinsi Maluku Utara dan Sulawesi Tengah mempunyai laju yang cukup tinggi. Kemudian dilakukan visualisasi hubungan faktor-faktor tersebut dengan IKLH melalui scatterplot:
## Scatterplot
Y_var <- "IKLH..Y."
X_vars <- c(
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5."
)
df_long <- indo_merged %>%
st_drop_geometry() %>%
select(all_of(c(Y_var, X_vars))) %>%
pivot_longer(
cols = all_of(X_vars),
names_to = "X_name",
values_to = "X_value"
)
ggplot(df_long, aes(x = X_value, y = !!sym(Y_var))) +
geom_point(color = "steelblue", size = 3) +
geom_smooth(method = "lm", se = TRUE, color = "red") +
facet_wrap(~ X_name, scales = "free", ncol = 2) +
labs(
x = "Nilai X",
y = "IKLH (Y)",
title = "Scatterplot IKLH (Y) terhadap Semua Variabel X",
subtitle = "Dengan garis regresi linear"
) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
Bisa dilihat dari grafik, hubungan Indeks Lingkungan dengan faktor-faktor sosio ekonomik. Dimana dilihat, sering kali perkembangan sektor pengolahan serta perhutanan-pertania-perikanan, dan juga RLS memiliki dampak negatif pada lingkungan. Namun, kemiskinan justru memiliki dampak positif dan sektor pertambangan memiliki dampak positif yang kecil. Kemudian dilakukan analisis deskriptif pada variabel dimana didapatkan:
## Deskriptif Data
data_desc <- indo_merged %>%
st_drop_geometry() %>%
select(
IKLH..Y.,
PDRB.Pertambangan..X2.,
PDRB.Pertambangan..X2.,
PDRB.Industri.Pengolahan..X3.,
Persentase.Kemiskinan..X4.,
RLS..X5.,
Pertumbuhan.Ekonomi..X6.
)
describe(data_desc)
## N mean Std.Dev. min Q1
## 204.000 18207.934 71605.378 1.800 7.805
## median Q3 max missing values
## 40.605 5103.315 683420.250 -198.000
Sebelumnya, diperlakukan penentuan bobot untuk analisis spasial. Pada analisis ini diberlakukan bobot Row-Standardized Weight dimana distribusi merata ke tetangga berdasarkan bobot. Kemudian dibergunakan sistem jarak Queen Contiguinity untuk mengukur jarak ke semua area di unit spasial. Sebelum dilakukan estimasi, dilakukan standarisasi karena perbedaannya variabel numerik antar X dan Y sehingga untuk mempermudah estimasi model. Pertama didapatkan jaringan hubungan tetangga seperti berikut:
nb <- poly2nb(as_Spatial(indo_merged), queen = TRUE)
## Warning in poly2nb(as_Spatial(indo_merged), queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(as_Spatial(indo_merged), queen = TRUE): neighbour object has 11 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
lwW <- nb2listw(nb, style = "W", zero.policy = TRUE)
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)
indo_sp <- as_Spatial(indo_merged)
coords <- coordinates(indo_sp)
plot(st_geometry(indo_merged),
col = "white",
border = "grey40",
main = "Jaringan Tetangga (Queen Contiguity) dengan Titik Centroid")
points(coords, pch = 19, col = "blue", cex = 1.1)
plot(nb, coords,
add = TRUE,
col = "red",
lwd = 2)
Didapatkan hasil pengujian dan pengecekan autokorelasi spasial pada IKLH sebagai berikut:
analyze_spatial <- function(data_sf, var_col, label) {
cat("\n\n==============================\n", label, "\n==============================\n")
x_raw <- data_sf[[var_col]]
# === Uji autokorelasi spasial global (asimptotik) ===
moran_res <- moran.test(x_raw, listw = lwW, zero.policy = TRUE)
geary_res <- geary.test(x_raw, listw = lwW, zero.policy = TRUE)
print(moran_res)
print(geary_res)
# === LISA (Local Moran’s I) ===
x <- scale(x_raw)[, 1]
lagx <- lag.listw(lwW, x, zero.policy = TRUE)
lisa <- localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
x >= 0 & lagx >= 0 ~ "High-High",
x < 0 & lagx < 0 ~ "Low-Low",
x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)
LISA_sf <- bind_cols(data_sf, lisa_df) |>
mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
# === Local Geary’s C ===
localC_vals <- localC(x_raw, lwW, zero.policy = TRUE)
Geary_sf <- mutate(data_sf, LocalC = as.numeric(localC_vals))
# === Getis–Ord Gi* ===
Wb <- listw2mat(lwB)
Wb_star <- Wb; diag(Wb_star) <- 1
G_star_raw <- as.numeric(Wb_star %*% x_raw) / sum(x_raw)
Gz <- localG(x_raw, listw = lwW, zero.policy = TRUE)
G_sf <- mutate(data_sf,
z_Gistar = as.numeric(Gz),
hotcold = case_when(
z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
TRUE ~ "Not significant"
))
# === Visualisasi ===
p1 <- ggplot(LISA_sf) +
geom_sf(aes(fill = quad), color = "white", size = 0.2) +
scale_fill_manual(values = c(
"High-High" = "#d73027",
"Low-Low" = "#4575b4",
"High-Low (Outlier)" = "#fdae61",
"Low-High (Outlier)" = "#74add1",
"Not significant" = "grey85"
)) +
labs(title = paste("Local Moran's I (LISA) —", label),
fill = "Kategori") +
theme_minimal()
p2 <- ggplot(Geary_sf) +
geom_sf(aes(fill = LocalC), color = "white", size = 0.2) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Local Geary’s C") +
labs(title = paste("Local Geary’s C —", label)) +
theme_minimal()
p3 <- ggplot(G_sf) +
geom_sf(aes(fill = hotcold), color = "white", size = 0.2) +
scale_fill_manual(values = c(
"Hot spot (p<0.05)" = "#b2182b",
"Cold spot (p<0.05)" = "#2166ac",
"Not significant" = "grey85"
)) +
labs(title = paste("Getis–Ord Gi* —", label), fill = NULL) +
theme_minimal()
print(p1); print(p2); print(p3)
}
analyze_spatial(indo_merged, "IKLH..Y.", "Indeks Kualitas Lingkungan Hidup")
##
##
## ==============================
## Indeks Kualitas Lingkungan Hidup
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 4.2155, p-value = 1.246e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.69131639 -0.03703704 0.02985230
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 4.636, p-value = 1.776e-06
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.21055896 1.00000000 0.02899745
Berdasarkan dari hasil analisis didapatkan bahwa Indeks Kualitas Lingkungan Hidup memiliki autokorelasi spasial yang positif dan signifikan pada Global Moran’s I dan Geary’s C, yang mengindikasi adanya ketergantungan spasial dimana wilayah berdekatan akan memiliki karakteristik yang sama. Pada analisis lokal Lisa dan Ger-Ord didapatkan cold-spot yaitu pada daerah DKI Jakarta dan provinsi sekitarnya, mengindikasi kualiats lingkungan rendah karena urbanisasi serta hot-spot pada Papua karena kurangnya pembangunan.
Didapatkan hasil pengujian dan pengecekan autokorelasi sebagai berikut:
analyze_spatial(indo_merged, "PDRB.Perhutanan.Pertanian.Perikanan..X1.", "PDRB Sektor Perhutanan/Pertanian/Perikanan")
##
##
## ==============================
## PDRB Sektor Perhutanan/Pertanian/Perikanan
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 2.0745, p-value = 0.01901
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.32179891 -0.03703704 0.02991899
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.3313, p-value = 0.009868
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.60345471 1.00000000 0.02893209
Berdasarkan dari hasil analisis didapatkan bahwa PDRB Sektor Perhutanan-Pertanian-Perikanan memiliki autokorelasi spasial yang positif dan signifikan pada Global Moran’s I dan Geary’s C, yang mengindikasi adanya ketergantungan spasial dimana wilayah berdekatan akan memiliki karakteristik yang sama. Pada analisis lokal Lisa dan Ger-Ord didapatkan hot-spot pada Jawa Tengah, Jawa Timur, dan DI Yogyakarta sebagai pusat Perhutanan-Pertanian-Perikanan Indonesia serta didapatkan hot-spot pada Aceh dan Sumatera Barat.
Didapatkan hasil pengujian dan pengecekan autokorelasi spasial sebagai berikut:
analyze_spatial(indo_merged, "PDRB.Pertambangan..X2.", "PDRB Sektor Pertambangan")
##
##
## ==============================
## PDRB Sektor Pertambangan
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.46348, p-value = 0.6785
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.08962129 -0.03703704 0.01287204
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = -2.0516, p-value = 0.9799
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 1.43828432 1.00000000 0.04563813
Berdasarkan dari hasil analisis didapatkan bahwa PDRB Sektor Pertambangan tidak memiliki autokorelasi spasial yang signifikan pada Global Moran’s maupun Geary C. Serta didapatkan hot-spot pada seluruh wilayah Kalimantan, namun hal tersebut bisa dikatakan efek random karena non-signifikansi autokorelasi spasial.
Didapatkan hasil pengujian dan pengecekan autokorelasi spasial sebagai berikut:
analyze_spatial(indo_merged, "PDRB.Industri.Pengolahan..X3.", "PDRB Sektor Pengolahan")
##
##
## ==============================
## PDRB Sektor Pengolahan
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 3.6037, p-value = 0.0001568
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.50583385 -0.03703704 0.02269291
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.7982, p-value = 0.002569
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.46897555 1.00000000 0.03601366
Berdasarkan dari hasil analisis didapatkan bahwa PDRB Sektor Pengolahan memiliki autokorelasi spasial yang positif dan signifikan pada Global Moran’s I dan Geary’s C, yang mengindikasi adanya ketergantungan spasial dimana wilayah berdekatan akan memiliki karakteristik yang sama. Pada analisis lokal Lisa dan Ger-Ord didapatkan hot-spot pada seluruh pulau Jawa menunjukkan majunya sektor industri pada pulau Jawa.
Didapatkan hasil pengujian dan pengecekan autokorelasi spasial sebagai berikut:
analyze_spatial(indo_merged, "Persentase.Kemiskinan..X4.", "Persentase Kemiskinan")
##
##
## ==============================
## Persentase Kemiskinan
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 3.401, p-value = 0.0003357
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.55317469 -0.03703704 0.03011588
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 4.7943, p-value = 8.161e-07
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.18723685 1.00000000 0.02873914
Berdasarkan dari hasil analisis didapatkan bahwa Persentase Kemiskinan memiliki autokorelasi spasial yang positif dan signifikan pada Global Moran’s I dan Geary’s C, yang mengindikasi adanya ketergantungan spasial dimana wilayah berdekatan akan memiliki karakteristik yang sama. Pada analisis lokal Lisa dan Ger-Ord didapatkan hot-spot pada pulau Papua yang menunjukkan betapa mundurnya pembangunan di Papua.
Didapatkan hasil pengujian dan pengecekan autokorelasi sebagai berikut:
analyze_spatial(indo_merged, "RLS..X5.", "Rata Lama Sekolah")
##
##
## ==============================
## Rata Lama Sekolah
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 0.9199, p-value = 0.1788
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12717596 -0.03703704 0.03186655
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.4023, p-value = 0.008145
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.60508470 1.00000000 0.02702347
Berdasarkan dari hasil analisis didapatkan bahwa Rata Lama Sekolah memiliki autokorelasi spasial yang positif dan signifikan pada Geary’s C namun bukan pada Global Moran’s I, yang mengindikasi bukti adanya ketergantungan spasial yang ada namun bisa dipertanyakan seberapa besar dampaknya. Pada analisis lokal Lisa dan Ger-Ord didapatkan cold-spot pada Papua Barat, yang menunjukkan mundurnya sistem pendidikan pada Papua.
Didapatkan hasil pengujian dan pengecekan autokorelasi sebagai berikut:
analyze_spatial(indo_merged, "Pertumbuhan.Ekonomi..X6.", "Laju Pertumbuhan Ekonomi")
##
##
## ==============================
## Laju Pertumbuhan Ekonomi
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 0.29567, p-value = 0.3837
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.009540027 -0.037037037 0.008648683
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 3.3751, p-value = 0.0003689
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.24698590 1.00000000 0.04977703
Berdasarkan dari hasil analisis didapatkan bahwa Laju Pertumbuhan PDRB tidak memiliki autokorelasi spasial yang signifikan pada Global Moran’s maupun Geary C. Serta tidak didapatkan autokorelasi spasial lokal.
Dilakukan standarisasi data sebagai berikut:
cols_to_scale <- c("IKLH..Y.", "PDRB.Perhutanan.Pertanian.Perikanan..X1.", "PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.", "Persentase.Kemiskinan..X4.", "RLS..X5.", "Pertumbuhan.Ekonomi..X6.")
indo_merged_scale <- indo_merged %>%
mutate(across(all_of(cols_to_scale), ~ scale(.)[,1]))
Dilakukan estimasi model OLS dengan hasil estimasi sebagai berikut:
model_ols <- lm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale)
summary(model_ols)
##
## Call:
## lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5380 -0.5864 0.0267 0.5529 1.5236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.021e-16 1.300e-01 0.000 1.00000
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 1.973e-01 1.908e-01 1.034 0.31045
## PDRB.Pertambangan..X2. 1.506e-01 1.359e-01 1.108 0.27767
## PDRB.Industri.Pengolahan..X3. -5.850e-01 1.845e-01 -3.170 0.00377
## Persentase.Kemiskinan..X4. 3.466e-01 1.556e-01 2.227 0.03444
## RLS..X5. -1.579e-01 1.564e-01 -1.009 0.32179
## Pertumbuhan.Ekonomi..X6. 3.095e-01 1.360e-01 2.275 0.03109
##
## (Intercept)
## PDRB.Perhutanan.Pertanian.Perikanan..X1.
## PDRB.Pertambangan..X2.
## PDRB.Industri.Pengolahan..X3. **
## Persentase.Kemiskinan..X4. *
## RLS..X5.
## Pertumbuhan.Ekonomi..X6. *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7579 on 27 degrees of freedom
## Multiple R-squared: 0.5301, Adjusted R-squared: 0.4256
## F-statistic: 5.076 on 6 and 27 DF, p-value: 0.001332
Dilihat dari hasil estimasi bahwa terdapat tiga variabel sosio-ekonomi yang memiliki pengaruh signifikan pada lingkungan yaitu PDRB sektor pengolahan, Persentase Kemiskinan, dan Pertumbuhan Ekonomi. Kemudian didapatkan R-Squared sebesar 42.56% yang mengindikasi model yang kurang bisa menangkap variasi pada Indeks Lingkungan.
Dilakukan model SAR dengan hasil estimasi sebagai berikut:
slag_model <- lagsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, method = "eigen", zero.policy = TRUE, tol.solve = 1e-10)
summary(slag_model)
##
## Call:lagsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, method = "eigen",
## zero.policy = TRUE, tol.solve = 1e-10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.557413 -0.400612 -0.078526 0.436975 1.361986
##
## Type: lag
## Regions with no neighbours included:
## 2 17 20 21 22 23
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.019015 0.097230 0.1956 0.84495
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.102892 0.143368 0.7177 0.47295
## PDRB.Pertambangan..X2. 0.061083 0.101992 0.5989 0.54924
## PDRB.Industri.Pengolahan..X3. -0.266735 0.145733 -1.8303 0.06721
## Persentase.Kemiskinan..X4. 0.241144 0.121011 1.9927 0.04629
## RLS..X5. -0.087162 0.117923 -0.7391 0.45982
## Pertumbuhan.Ekonomi..X6. 0.258281 0.102320 2.5243 0.01159
##
## Rho: 0.47384, LR test value: 8.7318, p-value: 0.003127
## Asymptotic standard error: 0.12219
## z-value: 3.8778, p-value: 0.0001054
## Wald statistic: 15.037, p-value: 0.0001054
##
## Log likelihood: -30.53278 for lag model
## ML residual variance (sigma squared): 0.32142, (sigma: 0.56694)
## Number of observations: 34
## Number of parameters estimated: 9
## AIC: 79.066, (AIC for lm: 85.797)
## LM test for residual autocorrelation
## test value: 5.8411, p-value: 0.015656
Bisa dilihat bahwa variabel Pertumbuhan Ekonomi dan Persentase Kemiskinan memberikan pengaruh signifikan pada lingkungan. Dengan Asymptotic Standard error sebesar 0.122 dan AIC 79.066.
Dilakukan model SEM dengan hasil estimasi sebagai berikut:
serr_model <- errorsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, method = "eigen", zero.policy = TRUE, tol.solve = 1e-10)
summary(serr_model)
##
## Call:
## errorsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, method = "eigen",
## zero.policy = TRUE, tol.solve = 1e-10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.70128 -0.40242 0.10328 0.38970 0.91223
##
## Type: error
## Regions with no neighbours included:
## 2 17 20 21 22 23
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.021379 0.195697 -0.1092 0.91301
## PDRB.Perhutanan.Pertanian.Perikanan..X1. -0.133515 0.166484 -0.8020 0.42257
## PDRB.Pertambangan..X2. -0.101602 0.082871 -1.2260 0.22019
## PDRB.Industri.Pengolahan..X3. 0.196081 0.174411 1.1242 0.26091
## Persentase.Kemiskinan..X4. 0.322943 0.149386 2.1618 0.03063
## RLS..X5. -0.091549 0.114604 -0.7988 0.42439
## Pertumbuhan.Ekonomi..X6. 0.235826 0.096776 2.4368 0.01482
##
## Lambda: 0.71075, LR test value: 7.424, p-value: 0.0064358
## Asymptotic standard error: 0.086904
## z-value: 8.1786, p-value: 2.2204e-16
## Wald statistic: 66.889, p-value: 3.3307e-16
##
## Log likelihood: -31.18668 for error model
## ML residual variance (sigma squared): 0.28413, (sigma: 0.53304)
## Number of observations: 34
## Number of parameters estimated: 9
## AIC: 80.373, (AIC for lm: 85.797)
Bisa dilihat bahwa variabel Pertumbuhan Ekonomi dan Persentase Kemiskinan memberikan pengaruh signifikan pada lingkungan. Dengan Asymptotic Standard error sebesar 0.087 dan AIC 80.373.
Dilakukan model SDM dengan hasil estimasi sebagai berikut:
sdm_model <- lagsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
summary(sdm_model)
##
## Call:lagsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, Durbin = TRUE, method = "eigen",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.916835 -0.180818 -0.018443 0.255421 0.962953
##
## Type: mixed
## Regions with no neighbours included:
## 2 17 20 21 22 23
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value
## (Intercept) -0.038493 0.075230 -0.5117
## PDRB.Perhutanan.Pertanian.Perikanan..X1. -0.035694 0.125507 -0.2844
## PDRB.Pertambangan..X2. 0.149804 0.084940 1.7637
## PDRB.Industri.Pengolahan..X3. 0.038511 0.138491 0.2781
## Persentase.Kemiskinan..X4. 0.272197 0.112286 2.4241
## RLS..X5. 0.019766 0.090050 0.2195
## Pertumbuhan.Ekonomi..X6. 0.285509 0.071659 3.9843
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.418237 0.146159 2.8615
## lag.PDRB.Pertambangan..X2. 0.322592 0.074451 4.3329
## lag.PDRB.Industri.Pengolahan..X3. -0.721602 0.155636 -4.6365
## lag.Persentase.Kemiskinan..X4. 0.408767 0.203857 2.0052
## lag.RLS..X5. 0.271630 0.183277 1.4821
## lag.Pertumbuhan.Ekonomi..X6. 0.593728 0.211418 2.8083
## Pr(>|z|)
## (Intercept) 0.608878
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.776104
## PDRB.Pertambangan..X2. 0.077791
## PDRB.Industri.Pengolahan..X3. 0.780953
## Persentase.Kemiskinan..X4. 0.015345
## RLS..X5. 0.826264
## Pertumbuhan.Ekonomi..X6. 6.768e-05
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.004216
## lag.PDRB.Pertambangan..X2. 1.471e-05
## lag.PDRB.Industri.Pengolahan..X3. 3.544e-06
## lag.Persentase.Kemiskinan..X4. 0.044946
## lag.RLS..X5. 0.138320
## lag.Pertumbuhan.Ekonomi..X6. 0.004980
##
## Rho: 0.12767, LR test value: 0.58138, p-value: 0.44577
## Asymptotic standard error: 0.16553
## z-value: 0.77129, p-value: 0.44053
## Wald statistic: 0.59489, p-value: 0.44053
##
## Log likelihood: -15.89823 for mixed model
## ML residual variance (sigma squared): 0.14829, (sigma: 0.38508)
## Number of observations: 34
## Number of parameters estimated: 15
## AIC: 61.796, (AIC for lm: 60.378)
## LM test for residual autocorrelation
## test value: 4.8544, p-value: 0.027576
Didapatkan model dengan Asymptotic Standard Error sebesar 0.16553 dengan AIC 61.796. Didapatkan variabel signifikan yaitu Persentase Kemiskinan, Laju Pertumbuhan Ekonomi, kemudian spatial lag PDRB sektor Perhutanan-Pertanian-Perikanan, Lag PDRB Pertambangan, Lag PDRB Pengolahan, Lag Persentase Kemiskinan, dan Lag Laju Pertumbuhan Ekonomi.
Dilakukan model SLX dengan hasil estimasi berikut:
slx_model <- lmSLX(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, Durbin = TRUE, zero.policy = TRUE)
summary(slx_model)
##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.037353 0.096487 -0.387129
## PDRB.Perhutanan.Pertanian.Perikanan..X1. -0.013971 0.157127 -0.088916
## PDRB.Pertambangan..X2. 0.167107 0.106768 1.565138
## PDRB.Industri.Pengolahan..X3. -0.001177 0.163856 -0.007182
## Persentase.Kemiskinan..X4. 0.276035 0.144405 1.911527
## RLS..X5. 0.023897 0.115759 0.206440
## Pertumbuhan.Ekonomi..X6. 0.290000 0.092041 3.150761
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.417086 0.187550 2.223871
## lag.PDRB.Pertambangan..X2. 0.341523 0.095033 3.593744
## lag.PDRB.Industri.Pengolahan..X3. -0.781796 0.173352 -4.509872
## lag.Persentase.Kemiskinan..X4. 0.465802 0.244575 1.904539
## lag.RLS..X5. 0.265661 0.235785 1.126709
## lag.Pertumbuhan.Ekonomi..X6. 0.655255 0.262510 2.496114
## Pr(>|t|)
## (Intercept) 0.702556
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.929991
## PDRB.Pertambangan..X2. 0.132496
## PDRB.Industri.Pengolahan..X3. 0.994338
## Persentase.Kemiskinan..X4. 0.069680
## RLS..X5. 0.838435
## Pertumbuhan.Ekonomi..X6. 0.004823
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.037260
## lag.PDRB.Pertambangan..X2. 0.001709
## lag.PDRB.Industri.Pengolahan..X3. 0.000192
## lag.Persentase.Kemiskinan..X4. 0.070630
## lag.RLS..X5. 0.272576
## lag.Pertumbuhan.Ekonomi..X6. 0.020953
Didapatkan variabel signifikan yaitu Persentase Kemiskinan, Laju Pertumbuhan Ekonomi, kemudian spatial lag PDRB sektor Perhutanan-Pertanian-Perikanan, Lag PDRB Pertambangan, Lag PDRB Pengolahan, dan Lag Laju Pertumbuhan Ekonomi.
Dilakukan model SDEM dengan hasil estimasi berikut:
sdem_model <- errorsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, Durbin = TRUE,
method = "eigen", zero.policy = TRUE)
summary(sdem_model)
##
## Call:
## errorsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, Durbin = TRUE, method = "eigen",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.94104 -0.19160 -0.03756 0.23533 0.94729
##
## Type: error
## Regions with no neighbours included:
## 2 17 20 21 22 23
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value
## (Intercept) -0.0427506 0.0711167 -0.6011
## PDRB.Perhutanan.Pertanian.Perikanan..X1. -0.0148052 0.1259030 -0.1176
## PDRB.Pertambangan..X2. 0.1702693 0.0843283 2.0191
## PDRB.Industri.Pengolahan..X3. 0.0037022 0.1318628 0.0281
## Persentase.Kemiskinan..X4. 0.2743728 0.1145525 2.3952
## RLS..X5. 0.0272468 0.0910443 0.2993
## Pertumbuhan.Ekonomi..X6. 0.2903196 0.0721282 4.0250
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.4276600 0.1508160 2.8356
## lag.PDRB.Pertambangan..X2. 0.3553650 0.0745292 4.7681
## lag.PDRB.Industri.Pengolahan..X3. -0.7839199 0.1375284 -5.7001
## lag.Persentase.Kemiskinan..X4. 0.4738373 0.1891247 2.5054
## lag.RLS..X5. 0.2589643 0.1831920 1.4136
## lag.Pertumbuhan.Ekonomi..X6. 0.6844063 0.2013032 3.3999
## Pr(>|z|)
## (Intercept) 0.5477511
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.9063908
## PDRB.Pertambangan..X2. 0.0434744
## PDRB.Industri.Pengolahan..X3. 0.9776015
## Persentase.Kemiskinan..X4. 0.0166126
## RLS..X5. 0.7647343
## Pertumbuhan.Ekonomi..X6. 5.696e-05
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.0045734
## lag.PDRB.Pertambangan..X2. 1.859e-06
## lag.PDRB.Industri.Pengolahan..X3. 1.198e-08
## lag.Persentase.Kemiskinan..X4. 0.0122305
## lag.RLS..X5. 0.1574729
## lag.Pertumbuhan.Ekonomi..X6. 0.0006742
##
## Lambda: -0.089296, LR test value: 0.15714, p-value: 0.69181
## Asymptotic standard error: 0.19171
## z-value: -0.46578, p-value: 0.64137
## Wald statistic: 0.21695, p-value: 0.64137
##
## Log likelihood: -16.11035 for error model
## ML residual variance (sigma squared): 0.15062, (sigma: 0.38809)
## Number of observations: 34
## Number of parameters estimated: 15
## AIC: 62.221, (AIC for lm: 60.378)
Asynptotic standard error sebesar 0.19171 dengan AIC sebesar 62.221. Didapatkan variabel signifikan yaitu PDRB Sektor Pertambanganm Persentase Kemiskinan, Laju Pertumbuhan Ekonomi, kemudian spatial lag PDRB sektor Perhutanan-Pertanian-Perikanan, Lag PDRB Pertambangan, Lag PDRB Pengolahan, Lag Persentase Kemiskinan, dan Lag Laju Pertumbuhan Ekonomi.
Dilakukan model SAC dengan hasil estimasi berikut:
sac_model <- sacsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, method = "eigen", zero.policy = TRUE)
summary(sac_model)
##
## Call:sacsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, method = "eigen",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.974580 -0.201381 -0.024682 0.254196 0.661781
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.013863 0.038430 -0.3607 0.7183055
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.156925 0.062351 2.5168 0.0118429
## PDRB.Pertambangan..X2. 0.182174 0.057577 3.1640 0.0015561
## PDRB.Industri.Pengolahan..X3. -0.218597 0.082075 -2.6634 0.0077362
## Persentase.Kemiskinan..X4. 0.216895 0.065584 3.3071 0.0009425
## RLS..X5. 0.067865 0.073171 0.9275 0.3536706
## Pertumbuhan.Ekonomi..X6. 0.276165 0.064709 4.2678 1.974e-05
##
## Rho: 0.72824
## Asymptotic standard error: 0.076356
## z-value: 9.5374, p-value: < 2.22e-16
## Lambda: -0.8444
## Asymptotic standard error: 0.090331
## z-value: -9.3479, p-value: < 2.22e-16
##
## LR test value: 21.345, p-value: 2.3169e-05
##
## Log likelihood: -24.22601 for sac model
## ML residual variance (sigma squared): 0.13465, (sigma: 0.36694)
## Number of observations: 34
## Number of parameters estimated: 10
## AIC: 68.452, (AIC for lm: 85.797)
AIC sebesar 68.452. Didapatkan variabel signifikan yaitu PDRB sektor perhutanan-pertanian-perikanan, PDRB Pertambangan, PDRB Pengolahan, Persentase Kemiskinan, dan Laju Pertumbuhan Ekonomi.
Dilakukan model GNS dengan hasil estimasi berikut:
gns_model <- sacsarlm(IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. + Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
data = indo_merged_scale, listw = lwW, Durbin = TRUE,
method = "eigen", zero.policy = TRUE)
summary(gns_model)
##
## Call:sacsarlm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1. +
## PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6.,
## data = indo_merged_scale, listw = lwW, Durbin = TRUE, method = "eigen",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.729428 -0.154775 -0.030504 0.149173 0.618683
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value
## (Intercept) -0.069621 0.037880 -1.8379
## PDRB.Perhutanan.Pertanian.Perikanan..X1. -0.084118 0.129529 -0.6494
## PDRB.Pertambangan..X2. 0.044150 0.083736 0.5273
## PDRB.Industri.Pengolahan..X3. 0.188855 0.138791 1.3607
## Persentase.Kemiskinan..X4. 0.252646 0.093271 2.7087
## RLS..X5. 0.034875 0.079319 0.4397
## Pertumbuhan.Ekonomi..X6. 0.268897 0.055763 4.8222
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.426700 0.163905 2.6033
## lag.PDRB.Pertambangan..X2. 0.362478 0.086341 4.1982
## lag.PDRB.Industri.Pengolahan..X3. -0.516840 0.156204 -3.3087
## lag.Persentase.Kemiskinan..X4. 0.228753 0.159017 1.4385
## lag.RLS..X5. 0.240083 0.155922 1.5398
## lag.Pertumbuhan.Ekonomi..X6. 0.421793 0.166834 2.5282
## Pr(>|z|)
## (Intercept) 0.0660728
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.5160716
## PDRB.Pertambangan..X2. 0.5980178
## PDRB.Industri.Pengolahan..X3. 0.1736053
## Persentase.Kemiskinan..X4. 0.0067542
## RLS..X5. 0.6601665
## Pertumbuhan.Ekonomi..X6. 1.42e-06
## lag.PDRB.Perhutanan.Pertanian.Perikanan..X1. 0.0092319
## lag.PDRB.Pertambangan..X2. 2.69e-05
## lag.PDRB.Industri.Pengolahan..X3. 0.0009371
## lag.Persentase.Kemiskinan..X4. 0.1502811
## lag.RLS..X5. 0.1236179
## lag.Pertumbuhan.Ekonomi..X6. 0.0114643
##
## Rho: 0.57094
## Asymptotic standard error: 0.13382
## z-value: 4.2663, p-value: 1.9871e-05
## Lambda: -0.76713
## Asymptotic standard error: 0.13212
## z-value: -5.8061, p-value: 6.3952e-09
##
## LR test value: 43.977, p-value: 5.7463e-07
##
## Log likelihood: -12.9101 for sacmixed model
## ML residual variance (sigma squared): 0.084704, (sigma: 0.29104)
## Number of observations: 34
## Number of parameters estimated: 16
## AIC: 57.82, (AIC for lm: 85.797)
Terdapat AIC sebesar 57.82. Didapatkan variabel signifikan yaitu Persentase Kemiskinan, Laju Pertumbuhan Ekonomi, kemudian spatial lag PDRB sektor perhutanan-pertanian-perikanan, Lag PDRB Pertambangan, Lag PDRB Pengolahan, dan Lag Laju Pertumbuhan Ekonomi.
Digunakan uji Lagrange Multiplier untuk melihat adanya pengaruh spasial serta adanya omitted variabel pada model yang tersederhana (OLS) sebagai alasan empiris menggunakan model spasial.
Didapatkan hasil uji Lagrange Multiplier OLS sebagai berikut:
lm_tests <- lm.LMtests(model_ols, listw = lwW, test=c("LMerr","LMlag","RLMerr","RLMlag","SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
print(lm_tests)
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1.
## + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6., data
## = indo_merged_scale)
## test weights: listw
##
## RSerr = 0.12297, df = 1, p-value = 0.7258
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1.
## + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6., data
## = indo_merged_scale)
## test weights: listw
##
## RSlag = 6.1276, df = 1, p-value = 0.01331
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1.
## + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6., data
## = indo_merged_scale)
## test weights: listw
##
## adjRSerr = 8.7641, df = 1, p-value = 0.003072
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1.
## + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6., data
## = indo_merged_scale)
## test weights: listw
##
## adjRSlag = 14.769, df = 1, p-value = 0.0001215
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = IKLH..Y. ~ PDRB.Perhutanan.Pertanian.Perikanan..X1.
## + PDRB.Pertambangan..X2. + PDRB.Industri.Pengolahan..X3. +
## Persentase.Kemiskinan..X4. + RLS..X5. + Pertumbuhan.Ekonomi..X6., data
## = indo_merged_scale)
## test weights: listw
##
## SARMA = 14.892, df = 2, p-value = 0.0005839
Didapatkan hasi analisis yang cukup menarik. Pada uji Lagrange Lag serta yang adjusted dibuktikan adanya dependensi spasial pada variabel dependen yaitu Indeks Kualitas Lingkungan Hidup dengan singifikansinya kedua uji. Namun, dibandingkan dengan model error didapatkan hasil yang cukup berbeda dimana tanpa adjust Uji Lagrangen akan keluar non-signifikan walaupun dengan adjust maka akan keluar signifikansi, menunjukkan adanya pengaruh omitted variable. Karena uji SARMA didapatkan singifikan, serta berdasarkan dari hasil riset sebelumnya maka untuk riset ini akan diasumsikam bahwa omitted variabel memang ada dan memengaruhi model secara signifikan.
Kemudian dilakukan Likelihood Ratio test pertama untuk membandingkan model OLS dengan Spasial Ekonometrika kemudian dari antara nested. Didapatkan hasil sebagai berikut:
spatialreg::LR.Sarlm(slag_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 8.7318, df = 1, p-value = 0.003127
## sample estimates:
## Log likelihood of slag_model Log likelihood of model_ols
## -30.53278 -34.89870
spatialreg::LR.Sarlm(serr_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 7.424, df = 1, p-value = 0.006436
## sample estimates:
## Log likelihood of serr_model Log likelihood of model_ols
## -31.18668 -34.89870
spatialreg::LR.Sarlm(sdm_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 38.001, df = 7, p-value = 3.029e-06
## sample estimates:
## Log likelihood of sdm_model Log likelihood of model_ols
## -15.89823 -34.89870
spatialreg::LR.Sarlm(slx_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 37.42, df = 6, p-value = 1.458e-06
## sample estimates:
## Log likelihood of slx_model Log likelihood of model_ols
## -16.18892 -34.89870
spatialreg::LR.Sarlm(sdem_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 37.577, df = 7, p-value = 3.647e-06
## sample estimates:
## Log likelihood of sdem_model Log likelihood of model_ols
## -16.11035 -34.89870
spatialreg::LR.Sarlm(sac_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 21.345, df = 2, p-value = 2.317e-05
## sample estimates:
## Log likelihood of sac_model Log likelihood of model_ols
## -24.22601 -34.89870
spatialreg::LR.Sarlm(gns_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 43.977, df = 8, p-value = 5.746e-07
## sample estimates:
## Log likelihood of gns_model Log likelihood of model_ols
## -12.9101 -34.8987
spatialreg::LR.Sarlm(sdm_model, slag_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 29.269, df = 6, p-value = 5.41e-05
## sample estimates:
## Log likelihood of sdm_model Log likelihood of slag_model
## -15.89823 -30.53278
spatialreg::LR.Sarlm(sdem_model, serr_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 30.153, df = 6, p-value = 3.677e-05
## sample estimates:
## Log likelihood of sdem_model Log likelihood of serr_model
## -16.11035 -31.18668
spatialreg::LR.Sarlm(gns_model, sdm_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 5.9763, df = 1, p-value = 0.0145
## sample estimates:
## Log likelihood of gns_model Log likelihood of sdm_model
## -12.91010 -15.89823
Didapatkan bahwa semua model spasial dapat menejelaskan dampak dependensi spasial pada variabel dependen, independen, error, atau semuanya dengan baik. Serta model OLS dibuktikan tidak seoptimal model spasial menunjukkan adanya interdependensi spasial yang signifikan. Kemudian berdasarkan pengujian antar nested bisa dilihat bahwa model GNS signifikan dari SDM yang signifikan dari SAR model. Hal tersebut menunjukkan adanya dependensi spasial pada variabel independen serta omitted variable ditambah pada dependensi spasial dependen.
JELASKAN BERDASARKAN SUMBER BAHWA GNS YANG TERBAIK
Terakhir adalah pemilihan model, menggunakan AIC. Dimana didapatkan hasil sebagai berikut:
model_comparison <- data.frame(
Model = c("OLS", "SAR", "SEM", "SDM", "SLX", "SDEM", "SAC", "GNS"),
AIC = c(
AIC(model_ols),
AIC(slag_model),
AIC(serr_model),
AIC(sdm_model),
AIC(slx_model),
AIC(sdem_model),
AIC(sac_model),
AIC(gns_model)
),
LogLik = c(
logLik(model_ols),
logLik(slag_model),
logLik(serr_model),
logLik(sdm_model),
logLik(slx_model),
logLik(sdem_model),
logLik(sac_model),
logLik(gns_model)
)
)
model_comparison <- model_comparison[order(model_comparison$AIC), ]
print(model_comparison)
## Model AIC LogLik
## 8 GNS 57.82021 -12.91010
## 5 SLX 60.37785 -16.18892
## 4 SDM 61.79647 -15.89823
## 6 SDEM 62.22071 -16.11035
## 7 SAC 68.45201 -24.22601
## 2 SAR 79.06556 -30.53278
## 3 SEM 80.37336 -31.18668
## 1 OLS 85.79740 -34.89870
Dilihat dari hasil menunjukkan model GNS dengan AIC terkecil. Dan dari hasil ditunjukkan bahwa model GNS merupakan yang teroptimal dibandingkan model yang lain.
Untuk pemilihan model, bisa dilihat dengan jelas berdasarkan bukti teoritis serta empiris. Dibuktikan secara teoritis bahwa kualitas lingkungan memang memiliki dependensi spasial karena lingkungan hidup tidak dibatasi oleh perbatasan administrasi sehingga kualitas lingkungan di sebuah provinsi pasti akan memberi dampak terhadap provinsi lain karena adanya polusi atau faktor lain. Secara teoritis juga dibuktikan bahwa PDRB memiliki dependensi spasial karena perkembangan ekonomi serta akses kapital juga tidak dibatasi batas administrasi terutama didalam sebuah negara serta hal seperti kemiskinan dan pertumbuhan ekonomi karena pembangunannya wilayah di daerah tersebut tidak sesuai batas administrasi. Secara empiris juga terbukti bahwa dependensi spasial pada variabel independen, dependen, dan omitted variabel ada dan signifikan. Dan juga berdasarkan AIC model bisa disimpulkan bahwa GNS yang melakukan estimasi dependensi spasial indpenden, memperhitungkan dependensi spasial pada variabel dependen, dan menangani omitted variabel.
Dengan didapatkannya model terbaik, dilakukan pengujian White Noise serta permasalahan data dalam bentuk multikolinearitas. Untuk pengujian white noise akan digunakan uji Moran’s I untuk menguji autokorelasi, Shapiro-Wilks untuk menguji normalitas, dan Breush Pagan untuk menguji homogenitas dimana didapatkan hasil sebagai berikut:
moran_rsd <- moran.test(as.numeric(resid(gns_model)), listw = lwW, zero.policy = T)
sprwlk <- shapiro.test(as.numeric(resid(gns_model)))
homo <- bptest.Sarlm(gns_model)
moran_rsd;sprwlk;homo
##
## Moran I test under randomisation
##
## data: as.numeric(resid(gns_model))
## weights: lwW
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.48235, p-value = 0.6852
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.12256773 -0.03703704 0.03144331
##
## Shapiro-Wilk normality test
##
## data: as.numeric(resid(gns_model))
## W = 0.96895, p-value = 0.433
##
## studentized Breusch-Pagan test
##
## data:
## BP = 18.265, df = 12, p-value = 0.1079
Didapatkan hasil pada ketiga uji white noise, pada Moran’s I karena (P-Value > 0.05) bisa disimpulkan bahwa tidak adanya autokorelasi spasial yang terjadi pada residual model. Karena Shapiro-Wilks memiliki (P-Value > 0.05) maka disimpulkan bahwa residual diambil dari populasi yang berdistribusi normal, dan Breusch-Pagan memiliki (P-Value > 0.05) maka bisa disimpulkan bahwa varians pada residual bersifat konstan. Bisa disimpulkan bahwa model telah memenuhi asumsi white noise dimana model sudah optimal.
Didapatkan hasil uji mutlikolinearitas menggunakan VIF (Variance Inflation Factor) sebagai berikut:
multikol <- vif(model_ols)
multikol
## PDRB.Perhutanan.Pertanian.Perikanan..X1.
## 2.092361
## PDRB.Pertambangan..X2.
## 1.061303
## PDRB.Industri.Pengolahan..X3.
## 1.956698
## Persentase.Kemiskinan..X4.
## 1.390896
## RLS..X5.
## 1.405674
## Pertumbuhan.Ekonomi..X6.
## 1.063453
Didapatkan nilai VIF pada semua variabel prediktor sebagai diatas dimana tidak ada nilai yang melebihi dari 5, menunjukkan tidak adanya multikolinearitas pada data sehingga variabel prediktor tidak saling berkorelasi satu sama lain.
Dilakukan analisis interpolasi, dimana data-data yang digunakan diperlakukan sebagai data titik, dimana titik tersebut akan menjadi sentorid atau titik tengah di area provinsi masing-masing data. Digunakan 5 metode yang kemudian dibandingkan berdasarkan RMSE serta diujikan pada asumsi yang dibutuhkan.
Pertama untuk melihat karakteristik dari pengaruh spasial variabel via jarak akan dilakukan analisis menggunakan variogram. Didapatkan hasil sebagai berikut:
variogram_plot <- function(indo_merged, var) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
sp_pts <- as(indo_centroid, "Spatial")
sp_pts@data$z <- as.numeric(sp_pts@data[[var]])
stopifnot(is.numeric(sp_pts@data$z))
vgm_emp <- variogram(z ~ 1, sp_pts)
vgm_fit <- fit.variogram(
vgm_emp,
model = vgm(psill = var(sp_pts@data$z),
model = "Sph",
range = max(vgm_emp$dist) / 2,
nugget = 0)
)
list(
variogram_empirical = vgm_emp,
variogram_model = vgm_fit
)
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6."
)
for (v in vars_interp) {
res_obj <- variogram_plot(indo_merged, v)
p <- plot(
res_obj$variogram_empirical,
res_obj$variogram_model,
main = v
)
print(p)
}
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, model = vgm(psill = var(sp_pts@data$z), : No
## convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: No convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, model = vgm(psill = var(sp_pts@data$z), :
## singular model in variogram fit
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, model = vgm(psill = var(sp_pts@data$z), : No
## convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: No convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
Dilihat secara jelas adanya trend global pada IKLH dan Rata Lama Sekolah dimana semakin jauh jarak maka semakin naik semivarians. Hal juga bisa dilihat pada PDRB Pertambangan, Persentase Kemiskinan, dan Pertumbuhan Ekonomi yang menunjukkan autokorelasi spasial yang kuat, menunjukkan tren lokal yang cukup kuat. Walaupun PDRB Perhutanan-Pertanian-Perikanan dan Industri Pengolahan memiliki indikasi autokroelasi spasial yang lemah.
Kemudian dilakukan analisis interpolasi spasial menggunakan nearest neighbour didapatkan hasil sebagai berikut:
nn_voronoi <- function(sf_poly, var, crs_proj = 3857) {
pts <- st_centroid(sf_poly)
pts <- st_transform(pts, crs_proj)
pts[[var]] <- as.numeric(pts[[var]])
boundary <- st_union(st_transform(sf_poly, crs_proj))
vor <- st_voronoi(
st_combine(st_geometry(pts)),
envelope = st_as_sfc(st_bbox(boundary))
)
vor_sf <- st_collection_extract(vor, "POLYGON") |>
st_sf() |>
st_intersection(boundary)
idx <- st_nearest_feature(vor_sf, pts)
vor_sf$value <- pts[[var]][idx]
p <- ggplot() +
geom_sf(data = vor_sf, aes(fill = value), color = "white", linewidth = 0.2) +
geom_sf(data = pts, color = "black", size = 1) +
scale_fill_viridis_c() +
labs(
title = paste("Nearest Neighbor (Voronoi) -", var),
fill = var
) +
theme_minimal()
list(
voronoi = vor_sf,
points = pts,
plot = p
)
}
vars <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6."
)
res_nn <- lapply(vars, \(v) nn_voronoi(indo_merged, v))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
res_nn[[1]]$plot
res_nn[[2]]$plot
res_nn[[3]]$plot
res_nn[[4]]$plot
res_nn[[5]]$plot
res_nn[[6]]$plot
res_nn[[7]]$plot
Bisa dilihat pendugaan variabel yang tidak bsia dilihat berdasarkan titik serta pembagian voronoi yang bisa mengikuti pola saat data berjenis area/polygon dengan akurasi yang cukup bagus.
Kemudian dilakukan analisis interpolasi spasial menggunakan metode inverted distance weight dengan hasil seperti berikut:
idw_plot <- function(indo_merged, var, power = 2) {
indo_centroid <- indo_merged |>
st_centroid() |>
dplyr::select(all_of(var)) |>
st_transform(4326)
sp_pts <- as(indo_centroid, "Spatial")
r <- raster(
extent(st_transform(indo_merged, 4326)),
res = 0.25,
crs = CRS("+proj=longlat +datum=WGS84")
)
idw_r <- raster::interpolate(
r,
gstat::gstat(
formula = as.formula(paste(var, "~ 1")),
locations = sp_pts,
set = list(idp = power)
)
)
indo_sp <- as(st_transform(indo_merged, 4326), "Spatial")
idw_r <- raster::mask(idw_r, indo_sp)
idw_df <- as.data.frame(idw_r, xy = TRUE, na.rm = TRUE)
colnames(idw_df) <- c("lon", "lat", "value")
p <- ggplot() +
geom_tile(
data = idw_df,
aes(lon, lat, fill = value)
) +
geom_sf(
data = st_transform(indo_merged, 4326),
fill = NA, color = "black", linewidth = 0.4
) +
scale_fill_viridis_c(name = var) +
coord_sf(expand = FALSE) +
labs(title = paste("IDW Interpolation -", var)) +
theme_minimal()
list(
raster = idw_r,
plot_idw = p
)
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6."
)
for(v in vars_interp){
var_name <- gsub("\\W+", "_", v)
res_obj <- idw_plot(indo_merged, v)
assign(paste0("res_", var_name, "_idw"), res_obj)
print(res_obj$plot_idw)
}
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
## Warning: st_centroid assumes attributes are constant over geometries
## [inverse distance weighted interpolation]
Dilihat penyebaran kontinu untuk semua variabel yang digunakan, semua mengikuti pola yang cukup dekat dengan data tersebut saat berjenis poligon/area.
Kemudian dilakukan analisis interpolasi spasial menggunakan thin-plate splines dengan hasil seperti berikut:
tps_plot <- function(indo_merged, var) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
coords <- st_coordinates(indo_centroid)[, 1:2]
z <- as.numeric(indo_centroid[[var]])
stopifnot(is.numeric(z))
stopifnot(is.matrix(coords))
stopifnot(ncol(coords) == 2)
tps_model <- fields::Tps(coords, z)
r <- raster(
extent(st_transform(indo_merged, 4326)),
res = 0.25,
crs = "+proj=longlat +datum=WGS84"
)
xy <- as.data.frame(coordinates(r))
colnames(xy) <- c("lon", "lat")
xy$value <- predict(tps_model, as.matrix(xy))
r[] <- xy$value
indo_sp <- as(st_transform(indo_merged, 4326), "Spatial")
r <- raster::mask(r, indo_sp)
df <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
colnames(df) <- c("lon", "lat", "value")
p <- ggplot() +
geom_tile(data = df, aes(lon, lat, fill = value)) +
geom_sf(
data = st_transform(indo_merged, 4326),
fill = NA, color = "black", linewidth = 0.4
) +
scale_fill_viridis_c(name = var) +
coord_sf(expand = FALSE) +
labs(title = paste("TPS Interpolation -", var)) +
theme_minimal()
list(
raster = r,
plot_tps = p
)
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
for(v in vars_interp){
var_name <- gsub("\\W+", "_", v)
res_obj <- tps_plot(indo_merged, v)
assign(paste0("res_tps_", var_name), res_obj)
print(res_obj$plot_tps)
}
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.178074e-05 (eff. df= 32.29999 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.178074e-05 (eff. df= 32.29999 )
Dibandingkan dengan IDW, pendugaan variabel yang tidak terobservasi memang lebih meluas dan tidak memusat pada satu area.
Kemudian dilakukan Ordinary Kriging untuk interpolasi dimana didapatkan hasil sebagai berikut:
ok_plot <- function(indo_merged, var, res = 0.25) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
sp_pts <- as(indo_centroid, "Spatial")
sp_pts@data$z <- as.numeric(as.character(sp_pts@data[[var]]))
vgm_fit <- fit.variogram(
variogram(z ~ 1, sp_pts),
vgm("Sph")
)
r <- raster(
extent(st_transform(indo_merged, 4326)),
res = res,
crs = "+proj=longlat +datum=WGS84"
)
grid <- as(r, "SpatialPixels")
ok <- krige(
z ~ 1,
sp_pts,
grid,
model = vgm_fit
)
r[] <- ok$var1.pred
indo_sp <- as(st_transform(indo_merged, 4326), "Spatial")
r <- mask(r, indo_sp)
df <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
colnames(df) <- c("lon", "lat", "value")
p <- ggplot() +
geom_tile(data = df, aes(lon, lat, fill = value)) +
geom_sf(data = st_transform(indo_merged, 4326),
fill = NA, color = "black", linewidth = 0.4) +
scale_fill_viridis_c(name = var) +
coord_sf(expand = FALSE) +
labs(title = paste("Ordinary Kriging -", var)) +
theme_minimal()
list(
raster = r,
plot_ok = p,
variogram = vgm_fit
)
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Pertambangan..X2.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6."
)
for (v in vars_interp) {
var_name <- gsub("\\W+", "_", v)
res_obj <- ok_plot(indo_merged, v)
print(res_obj$plot_ok)
}
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ 1, sp_pts), vgm("Sph")): No convergence
## after 200 iterations: try different initial values?
## [using ordinary kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: No convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## [using ordinary kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: No convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## [using ordinary kriging]
Karena lemahnya autokorelasi spasial pada PDRB Perhutanan-Pertanian-Perikanan serta PDRB Industri Pengolahan, melakukan analisis interpolasi pada kedua variabel tidak memungkinkan karena semivarians yang kecil.
Dan sebagai perbandingan akan dilakukan interpolasi spasial dengan Universal Kriging didapatkan hasil sebagai berikut:
uk_plot <- function(indo_merged, var, res = 0.25) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
coords <- st_coordinates(indo_centroid)
sp_pts <- as(indo_centroid, "Spatial")
sp_pts@data$z <- as.numeric(sp_pts@data[[var]])
sp_pts@data$lon <- coords[,1]
sp_pts@data$lat <- coords[,2]
vgm_emp <- variogram(z ~ lon + lat, sp_pts)
vgm_fit <- fit.variogram(
vgm_emp,
vgm(psill = var(sp_pts@data$z),
model = "Sph",
range = max(vgm_emp$dist)/2,
nugget = 0)
)
bbox <- st_bbox(st_transform(indo_merged, 4326))
grd <- expand.grid(
lon = seq(bbox["xmin"], bbox["xmax"], by = res),
lat = seq(bbox["ymin"], bbox["ymax"], by = res)
)
coordinates(grd) <- ~lon+lat
gridded(grd) <- TRUE
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84")
uk <- krige(
z ~ lon + lat,
sp_pts,
grd,
model = vgm_fit
)
r <- raster(uk["var1.pred"])
indo_sp <- as(st_transform(indo_merged, 4326), "Spatial")
r <- mask(r, indo_sp)
df <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
colnames(df) <- c("lon", "lat", "value")
p <- ggplot() +
geom_tile(data = df, aes(lon, lat, fill = value)) +
geom_sf(
data = st_transform(indo_merged, 4326),
fill = NA, color = "black", linewidth = 0.4
) +
scale_fill_viridis_c(name = var) +
coord_sf(expand = FALSE) +
labs(title = paste("Global Kriging (Universal) -", var)) +
theme_minimal()
list(
raster = r,
plot_uk = p,
variogram = vgm_fit
)
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Pertambangan..X2.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6."
)
for (v in vars_interp) {
var_name <- gsub("\\W+", "_", v)
res_obj <- uk_plot(indo_merged, v)
print(res_obj$plot_uk)
}
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, vgm(psill = var(sp_pts@data$z), model =
## "Sph", : linear model has singular covariance matrix
## Warning in fit.variogram(vgm_emp, vgm(psill = var(sp_pts@data$z), model =
## "Sph", : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## [using universal kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, vgm(psill = var(sp_pts@data$z), model =
## "Sph", : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial
## values?
## [using universal kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## [using universal kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(vgm_emp, vgm(psill = var(sp_pts@data$z), model =
## "Sph", : No convergence after 200 iterations: try different initial values?
## [using universal kriging]
## Warning: st_centroid assumes attributes are constant over geometries
## [using universal kriging]
Dilakukan perbandingan via RMSE untuk masing-masing metode pada masing-masing variabel serta konsiderasi lainnya seperti adanya autokorelasi spasial melalui variogram.
Didapatkan RMSE untuk Nearest Neighbour sebagai berikut:
rmse_nn <- function(indo_merged, var) {
coords <- st_coordinates(st_centroid(indo_merged))
z <- indo_merged[[var]]
valid <- !is.na(z)
coords <- coords[valid, , drop = FALSE]
z <- z[valid]
n <- length(z)
zhat <- numeric(n)
for (i in seq_len(n)) {
coords_train <- coords[-i, , drop = FALSE]
z_train <- z[-i]
nn_i <- get.knnx(
data = coords_train,
query = matrix(coords[i, ], nrow = 1),
k = 1
)$nn.index[1]
zhat[i] <- z_train[nn_i]
}
sqrt(mean((z - zhat)^2))
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
rmse_results <- data.frame(
variable = vars_interp,
RMSE = sapply(vars_interp, function(v) rmse_nn(indo_merged, v)))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
rmse_results
## variable
## IKLH..Y. IKLH..Y.
## PDRB.Perhutanan.Pertanian.Perikanan..X1. PDRB.Perhutanan.Pertanian.Perikanan..X1.
## PDRB.Pertambangan..X2. PDRB.Pertambangan..X2.
## PDRB.Industri.Pengolahan..X3. PDRB.Industri.Pengolahan..X3.
## Persentase.Kemiskinan..X4. Persentase.Kemiskinan..X4.
## RLS..X5. RLS..X5.
## Pertumbuhan.Ekonomi..X6. Pertumbuhan.Ekonomi..X6.
## RMSE
## IKLH..Y. 4.626195e+00
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 6.771807e+04
## PDRB.Pertambangan..X2. 6.435562e+04
## PDRB.Industri.Pengolahan..X3. 1.570269e+05
## Persentase.Kemiskinan..X4. 4.379690e+00
## RLS..X5. 1.261676e+00
## Pertumbuhan.Ekonomi..X6. 3.347488e+00
Kemudian RMSE untuk Inverse Distance Weight sebagai berikut:
idw_rmse_manual <- function(indo_merged, var, power = 2) {
coords <- st_coordinates(st_centroid(indo_merged))
z <- indo_merged[[var]]
n <- length(z)
zhat <- numeric(n)
for(i in seq_len(n)) {
d <- sqrt((coords[-i,1] - coords[i,1])^2 + (coords[-i,2] - coords[i,2])^2)
w <- 1 / (d^power)
zhat[i] <- sum(z[-i] * w) / sum(w)
}
sqrt(mean((z - zhat)^2, na.rm = TRUE))
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
rmse_results <- data.frame(
variable = vars_interp,
RMSE = sapply(vars_interp, function(v) idw_rmse_manual(indo_merged, v)))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
rmse_results
## variable
## IKLH..Y. IKLH..Y.
## PDRB.Perhutanan.Pertanian.Perikanan..X1. PDRB.Perhutanan.Pertanian.Perikanan..X1.
## PDRB.Pertambangan..X2. PDRB.Pertambangan..X2.
## PDRB.Industri.Pengolahan..X3. PDRB.Industri.Pengolahan..X3.
## Persentase.Kemiskinan..X4. Persentase.Kemiskinan..X4.
## RLS..X5. RLS..X5.
## Pertumbuhan.Ekonomi..X6. Pertumbuhan.Ekonomi..X6.
## RMSE
## IKLH..Y. 3.528029e+00
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 5.083597e+04
## PDRB.Pertambangan..X2. 4.741011e+04
## PDRB.Industri.Pengolahan..X3. 1.476636e+05
## Persentase.Kemiskinan..X4. 4.427066e+00
## RLS..X5. 1.033312e+00
## Pertumbuhan.Ekonomi..X6. 3.095713e+00
Kemudian dilakukan perhitungan RMSE untuk Thin-Spliced sebagai berikut:
tps_rmse <- function(sf_data, var) {
coords <- st_coordinates(st_centroid(sf_data))
z <- sf_data[[var]]
tps_model <- Tps(coords, z)
n <- length(z)
zhat <- numeric(n)
for(i in seq_len(n)) {
coords_train <- coords[-i, , drop = FALSE]
z_train <- z[-i]
tps_loocv <- Tps(coords_train, z_train)
zhat[i] <- predict(tps_loocv, matrix(coords[i, ], nrow = 1))
}
sqrt(mean((z - zhat)^2, na.rm = TRUE))
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Perhutanan.Pertanian.Perikanan..X1.",
"PDRB.Pertambangan..X2.",
"PDRB.Industri.Pengolahan..X3.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
rmse_results_tps <- data.frame(
variable = vars_interp,
RMSE = sapply(vars_interp, function(v) tps_rmse(indo_merged, v)))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001008 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001007 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000707 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000989 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001406 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.00071 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 130.8028 (eff. df= 3.001014 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001407 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 128.695 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000708 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000994 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 130.8028 (eff. df= 3.001018 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.285478e-05 (eff. df= 31.34997 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001008 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.315e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000707 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000989 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.243873e-05 (eff. df= 31.35002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001406 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.00071 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 130.8028 (eff. df= 3.001014 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001407 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 128.695 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000708 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000994 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001008 )
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning in golden.section.search(ax = starts[1], bx = starts[2], cx =
## starts[3], : Maximum iterations reached
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.178074e-05 (eff. df= 32.29999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.483004e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.329713e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.218869e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.315e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.172231e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.19412e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.184811e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.170543e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.13429e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.163609e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.202397e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.144812e-05 (eff. df= 31.34903 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.279828e-05 (eff. df= 31.34994 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.269459e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.488367e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.193969e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.213794e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.197325e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.219092e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.219937e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.148229e-05 (eff. df= 31.35001 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 132.2272 (eff. df= 3.000998 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 130.8028 (eff. df= 3.001018 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001004 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001008 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 129.3938 (eff. df= 3.001007 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000707 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000712 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000989 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.00071 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000714 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 130.8028 (eff. df= 3.001014 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001407 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 128.695 (eff. df= 3.001005 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000708 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000713 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000709 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 90.50967 (eff. df= 3.001411 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 181.0193 (eff. df= 3.000711 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 126.6212 (eff. df= 3.000994 )
## Warning: st_centroid assumes attributes are constant over geometries
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.178074e-05 (eff. df= 32.29999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.483004e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.329713e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.517872e-05 (eff. df= 31.35002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.285478e-05 (eff. df= 31.34997 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.275877e-05 (eff. df= 31.35012 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.478006e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.229125e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.218869e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.377695e-05 (eff. df= 31.35003 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.372269e-05 (eff. df= 31.35003 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.315e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.172231e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.19412e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.184811e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.170543e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.13429e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.243873e-05 (eff. df= 31.35002 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.163609e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.202397e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.144812e-05 (eff. df= 31.34903 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.279828e-05 (eff. df= 31.34994 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.269459e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.488367e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.132586e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.193969e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.213794e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.299845e-05 (eff. df= 31.35001 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.197325e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.275179e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.219092e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.219937e-05 (eff. df= 31.34999 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.320393e-05 (eff. df= 31.35 )
## Warning:
## Grid searches over lambda (nugget and sill variances) with minima at the endpoints:
## (GCV) Generalized Cross-Validation
## minimum at right endpoint lambda = 2.148229e-05 (eff. df= 31.35001 )
rmse_results_tps
## variable
## IKLH..Y. IKLH..Y.
## PDRB.Perhutanan.Pertanian.Perikanan..X1. PDRB.Perhutanan.Pertanian.Perikanan..X1.
## PDRB.Pertambangan..X2. PDRB.Pertambangan..X2.
## PDRB.Industri.Pengolahan..X3. PDRB.Industri.Pengolahan..X3.
## Persentase.Kemiskinan..X4. Persentase.Kemiskinan..X4.
## RLS..X5. RLS..X5.
## Pertumbuhan.Ekonomi..X6. Pertumbuhan.Ekonomi..X6.
## RMSE
## IKLH..Y. 3.332601e+00
## PDRB.Perhutanan.Pertanian.Perikanan..X1. 4.565695e+04
## PDRB.Pertambangan..X2. 4.879168e+04
## PDRB.Industri.Pengolahan..X3. 1.525437e+05
## Persentase.Kemiskinan..X4. 4.090673e+00
## RLS..X5. 9.393086e-01
## Pertumbuhan.Ekonomi..X6. 3.854590e+00
Kemudian dilakukan analisis RMSE untuk Ordinary Kriging dengan hasil sebagai berikut:
rmse_ok <- function(indo_merged, var) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
sp_pts <- as(indo_centroid, "Spatial")
sp_pts@data$z <- as.numeric(sp_pts@data[[var]])
vgm_fit <- fit.variogram(
variogram(z ~ 1, sp_pts),
vgm("Sph")
)
cv <- krige.cv(
z ~ 1,
sp_pts,
model = vgm_fit,
nfold = nrow(sp_pts)
)
sqrt(mean(cv$residual^2, na.rm = TRUE))
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Pertambangan..X2.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
rmse_results_ok <- data.frame(
variable = vars_interp,
RMSE = sapply(vars_interp, function(v) rmse_ok(indo_merged, v)))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ 1, sp_pts), vgm("Sph")): No convergence
## after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ 1, sp_pts), vgm("Sph")): No convergence
## after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ 1, sp_pts), vgm("Sph")): No convergence
## after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
rmse_results_ok
## variable RMSE
## IKLH..Y. IKLH..Y. 3.345214e+00
## PDRB.Pertambangan..X2. PDRB.Pertambangan..X2. 4.666437e+04
## Persentase.Kemiskinan..X4. Persentase.Kemiskinan..X4. 4.514019e+00
## RLS..X5. RLS..X5. 9.579574e-01
## Pertumbuhan.Ekonomi..X6. Pertumbuhan.Ekonomi..X6. 3.258747e+00
Didapatkan
Kemudian dilakukan analisis RMSE menggunakan Universal Kriging dimana didapatkan:
rmse_uk <- function(indo_merged, var) {
indo_centroid <- st_centroid(st_transform(indo_merged, 4326))
coords <- st_coordinates(indo_centroid)
sp_pts <- as(indo_centroid, "Spatial")
sp_pts@data$z <- as.numeric(sp_pts@data[[var]])
sp_pts@data$lon <- coords[,1]
sp_pts@data$lat <- coords[,2]
vgm_fit <- fit.variogram(
variogram(z ~ lon + lat, sp_pts),
vgm("Sph")
)
cv <- krige.cv(
z ~ lon + lat,
sp_pts,
model = vgm_fit,
nfold = nrow(sp_pts)
)
sqrt(mean(cv$residual^2, na.rm = TRUE))
}
vars_interp <- c(
"IKLH..Y.",
"PDRB.Pertambangan..X2.",
"Persentase.Kemiskinan..X4.",
"RLS..X5.",
"Pertumbuhan.Ekonomi..X6.")
rmse_results_uk <- data.frame(
variable = vars_interp,
RMSE = sapply(vars_interp, function(v) rmse_uk(indo_merged, v)))
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ lon + lat, sp_pts), vgm("Sph")): No
## convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ lon + lat, sp_pts), vgm("Sph")): No
## convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial
## values?
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning in fit.variogram(variogram(z ~ lon + lat, sp_pts), vgm("Sph")): No
## convergence after 200 iterations: try different initial values?
## Warning: st_centroid assumes attributes are constant over geometries
rmse_results_uk
## variable RMSE
## IKLH..Y. IKLH..Y. 3.590004e+00
## PDRB.Pertambangan..X2. PDRB.Pertambangan..X2. 4.860643e+04
## Persentase.Kemiskinan..X4. Persentase.Kemiskinan..X4. 4.369686e+00
## RLS..X5. RLS..X5. 9.513136e-01
## Pertumbuhan.Ekonomi..X6. Pertumbuhan.Ekonomi..X6. 3.333689e+00
Untuk IKLH yang didapatkan tren global serta pembuktian empiris dan teoritis adanya autokorelasi spasial didapatkan RMSE terkecil pada metode IDW sebesar 0.575. Namun karena kurangnya konsiderasi pada autokorelasi maka masih bisa dipertanyakan. Namun bisa memenuhi asumsi BLUE pada Kriging universal maka metode tersebut yang terpilih teroptimal, walaupun dengan RMSE terkecil ketiga dibawah TPS.
Untuk PDRB Perhutanan karena kurangnya pengaruh autokorelasi pada variogram, walaupun diuji signifikan, maka tidak bisa digunakan. Sehingga untuk interpolasi tersebut akan sangatlah susah dengan RMSE terkecil pada TPS sebesar 0.98 menunjukkan metode tersebut yang lebih optimal.
Untuk PDRB Pertambangan, secare empiris dibuktikan tidak ada autokorelasi spasial, sehingga seperti PDRB Perhutanan tidak bisa digunakan dalam Kriging walaupun variogram mengindikasi autokorelasi kuat dengan terendah pada IDW sebesar 1.08 menunjukkan metode tersebut yang lebih optimal.
Untuk PDRB Industri karena kurangnya pengaruh autokorelasi pada variogram, walaupun diuji signifikan, maka tidak bisa digunakan. Sehingga untuk interpolasi tersebut akan sangatlah susah dengan RMSE terkecil pada IDW sebesar 0.95 menunjukkan metode tersebut yang lebih optimal.
Untuk RLS atau Rata Lama Sekolah didapatkan RMSE terkecil pada TPS sebesar 1.01 menunjukkan metode tersebut adalah yang teroptimal.
Untuk Laju Pertumbuhan Ekonomi didapatkan RMSE terkecil pada TPS sebesar 1.02 menunjukkan metode tersebut yang lebih optimal, terutama dengan adanya autokorelasi signifikan pada data berdasarkan variogram.
Didapatkan bahwa model spasial ekonometrika General Nested Spatial Model dari konsiderasi empiris maupun teoritis dimana didapatkan model sebagai berikut dengan variabel yang telah distandarisasi:
\[ \begin{aligned} \text{IKLH} =\; -0.069621 - 0.084118 PDRBHutanTaniIkan_i + 0.044150 PDRBTambang_i + 0.188855 PDRBOlah_i + 0.252646 Miskin_i + 0.034875 RLS_i\\ + 0.268897 LajuPDRB_i + 0.4267 \sum_j \omega_{ij} PDRBHutanTaniIkan_j + 0.362478 \sum_j \omega_{ij} PDRBTambang_j\\ - 0.516840 \sum_j \omega_{ij} PDRBOlah_j + 0.228753 \sum_j \omega_{ij} Miskin_j + 0.240083 \sum_j \omega_{ij} RLS_j + 0.421793 \sum_j \omega_{ij} LajuPDRB_j \end{aligned} \] Bisa dilihat dari model beberapa kesimpulan mengenai determinasi faktor ekonomi pada Indeks Kualitas Lingkungan Hidup.
Bisa dilihat faktor PDRB sebagian besar memiliki pengaruh positif pada Lingkungan yaitu sektor Pertambangan dan Pengolahan serta PDRB sektor Perhutanan-Pertanian-Perikanan yang memiliki pengaruh negatif pada Lingkungan. Namun faktor aktivitas ekonomi dari ketiga sektor tersebut tidak memberikan pengaruh singifikan secara langsung pada lingkungan. Namun ketiga faktor memberikan pengaruh lag spasial signifikan pada tetangganya dimana pengaurh lag sektor pertanian serta pertambangan memberikan pengaruh positif secara lag pada lingkungan dan pertambangan memberikan pengaruh negatif.
Kemudian faktor ekonomi-sosial, dilihat bahwa faktor kemiskinan, RLS, dan pertumbuhan ekonomi ketiganya memiliki pengaruh positif pada lingkungan serta juga pengaruh spasial lag. Pengaruh langsung Persentase Kemiskinan dan Pertumbuhan Ekonomu bersifat signifikan pada lingkungan dimana pertumbuhan ekonomi juga memiliki pengaruh lag spasial yang signifikan.
Dilihat dari interpolasi, didapatkan penyebaran kualitas lingkungan juga dimana kualitas lingkungan menurun semakin dekat pada Pulau Jawa sebagai pusat perkembangan ekonomi Indonesia. Hal tersebut juga benar pada PDRB Perhutanan serta PDRB Industri pengolahan, dengan sedikit pembangunan di Sumatera namun tetap terpusat pada pulau Jawa.
Pertambangan juga memiliki pemusatan pada Kalimantan Timur sebagai pusat pertambangan di indonesia, sehingga sangatlah memusat di lokasi tersebut. Namun untuk RLS dan Kemiskinan bisa dilihat kebalikannya, dimana ada pemusatan kemunduran pendidikan serta standar hidup di Papua karena mundurnya pembangunan di lokasi.
Bisa dilihat dari model bahwa faktor ekonomi seperti pembangunan sektor-sektor serta faktor sosio-ekonomi seperti RLS memberikan pengaruh signifikan, secara langsung ataupun melalui spatial lag pada kualitas lingkungan. Namun ada beberapa pertanyaan yang muncul, terutama mengenai hubungan kausal tersebut seperti sektor pertambangan yang memberikan positif pada lingkungan walaupun secara dasar kegiatan tersebut merusak alam.
Karena itu dibutuhkan sebuah pendalaman teori mengenai hubungan kausal faktor ekonomi pada lingkungan serta pemilihan model yang lebih bisa menjelaskan hubungan kasual tersebut.
Listiyani, N. (2017). Dampak Pertambangan Terhadap Lingkungan Hidup Di Kalimantan Selatan Dan Implikasinya Bagi Hak-Hak Warga Negara. Al-Adl, 9(April), 67–86.
Perwithosuci, W., Anas, M., Hidayah, N., Putri, R., & Hadibasyir, H. (2025). Kualitas lingkungan, pertumbuhan ekonomi dan kepadatan penduduk: Studi panel di Indonesia. Optimal: Jurnal Ekonomi dan Pembangunan . https://doi.org/10.12928/optimum.v15i1.12559
Prasetyanto, P., & Sari, F. (2021). KURVA KUZNETS LINGKUNGAN: PERTUMBUHAN EKONOMI DENGAN DEGRADASI LINGKUNGAN DI INDONESIA. Jurnal Internasional Ekonomi dan Kebijakan Energi . https://doi.org/10.32479/ijeep.11609
Rahayu, H., & Handri, H. (2023). PENGARUH KUALITAS LINGKUNGAN TERHADAP PEMBANGUNAN BERKELANJUTAN DI INDONESIA. Jurnal Alwatzikhoebillah : Kajian Islam, Pendidikan, Ekonomi, Humaniora . https://doi.org/10.37567/alwatzikhoebillah.v9i1.1633 w
Setiawan, M., & Primandhana, W. (2022). Analisis pengaruh beberapa sektor PDRB terhadap indeks kualitas lingkungan hidup di Indonesia. KINERJA . https://doi.org/10.30872/jkin.v19i1.10830
Sumargo, B., & Haida, R. (2020). Keterkaitan Pertumbuhan Ekonomi, Kemiskinan dan Kualitas Lingkungan di Indonesia. Jurnal Ekonomi Pembangunan: Kajian Masalah Ekonomi dan Pembangunan , 21, 47-59. https://doi.org/10.23917/jep.v21i1.8262
Wafa, A. (2025). Analisis determinan kualitas lingkungan hidup di Indonesia. Tinjauan Ekonomi, Keuangan, dan Bisnis . https://doi.org/10.20885/efbr.vol1.iss1.art1
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Santa Barbara: Kluwer Academic Publishers.
Arbia, G. (2014). A Primer For Spatial Econometrics: With Applications in R, STATA and Python. Rome: Palgrave Macmillan.
Astaiza-Gomez, J. G. (2020). Lagrange Multiplier Tests in Applied Research. Journal de Ciencia e Ingeniería, 12(1), 13-19. doi:https://doi.org/10.46571/JCI.2020.1.2
Bai, J., Duan, J., & Han, X. (2024). The likelihood ratio test for structural changes in factor models. Journal of Econometrics(238).
Bivand, R. S., Pebesma, E. J., & Gómez-Rubio, V. (2008). Applied Spatial Data Analysis with R. 8 Springer Science. doi:10.1007/978-0-387-78171-6.
Chen, Y. (2020). New framework of Getis-Ord’s indexes associating spatial autocorrelation with interaction. PLoS ONE, 15(7).
Elhorst, J. P. (2022). The dynamic general nesting spatial econometric model for spatial panels with common factors: Further raising the bar. Rev Reg Res, 42, 249-267. doi:https://doi.org/10.1007/s10037-021-00163-w
Kholifia, N., Muksar, M., Atikah, N., & Afifah, D. L. (2021). Spatial analysis of factors influencing Gross Regional Domestic Product (GRDP) in East Java: a spatial durbin error model analysis. Journal of Physics: Conference Series. doi:10.1088/1742-6596/1918/4/042044
Verbeek, M. (2004). A Guide to Modern Econometrics. John Wiley & Sons, LTD.
Yamada, H. (2024). Geary’s c for Multivariate Spatial Data. Mathemathics, 1-12.
Comber, A., & Zeng, W. (2019). Spatial interpolation using areal features: A review of methods and opportunities using new forms of data with coded illustrations. Geography Compass, 13(10). https://doi.org/10.1111/gec3.12465
Mojtaba Zaresefat, Reza Derakhshani, & Griffioen, J. (2024). Empirical Bayesian Kriging, a Robust Method for Spatial Data Interpolation of a Large Groundwater Quality Dataset from the Western Netherlands. Water, 16(18), 2581–2581. https://doi.org/10.3390/w16182581
Ikechukwu, M. N., Ebinne, E., Idorenyin, U., & Raphael, N. I. (2017). Accuracy assessment and comparative analysis of IDW, spline and kriging in spatial interpolation of landform (topography): an experimental study. Journal of Geographic Information System, 9(3), 354-371.
FILE PETA https://drive.google.com/drive/folders/15nmgQnOc58jNLIaXt10To6aZyqdRigT2?usp=sharing