Pemodelan Regresi Spasial (Spatial Regression Model) dengan R
 |
| Regresi Spasial dengan R |
Pada unggahan sebelumnya, kita telah belajar bersama mengenai 2 jenis teknik pengambilan sampel, yaitu Simple Random Sampling dan Systematic Sampling. Kali ini, kita akan kembali mengulas mengenai pemodelan statistik, yaitu regresi spasial.
Regresi spasial pada dasarnya merupakan model hasil pengembangan regresi linier sederhana dan berganda atau yang biasa disebut sebagai regresi linier klasik. Istilah spasial sendiri mengacu pada pengaruh spasial, wilayah, atau tempat terhadap variabel yang dianalisis. Sebagai sebuah kelaziman, efek spasial sendiri demikian mudah kita temukan dalam keseharian, antara satu wilayah dengan wilayah lainnya. Dalam kondisi tertentu, suatu variabel bebas (independen) suatu wilayah amatan berkaitan erat dengan variabel bebas yang sama di wilayah lain, terlebih untuk wilayah-wilayah yang berdekatan atau diistilahkan dengan wilayah bertetangga (neighbour). Adanya keterkaitan spasial inilah yang mengakibatkan pendugaan (estimasi) menjadi tidak tepat sekaligus menengarai adanya pelanggaran asumsi keacakan. Untuk mengatasi hal tersebut, regresi linier perlu mengakomodir keterkaitan spasial antar wilayah yang kemudian disepakati sebagai regresi spasial.
 |
| Gambar 1. Alur kapan menggunakan regresi spasial |
Berdasarkan gambar 1, pertama kali yang perlu dilakukan untuk memastikan apakah dalam model regresi linier mengandung keterkaitan spasial adalah dengan mengecek uji asumsi homoskedastisitas dan uji asumsi non-autokorelasinya. Bila yang terlanggar adalah uji homoskedastisitas, maka arah pemodelan regresi linier berikutnya adalah model Geographically Weighted Regression (GWR). Sebaliknya, bila yang terlanggar adalah asumsi non-autokorelasinya maka model regresi spasial yang relevan adalah regresi spasial, SAR, SEM, GSM, Spasial Durbin dan variannya.
Apa saja tahapan pemodelan regresi spasial yang terlanggar asumsi non-autokorelasi (dependency spatial) dengan R? Berikut tahapannya:
1. Melakukan instalasi sejumlah package yang diperlukan sekaligus aktivasinya dengan fungsi library();
2. Melakukan import data;
3. Eksplorasi data pada peta. Pada bagian ini package visualisasi spasial sekaligus caranya dapat disimak pada tautan [1], [2], atau [3];
4. Membuat model regresi linier klasik dengan teknik estimasi Ordinary Least Square (OLS);
5. Melakukan uji asumsi regresi klasik. Pengujian asumsi dapat disimak lebih lanjut pada tautan [1],[2],[3],[4];
6. Apakah terdapat pelanggaran asumsi non-autokorelasi? Bila Ya, maka lanjut;
7. Mengecek nilai sekaligus uji signifikansi Indeks Moran serta Plotting Indeks Moran; Pada bagian ini, hasil pada praktik (sebagai contoh), menunjukkan bahwa Sidoarjo masuk pada kuadran 4 atau hot spot (high - low). Artinya, Sidoarjo masuk sebagai wilayah dengan Indeks Pembangunan Manusia (IPM) yang tinggi, sementara itu wilayah sekitarnya (tetangganya) justru memiliki IPM yang rendah; Pada bagian ini pula, hipotesis pengujian signifikansi Indeks Moran adalah sebagai berikut:
H0: Tidak terdapat autokorelasi spasial
H1: Terdapat autokorelasi spasial
Keputusan dari pengujian ini dapat dilihat berdasarkan nilai p-value. Apabila nilai p-value > alpha, maka dikatakan cukup bukti bahwa terdapat autokorelasi spasial. Pada tahapan ini pula, bila uji Indeks Moran signifikan, maka perlu dilakukan pengujian autokorelasi spasial lokalnya atau yang diistilahkan dengan Local Indicator of Spatial Assosiation (LISA). Kalau diumpamakan uji Indeks Moran mirip dengan uji F atau uji simultannya, sedangkan uji LISA adalah uji t atau uji parsialnya;
8. Melakukan pengujian efek spasial sekaligus memperoleh rekomendasi model spasial yang relevan digunakan. Adapun alur pengujian dan rekomendasi model spasial tertera pada bagan berikut:
 |
| Gambar 2. Alur pemodelan regresi spasial |
9. Membentuk model spasial;
10. Melakukan uji asumsi regresi spasial, mulai dari kenormalan residual, uji homoskedastisitas model spasial, serta uji asumsi non-autokorelasi model spasial; Pada pengujian ini, keputusan uji sama persis dengan pengujian asumsi regresi linier, terima H0 ketika nilai p-value > alpha penelitian;
11. Melakukan pemilihan model (bila model yang relevan lebih dari 1) terbaik. Dasar pemilihan model terbariknya adalah nilai AIC terkecil, Adjuster R Square terbesar, serta signifikansi nilai Rho dan Lambda;
12. Interpretasi model spasial.
Setelah mengetahui step-step pemodelan, berikut adalah langkah praktik pemodelan regresi spasial (dependency) menggunakan R. Adapun data yang digunakan dalam praktikum kali ini bersumber dari Badan Pusat Statistik Provinsi Jawa Timur tahun 2021. Untuk mempersiapkan datanya dapat mengunduhnya pada tautan [1]. Setelah datanya siap, praktik pemodelannya adalah sebagai berikut:
#Import Data
library(readxl)
modku <- read_excel("modku.xlsx")
modku
## # A tibble: 38 x 6
## Kako Kode IPM APK APM PDRBPK
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pacitan 3501 68.6 97.5 80.3 18855
## 2 Ponorogo 3502 71.1 96.4 84.6 15295
## 3 Trenggalek 3503 70.1 94.9 83.2 17634
## 4 Tulungagung 3504 73.2 91.8 82.1 24978
## 5 Lumajang 3508 66.1 89.3 74.1 20072
## 6 Bondowoso 3511 66.6 94.6 77.2 17882
## 7 Pasuruan 3514 68.9 94.4 75.4 66776
## 8 Jombang 3517 73.4 101. 86.8 21535
## 9 Nganjuk 3518 72.0 95.9 84.7 16798
## 10 Madiun 3519 71.9 97.6 86.2 17826
## # ... with 28 more rows
#Visualisasi peta jatim
library(rgdal)
ptjatim=readOGR(dsn=path.expand("E:/JOKO ADE/R/Modul UNUGIRI"),layer="jatim")
## OGR data source with driver: ESRI Shapefile
## Source: "E:\JOKO ADE\R\Modul UNUGIRI", layer: "jatim"
## with 38 features
## It has 17 fields
## Integer64 fields read as strings: KODE
#Menyatukan data dengan peta
colfunc <- colorRampPalette(c("green", "yellow", "red"))
color <- colfunc(16)
ptjatim$ipm <- modku$IPM
spplot(ptjatim, "ipm", col.regions=color, main="Peta Sebaran Indeks Pembangunan Manusia Kabkot Jawa Timur 2021")
#Mendapatkan nilai Indeks Moran
library(spdep)
w <- poly2nb(ptjatim)
ww <- nb2listw(w)
moran(modku$IPM, ww, n=length(ww$neighbours), S0=Szero(ww))
## $I
## [1] 0.4530106
##
## $K
## [1] 2.363749
#Uji signifikansi Indeks Moran
moran.test(modku$IPM, ww, alternative="greater")
##
## Moran I test under randomisation
##
## data: modku$IPM
## weights: ww
##
## Moran I statistic standard deviate = 3.6748, p-value = 0.000119
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.45301056 -0.02702703 0.01706421
#Plot Indeks Moran
moran.plot(modku$IPM, ww, labels= ptjatim$NAMA_KAB)
#Uji autokorelasi spasial lokal dengan LISA
local_M=localmoran(modku$IPM,ww)
local_M
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 0 0.243662199 -1.445514e-02 0.2631587129 0.50316254 0.61485000
## 1 0.027766405 -1.469038e-03 0.0079999269 0.32686368 0.74377099
## 2 0.112556751 -5.072316e-03 0.0603721383 0.47873631 0.63212623
## 3 -0.032770056 -9.251693e-04 0.0062442506 -0.40299490 0.68695198
## 4 0.972712174 -4.099002e-02 0.3423231992 1.73257571 0.08317111
## 5 0.911403079 -3.435685e-02 0.2889116303 1.75953697 0.07848635
## 6 -0.076378693 -1.174803e-02 0.0525214994 -0.28201359 0.77793309
## 7 -0.011125537 -1.622820e-03 0.0088360194 -0.10109262 0.91947694
## 8 0.002829727 -7.049547e-05 0.0003844343 0.14791773 0.88240770
## 9 -0.013149269 -1.289690e-04 0.0007032681 -0.49097632 0.62344320
## 10 0.121907237 -4.007999e-03 0.0347631022 0.67533463 0.49946319
## 11 0.016697732 -1.519898e-03 0.0181548496 0.13520587 0.89244910
## 12 0.052421047 -7.513337e-03 0.0406676532 0.29720184 0.76631242
## 13 0.115492888 -1.189107e-02 0.2170425563 0.27342756 0.78452457
## 14 0.010556413 -8.661150e-04 0.0058460202 0.14939363 0.88124304
## 15 2.967518116 -6.692855e-02 2.3730665895 1.96981095 0.04886004
## 16 1.703483969 -3.671267e-02 0.6532674243 2.15304360 0.03131525
## 17 0.085418593 -4.396524e-02 1.5972272635 0.10237564 0.91845851
## 18 -0.317782657 -4.936303e-02 1.7832000932 -0.20100841 0.84069200
## 19 -0.638541292 -1.042179e-01 3.5475474892 -0.28368743 0.77664993
## 20 -0.342600896 -2.227047e-03 0.0844393219 -1.17134227 0.24146123
## 21 -0.447798725 -1.246801e-02 0.4678772290 -0.63643397 0.52449360
## 22 0.478058850 -4.165160e-02 1.5168361354 0.42198024 0.67303944
## 23 0.285301189 -8.811584e-02 1.4842696949 0.30650511 0.75922009
## 24 2.563344240 -1.100307e-01 1.8088735168 1.98772262 0.04684238
## 25 -0.162083701 -1.778731e-02 0.2090039807 -0.31562994 0.75228341
## 26 -0.075143441 -1.494360e-03 0.0129939360 -0.64609588 0.51821726
## 27 0.012212395 -1.212249e-04 0.0006610447 0.47970600 0.63143646
## 28 0.194139176 -4.007999e-03 0.0152746964 1.60325126 0.10887919
## 29 0.169019553 -7.729810e-04 0.0092399946 1.76637427 0.07733308
## 30 0.912442386 -1.976986e-02 0.1687589347 2.26924516 0.02325342
## 31 0.913117104 -2.603212e-02 0.2207949897 1.99866495 0.04564462
## 32 -0.005512932 -2.857803e-03 0.0093580639 -0.02744685 0.97810333
## 33 0.840092156 -3.849853e-02 0.1674551046 2.14702647 0.03179117
## 34 2.582678039 -9.611073e-02 1.6047458326 2.11463464 0.03446110
## 35 1.095107986 -7.678875e-02 0.6173532409 1.49149732 0.13583097
## 36 0.738219152 -2.137863e-02 0.2502841330 1.51833299 0.12893047
## 37 1.209129976 -2.908885e-02 1.0732220893 1.19523307 0.23199598
## attr(,"call")
## localmoran(x = modku$IPM, listw = ww)
## attr(,"class")
## [1] "localmoran" "matrix" "array"
## attr(,"quadr")
## mean median pysal
## 1 Low-Low Low-Low Low-Low
## 2 Low-High Low-High Low-Low
## 3 Low-Low Low-Low Low-Low
## 4 High-High High-Low High-Low
## 5 Low-Low Low-Low Low-Low
## 6 Low-Low Low-Low Low-Low
## 7 Low-High Low-High Low-High
## 8 High-High High-High High-Low
## 9 Low-High High-High Low-Low
## 10 Low-High High-High Low-High
## 11 High-High High-High High-High
## 12 Low-High Low-High Low-Low
## 13 Low-High Low-High Low-Low
## 14 Low-High Low-Low Low-Low
## 15 High-High High-High High-High
## 16 Low-Low Low-Low Low-Low
## 17 Low-Low Low-Low Low-Low
## 18 High-High High-High High-High
## 19 High-High High-Low High-Low
## 20 High-Low High-Low High-Low
## 21 High-Low High-Low High-Low
## 22 High-Low High-Low High-Low
## 23 High-High High-High High-High
## 24 High-High High-High High-High
## 25 High-High High-High High-High
## 26 High-High High-Low High-Low
## 27 Low-High Low-High Low-High
## 28 High-High High-High High-High
## 29 High-High High-High High-High
## 30 Low-Low Low-Low Low-Low
## 31 High-High High-High High-High
## 32 Low-Low Low-Low Low-Low
## 33 Low-High Low-High Low-High
## 34 Low-Low Low-Low Low-Low
## 35 Low-Low Low-Low Low-Low
## 36 High-High High-High High-High
## 37 Low-Low Low-Low Low-Low
## 38 Low-Low Low-Low Low-Low
#Regresi OLS
regOLS <- lm(IPM~APK+APM+PDRBPK)
summary(regOLS)
##
## Call:
## lm(formula = IPM ~ APK + APM + PDRBPK)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1765 -1.3896 -0.6129 1.3433 5.9941
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.414e+00 9.826e+00 0.347 0.73042
## APK 4.242e-01 1.587e-01 2.672 0.01148 *
## APM 3.242e-01 1.481e-01 2.189 0.03552 *
## PDRBPK 3.138e-05 1.001e-05 3.134 0.00355 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.914 on 34 degrees of freedom
## Multiple R-squared: 0.6958, Adjusted R-squared: 0.669
## F-statistic: 25.93 on 3 and 34 DF, p-value: 6.547e-09
#Uji asumsi kenormalan residual
library(lmtest)
library(nortest)
ks.test(regOLS$residuals, ecdf(regOLS$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: regOLS$residuals
## D = 0.026316, p-value = 1
## alternative hypothesis: two-sided
#Uji asumsi homoskedastisitas residual
bptest(regOLS)
##
## studentized Breusch-Pagan test
##
## data: regOLS
## BP = 4.8658, df = 3, p-value = 0.1819
#Uji asumsi non-multikolinearitas variabel bebas
library(olsrr)
ols_vif_tol(regOLS)
## Variables Tolerance VIF
## 1 APK 0.3948246 2.532770
## 2 APM 0.3832948 2.608958
## 3 PDRBPK 0.9178622 1.089488
#Uji asumsi non-autokorelasi residual
dwtest(regOLS)
##
## Durbin-Watson test
##
## data: regOLS
## DW = 1.2098, p-value = 0.004027
## alternative hypothesis: true autocorrelation is greater than 0
#Uji efek spasial dan pemilihan model terbaik
LM <- lm.LMtests(regOLS, ww,test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
summary(LM)
## Lagrange multiplier diagnostics for spatial dependence
## data:
## model: lm(formula = IPM ~ APK + APM + PDRBPK)
## weights: ww
##
## statistic parameter p.value
## LMerr 0.02039 1 0.88645
## LMlag 2.96290 1 0.08520 .
## RLMerr 3.50454 1 0.06120 .
## RLMlag 6.44705 1 0.01111 *
## SARMA 6.46744 2 0.03941 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Regresi spasial
library(spatialreg)
#Model lag
SLX=lmSLX(regOLS,data=ptjatim,ww)
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
## -4.6012 -1.3976 -0.5537 0.9471 5.5108
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.553e+01 1.697e+01 -1.504 0.14267
## APK 2.982e-01 1.462e-01 2.040 0.04995 *
## APM 3.519e-01 1.495e-01 2.354 0.02506 *
## PDRBPK 2.611e-05 8.852e-06 2.949 0.00601 **
## lag.APK 5.695e-01 2.432e-01 2.342 0.02579 *
## lag.APM -2.057e-01 1.705e-01 -1.206 0.23675
## lag.PDRBPK 3.845e-05 2.693e-05 1.428 0.16342
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.534 on 31 degrees of freedom
## Multiple R-squared: 0.7903, Adjusted R-squared: 0.7497
## F-statistic: 19.47 on 6 and 31 DF, p-value: 2.853e-09
impacts(SLX, ww)
## Impact measures (SlX, estimable):
## Direct Indirect Total
## APK 2.981508e-01 5.695448e-01 8.676955e-01
## APM 3.519422e-01 -2.057335e-01 1.462087e-01
## PDRBPK 2.610947e-05 3.844804e-05 6.455751e-05
summary(impacts(SLX, ww,R=500),zstats=TRUE)
## Impact measures (SlX, estimable, n-k):
## Direct Indirect Total
## APK 2.981508e-01 5.695448e-01 8.676955e-01
## APM 3.519422e-01 -2.057335e-01 1.462087e-01
## PDRBPK 2.610947e-05 3.844804e-05 6.455751e-05
## ========================================================
## Standard errors:
## Direct Indirect Total
## APK 1.461501e-01 2.432101e-01 2.713632e-01
## APM 1.494764e-01 1.705211e-01 1.702789e-01
## PDRBPK 8.852439e-06 2.693294e-05 2.807575e-05
## ========================================================
## Z-values:
## Direct Indirect Total
## APK 2.040032 2.341781 3.1975433
## APM 2.354500 -1.206498 0.8586424
## PDRBPK 2.949410 1.427547 2.2994045
##
## p-values:
## Direct Indirect Total
## APK 0.0413472 0.019192 0.001386
## APM 0.0185477 0.227625 0.390538
## PDRBPK 0.0031838 0.153422 0.021482
AIC(SLX)
## [1] 186.778
#Uji asumsi kenormalan residual
ks.test(SLX$residuals, ecdf(SLX$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: SLX$residuals
## D = 0.026316, p-value = 1
## alternative hypothesis: two-sided
#Uji asumsi homoskedastisitas residual
bptest(SLX)
##
## studentized Breusch-Pagan test
##
## data: SLX
## BP = 3.9072, df = 6, p-value = 0.6892
#Uji asumsi non-autokorelasi residual
moran.test(residuals(SLX), ww,randomisation=T, alternative="greater")
##
## Moran I test under randomisation
##
## data: residuals(SLX)
## weights: ww
##
## Moran I statistic standard deviate = -0.18584, p-value = 0.5737
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.05102525 -0.02702703 0.01667492
#Model SAR
lag=lagsarlm(regOLS, listw=ww)
summary(lag)
##
## Call:lagsarlm(formula = regOLS, listw = ww)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.67246 -1.66639 -0.47512 1.14172 6.42440
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.8305e+00 1.0104e+01 -0.7750 0.438338
## APK 3.9866e-01 1.4081e-01 2.8312 0.004638
## APM 2.3834e-01 1.3604e-01 1.7520 0.079775
## PDRBPK 2.7219e-05 8.9613e-06 3.0374 0.002386
##
## Rho: 0.29492, LR test value: 4.0019, p-value: 0.045448
## Asymptotic standard error: 0.12232
## z-value: 2.411, p-value: 0.015908
## Wald statistic: 5.813, p-value: 0.015908
##
## Log likelihood: -90.44862 for lag model
## ML residual variance (sigma squared): 6.6665, (sigma: 2.582)
## Number of observations: 38
## Number of parameters estimated: 6
## AIC: 192.9, (AIC for lm: 194.9)
## LM test for residual autocorrelation
## test value: 2.648, p-value: 0.10368
impacts(lag,listw=ww)
## Impact measures (lag, exact):
## Direct Indirect Total
## APK 4.092854e-01 1.561251e-01 5.654105e-01
## APM 2.446947e-01 9.334068e-02 3.380354e-01
## PDRBPK 2.794459e-05 1.065968e-05 3.860427e-05
summary(impacts(lag, listw=ww,R=500),zstats=TRUE)
## Impact measures (lag, exact):
## Direct Indirect Total
## APK 4.092854e-01 1.561251e-01 5.654105e-01
## APM 2.446947e-01 9.334068e-02 3.380354e-01
## PDRBPK 2.794459e-05 1.065968e-05 3.860427e-05
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 4.109e-01 1.438e-01 6.430e-03 6.430e-03
## APM 2.389e-01 1.355e-01 6.060e-03 5.691e-03
## PDRBPK 2.849e-05 9.696e-06 4.336e-07 3.854e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 1.608e-01 3.190e-01 4.065e-01 5.021e-01 7.151e-01
## APM -4.393e-02 1.489e-01 2.470e-01 3.271e-01 5.034e-01
## PDRBPK 9.761e-06 2.179e-05 2.839e-05 3.504e-05 4.804e-05
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 1.690e-01 1.255e-01 5.611e-03 5.611e-03
## APM 8.754e-02 6.759e-02 3.023e-03 3.023e-03
## PDRBPK 1.122e-05 7.101e-06 3.176e-07 3.176e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 1.846e-02 8.376e-02 1.485e-01 2.226e-01 4.838e-01
## APM -1.676e-02 4.170e-02 7.657e-02 1.250e-01 2.443e-01
## PDRBPK 8.288e-07 6.096e-06 1.011e-05 1.478e-05 2.935e-05
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 5.799e-01 2.381e-01 1.065e-02 1.065e-02
## APM 3.265e-01 1.849e-01 8.267e-03 7.719e-03
## PDRBPK 3.971e-05 1.437e-05 6.428e-07 5.921e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 1.842e-01 4.291e-01 5.579e-01 7.036e-01 1.119e+00
## APM -6.306e-02 2.064e-01 3.397e-01 4.450e-01 6.780e-01
## PDRBPK 1.473e-05 3.028e-05 3.874e-05 4.893e-05 6.838e-05
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## APK 1.437816e-01 1.254704e-01 2.381096e-01
## APM 1.354954e-01 6.758523e-02 1.848592e-01
## PDRBPK 9.696448e-06 7.101429e-06 1.437442e-05
##
## Simulated z-values:
## Direct Indirect Total
## APK 2.857854 1.346730 2.435356
## APM 1.763331 1.295233 1.766003
## PDRBPK 2.938098 1.580345 2.762673
##
## Simulated p-values:
## Direct Indirect Total
## APK 0.0042652 0.17807 0.014877
## APM 0.0778446 0.19524 0.077395
## PDRBPK 0.0033023 0.11403 0.005733
AIC(lag)
## [1] 192.8972
#Uji asumsi kenormalan residual
ks.test(lag$residuals, ecdf(lag$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: lag$residuals
## D = 0.026316, p-value = 1
## alternative hypothesis: two-sided
#Uji asumsi homoskedastisitas residual
bptest.Sarlm(lag)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 3.0781, df = 3, p-value = 0.3797
#Uji asumsi non-autokorelasi residual
moran.test(residuals(lag), ww,randomisation=T, alternative="greater")
##
## Moran I test under randomisation
##
## data: residuals(lag)
## weights: ww
##
## Moran I statistic standard deviate = -1.023, p-value = 0.8468
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.15887332 -0.02702703 0.01661083
#Mencari R Square model lagSarlm
n_sar = length(lag$residuals)
dfres_sar = n_sar - (lag$parameters - 1)
R2Ajd_sar = 1- sum(sum(lag$residual^2)/dfres_sar)/var(lag$y)
R2Ajd_sar
## [1] 0.7007918
#Model SARMA
sarma=sacsarlm(regOLS, listw=ww, data=modku)
summary(sarma)
##
## Call:sacsarlm(formula = regOLS, data = modku, listw = ww)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.90060 -1.63451 -0.47752 0.87239 5.15425
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4828e+01 7.1031e+00 -2.0876 0.0368332
## APK 4.6081e-01 1.1329e-01 4.0675 4.751e-05
## APM 1.1967e-01 9.4542e-02 1.2658 0.2055982
## PDRBPK 2.8477e-05 7.9788e-06 3.5691 0.0003583
##
## Rho: 0.44197
## Asymptotic standard error: 0.12107
## z-value: 3.6504, p-value: 0.00026185
## Lambda: -0.59463
## Asymptotic standard error: 0.19445
## z-value: -3.0581, p-value: 0.0022276
##
## LR test value: 8.3444, p-value: 0.015418
##
## Log likelihood: -88.27738 for sac model
## ML residual variance (sigma squared): 5.1949, (sigma: 2.2792)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 190.55, (AIC for lm: 194.9)
impacts(sarma,listw=ww)
## Impact measures (sac, exact):
## Direct Indirect Total
## APK 4.914428e-01 3.343393e-01 8.257822e-01
## APM 1.276216e-01 8.682376e-02 2.144454e-01
## PDRBPK 3.036982e-05 2.066126e-05 5.103108e-05
summary(impacts(sarma, listw=ww,R=500),zstats=TRUE)
## Impact measures (sac, exact):
## Direct Indirect Total
## APK 4.914428e-01 3.343393e-01 8.257822e-01
## APM 1.276216e-01 8.682376e-02 2.144454e-01
## PDRBPK 3.036982e-05 2.066126e-05 5.103108e-05
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 5.005e-01 1.213e-01 5.423e-03 5.423e-03
## APM 1.236e-01 1.022e-01 4.572e-03 4.572e-03
## PDRBPK 3.048e-05 8.358e-06 3.738e-07 3.738e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 2.688e-01 4.159e-01 4.984e-01 5.857e-01 0.7379830
## APM -9.156e-02 5.544e-02 1.272e-01 1.903e-01 0.3071502
## PDRBPK 1.439e-05 2.483e-05 3.055e-05 3.601e-05 0.0000454
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 3.690e-01 1.999e-01 8.939e-03 8.939e-03
## APM 7.819e-02 8.690e-02 3.886e-03 3.886e-03
## PDRBPK 2.234e-05 1.227e-05 5.486e-07 5.955e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 1.206e-01 2.322e-01 3.240e-01 4.612e-01 8.304e-01
## APM -1.042e-01 3.214e-02 7.435e-02 1.224e-01 2.389e-01
## PDRBPK 5.511e-06 1.383e-05 2.057e-05 2.802e-05 5.236e-05
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## APK 8.695e-01 2.701e-01 1.208e-02 1.208e-02
## APM 2.018e-01 1.785e-01 7.982e-03 7.982e-03
## PDRBPK 5.282e-05 1.743e-05 7.793e-07 8.525e-07
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## APK 4.656e-01 6.810e-01 8.293e-01 1.027e+00 1.510e+00
## APM -1.918e-01 9.406e-02 2.174e-01 3.127e-01 5.251e-01
## PDRBPK 2.416e-05 4.055e-05 5.119e-05 6.317e-05 9.127e-05
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## APK 1.212614e-01 1.998930e-01 2.701247e-01
## APM 1.022417e-01 8.690405e-02 1.784817e-01
## PDRBPK 8.358077e-06 1.226652e-05 1.742674e-05
##
## Simulated z-values:
## Direct Indirect Total
## APK 4.127658 1.8458696 3.218891
## APM 1.208690 0.8996722 1.130444
## PDRBPK 3.647321 1.8210142 3.031094
##
## Simulated p-values:
## Direct Indirect Total
## APK 3.6648e-05 0.064911 0.0012869
## APM 0.22678206 0.368295 0.2582890
## PDRBPK 0.00026499 0.068605 0.0024367
AIC(sarma)
## [1] 190.5548
#Uji asumsi kenormalan residual
ks.test(sarma$residuals, ecdf(sarma$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: sarma$residuals
## D = 0.026316, p-value = 1
## alternative hypothesis: two-sided
#Uji asumsi homoskedastisitas residual
bptest.Sarlm(sarma)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 1.5535, df = 3, p-value = 0.67
#Uji asumsi non-autokorelasi residual
moran.test(residuals(sarma), ww,randomisation=T, alternative="greater")
##
## Moran I test under randomisation
##
## data: residuals(sarma)
## weights: ww
##
## Moran I statistic standard deviate = -0.38097, p-value = 0.6484
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.07653337 -0.02702703 0.01688667
#Mencari R Square model SARMA
n_sarma = length(sarma$residuals)
dfres_sarma = n_sarma - (sarma$parameters - 1)
R2Ajd_sarma = 1- sum(sum(sarma$residual^2)/dfres_sarma)/var(sarma$y)
R2Ajd_sarma
## [1] 0.759557
Dari ketiga jenis model tersebut, berikut ringkasan model untuk pemilihan model terbaik:
 |
| Ringkasan model |
Dari ketiga model tersebut, dapat disimpulkan bahwa model SARMA adalah model yang terbaik untuk digunakan dalam analisis regresi spasial pada kasus ini. Demikian sedikit sharing kita kali ini. Jangan lupa untuk terus menyimak setiap unggahan terbaru dan menarik pada blog sederhana ini. Semoga sedikit ini dapat bermanfaat dan membantu kebutuhan pembaca setiap blog ini. Selamat memahami dan mempraktikkan!
Tidak ada komentar: