library(readxl)
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
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, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(spdep)
## Loading required package: spData
## 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: '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.ZeroPolicyOption
library(CARBayes)
## Warning: package 'CARBayes' was built under R version 4.5.3
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: Rcpp
library(corrplot)
## corrplot 0.95 loaded
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
data_raw <- read_excel("D:/Persentase Penduduk Miskin Menurut Provinsi (3).xlsx", sheet = "Sheet2", skip = 1,
col_names = FALSE)
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
data <- data_raw[,1:8]
names(data) <- c(
"provinsi","kemiskinan","tpt","rls","ahh","ipm","pdrb","sanitasi")
# hapus baris kosong
data <- data %>% filter(!is.na(provinsi))
# fungsi cleaning numerik
clean_numeric <- function(x){
x <- as.character(x)
x[x == "-" | x == ""] <- NA
x <- gsub(",", ".", x)
x <- gsub("[^0-9.]", "", x)
as.numeric(x)
}
# konversi ke numerik
data$kemiskinan <- clean_numeric(data$kemiskinan)
data$tpt <- clean_numeric(data$tpt)
data$rls <- clean_numeric(data$rls)
data$ahh <- clean_numeric(data$ahh)
data$ipm <- clean_numeric(data$ipm)
data$pdrb <- clean_numeric(data$pdrb)
data$sanitasi <- clean_numeric(data$sanitasi)
# hapus NA
data <- na.omit(data)
# rapikan nama provinsi
data$provinsi <- toupper(trimws(data$provinsi))
summary(data)
## provinsi kemiskinan tpt rls
## Length:39 Min. : 3.100 Min. :1.180 Min. :63.98
## Class :character 1st Qu.: 4.725 1st Qu.:3.270 1st Qu.:67.11
## Mode :character Median : 6.760 Median :4.160 Median :69.64
## Mean : 7.045 Mean :4.380 Mean :69.92
## 3rd Qu.: 8.590 3rd Qu.:5.655 3rd Qu.:71.57
## Max. :16.560 Max. :7.020 Max. :87.53
## ahh ipm pdrb sanitasi
## Min. :66.68 Min. :53.42 Min. : 0.150 Min. :69.16
## 1st Qu.:70.95 1st Qu.:71.13 1st Qu.: 0.555 1st Qu.:82.31
## Median :72.95 Median :73.33 Median : 1.300 Median :86.31
## Mean :72.60 Mean :72.43 Mean : 3.401 Mean :86.05
## 3rd Qu.:74.31 3rd Qu.:74.31 3rd Qu.: 3.530 3rd Qu.:91.25
## Max. :77.40 Max. :83.08 Max. :26.560 Max. :97.90
map <- st_read("D:/gadm41_IDN_shp/gadm41_IDN_1.shp")
## Reading layer `gadm41_IDN_1' from data source `D:\gadm41_IDN_shp\gadm41_IDN_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84
map$NAME_1 <- toupper(trimws(map$NAME_1))
data_map <- left_join(map, data, by = c("NAME_1" = "provinsi"))
summary(data_map)
## GID_1 GID_0 COUNTRY NAME_1
## Length:34 Length:34 Length:34 Length:34
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## VARNAME_1 NL_NAME_1 TYPE_1 ENGTYPE_1
## Length:34 Length:34 Length:34 Length:34
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## CC_1 HASC_1 ISO_1 kemiskinan
## Length:34 Length:34 Length:34 Min. : 3.550
## Class :character Class :character Class :character 1st Qu.: 4.785
## Mode :character Mode :character Mode :character Median : 6.915
## Mean : 7.078
## 3rd Qu.: 8.600
## Max. :13.560
## NA's :4
## tpt rls ahh ipm
## Min. :1.870 Min. :63.98 Min. :68.27 Min. :67.02
## 1st Qu.:3.393 1st Qu.:66.50 1st Qu.:71.38 1st Qu.:71.56
## Median :4.140 Median :69.59 Median :72.70 Median :73.21
## Mean :4.383 Mean :69.38 Mean :72.77 Mean :72.84
## 3rd Qu.:5.445 3rd Qu.:71.51 3rd Qu.:74.32 3rd Qu.:74.08
## Max. :7.020 Max. :77.37 Max. :76.99 Max. :78.83
## NA's :4 NA's :4 NA's :4 NA's :4
## pdrb sanitasi geometry
## Min. : 0.2500 Min. :73.54 MULTIPOLYGON :34
## 1st Qu.: 0.6325 1st Qu.:82.19 epsg:4326 : 0
## Median : 1.3300 Median :86.31 +proj=long...: 0
## Mean : 2.6360 Mean :85.22
## 3rd Qu.: 3.1225 3rd Qu.:89.25
## Max. :14.3900 Max. :97.90
## NA's :4 NA's :4
data_map <- data_map %>%
filter(!is.na(kemiskinan),
!is.na(tpt),
!is.na(rls),
!is.na(ahh),
!is.na(ipm),
!is.na(pdrb),
!is.na(sanitasi))
data_nospasial <- st_drop_geometry(data_map)
cor_matrix <- cor(
data_nospasial[,c("kemiskinan","tpt","rls","ahh","ipm","pdrb","sanitasi")],
use = "complete.obs")
cor_matrix
## kemiskinan tpt rls ahh ipm
## kemiskinan 1.00000000 -0.22735299 0.288566729 -0.23933892 -0.3156050
## tpt -0.22735299 1.00000000 -0.580136168 0.20543209 0.3044948
## rls 0.28856673 -0.58013617 1.000000000 -0.28304055 -0.3483217
## ahh -0.23933892 0.20543209 -0.283040554 1.00000000 0.7652222
## ipm -0.31560504 0.30449478 -0.348321749 0.76522221 1.0000000
## pdrb 0.03471132 0.26937152 0.004337801 0.56156572 0.3499755
## sanitasi -0.29881657 -0.09988009 0.001967812 0.04667878 0.2800116
## pdrb sanitasi
## kemiskinan 0.034711318 -0.298816570
## tpt 0.269371524 -0.099880090
## rls 0.004337801 0.001967812
## ahh 0.561565725 0.046678783
## ipm 0.349975524 0.280011570
## pdrb 1.000000000 -0.092980531
## sanitasi -0.092980531 1.000000000
corrplot(cor_matrix, method = "color", type = "upper")
model_ols <- lm(
kemiskinan ~ tpt + rls + ahh + ipm + pdrb + sanitasi,
data = data_map
)
summary(model_ols)
##
## Call:
## lm(formula = kemiskinan ~ tpt + rls + ahh + ipm + pdrb + sanitasi,
## data = data_map)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0109 -1.6474 -0.1394 1.0541 4.9105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.259104 27.902417 1.335 0.195
## tpt -0.423969 0.518572 -0.818 0.422
## rls 0.062637 0.176742 0.354 0.726
## ahh -0.312088 0.392756 -0.795 0.435
## ipm -0.009302 0.328582 -0.028 0.978
## pdrb 0.163409 0.182309 0.896 0.379
## sanitasi -0.113977 0.081625 -1.396 0.176
##
## Residual standard error: 2.544 on 23 degrees of freedom
## Multiple R-squared: 0.2334, Adjusted R-squared: 0.0334
## F-statistic: 1.167 on 6 and 23 DF, p-value: 0.3576
library(lmtest)
# Heteroskedastisitas
bptest(model_ols)
##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 9.6659, df = 6, p-value = 0.1394
# Autokorelasi
dwtest(model_ols)
##
## Durbin-Watson test
##
## data: model_ols
## DW = 2.1803, p-value = 0.711
## alternative hypothesis: true autocorrelation is greater than 0
# Normalitas residual
shapiro.test(residuals(model_ols))
##
## Shapiro-Wilk normality test
##
## data: residuals(model_ols)
## W = 0.9619, p-value = 0.3461
ggplot(data_map) +
geom_sf(aes(fill = kemiskinan)) +
scale_fill_viridis_c() +
theme_minimal() +
labs(title = "Peta Persebaran Kemiskinan")
data_map$residual <- residuals(model_ols)
ggplot(data_map) +
geom_sf(aes(fill = residual)) +
scale_fill_viridis_c() +
theme_minimal() +
labs(title = "Peta Residual")
ggplot(data_map, aes(x = tpt, y = kemiskinan)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
ggplot(data_map, aes(x = kemiskinan)) +
geom_histogram(bins = 10) +
theme_minimal()
data_map <- data_map %>%
filter(!is.na(kemiskinan),
!is.na(tpt),
!is.na(rls),
!is.na(ahh),
!is.na(ipm),
!is.na(pdrb),
!is.na(sanitasi))
nb <- poly2nb(data_map)
## Warning in poly2nb(data_map): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(data_map): neighbour object has 9 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
# cek jumlah tetangga
card(nb)
## [1] 1 0 1 4 2 4 2 2 1 2 2 3 4 1 2 1 1 0 0 1 1 3 2 3 4 2 1 4 3 3
# hapus wilayah tanpa tetangga
data_map2 <- data_map[card(nb) > 0, ]
# buat neighbour baru
nb2 <- poly2nb(data_map2)
## Warning in poly2nb(data_map2): neighbour object has 6 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
# listw baru
lw2 <- nb2listw(nb2,
style = "W",
zero.policy = TRUE)
# matriks bobot baru
W2 <- nb2mat(nb2,
style = "B",
zero.policy = TRUE)
moran.test(
data_map2$kemiskinan,
lw2,
zero.policy = TRUE)
##
## Moran I test under randomisation
##
## data: data_map2$kemiskinan
## weights: lw2
##
## Moran I statistic standard deviate = 2.4286, p-value = 0.00758
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.42586839 -0.03846154 0.03655586
model_sar_bayes <- S.CARleroux(
formula = kemiskinan ~ tpt + rls + ahh + ipm + pdrb + sanitasi,family = "gaussian",
data = data_map2,
W = W2,
burnin = 20000,
n.sample = 100000)
## Setting up the model.
## Warning in mat2listw(W, style = "B"): neighbour object has 6 sub-graphs
##
## Markov chain 1 - generating 80000 post burnin and thinned samples.
## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
##
## Summarising results.
## Finished in 18.5 seconds.
summary(model_sar_bayes)
## Length Class Mode
## summary.results 70 -none- numeric
## samples 7 -none- list
## fitted.values 27 -none- numeric
## residuals 2 data.frame list
## modelfit 6 -none- numeric
## accept 5 -none- numeric
## localised.structure 0 -none- NULL
## formula 3 formula call
## model 2 -none- character
## mcmc.info 5 -none- numeric
## X 189 -none- numeric
# fitted value
data_map2$fitted_bayes <- model_sar_bayes$fitted.values
# residual
data_map2$residual_bayes <-
data_map2$kemiskinan - data_map2$fitted_bayes
ggplot(data_map2) +
geom_sf(aes(fill = fitted_bayes)) +
scale_fill_viridis_c() +
theme_minimal() +
labs(title = "Peta Fitted Value SAR Bayesian")
ggplot(data_map2) +
geom_sf(aes(fill = residual_bayes)) +
scale_fill_viridis_c() +
theme_minimal() +
labs(title = "Peta Residual SAR Bayesian")
ggplot(data_map2, aes(x = residual_bayes)) +
geom_histogram(bins = 10) +
theme_minimal() +
labs(title = "Histogram Residual SAR Bayesian")
