Pendahuluan

Dataset afcon dari paket spData berisi data Spatial patterns of conflict in Africa 1966-78 untuk 42 negara Afrika. Setiap baris adalah satu negara dengan koordinat titik (x, y). Variabel yang diinterpolasi adalah totcon (total konflik 1966-78).

Tujuan tugas ini adalah membuat permukaan interpolasi totcon di seluruh daratan Afrika dengan metode:

  1. Closest Observation (Voronoi / tetangga terdekat)
  2. Nearest Neighbors (NN)
  3. Inverse Distance Weighting (IDW)
  4. Kriging (Ordinary Kriging)
  5. Ensembel (gabungan Closest Observation, IDW, NN, dan Kriging)

Kinerja tiap metode dibandingkan memakai RMSE dari K-fold cross validation.

Persiapan Data

library(spData)
## Warning: package 'spData' was built under R version 4.6.1
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(sf)
## Warning: package 'sf' was built under R version 4.6.1
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(gstat)
## Warning: package 'gstat' was built under R version 4.6.1
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.6.1
## 
## 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(tidyr)
library(ggplot2)
library(viridis)
## Warning: package 'viridis' was built under R version 4.6.1
## Loading required package: viridisLite
library(patchwork)
data(afcon)
data(world)

# Struktur data
str(afcon)
## 'data.frame':    42 obs. of  5 variables:
##  $ x     : num  9.56 2.63 -6.32 18.02 29.77 ...
##  $ y     : num  34.1 28.2 31.9 27 26.6 ...
##  $ totcon: num  1363 1421 1861 2355 5246 ...
##  $ name  : Factor w/ 42 levels "ALGERIA","ANGOLA",..: 38 1 24 20 11 23 22 26 9 33 ...
##  $ id    : num  2040 2039 2038 2041 2043 ...
head(afcon[, c("name", "totcon", "x", "y")])
##          name totcon          x        y
## TS    TUNISIA   1363   9.556491 34.11685
## AG    ALGERIA   1421   2.631669 28.16258
## MO    MOROCCO   1861  -6.316814 31.88583
## LY      LIBYA   2355  18.022672 27.04399
## EG      EGYPT   5246  29.774473 26.56703
## MR MAURITANIA    811 -10.332211 20.25840
summary(afcon$totcon)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   147.0   620.8   989.0  1350.6  1784.2  5246.0
# Cek koordinat: harus dalam derajat lon/lat agar bisa diproyeksikan
stopifnot(all(c("x", "y", "totcon") %in% names(afcon)))
summary(afcon[, c("x", "y")])
##        x                 y          
##  Min.   :-15.412   Min.   :-29.574  
##  1st Qu.: -0.659   1st Qu.: -3.196  
##  Median : 17.798   Median :  7.100  
##  Mean   : 14.927   Mean   :  4.084  
##  3rd Qu.: 29.910   3rd Qu.: 13.741  
##  Max.   : 45.863   Max.   : 34.117
lonlat <- all(abs(afcon$x) <= 180 & abs(afcon$y) <= 90)
if (!lonlat) stop("Koordinat x,y tampaknya bukan lon/lat. Sesuaikan CRS pada chunk ini.")

# Cek lokasi duplikat (dapat mengganggu variogram)
sum(duplicated(afcon[, c("x", "y")]))
## [1] 0

Transformasi variabel

totcon sangat menjencong ke kanan (banyak negara dengan konflik rendah, beberapa sangat tinggi). Agar interpolasi (terutama Kriging) lebih stabil dan sama seperti contoh log(zinc) pada data meuse, digunakan log(1 + totcon). Seluruh RMSE di bawah ada pada skala log(1 + totcon).

par(mfrow = c(1, 2))
hist(afcon$totcon, main = "totcon", xlab = "totcon", col = "grey80")
hist(log1p(afcon$totcon), main = "log(1 + totcon)", xlab = "log(1 + totcon)", col = "steelblue")

par(mfrow = c(1, 1))

Membuat objek sf dan memproyeksikan ke CRS metrik

Interpolasi berbasis jarak lebih tepat memakai sistem koordinat proyeksi. Digunakan proyeksi Lambert Azimuthal Equal Area yang berpusat di Afrika dengan satuan kilometer.

crs_afr <- "+proj=laea +lat_0=0 +lon_0=20 +datum=WGS84 +units=km +no_defs"

dd <- st_as_sf(afcon, coords = c("x", "y"), crs = 4326) %>%
  st_transform(crs_afr)
dd$vble <- log1p(dd$totcon)

# Batas daratan Afrika
afr <- world %>%
  filter(continent == "Africa") %>%
  st_make_valid() %>%
  st_transform(crs_afr) %>%
  st_make_valid() %>%      # validasi lagi setelah proyeksi
  st_buffer(0)             # bereskan self-intersection tersisa

bnd <- st_sf(geometry = st_union(afr))

ggplot() +
  geom_sf(data = bnd, fill = "grey95") +
  geom_sf(data = dd, aes(color = vble), size = 3) +
  scale_color_viridis(name = "log(1+totcon)") +
  ggtitle("Sebaran data afcon") + theme_bw()

Grid prediksi

# Grid titik berjarak 75 km, hanya yang jatuh di daratan Afrika
grd <- st_make_grid(bnd, cellsize = 75, what = "centers") %>%
  st_sf(geometry = .) %>%
  st_filter(bnd)
nrow(grd)
## [1] 5326

Interpolasi pada Grid

Fungsi bantu untuk memetakan hasil:

plot_pred <- function(pred, judul) {
  ggplot() +
    geom_sf(data = bnd, fill = NA, color = "grey40", linewidth = 0.2) +
    geom_sf(data = pred, aes(color = pred), size = 0.6) +
    geom_sf(data = dd, color = "black", size = 0.8) +
    scale_color_viridis(name = "log(1+totcon)") +
    ggtitle(judul) + theme_bw()
}

1. Closest Observation

Setiap titik grid diberi nilai dari negara (titik data) terdekat. Hasilnya setara dengan poligon Voronoi.

# Ilustrasi diagram Voronoi
vor <- st_voronoi(st_union(dd)) %>%
  st_collection_extract("POLYGON") %>%
  st_sf(geometry = .) %>%
  st_intersection(bnd) %>%
  st_join(dd["vble"])

ggplot() + geom_sf(data = vor, aes(fill = vble), alpha = 0.7) +
  geom_sf(data = dd, size = 1) +
  scale_fill_viridis(name = "log(1+totcon)") +
  ggtitle("Diagram Voronoi") + theme_bw()

idx <- st_nearest_feature(grd, dd)
g1 <- grd; g1$pred <- dd$vble[idx]

2. IDW

Semua titik dipakai dengan pangkat jarak (power) \(p = 2\).

m_idw <- gstat(formula = vble ~ 1, locations = dd, nmax = nrow(dd),
               set = list(idp = 2))
g2 <- grd; g2$pred <- predict(m_idw, grd)$var1.pred
## [inverse distance weighted interpolation]

3. Nearest Neighbors (k = 5)

Rata-rata sederhana dari 5 titik terdekat (idp = 0 berarti bobot sama).

m_nn <- gstat(formula = vble ~ 1, locations = dd, nmax = 5,
              set = list(idp = 0))
g3 <- grd; g3$pred <- predict(m_nn, grd)$var1.pred
## [inverse distance weighted interpolation]

4. Ordinary Kriging

Variogram

vg <- variogram(vble ~ 1, dd)
vg
##    np      dist     gamma dir.hor dir.ver   id
## 1   3  163.4554 0.3161020       0       0 var1
## 2   5  307.4609 0.4325292       0       0 var1
## 3  17  525.0054 0.3722003       0       0 var1
## 4  26  729.8459 0.5571761       0       0 var1
## 5  26  971.8030 0.4181619       0       0 var1
## 6  38 1171.1253 0.3810122       0       0 var1
## 7  39 1380.4835 0.4421323       0       0 var1
## 8  32 1598.7367 0.4579631       0       0 var1
## 9  33 1815.0594 0.4255364       0       0 var1
## 10 39 2041.4281 0.4397625       0       0 var1
## 11 41 2239.8867 0.6858638       0       0 var1
## 12 32 2469.7816 0.6370000       0       0 var1
## 13 38 2687.2430 0.6336390       0       0 var1
## 14 39 2881.1626 0.4323512       0       0 var1
## 15 38 3095.1852 0.6761270       0       0 var1
plot(vg, main = "Variogram empiris log(1+totcon)")

Karena hanya ada 42 titik, variogram empiris cukup berisik. Fungsi berikut mencoba tiga model (Spherical, Exponential, Gaussian) dan memilih yang memiliki galat kuadrat terkecil. Bila semua gagal (singular), digunakan model Exponential dengan nilai awal sebagai cadangan.

fit_vgm <- function(dat) {
  v <- variogram(vble ~ 1, dat)
  s <- var(dat$vble)
  r <- max(v$dist) * 0.5

  best <- NULL; best_sse <- Inf
  for (m in c("Sph", "Exp", "Gau")) {
    f <- tryCatch(
      suppressWarnings(fit.variogram(v, vgm(psill = s * 0.7, model = m,
                                            range = r, nugget = s * 0.1))),
      error = function(e) NULL)
    if (is.null(f) || isTRUE(attr(f, "singular"))) next
    sse <- attr(f, "SSErr")
    if (sse < best_sse) { best <- f; best_sse <- sse }
  }
  if (is.null(best))   # cadangan
    best <- vgm(psill = s * 0.7, model = "Exp", range = r, nugget = s * 0.3)
  best
}

vgf <- fit_vgm(dd)
vgf
##   model     psill    range
## 1   Nug 0.1937366   0.0000
## 2   Exp 0.2868805 273.4003
plot(vg, vgf, cex = 1.5, main = "Variogram empiris dan model terpilih")

Prediksi

m_ok <- gstat(formula = vble ~ 1, locations = dd, model = vgf)
kp   <- predict(m_ok, grd)
## [using ordinary kriging]
g4 <- grd; g4$pred <- kp$var1.pred

# Varians kriging
ggplot() +
  geom_sf(data = bnd, fill = NA, color = "grey40", linewidth = 0.2) +
  geom_sf(data = kp, aes(color = var1.var), size = 0.6) +
  geom_sf(data = dd, size = 0.8) +
  scale_color_viridis(name = "varians", option = "magma") +
  ggtitle("Varians prediksi Ordinary Kriging") + theme_bw()

5. Ensembel

Ensembel menggabungkan empat metode: Closest Observation, IDW, NN, dan Kriging. Pada pemetaan ini digunakan bobot sama. Bobot berbasis kinerja (invers RMSE) dievaluasi pada tahap cross validation.

g5 <- grd; g5$pred <- (g1$pred + g2$pred + g3$pred + g4$pred) / 4

Perbandingan peta

(plot_pred(g1, "Closest Observation") | plot_pred(g2, "IDW (p = 2)")) /
(plot_pred(g3, "Nearest Neighbors (k = 5)") | plot_pred(g4, "Ordinary Kriging")) /
(plot_pred(g5, "Ensembel (bobot sama)") | plot_spacer())

Cross Validation (5-Fold)

Pada setiap fold, semua model dibangun ulang dari data latih (termasuk variogram untuk Kriging), lalu diuji pada data uji. Dua jenis ensembel dihitung:

set.seed(123)

RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))

# Pembagian fold acak
make_folds <- function(n, k) sample(rep(1:k, length.out = n))

# Empat prediksi dasar pada data uji `test` dari data latih `train`
pred_base <- function(train, test) {
  p1 <- train$vble[st_nearest_feature(test, train)]

  p2 <- predict(gstat(formula = vble ~ 1, locations = train,
                      nmax = nrow(train), set = list(idp = 2)),
                test)$var1.pred

  p3 <- predict(gstat(formula = vble ~ 1, locations = train,
                      nmax = 5, set = list(idp = 0)),
                test)$var1.pred

  p4 <- tryCatch(
    predict(gstat(formula = vble ~ 1, locations = train, model = fit_vgm(train)),
            test)$var1.pred,
    error = function(e) p2)   # cadangan jika kriging gagal

  cbind(closest = p1, IDW = p2, NN = p3, kriging = p4)
}

K  <- 5
kf <- make_folds(nrow(dd), K)

hasil <- data.frame(fold = 1:K, closest.obs = NA, IDW = NA, NN = NA,
                    kriging = NA, ens.equal = NA, ens.invRMSE = NA)

for (k in 1:K) {
  train <- dd[kf != k, ]
  test  <- dd[kf == k, ]

  P <- pred_base(train, test)
  hasil[k, c("closest.obs", "IDW", "NN", "kriging")] <-
    apply(P, 2, function(p) RMSE(test$vble, p))

  # Ensembel bobot sama
  hasil$ens.equal[k] <- RMSE(test$vble, rowMeans(P))

  # Ensembel bobot invers RMSE (CV internal pada data latih)
  inner <- make_folds(nrow(train), 5)
  Pin <- matrix(NA, nrow(train), 4)
  for (j in 1:5) {
    Pin[inner == j, ] <- pred_base(train[inner != j, ], train[inner == j, ])
  }
  w <- 1 / apply(Pin, 2, function(p) RMSE(train$vble, p))
  w <- w / sum(w)
  hasil$ens.invRMSE[k] <- RMSE(test$vble, as.vector(P %*% w))
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
# RMSE pada tiap fold
hasil
##   fold closest.obs       IDW        NN   kriging ens.equal ens.invRMSE
## 1    1   1.0016413 0.8924145 0.9356464 0.9409616 0.9174594   0.9169964
## 2    2   0.7501424 0.7113283 0.8380126 0.7478430 0.6991071   0.7034682
## 3    3   0.9862417 0.6154038 0.6486708 0.5559540 0.6555101   0.6412685
## 4    4   0.6568500 0.6156074 0.5842437 0.6592664 0.5858561   0.5860002
## 5    5   0.6483985 0.7661718 0.6769619 0.8830698 0.7199257   0.7240668
# Rata-rata RMSE
rmse_rata <- colMeans(hasil[, -1])
rmse_rata
## closest.obs         IDW          NN     kriging   ens.equal ens.invRMSE 
##   0.8086548   0.7201852   0.7367071   0.7574190   0.7155717   0.7143600

Perbandingan RMSE

# RMSE pada tiap fold

hasil_long <- data.frame(
  fold   = rep(hasil$fold, times = ncol(hasil) - 1),
  metode = rep(names(hasil)[-1], each = nrow(hasil)),
  RMSE   = unlist(hasil[, -1], use.names = FALSE)
)
hasil_long
##    fold      metode      RMSE
## 1     1 closest.obs 1.0016413
## 2     2 closest.obs 0.7501424
## 3     3 closest.obs 0.9862417
## 4     4 closest.obs 0.6568500
## 5     5 closest.obs 0.6483985
## 6     1         IDW 0.8924145
## 7     2         IDW 0.7113283
## 8     3         IDW 0.6154038
## 9     4         IDW 0.6156074
## 10    5         IDW 0.7661718
## 11    1          NN 0.9356464
## 12    2          NN 0.8380126
## 13    3          NN 0.6486708
## 14    4          NN 0.5842437
## 15    5          NN 0.6769619
## 16    1     kriging 0.9409616
## 17    2     kriging 0.7478430
## 18    3     kriging 0.5559540
## 19    4     kriging 0.6592664
## 20    5     kriging 0.8830698
## 21    1   ens.equal 0.9174594
## 22    2   ens.equal 0.6991071
## 23    3   ens.equal 0.6555101
## 24    4   ens.equal 0.5858561
## 25    5   ens.equal 0.7199257
## 26    1 ens.invRMSE 0.9169964
## 27    2 ens.invRMSE 0.7034682
## 28    3 ens.invRMSE 0.6412685
## 29    4 ens.invRMSE 0.5860002
## 30    5 ens.invRMSE 0.7240668
# Menghitung statistik RMSE setiap metode
df_rmse <- hasil_long %>%
  group_by(metode) %>%
  summarise(
    RMSE_mean = mean(RMSE, na.rm = TRUE),
    RMSE_sd = sd(RMSE, na.rm = TRUE),
    RMSE_min = min(RMSE, na.rm = TRUE),
    RMSE_max = max(RMSE, na.rm = TRUE)
  ) %>%
  arrange(RMSE_mean)

df_rmse
## # A tibble: 6 × 5
##   metode      RMSE_mean RMSE_sd RMSE_min RMSE_max
##   <chr>           <dbl>   <dbl>    <dbl>    <dbl>
## 1 ens.invRMSE     0.714   0.126    0.586    0.917
## 2 ens.equal       0.716   0.124    0.586    0.917
## 3 IDW             0.720   0.116    0.615    0.892
## 4 NN              0.737   0.145    0.584    0.936
## 5 kriging         0.757   0.158    0.556    0.941
## 6 closest.obs     0.809   0.174    0.648    1.00
# Visualisasi sebaran RMSE antar-fold
ggplot(hasil_long, aes(x = metode, y = RMSE)) +
  geom_boxplot(fill = "grey85") +
  geom_jitter(width = 0.1,
              color = "steelblue",
              size = 2) +
  labs(
    title = "Sebaran RMSE Antar-Fold",
    x = "Metode Interpolasi",
    y = "RMSE"
  ) +
  theme_bw()

# Metode dengan rata-rata RMSE terkecil
df_rmse %>%
  slice_min(RMSE_mean, n = 1)
## # A tibble: 1 × 5
##   metode      RMSE_mean RMSE_sd RMSE_min RMSE_max
##   <chr>           <dbl>   <dbl>    <dbl>    <dbl>
## 1 ens.invRMSE     0.714   0.126    0.586    0.917

Kesimpulan

Berdasarkan hasil perbandingan nilai RMSE dari enam metode interpolasi, Ensemble Inverse RMSE (ens.invRMSE) menghasilkan nilai RMSE rata-rata terendah sebesar 0,7144, diikuti oleh Ensemble Equal Weight (0,7156), IDW (0,7202), Nearest Neighbor (0,7367), dan Kriging (0,7574). Sementara itu, metode Closest Observation menghasilkan RMSE tertinggi sebesar 0,8087.

Hasil tersebut menunjukkan bahwa metode Ensemble Inverse RMSE memiliki kinerja interpolasi terbaik dalam menghasilkan prediksi dengan kesalahan paling kecil. Meskipun demikian, perbedaan nilai RMSE antara Ensemble Inverse RMSE dan Ensemble Equal Weight relatif kecil. Berdasarkan sebaran RMSE antar-fold, metode Ensemble Inverse RMSE juga menunjukkan variasi kesalahan yang perlu diperhatikan. Dengan demikian, metode Ensemble Inverse RMSE dapat dipilih sebagai metode interpolasi terbaik berdasarkan nilai RMSE rata-rata terkecil, dengan Ensemble Equal Weight sebagai alternatif yang memiliki performa hampir setara.