Analisis Autokorelasi Spasial Rata-Rata Lama Sekolah Kabupaten/Kota di Pulau Sumatera Tahun 2025 Menggunakan Moran’s I, Geary’s C, dan LISA

Studi Kasus

Pendidikan merupakan salah satu aspek penting dalam pembangunan kualitas sumber daya manusia. Salah satu indikator yang dapat digunakan untuk menggambarkan kondisi pendidikan suatu wilayah adalah Rata-Rata Lama Sekolah (RLS). RLS menunjukkan rata-rata jumlah tahun yang telah ditempuh oleh penduduk dalam mengikuti pendidikan formal. Perbedaan kondisi sosial, ekonomi, infrastruktur, serta akses terhadap fasilitas pendidikan dapat menyebabkan nilai RLS berbeda antarwilayah. Dalam analisis wilayah, nilai RLS pada suatu kabupaten/kota tidak selalu berdiri sendiri. Wilayah yang berdekatan secara geografis dapat memiliki karakteristik pendidikan yang serupa karena adanya kesamaan kondisi sosial ekonomi, aksesibilitas, fasilitas pendidikan, maupun hubungan antarwilayah. Oleh karena itu, diperlukan analisis spasial untuk mengetahui apakah terdapat pola keterkaitan geografis pada RLS kabupaten/kota di Pulau Sumatera.

Studi kasus ini menganalisis Rata-Rata Lama Sekolah kabupaten/kota di Pulau Sumatera tahun 2025. Analisis dilakukan untuk mengetahui apakah nilai RLS memiliki autokorelasi spasial, yaitu apakah kabupaten/kota dengan nilai RLS tinggi cenderung berdekatan dengan wilayah yang juga memiliki nilai RLS tinggi, maupun sebaliknya. Analisis dilakukan menggunakan Moran’s I Global dan Geary’s C untuk mengidentifikasi autokorelasi spasial secara keseluruhan. Selanjutnya, Local Moran’s I atau Local Indicators of Spatial Association (LISA) digunakan untuk mengetahui pola lokal pada setiap kabupaten/kota serta mengidentifikasi wilayah yang termasuk dalam kategori High-High, Low-Low, High-Low, dan Low-High. Sebelum analisis autokorelasi spasial dilakukan, data juga dibagi menjadi tiga zona berdasarkan posisi lintang, yaitu Selatan, Tengah, dan Utara, untuk menguji kehomogenan ragam RLS antarwilayah menggunakan Bartlett Test dan Levene Test/Brown-Forsythe.

Catatan: data yang digunakan adalah data simulasi untuk latihan (bukan data resmi BPS), sehingga seluruh interpretasi berlaku untuk data simulasi tersebut.

Deskripsi Data

Data yang digunakan merupakan data Rata-Rata Lama Sekolah (RLS) kabupaten/kota di Pulau Sumatera tahun 2025. Unit observasi dalam penelitian ini adalah kabupaten/kota yang berada di wilayah Pulau Sumatera. Setiap observasi dilengkapi dengan informasi lokasi geografis berupa koordinat bujur dan lintang sehingga dapat digunakan dalam pembentukan hubungan ketetanggaan spasial. Variabel yang digunakan dalam analisis terdiri atas:

# Membuat tabel deskripsi variabel
tabel_variabel <- data.frame(
  Variabel = c(
    "ID",
    "KabKot",
    "Long",
    "Lat",
    "RLS_2025",
    "Zona"
  ),
  Keterangan = c(
    "Kode identifikasi masing-masing kabupaten/kota yang digunakan untuk menghubungkan data dengan data spasial atau shapefile.",
    "Nama kabupaten atau kota yang menjadi unit observasi.",
    "Koordinat bujur (longitude) yang menunjukkan posisi geografis kabupaten/kota.",
    "Koordinat lintang (latitude) yang menunjukkan posisi geografis kabupaten/kota.",
    "Rata-Rata Lama Sekolah kabupaten/kota pada tahun 2025 dengan satuan tahun.",
    "Pengelompokan wilayah menjadi zona Selatan, Tengah, dan Utara berdasarkan tercile koordinat lintang."
  )
)

# Menampilkan tabel
knitr::kable(
  tabel_variabel,
  align = c("c", "l"),
  caption = "Deskripsi Variabel Data RLS Kabupaten/Kota di Pulau Sumatera Tahun 2025"
)

Deskripsi Variabel Data RLS Kabupaten/Kota di Pulau Sumatera Tahun 2025

Variabel	Keterangan
ID	Kode identifikasi masing-masing kabupaten/kota yang digunakan untuk menghubungkan data dengan data spasial atau shapefile.
KabKot	Nama kabupaten atau kota yang menjadi unit observasi.
Long	Koordinat bujur (longitude) yang menunjukkan posisi geografis kabupaten/kota.
Lat	Koordinat lintang (latitude) yang menunjukkan posisi geografis kabupaten/kota.
RLS_2025	Rata-Rata Lama Sekolah kabupaten/kota pada tahun 2025 dengan satuan tahun.
Zona	Pengelompokan wilayah menjadi zona Selatan, Tengah, dan Utara berdasarkan tercile koordinat lintang.

Variabel utama yang dianalisis adalah RLS_2025, sedangkan variabel Long dan Lat digunakan untuk membentuk struktur ketetanggaan spasial. Hubungan antarwilayah ditentukan menggunakan metode k-Nearest Neighbour (k-NN) dengan nilai \(k=6\), sehingga setiap kabupaten/kota dihubungkan dengan enam wilayah terdekat berdasarkan posisi geografisnya. Selanjutnya, bobot spasial distandarisasi menggunakan row-standardized weights atau style = “W”. Data RLS kemudian dianalisis untuk mengetahui pola spasial secara global maupun lokal. Moran’s I dan Geary’s C digunakan untuk mengevaluasi keberadaan autokorelasi spasial pada keseluruhan wilayah, sedangkan Local Moran’s I digunakan untuk mengidentifikasi klaster maupun pencilan spasial pada masing-masing kabupaten/kota.

Penyelesaian dan Analisis

Load package

library(readxl)
library(spdep)
library(ggplot2)
library(car)
library(sf)
library(tmap)
library(mapview)

Impor dan Pembersihan Data

# Impor Data
d <- read_excel(
  "data_simulasi_rls_sumatera.xlsx",
  sheet = "Data_RLS",
  skip = 3
)

# Melihat nama kolom asli
names(d)

## [1] "ID"                 "Kabupaten/Kota"     "Bujur (Longitude)" 
## [4] "Lintang (Latitude)" "RLS_2025_Tahun"

# Mengubah nama kolom agar lebih mudah digunakan
names(d) <- c(
  "ID",
  "KabKot",
  "Long",
  "Lat",
  "RLS_2025"
)

# Membersihkan data
# Menghapus baris yang bukan data kabupaten/kota,
d <- d[
  !is.na(d$ID) &
    !is.na(d$KabKot) &
    !is.na(d$Long) &
    !is.na(d$Lat) &
    !is.na(d$RLS_2025),
]

# Merapikan tipe variabel

d$ID       <- as.character(d$ID)
d$KabKot   <- trimws(as.character(d$KabKot))
d$Long     <- as.numeric(d$Long)
d$Lat      <- as.numeric(d$Lat)
d$RLS_2025 <- as.numeric(d$RLS_2025)

# Mengecek data
str(d)

## tibble [143 × 5] (S3: tbl_df/tbl/data.frame)
##  $ ID      : chr [1:143] "293" "277" "278" "279" ...
##  $ KabKot  : chr [1:143] "Aceh Barat" "Aceh Barat Daya" "Aceh Besar" "Aceh Jaya" ...
##  $ Long    : num [1:143] 96.1 96.8 95.3 95.8 96 ...
##  $ Lat     : num [1:143] 4.14 3.74 5.51 4.47 5.38 3.26 2.4 4.29 4.62 3.49 ...
##  $ RLS_2025: num [1:143] 8.2 8.3 8.5 8.1 8.7 8.6 8.5 8.8 8.5 8.6 ...

head(d)

## # A tibble: 6 × 5
##   ID    KabKot           Long   Lat RLS_2025
##   <chr> <chr>           <dbl> <dbl>    <dbl>
## 1 293   Aceh Barat       96.1  4.14      8.2
## 2 277   Aceh Barat Daya  96.8  3.74      8.3
## 3 278   Aceh Besar       95.4  5.51      8.5
## 4 279   Aceh Jaya        95.8  4.47      8.1
## 5 405   Aceh Pidie       96.0  5.38      8.7
## 6 365   Aceh Selatan     97.2  3.26      8.6

tail(d)

## # A tibble: 6 × 5
##   ID    KabKot           Long   Lat RLS_2025
##   <chr> <chr>           <dbl> <dbl>    <dbl>
## 1 270   Tapanuli Tengah  98.4  2.01      8.2
## 2 352   Tapanuli Utara   99.0  2.03      8.2
## 3 272   Tebo            102.  -1.49      8.5
## 4 271   Toba Samosir     99.1  2.33      8.3
## 5 243   Tulangbawang    105.  -4.49      8.6
## 6 312   Way Kanan       105.  -4.72      8.5

nrow(d)

## [1] 143

summary(d$RLS_2025)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.000   8.500   8.700   8.729   8.900   9.700

# Statistika deskriptif
summary(d$RLS_2025)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.000   8.500   8.700   8.729   8.900   9.700

cat(
  "\nRata-rata:",
  mean(d$RLS_2025, na.rm = TRUE)
)

## 
## Rata-rata: 8.729371

cat(
  "\nMedian:",
  median(d$RLS_2025, na.rm = TRUE)
)

## 
## Median: 8.7

cat(
  "\nStandar Deviasi:",
  sd(d$RLS_2025, na.rm = TRUE)
)

## 
## Standar Deviasi: 0.3431076

cat(
  "\nMinimum:",
  min(d$RLS_2025, na.rm = TRUE)
)

## 
## Minimum: 8

cat(
  "\nMaksimum:",
  max(d$RLS_2025, na.rm = TRUE)
)

## 
## Maksimum: 9.7

Setelah dibersihkan, data terdiri atas 143 kabupaten/kota dengan lima variabel (ID, KabKot, Long, Lat, RLS_2025) tanpa data hilang. Rata-rata RLS Pulau Sumatera sebesar 8,73 tahun dengan median 8,70 tahun. Nilai rata-rata yang hampir sama dengan median menunjukkan sebaran RLS cukup simetris. Standar deviasi sebesar 0,34 tahun tergolong kecil, artinya RLS antarwilayah relatif seragam. Nilai terendah adalah 8,0 tahun (Kaur) dan tertinggi 9,7 tahun (Kota Payakumbuh), sehingga rentangnya hanya 1,7 tahun. Kuartil 1 dan kuartil 3 masing-masing 8,5 dan 8,9 tahun, yang berarti separuh wilayah berada pada kisaran yang sempit.

Pembentukan zona

terc <- quantile(
  d$Lat,
  probs = c(1/3, 2/3),
  na.rm = TRUE
)

terc

## 33.33333% 66.66667% 
## -2.026667  1.500000

d$Zona <- cut(
  d$Lat,
  breaks = c(
    -Inf,
    terc[1],
    terc[2],
    Inf
  ),
  labels = c(
    "Selatan",
    "Tengah",
    "Utara"
  ),
  include.lowest = TRUE
)


# Mengecek jumlah kab/kota setiap zona
table(d$Zona)

## 
## Selatan  Tengah   Utara 
##      48      47      48

# Statistik RLS per zona
aggregate(
  RLS_2025 ~ Zona,
  data = d,
  FUN = mean
)

##      Zona RLS_2025
## 1 Selatan 8.589583
## 2  Tengah 8.948936
## 3   Utara 8.654167

aggregate(
  RLS_2025 ~ Zona,
  data = d,
  FUN = sd
)

##      Zona  RLS_2025
## 1 Selatan 0.2243094
## 2  Tengah 0.3810012
## 3   Utara 0.2989046

Batas tercile lintang berada di -2,03 dan 1,50. Wilayah dengan lintang di bawah -2,03 masuk zona Selatan, di antara -2,03 sampai 1,50 masuk zona Tengah, dan di atas 1,50 masuk zona Utara. Pembagian ini menghasilkan zona yang hampir seimbang: Selatan 48, Tengah 47, Utara 48 kabupaten/kota. Rata-rata RLS tertinggi ada di zona Tengah (8,95 tahun), diikuti Utara (8,65 tahun) dan Selatan (8,59 tahun). Simpangan baku zona Tengah juga paling besar (0,38) dibandingkan Utara (0,30) dan Selatan (0,22). Artinya, zona Tengah tidak hanya paling tinggi, tetapi juga paling beragam. Perbedaan keragaman ini diuji secara formal pada bagian berikutnya.

Uji kehomogenan ragam

Bartlett Test

Hipotesis

\[H_0 : \sigma^2_{Selatan} = \sigma^2_{Tengah} = \sigma^2_{Utara}\]

\[H_1 : \text{minimal ada satu } \sigma^2_i \text{ yang berbeda}\]

Statistik uji: Bartlett’s \(K^2\) yang berdistribusi \(\chi^2_{(k-1)}\) dengan \(k=3\) zona, sehingga derajat bebasnya 2.

Daerah penolakan (\(\alpha = 0{,}05\)): tolak \(H_0\) jika \(K^2 > \chi^2_{0{,}05;2} = 5{,}991\) atau \(p\text{-value} < 0{,}05\).

bartlett_res <- bartlett.test(
  RLS_2025 ~ Zona,
  data = d
)

bartlett_res

## 
##  Bartlett test of homogeneity of variances
## 
## data:  RLS_2025 by Zona
## Bartlett's K-squared = 12.517, df = 2, p-value = 0.001914

# Keputusan otomatis
if (bartlett_res$p.value < 0.05) {
  
  cat(
    "\nKeputusan Bartlett: Tolak H0",
    "\nKesimpulan: Ragam RLS antar zona tidak homogen.\n"
  )
  
} else {
  
  cat(
    "\nKeputusan Bartlett: Gagal menolak H0",
    "\nKesimpulan: Belum terdapat bukti bahwa ragam RLS antar zona berbeda.\n"
  )
}

## 
## Keputusan Bartlett: Tolak H0 
## Kesimpulan: Ragam RLS antar zona tidak homogen.

Diperoleh \(K^2 = 12{,}517\) dengan \(db = 2\) dan \(p\text{-value} = 0{,}0019\). Karena \(12{,}517 > 5{,}991\) (dan \(p < 0{,}05\)), \(H_0\) ditolak. Pada taraf signifikansi 5%, ragam RLS antarzona tidak homogen. Hasil ini sejalan dengan statistik deskriptif: ragam zona Tengah (0,145) sekitar tiga kali ragam zona Selatan (0,050).

Levene Test / Brown-Forsythe

Hipotesis

\[H_0 : \sigma^2_{Selatan} = \sigma^2_{Tengah} = \sigma^2_{Utara}\]

\[H_1 : \text{minimal ada satu } \sigma^2_i \text{ yang berbeda}\]

Statistik uji: \(F\) dari ANOVA satu arah terhadap simpangan absolut data dari pusatnya (median untuk Brown-Forsythe, yang menjadi default car::leveneTest), dengan \(db_1 = k-1 = 2\) dan \(db_2 = n-k = 140\).

Daerah penolakan: (\(\alpha = 0{,}05\)): tolak \(H_0\) jika \(F > F_{0{,}05;2;140} = 3{,}061\) atau \(p\text{-value} < 0{,}05\).

levene_res <- car::leveneTest(
  RLS_2025 ~ Zona,
  data = d
)

levene_res

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)    
## group   2  8.6384 0.00029 ***
##       140                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Mengambil p-value Levene
p_levene <- levene_res$`Pr(>F)`[1]


# Keputusan otomatis
if (p_levene < 0.05) {
  
  cat(
    "\nKeputusan Levene: Tolak H0",
    "\nKesimpulan: Ragam RLS antar zona tidak homogen.\n"
  )
  
} else {
  
  cat(
    "\nKeputusan Levene: Gagal menolak H0",
    "\nKesimpulan: Belum terdapat bukti bahwa ragam RLS antar zona berbeda.\n"
  )
}

## 
## Keputusan Levene: Tolak H0 
## Kesimpulan: Ragam RLS antar zona tidak homogen.

Diperoleh \(F = 8{,}638\) dengan \(p\text{-value} = 0{,}0003\). Karena \(8{,}638 > 3{,}061\), \(H_0\) ditolak. Kesimpulannya sama dengan Bartlett: ragam RLS antarzona tidak homogen. Levene/Brown-Forsythe lebih tahan terhadap pelanggaran normalitas dibanding Bartlett, sehingga konsistensi kedua uji memperkuat kesimpulan bahwa keragaman RLS memang berbeda antarzona (heterogenitas spasial). Implikasinya, bila RLS antarzona dibandingkan dengan metode yang mengasumsikan ragam sama (misalnya ANOVA klasik), asumsi tersebut tidak terpenuhi.

Pembentukan matriks pembobot spasial

# Membentuk matriks koordinat
koordinat <- cbind(
  d$Long,
  d$Lat
)

head(koordinat)

##       [,1] [,2]
## [1,] 96.13 4.14
## [2,] 96.83 3.74
## [3,] 95.35 5.51
## [4,] 95.79 4.47
## [5,] 95.95 5.38
## [6,] 97.18 3.26

# Membentuk k-nearest neighbour
knn <- knearneigh(
  koordinat,
  k = 6
)


# Mengubah menjadi neighbour list
nb <- knn2nb(
  knn,
  row.names = d$KabKot
)

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

# Membentuk spatial weight

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

W

## Characteristics of weights list object:
## 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
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 143 20449 143 41.88889 590.1667

Struktur ketetanggaan terbentuk untuk 143 wilayah dengan 858 hubungan (\(143 \times 6\)). Setiap wilayah memiliki tepat enam tetangga (rata-rata 6 tautan), sehingga tidak ada wilayah tanpa tetangga (isolated) dan pendekatan k-NN cocok untuk Sumatera yang berbentuk kepulauan dan memiliki wilayah dengan ukuran sangat beragam. Label non-symmetric berarti hubungan tidak selalu timbal balik: wilayah A bisa menjadi tetangga B tanpa B menjadi tetangga A. Hal ini wajar pada k-NN. Dengan style = "W", bobot tiap wilayah dibagi rata (\(w_{ij} = 1/6\)), sehingga spatial lag sama dengan rata-rata RLS dari enam tetangga terdekat.

Moran’s I Global

Hipotesis

\[H_0 : I = 0 \quad \text{(tidak ada autokorelasi spasial; RLS tersebar acak)}\]

\[H_1 : I > 0 \quad \text{(terdapat autokorelasi spasial positif)}\]

(moran.test secara default menggunakan uji satu arah greater.)

Statistik uji:

\[Z(I) = \frac{I - E(I)}{\sqrt{Var(I)}}, \qquad E(I) = -\frac{1}{n-1} = -\frac{1}{142} = -0{,}00704\]

Daerah penolakan: (\(\alpha = 0{,}05\)): tolak \(H_0\) jika \(Z(I) > Z_{0{,}05} = 1{,}645\) atau \(p\text{-value} < 0{,}05\).

moran_res <- moran.test(
  d$RLS_2025,
  W,
  zero.policy = TRUE
)

moran_res

## 
##  Moran I test under randomisation
## 
## data:  d$RLS_2025  
## weights: W    
## 
## Moran I statistic standard deviate = 16.815, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.733876321      -0.007042254       0.001941571

# Keputusan Moran's I
p_moran <- moran_res$p.value

I_moran <- unname(
  moran_res$estimate["Moran I statistic"]
)


if (p_moran < 0.05) {
  
  if (I_moran > 0) {
    
    cat(
      "\nKeputusan Moran: Tolak H0",
      "\nKesimpulan: Terdapat autokorelasi spasial positif yang signifikan.\n"
    )
    
  } else {
    
    cat(
      "\nKeputusan Moran: Tolak H0",
      "\nKesimpulan: Terdapat autokorelasi spasial negatif yang signifikan.\n"
    )
  }
  
} else {
  
  cat(
    "\nKeputusan Moran: Gagal menolak H0",
    "\nKesimpulan: Belum terdapat bukti autokorelasi spasial yang signifikan.\n"
  )
}

## 
## Keputusan Moran: Tolak H0 
## Kesimpulan: Terdapat autokorelasi spasial positif yang signifikan.

Nilai Moran’s I sebesar 0,7339, jauh di atas nilai harapannya (\(E(I) = -0{,}0070\)). Statistik uji \(Z = 16{,}815 > 1{,}645\) dan \(p\text{-value} < 2{,}2 \times 10^{-16}\), sehingga \(H_0\) ditolak. Terdapat autokorelasi spasial positif yang sangat signifikan pada RLS kabupaten/kota di Sumatera. Nilai I yang mendekati 1 menunjukkan keterkaitan yang kuat: wilayah dengan RLS tinggi cenderung dikelilingi wilayah ber-RLS tinggi, dan wilayah ber-RLS rendah cenderung dikelilingi wilayah ber-RLS rendah. Pola RLS di Sumatera jelas mengelompok (clustered), bukan acak.

Moran Monte Carlo

Hipotesis: sama dengan Moran’s I Global (\(H_0: I = 0\) vs \(H_1: I > 0\)). Uji ini tidak bergantung pada asumsi distribusi karena membandingkan \(I\) observasi dengan \(I\) dari 999 permutasi acak nilai RLS antarwilayah.

Daerah penolakan (\(\alpha = 0{,}05\)): tolak \(H_0\) jika \(p\text{-value} = \dfrac{R + 1}{nsim + 1} < 0{,}05\), dengan \(R\) = banyaknya \(I\) simulasi yang lebih besar atau sama dengan \(I\) observasi.

set.seed(123)

moran_mc <- moran.mc(
  d$RLS_2025,
  W,
  nsim = 999,
  zero.policy = TRUE
)

moran_mc

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  d$RLS_2025 
## weights: W  
## number of simulations + 1: 1000 
## 
## statistic = 0.73388, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Moran’s I observasi (0,7339) berada pada peringkat 1000 dari 1000 (lebih besar dari seluruh 999 hasil simulasi), sehingga \(p\text{-value} = 1/1000 = 0{,}001 < 0{,}05\) dan \(H_0\) ditolak. Nilai \(p = 0{,}001\) adalah nilai terkecil yang mungkin dengan 999 simulasi. Hasil ini mengonfirmasi uji Moran’s I sebelumnya secara nonparametrik: autokorelasi spasial positif bukan terjadi secara kebetulan.

Geary’s C

Hipotesis

\[H_0 : C = 1 \quad \text{(tidak ada autokorelasi spasial)}\]

\[H_1 : C < 1 \quad \text{(terdapat autokorelasi spasial positif)}\]

Statistik uji:

\[Z(C) = \frac{C - E(C)}{\sqrt{Var(C)}}, \qquad E(C) = 1\]

Daerah penolakan (\(\alpha = 0{,}05\)): tolak \(H_0\) jika \(Z(C) < -Z_{0{,}05} = -1{,}645\) atau \(p\text{-value} < 0{,}05\). Pada keluaran R, standard deviate ditampilkan sebagai \((E(C) - C)/\sqrt{Var(C)}\) (bertanda positif), sehingga kriteria tolak \(H_0\) menjadi nilai tersebut \(> 1{,}645\).

geary_res <- geary.test(
  d$RLS_2025,
  W,
  zero.policy = TRUE
)

geary_res

## 
##  Geary C test under randomisation
## 
## data:  d$RLS_2025 
## weights: W   
## 
## Geary C statistic standard deviate = 14.341, p-value < 2.2e-16
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.298100521       1.000000000       0.002395428

# Mengambil nilai Geary's C
C_geary <- unname(
  geary_res$estimate["Geary C statistic"]
)

p_geary <- geary_res$p.value

# Keputusan Geary's C
if (p_geary < 0.05) {
  
  if (C_geary < 1) {
    
    cat(
      "\nKeputusan Geary: Tolak H0",
      "\nKesimpulan: Terdapat autokorelasi spasial positif yang signifikan.\n"
    )
    
  } else if (C_geary > 1) {
    
    cat(
      "\nKeputusan Geary: Tolak H0",
      "\nKesimpulan: Terdapat autokorelasi spasial negatif yang signifikan.\n"
    )
    
  }
  
} else {
  
  cat(
    "\nKeputusan Geary: Gagal menolak H0",
    "\nKesimpulan: Belum terdapat bukti autokorelasi spasial yang signifikan.\n"
  )
}

## 
## Keputusan Geary: Tolak H0 
## Kesimpulan: Terdapat autokorelasi spasial positif yang signifikan.

Nilai Geary’s C sebesar 0,2981, jauh di bawah 1 (nilai harapan di bawah \(H_0\)). Statistik uji sebesar \(14{,}341 > 1{,}645\) dengan \(p\text{-value} < 2{,}2 \times 10^{-16}\), sehingga \(H_0\) ditolak. Karena \(C < 1\), terdapat autokorelasi spasial positif yang signifikan: selisih nilai RLS antarwilayah bertetangga jauh lebih kecil daripada yang diharapkan jika data tersebar acak. Kesimpulan ini konsisten dengan Moran’s I. Perlu diingat bahwa skala Geary’s C berkebalikan dengan Moran’s I (\(C\) kecil berarti autokorelasi positif kuat), dan Geary’s C lebih sensitif terhadap perbedaan nilai antar-tetangga langsung (pola lokal), sedangkan Moran’s I lebih menggambarkan pola global.

Moran Scatterplot

# Standardisasi RLS
z <- as.numeric(
  scale(d$RLS_2025)
)


# Menghitung spatial lag
wz <- lag.listw(
  W,
  z,
  zero.policy = TRUE
)


# Moran Scatterplot
moran.plot(
  d$RLS_2025,
  W,
  labels = d$KabKot,
  xlab = "Rata-Rata Lama Sekolah (Tahun)",
  ylab = "Spatial Lag RLS",
  main = "Moran Scatterplot RLS Kabupaten/Kota Pulau Sumatera",
  zero.policy = TRUE
)

Moran scatterplot memetakan RLS suatu wilayah (sumbu X) terhadap rata-rata RLS tetangganya (spatial lag, sumbu Y). Garis regresi memiliki kemiringan positif sebesar 0,7339, yang sama dengan nilai Moran’s I (karena bobot sudah row-standardized). Mayoritas titik berada di kuadran High-High (kanan atas) dan Low-Low (kiri bawah), yang menandakan kemiripan nilai antarwilayah bertetangga. Hanya sedikit titik di kuadran Low-High (kiri atas) dan High-Low (kanan bawah), yaitu wilayah yang berbeda dari tetangganya. Titik berlabel (simbol berlian) adalah pengamatan yang berpengaruh besar terhadap kemiringan garis.

Local Moran / LISA

Hipotesis (untuk setiap wilayah \(i\))

\[H_0 : I_i = 0 \quad \text{(tidak ada autokorelasi spasial lokal di wilayah } i)\]

\[H_1 : I_i \neq 0 \quad \text{(terdapat autokorelasi spasial lokal di wilayah } i)\]

Statistik uji:

\[Z(I_i) = \frac{I_i - E(I_i)}{\sqrt{Var(I_i)}}\]

Daerah penolakan (\(\alpha = 0{,}05\), dua arah sesuai keluaran localmoran): tolak \(H_0\) jika \(|Z(I_i)| > Z_{0{,}025} = 1{,}96\) atau \(p\text{-value} < 0{,}05\).

lisa <- localmoran(
  d$RLS_2025,
  W,
  zero.policy = TRUE
)

head(lisa)

##                         Ii          E.Ii      Var.Ii      Z.Ii Pr(z != E(Ii))
## Aceh Barat      1.71794995 -0.0168817825 0.381529942 2.8086204    0.004975428
## Aceh Barat Daya 1.08734090 -0.0111061436 0.252474514 2.1861016    0.028808170
## Aceh Besar      0.25384023 -0.0031693794 0.072627347 0.9536725    0.340249426
## Aceh Jaya       1.68355309 -0.0238622593 0.535460263 2.3333267    0.019631003
## Aceh Pidie      0.04087882 -0.0000519667 0.001194558 1.1842585    0.236310744
## Aceh Selatan    0.14317208 -0.0010082541 0.023154557 0.9475187    0.343374526

# Memasukkan hasil LISA ke data
d$Ii <- lisa[, 1]

d$EIi <- lisa[, 2]

d$VarIi <- lisa[, 3]

d$ZIi <- lisa[, 4]

d$PIi <- lisa[, 5]


# Mengecek hasil
head(
  d[
    ,
    c(
      "KabKot",
      "RLS_2025",
      "Ii",
      "ZIi",
      "PIi"
    )
  ]
)

## # A tibble: 6 × 5
##   KabKot          RLS_2025     Ii   ZIi     PIi
##   <chr>              <dbl>  <dbl> <dbl>   <dbl>
## 1 Aceh Barat           8.2 1.72   2.81  0.00498
## 2 Aceh Barat Daya      8.3 1.09   2.19  0.0288 
## 3 Aceh Besar           8.5 0.254  0.954 0.340  
## 4 Aceh Jaya            8.1 1.68   2.33  0.0196 
## 5 Aceh Pidie           8.7 0.0409 1.18  0.236  
## 6 Aceh Selatan         8.6 0.143  0.948 0.343

Setiap baris adalah hasil uji untuk satu kabupaten/kota. Contoh pada enam wilayah pertama: Aceh Barat (\(I_i = 1{,}718\); \(Z = 2{,}81\); \(p = 0{,}005\)), Aceh Barat Daya (\(I_i = 1{,}087\); \(Z = 2{,}19\); \(p = 0{,}029\)), dan Aceh Jaya (\(I_i = 1{,}684\); \(Z = 2{,}33\); \(p = 0{,}020\)) memiliki \(|Z| > 1{,}96\) sehingga \(H_0\) ditolak (autokorelasi lokal signifikan). Sebaliknya, Aceh Besar (\(p = 0{,}340\)), Aceh Pidie (\(p = 0{,}236\)), dan Aceh Selatan (\(p = 0{,}343\)) memiliki \(p > 0{,}05\) sehingga \(H_0\) gagal ditolak (tidak signifikan). Nilai \(I_i\) positif berarti wilayah tersebut mirip dengan tetangganya, sedangkan \(I_i\) negatif berarti berbeda dari tetangganya. Tanda (\(+\)/\(-\)) saja belum menunjukkan jenis klasternya (tinggi atau rendah), sehingga diperlukan klasifikasi kuadran pada langkah berikut.

Klasifikasi kuadran dan klaster

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"
    )
  )
)


# Melihat jumlah kabupaten/kota tiap kuadran
table(d$quadrant)

## 
## High-High  High-Low  Low-High   Low-Low 
##        48        10        14        71

# Signifikansi local moran
d$sig <- d$PIi < 0.05


# Jika tidak signifikan,
# kategori diubah menjadi "Not Significant"

d$cluster <- ifelse(
  d$sig,
  d$quadrant,
  "Not Significant"
)


# Distribusi klaster
table(d$cluster)

## 
##       High-High        Low-High         Low-Low Not Significant 
##              24               1              17             101

# Wilayah dengan local moran signifikan
hasil_lisa <- d[
  d$sig,
  c(
    "ID",
    "KabKot",
    "RLS_2025",
    "Ii",
    "ZIi",
    "PIi",
    "quadrant",
    "cluster"
  )
]


# Urutkan berdasarkan nilai Local Moran terbesar
hasil_lisa <- hasil_lisa[
  order(
    -hasil_lisa$Ii
  ),
]


hasil_lisa

## # A tibble: 42 × 8
##    ID    KabKot             RLS_2025    Ii   ZIi          PIi quadrant  cluster 
##    <chr> <chr>                 <dbl> <dbl> <dbl>        <dbl> <chr>     <chr>   
##  1 118   Kota Payahkumbuh        9.7  5.71  5.20 0.000000204  High-High High-Hi…
##  2 128   Kota Padangpanjang      9.6  5.49  5.53 0.0000000317 High-High High-Hi…
##  3 116   Kota Pariaman           9.6  5.24  5.28 0.000000126  High-High High-Hi…
##  4 296   Padang Pariaman         9.4  4.23  5.47 0.0000000454 High-High High-Hi…
##  5 1     Agam                    9.4  4.13  5.35 0.0000000901 High-High High-Hi…
##  6 95    Kota Bukittinggi        9.4  4.13  5.35 0.0000000901 High-High High-Hi…
##  7 123   Kota Sawahlunto         9.5  4.09  4.63 0.00000358   High-High High-Hi…
##  8 112   Kota Padang             9.5  3.87  4.39 0.0000114    High-High High-Hi…
##  9 393   Lima Puluh Kota         9.3  3.68  5.56 0.0000000263 High-High High-Hi…
## 10 228   Tanah Datar             9.3  3.60  5.44 0.0000000526 High-High High-Hi…
## # ℹ 32 more rows

Kuadran (sebelum uji signifikansi). Dari 143 wilayah, 48 High-High, 71 Low-Low, 14 Low-High, dan 10 High-Low. Sebanyak 119 wilayah (83,2%) berada pada kuadran searah (HH atau LL), sejalan dengan Moran’s I yang positif dan kuat. Klaster signifikan (\(\alpha = 5\%\)). Hanya 42 wilayah yang signifikan, sedangkan 101 wilayah tidak signifikan:

tabel_cluster <- data.frame(
  Klaster = c(
    "High-High",
    "Low-Low",
    "Low-High",
    "High-Low",
    "Not Significant"
  ),
  Jumlah = c(
    24,
    17,
    1,
    0,
    101
  ),
  Makna = c(
    "RLS tinggi, dikelilingi tetangga ber-RLS tinggi (hot spot)",
    "RLS rendah, dikelilingi tetangga ber-RLS rendah (cold spot)",
    "RLS rendah, dikelilingi tetangga ber-RLS tinggi (pencilan spasial)",
    "Tidak ada wilayah ber-RLS tinggi yang signifikan di tengah tetangga rendah",
    "Tidak ada bukti autokorelasi lokal"
  )
)

knitr::kable(
  tabel_cluster,
  align = c("c", "c", "l"),
  caption = "Distribusi Klaster LISA Rata-Rata Lama Sekolah Kabupaten/Kota di Pulau Sumatera Tahun 2025"
)

Distribusi Klaster LISA Rata-Rata Lama Sekolah Kabupaten/Kota di Pulau Sumatera Tahun 2025

Klaster	Jumlah	Makna
High-High	24	RLS tinggi, dikelilingi tetangga ber-RLS tinggi (hot spot)
Low-Low	17	RLS rendah, dikelilingi tetangga ber-RLS rendah (cold spot)
Low-High	1	RLS rendah, dikelilingi tetangga ber-RLS tinggi (pencilan spasial)
High-Low	0	Tidak ada wilayah ber-RLS tinggi yang signifikan di tengah tetangga rendah
Not Significant	101	Tidak ada bukti autokorelasi lokal

Hot spot (High-High). Klaster terkuat terpusat di Sumatera Barat, dengan nilai \(I_i\) tertinggi pada Kota Payakumbuh (9,7 tahun; \(I_i = 5{,}71\); \(Z = 5{,}20\)), diikuti Kota Padangpanjang, Kota Pariaman, Padang Pariaman, Agam, Kota Bukittinggi, Kota Sawahlunto, Kota Padang, Lima Puluh Kota, Tanah Datar, Pasaman, Kota Solok, Solok, Sijunjung, dan Pesisir Selatan. Klaster ini meluas ke Riau (Kampar, Kota Pekanbaru, Kuantan Singingi). Hot spot kecil lainnya muncul di Sumatera Utara pesisir timur (Batu Bara, Asahan, Serdang Bedagai, Kota Tebing Tinggi) dan Kepulauan Riau (Bintan, Kepulauan Anambas).

Cold spot (Low-Low). Terdapat tiga kelompok utama:

1. Aceh bagian barat (Aceh Barat, Aceh Jaya, Aceh Barat Daya, Nagan Raya)

2. Kawasan Tapanuli, Sumatera Utara (Tapanuli Utara, Tapanuli Tengah, Tapanuli Selatan, Humbang Hasundutan, Kota Gunung Sitoli, Kota Padangsidempuan)

3. Ujung selatan Sumatera (Kaur dan Bengkulu Selatan di Bengkulu; Lampung Barat dan Pesisir Barat di Lampung; OKU Selatan, OKU Timur, dan Kota Pagaralam di Sumatera Selatan). Kaur (RLS 8,0 tahun; \(I_i = 2{,}57\); \(Z = 3{,}09\)) merupakan cold spot dengan nilai \(I_i\) tertinggi.

Pencilan spasial. Hanya Pasaman Barat (RLS 8,7 tahun) yang signifikan sebagai Low-High: RLS-nya relatif rendah, padahal dikelilingi wilayah ber-RLS tinggi di Sumatera Barat (\(I_i = -0{,}08\); \(Z = -2{,}21\)).

Visualisasi peta

# Visualisasi sebaran RLS
ggplot(
  d,
  aes(
    x = Long,
    y = Lat
  )
) +
  
  geom_point(
    aes(
      size = RLS_2025,
      color = RLS_2025
    ),
    alpha = 0.85
  ) +
  
  scale_color_gradient(
    low = "#1C7293",
    high = "#C0392B"
  ) +
  
  labs(
    title =
      "Sebaran Rata-Rata Lama Sekolah Kabupaten/Kota Pulau Sumatera",
    subtitle =
      "Tahun 2025",
    x = "Bujur",
    y = "Lintang",
    color = "RLS (Tahun)",
    size = "RLS (Tahun)"
  ) +
  
  theme_minimal()

Titik berwarna merah dan berukuran besar menunjukkan RLS tinggi, sedangkan biru dan kecil menunjukkan RLS rendah. Terlihat pengelompokan warna secara geografis: kelompok titik merah menonjol di bagian tengah-barat pulau (sekitar Sumatera Barat dan Riau), sedangkan titik biru mengumpul di Aceh bagian barat, Tapanuli, dan ujung selatan pulau. Warna yang berdekatan cenderung sama, sesuai dengan hasil Moran’s I dan Geary’s C.

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


ggplot(
  d,
  aes(
    x = Long,
    y = Lat,
    color = cluster
  )
) +
  
  geom_point(
    size = 2.8,
    alpha = 0.9
  ) +
  
  scale_color_manual(
    values = cols
  ) +
  
  labs(
    title =
      "Klaster LISA Rata-Rata Lama Sekolah",
    subtitle =
      "Kabupaten/Kota Pulau Sumatera Tahun 2025",
    x = "Bujur",
    y = "Lintang",
    color = "Klaster LISA"
  ) +
  
  theme_minimal()

Titik merah tua (High-High) dan biru tua (Low-Low) tampak membentuk kantong-kantong yang terpisah: hot spot merah terkonsentrasi di sekitar Sumatera Barat–Riau, sedangkan cold spot biru tua tersebar di Aceh bagian barat, Tapanuli, dan Bengkulu–Lampung–Sumatera Selatan bagian barat daya. Sebagian besar titik berwarna abu-abu (tidak signifikan), yang menunjukkan bahwa autokorelasi spasial yang kuat secara global ditopang oleh beberapa klaster utama, bukan oleh seluruh wilayah.

# Membaca file SHP Sumatra sebagai peta dasar
shp <- st_read(
  "SHP Sumatera KabKot/Sumatera-KabKot.shp",
  quiet = TRUE
)

# Perbaiki geometri shapefile agar valid dan dapat digambar
sf_use_s2(FALSE)

shp_dasar <- st_make_valid(
  st_zm(shp, drop = TRUE, what = "ZM")
)

# Ubah data RLS menjadi titik berdasarkan koordinat
d_sf <- st_as_sf(
  d,
  coords = c("Long", "Lat"),
  crs = st_crs(shp),
  remove = FALSE
)

# Peta klaster LISA di atas peta dasar shapefile
ggplot() +
  
  geom_sf(
    data = shp_dasar,
    fill = "grey96",
    color = "grey70",
    linewidth = 0.15
  ) +
  
  geom_sf(
    data = d_sf,
    aes(color = cluster),
    size = 2.4,
    alpha = 0.95
  ) +
  
  scale_color_manual(
    values = cols
  ) +
  
  labs(
    title =
      "Peta Klaster LISA Rata-Rata Lama Sekolah",
    subtitle =
      "Kabupaten/Kota Pulau Sumatera Tahun 2025",
    color = "Klaster LISA"
  ) +
  
  theme_minimal()

Peta ini memperjelas sebaran spasial klaster pada peta Pulau Sumatera. Hot spot (High-High) terkonsentrasi di bagian tengah-barat Sumatera (Sumatera Barat dan sebagian Riau), dengan kantong kecil di pesisir timur Sumatera Utara dan Kepulauan Riau. Cold spot (Low-Low) muncul di Aceh bagian barat, kawasan Tapanuli, serta bagian barat daya pulau (Bengkulu bagian selatan, Lampung Barat–Pesisir Barat, dan OKU). Satu pencilan Low-High (Pasaman Barat) terlihat di tengah kelompok hot spot Sumatera Barat. Pola yang saling terpisah ini menunjukkan bahwa kesenjangan RLS di Sumatera tidak tersebar acak, melainkan mengikuti kedekatan geografis, sehingga kebijakan peningkatan pendidikan lebih efektif jika diprioritaskan pada klaster cold spot dan memperhatikan wilayah tetangganya.

Wilayah yang paling layak diprioritaskan untuk intervensi pendidikan adalah cold spot, terutama yang nilai \(I_i\)-nya besar seperti Kaur, Nagan Raya, Bengkulu Selatan, Tapanuli Utara, dan Aceh Barat.