Model Spasial Global

Tahapan Regresi Spasial

  1. Eksplorasi Data
  2. Regresi Klasik & Uji Asumsi
  3. Matriks Pembobot Spasial
  4. Uji Lagrange Multiplier
  5. Regresi Spasial & Uji Asumsi
  6. Kebaikan Model

Data Import

Silahkan download files yang ada pada link berikut ini:

https://github.com/raoy/SpatialReg

Data yang Anda download terdiri dari dua jenis data:

  1. Data polygon (peta Jawa Barat, dengan extension .shp)
  2. Dataframe (data pendidikan dan kemiskinan, diperoleh dari BPS)
library(openxlsx)
data.jabar = read.xlsx("Jabar Data (gabung).xlsx")
head(data.jabar)

Studi Kasus: Kemiskinan di Jawa Barat

  • Y : Persentase Penduduk Miskin Tahun 2016
  • X : Angka Melek Huruf Tahun 2016
library(spdep)
library(rgdal)
library(raster)

jabar2<-readOGR(dsn=“directory tempat folder utk file .shp", layer=“nama file shp")

Eksplorasi Data

Relationship of the data

plot(data.jabar$EYS2016, data.jabar$p.miskin16,
  xlab="Angka Melek Huruf Thn.2016", 
  ylab="Persentase Penduduk Miskin Thn.2016",
  pch=20, col="orange", cex=2)
reg.klasik = lm(p.miskin16~EYS2016, data = data.jabar)
library(DescTools)
lines.lm(reg.klasik, col=2, add=T)

Plot tersebut memperlihatkan adanya pola hubungan linear negatif antara angka melek huruf terhadap persentase penduduk miskin di Jawa Barat pada tahun 2016.

Plot Persentase Penduduk Miskin Tahun 2016

k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
library(sp)
jabar2$miskin2<- data.jabar$p.miskin16
spplot(jabar2, "miskin2", col.regions=color)

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol pada data persentase kemiskinan di kabupaten/kota di Jawa Barat. Hal ini tampak dari gradasi warna yang cenderung mengumpul, seperti pada warna merah dan oranye.

Moran Test

w<-poly2nb(jabar2)
ww<-nb2listw(w)
moran(data.jabar$p.miskin16, ww, n=length(ww$neighbours), 
      S0=Szero(ww))
## $I
## [1] 0.3932657
## 
## $K
## [1] 2.403804
moran.test(data.jabar$p.miskin16, ww,randomisation=T, 
           alternative="greater")
## 
##  Moran I test under randomisation
## 
## data:  data.jabar$p.miskin16  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.0168, p-value = 0.001277
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.3932657        -0.0400000         0.0206265

Moran Plot

moran.plot(data.jabar$p.miskin16, ww, labels=data.jabar$KABKOT)

Classical Regression modelling

OLS Method

reg.klasik = lm(p.miskin16~EYS2016, data = data.jabar)
err.regklasik<-residuals(reg.klasik)
summary(reg.klasik)
## 
## Call:
## lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3627 -1.9604 -0.3899  1.5238  8.1125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.5732     9.0973   4.350 0.000217 ***
## EYS2016      -2.3948     0.7209  -3.322 0.002854 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.788 on 24 degrees of freedom
## Multiple R-squared:  0.315,  Adjusted R-squared:  0.2865 
## F-statistic: 11.04 on 1 and 24 DF,  p-value: 0.002854

Model Diagnostics

Normality Test

library(nortest)
library(car)
library(DescTools)
library(lmtest)
ad.test(err.regklasik)
## 
##  Anderson-Darling normality test
## 
## data:  err.regklasik
## A = 0.39412, p-value = 0.3495
hist(err.regklasik)

qqnorm(err.regklasik,datax=T)
qqline(rnorm(length(err.regklasik),mean(err.regklasik),sd(err.regklasik)),datax=T, col="red")

H0: galat model menyebar normal

H1: galat model tidak menyebar normal

Heteroscedastics

bptest(reg.klasik)
## 
##  studentized Breusch-Pagan test
## 
## data:  reg.klasik
## BP = 1.0644, df = 1, p-value = 0.3022

H0: ragam galat homogen

H1: ragam galat tidak homogen

Explore for Spatial Autocorrelation

Uji kebebasan sisaan pada data spasial dapat dilakukan dengan uji moran menggunakan fungsi berikut.

w<-poly2nb(jabar2)
ww<-nb2listw(w)
lm.morantest(reg.klasik, ww, alternative="two.sided")
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## weights: ww
## 
## Moran I statistic standard deviate = 3.4736, p-value = 0.0005135
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.44799554      -0.05028275       0.02057696

Selain menggunakan fungsi lm.morantest, uji moran dapat dilakukan menggunakan fungsi moran.test seperti yang dibahas pada modul pertemuan sebelumnya. Perbedaannya adalah pada fungsi pertama, input yang digunakan adalah objek lm, sedangkan pada fungsi kedua, yang digunakan sebagai input adalah data sisaan model.

moran.test(err.regklasik, ww,randomisation=F, alternative="two.sided")
## 
##  Moran I test under normality
## 
## data:  err.regklasik  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.4266, p-value = 0.0006112
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.44799554       -0.04000000        0.02028183

Terlihat pada output bahwa hasil kedua tes menunjukkan kesimpulan yang sama, yaitu tolak H0 yang menyatakan bahwa tidak terdapat autokorelasi pada sisaan model regresi klasik pada taraf nyata 5%. Oleh karenanya, untuk mencari model yang lebih baik, kita dapat melakukan uji LM (lagrange multiplier) untuk mengidentifikasi model dependensi spasial yang dapat digunakan pada kasus ini.

Lagrange Multiplier Test

LM<-lm.LMtests(reg.klasik, nb2listw(w, style="W"),
               test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
summary(LM)
##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## weights: nb2listw(w, style = "W")
##  
##        statistic parameter  p.value   
## LMerr   8.375484         1 0.003803 **
## LMlag   6.642323         1 0.009958 **
## RLMerr  1.822499         1 0.177016   
## RLMlag  0.089338         1 0.765020   
## SARMA   8.464822         2 0.014517 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Output memperlihatkan bahwa hasil uji model SEM dan SAR sama-sama signifikan pada taraf 5%. Selanjutnya, hasil uji robust keduanya ternyata sama-sama tidak signifikan. Berdasarkan skema tersebut, kita dapat mencoba kandidat model SARMA atau GSM. Namun demikian, ada pula pendapat yang menyarankan agar kita mengambil kandidat model dengan p-value terkecil, dalam hal ini adalah model SEM ( p-value = 0.003803 ).

Pada modul ini, untuk kepentingan pembelajaran, kita akan mencoba ketiga model, SEM, SAR, dan SARMA, meskipun pada prakteknya, Anda hanya perlu memodelkan yang menurut Anda terbaik saja.

Spatial Regression Modelling

Model SEM

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

Output di atas menunjukkan bahwa koefisien Lambda signifikan pada taraf nyata 5% ( p-value = 0.0030663 ). AIC model SEM adalah sebesar 124.25. Selanjutnya kita akan coba memeriksa sisaan model SEM ini.

err.sem<-residuals(sem)
ad.test(err.sem)
## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.63011, p-value = 0.08968
bptest.sarlm(sem)
## Warning: Method bptest.sarlm moved to the spatialreg package
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.38918, df = 1, p-value = 0.5327
moran.test(err.sem, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.096995, p-value = 0.9227
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.05288671       -0.04000000        0.01765155

Terlihat pada output di atas bahwa sisaan telah memenuhi asumsi kenormalan, kehomogenan ragam, serta kebebasan.

Model SAR

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

Output di atas memperlihatkan bahwa koefisien Rho pada model SAR signifikan, dengan nilai AIC sebesar 125.08. Selain itu, terlihat pula hasil uji autokorelasi pada sisaan model memperlihatkan nilai p-value sebesar 0.8481, artinya tidak terdapat autokorelasi pada sisaan.

err.sar<-residuals(sar)
ad.test(err.sar)
## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.70904, p-value = 0.05644
bptest.sarlm(sar)
## Warning: Method bptest.sarlm moved to the spatialreg package
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.97078, df = 1, p-value = 0.3245

Berdasarkan output di atas, pada taraf 5% dapat disimpulkan bahwa sisaan model telah memenuhi asumsi kenormalan dan kehomogenan ragam.

Model GSM/SARMA

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

Output di atas memperlihatkan bahwa hanya salah satu koefisien dependensi spasial yang signifikan, yaitu Lambda. AIC model SARMA adalah sebesar 126.12.

err.gsm<-residuals(gsm)
ad.test(err.gsm)
## 
##  Anderson-Darling normality test
## 
## data:  err.gsm
## A = 0.55307, p-value = 0.1388
bptest.sarlm(gsm)
## Warning: Method bptest.sarlm moved to the spatialreg package
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.31097, df = 1, p-value = 0.5771
moran.test(err.gsm, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.gsm  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.0061095, p-value = 0.9951
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.04081914       -0.04000000        0.01797655

Berdasarkan output di atas, terlihat bahwa sisaan model SARMA telah memenuhi asumsi kenormalan, kehomogenan ragam, dan kebebasan.

Mohon diingat bahwa pada ilustrasi yang kita lakukan saat ini, kita hanya menggunakan satu peubah bebas sehingga kita tidak perlu mengkhawatirkan masalah multikolinieritas. Pada saat Anda memiliki lebih dari satu peubah bebas, pastikan Anda juga memperhatikan multikolinieritas pada model. Pemeriksaan dapat dilakukan dengan fungsi vif() pada package car.

Goodness of Fits

Akhirnya, kita akan coba merangkum hasil pemodelan yang telah dilakukan sepanjang ilustrasi pada modul ini.

Rincian OLS SEM SAR SARMA
AIC 131.020 124.250 125.080 126.120
p-value of Rho NA NA 0.005 0.706
p-value of Lambda NA 0.003 NA 0.003
Kenormalan(p-value) 0.349 0.090 0.056 0.139
Homoskedastisitas (p-value) 0.302 0.533 0.325 0.577
Kebebasan Sisaan(p-value) 0.001 0.923 0.848 0.995

Ilustrasi pada kasus ini memperlihatkan bahwa ternyata SEM merupakan model terbaik berdasarkan nilai AIC-nya. Hal ini ternyata sejalan dengan p-value nya yang juga terkecil pada uji LM.

Exercise

Sebagai latihan, silahkan Anda mencoba memodelkan dengan peubah lain yang terdapat pada data Jabar ini. Diskusikan hasilnya bersama rekan Anda. Anda juga dapat mencoba memodelkan lebih dari satu peubah penjelas.