1 Pendahuluan

Tugas ini menerapkan interpolasi spasial pada data meuse (kadar seng di tanah sekitar sungai Meuse, Belanda). Peubah yang diinterpolasi adalah log(zinc). Metode yang dibandingkan:

Closest Observation (Voronoi)
Inverse Distance Weighting (IDW)
Nearest Neighbors (NN)
Ensembel (gabungan Closest Observation, IDW, dan NN)
Ordinary Kriging

Kinerja setiap metode dibandingkan dengan RMSE dari 5-fold cross validation.

2 Persiapan

# install.packages(c("sp", "sf", "gstat", "ggplot2", "viridis", "knitr"))
library(sp)
library(sf)
library(gstat)
library(ggplot2)
library(viridis)
library(knitr)

2.1 Data

data(meuse)
data(meuse.grid)

# Titik pengamatan dan lokasi prediksi (grid)
d    <- st_as_sf(meuse,      coords = c("x", "y"), crs = 28992)
coop <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Peubah yang diinterpolasi
d$vble <- log(d$zinc)

summary(d$vble)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.727   5.288   5.787   5.886   6.514   7.517

ggplot() +
  geom_sf(data = coop, color = "grey85", size = 0.3) +
  geom_sf(data = d, aes(color = vble), size = 2) +
  scale_color_viridis(name = "log(zinc)") +
  theme_bw() +
  labs(title = "Sebaran titik pengamatan log(zinc)")

3 Fungsi Interpolasi

Setiap metode dibuat sebagai fungsi dengan argumen train (data latih) dan newloc (lokasi yang diprediksi), sehingga bisa dipakai ulang untuk peta maupun cross validation.

RMSE <- function(observed, predicted) sqrt(mean((observed - predicted)^2))

# 1. Closest observation: nilai titik terdekat.
#    Hasilnya identik dengan diagram Voronoi.
pred_closest <- function(train, newloc) {
  idx <- st_nearest_feature(newloc, train)
  train$vble[idx]
}

# 2. IDW (idp = 1, memakai semua titik)
pred_idw <- function(train, newloc) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = 1))
  predict(g, newloc, debug.level = 0)$var1.pred
}

# 3. Nearest neighbors (rata-rata k tetangga terdekat, bobot sama)
pred_nn <- function(train, newloc, k = 5) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = k, set = list(idp = 0))
  predict(g, newloc, debug.level = 0)$var1.pred
}

# 5. Ordinary Kriging (variogram dipasang pada data latih)
fit_vgm <- function(train) {
  v <- variogram(vble ~ 1, data = train)
  fit.variogram(v, model = vgm(psill = 0.5, model = "Sph",
                               range = 900, nugget = 0.1))
}
pred_kriging <- function(train, newloc, fv = fit_vgm(train)) {
  k <- gstat(formula = vble ~ 1, data = train, model = fv)
  predict(k, newloc, debug.level = 0)$var1.pred
}

4 Variogram dan Kriging

v      <- variogram(vble ~ 1, data = d)
fv_all <- fit_vgm(d)
fv_all

plot(v, fv_all, cex = 1.5)

Variogram eksperimental dipasang dengan model Spherical. Nilai nugget, sill, dan range pada tabel di atas menunjukkan struktur ketergantungan spasial: semivarian naik seiring jarak lalu mencapai sill pada jarak sekitar nilai range.

5 Peta Prediksi

Ensembel pada bagian ini memakai bobot sama (1/3 untuk Closest Observation, IDW, dan NN).

p1 <- pred_closest(d, coop)
p2 <- pred_idw(d, coop)
p3 <- pred_nn(d, coop)
p4 <- (p1 + p2 + p3) / 3
p5 <- pred_kriging(d, coop, fv_all)

nama_metode <- c("Closest Observation", "IDW", "Nearest Neighbors",
                 "Ensembel (bobot sama)", "Ordinary Kriging")
xy <- st_coordinates(coop)

hasil <- data.frame(
  x      = rep(xy[, 1], 5),
  y      = rep(xy[, 2], 5),
  pred   = c(p1, p2, p3, p4, p5),
  metode = factor(rep(nama_metode, each = nrow(xy)), levels = nama_metode)
)

ggplot(hasil, aes(x, y, fill = pred)) +
  geom_raster() +
  geom_point(data = data.frame(st_coordinates(d)), aes(X, Y),
             inherit.aes = FALSE, size = 0.4) +
  scale_fill_viridis(name = "log(zinc)") +
  coord_equal() +
  facet_wrap(~ metode, ncol = 3) +
  theme_bw() +
  labs(x = NULL, y = NULL,
       title = "Prediksi log(zinc) menurut metode interpolasi")

# Varians kriging (ketidakpastian prediksi)
k_all <- gstat(formula = vble ~ 1, data = d, model = fv_all)
kpred <- predict(k_all, coop, debug.level = 0)

ggplot() +
  geom_sf(data = kpred, aes(color = var1.var), size = 0.8) +
  geom_sf(data = d, size = 0.8) +
  scale_color_viridis(name = "Varians") +
  theme_bw() +
  labs(title = "Varians prediksi Ordinary Kriging")

6 Evaluasi dengan Cross Validation

Data dibagi menjadi 5 fold. Pada tiap putaran, satu fold menjadi data uji dan empat fold lainnya menjadi data latih. Dua versi ensembel dihitung:

Ensembel bobot sama: rata-rata sederhana ketiga metode.
Ensembel bobot kebalikan RMSE: metode dengan RMSE lebih kecil mendapat bobot lebih besar (seperti contoh di materi).

set.seed(123)
K  <- 5
kf <- sample(rep(1:K, length.out = nrow(d)))

rmse <- data.frame(fold = 1:K,
                   closest.obs = NA_real_, IDW = NA_real_,
                   nearest.neigh = NA_real_, ensemble_equal = NA_real_,
                   ensemble_invRMSE = NA_real_, kriging = NA_real_)

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

  q1 <- pred_closest(train, test)
  q2 <- pred_idw(train, test)
  q3 <- pred_nn(train, test)
  q5 <- pred_kriging(train, test)

  rmse$closest.obs[k]    <- RMSE(test$vble, q1)
  rmse$IDW[k]            <- RMSE(test$vble, q2)
  rmse$nearest.neigh[k]  <- RMSE(test$vble, q3)
  rmse$kriging[k]        <- RMSE(test$vble, q5)
  rmse$ensemble_equal[k] <- RMSE(test$vble, (q1 + q2 + q3) / 3)

  w <- 1 / c(rmse$closest.obs[k], rmse$IDW[k], rmse$nearest.neigh[k])
  w <- w / sum(w)
  rmse$ensemble_invRMSE[k] <- RMSE(test$vble, w[1] * q1 + w[2] * q2 + w[3] * q3)
}

6.1 RMSE tiap fold

kable(rmse, digits = 4, caption = "RMSE tiap fold (skala log(zinc))")

RMSE tiap fold (skala log(zinc))
fold	closest.obs	IDW	nearest.neigh	ensemble_equal	ensemble_invRMSE	kriging
1	0.5710	0.6870	0.5261	0.5251	0.5188	0.4213
2	0.4728	0.5705	0.3685	0.3981	0.3878	0.2919
3	0.7532	0.6398	0.5292	0.5656	0.5555	0.4411
4	0.7641	0.6585	0.5305	0.5671	0.5578	0.4685
5	0.6685	0.6404	0.5223	0.5283	0.5244	0.4157

6.2 Rata-rata RMSE

rata <- colMeans(rmse[, -1])
label <- c(closest.obs       = "Closest Observation",
           IDW               = "IDW",
           nearest.neigh     = "Nearest Neighbors",
           ensemble_equal    = "Ensembel (bobot sama)",
           ensemble_invRMSE  = "Ensembel (bobot kebalikan RMSE)",
           kriging           = "Ordinary Kriging")
rata_df <- data.frame(Metode = unname(label[names(rata)]),
                      RMSE = round(unname(rata), 4))
rata_df <- rata_df[order(rata_df$RMSE), ]
kable(rata_df, row.names = FALSE, caption = "Rata-rata RMSE 5 fold (urut terkecil)")

Rata-rata RMSE 5 fold (urut terkecil)
Metode	RMSE
Ordinary Kriging	0.4077
Nearest Neighbors	0.4953
Ensembel (bobot kebalikan RMSE)	0.5089
Ensembel (bobot sama)	0.5169
IDW	0.6393
Closest Observation	0.6459

ggplot(rata_df, aes(x = reorder(Metode, -RMSE), y = RMSE)) +
  geom_col(fill = "#3B6E8F") +
  geom_text(aes(label = RMSE), hjust = -0.1, size = 3.5) +
  coord_flip(ylim = c(0, max(rata_df$RMSE) * 1.15)) +
  theme_bw() +
  labs(x = NULL, y = "Rata-rata RMSE",
       title = "Perbandingan RMSE antar metode")

7 Kesimpulan

Berdasarkan rata-rata RMSE dari 5-fold cross validation, metode dengan kinerja terbaik adalah Ordinary Kriging (RMSE = 0.4077), sedangkan yang terendah adalah Closest Observation (RMSE = 0.6459).

7.1 Interpretasi

Kriging paling akurat. Ordinary Kriging memiliki RMSE terkecil dan unggul di kelima fold. Ini sesuai dengan sifat datanya: kadar seng memiliki struktur ketergantungan spasial yang jelas. Model Spherical yang dipasang memiliki nugget 0.051, partial sill 0.591 (sill total sekitar 0.641), dan range sekitar 897 meter. Nugget hanya sekitar 8% dari sill total, artinya ketergantungan spasialnya kuat: titik yang berdekatan memang bernilai mirip, dan pengaruh itu hilang pada jarak kurang lebih 900 meter. Kriging adalah satu-satunya metode yang memakai struktur itu secara eksplisit lewat model variogram. Metode deterministik hanya mengandalkan jarak, tanpa mempelajari pola ketergantungan dari data.
Di antara metode deterministik, Nearest Neighbors terbaik. Rata-rata 5 tetangga terdekat meredam pengaruh satu titik yang kebetulan ekstrem. Closest Observation menyalin nilai satu titik saja, sehingga paling sensitif terhadap noise. IDW (idp = 1) memakai seluruh titik, termasuk titik yang jauh, sehingga hasilnya cenderung terlalu halus.
Ensembel tidak mengalahkan metode terbaik penyusunnya. RMSE ensembel (baik bobot sama maupun bobot kebalikan RMSE) berada di antara metode penyusunnya dan masih lebih besar daripada Nearest Neighbors. Penyebabnya, ensembel mencampur dua metode yang lebih lemah (Closest Observation dan IDW) dengan satu metode yang lebih kuat, sehingga ketepatannya tertarik turun. Bobot kebalikan RMSE sedikit memperbaiki hasil karena memberi bobot lebih besar pada metode yang lebih baik, tetapi tetap tidak menyaingi Kriging. Ensembel lebih bermanfaat bila metode penyusunnya setara kinerjanya dan kesalahannya berbeda pola.
Keunggulan tambahan Kriging. Kriging juga menghasilkan varians prediksi, sehingga kita tahu di area mana prediksi kurang pasti (biasanya di tepi wilayah dan area yang jarang titik pengamatannya). Metode deterministik tidak memberi ukuran ketidakpastian ini.

Kesimpulan akhir: untuk data meuse, Ordinary Kriging dengan model Spherical adalah metode interpolasi yang paling direkomendasikan.

Tugas Pertemuan 9: Interpolasi Spasial pada Data Meuse

Darell Belva / M0401241071

10 November 2026