Email             :
RPubs            : https://rpubs.com/juliansalomo/
Department  : Business Statistics
Address         : ARA Center, Matana University Tower
                         Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.


Berikut merupakan library-library yang akan digunakan dalam pemodelan spasial

library(rgdal)            # untuk membuka file format .sdh (peta)
library(raster)
library(openxlsx)         # untuk membuka file excel
library(spdep)
library(sp)
library(raster)
library(rgeos)
library(rspatial)
library(latticeExtra)
library(RColorBrewer)
library(spatialreg)
library(DT)

1 Input Data

Disini saya akan memberikan contoh menggunakan data kemiskinan di Jawa Barat yang dapat di download di link ini Pertama kita membaca file data lengkap beserta peta Jawa Barat dengan fungsi berikut

data.jabar = read.xlsx("Jabar Data (gabung).xlsx")
jabar2 <- readOGR(dsn = "petaJabar2" , layer = "Jabar2")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Julian Salomo\Matana\Teknik Sampling\jabar1\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields
datatable(data.jabar, htmltools::em('Table data'),     
          extension = 'FixedColumns',
          options = list(scrollX = TRUE, fixedColums = TRUE))

2 Plot peta data

Sekarang kita akan memvisualisasikan kemiskinan di tahun 2016 di provinsi Jawa Barat. Kemiskinan dari dataset data.jabar berada di kolom p.miskin16. Kita akan menambahkan kolom ini ke dalam dataset jabar2 dengan script berikut

jabar2$miskin2<- data.jabar$p.miskin16

Selanjutnya kita akan memvisualisasikannya dengan package sp dengan fungsi berikut

kota <- aggregate(jabar2, "KABKOT")
grps <- 10           #mengelompokkan menjadi 10 interval
brks <- quantile(jabar2$miskin2, 0:(grps-1)/(grps-1), na.rm=TRUE)
p <- spplot(jabar2, "miskin2", at=brks, col.regions=rev(brewer.pal(grps, "RdBu")), col="transparent" )
p + layer(sp.polygons(kota))

Dengan plot di atas kita bisa melihat wilayah-wilayah yang terkonsentrasi dengan kemiskinan. Dan kita bisa lihat bahwa di bagian yang merah pekat, area sekitarnya pun memiliki tingkat kemiskinan yang cukup tinggi.

3 Pemodelan Regresi Spasial

Di bagian sini kita akan membangun model linier untuk data kita. Model spasial secara global terdiri dari:

  • SEM: \(y=XB+u\)
  • SAR: \(y=\rho Wy+XB+\varepsilon\)
  • SARMA: \(y=\rho Wy+XB+u\)

3.1 Model SEM

w<-poly2nb(jabar2)
sem<-errorsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(sem)
## 
## Call:errorsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, 
##     listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.24699 -1.36122 -0.13809  1.15579  7.03019 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 36.88515    7.80246  4.7274 2.274e-06
## EYS2016     -2.17498    0.60042 -3.6224 0.0002919
## 
## Lambda: 0.61793, LR test value: 8.7676, p-value: 0.0030663
## Asymptotic standard error: 0.1576
##     z-value: 3.9208, p-value: 8.8267e-05
## Wald statistic: 15.372, p-value: 8.8267e-05
## 
## Log likelihood: -58.12466 for error model
## ML residual variance (sigma squared): 4.5459, (sigma: 2.1321)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 124.25, (AIC for lm: 131.02)

3.2 Model SAR

sar<-lagsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(sar)
## 
## Call:
## lagsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.17670 -0.93185 -0.11318  0.91353  7.69131 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 28.22898    7.66250  3.6840 0.0002296
## EYS2016     -1.94997    0.57325 -3.4016 0.0006700
## 
## Rho: 0.59078, LR test value: 7.9343, p-value: 0.0048507
## Asymptotic standard error: 0.1559
##     z-value: 3.7894, p-value: 0.00015101
## Wald statistic: 14.36, p-value: 0.00015101
## 
## Log likelihood: -58.54132 for lag model
## ML residual variance (sigma squared): 4.7513, (sigma: 2.1798)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 125.08, (AIC for lm: 131.02)
## LM test for residual autocorrelation
## test value: 0.036687, p-value: 0.8481

3.3 Model SARMA

gsm<-sacsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(gsm)
## 
## Call:
## sacsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, listw = nb2listw(w))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0817 -1.2815 -0.1627  1.1317  6.7128 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 38.02312    8.19667  4.6388 3.504e-06
## EYS2016     -2.15311    0.59378 -3.6261 0.0002877
## 
## Rho: -0.14689
## Asymptotic standard error: 0.38983
##     z-value: -0.3768, p-value: 0.70633
## Lambda: 0.69352
## Asymptotic standard error: 0.23014
##     z-value: 3.0135, p-value: 0.0025829
## 
## LR test value: 8.8931, p-value: 0.011719
## 
## Log likelihood: -58.06189 for sac model
## ML residual variance (sigma squared): 4.3254, (sigma: 2.0797)
## Number of observations: 26 
## Number of parameters estimated: 5 
## AIC: 126.12, (AIC for lm: 131.02)