UTS SPASIAL

Load Library yang Dibutuhkan

# KODE R: ANALISIS SPASIAL IPM & KEMISKINAN BANTEN



library(sf)

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

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(ggplot2)
library(readxl)
library(tmap)

## Registered S3 method overwritten by 'stars':
##   method                  from
##   st_interpolate_aw.stars sf

Memuat data dan shp

# Membaca shapefile
banten_sf <- st_read("C:\\Users\\user\\Downloads\\UTS2\\2-SHPBANTEN\\Banten_ADMIN_BPS.shp")  

## Reading layer `Banten_ADMIN_BPS' from data source 
##   `C:\Users\user\Downloads\UTS2\2-SHPBANTEN\Banten_ADMIN_BPS.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -7.016779 xmax: 106.7799 ymax: -5.807418
## Geodetic CRS:  WGS 84

data_excel <- read_excel("C:\\Users\\user\\Downloads\\UTS2\\DataSpasial.xlsx")

# Samakan format teks 
banten_sf$Kabupaten_clean <- trimws(toupper(banten_sf$Kabupaten))
data_excel$Kabupaten_clean <- trimws(toupper(data_excel$`Kabupaten/Kota`))

# Lakukan penggabungan data (Left Join)
banten_merged <- banten_sf %>%
  left_join(data_excel, by = c("Kabupaten_clean" = "Kabupaten_clean"))

# Cek kembali apakah masih ada data yang NA (Missing) -> Target hasil harus 0
sum(is.na(banten_merged$`IPM 2025 (X)`)) 

## [1] 0

Statistik Deskriptif

summary(banten_merged$`IPM 2025 (X)`)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   69.24   73.53   77.12   77.02   80.26   84.81

summary(banten_merged$`Penduduk miskin Maret 2025 (Y, %)`)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.390   4.220   5.350   5.496   6.822   8.510

# Menghitung Standar Deviasi
sd_ipm <- sd(banten_merged$`IPM 2025 (X)`, na.rm = TRUE)
sd_miskin <- sd(banten_merged$`Penduduk miskin Maret 2025 (Y, %)`, na.rm = TRUE)

cat("Standar Deviasi IPM:", sd_ipm, "\n")

## Standar Deviasi IPM: 5.241082

cat("Standar Deviasi Penduduk Miskin:", sd_miskin, "\n")

## Standar Deviasi Penduduk Miskin: 2.116587

Berdasarkan delapan kabupaten/kota di Provinsi Banten, IPM memiliki rata-rata sebesar 77,02, dengan nilai minimum 69,24 dan maksimum 84,81. Sementara itu, persentase penduduk miskin memiliki rata-rata 5,496%, dengan nilai minimum 2,39% dan maksimum 8,51%.

Hasil tersebut menunjukkan adanya variasi kondisi pembangunan manusia dan kemiskinan antarwilayah. Kota Tangerang Selatan memiliki IPM tertinggi sekaligus persentase penduduk miskin terendah, sedangkan Kabupaten Lebak memiliki IPM terendah dan Kabupaten Pandeglang memiliki persentase penduduk miskin tertinggi.

Visualisasi Peta Tematik

tmap_mode("plot")

## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`

# Peta Tematik IPM
p_ipm <- tm_shape(banten_merged) +
  tm_polygons("IPM 2025 (X)", palette = "Blues", title = "IPM 2025") +
  tm_layout(main.title = "Peta Sebaran IPM Provinsi Banten", main.title.size = 1)

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

# Peta Tematik Penduduk Miskin
p_misk <- tm_shape(banten_merged) +
  tm_polygons("Penduduk miskin Maret 2025 (Y, %)", palette = "Reds", title = "Penduduk Miskin (%)") +
  tm_layout(main.title = "Peta Sebaran Kemiskinan Banten", main.title.size = 1)

## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

# Menampilkan peta berdampingan
tmap_arrange(p_ipm, p_misk)

## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"
## Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.
## 
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"
## Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

Gradasi warna biru menunjukkan perbedaan nilai IPM antarwilayah. Wilayah dengan warna lebih gelap memiliki IPM lebih tinggi, sedangkan warna lebih terang menunjukkan IPM lebih rendah. Kota Tangerang Selatan memiliki IPM tertinggi sebesar 84,81, sementara Kabupaten Lebak terendah sebesar 69,24.

Gradasi warna merah menunjukkan persentase penduduk miskin. Warna yang semakin gelap menunjukkan persentase yang semakin tinggi. Kabupaten Pandeglang memiliki persentase tertinggi sebesar 8,51%, sedangkan Kota Tangerang Selatan terendah sebesar 2,39%.

Peta tematik memperlihatkan bahwa tingkat pembangunan manusia dan kemiskinan bervariasi antarwilayah. Namun, pola warna pada peta saja belum cukup untuk membuktikan adanya autokorelasi spasial.

Analisis Aoutokorelasi spasial

# ANALISIS KORELASI DAN AUTOKORELASI SPASIAL
# IPM DAN PERSENTASE PENDUDUK MISKIN PROVINSI BANTEN


# Memastikan data lengkap dan urutannya sesuai
variabel <- c(
  "IPM 2025 (X)",
  "Penduduk miskin Maret 2025 (Y, %)"
)

data_analisis <- banten_merged

# Menghapus wilayah yang memiliki nilai hilang pada kedua variabel
data_analisis <- data_analisis[
  complete.cases(st_drop_geometry(data_analisis)[, variabel]),
]

x <- as.numeric(data_analisis[["IPM 2025 (X)"]])
y <- as.numeric(data_analisis[["Penduduk miskin Maret 2025 (Y, %)"]])

n <- length(x)

if (n < 3) {
  stop("Jumlah wilayah dengan data lengkap terlalu sedikit.")
}

if (sd(x) == 0 || sd(y) == 0) {
  stop("Variabel tidak memiliki variasi sehingga analisis tidak dapat dilakukan.")
}

cat("Jumlah wilayah yang dianalisis:", n, "\n")

## Jumlah wilayah yang dianalisis: 8

# A. KORELASI PEARSON

korelasi_pearson <- cor.test(x, y, method = "pearson")

cat("\nHASIL KORELASI PEARSON\n")

## 
## HASIL KORELASI PEARSON

print(korelasi_pearson)

## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = -3.5032, df = 6, p-value = 0.01278
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9662139 -0.2718522
## sample estimates:
##        cor 
## -0.8195303

# Interpretasi:
# Korelasi Pearson mengukur hubungan linear IPM dan kemiskinan
# tanpa memperhitungkan posisi atau ketetanggaan wilayah.

Hasil korelasi Pearson menunjukkan \(r=-0{,}8195\) dengan p-value sebesar 0,01278.

Interpretasi: Terdapat hubungan linear negatif yang kuat dan signifikan antara IPM dan persentase penduduk miskin pada wilayah yang sama. Artinya, wilayah dengan IPM lebih tinggi cenderung memiliki persentase penduduk miskin lebih rendah. Namun, korelasi ini tidak mengukur hubungan dengan wilayah tetangga.

# B. MENYIAPKAN MATRIKS PEMBOBOT SPASIAL

# Mengubah koordinat ke sistem proyeksi dalam satuan meter
banten_utm <- st_transform(data_analisis, 32748)

# B1. Queen Contiguity
# Dua wilayah bertetangga jika berbagi sisi atau titik sudut
nb_queen <- poly2nb(banten_utm, queen = TRUE)

lw_queen <- nb2listw(
  nb_queen,
  style = "W",
  zero.policy = TRUE
)

# B2. K-Nearest Neighbors (k = 3)
# Tetangga ditentukan berdasarkan tiga centroid terdekat
centroid <- st_centroid(st_geometry(banten_utm))
koordinat <- st_coordinates(centroid)

nb_knn <- knn2nb(
  knearneigh(koordinat, k = 3)
)

## Warning in knearneigh(koordinat, k = 3): k greater than one-third of the number
## of data points

# Membuat hubungan ketetanggaan simetris
nb_knn <- make.sym.nb(nb_knn)

lw_knn <- nb2listw(
  nb_knn,
  style = "W",
  zero.policy = TRUE
)

# B3. Inverse Distance
# Wilayah yang lebih dekat mendapatkan bobot lebih besar
jarak <- as.matrix(dist(koordinat))

bobot_idw <- 1 / pmax(jarak, 1e-10)
diag(bobot_idw) <- 0

lw_idw <- mat2listw(
  bobot_idw,
  style = "W",
  zero.policy = TRUE
)

# Memeriksa jumlah tetangga setiap wilayah
cat("\nJUMLAH TETANGGA QUEEN\n")

## 
## JUMLAH TETANGGA QUEEN

print(card(nb_queen))

## [1] 1 1 2 2 3 2 5 4

cat("\nJUMLAH TETANGGA KNN\n")

## 
## JUMLAH TETANGGA KNN

print(card(nb_knn))

## [1] 3 7 3 3 3 3 5 5

Metode Hasil jumlah tetangga Interpretasi Queen Contiguity 1–5 tetangga Ketetanggaan ditentukan berdasarkan persinggungan sisi atau titik sudut. KNN (\(k=3\)) 3–7 tetangga Berdasarkan tiga centroid terdekat, kemudian hubungan dibuat simetris. Inverse Distance Semua wilayah diberi bobot berdasarkan jarak Wilayah yang lebih dekat memiliki bobot lebih besar.

Ketiga metode membentuk struktur hubungan spasial yang berbeda. Struktur tersebut digunakan untuk menganalisis hubungan IPM dengan persentase penduduk miskin di wilayah tetangga.

# C. GLOBAL BIVARIATE MORAN'S I

# Standardisasi kedua variabel
zx <- as.numeric(scale(x))
zy <- as.numeric(scale(y))

# Fungsi menghitung Global Bivariate Moran's I
hitung_bivariat <- function(zx, zy, lw) {
  lag_zy <- lag.listw(
    lw, zy,
    zero.policy = TRUE
  )

  sum(zx * lag_zy) / sum(zx^2)
}

# Uji signifikansi menggunakan permutasi
uji_bivariat <- function(zx, zy, lw, nsim = 9999) {

  I_obs <- hitung_bivariat(zx, zy, lw)

  set.seed(123)
  I_sim <- replicate(nsim, {
    zy_acak <- sample(zy)
    hitung_bivariat(zx, zy_acak, lw)
  })

  # P-value dua sisi berdasarkan distribusi permutasi
  p_value <- (
    sum(abs(I_sim - mean(I_sim)) >=
          abs(I_obs - mean(I_sim))) + 1
  ) / (nsim + 1)

  data.frame(
    Moran_Bivariat = I_obs,
    Rata_rata_Permutasi = mean(I_sim),
    P_value = p_value,
    Signifikansi = ifelse(
      p_value < 0.05,
      "Signifikan",
      "Tidak signifikan"
    )
  )
}

# Menguji tiga matriks pembobot
hasil_global <- bind_rows(
  cbind(
    Bobot = "Inverse Distance",
    uji_bivariat(zx, zy, lw_idw)
  ),
  cbind(
    Bobot = "K-Nearest Neighbors (k=3)",
    uji_bivariat(zx, zy, lw_knn)
  ),
  cbind(
    Bobot = "Queen Contiguity",
    uji_bivariat(zx, zy, lw_queen)
  )
)

cat("\nHASIL GLOBAL BIVARIATE MORAN'S I\n")

## 
## HASIL GLOBAL BIVARIATE MORAN'S I

print(hasil_global)

##                       Bobot Moran_Bivariat Rata_rata_Permutasi P_value
## 1          Inverse Distance     0.01272686         0.001131906  0.9126
## 2 K-Nearest Neighbors (k=3)    -0.15979983         0.002435167  0.4674
## 3          Queen Contiguity    -0.26378551         0.002199056  0.2243
##       Signifikansi
## 1 Tidak signifikan
## 2 Tidak signifikan
## 3 Tidak signifikan

Matriks pembobot Moran’s I P-value Interpretasi Inverse Distance 0,0127 0,9126 Hubungan positif sangat lemah dan tidak signifikan. KNN (\(k=3\)) -0,1598 0,4674 Kecenderungan hubungan negatif, tidak signifikan. Queen Contiguity -0,2638 0,2243 Kecenderungan hubungan negatif, tidak signifikan.

Interpretasi: Ketiga metode menghasilkan p-value lebih besar dari 0,05. Dengan demikian, belum terdapat bukti statistik yang cukup untuk menyatakan bahwa IPM suatu wilayah berkaitan dengan persentase penduduk miskin di wilayah tetangganya. Nilai Moran’s I yang negatif menunjukkan kecenderungan hubungan berlawanan arah, tetapi kecenderungan tersebut tidak signifikan.

# D. GLOBAL MORAN'S I UNIVARIAT

# Analisis univariat tetap dilakukan sebagai analisis tambahan
cat("\nGLOBAL MORAN'S I UNIVARIAT: IPM\n")

## 
## GLOBAL MORAN'S I UNIVARIAT: IPM

print(moran.test(
  x, lw_queen,
  zero.policy = TRUE,
  randomisation = TRUE
))

## 
##  Moran I test under randomisation
## 
## data:  x  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 2.2867, p-value = 0.01111
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.45067712       -0.14285714        0.06736823

cat("\nGLOBAL MORAN'S I UNIVARIAT: KEMISKINAN\n")

## 
## GLOBAL MORAN'S I UNIVARIAT: KEMISKINAN

print(moran.test(
  y, lw_queen,
  zero.policy = TRUE,
  randomisation = TRUE
))

## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 1.1454, p-value = 0.126
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.15519921       -0.14285714        0.06771927

# E. LOCAL BIVARIATE MORAN (LISA)
# Menggunakan Queen Contiguity sebagai bobot utama


lw_best <- lw_queen

lag_zy <- lag.listw(
  lw_best, zy,
  zero.policy = TRUE
)

# Statistik lokal: IPM wilayah i x rata-rata tertimbang
# kemiskinan di wilayah tetangganya
I_lokal <- zx * lag_zy

# Menghitung p-value lokal melalui permutasi
set.seed(123)
nsim_local <- 9999
p_lokal <- rep(NA_real_, n)

for (i in seq_len(n)) {

  tetangga_i <- lw_best$neighbours[[i]]
  bobot_i <- lw_best$weights[[i]]

  # Wilayah tanpa tetangga tidak dapat diuji
  if (length(tetangga_i) == 0) next

  # Statistik lokal hasil observasi
  I_obs <- I_lokal[i]

  # Permutasi nilai Y pada lingkungan lokal
  I_sim <- replicate(nsim_local, {
    y_acak <- sample(zy, size = length(tetangga_i),
                     replace = FALSE)
    zx[i] * sum(bobot_i * y_acak)
  })

  p_lokal[i] <- (
    sum(abs(I_sim - mean(I_sim)) >=
          abs(I_obs - mean(I_sim))) + 1
  ) / (nsim_local + 1)
}

# Koreksi pengujian berganda menggunakan Benjamini-Hochberg
p_adj <- p.adjust(p_lokal, method = "BH")

# Menentukan kategori pola berdasarkan tanda nilai
# IPM wilayah dan lag kemiskinan tetangganya
kuadran <- case_when(
  zx >= 0 & lag_zy >= 0 ~ "High-High",
  zx <  0 & lag_zy <  0 ~ "Low-Low",
  zx <  0 & lag_zy >= 0 ~ "Low-High",
  zx >= 0 & lag_zy <  0 ~ "High-Low"
)

# Hanya menyebut pola sebagai signifikan jika p-value
# setelah koreksi lebih kecil dari 0,05
kategori_lisa <- ifelse(
  !is.na(p_adj) & p_adj < 0.05,
  kuadran,
  "Tidak signifikan"
)

data_analisis$Moran_Bivariat_Lokal <- I_lokal
data_analisis$P_value_Lokal <- p_lokal
data_analisis$P_value_BH <- p_adj
data_analisis$Kuadran <- kuadran
data_analisis$Kategori_LISA <- factor(
  kategori_lisa,
  levels = c(
    "High-High", "Low-Low",
    "High-Low", "Low-High",
    "Tidak signifikan"
  )
)

# Tabel hasil LISA
hasil_lisa <- st_drop_geometry(data_analisis) %>%
  select(
    Kabupaten,
    `IPM 2025 (X)`,
    `Penduduk miskin Maret 2025 (Y, %)`,
    Moran_Bivariat_Lokal,
    P_value_Lokal,
    P_value_BH,
    Kategori_LISA
  )

cat("\nHASIL LOCAL BIVARIATE MORAN (LISA)\n")

## 
## HASIL LOCAL BIVARIATE MORAN (LISA)

print(hasil_lisa)

##                Kabupaten IPM 2025 (X) Penduduk miskin Maret 2025 (Y, %)
## 1           Kota Cilegon        79.54                              3.44
## 2            Kota Serang        77.50                              5.51
## 3         Kota Tangerang        82.41                              5.19
## 4 Kota Tangerang Selatan        84.81                              2.39
## 5                  Lebak        69.24                              8.03
## 6             Pandeglang        71.84                              8.51
## 7                 Serang        74.10                              4.48
## 8              Tangerang        76.74                              6.42
##   Moran_Bivariat_Lokal P_value_Lokal P_value_BH    Kategori_LISA
## 1          -0.23062850        0.6233  0.7123429 Tidak signifikan
## 2          -0.04374384        0.6172  0.7123429 Tidak signifikan
## 3          -0.52997389        0.4986  0.7123429 Tidak signifikan
## 4           0.21674444        0.8219  0.8219000 Tidak signifikan
## 5          -0.68313971        0.3554  0.7123429 Tidak signifikan
## 6          -0.35447119        0.6025  0.7123429 Tidak signifikan
## 7          -0.23335043        0.1433  0.7123429 Tidak signifikan
## 8           0.01206454        0.5718  0.7123429 Tidak signifikan

Dari delapan wilayah yang dianalisis, Kabupaten Lebak memiliki nilai Moran bivariat lokal paling negatif, yaitu -0,6831. Sementara itu, Kota Tangerang Selatan memiliki nilai positif tertinggi, yaitu 0,2167. Akan tetapi, seluruh nilai P-value BH berada di atas 0,05 sehingga tidak ada wilayah yang signifikan setelah koreksi Benjamini–Hochberg.

# F. PETA LISA BIVARIAT


tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

peta_lisa <- tm_shape(data_analisis) +
  tm_polygons(
    "Kategori_LISA",
    palette = c(
      "red", "blue", "pink", "lightblue", "grey80"
    ),
    title = "Pola LISA Bivariat",
    colorNA = "grey80"
  ) +
  tm_layout(
    main.title = "Peta LISA Bivariat IPM dan Kemiskinan Banten",
    main.title.position = "center",
    legend.outside = TRUE
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values'), 'colorNA'
## (rename to 'value.na') to fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

peta_lisa

## The map variable `fill` of the layer "polygons" contains a unique value.
## Therefore a categorical scale is applied (tm_scale_categorical).

Peta LISA bivariat menunjukkan bahwa seluruh kabupaten/kota di Provinsi Banten masuk kategori tidak signifikan. Meskipun nilai Moran lokal berbeda antarwilayah, belum ditemukan bukti statistik yang cukup untuk mengidentifikasi pola hubungan lokal antara IPM suatu wilayah dan persentase penduduk miskin di wilayah tetangganya.

Warna abu-abu bukan berarti semua wilayah memiliki nilai IPM atau kemiskinan yang sama. Warna tersebut menunjukkan bahwa pola lokal tidak dinyatakan signifikan oleh prosedur pengujian yang digunakan.

Kesimpulan

Analisis menunjukkan bahwa IPM dan persentase penduduk miskin memiliki hubungan linear negatif yang kuat pada wilayah yang sama. Namun, Global Bivariate Moran’s I belum menunjukkan hubungan spasial bivariat global yang signifikan pada ketiga matriks pembobot. Hasil LISA juga menunjukkan seluruh wilayah tidak signifikan berdasarkan output yang diperoleh.