Spatial Regression Modelling

Eksplorasi Data

Import Data

#load package
#Packages Yang Digunakan 
library(spgwr)
library(readxl) 
library(openxlsx) 
library(spdep) 
library(rgdal) 
library(raster) 
library(corrplot)  
library(DescTools) 
library(nortest) 
library(car) 
library(spatialreg)
library(sf)
library(nortest) 
library(DescTools) 
library(lmtest)
library(tidyverse)


#input file excel
datauas <-read.csv("D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/UAS/data-uas-007.csv") 
head(datauas)
##   X ID_KAB              KAB        y     x1   x2   x3
## 1 0   3101 KEPULAUAN SERIBU 505.8809 267.44 2.18 1.58
## 2 1   3171  JAKARTA SELATAN 567.7238 244.14 1.68 2.69
## 3 2   3172    JAKARTA TIMUR 452.8043 137.81 2.68 1.46
## 4 3   3173    JAKARTA PUSAT 384.7793  38.78 0.55 2.33
## 5 4   3174    JAKARTA BARAT 470.8095 230.17 2.96 2.07
## 6 5   3175    JAKARTA UTARA 438.7692 113.31 0.77 1.82
dim(datauas)
## [1] 119   7

Data terdiri dari 119 baris/amatan dan 7 kolom/peubah. Terdapat satu peubah respon Y dan 3 peubah bebas X1 X2 dan X3. Dengan,

Y = PDRB kabupaten/kota

X1 = jumlah tenaga kerja (ribu orang)

X2 = pendapatan asli daerah (dalam triliun rupiah)

X3 = inflasi (persen)

#Input peta shp
peta<-readOGR(dsn="D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/UAS/Jawamap", layer="jawa") 
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\MATERI KULIAH S2 IPB\SEMESTER 2\SPASIAL\UAS\Jawamap", layer: "jawa"
## with 119 features
## It has 5 fields

Cek Missing Value Data

datauas[!complete.cases(datauas),]
## [1] X      ID_KAB KAB    y      x1     x2     x3    
## <0 rows> (or 0-length row.names)

Tidak terdapat missing value/nilai yang hilang pada data.

Summary Data

summary(datauas)
##        X             ID_KAB         KAB                  y        
##  Min.   :  0.0   Min.   :3101   Length:119         Min.   :361.1  
##  1st Qu.: 29.5   1st Qu.:3276   Class :character   1st Qu.:432.2  
##  Median : 59.0   Median :3327   Mode  :character   Median :490.5  
##  Mean   : 59.0   Mean   :3386                      Mean   :484.1  
##  3rd Qu.: 88.5   3rd Qu.:3516                      3rd Qu.:534.0  
##  Max.   :118.0   Max.   :3674                      Max.   :604.3  
##        x1               x2              x3       
##  Min.   : 19.13   Min.   :0.120   Min.   :1.170  
##  1st Qu.: 79.02   1st Qu.:1.660   1st Qu.:1.650  
##  Median :180.03   Median :3.310   Median :1.970  
##  Mean   :158.44   Mean   :3.500   Mean   :1.987  
##  3rd Qu.:240.86   3rd Qu.:5.305   3rd Qu.:2.330  
##  Max.   :297.68   Max.   :7.030   Max.   :2.740

Sebaran Data

Sebaran dari data akan dieksplorasi menggunakan boxplot dan diperoleh hasil sebagai berikut:

boxplot(datauas$y,main="Sebaran PDRB kabupaten/kota") 

boxplot(datauas$x1,main="Sebaran jumlah tenaga kerja") 

boxplot(datauas$x2,main="Sebaran pendapatan asli daerah") 

boxplot(datauas$x3,main="Sebaran inflasi")

Berdasrakan Boxplot di atas, terlihat bahwa tidak terdapat nilai pencilan pada peubah Y, X1, X2, dan X3.

Relationship of the data

  1. Scatter Plot X1 vs Y
plot(datauas$x1, datauas$y,
  xlab="X1", 
  ylab="Y",
  main="Scatter Plot of X1 and Y",
  pch=20, col="orange", cex=2)
reg.klasik = lm(y~x1, data=datauas)
lines.lm(reg.klasik, col=2, add=T)

Scatter plot X1 vs Y di atas menunjukkan bahwa terdapat hubungan linier positif antara X1 (jumlah tenaga kerja) dengan Y (PDRB).

  1. Scatter Plot X2 vs Y
plot(datauas$x2, datauas$y,
  xlab="X2", 
  ylab="Y",
  main="Scatter Plot of X2 and Y",
  pch=20, col="orange", cex=2)
reg.klasik = lm(y~x2, data=datauas)
lines.lm(reg.klasik, col=2, add=T)

Scatter plot X2 vs Y di atas menunjukkan bahwa terdapat hubungan linier positif antara X2 (pendapatan asli daerah) dengan Y(PDRB).

  1. Scatter Plot X3 vs Y
plot(datauas$x3, datauas$y,
  xlab="X3", 
  ylab="Y",
  main="Scatter Plot of X3 and Y",
  pch=20, col="orange", cex=2)
reg.klasik = lm(y~x3, data=datauas)
lines.lm(reg.klasik, col=2, add=T)

Scatter plot X3 vs Y di atas menunjukkan bahwa terdapat hubungan linier negatif antara X3 (inflasi) dengan Y (PDRB).

Correlation of the data

1. Matriks Korelasi

corrplot(cor(datauas[,4:7]), method = "number")

Berdasarkan matriks korelasi di atas terlihat bahwa Y memiliki korelasi positif dengan X1 dan X2, dan memiliki korelasi negatif dengan X3.

2. Pengujian Korelasi

  1. X1 (jumlah tenaga kerja) vs Y (PDRB)

H0: Tidak terdapat korelasi antara peubah X1 dan Y

H1: Terdapat korelasi antara peubah X1 dan Y

cor.test(datauas$y, datauas$x1)
## 
##  Pearson's product-moment correlation
## 
## data:  datauas$y and datauas$x1
## t = 23.692, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8725553 0.9363579
## sample estimates:
##     cor 
## 0.90968

Kesimpulan :

Karena diperoleh p-value < 0.05 maka dapat disimpulkan bahwa H0 ditolak yang berarti terdapat korelasi antara peubah X1 dan Y.

  1. X2 (pendapatan asli daerah) vs Y (PDRB)

H0: Tidak terdapat korelasi antara peubah X2 dan Y

H1: Terdapat korelasi antara peubah X2 dan Y

cor.test(datauas$y, datauas$x2)
## 
##  Pearson's product-moment correlation
## 
## data:  datauas$y and datauas$x2
## t = 2.0517, df = 117, p-value = 0.04243
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.006579527 0.354460279
## sample estimates:
##       cor 
## 0.1863544

Kesimpulan :

Karena diperoleh p-value < 0.05 maka dapat disimpulkan bahwa H0 ditolak yang berarti terdapat korelasi antara peubah X2 dan Y.

  1. X3 (inflasi) vs Y (PDRB)

H0: Tidak terdapat korelasi antara peubah X3 dan Y

H1: Terdapat korelasi antara peubah X3 dan Y

cor.test(datauas$y, datauas$x3)
## 
##  Pearson's product-moment correlation
## 
## data:  datauas$y and datauas$x3
## t = -1.6094, df = 117, p-value = 0.1102
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.31872560  0.03371542
## sample estimates:
##        cor 
## -0.1471733

Kesimpulan :

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa H0 gagal ditolak yang berarti tidak terdapat korelasi antara peubah X3 dan Y.

Sebaran Spasial Data

1. Sebaran Spasial Y

k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)

peta$y<- datauas$y
spplot(peta, "y", col.regions=color, main="Sebaran Spasial PDRB")

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol pada peubah Y. Hal ini tampak dari gradasi warna yang cenderung mengumpul, seperti pada warna hijau, kuning, dan oranye.

2. Sebaran Spasial X1

k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)

peta$x1<- datauas$x1
spplot(peta, "x1", col.regions=color, main="Sebaran Spasial jumlah tenaga kerja")

3. Sebaran Spasial X2

k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)

peta$x2<- datauas$x2
spplot(peta, "x2", col.regions=color, main="Sebaran Spasial pendapatan asli daerah")

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol pada pendapatan asli daerah. Hal ini tampak dari gradasi warna yang cenderung mengumpul, seperti pada warna hijau, kuning, dan merah

3. Sebaran Spasial X3

k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)

peta$x3<- datauas$x3
spplot(peta, "x3", col.regions=color, main="Sebaran Spasial Inflasi")

Matriks Pembobot Spasial

Matriks bobot berdasarkan jarak

longlat <- coordinates(peta) 
head(longlat)
##       [,1]      [,2]
## 0 106.4933 -5.796823
## 1 106.8100 -6.272549
## 2 106.9003 -6.255373
## 3 106.8351 -6.181230
## 4 106.7483 -6.165469
## 5 106.8695 -6.129185
plot(longlat)

peta$long <- longlat[,1]

peta$lat <- longlat[,2] 
coords <- peta[c("long","lat")] 
class(coords)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
# menghitung matriks jarak 
djarak<-dist(longlat) 
m.djarak<-as.matrix(djarak) #matriks jarak

1. K-Nearest Neighbor Weight (Bobot k-Tetangga Terdekat)

#k=3, 3 tetangga terdekat 
W.knn<-knn2nb(knearneigh(longlat,k=3,longlat=TRUE)) 
W.knn
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 357 
## Percentage nonzero weights: 2.521008 
## Average number of links: 3 
## Non-symmetric neighbours list
class(W.knn)
## [1] "nb"
W.knn1 <- nb2listw(W.knn,style='W') 
W.knn1
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 357 
## Percentage nonzero weights: 2.521008 
## Average number of links: 3 
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 68.77778 495.1111
plot(peta, col='gray', border='blue', main ="KNN dengan k=3") 

plot(W.knn1, longlat, col='red', lwd=2, add=TRUE)

2. Radial Distance Weight

#d=70 ditetapkan nilai ambang dmax adalah 70 km 
W.dmax<-dnearneigh(longlat,0,70,longlat=TRUE) 

class(W.dmax) #nb
## [1] "nb"
W.dmax.s <- nb2listw(W.dmax,style='W') 

W.dmax.s
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 1300 
## Percentage nonzero weights: 9.180143 
## Average number of links: 10.92437 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 25.99101 482.6102

3. Power Distance Weight / Invers Distance Weight (Bobot Jarak Invers)

#Alpha = 1 
alpha1=1 
W.idw <-1/(m.djarak^alpha1) 
class(W.idw)
## [1] "matrix" "array"
#dinormalisasi 
diag(W.idw)<-0 
rtot<-rowSums(W.idw,na.rm=TRUE)

W.idw.sd<-W.idw/rtot #row-normalized 
rowSums(W.idw.sd,na.rm=TRUE)
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
class(W.idw.sd)
## [1] "matrix" "array"
W.idw.1 = mat2listw(W.idw.sd,style='W') 

summary(W.idw.1)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 14042 
## Percentage nonzero weights: 99.15966 
## Average number of links: 118 
## Link number distribution:
## 
## 118 
## 119 
## 119 least connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 most connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 8.084636 479.8627
#Alpha = 2
alpha1=2 
W.idw2 <-2/(m.djarak^alpha1) 
class(W.idw2)
## [1] "matrix" "array"
#dinormalisasi 
diag(W.idw2)<-0 
rtot<-rowSums(W.idw2,na.rm=TRUE)

W.idw.sd2<-W.idw2/rtot #row-normalized 
rowSums(W.idw.sd2,na.rm=TRUE)
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
class(W.idw.sd2)
## [1] "matrix" "array"
W.idw.22 = mat2listw(W.idw.sd2,style='W') 

summary(W.idw.22)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 14042 
## Percentage nonzero weights: 99.15966 
## Average number of links: 118 
## Link number distribution:
## 
## 118 
## 119 
## 119 least connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 most connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 37.02077 490.9232

4. Exponential Distance Weight

alpha=1 

W.e<-exp((-alpha)*m.djarak) #dinormalisasi 
diag(W.e)<-0 

rtot<-rowSums(W.e,na.rm=TRUE) 

W.e.sd<-W.e/rtot #row-normalized 
rowSums(W.e.sd,na.rm=TRUE)
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
class(W.e.sd)
## [1] "matrix" "array"
W.ed1 = mat2listw(W.e.sd,style='W') 

summary(W.ed1)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 14042 
## Percentage nonzero weights: 99.15966 
## Average number of links: 118 
## Link number distribution:
## 
## 118 
## 119 
## 119 least connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 most connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 4.839428 478.0838
alpha=2

W.e2<-exp((-alpha)*m.djarak) #dinormalisasi 
diag(W.e2)<-0 

rtot<-rowSums(W.e2,na.rm=TRUE) 

W.e.sd2<-W.e2/rtot #row-normalized 
rowSums(W.e.sd2,na.rm=TRUE)
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
class(W.e.sd2)
## [1] "matrix" "array"
W.ed2 = mat2listw(W.e.sd2,style='W') 

summary(W.ed2)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 14042 
## Percentage nonzero weights: 99.15966 
## Average number of links: 118 
## Link number distribution:
## 
## 118 
## 119 
## 119 least connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 most connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 119 14161 119 8.900784 479.4497

List bobot berdasarkan jarak:

  1. W.knn1 (KNN)

  2. W.dmax.s (Radial Distance Weight)

  3. W.idw.1 (Power Distance Weight)

  4. W.idw.22 (Power Distance Weight)

  5. W.ed1 (Exponential Distance Weight)

  6. W.ed2 (Exponential Distance Weight)

Matriks Bobot berdasarkan ketetanggaan (Spatial Contiguity Weight)

1. Rook Contiguity

library(readr)
peta1 <- read_sf("D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/UAS/Jawamap/jawa.shp")
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
#Rook 
W.rook<-poly2nb(peta1,queen=FALSE) 
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
W.rook #matriks bobot 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
class(W.rook) #nb
## [1] "nb"
W.rook.s <- nb2listw(W.rook,style='W', zero.policy = TRUE) #W is row standardised (sums over all links to n). Standardisasi Baris 

print(W.rook.s, zero.policy=TRUE)
## Characteristics of weights list object:
## 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
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 118 13924 118 66.41255 515.0633

2. Queen Contiguity

# queen 
W.queen<-poly2nb(peta,queen=TRUE) 
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
W.queen #matriks bobot Rook
## 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:
## 0
class(W.queen) #nb
## [1] "nb"
W.queen.s <- nb2listw(W.queen,style='W', zero.policy = TRUE) #W is row standardised (sums over all links to n). Standardisasi Baris 

print(W.queen.s, zero.policy=TRUE)
## Characteristics of weights list object:
## 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:
## 0
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 118 13924 118 66.14482 515.6022

Pemilihan Matriks Bobot

1. KNN

MI.knn <- moran(datauas$y, W.knn1, n=length(W.knn1$neighbours), S0=Szero(W.knn1))

2. Radial Distance Weight

MI.radial <- moran(datauas$y, W.dmax.s, n=length(W.dmax.s$neighbours), S0=Szero(W.dmax.s))

3. Power Distance Weight

MI.power1 <- moran(datauas$y, W.idw.1, n=length(W.idw.1$neighbours), S0=Szero(W.idw.1))

MI.power2 <- moran(datauas$y, W.idw.22, n=length(W.idw.22$neighbours), S0=Szero(W.idw.22))

4. Exponential Distance Weight

MI.exp1 <- moran(datauas$y, W.ed1, n=length(W.ed1$neighbours), S0=Szero(W.ed1))

MI.exp2 <- moran(datauas$y, W.ed2, n=length(W.ed2$neighbours), S0=Szero(W.ed2))

5. Spatial Contiguity

MI.rook <- moran.test(datauas$y, W.rook.s, randomisation = TRUE, zero.policy = TRUE)
MI.rook$estimate
## Moran I statistic       Expectation          Variance 
##      -0.009409745      -0.008547009       0.004640627
MI.queen <- moran.test(datauas$y, W.queen.s, randomisation = TRUE, zero.policy = TRUE)
MI.queen$estimate
## Moran I statistic       Expectation          Variance 
##      -0.009671002      -0.008547009       0.004620888

Perbandingan Nilai Indeks Moran’s

moranindeks<-data.frame(
  "Matriks Bobot"=c("KNN (k=5)", "Radial distance weight", "Power distance weight (alpha=1)", "Power distance weight (alpha=2)", "Exponential distance weight (alpha=1)", "Exponential distance weight (alpha=2)", "Rook Contiguity", "Queen Contiguity"),
  "Nilai Indeks Moran"=c(MI.knn$I, MI.radial$I, MI.power1$I, MI.power2$I, MI.exp1$I, MI.exp2$I, -0.009409745, -0.009671002))

moranindeks
##                           Matriks.Bobot Nilai.Indeks.Moran
## 1                             KNN (k=5)        0.061699439
## 2                Radial distance weight       -0.059527412
## 3       Power distance weight (alpha=1)        0.049911441
## 4       Power distance weight (alpha=2)        0.102902435
## 5 Exponential distance weight (alpha=1)       -0.013118217
## 6 Exponential distance weight (alpha=2)       -0.008933659
## 7                       Rook Contiguity       -0.009409745
## 8                      Queen Contiguity       -0.009671002

Kesimpulan:

Matriks bobot yang menghasilkan nilai Indeks Moran terbesar adalah matriks bobot Power distance weight dengan alpha=2.

Selanjutnya akan dilakukan pengujian indeks moran dengan menggunakan matriks bobot Power distance weight dengan alpha=2:

H0: I=0 (tidak ada autokorelasi spasial antarlokasi) H1: I>0 (ada autokorelasi spasial antar lokasi)

Woptimum <- W.idw.22
moran.test(datauas$y, Woptimum, randomisation = TRUE, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  datauas$y  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = 2.2341, p-value = 0.01274
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.102902435      -0.008474576       0.002485374

Kesimpulan :

Berdasarkan hasil di atas, diperoleh p-value < 0.05 berarti H0 ditolak, dengan alternative hypothesis: greater berarti pada taraf nyata 5%, terdapat autokorelasi spasial positif pada peubah Y (PDRB).

Moran Plot

moran.plot(datauas$y, Woptimum, labels=datauas$KABKOT)

Plot di atas menunjukkan slope positif yang berarti nilai autokorelasi spasial pada Y (PDRB) bertanda positif.

Pengujian Efek Dependensi Spasial

Regresi Linier Y vs X1,X2, dan X3

reg.klasik1 = lm(y~x1+x2+x3, data = datauas)

summary(reg.klasik1)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datauas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.782 -14.051  -0.223  17.259  51.436 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 371.14937   12.25102  30.295  < 2e-16 ***
## x1            0.62775    0.02532  24.792  < 2e-16 ***
## x2            4.52374    1.06201   4.260 4.21e-05 ***
## x3           -1.18565    5.11064  -0.232    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.92 on 115 degrees of freedom
## Multiple R-squared:  0.8511, Adjusted R-squared:  0.8472 
## F-statistic: 219.1 on 3 and 115 DF,  p-value: < 2.2e-16

Pada tahap eksplorasi data saat dilakukan pengujian korelasi antara peubah Y dengan X3, diperoleh hasil bahwa tidak terdapat korelasi antara Y dengan X3. Berikutnya saat dilakukan pemodelan menggunakan regresi linier berganda koefisien X3 tidak berpengaruh signifikan/nyata terhadap model, sehingga X3 akan dikeluarkan dari model, dan diperoleh hasil regresi linier berganda sebagai berikut.

Regresi Linier Y vs X1 dan X2

reg.klasik = lm(y~x1+x2, data = datauas)

summary(reg.klasik)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = datauas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.122 -14.112  -0.248  17.382  51.546 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 368.64208    5.74597  64.157  < 2e-16 ***
## x1            0.62863    0.02493  25.214  < 2e-16 ***
## x2            4.52701    1.05757   4.281 3.86e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.83 on 116 degrees of freedom
## Multiple R-squared:  0.851,  Adjusted R-squared:  0.8485 
## F-statistic: 331.4 on 2 and 116 DF,  p-value: < 2.2e-16

Penduga parameter dalam regresi linier Y vs X1 dan X2 telah signifikan/nyata.

Uji Multikolinearitas

car::vif(reg.klasik)
##       x1       x2 
## 1.001323 1.001323

Kesimpulan:

Karena nilai VIF pada x1 dan x2 kurang dari 5, maka dapat disimpulkan bahwa tidak terdapat multikolinearitas antarpeubah penjelas.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.regklasik<-reg.klasik$residuals
ad.test(err.regklasik) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.regklasik
## A = 0.37095, p-value = 0.418
hist(err.regklasik)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest(reg.klasik) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  reg.klasik
## BP = 0.14401, df = 2, p-value = 0.9305

Kesimpulan:

Karena diperoleh p-value>0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(reg.klasik$residuals, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  reg.klasik$residuals  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = 3.8174, p-value = 0.0001349
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.180859997      -0.008474576       0.002459921

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value<0.05 yang berarti tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model regresi klasik pada taraf nyata 5%.

Karena terdapat autokorelasi spasial pada sisaan model regresi klasik maka berikutnya akan dilakukan pengujian menggunakan uji LM (lagrange multiplier) untuk mengidentifikasi model dependensi spasial yang dapat digunakan pada kasus ini.

Uji Lagrange Multiplier

ujilm <- lm.LMtests(reg.klasik,
                    listw=Woptimum,
                    zero.policy = TRUE,
                    test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
summary(ujilm)
##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = y ~ x1 + x2, data = datauas)
## weights: Woptimum
##  
##        statistic parameter   p.value    
## LMerr   12.51220         1 0.0004043 ***
## RLMerr   0.50671         1 0.4765672    
## LMlag   27.08995         1 1.942e-07 ***
## RLMlag  15.08446         1 0.0001028 ***
## SARMA   27.59666         2 1.017e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan hasil di atas, diperoleh hasil uji LM SEM, SAR, dan GSM signifikan pada taraf 5%. Sehingga, model yang akan dicoba adalah SEM, SAR, dan GSM.

Spatial Regression Modelling

1. Spatial Error Model (SEM)

SEM <- errorsarlm(y ~ x1 + x2,data=datauas,Woptimum, zero.policy = TRUE)
summary(SEM)
## 
## Call:errorsarlm(formula = y ~ x1 + x2, data = datauas, listw = Woptimum, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -56.9864 -14.1514  -1.4434  14.0970  56.1151 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 371.844477   6.102909 60.9291 < 2.2e-16
## x1            0.605783   0.022831 26.5337 < 2.2e-16
## x2            4.766685   0.960665  4.9619 6.982e-07
## 
## Lambda: 0.4323, LR test value: 10.15, p-value: 0.0014431
## Asymptotic standard error: 0.12997
##     z-value: 3.3261, p-value: 0.00088071
## Wald statistic: 11.063, p-value: 0.00088071
## 
## Log likelihood: -539.5869 for error model
## ML residual variance (sigma squared): 492.41, (sigma: 22.19)
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 1089.2, (AIC for lm: 1097.3)

Output di atas menunjukkan bahwa koefisien Lambda signifikan pada taraf nyata 5%. AIC model SEM adalah sebesar 1089.2.

Selanjutnya dilakukan pemeriksaan sisaan model SEM.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.sem<-SEM$residuals
ad.test(err.sem) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.39387, p-value = 0.3696
hist(err.sem)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest.Sarlm(SEM) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.4826, df = 2, p-value = 0.4765

Kesimpulan:

Karena diperoleh p-value >0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(SEM$residuals, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  SEM$residuals  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = -0.17489, p-value = 0.8612
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.017151992      -0.008474576       0.002461816

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SEM pada taraf nyata 5%.

2. Spatial Autoregressive Model (SAR)

sar<-lagsarlm(y ~ x1 + x2,data=datauas,Woptimum)
summary(sar)
## 
## Call:lagsarlm(formula = y ~ x1 + x2, data = datauas, listw = Woptimum)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -54.606088 -15.372896  -0.015188  12.157025  49.907866 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 214.030404  35.410499  6.0443 1.501e-09
## x1            0.620784   0.022517 27.5693 < 2.2e-16
## x2            4.459561   0.953652  4.6763 2.921e-06
## 
## Rho: 0.32395, LR test value: 19.711, p-value: 9.0061e-06
## Asymptotic standard error: 0.073413
##     z-value: 4.4127, p-value: 1.021e-05
## Wald statistic: 19.472, p-value: 1.021e-05
## 
## Log likelihood: -534.8061 for lag model
## ML residual variance (sigma squared): 461.02, (sigma: 21.471)
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 1079.6, (AIC for lm: 1097.3)
## LM test for residual autocorrelation
## test value: 0.98273, p-value: 0.32152

Output di atas menunjukkan bahwa koefisien Rho signifikan pada taraf nyata 5%. AIC model SAR adalah sebesar 1079.6.

Selanjutnya dilakukan pemeriksaan sisaan model SAR.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.sar<-sar$residuals
ad.test(err.sar) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.33555, p-value = 0.504
hist(err.sar)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest.Sarlm(sar) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.2736, df = 2, p-value = 0.529

Kesimpulan:

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(sar$residuals, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  sar$residuals  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = -0.71008, p-value = 0.4777
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.043747009      -0.008474576       0.002467472

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga dapat disimpulkan bahwa tidak terdapat autokorelasi pada sisaan model SAR pada taraf nyata 5%.

3. GSM/SARMA

SARMA <- sacsarlm(y ~ x1 + x2,data=datauas,Woptimum, zero.policy = TRUE)
summary(SARMA)
## 
## Call:sacsarlm(formula = y ~ x1 + x2, data = datauas, listw = Woptimum, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -52.74972 -15.49177  -0.75399  12.68545  49.18016 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 186.666013  32.753730  5.6991 1.205e-08
## x1            0.617955   0.022712 27.2079 < 2.2e-16
## x2            4.365375   0.939652  4.6457 3.389e-06
## 
## Rho: 0.3824
## Asymptotic standard error: 0.069011
##     z-value: 5.5412, p-value: 3.0049e-08
## Lambda: -0.25227
## Asymptotic standard error: 0.18169
##     z-value: -1.3885, p-value: 0.16499
## 
## LR test value: 21.086, p-value: 2.6381e-05
## 
## Log likelihood: -534.1189 for sac model
## ML residual variance (sigma squared): 448.18, (sigma: 21.17)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 1080.2, (AIC for lm: 1097.3)

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

Selanjutnya dilakukan pemeriksaan sisaan model SARMA.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.sarma<-SARMA$residuals
ad.test(err.sarma) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sarma
## A = 0.35112, p-value = 0.4645
hist(err.sarma)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest.Sarlm(SARMA) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.79348, df = 2, p-value = 0.6725

Kesimpulan:

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(err.sarma, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  err.sarma  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = 0.082221, p-value = 0.9345
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.004389037      -0.008474576       0.002469098

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SARMA pada taraf nyata 5%.

Berikutnya akan dicoba juga model dependensi spasial yang lain diantaranya yaitu SDM (Spatial Durbin Model), Spatial cross-regressive (SLX), dan Spatial Durbin Error Model (SDEM) yang kemudian dapat dibandingkan kebaikan modelnya dengan model SAR, SEM, dan SARMA.

SDM : Model regresi dengan peubah respon dan peubah bebas yang mempunyai dependensi spasial.

SLX : Model regresi dengan peubah bebas yang mempunyai dependensi spasial.

SDEM : Model regresi dengan peubah bebas dan galat yang mempunyai dependensi spasial.

1. SDM

sdm <- lagsarlm(y ~ x1 + x2, data=datauas, Woptimum, zero.policy = TRUE, type="mixed");
summary(sdm)
## 
## Call:lagsarlm(formula = y ~ x1 + x2, data = datauas, listw = Woptimum, 
##     type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -53.6224 -15.1502  -0.6623  12.0999  50.6339 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 232.130571  50.555955  4.5916 4.400e-06
## x1            0.622831   0.022918 27.1770 < 2.2e-16
## x2            4.402166   0.958490  4.5928 4.373e-06
## lag.x1        0.054829   0.124816  0.4393    0.6605
## lag.x2       -0.479943   2.504749 -0.1916    0.8480
## 
## Rho: 0.27193, LR test value: 3.7431, p-value: 0.053026
## Asymptotic standard error: 0.14095
##     z-value: 1.9292, p-value: 0.0537
## Wald statistic: 3.722, p-value: 0.0537
## 
## Log likelihood: -534.6382 for mixed model
## ML residual variance (sigma squared): 462.12, (sigma: 21.497)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1083.3, (AIC for lm: 1085)
## LM test for residual autocorrelation
## test value: 3.0185, p-value: 0.08232

Output di atas memperlihatkan bahwa koefisien lag.x1, lag.x2, dan rho tidak signifikan. AIC model SDM adalah sebesar 1083.3.

Selanjutnya dilakukan pemeriksaan sisaan model SDM.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.sdm<-sdm$residuals
ad.test(err.sdm) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sdm
## A = 0.33932, p-value = 0.4943
hist(err.sdm)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest.Sarlm(sdm) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.4029, df = 4, p-value = 0.8437

Kesimpulan:

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(err.sdm, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  err.sdm  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = -0.4292, p-value = 0.6678
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.029796999      -0.008474576       0.002468062

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SDM pada taraf nyata 5%.

2. SLX

slx <- lmSLX(y ~ x1 + x2, data=datauas, Woptimum, zero.policy = TRUE, Durbin=TRUE);
summary(slx)
## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -51.82 -15.73  -0.45  13.79  51.86 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 322.99522   16.23355  19.897  < 2e-16 ***
## x1            0.63141    0.02351  26.863  < 2e-16 ***
## x2            4.25604    0.99913   4.260 4.23e-05 ***
## lag.x1        0.27427    0.06800   4.033 9.98e-05 ***
## lag.x2        0.99483    2.48105   0.401    0.689    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.44 on 114 degrees of freedom
## Multiple R-squared:  0.8701, Adjusted R-squared:  0.8656 
## F-statistic: 190.9 on 4 and 114 DF,  p-value: < 2.2e-16

Output di atas memperlihatkan bahwa koefisien lag.x1 signifikan sedangkan koefisien lag.x2 tidak signifikan.

Selanjutnya dilakukan pemeriksaan sisaan model SLX.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.slx<-slx$residuals
ad.test(err.slx) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.slx
## A = 0.23261, p-value = 0.7949
hist(err.slx)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest(slx) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  slx
## BP = 1.3262, df = 4, p-value = 0.8569

Kesimpulan:

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(err.slx, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  err.slx  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = 1.5975, p-value = 0.1102
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.070941076      -0.008474576       0.002471316

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SLX pada taraf nyata 5%.

SDEM

sdem <- errorsarlm(y ~ x1 + x2, data=datauas, Woptimum, zero.policy = TRUE, etype="mixed");
summary(sdem)
## 
## Call:errorsarlm(formula = y ~ x1 + x2, data = datauas, listw = Woptimum, 
##     etype = "mixed", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -53.12119 -15.10542  -0.43288  12.33460  51.64130 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 327.916225  17.623874 18.6064 < 2.2e-16
## x1            0.629917   0.023025 27.3583 < 2.2e-16
## x2            4.333794   0.983342  4.4072 1.047e-05
## lag.x1        0.240324   0.072895  3.2969 0.0009778
## lag.x2        1.099304   2.456321  0.4475 0.6544847
## 
## Lambda: 0.22357, LR test value: 1.9239, p-value: 0.16543
## Asymptotic standard error: 0.15175
##     z-value: 1.4733, p-value: 0.14066
## Wald statistic: 2.1707, p-value: 0.14066
## 
## Log likelihood: -535.5478 for error model
## ML residual variance (sigma squared): 471.09, (sigma: 21.705)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1085.1, (AIC for lm: 1085)

Output di atas memperlihatkan bahwa koefisien lag.x1 signifikan sedangkan koefisien lag.x2 dan Lambda tidak signifikan.

Selanjutnya dilakukan pemeriksaan sisaan model SDEM.

Uji Normalitas

H0: sisaan model menyebar normal

H1: sisaan model tidak menyebar normal

err.sdem<-sdem$residuals
ad.test(err.sdem) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sdem
## A = 0.27735, p-value = 0.6473
hist(err.sdem)

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

Kesimpulan:

Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).

Uji Kehomogenan Ragam

H0: ragam sisaan homogen

H1: ragam sisaan tidak homogen

bptest.Sarlm(sdem) #inputnya adalah modelnya
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.2462, df = 4, p-value = 0.8704

Kesimpulan:

Karena diperoleh p-value > 0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.

Uji Autokorelasi Spasial

H0 : Tidak terdapat Autokorelasi pada sisaan

H1 : Terdapat Autokorelasi pada sisaan

moran.test(err.sdem, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  err.sdem  
## weights: Woptimum    
## 
## Moran I statistic standard deviate = 0.10012, p-value = 0.9203
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.003499385      -0.008474576       0.002469534

Kesimpulan:

Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SDEM pada taraf nyata 5%.

Pemilihan Model Dependensi Spasial Terbaik

#Pemilihan Model Berdasarkan Nilai AIC Terkecil
data.frame("MODEL" = c("OLS","SEM","SAR","SARMA", "SDM", "SLX", "SDEM"),
           "AIC" = c(AIC(reg.klasik),AIC(SEM), AIC(sar),AIC(SARMA), AIC(sdm), AIC(slx), AIC(sdem)),
           "Uji Diagnostik"=c("Terdapat autokorelasi spasial", "Seluruh asumsi terpenuhi", "Seluruh asumsi terpenuhi", "Seluruh asumsi terpenuhi", "Seluruh asumsi terpenuhi", "Seluruh asumsi terpenuh", "Seluruh asumsi terpenuhi"))%>% arrange(AIC)
##   MODEL      AIC                Uji.Diagnostik
## 1   SAR 1079.612      Seluruh asumsi terpenuhi
## 2 SARMA 1080.238      Seluruh asumsi terpenuhi
## 3   SDM 1083.276      Seluruh asumsi terpenuhi
## 4   SLX 1085.020       Seluruh asumsi terpenuh
## 5  SDEM 1085.096      Seluruh asumsi terpenuhi
## 6   SEM 1089.174      Seluruh asumsi terpenuhi
## 7   OLS 1097.324 Terdapat autokorelasi spasial

Model terbaik untuk data adalah Model SAR.

Interpretasi Efek Marginal (Spillover)

Model SAR merupakan salah satu model dependensi spasial yang memiliki efek marginal (spillover). Efek ini dapat dibedakan menjadi tiga, yaitu efek langsung (direct effect), efek tidak langsung (indirect effect), dan efek total (total effect): gabungan dari efek langsung dan efek tidak langsung.

Direct Effect : Pengaruh langsung pada peubah respon yang disebabkan oleh adanya perubahan prediktor di lokasi yang sama.

Indirect Effect : Pengaruh pada peubah respon yang disebabkan oleh perubahan prediktor di lokasi yang berbeda.

impacts(sar, listw = Woptimum) #Impact dari model SAR
## Impact measures (lag, exact):
##       Direct  Indirect     Total
## x1 0.6318858 0.2863646 0.9182504
## x2 4.5393117 2.0571725 6.5964842

Terlihat bahwa pengaruh langsung dari peubah x1 adalah sebesar 0.632 artinya saat jumlah tenaga kerja (x1) di wilayah-i meningkat 10 kali, maka PDRB (y) di wilayah tersebut akan bertambah sebesar 6.32 kali per KAB/KOTA, jika nilai x2 dan x3 tetap.

Sedangkan efek tak langsung dari peubah x1 bernilai 0.286. Artinya, jika nilai pendapatan asli daerah (x2) di wilayah-i meningkat 10 kali, maka nilai PDRB (y) di wilayah-j akan bertambah sebesar 2.86 per KAB/KOTA, jika nilai x2 dan x3 tetap.

Interpretasi serupa juga dapat dilakukan terhadap peubah x2 dan x3.