library(readr)
library(readxl)
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(spdep)
## Loading required package: spData
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(sp)
library(tmap)
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
data.q = read_excel("C:\\Users\\MUTHI'AH IFFA\\Downloads\\Data Kuis P2.xlsx")
str(data.q)
## tibble [119 × 8] (S3: tbl_df/tbl/data.frame)
## $ ID_KAB : num [1:119] 3101 3171 3172 3173 3174 ...
## $ PROPINSI: chr [1:119] "DKI JAKARTA" "DKI JAKARTA" "DKI JAKARTA" "DKI JAKARTA" ...
## $ KAB : chr [1:119] "KEPULAUAN SERIBU" "JAKARTA SELATAN" "JAKARTA TIMUR" "JAKARTA PUSAT" ...
## $ long : num [1:119] 106 107 107 107 107 ...
## $ lat : num [1:119] -5.8 -6.27 -6.26 -6.18 -6.17 ...
## $ y : num [1:119] 4.68 5.57 7.1 9.72 8.86 ...
## $ X1 : num [1:119] 0.586 0.294 0.814 0.724 0.949 ...
## $ X2 : num [1:119] 0.237 0.447 0.159 0.97 0.368 ...
data.sf = st_as_sf(data.q, coords = c("long","lat"), crs = 4326)
head(data.sf)
## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 106.4933 ymin: -6.272549 xmax: 106.9003 ymax: -5.796823
## Geodetic CRS: WGS 84
## # A tibble: 6 × 7
## ID_KAB PROPINSI KAB y X1 X2 geometry
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <POINT [°]>
## 1 3101 DKI JAKARTA KEPULAUAN SERI… 4.68 0.586 0.237 (106.4933 -5.796823)
## 2 3171 DKI JAKARTA JAKARTA SELATAN 5.57 0.294 0.447 (106.81 -6.272549)
## 3 3172 DKI JAKARTA JAKARTA TIMUR 7.10 0.814 0.159 (106.9003 -6.255373)
## 4 3173 DKI JAKARTA JAKARTA PUSAT 9.72 0.724 0.970 (106.8351 -6.18123)
## 5 3174 DKI JAKARTA JAKARTA BARAT 8.86 0.949 0.368 (106.7483 -6.165469)
## 6 3175 DKI JAKARTA JAKARTA UTARA 10.2 0.755 0.835 (106.8695 -6.129185)
str(data.sf)
## sf [119 × 7] (S3: sf/tbl_df/tbl/data.frame)
## $ ID_KAB : num [1:119] 3101 3171 3172 3173 3174 ...
## $ PROPINSI: chr [1:119] "DKI JAKARTA" "DKI JAKARTA" "DKI JAKARTA" "DKI JAKARTA" ...
## $ KAB : chr [1:119] "KEPULAUAN SERIBU" "JAKARTA SELATAN" "JAKARTA TIMUR" "JAKARTA PUSAT" ...
## $ y : num [1:119] 4.68 5.57 7.1 9.72 8.86 ...
## $ X1 : num [1:119] 0.586 0.294 0.814 0.724 0.949 ...
## $ X2 : num [1:119] 0.237 0.447 0.159 0.97 0.368 ...
## $ geometry:sfc_POINT of length 119; first list element: 'XY' num [1:2] 106.5 -5.8
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
## ..- attr(*, "names")= chr [1:6] "ID_KAB" "PROPINSI" "KAB" "y" ...
Eksplorasi Data
summary(data.sf$y)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.313 4.579 6.318 6.461 8.348 11.413
Histogram Distribusi
hist(data.sf$y, main="Distribusi Variabel y", col="skyblue", border="white")
tmap_mode("view")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(data.sf) +
tm_symbols(col="y", size=0.6, palette="YlOrRd", title.col="Nilai y") +
tm_layout(main.title="Sebaran Nilai y per Titik")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_symbols()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).[v3->v4] `symbols()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `symbols()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title.col' (rename to 'title') to 'fill.legend =
## tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.
data.sf = st_as_sf(data.q, coords = c("long", "lat"), crs = 4326)
data.sf_utm = st_transform(data.sf, crs = 32748)
coords = st_coordinates(data.sf_utm)
# Buat tetangga
knn5 = knearneigh(coords, k = 5)
nb_knn5 = knn2nb(knn5)
# Hitung jarak antar tetangga
dist_knn5 = nbdists(nb_knn5, coords)
# Ubah list
dist_vec = unlist(dist_knn5)
# Hitung radius
radius = median(dist_vec) * 3
radius
## [1] 103449.6
# Buat radial neighbors (0 ke radius)
nb_radial = dnearneigh(coords, 0, radius)
# Ubah ke matriks bobot spasial
lw_radial = nb2listw(nb_radial, style = "W")
# Cek hasil
summary(nb_radial)
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 2448
## Percentage nonzero weights: 17.28691
## Average number of links: 20.57143
## Link number distribution:
##
## 3 4 5 6 7 9 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 1 1 1 2 1 1 3 1 3 3 3 12 9 7 7 12 6 8 5 9 6 7 3 2 2 3
## 32
## 1
## 1 least connected region:
## 83 with 3 links
## 1 most connected region:
## 47 with 32 links
Eksplorasi Data Spasial
# Statistik deskriptif
summary(data.sf_utm$y)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.313 4.579 6.318 6.461 8.348 11.413
summary(data.sf_utm$X1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03614 0.28226 0.50627 0.52169 0.76510 0.99011
summary(data.sf_utm$X2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.007716 0.217327 0.445661 0.479496 0.740893 0.980693
# Buang geometry untuk eksplorasi numerik
data.num = st_drop_geometry(data.sf_utm)
# Ringkasan statistik
summary(data.num[, c("y", "X1", "X2")])
## y X1 X2
## Min. : 1.313 Min. :0.03614 Min. :0.007716
## 1st Qu.: 4.579 1st Qu.:0.28226 1st Qu.:0.217327
## Median : 6.318 Median :0.50627 Median :0.445661
## Mean : 6.461 Mean :0.52169 Mean :0.479496
## 3rd Qu.: 8.348 3rd Qu.:0.76510 3rd Qu.:0.740893
## Max. :11.413 Max. :0.99011 Max. :0.980693
# Korelasi antar variabel
cor(data.num[, c("y", "X1", "X2")])
## y X1 X2
## y 1.0000000 0.8505082 0.2961602
## X1 0.8505082 1.0000000 -0.0570887
## X2 0.2961602 -0.0570887 1.0000000
# Plot sebaran
pairs(data.num[, c("y", "X1", "X2")], main = "Plot Sebaran Antar Variabel")
# Uji korelasi signifikan (opsional)
cor.test(data.num$y, data.num$X1)
##
## Pearson's product-moment correlation
##
## data: data.num$y and data.num$X1
## t = 17.491, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7917144 0.8936907
## sample estimates:
## cor
## 0.8505082
cor.test(data.num$y, data.num$X2)
##
## Pearson's product-moment correlation
##
## data: data.num$y and data.num$X2
## t = 3.3539, df = 117, p-value = 0.001074
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1227058 0.4520577
## sample estimates:
## cor
## 0.2961602
Moran’s I
# Uji autokorelasi spasial global
moran_y = moran.test(data.num$y, lw_radial)
moran_y
##
## Moran I test under randomisation
##
## data: data.num$y
## weights: lw_radial
##
## Moran I statistic standard deviate = 3.3593, p-value = 0.0003907
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.0858215390 -0.0084745763 0.0007879286
moran.plot(data.num$y, lw_radial, main = "Moran Scatterplot untuk y")
Local Moran’s I
moran.plot(data.num$y, lw_radial, main = "Moran Scatterplot untuk y")
# Hitung LISA
lisa <- localmoran(data.num$y, lw_radial)
# Tambahkan hasil ke data spasial
data.sf_utm$Ii = lisa[,1]
data.sf_utm$Ei = lisa[,2]
data.sf_utm$Vi = lisa[,3]
data.sf_utm$Z_Ii = lisa[,4]
data.sf_utm$p_value = lisa[,5]
# Klasifikasi klaster LISA (opsional)
data.sf_utm$cluster_lisa <- factor(ifelse(data.sf_utm$p_value < 0.05,
ifelse(data.sf_utm$Z_Ii > 0, "High-High", "Low-Low"),
"Non-Signif"))
Uji Autokorelasi LISA
# Local Moran's I
local_moran = localmoran(data.num$y, lw_radial)
# Gabungkan hasil ke data spasial
data.sf_utm$Ii = local_moran[, 1]
data.sf_utm$P_value = local_moran[, 5]
# Visualisasi hasil LISA
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
tm_shape(data.sf_utm) +
tm_fill("Ii", palette = "RdBu", style = "pretty", title = "Local Moran's I") +
tm_borders()
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: instead of `style = "pretty"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
## [scale] tm_fill:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
# Statistik deskriptif dan korelasi
data.num <- st_drop_geometry(data.sf_utm)
summary(data.num[, c("y", "X1", "X2")])
## y X1 X2
## Min. : 1.313 Min. :0.03614 Min. :0.007716
## 1st Qu.: 4.579 1st Qu.:0.28226 1st Qu.:0.217327
## Median : 6.318 Median :0.50627 Median :0.445661
## Mean : 6.461 Mean :0.52169 Mean :0.479496
## 3rd Qu.: 8.348 3rd Qu.:0.76510 3rd Qu.:0.740893
## Max. :11.413 Max. :0.99011 Max. :0.980693
cor(data.num[, c("y", "X1", "X2")])
## y X1 X2
## y 1.0000000 0.8505082 0.2961602
## X1 0.8505082 1.0000000 -0.0570887
## X2 0.2961602 -0.0570887 1.0000000
# Uji Moran Global
moran_y <- moran.test(data.num$y, lw_radial)
print(moran_y)
##
## Moran I test under randomisation
##
## data: data.num$y
## weights: lw_radial
##
## Moran I statistic standard deviate = 3.3593, p-value = 0.0003907
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.0858215390 -0.0084745763 0.0007879286
# Uji LISA (Local Moran's I)
local_moran <- localmoran(data.num$y, lw_radial)
# Gabungkan hasil ke data spasial
data.sf_utm$Ii <- local_moran[, 1]
data.sf_utm$E.Ii <- local_moran[, 2]
data.sf_utm$Var.Ii <- local_moran[, 3]
data.sf_utm$Z.Ii <- local_moran[, 4]
data.sf_utm$p.value <- local_moran[, 5]
# Klasifikasi klaster LISA
data.sf_utm$cluster_lisa <- factor(
ifelse(data.sf_utm$p.value < 0.05,
ifelse(data.sf_utm$Z.Ii > 0, "High-High", "Low-Low"),
"Non-Signif"),
levels = c("High-High", "Low-Low", "Non-Signif")
)
# Cek hasil klasifikasi
table(data.sf_utm$cluster_lisa)
##
## High-High Low-Low Non-Signif
## 18 11 90
# Peta Local Moran's I
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
tm_shape(data.sf_utm) +
tm_fill(
"Ii",
fill.scale = tm_scale_intervals(values = "brewer.rd_bu"),
fill.legend = tm_legend(title = "Local Moran's I")
) +
tm_borders() +
tm_layout(title = "Peta Local Moran’s I untuk Variabel y")
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [scale] tm_fill:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
packageVersion("tmap")
## [1] '4.2'
# --- 8. Peta Klaster LISA (versi tmap 3.x)
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
tm_shape(data.sf_utm) +
tm_fill(
"cluster_lisa",
palette = c("red", "blue", "grey"),
title = "Klaster LISA"
) +
tm_borders() +
tm_layout(
title = "Peta Klaster LISA (High-High, Low-Low, Non-Signif)",
legend.outside = TRUE
)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
