Spatial Autocorrelation

Further reference: https://rpubs.com/quarcs-lab/tutorial-spatial-regression

library(cluster)  
library(stats)
library(sp)
library(spdep)
library(car)
library(corrplot)
library(lmtest)
library(tidyverse)
library(spatialreg)
library(vtable)
library(readxl)

Working Directory

setwd("D:/Kuliah/Tugas Statistik") #sesuaikan dengan komputer masing-masing

A. Spatial data

Banyumas_fix <- read_excel("Banyumas_fix.xlsx")

## New names:
## • `sekolah` -> `sekolah...6`
## • `sekolah` -> `sekolah...10`
## • `sekolah` -> `sekolah...14`

Banyumas_fix %>% dplyr::distinct(lat, lng, .keep_all = TRUE)

## # A tibble: 107 × 20
##    judul            alamat  harga   lat   lng sekolah...6 lat_SD long_SD  dis_SD
##    <chr>            <chr>   <dbl> <dbl> <dbl> <chr>        <dbl>   <dbl>   <dbl>
##  1 Rumah Dekat UNS… Jl.Ja… 5.76e8 -7.40  109. SDN 1 SUMA…  -7.40    109. 5.62e-3
##  2 Rumah Pusat Kot… Jl.Ve… 3.26e8 -7.43  109. SD MUHAMMA…  -7.42    109. 2.42e-3
##  3 Griya Satria Ba… Jalan… 3.25e8 -7.49  109. SDN 4 KALI…  -7.48    109. 6.10e-3
##  4 Rumah Murah Ban… Jl ra… 2.51e8 -7.48  109. SD NEGERI …  -7.48    109. 5.33e-3
##  5 Griya Satria In… Jalan… 5.36e8 -7.42  109. SD KARITAS   -7.42    109. 1.12e-3
##  6 Jual Perumahan … Jl. R… 2.05e8 -7.48  109. SD NEGERI …  -7.48    109. 1.68e-4
##  7 Jual Rumah Peru… Jl. M… 9.5 e8 -7.44  109. SD NEGERI …  -7.44    109. 1.73e-3
##  8 Jual Rumah Mura… Jl.Ra… 2.51e8 -7.49  109. SD IT NURU…  -7.49    109. 7.82e-3
##  9 Rumah Hook Kawa… Purwe… 6.90e8 -7.42  109. SD MUHAMMA…  -7.42    109. 2.78e-3
## 10 Rumah Bonus Kol… Purwo… 2.95e8 -7.44  109. SD NEGERI …  -7.43    109. 4.64e-3
## # ℹ 97 more rows
## # ℹ 11 more variables: sekolah...10 <chr>, lat_SMP <dbl>, long_SMP <dbl>,
## #   dis_SMP <dbl>, sekolah...14 <chr>, lat_SMA <dbl>, long_SMA <dbl>,
## #   dis_SMA <dbl>, kamar <dbl>, rumah <dbl>, lahan <dbl>

Banyumas_fix$lng <- as.numeric(Banyumas_fix$lng)
Banyumas_fix$lat <- as.numeric(Banyumas_fix$lat)
Spat.data<-Banyumas_fix
correlation<-as.data.frame(cbind( 
       Spat.data$kamar,
       Spat.data$rumah,
       Spat.data$lahan,
       Spat.data$dis_SD,
       Spat.data$dis_SMP,
       Spat.data$dis_SMA,
       Spat.data$harga))
M<-cor(correlation)
M

##            V1         V2         V3         V4         V5         V6         V7
## V1 1.00000000 0.74679743 0.55137424 0.01874103 0.01762895 0.01722088 0.74653718
## V2 0.74679743 1.00000000 0.82343723 0.09608296 0.09686859 0.10030143 0.82919883
## V3 0.55137424 0.82343723 1.00000000 0.04104797 0.04057898 0.04508849 0.73529856
## V4 0.01874103 0.09608296 0.04104797 1.00000000 0.99898026 0.99718266 0.04125475
## V5 0.01762895 0.09686859 0.04057898 0.99898026 1.00000000 0.99850012 0.04193665
## V6 0.01722088 0.10030143 0.04508849 0.99718266 0.99850012 1.00000000 0.04389572
## V7 0.74653718 0.82919883 0.73529856 0.04125475 0.04193665 0.04389572 1.00000000

corrplot(M, method="color")

coordinates(Spat.data) <- ~ lat + lng

plot(Spat.data$lat, Spat.data$lng)

B. Spatial Autocorrelation

reg = log(Spat.data$harga) ~ 
               log(Spat.data$kamar) + 
               log(Spat.data$rumah) +
               #log(Spat.data$lahan) + 
               log(Spat.data$dis_SD) +
               log(Spat.data$dis_SMP) +
               log(Spat.data$dis_SMA)

C. Ordinary least squares (OLS)

fdhhrgerfe

reg.OLS=lm(reg,data=Spat.data)
summary(reg.OLS)

## 
## Call:
## lm(formula = reg, data = Spat.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.48078 -0.27311 -0.01349  0.29328  1.58911 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            16.0900430  0.2352584  68.393   <2e-16 ***
## log(Spat.data$kamar)    0.1408423  0.0714125   1.972   0.0493 *  
## log(Spat.data$rumah)    0.8458456  0.0470287  17.986   <2e-16 ***
## log(Spat.data$dis_SD)  -0.0119304  0.0347144  -0.344   0.7313    
## log(Spat.data$dis_SMP) -0.0262800  0.0323503  -0.812   0.4171    
## log(Spat.data$dis_SMA)  0.0006635  0.0370628   0.018   0.9857    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4641 on 359 degrees of freedom
## Multiple R-squared:  0.7198, Adjusted R-squared:  0.7159 
## F-statistic: 184.5 on 5 and 359 DF,  p-value: < 2.2e-16

C.1 Spatial Weights

knea <- knearneigh(coordinates(Spat.data), longlat = TRUE, k = 5)

## Warning in knearneigh(coordinates(Spat.data), longlat = TRUE, k = 5):
## knearneigh: coordinates are not geographical: longlat argument wrong

## Warning in knearneigh(coordinates(Spat.data), longlat = TRUE, k = 5):
## knearneigh: identical points found

## Warning in knearneigh(coordinates(Spat.data), longlat = TRUE, k = 5):
## knearneigh: kd_tree not available for identical points

queen.nb <- knn2nb(knea)

## Warning in knn2nb(knea): neighbour object has 8 sub-graphs

queen.listw = nb2listw(queen.nb) 
listw1 = queen.listw

C.2 Detecting multicollinearity

car::vif(reg.OLS)

##   log(Spat.data$kamar)   log(Spat.data$rumah)  log(Spat.data$dis_SD) 
##               2.395075               2.420127               2.243739 
## log(Spat.data$dis_SMP) log(Spat.data$dis_SMA) 
##               1.878091               2.494791

C.3 Moran Test

Moran I = 0.086019920

Nilai Moran I yang dihasilhan positif, menunjukkan adanya autokorelasi spasial positif yang berarti adanya kecenderungan tempat-tempat yang berdekatan memiliki nilai yang serupa. Namun, meskipun nilai Moran I yang dihasilkan positif, nilainya medekati nol, yang menunjukkan pola ini tidak terlalu kuat. Sehingga, meskipun terdapat kecenderungan tempat-tempat yang berdekatan memiliki nilai serupa, sebagian besar datanya tersebar acak tanpa adanya pola yang jelas.

P-value = 8.001e-05 = 0.00008001 p-value kurang dari 0.05 (<0.05) menunjukkan hasil yang signifikan.

lm.morantest(reg.OLS,listw1)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## weights: listw1
## 
## Moran I statistic standard deviate = 3.775, p-value = 8.001e-05
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.086019920     -0.009400361      0.000638931

C.4 Lagrange Multiplier (LM) test

lm.RStests(reg.OLS,listw1,test=c(
  "LMerr", 
  "LMlag", 
  "RLMerr", 
  "RLMlag", 
  "SARMA"))

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## RSerr = 9.7218, df = 1, p-value = 0.001821
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## RSlag = 0.15479, df = 1, p-value = 0.694
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## adjRSerr = 11.011, df = 1, p-value = 0.0009058
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## adjRSlag = 1.4439, df = 1, p-value = 0.2295
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## SARMA = 11.166, df = 2, p-value = 0.003762

D. Spatial Autoregressive (SAR)

reg.SAR <- lagsarlm(reg, 
                    data = Spat.data, 
                    listw1)
summary(reg.SAR)

## 
## Call:lagsarlm(formula = reg, data = Spat.data, listw = listw1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.494878 -0.271454 -0.014936  0.296262  1.589421 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            15.7499539  0.9414745 16.7290  < 2e-16
## log(Spat.data$kamar)    0.1438895  0.0710630  2.0248  0.04289
## log(Spat.data$rumah)    0.8435511  0.0473834 17.8027  < 2e-16
## log(Spat.data$dis_SD)  -0.0132904  0.0344685 -0.3856  0.69981
## log(Spat.data$dis_SMP) -0.0269017  0.0324742 -0.8284  0.40744
## log(Spat.data$dis_SMA)  0.0014555  0.0370584  0.0393  0.96867
## 
## Rho: 0.016769, LR test value: 0.14232, p-value: 0.70598
## Asymptotic standard error: 0.045823
##     z-value: 0.36594, p-value: 0.71441
## Wald statistic: 0.13391, p-value: 0.71441
## 
## Log likelihood: -234.6057 for lag model
## ML residual variance (sigma squared): 0.21174, (sigma: 0.46015)
## Number of observations: 365 
## Number of parameters estimated: 8 
## AIC: 485.21, (AIC for lm: 483.35)
## LM test for residual autocorrelation
## test value: 10.573, p-value: 0.0011472

#impacts(reg.SAR, listw = listw1)

#summary(impacts(reg.SAR, listw=listw1, R=500), zstats=TRUE)

E. SLX Spatially Lagged X

reg.SLX <- lmSLX(reg, data = Spat.data, listw1)
summary(reg.SLX)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                             Estimate     Std. Error   t value      Pr(>|t|)   
## (Intercept)                   1.664e+01    4.506e-01    3.693e+01   1.863e-123
## log.Spat.data.kamar.          1.446e-01    7.224e-02    2.002e+00    4.606e-02
## log.Spat.data.rumah.          8.589e-01    4.805e-02    1.788e+01    2.202e-51
## log.Spat.data.dis_SD.        -4.815e-02    6.269e-02   -7.680e-01    4.430e-01
## log.Spat.data.dis_SMP.       -2.327e-02    7.120e-02   -3.269e-01    7.440e-01
## log.Spat.data.dis_SMA.       -1.094e-01    8.509e-02   -1.286e+00    1.994e-01
## lag.log.Spat.data.kamar.     -1.813e-01    2.325e-01   -7.799e-01    4.360e-01
## lag.log.Spat.data.rumah.     -1.637e-03    1.018e-01   -1.608e-02    9.872e-01
## lag.log.Spat.data.dis_SD.     6.801e-02    8.240e-02    8.253e-01    4.097e-01
## lag.log.Spat.data.dis_SMP.    2.814e-02    7.734e-02    3.639e-01    7.162e-01
## lag.log.Spat.data.dis_SMA.    1.177e-01    9.045e-02    1.301e+00    1.941e-01

F. Spatial Error Model (SEM)

Berdasarkan statistik derkiptif residual model yang dihasilkan (Min, 1Q, Median, 3Q, Max), dapat dilihat nilai median mendekati nol yang berarti bahwa residual tersebar secara seimbang disekitar nol.

Pada bagian coefficients dari output model, terdapat variabel-variabel yang digunakan dalam model beserta nilai-nilainya. Salah satu nilainya adalah nilai Pr(>|z|). Nilai Pr(>|z|) adalah nilai p-value untuk uji signifikansi koefisien.

log(Spat.data$kamar) = 0.03633

log(Spat.data$rumah) = < 2e-16 = 0.0000000000000002

log(Spat.data$dis_SD) = 0.61211

log(Spat.data$dis_SMP) = 0.79354

log(Spat.data$dis_SMA) = 0.67844

Dapat dilihat pada keterangan diatas merupakan Nilai Pr(>|z|) dari masing-masing variabel. variabel yang memiliki nilai p-value <0.05 adalah variabel jumlah kamar dan luas rumah, yang berarti kedua variabel ini merupakan variabel yang signifikan, yang menunjukkan bahwa variabel ini memiliki pengaruh penting terhadap variabel dependen.

reg.SEM <- errorsarlm(reg, data = Spat.data, listw1)
summary(reg.SEM)

## 
## Call:errorsarlm(formula = reg, data = Spat.data, listw = listw1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.5222157 -0.2800197 -0.0090969  0.2847529  1.5779226 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            15.9740171  0.2453254 65.1136  < 2e-16
## log(Spat.data$kamar)    0.1480746  0.0707421  2.0932  0.03633
## log(Spat.data$rumah)    0.8599675  0.0467105 18.4106  < 2e-16
## log(Spat.data$dis_SD)  -0.0195793  0.0386127 -0.5071  0.61211
## log(Spat.data$dis_SMP) -0.0098256  0.0375423 -0.2617  0.79354
## log(Spat.data$dis_SMA) -0.0177825  0.0428915 -0.4146  0.67844
## 
## Lambda: 0.19124, LR test value: 6.1233, p-value: 0.013341
## Asymptotic standard error: 0.084593
##     z-value: 2.2607, p-value: 0.023779
## Wald statistic: 5.1107, p-value: 0.023779
## 
## Log likelihood: -231.6152 for error model
## ML residual variance (sigma squared): 0.20765, (sigma: 0.45569)
## Number of observations: 365 
## Number of parameters estimated: 8 
## AIC: 479.23, (AIC for lm: 483.35)

G. Spatial Durbin Model (SDM)

reg.SEM <- errorsarlm(reg, data = Spat.data, listw1)
summary(reg.SEM)

## 
## Call:errorsarlm(formula = reg, data = Spat.data, listw = listw1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.5222157 -0.2800197 -0.0090969  0.2847529  1.5779226 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            15.9740171  0.2453254 65.1136  < 2e-16
## log(Spat.data$kamar)    0.1480746  0.0707421  2.0932  0.03633
## log(Spat.data$rumah)    0.8599675  0.0467105 18.4106  < 2e-16
## log(Spat.data$dis_SD)  -0.0195793  0.0386127 -0.5071  0.61211
## log(Spat.data$dis_SMP) -0.0098256  0.0375423 -0.2617  0.79354
## log(Spat.data$dis_SMA) -0.0177825  0.0428915 -0.4146  0.67844
## 
## Lambda: 0.19124, LR test value: 6.1233, p-value: 0.013341
## Asymptotic standard error: 0.084593
##     z-value: 2.2607, p-value: 0.023779
## Wald statistic: 5.1107, p-value: 0.023779
## 
## Log likelihood: -231.6152 for error model
## ML residual variance (sigma squared): 0.20765, (sigma: 0.45569)
## Number of observations: 365 
## Number of parameters estimated: 8 
## AIC: 479.23, (AIC for lm: 483.35)