1 Matriks Bobot Spasial

1.1 Definisi

Matriks bobot spasial W merupakan komponen penting dalam analisis spasial. Matriks ini berukuran n × n, bersifat simetris, non-stokastik, dan nilainya positif. Elemen \(w_{ij}\) menunjukkan kedekatan spasial antara lokasi i dan j.

Pembentukan matriks pembobot didasari oleh Hukum Pertama Geografi dari Waldo Tobler.

1.2 Jenis

Terdapat dua pendekatan umum dalam pembentukan matriks bobot spasial:

1.2.1 Pendekatan Persinggungan (Contiguity)

Pendekatan ini digunakan untuk data spasial berbasis poligon atau area. Dua lokasi \(i\) dan \(j\) memiliki hubungan spasial jika keduanya bersinggungan secara geografis. Matriks bobot spasialnya didefinisikan sebagai:

\[ w_{ij} = \begin{cases} 1 & \text{jika } i \text{ dan } j \text{ bersinggungan} \\ 0 & \text{lainnya} \end{cases} \]

Tiga jenis matriks bobot berdasarkan persinggungan:

  • Bishop Contiguity: Persinggungan sudut.
  • Rook Contiguity: Persinggungan sisi (atas, bawah, kiri, kanan).
  • Queen Contiguity: Persinggungan sisi dan sudut.

Misalkan lokasi amatan adalah lokasi E, maka bobot spasial digambarkan seperti berikut:

Sumber: Oyana, T.J. (2020). Spatial Analysis with R: Statistics, Visualization, and Computational Methods (2nd ed.). CRC Press. https://doi.org/10.1201/9781003021643

Penjelasan:

  • Bishop: Lokasi E bersinggungan dengan A, C, G, dan I.
  • Rook: Lokasi E bersinggungan dengan B, D, F, dan H.
  • Queen: Lokasi E bersinggungan dengan A, B, C, D, F, G, H, dan I.

1.2.2 Pendekatan Jarak (Distance)

  • K-Nearest Neighbors

    Jumlah tetangga terdekat \(k\) didefinisikan sebagai:

\[ w_{ij} = \begin{cases} 1 & \text{jika titik pusat lokasi } j \text{ adalah } k \text{ tetangga terdekat dari } i \\ 0 & \text{lainnya} \end{cases} \]

  • Jarak Ambang Batas :

    Bobot 1 diberikan jika jarak dua lokasi \(\leq d_{\text{max}}\) atau dapat ditulis sebagai berikut :

\[ w_{ij} = \begin{cases} 1 & \text{jika } 0 \leq d_{ij} \leq d_{\text{max}} \\ 0 & \text{lainnya} \end{cases} \]

Jika lokasi hanya memiliki satu tetangga, maka bobot otomatis 1.

Ukuran jarak yang sering digunakan adalah jarak Euclidean:

\[ d_{ij} = \sqrt{(u_i - u_j)^2 + (v_i - v_j)^2} \]

  • Inverse Distance :

    Jarak dalam pendekatan ini didefinisikan sebagai berikut:

\[ w_{ij} = \begin{cases} \frac{1}{d_{ij}^{\alpha}} & \text{jika } i \neq j \\ 0 & \text{jika } i = j \end{cases} \] Umumnya \(\alpha = \{1, 2\}\).

1.3 Visualisasi

Peta Kota Bandung

Indo_Kec<-readRDS('gadm36_IDN_3_sp.rds')
Bandung<-Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
plot(Bandung)

Matriks Bobot K-Nearest Neigbhors

CoordK<-coordinates(Bandung)
W <- knn2nb(knearneigh(CoordK, k = 5), row.names =rownames(Bandung))
W
## Neighbour list object:
## Number of regions: 30 
## Number of nonzero links: 150 
## Percentage nonzero weights: 16.66667 
## Average number of links: 5 
## Non-symmetric neighbours list
WB <- nb2mat(W, style='B', zero.policy
= TRUE) #menyajikan dalam bentuk matrix biner "B" 
WB[1:5,1:5]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    0    0    1    0
## 2    0    0    1    0    0
## 3    0    1    0    0    0
## 4    0    0    0    0    0
## 5    1    0    0    1    0
library(sp)
library(spdep)
plot(Bandung, border = "gray", main = "K-Nearest Neigbhors K=5")
plot(W, CoordK, add = TRUE, col = "blue", lwd = 1)
points(CoordK, pch = 20, col = "red")
text(CoordK, labels = 1:nrow(Bandung), cex = 0.6, pos = 3)

Matriks Bobot Inverse Jarak

Dist.mat<-as.matrix(dist(CoordK, method="euclidean"))
Dist.mat[1:5,1:5]
##            6199       6210       6221       6223       6224
## 6199 0.00000000 0.08105203 0.09725312 0.03125097 0.03277083
## 6210 0.08105203 0.00000000 0.01633011 0.06183450 0.08695008
## 6221 0.09725312 0.01633011 0.00000000 0.07664345 0.10141706
## 6223 0.03125097 0.06183450 0.07664345 0.00000000 0.02525779
## 6224 0.03277083 0.08695008 0.10141706 0.02525779 0.00000000
Dist.mat.inv<-1/Dist.mat
diag(Dist.mat.inv)<-0
Dist.mat.inv[1:5,1:5]
##          6199     6210      6221     6223      6224
## 6199  0.00000 12.33775 10.282447 31.99900 30.514940
## 6210 12.33775  0.00000 61.236564 16.17220 11.500852
## 6221 10.28245 61.23656  0.000000 13.04743  9.860274
## 6223 31.99900 16.17220 13.047429  0.00000 39.591741
## 6224 30.51494 11.50085  9.860274 39.59174  0.000000
Dist.mat.invs<-mat2listw(Dist.mat.inv, style="W")
summary(Dist.mat.invs)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 30 
## Number of nonzero links: 870 
## Percentage nonzero weights: 96.66667 
## Average number of links: 29 
## Link number distribution:
## 
## 29 
## 30 
## 30 least connected regions:
## 6199 6210 6221 6223 6224 6225 6226 6227 6228 6200 6201 6202 6203 6204 6205 6206 6207 6208 6209 6211 6212 6213 6214 6215 6216 6217 6218 6219 6220 6222 with 29 links
## 30 most connected regions:
## 6199 6210 6221 6223 6224 6225 6226 6227 6228 6200 6201 6202 6203 6204 6205 6206 6207 6208 6209 6211 6212 6213 6214 6215 6216 6217 6218 6219 6220 6222 with 29 links
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 30 900 30 2.901426 120.6623

Matriks Bobot Queen

W <- poly2nb(Bandung, row.names=rownames(Bandung), queen=TRUE) #Mendapatkan W
WB <- nb2mat(W, style='B', zero.policy = 
TRUE) #menyajikan dalam bentuk matrixbiner "B" 
WB[1:5,1:5]
##      [,1] [,2] [,3] [,4] [,5]
## 6199    0    0    0    1    1
## 6210    0    0    1    0    0
## 6221    0    1    0    0    0
## 6223    1    0    0    0    0
## 6224    1    0    0    0    0
WL<-nb2listw(W)
WL
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 30 
## Number of nonzero links: 140 
## Percentage nonzero weights: 15.55556 
## Average number of links: 4.666667 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 30 900 30 13.81692 123.3523
library(spdep)
library(sp)
plot(Bandung, border = "gray",main = "Queen")
coords <- coordinates(Bandung)
plot(W, coords, add = TRUE, col = "blue", lwd = 1)
points(coords, pch = 20, col = "red")
text(coords, labels = 1:nrow(Bandung), cex = 0.6, pos = 3)

2 Model Spasial Ekonometrika

Spasial ekonometrika merupakan cabang dari ilmu ekonometrika yang secara khusus mempelajari unit-unit geografis (seperti kota, wilayah, atau negara) beserta hubungan dan interaksi di antara unit-unit tersebut.

Dalam data spasial, sering ditemukan dua karakteristik penting: - Dependensi spasial, yaitu ketika suatu lokasi dipengaruhi oleh kondisi di lokasi sekitarnya. - Heterogenitas spasial, yaitu ketika hubungan antar variabel berbeda-beda antar wilayah.

Kedua efek tersebut dapat membuat teknik ekonometrika standar menjadi tidak tepat digunakan karena dapat menghasilkan estimasi yang bias, tidak konsisten, atau tidak efisien. Oleh karena itu, dibutuhkan pendekatan khusus dalam spasial ekonometrika untuk mengakomodasi struktur data yang bersifat spasial.

Regresi spasial merupakan salah satu bidang yang berkembang pesat dalam literatur ekonometrika dan statistika. Regresi spasial merupakan pengembangan dari regresi linier klasik yang dirancang untuk menangani efek spasial dalam data. Model ini secara khusus mempertimbangkan pengaruh lokasi melalui mekanisme spasial, seperti ketergantungan antar unit wilayah maupun variasi hubungan antar wilayah.

Menurut LeSage (1999), bentuk umum dari model regresi spasial dapat dituliskan sebagai berikut:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \] \[ \mathbf{u} = \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim N(0, \sigma_\varepsilon^2 \mathbf{I}_n) \]

Dengan:

  • \(\mathbf{y}\): vektor variabel respon berukuran \(n \times 1\)
  • \(\rho\): koefisien parameter spasial lag dari variabel respon
  • \(\mathbf{W}\): matriks pembobot spasial \(n \times n\)
  • \(\mathbf{X}\): matriks variabel prediktor \(n \times (p+1)\)
  • \(\boldsymbol{\beta}\): vektor parameter regresi \((p+1) \times 1\)
  • \(\lambda\): koefisien parameter spasial error
  • \(\mathbf{u}\): vektor error yang memiliki efek spasial berukuran \(n \times 1\)
  • \(\boldsymbol{\varepsilon}\): vektor error berukuran \(n \times 1\)

Struktur dari \(\mathbf{y}\), \(\mathbf{X}\), \(\boldsymbol{\beta}\), dan \(\boldsymbol{\varepsilon}\) adalah:

\[ \mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}, \quad \mathbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{pmatrix}, \quad \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{pmatrix} \]

Menurut Anselin (1988), dari model regresi spasial umum dapat diturunkan tiga jenis model berikut:

  1. Regresi Linier Klasik (jika \(\rho = 0\) dan \(\lambda = 0\)):

    \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

  2. Spatial Autoregressive Model (SAR) (jika \(\rho \neq 0\) dan \(\lambda = 0\)):

    \[ \mathbf{y} = \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

    Spatial Autoregressive Model (SAR), atau disebut juga Spatial Lag Model (SLM), merupakan model spasial berbasis area yang memperhitungkan pengaruh spasial lag pada variabel dependen. Model ini juga disebut Mixed Regressive–Autoregressive karena menggabungkan regresi biasa dengan regresi spasial lag pada variabel dependen.

  3. Spatial Error Model (SEM) (jika \(\rho = 0\) dan \(\lambda \neq 0\)):

    \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u}, \quad \mathbf{u} = \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon} \]

    Spatial Error Model (SEM) digunakan ketika terdapat korelasi spasial antar error, yaitu ketika nilai error pada suatu lokasi berkorelasi dengan nilai error di lokasi sekitarnya. Dengan kata lain, error pada lokasi \(i\) merupakan fungsi dari error pada lokasi \(j\) di sekitarnya.

Sumber:
Yasin, H., Warsito, B., & Hakim, A. R. (2020). Regresi Spasial (Aplikasi dengan R). WADE Publish.

3 Pengujian Hipotesis

3.1 Uji Dependensi Spasial

Salah satu cara untuk mengetahui adanya dependensi spasial antar lokasi adalah dengan melakukan uji autokorelasi spasial menggunakan statistik Moran’s I. Statistik ini mengukur taksiran dari korelasi nilai amatan yang berasal dari lokasi-lokasi pada variabel yang sama.

Hipotesis:

  • H₀: Tidak ada autokorelasi spasial antar lokasi
  • H₁: Ada autokorelasi spasial antar lokasi

Statistik Uji (Z):

\[ Z_{\text{Hitung}} = \frac{I - E(I)}{\sqrt{\text{Var}(I)}} \]

Dengan:

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

\[ E(I) = -\frac{1}{n-1}, \quad \text{Var}(I) = \frac{n^2 S_1 - n S_2 + 3 S_2^2}{(n-1)(n-2)(n-3)S_0^2} - [E(I)]^2 \]

Keterangan:

  • \(c_{ij} = (x_i - \bar{x})(x_j - \bar{x})\)
  • \(w_{ij}\): bobot spasial antar wilayah
  • \(S_0 = \sum_i \sum_j w_{ij}\)
  • \(S_1 = \frac{1}{2} \sum_i \sum_j (w_{ij} + w_{ji})^2\)
  • \(S_2 = \sum_i (w_{i+} + w_{+i})^2\)

Keputusan:
Tolak H₀ jika \(|Z_{\text{Hitung}}| > Z_{\alpha/2}\)

3.2 LM Test

Lagrange Multiplier (LM) Test merupakan uji statistik asimtotik yang sering digunakan karena berbasis pada estimasi model di bawah hipotesis nol (null hypothesis), yaitu bentuk model yang paling sederhana. Dalam banyak situasi, pendekatan ini menghindari prosedur estimasi nonlinier dan memungkinkan penggunaan metode Ordinary Least Squares (OLS).

Untuk menguji adanya ketergantungan spasial, uji LM menggunakan hipotesis spesifik dan menjadi dasar dalam pemilihan model regresi spasial yang sesuai. Dua jenis uji LM yang telah dikembangkan antara lain:

  • LM-LAG: untuk menguji ketergantungan spasial pada variabel dependen.
  • LM-ERR: untuk menguji ketergantungan spasial pada error.

Pada uji LM-Lag, model regresinya dinyatakan sebagai berikut:

\[ \mathbf{y} = \rho \mathbf{W}_\rho \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

Hipotesis:

  • \(H_0\): \(\rho = 0\), artinya tidak ada spatial lag dan model regresi standar:

    \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

    adalah valid.

  • \(H_1\): \(\rho \neq 0\)

Statistik uji LM-Lag dirumuskan sebagai berikut:

\[ \text{LM}_{\text{LAG}} = \frac{[(\mathbf{e}^T \mathbf{W}_A \mathbf{y}) / (\mathbf{e}^T \mathbf{e} / n)]^2}{[(\mathbf{W}_A \mathbf{X} \hat{\boldsymbol{\beta}})^T \mathbf{M} (\mathbf{W}_A \mathbf{X} \hat{\boldsymbol{\beta}}) / (\mathbf{e}^T \mathbf{e} / n)] + [\text{tr}(\mathbf{W}_A^T \mathbf{W}_A + \mathbf{W}_A^2)]} \sim \chi^2_{(1 - \alpha); df = 1} \]

Dengan: - \(\mathbf{e}\) : vektor residual OLS
- \(\mathbf{W}_A\) : matriks bobot spasial
- \(\mathbf{M} = (\mathbf{I}_n - \mathbf{X}(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T)\)

Pada uji LM-Error test, didasarkan pada estimasi model regresi:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

LM-ERR memperhitungkan dua tipe ketergantungan error spasial, yaitu:

  • Ketergantungan AR spasial (SpAR):

\[ \boldsymbol{\varepsilon} = \lambda \mathbf{W}_\lambda \boldsymbol{\varepsilon} + \boldsymbol{\zeta} \]

  • Ketergantungan MA spasial (SpMA):

\[ \boldsymbol{\varepsilon} = \boldsymbol{\zeta} + \theta \mathbf{W}_\theta \boldsymbol{\zeta} \]

3.3 Pengujian dalam Model SAR

Perbandingan dengan Model OLS

Hipotesis:

  • H₀: \(\rho = 0\) (Model OLS lebih baik)
  • H₁: \(\rho \neq 0\) (Model SAR lebih baik)

Statistik Uji:

\[ F_{\text{hitung}} = \frac{(SSE_C - SSE_U) / ((p+1) - p)}{SSE_U / (n - p)} \sim F_{\alpha/2; p; (n-p)} \]

Keterangan:

  • \(SSE_C\): SSE dari model constrained (OLS)
  • \(SSE_U\): SSE dari model unconstrained (SAR)

Uji Kecocokan Model Keseluruhan

Hipotesis:

  • H₀: \(\theta_i = 0, \forall i\) (Model SAR tidak cocok)
  • H₁: \(\exists \theta_i \neq 0\) (Model SAR cocok)

Statistik Uji:

\[ F_{\text{hitung}} = \frac{(SSE_C - SSE_U) / (p+1)}{SSE_U / (n - p)} \]

Uji Signifikansi Parameter Secara Individu (Uji Wald)

Uji ini digunakan untuk menguji signifikansi koefisien dalam model secara individu.

Hipotesis:

  • H₀: \(\beta_j = 0\) (parameter tidak signifikan)
  • H₁: \(\beta_j \neq 0\) (parameter signifikan)

Statistik Uji Wald (W):

\[ W = \left( \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \right)^2, \quad SE(\hat{\beta}_j) = \sqrt{\sigma^2(\hat{\beta}_j)} \]

Keputusan:
Tolak H₀ jika \(W > \chi^2_{\alpha, 1}\)

3.4 Pengujian dalam Model SEM

Perbandingan dengan Model Regresi OLS

Hipotesis:

\[ H_0 : \lambda = 0 \quad \text{(Model OLS lebih baik)} \\ H_1 : \lambda \neq 0 \quad \text{(Model SEM lebih baik)} \]

Statistik Uji:

\[ F_{\text{hitung}} = \frac{(SSE_C - SSE_U) / ((p+1) - p)}{SSE_U / (n - p)} \sim F_{\alpha/2; 1, (n-p)} \]

Keterangan:

  • \(SSE_C\): SSE dari model constrained (dengan asumsi \(H_0\) benar, yaitu model OLS)
  • \(SSE_U\): SSE dari model unconstrained (dengan asumsi \(H_0\) salah, yaitu model SEM)

Uji Signifikansi Parameter

Pengujian signifikansi parameter dilakukan menggunakan uji Wald (Anselin, 1988).

Hipotesis untuk Parameter \(\lambda\):

\[ H_0 : \lambda = 0 \\ H_1 : \lambda \neq 0 \]

Statistik Uji:

\[ Wald_{\lambda} = \frac{\hat{\lambda}^2}{\operatorname{var}(\hat{\lambda})} \]

Hipotesis untuk Parameter \(\beta\):

\[ H_0 : \beta_j = 0 \\ H_1 : \beta_j \neq 0, \quad j = 1, 2, \ldots, p \]

Statistik Uji:

\[ Wald_{\beta} = \frac{\hat{\beta}_j^2}{\operatorname{var}(\hat{\beta}_j)} \]

Keterangan:

  • \(\operatorname{var}(\hat{\lambda})\): elemen diagonal dari matriks varians untuk \(\lambda\)
  • \(\operatorname{var}(\hat{\beta}_j)\): elemen diagonal dari matriks varians untuk \(\beta_j\)

Kriteria Pengambilan Keputusan:
Tolak \(H_0\) jika:

\[ Wald > \chi^2_{\alpha;1} \]

4 Contoh Perhitungan SAR dan SEM dengan Data Simulasi

#30 Kecamatan di Kota Bandung
set.seed(123)
Indo_Kec <- readRDS('gadm36_IDN_3_sp.rds')
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
n <- nrow(Bandung)

4.1 Matriks Bobot

## Matriks Bobot Queen
coords <- coordinates(Bandung)
bandung_nb <- poly2nb(Bandung, queen=TRUE)
bandung_listw <- nb2listw(bandung_nb, style="W", zero.policy=TRUE)
WB <- nb2mat(bandung_nb, style='B', zero.policy=TRUE)

4.2 Data Simulasi

set.seed(123)
nama_kecamatan <- Bandung$NAME_3

# 1. Kepadatan Penduduk (jiwa/km²)
KEPADATAN <- runif(n, 5000, 30000)

# 2. Pendapatan Rata-rata Rumah Tangga per Bulan (Rp)
PENDAPATAN <- runif(n, 2500000, 8000000)

# 3. Persentase Luas Ruang Terbuka Hijau (%)
RTH <- runif(n, 5, 30)

# 4. Jumlah Fasilitas Pendidikan (unit)
PENDIDIKAN <- round(runif(n, 10, 50))

beta0 <- 50
beta1 <- -0.0003  # koefisien KEPADATAN (negatif, karena semakin padat, IKM cenderung lebih rendah)
beta2 <- 0.000005  # koefisien PENDAPATAN
beta3 <- 0.5      # koefisien RTH
beta4 <- 0.2      # koefisien PENDIDIKAN

epsilon <- rnorm(n, 0, 5)

X <- cbind(1, KEPADATAN, PENDAPATAN, RTH, PENDIDIKAN)
betas <- c(beta0, beta1, beta2, beta3, beta4)

rho <- 0.6
W_mat <- nb2mat(bandung_nb, style="W", zero.policy=TRUE)
I_mat <- diag(n)
IKM <- solve(I_mat - rho * W_mat) %*% (X %*% betas + epsilon)
IKM <- as.numeric(IKM)
IKM <- 40 + (IKM - min(IKM)) / (max(IKM) - min(IKM)) * 40 

data_simulasi <- data.frame(
  Kecamatan = nama_kecamatan,
  IKM = IKM,
  KEPADATAN = KEPADATAN,
  PENDAPATAN = PENDAPATAN,
  RTH = RTH,
  PENDIDIKAN = PENDIDIKAN
)
Bandung@data <- cbind(Bandung@data, data_simulasi[,-1])
head(data_simulasi, 5)
##         Kecamatan      IKM KEPADATAN PENDAPATAN       RTH PENDIDIKAN
## 1           Andir 79.27632  12189.44    7796633 21.627880         15
## 2        Antapani 57.98455  24707.63    7462645  7.371017         36
## 3       Arcamanik 59.84557  15224.42    6298879 14.599241         24
## 4     Astanaanyar 61.43071  27075.44    6875071 11.859591         36
## 5 Babakan Ciparay 43.76842  28511.68    2635375 25.366001         23
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
tmap_mode("plot")
## tmap mode set to plotting
# Plot IKM di peta Kota Bandung
tm_shape(Bandung) +
  tm_fill("IKM", 
          title = "Indeks Kesejahteraan Masyarakat",
          style = "quantile", 
          palette = "YlOrRd", 
          legend.hist = TRUE) +
  tm_borders() +
  tm_layout(main.title = "Peta Simulasi IKM Kota Bandung",
            legend.outside = TRUE)

4.3 Indeks Moran

# Uji Moran Global untuk IKM
moran_test <- moran.test(Bandung$IKM, bandung_listw)
print(moran_test)
## 
##  Moran I test under randomisation
## 
## data:  Bandung$IKM  
## weights: bandung_listw    
## 
## Moran I statistic standard deviate = 1.2302, p-value = 0.1093
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.10445805       -0.03448276        0.01275509
# Plot Moran Scatterplot
moran_plot <- moran.plot(Bandung$IKM, bandung_listw, 
                         labels=as.character(1:n), pch=16, col="blue",
                         xlab="Indeks Kesejahteraan Masyarakat", 
                         ylab="Spatial Lag IKM")

4.4 OLS

# Formula regresi
reg_formula <- IKM ~ KEPADATAN + PENDAPATAN + RTH + PENDIDIKAN

# Estimasi model OLS
ols_model <- lm(reg_formula, data=Bandung@data)
summary(ols_model)
## 
## Call:
## lm(formula = reg_formula, data = Bandung@data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.0408  -3.6049  -0.3424   3.3898  10.8444 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.977e+01  6.949e+00   2.845  0.00874 ** 
## KEPADATAN   -3.125e-04  1.509e-04  -2.071  0.04879 *  
## PENDAPATAN   4.907e-06  7.023e-07   6.987 2.53e-07 ***
## RTH          7.268e-01  1.563e-01   4.650 9.24e-05 ***
## PENDIDIKAN   2.261e-01  9.896e-02   2.285  0.03105 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.808 on 25 degrees of freedom
## Multiple R-squared:  0.7405, Adjusted R-squared:  0.699 
## F-statistic: 17.84 on 4 and 25 DF,  p-value: 4.868e-07
# AIC OLS
AIC_OLS <- AIC(ols_model)

# Diagnostic tests for spatial dependence
lm_tests <- lm.RStests(ols_model, bandung_listw, test=c("all"))
print(lm_tests)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg_formula, data = Bandung@data)
## test weights: bandung_listw
## 
## RSerr = 0.99418, df = 1, p-value = 0.3187
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg_formula, data = Bandung@data)
## test weights: bandung_listw
## 
## RSlag = 6.3945, df = 1, p-value = 0.01145
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg_formula, data = Bandung@data)
## test weights: bandung_listw
## 
## adjRSerr = 2.1514, df = 1, p-value = 0.1424
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg_formula, data = Bandung@data)
## test weights: bandung_listw
## 
## adjRSlag = 7.5517, df = 1, p-value = 0.005995
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = reg_formula, data = Bandung@data)
## test weights: bandung_listw
## 
## SARMA = 8.5459, df = 2, p-value = 0.01394

4.5 SAR

# Estimasi model SAR
sar_model <- lagsarlm(reg_formula, data=Bandung@data, bandung_listw)
## Warning in lagsarlm(reg_formula, data = Bandung@data, bandung_listw): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   reciprocal condition number = 2.17163e-16 - using numerical Hessian.
summary(sar_model)
## 
## Call:
## lagsarlm(formula = reg_formula, data = Bandung@data, listw = bandung_listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -13.66082  -1.84810  -0.05338   3.17892  10.84705 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -4.4356e+00  1.0603e+01 -0.4183  0.675703
## KEPADATAN   -3.3424e-04  1.2210e-04 -2.7374  0.006192
## PENDAPATAN   5.0435e-06  5.6939e-07  8.8577 < 2.2e-16
## RTH          6.5423e-01  1.2908e-01  5.0685 4.009e-07
## PENDIDIKAN   2.3652e-01  8.0006e-02  2.9562  0.003114
## 
## Rho: 0.43181, LR test value: 5.9404, p-value: 0.014797
## Approximate (numerical Hessian) standard error: 0.16056
##     z-value: 2.6894, p-value: 0.0071574
## Wald statistic: 7.233, p-value: 0.0071574
## 
## Log likelihood: -89.64262 for lag model
## ML residual variance (sigma squared): 22.005, (sigma: 4.691)
## Number of observations: 30 
## Number of parameters estimated: 7 
## AIC: 193.29, (AIC for lm: 197.23)
impacts(sar_model, listw=bandung_listw)
## Impact measures (lag, exact):
##                   Direct      Indirect         Total
## KEPADATAN  -3.516884e-04 -2.365661e-04 -5.882546e-04
## PENDAPATAN  5.306784e-06  3.569652e-06  8.876436e-06
## RTH         6.883881e-01  4.630500e-01  1.151438e+00
## PENDIDIKAN  2.488630e-01  1.673998e-01  4.162628e-01
# AIC untuk model SAR
AIC_SAR <- AIC(sar_model)

4.6 SEM

# Estimasi model SEM
sem_model <- errorsarlm(reg_formula, data=Bandung@data, bandung_listw)
summary(sem_model)
## 
## Call:
## errorsarlm(formula = reg_formula, data = Bandung@data, listw = bandung_listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -13.43640  -3.06471   0.18825   3.85941  10.90396 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  2.2023e+01  5.9608e+00  3.6947 0.0002201
## KEPADATAN   -2.8066e-04  1.2397e-04 -2.2640 0.0235751
## PENDAPATAN   4.9174e-06  6.0394e-07  8.1423 4.441e-16
## RTH          5.5662e-01  1.3626e-01  4.0848 4.411e-05
## PENDIDIKAN   2.2936e-01  8.1524e-02  2.8134 0.0049021
## 
## Lambda: 0.4115, LR test value: 1.8848, p-value: 0.16979
## Asymptotic standard error: 0.20824
##     z-value: 1.9761, p-value: 0.048143
## Wald statistic: 3.905, p-value: 0.048143
## 
## Log likelihood: -91.67042 for error model
## ML residual variance (sigma squared): 25.31, (sigma: 5.0309)
## Number of observations: 30 
## Number of parameters estimated: 7 
## AIC: 197.34, (AIC for lm: 197.23)
# AIC untuk model SEM
AIC_SEM <- AIC(sem_model)

4.7 MODEL OLS, SAR, dan SEM

# Penentuan model terbaik
LR_stat_sar <- 2 * (logLik(sar_model)[1] - logLik(ols_model)[1])
pval <- pchisq(LR_stat_sar, df = 1, lower.tail = FALSE)
cat("LR Statistic:", LR_stat_sar, "p-value:", pval, "\n")
## LR Statistic: 5.940432 p-value: 0.01479739
LR_stat_sem <- 2 * (logLik(sem_model)[1] - logLik(ols_model)[1])
pval <- pchisq(LR_stat_sem, df = 1, lower.tail = FALSE)
cat("LR Statistic:", LR_stat_sem, "p-value:", pval, "\n")
## LR Statistic: 1.884821 p-value: 0.1697872
best_model <- names(which.min(c("OLS"=AIC_OLS, "SAR"=AIC_SAR, "SEM"=AIC_SEM)))
print(paste("Model terbaik berdasarkan AIC adalah:", best_model))
## [1] "Model terbaik berdasarkan AIC adalah: SAR"
# Residual
res_ols <- residuals(ols_model)
res_sar <- residuals(sar_model)
res_sem <- residuals(sem_model)
Bandung@data$res_ols <- res_ols
Bandung@data$res_sar <- res_sar
Bandung@data$res_sem <- res_sem
library(RColorBrewer)

# Skema warna
col_scheme <- brewer.pal(8, "RdBu")

# Plot residual OLS
spplot(Bandung, "res_ols", 
       main = "Residual OLS Model", 
       col.regions = rev(col_scheme))

# Plot residual SAR
spplot(Bandung, "res_sar", 
       main = "Residual SAR Model", 
       col.regions = rev(col_scheme))

# Plot residual SEM
spplot(Bandung, "res_sem", 
       main = "Residual SEM Model", 
       col.regions = rev(col_scheme))

5 Contoh Perhitungan SAR dan SEM dengan Data Diare

Diare <- read.csv("Diare.csv", sep = ";")
head(Diare)
##   id       Kecamatan Diare  PHBS Kepadatan Penduduk Prevalensi    E   SIR Biner
## 1  1           Andir   643 48.58       235    99288    0.64761 1119 0.575     0
## 2  2        Antapani   885 72.12       188    79496    1.11326  896 0.988     0
## 3  3       Arcamanik   957 72.56       102    77750    1.23087  876 1.093     1
## 4  4     Astanaanyar  1859 41.53       274    73495    2.52942  828 2.245     1
## 5  5 Babakan Ciparay  1374 71.88       200   141196    0.97312 1591 0.864     0
## 6  6   Bandung Kidul  2898 54.29       112    60596    4.78249  683 4.245     1
# Mapping
Bandung$id<-c(1:30)
Bandung_sf<-st_as_sf(Bandung)
Bandung_merged <- Bandung_sf %>%
  left_join(Diare, by = "id") 

ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare),color=NA) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ 
theme(legend.position = "right",
axis.text.x = element_blank(),  # Remove x-axis labels
axis.text.y = element_blank())+  # Remove y-axis labels   
labs(title = "",
       fill = "Diare") 

summary(Diare$Diare)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   235.0   529.0   775.0   939.1  1119.2  2898.0
breaks <- c(-Inf, 550, 800, 1200,Inf)

# Define labels for each interval
labels <- c("Very Low", "Low", "High", "Very High")

# Create a new column with discretized Diare
Bandung_merged$Diare_Discrete <- cut(Bandung_merged$Diare, breaks = breaks, labels = labels, right = TRUE)
  

ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare_Discrete),color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow", 
                                "Low" = "orange",
                                "High" = "red",
                                "Very High" = "red3"))+
labs(fill = "Diare")+theme(legend.position = "right",
axis.text.x = element_blank(),  # Remove x-axis labels
axis.text.y = element_blank())+  # Remove y-axis labels   
labs(title = "",
       fill = "Diare")

5.1 Spatial Autocorrelation

5.1.1 Matriks Bobot Spasial

#Membuat matrix bobot
row.names(Bandung) <- as.character(1:30)
W <- poly2nb(Bandung, row.names(Bandung), queen=TRUE)  
WB <- nb2mat(W, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"  
WL<-nb2listw(W)#List neighbours 
CoordK<-coordinates(Bandung)
plot(Bandung, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Bandung), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)

5.1.2 Global Moran’s I

Global_Moran<-moran.test(Diare$Diare,WL)

Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  Diare$Diare  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.6592, p-value = 0.2549
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03685550       -0.03448276        0.01171163

5.1.3 Local Moran’s I

Local_Moran <- localmoran(Diare$Diare, WL)
   
quadr_data <- attr(Local_Moran, "quadr")
mean_values <- quadr_data$mean

Bandung_sf <- st_as_sf(Bandung)
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Bandung_sf$mean_values <- mean_values_char
ggplot() +
  geom_sf(data = Bandung_sf, aes(fill = mean_values)) +  # Fill based on 'mean_values'
  scale_fill_manual(values = c("Low-Low" = "blue", "High-Low" = "green", 
                               "Low-High" = "yellow", "High-High" = "red")) + 
  labs(title = "Local Moran's I Mean Values in Bandung", fill = "Mean Type") +
  theme_minimal()

5.2 Spatial Regression

5.2.1 OLS Regression

Diare.ols<-lm(Prevalensi~PHBS+Kepadatan, data=Diare)
summary(Diare.ols)
## 
## Call:
## lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0072 -0.5862 -0.1094  0.3048  3.2325 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.628067   0.998838   2.631   0.0139 *
## PHBS        -0.012999   0.011047  -1.177   0.2496  
## Kepadatan   -0.003325   0.002415  -1.377   0.1800  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.872 on 27 degrees of freedom
## Multiple R-squared:  0.07665,    Adjusted R-squared:  0.00825 
## F-statistic: 1.121 on 2 and 27 DF,  p-value: 0.3408
plot(Diare.ols)

5.2.2 LM Test

LM<-lm.LMtests(Diare.ols, WL, test="all")
## Please update scripts to use lm.RStests in place of lm.LMtests
print(LM)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
## 
## RSerr = 0.056764, df = 1, p-value = 0.8117
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
## 
## RSlag = 0.019686, df = 1, p-value = 0.8884
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
## 
## adjRSerr = 0.30273, df = 1, p-value = 0.5822
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
## 
## adjRSlag = 0.26565, df = 1, p-value = 0.6063
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
## 
## SARMA = 0.32242, df = 2, p-value = 0.8511

5.2.3 Spatial lag model (SAR)

sar.Diare<-lagsarlm(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(sar.Diare)
## 
## Call:lagsarlm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare, 
##     listw = WL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01483 -0.57146 -0.11372  0.30592  3.22645 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  2.7010837  1.0137718  2.6644 0.007713
## PHBS        -0.0133200  0.0104764 -1.2714 0.203579
## Kepadatan   -0.0033403  0.0022925 -1.4571 0.145095
## 
## Rho: -0.039432, LR test value: 0.02088, p-value: 0.88511
## Asymptotic standard error: 0.26773
##     z-value: -0.14728, p-value: 0.88291
## Wald statistic: 0.021692, p-value: 0.88291
## 
## Log likelihood: -36.86783 for lag model
## ML residual variance (sigma squared): 0.68362, (sigma: 0.82681)
## Number of observations: 30 
## Number of parameters estimated: 5 
## AIC: 83.736, (AIC for lm: 81.757)
## LM test for residual autocorrelation
## test value: 0.2879, p-value: 0.59157
sar2sls.Diare<-stsls(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(sar2sls.Diare)
## 
## Call:stsls(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare, 
##     listw = WL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.63976 -0.49031  0.13329  0.82213  3.64165 
## 
## Coefficients: 
##               Estimate Std. Error t value Pr(>|t|)
## Rho          2.6732720  3.8860330  0.6879   0.4915
## (Intercept) -2.3220679  7.3926166 -0.3141   0.7534
## PHBS         0.0087817  0.0367913  0.2387   0.8113
## Kepadatan   -0.0022674  0.0043759 -0.5182   0.6043
## 
## Residual variance (sigma squared): 2.188, (sigma: 1.4792)

5.2.4 Cek residual

Diare$Diare.ols.res<-resid(Diare.ols) #residuals ols

Diare$Diare.sar.res<-resid(sar.Diare) #residual sar


Bandung@data<-Diare

spplot(Bandung,"Diare.ols.res", at=seq(min(Bandung$Diare.ols.res,na.rm=TRUE),max(Bandung$Diare.ols.res,na.rm=TRUE),length=12),col.regions=rev(brewer.pal(11,"RdBu")))

spplot(Bandung,"Diare.sar.res", at=seq(min(Bandung$Diare.ols.res,na.rm=TRUE),max(Bandung$Diare.ols.res,na.rm=TRUE),length=12),col.regions=rev(brewer.pal(11,"RdBu")))

5.2.5 Menghitung Impact

impacts(sar.Diare, listw=WL)
## Impact measures (lag, exact):
##                 Direct     Indirect        Total
## PHBS      -0.013324391 0.0005097403 -0.012814651
## Kepadatan -0.003341392 0.0001278289 -0.003213563

5.2.6 Spatial Error model

errorsalm.Diare<-errorsarlm(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(errorsalm.Diare)
## 
## Call:errorsarlm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare, 
##     listw = WL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02631 -0.57603 -0.10839  0.30138  3.21357 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  2.6876712  0.9367007  2.8693 0.004114
## PHBS        -0.0139051  0.0103307 -1.3460 0.178301
## Kepadatan   -0.0033146  0.0022578 -1.4681 0.142086
## 
## Lambda: -0.070346, LR test value: 0.062348, p-value: 0.80282
## Asymptotic standard error: 0.27479
##     z-value: -0.256, p-value: 0.79795
## Wald statistic: 0.065535, p-value: 0.79795
## 
## Log likelihood: -36.8471 for error model
## ML residual variance (sigma squared): 0.68218, (sigma: 0.82594)
## Number of observations: 30 
## Number of parameters estimated: 5 
## AIC: 83.694, (AIC for lm: 81.757)
fgls.Diare<-GMerrorsar(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(fgls.Diare)
## 
## Call:GMerrorsar(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare, 
##     listw = WL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01180 -0.58178 -0.10406  0.30597  3.22663 
## 
## Type: GM SAR estimator
## Coefficients: (GM standard errors) 
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  2.6582711  0.9421772  2.8214 0.004781
## PHBS        -0.0134575  0.0104061 -1.2932 0.195930
## Kepadatan   -0.0033197  0.0022746 -1.4595 0.144438
## 
## Lambda: -0.034782 (standard error): 0.58547 (z-value): -0.059409
## Residual variance (sigma squared): 0.68315, (sigma: 0.82653)
## GM argmin sigma squared: 0.68955
## Number of observations: 30 
## Number of parameters estimated: 5