1. Pendahuluan

1.1 Latar Belakang

Penduduk dan fasilitas kesehatan merupakan dua aspek penting dalam perencanaan pembangunan wilayah. Jumlah penduduk dapat menggambarkan besarnya kebutuhan pelayanan, sedangkan jumlah fasilitas kesehatan memberikan gambaran awal mengenai ketersediaan sarana pelayanan. Persebaran kedua variabel tersebut mungkin berbeda antarwilayah dan dapat menunjukkan pola spasial.

Kabupaten Tangerang terdiri atas 29 kecamatan dengan karakteristik penduduk dan ketersediaan fasilitas yang beragam. Oleh karena itu, analisis statistik spasial digunakan untuk menggambarkan persebaran kedua variabel dan menilai apakah nilai suatu variabel berkaitan dengan nilai variabel di wilayah yang berdekatan.

Analisis ini menggunakan jumlah penduduk (ribu) sebagai variabel X dan jumlah fasilitas kesehatan tercatat sebagai variabel Y pada tahun 2024. Jumlah fasilitas merupakan penjumlahan kategori fasilitas yang tersedia pada dataset; hasilnya tidak langsung menggambarkan kecukupan pelayanan, aksesibilitas, atau mutu pelayanan kesehatan.

1.2 Tujuan

Mendeskripsikan jumlah penduduk dan fasilitas kesehatan per kecamatan.
Memvisualisasikan persebaran kedua variabel pada peta.
Membandingkan pembobot spasial invers jarak, k-nearest neighbours (k-NN), dan Queen Contiguity.
Menghitung dan menginterpretasikan Moran Bivariat Global serta LISA Bivariat.

2. Metode Penelitian

2.1 Data dan Wilayah Kajian

Unit analisis adalah 29 kecamatan di Kabupaten Tangerang, Provinsi Banten, tahun 2024. Data atribut dibaca dari file Excel yang telah disiapkan, sedangkan geometri wilayah dibaca dari shapefile batas kecamatan.

Variabel X: jumlah penduduk (ribu).
Variabel Y: jumlah fasilitas kesehatan tercatat.
Sumber atribut: BPS Kabupaten Tangerang (2025), tabel penduduk dan tabel jumlah fasilitas kesehatan menurut kecamatan tahun 2024.
Sumber geometri: Satu Data Indonesia, shapefile batas administrasi kecamatan.

Sebelum analisis, nama kecamatan diseragamkan dan data Excel digabungkan dengan shapefile berdasarkan nama kecamatan. Penggabungan menghasilkan 29 wilayah tanpa nilai kosong pada kedua variabel analisis.

2.2 Statistik Deskriptif dan Korelasi

Statistik deskriptif yang digunakan meliputi jumlah pengamatan, nilai minimum, maksimum, rata-rata, median, dan standar deviasi. Korelasi Pearson digunakan untuk menilai hubungan linear antara jumlah penduduk dan jumlah fasilitas kesehatan. Korelasi tidak membuktikan hubungan sebab-akibat.

2.3 Peta Tematik

Peta tematik digunakan untuk menampilkan perbedaan nilai jumlah penduduk dan jumlah fasilitas kesehatan antarkecamatan. Warna yang lebih tinggi pada legenda menunjukkan nilai variabel yang lebih besar.

2.4 Pembobot Spasial

Tiga pembobot dibandingkan:

Invers jarak: bobot antarwilayah ditentukan berdasarkan kebalikan jarak antartitik pusat kecamatan.
k-NN (k = 3): setiap kecamatan dihubungkan dengan tiga tetangga terdekat.
Queen Contiguity: dua kecamatan dianggap bertetangga apabila poligonnya berbagi sisi atau titik.

Pembobot distandardisasi menggunakan row-standardization (style W). Jarak dihitung setelah geometri ditransformasikan ke sistem koordinat proyeksi UTM zona 48S.

2.5 Moran Bivariat Global

Moran Bivariat Global digunakan untuk menilai hubungan antara nilai X di suatu wilayah dan nilai Y pada wilayah-wilayah tetangganya. Signifikansi diuji dengan 9.999 permutasi Monte Carlo dan taraf signifikansi 5%.

Statistik pada skrip dihitung menggunakan:

\[ I_{XY} = \frac{z_X^\mathsf{T}Wz_Y}{z_X^\mathsf{T}z_X} \]

dengan \(z_X\) dan \(z_Y\) adalah variabel yang distandardisasi, dan \(W\) adalah matriks bobot spasial. Karena definisi dan normalisasi Moran bivariat dapat berbeda antarliteratur dan perangkat lunak, nilai statistik pada laporan ini dibaca berdasarkan rumus di atas. Hasilnya juga divisualisasikan dalam Moran scatterplot, yaitu diagram pencar antara \(z_X\) dan \(Wz_Y\) yang kemiringan garis regresinya sama dengan nilai \(I_{XY}\).

2.6 LISA Bivariat

LISA Bivariat digunakan untuk mengenali pola lokal antara jumlah penduduk suatu kecamatan dan jumlah fasilitas kesehatan pada kecamatan tetangganya. Pada skrip ini, statistik lokal dihitung sebagai:

\[ I_i^{XY}=z_{X_i}\sum_j w_{ij}z_{Y_j} \]

Kategori kuadran ditentukan dari tanda nilai standar X dan spatial lag Y:

High–High: X tinggi dan lag Y tinggi.
Low–Low: X rendah dan lag Y rendah.
High–Low: X tinggi dan lag Y rendah.
Low–High: X rendah dan lag Y tinggi.

Kecamatan disebut signifikan jika p-value permutasi kurang dari 0,05. P-value lokal pada hasil ini belum dikoreksi untuk pengujian berganda; hasil LISA sebaiknya diperlakukan sebagai eksploratif.

3. Hasil dan Pembahasan

3.1 Persiapan Data

data_analisis <- read_excel(
  "C:/Users/rahay/Downloads/Data_Tangerang_2024_Rapih.xlsx"
)

peta_kecamatan <- st_read(
  "C:/Users/rahay/Downloads/batas_administrasi_kecamatan_ln/batas_administrasi_kecamatan_ln.shp",
  quiet = TRUE
)

names(data_analisis) <- trimws(names(data_analisis))

# Menyiapkan variabel analisis
data_analisis <- data_analisis %>%
  mutate(
    Kecamatan = trimws(as.character(Kecamatan)),
    Penduduk = as.numeric(`Jumlah Penduduk (Ribu)`),
    Fasilitas = as.numeric(`Jumlah Fasilitas Tercatat`)
  )

# Menyamakan nama kecamatan pada SHP
peta_kecamatan <- peta_kecamatan %>%
  mutate(
    Kecamatan = case_when(
      wadmkc == "Mekarbaru" ~ "Mekar Baru",
      wadmkc == "Sindangjaya" ~ "Sindang Jaya",
      TRUE ~ trimws(wadmkc)
    )
  )

# Memeriksa duplikasi
if (anyDuplicated(data_analisis$Kecamatan)) {
  stop("Terdapat nama kecamatan duplikat pada Excel.")
}

if (anyDuplicated(peta_kecamatan$Kecamatan)) {
  stop("Terdapat nama kecamatan duplikat pada SHP.")
}

# Menggabungkan data dengan peta
peta_analisis <- peta_kecamatan %>%
  left_join(
    data_analisis %>%
      select(Kecamatan, Penduduk, Fasilitas),
    by = "Kecamatan"
  )

# Validasi hasil penggabungan
cat("Jumlah kecamatan pada Excel:",
    nrow(data_analisis), "\n")

## Jumlah kecamatan pada Excel: 29

cat("Jumlah kecamatan pada SHP:",
    nrow(peta_kecamatan), "\n")

## Jumlah kecamatan pada SHP: 29

cat("Jumlah wilayah hasil penggabungan:",
    nrow(peta_analisis), "\n")

## Jumlah wilayah hasil penggabungan: 29

cat("Wilayah yang belum cocok:\n")

## Wilayah yang belum cocok:

print(setdiff(
  data_analisis$Kecamatan,
  peta_kecamatan$Kecamatan
))

## character(0)

if (nrow(peta_analisis) != 29 ||
    anyNA(peta_analisis$Penduduk) ||
    anyNA(peta_analisis$Fasilitas)) {
  stop("Periksa kembali pencocokan wilayah atau data yang kosong.")
}

3.2 Statistik Deskriptif

stat_deskriptif <- data.frame(
  Variabel = c(
    "Jumlah Penduduk (Ribu)",
    "Jumlah Fasilitas Kesehatan"
  ),
  N = c(
    sum(!is.na(peta_analisis$Penduduk)),
    sum(!is.na(peta_analisis$Fasilitas))
  ),
  Minimum = c(
    min(peta_analisis$Penduduk, na.rm = TRUE),
    min(peta_analisis$Fasilitas, na.rm = TRUE)
  ),
  Maksimum = c(
    max(peta_analisis$Penduduk, na.rm = TRUE),
    max(peta_analisis$Fasilitas, na.rm = TRUE)
  ),
  Rata_rata = c(
    mean(peta_analisis$Penduduk, na.rm = TRUE),
    mean(peta_analisis$Fasilitas, na.rm = TRUE)
  ),
  Median = c(
    median(peta_analisis$Penduduk, na.rm = TRUE),
    median(peta_analisis$Fasilitas, na.rm = TRUE)
  ),
  Standar_deviasi = c(
    sd(peta_analisis$Penduduk, na.rm = TRUE),
    sd(peta_analisis$Fasilitas, na.rm = TRUE)
  )
)

print(stat_deskriptif)

##                     Variabel  N Minimum Maksimum  Rata_rata Median
## 1     Jumlah Penduduk (Ribu) 29    44.9    258.5 116.310345  108.2
## 2 Jumlah Fasilitas Kesehatan 29     1.0      9.0   2.448276    2.0
##   Standar_deviasi
## 1       52.431597
## 2        1.938081

Analisis mencakup 29 kecamatan. Berdasarkan output statistik deskriptif yang diperoleh, ringkasannya adalah:

Statistik	Jumlah penduduk (ribu)	Jumlah fasilitas kesehatan
Jumlah wilayah	29	29
Minimum	44,9	1
Maksimum	258,5	9
Rata-rata	116,31	2,45
Median	108,2	2
Standar deviasi	52,43	1,94

Jumlah penduduk rata-rata adalah 116,31 ribu jiwa per kecamatan, dengan rentang 44,9–258,5 ribu jiwa. Standar deviasi sebesar 52,43 ribu menunjukkan variasi jumlah penduduk yang cukup besar antarwilayah. Jumlah fasilitas kesehatan tercatat rata-rata 2,45 per kecamatan, dengan rentang 1–9 dan standar deviasi 1,94. Jumlah fasilitas perlu ditafsirkan sesuai definisi kolom dataset dan tidak boleh langsung dianggap sebagai ukuran kecukupan pelayanan.

3.3 Korelasi Pearson

print(cor.test(
  peta_analisis$Penduduk,
  peta_analisis$Fasilitas,
  method = "pearson"
))

## 
##  Pearson's product-moment correlation
## 
## data:  peta_analisis$Penduduk and peta_analisis$Fasilitas
## t = 4.8778, df = 27, p-value = 4.229e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4243675 0.8401749
## sample estimates:
##       cor 
## 0.6844228

Hasil korelasi Pearson pada output adalah:

Koefisien korelasi: \(r = 0{,}6844\).
p-value: \(4{,}229 \times 10^{-5}\).
Interval kepercayaan 95%: 0,4244 sampai 0,8402.

Koefisien positif 0,6844 menunjukkan hubungan linear positif yang cukup kuat pada data ini: kecamatan dengan jumlah penduduk lebih besar cenderung memiliki jumlah fasilitas tercatat lebih banyak. Nilai p yang lebih kecil dari 0,05 menunjukkan hubungan linear tersebut signifikan secara statistik. Hasil ini tidak membuktikan bahwa kenaikan jumlah penduduk menyebabkan penambahan fasilitas.

3.4 Visualisasi Peta Tematik

peta_label <- peta_analisis %>%
  mutate(
    titik_label = st_point_on_surface(geometry),
    lon = st_coordinates(titik_label)[, 1],
    lat = st_coordinates(titik_label)[, 2]
  )

ggplot(peta_analisis) +
  geom_sf(aes(fill = Penduduk), color = "grey30", linewidth = 0.25) +
  geom_text_repel(
    data = peta_label,
    aes(x = lon, y = lat, label = Kecamatan),
    size = 2.3, color = "black", bg.color = "white", bg.r = 0.15,
    segment.size = 0.3, segment.color = "grey40",
    max.overlaps = 40, min.segment.length = 0
  ) +
  scale_fill_viridis_c(option = "C", name = "Penduduk\n(ribu)") +
  labs(
    title = "Jumlah Penduduk per Kecamatan",
    subtitle = "Kabupaten Tangerang, 2024",
    caption = "Sumber atribut: BPS Kabupaten Tangerang (2025)",
    x = NULL, y = NULL
  ) +
  theme_minimal()

Berdasarkan peta, jumlah penduduk di Kabupaten Tangerang tahun 2024 terlihat berbeda-beda di setiap kecamatan. Pasarkemis memiliki jumlah penduduk tertinggi, yaitu 258,5 ribu jiwa. Sementara itu, kecamatan mekar baru memiliki jumlah penduduk paling sedikit, yaitu 44,9 ribu jiwa. Perbedaan warna pada peta menunjukkan adanya perbedaan jumlah penduduk antarwilayah.

ggplot(peta_analisis) +
  geom_sf(aes(fill = Fasilitas), color = "grey30", linewidth = 0.25) +
  geom_text_repel(
    data = peta_label,
    aes(x = lon, y = lat, label = Kecamatan),
    size = 2.3, color = "black", bg.color = "white", bg.r = 0.15,
    segment.size = 0.3, segment.color = "grey40",
    max.overlaps = 40, min.segment.length = 0
  ) +
    scale_fill_viridis_c(option = "C", name = "Jumlah\nfasilitas") +
  labs(
    title = "Jumlah Fasilitas Kesehatan per Kecamatan",
    subtitle = "Kabupaten Tangerang, 2024",
    caption = "Sumber atribut: BPS Kabupaten Tangerang (2025)",
    x = NULL, y = NULL
  ) +
  theme_minimal()

Fasilitas kesehatan paling banyak berada di Kelapa Dua (9 fasilitas), sedangkan sebagian besar kecamatan di barat laut hanya memiliki 1-2 fasilitas. Pola ini mirip dengan peta penduduk, sejalan dengan korelasi Pearson sebesar 0,68.

3.5 Pembobot Spasial

# Proyeksi UTM zona 48S untuk perhitungan jarak
peta_m <- st_transform(peta_analisis, 32748)

# Titik pusat setiap kecamatan
titik <- st_centroid(st_geometry(peta_m))
koordinat <- st_coordinates(titik)

nama_wilayah <- peta_m$Kecamatan

Pembobot invers jarak

jarak <- as.matrix(dist(koordinat))
diag(jarak) <- Inf

W_invers <- 1 / jarak
diag(W_invers) <- 0
W_invers[!is.finite(W_invers)] <- 0

lw_invers <- mat2listw(
  W_invers,
  style = "W",
  zero.policy = TRUE
)

lw_invers

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 29 
## Number of nonzero links: 812 
## Percentage nonzero weights: 96.55172 
## Average number of links: 28 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 29 841 29 2.897078 116.4334

Pembobot k-NN (k = 3)

nb_knn <- knn2nb(
  knearneigh(koordinat, k = 3),
  row.names = nama_wilayah
)

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

cat("\nJumlah tetangga k-NN:\n")

## 
## Jumlah tetangga k-NN:

print(card(nb_knn))

##  [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

Pembobot Queen Contiguity

nb_queen <- poly2nb(
  peta_m,
  queen = TRUE,
  row.names = nama_wilayah
)

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

cat("\nJumlah tetangga Queen:\n")

## 
## Jumlah tetangga Queen:

print(card(nb_queen))

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

cat("\nKecamatan tanpa tetangga Queen:\n")

## 
## Kecamatan tanpa tetangga Queen:

print(nama_wilayah[card(nb_queen) == 0])

## character(0)

Pembobot Queen menghasilkan jumlah tetangga yang bervariasi antarwilayah, dari 1 sampai 6 tetangga. Tidak ada kecamatan tanpa tetangga pada hasil Queen. Pada pembobot k-NN dengan k = 3, setiap kecamatan memiliki tiga tetangga terdekat berdasarkan konstruksi yang digunakan.

Ketiga pembobot mewakili konsep kedekatan yang berbeda. Invers jarak memberikan bobot berdasarkan jarak, k-NN menetapkan jumlah tetangga terdekat, sedangkan Queen berdasarkan persinggungan batas poligon. Perbedaan definisi tersebut dapat menghasilkan nilai statistik dan signifikansi yang berbeda.

Struktur ketetanggaan Queen dan k-NN dapat dilihat pada peta konektivitas berikut. Setiap garis menghubungkan titik pusat dua kecamatan yang dianggap bertetangga.

buat_garis <- function(nb, koordinat, crs) {
  pasangan <- do.call(rbind, lapply(seq_along(nb), function(i) {
    if (nb[[i]][1] == 0) return(NULL)
    cbind(i, nb[[i]])
  }))
  garis <- lapply(seq_len(nrow(pasangan)), function(r) {
    st_linestring(rbind(
      koordinat[pasangan[r, 1], ],
      koordinat[pasangan[r, 2], ]
    ))
  })
  st_sf(geometry = st_sfc(garis, crs = crs))
}

peta_konektivitas <- function(nb, judul) {
  ggplot() +
    geom_sf(data = peta_m, fill = "grey95",
            color = "grey60", linewidth = 0.25) +
    geom_sf(data = buat_garis(nb, koordinat, st_crs(peta_m)),
            color = "firebrick", linewidth = 0.4) +
    geom_sf(data = st_sf(geometry = titik), size = 1.2) +
    labs(
      title = judul,
      subtitle = "Garis menghubungkan centroid kecamatan bertetangga",
      x = NULL, y = NULL
    ) +
    theme_minimal()
}

print(peta_konektivitas(nb_queen, "Konektivitas Queen Contiguity"))

print(peta_konektivitas(nb_knn, paste0("Konektivitas k-NN (k = 3)")))

3.6 Moran Bivariat Global

biv_moran_global <- function(
    x, y, listw, nsim = 9999, seed = 123) {

  set.seed(seed)

  zx <- as.numeric(scale(x))
  zy <- as.numeric(scale(y))

  W <- listw2mat(listw)

  # Statistik Moran Bivariat Global
  I_observed <- sum(zx * as.numeric(W %*% zy)) /
    sum(zx^2)

  # Permutasi variabel Y
  I_simulasi <- replicate(nsim, {
    y_perm <- sample(zy, length(zy), replace = FALSE)

    sum(zx * as.numeric(W %*% y_perm)) /
      sum(zx^2)
  })

  # Uji dua arah
  p_value <- (
    sum(abs(I_simulasi) >= abs(I_observed)) + 1
  ) / (nsim + 1)

  data.frame(
    Moran_Bivariat = I_observed,
    P_value = p_value,
    Jumlah_permutasi = nsim
  )
}

x <- peta_m$Penduduk
y <- peta_m$Fasilitas

# Hitung untuk ketiga pembobot
hasil_invers <- biv_moran_global(x, y, lw_invers)
hasil_invers$Pembobot <- "Invers Jarak"

hasil_knn <- biv_moran_global(x, y, lw_knn)
hasil_knn$Pembobot <- paste0("k-NN (k = 3)")

hasil_queen <- biv_moran_global(x, y, lw_queen)
hasil_queen$Pembobot <- "Queen Contiguity"

hasil_global <- bind_rows(
  hasil_invers,
  hasil_knn,
  hasil_queen
) %>%
  select(
    Pembobot,
    Moran_Bivariat,
    P_value,
    Jumlah_permutasi
  )

kable(hasil_global, digits = 4,
      caption = "Hasil Moran Bivariat Global (9.999 permutasi)")

Hasil Moran Bivariat Global (9.999 permutasi)
Pembobot	Moran_Bivariat	P_value	Jumlah_permutasi
Invers Jarak	0.0594	0.0446	9999
k-NN (k = 3)	0.1035	0.3985	9999
Queen Contiguity	0.2601	0.0212	9999

# Pembobot terpilih: p-value terkecil
metode_terbaik <- hasil_global$Pembobot[
  which.min(hasil_global$P_value)
]

cat("\nPembobot dengan p-value terkecil:",
    metode_terbaik, "\n")

## 
## Pembobot dengan p-value terkecil: Queen Contiguity

Hasil permutasi 9.999 kali dari skrip:

Pembobot	Statistik Moran bivariat	p-value	Keputusan pada 5%
Invers jarak	0,0594	0,0446	Signifikan
k-NN	0.1035	0.3985	Tidak signifikan
Queen Contiguity	0,2601	0,0212	Signifikan

Berdasarkan hasil Moran Bivariat Global, pembobot invers jarak dan Queen Contiguity menunjukkan hubungan spasial positif yang signifikan karena memiliki p-value di bawah 0,05. Sementara itu, pembobot k-NN dengan tiga tetangga terdekat tidak menunjukkan hubungan yang signifikan secara statistik. Dalam analisis lanjutan, Queen Contiguity dipilih karena menggunakan hubungan ketetanggaan berdasarkan batas wilayah dan menghasilkan p-value terkecil, yaitu 0,0212. Pemilihan ini berlaku untuk analisis yang dilakukan pada data Kabupaten Tangerang tahun 2024 dan tidak berarti bahwa Queen selalu menjadi pembobot terbaik untuk semua kasus.

3.6.1 Moran Scatterplot

# Pembobot terpilih (p-value Moran global terkecil)
listw_terbaik <- switch(
  metode_terbaik,
  "Invers Jarak" = lw_invers,
  "Queen Contiguity" = lw_queen,
  lw_knn
)

zx <- as.numeric(scale(x))
zy <- as.numeric(scale(y))
lag_zy <- as.numeric(listw2mat(listw_terbaik) %*% zy)

I_terpilih <- sum(zx * lag_zy) / sum(zx^2)

df_scatter <- data.frame(
  Kecamatan = nama_wilayah,
  ZX = zx,
  Lag_ZY = lag_zy
)

ggplot(df_scatter, aes(x = ZX, y = Lag_ZY)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "firebrick",
              linewidth = 0.8, formula = y ~ x) +
  geom_text(aes(label = Kecamatan), size = 2.8, vjust = -0.8,
            check_overlap = TRUE) +
  labs(
    title = "Moran Scatterplot Bivariat",
    subtitle = paste0("Pembobot: ", metode_terbaik,
                      " | I = ", round(I_terpilih, 4)),
    x = "Jumlah penduduk (terstandardisasi)",
    y = "Spatial lag fasilitas kesehatan (terstandardisasi)"
  ) +
  theme_minimal()

Moran scatterplot menampilkan hubungan antara jumlah penduduk terstandardisasi di suatu kecamatan (sumbu X) dan rata-rata berbobot fasilitas kesehatan terstandardisasi di kecamatan tetangganya (sumbu Y). Kemiringan garis merah sama dengan nilai Moran Bivariat Global untuk pembobot terpilih (Queen Contiguity), yaitu 0.2601. Titik di kuadran kanan atas (High–High) dan kiri bawah (Low–Low) mendukung asosiasi spasial positif, sedangkan titik di kuadran kanan bawah (High–Low) dan kiri atas (Low–High) menunjukkan kecamatan yang berbeda dari pola tetangganya.

3.7 LISA Bivariat

lisa_bivariat <- function(
    x, y, listw, wilayah,
    nsim = 9999, seed = 456) {

  set.seed(seed)

  zx <- as.numeric(scale(x))
  zy <- as.numeric(scale(y))

  W <- listw2mat(listw)

  # Spatial lag variabel Y
  lag_y <- as.numeric(W %*% zy)

  # Statistik lokal bivariat
  I_lokal <- zx * lag_y

  # Permutasi variabel Y
  simulasi <- replicate(nsim, {
    y_perm <- sample(zy, length(zy), replace = FALSE)
    zx * as.numeric(W %*% y_perm)
  })

  # P-value dua arah setiap wilayah
  p_value <- vapply(seq_along(I_lokal), function(i) {
    (sum(abs(simulasi[i, ]) >= abs(I_lokal[i])) + 1) /
      (nsim + 1)
  }, numeric(1))

  # Kuadran berdasarkan X dan lag Y
  kuadran <- ifelse(
    zx >= 0 & lag_y >= 0, "High-High",
    ifelse(
      zx < 0 & lag_y < 0, "Low-Low",
      ifelse(
        zx >= 0 & lag_y < 0,
        "High-Low",
        "Low-High"
      )
    )
  )

  data.frame(
    Kecamatan = wilayah,
    X_Standar = zx,
    Lag_Y_Standar = lag_y,
    LISA_Bivariat = I_lokal,
    P_value = p_value,
    Kuadran = kuadran,
    Signifikan = p_value < 0.05
  )
}

# listw_terbaik didefinisikan pada chunk moran-scatterplot
hasil_lisa <- lisa_bivariat(
  x,
  y,
  listw_terbaik,
  nama_wilayah
) %>%
  mutate(
    Kategori_LISA = ifelse(
      Signifikan,
      Kuadran,
      "Tidak Signifikan"
    )
  )

cat("Jumlah wilayah per kategori LISA:\n")

## Jumlah wilayah per kategori LISA:

print(table(hasil_lisa$Kategori_LISA))

## 
##        High-High Tidak Signifikan 
##                2               27

cat("\nWilayah dengan LISA signifikan:\n")

## 
## Wilayah dengan LISA signifikan:

kable(
  hasil_lisa %>%
    filter(Signifikan) %>%
    select(Kecamatan, Kuadran, LISA_Bivariat, P_value),
  digits = 4
)

Kecamatan	Kuadran	LISA_Bivariat	P_value
Curug	High-High	1.5968	0.0088
Legok	High-High	0.1617	0.0447

Dengan pembobot Queen, output LISA menunjukkan dua kecamatan signifikan pada taraf 5%:

Kecamatan	Kuadran	Statistik lokal	p-value
Curug	High–High	1,5968	0,0088
Legok	High–High	0,1617	0,0447

Curug dan Legok sama-sama berkategori High-High, artinya penduduk di kedua kecamatan berada di atas rata-rata dan rata-rata fasilitas kesehatan di kecamatan tetangganya juga di atas rata-rata. Asosiasi di Curug jauh lebih kuat (1,5968; p = 0,0088) dibandingkan Legok (0,1617; p = 0,0447) yang signifikansinya mendekati batas 5%. Tidak ada kecamatan Low-Low, High-Low, maupun Low-High yang signifikan.

3.8 Peta Klaster LISA Bivariat

peta_lisa <- peta_m %>%
  left_join(
    hasil_lisa %>% select(Kecamatan, Kategori_LISA),
    by = "Kecamatan"
  )

warna_lisa <- c(
  "High-High"        = "#d7191c",
  "Low-Low"          = "#2c7bb6",
  "High-Low"         = "#fdae61",
  "Low-High"         = "#abd9e9",
  "Tidak Signifikan" = "grey85"
)
ggplot(peta_lisa) +
  geom_sf(
    aes(fill = Kategori_LISA),
    color = "grey30",
    linewidth = 0.25
  ) +
  geom_text_repel(
    data = filter(peta_lisa, Kategori_LISA != "Tidak Signifikan"),
    aes(label = Kecamatan, geometry = geometry),
    stat = "sf_coordinates",
    size = 3, fontface = "bold", color = "black",
    bg.color = "white", bg.r = 0.15,
    min.segment.length = 0
  ) +
  scale_fill_manual(values = warna_lisa, name = "Kategori LISA") +
  labs(
    title = "Peta Klaster LISA Bivariat",
    subtitle = paste("Pembobot:", metode_terbaik),
    caption = paste(
      "X = jumlah penduduk; Y = jumlah fasilitas kesehatan.",
      "Uji permutasi 9.999 kali; p < 0,05."
    ),
    x = NULL, y = NULL
  ) +
  theme_minimal()

Peta klaster LISA Bivariat dengan pembobot Queen menunjukkan dua kecamatan yang signifikan pada taraf 5%, yaitu Curug dan Legok, keduanya berkategori High-High. Kedua kecamatan saling berbatasan dan membentuk satu kelompok di bagian selatan-tengah Kabupaten Tangerang. Artinya, jumlah penduduk di Curug dan Legok berada di atas rata-rata kabupaten (116,31 ribu jiwa), dan kecamatan-kecamatan di sekitarnya juga memiliki jumlah fasilitas kesehatan yang relatif tinggi.

4. Kesimpulan

Analisis mencakup 29 kecamatan di Kabupaten Tangerang tahun 2024. Jumlah penduduk memiliki rata-rata 116,31 ribu jiwa, sedangkan jumlah fasilitas kesehatan tercatat rata-rata 2,45 per kecamatan.
Korelasi Pearson menunjukkan hubungan linear positif dan signifikan antara jumlah penduduk dan jumlah fasilitas kesehatan (\(r=0{,}6844\); \(p<0{,}001\)).
Peta tematik memperlihatkan variasi geografis jumlah penduduk dan fasilitas kesehatan antarkecamatan.
Dari tiga pembobot yang diuji, Queen Contiguity menghasilkan p-value terkecil (0,0212), diikuti invers jarak (0,0446). k-NN dengan \(k=3\) tidak signifikan (0,3985), berdasarkan rumus statistik pada skrip.
LISA Bivariat dengan Queen mengidentifikasi Curug dan Legok sebagai klaster High–High yang signifikan pada taraf 5%. Hasil lokal ini perlu diperlakukan sebagai eksploratif karena belum dilakukan koreksi pengujian berganda.
Jumlah fasilitas tidak cukup untuk menyimpulkan kecukupan pelayanan kesehatan. Analisis lanjutan dapat mempertimbangkan rasio fasilitas terhadap penduduk, luas wilayah, kapasitas, dan aksesibilitas.

Daftar Pustaka

Badan Pusat Statistik Kabupaten Tangerang. (2025). Jumlah Rumah Sakit Umum, Rumah Sakit Khusus, Puskesmas, Klinik Pratama, dan Posyandu Menurut Kecamatan di Kabupaten Tangerang, 2024. https://tangerangkab.bps.go.id/id/statistics-table/3/YmlzemNGUkNVblZLVVhOblREWnZXbkEzWld0eVVUMDkjMw%3D%3D/jumlah-rumah-sakit-umum--rumah-sakit-khusus--puskesmas--klinik-pratama--dan-posyandu-menurut-kecamatan-di-kabupaten-tangerang--2024.html?year=2024

Badan Pusat Statistik Kabupaten Tangerang. (2025). Penduduk, Laju Pertumbuhan Penduduk, Distribusi Persentase Penduduk, Kepadatan Penduduk, dan Rasio Jenis Kelamin Penduduk Menurut Kecamatan di Kabupaten Tangerang, 2024. https://tangerangkab.bps.go.id/id/statistics-table/3/V1ZSbFRUY3lTbFpEYTNsVWNGcDZjek53YkhsNFFUMDkjMw%3D%3D/population--population-growth-rate--percentage-distribution-of-population--population-density--and-population-sex-ratio-by-district-in-tangerang-regency--2024.html?year=2024

Satu Data Indonesia. Batas Administrasi Kecamatan [Dataset]. https://data.go.id/dataset/dataset/batas_administrasi_kecamatan_ar1

Silalahi, S. D. C., Azharuadin, M., Sihombing, R. D. M., Khasanah, W. E., & Sianipar, R. J. (2024). Analisis hotspot cluster: Tinjauan pola spasial dari ketersediaan fasilitas kesehatan dan penduduk di Kalimantan Indonesia. Jurnal Pendidikan Geografi Undiksha, 12(3), 341–351. https://doi.org/10.23887/jjpg.v12i3.83398

Analisis Autokorelasi Spasial Bivariat Jumlah Penduduk dan Fasilitas Kesehatan di Kabupaten Tangerang Tahun 2024

Rahayu Puspita Sari