1 Library

packages <- c("sf", "spData", "spdep", "tmap", "ggplot2", "dplyr",
              "cartogram", "classInt", "RColorBrewer")

installed <- packages %in% rownames(installed.packages())
if (any(!installed)) install.packages(packages[!installed])

library(sf)
library(spData)
library(spdep)
library(tmap)
library(ggplot2)
library(dplyr)
library(cartogram)
library(classInt)

2 Import Data

nc_path <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(nc_path, quiet = TRUE)

nc <- nc %>%
  mutate(SID74_rate = (SID74 / BIR74) * 1000)

summary(nc$SID74_rate)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.084   1.855   2.046   2.604   9.554

3 Latihan 1 — Perbandingan Metode Pembagian Kelas

style_list <- c("equal", "quantile", "jenks", "sd")

maps <- lapply(style_list, function(s) {
  tm_shape(nc) +
    tm_polygons("SID74_rate",
                palette = "YlOrRd",
                style = s,
                n = 5,
                title = s) +
    tm_layout(main.title = paste("style =", s),
              main.title.size = 0.9,
              legend.outside = FALSE,
              legend.text.size = 0.5)
})

tmap_arrange(maps[[1]], maps[[2]], maps[[3]], maps[[4]], ncol = 2)

hist(nc$SID74_rate, breaks = 15, main = "Distribusi SID74_rate",
     xlab = "SID74 rate (per 1000 kelahiran)", col = "lightblue")

Jawaban:

3.0.1 1. Bentuk distribusi SID74_rate

Berdasarkan histogram, sebagian besar county memiliki nilai SID74_rate yang berada di sekitar 0–3 kasus per 1.000 kelahiran. Sementara itu, hanya beberapa county yang memiliki nilai cukup tinggi, yaitu sekitar 8–9,5. Hal ini menunjukkan bahwa distribusi SID74_rate menceng ke kanan (right-skewed). Kondisi ini juga dapat dilihat dari nilai rata-rata sebesar 2,046 yang lebih besar dibandingkan nilai median sebesar 1,855. Perbedaan tersebut menunjukkan bahwa terdapat beberapa county dengan nilai SID74_rate yang tinggi sehingga distribusi tertarik ke arah kanan.

3.0.2 2. Apakah keempat metode menghasilkan pola warna yang sama?

Tidak. Keempat metode menghasilkan pola warna yang berbeda karena masing-masing menggunakan cara yang berbeda dalam menentukan batas kelas, meskipun jumlah kelas yang digunakan sama, yaitu n = 5.

  • Equal membagi rentang nilai menjadi 5 interval dengan lebar yang sama. Karena data menceng ke kanan, sebagian besar county akan masuk ke kelas dengan nilai rendah, sedangkan kelas dengan nilai tinggi hanya berisi sedikit county.
  • Quantile membagi data berdasarkan jumlah county, sehingga jumlah county pada setiap kelas relatif sama. Oleh karena itu, warna pada peta cenderung lebih merata.
  • Jenks (Natural Breaks) menentukan batas kelas berdasarkan pengelompokan alami yang terdapat pada data. County dengan nilai yang berdekatan akan cenderung masuk ke kelas yang sama.
  • SD (Standard Deviation) menentukan kelas berdasarkan jarak nilai terhadap rata-rata dalam satuan simpangan baku. County yang memiliki nilai jauh dari rata-rata akan masuk ke kelas yang lebih ekstrem.

Jadi, karena dasar pembagian kelasnya berbeda, pola warna yang dihasilkan oleh keempat metode juga tidak sama.

3.0.3 3. Metode mana yang membuat sebagian besar county berada pada kelas warna yang sama?

Metode equal cenderung membuat sebagian besar county berada pada kelas dengan nilai rendah. Hal ini disebabkan oleh rentang data yang dibagi menjadi interval dengan lebar yang sama, sedangkan sebagian besar nilai SID74_rate memang berada pada nilai yang rendah. Akibatnya, beberapa kelas dengan nilai tinggi hanya diisi oleh sedikit county.

3.0.4 4. Metode mana yang paling sesuai untuk data ini?

Untuk data SID74_rate yang memiliki distribusi menceng ke kanan, metode quantile atau Jenks lebih sesuai dibandingkan metode equal dan SD. Metode quantile membuat jumlah county pada setiap kelas relatif seimbang sehingga perbedaan antarwilayah lebih mudah terlihat pada peta. Sementara itu, metode Jenks dapat digunakan karena batas kelasnya mengikuti pola pengelompokan alami pada data.

Metode equal kurang sesuai karena sebagian besar county akan terkumpul pada kelas rendah, sedangkan metode SD juga kurang sesuai karena distribusi data tidak simetris dan menceng ke kanan.

4 Latihan 2 — Cartogram Berdasarkan SID74

nc_proj <- st_transform(nc, crs = 32119)

# Cartogram berdasarkan jumlah kelahiran (BIR74) 
nc_cartogram_bir <- cartogram_cont(nc_proj, weight = "BIR74", itermax = 5)

# Cartogram berdasarkan jumlah kasus SIDS (SID74) 
nc_cartogram_sid <- cartogram_cont(nc_proj, weight = "SID74", itermax = 5)

p_bir <- tm_shape(nc_cartogram_bir) +
  tm_polygons("BIR74", palette = "Blues", style = "quantile", title = "BIR74") +
  tm_layout(main.title = "Cartogram: BIR74", main.title.size = 0.9)

p_sid <- tm_shape(nc_cartogram_sid) +
  tm_polygons("SID74", palette = "Reds", style = "quantile", title = "SID74") +
  tm_layout(main.title = "Cartogram: SID74", main.title.size = 0.9)

tmap_arrange(p_bir, p_sid, ncol = 2)

Untuk mengidentifikasi county mana yang paling menonjol/berubah ukurannya, kita bisa membandingkan langsung nilai BIR74 dan SID74 per county:

nc_proj %>%
  st_drop_geometry() %>%
  select(NAME, BIR74, SID74) %>%
  arrange(desc(SID74)) %>%
  head(10)
##           NAME BIR74 SID74
## 1  Mecklenburg 21588    44
## 2   Cumberland 20366    38
## 3      Robeson  7889    31
## 4       Onslow 11158    29
## 5     Guilford 16184    23
## 6      Halifax  3608    18
## 7        Wayne  6638    18
## 8   Rockingham  4449    16
## 9       Durham  7970    16
## 10        Wake 14484    16
nc_proj %>%
  st_drop_geometry() %>%
  select(NAME, BIR74, SID74) %>%
  arrange(desc(BIR74)) %>%
  head(10)
##           NAME BIR74 SID74
## 1  Mecklenburg 21588    44
## 2   Cumberland 20366    38
## 3     Guilford 16184    23
## 4         Wake 14484    16
## 5      Forsyth 11858    10
## 6       Onslow 11158    29
## 7       Gaston  9014    11
## 8       Durham  7970    16
## 9      Robeson  7889    31
## 10    Buncombe  7515     9

Jawaban:

4.0.1 1. Apakah county yang menonjol pada kedua cartogram sama?

Tidak selalu sama. County dengan nilai BIR74 tinggi umumnya merupakan county dengan jumlah kelahiran yang besar, sehingga biasanya county yang penduduknya juga relatif banyak. Sementara itu, county dengan nilai SID74 tinggi belum tentu memiliki jumlah kelahiran yang besar. Ada county dengan jumlah kelahiran sedang tetapi memiliki kasus SIDS yang cukup tinggi, sehingga tetap dapat terlihat menonjol pada cartogram SID74.

Untuk melihat perbedaannya, county yang berada pada peringkat tertinggi BIR74 dan SID74 dapat dibandingkan secara langsung.

4.0.2 2. County mana yang berubah ukuran paling jelas ketika menggunakan SID74?

County yang memiliki nilai SID74 paling tinggi akan mengalami perubahan ukuran yang paling jelas pada cartogram SID74. Hal ini karena ukuran wilayah pada cartogram dibuat mengikuti nilai peubah yang digunakan. Jadi, semakin besar nilai SID74, semakin besar pula wilayah county tersebut pada cartogram.

Berdasarkan hasil pengurutan arrange(desc(SID74)), county yang berada pada peringkat teratas merupakan county yang paling menonjol pada cartogram SID74.

4.0.3 3. Mengapa cartogram BIR74 dan SID74 dapat memberikan pola yang berbeda?

Cartogram BIR74 dan SID74 menggunakan peubah yang berbeda sehingga ukuran county yang ditampilkan juga dapat berbeda. BIR74 menunjukkan jumlah kelahiran pada suatu county, sehingga lebih menggambarkan besar kecilnya populasi atau aktivitas demografis di wilayah tersebut. Sementara itu, SID74 menunjukkan jumlah kasus SIDS, yang tidak hanya dipengaruhi oleh jumlah kelahiran tetapi juga oleh tingkat kejadian SIDS (SID74_rate).

Dengan demikian, county yang memiliki jumlah kelahiran besar belum tentu memiliki jumlah kasus SIDS paling tinggi. Sebaliknya, county dengan jumlah kelahiran yang lebih sedikit dapat memiliki SID74 yang cukup tinggi jika tingkat kejadian SIDS di county tersebut tinggi. Oleh karena itu, cartogram BIR74 dan SID74 dapat menonjolkan county yang berbeda.

5 Latihan 3 — Efek Zonasi (MAUP)

# 006 = 3 digit terakhir NIM saya
set.seed(006)

n_regions <- 10

nc$region_acak_baru <- sample(
  1:n_regions,
  nrow(nc),
  replace = TRUE
)

nc_region_baru <- nc %>%
  group_by(region_acak_baru) %>%
  summarise(
    SID74_total = sum(SID74),
    BIR74_total = sum(BIR74)
  ) %>%
  mutate(
    SID74_rate_region = (SID74_total / BIR74_total) * 1000
  )

# Bandingkan dengan pengelompokan seed = 123 dari modul
set.seed(123)
nc$region_acak_123 <- sample(1:n_regions, nrow(nc), replace = TRUE)
nc_region_123 <- nc %>%
  group_by(region_acak_123) %>%
  summarise(
    SID74_total = sum(SID74),
    BIR74_total = sum(BIR74)
  ) %>%
  mutate(SID74_rate_region = (SID74_total / BIR74_total) * 1000)

p_seed123 <- tm_shape(nc_region_123) +
  tm_polygons("SID74_rate_region", palette = "YlOrRd", style = "quantile",
              title = "Rate (seed=123)") +
  tm_layout(main.title = "Zonasi: seed = 123", main.title.size = 0.9)

p_seedbaru <- tm_shape(nc_region_baru) +
  tm_polygons("SID74_rate_region", palette = "YlOrRd", style = "quantile",
              title = "Rate (seed baru)") +
  tm_layout(main.title = "Zonasi: seed = NIM", main.title.size = 0.9)

tmap_arrange(p_seed123, p_seedbaru, ncol = 2)

summary(nc_region_123$SID74_rate_region)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.694   1.844   2.058   2.081   2.310   2.608
summary(nc_region_baru$SID74_rate_region)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.260   1.471   2.130   1.918   2.226   2.615

Jawaban:

Meskipun kedua percobaan sama-sama membentuk 10 region, pembagian county pada setiap region berbeda karena menggunakan seed yang berbeda. Hal ini dapat dilihat dari pola warna pada kedua peta yang juga berubah. Pada seed = 123, posisi county dengan rate tertentu berbeda dengan hasil pada seed = NIM.

Perbedaan tersebut terjadi karena proses pembentukan region dilakukan secara acak. Ketika seed diubah, county yang tergabung dalam masing-masing region juga dapat berubah. Akibatnya, nilai rate pada setiap region ikut berubah sehingga pola warna pada peta menjadi berbeda.

Rentang nilai rate pada kedua hasil juga berbeda. Pada seed = 123, rate berada pada rentang sekitar 1,694–2,608, sedangkan pada seed = NIM berada pada rentang sekitar 1,260–2,615. Jadi, meskipun jumlah region yang dibentuk tetap sama, perubahan seed dapat menghasilkan pembagian wilayah dan nilai rate yang berbeda.

Dengan demikian, hasil zonasi tidak hanya dipengaruhi oleh jumlah region yang digunakan, tetapi juga oleh bagaimana county dikelompokkan ke dalam setiap region.

6 Latihan 4 — K-Nearest Neighbors dengan k Berbeda

coords <- st_coordinates(st_centroid(st_geometry(nc)))

knn_k2 <- knearneigh(coords, k = 2) |> knn2nb()
knn_k4 <- knearneigh(coords, k = 4) |> knn2nb()
knn_k8 <- knearneigh(coords, k = 8) |> knn2nb()

knn_k2
## Neighbour list object:
## Number of regions: 100 
## Number of nonzero links: 200 
## Percentage nonzero weights: 2 
## Average number of links: 2 
## Non-symmetric neighbours list
knn_k4
## Neighbour list object:
## Number of regions: 100 
## Number of nonzero links: 400 
## Percentage nonzero weights: 4 
## Average number of links: 4 
## Non-symmetric neighbours list
knn_k8
## Neighbour list object:
## Number of regions: 100 
## Number of nonzero links: 800 
## Percentage nonzero weights: 8 
## Average number of links: 8 
## Non-symmetric neighbours list
par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
plot(st_geometry(nc), border = "grey60", main = "KNN k = 2")
plot(knn_k2, coords, add = TRUE, col = "darkgreen", lwd = 1.2)

plot(st_geometry(nc), border = "grey60", main = "KNN k = 4")
plot(knn_k4, coords, add = TRUE, col = "blue", lwd = 1.2)

plot(st_geometry(nc), border = "grey60", main = "KNN k = 8")
plot(knn_k8, coords, add = TRUE, col = "purple", lwd = 1.2)

par(mfrow = c(1, 1))

Jawaban:

6.0.1 1. Apa yang terjadi pada jumlah garis ketetanggaan ketika k diperbesar?

Semakin besar nilai k, semakin banyak juga garis ketetanggaan yang terbentuk. Pada n = 100 county, jumlah link kurang lebih mengikuti n × k, yaitu sekitar 200 link untuk k = 2, 400 link untuk k = 4, dan 800 link untuk k = 8. Jadi, semakin besar k, jaringan yang terbentuk akan semakin ramai dan padat.

6.0.2 2. Bagaimana pola jaringan k=2 dibandingkan k=8?

Pada k = 2, setiap county hanya memiliki 2 tetangga terdekat sehingga garis yang terbentuk masih sedikit dan jaringan terlihat lebih renggang. Sedangkan pada k = 8, setiap county memiliki lebih banyak tetangga sehingga garis yang muncul jauh lebih banyak dan jaringan terlihat lebih padat. Hubungannya juga bisa sampai ke county yang jaraknya lebih jauh dan tidak selalu berbatasan langsung.

6.0.3 3. Risiko k terlalu kecil?

Jika nilai k terlalu kecil, jumlah tetangga yang digunakan menjadi sangat terbatas. Akibatnya, ada kemungkinan beberapa hubungan antarwilayah yang sebenarnya masih berpengaruh tidak ikut masuk dalam analisis. Hal ini bisa membuat pola hubungan spasial yang terbentuk kurang menggambarkan kondisi sebenarnya dan dapat memengaruhi analisis selanjutnya, seperti Moran’s I.

6.0.4 4. Risiko k terlalu besar?

Sebaliknya, jika k terlalu besar, suatu county akan memiliki terlalu banyak tetangga, termasuk county yang letaknya cukup jauh. Akibatnya, hubungan yang terbentuk tidak lagi terlalu fokus pada kondisi di sekitar county tersebut. Pola spasial yang bersifat lokal juga bisa menjadi kurang terlihat karena terlalu banyak wilayah yang dianggap sebagai tetangga.

7 Latihan 5 — Analisis Data Tahun 1979

nc <- nc %>%
  mutate(SID79_rate = (SID79 / BIR79) * 1000)

summary(nc$SID74_rate)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.084   1.855   2.046   2.604   9.554
summary(nc$SID79_rate)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.249   2.075   2.039   2.539   6.114

7.1 Peta Choropleth Perbandingan 1974 vs 1979

p_1974 <- tm_shape(nc) +
  tm_polygons("SID74_rate",
              palette = "YlOrRd", style = "quantile", n = 5,
              title = "Rate 1974") +
  tm_layout(main.title = "SID74 Rate", main.title.size = 0.9)

p_1979 <- tm_shape(nc) +
  tm_polygons("SID79_rate",
              palette = "YlOrRd", style = "quantile", n = 5,
              title = "Rate 1979") +
  tm_layout(main.title = "SID79 Rate", main.title.size = 0.9)

tmap_arrange(p_1974, p_1979, ncol = 2)

7.2 Struktur Ketetanggaan Queen dan Bobot Spasial (style = “W”)

nb_queen <- poly2nb(nc, queen = TRUE)
listw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
listw_queen
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 100 
## Number of nonzero links: 490 
## Percentage nonzero weights: 4.9 
## Average number of links: 4.9 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 100 10000 100 44.65023 410.4746

7.3 Interpretasi Pola Spasial

cor(nc$SID74_rate, nc$SID79_rate, use = "complete.obs")
## [1] 0.2081463

Jawaban:

  1. Peubah SID79_rate telah dibentuk dengan rumus yang sama seperti SID74_rate, yaitu (SID79 / BIR79) * 1000.
  2. Peta choropleth SID79_rate menggunakan skema warna dan metode kelas yang sama (quantile, n = 5) dengan peta SID74_rate agar keduanya dapat dibandingkan secara adil (fair comparison).
  3. Perbandingan pola spasial 1974 vs 1979: Secara umum, county-county yang menunjukkan rate SIDS tinggi pada tahun 1974 cenderung tidak sepenuhnya sama dengan county yang menunjukkan rate tinggi pada tahun 1979. Beberapa county yang sebelumnya berwarna gelap (rate tinggi) di peta 1974 dapat berubah menjadi lebih terang di peta 1979, dan sebaliknya. Nilai korelasi antara SID74_rate dan SID79_rate (dihitung di atas) dapat digunakan untuk mengukur seberapa konsisten pola tersebut dari waktu ke waktu. Korelasi yang tidak terlalu tinggi mengindikasikan bahwa pola risiko SIDS antarcounty berubah dari tahun ke tahun, bukan pola yang sepenuhnya tetap/stabil secara spasial.
  4. Struktur ketetanggaan dibentuk menggunakan kriteria Queen (poly2nb(nc, queen = TRUE)), lalu bobot spasial dibentuk dengan row-standardization (style = "W") menggunakan nb2listw(). Struktur dan bobot ini nantinya menjadi dasar untuk pengujian autokorelasi spasial (Moran’s I) pada praktikum berikutnya.