Matriks Pembobot Spasial dan Autokorelasi Spasial

Author

Rosita Ria Rusesta, S.Stat - Statistika dan Sains Data, IPB

Review: Proksimitas & Ketetanggaan

1. Konsep dasar: Mengapa kedekatan penting dalam data spasial?

Dalam data spasial, setiap pengamatan memiliki lokasi. Karena itu, nilai suatu pengamatan tidak selalu dapat dianggap independen dari pengamatan lain.

Landasan konseptual yang sering digunakan adalah Tobler’s First Law of Geography, yaitu bahwa semua objek saling berkaitan, tetapi objek yang lebih dekat cenderung memiliki hubungan yang lebih kuat dibandingkan objek yang lebih jauh.

Misalnya:

Kecamatan yang berdekatan kemungkinan memiliki tingkat kemiskinan yang lebih mirip dibandingkan kecamatan yang berjauhan

Tetapi kata “kemungkinan” penting. Kedekatan geografis tidak otomatis berarti dua lokasi pasti memiliki nilai yang sama. Kedekatan adalah cara kita mendefinisikan struktur hubungan spasial yang kemudian perlu dianalisis.

2. Apa itu proksimitas?

Proksimitas (proximity) adalah kedekatan antar area pengamatan berdasarkan proses pembangkitan datanya

Dalam analisis geostatistik cross-sectional dengan data yang memiliki point support, jarak umumnya menjadi ukuran kedekatan yang masuk akal untuk merepresentasikan proses pembentukan data.

Sementara itu, pada data areal, kedekatan dapat ditentukan berdasarkan hubungan ketetanggaan antarwilayah, misalnya dengan melihat apakah dua wilayah berbatasan secara langsung(contiguity). Namun, jarak juga dapat digunakan sebagai ukuran kedekatan pada data areal, misalnya berdasarkan jarak antara titik pusat (centroid) dua wilayah. Dengan demikian, pemilihan ukuran proksimitas perlu disesuaikan dengan karakteristik dan proses spasial yang mendasari data.

Jadi, proksimitas bisa didasarkan pada:

  • jarak,

  • berbagi batas wilayah,

  • jumlah tetangga terdekat,

  • atau struktur hubungan lainnya.

Ketetanggan Spasial

Dalam analisis data spasial, salah satu konsep penting adalah ketetanggaan spasial (spatial neighborhood) yang digunakan untuk menentukan wilayah mana yang dianggap memiliki kedekatan atau hubungan spasial dengan wilayah lainnya. Dengan menentukan ketetanggaan, kita dapat mengeksplorasi apakah wilayah yang berdekatan cenderung memiliki nilai yang serupa maupun berbeda, serta menganalisis keberadaan autokorelasi spasial.

Wilayah yang dianggap sebagai tetangga dapat ditentukan dengan cara yang berbeda, bergantung pada karakteristik data dan proses spasial yang ingin direpresentasikan. Secara umum, dua pendekatan yang banyak digunakan adalah berdasarkan jarak (distance-based) dan berdasarkan contiguity 

1. Ketetanggaan berdasarkan Jarak (Distance-based)

Ketetanggaan juga dapat ditentukan berdasarkan jarak antarwilayah. Dalam pendekatan ini, dua wilayah tidak harus berbagi batas untuk dianggap sebagai tetangga. Suatu wilayah dapat dianggap bertetangga apabila jaraknya berada dalam batas tertentu.

Sebagai contoh, dua kecamatan dapat didefinisikan sebagai tetangga apabila jarak antara titik pusat (centroid) kedua kecamatan kurang dari atau sama dengan 70 km. Pendekatan ini memungkinkan wilayah yang tidak berbatasan secara langsung tetap memiliki hubungan ketetanggaan apabila secara geografis cukup dekat.

Selain menggunakan batas jarak tertentu (distance threshold), ketetanggaan berbasis jarak juga dapat ditentukan menggunakan K-nearest neighbors (KNN), yaitu dengan menetapkan sejumlah k wilayah terdekat sebagai tetangga bagi setiap wilayah.

Dengan demikian, ketetanggaan spasial merupakan konsep untuk mendefinisikan hubungan kedekatan antarunit spasial, sedangkan contiguity dan distance-based merupakan dua pendekatan yang dapat digunakan untuk membentuk hubungan tersebut. Pemilihan pendekatan sebaiknya disesuaikan dengan karakteristik fenomena dan asumsi mengenai bagaimana hubungan spasial terjadi dalam data.

2. Ketetanggaan berdasarkan Contiguity

Pada pendekatan contiguity, hubungan spasial dianggap berkaitan dengan interaksi antarwilayah yang berbatasan secara langsung.

Terdapat dua bentuk contiguity yang umum digunakan, yaitu Rook dan Queen.

  • Rook contiguity menganggap dua wilayah bertetangga apabila keduanya berbagi bagian sisi atau garis batas.

  • Queen contiguity menganggap dua wilayah bertetangga apabila keduanya berbagi sisi atau bertemu pada satu titik sudut.

Matriks Pembobot Spasial

Setelah menentukan ketetanggaan spasial (spatial neighborhood), kita dapat membangun matriks ketetanggaan spasial untuk merepresentasikan hubungan antarwilayah. Matriks ini kemudian digunakan dalam analisis autokorelasi spasial. Setiap elemen matriks menunjukkan bobot spasial (spatial weight) yang menggambarkan hubungan antara suatu wilayah dengan wilayah lainnya.

1. Berdasarkan Jarak (Distance Based)

1.1 Radial Distance Weight (Distance Band Weight)

Matriks bobot radial distance mendefinisikan tetangga berdasarkan jarak geografis. Setiap titik atau wilayah dianggap sebagai tetangga jika jaraknya dari titik/wilayah pusat berada dalam jarak radius tertentu.

Contoh: Jika radius yang ditentukan adalah 70 km, maka setiap titik atau wilayah yang berada dalam jarak 70 km dari titik pusat akan dianggap sebagai tetangga, dan bobot matriks akan diberikan 1.

library(sf)
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(spdep)
Loading required package: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
library(spatialreg)
Loading required package: Matrix

Attaching package: 'spatialreg'
The following objects are masked from 'package:spdep':

    get.ClusterOption, get.coresOption, get.mcOption,
    get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
    set.coresOption, set.mcOption, set.VerboseOption,
    set.ZeroPolicyOption
library(sp)
library(tmap)
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Columbus Data Set

data(columbus)
glimpse(columbus)
Rows: 49
Columns: 22
$ AREA       <dbl> 0.309441, 0.259329, 0.192468, 0.083841, 0.488888, 0.283079,…
$ PERIMETER  <dbl> 2.440629, 2.236939, 2.187547, 1.427635, 2.997133, 2.335634,…
$ COLUMBUS.  <int> 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
$ COLUMBUS.I <int> 5, 1, 6, 2, 7, 8, 4, 3, 18, 10, 38, 37, 39, 40, 9, 36, 11, …
$ POLYID     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ NEIG       <int> 5, 1, 6, 2, 7, 8, 4, 3, 18, 10, 38, 37, 39, 40, 9, 36, 11, …
$ HOVAL      <dbl> 80.467, 44.567, 26.350, 33.200, 23.225, 28.750, 75.000, 37.…
$ INC        <dbl> 19.531, 21.232, 15.956, 4.477, 11.252, 16.029, 8.438, 11.33…
$ CRIME      <dbl> 15.725980, 18.801754, 30.626781, 32.387760, 50.731510, 26.0…
$ OPEN       <dbl> 2.850747, 5.296720, 4.534649, 0.394427, 0.405664, 0.563075,…
$ PLUMB      <dbl> 0.217155, 0.320581, 0.374404, 1.186944, 0.624596, 0.254130,…
$ DISCBD     <dbl> 5.03, 4.27, 3.89, 3.70, 2.83, 3.78, 2.74, 2.89, 3.17, 4.33,…
$ X          <dbl> 38.80, 35.62, 39.82, 36.50, 40.01, 43.75, 33.36, 36.71, 43.…
$ Y          <dbl> 44.07, 42.38, 41.18, 40.52, 38.00, 39.28, 38.41, 38.71, 35.…
$ AREA       <dbl> 10.3910, 8.6210, 6.9810, 2.9080, 16.8270, 8.9290, 16.2200, …
$ NSA        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,…
$ NSB        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ EW         <dbl> 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,…
$ CP         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0,…
$ THOUS      <dbl> 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,…
$ NEIGNO     <dbl> 1005, 1001, 1006, 1002, 1007, 1008, 1004, 1003, 1018, 1010,…
$ PERIM      <dbl> 2.440629, 2.236939, 2.187547, 1.427635, 2.997133, 2.335634,…
columbus.map <- st_read(system.file("shapes/columbus.gpkg", package="spData"))
Reading layer `columbus' from data source 
  `/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/spData/shapes/columbus.gpkg' 
  using driver `GPKG'
Simple feature collection with 49 features and 20 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 5.874907 ymin: 10.78863 xmax: 11.28742 ymax: 14.74245
Projected CRS: Undefined Cartesian SRS with unknown unit
tm_shape(columbus.map) +
  tm_polygons("HOVAL") +
  tm_layout(title = "Map of HOVAL")
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

Ekstraksi Koordinat Latitude dan Longitude

koord <- st_coordinates(st_centroid(st_geometry(columbus.map)))
head(koord)
            X        Y
[1,] 8.827218 14.36908
[2,] 8.332658 14.03162
[3,] 9.012265 13.81972
[4,] 8.460801 13.71696
[5,] 9.007982 13.29637
[6,] 9.739926 13.47463
w.rdw <-dnearneigh(koord,0,70,longlat=TRUE)
Warning in dnearneigh(koord, 0, 70, longlat = TRUE): neighbour object has 2
sub-graphs
w.rdw.s <- nb2listw(w.rdw,style='W')

dnearneigh() digunakan untuk menentukan tetangga berdasarkan rentang jarak tertentu

nb2listw() mengubah neighbour list menjadi spatial weights list

style = "W"untuk bobot spasial yang distandardisasi berdasarkan baris.

# Plot poligon
plot(st_geometry(columbus.map), border="blue", col="gray",
     main="Radial Distance Weight dengan dmax=70")

# Tambahkan titik centroid
points(koord, col="black")

# Tambahkan garis ketetanggaan
plot(w.rdw, koord, add=TRUE, col="red")

1.2 Invers Distance Weight

Matriks pembobot Inverse Distance Weight (IDW) atau power distance weight menentukan bobot berdasarkan kebalikan dari jarak antar titik atau wilayah. Semakin dekat dua titik, semakin besar bobot yang diberikan.

Contoh: Jika jarak antara dua titik adalah 5 km, dan eksponen jarak adalah 2, maka bobot dihitung sebagai 0.04.

Ekstraksi Koordinat Latitude dan Longitude

koord <- st_coordinates(st_centroid(columbus.map))
Warning: st_centroid assumes attributes are constant over geometries

Menghitung Jarak Euclidean Distance dari Setiap Area

D<-dist(koord, method = "euclidean")
D<-as.matrix(D)
#inverse weight matrix
alpha1=1
w.idw=1/(D^alpha1)

# inverse weight matrix - row-normalized
diag(w.idw)<-0
rtot<-rowSums(w.idw, na.rm =T)
w.idw.sd<-w.idw/rtot
rowSums(w.idw.sd, na.rm=T)
 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 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
# matriks pembobot spasial
w.idw.s <- mat2listw(w.idw.sd,style='W')  
summary(w.idw.s)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 2352 
Percentage nonzero weights: 97.95918 
Average number of links: 48 
Link number distribution:

48 
49 
49 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 39 40 41 42 43 44 45 46 47 48 49 with 48 links
49 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 39 40 41 42 43 44 45 46 47 48 49 with 48 links

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 49 2401 49 3.236101 197.9036
#invers distance weight - plot
plot(st_geometry(columbus.map), border="blue",col='gray', main="Invers Distance Weight dengan alpha=1")
plot(w.idw.s, koord, add = TRUE, col = "red")

1.3 Exponential Distance Weight

Matriks exponential distance weight memberikan bobot yang menurun secara eksponensial seiring dengan peningkatan jarak antar unit spasial. Ini mirip dengan inverse distance weight, tetapi penurunan bobot terhadap jarak lebih tajam.

Alpha=1

alpha=1
w.edw<-exp((-alpha)*D)

#row-normalized
diag(w.edw)<-0
rtot<-rowSums(w.edw,na.rm=TRUE)
w.edw.sd<-w.edw/rtot #row-normalized
w.edw.s = mat2listw(w.edw.sd,style='W')
w.edw.s <- mat2listw(w.edw.sd,style='W') #untuk melihat matriks W
summary(w.edw.s)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 2352 
Percentage nonzero weights: 97.95918 
Average number of links: 48 
Link number distribution:

48 
49 
49 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 39 40 41 42 43 44 45 46 47 48 49 with 48 links
49 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 39 40 41 42 43 44 45 46 47 48 49 with 48 links

Weights style: W 
Weights constants summary:
   n   nn S0       S1      S2
W 49 2401 49 3.297375 198.343
#exponential distance weight - plot
plot(st_geometry(columbus.map), border="blue",col='gray', main="Exponential Distance Weight dengan alpha=1")
plot(w.edw.s, koord, add = TRUE, col = "red")

1.4 Jarak Pangkat Ganda (Double Power Distance)

# 1. Hitung koordinat centroid dan matriks jarak Euclidean (d_ij)
centroids <- st_centroid(columbus.map)
Warning: st_centroid assumes attributes are constant over geometries
dist_matrix <- as.matrix(st_distance(centroids)) # Matriks d_ij (dalam satuan peta)

# 2. Menentukan parameter d (threshold jarak) dan k (pangkat)
# Menghitung jarak minimum agar tidak ada unit yang terisolasi (no islands)
k_nn <- knearneigh(st_coordinates(centroids), k = 1)
nb_nn <- knn2nb(k_nn)
Warning in knn2nb(k_nn): neighbour object has 13 sub-graphs
d_min <- max(unlist(nbdists(nb_nn, st_coordinates(centroids))))

d_cutoff <- d_min * 1.5  # Batas threshold d (misal: 1.5x jarak tetangga terdekat)
k_param  <- 2            # Parameter pangkat k (dapat disesuaikan, misal k = 2 atau k = 3)

# 3. Fungsi Matriks Pembobot Spasial Jarak Pangkat Ganda (Double Power Distance)
calc_double_power_weight <- function(D, d, k) {
  W <- matrix(0, nrow = nrow(D), ncol = ncol(D))
  
  # Kondisi 0 <= d_ij <= d
  mask <- (D <= d)
  
  # Rumus: [1 - (d_ij / d)^k]^k
  W[mask] <- (1 - (D[mask] / d)^k)^k
  
  # Elemen diagonal utama w_ii = 0
  diag(W) <- 0
  
  return(W)
}

W_matrix <- calc_double_power_weight(D = dist_matrix, d = d_cutoff, k = k_param)

# 4. Standarisasi Baris
W_matrix_std <- W_matrix / rowSums(W_matrix)

# 5. Konversi ke objek listw untuk analisis spdep (Moran's I, Regresi Spasial, dll)
listw_dpd <- mat2listw(W_matrix_std, style = "W")

summary(listw_dpd)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 506 
Percentage nonzero weights: 21.07455 
Average number of links: 10.32653 
Link number distribution:

 3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 
 2  5  5  4  2  1  2  3  1  3  5  6  3  3  1  2  1 
2 least connected regions:
1 47 with 3 links
1 most connected region:
25 with 19 links

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 49 2401 49 19.00926 199.6745
plot(st_geometry(columbus.map), 
     col = "grey95", 
     border = "grey60", 
     main = "Jaringan Konektivitas Spasial\n(Double Power Distance)")

# Tambahkan garis hubungan antar-centroid (tetangga dengan bobot > 0)
plot(listw_dpd$neighbours, 
     st_coordinates(centroids), 
     add = TRUE, 
     col = "red", 
     lwd = 0.8, 
     pch = 19, 
     cex = 0.6)

1.5 K-Nearest Neighbor

Matriks bobot k-nearest neighbor menentukan tetangga berdasarkan sejumlah tetangga terdekat (k) yang dipilih untuk setiap titik atau wilayah.

Contoh: Misalkan k = 3, artinya setiap titik akan memiliki 3 tetangga terdekat. Matriks bobot akan memiliki 1 di elemen yang mewakili tetangga terdekat, dan 0 di elemen lainnya.

# k = 5
w.knn<-knn2nb(knearneigh(koord,k=5,longlat=TRUE))
w.knn.s <- nb2listw(w.knn,style='W')
summary(w.knn.s)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 245 
Percentage nonzero weights: 10.20408 
Average number of links: 5 
Non-symmetric neighbours list
Link number distribution:

 5 
49 
49 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 39 40 41 42 43 44 45 46 47 48 49 with 5 links
49 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 39 40 41 42 43 44 45 46 47 48 49 with 5 links

Weights style: W 
Weights constants summary:
   n   nn S0   S1     S2
W 49 2401 49 17.4 201.92
#k-nearest neighbor distance weight - plot
plot(st_geometry(columbus.map), border="blue",col='gray', main="K-Nearest Neigbor dengan k=5")
plot(w.knn.s, koord, add = TRUE, col = "red")

2. Berdasarkan Contiguity

Ilustrasi :

2.1 Rook

Matriks pembobot spasial rook contiguity didasarkan pada konsep bahwa dua wilayah dianggap bertetangga jika mereka berbagi sisi batas.

rook <- poly2nb(as(columbus.map, "Spatial"), queen=FALSE)
summary(rook)
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 200 
Percentage nonzero weights: 8.329863 
Average number of links: 4.081633 
Link number distribution:

 2  3  4  5  6  7  9 
 7 10 17  8  3  3  1 
7 least connected regions:
1 6 31 39 42 46 47 with 2 links
1 most connected region:
20 with 9 links
w.rook <- nb2mat(rook, style="B")
head(w.rook,3)
  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
1 0 1 1 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2 1 0 1 1 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3 1 1 0 1 1 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#standarisasi
diag(w.rook)<-0
rtot<-rowSums(w.rook, na.rm =T)
w.rook.sd<-w.rook/rtot
rowSums(w.rook.sd, na.rm=T)
 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 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
w.rook <- nb2mat(rook, style="B")
head(w.rook,3)
  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
1 0 1 1 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2 1 0 1 1 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3 1 1 0 1 1 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#standarisasi
diag(w.rook)<-0
rtot<-rowSums(w.rook, na.rm =T)
w.rook.sd<-w.rook/rtot
rowSums(w.rook.sd, na.rm=T)
 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 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
plot(st_geometry(columbus.map), border="blue",col='gray')
coords <- st_coordinates(st_centroid(st_geometry(columbus.map)))
plot(rook, koord, add = TRUE, col = "red")

2.2 Queen

Matriks pembobot spasial queen contiguity lebih luas daripada rook contiguity karena wilayah-wilayah dianggap bertetangga jika mereka berbagi sisi atau titik sudut.

queen <- poly2nb(as(columbus.map, "Spatial"),queen = TRUE)
summary(queen)
Neighbour list object:
Number of regions: 49 
Number of nonzero links: 236 
Percentage nonzero weights: 9.829238 
Average number of links: 4.816327 
Link number distribution:

 2  3  4  5  6  7  8  9 10 
 5  9 12  5  9  3  4  1  1 
5 least connected regions:
1 6 42 46 47 with 2 links
1 most connected region:
20 with 10 links
w.queen <- nb2mat(queen, style="B")
head(w.queen, 3)
  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
1 0 1 1 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2 1 0 1 1 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3 1 1 0 1 1 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#standarisasi
diag(w.queen)<-0
rtot<-rowSums(w.queen, na.rm =T)
w.queen.sd<-w.queen/rtot
rowSums(w.queen.sd, na.rm=T)
 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 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
plot(st_geometry(columbus.map), border="blue",col='gray')
coords <- st_coordinates(st_centroid(st_geometry(columbus.map)))
plot(queen, coords, add = TRUE, col = "red")

Autokorelasi Spasial

1. Autokorelasi Spasial Global

Autokorelasi spasial global mengukur kecenderungan seluruh wilayah geografis sebagai satu kesatuan. Hal ini membantu dalam mengevaluasi apakah ada pola spasial yang signifikan di seluruh dataset.

Autokorelasi spasial digunakan untuk menggambarkan sejauh mana suatu variabel memiliki keterkaitan dengan dirinya sendiri berdasarkan posisi geografisnya. Konsep ini berkaitan erat dengan Hukum Pertama Geografi Tobler (Tobler’s First Law of Geography), yang menyatakan bahwa:

“Segala sesuatu berhubungan dengan segala sesuatu yang lain, tetapi hal-hal yang berdekatan memiliki hubungan yang lebih kuat daripada hal-hal yang berjauhan.”
(Tobler, 1970)

Autokorelasi spasial positif terjadi ketika pengamatan yang memiliki nilai yang serupa berada lebih berdekatan satu sama lain, sehingga membentuk pola mengelompok (clustered). Sebaliknya, autokorelasi spasial negatif terjadi ketika pengamatan dengan nilai yang berbeda berada lebih berdekatan, sehingga membentuk pola menyebar atau terdispersi (dispersed).

1. Ruang Lattice (Lattice Space):

Menggunakan matriks pembobot spasial untuk menggambarkan hubungan antar lokasi. Korelasi spasial dapat dihitung dengan teknik seperti Moran’s I atau Geary’s C. Moran’s I mengukur keterkaitan nilai antar lokasi yang berdekatan.

2. Ruang Kontinu (Continuous Space):

Umumnya menggunakan semivariogram atau covariance function. Semivariogram mengukur perbedaan antara dua titik berdasarkan jaraknya, sehingga menunjukkan bagaimana hubungan spasial berubah seiring jarak.

A. Joint Count Statistics

Joint count merupakan metode paling dasar dalam menentukan autokorelasi spasial antar area. Statistik Joint Count adalah metode yang digunakan dalam analisis data spasial untuk mengukur autokorelasi spasial khususnya pada variabel kategorik biner.

Tipe Pasangan dalam Joint Count:

  1. BB (Black-Black): Pasangan tetangga di mana kedua nilai sama (misalnya, dua lokasi memiliki nilai 1 atau “hitam”).

  2. WW (White-White): Pasangan tetangga di mana kedua nilai sama (dua lokasi memiliki nilai 0 atau “putih”).

  3. BW (Black-White): Pasangan tetangga di mana nilai-nilai berbeda (satu lokasi bernilai 1 dan tetangganya bernilai 0).

Pola autokorelasi spasial ada tiga macam, yaitu pola mengelompok (cluster), pola menyebar (dispered), dan pola tidak beraturan (random).

  1. Pola mengelompok, banyakanya gabungan hitam-hitam (BB) atau putih-putih (WW) akan lebih banyak dibandingkan hitam-putih (BW), dan autokorelasi spasial bernilai positif.

  2. Pola menyebar, banyaknya gabungan hitam-putih (BW) lebih banyak dibandingkan gabungan hitam-hitam (BB) atau putih-putih (WW), dan mempunyai nilai autokorelasi negatif.

  3. Pola acak memiliki jumlah gabungan BB, WW, dan BW yang bersifat acak dengan nilai autokorelasi spasialnya mendekati nol.

pri <- rep(1,12)
seg <- rep(0,4)
ter <- rep(1,2)
cua <- rep(0,4)
qui <- rep(1,2)
sex <-  rep(0,12)

A <- matrix(c(pri, seg, ter, cua ,qui, sex), nrow=6, byrow=FALSE)

A
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    1    0    0    0    0
[2,]    1    1    0    0    0    0
[3,]    1    1    0    0    0    0
[4,]    1    1    0    0    0    0
[5,]    1    1    1    1    0    0
[6,]    1    1    1    1    0    0

Selanjutnya matriks A dikonversi menjadi raster menggunakan fungsi raster().

library(raster)

Attaching package: 'raster'
The following object is masked from 'package:dplyr':

    select
rA <- raster(A)

rA
class      : RasterLayer 
dimensions : 6, 6, 36  (nrow, ncol, ncell)
resolution : 0.1666667, 0.1666667  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 0, 1  (min, max)

Berikut ini adalah plot dari area A.

plot(rA)
text(coordinates(rA), labels=rA[ ], cex=1.5)

pA <- rasterToPolygons(rA, dissolve=FALSE)

pA
class       : SpatialPolygonsDataFrame 
features    : 36 
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs         : NA 
variables   : 1
names       : layer 
min values  :     0 
max values  :     1 
spA <- SpatialPolygons(pA@polygons)

nb1 <- poly2nb(spA, queen = T)

nb1
Neighbour list object:
Number of regions: 36 
Number of nonzero links: 220 
Percentage nonzero weights: 16.97531 
Average number of links: 6.111111 
par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb1, xy, col='red', lwd=2, add=TRUE)

nb2 <- poly2nb(spA, queen = FALSE)

nb2
Neighbour list object:
Number of regions: 36 
Number of nonzero links: 120 
Percentage nonzero weights: 9.259259 
Average number of links: 3.333333 
par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb2, xy, col='green', lwd=2, add=TRUE)

wl1 <- nb2listw(nb1, style='B')

wl2 <- nb2listw(nb2, style='B')

jc_test1 <- joincount.test(as.factor(pA$layer), wl1)

jc_test1

    Join count test under nonfree sampling

data:  as.factor(pA$layer) 
weights: wl1 

Std. deviate for 0 = 5.1529, p-value = 1.282e-07
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
             53.00000              33.17460              14.80263 


    Join count test under nonfree sampling

data:  as.factor(pA$layer) 
weights: wl1 

Std. deviate for 1 = 4.7634, p-value = 9.52e-07
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
             37.00000              20.95238              11.34999 
jc_test2 <- joincount.test(as.factor(pA$layer), wl2)

jc_test2

    Join count test under nonfree sampling

data:  as.factor(pA$layer) 
weights: wl2 

Std. deviate for 0 = 5.4677, p-value = 2.28e-08
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            30.000000             18.095238              4.740611 


    Join count test under nonfree sampling

data:  as.factor(pA$layer) 
weights: wl2 

Std. deviate for 1 = 5.1203, p-value = 1.525e-07
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            22.000000             11.428571              4.262554

Berdasarkan kedua output di atas, dengan p-value yang sangat kecil, artinya kita dapat menolak hipotesis nol yang menyatakan bahwa terdapat autokorelasi. Sesuai dengan output di atas, alternative hypothesis: greater, artinya kita dapat menyimpulkan bahwa terdapat cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf nyata 5%.

Pengujian hipotesis dapat pula dilakukan dengan melibatkan algoritma monte carlo seperti di bawah ini.

set.seed(123)
jc_test3 <- joincount.mc(as.factor(pA$layer), wl1, nsim=99)

jc_test3

    Monte-Carlo simulation of join-count statistic

data:  as.factor(pA$layer) 
weights: wl1 
number of simulations + 1: 100 

Join-count statistic for 0 = 53, rank of observed statistic = 100,
p-value = 0.01
alternative hypothesis: greater
sample estimates:
    mean of simulation variance of simulation 
              33.01010               15.15296 


    Monte-Carlo simulation of join-count statistic

data:  as.factor(pA$layer) 
weights: wl1 
number of simulations + 1: 100 

Join-count statistic for 1 = 37, rank of observed statistic = 100,
p-value = 0.01
alternative hypothesis: greater
sample estimates:
    mean of simulation variance of simulation 
              20.37374               12.37930

B. Indeks Moran

Nilai Indeks Moran dapat digunakan untuk menentukan pola dispersi/acak/gerombol.

  • Nilai positif (Moran’s I > 0): Menunjukkan bahwa lokasi yang berdekatan cenderung memiliki nilai yang mirip (pengelompokan spasial), misalnya wilayah yang berdekatan sama-sama memiliki tingkat kemiskinan yang tinggi.

  • Nilai negatif (Moran’s I < 0): Menunjukkan bahwa lokasi yang berdekatan cenderung memiliki nilai yang berbeda (penyebaran spasial), misalnya wilayah dengan tingkat kemiskinan tinggi cenderung berada di sebelah wilayah dengan tingkat kemiskinan rendah.

  • Nilai mendekati nol (Moran’s I ≈ 0): Menunjukkan bahwa nilai antar wilayah bersifat acak, tanpa adanya pola spasial yang jelas.

Sebelum melakukan pengujian perlu dicek terlebih dahulu apakah datanya mempunyai sebaran normal atau tidak.

Penerapan pada data Columbus

nb_queen <- poly2nb(columbus.map, queen = TRUE)
listw_queen <- nb2listw(nb_queen, style = "W")

#Uji Moran's I Analitis (misal pada variabel Pendapatan / INC)
moran_analitis <- moran.test(columbus$INC, listw_queen)
print(moran_analitis)

    Moran I test under randomisation

data:  columbus$INC  
weights: listw_queen    

Moran I statistic standard deviate = 4.7645, p-value = 9.467e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.415628778      -0.020833333       0.008391926

Hasil uji Moran’s I menunjukkan bahwa variabel pendapatan ($INC) pada wilayah Columbus memiliki nilai Moran’s I sebesar 0,4156 dengan p-value 9,467e-07. Nilai ini positif dan signifikan, sehingga dapat diinterpretasikan bahwa terdapat autokorelasi spasial positif yang kuat dalam distribusi pendapatan. Artinya, wilayah dengan pendapatan tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki pendapatan tinggi, sedangkan wilayah dengan pendapatan rendah cenderung berdekatan dengan wilayah berpendapatan rendah. Pola ini menandakan adanya kecenderungan pengelompokan (clustered)

Uji moran juga dapat dilakukan dengan melibatkan simulasi monte carlo.

set.seed(123) 
moran_mc <- moran.mc(columbus$INC, listw_queen, nsim = 999)
print(moran_mc)

    Monte-Carlo simulation of Moran I

data:  columbus$INC 
weights: listw_queen  
number of simulations + 1: 1000 

statistic = 0.41563, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
# 7. Visualisasi Moran Scatterplot
moran.plot(columbus$INC, listw_queen,
           xlab = "Income (INC)",
           ylab = "Spatially Lagged Income (W * INC)",
           main = "Moran Scatterplot - INC (Queen Contiguity)",
           pch = 19, col = "firebrick")

C. Global Geary’s C

Global Geary’s C adalah statistik yang digunakan untuk mengukur autokorelasi spasial global, serupa dengan Indeks Moran’s I, tetapi lebih sensitif terhadap perbedaan lokal di antara wilayah yang berdekatan. Geary C bervariasi pada skala dari 0 hingga 2.

  • C sekitar 1 menunjukkan tidak ada autokorelasi (pola acak)

  • C bernilai sekitar 0 menunjukkan autokorelasi positif sempurna (pola gerombol)

  • C bernilai sekitar 2 menunjukkan autokorelasi negatif sempurna (pola tersebar)

GC <- geary.test(columbus.map$INC, listw_queen, randomisation = T)
GC

    Geary C test under randomisation

data:  columbus.map$INC 
weights: listw_queen   

Geary C statistic standard deviate = 3.1756, p-value = 0.0007476
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic       Expectation          Variance 
       0.67588162        1.00000000        0.01041725

D. Getis-Ord G

Getis-Ord Global G menilai apakah nilai tinggi atau rendah dari variabel tertentu cenderung terkumpul di wilayah tertentu secara global (mencakup seluruh area studi). globalG.test(): Ini adalah fungsi yang digunakan untuk menjalankan Getis-Ord Global G test. Fungsi ini berasal dari paket spdep dalam R dan digunakan untuk menguji apakah ada pengelompokan nilai yang tinggi (hotspots) atau rendah (coldspots) di seluruh wilayah yang dianalisis.

GO <- globalG.test(columbus.map$INC,listw_queen, alternative = "greater")
Warning in globalG.test(columbus.map$INC, listw_queen, alternative =
"greater"): Binary weights recommended (especially for distance bands)
GO

    Getis-Ord global G statistic

data:  columbus.map$INC 
weights: listw_queen   

standard deviate = 3.4706, p-value = 0.0002596
alternative hypothesis: greater
sample estimates:
Global G statistic        Expectation           Variance 
      2.274862e-02       2.083333e-02       3.045420e-07

  • G tinggi (positif): Menunjukkan bahwa nilai tinggi dari variabel cenderung berkelompok (hotspots). Ini berarti area dengan nilai yang lebih tinggi dari variabel tertentu terkonsentrasi di beberapa wilayah.

  • G rendah (negatif): Menunjukkan bahwa nilai rendah cenderung berkelompok (coldspots), yaitu ada konsentrasi wilayah dengan nilai variabel yang rendah.

  • G mendekati 0: Menunjukkan tidak adanya pola yang signifikan atau distribusi acak dari nilai tinggi dan rendah.

2. Autokorelasi Spasial Lokal

Autokorelasi spasial lokal, di sisi lain, berfokus pada pola spasial di tingkat lokal. Ini membantu mengidentifikasi “hotspots” atau “coldspots” di wilayah tertentu, yaitu area yang secara signifikan memiliki nilai lebih tinggi atau lebih rendah dibandingkan dengan sekitarnya.

A. Local Moran’s I

Local Moran’s I (atau dikenal sebagai Local Indicators of Spatial Association, atau LISA) adalah statistik yang digunakan untuk mengukur autokorelasi spasial lokal. Ini merupakan variasi dari Moran’s I global, tetapi difokuskan pada pola spasial di tingkat lokal untuk mengidentifikasi hotspots (wilayah dengan nilai tinggi yang berdekatan) dan coldspots (wilayah dengan nilai rendah yang berdekatan), serta outliers (wilayah yang memiliki nilai yang berbeda dari tetangganya).

Nilai Local Moran’s I:

  • Nilai positif: Menunjukkan pengelompokan positif, di mana lokasi dengan nilai yang serupa (tinggi-tinggi atau rendah-rendah) cenderung berdekatan.

  • Nilai negatif: Menunjukkan outliers atau pengelompokan negatif, di mana lokasi dengan nilai yang berbeda (tinggi dikelilingi oleh rendah atau sebaliknya) berdekatan.

  • Nilai nol atau mendekati nol: Menunjukkan tidak ada hubungan spasial di tingkat lokal, atau data terdistribusi secara acak di lokasi tersebut.

LM <- localmoran(columbus.map$INC, listw_queen, alternative = "greater")
head(LM)
            Ii         E.Ii     Var.Ii        Z.Ii Pr(z > E(Ii))
1  0.682691373 -0.017381426 0.40954011  1.09394373    0.13698983
2 -0.226728724 -0.030741520 0.46596621 -0.28711144    0.61298650
3 -0.012500696 -0.001634356 0.01871230 -0.07943642    0.53165725
4 -0.176841847 -0.064052901 0.68751316 -0.13602729    0.55410014
5  0.301975916 -0.006376410 0.03302679  1.69673516    0.04487337
6  0.002287171 -0.001788759 0.04281544  0.01969820    0.49214206
columbus.map$lmI <- LM[, "Ii"] # local Moran's I
columbus.map$lmZ <- LM[, "Z.Ii"] # z-scores
# p-values corresponding to alternative greater
columbus.map$lmp <- LM[, "Pr(z > E(Ii))"]
p1 <- tm_shape(columbus.map) +
  tm_polygons(col = "lmI", title = "Local Moran's I",
              style = "quantile") +
  tm_layout(legend.outside = TRUE)
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
(instead of 'col'), and 'col' for the outlines (instead of 'border.col').
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
p1
[`tm_scale_intervals()`] Variable(s) "fill" contains positive and negative
values, so midpoint is set to 0. Set midpoint = NA to show the full range of
visual values.
This message is displayed once per session.

p2 <- tm_shape(columbus.map) +
  tm_polygons(col = "lmp", title = "p-value",
              breaks = c(-Inf, 0.05, Inf)) +
  tm_layout(legend.outside = TRUE)
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
the map variable `fill` namely 'breaks' to fill.scale = tm_scale(<HERE>).
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
p2

H0: Terdapat autokorelasi spasial

H1: Tidak terdapat autokorelasi spasial

Dalam uji dua sisi ini, nilai z-score yang lebih rendah dari –1,96 menunjukkan autokorelasi spasial negatif, dan nilai z-score yang lebih besar dari 1,96 menunjukkan autokorelasi spasial positif. Peta di bawah ini menunjukkan area dengan autokorelasi spasial negatif, tidak ada, dan positif, yang diperoleh dengan memecah legenda menurut nilai z-score.

tmap_mode("plot")
ℹ tmap modes "plot" - "view"
ℹ toggle with `tmap::ttm()`
tm_shape(columbus.map) + tm_polygons(col = "lmZ",
title = "Local Moran's I", style = "fixed",
breaks = c(-Inf, -1.96, 1.96, Inf),
labels = c("Negative SAC", "No SAC", "Positive SAC"),
palette =  c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
  'labels' to 'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
Multiple palettes called "blue" found: "kovesi.blue", "tableau.blue". The first one, "kovesi.blue", is returned.

breaks = c(-Inf, -1.96, 1.96, Inf): Rentang nilai untuk klasifikasi warna, dengan:

  • -Inf hingga -1.96 mewakili Negative Spatial Autocorrelation (SAC).

  • -1.96 hingga 1.96 mewakili No Spatial Autocorrelation (No SAC).

  • 1.96 hingga Inf mewakili Positive Spatial Autocorrelation (Positive SAC).

# Autokorelasi spasial positif (Z-score > 1.96)
positive_autocorr <- columbus.map[columbus.map$lmZ > 1.96, ]
# Autokorelasi spasial negatif (Z-score < -1.96)
negative_autocorr <- columbus.map[columbus.map$lmZ < -1.96, ]
# Area tanpa autokorelasi spasial (Z-score antara -1.96 dan 1.96)
no_autocorr <- columbus.map[columbus.map$lmZ >= -1.96 & columbus.map$lmZ <= 1.96, ]
moran.plot(columbus.map$INC, listw_queen)

  • Kuadran I (High-High): Lokasi dengan nilai tinggi yang dikelilingi oleh lokasi dengan nilai tinggi (indikasi hotspot).

  • Kuadran II (Low-High): Lokasi dengan nilai rendah yang dikelilingi oleh lokasi dengan nilai tinggi (indikasi outlier).

  • Kuadran III (Low-Low): Lokasi dengan nilai rendah yang dikelilingi oleh lokasi dengan nilai rendah (indikasi coldspot).

  • Kuadran IV (High-Low): Lokasi dengan nilai tinggi yang dikelilingi oleh lokasi dengan nilai rendah (indikasi outlier).

Interpretasi Moran Plot:

  • Kemiringan Positif: Jika garis regresi memiliki kemiringan positif, ini menunjukkan autokorelasi spasial positif, di mana lokasi dengan nilai tinggi cenderung berdekatan dengan lokasi bernilai tinggi, dan lokasi bernilai rendah cenderung berdekatan dengan lokasi bernilai rendah.

  • Kemiringan Negatif: Jika garis regresi memiliki kemiringan negatif, ini menunjukkan autokorelasi spasial negatif, di mana lokasi dengan nilai tinggi dikelilingi oleh lokasi dengan nilai rendah (dan sebaliknya).

B. Getis-Ord Gi

Getis-Ord 𝐺𝑖 adalah statistik spasial yang digunakan untuk mendeteksi adanya wilayah yang hotspot/coldspot.

Local_M <- localG(columbus.map$INC, listw_queen, alternative = "greater")
Local_M
 [1]  1.093943729 -0.287111443 -0.079436424  0.136027285 -1.696735164
 [6]  0.019698199 -1.567313000 -2.688524893 -0.295808952  1.141734152
[11] -1.870578711 -2.469392299 -1.654065363 -2.007162821 -2.007562814
[16] -2.158431063  2.334681175 -1.283306001 -0.613265761  1.084450000
[21] -0.027780662  0.184744185  1.797212143 -1.881269777 -1.912692805
[26] -1.754571636  0.250855179 -2.081122985 -1.490797599 -1.055423594
[31]  0.953784351  4.301158119  0.427808193  1.337045780  1.368568828
[36]  1.909049886 -1.279152870 -0.756001885  1.158311135  3.009217321
[41]  2.649313393  0.710148059 -0.191408976 -0.303598937  0.315883647
[46]  1.087157784  2.962968321  0.504744303  0.006299391
attr(,"internals")
              Gi      E(Gi)        V(Gi)        Z(Gi) Pr(z > E(Gi))
 [1,] 0.02715083 0.02083333 3.335030e-05  1.093943729  1.369898e-01
 [2,] 0.01950015 0.02083333 2.156145e-05 -0.287111443  6.129865e-01
 [3,] 0.02051521 0.02083333 1.603788e-05 -0.079436424  5.316573e-01
 [4,] 0.02135213 0.02083333 1.454604e-05  0.136027285  4.458999e-01
 [5,] 0.01629408 0.02083333 7.157173e-06 -1.696735164  9.551266e-01
 [6,] 0.02094741 0.02083333 3.353567e-05  0.019698199  4.921421e-01
 [7,] 0.01469141 0.02083333 1.535669e-05 -1.567313000  9.414792e-01
 [8,] 0.01232045 0.02083333 1.002595e-05 -2.688524893  9.964116e-01
 [9,] 0.02003481 0.02083333 7.287121e-06 -0.295808952  6.163120e-01
[10,] 0.02539289 0.02083333 1.594835e-05  1.141734152  1.267823e-01
[11,] 0.01438790 0.02083333 1.187277e-05 -1.870578711  9.692983e-01
[12,] 0.01305322 0.02083333 9.926387e-06 -2.469392299  9.932329e-01
[13,] 0.01431516 0.02083333 1.552912e-05 -1.654065363  9.509429e-01
[14,] 0.01451186 0.02083333 9.919078e-06 -2.007162821  9.776338e-01
[15,] 0.01451310 0.02083333 9.911239e-06 -2.007562814  9.776551e-01
[16,] 0.01515704 0.02083333 6.915963e-06 -2.158431063  9.845528e-01
[17,] 0.03158915 0.02083333 2.122420e-05  2.334681175  9.780041e-03
[18,] 0.01571282 0.02083333 1.592083e-05 -1.283306001  9.003076e-01
[19,] 0.01798811 0.02083333 2.152465e-05 -0.613265761  7.301497e-01
[20,] 0.02319539 0.02083333 4.744170e-06  1.084450000  1.390827e-01
[21,] 0.02070489 0.02083333 2.137696e-05 -0.027780662  5.110815e-01
[22,] 0.02141903 0.02083333 1.005073e-05  0.184744185  4.267148e-01
[23,] 0.02918058 0.02083333 2.157186e-05  1.797212143  3.615097e-02
[24,] 0.01534501 0.02083333 8.510938e-06 -1.881269777  9.700324e-01
[25,] 0.01577932 0.02083333 6.982049e-06 -1.912692805  9.721063e-01
[26,] 0.01535909 0.02083333 9.734335e-06 -1.754571636  9.603337e-01
[27,] 0.02182719 0.02083333 1.569656e-05  0.250855179  4.009630e-01
[28,] 0.01573131 0.02083333 6.010200e-06 -2.081122985  9.812887e-01
[29,] 0.01656483 0.02083333 8.198130e-06 -1.490797599  9.319927e-01
[30,] 0.01710439 0.02083333 1.248297e-05 -1.055423594  8.543842e-01
[31,] 0.02529414 0.02083333 2.187389e-05  0.953784351  1.700964e-01
[32,] 0.03802905 0.02083333 1.598342e-05  4.301158119  8.495388e-06
[33,] 0.02252201 0.02083333 1.558098e-05  0.427808193  3.343954e-01
[34,] 0.02618381 0.02083333 1.601375e-05  1.337045780  9.060385e-02
[35,] 0.02481457 0.02083333 8.462578e-06  1.368568828  8.556705e-02
[36,] 0.02758385 0.02083333 1.250374e-05  1.909049886  2.812783e-02
[37,] 0.01674656 0.02083333 1.020741e-05 -1.279152870  8.995784e-01
[38,] 0.01844148 0.02083333 1.000974e-05 -0.756001885  7.751760e-01
[39,] 0.02624454 0.02083333 2.182419e-05  1.158311135  1.233685e-01
[40,] 0.03083083 0.02083333 1.103763e-05  3.009217321  1.309608e-03
[41,] 0.03309317 0.02083333 2.141423e-05  2.649313393  4.032775e-03
[42,] 0.02482465 0.02083333 3.158890e-05  0.710148059  2.388062e-01
[43,] 0.02022382 0.02083333 1.013998e-05 -0.191408976  5.758974e-01
[44,] 0.01975819 0.02083333 1.254090e-05 -0.303598937  6.192833e-01
[45,] 0.02209603 0.02083333 1.597888e-05  0.315883647  3.760454e-01
[46,] 0.02712347 0.02083333 3.347607e-05  1.087157784  1.384835e-01
[47,] 0.03796202 0.02083333 3.341903e-05  2.962968321  1.523440e-03
[48,] 0.02283993 0.02083333 1.580440e-05  0.504744303  3.068692e-01
[49,] 0.02086275 0.02083333 2.180524e-05  0.006299391  4.974869e-01
attr(,"cluster")
 [1] High High High Low  Low  High Low  Low  High Low  Low  Low  Low  Low  Low 
[16] Low  Low  Low  Low  High Low  Low  High Low  Low  Low  Low  Low  Low  Low 
[31] High High Low  High Low  High High Low  High High High High Low  High Low 
[46] High High Low  High
Levels: Low High
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = columbus.map$INC, listw = listw_queen, alternative = "greater")
attr(,"class")
[1] "localG"
columbus.map$LocalGi <- as.numeric(Local_M)

columbus.map$Gi_cat <- cut(
  columbus.map$LocalGi,
  breaks = c(-Inf, -1.96, 1.96, Inf),
  labels = c("Coldspot", "No SAC", "Hotspot")
)

tmap_mode("plot")
ℹ tmap modes "plot" - "view"
tm_shape(columbus.map) +
  tm_polygons(
    fill = "Gi_cat",   # pakai string
    fill.scale = tm_scale_categorical(
      values = c("blue", "white", "red"),
      labels = c("Coldspot", "No SAC", "Hotspot")
    ),
    fill.legend = tm_legend(title = "Local Gi* (Getis-Ord)")
  ) +
  tm_layout(legend.outside = TRUE)
Multiple palettes called "blue" found: "kovesi.blue", "tableau.blue". The first one, "kovesi.blue", is returned.

Interpretasi Nilai Z(Gi)

  • Jika nilai z-score Z(Gi) suatu lokasi berada pada selang −Zα/2≤Z(Gi)≤Zα/2 maka lokasi tersebut bukan lokasi hotspot maupun coldspot

  • Suatu lokasi dikatakan hotspot jika nilai Z(Gi)>Zα/2 dan dikatakan sebagai coldspot jika nilai Z(Gi)<−Zα/2

Geostatistical Process

Proses geostatistik umumnya mencakup beberapa langkah utama:

  1. Analisis Data Eksploratif: Memahami distribusi dan karakteristik data, termasuk apakah data stasioner (mean dan varians tidak berubah di seluruh area) dan apakah ada tren spasial.

  2. Variografi: Menganalisis korelasi spasial dengan menggunakan variogram atau semivariogram. Langkah ini sangat penting untuk memodelkan struktur spasial data.

  3. Kriging (Interpolasi Spasial): Menggunakan model variogram yang telah dibuat untuk memprediksi nilai pada lokasi yang tidak disampel. Kriging adalah metode interpolasi geostatistik yang paling umum dan memberikan estimasi terbaik yang tidak bias. Kriging tidak hanya memberikan nilai prediksi, tetapi juga estimasi kesalahan (ketidakpastian) dari prediksi tersebut.

  4. Validasi: Memvalidasi hasil prediksi dengan membandingkan nilai prediksi dengan nilai yang disampel.

library(sp)
library(gstat)
library(spdep)
library(ncf)

#Meuse Dataset
data(meuse)
coordinates(meuse) <- ~x+y

# Inspect the first few rows of the data
head(meuse)
  cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse dist.m
1    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah     50
2     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1      Ah     30
3     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1      Ah    150
4     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0      Ga    270
5     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0      Ah    380
6     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0      Ga    470
plot(meuse)

x <- meuse@coords[,1]  # Koordinat X
y <- meuse@coords[,2]  # Koordinat Y
z <- meuse$zinc        # Nilai kandungan zinc

- Korelogram

Dalam geostatistik, korelogram adalah grafik yang menunjukkan autokorelasi spasial sebagai fungsi dari jarak antar titik data. Autokorelasi spasial adalah sejauh mana nilai-nilai data pada lokasi yang berdekatan memiliki kemiripan.

  • Sumbu x menunjukkan jarak (lag) antara pasangan titik data.

  • Sumbu y menunjukkan nilai koefisien korelasi.

Jika titik-titik data saling berdekatan (jarak kecil), koefisien korelasi akan tinggi (mendekati 1), menunjukkan kemiripan yang kuat. Seiring bertambahnya jarak, koefisien korelasi cenderung menurun dan mendekati nol. Korelogram sering digunakan sebagai alat bantu visual untuk memahami seberapa jauh pengaruh satu titik data terhadap titik data lainnya.

library(ncf)
# Menghitung correlogram dengan jarak (lag distance)
correlog <- correlog(x = x, y = y, z = z, increment = 100, resamp = 500)
50  of  500 
100  of  500 
150  of  500 
200  of  500 
250  of  500 
300  of  500 
350  of  500 
400  of  500 
450  of  500 
500  of  500 
plot(correlog)

  • Sumbu X menunjukkan kelas jarak (lag distance) antara titik-titik pengamatan dalam dataset. Jarak ini menunjukkan seberapa jauh titik-titik data yang dibandingkan satu sama lain untuk mengukur autokorelasi spasial.

  • Sumbu Y menunjukkan nilai autokorelasi spasial untuk setiap kelas jarak. Autokorelasi spasial positif berarti titik-titik dengan jarak tertentu memiliki nilai yang lebih mirip daripada yang diharapkan secara acak, sementara autokorelasi negatif berarti mereka lebih berbeda daripada yang diharapkan secara acak.

  • Titik hitam menunjukkan autokorelasi spasial yang signifikan secara statistik, sementara titik putih menunjukkan nilai autokorelasi spasial yang tidak signifikan.

  • Titik hitam positif di jarak dekat: Ada kecenderungan lokasi yang berdekatan memiliki nilai yang serupa (kluster/pengelompokan).

  • Titik hitam negatif di jarak menengah: Ada pola kontras (lokasi nilai tinggi dikelilingi lokasi nilai rendah, atau sebaliknya).

  • Titik putih di jarak jauh: Ketergantungan spasial sudah sepenuhnya hilang. Lokasi-lokasi tersebut sudah saling bebas (independent).

- Variogram dan Semivariogram

Apa itu Variogram dan Semivariogram?

Variogram dan semivariogram adalah alat fundamental dalam geostatistik untuk mengukur dan memvisualisasikan bagaimana variasi data berhubungan dengan jarak.

  • Secara Teknis: Semivariogram (γ(h)) dihitung sebagai setengah dari rata-rata kuadrat perbedaan antara nilai-nilai data pada dua lokasi yang dipisahkan oleh jarak (h). Variogram (2γ(h)) adalah dua kali semivariogram. Meskipun berbeda secara definisi matematis, dalam praktiknya, istilah “variogram” sering digunakan secara bergantian untuk merujuk pada plot semivariogram.

  • Tujuan Utama: Keduanya digunakan untuk mengukur dan memodelkan variabilitas spasial. Mereka menunjukkan seberapa besar perbedaan nilai data seiring dengan bertambahnya jarak antar titik. Ini membantu kita memahami apakah data Anda memiliki struktur spasial atau tidak.

- Komponen-komponen Penting pada Grafik

Grafik variogram atau semivariogram memiliki tiga komponen kunci yang sangat penting untuk interpretasi:

  1. Nugget Effect: Ini adalah nilai variabilitas pada jarak nol. Nugget mewakili variasi yang tidak dapat dijelaskan oleh korelasi spasial, seperti kesalahan pengukuran atau fluktuasi yang terjadi pada skala yang sangat kecil. Jika nilai nugget nol, berarti tidak ada variasi pada jarak yang sangat dekat.

  2. Range: Ini adalah jarak maksimum di mana titik-titik data masih memiliki korelasi spasial satu sama lain. Di luar jarak ini, nilai-nilai data dianggap tidak lagi saling terkait. Range sangat penting karena menentukan sejauh mana suatu titik data dapat digunakan untuk memprediksi nilai pada lokasi lain.

  3. Sill: Ini adalah total variabilitas dalam data. Sill adalah nilai variabilitas yang dicapai grafik saat mendatar, menunjukkan titik di mana korelasi spasial sudah tidak ada lagi.

Dengan memodelkan variogram/semivariogram, kita dapat memahami seberapa jauh pengaruh satu titik data terhadap titik lainnya, yang merupakan dasar untuk metode interpolasi geostatistik seperti kriging.

library(gstat)
# Membuat variogram untuk kandungan zinc di dataset meuse
vg <- variogram(log(zinc) ~ 1, meuse)

# Plot variogram
plot(vg)

  • variogram() menghitung semivariance (variabilitas spasial) berdasarkan jarak antar titik. Semivariance rendah pada jarak pendek menunjukkan bahwa nilai di titik yang berdekatan cenderung mirip, sementara semivariance yang lebih tinggi menunjukkan peningkatan variabilitas pada jarak yang lebih jauh.

  • Plot variogram menunjukkan semivariance sebagai fungsi jarak. Semivariance yang meningkat dengan jarak menunjukkan adanya struktur spasial.

  • Jika semivariance meningkat dengan jarak dan akhirnya mencapai nilai konstan (sill), ini menunjukkan adanya autokorelasi spasial positif pada jarak pendek, yang kemudian berkurang pada jarak yang lebih jauh (range).

Semivariogram adalah fungsi yang mirip dengan variogram, tetapi digunakan untuk memodelkanvariabilitas spasial dan memperkirakan hubungan spasial untuk prediksi lebih lanjut seperti kriging.

- Variogram Eksperimental vs. Fitted Variogram

  • Variogram Eksperimental: Dihitung langsung dari data Anda, menghasilkan plot titik-titik yang menunjukkan variabilitas pada jarak-jarak tertentu.

  • Fitted Variogram: Ini adalah model matematika (kurva halus) yang dibuat untuk menyesuaikan pola variogram eksperimental. Model inilah yang akhirnya digunakan dalam proses prediksi spasial seperti kriging.

# Fit a model to the semivariogram
semivg_model <- fit.variogram(vg, model = vgm(psill = 0.5, "Sph", range = 1000, nugget = 0.1))

# Plot semivariogram dengan model yang di-fit
plot(vg, model = semivg_model)

library(automap)
# Autofit variogram untuk data log-zinc
afvg <- autofitVariogram(log(zinc) ~ 1, meuse)

# Plot hasil
plot(afvg)

  • fit.variogram() digunakan untuk memodelkan variogram dengan model semivariogram tertentu. Dalam hal ini, kita menggunakan model Spherical.

  • Pada jarak pendek (di bawah 500 meter), nilai semivariance masih rendah, yang menunjukkan bahwa lokasi yang berdekatan memiliki nilai yang lebih mirip satu sama lain. Ini adalah indikasi adanya autokorelasi spasial positif. Titik-titik pengamatan yang berdekatan memiliki hubungan spasial.

  • Saat jarak meningkat (500 hingga sekitar 900 meter), semivariance juga meningkat. Ini menunjukkan bahwa semakin jauh jarak antar titik pengamatan, semakin besar perbedaan nilai di antara mereka, yang menunjukkan bahwa autokorelasi spasial menurun.

  • Angka yang ditampilkan di dekat titik semivariogram (misalnya 1314 atau 1349) menunjukkan jumlah pasangan titik (point pairs) yang digunakan untuk menghitung nilai semivariance pada kelas jarak tertentu. Semakin banyak jumlah pasangan titik, semakin stabil dan lebih dapat dipercaya nilai semivariance pada jarak tersebut. Sebaliknya, jika jumlah pasangan sedikit, nilai semivariance bisa lebih fluktuatif atau kurang representatif.