# Load Package
suppressPackageStartupMessages({
library(sf)
library(readxl)
library(tidyverse)
library(ggplot2)
library(spdep)
library(spatialreg)
library(MuMIn)
library(patchwork)
})Spatial Regression
A. Model Prediksi GRP Data China
Deskripsi Data china_join
Data yang digunakan dalam analisis ini bersumber dari file Data_China.xlsx yang memuat data utama 29 provinsi di China, mencakup koordinat geografis (longitude dan latitude) masing-masing provinsi, Gross Regional Product (GRP) sebagai ukuran output ekonomi wilayah, serta tingkat urbanisasi dalam persen yang menggambarkan proporsi penduduk yang tinggal di kawasan perkotaan. Selain itu, matriks jarak kereta api antar provinsi (km) yang menghubungkan 29 ibu kota provinsi satu sama lain, dan digunakan sebagai dasar pembentukan matriks bobot spasial dalam pemodelan regresi spasial. Secara keseluruhan, data ini digunakan untuk menganalisis sejauh mana tingkat urbanisasi berpengaruh terhadap output ekonomi regional di China dengan mempertimbangkan struktur keterkaitan spasial antar provinsi melalui konektivitas jaringan kereta api.
1. Persiapan Data
# Membaca data spasial dalam format shapefile (.shp)
shp_china <- st_read("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/CHINA/China29/China29.shp")Reading layer `China29' from data source
`/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/CHINA/China29/China29.shp'
using driver `ESRI Shapefile'
Simple feature collection with 29 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 73.5577 ymin: 20.22264 xmax: 134.7739 ymax: 53.56086
Geodetic CRS: WGS 84data_china <- read_excel("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/CHINA/Data_China.xlsx")
china_join <- shp_china %>%
left_join(data_china,
by = c("NAME_1" = "Region"))
head(china_join)Simple feature collection with 6 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.76416 ymin: 20.22264 xmax: 120.8376 ymax: 42.77393
Geodetic CRS: WGS 84
ID_0 ISO NAME_0 ID_1 NAME_1 HASC_1 CCN_1 CCA_1 TYPE_1 ENGTYPE_1
1 49 CHN China 1 Anhui CN.AH 0 <NA> Sh?ng Province
2 49 CHN China 2 Beijing CN.BJ 0 <NA> Zhíxiáshì Municipality
3 49 CHN China 3 Chongqing CN.CQ 0 <NA> Zhíxiáshì Municipality
4 49 CHN China 4 Fujian CN.FJ 0 <NA> Sh?ng Province
5 49 CHN China 5 Gansu CN.GS 0 <NA> Sh?ng Province
6 49 CHN China 6 Guangdong CN.GD 0 <NA> Sh?ng Province
NL_NAME_1 VARNAME_1 Long Lat GRP Level_of_Urbanization
1 ??|?? ?nhu? 117.2262 31.82579 28792 46.50
2 ??|?? B?ij?ng 116.4107 40.18491 87475 86.20
3 ??|?? Chóngqìng 107.8748 30.05865 38914 56.98
4 ?? Fújiàn 117.9815 26.08641 52763 59.60
5 ??|?? G?nsù 100.9339 37.82079 21978 38.75
6 ??|?? Gu?ngd?ng 113.4164 23.34225 54095 67.40
geometry
1 MULTIPOLYGON (((116.4263 34...
2 MULTIPOLYGON (((116.6669 40...
3 MULTIPOLYGON (((108.5419 32...
4 MULTIPOLYGON (((117.689 23....
5 MULTIPOLYGON (((97.18472 42...
6 MULTIPOLYGON (((110.1182 20...2. Eksplorasi Data
GRP <- ggplot(data = china_join) +
geom_sf(
aes(fill = GRP),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(
option = "magma",
name = "GRP"
) +
labs(
title = "Gross Regional Product",
subtitle = "Data kontinu per provinsi",
x = "Longitude",
y = "Latitude"
) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11)
)
Urbanization <- ggplot(data = china_join) +
geom_sf(
aes(fill = `Level_of_Urbanization`),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(
option = "magma",
name = "Level of Urbanization"
) +
labs(
title = "Level of Urbanization",
subtitle = "Data kontinu per provinsi",
x = "Longitude",
y = "Latitude"
) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11)
)
GRP | Urbanization
# Pengecekan struktur data
str(china_join)Classes 'sf' and 'data.frame': 29 obs. of 17 variables:
$ ID_0 : chr "49" "49" "49" "49" ...
$ ISO : chr "CHN" "CHN" "CHN" "CHN" ...
$ NAME_0 : chr "China" "China" "China" "China" ...
$ ID_1 : chr "1" "2" "3" "4" ...
$ NAME_1 : chr "Anhui" "Beijing" "Chongqing" "Fujian" ...
$ HASC_1 : chr "CN.AH" "CN.BJ" "CN.CQ" "CN.FJ" ...
$ CCN_1 : chr "0" "0" "0" "0" ...
$ CCA_1 : chr NA NA NA NA ...
$ TYPE_1 : chr "Sh?ng" "Zhíxiáshì" "Zhíxiáshì" "Sh?ng" ...
$ ENGTYPE_1 : chr "Province" "Municipality" "Municipality" "Province" ...
$ NL_NAME_1 : chr "??|??" "??|??" "??|??" "??" ...
$ VARNAME_1 : chr "?nhu?" "B?ij?ng" "Chóngqìng" "Fújiàn" ...
$ Long : num 117 116 108 118 101 ...
$ Lat : num 31.8 40.2 30.1 26.1 37.8 ...
$ GRP : num 28792 87475 38914 52763 21978 ...
$ Level_of_Urbanization: num 46.5 86.2 57 59.6 38.8 ...
$ geometry :sfc_MULTIPOLYGON of length 29; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:7291, 1:2] 116 116 116 116 116 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
..- attr(*, "names")= chr [1:16] "ID_0" "ISO" "NAME_0" "ID_1" ...# Pengecekan Missing Value
colSums(is.na(china_join)) ID_0 ISO NAME_0
0 0 0
ID_1 NAME_1 HASC_1
0 0 0
CCN_1 CCA_1 TYPE_1
0 29 0
ENGTYPE_1 NL_NAME_1 VARNAME_1
0 0 0
Long Lat GRP
0 0 0
Level_of_Urbanization geometry
0 0 3. Regresi OLS
china.lm<- lm(GRP ~ Level_of_Urbanization,data=china_join)
summary(china.lm)
Call:
lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
Residuals:
Min 1Q Median 3Q Max
-12050.9 -4410.4 -374.8 2961.6 14948.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -31895.5 5197.6 -6.137 1.48e-06 ***
Level_of_Urbanization 1400.0 92.6 15.118 1.06e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6595 on 27 degrees of freedom
Multiple R-squared: 0.8944, Adjusted R-squared: 0.8904
F-statistic: 228.6 on 1 and 27 DF, p-value: 1.064e-144. Matriks Pembobot Spasial
Membangun matriks bobot spasial berbasis railway distance antar provinsi di China yang akan digunakan dalam analisis regresi spasial. Data railway distance dibaca dari file Excel, kemudian dikonversi menjadi matriks numerik berukuran 29×29 yang merepresentasikan jarak antar setiap pasang provinsi. Karena dalam pemodelan spasial kedekatan diukur sebagai kebalikan dari jarak, semakin dekat suatu wilayah, semakin besar pengaruhnya maka setiap elemen matriks ditransformasi menggunakan invers jarak (1/d), dengan nilai diagonal diset nol agar suatu wilayah tidak dianggap bertetangga dengan dirinya sendiri. Matriks kemudian dinormalisasi dengan row-standardization sehingga total bobot untuk setiap provinsi bernilai satu, memastikan skala pengaruh yang sebanding antar provinsi. Terakhir, matriks bobot tersebut dikonversi menjadi objek listw yang siap digunakan dalam fungsi-fungsi analisis spasial seperti uji Moran’s I dan estimasi model regresi spasial.
# Baca jarak kereta dari sheet Railway Distance
rail <- read_excel("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/CHINA/Data_China.xlsx", sheet = "Railway Distance", skip = 1)New names:
• `` -> `...1`rail_matrix <- as.matrix(rail[, 3:31])
rownames(rail_matrix) <- colnames(rail_matrix)
# Inverse distance (w_ij = 1/d_ij, diagonal = 0)
W_rail <- 1 / rail_matrix
diag(W_rail) <- 0
# Row-standardize
W_rail_std <- W_rail / rowSums(W_rail)
# Buat listw object
W <- mat2listw(W_rail_std, style = "W")5. Uji Morans’I
Uji Diagnostik Autokorelasi Spasial Global
Untuk memahami pola hubungan spasial antara tingkat urbanisasi (Level of Urbanization) dan pertumbuhan ekonomi daerah (Gross Regional Product - GRP) di China, dilakukan pengujian autokorelasi spasial pada sisaan (residuals) model OLS dengan menggunakan Uji Moran.
col.moran <- lm.morantest(china.lm, W)
col.moran
Global Moran I for regression residuals
data:
model: lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
weights: W
Moran I statistic standard deviate = 1.1234, p-value = 0.1306
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
0.010629595 -0.033552435 0.001546711 Berdasarkan hasil uji Global Moran’s I terhadap sisaan model OLS, diperoleh nilai koefisien yang tidak signifikan (p-value>0.05). Hasil ini mengindikasikan bahwa sisaan dari model regresi menyebar secara acak dan tidak memperlihatkan adanya pola mengelompok (clustering) yang terstruktur secara spasial di sepanjang jaringan kereta api China.
6. Uji Lagrange Multiplier
Pengujian Lagrange Multiplier (LM Test). Baik uji LM-Lag maupun LM-Error menghasilkan nilai yang tidak signifikan (p-value>0.05). Tidak signifikannya uji LM-Lag menandakan tidak adanya efek ketergantungan spasial global pada variabel dependen (Global Spatial Lag Dependent). Artinya, tinggi rendahnya nilai GRP di suatu wilayah tidak didorong oleh efek penularan atau umpan balik sistemik dari nilai GRP wilayah lain di sepanjang koridor kereta api tersebut. Sementara itu, nilai LM-Error yang tidak signifikan mengonfirmasi tidak adanya bias dari variabel luar (omitted variables) yang polanya terstruktur secara spasial pada komponen error.
china.lagrange1 <- lm.LMtests(china.lm, W, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))Please update scripts to use lm.RStests in place of lm.LMtestssummary(china.lagrange1) Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
dependence
data:
model: lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
test weights: listw
statistic parameter p.value
RSerr 0.0277967 1 0.8676
adjRSerr 0.0726855 1 0.7875
RSlag 0.0035418 1 0.9525
adjRSlag 0.0484307 1 0.8258
SARMA 0.0762273 2 0.9626Oleh karena itu, jalur pemodelan spasial global yang mengasumsikan efek multi-arah tak terbatas (seperti Spatial Autoregressive [SAR], Spatial Error Model [SEM], maupun Spatial Durbin Model [SDM]) diputuskan tidak relevan untuk menjelaskan hubungan Urbanisasi dan GRP dalam penelitian ini.
7. SLX Model
Pemodelan Efek Limpahan Lokal: Spatially-Lagged X Model (SLX)
Meskipun interaksi spasial bersifat global tidak terdeteksi (Uji Morans’I tidak signifikan), analisis diarahkan pada model SLX (Spatially-lagged X). Model SLX ini mengasumsikan bahwa performa ekonomi suatu wilayah (Y: GRP) dipengaruhi secara langsung oleh tingkat urbanisasi internalnya sendiri, serta mendapat limpahan dari tingkat urbanisasi yang terjadi di 6 wilayah tetangga terdekat yang terhubung secara fungsional melalui jalur kereta api (WX: Lag Level of Urbanization).
china.SLX <- lmSLX(GRP ~ Level_of_Urbanization, data = china_join, listw = W, Durbin = TRUE);
summary(china.SLX)
Call:
lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
data = as.data.frame(x), weights = weights)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.335e+04 4.130e+04 -5.654e-01 5.766e-01
Level_of_Urbanization 1.387e+03 1.139e+02 1.217e+01 3.074e-12
lag.Level_of_Urbanization -1.418e+02 6.801e+02 -2.086e-01 8.364e-01Level_of_Urbanization berpengaruh positif dan signifikan secara statistik (p < 0,001). Artinya, setiap kenaikan satu satuan tingkat urbanisasi di suatu provinsi akan meningkatkan variabel dependen sebesar 1.387 satuan secara langsung (direct effect). Ini menunjukkan bahwa urbanisasi di wilayah itu sendiri memiliki peran nyata terhadap variabel yang dimodelkan.
Sebaliknya, lag.Level_of_Urbanization yang merepresentasikan rata-rata tingkat urbanisasi provinsi-provinsi tetangga memiliki koefisien negatif (-141,8) namun tidak signifikan (pvalue= 0,836>0.05). Artinya, tingkat urbanisasi di provinsi-provinsi sekitar tidak terbukti memberikan pengaruh spasial (spill-over) yang berarti terhadap variabel dependen.
Untuk membuktikan apakah penambahan komponen spasial lokal (WX) pada model SLX memberikan peningkatan performa yang bermakna dibandingkan model OLS biasa, dilakukan uji formal Likelihood Ratio (LR) Test
8. LR test : SLX vs OLS Regression
# log likelihood masing-masing model
logLik(china.SLX)'log Lik.' -295.1173 (df=4)logLik(china.lm)'log Lik.' -295.1415 (df=3)# nilai LR
LR <- 2 * (as.numeric(logLik(china.SLX)) -
as.numeric(logLik(china.lm)))
# derajat bebas
df <- attr(logLik(china.SLX), "df") -
attr(logLik(china.lm), "df")
# p-value
pchisq(LR, df = df, lower.tail = FALSE)[1] 0.8257453AICc(china.lm)[1] 597.2431AICc(china.SLX)[1] 599.9012Berdasarkan uji perbandingan tersebut, diperoleh nilai Likelihood Ratio yang tidak signifikan (p-value>0.05) yang disertai dengan nilai informasi kriteria AICc yang tidak lebih kecil dibandingkan OLS regression. Hasil statistik formal ini menyatakan gagal tolak H0, yang berarti penambahan variabel tingkat urbanisasi tetangga berbasis konektivitas kereta api (WX) tidak memberikan kontribusi signifikan dalam menurunkan AICc. Oleh karena itu, model regresi OLS menjadi model terbaik dalam memodelkan pengaruh Level_of_Urbanization terhadap GRP.
B. Model Dependensi Spasial (Jawa Barat)
Data Kemiskinan Jawa Barat
1. Data Preprocessing
shp_jabar <- st_read("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/JABAR/Jabar2/Jabar2.shp")Reading layer `Jabar2' from data source
`/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/JABAR/Jabar2/Jabar2.shp'
using driver `ESRI Shapefile'
Simple feature collection with 26 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
CRS: NAdata_jabar <- read_excel("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 11/JABAR/data_jabar.xlsx")data_jabar$KABKOTNO <- as.character(data_jabar$KABKOTNO)
jabar <- shp_jabar %>%
left_join(data_jabar, by = "KABKOTNO") %>%
na.omit()colnames(jabar) [1] "PROVNO.x" "KABKOTNO" "KODE2010.x" "PROVINSI.x" "KABKOT.x"
[6] "SUMBER" "IDSP2010.x" "PROVNO.y" "KODE2010.y" "PROVINSI.y"
[11] "KABKOT.y" "IDSP2010.y" "Long" "Lat" "p.miskin15"
[16] "p.miskin16" "j.miskin15" "j.miskin16" "AHH2015" "AHH2016"
[21] "EYS2015" "EYS2016" "MYS2015" "MYS2016" "EXP2015"
[26] "EXP2016" "APM.SD15" "APM.SMP15" "APM.SMA15" "APM.PT15"
[31] "APK.SD15" "APK.SMP15" "APK.SMA15" "APK.PT15" "APS.USIA15"
[36] "APS.USIA2" "APS.USIA3" "APS.USIA4" "geometry" 2. Exploratory Analysis
plot(jabar[,c(
"p.miskin15",
"AHH2015",
"EYS2015",
"MYS2015",
"EXP2015",
"APM.SD15",
"APM.SMP15",
"APM.SMA15",
"APM.PT15",
"APK.SD15",
"APK.SMP15",
"APK.SMA15",
"APK.PT15",
"APS.USIA15",
"APS.USIA2",
"APS.USIA3",
"APS.USIA4"
)])Warning: plotting the first 9 out of 17 attributes; use max.plot = 17 to plot
all
3. Matriks Pembobot Spasial
Membuat matriks bobot spasial berbasis ketetanggaan (contiguity) dari data poligon wilayah Jawa Barat. Fungsi poly2nb secara default menggunakan metode Queen dan dua wilayah dianggap bertetangga apabila berbagi batas sisi maupun titik sudut untuk mendeteksi secara otomatis wilayah-wilayah yang saling bersentuhan batas administratifnya. Kemudian dikonversi menjadi matriks bobot spasial menggunakan nb2listwdengan metode row standardization (style = "W"), sehingga total bobot dari semua tetangga untuk setiap wilayah bernilai satu dan pengaruh antar wilayah berada pada skala yang sebanding. Opsi zero.policy = TRUE memastikan proses tetap berjalan meskipun terdapat wilayah yang tidak memiliki tetangga sama sekali. Matriks bobot inilah yang selanjutnya digunakan sebagai dasar dalam uji autokorelasi spasial maupun estimasi model regresi spasial untuk data Jawa Barat.
# =====================================================
# SPATIAL WEIGHT
# =====================================================
nb_jabar <- suppressWarnings(poly2nb(jabar))
lw_jabar <- nb2listw(
nb_jabar,
style = "W",
zero.policy = TRUE
)4. Penentuan Variabel (Indeks Moran)
Sebelum dilakukan pemodelan, terlebih dulu mencari kombinasi variabel independen terbaik yang menghasilkan autokorelasi spasial paling kuat pada residual model regresi OLS, sebagai langkah awal dalam pemilihan variabel sebelum estimasi model spasial. Prosesnya dimulai dengan mendefinisikan 17 kandidat variabel yang mencakup indikator kemiskinan, kesehatan, pendidikan, dan pengeluaran di Jawa Barat tahun 2015, dengan p.miskin15 (persentase penduduk miskin) ditetapkan sebagai variabel dependen (Y) dan seluruh variabel lainnya sebagai kandidat prediktor (X). Selanjutnya, kode secara sistematis membentuk semua kemungkinan kombinasi variabel X mulai dari 1 hingga maksimum 6 variabel, lalu untuk setiap kombinasi diestimasi model OLS dan residualnya diuji menggunakan Moran’s I terhadap matriks bobot spasial lw_jabar. Hanya kombinasi yang menghasilkan p-value Moran’s I < 0,05 — artinya residualnya terbukti mengandung autokorelasi spasial yang signifikan — yang disimpan ke dalam dataframe hasil. Setelah seluruh kombinasi selesai dievaluasi, hasilnya diurutkan berdasarkan p-value terkecil, dan kombinasi dengan p-value terkecil dipilih sebagai model terbaik, karena mengindikasikan adanya dependensi spasial paling kuat yang perlu ditangani melalui pemodelan regresi spasial lebih lanjut.
# =========================================================
# ANALISIS MORAN'S I RESIDUAL
# MENCARI KOMBINASI X TERBAIK
# DENGAN P-VALUE TERKECIL
# =========================================================
# =========================================================
# 5. DAFTAR VARIABEL
# =========================================================
variabel_list <- c(
"p.miskin15",
"AHH2015",
"EYS2015",
"MYS2015",
"EXP2015",
"APM.SD15",
"APM.SMP15",
"APM.SMA15",
"APM.PT15",
"APK.SD15",
"APK.SMP15",
"APK.SMA15",
"APK.PT15",
"APS.USIA15",
"APS.USIA2",
"APS.USIA3",
"APS.USIA4"
)
# =========================================================
# 6. TENTUKAN VARIABEL Y
# =========================================================
variabel_Y <- "p.miskin15"
# variabel kandidat X
variabel_X <- setdiff(
variabel_list,
variabel_Y
)
# =========================================================
# 7. DATAFRAME HASIL
# =========================================================
hasil_model <- data.frame(
Formula = character(),
Jumlah_X = numeric(),
Moran_I = numeric(),
P_Value = numeric(),
stringsAsFactors = FALSE
)
# =========================================================
# 8. LOOP KOMBINASI VARIABEL
# =========================================================
maksimum_X <- 6
for (k in 1:maksimum_X) {
# membuat kombinasi
kombinasi_x <- combn(
variabel_X,
k,
simplify = FALSE
)
# loop tiap kombinasi
for (x_set in kombinasi_x) {
tryCatch({
# =====================================================
# MEMBUAT FORMULA
# =====================================================
formula_text <- paste(
variabel_Y,
"~",
paste(x_set, collapse = " + ")
)
rumus <- as.formula(formula_text)
# =====================================================
# MODEL OLS
# =====================================================
model_ols <- lm(
rumus,
data = jabar
)
# =====================================================
# MORAN TEST RESIDUAL
# =====================================================
uji_moran <- lm.morantest(
model_ols,
lw_jabar,
zero.policy = TRUE
)
# =====================================================
# AMBIL NILAI
# =====================================================
moran_i <- as.numeric(
uji_moran$estimate["Observed Moran I"]
)
p_val <- uji_moran$p.value
# =====================================================
# SIMPAN JIKA SIGNIFIKAN
# =====================================================
if (p_val < 0.05) {
hasil_model <- rbind(
hasil_model,
data.frame(
Formula = formula_text,
Jumlah_X = k,
Moran_I = moran_i,
P_Value = p_val
)
)
}
}, error = function(e) {
# abaikan error
NULL
})
}
}
# =========================================================
# 9. URUTKAN HASIL
# =========================================================
hasil_model <- hasil_model %>%
arrange(P_Value)
# =========================================================
# 11. MODEL TERBAIK
# P-VALUE PALING KECIL
# =========================================================
if (nrow(hasil_model) == 0) {
cat("Tidak ada model signifikan.\n")
} else {
# ambil model terbaik
model_terbaik <- hasil_model %>%
slice(1)
cat(
"Formula : ",
as.character(model_terbaik$Formula),
"\n",
sep = ""
)
cat(
"Jumlah Variabel X : ",
model_terbaik$Jumlah_X,
"\n",
sep = ""
)
cat(
"Moran's I : ",
round(model_terbaik$Moran_I, 5),
"\n",
sep = ""
)
cat(
"P-Value : ",
format(
model_terbaik$P_Value,
scientific = FALSE
),
"\n",
sep = ""
)
}Formula : p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15
Jumlah Variabel X : 5
Moran's I : 0.5012
P-Value : 0.005018003Berdasarkan hasil pencarian kombinasi variabel terbaik, terpilih model dengan formula p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15 yang terdiri dari 5 variabel independen. Model ini menghasilkan nilai Moran’s I sebesar 0,5012 dengan p-value 0,005 — signifikan pada taraf 1% — yang menunjukkan bahwa residual model OLS dengan kombinasi variabel ini mengandung autokorelasi spasial positif yang kuat. Artinya, wilayah-wilayah di Jawa Barat dengan pola kemiskinan yang serupa cenderung berkelompok secara geografis dan belum sepenuhnya dijelaskan oleh kelima variabel tersebut secara biasa. Kombinasi ini dipilih sebagai model terbaik karena menghasilkan p-value Moran’s I terkecil di antara seluruh kombinasi yang diuji, sehingga menjadi justifikasi terkuat untuk melanjutkan analisis menggunakan model regresi spasial yang mampu mengakomodasi struktur dependensi spasial tersebut secara eksplisit.
5. OLS Regression
jabar.lm <- lm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15, data=jabar)
summary(jabar.lm)
Call:
lm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 +
APS.USIA15, data = jabar)
Residuals:
Min 1Q Median 3Q Max
-6.2541 -1.5129 0.7468 1.7253 2.9766
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -263.5186 236.1804 -1.116 0.2883
EYS2015 2.7220 2.1113 1.289 0.2238
APM.SD15 -0.5289 0.6939 -0.762 0.4619
APK.SD15 0.4526 0.2946 1.536 0.1527
APK.PT15 -0.3228 0.1731 -1.865 0.0891 .
APS.USIA15 2.4849 2.5164 0.987 0.3446
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.114 on 11 degrees of freedom
Multiple R-squared: 0.5749, Adjusted R-squared: 0.3816
F-statistic: 2.975 on 5 and 11 DF, p-value: 0.061276. Uji Lagrange Multiplier
jabar.lagrange <- lm.LMtests(jabar.lm, lw_jabar, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))Please update scripts to use lm.RStests in place of lm.LMtestssummary(jabar.lagrange) Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
dependence
data:
model: lm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 +
APK.PT15 + APS.USIA15, data = jabar)
test weights: listw
statistic parameter p.value
RSerr 7.23955 1 0.007132 **
adjRSerr 5.48939 1 0.019132 *
RSlag 2.01129 1 0.156133
adjRSlag 0.26113 1 0.609346
SARMA 7.50067 2 0.023510 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hasil uji LM menunjukkan bahwa dependensi spasial dalam data Jawa Barat bersumber dari komponen error, bukan dari lag spasial variabel dependen. Uji RSerr signifikan (p = 0,007) mengindikasikan adanya autokorelasi spasial pada residual model, dan ketika dikontrol terhadap pengaruh lag, adjRSerr tetap signifikan (p = 0,019) — membuktikan bahwa spatial error hadir secara mandiri dalam model.
Sebaliknya, uji RSlag tidak signifikan (p = 0,156) dan adjRSlag pun tidak signifikan (p = 0,609), yang berarti tidak ada bukti bahwa tingkat kemiskinan di suatu wilayah dipengaruhi langsung oleh tingkat kemiskinan wilayah tetangganya.
5. Spatial Error Model (SEM)
jabar.err <- errorsarlm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar,data=jabar, zero.policy = TRUE)
summary(jabar.err)
Call:errorsarlm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 +
APK.PT15 + APS.USIA15, data = jabar, listw = lw_jabar, zero.policy = TRUE)
Residuals:
Min 1Q Median 3Q Max
-1.96720 -0.97706 0.17825 0.76731 1.45103
Type: error
Regions with no neighbours included:
18 19 21 25 26
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -425.458866 84.564213 -5.0312 4.874e-07
EYS2015 4.311729 1.035451 4.1641 3.126e-05
APM.SD15 -1.521207 0.226462 -6.7173 1.852e-11
APK.SD15 0.640932 0.101930 6.2880 3.216e-10
APK.PT15 -0.358322 0.069408 -5.1626 2.436e-07
APS.USIA15 4.680477 0.828208 5.6513 1.592e-08
Lambda: 0.9059, LR test value: 23.088, p-value: 1.5476e-06
Asymptotic standard error: 0.060606
z-value: 14.947, p-value: < 2.22e-16
Wald statistic: 223.43, p-value: < 2.22e-16
Log likelihood: -28.18863 for error model
ML residual variance (sigma squared): 1.0473, (sigma: 1.0234)
Number of observations: 17
Number of parameters estimated: 8
AIC: 72.377, (AIC for lm: 93.465)6. Model SLX
jabar.SLX <- lmSLX(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar,data=jabar, zero.policy = TRUE, Durbin = TRUE);
summary(jabar.SLX)
Call:
lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
data = as.data.frame(x), weights = weights)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -145.59982 307.35728 -0.47372 0.65245
EYS2015 2.63817 4.59362 0.57431 0.58663
APM.SD15 0.26488 0.97666 0.27122 0.79532
APK.SD15 0.59059 0.23882 2.47293 0.04827
APK.PT15 -0.20577 0.38277 -0.53759 0.61021
APS.USIA15 0.36388 3.17195 0.11472 0.91241
lag.EYS2015 -4.02200 9.31944 -0.43157 0.68113
lag.APM.SD15 3.42980 1.21392 2.82538 0.03014
lag.APK.SD15 -0.41090 0.69159 -0.59414 0.57412
lag.APK.PT15 0.38003 0.82137 0.46268 0.65990
lag.APS.USIA15 -2.48727 1.03245 -2.40910 0.05264Hasil regresi menggunakan model SLX, menunjukkan bahwa lag pada variabel X tidak berpengaruh signifikan. Model SEM tidak perlu ditambahakan depensi spasial variabel X, sehingga SDEM tidak relevan lagi.
7. Model SARMA
jabar.sarma <- sacsarlm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar, data=jabar, zero.policy = TRUE)
summary(jabar.sarma)
Call:sacsarlm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 +
APK.PT15 + APS.USIA15, data = jabar, listw = lw_jabar, zero.policy = TRUE)
Residuals:
Min 1Q Median 3Q Max
-1.837342 -0.900247 -0.050073 0.685523 1.606525
Type: sac
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -442.074285 83.148799 -5.3167 1.057e-07
EYS2015 4.793132 1.140231 4.2037 2.626e-05
APM.SD15 -1.707710 0.293038 -5.8276 5.623e-09
APK.SD15 0.672837 0.103147 6.5231 6.886e-11
APK.PT15 -0.392055 0.076727 -5.1097 3.226e-07
APS.USIA15 4.939024 0.838588 5.8897 3.869e-09
Rho: 0.20377
Asymptotic standard error: 0.19972
z-value: 1.0202, p-value: 0.30761
Lambda: 0.91256
Asymptotic standard error: 0.065071
z-value: 14.024, p-value: < 2.22e-16
LR test value: 23.918, p-value: 6.4015e-06
Log likelihood: -27.77358 for sac model
ML residual variance (sigma squared): 0.97176, (sigma: 0.98578)
Number of observations: 17
Number of parameters estimated: 9
AIC: 73.547, (AIC for lm: 93.465)Model SARMA (tipe sac) menghasilkan temuan yang menarik pada kedua parameter spasialnya. Nilai Rho (ρ) sebesar 0,204 tidak signifikan (p = 0,308), yang berarti tidak terbukti adanya pengaruh langsung dari tingkat kemiskinan wilayah tetangga terhadap kemiskinan suatu wilayah. Sebaliknya, Lambda (λ) sebesar 0,913 sangat signifikan (p < 2,22e-16), mengindikasikan adanya autokorelasi spasial yang sangat kuat pada komponen error. Faktor-faktor yang tidak teramati dalam model cenderung berkorelasi kuat antar wilayah yang berdekatan.
Pola ini ρ tidak signifikan, λ sangat signifikan sehingga mengonfirmasi bahwa SEM lebih tepat dibandingkan SARMA untuk data in. Hal ini juga akan dibuktikan dengan melakukan Uji Likelihood Ratio Test antara Model SEM dan Model SARMA
8. LR Test : SEM vs SARMA
anova(jabar.err,jabar.sarma) Model df AIC logLik Test L.Ratio p-value
jabar.err 1 8 72.377 -28.189 1
jabar.sarma 2 9 73.547 -27.774 2 0.83008 0.36225Hasil menunjukkan nilai L.Ratio sebesar 0,830 dengan p-value 0,362 — jauh di atas taraf signifikansi 0,05 — sehingga H₀ gagal ditolak. Artinya, penambahan parameter lag spasial (ρ) pada SARMA tidak memberikan peningkatan fit model yang signifikan dibandingkan SEM. Hal ini juga tercermin dari nilai AIC SEM (72,377) yang lebih kecil dari SARMA (73,547), mengonfirmasi bahwa SEM lebih parsimoni tanpa mengorbankan kualitas model.
9. Akaike Information Criterion Corrected (AICc)
AICc(jabar.lm)[1] 105.9096AICc(jabar.err)[1] 90.37725AICc(jabar.sarma)[1] 99.26145