1 Pendahuluan

1.1 Latar Belakang

Beras adalah komoditas pangan pokok masyarakat Indonesia, sehingga harganya menjadi indikator penting kesejahteraan dan stabilitas ekonomi daerah. Harga eceran beras tidak seragam di seluruh wilayah: ia dipengaruhi oleh lokasi sentra produksi, biaya transportasi dan distribusi, struktur pasar lokal, serta kondisi geografis (misalnya wilayah kepulauan). Karena faktor-faktor tersebut bersifat kewilayahan, harga beras di kabupaten/kota yang saling berdekatan diduga cenderung mirip.

Dugaan ini sejalan dengan Hukum I Tobler: segala sesuatu saling berhubungan, tetapi lokasi yang berdekatan lebih saling berhubungan dibandingkan lokasi yang berjauhan. Dalam statistika spasial, keterkaitan tersebut disebut autokorelasi spasial dan dapat diukur serta diuji secara statistik, misalnya dengan Indeks Moran, Geary’s C, dan Local Indicators of Spatial Association (LISA).

1.2 Rumusan Masalah

  1. Apakah ragam harga eceran beras homogen antar zona wilayah di Pulau Sumatera?
  2. Apakah terdapat autokorelasi spasial pada harga eceran beras kabupaten/kota di Pulau Sumatera, dan seberapa kuat?
  3. Di mana letak klaster harga tinggi (hot spot), klaster harga rendah (cold spot), dan pencilan spasial?

1.3 Tujuan

  1. Menguji kehomogenan ragam harga beras antar zona lintang menggunakan uji Bartlett dan Levene.
  2. Mengukur dan menguji autokorelasi spasial global dengan Moran’s I dan Geary’s C.
  3. Mengidentifikasi dan memetakan klaster lokal (High-High, Low-Low, High-Low, Low-High) dengan Local Moran (LISA).

2 Pedoman Interpretasi

Berikut pedoman membaca hasil uji pada bagian berikutnya.

  • Uji kehomogenan (Bartlett dan Levene). \(H_0\): ragam sama di semua zona. Jika p < 0,05, \(H_0\) ditolak sehingga ragam tidak homogen (heteroskedastisitas spasial) dan interpretasi autokorelasi perlu dilakukan dengan hati-hati. Levene (Brown–Forsythe) lebih tahan terhadap pelanggaran normalitas dibanding Bartlett.
  • Moran’s I. Dibandingkan dengan nilai harapan \(E(I)=-1/(n-1)\approx 0\). \(I\) jauh di atas \(E(I)\) berarti autokorelasi positif (nilai serupa mengelompok), \(I\approx E(I)\) berarti pola acak, dan \(I\) jauh di bawah \(E(I)\) berarti autokorelasi negatif (nilai berlawanan berselang-seling seperti papan catur). Signifikansi dinilai dari Z dan p-value (\(H_0\): tidak ada autokorelasi).
  • Geary’s C. Dibandingkan dengan \(E(C)=1\). Arahnya kebalikan Moran: \(C<1\) berarti autokorelasi positif, \(C=1\) acak, dan \(C>1\) negatif. Kedua indeks dipakai berdampingan sebagai pemeriksaan silang.
  • Local Moran (LISA). Menunjukkan di mana klaster berada. High-High: nilai tinggi dikelilingi tetangga bernilai tinggi (hot spot). Low-Low: nilai rendah dikelilingi tetangga rendah (cold spot). High-Low dan Low-High: pencilan spasial. Klaster dianggap signifikan jika p < 0,05.

Cara menjalankan: letakkan file .Rmd ini di folder yang sama dengan file data CSV dan folder SHP Sumatera KabKot. File dicari otomatis (folder .Rmd, lalu ~/Documents/intan, ~/Documents, ~/Desktop, ~/Downloads), jadi tidak perlu setwd(). Folder shapefile harus utuh (.shp, .shx, .dbf, .prj). Lalu klik Knit → Knit to HTML.

3 Hasil dan Pembahasan

3.1 Impor dan Eksplorasi Data

# Pencarian file otomatis. Urutan folder yang dicari (rekursif, berhenti di temuan pertama):
# folder .Rmd, lalu lokasi umum. Tambahkan folder Anda sendiri di sini bila perlu.
lokasi_tambahan <- c("~/Documents/intan", "~/Documents", "~/Desktop", "~/Downloads")
cari_dirs <- unique(c(".", lokasi_tambahan))
cari_dirs <- cari_dirs[dir.exists(cari_dirs)]

cari <- function(pola, wajib = TRUE) {
  for (dr in cari_dirs) {
    f <- list.files(dr, pattern = pola, recursive = TRUE,
                    full.names = TRUE, ignore.case = TRUE)
    if (length(f) > 0) return(sort(f)[1])
  }
  if (wajib)
    stop("File dengan pola '", pola, "' tidak ditemukan. Folder kerja: ", getwd(),
         ". Letakkan file data di folder yang sama dengan file .Rmd.")
  NA_character_
}
f_beras <- cari("^Beras Pulau Sumatera 2025.*\\.csv$")
f_shp   <- cari("^Sumatera-KabKot\\.shp$", wajib = FALSE)
ada_shp <- !is.na(f_shp)
cat("File CSV :", f_beras, "\n")
## File CSV : ./Beras Pulau Sumatera 2025.csv
cat("File SHP :", if (ada_shp) f_shp else "TIDAK DITEMUKAN (peta shapefile dilewati)", "\n")
## File SHP : /Users/user/Documents/intan/SP/SHP Sumatera KabKot/Sumatera-KabKot.shp
d <- read.csv(f_beras, sep = ";", stringsAsFactors = FALSE, check.names = FALSE,
              strip.white = TRUE, fileEncoding = "UTF-8-BOM")
names(d) <- trimws(names(d))
names(d)[names(d) %in% c("Haga Beras", "Harga Beras", "Haga.Beras")] <- "Harga_Beras"
d$KabKot      <- trimws(d$KabKot)
d$ID          <- as.integer(d$ID)
d$Long        <- as.numeric(d$Long)
d$Lat         <- as.numeric(d$Lat)
d$Harga_Beras <- as.numeric(d$Harga_Beras)
d <- d[!is.na(d$Harga_Beras), ]

head(d)
##    ID          KabKot  Long  Lat Harga_Beras
## 1 293      Aceh Barat 96.13 4.14       11333
## 2 277 Aceh Barat Daya 96.83 3.74       12875
## 3 278      Aceh Besar 95.35 5.51       10000
## 4 279       Aceh Jaya 95.79 4.47       11375
## 5 405      Aceh Pidie 95.95 5.38       11000
## 6 365    Aceh Selatan 97.18 3.26       13000
str(d)
## 'data.frame':    143 obs. of  5 variables:
##  $ ID         : int  293 277 278 279 405 365 415 280 281 282 ...
##  $ KabKot     : chr  "Aceh Barat" "Aceh Barat Daya" "Aceh Besar" "Aceh Jaya" ...
##  $ Long       : num  96.1 96.8 95.3 95.8 96 ...
##  $ Lat        : num  4.14 3.74 5.51 4.47 5.38 3.26 2.4 4.29 4.62 3.49 ...
##  $ Harga_Beras: num  11333 12875 10000 11375 11000 ...
summary(d$Harga_Beras)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    8875   10500   11044   11392   12106   16000

Data memuat 143 kabupaten/kota dengan harga beras berkisar Rp 8.875 hingga Rp 16.000 per kg (median Rp 11.044; rata-rata Rp 11.392).

3.2 Uji Kehomogenan Ragam

terc <- unname(quantile(d$Lat, probs = c(1/3, 2/3)))
d$Zona <- cut(d$Lat, breaks = c(-Inf, terc[1], terc[2], Inf),
              labels = c("Selatan", "Tengah", "Utara"))
table(d$Zona)
## 
## Selatan  Tengah   Utara 
##      48      47      48

Bartlett’s Test (asumsi data tiap kelompok menyebar normal):

bt <- bartlett.test(Harga_Beras ~ Zona, data = d)
bt
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Harga_Beras by Zona
## Bartlett's K-squared = 13.193, df = 2, p-value = 0.001365

Levene’s Test versi Brown–Forsythe (berbasis median, lebih robust terhadap pelanggaran normalitas):

levene_test <- function(y, g) {
  gm <- ave(y, g, FUN = median)
  z  <- abs(y - gm)
  summary(aov(z ~ g))
}
levene_test(d$Harga_Beras, d$Zona)
##              Df   Sum Sq Mean Sq F value Pr(>F)  
## g             2  3754211 1877105   4.571 0.0119 *
## Residuals   140 57486026  410614                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pembahasan. Bartlett: \(K^2\) = 13.193, db = 2, p 0.0014. Levene: F = 4.571, p 0.0119. Karena kedua uji signifikan (p < 0,05), sehingga \(H_0\) ditolak: ragam harga beras tidak homogen antar zona lintang. Ragam terbesar terdapat pada zona Utara. Hal ini wajar untuk data harga yang memang bervariasi menurut karakteristik wilayah (kepadatan penduduk, jalur distribusi, dan sebagainya). Ketidakhomogenan ini perlu diingat saat menafsirkan autokorelasi spasial.

3.3 Matriks Pembobot Spasial

koordinat <- cbind(d$Long, d$Lat)
nb <- knn2nb(knearneigh(koordinat, k = 6), row.names = d$KabKot)
W  <- nb2listw(nb, style = "W")
nb
## Neighbour list object:
## Number of regions: 143 
## Number of nonzero links: 858 
## Percentage nonzero weights: 4.195804 
## Average number of links: 6 
## Non-symmetric neighbours list

3.4 Autokorelasi Spasial Global

Indeks Moran’s I:

moran_res <- moran.test(d$Harga_Beras, W, zero.policy = TRUE)
moran_res
## 
##  Moran I test under randomisation
## 
## data:  d$Harga_Beras  
## weights: W    
## 
## Moran I statistic standard deviate = 13.173, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.572974950      -0.007042254       0.001938624

Uji permutasi (Monte Carlo) sebagai pembanding uji analitik:

set.seed(123)
moran_mc <- moran.mc(d$Harga_Beras, W, nsim = 999, zero.policy = TRUE)
moran_mc
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  d$Harga_Beras 
## weights: W  
## number of simulations + 1: 1000 
## 
## statistic = 0.57297, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Indeks Geary’s C:

geary_res <- geary.test(d$Harga_Beras, W, zero.policy = TRUE)
geary_res
## 
##  Geary C test under randomisation
## 
## data:  d$Harga_Beras 
## weights: W   
## 
## Geary C statistic standard deviate = 12.566, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.379248094       1.000000000       0.002440147
Ringkasan autokorelasi spasial global (k-NN, k = 6)
Indeks Nilai Nilai_Harapan Z p_value
Moran’s I 0.5730 -0.007 13.1733 0
Geary’s C 0.3792 1.000 12.5664 0

Pembahasan. Moran’s I = 0.573 jauh di atas nilai harapannya \(E(I)\) = -0.007 (Z = 13.17, p < 0.001; p uji permutasi = 0.001), dan Geary’s C = 0.379 berada di bawah 1 (Z = 12.57, p < 0.001). Kedua indeks konsisten dan menunjukkan autokorelasi spasial positif yang signifikan: harga beras di kab/kota yang berdekatan cenderung mirip (mengelompok).

3.5 Moran Scatterplot

z  <- scale(d$Harga_Beras)[, 1]
wz <- lag.listw(W, z, zero.policy = TRUE)

moran.plot(d$Harga_Beras, W, labels = d$KabKot, zero.policy = TRUE,
           xlab = "Harga Beras (Rp/kg)", ylab = "Rata-rata Tetangga (Spatial Lag)")

Kemiringan garis regresi pada scatterplot secara visual sebanding dengan nilai Moran’s I (0.573). Titik yang banyak berada di kuadran kanan-atas (High-High) dan kiri-bawah (Low-Low) mengindikasikan autokorelasi positif.

3.6 Local Moran (LISA)

lisa <- localmoran(d$Harga_Beras, W, zero.policy = TRUE)
d$Ii  <- lisa[, 1]   # nilai I_i
d$ZIi <- lisa[, 4]   # nilai Z
d$PIi <- lisa[, 5]   # p-value (analitik, tergantung argumen)

# klasifikasi kuadran berdasarkan z_i dan (Wz)_i
d$quadrant <- ifelse(z >= 0 & wz >= 0, "High-High",
               ifelse(z <  0 & wz <  0, "Low-Low",
               ifelse(z <  0 & wz >= 0, "Low-High", "High-Low")))
d$sig     <- d$PIi < 0.05
d$cluster <- ifelse(d$sig, d$quadrant, "Not Significant")

table(d$cluster)
## 
##       High-High        High-Low        Low-High         Low-Low Not Significant 
##              29               3               2              23              86

Wilayah dengan klaster signifikan, diurutkan berdasarkan kekuatan \(I_i\):

tab_sig <- d[d$sig, c("KabKot", "Harga_Beras", "Ii", "PIi", "cluster")]
tab_sig <- tab_sig[order(-tab_sig$Ii), ]
rownames(tab_sig) <- NULL
knitr::kable(tab_sig, digits = 4)
KabKot Harga_Beras Ii PIi cluster
Kep. Anambas 16000.00 4.0271 0.0045 High-High
Kota Padangpanjang 13812.50 2.7723 0.0004 High-High
Kota Bukittinggi 13500.00 2.6101 0.0002 High-High
Bengkalis 13750.00 2.5247 0.0010 High-High
Kuantan Sengingi 13875.00 2.4537 0.0023 High-High
Kota Dumai 13500.00 2.4377 0.0004 High-High
Kota Pekanbaru 13499.00 2.3766 0.0006 High-High
Lima Puluh Kota 13500.00 2.3390 0.0007 High-High
Agam 13500.00 2.2186 0.0013 High-High
Tanah Datar 13000.00 2.1517 0.0000 High-High
Kota Payahkumbuh 13125.00 1.9972 0.0004 High-High
Kota Solok 13500.00 1.9446 0.0047 High-High
Kota Pariaman 13000.00 1.8049 0.0006 High-High
Mesuji 9000.00 1.7372 0.0247 Low-Low
Kota Tanjung Pinang 13500.00 1.7217 0.0121 High-High
Palelawan 12875.00 1.5030 0.0021 High-High
Sijunjung 13375.00 1.4780 0.0222 High-High
Karimun 13000.00 1.4691 0.0054 High-High
Pasawaran 9750.00 1.4249 0.0082 Low-Low
Kota Sawahlunto 12625.00 1.4033 0.0006 High-High
OKU Timur 9812.50 1.3548 0.0090 Low-Low
Kampar 12500.00 1.2846 0.0005 High-High
Kepulauan Meranti 12500.00 1.2495 0.0006 High-High
Pringsewu 10000.00 1.2477 0.0065 Low-Low
Kota Batam 12750.00 1.2255 0.0061 High-High
Tanggamus 10000.00 1.2199 0.0078 Low-Low
Pesisir Selatan 13000.00 1.1961 0.0232 High-High
Lampung Barat 10000.00 1.1801 0.0100 Low-Low
Siak 12348.00 1.1220 0.0004 High-High
Kota Metro 10000.00 1.0886 0.0174 Low-Low
Bangka Selatan 10000.00 1.0572 0.0209 Low-Low
Pesisir Barat 10175.00 1.0561 0.0085 Low-Low
Padang Pariaman 12250.00 1.0364 0.0003 High-High
Lampung Utara 10000.00 1.0210 0.0256 Low-Low
Kota Bandar Lampung 10400.00 0.8779 0.0074 Low-Low
Lampung Tengah 10362.50 0.8478 0.0127 Low-Low
Bintan 12250.00 0.8232 0.0038 High-High
Belitung 10300.00 0.8169 0.0234 Low-Low
Indragiri Hulu 12500.00 0.8036 0.0278 High-High
Lampung Selatan 10500.00 0.7996 0.0068 Low-Low
Kota Padang 12112.50 0.7170 0.0027 High-High
Pasaman 12000.00 0.7005 0.0005 High-High
Kota Pagaralam 10500.00 0.6819 0.0208 Low-Low
Bangka Barat 10600.00 0.5450 0.0376 Low-Low
Way Kanan 10825.00 0.5286 0.0050 Low-Low
Empat Lawang 10666.67 0.4993 0.0377 Low-Low
Kota Lubuklinggau 10708.25 0.4837 0.0328 Low-Low
Tulangbawang 11000.00 0.3773 0.0038 Low-Low
Belitung Timur 11000.00 0.3247 0.0127 Low-Low
Ogan Komering Ulu (OKU) 11150.00 0.2275 0.0048 Low-Low
Kaur 11250.00 0.1431 0.0025 Low-Low
Solok 11500.00 0.1140 0.0015 High-High
Penukal Abab 11666.67 -0.2032 0.0265 High-Low
Dharmasraya 11125.00 -0.2664 0.0028 Low-High
Bengkulu Utara 12000.00 -0.4137 0.0418 High-Low
Aceh Tenggara 10700.00 -0.4610 0.0465 Low-High
Bangka Tengah 12000.00 -0.5310 0.0089 High-Low

Pembahasan. Pada p < 0,05 terdapat 29 wilayah High-High, 23 Low-Low, 3 High-Low, dan 2 Low-High; 86 wilayah lainnya tidak signifikan. High-High terkuat: Kep. Anambas, Kota Padangpanjang, Kota Bukittinggi, Bengkalis, Kuantan Sengingi, Kota Dumai. Low-Low terkuat: Mesuji, Pasawaran, OKU Timur, Pringsewu, Tanggamus, Lampung Barat.

Catatan: p-value di atas berasal dari pendekatan analitik localmoran(). Materi menggunakan uji permutasi, sehingga sebagai pembanding (jika spdep ≥ 1.3 terpasang):

if ("localmoran_perm" %in% getNamespaceExports("spdep")) {
  set.seed(123)
  lp   <- localmoran_perm(d$Harga_Beras, W, nsim = 999, zero.policy = TRUE)
  pcol <- grep("Sim", colnames(lp))[1]
  if (is.na(pcol)) pcol <- 5
  klaster_perm <- ifelse(lp[, pcol] < 0.05, d$quadrant, "Not Significant")
  print(table(klaster_perm))
} else {
  cat("localmoran_perm() tidak tersedia pada versi spdep ini; bagian pembanding dilewati.\n")
}
## klaster_perm
##       High-High        High-Low        Low-High         Low-Low Not Significant 
##              28               3               1              25              86

3.7 Visualisasi dan Pemetaan

cols <- c("High-High" = "#C0392B", "Low-Low" = "#2E86C1",
          "High-Low" = "#F1948A", "Low-High" = "#5DADE2",
          "Not Significant" = "#BFC9CA")

3.7.1 Peta titik

ggplot(d, aes(x = Long, y = Lat)) +
  geom_point(aes(size = Harga_Beras, color = Harga_Beras), alpha = 0.85) +
  scale_color_gradient(low = "#1C7293", high = "#C0392B") +
  labs(title = "Sebaran Harga Eceran Beras Kab/Kota Pulau Sumatera (Des 2025)",
       x = "Bujur", y = "Lintang", color = "Harga (Rp)", size = "Harga (Rp)") +
  theme_minimal()

ggplot(d, aes(x = Long, y = Lat, color = cluster)) +
  geom_point(size = 2.6, alpha = 0.9) +
  scale_color_manual(values = cols) +
  labs(title = "Peta Klaster LISA - Harga Beras Pulau Sumatera",
       x = "Bujur", y = "Lintang", color = "Klasifikasi") +
  theme_minimal()

3.7.2 Peta dengan shapefile (interaktif)

shp <- st_read(f_shp, quiet = TRUE)

# Bersihkan geometri agar bisa digambar oleh ggplot2/geom_sf:
# buang dimensi Z/M (campuran XY dan XYZ menyebabkan error "number of columns of
# matrices must match"), perbaiki geometri tidak valid, buang geometri kosong,
# dan seragamkan tipe menjadi MULTIPOLYGON.
bersihkan_geom <- function(x) {
  x <- sf::st_zm(x, drop = TRUE, what = "ZM")
  x <- sf::st_make_valid(x)
  x <- x[!sf::st_is_empty(x), ]
  if (any(as.character(sf::st_geometry_type(x)) == "GEOMETRYCOLLECTION")) {
    x <- sf::st_collection_extract(x, "POLYGON")
    x <- x[!sf::st_is_empty(x), ]
  }
  sf::st_cast(x, "MULTIPOLYGON")
}
shp <- bersihkan_geom(shp)

# gabungkan data spasial dan data harga beras
shp_data <- merge(shp, d, by = "ID")

# perbaiki geometri yang tidak valid, buang geometri kosong, pastikan objek sf
shp_data <- sf::st_make_valid(shp_data)
shp_data <- shp_data[!sf::st_is_empty(shp_data), ]
shp_data <- sf::st_as_sf(shp_data)
nrow(shp_data)
## [1] 166

Catatan: kolom ID tidak unik (SHP: 155 poligon dengan 142 ID; data: 143 baris dengan 136 ID), sehingga merge(..., by = "ID") menghasilkan 166 baris dengan sebagian poligon ganda yang saling bertumpuk. Karena itu peta statis yang sudah dikoreksi disajikan pada subbagian berikutnya.

library(mapview)
mapview(shp_data, zcol = "Harga_Beras")

Pemetaan kuadran Local Moran dan signifikansi autokorelasi spasial:

library(tmap)
tmap_mode("view")
tm_shape(shp_data) +
  tm_polygons("quadrant", fill.legend = tm_legend(
    title = "kuadran", show = TRUE)) +
  tm_layout(legend.outside = TRUE)
tm_shape(shp_data) +
  tm_polygons("cluster", fill.legend = tm_legend(
    title = "signifikansi", show = TRUE)) +
  tm_layout(legend.outside = TRUE)
library(tmap)
tmap_mode("view")
tm_shape(shp_data) + tm_polygons("quadrant", title = "kuadran")
tm_shape(shp_data) + tm_polygons("cluster", title = "signifikansi")

3.7.3 Peta statis dengan pencocokan poligon yang dikoreksi

Setiap baris data dicocokkan ke satu poligon yang ber-ID sama dan namanya paling mirip, sehingga tidak ada poligon ganda.

nm_shp <- toupper(gsub("[^A-Za-z]", "", as.character(shp$Kabupaten_)))
nm_d   <- toupper(gsub("[^A-Za-z]", "", d$KabKot))

idx_shp <- vapply(seq_len(nrow(d)), function(i) {
  cand <- which(shp$ID == d$ID[i])
  if (length(cand) == 0) return(NA_integer_)
  as.integer(cand[which.min(adist(nm_d[i], nm_shp[cand]))])
}, integer(1))

ok <- !is.na(idx_shp)
cat("Wilayah tanpa poligon di SHP:", paste(d$KabKot[!ok], collapse = ", "), "\n")
## Wilayah tanpa poligon di SHP: Pasaman
d_peta <- d[ok, ]
d_peta$Provinsi <- trimws(as.character(shp$Provinsi[idx_shp[ok]]))
d_peta$Provinsi[grepl("Aceh", d_peta$Provinsi, ignore.case = TRUE)] <- "Aceh"
d_peta$cluster_f  <- factor(d_peta$cluster,  levels = names(cols))
d_peta$quadrant_f <- factor(d_peta$quadrant, levels = names(cols))
shp_fix <- st_sf(d_peta, geometry = st_geometry(shp)[idx_shp[ok]])
# --- versi poligon (geom_sf)
gg_poligon_kat <- function(var, judul) {
  ggplot() +
    geom_sf(data = shp, fill = "grey95", color = "grey75") +
    geom_sf(data = shp_fix, aes(fill = .data[[var]]), color = "grey30") +
    scale_fill_manual(values = cols, drop = FALSE, name = "Klasifikasi") +
    labs(title = judul) + theme_minimal()
}
gg_poligon_num <- function(var, judul) {
  ggplot() +
    geom_sf(data = shp, fill = "grey95", color = "grey75") +
    geom_sf(data = shp_fix, aes(fill = .data[[var]]), color = "grey30") +
    scale_fill_gradient(low = "#FFF3B0", high = "#C0392B", name = "Rp/kg") +
    labs(title = judul) + theme_minimal()
}
# --- versi titik (cadangan bila poligon gagal digambar)
gg_titik_kat <- function(var, judul) {
  ggplot(d_peta, aes(x = Long, y = Lat, color = .data[[var]])) +
    geom_point(size = 2.6, alpha = 0.9) +
    scale_color_manual(values = cols, drop = FALSE, name = "Klasifikasi") +
    labs(title = judul, x = "Bujur", y = "Lintang") + theme_minimal()
}
gg_titik_num <- function(var, judul) {
  ggplot(d_peta, aes(x = Long, y = Lat, color = .data[[var]])) +
    geom_point(size = 2.6, alpha = 0.9) +
    scale_color_gradient(low = "#FFF3B0", high = "#C0392B", name = "Rp/kg") +
    labs(title = judul, x = "Bujur", y = "Lintang") + theme_minimal()
}
# --- menggambar dengan pengaman
peta_aman <- function(f_poligon, f_titik) {
  berhasil <- tryCatch({ print(f_poligon()); TRUE }, error = function(e) FALSE)
  if (!berhasil) {
    cat("Peta poligon tidak dapat digambar; peta titik ditampilkan sebagai pengganti.\n")
    print(f_titik())
  }
  invisible(NULL)
}

peta_aman(function() gg_poligon_num("Harga_Beras", "Sebaran Harga Eceran Beras (Des 2025)"),
          function() gg_titik_num("Harga_Beras", "Sebaran Harga Eceran Beras (Des 2025)"))

peta_aman(function() gg_poligon_kat("quadrant_f", "Kuadran Local Moran — Harga Beras"),
          function() gg_titik_kat("quadrant_f", "Kuadran Local Moran — Harga Beras"))

peta_aman(function() gg_poligon_kat("cluster_f", "Signifikansi Autokorelasi Spasial (LISA) — Harga Beras"),
          function() gg_titik_kat("cluster_f", "Signifikansi Autokorelasi Spasial (LISA) — Harga Beras"))

3.7.4 Sebaran klaster menurut provinsi

tb <- table(d_peta$Provinsi, d_peta$cluster_f)
knitr::kable(as.data.frame.matrix(tb), caption = "Jumlah wilayah per klaster LISA menurut provinsi")
Jumlah wilayah per klaster LISA menurut provinsi
High-High Low-Low High-Low Low-High Not Significant
Aceh 0 0 0 1 21
Bangka-Belitung 0 4 1 0 2
Bengkulu 0 1 1 0 7
Jambi 0 0 0 0 11
Kepulauan Riau 5 0 0 0 1
Lampung 0 12 0 0 0
Riau 9 0 0 0 3
Sumatera Barat 14 1 1 1 2
Sumatera Selatan 0 5 0 0 11
Sumatera Utara 0 0 0 0 28

Pembahasan. Klaster High-High (harga tinggi dikelilingi harga tinggi) terkonsentrasi di: Sumatera Barat (14), Riau (9), Kepulauan Riau (5). Klaster Low-Low (harga rendah dikelilingi harga rendah) terkonsentrasi di: Lampung (12), Sumatera Selatan (5), Bangka-Belitung (4) (angka dalam kurung adalah jumlah kab/kota). Pola ini mengindikasikan bahwa harga beras di Sumatera bersifat kewilayahan: wilayah yang bertetangga cenderung bernasib sama, misalnya karena berbagi jalur distribusi, sentra produksi, atau karakter pasar yang serupa. Wilayah High-Low/Low-High (5 wilayah) merupakan pencilan spasial yang harganya menyimpang dari tetangganya dan layak ditelaah lebih lanjut.

4 Kesimpulan

  1. Kehomogenan ragam. Ragam harga beras tidak homogen antar zona lintang (Bartlett p 0.0014; Levene p 0.0119). Hal ini wajar untuk data harga yang bervariasi menurut wilayah, dan perlu diperhatikan dalam pemodelan lanjutan.
  2. Autokorelasi global. Moran’s I = 0.573 dan Geary’s C = 0.379 sama-sama menunjukkan autokorelasi spasial positif yang signifikan: harga beras di kabupaten/kota yang berdekatan cenderung mirip.
  3. Autokorelasi lokal. LISA menemukan 29 wilayah High-High (antara lain Kep. Anambas, Kota Padangpanjang, Kota Bukittinggi, Bengkalis) dan 23 wilayah Low-Low (antara lain Mesuji, Pasawaran, OKU Timur, Pringsewu), serta 5 pencilan spasial (High-Low/Low-High).
  4. Catatan teknis. Kolom ID pada data dan shapefile tidak unik, sehingga merge(by = "ID") menggandakan sebagian poligon; peta statis yang dikoreksi (bila shapefile tersedia) disediakan untuk menghindari kekeliruan visual. Jumlah klaster LISA berbasis permutasi dapat bergeser sedikit antar-seed, tetapi kesimpulan umum tidak berubah.
  5. Saran. Analisis dapat dilanjutkan dengan autokorelasi spasial bivariat (Pertemuan 5) untuk melihat hubungan spasial harga beras dengan variabel lain, serta pemodelan regresi spasial bila diperlukan.

5 Informasi Sesi

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] tmap_4.4-1     mapview_2.11.4 ggplot2_4.0.2  spdep_1.4-2    sf_1.1-0      
## [6] spData_2.3.5  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        dplyr_1.1.4             farver_2.1.2           
##  [4] S7_0.2.1                fastmap_1.2.0           leaflegend_1.3.0       
##  [7] leaflet_2.2.3           XML_3.99-0.25           digest_0.6.39          
## [10] lifecycle_1.0.4         terra_1.9-50            magrittr_2.0.4         
## [13] dbscan_1.2.4            compiler_4.5.2          rlang_1.2.0            
## [16] sass_0.4.10             tools_4.5.2             leafpop_0.1.0          
## [19] yaml_2.3.11             data.table_1.17.8       knitr_1.50             
## [22] labeling_0.4.3          brew_1.0-10             htmlwidgets_1.6.4      
## [25] sp_2.2-1                classInt_0.4-11         RColorBrewer_1.1-3     
## [28] abind_1.4-8             KernSmooth_2.23-26      withr_3.0.2            
## [31] leafsync_0.1.0          grid_4.5.2              stats4_4.5.2           
## [34] cols4all_0.10           e1071_1.7-17            leafem_0.2.5           
## [37] colorspace_2.1-2        spacesXYZ_1.6-0         scales_1.4.0           
## [40] cli_3.6.5               rmarkdown_2.30          generics_0.1.4         
## [43] rstudioapi_0.17.1       tmaptools_3.3           DBI_1.2.3              
## [46] cachem_1.1.0            proxy_0.4-29            stars_0.7-3            
## [49] parallel_4.5.2          s2_1.1.9                base64enc_0.1-3        
## [52] vctrs_0.6.5             boot_1.3-32             jsonlite_2.0.0         
## [55] systemfonts_1.3.1       crosstalk_1.2.2         jquerylib_0.1.4        
## [58] units_1.0-0             maptiles_0.12.0         glue_1.8.0             
## [61] lwgeom_0.2-17           leaflet.providers_3.0.0 codetools_0.2-20       
## [64] gtable_0.3.6            deldir_2.0-4            raster_3.6-32          
## [67] tibble_3.3.0            logger_0.4.3            pillar_1.11.1          
## [70] htmltools_0.5.9         satellite_1.0.6         R6_2.6.1               
## [73] textshaping_1.0.4       wk_0.9.5                microbenchmark_1.5.0   
## [76] evaluate_1.0.5          lattice_0.22-7          png_0.1-9              
## [79] bslib_0.9.0             class_7.3-23            uuid_1.2-1             
## [82] Rcpp_1.1.0              svglite_2.2.2           xfun_0.61              
## [85] pkgconfig_2.0.3