Input Data

library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(sp)
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
library(ggplot2)
library(readxl)
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(terra)
## terra 1.7.78
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
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
petaju <- st_read("C:\\Users\\luthf\\Downloads\\DKI Jakarta\\commondata\\data1\\Indo_Districts.shp")
## Reading layer `Indo_Districts' from data source 
##   `C:\Users\luthf\Downloads\DKI Jakarta\commondata\data1\Indo_Districts.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 440 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.05965 ymin: -10.99741 xmax: 141.0072 ymax: 5.906884
## Geodetic CRS:  WGS 84
petaju <- st_make_valid(petaju)
peta<-petaju[c(113:117), ]
peta
## Simple feature collection with 5 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.6816 ymin: -6.376739 xmax: 106.9731 ymax: -6.088497
## Geodetic CRS:  WGS 84
##                 NAMA_KAB   NAMA_PROP KODE_KAB KODE_PROP        PROV
## 113 Kota Jakarta Selatan Dki Jakarta     3171        31 Dki Jakarta
## 114   Kota Jakarta Pusat Dki Jakarta     3173        31 Dki Jakarta
## 115   Kota Jakarta Timur Dki Jakarta     3172        31 Dki Jakarta
## 116   Kota Jakarta Barat Dki Jakarta     3174        31 Dki Jakarta
## 117   Kota Jakarta Utara Dki Jakarta     3175        31 Dki Jakarta
##                      KAB PROV_CODE KAB_CODE                       geometry
## 113 Kota Jakarta Selatan        31     3171 MULTIPOLYGON (((106.8584 -6...
## 114   Kota Jakarta Pusat        31     3173 MULTIPOLYGON (((106.8667 -6...
## 115   Kota Jakarta Timur        31     3172 MULTIPOLYGON (((106.8983 -6...
## 116   Kota Jakarta Barat        31     3174 MULTIPOLYGON (((106.7195 -6...
## 117   Kota Jakarta Utara        31     3175 MULTIPOLYGON (((106.9276 -6...
data <- readxl::read_excel("C:/Users/luthf/Downloads/Data Tugas Regspas Moran.xlsx")
str(data)
## tibble [6 Γ— 2] (S3: tbl_df/tbl/data.frame)
##  $ KODE_KAB: num [1:6] 3172 3171 3173 3174 3175 ...
##  $ Suara PG: num [1:6] 73181 49925 21463 44252 44179 ...
data
## # A tibble: 6 Γ— 2
##   KODE_KAB `Suara PG`
##      <dbl>      <dbl>
## 1     3172      73181
## 2     3171      49925
## 3     3173      21463
## 4     3174      44252
## 5     3175      44179
## 6     3101       2260
Data_merge<- merge(data, peta, by= "KODE_KAB")
str(Data_merge)
## 'data.frame':    5 obs. of  10 variables:
##  $ KODE_KAB : num  3171 3172 3173 3174 3175
##  $ Suara PG : num  49925 73181 21463 44252 44179
##  $ NAMA_KAB : chr  "Kota Jakarta Selatan" "Kota Jakarta Timur" "Kota Jakarta Pusat" "Kota Jakarta Barat" ...
##  $ NAMA_PROP: chr  "Dki Jakarta" "Dki Jakarta" "Dki Jakarta" "Dki Jakarta" ...
##  $ KODE_PROP: chr  "31" "31" "31" "31" ...
##  $ PROV     : chr  "Dki Jakarta" "Dki Jakarta" "Dki Jakarta" "Dki Jakarta" ...
##  $ KAB      : chr  "Kota Jakarta Selatan" "Kota Jakarta Timur" "Kota Jakarta Pusat" "Kota Jakarta Barat" ...
##  $ PROV_CODE: chr  "31" "31" "31" "31" ...
##  $ KAB_CODE : chr  "3171" "3172" "3173" "3174" ...
##  $ geometry :sfc_MULTIPOLYGON of length 5; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:400, 1:2] 107 107 107 107 107 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
# konversi ke sf
Data_merge <- st_sf(Data_merge)
Data_merge
## Simple feature collection with 5 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.6816 ymin: -6.376739 xmax: 106.9731 ymax: -6.088497
## Geodetic CRS:  WGS 84
##   KODE_KAB Suara PG             NAMA_KAB   NAMA_PROP KODE_PROP        PROV
## 1     3171    49925 Kota Jakarta Selatan Dki Jakarta        31 Dki Jakarta
## 2     3172    73181   Kota Jakarta Timur Dki Jakarta        31 Dki Jakarta
## 3     3173    21463   Kota Jakarta Pusat Dki Jakarta        31 Dki Jakarta
## 4     3174    44252   Kota Jakarta Barat Dki Jakarta        31 Dki Jakarta
## 5     3175    44179   Kota Jakarta Utara Dki Jakarta        31 Dki Jakarta
##                    KAB PROV_CODE KAB_CODE                       geometry
## 1 Kota Jakarta Selatan        31     3171 MULTIPOLYGON (((106.8584 -6...
## 2   Kota Jakarta Timur        31     3172 MULTIPOLYGON (((106.8983 -6...
## 3   Kota Jakarta Pusat        31     3173 MULTIPOLYGON (((106.8667 -6...
## 4   Kota Jakarta Barat        31     3174 MULTIPOLYGON (((106.7195 -6...
## 5   Kota Jakarta Utara        31     3175 MULTIPOLYGON (((106.9276 -6...
tm_shape(Data_merge) + tm_polygons("Suara PG") + tm_layout(title = "Map of Suara PG")

Matriks Pembobot Spasial

Queen Cognity

queen <- poly2nb(Data_merge, queen = TRUE)
summary(queen)
## Neighbour list object:
## Number of regions: 5 
## Number of nonzero links: 16 
## Percentage nonzero weights: 64 
## Average number of links: 3.2 
## Link number distribution:
## 
## 3 4 
## 4 1 
## 4 least connected regions:
## 1 2 4 5 with 3 links
## 1 most connected region:
## 3 with 4 links
W.queen <- nb2listw(queen, style='W',zero.policy=TRUE)
W.queen2 <- nb2listw(queen, style='B')
plot(st_geometry(Data_merge), border="blue",col='gray')
coords <- st_coordinates(st_centroid(st_geometry(Data_merge)))
plot(queen, coords, add = TRUE, col = "red")

Uji Normalitas

shapiro.test(Data_merge$"Suara PG")
## 
##  Shapiro-Wilk normality test
## 
## data:  Data_merge$"Suara PG"
## W = 0.93973, p-value = 0.664

\(H_0\): Data berasal dari distribusi normal.

\(H_1\): Data tidak berasal dari distribusi normal.

p-value>0.05, artinya data berdistribusi normal.

Plot QQ

qqnorm(Data_merge$"Suara PG")
qqline(Data_merge$"Suara PG", col = "red")

Autokorelasi Global

Indeks Moran Global

moran(Data_merge$"Suara PG", W.queen, n=length(W.queen$neighbours), S0=Szero(W.queen))
## $I
## [1] -0.2601213
## 
## $K
## [1] 2.426301

Nilai kurtosis K = 2.426 menunjukkan bahwa distribusi data memiliki bentuk yang sedikit lebih datar dibandingkan dengan distribusi normal, dengan lebih sedikit data yang berada di ekor distribusi (nilai ekstrem).

IM <- moran.test(Data_merge$"Suara PG", W.queen,
alternative="two.sided", zero.policy = FALSE)
IM
## 
##  Moran I test under randomisation
## 
## data:  Data_merge$"Suara PG"  
## weights: W.queen    
## 
## Moran I statistic standard deviate = -0.076793, p-value = 0.9388
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.26012135       -0.25000000        0.01737113
moran.mc(Data_merge$"Suara PG", listw = W.queen, zero.policy = FALSE, nsim = 99)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  Data_merge$"Suara PG" 
## weights: W.queen  
## number of simulations + 1: 100 
## 
## statistic = -0.26012, observed rank = 64.5, p-value = 0.355
## alternative hypothesis: greater

Moran’s I = -0.2601 menunjukkan kecenderungan spasial negatif, di mana lokasi-lokasi yang berdekatan cenderung memiliki nilai yang berbeda (penyebaran spasial).

Uji Moran dengan alternative two.sided menghasilkan nilai p-value sebesar 0.9388, lebih besar dari taraf nyata 5%. Sehingga dapat disimpulkan bahwa tidak terdapat bukti yang signifikan untuk menunjukkan adanya autocorrelation spasial dalam data β€˜Suara PG’.

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

GC <- geary.test(Data_merge$"Suara PG", W.queen, randomisation = T)
GC
## 
##  Geary C test under randomisation
## 
## data:  Data_merge$"Suara PG" 
## weights: W.queen   
## 
## Geary C statistic standard deviate = -0.42845, p-value = 0.6658
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.05215037        1.00000000        0.01481534

Hasil dari uji Global Geary’s C menunjukkan nilai c sebesar 1.0521, yang nilainya berada di sekitar 1. Hal ini mengindikasikan bahwa tidak terdapat autokorelasi(pola acak).

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(Data_merge$"Suara PG", W.queen2, alternative = "greater")
GO
## 
##  Getis-Ord global G statistic
## 
## data:  Data_merge$"Suara PG" 
## weights: W.queen2   
## 
## standard deviate = -1.5314, p-value = 0.9372
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##         0.74119341         0.80000000         0.00147466

Hasil Getis-Ord Global G sebesar 0.741 yang mendekati 0 Menunjukkan tidak adanya pola yang signifikan atau distribusi acak dari nilai tinggi dan rendah.

Autokorelasi Lokal

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).

LM <- localmoran(Data_merge$"Suara PG", W.queen, alternative = "greater")
head(LM)
##             Ii         E.Ii      Var.Ii        Z.Ii Pr(z > E(Ii))
## 1 -0.003681282 -0.010155086 0.005584422  0.08663041     0.4654826
## 2 -0.788891964 -0.648997358 0.126555437 -0.39324261     0.6529298
## 3 -0.580399712 -0.580399712 0.000000000         NaN           NaN
## 4  0.069685803 -0.005064032 0.002799104  1.41286507     0.0788477
## 5  0.002680416 -0.005383812 0.002974904  0.14785171     0.4412299
Data_merge$lmI <- LM[, "Ii"] # local Moran's I
Data_merge$lmZ <- LM[, "Z.Ii"] # z-scores
# p-values corresponding to alternative greater
Data_merge$lmp <- LM[, "Pr(z > E(Ii))"]

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.

p1 <- tm_shape(Data_merge) +
  tm_polygons(col = "lmI", title = "Local Moran's I",
              style = "quantile") +
  tm_layout(legend.outside = TRUE)
p1
## Variable(s) "lmI" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

\(H_0\): Tidak Terdapat autokorelasi spasial

\(H_1\): Terdapat autokorelasi spasial

p2 <- tm_shape(Data_merge) +
  tm_polygons(col = "lmp", title = "p-value",
              breaks = c(-Inf, 0.05, Inf)) +
  tm_layout(legend.outside = TRUE)
p2

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.

tm_shape(Data_merge) + 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)

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

# Autokorelasi spasial positif (Z-score > 1.96)
positive_autocorr <- Data_merge[Data_merge$lmZ > 1.96, ]
# Autokorelasi spasial negatif (Z-score < -1.96)
negative_autocorr <- Data_merge[Data_merge$lmZ < -1.96, ]
# Area tanpa autokorelasi spasial (Z-score antara -1.96 dan 1.96)
no_autocorr <-Data_merge[Data_merge$lmZ >= -1.96 & Data_merge$lmZ <= 1.96, ]
moran.plot(Data_merge$`Suara PG`, W.queen)

Moran plot diatas memiliki memiliki kemiringan negatif, ini menunjukkan autokorelasi spasial negatif, di mana lokasi dengan nilai tinggi dikelilingi oleh lokasi dengan nilai rendah (dan sebaliknya). Plot tersebut terlihat bahwa amatan cenderung menyebar di Kuadran II (Low-High) dan Kuadran IV (High-Low). Kuadran II (Low-High) merupakan lokasi dengan nilai rendah yang dikelilingi oleh lokasi dengan nilai tinggi (indikasi outlier) dan Kuadran IV (High-Low) merupakan lokasi dengan nilai tinggi yang dikelilingi oleh lokasi dengan nilai rendah (indikasi outlier).

Getis-Ord Gi

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

Local_M <- localG(Data_merge$"Suara PG", W.queen, alternative = "greater")
Local_M
## [1]  0.08663041 -0.39324261         NaN -1.41286507 -0.14785171
## attr(,"internals")
##             Gi E(Gi)        V(Gi)       Z(Gi) Pr(z > E(Gi))
## [1,] 0.2528945  0.25 0.0011163909  0.08663041     0.4654826
## [2,] 0.2410373  0.25 0.0005194704 -0.39324261     0.6529298
## [3,] 0.2500000  0.25 0.0000000000         NaN           NaN
## [4,] 0.2040940  0.25 0.0010556930 -1.41286507     0.9211523
## [5,] 0.2451987  0.25 0.0010545378 -0.14785171     0.5587701
## attr(,"cluster")
## [1] High High Low  Low  Low 
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = Data_merge$"Suara PG", listw = W.queen, alternative = "greater")
## attr(,"class")
## [1] "localG"
Data_merge$lmZ2 <- Local_M
tm_shape(Data_merge) + tm_polygons(col = "lmZ2",
title = "Local Gi", style = "fixed",
breaks = c(-Inf, -1.96, 1.96, Inf),
labels = c("Coldspot", "No SAC", "Hostpot"),
palette =  c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)

Hasil analisis menunjukkan bahwa lokasi-lokasi tersebut bukan lokasi hotspot maupun coldspot.

Kesimpulan

Berdasarkan analisis Indeks Moran Global dan Indeks Moran Lokal terkait Suara Partai Golkar di DKI Jakarta, dapat disimpulkan:

  1. Indeks Moran Global menunjukkan bahwa di tingkat Kabupaten/Kota di DKI Jakarta, tidak ditemukan adanya autokorelasi global. Artinya, suara di satu wilayah tidak memiliki keterkaitan secara umum dengan wilayah sekitarnya.

  2. Indeks Moran Lokal menunjukkan bahwa di tingkat Kabupaten/Kota di DKI Jakarta tidak ditemukan adanya autokorelasi lokal. Ini berarti bahwa pola distribusi suara Partai Golkar di setiap wilayah tidak memiliki keterkaitan signifikan dengan wilayah sekitarnya. Dengan kata lain, suara di satu wilayah tidak dipengaruhi atau tidak berkorelasi dengan suara di wilayah-wilayah terdekat, baik dalam bentuk hubungan positif maupun negatif.

  3. Hasil analisis menunjukkan bahwa di tingkat Kabupaten/Kota di DKI Jakarta, wilayah-wilayah tersebut tidak termasuk dalam wilayah hotspot maupun coldspot.