Spatial Econometrics Faktor Persentase Penduduk Miskin Jawa Barat tahun 2023
Aisya Putri Syanurli
2024-12-15
email: aisyasyanurli@gmail.com
Dosen pengampu: Dr. I Gede Nyoman Mindra Jaya, M.Si.
Set working directory
setwd("C:\\Users\\Aisya\\Downloads")Load Library
library(dplyr) ##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(scales)
library(magrittr)
library(sf)## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(ncdf4)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ✔ readr 2.1.5## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ purrr::discard() masks scales::discard()
## ✖ tidyr::extract() masks magrittr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::set_names() masks magrittr::set_names()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(reticulate)
library(gstat)
library(sp)
library(raster)## Warning: package 'raster' was built under R version 4.4.1##
## Attaching package: 'raster'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(Metrics)
library(gganimate)## No renderer backend detected. gganimate will default to writing frames to separate files
## Consider installing:
## - the `gifski` package for gif output
## - the `av` package for video output
## and restarting the R session
##
## Attaching package: 'gganimate'
##
## The following object is masked from 'package:raster':
##
## animatelibrary(spdep)## Loading required package: spData## Warning: package 'spData' was built under R version 4.4.1## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spatialreg)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'spatialreg'
##
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(olsrr)##
## Attaching package: 'olsrr'
##
## The following object is masked from 'package:datasets':
##
## riverslibrary(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:purrr':
##
## some
##
## The following object is masked from 'package:dplyr':
##
## recodelibrary(corrplot)## corrplot 0.92 loadedlibrary(lmtest)## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(tmaptools)Input the map
# Input the map
Jabar <- read_sf('C:\\Users\\Aisya\\Downloads\\[geosai.my.id]Jawa_Barat_Kab\\Jawa_Barat_ADMIN_BPS.shp')
Jabar <- Jabar[Jabar$Kabupaten != "Waduk Cirata",]
Jabar <- as_Spatial(Jabar)
class(Jabar)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"str(slot(Jabar,"data"))## 'data.frame': 27 obs. of 6 variables:
## $ ADM0_EN : chr "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...
## $ date : Date, format: "2019-12-20" "2019-12-20" ...
## $ validOn : Date, format: "2020-04-01" "2020-04-01" ...
## $ PROVINCE : chr "Jawa Barat" "Jawa Barat" "Jawa Barat" "Jawa Barat" ...
## $ Kabupaten: chr "Bandung" "Bandung Barat" "Bekasi" "Bogor" ...
## $ PRV2 : chr "Jawa_Barat" "Jawa_Barat" "Jawa_Barat" "Jawa_Barat" ...# Coordinate
IDs <- Jabar$Kabupaten
coords<-coordinates(Jabar)Input data
# Input data
data <- read.csv("Variabel UAS Spasial (Jabar).csv", sep = ";", dec = ".") ; data## Kab.Kota PDRB Tingkat.Partisipasi.Angkatan.Kerja
## 1 Bandung 92830.17 67.10
## 2 Bandung Barat 35062.19 67.01
## 3 Bekasi 279224.90 65.00
## 4 Bogor 176683.58 64.22
## 5 Ciamis 25111.71 66.26
## 6 Cianjur 36339.55 72.31
## 7 Cirebon 37246.57 66.16
## 8 Garut 44087.22 70.10
## 9 Indramayu 67237.69 63.88
## 10 Karawang 187051.65 63.40
## 11 Kota Bandung 221969.13 66.97
## 12 Kota Banjar 3679.08 67.44
## 13 Kota Bekasi 77241.79 64.65
## 14 Kota Bogor 37055.36 64.81
## 15 Kota Cimahi 25931.98 68.43
## 16 Kota Cirebon 18934.40 68.71
## 17 Kota Depok 55221.82 62.76
## 18 Kota Sukabumi 9801.88 62.57
## 19 Kota Tasikmalaya 17781.93 65.44
## 20 Kuningan 19418.73 61.95
## 21 Majalengka 25793.24 68.50
## 22 Pangandaran 8869.86 80.15
## 23 Purwakarta 51740.50 66.37
## 24 Subang 31604.98 70.03
## 25 Sukabumi 52993.98 67.75
## 26 Sumedang 26876.86 67.76
## 27 Tasikmalaya 27599.16 68.37
## Tingkat.Pengangguran.Terbuka Jumlah.Panjang.Jalan Laju.Pertumbuhan.penduduk
## 1 6.52 1387 0.92
## 2 8.11 794 1.39
## 3 8.87 1138 1.36
## 4 8.47 2019 1.27
## 5 3.52 1215 0.66
## 6 7.71 1729 1.12
## 7 7.65 1414 1.25
## 8 7.33 1211 1.31
## 9 6.46 1089 1.12
## 10 8.95 1941 1.19
## 11 8.83 1382 0.90
## 12 5.43 277 1.16
## 13 7.90 335 1.12
## 14 9.39 210 0.92
## 15 10.52 119 1.40
## 16 7.66 175 0.90
## 17 6.97 2003 1.50
## 18 8.53 160 1.47
## 19 6.55 499 1.27
## 20 9.49 923 1.07
## 21 4.12 1085 0.95
## 22 1.52 756 0.67
## 23 7.72 878 1.36
## 24 7.65 1366 0.87
## 25 7.32 1938 0.98
## 26 6.94 969 0.80
## 27 3.89 1493 0.78
## Jumlah.Tenaga.Kerja.Usaha.Mikro.dan.Kecil Persentase.Penduduk.Miskin IPM
## 1 205390 6.40 74.03
## 2 75317 10.52 70.33
## 3 53553 4.93 76.13
## 4 84589 7.27 73.02
## 5 67617 7.42 73.12
## 6 111340 10.22 68.18
## 7 67028 11.20 71.81
## 8 89786 9.77 69.22
## 9 31636 12.13 70.19
## 10 53214 7.87 73.25
## 11 118495 3.96 83.29
## 12 9862 6.14 74.45
## 13 40538 4.10 83.06
## 14 30860 6.67 78.36
## 15 23553 4.66 79.69
## 16 8872 9.16 77.45
## 17 31238 2.38 82.53
## 18 7569 7.50 77.16
## 19 40563 11.53 75.47
## 20 44148 12.12 70.99
## 21 50683 11.21 70.76
## 22 44367 8.98 70.57
## 23 26378 8.46 73.43
## 24 31563 9.52 71.42
## 25 85014 7.01 69.71
## 26 54347 9.36 74.02
## 27 91495 10.28 69.38
## Rata.rata.Lama.Sekolah
## 1 9.10
## 2 8.23
## 3 9.57
## 4 8.37
## 5 8.09
## 6 7.22
## 7 7.64
## 8 7.84
## 9 6.94
## 10 8.04
## 11 11.06
## 12 8.79
## 13 11.66
## 14 10.64
## 15 11.39
## 16 10.37
## 17 11.58
## 18 10.37
## 19 9.54
## 20 7.89
## 21 7.52
## 22 8.04
## 23 8.13
## 24 7.45
## 25 7.33
## 26 8.73
## 27 7.96attach(data)# Checking missing value
colSums(is.na(data))## Kab.Kota
## 0
## PDRB
## 0
## Tingkat.Partisipasi.Angkatan.Kerja
## 0
## Tingkat.Pengangguran.Terbuka
## 0
## Jumlah.Panjang.Jalan
## 0
## Laju.Pertumbuhan.penduduk
## 0
## Jumlah.Tenaga.Kerja.Usaha.Mikro.dan.Kecil
## 0
## Persentase.Penduduk.Miskin
## 0
## IPM
## 0
## Rata.rata.Lama.Sekolah
## 0#Define datasets
Jabar@data <- data[,-1]
summary(Jabar)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 106.370304 108.846881
## y -7.820979 -5.918148
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## PDRB Tingkat.Partisipasi.Angkatan.Kerja
## Min. : 3679 Min. :61.95
## 1st Qu.: 25452 1st Qu.:64.73
## Median : 36340 Median :66.97
## Mean : 62718 Mean :66.97
## 3rd Qu.: 61230 3rd Qu.:68.40
## Max. :279225 Max. :80.15
## Tingkat.Pengangguran.Terbuka Jumlah.Panjang.Jalan Laju.Pertumbuhan.penduduk
## Min. : 1.520 Min. : 119.0 Min. :0.66
## 1st Qu.: 6.535 1st Qu.: 627.5 1st Qu.:0.91
## Median : 7.650 Median :1089.0 Median :1.12
## Mean : 7.186 Mean :1055.7 Mean :1.10
## 3rd Qu.: 8.500 3rd Qu.:1400.5 3rd Qu.:1.29
## Max. :10.520 Max. :2019.0 Max. :1.50
## Jumlah.Tenaga.Kerja.Usaha.Mikro.dan.Kecil Persentase.Penduduk.Miskin
## Min. : 7569 Min. : 2.380
## 1st Qu.: 31401 1st Qu.: 6.535
## Median : 50683 Median : 8.460
## Mean : 58482 Mean : 8.177
## 3rd Qu.: 79953 3rd Qu.:10.250
## Max. :205390 Max. :12.130
## IPM Rata.rata.Lama.Sekolah
## Min. :68.18 Min. : 6.940
## 1st Qu.:70.67 1st Qu.: 7.865
## Median :73.25 Median : 8.230
## Mean :74.11 Mean : 8.870
## 3rd Qu.:76.64 3rd Qu.: 9.970
## Max. :83.29 Max. :11.660y <- Jabar$Persentase.Penduduk.Miskin
x1 <- Jabar$Tingkat.Partisipasi.Angkatan.Kerja
x2 <- Jabar$Laju.Pertumbuhan.penduduk
x3 <- Jabar$Jumlah.Panjang.Jalan
x4 <- Jabar$Tingkat.Pengangguran.Terbuka
x5 <- Jabar$Jumlah.Tenaga.Kerja.Usaha.Mikro.dan.KecilLinear model OLS
ppm.ols <- lm(y ~ x1+x2+x3+x4+x5)
summary(ppm.ols)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5861 -1.7251 0.2177 1.9518 4.6014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.476e+00 1.516e+01 0.559 0.582
## x1 3.880e-02 1.951e-01 0.199 0.844
## x2 -8.597e-01 2.966e+00 -0.290 0.775
## x3 9.869e-05 1.143e-03 0.086 0.932
## x4 -2.412e-01 3.750e-01 -0.643 0.527
## x5 -5.529e-06 1.681e-05 -0.329 0.745
##
## Residual standard error: 2.867 on 21 degrees of freedom
## Multiple R-squared: 0.07217, Adjusted R-squared: -0.1487
## F-statistic: 0.3267 on 5 and 21 DF, p-value: 0.8912Assumption Checking (OLS)
Multicolinearity
vif(ppm.ols)## x1 x2 x3 x4 x5
## 1.584058 1.658095 1.486464 1.819071 1.568941Homoskedastisitas
bptest(ppm.ols)##
## studentized Breusch-Pagan test
##
## data: ppm.ols
## BP = 5.7273, df = 5, p-value = 0.3337Normality
e <- residuals(ppm.ols)
shapiro.test(e)##
## Shapiro-Wilk normality test
##
## data: e
## W = 0.98336, p-value = 0.929Non-autocorrelation
durbinWatsonTest(ppm.ols)## lag Autocorrelation D-W Statistic p-value
## 1 0.2770028 1.428695 0.07
## Alternative hypothesis: rho != 0Testing for Spatial Autocorrelation
Spatial weight matrix
list.queen <- poly2nb(Jabar, row.names=IDs, queen=TRUE)
W <- nb2listw(list.queen, style="W", zero.policy=TRUE)
WB <- nb2mat(list.queen, style='W', zero.policy = TRUE)Moran’s I test
lm.morantest(ppm.ols, W, alternative="two.sided")##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## weights: W
##
## Moran I statistic standard deviate = 3.0859, p-value = 0.002029
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## 0.35679086 -0.06162300 0.01838444moran.plot(y, W, labels=IDs, xlim=c(0, 100), ylim=c(0,100))
Spatial Econometrics
Lagrange multiplier test
lm.LMtests(ppm.ols, W, test="all")## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## test weights: listw
##
## RSerr = 5.5452, df = 1, p-value = 0.01853
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## test weights: listw
##
## RSlag = 7.9141, df = 1, p-value = 0.004905
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## test weights: listw
##
## adjRSerr = 5.5514, df = 1, p-value = 0.01847
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## test weights: listw
##
## adjRSlag = 7.9204, df = 1, p-value = 0.004888
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## test weights: listw
##
## SARMA = 13.466, df = 2, p-value = 0.001191Spatial error model (SEM)
ppm.sem <- errorsarlm(y ~ x1+x2+x3+x4+x5, listw = W)
summary(ppm.sem)##
## Call:errorsarlm(formula = y ~ x1 + x2 + x3 + x4 + x5, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.146453 -0.765939 0.065764 1.511321 2.855214
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1111e+00 9.2759e+00 -0.3354 0.7373
## x1 7.5749e-02 1.0860e-01 0.6975 0.4855
## x2 2.1354e+00 1.3795e+00 1.5480 0.1216
## x3 -2.2436e-04 5.9130e-04 -0.3794 0.7044
## x4 2.5629e-01 2.7881e-01 0.9192 0.3580
## x5 1.1192e-05 1.1002e-05 1.0173 0.3090
##
## Lambda: 0.8596, LR test value: 14.88, p-value: 0.00011457
## Asymptotic standard error: 0.075499
## z-value: 11.386, p-value: < 2.22e-16
## Wald statistic: 129.63, p-value: < 2.22e-16
##
## Log likelihood: -55.91984 for error model
## ML residual variance (sigma squared): 2.7633, (sigma: 1.6623)
## Number of observations: 27
## Number of parameters estimated: 8
## AIC: 127.84, (AIC for lm: 140.72)Spatial lag model (SAR)
ppm.sar <- lagsarlm(y ~ x1+x2+x3+x4+x5, listw = W)
summary(ppm.sar)##
## Call:lagsarlm(formula = y ~ x1 + x2 + x3 + x4 + x5, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.10438 -1.29987 0.20385 1.01461 3.28779
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5529e+00 9.6326e+00 -0.1612 0.8719
## x1 4.2344e-02 1.2335e-01 0.3433 0.7314
## x2 1.3190e+00 1.8754e+00 0.7033 0.4819
## x3 -1.0528e-05 7.2278e-04 -0.0146 0.9884
## x4 -1.3438e-01 2.3715e-01 -0.5666 0.5710
## x5 -1.0872e-06 1.0627e-05 -0.1023 0.9185
##
## Rho: 0.77274, LR test value: 12.359, p-value: 0.00043896
## Asymptotic standard error: 0.10852
## z-value: 7.1207, p-value: 1.0736e-12
## Wald statistic: 50.705, p-value: 1.0735e-12
##
## Log likelihood: -57.18058 for lag model
## ML residual variance (sigma squared): 3.2869, (sigma: 1.813)
## Number of observations: 27
## Number of parameters estimated: 8
## AIC: 130.36, (AIC for lm: 140.72)
## LM test for residual autocorrelation
## test value: 0.82214, p-value: 0.36456Spatial Lag of X (SLX)
ppm.slx <- lmSLX(y ~ x1+x2+x3+x4+x5, data=Jabar@data, W)
summary(ppm.slx)##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.323e+01 3.725e+01 -2.503e+00 2.353e-02
## x1 4.572e-01 2.125e-01 2.152e+00 4.704e-02
## x2 -1.605e+00 3.093e+00 -5.191e-01 6.108e-01
## x3 6.249e-04 1.105e-03 5.658e-01 5.794e-01
## x4 1.189e+00 4.602e-01 2.585e+00 1.995e-02
## x5 -1.110e-05 1.555e-05 -7.141e-01 4.855e-01
## lag.x1 1.157e+00 4.181e-01 2.768e+00 1.372e-02
## lag.x2 -7.988e+00 8.157e+00 -9.793e-01 3.420e-01
## lag.x3 2.320e-03 1.993e-03 1.164e+00 2.614e-01
## lag.x4 -6.242e-02 9.156e-01 -6.817e-02 9.465e-01
## lag.x5 -1.019e-04 3.556e-05 -2.867e+00 1.118e-02Backward Method
lag_x1 <- WB%*%x1
lag_x2 <- WB%*%x2
lag_x3 <- WB%*%x3
lag_x4 <- WB%*%x4
lag_x5 <- WB%*%x5
ppm.slx2 <- lm(y ~ x1 + x2 + x3 + x4 +x5 + lag_x1 + lag_x2 +
lag_x3 + lag_x4 + lag_x5)
summary(ppm.slx2)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + lag_x1 + lag_x2 + lag_x3 +
## lag_x4 + lag_x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2375 -1.2056 -0.2542 1.0565 4.8281
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.323e+01 3.725e+01 -2.503 0.0235 *
## x1 4.572e-01 2.125e-01 2.152 0.0470 *
## x2 -1.605e+00 3.093e+00 -0.519 0.6108
## x3 6.249e-04 1.105e-03 0.566 0.5794
## x4 1.189e+00 4.602e-01 2.585 0.0200 *
## x5 -1.110e-05 1.555e-05 -0.714 0.4855
## lag_x1 1.157e+00 4.181e-01 2.768 0.0137 *
## lag_x2 -7.988e+00 8.157e+00 -0.979 0.3420
## lag_x3 2.320e-03 1.993e-03 1.164 0.2614
## lag_x4 -6.242e-02 9.156e-01 -0.068 0.9465
## lag_x5 -1.019e-04 3.556e-05 -2.867 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.196 on 16 degrees of freedom
## Multiple R-squared: 0.5853, Adjusted R-squared: 0.326
## F-statistic: 2.258 on 10 and 16 DF, p-value: 0.07073ols_step_backward_p(ppm.slx2)##
##
## Stepwise Summary
## ------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## 0 Full Model 128.980 144.530 66.537 0.58525 0.32604
## 1 lag_x4 126.988 141.242 63.161 0.58513 0.36550
## 2 x3 125.532 138.491 59.816 0.57668 0.38854
## 3 x2 123.960 135.623 56.588 0.56991 0.41146
## 4 x5 122.464 132.831 53.489 0.56181 0.43035
## ------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.750 RMSE 1.738
## R-Squared 0.562 MSE 4.077
## Adj. R-Squared 0.430 Coef. Var 24.694
## Pred R-Squared 0.041 AIC 122.464
## MAE 1.429 SBC 132.831
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 104.543 6 17.424 4.274 0.0063
## Residual 81.539 20 4.077
## Total 186.083 26
## ------------------------------------------------------------------
##
## Parameter Estimates
## --------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## --------------------------------------------------------------------------------------------
## (Intercept) -91.532 29.813 -3.070 0.006 -153.721 -29.343
## x1 0.490 0.185 0.664 2.647 0.015 0.104 0.876
## x4 1.108 0.403 0.838 2.752 0.012 0.268 1.948
## lag_x1 1.089 0.342 0.627 3.185 0.005 0.376 1.803
## lag_x2 -8.128 2.930 -0.553 -2.774 0.012 -14.241 -2.016
## lag_x3 0.002 0.002 0.239 1.392 0.179 -0.001 0.005
## lag_x5 0.000 0.000 -0.895 -3.589 0.002 0.000 0.000
## --------------------------------------------------------------------------------------------SLX (from backward Method)
ppm.slx_d <- lm(y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 + lag_x5)
summary(ppm.slx_d)##
## Call:
## lm(formula = y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 + lag_x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3738 -1.1128 0.0602 1.3308 4.6921
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.153e+01 2.981e+01 -3.070 0.00604 **
## x1 4.899e-01 1.851e-01 2.647 0.01548 *
## x4 1.108e+00 4.026e-01 2.752 0.01228 *
## lag_x1 1.089e+00 3.421e-01 3.185 0.00466 **
## lag_x2 -8.128e+00 2.930e+00 -2.774 0.01172 *
## lag_x3 2.184e-03 1.569e-03 1.392 0.17927
## lag_x5 -1.107e-04 3.083e-05 -3.589 0.00183 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.019 on 20 degrees of freedom
## Multiple R-squared: 0.5618, Adjusted R-squared: 0.4304
## F-statistic: 4.274 on 6 and 20 DF, p-value: 0.006258Assumption Checking
Homoskedastisitas
bptest(ppm.slx_d)##
## studentized Breusch-Pagan test
##
## data: ppm.slx_d
## BP = 8.2627, df = 6, p-value = 0.2195Normality
shapiro.test(residuals(ppm.slx_d))##
## Shapiro-Wilk normality test
##
## data: residuals(ppm.slx_d)
## W = 0.97332, p-value = 0.6911Non-autocorrelation
durbinWatsonTest(ppm.slx_d)## lag Autocorrelation D-W Statistic p-value
## 1 -0.03288401 1.981851 0.954
## Alternative hypothesis: rho != 0Spatial Durbin Model (SDM)
ppm.sdm = lagsarlm(y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 + lag_x5, listw = W)
summary(ppm.sdm)##
## Call:lagsarlm(formula = y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 + lag_x5,
## listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.40754 -1.06819 -0.26891 0.72857 3.20356
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.7330e+01 2.0328e+01 -3.3122 0.0009257
## x1 3.3413e-01 1.2763e-01 2.6179 0.0088484
## x4 8.6808e-01 2.7348e-01 3.1742 0.0015026
## lag_x1 7.6064e-01 2.3439e-01 3.2452 0.0011737
## lag_x2 -6.3922e+00 2.0062e+00 -3.1862 0.0014415
## lag_x3 2.3993e-03 1.0596e-03 2.2643 0.0235548
## lag_x5 -7.7071e-05 2.1258e-05 -3.6255 0.0002884
##
## Rho: 0.63946, LR test value: 9.707, p-value: 0.0018357
## Asymptotic standard error: 0.13568
## z-value: 4.7129, p-value: 2.4419e-06
## Wald statistic: 22.212, p-value: 2.4419e-06
##
## Log likelihood: -48.37872 for lag model
## ML residual variance (sigma squared): 1.8573, (sigma: 1.3628)
## Number of observations: 27
## Number of parameters estimated: 9
## AIC: 114.76, (AIC for lm: 122.46)
## LM test for residual autocorrelation
## test value: 0.23235, p-value: 0.62979Assumption Checking
Heteroskedasticity
bptest.Sarlm(ppm.sdm)##
## studentized Breusch-Pagan test
##
## data:
## BP = 4.5405, df = 6, p-value = 0.6039Normality
shapiro.test(residuals(ppm.sdm))##
## Shapiro-Wilk normality test
##
## data: residuals(ppm.sdm)
## W = 0.96261, p-value = 0.423Spatial Durbin Error Model (SDEM)
ppm.sdem = errorsarlm(y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 + lag_x5,
data = data, listw = W)
summary(ppm.sdem)##
## Call:errorsarlm(formula = y ~ x1 + x4 + lag_x1 + lag_x2 + lag_x3 +
## lag_x5, data = data, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.36461 -1.04032 -0.15574 1.01793 3.61950
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.4400e+01 2.5455e+01 -2.5300 0.011408
## x1 3.6270e-01 1.4045e-01 2.5824 0.009811
## x4 8.2600e-01 3.0989e-01 2.6655 0.007688
## lag_x1 7.4303e-01 2.7586e-01 2.6935 0.007071
## lag_x2 -5.0139e+00 2.0684e+00 -2.4241 0.015348
## lag_x3 2.2532e-03 1.1695e-03 1.9267 0.054023
## lag_x5 -7.2977e-05 2.5490e-05 -2.8629 0.004197
##
## Lambda: 0.68894, LR test value: 4.9348, p-value: 0.026322
## Asymptotic standard error: 0.1353
## z-value: 5.092, p-value: 3.544e-07
## Wald statistic: 25.928, p-value: 3.544e-07
##
## Log likelihood: -50.76483 for error model
## ML residual variance (sigma squared): 2.1594, (sigma: 1.4695)
## Number of observations: 27
## Number of parameters estimated: 9
## AIC: 119.53, (AIC for lm: 122.46)Assumption Checking
Heteroskedasticity
bptest.Sarlm(ppm.sdem)##
## studentized Breusch-Pagan test
##
## data:
## BP = 7.5509, df = 6, p-value = 0.2729Normality
shapiro.test(residuals(ppm.sdem))##
## Shapiro-Wilk normality test
##
## data: residuals(ppm.sdem)
## W = 0.96022, p-value = 0.3739Model selection
data.frame("MODEL" = c("OLS", "SAR", "SEM", "SLX", "SDM", "SDEM"),
"AIC" = c(AIC(ppm.ols), AIC(ppm.sar), AIC(ppm.sem), AIC(ppm.slx_d), AIC(ppm.sdm), AIC(ppm.sdem)))## MODEL AIC
## 1 OLS 140.7198
## 2 SAR 130.3612
## 3 SEM 127.8397
## 4 SLX 122.4644
## 5 SDM 114.7574
## 6 SDEM 119.5297FOV test
er <- e
eu <- residuals(ppm.sdm)
n <- length(y)
k <- 2 #x1 dan x3
q <- 4 #lag_x1, lag_x2, lag_x3, lag_x5
num <- t(er)%*%er - t(eu)%*%eu
dem <- t(eu)%*%eu / (n-k)
fov <- num/dem ; fov## [,1]
## [1,] 61.07283F_crit = qf(p=.05, df1=q, df2=(n-k), lower.tail=FALSE) ; F_crit## [1] 2.75871fov > F_crit #Reject if TRUE## [,1]
## [1,] TRUEData visualisation
tm_shape(Jabar) +
tm_polygons(col = "Persentase.Penduduk.Miskin")
Jabar$Kab.Kota <- as.character(Kab.Kota)
tm_shape(Jabar) +
tm_fill("Persentase.Penduduk.Miskin",
title = "ppm(Quantile)",
style="quantile", palette = "Greens") +
tm_borders(alpha = 0.1) +
tm_text("Kab.Kota", size = 0.6, col = "black", shadow = T) +
tm_layout(main.title = "Persentase Penduduk Miskin Jawa Barat, 2023",
main.title.size = 1.2 ,
legend.position = c("right", "bottom"),
legend.title.size = 0.8,
legend.outside = TRUE)
Mapping residuals
Jabar$residuals <- residuals(ppm.sdm)
summary(Jabar$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.4075 -1.0682 -0.2689 0.0000 0.7286 3.2036my_breaks <- c(-Inf, -1.0682, -0.2689, 0.7286,Inf)
tm_shape(Jabar) +
tm_fill("residuals", title = "Residuals", style = "fixed",
breaks = my_breaks, palette = "-RdBu", midpoint = NA) + # Set midpoint = NA
tm_borders(alpha = 0.1) +
tm_layout(main.title = "Residuals", main.title.size = 0.7,
legend.position = c("right", "bottom"), legend.title.size = 0.8)