Kesenjangan literasi digital antar kecamatan: analisis autokorelasi dan regresi spasial

# Contoh kode R: Spatial Autocorrelation and Spatial Regression (sf, spdep)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(sf)
  library(spdep)
  library(spatialreg)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Wilayah sintetis: sel grid 5 km (CRS UTM 49S) yang dipangkas menjadi bentuk wilayah tak beraturan ----
kotak <- sf::st_as_sfc(sf::st_bbox(c(xmin = 500000, ymin = 9100000, xmax = 570000, ymax = 9160000), crs = sf::st_crs(32749)))
sel <- sf::st_make_grid(kotak, cellsize = 5000, square = TRUE)
pusat_sel <- sf::st_coordinates(sf::st_centroid(sel))
u <- (pusat_sel[, 1] - 535000) / 33000
v <- (pusat_sel[, 2] - 9130000) / 28000
fase <- runif(1, 0, 2 * pi)
batas <- 1 + 0.22 * sin(3 * atan2(v, u) + fase) + 0.1 * cos(5 * atan2(v, u))
peta <- sf::st_sf(wilayah = NA_character_, geometry = sel[u^2 + v^2 < batas])
peta$wilayah <- sprintf("W%03d", seq_len(nrow(peta)))
xy <- sf::st_coordinates(sf::st_centroid(sf::st_geometry(peta)))
peta$koord_x <- round(xy[, 1]); peta$koord_y <- round(xy[, 2])

# ---- Matriks bobot spasial: tetangga queen, standarisasi baris ----
nb <- spdep::poly2nb(peta, queen = TRUE)
lw <- spdep::nb2listw(nb, style = "W", zero.policy = TRUE)
print(summary(nb))

# ---- Data sintetis: Y mengikuti proses SAR, Y = (I - rho W)^-1 (X b + e) ----
m <- nrow(peta)
W <- spdep::listw2mat(lw)
tren <- as.numeric(scale(xy[, 1] * runif(1, -1, 1) + xy[, 2] * runif(1, -1, 1)))
z1 <- 0.6 * tren + rnorm(m, 0, 0.8)                     # X1 sedikit berpola spasial
z2 <- rnorm(m)
rho_benar <- runif(1, 0.45, 0.7)
b1 <- runif(1, 0.5, 0.8); b2 <- runif(1, 0.2, 0.45) * sample(c(-1, 1), 1)
zy <- as.numeric(solve(diag(m) - rho_benar * W, b1 * z1 + b2 * z2 + rnorm(m, 0, runif(1, 0.5, 0.7))))
positif <- function(x, mu) if (mu > 0) pmax(x, 0.05 * mu) else x   # cegah nilai negatif yang tidak wajar
peta$digital_literacy_index <- round(positif(100 + 15 * as.numeric(scale(zy)), 100), 3)
peta$x1 <- round(positif(50 + 10 * as.numeric(scale(z1)), 50), 3)
peta$population_density <- round(positif(30 + 6 * z2, 30), 3)
d <- sf::st_drop_geometry(peta)                          # kolom geometri tidak ditulis ke CSV
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Autokorelasi spasial global (Moran's I) dan lokal (LISA) ----
mi <- spdep::moran.test(peta$digital_literacy_index, lw, zero.policy = TRUE)
print(mi)
lisa <- spdep::localmoran(peta$digital_literacy_index, lw, zero.policy = TRUE)
zy_std <- as.numeric(scale(peta$digital_literacy_index))
lag_y <- spdep::lag.listw(lw, zy_std, zero.policy = TRUE)
p_lisa <- lisa[, 5]
peta$lisa <- ifelse(p_lisa >= 0.05, "Not significant",
                    ifelse(zy_std > 0 & lag_y > 0, "High-High", ifelse(zy_std < 0 & lag_y < 0, "Low-Low",
                           ifelse(zy_std > 0, "High-Low", "Low-High"))))
print(table(peta$lisa))

# ---- Regresi OLS, uji pengganda Lagrange, dan model spasial lag (SAR) ----
ols <- lm(digital_literacy_index ~ x1 + population_density, data = peta)
print(summary(ols))
print(spdep::lm.morantest(ols, lw, zero.policy = TRUE))
rs <- spdep::lm.RStests(ols, lw, test = c("RSerr", "RSlag", "adjRSerr", "adjRSlag"), zero.policy = TRUE)
print(rs)
sar <- spatialreg::lagsarlm(digital_literacy_index ~ x1 + population_density, data = peta, listw = lw, zero.policy = TRUE)
print(summary(sar))
rho <- as.numeric(sar$rho)
lr <- 2 * (as.numeric(logLik(sar)) - as.numeric(logLik(ols)))
p_lr <- pchisq(lr, df = 1, lower.tail = FALSE)
aic_ols <- AIC(ols); aic_sar <- AIC(sar)


# ---- Grafik ----
tema_peta <- theme(axis.text = element_blank(), axis.title = element_blank(), panel.grid = element_blank())
p1 <- ggplot(peta) +
  geom_sf(aes(fill = digital_literacy_index), colour = warna$panel, linewidth = 0.25) +
  scale_fill_gradient(low = warna$grid, high = pal[1], name = NULL) +
  coord_sf(datum = NA) +
  labs(title = "Choropleth map", subtitle = paste0("digital_literacy_index", " by ", tolower("Group"))) +
  tema_peta + theme(legend.key.width = unit(0.5, "lines"))
warna_lisa <- c("High-High" = pal[2], "Low-Low" = pal[1], "High-Low" = pal[5], "Low-High" = pal[4], "Not significant" = warna$grid)
peta$lisa <- factor(peta$lisa, levels = names(warna_lisa))
p2 <- ggplot(peta) +
  geom_sf(aes(fill = lisa), colour = warna$panel, linewidth = 0.25, show.legend = TRUE) +
  scale_fill_manual(values = warna_lisa, drop = FALSE, name = NULL) +
  coord_sf(datum = NA) +
  labs(title = "LISA cluster map (p < 0.05)", subtitle = "Queen contiguity, row-standardised W") +
  tema_peta + theme(legend.key.size = unit(0.6, "lines"))
dm <- data.frame(z = zy_std, lag = lag_y, lisa = peta$lisa)
p3 <- ggplot(dm, aes(z, lag)) +
  geom_hline(yintercept = 0, colour = warna$redup, linetype = "dashed") +
  geom_vline(xintercept = 0, colour = warna$redup, linetype = "dashed") +
  geom_point(aes(colour = lisa), size = 1.8, alpha = 0.85, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, colour = warna$teks, linewidth = 0.9) +
  scale_colour_manual(values = replace(warna_lisa, "Not significant", warna$redup)) +
  labs(title = sprintf("Moran scatterplot (I = %.3f)", mi$estimate[1]), x = "digital_literacy_index (z-score)", y = "Spatial lag of z")
p1 <- wrap_elements(full = p1)   # peta berasio tetap: isi penuh sel tata letak
p2 <- wrap_elements(full = p2)
p_hasil <- if (tata == 1) ((p1 | p2) / p3) + plot_layout(heights = c(1.25, 1)) else
  (p3 | (p1 / p2)) + plot_layout(widths = c(1, 1.1))

print(p_hasil)