Pembentukan Matriks Pembobot Spasial
dan Moran’s Index
Sapta Hastho Ponco
04-08-2024
Regresi spasial adalah pendekatan regresi yang dirancang untuk data spasial atau data yang dipengaruhi oleh lokasi (spatial effect). Ada dua jenis yang perlu diperhatikan: dependensi spasial dan heterogenitas (keragaman) spasial.
Dependensi spasial berarti pengamatan pada lokasi \(i\) dipengaruhi oleh pengamatan di lokasi \(j\), dimana \(j ≠ i\).
Heterogenitas spasial disebabkan oleh efek lokasi yang acak, yakni perbedaan antara satu lokasi dengan lokasi lainnya.
Metode regresi spasial berkembang dari regresi linier klasik dengan penyesuaian pertimbangan pengaruh lokasi. Pada tahun 1979, Tobler menyatakan dalam Hukum Geografi Pertama bahwa segala sesuatu saling berhubungan, namun sesuatu yang berada lebih dekat cenderung memiliki pengaruh yang lebih besar dibandingkan dengan sesuatu yang lebih jauh (Anselin, 1988).
library(sf)
library(ggplot2)
library(ggspatial)Membaca shapefile data
library(dplyr)
spJatim <- st_read("Jatim.shp")## Reading layer `Jatim' from data source
## `D:\2. Pengembangan diri\1 Exercise Bagus\Spasial in R\Dependensi Spasial\Jatim.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 110.8988 ymin: -8.78036 xmax: 116.2701 ymax: -5.043025
## Geodetic CRS: WGS 84spJatim <- spJatim %>%
dplyr::select(-c("putus_kode","putus_kabu","putus_putu","putus_X1",
"putus_X2","putus_X3","putus_X4","putus_X5","putus_X6","putus_U","putus_V"))
# Menyimpan shapefile yang telah dimodifikasi
st_write(spJatim, "Jatim_cleaned.shp", delete_layer = TRUE)## Deleting layer `Jatim_cleaned' using driver `ESRI Shapefile'
## Writing layer `Jatim_cleaned' to data source
## `Jatim_cleaned.shp' using driver `ESRI Shapefile'
## Writing 38 features with 8 fields and geometry type Multi Polygon.names(spJatim)## [1] "fid" "idkab" "nmprov" "nmkab" "kdprov" "kdkab" "sumber"
## [8] "periode" "geometry"1 Visualisasi Peta
Berikut merupakan plot peta dasar dengan menampilkan koordinat dengan sistem proyeksi
ggplot(data = spJatim) +
geom_sf() +
theme_minimal() +
labs(title = "Peta Jawa Timur",
subtitle = "Visualisasi data shapefile",
caption = "Sumber: Data Shapefile") +
coord_sf() # Menampilkan koordinat dengan sistem proyeksi
Membuat Layout Peta
ggplot(data = spJatim) +
geom_sf() +
theme_minimal() +
labs(title = "Peta Jawa Timur",
subtitle = "Visualisasi data shapefile",
caption = "Sumber: Data Shapefile") +
coord_sf(crs = 4326) + #menggunakan proyeksi EPSG:4326
annotation_scale(location = "bl", width_hint = 0.5) + # Skala bar
annotation_north_arrow(location = "tr", which_north = "true",
pad_x = unit(0.5, "in"), pad_y = unit(0.5, "in")) # Kompas arah
Berikut merupakan variabel yang akan dijadikan sebagai atribut pada data shapefile.
data=read.csv("putus.csv", header=T, sep=",")
data[, 5:12] <- round(data[, 5:12], 2)
head(data,5)## kode idkab kabupaten y x1 x2 x3 x4 x5 x6 u v
## 1 1 3501 Pacitan 109 22.31 1.07 15.66 88.88 17.22 7.74 111.11 -8.20
## 2 2 3502 Ponorogo 152 29.13 1.21 15.56 91.72 11.72 11.03 111.47 -7.86
## 3 3 3503 Trenggalek 215 16.54 1.22 15.80 93.38 14.15 11.09 111.71 -8.05
## 4 4 3504 Tulungagung 160 23.87 1.83 15.38 95.10 9.37 20.64 111.91 -8.08
## 5 5 3505 Blitar 271 16.31 1.54 15.45 92.45 10.70 13.65 112.24 -8.09Menambahkan kolom atribut pada data shapefile
library(sf)
library(tmap)
spJatim$y=data$y
spJatim$x1=data$x1
spJatim$x2=data$x2
spJatim$x3=data$x3
spJatim$x4=data$x4
spJatim$x5=data$x5
spJatim$x6=data$x6
spJatim$u=data$u
spJatim$v=data$vMenampilkan peta dengan variabel tertentu
library(tmap)
map1 <- tm_shape(spJatim) +
tm_polygons(col = "x1") +
tm_layout(title = " ", legend.outside = TRUE)
map2 <- tm_shape(spJatim) +
tm_polygons(col = "x2") +
tm_layout(title = " ", legend.outside = TRUE)
map3 <- tm_shape(spJatim) +
tm_polygons(col = "x3") +
tm_layout(title = " ", legend.outside = TRUE)
map4 <- tm_shape(spJatim) +
tm_polygons(col = "x4") +
tm_layout(title = " ", legend.outside = TRUE)
map5 <- tm_shape(spJatim) +
tm_polygons(col = "x5") +
tm_layout(title = " ", legend.outside = TRUE)
map6 <- tm_shape(spJatim) +
tm_polygons(col = "x6") +
tm_layout(title = " ", legend.outside = TRUE)
# Menggabungkan semua peta menjadi satu panel
tmap_arrange(map1, map2, map3, map4, map5, map6)
names(spJatim)## [1] "fid" "idkab" "nmprov" "nmkab" "kdprov" "kdkab"
## [7] "sumber" "periode" "geometry" "y" "x1" "x2"
## [13] "x3" "x4" "x5" "x6" "u" "v"Mengklasifikasikan simbologi
Misalkan ingin menampilkan peta dengan 5 kelas secara klasifikasi equal interval dengan warna YlGnBu
st_crs(spJatim) <- 3395 # Mengatur proyeksi peta
library(classInt)
breaks1 <- classIntervals(spJatim$x1, n = 5, style = "jenks") #klasifikasi Natural Breaks
breaks2 <- classIntervals(spJatim$x1, n = 5, style = "quantile") #klasifikasi Kuadran
breaks3 <- classIntervals(spJatim$x1, n = 5, style = "equal") #klasifikasi Equal Interval
library(RColorBrewer)
colors1 <- brewer.pal(n = 5, name = "BuGn") # palet warna "BuGn" dengan 5 warna
colors2 <- brewer.pal(n = 5, name = "YlOrRd")
colors3 <- brewer.pal(n = 5, name = "YlGnBu")
breaks3x2 <- classIntervals(spJatim$x2, n = 5, style = "equal")
breaks3x3 <- classIntervals(spJatim$x3, n = 5, style = "equal")
breaks3x4 <- classIntervals(spJatim$x4, n = 5, style = "equal")
breaks3x5 <- classIntervals(spJatim$x5, n = 5, style = "equal")
breaks3x6 <- classIntervals(spJatim$x6, n = 5, style = "equal")
# Menampilkan peta dengan 5 kelas
sim1 <- tm_shape(spJatim) +
tm_polygons(col = "x1",
style = "fixed",
breaks = breaks3$brks,
palette = colors3,
title = "Variabel x1")
sim2 <- tm_shape(spJatim) +
tm_polygons(col = "x2",
style = "fixed",
breaks = breaks3x2$brks,
palette = colors3,
title = "Variabel x2")
sim3 <- tm_shape(spJatim) +
tm_polygons(col = "x3",
style = "fixed",
breaks = breaks3x3$brks,
palette = colors3,
title = "Variabel x3")
sim4 <- tm_shape(spJatim) +
tm_polygons(col = "x4",
style = "fixed",
breaks = breaks3x4$brks,
palette = colors3,
title = "Variabel x4")
sim5 <- tm_shape(spJatim) +
tm_polygons(col = "x5",
style = "fixed",
breaks = breaks3x5$brks,
palette = colors3,
title = "Variabel x5")
sim6 <- tm_shape(spJatim) +
tm_polygons(col = "x6",
style = "fixed",
breaks = breaks3x6$brks,
palette = colors3,
title = "Variabel x6")
tmap_arrange(sim1, sim2, sim3, sim4, sim5, sim6)
2 Pembobot Spasial
Pembobot spasial digunakan untuk menentukan bobot antar lokasi yang diamati berdasarkan hubungan ketetanggaan antar lokasi.
Pembobot spasial terdiri dari dua jenis utama, yaitu pembobot titik dan pembobot area atau poligon. Pembobot titik berfokus pada kedekatan antara titik-titik tertentu dalam ruang, dimana jarak antara titik-titik tersebut dapat dihitung menggunakan berbagai pendekatan. Beberapa metode yang umum digunakan untuk menghitung pembobot titik termasuk K-NN (K-Nearest Neighbors), Inverse Distance Weighting, dan Exponential Decay. Metode K-NN mempertimbangkan sejumlah titik terdekat dari suatu titik sebagai tetangganya, sementara metode Inverse Distance memberikan bobot yang lebih besar pada tetangga yang lebih dekat dan semakin kecil seiring bertambahnya jarak. Metode Exponential Decay juga menekankan pada kedekatan tetapi dengan penurunan yang lebih tajam dibandingkan dengan Inverse.
Di sisi lain, pembobot area/poligon digunakan untuk mengukur keterhubungan antar area yang berbatasan, sering disebut juga sebagai contiguity. Ada beberapa pendekatan dalam membentuk hubungan tetangga antar area. Tiga metode yang paling dikenal adalah Bishop, Rook, dan Queen.
Menurut Kosfeld dalam Wuryandari dkk. (2014), ketetanggaan didefinisikan:
Bishop contiguity
Daerah pengamatannya ditentukan berdasarkan sudut-sudut yang saling bersinggungan dan sisi tidak diperhitungkan.
Rook contiguity
Daerah pengamatannya ditentukan berdasarkan sisi-sisi yang saling bersinggungan dan sudut tidak diperhitungkan.
Queen contiguity
Daerah pengamatannya ditentukan berdasarkan sisi-sisi yang saling bersinggungan dan sudut juga diperhitungkan. Menurut Kosfeld seperti yang dikutip oleh Wuryandari dkk. (2014), terdapat dua metode untuk memperoleh matriks pembobot spasial \(W\): matriks pembobot terstandarisasi (standardized contiguity matrix, \(W\)) dan matriks pembobot tak terstandarisasi (unstandardized contiguity matrix, \(W∗\)). Matriks pembobot terstandarisasi \(W\) adalah matriks yang memberikan bobot yang sama rata untuk semua tetangga terdekat, sementara bobot untuk lokasi lainnya adalah nol. Sebaliknya, matriks pembobot tak terstandarisasi \(W∗\) memberikan bobot satu untuk tetangga terdekat dan nol untuk lokasi lainnya.
2.1 Pembobot titik
library(spdep)
library(sf)
dist_matrix <- st_distance(spJatim,method="euclidean")
head(dist_matrix) # Jarak antara semua pasangan titik## Units: [m]
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.0000000 0.00000000 0.00000000 0.2880656 0.5556593 0.46609393 0.9102056
## [2,] 0.0000000 0.00000000 0.00000000 0.0000000 0.2155658 0.01345717 0.5051232
## [3,] 0.0000000 0.00000000 0.00000000 0.0000000 0.1387075 0.05738263 0.4878544
## [4,] 0.2880656 0.00000000 0.00000000 0.0000000 0.0000000 0.00000000 0.2490882
## [5,] 0.5556593 0.21556584 0.13870752 0.0000000 0.0000000 0.00000000 0.0000000
## [6,] 0.4660939 0.01345717 0.05738263 0.0000000 0.0000000 0.00000000 0.0000000
## [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1.4737908 1.8589155 2.433741 2.206950 2.189259 1.5155593 1.2281414
## [2,] 1.1372097 1.5767212 2.150812 1.842029 1.789108 1.1495789 0.7949769
## [3,] 1.0538529 1.4649532 2.030932 1.782173 1.763669 1.0898082 0.8117514
## [4,] 0.7802590 1.1916557 1.756159 1.510022 1.499300 0.8210907 0.5871237
## [5,] 0.4304001 0.8606537 1.428534 1.144337 1.117336 0.4488264 0.2019538
## [6,] 0.5540529 0.9694937 1.577326 1.216589 1.150579 0.5416215 0.1558736
## [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 1.2292769 1.09606939 0.8533111 0.42885056 0.20521158 0.1262471 0.2939509
## [2,] 0.7733256 0.63542478 0.4017857 0.00000000 0.00000000 0.0000000 0.1841252
## [3,] 0.8339027 0.68819235 0.4664553 0.06763376 0.07780637 0.2606199 0.3763120
## [4,] 0.7198754 0.48080011 0.3779315 0.00000000 0.04656169 0.3143149 0.3638065
## [5,] 0.4129809 0.13207722 0.1318817 0.22125244 0.28331001 0.5467676 0.5992029
## [6,] 0.3132107 0.02555352 0.0000000 0.00000000 0.05894272 0.3310445 0.3669122
## [,22] [,23] [,24] [,25] [,26] [,27] [,28]
## [1,] 0.5840774 0.9415266 0.9445246 1.2129139 1.5755432 1.8570210 2.169537
## [2,] 0.3159278 0.6633919 0.5320498 0.7723283 1.1303819 1.3986642 1.708740
## [3,] 0.4197081 0.7524694 0.6067804 0.8392842 1.1949325 1.4560380 1.762333
## [4,] 0.3639805 0.6889934 0.5346355 0.7581082 1.0929094 1.2938774 1.554227
## [5,] 0.5522216 0.8131130 0.5604740 0.5051753 0.7854025 0.9160317 1.165785
## [6,] 0.2008448 0.4475923 0.2089945 0.3529397 0.6628561 0.8295935 1.106322
## [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 2.408819 0.62654291 0.71076327 1.16673638 1.7819196 1.5337001 1.1668069
## [2,] 1.947821 0.17014584 0.41250755 0.79038866 1.3785354 1.0996088 0.7097708
## [3,] 2.000734 0.21635553 0.28581662 0.74124221 1.3570310 1.1159200 0.7697221
## [4,] 1.771223 0.11445169 0.02121788 0.47712665 1.0978099 0.8810969 0.6551543
## [5,] 1.382091 0.09660213 0.00000000 0.09578288 0.7113539 0.4914709 0.4175966
## [6,] 1.341452 0.00000000 0.08240473 0.22208257 0.7406851 0.4618274 0.2586403
## [,36] [,37] [,38]
## [1,] 0.3379983 1.4430309 1.073112896
## [2,] 0.1237373 0.9894892 0.693195609
## [3,] 0.2846634 1.0510226 0.648284726
## [4,] 0.2974412 0.9362193 0.396721600
## [5,] 0.5334359 0.6007022 0.009135827
## [6,] 0.3108889 0.4810222 0.055912610Invers Distance Weight
Inverse Distance Weight biasanya dihitung sebagai kebalikan dari jarak antar titik, dengan bobot yang lebih besar diberikan kepada titik-titik yang lebih dekat.
Langkah-langkah:
Mengambil koordinat centroid dari setiap poligon.
Menghitung jarak antar semua centroid.
Menghitung bobot jarak terbalik (invers distance).
Membuat objek listw yaitu menggunakan bobot jarak terbalik untuk membuat objek listw.
library(spdep)
library(sp)
# 1. Mengambil koordinat centroid dari setiap poligon
coords <- st_centroid(st_geometry(spJatim)) # Mengambil koordinat centroid dari tiap polygon
coords <- st_coordinates(coords) # Mengubahnya menjadi matriks koordinat
# 2. Menghitung matriks jarak
dist_mat <- spDists(coords, longlat = TRUE) # Menghitung matriks jarak
# 4. Menghitung bobot invers
inv_dist_weights <- 1 / dist_mat # Membuat inverse distance weights
diag(inv_dist_weights) <- 0 # Mengatur diagonal menjadi 0 karena jarak dengan diri sendiri adalah 0
# 5. Membuat listw object (spatial weights list)
nb <- dnearneigh(coords, 0, Inf) # Neighbor list dengan jarak tak terbatas
# Membuat daftar bobot jarak terbalik untuk setiap tetangga
inv_dist_weights_list <- lapply(1:length(nb), function(i) {
neighbors <- nb[[i]]
if (length(neighbors) > 0) {
1 / dist_mat[i, neighbors]
} else {
NULL
}
})
listw <- nb2listw(nb, glist = inv_dist_weights_list, style = "W") # Ubah menjadi listw object (spatial weights list)
summary(listw)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 1406
## Percentage nonzero weights: 97.36842
## Average number of links: 37
## Link number distribution:
##
## 37
## 38
## 38 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 37 links
## 38 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 37 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 3.369623 153.6384spJatim$inv_dist_sum <- sapply(listw$weights, sum)
plot(st_geometry(spJatim), border="white", col='gray', sub="Inverse Distance Weights")
plot(listw, coords, add = TRUE, col = "red") # Menambahkan garis yang mewakili bobot jarak invers
Exponential Distance Weight
Untuk membuat Exponential Distance Weights dalam analisis spasial, dapat digunakan fungsi eksponensial untuk menghitung bobot dari jarak antar titik. Bobot ini bisa dihitung sebagai exp(\(-α\) * distance), dimana α adalah parameter yang mengontrol seberapa cepat bobot menurun seiring dengan meningkatnya jarak.
Berikut adalah langkah-langkah untuk membuat Exponential Distance Weights:
Menghitung jarak antar centroid: Gunakan metode yang sama seperti sebelumnya untuk menghitung jarak antara centroid dari poligon.
Menghitung bobot eksponensial: Terapkan fungsi eksponensial pada jarak yang telah dihitung.
Membuat objek listw: Menggunakan bobot yang telah dihitung dalam analisis spasial.
library(spdep)
# 1. Mengambil koordinat centroid dari setiap poligon
coords <- st_centroid(st_geometry(spJatim))
coords <- st_coordinates(coords)
# 2. Menghitung matriks jarak
dist_mat <- spDists(coords, longlat = TRUE)
# 3. Menentukan parameter α
alpha <- 0.1 # Contoh nilai α, bisa disesuaikan
# 4. Menghitung bobot eksponensial
exp_dist_weights <- exp(-alpha * dist_mat)
diag(exp_dist_weights) <- 0 # Mengatur diagonal menjadi 0 karena jarak dengan diri sendiri
# 5. Membuat listw object (spatial weights list)
nb <- dnearneigh(coords, 0, Inf) # Membuat neighbor list tanpa batasan jarak
exp_dist_weights_list <- lapply(1:length(nb), function(i) {
neighbors <- nb[[i]]
if (length(neighbors) > 0) {
exp_dist_weights[i, neighbors]
} else {
NULL
}
})
listw_exp <- nb2listw(nb, glist = exp_dist_weights_list, style = "W")
summary(listw_exp)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 1406
## Percentage nonzero weights: 97.36842
## Average number of links: 37
## Link number distribution:
##
## 37
## 38
## 38 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 37 links
## 38 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 37 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 29.58895 160.2562plot(st_geometry(spJatim), border="white", col='gray', sub="Exponential Distance Weights")
plot(listw_exp, coords, add = TRUE, col = "red") # Menambahkan garis
Radial Distance Weight
Untuk Radial Distance Weights, bobot dihitung berdasarkan jarak radial antara titik-titik atau poligon, dimana bobotnya bisa mengikuti pola Gaussian atau fungsi lain yang menurun seiring bertambahnya jarak dari titik pusat. Bobot radial sering digunakan dalam konteks analisis spasial untuk memodelkan pengaruh yang menyebar dari satu titik, misalnya dalam bentuk kerucut atau lingkaran.
Langkah-langkah untuk menghitungnya:
Menghitung jarak antar centroid poligon.
Menghitung bobot radial dengan fungsi tertentu, misalnya Gaussian radial atau bentuk lain seperti linear decay.
Membuat neighbor list untuk mendefinisikan tetangga dengan kriteria tertentu.
Membuat objek listw untuk menerapkan bobot tersebut dalam analisis spasial.
library(spdep)
library(sp)
# 1. Mengambil koordinat centroid dari setiap poligon
coords <- st_centroid(st_geometry(spJatim))
coords <- st_coordinates(coords)
# 2. Menghitung matriks jarak
dist_mat <- spDists(coords, longlat = TRUE)
# 3. Menghitung bobot radial distance menggunakan fungsi Gaussian
# Misalnya kita menggunakan Gaussian kernel dengan sigma (variance) tertentu
sigma <- 1 # Misalkan kita menetapkan sigma = 1 (Anda dapat mengatur sesuai kebutuhan)
radial_dist_weights <- exp(- (dist_mat^2) / (2 * sigma^2)) # Fungsi Gaussian
# 4. Menyetel diagonal menjadi 0 karena tidak ada bobot untuk jarak dengan diri sendiri
diag(radial_dist_weights) <- 0
# 5. Membuat listw object (spatial weights list)
nb <- dnearneigh(coords, 0, Inf)
# Membuat daftar bobot radial distance untuk setiap tetangga
radial_dist_weights_list <- lapply(1:length(nb), function(i) {
neighbors <- nb[[i]]
if (length(neighbors) > 0) {
exp(- (dist_mat[i, neighbors]^2) / (2 * sigma^2))
} else {
NULL
}
})
listw_radial <- nb2listw(nb, glist = radial_dist_weights_list, style = "W")
# 6. Menambahkan ke data dan visualisasi
spJatim$radial_dist_sum <- sapply(listw_radial$weights, sum)
plot(st_geometry(spJatim), border = "white", col = 'gray', sub = "Radial Distance Weights")
plot(listw_radial, coords, add = TRUE, col = "red") # Garis mewakili bobot radial
K-Nearest Neighbor Weight
Untuk membuat K-Nearest Neighbor (K-NN) Weights digunakan pendekatan berikut:
Menentukan tetangga terdekat: menggunakan metode k-nearest neighbors untuk menentukan k tetangga terdekat dari setiap titik.
Membuat bobot berbasis k-nearest neighbors: Setelah menentukan tetangga, dapat dibuat bobot berbasis jarak atau menggunakan bobot biner (1 untuk tetangga, 0 untuk yang lainnya).
library(spdep)
# 1. Mengambil koordinat centroid dari setiap poligon
coords <- st_centroid(st_geometry(spJatim))
coords <- st_coordinates(coords)
# 2. Menentukan nilai k (jumlah tetangga terdekat)
k <- 5 # Contoh nilai k, bisa disesuaikan
# 3. Menghitung k-nearest neighbors
knn <- knn2nb(knearneigh(coords, k = k))
# 4. Membuat listw object (spatial weights list) dengan bobot biner
listw_knn <- nb2listw(knn, style = "W")
summary(listw_knn)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 190
## Percentage nonzero weights: 13.15789
## Average number of links: 5
## Non-symmetric neighbours list
## Link number distribution:
##
## 5
## 38
## 38 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 5 links
## 38 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 with 5 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 13.2 158.24spJatim$knn_sum <- sapply(listw_knn$weights, sum)
plot(st_geometry(spJatim), border="white", col='gray', sub="K-NN Distance Weights")
plot(listw_knn, coords, add = TRUE, col = "red")
2.2 Pembobot area
Queen Contiguity
library(spdep)
# Mengkonversi data sf ke objek polygon neighbors
queen.nb <- poly2nb(spJatim)
summary(queen.nb)## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 130
## Percentage nonzero weights: 9.00277
## Average number of links: 3.421053
## 4 regions with no links:
## 30 31 32 35
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8
## 4 4 6 7 5 3 6 2 1
## 4 least connected regions:
## 26 29 33 34 with 1 link
## 1 most connected region:
## 7 with 8 linksqueen.listw <- nb2listw(queen.nb, zero.policy = TRUE) #area tanpa tetangga diberikan nilai nol
queen.jatim= queen.listw
bobot.queen = listw2mat(queen.listw)
#bobot.queenDari output terlihat banyaknya lokasi ada 38 yang menunjukkan banyaknya kabupaten/kota di Jawa Timur. Banyaknya Nonzero link menunjukkan elemen matriks pembobot yang tidak nol (persinggungan antar kabupaten/kota).
# Visualisasi tetangga berdasarkan aturan queen
plot(st_geometry(spJatim), border="gray", main="Tetangga Queen")
plot(queen.nb, st_coordinates(st_centroid(st_geometry(spJatim))), add=TRUE, col="blue") 
Rook Contiguity
rook.nb <- poly2nb(spJatim,queen=FALSE)
summary(rook.nb)## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 130
## Percentage nonzero weights: 9.00277
## Average number of links: 3.421053
## 4 regions with no links:
## 30 31 32 35
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8
## 4 4 6 7 5 3 6 2 1
## 4 least connected regions:
## 26 29 33 34 with 1 link
## 1 most connected region:
## 7 with 8 linksrook.listw <- nb2listw(rook.nb, zero.policy = TRUE) #convert nb to listw type
rook.jatim = rook.listw
bobot.rook = listw2mat(rook.listw)
#bobot.rook# Visualisasi tetangga berdasarkan aturan rook
plot(st_geometry(spJatim), border="gray", main="Tetangga Rook")
plot(rook.nb, st_coordinates(st_centroid(st_geometry(spJatim))), add=TRUE, col="blue") 
Bishop Contiguity
Dalam konteks analisis spasial, aturan bishop menganggap dua poligon
sebagai tetangga jika mereka hanya berbagi titik sudut (vertex), tanpa
berbagi sisi yang sama. Ini berlawanan dengan aturan rook (berbagi sisi)
dan queen (berbagi sisi atau titik sudut). Namun, secara default,
spdep::poly2nb hanya mendukung aturan rook dan queen.
Untuk aturan bishop, harus memodifikasi secara manual tetangga yang hanya berbagi titik sudut (vertex). Berikut adalah langkah-langkah untuk membuat tetangga bishop:
Menggunakan poly2nb untuk menemukan tetangga queen terlebih dahulu, karena queen mencakup tetangga yang berbagi titik sudut.
Menghapus tetangga rook dari tetangga queen untuk mendapatkan tetangga bishop.
Membuat neighbor list berdasarkan aturan bishop.
library(spdep)
# 1. Membuat neighbor list queen (tetangga berbagi sisi atau titik)
queen.nb <- poly2nb(spJatim, queen = TRUE)
# 2. Membuat neighbor list rook (tetangga berbagi sisi)
rook.nb <- poly2nb(spJatim, queen = FALSE)
# 3. Mengidentifikasi tetangga bishop dengan menghapus tetangga rook dari tetangga queen
bishop.nb <- lapply(1:length(queen.nb), function(i) {
setdiff(queen.nb[[i]], rook.nb[[i]]) # Tetangga bishop adalah tetangga queen minus tetangga rook
})
bishop.nb
# Konversi ke nb object
bishop.jatim <- nb2listw(nb = nb(bishop.nb))
# 4. Visualisasi tetangga bishop
plot(st_geometry(spJatim), border = "gray") # Plot batas poligon
plot(bishop.jatim, st_coordinates(st_centroid(st_geometry(spJatim))), add = TRUE, col = "green")Menyimpan Matriks Pembobot
Untuk menyimpan matriks pembobot, bobot harus diubah dalam format matriks, kemudian disimpan dalam file *.csv.
bobot.queen = listw2mat(queen.listw) #convert listw to matrix
write.csv(bobot.queen, "Matriks Bobot Queen.csv")
bobot.rook = listw2mat(rook.listw) #convert listw to matrix
write.csv(bobot.rook, "Matriks Bobot Rook.csv")3 Moran’s Index
Cara lain untuk mengetahui adanya dependensi spasial antar lokasi adalah dengan melakukan uji autokorelasi spasial dengan menggunakan statistik Moran’s I. Autokorelasi spasial adalah taksiran dari korelasi antar nilai amatan yang berkaitan dengan lokasi pada variabel yang sama. Jika terdapat pola sistematik dalam penyebaran sebuah variabel, maka terdapat autokorelasi spasial.
Indeks Moran \(I\) adalah ukuran autokorelasi spasial yang mengukur sejauh mana nilai-nilai variabel di lokasi tertentu mirip dengan nilai-nilai di lokasi-lokasi tetangganya. Rumus Indeks Moran didefinisikan sebagai berikut:
\[ I = \frac{N}{W} \cdot \frac{\sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \]
dimana, - \(N\) adalah jumlah lokasi atau wilayah.
\(x_i\) dan \(x_j\) adalah nilai dari variabel di lokasi \(i\) dan \(j\).
\(\bar{x}\) adalah rata-rata dari variabel.
\(w_{ij}\) adalah bobot spasial yang menunjukkan hubungan antara lokasi \(i\) dan \(j\).
\(W\) adalah jumlah dari semua bobot spasial \(w_{ij}\).
Penjelasan Rumus
Term di Numerator:
\[ \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x}) \]
Ini menghitung kontribusi dari semua pasangan lokasi \(i\) dan \(j\) dengan memperhitungkan bobot spasial \(w_{ij}\). Jika \(w_{ij}\) besar, kontribusi dari pasangan ini lebih besar dalam perhitungan autokorelasi.
Term di Denominator:
\[ \sum_{i=1}^N (x_i - \bar{x})^2 \]
Ini menghitung varians dari nilai-nilai variabel di seluruh lokasi, yang berfungsi sebagai skala untuk ukuran autokorelasi.
Kalkulasi Rasio:
\[ \frac{N}{W} \]
Mengatur skala dari hasil yang didapat agar indeks Moran berada dalam rentang yang dapat dibandingkan dengan hasil lain.
Hipotesis Moran’s I untuk melihat analisis autokorelasi spasial adalah sebagai berikut:
Hipotesis
\(H_0\): Tidak ada autokorelasi spasial dalam data.
\(H_1\): Ada autokorelasi spasial dalam data.
Interpretasi Hasil
- dekat 1: menunjukkan autokorelasi positif, dimana nilai-nilai tinggi atau rendah cenderung dikelompokkan.
- dekat 0: menunjukkan tidak adanya autokorelasi spasial, dimana nilai-nilai tersebar secara acak.
- dekat -1: menunjukkan autokorelasi negatif, dimana nilai-nilai tinggi dan rendah cenderung tersebar.
Berikut dicontohkan Indeks Moran dengan menggunakan Pembobot Queen.
moran.test(spJatim$x1, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x1
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 1.2235, p-value = 0.1106
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12964485 -0.03030303 0.01709022moran.test(spJatim$x2, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x2
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 3.2128, p-value = 0.0006573
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.38969935 -0.03030303 0.01709022moran.test(spJatim$x3, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x3
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 3.5805, p-value = 0.0001714
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.43777908 -0.03030303 0.01709022moran.test(spJatim$x4, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x4
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 4.0774, p-value = 2.277e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.50272901 -0.03030303 0.01709022moran.test(spJatim$x5, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x5
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 4.2339, p-value = 1.148e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.52319514 -0.03030303 0.01709022moran.test(spJatim$x6, queen.jatim, randomisation = FALSE, zero.policy = TRUE)##
## Moran I test under normality
##
## data: spJatim$x6
## weights: queen.jatim n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 2.3833, p-value = 0.008578
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.28127096 -0.03030303 0.01709022Dari output terlihat:
Variabel x1: nilai p-value > alpha sehingga H0 tidak ditolak. Hal ini menunjukkan tidak ada autokorelasi spasial antar lokasi.
Variabel x2: nilai p-value < alpha sehingga H0 ditolak. Hal ini menunjukkan ada autokorelasi spasial antar lokasi.
Variabel x3: nilai p-value < alpha sehingga H0 ditolak. Hal ini menunjukkan ada autokorelasi spasial antar lokasi.
Variabel x4: nilai p-value < alpha sehingga H0 ditolak. Hal ini menunjukkan ada autokorelasi spasial antar lokasi.
Variabel x5: nilai p-value < alpha sehingga H0 ditolak. Hal ini menunjukkan ada autokorelasi spasial antar lokasi.
Variabel x6: nilai p-value < alpha sehingga H0 ditolak. Hal ini menunjukkan ada autokorelasi spasial antar lokasi.
par(mfrow = c(2, 3))
moran.plot(spJatim$x1, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x1")
moran.plot(spJatim$x2, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x2")
moran.plot(spJatim$x3, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x3")
moran.plot(spJatim$x4, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x4")
moran.plot(spJatim$x5, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x5")
moran.plot(spJatim$x6, queen.jatim, zero.policy = TRUE, main = "Moran Plot of x6")
par(mfrow = c(1, 1))4 References
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers.
Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37(1-2), 17-23.
T. Wuryandari, A. Hoyyi, D. S. Kusumawardani, and D. Rahmawati. 2014. Identifikasi Autokorelasi Spasial pada Jumlah Pengangguran di Jawa Tengah Menggunakan Indeks Moran.
Direktorat Statistik Kesejahteraan Rakyat, BPS, saptahas@bps.go.id