Cek Missing Value

data[!complete.cases(data),]

## [1] V1     ID_KAB KAB    y      x1     x2     x3    
## <0 rows> (or 0-length row.names)

Data yang digunakan memiliki 119 baris, dan tidak ada missing data.

Ringkasan Data

summary(data[,-c(1,2)])

##      KAB                  y               x1               x2       
##  Length:119         Min.   :207.0   Min.   : 16.25   Min.   :0.110  
##  Class :character   1st Qu.:265.6   1st Qu.: 69.83   1st Qu.:2.010  
##  Mode  :character   Median :288.7   Median :136.01   Median :3.700  
##                     Mean   :292.6   Mean   :144.65   Mean   :3.802  
##                     3rd Qu.:315.3   3rd Qu.:217.48   3rd Qu.:5.590  
##                     Max.   :402.9   Max.   :298.31   Max.   :7.120  
##        x3       
##  Min.   :1.180  
##  1st Qu.:1.595  
##  Median :2.020  
##  Mean   :2.000  
##  3rd Qu.:2.380  
##  Max.   :2.780

Boxplot

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

Dari Boxplot keempat peubah, terdapat data yang merupakan pencilan atas dari peubah y.

Scatter Plot

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

Dari scatter plot di atas, tampak 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

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

corrplot(cor(data[,4:7]), method = "number", is.cor = T, type = "lower", diag=F)

corrplot(cor(data[,4:7]), method = "square", bg = "#F7EFE5",tl.col = "brown",rect.col = "black", pch.col = "black", type = "lower")
corrplot(cor(data[,4:7]), method = "number", bg = "#F7EFE5",tl.col = "brown",rect.col = "black", pch.col = "black", type = "upper", add = TRUE,tl.pos="lt")

Terlihat bahwa peubah y dan x1 memiliki korelasi positif yang cukup kuat dengan nilai 0.53. Sedangkan peubah Y dengan x2 dan x3 memiliki korelasi yang kecil dengan masing-masing 0.23 dan -0.02

Sebaran Spasial Data Pengamatan

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

Kabupaten/Kota yang memiliki PDRB yang tinggi saling berdekatan dengan Kabupaten/Kota yang memiliki PDRB yang tinggi juga. Hal ini menunjukkan ada pola spasial pada variabel tersebut.

Peta tematik peubah respon Y (PDRB) menunjukkan kemiripan antar kabupaten/kota terhadap sebaran peubah Y (PDRB). Kemiripan ini ditunjukkan dengan pola kecenderungan untuk mengelompok dengan wilayah tetangganya sehingga terindikasi adanya autokorelasi spasial positif.

Matriks Ketetanggaan

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

W.qc <- nb2listw(qc, style='W',zero.policy=TRUE)
ols <- lm(y~x1+x2+x3, data=data)
qct = lm.morantest(ols, W.qc, alternative="greater", zero.policy = TRUE)
qct

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.qc
## 
## Moran I statistic standard deviate = 4.8151, p-value = 7.355e-07
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.317462275     -0.008104406      0.004571536

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

W.rc <- nb2listw(rc, style='W',zero.policy=TRUE)
rct = lm.morantest(ols, W.rc, alternative="greater", zero.policy = TRUE)
rct

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.rc
## 
## Moran I statistic standard deviate = 4.7796, p-value = 8.781e-07
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.315616525     -0.008156051      0.004588711

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)

# k = 5
W.knn<-knn2nb(knearneigh(longlat,k=5,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')
mt1 = lm.morantest(ols,W.knn.s,alternative = "greater")
mt1

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.knn.s
## 
## Moran I statistic standard deviate = 7.4153, p-value = 6.066e-14
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.385620363     -0.008836269      0.002829678

W.dmax<-dnearneigh(longlat,0,60,longlat=TRUE)
W.dmax.s <- nb2listw(W.dmax,style='W')
mt2 = lm.morantest(ols,W.dmax.s,alternative = "greater")
mt2

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.dmax.s
## 
## Moran I statistic standard deviate = 8.8702, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.437564040     -0.008595749      0.002529937

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')
mt3 = lm.morantest(ols,W.e.s,alternative ="greater")
mt3

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.e.s
## 
## Moran I statistic standard deviate = 24.598, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.3345470268    -0.0085773304     0.0001945776

alpha2=2
W.e2 <- exp((-alpha2)*m.djarak)
diag(W.e2) <- 0
rtot2 <- rowSums(W.e2,na.rm=TRUE)
W.e2.sd <- W.e2/rtot2 #row-normalized
W.e2.s = mat2listw(W.e2.sd,style='W')
mt4 = lm.morantest(ols,W.e2.s,alternative ="greater")
mt4

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.e2.s
## 
## Moran I statistic standard deviate = 18.159, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.3898383209    -0.0085581075     0.0004813406

Matriks Bobot Terpilih

MatrikBobot <- c("KNN","dmax","Exp1","Exp2","Rook","Queen")
IndeksMoran <- c(mt1$estimate[1],mt2$estimate[1],mt3$estimate[1],mt4$estimate[1],rct$estimate[1],qct$estimate[1])
pv = c(mt1$p.value,mt2$p.value,mt3$p.value,mt4$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.3856204  6.065876e-14
## 2        dmax   0.4375640  3.649261e-19
## 3        Exp1   0.3345470 6.586173e-134
## 4        Exp2   0.3898383  5.460913e-74
## 5        Rook   0.3156165  8.780689e-07
## 6       Queen   0.3174623  7.354846e-07

Matriks bobot yang memiliki nilai Indeks Moran yang tertinggi adalah Radial Distance Weigth. Selanjutnya, matriks bobot ini yang akan digunakan dalam analisis.

W.dmax<-dnearneigh(longlat,0,60,longlat=TRUE)
W.dmax.s <- nb2listw(W.dmax,style='W')
mt2 = lm.morantest(ols,W.dmax.s,alternative = "greater")
mt2

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.dmax.s
## 
## Moran I statistic standard deviate = 8.8702, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.437564040     -0.008595749      0.002529937

ww = W.dmax.s
ww

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 994 
## Percentage nonzero weights: 7.019278 
## Average number of links: 8.352941 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1      S2
## W 119 14161 119 38.00293 485.358

Tahap kedua dilakukan pemodelan OLS untuk memeriksa keterpenuhan persyaratan multikolinieritas serta menguji asumsi yang ada di regresi OLS yaitu autokorelasi spasial, ragam galat homogen, kenormalan galat.

Matriks Bobot yang digunakan adalah Radial Distance Weigth sebagaimana hasil dari tahap sebelumnya.

Pemeriksaan multikolinieritas

Pada regresi OLS 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.001146 1.000421 1.000836

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

Uji Asumsi

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.3

## 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.61211, p-value = 0.109

ww = W.dmax.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 = 8.8702, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.437564040     -0.008595749      0.002529937

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 1.6864, df = 3, p-value = 0.64

Kesimpulan

Dari analsis terhadap model OLS, disimpulkan bahwa terjadi pelanggaran Asumsi kebebasan sisaan pada model regresi global di mana terjadi Autokorelasi Spasial pada sisaan. 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. Matriks Bobot yang digunakan adalah Radial Distance Weigth sebagaimana hasil dari tahap sebelumnya.

Tahap 3 dilakukan Pemodelan Dependensi Spasial, Pemilihan model terbaik berdasarkan AIC, Pengujian asumsi kenormalan sisaan, kehomogenan ragam, autokorelasi spasial, serta Interpretasi model terbaik.

Spatial Lag X (SLX)

Jika semua variabel x signifikan, maka bisa jadi model terbaiknya SLX atau SDM atau SDEM. Tapi jika tidak semua variabel x signifikan, maka lanjut uji LM untuk menentukan model lain (SAR, SEM, atau SARMA).

# Matriks pembobot: Radial Distance Weigth

library(spdep)
W.dmax<-dnearneigh(longlat,0,60,longlat=TRUE)
W.dmax.s <- nb2listw(W.dmax,style='W')
mt2 = lm.morantest(ols,W.dmax.s,alternative = "greater")
mt2

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: W.dmax.s
## 
## Moran I statistic standard deviate = 8.8702, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.437564040     -0.008595749      0.002529937

Didapatkan indeks moran sebesar 0.437564040 dengan p-value sebesar < 2.2e-16.

W.qc <- nb2listw(qc, style='W',zero.policy=TRUE)
qcty = moran.test(data$y, W.qc,alternative="greater", zero.policy = TRUE)
qcty

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.qc  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 5.2015, p-value = 9.886e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.343031012      -0.008547009       0.004568695

qctx1 = moran.test(data$x1, W.qc,alternative="greater", zero.policy = TRUE)
qctx1

## 
##  Moran I test under randomisation
## 
## data:  data$x1  
## weights: W.qc  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = -0.39643, p-value = 0.6541
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.035507683      -0.008547009       0.004625208

qctx2 = moran.test(data$x2, W.qc,alternative="greater", zero.policy = TRUE)
qctx2

## 
##  Moran I test under randomisation
## 
## data:  data$x2  
## weights: W.qc  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 0.71224, p-value = 0.2382
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.039909439      -0.008547009       0.004628657

qctx3 = moran.test(data$x3, W.qc,alternative="greater", zero.policy = TRUE)
qctx3

## 
##  Moran I test under randomisation
## 
## data:  data$x3  
## weights: W.qc  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = -0.80605, p-value = 0.7899
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.063345032      -0.008547009       0.004621696

SLX <- lmSLX(y ~ x1+x2+x3, data=data, listw=ww, 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 
## -69.299 -17.763  -0.666  18.579  71.004 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 192.67974   40.16175   4.798 4.99e-06 ***
## x1            0.23518    0.03277   7.176 8.30e-11 ***
## x2            3.50935    1.35192   2.596  0.01070 *  
## x3           -1.41065    6.22932  -0.226  0.82126    
## lag.x1        0.20295    0.08662   2.343  0.02090 *  
## lag.x2        8.81718    3.31131   2.663  0.00889 ** 
## lag.x3       -3.98719   15.77224  -0.253  0.80089    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.27 on 112 degrees of freedom
## Multiple R-squared:  0.414,  Adjusted R-squared:  0.3826 
## F-statistic: 13.19 on 6 and 112 DF,  p-value: 2.992e-11

Variabel <- c("Y","X1","X2","X3")
IndeksMoran <- c(qcty$estimate[1],qctx1$estimate[1],qctx2$estimate[1],qctx3$estimate[1])
pv = c(qcty$p.value,qctx1$p.value,qctx2$p.value,qctx3$p.value)
M = cbind.data.frame(Variabel,IndeksMoran,"p-value"=pv)
M

##   Variabel IndeksMoran      p-value
## 1        Y  0.34303101 9.886209e-08
## 2       X1 -0.03550768 6.541056e-01
## 3       X2  0.03990944 2.381593e-01
## 4       X3 -0.06334503 7.898942e-01

Pada uji Morans’I tidak ada variabel x yang terdapat autokorelasi atau keterkaitan antar wilayah Kabupaten/Kota di Pulau Jawa. Berdasarkan nilai p-value, tidak ada p-value yang kurang dari \(\alpha\). Hal ini menandakan bahwa dependensi spasial tidak ada di x, melainkan ada di lag (y) atau error. Karena itulah, akan dilanjutkan dengan uji LM.

Spatial Durbin Model (SDM)

Tambahan aja sebelum ke uji LM

SDM <- lagsarlm(y ~ x1+x2+x3, data=data, listw=ww, type="mixed");
summary(SDM)

## 
## Call:lagsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww, 
##     type = "mixed")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -52.5558 -15.1988  -2.6865  15.1846  58.5746 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 56.874320  38.687351  1.4701   0.14153
## x1           0.212545   0.026843  7.9180 2.442e-15
## x2           2.809180   1.103939  2.5447   0.01094
## x3          -1.190135   5.073011 -0.2346   0.81452
## lag.x1       0.015074   0.077056  0.1956   0.84490
## lag.x2       1.828848   2.768484  0.6606   0.50887
## lag.x3       6.918889  12.845126  0.5386   0.59014
## 
## Rho: 0.59238, LR test value: 33.705, p-value: 6.4122e-09
## Asymptotic standard error: 0.095764
##     z-value: 6.1858, p-value: 6.1772e-10
## Wald statistic: 38.265, p-value: 6.1772e-10
## 
## Log likelihood: -554.21 for mixed model
## ML residual variance (sigma squared): 607.55, (sigma: 24.648)
## Number of observations: 119 
## Number of parameters estimated: 9 
## AIC: 1126.4, (AIC for lm: 1158.1)
## LM test for residual autocorrelation
## test value: 8.9152, p-value: 0.0028281

Output di atas menunjukkan bahwa koefisien rho signifikan pada taraf nyata 5% (p-value = 6.4122e-09). AIC model SDM adalah sebesar 1126.4

Spatial Durbin Error Model (SDEM)

SDEM <- errorsarlm(y ~ x1+x2+x3, data=data, listw=ww, etype="mixed");
summary(SDEM)

## 
## Call:errorsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww, 
##     etype = "mixed")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -55.8418 -16.5728  -1.2647  13.8470  64.3167 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 187.162859  41.142822  4.5491 5.388e-06
## x1            0.220925   0.029045  7.6063 2.820e-14
## x2            3.057088   1.187782  2.5738   0.01006
## x3            0.774772   5.579618  0.1389   0.88956
## lag.x1        0.134598   0.084275  1.5971   0.11024
## lag.x2        3.358980   3.544103  0.9478   0.34325
## lag.x3       13.506781  15.254057  0.8855   0.37591
## 
## Lambda: 0.62641, LR test value: 32.59, p-value: 1.138e-08
## Asymptotic standard error: 0.092574
##     z-value: 6.7665, p-value: 1.3189e-11
## Wald statistic: 45.786, p-value: 1.3189e-11
## 
## Log likelihood: -554.7678 for error model
## ML residual variance (sigma squared): 607.23, (sigma: 24.642)
## Number of observations: 119 
## Number of parameters estimated: 9 
## AIC: NA (not available for weighted model), (AIC for lm: 1158.1)

Output di atas menunjukkan bahwa koefisien lambda signifikan pada taraf nyata 5% (p-value = 1.138e-08). AIC model SDM adalah sebesar 1158.1.

Uji LM (LM dan robust LM)

Uji LM untuk mengetahui model dependensi mana yang mungkin dipilih.

ols <- lm(y~x1+x2+x3, data=data)
ww = W.dmax.s
model <- lm.LMtests(ols,listw=ww,zero.policy = TRUE, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
summary(model)

##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = y ~ x1 + x2 + x3, data = data)
## weights: ww
##  
##        statistic parameter   p.value    
## LMerr   71.34443         1 < 2.2e-16 ***
## RLMerr   0.48712         1 0.4852150    
## LMlag   83.74747         1 < 2.2e-16 ***
## RLMlag  12.89016         1 0.0003303 ***
## SARMA   84.23459         2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keterangan: - LMerr: LM error - RLMerr: Robust LM error - LMlag: LM lag (respon) - RLMlag: Robust LM lag (respon) - SARMA

Pada uji LM pada taraf nyata 5%, diketahui signifikan pada LM error dan LM lag. Pada uji Robust LM, diketahui signifikan pada RLM lag. Selain itu, uji SARMA juga signifikan.

Sehingga kemungkinan modelnya adalah SAR atau SARMA.SARMA dimasukkan karena di uji LM lag dan LM errornya signifikan. Namun, akan dilakukan juga SEM untuk lebih memastikannya.

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 
## -52.3445 -16.6792  -2.9435  14.5660  57.1540 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 74.877606  26.558139  2.8194  0.004812
## x1           0.213258   0.026652  8.0015 1.332e-15
## x2           2.754510   1.084863  2.5390  0.011116
## x3          -1.930412   4.959795 -0.3892  0.697119
## 
## Rho: 0.61526, LR test value: 48.995, p-value: 2.5659e-12
## Asymptotic standard error: 0.083287
##     z-value: 7.3872, p-value: 1.4988e-13
## Wald statistic: 54.571, p-value: 1.4999e-13
## 
## Log likelihood: -554.5723 for lag model
## ML residual variance (sigma squared): 607.27, (sigma: 24.643)
## Nagelkerke pseudo-R-squared: 0.55585 
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 1121.1, (AIC for lm: 1168.1)
## LM test for residual autocorrelation
## test value: 6.894, p-value: 0.0086487

Output di atas menunjukkan bahwa koefisien Rho signifikan pada taraf nyata 5% (p-value = 2.5659e-12). AIC model SAR adalah sebesar 1121.1

SEM

SEM <- errorsarlm(y ~ x1+x2+x3, data= data, listw = ww)
summary(SEM)

## 
## Call:errorsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -58.3396 -15.9169  -2.2299  15.0532  68.3221 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 260.332754  12.483414 20.8543 < 2.2e-16
## x1            0.197085   0.026087  7.5550 4.197e-14
## x2            2.391002   1.099753  2.1741    0.0297
## x3           -2.910802   4.729261 -0.6155    0.5382
## 
## Lambda: 0.64855, LR test value: 43.163, p-value: 5.0357e-11
## Asymptotic standard error: 0.088781
##     z-value: 7.305, p-value: 2.7733e-13
## Wald statistic: 53.363, p-value: 2.7722e-13
## 
## Log likelihood: -557.4882 for error model
## ML residual variance (sigma squared): 631.18, (sigma: 25.123)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model), (AIC for lm: 1168.1)

Output di atas menunjukkan bahwa koefisien lambda signifikan pada taraf nyata 5% (p-value = 5.03573e-11). AIC model SEM adalah sebesar 1127

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 
## -50.1070 -16.0771  -1.0449  15.1752  56.9846 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 22.403963  18.249833  1.2276 0.219588
## x1           0.202330   0.024639  8.2116 2.22e-16
## x2           2.448615   0.946076  2.5882 0.009648
## x3          -1.202825   4.515978 -0.2663 0.789971
## 
## Rho: 0.79951
## Asymptotic standard error: 0.056396
##     z-value: 14.177, p-value: < 2.22e-16
## Lambda: -0.62969
## Asymptotic standard error: 0.17554
##     z-value: -3.5871, p-value: 0.00033442
## 
## LR test value: 56.665, p-value: 4.9583e-13
## 
## Log likelihood: -550.7373 for sac model
## ML residual variance (sigma squared): 499.41, (sigma: 22.347)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1115.5, (AIC for lm: 1168.1)

Output di atas menunjukkan bahwa koefisien Rho (p-value < 2.22e-16) dan Lambda (p-value = 0.00033442) signifikan pada taraf nyata 5%. AIC model SARMA adalah sebesar 1115.5.

Uji Asumsi Model Terbaik

library(nortest)
residual.SARMA <- residuals(SARMA)

ad.test(residual.SARMA)

## 
##  Anderson-Darling normality test
## 
## data:  residual.SARMA
## A = 0.19481, p-value = 0.8899

bptest.Sarlm(SARMA)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 5.3496, df = 3, p-value = 0.1479

moran.test(residual.SARMA, ww, alternative="two.sided", zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  residual.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.18252, p-value = 0.8552
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0007282153     -0.0084745763      0.0025423833

Interpretasi Model Terbaik

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

catatan: Jika model terbaiknya SEM maka tidak memiliki spillover, karena SEM mirip seperti regresi biasa.

S = summary(SARMA)
S

## 
## Call:sacsarlm(formula = y ~ x1 + x2 + x3, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -50.1070 -16.0771  -1.0449  15.1752  56.9846 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 22.403963  18.249833  1.2276 0.219588
## x1           0.202330   0.024639  8.2116 2.22e-16
## x2           2.448615   0.946076  2.5882 0.009648
## x3          -1.202825   4.515978 -0.2663 0.789971
## 
## Rho: 0.79951
## Asymptotic standard error: 0.056396
##     z-value: 14.177, p-value: < 2.22e-16
## Lambda: -0.62969
## Asymptotic standard error: 0.17554
##     z-value: -3.5871, p-value: 0.00033442
## 
## LR test value: 56.665, p-value: 4.9583e-13
## 
## Log likelihood: -550.7373 for sac model
## ML residual variance (sigma squared): 499.41, (sigma: 22.347)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 1115.5, (AIC for lm: 1168.1)

# Spill over
Im = impacts(SARMA, listw = ww)
Im

## Impact measures (sac, exact):
##        Direct   Indirect     Total
## x1  0.2480898  0.7610841  1.009174
## x2  3.0023986  9.2106887 12.213087
## x3 -1.4748583 -4.5245362 -5.999394

# 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.2023304    0.2480898  0.0457594
## x2  2.4486154    3.0023986  0.5537832
## x3 -1.2028253   -1.4748583 -0.2720331

Terlihat bahwa pengaruh total dari peubah x1 adalah sebesar 1.009174, artinya jika x1 (jumlah tenaga kerja) naik 1 satuan (1000 orang), maka secara rata-rata y (PDRB kabupaten/kota) meningkat sebesar 1.009174 satuan jika peubah lain tetap. Efek langsung x1 sebesar 0.240898, sedangkan koefisien SARMA x1 sebesar 0.202330, artinya efek umpan balik adalah sebesar 0.240898 - 0.202330 = 0.0456594.

Terlihat bahwa pengaruh total dari peubah x2 adalah sebesar 12.213087, artinya jika x2 (pendapatan asli daerah) naik 1 satuan (1 triliun), maka secara rata-rata y (PDRB kabupaten/kota) meningkat sebesar 12.213087 satuan jika peubah lain tetap. Efek langsung x2 sebesar 3.0023986, sedangkan koefisien SARMA x2 sebesar 2.4486154, artinya efek umpan balik adalah sebesar 3.0023986 - 2.4486154 = 0.5537832.

Terlihat bahwa pengaruh total dari peubah x3 adalah sebesar -5.999394, artinya jika x3 (inflasi) naik 1 satuan (1 persen), maka secara rata-rata y (PDRB kabupaten/kota) turun sebesar 5.999394 satuan jika peubah lain tetap. Efek langsung x3 sebesar -1.4748583, sedangkan koefisien SARMA x3 sebesar -1.2028253, artinya efek umpan balik adalah sebesar -1.4748583 - (-1.2028253) = -0.2720331.