1 Latar Belakang dan Tujuan

Pendidikan merupakan salah satu dimensi utama pembangunan manusia. Dua indikator yang lazim dipakai untuk menggambarkannya adalah Harapan Lama Sekolah (HLS), yaitu lamanya sekolah (tahun) yang diharapkan akan dirasakan anak usia 7 tahun pada masa mendatang, dan Rata-rata Lama Sekolah (RLS), yaitu rata-rata lama sekolah penduduk usia 25 tahun ke atas. HLS mencerminkan harapan untuk generasi muda, sedangkan RLS mencerminkan capaian generasi dewasa saat ini.

Pertanyaan yang diajukan pada laporan ini: apakah kabupaten/kota yang HLS-nya tinggi cenderung bertetangga dengan kabupaten/kota yang RLS-nya tinggi? Pertanyaan semacam ini tidak dapat dijawab dengan korelasi biasa karena korelasi mengabaikan letak wilayah. Karena itu digunakan Moran’s I bivariat (ukuran global) dan Local Indicator of Spatial Association (LISA) bivariat (ukuran lokal).

Tujuan analisis:

  1. Menyajikan statistik deskriptif HLS dan RLS dan menafsirkannya.
  2. Memvisualisasikan sebaran kedua indikator dalam bentuk peta tematik.
  3. Menyusun tiga pembobot spasial (invers jarak, k-nearest neighbors, dan Queen contiguity), menghitung Moran’s I bivariat global, lalu memilih pembobot terbaik.
  4. Mengidentifikasi klaster lokal dengan LISA bivariat.

2 Data dan Persiapan

2.1 Sumber data

Data berasal dari BPS Provinsi Banten tahun 2024 untuk 8 kabupaten/kota (file Harapan_Lama_Sekolah_...xlsx dan Mean_Years_Schooling_...xlsx). Kode wilayah mengikuti kode BPS.

data_banten <- tibble::tribble(
  ~kode,  ~wilayah,                 ~hls,  ~rls,
  "3601", "Kab. Pandeglang",        13.74,  7.16,
  "3602", "Kab. Lebak",             12.11,  6.61,
  "3603", "Kab. Tangerang",         12.88,  9.06,
  "3604", "Kab. Serang",            12.87,  7.88,
  "3671", "Kota Tangerang",         13.90, 11.14,
  "3672", "Kota Cilegon",           13.23, 10.39,
  "3673", "Kota Serang",            12.83,  8.93,
  "3674", "Kota Tangerang Selatan", 14.70, 11.86
)
n <- nrow(data_banten)
kable(data_banten, col.names = c("Kode", "Kabupaten/Kota", "HLS (tahun)", "RLS (tahun)"),
      caption = "Data HLS dan RLS menurut kabupaten/kota, Provinsi Banten 2024")
Data HLS dan RLS menurut kabupaten/kota, Provinsi Banten 2024
Kode Kabupaten/Kota HLS (tahun) RLS (tahun)
3601 Kab. Pandeglang 13.74 7.16
3602 Kab. Lebak 12.11 6.61
3603 Kab. Tangerang 12.88 9.06
3604 Kab. Serang 12.87 7.88
3671 Kota Tangerang 13.90 11.14
3672 Kota Cilegon 13.23 10.39
3673 Kota Serang 12.83 8.93
3674 Kota Tangerang Selatan 14.70 11.86

2.2 Data spasial

Peta dibuat tanpa shapefile: bentuk wilayah digambar manual (skematik) dari koordinat simpul pada chunk berikut, dengan simpul batas yang dipakai bersama sehingga ketetanggaan Queen tetap dapat dihitung. Bentuknya mendekati batas sebenarnya, tetapi bukan batas resmi.

# Peta skematik Banten (tanpa shapefile): poligon digambar manual dari koordinat piksel
# lalu dikonversi ke lon/lat. Simpul batas dipakai bersama sehingga Queen contiguity valid.
N <- list(
  a = c(322,182),
  b = c(335,168),
  d1 = c(350,148),
  d2 = c(362,118),
  d3 = c(378,108),
  d4 = c(400,120),
  k1 = c(410,166),
  k2 = c(392,200),
  c1 = c(355,176),
  c2 = c(372,186),
  k3 = c(392,216),
  k4 = c(402,227),
  k5 = c(420,236),
  k6 = c(440,236),
  k7 = c(458,206),
  k8 = c(448,185),
  k9 = c(440,160),
  s1 = c(446,150),
  s2 = c(480,148),
  s3 = c(500,163),
  t1 = c(493,190),
  t2 = c(510,280),
  m2 = c(522,290),
  j1 = c(470,272),
  j2 = c(445,264),
  l1 = c(420,265),
  p1 = c(300,245),
  p2 = c(360,251),
  p3 = c(405,245),
  q1 = c(550,158),
  q2 = c(600,155),
  q3 = c(628,172),
  q4 = c(632,190),
  u3 = c(612,205),
  u2 = c(590,218),
  u1 = c(570,230),
  v1 = c(603,240),
  e1 = c(645,240),
  e2 = c(650,270),
  e3 = c(642,288),
  w1 = c(610,288),
  w2 = c(603,265),
  m1 = c(560,288),
  z1 = c(530,320),
  z2 = c(515,345),
  z3 = c(525,385),
  z4 = c(555,435),
  z5 = c(540,450),
  z6 = c(492,455),
  z7 = c(482,500),
  z8 = c(470,520),
  z9 = c(430,500),
  z10 = c(400,485),
  z11 = c(370,478),
  z12 = c(335,466),
  pl1 = c(322,462),
  pl2 = c(330,410),
  pl3 = c(342,372),
  pl4 = c(376,318),
  pl5 = c(405,298)
)
cincin <- list(
  pandeglang = list(N$p1, N$p2, N$p3, N$l1, N$pl5, N$pl4, N$pl3, N$pl2, N$pl1, c(285,470), c(245,472), c(205,468), c(165,475), c(158,462), c(150,442), c(125,432), c(100,436), c(72,440), c(78,428), c(100,418), c(125,398), c(138,395), c(150,412), c(165,420), c(190,415), c(205,395), c(220,355), c(232,335), c(262,338), c(290,312), c(302,280)),
  lebak = list(N$l1, N$j2, N$j1, N$t2, N$m2, N$z1, N$z2, N$z3, N$z4, N$z5, N$z6, N$z7, N$z8, N$z9, N$z10, N$z11, N$z12, N$pl1, N$pl2, N$pl3, N$pl4, N$pl5),
  kab_serang = list(N$a, N$b, N$c1, N$c2, N$k2, N$k3, N$k4, N$k5, N$k6, N$k7, N$k8, N$k9, N$s1, N$s2, N$s3, N$t1, N$t2, N$j1, N$j2, N$l1, N$p3, N$p2, N$p1),
  kab_tangerang = list(N$s3, N$q1, N$q2, N$q3, N$q4, N$u3, N$u2, N$u1, N$v1, N$w2, N$w1, N$m1, N$m2, N$t2, N$t1),
  kota_tangerang = list(N$q4, N$u3, N$u2, N$u1, N$v1, N$e1),
  kota_cilegon = list(N$b, N$d1, N$d2, N$d3, N$d4, N$k1, N$k2, N$c2, N$c1),
  kota_serang = list(N$k1, N$k2, N$k3, N$k4, N$k5, N$k6, N$k7, N$k8, N$k9),
  tangsel = list(N$v1, N$e1, N$e2, N$e3, N$w1, N$w2)
)
pulau_panaitan <- list(c(32,380), c(60,358), c(90,348), c(97,362), c(76,384), c(62,402), c(52,402), c(55,390))   # P. Panaitan (bagian Kab. Pandeglang)

px2ll <- function(p) c(105.13 + (p[1] - 30) * 0.00266, -5.92 - (p[2] - 108) * 0.00262)
buat_poly <- function(titik) {
  m <- t(sapply(titik, px2ll)); m <- rbind(m, m[1, ])
  st_polygon(list(m))
}
kode_peta <- c(pandeglang = "3601", lebak = "3602", kab_tangerang = "3603", kab_serang = "3604",
               kota_tangerang = "3671", kota_cilegon = "3672", kota_serang = "3673", tangsel = "3674")

geom <- lapply(names(kode_peta), function(k) {
  mp <- list(buat_poly(cincin[[k]]))
  if (k == "pandeglang") mp <- c(mp, list(buat_poly(pulau_panaitan)))
  st_multipolygon(lapply(mp, unclass))
})
peta <- st_sf(kode = unname(kode_peta), geometry = st_sfc(geom, crs = 4326))
peta <- left_join(peta, data_banten, by = "kode")
peta <- peta[match(data_banten$kode, peta$kode), ]

ada_peta <- TRUE
titik <- st_centroid(st_geometry(peta))
koord <- st_coordinates(titik)
cat("Peta skematik dimuat:", nrow(peta), "poligon (tanpa shapefile).")
## Peta skematik dimuat: 8 poligon (tanpa shapefile).

2.2.1 Memeriksa data spasial

print(peta)                                              # data wilayah
## Simple feature collection with 8 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.1 ymin: -6.999 xmax: 106.8 ymax: -5.92
## Geodetic CRS:  WGS 84
##   kode                wilayah   hls   rls                       geometry
## 1 3601        Kab. Pandeglang 13.74  7.16 MULTIPOLYGON (((105.8 -6.27...
## 2 3602             Kab. Lebak 12.11  6.61 MULTIPOLYGON (((106.2 -6.33...
## 3 3603         Kab. Tangerang 12.88  9.06 MULTIPOLYGON (((106.4 -6.06...
## 4 3604            Kab. Serang 12.87  7.88 MULTIPOLYGON (((105.9 -6.11...
## 5 3671         Kota Tangerang 13.90 11.14 MULTIPOLYGON (((106.7 -6.13...
## 6 3672           Kota Cilegon 13.23 10.39 MULTIPOLYGON (((105.9 -6.07...
## 7 3673            Kota Serang 12.83  8.93 MULTIPOLYGON (((106.1 -6.07...
## 8 3674 Kota Tangerang Selatan 14.70 11.86 MULTIPOLYGON (((106.7 -6.26...
plot(st_geometry(peta), main = "Batas kabupaten/kota, Provinsi Banten")  # bentuk peta

3 Statistik Deskriptif

skew <- function(x) { m <- length(x); mu <- mean(x); s <- sd(x)
  m / ((m - 1) * (m - 2)) * sum(((x - mu) / s)^3) }
kurt <- function(x) { m <- length(x); mu <- mean(x); s <- sd(x)
  m * (m + 1) / ((m - 1) * (m - 2) * (m - 3)) * sum(((x - mu) / s)^4) -
    3 * (m - 1)^2 / ((m - 2) * (m - 3)) }
ringkas <- function(x) c(
  "n" = length(x), "Rata-rata" = mean(x), "Median" = median(x),
  "Simpangan baku" = sd(x), "Varians" = var(x),
  "Minimum" = min(x), "Maksimum" = max(x), "Range" = diff(range(x)),
  "Q1" = unname(quantile(x, .25)), "Q3" = unname(quantile(x, .75)),
  "IQR" = IQR(x), "Koef. variasi (%)" = sd(x) / mean(x) * 100,
  "Skewness" = skew(x), "Kurtosis (excess)" = kurt(x))
tab_desk <- data.frame(Statistik = names(ringkas(data_banten$hls)),
                       HLS = ringkas(data_banten$hls),
                       RLS = ringkas(data_banten$rls), row.names = NULL)
kable(tab_desk, digits = 2, caption = "Statistik deskriptif HLS dan RLS (n = 8 kabupaten/kota)")
Statistik deskriptif HLS dan RLS (n = 8 kabupaten/kota)
Statistik HLS RLS
n 8.00 8.00
Rata-rata 13.28 9.13
Median 13.05 9.00
Simpangan baku 0.80 1.89
Varians 0.64 3.56
Minimum 12.11 6.61
Maksimum 14.70 11.86
Range 2.59 5.25
Q1 12.86 7.70
Q3 13.78 10.58
IQR 0.92 2.88
Koef. variasi (%) 6.05 20.67
Skewness 0.49 0.14
Kurtosis (excess) 0.22 -1.32 
h <- data_banten$hls; r <- data_banten$rls
cat(sprintf(
"**Tafsiran.** Rata-rata HLS sebesar %.2f tahun, sedangkan rata-rata RLS hanya %.2f tahun; selisih sekitar %.1f tahun menunjukkan bahwa generasi muda diharapkan bersekolah lebih lama daripada capaian penduduk dewasa saat ini. Koefisien variasi RLS (%.1f%%) jauh lebih besar daripada HLS (%.1f%%), artinya **kesenjangan antarwilayah lebih terasa pada capaian aktual (RLS) daripada pada harapan bersekolah (HLS)**. %s memiliki nilai terendah pada kedua indikator (HLS %.2f; RLS %.2f), sementara %s tertinggi (HLS %.2f; RLS %.2f). Median dan rata-rata yang berdekatan serta skewness yang kecil menunjukkan sebaran yang cukup simetris. Korelasi Pearson HLS–RLS adalah %.2f, yaitu hubungan positif yang kuat namun belum memperhitungkan lokasi.",
  mean(h), mean(r), mean(h) - mean(r), sd(r)/mean(r)*100, sd(h)/mean(h)*100,
  data_banten$wilayah[which.min(r)], min(h), min(r),
  data_banten$wilayah[which.max(r)], max(h), max(r), cor(h, r)))

Tafsiran. Rata-rata HLS sebesar 13.28 tahun, sedangkan rata-rata RLS hanya 9.13 tahun; selisih sekitar 4.2 tahun menunjukkan bahwa generasi muda diharapkan bersekolah lebih lama daripada capaian penduduk dewasa saat ini. Koefisien variasi RLS (20.7%) jauh lebih besar daripada HLS (6.0%), artinya kesenjangan antarwilayah lebih terasa pada capaian aktual (RLS) daripada pada harapan bersekolah (HLS). Kab. Lebak memiliki nilai terendah pada kedua indikator (HLS 12.11; RLS 6.61), sementara Kota Tangerang Selatan tertinggi (HLS 14.70; RLS 11.86). Median dan rata-rata yang berdekatan serta skewness yang kecil menunjukkan sebaran yang cukup simetris. Korelasi Pearson HLS–RLS adalah 0.73, yaitu hubungan positif yang kuat namun belum memperhitungkan lokasi.

4 Peta Tematik

# Bagi nilai menjadi 4 kelas interval sama lebar, label seperti legenda peta tematik
kelas <- function(x, k = 4) {
  br  <- seq(floor(min(x)), ceiling(max(x)), length.out = k + 1)
  lab <- sprintf("%.1f - %.1f", br[-(k + 1)], br[-1])
  factor(lab[cut(x, br, labels = FALSE, include.lowest = TRUE)], levels = lab)
}

buat_peta <- function(dat, var, judul, legenda, palet) {
  dat$kelas <- kelas(dat[[var]])
  ggplot(dat) +
    geom_sf(aes(fill = kelas), colour = "grey25", linewidth = .35) +
    geom_sf_text(aes(label = sprintf("%.2f", .data[[var]])), size = 2.8,
                 fontface = "bold", check_overlap = TRUE) +
    scale_fill_brewer(palette = palet, drop = FALSE, name = legenda) +
    guides(fill = guide_legend(title.position = "top", nrow = 2)) +
    labs(title = judul) +
    theme_void(base_size = 11) +
    theme(plot.title = element_text(size = 13, margin = margin(b = 6)),
          panel.border = element_rect(colour = "black", fill = NA, linewidth = .6),
          legend.position = "bottom", legend.justification = "left",
          legend.box.background = element_rect(colour = "black", fill = NA),
          legend.margin = margin(5, 8, 5, 8), plot.margin = margin(8, 8, 8, 8))
}

if (ada_peta) {
  p_hls <- buat_peta(peta, "hls", "Peta Sebaran Harapan Lama Sekolah Provinsi Banten", "HLS 2024 (tahun)", "Blues")
  p_rls <- buat_peta(peta, "rls", "Peta Sebaran Rata-rata Lama Sekolah Provinsi Banten", "RLS 2024 (tahun)", "OrRd")
  p_hls + p_rls
} else {
  data_banten |> mutate(X = koord[, 1], Y = koord[, 2]) |>
    pivot_longer(c(hls, rls), names_to = "indikator", values_to = "nilai") |>
    ggplot(aes(X, Y, size = nilai, colour = nilai)) +
    geom_point() + facet_wrap(~indikator) + scale_colour_viridis_c(direction = -1) +
    labs(title = "Sebaran HLS dan RLS (titik pusat perkiraan)", x = "Bujur", y = "Lintang") +
    theme_minimal()
}

urut_r <- data_banten$wilayah[order(-data_banten$rls)]
urut_h <- data_banten$wilayah[order(-data_banten$hls)]
rank_h <- rank(-data_banten$hls); rank_r <- rank(-data_banten$rls)
selisih <- rank_r - rank_h            # positif: peringkat HLS jauh lebih baik daripada RLS
kota <- grepl("^Kota", data_banten$wilayah)
cat(sprintf(
"**Tafsiran.** Pada peta RLS, tiga wilayah tertinggi adalah %s, sedangkan tiga terendah adalah %s. Pada peta HLS, tiga tertinggi adalah %s dan tiga terendah %s. Rata-rata RLS wilayah berstatus kota (%.2f tahun) lebih tinggi %.2f tahun daripada kabupaten (%.2f tahun), sedangkan selisih serupa pada HLS hanya %.2f tahun; jadi kontras kota–kabupaten **lebih tegas pada capaian (RLS) daripada pada harapan (HLS)**. Wilayah dengan peringkat HLS jauh lebih baik daripada peringkat RLS-nya adalah **%s** (peringkat HLS %d, peringkat RLS %d), menandakan perbaikan akses pendidikan pada generasi muda yang belum tercermin pada penduduk dewasa. Secara spasial, nilai tinggi cenderung bertumpu pada wilayah yang saling berdekatan, yang diuji secara formal pada bagian berikutnya.",
  paste(head(urut_r, 3), collapse = ", "), paste(tail(urut_r, 3), collapse = ", "),
  paste(head(urut_h, 3), collapse = ", "), paste(tail(urut_h, 3), collapse = ", "),
  mean(data_banten$rls[kota]), mean(data_banten$rls[kota]) - mean(data_banten$rls[!kota]),
  mean(data_banten$rls[!kota]),
  mean(data_banten$hls[kota]) - mean(data_banten$hls[!kota]),
  data_banten$wilayah[which.max(selisih)], as.integer(rank_h[which.max(selisih)]),
  as.integer(rank_r[which.max(selisih)])))

Tafsiran. Pada peta RLS, tiga wilayah tertinggi adalah Kota Tangerang Selatan, Kota Tangerang, Kota Cilegon, sedangkan tiga terendah adalah Kab. Serang, Kab. Pandeglang, Kab. Lebak. Pada peta HLS, tiga tertinggi adalah Kota Tangerang Selatan, Kota Tangerang, Kab. Pandeglang dan tiga terendah Kab. Serang, Kota Serang, Kab. Lebak. Rata-rata RLS wilayah berstatus kota (10.58 tahun) lebih tinggi 2.90 tahun daripada kabupaten (7.68 tahun), sedangkan selisih serupa pada HLS hanya 0.76 tahun; jadi kontras kota–kabupaten lebih tegas pada capaian (RLS) daripada pada harapan (HLS). Wilayah dengan peringkat HLS jauh lebih baik daripada peringkat RLS-nya adalah Kab. Pandeglang (peringkat HLS 3, peringkat RLS 7), menandakan perbaikan akses pendidikan pada generasi muda yang belum tercermin pada penduduk dewasa. Secara spasial, nilai tinggi cenderung bertumpu pada wilayah yang saling berdekatan, yang diuji secara formal pada bagian berikutnya.

5 Pembobot Spasial

Misalkan \(w_{ij}\) adalah bobot hubungan wilayah \(i\) dan \(j\). Seluruh matriks dinormalisasi per baris sehingga \(\sum_j w_{ij}=1\).

5.1 Invers jarak

\[w_{ij}=\frac{1/d_{ij}}{\sum_{j}1/d_{ij}},\quad i\neq j\]

dengan \(d_{ij}\) jarak antartitik pusat wilayah (km). Semua wilayah terhubung, tetapi wilayah yang lebih dekat berbobot lebih besar.

d_km <- matrix(as.numeric(st_distance(titik)), n, n) / 1000
w_id <- 1 / d_km
diag(w_id) <- 0
w_id <- w_id / rowSums(w_id)
dimnames(w_id) <- list(data_banten$kode, data_banten$kode)
lw_id <- mat2listw(w_id, style = "W")
kable(round(w_id, 3), caption = "Matriks pembobot invers jarak (row-standardized)")
Matriks pembobot invers jarak (row-standardized)
3601 3602 3603 3604 3671 3672 3673 3674
3601 0.000 0.213 0.112 0.171 0.097 0.150 0.158 0.098
3602 0.168 0.000 0.141 0.169 0.120 0.120 0.154 0.129
3603 0.054 0.085 0.000 0.126 0.296 0.093 0.132 0.213
3604 0.050 0.062 0.077 0.000 0.054 0.143 0.564 0.050
3671 0.043 0.067 0.273 0.082 0.000 0.066 0.084 0.385
3672 0.082 0.083 0.106 0.268 0.081 0.000 0.306 0.074
3673 0.045 0.056 0.079 0.555 0.055 0.160 0.000 0.050
3674 0.048 0.079 0.216 0.083 0.424 0.066 0.085 0.000

5.2 k-Nearest Neighbors (k = 3)

Setiap wilayah dihubungkan dengan \(k=3\) wilayah terdekat berdasarkan jarak titik pusat. Nilai \(k=3\) dipilih karena mendekati rata-rata jumlah tetangga pada pembobot Queen. Sensitivitas terhadap \(k\) diperiksa pada bagian 6.3.

nb_knn <- knn2nb(knearneigh(koord, k = 3, longlat = TRUE))
lw_knn <- nb2listw(nb_knn, style = "W")
data.frame(Wilayah = data_banten$wilayah,
           `Tetangga terdekat (k = 3)` = sapply(nb_knn, function(i) paste(data_banten$wilayah[i], collapse = "; ")),
           check.names = FALSE) |>
  kable(caption = "Daftar tetangga menurut k-NN (k = 3)")
Daftar tetangga menurut k-NN (k = 3)
Wilayah Tetangga terdekat (k = 3)
Kab. Pandeglang Kab. Lebak; Kab. Serang; Kota Serang
Kab. Lebak Kab. Pandeglang; Kab. Serang; Kota Serang
Kab. Tangerang Kota Tangerang; Kota Serang; Kota Tangerang Selatan
Kab. Serang Kab. Tangerang; Kota Cilegon; Kota Serang
Kota Tangerang Kab. Tangerang; Kota Serang; Kota Tangerang Selatan
Kota Cilegon Kab. Tangerang; Kab. Serang; Kota Serang
Kota Serang Kab. Tangerang; Kab. Serang; Kota Cilegon
Kota Tangerang Selatan Kab. Tangerang; Kota Tangerang; Kota Serang 
w_knn <- nb2mat(nb_knn, style = "W")
dimnames(w_knn) <- list(data_banten$kode, data_banten$kode)
kable(round(w_knn, 3), caption = "Matriks pembobot k-NN (row-standardized)")
Matriks pembobot k-NN (row-standardized)
3601 3602 3603 3604 3671 3672 3673 3674
3601 0.000 0.333 0.000 0.333 0.000 0.000 0.333 0.000
3602 0.333 0.000 0.000 0.333 0.000 0.000 0.333 0.000
3603 0.000 0.000 0.000 0.000 0.333 0.000 0.333 0.333
3604 0.000 0.000 0.333 0.000 0.000 0.333 0.333 0.000
3671 0.000 0.000 0.333 0.000 0.000 0.000 0.333 0.333
3672 0.000 0.000 0.333 0.333 0.000 0.000 0.333 0.000
3673 0.000 0.000 0.333 0.333 0.000 0.333 0.000 0.000
3674 0.000 0.000 0.333 0.000 0.333 0.000 0.333 0.000

Perhatikan bahwa k-NN bersifat asimetris: wilayah A dapat menjadi tetangga B tanpa B menjadi tetangga A.

5.3 Queen contiguity

Dua wilayah bertetangga bila berbagi sisi atau titik batas.

if (ada_peta) {
  nb_queen <- poly2nb(peta, queen = TRUE)
} else {
  # Ketetanggaan dari pengetahuan geografis (cadangan)
  pasangan <- list(c("3601","3602"), c("3601","3604"), c("3602","3604"), c("3602","3603"),
                   c("3604","3673"), c("3604","3672"), c("3604","3603"), c("3673","3672"),
                   c("3603","3671"), c("3603","3674"), c("3671","3674"))
  A <- matrix(0, n, n, dimnames = list(data_banten$kode, data_banten$kode))
  for (p in pasangan) A[p[1], p[2]] <- A[p[2], p[1]] <- 1
  nb_queen <- mat2listw(A, style = "B")$neighbours
}
lw_q <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
data.frame(Wilayah = data_banten$wilayah,
           `Jumlah tetangga` = card(nb_queen),
           Tetangga = sapply(nb_queen, function(i) paste(data_banten$wilayah[i], collapse = "; ")),
           check.names = FALSE) |>
  kable(caption = "Daftar tetangga menurut Queen contiguity")
Daftar tetangga menurut Queen contiguity
Wilayah Jumlah tetangga Tetangga
Kab. Pandeglang 2 Kab. Lebak; Kab. Serang
Kab. Lebak 3 Kab. Pandeglang; Kab. Tangerang; Kab. Serang
Kab. Tangerang 4 Kab. Lebak; Kab. Serang; Kota Tangerang; Kota Tangerang Selatan
Kab. Serang 5 Kab. Pandeglang; Kab. Lebak; Kab. Tangerang; Kota Cilegon; Kota Serang
Kota Tangerang 2 Kab. Tangerang; Kota Tangerang Selatan
Kota Cilegon 2 Kab. Serang; Kota Serang
Kota Serang 2 Kab. Serang; Kota Cilegon
Kota Tangerang Selatan 2 Kab. Tangerang; Kota Tangerang 
w_q <- nb2mat(nb_queen, style = "W", zero.policy = TRUE)
dimnames(w_q) <- list(data_banten$kode, data_banten$kode)
kable(round(w_q, 3), caption = "Matriks pembobot Queen (row-standardized)")
Matriks pembobot Queen (row-standardized)
3601 3602 3603 3604 3671 3672 3673 3674
3601 0.000 0.50 0.000 0.500 0.00 0.0 0.0 0.00
3602 0.333 0.00 0.333 0.333 0.00 0.0 0.0 0.00
3603 0.000 0.25 0.000 0.250 0.25 0.0 0.0 0.25
3604 0.200 0.20 0.200 0.000 0.00 0.2 0.2 0.00
3671 0.000 0.00 0.500 0.000 0.00 0.0 0.0 0.50
3672 0.000 0.00 0.000 0.500 0.00 0.0 0.5 0.00
3673 0.000 0.00 0.000 0.500 0.00 0.5 0.0 0.00
3674 0.000 0.00 0.500 0.000 0.50 0.0 0.0 0.00

6 Moran’s I Bivariat Global

6.1 Formula

Misalkan \(z^{X}\) dan \(z^{Y}\) adalah nilai terstandarisasi HLS dan RLS (simpangan baku populasi). Dengan \(W\) ter-normalisasi baris:

\[I_{XY}=\frac{\sum_{i}z^{X}_{i}\left(\sum_{j}w_{ij}z^{Y}_{j}\right)}{n}\]

\(I_{XY}>0\) berarti nilai HLS yang tinggi cenderung berdampingan dengan nilai RLS yang tinggi (dan sebaliknya). Nilai harapan di bawah \(H_0\) (acak spasial) mendekati \(-1/(n-1)\). Signifikansi diuji dengan permutasi nilai \(Y\) sebanyak 9.999 kali, karena asumsi normalitas sulit dipenuhi untuk \(n=8\).

zs <- function(v) (v - mean(v)) / sqrt(mean((v - mean(v))^2))
zx <- zs(data_banten$hls)
zy <- zs(data_banten$rls)

moran_bv <- function(zx, zy, lw) sum(zx * lag.listw(lw, zy, zero.policy = TRUE)) / length(zx)

uji_bv <- function(lw, nsim = 9999) {
  I   <- moran_bv(zx, zy, lw)
  sim <- replicate(nsim, moran_bv(zx, sample(zy), lw))
  data.frame(I = I, `E[I]` = mean(sim), SD = sd(sim), z = (I - mean(sim)) / sd(sim),
             `p (1 sisi)` = (sum(sim >= I) + 1) / (nsim + 1),
             `p (2 sisi)` = (sum(abs(sim - mean(sim)) >= abs(I - mean(sim))) + 1) / (nsim + 1),
             `I (RLS,HLS)` = moran_bv(zy, zx, lw), check.names = FALSE)
}

6.2 Hasil dan pemilihan pembobot terbaik

daftar_w <- list("Invers jarak" = lw_id, "k-NN (k=3)" = lw_knn, "Queen" = lw_q)
hasil <- bind_rows(lapply(daftar_w, uji_bv), .id = "Pembobot")
kable(hasil, digits = 4, caption = "Moran's I bivariat global (X = HLS, Y = RLS)")
Moran’s I bivariat global (X = HLS, Y = RLS)
Pembobot I E[I] SD z p (1 sisi) p (2 sisi) I (RLS,HLS)
Invers jarak 0.0779 -0.0007 0.1110 0.7079 0.2506 0.4998 0.0860
k-NN (k=3) 0.1273 -0.0018 0.1503 0.8591 0.2082 0.4158 0.1609
Queen 0.2684 0.0005 0.2053 1.3048 0.1072 0.2137 0.2897

Pembobot terbaik dipilih berdasarkan nilai \(p\) terkecil (bukti autokorelasi paling kuat); bila seri, dipilih \(I\) terbesar.

terbaik <- hasil |> arrange(`p (1 sisi)`, desc(I)) |> slice(1)
lw_terbaik <- daftar_w[[terbaik$Pembobot]]
terbaik |> kable(digits = 4, caption = "Pembobot terbaik")
Pembobot terbaik
Pembobot I E[I] SD z p (1 sisi) p (2 sisi) I (RLS,HLS)
Queen 0.2684 5e-04 0.2053 1.305 0.1072 0.2137 0.2897 
sc <- data.frame(wilayah = data_banten$wilayah, zx = zx,
                 lag = lag.listw(lw_terbaik, zy, zero.policy = TRUE))
ggplot(sc, aes(zx, lag)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
  geom_smooth(method = "lm", se = FALSE, colour = "firebrick", linetype = "dashed", formula = y ~ x) +
  geom_point(size = 3) +
  geom_text_repel(aes(label = wilayah), size = 3.2) +
  labs(x = "z(HLS)", y = "Lag spasial z(RLS)",
       title = sprintf("Moran scatterplot bivariat — pembobot %s (I = %.3f)",
                       terbaik$Pembobot, terbaik$I)) +
  theme_minimal()

sig <- terbaik$`p (1 sisi)` < 0.05
cat(sprintf(
"**Tafsiran.** Pembobot **%s** memberikan Moran's I bivariat sebesar **%.3f** (nilai harapan %.3f; $z$ = %.2f; $p$ satu sisi = %.3f). Tanda %s menunjukkan kecenderungan %s. Pada taraf nyata 5%% hasil ini **%s**; %s Pembobot lain menghasilkan I yang jauh lebih kecil karena mengikutkan wilayah yang tidak berbatasan langsung sehingga pola kedekatan fisik terdilusi. Perlu diingat bahwa dengan hanya 8 wilayah, daya uji sangat terbatas sehingga hasil ini bersifat indikatif.",
  terbaik$Pembobot, terbaik$I, terbaik$`E[I]`, terbaik$z, terbaik$`p (1 sisi)`,
  ifelse(terbaik$I > 0, "positif", "negatif"),
  ifelse(terbaik$I > 0,
         "wilayah ber-HLS tinggi berdekatan dengan wilayah ber-RLS tinggi (dan sebaliknya)",
         "wilayah ber-HLS tinggi berdekatan dengan wilayah ber-RLS rendah"),
  ifelse(sig, "signifikan", "belum signifikan"),
  ifelse(sig, "terdapat bukti autokorelasi spasial bivariat.",
              "belum cukup bukti untuk menyatakan adanya autokorelasi spasial bivariat, meskipun arahnya konsisten dengan dugaan.")))

Tafsiran. Pembobot Queen memberikan Moran’s I bivariat sebesar 0.268 (nilai harapan 0.000; \(z\) = 1.30; \(p\) satu sisi = 0.107). Tanda positif menunjukkan kecenderungan wilayah ber-HLS tinggi berdekatan dengan wilayah ber-RLS tinggi (dan sebaliknya). Pada taraf nyata 5% hasil ini belum signifikan; belum cukup bukti untuk menyatakan adanya autokorelasi spasial bivariat, meskipun arahnya konsisten dengan dugaan. Pembobot lain menghasilkan I yang jauh lebih kecil karena mengikutkan wilayah yang tidak berbatasan langsung sehingga pola kedekatan fisik terdilusi. Perlu diingat bahwa dengan hanya 8 wilayah, daya uji sangat terbatas sehingga hasil ini bersifat indikatif.

6.3 Uji sensitivitas terhadap k

sens <- bind_rows(lapply(2:5, function(k) {
  lw_k <- nb2listw(knn2nb(knearneigh(koord, k = k, longlat = TRUE)), style = "W")
  cbind(k = k, uji_bv(lw_k, nsim = 4999))
}))
kable(sens[, c("k", "I", "z", "p (1 sisi)")], digits = 4,
      caption = "Sensitivitas Moran's I bivariat terhadap jumlah tetangga k-NN")
Sensitivitas Moran’s I bivariat terhadap jumlah tetangga k-NN
k I z p (1 sisi)
2 0.1966 0.7227 0.2574
3 0.1273 0.8621 0.2046
4 0.1174 1.4890 0.0754
5 -0.1041 -1.2035 0.8778

Perubahan nilai \(I\) dan tanda yang cukup besar antar-\(k\) menunjukkan bahwa kesimpulan sangat bergantung pada definisi ketetanggaan, sehingga pembobot berbasis batas wilayah (Queen) lebih mudah dipertanggungjawabkan secara substantif.

7 LISA Bivariat

7.1 Formula

\[I_{i}^{XY}=z^{X}_{i}\sum_{j}w_{ij}z^{Y}_{j}\]

Kuadran ditentukan oleh tanda \(z^{X}_{i}\) dan lag spasial \(\sum_j w_{ij}z^Y_j\):

Kategori Arti
High-High HLS tinggi, tetangga ber-RLS tinggi
Low-Low HLS rendah, tetangga ber-RLS rendah
High-Low HLS tinggi, tetangga ber-RLS rendah (pencilan spasial)
Low-High HLS rendah, tetangga ber-RLS tinggi (pencilan spasial)

Signifikansi memakai permutasi bersyarat (nilai \(Y\) tetangga diacak dari wilayah lain dengan \(z^X_i\) tetap) sebanyak 9.999 kali.

lisa_bv <- function(lw, nsim = 9999) {
  lag_y <- lag.listw(lw, zy, zero.policy = TRUE)
  Ii <- zx * lag_y
  p <- sapply(seq_len(n), function(i) {
    nb <- lw$neighbours[[i]]; w <- lw$weights[[i]]
    if (length(nb) == 1 && nb[1] == 0) return(NA_real_)
    lain <- setdiff(seq_len(n), i)
    sim <- replicate(nsim, zx[i] * sum(w * zy[sample(lain, length(nb))]))
    if (Ii[i] >= 0) (sum(sim >= Ii[i]) + 1) / (nsim + 1) else (sum(sim <= Ii[i]) + 1) / (nsim + 1)
  })
  kuadran <- case_when(zx > 0 & lag_y > 0 ~ "High-High",
                       zx < 0 & lag_y < 0 ~ "Low-Low",
                       zx > 0 & lag_y < 0 ~ "High-Low",
                       TRUE               ~ "Low-High")
  data.frame(kode = data_banten$kode, wilayah = data_banten$wilayah,
             zHLS = zx, lagRLS = lag_y, Ii = Ii, p = p, kuadran = kuadran,
             kategori = ifelse(!is.na(p) & p < 0.05, kuadran, "Tidak signifikan"),
             kategori10 = ifelse(!is.na(p) & p < 0.10, kuadran, "Tidak signifikan"))
}

7.2 Hasil

lisa <- lisa_bv(lw_terbaik)
kable(lisa[, c("wilayah", "zHLS", "lagRLS", "Ii", "p", "kuadran", "kategori")], digits = 3,
      col.names = c("Wilayah", "z(HLS)", "Lag z(RLS)", "Ii", "Pseudo p", "Kuadran", "Kategori (p < 0,05)"),
      caption = sprintf("LISA bivariat dengan pembobot %s", terbaik$Pembobot))
LISA bivariat dengan pembobot Queen
Wilayah z(HLS) Lag z(RLS) Ii Pseudo p Kuadran Kategori (p < 0,05)
Kab. Pandeglang 0.609 -1.067 -0.650 0.051 High-Low Tidak signifikan
Kab. Lebak -1.561 -0.621 0.969 0.057 Low-Low Tidak signifikan
Kab. Tangerang -0.536 0.138 -0.074 0.404 Low-High Tidak signifikan
Kab. Serang -0.549 -0.396 0.217 0.049 Low-Low Low-Low
Kota Tangerang 0.822 0.754 0.620 0.102 High-High Tidak signifikan
Kota Cilegon -0.070 -0.410 0.029 0.373 Low-Low Tidak signifikan
Kota Serang -0.602 0.004 -0.002 0.567 Low-High Tidak signifikan
Kota Tangerang Selatan 1.887 0.550 1.038 0.093 High-High Tidak signifikan 
warna <- c("High-High" = "#d7191c", "Low-Low" = "#2c7bb6", "High-Low" = "#fdae61",
           "Low-High" = "#abd9e9", "Tidak signifikan" = "grey88")
if (ada_peta) {
  peta_l <- left_join(peta[, "kode"], lisa, by = "kode")
  ggplot(peta_l) +
    geom_sf(aes(fill = kategori), colour = "grey30") +
    geom_sf_text(aes(label = wilayah), size = 2.8) +
    scale_fill_manual(values = warna, name = "Klaster (p < 0,05)", drop = FALSE) +
    labs(title = "LISA bivariat: HLS (X) terhadap RLS tetangga (Y)") +
    theme_minimal()
} else {
  ggplot(lisa |> mutate(X = koord[, 1], Y = koord[, 2]), aes(X, Y, fill = kategori)) +
    geom_point(shape = 21, size = 6) + geom_text_repel(aes(label = wilayah), size = 3) +
    scale_fill_manual(values = warna, name = "Klaster (p < 0,05)", drop = FALSE) +
    labs(title = "LISA bivariat (titik pusat perkiraan)", x = "Bujur", y = "Lintang") +
    theme_minimal()
}

7.3 Interpretasi

arti <- c("High-High" = "HLS tinggi dan dikelilingi wilayah ber-RLS tinggi",
          "Low-Low"   = "HLS rendah dan dikelilingi wilayah ber-RLS rendah",
          "High-Low"  = "HLS tinggi tetapi dikelilingi wilayah ber-RLS rendah (pencilan spasial)",
          "Low-High"  = "HLS rendah tetapi dikelilingi wilayah ber-RLS tinggi (pencilan spasial)")
s5 <- lisa[lisa$kategori != "Tidak signifikan", ]
s10 <- lisa[lisa$kategori == "Tidak signifikan" & lisa$kategori10 != "Tidak signifikan", ]
if (nrow(s5) > 0) {
  cat("**Signifikan pada taraf 5%:**\n\n")
  for (i in seq_len(nrow(s5)))
    cat(sprintf("- **%s** (%s, $p$ = %.3f): %s.\n", s5$wilayah[i], s5$kategori[i], s5$p[i], arti[[s5$kategori[i]]]))
} else cat("Tidak ada wilayah yang signifikan pada taraf 5%.\n")

Signifikan pada taraf 5%:

  • Kab. Serang (Low-Low, \(p\) = 0.049): HLS rendah dan dikelilingi wilayah ber-RLS rendah.
if (nrow(s10) > 0) {
  cat("\n**Hanya signifikan pada taraf 10% (bukti lemah):**\n\n")
  for (i in seq_len(nrow(s10)))
    cat(sprintf("- %s (%s, $p$ = %.3f)\n", s10$wilayah[i], s10$kategori10[i], s10$p[i]))
}

Hanya signifikan pada taraf 10% (bukti lemah):

  • Kab. Pandeglang (High-Low, \(p\) = 0.051)
  • Kab. Lebak (Low-Low, \(p\) = 0.057)
  • Kota Tangerang Selatan (High-High, \(p\) = 0.093)
cat(sprintf("\nWilayah lainnya (%d) tidak menunjukkan klaster maupun pencilan yang signifikan.\n",
            n - nrow(s5) - nrow(s10)))

Wilayah lainnya (4) tidak menunjukkan klaster maupun pencilan yang signifikan.

Catatan kehati-hatian. Dengan \(n=8\) dan hanya 2–5 tetangga per wilayah, jumlah kombinasi permutasi sangat kecil sehingga pseudo p-value bersifat diskret dan kurang stabil. Selain itu, pengujian dilakukan serentak pada delapan wilayah tanpa koreksi multiple testing (mis. Bonferroni, \(\alpha = 0{,}05/8 \approx 0{,}006\)). Hasil LISA sebaiknya diperlakukan sebagai temuan eksploratif.

8 Kesimpulan

  1. HLS rata-rata jauh di atas RLS (generasi muda diharapkan bersekolah lebih lama). Variasi antarwilayah jauh lebih besar pada RLS daripada HLS.
  2. Dari tiga pembobot, Queen menghasilkan autokorelasi spasial bivariat terkuat dengan Moran’s I = 0.268 (p satu sisi = 0.107). Arah hubungan positif, tetapi belum signifikan pada taraf 5%.
  3. LISA bivariat menunjukkan kontras antara kawasan Tangerang Raya (cenderung klaster tinggi) dan kawasan selatan serta barat daya (cenderung klaster rendah), dengan Kab. Pandeglang sebagai pencilan spasial ber-HLS tinggi di tengah wilayah ber-RLS rendah.
  4. Secara kebijakan, peningkatan akses pendidikan menengah dan pendidikan kesetaraan bagi penduduk dewasa perlu diprioritaskan di Lebak, Pandeglang, dan Kab. Serang.

9 Informasi Sesi dan Referensi

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.50      patchwork_1.3.2 ggrepel_0.9.7   spdep_1.4-2    
## [5] spData_2.3.5    sf_1.1-0        ggplot2_4.0.2   tidyr_1.3.1    
## [9] dplyr_1.1.4    
## 
## loaded via a namespace (and not attached):
##  [1] s2_1.1.9           sass_0.4.10        generics_0.1.4     class_7.3-23      
##  [5] KernSmooth_2.23-26 lattice_0.22-7     digest_0.6.39      magrittr_2.0.4    
##  [9] evaluate_1.0.5     grid_4.5.2         RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] Matrix_1.7-4       jsonlite_2.0.0     e1071_1.7-17       DBI_1.2.3         
## [17] mgcv_1.9-3         purrr_1.2.0        scales_1.4.0       jquerylib_0.1.4   
## [21] cli_3.6.5          rlang_1.2.0        units_1.0-0        splines_4.5.2     
## [25] withr_3.0.2        cachem_1.1.0       yaml_2.3.11        tools_4.5.2       
## [29] deldir_2.0-4       boot_1.3-32        vctrs_0.6.5        R6_2.6.1          
## [33] proxy_0.4-29       lifecycle_1.0.4    classInt_0.4-11    pkgconfig_2.0.3   
## [37] pillar_1.11.1      bslib_0.9.0        gtable_0.3.6       glue_1.8.0        
## [41] Rcpp_1.1.0         xfun_0.61          tibble_3.3.0       tidyselect_1.2.1  
## [45] rstudioapi_0.17.1  farver_2.1.2       nlme_3.1-168       htmltools_0.5.9   
## [49] labeling_0.4.3     rmarkdown_2.30     wk_0.9.5           compiler_4.5.2    
## [53] S7_0.2.1           sp_2.2-1
  • Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
  • Wartenberg, D. (1985). Multivariate spatial correlation: a method for exploratory geographical analysis. Geographical Analysis, 17(4), 263–283.
  • Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer.
  • BPS Provinsi Banten (2024). Harapan Lama Sekolah dan Rata-rata Lama Sekolah menurut Kabupaten/Kota.