Regresi Spasial

Tahapan analisis yang akan digunakan: 1.Eksplorasi Data 2.Ordinary Least Square atau Regresi Klasik dan Uji Asumsi 3.Menentukan Matriks Pembobot Spasial 4.Uji Efek Dependensi Spasial dan Heterogenitas Spasial 5.Regresi Spasial dan Uji Asumsi 6.Kebaikan Model

Library

library(spdep)        
library(sp)           
library(sf)           
library(readxl)       
library(corrplot)     
library(DescTools)    
library(nortest)      
library(car)          
library(spatialreg)   
library(lmtest)       
library(dplyr)
library(RColorBrewer)
library(ggplot2)

menyiapkan library yang akan digunakan selama analisis

Import Data

#Melakukan import data dan menampilkan data yang akan digunakan
dt<-read_xlsx("C:\\Users\\Acer\\Downloads\\UAS REGSPAS\\data-38.xlsx")

#Melakukan import data peta Jawa
jawa<-st_read("C:\\Users\\Acer\\Downloads\\UAS REGSPAS\\jawa.shp")

## Reading layer `jawa' from data source 
##   `C:\Users\Acer\Downloads\UAS REGSPAS\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

1.Eksplorasi Data

par(mfrow=c(1,3))
type <- names(dt)[1:ncol(dt)]
for (i in type){
  boxplot(dt[,i],
       main = paste(i))
}

Berdasarkan hasil eksplorasi sebaran menggunakan boxplot, tidak terdapat pencilan pada setiap peubah.

cor(dt)

##            Y         X1         X2
## Y  1.0000000  0.9256759  0.2431457
## X1 0.9256759  1.0000000 -0.1115304
## X2 0.2431457 -0.1115304  1.0000000

Berdasarkan eksplorasi korelasi terdapat korelasi besar antara Y dengan Peubah X1 yang selnajutnya akan diuji. Diketahui juga bahwa peubah X1 dan X2 memiliki korelasi yang kecil yang menandakan tidak adanya multikolinearitas

cor.test(dt$Y, dt$X1)

## 
##  Pearson's product-moment correlation
## 
## data:  dt$Y and dt$X1
## t = 26.466, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8947636 0.9477583
## sample estimates:
##       cor 
## 0.9256759

cor.test(dt$Y, dt$X2)

## 
##  Pearson's product-moment correlation
## 
## data:  dt$Y and dt$X2
## t = 2.7114, df = 117, p-value = 0.00771
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06604039 0.40539888
## sample estimates:
##       cor 
## 0.2431457

Dari hasil uji kedua peubah penjelas, didapatkan tolak H0 dengan taraf 5% yang berarti diketahui bahwa terdapat korelasi yang signifikan antara Y dengan kedua peubah X

jawa$y <- dt$Y

# Plot menggunakan ggplot2
ggplot(jawa) +
  geom_sf(aes(fill = y)) +
  scale_fill_viridis_c(option = "magma")

  labs(title = "Peta Sebaran Spasial Peubah Respon Y") +
  theme_minimal()

## NULL

Memvisualisasi data peubah Y pada peta

par(mfrow=c(1,2))
jawa$x1 <- dt$X1
jawa$x2 <- dt$X2

# Plot menggunakan ggplot2
ggplot(jawa) +
  geom_sf(aes(fill = x1)) +
  scale_fill_viridis_c(option = "magma")

  labs(title = "Peta Sebaran Spasial Peubah Respon X1") +
  theme_minimal()

## NULL

ggplot(jawa) +
  geom_sf(aes(fill = x2)) +
  scale_fill_viridis_c(option = "magma")

  labs(title = "Peta Sebaran Spasial Peubah Respon X2") +
  theme_minimal()

## NULL

2.Ordinary Least Square atau Regresi Klasik dan Uji Asumsi

reg.ols<-lm(y~x1+x2,data=jawa)
summary(reg.ols)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = jawa)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.38063 -0.18945  0.00517  0.20844  1.04850 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.38040    0.08588   16.07   <2e-16 ***
## x1           7.47094    0.10640   70.21   <2e-16 ***
## x2           2.79758    0.10960   25.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3417 on 116 degrees of freedom
## Multiple R-squared:  0.9784, Adjusted R-squared:  0.978 
## F-statistic:  2624 on 2 and 116 DF,  p-value: < 2.2e-16

Dari pemodelan dengan manggunakan OLS didapatkan model yang memiliki peubah yang signifikan secara keseluruhan, selain itu didapatkan juga nilai R-sq dan adj R-sq yang berada pada kisaran 0.978. Hal tersebut menunjukkan bahwa pemodelan linear dengan menggunakan metode OLS sudah bagus.

Asumsi sisaan pemodelan OLS

e<-residuals(reg.ols)

#Uji Normalitas
ad.test(e)

## 
##  Anderson-Darling normality test
## 
## data:  e
## A = 0.37361, p-value = 0.4121

#Uji Homogenitas
lmtest::bptest(reg.ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg.ols
## BP = 1.7961, df = 2, p-value = 0.4074

#Mutikolinearitas
vif(reg.ols)

##       x1       x2 
## 1.012596 1.012596

Didapatkan dari semua uji sisaan tidak terdapat asumsi sisaan yang tidak terpenuhi

3.Menentukan Matriks Pembobot Spasial

Berdasarkan ketetanggaan

Queen

peta <- st_make_valid(jawa)
sp.peta <- as(peta, "Spatial")


sp.polygons <- SpatialPolygons(sp.peta@polygons)
queen <- poly2nb(sp.peta, queen = T)

queen

## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 522 
## Percentage nonzero weights: 3.68618 
## Average number of links: 4.386555 
## 1 region with no links:
## 1
## 3 disjoint connected subgraphs

W.queen <- nb2listw(queen, style='W',zero.policy=TRUE)
qc = moran.test(dt$Y, W.queen,randomisation=T,
alternative="greater", zero.policy = TRUE)
qc

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.queen  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.95027, p-value = 0.171
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.056029700      -0.008547009       0.004618057

Didapatkan nilai moran dan p-value dari matriks pembobot queen tidak cukup signifikan untuk digunakan sebagai pembobot spasial

Rook

rook <- poly2nb(sp.peta, queen=F)
rook

## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 518 
## Percentage nonzero weights: 3.657934 
## Average number of links: 4.352941 
## 1 region with no links:
## 1
## 3 disjoint connected subgraphs

W.rook <- nb2listw(rook, style='W',zero.policy=TRUE)
qc = moran.test(dt$Y, W.rook,randomisation=T,
alternative="greater", zero.policy = TRUE)
qc

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.rook  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.87516, p-value = 0.1907
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.051052843      -0.008547009       0.004637784

Didapatkan nilai moran dan p-value dari matriks pembobot rook tidak cukup signifikan untuk digunakan sebagai pembobot spasial

Berdasarkan Jarak

centroids <- st_centroid(peta)
longlat <- st_coordinates(centroids)

# Tambahkan koordinat long dan lat ke objek peta
peta$long <- longlat[, 1]  # Kolom X (long)
peta$lat  <- longlat[, 2]  # Kolom Y (lat)
coords <- peta[c("long","lat")]
#class(coords)
koord <- as.data.frame(coords)
jarak<-dist(longlat, method = "euclidean")
m.jarak<-as.matrix(jarak)

K-Means

W.knn<-knn2nb(knearneigh(longlat,k=7,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')

MI.knn <- moran(dt$Y,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(dt$Y,W.knn.s,randomisation = T,alternative = "greater")
mt1

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 1.9159, p-value = 0.02769
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.076934967      -0.008474576       0.001987228

Metode K-Means dengan jumlah k-7 menghasilkan nilai moran dengan p-value yang cukup signifikan untuk digunakan sebagai pembobot spasial

Radial Distance Weight

W.rdw <-dnearneigh(longlat,0,105,longlat=TRUE)
W.rdw.s <- nb2listw(W.rdw,style='W')

MI.rdw <- moran(dt$Y,W.rdw.s,n=length(W.rdw.s$neighbours),S0=Szero(W.rdw.s))
mt2 = moran.test(dt$Y,W.rdw.s,randomisation = T,alternative = "greater")
mt2

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.rdw.s    
## 
## Moran I statistic standard deviate = 0.29537, p-value = 0.3839
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0001566656     -0.0084745763      0.0007930208

Metode Radial Distance Weight menghasilkan nilai moran dengan p-value yang tidak signifikan

Inver Distance Weight

alpha1=1
W.idw <-1/(m.jarak^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')

MI.idw <- moran(dt$Y,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(dt$Y,W.idw.s,randomisation = T,alternative = "greater")
mt3

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 2.3222, p-value = 0.01011
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0396068097     -0.0084745763      0.0004286893

Metode IDW menghasilkan nilai moran dengan p-value yang signifikan, tetapi nilai pembobot moran yang lebih kecil dibandingkan dengan metode K-Means

Eksponential Distance Weight

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

MI.edw <- moran(dt$Y,W.edw.s,n=length(W.edw.s$neighbours),S0=Szero(W.edw.s))
mt5 = moran.test(dt$Y,W.edw.s,randomisation = T,alternative = "greater")
mt5

## 
##  Moran I test under randomisation
## 
## data:  dt$Y  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 1.187, p-value = 0.1176
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0082930975     -0.0084745763      0.0001995308

Tidak didapatkan pembobot moran yang signifikan

Dari tahapan pennetuan pembobot spasial maka akan digunakan pembobot spasial dari metode K-Means dengan nilai pembobot moran sebesar 0.0769.

4.Uji Efek Dependensi Spasial dan Heterogenitas Spasial

Efek Dependensi spasial

ww = W.knn.s
lm.morantest(reg.ols, listw=ww, alternative="two.sided", zero.policy = TRUE)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2, data = jawa)
## weights: ww
## 
## Moran I statistic standard deviate = 10.375, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.452503783     -0.007967164      0.001969803

merr <- moran.test(e, ww,randomisation=T, 
           alternative="two.sided")
merr

## 
##  Moran I test under randomisation
## 
## data:  e  
## weights: ww    
## 
## Moran I statistic standard deviate = 10.466, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.452503783      -0.008474576       0.001939819

Peubah = c("Y", "X1", "X2", "Sisaan")
Indeks_Moran = c(my$estimate[1], mx1$estimate[1], mx2$estimate[1], merr$estimate[1])
p_value = c(my$p.value, mx1$p.value, mx2$p.value, merr$p.value)

df1 = cbind.data.frame(Peubah, Indeks_Moran, p_value)
colnames(df1) <- c("Peubah", "Indeks Moran", "p-value")
df1

##   Peubah Indeks Moran      p-value
## 1      Y   0.07693497 5.537238e-02
## 2     X1   0.01497672 5.992705e-01
## 3     X2  -0.08833915 7.340062e-02
## 4 Sisaan   0.45250378 1.231583e-25

Dari uji indeks moran diketahui bahwa terdapat dipendensi spasial pada Sisaan model OLS saja pada taraf 5%, sedangkan peubah Y,X1 dan X2 tidak terdapat pengaruh spasial yang signifikan pada taraf 5%. Tahapan selanjutnya akan dilakukan pengujian dengan langrange multiple

Uji langrange Multiple

#Uji Lagrange Multiplier
model <- lm.LMtests(reg.ols,listw=ww,zero.policy = TRUE, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

## Please update scripts to use lm.RStests in place of lm.LMtests

## Please update scripts to use lm.RStests in place of lm.LMtests
summary(model)

##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = y ~ x1 + x2, data = jawa)
## test weights: listw
##  
##          statistic parameter   p.value    
## RSerr       96.457         1 < 2.2e-16 ***
## adjRSerr    67.745         1 2.220e-16 ***
## RSlag       82.533         1 < 2.2e-16 ***
## adjRSlag    53.821         1 2.196e-13 ***
## SARMA      150.278         2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Output menunjukkan bahwa hasil uji LM model SAR, SEM, dan GSM serta uji robust LM SAR dan LM-SEM signifikan pada taraf 5%. Berdasarkan skema tersebut, maka kita dapat mencoba kandidat model SAR, SEM, dan SARMA.

Efek heterogenitas Spasial

lmtest::bptest(reg.ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg.ols
## BP = 1.7961, df = 2, p-value = 0.4074

Diketahui karena p-value > dari 0.05 maka tidak terdapat heterogenitas spasial pada model

5.Model Regresi Spasial

SAR

sar <- lagsarlm(y ~ x1 + x2, data = jawa, listw = ww)
summary(sar, Nagelkerke = T)

## 
## Call:lagsarlm(formula = y ~ x1 + x2, data = jawa, listw = ww)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.0343691 -0.1275685 -0.0092897  0.1324040  0.8009187 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value Pr(>|z|)
## (Intercept) -0.141807   0.150603  -0.9416   0.3464
## x1           7.393296   0.073507 100.5791   <2e-16
## x2           2.873339   0.075694  37.9601   <2e-16
## 
## Rho: 0.23428, LR test value: 85.167, p-value: < 2.22e-16
## Asymptotic standard error: 0.021324
##     z-value: 10.987, p-value: < 2.22e-16
## Wald statistic: 120.71, p-value: < 2.22e-16
## 
## Log likelihood: 3.029736 for lag model
## ML residual variance (sigma squared): 0.055275, (sigma: 0.23511)
## Nagelkerke pseudo-R-squared: 0.98943 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 3.9405, (AIC for lm: 87.107)
## LM test for residual autocorrelation
## test value: 4.0226, p-value: 0.044893

Asumsi Model SAR

#normalitas
err.sar<-residuals(sar)
ad.test(err.sar)

## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 1.0332, p-value = 0.009825

#Homogenitas
bptest.Sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.33291, df = 2, p-value = 0.8467

#Kebebasan Sisaan
moran.test(err.sar, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: ww    
## 
## Moran I statistic standard deviate = 2.2882, p-value = 0.02212
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.091563026      -0.008474576       0.001911302

Dari uji asumsi sar, asumsi normalitas dan homogenitas tidak ada yang memenuhi pada taraf 5% tetapi tidak terdapat autokorelasi spasial pada sisaan model pada taraf 5%.

SEM

sem<-errorsarlm(y ~ x1 + x2,data=jawa, listw= ww)
summary(sem)

## 
## Call:errorsarlm(formula = y ~ x1 + x2, data = jawa, listw = ww)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.1092745 -0.1314719 -0.0085899  0.1158370  0.8189233 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 1.552772   0.114735  13.534 < 2.2e-16
## x1          7.189591   0.074150  96.961 < 2.2e-16
## x2          2.727190   0.071325  38.236 < 2.2e-16
## 
## Lambda: 0.77584, LR test value: 62.207, p-value: 3.1086e-15
## Asymptotic standard error: 0.065777
##     z-value: 11.795, p-value: < 2.22e-16
## Wald statistic: 139.12, p-value: < 2.22e-16
## 
## Log likelihood: -8.449967 for error model
## ML residual variance (sigma squared): 0.060483, (sigma: 0.24593)
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 26.9, (AIC for lm: 87.107)

Asumsi Model SEM

#normalitas
err.sem<-residuals(sem)
ad.test(err.sem)

## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 1.7208, p-value = 0.0001959

#Homogenitas
bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.11977, df = 2, p-value = 0.9419

#Kebebasan Sisaan
moran.test(err.sem, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.56273, p-value = 0.5736
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.016074815      -0.008474576       0.001903178

Model SEM memenuhi asumsi sisaan homogenitas dan kebebasan pada taraf 5% namun model SEM tidak memnuhi asumsi normalitas sisaan

SARMA

sarma <- sacsarlm(y ~ x1+x2,data=jawa,ww, zero.policy = TRUE)
summary(sarma)

## 
## Call:
## sacsarlm(formula = y ~ x1 + x2, data = jawa, listw = ww, zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.0598031 -0.1118129 -0.0072247  0.1328563  0.7886158 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value Pr(>|z|)
## (Intercept) -0.012409   0.192859  -0.0643   0.9487
## x1           7.390034   0.072046 102.5741   <2e-16
## x2           2.862881   0.072564  39.4530   <2e-16
## 
## Rho: 0.21516
## Asymptotic standard error: 0.026465
##     z-value: 8.1298, p-value: 4.4409e-16
## Lambda: 0.31269
## Asymptotic standard error: 0.14915
##     z-value: 2.0965, p-value: 0.036041
## 
## LR test value: 88.675, p-value: < 2.22e-16
## 
## Log likelihood: 4.784086 for sac model
## ML residual variance (sigma squared): 0.053071, (sigma: 0.23037)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 2.4318, (AIC for lm: 87.107)

Asumsi Model SARMA

#normalitas
err.sarma<-residuals(sarma)
ad.test(err.sarma)

## 
##  Anderson-Darling normality test
## 
## data:  err.sarma
## A = 1.4839, p-value = 0.0007531

#Homogenitas
bptest.Sarlm(sarma)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.20947, df = 2, p-value = 0.9006

#Kebebasan Sisaan
moran.test(err.sarma, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sarma  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.23122, p-value = 0.8171
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.001602998      -0.008474576       0.001899591

Model SEM memenuhi asumsi sisaan homogenitas dan kebebasan pada taraf 5% namun model SEM tidak memnuhi asumsi normalitas sisaan

6.Menentukan Model Terbaik

df <- data.frame("Model" = c("OLS (Regresi KlasiK)","SAR","SEM","SARMA"),
           "AIC" = c(AIC(reg.ols),AIC(sar), AIC(sem),AIC(sarma)))

df

##                  Model       AIC
## 1 OLS (Regresi KlasiK) 87.107070
## 2                  SAR  3.940528
## 3                  SEM 26.899933
## 4                SARMA  2.431829

Dari hasil perbandingan AIC model terbaik yang didapatkan adalah model SARMA