#Pemetaan Variabel Sanitasi
mapview(shp_joined, zcol = "Sanitasi")
#Pemetaan Variabel Jumlah Air Bersih
mapview(shp_joined, zcol = "Air")
# 2. Membuat matriks bobot spasial
# Queen contiguity
nb1 <- poly2nb(shp_joined, queen = TRUE)
lw1 <- nb2listw(nb1, style = "W", zero.policy = TRUE)
moranbiv_queen <- moranbi.test(shp_joined$Sanitasi,shp_joined$Air,lw1,rank=TRUE,
zero.policy=TRUE,alternative="two.sided")
moranbiv_queen
##
## Bivariate Moran I_{xy} test under randomisation
##
## data: shp_joined$Sanitasi using rank correction
## weights: lw1
##
## Bivariate Moran Z(I_{xy}) statistic = 2.1347, p-value = 0.03279
## alternative hypothesis: two.sided
## sample estimates:
## Bivariate Moran I_{xy} statistic Expectation
## 0.39866604 -0.13660229
## Variance
## 0.06287551
#Pembobot 1/jarak
centroid<-st_centroid(shp_joined)
longlat<-st_coordinates(centroid)
jarak<-dist(longlat)
w<-as.matrix(1/jarak)
lw<-mat2listw(w)
moranbiv_jarak <- moranbi.test(shp_joined$Sanitasi,shp_joined$Air,lw,rank=TRUE,
zero.policy=TRUE,alternative="two.sided")
moranbiv_jarak
##
## Bivariate Moran I_{xy} test under randomisation
##
## data: shp_joined$Sanitasi using rank correction
## weights: lw
##
## Bivariate Moran Z(I_{xy}) statistic = 0.63864, p-value = 0.5231
## alternative hypothesis: two.sided
## sample estimates:
## Bivariate Moran I_{xy} statistic Expectation
## 0.05087508 -0.13660229
## Variance
## 0.08617615
#PEMBOBOT KNN
longlat<- as.matrix(longlat)
k1<-knn2nb(knearneigh(longlat,k=2))
listw<-nb2listw(k1)
moranbiv_knn <- moranbi.test(shp_joined$Sanitasi,shp_joined$Air,listw,rank=TRUE,
zero.policy=TRUE,alternative="two.sided")
moranbiv_knn
##
## Bivariate Moran I_{xy} test under randomisation
##
## data: shp_joined$Sanitasi using rank correction
## weights: listw
##
## Bivariate Moran Z(I_{xy}) statistic = 2.3106, p-value = 0.02086
## alternative hypothesis: two.sided
## sample estimates:
## Bivariate Moran I_{xy} statistic Expectation
## 0.51659579 -0.13660229
## Variance
## 0.07991865
# 3. Global Bivariate Moran’s I
# Jumlah Sampah yang masuk TPA (Y) dengan Sanitasi (X)
moran_biv_as <- moranbi.test(shp_joined$Sanitasi,shp_joined$Air,lw1,rank=TRUE,
zero.policy=TRUE,alternative="two.sided")
moran_biv_as
##
## Bivariate Moran I_{xy} test under randomisation
##
## data: shp_joined$Sanitasi using rank correction
## weights: lw1
##
## Bivariate Moran Z(I_{xy}) statistic = 2.1347, p-value = 0.03279
## alternative hypothesis: two.sided
## sample estimates:
## Bivariate Moran I_{xy} statistic Expectation
## 0.39866604 -0.13660229
## Variance
## 0.06287551
# standarisasi peubah
scale.sanitasi <- scale(shp_joined$Sanitasi) %>% as.vector()
scale.air <- scale(shp_joined$Air) %>% as.vector()
# moran sccaterplot dasar
moran.plot(
x = scale.sanitasi,
listw = lw1, # Spatial weight
labels=as.character(shp_joined$WADMKK),
y = scale.air,
zero.policy = TRUE,
xlab = "Persentase Sanitasi Layak",
ylab = "Persentase Air Bersih",
main = "Moran Scatter Plot Bivariat: sanitasi layak Vs air bersih"
)
# Dark mode untuk poster
par(col.lab ="black", col.axis = "black", col.main = "black")
# 4. Local Bivariate Moran’s I (LISA)
local_biv_ss<- localmoran.bi(shp_joined$Sanitasi,shp_joined$Air,lw1,
zero.policy=TRUE,alternative="two.sided")
# Ubah ke data frame
local_biv_ss_df <- as.data.frame(local_biv_ss)
# Tambahkan nama daerah
local_biv_ss_df$Wilayah <- shp_joined$WADMKK
# Gabungkan hasil ke shapefile
shp_joined$Ii <- local_biv_ss[,1]
shp_joined$E.Ii <- local_biv_ss[,2]
shp_joined$Var.Ii <- local_biv_ss[,3]
shp_joined$Z.Ii <- local_biv_ss[,4]
shp_joined$Pvalue <- local_biv_ss[,5]
# Buat kolom signifikan / tidak
shp_joined$Signif <- ifelse(shp_joined$Pvalue < 0.05, "Signifikan", "Tidak Signifikan")
# Lihat tabel ringkas
# Hilangkan geometry
hasil_LISA <- st_drop_geometry(shp_joined[, c("WADMKK","Ii","Z.Ii","Pvalue","Signif")])
print(hasil_LISA)
## WADMKK Ii Z.Ii Pvalue Signif
## 1 Kota Cilegon 0.11565655 0.24897453 0.8033805 Tidak Signifikan
## 2 Kota Serang 0.06620612 0.19409708 0.8460999 Tidak Signifikan
## 3 Kota Tangerang 0.33427121 0.79307541 0.4277339 Tidak Signifikan
## 4 Kota Tangerang Selatan 0.40051754 0.83730845 0.4024192 Tidak Signifikan
## 5 Lebak 0.29209357 0.92370370 0.3556406 Tidak Signifikan
## 6 Pandeglang 1.01787488 1.64204863 0.1005799 Tidak Signifikan
## 7 Serang 0.02408981 0.76967347 0.4414936 Tidak Signifikan
## 8 Tangerang -0.01273981 -0.08908576 0.9290138 Tidak Signifikan
# Ubah ke data frame
hasil_lisa_df <- as.data.frame(hasil_LISA)
# Setting sumbu X dan Y pada kuadran
shp_joined$scale.Y <- scale(shp_joined$Air) %>% as.vector()
shp_joined$lag.X <- lag.listw(lw1, scale.sanitasi)
shp_joined$quad_sig <- NA
# high-high quadrant
shp_joined[(shp_joined$scale.Y >= 0 &
shp_joined$lag.X >= 0), "quad_sig"] <- "high-high"
# low-low quadrant
shp_joined[(shp_joined$scale.Y <= 0 &
shp_joined$lag.X <= 0), "quad_sig"] <- "low-low"
# high-low quadrant
shp_joined[(shp_joined$scale.Y >= 0 &
shp_joined$lag.X <= 0), "quad_sig"] <- "high-low"
# low-high quadrant
shp_joined[(shp_joined$scale.Y <= 0 &
shp_joined$lag.X >= 0), "quad_sig"] <- "low-high"
#Menghitung signifikansi
shp_joined$sig <- NA
# high-high quadrant
shp_joined[(local_biv_ss[,5] <= 0.05), "sig"] <- "Signifikan"
# low-low quadrant
shp_joined[(local_biv_ss[,5] > 0.05), "sig"] <- "Tidak Signifikan"
# 5. Visualisasi Hasil
unique(shp_joined$quad_sig)
table(shp_joined$quad_sig, useNA = "ifany")
#PEMETAAN MORAN SCATTERPLOT
#Pemetaan Kuadran Local Moran
tm_shape(shp_joined) +
tm_polygons(
fill = "quad_sig",
fill.scale = tm_scale(
values = c("high-high" = "royalblue4",
"low-low" = "indianred1",
"low-high" = "orchid4",
"high-low" = "hotpink2"),
labels = c("high-high", "low-low", "low-high", "high-low")
),
fill.legend = tm_legend(
title = "Kuadran LISA Sanitasi Layak Vs Air bersih"
)
) +
tm_layout(legend.outside = TRUE)
#PEMETAAN SIGNIFIKANSI AUTOKORELASI SPASIAL
tmap_mode("plot")
tm_shape(shp_joined) +
tm_polygons(
fill = "sig",
fill.scale = tm_scale(
values = c("Signifikan" = "steelblue4",
"Tidak Signifikan" = "red"),
labels = c("Signifikan", "Tidak Signifikan")
),
fill.legend = tm_legend(title = "Signifikansi")
) +
tm_borders(col = "black", lwd = 0.5, alpha = 0.5) +
tm_layout(frame = FALSE, legend.outside = TRUE) +
tm_title("Clusters")