Pendekatan Ekonometrika Spasial dan Pemodelan Interpolasi untuk Menganalisis Pengaruh PDRB, Luas Wilayah, dan Jumlah Penduduk Usia Sekolah terhadap Jumlah Sekolah Provinsi di Indonesia Tahun 2024

1 Pendahuluan

1.1 Latar Belakang

Jumlah sekolah merupakan indikator penting dalam menggambarkan ketersediaan dan pemerataan layanan pendidikan antarwilayah. Di Indonesia, distribusi jumlah sekolah pada 38 provinsi menunjukkan variasi yang dipengaruhi oleh perbedaan kondisi ekonomi, demografis, dan geografis. Selain itu, karakteristik wilayah yang saling berdekatan secara geografis memungkinkan terjadinya keterkaitan spasial, di mana kondisi suatu provinsi dapat dipengaruhi oleh provinsi di sekitarnya. Oleh karena itu, pendekatan analisis spasial menjadi relevan untuk mengkaji pola persebaran jumlah sekolah serta mengidentifikasi adanya ketergantungan antarwilayah dalam penyediaan fasilitas pendidikan.

Faktor ekonomi yang direpresentasikan oleh Produk Domestik Regional Bruto (PDRB), luas wilayah, dan jumlah penduduk usia sekolah diduga memiliki hubungan yang tidak hanya bersifat langsung, tetapi juga dipengaruhi oleh efek spasial. Provinsi dengan PDRB tinggi cenderung memiliki kapasitas fiskal yang lebih baik dalam pembangunan infrastruktur pendidikan, sementara luas wilayah menentukan kebutuhan sekolah dari sisi jangkauan pelayanan. Di sisi demografis, jumlah penduduk usia sekolah menjadi penentu utama kebutuhan fasilitas pendidikan formal. Namun demikian, pengaruh variabel-variabel tersebut dapat melampaui batas administrasi provinsi, sehingga penting untuk menguji adanya autokorelasi spasial baik secara global maupun lokal guna memahami pola klaster dan wilayah dengan karakteristik serupa.

Melalui pengujian autokorelasi spasial global dan lokal, penelitian ini bertujuan untuk mengidentifikasi pola pengelompokan (cluster) dan wilayah yang menjadi pusat atau kantong kekurangan fasilitas pendidikan. Selanjutnya, pendekatan ekonometrika spasial digunakan untuk mengestimasi pengaruh PDRB, luas wilayah, dan jumlah penduduk usia sekolah terhadap jumlah sekolah dengan mempertimbangkan efek ketergantungan spasial antarprovinsi. Selain itu, pemodelan interpolasi spasial diterapkan untuk memvisualisasikan dan memperkirakan variasi jumlah sekolah secara kontinu di wilayah Indonesia. Dengan pendekatan ini, diharapkan penelitian dapat memberikan gambaran yang lebih komprehensif mengenai distribusi dan determinan jumlah sekolah, serta menjadi dasar perumusan kebijakan pendidikan yang lebih tepat sasaran dan berbasis wilayah.

1.2 Identifikasi Masalah

Berdasarkan latar belakang tersebut, permasalahan utama dalam penelitian ini adalah belum diketahuinya secara jelas pola persebaran spasial jumlah sekolah antarprovinsi di Indonesia serta sejauh mana ketergantungan spasial memengaruhi hubungan antara jumlah sekolah dengan PDRB, luas wilayah, dan jumlah penduduk usia sekolah. Perbedaan karakteristik ekonomi, demografis, dan geografis antarprovinsi diduga menimbulkan ketimpangan penyediaan fasilitas pendidikan, yang tidak dapat dijelaskan secara memadai hanya dengan pendekatan nonspasial. Oleh karena itu, diperlukan identifikasi adanya autokorelasi spasial, baik secara global maupun lokal, serta pemodelan ekonometrika dan interpolasi spasial untuk memahami pola klaster, pengaruh antarwilayah, dan variasi jumlah sekolah secara lebih komprehensif

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk:

1. Mengidentifikasi pola persebaran spasial jumlah sekolah pada 38 provinsi di Indonesia melalui pengujian autokorelasi spasial global.
2. Menentukan wilayah-wilayah yang membentuk klaster dan pencilan (outlier) jumlah sekolah menggunakan analisis autokorelasi spasial lokal.
3. Menganalisis pengaruh Produk Domestik Regional Bruto (PDRB), luas wilayah, dan jumlah penduduk usia sekolah terhadap jumlah sekolah dengan menggunakan pendekatan ekonometrika spasial.
4. Membandingkan hasil model ekonometrika spasial dengan model nonspasial untuk menilai keberadaan dan besarnya efek spasial antarprovinsi.
5. Memvisualisasikan dan memperkirakan variasi jumlah sekolah secara spasial melalui pemodelan interpolasi sebagai dasar pemahaman distribusi fasilitas pendidikan di Indonesia.

1.4 Batasan Penelitian

Penelitian ini dibatasi pada analisis jumlah sekolah sebagai variabel dependen dengan Produk Domestik Regional Bruto (PDRB), luas wilayah, dan jumlah penduduk usia sekolah sebagai variabel independen pada 38 provinsi di Indonesia. Data yang digunakan merupakan data sekunder pada satu periode waktu tertentu sehingga analisis bersifat cross-section dan tidak menangkap dinamika perubahan antarwaktu. Pendekatan spasial yang digunakan meliputi autokorelasi spasial global dan lokal, ekonometrika spasial, serta pemodelan interpolasi, dengan unit analisis pada tingkat provinsi dan keterbatasan pada ketelitian batas administrasi wilayah. Faktor lain di luar variabel penelitian, seperti kualitas pendidikan, kebijakan daerah, dan kondisi sosial budaya, tidak dimasukkan dalam model sehingga berada di luar cakupan penelitian ini.

2 Tinjauan Pustaka

2.1 Konsep Spatial dan Spatial Heterogenity

Keterkaitan spasial (spatial dependence) menunjukkan bahwa nilai suatu variabel di satu wilayah dipengaruhi oleh wilayah lain yang berdekatan secara geografis. Dalam konteks PDRB antarprovinsi di Indonesia, hal ini berarti pertumbuhan ekonomi suatu provinsi dapat terdorong oleh interaksi ekonomi, mobilitas tenaga kerja, keterhubungan infrastruktur, maupun kesamaan kebijakan dengan provinsi sekitarnya. Cartone et al. (2022) menyebutnya sebagai efek limpahan spasial (spatial spillover), yaitu penularan dinamika ekonomi antarwilayah yang membentuk pola pengelompokan (clustering) wilayah ber-PDRB tinggi maupun rendah.

Menurut Jankiewicz et al. (2025), ketergantungan spasial tidak hanya dipicu oleh kedekatan fisik, tetapi juga oleh keterhubungan fungsional seperti arus produksi dan pasar tenaga kerja. Deng et al. (2024) menambahkan bahwa pengaruh spasial menurun seiring jarak, namun tetap signifikan di wilayah yang memiliki hubungan ekonomi erat.

Selain itu, keragaman spasial (spatial heterogeneity) mencerminkan perbedaan karakteristik antarwilayah yang menyebabkan respons ekonomi tidak seragam. Trinh et al. (2024) dan Wu et al. (2024) menegaskan bahwa variasi kapasitas ekonomi, infrastruktur, dan aksesibilitas turut menentukan perbedaan tersebut.

Mengabaikan dependensi dan heterogenitas spasial dalam analisis, seperti dengan menggunakan metode OLS, dapat menyebabkan estimasi bias dan kesimpulan kebijakan keliru. Karena itu, pengujian autokorelasi spasial dan pengakuan terhadap perbedaan struktural antarwilayah menjadi langkah penting dalam membangun model spatial econometrics yang valid dan kontekstual.

2.2 Autokorelasi Spasial

Keterkaitan spasial tidak hanya bersifat teoritis, tetapi juga dapat diuji secara empiris melalui analisis autokorelasi spasial. Autokorelasi spasial menggambarkan hubungan antara nilai suatu variabel di satu wilayah dengan nilai variabel serupa di wilayah sekitarnya. Dalam konteks PDRB antarprovinsi di Indonesia, analisis ini berfungsi untuk mengidentifikasi apakah provinsi dengan PDRB tinggi membentuk klaster pertumbuhan di wilayah tertentu atau tersebar secara acak. Noviyanti et al. (2023) dalam Jurnal Artha Info menemukan bahwa provinsi dengan tingkat kemiskinan rendah cenderung dikelilingi oleh provinsi dengan karakteristik serupa, yang menandakan adanya pola spasial positif.

Dalam analisis spasial, autokorelasi dibedakan menjadi dua pendekatan, yaitu autokorelasi global dan autokorelasi lokal. Autokorelasi global digunakan untuk menilai apakah secara keseluruhan terdapat pola spasial dalam data. Jika nilai statistik global signifikan, maka terdapat indikasi pengelompokan antarwilayah dan pola data tidak terjadi secara acak. Namun, pendekatan ini tidak menunjukkan lokasi spesifik dari klaster atau anomali. Sebaliknya, autokorelasi lokal digunakan untuk mengidentifikasi area tertentu yang memiliki kemiripan nilai, seperti hotspot (wilayah bernilai tinggi yang berdekatan), coldspot (wilayah bernilai rendah yang membentuk klaster), maupun anomali spasial.

Beberapa ukuran yang umum digunakan dalam pengujian autokorelasi spasial antara lain Moran’s I Global, Geary’s C, Local Moran’s I (LISA), serta Getis–Ord Gi. Moran’s I dan Geary’s C mengukur keterkaitan spasial secara keseluruhan, sementara LISA dan Getis–Ord Gi digunakan untuk mendeteksi pola spasial lokal secara lebih detail.

2.2.1 Moran’s I Global

Moran’s I Global adalah ukuran yang paling umum digunakan untuk mendeteksi autokorelasi spasial global pada data. Indeks ini mengukur apakah terdapat pola pengelompokan nilai yang tinggi atau rendah antarwilayah. Nilai Moran’s I yang positif menunjukkan adanya klaster wilayah dengan nilai yang mirip (baik tinggi maupun rendah), sedangkan nilai negatif mengindikasikan pola penyebaran yang bersifat menyebar atau kontras antarwilayah.

Secara matematis, Moran’s I dirumuskan sebagai berikut:

\[ I = \frac{n}{\sum_{i=1}^n \sum_{j=1}^n w_{ij}} \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} \]

di mana \(n\) adalah jumlah wilayah, \(w_{ij}\) adalah elemen matriks bobot spasial yang menunjukkan kedekatan antara wilayah \(i\) dan \(j\), \(x_i\) adalah nilai variabel pada wilayah \(i\), dan \(\bar{x}\) adalah nilai rata-rata variabel pada seluruh wilayah.

2.2.2 Geary’s C

Geary’s C merupakan ukuran autokorelasi spasial global yang lebih sensitif terhadap perbedaan lokal antarwilayah tetangga dibandingkan Moran’s I. Jika Moran’s I lebih memperhatikan kovarians antarwilayah, Geary’s C lebih menekankan perbedaan kuadrat antarwilayah yang berdekatan.

Rumus Geary’s C dituliskan sebagai berikut:

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

---

2.2.3 Local Moran’s I (LISA)

Local Moran’s I atau LISA (Local Indicators of Spatial Association) digunakan untuk mengidentifikasi klaster spasial lokal dan mendeteksi wilayah yang berpotensi menjadi hotspot, coldspot, atau outlier spasial

Formulasinya adalah sebagai berikut:

\[ I_i = \frac{(x_i - \bar{x})}{m_2} \sum_{j=1}^n w_{ij} (x_j - \bar{x}) \]

dengan \[ m_2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n} \]

2.2.4 Getis–Ord Gi*

Indeks Getis–Ord Gi* digunakan untuk mendeteksi hotspot atau coldspot dengan melihat intensitas nilai tinggi atau rendah dalam radius spasial tertentu. Berbeda dengan Local Moran yang mengukur korelasi lokal, Getis–Ord Gi* fokus pada tingkat konsentrasi nilai tinggi atau rendah dalam suatu area

Rumus Getis–Ord Gi* adalah:

\[ G_i^* = \frac{\sum_{j=1}^n w_{ij} x_j - \bar{x} \sum_{j=1}^n w_{ij}} {s \sqrt{\frac{n \sum_{j=1}^n w_{ij}^2 - (\sum_{j=1}^n w_{ij})^2}{n - 1}}} \] ## 2.3 Spatial Econometrics

Spatial econometrics merupakan cabang dari ilmu ekonometrika yang memasukkan dimensi spasial ke dalam model regresi guna menangkap hubungan antarwilayah. Pendekatan ini penting digunakan ketika asumsi independensi antarunit geografis tidak terpenuhi, seperti pada kasus PDRB antarprovinsi di Indonesia, di mana aktivitas ekonomi suatu provinsi kerap dipengaruhi oleh kondisi ekonomi wilayah di sekitarnya. Jika faktor spasial diabaikan, maka hasil estimasi menggunakan Ordinary Least Squares (OLS) berpotensi menjadi bias dan tidak efisien akibat adanya autokorelasi spasial pada residual maupun variabel dependennya (Anselin, 2010).

Model dalam spatial econometrics secara umum terbagi menjadi dua kelompok utama, yakni model yang menekankan pada ketergantungan spasial (spatial dependence) dan model yang menangkap keragaman spasial (spatial heterogeneity). Beberapa model yang umum diterapkan meliputi OLS, Spatial Autoregressive (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Autoregressive Combined (SAC), serta General Nesting Spatial (GNS) model.

2.3 Ekonometrika Spasial

2.3.1 Ordinary Least Squares (OLS)

Model OLS digunakan sebagai pendekatan dasar sebelum mempertimbangkan komponen spasial. Model ini mengasumsikan bahwa antarobservasi bersifat independen, sehingga hubungan antara variabel dependen dan independen dapat dijelaskan melalui persamaan linear klasik:

\[ y = X\beta + \varepsilon \]

dengan \(y\) adalah vektor variabel dependen (misalnya PDRB per provinsi), \(X\) matriks variabel independen (seperti TPAK dan IPM), \(\beta\) parameter yang diestimasi, dan \(\varepsilon\) adalah komponen error yang diasumsikan berdistribusi normal dengan varian homogen dan tanpa autokorelasi spasial. Namun, jika terdapat korelasi spasial antarwilayah, maka model OLS tidak lagi menghasilkan estimasi BLUE (Best Linear Unbiased Estimator).

2.3.2 Spatial Autoregressive Model (SAR)

Model SAR menambahkan efek ketergantungan spasial langsung pada variabel dependen melalui matriks bobot spasial \(W\):

\[ y = \rho W y + X\beta + \varepsilon \]

dengan \(\rho\) adalah koefisien autokorelasi spasial, \(W y\) menunjukkan pengaruh nilai \(y\) dari wilayah sekitar terhadap nilai \(y\) pada suatu wilayah, dan \(\varepsilon\) adalah residual acak. Model ini relevan ketika interaksi antarwilayah bersifat langsung, misalnya pertumbuhan ekonomi suatu provinsi dipengaruhi oleh pertumbuhan ekonomi provinsi tetangga.

2.3.2.1 Estimasi dengan Maximum Likelihood (ML)

Pendekatan estimasi yang paling umum digunakan dalam model SAR adalah Maximum Likelihood Estimation (MLE). Model dapat ditulis ulang menjadi bentuk setara dengan mendefinisikan \(A(\rho) = I - \rho W\):

\[ A(\rho)y = X\beta + \varepsilon \]

dengan \(\varepsilon \sim N(0, \sigma^2 I)\).

Fungsi log-likelihood untuk model SAR dirumuskan sebagai berikut:

\[ \ell(\beta, \rho, \sigma^2) = \log|A(\rho)| - \frac{N}{2}\log\sigma^2 - \frac{1}{2\sigma^2}(A(\rho)y - X\beta)^\top(A(\rho)y - X\beta) \]

Untuk setiap nilai \(\rho\), estimasi parameter \(\beta\) dan \(\sigma^2\) diperoleh sebagai berikut:

\[ \hat{\beta}(\rho) = (X^\top X)^{-1}X^\top A(\rho)y \]

\[ \hat{\sigma}^2(\rho) = \frac{1}{N}\|A(\rho)y - X\hat{\beta}(\rho)\|^2 \]

Sehingga fungsi likelihood terkonsentrasi (concentrated likelihood) dapat ditulis sebagai:

\[ \ell_c(\rho) = \log|A(\rho)| - \frac{N}{2}\log\hat{\sigma}^2(\rho) \]

Nilai \(\rho\) yang memaksimalkan \(\ell_c(\rho)\) memberikan estimasi maximum likelihood untuk parameter autokorelasi spasial.

Interpretasi koefisien \(\rho\) menunjukkan tingkat pengaruh langsung antara nilai variabel dependen suatu wilayah dengan nilai variabel dependen di wilayah sekitarnya. Semakin tinggi \(\rho\), semakin kuat pengaruh antarwilayah yang terjadi.

2.3.3 Spatial Error Model (SEM)

Model SEM tepat digunakan jika hasil uji Moran’s I pada residual OLS menunjukkan adanya autokorelasi spasial, namun teori ekonomi atau sosial tidak mendukung adanya interaksi langsung antarwilayah pada variabel dependen.

Rumus Model:

\[ y = X\beta + u, \quad u = \lambda W u + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2 I) \]

dengan \(y\) sebagai variabel dependen (misalnya PDRB), \(X\) sebagai matriks kovariat (misalnya TPAK dan IPM), \(W\) adalah matriks bobot spasial, \(\lambda\) parameter autokorelasi error, dan \(\varepsilon\) adalah error acak.

2.3.4 Spatial Durbin Model (SDM)

Spatial Durbin Model digunakan ketika ada interaksi endogen pada y dan spillover dari kovariat. Model SDM merupakan perluasan dari SAR, dengan menambahkan pengaruh variabel independen dari wilayah sekitar ke dalam model.

\[ y = \rho W y + X\beta + W X\theta + \varepsilon \]

Parameter \(\theta\) merepresentasikan efek spasial dari variabel independen di wilayah tetangga. Model ini lebih fleksibel karena dapat menangkap efek spillover dari variabel-variabel penjelas, misalnya peningkatan IPM di suatu provinsi dapat berdampak pada peningkatan PDRB di provinsi sekitarnya.

2.3.5 Spatial Autoregressive Combined Model (SAC)

Model SAC menggabungkan unsur SAR dan SEM sekaligus, sehingga dapat menangkap ketergantungan spasial baik pada variabel dependen maupun pada error. SAC digunakan ketika ada interaksi endogen dan juga omitted variables berpola spasial. Berikut model SAC:

\[ y = \rho W y + X\beta + u, \quad u = \lambda W u + \varepsilon \]

2.3.6 General Nesting Spatial Model (GNS)

Model GNS merupakan bentuk paling umum dan komprehensif dalam ekonometrika spasial karena menggabungkan seluruh komponen ketergantungan spasial, baik pada variabel dependen, independen, maupun error:

\[ y = \rho W y + X\beta + W X\theta + u, \quad u = \lambda W u + \varepsilon \]

Model ini memungkinkan peneliti untuk menguji berbagai bentuk hubungan spasial yang kompleks dan saling tumpang tindih. Dalam konteks PDRB antarprovinsi, model ini dapat mengungkapkan apakah pengaruh IPM dan TPAK suatu provinsi tidak hanya berdampak pada wilayahnya sendiri tetapi juga menyebar ke wilayah sekitar.

2.4 Spatial Interpolation

Interpolasi spasial merupakan teknik analisis untuk mengestimasi nilai suatu variabel pada lokasi yang tidak teramati berdasarkan nilai pada lokasi yang diketahui, dengan memanfaatkan prinsip kedekatan spasial. Konsep ini sejalan dengan First Law of Geography yang menyatakan bahwa lokasi yang berdekatan cenderung memiliki karakteristik yang lebih mirip dibandingkan lokasi yang berjauhan (Tobler, 1970). Dalam kerangka analisis spasial, interpolasi berfungsi sebagai alat pendukung untuk memvisualisasikan pola distribusi suatu fenomena secara kontinu, melengkapi analisis autokorelasi spasial global (Moran’s I) dan lokal (LISA).

Dalam penelitian ini, interpolasi spasial digunakan untuk menggambarkan variasi jumlah sekolah antarprovinsi di Indonesia. Meskipun unit analisis bersifat areal, pendekatan interpolasi memungkinkan pemahaman pola spasial yang lebih intuitif dan membantu mengidentifikasi wilayah dengan kecenderungan nilai tinggi atau rendah, yang konsisten dengan hasil pengujian ketergantungan spasial.

2.4.1 Metode Inverse Distance Weighting (IDW)

Inverse Distance Weighting (IDW) merupakan metode interpolasi deterministik yang mengasumsikan bahwa pengaruh suatu lokasi terhadap lokasi lain menurun seiring bertambahnya jarak (Shepard, 1968). Metode ini tidak secara eksplisit memodelkan autokorelasi spasial, sehingga lebih sesuai digunakan sebagai visualisasi awal pola spasial.

Secara matematis, estimasi IDW dirumuskan sebagai: \[ \hat{Z}(x_0) = \frac{\sum_{i=1}^{n} \frac{Z(x_i)}{d_{i0}^p}}{\sum_{i=1}^{n} \frac{1}{d_{i0}^p}} \] dengan \(\hat{Z}(x_0)\) merupakan nilai estimasi pada lokasi \(x_0\), \(Z(x_i)\) adalah nilai variabel pada lokasi ke-\(i\), \(d_{i0}\) adalah jarak antara lokasi \(x_i\) dan \(x_0\), serta \(p\) merupakan parameter pangkat yang mengatur tingkat pengaruh jarak.

3 Metodologi Penelitian

3.1 Data dan Sumber Data

Pada penelitian ini digunakan data sekunder yang bersumber dari situs resmi Badan Pusat Statistik (BPS) Indonesia (https://www.bps.go.id ). Data yang diperoleh dari BPS dianggap memiliki tingkat keandalan yang tinggi karena dihimpun melalui prosedur statistik yang baku dan terverifikasi, sehingga layak dijadikan dasar dalam analisis spasial. Jenis data yang digunakan berupa data cross-sectional tahun 2024 dengan unit analisis mencakup 38 provinsi di Indonesia. Pemilihan jenis data tersebut dilakukan karena penelitian ini berfokus pada pengujian keterkaitan dan ketergantungan spasial antarprovinsi dalam menganalisis faktor-faktor yang memengaruhi jumlah sekolah.

Variabel yang digunakan dalam penelitian ini terdiri atas empat komponen utama. Jumlah sekolah digunakan sebagai variabel dependen yang merepresentasikan ketersediaan fasilitas pendidikan formal di masing-masing provinsi. Produk Domestik Regional Bruto (PDRB) digunakan sebagai variabel independen yang mencerminkan kondisi dan kapasitas ekonomi daerah. Luas wilayah digunakan untuk merepresentasikan karakteristik geografis provinsi yang memengaruhi kebutuhan dan persebaran sarana pendidikan. Sementara itu, jumlah penduduk usia sekolah digunakan sebagai variabel demografis yang menggambarkan besarnya potensi kebutuhan layanan pendidikan di setiap provinsi. Keempat variabel tersebut dianalisis menggunakan pendekatan spasial untuk menangkap pengaruh langsung maupun ketergantungan antarwilayah dalam distribusi jumlah sekolah di Indonesia.

library(sf)
library(dplyr)

shp <- st_read("C:/Users/Johanes Credo/OneDrive/Documents/KULIAH/SEM 5/Spatial/indonesia-38-provinces.geojson", quiet = TRUE)
class(shp)     # harus "sf"

## [1] "sf"         "data.frame"

nrow(shp)      # harus 38

## [1] 38

st_geometry_type(shp)

##  [1] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
##  [6] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [11] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [16] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [21] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [26] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [31] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## [36] MULTIPOLYGON MULTIPOLYGON MULTIPOLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

st_crs(shp)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

if ("Shape" %in% names(shp)) {
  shp <- st_as_sf(shp, wkt = "Shape")
} else {
  shp <- st_as_sf(shp)
}
names(shp)

## [1] "id"        "KODE_PROV" "PROVINSI"  "geometry"

unique(shp$PROVINSI)

##  [1] "Sulawesi Tengah"            "Sulawesi Barat"            
##  [3] "Sulawesi Selatan"           "Papua Tengah"              
##  [5] "Papua Barat"                "Gorontalo"                 
##  [7] "Riau"                       "Papua Selatan"             
##  [9] "Daerah Istimewa Yogyakarta" "Sumatera Barat"            
## [11] "DKI Jakarta"                "Maluku"                    
## [13] "Bengkulu"                   "Lampung"                   
## [15] "Papua"                      "Kepulauan Riau"            
## [17] "Nusa Tenggara Barat"        "Jambi"                     
## [19] "Bali"                       "Jawa Timur"                
## [21] "Papua Barat Daya"           "Sumatera Utara"            
## [23] "Sulawesi Tenggara"          "Nusa Tenggara Timur"       
## [25] "Kalimantan Selatan"         "Aceh"                      
## [27] "Kalimantan Tengah"          "Papua Pegunungan"          
## [29] "Kepulauan Bangka Belitung"  "Sumatera Selatan"          
## [31] "Banten"                     "Sulawesi Utara"            
## [33] "Kalimantan Utara"           "Kalimantan Timur"          
## [35] "Jawa Tengah"                "Maluku Utara"              
## [37] "Kalimantan Barat"           "Jawa Barat"

data <- read.csv(file.choose(), sep=";")
data

##                      Provinsi Jumlah.Sekolah       PDRB Luas.Wilayah
## 1                        Aceh           5414  153780.44     56834.75
## 2              Sumatera Utara          13531  632534.73     72460.74
## 3              Sumatera Barat           5429  199407.38     42119.54
## 4                        Riau           5627  571233.59     89935.90
## 5                       Jambi           3446  176906.50     49026.58
## 6            Sumatera Selatan           6822  379119.63     86771.68
## 7                    Bengkulu           2021   54454.65     20128.34
## 8                     Lampung           6763  281557.20     33570.26
## 9   Kepulauan Bangka Belitung           1147   60802.64     16690.13
## 10             Kepulauan Riau           1590  209939.07      8269.71
## 11                DKI Jakarta           3795 2151041.33       660.98
## 12                 Jawa Barat          27409 1752071.20     37044.86
## 13                Jawa Tengah          22921 1157025.94     34337.49
## 14 Daerah Istimewa Yogyakarta           2469  124590.45      3170.65
## 15                 Jawa Timur          25539 1935810.15     48036.84
## 16                     Banten           6884  531735.25      9352.77
## 17                       Bali           2970  168186.03      5590.15
## 18        Nusa Tenggara Barat           4823  109414.97     19675.89
## 19        Nusa Tenggara Timur           7723   78044.57     46446.64
## 20           Kalimantan Barat           6326  162574.47    147037.04
## 21          Kalimantan Tengah           3778  118682.32    153443.91
## 22         Kalimantan Selatan           3766  156756.94     37135.05
## 23           Kalimantan Timur           2902  570824.01    126981.28
## 24           Kalimantan Utara            755   73007.99     70101.18
## 25             Sulawesi Utara           3173  107575.01     14500.28
## 26            Sulawesi Tengah           4070  212281.54     61605.72
## 27           Sulawesi Selatan           8794  396141.74     45330.55
## 28          Sulawesi Tenggara           3465  113989.49     36159.71
## 29                  Gorontalo           1348   32949.62     12025.15
## 30             Sulawesi Barat           1803   37088.07     16594.75
## 31                     Maluku           2833   37209.36     46158.27
## 32               Maluku Utara           2066   55152.33     32998.70
## 33                Papua Barat            805   49486.23     60275.31
## 34           Papua Barat Daya            821   24873.72     39122.95
## 35                      Papua           1244   51587.16     82680.96
## 36              Papua Selatan            790   18963.85    117849.16
## 37               Papua Tengah            791  105471.48     61072.91
## 38           Papua Pegunungan            954   13798.74     51213.33
##    Usia.Sekolah
## 1       1304.61
## 2       3661.14
## 3       1370.70
## 4       1580.17
## 5        874.70
## 6       2075.54
## 7        496.08
## 8       2212.31
## 9        359.69
## 10       512.77
## 11      2509.48
## 12     11824.17
## 13      8899.46
## 14       882.96
## 15      9820.63
## 16      2919.66
## 17      1041.21
## 18      1326.03
## 19      1328.38
## 20      1337.66
## 21       659.89
## 22      1003.66
## 23       950.23
## 24       173.75
## 25       634.55
## 26       733.19
## 27      2222.59
## 28       655.99
## 29       288.36
## 30       353.04
## 31       456.95
## 32       318.38
## 33       135.91
## 34       147.28
## 35       249.09
## 36       127.32
## 37       345.93
## 38       344.54

shp <- shp %>%
  mutate(
    prov_bps = toupper(PROVINSI)
  )

data <- data %>%
  mutate(
    Provinsi = toupper(Provinsi)
  )

shp_join <- shp %>%
  left_join(data, by = c("prov_bps" = "Provinsi"))
shp_join %>%
  filter(is.na(Jumlah.Sekolah)) %>%
  .[, c("PROVINSI", "prov_bps")]

## Simple feature collection with 0 features and 2 fields
## Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
## Geodetic CRS:  WGS 84
## [1] PROVINSI prov_bps geometry
## <0 rows> (or 0-length row.names)

3.2 Populasi dan Sampel Penelitian

Populasi dalam penelitian ini mencakup seluruh provinsi di Indonesia yang berjumlah 38 provinsi. Setiap provinsi merepresentasikan satuan wilayah administratif yang menjadi unit analisis dalam studi ini. Karena penelitian ini menggunakan pendekatan data cross-sectional tahun 2024, maka seluruh provinsi dilibatkan dalam analisis tanpa melakukan proses pengambilan sampel. Dengan demikian, penelitian ini menggunakan metode sensus penuh (census), di mana seluruh anggota populasi dijadikan sebagai sampel penelitian. Pendekatan ini dipilih untuk memperoleh gambaran yang lebih komprehensif mengenai pola dan ketergantungan spasial antarprovinsi dalam menganalisis pengaruh Produk Domestik Regional Bruto (PDRB), luas wilayah, dan jumlah penduduk usia sekolah terhadap jumlah sekolah di Indonesia.

3.3 Variabel Penelitian dan Definisi Operasional

Variabel penelitian yang digunakan dalam penelitian ini dapat dijelaskan pada tabel berikut.

Jenis Variabel	Nama Variabel	Definisi Operasional	Satuan	Sumber Data
Variabel Dependen (Y)	Jumlah Sekolah	Jumlah total sekolah formal (SD, SMP, dan SMA) yang terdapat di masing-masing provinsi di Indonesia pada tahun 2024.	Unit Sekolah	BPS
Variabel Independen (X₁)	Produk Domestik Regional Bruto (PDRB)	Nilai total output ekonomi dari seluruh aktivitas produksi di suatu provinsi berdasarkan harga konstan tahun 2024.	Miliar Rupiah	BPS
Variabel Independen (X₂)	Luas Wilayah	Luas administratif masing-masing provinsi di Indonesia yang mencerminkan karakteristik geografis wilayah.	Kilometer Persegi (km²)	BPS
Variabel Independen (X₃)	Jumlah Penduduk Usia Sekolah	Jumlah penduduk yang berada pada usia sekolah formal di setiap provinsi pada tahun 2024.	Jiwa	BPS

3.4 Metode Analisis

Proses analisis dalam penelitian ini dilakukan melalui beberapa tahapan utama sebagai berikut:

1. Analisis Deskriptif Spasial
   • Tahap awal dilakukan dengan menyajikan peta tematik (choropleth map) untuk menggambarkan distribusi jumlah sekolah, Produk Domestik Regional Bruto (PDRB), luas wilayah, dan jumlah penduduk usia sekolah di seluruh provinsi di Indonesia.
     
   • Tujuan tahap ini adalah untuk memberikan gambaran awal mengenai pola persebaran spasial, tingkat ketimpangan antarprovinsi, serta indikasi visual adanya pengelompokan wilayah (spatial clustering).
2. Uji Autokorelasi Spasial Global dan Lokal
   • Uji autokorelasi spasial dilakukan untuk mengidentifikasi adanya ketergantungan spasial antarprovinsi dalam variabel jumlah sekolah.
     
   • Pengujian dilakukan menggunakan Moran’s I untuk mendeteksi autokorelasi spasial secara global dan Local Indicators of Spatial Association (LISA) untuk mengidentifikasi pola lokal seperti high-high, low-low, serta wilayah pencilan (spatial outliers).
     
   • Matriks pembobot spasial (\(W\)) dibentuk menggunakan pendekatan queen contiguity, di mana dua provinsi dianggap bertetangga apabila memiliki sisi atau titik batas wilayah yang bersinggungan.
3. Pemodelan Regresi dan Ekonometrika Spasial
   • Analisis diawali dengan estimasi model regresi linier klasik (Ordinary Least Squares / OLS) untuk menganalisis pengaruh PDRB, luas wilayah, dan jumlah penduduk usia sekolah terhadap jumlah sekolah.
     
   • Selanjutnya dilakukan pengujian autokorelasi spasial pada residual model OLS untuk mendeteksi keberadaan efek spasial yang tidak tertangkap oleh model nonspasial.
     
   • Apabila ditemukan indikasi ketergantungan spasial, maka analisis dilanjutkan dengan estimasi beberapa model ekonometrika spasial, yaitu:
     • Spatial Lag Model (SAR), yang mempertimbangkan pengaruh ketergantungan antarprovinsi pada variabel dependen (jumlah sekolah).
       
     • Spatial Error Model (SEM), yang menangkap efek spasial melalui komponen galat.
       
     • Spatial Durbin Model (SDM), yang mengakomodasi efek lag spasial pada variabel dependen dan variabel independen.
       
     • Spatial Autoregressive Combined Model (SAC), yang menggabungkan efek spasial pada lag dependen dan error.
       
     • General Nesting Spatial Model (GNS), sebagai model paling umum yang mencakup seluruh bentuk ketergantungan spasial.
       
   • Estimasi model dilakukan menggunakan fungsi lm() untuk OLS serta fungsi lagsarlm(), errorsarlm(), dan lagsarlm(..., type = "Durbin") dari paket spdep atau spatialreg di R.
4. Pemilihan Model Spasial Terbaik
   • Tahap ini bertujuan untuk menentukan model ekonometrika spasial yang paling sesuai dengan karakteristik data.
     
   • Pemilihan model dilakukan melalui Lagrange Multiplier (LM) test, baik LM-lag maupun LM-error, serta dengan membandingkan nilai Akaike Information Criterion (AIC) antar model.
     
   • Model dengan nilai AIC terendah dan hasil pengujian yang signifikan dipilih sebagai model terbaik dalam merepresentasikan hubungan spasial jumlah sekolah antarprovinsi.
5. Analisis Efek Langsung dan Efek Tidak Langsung (Spillover Effect)
   • Setelah model terbaik diperoleh, dilakukan analisis dampak menggunakan fungsi impacts() untuk memisahkan efek langsung (pengaruh variabel independen dalam provinsi itu sendiri) dan efek tidak langsung atau spillover effect ke provinsi lain yang bertetangga.
     
   • Analisis ini bertujuan untuk memahami bagaimana perubahan PDRB, luas wilayah, dan jumlah penduduk usia sekolah di suatu provinsi memengaruhi jumlah sekolah baik secara lokal maupun lintas wilayah.
6. Pemodelan Interpolasi Spasial
   • Tahap akhir dilakukan pemodelan interpolasi spasial untuk memvisualisasikan variasi jumlah sekolah secara kontinu di wilayah Indonesia.
     
   • Metode interpolasi yang digunakan meliputi Inverse Distance Weighting (IDW) sebagai pendekatan deterministik dan Kriging sebagai pendekatan geostatistik yang mempertimbangkan autokorelasi spasial.
     
   • Hasil interpolasi digunakan sebagai alat pendukung untuk memperkuat interpretasi pola spasial yang diperoleh dari analisis Moran’s I, LISA, dan model ekonometrika spasial.

Metode analisis ini dipilih karena mampu menangkap hubungan spasial yang kompleks, baik secara global maupun lokal, serta memberikan pemahaman yang lebih komprehensif mengenai distribusi dan faktor-faktor yang memengaruhi jumlah sekolah antarprovinsi di Indonesia pada tahun 2024.

3.5 Alur Penelitian

library(DiagrammeR)

grViz("
digraph flowchart {
graph [layout = dot, rankdir = TB]  # TB = top to bottom (vertikal)

node [shape = box, style = filled, color = black, fontname = Helvetica, fontsize = 11, width = 3]

A [label = 'Pengumpulan Data\n(BPS & Data Spasial Provinsi)', fillcolor = '#a3c4f3']
B [label = 'Pengolahan Data Awal\n(Cleaning, Transformasi,\nStandarisasi)', fillcolor = '#b9e0a5']
C [label = 'Analisis Deskriptif Spasial\n(Peta Tematik & Statistik Ringkas)', fillcolor = '#fbe6a2']
D [label = 'Uji Autokorelasi Spasial\n(Moran’s I dan LISA)', fillcolor = '#f8c0a3']
E [label = 'Estimasi Model\nRegresi & Ekonometrika Spasial\n(OLS, SAR, SEM, SDM, SAC, GNS)', fillcolor = '#d3b3f2']
F [label = 'Pemilihan Model Terbaik\n(LM Test & AIC)', fillcolor = '#f4b4e1']
G [label = 'Analisis Efek Spasial\n(Efek Langsung & Spillover)', fillcolor = '#f4b4e1']
H [label = 'Pemodelan Interpolasi Spasial\n(IDW & Kriging)', fillcolor = '#cdb4db']
I [label = 'Kesimpulan &\nImplikasi Kebijakan Pendidikan', fillcolor = '#a8dee0']

A -> B -> C -> D -> E -> F -> G -> H -> I
}
")

4 Hasil dan Pembahasan

4.1 Peta Deskriptif

library(tmap)
tmap_mode("plot")

4.1.1 Peta Jumlah Sekolah

tm_shape(shp_join) +
  tm_polygons("Jumlah.Sekolah", palette = "YlGnBu", style = "quantile", 
              title = "Jumlah Sekolah (SD-SMA)") +
  tm_layout(title = "Peta Sebaran Jumlah Sekolah Indonesia 2024", legend.outside = TRUE)

Berdasarkan visualisasi peta di atas, dapat dilihat bahwa persebaran jumlah sekolah belum sepenuhnya merata di Indonesia. Pada Indonesia bagian barat sebenarnya penyebaran sekolah sudah cukup baik dimana secara jumlah didominan menegah ke atas. Namun, jika melihat ke Indonesia bagian timur terdapat ketimpangan persebaran jumlah sekolah. Jumlah sekolah yang ada dominan menengah ke bawah. ### Peta PDRB

tm_shape(shp_join) +
  tm_polygons("PDRB", palette = "YlGnBu", style = "quantile", 
              title = "PDRB") +
  tm_layout(title = "Peta Sebaran PDRB Indonesia 2024", legend.outside = TRUE)

Berdasarkan peta sebaran PDRB di Indonesia, dapat dilihat bahwa juga terdapat ketimpangan PDRB antar provinsi di Indonesia. Pada pulau Jawa dan Sumatra, PDRB tiap provinsi nya dominan cukup besar karena memang banyak dijadikan daerah industrial dan pusat perekonomian. Namun, pada pulau lainnya, PDRB cenderung cukup kecil. ### Peta Luas Wilayah

tm_shape(shp_join) +
  tm_polygons("Luas.Wilayah", palette = "YlGnBu", style = "quantile", 
              title = "Luas Wilayah") +
  tm_layout(title = "Peta Luas Wilayah Indonesia 2024", legend.outside = TRUE)

Peta ini menggambarkan luas wilayah tiap provinsi di Indonesia. Perbedaan luas merupakan hal yang sangat wajar karena merupakan kondisi geografis faktual. ### Peta Penduduk Usia Sekolah

tm_shape(shp_join) +
  tm_polygons("Usia.Sekolah", palette = "YlGnBu", style = "quantile", 
              title = "Penduduk Usia Sekolah") +
  tm_layout(title = "Peta Sebaran Jumlah Penduduk Usia Sekolah Indonesia 2024", legend.outside = TRUE)

Peta ini menggambarkan kondisi jumlah penduduk yang termasuk dalam usia sekolah yaitu 6-18 tahun. Dimana ketimpangan jumlah merupakan hal yang cukup wajar karena juga dipengaruhi oleh jumlah penduduk secara umum di wilayah tersebut, yang mana jumlah penduduk antar provinsi nya sangat beragam.

4.2 Autokorelasi Spasial

4.2.1 Pembobotan dengan k-NN

library(spdep)
crs_indonesia <- "+proj=aea +lat_1=-5 +lat_2=10 +lat_0=0 +lon_0=120 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
shp_proj <- st_transform(shp_join, crs_indonesia)
coords <- st_coordinates(st_centroid(shp_proj))

k_vals <- 2:9
I_vals <- numeric(length(k_vals))

for (i in seq_along(k_vals)) {
  k <- k_vals[i]
  nb <- knn2nb(knearneigh(coords, k))
  lw <- nb2listw(nb, style = "W")
  I_vals[i] <- moran.test(shp_join$Jumlah.Sekolah, lw)$estimate[1]
  
  cat(sprintf("k = %d: Moran's I = %.4f, p-value = %.4f\n", 
              k, moran.test(shp_join$Jumlah.Sekolah, lw)$estimate[1], moran.test(shp_join$Jumlah.Sekolah,lw)$p.value))
}

## k = 2: Moran's I = 0.1161, p-value = 0.1443
## k = 3: Moran's I = 0.3378, p-value = 0.0004
## k = 4: Moran's I = 0.2504, p-value = 0.0015
## k = 5: Moran's I = 0.2771, p-value = 0.0001
## k = 6: Moran's I = 0.2273, p-value = 0.0003
## k = 7: Moran's I = 0.2076, p-value = 0.0002
## k = 8: Moran's I = 0.2100, p-value = 0.0001
## k = 9: Moran's I = 0.1927, p-value = 0.0001

delta_I <- diff(I_vals)
k_opt <- k_vals[which.min(abs(delta_I)) + 1]
cat("\nPerubahan Moran's I (delta):\n")

## 
## Perubahan Moran's I (delta):

for (i in 1:length(delta_I)) {
  cat(sprintf("ΔI(k%d→k%d) = %.4f\n", k_vals[i], k_vals[i+1], delta_I[i]))
}

## ΔI(k2→k3) = 0.2217
## ΔI(k3→k4) = -0.0874
## ΔI(k4→k5) = 0.0267
## ΔI(k5→k6) = -0.0497
## ΔI(k6→k7) = -0.0197
## ΔI(k7→k8) = 0.0024
## ΔI(k8→k9) = -0.0173

# Pastikan k dalam range
k_opt <- min(max(k_opt, 3), 9)

k <- k_opt
k

## [1] 8

knn <- knearneigh(coords, k = k)
nb  <- knn2nb(knn)

#nb <- poly2nb(shp_join, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

Dari hasil di atas, didapatkan bahwa k optimum adalah 8. Hal ini dikarenakan saat k=8 perubahan nilai Moran’s I yang terjadi memiliki nilai paling kecil, sehingga Moran’s I paling stabil di saat k=8.

4.2.2 Autokorelasi Global

Moran’s I

safe_try <- function(expr) tryCatch(eval(expr), error = function(e) { warning(e$message); NULL })

moran <- safe_try(quote(moran.test(shp_join$Jumlah.Sekolah, lw, zero.policy = TRUE)))
moran

## 
##  Moran I test under randomisation
## 
## data:  shp_join$Jumlah.Sekolah  
## weights: lw    
## 
## Moran I statistic standard deviate = 3.8204, p-value = 6.661e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.210036119      -0.027027027       0.003850348

Berdasarkan hasil Moran’s I di atas, didapatkan bahwa nilai statistik Moran menunjukkan nilai positf yang artinya menunjukkan terjadi autokorelasi positif secara global pada persebaran Jumlah Sekolah di Indonesia. Dari nilai p-value = 6.661e-05 menunjukkan bahwa autokorelasi global Moran’s I signifikan secara statistik. Geary C

geary <- safe_try(quote(geary.test(shp_join$Jumlah.Sekolah, lw, zero.policy = TRUE)))
geary

## 
##  Geary C test under randomisation
## 
## data:  shp_join$Jumlah.Sekolah 
## weights: lw   
## 
## Geary C statistic standard deviate = 2.8352, p-value = 0.00229
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.736174651       1.000000000       0.008659047

Berdasarkan hasil Geary C di atas, didapatkan bahwa nilai statistik Geary C < 1 yang menunjukkan terdapat autokorelasi spasial positif pada sebaran Jumlah Sekolah di Indonesia. Nilai p-value = 0.00229 juga menunjukkan signifikansi Geary C.

4.2.3 Autokorelasi Lokal

lisa <- tryCatch(
  localmoran(shp_join$Jumlah.Sekolah, lw, zero.policy = TRUE),
  error = function(e) NULL
)

# Tambahkan hasil ke data
if(!is.null(lisa)) {
  shp_join$LISA_I <- lisa[, 1]   # Local Moran's I
  shp_join$LISA_p <- lisa[, 5]   # p-value
  shp_join$lag_Y  <- lag.listw(
    lw,
    shp_join$Jumlah.Sekolah,
    zero.policy = TRUE
  )
} else {
  shp_join$LISA_I <- NA
  shp_join$LISA_p <- NA
  shp_join$lag_Y  <- NA
}

shp_join$LISA_cluster <- NA

# Signifikansi
sig <- shp_join$LISA_p < 0.05

# Rata-rata global
mean_y   <- mean(shp_join$Jumlah.Sekolah, na.rm = TRUE)
mean_lag <- mean(shp_join$lag_Y, na.rm = TRUE)

# High-High
shp_join$LISA_cluster[
  sig & shp_join$Jumlah.Sekolah > mean_y & shp_join$lag_Y > mean_lag
] <- "High-High"

# Low-Low
shp_join$LISA_cluster[
  sig & shp_join$Jumlah.Sekolah < mean_y & shp_join$lag_Y < mean_lag
] <- "Low-Low"

# High-Low
shp_join$LISA_cluster[
  sig & shp_join$Jumlah.Sekolah > mean_y & shp_join$lag_Y < mean_lag
] <- "High-Low"

# Low-High
shp_join$LISA_cluster[
  sig & shp_join$Jumlah.Sekolah < mean_y & shp_join$lag_Y > mean_lag
] <- "Low-High"

# Tidak signifikan
shp_join$LISA_cluster[!sig] <- "Not significant"

# Jadikan faktor
shp_join$LISA_cluster <- factor(
  shp_join$LISA_cluster,
  levels = c(
    "High-High",
    "Low-Low",
    "High-Low",
    "Low-High",
    "Not significant"
  )
)

if (!inherits(shp_join, "sf")) {
  shp_join <- st_as_sf(shp_join)
}

Visualisasi LISA

tmap_mode("plot")

tm_shape(shp_join) +
  tm_polygons(
    "LISA_cluster",
    palette = c("red", "blue", "pink", "lightblue", "grey80"),
    title = "LISA Cluster (Local Moran's I)"
  ) +
  tm_layout(
    main.title = "Peta Klaster LISA – Jumlah Sekolah Provinsi Indonesia",
    legend.outside = TRUE
  )

Berdasarkan hasil visualisasi LISA di atas, didapatkan bahwa sebaran jumlah sekolah di Indonesia memiliki autokorelasi spasial lokal yang signifikan. Terdapat beberapa provinsi di Pulau Jawa yang menunjukkan “High-High” artinya provinsi dengan jumlah sekolah yang tinggi dikelilingi dengan provinsi jumlah sekolah tinggi juga. Pada Pulau Jawa juga terdapat beberapa provinsi yang menunjukkan “Low-High” artinya provinsi dengan jumlah sekolah rendah dikelilingi dengan provinsi jumlah sekolah tinggi. Sementara itu, pada pulau Papua terdapat wilayah yang terindikasi “Low-Low” artinya provinsi dengan jumlah sekolah rendah dikelilingi dengan provnsi jumlah sekolah rendah juga.

4.3 Ekonometrika Spasial

4.3.1 Pemodelan OLS

ols  <- safe_try(quote(lm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join)))

if(!is.null(ols)) {
  cat("\n### 🔹 Model OLS\n")
  print(summary(ols))
} else cat("\nModel OLS gagal dijalankan.\n")

## 
## ### 🔹 Model OLS
## 
## Call:
## lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2523.0  -764.1  -109.2   602.8  3189.7 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.279e+02  3.618e+02   2.288   0.0285 *  
## PDRB         -1.686e-03  6.463e-04  -2.609   0.0134 *  
## Luas.Wilayah  1.072e-02  5.197e-03   2.062   0.0469 *  
## Usia.Sekolah  2.678e+00  1.282e-01  20.880   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1203 on 34 degrees of freedom
## Multiple R-squared:  0.9686, Adjusted R-squared:  0.9658 
## F-statistic: 349.1 on 3 and 34 DF,  p-value: < 2.2e-16

Berdasarkan hasil di atas, didapatkan bahwa model OLS telah dibentuk, dimana semua variabel prediktor berpengaruh signifikan terhadap Y (Jumlah Sekolah), dan juga dengan R-Squared yang sangat baik yaitu 96.58%. Namun, karena Jumlah Sekolah terindikasi autokorelasi spasial positif perlu dicek apakah model OLS sudah cukup atau membutuhkan model spasial untuk memodelkannya. Pengecekan dilakukan dengan melakukan Moran Test pada residual OLS. ### Pengecekan OLS

moran.test(residuals(ols), lw)

## 
##  Moran I test under randomisation
## 
## data:  residuals(ols)  
## weights: lw    
## 
## Moran I statistic standard deviate = 3.3406, p-value = 0.000418
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.192471795      -0.027027027       0.004317431

Berdasarkan hasil Moran’s I Test pada residual OLS, didapatkan bahwa nilai Moran’s I yang positif dan p-value = 0.0004 < 0.05 menunjukkan bahwa model OLS tidaklah cukup untuk memodelkan, dan diperlukan model spasial lainnya untuk memodelkannya. ### LM-Test

library(spdep)
lm_test <- lm.LMtests(ols, lw, test = "all", zero.policy = TRUE)
print(lm_test)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah +
## Usia.Sekolah, data = shp_join)
## test weights: listw
## 
## RSerr = 6.1798, df = 1, p-value = 0.01292
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah +
## Usia.Sekolah, data = shp_join)
## test weights: listw
## 
## RSlag = 0.078617, df = 1, p-value = 0.7792
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah +
## Usia.Sekolah, data = shp_join)
## test weights: listw
## 
## adjRSerr = 6.8645, df = 1, p-value = 0.008793
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah +
## Usia.Sekolah, data = shp_join)
## test weights: listw
## 
## adjRSlag = 0.76333, df = 1, p-value = 0.3823
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah +
## Usia.Sekolah, data = shp_join)
## test weights: listw
## 
## SARMA = 6.9431, df = 2, p-value = 0.03107

Jenis Uji	Statistik Uji	df	p-value	Keterangan
LM-Error (RSerr)	6.1798	1	0.01292	Signifikan (efek error)
LM-Lag (RSlag)	0.0786	1	0.77920	Tidak signifikan
Robust LM-Error (adjRSerr)	6.8645	1	0.00879	Signifikan (efek error)
Robust LM-Lag (adjRSlag)	0.7633	1	0.38230	Tidak signifikan
LM-SARMA	6.9431	2	0.03107	Signifikan (dependensi spasial)

Berdasarkan hasil uji Lagrange Multiplier, LM-Error dan Robust LM-Error menunjukkan nilai p-value yang signifikan, sedangkan LM-Lag dan Robust LM-Lag tidak signifikan. Hasil ini mengindikasikan bahwa ketergantungan spasial dalam data jumlah sekolah antarprovinsi lebih dominan terjadi pada komponen galat. Dengan demikian, model Spatial Error Model (SEM) lebih tepat digunakan dibandingkan Spatial Lag Model (SAR).

4.3.2 Pembentukan Model Spasial

# SAR
library(spatialreg)
sar  <- safe_try(quote(lagsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, zero.policy = TRUE)))

if(!is.null(sar)) {
  cat("\n### 🔹 Model SAR (Spatial Autoregressive)\n")
  print(summary(sar))
} else cat("\nModel SAR gagal dijalankan.\n")

## 
## ### 🔹 Model SAR (Spatial Autoregressive)
## 
## Call:lagsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2549.61  -802.85  -119.52   609.22  3145.82 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   9.4988e+02  5.3395e+02  1.7790  0.075245
## PDRB         -1.6547e-03  6.1999e-04 -2.6690  0.007608
## Luas.Wilayah  1.0248e-02  5.1540e-03  1.9883  0.046777
## Usia.Sekolah  2.6851e+00  1.2373e-01 21.7007 < 2.2e-16
## 
## Rho: -0.02317, LR test value: 0.083438, p-value: 0.77269
## Approximate (numerical Hessian) standard error: 0.07817
##     z-value: -0.2964, p-value: 0.76692
## Wald statistic: 0.087854, p-value: 0.76692
## 
## Log likelihood: -321.275 for lag model
## ML residual variance (sigma squared): 1291500, (sigma: 1136.4)
## Number of observations: 38 
## Number of parameters estimated: 6 
## AIC: 654.55, (AIC for lm: 652.63)

# SEM
sem  <- safe_try(quote(errorsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, zero.policy = TRUE)))

if(!is.null(sem)) {
  cat("\n### 🔹 Model SEM (Spatial Error)\n")
  print(summary(sem))
} else cat("\nModel SEM gagal dijalankan.\n")

## 
## ### 🔹 Model SEM (Spatial Error)
## 
## Call:errorsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2665.289  -493.706   -23.883   437.530  2969.634 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   5.9836e+02  5.4721e+02  1.0935  0.274188
## PDRB         -1.4732e-03  5.1389e-04 -2.8668  0.004146
## Luas.Wilayah  1.3505e-02  4.7615e-03  2.8362  0.004566
## Usia.Sekolah  2.6777e+00  1.0528e-01 25.4340 < 2.2e-16
## 
## Lambda: 0.63852, LR test value: 5.7152, p-value: 0.016818
## Approximate (numerical Hessian) standard error: 0.20773
##     z-value: 3.0738, p-value: 0.0021135
## Wald statistic: 9.4482, p-value: 0.0021135
## 
## Log likelihood: -318.4591 for error model
## ML residual variance (sigma squared): 1045500, (sigma: 1022.5)
## Number of observations: 38 
## Number of parameters estimated: 6 
## AIC: 648.92, (AIC for lm: 652.63)

# SDM
sdm  <- safe_try(quote(lagsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, type = "mixed", zero.policy = TRUE)))

if(!is.null(sdm)) {
  cat("\n### 🔹 Model SDM (Spatial Durbin)\n")
  print(summary(sdm))
} else cat("\nModel SDM gagal dijalankan.\n")

## 
## ### 🔹 Model SDM (Spatial Durbin)
## 
## Call:lagsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2493.21  -406.76  -147.00   496.25  2908.76 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                     Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)       1.7909e+03  8.5032e+02  2.1062  0.035188
## PDRB             -1.4517e-03  6.1818e-04 -2.3483  0.018858
## Luas.Wilayah      1.4630e-02  4.8478e-03  3.0179  0.002545
## Usia.Sekolah      2.6502e+00  1.2485e-01 21.2274 < 2.2e-16
## lag.PDRB          5.5481e-04  2.6067e-03  0.2128  0.831453
## lag.Luas.Wilayah -3.6641e-02  1.5237e-02 -2.4046  0.016188
## lag.Usia.Sekolah -1.7651e+00  8.3001e-01 -2.1266  0.033453
## 
## Rho: 0.62147, LR test value: 5.4584, p-value: 0.019475
## Approximate (numerical Hessian) standard error: 0.21279
##     z-value: 2.9206, p-value: 0.0034932
## Wald statistic: 8.5301, p-value: 0.0034932
## 
## Log likelihood: -316.609 for mixed model
## ML residual variance (sigma squared): 952590, (sigma: 976)
## Number of observations: 38 
## Number of parameters estimated: 9 
## AIC: 651.22, (AIC for lm: 654.68)

# SAC
sac  <- safe_try(quote(sacsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, zero.policy = TRUE)))

if(!is.null(sac)) {
  cat("\n### 🔹 Model SAC (Spatial Autoregressive Combined)\n")
  print(summary(sac))
} else cat("\nModel SAC gagal dijalankan.\n")

## 
## ### 🔹 Model SAC (Spatial Autoregressive Combined)
## 
## Call:sacsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2557.412  -447.668   -60.259   410.492  2865.895 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   1.0260e+03  8.2785e+02  1.2393  0.215231
## PDRB         -1.4340e-03  5.1046e-04 -2.8092  0.004967
## Luas.Wilayah  1.2690e-02  4.9669e-03  2.5548  0.010624
## Usia.Sekolah  2.6790e+00  1.0431e-01 25.6846 < 2.2e-16
## 
## Rho: -0.079583
## Approximate (numerical Hessian) standard error: 0.11153
##     z-value: -0.71353, p-value: 0.47552
## Lambda: 0.6658
## Approximate (numerical Hessian) standard error: 0.20679
##     z-value: 3.2198, p-value: 0.001283
## 
## LR test value: 6.2514, p-value: 0.043907
## 
## Log likelihood: -318.1911 for sac model
## ML residual variance (sigma squared): 1022500, (sigma: 1011.2)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 650.38, (AIC for lm: 652.63)

# GNS
gns  <- safe_try(quote(sacsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, type = "sacmixed", zero.policy = TRUE)))

if(!is.null(gns)) {
  cat("\n### 🔹 Model GNS (General Nesting Spatial Model)\n")
  print(summary(gns))
} else cat("\nModel GNS gagal dijalankan.\n")

## 
## ### 🔹 Model GNS (General Nesting Spatial Model)
## 
## Call:sacsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, type = "sacmixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2298.54  -456.67  -126.29   401.04  2817.75 
## 
## Type: sacmixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                     Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)       3.2661e+03  1.3954e+03  2.3406  0.019255
## PDRB             -1.4834e-03  6.5582e-04 -2.2619  0.023703
## Luas.Wilayah      1.3546e-02  4.5751e-03  2.9609  0.003068
## Usia.Sekolah      2.6537e+00  1.3043e-01 20.3450 < 2.2e-16
## lag.PDRB         -9.3915e-04  3.2199e-03 -0.2917  0.770541
## lag.Luas.Wilayah -4.2125e-02  1.6575e-02 -2.5414  0.011041
## lag.Usia.Sekolah -1.5720e-01  1.3843e+00 -0.1136  0.909587
## 
## Rho: -0.019792
## Approximate (numerical Hessian) standard error: 0.43856
##     z-value: -0.045129, p-value: 0.964
## Lambda: 0.68141
## Approximate (numerical Hessian) standard error: 0.2329
##     z-value: 2.9258, p-value: 0.0034357
## 
## LR test value: 12.253, p-value: 0.031475
## 
## Log likelihood: -315.19 for sacmixed model
## ML residual variance (sigma squared): 869540, (sigma: 932.49)
## Number of observations: 38 
## Number of parameters estimated: 10 
## AIC: 650.38, (AIC for lm: 652.63)

Pembentukan model spasial di atas digunakan untuk memilih model terbaik berdasarkan nilai AIC nya. Model terbaik dipilih berdasarkan nilai AIC terkecil. ### Pemilihan Model Spasial Terbaik

safe_AIC <- function(m) { if(is.null(m)) return(NA); tryCatch(AIC(m), error = function(e) NA) }

models <- list(
  OLS = if (exists("ols")) ols else NULL,
  SAR = if (exists("sar")) sar else NULL,
  SEM = if (exists("sem")) sem else NULL,
  SDM = if (exists("sdm")) sdm else NULL,
  SAC = if (exists("sac")) sac else NULL,
  GNS = if (exists("gns")) gns else NULL
)

model_table <- data.frame(
  Model = names(models),
  AIC = sapply(models, safe_AIC),
  row.names = NULL
)
cat("### 🔹 Perbandingan AIC dan BIC antar model\n\n")

## ### 🔹 Perbandingan AIC dan BIC antar model

print(model_table)

##   Model      AIC
## 1   OLS 652.6335
## 2   SAR 654.5501
## 3   SEM 648.9183
## 4   SDM 651.2180
## 5   SAC 650.3822
## 6   GNS 650.3801

Berdasarakan hasil perbandingan nilai AIC, didapatkan bahwa model spasial terbaik adalah SEM dengan AIC terkecil yaitu 648.9183.

Model SEM

# SEM
sem  <- safe_try(quote(errorsarlm(Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, data = shp_join, listw = lw, zero.policy = TRUE)))

if(!is.null(sem)) {
  cat("\n### 🔹 Model SEM (Spatial Error)\n")
  print(summary(sem))
} else cat("\nModel SEM gagal dijalankan.\n")

## 
## ### 🔹 Model SEM (Spatial Error)
## 
## Call:errorsarlm(formula = Jumlah.Sekolah ~ PDRB + Luas.Wilayah + Usia.Sekolah, 
##     data = shp_join, listw = lw, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2665.289  -493.706   -23.883   437.530  2969.634 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   5.9836e+02  5.4721e+02  1.0935  0.274188
## PDRB         -1.4732e-03  5.1389e-04 -2.8668  0.004146
## Luas.Wilayah  1.3505e-02  4.7615e-03  2.8362  0.004566
## Usia.Sekolah  2.6777e+00  1.0528e-01 25.4340 < 2.2e-16
## 
## Lambda: 0.63852, LR test value: 5.7152, p-value: 0.016818
## Approximate (numerical Hessian) standard error: 0.20773
##     z-value: 3.0738, p-value: 0.0021135
## Wald statistic: 9.4482, p-value: 0.0021135
## 
## Log likelihood: -318.4591 for error model
## ML residual variance (sigma squared): 1045500, (sigma: 1022.5)
## Number of observations: 38 
## Number of parameters estimated: 6 
## AIC: 648.92, (AIC for lm: 652.63)

Berdasarkan hasil estimasi Spatial Error Model (SEM) dengan variabel dependen Jumlah Sekolah dan variabel independen PDRB, Luas Wilayah, serta Jumlah Penduduk Usia Sekolah, diperoleh bukti kuat adanya ketergantungan spasial pada komponen galat. Hal ini ditunjukkan oleh nilai parameter Lambda sebesar 0,63852 yang signifikan secara statistik (p-value = 0,00211). Nilai lambda yang positif dan signifikan mengindikasikan bahwa faktor-faktor tidak teramati yang memengaruhi jumlah sekolah di suatu provinsi memiliki keterkaitan spasial dengan provinsi-provinsi sekitarnya.

Secara parsial, variabel PDRB berpengaruh negatif dan signifikan terhadap jumlah sekolah, dengan koefisien sebesar -0,00147 dan p-value 0,00415. Temuan ini menunjukkan bahwa peningkatan PDRB tidak selalu diikuti oleh peningkatan jumlah sekolah, yang dapat mengindikasikan bahwa daerah dengan tingkat aktivitas ekonomi tinggi cenderung memiliki efisiensi atau konsolidasi fasilitas pendidikan. Sementara itu, Luas Wilayah berpengaruh positif dan signifikan terhadap jumlah sekolah (koefisien = 0,01350, p-value = 0,00457), yang menunjukkan bahwa provinsi dengan wilayah geografis yang lebih luas membutuhkan lebih banyak unit sekolah untuk menjangkau penduduknya.

Variabel Jumlah Penduduk Usia Sekolah memiliki pengaruh positif dan sangat signifikan terhadap jumlah sekolah, dengan koefisien sebesar 2,67765 dan p-value jauh di bawah 0,01. Hasil ini sejalan dengan teori permintaan layanan pendidikan, di mana semakin besar jumlah penduduk usia sekolah di suatu wilayah, maka kebutuhan akan fasilitas pendidikan juga semakin meningkat.

Dari sisi kinerja model, nilai AIC SEM sebesar 648,92 lebih rendah dibandingkan dengan AIC model OLS sebesar 652,63, yang menandakan bahwa SEM memberikan kecocokan model yang lebih baik. Selain itu, hasil Likelihood Ratio (LR) test yang signifikan (p-value = 0,01682) memperkuat bukti bahwa penggunaan model spasial lebih tepat dibandingkan regresi linier klasik. Dengan demikian, dapat disimpulkan bahwa Spatial Error Model (SEM) merupakan model yang paling sesuai untuk menganalisis pola jumlah sekolah antarprovinsi di Indonesia, karena mampu menangkap efek spasial yang tidak terobservasi secara langsung. ## Pemodelan Interpolasi ### Pembentukan Centroid Data Area

# Cek sistem koordinat awal
st_crs(shp_join)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

# Definisikan CRS untuk Indonesia (baik untuk analisis area)
crs_indonesia <- "+proj=aea +lat_1=-5 +lat_2=10 +lat_0=0 +lon_0=120 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# Transformasi ke CRS yang sesuai
shp_proj <- st_transform(shp_join, crs_indonesia)

# Ambil centroid dari setiap polygon
centroid <- st_point_on_surface(shp_proj)
nrow(centroid) # Cek jumlah titik

## [1] 38

centroid_sp <- as(centroid, "Spatial")

Visualisasi Titik Centroid

tm_shape(centroid) +
  tm_dots(col = "Jumlah.Sekolah", 
          palette = "viridis", 
          size = 0.2,
          title = "Jumlah Sekolah per Wilayah") +
  tm_layout(legend.outside = TRUE)

Dari visualisasi tersebut, didapatkan bahwa tiap provinsi sekarang diwakilkan oleh satu titik representatif dari masing-masing wilayah provinsi. ### Pembentukan Grid dan Centroid

# Buat grid untuk interpolasi (50 km x 50 km)
grid <- st_make_grid(shp_proj, cellsize = 50000)  # 50 km
grid_sf <- st_sf(geometry = grid)

# Konversi ke Spatial* object untuk analisis lebih lanjut
grid_sp <- as(grid_sf, "Spatial")

# OPSIONAL: Tambahkan bounding box dari shapefile asli
# untuk memastikan grid hanya di area yang relevan
# Fix the error with this minimal code
shp_proj <- st_make_valid(shp_proj)
grid_clipped <- st_make_grid(shp_proj, cellsize = 50000) %>%
  st_intersection(shp_proj)
bbox <- st_bbox(shp_proj)
grid_clipped <- st_make_grid(shp_proj, 
                             cellsize = 50000,
                             what = "polygons") %>%
  st_intersection(shp_proj)

4.3.3 Metode Invers Distance Weighting (IDW)

Peta Hasil Interpolasi IDW p=1

Peta interpolasi IDW (Inverse Distance Weighting) dengan parameter idp = 1 menggambarkan pola prediksi spasial jumlah sekolah di Indonesia dengan menekankan pengaruh kedekatan geografis antarprovinsi. Dalam metode ini, provinsi yang berlokasi lebih dekat memberikan kontribusi bobot yang lebih besar terhadap nilai prediksi suatu lokasi, sehingga pola spasial yang terbentuk sangat mencerminkan struktur ketetanggaan wilayah.

Hasil interpolasi menunjukkan bahwa wilayah dengan prediksi jumlah sekolah tinggi, yang ditunjukkan oleh warna kuning hingga merah tua, terkonsentrasi di Pulau Jawa serta sebagian wilayah Sumatra bagian selatan. Pola ini mencerminkan tingginya konsentrasi penduduk usia sekolah, tingkat urbanisasi, serta aktivitas ekonomi yang relatif lebih maju di wilayah tersebut. Temuan ini sejalan dengan hasil analisis deskriptif spasial dan uji autokorelasi sebelumnya, yang mengindikasikan adanya pengelompokan wilayah dengan karakteristik pendidikan yang serupa.

Sebaliknya, wilayah dengan prediksi jumlah sekolah rendah, yang ditunjukkan oleh warna lebih terang hingga kekuningan pucat, cenderung berada di Indonesia Timur, khususnya di wilayah Papua, Maluku, dan sebagian Nusa Tenggara. Pola ini mengindikasikan keterbatasan penyediaan fasilitas pendidikan yang kemungkinan dipengaruhi oleh faktor geografis, luas wilayah yang besar, kondisi topografi yang menantang, serta persebaran penduduk yang relatif jarang.

Gradasi warna yang relatif halus dan kontinu pada peta menunjukkan bahwa penggunaan idp = 1 menghasilkan permukaan interpolasi yang cukup smooth, sehingga transisi nilai prediksi antarwilayah terlihat gradual. Hal ini menandakan bahwa efek kedekatan spasial memiliki peran yang kuat dalam membentuk pola distribusi jumlah sekolah antarprovinsi.

Secara keseluruhan, peta interpolasi IDW ini memperkuat temuan dari uji autokorelasi spasial global dan lokal (Moran’s I dan LISA), bahwa distribusi jumlah sekolah antarprovinsi di Indonesia menunjukkan pola spasial yang terstruktur dan saling berkaitan. Oleh karena itu, kebijakan pemerataan pendidikan sebaiknya tidak dirancang secara parsial per provinsi, melainkan mempertimbangkan pendekatan spasial regional, khususnya untuk meningkatkan ketersediaan sekolah di wilayah dengan prediksi jumlah sekolah yang relatif rendah.

5 Kesimpulan dan Saran

5.1 Kesimpulan

Berdasarkan hasil analisis autokorelasi spasial, baik menggunakan Moran’s I, Geary’s C, maupun analisis LISA, dapat disimpulkan bahwa jumlah sekolah antarprovinsi di Indonesia menunjukkan adanya ketergantungan spasial yang signifikan. Provinsi-provinsi yang berdekatan secara geografis cenderung memiliki jumlah sekolah yang relatif mirip, yang ditunjukkan oleh pola klaster high–high dan low–low pada peta LISA. Temuan ini menegaskan bahwa distribusi fasilitas pendidikan tidak bersifat acak secara spasial.

Hasil pengujian lanjutan melalui regresi spasial menunjukkan bahwa model Spatial Error Model (SEM) merupakan model yang paling sesuai untuk menjelaskan hubungan antara jumlah sekolah dengan PDRB, luas wilayah, dan jumlah penduduk usia sekolah. Nilai parameter lambda yang signifikan mengindikasikan bahwa pengaruh spasial lebih dominan muncul melalui komponen galat, yang mencerminkan adanya faktor-faktor lain yang bersifat spasial dan tidak terobservasi dalam model. Secara parsial, PDRB berpengaruh negatif terhadap jumlah sekolah, sementara luas wilayah dan jumlah penduduk usia sekolah berpengaruh positif dan signifikan. Selain itu, nilai AIC SEM yang lebih kecil dibandingkan model OLS menunjukkan bahwa penggunaan model spasial mampu meningkatkan kualitas estimasi.

Pada tahap interpolasi spasial, metode Inverse Distance Weighting (IDW) mampu menggambarkan pola persebaran jumlah sekolah antarprovinsi secara kontinu berdasarkan kedekatan geografis. Hasil interpolasi menunjukkan bahwa wilayah dengan jumlah sekolah tinggi cenderung mengelompok dan memengaruhi wilayah di sekitarnya, sementara daerah dengan nilai rendah berada pada kawasan yang relatif terisolasi. Temuan ini memperkuat adanya keterkaitan spasial dalam distribusi jumlah sekolah dan mendukung penggunaan pendekatan spasial dalam analisis pemerataan pendidikan.

Secara keseluruhan, penelitian ini membuktikan bahwa analisis spasial memberikan pemahaman yang lebih komprehensif dalam mengkaji distribusi jumlah sekolah antarprovinsi dibandingkan pendekatan nonspasial semata.

5.2 Saran

Berdasarkan hasil penelitian yang telah diperoleh, beberapa saran yang dapat diajukan adalah sebagai berikut. 1. Bagi pembuat kebijakan, hasil penelitian ini dapat dijadikan dasar dalam perencanaan pembangunan pendidikan berbasis wilayah, dengan memperhatikan keterkaitan antarprovinsi. Intervensi kebijakan sebaiknya tidak dilakukan secara parsial, tetapi mempertimbangkan efek spasial dan kondisi wilayah sekitar, terutama pada provinsi yang tergolong low–low cluster. 2. Untuk penelitian selanjutnya, disarankan untuk menambahkan variabel penjelas lain yang berpotensi memengaruhi jumlah sekolah, seperti tingkat kemiskinan, belanja pendidikan daerah, atau aksesibilitas infrastruktur, guna mengurangi efek galat spasial yang masih tertangkap dalam SEM. 3. Pengembangan metode interpolasi dapat dilakukan dengan membandingkan lebih lanjut beberapa model variogram dan jenis kriging lain, seperti Universal Kriging atau Kriging with External Drift, agar diperoleh hasil prediksi yang lebih akurat. 4. Penggunaan data panel spasial antarwaktu direkomendasikan agar penelitian selanjutnya tidak hanya mampu menangkap pola spasial, tetapi juga dinamika temporal distribusi fasilitas pendidikan di Indonesia.

6 Lampiran

Dashboard: https://credo13.shinyapps.io/DashboardAnalisisSpasialUAS/

Video Presentasi: https://youtu.be/Cs_fDFBnhEM