PENDEKATAN SPASIAL UNTUK MENGANALISIS INCIDENCE RATE DBD DI SUMATERA UTARA TAHUN 2024
Ghisela Nisrina Nazifah (140610230032)
Departemen Statistika, FMIPA Universitas PadjadjaranABSTRAK
Demam Berdarah Dengue (DBD) merupakan masalah kesehatan masyarakat yang penyebarannya dipengaruhi oleh kondisi lingkungan dan keterkaitan antarwilayah. Penelitian ini bertujuan untuk mengidentifikasi ketergantungan spasial serta memodelkan sebaran Incidence Rate (IR) DBD antar kabupaten/kota di Provinsi Sumatera Utara. Data yang digunakan meliputi data DBD serta data lingkungan dan sosial ekonomi pada tingkat kabupaten/kota. Metode penelitian mencakup uji autokorelasi spasial menggunakan Moran’s I, pemodelan regresi spasial dengan model Spatial Autoregressive (SAR), serta interpolasi spasial menggunakan metode Inverse Distance Weighting (IDW) dan Ordinary Kriging. Hasil penelitian menunjukkan adanya autokorelasi spasial positif dan signifikan pada DBD IR (Moran’s I = 0,397; p-value = 0,001) serta parameter autoregresif spasial yang signifikan pada model SAR (ρ = 0,591; p-value = 0,006). Variabel kejadian banjir berpengaruh signifikan terhadap DBD IR (p-value = 0,032), sementara variabel lainnya tidak signifikan. Hasil interpolasi menunjukkan bahwa Kriging menghasilkan permukaan yang lebih halus dan terstruktur dibandingkan IDW. Secara keseluruhan, penelitian ini menegaskan pentingnya pendekatan spasial dalam analisis dan pengendalian DBD antarwilayah, dengan menekankan peran kondisi lingkungan dalam upaya mitigasi risiko DBD.
BAB 1. PENDAHULUAN
1.1 Latar Belakang
Demam Berdarah Dengue (DBD) merupakan salah satu penyakit menular yang masih menjadi permasalahan kesehatan masyarakat di Indonesia dan sangat dipengaruhi oleh kondisi lingkungan. Salah satu faktor lingkungan yang berperan penting dalam peningkatan risiko DBD adalah kejadian banjir, karena genangan air yang terbentuk pascabanjir dapat menjadi tempat berkembang biaknya nyamuk Aedes aegypti sebagai vektor utama DBD (WHO, 2012; Kemenkes RI, 2023). Wilayah Sumatera, termasuk Provinsi Sumatera Utara, senantiasa mengalami kejadian banjir setiap tahunnya dan setiap tahun menunjukkan peningkatan frekuensi banjir. Pada tahun 2024, terdapat 1.420 kejadian banjir dan tahun ini (per tanggal 21 Desember) terdapat 1.598 kejadian yang masih terus bertambah (Pusdatinkom BNPB, 2025). Kondisi ini menunjukkan meningkatnya tekanan lingkungan yang berpotensi memperbesar risiko penularan DBD apabila tidak diantisipasi secara tepat.
Peningkatan kejadian banjir tersebut berpotensi menyebabkan perbedaan tingkat kejadian DBD antarwilayah di Provinsi Sumatera Utara. Oleh karena itu, analisis terhadap incidence rate DBD per 100.000 penduduk menjadi penting untuk menggambarkan tingkat kejadian DBD secara proporsional dan dapat dibandingkan antarwilayah. Pemahaman yang baik mengenai variasi tingkat kejadian DBD antarwilayah diharapkan dapat menjadi dasar dalam upaya pencegahan dan pengendalian DBD yang lebih efektif, tepat sasaran, dan komprehensif.
1.2 Identifikasi Masalah
Peta Incidence Rate DBD menunjukkan adanya variasi spasial antar kabupaten/kota di Sumatera Utara, di mana beberapa wilayah memiliki angka kesakitan yang jauh lebih tinggi dibandingkan wilayah lainnya. Hal ini mengindikasikan bahwa distribusi DBD tidak bersifat homogen secara spasial. Selanjutnya secara visual dapat dilihat bahwa beberapa daerah dengan nilai Incidence Rate rendah berdekatan dengan daerah nilai rendah lainnya, begitu juga dengan daerah denga Incidence Rate tinggi. Sehingga dalam hal ini diperlukan analisis spasial lanjutan untuk menangkap pola penyebaran yang lebih akurat.
1.3 Maksud dan Tujuan Penelitian
Maksud dari penelitian ini adalah untuk menerapkan pendekatan spasial dalam menganalisis pola penyebaran Incidence Rate Demam Berdarah Dengue (DBD) di Provinsi Sumatera Utara tahun 2024. Penelitian ini bertujuan untuk mengidentifikasi ketergantungan spasial antarwilayah serta memodelkan sebaran Incidence Rate DBD melalui metode pemodelan interpolasi spasial guna memperoleh gambaran permukaan spasial DBD secara lebih komprehensif, sehingga dapat menggambarkan variasi spasial DBD secara lebih detail antarwilayah kabupaten/kota.
1.4 Batasan Penelitian
Penelitian ini dibatasi pada analisis spasial angka kesakitan Demam Berdarah Dengue (DBD) per 100.000 penduduk di Provinsi Sumatera Utara tahun 2024 pada tingkat kabupaten/kota. Analisis yang dilakukan meliputi pemodelan ekonometrik spasial untuk mengkaji pengaruh faktor-faktor yang berkaitan dengan kejadian DBD, serta pemodelan interpolasi spasial menggunakan metode IDW dan Ordinary Kriging. Variabel curah hujan yang mungkin menjadi penyebab DBD tidak digunakan secara langsung dan digantikan dengan variabel cuaca ekstrem sebagai proksi kondisi iklim, sesuai dengan ketersediaan data.
BAB 2. TINJAUAN PUSTAKA
2.1 Spatial Autocorrellation
Autokorelasi spasial menggambarkan adanya keterkaitan nilai suatu variabel antarwilayah yang berdekatan secara geografis. Keberadaan autokorelasi spasial menunjukkan bahwa data tidak terdistribusi secara acak di ruang, sehingga pelanggaran asumsi independensi pada model regresi klasik dapat terjadi dan diperlukan pendekatan analisis spasial lanjutan (Anselin, 1988).
2.2 Ordinary Least Square
Ordinary Least Squares (OLS) merupakan metode estimasi parameter yang paling umum digunakan dalam analisis regresi linier untuk mengkaji hubungan antara variabel dependen dan variabel independen.
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
Metode ini mengestimasi parameter dengan meminimalkan jumlah kuadrat galat. Namun, dalam konteks data spasial, OLS sering menghadapi permasalahan pelanggaran asumsi klasik, khususnya autokorelasi spasial, yang dapat menyebabkan estimasi menjadi tidak efisien dan inferensi statistik menjadi bias (Anselin, 1988)
2.3 Uji Lagrange Multiplier
Uji Lagrange Multiplier (LM) digunakan untuk mendeteksi keberadaan ketergantungan spasial dalam model regresi OLS. Uji ini membantu menentukan apakah model spasial diperlukan dan jenis model spasial yang paling sesuai, seperti Spatial Autoregressive (SAR) atau Spatial Error Model (SEM). LM test dan robust LM test sering digunakan sebagai dasar pemilihan model spasial lanjutan dalam analisis spasial (Anselin et al., 1996).
2.4 Spatial Autoregressive
Spatial Autoregressive Model (SAR) merupakan model regresi spasial yang memasukkan pengaruh ketergantungan spasial pada variabel dependen melalui lag spasial.
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
Model ini mengasumsikan bahwa nilai suatu wilayah dipengaruhi oleh nilai wilayah sekitarnya. SAR banyak digunakan untuk menangkap efek penyebaran (spillover effects) antarwilayah, terutama ketika terdapat interaksi spasial yang kuat (LeSage & Pace, 2009).
2.5 Spatial Durbin Model
Spatial Durbin Model (SDM) merupakan pengembangan dari model SAR dengan memasukkan lag spasial baik pada variabel dependen maupun variabel independen.
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \]
Model ini mampu menangkap efek langsung dan tidak langsung (spillover) dari variabel independen antarwilayah. SDM dianggap sebagai model yang lebih fleksibel dan komprehensif dalam analisis spasial, khususnya ketika terdapat kemungkinan interaksi kompleks antarwilayah (LeSage & Pace, 2009).
2.6 Metode Evaluasi Model
2.6.1 Akaike Information Criterion (AIC)
Akaike Information Criterion (AIC) merupakan ukuran evaluasi model yang mempertimbangkan kecocokan model terhadap data dan jumlah parameter yang digunakan. Model dengan nilai AIC yang lebih kecil dianggap lebih baik. AIC banyak digunakan dalam pemilihan model ekonometrika, termasuk dalam analisis spasial (Akaike, 1974).
2.6.2 Bayesian Information Criterion (BIC)
Bayesian Information Criterion (BIC) merupakan kriteria evaluasi model yang memberikan penalti lebih besar terhadap kompleksitas model dibandingkan AIC. BIC cenderung memilih model yang lebih sederhana dan sering digunakan sebagai pembanding dalam pemilihan model ekonometrika spasial (Schwarz, 1978).
2.7 Inverse Disctance Weighting
Inverse Distance Weighting (IDW) merupakan metode interpolasi deterministik yang mengestimasi nilai pada lokasi tidak terobservasi berdasarkan nilai titik di sekitarnya dengan bobot yang berbanding terbalik dengan jarak.
\[ \hat{Z}(x_0) = \frac{\displaystyle\sum_{i=1}^{n} Z(x_i)\, d_{i0}^{-p}} {\displaystyle\sum_{i=1}^{n} d_{i0}^{-p}} \]
Metode ini sederhana, mudah diimplementasikan, dan sering digunakan sebagai pendekatan awal (baseline) dalam analisis interpolasi spasial. Namun, IDW tidak mempertimbangkan struktur autokorelasi spasial secara eksplisit (Burrough & McDonnell, 1998).
2.8 Kriging
Kriging merupakan metode interpolasi geostatistik yang memanfaatkan struktur autokorelasi spasial melalui variogram. Berbeda dengan IDW, kriging menghasilkan estimasi yang bersifat optimal secara statistik serta menyediakan informasi mengenai ketidakpastian prediksi.
\[ \hat{Z}(x_0) = \sum_{i=1}^{n} \lambda_i Z(x_i), \quad \text{dengan} \quad \sum_{i=1}^{n} \lambda_i = 1 \]
Ordinary Kriging merupakan jenis kriging yang paling umum digunakan ketika nilai rata-rata tidak diketahui dan diasumsikan konstan secara lokal (Cressie, 1993).
BAB 3. METODOLOGI PENELITIAN
3.1 Data
Penelitian ini menggunakan data sekunder yang diperoleh dari Badan Pusat Statistik (BPS) Provinsi Sumatera Utara dan Badan Nasional Penanggulangan Bencana (BNPB). Unit analisis penelitian adalah kabupaten/kota di Provinsi Sumatera Utara, dengan data bersifat cross section yang merepresentasikan kondisi wilayah pada tahun 2024.
Variabel dependen dalam penelitian ini adalah angka kesakitan Demam Berdarah Dengue (DBD) yang diukur menggunakan incidence rate (IR) per 100.000 penduduk menurut kabupaten/kota, yang bersumber dari publikasi BPS Provinsi Sumatera Utara. Variabel independen yang digunakan mencakup faktor lingkungan dan sosial ekonomi, yaitu kejadian banjir dan kejadian cuaca ekstrem yang diperoleh dari BNPB, serta kepadatan penduduk, persentase penduduk miskin, dan persentase rumah tangga dengan akses sanitasi layak yang bersumber dari BPS Provinsi Sumatera Utara dan publikasi Provinsi Sumatera Utara dalam Angka 2025.
Selain data atribut, penelitian ini juga menggunakan data spasial berupa peta batas administrasi kabupaten/kota di Provinsi Sumatera Utara dalam format shapefile (.shp) yang digunakan untuk membentuk hubungan antarwilayah dalam analisis spasial dan pemodelan interpolasi.
| Variabel | Simbol | Definisi.Operasional | Satuan |
|---|---|---|---|
| Angka Kesakitan DBD | Y | Jumlah kasus Demam Berdarah Dengue per 100.000 penduduk menurut kabupaten/kota. | Kasus per 100.000 penduduk |
| Banjir | X₁ | Jumlah kejadian banjir yang terjadi dalam satu tahun menurut kabupaten/kota. | Kejadian |
| Cuaca Ekstrem | X₂ | Jumlah kejadian cuaca ekstrem yang terjadi dalam satu tahun menurut kabupaten/kota. | Kejadian |
| Kepadatan Penduduk | X₃ | Jumlah penduduk per satuan luas wilayah menurut kabupaten/kota. | Jiwa/km² |
| Persentase Penduduk Miskin | X₄ | Persentase penduduk yang berada di bawah garis kemiskinan menurut kabupaten/kota. | % |
| Akses Sanitasi Layak | X₅ | Persentase rumah tangga yang memiliki akses terhadap fasilitas sanitasi layak. | % |
3.2 Metode Analisis
Analisis dalam penelitian ini dilakukan menggunakan pendekatan statistik spasial dengan bantuan perangkat lunak R. Tahapan analisis diawali dengan analisis autokorelasi spasial menggunakan Moran’s I, Geary’s C, Local Indicators of Spatial Association (LISA), dan Getis-Ord Gi* untuk mengidentifikasi pola keterkaitan spasial angka kesakitan Demam Berdarah Dengue (DBD) antar kabupaten/kota di Provinsi Sumatera Utara.
Selanjutnya, dilakukan pemodelan regresi Ordinary Least Squares (OLS) sebagai model dasar yang diikuti dengan uji Lagrange Multiplier (LM) untuk menentukan model ekonometrika spasial yang sesuai. Berdasarkan hasil uji tersebut, analisis dilanjutkan dengan penerapan model ekonometrika spasial dan pemilihan model terbaik dilakukan berdasarkan kriteria informasi dan signifikansi parameter spasial.
Selain itu, penelitian ini juga menerapkan pemodelan interpolasi spasial menggunakan metode Inverse Distance Weighting (IDW) dan Ordinary Kriging untuk menggambarkan sebaran angka kesakitan DBD secara kontinu di wilayah penelitian.
3.3 Alur Kerja Penelitian
Alur kerja pada penelitian ini disususn dalam Flowchart sebagai berikut:
BAB 4. HASIL DAN PEMBAHASAN
4.1 Analisis Statistik Deskriptif
Analisis deskriptif dilakukan untuk memberikan gambaran umum mengenai karakteristik data penelitian. Hasil analisis deskriptif data dapat dilihat pada tabel berikut:
| Variabel | Mean | Minimum | Maximum | Simpangan Baku |
|---|---|---|---|---|
| Banjir | 3.64 | 0.00 | 14.00 | 3.30 |
| Cuaca Ekstrem | 1.30 | 0.00 | 15.00 | 2.72 |
| Incidence Rate DBD | 82.91 | 2.30 | 276.90 | 76.49 |
| Kepadatan Penduduk | 1159.93 | 41.16 | 8902.16 | 2221.83 |
| Persentase Penduduk Miskin | 9.73 | 3.44 | 22.68 | 4.22 |
| Akses Sanitasi Layak | 76.69 | 15.33 | 97.92 | 24.12 |
Selanjutnya, untuk memberikan gambaran awal terkait penyebaran masing-masing variabel, dilakukan pemetaan dengan peta choropleth sebagai berikut:
4.2 Ordinary Least Square
Sebagai baseline dalam pemodelan, dilakukan pemodelan dengan OLS, didapatkan persamaan model regresi linear berikut:
\[ Y = 100.11 - 10.15X_1 + 5.98X_2 + 0.001X_3 + 2.60X_4 - 0.19X_5 \]
Selanjutnya, akan dilakukan uji asumsi klasik untuk model tersebut demi menghasilkan model yang hasilnya valid dan reliable. Hasil uji asumsi klasik dapat dilihat pada tabel berikut.
| Jenis Uji | Statistik | Statistik Hitung | p-value | Keterangan |
|---|---|---|---|---|
| Linearitas (Ramsey RESET) | RESET | 1.40906 | 0.26312 | Variabel independen terhadap variabel dependen membentuk model linier |
| Normalitas (Shapiro Wilk) | W | 0.92365 | 0.02316 | Residual tidak berdistribusi normal |
| Homoskedastisitas (Breusch–Pagan) | BP | 2.99327 | 0.70102 | Tidak terdapat gejala heteroskedastisitas |
Meskipun uji Shapiro–Wilk menunjukkan residual tidak berdistribusi normal, transformasi variabel tidak dilakukan karena normalitas residual bukan asumsi utama dalam estimasi OLS, terutama pada ukuran sampel yang relatif besar (>30). Selain itu, hasil uji linearitas dan homoskedastisitas menunjukkan bahwa model telah memenuhi asumsi dasar regresi linear.
| X1 | X2 | X3 | X4 | X5 | Keterangan |
|---|---|---|---|---|---|
| 1.7657 | 1.7275 | 1.0202 | 1.7886 | 1.6084 | Tidak terdapat multikolinearitas |
Berdasarkan hasil uji multikolinearitas data dengan VIF juga diketahui bahwa data terbebas dari multikolinearitas. Selanjutnya dilakukan uji Moran’s I residual untuk melihat apakah residual OLS bersifat acak atau mengalami dependensi spasial:
| I | p(Monte Carlo) | Keterangan |
|---|---|---|
| 0.2301 | 0.0330 | Masih terdapat autokorelasi spasial residual |
Dapat dilihat bahwa permasalahan utama yang teridentifikasi justru adanya autokorelasi spasial residual berdasarkan uji Moran’s I, sehingga pendekatan yang lebih tepat adalah penggunaan model regresi spasial.
4.3 Spatial Autocorrellation
Untuk mengetahui ada atau tidaknya dependensi spasial pada setiap variabel penelitian, dilakukan uji autokorelasi spasial menggunakan Moran’s I, Geary’s C, LISA, dan Getis-Ord Gi*.
4.3.1 Moran’s I
| Variabel | Moran.s.I | Z | p.value | Keterangan |
|---|---|---|---|---|
| DBD_IR | 0.397 | 3.192 | 0.001 | Signifikan |
| Banjir | 0.459 | 3.723 | 0.000 | Signifikan |
| Cuaca_Ekstreme | 0.204 | 2.836 | 0.002 | Signifikan |
| Kepadatan_Penduduk | 0.005 | 0.298 | 0.383 | Tidak Signifikan |
| Persentase_Miskin | 0.763 | 6.161 | 0.000 | Signifikan |
| Persentase_Sanitasi | 0.778 | 5.980 | 0.000 | Signifikan |
Hasil uji Moran’s I menunjukkan bahwa variabel DBD_IR, Banjir, Cuaca_Ekstreme, Persentase_Miskin, dan Persentase_Sanitasi memiliki nilai Moran’s I positif dan signifikan (p-value < 0,05), yang mengindikasikan adanya autokorelasi spasial berupa pola pengelompokan antarwilayah, sedangkan variabel Kepadatan_Penduduk tidak menunjukkan autokorelasi spasial yang signifikan (p-value = 0,383) sehingga pola sebarannya cenderung acak; temuan ini menegaskan bahwa sebagian besar variabel dalam penelitian memiliki dependensi spasial.
4.3.2 Geary’s C
| Variabel | Geary.s.C | Z | p.value | Keterangan |
|---|---|---|---|---|
| DBD_IR | 0.577 | 2.639 | 0.004 | Signifikan |
| Banjir | 0.732 | 1.566 | 0.059 | Tidak Signifikan |
| Cuaca_Ekstreme | 1.368 | -1.255 | 0.895 | Tidak Signifikan |
| Kepadatan_Penduduk | 0.649 | 1.713 | 0.043 | Signifikan |
| Persentase_Miskin | 0.278 | 3.951 | 0.000 | Signifikan |
| Persentase_Sanitasi | 0.170 | 5.332 | 0.000 | Signifikan |
Hasil uji Geary’s C menunjukkan bahwa variabel DBD_IR, Kepadatan_Penduduk, Persentase_Miskin, dan Persentase_Sanitasi memiliki autokorelasi spasial yang signifikan (p-value < 0,05), sedangkan variabel Banjir dan Cuaca_Ekstreme tidak menunjukkan autokorelasi spasial yang signifikan. Perbedaan hasil antara Geary’s C dan Moran’s I terjadi karena Geary’s C lebih sensitif terhadap perbedaan nilai antarwilayah yang bertetangga (lokal), sementara Moran’s I lebih menekankan pola pengelompokan secara global, sehingga variabel tertentu dapat terdeteksi signifikan pada Moran’s I namun tidak pada Geary’s C.
4.3.3 Local Moran’s I (LISA)
Berikut adalah hasil analisis Local Moran’s I yang ditampilkan menggunakan peta untuk setiap variabel.
Peta Local Moran’s I (LISA)
4.3.4 Getis-Ord Gi*
Berikut adalah hasil analisis Getis-Ord Gi* yang ditampilkan menggunakan peta untuk setiap variabel.
Peta Getis-Ord Gi*
4.4 Estimasi Model Spatial Econometrics
Adanya indikasi autokorelasi spasial pada data menunjukkan bahwa model OLS belum sepenuhnya mampu menangkap struktur spasial yang ada. Oleh karena itu, uji Lagrange Multiplier dilakukan untuk menentukan bentuk ketergantungan spasial yang tepat sebelum dilakukan estimasi model regresi spasial.
4.4.1 Uji Lagrange Multiplier
| Statistik.Uji | Value | p.value | Keputusan | |
|---|---|---|---|---|
| RSerr | LM Error (RSer) | 2.5720 | 0.1088 | Tidak Signifikan |
| RSlag | LM Lag (RSlag) | 4.5570 | 0.0328 | Signifikan |
| adjRSerr | Robust LM Error (AdjRSer) | 1.5892 | 0.2074 | Tidak Signifikan |
| adjRSlag | Robust LM Lag (AdjRSlag) | 3.5742 | 0.0587 | Tidak Signifikan |
| SARMA | LM SARMA (Gabungan) | 6.1463 | 0.0463 | Signifikan |
Hasil uji Lagrange Multiplier menunjukkan bahwa LM Lag (RSlag) signifikan (p-value = 0,0328), sedangkan LM Error (RSerr) tidak signifikan, yang mengindikasikan bahwa ketergantungan spasial lebih dominan terjadi pada lag variabel dependen dibandingkan pada error. Meskipun uji Robust LM Lag tidak signifikan, hasil uji LM SARMA yang signifikan (p-value = 0,0463) menguatkan adanya ketergantungan spasial secara umum. Berdasarkan temuan ini, model Spatial Autoregressive (SAR) dipilih sebagai model utama untuk tahap estimasi selanjutnya. Sebagai pembanding, estimasi Spatial Durbin Model (SDM) juga dilakukan untuk mengakomodasi kemungkinan adanya pengaruh spasial tidak hanya pada variabel dependen, tetapi juga pada variabel independen, sehingga memungkinkan evaluasi yang lebih komprehensif terhadap mekanisme ketergantungan spasial dalam data.
4.4.2 Estimasi Model
Estimasi model regresi spasial dalam penelitian ini dilakukan
menggunakan metode maximum likelihood dengan fungsi
lagsarlm. Model Spatial Autoregressive (SAR) diestimasi
untuk menangkap pengaruh ketergantungan spasial pada variabel dependen
melalui parameter autoregresif spasial (ρ), yang merepresentasikan
pengaruh nilai variabel dependen di wilayah sekitar terhadap wilayah
yang diamati. Selain itu, Spatial Durbin Model (SDM) juga diestimasi
dengan mengaktifkan komponen spatial lag pada variabel independen,
sehingga memungkinkan adanya pengaruh tidak langsung dari variabel
penjelas di wilayah tetangga. Estimasi kedua model dilakukan menggunakan
matriks bobot spasial dengan tujuan memperoleh parameter model yang
mencerminkan struktur spasial pada data penelitian.
4.5 Pemilihan Model Terbaik
| Model | AIC | BIC | p uji ρ | p X1 | p X2 | p X3 | p X4 | p X5 |
|---|---|---|---|---|---|---|---|---|
| SAR | 380.527 | 392.499 | 1e-05*** | 0.0293* | 0.461 | 0.298 | 0.398 | 0.462 |
| SDM | 380.331 | 399.785 | 0.00877** | 0.53 | 0.951 | 0.363 | 0.0426* | 0.889 |
Berdasarkan nilai AIC, model SDM menunjukkan kinerja sedikit lebih baik dibandingkan model SAR. Namun, berdasarkan nilai BIC yang memberikan penalti lebih besar terhadap kompleksitas model, model SAR dipilih sebagai model terbaik karena lebih parsimonious. Perbedaan ini menunjukkan adanya trade-off antara goodness-of-fit dan kompleksitas model.
4.5.1 Estimasi Model Terbaik
Berdasarkan hasil perbandingan model sebelumnya, model Spatial Autoregressive (SAR) dipilih menjadi model terbaik. Dengan menggunakan Maximum Likelihood Estimation (MLE), didapatkan model SAR untuk incidence rate DBD di Sumatera Utara sebagai berikut:
\[ \hat{Y}_i = 118.27 - 8.94^{*} X_{1i} + 3.61 X_{2i} + 0.0048 X_{3i} - 2.72 X_{4i} - 0.39 X_{5i} + 0.59^{***} \sum_j w_{ij} Y_j \]
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 7.5003, df = 1, p-value = 0.006169
## sample estimates:
## Log likelihood of sar Log likelihood of ols
## -182.2634 -186.0136| rho | LR.test.value | p.value | Keterangan |
|---|---|---|---|
| 0.59143 | 7.5 | 0.00617 | Signifikan |
Hasil uji Likelihood Ratio menunjukkan bahwa nilai parameter autoregresif spasial (ρ) sebesar 0,59143 dengan p-value 0,00617, yang berarti signifikan pada tingkat kepercayaan 95%. Temuan ini mengindikasikan adanya ketergantungan spasial positif pada variabel DBD IR, di mana peningkatan nilai DBD IR di suatu kabupaten/kota cenderung diikuti oleh peningkatan nilai DBD IR di wilayah sekitarnya. Selanjutnya hasil perhitungan impacts pada model SAR ditampilkan dalam tabel berikut:
| Variabel | Jenis_Impact | Impact Measures | p-value | Keputusan |
|---|---|---|---|---|
| Banjir dy/dx | Direct | -10.1734 | 0.035 | Signifikan |
| Cuaca_Ekstreme dy/dx | Direct | 4.1026 | 0.466 | Tidak Signifikan |
| Kepadatan_Penduduk dy/dx | Direct | 0.0055 | 0.275 | Tidak Signifikan |
| Persentase_Miskin dy/dx | Direct | -3.0953 | 0.454 | Tidak Signifikan |
| Persentase_Sanitasi dy/dx | Direct | -0.4462 | 0.481 | Tidak Signifikan |
| Banjir dy/dx | Indirect | -11.7097 | 0.364 | Tidak Signifikan |
| Cuaca_Ekstreme dy/dx | Indirect | 4.7221 | 0.626 | Tidak Signifikan |
| Kepadatan_Penduduk dy/dx | Indirect | 0.0063 | 0.626 | Tidak Signifikan |
| Persentase_Miskin dy/dx | Indirect | -3.5627 | 0.719 | Tidak Signifikan |
| Persentase_Sanitasi dy/dx | Indirect | -0.5136 | 0.689 | Tidak Signifikan |
Berdasarkan hasil estimasi direct dan indirect impact pada model SAR, variabel Banjir memiliki pengaruh langsung (direct effect) yang signifikan terhadap DBD IR, dengan nilai impact measures sebesar −10,1734 dan p-value 0,032, yang menunjukkan bahwa peningkatan kejadian banjir di suatu wilayah berkaitan dengan penurunan nilai DBD IR di wilayah tersebut. Sementara itu, seluruh pengaruh tidak langsung (indirect effect) dari variabel independen, termasuk banjir, cuaca ekstrem, kepadatan penduduk, persentase penduduk miskin, dan persentase sanitasi, tidak menunjukkan signifikansi statistik, yang mengindikasikan bahwa perubahan variabel-variabel tersebut di wilayah sekitar belum memberikan efek spillover yang signifikan terhadap DBD IR di wilayah lain. Hasil ini menunjukkan bahwa pengaruh variabel penjelas dalam model SAR lebih dominan terjadi secara lokal (langsung) dibandingkan melalui mekanisme spasial tidak langsung.
4.5.2 Uji Asumsi Model
| Jenis Uji | Statistik | Statistik Hitung | p-value | Kesimpulan |
|---|---|---|---|---|
| Linearitas (Ramsey RESET) | RESET | 1.82329 | 0.17953 | Hubungan variabel independen terhadap variabel dependen membentuk model linier |
| Normalitas (Shapiro Wilk) | W | 0.94267 | 0.08127 | Residual SAR berdistribusi normal |
| Homoskedastisitas (Breusch–Pagan) | BP | 1.32434 | 0.24981 | Tidak terdapat gejala heteroskedastisitas pada residual SAR |
Berdasarkan tabel di atas, hasil uji asumsi klasik pada model SAR menunjukkan bahwa model telah memenuhi asumsi linearitas, normalitas residual, dan homoskedastisitas, yang ditunjukkan oleh nilai p-value > 0,05 pada seluruh pengujian.
| I | p (Monte Carlo) | Keterangan |
|---|---|---|
| 0.0576 | 0.2460 | Tidak terdapat autokorelasi spasial residual |
Selanjutnya berdasarkan uji Moran’s I di atas, residual model SAR menunjukkan tidak adanya autokorelasi spasial residual, sehingga dapat disimpulkan bahwa model SAR telah mampu menangkap ketergantungan spasial pada data secara efektif.
4.5.3 Uji Sensitivitas Bobot
Uji sensitivitas bobot dilakukan untuk menilai kestabilan hasil estimasi model SAR terhadap perbedaan definisi bobot spasial yang digunakan.
| Bobot Spasial | rho | AIC | Log-Likelihood |
|---|---|---|---|
| Queen | 0.591 | 380.53 | -182.26 |
| Rook | 0.591 | 380.53 | -182.26 |
| KNN-5 | 0.484 | 384.90 | -184.45 |
Hasil uji sensitivitas bobot spasial menunjukkan bahwa model SAR dengan bobot Queen dan Rook menghasilkan nilai parameter autoregresif spasial (ρ), AIC, dan log-likelihood yang sama, yang mengindikasikan bahwa struktur ketetanggaan berbasis kontiguitas memberikan hasil estimasi yang konsisten. Sebaliknya, penggunaan bobot KNN-5 menghasilkan nilai ρ yang lebih kecil serta AIC yang lebih besar dan log-likelihood yang lebih rendah, yang menunjukkan kinerja model yang relatif lebih lemah. Berdasarkan kriteria kelayakan model, bobot spasial Queen dipilih sebagai bobot yang paling sesuai karena memberikan model dengan kinerja terbaik dan stabil.
4.6 Pemodelan Interpolasi
4.6.1 IDW
Hasil interpolasi menggunakan metode Inverse Distance Weighting (IDW) menunjukkan bahwa nilai DBD IR memiliki pola spasial yang tidak merata, dengan konsentrasi nilai tinggi pada beberapa wilayah tertentu dan nilai yang lebih rendah di wilayah lainnya. Pola ini mengindikasikan bahwa daerah dengan nilai DBD IR tinggi cenderung memengaruhi wilayah sekitarnya, sesuai dengan prinsip dasar IDW yang memberikan bobot lebih besar pada titik pengamatan terdekat. Permukaan interpolasi yang dihasilkan relatif halus namun masih menunjukkan variasi lokal, sehingga metode IDW mampu merepresentasikan pola sebaran spasial DBD IR berdasarkan kedekatan antarwilayah.
4.6.2 Kriging
Hasil variogram empiris DBD IR menunjukkan bahwa nilai semivarian meningkat seiring bertambahnya jarak hingga mencapai kondisi relatif stabil, yang mengindikasikan adanya autokorelasi spasial pada jarak tertentu. Variogram tersebut dapat direpresentasikan dengan baik menggunakan model spherical, yang menunjukkan bahwa pengaruh spasial melemah setelah jarak tertentu.
| Parameter | Model | Nilai | Range (m) |
|---|---|---|---|
| Partial Sill | Sph | 4131.582 | 17847.82 |
Parameter variogram menunjukkan nilai partial sill sebesar 4.131,582 dan range sekitar 17.847,82 meter, yang berarti bahwa pengaruh spasial DBD IR masih kuat hingga jarak tersebut dan menjadi relatif tidak berkorelasi setelahnya. Selanjutnya metode Ordinary Kriging digunakan karena mengasumsikan nilai rata-rata yang konstan namun tidak diketahui, sehingga sesuai untuk data DBD IR yang tidak menunjukkan tren spasial yang jelas dan menghasilkan estimasi yang fleksibel serta stabil.
Hasil interpolasi menggunakan Ordinary Kriging menunjukkan bahwa sebaran DBD IR memiliki pola spasial yang relatif halus dengan konsentrasi nilai tinggi pada beberapa wilayah tertentu, sementara sebagian besar wilayah lainnya memiliki nilai yang lebih rendah. Pola ini mencerminkan kemampuan Kriging dalam menangkap struktur autokorelasi spasial berdasarkan variogram, sehingga menghasilkan permukaan prediksi yang lebih terstruktur dibandingkan metode deterministik. Hasil Kriging ini memberikan gambaran sebaran DBD IR yang mempertimbangkan hubungan spasial antarwilayah secara lebih sistematis.
4.7 Interpretasi dan Pembahasan Hasil
Berdasarkan peta prediksi model Spatial Autoregressive (SAR), terlihat bahwa nilai DBD IR di Provinsi Sumatera Utara menunjukkan variasi spasial yang cukup jelas. Beberapa kabupaten/kota diprediksi memiliki nilai DBD IR yang relatif tinggi, terutama pada wilayah tertentu di bagian tengah dan timur provinsi, sementara wilayah lainnya menunjukkan nilai yang lebih rendah. Pola ini mengindikasikan adanya pengelompokan spasial, di mana wilayah dengan tingkat DBD IR tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki nilai tinggi.
Selanjutnya, peta residual standar model SAR menunjukkan bahwa sebagian besar kabupaten/kota memiliki residual yang relatif kecil dan tersebar secara acak. Residual bernilai positif (biru) mengindikasikan bahwa model cenderung mengunderestimasi nilai DBD IR di wilayah tersebut, sedangkan residual bernilai negatif (coklat) menunjukkan adanya overestimasi. Meskipun masih terdapat beberapa wilayah dengan residual yang cukup besar, secara umum tidak terlihat pola spasial tertentu pada residual, yang mengindikasikan bahwa model SAR telah mampu menangkap struktur spasial DBD IR di Provinsi Sumatera Utara dengan baik.
Secara matematis, model SAR yang terbentuk untuk DBD IR di Provinsi Sumatera Utara dapat dituliskan sebagai berikut:
\[ \hat{Y}_i = 118.27 - 8.94^{*} X_{1i} + 3.61 X_{2i} + 0.0048 X_{3i} - 2.72 X_{4i} - 0.39 X_{5i} + 0.59^{***} \sum_{j=1}^{n} w_{ij} Y_j \]
dengan \(Y_i\) menyatakan nilai incidence rate DBD pada kabupaten/kota ke-\(i\), \(X_1\) adalah kejadian banjir, \(X_2\) cuaca ekstrem, \(X_3\) kepadatan penduduk, \(X_4\) persentase penduduk miskin, dan \(X_5\) persentase sanitasi. Parameter \(\rho\) yang signifikan menunjukkan adanya ketergantungan spasial positif, di mana nilai incidence rate DBD suatu wilayah dipengaruhi oleh nilai incidence rate DBD di wilayah sekitarnya.
Hasil estimasi menunjukkan bahwa variabel banjir berpengaruh signifikan terhadap incidence rate DBD, sementara variabel cuaca ekstrem, kepadatan penduduk, persentase penduduk miskin, dan persentase sanitasi tidak menunjukkan pengaruh signifikan secara statistik. Temuan ini mengindikasikan bahwa variasi incidence rate DBD di Provinsi Sumatera Utara tidak hanya dipengaruhi oleh karakteristik wilayah itu sendiri, tetapi juga oleh kondisi wilayah tetangga melalui mekanisme spasial.
Perbandingan hasil interpolasi menunjukkan bahwa metode IDW menghasilkan permukaan prediksi dengan variasi nilai yang lebih kontras dan dipengaruhi kuat oleh titik pengamatan terdekat, sehingga membentuk pola lokal yang lebih menonjol. Sebaliknya, metode Kriging menghasilkan permukaan yang lebih halus dan terstruktur karena mempertimbangkan autokorelasi spasial melalui variogram, sehingga variasi ekstrem menjadi lebih teredam. Meskipun kedua metode menunjukkan pola sebaran DBD IR yang relatif serupa, Kriging memberikan representasi spasial yang lebih stabil dan konsisten, sedangkan IDW lebih sensitif terhadap sebaran dan kerapatan titik data.
BAB 5. KESIMPULAN DAN SARAN
5.1 Kesimpulan
Penelitian ini bertujuan untuk mengidentifikasi ketergantungan spasial antarwilayah serta memodelkan sebaran Incidence Rate (IR) DBD melalui pendekatan interpolasi spasial guna memperoleh gambaran permukaan spasial DBD yang lebih komprehensif. Hasil analisis menunjukkan bahwa DBD IR di Provinsi Sumatera Utara memiliki ketergantungan spasial yang signifikan, yang mengindikasikan bahwa kondisi DBD di suatu kabupaten/kota dipengaruhi oleh kondisi wilayah sekitarnya. Hal ini dibuktikan melalui uji autokorelasi spasial dan estimasi model Spatial Autoregressive (SAR) yang menghasilkan parameter spasial (ρ) signifikan.
Berdasarkan hasil estimasi model SAR, kejadian banjir terbukti berpengaruh signifikan terhadap DBD IR, sementara variabel lain tidak menunjukkan pengaruh signifikan secara statistik. Temuan ini menunjukkan bahwa faktor lingkungan memiliki peran penting dalam meningkatkan risiko penyebaran DBD. Selain itu, hasil interpolasi spasial menggunakan metode IDW dan Kriging mampu menggambarkan variasi spasial DBD IR secara lebih detail, dengan Kriging menghasilkan permukaan yang lebih halus dan terstruktur, sedangkan IDW menonjolkan variasi lokal berdasarkan kedekatan titik pengamatan. Secara keseluruhan, kombinasi analisis spasial dan interpolasi memberikan pemahaman yang lebih mendalam mengenai pola sebaran DBD antarwilayah kabupaten/kota di Provinsi Sumatera Utara.
5.2 Saran
Berdasarkan hasil penelitian, pemerintah daerah dan pemangku kebijakan diharapkan dapat menyusun kebijakan yang lebih bijak dan terintegrasi terkait pengelolaan lingkungan, khususnya dalam pengendalian deforestasi dan kebakaran hutan/lahan (karhutla). Secara teoritis, kerusakan lingkungan tersebut dapat meningkatkan risiko banjir, yang selanjutnya berpotensi memperbesar kejadian DBD melalui peningkatan genangan air sebagai tempat berkembang biaknya vektor penyakit.
Selain itu, hasil penelitian ini dapat digunakan sebagai dasar pertimbangan dalam perencanaan mitigasi dan pengendalian DBD berbasis wilayah, dengan memperhatikan keterkaitan antarwilayah dan pola spasial yang terbentuk. Untuk penelitian selanjutnya, disarankan untuk menambahkan variabel lingkungan dan iklim yang lebih rinci, serta menggunakan data dengan resolusi spasial dan temporal yang lebih tinggi guna memperoleh hasil pemodelan yang lebih akurat.
Referensi
Kementerian Kesehatan Republik Indonesia. (2023). Profil Kesehatan Indonesia.
World Health Organization (WHO). (2012). Global Strategy for Dengue Prevention and Control 2012–2020.
Badan Nasional Penanggulangan Bencana (BNPB). (2025). Geoportal data bencana Indonesia: Data kejadian banjir dan cuaca ekstrem. Diakses pada 15 Desember 2025 dari https://gis.bnpb.go.id
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic Publishers.
Anselin, L., Bera, A. K., Florax, R., & Yoon, M. J. (1996). Simple diagnostic tests for spatial dependence. Regional Science and Urban Economics.
LeSage, J. P., & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press.
Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control.
Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics.
Burrough, P. A., & McDonnell, R. A. (1998). Principles of Geographical Information Systems. Oxford University Press.
Cressie, N. (1993). Statistics for Spatial Data. Wiley.
Badan Pusat Statistik Provinsi Sumatera Utara. (2024). Kasus penyakit menurut kabupaten/kota dan jenis penyakit di Provinsi Sumatera Utara. Diakses pada 15 Desember 2025 dari https://sumut.bps.go.id
Badan Nasional Penanggulangan Bencana. (2024). Data kejadian bencana Indonesia menurut kabupaten/kota. Diakses pada 15 Desember 2025 dari https://data.bnpb.go.id
Badan Pusat Statistik Provinsi Sumatera Utara. (2024). Kepadatan penduduk dan rasio jenis kelamin menurut kabupaten/kota di Provinsi Sumatera Utara. Diakses pada 15 Desember 2025 dari https://sumut.bps.go.id
Badan Pusat Statistik Provinsi Sumatera Utara. (2024). Persentase penduduk miskin menurut kabupaten/kota di Provinsi Sumatera Utara. Diakses pada 15 Desember 2025 dari https://sumut.bps.go.id
Badan Pusat Statistik Provinsi Sumatera Utara. (2025). Provinsi Sumatera Utara dalam Angka 2025. Medan: BPS Provinsi Sumatera Utara.
Lampiran
Dokumen pendukung lain yang diperlukan adalah sebagai berikut:
Syntax
# 1. SET WORKING DIRECTORY ----------------------------------------------------------
setwd("C:/Users/ACER/Documents/SEMESTER 5/SPATIAL/UAS")
getwd()
# 2. LOAD LIBRARY ----------------------------------------------------------
library(sf); library(dplyr); library(stringr)
library(spdep); library(spatialreg); library(leaflet)
library(tmap); library(broom); library(ggplot2)
library(psych); library(readxl); library(viridis)
library(scales); library(tibble); library(purrr)
library(forcats); library(htmltools);library(knitr)
library(kableExtra); library(gstat); library(sp)
library(raster); library(patchwork)
# 3. BACA SHAPEFILE KAB/KOTA ----------------------------------------------------------
shp_kabkota <- st_read(
"batas-administrasi-indonesia-master/batas-administrasi-indonesia-master/Kab_Kota/Kab_Kota.shp"
)
# 4. FILTER KAB/KOTA SUMATERA UTARA ---------------------------------------
kab_sumut <- c(
"Asahan","Batu Bara","Dairi","Deli Serdang","Humbang Hasundutan","Karo",
"Labuhanbatu","Labuhanbatu Selatan","Labuhanbatu Utara","Langkat",
"Mandailing Natal","Nias","Nias Barat","Nias Selatan","Nias Utara",
"Padang Lawas","Padang Lawas Utara","Pakpak Bharat","Samosir",
"Serdang Bedagai","Simalungun","Tapanuli Selatan","Tapanuli Tengah",
"Tapanuli Utara","Toba","Kota Binjai","Kota Gunungsitoli","Kota Medan",
"Kota Padang Sidempuan","Kota Pematangsiantar","Kota Sibolga",
"Kota Tanjung Balai","Kota Tebing Tinggi"
)
shp_sumut <- shp_kabkota %>%
filter(KAB_KOTA %in% kab_sumut)
nrow(shp_sumut)
# 5. BACA DATA BANJIR ----------------------------------------------------------
data_sumut <- read_excel("data_sumut.xlsx")
# Cek kolom
names(data_sumut)
# 6. NORMALISASI NAMA KAB/KOTA --------------------------------------------------
shp_sumut$KAB_KOTA <- str_to_title(shp_sumut$KAB_KOTA)
data_sumut$KAB_KOTA <- str_to_title(data_sumut$KAB_KOTA)
# 7. JOIN SHAPEFILE + DATA BANJIR -----------------------------------------
shp_sumut_join <- shp_sumut %>%
left_join(
data_sumut,
by = "KAB_KOTA"
)
# 8. CEK HASIL JOIN -------------------------------------------------------
# Pastikan tidak ada NA di variabel utama
summary(shp_sumut_join)
head(shp_sumut_join)
# 11. AUTOKORELASI SPASIAL -------------------------------------------------------
# PASTIKAN SHP AMAN
shp_sumut_join <- st_make_valid(shp_sumut_join)
shp_sumut_join <- shp_sumut_join[!st_is_empty(shp_sumut_join), ]
shp_sumut_join <- st_cast(shp_sumut_join, "MULTIPOLYGON", warn = FALSE)
# Bangun Queen contiguity
nb_queen <- poly2nb(shp_sumut_join, queen = TRUE)
# Ringkasan tetangga
summary(nb_queen)
# Matriks bobot spasial (row-standardized)
lw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
# Moran's I untuk BANJIR
moran_banjir <- moran.test(
shp_sumut_join$Banjir,
lw_queen,
zero.policy = TRUE
)
moran_banjir
# Queen Contiguity
make_queen_listw <- function(sf_obj, queen = TRUE, style = "W", zero.policy = TRUE){
nb <- spdep::poly2nb(sf_obj, queen = queen)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_queen <- make_queen_listw(shp_sumut_join)
lw_queen <- w_queen$lw
# KNN Contiguity
make_knn_listw <- function(sf_obj, k = 5, style = "W", zero.policy = TRUE){
pts <- sf::st_point_on_surface(sf_obj) # lebih aman dari centroid yang bisa keluar polygon
coords <- sf::st_coordinates(pts)
knn <- spdep::knearneigh(coords, k = k)
nb <- spdep::knn2nb(knn)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_knn5 <- make_knn_listw(shp_sumut_join, k = 5)
lw_knn5 <- w_knn5$lw
# Rook Contiguity
make_rook_listw <- function(sf_obj, style = "W", zero.policy = TRUE){
nb <- spdep::poly2nb(sf_obj, queen = FALSE)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_rook <- make_rook_listw(shp_sumut_join)
lw_rook <- w_rook$lw
# Identifikasi Masalah
choropleth_map <- function(sf_obj, var,
title = NULL,
subtitle = NULL,
pal = "viridis",
direction = 1){
stopifnot(var %in% names(sf_obj))
if (!is.numeric(sf_obj[[var]]))
stop("Variabel harus numeric: ", var)
if (is.null(title)) title <- var
ggplot2::ggplot(sf_obj) +
ggplot2::geom_sf(
ggplot2::aes(fill = .data[[var]]),
color = "white",
linewidth = 0.25
) +
ggplot2::scale_fill_viridis_c(
option = pal,
direction = direction,
na.value = "grey90"
) +
ggplot2::labs(
title = title,
subtitle = subtitle,
fill = NULL
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(
axis.text = ggplot2::element_blank(),
axis.title = ggplot2::element_blank(),
panel.grid = ggplot2::element_blank(),
legend.position = "right"
)
}
P_DBD_Awal <- choropleth_map(
shp_sumut_join, "DBD_IR",
title = expression(paste("Peta ", italic("Incidence Rate"), " DBD Sumatera Utara 2024"))
)
P_DBD_Awal <- P_DBD_Awal +
labs(fill = "DBD IR\n(per 100.000 penduduk)")
P_DBD_Awal
# Data
## Tabel Definisi Operasional
var_df <- data.frame(
Variabel = c(
"Angka Kesakitan DBD",
"Banjir",
"Cuaca Ekstrem",
"Kepadatan Penduduk",
"Persentase Penduduk Miskin",
"Akses Sanitasi Layak"
),
Simbol = c("Y", "X₁", "X₂", "X₃", "X₄", "X₅"),
`Definisi Operasional` = c(
"Jumlah kasus Demam Berdarah Dengue per 100.000 penduduk menurut kabupaten/kota.",
"Jumlah kejadian banjir yang terjadi dalam satu tahun menurut kabupaten/kota.",
"Jumlah kejadian cuaca ekstrem yang terjadi dalam satu tahun menurut kabupaten/kota.",
"Jumlah penduduk per satuan luas wilayah menurut kabupaten/kota.",
"Persentase penduduk yang berada di bawah garis kemiskinan menurut kabupaten/kota.",
"Persentase rumah tangga yang memiliki akses terhadap fasilitas sanitasi layak."
),
Satuan = c(
"Kasus per 100.000 penduduk",
"Kejadian",
"Kejadian",
"Jiwa/km²",
"%",
"%"
)
)
var_df |>
kable(
caption = "Tabel 1. Definisi Operasional Variabel",
align = c("l", "c", "l", "c"),
booktabs = TRUE
) |>
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
column_spec(1, width = "22%") |>
column_spec(2, width = "8%") |>
column_spec(3, width = "50%") |>
column_spec(4, width = "20%")
# Analisis Statistik Deskriptif
# Statistik Deskriptif
vars_desc <- c(
"DBD_IR",
"Banjir",
"Cuaca_Ekstreme",
"Kepadatan_Penduduk",
"Persentase_Miskin",
"Persentase_Sanitasi"
)
# 2) ubah data jadi long (AMAN meski nama kolom ada underscore)
dat_long <- shp_sumut_join |>
st_drop_geometry() |>
dplyr::select(dplyr::all_of(vars_desc)) |>
pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")
# 3) hitung statistik deskriptif per variabel (1 baris = 1 variabel)
desc_table <- dat_long |>
group_by(Variabel) |>
summarise(
Mean = mean(Nilai, na.rm = TRUE),
Minimum = min(Nilai, na.rm = TRUE),
Maximum = max(Nilai, na.rm = TRUE),
Simpangan_Baku = sd(Nilai, na.rm = TRUE),
.groups = "drop"
) |>
mutate(across(where(is.numeric), ~round(.x, 2)))
# 4) (opsional tapi rapi) ganti label variabel biar enak dibaca
desc_table <- desc_table |>
mutate(Variabel = case_when(
Variabel == "DBD_IR" ~ "Incidence Rate DBD",
Variabel == "Cuaca_Ekstreme" ~ "Cuaca Ekstrem",
Variabel == "Kepadatan_Penduduk" ~ "Kepadatan Penduduk",
Variabel == "Persentase_Miskin" ~ "Persentase Penduduk Miskin",
Variabel == "Persentase_Sanitasi" ~ "Akses Sanitasi Layak",
TRUE ~ Variabel
))
# 5) tampilkan tabel dengan kable styling
desc_table |>
kable(
caption = "Tabel 2. Hasil Analisis Deskriptif",
col.names = c("Variabel", "Mean", "Minimum", "Maximum", "Simpangan Baku"),
align = c("l", "c", "c", "c", "c"),
booktabs = TRUE
) |>
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
)
# Gambaran Awal Peta
choropleth_map <- function(sf_obj, var,
title = NULL,
subtitle = NULL,
pal = "viridis",
direction = 1){
stopifnot(var %in% names(sf_obj))
if (!is.numeric(sf_obj[[var]]))
stop("Variabel harus numeric: ", var)
if (is.null(title)) title <- var
ggplot2::ggplot(sf_obj) +
ggplot2::geom_sf(
ggplot2::aes(fill = .data[[var]]),
color = "white",
linewidth = 0.25
) +
ggplot2::scale_fill_viridis_c(
option = pal,
direction = direction,
na.value = "grey90"
) +
ggplot2::labs(
title = title,
subtitle = subtitle,
fill = NULL
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(
axis.text = ggplot2::element_blank(),
axis.title = ggplot2::element_blank(),
panel.grid = ggplot2::element_blank(),
legend.position = "right"
)
}
P_DBD <- choropleth_map(
shp_sumut_join, "DBD_IR",
title = "DBD IR"
)
P_Banjir <- choropleth_map(
shp_sumut_join, "Banjir",
title = "Banjir"
)
P_Cuaca <- choropleth_map(
shp_sumut_join, "Cuaca_Ekstreme",
title = "Cuaca Ekstreme"
)
P_Kepadatan <- choropleth_map(
shp_sumut_join, "Kepadatan_Penduduk",
title = "Kepadatan Penduduk"
)
P_Miskin <- choropleth_map(
shp_sumut_join, "Persentase_Miskin",
title = "Penduduk Miskin (%)"
)
P_Sanitasi <- choropleth_map(
shp_sumut_join, "Persentase_Sanitasi",
title = "Sanitasi (%)"
)
fix_panel <- function(p){
p +
ggplot2::theme(
plot.title = ggplot2::element_text(size = 12, face = "bold"),
plot.subtitle = ggplot2::element_text(size = 9),
legend.title = ggplot2::element_text(size = 10),
legend.text = ggplot2::element_text(size = 9),
plot.margin = ggplot2::margin(2, 2, 2, 2)
)
}
P_DBD2 <- fix_panel(P_DBD)
P_Banjir2 <- fix_panel(P_Banjir)
P_Cuaca2 <- fix_panel(P_Cuaca)
P_Kepadatan2 <- fix_panel(P_Kepadatan)
P_Miskin2 <- fix_panel(P_Miskin)
P_Sanitasi2 <- fix_panel(P_Sanitasi)
peta6 <- gridExtra::grid.arrange(
P_DBD2, P_Banjir2, P_Cuaca2,
P_Kepadatan2, P_Miskin2, P_Sanitasi2,
ncol = 3
)
ggplot2::ggsave("Peta_6panel.png", peta6,
width = 14, height = 8, units = "in", dpi = 300, bg = "white")
# Ordinary Least Square
## Uji Asumsi Klasik
ols <- lm(DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi, data = shp_sumut_join)
summary(ols)
# Moran's I residual OLS
shp_sumut_join$res_ols <- residuals(ols)
# Moran's I Residual model OLS
moran_res_ols <- moran.test(
shp_sumut_join$res_ols,
lw_queen,
zero.policy = TRUE
)
moran_res_ols
# Uji Asumsi Klasik
# Normalitas
res <- residuals(ols)
shapiro.test(res)
# Multikolinearitas
library(car)
vif(ols)
# Heteroskedastisitas
library(lmtest)
bptest(ols)
# Linearitas
library(lmtest)
resettest(ols)
uji_asumsi_klasik_kable <- function(model,
reset_power = 2:3,
alpha = 0.05,
caption = "Hasil Uji Asumsi Klasik",
digits = 5) {
# --- validasi ---
if (!inherits(model, "lm")) stop("model harus objek 'lm' (OLS).")
# paket wajib
if (!requireNamespace("lmtest", quietly = TRUE)) stop("Paket 'lmtest' belum terpasang.")
if (!requireNamespace("knitr", quietly = TRUE)) stop("Paket 'knitr' belum terpasang.")
if (!requireNamespace("kableExtra", quietly = TRUE)) stop("Paket 'kableExtra' belum terpasang.")
# residual
res <- residuals(model)
# --- 1) Linearitas: Ramsey RESET ---
rst <- lmtest::resettest(model, power = reset_power, type = "fitted")
rst_stat <- unname(as.numeric(rst$statistic[[1]]))
rst_p <- rst$p.value
kes_reset <- if (rst_p > alpha) {
"Variabel independen terhadap variabel dependen membentuk model linier"
} else {
"Indikasi ketidaklinieran / spesifikasi model belum tepat"
}
# --- 2) Normalitas: Shapiro-Wilk ---
shp <- shapiro.test(res)
shp_stat <- unname(as.numeric(shp$statistic[["W"]]))
shp_p <- shp$p.value
kes_shp <- if (shp_p > alpha) {
"Residual berdistribusi normal"
} else {
"Residual tidak berdistribusi normal"
}
# --- 3) Homoskedastisitas: Breusch-Pagan ---
bp <- lmtest::bptest(model)
bp_stat <- unname(as.numeric(bp$statistic[["BP"]]))
bp_p <- bp$p.value
kes_bp <- if (bp_p > alpha) {
"Tidak terdapat gejala heteroskedastisitas"
} else {
"Terdapat gejala heteroskedastisitas"
}
# helper format angka
fmt_num <- function(x) formatC(x, digits = digits, format = "f")
fmt_p <- function(p) formatC(p, digits = digits, format = "f")
# --- tabel ringkas ---
tbl <- data.frame(
`Jenis Uji` = c(
"Linearitas (Ramsey RESET)",
"Normalitas (Shapiro Wilk)",
"Homoskedastisitas (Breusch–Pagan)"
),
Statistik = c("RESET", "W", "BP"),
`Statistik Hitung` = c(
fmt_num(rst_stat),
fmt_num(shp_stat),
fmt_num(bp_stat)
),
`p-value` = c(
fmt_p(rst_p),
fmt_p(shp_p),
fmt_p(bp_p)
),
Keterangan = c(kes_reset, kes_shp, kes_bp),
check.names = FALSE
)
# --- output kable ---
knitr::kable(
tbl,
format = "html",
escape = FALSE,
caption = caption,
align = c("l", "c", "c", "c", "l")
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
kableExtra::column_spec(1, width = "28%") |>
kableExtra::column_spec(2, width = "10%") |>
kableExtra::column_spec(3, width = "14%") |>
kableExtra::column_spec(4, width = "10%") |>
kableExtra::column_spec(5, width = "38%")
}
hasil <- uji_asumsi_klasik_kable(
ols,
caption = "Tabel 3. Hasil Uji Asumsi Klasik"
)
hasil
## Tabel Multikolinearitas
tabel_vif_horizontal_kable <- function(model,
cutoff = 10,
caption = "Hasil Uji Multikolinearitas (VIF)",
digits = 4) {
if (!inherits(model, "lm")) stop("model harus objek 'lm' (OLS).")
if (!requireNamespace("car", quietly = TRUE)) stop("Paket 'car' belum terpasang.")
if (!requireNamespace("knitr", quietly = TRUE)) stop("Paket 'knitr' belum terpasang.")
if (!requireNamespace("kableExtra", quietly = TRUE)) stop("Paket 'kableExtra' belum terpasang.")
vif_raw <- car::vif(model)
# handle GVIF jika ada faktor
if (is.matrix(vif_raw)) {
vif_val <- vif_raw[, "GVIF"]^(1/(2 * vif_raw[, "Df"]))
} else {
vif_val <- as.numeric(vif_raw)
}
n <- length(vif_val)
xnames <- paste0("X", seq_len(n))
fmt <- function(x) formatC(x, digits = digits, format = "f")
# kesimpulan umum berdasarkan max VIF
max_vif <- max(vif_val, na.rm = TRUE)
ket <- if (max_vif < cutoff) "Tidak terdapat multikolinearitas" else "Terdapat multikolinearitas"
# buat 1 baris, kolomnya X1..Xn + Keterangan
tbl <- as.data.frame(t(setNames(fmt(vif_val), xnames)), check.names = FALSE)
tbl$Keterangan <- ket
knitr::kable(
tbl,
format = "html",
caption = caption,
align = c(rep("c", n), "l"),
escape = FALSE
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
# kolom X1..Xn dipersempit
kableExtra::column_spec(1:n, width = "7%") |>
# kolom Keterangan dibatasi + boleh wrap
kableExtra::column_spec(n + 1, width = "26%", extra_css = "white-space: normal;")
}
tabel_vif_horizontal_kable(
model = ols,
cutoff = 10,
caption = "Tabel 4. Hasil Uji Multikolinearitas"
)
## Morans'I Residual
moran_residual_ols_kable <- function(model_ols,
lw,
nsim = 999,
alpha = 0.05,
seed = 123,
caption = "Moran's I Residual Model OLS (Monte Carlo)",
digits = 4) {
# validasi
if (!inherits(model_ols, "lm")) stop("model_ols harus objek OLS dari lm().")
if (!inherits(lw, "listw")) stop("lw harus objek listw.")
if (!requireNamespace("spdep", quietly = TRUE)) stop("Paket 'spdep' belum terpasang.")
if (!requireNamespace("knitr", quietly = TRUE)) stop("Paket 'knitr' belum terpasang.")
if (!requireNamespace("kableExtra", quietly = TRUE)) stop("Paket 'kableExtra' belum terpasang.")
# residual OLS
res_ols <- as.numeric(residuals(model_ols))
# Moran MC
set.seed(seed)
mor_mc <- spdep::moran.mc(res_ols, lw, nsim = nsim, zero.policy = TRUE)
I_val <- unname(as.numeric(mor_mc$statistic))
p_val <- mor_mc$p.value
ket <- if (p_val > alpha) {
"Tidak terdapat autokorelasi spasial residual"
} else {
"Masih terdapat autokorelasi spasial residual"
}
# tabel (1 baris)
tbl <- data.frame(
`I` = formatC(I_val, digits = digits, format = "f"),
`p(Monte Carlo)` = formatC(p_val, digits = digits, format = "f"),
`Keterangan` = ket,
check.names = FALSE
)
row.names(tbl) <- NULL
# kable style (ikut style kamu)
knitr::kable(
tbl,
format = "html",
escape = FALSE,
caption = caption,
align = c("c", "c", "l")
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
kableExtra::column_spec(1, width = "20%") |>
kableExtra::column_spec(2, width = "35%") |>
kableExtra::column_spec(3, width = "45%")
}
moran_residual_ols_kable(
model_ols = ols, # ganti sesuai nama model OLS kamu
lw = lw_queen,
nsim = 999,
alpha = 0.05,
caption = "Moran's I Residual Model OLS (Queen, Monte Carlo)"
)
# Autokorelasi Spasial
table_kable <- function(df,
caption = NULL,
digits = 2,
align = NULL,
striped = TRUE,
full_width = TRUE,
bold_header = TRUE){
library(knitr)
library(kableExtra)
if (is.null(align)) {
align <- paste0(
"l",
paste(rep("c", ncol(df) - 1), collapse = "")
)
}
kable(
df,
caption = caption,
digits = digits,
align = align,
booktabs = TRUE
) %>%
kable_styling(
full_width = full_width,
position = "center",
bootstrap_options = if (striped) "striped" else NULL
) %>%
{ if (bold_header) row_spec(., 0, bold = TRUE) else . }
}
## Function Moran's I
moran_table <- function(sf_obj, vars, lw, alpha = 0.05, zero.policy = TRUE){
results <- lapply(vars, function(var){
x <- sf_obj[[var]]
if (!is.numeric(x)) stop(paste("Variabel harus numeric:", var))
ok <- !is.na(x)
x_ok <- x[ok]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
moran_res <- spdep::moran.test(x_ok, lw_ok, zero.policy = zero.policy)
data.frame(
Variabel = var,
`Moran's I` = as.numeric(moran_res$estimate["Moran I statistic"]),
Z = as.numeric(moran_res$statistic),
`p.value` = moran_res$p.value,
Keterangan = ifelse(moran_res$p.value < alpha,
"Signifikan",
"Tidak Signifikan")
)
})
do.call(rbind, results)
}
vars_uji <- c("DBD_IR","Banjir", "Cuaca_Ekstreme", "Kepadatan_Penduduk","Persentase_Miskin", "Persentase_Sanitasi")
tabel_moran <- moran_table(
sf_obj = shp_sumut_join,
vars = vars_uji,
lw = lw_queen
)
## Geary's C
geary_table <- function(sf_obj, vars, lw, alpha = 0.05, zero.policy = TRUE){
one_var <- function(v){
x <- sf_obj[[v]]
if (!is.numeric(x)) stop("Variabel harus numeric: ", v)
ok <- is.finite(x) # buang NA/Inf
x <- x[ok]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
# Geary's C (asymptotic) - TANPA Monte Carlo
gt <- spdep::geary.test(x, lw_ok, zero.policy = zero.policy)
C <- unname(gt$estimate[["Geary C statistic"]])
Z <- unname(gt$statistic)
p <- gt$p.value
data.frame(
Variabel = v,
`Geary's C` = as.numeric(C),
Z = as.numeric(Z),
`p.value` = as.numeric(p),
Keterangan = ifelse(p < alpha, "Signifikan", "Tidak Signifikan"),
row.names = NULL
)
}
do.call(rbind, lapply(vars, one_var))
}
tab_geary <- geary_table(
sf_obj = shp_sumut_join,
vars = vars_uji,
lw = lw_queen
)
table_kable(
tab_geary,
caption = "Hasil Uji Autokorelasi Spasial (Geary's C)",
digits = 3)
## LISA
lisa_localmoran_map <- function(sf_obj, var, lw,
alpha = 0.05, nsim = 999,
zero.policy = TRUE, seed = 123,
title = NULL){
stopifnot(var %in% names(sf_obj))
x <- sf_obj[[var]]
if (!is.numeric(x)) stop("Variabel harus numeric: ", var)
# Buang NA/Inf dan subset listw biar konsisten
ok <- is.finite(x)
sf_ok <- sf_obj[ok, ]
x_ok <- x[ok]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
# Standardize (z-score)
z <- as.numeric(scale(x_ok))
wz <- as.numeric(spdep::lag.listw(lw_ok, z, zero.policy = zero.policy))
set.seed(seed)
# Prefer permutation-based local Moran jika tersedia
if ("localmoran_perm" %in% getNamespaceExports("spdep")) {
lm <- spdep::localmoran_perm(z, lw_ok, nsim = nsim, zero.policy = zero.policy)
Ii <- lm[, "Ii"]
p <- lm[, "Pr(z != E(Ii))"]
} else {
lm <- spdep::localmoran(z, lw_ok, zero.policy = zero.policy)
Ii <- lm[, "Ii"]
# p-value approx normal (kurang kuat daripada permutation)
p <- 2 * pnorm(abs(lm[, "Z.Ii"]), lower.tail = FALSE)
}
sig <- p < alpha
# Quadrant / Cluster
quad <- ifelse(z >= 0 & wz >= 0, "HH",
ifelse(z < 0 & wz < 0, "LL",
ifelse(z >= 0 & wz < 0, "HL", "LH")))
cluster <- ifelse(sig, quad, "Not sig")
cluster <- factor(cluster, levels = c("HH","LL","HL","LH","Not sig"))
# Simpan hasil ke sf
sf_ok$z <- z
sf_ok$wz <- wz
sf_ok$Ii <- as.numeric(Ii)
sf_ok$pvalue <- as.numeric(p)
sf_ok$cluster <- cluster
# Warna yang umum dipakai di LISA maps
pal <- c(
"HH" = "#d7191c", # hot spot
"LL" = "#2c7bb6", # cold spot
"HL" = "#fdae61", # high-low (outlier)
"LH" = "#abd9e9", # low-high (outlier)
"Not sig" = "grey85"
)
# Ringkasan count untuk anotasi
tab <- table(sf_ok$cluster)
ann <- paste0(names(tab), ": ", as.integer(tab), collapse = " | ")
if (is.null(title)) title <- paste0(var)
p_map <- ggplot2::ggplot(sf_ok) +
ggplot2::geom_sf(ggplot2::aes(fill = cluster), color = "white", linewidth = 0.25) +
ggplot2::scale_fill_manual(values = pal, drop = FALSE) +
ggplot2::labs(
title = title,
subtitle = paste0(ann),
fill = "Cluster"
) +
ggplot2::theme_void(base_size = 12) +
ggplot2::theme(
plot.title = ggplot2::element_text(face = "bold"),
legend.position = "right",
panel.grid = ggplot2::element_blank()
)
list(data = sf_ok, plot = p_map)
}
res_lisa_DBD <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "DBD_IR", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
res_lisa_Banjir <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "Banjir", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
res_lisa_Cuaca <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "Cuaca_Ekstreme", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
res_lisa_Kepadatan <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "Kepadatan_Penduduk", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
res_lisa_Kemiskinan <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "Persentase_Miskin", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
res_lisa_Sanitasi <- lisa_localmoran_map(
sf_obj = shp_sumut_join,
var = "Persentase_Sanitasi", # ganti sesuai nama kolom Y kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999
)
LISA_DBD <- res_lisa_DBD$plot
LISA_Banjir <- res_lisa_Banjir$plot
LISA_Cuaca <- res_lisa_Cuaca$plot
LISA_Kepadatan <- res_lisa_Kepadatan$plot
LISA_Kemiskinan <- res_lisa_Kemiskinan$plot
LISA_Sanitasi <- res_lisa_Sanitasi$plot
LISA_DBD
LISA_Banjir
LISA_Cuaca
LISA_Kepadatan
LISA_Kemiskinan
LISA_Sanitasi
fix_panel <- function(p){
p +
ggplot2::theme(
plot.title = ggplot2::element_text(size = 12, face = "bold"),
plot.subtitle = ggplot2::element_text(size = 9),
legend.title = ggplot2::element_text(size = 10),
legend.text = ggplot2::element_text(size = 9),
plot.margin = ggplot2::margin(2, 2, 2, 2)
)
}
LISA_DBD2 <- fix_panel(LISA_DBD)
LISA_Banjir2 <- fix_panel(LISA_Banjir)
LISA_Cuaca2 <- fix_panel(LISA_Cuaca)
LISA_Kepadatan2 <- fix_panel(LISA_Kepadatan)
LISA_Kemiskinan2 <- fix_panel(LISA_Kemiskinan)
LISA_Sanitasi2 <- fix_panel(LISA_Sanitasi)
p6 <- gridExtra::grid.arrange(
LISA_DBD2, LISA_Banjir2, LISA_Cuaca2,
LISA_Kepadatan2, LISA_Kemiskinan2, LISA_Sanitasi2,
ncol = 3
)
ggplot2::ggsave("LISA_6panel.png", p6,
width = 14, height = 8, units = "in", dpi = 300, bg = "white")
### LISA LEAFLET
lisa_leaflet <- function(sf_obj, var, lw,
alpha = 0.05, nsim = 999,
zero.policy = TRUE, seed = 123,
id_col = NULL, title = NULL){
stopifnot(var %in% names(sf_obj))
if (!is.numeric(sf_obj[[var]])) stop("Variabel harus numeric: ", var)
# leaflet perlu lon/lat
sf_ll <- sf::st_transform(sf_obj, 4326)
x <- sf_ll[[var]]
ok <- is.finite(x)
sf_ok <- sf_ll[ok, ]
x_ok <- x[ok]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
# z-score + spatial lag of z
z <- as.numeric(scale(x_ok))
wz <- as.numeric(spdep::lag.listw(lw_ok, z, zero.policy = zero.policy))
set.seed(seed)
# Permutation-based local Moran jika tersedia, kalau tidak fallback
if ("localmoran_perm" %in% getNamespaceExports("spdep")) {
lm <- spdep::localmoran_perm(z, lw_ok, nsim = nsim, zero.policy = zero.policy)
Ii <- lm[, "Ii"]
p <- lm[, "Pr(z != E(Ii))"]
} else {
lm <- spdep::localmoran(z, lw_ok, zero.policy = zero.policy)
Ii <- lm[, "Ii"]
p <- 2 * pnorm(abs(lm[, "Z.Ii"]), lower.tail = FALSE)
}
sig <- p < alpha
quad <- ifelse(z >= 0 & wz >= 0, "HH",
ifelse(z < 0 & wz < 0, "LL",
ifelse(z >= 0 & wz < 0, "HL", "LH")))
cluster <- ifelse(sig, quad, "Not sig")
cluster <- factor(cluster, levels = c("HH","LL","HL","LH","Not sig"))
# simpan hasil
sf_ok$z <- z
sf_ok$wz <- wz
sf_ok$Ii <- as.numeric(Ii)
sf_ok$pvalue <- as.numeric(p)
sf_ok$cluster <- cluster
# nama wilayah
if (!is.null(id_col) && id_col %in% names(sf_ok)) {
nm <- sf_ok[[id_col]]
} else {
nm <- paste("Area", seq_len(nrow(sf_ok)))
}
sf_ok$._name <- nm
# palette LISA umum
pal <- leaflet::colorFactor(
palette = c(
"HH" = "#d7191c", # hotspot
"LL" = "#2c7bb6", # coldspot
"HL" = "#fdae61", # high surrounded by low
"LH" = "#abd9e9", # low surrounded by high
"Not sig" = "#d9d9d9"
),
domain = sf_ok$cluster
)
fmt3 <- function(v) ifelse(is.na(v), "-", format(round(v, 3), nsmall = 3))
sf_ok$._tooltip <- paste0(
"<b>", sf_ok$._name, "</b><br/>",
var, ": ", sf_ok[[var]], "<br/>",
"Local I: ", fmt3(sf_ok$Ii), "<br/>",
"p-value: ", fmt3(sf_ok$pvalue), "<br/>",
"Cluster: <b>", sf_ok$cluster, "</b>"
)
sf_ok$._popup <- paste0(
"<div style='font-size:14px'>",
"<h4 style='margin:0 0 6px 0'>", sf_ok$._name, "</h4>",
"<b>", var, ":</b> ", sf_ok[[var]], "<br/>",
"<b>z:</b> ", fmt3(sf_ok$z), " | <b>Wz:</b> ", fmt3(sf_ok$wz), "<br/>",
"<b>Local Moran's I:</b> ", fmt3(sf_ok$Ii), "<br/>",
"<b>p-value:</b> ", fmt3(sf_ok$pvalue), "<br/>",
"<b>Cluster:</b> ", sf_ok$cluster,
"</div>"
)
if (is.null(title)) {
title <- paste0("LISA - Local Moran's I (",var,")")
}
m <- leaflet::leaflet(sf_ok, options = leaflet::leafletOptions(preferCanvas = TRUE)) |>
leaflet::addProviderTiles(leaflet::providers$CartoDB.Positron) |>
leaflet::addPolygons(
fillColor = ~pal(cluster),
fillOpacity = 0.75,
color = "white",
weight = 1,
opacity = 1,
smoothFactor = 0.2,
highlightOptions = leaflet::highlightOptions(weight = 3, color = "#333", bringToFront = TRUE),
label = lapply(sf_ok$._tooltip, htmltools::HTML),
popup = lapply(sf_ok$._popup, htmltools::HTML)
) |>
leaflet::addLegend(
position = "bottomright",
pal = pal,
values = sf_ok$cluster,
title = title,
opacity = 1
)
list(data = sf_ok, map = m)
}
res_lisa_leaf <- lisa_leaflet(
sf_obj = shp_sumut_join,
var = "DBD_IR", # ganti sesuai nama kolommu
lw = lw_queen,
alpha = 0.05,
nsim = 999,
id_col = "KAB_KOTA" # opsional: ganti sesuai kolom nama kab/kota
)
res_lisa_leaf$map
## Getis Ord
getis_ord_map2 <- function(sf_obj, var, lw,
alpha = 0.05, nsim = 999,
alternative = c("two.sided", "greater", "less"),
zero.policy = TRUE, seed = 123,
title = NULL,
show_axes = FALSE,
add_subtitle = TRUE){
alternative <- match.arg(alternative)
stopifnot(var %in% names(sf_obj))
x <- sf_obj[[var]]
if (!is.numeric(x)) stop("Variabel harus numeric: ", var)
# Buang NA/Inf dan subset listw biar konsisten
ok <- is.finite(x)
sf_ok <- sf_obj[ok, ]
x_ok <- x[ok]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
set.seed(seed)
# Gi* permutation
gi <- spdep::localG_perm(x_ok, lw_ok, nsim = nsim, zero.policy = zero.policy)
# ---- Ambil Z dan p-value se-robust mungkin ----
# localG_perm biasanya return matrix/data.frame dengan kolom "Z.Ii" / "Pr(...)" tergantung versi,
# atau kadang hanya numeric vector z.
if (is.matrix(gi) || is.data.frame(gi)) {
gi_df <- as.data.frame(gi)
# coba tebak kolom z
z_col <- intersect(names(gi_df), c("Z.Ii","Z.Gi","Z", "GiZ", "statistic"))
if (length(z_col) == 0) {
# fallback: ambil kolom pertama
z <- as.numeric(gi_df[[1]])
} else {
z <- as.numeric(gi_df[[z_col[1]]])
}
# coba tebak kolom p
p_col <- intersect(names(gi_df), c("Pr(z != E(Ii))","Pr(z > E(Ii))","Pr(z < E(Ii))","p.value","pvalue"))
if (length(p_col) > 0) {
p_raw <- as.numeric(gi_df[[p_col[1]]])
# kalau p yang tersedia bukan two-sided, kita tetap pakai sesuai alternative jika bisa
p <- p_raw
} else {
# fallback normal approx
if (alternative == "two.sided") p <- 2 * pnorm(abs(z), lower.tail = FALSE)
if (alternative == "greater") p <- pnorm(z, lower.tail = FALSE)
if (alternative == "less") p <- pnorm(z, lower.tail = TRUE)
}
} else {
# gi hanya numeric vector (z-score)
z <- as.numeric(gi)
if (alternative == "two.sided") p <- 2 * pnorm(abs(z), lower.tail = FALSE)
if (alternative == "greater") p <- pnorm(z, lower.tail = FALSE)
if (alternative == "less") p <- pnorm(z, lower.tail = TRUE)
}
sig <- is.finite(p) & (p < alpha)
# Klasifikasi mirip LISA
cls <- ifelse(sig & z > 0, "Hotspot",
ifelse(sig & z < 0, "Coldspot", "Not sig"))
cls <- factor(cls, levels = c("Hotspot","Coldspot","Not sig"))
# Simpan ke sf
sf_ok$GiZ <- z
sf_ok$pvalue <- p
sf_ok$cluster <- cls
pal <- c(
"Hotspot" = "#d73027",
"Coldspot" = "#4575b4",
"Not sig" = "grey85"
)
tab <- table(sf_ok$cluster)
ann <- paste0(names(tab), ": ", as.integer(tab), collapse = " | ")
if (is.null(title)) title <- var
subtitle_txt <- if (add_subtitle) {
paste0(ann)
} else {
NULL
}
p_map <- ggplot2::ggplot(sf_ok) +
ggplot2::geom_sf(ggplot2::aes(fill = cluster),
color = "white", linewidth = 0.25) +
ggplot2::scale_fill_manual(values = pal, drop = FALSE) +
ggplot2::labs(
title = title,
subtitle = subtitle_txt,
fill = "Cluster"
) +
ggplot2::theme_minimal(base_size = 12) +
ggplot2::theme(
plot.title = ggplot2::element_text(face = "bold"),
legend.position = "right",
panel.grid = ggplot2::element_blank()
)
# Hilangkan sumbu koordinat kalau mau “clean map”
if (!show_axes) {
p_map <- p_map + ggplot2::theme_void(base_size = 12) +
ggplot2::theme(
plot.title = ggplot2::element_text(face = "bold"),
legend.position = "right"
)
# theme_void ngilangin subtitle juga? enggak, subtitle tetap ada.
# Tapi kalau kamu ingin subtitle tetap, theme_void aman.
}
list(data = sf_ok, plot = p_map)
}
GI_DBD <- getis_ord_map2(shp_sumut_join, "DBD_IR", lw_queen, title="DBD_IR")$plot
GI_Banjir <- getis_ord_map2(shp_sumut_join, "Banjir", lw_queen, title="Banjir")$plot
GI_Cuaca <- getis_ord_map2(shp_sumut_join, "Cuaca_Ekstreme", lw_queen, title="Cuaca_Ekstreme")$plot
GI_Kepadatan <- getis_ord_map2(shp_sumut_join, "Kepadatan_Penduduk", lw_queen, title="Kepadatan_Penduduk")$plot
GI_Miskin <- getis_ord_map2(shp_sumut_join, "Persentase_Miskin", lw_queen, title="Persentase_Miskin")$plot
GI_Sanitasi <- getis_ord_map2(shp_sumut_join, "Persentase_Sanitasi", lw_queen, title="Persentase_Sanitasi")$plot
fix_panel <- function(p){
p +
ggplot2::theme(
plot.title = ggplot2::element_text(size = 12, face = "bold"),
plot.subtitle = ggplot2::element_text(size = 9),
legend.title = ggplot2::element_text(size = 10),
legend.text = ggplot2::element_text(size = 9),
plot.margin = ggplot2::margin(2, 2, 2, 2)
)
}
GI_DBD2 <- fix_panel(GI_DBD)
GI_Banjir2 <- fix_panel(GI_Banjir)
GI_Cuaca2 <- fix_panel(GI_Cuaca)
GI_Kepadatan2 <- fix_panel(GI_Kepadatan)
GI_Miskin2 <- fix_panel(GI_Miskin)
GI_Sanitasi2 <- fix_panel(GI_Sanitasi)
p6_gi <- gridExtra::arrangeGrob(
GI_DBD2, GI_Banjir2, GI_Cuaca2,
GI_Kepadatan2, GI_Miskin2, GI_Sanitasi2,
ncol = 3
)
grid::grid.draw(p6_gi)
ggplot2::ggsave("GI_6panel.png", p6_gi,
width = 16, height = 9, units = "in", dpi = 300, bg = "white")
### Getis Ord Leaflet
getis_ord_leaflet <- function(sf_obj, var, lw,
alpha = 0.05, nsim = 999,
zero.policy = TRUE, seed = 123,
id_col = NULL, title = NULL){
stopifnot(var %in% names(sf_obj))
x <- sf_obj[[var]]
if (!is.numeric(x)) stop("Variabel harus numeric: ", var)
# pastikan leaflet pakai lon/lat (EPSG:4326)
sf_ll <- sf::st_transform(sf_obj, 4326)
# buang NA/Inf dan subset listw
ok <- is.finite(sf_ll[[var]])
sf_ok <- sf_ll[ok, ]
x_ok <- sf_ok[[var]]
lw_ok <- spdep::subset.listw(lw, subset = ok, zero.policy = zero.policy)
set.seed(seed)
# Hitung Gi* (permutation) -> menghasilkan Z-score
gi <- spdep::localG_perm(x_ok, lw_ok, nsim = nsim, zero.policy = zero.policy)
z <- as.numeric(gi)
p <- 2 * pnorm(abs(z), lower.tail = FALSE)
cls <- ifelse(p < alpha & z > 0, "Hotspot",
ifelse(p < alpha & z < 0, "Coldspot", "Not sig"))
sf_ok$GiZ <- z
sf_ok$pvalue <- p
sf_ok$cluster <- factor(cls, levels = c("Hotspot","Coldspot","Not sig"))
# label wilayah untuk tooltip/popup
if (!is.null(id_col) && id_col %in% names(sf_ok)) {
nm <- sf_ok[[id_col]]
} else {
# fallback: pakai row number
nm <- paste("Area", seq_len(nrow(sf_ok)))
}
sf_ok$._name <- nm
# palet warna
pal <- leaflet::colorFactor(
palette = c("Hotspot" = "#d73027", "Coldspot" = "#4575b4", "Not sig" = "#d9d9d9"),
domain = sf_ok$cluster
)
# tooltip & popup
fmt_num <- function(v) ifelse(is.na(v), "-", format(round(v, 3), nsmall = 3))
sf_ok$._tooltip <- paste0(
"<b>", sf_ok$._name, "</b><br/>",
var, ": ", sf_ok[[var]], "<br/>",
"Gi* Z: ", fmt_num(sf_ok$GiZ), "<br/>",
"p-value: ", fmt_num(sf_ok$pvalue), "<br/>",
"Cluster: <b>", sf_ok$cluster, "</b>"
)
sf_ok$._popup <- paste0(
"<div style='font-size:14px'>",
"<h4 style='margin:0 0 6px 0'>", sf_ok$._name, "</h4>",
"<b>Variabel:</b> ", var, " = ", sf_ok[[var]], "<br/>",
"<b>Getis-Ord Gi* (Z):</b> ", fmt_num(sf_ok$GiZ), "<br/>",
"<b>p-value:</b> ", fmt_num(sf_ok$pvalue), "<br/>",
"<b>Cluster:</b> ", sf_ok$cluster,
"</div>"
)
if (is.null(title)) title <- paste0("Getis–Ord Gi* (", var, ")")
m <- leaflet::leaflet(sf_ok, options = leaflet::leafletOptions(preferCanvas = TRUE)) |>
leaflet::addProviderTiles(leaflet::providers$CartoDB.Positron) |>
leaflet::addPolygons(
fillColor = ~pal(cluster),
fillOpacity = 0.75,
color = "white",
weight = 1,
opacity = 1,
smoothFactor = 0.2,
highlightOptions = leaflet::highlightOptions(weight = 3, color = "#333", bringToFront = TRUE),
label = lapply(sf_ok$._tooltip, htmltools::HTML),
popup = lapply(sf_ok$._popup, htmltools::HTML)
) |>
leaflet::addLegend(
position = "bottomright",
pal = pal,
values = sf_ok$cluster,
title = title,
opacity = 1
)
list(data = sf_ok, map = m)
}
res_gi_leaf <- getis_ord_leaflet(
sf_obj = shp_sumut_join,
var = "DBD_IR", # ganti sesuai nama kolom kamu
lw = lw_queen,
alpha = 0.05,
nsim = 999,
id_col = "KAB_KOTA" # ganti sesuai kolom nama kab/kota (opsional)
)
res_gi_leaf$map
# Uji Lagrange Multiplier
lm_tests <- lm.RStests(ols, listw = lw_queen, test = "all", zero.policy = TRUE)
lm_table <- data.frame(
`Statistik Uji` = c("LM Error (RSer)",
"LM Lag (RSlag)",
"Robust LM Error (AdjRSer)",
"Robust LM Lag (AdjRSlag)",
"LM SARMA (Gabungan)"),
Value = c(lm_tests$RSer$statistic,
lm_tests$RSlag$statistic,
lm_tests$adjRSer$statistic,
lm_tests$adjRSlag$statistic,
lm_tests$SARMA$statistic),
`p-value` = c(lm_tests$RSer$p.value,
lm_tests$RSlag$p.value,
lm_tests$adjRSer$p.value,
lm_tests$adjRSlag$p.value,
lm_tests$SARMA$p.value)
)
alpha <- 0.05
lm_table$Keputusan <- ifelse(lm_table$`p.value` < alpha, "Signifikan", "Tidak Signifikan")
names(lm_table)
str(lm_table)
# pembulatan biar mirip contoh
lm_table$Value <- round(lm_table$Value, 4)
lm_table$`p.value` <- round(lm_table$`p.value`, 4)
kable(lm_table, align = "lccc", caption = "Hasil Uji Lagrange Multiplier") |>
kable_styling(full_width = TRUE, bootstrap_options = c("striped")) |>
row_spec(0, bold = TRUE)
# Estimasi Model
sar <- lagsarlm(
DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi,
data = shp_sumut_join,
listw = lw_queen,
zero.policy = TRUE
)
summary(sar)
sdm <- lagsarlm(
DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi,
data = shp_sumut_join,
listw = lw_queen,
Durbin = TRUE, # aktifkan Spatial Lag of X
method = "eigen", # metode umum untuk estimasi SDM
zero.policy = TRUE
)
summary(sdm)
# Pemilihan Model Terbaik
make_model_compare_table <- function(models, model_names,
x_vars,
digits = 3,
p_digits = 3,
star = TRUE) {
stopifnot(length(models) == length(model_names))
# helper format
fnum <- function(x) ifelse(is.na(x), ".", formatC(x, digits = digits, format = "f"))
fp <- function(p) {
ifelse(is.na(p), ".",
formatC(p, digits = p_digits, format = "g"))
}
fp_star <- function(p){
if (is.na(p)) return(".")
txt <- fp(p)
if (!star) return(txt)
if (p < 0.001) return(paste0(txt, "***"))
if (p < 0.01) return(paste0(txt, "**"))
if (p < 0.05) return(paste0(txt, "*"))
txt
}
# ambil p-value koef dari summary.sarlm
get_pvals <- function(m){
s <- summary(m)
coefs <- as.data.frame(s$Coef)
pcol <- intersect(colnames(coefs), c("Pr(>|z|)", "Pr(>|t|)", "p-value"))
if (length(pcol) == 0) pcol <- tail(colnames(coefs), 1)
pvals <- coefs[[pcol[1]]]
names(pvals) <- rownames(coefs)
pvals
}
rows <- lapply(seq_along(models), function(i){
m <- models[[i]]
nm <- model_names[[i]]
aic <- AIC(m)
bic <- BIC(m)
# p uji rho (Wald z approx)
s <- summary(m)
p_rho <- NA_real_
if (!is.null(s$rho) && !is.null(s$rho.se)) {
z <- as.numeric(s$rho) / as.numeric(s$rho.se)
p_rho <- 2 * pnorm(abs(z), lower.tail = FALSE)
}
pv <- get_pvals(m)
# isi row (TANPA lambda)
out <- c(
Model = nm,
AIC = fnum(aic),
BIC = fnum(bic),
`p uji ρ` = fp_star(p_rho)
)
# p untuk X1..Xk (sesuai urutan x_vars)
for (j in seq_along(x_vars)){
v <- x_vars[j]
p_here <- pv[v]
if (is.null(p_here) || length(p_here) == 0 || is.na(p_here)) p_here <- NA_real_
out[paste0("p X", j)] <- fp_star(as.numeric(p_here))
}
out
})
df <- do.call(rbind, rows)
df <- as.data.frame(df, stringsAsFactors = FALSE, check.names = FALSE)
rownames(df) <- NULL
df
}
tbl_compare <- make_model_compare_table(
models = list(sar, sdm),
model_names = c("SAR", "SDM"),
x_vars = c("Banjir","Cuaca_Ekstreme","Kepadatan_Penduduk","Persentase_Miskin","Persentase_Sanitasi")
)
knitr::kable(
tbl_compare,
format = "html",
escape = FALSE,
caption = "Tabel. Perbandingan Model Spasial",
align = "c"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
)
# Uji Asumsi Klasik SAR (Model Terbaik)
uji_asumsi_klasik_SAR_kable <- function(model_sar,
lw,
reset_power = 2:3,
alpha = 0.05,
caption = "Hasil Uji Asumsi Klasik Model SAR",
digits = 5) {
# ========== VALIDASI ==========
if (inherits(model_sar, "summary.sarlm")) {
model_sar <- model_sar$object
}
if (!inherits(model_sar, "Sarlm"))
stop("model_sar harus objek SAR hasil lagsarlm().")
if (!inherits(lw, "listw"))
stop("lw harus objek listw.")
if (!requireNamespace("lmtest", quietly = TRUE)) stop("Paket 'lmtest' belum terpasang.")
if (!requireNamespace("knitr", quietly = TRUE)) stop("Paket 'knitr' belum terpasang.")
if (!requireNamespace("kableExtra", quietly = TRUE)) stop("Paket 'kableExtra' belum terpasang.")
# ========== RESIDUAL & FITTED (KHUSUS SAR) ==========
res_sar <- as.numeric(residuals(model_sar))
fit_sar <- as.numeric(fitted(model_sar))
fmt_num <- function(x) formatC(x, digits = digits, format = "f")
fmt_p <- function(p) formatC(p, digits = digits, format = "f")
# ========== 1) RESET (auxiliary – SAR) ==========
aux_reset <- data.frame(res_sar = res_sar, fit_sar = fit_sar)
for (p in reset_power) {
aux_reset[[paste0("fit_sar", p)]] <- fit_sar^p
}
m_restricted <- lm(res_sar ~ fit_sar, data = aux_reset)
m_full <- lm(
as.formula(
paste(
"res_sar ~",
paste(c("fit_sar", paste0("fit_sar", reset_power)), collapse = " + ")
)
),
data = aux_reset
)
rst <- lmtest::waldtest(m_restricted, m_full)
rst_stat <- unname(as.numeric(rst$F[2]))
rst_p <- unname(as.numeric(rst$`Pr(>F)`[2]))
kes_reset <- if (rst_p > alpha) {
"Hubungan variabel independen terhadap variabel dependen membentuk model linier"
} else {
"Indikasi ketidaklinieran / spesifikasi model belum tepat"
}
# ========== 2) NORMALITAS (Shapiro–Wilk) ==========
shp <- shapiro.test(res_sar)
shp_stat <- unname(as.numeric(shp$statistic[["W"]]))
shp_p <- shp$p.value
kes_shp <- if (shp_p > alpha) {
"Residual SAR berdistribusi normal"
} else {
"Residual SAR tidak berdistribusi normal"
}
# ========== 3) HOMOSKEDASTISITAS (Breusch–Pagan) ==========
aux_bp <- data.frame(e2_sar = res_sar^2, fit_sar = fit_sar)
m_bp <- lm(e2_sar ~ fit_sar, data = aux_bp)
bp <- lmtest::bptest(m_bp)
bp_stat <- unname(as.numeric(bp$statistic))
bp_p <- bp$p.value
kes_bp <- if (bp_p > alpha) {
"Tidak terdapat gejala heteroskedastisitas pada residual SAR"
} else {
"Terdapat gejala heteroskedastisitas pada residual SAR"
}
# ========== TABEL ==========
tbl <- data.frame(
`Jenis Uji` = c(
"Linearitas (Ramsey RESET)",
"Normalitas (Shapiro Wilk)",
"Homoskedastisitas (Breusch–Pagan)"
),
Statistik = c("RESET", "W", "BP"),
`Statistik Hitung` = c(
fmt_num(rst_stat),
fmt_num(shp_stat),
fmt_num(bp_stat)
),
`p-value` = c(
fmt_p(rst_p),
fmt_p(shp_p),
fmt_p(bp_p)
),
Kesimpulan = c(kes_reset, kes_shp, kes_bp),
check.names = FALSE
)
row.names(tbl) <- NULL
# --- OUTPUT KABLE (STYLE TETAP) ---
knitr::kable(
tbl,
format = "html",
escape = FALSE,
caption = caption,
align = c("l", "c", "c", "c", "l")
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
kableExtra::column_spec(1, width = "28%") |>
kableExtra::column_spec(2, width = "10%") |>
kableExtra::column_spec(3, width = "14%") |>
kableExtra::column_spec(4, width = "10%") |>
kableExtra::column_spec(5, width = "38%")
}
uji_asumsi_klasik_SAR_kable(
model_sar = sar,
lw = lw_queen,
caption = "Uji Asumsi Klasik Model SAR (Queen Contiguity)"
)
# Moran's I SAR
moran_residual_sar_kable <- function(model_sar,
lw,
nsim = 999,
alpha = 0.05,
seed = 123,
caption = "Moran's I Residual Model SAR (Monte Carlo)",
digits = 4) {
# validasi
if (inherits(model_sar, "summary.sarlm")) model_sar <- model_sar$object
if (!inherits(model_sar, "Sarlm")) stop("model_sar harus objek SAR hasil lagsarlm().")
if (!inherits(lw, "listw")) stop("lw harus objek listw.")
if (!requireNamespace("spdep", quietly = TRUE)) stop("Paket 'spdep' belum terpasang.")
if (!requireNamespace("knitr", quietly = TRUE)) stop("Paket 'knitr' belum terpasang.")
if (!requireNamespace("kableExtra", quietly = TRUE)) stop("Paket 'kableExtra' belum terpasang.")
# residual SAR
res_sar <- as.numeric(residuals(model_sar))
# Moran MC
set.seed(seed)
mor_mc <- spdep::moran.mc(res_sar, lw, nsim = nsim, zero.policy = TRUE)
I_val <- unname(as.numeric(mor_mc$statistic))
p_val <- mor_mc$p.value
ket <- if (p_val > alpha) {
"Tidak terdapat autokorelasi spasial residual"
} else {
"Masih terdapat autokorelasi spasial residual"
}
# tabel (1 baris)
tbl <- data.frame(
`I` = formatC(I_val, digits = digits, format = "f"),
`p (Monte Carlo)` = formatC(p_val, digits = digits, format = "f"),
`Keterangan` = ket,
check.names = FALSE
)
row.names(tbl) <- NULL
# kable style (ikuti style kamu)
knitr::kable(
tbl,
format = "html",
escape = FALSE,
caption = caption,
align = c("c", "c", "l")
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
) |>
kableExtra::column_spec(1, width = "20%") |>
kableExtra::column_spec(2, width = "35%") |>
kableExtra::column_spec(3, width = "45%")
}
moran_residual_sar_kable(
model_sar = sar,
lw = lw_queen,
nsim = 999,
alpha = 0.05,
caption = "Moran's I Residual Model SAR (Queen, Monte Carlo)"
)
# Uji Sensitivitas Bobot
# Queen Contiguity
make_queen_listw <- function(sf_obj, queen = TRUE, style = "W", zero.policy = TRUE){
nb <- spdep::poly2nb(sf_obj, queen = queen)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_queen <- make_queen_listw(shp_sumut_join)
lw_queen <- w_queen$lw
# KNN Contiguity
make_knn_listw <- function(sf_obj, k = 5, style = "W", zero.policy = TRUE){
pts <- sf::st_point_on_surface(sf_obj) # lebih aman dari centroid yang bisa keluar polygon
coords <- sf::st_coordinates(pts)
knn <- spdep::knearneigh(coords, k = k)
nb <- spdep::knn2nb(knn)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_knn5 <- make_knn_listw(shp_sumut_join, k = 5)
lw_knn5 <- w_knn5$lw
# Rook Contiguity
make_rook_listw <- function(sf_obj, style = "W", zero.policy = TRUE){
nb <- spdep::poly2nb(sf_obj, queen = FALSE)
lw <- spdep::nb2listw(nb, style = style, zero.policy = zero.policy)
list(nb = nb, lw = lw)
}
w_rook <- make_rook_listw(shp_sumut_join)
lw_rook <- w_rook$lw
# --- 2. Jalankan model SAR untuk setiap W ---
sar_queen <- lagsarlm(DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi,
data = shp_sumut_join, listw = lw_queen)
sar_rook <- lagsarlm(DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi,
data = shp_sumut_join, listw = lw_rook)
sar_knn5 <- lagsarlm(DBD_IR ~ Banjir + Cuaca_Ekstreme + Kepadatan_Penduduk + Persentase_Miskin + Persentase_Sanitasi,
data = shp_sumut_join, listw = lw_knn5)
# --- 3. Bandingkan hasil utama ---
data.frame(
Bobot = c("Queen", "Rook", "KNN-5"),
Rho = c(sar_queen$rho, sar_rook$rho, sar_knn5$rho),
AIC = c(AIC(sar_queen), AIC(sar_rook), AIC(sar_knn5)),
LogLik = c(logLik(sar_queen), logLik(sar_rook), logLik(sar_knn5))
)
df_sar <- data.frame(
Bobot = c("Queen", "Rook", "KNN-5"),
Rho = c(sar_queen$rho, sar_rook$rho, sar_knn5$rho),
AIC = c(AIC(sar_queen), AIC(sar_rook), AIC(sar_knn5)),
LogLik = c(logLik(sar_queen), logLik(sar_rook), logLik(sar_knn5))
)
df_sar <- transform(
df_sar,
Rho = round(Rho, 3),
AIC = round(AIC, 2),
LogLik = round(as.numeric(LogLik), 2)
)
df_sar |>
knitr::kable(
caption = "Perbandingan Model SAR berdasarkan Definisi Bobot Spasial",
align = "lccc",
col.names = c("Bobot Spasial", expression(rho), "AIC", "Log-Likelihood"),
format = "html",
escape = FALSE
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE,
position = "center"
) |>
kableExtra::row_spec(0, bold = TRUE) |>
kableExtra::column_spec(1, bold = TRUE)
# IDW
# 1) CEK DATA & KOLOM TARGET ----------------------------------
target_var <- "DBD_IR" # <- ganti kalau nama kolom berbeda
if (!exists("shp_sumut_join")) {
stop("Objek 'shp_sumut_join' belum ada. Pastikan kamu sudah join shapefile + data.")
}
if (!target_var %in% names(shp_sumut_join)) {
stop(paste0("Kolom '", target_var, "' tidak ditemukan. Cek: names(shp_sumut_join)"))
}
# 2) SIAPKAN POLYGON (VALID + HAPUS GEOMETRI KOSONG) -----------
shp_int <- shp_sumut_join |>
st_make_valid() |>
dplyr::filter(!sf::st_is_empty(geometry))
# 3) PROYEKSI KE UTM (METER) -----------------------------------
# Sumut umumnya cocok UTM Zone 47N: EPSG 32647
epsg_utm <- 32647
shp_int_utm <- st_transform(shp_int, epsg_utm)
# 4) POLYGON -> TITIK OBSERVASI (POINT ON SURFACE) -------------
pts_int <- shp_int_utm |>
st_point_on_surface() |>
dplyr::select(dplyr::all_of(target_var)) |>
dplyr::filter(!is.na(.data[[target_var]]))
# sf -> sp (untuk gstat)
pts_sp <- as(pts_int, "Spatial")
# 5) BUAT GRID PREDIKSI (CEPAT) --------------------------------
# Resolusi grid (meter). 5000 = 5 km (baseline cepat).
grid_size_m <- 5000
grid_sf <- st_make_grid(
st_as_sf(pts_int), # bbox titik -> cepat
cellsize = grid_size_m,
what = "centers"
) |>
st_sf(crs = st_crs(shp_int_utm))
grid_sp <- as(grid_sf, "Spatial")
# 6) IDW BASELINE ----------------------------------------------
p_idw <- 2 # power
k_idw <- 12 # jumlah tetangga terdekat
idw_res <- gstat::idw(
formula = as.formula(paste0(target_var, " ~ 1")),
locations = pts_sp,
newdata = grid_sp,
idp = p_idw,
nmax = k_idw
)
# 7) HASIL IDW (POINT) -> RASTER -> CLIP KE WILAYAH SUMUT -------
idw_ras <- rasterFromXYZ(
as.data.frame(idw_res)[, c("coords.x1", "coords.x2", "var1.pred")]
)
# set CRS raster agar sama dengan shp_int_utm
crs(idw_ras) <- CRS(st_crs(shp_int_utm)$wkt)
# clip pakai mask
sumut_sp <- as(shp_int_utm, "Spatial")
idw_clip <- raster::mask(idw_ras, sumut_sp)
# 8) SIAPKAN DATAFRAME UNTUK PLOT ------------------------------
# Raster -> data frame (permukaan)
idw_df <- as.data.frame(idw_clip, xy = TRUE, na.rm = TRUE)
colnames(idw_df) <- c("x", "y", "pred")
# Titik observasi (untuk plotting)
pts_df <- sf::st_coordinates(pts_int) |> as.data.frame()
colnames(pts_df) <- c("x", "y")
## Dengan Contour
p_idw_contour_noborder <- ggplot() +
geom_raster(data = idw_df, aes(x = x, y = y, fill = pred)) +
# garis kontur (lengkung-lengkung)
geom_contour(
data = idw_df,
aes(x = x, y = y, z = pred),
bins = 20,
color = "white",
linewidth = 0.45
) +
# titik observasi (opsional)
geom_point(data = pts_df, aes(x = x, y = y), size = 1.4) +
scale_fill_viridis_c(
name = "DBD IR (pred)",
option = "C"
) +
coord_sf(crs = sf::st_crs(shp_int_utm)) +
labs(
title = "Contour Interpolasi IDW (DBD IR)",
x = NULL,
y = NULL
) +
theme_minimal() +
theme(
axis.title = element_blank(), # <<< ini kuncinya
plot.title = element_text(face = "bold", hjust = 0.5)
)
p_idw_contour_noborder
ggsave("Permukaan_IDW_DBD_IR_Sumut.png", p_idw_contour_noborder, width = 7, height = 6, dpi = 300)
## Tanpa Contour
p_idw_grid <- ggplot() +
geom_raster(data = idw_df, aes(x = x, y = y, fill = pred)) +
scale_fill_viridis_c(
name = "DBD IR (pred)",
option = "C"
) +
coord_sf(
crs = sf::st_crs(shp_int_utm),
datum = sf::st_crs(4326) # supaya koordinat tampil derajat (E, N)
) +
labs(
title = "Interpolasi IDW (DBD IR)",
subtitle = paste0("p = ", p_idw, ", k = ", k_idw, ", grid = ", grid_size_m/1000, " km"),
x = NULL,
y = NULL
) +
theme_minimal() +
theme(
# grid belakang
panel.grid.major = element_line(color = "grey80", linewidth = 0.3),
panel.grid.minor = element_line(color = "grey90", linewidth = 0.2),
# hilangkan judul sumbu tapi angka tetap ada
axis.title = element_blank(),
# styling teks
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
print(p_idw_grid)
# 9) SIMPAN OUTPUT ----------------------------------
ggsave("IDW_DBD_IR_Sumut.png", p_idw_grid, width = 7, height = 6, dpi = 300)
# Kriging
# sf -> sp
pts_sp <- as(pts_int, "Spatial")
# Variogram empiris
vgm_emp <- gstat::variogram(DBD_IR ~ 1, pts_sp)
# Coba fit model variogram (pilih salah satu model: "Sph", "Exp", "Gau")
vgm_fit <- gstat::fit.variogram(vgm_emp, model = gstat::vgm(model = "Sph"))
# Plot variogram
plot(vgm_emp, vgm_fit, main = "Variogram DBD_IR (Empiris + Model)")
## Tabel Parameter Kriging
vgm_df <- as.data.frame(vgm_fit) |>
dplyr::mutate(
Parameter = dplyr::case_when(
model == "Nug" ~ "Nugget",
TRUE ~ "Partial Sill"
)
) |>
dplyr::select(Parameter, model, psill, range) |>
dplyr::mutate(
range = ifelse(Parameter == "Nugget", NA, range)
)
vgm_df |>
knitr::kable(
caption = "Parameter Variogram DBD IR (Ordinary Kriging)",
col.names = c("Parameter", "Model", "Nilai", "Range (m)"),
digits = 3,
align = "c"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
position = "center"
)
# sf -> sp
pts_sp <- as(pts_int, "Spatial")
grid_sp <- as(grid_sf, "Spatial")
# Ordinary Kriging
krig_res <- gstat::krige(
formula = DBD_IR ~ 1,
locations = pts_sp,
newdata = grid_sp,
model = vgm_fit
)
# hasil kriging (sp)
krig_res
# point -> raster
krig_ras <- rasterFromXYZ(
as.data.frame(krig_res)[, c("coords.x1","coords.x2","var1.pred")]
)
# set CRS
crs(krig_ras) <- CRS(st_crs(shp_int_utm)$wkt)
# clip ke wilayah Sumut
krig_clip <- raster::mask(krig_ras, as(shp_int_utm, "Spatial"))
## Dengan Contour
# 1) sf -> sp
pts_sp <- as(pts_int, "Spatial")
grid_sp <- as(grid_sf, "Spatial")
# 2) Ordinary Kriging (pakai variogram ter-fit)
krig_res <- gstat::krige(
formula = DBD_IR ~ 1,
locations = pts_sp,
newdata = grid_sp,
model = vgm_fit
)
# 3) Kriging -> raster -> clip/mask ke Sumut
krig_ras <- raster::rasterFromXYZ(
as.data.frame(krig_res)[, c("coords.x1", "coords.x2", "var1.pred")]
)
raster::crs(krig_ras) <- sp::CRS(sf::st_crs(shp_int_utm)$wkt)
krig_clip <- raster::mask(krig_ras, as(shp_int_utm, "Spatial"))
# 4) Raster -> data frame
krig_df <- as.data.frame(krig_clip, xy = TRUE, na.rm = TRUE)
colnames(krig_df) <- c("x", "y", "z")
# (opsional) titik observasi
pts_df <- sf::st_coordinates(pts_int) |> as.data.frame()
colnames(pts_df) <- c("x", "y")
# 5) Plot: surface + contour + titik (STYLE sama IDW)
p_krig_contour <- ggplot() +
geom_raster(data = krig_df, aes(x = x, y = y, fill = z)) +
geom_contour(
data = krig_df,
aes(x = x, y = y, z = z),
bins = 20,
color = "white",
linewidth = 0.45
) +
geom_point(data = pts_df, aes(x = x, y = y), size = 1.4) +
scale_fill_viridis_c(name = "Z", option = "C") +
coord_sf(crs = sf::st_crs(shp_int_utm)) +
labs(
title = "Contour Kriging (DBD IR)",
subtitle = "Ordinary Kriging",
x = NULL, y = NULL
) +
theme_minimal() +
theme(
axis.title = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
# raster -> df
krig_df <- as.data.frame(krig_clip, xy = TRUE, na.rm = TRUE)
colnames(krig_df) <- c("x", "y", "z")
# titik observasi (opsional)
pts_df <- sf::st_coordinates(pts_int) |> as.data.frame()
colnames(pts_df) <- c("x", "y")
## Tanpa Contour
p_krig_surface <- ggplot() +
geom_raster(data = krig_df, aes(x = x, y = y, fill = z)) +
scale_fill_viridis_c(name = "DBD IR (pred)", option = "C") +
coord_sf(crs = sf::st_crs(shp_int_utm)) +
labs(
title = "Kriging (DBD IR)",
subtitle = "Ordinary Kriging",
x = NULL, y = NULL
) +
theme_minimal() +
theme(
axis.title = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
print(p_krig_surface)
# Interpretasi Hasil
## Pemodelan Ekonometrik
# ============================================================
# 1) Tambahkan nilai prediksi & residual ke data spasial
# ============================================================
shp_sar_map <- shp_sumut_join |>
mutate(
pred_sar = as.numeric(fitted(sar)),
resid_sar = as.numeric(residuals(sar)),
resid_z = as.numeric(scale(resid_sar))
)
# ============================================================
# 2) Theme peta bersih (tanpa sumbu X-Y)
# ============================================================
theme_map_clean <- theme_minimal() +
theme(
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "right"
)
# ============================================================
# PETA PREDIKSI & RESIDUAL SAR (LAYOUT KONSISTEN)
# ============================================================
# --- Peta Prediksi (SAR) ---
p_pred_clean <- ggplot(shp_sar_map) +
geom_sf(aes(fill = pred_sar), color = "white", linewidth = 0.25) +
scale_fill_viridis_c(
name = "Prediksi\nDBD IR",
option = "C"
) +
labs(title = "Peta Prediksi (SAR)") +
coord_sf() +
theme_map_clean +
theme(
legend.position = "right",
plot.margin = margin(t = 5, r = 20, b = 5, l = 5)
)
# --- Peta Residual Standar (SAR) ---
p_resid_clean <- ggplot(shp_sar_map) +
geom_sf(aes(fill = resid_z), color = "white", linewidth = 0.25) +
scale_fill_gradient2(
name = "Residual (z)",
midpoint = 0
) +
labs(
title = "Peta Residual Standar (SAR)",
caption = "• Coklat: terprediksi terlalu rendah (underfit)\n• Biru: terprediksi terlalu tinggi (overfit)"
) +
coord_sf() +
theme_map_clean +
theme(
legend.position = "right",
plot.caption = element_text(hjust = 0.5, size = 9),
plot.margin = margin(t = 5, r = 5, b = 5, l = 20)
)
# --- Gabungkan kiri–kanan (legend per peta, map besar) ---
final_map <- p_pred_clean + p_resid_clean +
plot_layout(ncol = 2)
final_map
## Pemodelan Interpolasi
# range global dari dua metode
fill_lim <- range(
idw_df$pred,
krig_df$z,
na.rm = TRUE
)
p_idw_clean <- ggplot() +
geom_raster(data = idw_df, aes(x = x, y = y, fill = pred)) +
scale_fill_viridis_c(
name = "DBD IR (pred)",
option = "C",
limits = fill_lim
) +
coord_sf(crs = sf::st_crs(shp_int_utm)) +
labs(
title = "Interpolasi IDW (DBD IR)",
subtitle = "p = 2, k = 12, grid = 5 km"
) +
theme_void() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
p_krig_clean <- ggplot() +
geom_raster(data = krig_df, aes(x = x, y = y, fill = z)) +
scale_fill_viridis_c(
name = "DBD IR (pred)",
option = "C",
limits = fill_lim
) +
coord_sf(crs = sf::st_crs(shp_int_utm)) +
labs(
title = "Kriging (DBD IR)",
subtitle = "Ordinary Kriging"
) +
theme_void() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
p_idw_clean <- p_idw_clean +
theme(
plot.margin = margin(t = 5, r = 25, b = 5, l = 5) # kanan diperlebar
)
p_krig_clean <- p_krig_clean +
theme(
plot.margin = margin(t = 5, r = 5, b = 5, l = 25) # kiri diperlebar
)
p_final <- p_idw_clean + p_krig_clean +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "right")
p_final