Pemodelan Spatial LAG Model (SAR) pada Rumah dengan sanitasi layak di Provinsi Jawa Barat
Azmi
2024-12-11
INPUT DATA DAN PLOT PETA
setwd("c:/Materi Kuliah/Spatial/Peta Jabar")
jabar1=readShapePoly("Jabar.shp")## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among othersplot(jabar1)
KAB = factor(c("Kab Bogor", "Kab Sukabumi", "Kab Cianjur", "Kab Bandung", "Kab Garut",
"Kab Tasikmalaya", "Kab Ciamis", "Kab Kuningan", "Kab Cirebon", "Kab Majalengka",
"Kab Sumedang", "Kab Indramayu", "Kab Subang", "Kab Purwakarta", "Kab Karawang",
"Kab Bekasi", "Kab Bandung Barat", "Kab Pangandaran", "Kota Bogor", "Kota Sukabumi",
"Kota Bandung", "Kota Cirebon", "Kota Bekasi", "Kota Depok", "Kota Cimahi",
"Kota Tasikmalaya", "Kota Banjar"))
jabar1$KAB = KAB
ID = c(1:27)PLOT PERSEBARAN TIAP VARIABEL
setwd("C:/Materi Kuliah/Spatial")
sani <- read.csv("sanitasi.csv", sep = ",", dec =",")
head(sani)## ID Kabupaten.Kota Rumah.dengan.sanitasi.layak Indeks.Pendidikan
## 1 1 Bogor 49.98 62.52
## 2 2 Sukabumi 42.48 57.73
## 3 3 Cianjur 36.57 57.36
## 4 4 Bandung 52.58 65.57
## 5 5 Garut 34.24 59.85
## 6 6 Tasikmalaya 40.54 60.74
## Indeks.Kesehatan Rata.rata.pengeluaran.perkapita
## 1 79.46 709091
## 2 79.29 629476
## 3 77.82 508099
## 4 83.09 570272
## 5 79.77 522212
## 6 76.85 588690
## Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak
## 1 92.28
## 2 85.70
## 3 80.79
## 4 93.79
## 5 78.15
## 6 89.35str(sani)## 'data.frame': 27 obs. of 7 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Kabupaten.Kota : chr "Bogor" "Sukabumi" "Cianjur" "Bandung" ...
## $ Rumah.dengan.sanitasi.layak : num 50 42.5 36.6 52.6 34.2 ...
## $ Indeks.Pendidikan : num 62.5 57.7 57.4 65.6 59.9 ...
## $ Indeks.Kesehatan : num 79.5 79.3 77.8 83.1 79.8 ...
## $ Rata.rata.pengeluaran.perkapita : int 709091 629476 508099 570272 522212 588690 653024 621212 636984 555269 ...
## $ Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak: num 92.3 85.7 80.8 93.8 78.2 ...jabar1$sanitasi=sani$Rumah.dengan.sanitasi.layak
jabar1$pendidikan=sani$Indeks.Pendidikan
jabar1$kesehatan=sani$Indeks.Kesehatan
jabar1$pengeluaran=sani$Rata.rata.pengeluaran.perkapita
jabar1$airbersih=sani$Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak
tmap_options(check.and.fix = TRUE)
par(mfrow=c(2,2))1. Persebaran Jumlah Rumah dengan Sanitasi Layak
tm_shape(jabar1) + tm_polygons("sanitasi",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)## Warning: Currect projection of shape jabar1 unknown. Long-lat (WGS84) is
## assumed.## Warning: The shape jabar1 is invalid. See sf::st_is_valid## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
## 2. Persebaran Jumlah Pendidikan
tm_shape(jabar1) + tm_polygons("pendidikan",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)## Warning: Currect projection of shape jabar1 unknown. Long-lat (WGS84) is
## assumed.## Warning: The shape jabar1 is invalid. See sf::st_is_valid## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
## 3. Persebaran Jumlah Kesehatan
tm_shape(jabar1) + tm_polygons("kesehatan",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)## Warning: Currect projection of shape jabar1 unknown. Long-lat (WGS84) is
## assumed.## Warning: The shape jabar1 is invalid. See sf::st_is_valid## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
## 3. Persebaran Jumlah Kesehatan
tm_shape(jabar1) + tm_polygons("pengeluaran",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)## Warning: Currect projection of shape jabar1 unknown. Long-lat (WGS84) is
## assumed.## Warning: The shape jabar1 is invalid. See sf::st_is_valid## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
4. Persebaran Jumlah Rumah dengan ketersediaan air bersih
tm_shape(jabar1) + tm_polygons("airbersih",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)## Warning: Currect projection of shape jabar1 unknown. Long-lat (WGS84) is
## assumed.## Warning: The shape jabar1 is invalid. See sf::st_is_valid## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
Mendapatkan koordinat
CoordK = coordinates(jabar1)Mendapatkan Matriks Pembobot Spasial
W = poly2nb(jabar1, row.names = ID, queen = F)
W## Neighbour list object:
## Number of regions: 27
## Number of nonzero links: 106
## Percentage nonzero weights: 14.54047
## Average number of links: 3.925926WB = nb2mat(W, style = "B", zero.policy = T)
WB## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## 1 0 1 1 0 0 0 0 0 0 0 0 0 0 1
## 2 1 0 1 0 0 0 0 0 0 0 0 0 0 0
## 3 1 1 0 1 1 0 0 0 0 0 0 0 0 1
## 4 0 0 1 0 1 0 0 0 0 0 1 0 1 0
## 5 0 0 1 1 0 1 0 0 0 0 1 0 0 0
## 6 0 0 0 0 1 0 1 0 0 1 1 0 0 0
## 7 0 0 0 0 0 1 0 1 0 1 0 0 0 0
## 8 0 0 0 0 0 0 1 0 1 1 0 0 0 0
## 9 0 0 0 0 0 0 0 1 0 1 0 1 0 0
## 10 0 0 0 0 0 1 1 1 1 0 1 1 0 0
## 11 0 0 0 1 1 1 0 0 0 1 0 1 1 0
## 12 0 0 0 0 0 0 0 0 1 1 1 0 1 0
## 13 0 0 0 1 0 0 0 0 0 0 1 1 0 1
## 14 1 0 1 0 0 0 0 0 0 0 0 0 1 0
## 15 1 0 0 0 0 0 0 0 0 0 0 0 1 1
## 16 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 17 0 0 1 1 0 0 0 0 0 0 0 0 1 1
## 18 0 0 0 0 0 1 1 0 0 0 0 0 0 0
## 19 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 20 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## 21 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## 22 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 23 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 24 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 25 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## 26 0 0 0 0 0 1 1 0 0 0 0 0 0 0
## 27 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## 1 1 1 0 0 1 0 0 0 1 1 0 0
## 2 0 0 0 0 0 1 0 0 0 0 0 0
## 3 0 0 1 0 0 0 0 0 0 0 0 0
## 4 0 0 1 0 0 0 1 0 0 0 1 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 1 0 0 0 0 0 0 0 1
## 7 0 0 0 1 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 1 0 0 0 0
## 10 0 0 0 0 0 0 0 0 0 0 0 0
## 11 0 0 0 0 0 0 0 0 0 0 0 0
## 12 0 0 0 0 0 0 0 0 0 0 0 0
## 13 1 0 1 0 0 0 0 0 0 0 0 0
## 14 1 0 1 0 0 0 0 0 0 0 0 0
## 15 0 1 0 0 0 0 0 0 0 0 0 0
## 16 1 0 0 0 0 0 0 0 1 0 0 0
## 17 0 0 0 0 0 0 1 0 0 0 1 0
## 18 0 0 0 0 0 0 0 0 0 0 0 0
## 19 0 0 0 0 0 0 0 0 0 0 0 0
## 20 0 0 0 0 0 0 0 0 0 0 0 0
## 21 0 0 1 0 0 0 0 0 0 0 1 0
## 22 0 0 0 0 0 0 0 0 0 0 0 0
## 23 0 1 0 0 0 0 0 0 0 1 0 0
## 24 0 0 0 0 0 0 0 0 1 0 0 0
## 25 0 0 1 0 0 0 1 0 0 0 0 0
## 26 0 0 0 0 0 0 0 0 0 0 0 0
## 27 0 0 0 0 0 0 0 0 0 0 0 0
## [,27]
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
## 7 1
## 8 0
## 9 0
## 10 0
## 11 0
## 12 0
## 13 0
## 14 0
## 15 0
## 16 0
## 17 0
## 18 0
## 19 0
## 20 0
## 21 0
## 22 0
## 23 0
## 24 0
## 25 0
## 26 0
## 27 0
## attr(,"call")
## nb2mat(neighbours = W, style = "B", zero.policy = T)WL = nb2listw(W)
WL## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27
## Number of nonzero links: 106
## Percentage nonzero weights: 14.54047
## Average number of links: 3.925926
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 27 729 27 16.73552 121.5371Plot Hubungan Lokasi Peta
plot(jabar1, axes=T, col="gray90")
text(CoordK[,1],CoordK[,2],jabar1$KAB,
col="black",cex=.5, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot(W,coordinates(jabar1), add=T, col='red')
## Moran’s Test
setwd("C:/Materi Kuliah/Spatial")
sani <- read.csv("sanitasi.csv",sep = ",", dec =",")
kasus <- sani$Rumah.dengan.sanitasi.layak
moran.test(kasus, WL)##
## Moran I test under randomisation
##
## data: kasus
## weights: WL
##
## Moran I statistic standard deviate = 3.7589, p-value = 8.532e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.49820654 -0.03846154 0.02038358moran.plot(kasus, WL)
### Local Moran’s I
local <- localmoran(kasus,WL)
local## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.06211892 -0.0055307243 0.013365329 0.58516147 0.55843915
## 2 0.97522375 -0.0326145725 0.261241139 1.97183279 0.04862870
## 3 0.81666396 -0.0698551803 0.233911563 1.83299918 0.06680269
## 4 0.10810605 -0.0014095334 0.004126122 1.70492345 0.08820870
## 5 0.67010978 -0.0883886272 0.478621902 1.09637314 0.27291551
## 6 -0.03174170 -0.0432940534 0.149110842 0.02991686 0.97613336
## 7 0.06253742 -0.0028956010 0.010393979 0.64180933 0.52099699
## 8 0.39939399 -0.0062892957 0.051747851 1.78336794 0.07452643
## 9 1.62857859 -0.0813631671 0.443974620 2.56626941 0.01027990
## 10 0.67778876 -0.0585851217 0.198550459 1.65258166 0.09841603
## 11 0.11080535 -0.0266188915 0.093277174 0.44996182 0.65273797
## 12 1.55189184 -0.0744080663 0.409096746 2.54265609 0.01100135
## 13 0.38914663 -0.0341595845 0.118773746 1.22827162 0.21934502
## 14 -0.18023028 -0.0243140361 0.107606910 -0.47530379 0.63457046
## 15 0.47392140 -0.0444123222 0.252092815 1.03235539 0.30190566
## 16 0.08555905 -0.0063343054 0.052115827 0.40253096 0.68729331
## 17 0.29461445 -0.0779375503 0.258707839 0.73245696 0.46388970
## 18 -0.02554207 -0.0001622224 0.002102061 -0.55356215 0.57987854
## 19 0.25964326 -0.0180312320 0.478064880 0.40159895 0.68797920
## 20 1.33545529 -0.0808910211 2.007386923 0.99966367 0.31747330
## 21 0.76541082 -0.0685140450 0.528428529 1.14718589 0.25130480
## 22 2.10808949 -0.0807985461 2.005293809 1.54573322 0.12216901
## 23 0.01138749 -0.0001520963 0.001259166 0.32519868 0.74503071
## 24 0.12596750 -0.0124904018 0.159853716 0.34630311 0.72911493
## 25 0.10372309 -0.0004280427 0.003542677 1.74983959 0.08014600
## 26 0.27922700 -0.0193871221 0.246385950 0.60159245 0.54744546
## 27 0.39372665 -0.0791961761 1.968951829 0.33703324 0.73609184
## attr(,"call")
## localmoran(x = kasus, listw = WL)
## attr(,"class")
## [1] "localmoran" "matrix" "array"
## attr(,"quadr")
## mean median pysal
## 1 Low-Low Low-High Low-Low
## 2 Low-Low Low-Low Low-Low
## 3 Low-Low Low-Low Low-Low
## 4 Low-Low Low-Low Low-Low
## 5 Low-Low Low-Low Low-Low
## 6 Low-High Low-High Low-High
## 7 High-High High-High High-High
## 8 High-High High-High High-High
## 9 High-High High-High High-High
## 10 High-High High-High High-High
## 11 High-High High-High High-High
## 12 High-High High-High High-High
## 13 High-High High-High High-High
## 14 High-Low High-Low High-Low
## 15 High-High High-High High-High
## 16 High-High High-High High-High
## 17 Low-Low Low-Low Low-Low
## 18 High-Low High-Low High-Low
## 19 Low-Low Low-Low Low-Low
## 20 Low-Low Low-Low Low-Low
## 21 Low-Low Low-Low Low-Low
## 22 High-High High-High High-High
## 23 Low-Low Low-Low Low-Low
## 24 Low-Low Low-Low Low-Low
## 25 Low-Low Low-Low Low-Low
## 26 Low-Low Low-Low Low-Low
## 27 High-High High-High High-HighRegresi linier OLS
setwd("C:/Materi Kuliah/Spatial")
sani <- read.csv("sanitasi.csv", sep = ",", dec =",")
head(sani)## ID Kabupaten.Kota Rumah.dengan.sanitasi.layak Indeks.Pendidikan
## 1 1 Bogor 49.98 62.52
## 2 2 Sukabumi 42.48 57.73
## 3 3 Cianjur 36.57 57.36
## 4 4 Bandung 52.58 65.57
## 5 5 Garut 34.24 59.85
## 6 6 Tasikmalaya 40.54 60.74
## Indeks.Kesehatan Rata.rata.pengeluaran.perkapita
## 1 79.46 709091
## 2 79.29 629476
## 3 77.82 508099
## 4 83.09 570272
## 5 79.77 522212
## 6 76.85 588690
## Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak
## 1 92.28
## 2 85.70
## 3 80.79
## 4 93.79
## 5 78.15
## 6 89.35summary(sani)## ID Kabupaten.Kota Rumah.dengan.sanitasi.layak
## Min. : 1.0 Length:27 Min. :34.24
## 1st Qu.: 7.5 Class :character 1st Qu.:43.94
## Median :14.0 Mode :character Median :54.36
## Mean :14.0 Mean :55.23
## 3rd Qu.:20.5 3rd Qu.:67.52
## Max. :27.0 Max. :75.37
## Indeks.Pendidikan Indeks.Kesehatan Rata.rata.pengeluaran.perkapita
## Min. :56.72 Min. :76.85 Min. : 508099
## 1st Qu.:60.06 1st Qu.:79.61 1st Qu.: 627210
## Median :62.52 Median :80.97 Median : 671250
## Mean :64.95 Mean :81.11 Mean : 695377
## 3rd Qu.:70.06 3rd Qu.:83.11 3rd Qu.: 730263
## Max. :77.33 Max. :85.35 Max. :1053945
## Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak
## Min. :78.15
## 1st Qu.:92.44
## Median :95.51
## Mean :93.71
## 3rd Qu.:97.92
## Max. :99.44y=sani$Rumah.dengan.sanitasi.layak
x1=sani$Indeks.Pendidikan
x2=sani$Indeks.Kesehatan
x3=sani$Rata.rata.pengeluaran.perkapita
x4=sani$Rumah.Tangga.yang.Memiliki.Akses.Terhadap.Sumber.Air.Minum.Layak
modelols=lm(y~x1+x2+x3+x4, data=sani)
summary(modelols)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.7388 -8.0378 -0.2753 6.7542 18.6617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.200e+01 1.009e+02 0.218 0.82936
## x1 -1.666e+00 4.854e-01 -3.431 0.00239 **
## x2 -8.763e-01 1.480e+00 -0.592 0.55981
## x3 2.667e-05 2.502e-05 1.066 0.29810
## x4 2.070e+00 4.241e-01 4.880 7.05e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.889 on 22 degrees of freedom
## Multiple R-squared: 0.5844, Adjusted R-squared: 0.5088
## F-statistic: 7.734 on 4 and 22 DF, p-value: 0.0004746AIC(modelols)## [1] 206.8302plot(modelols)
## Pengujian Asumsi OLS
1. Normalitas Residual
shapiro.test(modelols$residuals)##
## Shapiro-Wilk normality test
##
## data: modelols$residuals
## W = 0.96656, p-value = 0.514hist(modelols$residuals)
### 2. Homoskedastisitas
bptest(modelols)##
## studentized Breusch-Pagan test
##
## data: modelols
## BP = 8.7518, df = 4, p-value = 0.067613. Multikolinearitas
library(car)## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:maptools':
##
## pointLabelvif(modelols)## x1 x2 x3 x4
## 2.845256 2.685698 2.747865 1.6332274. Autokorelasi spasial
moran.test(modelols$residuals,listw=WL, alternative="two.sided")##
## Moran I test under randomisation
##
## data: modelols$residuals
## weights: WL
##
## Moran I statistic standard deviate = 1.4874, p-value = 0.1369
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.17203920 -0.03846154 0.02002882durbinWatsonTest(modelols)## lag Autocorrelation D-W Statistic p-value
## 1 -0.2064719 2.303557 0.64
## Alternative hypothesis: rho != 0LM Test
LM<-lm.LMtests(modelols, WL, test="all")## Please update scripts to use lm.RStests in place of lm.LMtestsprint(LM)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
## test weights: listw
##
## RSerr = 1.2893, df = 1, p-value = 0.2562
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
## test weights: listw
##
## RSlag = 5.6858, df = 1, p-value = 0.0171
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
## test weights: listw
##
## adjRSerr = 3.7924, df = 1, p-value = 0.05149
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
## test weights: listw
##
## adjRSlag = 8.1889, df = 1, p-value = 0.004215
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sani)
## test weights: listw
##
## SARMA = 9.4782, df = 2, p-value = 0.008747Spatial lag model (SAR)
sar.sani<-lagsarlm(y~x1+x2+x3+x4, data=sani, WL)
summary(sar.sani)##
## Call:lagsarlm(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.45232 -6.31053 0.25799 5.55722 12.36659
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0709e+01 7.7215e+01 0.1387 0.889691
## x1 -1.0921e+00 3.9424e-01 -2.7701 0.005603
## x2 -8.3562e-01 1.1308e+00 -0.7390 0.459929
## x3 2.0631e-05 1.9258e-05 1.0713 0.284041
## x4 1.4927e+00 3.4736e-01 4.2973 1.729e-05
##
## Rho: 0.52688, LR test value: 6.8378, p-value: 0.008925
## Asymptotic standard error: 0.15413
## z-value: 3.4184, p-value: 0.00062986
## Wald statistic: 11.686, p-value: 0.00062986
##
## Log likelihood: -93.99621 for lag model
## ML residual variance (sigma squared): 57.093, (sigma: 7.556)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 201.99, (AIC for lm: 206.83)
## LM test for residual autocorrelation
## test value: 1.3981, p-value: 0.23703sar2sls.sani<-stsls(y~x1+x2+x3+x4, data=sani, WL)
summary(sar2sls.sani)##
## Call:stsls(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.10597 -5.10108 -0.41626 3.74788 13.82514
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Rho 9.1553e-01 2.5375e-01 3.6081 0.0003085
## (Intercept) 2.3800e+00 8.6712e+01 0.0274 0.9781029
## x1 -6.6914e-01 4.9971e-01 -1.3390 0.1805557
## x2 -8.0562e-01 1.2698e+00 -0.6344 0.5257905
## x3 1.6178e-05 2.1663e-05 0.7468 0.4551718
## x4 1.0672e+00 4.5777e-01 2.3312 0.0197414
##
## Residual variance (sigma squared): 71.978, (sigma: 8.484)Uji Asumsi SAR
1. Normalitas Residual
shapiro.test(sar.sani$residuals)##
## Shapiro-Wilk normality test
##
## data: sar.sani$residuals
## W = 0.96166, p-value = 0.4029hist(sar.sani$residuals)
### 2. Homoskedastisitas
# Ekstraksi residual dari model SAR
residuals_sar <- residuals(sar.sani)
# Buat model linier biasa untuk variabel independen
lm_model <- lm(residuals_sar ~ x1 + x2 + x3 + x4, data = sani)
# Uji Breusch-Pagan
bptest(lm_model)##
## studentized Breusch-Pagan test
##
## data: lm_model
## BP = 6.3517, df = 4, p-value = 0.17443. Autokorelasi spasial
moran.test(sar.sani$residuals,listw=WL, alternative="two.sided")##
## Moran I test under randomisation
##
## data: sar.sani$residuals
## weights: WL
##
## Moran I statistic standard deviate = -0.42185, p-value = 0.6731
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## -0.09833927 -0.03846154 0.02014680Sebaran Parameter SAR
sar.sani <- lagsarlm(y ~ x1 + x2 + x3 + x4, data = sani, WL)
summary(sar.sani)##
## Call:lagsarlm(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.45232 -6.31053 0.25799 5.55722 12.36659
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0709e+01 7.7215e+01 0.1387 0.889691
## x1 -1.0921e+00 3.9424e-01 -2.7701 0.005603
## x2 -8.3562e-01 1.1308e+00 -0.7390 0.459929
## x3 2.0631e-05 1.9258e-05 1.0713 0.284041
## x4 1.4927e+00 3.4736e-01 4.2973 1.729e-05
##
## Rho: 0.52688, LR test value: 6.8378, p-value: 0.008925
## Asymptotic standard error: 0.15413
## z-value: 3.4184, p-value: 0.00062986
## Wald statistic: 11.686, p-value: 0.00062986
##
## Log likelihood: -93.99621 for lag model
## ML residual variance (sigma squared): 57.093, (sigma: 7.556)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 201.99, (AIC for lm: 206.83)
## LM test for residual autocorrelation
## test value: 1.3981, p-value: 0.23703# Tambahkan koefisien ke data spatial
jabar1$coef_x1 <- rep(sar.sani$coefficients["x1"], nrow(jabar1))
jabar1$coef_x2 <- rep(sar.sani$coefficients["x2"], nrow(jabar1))
jabar1$coef_x3 <- rep(sar.sani$coefficients["x3"], nrow(jabar1))
jabar1$coef_x4 <- rep(sar.sani$coefficients["x4"], nrow(jabar1))
# Visualisasi peta untuk masing-masing koefisien
library(RColorBrewer)
library(sp)
spplot(jabar1, "coef_x1", main = "Koefisien SAR untuk x1", col.regions = terrain.colors(100))
spplot(jabar1, "coef_x2", main = "Koefisien SAR untuk x2", col.regions = terrain.colors(100))
spplot(jabar1, "coef_x3", main = "Koefisien SAR untuk x3", col.regions = terrain.colors(100))
spplot(jabar1, "coef_x4", main = "Koefisien SAR untuk x4", col.regions = terrain.colors(100))
## Sebaran Residu dan Predict SAR
sar.sani <- lagsarlm(y ~ x1 + x2 + x3 + x4, data = sani, WL)
summary(sar.sani)##
## Call:lagsarlm(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.45232 -6.31053 0.25799 5.55722 12.36659
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0709e+01 7.7215e+01 0.1387 0.889691
## x1 -1.0921e+00 3.9424e-01 -2.7701 0.005603
## x2 -8.3562e-01 1.1308e+00 -0.7390 0.459929
## x3 2.0631e-05 1.9258e-05 1.0713 0.284041
## x4 1.4927e+00 3.4736e-01 4.2973 1.729e-05
##
## Rho: 0.52688, LR test value: 6.8378, p-value: 0.008925
## Asymptotic standard error: 0.15413
## z-value: 3.4184, p-value: 0.00062986
## Wald statistic: 11.686, p-value: 0.00062986
##
## Log likelihood: -93.99621 for lag model
## ML residual variance (sigma squared): 57.093, (sigma: 7.556)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 201.99, (AIC for lm: 206.83)
## LM test for residual autocorrelation
## test value: 1.3981, p-value: 0.23703# Menambahkan prediksi dan residu ke shapefile
sani$predicted <- sar.sani$fitted.values
sani$residuals <- residuals(sar.sani)
# Menambahkan prediksi dan residu ke shapefile spatial
jabar1$predicted <- sani$predicted
jabar1$residuals <- sani$residuals
# Plot prediksi pada peta
library(RColorBrewer)
library(sp)
spplot(jabar1, "predicted", main = "Prediksi SAR", col.regions = brewer.pal(9, "Blues"))
# Plot residu pada peta
spplot(jabar1, "residuals", main = "Residu SAR", col.regions = brewer.pal(9, "Reds"))
Spatial Error model
errorsalm.sani<-errorsarlm(y ~ x1 + x2 + x3+x4, data=sani, WL)
summary(errorsalm.sani)##
## Call:errorsarlm(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.6471 -5.9460 1.5220 6.6831 13.4919
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.9840e+01 8.7466e+01 0.3412 0.732983
## x1 -1.3732e+00 4.4695e-01 -3.0723 0.002124
## x2 -7.1556e-01 1.2471e+00 -0.5738 0.566109
## x3 2.9474e-05 2.3131e-05 1.2742 0.202581
## x4 1.6297e+00 4.0087e-01 4.0655 4.793e-05
##
## Lambda: 0.49205, LR test value: 2.5957, p-value: 0.10716
## Asymptotic standard error: 0.18403
## z-value: 2.6737, p-value: 0.0075022
## Wald statistic: 7.1486, p-value: 0.0075022
##
## Log likelihood: -96.11727 for error model
## ML residual variance (sigma squared): 67.579, (sigma: 8.2206)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 206.23, (AIC for lm: 206.83)fgls.sani<-GMerrorsar(y ~ x1 + x2 + x3+x4, data=sani, WL)
summary(fgls.sani)##
## Call:GMerrorsar(formula = y ~ x1 + x2 + x3 + x4, data = sani, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.0521 -8.3841 1.2029 6.0168 18.7543
##
## Type: GM SAR estimator
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.8281e+01 9.0235e+01 0.3134 0.7539635
## x1 -1.4825e+00 4.4930e-01 -3.2995 0.0009685
## x2 -7.7243e-01 1.2997e+00 -0.5943 0.5522995
## x3 2.9044e-05 2.3281e-05 1.2475 0.2122139
## x4 1.7720e+00 3.9930e-01 4.4379 9.085e-06
##
## Lambda: 0.34192 (standard error): 0.79228 (z-value): 0.43157
## Residual variance (sigma squared): 72.709, (sigma: 8.5269)
## GM argmin sigma squared: 69.006
## Number of observations: 27
## Number of parameters estimated: 7Menghitung Impact
impacts(sar.sani, listw=WL)## Impact measures (lag, exact):
## Direct Indirect Total
## x1 -1.192100e+00 -1.116218e+00 -2.308318e+00
## x2 -9.121319e-01 -8.540714e-01 -1.766203e+00
## x3 2.251988e-05 2.108641e-05 4.360629e-05
## x4 1.629356e+00 1.525642e+00 3.154998e+00#Perbandingan AIC untuk model
AIC(modelols)## [1] 206.8302AIC(sar.sani)## [1] 201.9924AIC(errorsalm.sani)## [1] 206.2345