# Load Package
suppressPackageStartupMessages({
library(sf)
library(spdep)
library(spatialreg)
library(tidyverse)
library(MuMIn)
library(dplyr)
library(car)
})Marginal Effects (Spill-over) on the Spatial Regression Modeling
Studi Kasus 1 : Data chi.poly
Data chi.poly akan digunakan sebagai bahan analisis dalam latihan ini. Proses pemodelan dimulai dengan tahap diagnostik, yakni pengujian multikolinieritas antar peubah bebas melalui pendekatan VIF guna memastikan tidak terdapat hubungan linear yang kuat di antara prediktor yang digunakan. Di samping itu, autokorelasi spasial pada model juga akan ditelusuri dengan memanfaatkan matriks pembobot berbasis jarak W_dist, sehingga dapat diketahui apakah terdapat ketergantungan spasial yang perlu diperhitungkan dalam pemilihan model. Berdasarkan hasil kedua pengujian tersebut, model terbaik akan ditentukan dan diinterpretasikan secara komprehensif.
1. Persiapan Data
# Membaca data spasial foreclosures dalam format shapefile (.shp)
chi.poly <- st_read("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 12/poly/foreclosures.shp")Reading layer `foreclosures' from data source
`/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 12/poly/foreclosures.shp'
using driver `ESRI Shapefile'
Simple feature collection with 897 features and 16 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -87.94027 ymin: 41.64423 xmax: -87.52404 ymax: 42.03793
CRS: NA2. Eksplorasi dan Pra-pemrosesan Data
# Melihat Struktur Data
str(chi.poly)Classes 'sf' and 'data.frame': 897 obs. of 17 variables:
$ SP_ID : chr "1" "2" "3" "4" ...
$ fips : chr "17031010100" "17031010200" "17031010300" "17031010400" ...
$ est_fcs : int 43 129 55 21 64 56 107 43 7 51 ...
$ est_mtgs : int 904 2122 1151 574 1427 1241 1959 830 208 928 ...
$ est_fcs_rt: num 4.76 6.08 4.78 3.66 4.48 4.51 5.46 5.18 3.37 5.5 ...
$ res_addr : int 2530 3947 3204 2306 5485 2994 3701 1694 443 1552 ...
$ est_90d_va: num 12.61 12.36 10.46 5.03 8.44 ...
$ bls_unemp : num 8.16 8.16 8.16 8.16 8.16 8.16 8.16 8.16 8.16 8.16 ...
$ county : chr "Cook County" "Cook County" "Cook County" "Cook County" ...
$ fips_num : num 1.7e+10 1.7e+10 1.7e+10 1.7e+10 1.7e+10 ...
$ totpop : int 5391 10706 6649 5325 10944 7178 10799 5403 1089 3634 ...
$ tothu : int 2557 3981 3281 2464 5843 3136 3875 1768 453 1555 ...
$ huage : int 61 53 56 60 54 58 48 57 61 48 ...
$ oomedval : int 169900 147000 119800 151500 143600 145900 153400 170500 215900 114700 ...
$ property : num 646 914 478 509 641 612 678 332 147 351 ...
$ violent : num 433 421 235 159 240 266 272 146 78 84 ...
$ geometry :sfc_POLYGON of length 897; first list element: List of 1
..$ : num [1:15, 1:2] -87.7 -87.7 -87.7 -87.7 -87.7 ...
..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "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] "SP_ID" "fips" "est_fcs" "est_mtgs" ...2.1 Deskripsi Data chi.poly
Data chi.poly merupakan data spasial berbentuk poligon (sf/spatial data frame) yang menggambarkan kondisi sosial-ekonomi di wilayah Cook County, Chicago, Illinois, Amerika Serikat. Data ini terdiri dari 897 observasi
Pada studi kasus ini, bertujuan untuk membuat model regresi dependensi spasial dan efek marginal dari faktor-faktor sosial-ekonomi yang mempengaruhi tingkat foreclosure di wilayah Cook County, Chicago. Tingkat foreclosure yang tinggi pada suatu wilayah tidak terjadi secara acak, melainkan diduga dipengaruhi oleh berbagai kondisi lingkungan sekitarnya, seperti tingkat pengangguran, nilai properti, kepadatan hunian, usia bangunan, serta tingkat kriminalitas. Dengan memanfaatkan data spasial berbasis unit sensus tract, analisis ini tidak hanya menelaah hubungan antara peubah bebas dan peubah respon secara klasik, tetapi juga mempertimbangkan struktur ketergantungan spasial yang mungkin terbentuk akibat kedekatan geografis antar wilayah. Oleh karena itu, dilakukan eksplorasi data pada est_fcs_rt berikut:
2.2 Pengecekan Missing Values dan Imputasi
colSums(is.na(chi.poly)) SP_ID fips est_fcs est_mtgs est_fcs_rt res_addr est_90d_va
0 0 0 0 0 0 0
bls_unemp county fips_num totpop tothu huage oomedval
0 0 0 0 0 0 0
property violent geometry
3 0 0 # Buat nb dulu dari chi.poly
nb_queen <- poly2nb(chi.poly, queen = TRUE)
# Baru imputasi
idx_na <- which(is.na(chi.poly$property))
for (i in idx_na) {
tetangga <- nb_queen[[i]]
chi.poly$property[i] <- mean(chi.poly$property[tetangga],
na.rm = TRUE)
}
sum(is.na(chi.poly$property)) # harus 0[1] 02.3 Visualisasi Data
# Visualisasi Rasio Tingkat Penyitaan Rumah
plot(chi.poly["est_fcs_rt"])
Eksplorasi awal dilakukan dengan memvisualisasikan sebaran spasial variabel tingkat foreclosure (est_fcs_rt) pada 897 unit sensus tract di wilayah Cook County, Chicago. Berdasarkan peta yang dihasilkan, sebagian besar wilayah menunjukkan tingkat foreclosure yang rendah (0–5), tercermin dari dominasi warna biru tua pada hampir seluruh area. Namun demikian, terdapat pengelompokan wilayah dengan nilai yang lebih tinggi (10–20) yang terkonsentrasi di bagian tengah hingga selatan kota, serta beberapa titik ekstrem dengan nilai di atas 30 yang tersebar secara terisolasi sebagai pencilan spasial.
Pola sebaran yang tidak merata ini memberikan indikasi awal bahwa terdapat ketergantungan spasial dalam data, sehingga pemodelan yang mengabaikan efek spasial berpotensi menghasilkan estimasi yang bias. Oleh karena itu, tahap selanjutnya akan melibatkan pengujian autokorelasi spasial melalui statistik Moran’s I serta pemilihan model regresi spasial yang paling sesuai dengan karakteristik data.
2.4 Standarisasi Data
Standardisasi variabel dalam regresi spasial penting dilakukan karena variabel-variabel yang digunakan dalam model seringkali memiliki satuan dan skala pengukuran yang sangat berbeda — misalnya tingkat pengangguran dalam persen, jumlah unit hunian dalam ribuan, dan tingkat kejahatan dalam angka absolut. Perbedaan skala ini dapat menyebabkan koefisien regresi tidak bisa dibandingkan secara langsung, sehingga sulit untuk menilai variabel mana yang memiliki pengaruh relatif lebih besar terhadap variabel dependen. Dengan melakukan standardisasi (misalnya menjadi z-score), seluruh variabel berada pada skala yang seragam sehingga koefisien yang dihasilkan mencerminkan pengaruh relatif yang sesungguhnya dan dapat dibandingkan satu sama lain secara adil. Selain itu, dalam konteks regresi spasial yang melibatkan matriks bobot dan perhitungan autokorelasi, perbedaan skala yang ekstrem dapat mempengaruhi stabilitas numerik dalam proses estimasi. Oleh karena itu, standardisasi tidak hanya membantu interpretasi, tetapi juga berkontribusi pada kestabilan dan keandalan hasil estimasi model spasial secara keseluruhan.
# Drop geometry dulu, lalu pilih kolom yang ada
chi_df <- st_drop_geometry(chi.poly)
# Pastikan vars_std hanya yang benar-benar ada di data
vars_std <- c("est_fcs_rt", "est_mtgs",
"oomedval", "bls_unemp", "est_90d_va", "res_addr",
"tothu", "totpop", "violent", "property", "huage")
# Cek variabel mana yang ada
vars_std <- vars_std[vars_std %in% names(chi_df)]
# Ambil mean dan SD
stats_asli <- data.frame(
variabel = vars_std,
mean = sapply(chi_df[, vars_std], function(x) mean(as.numeric(x), na.rm = TRUE)),
sd = sapply(chi_df[, vars_std], function(x) sd(as.numeric(x), na.rm = TRUE))
)
print(stats_asli) variabel mean sd
est_fcs_rt est_fcs_rt 6.729710e+00 3.760875e+00
est_mtgs est_mtgs 6.606198e+02 6.239358e+02
oomedval oomedval 1.596818e+05 9.293880e+04
bls_unemp bls_unemp 8.015942e+00 9.352458e-01
est_90d_va est_90d_va 6.766232e+00 5.936992e+00
res_addr res_addr 1.197284e+03 1.152136e+03
tothu tothu 1.353851e+03 1.113852e+03
totpop totpop 3.402415e+03 2.584909e+03
violent violent 1.693389e+02 1.897759e+02
property property 4.199578e+02 3.441708e+02
huage huage 5.114270e+01 1.159278e+01# Baru jalankan standardisasi
chi.poly <- chi.poly %>%
mutate(across(
c(est_fcs, est_mtgs, est_fcs_rt,
oomedval, bls_unemp, est_90d_va, res_addr,
totpop, violent, property, huage),
~ as.numeric(scale(.))
))2.5 Pengecekan Multikolinieritas
df_vif <- st_drop_geometry(chi.poly)
model_ols <- lm(est_fcs_rt ~ est_mtgs + est_90d_va + oomedval+ res_addr + bls_unemp + totpop + tothu + huage + property + violent, data = df_vif)
# Menghitung nilai VIF
vif_values <- vif(model_ols)
# Menampilkan hasil VIF
print(vif_values) est_mtgs est_90d_va oomedval res_addr bls_unemp totpop tothu
3.719782 1.426790 1.390321 13.406935 1.429911 8.295513 20.308236
huage property violent
1.536820 3.374562 3.727998 vif_df <- data.frame(Variable = names(vif_values), VIF = vif_values)
vif_df <- vif_df[order(-vif_df$VIF), ] # Urutkan dari VIF tertinggi
print(vif_df) Variable VIF
tothu tothu 20.308236
res_addr res_addr 13.406935
totpop totpop 8.295513
violent violent 3.727998
est_mtgs est_mtgs 3.719782
property property 3.374562
huage huage 1.536820
bls_unemp bls_unemp 1.429911
est_90d_va est_90d_va 1.426790
oomedval oomedval 1.390321VIF digunakan untuk mendeteksi adanya multikolinearitas antar variabel independen dalam model. Dengan menggunakan batas kritis VIF > 5, hasil uji menunjukkan bahwa terdapat tiga variabel yang bermasalah. Variabel tothu (total unit hunian) memiliki nilai VIF tertinggi sebesar 20,31, diikuti oleh res_addr dengan VIF 13,41, dan totpop dengan VIF 8,30 — ketiganya melampaui batas 5 dan mengindikasikan adanya multikolinearitas yang cukup serius. Artinya, ketiga variabel tersebut berkorelasi kuat dengan variabel lain dalam model. Sebaliknya, tujuh variabel lainnya yakni violent, est_mtgs, property, huage, bls_unemp, est_90d_va, dan oomedvalseluruhnya memiliki nilai VIF di bawah 5, yang berarti bebas dari masalah multikolinearitas.
Berdasarkan hal ini, variabel tothu, res_addr, dan totpop perlu dipertimbangkan untuk dipilih salah satu diantaranya, dalam hal ini tothu yang dimasukkan ke dalam model.
3. Matriks Pembobot Spasial
coords <- st_coordinates(st_centroid(st_geometry(chi.poly)))
W_dist<-dnearneigh(coords,0,5,longlat = TRUE)
W_dist_listw <- nb2listw(W_dist, style = "W", zero.policy = TRUE)Matriks bobot spasial digunakan sebagai dasar dalam analisis dependensi spasial. Prosesnya dimulai dengan mengekstrak titik pusat (centroid) dari setiap poligon wilayah di Chicago, kemudian menentukan hubungan ketetanggaan antar wilayah berdasarkan jarak. Dua wilayah dianggap bertetangga apabila jarak antar centroid-nya tidak melebihi 5 km (radius agar semua mempunyai tetangga). Hasil daftar tetangga tersebut selanjutnya dikonversi menjadi matriks bobot spasial yang dinormalisasi per baris (row-standardized), sehingga total pengaruh dari semua tetangga untuk setiap wilayah bernilai satu. Matriks inilah yang kemudian digunakan dalam uji Moran’s I maupun pemodelan regresi spasial untuk menangkap sejauh mana suatu wilayah dipengaruhi oleh wilayah-wilayah di sekitarnya.
4. OLS Regression
foreclosure.lm <- lm(est_fcs_rt~ est_mtgs + bls_unemp+ est_90d_va + tothu+ violent + property, data = chi.poly)
summary(foreclosure.lm)
Call:
lm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
tothu + violent + property, data = chi.poly)
Residuals:
Min 1Q Median 3Q Max
-3.1288 -0.4272 -0.0082 0.4192 7.7340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.286e-01 5.036e-02 6.524 1.15e-10 ***
est_mtgs 1.347e-01 3.839e-02 3.508 0.000475 ***
bls_unemp 1.600e-01 2.406e-02 6.649 5.14e-11 ***
est_90d_va 2.340e-01 2.708e-02 8.642 < 2e-16 ***
tothu -2.427e-04 3.291e-05 -7.374 3.80e-13 ***
violent 7.267e-01 3.852e-02 18.866 < 2e-16 ***
property -2.844e-01 4.123e-02 -6.898 1.00e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.703 on 890 degrees of freedom
Multiple R-squared: 0.5091, Adjusted R-squared: 0.5058
F-statistic: 153.8 on 6 and 890 DF, p-value: < 2.2e-165. Uji Morans’I
lm.morantest(foreclosure.lm,W_dist_listw)
Global Moran I for regression residuals
data:
model: lm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
tothu + violent + property, data = chi.poly)
weights: W_dist_listw
Moran I statistic standard deviate = 43.543, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
2.034396e-01 -2.129529e-03 2.228843e-05 Uji Global Moran’s I bertujuan untuk mengetahui apakah residual dari model OLS masih mengandung pola spasial atau sudah tersebar secara acak. Jika residual masih berpola, berarti model OLS belum mampu menjelaskan pengaruh lokasi dalam data.
Hasil uji menunjukkan nilai Moran’s I sebesar 0,2034 yang signifikan secara statistik (p-value < 2,2e-16). Nilai ini jauh di atas ekspektasi ketika tidak ada autokorelasi (≈ 0), yang berarti H₀ ditolak — residual tidak tersebar acak, melainkan masih mengelompok berdasarkan lokasi geografis. Wilayah yang residualnya tinggi cenderung berdekatan dengan wilayah yang residualnya juga tinggi, dan begitu sebaliknya.
Hal ini menjadi dasar yang kuat untuk melanjutkan analisis menggunakan model regresi spasial, yang mampu mengakomodasi pengaruh kedekatan lokasi secara eksplisit.
6. Uji Lagrange Multiplier
chipoly.lagrange <- lm.RStests(foreclosure.lm, W_dist_listw, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
summary(chipoly.lagrange) Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
dependence
data:
model: lm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
tothu + violent + property, data = chi.poly)
test weights: W_dist_listw
statistic parameter p.value
RSerr 1502.14 1 < 2.2e-16 ***
adjRSerr 885.17 1 < 2.2e-16 ***
RSlag 873.99 1 < 2.2e-16 ***
adjRSlag 257.02 1 < 2.2e-16 ***
SARMA 1759.16 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Berdasarkan hasil uji LM, model SARMA menjadi kandidat utama karena mampu menangkap dependensi spasial dari dua sumber sekaligus. Namun demikian, kompleksitas SARMA perlu diimbangi dengan bukti empiris berupa nilai AIC yang lebih rendah dan residual yang terbebas dari autokorelasi spasial dibandingkan SAR maupun SEM. Jika perbedaan AIC antara SARMA dan salah satu model yang lebih sederhana tidak substansial, model yang lebih parsimoni (SAR atau SEM) dapat dipertimbangkan demi kemudahan interpretasi.
| Uji | P-Value | Interpretasi |
|---|---|---|
| RSerr (LMerr) | < 2.2e-16 | Indikasi kuat adanya spatial error dependence — error model berkorelasi antar wilayah tetangga |
| adjRSerr (RLMerr) | < 2.2e-16 | Setelah dikontrol lag, spatial error tetap signifikan secara mandiri |
| RSlag (LMlag) | < 2.2e-16 | Indikasi kuat adanya spatial lag dependence — nilai variabel dependen di suatu wilayah dipengaruhi nilai tetangganya |
| adjRSlag (RLMlag) | < 2.2e-16 | Setelah dikontrol error, spatial lag tetap signifikan secara mandiri |
| SARMA | < 2.2e-16 | Kombinasi keduanya menghasilkan statistik tertinggi, mengindikasikan model yang paling sesuai |
7. Model Regresi Dependensi Spasial
7.1 Spatial Error Model (SEM)
# Menjalankan Model SEM
chipoly.err <- errorsarlm(est_fcs_rt~ est_mtgs + bls_unemp+ est_90d_va + tothu+ violent + property, data = chi.poly, W_dist_listw, zero.policy = TRUE)
summary(chipoly.err)
Call:errorsarlm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
tothu + violent + property, data = chi.poly, listw = W_dist_listw,
zero.policy = TRUE)
Residuals:
Min 1Q Median 3Q Max
-3.064871 -0.253228 -0.017144 0.279277 7.544985
Type: error
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.0563e-01 1.1372e+00 -0.1808 0.8565019
est_mtgs 6.1451e-02 3.3611e-02 1.8283 0.0675045
bls_unemp 1.8196e-01 1.9539e-02 9.3124 < 2.2e-16
est_90d_va 2.1277e-01 2.2984e-02 9.2574 < 2.2e-16
tothu -4.3159e-05 2.9484e-05 -1.4638 0.1432365
violent 3.0197e-01 3.7576e-02 8.0363 8.882e-16
property -1.3566e-01 3.6007e-02 -3.7674 0.0001649
Lambda: 0.98326, LR test value: 354.09, p-value: < 2.22e-16
Asymptotic standard error: 0.010847
z-value: 90.651, p-value: < 2.22e-16
Wald statistic: 8217.6, p-value: < 2.22e-16
Log likelihood: -776.1307 for error model
ML residual variance (sigma squared): 0.32166, (sigma: 0.56715)
Number of observations: 897
Number of parameters estimated: 9
AIC: 1570.3, (AIC for lm: 1922.4)7.2 Spatial Lag Model (SAR)
# Menjalankan Model SAR
chipoly.sar <- lagsarlm(est_fcs_rt~ est_mtgs + bls_unemp+ est_90d_va + tothu+ violent + property, data=chi.poly, listw = W_dist_listw, zero.policy = TRUE)
# Tampilkan hasil estimasi
summary(chipoly.sar)
Call:lagsarlm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
tothu + violent + property, data = chi.poly, listw = W_dist_listw,
zero.policy = TRUE)
Residuals:
Min 1Q Median 3Q Max
-2.85852661 -0.24217208 0.00035621 0.26626261 7.52079127
Type: lag
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.4516e-02 4.0869e-02 1.5786 0.114421
est_mtgs 9.5498e-02 3.0726e-02 3.1081 0.001883
bls_unemp 1.7423e-01 1.9242e-02 9.0549 < 2.2e-16
est_90d_va 1.8867e-01 2.1835e-02 8.6408 < 2.2e-16
tothu -5.6106e-05 2.6738e-05 -2.0983 0.035875
violent 3.1349e-01 3.3040e-02 9.4881 < 2.2e-16
property -1.5873e-01 3.3036e-02 -4.8047 1.55e-06
Rho: 0.83679, LR test value: 381.67, p-value: < 2.22e-16
Asymptotic standard error: 0.02953
z-value: 28.337, p-value: < 2.22e-16
Wald statistic: 803, p-value: < 2.22e-16
Log likelihood: -762.3397 for lag model
ML residual variance (sigma squared): 0.31603, (sigma: 0.56216)
Number of observations: 897
Number of parameters estimated: 9
AIC: 1542.7, (AIC for lm: 1922.4)
LM test for residual autocorrelation
test value: 29.049, p-value: 7.0559e-087.3 Spatial Autocorrelation Model (SAC)
chipoly.sarma <- sacsarlm(est_fcs_rt~ est_mtgs + bls_unemp+ est_90d_va + violent + property , data = chi.poly, listw = W_dist_listw, zero.policy = TRUE)
summary(chipoly.sarma)
Call:sacsarlm(formula = est_fcs_rt ~ est_mtgs + bls_unemp + est_90d_va +
violent + property, data = chi.poly, listw = W_dist_listw,
zero.policy = TRUE)
Residuals:
Min 1Q Median 3Q Max
-3.0146293 -0.2251465 0.0063244 0.2527593 7.5750086
Type: sac
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.012803 0.049773 -0.2572 0.7970099
est_mtgs 0.051543 0.026364 1.9550 0.0505800
bls_unemp 0.182907 0.019182 9.5353 < 2.2e-16
est_90d_va 0.204964 0.022363 9.1652 < 2.2e-16
violent 0.259874 0.034796 7.4685 8.105e-14
property -0.134602 0.034828 -3.8648 0.0001112
Rho: 0.79874
Asymptotic standard error: 0.077982
z-value: 10.243, p-value: < 2.22e-16
Lambda: 0.62618
Asymptotic standard error: 0.16424
z-value: 3.8126, p-value: 0.00013751
LR test value: 448.42, p-value: < 2.22e-16
Log likelihood: -755.5611 for sac model
ML residual variance (sigma squared): 0.30991, (sigma: 0.5567)
Number of observations: 897
Number of parameters estimated: 9
AIC: 1529.1, (AIC for lm: 1973.5)10. Pemilihan Model Terbaik
AICc(foreclosure.lm)[1] 1922.515AICc(chipoly.sar)[1] 1542.882AICc(chipoly.err)[1] 1570.464AICc(chipoly.sarma)[1] 1529.325Berdasarkan nilai AICc, seluruh model spasial menunjukkan penurunan AICc yang sangat drastis dibandingkan model OLS. Hal ini mengonfirmasi bahwa OLS memang tidak layak digunakan karena mengabaikan struktur dependensi spasial dalam data. Di antara model spasial, SARMA unggul atas SAR maupun SEM karena mampu menangkap dua sumber dependensi spasial sekaligus, yakni pengaruh lag spasial pada variabel dependen (spatial lag) dan autokorelasi pada komponen error (spatial error).
11. Interpretasi Efek Marginal
impacts(chipoly.sarma, listw = W_dist_listw)Impact measures (sac, exact):
Direct Indirect Total
est_mtgs dy/dx 0.05245476 0.2036477 0.2561024
bls_unemp dy/dx 0.18614185 0.7226676 0.9088094
est_90d_va dy/dx 0.20858933 0.8098165 1.0184059
violent dy/dx 0.26447046 1.0267666 1.2912370
property dy/dx -0.13698310 -0.5318162 -0.6687993# Impact dari output impacts()
impacts_raw <- data.frame(
variabel = c("est_mtgs",
"bls_unemp",
"est_90d_va",
"violent",
"property"),
direct = c(0.05245476,
0.18614185,
0.20858933,
0.26447046,
-0.13698310),
indirect = c(0.2036477,
0.7226676,
0.8098165,
1.0267666,
-0.5318162),
total = c(0.2561024,
0.9088094,
1.0184059,
1.2912370,
-0.6687993)
)
# Statistik asli manual
stats_asli <- data.frame(
variabel = c("est_mtgs",
"bls_unemp",
"est_90d_va",
"violent",
"property",
"est_fcs"),
mean = c(
mean(chi_df$est_mtgs, na.rm = TRUE),
mean(chi_df$bls_unemp, na.rm = TRUE),
mean(chi_df$est_90d_va, na.rm = TRUE),
mean(chi_df$violent, na.rm = TRUE),
mean(chi_df$property, na.rm = TRUE),
mean(chi_df$est_fcs, na.rm = TRUE)
),
sd = c(
sd(chi_df$est_mtgs, na.rm = TRUE),
sd(chi_df$bls_unemp, na.rm = TRUE),
sd(chi_df$est_90d_va, na.rm = TRUE),
sd(chi_df$violent, na.rm = TRUE),
sd(chi_df$property, na.rm = TRUE),
sd(chi_df$est_fcs, na.rm = TRUE)
)
)
# SD variabel dependen
sd_Y <- stats_asli$sd[
stats_asli$variabel == "est_fcs"
]
# Konversi ke satuan asli
impacts_asli <- impacts_raw %>%
left_join(stats_asli, by = "variabel") %>%
mutate(
direct_asli = direct * (sd_Y / sd),
indirect_asli = indirect * (sd_Y / sd),
total_asli = total * (sd_Y / sd)
) %>%
select(
variabel,
mean,
sd,
direct_asli,
indirect_asli,
total_asli
)
print(impacts_asli) variabel mean sd direct_asli indirect_asli total_asli
1 est_mtgs 660.619844 623.9357845 0.003532374 0.01371391 0.01724628
2 bls_unemp 8.015942 0.9352458 8.362574000 32.46642968 40.82900143
3 est_90d_va 6.766232 5.9369915 1.476207329 5.73115150 7.20735933
4 violent 169.338907 189.7759281 0.058554158 0.22732767 0.28588182
5 property 419.957822 344.1708470 -0.016723016 -0.06492458 -0.08164760Hasil impact menunjukkan adanya pengaruh langsung (direct), tidak langsung/spillover (indirect), dan total antarwilayah pada model spasial.
est_mtgsmemiliki pengaruh positif tetapi sangat kecil terhadap variabel dependen.bls_unempmemiliki pengaruh paling besar dan positif. Kenaikan tingkat pengangguran meningkatkan variabel dependen sebesar 8.36 pada wilayah sendiri dan 32.47 pada wilayah tetangga, dengan total pengaruh 40.83. Ini menunjukkan efek spasial yang sangat kuat.est_90d_vajuga berpengaruh positif cukup besar, dengan total effect sebesar 7.21. Artinya peningkatan variabel ini meningkatkan variabel dependen baik secara lokal maupun pada wilayah sekitar.violentberpengaruh positif namun relatif kecil. Peningkatan kriminalitas kekerasan meningkatkan variabel dependen dengan total effect sebesar 0.29.propertymemiliki pengaruh negatif. Kenaikan kriminalitas properti menurunkan variabel dependen baik secara langsung maupun melalui spillover, dengan total effect sebesar −0.08.
Secara umum, nilai indirect effect yang lebih besar daripada direct effect menunjukkan bahwa hubungan spasial antarwilayah cukup kuat, sehingga model spasial relevan digunakan.
Studi Kasus 2: Data Artikel China
1. Persiapan Data
suppressPackageStartupMessages({
library(sf)
library(readxl)
library(tidyverse)
library(ggplot2)
library(spdep)
library(spatialreg)
library(MuMIn)
library(patchwork)
})shp_china <- st_read("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 12/china/China29/China29.shp")Reading layer `China29' from data source
`/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 12/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 84covid_china <- read_excel("/Users/rositariarusesta/Desktop/Tugas Spasial Pertemuan 12/china/COVID19_China.xlsx")
china_join <- shp_china %>%
left_join(covid_china,
by = c("NAME_1" = "Region"))
head(china_join)Simple feature collection with 6 features and 17 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 Rc Rr Rd
1 ??|?? ?nhu? 117.2262 31.82579 0.1122 0.0610 6e-04
2 ??|?? B?ij?ng 116.4107 40.18491 0.0481 0.0201 6e-04
3 ??|?? Chóngqìng 107.8748 30.05865 0.0683 0.0342 6e-04
4 ?? Fújiàn 117.9815 26.08641 0.0356 0.0177 1e-04
5 ??|?? G?nsù 100.9339 37.82079 0.0115 0.0069 2e-04
6 ??|?? Gu?ngd?ng 113.4164 23.34225 0.1586 0.0784 5e-04
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 dan Pra-pemrosesan data
Rc <- ggplot(data = china_join) +
geom_sf(
aes(fill = Rc),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(
option = "magma",
name = "Rc"
) +
labs(
title = "Rc di China",
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)
)
Rr <- ggplot(data = china_join) +
geom_sf(
aes(fill = Rr),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(
option = "magma",
name = "Rr"
) +
labs(
title = "Rr di China",
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)
)
Rd <- ggplot(data = china_join) +
geom_sf(
aes(fill = Rd),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(
option = "magma",
name = "Rd"
) +
labs(
title = "Rd di China",
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)
)Rc | Rr | Rd
3. Matriks Pembobot Spasial
Untuk menangkap dependensi spasial antar wilayah, analisis ini menggunakan pendekatan k-nearest neighbor (KNN)dalam membangun matriks bobot spasial. Pertama, titik sentroid dari setiap unit wilayah diekstraksi menggunakan fungsi st_centroid(), kemudian koordinatnya disimpan melalui st_coordinates() sebagai representasi lokasi geografis masing-masing observasi.
Berdasarkan koordinat sentroid tersebut, algoritma KNN diterapkan dengan nilai k = 6, artinya setiap wilayah didefinisikan memiliki tepat enam tetangga terdekat berdasarkan jarak Euclidean. Pendekatan ini dipilih karena bersifat adaptif terhadap keragaman kepadatan spasial — wilayah yang secara geografis terisolasi tetap memperoleh tetangga, sehingga tidak ada unit observasi yang terisolir dari struktur ketetanggaan (no islands problem).
Objek ketetanggaan kemudian dikonversi ke dalam bentuk spatial weights matrix menggunakan skema row-standardization (style = "W"). Dalam skema ini, bobot setiap tetangga dinormalisasi sehingga jumlah bobot untuk setiap wilayah bernilai satu. Hal ini memastikan bahwa rata-rata tertimbang spasial (spatially lagged variable) dapat diinterpretasikan secara konsisten tanpa dipengaruhi oleh perbedaan jumlah tetangga antar wilayah.
Matriks bobot spasial W yang dihasilkan selanjutnya digunakan sebagai dasar dalam uji autokorelasi spasial maupun pemodelan regresi spasial pada tahapan analisis berikutnya.
coords <- st_coordinates(st_centroid(china_join))Warning: st_centroid assumes attributes are constant over geometriesknn <- knearneigh(coords, k = 6)
nb_knn <- knn2nb(knn)
W <- nb2listw(nb_knn, style = "W")4. OLS Regression
covid.lm<-lm(Rc~Rd+Rr, data=china_join)
summary(covid.lm)
Call:
lm(formula = Rc ~ Rd + Rr, data = china_join)
Residuals:
Min 1Q Median 3Q Max
-0.0227795 -0.0021997 -0.0000442 0.0028980 0.0142696
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.415e-04 1.782e-03 0.192 0.849
Rd 1.142e+01 4.536e-01 25.180 <2e-16 ***
Rr 1.764e+00 5.183e-02 34.027 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.006583 on 26 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 4.635e+05 on 2 and 26 DF, p-value: < 2.2e-165. Uji Morans’I
lm.morantest(covid.lm,W)
Global Moran I for regression residuals
data:
model: lm(formula = Rc ~ Rd + Rr, data = china_join)
weights: W
Moran I statistic standard deviate = -0.16413, p-value = 0.5652
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
-0.061455011 -0.048414285 0.006313236 Uji Global Moran I pada residual model regresi OLS (Rc ~ Rd + Rr) menghasilkan nilai Moran I sebesar -0.0615 dengan p-value = 0.5652. Karena p-value jauh di atas ambang signifikansi α = 0.05, maka H₀ tidak dapat ditolak — artinya tidak terdapat autokorelasi spasial yang signifikan pada residual model.
Nilai Moran I yang sangat kecil dan mendekati nilai ekspektasi teoritis (-0.0484) menunjukkan bahwa residual regresi terdistribusi secara acak di ruang geografis, tanpa pola pengelompokan (clustering) maupun dispersi spasial yang sistematis.
6. Uji Lagrange Multiplier
LM<-lm.LMtests(covid.lm, W, test="all")Please update scripts to use lm.RStests in place of lm.LMtestssummary(LM) Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
dependence
data:
model: lm(formula = Rc ~ Rd + Rr, data = china_join)
test weights: listw
statistic parameter p.value
RSerr 0.37613 1 0.5397
RSlag 0.72932 1 0.3931
adjRSerr 0.36732 1 0.5445
adjRSlag 0.72051 1 0.3960
SARMA 1.09664 2 0.5779Seluruh statistik uji LM menghasilkan p-value jauh di atas α = 0.05, sehingga H₀ tidak dapat ditolak pada semua spesifikasi uji.
RSerr & adjRSerr — tidak terdapat dependensi spasial pada komponen error, artinya tidak ada indikasi bahwa gangguan model menyebar secara sistematis mengikuti struktur ketetanggaan spasial. Model Spatial Error (SEM) tidak diperlukan.
RSlag & adjRSlag — tidak terdapat dependensi spasial pada variabel dependen (
Rc), artinya nilai konsumsi rokok di suatu provinsi tidak secara signifikan dipengaruhi oleh nilai konsumsi rokok di provinsi-provinsi tetangganya. Model Spatial Lag (SLM) tidak diperlukan.SARMA — kombinasi keduanya pun tidak signifikan, mengonfirmasi bahwa tidak ada struktur autokorelasi spasial baik pada lag maupun error secara bersamaan.
Hasil ini konsisten dengan temuan Global Moran I pada residual yang juga tidak signifikan. Secara diagnostik, model OLS sudah memadai — namun demikian, ketiadaan dependensi spasial pada Y dan error tidak berarti bahwa dependensi spasial pada variabel X (spillover effect) juga absen. Oleh karena itu, estimasi model SLX tetap relevan dan perlu dilakukan
7. Model SLX
Uji Lagrange Multiplier (LM) yang umum digunakan dalam konteks pemilihan model spasial — yakni LM-Lag dan LM-Error — dirancang secara spesifik untuk mendeteksi dependensi spasial pada variabel dependen (Y) melalui mekanisme spatial lag, serta dependensi spasial pada komponen error, bukan pada variabel independen (X). Dengan demikian, uji LM tidak mengevaluasi apakah variabel penjelas seperti Rd (rasio ketergantungan) dan Rr (rasio ruralitas) itu sendiri mengandung struktur autokorelasi spasial.
Justru di sinilah letak urgensi penggunaan Spatially Lagged X (SLX) model. Berbeda dengan SAR maupun SEM, model SLX secara eksplisit mengakomodasi kemungkinan bahwa nilai variabel independen di wilayah-wilayah tetangga turut memengaruhi variabel dependen di wilayah yang bersangkutan — suatu fenomena yang dikenal sebagai spillover effect atau efek limpahan spasial.
china.SLX <- lmSLX(Rc ~ Rd+Rr , 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) -3.059e-04 2.808e-03 -1.089e-01 9.142e-01
Rd 1.131e+01 6.068e-01 1.864e+01 8.808e-16
Rr 1.776e+00 6.970e-02 2.548e+01 6.914e-19
lag.Rd -3.893e-01 9.052e-01 -4.301e-01 6.710e-01
lag.Rr 3.709e-02 1.040e-01 3.565e-01 7.246e-01Efek Lokal (Direct Effect)
Rd — rasio ketergantungan di suatu provinsi berpengaruh positif dan sangat signifikan terhadap tingkat konsumsi rokok (
Rc). Setiap kenaikan satu satuanRddi provinsi tersebut, konsumsi rokok meningkat sebesar 11.31 satuan, ceteris paribus.Rr — rasio ruralitas juga berpengaruh positif dan sangat signifikan. Setiap kenaikan satu satuan
Rr, konsumsi rokok meningkat sebesar 1.776 satuan, mengindikasikan bahwa wilayah dengan tingkat ruralitas lebih tinggi cenderung memiliki konsumsi rokok yang lebih besar.
Efek Limpahan Spasial (Spatial Spillover Effect)
lag.Rd — rasio ketergantungan di provinsi-provinsi tetangga tidak berpengaruh signifikan (p = 0.671) terhadap konsumsi rokok di provinsi yang bersangkutan. Tidak terdapat bukti adanya efek limpahan spasial dari variabel ini.
lag.Rr — demikian pula, rasio ruralitas di wilayah tetangga tidak berpengaruh signifikan (p = 0.725) terhadap konsumsi rokok lokal. Efek spillover dari ruralitas tetangga tidak terbukti secara statistik.
AICc(covid.lm)[1] -202.5489AICc(china.SLX)[1] -197.3188Perbandingan nilai AICc (Akaike Information Criterion corrected) digunakan untuk mengevaluasi keseimbangan antara goodness-of-fit dan kompleksitas model. Semakin kecil (lebih negatif) nilai AICc, semakin baik model tersebut.
Hasil menunjukkan bahwa model OLS memiliki AICc sebesar -202.55, sedangkan model SLX menghasilkan AICc sebesar -197.32. Selisih keduanya adalah ΔAICc = 5.23, di mana model OLS memiliki nilai yang lebih kecil (lebih baik).
Berdasarkan kaidah umum interpretasi AICc, selisih lebih dari 2 poin sudah cukup menjadi bukti bahwa model dengan AICc lebih rendah lebih didukung oleh data. Dengan demikian, model OLS lebih disukai dibandingkan model SLX secara kriteria informasi.