Regresi Spasial
Regresi spasial merupakan hasil pengembangan dari metode regresi klasik. Pengembangan itu berdasarkan adanya pengaruh tempat atau pengaruh spasial pada data yang dianalisis. Adanya efek spasial merupakan hal yang lazim terjadi antara satu wilayah dengan wilayah yang lain. Pada beberapa kasus, peubah tak bebas yang diamati memiliki keterkaitan dengan hasil pengamatan di wilayah yang berbeda, terutama wilayah yang berdekatan. Adanya hubungan spasial dalam peubah tak bebas akan menyebabkan pendugaan menjadi tidak tepat karena asumsi keacakan galat dilanggar. Untuk mengatasi permasalahan di atas diperlukan suatu model regresi yang memasukkan hubungan spasial antar wilayah ke dalam model. Adanya informasi hubungan spasial antar wilayah menyebabkan perlu mengakomodir keragaman spasial ke dalam model, sehingga model yang digunakan adalah model regresi spasial.
Berikut adalah gambaran umum mengenai Analisis Regresi Spasial
Package yang Digunakan
library(readxl) #membaca file excel
library(rgdal) #membaca file shp, fs readOGR
library(raster) #eksplorasi data
library(spdep) #pembobot spasial (fungsi poly2nb)
library(lmtest) #regresi klasik
library(nortest) #uji kenormalan galat
library(DescTools) #uji kebebasan galat
library(spatialreg) #model SAR
library(RColorBrewer) #memberi warna peta
library(car) #uji multikolinieritas
library(GWmodel) #GWR Jarak SpasialInput Data
Data yang akan diinput berupa data frame dan data polygon Jawa Barat
Data Frame
Berikut adalah data yang digunakan :
## # A tibble: 6 x 32
## PROVNO KABKOTNO KODE2010 PROVINSI KABKOT IDSP2010 Long Lat p.miskin15
## <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 32 1 3201 JAWA BA~ BOGOR 3201 107. -6.56 8.96
## 2 32 2 3202 JAWA BA~ SUKAB~ 3202 107. -7.07 8.96
## 3 32 3 3203 JAWA BA~ CIANJ~ 3203 107. -7.13 12.2
## 4 32 4 3204 JAWA BA~ BANDU~ 3204 108. -7.10 8.00
## 5 32 5 3205 JAWA BA~ GARUT 3205 108. -7.36 12.8
## 6 32 6 3206 JAWA BA~ TASIK~ 3206 108. -7.50 12.0
## # ... with 23 more variables: p.miskin16 <dbl>, j.miskin15 <dbl>,
## # j.miskin16 <dbl>, AHH2015 <dbl>, AHH2016 <dbl>, EYS2015 <dbl>,
## # EYS2016 <dbl>, MYS2015 <dbl>, MYS2016 <dbl>, EXP2015 <dbl>, EXP2016 <dbl>,
## # APM.SD15 <dbl>, APM.SMP15 <dbl>, APM.SMA15 <dbl>, APM.PT15 <dbl>,
## # APK.SD15 <dbl>, APK.SMP15 <dbl>, APK.SMA15 <dbl>, APK.PT15 <dbl>,
## # APS.USIA15 <dbl>, APS.USIA2 <dbl>, APS.USIA3 <dbl>, APS.USIA4 <dbl>Data Polygon
setwd("D:\\Spasial\\")
petajabar <- readOGR(dsn="petaJabar2", layer="Jabar2") #dsn:nama folder,layer:nama file shp## OGR data source with driver: ESRI Shapefile
## Source: "D:\Spasial\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields
Eksplorasi Data
Berikut adalah sebaran persentase penduduk miskin tahun 2016 di Jawa Barat. Terdapat 3 pembagian warna yaitu hijau, kuning, dan merah yang didasarkan pada persentase penduduk miskin. Semakin wilayah berwarna merah, maka semakin tinggi persentase penduduk miskin di wilayah tersebut.
colfunc <- colorRampPalette(c("green", "yellow", "red"))
color <- colfunc(16)
petajabar$miskin<-dataJabar$p.miskin16
spplot(petajabar, "miskin", col.regions=color, main="Peta Sebaran Persentase Penduduk Miskin 2016")
Indeks Moran
Indeks Moran digunakan untuk melihat ada tidaknya autokorelasi spasial.
Autokorelasi spasial merupakan teknik dalam analisis spasial untuk mengukur kemiripan nilai atribut dalam suatu ruang (jarak, waktu dan area). Jika terdapat pola sistematik dalam nilai atribut tersebut, maka terdapat autokorelasi spasial. Adanya autokorelasi spasial mengindikasikan bahwa nilai atribut pada area tertentu terkait oleh nilai atribut tersebut pada area lain yang letaknya saling berdekatan (bertetangga). Ketetanggaan tersebut diharapkan dapat mencerminkan derajat ketergantungan area (spasial) yang tinggi apabila dibandingkan dengan area lain yang letaknya terpisah jauh.
w <- poly2nb(petajabar)
ww <- nb2listw(w)
moran(dataJabar$p.miskin16, ww, n=length(ww$neighbours), S0=Szero(ww))## $I
## [1] 0.3932657
##
## $K
## [1] 2.403804##
## Moran I test under randomisation
##
## data: dataJabar$p.miskin16
## weights: ww
##
## Moran I statistic standard deviate = 3.0168, p-value = 0.001277
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3932657 -0.0400000 0.0206265H0 : Tidak terdapat autokorelasi spasial
H1 : Terdapat autokorelasi spasial
Statistik uji : p-value=0.001277
Keputusan : karena p-value < alpha (0.05), maka TOLAK H0
Kesimpulan : cukup bukti untuk mengatakan bahwa terdapat autokorelasi spasial
Terlihat pula dari nilai indeks moran yaitu 0.3932657. Nilai indeks moran berkisar antara -1 sampai dengan 1. Dimana jika nilai indeks moran positif, maka menunjukkan adanya autokorealasi positif.
Plot Moran
Plot Moran terbagi ke dalam 4 kuadran. Dalam plot terlihat kota Cirebon masuk ke dalam kuadran IV yaitu hotspot - high low. Artinya, Cirebon termasuk ke dalam daerah yang memiliki nilai persentase penduduk miskin yang tinggi, sedangkan wilayah sekitarnya memiliki persentase yang rendah.
Jika daerah masuk ke dalam kuadran low-low, artinya daerah tersebut memiliki persentase kemiskinan yang rendah dan daerah sekitarnya juga memiliki persentase kemiskinan yang rendah pula.
Model Regresi Klasik
Selanjutnya dibuat model regresi klasik penduduk miskin dengan variabel penjelas EYS2016.
reg.klasik <- lm(p.miskin16~EYS2016 , dataJabar)
err.regklasik <- residuals(reg.klasik)
summary(reg.klasik)##
## Call:
## lm(formula = p.miskin16 ~ EYS2016, data = dataJabar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3627 -1.9604 -0.3899 1.5238 8.1125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.5732 9.0973 4.350 0.000217 ***
## EYS2016 -2.3948 0.7209 -3.322 0.002854 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.788 on 24 degrees of freedom
## Multiple R-squared: 0.315, Adjusted R-squared: 0.2865
## F-statistic: 11.04 on 1 and 24 DF, p-value: 0.002854Uji Asumsi klasik
Akan dilakukan pengecekan asumsi kenormalan, kehomogenan ragam, dan kebebasan galat.
Uji Kenormalan galat
##
## Anderson-Darling normality test
##
## data: err.regklasik
## A = 0.39412, p-value = 0.3495H0 : galat model menyebar normal
H1 : galat model tidak menyebar normal
Statistik uji : p-value = 0.3495
Keputusan : karena p-value > alpha (0.05) maka Tidak tolak H0
Kesimpulan : Tidak cukup bukti untuk menyatakan galat model tidak menyebar normal
Q-Q Plot
qqnorm(err.regklasik,datax=T)
qqline(rnorm(length(err.regklasik),mean(err.regklasik),sd(err.regklasik)),datax=T, col="red")
Dari Q-Q plot juga terlihat bahwa galat berada di sekitar garis 45 derajat, sehingga makin terlihat bahwa asumsi galat menyebar normal terpenuhi.
Uji kehomogenan ragam
##
## studentized Breusch-Pagan test
##
## data: reg.klasik
## BP = 1.0644, df = 1, p-value = 0.3022H0 : ragam galat homogen
H1 : ragam galat tidak homogen
Statistik uji : p-value = 0.3022
Keputusan : karena p-value > alpha (0.05) maka Tidak tolak H0
Kesimpulan : Tidak cukup bukti untuk menyatakan ragam galat tidak homogen
Uji Kebebasan Galat
##
## Runs Test for Randomness
##
## data: err.regklasik
## runs = 13, m = 13, n = 13, p-value = 0.8358
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x)
## -0.3898635H0 : galat model saling bebas
H1 : galat model tidak saling bebas
Statistik uji : p-value = 0.8358
Keputusan : karena p-value > alpha (0.05) maka Tidak tolak H0
Kesimpulan : Tidak cukup bukti untuk menyatakan galat model tidak saling bebas.
Dari ketiga uji diatas, dapat disimpulkan bahwa model yang dibuat memenuhi seluruh asumsi klasik yaitu asumsi kenormalan, kehomogenan ragam, dan kebebasan galat.
Spatial Autocorrelation : Moran’s Index
Jika sebelumnya moran index digunakan untuk melihat autokorelasi pada respon penduduk miskin, maka selanjutnya akan dilihat moran index atau autokorelasi pada error dari regresi klasik yang telah dibuat.
w <- poly2nb(petajabar)
ww <- nb2listw(w)
moran.test(err.regklasik, ww, randomisation=T, alternative="greater")##
## Moran I test under randomisation
##
## data: err.regklasik
## weights: ww
##
## Moran I statistic standard deviate = 3.5604, p-value = 0.0001851
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.44799554 -0.04000000 0.01878561H0 : Tidak terdapat autokorelasi spasial
H1 : Terdapat autokorelasi spasial
Statistik uji : p-value=0.0001851
Keputusan : karena p-value < alpha (0.05), maka TOLAK H0
Kesimpulan : cukup bukti untuk mengatakan bahwa terdapat autokorelasi spasial
Terlihat pula dari nilai indeks moran yaitu 0.4479. Nilai indeks moran berkisar antara -1 sampai dengan 1. Dimana jika nilai indeks moran positif, maka menunjukkan adanya autokorealasi positif.
Dapat disimpulkan bahwa dari model regresi klasik yang dibuat, terdapat autokorelasi spasial, sehingga akan dilakukan uji efek spasial terhadap model tersebut.
Uji Efek Spasial
Berikut adalah gambaran umum tentang alur uji efek spasial
LM <- lm.LMtests(reg.klasik, nb2listw(w, style="W"),test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
summary(LM)## Lagrange multiplier diagnostics for spatial dependence
## data:
## model: lm(formula = p.miskin16 ~ EYS2016, data = dataJabar)
## weights: nb2listw(w, style = "W")
##
## statistic parameter p.value
## LMerr 8.375484 1 0.003803 **
## LMlag 6.642323 1 0.009958 **
## RLMerr 1.822499 1 0.177016
## RLMlag 0.089338 1 0.765020
## SARMA 8.464822 2 0.014517 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Uji diatas memberikan hasil dari berbagai tes model efek spasial. Terlihat dari ke-5 model yang diujikan, terdapat 3 model yang signifikan berdasarkan p-value nya yaitu LMerr, LMlag, dan SARMA.
Model Spasial
Model SAR
Spatial Autoregresive Model (SAR) merupakan model spasial yang terjadi akibat adanya pengaruh spasial pada variabel dependen. Apabila data yang diperoleh menghasilkan dependensi lag maka data dimodelkan dengan SAR.
##
## Call:
## lagsarlm(formula = p.miskin16 ~ EYS2016, data = dataJabar, listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.17670 -0.93185 -0.11318 0.91353 7.69131
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 28.22898 7.66250 3.6840 0.0002296
## EYS2016 -1.94997 0.57325 -3.4016 0.0006700
##
## Rho: 0.59078, LR test value: 7.9343, p-value: 0.0048507
## Asymptotic standard error: 0.1559
## z-value: 3.7894, p-value: 0.00015101
## Wald statistic: 14.36, p-value: 0.00015101
##
## Log likelihood: -58.54132 for lag model
## ML residual variance (sigma squared): 4.7513, (sigma: 2.1798)
## Number of observations: 26
## Number of parameters estimated: 4
## AIC: 125.08, (AIC for lm: 131.02)
## LM test for residual autocorrelation
## test value: 0.036687, p-value: 0.8481Dari hasil Model SAR, dapat dilihat nilai p-value (0.0048507) yang signifikan, sehingga akan dilakukan pengecekan asumsi dari model SAR.
Pengecekan Asumsi Model SAR
Akan dilakukan pengecekan asumsi kenormalan, kehomogenan ragam, dan autokorelasi (dengan indeks moran) pada model SAR.
##
## Anderson-Darling normality test
##
## data: err.sar
## A = 0.70904, p-value = 0.05644##
## studentized Breusch-Pagan test
##
## data:
## BP = 0.97078, df = 1, p-value = 0.3245##
## Moran I test under randomisation
##
## data: err.sar
## weights: ww
##
## Moran I statistic standard deviate = 0.20765, p-value = 0.4178
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.01331736 -0.04000000 0.01651197Dari uji asumsi model SAR terlihat bahwa dari uji kenormalan, kehomogenan ragam, dan autokorelasi, p-value > alpha (0.05) tidak tolak H0. Artinya, galat menyebar normal, ragam galat homogen, dan tidak terdapat autokorelasi spasial pada model SAR sehingga dapat dikatakan bahwa semua asumsi klasik terpenuhi.
Model SEM
Spatial Error Model (SEM) merupakan model spasial yang terjadi akibat adanya pengaruh spasial pada error. Apabila data yang diperoleh menghasilkan dependensi error maka data dimodelkan dengan SEM.
##
## Call:errorsarlm(formula = p.miskin16 ~ EYS2016, data = dataJabar,
## listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.24699 -1.36122 -0.13809 1.15579 7.03019
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 36.88515 7.80246 4.7274 2.274e-06
## EYS2016 -2.17498 0.60042 -3.6224 0.0002919
##
## Lambda: 0.61793, LR test value: 8.7676, p-value: 0.0030663
## Asymptotic standard error: 0.1576
## z-value: 3.9208, p-value: 8.8267e-05
## Wald statistic: 15.372, p-value: 8.8267e-05
##
## Log likelihood: -58.12466 for error model
## ML residual variance (sigma squared): 4.5459, (sigma: 2.1321)
## Number of observations: 26
## Number of parameters estimated: 4
## AIC: 124.25, (AIC for lm: 131.02)Dari hasil Model SEM, dapat dilihat nilai p-value (0.0030663) yang signifikan, sehingga akan dilakukan pengecekan asumsi dari model SEM.
Pengecekan Asumsi Model SEM
Akan dilakukan pengecekan asumsi kenormalan, kehomogenan ragam, dan autokorelasi (dengan indeks moran) pada model SEM.
##
## Anderson-Darling normality test
##
## data: err.sem
## A = 0.63011, p-value = 0.08968##
## studentized Breusch-Pagan test
##
## data:
## BP = 0.38918, df = 1, p-value = 0.5327##
## Moran I test under randomisation
##
## data: err.sem
## weights: ww
##
## Moran I statistic standard deviate = -0.096995, p-value = 0.5386
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.05288671 -0.04000000 0.01765155Dari uji asumsi model SEM terlihat bahwa dari uji kenormalan, kehomogenan ragam, dan autokorelasi, p-value > alpha (0.05) tidak tolak H0. Artinya, galat menyebar normal, ragam galat homogen, dan tidak terdapat autokorelasi spasial pada model SEM sehingga dapat dikatakan bahwa semua asumsi klasik terpenuhi.
Model SARMA/GSM
Jika data menghasilkan dependensi lag dan dependensi error maka data dimodelkan dengan Spatial Autoregressive Moving Average (SARMA).
##
## Call:
## sacsarlm(formula = p.miskin16 ~ EYS2016, data = dataJabar, listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0817 -1.2815 -0.1627 1.1317 6.7128
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 38.02312 8.19667 4.6388 3.504e-06
## EYS2016 -2.15311 0.59378 -3.6261 0.0002877
##
## Rho: -0.14689
## Asymptotic standard error: 0.38983
## z-value: -0.3768, p-value: 0.70633
## Lambda: 0.69352
## Asymptotic standard error: 0.23014
## z-value: 3.0135, p-value: 0.0025829
##
## LR test value: 8.8931, p-value: 0.011719
##
## Log likelihood: -58.06189 for sac model
## ML residual variance (sigma squared): 4.3254, (sigma: 2.0797)
## Number of observations: 26
## Number of parameters estimated: 5
## AIC: 126.12, (AIC for lm: 131.02)Dari hasil Model SARMA/GSM, dapat dilihat nilai p-value yang signifikan, sehingga akan dilakukan pengecekan asumsi dari model SARMA/GSM.
Pengecekan Asumsi Model SARMA/GSM
Akan dilakukan pengecekan asumsi kenormalan, kehomogenan ragam, dan autokorelasi (dengan indeks moran) pada model SARMA/GSM.
##
## Anderson-Darling normality test
##
## data: err.gsm
## A = 0.55307, p-value = 0.1388##
## studentized Breusch-Pagan test
##
## data:
## BP = 0.31097, df = 1, p-value = 0.5771##
## Moran I test under randomisation
##
## data: err.gsm
## weights: ww
##
## Moran I statistic standard deviate = -0.0061095, p-value = 0.5024
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04081914 -0.04000000 0.01797655Dari uji asumsi model SARMA/GSM terlihat bahwa dari uji kenormalan, kehomogenan ragam, dan autokorelasi, p-value > alpha (0.05) tidak tolak H0. Artinya, galat menyebar normal, ragam galat homogen, dan tidak terdapat autokorelasi spasial pada model SARMA/GSM sehingga dapat dikatakan bahwa semua asumsi klasik terpenuhi.
Goodness of Fit
Tabel diatas membandingkan model regresi klasik dan ketiga model spasial yang sudah dibuat sebelumnya yaitu SAR, SEM, dan SARMA.
Variabel yang dibandingkan adalah nilai AIC, Rho, Lambda, dan asumsi (kenormalan, kehomogenan ragam, dan kebebasan galat).
- Asumsi
Seperti sudah diuji sebelumnya, model regresi klasik dan ketiga model spasial semuanya memenuhi ketiga asumsi.
- Signifikansi Rho dan Lambda
Model SAR hanya memiliki Rho, model SEM hanya memiliki lambda, dan Model SARMA memiliki Rho dan Lambda. Model SARMA baik jika Rho dan Lambda nya keduanya signifikan (Jika Rho tidak signifikan akan sama dengan SEM, jika Lambda tidak signifikan akan sama dengan SAR, jika keduanya tidak signifikan sama dengan Model Regresi Klasik).
Terlihat dari model SARMA, Rho tidak signifikan, sehingga model SARMA sama saja dengan model SEM. Sehingga hanya tersisa 3 model saja yaitu Model Regresi Klasik, Model SAR, dan Model SEM.
- AIC
Model terbaik adalah model yang memiliki nilai AIC terkecil. Sehingga jika dilihat dari nilai AIC ketiga kandidat model, Model SEM yang terbaik.
Heterogenity Spasial
Contoh kedua, akan dilakukan analisis regresi spasial tentang PDRB di Jawa Timur.
Input Data
setwd("D:\\Spasial\\")
datajatim <- read_xlsx("data gwr jatim.xlsx")
jatim <- readOGR(dsn="Poligon Jatim", layer="poly3500")## OGR data source with driver: ESRI Shapefile
## Source: "D:\Spasial\Poligon Jatim", layer: "poly3500"
## with 38 features
## It has 8 fields## # A tibble: 38 x 6
## Kabupaten Long Lat Y x1 x2
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Pacitan 111. -8.1 9.02 350 0.97
## 2 Kab. Ponorogo 112. -7.97 11.7 467 3.68
## 3 Kab. Trenggalek 112. -8.14 10.5 395 2.46
## 4 Kab. Tulungagung 112. -8.08 22.3 526 3.95
## 5 Kab. Blitar 112. -8.14 20.9 581 2.79
## 6 Kab. Kediri 112. -7.8 24.0 761 5.02
## 7 Kab. Malang 113. -8.09 55.3 1228 4.95
## 8 Kab. Lumajang 113. -8.12 18.7 518 2.6
## 9 Kab. Jember 114. -8.26 44.2 1117 4.77
## 10 Kab. Banyuwangi 114. -8.33 44.5 871 2.55
## # ... with 28 more rows
Eksplorasi Data
Akan dilakukan eksplorasi data yaitu data PDRB di Provinsi Jawa Timur. Semakin pekat warna, maka nilai PDRB di daerah tersebut semakin tinggi.
jatim@data<-cbind(pdrb=datajatim[,2],jatim@data)
my.palette.1<-brewer.pal(n=9,name="YlOrRd")
spplot(jatim[,1], cuts=8, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="PDRB jatim")
Model Regresi Klasik
Akan dilakukan permodelan antara PDRB (Y) dengan 2 variabel penjelas yaitu X1 dan X2
##
## Call:
## lm(formula = Y ~ x1 + x2, data = datajatim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.616 -13.475 -0.679 10.545 164.572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -72.45320 18.54899 -3.906 0.000409 ***
## x1 0.10788 0.01952 5.528 3.24e-06 ***
## x2 12.10245 3.51225 3.446 0.001497 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.82 on 35 degrees of freedom
## Multiple R-squared: 0.5671, Adjusted R-squared: 0.5424
## F-statistic: 22.93 on 2 and 35 DF, p-value: 4.332e-07Uji Asumsi Model Regresi Klasik
Uji asumsi yang akan dilakukan adalah uji multikolinieritas (karena terdapat lebih dari 1 variabel penjelas), uji kenormalan, dan uji kehomogenan ragam.
## x1 x2
## 1.006723 1.006723##
## Anderson-Darling normality test
##
## data: err.ols
## A = 1.9608, p-value = 4.258e-05##
## studentized Breusch-Pagan test
##
## data: reglin
## BP = 13.836, df = 2, p-value = 0.0009896Uji Multikolinieritas
Jika nilai VIF > 5, maka terjadi multikolinieritas. Terlihat dalam hasil uji VIF, nilai VIF x1 dan x2 adalah 1.006723 sehingga dapat dikatakan bahwa tidak ada multikolinieritas
Uji kenormalan galat
H0 : galat menyebar normal
H1 : galat tidak menyebar normal
Statistik Uji : p-value=4.258e-05
Keputusan : karena p-value < alpha (0.05) maka TOLAK H0
Kesimpulan : Cukup bukti untuk menyatakan bahwa galat tidak menyebar normal
Uji Kehomogenan Ragam
H0 : Ragam galat homogen
H1 : Ragam galat tidak homogen
Statistik uji : p-value = 0.0009896
Keputusan : karena p-value < alpha (0.05) maka TOLAK H0
Kesimpulan : Cukup bukti untuk menyatakan bahwa ragam galat tidak homogen.
Uji Efek Spasial
Dari uji Breusch-Pagan yang telah dilakukan, ditemukan adanya ragam galat yang tidak homogen. Sehingga akan dilakukan mode GWR.
Geographically Weighted Regression (GWR)
Model GWR adalah model regresi yang dikembangkan untuk memodelkan data dengan variabel respon yang bersifat kontinu dan mempertimbangkan aspek spasial atau lokasi.
Pendekatan yang dilakukan dalam GWR adalah pendekatan titik. Setiap nilai parameter ditaksir pada setiap titik lokasi pengamatan, sehingga setiap titik lokasi pengamatan mempunyai nilai parameter yang berbeda-beda.
GWR - Jarak Spasial
Matriks pembobot pada GWR merupakan matriks pembobot berbasis pada kedekatan titik lokasi pengamatan ke-i dengan titik lokasi pengamatan lainnya. Pengamatan terdekat ke titik lokasi pengamatan ke-i umumnya diasumsikan memiliki pengaruh paling besar terhadap penaksiran parameter di titik lokasi pengamatan ke-i. Oleh karena itu, matriks pembobot W (ui,vj) akan semakin besar seperti jarak yang semakin dekat.
cord<-cbind(datajatim$Long,datajatim$Lat)
DM <-gw.dist(dp.locat=cord)
datajatim.sp<-datajatim
coordinates(datajatim.sp)<-~Long+Lat
class(datajatim.sp)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"Berikut adalah datajatim.sp
## class : SpatialPointsDataFrame
## features : 38
## extent : 111.16, 114.26, -8.33, -4.55 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 4
## names : Kabupaten, Y, x1, x2
## min values : Kab. Bangkalan, 3.86, 64, 0.97
## max values : Kota Surabaya, 324.22, 1365, 8.46GWR - Bandwitch
Bandwidth dapat dianalogikan sebagai radius (b) suatu lingkaran, sehingga sebuah titik lokasi pengamatan yang berada dalam radius (b) lingkaran masih dianggap berpengaruh dalam membentuk parameter di titik lokasi pengamatan ke-i. Pemilihan bandwidth optimum dalam GWR merupakan hal yang penting karena akan mempengaruhi ketepatan model terhadap data.
Nilai bandwidth yang sangat kecil akan mengakibatkan penaksiran parameter di lokasi pengamatan ke-i semakin bergantung pada titik lokasi pengamatan lain yang memiliki jarak terdekat dengan lokasi pengamatan ke-i, sehingga varians yang dihasilkan akan semakin besar. Sebaliknya, jika nilai bandwidth sangat besar makan akan mengakibatkan bias semakin besar, sehingga model yang diperoleh terlalu halus (Dwinata, 2012).
Metode yang digunakan untuk menentukan bandwidth optimum salah satunya adalah metode validasi silang atau cross validation (CV). Nilai bandwidth optimum diperoleh ketika CV minimum (Fotheringham, 2002).
Perhitungan bandwith akan dipilih nilai CV terkecil, sehingga akan dilakukan iterasi. Ketika telah ditemukan CV terkecil, maka bandwith tersebut yang akan digunakan untuk menghitung pembobot.
## Fixed bandwidth: 2.568552 CV score: 77431.75
## Fixed bandwidth: 1.58777 CV score: 76850.32
## Fixed bandwidth: 0.9816131 CV score: 74446.36
## Fixed bandwidth: 0.6069877 CV score: 67449.85
## Fixed bandwidth: 0.3754565 CV score: 57557.27
## Fixed bandwidth: 0.2323623 CV score: 58990.41
## Fixed bandwidth: 0.4638935 CV score: 61547.68
## Fixed bandwidth: 0.3207994 CV score: 55973.74
## Fixed bandwidth: 0.2870194 CV score: 56026.58
## Fixed bandwidth: 0.3416765 CV score: 56403.23
## Fixed bandwidth: 0.3078966 CV score: 55867.61
## Fixed bandwidth: 0.2999222 CV score: 55876.05
## Fixed bandwidth: 0.312825 CV score: 55891.64
## Fixed bandwidth: 0.3048507 CV score: 55863.7
## Fixed bandwidth: 0.3029682 CV score: 55865.64
## Fixed bandwidth: 0.3060141 CV score: 55864.17
## Fixed bandwidth: 0.3041316 CV score: 55864.04
## Fixed bandwidth: 0.305295 CV score: 55863.73
## Fixed bandwidth: 0.304576 CV score: 55863.77
## Fixed bandwidth: 0.3050204 CV score: 55863.69
## Fixed bandwidth: 0.3051253 CV score: 55863.7
## Fixed bandwidth: 0.3049556 CV score: 55863.69
## Fixed bandwidth: 0.3050605 CV score: 55863.69
## Fixed bandwidth: 0.3049956 CV score: 55863.69GWR - Model
GWR menghasilkan penduga parameter yang berbeda untuk tiap wilayah. Data yang ada dalam datajatim.sp adalah sebanyak 38 wilayah, maka akan ada 38 model yang berbeda dengan nilai penduga parameter (beta) yang berbeda pula.
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2021-03-22 23:14:38
## Call:
## gwr.basic(formula = Y ~ x1 + x2, data = datajatim.sp, bw = bw.fg,
## kernel = "gaussian", dMat = DM)
##
## Dependent (y) variable: Y
## Independent variables: x1 x2
## Number of data points: 38
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.616 -13.475 -0.679 10.545 164.572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -72.45320 18.54899 -3.906 0.000409 ***
## x1 0.10788 0.01952 5.528 3.24e-06 ***
## x2 12.10245 3.51225 3.446 0.001497 **
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 36.82 on 35 degrees of freedom
## Multiple R-squared: 0.5671
## Adjusted R-squared: 0.5424
## F-statistic: 22.93 on 2 and 35 DF, p-value: 4.332e-07
## ***Extra Diagnostic information
## Residual sum of squares: 47462.13
## Sigma(hat): 36.30967
## AIC: 386.7832
## AICc: 387.9953
## BIC: 369.8839
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: gaussian
## Fixed bandwidth: 0.3049956
## Regression points: the same locations as observations are used.
## Distance metric: A distance matrix is specified for this model calibration.
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept -185.220437 -80.448573 -34.683457 -13.604522 4.2308
## x1 0.010088 0.033628 0.062080 0.122049 0.2896
## x2 -3.672824 3.685552 8.853189 13.036107 27.4911
## ************************Diagnostic information*************************
## Number of data points: 38
## Effective number of parameters (2trace(S) - trace(S'S)): 24.9044
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 13.0956
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 430.7659
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 345.5029
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 361.6579
## Residual sum of squares: 11516.94
## R-square value: 0.894956
## Adjusted R-square value: 0.6786741
##
## ***********************************************************************
## Program stops at: 2021-03-22 23:14:38Tampilan model GWR dalam boxplot
hasil=as.data.frame(gwr1$SDF)
par(mfrow=c(1,3))
boxplot(hasil[,1],xlab="Intercept")
boxplot(hasil[,2],xlab="X1")
boxplot(hasil[,3],xlab="X2")
GWR - Uji t
Akan dilakukan ji keberartian model regresi linear berganda dengan menggunakan uji t. Akan diuji tiap model sehingga akan ada 38 hasil dari uji t.
## Intercept_t x1_t x2_t
## 1 0.07500666 0.1013885 0.130699595
## 2 -0.14929048 0.2157656 0.812943235
## 3 -0.24382033 0.2121404 1.100402805
## 4 -0.77809332 0.5073455 1.817783449
## 5 -1.03150976 1.1678547 1.901986651
## 6 -1.28452274 1.1962923 2.077186689
## 7 -1.30932769 1.8582087 1.659964084
## 8 -0.94327543 1.9873844 1.284157922
## 9 -0.12244755 1.1741318 0.071774007
## 10 -0.04982797 1.0021626 -0.222302493
## 11 -0.15878387 1.0351021 0.006655903
## 12 -0.21904872 1.0101595 0.038605042
## 13 -1.48059172 2.2312068 1.622192589
## 14 -2.90120720 5.1851252 2.495743406
## 15 -3.57748166 7.1541738 2.833026493
## 16 -3.57084305 6.6247468 3.180363650
## 17 -3.10016529 4.8917619 2.976466370
## 18 -1.05069180 1.2991928 1.681710514
## 19 -0.43015507 0.6148005 0.919913579
## 20 -0.09618664 0.3231631 0.333303056
## 21 -0.23681017 0.6419804 0.325916469
## 22 -0.87810784 1.7396100 0.928864528
## 23 -1.33830841 3.7635588 0.802950342
## 24 -4.26880431 6.3071611 2.974985041
## 25 -4.09884082 6.1570722 2.909876338
## 26 -4.36032747 4.4683034 0.711817195
## 27 -0.14320952 0.2399921 0.109873140
## 28 -0.11969202 0.2650994 -0.015039046
## 29 -0.36947655 0.4607246 0.276731865
## 30 -1.11517121 0.9381891 1.981325867
## 31 -1.06567257 1.0890099 1.983133030
## 32 -1.69965880 2.5145669 1.951437978
## 33 -1.95150401 3.1030348 1.925933418
## 34 -3.19647847 5.8304793 2.534254030
## 35 -3.81335593 6.6067751 3.293334465
## 36 -0.28338527 0.4823571 0.660346626
## 37 -3.85246507 7.0155083 2.869742511
## 38 -2.13110280 3.3716669 2.424206100Referensi
Djuraidah, Anik. (2021). STA500 Metode Penelitan Kuantitatif: Spasial. Retrieved from https://newlms.ipb.ac.id/
Hujjah A, Al. (2021). Spatial Regression. Retrieved from https://newlms.ipb.ac.id/
https://www.mobilestatistik.com/regresi-spasial-gwr/
Nurul Muthiah, Raupong, Anisa. ESTIMASI PARAMETER REGRESI SPATIAL AUTOREGRESSIVE MODEL. Retrieved from https://core.ac.uk/download/pdf/77620924.pdf