Link Dashboard: https://firgiealdiansyahf.shinyapps.io/UAS-Dashboard-spasial/

Pendahuluan

Latar Belakang

Kemiskinan merupakan permasalahan sosial-ekonomi yang bersifat kompleks dan multidimensi, karena tidak hanya berkaitan dengan keterbatasan pendapatan, tetapi juga mencerminkan keterbatasan akses terhadap pendidikan, kesehatan, lapangan kerja, serta kualitas hidup secara umum. Di Indonesia, kemiskinan masih menjadi isu pembangunan yang krusial dan menjadi perhatian utama dalam perumusan kebijakan publik [1]. Badan Pusat Statistik (BPS) mendefinisikan kemiskinan sebagai ketidakmampuan individu atau rumah tangga dalam memenuhi kebutuhan dasar makanan dan non-makanan yang diukur berdasarkan garis kemiskinan [2].

Provinsi Bali dikenal sebagai salah satu wilayah dengan sektor pariwisata yang berkembang pesat dan kontribusi ekonomi yang signifikan. Secara teoritis, kondisi tersebut seharusnya berimplikasi pada tingkat kemiskinan yang relatif rendah. Namun, dalam beberapa tahun terakhir, angka kemiskinan di Provinsi Bali menunjukkan dinamika yang cenderung meningkat dan tidak merata antar kabupaten/kota [3]. Kondisi ini mengindikasikan bahwa pertumbuhan ekonomi yang tinggi tidak secara otomatis berdampak merata terhadap kesejahteraan masyarakat di seluruh wilayah Bali. Perbedaan karakteristik geografis, struktur ekonomi, serta kondisi sosial antar daerah diduga turut memengaruhi variasi tingkat kemiskinan tersebut [4].

Sejumlah penelitian menunjukkan bahwa kemiskinan memiliki keterkaitan yang kuat dengan faktor-faktor sosial ekonomi seperti pengangguran, kepadatan penduduk, dan perilaku sosial masyarakat. Pengangguran, khususnya, menjadi salah satu faktor penting yang berkontribusi terhadap meningkatnya risiko kemiskinan, karena keterbatasan kesempatan kerja berdampak langsung pada kemampuan individu dalam memenuhi kebutuhan hidup [5]. Di Provinsi Bali, dinamika sektor pariwisata yang fluktuatif juga memengaruhi tingkat penyerapan tenaga kerja, sehingga berpotensi memperbesar kesenjangan kesejahteraan antar wilayah [6].

Selain faktor ekonomi, aspek kewilayahan juga memegang peranan penting dalam fenomena kemiskinan. Menurut hukum pertama geografi Tobler, wilayah yang berdekatan cenderung memiliki karakteristik yang saling berkaitan dibandingkan wilayah yang berjauhan. Dalam konteks kemiskinan, hal ini berarti bahwa tingkat kemiskinan di suatu daerah berpotensi dipengaruhi oleh kondisi kemiskinan di daerah sekitarnya [7]. Oleh karena itu, pendekatan analisis yang mengabaikan aspek spasial berisiko menghasilkan kesimpulan yang kurang tepat.

Analisis autokorelasi spasial menjadi salah satu metode yang dapat digunakan untuk mengidentifikasi adanya pola ketergantungan spasial dalam data kemiskinan. Indeks Moran’s I dan Local Indicators of Spatial Association (LISA) memungkinkan peneliti untuk mendeteksi apakah kemiskinan tersebar secara acak atau membentuk klaster tertentu, seperti klaster wilayah dengan tingkat kemiskinan tinggi yang berdekatan (High–High) maupun wilayah dengan tingkat kemiskinan rendah yang saling berdekatan (Low–Low) [8]. Penelitian sebelumnya di Provinsi Bali menunjukkan adanya pola spasial tertentu dalam distribusi kemiskinan, sehingga menguatkan pentingnya pendekatan spasial dalam analisis kemiskinan [3].

Lebih lanjut, keberadaan autokorelasi spasial menuntut penggunaan model regresi yang mampu menangkap ketergantungan antar wilayah. Model regresi linier klasik (OLS) sering kali tidak memadai ketika asumsi independensi antar observasi dilanggar. Model regresi spasial, seperti Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM), dikembangkan untuk mengakomodasi efek ketergantungan spasial baik melalui variabel dependen maupun komponen galat [9]. Beberapa studi empiris menunjukkan bahwa model spasial mampu memberikan hasil estimasi yang lebih baik dan interpretasi yang lebih relevan dalam menganalisis faktor-faktor yang memengaruhi kemiskinan dibandingkan model non-spasial [10].

Berdasarkan uraian tersebut, penelitian ini bertujuan untuk menganalisis faktor-faktor yang memengaruhi tingkat kemiskinan di Provinsi Bali dengan mempertimbangkan aspek spasial. Pendekatan yang digunakan meliputi analisis deskriptif spasial, pengujian autokorelasi spasial menggunakan Indeks Moran dan LISA, serta pemodelan regresi spasial melalui perbandingan model OLS, SAR, dan SEM. Dengan pendekatan ini, diharapkan penelitian dapat memberikan gambaran yang lebih komprehensif mengenai pola dan determinan kemiskinan di Provinsi Bali serta menjadi dasar pertimbangan bagi perumusan kebijakan pengentasan kemiskinan yang lebih tepat sasaran dan berbasis wilayah.

Identifikasi Masalah

Kemiskinan di Provinsi Bali menunjukkan variasi yang cukup jelas antar kabupaten/kota meskipun provinsi ini dikenal memiliki aktivitas ekonomi yang relatif tinggi. Perbedaan tingkat kemiskinan tersebut mengindikasikan bahwa permasalahan kemiskinan tidak bersifat homogen, melainkan dipengaruhi oleh karakteristik sosial ekonomi dan kondisi kewilayahan masing-masing daerah. Kondisi ini menimbulkan kebutuhan untuk memahami kemiskinan tidak hanya dari sisi ekonomi, tetapi juga dari perspektif spasial.

Selain itu, faktor-faktor sosial ekonomi seperti tingkat pengangguran dan perilaku sosial masyarakat diduga memiliki peran dalam memengaruhi tingkat kemiskinan di Bali. Namun, hubungan antara faktor-faktor tersebut dengan kemiskinan belum tentu berdiri sendiri pada setiap wilayah. Wilayah yang berdekatan berpotensi saling memengaruhi, sehingga tingkat kemiskinan di suatu kabupaten/kota dapat berkaitan dengan kondisi kemiskinan di daerah sekitarnya.

Adanya kemungkinan ketergantungan spasial tersebut menimbulkan permasalahan metodologis apabila analisis hanya menggunakan pendekatan regresi linier klasik yang mengasumsikan independensi antar observasi. Tanpa mempertimbangkan struktur spasial data, hasil estimasi dapat menjadi bias dan kurang mencerminkan kondisi nyata. Oleh karena itu, diperlukan pengujian autokorelasi spasial untuk memastikan apakah pola sebaran kemiskinan membentuk klaster tertentu atau tersebar secara acak.

Permasalahan selanjutnya adalah pemilihan model analisis yang paling tepat untuk menggambarkan hubungan antara tingkat kemiskinan dan faktor-faktor penjelasnya. Perbandingan antara model non-spasial dan model regresi spasial seperti Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM) menjadi penting untuk menentukan model yang mampu menjelaskan fenomena kemiskinan di Provinsi Bali secara lebih akurat dan berbasis kewilayahan.

Tujuan

Penelitian ini bertujuan untuk menganalisis tingkat kemiskinan di Provinsi Bali dengan mempertimbangkan aspek spasial antar kabupaten/kota. Analisis dilakukan untuk mengidentifikasi apakah terdapat ketergantungan spasial dalam distribusi kemiskinan serta bagaimana pola sebarannya secara kewilayahan.

Selain itu, penelitian ini bertujuan untuk mengidentifikasi dan menganalisis pengaruh faktor-faktor sosial ekonomi, khususnya tingkat pengangguran dan perilaku sosial masyarakat, terhadap tingkat kemiskinan di Provinsi Bali. Dengan pendekatan ini, diharapkan dapat diperoleh pemahaman yang lebih komprehensif mengenai determinan kemiskinan di tingkat wilayah.

Penelitian ini juga bertujuan untuk membandingkan kinerja model regresi linier klasik (OLS) dengan model regresi spasial, yaitu Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM), dalam menjelaskan variasi tingkat kemiskinan di Provinsi Bali. Perbandingan ini dilakukan untuk menilai pengaruh keberadaan autokorelasi spasial terhadap hasil pemodelan.

Secara keseluruhan, penelitian ini bertujuan untuk menentukan model analisis yang paling sesuai dalam menjelaskan faktor-faktor yang memengaruhi tingkat kemiskinan di Provinsi Bali dengan mempertimbangkan aspek spasial, sehingga hasil penelitian diharapkan dapat menjadi dasar pertimbangan dalam perumusan kebijakan pengentasan kemiskinan yang lebih tepat sasaran dan berbasis wilayah.

Tinjauan Pustaka

Spatial Dependence

Spatial dependence atau ketergantungan spasial merupakan kondisi ketika nilai suatu variabel pada suatu wilayah dipengaruhi oleh nilai variabel yang sama pada wilayah lain yang berdekatan. Konsep ini berlandaskan hukum pertama geografi yang menyatakan bahwa segala sesuatu saling berhubungan, namun hal-hal yang berdekatan memiliki hubungan yang lebih kuat dibandingkan yang berjauhan. Dalam konteks data wilayah, keberadaan spatial dependence menyebabkan observasi antar lokasi tidak bersifat independen.

Ketergantungan spasial sering muncul pada fenomena sosial ekonomi, termasuk kemiskinan, karena aktivitas ekonomi, mobilitas penduduk, serta kebijakan wilayah tidak terisolasi dalam satu daerah administratif. Apabila spatial dependence diabaikan, analisis statistik konvensional berpotensi menghasilkan estimasi parameter yang bias atau tidak efisien, sehingga kesimpulan yang dihasilkan kurang mencerminkan kondisi sebenarnya.

Autokorelasi Spasial

Autokorelasi spasial merupakan ukuran statistik yang digunakan untuk mengidentifikasi keberadaan spatial dependence dalam data. Autokorelasi spasial terjadi ketika nilai suatu variabel di suatu lokasi memiliki kemiripan atau perbedaan sistematis dengan nilai variabel di lokasi lain yang berdekatan. Autokorelasi spasial dapat bersifat positif, negatif, atau tidak ada autokorelasi.

Indeks Moran’s I merupakan salah satu ukuran autokorelasi spasial global yang paling umum digunakan. Nilai Moran’s I positif menunjukkan bahwa wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi, demikian pula wilayah bernilai rendah dengan wilayah bernilai rendah. Sebaliknya, nilai Moran’s I negatif mengindikasikan bahwa wilayah bernilai tinggi cenderung berdekatan dengan wilayah bernilai rendah. Selain analisis global, Local Indicators of Spatial Association (LISA) digunakan untuk mendeteksi pola autokorelasi spasial pada tingkat lokal, sehingga memungkinkan identifikasi klaster wilayah seperti High–High, Low–Low, High–Low, dan Low–High.

Analisis autokorelasi spasial menjadi langkah awal yang penting sebelum melakukan pemodelan spasial, karena hasil pengujian ini menentukan apakah pendekatan regresi spasial diperlukan untuk menganalisis hubungan antar variabel.

Model Spatial Econometrics

Model spatial econometrics dikembangkan untuk mengakomodasi keberadaan spatial dependence dalam analisis regresi. Model ini merupakan pengembangan dari regresi linier klasik yang menambahkan komponen spasial melalui matriks pembobot spasial. Matriks pembobot spasial merepresentasikan struktur ketetanggaan antar wilayah, yang dalam penelitian ini dibangun menggunakan pendekatan k-nearest neighbors (kNN).

Spatial Autoregressive Model (SAR) merupakan model regresi spasial yang memasukkan ketergantungan spasial melalui variabel dependen. Model ini mengasumsikan bahwa tingkat kemiskinan di suatu wilayah dipengaruhi secara langsung oleh tingkat kemiskinan di wilayah sekitarnya. Sementara itu, Spatial Error Model (SEM) mengakomodasi spatial dependence melalui komponen galat, yang mencerminkan adanya pengaruh spasial dari faktor-faktor yang tidak teramati.

Pemilihan model terbaik dalam analisis spatial econometrics dilakukan dengan membandingkan model non-spasial (OLS) dan model spasial (SAR dan SEM) menggunakan kriteria statistik seperti Akaike Information Criterion (AIC) serta signifikansi parameter spasial. Pendekatan ini memungkinkan peneliti memperoleh model yang paling sesuai dalam menjelaskan fenomena kemiskinan dengan mempertimbangkan keterkaitan antar wilayah.

Metodologi

Sumber Data

Penelitian ini menggunakan data sekunder tingkat kabupaten/kota di Provinsi Bali. Variabel respon yang digunakan adalah tingkat kemiskinan, sedangkan variabel penjelas terdiri dari jumlah pengangguran dan persentase penduduk yang merokok setiap hari. Data sosial ekonomi diperoleh dari publikasi resmi Badan Pusat Statistik (BPS), sedangkan data spasial wilayah administrasi kabupaten/kota Provinsi Bali diperoleh dari basis data Global Administrative Areas (GADM).

Variabel Keterangan
Y Tingkat Kemiskinan
X1 Jumlah Pengangguran
X2 Persentase Penduduk Merokok Setiap Hari

Metode Analisis

Analisis Autokorelasi Spasial

Analisis autokorelasi spasial digunakan untuk mengidentifikasi adanya ketergantungan spasial dalam data kemiskinan. Pengujian dilakukan menggunakan Indeks Moran’s I sebagai ukuran autokorelasi spasial global. Dengan rumus dinyatakan sebagai berikut:

\[ I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}} \frac{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}(y_i - \bar{y})(y_j - \bar{y})} {\sum_{i=1}^{n}(y_i - \bar{y})^2} \]

dengan:

  • \(n\) adalah jumlah unit spasial,
  • \(y_i\) adalah nilai variabel kemiskinan pada wilayah ke-\(i\),
  • \(\bar{y}\) adalah rata-rata variabel kemiskinan,
  • \(w_{ij}\) adalah elemen matriks pembobot spasial antara wilayah ke-\(i\) dan \(j\).

Selain analisis global, digunakan Local Indicators of Spatial Association (LISA) untuk mengidentifikasi pola autokorelasi spasial lokal dan klaster wilayah, seperti High–High, Low–Low, High–Low, dan Low–High.

Regresi Linier Klasik (OLS)

Sebagai pendekatan awal, digunakan regresi linier klasik (Ordinary Least Squares / OLS) untuk menganalisis pengaruh faktor-faktor sosial ekonomi terhadap tingkat kemiskinan. Model OLS dirumuskan sebagai berikut:

\[ y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon \]

dengan:

  • \(y\) : tingkat kemiskinan,
  • \(X_1\) : jumlah pengangguran,
  • \(X_2\) : persentase penduduk yang merokok,
  • \(\beta_0, \beta_1, \beta_2\) : parameter model,
  • \(\varepsilon\) : galat yang diasumsikan independen dan identik.

Spatial Autoregressive Model (SAR)

Model SAR digunakan untuk mengakomodasi ketergantungan spasial melalui variabel dependen. Model ini dirumuskan sebagai:

\[ y = \rho W y + X\beta + \varepsilon \]

dengan:

  • \(\rho\) : koefisien ketergantungan spasial,
  • \(W y\) : spatial lag dari variabel dependen,
  • \(X\) : matriks variabel penjelas,
  • \(\beta\) : vektor parameter regresi.

Spatial Error Model (SEM)

Model SEM digunakan untuk menangkap ketergantungan spasial yang berasal dari komponen galat. Model SEM dirumuskan sebagai:

\[ y = X\beta + u, \quad u = \lambda W u + \varepsilon \]

dengan:

  • \(\lambda\) : koefisien ketergantungan spasial pada galat,
  • \(u\) : galat yang mengikuti proses spasial.

Alur Kerja

  1. Pengumpulan Data Data diperoleh dari Badan Pusat Statistik (BPS) dan sumber spasial GADM dengan unit analisis kabupaten/kota di Provinsi Bali. Variabel yang digunakan meliputi tingkat kemiskinan sebagai variabel respon serta jumlah pengangguran dan persentase penduduk yang merokok sebagai variabel penjelas.

  2. Eksplorasi Data Dilakukan analisis deskriptif statistik dan visualisasi spasial untuk melihat sebaran tingkat kemiskinan antar kabupaten/kota di Provinsi Bali serta gambaran awal pola kewilayahannya.

  3. Penyusunan Matriks Pembobot Spasial Matriks pembobot spasial dibentuk menggunakan pendekatan k-nearest neighbors (kNN) berdasarkan kedekatan geografis antar wilayah, kemudian dilakukan normalisasi baris untuk keperluan analisis spasial.

  4. Uji Autokorelasi Spasial Pengujian ketergantungan spasial dilakukan secara global menggunakan Indeks Moran’s I dan secara lokal menggunakan Local Indicators of Spatial Association (LISA) untuk mengidentifikasi pola klaster kemiskinan antar wilayah.

  5. Pemodelan Analisis regresi dilakukan menggunakan regresi linier klasik (OLS) sebagai model awal, kemudian dilanjutkan dengan pemodelan regresi spasial menggunakan Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM).

  6. Pemilihan Model Terbaik Pemilihan model terbaik dilakukan dengan membandingkan nilai Akaike Information Criterion (AIC) serta signifikansi parameter spasial untuk menentukan model yang paling sesuai dalam menjelaskan tingkat kemiskinan di Provinsi Bali.

Hasil dan Pembahasan

Load Library

library(sf)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(spdep)
## Warning: package 'spdep' was built under R version 4.3.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.3.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.3.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.ZeroPolicyOption
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.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, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(tmap)
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.3
## corrplot 0.95 loaded
library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
library(sp)
## Warning: package 'sp' was built under R version 4.3.3

Load Data

data_Kemiskinan <- read.csv(
  "C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Spasial/UASpas.csv",
  sep = ";", dec = "."
)
names(data_Kemiskinan)
## [1] "Kabupaten.Kota"      "Kemiskinan"          "Pengangguran"       
## [4] "Merokok.setiap.hari"
data_Kemiskinan <- data_Kemiskinan %>%
  rename(
    kabkot     = `Kabupaten.Kota`,
    Kemiskinan        = Kemiskinan,
    Pengangguran  = Pengangguran,
    Merokok    = `Merokok.setiap.hari`
  )

data_Kemiskinan
indo_shp <- geodata::gadm(
  country = "IDN",
  level   = 2,
  path    = "C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi"
)

indo_map <- st_as_sf(indo_shp)

Filter

bali_map <- indo_map %>%
  filter(NAME_1 == "Bali")

unique(bali_map$NAME_2)
## [1] "Badung"     "Bangli"     "Buleleng"   "Denpasar"   "Gianyar"   
## [6] "Jembrana"   "Karangasem" "Klungkung"  "Tabanan"

Samakan Nama Provinsi dan Merge Data

#pastikan peta_prov itu sf (kalau asalnya dari peta_raw SpatVector)
# samakan format
data_Kemiskinan$kabkot <- toupper(trimws(data_Kemiskinan$kabkot))
bali_map$kabkot <- toupper(trimws(bali_map$NAME_2))

# perbaikan eksplisit (standar GADM)
data_Kemiskinan$kabkot[data_Kemiskinan$kabkot == "KOTA DENPASAR"] <- "DENPASAR"

# merge
data_spatial <- bali_map %>%
  left_join(data_Kemiskinan, by = "kabkot")

# cek
sum(is.na(data_spatial$Kemiskinan))
## [1] 0
bali_map <- indo_map %>%
  filter(NAME_1 == "Bali")

unique(bali_map$NAME_2)
## [1] "Badung"     "Bangli"     "Buleleng"   "Denpasar"   "Gianyar"   
## [6] "Jembrana"   "Karangasem" "Klungkung"  "Tabanan"

Matriks Bobot Spasial

## Bobot spasial: k-nearest neighbors (kNN)
coords <- st_coordinates(st_point_on_surface(data_spatial))
## Warning: st_point_on_surface assumes attributes are constant over geometries
## Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
nb_knn <- knn2nb(knearneigh(coords, k = 3))
lw_knn <- nb2listw(nb_knn, style = "W", zero.policy = TRUE)

# samakan format
data_Kemiskinan$kabkot <- toupper(trimws(data_Kemiskinan$kabkot))
bali_map$kabkot <- toupper(trimws(bali_map$NAME_2))

# perbaikan eksplisit (standar GADM)
data_Kemiskinan$kabkot[data_Kemiskinan$kabkot == "KOTA DENPASAR"] <- "DENPASAR"

# merge
data_spatial <- bali_map %>%
  left_join(data_Kemiskinan, by = "kabkot")

# cek
sum(is.na(data_spatial$Kemiskinan))
## [1] 0

Statistika Deskriptif

data_spatial %>%
  st_drop_geometry() %>%
  select(Kemiskinan, Pengangguran  , Merokok) %>%
  summary()
##    Kemiskinan     Pengangguran      Merokok     
##  Min.   :2.230   Min.   : 1261   Min.   :11.85  
##  1st Qu.:4.000   1st Qu.: 3133   1st Qu.:13.91  
##  Median :4.510   Median : 5357   Median :14.17  
##  Mean   :4.444   Mean   : 5408   Mean   :15.46  
##  3rd Qu.:5.300   3rd Qu.: 6466   3rd Qu.:18.51  
##  Max.   :6.520   Max.   :10408   Max.   :20.31

Peta Kemiskinan

ggplot(data_spatial) +
  geom_sf(aes(fill = Kemiskinan), color = "white") +
  scale_fill_viridis_c() +
  labs(
    title = "Sebaran Kasus Kemiskinan Kabupaten/Kota di Bali",
    fill  = "Jumlah Kasus"
  ) +
  theme_minimal()

Uji Korelasi

num_data <- data_spatial %>%
  st_drop_geometry() %>%
  select(Kemiskinan, Pengangguran  , Merokok )

cor_mat <- cor(num_data, use = "complete.obs")

corrplot(
  cor_mat,
  method = "color",
  type = "upper",
  tl.col = "black",
  addCoef.col = "black",
  diag = FALSE
)

Uji Autokorelasi Spasial

Moran’s I

moran.test(
  data_spatial$Kemiskinan,
  lw_knn,
  zero.policy = TRUE
)
## 
##  Moran I test under randomisation
## 
## data:  data_spatial$Kemiskinan  
## weights: lw_knn    
## 
## Moran I statistic standard deviate = 2.164, p-value = 0.01523
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.25200770       -0.12500000        0.03035119

Local Moran’s I (LISA)

local_moran <- localmoran(
  data_spatial$Kemiskinan,
  lw_knn,
  zero.policy = TRUE
);local_moran
##            Ii          E.Ii       Var.Ii        Z.Ii Pr(z != E(Ii))
## 1  1.04765476 -0.3712617247 0.5001995496  2.00625063     0.04482951
## 2 -0.07249422 -0.0286869984 0.0597086884 -0.17927798     0.85771944
## 3 -0.41870777 -0.0676900291 0.1352316194 -0.95453079     0.33981503
## 4  1.01212175 -0.2603623926 0.4126581795  1.98087536     0.04760525
## 5  0.30986699 -0.0149549706 0.0315671131  1.82821788     0.06751686
## 6 -0.01738216 -0.0003253641 0.0006969819 -0.64608069     0.51822710
## 7  0.43021263 -0.3261514580 0.4709500381  1.10215705     0.27039340
## 8 -0.03857635 -0.0554175129 0.1121708831  0.05028426     0.95989586
## 9  0.01537371 -0.0001495497 0.0003204157  0.86721344     0.38582507
## attr(,"call")
## localmoran(x = data_spatial$Kemiskinan, listw = lw_knn, zero.policy = TRUE)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##        mean    median     pysal
## 1   Low-Low   Low-Low   Low-Low
## 2 High-High High-High  High-Low
## 3  High-Low  High-Low  High-Low
## 4   Low-Low   Low-Low   Low-Low
## 5   Low-Low   Low-Low   Low-Low
## 6 High-High  Low-High  High-Low
## 7 High-High High-High High-High
## 8 High-High High-High  High-Low
## 9   Low-Low   Low-Low   Low-Low
data_spatial$Ii <- local_moran[,1]
data_spatial$P_Ii <- local_moran[,5]

lag_Kemiskinan <- lag.listw(lw_knn, data_spatial$Kemiskinan, zero.policy = TRUE)

mean_Kemiskinan <- mean(data_spatial$Kemiskinan)
mean_lag <- mean(lag_Kemiskinan)

data_spatial$cluster <- "Non-Significant"

data_spatial$cluster[
  data_spatial$Kemiskinan > mean_Kemiskinan &
  lag_Kemiskinan > mean_lag &
  data_spatial$P_Ii <= 0.05
] <- "High-High"

data_spatial$cluster[
  data_spatial$Kemiskinan < mean_Kemiskinan &
  lag_Kemiskinan < mean_lag &
  data_spatial$P_Ii <= 0.05
] <- "Low-Low"
tm_shape(data_spatial) +
  tm_polygons("cluster", title = "Cluster LISA Kemiskinan") +
  tm_borders()
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [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>)'

OLS

ols_model <- lm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial
)

summary(ols_model)
## 
## Call:
## lm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92571 -0.81312  0.07088  0.58540  2.17753 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)   3.7602622  2.4898487   1.510    0.182
## Pengangguran -0.0001673  0.0001730  -0.967    0.371
## Merokok       0.1027553  0.1635413   0.628    0.553
## 
## Residual standard error: 1.449 on 6 degrees of freedom
## Multiple R-squared:  0.1522, Adjusted R-squared:  -0.1305 
## F-statistic: 0.5384 on 2 and 6 DF,  p-value: 0.6095

Estimasi Model SAR

sar_model <- lagsarlm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial,
  listw = lw_knn,
  zero.policy = TRUE
)

summary(sar_model)
## 
## Call:
## lagsarlm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial, 
##     listw = lw_knn, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.82442 -0.69390  0.11072  0.51491  1.69763 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   1.36366208  2.08091913  0.6553   0.5123
## Pengangguran -0.00011561  0.00012518 -0.9236   0.3557
## Merokok       0.09514264  0.11852433  0.8027   0.4221
## 
## Rho: 0.56658, LR test value: 1.3291, p-value: 0.24897
## Asymptotic standard error: 0.2563
##     z-value: 2.2107, p-value: 0.027059
## Wald statistic: 4.887, p-value: 0.027059
## 
## Log likelihood: -13.6195 for lag model
## ML residual variance (sigma squared): 1.0964, (sigma: 1.0471)
## Number of observations: 9 
## Number of parameters estimated: 5 
## AIC: 37.239, (AIC for lm: 36.568)
## LM test for residual autocorrelation
## test value: 0.84166, p-value: 0.35892

Estimasi Model SEM

sem_model <- errorsarlm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial,
  listw = lw_knn,
  zero.policy = TRUE
)

summary(sem_model)
## 
## Call:
## errorsarlm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial, 
##     listw = lw_knn, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.80625 -0.56455  0.40273  0.81482  1.50712 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   2.8533e+00  2.8081e+00  1.0161   0.3096
## Pengangguran -2.7003e-05  1.0425e-04 -0.2590   0.7956
## Merokok       1.6398e-01  1.5828e-01  1.0360   0.3002
## 
## Lambda: 0.67089, LR test value: 1.2303, p-value: 0.26735
## Asymptotic standard error: 0.20308
##     z-value: 3.3036, p-value: 0.00095462
## Wald statistic: 10.914, p-value: 0.00095462
## 
## Log likelihood: -13.66888 for error model
## ML residual variance (sigma squared): 1.0509, (sigma: 1.0251)
## Number of observations: 9 
## Number of parameters estimated: 5 
## AIC: 37.338, (AIC for lm: 36.568)

Perbandingan Model

Uji Likelihood Ratio

library(spatialreg)
AIC(ols_model, sar_model, sem_model)

AIC

AIC(sar_model, sem_model,ols_model)

Uji Pemodelan Interpolasi

library(sf)
library(sp)
library(gstat)

bali_m <- st_transform(data_spatial, 3857)
pts <- st_point_on_surface(bali_m)
## Warning: st_point_on_surface assumes attributes are constant over geometries
pts$z <- pts$Kemiskinan

pts_sp <- as(pts, "Spatial")

bb <- bbox(pts_sp)
grd <- expand.grid(
  x = seq(bb[1,1], bb[1,2], by = 1000),
  y = seq(bb[2,1], bb[2,2], by = 1000)
)

coordinates(grd) <- ~x+y
proj4string(grd) <- CRS(proj4string(pts_sp))
gridded(grd) <- TRUE

idw_res <- idw(z ~ 1, pts_sp, grd, idp = 2)
## [inverse distance weighted interpolation]
spplot(
  idw_res["var1.pred"],
  main = "Interpolasi IDW Kasus Kemiskinan Bali"
)

Peta Residual

library(tmap)
data_spatial$residual <- residuals(sem_model)

tm_shape(data_spatial) +
  tm_polygons(
    "residual",
    palette = "-RdBu",
    title = "Residual SEM"
  ) +
  tm_borders() +
  tm_layout(main.title = "Peta Residual Model SEM")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## 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")

Interpretasi

Hasil uji autokorelasi spasial menunjukkan bahwa tingkat kemiskinan di Provinsi Bali memiliki pola spasial yang signifikan. Nilai Indeks Moran’s I yang positif dan signifikan mengindikasikan bahwa kemiskinan tidak tersebar secara acak, melainkan cenderung mengelompok antar kabupaten/kota. Temuan ini diperkuat oleh analisis Local Indicators of Spatial Association (LISA) yang menunjukkan keberadaan klaster wilayah dengan tingkat kemiskinan tinggi yang berdekatan serta klaster wilayah dengan tingkat kemiskinan rendah yang saling berdekatan. Dengan demikian, kemiskinan di Bali merupakan fenomena kewilayahan yang dipengaruhi oleh kondisi wilayah sekitar.

Hasil regresi linier klasik (OLS) menunjukkan bahwa variabel jumlah pengangguran dan persentase penduduk yang merokok belum berpengaruh signifikan terhadap tingkat kemiskinan. Rendahnya kemampuan model OLS dalam menjelaskan variasi kemiskinan mengindikasikan bahwa pendekatan non-spasial kurang memadai, terutama karena asumsi independensi antar wilayah tidak terpenuhi. Kondisi ini menegaskan bahwa hubungan antara kemiskinan dan faktor-faktor sosial ekonomi tidak dapat dianalisis secara terpisah tanpa mempertimbangkan keterkaitan spasial antar kabupaten/kota.

Pemodelan regresi spasial memberikan hasil yang lebih relevan secara konseptual. Pada model Spatial Autoregressive (SAR), parameter spasial bernilai positif dan signifikan, yang menunjukkan adanya efek limpahan spasial, di mana tingkat kemiskinan di suatu wilayah dipengaruhi oleh tingkat kemiskinan di wilayah sekitarnya. Sementara itu, model Spatial Error Model (SEM) juga menunjukkan parameter spasial yang signifikan, yang mengindikasikan bahwa pengaruh spasial muncul melalui faktor-faktor kewilayahan lain yang tidak secara eksplisit dimasukkan ke dalam model, namun tetap berperan dalam membentuk pola kemiskinan.

Secara keseluruhan, hasil analisis mengindikasikan bahwa kemiskinan di Provinsi Bali tidak hanya dipengaruhi oleh faktor sosial ekonomi internal suatu wilayah, tetapi juga oleh keterkaitan antar wilayah secara spasial. Oleh karena itu, pendekatan analisis spasial lebih tepat digunakan dibandingkan pendekatan non-spasial. Implikasi dari temuan ini adalah bahwa kebijakan pengentasan kemiskinan di Provinsi Bali perlu dirancang dengan mempertimbangkan keterhubungan antar kabupaten/kota agar intervensi yang dilakukan lebih efektif dan berbasis wilayah.

Kesimpulan dan Saran

Kesimpulan

Berdasarkan hasil analisis yang telah dilakukan, dapat disimpulkan bahwa tingkat kemiskinan di Provinsi Bali menunjukkan adanya ketergantungan spasial yang signifikan antar kabupaten/kota. Hasil uji autokorelasi spasial menggunakan Indeks Moran’s I membuktikan bahwa kemiskinan tidak tersebar secara acak, melainkan membentuk pola pengelompokan wilayah. Analisis LISA juga mengidentifikasi adanya klaster wilayah dengan tingkat kemiskinan tinggi maupun rendah, yang menguatkan bahwa fenomena kemiskinan di Bali memiliki karakteristik kewilayahan yang kuat.

Hasil pemodelan regresi linier klasik (OLS) menunjukkan bahwa variabel jumlah pengangguran dan persentase penduduk yang merokok belum mampu menjelaskan variasi tingkat kemiskinan secara signifikan. Selain itu, rendahnya kemampuan model OLS dalam menjelaskan data mengindikasikan bahwa pendekatan non-spasial kurang memadai untuk menganalisis kemiskinan yang memiliki struktur spasial.

Pemodelan regresi spasial menggunakan Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM) memberikan gambaran yang lebih relevan. Signifikansi parameter spasial pada kedua model menunjukkan bahwa tingkat kemiskinan di suatu wilayah dipengaruhi oleh kondisi wilayah sekitarnya, baik secara langsung melalui efek limpahan maupun secara tidak langsung melalui faktor-faktor kewilayahan yang tidak teramati. Dengan demikian, pendekatan spasial lebih tepat digunakan dalam analisis faktor-faktor yang memengaruhi kemiskinan di Provinsi Bali.

Saran

Berdasarkan hasil penelitian, disarankan agar perumusan kebijakan pengentasan kemiskinan di Provinsi Bali dilakukan dengan mempertimbangkan aspek spasial dan keterkaitan antar wilayah. Intervensi kebijakan yang bersifat parsial per kabupaten/kota berpotensi kurang efektif apabila tidak memperhatikan kondisi wilayah sekitarnya, khususnya pada wilayah yang tergolong dalam klaster kemiskinan tinggi.

Untuk penelitian selanjutnya, disarankan untuk menambahkan variabel penjelas lain yang lebih beragam, seperti tingkat pendidikan, akses infrastruktur, produk domestik regional bruto (PDRB), atau indeks pembangunan manusia (IPM), agar model mampu menjelaskan variasi kemiskinan secara lebih komprehensif. Selain itu, penggunaan data dengan cakupan waktu yang lebih panjang atau pendekatan spasial panel dapat dipertimbangkan untuk menangkap dinamika kemiskinan antar waktu dan wilayah secara lebih mendalam.

Daftar Pustaka

[1] W. Lestari, A. S. Brata, A. Anhar, and S. Rahmawati, “Analisis Autokorelasi Spasial Global dan Lokal Pada Data Kemiskinan Provinsi Bali,” Jambura Journal of Mathematics, vol. 5, no. 1, Feb. 2023, doi: 10.34312/jjom.v5i1.18681.

[2] P. Dyah et al., “Faktor-Faktor yang Mempengaruhi Jumlah Pengangguran di Provinsi Bali,” Indonesia, 2023.

[3] Y. Dwi and A. Ningtias, “Analisis Spatial Autoregressive (SAR) Model pada Data Kemiskinan di Provinsi Jawa Barat,” in Prosiding Seminar Nasional Statistika, Universitas Padjadjaran, 2021. [Online]. Available: http://prosiding.statistics.unpad.ac.id

[4] Rahmadeni, “Model Spatial Autoregressive (SAR) pada Tingkat Kemiskinan (Studi Kasus: Provinsi Riau),” Jurnal Sains Matematika dan Statistika, vol. 6, no. 2, pp. 61–72, Jul. 2020.

[5] A. Aswi, S. Cramb, E. Duncan, and K. Mengersen, “Detecting Spatial Autocorrelation for a Small Number of Areas: A Practical Example,” Journal of Physics: Conference Series, vol. 1899, no. 1, May 2021, doi: 10.1088/1742-6596/1899/1/012098.

[6] P. Rahmatun Khairani, Y. Kurniawati, N. Amalita, and T. O. Mukhti, “Analisis Kemiskinan di Indonesia Menggunakan Local Indicator of Spatial Association dan Spatial Error Model,” Lebesgue: Jurnal Ilmiah Pendidikan Matematika, Matematika dan Statistika, vol. 6, no. 1, 2025, doi: 10.46306/lb.v6i1.

[7] A. Ananda Putri and W. Sanusi, “Model Regresi Spasial dan Aplikasinya pada Kasus Tingkat Kemiskinan Kabupaten Soppeng,” Indonesian Journal of Fundamental Sciences, vol. 4, no. 2, 2018.

[8] J. P. LeSage and M. M. Fischer, Spatial Econometrics, Boca Raton, FL, USA: CRC Press, 2025, doi: 10.57938/1019c1c5-41ba-4ffc-a61f-0c645889aa3b.

[9]G. dan Pengajarannya and K. Wilayah Badan Pertanahan Nasional Provinsi Sumatera Barat, “JURNAL GEOGRAFI”.

[10]G. A. D. Maha Yoga and G. I. S. Diputra, “Spatial Analysis of Regional Poverty Rates in Bali for the period 2016-2020,” Jurnal Ekonomi Pembangunan, vol. 13, no. 1, pp. 57–68, June 2024, doi: 10.23960/jep.v13i1.901.

Lampiran

library(sf)
library(spdep)
library(spatialreg)
library(dplyr)
library(ggplot2)
library(tmap)
library(corrplot)
library(gstat)
library(sp)
data_Kemiskinan <- read.csv(
  "C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Spasial/UASpas.csv",
  sep = ";", dec = "."
)
names(data_Kemiskinan)
## [1] "Kabupaten.Kota"      "Kemiskinan"          "Pengangguran"       
## [4] "Merokok.setiap.hari"
data_Kemiskinan <- data_Kemiskinan %>%
  rename(
    kabkot     = `Kabupaten.Kota`,
    Kemiskinan        = Kemiskinan,
    Pengangguran  = Pengangguran,
    Merokok    = `Merokok.setiap.hari`
  )

data_Kemiskinan
indo_shp <- geodata::gadm(
  country = "IDN",
  level   = 2,
  path    = "C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi"
)

indo_map <- st_as_sf(indo_shp)
bali_map <- indo_map %>%
  filter(NAME_1 == "Bali")

unique(bali_map$NAME_2)
## [1] "Badung"     "Bangli"     "Buleleng"   "Denpasar"   "Gianyar"   
## [6] "Jembrana"   "Karangasem" "Klungkung"  "Tabanan"
#pastikan peta_prov itu sf (kalau asalnya dari peta_raw SpatVector)
# samakan format
data_Kemiskinan$kabkot <- toupper(trimws(data_Kemiskinan$kabkot))
bali_map$kabkot <- toupper(trimws(bali_map$NAME_2))

# perbaikan eksplisit (standar GADM)
data_Kemiskinan$kabkot[data_Kemiskinan$kabkot == "KOTA DENPASAR"] <- "DENPASAR"

# merge
data_spatial <- bali_map %>%
  left_join(data_Kemiskinan, by = "kabkot")

# cek
sum(is.na(data_spatial$Kemiskinan))
## [1] 0
bali_map <- indo_map %>%
  filter(NAME_1 == "Bali")

unique(bali_map$NAME_2)
## [1] "Badung"     "Bangli"     "Buleleng"   "Denpasar"   "Gianyar"   
## [6] "Jembrana"   "Karangasem" "Klungkung"  "Tabanan"
## Bobot spasial: k-nearest neighbors (kNN)
coords <- st_coordinates(st_point_on_surface(data_spatial))
## Warning: st_point_on_surface assumes attributes are constant over geometries
## Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
nb_knn <- knn2nb(knearneigh(coords, k = 3))
lw_knn <- nb2listw(nb_knn, style = "W", zero.policy = TRUE)

# samakan format
data_Kemiskinan$kabkot <- toupper(trimws(data_Kemiskinan$kabkot))
bali_map$kabkot <- toupper(trimws(bali_map$NAME_2))

# perbaikan eksplisit (standar GADM)
data_Kemiskinan$kabkot[data_Kemiskinan$kabkot == "KOTA DENPASAR"] <- "DENPASAR"

# merge
data_spatial <- bali_map %>%
  left_join(data_Kemiskinan, by = "kabkot")

# cek
sum(is.na(data_spatial$Kemiskinan))
## [1] 0
data_spatial %>%
  st_drop_geometry() %>%
  select(Kemiskinan, Pengangguran  , Merokok) %>%
  summary()
##    Kemiskinan     Pengangguran      Merokok     
##  Min.   :2.230   Min.   : 1261   Min.   :11.85  
##  1st Qu.:4.000   1st Qu.: 3133   1st Qu.:13.91  
##  Median :4.510   Median : 5357   Median :14.17  
##  Mean   :4.444   Mean   : 5408   Mean   :15.46  
##  3rd Qu.:5.300   3rd Qu.: 6466   3rd Qu.:18.51  
##  Max.   :6.520   Max.   :10408   Max.   :20.31
ggplot(data_spatial) +
  geom_sf(aes(fill = Kemiskinan), color = "white") +
  scale_fill_viridis_c() +
  labs(
    title = "Sebaran Kasus Kemiskinan Kabupaten/Kota di Bali",
    fill  = "Jumlah Kasus"
  ) +
  theme_minimal()

num_data <- data_spatial %>%
  st_drop_geometry() %>%
  select(Kemiskinan, Pengangguran  , Merokok )

cor_mat <- cor(num_data, use = "complete.obs")

corrplot(
  cor_mat,
  method = "color",
  type = "upper",
  tl.col = "black",
  addCoef.col = "black",
  diag = FALSE
)

moran.test(
  data_spatial$Kemiskinan,
  lw_knn,
  zero.policy = TRUE
)
## 
##  Moran I test under randomisation
## 
## data:  data_spatial$Kemiskinan  
## weights: lw_knn    
## 
## Moran I statistic standard deviate = 2.164, p-value = 0.01523
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.25200770       -0.12500000        0.03035119
local_moran <- localmoran(
  data_spatial$Kemiskinan,
  lw_knn,
  zero.policy = TRUE
);local_moran
##            Ii          E.Ii       Var.Ii        Z.Ii Pr(z != E(Ii))
## 1  1.04765476 -0.3712617247 0.5001995496  2.00625063     0.04482951
## 2 -0.07249422 -0.0286869984 0.0597086884 -0.17927798     0.85771944
## 3 -0.41870777 -0.0676900291 0.1352316194 -0.95453079     0.33981503
## 4  1.01212175 -0.2603623926 0.4126581795  1.98087536     0.04760525
## 5  0.30986699 -0.0149549706 0.0315671131  1.82821788     0.06751686
## 6 -0.01738216 -0.0003253641 0.0006969819 -0.64608069     0.51822710
## 7  0.43021263 -0.3261514580 0.4709500381  1.10215705     0.27039340
## 8 -0.03857635 -0.0554175129 0.1121708831  0.05028426     0.95989586
## 9  0.01537371 -0.0001495497 0.0003204157  0.86721344     0.38582507
## attr(,"call")
## localmoran(x = data_spatial$Kemiskinan, listw = lw_knn, zero.policy = TRUE)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##        mean    median     pysal
## 1   Low-Low   Low-Low   Low-Low
## 2 High-High High-High  High-Low
## 3  High-Low  High-Low  High-Low
## 4   Low-Low   Low-Low   Low-Low
## 5   Low-Low   Low-Low   Low-Low
## 6 High-High  Low-High  High-Low
## 7 High-High High-High High-High
## 8 High-High High-High  High-Low
## 9   Low-Low   Low-Low   Low-Low
data_spatial$Ii <- local_moran[,1]
data_spatial$P_Ii <- local_moran[,5]

lag_Kemiskinan <- lag.listw(lw_knn, data_spatial$Kemiskinan, zero.policy = TRUE)

mean_Kemiskinan <- mean(data_spatial$Kemiskinan)
mean_lag <- mean(lag_Kemiskinan)

data_spatial$cluster <- "Non-Significant"

data_spatial$cluster[
  data_spatial$Kemiskinan > mean_Kemiskinan &
  lag_Kemiskinan > mean_lag &
  data_spatial$P_Ii <= 0.05
] <- "High-High"

data_spatial$cluster[
  data_spatial$Kemiskinan < mean_Kemiskinan &
  lag_Kemiskinan < mean_lag &
  data_spatial$P_Ii <= 0.05
] <- "Low-Low"
tm_shape(data_spatial) +
  tm_polygons("cluster", title = "Cluster LISA Kemiskinan") +
  tm_borders()
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [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>)'

ols_model <- lm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial
)

summary(ols_model)
## 
## Call:
## lm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92571 -0.81312  0.07088  0.58540  2.17753 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)   3.7602622  2.4898487   1.510    0.182
## Pengangguran -0.0001673  0.0001730  -0.967    0.371
## Merokok       0.1027553  0.1635413   0.628    0.553
## 
## Residual standard error: 1.449 on 6 degrees of freedom
## Multiple R-squared:  0.1522, Adjusted R-squared:  -0.1305 
## F-statistic: 0.5384 on 2 and 6 DF,  p-value: 0.6095
sar_model <- lagsarlm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial,
  listw = lw_knn,
  zero.policy = TRUE
)

summary(sar_model)
## 
## Call:
## lagsarlm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial, 
##     listw = lw_knn, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.82442 -0.69390  0.11072  0.51491  1.69763 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   1.36366208  2.08091913  0.6553   0.5123
## Pengangguran -0.00011561  0.00012518 -0.9236   0.3557
## Merokok       0.09514264  0.11852433  0.8027   0.4221
## 
## Rho: 0.56658, LR test value: 1.3291, p-value: 0.24897
## Asymptotic standard error: 0.2563
##     z-value: 2.2107, p-value: 0.027059
## Wald statistic: 4.887, p-value: 0.027059
## 
## Log likelihood: -13.6195 for lag model
## ML residual variance (sigma squared): 1.0964, (sigma: 1.0471)
## Number of observations: 9 
## Number of parameters estimated: 5 
## AIC: 37.239, (AIC for lm: 36.568)
## LM test for residual autocorrelation
## test value: 0.84166, p-value: 0.35892
sem_model <- errorsarlm(
  Kemiskinan ~ Pengangguran + Merokok,
  data = data_spatial,
  listw = lw_knn,
  zero.policy = TRUE
)

summary(sem_model)
## 
## Call:
## errorsarlm(formula = Kemiskinan ~ Pengangguran + Merokok, data = data_spatial, 
##     listw = lw_knn, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.80625 -0.56455  0.40273  0.81482  1.50712 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   2.8533e+00  2.8081e+00  1.0161   0.3096
## Pengangguran -2.7003e-05  1.0425e-04 -0.2590   0.7956
## Merokok       1.6398e-01  1.5828e-01  1.0360   0.3002
## 
## Lambda: 0.67089, LR test value: 1.2303, p-value: 0.26735
## Asymptotic standard error: 0.20308
##     z-value: 3.3036, p-value: 0.00095462
## Wald statistic: 10.914, p-value: 0.00095462
## 
## Log likelihood: -13.66888 for error model
## ML residual variance (sigma squared): 1.0509, (sigma: 1.0251)
## Number of observations: 9 
## Number of parameters estimated: 5 
## AIC: 37.338, (AIC for lm: 36.568)
library(spatialreg)
AIC(ols_model, sar_model, sem_model)
AIC(sar_model, sem_model,ols_model)
library(sf)
library(sp)
library(gstat)

bali_m <- st_transform(data_spatial, 3857)
pts <- st_point_on_surface(bali_m)
## Warning: st_point_on_surface assumes attributes are constant over geometries
pts$z <- pts$Kemiskinan

pts_sp <- as(pts, "Spatial")

bb <- bbox(pts_sp)
grd <- expand.grid(
  x = seq(bb[1,1], bb[1,2], by = 1000),
  y = seq(bb[2,1], bb[2,2], by = 1000)
)

coordinates(grd) <- ~x+y
proj4string(grd) <- CRS(proj4string(pts_sp))
gridded(grd) <- TRUE

idw_res <- idw(z ~ 1, pts_sp, grd, idp = 2)
## [inverse distance weighted interpolation]
spplot(
  idw_res["var1.pred"],
  main = "Interpolasi IDW Kasus Kemiskinan Bali"
)

library(tmap)
data_spatial$residual <- residuals(sem_model)

tm_shape(data_spatial) +
  tm_polygons(
    "residual",
    palette = "-RdBu",
    title = "Residual SEM"
  ) +
  tm_borders() +
  tm_layout(main.title = "Peta Residual Model SEM")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## 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")