Spatial Autoregression Model

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("C:/Users/pf2tv/Documents/DATA") #sesuaikan dengan komputer masing-masing

A. Spatial data

Semarang_fix <- read_excel("Semarang_fix.xlsx")

## New names:
## • `sekolah` -> `sekolah...9`
## • `sekolah` -> `sekolah...13`
## • `sekolah` -> `sekolah...17`

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

## # A tibble: 3,000 × 20
##    judul  alamat kamar rumah lahan  harga   lat   lng sekolah...9 lat_SD long_SD
##    <chr>  <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <chr>        <dbl>   <dbl>
##  1 Rumah… Pedur…     2    80    80 5.25e8 -7.01  110. SD NEGERI …  -7.01    110.
##  2 Rumah… Pedur…     3   408   408 3.5 e9 -7.01  110. MIS Atthoh…  -7.01    110.
##  3 Rumah… Mijen…     3   240   240 3.20e9 -7.05  110. SD NEGERI …  -7.05    110.
##  4 Rumah… Banyu…     3    84    84 6.47e8 -7.08  110. SD NEGERI …  -7.07    110.
##  5 DE VI… Kedun…     3    65    65 8.11e8 -7.03  110. SD NEGERI …  -7.03    110.
##  6 RUMAH… Jatin…     4   105   105 1.83e9 -7.03  110. SD NEGERI …  -7.03    110.
##  7 STRAT… Semar…    10   219   219 3.60e9 -7.00  110. SD NEGERI …  -7.00    110.
##  8 STRAT… Semar…    10   120   120 2.65e9 -6.98  110. SD NEGERI …  -6.98    110.
##  9 RUMAH… Pedur…     3   195   195 1.10e9 -7.00  110. SD KEMALA …  -7.00    110.
## 10 DIJUA… Mijen…     2   120   120 5.50e8 -6.97  110. SD NEGERI …  -6.97    110.
## # ℹ 2,990 more rows
## # ℹ 9 more variables: dis_SD <dbl>, sekolah...13 <chr>, lat_SMP <dbl>,
## #   long_SMP <dbl>, dis_SMP <dbl>, sekolah...17 <chr>, lat_SMA <dbl>,
## #   long_SMA <dbl>, dis_SMA <dbl>

Semarang_fix$lng <- as.numeric(Semarang_fix$lng)
Semarang_fix$lat <- as.numeric(Semarang_fix$lat)
Spat.data<-Semarang_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.0000000  0.3603791  0.3603791 -0.1097817 -0.1099915 -0.1058580  0.6357817
## V2  0.3603791  1.0000000  1.0000000 -0.1847086 -0.1886411 -0.1864979  0.8000423
## V3  0.3603791  1.0000000  1.0000000 -0.1847086 -0.1886411 -0.1864979  0.8000423
## V4 -0.1097817 -0.1847086 -0.1847086  1.0000000  0.9999257  0.9997450 -0.1713952
## V5 -0.1099915 -0.1886411 -0.1886411  0.9999257  1.0000000  0.9997580 -0.1763965
## V6 -0.1058580 -0.1864979 -0.1864979  0.9997450  0.9997580  1.0000000 -0.1738149
## V7  0.6357817  0.8000423  0.8000423 -0.1713952 -0.1763965 -0.1738149  1.0000000

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)

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

## 
## Call:
## lm(formula = reg, data = Spat.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.40121 -0.27088  0.07124  0.35131  1.34927 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            15.07499    0.09749 154.630  < 2e-16 ***
## log(Spat.data$kamar)    0.50193    0.01574  31.897  < 2e-16 ***
## log(Spat.data$rumah)    0.98598    0.02461  40.060  < 2e-16 ***
## log(Spat.data$dis_SD)   0.27728    0.01532  18.096  < 2e-16 ***
## log(Spat.data$dis_SMP) -0.34170    0.01898 -18.007  < 2e-16 ***
## log(Spat.data$dis_SMA) -0.02904    0.01023  -2.839  0.00455 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4975 on 2994 degrees of freedom
## Multiple R-squared:  0.7337, Adjusted R-squared:  0.7332 
## F-statistic:  1650 on 5 and 2994 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

queen.nb <- knn2nb(knea)

## Warning in knn2nb(knea): neighbour object has 21 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) 
##               1.531817               1.935913               4.918861 
## log(Spat.data$dis_SMP) log(Spat.data$dis_SMA) 
##               6.763110               2.174723

C.3 Moran Test

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 = 85.509, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.9274817855    -0.0018269507     0.0001181119

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 = 7220.9, df = 1, p-value < 2.2e-16
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## RSlag = 3001.8, df = 1, p-value < 2.2e-16
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## adjRSerr = 4245.3, df = 1, p-value < 2.2e-16
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## adjRSlag = 26.118, df = 1, p-value = 3.212e-07
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg, data = Spat.data)
## test weights: listw1
## 
## SARMA = 7247.1, df = 2, p-value < 2.2e-16

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.28929 -0.19928  0.05078  0.18274  0.91146 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value  Pr(>|z|)
## (Intercept)             5.4573932  0.1492223 36.5722 < 2.2e-16
## log(Spat.data$kamar)    0.3511931  0.0105843 33.1805 < 2.2e-16
## log(Spat.data$rumah)    0.5785011  0.0172016 33.6306 < 2.2e-16
## log(Spat.data$dis_SD)   0.1192800  0.0097325 12.2559 < 2.2e-16
## log(Spat.data$dis_SMP) -0.1004224  0.0120969 -8.3015 < 2.2e-16
## log(Spat.data$dis_SMA) -0.0273541  0.0063028 -4.3400 1.425e-05
## 
## Rho: 0.58354, LR test value: 2680.4, p-value: < 2.22e-16
## Asymptotic standard error: 0.0088834
##     z-value: 65.689, p-value: < 2.22e-16
## Wald statistic: 4315, p-value: < 2.22e-16
## 
## Log likelihood: -818.8653 for lag model
## ML residual variance (sigma squared): 0.093915, (sigma: 0.30646)
## Number of observations: 3000 
## Number of parameters estimated: 8 
## AIC: 1653.7, (AIC for lm: 4332.1)
## LM test for residual autocorrelation
## test value: 2280.3, p-value: < 2.22e-16

#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.487e+01   1.058e-01   1.406e+02   0.000e+00
## log.Spat.data.kamar.         6.655e-01   4.731e-02   1.407e+01   1.462e-43
## log.Spat.data.rumah.         6.060e-01   7.675e-02   7.896e+00   4.012e-15
## log.Spat.data.dis_SD.       -1.018e-01   1.139e-01  -8.937e-01   3.715e-01
## log.Spat.data.dis_SMP.       2.780e-02   1.724e-01   1.612e-01   8.719e-01
## log.Spat.data.dis_SMA.       1.218e-02   1.669e-01   7.302e-02   9.418e-01
## lag.log.Spat.data.kamar.    -1.448e-01   4.944e-02  -2.930e+00   3.417e-03
## lag.log.Spat.data.rumah.     4.268e-01   7.977e-02   5.350e+00   9.478e-08
## lag.log.Spat.data.dis_SD.    3.891e-01   1.166e-01   3.338e+00   8.547e-04
## lag.log.Spat.data.dis_SMP.  -3.673e-01   1.745e-01  -2.106e+00   3.532e-02
## lag.log.Spat.data.dis_SMA.  -4.514e-02   1.663e-01  -2.715e-01   7.860e-01

F. Spatial Error Model (SEM)

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 
## -0.5404110 -0.0259135 -0.0042324  0.0401417  0.4985802 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                         Estimate Std. Error z value  Pr(>|z|)
## (Intercept)            16.688703   0.176206 94.7115 < 2.2e-16
## log(Spat.data$kamar)    0.681557   0.012496 54.5410 < 2.2e-16
## log(Spat.data$rumah)    0.544209   0.020257 26.8657 < 2.2e-16
## log(Spat.data$dis_SD)   0.031528   0.025057  1.2583 0.2083010
## log(Spat.data$dis_SMP) -0.118563   0.034819 -3.4051 0.0006613
## log(Spat.data$dis_SMA) -0.043190   0.031256 -1.3818 0.1670286
## 
## Lambda: 0.93982, LR test value: 6962.3, p-value: < 2.22e-16
## Asymptotic standard error: 0.0036504
##     z-value: 257.46, p-value: < 2.22e-16
## Wald statistic: 66284, p-value: < 2.22e-16
## 
## Log likelihood: 1322.059 for error model
## ML residual variance (sigma squared): 0.01749, (sigma: 0.13225)
## Number of observations: 3000 
## Number of parameters estimated: 8 
## AIC: -2628.1, (AIC for lm: 4332.1)

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 
## -0.5404110 -0.0259135 -0.0042324  0.0401417  0.4985802 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                         Estimate Std. Error z value  Pr(>|z|)
## (Intercept)            16.688703   0.176206 94.7115 < 2.2e-16
## log(Spat.data$kamar)    0.681557   0.012496 54.5410 < 2.2e-16
## log(Spat.data$rumah)    0.544209   0.020257 26.8657 < 2.2e-16
## log(Spat.data$dis_SD)   0.031528   0.025057  1.2583 0.2083010
## log(Spat.data$dis_SMP) -0.118563   0.034819 -3.4051 0.0006613
## log(Spat.data$dis_SMA) -0.043190   0.031256 -1.3818 0.1670286
## 
## Lambda: 0.93982, LR test value: 6962.3, p-value: < 2.22e-16
## Asymptotic standard error: 0.0036504
##     z-value: 257.46, p-value: < 2.22e-16
## Wald statistic: 66284, p-value: < 2.22e-16
## 
## Log likelihood: 1322.059 for error model
## ML residual variance (sigma squared): 0.01749, (sigma: 0.13225)
## Number of observations: 3000 
## Number of parameters estimated: 8 
## AIC: -2628.1, (AIC for lm: 4332.1)