Mata Kuliah: Spatial Statistics

Dosen Pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.

Library

# Initial setup, loading libraries and setting working directory

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(maps)

## Warning: package 'maps' was built under R version 4.2.3

library(raster)

## Warning: package 'raster' was built under R version 4.2.3

## Loading required package: sp

## Warning: package 'sp' was built under R version 4.2.3

library(maptools)

## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: FALSE

## 
## Attaching package: 'maptools'

## The following object is masked from 'package:sp':
## 
##     sp2Mondrian

library(sp)
library(spdep)

## Warning: package 'spdep' was built under R version 4.2.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.2.3

## 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

## Warning: package 'sf' was built under R version 4.2.3

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.2.3

## 
## 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

library(sf)
library(spdep)
library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.2.3

## Loading required package: Matrix

## Warning: package 'Matrix' was built under R version 4.2.3

## 
## 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(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'purrr' was built under R version 4.2.3

## Warning: package 'stringr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()  masks Matrix::expand()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ purrr::map()     masks maps::map()
## ✖ tidyr::pack()    masks Matrix::pack()
## ✖ dplyr::select()  masks raster::select()
## ✖ tidyr::unpack()  masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(RColorBrewer)
library(car)

## Warning: package 'car' was built under R version 4.2.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.2.2

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:maptools':
## 
##     pointLabel

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.2.3

## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

setwd("C:/Users/User/Documents/Stat/sem 5/Spatial Statistics")

Data

Input Data

Data yang digunakan adalah jumlah kejahatan sebagai variabel y, jumlah penduduk sebagai variabel x1, dan rasio jenis kelamin sebagai variabel x2 dari Provinsi Jawa Timur pada tahun 2022.

# Load the kejahatan dataset
kejahatan <- read.csv("data uas.csv")
kejahatan

##    id   Kabupaten.Kota Jumlah.Kejahatan Jumlah.Penduduk Rasio.JK
## 1   1        Bangkalan              976         1086620    97.17
## 2   2       Banyuwangi             2545         1731731   100.20
## 3   3             Batu              413          216735   101.40
## 4   4           Blitar             2238         1240322   101.45
## 5   5       Bojonegoro             1314         1315125   100.81
## 6   6        Bondowoso              947          781417    96.96
## 7   7           Gresik             2184         1332664   101.36
## 8   8           Jember             2773         2567718    99.39
## 9   9          Jombang             2243         1335972   101.63
## 10 10           Kediri             1710         1656020   101.95
## 11 11      Kota Blitar             1217          151960    98.97
## 12 12      Kota Kediri             1047          289418   100.23
## 13 13      Kota Madiun             1132          199192    95.63
## 14 14      Kota Malang             1603          846126    98.98
## 15 15   Kota Mojokerto             1353          134350    98.33
## 16 16    Kota Pasuruan              893          211497   100.01
## 17 17 Kota Probolinggo              886          243200    98.28
## 18 18         Lamongan             1509         1371509   100.11
## 19 19         Lumajang             1378         1137227    97.79
## 20 20           Madiun             1293          757665    98.19
## 21 21          Magetan             1304          678343    96.77
## 22 22           Malang             2010         2685900   101.53
## 23 23        Mojokerto             1315         1133584   101.13
## 24 24          Nganjuk             1404         1117033   101.14
## 25 25            Ngawi             1292          877432    98.53
## 26 26          Pacitan              300          592916   100.38
## 27 27        Pamekasan             1307          857818    96.20
## 28 28         Pasuruan             1480         1619035   100.11
## 29 29         Ponorogo              911          964253    99.76
## 30 30      Probolinggo             1057         1159965    97.33
## 31 31          Sampang              725          984162    98.62
## 32 32         Sidoarjo             3344         2103401   101.31
## 33 33        Situbondo              911          691260    96.36
## 34 34          Sumenep              759         1136632    93.23
## 35 35         Surabaya             8759         2887223    98.27
## 36 36       Trenggalek              709          739669   100.92
## 37 37            Tuban             1208         1209543    99.70
## 38 38      Tulungagung             1787         1105337   100.06

Statistik Deskriptif

Deskriptif data dilakukan untuk melihat karakteristik dari data yang digunakan.

summary(kejahatan)

##        id        Kabupaten.Kota     Jumlah.Kejahatan Jumlah.Penduduk  
##  Min.   : 1.00   Length:38          Min.   : 300.0   Min.   : 134350  
##  1st Qu.:10.25   Class :character   1st Qu.: 954.2   1st Qu.: 703362  
##  Median :19.50   Mode  :character   Median :1305.5   Median :1095979  
##  Mean   :19.50                      Mean   :1585.2   Mean   :1082894  
##  3rd Qu.:28.75                      3rd Qu.:1683.2   3rd Qu.:1328279  
##  Max.   :38.00                      Max.   :8759.0   Max.   :2887223  
##     Rasio.JK     
##  Min.   : 93.23  
##  1st Qu.: 98.21  
##  Median : 99.73  
##  Mean   : 99.22  
##  3rd Qu.:100.89  
##  Max.   :101.95

Peta

Peta Jawa Timur

# Mapping
Indo_Kab<-readRDS('gadm36_IDN_2_sp.rds')
Jatim<-Indo_Kab[Indo_Kab$NAME_1 == "Jawa Timur",]
plot(Jatim)

Peta Persebaran Jumlah Kejahatan

Jatim$id<-c(1:38)
Jatim_sf<-st_as_sf(Jatim)
Jatim_merged <- Jatim_sf %>%
  left_join(kejahatan, by = "id")

Kontinu

ggplot() +
  geom_sf(data=Jatim_merged, aes(fill = Jumlah.Kejahatan),color=NA) +
  theme_bw() +
  scale_fill_gradient(low = "yellow", high = "red") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ 
  theme(legend.position = "right",
        axis.text.x = element_blank(),  # Menghapus label x-axis
        axis.text.y = element_blank())+  # Menghapus label y-axis   
  labs(title = "Peta Persebaran Jumlah Kejahatan Jawa Timur 2022",
       fill = "Jumlah Kejahatan")

Dari peta persebaran secara kontinu, terlihat bahwa terdapat satu kabupaten/kota yang memiliki jumlah kejahatan yang berbeda signifikan dari kabupaten/kota lainnya di Jawa Timur, yaitu Kota Surabaya sebanyak sekitar 8000 jumlah kejahatan.

Diskrit

summary(kejahatan$Jumlah.Kejahatan)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   300.0   954.2  1305.5  1585.2  1683.2  8759.0

breaks <- c(-Inf, 954.2, 1305.5, 1683.2,Inf)

# Memberi label pada setiap interval
labels <- c("Very Low", "Low", "High", "Very High")

# Membuat kolom baru dengan kejahatan diskrit
Jatim_merged$Kejahatan_Discrete <- cut(Jatim_merged$Jumlah.Kejahatan, 
                                  breaks = breaks, labels = labels, 
                                  right = TRUE)

ggplot() +
  geom_sf(data=Jatim_merged, aes(fill = Kejahatan_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(fill = "kejahatan")+theme(legend.position = "right",
                          axis.text.x = element_blank(),  # Menghapus label x-axis
                          axis.text.y = element_blank())+  # Menghapus label y-axis
  labs(title = "Peta Persebaran Jumlah Kejahatan Jawa Timur 2022",
       fill = "Jumlah Kejahatan")

Secara diskrit, terlihat bahwa persebaran dari jumlah kejahatan di kabupaten/kota di Jawa Timur condong menjadi sangat besar. Hal ini kemungkinan terjadi akibat jumlah kejahatan di Kota Surabaya yang jauh lebih besar dari kabupaten/kota lainnya di Jawa Timur.

Spatial Autocorrelation

Matriks Bobot Spasial

row.names(Jatim) <- as.character(1:38)
CoordK<-coordinates(Jatim)

Rook

WR <- poly2nb(Jatim, row.names(Jatim), queen=FALSE)  
WBR <- nb2mat(WR, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"  
WLR<-nb2listw(WR)#List neighbours 
plot(Jatim, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Jatim), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot.nb(WR, CoordK, add=TRUE, col="red", lwd=2)

Queen

WQ <- poly2nb(Jatim, row.names(Jatim), queen=TRUE)  
WBQ <- nb2mat(WQ, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"  
WLQ<-nb2listw(WQ)#List neighbours 
plot(Jatim, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Jatim), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot.nb(WQ, CoordK, add=TRUE, col="red", lwd=2)

KNN

#k=5
WK <- knn2nb(knearneigh(CoordK, k = 5), row.names(Jatim))
WBK <- nb2mat(WK, style='B', 
              zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WLK<-nb2listw(WK)
plot(Jatim, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Jatim), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot.nb(WK, CoordK, add=TRUE, col="red", lwd=2)

Moran’s I

moran.test(kejahatan$Jumlah.Kejahatan, WLR) #Rook

## 
##  Moran I test under randomisation
## 
## data:  kejahatan$Jumlah.Kejahatan  
## weights: WLR    
## 
## Moran I statistic standard deviate = 3.1535, p-value = 0.0008067
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.243582388      -0.027027027       0.007363864

moran.test(kejahatan$Jumlah.Kejahatan, WLQ) #Queen

## 
##  Moran I test under randomisation
## 
## data:  kejahatan$Jumlah.Kejahatan  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 3.1535, p-value = 0.0008067
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.243582388      -0.027027027       0.007363864

moran.test(kejahatan$Jumlah.Kejahatan, WLK) #KNN

## 
##  Moran I test under randomisation
## 
## data:  kejahatan$Jumlah.Kejahatan  
## weights: WLK    
## 
## Moran I statistic standard deviate = 1.0619, p-value = 0.1441
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.034440515      -0.027027027       0.003350628

Dari hasil diatas, didapatkan nilai Moran’s I terbesar oleh matriks bobot Rook dan Queen sehingga dipilih matriks bobot Queen sebagai matriks bobot optimum.

Global Moran’s I

Global_Moran<-moran.test(kejahatan$Jumlah.Kejahatan, WLQ)
Global_Moran

## 
##  Moran I test under randomisation
## 
## data:  kejahatan$Jumlah.Kejahatan  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 3.1535, p-value = 0.0008067
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.243582388      -0.027027027       0.007363864

Hasil uji Moran’s I dari matriks bobot Queen menunjukkan hasil p-value yang signifikan, sehingga terdapat dependensi spasial pada jumlah kejahatan di kabupaten/kota Jawa Timur.

#Membuat plot Moran's I
moran.plot(kejahatan$Jumlah.Kejahatan,WLQ)

Local Moran’s I

Local_Moran <- localmoran(kejahatan$Jumlah.Kejahatan, WLQ)
   
quadr_data <- attr(Local_Moran, "quadr")
mean_values <- quadr_data$mean

Jatim_sf <- st_as_sf(Jatim)
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Jatim_sf$mean_values <- mean_values_char
ggplot() +
  geom_sf(data = Jatim_sf, aes(fill = mean_values)) +  # Fill based on 'mean_values'
  scale_fill_manual(values = c("Low-Low" = "blue", "High-Low" = "green", 
                               "Low-High" = "yellow", "High-High" = "red")) + 
  labs(title = "Local Moran's I Mean Values in Jatim", fill = "Mean Type") +
  theme_minimal()

Dari hasil uji Local Moran’s I, terlihat bahwa pada peta di atas, terdapat lebih banyak kabupaten/kota yang memiliki jumlah kejahatan sedikit dan dikelilingi oleh kabupaten/kota yang memiliki jumlah kejahatan sedikit juga.

Spatial Regression

OLS Regression

kejahatan.ols<-lm(Jumlah.Kejahatan~Jumlah.Penduduk+Rasio.JK, data=kejahatan) 
summary(kejahatan.ols)

## 
## Call:
## lm(formula = Jumlah.Kejahatan ~ Jumlah.Penduduk + Rasio.JK, data = kejahatan)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1749.6  -519.3  -122.1   395.1  4534.8 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      5.381e+03  8.364e+03   0.643    0.524    
## Jumlah.Penduduk  1.434e-03  2.521e-04   5.690 1.98e-06 ***
## Rasio.JK        -5.391e+01  8.498e+01  -0.634    0.530    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 997.5 on 35 degrees of freedom
## Multiple R-squared:  0.487,  Adjusted R-squared:  0.4577 
## F-statistic: 16.61 on 2 and 35 DF,  p-value: 8.462e-06

plot(kejahatan.ols)

Dari hasil analisis OLS, didapatkan bahwa hany variabel x1 yang signifikan. Selain itu, nilai R-square sebesar 0.4577 menunjukkan bahwa sekitar 45.77% varians jumlah kejahatan dapat dijelaskan oleh variabel prediktor. Kemudian dilakukan uji asumsi untuk residual dari model OLS.

Uji Asumsi

err.ols <- residuals(kejahatan.ols) 

shapiro.test(err.ols)           # Saphiro test for normality of residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  err.ols
## W = 0.77175, p-value = 2.92e-06

vif(kejahatan.ols)              # Multicollinearity test

## Jumlah.Penduduk        Rasio.JK 
##        1.076945        1.076945

bptest(kejahatan.ols)           # Breusch-Pagan test for heteroskedasticity

## 
##  studentized Breusch-Pagan test
## 
## data:  kejahatan.ols
## BP = 11.238, df = 2, p-value = 0.003629

durbinWatsonTest(kejahatan.ols) # Non-autocorrelation

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1285231      2.240634    0.45
##  Alternative hypothesis: rho != 0

moran.test(err.ols, WLQ)        # Spatial autocorrelation

## 
##  Moran I test under randomisation
## 
## data:  err.ols  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 1.3377, p-value = 0.09049
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.11634136       -0.02702703        0.01148600

Dari hasil uji asumsi di atas, terlihat bahwa pada uji normalitas, didapatkan p-value sebesar 0.001219 yang artinya data tidak berdistribusi normal. Oleh karena itu, dilakukan transformasi data menggunakan log.

Transformasi Data

# Transformasi log pada variabel 
y_trans = log(kejahatan$Jumlah.Kejahatan) 
x1_trans = log(kejahatan$Jumlah.Penduduk) 
x2_trans = log(kejahatan$Rasio.JK) 
data_trans = data.frame(y_trans, x1_trans,x2_trans)  

# Model regresi dengan variabel yang ditransformasi 
ols_trans <- lm(y_trans ~ x1_trans+x2_trans,                   
                data = data_trans) 
summary(ols_trans)

## 
## Call:
## lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35972 -0.23446 -0.02998  0.26398  1.42784 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.3958    18.2130  -0.406  0.68716    
## x1_trans      0.3972     0.1034   3.843  0.00049 ***
## x2_trans      1.9916     4.0173   0.496  0.62317    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4857 on 35 degrees of freedom
## Multiple R-squared:  0.3222, Adjusted R-squared:  0.2834 
## F-statistic: 8.317 on 2 and 35 DF,  p-value: 0.001108

plot(ols_trans)

Dari hasil OLS pada data yang sudah ditransformasi, didapatkan bahwa hanya variabel x1 juga yang signifikan. Selain itu, nilai R-square sebesar 0.2834 menunjukkan bahwa hanya sekitar 28.34% variasi jumlah kejahatan dapat dijelaskan oleh variabel prediktor. Kemudian, dilakukan kembali uji asumsi pada residual model.

Uji Asumsi

err.trans <- residuals(ols_trans)  

# Uji Asumsi 
shapiro.test(err.trans)     # Saphiro test for normality of residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  err.trans
## W = 0.96211, p-value = 0.2222

vif(ols_trans)              # Multicollinearity test

## x1_trans x2_trans 
## 1.049818 1.049818

bptest(ols_trans)           # Breusch-Pagan test for heteroskedasticity

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_trans
## BP = 0.021839, df = 2, p-value = 0.9891

durbinWatsonTest(ols_trans) # Non-autocorrelation

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1780184        2.3366   0.306
##  Alternative hypothesis: rho != 0

moran.test(err.trans, WLQ)  # Spatial autocorrelation

## 
##  Moran I test under randomisation
## 
## data:  err.trans  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 2.5548, p-value = 0.005312
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.29512746       -0.02702703        0.01590031

Dari hasil uji asumsi diatas, didapatkan bahwa seluruh variabel x memiliki nilai VIF kurang dari 5. Selain itu, didapatkan p-value untuk uji normalitas, heteroskedastisitas, dan autokorelasi lebih besar dari alpha (0.05) sehingga data sudah berdistribusi normal, homogen, dan tidak berkorelasi. Namun, nilai p-value dari uji autokorelasi spasial sebesar 0.005312 yang dibawah 0.05 menandakan bahwa terdapat autokorelasi dari kabupaten/kota di Jawa Timur secara spasial yang mendukung hasil uji Global Moran’s I, sehingga dilakukan analisis yang mempertimbangkan efek spasial.

LM Test

LM<-lm.RStests(ols_trans, WLQ, test="all")
print(LM)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## test weights: WLQ
## 
## RSerr = 4.6582, df = 1, p-value = 0.03091
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## test weights: WLQ
## 
## RSlag = 7.3295, df = 1, p-value = 0.006783
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## test weights: WLQ
## 
## adjRSerr = 1.8104, df = 1, p-value = 0.1785
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## test weights: WLQ
## 
## adjRSlag = 4.4817, df = 1, p-value = 0.03426
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans)
## test weights: WLQ
## 
## SARMA = 9.1398, df = 2, p-value = 0.01036

Dari hasil uji di atas, didapatkan bahwa hasil uji yang signifikan adalah LM Lag, LM Error, Robust LM Lag, dan SARMA yang menunjukkan terdapat dependensi spasial pada lag dan error. Model regresi spasial ditentukan dari hasil uji LM yang signifikan, dan pada hasil diatas meskipun LM error signifikan, Robust LM Error tidak signifikan yang menunjukkan tidak terdapat dependensi spasial kemungkinan besar berasal dari lag. Sehingga, model yang paling baik digunakan adalah Spatial Autoregressive Model(SAR).

Spatial Lag Model (SAR)

sar.kejahatan<-lagsarlm(y_trans~x1_trans+x2_trans, data=data_trans, WLQ)
summary(sar.kejahatan)

## 
## Call:lagsarlm(formula = y_trans ~ x1_trans + x2_trans, data = data_trans, 
##     listw = WLQ)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.019944 -0.213071 -0.014215  0.267218  1.052603 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.865839  14.496234 -0.0597    0.9524
## x1_trans     0.379005   0.082332  4.6034 4.157e-06
## x2_trans    -0.278895   3.206824 -0.0870    0.9307
## 
## Rho: 0.5736, LR test value: 10.05, p-value: 0.0015239
## Asymptotic standard error: 0.12486
##     z-value: 4.5941, p-value: 4.3465e-06
## Wald statistic: 21.106, p-value: 4.3465e-06
## 
## Log likelihood: -19.88623 for lag model
## ML residual variance (sigma squared): 0.14935, (sigma: 0.38646)
## Number of observations: 38 
## Number of parameters estimated: 5 
## AIC: 49.772, (AIC for lm: 57.822)
## LM test for residual autocorrelation
## test value: 0.046008, p-value: 0.83016

Hasil di atas menunjukkan bahwa dari data yang sudah ditransformasi pada model SAR, jumlah penduduk berpengaruh secara signifikan terhadap jumlah kejahatan pada kabupaten/kota di Provinsi Jawa Timur.

Cek Residual

err.sar <- residuals(sar.kejahatan)

shapiro.test(err.sar)         # Saphiro test for normality of residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  err.sar
## W = 0.97632, p-value = 0.5875

bptest.Sarlm(sar.kejahatan)   # Breusch-Pagan test for heteroskedasticity

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.89212, df = 2, p-value = 0.6401

moran.test(err.sar, WLQ)      # Spatial autocorrelation

## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 0.10804, p-value = 0.457
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.01317811       -0.02702703        0.01643190

Kemudian setelah dilakukan uji asumsi residual pada model SAR, didapatkan p-value pada uji normalitas, heteroskedastisitas, dan autokorelasi spasial berada di atas 0.05 (alpha) sehingga model SAR memenuhi asumsi residual dan sudah baik untuk digunakan.

# Menambahkan residual OLS dan SAR ke data
kejahatan$kejahatan.ols.res<-resid(ols_trans)     #Residuals OLS
kejahatan$kejahatan.sar.res<-resid(sar.kejahatan) #Residual SAR

# Gabungkan data residual dengan SpatialPolygonsDataFrame
Jatim@data<-kejahatan

# Visualisasi residual OLS
spplot(Jatim, "kejahatan.ols.res", 
       at = seq(min(Jatim$kejahatan.ols.res, na.rm = TRUE), 
                max(Jatim$kejahatan.ols.res, na.rm = TRUE), 
                length = 12),
       col.regions = rev(brewer.pal(11, "RdBu")),
       main = "Residual Model OLS")

# Visualisasi residual SAR
spplot(Jatim, "kejahatan.sar.res", 
       at = seq(min(Jatim$kejahatan.sar.res, na.rm = TRUE), 
                max(Jatim$kejahatan.sar.res, na.rm = TRUE), 
                length = 12),
       col.regions = rev(brewer.pal(11, "RdBu")),
       main = "Residual Spatial Lag Model (SAR)")

Menghitung Impact

impacts(sar.kejahatan, listw=WLQ)

## Impact measures (lag, exact):
##              Direct   Indirect      Total
## x1_trans  0.4281708  0.4606828  0.8888536
## x2_trans -0.3150743 -0.3389986 -0.6540730

Dari hasil di atas, didapatkan hasil bahwa setiap peningkatan satu jumlah penduduk (x1) di Jawa Timur, maka jumlah kejahatan meningkat sebesar 0.428 di kabupaten/kota tersebut dan meningkat sebesar 0.461 di kabupaten/kota sekitarnya. Sedangkan setiap peningkatan satu rasio jenis kelamin (x2) di Jawa Timur, maka jumlah kejahatan akan menurun sebesar 0.315 di kabupaten/kota tersebut dan menurun sebesar 0.339 pada kabupaten/kota di sekitarnya.

Pemilihan Model Terbaik

Untuk menentukan model terbaik yang digunakan dapat dilihat dari nilai AIC. Nilai AIC paling kecil menandakan bahwa model tersebut adalah model terbaik. Berikut merupakan nilai AIC dari masing-masing model.

data.frame("Model" = c("OLS","SAR"),
           "AIC" = c(AIC(ols_trans),AIC(sar.kejahatan)))%>% arrange(AIC)

##   Model      AIC
## 1   SAR 49.77247
## 2   OLS 57.82201

Dari hasil di atas, dapat disimpulkan bahwa nilai AIC paling kecil berasal dari model SAR yaitu sebesar 49.77247, sehingga model terbaik yang digunakan pada penelitian ini adalah model Spatial Autoregressive (SAR).

Pemodelan Spatial Autoregressive (SAR) untuk Menentukan Faktor-Faktor yang Mempengaruhi Jumlah Kejahatan di Jawa Timur

Hotlas Septiani Silaen

2024-12-04

Library

Data

Input Data

Statistik Deskriptif

Peta

Peta Jawa Timur

Peta Persebaran Jumlah Kejahatan

Kontinu

Diskrit

Spatial Autocorrelation

Matriks Bobot Spasial

Rook

Queen

KNN

Moran’s I

Global Moran’s I

Local Moran’s I

Spatial Regression

OLS Regression

Uji Asumsi

Transformasi Data

Uji Asumsi

LM Test

Spatial Lag Model (SAR)

Cek Residual

Menghitung Impact

Pemilihan Model Terbaik