Langsung ke konten utama

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")

plot of chunk unnamed-chunk-34

#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)

plot of chunk unnamed-chunk-37

#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!

Komentar

Postingan populer dari blog ini

Mencari P - Value dan Titik Kritis Uji F, Uji t, Uji Chi Square, dan Uji Z Normal dengan R

Mencari nilai p-value dan titik kritis Bagi teman-teman yang pernah mengenal statistika, pasti familier dengan istilah p-value dan titik kritis. P-value biasanya didefinisikan sebagai probabiltas atau peluang maksimal yang diamati dari hasil uji statistik, bahasa gampangnya adalah besarnya kesalahan penelitian berdasarkan uji statistik. Sebagai contoh sederhana, dari 100 orang dengan nama masing-masing dan diklasifikasikan ke dalam gender nama perempuan dan nama laki-laki, didapatkan nilai p-value uji statistiknya sebesar 0,05 atau 5%. Itu artinya, dari 100 orang, ada kemungkinan sebanyak 5 orang yang namanya salah klasifikasi. Dari namanya terdeteksi sebagai nama perempuan, padahal aktualnya yang bersangkutan bergender laki-laki. Sedangkan titik kritis atau titik uji adalah nilai batas pengujian hipotesis statistik, apakah masuk dalam wilayah tolak hipotesis, ataukah gagal menolaknya. Titik ini berkaitan erat dengan nilai p-value . Kalau biasanya kita mendapatkan kedua nilai ini da...

Cara Mendowload dan Install R serta RStudio di Windows (Step by Step)

Cara Download dan Install R serta R Studio di Windows Halo teman-teman, mohon maaf karena beberapa waktu ini, blog ini sempat vakum dari unggahan. Kali ini saya akan coba berbagai mengenai bagaimana cara mengunduh ( download ) dan menginstal ( install ) program R sekaligus R Studio khususnya di Windows. Unggahan kali ini sedikit terbalik karena semestinya saya unggah terlebih dahulu pertama kali di blog ini, namun bukan masalah, mengingat kemarin ada beberapa pihak yang meminta untuk menerangkap bagaimana tahapan mengunduh dan instalasi R dan R Studio, jadinya saya dahulukan pada unggahan ini sebelum pembahasan mengenai Data Mining , Data Science , atau bahasan Big Data kita terlampau jauh. Baik, kita akan mulai dengan bagaimana mengunduh R dan R Studio melalui mesin pencari Google. R dan R Studio ini memang beberapa waktu terakhir ini booming , apalagi dengan munculnya konsep mengenai Big Data , Data Modelling, Data Mining, dan Data Science serta Data Visualization . Sebenarnya, men...

Analisis Tipologi Klassen (Klassen Typology) dan Visualisasi Spasialnya dengan R

Tipologi Klassen dan visualisasinya dengan R Halo teman-teman, sebelumnya kita telah membahas tentang analisis Shift Share dan Location Quotient (LQ) dengan menggunakan R. Kali ini, kita akan membahas mengenai satu lagi alat analisis yang sebenarnya merupakan alat analisis tiga serangkai dari SS dan LQ, yaitu analisis Tipologi Klassen. Dalam penelitian ekonomi kewilayahan, ketiga analisis ini seringkali digunakan, baik dalam rangka melihat perkembangan dan transformasi struktur ekonomi suatu wilayah maupun melihat keunggulan kompetitif dan keunggulan komparatif wilayah satu dengan wilayah lainnya dengan mengacu wilayah referensi. Terlebih dulu, sebelum melakukan visualisasi spasial menggunakan fungsi plot(), ada baiknya kita bahas terlebih dahulu mengenai Tipologi Klassen itu sendiri. Tipologi Klassen merupakan teknik pengelompokan sektor, subsektor, lapangan usaha, atau komoditas tertentu di wilayah analisis berdasarkan pertumbuhan nilai tambah wilayah analisis terhadap nasional atau...