Spatial Lag Model
Ardiana
2024-12-11
Import Packages
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')`## Loading required package: sf## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(sp)
library(sf)
library(ggplot2)
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(spdep)
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.ZeroPolicyOptionlibrary(RColorBrewer)Import Data
setwd("D:/DIAN/unpad/studi/SEMESTER 5/Spatial Statistics/Tugas/uas")
# Data Jawa Timur
miskin <- readxl::read_xlsx("data uas spasial.xlsx")
miskin$id = c(1:38)
miskin## # A tibble: 38 × 6
## `Kab/Kota` `miskin %` pop `pdb rate` `Penempatan Tenaga Kerja` id
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Bangkalan 19.4 1092. 0.18 1.98 1
## 2 Banyuwangi 7.34 1744. 4.38 2.25 2
## 3 Batu 3.31 220. 4.94 5.40 3
## 4 Blitar 8.69 1254. 3.58 2.32 4
## 5 Bojonegoro 12.2 1320. 2 1.33 5
## 6 Bondowoso 13.3 788. 4.05 0.947 6
## 7 Gresik 11.0 1350. 3.54 11.0 7
## 8 Jember 9.51 2587. 4.21 0.676 8
## 9 Jombang 9.15 1351. 4.13 6.01 9
## 10 Kediri 10.7 1677. 3.71 5.83 10
## # ℹ 28 more rows# Data Persentase Penduduk Miskin Pulau Jawa
miskin_jawa = readxl::read_xlsx("Persentase Penduduk Miskin (P0) Menurut Provinsi dan Daerah, 2023.xlsx")
miskin_jawa## # A tibble: 6 × 2
## prov `miskin %`
## <chr> <dbl>
## 1 BANTEN 6.17
## 2 Jakarta raya 4.44
## 3 JAWA BARAT 7.62
## 4 JAWA TIMUR 10.4
## 5 JAWA TENGAH 10.8
## 6 Yogyakarta 11.0miskin_jawa$id = c(1:6)Indo<-st_read("gadm41_IDN.gpkg", layer = "ADM_ADM_1") #sf package, reading .gpkg file## Reading layer `ADM_ADM_1' from data source
## `D:\DIAN\unpad\studi\SEMESTER 5\Spatial Statistics\Tugas\uas\gadm41_IDN.gpkg'
## using driver `GPKG'
## 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 84head(Indo)## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -8.849262 xmax: 123.5523 ymax: 6.076941
## Geodetic CRS: WGS 84
## GID_1 GID_0 COUNTRY NAME_1 VARNAME_1 NL_NAME_1 TYPE_1
## 1 IDN.1_1 IDN Indonesia Aceh NA NA Propinisi
## 2 IDN.2_1 IDN Indonesia Bali NA NA Propinisi
## 3 IDN.3_1 IDN Indonesia Bangka Belitung NA NA Propinisi
## 4 IDN.4_1 IDN Indonesia Banten NA NA Propinisi
## 5 IDN.5_1 IDN Indonesia Bengkulu NA NA Propinisi
## 6 IDN.6_1 IDN Indonesia Gorontalo NA NA Propinisi
## ENGTYPE_1 CC_1 HASC_1 ISO_1 geom
## 1 Province 11 ID.AC ID-AC MULTIPOLYGON (((98.14903 2....
## 2 Province 51 ID.BA ID-BA MULTIPOLYGON (((115.5257 -8...
## 3 Province 19 ID.BB NA MULTIPOLYGON (((107.9614 -3...
## 4 Province 36 ID.BT ID-BT MULTIPOLYGON (((106.3868 -6...
## 5 Province 17 ID.BE ID-BE MULTIPOLYGON (((103.573 -4....
## 6 Province 75 ID.GO ID-GO MULTIPOLYGON (((123.5491 0....jawa<-Indo[c(4,7,9,10,11,34),]
jawa## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS: WGS 84
## GID_1 GID_0 COUNTRY NAME_1 VARNAME_1 NL_NAME_1 TYPE_1
## 4 IDN.4_1 IDN Indonesia Banten NA NA Propinisi
## 7 IDN.7_1 IDN Indonesia Jakarta Raya NA NA Propinisi
## 9 IDN.9_1 IDN Indonesia Jawa Barat NA NA Propinisi
## 10 IDN.10_1 IDN Indonesia Jawa Tengah NA NA Propinisi
## 11 IDN.11_1 IDN Indonesia Jawa Timur NA NA Propinisi
## 34 IDN.33_1 IDN Indonesia Yogyakarta NA NA Propinisi
## ENGTYPE_1 CC_1 HASC_1 ISO_1 geom
## 4 Province 36 ID.BT ID-BT MULTIPOLYGON (((106.3868 -6...
## 7 Province 31 ID.JK ID-JK MULTIPOLYGON (((106.8672 -6...
## 9 Province 32 ID.JR ID-JB MULTIPOLYGON (((108.1331 -7...
## 10 Province 33 ID.JT ID-JT MULTIPOLYGON (((110.0044 -7...
## 11 Province 35 ID.JI ID-JI MULTIPOLYGON (((114.0286 -8...
## 34 Province 34 ID.YO ID-YO MULTIPOLYGON (((110.7128 -8...jawa$NAME_1## [1] "Banten" "Jakarta Raya" "Jawa Barat" "Jawa Tengah" "Jawa Timur"
## [6] "Yogyakarta"jawa$id = c(1:6)
indo_sf<-st_as_sf(Indo)
jawa_sf = st_as_sf(jawa)
####
Indo_KabKot<-readRDS('gadm36_IDN_2_sp.rds')
jatim<-Indo_KabKot[Indo_KabKot$NAME_1 == "Jawa Timur",]
jatim$NAME_2## [1] "Bangkalan" "Banyuwangi" "Batu" "Blitar"
## [5] "Bojonegoro" "Bondowoso" "Gresik" "Jember"
## [9] "Jombang" "Kediri" "Kota Blitar" "Kota Kediri"
## [13] "Kota Madiun" "Kota Malang" "Kota Mojokerto" "Kota Pasuruan"
## [17] "Kota Probolinggo" "Lamongan" "Lumajang" "Madiun"
## [21] "Magetan" "Malang" "Mojokerto" "Nganjuk"
## [25] "Ngawi" "Pacitan" "Pamekasan" "Pasuruan"
## [29] "Ponorogo" "Probolinggo" "Sampang" "Sidoarjo"
## [33] "Situbondo" "Sumenep" "Surabaya" "Trenggalek"
## [37] "Tuban" "Tulungagung"plot(jatim)
Analisis Deskriptif
Ringkasan
summary(miskin)## Kab/Kota miskin % pop pdb rate
## Length:38 Min. : 3.310 Min. : 136.1 Min. :0.180
## Class :character 1st Qu.: 7.188 1st Qu.: 708.0 1st Qu.:3.785
## Mode :character Median : 9.665 Median :1099.8 Median :4.230
## Mean :10.293 Mean :1092.8 Mean :3.915
## 3rd Qu.:12.360 3rd Qu.:1342.7 3rd Qu.:4.577
## Max. :21.760 Max. :2911.4 Max. :5.280
## Penempatan Tenaga Kerja id
## Min. : 0.676 Min. : 1.00
## 1st Qu.: 1.535 1st Qu.:10.25
## Median : 1.986 Median :19.50
## Mean : 4.068 Mean :19.50
## 3rd Qu.: 5.152 3rd Qu.:28.75
## Max. :24.269 Max. :38.00data.frame(sd(miskin$`miskin %`),sd(miskin$pop),sd(miskin$`pdb rate`),sd(miskin$`Penempatan Tenaga Kerja`))## sd.miskin..miskin.... sd.miskin.pop. sd.miskin..pdb.rate..
## 1 4.321289 682.4675 1.164824
## sd.miskin..Penempatan.Tenaga.Kerja..
## 1 4.805769Korelasi
cor = cor(miskin[,-c(1,6)]); cor## miskin % pop pdb rate
## miskin % 1.00000000 0.03683401 -0.4426950
## pop 0.03683401 1.00000000 0.1622698
## pdb rate -0.44269505 0.16226985 1.0000000
## Penempatan Tenaga Kerja -0.35831560 0.56189991 0.3133745
## Penempatan Tenaga Kerja
## miskin % -0.3583156
## pop 0.5618999
## pdb rate 0.3133745
## Penempatan Tenaga Kerja 1.0000000heatmap(cor,
col = colorRampPalette(c("blue", "white", "red"))(256),
cexRow = 1,
cexCol = 1)
Bar chart
data_sorted <- miskin[order(miskin$`miskin %`, decreasing = TRUE), ]
# Create a barplot
barplot(data_sorted$`miskin %`,
names.arg = data_sorted$`Kab/Kota`,
col = "lightblue",
main = "Persentase Penduduk Miskin di Jawa Timur 2023",
xlab = "Kabupaten/Kota",
ylab = "Persentase Penduduk Miskin %",
cex.names = 0.5,
las = 2) 
Peta Persebaran
Pulau Jawa
jawa_merged <- jawa_sf %>% left_join(miskin_jawa, by = "id")
head(jawa_merged)## Simple feature collection with 6 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS: WGS 84
## GID_1 GID_0 COUNTRY NAME_1 VARNAME_1 NL_NAME_1 TYPE_1 ENGTYPE_1
## 1 IDN.4_1 IDN Indonesia Banten NA NA Propinisi Province
## 2 IDN.7_1 IDN Indonesia Jakarta Raya NA NA Propinisi Province
## 3 IDN.9_1 IDN Indonesia Jawa Barat NA NA Propinisi Province
## 4 IDN.10_1 IDN Indonesia Jawa Tengah NA NA Propinisi Province
## 5 IDN.11_1 IDN Indonesia Jawa Timur NA NA Propinisi Province
## 6 IDN.33_1 IDN Indonesia Yogyakarta NA NA Propinisi Province
## CC_1 HASC_1 ISO_1 id prov miskin % geom
## 1 36 ID.BT ID-BT 1 BANTEN 6.17 MULTIPOLYGON (((106.3868 -6...
## 2 31 ID.JK ID-JK 2 Jakarta raya 4.44 MULTIPOLYGON (((106.8672 -6...
## 3 32 ID.JR ID-JB 3 JAWA BARAT 7.62 MULTIPOLYGON (((108.1331 -7...
## 4 33 ID.JT ID-JT 4 JAWA TIMUR 10.35 MULTIPOLYGON (((110.0044 -7...
## 5 35 ID.JI ID-JI 5 JAWA TENGAH 10.77 MULTIPOLYGON (((114.0286 -8...
## 6 34 ID.YO ID-YO 6 Yogyakarta 11.04 MULTIPOLYGON (((110.7128 -8...breaks <- c(0, 6.17, 4.44, 7.62, 10.35, 10.77, 11.04)
labels <- c("6.17", "4.44", "7.62", "10.35", "10.77", "11.04")
jawa_merged$`miskin_bin` <- cut(jawa_merged$`miskin %`,
breaks = breaks,
labels = labels,
right = TRUE)
ggplot() +
geom_sf(data=jawa_merged, aes(fill = jawa_merged$miskin_bin),color="white") +
theme_bw() +
scale_fill_manual(values = c("6.17" = "lightyellow",
"4.44" = "yellow",
"7.62" = "orange",
"10.35" = "red",
"10.77" = "red3",
"11.04"="red4"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank()) +
theme(legend.position = "right",
legend.direction = "vertical",
legend.key.width = unit(0.5,"cm"),
legend.key.height = unit(0.3,"cm"),
legend.text = element_text(size = 6),
plot.title = element_text(hjust = 0.5, size = 10),
legend.title = element_text(size=8),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank()) +
labs(title = "Persebaran Persentase Penduduk Miskin di Pulau Jawa 2023",fill = "% Penduduk Miskin") +
geom_sf_text(data = jawa_merged, aes(label = jawa_merged$NAME_1), size = 2.5, color = "black")## Warning: Use of `jawa_merged$miskin_bin` is discouraged.
## ℹ Use `miskin_bin` instead.## Warning: Use of `jawa_merged$NAME_1` is discouraged.
## ℹ Use `NAME_1` instead.## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Jawa Timur
jatim$id<-c(1:38)
jatim_sf<-st_as_sf(jatim)
jatim_merged <- jatim_sf %>%
left_join(miskin, by = "id")
head(jatim_merged)## Simple feature collection with 6 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 111.4236 ymin: -8.78036 xmax: 114.605 ymax: -6.879476
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## GID_0 NAME_0 GID_1 NAME_1 NL_NAME_1 GID_2 NAME_2 VARNAME_2
## 1 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.1_1 Bangkalan <NA>
## 2 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.2_1 Banyuwangi <NA>
## 3 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.3_1 Batu <NA>
## 4 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.4_1 Blitar <NA>
## 5 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.5_1 Bojonegoro <NA>
## 6 IDN Indonesia IDN.11_1 Jawa Timur <NA> IDN.11.6_1 Bondowoso <NA>
## NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2 id Kab/Kota miskin % pop
## 1 <NA> Kabupaten Regency 3526 ID.JI.BK 1 Bangkalan 19.35 1091.8
## 2 <NA> Kabupaten Regency 3510 ID.JI.BW 2 Banyuwangi 7.34 1743.9
## 3 <NA> Kota City 3579 ID.JI.BA 3 Batu 3.31 220.2
## 4 <NA> Kabupaten Regency 3505 ID.JI.BR 4 Blitar 8.69 1253.6
## 5 <NA> Kabupaten Regency 3522 ID.JI.BJ 5 Bojonegoro 12.18 1319.6
## 6 <NA> Kabupaten Regency 3511 ID.JI.BD 6 Bondowoso 13.34 788.2
## pdb rate Penempatan Tenaga Kerja geometry
## 1 0.18 1.977 MULTIPOLYGON (((112.695 -7....
## 2 4.38 2.253 MULTIPOLYGON (((114.0247 -8...
## 3 4.94 5.402 MULTIPOLYGON (((112.5227 -7...
## 4 3.58 2.317 MULTIPOLYGON (((112.237 -8....
## 5 2.00 1.332 MULTIPOLYGON (((111.5608 -7...
## 6 4.05 0.947 MULTIPOLYGON (((113.7905 -8...ggplot() +
geom_sf(data=jatim_merged, aes(fill = jatim_merged$`miskin %`),color="grey",size=0.1) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank()) +
theme(legend.position = "right",
legend.direction = "vertical",
legend.key.width = unit(0.3,"cm"),
legend.key.height = unit(1,"cm"),
legend.text = element_text(size = 6),
plot.title = element_text(hjust = 0.5, size = 10),
legend.title = element_text(size=8),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank()) +
labs(title = "Persebaran Persentase Kemiskinan di Jawa Timur 2023",fill = "% Penduduk Miskin") +
geom_sf_text(data = jatim_merged, aes(label = id), size = 2, color = "black")## Warning: Use of `` jatim_merged$`miskin %` `` is discouraged.
## ℹ Use `miskin %` instead.## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
jatim_merged$NAME_2## [1] "Bangkalan" "Banyuwangi" "Batu" "Blitar"
## [5] "Bojonegoro" "Bondowoso" "Gresik" "Jember"
## [9] "Jombang" "Kediri" "Kota Blitar" "Kota Kediri"
## [13] "Kota Madiun" "Kota Malang" "Kota Mojokerto" "Kota Pasuruan"
## [17] "Kota Probolinggo" "Lamongan" "Lumajang" "Madiun"
## [21] "Magetan" "Malang" "Mojokerto" "Nganjuk"
## [25] "Ngawi" "Pacitan" "Pamekasan" "Pasuruan"
## [29] "Ponorogo" "Probolinggo" "Sampang" "Sidoarjo"
## [33] "Situbondo" "Sumenep" "Surabaya" "Trenggalek"
## [37] "Tuban" "Tulungagung"Jawa Timur (Quantile)
summary(miskin$`miskin %`)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.310 7.188 9.665 10.293 12.360 21.760breaks <- c(-Inf, 7.2, 9.7, 12.4,Inf)
# Define labels for each interval
labels <- c("Very Low", "Low", "High", "Very High")
# Create a new column with discretized Diare
jatim_merged$miskin_discrete <- cut(jatim_merged$`miskin %`, breaks = breaks, labels = labels, right = TRUE)
ggplot() +
geom_sf(data=jatim_merged, aes(fill = miskin_discrete),color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow",
"Low" = "orange",
"High" = "red",
"Very High" = "red3"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank()) +
theme(legend.position = "right",
legend.direction = "vertical",
legend.key.width = unit(0.5,"cm"),
legend.key.height = unit(0.3,"cm"),
legend.text = element_text(size = 6),
plot.title = element_text(hjust = 0.5, size = 10),
legend.title = element_text(size=8),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank()) +
labs(title = "Persebaran Persentase Penduduk Miskin di Pulau Jawa 2023",fill = "% Penduduk Miskin") +
geom_sf_text(data = jatim_merged, aes(label = id), size = 2.5, color = "black")## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Autokorelasi Spatial
Matriks Bobot Spasial
#Membuat matriks bobot
row.names(jatim) <- as.character(1:38)
W <- poly2nb(jatim, row.names(jatim), queen=TRUE) ## Warning in poly2nb(jatim, row.names(jatim), queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.summary(W)## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 138
## Percentage nonzero weights: 9.556787
## Average number of links: 3.631579
## 2 disjoint connected subgraphs
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9
## 8 6 6 6 2 7 1 1 1
## 8 least connected regions:
## 1 11 12 14 15 16 17 34 with 1 link
## 1 most connected region:
## 22 with 9 linksis.symmetric.nb(W)## [1] TRUEWB <- nb2mat(W, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matriks biner "B"
WB[0:10, 0:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## 1 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 1 0 1 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 1
## 5 0 0 0 0 0 0 0 0 1 0
## 6 0 1 0 0 0 0 0 1 0 0
## 7 0 0 0 0 0 0 0 0 0 0
## 8 0 1 0 0 0 1 0 0 0 0
## 9 0 0 0 0 1 0 0 0 0 1
## 10 0 0 0 1 0 0 0 0 1 0WB## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 1 0 1 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 1 1 0 0 0
## 5 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 6 0 1 0 0 0 0 0 1 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 8 0 1 0 0 0 1 0 0 0 0 0 0 0 0
## 9 0 0 0 0 1 0 0 0 0 1 0 0 0 0
## 10 0 0 0 1 0 0 0 0 1 0 0 1 0 0
## 11 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## 12 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 17 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 18 0 0 0 0 1 0 1 0 1 0 0 0 0 0
## 19 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## 20 0 0 0 0 1 0 0 0 0 0 0 0 1 0
## 21 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## 22 0 0 1 1 0 0 0 0 1 1 0 0 0 1
## 23 0 0 1 0 0 0 1 0 1 0 0 0 0 0
## 24 0 0 0 0 1 0 0 0 1 1 0 0 0 0
## 25 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 26 0 0 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 0 0
## 28 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## 29 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 30 0 0 0 0 0 1 0 1 0 0 0 0 0 0
## 31 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 32 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## 33 0 1 0 0 0 1 0 0 0 0 0 0 0 0
## 34 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 35 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## 36 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 37 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 38 0 0 0 1 0 0 0 0 0 1 0 0 0 0
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## 1 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 1 1 0 0 0
## 4 0 0 0 0 0 0 0 1 0 0 0 0
## 5 0 0 0 1 0 1 0 0 0 1 1 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0
## 7 0 0 0 1 0 0 0 0 1 0 0 0
## 8 0 0 0 0 1 0 0 0 0 0 0 0
## 9 0 0 0 1 0 0 0 1 1 1 0 0
## 10 0 0 0 0 0 0 0 1 0 1 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 0 0 0 0 0 1 1 0 0 0 0 0
## 14 0 0 0 0 0 0 0 1 0 0 0 0
## 15 0 0 0 0 0 0 0 0 1 0 0 0
## 16 0 0 0 0 0 0 0 0 0 0 0 0
## 17 0 0 0 0 0 0 0 0 0 0 0 0
## 18 0 0 0 0 0 0 0 0 1 0 0 0
## 19 0 0 0 0 0 0 0 1 0 0 0 0
## 20 0 0 0 0 0 0 1 0 0 1 1 0
## 21 0 0 0 0 0 1 0 0 0 0 1 0
## 22 0 0 0 0 1 0 0 0 1 0 0 0
## 23 1 0 0 1 0 0 0 1 0 0 0 0
## 24 0 0 0 0 0 1 0 0 0 0 0 0
## 25 0 0 0 0 0 1 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
## 28 0 1 0 0 0 0 0 1 1 0 0 0
## 29 0 0 0 0 0 1 1 0 0 1 0 1
## 30 0 0 1 0 1 0 0 1 0 0 0 0
## 31 0 0 0 0 0 0 0 0 0 0 0 0
## 32 0 0 0 0 0 0 0 0 1 0 0 0
## 33 0 0 0 0 0 0 0 0 0 0 0 0
## 34 0 0 0 0 0 0 0 0 0 0 0 0
## 35 0 0 0 0 0 0 0 0 0 0 0 0
## 36 0 0 0 0 0 0 0 0 0 0 0 1
## 37 0 0 0 1 0 0 0 0 0 0 0 0
## 38 0 0 0 0 0 0 0 0 0 1 0 0
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## 1 0 0 0 0 1 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 1 0 0 0 0 0
## 3 0 1 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 1
## 5 0 0 0 0 0 0 0 0 0 0 1 0
## 6 0 0 0 1 0 0 1 0 0 0 0 0
## 7 0 0 0 0 0 1 0 0 1 0 0 0
## 8 0 0 0 1 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0 0 0 0
## 10 0 0 0 0 0 0 0 0 0 0 0 1
## 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 0 0 0 0 0 0 0 0 0 0 0 0
## 14 0 0 0 0 0 0 0 0 0 0 0 0
## 15 0 0 0 0 0 0 0 0 0 0 0 0
## 16 0 1 0 0 0 0 0 0 0 0 0 0
## 17 0 0 0 1 0 0 0 0 0 0 0 0
## 18 0 0 0 0 0 0 0 0 0 0 1 0
## 19 0 0 0 1 0 0 0 0 0 0 0 0
## 20 0 0 1 0 0 0 0 0 0 0 0 0
## 21 0 0 1 0 0 0 0 0 0 0 0 0
## 22 0 1 0 1 0 0 0 0 0 0 0 0
## 23 0 1 0 0 0 1 0 0 0 0 0 0
## 24 0 0 1 0 0 0 0 0 0 0 0 1
## 25 0 0 0 0 0 0 0 0 0 0 0 0
## 26 0 0 1 0 0 0 0 0 0 1 0 0
## 27 0 0 0 0 1 0 0 1 0 0 0 0
## 28 0 0 0 1 0 1 0 0 0 0 0 0
## 29 0 0 0 0 0 0 0 0 0 1 0 1
## 30 0 1 0 0 0 0 1 0 0 0 0 0
## 31 1 0 0 0 0 0 0 0 0 0 0 0
## 32 0 1 0 0 0 0 0 0 1 0 0 0
## 33 0 0 0 1 0 0 0 0 0 0 0 0
## 34 1 0 0 0 0 0 0 0 0 0 0 0
## 35 0 0 0 0 0 1 0 0 0 0 0 0
## 36 0 0 1 0 0 0 0 0 0 0 0 1
## 37 0 0 0 0 0 0 0 0 0 0 0 0
## 38 0 0 1 0 0 0 0 0 0 1 0 0
## attr(,"call")
## nb2mat(neighbours = W, style = "B", zero.policy = TRUE)WL<-nb2listw(W)#List neighbours
WL## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 138
## Percentage nonzero weights: 9.556787
## Average number of links: 3.631579
## 2 disjoint connected subgraphs
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 27.00053 169.545CoordK<-coordinates(jatim)
CoordK## [,1] [,2]
## 1 112.9302 -7.043867
## 2 114.2054 -8.364818
## 3 112.5311 -7.833109
## 4 112.2368 -8.129433
## 5 111.8098 -7.255433
## 6 113.9476 -7.944023
## 7 112.5386 -7.146270
## 8 113.6635 -8.230025
## 9 112.2651 -7.545102
## 10 112.0863 -7.828760
## 11 112.1670 -8.095024
## 12 112.0138 -7.826305
## 13 111.5284 -7.627492
## 14 112.6364 -7.978964
## 15 112.4374 -7.471312
## 16 112.9097 -7.649963
## 17 113.2053 -7.775074
## 18 112.3008 -7.131243
## 19 113.1376 -8.125301
## 20 111.6457 -7.624239
## 21 111.3579 -7.662881
## 22 112.6407 -8.121300
## 23 112.4843 -7.548270
## 24 111.9384 -7.597423
## 25 111.3431 -7.439238
## 26 111.1791 -8.126201
## 27 113.5039 -7.064970
## 28 112.8316 -7.743180
## 29 111.4994 -7.931320
## 30 113.3210 -7.866758
## 31 113.2561 -7.051974
## 32 112.7003 -7.451632
## 33 114.0524 -7.801377
## 34 113.8205 -6.990251
## 35 112.7227 -7.275410
## 36 111.6263 -8.161720
## 37 111.8915 -6.953156
## 38 111.8872 -8.113254plot(jatim, axes=T, col="gray90")+
text(CoordK[,1], CoordK[,2], row.names(jatim), col="black", cex=0.8, pos=1.5)+
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")+
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
## integer(0)jatim$NAME_2## [1] "Bangkalan" "Banyuwangi" "Batu" "Blitar"
## [5] "Bojonegoro" "Bondowoso" "Gresik" "Jember"
## [9] "Jombang" "Kediri" "Kota Blitar" "Kota Kediri"
## [13] "Kota Madiun" "Kota Malang" "Kota Mojokerto" "Kota Pasuruan"
## [17] "Kota Probolinggo" "Lamongan" "Lumajang" "Madiun"
## [21] "Magetan" "Malang" "Mojokerto" "Nganjuk"
## [25] "Ngawi" "Pacitan" "Pamekasan" "Pasuruan"
## [29] "Ponorogo" "Probolinggo" "Sampang" "Sidoarjo"
## [33] "Situbondo" "Sumenep" "Surabaya" "Trenggalek"
## [37] "Tuban" "Tulungagung"GloNAME_2## Global Moran’s I
Global_Moran<-moran.test(miskin$`miskin %`,WL)
Global_Moran##
## Moran I test under randomisation
##
## data: miskin$`miskin %`
## weights: WL
##
## Moran I statistic standard deviate = 3.2783, p-value = 0.0005221
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.39778983 -0.02702703 0.01679181Local Moran’s I
Local_Moran <- localmoran(miskin$`miskin %`, WL)
quadr_data <- attr(Local_Moran, "quadr")
mean_values <- quadr_data$mean
jatim_sf <- st_as_sf(jatim)
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
jatim_sf$mean_values <- mean_values_char
ggplot() +
geom_sf(data = jatim_sf, aes(fill = mean_values)) + # Fill based on 'mean_values'
scale_fill_manual(values = c("Low-Low" = "blue", "High-Low" = "green",
"Low-High" = "yellow", "High-High" = "red")) +
labs(title = "Local Moran's I Mean Values in East Java", fill = "Mean Type") +
theme_minimal()
Dataframe
Y = miskin$`miskin %`
X1 = miskin$`pop`
X2 = miskin$`pdb rate`
X3 = miskin$`Penempatan Tenaga Kerja`
df = data.frame(Y,X1,X2,X3)Uji Multikolinieritas
vif=diag(solve(cor(df[,-1])))
vif## X1 X2 X3
## 1.461867 1.109241 1.578376OLS Regression
miskin.ols<-lm(Y~., data=df);summary(miskin.ols)##
## Call:
## lm(formula = Y ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4933 -2.9288 -0.1003 2.5381 8.0818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.776841 2.298066 6.430 2.39e-07 ***
## X1 0.002158 0.001076 2.006 0.0529 .
## X2 -1.339870 0.548940 -2.441 0.0200 *
## X3 -0.392600 0.158714 -2.474 0.0185 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.693 on 34 degrees of freedom
## Multiple R-squared: 0.3289, Adjusted R-squared: 0.2697
## F-statistic: 5.554 on 3 and 34 DF, p-value: 0.003266plot(miskin.ols)
# Uji Asumsi Model OLS ## Uji Asumsi: Normalitas
shapiro.test(miskin.ols$residuals)##
## Shapiro-Wilk normality test
##
## data: miskin.ols$residuals
## W = 0.9813, p-value = 0.7632Uji Asumsi: Homoskedastisitas
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(miskin.ols)##
## studentized Breusch-Pagan test
##
## data: miskin.ols
## BP = 7.6416, df = 3, p-value = 0.05403Uji Asums: Autokorelasi
moran.test(miskin.ols$residuals, WL)##
## Moran I test under randomisation
##
## data: miskin.ols$residuals
## weights: WL
##
## Moran I statistic standard deviate = 1.3512, p-value = 0.08831
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.14971527 -0.02702703 0.01710883LM Test
LM<-lm.RStests(miskin.ols, WL, test="all")
print(LM)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ ., data = df)
## test weights: WL
##
## RSerr = 1.1987, df = 1, p-value = 0.2736
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ ., data = df)
## test weights: WL
##
## RSlag = 3.6262, df = 1, p-value = 0.05688
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ ., data = df)
## test weights: WL
##
## adjRSerr = 2.1286, df = 1, p-value = 0.1446
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ ., data = df)
## test weights: WL
##
## adjRSlag = 4.5561, df = 1, p-value = 0.0328
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ ., data = df)
## test weights: WL
##
## SARMA = 5.7548, df = 2, p-value = 0.05628Spatial lag model (SAR)
sar.miskin<-lagsarlm(Y~., data=df, WL)
summary(sar.miskin)##
## Call:lagsarlm(formula = Y ~ ., data = df, listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.26105 -2.45834 -0.13609 2.36927 7.15162
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.53862758 2.70583747 3.5252 0.0004232
## X1 0.00201508 0.00095418 2.1118 0.0347011
## X2 -1.02043095 0.48651750 -2.0974 0.0359565
## X3 -0.31398043 0.14428409 -2.1761 0.0295458
##
## Rho: 0.3439, LR test value: 4.0222, p-value: 0.044905
## Asymptotic standard error: 0.15364
## z-value: 2.2383, p-value: 0.025198
## Wald statistic: 5.0102, p-value: 0.025198
##
## Log likelihood: -99.43941 for lag model
## ML residual variance (sigma squared): 10.593, (sigma: 3.2547)
## Number of observations: 38
## Number of parameters estimated: 6
## AIC: 210.88, (AIC for lm: 212.9)
## LM test for residual autocorrelation
## test value: 4.0379, p-value: 0.044489Koefisien Model SAR
koef <- impacts(sar.miskin, listw = WL)
print(koef)## Impact measures (lag, exact):
## Direct Indirect Total
## X1 0.002091061 0.000980213 0.003071274
## X2 -1.058910328 -0.496378372 -1.555288700
## X3 -0.325820303 -0.152732623 -0.478552925Uji Asumsi Model SAR
Uji Asumsi: Homoskedastisitas
bptest.Sarlm(sar.miskin)##
## studentized Breusch-Pagan test
##
## data:
## BP = 3.8375, df = 3, p-value = 0.2796Uji Asumsi: Normalitas Residu
ResSAR <- residuals(sar.miskin)
shapiro.test(ResSAR)##
## Shapiro-Wilk normality test
##
## data: ResSAR
## W = 0.97011, p-value = 0.3947Uji Asumsi: NonAutokorelasi
moran.test(ResSAR, WL)##
## Moran I test under randomisation
##
## data: ResSAR
## weights: WL
##
## Moran I statistic standard deviate = -0.69248, p-value = 0.7557
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.11777566 -0.02702703 0.01717374Evaluasi Metrik
AIC
data.frame(Model = c("OLS","SAR"), AIC = c(AIC(miskin.ols),AIC(sar.miskin)))## Model AIC
## 1 OLS 212.9010
## 2 SAR 210.8788Koefisien Determinasi
summary(miskin.ols)##
## Call:
## lm(formula = Y ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4933 -2.9288 -0.1003 2.5381 8.0818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.776841 2.298066 6.430 2.39e-07 ***
## X1 0.002158 0.001076 2.006 0.0529 .
## X2 -1.339870 0.548940 -2.441 0.0200 *
## X3 -0.392600 0.158714 -2.474 0.0185 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.693 on 34 degrees of freedom
## Multiple R-squared: 0.3289, Adjusted R-squared: 0.2697
## F-statistic: 5.554 on 3 and 34 DF, p-value: 0.003266cat("R-squared Model SAR:",1-(sar.miskin$SSE/(var(df$Y)*(length(df$Y)-1))),"\n",
"R-squared Model OLS:",0.2697)## R-squared Model SAR: 0.4173767
## R-squared Model OLS: 0.2697Peta Residual
jatim_merged$resid <- ResSAR
library(ggplot2)
ggplot(data = jatim_merged) +
geom_sf(aes(fill = resid)) +
scale_fill_viridis_c(option = "plasma") +
labs(title = "Residuals Map", fill = "resid") +
theme_minimal()