Pemetaan Pola Sebaran dan Pemodelan Faktor Risiko Stunting di Indonesia Menggunakan Pendekatan Spatial Econometrics
Ara Figura Bilbina - 140610230004 - Statistika Universitas Padjadjaran 2025
2025-12-22
BAB I : Pendahuluan
1.1 Latar Belakang
Stunting masih menjadi isu strategis di Indonesia yang ditandai dengan tingginya disparitas prevalensi antarprovinsi. Secara substantif, kasus ini dipengaruhi oleh faktor multidimensi yang kompleks, meliputi determinan lingkungan (sanitasi layak, akses air bersih), akses layanan kesehatan (kepemilikan JKN, cakupan imunisasi, pemberian TTD), serta kondisi ekonomi makro (tingkat kemiskinan). Heterogenitas karakteristik antarwilayah ini mengindikasikan bahwa masalah stunting tidak hanya bersifat lokal, tetapi juga memiliki keterkaitan geografis.
Dalam analisis statistika, penggunaan regresi linier klasik (Ordinary Least Square/OLS) sering kali kurang tepat untuk data kewilayahan karena mengasumsikan independensi antar-pengamatan. Padahal, wilayah yang bertetangga cenderung memiliki kemiripan pola (spatial dependence), baik karena kesamaan karakteristik sosial-budaya maupun dampak kebijakan lintas daerah (spillover effect). Mengabaikan efek spasial ini dapat menyebabkan estimator menjadi bias dan kesimpulan yang tidak valid.
Oleh karena itu, diperlukan pendekatan Spatial Econometrics untuk mengakomodasi ketergantungan spasial tersebut. Penelitian ini bertujuan memetakan pola sebaran stunting di Indonesia tahun 2024 serta membandingkan performa model spasial, yaitu Spatial Error Model (SEM) dan Spatial Durbin Error Model (SDEM). Melalui komparasi ini, diharapkan diperoleh model terbaik yang mampu menjelaskan determinan stunting secara presisi dan menghasilkan rekomendasi kebijakan berbasis wilayah yang lebih akurat.
1.2 Identifikasi Masalah
Pemetaan persentase stunting di Indonesia pada tahun 2024 menunjukkan adanya indikasi tidak independen secara spasial, dimana nilai stunting di satu provinsi memiliki keterkaitan dengan nilai observasi di provinsi lain yang berdekatan. Kondisi ini menunjukkan adanya spatial dependence. Adanya dependensi spasial ini mengindikasikan bahwa asumsi independensi antar observasi yang digunakan pada model regresi klasik (OLS) tidak terpenuhi. Jika kondisi tersebut diabaikan, maka hasil estimasi dapat menjadi bias dan tidak efisien. Oleh karena itu, perlu digunakan pendekatan ekonometrika spasial untuk menangkap pengaruh spasial secara lebih akurat.
1.3 Maksud Tujuan Penelitian
Maksud dari penelitian ini adalah untuk menerapkan pendekatan spasial untuk menganalisis faktor yang mempengaruhi stunting di Indonesia pada tahun 2024. Melalui analisis spasial, penelitian ini bertujuan untuk memperoleh model terbaik yang dapat menjelaskan hubungan antarwilayah serta faktor-faktor utama yang berpengaruh terhadap stunting di Indonesia.
1.4 Batasan Penelitian
Penelitian ini dibatasi pada wilayah Repiblik Indonesia dengan unit analisis Provinsi. Analisis ini menggunakan data 34 provinsi, dikarenakan data dari provinsi baru seperti Papua Pegunungan, Papua Tengah , Papua Selatan, dan Papua Barat Daya belum tersedia. Analisis difokuskan pada pola spasial stunting tanpa membahas aspek penyebab non-spasial secara mendalam. Selain itu, penelitian ini hanya menggunakan data tahun 2024 sehingga hasilnya menggambarkan kondisi pada periode tertentu dan tidak mencakup dinamika perubahan antarwaktu.
BAB II : Tinjauan Pustaka
2.1 Spatial Autocorrelation
Dalam analisis spasial, keterkaitan antar wilayah secara spasial diukur menggunakan Moran’s I. Nilai I > 0 menunjukkan adanya autokorelasi spasial positif, yaitu wilayah dengan nilai mirip (high-high atau low-low) cenderung berdekatan. Sebaliknya, I < 0 menandakan pola checkerboard (high-low), sedangkan I ≈ 0 menunjukkan pola yang acak. Dalam perhitungannya, Moran’s I memanfaatkanmatriks bobot spasial (W) yang menentukan seberapa kuat pengaruh satu wilayah terhadap wilayah lain berdasarkan kriteria kedekatan tertentu, yakni Rook, Queen, ataupun Jarak. Selain Moran’s I, terdapat metode pengukuran lain seperti Greary’s C untuk ukuran Global, serta Local Moran’s I dan Getis-Ord G* untuk ukuran lokal.
2.2 OLS
OLS adalah metode estimasi parameter untuk regresi klasik. Metode ini membutuhkan pemenuhan pengujian asumsi klasik, diantaranya pengujian multikolinearitas, residual berdistribusi normal, dan homoskedastisitas. Metode ini mengasumsikan bahwa error antarwilayah bersifat independen dan homogen, sehingga tidak dapat menangkap spillover atau pengaruh kemiskinan dari satu daerah ke daerah lain. Hal ini menyebabkan hasil estimasi menjadi bias dan kurang akurat. Untuk mengatasi keterbatasan tersebut, digunakan model ekonometrika spasial yang mampu memperhitungkan spatial dependence antarwilayah. Adapun model matematis OLS adalah sebagai berikut.
\[\begin{equation} Y_i = \beta_0 + \sum_{k=1}^{K} \beta_k X_{ki} + \varepsilon_i, \quad \varepsilon_i \sim N(0,\sigma^2) \end{equation}\]
2.3 Uji Lagrange Multipier
Uji ini digunakan untuk mengidentifikasi ada tidaknya dependensi spasial pada model. Jika hasil uji menunjukkan adanya dependensi lag, maka digunakan model Spatial Autoregressive (SAR). Jika yang muncul adalah dependensi galat, maka digunakan model Spatial Error Model (SEM). Uji ini membantu menentukan model spasial yang paling sesuai sebelum analisis lanjut.
2.5 Spatial Error Model (SEM)
Spatial Error Model (SEM) merupakan model ekonometrika spasial yang mengasumsikan bahwa dependensi spasial tidak terjadi antarvariabel dependen secara langsung, melainkan melalui komponen galat (error) akibat adanya variabel-variabel yang tidak teramati (unobserved variables) yang memiliki pola spasial. Kondisi ini sering dijumpai ketika terdapat faktor-faktor penting yang bersifat lokal (misalnya budaya atau kondisi geografis mikro) yang mempengaruhi variabel \(Y\) namun datanya tidak tersedia dalam model, dan faktor-faktor tersebut berkorelasi antar wilayah.
Secara matematis, SEM dirumuskan sebagai:
\[Y = X\beta + u\] \[u = \lambda W u + \epsilon\]
Di mana \(\lambda\) adalah koefisien error spasial yang mengukur kekuatan korelasi antar-galat wilayah.
Model SEM digunakan ketika hasil Uji Lagrange Multiplier (LM) menunjukkan adanya dependensi pada error (\(LM_{error}\) signifikan). Mengabaikan dependensi ini dapat menyebabkan estimator varians menjadi bias dan pengujian hipotesis (Uji-t) menjadi tidak valid, meskipun estimasi koefisien regresi (\(\beta\)) tetap tidak bias (unbiased). Dengan demikian, SEM bertujuan untuk memperbaiki efisiensi dan validitas inferensi statistik.
2.6 Spatial Durbin Error Model (SDEM)
Spatial Durbin Error Model (SDEM) merupakan pengembangan dari model SEM yang memperluas analisis dengan memasukkan pengaruh lag spasial dari variabel independen (\(WX\)). Model ini mengasumsikan bahwa variabel dependen (\(Y\)) dipengaruhi oleh variabel independen wilayah itu sendiri (\(X\)), variabel independen dari wilayah tetangga (\(WX\)), serta adanya dependensi spasial pada komponen galat (error). Berbeda dengan model Spatial Durbin Model (SDM) yang memiliki efek spillover global, SDEM hanya menangkap efek spillover lokal, di mana perubahan pada suatu wilayah hanya memengaruhi tetangga langsungnya tanpa efek umpan balik yang berulang (feedback loop).
Secara matematis, persamaan model SDEM dapat dituliskan sebagai berikut:
\[Y = X\beta + WX\theta + u\] \[u = \lambda W u + \epsilon\]
Di mana: * \(\theta\) adalah vektor koefisien regresi untuk variabel independen yang di-lag secara spasial (pengaruh karakteristik tetangga terhadap \(Y\) amatan). * \(\lambda\) adalah koefisien error spasial.
SDEM dianggap sebagai model yang kuat karena mampu mengakomodasi heterogenitas spasial melalui variabel independen tetangga sekaligus menangani autokorelasi spasial yang muncul dari faktor-faktor yang tidak teramati (unobserved variables) pada komponen galat. Dalam konteks stunting, model ini memungkinkan peneliti untuk melihat apakah perbaikan sanitasi atau akses kesehatan di provinsi tetangga (\(WX\)) turut berdampak pada penurunan stunting di provinsi amatan, sambil tetap mengontrol bias akibat variabel pengganggu yang tidak terukur.
2.7 Akaike’s Information Criterion (AIC)
Menurut Suryaningrat, S. H. (2021) Akaike Information Criterion (AIC) merupakan metode analisis yang digunakan untuk memperoleh model yang terbaik dengan menggunakan estimasi maximum likelihood sebagai perhitungan yang sesuai. AIC mempertimbangkan keseimbangan antara kompleksitas model (jumlah parameter) dan kecocokan model terhadap data (likelihood). AIC dihitung dengan rumus:
\[\begin{equation} \text{AIC} = -2 \ln(\hat{L}) + 2k, \end{equation}\]
Model dengan nilai AIC yang lebih rendah dianggap lebih baik karena menunjukkan bahwa model tersebut memiliki kompleksitas yang wajar tetapi tetap mampu menjelaskan data dengan baik.
BAB III : Metodologi Penelitian
3.1 Data
Penelitian ini menggunakan data sekunder yang diperoleh dari Badan Pusat Statistik (BPS) tahun 2024 dengan unit analisis sebanyak 34 Provinsi di Indonesia. Data yang digunakan bersifat cross section dan mencerminkan kondisi kesehatan dan sosial yang diduga berpengaruh terhadap tingkat stunting di Indonesia. Berikut variabel yang diamati dalam analisis ini:
| Variabel | Satuan |
|---|---|
| Prevalensi balita stunting berdasarkan provinsi di indonesia (Y) | (%) |
| Persentase penduduk miskin (X1) | (%) |
| Persentase Rumah Tangga yang Memiliki Akses terhadap Sanitasi Layak (X2) | (%) |
| Persentase Rumah Tangga yang Memiliki Akses terhadap Sumber Air Minum Layak Menurut Provinsi (X3) | (%) |
| Persentase Anak Umur 0-59 Bulan (Balita) yang Pernah Mendapat Imunisasi Campak (X4) | (%) |
| Proporsi Balita Mengalami Diare dalam 1 Bulan Terakhir (X5) | (%) |
| Proporsi Balita Mengalami ISPA dalam 1 Bulan Terakhir (X6) | |
| Proporsi Pemanfaatan Jaminan Kesehatan (X7) | (%) |
| Proporsi Ibu Balita Berdasarkan Usia Kehamilan Pertama Kali Mendapat/Diberi/Membeli Tablet Tambah Darah di trisemester 1 | (%) |
Selain data atribut, penelitian ini juga menggunakan data spasial berupa peta batas administrasi provinsi di Indonesia yang diperoleh melalui situs https://gadm.org/download_country.html dalam format shapefile (.shp). Data ini digunakan untuk membentuk hubungan antarwilayah dalam analisis spasial.
3.2 Metode Analisis
Analisis dilakukan menggunakan pendekatan statistik spasial dengan bantuan perangkat lunak R. Pendekatan ini bertujuan untuk mengidentifikasi pengaruh faktor kesehatan-sosial terhadap stunting serta mendeteksi adanya ketergantungan spasial antarwilayah di Indonesia. Analisis diawali dengan uji autokorelasi spasial global (Moran’s I, Geary’s C) dan lokal (LISA serta Getis-Ord Gi*) untuk melihat pola keterkaitan antarwilayah.
Selanjutnya digunakan model regresi OLS sebagai dasar analisis hubungan antarvariabel, diikuti uji Lagrange Multiplier (LM) untuk menentukan model spasial yang tepat. Beberapa model yang dipertimbangkan mencakup SEM dan SDEM. Pemilihan model terbaik dilakukan berdasarkan nilai AIC, signifikansi parameter spasial (ρ dan λ), serta kesesuaian teori yang mendukung hubungan spasial antarwilayah.
3.3 Alur Kerja Penelitian
- Persiapan data
- Analisis statistika deskriptif
- Estimasi parameter menggunakan metode OLS
- Uji asumsi klasik model OLS
- Uji autokorelasi spasial
- Uji lagrange multiplier
- Estimasi parameter menggunakan model terpilih ( SEM dan SDEM)
- Uji asumsi model spasial
- Pemilihan model terbaik berdasarkan AIC
- Uji sensitivitas bobot spasial
- Interpolasi spasial
- Interpretasi hasil
BAB IV : Hasil dan Pembahasan
4.1 Analisis Deskriptif
Analisis deskriptif dilakukan untuk memberikan gambaran umum mengenai karakteristik data penelitian. Hasil analisis deskriptif data dapat dilihat pada tabel berikut:
| Variabel | Rata-rata | Minimum | Maximum | Simpangan Baku |
|---|---|---|---|---|
| Stunting (Y) | 21.92 | 22.05 | 37.00 | 5.92 |
| Kemiskinan (X1) | 9.56 | 4.00 | 21.66 | 4.56 |
| Akses Sanitasi Layak (X2) | 85.02 | 72.82 | 96.83 | 6.02 |
| Sumber Air Minum Layak (X3) | 89.42 | 72.10 | 99.96 | 6.75 |
| Balita yang Pernah Mendapat Imunisasi Campak (X4) | 72.55 | 41.44 | 82.93 | 7.48 |
| Balita Mengalami Diare (X5) | 4.09 | 1.30 | 8.50 | 1.97 |
| Balita Mengalami ISPA (X6) | 30.24 | 16.80 | 49.10 | 8.89 |
| Penduduk yang memiliki dan memanfaatkan JKN/Jamkesda (X7) | 44.09 | 0.00 | 71.00 | 16.00 |
| Ibu balita yang mendapatkan TTD sejak trisemester 1 (X8) | 8.57 | 4.00 | 13.10 | 2.65 |
Hasil statistik deskriptif menunjukkan bahwa prevalensi stunting memiliki rata-rata sebesar 21,92% dengan variasi yang cukup tinggi (SD 5,92) dan nilai maksimum mencapai 37,00%, menandakan masih adanya wilayah dengan kondisi stunting yang tergolong kritis. Dari sisi sosial-ekonomi, tingkat kemiskinan rata-rata sebesar 9,56% dengan disparitas yang lebar antar wilayah, sementara capaian infrastruktur dasar relatif lebih baik dengan rata-rata akses sanitasi layak dan sumber air minum layak masing-masing sebesar 85,02% dan 89,42%, meskipun masih terdapat wilayah dengan akses yang rendah. Pada aspek kesehatan balita, cakupan imunisasi campak menunjukkan ketimpangan yang besar, prevalensi diare relatif rendah, namun angka ISPA tergolong tinggi dengan rata-rata 30,24%, yang berpotensi menghambat pertumbuhan anak. Selain itu, indikator akses layanan kesehatan memperlihatkan kesenjangan yang signifikan, terutama pada kepemilikan jaminan kesehatan dan cakupan pemberian tablet tambah darah pada ibu sejak trimester pertama yang masih sangat rendah (rata-rata 8,57%), sehingga berpotensi menjadi faktor risiko utama dalam kejadian stunting di Indonesia.
Prevalensi stunting di
Indonesia masih tergolong tinggi dengan rata-rata 21,92% dan menunjukkan
ketimpangan spasial yang tajam, di mana wilayah Indonesia
Timur—khususnya Papua dan Nusa Tenggara Timur—mendominasi zona
prevalensi tertinggi, sementara wilayah Barat cenderung lebih rendah.
Pola ini sejalan dengan sebaran kemiskinan yang terkonsentrasi di
wilayah timur, menegaskan kemiskinan sebagai determinan struktural utama
stunting. Ketimpangan infrastruktur dasar juga terlihat jelas, di mana
akses sanitasi dan air minum layak relatif tinggi di Jawa dan Bali,
namun masih rendah di Papua, Kalimantan, dan Bengkulu. Dari sisi
kesehatan balita, cakupan imunisasi campak menunjukkan variasi antar
wilayah, sementara hotspot diare teridentifikasi di Aceh dan Nusa
Tenggara Barat serta beban ISPA yang tinggi tersebar di Nusa Tenggara
Timur, Jawa Timur, dan Bengkulu. Pola kepemilikan jaminan kesehatan
menunjukkan anomali spasial dengan cakupan tinggi di wilayah ujung barat
dan timur, namun rendah di beberapa wilayah Sumatera dan Maluku,
sedangkan konsumsi Tablet Tambah Darah pada ibu hamil secara umum masih
rendah secara nasional dan belum merata, bahkan di wilayah dengan
fasilitas kesehatan yang relatif baik.
4.2 Model OLS
\[ \begin{aligned} \hat{Y}_i =\;& 51.196^{*} + 0.441\,X_{1}^{*} - 0.036\,X_{2} - 0.163\,X_{3} - 0.244\,X_{4} + 0.924\,X_{5} \\ & - 0.335\,X_{6} + 0.135\,X_{7}^{*} + 0.263\,X_{8} \end{aligned} \]
Model regresi OLS yang terbentuk secara simultan berpengaruh signifikan terhadap prevalensi stunting (\(F\text{-statistic } p\text{-value} < 0.01\)) dengan kemampuan menjelaskan variabilitas data sebesar 53,99% (\(R^2 = 0.5399\)). Secara parsial, variabel Kemiskinan (\(X_1\)) dan Kepemilikan JKN (\(X_7\)) terbukti memiliki pengaruh positif yang signifikan pada taraf nyata 5%, di mana setiap peningkatan 1% kemiskinan dan kepemilikan JKN diasosiasikan dengan kenaikan stunting masing-masing sebesar 0,44% dan 0,14%, sementara variabel Imunisasi Campak (\(X_4\)) menunjukkan indikasi pengaruh negatif (menurunkan stunting) yang signifikan pada taraf 10%.
4.2.1 Uji Asumsi OLS
| Uji Asumsi | p-value | Ket |
|---|---|---|
| Normalitas | 0.94 | Normalitas terpenuhi |
| Homoskedastisitas | 0.5844 | Homoskedastisitas terpenuhi |
| Linearitas | 0.2219 | Linearitas terpenuhi |
| Uji Moran’s I Residual | 0.01553 | Terdapat autokorelasi spasial pada error |
| Variabel | VIF |
|---|---|
| Kemiskinan (X1) | 1.34 |
| Air layak (X2) | 2.10 |
| Sumber air minum layak (X3) | 1.73 |
| Imunisasi campak (X4) | 1.71 |
| Diare (X5) | 5.05 |
| ISPA (X6) | 7.27 |
| JKN/Jamkesda (X7) | 1.68 |
| Tablet tambah darah (X8) | 2.45 |
Berdasarkan hasil uji asumsi klasik yang disajikan pada Tabel Uji Asumsi, dapat diketahui bahwa model regresi OLS telah memenuhi sebagian besar asumsi dasar regresi linier. Uji normalitas menunjukkan nilai p-value sebesar 0,94, yang lebih besar dari tingkat signifikansi 5%, sehingga residual berdistribusi normal. Selanjutnya, uji homoskedastisitas menghasilkan p-value sebesar 0,5844, yang mengindikasikan bahwa varians residual bersifat konstan. Uji linearitas juga menunjukkan p-value sebesar 0,2219, sehingga hubungan antara variabel independen dan variabel dependen dapat dianggap linier.
Namun demikian, hasil uji Moran’s I pada residual menunjukkan p-value sebesar 0,01553, yang lebih kecil dari tingkat signifikansi 5%. Hal ini menandakan adanya autokorelasi spasial yang signifikan pada error model, sehingga asumsi independensi residual tidak terpenuhi. Dengan demikian, meskipun asumsi klasik regresi linier telah terpenuhi, model OLS belum sepenuhnya memadai dan perlu dikembangkan menggunakan model regresi spasial untuk mengakomodasi ketergantungan spasial antar wilayah.
4.3 Autokorelasi Spasial
4.3.1 Moran’s I
| Variabel | Moran’s I | p-permutasi | Ket |
|---|---|---|---|
| Stunting (Y) | 0.30 | 0.023 | Ada autokorelasi spasial |
| Kemiskinan (X1) | 0.39 | 0.007 | Ada autokorelasi spasial |
| Air layak (X2) | -0.02 | 0.483 | Tidak ada autokorelasi spasial |
| Sumber air minum layak (X3) | 0.43 | 0.003 | Ada autokorelasi spasial |
| Imunisasi campak (X4) | 0.17 | 0.071 | Tidak ada autokorelasi spasial |
| Diare (X5) | 0.60 | 0.001 | Tidak ada autokorelasi spasial |
| ISPA (X6) | 0.68 | 0.001 | Ada autokorelasi spasial |
| JKN/Jamkesda (X7) | 0.09 | 0.209 | Tidak ada autokorelasi spasial |
| Tablet tambah darah (X8) | 0.61 | 0.001 | Ada autokorelasi spasial |
4.3.2 Geary’s C
| Variabel | Geary’s C | p-value | Ket |
|---|---|---|---|
| Stunting (Y) | 0.22 | 0.001 | Ada autokorelasi spasial |
| Kemiskinan (X1) | 0.34 | 0.002 | Ada autokorelasi spasial |
| Air layak (X2) | 0.77 | 0.302 | Tidak ada autokorelasi spasial |
| Sumber air minum layak (X3) | 0.43 | 0.004 | Ada autokorelasi spasial |
| Imunisasi campak (X4) | 0.53 | 0.127 | Tidak ada autokorelasi spasial |
| Diare (X5) | 0.22 | 0.001 | Tidak ada autokorelasi spasial |
| ISPA (X6) | 0.24 | 0.001 | Ada autokorelasi spasial |
| JKN/Jamkesda (X7) | 0.78 | 0.330 | Tidak ada autokorelasi spasial |
| Tablet tambah darah (X8) | 0.25 | 0.001 | Ada autokorelasi spasial |
4.3.3 Local Moran’s I (LISA)
4.3.4 Getis-Ord*
4.4 Uji Lagrange Multiplier
| Statistik Uji | Value | p-value | Ket |
|---|---|---|---|
| LM | |||
| Error (RSerr) | 4.58 | 0.03 | Sig |
| LM | |||
| Lag (RSlag) | 0.03 | 0.85 | not.sig |
| Robust | |||
| LM Error (AdjRSerr) | 6.55 | 0.01 | sig |
| Robust | |||
| LM Lag (AdjRSlag) | 2.01 | 0.15 | not.sig |
Hasil uji Lagrange Multiplier menunjukkan bahwa LM Error (RSerr) dan Robust LM Error (AdjRSerr) signifikan, masing-masing dengan p-value sebesar 0,03 dan 0,01, sedangkan LM Lag (RSlag) dan Robust LM Lag (AdjRSlag) tidak signifikan. Temuan ini mengindikasikan bahwa ketergantungan spasial dalam data terutama berasal dari komponen error, bukan dari pengaruh langsung variabel dependen antar wilayah.
Berdasarkan hasil tersebut, Spatial Error Model (SEM) dipilih sebagai model utama, karena SEM secara khusus dirancang untuk menangkap autokorelasi spasial yang muncul pada error term. Penggunaan SEM diharapkan mampu mengoreksi pelanggaran asumsi independensi residual yang terdeteksi pada model OLS melalui uji Moran’s I.
Selain SEM, juga dipertimbangkan Spatial Durbin Error Model (SDEM) sebagai pengembangan model. Pemilihan SDEM didasarkan pada pertimbangan bahwa meskipun ketergantungan spasial utama terletak pada error, variabel independen di wilayah sekitar (spatially lagged independent variables) berpotensi memengaruhi variabel dependen suatu wilayah. Dengan demikian, SDEM memungkinkan analisis yang lebih komprehensif karena mengakomodasi spillover effect dari variabel penjelas antar wilayah, tanpa mengasumsikan adanya lag pada variabel dependen.
4.5 Model SEM
\[ \begin{aligned} Y_i =\;& 53.5004^{*} + 0.6475\,X_{1i}^{*} - 0.0540\,X_{2i} - 0.2400\,X_{3i}^{*} - 0.0580\,X_{4i} + 1.4597\,X_{5i}^{*} \\ & - 0.4233\,X_{6i}^{*} + 0.0915\,X_{7i}^{*} - 0.4643\,X_{8i} \end{aligned} \]
Secara keseluruhan, hasil estimasi menunjukkan bahwa SEM mampu mengakomodasi autokorelasi spasial pada residual, ditunjukkan oleh parameter λ yang signifikan serta peningkatan goodness-of-fit model. Dengan demikian, SEM lebih tepat digunakan dibandingkan OLS untuk menganalisis hubungan antara variabel independen dan variabel dependen dalam konteks data spasial.
4.6 Model SDEM
\[ \begin{aligned} Y_i =\;& 70.6776^{*} + 0.5248\,X_{1i}^{*} - 0.1548\,X_{2i} - 0.2702\,X_{3i}^{*} - 0.1640\,X_{4i} + 0.9739\,X_{5i} \\ & - 0.2588\,X_{6i} + 0.0915\,X_{7i}^{*} - 0.2375\,X_{8i} - 0.5808\,(W X_{1})_i^{*} - 0.1836\,(W X_{2})_i \\ & + 0.1611\,(W X_{3})_i - 0.3483\,(W X_{4})_i - 0.0444\,(W X_{5})_i + 0.2258\,(W X_{6})_i \\ & + 0.0079\,(W X_{7})_i + 1.5680\,(W X_{8})_i^{*} \end{aligned} \]
Hasil estimasi Spatial Durbin Error Model (SDEM) menunjukkan bahwa variabel internal wilayah tetap berperan penting dalam memengaruhi variabel dependen, sekaligus mengungkap adanya spillover spasial antar wilayah. Beberapa variabel independen berpengaruh signifikan secara langsung, sementara variabel di wilayah sekitar juga memberikan pengaruh tidak langsung yang signifikan, baik positif maupun negatif. Parameter spasial error yang signifikan mengindikasikan bahwa ketergantungan spasial masih ada, namun sebagian telah dijelaskan oleh keberadaan lag variabel independen. Dengan demikian, SDEM memberikan pemahaman yang lebih komprehensif terhadap mekanisme pengaruh spasial, meskipun digunakan sebagai model pelengkap dibandingkan SEM.
4.7 Uji Asumsi Model
Uji Asumsi Model SEM
| Uji | Statistik | p-value | Ket |
|---|---|---|---|
| Residual autokorelasi | 0.052 | 0.314 | Tidak ada autokorelasi residual |
| Normalitas | 0.975 | 0.612 | Normalitas terpenuhi |
| Homoskedastisitas | 1.7004e - 05 | 0.996 | Homoskedastisitas terpenuhi |
Uji Asumsi Model SDEM
| Uji | Statistik | p-value | Ket |
|---|---|---|---|
| Residual autokorelasi | -0.077 | 0.592 | Tidak ada autokorelasi residual |
| Normalitas | 0.977 | 0.687 | Normalitas terpenuhi |
| Homoskedastisitas | 1.804 | 0.179 | Homoskedastisitas terpenuhi |
Dengan terpenuhinya seluruh asumsi-asumsi diatas, baik SEM maupun SDEM layak digunakan secara statistik untuk analisis lebih lanjut. Pemilihan model akhir selanjutnya dapat didasarkan pada kecocokan model (AIC).
4.8 Perbandingan Model
| Model | df | LogLik | AIC |
|---|---|---|---|
| OLS | 10 | -95.012 | 210.025 |
| SEM | 11 | -89.840 | 201.680 |
| SDEM | 19 | -84.246 | 206.493 |
Berdasarkan perbandingan nilai Log-Likelihood dan Akaike Information Criterion (AIC), model Spatial Error Model (SEM) dipilih sebagai model terbaik dalam penelitian ini. Meskipun model SDEM memiliki nilai Log-Likelihood tertinggi, jumlah parameter yang lebih besar menyebabkan nilai AIC SDEM lebih tinggi dibandingkan SEM. Sementara itu, model OLS menunjukkan nilai Log-Likelihood terendah dan AIC tertinggi, sehingga memiliki kecocokan model yang paling rendah.
Dengan nilai AIC terendah sebesar 201,68, SEM memberikan keseimbangan terbaik antara kecocokan model dan kompleksitas parameter, serta telah terbukti mampu mengakomodasi autokorelasi spasial pada residual. Oleh karena itu, SEM dipilih sebagai model akhir yang paling tepat untuk menjelaskan hubungan antar variabel dalam konteks data spasial penelitian ini.
4.9 Uji Sensitivitas Bobot
| Bobot | rho | AIC | LogLik |
|---|---|---|---|
| Queen | 0.6147 | 201.6806 | -89.840 |
| Rook | 0.6147 | 201.6806 | -89.840 |
| KNN-4 | 0.6884 | 204.6052 | -91.302 |
Berdasarkan hasil uji sensitivitas bobot spasial pada model Spatial Error Model (SEM), diperoleh bahwa penggunaan bobot Queen contiguity dan Rook contiguity menghasilkan nilai parameter spasial (ρ), AIC, dan Log-Likelihood yang identik. Hal ini menunjukkan bahwa hasil estimasi model stabil dan tidak sensitif terhadap perbedaan jenis bobot kontiguitas berbasis batas wilayah.
Sebaliknya, penggunaan bobot K-nearest neighbors (KNN-4) menghasilkan nilai ρ yang lebih besar, namun disertai dengan nilai AIC yang lebih tinggi dan Log-Likelihood yang lebih rendah dibandingkan bobot Queen dan Rook. Temuan ini mengindikasikan bahwa meskipun ketergantungan spasial terdeteksi lebih kuat pada bobot KNN-4, kecocokan model secara keseluruhan menjadi lebih rendah.
Dengan demikian, bobot spasial Queen contiguity (atau Rook contiguity) dipilih sebagai bobot terbaik dalam penelitian ini karena memberikan kecocokan model paling optimal dan hasil estimasi yang paling stabil.
4.10 Spasial Interpolasi
Spasial interpolasi merupakan metode analisis spasial yang digunakan untuk memperkirakan nilai suatu variabel pada lokasi yang tidak memiliki data pengamatan, berdasarkan nilai dari lokasi-lokasi di sekitarnya. Metode ini berangkat dari prinsip kedekatan geografis, yaitu bahwa lokasi yang berdekatan cenderung memiliki karakteristik yang lebih mirip dibandingkan lokasi yang berjauhan. Dalam analisis ini, interpolasi digunakan untuk menggambarkan pola sebaran spasial variabel penelitian secara kontinu, sehingga informasi spasial dapat disajikan tidak hanya pada titik atau wilayah pengamatan, tetapi juga pada area di antaranya.
Spasial interpolasi dilakukan untuk memvisualisasikan sebaran nilai prediksi hasil Spatial Error Model (SEM) secara kontinu antar wilayah, sehingga pola spasial yang tidak terlihat pada data diskrit dapat ditangkap dengan lebih jelas. Metode Inverse Distance Weighting (IDW) digunakan dengan asumsi bahwa wilayah yang berdekatan memiliki karakteristik yang serupa, sehingga nilai prediksi pada suatu lokasi dipengaruhi oleh wilayah di sekitarnya. Hasil interpolasi prediksi SEM disajikan pada Gambar berikut, yang menunjukkan variasi spasial nilai prediksi antar wilayah di Indonesia.
Peta interpolasi prediksi SEM menunjukkan adanya ketimpangan spasial nilai prediksi antar wilayah di Indonesia. Wilayah Indonesia bagian barat, khususnya sebagian Sumatra dan Jawa, didominasi oleh nilai prediksi rendah hingga menengah, sementara wilayah Indonesia bagian timur—terutama Papua, Maluku, dan Nusa Tenggara—menunjukkan nilai prediksi yang lebih tinggi. Pola sebaran yang halus dan berkelanjutan mengindikasikan bahwa ketergantungan spasial telah berhasil ditangkap oleh model SEM, sejalan dengan hasil uji Moran’s I residual yang menunjukkan tidak adanya autokorelasi spasial pada error model.
4.10 Interpretasi Hasil
Berdasarkan visualisasi hasil estimasi Spatial Error Model (SEM), terlihat adanya pola ketimpangan spasial yang signifikan terkait persentase stunting di Indonesia antara wilayah Barat dan Timur. Peta prediksi memperlihatkan bahwa prevalensi stunting tertinggi terkonsentrasi di wilayah Papua dan Papua Barat dengan kisaran 26,49% hingga 31,36%, sedangkan wilayah Sumatera dan sebagian besar Jawa cenderung memiliki tingkat stunting yang lebih rendah (12,20% – 18,98%). Dari sisi evaluasi model, peta residual standar mengindikasikan bahwa model SEM memiliki performa yang cukup baik (goodness of fit) dalam menjelaskan data di mayoritas provinsi, ditandai dengan dominasi warna pudar yang berarti selisih antara nilai riil dan prediksi relatif kecil. Meskipun demikian, teridentifikasi adanya outlier spasial yang perlu perhatian khusus; wilayah Sulawesi menunjukkan residual positif ekstrem (merah) yang mengindikasikan angka stunting aktual jauh lebih tinggi dibandingkan prediksi model (underestimation), sementara beberapa wilayah di Sumatera bagian selatan menunjukkan residual negatif ekstrem (biru tua) di mana kondisi aktualnya justru lebih baik daripada yang diperkirakan oleh model (overestimation).
Hasil estimasi menggunakan Spatial Durbin Error Model (SDEM) kembali mempertegas adanya ketimpangan spasial yang signifikan terkait prevalensi stunting di Indonesia. Berdasarkan peta prediksi, wilayah dengan tingkat kerentanan stunting tertinggi (zona kuning dan oranye) terkonsentrasi secara mencolok di wilayah Nusa Tenggara, sebagian besar Sulawesi, dan Papua dengan estimasi persentase stunting mencapai kisaran 25,93% hingga 31,69%. Sebaliknya, wilayah Indonesia bagian barat, khususnya Sumatera dan Jawa, didominasi oleh warna gelap yang mengindikasikan angka prediksi stunting terendah, berkisar antara 13,55% hingga 17,11%. Dari sisi diagnostik, peta residual standar SDEM memperlihatkan performa model yang cukup baik, ditandai dengan dominasi warna pudar (krem dan biru muda) di sebagian besar wilayah, yang berarti selisih antara nilai aktual dan prediksi relatif minim. Namun, model mendeteksi adanya anomali spasial di beberapa titik; terdapat wilayah di Sumatera bagian selatan dengan residual negatif (biru) yang menunjukkan kondisi lapangan lebih baik daripada prediksi, serta beberapa area di Kalimantan dan Sulawesi dengan residual positif (oranye) yang menandakan angka stunting aktual justru lebih tinggi dibandingkan yang diekspektasikan oleh variabel-variabel dalam model.
BAB V : Kesimpulan dan Saran
Hasil analisis ini menunjukkan bahwa prevalensi stunting di Indonesia memiliki pola sebaran spasial yang tidak merata, dengan ketimpangan yang jelas antara wilayah Indonesia bagian Barat dan Timur, di mana Papua, Nusa Tenggara, dan sebagian Sulawesi cenderung memiliki tingkat stunting lebih tinggi dibandingkan Sumatera dan Jawa. Hasil uji Moran’s I dan Lagrange Multiplier mengonfirmasi adanya ketergantungan spasial yang dominan pada komponen error, sehingga Spatial Error Model (SEM) dipilih sebagai model terbaik dengan kinerja yang lebih optimal dibandingkan OLS dan SDEM, ditunjukkan oleh nilai AIC yang lebih rendah dan Log Likelihood yang lebih tinggi.
Model SEM terbukti mampu menangkap pengaruh faktor risiko stunting sekaligus mengakomodasi dependensi spasial yang tidak teramati, serta telah memenuhi seluruh asumsi diagnostik model. Berdasarkan temuan ini, disarankan agar kebijakan penanggulangan stunting di Indonesia dirancang secara berbasis wilayah dan kontekstual, dengan prioritas pada daerah berisiko tinggi serta memperhatikan keterkaitan spasial antar provinsi, sehingga intervensi yang dilakukan dapat lebih tepat sasaran dan efektif.
REFERENSI
Revildy, W. D., Lestari, S. S. S., & Nalita, Y. (2020). Pemodelan spatial error model (sem) angka prevalensi balita pendek (stunting) di indonesia tahun 2018. In Seminar Nasional Official Statistics (Vol. 2020, No. 1, pp. 1224-1231).
Singrapati, L. R., & Astuti, E. T. (2024, November). Determinan Prevalensi Stunting di Nusa Tenggara Tahun 2023. In Seminar Nasional Official Statistics (Vol. 2024, No. 1, pp. 183-192).
Eliezer, W. R., & Alwanti, N. (2025). Pengaruh Variabel Sosioekonomi dan Lingkungan Terhadap Prevalensi Stunting Secara Spasial di Provinsi Nusa Tenggara Timur Tahun 2023. Jurnal Statistika Terapan (ISSN 2807-6214), 5(1), 107-126.
Pebrianty, P., Lalli, L., & Embong, M. (2023). Percepatan pencegahan stunting pada anak usia dini dengan pendekatan analisis spasial. Murhum: Jurnal Pendidikan Anak Usia Dini, 4(2), 259-271.
LINK DASHBOARD
LINK VIDEO PRESENTASI
LAMPIRAN
### 1. PACKAGES ###
library(sf)## Warning: package 'sf' was built under R version 4.4.3## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(spdep)## Warning: package 'spdep' was built under R version 4.4.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.4.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.4.3## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(spgwr)## Warning: package 'spgwr' was built under R version 4.4.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.4.3## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)library(tmap)## Warning: package 'tmap' was built under R version 4.4.3library(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(gstat)## Warning: package 'gstat' was built under R version 4.4.3library(raster)## Warning: package 'raster' was built under R version 4.4.3##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## select### 2. MENYIAPKAN DATA ###
# Data atribut
data <- read_excel("~/projekstunting/data/metopen.xlsx")
data## # A tibble: 34 × 11
## No provinsi Y X1 X2 X3 X4 X5 X6 X7 X8
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 Aceh 28.6 14.2 81.1 90.1 41.4 3.4 26.1 26.5 4.5
## 2 2 Sumatera Utara 22 7.99 85.7 92.9 66.3 4.4 34.3 31.4 6.2
## 3 3 Sumatera Barat 24.9 5.97 72.8 86.8 61.2 5.5 42.6 68.6 4.4
## 4 4 Riau 20.1 6.67 86.3 92.2 68.6 5 36 46.3 6
## 5 5 Jambi 17.1 7.1 84.0 82.2 63.6 2.6 25.2 42.8 7.1
## 6 6 Sumatera Selatan 15.9 11.0 82.4 87.2 75.7 4.4 30.6 52.7 9.8
## 7 7 Bengkulu 18.8 13.6 83.0 72.1 78.1 4.7 36.3 29.5 8.7
## 8 8 Lampung 15.9 10.7 85.4 84.5 78.5 3.3 25.6 28.5 6.7
## 9 9 Bangka Belitung 20.1 4.55 94.2 83.0 75.1 3.4 23 41.8 8.8
## 10 10 Kepulauan Riau 15 5.37 91.2 94.7 75.2 1.3 23 44.2 4
## # ℹ 24 more rows# Data spasial
peta_prov <- st_read("~/projekstunting/data/gadm41_IDN_1.shp")## Reading layer `gadm41_IDN_1' from data source
## `C:\Users\Owner\Documents\projekstunting\data\gadm41_IDN_1.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84# Standarisasi nama provinsi
data <- data %>%
mutate(provinsi = toupper(trimws(provinsi)))
peta_prov <- peta_prov %>%
mutate(NAME_1 = toupper(trimws(NAME_1)))
# Join data ke peta
peta_sf <- peta_prov %>%
left_join(data, by = c("NAME_1" = "provinsi"))
# Pastikan tidak ada NA pada Y
summary(peta_sf$Y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.70 17.32 22.05 21.92 24.85 37.00### 3. STATISTIKA DESKRIPTIF & PETA SEBARAN ###
summary(data)## No provinsi Y X1
## Min. : 1.00 Length:34 Min. : 8.70 Min. : 4.000
## 1st Qu.: 9.25 Class :character 1st Qu.:17.32 1st Qu.: 6.058
## Median :17.50 Mode :character Median :22.05 Median : 8.025
## Mean :17.53 Mean :21.92 Mean : 9.566
## 3rd Qu.:25.75 3rd Qu.:24.85 3rd Qu.:11.630
## Max. :35.00 Max. :37.00 Max. :21.660
## X2 X3 X4 X5
## Min. :72.82 Min. :72.10 Min. :41.44 Min. :1.300
## 1st Qu.:81.22 1st Qu.:84.84 1st Qu.:69.46 1st Qu.:2.525
## Median :84.57 Median :89.75 Median :74.60 Median :3.500
## Mean :85.02 Mean :89.42 Mean :72.55 Mean :4.091
## 3rd Qu.:87.78 3rd Qu.:94.85 3rd Qu.:77.48 3rd Qu.:5.250
## Max. :96.83 Max. :99.96 Max. :82.93 Max. :8.500
## X6 X7 X8
## Min. :16.80 Min. : 0.00 Min. : 4.000
## 1st Qu.:23.18 1st Qu.:34.52 1st Qu.: 6.800
## Median :29.25 Median :43.50 Median : 8.600
## Mean :30.24 Mean :44.09 Mean : 8.574
## 3rd Qu.:36.23 3rd Qu.:56.70 3rd Qu.:10.675
## Max. :49.10 Max. :71.00 Max. :13.100sd(peta_sf$Y, na.rm = TRUE)## [1] 5.921766sapply(data[, c("Y", "X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")], sd, na.rm = TRUE)## Y X1 X2 X3 X4 X5 X6 X7
## 5.921766 4.566111 6.021335 6.759146 7.482184 1.976279 8.890960 16.000708
## X8
## 2.659381### HISTOGRAM
hist(
peta_sf$Y,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Prevalensi Stunting (Y)",
xlab = "Prevalensi Stunting (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X1,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Presentase Kemiskinan (X1)",
xlab = "Presentase Kemiskinan (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X2,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Akses Sanitasi Layak (X2)",
xlab = "Presentase Akses Sanitasi Layak (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X3,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Sumber Air Minum Layak (X3)",
xlab = "Presentase Sumber Air Minum Layak (X3) (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X4,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Balita diimunisasi campak (X4)",
xlab = "Presentase Balita diimunisasi campak (X4) (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X5,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Balita mengalami diare (X5)",
xlab = "Presentase Balita mengalami diare 1 bulan terakhir (X5) (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X6,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Balita mengalami ISPA (X6)",
xlab = "Presentase Balita mengalami ISPA 1 bulan terakhir (X6) (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X7,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Penduduk yang memanfaatkan JKN/Jamkesda (X7)",
xlab = "Presentase Penduduk yang memanfaatkan JKN/Jamkesda (X7) (%)",
ylab = "Frekuensi"
)
hist(
peta_sf$X8,
breaks = 8,
col = "lightblue",
border = "white",
main = "Histogram Ibu Balita Berdasarkan Usia Kehamilan Pertama Kali mendapat TTD (Trisemester 1/ minggu 1-12)",
xlab = "Presentase Ibu balita mendapat TTD (X8) (%)",
ylab = "Frekuensi"
)
### BOXPLOT
boxplot(data$Y,
main = "Boxplot Y (Stunting)",
ylab = "Persentase",
col = "orange")
boxplot(data$X1,
main = "Boxplot X1 (Kemiskinan)",
ylab = "Persentase",
col = "orange")
boxplot(data$X2,
main = "Boxplot X2 (Sanitasi Layak)",
ylab = "Persentase",
col = "orange")
boxplot(data$X3,
main = "Boxplot X3 (Air Minum Layak)",
ylab = "Persentase",
col = "orange")
boxplot(data$X4,
main = "Boxplot X4 (Imunisasi Campak)",
ylab = "Persentase",
col = "orange")
boxplot(data$X5,
main = "Boxplot X5 (Diare Balita)",
ylab = "Persentase",
col = "orange")
boxplot(data$X6,
main = "Boxplot X6 (ISPA Balita)",
ylab = "Persentase",
col = "orange")
boxplot(data$X7,
main = "Boxplot X7 (Kepemilikan JKN)",
ylab = "Persentase",
col = "orange")
boxplot(data$X8,
main = "Boxplot X8 (TTD Ibu Hamil)",
ylab = "Persentase",
col = "orange")
# --- 1. Peta Sebaran Variabel Dependen (Y) ---
tm_shape(peta_sf) +
tm_polygons("Y",
palette = "YlOrRd", # Warna kuning ke merah (merah menandakan nilai Y tinggi/buruk)
style = "quantile",
title = "Prevalensi Stunting (%)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Prevalensi Stunting (Y)",
frame = FALSE)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"
## Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.
# --- 2. Peta Sebaran Kemiskinan (X1) ---
tm_shape(peta_sf) +
tm_polygons("X1",
palette = "Oranges", # Oranye untuk Kemiskinan
style = "quantile",
title = "Persentase Kemiskinan (X1)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran Variabel X1",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Oranges" is named
## "brewer.oranges"Multiple palettes called "oranges" found: "brewer.oranges", "matplotlib.oranges". The first one, "brewer.oranges", is returned.
# --- 3. Peta Sebaran Sanitasi Layak (X2) ---
tm_shape(peta_sf) +
tm_polygons("X2",
palette = "Blues", # Biru untuk Hal Positif (Akses Layak)
style = "quantile",
title = "Akses Sanitasi Layak (X2)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran Variabel X2",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.
# --- 4. Peta Sebaran Air Minum Layak (X3) ---
tm_shape(peta_sf) +
tm_polygons("X3",
palette = "BuGn", # Biru-Hijau untuk Air Layak
style = "quantile",
title = "Akses Air Minum Layak (X3)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran Variabel X3",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuGn" is named
## "brewer.bu_gn"Multiple palettes called "bu_gn" found: "brewer.bu_gn", "matplotlib.bu_gn". The first one, "brewer.bu_gn", is returned.
# --- 5. Peta Sebaran Imunisasi Campak (X4) ---
tm_shape(peta_sf) +
tm_polygons("X4",
palette = "Greens", # Hijau untuk Cakupan Kesehatan
style = "quantile",
title = "Cakupan Imunisasi Campak (X4)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran Variabel X4",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
## "brewer.greens"Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.
# --- 5. Peta Sebaran Balita Mengalami Diare dalam 1 Bulan Terakhir (X5) ---
tm_shape(peta_sf) +
tm_polygons("X5",
palette = "Reds",
style = "quantile",
title = "Prevalensi Diare Balita (%)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran Diare (X5)",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.
# --- 6. Peta Sebaran Balita Mengalami ISPA dalam 1 Bulan Terakhir (X6) ---
tm_shape(peta_sf) +
tm_polygons("X6",
palette = "PuRd",
style = "quantile",
title = "Prevalensi ISPA Balita (%)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran ISPA (X6)",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "PuRd" is named
## "brewer.pu_rd"Multiple palettes called "pu_rd" found: "brewer.pu_rd", "matplotlib.pu_rd". The first one, "brewer.pu_rd", is returned.
# --- 6. Peta Sebaran Proporsi penduduk yang memiliki dan memanfaatkan JKN/Jamkesda (X7) ---
tm_shape(peta_sf) +
tm_polygons("X7",
palette = "Blues",
style = "quantile",
title = "Kepemilikan JKN/Jamkesda (%)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran JKN (X7)",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.
# --- 6. Peta Sebaran Ibu Balita Berdasarkan Usia Kehamilan Pertama Kali Mendapat/Diberi/Membeli TTD Menurut Provinsi (Trisemester 1/ minggu 1-12) (X8) ---
tm_shape(peta_sf) +
tm_polygons("X8",
palette = "Greens",
style = "quantile",
title = "Cakupan TTD Ibu Hamil (%)") +
tm_borders(col = "black", lwd = 0.5) +
tm_layout(main.title = "Peta Sebaran TTD (X8)",
frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
## "brewer.greens"Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.
### 4. UJI AUTOKORELASI SPASIAL (GLOBAL & LOKAL) ###
# Konversi ke Spatial #
peta_sp <- as(peta_sf, "Spatial")
# Uji sensitivitas Bobot Spasial #
library(spdep)
library(spatialreg)
# Queen
nb_queen <- poly2nb(peta_sp, queen = TRUE)## Warning in poly2nb(peta_sp, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(peta_sp, queen = TRUE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.listw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
# Rook
nb_rook <- poly2nb(peta_sp, queen = FALSE)## Warning in poly2nb(peta_sp, queen = FALSE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(peta_sp, queen = FALSE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.listw_rook <- nb2listw(nb_rook, style = "W", zero.policy = TRUE)
# Koordinat
library(sf)
# Konversi sp -> sf
peta_sf <- st_as_sf(peta_sp)
# Ambil centroid
coords <- st_coordinates(st_centroid(peta_sf))## Warning: st_centroid assumes attributes are constant over geometries# KNN-4
knn4 <- knearneigh(coords, k = 4)
nb_knn4 <- knn2nb(knn4)
listw_knn4 <- nb2listw(nb_knn4, style = "W", zero.policy = TRUE)
# SEM - Queen
sem_queen <- errorsarlm(
Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp,
listw = listw_queen,
zero.policy = TRUE
)
# SEM - Rook
sem_rook <- errorsarlm(
Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp,
listw = listw_rook,
zero.policy = TRUE
)
# SEM - KNN-4
sem_knn4 <- errorsarlm(
Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp,
listw = listw_knn4,
zero.policy = TRUE
)
sens_bobot <- data.frame(
Bobot = c("Queen", "Rook", "KNN-4"),
Rho = c(
sem_queen$lambda,
sem_rook$lambda,
sem_knn4$lambda
),
AIC = c(
AIC(sem_queen),
AIC(sem_rook),
AIC(sem_knn4)
),
LogLik = c(
as.numeric(logLik(sem_queen)),
as.numeric(logLik(sem_rook)),
as.numeric(logLik(sem_knn4))
)
)
nb <- poly2nb(peta_sp)## Warning in poly2nb(peta_sp): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(peta_sp): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.listw <- nb2listw(nb, style = "W", zero.policy = TRUE)
sens_bobot## Bobot Rho AIC LogLik
## 1 Queen 0.6147611 201.6806 -89.84032
## 2 Rook 0.6147611 201.6806 -89.84032
## 3 KNN-4 0.6884452 204.6052 -91.30260# Global Moran’s I
vars <- c("Y","X1","X2","X3","X4","X5","X6","X7","X8")
moran_perm <- lapply(vars, function(v){
mc <- moran.mc(peta_sp[[v]],
listw,
nsim = 999,
zero.policy = TRUE)
data.frame(
Variabel = v,
Moran_I = mc$statistic,
Mean_MC = mean(mc$res),
SD_MC = sd(mc$res),
P_value = mc$p.value
)
})
moran_perm_table <- do.call(rbind, moran_perm)
moran_perm_table## Variabel Moran_I Mean_MC SD_MC P_value
## statistic Y 0.30114919 -0.02522324 0.1551077 0.017
## statistic1 X1 0.39123982 -0.02124880 0.1613032 0.011
## statistic2 X2 -0.02745353 -0.02225688 0.1638811 0.502
## statistic3 X3 0.43481676 -0.02545351 0.1571400 0.001
## statistic4 X4 0.17705169 -0.03269849 0.1381662 0.075
## statistic5 X5 0.60238422 -0.01152743 0.1595558 0.001
## statistic6 X6 0.68829155 -0.02235610 0.1626601 0.001
## statistic7 X7 0.09693227 -0.02659225 0.1601075 0.228
## statistic8 X8 0.61683899 -0.02833169 0.1603854 0.001moran_perm_table <- moran_perm_table %>%
mutate(
Keputusan = ifelse(P_value < 0.05,
"Ada autokorelasi spasial",
"Tidak ada autokorelasi spasial")
)
moran_perm_table## Variabel Moran_I Mean_MC SD_MC P_value
## statistic Y 0.30114919 -0.02522324 0.1551077 0.017
## statistic1 X1 0.39123982 -0.02124880 0.1613032 0.011
## statistic2 X2 -0.02745353 -0.02225688 0.1638811 0.502
## statistic3 X3 0.43481676 -0.02545351 0.1571400 0.001
## statistic4 X4 0.17705169 -0.03269849 0.1381662 0.075
## statistic5 X5 0.60238422 -0.01152743 0.1595558 0.001
## statistic6 X6 0.68829155 -0.02235610 0.1626601 0.001
## statistic7 X7 0.09693227 -0.02659225 0.1601075 0.228
## statistic8 X8 0.61683899 -0.02833169 0.1603854 0.001
## Keputusan
## statistic Ada autokorelasi spasial
## statistic1 Ada autokorelasi spasial
## statistic2 Tidak ada autokorelasi spasial
## statistic3 Ada autokorelasi spasial
## statistic4 Tidak ada autokorelasi spasial
## statistic5 Ada autokorelasi spasial
## statistic6 Ada autokorelasi spasial
## statistic7 Tidak ada autokorelasi spasial
## statistic8 Ada autokorelasi spasial# Geary's C
vars <- c("Y","X1","X2","X3","X4","X5","X6","X7","X8")
geary_mc_results <- lapply(vars, function(v){
gc <- geary.mc(peta_sp[[v]],
listw,
nsim = 999,
zero.policy = TRUE)
data.frame(
Variabel = v,
Geary_C = gc$statistic,
P_value_MC = gc$p.value
)
})
geary_mc_table <- do.call(rbind, geary_mc_results)
geary_mc_table## Variabel Geary_C P_value_MC
## statistic Y 0.2200024 0.001
## statistic1 X1 0.3440318 0.001
## statistic2 X2 0.7764312 0.300
## statistic3 X3 0.4327746 0.009
## statistic4 X4 0.5313374 0.123
## statistic5 X5 0.2274290 0.001
## statistic6 X6 0.2467170 0.001
## statistic7 X7 0.7858651 0.338
## statistic8 X8 0.2586480 0.001# Getis Ord
for(v in vars){
gi <- localG(peta_sp@data[[v]], listw, zero.policy = TRUE)
peta_sf[[paste0("GiZ_",v)]] <- as.numeric(gi)
peta_sf[[paste0("GiCat_",v)]] <- "Not Significant"
peta_sf[[paste0("GiCat_",v)]][gi > 1.96] <- "Hot Spot"
peta_sf[[paste0("GiCat_",v)]][gi < -1.96] <- "Cold Spot"
}
tm_shape(peta_sf) +
tm_polygons("GiCat_Y",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (Y)") +
tm_borders() +
tm_layout(frame = FALSE)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X1",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X1)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X2",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X2)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'The visual variable `fill` of the layer "polygons" contains a unique value.
## Therefore a categorical scale is applied (tm_scale_categorical).
tm_shape(peta_sf) +
tm_polygons("GiCat_X3",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X3)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X4",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X4)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X5",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X5)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X6",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X6)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X7",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X7)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("GiCat_X8",
palette = c("red","blue","grey80"),
title = "Getis-Ord Gi* (X8)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
#LISA
for(v in vars){
lisa <- localmoran(peta_sp@data[[v]],
listw,
zero.policy = TRUE)
peta_sf[[paste0("Ii_",v)]] <- lisa[,1]
peta_sf[[paste0("P_",v)]] <- lisa[,5]
# Spatial lag
lagv <- lag.listw(listw, peta_sp@data[[v]], zero.policy = TRUE)
# Klasifikasi LISA
peta_sf[[paste0("LISA_",v)]] <- "Not Significant"
sig <- peta_sf[[paste0("P_",v)]] < 0.05
peta_sf[[paste0("LISA_",v)]][ sig & peta_sp@data[[v]] > mean(peta_sp@data[[v]]) & lagv > mean(lagv) ] <- "High-High"
peta_sf[[paste0("LISA_",v)]][ sig & peta_sp@data[[v]] < mean(peta_sp@data[[v]]) & lagv < mean(lagv) ] <- "Low-Low"
peta_sf[[paste0("LISA_",v)]][ sig & peta_sp@data[[v]] > mean(peta_sp@data[[v]]) & lagv < mean(lagv) ] <- "High-Low"
peta_sf[[paste0("LISA_",v)]][ sig & peta_sp@data[[v]] < mean(peta_sp@data[[v]]) & lagv > mean(lagv) ] <- "Low-High"
}
tm_shape(peta_sf) +
tm_polygons("LISA_Y",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Stunting (Y)") +
tm_borders() +
tm_layout(frame = FALSE)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X1",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Kemsikinan (X1)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X2",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Sanitasi Layak (X2)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'The visual variable `fill` of the layer "polygons" contains a unique value.
## Therefore a categorical scale is applied (tm_scale_categorical).
tm_shape(peta_sf) +
tm_polygons("LISA_X3",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Sumber Air Minum Layak (X3)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X4",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Imunisasi Campak (X4)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X5",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Balita Diare (X5)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X6",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA ISPA (X6)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X7",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA Pemanfaatkan Jamkesda (X7)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
tm_shape(peta_sf) +
tm_polygons("LISA_X8",
palette = c("red","blue","orange","lightblue","grey80"),
title = "LISA TTD pada Ibu Hamil (X8)") +
tm_borders() +
tm_layout(frame = FALSE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
### MODEL OLS DAN UJI ASUMSI
model_ols <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp)
summary(model_ols)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data = peta_sp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.9783 -3.0011 -0.1851 2.7768 9.8035
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.19596 16.88791 3.032 0.0056 **
## X1 0.44128 0.20439 2.159 0.0407 *
## X2 -0.03630 0.19359 -0.187 0.8528
## X3 -0.16289 0.15643 -1.041 0.3077
## X4 -0.24448 0.14071 -1.737 0.0946 .
## X5 0.92409 0.91384 1.011 0.3216
## X6 -0.33459 0.24374 -1.373 0.1820
## X7 0.13549 0.06526 2.076 0.0483 *
## X8 0.26343 0.47314 0.557 0.5826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.615 on 25 degrees of freedom
## Multiple R-squared: 0.5399, Adjusted R-squared: 0.3927
## F-statistic: 3.667 on 8 and 25 DF, p-value: 0.005892shapiro.test(residuals(model_ols))##
## Shapiro-Wilk normality test
##
## data: residuals(model_ols)
## W = 0.98697, p-value = 0.9495library(car)## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(model_ols)## X1 X2 X3 X4 X5 X6 X7 X8
## 1.349650 2.105410 1.732166 1.717528 5.053915 7.276874 1.689753 2.453203moran.test(residuals(model_ols), listw, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: residuals(model_ols)
## weights: listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 2.1563, p-value = 0.01553
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.35409264 -0.03448276 0.03247513library(lmtest)## Warning: package 'lmtest' was built under R version 4.4.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.4.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(model_ols)##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 6.5636, df = 8, p-value = 0.5844resettest(model_ols, power = 2:3, type = "fitted")##
## RESET test
##
## data: model_ols
## RESET = 1.6085, df1 = 2, df2 = 23, p-value = 0.2219### UJI LAGRANGE MULTIPIER
lm_tests <- lm.LMtests(model_ols, listw,
test = c("LMerr","LMlag","RLMerr","RLMlag"))## Please update scripts to use lm.RStests in place of lm.LMtestsprint(lm_tests)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data =
## peta_sp)
## test weights: listw
##
## RSerr = 4.5815, df = 1, p-value = 0.03232
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data =
## peta_sp)
## test weights: listw
##
## RSlag = 0.034662, df = 1, p-value = 0.8523
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data =
## peta_sp)
## test weights: listw
##
## adjRSerr = 6.5598, df = 1, p-value = 0.01043
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data =
## peta_sp)
## test weights: listw
##
## adjRSlag = 2.0129, df = 1, p-value = 0.156### MODEL SEM
model_sem <- errorsarlm(
Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp,
listw = listw,
zero.policy = TRUE
)
summary(model_sem)##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
## data = peta_sp, listw = listw, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.79176 -1.97548 -0.48389 2.15162 7.68875
##
## Type: error
## Regions with no neighbours included:
## 2 3 21 22
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 53.500420 14.178007 3.7735 0.000161
## X1 0.647483 0.152695 4.2404 2.231e-05
## X2 -0.054049 0.120969 -0.4468 0.655016
## X3 -0.240040 0.120001 -2.0003 0.045467
## X4 -0.057952 0.099201 -0.5842 0.559093
## X5 1.459686 0.653286 2.2344 0.025458
## X6 -0.423315 0.161225 -2.6256 0.008649
## X7 0.091548 0.041934 2.1831 0.029027
## X8 -0.464273 0.416227 -1.1154 0.264665
##
## Lambda: 0.61476, LR test value: 10.344, p-value: 0.0012986
## Asymptotic standard error: 0.10335
## z-value: 5.9483, p-value: 2.7091e-09
## Wald statistic: 35.383, p-value: 2.7091e-09
##
## Log likelihood: -89.84032 for error model
## ML residual variance (sigma squared): 9.4459, (sigma: 3.0734)
## Number of observations: 34
## Number of parameters estimated: 11
## AIC: 201.68, (AIC for lm: 210.03)### MODEL SDEM
model_sdem <- lagsarlm(
Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
data = peta_sp,
listw = listw,
type = "mixed",
zero.policy = TRUE
)
summary(model_sdem)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8,
## data = peta_sp, listw = listw, type = "mixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.09751 -1.69239 -0.16083 2.43934 5.58129
##
## Type: mixed
## Regions with no neighbours included:
## 2 3 21 22
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 70.6776301 15.3623566 4.6007 4.211e-06
## X1 0.5247868 0.1553907 3.3772 0.0007323
## X2 -0.1548318 0.1294267 -1.1963 0.2315837
## X3 -0.2701861 0.1091946 -2.4744 0.0133478
## X4 -0.1640224 0.1039014 -1.5786 0.1144198
## X5 0.9739281 0.6711168 1.4512 0.1467228
## X6 -0.2588096 0.1625377 -1.5923 0.1113163
## X7 0.0915262 0.0454409 2.0142 0.0439905
## X8 -0.2375314 0.4550365 -0.5220 0.6016668
## lag.X1 -0.5807532 0.1977490 -2.9368 0.0033160
## lag.X2 -0.1835913 0.2143989 -0.8563 0.3918278
## lag.X3 0.1611125 0.1567840 1.0276 0.3041342
## lag.X4 -0.3482804 0.1932158 -1.8025 0.0714596
## lag.X5 -0.0444022 1.3441293 -0.0330 0.9736473
## lag.X6 0.2258289 0.3379789 0.6682 0.5040221
## lag.X7 0.0079314 0.0620362 0.1279 0.8982667
## lag.X8 1.5680448 0.5931986 2.6434 0.0082085
##
## Rho: 0.36314, LR test value: 3.6136, p-value: 0.05731
## Asymptotic standard error: 0.13831
## z-value: 2.6255, p-value: 0.0086512
## Wald statistic: 6.8934, p-value: 0.0086512
##
## Log likelihood: -84.24684 for mixed model
## ML residual variance (sigma squared): 7.8285, (sigma: 2.7979)
## Number of observations: 34
## Number of parameters estimated: 19
## AIC: 206.49, (AIC for lm: 208.11)
## LM test for residual autocorrelation
## test value: 1.3553, p-value: 0.24435### Uji Asumsi Model
# Residual
res_sem <- residuals(model_sem)
res_sdem <- residuals(model_sdem)
library(spdep)
# SEM
moran.test(res_sem, listw, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: res_sem
## weights: listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 0.48278, p-value = 0.3146
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.05265062 -0.03448276 0.03257430# SDEM
moran.test(res_sdem, listw, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: res_sdem
## weights: listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.23414, p-value = 0.5926
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.07712971 -0.03448276 0.03317696### Normalitas
# SEM
shapiro.test(res_sem)##
## Shapiro-Wilk normality test
##
## data: res_sem
## W = 0.97502, p-value = 0.612# SDEM
shapiro.test(res_sdem)##
## Shapiro-Wilk normality test
##
## data: res_sdem
## W = 0.97735, p-value = 0.6874### Homoskedastisitas
library(lmtest)
# SEM
bptest(res_sem ~ fitted(model_sem))## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page##
## studentized Breusch-Pagan test
##
## data: res_sem ~ fitted(model_sem)
## BP = 1.7004e-05, df = 1, p-value = 0.9967# SDEM
bptest(res_sdem ~ fitted(model_sdem))## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page##
## studentized Breusch-Pagan test
##
## data: res_sdem ~ fitted(model_sdem)
## BP = 1.8049, df = 1, p-value = 0.1791### PEMILIHAN MODEL TERBAIK
# Bandingkan AIC dan Log-Likelihood
data.frame(
Model = c("OLS", "SEM", "SDEM"),
df = c(
attr(logLik(model_ols), "df"),
attr(logLik(model_sem), "df"),
attr(logLik(model_sdem), "df")
),
LogLik = c(
as.numeric(logLik(model_ols)),
as.numeric(logLik(model_sem)),
as.numeric(logLik(model_sdem))
),
AIC = c(
AIC(model_ols),
AIC(model_sem),
AIC(model_sdem)
)
)## Model df LogLik AIC
## 1 OLS 10 -95.01255 210.0251
## 2 SEM 11 -89.84032 201.6806
## 3 SDEM 19 -84.24684 206.4937# Prediksi dari SEM dan SDEM
# Prediksi
peta_sf$Y_pred_SEM <- as.numeric(predict(model_sem))## This method assumes the response is known - see manual pagepeta_sf$Y_pred_SDEM <- as.numeric(predict(model_sdem))## This method assumes the response is known - see manual page# Residual
peta_sf$res_SEM <- residuals(model_sem)
peta_sf$res_SDEM <- residuals(model_sdem)
# Residual standar (z-score)
peta_sf$z_res_SEM <- scale(peta_sf$res_SEM)
peta_sf$z_res_SDEM <- scale(peta_sf$res_SDEM)
### PETA PREDIKSI SEM
peta_pred_sem <-
tm_shape(peta_sf) +
tm_polygons(
"Y_pred_SEM",
palette = "plasma",
style = "quantile",
title = "Prediksi Y (%)"
) +
tm_borders(col = "white", lwd = 0.7) +
tm_layout(
main.title = "Peta Prediksi SEM",
frame = FALSE
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`peta_pred_sem
### PETA RESIDUAL STANDAR SEM
peta_res_sem <-
tm_shape(peta_sf) +
tm_polygons(
"z_res_SEM",
palette = "-RdBu",
style = "fixed",
breaks = c(-3,-2,-1,0,1,2,3),
title = "Residual (z)"
) +
tm_borders(col = "white", lwd = 0.7) +
tm_layout(
main.title = "Peta Residual Standar SEM",
frame = FALSE
)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`peta_res_sem## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
## "rd_bu" (in long format "brewer.rd_bu")
### PETA PREDIKSI SDEM
peta_pred_sdem <-
tm_shape(peta_sf) +
tm_polygons(
"Y_pred_SDEM",
palette = "plasma",
style = "quantile",
title = "Prediksi Y (%)"
) +
tm_borders(col = "white", lwd = 0.7) +
tm_layout(
main.title = "Peta Prediksi SDEM",
frame = FALSE
)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`peta_pred_sdem
### PETA RESIDUAL STANDAR SDEM
peta_res_sdem <-
tm_shape(peta_sf) +
tm_polygons(
"z_res_SDEM",
palette = "-RdBu",
style = "fixed",
breaks = c(-3,-2,-1,0,1,2,3),
title = "Residual (z)"
) +
tm_borders(col = "white", lwd = 0.7) +
tm_layout(
main.title = "Peta Residual Standar SDEM",
frame = FALSE
)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`peta_res_sdem## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
## "rd_bu" (in long format "brewer.rd_bu")
### SPASIAL INTERPOLATION
# Siapkan data titik (centroid)
# Ambil centroid provinsi
coords <- coordinates(peta_sp)
titik_prov <- SpatialPointsDataFrame(
coords = coords,
data = peta_sp@data,
proj4string = CRS(proj4string(peta_sp))
)
# Buat grid interpolasi
# Batas wilayah
x.range <- bbox(peta_sp)[1,]
y.range <- bbox(peta_sp)[2,]
# Grid (atur resolusi: makin kecil makin halus)
grd <- expand.grid(
x = seq(x.range[1], x.range[2], by = 0.1),
y = seq(y.range[1], y.range[2], by = 0.1)
)
coordinates(grd) <- ~x+y
gridded(grd) <- TRUE
proj4string(grd) <- CRS(proj4string(peta_sp))
# Interpolasi IDW
# Untuk variabel observasi (Y)
idw_Y <- idw(
formula = Y ~ 1,
locations = titik_prov,
newdata = grd,
idp = 2
)## [inverse distance weighted interpolation]# Untuk prediksi SEM
peta_sp$Y_pred_SEM <- predict(model_sem)## This method assumes the response is known - see manual pagetitik_prov$Y_pred_SEM <- peta_sp$Y_pred_SEM
idw_SEM <- idw(
formula = Y_pred_SEM ~ 1,
locations = titik_prov,
newdata = grd,
idp = 2
)## [inverse distance weighted interpolation]# Konversi ke raster & mask peta
r_idw_sem <- raster(idw_SEM)
r_idw_sem <- mask(r_idw_sem, peta_sp)
# Plot peta interpolasi
tm_shape(r_idw_sem) +
tm_raster(
col = "var1.pred",
palette = "-Spectral",
style = "cont",
title = "Interpolasi Prediksi SEM"
) +
tm_shape(peta_sp) +
tm_borders(col = "black", alpha = 0.3) +
tm_layout(frame = FALSE)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
## `tm_scale_continuous()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values') to
## 'tm_scale_continuous(<HERE>)'
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "-Spectral" is
## named "spectral" (in long format "brewer.spectral")