1 Pendahuluan

1.1 Latar belakang

Padi adalah komoditas pangan strategis, dan besarnya produksi suatu daerah ditentukan terutama oleh luas lahan yang dipanen dan produktivitasnya. Menurut BPS Provinsi Banten, luas panen padi di Banten pada 2025 mencapai sekitar 345,42 ribu hektare, naik 46,33 ribu hektare (15,49%) dari 299,09 ribu hektare pada 2024. Produksinya diperkirakan sekitar 1,77 juta ton gabah kering giling (GKG), naik 14,41% dari sekitar 1,55 juta ton pada 2024. Pandeglang, Lebak, dan Serang disebut sebagai tiga penghasil padi terbesar, sedangkan Kota Tangerang dan Kota Cilegon mengalami penurunan produksi yang cukup besar (Banten Ekspres, 2026, yang mengutip rilis BPS Provinsi Banten). Data kabupaten/kota yang dipakai dalam laporan ini menunjukkan ketimpangan yang sama: tiga kabupaten terbesar menyumbang sekitar 79% luas panen dari tujuh wilayah yang dianalisis, sedangkan kota-kota kecil di utara hampir tidak berproduksi.

Ketimpangan antardaerah seperti ini memunculkan pertanyaan apakah daerah dengan luas panen besar cenderung bertetangga dengan daerah yang produksinya juga besar (dan sebaliknya), atau pola itu hanya mencerminkan ciri masing-masing daerah? Korelasi biasa tidak dapat menjawabnya karena tidak memperhitungkan letak dan ketetanggaan wilayah. Dasarnya adalah Hukum Pertama Geografi Tobler: “segala sesuatu saling berhubungan, tetapi yang berdekatan lebih berhubungan daripada yang berjauhan” (Tobler, 1970). Kekuatan hubungan seperti itu diukur dengan autokorelasi spasial, yaitu indeks Moran (Moran, 1950) dan turunan lokalnya, LISA (Anselin, 1995). Untuk dua variabel sekaligus digunakan Moran Bivariat, yang membandingkan nilai suatu variabel di sebuah wilayah dengan rata-rata tertimbang variabel lain di wilayah tetangganya (Wartenberg, 1985; Anselin, Syabri, & Smirnov, 2002).

Hasil ukuran autokorelasi bergantung pada cara ketetanggaan didefinisikan (pembobot spasial), pengaturan uji, dan implementasi perangkat lunak (Bivand & Wong, 2018). Karena itu laporan ini membandingkan tiga pembobot, yaitu invers jarak, k-Nearest Neighbors, dan Queen Contiguity, lalu memilih yang terbaik sebelum menghitung LISA Bivariat. Analisis dilakukan pada tingkat kabupaten/kota Provinsi Banten tahun 2025 dengan data resmi BPS.

1.2 Rumusan masalah

  1. Bagaimana gambaran (statistik deskriptif) luas panen dan produksi padi kabupaten/kota di Provinsi Banten tahun 2025, dan apakah keduanya saling berkorelasi?
  2. Bagaimana sebaran spasial luas panen dan produksi padi di Provinsi Banten?
  3. Apakah terdapat autokorelasi spasial bivariat global antara luas panen dan produksi padi, dan pembobot spasial mana (invers jarak, k-NN, Queen Contiguity) yang terbaik?
  4. Wilayah mana yang membentuk klaster (High-High, Low-Low) atau pencilan spasial (High-Low, Low-High) menurut LISA Bivariat?

1.3 Tujuan

  1. Mencari dan menjelaskan data spasial dua variabel numerik yang saling berkorelasi untuk kabupaten/kota di Provinsi Banten dari sumber resmi (BPS).
  2. Menyajikan statistik deskriptif beserta interpretasinya.
  3. Memvisualisasikan data dalam peta spasial dan menginterpretasikannya.
  4. Menghitung tiga pembobot spasial, Moran Bivariat Global dari ketiganya, memilih pembobot terbaik, dan menafsirkan hasilnya.
  5. Menghitung dan menafsirkan Local Moran Bivariat (LISA Bivariat).

1.4 Manfaat

  • Akademik: contoh penerapan Moran Bivariat Global dan LISA Bivariat lengkap dengan perbandingan pembobot, uji permutasi, dan koreksi uji berganda, yang seluruh kodenya dapat direproduksi (R Markdown).
  • Praktis: gambaran awal bagi pemangku kepentingan pertanian tentang apakah produksi padi Banten membentuk klaster spasial atau lebih ditentukan oleh ciri tiap daerah. Hasil ini bersifat eksploratif dan bukan dasar kebijakan tunggal.

1.5 Batasan penelitian

  • Data satu tahun (2025) pada tingkat kabupaten/kota; produksi dalam bentuk gabah kering giling (GKG).
  • Kota Tangerang Selatan tidak memiliki data pada tabel sumber sehingga analisis memakai 7 dari 8 kabupaten/kota.
  • Analisis menggambarkan keterkaitan spasial, bukan hubungan sebab-akibat.

1.6 Keterkaitan dengan soal tugas

No. soal Isi Bagian laporan
1 Data dua variabel berkorelasi dan sumbernya 2.1, 3.1, dan 3.2 (korelasi)
2 Statistik deskriptif dan interpretasi 3.2
3 Peta spasial dan interpretasi 3.3
4a Pembobot invers jarak, k-NN, Queen 3.4
4b Moran Bivariat Global, pembobot terbaik, interpretasi 3.5
4c LISA Bivariat dan interpretasi 3.6

2 Metode Penelitian

2.1 Sumber data dan variabel

Data yang digunakan adalah tabel Luas Panen, Produktivitas, dan Produksi Padi Menurut Kabupaten/Kota di Provinsi Banten, 2025 dari Badan Pusat Statistik (BPS) Provinsi Banten (berkas data_padi_banten_2025.xlsx). Data BPS berasal dari pendataan Kerangka Sampel Area (KSA) sepanjang 2025 (BPS Provinsi Banten, 2026). Total Banten pada tabel (345.421 ha dan 1.774.017 ton) sesuai dengan angka yang diberitakan dari rilis BPS Banten (Banten Ekspres, 2026). Batas wilayah kabupaten/kota berasal dari GADM versi 4.1 (GADM, 2022; berkas gadm41_IDN_2_pk.rds). Pengolahan dilakukan dengan R (R Core Team), paket sf (Pebesma, 2018) dan spdep (Bivand & Wong, 2018; Bivand, Pebesma, & Gómez-Rubio, 2013).

Variabel Simbol Satuan Peran
Luas panen padi X hektare (ha) variabel pertama
Produksi padi (gabah kering giling, GKG) Y ton variabel kedua
Produktivitas padi — kuintal/hektare (ku/ha) variabel pembanding (analisis sensitivitas)

Alasan pemilihan pasangan variabel. Penelitian ini memerlukan dua variabel yang saling berkorelasi. Luas panen dan produksi secara substantif saling berkaitan (produksi = luas panen × produktivitas), dan korelasinya sangat kuat (dibuktikan pada Bagian 3.2). Pasangan alternatif, luas panen dan produktivitas, hanya berkorelasi lemah dan tidak signifikan sehingga tidak dipakai sebagai analisis utama. Pasangan itu tetap disajikan sebagai analisis sensitivitas pada Bagian 3.7.

Unit analisis. Provinsi Banten memiliki 8 kabupaten/kota. Baris agregat “Banten” tidak digunakan. Kota Tangerang Selatan tidak memiliki nilai luas panen dan produksi pada tabel BPS (tertulis “–”), sehingga tidak dapat masuk analisis bivariat. Nilainya dibiarkan hilang (NA) dan tidak diganti nol. Akibatnya analisis memakai 7 wilayah, dan Tangerang Selatan hanya ditampilkan sebagai wilayah tanpa data pada peta.

2.2 Tahapan analisis

  1. Membaca dan membersihkan data, lalu menggabungkannya dengan peta.
  2. Statistik deskriptif, diagram pencar, dan korelasi Pearson.
  3. Peta tematik.
  4. Pembobot spasial W (matriks baku-baris): invers jarak, k-NN, dan Queen Contiguity.
  5. Moran Bivariat Global untuk ketiga pembobot, lalu memilih pembobot terbaik.
  6. LISA Bivariat dengan pembobot terpilih, pemetaan kuadran dan signifikansi.
  7. Analisis sensitivitas dengan pasangan luas panen dan produktivitas.

2.2.1 Diagram alir penelitian

library(ggplot2)
kotak <- function(id, x, y, w, h, teks, isi = "#E4EFDF", tebal = FALSE)
  data.frame(id, x, y, w, h, teks, isi, tebal)
bx <- rbind(
  kotak("a1", 2.6, 11.2, 3.6, 0.9, "Data BPS Banten 2025\n(luas panen, produksi)", "#D6E6F2"),
  kotak("a2", 7.4, 11.2, 3.6, 0.9, "Peta GADM 4.1\n(batas kabupaten/kota)", "#D6E6F2"),
  kotak("b",  5.0, 9.7, 7.6, 1.0, "Pembersihan dan pencocokan nama wilayah\n7 wilayah dianalisis (Tangsel tanpa data = NA)"),
  kotak("c",  5.0, 8.1, 7.6, 1.0, "Eksplorasi data: statistik deskriptif, scatterplot, korelasi Pearson,\npeta tematik luas panen dan produksi"),
  kotak("d1", 1.7, 6.4, 2.6, 0.9, "Queen\nContiguity", "#FCE9C8"),
  kotak("d2", 5.0, 6.4, 2.6, 0.9, "k-NN\n(k = 2)", "#FCE9C8"),
  kotak("d3", 8.3, 6.4, 2.6, 0.9, "Invers jarak\n(1/d)", "#FCE9C8"),
  kotak("e",  5.0, 4.9, 7.6, 0.9, "Moran Bivariat Global tiap pembobot\n(uji permutasi 9.999 kali, taraf 5%)"),
  kotak("f",  5.0, 3.5, 7.6, 0.9, "Pilih pembobot terbaik\n(p-value permutasi terkecil; seri: |I| terbesar)", "#F6D5D0"),
  kotak("g",  5.0, 2.1, 7.6, 0.9, "LISA Bivariat dengan pembobot terpilih (999 permutasi, koreksi BH)\nKuadran HH, LL, HL, LH dan peta signifikansi"),
  kotak("h",  5.0, 0.7, 7.6, 0.9, "Interpretasi, sensitivitas (luas panen vs produktivitas),\nkesimpulan dan keterbatasan", "#D6E6F2")
)
pan <- function(a, b, ax = NULL, bx_ = NULL) {
  A <- bx[bx$id == a, ]; B <- bx[bx$id == b, ]
  data.frame(x = ifelse(is.null(ax), A$x, ax), y = A$y - A$h / 2,
             xend = ifelse(is.null(bx_), B$x, bx_), yend = B$y + B$h / 2)
}
ar <- rbind(pan("a1","b",2.6,3.4), pan("a2","b",7.4,6.6),
            pan("b","c"),
            pan("c","d1",3.4,1.7), pan("c","d2"), pan("c","d3",6.6,8.3),
            pan("d1","e",1.7,3.2), pan("d2","e"), pan("d3","e",8.3,6.8),
            pan("e","f"), pan("f","g"), pan("g","h"))
ggplot() +
  geom_rect(data = bx, aes(xmin = x - w/2, xmax = x + w/2, ymin = y - h/2, ymax = y + h/2, fill = isi),
            color = "#1F4E3D", linewidth = 0.5) +
  geom_text(data = bx, aes(x, y, label = teks), size = 3.4, lineheight = 0.95) +
  geom_segment(data = ar, aes(x = x, y = y, xend = xend, yend = yend),
               arrow = arrow(length = unit(0.18, "cm"), type = "closed"), color = "#1F4E3D") +
  scale_fill_identity() + coord_cartesian(xlim = c(0, 10), ylim = c(0, 11.9)) +
  labs(title = "Diagram alir penelitian") +
  theme_void(base_size = 12) + theme(plot.title = element_text(face = "bold", hjust = 0.5))

Diagram menunjukkan bahwa tahap pembobot dan Moran Global dilakukan untuk ketiga pembobot sekaligus, lalu hanya satu pembobot (hasil seleksi) yang diteruskan ke LISA.

2.3 Rumus

Notasi mengikuti materi kuliah Statistika Spasial Pertemuan 4 (Moran’s I univariat; Moran, 1950), diperluas ke dua variabel (Wartenberg, 1985). Moran’s I univariat adalah

\[I=\frac{n}{S_0}\cdot\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2},\qquad S_0=\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}\]

Untuk Moran Bivariat Global, suku \((x_j-\bar{x})\) diganti dengan nilai variabel kedua di wilayah tetangga, dan penyebut memakai ragam kedua variabel:

\[I_{xy}=\frac{n}{S_0}\cdot\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(y_j-\bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2\;\sum_{i=1}^{n}(y_i-\bar{y})^2}}\]

Pembobot yang dipakai baku-baris, sehingga \(S_0=n\). Jika X dan Y distandardisasi menjadi \(z_x\) dan \(z_y\) (simpangan baku dengan pembagi \(n-1\)), rumus di atas menjadi bentuk yang dihitung spdep::moran_bv():

\[I_{xy}=\frac{\sum_{i}\sum_{j}w_{ij}\,z_{x,i}\,z_{y,j}}{n-1}\]

Anselin, Syabri, dan Smirnov (2002) menuliskan bentuk yang sama sebagai \(z_x'Wz_y/n\) dengan standardisasi memakai simpangan baku populasi (pembagi \(n\)); kedua penulisan menghasilkan nilai yang sama.

Nilai \(I_{xy}\) positif berarti wilayah dengan X tinggi cenderung bertetangga dengan wilayah dengan Y tinggi (dan X rendah dengan Y rendah), nilai negatif berarti sebaliknya. Nilai harapan di bawah \(H_0\) (tidak ada autokorelasi) adalah \(E(I)=-1/(n-1)\), dan statistik ujinya \(Z(I)=\dfrac{I-E(I)}{\sqrt{Var(I)}}\). Karena \(n\) hanya 7, ragam analitik tidak diandalkan; signifikansi diuji dengan uji permutasi (Monte Carlo), yaitu nilai Y diacak berulang kali (9.999 kali) untuk membentuk distribusi \(I\) di bawah \(H_0\). Rata-rata hasil permutasi pada Bagian 3.5 mendekati \(E(I)=-1/6\approx-0{,}167\) sebagaimana teori.

Local Moran Bivariat (LISA Bivariat; Anselin, 1995; Anselin et al., 2002). Pada materi, \(I_i=\dfrac{x_i-\bar{x}}{S^2}\sum_{j\neq i}w_{ij}(x_j-\bar{x})\). Dengan variabel yang sudah distandardisasi (ragam = 1) dan Y pada tetangga, bentuknya menjadi

\[I_i^{xy}=z_{x,i}\sum_{j}w_{ij}\,z_{y,j}=z_{x,i}\,(Wz_y)_i\]

dan rata-rata \(I_i^{xy}\) sebanding dengan \(I_{xy}\) (LISA adalah dekomposisi indeks global). Kuadran ditentukan oleh tanda \(z_{x,i}\) dan \((Wz_y)_i\) seperti pada materi: High-High (\(+,+\)), Low-Low (\(-,-\)), High-Low (\(+,-\)), dan Low-High (\(-,+\)). Signifikansi tiap \(I_i^{xy}\) diuji dengan permutasi 999 kali, lalu p-value dikoreksi Benjamini–Hochberg (Benjamini & Hochberg, 1995) karena menguji tujuh wilayah sekaligus. Taraf signifikansi 5%.

2.4 Kriteria keputusan

  • Moran Bivariat Global: \(H_0\) ditolak bila p-value permutasi < 0,05.
  • Pembobot terbaik: pembobot dengan p-value permutasi terkecil (bukti autokorelasi terkuat); bila seri, dipilih |\(I\)| terbesar. Pertimbangan substantif (struktur tetangga yang wajar) dibahas bersama hasilnya karena hasil bergantung pada pembobot (Bivand & Wong, 2018).
  • LISA Bivariat: wilayah dinyatakan signifikan bila p-value permutasi yang dikoreksi BH < 0,05. Kuadran tanpa signifikansi hanya dibaca sebagai arah pola, bukan klaster.

3 Hasil dan Pembahasan

3.1 Persiapan perangkat dan data

# install.packages(c("sf", "spdep", "dplyr", "ggplot2", "readxl", "knitr", "scales", "terra"))
library(terra)   # dimuat lebih dulu: berkas .rds GADM berformat terra (SpatVector)
library(sf)
library(spdep)
library(dplyr)
library(ggplot2)
library(readxl)
library(knitr)
library(scales)

Letakkan data_padi_banten_2025.xlsx dan gadm41_IDN_2_pk.rds di folder yang sama dengan berkas Rmd ini.

stopifnot(file.exists("data_padi_banten_2025.xlsx"))

nama_kabkota <- c("Pandeglang", "Lebak", "Tangerang", "Serang", "Kota Tangerang",
                  "Kota Cilegon", "Kota Serang", "Kota Tangerang Selatan")

mentah <- read_excel("data_padi_banten_2025.xlsx", col_types = "text")
names(mentah) <- c("wilayah", "luas_panen", "produktivitas", "produksi")

# as.numeric() mengubah tanda "–" dan "..." (data tidak tersedia) menjadi NA.
bersih <- mentah %>%
  filter(wilayah %in% c(nama_kabkota, "Banten")) %>%
  mutate(across(-wilayah, ~ suppressWarnings(as.numeric(.x))))

total_banten <- filter(bersih, wilayah == "Banten")

data_padi <- bersih %>%
  filter(wilayah != "Banten") %>%
  mutate(wilayah = case_when(
    wilayah == "Tangerang" ~ "Kabupaten Tangerang",
    wilayah == "Serang"    ~ "Kabupaten Serang",
    TRUE ~ wilayah
  ))

kable(data_padi, digits = 2, format.args = list(big.mark = ","),
      caption = "Data padi kabupaten/kota di Provinsi Banten, 2025 (sumber: BPS Provinsi Banten)")
Data padi kabupaten/kota di Provinsi Banten, 2025 (sumber: BPS Provinsi Banten)
wilayah luas_panen produktivitas produksi
Pandeglang 105,490 52.42 552,983
Lebak 88,819 47.62 422,943
Kabupaten Tangerang 58,298 52.80 307,816
Kabupaten Serang 79,503 52.80 419,792
Kota Tangerang 164 62.36 1,025
Kota Cilegon 1,664 50.57 8,413
Kota Serang 11,482 53.17 61,045
Kota Tangerang Selatan NA NA NA 
data_analisis <- filter(data_padi, !is.na(luas_panen), !is.na(produksi))

cat("Jumlah kabupaten/kota        :", nrow(data_padi), "\n")
## Jumlah kabupaten/kota        : 8
cat("Jumlah wilayah dianalisis    :", nrow(data_analisis), "\n")
## Jumlah wilayah dianalisis    : 7
cat("Wilayah tanpa data           :", paste(setdiff(data_padi$wilayah, data_analisis$wilayah), collapse = ", "), "\n")
## Wilayah tanpa data           : Kota Tangerang Selatan
cat("Jumlah luas panen 7 wilayah  :", sum(data_analisis$luas_panen), "ha (BPS Banten:", total_banten$luas_panen, "ha; selisih pembulatan)\n")
## Jumlah luas panen 7 wilayah  : 345420 ha (BPS Banten: 345421 ha; selisih pembulatan)
cat("Jumlah produksi 7 wilayah    :", sum(data_analisis$produksi), "ton (BPS Banten:", total_banten$produksi, "ton)\n")
## Jumlah produksi 7 wilayah    : 1774017 ton (BPS Banten: 1774017 ton)

Jumlah produksi tujuh wilayah sama persis dengan total Banten, dan jumlah luas panen hanya berbeda 1 ha (pembulatan). Ini memastikan data terbaca dengan benar dan Tangerang Selatan memang tidak menyumbang produksi pada tabel sumber.

3.2 Statistik deskriptif

ringkas <- function(v) {
  c(n = length(v), Rata2 = mean(v), Median = median(v), SD = sd(v),
    Min = min(v), Q1 = unname(quantile(v, .25)), Q3 = unname(quantile(v, .75)),
    Maks = max(v), KV_persen = 100 * sd(v) / mean(v))
}

tabel_desk <- rbind(
  "Luas panen (ha)"            = ringkas(data_analisis$luas_panen),
  "Produksi (ton)"             = ringkas(data_analisis$produksi),
  "Produktivitas (ku/ha)"      = ringkas(data_analisis$produktivitas)
)
kable(tabel_desk, digits = 2, format.args = list(big.mark = ","),
      caption = "Statistik deskriptif (n = 7 kabupaten/kota)")
Statistik deskriptif (n = 7 kabupaten/kota)
n Rata2 Median SD Min Q1 Q3 Maks KV_persen
Luas panen (ha) 7 49,345.71 58,298.0 44,395.90 164.00 6,573.0 84,161.00 105,490.00 89.97
Produksi (ton) 7 253,431.00 307,816.0 227,249.65 1,025.00 34,729.0 421,367.50 552,983.00 89.67
Produktivitas (ku/ha) 7 53.11 52.8 4.52 47.62 51.5 52.98 62.36 8.52 
# Nilai-nilai yang dipakai pada narasi interpretasi.
w_max_luas <- data_analisis$wilayah[which.max(data_analisis$luas_panen)]
w_min_luas <- data_analisis$wilayah[which.min(data_analisis$luas_panen)]
w_max_prod <- data_analisis$wilayah[which.max(data_analisis$produksi)]
w_min_prod <- data_analisis$wilayah[which.min(data_analisis$produksi)]
w_max_pvt  <- data_analisis$wilayah[which.max(data_analisis$produktivitas)]
w_min_pvt  <- data_analisis$wilayah[which.min(data_analisis$produktivitas)]
kv_luas <- tabel_desk["Luas panen (ha)", "KV_persen"]
kv_prod <- tabel_desk["Produksi (ton)", "KV_persen"]
kv_pvt  <- tabel_desk["Produktivitas (ku/ha)", "KV_persen"]
top3_luas <- data_analisis %>% arrange(desc(luas_panen)) %>% slice(1:3)
share_top3 <- 100 * sum(top3_luas$luas_panen) / sum(data_analisis$luas_panen)
rp <- function(v, d = 0) formatC(v, format = "f", digits = d, big.mark = ".", decimal.mark = ",")

Interpretasi.

  • Luas panen. Rata-rata luas panen 49.346 ha per kabupaten/kota, lebih rendah dari median 58.298 ha. Artinya sebaran sedikit menjulur ke kiri karena ada daerah dengan luas panen sangat kecil. Rentangnya sangat lebar: dari 164 ha di Kota Tangerang hingga 105.490 ha di Pandeglang. Koefisien variasi 90,0% menunjukkan keragaman yang sangat tinggi. Tiga daerah terbesar (Pandeglang, Lebak, Kabupaten Serang) menyumbang 79,3% luas panen dari ketujuh wilayah.
  • Produksi. Pola produksi mengikuti luas panen: tertinggi di Pandeglang (552.983 ton) dan terendah di Kota Tangerang (1.025 ton), dengan koefisien variasi 89,7%. Produksi terkonsentrasi di daerah agraris, sedangkan kota-kota kecil hampir tidak berkontribusi.
  • Produktivitas (pembanding). Produktivitas jauh lebih seragam (koefisien variasi hanya 8,5%), rata-rata 53,11 ku/ha. Tertinggi di Kota Tangerang (62,36 ku/ha) dan terendah di Lebak (47,62 ku/ha). Nilai tertinggi justru dimiliki daerah dengan luas panen paling kecil, sehingga produktivitas tinggi tidak otomatis berarti produksi besar.
  • Kesimpulan deskriptif. Perbedaan antardaerah di Banten terutama berasal dari skala lahan (luas panen), bukan dari produktivitas per hektare. Karena sebarannya sangat tidak merata, rata-rata saja tidak cukup untuk menggambarkan kondisi daerah.

3.2.1 Hubungan antarvariabel

ggplot(data_analisis, aes(luas_panen, produksi, label = wilayah)) +
  geom_smooth(method = "lm", se = FALSE, color = "grey60", linewidth = 0.6, formula = y ~ x) +
  geom_point(size = 3, color = "#1b6ca8") +
  geom_text(vjust = -0.9, size = 3.2, check_overlap = TRUE) +
  scale_x_continuous(labels = label_number(big.mark = ".", decimal.mark = ",")) +
  scale_y_continuous(labels = label_number(big.mark = ".", decimal.mark = ",")) +
  labs(title = "Luas panen vs produksi padi, kabupaten/kota Banten 2025",
       x = "Luas panen (ha)", y = "Produksi (ton)") +
  theme_minimal(base_size = 12)

uji_xy  <- cor.test(data_analisis$luas_panen, data_analisis$produksi, method = "pearson")
uji_xp  <- cor.test(data_analisis$luas_panen, data_analisis$produktivitas, method = "pearson")
uji_xy
## 
##  Pearson's product-moment correlation
## 
## data:  data_analisis$luas_panen and data_analisis$produksi
## t = 32.888, df = 5, p-value = 0.0000004884
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9837624 0.9996752
## sample estimates:
##       cor 
## 0.9976967

Korelasi Pearson luas panen dan produksi sebesar r = 0,998 (p < 0,001; selang kepercayaan 95%: 0,984 sampai 1,000). Kedua variabel berkorelasi positif sangat kuat dan signifikan, sehingga memenuhi syarat soal. Sebagai perbandingan, luas panen dan produktivitas hanya berkorelasi r = -0,501 (p = 0,252, tidak signifikan). Perlu dicatat bahwa korelasi ini bersifat nonspasial dan tidak menggantikan Moran Bivariat, yang menilai hubungan X di suatu wilayah dengan Y di wilayah tetangganya.

3.3 Peta spasial

peta_mentah <- readRDS("gadm41_IDN_2_pk.rds")
# Berkas dapat berupa PackedSpatVector (terra) atau sf; ubah ke sf.
if (inherits(peta_mentah, "PackedSpatVector")) peta_mentah <- terra::vect(peta_mentah)
peta_admin <- if (inherits(peta_mentah, "sf")) peta_mentah else sf::st_as_sf(peta_mentah)

peta_banten <- peta_admin %>%
  filter(tolower(trimws(NAME_1)) == "banten") %>%
  mutate(wilayah = case_when(
    NAME_2 == "Pandeglang"        ~ "Pandeglang",
    NAME_2 == "Lebak"             ~ "Lebak",
    NAME_2 == "Cilegon"           ~ "Kota Cilegon",
    NAME_2 == "Kota Serang"       ~ "Kota Serang",
    NAME_2 == "Kota Tangerang"    ~ "Kota Tangerang",
    NAME_2 == "Tangerang Selatan" ~ "Kota Tangerang Selatan",
    NAME_2 == "Serang"            ~ "Kabupaten Serang",
    NAME_2 == "Tangerang"         ~ "Kabupaten Tangerang",
    TRUE ~ NA_character_
  ))

stopifnot(!anyNA(peta_banten$wilayah), !anyDuplicated(peta_banten$wilayah),
          setequal(peta_banten$wilayah, data_padi$wilayah))

# Peta 8 wilayah (untuk tampilan) dan peta 7 wilayah (untuk analisis).
peta_8 <- peta_banten %>% left_join(data_padi, by = "wilayah") %>% st_make_valid()
peta_7 <- peta_8 %>% filter(!is.na(luas_panen), !is.na(produksi)) %>% arrange(wilayah)

# Titik label peta: di dalam poligon terbesar tiap wilayah; kota kecil digeser agar tidak menimpa kabupaten.
label_peta <- peta_8 %>%
  st_cast("POLYGON", warn = FALSE) %>%
  group_by(wilayah) %>% slice_max(as.numeric(st_area(geometry)), n = 1, with_ties = FALSE) %>% ungroup() %>%
  st_point_on_surface() %>%
  mutate(nama = sub("Kota Tangerang Selatan", "Tangsel", sub("Kabupaten ", "Kab. ", wilayah)))
xy_lab <- st_coordinates(label_peta)
label_peta$lx <- xy_lab[, 1]; label_peta$ly <- xy_lab[, 2]
geser <- c("Kota Cilegon" = -0.12, "Kota Serang" = 0.07, "Kota Tangerang" = 0.10, "Kota Tangerang Selatan" = 0.02)
label_peta$lx <- label_peta$lx + ifelse(label_peta$wilayah %in% names(geser), geser[label_peta$wilayah], 0)
label_peta$ly <- label_peta$ly + ifelse(label_peta$wilayah %in% c("Kota Serang", "Kota Tangerang"), 0.07, 0) -
                 ifelse(label_peta$wilayah == "Kota Tangerang Selatan", 0.04, 0)

# Proyeksi UTM 48S (meter) untuk jarak dan titik pusat wilayah.
peta_utm <- st_transform(peta_7, 32748)
titik    <- st_centroid(st_geometry(peta_utm))
koord    <- st_coordinates(titik)

# x, y diambil dari objek peta yang sama sehingga urutannya pasti selaras dengan pembobot.
nama_w <- peta_7$wilayah
x  <- peta_7$luas_panen
y  <- peta_7$produksi
yp <- peta_7$produktivitas

kable(st_drop_geometry(peta_7)[, c("wilayah", "luas_panen", "produksi", "produktivitas")],
      format.args = list(big.mark = ","), caption = "Data analisis yang telah dicocokkan dengan peta (7 wilayah)")
Data analisis yang telah dicocokkan dengan peta (7 wilayah)
wilayah luas_panen produksi produktivitas
Kabupaten Serang 79,503 419,792 52.80
Kabupaten Tangerang 58,298 307,816 52.80
Kota Cilegon 1,664 8,413 50.57
Kota Serang 11,482 61,045 53.17
Kota Tangerang 164 1,025 62.36
Lebak 88,819 422,943 47.62
Pandeglang 105,490 552,983 52.42

Seluruh 8 nama wilayah pada GADM berhasil dicocokkan dengan tabel BPS tanpa duplikasi.

peta_tematik <- function(peta, kolom, judul, nama_legenda, palet = "viridis", arah = -1) {
  ggplot(peta) +
    geom_sf(aes(fill = .data[[kolom]]), color = "white", linewidth = 0.4) +
    geom_label(data = label_peta, aes(lx, ly, label = nama), size = 2.6, fontface = "bold",
               fill = alpha("white", 0.75), label.size = 0, label.padding = unit(0.12, "lines")) +
    scale_fill_viridis_c(option = palet, direction = arah, na.value = "grey85",
                         labels = label_number(big.mark = ".", decimal.mark = ","),
                         name = nama_legenda) +
    labs(title = judul,
         subtitle = "Abu-abu = tidak ada data (Tangsel = Kota Tangerang Selatan)",
         caption = "Sumber: BPS Provinsi Banten; batas wilayah: GADM 4.1") +
    labs(x = NULL, y = NULL) +
    theme_minimal(base_size = 11) + theme(panel.grid = element_blank(),
                                          axis.text = element_blank())
}
peta_tematik(peta_8, "luas_panen", "Luas panen padi, Provinsi Banten 2025", "Luas panen\n(ha)", "viridis")

peta_tematik(peta_8, "produksi", "Produksi padi, Provinsi Banten 2025", "Produksi\n(ton)", "magma")

peta_tematik(peta_8, "produktivitas", "Produktivitas padi, Provinsi Banten 2025", "Produktivitas\n(ku/ha)", "cividis")

Interpretasi peta.

  • Luas panen dan produksi memiliki pola spasial yang hampir identik. Nilai tertinggi terkonsentrasi di bagian selatan dan barat: Pandeglang dan Lebak, disusul Kabupaten Serang dan Kabupaten Tangerang. Wilayah pesisir utara yang berkarakter perkotaan (Kota Tangerang, Kota Cilegon, Kota Serang) berwarna paling rendah. Pola ini kontras antara kabupaten agraris yang luas dan kota-kota yang kecil dan padat bangunan.
  • Produktivitas memiliki pola berbeda. Nilai tertinggi ada di Kota Tangerang, tetapi dengan luas panen hanya 164 ha. Produktivitas terendah di Lebak, padahal luas panen dan produksinya termasuk besar. Jadi sebaran produktivitas tidak mengikuti sebaran luas panen.
  • Peta hanya menunjukkan sebaran, bukan bukti autokorelasi. Daerah bernilai tinggi tampak mengelompok di selatan dan barat, tetapi kesan visual ini harus diuji secara statistik. Pengujiannya dilakukan melalui Moran Bivariat Global dan LISA pada bagian berikut.

3.4 Pembobot spasial

Tiga pembobot dibentuk dari 7 poligon yang dianalisis, semuanya baku-baris (jumlah bobot tiap wilayah = 1):

  1. Queen Contiguity: dua wilayah bertetangga bila berbagi batas atau titik sudut.
  2. k-NN dengan k = 2: tiap wilayah bertetangga dengan 2 titik pusat terdekat. Nilai k = 2 dipilih karena rata-rata jumlah tetangga Queen sekitar 2,3. Dengan demikian kepadatan kedua pembobot setara.
  3. Invers jarak: \(w_{ij}=1/d_{ij}\) dengan \(d_{ij}\) jarak antar titik pusat (km), sehingga semua wilayah saling berpengaruh dan yang lebih dekat berpengaruh lebih besar.
# 1. Queen Contiguity
nb_queen <- poly2nb(peta_7, queen = TRUE, row.names = nama_w)
lw_queen <- nb2listw(nb_queen, style = "W")

# 2. k-NN (k = 2)
nb_knn <- knn2nb(knearneigh(koord, k = 2), row.names = nama_w)
lw_knn <- nb2listw(nb_knn, style = "W")

# 3. Invers jarak (jarak dalam km)
jarak_km <- as.matrix(dist(koord)) / 1000
dimnames(jarak_km) <- list(nama_w, nama_w)
w_inv <- 1 / jarak_km
diag(w_inv) <- 0
lw_inv <- mat2listw(w_inv, row.names = nama_w, style = "W")

daftar_lw <- list("Invers jarak" = lw_inv, "k-NN (k = 2)" = lw_knn, "Queen Contiguity" = lw_queen)

stopifnot(all(card(nb_queen) > 0))   # tidak ada wilayah terisolasi

nama_pendek <- sub("Kabupaten ", "Kab. ", sub("Kota ", "Kota ", nama_w))
daftar_nb <- function(nb) sapply(nb, function(j) paste(nama_pendek[j], collapse = "; "))

kable(tibble(
  Wilayah = nama_w,
  `Jumlah tetangga Queen` = card(nb_queen),
  `Tetangga Queen` = daftar_nb(nb_queen),
  `Tetangga k-NN` = daftar_nb(nb_knn)
), caption = "Struktur ketetanggaan Queen dan k-NN")
Struktur ketetanggaan Queen dan k-NN
Wilayah Jumlah tetangga Queen Tetangga Queen Tetangga k-NN
Kabupaten Serang 5 Kab. Tangerang; Kota Cilegon; Kota Serang; Lebak; Pandeglang Kota Cilegon; Kota Serang
Kabupaten Tangerang 3 Kab. Serang; Kota Tangerang; Lebak Kota Serang; Kota Tangerang
Kota Cilegon 1 Kab. Serang Kab. Serang; Kota Serang
Kota Serang 1 Kab. Serang Kab. Serang; Kota Cilegon
Kota Tangerang 1 Kab. Tangerang Kab. Tangerang; Kota Serang
Lebak 3 Kab. Serang; Kab. Tangerang; Pandeglang Kab. Serang; Pandeglang
Pandeglang 2 Kab. Serang; Lebak Kab. Serang; Lebak 
tepi <- function(lw, nama_bobot) {
  bind_rows(lapply(seq_along(lw$neighbours), function(a) {
    tibble(i = a, j = lw$neighbours[[a]], w = lw$weights[[a]])
  })) %>%
    mutate(x = koord[i, 1], y = koord[i, 2], xend = koord[j, 1], yend = koord[j, 2],
           pembobot = nama_bobot)
}
tepi_semua <- bind_rows(lapply(names(daftar_lw), function(k) tepi(daftar_lw[[k]], k))) %>%
  mutate(pembobot = factor(pembobot, levels = names(daftar_lw)))

ggplot() +
  geom_sf(data = peta_utm, fill = "grey95", color = "grey60", linewidth = 0.3) +
  geom_segment(data = tepi_semua, aes(x = x, y = y, xend = xend, yend = yend, linewidth = w),
               color = "#c0392b", alpha = 0.6, lineend = "round") +
  geom_sf(data = titik, color = "black", size = 1.5) +
  scale_linewidth_continuous(range = c(0.15, 1.6), guide = "none") +
  facet_wrap(~ pembobot) +
  labs(title = "Jaringan ketetanggaan menurut pembobot",
       subtitle = "Ketebalan garis = besar bobot baku-baris; titik = pusat wilayah") +
  theme_void(base_size = 11) +
  theme(strip.text = element_text(face = "bold", margin = margin(t = 12, b = 4)),
        plot.title = element_text(size = 13), plot.subtitle = element_text(margin = margin(b = 6)),
        plot.margin = margin(10, 10, 10, 10))

Berikut matriks pembobot baku-baris yang dihasilkan (baris = wilayah ke-\(i\), kolom = tetangga ke-\(j\)).

tampil_w <- function(lw) {
  m <- round(listw2mat(lw), 3); dimnames(m) <- list(nama_pendek, nama_pendek); m
}
kable(tampil_w(lw_inv),   caption = "Matriks W — Invers jarak (baku-baris)")
Matriks W — Invers jarak (baku-baris)
Kab. Serang Kab. Tangerang Kota Cilegon Kota Serang Kota Tangerang Lebak Pandeglang
Kab. Serang 0.000 0.052 0.100 0.738 0.039 0.039 0.032
Kab. Tangerang 0.151 0.000 0.106 0.158 0.421 0.101 0.064
Kota Cilegon 0.307 0.113 0.000 0.309 0.091 0.089 0.091
Kota Serang 0.737 0.055 0.101 0.000 0.040 0.037 0.031
Kota Tangerang 0.127 0.480 0.097 0.131 0.000 0.101 0.064
Lebak 0.183 0.165 0.137 0.177 0.145 0.000 0.195
Pandeglang 0.183 0.126 0.170 0.176 0.111 0.234 0.000 
kable(tampil_w(lw_knn),   caption = "Matriks W — k-NN, k = 2 (baku-baris)")
Matriks W — k-NN, k = 2 (baku-baris)
Kab. Serang Kab. Tangerang Kota Cilegon Kota Serang Kota Tangerang Lebak Pandeglang
Kab. Serang 0.0 0.0 0.5 0.5 0.0 0.0 0.0
Kab. Tangerang 0.0 0.0 0.0 0.5 0.5 0.0 0.0
Kota Cilegon 0.5 0.0 0.0 0.5 0.0 0.0 0.0
Kota Serang 0.5 0.0 0.5 0.0 0.0 0.0 0.0
Kota Tangerang 0.0 0.5 0.0 0.5 0.0 0.0 0.0
Lebak 0.5 0.0 0.0 0.0 0.0 0.0 0.5
Pandeglang 0.5 0.0 0.0 0.0 0.0 0.5 0.0 
kable(tampil_w(lw_queen), caption = "Matriks W — Queen Contiguity (baku-baris)")
Matriks W — Queen Contiguity (baku-baris)
Kab. Serang Kab. Tangerang Kota Cilegon Kota Serang Kota Tangerang Lebak Pandeglang
Kab. Serang 0.000 0.200 0.2 0.2 0.000 0.200 0.200
Kab. Tangerang 0.333 0.000 0.0 0.0 0.333 0.333 0.000
Kota Cilegon 1.000 0.000 0.0 0.0 0.000 0.000 0.000
Kota Serang 1.000 0.000 0.0 0.0 0.000 0.000 0.000
Kota Tangerang 0.000 1.000 0.0 0.0 0.000 0.000 0.000
Lebak 0.333 0.333 0.0 0.0 0.000 0.000 0.333
Pandeglang 0.500 0.000 0.0 0.0 0.000 0.500 0.000

Pada Queen, Kota Tangerang, Kota Cilegon, dan Kota Serang hanya memiliki satu tetangga, sedangkan Kabupaten Serang memiliki lima. Akibatnya nilai lag di kota-kota kecil itu hanya ditentukan satu wilayah. Pada k-NN setiap wilayah memiliki dua tetangga, tetapi ketetanggaannya bergantung pada jarak titik pusat dan tidak selalu berbagi batas. Invers jarak melibatkan seluruh wilayah dengan bobot yang menurun menurut jarak.

3.5 Moran Bivariat Global

Hipotesis untuk tiap pembobot:

  • \(H_0\): tidak ada autokorelasi spasial bivariat (luas panen suatu wilayah tidak berkaitan dengan produksi di wilayah tetangganya).
  • \(H_1\): ada autokorelasi spasial bivariat.
hitung_global <- function(x, y, daftar_lw, nsim = 9999) {
  hasil <- lapply(names(daftar_lw), function(k) {
    set.seed(123)
    r  <- spdep::moran_bv(x, y, daftar_lw[[k]], nsim = nsim, scale = TRUE)
    I  <- as.numeric(r$t0)
    Ip <- as.numeric(r$t); Ip <- Ip[is.finite(Ip)]
    EI <- mean(Ip)
    p  <- (sum(abs(Ip - EI) >= abs(I - EI)) + 1) / (length(Ip) + 1)   # dua arah
    list(tab = tibble(pembobot = k, I = I, E_I = EI, p_permutasi = p),
         perm = tibble(pembobot = k, Ip = Ip, I_obs = I))
  })
  tab <- bind_rows(lapply(hasil, `[[`, "tab"))
  # Pembobot terbaik: bukti terkuat (p-value permutasi terkecil); bila seri, |I| terbesar.
  tab$terbaik <- seq_len(nrow(tab)) == order(tab$p_permutasi, -abs(tab$I))[1]
  tab$keputusan <- ifelse(tab$p_permutasi < 0.05, "Tolak H0", "Gagal tolak H0")
  list(tabel = tab, perm = bind_rows(lapply(hasil, `[[`, "perm")))
}

glob <- hitung_global(x, y, daftar_lw)
tabel_global <- glob$tabel
kable(tabel_global %>% mutate(terbaik = ifelse(terbaik, "✔", "")),
      digits = 4, col.names = c("Pembobot", "Moran Bivariat (I)", "E[I] permutasi",
                                "p-value permutasi", "Terbaik", "Keputusan (α = 5%)"),
      caption = "Moran Bivariat Global: luas panen (X) dan produksi (Y), 9.999 permutasi")
Moran Bivariat Global: luas panen (X) dan produksi (Y), 9.999 permutasi
Pembobot Moran Bivariat (I) E[I] permutasi p-value permutasi Terbaik Keputusan (α = 5%)
Invers jarak -0.2302 -0.1666 0.8538 Gagal tolak H0
k-NN (k = 2) 0.2567 -0.1692 0.1214 ✔ Gagal tolak H0
Queen Contiguity 0.0023 -0.1681 0.5724 Gagal tolak H0 
glob$perm %>%
  mutate(pembobot = factor(pembobot, levels = names(daftar_lw))) %>%
  ggplot(aes(Ip)) +
  geom_histogram(bins = 40, fill = "grey75", color = "white") +
  geom_vline(aes(xintercept = I_obs), color = "#c0392b", linewidth = 0.9) +
  facet_wrap(~ pembobot) +
  labs(title = "Distribusi permutasi Moran Bivariat dan nilai observasi (garis merah)",
       x = "Moran Bivariat (I)", y = "Frekuensi") +
  theme_minimal(base_size = 11)

b <- tabel_global[tabel_global$terbaik, ]
nama_pilihan <- b$pembobot
lw_pilihan   <- daftar_lw[[nama_pilihan]]
ada_signif   <- any(tabel_global$p_permutasi < 0.05)
fp <- function(p) ifelse(p < 0.001, "< 0,001", paste("=", rp(p, 3)))
alasan <- list(
  "Invers jarak" = "Secara substantif pembobot ini masuk akal karena memakai seluruh informasi jarak antarwilayah dan tidak ada wilayah yang lag-nya ditentukan oleh satu tetangga saja.",
  "k-NN (k = 2)" = "Secara substantif pembobot ini masuk akal karena setiap wilayah dihitung dari jumlah tetangga yang sama (dua), sehingga tidak ada wilayah yang lag-nya ditentukan satu tetangga saja seperti pada Queen (tiga kota kecil hanya punya satu tetangga). Kelemahannya, tetangga ditentukan oleh jarak titik pusat, sehingga kabupaten yang luas dapat terhubung dengan kota kecil yang ada di dalam atau di sekitarnya.",
  "Queen Contiguity" = "Secara substantif pembobot ini mengikuti batas administrasi yang nyata. Kelemahannya, tiga kota kecil hanya memiliki satu tetangga."
)
arah <- function(I) ifelse(I > 0, "positif", "negatif")
ket <- apply(tabel_global, 1, function(r) sprintf("%s (I = %s; p %s)", r["pembobot"],
        rp(as.numeric(r["I"]), 3), fp(as.numeric(r["p_permutasi"]))))

Pemilihan pembobot terbaik. Kriteria yang dipakai: pembobot dengan bukti autokorelasi terkuat, yaitu p-value permutasi terkecil (jika seri, nilai |I| terbesar). Hasilnya: Invers jarak (I = -0,230; p = 0,854); k-NN (k = 2) (I = 0,257; p = 0,121); Queen Contiguity (I = 0,002; p = 0,572). Pembobot terbaik adalah k-NN (k = 2) (I = 0,257, p = 0,121). Secara substantif pembobot ini masuk akal karena setiap wilayah dihitung dari jumlah tetangga yang sama (dua), sehingga tidak ada wilayah yang lag-nya ditentukan satu tetangga saja seperti pada Queen (tiga kota kecil hanya punya satu tetangga). Kelemahannya, tetangga ditentukan oleh jarak titik pusat, sehingga kabupaten yang luas dapat terhubung dengan kota kecil yang ada di dalam atau di sekitarnya.

Interpretasi.

  • Nilai Moran Bivariat dengan k-NN (k = 2) bernilai positif (I = 0,257), tetapi tidak signifikan pada taraf 5% (p = 0,121). Ketiga pembobot sama-sama gagal menolak H0.
  • Dengan demikian tidak ada bukti autokorelasi spasial bivariat: luas panen suatu kabupaten/kota tidak berhubungan secara sistematis dengan produksi padi di wilayah tetangganya. Besar dan arah I yang berubah-ubah antarpembobot (Invers jarak: negatif; k-NN (k = 2): positif; Queen Contiguity: positif) menunjukkan bahwa dugaan pola pada peta tidak stabil dan tidak cukup kuat untuk disebut klaster spasial.
  • Keputusan gagal menolak H0 berlaku untuk ketiga pembobot, sehingga kesimpulan tidak bergantung pada pilihan pembobot (yang berubah hanya arah dan besar I). Besarnya produksi di Banten lebih ditentukan oleh karakteristik tiap daerah (kabupaten agraris vs kota) daripada oleh pengaruh wilayah sekitarnya.
  • Karena korelasi X dan Y sangat tinggi (r = 0,998), Moran Bivariat di sini praktis setara dengan autokorelasi spasial univariat luas panen. Hasilnya tidak menunjukkan penularan spasial luas panen.
  • Catatan kekuatan uji. Dengan hanya 7 wilayah, jumlah susunan permutasi hanya 7! = 5.040, sehingga uji sangat sulit mendeteksi pola yang lemah. Hasil “tidak signifikan” berarti tidak ada bukti, bukan bukti bahwa pola spasial pasti tidak ada.

3.6 Local Moran Bivariat — LISA Bivariat

LISA Bivariat menunjukkan wilayah mana yang berkontribusi pada pola global. Kuadran ditentukan oleh \(z_x\) (luas panen wilayah itu) dan rata-rata tertimbang \(z_y\) (produksi) di wilayah tetangganya:

  • High-High (HH): luas panen tinggi, produksi tetangga tinggi.
  • Low-Low (LL): luas panen rendah, produksi tetangga rendah.
  • High-Low (HL): luas panen tinggi, produksi tetangga rendah.
  • Low-High (LH): luas panen rendah, produksi tetangga tinggi.

Pembobot yang dipakai adalah pembobot terbaik dari tahap global: k-NN (k = 2).

hitung_lisa <- function(x, y, lw, nama, nsim = 999) {
  l <- as.data.frame(spdep::localmoran_bv(x, y, lw, nsim = nsim, scale = TRUE,
                                          alternative = "two.sided", iseed = 123))
  # Nama kolom p-value berbeda antarversi spdep ("p_sim" pada versi baru,
  # "Pr(z != E(Ibvi)) Sim" pada versi lama). Keduanya adalah p-value permutasi dua arah.
  kol_p <- intersect(c("p_sim", "Pr(z != E(Ibvi)) Sim"), names(l))
  stopifnot(length(kol_p) == 1)
  zx     <- as.numeric(scale(x))
  lag_zy <- as.numeric(spdep::lag.listw(lw, as.numeric(scale(y))))
  data.frame(wilayah = nama, zx = zx, lag_zy = lag_zy,
             Ibvi = zx * lag_zy, p_perm = l[[kol_p]]) %>%
    mutate(kuadran = case_when(zx >= 0 & lag_zy >= 0 ~ "High-High",
                               zx <  0 & lag_zy <  0 ~ "Low-Low",
                               zx >= 0 & lag_zy <  0 ~ "High-Low",
                               TRUE ~ "Low-High"),
           p_BH   = p.adjust(p_perm, method = "BH"),
           status = ifelse(p_BH < 0.05, "Signifikan", "Tidak signifikan")) %>%
    select(wilayah, zx, lag_zy, Ibvi, p_perm, p_BH, kuadran, status)
}

lisa <- hitung_lisa(x, y, lw_pilihan, nama_w)
kable(lisa, digits = 3,
      col.names = c("Wilayah", "z(X)", "lag z(Y)", "Ii (bivariat)",
                    "p permutasi", "p BH", "Kuadran", "Status (BH 5%)"),
      caption = paste("LISA Bivariat, pembobot", nama_pilihan, "(999 permutasi, dua arah)"))
LISA Bivariat, pembobot k-NN (k = 2) (999 permutasi, dua arah)
Wilayah z(X) lag z(Y) Ii (bivariat) p permutasi p BH Kuadran Status (BH 5%)
Kabupaten Serang 0.679 -0.962 -0.654 0.372 0.521 High-Low Tidak signifikan
Kabupaten Tangerang 0.202 -0.979 -0.197 0.344 0.521 High-Low Tidak signifikan
Kota Cilegon -1.074 -0.057 0.061 0.704 0.740 Low-Low Tidak signifikan
Kota Serang -0.853 -0.173 0.148 0.740 0.740 Low-Low Tidak signifikan
Kota Tangerang -1.108 -0.304 0.336 0.332 0.521 Low-Low Tidak signifikan
Lebak 0.889 1.025 0.911 0.124 0.434 High-High Tidak signifikan
Pandeglang 1.265 0.739 0.935 0.114 0.434 High-High Tidak signifikan 
warna_kuadran <- c("High-High" = "#c0392b", "Low-Low" = "#2471a3",
                   "High-Low" = "#f0a35e", "Low-High" = "#8fc1e3")

# Moran scatterplot bivariat
ggplot(lisa, aes(zx, lag_zy, color = kuadran, label = sub("Kabupaten ", "Kab. ", wilayah))) +
  geom_hline(yintercept = 0, color = "grey50") + geom_vline(xintercept = 0, color = "grey50") +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE, color = "grey35",
              linewidth = 0.6, linetype = "dashed", formula = y ~ x) +
  geom_point(size = 3.5) +
  geom_text(vjust = -1, size = 3.2, show.legend = FALSE, check_overlap = TRUE) +
  scale_color_manual(values = warna_kuadran, name = "Kuadran") +
  labs(title = "Moran scatterplot bivariat",
       subtitle = paste("X = luas panen (standar), sumbu tegak = lag produksi (standar);", nama_pilihan),
       x = "Luas panen, z(X)", y = "Rata-rata tertimbang produksi tetangga, lag z(Y)") +
  theme_minimal(base_size = 12)

peta_lisa <- peta_8 %>%
  left_join(lisa %>% select(wilayah, kuadran, status, p_BH), by = "wilayah") %>%
  mutate(
    kuadran_peta = ifelse(is.na(kuadran), "Tidak ada data", kuadran),
    sig_peta = case_when(is.na(kuadran) ~ "Tidak ada data",
                         status == "Signifikan" ~ kuadran,
                         TRUE ~ "Tidak signifikan"))

pal_peta <- c(warna_kuadran, "Tidak signifikan" = "grey80", "Tidak ada data" = "white")

peta_kelas <- function(kolom, judul, sub) {
  ggplot(peta_lisa) +
    geom_sf(aes(fill = .data[[kolom]]), color = "grey30", linewidth = 0.3) +
    geom_label(data = label_peta, aes(lx, ly, label = nama), size = 2.6, fontface = "bold",
               fill = alpha("white", 0.75), label.size = 0, label.padding = unit(0.12, "lines")) +
    scale_fill_manual(values = pal_peta, name = "Kategori", drop = TRUE) +
    labs(title = judul, subtitle = sub, x = NULL, y = NULL) +
    theme_minimal(base_size = 11) + theme(panel.grid = element_blank(), axis.text = element_blank())
}
peta_kelas("kuadran_peta", "Peta kuadran LISA Bivariat (tanpa uji signifikansi)",
           paste("Luas panen (X) vs produksi tetangga (Y); pembobot:", nama_pilihan))

peta_kelas("sig_peta", "Peta signifikansi LISA Bivariat",
           paste("Signifikan bila p BH < 0,05; pembobot:", nama_pilihan))

n_sig <- sum(lisa$status == "Signifikan")
kel_kuadran <- lisa %>% group_by(kuadran) %>% summarise(w = paste(sub("Kabupaten ", "Kab. ", wilayah), collapse = ", "), .groups = "drop")
tulis_kuadran <- function(k) { v <- kel_kuadran$w[kel_kuadran$kuadran == k]; ifelse(length(v) == 0, "tidak ada", v) }
p_min_lokal <- min(lisa$p_perm)
w_p_min <- lisa$wilayah[which.min(lisa$p_perm)]

Interpretasi.

  • Signifikansi. Setelah koreksi Benjamini–Hochberg, 0 dari 7 wilayah memiliki LISA signifikan. p-value permutasi terkecil adalah 0,114 (Pandeglang), yang bahkan tanpa koreksi belum mencapai 0,05. Tidak ada klaster (HH/LL) maupun pencilan spasial (HL/LH) yang dapat dipastikan secara statistik. Seluruh wilayah pada peta signifikansi berwarna abu-abu.
  • Arah pola lokal (deskriptif saja). Menurut kuadran Moran scatterplot, wilayah High-High: Lebak, Pandeglang; Low-Low: Kota Cilegon, Kota Serang, Kota Tangerang; High-Low: Kab. Serang, Kab. Tangerang; Low-High: tidak ada. Kuadran ini hanya menggambarkan arah hubungan lokal dan tidak boleh dibaca sebagai klaster karena tidak signifikan.
  • Pola lokal. Pandeglang (Ii = 0,93) dan Lebak (Ii = 0,91) adalah dua wilayah berluas panen terbesar yang saling bertetangga pada pembobot terpilih, sehingga keduanya membentuk blok High-High di selatan-barat dengan Ii terbesar. Kota Cilegon, Kota Serang, dan Kota Tangerang berada di sisi Low-Low. Kabupaten Serang dan Kabupaten Tangerang berluas panen besar, tetapi tetangga terdekat menurut titik pusat adalah kota-kota kecil berproduksi rendah, sehingga keduanya High-Low dengan Ii negatif. Kontras High-Low ini mengurangi kekuatan hubungan global.
  • Konsisten dengan hasil global. Tidak signifikannya seluruh wilayah sejalan dengan Moran Bivariat Global yang tidak signifikan.

3.7 Analisis sensitivitas: luas panen dan produktivitas

Pasangan variabel awal (luas panen dan produktivitas) dianalisis dengan prosedur yang sama untuk memeriksa apakah kesimpulan berubah.

glob_p <- hitung_global(x, yp, daftar_lw)
kable(glob_p$tabel %>% mutate(terbaik = ifelse(terbaik, "✔", "")),
      digits = 4, col.names = c("Pembobot", "Moran Bivariat (I)", "E[I] permutasi",
                                "p-value permutasi", "Terbaik", "Keputusan (α = 5%)"),
      caption = "Sensitivitas: Moran Bivariat Global luas panen (X) dan produktivitas (Y)")
Sensitivitas: Moran Bivariat Global luas panen (X) dan produktivitas (Y)
Pembobot Moran Bivariat (I) E[I] permutasi p-value permutasi Terbaik Keputusan (α = 5%)
Invers jarak 0.0546 0.0844 0.8273 Gagal tolak H0
k-NN (k = 2) -0.0931 0.0855 0.4743 Gagal tolak H0
Queen Contiguity -0.1511 0.0847 0.2980 ✔ Gagal tolak H0 
nama_pilihan_p <- glob_p$tabel$pembobot[glob_p$tabel$terbaik]
lisa_p <- hitung_lisa(x, yp, daftar_lw[[nama_pilihan_p]], nama_w)
n_sig_p <- sum(lisa_p$status == "Signifikan")

Pada pasangan luas panen dan produktivitas, tidak satu pun pembobot menghasilkan Moran Bivariat Global yang signifikan (semua p-value > 0,05), dan LISA dengan pembobot terbaik (Queen Contiguity) memiliki 0 wilayah signifikan. Kesimpulan tidak berubah: tidak ditemukan autokorelasi spasial bivariat pada data padi Banten 2025, baik untuk produksi maupun produktivitas.

4 Kesimpulan

4.1 Jawaban atas rumusan masalah

kable(tibble(
  `Rumusan masalah` = c("1. Gambaran dan korelasi", "2. Sebaran spasial", "3. Moran Bivariat Global dan pembobot terbaik", "4. LISA Bivariat"),
  Jawaban = c(
    sprintf("Luas panen sangat beragam (KV %s%%; %s ha di %s sampai %s ha di %s), produktivitas relatif seragam (KV %s%%). Luas panen dan produksi berkorelasi r = %s (signifikan).",
            rp(kv_luas, 1), rp(tabel_desk[1, "Min"]), w_min_luas, rp(tabel_desk[1, "Maks"]), w_max_luas, rp(kv_pvt, 1), rp(uji_xy$estimate, 3)),
    "Luas panen dan produksi tinggi di kabupaten agraris (Pandeglang, Lebak, Kab. Serang, Kab. Tangerang) dan rendah di kota-kota utara.",
    sprintf("Tidak ada autokorelasi spasial bivariat yang signifikan pada ketiga pembobot. Pembobot terbaik: %s (I = %s; p %s).", nama_pilihan, rp(b$I, 3), fp(b$p_permutasi)),
    sprintf("%s dari 7 wilayah signifikan (BH 5%%): tidak ada klaster atau pencilan spasial yang terbukti secara statistik.", n_sig))
), caption = "Ringkasan jawaban")
Ringkasan jawaban
Rumusan masalah Jawaban
1. Gambaran dan korelasi Luas panen sangat beragam (KV 90,0%; 164 ha di Kota Tangerang sampai 105.490 ha di Pandeglang), produktivitas relatif seragam (KV 8,5%). Luas panen dan produksi berkorelasi r = 0,998 (signifikan).
2. Sebaran spasial Luas panen dan produksi tinggi di kabupaten agraris (Pandeglang, Lebak, Kab. Serang, Kab. Tangerang) dan rendah di kota-kota utara.
3. Moran Bivariat Global dan pembobot terbaik Tidak ada autokorelasi spasial bivariat yang signifikan pada ketiga pembobot. Pembobot terbaik: k-NN (k = 2) (I = 0,257; p = 0,121).
4. LISA Bivariat 0 dari 7 wilayah signifikan (BH 5%): tidak ada klaster atau pencilan spasial yang terbukti secara statistik.

4.2 Kesimpulan rinci

  1. Data. Dua variabel yang dipakai adalah luas panen (ha) dan produksi padi (ton) kabupaten/kota di Provinsi Banten 2025 dari BPS, dengan korelasi Pearson r = 0,998 (signifikan). Analisis memakai 7 wilayah karena Kota Tangerang Selatan tidak memiliki data.
  2. Deskriptif. Luas panen sangat beragam (koefisien variasi 90,0%, dari 164 ha di Kota Tangerang sampai 105.490 ha di Pandeglang), sedangkan produktivitas relatif seragam (koefisien variasi 8,5%). Perbedaan produksi antardaerah terutama dipengaruhi skala lahan.
  3. Peta. Luas panen dan produksi terkonsentrasi di kabupaten agraris (Pandeglang, Lebak, Kabupaten Serang, Kabupaten Tangerang), sedangkan kota-kota di utara rendah.
  4. Pembobot. Invers jarak, k-NN (k = 2), dan Queen Contiguity dibentuk dan dibakukan per baris. Moran Bivariat Global bernilai Invers jarak: -0,230; k-NN (k = 2): 0,257; Queen Contiguity: 0,002, seluruhnya tidak signifikan. Pembobot terbaik menurut kriteria p-value terkecil adalah k-NN (k = 2).
  5. LISA Bivariat. Dengan k-NN (k = 2), 0 dari 7 wilayah signifikan (BH 5%). Tidak ada klaster HH/LL maupun pencilan HL/LH yang terbukti secara statistik.
  6. Makna. Tidak ditemukan bukti bahwa luas panen suatu wilayah di Banten terkait dengan produksi padi wilayah tetangganya. Produksi lebih ditentukan oleh ciri tiap daerah (agraris vs perkotaan). Hasil yang sama diperoleh pada analisis sensitivitas dengan produktivitas.
  7. Keterbatasan. (a) Hanya 7 wilayah, sehingga kekuatan uji rendah dan hasil “tidak signifikan” tidak membuktikan ketiadaan pola. (b) Kota Tangerang Selatan dikeluarkan sehingga struktur ketetanggaan berubah. (c) Data tahun tunggal tingkat kabupaten/kota. (d) Besar dan arah Moran Bivariat berubah antarpembobot dan antara titik pusat berbeda, sehingga estimasinya tidak stabil pada n sekecil ini. (e) Luas panen dan produksi hampir kolinear, sehingga Moran Bivariat praktis setara dengan autokorelasi luas panen. (f) Analisis bersifat eksploratif dan tidak menyatakan hubungan sebab-akibat. Penelitian lanjutan dapat memakai data tingkat kecamatan agar jumlah unit lebih besar.

5 Daftar pustaka

Tanggal akses seluruh tautan: 27 Oktober 2026. DOI dapat dibuka melalui https://doi.org/<DOI>.

Sumber data

Teori dan metode

Perangkat lunak

Bahan kuliah

  • Sholihin, M., & Sari, P. D. (2026). Uji kehomogenan dan autokorelasi spasial univariat (Pertemuan 4, Statistika Spasial STL623303, Semester Gasal 2026/2027) [Slide kuliah]. Program Studi Statistika, Universitas Sultan Ageng Tirtayasa.
sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] scales_1.4.0  knitr_1.51    readxl_1.4.5  dplyr_1.2.0   spdep_1.4-2  
## [6] spData_2.3.5  sf_1.1-0      terra_1.9-50  ggplot2_4.0.2
## 
## 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       cellranger_1.1.0   jsonlite_2.0.0     e1071_1.7-17      
## [17] DBI_1.3.0          mgcv_1.9-3         viridisLite_0.4.3  codetools_0.2-20  
## [21] jquerylib_0.1.4    cli_3.6.5          rlang_1.1.7        units_1.0-0       
## [25] splines_4.5.2      withr_3.0.2        cachem_1.1.0       yaml_2.3.12       
## [29] otel_0.2.0         tools_4.5.2        deldir_2.0-4       boot_1.3-32       
## [33] vctrs_0.7.1        R6_2.6.1           proxy_0.4-29       lifecycle_1.0.5   
## [37] classInt_0.4-11    dbscan_1.2.4       pkgconfig_2.0.3    pillar_1.11.1     
## [41] bslib_0.10.0       gtable_0.3.6       glue_1.8.0         Rcpp_1.1.1        
## [45] xfun_0.61          tibble_3.3.1       tidyselect_1.2.1   rstudioapi_0.18.0 
## [49] farver_2.1.2       igraph_2.2.3       nlme_3.1-168       htmltools_0.5.9   
## [53] rmarkdown_2.30     labeling_0.4.3     wk_0.9.5           compiler_4.5.2    
## [57] S7_0.2.1           sp_2.2-3