PROJECT SPASIAL
LOAD & PERSIAPAN DATA
# STEP 1: LOAD DATA
sf_use_s2(FALSE)
jatim_shp <- st_read(
"D:/KULIAH/SEMESTER 6/Statistika Spasial/jatim"
)## Reading layer `jatim' from data source
## `D:\KULIAH\SEMESTER 6\Statistika Spasial\jatim' using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 110.8953 ymin: -8.780223 xmax: 116.2702 ymax: -5.042965
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defsdata_jatim <- read.csv(
"D:/KULIAH/SEMESTER 6/Statistika Spasial/datajatim.csv",
stringsAsFactors = FALSE
)
colnames(data_jatim) <- c("Kab_Kota", "Kemiskinan", "RLS", "Jumlah_Sekolah",
"PDRB", "APBD", "TPT", "AHH", "Nakes", "Faskes", "Kepadatan")Interpretasi: Data spasial Jawa Timur berhasil dimuat dengan 38 wilayah (kabupaten/kota) dan 8 field atribut. CRS yang digunakan adalah WGS84 dengan tipe geometri MULTIPOLYGON. Data CSV memiliki 11 variabel yang mencakup indikator sosial-ekonomi dan pembangunan.
CEK KECOCOKAN NAMA SHP vs CSV
# STEP 2: CEK KECOCOKAN NAMA SHP vs CSV
nama_shp <- sort(jatim_shp$WADMKK)
nama_csv <- sort(data_jatim$Kab_Kota)
tidak_cocok_shp <- setdiff(nama_shp, nama_csv)
tidak_cocok_csv <- setdiff(nama_csv, nama_shp)
cat("\n NAMA DI SHP TAPI TIDAK ADA DI CSV \n")##
## NAMA DI SHP TAPI TIDAK ADA DI CSV## character(0)##
## NAMA DI CSV TAPI TIDAK ADA DI SHP## character(0)Interpretasi: Hasil pengecekan menunjukkan
character(0)pada kedua sisi, artinya tidak ada ketidakcocokan nama antara shapefile dan CSV. Seluruh 38 nama kabupaten/kota pada shapefile berhasil dicocokkan sempurna dengan data CSV, sehingga proses merge dapat dilanjutkan tanpa perlu perbaikan nama.
GABUNGKAN SHP + CSV
# STEP 3: GABUNGKAN SHP + CSV
jatim_merged <- jatim_shp %>%
left_join(data_jatim, by = c("WADMKK" = "Kab_Kota")) %>%
filter(!is.na(Kemiskinan))
cat("\n CEK NA SETELAH MERGE \n")##
## CEK NA SETELAH MERGE## KDPKAB KDPPUM WADMKK WADMPR METADATA
## 0 0 0 0 0
## UPDATED sumut jatim Kemiskinan RLS
## 0 0 0 0 0
## Jumlah_Sekolah PDRB APBD TPT AHH
## 0 0 0 0 0
## Nakes Faskes Kepadatan
## 0 0 0##
## JUMLAH WILAYAH SETELAH MERGE: 38Interpretasi: Proses merge berhasil dengan sempurna. Seluruh kolom menunjukkan nilai NA = 0, artinya tidak ada data yang hilang pada semua variabel (Kemiskinan, RLS, Jumlah_Sekolah, PDRB, APBD, TPT, AHH, Nakes, Faskes, Kepadatan). Jumlah wilayah setelah merge tetap 38 wilayah, sesuai dengan jumlah kabupaten/kota di Jawa Timur.
MATRIKS PEMBOBOT SPASIAL (Queen Contiguity)
# STEP 4: BUAT MATRIKS PEMBOBOT SPASIAL (Queen Contiguity)
nb <- poly2nb(jatim_merged, queen = TRUE)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)
cat("\n RINGKASAN BOBOT SPASIAL \n")##
## RINGKASAN BOBOT SPASIAL## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 138
## Percentage nonzero weights: 9.556787
## Average number of links: 3.631579
## 2 disjoint connected subgraphs
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9
## 8 6 6 6 2 7 1 1 1
## 8 least connected regions:
## 26 29 30 31 32 33 34 35 with 1 link
## 1 most connected region:
## 7 with 9 linksInterpretasi: Matriks pembobot spasial berhasil dibentuk menggunakan metode Queen Contiguity (berbagi sisi maupun titik sudut). Ringkasan menunjukkan:
- Jumlah wilayah: 38
- Total tautan (links): 138, dengan persentase nonzero weights sebesar 9,56%
- Rata-rata jumlah tetangga: 3,63 per wilayah
- Terdapat 2 subgraf disjoint, mengindikasikan adanya wilayah kepulauan (Madura) yang tidak berbagi batas darat dengan wilayah lain
- Wilayah paling terisolasi: 8 wilayah (indeks 26, 29–35) hanya memiliki 1 tetangga — kemungkinan besar merupakan wilayah kepulauan seperti Sumenep
- Wilayah paling terhubung: wilayah indeks 7 memiliki 9 tetangga, menunjukkan posisi sentral secara geografis
REGRESI OLS
# STEP 5: REGRESI OLS
# Variabel dependen : Kemiskinan
# Variabel independen: RLS, Jumlah_Sekolah, PDRB, APBD, TPT, AHH, Nakes, Faskes
ols_model <- lm(
Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD + TPT + AHH + Kepadatan,
data = jatim_merged
)
summary(ols_model)##
## Call:
## lm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD +
## TPT + AHH + Kepadatan, data = jatim_merged)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5885 -1.5133 0.2031 1.4953 5.1110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.0534524 69.2656755 0.925 0.3625
## RLS -1.8510476 0.8587549 -2.156 0.0393 *
## Jumlah_Sekolah 0.0015425 0.0013634 1.131 0.2668
## PDRB 0.0022022 0.0056443 0.390 0.6992
## APBD -0.0003098 0.0003251 -0.953 0.3483
## TPT -0.4703897 0.5229267 -0.900 0.3755
## AHH -0.4995198 0.9781347 -0.511 0.6133
## Kepadatan 0.0002243 0.0003560 0.630 0.5334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.713 on 30 degrees of freedom
## Multiple R-squared: 0.6626, Adjusted R-squared: 0.5839
## F-statistic: 8.417 on 7 and 30 DF, p-value: 1.115e-05##
## Durbin-Watson test
##
## data: ols_model
## DW = 0.7993, p-value = 3.899e-06
## alternative hypothesis: true autocorrelation is greater than 0Interpretasi: Model OLS menghasilkan R² = 0,6626 (Adjusted R² = 0,5839), artinya sekitar 66,26% variasi tingkat kemiskinan dapat dijelaskan oleh tujuh variabel prediktor. Model secara keseluruhan signifikan (F = 8,417; p < 0,001).
Secara parsial, hanya variabel RLS (Rata-rata Lama Sekolah) yang signifikan pada α = 5% (t = -2,156; p = 0,039), dengan koefisien -1,851 — setiap kenaikan 1 tahun RLS dikaitkan dengan penurunan tingkat kemiskinan sebesar 1,85 persen, ceteris paribus.
Variabel lainnya (Jumlah_Sekolah, PDRB, APBD, TPT, AHH, Kepadatan) tidak signifikan secara parsial pada α = 5%.
Uji Durbin-Watson menghasilkan DW = 0,7993 (p = 3,899e-06), jauh di bawah nilai 2, sehingga H₀ ditolak — terdapat autokorelasi positif yang kuat pada residual OLS. Ini menjadi dasar kuat untuk melanjutkan ke uji spasial.
UJI MORAN’S I PADA RESIDUAL OLS
# STEP 6: UJI MORAN'S I PADA RESIDUAL OLS
residual_ols <- residuals(ols_model)
moran_test <- moran.test(
residual_ols,
listw = listw,
zero.policy = TRUE
)
print(moran_test)##
## Moran I test under randomisation
##
## data: residual_ols
## weights: listw
##
## Moran I statistic standard deviate = 3.3418, p-value = 0.0004162
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.40874658 -0.02702703 0.01700458Interpretasi: Uji Moran’s I pada residual OLS menghasilkan Moran I = 0,4087 dengan standar deviat = 3,3418 dan p-value = 0,0004162. Karena p-value jauh di bawah α = 5%, maka H₀ ditolak — terdapat autokorelasi spasial yang signifikan pada residual model OLS.
Nilai Moran I positif (0,41) mengindikasikan bahwa wilayah-wilayah dengan tingkat kemiskinan serupa cenderung berkelompok secara spasial (spatial clustering). Kondisi ini melanggar asumsi independensi residual OLS, sehingga model OLS tidak memadai dan pemodelan spasial diperlukan.
UJI LAGRANGE MULTIPLIER (LM)
# STEP 7: UJI LAGRANGE MULTIPLIER (LM)
lm_tests <- lm.LMtests(
ols_model,
listw = listw,
test = c("LMerr", "LMlag", "RLMerr", "RLMlag"),
zero.policy = TRUE
)
print(lm_tests)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD +
## TPT + AHH + Kepadatan, data = jatim_merged)
## test weights: listw
##
## RSerr = 8.9352, df = 1, p-value = 0.002797
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD +
## TPT + AHH + Kepadatan, data = jatim_merged)
## test weights: listw
##
## RSlag = 2.8659, df = 1, p-value = 0.09047
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD +
## TPT + AHH + Kepadatan, data = jatim_merged)
## test weights: listw
##
## adjRSerr = 8.1623, df = 1, p-value = 0.004277
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD +
## TPT + AHH + Kepadatan, data = jatim_merged)
## test weights: listw
##
## adjRSlag = 2.093, df = 1, p-value = 0.148Interpretasi: Hasil uji Lagrange Multiplier dirangkum sebagai berikut:
Uji Statistik p-value Keputusan LM Error (LMerr) 8,9352 0,0028 Signifikan ✓ LM Lag (LMlag) 2,8659 0,0905 Tidak Signifikan ✗ Robust LM Error (RLMerr) 8,1623 0,0043 Signifikan ✓ Robust LM Lag (RLMlag) 2,0930 0,1480 Tidak Signifikan ✗ Berdasarkan prosedur pemilihan model Anselin (1988): LM Error signifikan sedangkan LM Lag tidak, dan Robust LM Error juga signifikan sedangkan Robust LM Lag tidak. Ini secara konsisten mengarah pada Spatial Error Model (SEM) sebagai model yang tepat. Dengan menambahkan efek spasial pada lag variabel independen, pilihan jatuh pada SDEM (Spatial Durbin Error Model).
MODEL SDEM (Spatial Durbin Error Model)
# STEP 8: MODEL SDEM (Spatial Durbin Error Model)
sdem_model <- errorsarlm(
Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB + APBD + TPT + AHH + Kepadatan,
data = jatim_merged,
listw = listw,
etype = "emixed",
zero.policy = TRUE
)
summary(sdem_model)##
## Call:errorsarlm(formula = Kemiskinan ~ RLS + Jumlah_Sekolah + PDRB +
## APBD + TPT + AHH + Kepadatan, data = jatim_merged, listw = listw,
## etype = "emixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4947274 -0.9202136 -0.0060823 0.9821162 3.0313167
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.5941e+01 1.1479e+02 0.2260 0.821212
## RLS -2.6559e+00 6.2940e-01 -4.2197 2.446e-05
## Jumlah_Sekolah 1.7990e-03 1.0025e-03 1.7945 0.072731
## PDRB 7.0504e-03 4.8979e-03 1.4395 0.150011
## APBD -5.2858e-04 2.4610e-04 -2.1478 0.031726
## TPT -3.7935e-01 4.1747e-01 -0.9087 0.363508
## AHH -1.9616e+00 8.8650e-01 -2.2128 0.026911
## Kepadatan 8.4266e-04 3.0846e-04 2.7318 0.006299
## lag.RLS -1.0612e+00 9.9421e-01 -1.0674 0.285789
## lag.Jumlah_Sekolah -4.1386e-03 1.7215e-03 -2.4040 0.016215
## lag.PDRB 1.8447e-02 1.6323e-02 1.1301 0.258426
## lag.APBD 1.1443e-03 5.9135e-04 1.9350 0.052988
## lag.TPT -8.0544e-01 8.3370e-01 -0.9661 0.333993
## lag.AHH 2.2078e+00 1.3351e+00 1.6536 0.098205
## lag.Kepadatan 2.4069e-04 1.0037e-03 0.2398 0.810478
##
## Lambda: 0.66163, LR test value: 12.022, p-value: 0.00052585
## Asymptotic standard error: 0.11316
## z-value: 5.8467, p-value: 5.0156e-09
## Wald statistic: 34.183, p-value: 5.0156e-09
##
## Log likelihood: -75.20881 for error model
## ML residual variance (sigma squared): 2.6237, (sigma: 1.6198)
## Number of observations: 38
## Number of parameters estimated: 17
## AIC: 184.42, (AIC for lm: 194.44)Interpretasi: Model SDEM berhasil diestimasi. Berikut temuan utama:
Variabel yang signifikan (efek lokal/direct):
- RLS (z = -4,220; p < 0,001): Setiap kenaikan 1 tahun rata-rata lama sekolah menurunkan kemiskinan sebesar 2,656%. Efek ini lebih kuat dibanding OLS, menunjukkan pentingnya pendidikan dalam mengurangi kemiskinan.
- APBD (z = -2,148; p = 0,032): Peningkatan belanja daerah berpengaruh negatif signifikan terhadap kemiskinan, meski koefisiennya sangat kecil (-0,00053).
- AHH (z = -2,213; p = 0,027): Angka Harapan Hidup yang lebih tinggi dikaitkan dengan kemiskinan yang lebih rendah (-1,962).
- Kepadatan (z = 2,732; p = 0,006): Kepadatan penduduk justru berpengaruh positif terhadap kemiskinan (0,00084), kemungkinan mencerminkan tekanan urban di wilayah padat.
Variabel lag spasial yang signifikan (efek spillover):
- lag.Jumlah_Sekolah (z = -2,404; p = 0,016): Jumlah sekolah SMA di wilayah tetangga berpengaruh negatif signifikan terhadap kemiskinan suatu wilayah, mengindikasikan adanya efek limpahan (spillover) dari ketersediaan fasilitas pendidikan antarwilayah.
Parameter Lambda (λ): λ = 0,6616 dengan p-value = 5,02e-09 — sangat signifikan, mengkonfirmasi bahwa terdapat dependensi spasial pada error yang kuat antar wilayah di Jawa Timur.
PERBANDINGAN R² DAN AIC
# STEP 9: PERBANDINGAN R² DAN AIC
y_actual <- jatim_merged$Kemiskinan
y_fitted <- fitted(sdem_model)
ss_res <- sum((y_actual - y_fitted)^2)
ss_tot <- sum((y_actual - mean(y_actual))^2)
r2_sdem <- 1 - (ss_res / ss_tot)
cat("R² OLS :", round(summary(ols_model)$r.squared, 4), "\n")## R² OLS : 0.6626## R² SDEM : 0.8477## AIC OLS : 192.72## AIC SDEM : 184.42Interpretasi: Perbandingan performa model OLS vs SDEM menunjukkan hasil yang jelas:
Metrik OLS SDEM Selisih R² 0,6626 0,8477 +0,1851 AIC 192,72 184,42 -8,30 Model SDEM unggul pada kedua kriteria: R² meningkat signifikan dari 0,6626 menjadi 0,8477 (+18,51 poin persentase), artinya SDEM mampu menjelaskan 84,77% variasi kemiskinan di Jawa Timur. Nilai AIC SDEM (184,42) juga lebih kecil dari OLS (192,72), mengindikasikan model yang lebih parsimoni dan lebih baik fit-nya.
Kesimpulan: Model SDEM terbukti lebih unggul dibandingkan OLS dalam memodelkan kemiskinan di Jawa Timur. Peningkatan performa ini dikaitkan dengan kemampuan SDEM mengakomodasi dependensi spasial pada error (λ = 0,66) dan efek spillover dari variabel di wilayah tetangga, yang diabaikan oleh OLS biasa.