Konsep Dasar

Regresi spasial merupakan metode regresi untuk tipe data spasial atau data yang memiliki efek lokasi (spatial effect). Spatial effect terdiri dari dua jenis yaitu dependensi spasial dan heterogenitas spasial. Dependensi spasial dapat diartikan bahwa pengamatan pada lokasi i bergantung pada pengamatan lain di lokasi j, j≠i. Sedangkan heterogenitas spasial terjadi akibat adanya efek lokasi random, yaitu perbedaan antara satu lokasi dengan lokasi lainnya. Dasar berkembangnya metode regresi spasial adalah metode regresi linear berganda. Pengembangan itu berdasarkan adanya pengaruh tempat atau spasial pada data yang dianalisis. Menurut Tobler (1979), di dalam Hukum Geografi pertama, segala sesuatu saling berhubungan satu dengan yang lainnya, tetapi sesuatu yang dekat lebih mempunyai pengaruh daripada sesuatu yang jauh.

Model regresi linear dalam bentuk matriks: \[ y = X\beta + \epsilon \] dengan: \[ y = \begin{bmatrix} y_1 \\ y_2 \\ . \\ . \\y_n \end{bmatrix} , \quad X = \begin{bmatrix} 1 & x_{11} & x_{12} & ... & x_{1p} \\ 1 & x_{21} & x_{22} & ... & x_{2p} \\ ... &... &... &... &... \\ 1 & x_{n1} & x_{n2} & ... & x_{np} \end{bmatrix} , \quad \beta = \begin{bmatrix} y_1 \\ y_2 \\ . \\ . \\y_n \end{bmatrix} \] Menurut LeSage (1999), model umum regresi spasial dapat dituliskan sebagai berikut: \[ y \; = \;\rho Wy \; + \; X\beta \; + \; u \; , \quad u \; = \; \lambda Wu \; + \; \epsilon \; , \quad \epsilon \sim N(0, \sigma^2_\epsilon I_n) \] dengan:
\[ \begin{alignat*}{2} y &= \text{vektor variabel respon berukuran } n \times 1 \\ \rho &= \text{koefisien parameter spasial lag dari variabel respon} \\ W &= \text{matriks pembobot spasial berukuran } n \times n \\ X &= \text{matriks variabel prediktor berukuran } n \times (p+1) \\ \beta &= \text{vektor koefisien regresi berukuran } (p+1) \times 1 \\ \lambda &= \text{koefisien parameter spasial error} \\ u &= \text{vektor error yang mempunyai efek spasial} \end{alignat*} \]

Aktifkan Library

Library yang digunakan yaitu tmap, raster, stats, spdep, dan spatialreg.

library(tmap)
library(raster)
library(stats)
library(spdep)
library(spatialreg)
library(sf)
library(spData)

Load Peta dan Data

Import data peta dengan tipe shp dan tabel data pendukung dengan tipe csv. Contoh di bawah ini menggunakan dataset sekunder columbus yang tersedia dalam paket spdep pada software R.

# Load data Columbus
columbus <- st_read(system.file("shapes/columbus.gpkg",
                                package="spData"), quiet=TRUE)
col.gal.nb <- read.gal(system.file("weights/columbus.gal",
                                   package="spData"))
# Buat listw
col.listw <- nb2listw(col.gal.nb, style="W")
# Plot map}
plot(st_geometry(columbus),
     border="black",
     col="gray",
     main="Peta Columbus")

names(columbus)
##  [1] "AREA"       "PERIMETER"  "COLUMBUS_"  "COLUMBUS_I" "POLYID"    
##  [6] "NEIG"       "HOVAL"      "INC"        "CRIME"      "OPEN"      
## [11] "PLUMB"      "DISCBD"     "X"          "Y"          "NSA"       
## [16] "NSB"        "EW"         "CP"         "THOUS"      "NEIGNO"    
## [21] "geom"
head(columbus)
## Simple feature collection with 6 features and 20 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 7.950089 ymin: 12.86109 xmax: 10.1806 ymax: 14.74245
## Projected CRS: Undefined Cartesian SRS with unknown unit
##       AREA PERIMETER COLUMBUS_ COLUMBUS_I POLYID NEIG  HOVAL    INC    CRIME
## 1 0.309441  2.440629         2          5      1    5 80.467 19.531 15.72598
## 2 0.259329  2.236939         3          1      2    1 44.567 21.232 18.80175
## 3 0.192468  2.187547         4          6      3    6 26.350 15.956 30.62678
## 4 0.083841  1.427635         5          2      4    2 33.200  4.477 32.38776
## 5 0.488888  2.997133         6          7      5    7 23.225 11.252 50.73151
## 6 0.283079  2.335634         7          8      6    8 28.750 16.029 26.06666
##       OPEN    PLUMB DISCBD     X     Y NSA NSB EW CP THOUS NEIGNO
## 1 2.850747 0.217155   5.03 38.80 44.07   1   1  1  0  1000   1005
## 2 5.296720 0.320581   4.27 35.62 42.38   1   1  0  0  1000   1001
## 3 4.534649 0.374404   3.89 39.82 41.18   1   1  1  0  1000   1006
## 4 0.394427 1.186944   3.70 36.50 40.52   1   1  0  0  1000   1002
## 5 0.405664 0.624596   2.83 40.01 38.00   1   1  1  0  1000   1007
## 6 0.563075 0.254130   3.78 43.75 39.28   1   1  1  0  1000   1008
##                             geom
## 1 POLYGON ((8.624129 14.23698...
## 2 POLYGON ((8.25279 14.23694,...
## 3 POLYGON ((8.653305 14.00809...
## 4 POLYGON ((8.459499 13.82035...
## 5 POLYGON ((8.685274 13.63952...
## 6 POLYGON ((9.401384 13.5504,...

Membuat Matriks Pembobot

Matriks Pembobot pada model regresi spasial disusun berdasarkan persinggungan antar wilayah.

#MATRIKS PEMBOBOT QUEEN
queen.nb <- poly2nb(as(columbus, "Spatial"), queen = TRUE) #create nb
queen.listw=nb2listw(queen.nb) #convert nb to listw type
queen.columbus= queen.listw
summary(queen.columbus)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49 
## Number of nonzero links: 236 
## Percentage nonzero weights: 9.829238 
## Average number of links: 4.816327 
## Link number distribution:
## 
##  2  3  4  5  6  7  8  9 10 
##  5  9 12  5  9  3  4  1  1 
## 5 least connected regions:
## 1 6 42 46 47 with 2 links
## 1 most connected region:
## 20 with 10 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 49 2401 49 22.75119 203.7091
#Ilustrasi Matriks Pembobot
plot(st_geometry(columbus), border="white",col='gray', sub="Pembobot Queen Pada Peta Columbus")
coords <- st_coordinates(st_centroid(st_geometry(columbus)))
plot(queen.columbus, coords, add = TRUE, col = "red")

Uji Dependensi Spasial

Salah satu cara untuk mengetahui adanya dependensi spasial antar lokasi adalah dengan melakukan uji autokorelasi spasial dengan menggunakan statistik Moran’s I. Autokorelasi spasial adalah taksiran dari korelasi antar nilai amatan yang berkaitan dengan lokasi pada variabel yang sama (Yasin dan Saputra, 2013). Jika terdapat pola sistematik dalam penyebaran sebuah variabel, maka terdapat autokorelasi spasial.

#Moran Test: Pembobot Queen
#Variabel Respon CRIME
moran.test(columbus$CRIME,queen.columbus,randomisation=FALSE)
## 
##  Moran I test under normality
## 
## data:  columbus$CRIME  
## weights: queen.columbus    
## 
## Moran I statistic standard deviate = 5.6303, p-value = 8.994e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.500188557      -0.020833333       0.008563413
moran.plot(columbus$CRIME,queen.columbus)

#Variabel Prediktor HOVAL 
moran.test(columbus$HOVAL,queen.columbus,randomisation=FALSE)
## 
##  Moran I test under normality
## 
## data:  columbus$HOVAL  
## weights: queen.columbus    
## 
## Moran I statistic standard deviate = 2.1713, p-value = 0.01496
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.180093114      -0.020833333       0.008563413
moran.plot(columbus$HOVAL,queen.columbus)

#Variabel Respon INC
moran.test(columbus$INC,queen.columbus,randomisation=FALSE) 
## 
##  Moran I test under normality
## 
## data:  columbus$INC  
## weights: queen.columbus    
## 
## Moran I statistic standard deviate = 4.7165, p-value = 1.199e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.415628778      -0.020833333       0.008563413
moran.plot(columbus$INC,queen.columbus)

Pemodelan Regresi Spasial

Model Regresi Linear Berganda

Persamaan: \[ y\; = \; X\beta \; + \; \epsilon \]

#PERSAMAAN REGRESI
reg.eq = CRIME ~ INC + HOVAL
#OLS
reg.OLS=lm(reg.eq,data=columbus)
summary(reg.OLS)
## 
## Call:
## lm(formula = reg.eq, data = columbus)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.418  -6.388  -1.580   9.052  28.649 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  68.6190     4.7355  14.490  < 2e-16 ***
## INC          -1.5973     0.3341  -4.780 1.83e-05 ***
## HOVAL        -0.2739     0.1032  -2.654   0.0109 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.43 on 46 degrees of freedom
## Multiple R-squared:  0.5524, Adjusted R-squared:  0.5329 
## F-statistic: 28.39 on 2 and 46 DF,  p-value: 9.341e-09
AIC(reg.OLS)
## [1] 382.7545

Uji LM Test

LM Test digunakan untuk deteksi awal efek spasial.

lm.LMtests(reg.OLS,queen.columbus,test=c("LMlag", "LMerr",  "RLMlag", "RLMerr", "SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg.eq, data = columbus)
## test weights: listw
## 
## RSlag = 8.898, df = 1, p-value = 0.002855
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg.eq, data = columbus)
## test weights: listw
## 
## RSerr = 5.2062, df = 1, p-value = 0.02251
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg.eq, data = columbus)
## test weights: listw
## 
## adjRSlag = 3.7357, df = 1, p-value = 0.05326
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg.eq, data = columbus)
## test weights: listw
## 
## adjRSerr = 0.043906, df = 1, p-value = 0.834
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg.eq, data = columbus)
## test weights: listw
## 
## SARMA = 8.9419, df = 2, p-value = 0.01144

Model SCR

Spatial Cross Regressive (SCR) merupakan model regresi spasial, dimana efek spasial melekat pada variabel independen (). Menurut LeSage & Pace (2009) dalam hal pemodelan, variabel lag spasial dari prediktor () dapat berperan langsung dalam menentukan nilai dari variabel dependen (). Persamaan Model SCR: \[y\; = \;X\beta \; + \; WX\gamma \; + \; \epsilon\]

#SLX (SCR)
reg.SLX=lmSLX(reg.eq,data=columbus, queen.columbus)
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)   7.455e+01   6.716e+00   1.110e+01   2.416e-14
## INC          -1.097e+00   3.738e-01  -2.936e+00   5.277e-03
## HOVAL        -2.944e-01   1.017e-01  -2.896e+00   5.868e-03
## lag.INC      -1.399e+00   5.601e-01  -2.497e+00   1.633e-02
## lag.HOVAL     2.148e-01   2.079e-01   1.033e+00   3.071e-01
#Perhitungan AIC
AIC(reg.SLX)
## [1] 379.9412
AIC(reg.OLS)
## [1] 382.7545
#Perbandingan Model dengan OLS
LR.Sarlm(reg.SLX,reg.OLS) # OLS vs SCR
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 6.8133, df = 2, p-value = 0.03315
## sample estimates:
## Log likelihood of reg.SLX Log likelihood of reg.OLS 
##                 -183.9706                 -187.3772

Model SAR

Spatial Autoregressive Model (SAR) disebut juga Spatial Lag Model (SLM) adalah salah satu model spasial dengan pendekatan area dengan memperhitungkan pengaruh spasial lag pada variabel dependen saja. Model ini dinamakan juga Mixed Regressive-Autoregressive karena mengkombinasikan model regresi biasa dengan model regresi spasial lag pada variabel dependen (Anselin, 1988). Persamaan Model SAR:

#SAR
reg.SAR=lagsarlm(reg.eq,data=columbus, queen.columbus)
summary(reg.SAR)
## 
## Call:lagsarlm(formula = reg.eq, data = columbus, listw = queen.columbus)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -37.652017  -5.334611   0.071473   6.107196  23.302617 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.603248   7.257404  6.2837 3.306e-10
## INC         -1.048728   0.307406 -3.4115  0.000646
## HOVAL       -0.266335   0.089096 -2.9893  0.002796
## 
## Rho: 0.42333, LR test value: 9.4065, p-value: 0.0021621
## Asymptotic standard error: 0.11951
##     z-value: 3.5422, p-value: 0.00039686
## Wald statistic: 12.547, p-value: 0.00039686
## 
## Log likelihood: -182.674 for lag model
## ML residual variance (sigma squared): 96.857, (sigma: 9.8416)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 375.35, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 0.24703, p-value: 0.61917
#Perhitungan AIC
AIC(reg.SAR)
## [1] 375.3479
AIC(reg.OLS)
## [1] 382.7545
#Perbandingan Model dengan OLS
LR.Sarlm(reg.SAR,reg.OLS) # OLS vs SAR
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 9.4065, df = 1, p-value = 0.002162
## sample estimates:
## Log likelihood of reg.SAR Log likelihood of reg.OLS 
##                 -182.6740                 -187.3772

Model SEM

Spatial Error Model (SEM) dapat digunakan saat nilai error pada suatu lokasi berkorelasi dengan nilai error dengan lokasi sekitarnya atau bisa dikatakan terdapat korelasi spasial antar error. Pada model SEM, bentuk error pada lokasi i merupakan fungsi dari error pada lokasi j dimana j merupakan suatu lokasi yang terletak disekitar lokasi i. Persamaan model SEM: \[y\; = \;X\beta \; + \; u \; , \quad u\; = \;\lambda Wu \; + \; \epsilon\]

#SEM
reg.SEM=errorsarlm(reg.eq,data=columbus, queen.columbus)
summary(reg.SEM)
## 
## Call:errorsarlm(formula = reg.eq, data = columbus, listw = queen.columbus)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -34.65998  -6.16943  -0.70623   7.75392  23.43878 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 60.279470   5.365594 11.2344 < 2.2e-16
## INC         -0.957305   0.334231 -2.8642 0.0041806
## HOVAL       -0.304559   0.092047 -3.3087 0.0009372
## 
## Lambda: 0.54675, LR test value: 7.2556, p-value: 0.0070679
## Asymptotic standard error: 0.13805
##     z-value: 3.9605, p-value: 7.4786e-05
## Wald statistic: 15.686, p-value: 7.4786e-05
## 
## Log likelihood: -183.7494 for error model
## ML residual variance (sigma squared): 97.674, (sigma: 9.883)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 377.5, (AIC for lm: 382.75)
#Perhitungan AIC
AIC(reg.SEM)
## [1] 377.4989
AIC(reg.OLS)
## [1] 382.7545
#Perbandingan Model dengan OLS
LR.Sarlm(reg.SEM,reg.OLS) # SEM vs OLS
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 7.2556, df = 1, p-value = 0.007068
## sample estimates:
## Log likelihood of reg.SEM Log likelihood of reg.OLS 
##                 -183.7494                 -187.3772

Model Manski

Model ini merupakan model Regresi Spasial paling lengkap. Model ini memuat semua efek spasial pada seluruh komponen model regresi. Persamaan Model Manski: \[y\; = \;\rho Wy \; +\; X\beta \; + \; WX\gamma \; + \; u \; , \quad u \; = \;\lambda Wu \; + \; \epsilon\]

#Manski
reg.Manski=sacsarlm(reg.eq,data=columbus, queen.columbus, type="sacmixed") 
summary(reg.Manski)
## 
## Call:sacsarlm(formula = reg.eq, data = columbus, listw = queen.columbus, 
##     type = "sacmixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -37.64377  -6.54301  -0.36732   5.95810  22.92013 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) 53.04801   71.85207  0.7383 0.460335
## INC         -0.95880    0.45677 -2.0991 0.035809
## HOVAL       -0.28903    0.10277 -2.8126 0.004915
## lag.INC     -0.77341    1.78153 -0.4341 0.664196
## lag.HOVAL    0.21846    0.29863  0.7316 0.464436
## 
## Rho: 0.28438
## Asymptotic standard error: 0.99595
##     z-value: 0.28553, p-value: 0.77524
## Lambda: 0.16325
## Asymptotic standard error: 1.0908
##     z-value: 0.14966, p-value: 0.88103
## 
## LR test value: 11.594, p-value: 0.020638
## 
## Log likelihood: -181.5801 for sacmixed model
## ML residual variance (sigma squared): 94.485, (sigma: 9.7204)
## Number of observations: 49 
## Number of parameters estimated: 8 
## AIC: 379.16, (AIC for lm: 382.75)

Perbandingan Kinerja Model

Pemilihan model terbaik dapat dilakukan dengan kriteria AIC terkecil.

#Perhitungan AIC
AIC(reg.SLX)
## [1] 379.9412
AIC(reg.SAR)
## [1] 375.3479
AIC(reg.SEM)
## [1] 377.4989
AIC(reg.Manski)
## [1] 379.1603
#Perbandingan dengan Model Regresi Spasial Lainnya
LR.Sarlm(reg.Manski,reg.OLS) # vs OLS
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 11.594, df = 4, p-value = 0.02064
## sample estimates:
## Log likelihood of reg.Manski    Log likelihood of reg.OLS 
##                    -181.5801                    -187.3772
#Perbandingan dengan Model Regresi Spasial Lainnya
LR.Sarlm(reg.Manski,reg.SLX) # vs SCR
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 4.7809, df = 2, p-value = 0.09159
## sample estimates:
## Log likelihood of reg.Manski    Log likelihood of reg.SLX 
##                    -181.5801                    -183.9706
#Perbandingan dengan Model Regresi Spasial Lainnya
LR.Sarlm(reg.Manski,reg.SAR) # vs SAR
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 2.1877, df = 3, p-value = 0.5344
## sample estimates:
## Log likelihood of reg.Manski    Log likelihood of reg.SAR 
##                    -181.5801                    -182.6740
#Perbandingan dengan Model Regresi Spasial Lainnya
LR.Sarlm(reg.Manski,reg.SEM) # vs SEM
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 4.3386, df = 3, p-value = 0.2271
## sample estimates:
## Log likelihood of reg.Manski    Log likelihood of reg.SEM 
##                    -181.5801                    -183.7494

Berdasarkan output yg dihasilkan, dapat disimpulkan bahwa model terbaiknya adalah model SAR karena memiliki nilai AIC terkecil.