Spatial Autocorrelation

Joint Count Statistics

Joint count merupakan metode paling dasar dalam menentukan autokorelasi spasial antar area. Penjelasan lebih lengkap mengenai metode ini dapat dilihat pada Lizazaro (2016). Ilustrasi berikut ini diadaptasi dari Lizazaro (2016). Untuk melakukan analisis joint count, terlebih dulu kita load beberapa package pada R.

1. Loading Packages

library(raster) #membuat dan memanipulasi data raster (grid)
## Warning: package 'raster' was built under R version 4.5.2
## Loading required package: sp
library(sp) #kelas dasar untuk data spasial (Spatial* objects) di R
library(spdep) #paket utama untuk analisis autokorelasi spasial (Moran’s I, Geary’s C, LISA, dll).
## Warning: package 'spdep' was built under R version 4.5.2
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.5.2
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE

2. Membuat Data Contoh (Matriks 6×6)

Sebagai ilustrasi, kita akan membuat data contoh dengan menggunakan syntax berikut.

pri <- rep(1,12)
seg <- rep(0,4)
ter <- rep(1,2)
cua <- rep(0,4)
qui <- rep(1,2)
sex <-  rep(0,12)

A <- matrix(c(pri, seg, ter, cua ,qui, sex), nrow=6, byrow=FALSE)

A
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    0    0    0    0
## [2,]    1    1    0    0    0    0
## [3,]    1    1    0    0    0    0
## [4,]    1    1    0    0    0    0
## [5,]    1    1    1    1    0    0
## [6,]    1    1    1    1    0    0

Hasilnya adalah grid 6×6 dengan nilai 0 dan 1 yang membentuk pola huruf “A” dari nilai 1 (hitam) dan latar belakang 0 (putih). Selanjutnya matriks A dikonversi menjadi raster menggunakan fungsi raster().

3. Mengubah Matriks menjadi Raster

rA <- raster(A)

rA
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

Berikut ini adalah plot dari area A.

plot(rA)
text(coordinates(rA), labels=rA[], cex=1.5)

4. Mengubah Raster menjadi Polygon (untuk analisis tetangga)

Kriteria ketetanggaan selanjutnya digunakan untuk mengkuantifikasi pola kedekatan data tersebut.

pA <- rasterToPolygons(rA, dissolve=FALSE)

pA
## class       : SpatialPolygonsDataFrame 
## features    : 36 
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 1
## names       : layer 
## min values  :     0 
## max values  :     1

Setiap sel raster diubah menjadi polygon persegi (grid kotak-kotak). Ini diperlukan karena fungsi poly2nb() hanya bekerja pada polygon, bukan raster langsung.

5. Membuat Matriks Tetangga (Spatial Weights)

a. Queen contiguity (8 tetangga: horizontal, vertikal, diagonal)

Seandainya kita akan menggunakan kriteria queen contiguity, maka dapat dilakukan dengan syntax berikut.

spA <- SpatialPolygons(pA@polygons)

nb1 <- poly2nb(spA, queen = T)

nb1
## Neighbour list object:
## Number of regions: 36 
## Number of nonzero links: 220 
## Percentage nonzero weights: 16.97531 
## Average number of links: 6.111111

Selanjutnya kita dapat memvisualisasikan link dari ketetanggannya.

par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb1, xy, col='red', lwd=2, add=TRUE)

b. Rook contiguity (hanya 4 tetangga: horizontal & vertikal)

Seandainya kriteria yang digunakan adalah rook contiguity.

nb2 <- poly2nb(spA, queen = F)

nb2
## Neighbour list object:
## Number of regions: 36 
## Number of nonzero links: 120 
## Percentage nonzero weights: 9.259259 
## Average number of links: 3.333333
par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb2, xy, col='green', lwd=2, add=TRUE)

Perbedaan: queen lebih banyak menghubungkan sel (termasuk diagonal), rook lebih ketat.

6. Membuat Spatial Weights List (matriks bobot)

Agar dapat mengidentifikasi pola autokorelasi pada data ini, maka matriks bobot perlu disusun dengan tipe biner.

wl1 <- nb2listw(nb1, style='B')# Queen, biner (1 jika tetangga, 0 jika tidak)

wl2 <- nb2listw(nb2, style='B')# Rook, biner

Style = “B” → biner (paling umum untuk joint count dan Moran’s I dasar).

7. Joint Count Statistics (khusus data biner 0/1)

jc_test1 <- joincount.test(as.factor(pA$layer), wl1)# Queen

jc_test1
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## 
## Std. deviate for 0 = 5.1529, p-value = 1.282e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              53.00000              33.17460              14.80263 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## 
## Std. deviate for 1 = 4.7634, p-value = 9.52e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              37.00000              20.95238              11.34999

Kesimpulan : Karena p-value sangat kecil (< 0.001) → tolak H0 → ada autokorelasi spasial positif (clustering nilai sama).

jc_test2 <- joincount.test(as.factor(pA$layer), wl2) # Rook

jc_test2
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl2 
## 
## Std. deviate for 0 = 5.4677, p-value = 2.28e-08
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             30.000000             18.095238              4.740611 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl2 
## 
## Std. deviate for 1 = 5.1203, p-value = 1.525e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             22.000000             11.428571              4.262554

Kesimpulan : Karena p-value sangat kecil (< 0.001) → tolak H0 → ada autokorelasi spasial positif (clustering nilai sama).

Berdasarkan kedua output di atas, dengan p-value yang sangat kecil, artinya kita dapat menolak hipotesis nol yang menyatakan bahwa tidak terdapat autokorelasi. Sesuai dengan output di atas, alternative hypothesis: greater, artinya kita dapat menyimpulkan bahwa terdapat cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf nyata 5%.

Pengujian hipotesis dapat pula dilakukan dengan melibatkan algoritma monte carlo seperti di bawah ini.

set.seed(123)
jc_test3 <- joincount.mc(as.factor(pA$layer), wl1, nsim=99) #Mengacak nilai 99 kali untuk simulasi distribusi di bawah H0.

jc_test3
## 
##  Monte-Carlo simulation of join-count statistic
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## number of simulations + 1: 100 
## 
## Join-count statistic for 0 = 53, rank of observed statistic = 100,
## p-value = 0.01
## alternative hypothesis: greater
## sample estimates:
##     mean of simulation variance of simulation 
##               33.01010               15.15296 
## 
## 
##  Monte-Carlo simulation of join-count statistic
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## number of simulations + 1: 100 
## 
## Join-count statistic for 1 = 37, rank of observed statistic = 100,
## p-value = 0.01
## alternative hypothesis: greater
## sample estimates:
##     mean of simulation variance of simulation 
##               20.37374               12.37930

Global Autocorrelation

Moran’s I (bisa untuk data kontinu maupun biner)

Sebagai ilustrasi, kita masih akan menggunakan data yang sama. Seandainya kita ingin menguji autokorelasi menggunakan pendekatan indeks moran, maka kita dapat menggunakan fungsi moran.test().

I1 <- moran.test(pA$layer,wl1)

I1
## 
##  Moran I test under randomisation
## 
## data:  pA$layer  
## weights: wl1    
## 
## Moran I statistic standard deviate = 7.4654, p-value = 4.151e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.624090909      -0.028571429       0.007643049

Hasil uji Moran : Nilai I positif (0.62409090) dan signifikan (pvalue < 5%) → autokorelasi spasial positif global (nilai tinggi cenderung berdekatan).

Sama halnya seperti yang telah diperlihatkan pada uji dengan metode joint count, uji moran juga dapat dilakukan dengan melibatkan simulasi monte carlo.

set.seed(123)
MC<- moran.mc(pA$layer, wl1, nsim=599)

# View results (including p-value)
MC
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  pA$layer 
## weights: wl1  
## number of simulations + 1: 600 
## 
## statistic = 0.62409, observed rank = 600, p-value = 0.001667
## alternative hypothesis: greater

Global Geary’s C

Geary’s C merupakan alternatif dari indeks Moran, yang memiliki nilai antara 0 s.d 2. Nilai 0 menunjukkan autokorelasi positif, 1 menunjukkan tidak ada autokorelasi, dan 2 menunjukkan autokorelasi negatif.

C1 <- geary.test(pA$layer,wl1)

C1
## 
##  Geary C test under randomisation
## 
## data:  pA$layer 
## weights: wl1   
## 
## Geary C statistic standard deviate = 7.4869, p-value = 3.526e-14
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.357954545       1.000000000       0.007354046

Kesimpulan : Nilai C < 1 dan signifikan → autokorelasi positif (sama seperti Moran’s I positif).

GS1 <- geary.mc(pA$layer, wl1, nsim=599)

GS1
## 
##  Monte-Carlo simulation of Geary C
## 
## data:  pA$layer 
## weights: wl1  
## number of simulations + 1: 600 
## 
## statistic = 0.35795, observed rank = 1, p-value = 0.001667
## alternative hypothesis: greater

Local Autocorrelation (LISA)

Local Moran’s I

Pendekatan ini termasuk ke dalam Local Indicators for Spatial Association (LISA), yang mengindentifikasi autokorelasi pada tingkat lokal.

oid <- order(pA$layer)
resI <- localmoran(pA$layer, wl1)
head(resI)
##     Ii        E.Ii   Var.Ii      Z.Ii Pr(z != E(Ii))
## 1 3.75 -0.10714286 3.500600 2.0615528     0.03925033
## 2 1.75 -0.17857143 5.469688 0.8246211     0.40958672
## 3 0.40 -0.11428571 3.547275 0.2730593     0.78480760
## 4 4.00 -0.11428571 3.547275 2.1844747     0.02892738
## 5 4.00 -0.11428571 3.547275 2.1844747     0.02892738
## 6 2.40 -0.06857143 2.270256 1.6383560     0.10134744
pA$z.li <- resI[,4]
pA$pvalue <- resI[,5]
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="z.li", col.regions=lm.palette(20), main="Local Moran")

Intepretasi :

Pola yang Terlihat pada Heatmap Ini

  1. Bagian huruf “A” yang terbentuk dari nilai 1 → Hampir seluruh sel yang membentuk huruf A berwarna merah tua hingga merah cerah (z > +1.5, bahkan banyak yang > +2.5). Artinya: Sel-sel bernilai tinggi (1) dikelilingi oleh sel-sel bernilai tinggi lainnya → kluster hotspot yang sangat signifikan secara statistik (p < 0.05, bahkan p < 0.01).

  2. Latar belakang (sel bernilai 0) → Sebagian besar berwarna kuning muda hingga putih, artinya z-value mendekati 0 → tidak ada autokorelasi lokal yang signifikan. Ini wajar karena nilai 0 tersebar cukup acak di luar huruf A, sehingga tidak membentuk kluster Low-Low yang kuat.

  3. Beberapa sel di dalam atau di tepi huruf A berwarna oranye/kuning muda Ini biasanya terjadi pada sel-sel di ujung-ujung atau sudut huruf A yang memiliki sedikit tetangga bernilai tinggi (misalnya hanya 2–3 tetangga 1), sehingga nilai-z-nya tidak sekuat sel-sel di tengah-tengah huruf A.

moran.plot(pA$layer,wl1)

Intepretasi :

  • Sumbu X: nilai sel itu sendiri
  • Sumbu Y: rata-rata nilai tetangganya (spatially lagged)
  • Kuadran kanan atas = High-High → cluster positif

Getis-Ord Gi

Menurut Mendez (2020), pendekatan Getis-ord Gi dapat membantu mengidentifikasi pola penggerombolan berdasarkan ukuran autokorelasi pada level lokal.

local_g <- localG(pA$layer, wl1)

local_g
##  [1]  2.0615528  0.8246211 -0.2730593 -2.1844747 -2.1844747 -1.6383560
##  [7]  2.7487371  1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [13]  2.7487371  1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [19]  2.7487371  2.0615528  1.0694824 -1.3197868 -2.1162099 -2.1844747
## [25]  2.7487371  2.8632678  2.0615528 -0.3435921 -1.3197868 -2.1844747
## [31]  2.0615528  2.7487371  2.7487371  0.8246211 -0.2730593 -1.6383560
## attr(,"internals")
##              Gi      E(Gi)       V(Gi)      Z(Gi) Pr(z != E(Gi))
##  [1,] 0.2000000 0.08571429 0.003073229  2.0615528    0.039250330
##  [2,] 0.2000000 0.14285714 0.004801921  0.8246211    0.409586724
##  [3,] 0.1250000 0.14285714 0.004276711 -0.2730593    0.784807601
##  [4,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
##  [5,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
##  [6,] 0.0000000 0.08571429 0.002737095 -1.6383560    0.101347443
##  [7,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
##  [8,] 0.3333333 0.22857143 0.006914766  1.2598378    0.207727870
##  [9,] 0.1875000 0.22857143 0.006158463 -0.5233637    0.600721155
## [10,] 0.0000000 0.22857143 0.006158463 -2.9126330    0.003583956
## [11,] 0.0000000 0.22857143 0.006158463 -2.9126330    0.003583956
## [12,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
## [13,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
## [14,] 0.3333333 0.22857143 0.006914766  1.2598378    0.207727870
## [15,] 0.1875000 0.22857143 0.006158463 -0.5233637    0.600721155
## [16,] 0.0000000 0.22857143 0.006158463 -2.9126330    0.003583956
## [17,] 0.0000000 0.22857143 0.006158463 -2.9126330    0.003583956
## [18,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
## [19,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
## [20,] 0.4000000 0.22857143 0.006914766  2.0615528    0.039250330
## [21,] 0.3125000 0.22857143 0.006158463  1.0694824    0.284852347
## [22,] 0.1250000 0.22857143 0.006158463 -1.3197868    0.186906206
## [23,] 0.0625000 0.22857143 0.006158463 -2.1162099    0.034326960
## [24,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
## [25,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
## [26,] 0.4666667 0.22857143 0.006914766  2.8632678    0.004192960
## [27,] 0.4000000 0.22857143 0.006914766  2.0615528    0.039250330
## [28,] 0.2000000 0.22857143 0.006914766 -0.3435921    0.731153040
## [29,] 0.1250000 0.22857143 0.006158463 -1.3197868    0.186906206
## [30,] 0.0000000 0.14285714 0.004276711 -2.1844747    0.028927383
## [31,] 0.2000000 0.08571429 0.003073229  2.0615528    0.039250330
## [32,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
## [33,] 0.3333333 0.14285714 0.004801921  2.7487371    0.005982535
## [34,] 0.2000000 0.14285714 0.004801921  0.8246211    0.409586724
## [35,] 0.1250000 0.14285714 0.004276711 -0.2730593    0.784807601
## [36,] 0.0000000 0.08571429 0.002737095 -1.6383560    0.101347443
## attr(,"cluster")
##  [1] High High Low  Low  Low  Low  High High Low  Low  Low  Low  High High Low 
## [16] Low  Low  Low  High High Low  Low  Low  Low  High High High High Low  Low 
## [31] High High High High Low  Low 
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = pA$layer, listw = wl1)
## attr(,"class")
## [1] "localG"

Output di atas menghasilkan z-score, yang biasanya disajikan secara visual untuk mengidentifikasi cluster maupun hotspot.

pA$localg <- as.numeric(local_g)
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="localg", col.regions=lm.palette(20), main="Local Gi")

Exercise

Sebagai latihan, Anda dipersilahkan menggunakan data yang tersedia pada: https://github.com/raoy/SpatialReg . Terdapat dua data yang harus Anda download, yaitu:

  1. Jabar Data (gabung).xlsx

  2. petaJabar2.zip

Data pertama (dengan format Excel) menyimpan data kependudukan yang diperoleh dari BPS. Sedangkan data kedua merupakan data shapefile berisi peta Provinsi Jawa Barat. Silahkan manfaatkan kedua data tersebut untuk mengeksplorasi pola depedensi spasial untuk peubah kemiskinan antar kota/kabupaten di Jawa Barat pada tahun 2015. Data tersebut terdapat pada kolom I dengan nama kolom p.miskin15 pada file Excel.

1. Uji Autokorelasi Spasial

# ==============================
# ANALISIS AUTOKORELASI SPASIAL 
# Provinsi Jawa Barat (Data Kemiskinan 2015–2016)
# ==============================

# 1. Load library yang benar-benar dibutuhkan
library(sf)         # untuk data spasial modern (geopackage/shapefile)
library(spdep)      # untuk autokorelasi spasial (Moran, LISA, dll)
library(readxl)     # membaca file Excel
library(ggplot2)    # opsional, untuk plot yang lebih bagus
library(classInt)   # untuk breaks warna yang bagus (opsional)
# 2. Baca shapefile Jawa Barat
jabar <- st_read("D:/Kuliah/IPB 2025 Semester 3/Analisis Spasial/Praktikum/Sesi Uas/Praktikum 9/Regresi Spasial/petaJabar2/petaJabar2/Jabar2.shp")
## Reading layer `Jabar2' from data source 
##   `D:\Kuliah\IPB 2025 Semester 3\Analisis Spasial\Praktikum\Sesi Uas\Praktikum 9\Regresi Spasial\petaJabar2\petaJabar2\Jabar2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS:           NA
# Cek dan pastikan CRS (jika belum ada, set ke WGS84)
if (is.na(st_crs(jabar))) {
  jabar <- st_set_crs(jabar, 4326)
}
st_crs(jabar)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
names(jabar)
## [1] "PROVNO"   "KABKOTNO" "KODE2010" "PROVINSI" "KABKOT"   "SUMBER"   "IDSP2010"
## [8] "geometry"
# Plot peta kosong
plot(st_geometry(jabar), main = "Batas Administrasi Jawa Barat")

# 3. Baca data Excel
data.jabar <- read_excel("D:/Kuliah/IPB 2025 Semester 3/Analisis Spasial/Praktikum/Sesi Uas/Praktikum 9/Regresi Spasial/Jabar Data.xlsx", 
                         sheet = "data")
names(data.jabar)
##  [1] "PROVNO"     "KABKOTNO"   "KODE2010"   "PROVINSI"   "KABKOT"    
##  [6] "IDSP2010"   "Long"       "Lat"        "p.miskin15" "p.miskin16"
## [11] "j.miskin15" "j.miskin16" "AHH2015"    "AHH2016"    "EYS2015"   
## [16] "EYS2016"    "MYS2015"    "MYS2016"    "EXP2015"    "EXP2016"   
## [21] "APM.SD15"   "APM.SMP15"  "APM.SMA15"  "APM.PT15"   "APK.SD15"  
## [26] "APK.SMP15"  "APK.SMA15"  "APK.PT15"   "APS.USIA15" "APS.USIA2" 
## [31] "APS.USIA3"  "APS.USIA4"
# Pastikan nama kolom kunci sama persis (huruf besar/kecil & spasi)
# Misalnya di shapefile kolomnya "KABKOT", di Excel juga harus "KABKOT"
# Jika beda, samakan dulu:
#colnames(data.jabar)[colnames(data.jabar) == "Kabupaten"] <- "KABKOT"
# 4. Gabungkan data atribut Excel ke shapefile
jabar_join <- merge(jabar, data.jabar, by = "KABKOT", all.x = TRUE)

# Cek apakah penggabungan berhasil
head(jabar_join)
## Simple feature collection with 6 features and 38 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 107.1819 ymin: -7.437613 xmax: 108.6648 ymax: -6.688515
## Geodetic CRS:  WGS 84
##          KABKOT PROVNO.x KABKOTNO.x KODE2010.x PROVINSI.x
## 1       BANDUNG       32         04       3204 JAWA BARAT
## 2       BANDUNG       32         04       3204 JAWA BARAT
## 3       BANDUNG       32         73       3273 JAWA BARAT
## 4       BANDUNG       32         73       3273 JAWA BARAT
## 5 BANDUNG BARAT       32         17       3217 JAWA BARAT
## 6        BANJAR       32         79       3279 JAWA BARAT
##                         SUMBER IDSP2010.x PROVNO.y KABKOTNO.y KODE2010.y
## 1 SP2010_BADAN PUSAT STATISTIK       3204       32          4       3204
## 2 SP2010_BADAN PUSAT STATISTIK       3204       32         73       3273
## 3 SP2010_BADAN PUSAT STATISTIK       3273       32          4       3204
## 4 SP2010_BADAN PUSAT STATISTIK       3273       32         73       3273
## 5 SP2010_BADAN PUSAT STATISTIK       3217       32         17       3217
## 6 SP2010_BADAN PUSAT STATISTIK       3279       32         79       3279
##   PROVINSI.y IDSP2010.y     Long       Lat p.miskin15 p.miskin16 j.miskin15
## 1 JAWA BARAT       3204 107.6108 -7.099969   7.998184   7.613293     281.04
## 2 JAWA BARAT       3273 107.6366 -6.919135   4.611404   4.323270     114.12
## 3 JAWA BARAT       3204 107.6108 -7.099969   7.998184   7.613293     281.04
## 4 JAWA BARAT       3273 107.6366 -6.919135   4.611404   4.323270     114.12
## 5 JAWA BARAT       3217 107.4150 -6.897056  12.669240  11.710438     205.69
## 6 JAWA BARAT       3279 108.5665 -7.376302   7.409971   7.014890      13.42
##   j.miskin16 AHH2015 AHH2016 EYS2015 EYS2016 MYS2015 MYS2016 EXP2015 EXP2016
## 1     272.65   73.07   73.10   12.13   12.42    8.41    8.50    9375    9580
## 2     107.58   73.82   73.84   13.63   13.89   10.52   10.58   15609   15805
## 3     272.65   73.07   73.10   12.13   12.42    8.41    8.50    9375    9580
## 4     107.58   73.82   73.84   13.63   13.89   10.52   10.58   15609   15805
## 5     192.48   71.76   71.82   11.39   11.56    7.53    7.63    7522    7698
## 6      12.74   70.26   70.33   12.95   13.18    8.06    8.19    9476    9815
##   APM.SD15 APM.SMP15 APM.SMA15  APM.PT15 APK.SD15 APK.SMP15 APK.SMA15  APK.PT15
## 1 98.00367  82.13540  55.49216 17.042698 110.0396  89.63486  66.20417 19.736955
## 2 96.24954  87.22698  74.06181 39.604943 103.7611  98.32964  92.70930 45.223528
## 3 98.00367  82.13540  55.49216 17.042698 110.0396  89.63486  66.20417 19.736955
## 4 96.24954  87.22698  74.06181 39.604943 103.7611  98.32964  92.70930 45.223528
## 5 96.56357  79.39350  47.25078  6.934742 108.5784  93.17800  58.30069  7.943457
## 6 98.06811  86.53822  71.15240 15.854381 104.6057  97.17299  83.98661 19.289534
##   APS.USIA15 APS.USIA2 APS.USIA3 APS.USIA4                       geometry
## 1   99.91440  95.00882  60.55318  19.99148 MULTIPOLYGON (((107.75 -6.8...
## 2   99.26582  97.08581  85.68069  43.23437 MULTIPOLYGON (((107.75 -6.8...
## 3   99.91440  95.00882  60.55318  19.99148 MULTIPOLYGON (((107.6286 -6...
## 4   99.26582  97.08581  85.68069  43.23437 MULTIPOLYGON (((107.6286 -6...
## 5   99.62551  90.27611  57.73897  11.58464 MULTIPOLYGON (((107.5955 -6...
## 6  100.00000  95.77732  80.61464  21.27231 MULTIPOLYGON (((108.5594 -7...
summary(jabar_join$p.miskin15)   # pastikan tidak ada NA yang banyak
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.397   7.597   8.971   9.755  12.253  16.280
summary(jabar_join$p.miskin16)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.339   7.288   8.835   9.194  11.525  15.596
summary(jabar_join$EYS2016)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.56   12.04   12.44   12.66   13.38   13.89
# 5. Konversi ke Spatial (diperlukan untuk fungsi spdep klasik)
jabar.sp <- as(jabar_join, "Spatial")
jabar.sp
## class       : SpatialPolygonsDataFrame 
## features    : 38 
## extent      : 106.3705, 108.8338, -7.823398, -5.91377  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 38
## names       :      KABKOT, PROVNO.x, KABKOTNO.x, KODE2010.x, PROVINSI.x,                       SUMBER, IDSP2010.x, PROVNO.y, KABKOTNO.y, KODE2010.y, PROVINSI.y, IDSP2010.y,     Long,       Lat,       p.miskin15, ... 
## min values  :     BANDUNG,       32,         01,       3201, JAWA BARAT, SP2010_BADAN PUSAT STATISTIK,       3201,       32,          1,       3201, JAWA BARAT,       3201, 106.7101, -7.496892, 2.39749571489967, ... 
## max values  : TASIKMALAYA,       32,         79,       3279, JAWA BARAT, SP2010_BADAN PUSAT STATISTIK,       3279,       32,         79,       3279, JAWA BARAT,       3279, 108.5665, -6.215149, 16.2801622356294, ...
# 6. Buat matriks tetangga (Queen contiguity) dan spatial weights
nb_q <- poly2nb(jabar.sp, queen = TRUE)        # tetangga
summary(nb_q)                                  # lihat ada berapa link
## Neighbour list object:
## Number of regions: 38 
## Number of nonzero links: 226 
## Percentage nonzero weights: 15.65097 
## Average number of links: 5.947368 
## Link number distribution:
## 
##  1  3  4  5  6  7  8  9 13 
##  1  6  4  4 11  4  3  3  2 
## 1 least connected region:
## 6 with 1 link
## 2 most connected regions:
## 11 12 with 13 links
# Row-standardized weights (style = "W") → paling umum digunakan
listw <- nb2listw(nb_q, style = "W", zero.policy = TRUE)
listw
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38 
## Number of nonzero links: 226 
## Percentage nonzero weights: 15.65097 
## Average number of links: 5.947368 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 38 1444 38 14.48515 161.3108
# 7. Uji Autokorelasi Global: Moran’s I
# Kita gunakan variabel persentase penduduk miskin tahun 2015
moran_result <- moran.test(jabar_join$p.miskin16, listw, 
                           randomisation = TRUE, alternative = "greater")
print(moran_result)
## 
##  Moran I test under randomisation
## 
## data:  jabar_join$p.miskin16  
## weights: listw    
## 
## Moran I statistic standard deviate = 5.4849, p-value = 2.068e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.480818973      -0.027027027       0.008572773
# Monte Carlo (lebih akurat untuk p-value)
set.seed(123)
moran_mc <- moran.mc(jabar_join$p.miskin16, listw, nsim = 999, 
                     alternative = "greater")
print(moran_mc)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar_join$p.miskin16 
## weights: listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.48082, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
plot(moran_mc, main = "Monte Carlo Moran’s I")

# 8. Moran Scatterplot (untuk melihat outlier dan kluster)
moran.plot(jabar_join$p.miskin16, listw, 
           labels = jabar_join$KABKOT,
           xlab = "Persentase Penduduk Miskin 2016",
           ylab = "Spatially Lagged Persentase Miskin",
           main = "Moran Scatterplot - Kemiskinan Jawa Barat 2016")

Kesimpulan :

  1. Terdapat autokorelasi spasial positif yang sangat kuat. Artinya: kabupaten/kota dengan tingkat kemiskinan tinggi cenderung berdekatan dengan kabupaten/kota miskin lainnya (kluster kemiskinan).

  2. Pada morans plot terlihat yg menujukan diamond artinya signifikan yaitu Tasikmalaya dan cirebon (Kuadran I) dan ciamis (Kuadran II)

# 9. Peta choropleth kemiskinan (opsional, tapi bagus)
library(tmap)
library(tmaptools)
library(dplyr)# untuk smooth_map (opsional)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.
# Hitung titik centroid untuk tempat menaruh label
jabar_label <- jabar_join %>%
  mutate(centroid = st_centroid(geometry)) %>%          # centroid tiap kab/kota
  st_sf()

tm_shape(jabar_join) +
  tm_polygons(col = "p.miskin16",
              palette = c("#2C7BB6", "#ABD9E9", "#FFFFBF", "#FDAE61", "#F46D43", "#D7191C"),
              style = "quantile",
              n = 6,
              title = "% Penduduk Miskin 2016",
              border.col = "grey40",
              lwd = 0.7) +
  
  # Tambahkan label nama kab/kota di centroid
  tm_shape(jabar_label) +
  tm_text(text = "KABKOT",          # ganti dengan kolom yang berisi nama kabkot
          size = 0.4,               # ukuran teks (sesuaikan)
          col = "black", 
          fontface = "bold",
          shadow = TRUE,            # kasih bayangan biar lebih jelas
          bg.color = "white",       # background putih biar kontras
          bg.alpha = 0.7,           # sedikit transparan
          auto.placement = TRUE,    # otomatis hindari tumpang tindih (penting!)
          remove.overlap = TRUE) +  # kalau masih bertabrakan, hapus yang kecil
  
  tm_layout(main.title = "Sebaran Persentase Penduduk Miskin Jawa Barat 2016",
            main.title.size = 1.2,
            legend.position = c("left", "bottom"),
            legend.title.size = 0.5,
            legend.text.size = 0.5,
            frame = FALSE) +
  
  tm_credits("Sumber: BPS Jawa Barat 2016", 
             position = c("RIGHT", "BOTTOM"), size = 0.7,fontface = "bold",)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[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_text()`: migrate the layer options 'shadow', 'auto.placement'
## (rename to 'point.label') to 'options = opt_tm_text(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

2. Regresi Spasial

library(sf)
library(spdep)
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.5.2
## 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(readxl)
library(nortest)     # ad.test
library(lmtest)      # bptest, dwtest, dll
## Warning: package 'lmtest' was built under R version 4.5.2
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)         # untuk beberapa uji (opsional)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(DescTools)   # RunsTest
## Warning: package 'DescTools' was built under R version 4.5.2
## 
## Attaching package: 'DescTools'
## The following object is masked from 'package:car':
## 
##     Recode
library(arules)      # discretize
## Warning: package 'arules' was built under R version 4.5.2
## 
## Attaching package: 'arules'
## The following object is masked from 'package:car':
## 
##     recode
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following objects are masked from 'package:base':
## 
##     abbreviate, write
library(tmap)        # untuk plot bagus (opsional)
# ----------------------------------------------------
# 1. Baca dan gabungkan data (pastikan sudah benar seperti sebelumnya)
# ----------------------------------------------------
jabar <- st_read("D:/Kuliah/IPB 2025 Semester 3/Analisis Spasial/Praktikum/Sesi Uas/Praktikum 9/Regresi Spasial/petaJabar2/petaJabar2/Jabar2.shp")
## Reading layer `Jabar2' from data source 
##   `D:\Kuliah\IPB 2025 Semester 3\Analisis Spasial\Praktikum\Sesi Uas\Praktikum 9\Regresi Spasial\petaJabar2\petaJabar2\Jabar2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS:           NA
data.jabar <- read_excel("D:/Kuliah/IPB 2025 Semester 3/Analisis Spasial/Praktikum/Sesi Uas/Praktikum 9/Regresi Spasial/Jabar Data.xlsx", sheet = "data")

# Gabungkan berdasarkan kolom kunci yang benar (ganti "KABKOT" jika beda)
jabar_join <- merge(jabar, data.jabar, by = "KABKOT", all.x = TRUE)

# Konversi ke Spatial untuk spdep
jabar.sp <- as(jabar_join, "Spatial")
# Buat matriks tetangga dan weights (Queen + Row-standardized)
nb <- poly2nb(jabar.sp, queen = TRUE)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# 3. Regresi Klasik (OLS)
# ----------------------------------------------------
ols <- lm(p.miskin16 ~ EYS2016, data = jabar_join)
summary(ols)
## 
## Call:
## lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1785 -2.0715 -0.5418  1.5651  7.9577 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  35.8463     8.4317   4.251 0.000144 ***
## EYS2016      -2.1051     0.6649  -3.166 0.003140 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.963 on 36 degrees of freedom
## Multiple R-squared:  0.2178, Adjusted R-squared:  0.1961 
## F-statistic: 10.02 on 1 and 36 DF,  p-value: 0.00314
err_ols <- residuals(ols)
# Uji Asumsi Regresi Klasik
ad.test(err_ols)                    # Normalitas (Anderson-Darling)
## 
##  Anderson-Darling normality test
## 
## data:  err_ols
## A = 0.72137, p-value = 0.0549
hist(err_ols, main = "Histogram Residual OLS", col = "lightblue", breaks = 15)

qqnorm(err_ols); qqline(err_ols, col = "red", lwd = 2)

durbinWatsonTest(ols)               # Autokorelasi temporal (kurang relevan di spasial)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.3475851      1.282738   0.034
##  Alternative hypothesis: rho != 0
RunsTest(err_ols)                   # Uji acak residual
## 
##  Runs Test for Randomness
## 
## data:  err_ols
## z = -1.1512, runs = 16, m = 19, n = 19, p-value = 0.2496
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
##  median(x) 
## -0.5417598
bptest(ols)                         # Heteroskedastisitas (Breusch-Pagan)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 1.5956, df = 1, p-value = 0.2065
# Uji Autokorelasi Spasial pada Residual OLS (PALING PENTING!)
moran.test(err_ols, listw, randomisation = TRUE, alternative = "greater")
## 
##  Moran I test under randomisation
## 
## data:  err_ols  
## weights: listw    
## 
## Moran I statistic standard deviate = 6.4291, p-value = 6.419e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.553529086      -0.027027027       0.008154371
# ----------------------------------------------------
# 4. Lagrange Multiplier Tests (untuk pilih model terbaik)
# ----------------------------------------------------
LMtests <- lm.LMtests(ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
print(LMtests)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## test weights: listw
## 
## RSerr = 30.544, df = 1, p-value = 3.264e-08
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## test weights: listw
## 
## RSlag = 28.96, df = 1, p-value = 7.388e-08
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## test weights: listw
## 
## adjRSerr = 1.8286, df = 1, p-value = 0.1763
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## test weights: listw
## 
## adjRSlag = 0.24482, df = 1, p-value = 0.6207
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = jabar_join)
## test weights: listw
## 
## SARMA = 30.789, df = 2, p-value = 2.062e-07
# ----------------------------------------------------
# 5. Pemodelan Spasial
# ----------------------------------------------------
library(lmtest)
library(spatialreg)   # pastikan sudah di-load
# A. Spatial Lag Model (SAR / SLM)
sar <- lagsarlm(p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, zero.policy = TRUE)
summary(sar)
## 
## Call:
## lagsarlm(formula = p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.20154 -1.26555 -0.16793  0.78134  5.94654 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 25.38595    5.64592  4.4963 6.913e-06
## EYS2016     -1.84750    0.43429 -4.2540 2.099e-05
## 
## Rho: 0.77615, LR test value: 24.337, p-value: 8.0878e-07
## Asymptotic standard error: 0.09672
##     z-value: 8.0246, p-value: 1.1102e-15
## Wald statistic: 64.395, p-value: 9.992e-16
## 
## Log likelihood: -82.00512 for lag model
## ML residual variance (sigma squared): 3.7246, (sigma: 1.9299)
## Number of observations: 38 
## Number of parameters estimated: 4 
## AIC: 172.01, (AIC for lm: 194.35)
## LM test for residual autocorrelation
## test value: 0.25918, p-value: 0.61068
err.sar<-residuals(sar)
ad.test(err.sar) #uji kenormalan,hrs menggunakan package "nortest"
## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.86538, p-value = 0.02382
bptest(residuals(sar) ~ fitted(sar))
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## 
##  studentized Breusch-Pagan test
## 
## data:  residuals(sar) ~ fitted(sar)
## BP = 0.10831, df = 1, p-value = 0.7421
RunsTest(err.sar)
## 
##  Runs Test for Randomness
## 
## data:  err.sar
## z = 0.16446, runs = 21, m = 19, n = 19, p-value = 0.8694
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x) 
## -0.167926
moran.test(err.sar, listw, randomisation = TRUE)  # Autokorelasi spasial residual (harusnya tidak signifikan)
## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.019044, p-value = 0.4924
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.025342568      -0.027027027       0.007823602
# Impacts SAR (Direct, Indirect, Total)
imp_sar <- impacts(sar, listw = listw)
print(imp_sar, zstats = TRUE)   # Tampilkan impacts lengkap + z-stats
## Impact measures (lag, exact):
##                  Direct  Indirect     Total
## EYS2016 dy/dx -2.298508 -5.954681 -8.253188
# B. Spatial Error Model (SEM)
sem <- errorsarlm(p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, zero.policy = TRUE)
summary(sem)
## 
## Call:errorsarlm(formula = p.miskin16 ~ EYS2016, data = jabar_join, 
##     listw = listw, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.20331 -1.51842 -0.16988  1.22079  5.89772 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 31.26867    5.67174  5.5131 3.526e-08
## EYS2016     -1.76389    0.43366 -4.0674 4.754e-05
## 
## Lambda: 0.77312, LR test value: 23.206, p-value: 1.4556e-06
## Asymptotic standard error: 0.10111
##     z-value: 7.6463, p-value: 2.065e-14
## Wald statistic: 58.466, p-value: 2.065e-14
## 
## Log likelihood: -82.57064 for error model
## ML residual variance (sigma squared): 3.8444, (sigma: 1.9607)
## Number of observations: 38 
## Number of parameters estimated: 4 
## AIC: 173.14, (AIC for lm: 194.35)
err_sem <- residuals(sem)
ad.test(err_sem)
## 
##  Anderson-Darling normality test
## 
## data:  err_sem
## A = 0.6705, p-value = 0.07376
bptest(residuals(sem) ~ fitted(sem))
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## 
##  studentized Breusch-Pagan test
## 
## data:  residuals(sem) ~ fitted(sem)
## BP = 0.32238, df = 1, p-value = 0.5702
RunsTest(err_sem)
## 
##  Runs Test for Randomness
## 
## data:  err_sem
## z = 0.16446, runs = 21, m = 19, n = 19, p-value = 0.8694
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
##  median(x) 
## -0.1698759
moran.test(err_sem, listw)   # harusnya tidak signifikan
## 
##  Moran I test under randomisation
## 
## data:  err_sem  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.29144, p-value = 0.3854
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0008751173     -0.0270270270      0.0080520512
cat("Model SEM: Hanya efek DIRECT dari EYS2016 = ", coef(sem)["EYS2016"], "\n")
## Model SEM: Hanya efek DIRECT dari EYS2016 =  -1.763887
cat("Tidak ada efek indirect (spillover) karena model error.\n")
## Tidak ada efek indirect (spillover) karena model error.
# C. SARAR (SAC = SAR + Error) – sering jadi model terbaik
sarar <- sacsarlm(p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, zero.policy = TRUE)
summary(sarar)
## 
## Call:
## sacsarlm(formula = p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.02036 -1.17876 -0.10730  0.68317  5.88012 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  24.2541     5.5167  4.3965 1.100e-05
## EYS2016      -1.8147     0.4144 -4.3790 1.192e-05
## 
## Rho: 0.85293
## Asymptotic standard error: 0.11164
##     z-value: 7.6402, p-value: 2.176e-14
## Lambda: -0.2798
## Asymptotic standard error: 0.45202
##     z-value: -0.619, p-value: 0.53592
## 
## LR test value: 24.649, p-value: 4.4422e-06
## 
## Log likelihood: -81.84917 for sac model
## ML residual variance (sigma squared): 3.4444, (sigma: 1.8559)
## Number of observations: 38 
## Number of parameters estimated: 5 
## AIC: 173.7, (AIC for lm: 194.35)
err_sarar <- residuals(sarar)
ad.test(err_sarar)
## 
##  Anderson-Darling normality test
## 
## data:  err_sarar
## A = 1.0029, p-value = 0.01073
bptest(residuals(sarar) ~ fitted(sarar))
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## 
##  studentized Breusch-Pagan test
## 
## data:  residuals(sarar) ~ fitted(sarar)
## BP = 0.085311, df = 1, p-value = 0.7702
moran.test(err_sarar, listw)
## 
##  Moran I test under randomisation
## 
## data:  err_sarar  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.35656, p-value = 0.3607
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.004303692      -0.027027027       0.007720965
imp_sarar <- impacts(sarar, listw = listw)
print(imp_sarar, zstats = TRUE)   # JALAN! Ada Direct, Indirect, Total
## Impact measures (sac, exact):
##                  Direct  Indirect     Total
## EYS2016 dy/dx -2.512785 -9.825628 -12.33841
#pemodelan GSM
gsm<-sacsarlm(p.miskin16 ~ EYS2016, data = jabar_join, listw = listw, zero.policy = TRUE)
err.gsm<-residuals(gsm)
ad.test(err.gsm) #uji kenormalan,hrs menggunakan package "nortest"
## 
##  Anderson-Darling normality test
## 
## data:  err.gsm
## A = 1.0029, p-value = 0.01073
bptest(residuals(gsm) ~ fitted(gsm))
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## 
##  studentized Breusch-Pagan test
## 
## data:  residuals(gsm) ~ fitted(gsm)
## BP = 0.085311, df = 1, p-value = 0.7702
    #uji kehomogenan residual
RunsTest(err.gsm)
## 
##  Runs Test for Randomness
## 
## data:  err.gsm
## z = 0.16446, runs = 21, m = 19, n = 19, p-value = 0.8694
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
##  median(x) 
## -0.1072957
imp_gsm <- impacts(gsm, listw = listw)
print(imp_gsm, zstats = TRUE)   # JALAN! Ada Direct, Indirect, Total
## Impact measures (sac, exact):
##                  Direct  Indirect     Total
## EYS2016 dy/dx -2.512785 -9.825628 -12.33841
# ----------------------------------------------------
# 6. Bandingkan Model dengan AIC & LogLik
# ----------------------------------------------------
AIC(ols, sar, sem, sarar,gsm)
##       df      AIC
## ols    3 194.3471
## sar    4 172.0102
## sem    4 173.1413
## sarar  5 173.6983
## gsm    5 173.6983
logLik(ols); logLik(sar); logLik(sem); logLik(sarar); logLik(gsm)
## 'log Lik.' -94.17354 (df=3)
## 'log Lik.' -82.00512 (df=4)
## 'log Lik.' -82.57064 (df=4)
## 'log Lik.' -81.84917 (df=5)
## 'log Lik.' -81.84917 (df=5)
# Model dengan AIC terkecil = model terbaik
# MOdel TErbaik --> sar
# ----------------------------------------------------
# 7. Peta Residual Model Terbaik (contoh: SAR)
# ----------------------------------------------------
jabar_join$resid_sar <- residuals(sar)

tm_shape(jabar_join) +
  tm_polygons(col = "resid_sar",
              palette = "-RdBu",
              style = "quantile",
              n = 7,
              title = "Residual SAR") +
  tm_text("KABKOT", size = 0.5, col = "black") +
  tm_layout(main.title = "Peta Residual Model SAR \n(Kemiskinan Jabar 2016)",
            legend.position = c("left", "bottom"))
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## 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")

Kesimpulan:

“Model regresi spasial terbaik adalah Spatial Lag Model (SAR) dengan AIC = 172,01. Setiap kenaikan satu tahun Harapan Lama Sekolah berpotensi menurunkan angka kemiskinan hingga 8,25% poin secara total (direct + spillover), sehingga kebijakan pendidikan harus bersifat regional dan tidak terbatas pada batas administrasi kabupaten/kota.”

Impact measures (lag, exact):                  Direct  Indirect     Total EYS2016 dy/dx -2.298508 -5.954681 -8.253188