Analisis Spasial Rata-rata Lama Sekolah di Provinsi Jawa Timur Tahun 2024
Disusun oleh: Nafalla Afftanur Rismawanti (140610230044)
Dosen Pengampu: Dr. I Gede Nyoman Mindra Jaya, M.Si
ABSTRAK
Rata-rata Lama Sekolah (RLS) merupakan salah satu indikator penting dalam mengukur capaian pembangunan pendidikan. Namun, distribusi RLS antarwilayah sering kali menunjukkan ketimpangan yang bersifat spasial, sehingga analisis yang mengabaikan ketergantungan antarwilayah berpotensi menghasilkan estimasi yang bias. Penelitian ini bertujuan untuk menganalisis pola distribusi spasial RLS, menguji keberadaan autokorelasi spasial, mengidentifikasi faktor-faktor sosial ekonomi dan pendidikan yang memengaruhi RLS, menentukan model spatial econometrics yang paling sesuai, serta membandingkan kinerja metode interpolasi spasial dalam memprediksi RLS secara kontinu di Provinsi Jawa Timur tahun 2024.
Data yang digunakan merupakan data cross-sectional tingkat kabupaten/kota yang bersumber dari Badan Pusat Statistik (BPS) dan Kementerian Keuangan, dengan unit analisis sebanyak 38 kabupaten/kota. Analisis autokorelasi spasial dilakukan menggunakan Moran’s I, Geary’s C, dan Getis–Ord G, serta analisis lokal menggunakan LISA dan Getis–Ord Gi*. Pemodelan dilakukan dengan pendekatan spatial econometrics, meliputi Ordinary Least Squares (OLS), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS), menggunakan matriks bobot spasial berbasis k-nearest neighbors (KNN) dengan k = 4. Selain itu, interpolasi spasial dilakukan menggunakan Universal Kriging (UK) dan Inverse Distance Weighting (IDW), dengan evaluasi akurasi melalui Leave-One-Out Cross-Validation.
Hasil penelitian menunjukkan bahwa Rata-Rata Lama Sekolah (RLS) di Provinsi Jawa Timur memiliki pola spasial yang signifikan dan tidak acak, ditandai dengan keberadaan wilayah hot spot dan cold spot. Variabel Tingkat Kemiskinan berpengaruh negatif dan signifikan terhadap RLS, sedangkan Angka Partisipasi Murni (APM) SMA berpengaruh positif dan signifikan. Berdasarkan perbandingan kriteria kecocokan model, General Nested Spatial Model (GNS) terpilih sebagai model terbaik dengan nilai Akaike Information Criterion (AIC) terendah sebesar 104,48, serta residual yang memenuhi asumsi normalitas, homoskedastisitas, dan bebas dari autokorelasi spasial. Pada analisis interpolasi, Universal Kriging menunjukkan kinerja yang lebih baik dibandingkan IDW, dengan nilai RMSE sebesar 1,6721 dan MAE sebesar 1,4457, serta mampu menyediakan informasi ketidakpastian prediksi melalui varians kriging.
Penelitian ini menegaskan pentingnya pendekatan spasial dalam analisis pendidikan serta memberikan implikasi kebijakan berbasis wilayah untuk mengurangi ketimpangan capaian pendidikan antar kabupaten/kota di Provinsi Jawa Timur.
Kata kunci: rata-rata lama sekolah, autokorelasi spasial, spatial econometrics, SDEM, Universal Kriging, Provinsi Jawa Timur.
BAB I PENDAHULUAN
1.1 Latar Belakang
Pendidikan merupakan salah satu pilar utama pembangunan manusia yang berperan strategis dalam meningkatkan kualitas sumber daya manusia, produktivitas ekonomi, serta kesejahteraan sosial masyarakat. Dalam kerangka pembangunan berkelanjutan, pendidikan tidak hanya dipandang sebagai tujuan akhir, tetapi juga sebagai instrumen penting untuk mengurangi ketimpangan sosial dan meningkatkan mobilitas ekonomi antargenerasi (UNDP, 2023). Oleh karena itu, pengukuran dan evaluasi capaian pendidikan menjadi aspek krusial dalam perumusan kebijakan pembangunan yang berbasis bukti.
Salah satu indikator pendidikan yang digunakan secara luas adalah Rata-rata Lama Sekolah (RLS), yaitu rata-rata jumlah tahun pendidikan formal yang telah ditempuh oleh penduduk usia 25 tahun ke atas. Indikator ini merupakan komponen utama dalam penyusunan Indeks Pembangunan Manusia (IPM) dan mencerminkan akumulasi investasi pendidikan jangka panjang di suatu wilayah (Badan Pusat Statistik, 2024). Nilai RLS yang tinggi mengindikasikan akses pendidikan yang lebih baik dan berkelanjutan, sedangkan nilai yang rendah dapat menjadi sinyal adanya keterbatasan struktural dalam sistem pendidikan maupun kondisi sosial ekonomi masyarakat.
Secara nasional, capaian RLS di Indonesia menunjukkan tren peningkatan dari tahun ke tahun. Namun demikian, peningkatan tersebut belum sepenuhnya merata antarprovinsi maupun antar kabupaten/kota. Ketimpangan spasial dalam capaian pendidikan masih menjadi permasalahan yang signifikan, terutama pada wilayah dengan karakteristik sosial ekonomi dan geografis yang heterogen (Badan Pusat Statistik, 2024). Faktor-faktor seperti tingkat kemiskinan, kapasitas fiskal daerah, ketersediaan tenaga pendidik, serta akses terhadap fasilitas pendidikan berkontribusi terhadap variasi capaian RLS antarwilayah.
Provinsi Jawa Timur merupakan salah satu provinsi dengan jumlah penduduk terbesar di Indonesia dan memiliki struktur wilayah yang sangat beragam, mulai dari kawasan metropolitan, wilayah industri, hingga daerah perdesaan dan kepulauan. Kondisi ini berpotensi menimbulkan disparitas capaian pendidikan antar kabupaten/kota. Data tahun 2024 menunjukkan bahwa nilai RLS di Jawa Timur berkisar antara 5,08 tahun hingga 12,11 tahun, dengan capaian terendah berada di Kabupaten Sampang dan capaian tertinggi di Kota Madiun. Perbedaan ini mencerminkan adanya ketimpangan pembangunan pendidikan antara wilayah perkotaan dan perdesaan, serta antara kawasan daratan utama dan wilayah kepulauan. Meskipun secara agregat Jawa Timur memiliki capaian IPM yang relatif baik, variasi internal antarwilayah mengindikasikan bahwa manfaat pembangunan pendidikan belum terdistribusi secara merata (BPS Provinsi Jawa Timur, 2024).
Selain faktor sosial ekonomi, capaian pendidikan juga dipengaruhi oleh interaksi spasial antarwilayah. Wilayah yang berdekatan secara geografis cenderung memiliki karakteristik pendidikan yang saling berkaitan akibat adanya efek limpahan (spatial spillover), seperti mobilitas penduduk, distribusi fasilitas pendidikan, dan keterkaitan kebijakan regional. Fenomena ini sejalan dengan First Law of Geography yang menyatakan bahwa segala sesuatu saling berhubungan, tetapi hal-hal yang berdekatan memiliki hubungan yang lebih kuat dibandingkan yang berjauhan (Tobler, 1970). Oleh karena itu, analisis yang mengabaikan dimensi spasial berpotensi menghasilkan estimasi yang bias dan kurang mencerminkan kondisi empiris yang sebenarnya.
Pendekatan analisis spasial menjadi penting untuk mengidentifikasi pola distribusi RLS, menguji keberadaan autokorelasi spasial, serta memodelkan hubungan antara RLS dan faktor-faktor penentunya dengan mempertimbangkan ketergantungan antarwilayah. Penggunaan metode spatial econometrics, seperti Spatial Lag Model (SLM), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM), memungkinkan peneliti untuk menangkap pengaruh langsung maupun tidak langsung antarwilayah yang tidak dapat dijelaskan oleh model regresi klasik (Anselin, 1988; LeSage & Pace, 2009).
Di sisi lain, data capaian RLS umumnya tersedia dalam bentuk diskret pada unit administrasi kabupaten/kota, sehingga menyulitkan analisis pada wilayah yang tidak memiliki observasi langsung. Untuk mengatasi keterbatasan tersebut, penelitian ini juga menerapkan metode interpolasi spasial, yaitu Universal Kriging (UK) dan Inverse Distance Weighting (IDW), guna memprediksi nilai RLS secara kontinu pada seluruh wilayah Provinsi Jawa Timur. Metode kriging secara khusus mampu mempertimbangkan struktur ketergantungan spasial melalui variogram serta menyediakan ukuran ketidakpastian prediksi, sehingga banyak direkomendasikan dalam analisis geostatistik (Cressie, 1993).
Berdasarkan latar belakang tersebut, penelitian ini bertujuan untuk menganalisis distribusi spasial Rata-rata Lama Sekolah di Provinsi Jawa Timur tahun 2024 dengan mengintegrasikan analisis autokorelasi spasial, pemodelan spatial econometrics, dan interpolasi spasial. Hasil penelitian ini diharapkan dapat memberikan pemahaman yang lebih komprehensif mengenai pola ketimpangan pendidikan antarwilayah serta menjadi dasar perumusan kebijakan pendidikan yang lebih tepat sasaran dan berbasis wilayah.
1.2 Rumusan Masalah
- Bagaimana distribusi spasial rata-rata lama sekolah di Jawa Timur tahun 2024?
- Apakah terdapat autokorelasi spasial pada distribusi RLS?
- Variabel apa saja yang signifikan mempengaruhi RLS di Jawa Timur?
- Model spatial econometrics mana yang paling sesuai?
- Bagaimana hasil interpolasi spasial (UK vs IDW) untuk prediksi RLS?
1.3 Tujuan Penelitian
Penelitian ini bertujuan untuk menganalisis pola dan faktor penentu rata-rata lama sekolah (RLS) di Provinsi Jawa Timur tahun 2024 dengan pendekatan analisis spasial. Secara khusus, tujuan penelitian ini adalah sebagai berikut:
Menganalisis dan memetakan distribusi spasial rata-rata lama sekolah pada tingkat kabupaten/kota di Provinsi Jawa Timur tahun 2024 guna mengidentifikasi variasi dan ketimpangan antarwilayah.
Menguji keberadaan dan tingkat autokorelasi spasial pada distribusi RLS menggunakan ukuran autokorelasi spasial global dan lokal, sehingga dapat diketahui apakah nilai RLS di suatu wilayah dipengaruhi oleh wilayah sekitarnya.
Mengidentifikasi variabel-variabel sosial ekonomi dan pendidikan yang berpengaruh signifikan terhadap RLS di Provinsi Jawa Timur dengan menggunakan model regresi sebagai pendekatan awal.
Menentukan model spatial econometrics yang paling sesuai dalam menjelaskan hubungan antara RLS dan variabel penentunya dengan mempertimbangkan ketergantungan spasial antarwilayah, berdasarkan uji diagnostik dan kriteria pemilihan model.
Membandingkan kinerja metode interpolasi spasial Universal Kriging (UK) dan Inverse Distance Weighting (IDW) dalam memprediksi nilai RLS secara kontinu di seluruh wilayah Provinsi Jawa Timur, serta mengevaluasi tingkat akurasi dan ketidakpastian prediksi yang dihasilkan.
1.4 Manfaat Penelitian
Penelitian ini diharapkan dapat memberikan manfaat sebagai berikut.
Manfaat Akademis
Menambah referensi empiris mengenai penerapan analisis spasial dalam kajian pendidikan, khususnya pada analisis rata-rata lama sekolah di tingkat kabupaten/kota.
Memberikan contoh penggunaan statistika spasial dan interpolasi spasial dalam mengidentifikasi pola distribusi dan ketimpangan capaian pendidikan antarwilayah.
Manfaat Praktis
Menyediakan informasi berbasis wilayah mengenai distribusi dan ketimpangan rata-rata lama sekolah di Provinsi Jawa Timur yang dapat digunakan sebagai dasar perencanaan pembangunan pendidikan.
Memberikan rekomendasi kebijakan pendidikan berbasis spasial bagi pemerintah daerah dalam menentukan wilayah prioritas intervensi pendidikan.
1.5 Keterbatasan Penelitian
- Penelitian ini memiliki beberapa keterbatasan yang perlu diperhatikan dalam interpretasi hasil. Pertama, data yang digunakan bersifat cross-sectional, sehingga analisis yang dilakukan hanya menggambarkan kondisi pada satu periode waktu dan belum mampu menangkap dinamika perubahan rata-rata lama sekolah antarwilayah dari waktu ke waktu. Oleh karena itu, hasil penelitian ini tidak dapat digunakan untuk menarik kesimpulan kausalitas jangka panjang. Kedua, unit analisis yang digunakan adalah kabupaten/kota, sehingga penelitian ini berpotensi mengandung ecological fallacy, yaitu ketidaksesuaian antara hubungan yang diamati pada tingkat agregat wilayah dengan kondisi pada tingkat individu. Interpretasi hasil penelitian ini dengan demikian terbatas pada konteks wilayah, bukan pada perilaku atau keputusan individu.
Ketiga, variabel independen yang digunakan dalam penelitian ini terbatas pada variabel sosial ekonomi dan pendidikan yang tersedia secara konsisten pada tingkat kabupaten/kota. Faktor-faktor lain yang berpotensi memengaruhi rata-rata lama sekolah, seperti kualitas fasilitas pendidikan, karakteristik rumah tangga, dan faktor budaya, belum dapat dimasukkan ke dalam model dan menjadi peluang pengembangan pada penelitian selanjutnya.
BAB II TINJAUAN PUSTAKA
2.1 Teori Spatial Dependence
Ketergantungan spasial merupakan fenomena di mana nilai suatu variabel pada suatu wilayah tidak sepenuhnya independen, melainkan dipengaruhi oleh nilai variabel yang terdapat pada wilayah tetangganya (Wiguna et al., 2022). Fenomena ini sangat penting dalam analisis data geospasial karena dapat memengaruhi hasil estimasi dan interpretasi model statistik. Secara umum, terdapat dua bentuk utama ketergantungan spasial, yaitu spatial lag dan spatial error (Novitasari et al., 2022).
Spatial lag muncul ketika nilai variabel dependen pada suatu wilayah secara langsung dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Dalam model regresi, efek ini direpresentasikan dengan menambahkan komponen WY ke persamaan regresi, sehingga persamaannya menjadi:
\[Y = \rho WY + X\beta + \varepsilon.\]
Di mana W adalah matriks bobot ruang yang menunjukkan kedekatan atau keterhubungan antarwilayah, mengukur kekuatan pengaruh lag, dan adalah kesalahan acak. Model spatial lag memungkinkan peneliti menangkap efek spillover, yaitu pengaruh lintas wilayah yang dapat terjadi akibat interaksi sosial, ekonomi, atau faktor pendidikan.
Sementara itu, spatial error mengacu pada autokorelasi yang terdapat pada komponen kesalahan model, bukan pada variabel dependen itu sendiri. Model ini menuliskan kesalahan sebagai:
\[\varepsilon = \lambda W \varepsilon + u,\]
dengan mengukur intensitas autokorelasi pada residual, u adalag error i.i.d., dan W kembali matriks bobot ruang (spatial weight matrix). Spatial error digunakan untuk menangkap penyebaran faktor-faktor tak terobservasi yang memiliki pola spasial. Jika faktor-faktor ini diabaikan, residual OLS akan menunjukkan pola yang terstruktur, dan model spatial error membantu memperbaiki bias tersebut.
Kedua bentuk ketergantungan ini sangat penting karena mengabaikannya dapat menghasilkan koefisien regresi yang bias, standar error yang terlalu kecil, dan keputusan kebijakan yang menyesatkan. Pemilihan antara spatial lag, spatial error, atau kombinasi keduanya bergantung pada apakah penyebaran terjadi pada variabel dependen secara langsung atau melalui faktor-faktor tak terukur yang tercermin pada residual. Jika spatial dependence diabaikan, estimasi koefisien OLS menjadi tidak efisien dan standar error dapat direduksi secara artifisial, sehingga uji signifikansi menjadi menyesatkan (Mar’ah et al., 2025). Model OLS yang mengabaikan autokorelasi spasial mungkin menampilkan koefisien X yang signifikan, padahal efek tersebut dipengaruhi oleh pola spasial yang tidak terkontrol. Sebaliknya, model spatial error mengoreksi korelasi residual dan memberikan estimasi yang lebih realistis (Eryando et al., 2022)
Dengan demikian, identifikasi dan penanganan ketergantungan spasial menjadi prasyarat penting sebelum menafsirkan hubungan antara variabel independen (X) dan rata-rata lama sekolah (Y).
2.2 Autokorelasi Spasial
Autokorelasi spasial merujuk pada korelasi yang muncul antara nilai suatu variabel dengan nilai variabel yang sama di lokasi-lokasi lain, yang disebabkan oleh kedekatan atau hubungan spasial antar pengamatan. Dengan kata lain, autokorelasi spasial menunjukkan sejauh mana nilai suatu variabel pada satu wilayah dipengaruhi oleh nilai variabel di wilayah sekitarnya. Secara matematis, autokorelasi spasial dapat dinyatakan melalui momen kovarians:
\[cov[y_i, y_j] = E[y_i y_j] – E[y_i] \cdot E[y_j] \neq 0\] di mana i dan j merujuk pada pengamatan individual (lokasi), dan yi adalah nilai variabel acak yang diamati pada lokasi i. Salah satu indeks yang paling umum digunakan untuk mengukur autokorelasi spasial global adalah Moran’s I, yang dirumuskan sebagai:
\[I = \frac{n \sum_{i} \sum_{j} w_{ij} (Y_i - \bar{Y})(Y_j - \bar{Y})}{(\sum_{i \neq j} w_{ij}) \sum_{i} (Y_i - \bar{Y})^2}\] dimana
n adalah jumlah wilayah
\(Y_{i}\) adalah nilai variabel yang diamati di wilayah i,
\(\overline{Y}\) adalah rata-rata dari nilai variabel, dan
\(w_{ij}\) adalah elemen matriks pembobot spasial (W) yang menunjukkan kedekatan spasial antara wilayah i dan j dengan \(w_{ii}\)= 0.
Matriks pembobot spasial (W) merupakan komponen penting dalam analisis autokorelasi, karena menentukan struktur ketetanggaan antarwilayah. Elemen \(w_{ij}\) bernilai 1 jika wilayah i dan j bertetangga, dan 0 jika tidak bertetangga. Matriks ini biasanya distandarisasi baris sehingga \[\sum{j} w{ij} = 1\] untuk setiap i, sehingga efek spatial lag dapat diinterpretasikan sebagai rata-rata tertimbang dari nilai-nilai tetangga. Nilai Moran’s I umumnya berada di kisaran -1 hingga +1. Nilai positif (I > 0) menandakan autokorelasi spasial positif, yaitu wilayah dengan nilai tinggi dikelilingi oleh wilayah dengan nilai tinggi, dan wilayah dengan nilai rendah dikelilingi oleh wilayah dengan nilai rendah. Nilai Moran’s I mendekati 0 (I0) mengindikasikan pola spasial yang acak (random) tanpa keteraturan spasial, di mana tidak ada pola sistematis dalam distribusi spasial variabel (tidak ada autokorelasi spasial). Sementara nilai negatif (I < 0) menunjukkan autokorelasi spasial negatif, di mana wilayah dengan nilai tinggi cenderung dikelilingi wilayah dengan nilai rendah, atau sebaliknya, yang dapat diartikan sebagai pola dispersi atau kompetisi spasial.
Autokorelasi spasial dapat terjadi pada variabel dependen maupun residual model regresi. Jika terdapat pada variabel dependen atau independen, hal ini menunjukkan adanya dependensi spasial substantif dalam data, sehingga model yang digunakan perlu mengakomodasi efek ini, misalnya melalui spatial lag model atau Spatial Durbin Model. Mengabaikannya dapat menyebabkan estimasi koefisien bias dan inkonsisten. Autokorelasi pada residual regresi merupakan indikator yang lebih kritis karena menandakan pelanggaran asumsi independensi error dalam model OLS. Dalam konteks regresi, kondisi ini dapat dinyatakan sebagai:
\[cov[\varepsilon_i, \varepsilon_j] \neq 0, untuk \ i \neq j\] Dalam praktiknya, pengujian autokorelasi spasial residual dilakukan sebagai bagian dari diagnostik spesifikasi model setelah estimasi OLS. Jika ditemukan autokorelasi yang signifikan, model perlu direspesifikasi, misalnya dengan menggunakan Spatial Error Model (SEM) atau model spasial lain yang lebih kompleks, agar estimasi menjadi konsisten dan efisien.
2.3 Teori Spatial Dependence
Dalam analisis ekonometrika spasial, model dasarnya adalah Ordinary Least Squares (OLS), yang memiliki asumsi bahwa setiap observasi bersifat independen satu sama lain. Pada konteks spasial, OLS memodelkan hubungan antara variabel dependen Y dan sejumlah variabel independen X1, X2, …, Xk tanpa mempertimbangkanposisi geografis masing-masing unit analisis. Bentuk umum model OLS adalah:
\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... + \beta_kX_k + \varepsilon\] Namun, asumsi independensi sering tidak terpenuhi dalam data wilayah. Fenomena ekonomi, sosial, maupun pendidikan cenderung menunjukkan pola non-acak di ruang, sehingga wilayah yang berdekatan sering memiliki karakteristik serupa akibat interaksi sosial, kebijakan publik, mobilitas tenaga kerja, dan kesamaan kondisi geografis. Apabila terdapat dependensi spasial (spatial dependence) baik pada variabel maupun pada komponen error maka penggunaan OLS menjadi tidak tepat. Dua sumber utama dependensi spasial adalah efek lag spasial dan efek error spasial. Efek lag spasial terjadi ketika nilai variabel dependen di suatu wilayah dipengaruhi secara langsung oleh nilai dependen di wilayah sekitarnya. Sebaliknya, efek error spasial muncul ketika kesalahan atau faktor tak teramati menunjukkan korelasi spasial, menandakan adanya variabel yang tidak terukur namun memengaruhi beberapa wilayah secara simultan. Ketika kedua jenis efek ini muncul secara bersamaan, estimasi OLS akan menghasilkan koefisien yang bias dan tidak efisien, serta residual yang menunjukkan autokorelasi spasial positif. Hal ini menandakan bahwa OLS gagal menangkap interaksi antarwilayah secara akurat. Oleh karena itu, diperlukan model spatial econometrics, yang secara eksplisit memperhitungkan ketergantungan spasial untuk memperoleh estimasi yang lebih tepat.
Beberapa model spasial yang umum digunakan antara lain: Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nesting Spatial Model (GNS). Model-model ini memungkinkan analisis yang lebih komprehensif, baik untuk efek langsung antarwilayah maupun efek spillover dan distribusi error yang bersifat spasial. Berikut adalah penjelasan untuk masing-masing modelnya.
- Spatial Autoregressive (SAR) Model SAR memasukkan lag spasial dari variabel dependen ke dalam persamaan regresi. Dengan kata lain, nilai dependen di suatu wilayah dipengaruhi secara langsung oleh nilai dependen di wilayah tetangganya. Persamaan umumnya dapat dituliskan sebagai:
\[Y = \rho WY + X\beta + \varepsilon\]
Keterangan:
Y : Variabel dependen, yaitu ata-rata Lama Sekolah (RLS) di masing-masing kabupaten/kota.
W : Matriks bobot spasial (spatial weight matrix), yang menunjukkan hubungan kedekatan antarwilayah (misalnya berbatasan langsung, atau berjarak dalam radius tertentu).
WY: Lag spasial dari variabel dependen, yaitu rata-rata tertimbang RLS di wilayah-wilayah tetangga.
\(\rho\): Koefisien lag spasial, mengukur seberapa besar pengaruh rata-rata lama sekolah di daerah tetangga terhadap daerah bersangkutan.
X: Matriks variabel independen.
\(\beta\):Vektor koefisien regresi, menunjukkan besarnya pengaruh masing-masing variabel X terhadap RLS.
\(\varepsilon\): Komponen error acak, diasumsikan identik dan independen.
Model SAR cocok digunakan ketika terdapat pengaruh langsung antarwilayah, misalnya daerah dengan RLS tinggi memengaruhi wilayah tetangganya melalui pertukaran sumber daya pendidikan atau kebijakan lokal.
- Spatial Error Model (SEM) Model SEM digunakan ketika pengaruh spasial muncul dari variabel yang tidak teramati tetapi memiliki pola spasial. Fokus model ini adalah menangkap autokorelasi pada komponen error, bukan variabel dependen. Persamaan model SEM adalah:
\[Y = X\beta + u, u = \lambda Wu + \varepsilon\] Keterangan:
u : Error yang memiliki dependensi spasial.
λ : Koefisien autokorelasi spasial pada error, menunjukkan seberapa kuat keterkaitan error antarwilayah.
Wu : Lag spasial dari error, yaitu rata-rata kesalahan prediksi di wilayah tetangga.
ε : Komponen error acak yang bersifat independen.
SEM relevan ketika faktor tak teramati yang memengaruhi RLS tersebar secara spasial, sehingga residual OLS menunjukkan pola yang terstruktur.
- Spatial Durbin Model (SDM)
SDM merupakan pengembangan dari SAR karena tidak hanya memasukkan lag spasial pada variabel dependen, tetapi juga pada variabel independen. Model ini mampu menangkap efek langsung dan efek tidak langsung (spillover) antarwilayah. Persamaannya:
\[Y = \rho WY + WX\theta + X\beta + \varepsilon\]
Keterangan: θ : Koefisien efek tidak langsung (spillover) dari variabel independen di wilayah sekitar.
Contohnya, peningkatan PAD atau DAU Pendidikan di satu daerah dapat berpengaruh pada RLS daerah tersebut sekaligus wilayah sekitarnya melalui pergerakan tenaga pengajar atau fasilitas pendidikan.
Spatial Durbin Error Model (SDEM)
SDEM memperluas SEM dengan menambahkan lag spasial pada variabel independen. Model ini sesuai ketika pengaruh antarwilayah berasal dari faktor penjelas (X) dan error secara bersamaan, tanpa ketergantungan langsung antarvariabel dependen. Persamaannya:
\[Y = X\beta + WX\Theta + u, u = \lambda Wu + \varepsilon\] 5. Spatial Autoregressive Combined Model (SAC)
Model SAC menggabungkan efek lag spasial (SAR) dan error spasial (SEM). Model ini digunakan ketika terdapat interaksi langsung antarwilayah sekaligus variabel tak teramati yang berpola spasial. Persamaannya:
\[Y = \rho WY + X\beta + u, u = \lambda Wu + \varepsilon\]
6.General Nesting Spatial (GNS)
Model GNS merupakan model paling umum karena menggabungkan semua kemungkinan dependensi spasial. Model ini bersifat nested terhadap model lain (SAR, SEM, SDM, SDEM, SAC) sehingga dapat disederhanakan jika beberapa parameter tidak signifikan. Persamaannya:
\[Y = \rho WY + WX\theta + X\beta + u, u = \lambda Wu + \varepsilon\]
GNS memungkinkan estimasi efek langsung, tidak langsung, dan total dari variabel independen terhadap RLS antarwilayah, sehingga memberikan gambaran paling komprehensif mengenai interaksi spasial di wilayah penelitian.
2.4 Interpolasi Spasial
Interpolasi spasial merupakan metode statistik yang digunakan untuk mengestimasi nilai suatu variabel pada lokasi yang tidak terobservasi berdasarkan nilai yang diketahui pada lokasi sekitarnya.
Pemilihan metode interpolasi bergantung pada karakteristik data, keberadaan tren spasial, serta asumsi mengenai struktur ketergantungan spasial. Dalam penelitian ini, dua metode interpolasi yang digunakan adalah Universal Kriging (UK) dan Inverse Distance Weighting (IDW).
2.4.1 Universal Kriging (UK)
Universal Kriging merupakan salah satu teknik interpolasi geostatistik yang tergolong dalam metode kriging. Berbeda dengan Ordinary Kriging yang mengasumsikan rata-rata lokal konstan (stasioneritas orde kedua), Universal Kriging memperbolehkan adanya tren deterministik dalam data, sehingga cocok digunakan ketika terdapat variasi sistematis berdasarkan posisi geografis (Cressie, 1993).
Model Universal Kriging dapat dinyatakan sebagai:
\[ Z(s) = \mu(s) + \delta(s), \quad s \in D \]
Variogram Salah satu komponen utama dalam metode kriging adalah variogram (atau semivariogram), yang digunakan untuk mengukur ketergantungan spasial antara nilai-nilai yang diamati pada berbagai jarak. Variogram empiris didefinisikan sebagai:
\[ \gamma(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} \left[ Z(s_i) - Z(s_i + h) \right]^2 \]
Variogram empiris kemudian dimodelkan menggunakan fungsi teoritis, seperti model Spherical, Exponential, atau Gaussian. Parameter utama dari model variogram adalah:
Nugget: Variasi pada jarak nol, yang mencerminkan kesalahan pengukuran atau variasi skala mikro.
Sill: Nilai variogram pada saat mencapai titik jenuh, menunjukkan total variabilitas data.
Range: Jarak di mana variogram mencapai sill, menandakan batas pengaruh spasial.
Universal Kriging lebih sesuai digunakan ketika:
Terdapat tren spasial yang signifikan dalam data (misalnya, nilai RLS cenderung meningkat atau menurun berdasarkan posisi geografis).
Asumsi rata-rata konstan tidak terpenuhi.
Dibutuhkan estimasi yang memperhitungkan struktur ketergantungan spasial melalui variogram.
Diperlukan ukuran ketidakpastian prediksi melalui kriging variance.
Dalam penelitian ini, hasil uji tren spasial menunjukkan adanya pengaruh signifikan dari koordinat X (barat–timur), sehingga Universal Kriging dipilih untuk mengakomodasi tren tersebut.
2.4.2 Inverse Distance Weighting (IDW)
Inverse Distance Weighting (IDW) adalah metode interpolasi deterministik yang mengestimasi nilai pada lokasi yang tidak terobservasi berdasarkan rata-rata tertimbang dari nilai-nilai di lokasi sekitarnya, dengan bobot yang berbanding terbalik terhadap jarak (Shepard, 1968).
\[ \hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i \, Z(s_i)}{\sum_{i=1}^{n} w_i} \]
dengan:
\(\hat{Z}(s_0)\) adalah nilai prediksi pada lokasi \(s_0\)
\(Z(s_i)\) adalah nilai observasi pada lokasi \(s_i\),
\(w_i\) adalah bobot pengamatan ke-\(i\) yang ditentukan berdasarkan jarak,
\(d_{i0}\) adalah jarak antara lokasi prediksi \(s_0\) dan lokasi pengamatan \(s_i\),
\(p\) adalah parameter power yang mengontrol tingkat pengaruh jarak terhadap bobot.
Parameter Power
Parameter power (p) menentukan seberapa cepat pengaruh titik pengamatan menurun terhadap jarak. Nilai p yang umum digunakan adalah:
Nilai parameter \(p\) pada metode Inverse Distance Weighting (IDW) menentukan seberapa cepat pengaruh suatu titik pengamatan menurun seiring bertambahnya jarak. Interpretasi nilai \(p\) adalah sebagai berikut:
\(p = 1\): Pengaruh menurun secara linier terhadap jarak, menghasilkan permukaan interpolasi yang lebih halus.
\(p = 2\): Pengaruh menurun secara kuadratik terhadap jarak, merupakan nilai standar yang paling umum digunakan.
\(p > 2\): Pengaruh menurun lebih cepat, sehingga interpolasi menjadi lebih bersifat lokal dan lebih sensitif terhadap variasi spasial di sekitar titik pengamatan.
Dalam penelitian ini, nilai p=2 digunakan sebagai standar untuk memberikan keseimbangan antara pengaruh lokal dan global.
Kelebihan dan Kekurangan IDW
Kelebihan:
Sederhana dan mudah dipahami, tidak memerlukan asumsi distribusi data. Cepat secara komputasi, cocok untuk dataset besar. Tidak memerlukan model variogram atau estimasi parameter struktural.
Kekurangan:
Tidak memperhitungkan struktur ketergantungan spasial secara eksplisit. Rentan terhadap efek “bull’s eye”, yaitu hasil interpolasi yang terlalu ekstrem di sekitar titik pengamatan. Tidak menyediakan ukuran ketidakpastian prediksi. Sensitif terhadap nilai ekstrem (outliers) dan pemilihan parameter power.
Meskipun memiliki keterbatasan, IDW tetap digunakan sebagai metode pembanding karena kesederhanaannya dan efektivitasnya dalam interpolasi berbasis jarak.
2.4.3 Perbandingan UK vs IDW
Berdasarkan karakteristik data dan hasil analisis, metode interpolasi yang paling sesuai dipilih dengan mempertimbangkan hal-hal berikut:
Keberadaan Tren Spasial Jika terdapat tren spasial yang signifikan (misalnya, pengaruh koordinat geografis terhadap RLS), maka Universal Kriging lebih sesuai karena mampu mengakomodasi tren tersebut. Sebaliknya, jika data tidak memiliki tren sistematis, Ordinary Kriging atau bahkan IDW dapat dipertimbangkan. Struktur Ketergantungan Spasial Jika variogram menunjukkan struktur ketergantungan spasial yang kuat dan dapat dimodelkan dengan baik, maka kriging lebih direkomendasikan. IDW lebih sesuai jika ketergantungan spasial lemah atau tidak jelas. Kebutuhan Ukuran Ketidakpastian Jika diperlukan informasi mengenai ketidakpastian prediksi (misalnya untuk analisis risiko atau perencanaan kebijakan), maka kriging menjadi pilihan utama karena menyediakan kriging variance. IDW tidak menyediakan ukuran ketidakpastian. Validasi Empiris Kinerja kedua metode dibandingkan menggunakan teknik Leave-One-Out Cross-Validation (LOOCV), dengan menghitung metrik kesalahan prediksi seperti:
Kinerja metode interpolasi dievaluasi menggunakan beberapa ukuran kesalahan, yaitu Root Mean Square Error (RMSE), Mean Absolute Error (MAE), dan Koefisien Determinasi (\(R^2\)), yang dirumuskan sebagai berikut.
- Root Mean Square Error (RMSE)
\[ \mathrm{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( Z_i - \hat{Z}_i \right)^2} \]
Mean Absolute Error (MAE) \[ \mathrm{MAE} = \frac{1}{n} \sum_{i=1}^{n} \left| Z_i - \hat{Z}_i \right| \]
Koefisien Determinasi \[ R^2 = 1 - \frac{\sum_{i=1}^{n} \left( Z_i - \hat{Z}_i \right)^2} {\sum_{i=1}^{n} \left( Z_i - \bar{Z} \right)^2} \]
Metode interpolasi terbaik ditentukan berdasarkan nilai RMSE dan MAE yang paling kecil serta nilai \(R^2\) yang paling besar, yang menunjukkan tingkat akurasi dan kemampuan model dalam menjelaskan variasi data.
BAB III METODOLOGI PENELITIAN
3.1 Data
Penelitian ini menggunakan data sekunder yang bersumber dari instansi resmi pemerintah, khususnya Badan Pusat Statistik (BPS) dan Kementerian Keuangan Republik Indonesia, dengan unit analisis berupa 38 kabupaten/kota di Provinsi Jawa Timur. Data yang digunakan bersifat cross-sectional untuk tahun 2024, sehingga merepresentasikan kondisi pendidikan dan sosial ekonomi pada satu periode waktu tertentu.
3.1.1 Variabel Penelitian
Variabel yang digunakan dalam penelitian ini terdiri dari satu variabel dependen (Y) dan empat variabel independen (X), yang dipilih berdasarkan relevansinya terhadap capaian pendidikan serta ketersediaan data pada tingkat kabupaten/kota.
Variabel Dependen (Y)
Rata-rata Lama Sekolah (RLS) Rata-rata Lama Sekolah merupakan indikator pendidikan yang menggambarkan rata-rata jumlah tahun pendidikan formal yang telah ditempuh oleh penduduk berusia 15 tahun ke atas. Variabel ini diukur dalam satuan tahun dan digunakan sebagai indikator utama untuk menilai tingkat capaian pendidikan penduduk di suatu wilayah (Badan Pusat Statistik, 2024).
RLS dipilih sebagai variabel dependen karena merupakan salah satu komponen utama dalam penyusunan Indeks Pembangunan Manusia (IPM) dan mencerminkan akumulasi investasi pendidikan jangka panjang di suatu daerah (UNDP, 2023). Variasi nilai RLS antarwilayah juga sering digunakan untuk mengidentifikasi ketimpangan akses dan kualitas pendidikan secara spasial.
Variabel Independen (X)
Pendapatan Asli Daerah (PAD)
Pendapatan Asli Daerah merupakan total penerimaan daerah yang berasal dari pajak daerah, retribusi daerah, hasil pengelolaan kekayaan daerah yang dipisahkan, serta sumber pendapatan sah lainnya. PAD diukur dalam satuan miliar rupiah. Secara teoritis, PAD mencerminkan kapasitas fiskal daerah dalam membiayai pembangunan, termasuk penyediaan infrastruktur dan layanan pendidikan. Daerah dengan PAD yang lebih tinggi memiliki peluang lebih besar untuk meningkatkan akses dan kualitas pendidikan, sehingga diharapkan berdampak positif terhadap RLS (Farizki, 2022).
Tingkat Kemiskinan Tingkat kemiskinan menunjukkan persentase penduduk yang memiliki pengeluaran di bawah garis kemiskinan yang ditetapkan oleh Badan Pusat Statistik. Variabel ini diukur dalam satuan persentase (%). Tingkat kemiskinan yang tinggi mencerminkan keterbatasan kemampuan ekonomi rumah tangga dalam membiayai pendidikan, yang berpotensi menurunkan keberlanjutan pendidikan dan berdampak negatif terhadap RLS (Efficient Journal, 2020).
Rasio Guru–Murid SMA
Rasio guru–murid SMA merupakan perbandingan antara jumlah murid dan jumlah guru pada jenjang Sekolah Menengah Atas (SMA) di setiap kabupaten/kota. Variabel ini diukur dalam bentuk rasio murid per guru. Rasio yang tinggi menunjukkan beban pengajaran guru yang lebih besar, yang berpotensi menurunkan kualitas pembelajaran dan berdampak negatif terhadap capaian pendidikan, termasuk rata-rata lama sekolah (Stochastic Frontier Analysis, 2013).
Angka Partisipasi Murni (APM) SMA Angka Partisipasi Murni SMA menunjukkan proporsi penduduk usia sekolah menengah atas yang sedang bersekolah pada jenjang SMA sesuai dengan kelompok usia resminya. Variabel ini diukur dalam satuan persentase (%). Tingginya APM SMA mencerminkan keberlanjutan pendidikan pada usia sekolah menengah, yang secara langsung berkontribusi terhadap peningkatan rata-rata lama sekolah di suatu wilayah (Badan Pusat Statistik, 2024).
Tabel Daftar Variabel
| Jenis Variabel | Nama Variabel | Keterangan | Satuan |
|---|---|---|---|
| Dependen | Rata-rata Lama Sekolah (RLS) | Rata-rata jumlah tahun pendidikan penduduk usia ≥15 tahun | Tahun |
| Independen | Pendapatan Asli Daerah (PAD) | Penerimaan asli pemerintah daerah | Miliar Rupiah |
| Independen | Tingkat Kemiskinan | Persentase penduduk miskin | Persen (%) |
| Independen | Rasio Guru–Murid SMA | Jumlah murid per guru SMA | Rasio |
| Independen | Angka Partisipasi Murni (APM) SMA | Persentase penduduk usia SMA yang bersekolah | Persen (%) |
3.1.2 Sumber Data
Sumber data dalam penelitian ini dirangkum pada Tabel berikut.
| Variabel | Sumber Data |
|---|---|
| Rata-rata Lama Sekolah (RLS) | Badan Pusat Statistik Provinsi Jawa Timur – Tabel Statistik RLS Kabupaten/Kota 2024 |
| Pendapatan Asli Daerah (PAD) | Direktorat Jenderal Perimbangan Keuangan (DJPK), Kementerian Keuangan RI – Data APBD Tahun 2024 |
| Tingkat Kemiskinan | Badan Pusat Statistik Provinsi Jawa Timur – Jumlah dan Persentase Penduduk Miskin Kabupaten/Kota 2024 |
| Jumlah Guru dan Murid SMA | Badan Pusat Statistik Provinsi Jawa Timur – Statistik Pendidikan SMA |
| Angka Partisipasi Murni (APM) SMA | Badan Pusat Statistik Provinsi Jawa Timur – Statistik Pendidikan Provinsi Jawa Timur 2024 |
Seluruh data dikumpulkan dari situs resmi lembaga terkait, yaitu:
BPS Provinsi Jawa Timur (https://jatim.bps.go.id )
Direktorat Jenderal Perimbangan Keuangan, Kementerian Keuangan RI (https://djpk.kemenkeu.go.id )
3.1.3 Karakteristik Data
Data yang digunakan bersifat:
Cross-sectional, sehingga tidak menangkap dinamika perubahan antarwaktu,
Spasial, karena setiap observasi memiliki keterkaitan geografis antarwilayah,
Agregat wilayah, dengan unit kabupaten/kota, sehingga hasil analisis tidak merepresentasikan perilaku individu (ecological data).
Karakteristik ini menjadi dasar dalam pemilihan metode spatial econometrics dan interpolasi spasial pada tahapan analisis selanjutnya.
3.2 Unit Spasial
Unit spasial dalam penelitian ini adalah kabupaten/kota di Provinsi Jawa Timur, yang berjumlah 38 wilayah administratif, terdiri atas 29 kabupaten dan 9 kota. Setiap kabupaten/kota diperlakukan sebagai satu unit observasi yang memiliki karakteristik sosial, ekonomi, dan pendidikan yang berbeda, serta berpotensi saling memengaruhi secara spasial.
Representasi spasial wilayah penelitian menggunakan data batas administrasi tingkat kabupaten/kota dalam bentuk shapefile, yang memuat informasi geometris berupa poligon wilayah. Data spasial ini digunakan dalam penyusunan peta tematik, analisis autokorelasi spasial, pemodelan spatial econometrics, serta interpolasi spasial.
Hubungan spasial antarwilayah ditentukan menggunakan matriks bobot spasial (spatial weight matrix) berbasis k-nearest neighbors (KNN) dengan k = 4. Dalam pendekatan ini, setiap kabupaten/kota diasumsikan memiliki hubungan spasial dengan empat wilayah terdekat berdasarkan jarak geografis antar pusat wilayah (centroid). Pendekatan KNN memastikan bahwa setiap unit spasial memiliki jumlah tetangga yang sama, sehingga matriks bobot menjadi lebih seimbang dan terhindar dari masalah wilayah terisolasi (isolated regions).
Matriks bobot spasial kemudian dilakukan standardisasi baris (row-standardized), sehingga jumlah bobot pada setiap baris bernilai satu. Dengan standardisasi ini, pengaruh spasial dapat diinterpretasikan sebagai rata-rata tertimbang dari empat wilayah terdekat, yang mencerminkan intensitas keterkaitan spasial lokal antar kabupaten/kota.
Pemilihan pendekatan KNN (k = 4) dinilai sesuai dengan karakteristik wilayah Provinsi Jawa Timur yang memiliki bentuk administratif dan ukuran wilayah yang beragam. Pendekatan ini lebih fleksibel dibandingkan metode kontiguitas karena tidak bergantung pada batas administratif semata, tetapi pada kedekatan geografis aktual, sehingga lebih mampu menangkap interaksi spasial pendidikan dan sosial ekonomi antarwilayah. Oleh karena itu, unit spasial dan matriks bobot yang digunakan diharapkan mampu merepresentasikan pola ketergantungan spasial Rata-rata Lama Sekolah (RLS) secara lebih akurat.
3.3 Metode Analisis
Analisis dalam penelitian ini dilakukan melalui beberapa tahapan yang sistematis untuk mengidentifikasi pola spasial, menguji keberadaan ketergantungan spasial, serta menentukan model spasial dan metode interpolasi yang paling sesuai dalam menjelaskan variasi Rata-rata Lama Sekolah (RLS) antarwilayah di Provinsi Jawa Timur. Tahapan metode analisis yang digunakan dijelaskan sebagai berikut.
- Analisis Eksploratif dan Visualisasi Spasial
Tahap awal penelitian dilakukan dengan analisis eksploratif untuk memahami karakteristik data dan distribusi spasial variabel yang diteliti. Analisis ini meliputi statistik deskriptif untuk seluruh variabel penelitian serta visualisasi spasial dalam bentuk peta tematik (choropleth map) pada tingkat kabupaten/kota di Provinsi Jawa Timur. Visualisasi ini bertujuan untuk mengidentifikasi pola awal, variasi antarwilayah, serta indikasi ketimpangan capaian RLS secara geografis.
- Analisis Autokorelasi Spasial Global dan Lokal
Setelah tahap eksplorasi awal, dilakukan pengujian autokorelasi spasial global tidak hanya pada variabel dependen Rata-rata Lama Sekolah (RLS), tetapi juga pada seluruh variabel independen yang digunakan dalam penelitian. Pengujian ini dilakukan menggunakan beberapa ukuran autokorelasi spasial global, yaitu Moran’s I, Geary’s C, dan Getis–Ord General G, dengan tujuan untuk mengetahui apakah masing-masing variabel tersebar secara acak atau menunjukkan pola pengelompokan spasial yang signifikan antarwilayah.
Pengujian autokorelasi spasial pada variabel independen penting dilakukan untuk mengidentifikasi potensi ketergantungan spasial pada faktor-faktor sosial ekonomi dan pendidikan yang memengaruhi RLS. Keberadaan autokorelasi spasial pada variabel independen mengindikasikan adanya pola spasial yang dapat berkontribusi terhadap munculnya efek limpahan (spillover effects) dalam pemodelan regresi spasial.
- Estimasi Model Regresi OLS sebagai Model Dasar
Sebelum menerapkan model spasial, dilakukan estimasi regresi linier berganda (Ordinary Least Squares/OLS) sebagai model dasar untuk menganalisis pengaruh variabel independen terhadap RLS. Model OLS digunakan sebagai pembanding awal dan untuk menguji apakah hubungan antarvariabel dapat dijelaskan tanpa mempertimbangkan ketergantungan spasial.
Residual dari model OLS kemudian diuji menggunakan Moran’s I untuk mengetahui apakah masih terdapat autokorelasi spasial yang signifikan pada residual. Jika residual menunjukkan ketergantungan spasial, maka model OLS dianggap tidak memadai dan perlu dikembangkan ke model regresi spasial.
- Uji Diagnostik Lagrange Multiplier (LM)
Untuk menentukan bentuk ketergantungan spasial yang paling sesuai, dilakukan uji Lagrange Multiplier (LM) dan Robust LM. Uji ini digunakan untuk mengidentifikasi apakah ketergantungan spasial lebih dominan terjadi pada variabel dependen (spatial lag) atau pada komponen error (spatial error).
Apabila LM-Lag signifikan, maka model dengan ketergantungan spasial pada variabel dependen lebih sesuai. Sebaliknya, jika LM-Error signifikan, maka model dengan ketergantungan spasial pada error lebih tepat. Jika kedua uji signifikan, maka diperlukan penggunaan model spasial lanjutan yang mampu menangkap efek spasial secara lebih komprehensif.
- Estimasi Model Spatial Econometrics
Berdasarkan hasil uji LM, penelitian ini mengestimasi beberapa model regresi spasial lanjutan untuk menangkap berbagai bentuk ketergantungan spasial, yaitu:
Spatial Durbin Model (SDM)
Spatial Durbin Error Model (SDEM)
Spatial Autoregressive Combined Model (SAC)
General Nested Spatial Model (GNS)
Model-model tersebut dipilih karena mampu menangkap pengaruh spasial baik melalui variabel dependen, variabel independen, maupun komponen error, serta memungkinkan analisis efek langsung dan tidak langsung (spillover effects) antarwilayah.
- Pemilihan dan Evaluasi Model Terbaik
Kinerja masing-masing model spasial dievaluasi menggunakan beberapa kriteria, yaitu Akaike Information Criterion (AIC), log-likelihood, serta signifikansi parameter spasial. Model dengan nilai AIC terendah dan nilai log-likelihood tertinggi dipilih sebagai model terbaik. Selain itu, dilakukan analisis Delta AIC, Akaike Weights, dan Likelihood Ratio Test untuk memperkuat keputusan pemilihan model.
- Interpretasi Efek Spasial
Setelah model terbaik ditentukan, dilakukan interpretasi terhadap hasil estimasi dengan menitikberatkan pada efek langsung, efek tidak langsung (spillover), dan efek total dari variabel independen terhadap RLS. Interpretasi ini bertujuan untuk memahami sejauh mana faktor sosial ekonomi dan pendidikan di suatu wilayah tidak hanya memengaruhi wilayah tersebut, tetapi juga wilayah sekitarnya di Provinsi Jawa Timur.
- Analisis Interpolasi Spasial
Selain pemodelan regresi spasial, penelitian ini juga melakukan analisis interpolasi spasial untuk memprediksi nilai RLS secara kontinu di seluruh wilayah Provinsi Jawa Timur. Dua metode interpolasi yang digunakan adalah Universal Kriging (UK) dan Inverse Distance Weighting (IDW).
Metode Universal Kriging diterapkan dengan mempertimbangkan adanya tren spasial melalui koordinat geografis, serta struktur ketergantungan spasial yang direpresentasikan oleh variogram. Sementara itu, metode IDW digunakan sebagai metode pembanding berbasis jarak. Kinerja kedua metode interpolasi dievaluasi menggunakan cross-validation, dengan membandingkan nilai Root Mean Square Error (RMSE), Mean Absolute Error (MAE), dan koefisien determinasi (R²). Metode dengan tingkat kesalahan prediksi yang lebih rendah dianggap memberikan hasil interpolasi yang lebih akurat.
3.5 Alur Penelitian
Penelitian ini dilaksanakan melalui beberapa tahapan analisis yang sistematis untuk memperoleh hasil yang komprehensif. Alur penelitian disusun agar setiap tahap saling berkaitan dan menghasilkan model spasial terbaik untuk menjelaskan variasi Rata-rata Lama Sekolah (RLS) di Provinsi Jawa Barat. Alur penelitian dituangkan kedalam diagram di bawah ini.
## Warning: package 'DiagrammeR' was built under R version 4.5.1grViz("
digraph flowchart {
rankdir=TB
Start [shape=ellipse, label='Mulai']
Data [label='Input Data']
Preprocess [label='Pra-pemrosesan Data\\n']
EDA [label='Analisis Eksploratif\\n& Peta Tematik']
GlobalAuto [label='Autokorelasi Spasial Global\\n(Moran I, Geary C, Getis–Ord G)']
LocalAuto [label='Autokorelasi Spasial Lokal\\n(LISA & Getis–Ord Gi*)']
OLS [label='Estimasi Model Dasar\\n(OLS)']
MoranRes [shape=diamond, label='Autokorelasi Spasial\\npada Residual OLS?']
LM [label='Uji Lagrange Multiplier\\n(LM & Robust LM)']
Spasial [label='Model Spatial Econometrics\\n(SDM, SDEM, SAC, GNS)']
Select [label='Pemilihan Model Terbaik\\n(AIC, LogLik, LR Test)']
Impact [label='Interpretasi Efek Spasial\\n(Direct, Indirect, Total)']
FinishOLS [shape=ellipse, label='Selesai\\n(OLS Memadai)']
Finish [shape=ellipse, label='Kesimpulan & Rekomendasi']
Start -> Data -> Preprocess -> EDA -> GlobalAuto -> LocalAuto -> OLS -> MoranRes
MoranRes -> FinishOLS [label='Tidak']
MoranRes -> LM [label='Ya']
LM -> Spasial -> Select -> Impact -> Finish
}
")BAB IV HASIL DAN PEMBAHASAN
4.1 Statistika Deskriptif
Analisis deskriptif dilakukan untuk memberikan gambaran awal mengenai karakteristik data serta variasi spasial Rata-rata Lama Sekolah (RLS) dan variabel-variabel ekonomi serta pendidikan yang digunakan dalam penelitian ini. Analisis ini mencakup penyajian statistik ringkas dan visualisasi spasial dalam bentuk peta tematik pada tingkat kabupaten/kota di Provinsi Jawa Timur tahun 2024. Tahap ini bertujuan untuk mengidentifikasi pola awal, ketimpangan antarwilayah, serta indikasi awal adanya keterkaitan spasial sebelum dilakukan analisis autokorelasi dan pemodelan spasial lanjutan.
4.1.1 Peta Distribusi Rata-Rata Lama Sekolah
## Warning: package 'spdep' was built under R version 4.5.1## Loading required package: spData## Warning: package 'spData' was built under R version 4.5.1## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Warning: package 'sf' was built under R version 4.5.1## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE## Warning: package 'sp' was built under R version 4.5.1## Warning: package 'ggplot2' was built under R version 4.5.1## Warning: package 'dplyr' was built under R version 4.5.1##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'openxlsx' was built under R version 4.5.1## Warning: package 'viridis' was built under R version 4.5.1## Loading required package: viridisLite## Warning: package 'spatialreg' was built under R version 4.5.1## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Warning: package 'car' was built under R version 4.5.1## Loading required package: carData## Warning: package 'carData' was built under R version 4.5.1##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeIndo2 <- readRDS("E:/SMT 5/Analisis Data Spasial/UAS/gadm36_IDN_2_sp.rds")
# Filter Jawa Timur & Konversi ke sf
jatim_sf <- st_as_sf(Indo2[Indo2$NAME_1 == "Jawa Timur", ])
jatim_sf <- st_zm(jatim_sf, drop = TRUE, what = "ZM")
# Load Data Excel Jatim
data <- read.xlsx("E:/SMT 5/Analisis Data Spasial/UAS/DATA UAS SPASIAL.xlsx")
# Penyeragaman Nama Daerah agar Join Berhasil
jatim_sf_clean <- jatim_sf %>%
mutate(NAME_FIX = case_when(
NAME_2 == "Surabaya" ~ "Kota Surabaya",
NAME_2 == "Batu" ~ "Kota Batu",
TYPE_2 == "Kota" & !grepl("Kota", NAME_2) ~ paste("Kota", NAME_2),
TRUE ~ NAME_2
))
# Join Data Excel dengan Shapefile
jatim_merged <- jatim_sf_clean %>%
left_join(data, by = c("NAME_FIX" = "Kabupaten/Kota"))
library(sf)
library(spdep)
library(ggplot2)
library(dplyr)
library(viridis)
cat("\n=== PETA SPASIAL RATA-RATA LAMA SEKOLAH ===\n")##
## === PETA SPASIAL RATA-RATA LAMA SEKOLAH ===# Peta Choropleth dengan skala warna kontinu
p1 = ggplot(jatim_merged) +
geom_sf(aes(fill = Rata.rata.Lama.Sekolah), color = "white", size = 0.3) +
scale_fill_viridis_c(option = "C",
name = "RLS (Tahun)",
breaks = seq(6, 10, by = 0.5)) +
labs(title = "Sebaran Rata-rata Lama Sekolah di Jawa Timur Tahun 2024",
subtitle = "Berdasarkan Kabupaten/Kota",
caption = "Sumber: BPS Jawa Timur, 2024") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5),
legend.position = "right",
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
# Peta Choropleth dengan kategori (quartile)
jatim_merged <- jatim_merged %>%
mutate(RLS_Kategori = cut(Rata.rata.Lama.Sekolah,
breaks = quantile(Rata.rata.Lama.Sekolah,
probs = c(0, 0.25, 0.5, 0.75, 1)),
labels = c("Rendah (Q1)", "Sedang (Q2)",
"Tinggi (Q3)", "Sangat Tinggi (Q4)"),
include.lowest = TRUE))
p2 = ggplot(jatim_merged) +
geom_sf(aes(fill = RLS_Kategori), color = "black", size = 0.3) +
scale_fill_manual(values = c("Rendah (Q1)" = "#ffffb2",
"Sedang (Q2)" = "#fecc5c",
"Tinggi (Q3)" = "#fd8d3c",
"Sangat Tinggi (Q4)" = "#e31a1c"),
name = "Kategori RLS") +
labs(title = "Klasifikasi Rata-rata Lama Sekolah di Jawa Timur",
subtitle = "Berdasarkan Kuartil",
caption = "Sumber: BPS Jawa Timur, 2024") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5),
legend.position = "right",
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
print(p1)
Nilai Rata-Rata Lama Sekolah (RLS) tertinggi di Provinsi Jawa Timur tahun 2024 terkonsentrasi pada wilayah perkotaan, khususnya Kota Madiun, Kota Mojokerto, Kota Malang, Kota Kediri, dan Kota Surabaya, serta kabupaten penyangga seperti Sidoarjo dan Gresik.
Sebaliknya, nilai RLS terendah ditemukan di wilayah Pulau Madura, terutama Kabupaten Sampang dan Bangkalan, serta beberapa kabupaten di kawasan Tapal Kuda seperti Probolinggo, Bondowoso, Jember, dan Situbondo.
Pola ini menunjukkan adanya kesenjangan pendidikan antarwilayah yang secara spasial membentuk pengelompokan antara daerah dengan RLS tinggi dan rendah.
4.1.2 Peta Distribusi Variabel Ekonomi dan Sosial
Selain Rata-rata Lama Sekolah, penelitian ini juga memetakan distribusi spasial variabel-variabel independen yang meliputi Pendapatan Asli Daerah (PAD), Tingkat Kemiskinan, Rasio Guru–Murid SMA, serta Angka Partisipasi Murni (APM) SMA. Peta tematik untuk masing-masing variabel digunakan untuk mengamati pola persebaran dan keterkaitannya secara visual dengan capaian RLS.
- Distribusi Pendapatan Asli Daerah (PAD)
# Peta PAD
ggplot(jatim_merged) +
geom_sf(aes(fill = PAD), color = "white", size = 0.3) +
scale_fill_viridis_c(option = "D",
name = "PAD\n(Miliar Rp)",
trans = "log10",
labels = scales::comma) +
labs(title = "Sebaran Pendapatan Asli Daerah (PAD) di Jawa Timur",
subtitle = "Skala logaritmik") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
Peta PAD menunjukkan ketimpangan ekonomi yang sangat kontras antarwilayah di Provinsi Jawa Timur. Nilai PAD tertinggi terkonsentrasi di wilayah perkotaan dan kawasan industri, terutama:
Kota Surabaya
Kabupaten Sidoarjo
Kabupaten Gresik
Kota dan Kabupaten Mojokerto
Kota Malang
Wilayah-wilayah tersebut tampak berwarna paling terang pada peta, mencerminkan kapasitas fiskal daerah yang jauh lebih besar dibandingkan wilayah lain. Sebaliknya, kabupaten di Pulau Madura (Bangkalan, Sampang, Pamekasan, Sumenep) serta beberapa wilayah tapal kuda memiliki PAD relatif rendah, yang terlihat dari warna lebih gelap.
- Distribusi Tingkat Kemiskinan
# Peta Tingkat Kemiskinan
ggplot(jatim_merged) +
geom_sf(aes(fill = Tingkat.Kemiskinan), color = "white", size = 0.3) +
scale_fill_viridis_c(option = "A",
name = "Kemiskinan\n(%)",
direction = -1) +
labs(title = "Sebaran Tingkat Kemiskinan di Jawa Timur") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
Peta tingkat kemiskinan memperlihatkan pola spasial yang hampir berlawanan dengan PAD. Wilayah dengan tingkat kemiskinan tinggi terkonsentrasi di:
Pulau Madura, terutama Sampang
Wilayah tapal kuda seperti Probolinggo, Bondowoso, dan Situbondo
Sebaliknya, wilayah perkotaan seperti Surabaya, Sidoarjo, Gresik, dan Kota Malang memiliki tingkat kemiskinan yang relatif lebih rendah. Pola ini menunjukkan adanya hubungan spasial negatif antara kapasitas ekonomi daerah dan tingkat kesejahteraan penduduk, yang berpotensi berdampak langsung pada capaian pendidikan.
- Distribusi Angka Partisipasi Murni (APM) SMA
# Peta APM SMA
ggplot(jatim_merged) +
geom_sf(aes(fill = APM.SMA), color = "white", size = 0.3) +
scale_fill_viridis_c(option = "B",
name = "APM SMA\n(%)") +
labs(title = "Sebaran Angka Partisipasi Murni (APM) SMA di Jawa Timur") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
Peta APM SMA menunjukkan bahwa partisipasi pendidikan menengah atas lebih tinggi di wilayah perkotaan dan daerah dengan infrastruktur pendidikan yang baik, seperti:
Surabaya
Sidoarjo
Gresik
Malang Raya
Sebagian wilayah Jawa Timur bagian barat
Sebaliknya, wilayah Madura dan beberapa kabupaten di bagian timur menunjukkan APM SMA yang lebih rendah, menandakan masih adanya hambatan dalam akses atau keberlanjutan pendidikan SMA, baik dari sisi ekonomi, jarak sekolah, maupun faktor sosial.
- Distribusi Rasio Guru dan Murid SMA
# Peta Rasio Guru-Murid SMA
ggplot(jatim_merged) +
geom_sf(aes(fill = Rasio.Guru.dan.Murid.SMA), color = "white", size = 0.3) +
scale_fill_viridis_c(option = "E",
name = "Rasio\nGuru:Murid",
direction = -1) +
labs(title = "Sebaran Rasio Guru dan Murid SMA di Jawa Timur",
subtitle = "Nilai lebih rendah = Rasio lebih ideal") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank()
)
Peta rasio guru–murid SMA menunjukkan ketimpangan ketersediaan tenaga pendidik antarwilayah. Wilayah dengan rasio lebih rendah (lebih ideal) umumnya berada di:
Kota-kota besar dan pusat pendidikan
Beberapa kabupaten di Jawa Timur bagian tengah
Sebaliknya, wilayah dengan rasio tinggi (kurang ideal) cenderung berada di:
Pulau Madura
Wilayah timur dan pinggiran provinsi
Rasio yang tinggi mengindikasikan beban guru yang lebih besar, yang berpotensi menurunkan kualitas pembelajaran dan berdampak pada capaian pendidikan jangka panjang.
4.1.3 Statistik Deskriptif Variabel Penelitian
Statistik deskriptif digunakan untuk memberikan gambaran umum mengenai sebaran, kecenderungan pusat, dan variasi dari setiap variabel penelitian yang dianalisis pada 38 kabupaten/kota di Provinsi Jawa Timur tahun 2024. Variabel yang dikaji meliputi Rata-rata Lama Sekolah (RLS) sebagai variabel dependen, serta Pendapatan Asli Daerah (PAD), Tingkat Kemiskinan, Rasio Guru dan Murid SMA, dan Angka Partisipasi Murni (APM) SMA sebagai variabel independen.
Berikut adalah tabel data.
vars <- c("Rata.rata.Lama.Sekolah",
"PAD", "Tingkat.Kemiskinan", "Rasio.Guru.dan.Murid.SMA", "APM.SMA"
)
data_vars = data[, vars]
data_vars## Rata.rata.Lama.Sekolah PAD Tingkat.Kemiskinan Rasio.Guru.dan.Murid.SMA
## 1 7.90 79.39 13.08 0.07487521
## 2 7.80 137.74 9.11 0.07062885
## 3 7.92 98.24 10.50 0.06319837
## 4 8.68 243.77 6.28 0.05679428
## 5 7.87 144.28 8.16 0.05988220
## 6 8.26 309.73 9.95 0.05536274
## 7 7.80 305.27 8.98 0.05670424
## 8 7.27 160.81 8.65 0.05977804
## 9 6.54 289.24 9.01 0.05886180
## 10 7.78 234.45 6.54 0.06524725
## 11 6.53 73.36 12.60 0.07548953
## 12 6.93 109.26 11.51 0.04929012
## 13 6.31 124.90 16.45 0.07096690
## 14 7.46 456.16 8.63 0.06131838
## 15 10.91 857.52 4.53 0.05608574
## 16 9.13 343.62 9.37 0.06171004
## 17 8.78 231.38 8.60 0.06390146
## 18 8.25 173.11 10.17 0.05890716
## 19 8.20 126.79 10.63 0.05847210
## 20 8.69 92.03 9.32 0.05779315
## 21 7.84 48.54 13.81 0.06523460
## 22 7.59 456.77 11.69 0.06382499
## 23 7.53 287.88 14.36 0.06276276
## 24 8.48 186.32 12.16 0.07695324
## 25 10.03 511.99 10.32 0.06242611
## 26 6.01 138.50 18.66 0.05833477
## 27 5.08 143.28 20.83 0.08087496
## 28 7.17 93.18 13.41 0.08318345
## 29 6.10 130.59 17.78 0.08626153
## 30 10.92 150.13 6.51 0.05787886
## 31 10.82 71.46 6.75 0.05335157
## 32 11.14 325.68 3.91 0.05876934
## 33 9.72 64.59 6.18 0.05345843
## 34 9.94 75.76 6.32 0.06333108
## 35 11.38 112.42 5.57 0.05819991
## 36 12.11 108.82 4.38 0.06050641
## 37 10.89 2686.57 3.96 0.05864813
## 38 9.87 90.01 3.06 0.06223565
## APM.SMA
## 1 70.62
## 2 69.24
## 3 67.86
## 4 62.73
## 5 61.20
## 6 71.08
## 7 56.25
## 8 55.18
## 9 66.56
## 10 68.11
## 11 59.75
## 12 60.77
## 13 54.20
## 14 51.01
## 15 74.93
## 16 73.87
## 17 76.01
## 18 73.30
## 19 68.28
## 20 75.78
## 21 67.93
## 22 72.41
## 23 58.68
## 24 64.38
## 25 74.89
## 26 44.25
## 27 46.89
## 28 56.76
## 29 60.20
## 30 85.68
## 31 85.84
## 32 76.48
## 33 71.45
## 34 65.85
## 35 74.97
## 36 60.81
## 37 62.48
## 38 72.37Tabel berikut menampilkan ringkasan statistik deskriptif untuk setiap variabel penelitian.
## Rata.rata.Lama.Sekolah PAD Tingkat.Kemiskinan
## Min. : 5.080 Min. : 48.54 Min. : 3.060
## 1st Qu.: 7.478 1st Qu.: 100.89 1st Qu.: 6.518
## Median : 8.060 Median : 143.78 Median : 9.215
## Mean : 8.464 Mean : 270.36 Mean : 9.782
## 3rd Qu.: 9.832 3rd Qu.: 288.90 3rd Qu.:12.043
## Max. :12.110 Max. :2686.57 Max. :20.830
## Rasio.Guru.dan.Murid.SMA APM.SMA
## Min. :0.04929 Min. :44.25
## 1st Qu.:0.05823 1st Qu.:60.34
## Median :0.06091 Median :67.89
## Mean :0.06320 Mean :66.29
## 3rd Qu.:0.06490 3rd Qu.:73.08
## Max. :0.08626 Max. :85.84Berdasarkan nilai statistik deskriptif di atas didapatkan hasil yakni.
1. Nilai Rata-rata Lama Sekolah di Jawa Timur berkisar antara 5,08 tahun hingga 12,11 tahun, dengan nilai rata-rata sebesar 8,46 tahun dan median 8,06 tahun. Kuartil pertama (Q1) berada pada 7,48 tahun, sedangkan kuartil ketiga (Q3) pada 9,83 tahun, menunjukkan adanya variasi capaian pendidikan yang cukup lebar antarwilayah.
2. Pendapatan Asli Daerah memiliki rentang nilai yang sangat lebar, dengan nilai minimum sebesar 48,54 miliar rupiah dan maksimum mencapai 2.686,57 miliar rupiah. Nilai rata-rata PAD (270,36 miliar rupiah) jauh lebih besar dibandingkan median (143,78 miliar rupiah), yang mengindikasikan distribusi PAD yang sangat menceng ke kanan (right-skewed).
3. Tingkat kemiskinan di Jawa Timur berada pada rentang 3,06% hingga 20,83%, dengan rata-rata 9,78% dan median 9,22%. Nilai kuartil pertama sebesar 6,52% dan kuartil ketiga sebesar 12,04% menunjukkan bahwa sebagian besar wilayah memiliki tingkat kemiskinan menengah, namun masih terdapat beberapa daerah dengan tingkat kemiskinan yang sangat tinggi.Sebaran ini mencerminkan adanya ketimpangan kesejahteraan ekonomi antarwilayah, yang berpotensi memengaruhi capaian pendidikan dan keberlanjutan sekolah di tingkat menengah.
4. Rasio guru dan murid SMA memiliki nilai minimum 0,049 dan maksimum 0,086, dengan rata-rata 0,063 dan median 0,061. Rentang nilai yang relatif sempit menunjukkan bahwa variasi rasio antarwilayah tidak terlalu ekstrem, namun tetap mencerminkan perbedaan dalam ketersediaan dan distribusi tenaga pendidik. Nilai rasio yang lebih tinggi mengindikasikan beban guru yang lebih besar, yang berpotensi berdampak pada kualitas proses pembelajaran di wilayah tertentu.
5. Nilai APM SMA berkisar antara 44,25% hingga 85,84%, dengan rata-rata 66,29% dan median 67,89%. Perbedaan yang cukup besar antara nilai minimum dan maksimum menunjukkan adanya kesenjangan partisipasi pendidikan SMA antar kabupaten/kota. Wilayah dengan APM tinggi umumnya merupakan daerah perkotaan dengan akses pendidikan yang lebih baik, sementara wilayah dengan APM rendah masih menghadapi tantangan dalam mempertahankan partisipasi siswa pada jenjang pendidikan menengah atas.
# Histogram
par(mfrow = c(1, 2), mar = c(5, 4, 4, 2))
hist(jatim_merged$Rata.rata.Lama.Sekolah,
breaks = 15,
main = "Distribusi Rata-rata Lama Sekolah\ndi Jawa Timur",
xlab = "Rata-rata Lama Sekolah (Tahun)",
ylab = "Frekuensi",
col = "steelblue",
border = "white")
grid()
# Boxplot
boxplot(jatim_merged$Rata.rata.Lama.Sekolah,
main = "Boxplot Rata-rata Lama Sekolah\ndi Jawa Timur",
ylab = "Rata-rata Lama Sekolah (Tahun)",
col = "lightgreen",
border = "darkgreen",
horizontal = FALSE)
grid(nx = NA, ny = NULL)
Histogram RLS menunjukkan sebaran data yang relatif tidak simetris sempurna, dengan konsentrasi utama pada rentang 7–9 tahun, sementara boxplot memperlihatkan keberadaan nilai ekstrem di wilayah dengan RLS sangat rendah maupun sangat tinggi. Hal ini mengindikasikan ketimpangan pendidikan antar kabupaten/kota, terutama antara wilayah perkotaan dan wilayah tertinggal.
Secara keseluruhan, statistik deskriptif menunjukkan bahwa seluruh variabel penelitian memiliki variasi yang cukup besar antarwilayah, terutama pada PAD, tingkat kemiskinan, dan Rata-rata Lama Sekolah. Temuan ini mengindikasikan adanya ketimpangan struktural antar kabupaten/kota di Provinsi Jawa Timur, yang memperkuat urgensi penggunaan pendekatan spasial dalam analisis selanjutnya untuk menangkap pengaruh ketergantungan antarwilayah.
4.2 Hasil Uji Autokorelasi Spasial
Sebelum dilakukan pengujian autokorelasi spasial secara formal, terlebih dahulu divisualisasikan pola keterkaitan spasial antar kabupaten/kota di Provinsi Jawa Timur melalui jaringan ketetanggaan (spatial neighborhood). Visualisasi ini bertujuan untuk memberikan gambaran awal mengenai struktur hubungan spasial yang digunakan dalam analisis selanjutnya.
Dalam penelitian ini, matriks bobot spasial dibangun menggunakan pendekatan K-Nearest Neighbors (KNN) dengan nilai k = 4, di mana setiap kabupaten/kota dihubungkan dengan empat wilayah terdekat berdasarkan jarak antar centroid. Pendekatan KNN dipilih untuk memastikan bahwa setiap unit wilayah memiliki jumlah tetangga yang sama dan tidak terdapat wilayah yang terisolasi, khususnya pada daerah pesisir dan kepulauan yang sering mengalami permasalahan isolated polygons apabila menggunakan pendekatan berbasis kontiguitas.
4.2.1 Pola Keterkaitan Spasial Antarwilayah
#Plot Ketetanggaan Spasial
# Konversi ke Spatial untuk analisis spdep
jatim_sp <- as_Spatial(jatim_merged)
row.names(jatim_sp) <- jatim_sp$NAME_FIX
# Buat matriks ketetanggaan KNN (k=4)
coords <- coordinates(jatim_sp)
knn_jatim <- knearneigh(coords, k = 4)
W <- knn2nb(knn_jatim, row.names = row.names(jatim_sp))
WL <- nb2listw(W, style = "W", zero.policy = TRUE)
# Visualisasi Jaringan KNN
par(mar = c(2, 2, 3, 2))
p3 = plot(jatim_sp, border = "gray80", col = "lightblue", axes = TRUE,
main = "Jaringan Konektivitas KNN (k=4)\nKabupaten/Kota di Jawa Timur")
# Plot edges (garis penghubung)
plot(W, coords, add = TRUE, col = "red", lwd = 1.5)
# Plot nodes (titik centroid)
points(coords, pch = 20, cex = 1, col = "darkblue")
legend("bottomright",
legend = c("Kabupaten/Kota", "Koneksi KNN", "Centroid"),
col = c("lightblue", "red", "darkblue"),
pch = c(15, NA, 20),
lty = c(NA, 1, NA),
lwd = c(NA, 1.5, NA),
bg = "white",
cex = 0.8)
## NULLGambar jaringan konektivitas KNN menunjukkan bahwa setiap kabupaten/kota di Jawa Timur memiliki hubungan spasial yang relatif merata dengan wilayah di sekitarnya. Titik biru pada peta merepresentasikan centroid masing-masing kabupaten/kota, sedangkan garis merah menunjukkan hubungan ketetanggaan berdasarkan KNN. Pola jaringan ini mencerminkan struktur ketergantungan spasial yang bersifat jarak-based, di mana interaksi antarwilayah ditentukan oleh kedekatan geografis, bukan semata-mata oleh batas administrasi.
Secara visual, jaringan KNN memperlihatkan bahwa wilayah-wilayah di Pulau Jawa memiliki konektivitas yang lebih rapat dibandingkan wilayah kepulauan di bagian timur Jawa Timur. Meskipun demikian, pendekatan KNN tetap menjamin bahwa wilayah kepulauan tetap terhubung dengan wilayah lain melalui tetangga terdekatnya. Struktur jaringan ini menjadi dasar dalam pengujian autokorelasi spasial serta pemodelan spatial econometrics pada tahap selanjutnya, karena bobot spasial yang digunakan secara langsung memengaruhi hasil uji Moran’s I, LISA, dan estimasi model spasial.
4.2.2 Uji Autokorelasi Spasial pada Variabel Rata-Rata Lama Sekolah
4.2.2.1 Autokorelasi Global
Uji autokorelasi spasial global dilakukan untuk mengetahui apakah distribusi Rata-Rata Lama Sekolah (RLS) antar kabupaten/kota di Provinsi Jawa Timur menunjukkan pola spasial yang signifikan atau tersebar secara acak. Pengujian dilakukan menggunakan dua ukuran autokorelasi spasial global, yaitu Moran’s I dan Geary’s C, dengan matriks bobot spasial berbasis K-Nearest Neighbors (KNN) dengan k = 4.
Hasil Uji Autokorelasi Global
# Moran's I
Global_Moran <- moran.test(jatim_merged$Rata.rata.Lama.Sekolah, WL, zero.policy = TRUE)
print(Global_Moran)##
## Moran I test under randomisation
##
## data: jatim_merged$Rata.rata.Lama.Sekolah
## weights: WL
##
## Moran I statistic standard deviate = 3.216, p-value = 0.0006499
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.293862921 -0.027027027 0.009955761# Geary's C
Global_Geary <- geary.test(jatim_merged$Rata.rata.Lama.Sekolah, WL, zero.policy = TRUE)
print(Global_Geary)##
## Geary C test under randomisation
##
## data: jatim_merged$Rata.rata.Lama.Sekolah
## weights: WL
##
## Geary C statistic standard deviate = 2.4206, p-value = 0.007749
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.7464504 1.0000000 0.0109723# Getis-Ord General G
Global_Getis <- globalG.test(jatim_merged$Rata.rata.Lama.Sekolah, WL, zero.policy = TRUE)## Warning in globalG.test(jatim_merged$Rata.rata.Lama.Sekolah, WL, zero.policy =
## TRUE): Binary weights recommended (especially for distance bands)##
## Getis-Ord global G statistic
##
## data: jatim_merged$Rata.rata.Lama.Sekolah
## weights: WL
##
## standard deviate = 3.6127, p-value = 0.0001515
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 2.834979e-02 2.702703e-02 1.340632e-07Hasil uji Moran’s I menunjukkan nilai statistik sebesar 0,294 dengan nilai z-score 3,216 dan p-value 0,00065. Nilai Moran’s I yang positif dan signifikan pada tingkat signifikansi 5 persen mengindikasikan adanya autokorelasi spasial positif pada variabel RLS. Hal ini berarti kabupaten/kota dengan nilai RLS tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki RLS tinggi, demikian pula wilayah dengan RLS rendah cenderung berdekatan dengan wilayah dengan RLS rendah.
Hasil ini diperkuat oleh uji Geary’s C yang menghasilkan nilai statistik sebesar 0,746 dengan p-value 0,00775. Nilai Geary’s C yang lebih kecil dari 1 dan signifikan menunjukkan adanya kemiripan nilai antarwilayah yang bertetangga, sehingga mengonfirmasi keberadaan ketergantungan spasial pada distribusi RLS. Secara konseptual, Geary’s C yang signifikan menegaskan bahwa perbedaan RLS antarwilayah tidak bersifat acak, melainkan dipengaruhi oleh kondisi wilayah sekitarnya.
Secara keseluruhan, hasil uji autokorelasi spasial global menunjukkan bahwa distribusi Rata-Rata Lama Sekolah di Provinsi Jawa Timur tahun 2024 memiliki pola pengelompokan spasial yang signifikan. Temuan ini mengindikasikan bahwa asumsi independensi spasial pada model regresi klasik berpotensi dilanggar, sehingga penggunaan analisis autokorelasi lokal dan pemodelan spatial econometrics pada tahap selanjutnya menjadi relevan dan diperlukan.
4.2.2.2 Analisis LISA
Peta Klaster LISA
Local_Moran <- localmoran(jatim_merged$Rata.rata.Lama.Sekolah, WL, zero.policy = TRUE)
colnames(Local_Moran) <- c("Ii", "E.Ii", "Var.Ii", "Z.Ii", "P.value")
jatim_merged$Ii <- Local_Moran[, "Ii"]
jatim_merged$Z.Ii <- Local_Moran[, "Z.Ii"]
jatim_merged$P.value <- Local_Moran[, "P.value"]
# Klasifikasi Klaster LISA
scaled_x <- as.numeric(scale(jatim_merged$Rata.rata.Lama.Sekolah))
lag_x <- lag.listw(WL, scaled_x)
jatim_merged$lisa_cluster <- "Non-signifikan"
jatim_merged$lisa_cluster[jatim_merged$P.value <= 0.05 & scaled_x > 0 & lag_x > 0] <- "High-High"
jatim_merged$lisa_cluster[jatim_merged$P.value <= 0.05 & scaled_x < 0 & lag_x < 0] <- "Low-Low"
jatim_merged$lisa_cluster[jatim_merged$P.value <= 0.05 & scaled_x > 0 & lag_x < 0] <- "High-Low"
jatim_merged$lisa_cluster[jatim_merged$P.value <= 0.05 & scaled_x < 0 & lag_x > 0] <- "Low-High"
# Plot LISA Map
ggplot(jatim_merged) +
geom_sf(aes(fill = lisa_cluster), color = "white", size = 0.2) +
scale_fill_manual(values = c(
"High-High" = "#ff0000",
"Low-Low" = "#0000ff",
"High-Low" = "#ff9999",
"Low-High" = "#9999ff",
"Non-signifikan" = "#eeeeee"
)) +
labs(title = "Peta Klaster Autokorelasi Lokal (LISA)",
subtitle = "Variabel: Rata.rata.Lama.Sekolah - Jawa Timur",
fill = "Tipe Klaster") +
theme_minimal()
Analisis autokorelasi spasial lokal dilakukan menggunakan Local Indicators of Spatial Association (LISA) untuk mengidentifikasi pola keterkaitan spasial pada tingkat kabupaten/kota secara lebih rinci. Berbeda dengan autokorelasi global yang hanya menunjukkan keberadaan pola spasial secara umum, analisis LISA memungkinkan identifikasi lokasi spesifik klaster Rata-Rata Lama Sekolah (RLS) yang signifikan secara statistik.
Berdasarkan peta klaster LISA, teridentifikasi beberapa tipe klaster spasial, yaitu High–High (HH), Low–Low (LL), Low–High (LH), serta wilayah non-signifikan. Klaster High–High menunjukkan wilayah dengan nilai RLS tinggi yang dikelilingi oleh wilayah lain dengan RLS tinggi, sedangkan klaster Low–Low menunjukkan wilayah dengan RLS rendah yang berdekatan dengan wilayah lain yang juga memiliki RLS rendah. Kedua klaster ini mencerminkan adanya pengelompokan spasial yang kuat dan konsisten secara lokal.
Selain itu, terdapat klaster Low–High, yaitu wilayah dengan RLS relatif rendah yang dikelilingi oleh wilayah dengan RLS tinggi. Keberadaan klaster ini mengindikasikan adanya ketimpangan lokal atau kondisi wilayah yang tertinggal di tengah lingkungan dengan capaian pendidikan yang lebih baik. Sementara itu, sebagian besar wilayah lainnya tergolong ke dalam kategori non-signifikan, yang menunjukkan tidak adanya keterkaitan spasial lokal yang signifikan pada tingkat signifikansi yang digunakan.
Secara keseluruhan, hasil analisis LISA memperkuat temuan autokorelasi spasial global sebelumnya dengan menunjukkan bahwa ketergantungan spasial Rata-Rata Lama Sekolah tidak hanya bersifat umum, tetapi juga muncul secara spesifik pada wilayah-wilayah tertentu. Temuan ini penting sebagai dasar dalam perumusan kebijakan pendidikan berbasis wilayah, karena menunjukkan bahwa intervensi kebijakan perlu difokuskan pada klaster wilayah dengan RLS rendah, khususnya yang tergolong dalam klaster Low–Low dan Low–High.
4.2.2.3 Hotspot Analysis (Getis-Ord Gi*)
Peta Hotspot dan Coldspot
x_val <- jatim_merged$Rata.rata.Lama.Sekolah
x_val[is.na(x_val)] <- 0
lwB <- nb2listw(W, style = "B", zero.policy = TRUE)
Wb_mat <- listw2mat(lwB)
# Hitung G dan G*
sum_x <- sum(x_val)
num_G <- as.numeric(Wb_mat %*% x_val)
den_G <- (sum_x - x_val)
G_raw <- num_G / den_G
Wb_star <- Wb_mat
diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_val)
den_Gs <- sum_x
G_star_raw <- num_Gs / den_Gs
# Hitung z-score
Gz_res <- as.numeric(localG(x_val, WL, zero.policy = TRUE))
jatim_merged <- jatim_merged %>%
mutate(
G_raw = G_raw,
G_star_raw = G_star_raw,
z_Gistar = Gz_res,
hotspot_cat = case_when(
z_Gistar >= 1.96 ~ "Hot spot (p<=0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p<=0.05)",
TRUE ~ "Tidak Signifikan"
)
)
# Ringkasan
summary(select(jatim_merged, G_raw, G_star_raw, z_Gistar))## G_raw G_star_raw z_Gistar geometry
## Min. :0.07983 Min. :0.09729 Min. :-2.8116 MULTIPOLYGON :38
## 1st Qu.:0.10536 1st Qu.:0.12746 1st Qu.:-0.2778 epsg:NA : 0
## Median :0.11559 Median :0.13847 Median : 0.7204 +proj=long...: 0
## Mean :0.11206 Mean :0.13539 Mean : 0.3767
## 3rd Qu.:0.12275 3rd Qu.:0.14566 3rd Qu.: 1.4098
## Max. :0.13370 Max. :0.16245 Max. : 2.4740##
## Cold spot (p<=0.05) Hot spot (p<=0.05) Tidak Signifikan
## 5 5 28# Plot Getis-Ord G* Raw
p1 <- ggplot(jatim_merged) +
geom_sf(aes(fill = G_star_raw), color = "white", size = 0.2) +
scale_fill_viridis_c(option = "C", na.value = "grey90") +
labs(title = "Raw Getis–Ord G* (proporsi massa tetangga)",
subtitle = "Variabel: Rata.rata.Lama.Sekolah - Jawa Timur",
fill = "G*_raw") +
theme_minimal()
# Plot Hot/Cold Spots
jatim_merged$hotspot_cat <- factor(jatim_merged$hotspot_cat,
levels = c("Hot spot (p<=0.05)",
"Cold spot (p<=0.05)",
"Tidak Signifikan"))
p2 <- ggplot(jatim_merged) +
geom_sf(aes(fill = hotspot_cat), color = "white", size = 0.2) +
scale_fill_manual(values = c("Hot spot (p<=0.05)" = "#d73027",
"Cold spot (p<=0.05)" = "#4575b4",
"Tidak Signifikan" = "#e0e0e0")) +
labs(title = "Getis–Ord Gi* — Hot/Cold Spots",
subtitle = "Rata.rata.Lama.Sekolah - Jawa Timur (alpha = 0.05)",
fill = NULL) +
theme_minimal()
print(p1)
Analisis hotspot dan coldspot dilakukan menggunakan statistik Getis–Ord Gi* untuk mengidentifikasi wilayah yang menjadi pusat konsentrasi nilai Rata-Rata Lama Sekolah (RLS) tinggi maupun rendah secara spasial. Berbeda dengan analisis LISA yang menekankan kesamaan lokal antarwilayah, metode Gi* berfokus pada intensitas pengelompokan nilai tinggi atau rendah dalam suatu lingkungan spasial.
Peta Getis–Ord G*_raw menunjukkan variasi proporsi massa spasial RLS pada masing-masing wilayah, di mana nilai G*_raw yang lebih tinggi mengindikasikan kecenderungan wilayah tersebut berada di lingkungan dengan nilai RLS yang relatif tinggi. Sebaliknya, nilai G*_raw yang lebih rendah mencerminkan lingkungan spasial dengan konsentrasi RLS yang lebih rendah. Visualisasi ini memberikan gambaran awal mengenai gradien spasial capaian pendidikan di Provinsi Jawa Timur.
Selanjutnya, peta Getis–Ord Gi* dengan tingkat signifikansi α = 0,05 mengklasifikasikan wilayah ke dalam tiga kategori, yaitu hot spot, cold spot, dan wilayah tidak signifikan. Wilayah yang tergolong sebagai hot spot merupakan wilayah dengan nilai RLS tinggi yang dikelilingi oleh wilayah lain dengan nilai RLS tinggi secara signifikan. Sebaliknya, cold spot menunjukkan wilayah dengan nilai RLS rendah yang terkonsentrasi secara spasial dengan wilayah lain yang juga memiliki nilai RLS rendah.
Hasil analisis Gi* menunjukkan bahwa hot spot dan cold spot RLS di Provinsi Jawa Timur tidak tersebar secara acak, melainkan membentuk pola spasial yang terlokalisasi. Temuan ini sejalan dengan hasil uji autokorelasi spasial global dan analisis LISA sebelumnya, yang sama-sama mengindikasikan adanya ketergantungan spasial pada capaian pendidikan antar kabupaten/kota.
Secara keseluruhan, analisis hotspot dan coldspot memperkuat bukti bahwa disparitas Rata-Rata Lama Sekolah di Provinsi Jawa Timur memiliki dimensi spasial yang kuat. Informasi ini penting dalam konteks kebijakan pendidikan berbasis wilayah, karena memungkinkan identifikasi area prioritas yang memerlukan intervensi khusus, baik untuk peningkatan capaian pendidikan pada wilayah cold spot maupun untuk mempertahankan dan memperkuat capaian pada wilayah hot spot.
Setelah mengidentifikasi pola pengelompokan spasial RLS melalui analisis autokorelasi global dan lokal, tahap selanjutnya adalah mengevaluasi apakah ketergantungan spasial tersebut juga muncul dalam hubungan antara RLS dan variabel-variabel penentunya melalui pemodelan spatial econometrics.
4.2.2.4 Autokorelasi Variabel X
Selain variabel dependen Rata-Rata Lama Sekolah (RLS), pengujian autokorelasi spasial juga dilakukan terhadap seluruh variabel independen untuk mengetahui apakah faktor-faktor ekonomi dan pendidikan tersebut memiliki pola keterkaitan spasial antar kabupaten/kota di Provinsi Jawa Timur. Pengujian dilakukan menggunakan ukuran autokorelasi spasial global, yaitu Moran’s I, Geary’s C, dan Getis–Ord General G, dengan matriks bobot spasial berbasis KNN (k = 4).
Ringkasan Autokorelasi Variabel Independen
# Daftar variabel X
x_vars <- c(
"PAD", "Tingkat.Kemiskinan", "Rasio.Guru.dan.Murid.SMA", "APM.SMA"
)
# Buat list kosong untuk menyimpan hasil global
hasil_moran <- list()
hasil_geary <- list()
hasil_getis <- list()
# Loop semua variabel
for (v in x_vars) {
cat("\n=====================================\n")
cat("🔹 Variabel:", v, "\n")
cat("=====================================\n")
# Skip kalau variabel tidak ada
if (!v %in% names(jatim_merged)) {
cat("❌ Variabel", v, "tidak ditemukan di data.\n")
next
}
# -------------------- Global Tests --------------------
moran_res <- moran.test(jatim_merged[[v]], WL, zero.policy = TRUE)
geary_res <- geary.test(jatim_merged[[v]], WL, zero.policy = TRUE)
getis_res <- globalG.test(jatim_merged[[v]], WL, zero.policy = TRUE)
hasil_moran[[v]] <- moran_res
hasil_geary[[v]] <- geary_res
hasil_getis[[v]] <- getis_res
cat("Moran's I :", round(moran_res$estimate[1], 4),
"| p =", round(moran_res$p.value, 4), "\n")
cat("Geary's C :", round(geary_res$estimate[1], 4),
"| p =", round(geary_res$p.value, 4), "\n")
cat("Getis-Ord G :", round(getis_res$estimate[1], 4),
"| p =", round(getis_res$p.value, 4), "\n")
# -------------------- LISA (Local Moran's I) --------------------
lisa <- localmoran(jatim_merged[[v]], WL, zero.policy = TRUE)
colnames(lisa) <- c("Ii", "E.Ii", "Var.Ii", "Z.Ii", "P.value")
# Salin data agar tidak menimpa variabel lain
temp_map <- jatim_merged
temp_map$Ii <- lisa[, "Ii"]
temp_map$P.value <- lisa[, "P.value"]
temp_map$lag_var <- lag.listw(WL, jatim_merged[[v]])
mean_var <- mean(jatim_merged[[v]], na.rm = TRUE)
mean_lag <- mean(temp_map$lag_var, na.rm = TRUE)
temp_map$cluster <- NA
temp_map$cluster[jatim_merged[[v]] >= mean_var & temp_map$lag_var >= mean_lag] <- "High-High"
temp_map$cluster[jatim_merged[[v]] <= mean_var & temp_map$lag_var <= mean_lag] <- "Low-Low"
temp_map$cluster[jatim_merged[[v]] >= mean_var & temp_map$lag_var <= mean_lag] <- "High-Low"
temp_map$cluster[jatim_merged[[v]] <= mean_var & temp_map$lag_var >= mean_lag] <- "Low-High"
temp_map$cluster[temp_map$P.value > 0.05] <- "Non-signifikan"
# -------------------- Plot LISA --------------------
print(
ggplot(temp_map) +
geom_sf(aes(fill = cluster)) +
scale_fill_manual(values = c(
"High-High" = "red",
"Low-Low" = "blue",
"High-Low" = "orange",
"Low-High" = "green",
"Non-signifikan" = "grey80"
)) +
labs(
title = paste("Peta Klaster LISA —", v),
subtitle = "Autokorelasi Lokal (Local Moran's I) - Jawa Timur",
fill = "Tipe Klaster"
) +
theme_minimal()
)
}##
## =====================================
## 🔹 Variabel: PAD
## =====================================## Warning in globalG.test(jatim_merged[[v]], WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)## Moran's I : 0.1022 | p = 0.007
## Geary's C : 1.0948 | p = 0.7026
## Getis-Ord G : 0.0417 | p = 5e-04
##
## =====================================
## 🔹 Variabel: Tingkat.Kemiskinan
## =====================================## Warning in globalG.test(jatim_merged[[v]], WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)## Moran's I : 0.2696 | p = 0.0013
## Geary's C : 0.6885 | p = 0.0019
## Getis-Ord G : 0.0265 | p = 0.721
##
## =====================================
## 🔹 Variabel: Rasio.Guru.dan.Murid.SMA
## =====================================## Warning in globalG.test(jatim_merged[[v]], WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)## Moran's I : 0.1804 | p = 0.017
## Geary's C : 0.6761 | p = 0.0016
## Getis-Ord G : 0.0265 | p = 0.9893
##
## =====================================
## 🔹 Variabel: APM.SMA
## =====================================## Warning in globalG.test(jatim_merged[[v]], WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)## Moran's I : 0.245 | p = 0.003
## Geary's C : 0.7719 | p = 0.0162
## Getis-Ord G : 0.0276 | p = 0.0202
# ============================================================
# RINGKASAN HASIL GLOBAL (opsional)
# ============================================================
summary_autokorelasi <- data.frame(
Variabel = names(hasil_moran),
Moran_I = sapply(hasil_moran, \(x) x$estimate[1]),
Moran_p = sapply(hasil_moran, \(x) x$p.value),
Geary_C = sapply(hasil_geary, \(x) x$estimate[1]),
Geary_p = sapply(hasil_geary, \(x) x$p.value),
Getis_G = sapply(hasil_getis, \(x) x$estimate[1]),
Getis_p = sapply(hasil_getis, \(x) x$p.value)
)
print(summary_autokorelasi)## Variabel Moran_I
## PAD.Moran I statistic PAD 0.1021606
## Tingkat.Kemiskinan.Moran I statistic Tingkat.Kemiskinan 0.2696372
## Rasio.Guru.dan.Murid.SMA.Moran I statistic Rasio.Guru.dan.Murid.SMA 0.1803834
## APM.SMA.Moran I statistic APM.SMA 0.2449683
## Moran_p Geary_C Geary_p
## PAD.Moran I statistic 0.006978270 1.0948400 0.702552881
## Tingkat.Kemiskinan.Moran I statistic 0.001321039 0.6884914 0.001918682
## Rasio.Guru.dan.Murid.SMA.Moran I statistic 0.016964528 0.6761142 0.001623255
## APM.SMA.Moran I statistic 0.003025488 0.7718704 0.016201430
## Getis_G Getis_p
## PAD.Moran I statistic 0.04168956 0.0004758545
## Tingkat.Kemiskinan.Moran I statistic 0.02650502 0.7210072145
## Rasio.Guru.dan.Murid.SMA.Moran I statistic 0.02647954 0.9893281457
## APM.SMA.Moran I statistic 0.02755199 0.0202478064Hasil uji menunjukkan bahwa seluruh variabel independen memiliki indikasi autokorelasi spasial, meskipun dengan tingkat dan karakteristik yang berbeda. Variabel Pendapatan Asli Daerah (PAD) menunjukkan nilai Moran’s I positif dan signifikan (I = 0,102; p < 0,01), yang mengindikasikan adanya pengelompokan spasial PAD antarwilayah. Namun, uji Geary’s C tidak signifikan, sementara Getis–Ord G signifikan. Hal ini menunjukkan bahwa pola spasial PAD lebih terlihat pada skala agregat (global clustering), tetapi tidak terlalu kuat pada tingkat perbedaan lokal antarwilayah yang berdekatan.
Variabel Tingkat Kemiskinan menunjukkan autokorelasi spasial yang kuat dan konsisten. Nilai Moran’s I positif dan signifikan (I = 0,270; p < 0,01) serta nilai Geary’s C yang signifikan (< 1) mengindikasikan adanya kemiripan nilai kemiskinan antarwilayah yang bertetangga. Temuan ini menunjukkan bahwa wilayah dengan tingkat kemiskinan tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki tingkat kemiskinan tinggi, demikian pula sebaliknya.
Autokorelasi spasial juga teridentifikasi pada variabel Rasio Guru–Murid SMA, dengan nilai Moran’s I dan Geary’s C yang signifikan. Hal ini mengindikasikan bahwa distribusi tenaga pendidik pada jenjang SMA tidak bersifat acak secara spasial, melainkan membentuk pola pengelompokan antarwilayah. Namun, uji Getis–Ord G yang tidak signifikan menunjukkan bahwa pengelompokan tersebut tidak membentuk pusat konsentrasi nilai tinggi atau rendah yang kuat secara regional.
Variabel Angka Partisipasi Murni (APM) SMA menunjukkan autokorelasi spasial yang relatif konsisten pada ketiga ukuran yang digunakan. Nilai Moran’s I dan Geary’s C yang signifikan menunjukkan adanya keterkaitan lokal antarwilayah, sementara signifikansi Getis–Ord G mengindikasikan adanya konsentrasi wilayah dengan tingkat partisipasi sekolah yang relatif tinggi maupun rendah.
Secara keseluruhan, hasil uji autokorelasi spasial pada variabel independen menunjukkan bahwa faktor-faktor ekonomi dan pendidikan di Provinsi Jawa Timur tidak terdistribusi secara acak, melainkan memiliki pola spasial tertentu. Temuan ini menguatkan justifikasi penggunaan model regresi spasial, karena variabel penentu RLS tidak hanya dipengaruhi oleh karakteristik internal suatu wilayah, tetapi juga oleh kondisi wilayah sekitarnya.
Berdasarkan hasil uji autokorelasi spasial pada variabel dependen dan independen, tahap selanjutnya adalah mengevaluasi hubungan kausal antara RLS dan variabel penentunya menggunakan pemodelan spatial econometrics.
4.3 Pemodelan Spatial Econometrics
4.3.1 Model OLS
Tahap awal pemodelan dilakukan dengan mengestimasi model regresi linier berganda menggunakan metode Ordinary Least Squares (OLS). Model OLS digunakan sebagai model dasar untuk mengidentifikasi pengaruh variabel independen, yaitu Pendapatan Asli Daerah (PAD), Tingkat Kemiskinan, Rasio Guru–Murid SMA, dan Angka Partisipasi Murni (APM) SMA terhadap Rata-Rata Lama Sekolah (RLS) pada kabupaten/kota di Provinsi Jawa Timur.
Hasil Estimasi OLS
ols_model <- lm(Rata.rata.Lama.Sekolah ~
PAD +
Tingkat.Kemiskinan +
Rasio.Guru.dan.Murid.SMA +
APM.SMA,
data = jatim_merged
)
summary(ols_model)##
## Call:
## lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1222 -0.5166 -0.0056 0.5606 2.7075
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.6398873 2.0136712 3.297 0.00234 **
## PAD 0.0004552 0.0003716 1.225 0.22921
## Tingkat.Kemiskinan -0.2486402 0.0531182 -4.681 4.71e-05 ***
## Rasio.Guru.dan.Murid.SMA 4.1776684 22.0245809 0.190 0.85072
## APM.SMA 0.0583680 0.0196086 2.977 0.00542 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9133 on 33 degrees of freedom
## Multiple R-squared: 0.7436, Adjusted R-squared: 0.7125
## F-statistic: 23.92 on 4 and 33 DF, p-value: 2.346e-09# -------------------- Uji Asumsi OLS --------------------
# Uji Multikolinearitas
cat("\n========== Uji Multikolinearitas (VIF) ==========\n")##
## ========== Uji Multikolinearitas (VIF) ==========vif_check <- vif(lm((Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged))
print(vif_check)## PAD Tingkat.Kemiskinan APM.SMA
## 1.149498 2.222186 1.545794
## Rasio.Guru.dan.Murid.SMA
## 1.552363## Warning: package 'nortest' was built under R version 4.5.2##
## ========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(ols_model)
## D = 0.093751, p-value = 0.5448##
## ========== Uji Homoskedastisitas (Breusch-Pagan) ==========##
## studentized Breusch-Pagan test
##
## data: ols_model
## BP = 5.8286, df = 4, p-value = 0.2123##
## ========== Uji Linearitas (RESET Test) ==========##
## RESET test
##
## data: ols_model
## RESET = 2.9753, df1 = 2, df2 = 31, p-value = 0.06576# 4. Uji Autokorelasi Spasial Residual
cat("\n========== Uji Autokorelasi Spasial Residual (Moran's I) ==========\n")##
## ========== Uji Autokorelasi Spasial Residual (Moran's I) ==========jatim_merged$residual_ols <- residuals(ols_model)
moran_res <- moran.test(jatim_merged$residual_ols, WL, zero.policy = TRUE)
print(moran_res)##
## Moran I test under randomisation
##
## data: jatim_merged$residual_ols
## weights: WL
##
## Moran I statistic standard deviate = 2.0352, p-value = 0.02091
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.169235399 -0.027027027 0.009299124Berdasarkan hasil estimasi, persamaan regresi OLS yang diperoleh adalah sebagai berikut:
\[\hat{y} = 6.64 + 0.000468 PAD - 2.46Tingkat Kemiskinan + 4.18Rasio Guru Murid SMA + 0.0583 APM.SMA \]
Hasil estimasi menunjukkan bahwa Tingkat Kemiskinan dan APM SMA berpengaruh signifikan terhadap RLS pada taraf signifikansi 5 persen. Variabel Tingkat Kemiskinan memiliki koefisien negatif dan signifikan, yang mengindikasikan bahwa peningkatan persentase penduduk miskin cenderung menurunkan capaian rata-rata lama sekolah. Sebaliknya, APM SMA memiliki pengaruh positif dan signifikan, yang menunjukkan bahwa semakin tinggi tingkat partisipasi pendidikan menengah, semakin tinggi pula RLS di suatu wilayah.
Sementara itu, variabel PAD dan Rasio Guru–Murid SMA tidak menunjukkan pengaruh yang signifikan secara statistik dalam model OLS. Hal ini mengindikasikan bahwa kapasitas fiskal daerah serta rasio tenaga pendidik SMA belum secara langsung menjelaskan variasi RLS antarwilayah apabila pengaruh spasial belum diperhitungkan.
Secara keseluruhan, model OLS memiliki nilai Adjusted R-squared sebesar 0,7125, yang menunjukkan bahwa sekitar 71,25 persen variasi RLS dapat dijelaskan oleh variabel independen dalam model. Uji simultan (F-test) juga menunjukkan bahwa model signifikan secara keseluruhan dengan p-value < 0,001.
Uji Asumsi Klasik Model OLS
Pengujian asumsi klasik dilakukan untuk menilai kelayakan model OLS secara statistik. Hasil uji multikolinearitas berdasarkan nilai Variance Inflation Factor (VIF) menunjukkan bahwa seluruh variabel memiliki nilai VIF di bawah 10, sehingga tidak terdapat indikasi multikolinearitas serius.
Uji normalitas residual dilakukan menggunakan uji Kolmogorov–Smirnov dengan koreksi Lilliefors, yang menghasilkan p-value sebesar 0,5448. Nilai ini lebih besar dari taraf signifikansi 5 persen, sehingga residual model dapat dianggap berdistribusi normal.
Uji homoskedastisitas menggunakan Breusch–Pagan Test menghasilkan p-value sebesar 0,2123, yang menunjukkan tidak terdapat indikasi heteroskedastisitas. Selain itu, uji linearitas model menggunakan RESET Test menghasilkan p-value sebesar 0,0658, sehingga tidak terdapat bukti kuat adanya kesalahan spesifikasi nonlinier pada model.
Berdasarkan hasil tersebut, dapat disimpulkan bahwa model OLS memenuhi sebagian besar asumsi klasik regresi linear.
Uji Autokorelasi Spasial Residual OLS
Meskipun asumsi klasik terpenuhi, hasil uji Moran’s I terhadap residual model OLS menunjukkan adanya autokorelasi spasial yang signifikan. Nilai Moran’s I sebesar 0,169 dengan p-value 0,0209 mengindikasikan adanya autokorelasi spasial positif pada residual, di mana kesalahan prediksi di suatu wilayah cenderung berkorelasi dengan wilayah sekitarnya.
Temuan ini menunjukkan bahwa asumsi independensi spasial pada model OLS dilanggar. Dengan demikian, model OLS belum sepenuhnya mampu menangkap ketergantungan spasial yang terdapat dalam data RLS antarwilayah di Provinsi Jawa Timur.
Oleh karena itu, diperlukan estimasi model spatial econometrics pada tahap selanjutnya untuk mengakomodasi ketergantungan spasial, baik yang bersumber dari variabel dependen maupun dari komponen error. Sebagai dasar pemilihan model spasial yang tepat, tahap berikutnya dilakukan uji Lagrange Multiplier (LM) dan Robust LM.
4.3.2 Lagrange Multiplier Tests
Uji Lagrange Multiplier (LM) dilakukan untuk menentukan bentuk ketergantungan spasial yang paling sesuai dalam pemodelan, yaitu apakah ketergantungan spasial muncul pada variabel dependen (Spatial Lag) atau pada komponen error (Spatial Error).
Hasil LM Tests
##
## ========== Lagrange Multiplier Tests ==========##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
## test weights: WL
##
## RSerr = 2.5256, df = 1, p-value = 0.112
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
## test weights: WL
##
## RSlag = 4.5358, df = 1, p-value = 0.03319
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
## test weights: WL
##
## adjRSerr = 0.20342, df = 1, p-value = 0.652
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
## test weights: WL
##
## adjRSlag = 2.2136, df = 1, p-value = 0.1368
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Rata.rata.Lama.Sekolah ~ PAD + Tingkat.Kemiskinan +
## Rasio.Guru.dan.Murid.SMA + APM.SMA, data = jatim_merged)
## test weights: WL
##
## SARMA = 4.7392, df = 2, p-value = 0.09352Uji Lagrange Multiplier (LM) dilakukan untuk menentukan bentuk ketergantungan spasial yang paling sesuai dalam pemodelan, yaitu apakah ketergantungan spasial muncul pada variabel dependen (spatial lag) atau pada komponen error (spatial error). Pengujian ini dilakukan berdasarkan residual dari model Ordinary Least Squares (OLS) dengan menggunakan matriks bobot spasial berbasis k-nearest neighbors (KNN) dengan k = 4.
Hasil uji LM standar menunjukkan bahwa LM-Lag (RSlag) memiliki nilai statistik sebesar 4,5358 dengan p-value = 0,0332, yang signifikan pada taraf 5 persen. Sementara itu, LM-Error (RSerr) memiliki nilai statistik sebesar 2,5256 dengan p-value = 0,112, sehingga tidak signifikan secara statistik. Temuan ini mengindikasikan bahwa ketergantungan spasial lebih dominan muncul melalui variabel dependen, yaitu adanya pengaruh nilai Rata-rata Lama Sekolah (RLS) di suatu wilayah terhadap RLS di wilayah tetangganya.
Selanjutnya, hasil Robust Lagrange Multiplier menunjukkan bahwa baik Robust LM-Lag (adjRSlag) maupun Robust LM-Error (adjRSerr) tidak signifikan pada taraf 5 persen, dengan p-value masing-masing sebesar 0,1368 dan 0,6520. Selain itu, uji SARMA yang menguji kemungkinan ketergantungan spasial gabungan (lag dan error) juga tidak signifikan dengan p-value sebesar 0,0935. Ketidaksignifikanan uji robust ini mengindikasikan bahwa bentuk ketergantungan spasial yang terdeteksi relatif lemah dan sensitif terhadap spesifikasi model.
Meskipun demikian, hasil uji Moran’s I pada residual model OLS pada bagian sebelumnya menunjukkan adanya autokorelasi spasial positif yang signifikan, yang menandakan bahwa asumsi independensi residual masih dilanggar. Selain itu, analisis autokorelasi global dan lokal pada variabel dependen serta beberapa variabel independen juga memperlihatkan adanya pola spasial yang konsisten antarwilayah.
Oleh karena itu, penelitian ini tidak hanya bergantung pada hasil uji LM dalam menentukan spesifikasi model, tetapi tetap melanjutkan estimasi model spatial econometrics yang lebih umum, khususnya Spatial Durbin Model (SDM) dan turunannya. Pendekatan ini dipilih untuk mengakomodasi kemungkinan adanya efek limpahan (spillover effects) baik melalui variabel dependen maupun variabel independen, serta untuk menghindari potensi bias spesifikasi model akibat keterbatasan daya uji LM pada ukuran sampel yang relatif kecil.
Berdasarkan hasil uji Lagrange Multiplier dan indikasi autokorelasi spasial pada residual OLS, tahap selanjutnya dilakukan estimasi dan perbandingan beberapa model spatial econometrics guna menentukan spesifikasi model yang paling sesuai.
4.3.3 Estimasi Model Spasial
4.3.3.1 Spatial Durbin Model (SDM)
##
## ========== SPATIAL DURBIN MODEL (SDM) ==========sdm_model <- lagsarlm(
(Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged,
listw = WL,
Durbin = TRUE,
method = "eigen",
zero.policy = TRUE
)
summary(sdm_model)##
## Call:lagsarlm(formula = (Rata.rata.Lama.Sekolah) ~ PAD + Tingkat.Kemiskinan +
## APM.SMA + Rasio.Guru.dan.Murid.SMA, data = jatim_merged,
## listw = WL, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.642366 -0.466435 -0.074493 0.357917 2.590046
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.30147083 4.61454491 0.2820 0.7779154
## PAD 0.00012333 0.00031601 0.3903 0.6963412
## Tingkat.Kemiskinan -0.28794396 0.04610414 -6.2455 4.224e-10
## APM.SMA 0.05691455 0.01624299 3.5039 0.0004584
## Rasio.Guru.dan.Murid.SMA 7.35174642 19.15947754 0.3837 0.7011910
## lag.PAD 0.00144444 0.00054392 2.6556 0.0079166
## lag.Tingkat.Kemiskinan 0.12478513 0.10051713 1.2414 0.2144464
## lag.APM.SMA 0.02468830 0.03043072 0.8113 0.4171961
## lag.Rasio.Guru.dan.Murid.SMA 11.05072542 55.76540019 0.1982 0.8429163
##
## Rho: 0.20221, LR test value: 0.66987, p-value: 0.4131
## Asymptotic standard error: 0.21239
## z-value: 0.95206, p-value: 0.34107
## Wald statistic: 0.90642, p-value: 0.34107
##
## Log likelihood: -41.85355 for mixed model
## ML residual variance (sigma squared): 0.52572, (sigma: 0.72506)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 105.71, (AIC for lm: 104.38)
## LM test for residual autocorrelation
## test value: 4.5043, p-value: 0.033809##
## --- Uji Multikolinearitas SDM ---print(vif(lm((Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged)))## PAD Tingkat.Kemiskinan APM.SMA
## 1.149498 2.222186 1.545794
## Rasio.Guru.dan.Murid.SMA
## 1.552363library(nortest)
cat("\n========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========\n")##
## ========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(sdm_model)
## D = 0.098443, p-value = 0.4649##
## --- Uji Homoskedastisitas SDM ---##
## studentized Breusch-Pagan test
##
## data:
## BP = 11.973, df = 8, p-value = 0.1524##
## --- Uji Autokorelasi Spasial Residual SDM ---##
## Moran I test under randomisation
##
## data: residuals(sdm_model)
## weights: WL
##
## Moran I statistic standard deviate = -0.44701, p-value = 0.6726
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.069379061 -0.027027027 0.008976726Spatial Durbin Model (SDM) diestimasi untuk mengakomodasi ketergantungan spasial yang teridentifikasi pada residual model OLS, sekaligus untuk menangkap kemungkinan adanya efek limpahan (spillover effects) dari variabel independen antarwilayah. Model ini memasukkan lag spasial pada variabel dependen maupun variabel independen.
Hasil estimasi SDM menunjukkan bahwa model memiliki nilai Akaike Information Criterion (AIC) sebesar 105,71, yang lebih rendah dibandingkan model OLS (AIC = 104,38 pada data tanpa komponen spasial), sehingga menunjukkan adanya peningkatan kecocokan model setelah memperhitungkan struktur spasial.
Dari sisi diagnostik residual, uji normalitas residual menggunakan Kolmogorov–Smirnov dengan koreksi Lilliefors menghasilkan p-value sebesar 0,4649, sehingga residual dapat dianggap berdistribusi normal. Uji homoskedastisitas menggunakan Breusch–Pagan Test memberikan p-value sebesar 0,1524, yang menunjukkan tidak terdapat indikasi heteroskedastisitas.
Uji Moran’s I terhadap residual SDM menghasilkan p-value sebesar 0,6726, yang menunjukkan tidak adanya autokorelasi spasial residual yang signifikan. Hasil ini mengindikasikan bahwa SDM berhasil mengakomodasi ketergantungan spasial yang sebelumnya terdeteksi pada residual model OLS.
4.3.3.2 Spatial Durbin Error Model (SDEM)
##
## ========== SPATIAL DURBIN ERROR MODEL (SDEM) ==========sdem_model <- errorsarlm(
(Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged,
listw = WL,
Durbin = TRUE,
method = "eigen",
zero.policy = TRUE
)
summary(sdem_model)##
## Call:
## errorsarlm(formula = (Rata.rata.Lama.Sekolah) ~ PAD + Tingkat.Kemiskinan +
## APM.SMA + Rasio.Guru.dan.Murid.SMA, data = jatim_merged,
## listw = WL, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.736257 -0.497481 -0.023573 0.379901 2.591738
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.85994866 4.61312818 0.4032 0.6868114
## PAD 0.00022141 0.00030948 0.7154 0.4743435
## Tingkat.Kemiskinan -0.28378719 0.04654213 -6.0974 1.078e-09
## APM.SMA 0.05885337 0.01639248 3.5903 0.0003303
## Rasio.Guru.dan.Murid.SMA 10.55642650 19.35951786 0.5453 0.5855585
## lag.PAD 0.00164316 0.00051670 3.1801 0.0014722
## lag.Tingkat.Kemiskinan 0.07502788 0.08723980 0.8600 0.3897786
## lag.APM.SMA 0.04086588 0.02624398 1.5572 0.1194342
## lag.Rasio.Guru.dan.Murid.SMA 12.94177513 56.43387642 0.2293 0.8186153
##
## Lambda: 0.014729, LR test value: 0.0027334, p-value: 0.9583
## Asymptotic standard error: 0.24559
## z-value: 0.059971, p-value: 0.95218
## Wald statistic: 0.0035965, p-value: 0.95218
##
## Log likelihood: -42.18711 for error model
## ML residual variance (sigma squared): 0.53927, (sigma: 0.73435)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 106.37, (AIC for lm: 104.38)# Uji Asumsi SDEM
library(nortest)
cat("\n========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========\n")##
## ========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(sdem_model)
## D = 0.11213, p-value = 0.2652##
## --- Uji Homoskedastisitas SDEM ---##
## studentized Breusch-Pagan test
##
## data:
## BP = 12.642, df = 8, p-value = 0.1248##
## --- Uji Autokorelasi Spasial Residual SDEM ---##
## Moran I test under randomisation
##
## data: residuals(sdem_model)
## weights: WL
##
## Moran I statistic standard deviate = 0.28365, p-value = 0.3883
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0001292606 -0.0270270270 0.0089922234Spatial Durbin Error Model (SDEM) diestimasi untuk mengakomodasi ketergantungan spasial yang bersumber dari komponen error serta kemungkinan adanya efek limpahan (spillover effects) dari variabel independen antarwilayah. Berbeda dengan SDM, SDEM tidak memasukkan lag spasial dari variabel dependen, melainkan memfokuskan pada autokorelasi spasial pada residual dan pengaruh tidak langsung dari variabel penjelas.
Hasil estimasi menunjukkan bahwa model SDEM memiliki nilai Akaike Information Criterion (AIC) sebesar 106,37, yang lebih rendah dibandingkan model OLS (AIC = 104,38), namun lebih tinggi dibandingkan model SDM (AIC = 105,71). Hal ini menunjukkan bahwa SDEM mampu memperbaiki model OLS dengan mempertimbangkan struktur spasial, tetapi peningkatan kecocokan model masih lebih terbatas dibandingkan SDM.
Dari sisi diagnostik residual, uji normalitas residual menggunakan Kolmogorov–Smirnov dengan koreksi Lilliefors menghasilkan p-value sebesar 0,2652, sehingga residual dapat dianggap berdistribusi normal. Uji homoskedastisitas menggunakan Breusch–Pagan Test memberikan p-value sebesar 0,1248, yang menunjukkan tidak terdapat indikasi heteroskedastisitas.
Uji Moran’s I terhadap residual SDEM menghasilkan p-value sebesar 0,3883, yang menunjukkan tidak adanya autokorelasi spasial residual yang signifikan. Temuan ini mengindikasikan bahwa SDEM berhasil menghilangkan autokorelasi spasial pada residual yang sebelumnya teridentifikasi pada model OLS.
4.3.3.3 Spatial Autoregressive Combined Model (SAC)
##
## ========== SPATIAL AUTOREGRESSIVE COMBINED (SAC) ==========sac_model <- sacsarlm(
(Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged,
listw = WL,
method = "eigen",
zero.policy = TRUE
)
summary(sac_model)##
## Call:sacsarlm(formula = (Rata.rata.Lama.Sekolah) ~ PAD + Tingkat.Kemiskinan +
## APM.SMA + Rasio.Guru.dan.Murid.SMA, data = jatim_merged,
## listw = WL, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.682902 -0.542240 0.060242 0.564194 2.574586
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.6781e+00 2.5667e+00 2.9914 0.002777
## PAD -8.1208e-05 3.2511e-04 -0.2498 0.802751
## Tingkat.Kemiskinan -2.7997e-01 4.6696e-02 -5.9957 2.026e-09
## APM.SMA 4.7329e-02 1.7131e-02 2.7627 0.005732
## Rasio.Guru.dan.Murid.SMA -2.1857e-01 1.9006e+01 -0.0115 0.990825
##
## Rho: 0.041867
## Asymptotic standard error: 0.18415
## z-value: 0.22736, p-value: 0.82015
## Lambda: 0.48362
## Asymptotic standard error: 0.21448
## z-value: 2.2549, p-value: 0.024141
##
## LR test value: 4.5483, p-value: 0.10289
##
## Log likelihood: -45.51688 for sac model
## ML residual variance (sigma squared): 0.60977, (sigma: 0.78088)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 107.03, (AIC for lm: 107.58)# Uji Asumsi SAC
library(nortest)
cat("\n========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========\n")##
## ========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(sac_model)
## D = 0.067884, p-value = 0.9297##
## --- Uji Homoskedastisitas SAC ---##
## studentized Breusch-Pagan test
##
## data:
## BP = 9.574, df = 4, p-value = 0.04825##
## --- Uji Autokorelasi Spasial Residual SAC ---##
## Moran I test under randomisation
##
## data: residuals(sac_model)
## weights: WL
##
## Moran I statistic standard deviate = 0.084411, p-value = 0.4664
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.018862722 -0.027027027 0.009354959Spatial Autoregressive Combined Model (SAC) diestimasi untuk menangkap dua bentuk ketergantungan spasial secara simultan, yaitu ketergantungan spasial pada variabel dependen (spatial lag) dan ketergantungan spasial pada komponen error (spatial error). Dengan demikian, model SAC bersifat lebih umum dibandingkan SAR dan SEM karena mengombinasikan kedua mekanisme tersebut dalam satu spesifikasi model.
Hasil estimasi menunjukkan bahwa model SAC memiliki nilai Akaike Information Criterion (AIC) sebesar 107,03, yang lebih rendah dibandingkan model OLS (AIC = 107,58), sehingga mengindikasikan adanya perbaikan kecocokan model setelah memperhitungkan struktur spasial. Namun demikian, nilai AIC SAC masih lebih tinggi dibandingkan model SDM dan SDEM, sehingga dari sisi kecocokan model, SAC belum menjadi spesifikasi yang paling optimal.
Hasil uji normalitas residual menggunakan Kolmogorov–Smirnov dengan koreksi Lilliefors menghasilkan p-value sebesar 0,9297, yang menunjukkan bahwa residual berdistribusi normal. Namun, uji homoskedastisitas menggunakan Breusch–Pagan Test menghasilkan p-value sebesar 0,0483, yang mengindikasikan adanya gejala heteroskedastisitas pada taraf signifikansi 5 persen.
Uji Moran’s I terhadap residual model SAC menghasilkan p-value sebesar 0,4664, yang menunjukkan tidak terdapat autokorelasi spasial residual yang signifikan. Hal ini mengindikasikan bahwa model SAC berhasil menghilangkan autokorelasi spasial residual yang sebelumnya terdeteksi pada model OLS.
4.3.3.4 General Nested Spatial Model (GNS)
##
## ========== GENERAL NESTED SPATIAL MODEL (GNS) ==========gns_model <- sacsarlm(
(Rata.rata.Lama.Sekolah) ~
PAD +
Tingkat.Kemiskinan +
APM.SMA +
Rasio.Guru.dan.Murid.SMA,
data = jatim_merged,
listw = WL,
Durbin = TRUE,
method = "eigen",
zero.policy = TRUE
)
summary(gns_model)##
## Call:sacsarlm(formula = (Rata.rata.Lama.Sekolah) ~ PAD + Tingkat.Kemiskinan +
## APM.SMA + Rasio.Guru.dan.Murid.SMA, data = jatim_merged,
## listw = WL, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44043 -0.34053 -0.06858 0.29107 2.29880
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.85896571 3.28712083 -0.2613 0.7938515
## PAD -0.00012397 0.00031249 -0.3967 0.6915721
## Tingkat.Kemiskinan -0.29976368 0.04505967 -6.6526 2.88e-11
## APM.SMA 0.05530039 0.01641055 3.3698 0.0007522
## Rasio.Guru.dan.Murid.SMA 1.05276393 17.53700374 0.0600 0.9521309
## lag.PAD 0.00107962 0.00049428 2.1842 0.0289462
## lag.Tingkat.Kemiskinan 0.25653440 0.07915741 3.2408 0.0011919
## lag.APM.SMA -0.01457970 0.02395561 -0.6086 0.5427809
## lag.Rasio.Guru.dan.Murid.SMA 11.05623435 47.75691867 0.2315 0.8169181
##
## Rho: 0.70681
## Asymptotic standard error: 0.15225
## z-value: 4.6424, p-value: 3.4433e-06
## Lambda: -0.97166
## Asymptotic standard error: 0.34168
## z-value: -2.8438, p-value: 0.0044584
##
## LR test value: 15.107, p-value: 0.019442
##
## Log likelihood: -40.23758 for sacmixed model
## ML residual variance (sigma squared): 0.36781, (sigma: 0.60648)
## Number of observations: 38
## Number of parameters estimated: 12
## AIC: 104.48, (AIC for lm: 107.58)# Uji Asumsi GNS
library(nortest)
cat("\n========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========\n")##
## ========== Uji Normalitas Residual (Kolmogorov–Smirnov / Lilliefors) ==========##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(gns_model)
## D = 0.11858, p-value = 0.1947##
## --- Uji Homoskedastisitas GNS ---##
## studentized Breusch-Pagan test
##
## data:
## BP = 11.501, df = 8, p-value = 0.1749##
## --- Uji Autokorelasi Spasial Residual GNS ---##
## Moran I test under randomisation
##
## data: residuals(gns_model)
## weights: WL
##
## Moran I statistic standard deviate = 0.01589, p-value = 0.4937
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.025549744 -0.027027027 0.008643697General Nested Spatial Model (GNS) diestimasi sebagai spesifikasi spasial paling umum karena mampu mengakomodasi seluruh bentuk ketergantungan spasial secara simultan, yaitu spatial lag pada variabel dependen, spatial error pada komponen galat, serta spillover effects dari variabel independen melalui lag spasial. Model ini mencakup dan menurunkan model SAR, SEM, SDM, dan SAC sebagai kasus khusus, sehingga digunakan untuk mengevaluasi apakah struktur spasial yang lebih kompleks mampu memberikan peningkatan kecocokan model yang signifikan.
Hasil estimasi menunjukkan bahwa model GNS memiliki nilai Akaike Information Criterion (AIC) sebesar 104,48, yang merupakan nilai AIC terendah dibandingkan seluruh model yang telah diestimasi sebelumnya, termasuk OLS, SDM, SDEM, dan SAC. Penurunan nilai AIC yang cukup signifikan ini menunjukkan bahwa model GNS memberikan kecocokan terbaik terhadap data Rata-rata Lama Sekolah (RLS) antar kabupaten/kota di Provinsi Jawa Timur, meskipun dengan jumlah parameter yang lebih banyak.
Selain itu, nilai log-likelihood model GNS (−40,24) lebih tinggi dibandingkan model spasial lainnya, yang semakin menguatkan bahwa spesifikasi ini mampu menangkap struktur variasi data secara lebih komprehensif.
Dari sisi diagnostik residual, uji normalitas residual menggunakan Kolmogorov–Smirnov dengan koreksi Lilliefors menghasilkan p-value sebesar 0,1947, yang menunjukkan bahwa residual berdistribusi normal. Uji homoskedastisitas menggunakan Breusch–Pagan Test memberikan p-value sebesar 0,1749, sehingga tidak terdapat indikasi heteroskedastisitas.
Selanjutnya, uji Moran’s I terhadap residual model GNS menghasilkan p-value sebesar 0,4937, yang menunjukkan tidak adanya autokorelasi spasial residual yang signifikan. Hal ini mengindikasikan bahwa model GNS berhasil mengakomodasi seluruh bentuk ketergantungan spasial yang sebelumnya terdeteksi pada model OLS dan model spasial lainnya.
4.3.4 Pemilihan Model Terbaik
Pemilihan model spasial terbaik dilakukan dengan membandingkan kinerja beberapa model, yaitu Ordinary Least Squares (OLS), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Autoregressive Combined Model (SAC), dan General Nested Spatial Model (GNS).
Tabel berikut menyajikan perbandingan nilai AIC dari masing-masing model:
| Model | AIC |
|---|---|
| OLS | 107.58 |
| SDM | 105.71 |
| SDEM | 106.37 |
| SAC | 107.03 |
| GNS | 104.48 |
Kriteria pemilihan model meliputi Akaike Information Criterion (AIC), keberadaan autokorelasi spasial residual, serta pemenuhan asumsi statistik model.
Berdasarkan Tabel perbandingan AIC, model OLS memiliki nilai AIC terbesar, yaitu 107,58, serta pada tahap sebelumnya terbukti masih mengandung autokorelasi spasial pada residual. Hal ini menunjukkan bahwa model non-spasial belum mampu menangkap struktur ketergantungan antarwilayah secara memadai, sehingga kurang tepat digunakan sebagai model akhir.
Model SDM, SDEM, dan SAC menunjukkan peningkatan kinerja dibandingkan OLS, yang tercermin dari nilai AIC yang lebih kecil. Ketiga model tersebut juga berhasil menghilangkan autokorelasi spasial pada residual serta memenuhi asumsi normalitas dan homoskedastisitas. Namun, jika dibandingkan satu sama lain, masih terdapat perbedaan tingkat kecocokan model.
Model General Nested Spatial Model (GNS) menghasilkan nilai AIC terendah, yaitu 104,48, yang menunjukkan bahwa model ini memiliki kecocokan terbaik di antara seluruh model yang diuji. Selain itu, hasil diagnostik residual menunjukkan bahwa residual model GNS tidak mengandung autokorelasi spasial serta memenuhi asumsi statistik yang diperlukan.
Dengan demikian, berdasarkan perbandingan nilai AIC, hasil uji diagnostik residual, dan kemampuan model dalam menangkap ketergantungan spasial yang kompleks, General Nested Spatial Model (GNS) dipilih sebagai model spasial terbaik untuk menganalisis faktor-faktor yang memengaruhi Rata-Rata Lama Sekolah di Provinsi Jawa Timur. Model ini selanjutnya digunakan sebagai dasar dalam interpretasi hasil dan pembahasan implikasi kebijakan.
4.3.5 Interpretasi Model Terbaik
General Nested Spatial Model (GNS) dipilih sebagai model terbaik karena mampu menangkap seluruh bentuk ketergantungan spasial yang mungkin terjadi, baik melalui lag variabel dependen, lag variabel independen (spillover effects), maupun autokorelasi pada komponen error. Model ini memberikan representasi paling komprehensif terhadap struktur spasial Rata-Rata Lama Sekolah (RLS) di Provinsi Jawa Timur.
Pengaruh Variabel Independen (Efek Langsung)
Hasil estimasi menunjukkan bahwa Tingkat Kemiskinan memiliki koefisien negatif dan signifikan (β = −0,2998; p-value < 0,001). Temuan ini mengindikasikan bahwa peningkatan persentase penduduk miskin di suatu kabupaten/kota secara signifikan menurunkan capaian rata-rata lama sekolah. Secara substantif, kondisi ekonomi yang lemah membatasi akses pendidikan melalui keterbatasan biaya, kebutuhan bekerja dini, serta rendahnya daya dukung rumah tangga terhadap keberlanjutan pendidikan.
Variabel Angka Partisipasi Murni (APM) SMA berpengaruh positif dan signifikan terhadap RLS (β = 0,0553; p-value < 0,01). Hal ini menunjukkan bahwa semakin tinggi tingkat partisipasi penduduk usia sekolah pada jenjang SMA, semakin tinggi pula rata-rata lama sekolah yang dicapai di wilayah tersebut. Temuan ini menegaskan pentingnya pendidikan menengah sebagai determinan utama peningkatan capaian pendidikan formal.
Sementara itu, Pendapatan Asli Daerah (PAD) dan Rasio Guru–Murid SMA tidak menunjukkan pengaruh signifikan secara langsung terhadap RLS. Hal ini mengindikasikan bahwa kapasitas fiskal daerah dan rasio tenaga pendidik belum secara otomatis berkontribusi pada peningkatan lama sekolah tanpa diiringi oleh efektivitas kebijakan dan pemerataan akses pendidikan.
Efek Limpahan Spasial (Spillover Effects)
Model GNS memungkinkan identifikasi efek tidak langsung melalui lag spasial variabel independen. Hasil estimasi menunjukkan bahwa lag PAD memiliki koefisien positif dan signifikan (θ = 0,00108; p-value < 0,05). Artinya, peningkatan PAD di wilayah sekitar berkontribusi positif terhadap RLS di suatu daerah. Temuan ini mengindikasikan adanya spillover fiskal, di mana daerah dengan kapasitas ekonomi tinggi dapat memberikan dampak tidak langsung bagi wilayah tetangganya melalui interaksi ekonomi, infrastruktur bersama, dan mobilitas penduduk.
Selain itu, lag Tingkat Kemiskinan juga berpengaruh positif dan signifikan (θ = 0,2565; p-value < 0,01). Interpretasi temuan ini menunjukkan bahwa tingkat kemiskinan di wilayah sekitar turut memengaruhi kondisi pendidikan suatu daerah, yang mencerminkan keterkaitan regional dalam dinamika sosial ekonomi. Wilayah dengan konsentrasi kemiskinan tinggi dapat membentuk klaster spasial yang memengaruhi kualitas pendidikan lintas batas administratif.
Sementara itu, lag APM SMA dan lag Rasio Guru–Murid SMA tidak signifikan, yang menunjukkan bahwa efek limpahan pendidikan menengah dan tenaga pendidik lebih bersifat lokal dibandingkan regional.
Ketergantungan Spasial Global
Parameter spasial ρ (rho) bernilai positif dan signifikan (ρ = 0,7068; p-value < 0,001), yang menunjukkan adanya ketergantungan spasial kuat pada variabel dependen. Artinya, capaian RLS di suatu wilayah sangat dipengaruhi oleh capaian RLS di wilayah sekitarnya, sehingga peningkatan pendidikan di satu daerah berpotensi mendorong peningkatan pendidikan di daerah lain.
Parameter λ (lambda) bernilai negatif dan signifikan (λ = −0,9717; p-value < 0,01), yang mengindikasikan adanya korelasi spasial pada komponen error. Hal ini menunjukkan bahwa terdapat faktor-faktor tak teramati yang berpola spasial dan memengaruhi RLS, seperti budaya pendidikan, karakteristik geografis, atau kualitas institusi pendidikan yang belum sepenuhnya tertangkap oleh variabel dalam model.
Diagnostik Model
Hasil uji diagnostik menunjukkan bahwa model GNS memenuhi asumsi statistik yang diperlukan. Uji normalitas residual menggunakan Kolmogorov–Smirnov dengan koreksi Lilliefors menghasilkan p-value sebesar 0,1947, sehingga residual dapat dianggap berdistribusi normal. Uji homoskedastisitas Breusch–Pagan juga tidak signifikan (p-value = 0,1749), yang menunjukkan tidak adanya heteroskedastisitas.
Selain itu, uji Moran’s I terhadap residual menghasilkan p-value sebesar 0,4937, yang mengindikasikan tidak terdapat autokorelasi spasial residual yang signifikan. Temuan ini menegaskan bahwa model GNS telah berhasil mengakomodasi ketergantungan spasial yang sebelumnya terdeteksi pada model OLS.
Secara keseluruhan, hasil estimasi General Nested Spatial Model menunjukkan bahwa capaian rata-rata lama sekolah di Provinsi Jawa Timur dipengaruhi oleh faktor ekonomi dan pendidikan secara langsung, serta oleh interaksi spasial antarwilayah. Oleh karena itu, kebijakan peningkatan pendidikan tidak dapat dirancang secara terisolasi per kabupaten/kota, melainkan perlu mempertimbangkan keterkaitan regional dan efek limpahan antarwilayah.
4.4 Interpolasi Spasial
Interpolasi spasial dilakukan untuk mengestimasi nilai Rata-Rata Lama Sekolah (RLS) pada lokasi yang tidak terobservasi secara langsung, sehingga diperoleh gambaran permukaan spasial RLS yang lebih kontinu di seluruh wilayah Provinsi Jawa Timur. Metode interpolasi yang digunakan ditentukan berdasarkan karakteristik tren spasial pada data.
4.4.1 Uji Tren Spasial
Sebelum dilakukan interpolasi, terlebih dahulu dilakukan uji tren spasial untuk menentukan apakah terdapat kecenderungan tren global pada data RLS. Uji ini dilakukan dengan meregresikan nilai RLS terhadap koordinat spasial (X dan Y) dari titik centroid kabupaten/kota.
## Warning: package 'gstat' was built under R version 4.5.2## Warning: package 'raster' was built under R version 4.5.1##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(sp)
library(sf)
# ------------------------------------------------------------------------------
# PERSIAPAN DATA TITIK (CENTROID)
# ------------------------------------------------------------------------------
# Transformasi ke CRS proyeksi (UTM 49S untuk Jatim)
jatim_merged_proj <- st_transform(jatim_merged, 32749)
# Buat centroid dari polygon
jatim_centroid <- st_centroid(jatim_merged_proj)## Warning: st_centroid assumes attributes are constant over geometries# Konversi ke objek Spatial (untuk kompatibilitas gstat)
jatim_centroid_sp <- as(jatim_centroid, "Spatial")
# Ekstrak koordinat dan tambahkan sebagai atribut
coords_krig <- coordinates(jatim_centroid_sp)
jatim_centroid_sp$X <- coords_krig[, 1]
jatim_centroid_sp$Y <- coords_krig[, 2]
# Ringkasan variabel target
cat("\n=== RINGKASAN DATA ===\n")##
## === RINGKASAN DATA ===## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.080 7.478 8.060 8.464 9.832 12.110# ------------------------------------------------------------------------------
# UJI TREN SPASIAL (PENENTU ORDINARY vs UNIVERSAL KRIGING)
# ------------------------------------------------------------------------------
cat("\n=== UJI TREN SPASIAL ===\n")##
## === UJI TREN SPASIAL ===trend_test <- lm(
Rata.rata.Lama.Sekolah ~ X + Y,
data = jatim_centroid_sp@data
)
summary(trend_test)##
## Call:
## lm(formula = Rata.rata.Lama.Sekolah ~ X + Y, data = jatim_centroid_sp@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7119 -1.1514 -0.6224 1.6330 2.8857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.152e+01 5.415e+01 0.582 0.5642
## X -6.834e-06 2.879e-06 -2.374 0.0232 *
## Y -2.019e-06 5.919e-06 -0.341 0.7351
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.622 on 35 degrees of freedom
## Multiple R-squared: 0.1423, Adjusted R-squared: 0.09325
## F-statistic: 2.902 on 2 and 35 DF, p-value: 0.06819Hasil uji tren menunjukkan bahwa koordinat X (arah barat–timur) berpengaruh signifikan terhadap RLS (p-value = 0,023), sedangkan koordinat Y (arah selatan–utara) tidak signifikan (p-value = 0,735). Temuan ini mengindikasikan adanya tren spasial lemah satu arah, di mana nilai RLS cenderung berubah mengikuti posisi geografis barat–timur.
Meskipun nilai koefisien determinasi relatif kecil (R² = 0,142), keberadaan koefisien koordinat yang signifikan menunjukkan bahwa asumsi rata-rata konstan secara global tidak sepenuhnya terpenuhi. Oleh karena itu, metode Universal Kriging (UK) lebih sesuai digunakan dibandingkan Ordinary Kriging, karena UK mampu mengakomodasi tren spasial melalui fungsi deterministik berbasis koordinat.
Sebagai pembanding, metode Inverse Distance Weighting (IDW) juga diterapkan untuk mengevaluasi kinerja interpolasi non-parametrik berbasis jarak.
4.4.2 Analisis Variogram
Plot Variogram Empiris dan Fitted
##
## === VARIOGRAM ANALYSIS ===# Variogram empiris dengan tren (untuk Universal Kriging)
vgm_uk <- variogram(
Rata.rata.Lama.Sekolah ~ X + Y,
jatim_centroid_sp
)
# Plot variogram empiris
plot(vgm_uk, main = "Empirical Variogram")
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit##
## Parameter Variogram:## model psill range
## 1 Nug 6.128128 0.00
## 2 Sph 0.000000 41847.04
Variogram empiris menunjukkan bahwa semivarians Rata-Rata Lama Sekolah meningkat pada jarak antarwilayah yang relatif dekat, menandakan adanya ketergantungan spasial lokal. Setelah jarak tertentu, peningkatan semivarians cenderung melandai, yang mengindikasikan tercapainya range pengaruh spasial.
Model variogram spherical dipilih sebagai model terbaik. Hasil estimasi parameter menunjukkan adanya komponen nugget sebesar 6,13 yang merepresentasikan variasi mikro atau kesalahan pengukuran pada jarak sangat dekat. Sementara itu, komponen sill parsial (psill) pada struktur spherical bernilai 0 dengan range sekitar 41.847 meter, yang menunjukkan bahwa pengaruh spasial Rata-Rata Lama Sekolah masih terdeteksi hingga jarak tersebut sebelum menjadi relatif konstan.
Kesesuaian antara variogram empiris dan variogram fitted mengindikasikan bahwa model spherical mampu menangkap struktur spasial data secara memadai, sehingga layak digunakan sebagai dasar dalam proses interpolasi spasial menggunakan metode Kriging.
4.4.3 Hasil Interpolasi
# Konversi boundary ke Spatial
jatim_poly_sp <- as(jatim_merged_proj, "Spatial")
# Buat grid reguler dalam bounding box Jatim
bbox_jatim <- bbox(jatim_poly_sp)
grid_krig <- expand.grid(
X = seq(bbox_jatim[1, 1], bbox_jatim[1, 2], length.out = 200),
Y = seq(bbox_jatim[2, 1], bbox_jatim[2, 2], length.out = 200)
)
# Konversi ke SpatialPointsDataFrame
coordinates(grid_krig) <- ~ X + Y
gridded(grid_krig) <- TRUE
proj4string(grid_krig) <- CRS(proj4string(jatim_poly_sp))Universal Kriging menghasilkan permukaan yang lebih halus dan gradual karena mempertimbangkan tren spasial serta struktur variogram. Pola ini lebih stabil untuk merepresentasikan kecenderungan regional Rata-rata Lama Sekolah.
Peta Interpolasi UK
##
## === UNIVERSAL KRIGING ===uk_result <- krige(
Rata.rata.Lama.Sekolah ~ X + Y,
locations = jatim_centroid_sp,
newdata = grid_krig,
model = vgm_uk_fit
)## [using universal kriging]# Konversi ke raster dan clip ke batas Jatim
r_uk <- raster(uk_result["var1.pred"])
r_uk <- crop(r_uk, extent(jatim_poly_sp))
r_uk <- mask(r_uk, jatim_poly_sp)
par(mar = c(4, 4, 3, 6))
plot(r_uk,
main = "Interpolasi Rata-rata Lama Sekolah\nUniversal Kriging",
col = terrain.colors(50),
xlab = "Longitude",
ylab = "Latitude",
legend.args = list(text = "Tahun", side = 4, line = 2.5),
cex.main = 1.3)
Peta Interpolasi IDW
##
## === INVERSE DISTANCE WEIGHTING ===idw_result <- idw(
Rata.rata.Lama.Sekolah ~ 1,
locations = jatim_centroid_sp,
newdata = grid_krig,
idp = 2 # power parameter
)## [inverse distance weighted interpolation]# Konversi ke raster dan clip
r_idw <- raster(idw_result["var1.pred"])
r_idw <- crop(r_idw, extent(jatim_poly_sp))
r_idw <- mask(r_idw, jatim_poly_sp)
# =============================
par(mar = c(4, 4, 3, 6))
plot(r_idw,
main = "Interpolasi Rata-rata Lama Sekolah\nInverse Distance Weighting (IDW)",
col = terrain.colors(50),
xlab = "Longitude",
ylab = "Latitude",
legend.args = list(text = "Tahun", side = 4, line = 2.5),
cex.main = 1.3)
IDW menghasilkan permukaan interpolasi yang lebih lokal dan mengikuti pola titik pengamatan, sehingga terlihat variasi yang lebih tajam antarwilayah. Metode ini menekankan pengaruh jarak tanpa mempertimbangkan struktur spasial global. Peta Ketidakpastian (Kriging Variance)
Secara visual, Universal Kriging memberikan hasil interpolasi yang lebih konsisten secara spasial, sedangkan IDW lebih sensitif terhadap nilai ekstrem di sekitar titik pengamatan.
##
## === PETA UNCERTAINTY (KRIGING VARIANCE) ===par(mar = c(4, 4, 3, 6))
# Konversi variance ke raster
r_uk_var <- raster(uk_result["var1.var"])
r_uk_var <- crop(r_uk_var, extent(jatim_poly_sp))
r_uk_var <- mask(r_uk_var, jatim_poly_sp)
# Plot variance (uncertainty)
plot(r_uk_var,
main = "Peta Ketidakpastian (Kriging Variance)\nUniversal Kriging",
col = heat.colors(50),
xlab = "Longitude",
ylab = "Latitude",
legend.args = list(text = "Variance", side = 4, line = 2.5),
cex.main = 1.3)
Peta Ketidakpastian (Kriging Variance) pada metode Universal Kriging menunjukkan tingkat ketelitian hasil interpolasi di setiap lokasi. Nilai varians kriging yang lebih rendah (warna lebih gelap) umumnya terdapat pada wilayah yang dekat dengan titik pengamatan, menandakan estimasi yang lebih andal. Sebaliknya, varians yang lebih tinggi muncul pada area yang relatif jauh dari titik observasi atau memiliki kepadatan data rendah, sehingga tingkat ketidakpastiannya lebih besar. Pola ini menunjukkan bahwa keandalan hasil interpolasi sangat dipengaruhi oleh distribusi spasial data Rata-rata Lama Sekolah.
4.4.4 Validasi Model Interpolasi
##
## === VALIDASI MODEL (LEAVE-ONE-OUT CROSS-VALIDATION) ===# Cross-validation Universal Kriging
cv_uk <- krige.cv(
Rata.rata.Lama.Sekolah ~ X + Y,
jatim_centroid_sp,
model = vgm_uk_fit
)
# Cross-validation IDW
cv_idw <- krige.cv(
Rata.rata.Lama.Sekolah ~ 1,
locations = jatim_centroid_sp,
nmax = 12,
set = list(idp = 2)
)
# Hitung metrik validasi untuk UK
rmse_uk <- sqrt(mean(cv_uk$residual^2))
mae_uk <- mean(abs(cv_uk$residual))
# Hitung metrik validasi untuk IDW
rmse_idw <- sqrt(mean(cv_idw$residual^2))
mae_idw <- mean(abs(cv_idw$residual))
# Tampilkan hasil
cat("\n--- METRIK VALIDASI ---\n")##
## --- METRIK VALIDASI ---## UK -> RMSE: 1.6721 | MAE: 1.4457## IDW -> RMSE: 1.7631 | MAE: 1.5407Hasil validasi menunjukkan bahwa Universal Kriging memiliki kinerja yang lebih baik dibandingkan IDW, ditunjukkan oleh nilai kesalahan yang lebih rendah.
Tabel Metrik Validasi UK vs IDW
| Metode | RMSE | MAE |
|---|---|---|
| UK | 1.6721 | 1.4457 |
| IDW | 1.7631 | 1.5407 |
Nilai RMSE dan MAE yang lebih kecil pada UK mengindikasikan bahwa metode ini lebih akurat dalam memprediksi Rata-rata Lama Sekolah dibandingkan IDW.
Residual Plot
# =============================
par(mar = c(5, 5, 4, 2))
plot(cv_uk$observed, cv_uk$residual,
main = "Residual Plot - Universal Kriging",
xlab = "Nilai Observasi (Tahun)",
ylab = "Residual (Tahun)",
pch = 19,
col = "steelblue",
cex = 1.5,
cex.main = 1.3,
cex.lab = 1.1)
abline(h = 0, col = "red", lty = 2, lwd = 2)
grid()
# Tambahkan info RMSE
text(min(cv_uk$observed), max(cv_uk$residual) * 0.9,
sprintf("RMSE = %.4f", rmse_uk),
pos = 4, col = "darkblue", font = 2, cex = 1.1)
par(mar = c(5, 5, 4, 2))
plot(cv_idw$observed, cv_idw$residual,
main = "Residual Plot - IDW",
xlab = "Nilai Observasi (Tahun)",
ylab = "Residual (Tahun)",
pch = 19,
col = "darkgreen",
cex = 1.5,
cex.main = 1.3,
cex.lab = 1.1)
abline(h = 0, col = "red", lty = 2, lwd = 2)
grid()
# Tambahkan info RMSE
text(min(cv_idw$observed), max(cv_idw$residual) * 0.9,
sprintf("RMSE = %.4f", rmse_idw),
pos = 4, col = "darkred", font = 2, cex = 1.1)
Plot residual menunjukkan bahwa:
Residual Universal Kriging cenderung lebih terkonsentrasi di sekitar nol dengan sebaran yang relatif lebih stabil. Sedangkan residual IDW menunjukkan sebaran yang lebih lebar, terutama pada nilai observasi rendah dan tinggi, yang mengindikasikan potensi bias lokal.
Observed vs Predicted
# UK
par(mar = c(5, 5, 4, 2))
plot(cv_uk$observed, cv_uk$var1.pred,
main = "Observed vs Predicted - Universal Kriging",
xlab = "Nilai Observasi (Tahun)",
ylab = "Nilai Prediksi (Tahun)",
pch = 19,
col = "steelblue",
cex = 1.5,
cex.main = 1.3,
cex.lab = 1.1,
xlim = range(cv_uk$observed),
ylim = range(cv_uk$observed))
abline(0, 1, col = "red", lwd = 2, lty = 2)
grid()
# IDW
par(mar = c(5, 5, 4, 2))
plot(cv_idw$observed, cv_idw$var1.pred,
main = "Observed vs Predicted - IDW",
xlab = "Nilai Observasi (Tahun)",
ylab = "Nilai Prediksi (Tahun)",
pch = 19,
col = "darkgreen",
cex = 1.5,
cex.main = 1.3,
cex.lab = 1.1,
xlim = range(cv_idw$observed),
ylim = range(cv_idw$observed))
abline(0, 1, col = "red", lwd = 2, lty = 2)
grid()
Grafik perbandingan nilai observasi dan prediksi memperlihatkan bahwa titik hasil Universal Kriging lebih dekat terhadap garis diagonal (garis ideal 1:1), menandakan kesesuaian prediksi yang lebih baik.
Pada metode IDW, deviasi dari garis diagonal lebih besar, terutama pada wilayah dengan nilai Rata-rata Lama Sekolah yang ekstrem
Berdasarkan hasil LOOCV, analisis residual, dan kesesuaian antara nilai observasi dan prediksi, Universal Kriging dipilih sebagai metode interpolasi terbaik dalam penelitian ini. Metode ini lebih mampu menangkap struktur spasial Rata-rata Lama Sekolah di Provinsi Jawa Timur serta menghasilkan prediksi yang lebih stabil dan akurat dibandingkan IDW.
4.4.5 Distribusi Nilai Prediksi
Histogram Perbandingan UK vs IDW
par(mar = c(5, 4, 4, 2))
# Ekstrak nilai dari raster
vals_uk <- values(r_uk)
vals_idw <- values(r_idw)
# Hapus NA
vals_uk <- vals_uk[!is.na(vals_uk)]
vals_idw <- vals_idw[!is.na(vals_idw)]
# Plot histogram
hist(vals_uk,
col = rgb(0, 0, 1, 0.5),
main = "Distribusi Nilai Prediksi UK vs IDW",
xlab = "Rata-rata Lama Sekolah (Tahun)",
ylab = "Frekuensi",
xlim = range(c(vals_uk, vals_idw)),
breaks = 30,
border = "blue",
lwd = 2)
hist(vals_idw,
col = rgb(1, 0, 0, 0.5),
add = TRUE,
breaks = 30,
border = "red",
lwd = 2)
legend("topright",
legend = c("Universal Kriging", "IDW"),
fill = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
border = c("blue", "red"),
lwd = 2,
bg = "white",
box.lwd = 1.5)
Distribusi nilai prediksi Rata-rata Lama Sekolah dianalisis untuk membandingkan karakteristik hasil interpolasi antara metode Universal Kriging (UK) dan Inverse Distance Weighting (IDW). Histogram menunjukkan bahwa kedua metode menghasilkan pola distribusi yang relatif serupa, dengan konsentrasi nilai prediksi utama berada pada kisaran 7–9 tahun, yang sejalan dengan distribusi data observasi.
Namun demikian, terlihat perbedaan pada sebaran ekor distribusi. Metode IDW cenderung menghasilkan distribusi yang lebih menyebar dengan ekor kanan yang lebih panjang, terutama pada nilai prediksi tinggi (>9,5 tahun). Hal ini mengindikasikan bahwa IDW lebih sensitif terhadap pengaruh titik pengamatan bernilai ekstrem di sekitarnya.
Sebaliknya, Universal Kriging menghasilkan distribusi prediksi yang lebih terkonsentrasi dan halus, mencerminkan peran model variogram dan struktur spasial dalam mengendalikan variasi nilai prediksi. Karakteristik ini menunjukkan bahwa UK lebih stabil dalam menghasilkan estimasi Rata-rata Lama Sekolah di seluruh wilayah studi.
Secara keseluruhan, analisis distribusi nilai prediksi memperkuat hasil validasi sebelumnya, bahwa Universal Kriging memberikan estimasi yang lebih konsisten dan representatif secara spasial dibandingkan IDW, sehingga lebih sesuai digunakan dalam pemetaan Rata-rata Lama Sekolah di Provinsi Jawa Timur.
BAB V KESIMPULAN DAN SARAN
5.1 Kesimpulan
Pola Spasial Rata-rata Lama Sekolah (RLS) Hasil uji autokorelasi spasial (Moran’s I, Geary’s C, dan Getis–Ord G) menunjukkan bahwa Rata-rata Lama Sekolah di Provinsi Jawa Timur memiliki pola spasial yang signifikan dan tidak acak. Analisis Getis–Ord Gi* mengidentifikasi keberadaan wilayah hot spot dan cold spot, yang menandakan adanya pengelompokan daerah dengan RLS tinggi dan rendah secara geografis.
Faktor Penentu RLS Berdasarkan model spasial, variabel Tingkat Kemiskinan berpengaruh negatif dan signifikan terhadap RLS, sedangkan Angka Partisipasi Murni (APM) SMA berpengaruh positif dan signifikan. Selain itu, terdapat indikasi spillover spasial pada beberapa variabel, yang menunjukkan bahwa kondisi pendidikan di suatu daerah dipengaruhi pula oleh daerah sekitarnya.
Perbandingan nilai Akaike Information Criterion (AIC) menunjukkan bahwa General Nested Spatial Model (GNS) merupakan model terbaik dengan nilai AIC terendah (AIC = 104,48) dibandingkan model OLS, SDM, SDEM, dan SAC. Selain memiliki kecocokan model terbaik, GNS juga memenuhi seluruh asumsi diagnostik residual, yaitu residual berdistribusi normal, bersifat homoskedastis, dan tidak menunjukkan autokorelasi spasial.
Metode Interpolasi Terbaik Berdasarkan validasi Leave-One-Out Cross-Validation, Universal Kriging (UK) memberikan kinerja yang lebih baik dibandingkan IDW, ditunjukkan oleh nilai RMSE (1,6721) dan MAE (1,4457) yang lebih rendah. Selain itu, UK mampu menghasilkan permukaan spasial yang lebih halus serta menyediakan informasi ketidakpastian melalui peta kriging variance.
Implikasi Kebijakan Wilayah cold spot RLS perlu menjadi prioritas kebijakan pendidikan, khususnya melalui penurunan tingkat kemiskinan, peningkatan akses dan partisipasi pendidikan menengah, serta pemerataan sumber daya pendidikan. Pendekatan kebijakan berbasis wilayah menjadi penting mengingat adanya keterkaitan spasial antar kabupaten/kota.
5.2 Saran Bagi Pemerintah Daerah
Memfokuskan intervensi pendidikan pada wilayah cold spot RLS, terutama melalui program peningkatan partisipasi sekolah menengah dan dukungan sosial bagi rumah tangga miskin.
Mengembangkan kebijakan pendidikan berbasis regional dengan memperhatikan efek spillover antar wilayah, sehingga perencanaan tidak dilakukan secara terisolasi per daerah.
Menambahkan variabel yang mencerminkan kualitas sekolah, seperti rasio fasilitas pendidikan atau kualifikasi tenaga pendidik.
Menggunakan data panel (multi-tahun) untuk menangkap dinamika spasial dan temporal RLS.
Mengembangkan analisis spasio-temporal untuk melihat perubahan pola RLS dan efektivitas kebijakan dari waktu ke waktu.
DAFTAR PUSTAKA
Eryando, T., Nurhayati, S., & Sari, R. K. (2022). Spatial error modeling for regional socioeconomic analysis in Indonesia. Journal of Applied Statistics, 49(12), 3121–3138. https://doi.org/10.1080/02664763.2021.1962263
Mar’ah, S. N., Wibowo, A., & Rahmawati, D. (2025). Spatial dependence and bias correction in regional education models. Journal of Regional Science and Policy, 17(1), 45–63.
Novitasari, D., Pratiwi, D., & Kusuma, H. (2022). Spatial lag and spatial error models in regional development analysis. Jurnal Statistika dan Aplikasinya, 6(2), 85–98. https://doi.org/10.21009/JSA.06207
Wiguna, A. R., Putri, L. S., & Hidayat, R. (2022). Analisis ketergantungan spasial pada indikator pendidikan di Indonesia. Jurnal Ekonomi dan Pembangunan Indonesia, 22(1), 67–82. https://doi.org/10.21002/jepi.v22i1.1356
Anselin, L. (1988). Spatial econometrics: Methods and models. Dordrecht, Netherlands: Kluwer Academic Publishers. https://doi.org/10.1007/978-94-015-7799-1
Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x
Anselin, L., Bera, A. K., Florax, R., & Yoon, M. J. (1996). Simple diagnostic tests for spatial dependence. Regional Science and Urban Economics, 26(1), 77–104. https://doi.org/10.1016/0166-0462(95)02111-6
Badan Pusat Statistik Provinsi Jawa Timur. (2024). Provinsi Jawa Timur dalam angka 2024. Surabaya: BPS Provinsi Jawa Timur. https://jatim.bps.go.id
Badan Pusat Statistik Provinsi Jawa Timur. (2024). Rata-rata lama sekolah menurut kabupaten/kota di Provinsi Jawa Timur. Surabaya: BPS Provinsi Jawa Timur.
Cressie, N. (1993). Statistics for spatial data (Revised ed.). New York, NY: John Wiley & Sons. https://doi.org/10.1002/9780470525212
Getis, A., & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189–206. https://doi.org/10.1111/j.1538-4632.1992.tb00261.x
Getis, A., & Ord, J. K. (1996). Local spatial statistics: An overview. In P. Longley & M. Batty (Eds.), Spatial analysis: Modelling in a GIS environment (pp. 261–277). Cambridge: GeoInformation International.
LeSage, J. P., & Pace, R. K. (2009). Introduction to spatial econometrics. Boca Raton, FL: CRC Press. https://doi.org/10.1201/9781420064254
Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Computers & Geosciences, 30(7), 683–691. https://doi.org/10.1016/j.cageo.2004.03.012
Pebesma, E., & Bivand, R. (2023). Spatial data science: With applications in R. Boca Raton, FL: CRC Press. https://doi.org/10.1201/9780429459016
Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(Suppl.), 234–240. https://doi.org/10.2307/143141
Waller, L. A., & Gotway, C. A. (2004). Applied spatial statistics for public health data. Hoboken, NJ: John Wiley & Sons. https://doi.org/10.1002/0471662682
LAMPIRAN
Lampiran A: Syntax Lengkap
Lampiran B: Output Tambahan
Link Dashboard: https://nafallaffta.shinyapps.io/dashboard/
Link Video Presentasi: https://www.youtube.com/watch?v=s48HI0T2_Ug
Link RPubs: https://rpubs.com/nafallaffta/AnalisisSpasialRatarataLamaSekolahJawaTimur2024_NafallaAfftanurRismawanti