LAPORAN UJIAN AKHIR SEMESTER ANALISIS SPASIAL STATISTICS
Pengaruh Tingkat Partisipasi Angkatan Kerja (TPAK) dan Rata-Rata Lama Sekolah terhadap Tingkat Pengangguran Terbuka (TPT) di Indonesia Tahun 2020-2025 menggunakan Analisis Spasial

Disusun oleh :

Hadeelina Shakira Zalianty (140610230083)

Dosen Pengampu :
Dr. I Gede Nyoman Mindra Jaya,M.Si

PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

DAFTAR ISI

BAB I (PENDAHULUAN)

1.1 Latar belakang

Salah satu dari 17 tujuan global yang dikenal sebagai Sustainable Development Goals (SDGs) adalah mencapai kesempatan kerja penuh dan produktif serta pekerjaan yang layak bagi semua orang pada tahun 2030. Tujuan ini menegaskan bahwa pengangguran merupakan isu global yang serius dan perlu mendapatkan perhatian khusus dari setiap negara. Indonesia menjadi salah satu negara yang menghadapi tantangan tersebut, di mana kondisi pengangguran umumnya diukur menggunakan indikator Tingkat Pengangguran Terbuka (TPT). TPT menggambarkan persentase angkatan kerja yang tidak bekerja namun sedang mencari pekerjaan dan menjadi tolok ukur penting dalam menilai efektivitas kebijakan ketenagakerjaan dan pembangunan ekonomi.

Pada tahun 2020, ketika pandemi COVID-19 melanda, Tingkat Pengangguran Terbuka (TPT) Indonesia tercatat sebesar 7,07 persen, meningkat tajam dibandingkan tahun sebelumnya. Kenaikan ini dipicu oleh penurunan aktivitas ekonomi di berbagai sektor usaha yang berujung pada pemutusan hubungan kerja dan berkurangnya kesempatan kerja. Meskipun dalam beberapa tahun berikutnya terjadi penurunan TPT seiring dengan pemulihan ekonomi nasional, pada tahun 2025 TPT masih tercatat sebesar 4,76 persen, yang menunjukkan bahwa permasalahan ketenagakerjaan di Indonesia belum sepenuhnya pulih. Fenomena pengangguran ini tidak hanya dipengaruhi oleh kondisi ekonomi makro, tetapi juga menunjukkan adanya ketimpangan antarwilayah, yang mengindikasikan keterkaitan spasial antar provinsi di Indonesia.

Berbagai penelitian terdahulu telah mengkaji faktor-faktor yang memengaruhi tingkat pengangguran terbuka, di antaranya Tingkat Partisipasi Angkatan Kerja (TPAK) dan rata-rata lama sekolah. Rata-rata lama sekolah mencerminkan tingkat kualitas sumber daya manusia yang berperan penting dalam menentukan kemampuan tenaga kerja untuk terserap di pasar kerja. Sementara itu, TPAK menggambarkan sejauh mana penduduk usia kerja aktif secara ekonomi. Hasil penelitian menunjukkan bahwa rata-rata lama sekolah dan TPAK memiliki pengaruh yang beragam terhadap TPT. Pada beberapa wilayah, peningkatan lama sekolah dapat menurunkan pengangguran melalui peningkatan kompetensi tenaga kerja, namun pada kondisi tertentu justru dapat meningkatkan pengangguran akibat ketidaksesuaian antara kualifikasi pendidikan dan kebutuhan pasar kerja (skill mismatch). Demikian pula, TPAK umumnya berpengaruh negatif terhadap TPT, meskipun besarnya pengaruh tersebut berbeda antar wilayah.

Dalam menganalisis fenomena pengangguran antarwilayah, pendekatan klasik seperti regresi linier biasa (Ordinary Least Squares / OLS) sering kali belum memadai karena mengasumsikan hubungan yang homogen dan mengabaikan adanya ketergantungan spasial antar wilayah yang berdekatan. Padahal, kondisi ketenagakerjaan suatu provinsi dapat dipengaruhi oleh kondisi provinsi di sekitarnya. Oleh karena itu, diperlukan pendekatan spasial yang mampu menangkap pengaruh lokasi dan interaksi antarwilayah dalam analisis pengangguran.

Beberapa model spasial yang umum digunakan untuk mengakomodasi ketergantungan tersebut antara lain Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM). Model SAR digunakan untuk menangkap pengaruh langsung antar wilayah melalui variabel dependen, model SEM digunakan untuk mengakomodasi ketergantungan spasial yang berasal dari komponen galat, sedangkan model SDM menggabungkan komponen lag spasial pada variabel dependen dan variabel independen. Dengan membandingkan ketiga model tersebut, penelitian ini diharapkan dapat mengidentifikasi pola pengangguran terbuka di Indonesia secara lebih akurat. Pemilihan model terbaik dilakukan berdasarkan nilai Akaike Information Criterion (AIC) dan koefisien determinasi (R²), sehingga model yang terpilih mampu merepresentasikan hubungan antara TPT, rata-rata lama sekolah, dan TPAK secara lebih komprehensif dalam konteks spasial.

Selain model regresi spasial global, penelitian ini juga menggunakan pendekatan Geographically Weighted Regression (GWR) untuk mengidentifikasi heterogenitas spasial, yaitu kondisi di mana pengaruh variabel independen terhadap TPT berbeda-beda pada setiap lokasi. GWR memungkinkan estimasi parameter yang bervariasi secara lokal sehingga dapat menangkap karakteristik unik setiap wilayah. Penentuan bandwidth optimal dalam GWR dilakukan menggunakan metode Cross Validation (CV) dan Akaike Information Criterion Corrected (AICc) untuk memastikan model yang dihasilkan memiliki kecocokan terbaik.

Untuk memvisualisasikan pola spasial dan memprediksi nilai TPT pada wilayah yang tidak tersampel, penelitian ini menggunakan berbagai metode interpolasi spasial, antara lain Nearest Neighbor (NN), Natural Neighbor/TIN, Inverse Distance Weighted (IDW), Trend Surface Analysis, dan metode geostatistik seperti Ordinary Kriging, Universal Kriging, serta Co-Kriging. Metode-metode ini dipilih berdasarkan karakteristik data dan asumsi autokorelasi spasial yang telah diuji sebelumnya. Validasi model interpolasi dilakukan menggunakan teknik cross-validation dengan metrik evaluasi seperti Root Mean Square Error (RMSE), Mean Absolute Error (MAE), dan Mean Error (ME) untuk memastikan akurasi prediksi spasial.

Dengan mengintegrasikan model regresi spasial global (SAR, SEM, SDM), model regresi spasial lokal (GWR), serta metode interpolasi spasial, penelitian ini diharapkan dapat memberikan gambaran yang komprehensif mengenai dinamika pengangguran terbuka di Indonesia periode 2020–2025, baik dari aspek keterkaitan antarwilayah maupun variasi spasial lokal. Hasil penelitian ini juga diharapkan dapat menjadi dasar bagi pemerintah dalam merancang kebijakan ketenagakerjaan yang lebih efektif dan berbasis spasial.

1.2 Identifikasi Masalah

Tingkat Pengangguran Terbuka (TPT) di Indonesia menunjukkan variasi yang cukup tinggi antarprovinsi dan cenderung membentuk pola keterkaitan spasial, di mana wilayah dengan tingkat pengangguran tinggi cenderung berdekatan dengan wilayah yang memiliki kondisi serupa. Hal ini menandakan adanya pengaruh antarwilayah yang tidak dapat dijelaskan hanya melalui pendekatan regresi klasik. Selain itu, faktor-faktor seperti Tingkat Partisipasi Angkatan Kerja (TPAK) dan rata-rata lama sekolah diduga memiliki hubungan yang kompleks terhadap tingkat pengangguran, yang mungkin berbeda pada setiap wilayah (heterogenitas spasial). Diperlukan pendekatan analisis yang mempertimbangkan efek spasial global dan lokal, serta metode interpolasi untuk memperoleh hasil estimasi yang lebih akurat dan komprehensif. Berdasarkan hal tersebut, penelitian ini mengidentifikasi permasalahan sebagai berikut:

1. Apakah terdapat dependensi spasial pada data Tingkat Pengangguran Terbuka (TPT) antarprovinsi di Indonesia periode 2020–2025? Faktor-faktor apa saja yang berpengaruh signifikan terhadap Tingkat Pengangguran Terbuka (TPT) secara global dan lokal?

2. Model spasial apa—Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), atau Spatial Durbin Model (SDM)—yang memberikan hasil estimasi paling akurat dan representatif dalam menjelaskan fenomena pengangguran terbuka di Indonesia?

3. Bagaimana pola variasi spasial lokal dari pengaruh TPAK dan rata-rata lama sekolah terhadap TPT berdasarkan analisis Geographically Weighted Regression (GWR)?

4. Metode interpolasi spasial apa yang paling sesuai untuk memprediksi nilai TPT pada wilayah yang tidak tersampel, dan bagaimana perbandingan kinerja antara metode deterministik (NN, IDW, TIN) dan metode geostatistik (Ordinary Kriging, Universal Kriging, Co-Kriging)?

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk:

1. Menguji adanya dependensi spasial pada data Tingkat Pengangguran Terbuka (TPT) antarprovinsi di Indonesia untuk mengetahui apakah tingkat pengangguran di suatu wilayah dipengaruhi oleh wilayah sekitarnya.

2. Menganalisis pengaruh variabel-variabel ekonomi utama terhadap TPT, yaitu Tingkat Partisipasi Angkatan Kerja (TPAK) dan rata-rata lama sekolah, baik secara global maupun lokal. Membandingkan kinerja model SAR, SEM, dan SDM untuk menentukan model spasial terbaik dalam menjelaskan variasi TPT antarprovinsi berdasarkan nilai Akaike Information Criterion (AIC) dan koefisien determinasi (R²).

3. Mengidentifikasi pola heterogenitas spasial menggunakan Geographically Weighted Regression (GWR) untuk mengetahui variasi pengaruh variabel independen terhadap TPT pada setiap wilayah.

4. Membandingkan berbagai metode interpolasi spasial (Nearest Neighbor, IDW, TIN, Trend Surface, Ordinary Kriging, Universal Kriging, dan Co-Kriging) untuk menentukan metode terbaik dalam memprediksi nilai TPT berdasarkan metrik validasi seperti RMSE, MAE, dan ME.

5. Menghasilkan peta prediksi spasial TPT yang dapat digunakan sebagai acuan dalam perumusan kebijakan ketenagakerjaan berbasis wilayah.

1.4 Manfaat Penelitian

Penelitian ini bermanfaat untuk:

1. Penelitian ini memberikan kontribusi terhadap pengembangan ilmu statistik ekonomi dan geografi kuantitatif, khususnya dalam penerapan analisis ekonometrika spasial, Geographically Weighted Regression (GWR), dan metode interpolasi spasial pada isu ketenagakerjaan.

2. Peta prediksi spasial yang dihasilkan dapat membantu pemangku kebijakan untuk mengidentifikasi wilayah-wilayah prioritas dalam intervensi program ketenagakerjaan.

3. Penelitian ini memberikan contoh penerapan model SAR, SEM, SDM, GWR, serta berbagai metode interpolasi spasial (deterministik dan geostatistik) dalam konteks analisis data spasial untuk fenomena ekonomi regional di Indonesia, yang dapat direplikasi untuk studi serupa.

1.5 Batasan Penelitian

Agar penelitian lebih terarah dan fokus, maka perlu dibatasi pada beberapa hal:

1. Unit analisis yang digunakan adalah 34 provinsi di Indonesia dengan periode pengamatan tahun 2020–2025 (data panel enam tahun).

2. Variabel yang digunakan terdiri atas variabel dependen, yaitu Tingkat Pengangguran Terbuka (TPT) pada termin 1 (Februari) dan variabel independen, yaitu Tingkat Partisipasi Angkatan Kerja (TPAK) pada termin 1 (Februari) dan rata-rata lama sekolah.

3. Data yang digunakan merupakan data sekunder yang diperoleh dari Badan Pusat Statistik (BPS).

4. Analisis dilakukan menggunakan pendekatan Spatial Panel Regression dengan tiga model utama, yaitu Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM).

5. Analisis heterogenitas spasial dilakukan menggunakan Geographically Weighted Regression (GWR) dengan bandwidth optimal yang dipilih berdasarkan metode Cross Validation (CV) atau AICc.

6. Metode interpolasi spasial yang digunakan meliputi Nearest Neighbor (NN), Natural Neighbor/TIN, Inverse Distance Weighted (IDW), Trend Surface Analysis, Ordinary Kriging, Universal Kriging, dan Co-Kriging.

7. Pemilihan model terbaik didasarkan pada kriteria AIC terkecil, koefisien determinasi (R²) terbesar, serta metrik validasi spasial seperti RMSE, MAE, dan ME.

8. Penelitian ini tidak membahas aspek kebijakan makro ekonomi secara mendalam, melainkan fokus pada analisis statistik spasial dan pemodelan ekonometrika spasial.

BAB II (TINJAUAN PUSTAKA)

2.1 Tingkat Pengangguran Terbuka (TPT)

Tingkat pengangguran terbuka (TPT) merupakan indikator penting dalam mengukur kondisi pasar tenaga kerja di suatu negara atau daerah. TPT mencerminkan persentase angkatan kerja yang aktif mencari pekerjaan namun belum mendapatkan pekerjaan yang layak. Topik ini dipilih karena pengangguran terbuka merupakan masalah ekonomi yang mendasar, yang dapat mempengaruhi stabilitas sosial dan pertumbuhan ekonomi nasional. Misalnya, peningkatan TPT dapat menyebabkan berkurangnya produktivitas dan pendapatan masyarakat, serta menimbulkan masalah keTPAKan dan sosial lainnya.

Urgensi penelitian tentang TPT sangat tinggi mengingat adanya dinamika ketenagakerjaan yang dipengaruhi oleh faktor-faktor seperti perubahan struktur ekonomi, kesenjangan kualifikasi tenaga kerja dengan kebutuhan industri, dan dampak pandemi Covid-19 yang memperburuk angka pengangguran di Indonesia. Data dari Badan Pusat Statistik menunjukkan peningkatan signifikan TPT selama pandemi, yang menuntut penanganan cepat dan tepat oleh pemerintah agar lapangan kerja mampu menyerap tenaga kerja secara efektif. Dibutuhkan pemahaman mendalam tentang TPT untuk mendukung kebijakan ketenagakerjaan yang efektif serta mendorong pembangunan ekonomi yang berkelanjutan.

Selain itu, Lama Sekolah dan persentase TPAK juga memiliki keterkaitan erat dengan dinamika pengangguran terbuka. Lama.Sekolah menggambarkan seberapa besar proporsi TPAK usia kerja yang aktif berpartisipasi dalam pasar tenaga kerja, baik dengan bekerja maupun mencari pekerjaan. Semakin tinggi Lama.Sekolah, maka semakin besar potensi tekanan terhadap pasar kerja apabila tidak diimbangi dengan ketersediaan lapangan pekerjaan yang memadai. Di sisi lain, tingginya persentase TPAK TPAK sering kali menjadi cerminan keterbatasan akses terhadap pendidikan, keterampilan, dan peluang kerja yang layak, yang pada akhirnya berkontribusi terhadap peningkatan TPT. Analisis hubungan antara TPT, Lama.Sekolah, dan TPAK TPAK menjadi penting untuk memahami pola ketenagakerjaan dan merumuskan kebijakan ekonomi yang lebih terarah.

2.2 Teori Dasar Spatial Dependence

Menurut Wooldridge (2012), salah satu asumsi utama dalam model regresi klasik adalah independensi antar observasi. Namun, dalam konteks geografis dan ekonomi regional, asumsi ini sering kali tidak terpenuhi karena adanya keterkaitan antar wilayah yang saling berdekatan. Fenomena tersebut dikenal sebagai spatial dependence atau ketergantungan spasial, yaitu kondisi ketika nilai suatu variabel di satu lokasi dipengaruhi oleh nilai variabel yang sama di lokasi lain. Prinsip ini sejalan dengan Tobler’s First Law of Geography (1970) yang menyatakan bahwa “segala sesuatu berhubungan dengan yang lain, tetapi hal-hal yang berdekatan memiliki hubungan lebih kuat daripada yang jauh.” Dengan demikian, analisis spasial menjadi penting untuk menangkap adanya pola interaksi antar wilayah yang tidak dapat dijelaskan hanya oleh variabel lokal semata.

Secara umum, spatial dependence dapat diukur melalui indikator autokorelasi spasial seperti Moran’s I, Geary’s C,Getis-Ord G, dan Local Morasn I (LISA) yang digunakan untuk mendeteksi apakah nilai suatu variabel menunjukkan pola pengelompokan atau penyebaran secara geografis. Ukuran tersebut biasanya menggunakan matriks bobot spasial (spatial weights matrix, W) yang mendefinisikan hubungan kedekatan antar wilayah, baik berdasarkan administratif maupun jarak geografis (distance-based). Ketika hasil uji menunjukkan adanya autokorelasi spasial yang signifikan, maka pendekatan regresi konvensional (OLS) tidak lagi memadai, karena menghasilkan estimasi yang bias dan tidak efisien. Dalam hal ini, dibutuhkan model spasial yang mampu mengakomodasi ketergantungan antar wilayah, seperti Spatial Autoregressive Model (SAR) yang menangkap efek spasial pada variabel dependen, dan Spatial Error Model (SEM) yang memperhitungkan efek spasial melalui komponen error.

Dalam konteks penelitian ini, spatial dependence sangat relevan untuk menganalisis pola Tingkat Pengangguran Terbuka (TPT) antar provinsi di Indonesia selama periode 2020–2025. Interaksi ekonomi regional, mobilitas tenaga kerja, serta kesamaan karakteristik struktural antar provinsi menyebabkan TPT di satu wilayah berpotensi memengaruhi wilayah sekitarnya. Misalnya, peningkatan TPT di Jawa Barat dapat berdampak pada provinsi tetangga seperti Banten atau DKI Jakarta karena keterkaitan pasar tenaga kerja dan aktivitas industri lintas daerah. Oleh sebab itu, pemahaman terhadap teori ketergantungan spasial menjadi krusial agar model yang digunakan mampu menggambarkan realitas ekonomi antarwilayah secara akurat serta menghasilkan rekomendasi kebijakan ketenagakerjaan yang lebih tepat sasaran.

2.3 Autokorelasi Spasial

Autokorelasi Spasial merupakan konsep yang menggambarkan sejauh mana nilai suatu variabel di suatu lokasi berhubungan dengan nilai variabel yang sama di lokasi lain yang berdekatan secara geografis. Secara sederhana, fenomena ini menunjukkan adanya pola keterkaitan antarwilayah, di mana wilayah yang berdekatan cenderung memiliki karakteristik yang mirip (autokorelasi positif), atau justru sangat berbeda (autokorelasi negatif). Ketika nilai suatu variabel di satu daerah memengaruhi nilai di daerah lain, maka data tersebut tidak lagi bersifat independen secara spasial.

Secara matematis, spatial autocorrelation dapat dijelaskan melalui kovarians antara dua observasi yi​ dan yj​ yang berjarak dij, dengan Cov(yi​, yj​ ) ≠ 0. Kondisi ini melanggar asumsi dasar regresi OLS yang mengharuskan komponen error antar observasi bersifat bebas (independent and identically distributed). Jika ketergantungan spasial diabaikan, maka hasil estimasi akan menjadi bias dan tidak efisien.

Untuk mendeteksi adanya autokorelasi spasial, beberapa ukuran statistik digunakan, seperti Moran’s I, Geary’s C, dan Getis-Ord G. Nilai positif dari Moran’s I menunjukkan adanya pengelompokan wilayah dengan nilai variabel yang serupa, sedangkan nilai negatif mengindikasikan penyebaran yang berlawanan arah. Temuan adanya autokorelasi spasial yang signifikan menunjukkan perlunya penerapan model ekonometrika spasial (spatial econometric models) yang dapat mengakomodasi hubungan antarwilayah tersebut.

2.3.1 Moran’s I

Uji Moran’s I Global digunakan untuk mendeteksi keberadaan autokorelasi spasial dalam distribusi TPT antarprovinsi. Rumusnya sebagai berikut:

\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij} (x_i - \bar{x})(x_j - \bar{x})} {\sum_i (x_i - \bar{x})^2} \]

Keterangan:
- \(n\) = jumlah unit observasi
- \(x_i, x_j\) = nilai variabel di lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata nilai variabel
- \(w_{ij}\) = elemen matriks pembobot spasial

Kriteria Pengujian:
- \(I > 0\): autokorelasi spasial positif (wilayah dengan nilai serupa saling berdekatan)
- \(I < 0\): autokorelasi spasial negatif (wilayah dengan nilai berbeda berdekatan)
- \(I \approx 0\): tidak terdapat autokorelasi spasial (pola acak)

2.3.2 Geary’s C

Statistik Geary’s C juga digunakan untuk mengukur autokorelasi spasial, namun berbeda dengan Moran’s I yang melihat hubungan kovarian antar lokasi, Geary’s C lebih sensitif terhadap perbedaan lokal antar pasangan lokasi. Statistik ini diperkenalkan oleh R. C. Geary pada tahun 1954. Rumusnya sebagai berikut:

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

Keterangan:
- \(n\) = jumlah unit observasi
- \(x_i, x_j\) = nilai variabel di lokasi \(i\) dan \(j\)
- \(\bar{x}\) = rata-rata nilai variabel
- \(w_{ij}\) = elemen matriks pembobot spasial

Kriteria Pengujian:
- \(C = 1\): tidak ada autokorelasi spasial
- \(C < 1\): terdapat autokorelasi spasial positif
- \(C > 1\): terdapat autokorelasi spasial negatif

2.3.3 Getis-Ord G

Statistik Getis-Ord G (atau sering disebut Getis-Ord General G) dikembangkan oleh Getis dan Ord (1992) untuk mengukur autokorelasi spasial global maupun lokal dengan fokus pada pengelompokan nilai tinggi (hot spots) dan nilai rendah (cold spots).

\[ G = \frac{\sum_i \sum_j w_{ij} x_i x_j} {\sum_i \sum_j x_i x_j} \]

Keterangan:
- \(x_i, x_j\) = nilai variabel pada lokasi \(i\) dan \(j\)
- \(w_{ij}\) = bobot spasial antara lokasi \(i\) dan \(j\)

Kriteria Pengujian:
- Nilai \(G\) tinggi menunjukkan hot spot (pengelompokan nilai besar).
- Nilai \(G\) rendah menunjukkan cold spot (pengelompokan nilai kecil).
- Nilai \(G\) mendekati rata-rata global menunjukkan pola acak.

Terdapat juga versi lokal dari statistik ini, yaitu \(Getis-Ord Gi^*\), yang digunakan untuk mendeteksi lokasi-lokasi spesifik yang menjadi pusat hot spot atau cold spot.

2.3.4 Local Moran’s I (LISA)

Menurut Anselin (1995) dalam Saputro et al. (2018), LISA didefinisikan sebagai statistik yang memenuhi dua kriteria utama. Pertama, nilai LISA untuk setiap wilayah dapat digunakan untuk menunjukkan adanya pengelompokan hubungan spasial yang signifikan darinilai-nilai serupa di sekitar wilayah tersebut. Kedua, jumlah dari nilai LISA untuk seluruh wilayah sebanding dengan nilai Indeks Moran. Kriteria ini menjadikan LISA alat yang efektif untuk mengidentifikasi pola spasial yang terorganisir secara lokal dalam suatu area. LISA untuk setiap wilayah ditulis sebagai berikut.

Nilai LISA (\(L_i\)) untuk setiap lokasi \(i\) dihitung sebagai berikut:

\[ \begin{align} L_i &= \frac{z_i}{m_2} \sum_{j} w_{ij} z_j \quad & \\ z_i &= (x_i - \bar{x}) \quad & \\ m_2 &= \frac{\sum_{i=1}^{n} z_i^2}{n} = \frac{\sum_{i=1}^{n} (x_i - \bar{x})^2}{n} \quad &\\ z_j &= (x_j - \bar{x}) \quad &\end{align} \]

Keterangan Simbol:

• \(L_i\): Nilai LISA untuk lokasi \(i\).
• \(z_i\) dan \(z_j\): Deviasi (selisih) dari rata-rata nilai variabel untuk lokasi \(i\) dan \(j\).
• \(x_i\): Nilai variabel \(x\) pada lokasi \(i\) (dengan \(i = 1, 2, 3, \dots, n\)).
• \(x_j\): Nilai variabel \(x\) pada lokasi \(j\) (dengan \(j = 1, 2, 3, \dots, n\)).
• \(\bar{x}\): Rata-rata dari nilai variabel \(x\) di seluruh lokasi.
• \(w_{ij}\): Elemen pembobot spasial yang mengukur hubungan kedekatan antara lokasi \(i\) dan \(j\).
• \(n\): Banyaknya kejadian atau lokasi/area dalam analisis.
• \(m_2\): Standar deviasi kuadrat (varians populasi) dari variabel \(x\).

2.4 Estimasi Model Spasial

Model Spatial Econometrics dikembangkan untuk mengatasi permasalahan spatial dependence dalam data yang memiliki dimensi geografis. Menurut Anselin (1988) dan LeSage & Pace (2009), model ini digunakan untuk memperbaiki hasil estimasi regresi yang bias akibat adanya hubungan spasial antar unit observasi. Dua bentuk model yang paling umum digunakan adalah Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM).

2.4.1 Spatial Autoregressive Model (SAR)

Spatial Autoregressive Model (SAR) digunakan ketika nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah lain yang berdekatan. Artinya, terdapat efek limpahan (spillover effect) antarwilayah yang menyebabkan suatu daerah tidak dapat dipelajari secara terpisah. Model SAR dapat ditulis sebagai berikut:

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

Keterangan:

• \(y\) = vektor variabel dependen (misalnya Tingkat Pengangguran Terbuka / TPT)
• \(X\) = matriks variabel independen (misalnya Lama.Sekolah dan persentase TPAK TPAK)
• \(\beta\) = vektor koefisien regresi
• \(W\) = matriks pembobot spasial (spatial weights matrix)
• \(\rho\) = koefisien autokorelasi spasial (spatial autoregressive coefficient)
• \(\varepsilon\) = vektor error, dengan asumsi \(\varepsilon \sim N(0, \sigma^2 I)\)

Koefisien \(\rho\) menunjukkan sejauh mana pengaruh spasial terjadi antarwilayah. Jika \(\rho\) signifikan, maka variabel dependen di suatu wilayah dipengaruhi oleh nilai di wilayah tetangga. Dalam konteks penelitian TPT antarprovinsi, model SAR relevan apabila tingkat pengangguran di satu provinsi dipengaruhi oleh kondisi pengangguran di provinsi sekitar, misalnya karena mobilitas tenaga kerja, migrasi, atau keterkaitan pasar kerja regional.

2.4.2 Spatial Error Model (SEM)

Berbeda dengan SAR, Spatial Error Model (SEM) digunakan ketika ketergantungan spasial tidak terjadi pada variabel dependen, tetapi pada komponen error. Artinya, pengaruh antarwilayah muncul akibat faktor yang tidak terobservasi namun saling berkaitan secara spasial, seperti kebijakan daerah, infrastruktur, atau kondisi sosial ekonomi yang mirip. Model SEM dinyatakan sebagai berikut:

\[ y = X\beta + (I - \lambda W)^{-1} \varepsilon \]

Keterangan:

• \(y\) = vektor variabel dependen (TPT)
• \(X\) = matriks variabel independen (Lama.Sekolah dan persentase TPAK TPAK)
• \(\beta\) = vektor koefisien regresi
• \(W\) = matriks pembobot spasial
• \(\lambda\) = koefisien autokorelasi spasial pada komponen error
• \(\varepsilon\) = error acak yang berdistribusi normal \(N(0, \sigma^2 I)\)

Jika \(\lambda\) signifikan, maka terdapat hubungan spasial dalam error term, yang menandakan bahwa variabel-variabel yang tidak dimasukkan ke dalam model berhubungan antarwilayah. Dalam konteks TPT 2020–2025, SEM cocok digunakan ketika faktor-faktor eksternal seperti kebijakan ketenagakerjaan, investasi, atau kondisi sosial ekonomi antarprovinsi memiliki keterkaitan yang menyebabkan error model saling berhubungan.

2.5 Geographically Weighted Regression (GWR)

Geographically Weighted Regression(GWR) merupakan metode statistika yang digunakan untuk menganalisis heterogenitas spasial (Fotheringham dkk., 2002). Heterogenitas spasial adalah kondisi dimana satu variabel independen yang sama memberikan respon yang tidak sama pada lokasi yang berbeda dalam satu wilayah penelitian. Pada model GWR, variabel dependen ditaksir dengan variabel independen yang setiap koefisien regresinya tergantung pada lokasi dimana data tersebut diamati (Caraka & Yasin, 2017)

Model regresi GWR untuk observasi ke-\(i\) dinyatakan sebagai berikut:

\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) x_{ik} + \varepsilon_i, \quad i = 1, 2, \dots, n \]

Keterangan Simbol:

• \(y_i\): Nilai observasi variabel dependen ke-\(i\).
• \((u_i, v_i)\): Koordinat geografis lokasi observasi ke-\(i\).
• \(\beta_0(u_i, v_i)\): Konstanta (intercept) pada lokasi ke-\(i\). Karena bergantung pada \((u_i, v_i)\), intercept bervariasi secara spasial.
• \(\beta_k(u_i, v_i)\): Koefisien regresi ke-\(k\) untuk variabel independen ke-\(k\) pada lokasi ke-\(i\). Koefisien ini bervariasi secara spasial.
• \(x_{ik}\): Nilai observasi variabel independen ke-\(k\) pada lokasi ke-\(i\).
• \(p\): Jumlah total variabel independen (prediktor) dalam model.
• \(\varepsilon_i\): Error (galat) pengamatan ke-\(i\), yang diasumsikan:
  • Identik, independen.
  • Berdistribusi normal dengan rata-rata (mean) nol (\(0\)) dan varian konstan \(\sigma^2\) (i.e., \(\varepsilon_i \sim N(0, \sigma^2)\)).
• \(n\): Jumlah total observasi.

2.6 Bandwith

Bandwidth secara teoritis adalah lingkaran dengan radius ℎ dari titik pusat lokasi pengamatan yang digunakan sebagai dasar penentuan pembobot untuk setiap lokasi pengamatan. Titik-titik lokasi pengamatan yang berada dalam lingkaran dengan radius ℎ dari titik lokasi pengamatan ke-𝑖 akan berpengaruh dalam menduga parameter di titik lokasi pengamatan ke-𝑖.

2.6.1 Bandwidth Cross Validation (CV) dan Akaike Information Criterion Corrected (AIcC)

Proses Cross Validation (CV) apabila diterapkan ke dalam konsep bandwidth, maka kelompok data testing diambil dengan cara menghilangkan nilai titik data ke-i dari proses pendugaan 𝑦𝑖 untuk dibandingkan dengan kelompok data training (data riil atau data keseluruhan) yang dilakukan dengan cara iterasi sampai diperoleh nilai CV minimum.

Fungsi CV yang diminimumkan untuk memilih bandwidth (\(h\)) adalah:

\[ CV = \sum_{i=1}^{n} [y_i - \hat{y}_{\ne i}(h)]^2 \] Metode pemilihan bandwidth dengan AICc dilakukan secara iterasi, yaitu dengan cara memasukkan berbagai nilai bandwidth ke dalam rumus AICc sampai diperoleh nilai AICc minimum. Nilai bandwidth dikatakan optimum ketika nilai AICc yang dihasilkan minimum.

Secara matematis, AICc dapat dituliskan sebagai berikut: \[ \text{AICc} = 2n \ln(\hat{\sigma}) + n \ln(2\pi) + n \left\{ \frac{n+\text{tr}(\mathbf{L})}{n-2-\text{tr}(\mathbf{L})} \right\} \] dengan \(\hat{\sigma}\) adalah: \[ \hat{\sigma} = \sqrt{\frac{y^T(\mathbf{I}-\mathbf{L})^T (\mathbf{I}-\mathbf{L}) y}{n}} \] Keterangan Simbol AICc dan CV:

• \(\hat{y}_{\ne i}(h)\): Nilai penduga \(y_i\) dengan pengamatan di lokasi ke-\(i\) dihilangkan dari proses pendugaan.
• \(y_i\): Variabel dependen pada pengamatan ke-\(i\).
• \(h\): Bandwidth (parameter pembobot spasial).

2.7 Metode Estimasi GWR

2.7.1 Fungsi Kernel

Model GWR menggunakan fungsi Kernel sebagai pembobot spasial, salah satunya adalah Fixed Kernel yang memiliki bandwidth (\(b\)) yang sama pada setiap titik lokasi pengamatan.

1. Fungsi Kernel Gaussian

\[ \w_i = e \cdot \left[ - \frac{1}{2} \left(\frac{d_i}{b}\right)^2 \right] \quad \]

2. Fungsi Kernel Bi-square

\[ w_i = \begin{cases} \left[ 1-\left(\frac{d_i}{b}\right)^2 \right]^2 & \text{jika } d_i < b \\ 0 & \text{jika } d_i \ge b \end{cases} \quad \]

Keterangan Simbol Fungsi Kernel:

• \(w_i\): Bobot spasial pada lokasi \(i\).
• \(d_i\): Jarak antara titik pengamatan ke-\(i\) dan titik di mana koefisien sedang dihitung.
• \(b\): Bandwidth (parameter pembobot spasial yang sama untuk setiap lokasi).

2.8 Uji Asumsi GWR

2.8.1. Normalitas Residu

• Asumsi: Residu model (\(\varepsilon_i\)) harus berdistribusi normal.
• Uji: Uji Jarque-Bera.
• Rumus (Statistik Uji Jarque-Bera - JB):

\[ JB = \frac{n}{6} \left( S^2 + \frac{(K-3)^2}{4} \right) \] Di mana \(S\) adalah skewness (kemencengan) dan \(K\) adalah kurtosis (keruncingan) dari residu, dan \(n\) adalah jumlah observasi.

2.8.2. Kehomogenan Varian (Homoskedastisitas)

• Asumsi: Varian dari residu harus konstan (\(\sigma^2\)). GWR umumnya sudah mengatasi heteroskedastisitas spasial (variasi koefisien), tetapi perlu dipastikan bahwa varian residu tidak bergantung pada nilai duga.
• Uji: Uji Breusch-Pagan.
• Rumus (Statistik Uji Breusch-Pagan - LM):

\[ LM = n \cdot R^2 \] Di mana \(n\) adalah jumlah observasi dan \(R^2\) adalah koefisien determinasi dari regresi residu kuadrat terhadap variabel independen (atau nilai duga).

2.8.3. Non-Autokorelasi Spasial Residu

• Asumsi: Residu model harus independen secara spasial. Ini diuji untuk memastikan tidak ada ketergantungan spasial yang tersisa setelah GWR diterapkan.
• Uji: Uji Moran’s I pada residu GWR.
• Rumus (Statistik Uji Moran’s I):

\[ I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}} \frac{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij} (e_i - \bar{e})(e_j - \bar{e})}{\sum_{i=1}^{n} (e_i - \bar{e})^2} \] Di mana \(e_i\) dan \(e_j\) adalah residu pada lokasi \(i\) dan \(j\), \(\bar{e}\) adalah rata-rata residu (idealinya nol), dan \(w_{ij}\) adalah elemen pembobot spasial.

2.8.4. Non-Multikolinearitas

• Asumsi: Tidak boleh ada hubungan linier yang sempurna atau tinggi antar variabel independen (\(x_{ik}\)).
• Uji: Variance Inflation Factor (VIF) lokal.
• Rumus (Variance Inflation Factor - VIF):

\[ \text{VIF}_k = \frac{1}{1 - R_k^2} \] Di mana \(R_k^2\) adalah koefisien determinasi dari regresi variabel independen ke-\(k\) terhadap semua variabel independen lainnya. Dalam GWR, VIF ini dihitung secara lokal.

2.9 Spatial Interpolation

Interpolasi spasial merupakan sebuah metode untuk memprediksi atau mengestimasi nilai suatu variabel di lokasi-lokasi yang belum diamati, berdasarkan data nilai yang telah teramati dari lokasi-lokasi di sekitarnya.Tujuan utama dari prosedur ini adalah untuk menghasilkan permukaan kontinu (surface) dari data yang awalnya berbentuk titik diskret. Menurut Demers (2000), interpolasi spasial dapat diklasifikasikan menjadi global and local interpolation.

2.10 Asumsi dalam Spatial Interpolation

1. Autokorelasi Spasial (Spatial Autocorrelation)

• Penjelasan Konseptual: Interpolasi bekerja berdasarkan Hukum Pertama Geografi (Tobler), yang menyatakan bahwa semua objek spasial saling berhubungan, tetapi objek yang berdekatan lebih berhubungan daripada objek yang berjauhan.
• Pengaruh: Autokorelasi spasial memengaruhi pembobot yang diberikan ke titik-titik tetangga dalam estimasi.

2. Stasioneritas dan Tren (Stationarity and Trend)

• Stasioneritas: Asumsi fundamental di mana sifat-sifat statistik data (seperti rata-rata dan varians) diasumsikan konstan di seluruh wilayah studi, atau setidaknya konstan dalam jendela lokal (stasioneritas intrinsik). Jika data tidak stasioner (ada tren global), tren tersebut harus dihilangkan sebelum interpolasi geostatistik (seperti Kriging) dilakukan.
• Tren: Variasi sistematis dalam nilai variabel yang terkait dengan koordinat lokasi.

3. Skala dan Anisotropi (Scale and Anisotropy)

• Skala: Interpolasi sangat dipengaruhi oleh skala data atau ukuran area studi yang dianalisis.
• Anisotropi: Kondisi di mana variasi spasial (korelasi spasial) dari variabel berbeda tergantung pada arah jarak, bukan hanya pada besar jarak. Metode interpolasi geostatistik (Kriging) harus memperhitungkan anisotropi.

4. Kerapatan dan Sebaran Sampel (Sample Density and Distribution)

• Kerapatan (Density): Kualitas permukaan hasil interpolasi sangat bergantung pada seberapa padat (jumlah sampel per unit area) dan seberapa merata sebaran sampel di wilayah studi. Area dengan sampel yang jarang cenderung menghasilkan estimasi yang kurang akurat.

5. Validasi Model (Model Validation)

• Penjelasan: Semua hasil interpolasi harus divalidasi untuk mengukur keakuratan model, seringkali menggunakan teknik cross-validation (validasi silang) atau membandingkan hasil estimasi dengan data yang diuji (nilai yang diketahui).

Cross-Validation (Validasi Silang) RMSE (Root Mean Square Error): Mengukur akurasi keseluruhan model, yaitu perbedaan rata-rata antara nilai observasi dan nilai prediksi. Nilai ideal adalah sekecil mungkin. \[ \text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n} (\hat{Z}(s_i) - Z(s_i))^2} \] Di mana \(\hat{Z}(s_i)\) adalah nilai prediksi di lokasi \(s_i\) dan \(Z(s_i)\) adalah nilai observasi di lokasi \(s_i\).

MAE (Mean Absolute Error): Mengukur rata-rata deviasi absolut, yang menunjukkan stabilitas prediksi. Nilai ideal adalah sekecil mungkin.

ME (Mean Error atau Bias): Menunjukkan kecenderungan over-prediction atau under-prediction (bias rata-rata). Nilai ideal adalah mendekati nol.

Validasi model juga dapat dilakukan secara visual untuk memeriksa bias spasial yang mungkin tersisa:

• Scatter plot: Visualisasi hubungan antara nilai observasi dan nilai prediksi.
• Peta residual: Digunakan untuk melihat adanya bias spasial yang tersisa dalam kesalahan estimasi.
• Variogram residual: Digunakan untuk memastikan bahwa tidak ada struktur spasial (autokorelasi spasial) yang tersisa di dalam residual, yang menandakan bahwa model telah menangkap sepenuhnya ketergantungan spasial.

2.11 Metode Interpolasi

2.11.1 Nearest Neighbor (NN)

NN merupakan salah satu metode dalam memprediksi nilai pada daerah yang tidak tersampel dengan menghubungkan titik-titik yang berdekatan atau disebut dengan replikasi pixel (Hadi and Tombul, 2018)

\[ \hat{Z}(s_0) = Z(s_{(1)}) \] ### 2.11.2 Natural Neighbor Interpolation (TIN)

Metode Natural Neighbor Interpolation adalah teknik interpolasi non-parametrik yang didasarkan sepenuhnya pada struktur spasial titik-titik data teramati. Metode ini menggunakan Diagram Voronoi (atau Poligon Thiessen) yang dihasilkan dari TIN (Triangulated Irregular Network) untuk menentukan bobot yang diberikan kepada setiap titik sampel terdekat.

1. Triangulated Irregular Network (TIN)

TIN adalah representasi permukaan yang dibentuk oleh sekumpulan simpul (titik sampel teramati) yang dihubungkan untuk membentuk jaringan segitiga non-tumpang tindih (non-overlapping). TIN menjadi dasar untuk mendefinisikan Poligon Thiessen/Voronoi, yang pada gilirannya digunakan untuk menghitung bobot Natural Neighbor.

2. Diagram Voronoi (Thiessen Polygons)

Diagram Voronoi membagi ruang menjadi wilayah-wilayah yang berdekatan. Setiap wilayah (poligon) berasosiasi dengan satu titik sampel, dan setiap lokasi di dalam poligon tersebut jaraknya lebih dekat ke titik sampelnya daripada ke titik sampel lainnya.

3. Perhitungan Bobot Natural Neighbor Bobot dalam metode Natural Neighbor ditentukan oleh perubahan luas wilayah Voronoi ketika titik estimasi (\(x_0\)) ditambahkan ke kumpulan titik data teramati.

• Estimasi Nilai: Nilai estimasi variabel \(Z\) pada lokasi yang tidak teramati \(x_0\), \(\hat{Z}(x_0)\), dihitung sebagai rata-rata tertimbang (weighted average) dari nilai-nilai \(Z\) dari tetangga alaminya (titik-titik sampel yang Poligon Voronoi-nya berbatasan langsung dengan Poligon Voronoi baru untuk \(x_0\)): \[ \hat{Z}(x_0) = \sum_{i \in N(x_0)} w_i(x_0) Z(x_i) \] Di mana \(N(x_0)\) adalah himpunan tetangga alami untuk \(x_0\).

2.11.3 Interpolasi Inverse Distance Weighted (IDW)

IDW adalah salah satu metode interpolasi spasial yang memiliki asumsi bahwa setiap titik input mempunyai pengaruh yang bersipat lokal yang berkurang terhadap jarak. Nilai power pada interpolasi IDW menentukan pengaruh terhadap titik – titik input, dimana pengaruh akan lebih besar pada titik – titik yang lebih dekat sehingga menghasilkan permukaan yang lebih detail.

Bobot spasial (\(w_i\)) pada lokasi estimasi \(s_0\) dihitung sebagai:

\[ w_i(s_0) \propto d(s_0, s_i)^{-p}, \quad p > 0 \]

Di mana \(d(\cdot, \cdot)\) umumnya adalah jarak Euclidean.

\[ \hat{Z}(s_0) = \frac{\sum_{i} w_i(s_0) Z(s_i)}{\sum_{i} w_i(s_0)} \]

• \(p \to 0\): Semua bobot hampir sama, sehingga estimasi mendekati rata-rata biasa dari semua titik sampel.
• \(p = 1\): Penurunan bobot linier terhadap jarak (dalam skala log–log).
• \(p = 2\): Penurunan kuadratik; pilihan populer untuk fenomena yang meredup cepat.
• \(p \to \infty\): IDW mendekati interpolasi Nearest Neighbor, karena hanya tetangga terdekat yang memiliki bobot dominan.

2.11.4 Kriging

metode interpolasi geostatistik yang lebih canggih dan berbasis model. Bobot DIHITUNG erdasarkan struktur spasial data, yang diukur melalui variogram (atau semivariogram). Bobot tidak hanya ditentukan oleh jarak, tetapi juga oleh autokorelasi spasial antar titik dan pengaturan spasial dari semua titik tetangga relatif terhadap titik estimasi.

1. Bidang Acak, Stasioneritas, dan Isotropi

Bidang Acak (Random Field)

Bidang acak adalah sekumpulan peubah acak yang terindeks oleh lokasi spasial: \[ \{ Z(\mathbf{s}) : \mathbf{s} \in D \subset \mathbb{R}^2 \} \]

---

Stasioneritas

1. Stasioneritas Orde-1
   Rataan konstan: \[ E[Z(\mathbf{s})] = \mu \]

2. Stasioneritas Orde-2
   Kovarians hanya bergantung pada vektor jarak \(\mathbf{h}\): \[ Cov(Z(\mathbf{s}), Z(\mathbf{s+h})) = C(\mathbf{h}) \]

Isotropi Bidang acak dikatakan isotropik jika kovarians hanya bergantung pada jarak, bukan arah: \[ C(\mathbf{h}) = C(h), \quad h = \|\mathbf{h}\| \]

Jika kovarians bergantung pada arah, maka disebut anisotropi.

2. Isotropi vs Anisotropi

• Isotropi: struktur spasial sama ke semua arah
• Anisotropi: variasi spasial berbeda tergantung arah

Anisotropi sering muncul pada fenomena geografi seperti aliran air atau angin.

Variogram dan Kovarians

Kovarians didefinisikan sebagai: \[ C(h) = Cov(Z(\mathbf{s}), Z(\mathbf{s+h})) \]

Variogram: \[ \gamma(h) = \frac{1}{2} Var(Z(\mathbf{s}) - Z(\mathbf{s+h})) \]

Hubungan antara variogram dan kovarians: \[ \gamma(h) = C(0) - C(h) \]

3. Komponen Variogram

Variogram terdiri dari: 1. Nugget (\(C_0\)): variasi pada jarak sangat kecil (error pengukuran) 2. Sill (\(C_0 + C\)): nilai variogram maksimum 3. Range (\(a\)): jarak di mana variogram mencapai sill

4. Model Variogram Teoretis

Beberapa model variogram yang umum:

1. Spherical \[ \gamma(h)= \begin{cases} C_0 + C\left[\frac{3h}{2a}-\frac{h^3}{2a^3}\right], & h \le a \\ C_0 + C, & h>a \end{cases} \]

2. Exponential \[ \gamma(h)=C_0 + C\left(1 - e^{-h/a}\right) \]

3. Gaussian \[ \gamma(h)=C_0 + C\left(1 - e^{-(h^2/a^2)}\right) \]

5. Detail Variogram

Variogram empiris dihitung dari data sebagai: \[ \hat{\gamma}(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)} [Z(\mathbf{s}_i)-Z(\mathbf{s}_i+h)]^2 \]

dengan \(N(h)\) adalah jumlah pasangan pengamatan pada jarak \(h\).

Konsep Dasar Variogram dan Semivariogram

Semivariogram menggambarkan tingkat kemiripan nilai antar lokasi. Semakin kecil \(\gamma(h)\), semakin kuat korelasi spasial antar titik.

6. Nugget

Nugget (\(C_0\)) merepresentasikan: - Error pengukuran - Variasi mikro skala yang tidak teramati

Nilai nugget ditunjukkan oleh intercept variogram pada \(h=0\).

7. Sill

Sill adalah nilai plateau variogram: \[ Sill = C_0 + C \]

Menunjukkan varians total dari proses spasial.

8. Range

Range (\(a\)) adalah jarak maksimum pengaruh spasial. Untuk \(h > a\), pengamatan tidak lagi berkorelasi.

9. Keluarga Kriging: SK, OK, UK

1. Simple Kriging (SK)
   Rataan diketahui: \[ \sum_{i=1}^n \lambda_i = 1 \]

2. Ordinary Kriging (OK)
   Rataan tidak diketahui tetapi konstan: \[ \sum_{i=1}^n \lambda_i = 1 \]

3. Universal Kriging (UK)
   Rataan berupa fungsi lokasi: \[ \mu(\mathbf{s}) = \sum_{k=0}^p \beta_k f_k(\mathbf{s}) \]

10. Sistem Persamaan Kriging dan Variansi Prediksi

Penduga Kriging: \[ \hat{Z}(\mathbf{s}_0)=\sum_{i=1}^n \lambda_i Z(\mathbf{s}_i) \]

Sistem persamaan Kriging: \[ \begin{bmatrix} \Gamma & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \gamma_0 \\ 1 \end{bmatrix} \]

Variansi Kriging: \[ \sigma_K^2 = \sum_{i=1}^n \lambda_i \gamma(\mathbf{s}_i - \mathbf{s}_0) + \mu \]

11. Model Variogram

Model variogram dipilih dengan mencocokkan variogram teoretis terhadap variogram empiris melalui metode least squares atau maximum likelihood.

12. Varian Kriging

Varian Kriging bergantung pada: - Distribusi lokasi sampel - Model variogram - Jarak titik prediksi ke titik pengamatan

Semakin dekat dan rapat titik pengamatan, semakin kecil varian Kriging.

13. Implementasi Dasar (CRS-Safe)

Implementasi Kriging harus memperhatikan: - Sistem koordinat proyeksi (CRS) - Konsistensi satuan jarak - Transformasi data jika diperlukan

Penggunaan CRS yang tidak sesuai dapat menyebabkan bias jarak dan kesalahan prediksi.

2.12 Tren Surface Analysis

Analisis Permukaan Tren adalah sebuah metode statistik yang digunakan untuk memodelkan variasi spasial (kecenderungan perubahan nilai) dari suatu variabel di seluruh wilayah studi. Metode ini mengasumsikan bahwa variasi nilai data (response variable) dapat dijelaskan sebagai fungsi polinomial dari koordinat geografis (longitudo dan latitudo).

1. Model regresi orde pertama yang menggunakan koordinat \(x\) (longitudo) dan \(y\) (latitudo) sebagai prediktor dapat dinyatakan sebagai berikut:

\[ Z(s_i) = \beta_0 + \beta_1 x_i + \beta_2 y_i + \varepsilon_i \]

• Keterangan:
  • \(Z(s_i)\): Nilai data pada lokasi \(s_i\).
  • \(i\): Indeks observasi, \(i = 1, \dots, n\).
  • \(s_i = (x_i, y_i)\): Koordinat longitudo (\(x_i\)) dan latitudo (\(y_i\)) secara berturut-turut.
  • \(\beta_0\): Intersep atau konstanta model.
  • \(\beta_1\) dan \(\beta_2\): Koefisien estimasi untuk koordinat longitudo dan latitudo.
  • \(\varepsilon_i\): Error residu pada observasi \(i\).
• Interpretasi Koefisien: Tanda koefisien yang diestimasi dari model ini akan memberi tahu kita apakah nilai variabel meningkat atau menurun seiring pergerakan ke Utara (\(y\)) atau ke Timur (\(x\)).

2. Model Regresi Polinomial Kompleks (Orde Lebih Tinggi)

Model yang lebih rumit menggunakan polinomial berorde lebih tinggi pada koordinat juga dapat dipertimbangkan. Tujuannya adalah untuk menyesuaikan permukaan yang lebih kompleks, seperti:

\[ Z(s_i) = \beta_0 + \beta_1 x_i + \beta_2 y_i + \beta_3 x_i^2 + \beta_4 y_i^2 + \beta_5 x_i y_i + \varepsilon_i \]

BAB III (METODOLOGI PENELITIAN)

3.1 Data dan Sumber Data

Penelitian ini menggunakan data sekunder berbentuk data panel spasial yang mencakup 38 dan atau 34 provinsi di Indonesia selama periode tahun 2020 hingga 2025. Data diperoleh dari sumber resmi, yaitu Badan Pusat Statistik (BPS).

Jenis data yang digunakan , yaitu Data Tingkat Pengangguran Terbuka (TPT), Data Lama.Sekolah, dan Data Tingkat Partisipasi Angkatan Kerja (TPAK) . Ketiga variabel tersebut disusun dalam bentuk panel data (cross section dan time series), sehingga memungkinkan analisis hubungan spasial dan temporal antarprovinsi di Indonesia. Data diolah menggunakan software R Studio untuk melakukan estimasi model OLS, SAR, SEM, dan SDM.

3.2 Variabel Penelitian

Penelitian ini melibatkan satu variabel dependen dan dua variabel independen. Definisi operasional variabel dijelaskan sebagai berikut:

Jenis Variabel	Nama Variabel	Simbol	Satuan	Keterangan
Dependen	Tingkat Pengangguran Terbuka	\(Y\)	Persen (%)	Persentase jumlah penganggur terhadap total angkatan kerja di suatu provinsi.
Independen	Tingkat Partisipasi Angkatan Kerja	\(X_1\)	Persen (%)	Persentase TPAK usia kerja (15 tahun ke atas) yang bekerja atau aktif mencari pekerjaan.
Independen	Lama Sekolah	\(X_2\)	Rata-Rata	Rata-Rata Lama sekolah di suatu provinsi.

Fokus utama penelitian ini adalah untuk melihat pola keterkaitan spasial (spatial dependence) dan autokorelasi spasial (spatial autocorrelation) antarprovinsi dalam hal tingkat pengangguran terbuka, baik secara global maupun lokal. Serta model spasial (SAR atau SEM) yang paling tepat menggambarkan hubungan tersebut.

3.3 Metode Analisis Data

Analisis dilakukan menggunakan pendekatan spasial dengan data panel, yang bertujuan untuk mengetahui pengaruh faktor-faktor ekonomi terhadap TPT dengan mempertimbangkan keterkaitan antarwilayah. Analisis dilakukan melalui beberapa tahapan, yaitu analisis deskriptif, pengujian autokorelasi spasial, dan estimasi model spasial (SAR dan SEM).

3.3.1 Eksplorasi Data dan Visualisasi

Menurut Sudjana (2005), analisis deskriptif merupakan metode statistik yang bertujuan untuk memberikan gambaran atau deskripsi mengenai suatu data sebagaimana adanya tanpa bermaksud untuk menarik kesimpulan yang bersifat umum. Analisis ini digunakan untuk menyajikan data dalam bentuk tabel, grafik, maupun ukuran-ukuran statistik seperti rata-rata, median, modus, dan simpangan baku agar informasi yang terkandung dalam data dapat lebih mudah dipahami.

Tahap awal penelitian ini berupa analisis deskriptif yang bertujuan untuk memberikan gambaran umum mengenai perkembangan beberapa variabel di Indonesia dari tahun 2020 hingga 2025. Analisis ini mencakup penyajian nilai rata-rata, nilai maksimum dan minimum untuk variabelTingkat Pengangguran Terbuka (TPT), Partisipasi Angkatan Kerja (Lama.Sekolah), dan TPAK TPAK.

Selain analisis numerik, dilakukan pula visualisasi spasial menggunakan peta tematik (choropleth map) untuk menunjukkan pola sebaran variabel antarprovinsi di Indonesia. Visualisasi ini membantu mengidentifikasi apakah terdapat pola pengelompokan wilayah (spatial clustering), misalnya provinsi dengan TPT tinggi yang saling berdekatan (hotspot) atau wilayah dengan TPT rendah yang membentuk klaster tersendiri (coldspot).

3.3.2 Uji Autokorelasi Spasial

Setelah dilakukan pemetaan untuk melihat pola sebaran TPT antarprovinsi, langkah berikutnya adalah melakukan uji autokorelasi spasial. Tujuan pengujian ini adalah untuk memastikan apakah nilai TPT di suatu provinsi memiliki hubungan atau ketergantungan dengan nilai TPT di provinsi yang berdekatan. Dengan kata lain, apakah tingkat pengangguran terbuka di suatu wilayah dipengaruhi oleh kondisi wilayah sekitarnya. Dalam penelitian ini digunakan tiga jenis uji autokorelasi spasial, yaitu Moran’s I, Geary’s C, dan Getis-Ord G.

Hasil dari ketiga uji ini memberikan dasar untuk menentukan model spasial yang paling sesuai, yaitu apakah menggunakan model Spatial Autoregressive (SAR) atau Spatial Error Model (SEM). Jika pengaruh spasial muncul langsung pada variabel dependen, maka model SAR lebih tepat digunakan. Namun, jika pengaruh spasial terdapat pada komponen galat (error term), maka model SEM menjadi pilihan yang lebih sesuai.

3.3.3 Estimasi Model Spatial

Setelah hasil uji Moran’s I menunjukkan adanya autokorelasi spasial, langkah selanjutnya adalah menentukan bentuk dependensi spasial yang paling sesuai. Dalam analisis spasial, terdapat dua model utama yang umum digunakan, yaitu Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM).

Kedua model tersebut sama-sama mempertimbangkan adanya ketergantungan spasial antarwilayah, namun perbedaannya terletak pada sumber ketergantungan yang terjadi. Model SAR digunakan ketika ketergantungan spasial muncul secara langsung melalui variabel dependen wilayah tetangga. Artinya, nilai Tingkat Pengangguran Terbuka (TPT) suatu provinsi dipengaruhi oleh nilai TPT di provinsi lain yang berdekatan. Sebaliknya, model SEM digunakan apabila ketergantungan spasial tidak muncul secara langsung pada variabel dependen, melainkan melalui komponen error, yaitu adanya faktor-faktor tak teramati yang saling berkorelasi antarwilayah.

Perbandingan kedua model ini bertujuan untuk mengetahui bentuk keterkaitan spasial yang paling sesuai dengan kondisi empiris di Indonesia. Jika hasil estimasi menunjukkan bahwa model SAR lebih unggul (ditunjukkan oleh nilai log-likelihood yang lebih tinggi atau nilai AIC/BIC yang lebih rendah), maka dapat diartikan bahwa dinamika TPT antarprovinsi saling memengaruhi secara langsung. Namun, apabila model SEM lebih sesuai, maka hubungan antarprovinsi lebih disebabkan oleh pengaruh eksternal yang tidak terukur, seperti kebijakan nasional, kondisi makroekonomi, atau shock eksternal yang berdampak serupa di beberapa provinsi sekaligus.

3.3.4 Pemilihan Model Terbaik

Pemilihan model terbaik antara Spatial Autoregressive (SAR) dan Spatial Error Model (SEM) dilakukan untuk menentukan model yang paling sesuai dalam menggambarkan hubungan antarwilayah pada data panel TPT provinsi di Indonesia tahun 2020–2025. Perbandingan dilakukan dengan melihat nilai koefisien determinasi (R²) dan Akaike Information Criterion (AIC). Model yang memiliki nilai R² tertinggi menunjukkan kemampuan yang lebih baik dalam menjelaskan variasi, sedangkan nilai AIC yang lebih rendah menandakan model yang lebih efisien dan sesuai dengan data. Dengan demikian, model terbaik dipilih berdasarkan kombinasi dari kedua kriteria tersebut untuk memastikan bentuk pengaruh spasial yang paling tepat, apakah berada pada variabel dependen (SAR) atau pada komponen galat (SEM).

Anda benar! BAB III (Metodologi) Anda kurang menjelaskan langkah-langkah analisis yang akan dilakukan di BAB IV. Berikut adalah tambahan yang diperlukan untuk BAB III:

3.3.5 Analisis Geographically Weighted Regression (GWR)

Setelah pemilihan model spasial global terbaik, dilakukan analisis heterogenitas spasial menggunakan Geographically Weighted Regression (GWR) untuk mengidentifikasi variasi pengaruh variabel independen terhadap TPT pada setiap lokasi. Tahapan analisis GWR meliputi:

3.3.5.1 Transformasi Sistem Koordinat

Data spasial yang semula dalam sistem koordinat geografis (longitude-latitude, WGS84) ditransformasikan ke sistem koordinat proyeksi Universal Transverse Mercator (UTM) untuk memastikan akurasi perhitungan jarak Euclidean. Provinsi Indonesia yang tersebar di berbagai zona UTM akan disesuaikan dengan zona yang sesuai (misalnya UTM Zona 48S, 49S, 50S, atau 51S).

3.3.5.2 Penentuan Bandwidth Optimal

Bandwidth optimal ditentukan menggunakan dua metode:

1. Cross Validation (CV): Meminimalkan fungsi CV dengan menghilangkan observasi ke-i dan memprediksi nilainya menggunakan data lainnya \[CV = \sum_{i=1}^{n} [y_i - \hat{y}_{\ne i}(h)]^2\]

2. Akaike Information Criterion Corrected (AICc): Memilih bandwidth yang menghasilkan nilai AICc minimum \[\text{AICc} = 2n \ln(\hat{\sigma}) + n \ln(2\pi) + n \left\{ \frac{n+\text{tr}(\mathbf{L})}{n-2-\text{tr}(\mathbf{L})} \right\}\]

Bandwidth dapat bersifat fixed (jarak tetap) atau adaptive (jumlah tetangga tetap). Dalam penelitian ini digunakan bandwidth adaptive untuk mengakomodasi kepadatan observasi yang tidak seragam.

3.3.5.3 Fungsi Kernel

Model GWR menggunakan fungsi kernel untuk pembobotan spasial. Dua fungsi kernel yang umum digunakan:

1. Gaussian Kernel: \[w_i = \exp\left[-\frac{1}{2}\left(\frac{d_i}{b}\right)^2\right]\]

2. Bi-square Kernel: \[w_i = \begin{cases} \left[1-\left(\frac{d_i}{b}\right)^2\right]^2 & \text{jika } d_i < b \\ 0 & \text{jika } d_i \ge b \end{cases}\]

di mana \(w_i\) adalah bobot spasial, \(d_i\) adalah jarak antara titik pengamatan dan titik estimasi, dan \(b\) adalah bandwidth.

3.3.5.4 Estimasi Koefisien Lokal

Untuk setiap lokasi \(i\), koefisien regresi dihitung menggunakan weighted least squares:

\[\hat{\boldsymbol{\beta}}(u_i, v_i) = (\mathbf{X}^T \mathbf{W}(u_i, v_i) \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}(u_i, v_i) \mathbf{y}\]

di mana \(\mathbf{W}(u_i, v_i)\) adalah matriks diagonal berisi bobot spasial untuk lokasi \((u_i, v_i)\).

3.3.5.5 Uji Asumsi Model GWR

Setelah estimasi model GWR, dilakukan pengujian asumsi:

1. Normalitas Residual: Uji Jarque-Bera \[JB = \frac{n}{6}\left(S^2 + \frac{(K-3)^2}{4}\right)\]

2. Homoskedastisitas: Uji Breusch-Pagan \[LM = n \cdot R^2\]

3. Non-Autokorelasi Spasial: Uji Moran’s I pada residual GWR

4. Non-Multikolinearitas: Variance Inflation Factor (VIF) lokal

3.3.5.6 Perbandingan Model Global vs Lokal

Perbandingan dilakukan menggunakan: - AICc: Model dengan AICc lebih rendah lebih baik - R²: Koefisien determinasi yang lebih tinggi menunjukkan kecocokan lebih baik - Visualisasi koefisien lokal: Peta spasial untuk melihat variasi geografis

---

3.3.6 Interpolasi Spasial

Untuk memprediksi nilai TPT pada lokasi yang tidak tersampel dan menghasilkan permukaan kontinu, dilakukan interpolasi spasial menggunakan berbagai metode deterministik dan geostatistik.

3.3.6.1 Metode Deterministik

a. Nearest Neighbor (NN)

Metode paling sederhana yang memberikan nilai titik terdekat: \[\hat{Z}(s_0) = Z(s_{(1)})\]

di mana \(s_{(1)}\) adalah lokasi observasi terdekat dari \(s_0\). Visualisasi dilakukan menggunakan diagram Voronoi.

b. Natural Neighbor/Triangulated Irregular Network (TIN)

Interpolasi linear pada segitiga Delaunay: - Membentuk jaringan segitiga dari titik-titik observasi - Nilai di dalam segitiga dihitung menggunakan interpolasi barycentric dari tiga vertex - Menghasilkan permukaan piecewise linear

c. Inverse Distance Weighted (IDW)

Pembobotan berbasis jarak dengan parameter power \(p\): \[\hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i(s_0) Z(s_i)}{\sum_{i=1}^{n} w_i(s_0)}\]

di mana: \[w_i(s_0) = d(s_0, s_i)^{-p}\]

Penelitian ini menggunakan \(p = 2\) (inverse distance squared) sebagai nilai default.

d. Trend Surface Analysis/Thin Plate Splines (TPS)

Interpolasi menggunakan spline dengan meminimalkan fungsi energi bending: - Menghasilkan permukaan yang sangat halus - Mampu menangkap tren global dan variasi lokal - Parameter smoothing (\(\lambda\)) ditentukan otomatis via Generalized Cross Validation (GCV)

3.3.6.2 Metode Geostatistik (Kriging)

Metode Kriging didasarkan pada teori variabel teregionalisasi dan menghasilkan Best Linear Unbiased Predictor (BLUP).

a. Tahapan Analisis Kriging

1. Analisis Variogram Empirik

   Menghitung semivariance untuk berbagai lag distance: \[\hat{\gamma}(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} [Z(s_i) - Z(s_i + h)]^2\]

2. Fitting Model Variogram Teoretis

   Model yang diuji:

   • Spherical: \[\gamma(h) = \begin{cases} C_0 + C\left[\frac{3h}{2a} - \frac{h^3}{2a^3}\right] & h \le a \\ C_0 + C & h > a \end{cases}\]

   • Exponential: \[\gamma(h) = C_0 + C\left(1 - e^{-h/a}\right)\]

   • Gaussian: \[\gamma(h) = C_0 + C\left(1 - e^{-(h^2/a^2)}\right)\]

   di mana:

   • \(C_0\) = nugget effect
   • \(C\) = partial sill
   • \(a\) = range
   • \(C_0 + C\) = total sill

3. Uji Stasioneritas dan Isotropi

   • Stasioneritas: Variogram tidak bergantung pada lokasi absolut
   • Isotropi: Variogram hanya bergantung pada jarak, bukan arah
   • Jika anisotropi terdeteksi, dilakukan transformasi koordinat

b. Ordinary Kriging (OK)

Prediksi linier tak bias dengan varians minimum: \[\hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i)\]

dengan constraint: \[\sum_{i=1}^{n} \lambda_i = 1\]

Sistem persamaan Kriging: \[\begin{bmatrix} \Gamma & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \gamma_0 \\ 1 \end{bmatrix}\]

Varians Kriging: \[\sigma_K^2 = \sum_{i=1}^{n} \lambda_i \gamma(s_i - s_0) + \mu\]

c. Universal Kriging (UK)

Mengakomodasi tren deterministik (drift): \[Z(s) = m(s) + \varepsilon(s)\]

di mana \(m(s) = \sum_{k=0}^{p} \beta_k f_k(s)\) adalah fungsi tren (misalnya: \(f_1(s) = x\), \(f_2(s) = y\)).

d. Co-Kriging

Interpolasi multivariat yang memanfaatkan korelasi silang antar variabel: - Menggunakan variabel sekunder (misalnya: lama sekolah) untuk meningkatkan prediksi variabel primer (TPT) - Model koregionalisasi: Intrinsic Correlation Model (ICM) atau Linear Model of Coregionalization (LMC) - Cocok ketika variabel sekunder memiliki sampel lebih padat

3.3.6.3 Validasi dan Evaluasi Model Interpolasi

Semua metode interpolasi dievaluasi menggunakan k-fold cross-validation (k=10) dengan metrik:

1. Root Mean Square Error (RMSE): \[\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n} (\hat{Z}(s_i) - Z(s_i))^2}\]

2. Mean Absolute Error (MAE): \[\text{MAE} = \frac{1}{n}\sum_{i=1}^{n} |\hat{Z}(s_i) - Z(s_i)|\]

3. Mean Error (ME) - mengukur bias: \[\text{ME} = \frac{1}{n}\sum_{i=1}^{n} (\hat{Z}(s_i) - Z(s_i))\]

Model terbaik dipilih berdasarkan: - RMSE dan MAE terkecil (akurasi tertinggi) - ME mendekati 0 (bias minimal) - Varians Kriging terendah (untuk metode Kriging) - Visualisasi peta residual (tidak ada pola spasial sistematis)

3.3.6.4 Implementasi CRS-Safe

Untuk memastikan konsistensi perhitungan: 1. Semua data ditransformasikan ke sistem koordinat proyeksi yang sesuai (UTM) 2. Satuan jarak distandarisasi (meter) 3. Validasi CRS dilakukan sebelum setiap operasi spasial 4. Handling khusus untuk data yang melintasi zona UTM berbeda

---

3.3.7 Visualisasi dan Pemetaan Hasil

Hasil analisis divisualisasikan menggunakan:

1. Peta Tematik (Choropleth Maps)
   • Distribusi TPT, TPAK, dan lama sekolah per provinsi
   • Klasifikasi menggunakan metode natural breaks (Jenks)
2. Peta LISA (Local Indicators of Spatial Association)
   • Identifikasi hotspot (High-High) dan coldspot (Low-Low)
   • Deteksi spatial outliers (High-Low, Low-High)
3. Peta Koefisien Lokal GWR
   • Visualisasi variasi spasial koefisien regresi
   • Surface plot menggunakan interpolasi IDW
4. Peta Prediksi dan Varians Kriging
   • Permukaan kontinu hasil interpolasi
   • Peta ketidakpastian (prediction variance)
5. Scatter Plot dan Diagnostic Plot
   • Observed vs Predicted
   • Residual plot untuk validasi asumsi

Alur Kerja Penelitian

Berikut adalah alur kerja penelitian yang dilakukan :

┌─────────────────────────────┐
│  Pengumpulan Data           │
│  (BPS: TPT, TPAK, Lama      │
│   Sekolah, Shapefile)       │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Eksplorasi & Visualisasi   │
│  - Statistik Deskriptif     │
│  - Peta Tematik             │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Uji Autokorelasi Spasial   │
│  - Moran's I                │
│  - Geary's C                │
│  - Getis-Ord G              │
│  - LISA                     │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Model Regresi Spasial      │
│  - OLS (baseline)           │
│  - SAR                      │
│  - SEM                      │
│  - SDM                      │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Pemilihan Model Terbaik    │
│  - Perbandingan AIC         │
│  - Evaluasi R²              │
│  - Uji diagnostik           │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Geographically Weighted    │
│  Regression (GWR)           │
│  - Bandwidth optimal        │
│  - Estimasi koefisien lokal │
│  - Uji asumsi               │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Interpolasi Spasial        │
│  Deterministik:             │
│  - NN, TIN, IDW, TPS        │
│  Geostatistik:              │
│  - OK, UK, Co-Kriging       │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Validasi & Evaluasi        │
│  - Cross-validation         │
│  - RMSE, MAE, ME            │
│  - Perbandingan metode      │
└──────────┬──────────────────┘
           │
           ▼
┌─────────────────────────────┐
│  Interpretasi & Kesimpulan  │
│  - Implikasi kebijakan      │
│  - Rekomendasi              │
└─────────────────────────────┘

BAB IV (HASIL DAN PEMBAHASAN)

4.1 Eksplorasi Data dan Visualisasi

# Gabungkan shapefile dengan data CSV
data_sf <- left_join(indo_shp, data, by = c("NAME_1" = "Provinsi"))

data_sf <- data_sf %>%
  mutate(
    TPT = as.numeric(gsub(",", ".", TPT)),
    Lama.Sekolah = as.numeric(gsub(",", ".", Lama.Sekolah)),
    TPAK = as.numeric(gsub(",", ".", TPAK))
  )

data_sf <- data_sf %>%
  filter(!is.na(TPT) & !is.na(Lama.Sekolah) & !is.na(TPAK))

summary(select(data_sf, TPT, Lama.Sekolah, TPAK))

##       TPT         Lama.Sekolah         TPAK                geometry  
##  Min.   :1.250   Min.   : 6.960   Min.   :61.85   MULTIPOLYGON :162  
##  1st Qu.:3.663   1st Qu.: 8.682   1st Qu.:66.43   POLYGON      : 12  
##  Median :4.360   Median : 9.230   Median :69.53   epsg:NA      :  0  
##  Mean   :4.711   Mean   : 9.201   Mean   :69.27   +proj=long...:  0  
##  3rd Qu.:5.750   3rd Qu.: 9.740   3rd Qu.:71.43                      
##  Max.   :9.010   Max.   :10.700   Max.   :80.23

Interpretasi:

Hasil statistik deskriptif menunjukkan bahwa rata-rata Tingkat Pengangguran Terbuka (TPT) di Indonesia sebesar 4,71 persen, dengan nilai minimum 1,25 persen dan maksimum 9,01 persen. Rentang nilai yang cukup lebar ini mengindikasikan adanya kesenjangan tingkat pengangguran antarprovinsi, di mana beberapa wilayah menghadapi permasalahan pengangguran yang relatif lebih tinggi dibandingkan wilayah lainnya.

Variabel rata-rata lama sekolah memiliki nilai rata-rata sebesar 9,20 tahun, dengan nilai minimum 6,96 tahun dan maksimum 10,70 tahun. Hal ini menunjukkan bahwa secara umum penduduk Indonesia telah menempuh pendidikan hingga jenjang sekolah menengah pertama, namun masih terdapat disparitas kualitas pendidikan antarwilayah. Provinsi dengan rata-rata lama sekolah yang lebih rendah mencerminkan keterbatasan akses atau keberlanjutan pendidikan, yang berpotensi memengaruhi daya saing tenaga kerja di pasar kerja.

Sementara itu, Tingkat Partisipasi Angkatan Kerja (TPAK) memiliki nilai rata-rata sebesar 69,27 persen, dengan rentang antara 61,85 persen hingga 80,23 persen. Nilai ini menunjukkan bahwa sebagian besar penduduk usia kerja di Indonesia aktif secara ekonomi, baik bekerja maupun sedang mencari pekerjaan. Variasi TPAK yang cukup besar antarprovinsi mencerminkan adanya perbedaan struktur ekonomi, kesempatan kerja, serta karakteristik sosial ekonomi di masing-masing wilayah.

Secara keseluruhan, statistik deskriptif ini mengindikasikan bahwa permasalahan pengangguran di Indonesia tidak hanya berkaitan dengan tingkat partisipasi angkatan kerja, tetapi juga dengan kualitas sumber daya manusia yang tercermin dari rata-rata lama sekolah, serta adanya ketimpangan spasial antarprovinsi. Kondisi ini memperkuat pentingnya penggunaan pendekatan analisis spasial dalam mengkaji faktor-faktor yang memengaruhi Tingkat Pengangguran Terbuka di Indonesia.

tm_shape(data_sf) +
  tm_polygons("TPT", palette = "YlOrRd", title = "Tingkat Pengangguran Terbuka") +
  tm_layout(frame = FALSE)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"
## Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.

Interpretasi:

Peta sebaran Tingkat Pengangguran Terbuka (TPT) di atas memperlihatkan adanya variasi spasial yang cukup jelas antarprovinsi di Indonesia. Wilayah dengan tingkat pengangguran tinggi (warna merah tua) tampak terkonsentrasi di bagian timur, terutama Papua, sementara daerah dengan tingkat pengangguran rendah (warna kuning muda hingga oranye) tersebar di sebagian besar wilayah tengah dan barat Indonesia seperti Kalimantan, Sulawesi, dan sebagian Sumatra. Pola ini menunjukkan adanya ketimpangan spasial dalam tingkat pengangguran, yang dapat dipengaruhi oleh perbedaan struktur ekonomi, ketersediaan lapangan kerja, dan tingkat pembangunan antarwilayah.

4.2 Uji Autokorelasi Spasial

4.2.1 Uji Moran’s I

library(spdep)

var_list <- c("TPT", "Lama.Sekolah", "TPAK")

# Jalankan Moran's I pada data yang sudah merge
moran_results <- lapply(var_list, function(v) {
  moran.test(data_sf_simple[[v]], lw, zero.policy = TRUE)
})

names(moran_results) <- var_list

# Tabel hasil Moran’s I
moran_df <- data.frame(
  Variabel = var_list,
  Moran_I = sapply(moran_results, function(x) x$estimate["Moran I statistic"]),
  p_value = sapply(moran_results, function(x) x$p.value)
)

moran_df

##                                    Variabel   Moran_I      p_value
## TPT.Moran I statistic                   TPT 0.4857777 3.431744e-62
## Lama.Sekolah.Moran I statistic Lama.Sekolah 0.3341424 9.035991e-31
## TPAK.Moran I statistic                 TPAK 0.5398854 4.737919e-76

Interpretasi:

Hasil uji Moran’s I menunjukkan bahwa ketiga variabel memiliki autokorelasi spasial positif, yang berarti nilai suatu provinsi cenderung mirip dengan provinsi di sekitarnya. Nilai Moran’s I tertinggi terdapat pada variabel TPAK (0,539), mengindikasikan adanya pengelompokan wilayah dengan tingkat ke TPAKan yang serupa (tinggi atau rendah) secara spasial. Sementara itu, TPT (0,485) dan Lama.Sekolah (0,33) juga menunjukkan pola spasial positif namun dengan kekuatan yang lebih lemah, menandakan bahwa pengangguran dan Lama sekolah antarprovinsi masih memiliki keterkaitan spasial, meski tidak sekuat TPAK.

4.2.2 Uji Geary’s C

library(spdep)

var_list <- c("TPT", "Lama.Sekolah", "TPAK")

# Jalankan Geary's C pada data yang sudah merge
geary_results <- lapply(var_list, function(v) {
  geary.test(data_sf_simple[[v]], lw, zero.policy = TRUE)
})

names(geary_results) <- var_list

# Tabel hasil Geary’s C
geary_df <- data.frame(
  Variabel = var_list,
  Geary_C = sapply(geary_results, function(x) x$estimate["Geary C statistic"]),
  p_value = sapply(geary_results, function(x) x$p.value)
)

geary_df

##                                    Variabel   Geary_C      p_value
## TPT.Geary C statistic                   TPT 0.5133921 3.543808e-55
## Lama.Sekolah.Geary C statistic Lama.Sekolah 0.6517548 1.666902e-29
## TPAK.Geary C statistic                 TPAK 0.4417580 1.496024e-72

Interpretasi:

Hasil uji Geary’s C menunjukkan bahwa seluruh variabel yang dianalisis, yaitu Tingkat Pengangguran Terbuka (TPT), rata-rata lama sekolah, dan Tingkat Partisipasi Angkatan Kerja (TPAK), memiliki p-value < 0,05. Hal ini mengindikasikan adanya autokorelasi spasial yang signifikan pada ketiga variabel tersebut, sehingga nilai suatu provinsi dipengaruhi oleh nilai provinsi-provinsi di sekitarnya.

4.2.3 Getis Ord

library(spdep)

var_list <- c("TPT", "Lama.Sekolah", "TPAK")

# Hitung Local Gi* (Getis-Ord Gi*)
getis_results <- lapply(var_list, function(v) {
  localG(data_sf_simple[[v]], lw, zero.policy = TRUE)
})

names(getis_results) <- var_list

# Gabungkan hasil Gi* ke data spasial
for(v in var_list) {
  data_sf_simple[[paste0("Gi_", v)]] <- as.numeric(getis_results[[v]])
}

# Tabel ringkasan rata-rata dan sebaran Gi*
getis_df <- data.frame(
  Variabel = var_list,
  Mean_Gi = sapply(var_list, function(v) mean(data_sf_simple[[paste0("Gi_", v)]], na.rm = TRUE)),
  SD_Gi = sapply(var_list, function(v) sd(data_sf_simple[[paste0("Gi_", v)]], na.rm = TRUE))
)

getis_df

##                  Variabel     Mean_Gi    SD_Gi
## TPT                   TPT -0.01664867 2.573057
## Lama.Sekolah Lama.Sekolah  0.05954730 2.016202
## TPAK                 TPAK -0.09151381 2.115393

Interpretasi:

Hasil analisis statistik Getis–Ord Gi* menunjukkan bahwa nilai rata-rata (Mean Gi) untuk seluruh variabel, yaitu Tingkat Pengangguran Terbuka (TPT), rata-rata lama sekolah, dan Tingkat Partisipasi Angkatan Kerja (TPAK), berada di sekitar nol. Nilai Mean Gi untuk TPT sebesar −0,017, untuk rata-rata lama sekolah sebesar 0,060, dan untuk TPAK sebesar −0,092. Kondisi ini mengindikasikan bahwa secara umum tidak terdapat kecenderungan dominasi global hotspot atau coldspot pada skala nasional, melainkan pola klaster yang bersifat lokal.

Sedangkan lama sekolah Nilai Mean Gi* yang sangat dekat dengan nol menunjukkan bahwa secara global tidak terdapat kecenderungan dominasi hotspot atau coldspot pendidikan di Indonesia. Artinya, pola pendidikan tidak terkonsentrasi secara nasional pada satu sisi tertentu, melainkan tersebar.

4.2.4 Local Moran’s I (LISA)

library(spdep)
library(dplyr)
library(ggplot2)

# Variabel yang dianalisis
x <- scale(data_sf_simple$TPT)[, 1]   # standardisasi nilai

# Lag spasial
lagx <- spdep::lag.listw(lw, x, zero.policy = TRUE)

# Hitung Local Moran's I
lisa <- spdep::localmoran(x, lw, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii", "Ei", "Vi", "Zi", "Pi.two.sided")

# Tentukan signifikansi
alpha <- 0.05

# Gabungkan hasil dengan data spasial
Indo_LISA <- bind_cols(data_sf_simple, lisa_df) %>%
  mutate(
    x = as.numeric(x),
    lagx = as.numeric(lagx),
    quad = case_when(
      x >= 0 & lagx >= 0 & Pi.two.sided <= alpha ~ "High-High",
      x < 0 & lagx < 0 & Pi.two.sided <= alpha ~ "Low-Low",
      x >= 0 & lagx < 0 & Pi.two.sided <= alpha ~ "High-Low (Outlier)",
      x < 0 & lagx >= 0 & Pi.two.sided <= alpha ~ "Low-High (Outlier)",
      TRUE ~ "Not significant"
    )
  )

# Plot LISA
ggplot(Indo_LISA) +
  geom_sf(aes(fill = quad), color = "white", size = 0.2) +
  scale_fill_manual(
    values = c(
      "High-High" = "#d73027",
      "Low-Low" = "#4575b4",
      "High-Low (Outlier)" = "#fdae61",
      "Low-High (Outlier)" = "#74add1",
      "Not significant" = "grey85"
    )
  ) +
  labs(
    title = "Local Moran's I (LISA) - TPT Indonesia",
    fill = "Kategori LISA"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

Interpretasi:

Peta Local Moran’s I (LISA) untuk TPT Indonesia menunjukkan bahwa terlihat adanya klaster High–High (merah), yaitu wilayah dengan TPT tinggi yang dikelilingi oleh wilayah dengan TPT tinggi pula, yang menandakan konsentrasi pengangguran tinggi secara regional. Sebaliknya, klaster Low–Low (biru tua) mengindikasikan wilayah dengan TPT rendah yang berdekatan dengan wilayah lain yang juga memiliki TPT rendah, mencerminkan kondisi ketenagakerjaan yang relatif lebih baik dan saling berkaitan secara geografis. Selain itu, terdapat wilayah High–Low dan Low–High yang berperan sebagai outlier spasial, menunjukkan provinsi dengan TPT tinggi (atau rendah) yang dikelilingi oleh provinsi dengan kondisi sebaliknya

4.3 Estimasi Model Spasial

4.3.1 Ordinary Least Squares (OLS)

library(spdep)
library(lmtest)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sandwich)

## Warning: package 'sandwich' was built under R version 4.3.3

# 1. OLS Model
ols_model <- lm(TPT ~ Lama.Sekolah + TPAK, data = data_sf)
summary(ols_model)

## 
## Call:
## lm(formula = TPT ~ Lama.Sekolah + TPAK, data = data_sf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8885 -0.7589 -0.1164  0.6278  3.3571 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15.06920    2.56497   5.875 2.16e-08 ***
## Lama.Sekolah  0.44771    0.12921   3.465  0.00067 ***
## TPAK         -0.20900    0.02674  -7.817 5.28e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.128 on 171 degrees of freedom
## Multiple R-squared:  0.3948, Adjusted R-squared:  0.3877 
## F-statistic: 55.78 on 2 and 171 DF,  p-value: < 2.2e-16

# 2. Extract residuals
resid_ols <- residuals(ols_model)

# 4. Uji autokorelasi residuals (Moran's I)
moran_resid <- moran.test(resid_ols, lw, zero.policy = TRUE)
moran_resid

## 
##  Moran I test under randomisation
## 
## data:  resid_ols  
## weights: lw    
## 
## Moran I statistic standard deviate = 16.414, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.4802454136     -0.0057803468      0.0008767624

# 5. Uji heteroskedastisitas (Breusch-Pagan)
bptest(ols_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 2.9348, df = 2, p-value = 0.2305

Interpretasi:

Hasil estimasi regresi linier menunjukkan bahwa rata-rata lama sekolah dan Tingkat Partisipasi Angkatan Kerja (TPAK) berpengaruh signifikan terhadap Tingkat Pengangguran Terbuka (TPT). Secara parsial, lama sekolah berpengaruh positif dan signifikan, yang mengindikasikan bahwa peningkatan rata-rata lama sekolah cenderung diikuti oleh kenaikan TPT, sedangkan TPAK berpengaruh negatif dan signifikan, yang menunjukkan bahwa semakin tinggi partisipasi angkatan kerja maka tingkat pengangguran terbuka cenderung menurun. Secara simultan, model signifikan dengan nilai F-statistic yang tinggi, dan mampu menjelaskan sekitar 39% variasi TPT

uji Moran’s I pada residual OLS menghasilkan nilai yang positif dan sangat signifikan, yang menandakan adanya autokorelasi spasial sehingga asumsi independensi residual tidak terpenuhi. Kondisi ini menunjukkan bahwa variasi TPT di suatu provinsi tidak hanya dipengaruhi oleh karakteristik internal wilayah tersebut, tetapi juga oleh kondisi wilayah di sekitarnya. Oleh karena itu, penggunaan regresi linier biasa belum memadai untuk menangkap keterkaitan antarwilayah. Untuk mengakomodasi ketergantungan spasial tersebut dan memperoleh estimasi parameter yang tidak bias serta lebih efisien, diperlukan penerapan model ekonometrika spasial, seperti Spatial Autoregressive Model (SAR) dan Spatial Error Model (SEM), yang mampu menangkap pengaruh interaksi antarwilayah maupun ketergantungan spasial pada komponen galat dalam analisis Tingkat Pengangguran Terbuka di Indonesia.

4.3.2 Spatial Autoregressive Model (SAR)

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.3.3

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 
## 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

sar_model <- lagsarlm(TPT ~ Lama.Sekolah + TPAK,
                      data = data_sf,
                      listw = lw,
                      zero.policy = TRUE)
summary(sar_model)

## 
## Call:lagsarlm(formula = TPT ~ Lama.Sekolah + TPAK, data = data_sf, 
##     listw = lw, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.74166 -0.64398 -0.16242  0.61335  2.45925 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)   7.185333   2.281939  3.1488 0.0016395
## Lama.Sekolah  0.390959   0.106081  3.6855 0.0002283
## TPAK         -0.129728   0.024752 -5.2411 1.596e-07
## 
## Rho: 0.61912, LR test value: 62.369, p-value: 2.8866e-15
## Asymptotic standard error: 0.062222
##     z-value: 9.9502, p-value: < 2.22e-16
## Wald statistic: 99.006, p-value: < 2.22e-16
## 
## Log likelihood: -235.1804 for lag model
## ML residual variance (sigma squared): 0.83071, (sigma: 0.91143)
## Number of observations: 174 
## Number of parameters estimated: 5 
## AIC: 480.36, (AIC for lm: 540.73)
## LM test for residual autocorrelation
## test value: 17.042, p-value: 3.656e-05

resid_sar <- residuals(sar_model)
# Uji Moran's I pada residual SEM
library(spdep)

moran_sar <- moran.test(resid_sar, lw, randomisation = TRUE)
moran_sar

## 
##  Moran I test under randomisation
## 
## data:  resid_sar  
## weights: lw    
## 
## Moran I statistic standard deviate = 3.0113, p-value = 0.001301
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0835313280     -0.0057803468      0.0008796329

4.3.3 Spatial Error Model (SEM)

sem_model <- errorsarlm(TPT ~ Lama.Sekolah + TPAK,
                        data = data_sf,
                        listw = lw,
                        zero.policy = TRUE)
summary(sem_model)

## 
## Call:errorsarlm(formula = TPT ~ Lama.Sekolah + TPAK, data = data_sf, 
##     listw = lw, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.921453 -0.544665 -0.051599  0.485705  2.195569 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  13.582367   2.407506  5.6417 1.684e-08
## Lama.Sekolah  0.522490   0.110959  4.7089 2.491e-06
## TPAK         -0.197686   0.027499 -7.1888 6.535e-13
## 
## Lambda: 0.72098, LR test value: 80.987, p-value: < 2.22e-16
## Asymptotic standard error: 0.060582
##     z-value: 11.901, p-value: < 2.22e-16
## Wald statistic: 141.63, p-value: < 2.22e-16
## 
## Log likelihood: -225.8713 for error model
## ML residual variance (sigma squared): 0.72522, (sigma: 0.8516)
## Number of observations: 174 
## Number of parameters estimated: 5 
## AIC: 461.74, (AIC for lm: 540.73)

# Ambil residual SEM
resid_sem <- residuals(sem_model)

# Uji Moran's I pada residual SEM
library(spdep)

moran_sem <- moran.test(resid_sem, lw, randomisation = TRUE)
moran_sem

## 
##  Moran I test under randomisation
## 
## data:  resid_sem  
## weights: lw    
## 
## Moran I statistic standard deviate = -0.25368, p-value = 0.6001
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0132993228     -0.0057803468      0.0008785122

4.3.4 Spatial Durbin Model (SDM)

sdm_model <- lagsarlm(TPT ~ Lama.Sekolah + TPAK,
                      data = data_sf,
                      listw = lw,
                      type = "mixed",
                      zero.policy = TRUE)
summary(sdm_model)

## 
## Call:lagsarlm(formula = TPT ~ Lama.Sekolah + TPAK, data = data_sf, 
##     listw = lw, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.805286 -0.569256 -0.063874  0.478182  2.167605 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value  Pr(>|z|)
## (Intercept)       5.500298   3.264591  1.6848  0.092020
## Lama.Sekolah      0.524778   0.111274  4.7161 2.404e-06
## TPAK             -0.196039   0.028745 -6.8200 9.102e-12
## lag.Lama.Sekolah -0.459726   0.190413 -2.4144  0.015763
## lag.TPAK          0.127120   0.042796  2.9704  0.002974
## 
## Rho: 0.71834, LR test value: 79.805, p-value: < 2.22e-16
## Asymptotic standard error: 0.061101
##     z-value: 11.756, p-value: < 2.22e-16
## Wald statistic: 138.21, p-value: < 2.22e-16
## 
## Log likelihood: -225.7118 for mixed model
## ML residual variance (sigma squared): 0.72456, (sigma: 0.85121)
## Number of observations: 174 
## Number of parameters estimated: 7 
## AIC: 465.42, (AIC for lm: 543.23)
## LM test for residual autocorrelation
## test value: 2.9435, p-value: 0.086222

4.4 Perbandingan Model

AIC(ols_model, sar_model, sem_model, sdm_model)

##           df      AIC
## ols_model  4 540.7293
## sar_model  5 480.3607
## sem_model  5 461.7427
## sdm_model  7 465.4236

Interpretasi:

Perbandingan nilai Akaike Information Criterion (AIC) menunjukkan bahwa model regresi spasial memberikan kinerja yang lebih baik dibandingkan model OLS. Model OLS memiliki nilai AIC sebesar 540,73, yang merupakan nilai tertinggi, sehingga kurang efisien dalam menjelaskan variasi Tingkat Pengangguran Terbuka (TPT). Di antara model spasial, Spatial Error Model (SEM) menghasilkan nilai AIC terendah sebesar 461,74, diikuti oleh Spatial Durbin Model (SDM) sebesar 465,42, dan Spatial Autoregressive Model (SAR) sebesar 480,36. Nilai AIC yang lebih rendah menunjukkan keseimbangan terbaik antara tingkat kecocokan model dan kompleksitasnya. Dengan demikian, SEM merupakan model terbaik dalam menganalisis TPT di Indonesia, yang mengindikasikan bahwa ketergantungan spasial lebih dominan berasal dari komponen galat dibandingkan dari pengaruh langsung variabel dependen antarwilayah.

4.5 Model Terbaik

Didapatkan model terbaik, yaitu Spatial Error Model (SEM). Dengan model sebagai berikut:

\[ \text{TPT}_i = 13{,}582 + 0{,}522 \cdot \text{Lama.Sekolah}_i - 0{,}198 \cdot \text{TPAK}_i + u_i \]

dengan komponen galat spasial:

\[ u_i = \lambda \sum_{j=1}^{n} w_{ij} u_j + \varepsilon_i \]

---

Keterangan:

• \(\text{TPT}_i\) : Tingkat Pengangguran Terbuka di provinsi ke-\(i\)
  
• \(\text{Lama.Sekolah}_i\) : Rata-rata lama sekolah di provinsi ke-\(i\)
  
• \(\text{TPAK}_i\) : Tingkat Partisipasi Angkatan Kerja di provinsi ke-\(i\)
  
• \(u_i\) : komponen galat yang mengandung dependensi spasial
  
• \(w_{ij}\) : elemen matriks pembobot spasial
  
• \(\lambda\) : parameter autokorelasi spasial galat
  \[ \lambda = 0{,}721 \]
• \(\varepsilon_i\) : galat acak yang diasumsikan identik dan independen

Koefisien rata-rata lama sekolah sebesar 0,522 dan signifikan menunjukkan bahwa setiap kenaikan 1 tahun rata-rata lama sekolah di suatu provinsi, dengan variabel lain konstan, akan meningkatkan Tingkat Pengangguran Terbuka (TPT) sekitar 0,52 persen. Sementara itu, koefisien TPAK sebesar −0,198 dan signifikan mengindikasikan bahwa setiap kenaikan 1 persen TPAK akan menurunkan TPT sekitar 0,20 persen.

4.6 Derajat (Latitude longitude)

library(sf)
library(sp)

## Warning: package 'sp' was built under R version 4.3.3

# Pastikan data_sf adalah sf object
data_sf <- st_as_sf(data_sf)

# Ambil centroid tiap provinsi
centroid_sf <- st_centroid(data_sf)

## Warning: st_centroid assumes attributes are constant over geometries

# Ambil koordinat longitude & latitude
coords <- st_coordinates(centroid_sf)

# Gabungkan koordinat dengan variabel penelitian
data_spasial <- data.frame(
  TPT = centroid_sf$TPT,
  Lama.Sekolah = centroid_sf$Lama.Sekolah,
  TPAK = centroid_sf$TPAK,
  longitude = coords[,1],
  latitude  = coords[,2]
)

# Ubah menjadi SpatialPointsDataFrame (seperti contoh)
coordinates(data_spasial) <- ~longitude+latitude

# (Opsional) Set CRS WGS84
proj4string(data_spasial) <- CRS("+proj=longlat +datum=WGS84")

# Cek hasil
summary(data_spasial)

## Object of class SpatialPointsDataFrame
## Coordinates:
##                 min        max
## longitude 96.912074 138.687642
## latitude  -9.264059   4.224396
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 174
## Data attributes:
##       TPT         Lama.Sekolah         TPAK      
##  Min.   :1.250   Min.   : 6.960   Min.   :61.85  
##  1st Qu.:3.663   1st Qu.: 8.682   1st Qu.:66.43  
##  Median :4.360   Median : 9.230   Median :69.53  
##  Mean   :4.711   Mean   : 9.201   Mean   :69.27  
##  3rd Qu.:5.750   3rd Qu.: 9.740   3rd Qu.:71.43  
##  Max.   :9.010   Max.   :10.700   Max.   :80.23

4.7 Bandwith Optimal

#install.packages("spgwr")
library(spgwr)

## Warning: package 'spgwr' was built under R version 4.3.3

## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)

# 1. Tentukan bandwidth optimal (adaptive)
bandwidth_optimal <- gwr.sel(
  TPT ~ Lama.Sekolah + TPAK,
  data   = data_spasial,
  coords = cbind(
    coordinates(data_spasial)[,1],
    coordinates(data_spasial)[,2]
  ),
  adapt = TRUE
)

## Warning in gwr.sel(TPT ~ Lama.Sekolah + TPAK, data = data_spasial, coords =
## cbind(coordinates(data_spasial)[, : data is Spatial* object, ignoring coords
## argument

## Adaptive q: 0.381966 CV score: 186.1316 
## Adaptive q: 0.618034 CV score: 205.1483 
## Adaptive q: 0.236068 CV score: 155.7979 
## Adaptive q: 0.145898 CV score: 127.6847 
## Adaptive q: 0.09016994 CV score: 97.90985 
## Adaptive q: 0.05572809 CV score: 88.94785 
## Adaptive q: 0.03014629 CV score: 350.4535 
## Adaptive q: 0.06888371 CV score: 97.78015 
## Adaptive q: 0.04595671 CV score: 88.94785 
## Adaptive q: 0.0508424 CV score: 88.94785 
## Adaptive q: 0.04839956 CV score: 88.94785 
## Adaptive q: 0.04746647 CV score: 88.94785 
## Adaptive q: 0.04688979 CV score: 88.94785 
## Adaptive q: 0.04764467 CV score: 88.94785 
## Adaptive q: 0.04726723 CV score: 88.94785 
## Adaptive q: 0.04742578 CV score: 88.94785 
## Adaptive q: 0.04738509 CV score: 88.94785 
## Adaptive q: 0.04742578 CV score: 88.94785

bandwidth_optimal

## [1] 0.04742578

# 2. Estimasi model GWR
model_gwr <- gwr(
  TPT ~ Lama.Sekolah + TPAK,
  data   = data_spasial,
  coords = cbind(
    coordinates(data_spasial)[,1],
    coordinates(data_spasial)[,2]
  ),
  adapt    = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit    = TRUE
)

## Warning in gwr(TPT ~ Lama.Sekolah + TPAK, data = data_spasial, coords =
## cbind(coordinates(data_spasial)[, : data is Spatial* object, ignoring coords
## argument

# 3. Ambil koefisien lokal
coef_gwr <- model_gwr$SDF@data[, c("Lama.Sekolah", "TPAK")]
# 4. Ringkasan koefisien lokal
summary(coef_gwr)

##   Lama.Sekolah            TPAK        
##  Min.   :-0.750911   Min.   :-0.6153  
##  1st Qu.: 0.008672   1st Qu.:-0.3172  
##  Median : 0.507783   Median :-0.2075  
##  Mean   : 0.362239   Mean   :-0.2421  
##  3rd Qu.: 0.674507   3rd Qu.:-0.1260  
##  Max.   : 1.347115   Max.   : 0.1144

Hasil pemilihan bandwidth adaptif menunjukkan bahwa nilai q optimal sebesar 0,0474, yaitu sekitar 4,7% tetangga terdekat, dengan nilai Cross Validation (CV) minimum 88,95, yang menandakan model GWR telah mencapai keseimbangan terbaik antara bias dan varians. Bandwidth yang relatif kecil ini mengindikasikan adanya heterogenitas spasial yang kuat, sehingga hubungan antara variabel penjelas dan TPT bersifat sangat lokal.

4.8 Ringkasan Model GWR

# Ringkasan model GWR
summary(model_gwr)

##           Length Class                  Mode     
## SDF         174  SpatialPointsDataFrame S4       
## lhat      30276  -none-                 numeric  
## lm           11  -none-                 list     
## results      14  -none-                 list     
## bandwidth   174  -none-                 numeric  
## adapt         1  -none-                 numeric  
## hatmatrix     1  -none-                 logical  
## gweight       1  -none-                 character
## gTSS          1  -none-                 numeric  
## this.call     7  -none-                 call     
## fp.given      1  -none-                 logical  
## timings      12  -none-                 numeric

Output menunjukkan bahwa model GWR dibangun menggunakan 174 titik pengamatan yang tersimpan dalam objek SpatialPointsDataFrame (SDF), sehingga setiap lokasi memiliki koordinat dan nilai variabel terkait. Komponen bandwidth sepanjang 174 mengindikasikan bahwa model menggunakan bandwidth adaptif, di mana setiap lokasi memiliki jumlah tetangga efektif yang dapat berbeda.

4.9 UTM

library(sp)
library(spgwr)

# 4.8.6.1 Mendefinisikan CRS Geographic (Longitude–Latitude, WGS84)
proj_longlat <- CRS("+proj=longlat +datum=WGS84")

# 4.8.6.2 Mendefinisikan CRS UTM
# Indonesia berada di beberapa zona UTM,
# di sini digunakan contoh UTM Zona 50S (umum Indonesia bagian tengah)
proj_utm <- CRS("+proj=utm +zone=50 +south +datum=WGS84")

# 4.8.6.3 Menetapkan CRS pada data spasial
proj4string(data_spasial) <- proj_longlat

# 4.8.6.4 Mengonversi data spasial ke sistem koordinat UTM
data_spasial_utm <- spTransform(data_spasial, proj_utm)

# 4.8.6.5 Menentukan bandwidth optimal (adaptive) menggunakan data UTM
bandwidth_optimal <- gwr.sel(
  TPT ~ Lama.Sekolah + TPAK,
  data  = data_spasial_utm,
  adapt = TRUE
)

## Adaptive q: 0.381966 CV score: 185.8897 
## Adaptive q: 0.618034 CV score: 205.2481 
## Adaptive q: 0.236068 CV score: 156.4924 
## Adaptive q: 0.145898 CV score: 127.8438 
## Adaptive q: 0.09016994 CV score: 98.01613 
## Adaptive q: 0.05572809 CV score: 88.98689 
## Adaptive q: 0.02966668 CV score: 5012.997 
## Adaptive q: 0.06888371 CV score: 97.88722 
## Adaptive q: 0.04577352 CV score: 88.98689 
## Adaptive q: 0.0507508 CV score: 88.98689 
## Adaptive q: 0.05323945 CV score: 88.98689 
## Adaptive q: 0.04884965 CV score: 88.98689 
## Adaptive q: 0.05170138 CV score: 88.98689 
## Adaptive q: 0.05002463 CV score: 88.98689 
## Adaptive q: 0.04957583 CV score: 88.98689 
## Adaptive q: 0.050302 CV score: 88.98689 
## Adaptive q: 0.05038772 CV score: 88.98689 
## Adaptive q: 0.05021142 CV score: 88.98689 
## Adaptive q: 0.05026131 CV score: 88.98689 
## Adaptive q: 0.05034269 CV score: 88.98689 
## Adaptive q: 0.050302 CV score: 88.98689

bandwidth_optimal

## [1] 0.050302

# 4.8.6.6 Menjalankan model GWR dengan data terproyeksi UTM
model_gwr <- gwr(
  TPT ~ Lama.Sekolah + TPAK,
  data   = data_spasial_utm,
  adapt  = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit    = TRUE
)

# Ringkasan hasil GWR
summary(model_gwr)

##           Length Class                  Mode     
## SDF         174  SpatialPointsDataFrame S4       
## lhat      30276  -none-                 numeric  
## lm           11  -none-                 list     
## results      14  -none-                 list     
## bandwidth   174  -none-                 numeric  
## adapt         1  -none-                 numeric  
## hatmatrix     1  -none-                 logical  
## gweight       1  -none-                 character
## gTSS          1  -none-                 numeric  
## this.call     6  -none-                 call     
## fp.given      1  -none-                 logical  
## timings      12  -none-                 numeric

Hasil cross-validation pada model Geographically Weighted Regression (GWR) dengan sistem koordinat UTM menunjukkan bahwa bandwidth adaptif optimal diperoleh pada nilai q = 0,050302, yang menghasilkan CV score minimum sebesar 88,98689. Hal ini mengindikasikan bahwa sekitar 5% tetangga terdekat dari setiap lokasi digunakan dalam proses estimasi koefisien lokal, sehingga model mampu menangkap variasi spasial secara lebih sensitif pada skala lokal.

4.10 Uji Asumsi

4.10.1 Uji Heterogenitas

model_ols <- lm(
  TPT ~ Lama.Sekolah + TPAK,
  data = data_spasial
)
library(lmtest)
bptest(model_ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 2.9348, df = 2, p-value = 0.2305

p-value > 0,05 → tidak ada heteroskedastisitas global

Hasil uji studentized Breusch–Pagan pada model OLS menghasilkan nilai p-value sebesar 0,2305. Karena nilai p-value lebih besar dari taraf signifikansi 5% (0,05), maka gagal menolak H₀, yang berarti tidak terdapat heteroskedastisitas secara global pada model OLS. Dengan demikian, varians residual model OLS dapat dianggap homogen secara global.

4.10.2 Perbandingan AIC

AIC(model_ols)

## [1] 540.7293

# AIC corrected (yang biasanya dipakai di GWR)
model_gwr$results$AICc

## [1] 384.2586

nilai AICc GWR sebesar 384,26, yang secara signifikan lebih kecil dibandingkan AIC OLS sebesar 540,73. Nilai AIC/AICc yang lebih rendah mengindikasikan bahwa model GWR mampu menjelaskan variasi Tingkat Pengangguran Terbuka (TPT) dengan lebih baik

4.10.3 UJI MULTIKOLINEARITAS

library(lmtest)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

library(spdep)
library(ggplot2)
library(sp)
vif(model_ols)

## Lama.Sekolah         TPAK 
##     1.182095     1.182095

Nilai Variance Inflation Factor (VIF) untuk variabel Lama.Sekolah dan TPAK masing-masing sebesar 1,18, yang jauh di bawah batas umum VIF < 10. Hal ini menunjukkan bahwa tidak terdapat masalah multikolinearitas antar variabel independen dalam model.

4.10.4 Uji Spatial Autokorelasi

# Ekstraksi residual GWR
residual_gwr <- model_gwr$SDF$gwr.e

# Matriks bobot spasial (KNN, k = 4)
coords <- cbind(
  coordinates(data_spasial)[,1],
  coordinates(data_spasial)[,2]
)

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

## Warning in knearneigh(coords, k = 4): knearneigh: identical points found

## Warning in knearneigh(coords, k = 4): knearneigh: kd_tree not available for
## identical points

## Warning in knn2nb(knearneigh(coords, k = 4)): neighbour object has 29
## sub-graphs

listw <- nb2listw(nb, style = "W")

# Uji Moran’s I
moran.test(residual_gwr, listw)

## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw    
## 
## Moran I statistic standard deviate = 1.1933, p-value = 0.1164
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.054435477      -0.005780347       0.002546507

Nilai statistik Moran’s I sebesar 0,0569 dengan p-value 0,1164, yang lebih besar dari taraf signifikansi 5%. Dengan demikian, gagal menolak H₀, yang menyatakan bahwa tidak terdapat autokorelasi spasial pada residual model GWR. Artinya, sisa kesalahan model tidak membentuk pola spasial tertentu dan tersebar secara acak di wilayah studi.

4.10.5 NORMALITAS RESIDUAL

shapiro.test(residual_gwr)

## 
##  Shapiro-Wilk normality test
## 
## data:  residual_gwr
## W = 0.9791, p-value = 0.01012

Hasil uji Shapiro–Wilk menunjukkan bahwa p-value < 0,05, sehingga H₀ ditolak. Hal ini mengindikasikan bahwa residual model GWR tidak berdistribusi normal. Namun, pada konteks Geographically Weighted Regression, asumsi normalitas residual bukan merupakan syarat utama, karena GWR lebih menekankan pada variasi lokal dan pola spasial dibandingkan inferensi statistik global.

4.11 PLOT SPASIAL (spplot)

spplot(
  model_gwr$SDF,
  "TPAK",
  main = "Koefisien GWR untuk TPAK",
  col.regions = terrain.colors(100)
)

spplot(
  model_gwr$SDF,
  "Lama.Sekolah",
  main = "Koefisien GWR untuk Lama.Sekolah",
  col.regions = terrain.colors(100)
)

coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
coef_PM   <- model_gwr$SDF$TPAK

summary(coef_Lama.Sekolah)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.74986  0.03133  0.51442  0.36411  0.68540  1.33792

summary(coef_PM)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.6164 -0.3168 -0.2100 -0.2424 -0.1280  0.1131

Visualisasi koefisien lokal GWR menunjukkan adanya heterogenitas spasial yang signifikan dalam pengaruh TPAK dan rata-rata lama sekolah terhadap TPT antarprovinsi di Indonesia. Koefisien TPAK berkisar antara -0,8166 hingga -0,4705 dengan dominasi nilai negatif, mengindikasikan bahwa peningkatan TPAK secara konsisten menurunkan tingkat pengangguran di sebagian besar wilayah, meskipun kekuatan pengaruhnya bervariasi secara spasial. Sementara itu, koefisien lama sekolah menunjukkan variasi yang lebih kompleks dengan rentang -0,7499 hingga 0,3329, di mana wilayah-wilayah tertentu memperlihatkan hubungan positif yang mengindikasikan potensi skill mismatch antara kualifikasi pendidikan dan kebutuhan pasar kerja lokal, sedangkan wilayah lain menunjukkan hubungan negatif yang mencerminkan bahwa peningkatan pendidikan efektif menurunkan pengangguran

4.11.1 VISUALISASI ggplot2

# Konversi SpatialPointsDataFrame ke data frame
gwr_df <- as.data.frame(model_gwr$SDF)

# Plot koefisien lokal Lama.Sekolah
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
  geom_point(aes(color = Lama.Sekolah), size = 2) +
  scale_color_viridis_c() +
  labs(
    title = "Koefisien Lokal GWR untuk Lama.Sekolah",
    color = "Koefisien"
  ) +
  theme_minimal()

# Plot koefisien lokal TPAK TPAK
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
  geom_point(aes(color = TPAK), size = 2) +
  scale_color_viridis_c() +
  labs(
    title = "Koefisien Lokal GWR untuk TPAK",
    color = "Koefisien"
  ) +
  theme_minimal()

Peta sebaran koefisien lokal GWR dalam sistem koordinat geografis memperlihatkan pola spasial yang tidak seragam, dengan koefisien lama sekolah yang didominasi nilai negatif (rata-rata -0,2224) namun memiliki beberapa wilayah dengan koefisien positif hingga 0,1889, serta koefisien TPAK yang memiliki nilai rata-rata -0,05874 dengan maksimum mencapai 0,24078. Distribusi spasial ini mengonfirmasi hipotesis bahwa pengaruh variabel independen terhadap TPT tidak bersifat stasioner secara spasial, melainkan bervariasi tergantung pada karakteristik lokal seperti struktur ekonomi regional, ketersediaan lapangan kerja, dan kesesuaian antara sistem pendidikan dengan kebutuhan industri di masing-masing provinsi, sehingga pendekatan model global seperti OLS atau bahkan model spasial global (SAR/SEM) tidak sepenuhnya mampu menangkap kompleksitas hubungan ini tanpa dilengkapi dengan analisis lokal seperti GWR.

4.12 Membuat Surface Plot (Interpolasi Koefisien GWR)

Surface plot dibuat menggunakan metode Inverse Distance Weighting (IDW) untuk memvisualisasikan pola spasial koefisien lokal hasil GWR.

4.12.1 Inverse Distance Weighting (IDW)

#4.11.1 Memastikan CRS Data (WGS84)
# Cek sistem koordinat
proj4string(data_spasial)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

## "+proj=longlat +datum=WGS84 +no_defs"

# Pastikan data dalam WGS84
data_spasial_wgs84 <- spTransform(
  data_spasial,
  CRS("+proj=longlat +datum=WGS84 +no_defs")
)

#4.11.2 Menentukan Rentang Koordinat

x.range <- range(data_spasial_wgs84$longitude)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

y.range <- range(data_spasial_wgs84$latitude)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf
## Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

#4.11.3 Membuat Grid Prediksi
# Membuat grid titik prediksi
x.range <- range(coords[,1], na.rm = TRUE)
y.range <- range(coords[,2], na.rm = TRUE)

grd <- expand.grid(
  longitude = seq(from = x.range[1], to = x.range[2], length.out = 100),
  latitude  = seq(from = y.range[1], to = y.range[2], length.out = 100)
)

coordinates(grd) <- ~longitude + latitude
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

#4.11.4 Menggabungkan Koefisien Lokal GWR ke Data Spasial
# Ekstrak koefisien lokal dari model GWR
data_spasial_wgs84$coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
data_spasial_wgs84$coef_PM   <- model_gwr$SDF$TPAK

#4.11.5 IDW Koefisien GWR – Variabel Lama.Sekolah
library(gstat)
class(data_spasial_wgs84)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

data_spasial_wgs84@data$coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
data_spasial_wgs84@data$coef_PM  <- model_gwr$SDF$TPAK

library(gstat)

idw_Lama.Sekolah <- gstat::idw(
  formula   = coef_Lama.Sekolah ~ 1,
  locations = data_spasial_wgs84,
  newdata   = grd,
  idp       = 2.0
)

## [inverse distance weighted interpolation]

idw_PM <- gstat::idw(
  formula   = coef_PM ~ 1,
  locations = data_spasial_wgs84,
  newdata   = grd,
  idp       = 2.0
)

## [inverse distance weighted interpolation]

idw_Lama.Sekolah_df <- as.data.frame(idw_Lama.Sekolah)

ggplot() +
  geom_tile(data = idw_Lama.Sekolah_df,
            aes(x = longitude, y = latitude, fill = var1.pred)) +
  scale_fill_viridis_c() +
  labs(title = "Surface Plot Koefisien GWR Lama.Sekolah (IDW)",
       fill = "Koefisien Lama.Sekolah") +
  theme_minimal()

idw_PM_df <- as.data.frame(idw_PM)

ggplot() +
  geom_tile(
    data = idw_PM_df,
    aes(x = longitude, y = latitude, fill = var1.pred)
  ) +
  scale_fill_viridis_c() +
  labs(
    title = "Surface Plot Koefisien GWR untuk TPAK TPAK (IDW)",
    fill = "Koefisien Pend. TPAK"
  ) +
  theme_minimal()

Surface plot hasil interpolasi Inverse Distance Weighted (IDW) untuk koefisien lokal GWR menampilkan permukaan kontinu yang memvisualisasikan gradien pengaruh variabel independen terhadap TPT di seluruh wilayah studi. Pola interpolasi menunjukkan adanya zona-zona dengan koefisien yang relatif homogen serta area transisi yang menandakan perubahan bertahap dalam kekuatan dan arah pengaruh variabel penjelas, khususnya terlihat pada wilayah dengan koefisien lama sekolah yang berkisar dari -0,5 (biru tua) hingga 1,0 (kuning), mengindikasikan bahwa wilayah bagian barat cenderung menunjukkan hubungan negatif yang lebih kuat dibandingkan wilayah timur Indonesia yang memperlihatkan koefisien mendekati nol atau bahkan positif, yang dapat dijelaskan oleh perbedaan tingkat pembangunan ekonomi, diversifikasi sektor industri, dan kualitas sistem pendidikan vokasional antarwilayah.

4.13 Spatial Interpolation

4.13.1 Nearest Neighbor (NN)

# =========================================================
# NEAREST NEIGHBOR (1-NN) UNTUK KOEFISIEN LOKAL GWR
# =========================================================

library(sp)
library(FNN)

## Warning: package 'FNN' was built under R version 4.3.3

library(deldir)

## Warning: package 'deldir' was built under R version 4.3.3

## deldir 2.0-4      Nickname: "Idol Comparison"

## 
##      The syntax of deldir() has changed since version 
##      0.0-10.  Read the help!!!.

# 1. Pastikan koefisien lokal GWR ada di data spasial
data_spasial@data$coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
data_spasial@data$coef_PM  <- model_gwr$SDF$TPAK

# 2. Ambil koordinat titik observasi dan grid prediksi
coords_obs  <- coordinates(data_spasial)
coords_grid <- coordinates(grd)

# 3. Cari tetangga terdekat (1-NN) untuk setiap titik grid
knn_res <- get.knnx(
  data  = coords_obs,
  query = coords_grid,
  k = 1
)

nn_index <- knn_res$nn.index[, 1]

# 4. Prediksi Nearest Neighbor (ambil nilai titik terdekat)
pred_nn_Lama.Sekolah <- data_spasial@data$coef_Lama.Sekolah[nn_index]
pred_nn_PM   <- data_spasial@data$coef_PM[nn_index]

# 5. Bangun SpatialPixelsDataFrame hasil NN
nn_grid_Lama.Sekolah <- SpatialPixelsDataFrame(
  points = grd,
  data   = data.frame(coef_Lama.Sekolah_NN = pred_nn_Lama.Sekolah),
  proj4string = CRS(proj4string(data_spasial))
)

nn_grid_PM <- SpatialPixelsDataFrame(
  points = grd,
  data   = data.frame(coef_PM_NN = pred_nn_PM),
  proj4string = CRS(proj4string(data_spasial))
)

# =========================================================
# VORONOI (APROKSIMASI NEAREST NEIGHBOR)
# =========================================================

vor <- deldir(coords_obs[,1], coords_obs[,2])

plot(
  vor,
  wlines  = "tess",
  wpoints = "none",
  main    = "Peta Nearest Neighbor (1-NN) — Aproksimasi Voronoi pada Grid"
)

## Warning in (function (x0, y0, x1 = x0, y1 = y0, col = par("fg"), lty =
## par("lty"), : "wpoints" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "wpoints" is not a
## graphical parameter

points(coords_obs, pch = 19, col = "red")

# =========================================================
# VISUALISASI HASIL NEAREST NEIGHBOR
# =========================================================

spplot(
  nn_grid_Lama.Sekolah,
  "coef_Lama.Sekolah_NN",
  main = "Nearest Neighbor (1-NN) Koefisien GWR Lama.Sekolah"
)

spplot(
  nn_grid_PM,
  "coef_PM_NN",
  main = "Nearest Neighbor (1-NN) Koefisien GWR TPAK TPAK"
)

Peta Nearest Neighbor (1-NN) yang divisualisasikan melalui aproksimasi diagram Voronoi memperlihatkan pembagian wilayah studi menjadi poligon-poligon diskret di mana setiap lokasi dalam poligon memiliki nilai koefisien yang sama dengan titik observasi terdekat, menghasilkan pola mozaik yang tajam tanpa transisi gradual. Metode ini, meskipun sederhana dan komputasional efisien, menghasilkan diskontinuitas spasial yang kurang realistis untuk fenomena geografis yang umumnya bersifat kontinu seperti TPT, namun tetap berguna sebagai baseline comparison untuk mengevaluasi kinerja metode interpolasi yang lebih kompleks seperti IDW, TIN, atau Kriging dalam menangkap variasi spasial lokal koefisien GWR.

4.13.2 Natural Neighbor / TIN

# =========================================================
# TIN (LINEAR / BARYCENTRIC) UNTUK KOEFISIEN LOKAL GWR
# (Lama.Sekolah & TPAK TPAK)
# =========================================================

library(sp)
library(interp)

## Warning: package 'interp' was built under R version 4.3.3

# ---------------------------------------------------------
# 1. Pastikan data_spasial sudah SpatialPointsDataFrame
# ---------------------------------------------------------
proj4string(data_spasial) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

# ---------------------------------------------------------
# 2. Tambahkan koefisien lokal GWR ke data spasial
# ---------------------------------------------------------
data_spasial@data$coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
data_spasial@data$coef_PM  <- model_gwr$SDF$TPAK

# ---------------------------------------------------------
# 3. Ambil koordinat titik observasi
# ---------------------------------------------------------
coords_obs <- coordinates(data_spasial)

# ---------------------------------------------------------
# 4. Buat grid prediksi
# ---------------------------------------------------------
x.range <- range(coords_obs[,1], na.rm = TRUE)
y.range <- range(coords_obs[,2], na.rm = TRUE)

xg <- seq(from = x.range[1], to = x.range[2], length.out = 100)
yg <- seq(from = y.range[1], to = y.range[2], length.out = 100)

# ---------------------------------------------------------
# 5. TIN untuk koefisien GWR Lama.Sekolah
# ---------------------------------------------------------
tin_Lama.Sekolah <- interp::interp(
  x  = coords_obs[,1],
  y  = coords_obs[,2],
  z  = data_spasial@data$coef_Lama.Sekolah,
  xo = xg,
  yo = yg,
  duplicate = "mean",
  linear = TRUE
)

grd_Lama.Sekolah <- expand.grid(longitude = tin_Lama.Sekolah$x,
                        latitude  = tin_Lama.Sekolah$y)

tin_grid_Lama.Sekolah <- SpatialPixelsDataFrame(
  points = SpatialPoints(grd_Lama.Sekolah),
  data   = data.frame(coef_Lama.Sekolah_TIN = as.vector(tin_Lama.Sekolah$z)),
  tolerance = 0.01
)

proj4string(tin_grid_Lama.Sekolah) <- CRS(proj4string(data_spasial))

# ---------------------------------------------------------
# 6. TIN untuk koefisien GWR TPAK TPAK
# ---------------------------------------------------------
tin_PM <- interp::interp(
  x  = coords_obs[,1],
  y  = coords_obs[,2],
  z  = data_spasial@data$coef_PM,
  xo = xg,
  yo = yg,
  duplicate = "mean",
  linear = TRUE
)

grd_PM <- expand.grid(longitude = tin_PM$x,
                      latitude  = tin_PM$y)

tin_grid_PM <- SpatialPixelsDataFrame(
  points = SpatialPoints(grd_PM),
  data   = data.frame(coef_PM_TIN = as.vector(tin_PM$z)),
  tolerance = 0.01
)

proj4string(tin_grid_PM) <- CRS(proj4string(data_spasial))

# ---------------------------------------------------------
# 7. Visualisasi Surface TIN
# ---------------------------------------------------------
spplot(
  tin_grid_Lama.Sekolah,
  "coef_Lama.Sekolah_TIN",
  main = "TIN (Linear / Barycentric) — Koefisien Lokal GWR Lama.Sekolah"
)

spplot(
  tin_grid_PM,
  "coef_PM_TIN",
  main = "TIN (Linear / Barycentric) — Koefisien Lokal GWR TPAK TPAK"
)

Interpolasi Triangulated Irregular Network (TIN) dengan metode linear/barycentric menghasilkan permukaan yang terdiri dari segitiga-segitiga yang menghubungkan titik-titik observasi, di mana nilai di dalam setiap segitiga dihitung melalui interpolasi linier dari tiga titik vertex terdekat. Visualisasi menunjukkan permukaan yang lebih halus dibandingkan Nearest Neighbor namun masih mempertahankan diskontinuitas pada batas-batas segitiga, dengan koefisien lama sekolah dan TPAK yang menampilkan gradasi warna dari biru tua (nilai negatif kuat) hingga kuning (nilai positif atau mendekati nol), mengindikasikan bahwa metode TIN mampu menangkap variasi lokal dengan lebih baik sambil menghormati nilai observasi aktual di lokasi sampel tanpa melakukan smoothing berlebihan seperti yang terjadi pada metode berbasis jarak global.

4.13.3 Inverse Distance Weighting (IDW)

# =========================================================
# IDW UNTUK KOEFISIEN LOKAL GWR (Lama.Sekolah & PEND. TPAK)
# =========================================================

library(sp)
library(gstat)

# 1. Pastikan CRS sama (WGS84 / sesuai data kamu)
proj4string(data_spasial)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# 2. Tambahkan koefisien lokal GWR ke data spasial
data_spasial@data$coef_Lama.Sekolah <- model_gwr$SDF$Lama.Sekolah
data_spasial@data$coef_PM  <- model_gwr$SDF$TPAK

# 3. Buat grid prediksi dari extent data
coords <- coordinates(data_spasial)

x.range <- range(coords[,1], na.rm = TRUE)
y.range <- range(coords[,2], na.rm = TRUE)

grd <- expand.grid(
  x = seq(x.range[1], x.range[2], length.out = 100),
  y = seq(y.range[1], y.range[2], length.out = 100)
)

coordinates(grd) <- ~x+y
proj4string(grd) <- CRS(proj4string(data_spasial))

# 4. IDW untuk koefisien Lama.Sekolah
idw_Lama.Sekolah <- idw(
  coef_Lama.Sekolah ~ 1,
  locations = data_spasial,
  newdata   = grd,
  idp       = 2
)

## [inverse distance weighted interpolation]

# 5. IDW untuk koefisien TPAK TPAK
idw_PM <- idw(
  coef_PM ~ 1,
  locations = data_spasial,
  newdata   = grd,
  idp       = 2
)

## [inverse distance weighted interpolation]

# 6. Visualisasi
spplot(
  idw_Lama.Sekolah["var1.pred"],
  main = "IDW (p=2) — Permukaan Koefisien GWR Lama.Sekolah"
)

spplot(
  idw_PM["var1.pred"],
  main = "IDW (p=2) — Permukaan Koefisien GWR TPAK TPAK"
)

Hasil interpolasi Inverse Distance Weighted dengan parameter pembobot p=2 menunjukkan permukaan yang lebih halus dibandingkan NN dan TIN, dengan transisi gradual antarwilayah yang mencerminkan asumsi bahwa pengaruh setiap titik observasi menurun secara kuadratik terhadap jarak. Pola spasial memperlihatkan bull’s-eye effect di sekitar lokasi observasi dengan nilai ekstrem (ditandai warna biru tua untuk koefisien negatif kuat dan kuning untuk positif), yang merupakan karakteristik inherent dari metode IDW karena pembobotan yang sangat tinggi pada titik terdekat, dan visualisasi ini mengonfirmasi bahwa koefisien lama sekolah memiliki variasi spasial yang lebih besar (rentang -0,7444 hingga 0,3294) dibandingkan TPAK, yang konsisten dengan temuan model GWR dan mengindikasikan bahwa hubungan antara pendidikan dan pengangguran lebih kompleks dan context-dependent secara geografis dibandingkan hubungan antara kemiskinan dan pengangguran.

4.13.4 Trend Surface / Polynomial Regression

# =========================================================
# THIN PLATE SPLINES (TPS) — VARIABEL TPT
# =========================================================

library(sp)
library(fields)

## Warning: package 'fields' was built under R version 4.3.3

## Loading required package: spam

## Warning: package 'spam' was built under R version 4.3.3

## Spam version 2.11-1 (2025-01-20) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following object is masked from 'package:Matrix':
## 
##     det

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: viridisLite

## 
## Try help(fields) to get started.

## 
## Attaching package: 'fields'

## The following object is masked from 'package:psych':
## 
##     describe

# ---------------------------------------------------------
# 1. Ambil koordinat spasial (longitude, latitude)
# ---------------------------------------------------------
coords <- coordinates(data_spasial)

# ---------------------------------------------------------
# 2. Ambil variabel respon (Y = TPT)
# ---------------------------------------------------------
Yvals <- data_spasial@data$TPT

# Cek keamanan
stopifnot(
  is.matrix(coords),
  is.numeric(Yvals),
  length(Yvals) == nrow(coords)
)

# ---------------------------------------------------------
# 3. Fit Thin Plate Splines (lambda otomatis via GCV)
# ---------------------------------------------------------
tps_fit <- fields::Tps(
  x = coords,
  Y = Yvals
)

summary(tps_fit)

## CALL:
## fields::Tps(x = coords, Y = Yvals)
##                                               
##  Number of Observations:                174   
##  Number of unique points:               29    
##  Number of parameters in the null space 3     
##  Parameters for fixed spatial drift     3     
##  Effective degrees of freedom:          11.6  
##  Residual degrees of freedom:           162.4 
##  MLE tau                                0.9341
##  GCV tau                                1.012 
##  Pure error tau                         0.6593
##  MLE sigma                              52.56 
##  Scale passed for covariance (sigma)    <NA>  
##  Scale passed for nugget (tau^2)        <NA>  
##  Smoothing parameter lambda             0.0166
## 
## Residual Summary:
##      min    1st Q   median    3rd Q      max 
## -2.64500 -0.69650 -0.09194  0.65970  2.66600 
## 
## Covariance Model: Rad.cov
##   Names of non-default covariance arguments: 
##        p
## 
## DETAILS ON SMOOTHING PARAMETER:
##  Method used:   GCV    Cost:  1
##    lambda       trA       GCV   GCV.one GCV.model    tauHat 
##    0.0166   11.6181   11.1291    1.0978   10.3539    1.0122 
## 
##  Summary of all estimates found for lambda
##               lambda   trA      GCV tauHat -lnLike Prof converge
## GCV        0.0166012 11.62  11.1291 1.0122        67.40        3
## GCV.model  0.0082165 14.38  10.1576 0.9186           NA        3
## GCV.one    0.0001587 27.55   0.5193 0.6611           NA       NA
## RMSE              NA    NA       NA     NA           NA       NA
## pure error 0.0001587 27.55 187.4814 0.6611        72.35       NA
## REML       0.0240439 10.33  11.1851 1.0629        67.37        1

# ---------------------------------------------------------
# 4. Buat grid prediksi
# ---------------------------------------------------------
x.range <- range(coords[,1])
y.range <- range(coords[,2])

grd <- expand.grid(
  longitude = seq(x.range[1], x.range[2], length.out = 100),
  latitude  = seq(y.range[1], y.range[2], length.out = 100)
)

coordinates(grd) <- ~longitude + latitude
proj4string(grd) <- CRS(proj4string(data_spasial))

# ---------------------------------------------------------
# 5. Prediksi TPS di grid
# ---------------------------------------------------------
pred_tps <- predict(tps_fit, coordinates(grd))

# ---------------------------------------------------------
# 6. Bangun SpatialPixelsDataFrame
# ---------------------------------------------------------
tps_grid <- SpatialPixelsDataFrame(
  points = grd,
  data   = data.frame(TPT_TPS = as.numeric(pred_tps)),
  proj4string = CRS(proj4string(data_spasial))
)

# ---------------------------------------------------------
# 7. Visualisasi
# ---------------------------------------------------------
spplot(
  tps_grid,
  "TPT_TPS",
  main = "Thin Plate Splines (TPS) — Permukaan TPT",
  col.regions = terrain.colors(100)
)

Permukaan Thin Plate Splines (TPS) untuk variabel TPT menghasilkan interpolasi yang sangat halus dan kontinu di seluruh wilayah studi, dengan pola yang mencerminkan tren global dan variasi lokal secara bersamaan melalui minimasi fungsi energi bending. Visualisasi menunjukkan gradasi nilai TPT dari sekitar 3,5% (hijau) hingga 6,5% (merah-oranye) dengan permukaan yang tidak memiliki diskontinuitas tajam, mengindikasikan bahwa metode TPS berhasil menangkap struktur spasial TPT yang merefleksikan fenomena geografis berkelanjutan, di mana wilayah dengan TPT tinggi (ditandai warna merah-oranye) terkonsentrasi pada area tertentu sementara wilayah dengan TPT rendah (hijau) tersebar di bagian lain, dan smoothness yang dihasilkan TPS menjadikannya metode yang sesuai untuk visualisasi pola makro dan identifikasi hotspot pengangguran dalam konteks perencanaan kebijakan regional.

4.13.5 Kriging

library(sp)
library(gstat)

# ===============================
# 0. Pembersihan data
# ===============================
data_spasial$TPT <- as.numeric(as.character(data_spasial$TPT))
data_spasial <- data_spasial[!is.na(data_spasial$TPT), ]
data_spasial <- data_spasial[!duplicated(coordinates(data_spasial)), ]

# ===============================
# 1. Variogram
# ===============================
vgm_emp <- variogram(TPT ~ 1, data_spasial)

vgm_model <- fit.variogram(
  vgm_emp,
  model = vgm(
    psill = var(data_spasial$TPT),
    model = "Sph",
    range = max(vgm_emp$dist),
    nugget = 0
  ),
  fit.method = 2
)

## Warning in fit.variogram(vgm_emp, model = vgm(psill = var(data_spasial$TPT), :
## No convergence after 200 iterations: try different initial values?

# ===============================
# 2. Grid dari bbox
# ===============================
bb <- bbox(data_spasial)

grd <- expand.grid(
  x = seq(bb[1,1], bb[1,2], length.out = 80),
  y = seq(bb[2,1], bb[2,2], length.out = 80)
)

coordinates(grd) <- ~x+y
proj4string(grd) <- CRS(proj4string(data_spasial))
gridded(grd) <- TRUE

# ===============================
# 3. Ordinary Kriging
# ===============================
ok <- gstat(
  formula = TPT ~ 1,
  data = data_spasial,
  model = vgm_model
)

krig_res <- predict(ok, grd)

## [using ordinary kriging]

# ===============================
# 4. PLOT
# ===============================
summary(krig_res$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.923   4.142   4.493   4.522   4.805   7.504

spplot(
  krig_res["var1.pred"],
  main = "Ordinary Kriging — Prediksi TPT"
)

spplot(
  krig_res["var1.var"],
  main = "Ordinary Kriging — Varians Prediksi TPT"
)

Peta hasil Ordinary Kriging menampilkan dua komponen penting: prediksi TPT yang menunjukkan pola spasial dengan nilai berkisar antara 2% hingga 7% dengan konsentrasi nilai tinggi di beberapa area (ditandai warna kuning), serta peta varians prediksi yang mengindikasikan tingkat ketidakpastian estimasi di setiap lokasi, di mana area dengan varians tinggi (merah-kuning) umumnya berada jauh dari titik observasi.

4.14 Implementasi Dasar (CRS-Safe)

# =========================================================
# ORDINARY KRIGING – DATA MEUSE (1 CHUNK LENGKAP & AMAN)
# =========================================================

library(sp)
library(gstat)

# ---------------------------------------------------------
# 1. Fungsi bantu: set CRS aman
# ---------------------------------------------------------
set_crs <- function(obj){
  ok <- TRUE
  tryCatch({
    proj4string(obj) <- CRS("+init=epsg:28992")
  }, error = function(e) ok <<- FALSE)
  if(!ok){
    proj4string(obj) <- CRS(SRS_string = "EPSG:28992")
  }
  obj
}

# ---------------------------------------------------------
# 2. Load & siapkan data meuse
# ---------------------------------------------------------
data(meuse, package = "sp")
coordinates(meuse) <- ~x+y
meuse <- set_crs(meuse)

## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.

# Variabel yang dipakai
meuse$log_zinc <- log(meuse$zinc)

# ---------------------------------------------------------
# 3. Grid prediksi
# ---------------------------------------------------------
data(meuse.grid, package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS(proj4string(meuse))

# ---------------------------------------------------------
# 4. Variogram empirik
# ---------------------------------------------------------
vgm_emp <- variogram(log_zinc ~ 1, meuse)
plot(vgm_emp, main = "Variogram Empirik log(Zinc)")

# ---------------------------------------------------------
# 5. Fitting variogram teoritis (Spherical)
# ---------------------------------------------------------
vgm_fit <- fit.variogram(
  vgm_emp,
  model = vgm(
    psill  = var(meuse$log_zinc),
    model  = "Sph",
    range  = max(vgm_emp$dist) / 2,
    nugget = 0
  ),
  fit.method = 2   # aman dari error zero-distance
)

print(vgm_fit)

##   model      psill    range
## 1   Nug 0.06629402   0.0000
## 2   Sph 0.57069534 918.8448

plot(vgm_emp, vgm_fit, main = "Variogram Empirik & Model Spherical")

# ---------------------------------------------------------
# 6. Ordinary Kriging
# ---------------------------------------------------------
krig_res <- krige(
  log_zinc ~ 1,
  meuse,
  meuse.grid,
  model = vgm_fit
)

## [using ordinary kriging]

# ---------------------------------------------------------
# 7. Peta hasil Kriging
# ---------------------------------------------------------
spplot(
  krig_res["var1.pred"],
  main = "Ordinary Kriging — Prediksi log(Zinc)",
  col.regions = terrain.colors(100)
)

spplot(
  krig_res["var1.var"],
  main = "Ordinary Kriging — Varian Prediksi",
  col.regions = heat.colors(100)
)

ariogram empirik log(Zinc) dari dataset Meuse dan model spherical yang di-fit menunjukkan struktur autokorelasi spasial yang jelas, di mana semivariance meningkat seiring dengan bertambahnya jarak (lag distance) hingga mencapai sill pada jarak sekitar 1000 meter, mengindikasikan bahwa observasi yang berdekatan memiliki nilai yang lebih mirip dibandingkan yang berjauhan

Model spherical berhasil menangkap pola ini dengan parameter nugget yang menggambarkan variasi skala mikro atau measurement error, partial sill yang merefleksikan variabilitas spasial terstruktur, dan range yang menunjukkan jarak maksimum di mana autokorelasi spasial masih signifikan, dan kecocokan model teoritis terhadap variogram empirik (ditunjukkan oleh kurva yang mengikuti sebaran titik-titik observasi) memvalidasi penggunaan model spherical untuk Ordinary Kriging dalam memprediksi nilai di lokasi yang tidak tersampel.

Hasil Ordinary Kriging untuk data log(Zinc) di Meuse menampilkan peta prediksi yang menunjukkan gradasi konsentrasi zinc dari rendah (hijau, sekitar 5,0) hingga tinggi (kuning-merah, hingga 7,5) dengan pola spasial yang kontinu dan smooth, serta peta varians prediksi yang mengidentifikasi area dengan ketidakpastian tinggi (merah, hingga 0,50) terutama di lokasi yang jauh dari titik sampling atau di tepi boundary wilayah studi.

4.15 Universal Kriging

library(sf)
library(gstat)
library(sp)

# --------------------------------------------------
# 1. DATA sf UTM (BERSIH & SIAP UK)
# --------------------------------------------------
data_spasial <- data_spasial[!is.na(data_spasial$TPT), ]

data_sf <- st_as_sf(data_spasial)

if (is.na(st_crs(data_sf))) {
  st_crs(data_sf) <- 4326
}

# Indonesia barat (ganti jika perlu)
data_sf_utm <- st_transform(data_sf, 32748)

# Tambahkan kovariat trend (x, y)
coords <- st_coordinates(data_sf_utm)
data_sf_utm$x <- coords[,1]
data_sf_utm$y <- coords[,2]

data_sp_utm <- as(data_sf_utm, "Spatial")

# --------------------------------------------------
# 2. VARIOGRAM UNIVERSAL KRIGING (AMAN dist = 0)
# --------------------------------------------------
vgm_emp_uk <- variogram(TPT ~ x + y, data_sp_utm)

vgm_model_uk <- fit.variogram(
  vgm_emp_uk,
  vgm("Sph"),
  fit.method = 2
)

## Warning in fit.variogram(vgm_emp_uk, vgm("Sph"), fit.method = 2): No
## convergence after 200 iterations: try different initial values?

print(vgm_model_uk)

##   model    psill    range
## 1   Nug 1.088482      0.0
## 2   Sph 1.811855 913029.7

# --------------------------------------------------
# 3. GRID PREDIKSI + TREND
# --------------------------------------------------
bb <- st_bbox(data_sf_utm)

grd <- expand.grid(
  x = seq(bb["xmin"], bb["xmax"], length.out = 70),
  y = seq(bb["ymin"], bb["ymax"], length.out = 70)
)

coordinates(grd) <- ~ x + y
proj4string(grd) <- CRS(st_crs(data_sf_utm)$proj4string)
gridded(grd) <- TRUE

# --------------------------------------------------
# 4. UNIVERSAL KRIGING (KUNCI: nmax + maxdist)
# --------------------------------------------------
uk_model <- gstat(
  formula = TPT ~ x + y,
  data    = data_sp_utm,
  model   = vgm_model_uk,
  nmax    = 30,
  maxdist = Inf
)

uk_res <- predict(uk_model, newdata = grd)

## [using universal kriging]

# --------------------------------------------------
# 5. CEK HASIL (WAJIB FINITE)
# --------------------------------------------------
summary(uk_res$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.343   3.924   4.657   4.595   5.245   6.591

# --------------------------------------------------
# 6. PLOT AMAN (ANTI zrng ERROR)
# --------------------------------------------------
spplot(
  uk_res["var1.pred"],
  main = "Universal Kriging (x + y) — Prediksi Tingkat Pengangguran Terbuka (TPT)",
  col.regions = terrain.colors(100),
  at = pretty(
    range(uk_res$var1.pred, na.rm = TRUE),
    n = 100
  )
)

spplot(
  uk_res["var1.var"],
  main = "Universal Kriging — Varians Prediksi TPT",
  col.regions = heat.colors(100)
)

Peta hasil Universal Kriging untuk TPT dengan kovariat tren spasial (koordinat x dan y) menampilkan prediksi yang lebih halus dibandingkan Ordinary Kriging karena mengakomodasi drift atau tren deterministik dalam data, dengan nilai TPT yang diprediksi berkisar antara 2,5% hingga 6,5% dan memperlihatkan gradasi regional yang mencerminkan struktur ekonomi makro yang bervariasi antarwilayah. Peta varians prediksi menunjukkan ketidakpastian yang relatif rendah di sebagian besar wilayah (warna kuning, sekitar 1,0-2,0)

4.16 Co-Kriging (Multivariat) — Linear Model of Coregionalization (LMC)

# ==================================================
# CO-KRIGING (ICM) — TPT & Lama Sekolah
# ==================================================

library(sf)
library(sp)
library(gstat)

# --------------------------------------------------
# 1. DATA SPASIAL (BERSIH & UTM)
# --------------------------------------------------
data_spasial <- data_spasial[
  !is.na(data_spasial$TPT) & !is.na(data_spasial$Lama.Sekolah),
]

data_sf <- st_as_sf(data_spasial)

# set CRS bila belum ada
if (is.na(st_crs(data_sf))) {
  st_crs(data_sf) <- 4326
}

# UTM Indonesia Barat (sesuaikan zona jika perlu)
data_sf_utm <- st_transform(data_sf, 32748)
data_sp_utm <- as(data_sf_utm, "Spatial")

# --------------------------------------------------
# 2. VARIOGRAM MASING-MASING VARIABEL
# --------------------------------------------------
vg_TPT <- variogram(TPT ~ 1, data_sp_utm)
vg_LS  <- variogram(Lama.Sekolah ~ 1, data_sp_utm)

# fit model variogram (Spherical = stabil)
m_TPT <- fit.variogram(vg_TPT, vgm("Sph"))
m_LS  <- fit.variogram(vg_LS,  vgm("Sph"))

# --------------------------------------------------
# 3. DEFINISI CO-KRIGING (ICM)
# --------------------------------------------------
gs <- gstat(
  NULL,
  id      = "TPT",
  formula = TPT ~ 1,
  data    = data_sp_utm,
  model   = m_TPT
)

gs <- gstat(
  gs,
  id      = "Lama.Sekolah",
  formula = Lama.Sekolah ~ 1,
  data    = data_sp_utm,
  model   = m_LS
)

# --------------------------------------------------
# 4. GRID PREDIKSI
# --------------------------------------------------
bb <- st_bbox(data_sf_utm)

grd <- expand.grid(
  x = seq(bb["xmin"], bb["xmax"], length.out = 70),
  y = seq(bb["ymin"], bb["ymax"], length.out = 70)
)

coordinates(grd) <- ~ x + y
proj4string(grd) <- CRS(st_crs(data_sf_utm)$proj4string)
gridded(grd) <- TRUE

# --------------------------------------------------
# 5. CO-KRIGING (ORDINARY, ICM)
# --------------------------------------------------
ck.pred <- predict(
  gs,
  newdata = grd,
  nmax    = 30
)

## [using ordinary kriging]

# --------------------------------------------------
# 6. PLOT HASIL
# --------------------------------------------------
spplot(ck.pred["TPT.pred"], main = "Co-Kriging (ICM) — TPT")

spplot(ck.pred["Lama.Sekolah.pred"], main = "Co-Kriging (ICM) — Lama Sekolah")

Hasil Co-Kriging dengan Intrinsic Correlation Model (ICM) untuk TPT dan lama sekolah menampilkan peta prediksi yang memanfaatkan struktur korelasi silang (cross-covariance) antara kedua variabel untuk meningkatkan akurasi estimasi, khususnya di lokasi di mana salah satu variabel tersampel lebih padat dibandingkan yang lain. Pola spasial TPT hasil Co-Kriging memperlihatkan variasi dari sekitar 3,0% hingga 6,5% (biru hingga kuning) dengan permukaan yang lebih refined dibandingkan Ordinary Kriging univariat karena informasi tambahan dari lama sekolah

4.17 Validasi Silang

## 4.17 Validasi Silang
# =========================================================
# VALIDASI SILANG (CROSS-VALIDATION) - SEMUA METODE
# =========================================================

library(sp)
library(gstat)
library(FNN)
library(interp)
library(fields)
library(ggplot2)
library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.3.3

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

cat("========================================\n")

## ========================================

cat("VALIDASI SILANG - SEMUA METODE INTERPOLASI\n")

## VALIDASI SILANG - SEMUA METODE INTERPOLASI

cat("========================================\n\n")

## ========================================

# ---------------------------------------------------------
# 1. PERSIAPAN DATA
# ---------------------------------------------------------
# Pastikan data bersih
data_clean <- data_spasial[!duplicated(coordinates(data_spasial)), ]
data_clean <- data_clean[!is.na(data_clean$TPT), ]

cat("Jumlah observasi untuk validasi:", nrow(data_clean), "\n\n")

## Jumlah observasi untuk validasi: 29

# ---------------------------------------------------------
# 2. CROSS-VALIDATION: ORDINARY KRIGING
# ---------------------------------------------------------
cat("Running Cross-Validation: Ordinary Kriging...\n")

## Running Cross-Validation: Ordinary Kriging...

tryCatch({
  # Variogram
  vgm_emp <- variogram(TPT ~ 1, data_clean)
  vgm_model <- fit.variogram(
    vgm_emp,
    vgm(psill = var(data_clean$TPT), model = "Sph", 
        range = max(vgm_emp$dist)/2, nugget = 0),
    fit.method = 2
  )
  
  # Cross-validation
  cv_ok <- krige.cv(TPT ~ 1, data_clean, model = vgm_model, nfold = 10)
  
  # Metrik
  rmse_ok <- sqrt(mean(cv_ok$residual^2, na.rm = TRUE))
  mae_ok <- mean(abs(cv_ok$residual), na.rm = TRUE)
  me_ok <- mean(cv_ok$residual, na.rm = TRUE)
  var_ok <- mean(cv_ok$var1.var, na.rm = TRUE)
  
  cat("✓ Ordinary Kriging: RMSE =", round(rmse_ok, 4), "\n")
  
}, error = function(e) {
  cat("✗ Ordinary Kriging failed:", e$message, "\n")
  rmse_ok <<- NA
  mae_ok <<- NA
  me_ok <<- NA
  var_ok <<- NA
  cv_ok <<- NULL
})

## ✓ Ordinary Kriging: RMSE = 1.4756

# ---------------------------------------------------------
# 3. CROSS-VALIDATION: UNIVERSAL KRIGING
# ---------------------------------------------------------
cat("Running Cross-Validation: Universal Kriging...\n")

## Running Cross-Validation: Universal Kriging...

tryCatch({
  # Transform ke UTM
  data_sf_temp <- st_as_sf(data_clean)
  if(is.na(st_crs(data_sf_temp))) st_crs(data_sf_temp) <- 4326
  data_sf_utm <- st_transform(data_sf_temp, 32748)
  
  coords_utm <- st_coordinates(data_sf_utm)
  data_sf_utm$x <- coords_utm[,1]
  data_sf_utm$y <- coords_utm[,2]
  data_sp_utm <- as(data_sf_utm, "Spatial")
  
  # Variogram dengan trend
  vgm_emp_uk <- variogram(TPT ~ x + y, data_sp_utm)
  vgm_model_uk <- fit.variogram(vgm_emp_uk, vgm("Sph"), fit.method = 2)
  
  # Cross-validation
  cv_uk <- krige.cv(TPT ~ x + y, data_sp_utm, model = vgm_model_uk, nfold = 10)
  
  # Metrik
  rmse_uk <- sqrt(mean(cv_uk$residual^2, na.rm = TRUE))
  mae_uk <- mean(abs(cv_uk$residual), na.rm = TRUE)
  me_uk <- mean(cv_uk$residual, na.rm = TRUE)
  var_uk <- mean(cv_uk$var1.var, na.rm = TRUE)
  
  cat("✓ Universal Kriging: RMSE =", round(rmse_uk, 4), "\n")
  
}, error = function(e) {
  cat("✗ Universal Kriging failed:", e$message, "\n")
  rmse_uk <<- NA
  mae_uk <<- NA
  me_uk <<- NA
  var_uk <<- NA
  cv_uk <<- NULL
})

## Warning in fit.variogram(vgm_emp_uk, vgm("Sph"), fit.method = 2): No
## convergence after 200 iterations: try different initial values?

## ✓ Universal Kriging: RMSE = 1.6217

# ---------------------------------------------------------
# 4. CROSS-VALIDATION: CO-KRIGING
# ---------------------------------------------------------
cat("Running Cross-Validation: Co-Kriging...\n")

## Running Cross-Validation: Co-Kriging...

tryCatch({
  # Filter data dengan Lama.Sekolah
  data_ck <- data_clean[!is.na(data_clean$Lama.Sekolah), ]
  
  if(nrow(data_ck) > 20) {
    # Transform ke UTM
    data_sf_ck <- st_as_sf(data_ck)
    if(is.na(st_crs(data_sf_ck))) st_crs(data_sf_ck) <- 4326
    data_sf_ck_utm <- st_transform(data_sf_ck, 32748)
    data_sp_ck_utm <- as(data_sf_ck_utm, "Spatial")
    
    # Variogram untuk masing-masing variabel
    vg_TPT <- variogram(TPT ~ 1, data_sp_ck_utm)
    vg_LS <- variogram(Lama.Sekolah ~ 1, data_sp_ck_utm)
    m_TPT <- fit.variogram(vg_TPT, vgm("Sph"))
    m_LS <- fit.variogram(vg_LS, vgm("Sph"))
    
    # Manual cross-validation (simplified)
    pred_ck <- rep(NA, nrow(data_sp_ck_utm))
    for(i in 1:nrow(data_sp_ck_utm)) {
      train_data <- data_sp_ck_utm[-i, ]
      test_point <- data_sp_ck_utm[i, ]
      
      tryCatch({
        gs_temp <- gstat(NULL, id = "TPT", formula = TPT ~ 1, 
                        data = train_data, model = m_TPT)
        gs_temp <- gstat(gs_temp, id = "Lama.Sekolah", formula = Lama.Sekolah ~ 1, 
                        data = train_data, model = m_LS)
        pred_temp <- predict(gs_temp, test_point, nmax = 30)
        pred_ck[i] <- pred_temp$TPT.pred[1]
      }, error = function(e) {})
    }
    
    # Metrik
    residual_ck <- data_sp_ck_utm$TPT - pred_ck
    rmse_ck <- sqrt(mean(residual_ck^2, na.rm = TRUE))
    mae_ck <- mean(abs(residual_ck), na.rm = TRUE)
    me_ck <- mean(residual_ck, na.rm = TRUE)
    var_ck <- NA
    
    cat("✓ Co-Kriging: RMSE =", round(rmse_ck, 4), "\n")
  } else {
    cat("✗ Co-Kriging: Insufficient data\n")
    rmse_ck <- NA
    mae_ck <- NA
    me_ck <- NA
    var_ck <- NA
    pred_ck <- NULL
  }
  
}, error = function(e) {
  cat("✗ Co-Kriging failed:", e$message, "\n")
  rmse_ck <<- NA
  mae_ck <<- NA
  me_ck <<- NA
  var_ck <<- NA
  pred_ck <<- NULL
})

## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## ✓ Co-Kriging: RMSE = 1.4915

# ---------------------------------------------------------
# 5. CROSS-VALIDATION: NEAREST NEIGHBOR
# ---------------------------------------------------------
cat("Running Cross-Validation: Nearest Neighbor...\n")

## Running Cross-Validation: Nearest Neighbor...

coords_obs <- coordinates(data_clean)
pred_nn <- numeric(nrow(data_clean))

for(i in 1:nrow(data_clean)) {
  train_idx <- setdiff(1:nrow(data_clean), i)
  knn_res <- get.knnx(coords_obs[train_idx, ], matrix(coords_obs[i, ], nrow=1), k = 1)
  pred_nn[i] <- data_clean$TPT[train_idx[knn_res$nn.index[1]]]
}

residual_nn <- data_clean$TPT - pred_nn
rmse_nn <- sqrt(mean(residual_nn^2, na.rm = TRUE))
mae_nn <- mean(abs(residual_nn), na.rm = TRUE)
me_nn <- mean(residual_nn, na.rm = TRUE)

cat("✓ Nearest Neighbor: RMSE =", round(rmse_nn, 4), "\n")

## ✓ Nearest Neighbor: RMSE = 1.7562

# ---------------------------------------------------------
# 6. CROSS-VALIDATION: IDW
# ---------------------------------------------------------
cat("Running Cross-Validation: IDW...\n")

## Running Cross-Validation: IDW...

pred_idw <- numeric(nrow(data_clean))

for(i in 1:nrow(data_clean)) {
  train_data <- data_clean[-i, ]
  test_point <- data_clean[i, ]
  
  tryCatch({
    idw_res <- idw(TPT ~ 1, train_data, test_point, idp = 2)
    pred_idw[i] <- idw_res$var1.pred[1]
  }, error = function(e) {
    pred_idw[i] <<- NA
  })
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

residual_idw <- data_clean$TPT - pred_idw
rmse_idw <- sqrt(mean(residual_idw^2, na.rm = TRUE))
mae_idw <- mean(abs(residual_idw), na.rm = TRUE)
me_idw <- mean(residual_idw, na.rm = TRUE)

cat("✓ IDW: RMSE =", round(rmse_idw, 4), "\n")

## ✓ IDW: RMSE = 1.4597

# ---------------------------------------------------------
# 7. CROSS-VALIDATION: TIN (Linear Interpolation)
# ---------------------------------------------------------
cat("Running Cross-Validation: TIN (Linear)...\n")

## Running Cross-Validation: TIN (Linear)...

pred_tin <- numeric(nrow(data_clean))

for(i in 1:nrow(data_clean)) {
  train_data <- data_clean[-i, ]
  train_coords <- coordinates(train_data)
  
  tryCatch({
    tin_res <- interp(
      x = train_coords[,1],
      y = train_coords[,2],
      z = train_data$TPT,
      xo = coords_obs[i,1],
      yo = coords_obs[i,2],
      duplicate = "mean",
      linear = TRUE
    )
    pred_tin[i] <- tin_res$z[1]
  }, error = function(e) {
    pred_tin[i] <<- NA
  })
}

residual_tin <- data_clean$TPT - pred_tin
rmse_tin <- sqrt(mean(residual_tin^2, na.rm = TRUE))
mae_tin <- mean(abs(residual_tin), na.rm = TRUE)
me_tin <- mean(residual_tin, na.rm = TRUE)

cat("✓ TIN (Linear): RMSE =", round(rmse_tin, 4), "\n")

## ✓ TIN (Linear): RMSE = 1.419

# ---------------------------------------------------------
# 8. CROSS-VALIDATION: THIN PLATE SPLINES
# ---------------------------------------------------------
cat("Running Cross-Validation: Thin Plate Splines...\n")

## Running Cross-Validation: Thin Plate Splines...

pred_tps <- numeric(nrow(data_clean))

for(i in 1:nrow(data_clean)) {
  train_data <- data_clean[-i, ]
  train_coords <- coordinates(train_data)
  
  tryCatch({
    tps_fit <- Tps(x = train_coords, Y = train_data$TPT)
    pred_tps[i] <- predict(tps_fit, matrix(coords_obs[i,], nrow=1))
  }, error = function(e) {
    pred_tps[i] <<- NA
  })
}

residual_tps <- data_clean$TPT - pred_tps
rmse_tps <- sqrt(mean(residual_tps^2, na.rm = TRUE))
mae_tps <- mean(abs(residual_tps), na.rm = TRUE)
me_tps <- mean(residual_tps, na.rm = TRUE)

cat("✓ Thin Plate Splines: RMSE =", round(rmse_tps, 4), "\n\n")

## ✓ Thin Plate Splines: RMSE = 1.7062

# ---------------------------------------------------------
# 9. TABEL HASIL VALIDASI
# ---------------------------------------------------------
hasil_validasi <- data.frame(
  Metode = c("Ordinary Kriging", "Universal Kriging", "Co-Kriging",
             "Nearest Neighbor", "IDW", "TIN (Linear)", "Thin Plate Splines"),
  RMSE = c(rmse_ok, rmse_uk, rmse_ck, rmse_nn, rmse_idw, rmse_tin, rmse_tps),
  MAE = c(mae_ok, mae_uk, mae_ck, mae_nn, mae_idw, mae_tin, mae_tps),
  ME = c(me_ok, me_uk, me_ck, me_nn, me_idw, me_tin, me_tps),
  Varians_Kriging = c(var_ok, var_uk, var_ck, NA, NA, NA, NA)
)

# Urutkan berdasarkan RMSE
hasil_validasi <- hasil_validasi[order(hasil_validasi$RMSE, na.last = TRUE), ]

cat("========================================\n")

## ========================================

cat("HASIL VALIDASI CROSS-VALIDATION\n")

## HASIL VALIDASI CROSS-VALIDATION

cat("========================================\n\n")

## ========================================

print(hasil_validasi, row.names = FALSE)

##              Metode     RMSE      MAE           ME Varians_Kriging
##        TIN (Linear) 1.418993 1.176906  0.006022098              NA
##                 IDW 1.459661 1.180494  0.042583074              NA
##    Ordinary Kriging 1.475591 1.230538 -0.007289128        2.370894
##          Co-Kriging 1.491515 1.242707  0.015895032              NA
##   Universal Kriging 1.621664 1.353181  0.069037511        2.606101
##  Thin Plate Splines 1.706161 1.313684 -0.009514328              NA
##    Nearest Neighbor 1.756233 1.311379 -0.021724138              NA

# ---------------------------------------------------------
# 10. IDENTIFIKASI MODEL TERBAIK
# ---------------------------------------------------------
model_terbaik <- as.character(hasil_validasi$Metode[1])

cat("\n========================================\n")

## 
## ========================================

cat("MODEL TERBAIK:", model_terbaik, "\n")

## MODEL TERBAIK: TIN (Linear)

cat("========================================\n\n")

## ========================================

# Pilih prediksi berdasarkan model terbaik
pred_best <- NULL
obs_best <- NULL

if(grepl("Ordinary Kriging", model_terbaik, ignore.case = TRUE) && !is.null(cv_ok)) {
  pred_best <- cv_ok$var1.pred
  obs_best <- cv_ok$observed
  cat("✓ Menggunakan hasil Ordinary Kriging\n")
} else if(grepl("Universal Kriging", model_terbaik, ignore.case = TRUE) && !is.null(cv_uk)) {
  pred_best <- cv_uk$var1.pred
  obs_best <- cv_uk$observed
  cat("✓ Menggunakan hasil Universal Kriging\n")
} else if(grepl("Co-Kriging", model_terbaik, ignore.case = TRUE) && !is.null(pred_ck)) {
  pred_best <- pred_ck
  obs_best <- data_sp_ck_utm$TPT
  cat("✓ Menggunakan hasil Co-Kriging\n")
} else if(grepl("Nearest Neighbor", model_terbaik, ignore.case = TRUE)) {
  pred_best <- pred_nn
  obs_best <- data_clean$TPT
  cat("✓ Menggunakan hasil Nearest Neighbor\n")
} else if(grepl("TIN", model_terbaik, ignore.case = TRUE)) {
  pred_best <- pred_tin
  obs_best <- data_clean$TPT
  cat("✓ Menggunakan hasil TIN (Linear)\n")
} else if(grepl("IDW", model_terbaik, ignore.case = TRUE)) {
  pred_best <- pred_idw
  obs_best <- data_clean$TPT
  cat("✓ Menggunakan hasil IDW\n")
} else if(grepl("Thin Plate", model_terbaik, ignore.case = TRUE)) {
  pred_best <- pred_tps
  obs_best <- data_clean$TPT
  cat("✓ Menggunakan hasil Thin Plate Splines\n")
}

## ✓ Menggunakan hasil TIN (Linear)

# ---------------------------------------------------------
# 11. VISUALISASI (JIKA DATA VALID)
# ---------------------------------------------------------
if(!is.null(pred_best) && !is.null(obs_best)) {
  
  # Buat dataframe
  scatter_df <- data.frame(
    Observed = obs_best,
    Predicted = pred_best
  )
  
  # Hapus NA
  scatter_df <- scatter_df[complete.cases(scatter_df), ]
  
  if(nrow(scatter_df) > 0) {
    # Hitung R-squared
    lm_fit <- lm(Predicted ~ Observed, data = scatter_df)
    r_squared <- summary(lm_fit)$r.squared
    
    # Hitung residual
    scatter_df$Residual <- scatter_df$Observed - scatter_df$Predicted
    
    # PLOT 1: Scatter Plot
    p1 <- ggplot(scatter_df, aes(x = Observed, y = Predicted)) +
      geom_point(alpha = 0.6, color = "#3c8dbc", size = 3) +
      geom_abline(slope = 1, intercept = 0, linetype = "dashed", 
                  color = "red", size = 1) +
      geom_smooth(method = "lm", se = TRUE, color = "darkgreen", 
                  fill = "lightgreen", alpha = 0.2) +
      annotate("text", 
               x = min(scatter_df$Observed, na.rm = TRUE) + 0.5,
               y = max(scatter_df$Predicted, na.rm = TRUE) - 0.5,
               label = paste0("R² = ", round(r_squared, 3)),
               size = 5, fontface = "bold") +
      labs(
        title = paste("Observed vs Predicted TPT -", model_terbaik),
        subtitle = paste("RMSE:", round(hasil_validasi$RMSE[1], 3), 
                         "| MAE:", round(hasil_validasi$MAE[1], 3),
                         "| ME:", round(hasil_validasi$ME[1], 3)),
        x = "TPT Observed (%)",
        y = "TPT Predicted (%)"
      ) +
      theme_minimal(base_size = 12) +
      theme(
        plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        plot.subtitle = element_text(size = 11, hjust = 0.5),
        panel.grid.major = element_line(color = "gray90"),
        panel.grid.minor = element_line(color = "gray95")
      )
    
    print(p1)
    
    # PLOT 2: Residual Plot
    p2 <- ggplot(scatter_df, aes(x = Predicted, y = Residual)) +
      geom_point(alpha = 0.6, color = "purple", size = 3) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
      geom_smooth(method = "loess", se = TRUE, color = "blue", 
                  fill = "lightblue", alpha = 0.2) +
      labs(
        title = paste("Residual Plot -", model_terbaik),
        subtitle = "Residual seharusnya tersebar acak di sekitar 0",
        x = "TPT Predicted (%)",
        y = "Residual (Observed - Predicted)"
      ) +
      theme_minimal(base_size = 12) +
      theme(
        plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        plot.subtitle = element_text(size = 11, hjust = 0.5)
      )
    
    print(p2)
    
    # PLOT 3: Histogram Residual
    p3 <- ggplot(scatter_df, aes(x = Residual)) +
      geom_histogram(aes(y = ..density..), bins = 15, 
                     fill = "steelblue", color = "black", alpha = 0.7) +
      geom_density(color = "red", size = 1) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "darkred", size = 1) +
      labs(
        title = paste("Distribusi Residual -", model_terbaik),
        subtitle = paste("Mean:", round(mean(scatter_df$Residual), 3),
                         "| SD:", round(sd(scatter_df$Residual), 3)),
        x = "Residual",
        y = "Density"
      ) +
      theme_minimal(base_size = 12) +
      theme(
        plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        plot.subtitle = element_text(size = 11, hjust = 0.5)
      )
    
    print(p3)
    
    # ---------------------------------------------------------
    # 12. INTERPRETASI DETAIL
    # ---------------------------------------------------------
    cat("\n========================================\n")
    cat("INTERPRETASI HASIL VALIDASI\n")
    cat("========================================\n\n")
    
    cat("METRIK AKURASI:\n")
    cat("  RMSE         :", round(hasil_validasi$RMSE[1], 4), 
        "(semakin kecil semakin baik)\n")
    cat("  MAE          :", round(hasil_validasi$MAE[1], 4), 
        "(semakin kecil semakin baik)\n")
    cat("  ME (Bias)    :", round(hasil_validasi$ME[1], 4), 
        "(mendekati 0 = tidak bias)\n")
    cat("  R-squared    :", round(r_squared, 4), 
        "(mendekati 1 = model baik)\n\n")
    
    # Interpretasi RMSE
    cat("EVALUASI AKURASI:\n")
    if(hasil_validasi$RMSE[1] < 0.5) {
      cat("✓ RMSE sangat rendah - Model sangat akurat\n")
    } else if(hasil_validasi$RMSE[1] < 1.0) {
      cat("✓ RMSE rendah - Model cukup akurat\n")
    } else if(hasil_validasi$RMSE[1] < 1.5) {
      cat("⚠ RMSE sedang - Model acceptable\n")
    } else {
      cat("✗ RMSE tinggi - Model perlu perbaikan\n")
    }
    
    # Interpretasi Bias
    if(abs(hasil_validasi$ME[1]) < 0.1) {
      cat("✓ Bias sangat kecil - Prediksi seimbang\n")
    } else if(abs(hasil_validasi$ME[1]) < 0.3) {
      cat("✓ Bias kecil - Prediksi cukup balanced\n")
    } else {
      cat("⚠ Bias cukup besar - Ada kecenderungan sistematis\n")
    }
    
    # Interpretasi R-squared
    if(r_squared > 0.8) {
      cat("✓ R² tinggi - Model menjelaskan variasi sangat baik\n")
    } else if(r_squared > 0.6) {
      cat("✓ R² cukup - Model menjelaskan variasi dengan baik\n")
    } else if(r_squared > 0.3) {
      cat("⚠ R² sedang - Variasi data cukup terjelaskan\n")
    } else {
      cat("✗ R² rendah - Variasi data kurang terjelaskan\n")
    }
    
    cat("\n========================================\n")
    cat("KESIMPULAN:\n")
    cat("========================================\n")
    cat("Hasil validasi menunjukkan bahwa metode", model_terbaik, 
        "merupakan interpolasi terbaik dengan RMSE", round(hasil_validasi$RMSE[1], 4), 
        "dan MAE", round(hasil_validasi$MAE[1], 4), ".\n")
    cat("Rendahnya performa seluruh metode interpolasi menunjukkan bahwa pola TPT\n")
    cat("antarprovinsi di Indonesia sangat kompleks dan dipengaruhi oleh faktor\n")
    cat("non-geografis seperti struktur ekonomi lokal, kebijakan daerah, dan\n")
    cat("karakteristik industri yang tidak dapat sepenuhnya ditangkap hanya\n")
    cat("melalui posisi spasial.\n")
    cat("========================================\n\n")
    
  } else {
    cat("⚠ Warning: Data tidak cukup untuk membuat visualisasi\n")
  }
  
} else {
  cat("⚠ Warning: Tidak dapat membuat visualisasi untuk model terbaik\n")
}

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## 
## ========================================
## INTERPRETASI HASIL VALIDASI
## ========================================
## 
## METRIK AKURASI:
##   RMSE         : 1.419 (semakin kecil semakin baik)
##   MAE          : 1.1769 (semakin kecil semakin baik)
##   ME (Bias)    : 0.006 (mendekati 0 = tidak bias)
##   R-squared    : 0.0707 (mendekati 1 = model baik)
## 
## EVALUASI AKURASI:
## ⚠ RMSE sedang - Model acceptable
## ✓ Bias sangat kecil - Prediksi seimbang
## ✗ R² rendah - Variasi data kurang terjelaskan
## 
## ========================================
## KESIMPULAN:
## ========================================
## Hasil validasi menunjukkan bahwa metode TIN (Linear) merupakan interpolasi terbaik dengan RMSE 1.419 dan MAE 1.1769 .
## Rendahnya performa seluruh metode interpolasi menunjukkan bahwa pola TPT
## antarprovinsi di Indonesia sangat kompleks dan dipengaruhi oleh faktor
## non-geografis seperti struktur ekonomi lokal, kebijakan daerah, dan
## karakteristik industri yang tidak dapat sepenuhnya ditangkap hanya
## melalui posisi spasial.
## ========================================

Hasil validasi menunjukkan bahwa metode TIN (Linear) merupakan interpolasi terbaik dengan RMSE 1,4359 dan MAE 1,1737, meskipun akurasinya masih kurang optimal karena hanya mampu menjelaskan 4,33% variasi data (R² = 0,0433). Bias yang sangat kecil (ME = -0,1074) mengindikasikan prediksi cukup seimbang tanpa kecenderungan sistematis overestimate atau underestimate. Rendahnya performa seluruh metode interpolasi, baik deterministik (NN, TIN, IDW, TPS) maupun geostatistik (OK, UK, Co-Kriging)—menunjukkan bahwa pola TPT antarprovinsi di Indonesia sangat kompleks dan dipengaruhi oleh faktor non-geografis seperti struktur ekonomi lokal, kebijakan daerah, dan karakteristik industri yang tidak dapat sepenuhnya ditangkap hanya melalui posisi spasial.

4.18 Spatial Analysis Dashboard

Dashboard menyajikan analisis regresi spasial untuk membandingkan berbagai model ekonometrika dalam memahami dan memprediksi tingkat pengangguran terbuka (TPT) di seluruh provinsi Indonesia. Model yang digunakan meliputi OLS (Ordinary Least Squares) sebagai model regresi klasik, SAR (Spatial Autoregressive Model) yang menangkap efek limpahan antarwilayah, SEM (Spatial Error Model) yang memperhitungkan korelasi spasial pada galat, serta SDM (Spatial Durbin Model) yang menggabungkan komponen lag dan error.

Dashboard juga menyajikan analisisi GWR dan visualisasinya. Selain itu dashboard juga menampilkan visualisasi analisis Spasial interpolasi. Validasi pada analisis sapial interpolasi juga dilakukan untuk mengetahui mana metode terbaik dengan melihat nilai RMSE terkecil.

Variabel dependen yang dianalisis adalah TPT, dengan variabel independen meliputi Lama.Sekolah dan TPAK. Dashboard ini dilengkapi fitur interaktif seperti visualisasi spasial dengan label provinsi, perbandingan model berdasarkan metrik statistik, analisis tren deret waktu, analisis korelasi global dan lokal, serta opsi untuk mengunduh hasil dan visualisasi. Untuk menjalankan dashboard ini, pengguna perlu mengunggah file CSV yang berisi kolom Provinsi, Tahun, TPT, Lama.Sekolah, dan Pend_TPAK, serta file shapefile (.rds) untuk batas wilayah provinsi Indonesia.

Link Dashboard : https://hadeelina23001.shinyapps.io/spasial_140610230083UAS/

Link Video : Link Youtube : https://youtu.be/2FO4vdExSxM

BAB V (HASIL DAN PEMBAHASAN)

5.1 Kesimpulan

Berdasarkan hasil analisis spasial terhadap Tingkat Pengangguran Terbuka (TPT) di Indonesia periode 2020–2025, penelitian ini berhasil mengidentifikasi adanya pola keterkaitan spasial yang signifikan antarprovinsi. Hasil uji autokorelasi spasial menggunakan Moran’s I, Geary’s C, dan Getis-Ord G secara konsisten menunjukkan bahwa distribusi TPT tidak bersifat acak, melainkan membentuk pola pengelompokan (clustering) di mana provinsi dengan tingkat pengangguran tinggi cenderung berdekatan dengan provinsi lain yang memiliki kondisi serupa, dan sebaliknya. Temuan ini mengonfirmasi hipotesis bahwa tingkat pengangguran di suatu wilayah tidak hanya dipengaruhi oleh faktor internal wilayah tersebut, tetapi juga oleh kondisi dan karakteristik wilayah sekitarnya, yang dapat dijelaskan melalui mekanisme spillover effect seperti mobilitas tenaga kerja antarprovinsi, difusi kebijakan ketenagakerjaan regional, dan keterkaitan struktur ekonomi antarwilayah.

Perbandingan model ekonometrika spasial menunjukkan bahwa Spatial Error Model (SEM) merupakan model terbaik dalam menjelaskan variasi TPT di Indonesia, dengan nilai AIC terendah (461,74) dibandingkan model OLS (540,73), SAR (480,36), dan SDM (465,42). Koefisien lambda (λ) yang signifikan sebesar 0,721 pada model SEM mengindikasikan bahwa ketergantungan spasial lebih dominan berasal dari komponen galat (error term) dibandingkan dari pengaruh langsung variabel dependen antarwilayah. Hal ini mengimplikasikan bahwa faktor-faktor tidak teramati yang memengaruhi pengangguran seperti kebijakan daerah, iklim investasi, kualitas infrastruktur, dan kondisi kelembagaan memiliki korelasi spasial yang kuat antarprovinsi. Secara substantif, variabel Tingkat Partisipasi Angkatan Kerja (TPAK) terbukti berpengaruh negatif dan signifikan terhadap TPT, yang mengindikasikan bahwa peningkatan partisipasi aktif penduduk usia kerja dalam pasar tenaga kerja dapat mengurangi tingkat pengangguran. Sementara itu, variabel rata-rata lama sekolah menunjukkan pengaruh positif dan signifikan, yang mencerminkan fenomena skill mismatch di mana peningkatan pendidikan tidak selalu diikuti oleh penyerapan tenaga kerja yang memadai jika tidak sejalan dengan kebutuhan pasar kerja lokal.

Analisis Geographically Weighted Regression (GWR) berhasil mengungkap adanya heterogenitas spasial yang signifikan dalam hubungan antara variabel independen dan TPT. Koefisien lokal untuk variabel rata-rata lama sekolah berkisar antara -0,6141 hingga 0,1889, sedangkan koefisien TPAK berkisar antara -0,2928 hingga 0,2408, menunjukkan bahwa pengaruh kedua variabel tersebut tidak seragam antarprovinsi melainkan bervariasi secara geografis. Bandwidth adaptif optimal sebesar 0,050 menandakan bahwa model GWR menggunakan sekitar 5% tetangga terdekat dalam estimasi lokal, yang mengindikasikan variasi spasial yang cukup kuat dan perlunya pendekatan lokal dalam analisis. Model GWR juga menunjukkan performa superior dibandingkan model OLS dengan nilai AICc yang jauh lebih rendah (384,26 vs 540,73), mengonfirmasi bahwa pendekatan regresi global tidak memadai untuk menangkap kompleksitas hubungan spasial yang bervariasi antarwilayah. Uji asumsi menunjukkan bahwa model GWR tidak memiliki masalah multikolinearitas (VIF < 2) dan autokorelasi spasial residual (Moran’s I tidak signifikan), meskipun asumsi normalitas residual tidak terpenuhi, namun hal ini dapat ditoleransi dalam konteks GWR yang lebih menekankan pada estimasi parameter lokal dibandingkan inferensi statistik global.

Perbandingan berbagai metode interpolasi spasial menunjukkan bahwa metode geostatistik (Ordinary Kriging, Universal Kriging, dan Co-Kriging) memberikan hasil yang lebih superior dibandingkan metode deterministik (Nearest Neighbor, TIN, IDW, dan Thin Plate Splines). Hasil validasi menunjukkan bahwa metode TIN (Linear) merupakan interpolasi terbaik dengan RMSE 1,4359 dan MAE 1,1737, meskipun akurasinya masih kurang optimal karena hanya mampu menjelaskan 4,33% variasi data (R² = 0,0433). Bias yang sangat kecil (ME = -0,1074) mengindikasikan prediksi cukup seimbang tanpa kecenderungan sistematis overestimate atau underestimate. Rendahnya performa seluruh metode interpolasi, baik deterministik (NN, TIN, IDW, TPS) maupun geostatistik (OK, UK, Co-Kriging) menunjukkan bahwa pola TPT antarprovinsi di Indonesia sangat kompleks dan dipengaruhi oleh faktor non-geografis seperti struktur ekonomi lokal, kebijakan daerah, dan karakteristik industri yang tidak dapat sepenuhnya ditangkap hanya melalui posisi spasial.

Secara keseluruhan, penelitian ini mendemonstrasikan pentingnya pendekatan analisis spasial dalam memahami fenomena pengangguran di Indonesia, di mana ketergantungan spasial (spatial dependence) dan heterogenitas spasial (spatial heterogeneity) menjadi karakteristik inheren yang tidak dapat diabaikan. Integrasi antara model regresi spasial global (SEM), model regresi spasial lokal (GWR), dan metode interpolasi geostatistik (Kriging) memberikan pemahaman yang komprehensif tentang dinamika TPT dari berbagai dimensi, keterkaitan antarwilayah melalui komponen error (SEM), variasi lokal pengaruh variabel penjelas (GWR), dan prediksi nilai di lokasi tidak tersampel dengan ukuran ketidakpastian (Kriging). Temuan ini memiliki implikasi penting bagi perumusan kebijakan ketenagakerjaan di Indonesia, yang seharusnya tidak bersifat uniform untuk seluruh wilayah tetapi disesuaikan dengan konteks lokal dan mempertimbangkan efek spillover antarprovinsi. Kebijakan yang efektif perlu memperhatikan tidak hanya karakteristik internal suatu provinsi, tetapi juga kondisi provinsi-provinsi di sekitarnya, serta mengakui bahwa hubungan antara pendidikan, partisipasi angkatan kerja, dan pengangguran bervariasi secara geografis sehingga memerlukan strategi intervensi yang berbeda untuk wilayah yang berbeda. Penelitian lebih lanjut disarankan untuk menambahkan variabel penjelas lain seperti pertumbuhan ekonomi, upah minimum, investasi daerah, dan sektor dominan ekonomi, serta menggunakan data dengan resolusi temporal dan spasial yang lebih tinggi untuk meningkatkan akurasi dan relevansi analisis bagi pengambilan kebijakan regional.

5.2 Saran

Berdasarkan hasil penelitian ini, terdapat beberapa saran yang dapat diberikan:

1. Pemerintah daerah perlu memperhatikan keterkaitan spasial antarprovinsi dalam menyusun kebijakan ketenagakerjaan. Upaya penurunan pengangguran sebaiknya dilakukan secara terpadu antarwilayah, terutama pada provinsi yang membentuk klaster pengangguran tinggi.
2. Peningkatan partisipasi angkatan kerja perlu difokuskan melalui program pelatihan keterampilan, pendidikan vokasi, dan penciptaan lapangan kerja produktif agar dapat menekan tingkat pengangguran secara berkelanjutan.
3. Penelitian selanjutnya disarankan untuk menambahkan variabel lain seperti pertumbuhan ekonomi, pendidikan, dan investasi agar model dapat menangkap lebih banyak faktor penyebab pengangguran.
4. Penggunaan data dengan resolusi spasial dan temporal yang lebih tinggi juga direkomendasikan agar hasil analisis lebih akurat dan relevan untuk pengambilan kebijakan di tingkat regional.
5. Penggunaan data dengan resolusi spasial lebih tinggi (kabupaten/kota) dan temporal lebih detail (triwulanan/bulanan) direkomendasikan untuk meningkatkan akurasi metode interpolasi spasial yang saat ini menunjukkan R² hanya 4,33% dan memberikan hasil yang lebih relevan untuk kebijakan tingkat regional.

REFERENSI

1. Arbia, G. (2014). A primer for spatial econometrics: With applications in R. Palgrave Macmillan.

2. Anselin, L. (1988). Spatial econometrics: Methods and models. Kluwer Academic Publishers.

3. Elhorst, J. P. (2014). Spatial econometrics: From cross-sectional data to spatial panels. SpringerBriefs in Regional Science.

4. Lee, L. F. (2021). Spatial econometrics: Spatial autoregressive models. World Scientific.

5. Oktafianto, E. K., et al. (2019). The determinant of regional unemployment in Indonesia: The spatial Durbin models. Signifikan: Jurnal Ilmu Ekonomi

6. Yanuar, F. (2023). Spatial autoregressive quantile regression with open unemployment level in all provinces in Indonesia. Jurnal Sains Teknologi Indonesia

7. Fauzi, F., Wenur, G. H., & Wasono, R. (2023). Spatial Durbin Model of Unemployment Rate in Central Java. Parameter: Journal of Statistics

8. Rosmanah, R., Damu Djara, V. A., & Jaya, I. G. N. M. (2022). The spatial econometrics of economic growth in Sumatera Utara Province. J. Math. Comput. Sci.

9. Hariyatih, S. (2024). Analisis Pengaruh IPM, UMP dan Tingkat Partisipasi Angkatan Kerja terhadap Pengangguran Terbuka di Indonesia. Journal of Social and Economics Research

10. Adianita, H., Susilowati, D., & Putri Karisma, D. A. (2024). Analisis Tingkat Partisipasi Angkatan Kerja dan Jumlah Lapangan Kerja melalui Pendidikan terhadap Tingkat Pengangguran di Indonesia. Development Review

11. Ade Kantari, D. P., Kamila, N., & Safitri, O. (2024). Analisis Pengaruh Indeks Pendidikan, Tingkat Partisipasi Angkatan Kerja dan Jumlah TPAK terhadap Tingkat Pengangguran Terbuka pada Kabupaten/Kota Provinsi Jawa Timur Tahun 2021. Acmatics Journal: Actuarial, Mathematics, and Statistics Journal

12. Safrina, Y., & Ratna, R. (2023). Analisis Pengaruh Tingkat Partisipasi Angkatan Kerja, Indeks Pembangunan Manusia, dan Tingkat Pengangguran Terbuka terhadap Pertumbuhan Ekonomi di Indonesia. Jurnal Ekonomi Regional Unimal

13. Badan Pusat Statistik. Tingkat Pengangguran Terbuka menurut Provinsi. Tabel Statistik.[https://www.bps.go.id/id/statistics-table/2/NTQzIzI=/tingkat-pengangguran-terbuka-menurut-provinsi.html]

14. Badan Pusat Statistik. Persentase TPAK TPAK (P0) menurut Provinsi dan Daerah - persen. Tabel Statistik.[https://www.bps.go.id/id/statistics-table/2/MTkyIzI=/persentase-TPAK-TPAK--p0--menurut-provinsi-dan-daerah--persen-.html]

15. Badan Pusat Statistik Persentase Angkatan Kerja terhadap TPAK Usia Kerja (Lama.Sekolah) menurut Provinsi - persen. Tabel Statistik.[https://www.bps.go.id/id/statistics-table/2/MjM5NiMy/persentase-angkatan-kerja-terhadap-TPAK-usia-kerja--Lama.Sekolah--menurut-provinsi--persen-.html]

16. United Nations. (2015). Transforming our world: The 2030 Agenda for Sustainable Development. United Nations Department of Economic and Social Affairs. https://sdgs.un.org/2030agenda

17. Anselin, L. 1995. Local Indicator of Spatial Association-LISA. Geographical Analysis, 27(2): 93-115

18. Saputro, D. R., Widyaningsih, P., Kurdi, N. A., & Susanti, A. 2018. Proporsionalitas Autokorelasi Spasial dengan Indeks Global (Indeks Moran) dan Indeks Lokal (Local Indicator of Spatial Association (LISA)). KNPMP III. Hal : 701-709. ISSN : 2502-6526

19. Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. John Wiley & Sons, Inc

20. Caraka, R. E., & Yasin, H. (2017). Geographically Weighted Regression (GWR) ; Sebuah Pendekatan Regresi Geografis. Mobius

21. Carisa Putri Salsabila Purnamasari, & Yekti Widyaningsih. (2022). Perbandingan Performa Bandwidth CV, AICc, dan BIC pada Model Geographically Weighted Regression

22. Hadi, B. S. (2013). Metode interpolasi spasial dalam studi geografi.

23. Yudanegara, R. A., Astutik, D., Hernandi, A., Soedarmodjo, T. P., & Alexander, E. (2021). PENGGUNAAN METODE INVERSE DISTANCE WEIGHTED (IDW) UNTUK PEMETAAN ZONA NILAI TANAH (STUDI KASUS: KELURAHAN GEDONG MENENG, BANDAR LAMPUNG)

24. Christou, N. (2021). Spatial Data in Undergraduate Statistics Curriculum. JOURNAL OF STATISTICS AND DATA SCIENCE EDUCATION

25. Cressie, N. (1993). Statistics for Spatial Data. New York: Wiley.

26. Chilès, J. P., & Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty. Hoboken, NJ: Wiley.

27. Isaaks, E. H., & Srivastava, R. M. (1989). An Introduction to Applied Geostatistics. New York: Oxford University Press.

28. Wackernagel, H. (2003). Multivariate Geostatistics. Berlin: Springer.