Tugas Analisis Spasial

Berikut ini adalah data dan peta Jawa dengan peubah-peubah sebagai berikut: Y = PDRB kabupaten/kota X1= jumlah tenaga kerja (ribu orang) X2= pendapatan asli daerah (dalam triliun rupiah) X3= inflasi (persen) Akan dilakukan pemodelan spasial.

Eksplorasi Data

Eksplorasi Data yang dilakukan meliputi Distribusi Data Pengamatan, Scatter Plot antara Y dan X, sebaran data Y secara Spasial, jika peubah prediktor lebih dari satu, maka di cek hubungan antar peubah prediktor, penentuan matrix bobot yang digunakan dengan mempertimbangkan besaran nilai Autokorelasi). Berikut adalah eksplorasi yang dilakukan:

Import Data

library(raster)
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.3
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
library(sp)
library(readxl)
library(openxlsx)
library(rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.8.2, released 2023/16/12
## Path to GDAL shared files: C:/Users/DELL/AppData/Local/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.3.1, December 1st, 2023, [PJ_VERSION: 931]
## Path to PROJ shared files: C:/Users/DELL/AppData/Local/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.1-3
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(corrplot)
## corrplot 0.92 loaded
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
library(nortest)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.3.3
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.3.3
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

Data CSV

data <- read.csv("data-tugas-praktikum.csv")
head(data)
##   X ID_KAB              KAB        y     x1   x2   x3
## 1 0   3101 KEPULAUAN SERIBU 155.4070 177.93 7.19 1.18
## 2 1   3171  JAKARTA SELATAN 131.2311 172.59 0.13 1.80
## 3 2   3172    JAKARTA TIMUR 221.7079 267.46 6.30 1.79
## 4 3   3173    JAKARTA PUSAT 109.6782 114.52 0.71 1.88
## 5 4   3174    JAKARTA BARAT 147.1425 129.90 5.24 2.41
## 6 5   3175    JAKARTA UTARA 156.0865 207.88 4.89 1.97

SHP

peta <- readOGR(dsn ="Jawamap", layer = "jawa")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\DELL\Downloads\Jawamap", layer: "jawa"
## with 119 features
## It has 5 fields

Distribusi Data Pengamatan

Cek Missing Value

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

Tidak ada missing data pada data ini.

Summary Data

Berikut summary dari data

summary(data[,-c(1,2)])
##      KAB                  y               x1               x2       
##  Length:119         Min.   : 13.8   Min.   : 16.60   Min.   :0.130  
##  Class :character   1st Qu.:110.2   1st Qu.: 89.35   1st Qu.:2.095  
##  Mode  :character   Median :150.4   Median :158.19   Median :3.930  
##                     Mean   :149.4   Mean   :156.82   Mean   :3.791  
##                     3rd Qu.:192.3   3rd Qu.:226.07   3rd Qu.:5.555  
##                     Max.   :275.3   Max.   :295.50   Max.   :7.190  
##        x3       
##  Min.   :1.180  
##  1st Qu.:1.625  
##  Median :1.960  
##  Mean   :1.977  
##  3rd Qu.:2.345  
##  Max.   :2.790

Boxplot

par(mfrow=c(1,4))
boxplot(data$y, main = "Boxplot y", col = "orange")
boxplot(data$x1, main = "Boxplot x1", col = "orange")
boxplot(data$x2, main = "Boxplot x2", col = "orange")
boxplot(data$x3, main = "Boxplot x3", col = "orange")

Dari Boxplot keempat peubah, tidak ada data yang merupakan pencilan

Scatter Plot

Scatter plot ini untuk melihat pola hubungan y dengan masing-masing peubah x1, x2, dan x3, serta pola hubungan antar peubah x

library(DescTools)
pairs(data[,4:7], pch = 19, col = "purple", lower.panel = NULL)

Dari scatter plot di atas, terlihat bahwa adanya hubungan linear positif antara peubah respon y dan x1 maupun x2, dimana x1 tampaknya hubungannya lebih kuat. Sedangkan peubah y dengan x3 tampaknya memiliki hubungan linear negatif yang kecil. Serta antar peubah x tampaknya tidak memiliki hubungan karena tidak memiliki pola tertentu.

Hubungan antar Peubah

Nilai korelasi antar peubah adalah sebagai berikut.

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

Korelasi ini menunjukan bahwa y dan x1 memiliki korelasi positif yang sangat kuat dgn nilai 0.9. Sedangkan Y dengan x2 dan x3 memiliki korelasi yang kecil dengan masing-masing 0.31 dan -0.23

Sebaran Spasial Data Pengamatan

Sebaran Spasial y

k=101
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$y<- data$y
spplot(peta, "y", col.regions=color)

Peta peubah respon Y menunjukkan kemiripan antarkabupaten/kota terhadap sebaran peubah Y. Kemiripan ini ditunjukkan dengan pola kecenderungan untuk mengelompok dengan wilayah tetangganya sehingga terindikasi adanya autokorelasi spasial positif.

Matriks Bobot Spasial

Pemilihan matriks bobot didasarkan pada matriks bobot dengan nilai indeks moran yang signifikan pada taraf 5%

Matriks Kebertetanggaan

Pada Pemeriksaan Bobot Spasial berdasarkan matriks ketetanggaan terdapat masalah, yaitu ada 1 daerah yang tidak memiliki tetangga, sehingga diaktifkan opsi zero.policy = TRUE.

1 Queen Contiguity
# Queen Contiguity
sp.peta <- SpatialPolygons(peta@polygons)
qc <- poly2nb(sp.peta, queen = T)
qc
## 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
## 3 disjoint connected subgraphs
W.qc <- nb2listw(qc, style='W',zero.policy=TRUE)
qct = moran.test(data$y, W.qc,randomisation=T,
alternative="greater", zero.policy = TRUE)
qct
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.qc  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.85759, p-value = 0.1956
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.049643472      -0.008547009       0.004604072
moran.mc(data$y, listw = W.qc, zero.policy = TRUE, nsim = 99)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$y 
## weights: W.qc  
## number of simulations + 1: 100 
## 
## statistic = 0.049643, observed rank = 78, p-value = 0.22
## alternative hypothesis: greater
2 Rook Contiguity
# Rook Contiguity
sp.peta <- SpatialPolygons(peta@polygons)
rc <- poly2nb(sp.peta, queen = F)
rc
## 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:
## 0
## 3 disjoint connected subgraphs
W.rc <- nb2listw(rc, style='W',zero.policy=TRUE)
rct = moran.test(data$y, W.rc,randomisation=T,
alternative="greater", zero.policy = TRUE)
rct
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.rc  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 0.7872, p-value = 0.2156
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.044981054      -0.008547009       0.004623737
moran.mc(data$y, listw = W.rc, zero.policy = TRUE, nsim = 99)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$y 
## weights: W.rc  
## number of simulations + 1: 100 
## 
## statistic = 0.044981, observed rank = 79, p-value = 0.21
## alternative hypothesis: greater

Matriks 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
peta$long <- longlat[,1]
peta$lat <- longlat[,2]
coords <- peta[c("long","lat")]
#class(coords)
koord <- as.data.frame(coords)
djarak<-dist(longlat)
m.djarak<-as.matrix(djarak)
1. K-NEAREST NEIGHBOR (KNN)
# k = 5 (Digunakan k tetangga terdekat adalah 5)
W.knn<-knn2nb(knearneigh(longlat,k=5,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')
MI.knn <- moran(data$y,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(data$y,W.knn.s,randomisation = T,alternative = "greater")
mt1
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = -0.2758, p-value = 0.6086
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.023212662      -0.008474576       0.002855605
2. Radial Distance Weigth (Bobot Jarak Ambang)
# d=60 ditetapkan nilai ambang dmax adalah 60 km
W.dmax<-dnearneigh(longlat,0,60,longlat=TRUE)
W.dmax.s <- nb2listw(W.dmax,style='W')
MI.dmax <- moran(data$y,W.dmax.s,n=length(W.dmax.s$neighbours),S0=Szero(W.dmax.s))
mt2 = moran.test(data$y,W.dmax.s,randomisation = T,alternative = "greater")
mt2
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.dmax.s    
## 
## Moran I statistic standard deviate = 0.15641, p-value = 0.4379
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0005773405     -0.0084745763      0.0025494618
3. Power Distance Weigth / Invers Distance Weight (Bobot Jarak Invers)
# Untuk Alpha = 1
alpha1=1
alpha1 <- 1
W.idw <- 1/(m.djarak^alpha1)
diag(W.idw)<-0
rtot<-rowSums(W.idw,na.rm=TRUE)
W.idw.sd<-W.idw/rtot #row-normalized
W.idw.s = mat2listw(W.idw.sd,style='W')
MI.idw <- moran(data$y,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(data$y,W.idw.s,randomisation = T,alternative = "greater")
mt3
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 1.6573, p-value = 0.04873
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0258327380     -0.0084745763      0.0004285004
# Untuk alpha = 2
alpha2=2
W.idw2<-1/(m.djarak^alpha2)
diag(W.idw2)<-0
rtot<-rowSums(W.idw2,na.rm=TRUE)
W.idw2.sd<-W.idw2/rtot #row-normalized
W.idw2.s = mat2listw(W.idw2.sd,style='W')
MI.idw2 <- moran(data$y,W.idw2.s,n=length(W.idw2.s$neighbours),S0=Szero(W.idw2.s))
mt4 = moran.test(data$y,W.idw2.s,randomisation = T,alternative = "greater")
mt4
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 0.49133, p-value = 0.3116
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.015975652      -0.008474576       0.002476387
4. Exponential Distance Weigth (Bobot Jarak Eksponensial)
# Untuk alpha = 1
alpha=1
W.e<-exp((-alpha)*m.djarak)
diag(W.e)<-0
rtot<-rowSums(W.e,na.rm=TRUE)
W.e.sd<-W.e/rtot #row-normalized
W.e.s = mat2listw(W.e.sd,style='W')
MI.e <- moran(data$y,W.e.s,n=length(W.e.s$neighbours),S0=Szero(W.e.s))
mt5 = moran.test(data$y,W.e.s,randomisation = T,alternative = "greater")
mt5
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.e.s    
## 
## Moran I statistic standard deviate = -0.19817, p-value = 0.5785
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0112712359     -0.0084745763      0.0001991515
# Untuk alpha = 2
alpha2=2
W.e2<-exp((-alpha2)*m.djarak)
diag(W.e2)<-0
rtot<-rowSums(W.e2,na.rm=TRUE)
W.e2.sd<-W.e2/rtot #row-normalized
W.e2.s = mat2listw(W.e2.sd,style='W')
MI.e2 <- moran(data$y,W.e2.s,n=length(W.e2.s$neighbours),S0=Szero(W.e2.s))
mt6 = moran.test(data$y,W.e2.s,randomisation = T,alternative = "greater")
mt6
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.e2.s    
## 
## Moran I statistic standard deviate = -0.20904, p-value = 0.5828
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0130862592     -0.0084745763      0.0004866963

Matriks Bobot Terpilih

MatrikBobot <- c("KNN","dmax","Idw1","Idw2","Exp1","Exp2","Rook","Queen")
IndeksMoran <- c(MI.knn$I,MI.dmax$I,MI.idw$I,MI.idw2$I,MI.e$I,MI.e2$I,rct$estimate[1],qct$estimate[1])
pv = c(mt1$p.value,mt2$p.value,mt3$p.value,mt4$p.value,mt5$p.value,mt6$p.value,rct$p.value,qct$p.value)
M = cbind.data.frame(MatrikBobot,IndeksMoran,"p-value"=pv)
M
##   MatrikBobot   IndeksMoran    p-value
## 1         KNN -0.0232126615 0.60864860
## 2        dmax -0.0005773405 0.43785687
## 3        Idw1  0.0258327380 0.04872558
## 4        Idw2  0.0159756520 0.31159640
## 5        Exp1 -0.0112712359 0.57854575
## 6        Exp2 -0.0130862592 0.58279169
## 7        Rook  0.0449810541 0.21558247
## 8       Queen  0.0496434720 0.19555880
M[M$`p-value`< 0.05,]
##   MatrikBobot IndeksMoran    p-value
## 3        Idw1  0.02583274 0.04872558

Matriks bobot yang memiliki nilai Indeks Moran yang signifikan pada taraf 5% adalah Invers Distance Weight (Bobot Jarak Invers) dengan alpha =1. Selanjutnya, matriks bobot ini yang akan digunakan dalam analisis.

Power Distance Weigth / Invers Distance Weight (Bobot Jarak Invers)

#Alpha = 1
alpha1=1
W.idw <-1/(m.djarak^alpha1)
diag(W.idw)<-0
rtot<-rowSums(W.idw,na.rm=TRUE)
W.idw.sd<-W.idw/rtot #row-normalized
W.idw.s = mat2listw(W.idw.sd,style='W')
MI.idw <- moran(data$y,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(data$y,W.idw.s,randomisation = T,alternative = "greater")
mt3
## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 1.6573, p-value = 0.04873
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0258327380     -0.0084745763      0.0004285004
ww = W.idw.s

Pemodelan Ordinary Least Square

Memeriksa keterpenuhan persyaratan multikolinearitas, serta menguji asumsi autokorelasi spasial, kehomogenan ragam, kenormalan sisaan, dan kesimpulan Matriks Bobot yang digunakan adalah Power Distance Weigth / Invers Distance Weight (Bobot Jarak Invers) dengan alpha = 1 sebagaimana hasil dari tahap sebelumnya.

Model Ordinary Least Square

ols <- lm(y~x1+x2+x3, data=data)

Pemeriksaan Multikolinearitas

Pada regresi linear berganda perlu diperiksa persyaratan tidak adanya kolinearitas antar peubah penjelas agar koefisien regresi yang diperoleh dapat dikatakan valid. Pemeriksaan multikolinearitas dapat dilakukan dengan menggunakan nilai variance inflation factor (VIF). Multikolinearitas terjadi jika nilai VIF > 5.

car::vif(ols)
##       x1       x2       x3 
## 1.051166 1.077603 1.045575

Berdasarkan perhitungan, diperoleh nilai VIF < 5 yang menunjukkan tidak terjadi multikolinearitas antar peubah penjelas.

ols <- lm(y~x1+x2+x3, data=data)
summary(ols)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.469 -14.140  -3.913  14.493  70.459 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  70.94445   11.86595   5.979 2.59e-08 ***
## x1            0.60524    0.02682  22.570  < 2e-16 ***
## x2            2.97143    1.07142   2.773  0.00648 ** 
## x3          -14.01285    4.76881  -2.938  0.00399 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.24 on 115 degrees of freedom
## Multiple R-squared:  0.8391, Adjusted R-squared:  0.8349 
## F-statistic: 199.9 on 3 and 115 DF,  p-value: < 2.2e-16
AIC(ols)
## [1] 1092.336

Uji Asumsi

Uji Normalitas Sisaan

H0 : galat model menyebar normal H1 : galat model tidak menyebar normal

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
err.ols <- residuals(ols)
library(nortest)
ad.test(err.ols)
## 
##  Anderson-Darling normality test
## 
## data:  err.ols
## A = 0.47582, p-value = 0.235

Karena p-value > 5%, maka GAGAL TOLAK H0. Artinya sisaan menyebar NORMAL.

Uji Autokorelasi Spasial

Pengujian autokorelasi spasial menggunakan Indeks Moran yang memerlukan matriks pembobot spasial. Dalam kasus ini, digunakan matriks bobot terbaik yang telah diperoleh sebelumnya, yaitu Invers Distance Weight (Bobot Jarak Invers) dengan alpha = 1. H0 : Tidak Ada Autokorelasi H1 : Ada Autokorelasi

ww = W.idw.s
lm.morantest(ols, listw=ww, alternative="two.sided", zero.policy = TRUE)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 2.1803, p-value = 0.02924
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.036234217     -0.008645290      0.000423707

Karena p-value < 5%, maka TOLAK H0. Artinya terdapat AUTOKORELASI SPASIAL pada data.

Uji Kehomogenan Ragam

H0 : Ragam Sisaan Homogen H1 : Ragam Sisaan Tidak Homogen

lmtest::bptest(ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 0.19369, df = 3, p-value = 0.9786

Karena p-value > 5%, maka GAGAL TOLAK H0. Artinya ragam sisaan homogen.

Kesimpulan

Dari analsis terhadap OLS di atas, disimpulkan bahwa terdapat Autokorelasi Spasial, artinya model OLS masih belum baik karena terdapat asumsi yang tidak terpenuhi, sehingga perlu dilanjutkan dengan model dependensi spasial untuk mencari model yang lebih baik.

Pemodelan Dependensi Spasial

Pemilihan model terbaik berdasarkan AIC, menguji asumsi kenormalan sisaan, kehomogenan ragam, autokorelasi spasial, serta interpretasi model. Pada tahap sebelumnya, diketahui bahwa terjadi pelanggaran Asumsi kebebasan sisaan pada model regresi global di mana terjadi Autokorelasi Spasial pada sisaan. Untuk mengatasi hal ini, perlu dilakukan pemodelan dependensi spasial. Matriks Bobot yang digunakan adalah Power Distance Weigth / Invers Distance Weight (Bobot Jarak Invers) dengan alpha =

Model Dependensi Spasial

ols <- lm(y~x1+x2+x3, data=data)

UJI LM

Uji LM untuk mengetahui model dependensi mana yang mungkin dipilih.

ww = W.idw.s
model <- lm.LMtests(ols,listw=ww,zero.policy = TRUE, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
summary(model)
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## test weights: listw
##  
##          statistic parameter  p.value   
## RSerr     2.299700         1 0.129399   
## RSlag     7.793113         1 0.005245 **
## adjRSerr  0.050063         1 0.822953   
## adjRSlag  5.543477         1 0.018550 * 
## SARMA     7.843177         2 0.019810 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pada Uji LM pada taraf 5%, diketahui nyata pada LMlag, RLlag, dan SARMA, sehingga kemungkinan modelnya adalah SAR atau SARMA.

SAR
sar <- lagsarlm(y ~ x1 + x2+x3, data = data, listw = ww)
summary(sar, Nagelkerke = T)
## 
## Call:lagsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -66.2865 -13.5852  -4.4545  13.7061  60.9695 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  19.87310   28.54442  0.6962  0.486293
## x1            0.60001    0.02588 23.1844 < 2.2e-16
## x2            3.36011    1.03950  3.2324  0.001227
## x3          -13.83058    4.59155 -3.0122  0.002594
## 
## Rho: 0.33585, LR test value: 4.5098, p-value: 0.033702
## Asymptotic standard error: 0.17557
##     z-value: 1.9129, p-value: 0.055766
## Wald statistic: 3.659, p-value: 0.055766
## 
## Log likelihood: -538.9132 for lag model
## ML residual variance (sigma squared): 500.32, (sigma: 22.368)
## Nagelkerke pseudo-R-squared: 0.8451 
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 1089.8, (AIC for lm: 1092.3)
## LM test for residual autocorrelation
## test value: 2.0286, p-value: 0.15436

Output di atas menunjukkan bahwa koefisien Rho signifikan pada taraf nyata 5% (p-value 0.033702). AIC model SAR adalah sebesar 1089.8.

SARMA
SARMA <- sacsarlm(y ~ x1 + x2+x3,data=data,ww, zero.policy = TRUE)
summary(SARMA)
## 
## Call:sacsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -68.04274 -13.42963  -0.67171  13.85824  59.32255 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  -2.244431  20.811651 -0.1078 0.9141187
## x1            0.595749   0.025549 23.3183 < 2.2e-16
## x2            3.120825   1.006513  3.1006 0.0019311
## x3          -15.710042   4.543932 -3.4574 0.0005455
## 
## Rho: 0.51918
## Asymptotic standard error: 0.12137
##     z-value: 4.2777, p-value: 1.8884e-05
## Lambda: -0.78421
## Asymptotic standard error: 0.31481
##     z-value: -2.4911, p-value: 0.012735
## 
## LR test value: 7.7221, p-value: 0.021046
## 
## Log likelihood: -537.3071 for sac model
## ML residual variance (sigma squared): 474.25, (sigma: 21.777)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1088.6, (AIC for lm: 1092.3)

Output di atas menunjukkan bahwa koefisien Rho dan Lambda signifikan pada taraf nyata 5%. AIC model SARMA adalah sebesar 1088.6.

Model Terbaik
cbind.data.frame(AIC(ols),AIC(sar),AIC(SARMA))
##   AIC(ols) AIC(sar) AIC(SARMA)
## 1 1092.336 1089.826   1088.614

Dari kedua model dependensi di atas model SARMA memiliki AIC yang lebih kecil, sehingga model SARMA dipilih menjadi model terbaik.

Uji Asumsi Model
library(nortest)
err.SARMA<-residuals(SARMA)
Kenormalan Sisaan

H0 : galat model menyebar normal H1 : galat model tidak menyebar normal

ad.test(err.SARMA)
## 
##  Anderson-Darling normality test
## 
## data:  err.SARMA
## A = 0.48107, p-value = 0.2282

Karena p-value > 5%, maka GAGAL TOLAK H0. Artinya sisaan menyebar normal.

Kehomogenan Ragam

H0 : Ragam Sisaan Homogen H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sar)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.8048, df = 3, p-value = 0.6139

Karena p-value > 5%, maka GAGAL TOLAK H0. Artinya ragam sisaan homogen.

Uji Kebebasan Sisaan

H0 : Tidak Ada Autokorelasi H1 : Ada Autokorelasi

moran.test(err.SARMA, ww, alternative="two.sided", zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  err.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.037194, p-value = 0.9703
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0077079920     -0.0084745763      0.0004247816

Karena p-value > 5%, maka GAGAL TOLAK H0. Artinya Tidak ada korelasi Spasial. Karena semua asumsi telah terpenuhi, maka model ini telah valid. Selanjutnya, dilakukan interpretasi.

Interpretasi

Karena Model SARMA memiliki spillover, maka perlu diinterpretasi dengan memperhatikan spillovernya.

# Summary Model
S = summary(SARMA)
S
## 
## Call:sacsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -68.04274 -13.42963  -0.67171  13.85824  59.32255 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  -2.244431  20.811651 -0.1078 0.9141187
## x1            0.595749   0.025549 23.3183 < 2.2e-16
## x2            3.120825   1.006513  3.1006 0.0019311
## x3          -15.710042   4.543932 -3.4574 0.0005455
## 
## Rho: 0.51918
## Asymptotic standard error: 0.12137
##     z-value: 4.2777, p-value: 1.8884e-05
## Lambda: -0.78421
## Asymptotic standard error: 0.31481
##     z-value: -2.4911, p-value: 0.012735
## 
## LR test value: 7.7221, p-value: 0.021046
## 
## Log likelihood: -537.3071 for sac model
## ML residual variance (sigma squared): 474.25, (sigma: 21.777)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1088.6, (AIC for lm: 1092.3)
# Spill over
Im = impacts(SARMA, listw = ww)
Im
## Impact measures (sac, exact):
##         Direct   Indirect      Total
## x1   0.6033648   0.635669   1.239034
## x2   3.1607177   3.329942   6.490660
## x3 -15.9108610 -16.762728 -32.673589
# Efek Umpan Balik
koef = S$coefficients[-1]
diref = Im$direct
umbal = diref-koef
cbind.data.frame(Koefisien=koef,EfekLangsung=diref,UmpanBalik=umbal)
##      Koefisien EfekLangsung   UmpanBalik
## x1   0.5957494    0.6033648  0.007615374
## x2   3.1208246    3.1607177  0.039893023
## x3 -15.7100419  -15.9108610 -0.200819057

Berdasarkan hasil di atas menunjukan bahwa pengaruh total dari peubah x1 adalah sebesar 1.239034, artinya jika jumlah tenaga kerja naik 1 satuan (1000 orang), maka secara rata-rata y (PDRB kabupaten/kota) meningkat sebesar 1.239034 satuan jika peubah lain tetap. Efek langsung x1 sebesar 0.6033648, sedangkan koefisien SARMA x1 sebesar 0.595749, artinya efek umpan balik adalah sebesar 0.6033648 - 0.595749 = 0.0076158. Selain itu, pengaruh total dari peubah x2 adalah sebesar 6.490660, artinya jika pendapatan asli daerah naik 1 satuan (1 triliun), maka secara rata-rata y (PDRB kabupaten/kota) meningkat sebesar 6.490660 satuan jika peubah lain tetap. Efek langsung x2 sebesar 3.1607177, sedangkan koefisien SARMA x2 sebesar 3.120825, artinya efek umpan balik adalah sebesar 3.1607177 - 3.120825 = 0.0398927. Pengaruh total dari peubah x3 adalah sebesar -32.673589, artinya jika inflasi naik 1 satuan (1 persen), maka secara rata-rata y (PDRB kabupaten/kota) turun sebesar 32.673589 satuan jika peubah lain tetap. Efek langsung x3 sebesar -15.9108607, sedangkan koefisien SARMA x3 sebesar -15.710042, artinya efek umpan balik adalah sebesar -15.9108607 - (-15.710042) = -0.2008187.