Regresi spasial adalah metode statistika yang secara khusus memasukkan aspek spasial ke dalam kerangka model regresi. Terdapat dua efek inti pada regresi spasial, yaitu keragaman spasial (spatial heterogeneity) atau ketergantungan spasial (spatial dependence).

1 Library

library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(spatialreg)
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library(performance)
library(rsconnect)

2 DATA

2.1 Import SHP

library(sf)
jawa <- st_read("D:/Jawamap/jawa.shp")
## Reading layer `jawa' from data source `D:\Jawamap\jawa.shp' using driver `ESRI Shapefile'
## Simple feature collection with 119 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84

2.2 Plot per Provinsi

Pulau Jawa Barat terdiri dari 6 Provinsi, yaitu Provinsi Banten, DKI Jakarta, Jawa Barat, Jawa Tengah, DI Yogyakarta dan Jawa Timur dengan jumlah keseluruhan 119 kabupaten/kota.

plot(jawa["PROVINSI"])

2.3 Import data amatan

Data yang digunakan diperoleh dari publikasi BPS. Peubah yang digunakan yaitu:

Peubah Keterangan Satuan
Y Persentase Penduduk Miskin (PPM) %
X1 Persentase pengguna BPJS Penerima Bantuan Iuran (BPJSPBI) %
X2 Persentase penduduk usia 15 tahun ke atas dengan status tidak bekerja (PPTB) %
X3 Persentase Rumah Tangga Penerima Bantuan Pangan Non Tunai (PRTPBPNT) %
X4 Persentase pengeluaran perkapita kelompok makanan (PPKM) %
X6 Indeks Pembangunan Manusia (IPM) % 
library(readxl)
data_jawa <- read_excel("D:/jawa.xlsx", col_names = TRUE)
head(data_jawa)
## # A tibble: 6 × 11
##    kode NAMA_KAB        TK BPJSPBI  PPTB PRTPBPNT  PPKM   IPM  long   lat ID2013
##   <dbl> <chr>        <dbl>   <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1 kepulauan S… 14.1     98.4  52.8     16.8  62.1  72.8  106. -5.80   3101
## 2     2 Jakarta Sel…  3.52    87.1  63.7     15.5  38.4  85.2  107. -6.27   3171
## 3     3 Jakarta Tim…  4.9     84.5  51.6      1.2  43.2  83.4  107. -6.26   3172
## 4     4 Jakarta Pus…  4.3     86.5  50.6      0    42.3  82.1  107. -6.18   3173
## 5     5 Jakarta Bar…  4.22    87.3  55.5      0    35.8  82.5  107. -6.17   3174
## 6     6 Jakarta Uta…  7.24    82.2  54.1     10.5  32.1  80.8  107. -6.13   3175

2.4 Menggabungkan SHP dan data amatan

library(dplyr)
## 
## 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
jawa$ID2013 <- as.character(jawa$ID2013) 
data_jawa$ID2013 <- as.character(data_jawa$ID2013)
# Join data
jawa <- jawa %>%
  left_join(data_jawa, by = c("ID2013" = "ID2013"))

3 Eksplorasi Data

3.1 Plot peta peubah respon

library(ggplot2)
library(ggspatial)
library(sf)

# Buat kelas kategori manual
jawa$Kategori <- cut(
  jawa$TK,
  breaks = c(-Inf, 5, 10, 15, 20, Inf),
  labels = c("<5% sangat rendah", "5% - 10% rendah", "10% - 15% sedang",
             "15% - 20% tinggi", ">20% sangat tinggi"),
  right = FALSE
)

ggplot(jawa) +
  geom_sf(aes(fill = Kategori), color = "grey40", linewidth = 0.2) +
  ggtitle("Persentase Penduduk Miskin") +
  scale_fill_manual(values = rev(hcl.colors(5, "ArmyRose"))) +
  annotation_north_arrow(
    location = "bl",       # bottom-left (kiri bawah)
    which_north = "true",
    pad_x = unit(0.5, "cm"),
    pad_y = unit(0.5, "cm"),
    style = north_arrow_fancy_orienteering
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold"),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    panel.grid.major = element_blank(),  # hilangkan grid
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  ) +
  labs(fill = "Kategori")

Terlihat bahwa Sampang memiliki persentase penduduk miskin tertinggi, dengan nilai persentase penduduk miskin lebih dari 20%, sedangkan Cilegon merupakan wilayah dengan persentase penduduk miskin terendah, yaitu kurang dari 5%. Kemudian terlihat terbentuknya kelompok berdasarkan nilai persentase kelompok miskin, contohnya Kabupaten Kebumen, Banjarnegara, Wonosobo, Purbalingga, pemalang, dan Brebes dengan kategori persentase penduduk miskin tinggi.

3.2 Plot peta peubah penjelas

library(ggplot2)
library(sf)
library(patchwork)

# X1
X1 <- ggplot(jawa) +
  geom_sf(aes(fill = BPJSPBI), color = "grey40", linewidth = 0.2) +
  ggtitle("Persentase pengguna BPJSPBI") +
  scale_fill_gradientn(colours = rev(hcl.colors(7, "Fall"))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(fill = "%")
X1

# X2
X2 <- ggplot(jawa) +
  geom_sf(aes(fill = PPTB), color = "grey40", linewidth = 0.2) +
  ggtitle("Persentase penduduk usia 15 tahun ke atas dengan status tidak bekerja") +
  scale_fill_gradientn(colours = rev(hcl.colors(7, "Geyser"))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(fill = "%")
X2

# X3
X3 <- ggplot(jawa) +
  geom_sf(aes(fill = PRTPBPNT), color = "grey40", linewidth = 0.2) +
  ggtitle("Persentase RT penerima BPNT") +
  scale_fill_gradientn(colours = rev(hcl.colors(7, "Temps"))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(fill = "%")
X3

# X4
X4 <- ggplot(jawa) +
  geom_sf(aes(fill = PPKM), color = "grey40", linewidth = 0.2) +
  ggtitle("Persentase pengeluaran per kapita kelompok makanan") +
  scale_fill_gradientn(colours = rev(hcl.colors(7, "TealRose"))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(fill = "%")
X4

# X5
X5 <- ggplot(jawa) +
  geom_sf(aes(fill = IPM), color = "grey40", linewidth = 0.2) +
  ggtitle("Indeks pembangunan Manusia") +
  scale_fill_gradientn(colours = rev(hcl.colors(7, "Tropic"))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(fill = "Indeks")
X5

3.3 Plot Peubah penjelas

library(ggplot2)
library(dplyr)
library(ggpubr)

df <- data.frame(
  Peubah = rep(c('X1', 'X2', 'X3', 'X4', 'X5'), each = 119),
  kabkot = data_jawa$kode,
  nilai_peubah = c(data_jawa$BPJSPBI, data_jawa$PPTB, data_jawa$PRTPBPNT, data_jawa$PPKM, data_jawa$IPM)
)

findoutlier <- function(x) {
  x < quantile(x, .25) - 1.5 * IQR(x) |
    x > quantile(x, .75) + 1.5 * IQR(x)
}

df <- df %>%
  group_by(Peubah) %>%
  mutate(outlier = ifelse(findoutlier(nilai_peubah), kabkot, NA)) %>%
  ungroup()

warna_peubah <- c(
  "X1"  = "#A7C6DA",
  "X2" = "#C0C7BD",
  "X3" = "#EAD3C6",
  "X4" = "#D7A1AA",
  "X5" = "#AC4D61"
)

plot_list <- lapply(unique(df$Peubah), function(v) {
  ggplot(df %>% filter(Peubah == v),
         aes(x = Peubah, y = nilai_peubah, fill = Peubah)) +
    geom_boxplot(color = "black", notch = TRUE) +
    geom_text(aes(label = outlier),
              na.rm = TRUE, hjust = -0.2,
              size = 3, color = "black") +
    theme_minimal() +
    theme(
      panel.border = element_rect(color = "black", fill = NA),
      panel.grid = element_blank(),
      legend.position = "none",
      text = element_text(size = 12, color = "black"),
      axis.title.x = element_blank(),
      axis.text.x = element_text(face = "bold")
    ) +
    scale_fill_manual(values = warna_peubah) +
    labs(y = paste("Nilai", v))
})

ggarrange(plotlist = plot_list, ncol = 3, nrow = 2)

3.4 Korelasi Antar Peubah

corrplot::corrplot(cor(data_jawa[,c(3:8)]), method = "color", type = "upper", tl.cex = 0.8,
  tl.col = "black",
  tl.srt = 45,
  addCoef.col = "black")

peubah PPTB dan IPM memiliki hubungan negatif dengan persentase penduduk miskin, sedangkan BPJSPBI, PRTPBPNT, dan PPKM memiliki hubungan positif dengan persentase penduduk miskin.

4 Matriks pembobot spasial

library(sf)
library(spdep)

# Perbaikan geometri shapefile
jawa <- st_make_valid(jawa)
jawa <- st_collection_extract(jawa, "POLYGON") # memastikan poligon
## Warning in st_collection_extract.sf(jawa, "POLYGON"): x is already of type
## POLYGON.
# Menghapus vertex duplikat agar adjacency valid
jawa <- st_buffer(jawa, dist = 0)
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
coord <- st_coordinates(st_centroid(jawa$geometry))
## Warning in st_centroid.sfc(jawa$geometry): st_centroid does not give correct
## centroids for longitude/latitude data

4.1 Persinggungan

4.1.1 Queen Contiguity

#Queen Contiguity
queen <- poly2nb(jawa, queen = TRUE)
## although coordinates are longitude/latitude, st_intersects assumes that they
## are planar
W.queen <- nb2listw(queen, style = 'W', zero.policy = TRUE)

#uji index moran
mi_q <- moran.test(jawa$TK, W.queen, randomisation = TRUE, alternative = "greater", zero.policy = TRUE); mi_q
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: W.queen  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 6.8641, p-value = 3.346e-12
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.475128985      -0.009433962       0.004983496
plot(st_geometry(jawa), col = "white", border = "black")
plot(queen, st_coordinates(st_centroid(st_geometry(jawa))), add = TRUE, col = "red")
## Warning in st_centroid.sfc(st_geometry(jawa)): st_centroid does not give
## correct centroids for longitude/latitude data

4.1.2 Rook Contiguity

#Rook Contiguity
rook <- poly2nb(jawa, queen = FALSE)
## although coordinates are longitude/latitude, st_intersects assumes that they
## are planar
W.rook <- nb2listw(rook, style = 'W', zero.policy = TRUE)

#uji index moran
mi_r <- moran.test(jawa$TK, W.rook, randomisation = TRUE, alternative = "greater", zero.policy = TRUE); mi_r
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: W.rook  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 5.367, p-value = 4.002e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.432425770      -0.010101010       0.006798488
plot(st_geometry(jawa), col = "white", border = "black")
plot(rook, st_coordinates(st_centroid(st_geometry(jawa))), add = TRUE, col = "red")
## Warning in st_centroid.sfc(st_geometry(jawa)): st_centroid does not give
## correct centroids for longitude/latitude data

4.2 Jarak

4.2.1 k-NN

#k-Nearest Neighbors dengan k=4
knn <- knn2nb(knearneigh(coord, k = 4, longlat = TRUE))
knn.s <- nb2listw(knn, style = 'W')

#uji index moran
m.knn <- moran(jawa$TK, knn.s, n = length(knn.s$neighbours), S0 = Szero(knn.s))
mi_knn <- moran.test(jawa$TK, knn.s, randomisation = TRUE, alternative = "greater");mi_knn
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: knn.s    
## 
## Moran I statistic standard deviate = 6.2159, p-value = 2.551e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.360051793      -0.008474576       0.003515000
plot(st_geometry(jawa), col = "white", border = "black")
plot(knn, coord, add = TRUE, col = "red")
points(coord, pch = 10)

4.2.2 RDW

#Radial Distance Weigth
W.dmax <- dnearneigh(coord,0,60,longlat=TRUE)
W.dmax.s <- nb2listw(W.dmax,style='W', zero.policy = TRUE)

#uji index moran
m.dmax <- moran(jawa$TK, W.dmax.s, n = length(W.dmax.s$neighbours), S0 = Szero(W.dmax.s))
mi_dmax <- moran.test(jawa$TK, W.dmax.s, randomisation = TRUE, alternative = "greater")
mi_dmax 
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: W.dmax.s  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 7.6029, p-value = 1.448e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.376044794      -0.008547009       0.002558828
plot(st_geometry(jawa), col = "white", border = "black")
plot(W.dmax, coord, add = TRUE, col = "red")
points(coord, pch = 10)

4.2.3 IDW

#Invers Distance Weight alpha 1
djarak <- dist(coord)
m.djarak <- as.matrix(djarak)
alpha1 = 1
W.idw <- 1/(m.djarak^alpha1)
diag(W.idw) <- 0
rtot <- rowSums(W.idw, na.rm = TRUE)
W.idw.sd <- W.idw/rtot #row-normalized
W.idw.s <- mat2listw(W.idw.sd, style = 'W')

#uji index moran
m.idw <- moran(jawa$TK, W.idw.s, n = length(W.idw.s$neighbours), S0 = Szero(W.idw.s))
mi_idw <- moran.test(jawa$TK, W.idw.s, randomisation = TRUE, alternative = "greater")
mi_idw 
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 7.7709, p-value = 3.896e-15
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1497900058     -0.0084745763      0.0004147859

4.2.4 EDW

#Exponential Distance Weigth alpha 1
alpha <- 1
W.e <- exp((-alpha)*m.djarak)
diag(W.e) <- 0
rtot <- rowSums(W.e, na.rm = TRUE)
W.e.sd <- W.e/rtot #row-normalized
W.e.s <- mat2listw(W.e.sd, style = 'W')

#uji index moran
m.edw <- moran(jawa$TK, W.e.s, n = length(W.e.s$neighbours), S0 = Szero(W.e.s))
mi_edw <- moran.test(jawa$TK, W.e.s, randomisation = TRUE, alternative = "greater")
mi_edw 
## 
##  Moran I test under randomisation
## 
## data:  jawa$TK  
## weights: W.e.s    
## 
## Moran I statistic standard deviate = 10.553, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1398903194     -0.0084745763      0.0001976647

4.3 Rekap

# Rekapan
MatrikBobot <- c("Queen", "Rook", "KNN (k=4)", "Radial distance weight", "Invers distance weight", "Exponential distance weight")

IndeksMoran <- c(mi_q$estimate[1], mi_r$estimate[1], mi_knn$estimate[1], mi_dmax$estimate[1], mi_idw$estimate[1], mi_edw$estimate[1])

pv <- c(mi_q$p.value, mi_r$p.value, mi_knn$p.value, mi_dmax$p.value, mi_idw$p.value, mi_edw$p.value)

M <- cbind.data.frame(MatrikBobot,IndeksMoran,"p-value"=pv);M
##                   MatrikBobot IndeksMoran      p-value
## 1                       Queen   0.4751290 3.345763e-12
## 2                        Rook   0.4324258 4.002361e-08
## 3                   KNN (k=4)   0.3600518 2.551131e-10
## 4      Radial distance weight   0.3760448 1.447796e-14
## 5      Invers distance weight   0.1497900 3.896293e-15
## 6 Exponential distance weight   0.1398903 2.465956e-26

Matriks pembobot Queen Contiguity dipilih karena menghasilkan nilai Indeks Moran paling optimum dibandingkan matriks pembobot spasial lainnya.

5 Model Regresi Linear

5.1 MKT

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
mkt <- lm (TK ~ BPJSPBI + PPTB + PRTPBPNT + PPKM + IPM, data = data_jawa)
summary(mkt)
## 
## Call:
## lm(formula = TK ~ BPJSPBI + PPTB + PRTPBPNT + PPKM + IPM, data = data_jawa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8068 -1.5374  0.0201  1.2346  7.2854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.04523   10.44407   2.207 0.029369 *  
## BPJSPBI      0.05653    0.01374   4.115 7.37e-05 ***
## PPTB        -0.13364    0.03085  -4.331 3.22e-05 ***
## PRTPBPNT     0.05076    0.01445   3.513 0.000638 ***
## PPKM         0.14716    0.07061   2.084 0.039393 *  
## IPM         -0.26230    0.09301  -2.820 0.005672 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.278 on 113 degrees of freedom
## Multiple R-squared:  0.6581, Adjusted R-squared:  0.6429 
## F-statistic: 43.49 on 5 and 113 DF,  p-value: < 2.2e-16

Berdasarkan penduga parameter menggunakan MKT, diperoleh seluruh peubah berpengaruh terhadap persentase penduduk miskin. peubah BPJSPBI, PRTPBPNT, dan PPKM berpengaruh positif terhadap persentase penduduk miskin, sedangkan peubah PPTB dan IPM berpengaruh negatif terhadap persentase penduduk miskin. Berikut merupakan model MKT: \[ \hat{y}_i = 23.0452 + 0.0565 X_1 - 0.1336 X_2 + 0.0508 X_3 + 0.1472 X_4 - 0.2623 X_5 \]

rmse_mkt <- sqrt(mean((mkt$residuals)^2)); rmse_mkt
## [1] 2.219427

5.2 Uji Multikolinearitas

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
car::vif(mkt)
##  BPJSPBI     PPTB PRTPBPNT     PPKM      IPM 
## 1.076847 1.180005 1.253008 5.322379 5.555501

Uji multikolinearitas pada peubah penjelas menghasilkan nilai VIF < 10, artinya setiap pengaruh dari masing-masing peubah penjelas terhadap peubah respon dapat dilihat secara terpisah dan lebih akurat.

6 Deteksi ketergantungan spasial

6.1 Uji Moran’s I

lm.morantest(mkt, W.queen, alternative="greater")
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = TK ~ BPJSPBI + PPTB + PRTPBPNT + PPKM + IPM, data =
## data_jawa)
## weights: W.queen
## 
## Moran I statistic standard deviate = 3.9028, p-value = 4.754e-05
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.247784198     -0.023502296      0.004831661

6.2 Uji Lagrange Multiplier

lm_test <- lm.LMtests(mkt, W.queen, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
## Please update scripts to use lm.RStests in place of lm.LMtests
summary(lm_test)
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = TK ~ BPJSPBI + PPTB + PRTPBPNT + PPKM + IPM, data =
## data_jawa)
## test weights: listw
##  
##          statistic parameter   p.value    
## RSerr    14.637486         1 0.0001303 ***
## RSlag     4.044323         1 0.0443202 *  
## adjRSerr 10.648362         1 0.0011017 ** 
## adjRSlag  0.055199         1 0.8142513    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7 Model SAR

library(spatialreg)
sar <- lagsarlm(mkt, data = data_jawa, W.queen, zero.policy=TRUE)
summary(sar)
## 
## Call:
## lagsarlm(formula = mkt, data = data_jawa, listw = W.queen, zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.0262312 -1.3491877  0.0025808  1.2632787  6.9210092 
## 
## Type: lag 
## Regions with no neighbours included:
##  1 25 26 63 65 68 102 103 104 105 107 108 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 16.255695  10.580562  1.5364 0.1244468
## BPJSPBI      0.057137   0.013178  4.3357 1.453e-05
## PPTB        -0.107612   0.031149 -3.4547 0.0005509
## PRTPBPNT     0.051390   0.013856  3.7088 0.0002083
## PPKM         0.161635   0.067902  2.3804 0.0172925
## IPM         -0.210365   0.093528 -2.2492 0.0244994
## 
## Rho: 0.11259, LR test value: 3.7637, p-value: 0.052378
## Asymptotic standard error: 0.059834
##     z-value: 1.8818, p-value: 0.059866
## Wald statistic: 3.5411, p-value: 0.059866
## 
## Log likelihood: -261.8445 for lag model
## ML residual variance (sigma squared): 4.7593, (sigma: 2.1816)
## Number of observations: 119 
## Number of parameters estimated: 8 
## AIC: 539.69, (AIC for lm: 541.45)
## LM test for residual autocorrelation
## test value: 8.73, p-value: 0.0031301

Berdasarkan penduga parameter menggunakan SAR, diperoleh seluruh peubah berpengaruh terhadap persentase penduduk miskin. peubah BPJSPBI, PRTPBPNT, dan PPKM berpengaruh positif terhadap persentase penduduk miskin, sedangkan peubah PPTB dan IPM berpengaruh negatif terhadap persentase penduduk miskin. Berikut merupakan model SAR:

\[ \hat{y}_i = 0.1126 \sum_{j=1,\, j\neq i}^{n} w_{ij} y_j + 16.2557 + 0.0571 X_1 - 0.1076 X_2 + 0.0514 X_3 + 0.1616 X_4 - 0.2104 \]

8 Model SEM

sem <- errorsarlm(mkt, data = data_jawa, W.queen, zero.policy=TRUE)
summary(sem)
## 
## Call:errorsarlm(formula = mkt, data = data_jawa, listw = W.queen, 
##     zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.1905378 -1.2884992 -0.0013419  0.9056905  5.7876878 
## 
## Type: error 
## Regions with no neighbours included:
##  1 25 26 63 65 68 102 103 104 105 107 108 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 25.723107   9.350795  2.7509 0.0059432
## BPJSPBI      0.049559   0.013441  3.6871 0.0002268
## PPTB        -0.066182   0.030506 -2.1695 0.0300449
## PRTPBPNT     0.039859   0.014114  2.8241 0.0047421
## PPKM         0.145388   0.061062  2.3810 0.0172670
## IPM         -0.327985   0.087035 -3.7684 0.0001643
## 
## Lambda: 0.48062, LR test value: 16.301, p-value: 5.4039e-05
## Asymptotic standard error: 0.097597
##     z-value: 4.9246, p-value: 8.4546e-07
## Wald statistic: 24.251, p-value: 8.4546e-07
## 
## Log likelihood: -255.5759 for error model
## ML residual variance (sigma squared): 4.0559, (sigma: 2.0139)
## Number of observations: 119 
## Number of parameters estimated: 8 
## AIC: 527.15, (AIC for lm: 541.45)

Berdasarkan penduga parameter menggunakan SEM, diperoleh seluruh peubah berpengaruh terhadap persentase penduduk miskin. peubah BPJSPBI, PRTPBPNT, dan PPKM berpengaruh positif terhadap persentase penduduk miskin, sedangkan peubah PPTB dan IPM berpengaruh negatif terhadap persentase penduduk miskin. Berikut merupakan model SEM:

\[ \hat{y}_i = 25.7231 + 0.0496 X_1 - 0.0662 X_2 + 0.0399 X_3 + 0.1454 X_4 - 0.3280 X_5 + \hat{u}_i \]

\[ \hat{u}_i = 0.4806 \sum_{j=1,\, j\neq i}^{n} w_{ij} y_j + \varepsilon_i \]

9 Perbandingan Model Terbaik

#SAR
sar_sst <- sum((jawa$TK - mean(jawa$TK))^2)

rmse_sar <- sqrt(mean((sar$residuals)^2)); rmse_sar
## [1] 2.181577
#SEM
sem_sst <- sum((jawa$TK - mean(jawa$TK))^2)

rmse_sem <- sqrt(mean((sem$residuals)^2)); rmse_sem
## [1] 2.013931
model <- c("mkt", "sar", "sem")

AIC <- c(541.45, 539.69, 527.15)
RMSE <- c(rmse_mkt, rmse_sar, rmse_sem)

rekap <- cbind.data.frame(model, AIC, RMSE); rekap
##   model    AIC     RMSE
## 1   mkt 541.45 2.219427
## 2   sar 539.69 2.181577
## 3   sem 527.15 2.013931

perbandingan model menunjukkan bahwa sem merupakan model terbaik dalam memodelkan persentase penduduk miskin di kabupaten/kota Pulau Jawa.