1 Persiapan

library(sp)        
library(sf)
library(gstat)
library(terra)     
library(dismo)     
library(ggplot2)
library(viridis)
library(knitr)

1.1 Data

data(meuse)
data(meuse.grid)

meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse$logzinc <- log(meuse$zinc)

grid_df <- meuse.grid[, c("x", "y")]
grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)
ggplot() +
  geom_tile(data = grid_df, aes(x, y), fill = "grey90") +
  geom_sf(data = meuse, aes(color = logzinc), size = 2.2) +
  scale_color_viridis_c(name = "log(zinc)") +
  coord_sf(datum = NA) +
  labs(title = "Sebaran titik pengamatan log(zinc) dan grid prediksi") +
  theme_minimal()

2 Metode Deterministik

2.1 Closest Observation (Voronoi)

Diagram Voronoi membagi ruang menjadi poligon, setiap poligon berisi satu titik pengamatan.

v <- voronoi(vect(meuse))

plot(v, main = "Diagram Voronoi data meuse")
points(vect(meuse), cex = 0.5)

v_sf <- st_as_sf(v)

ggplot(v_sf) +
  geom_sf(aes(fill = logzinc), color = "grey40", linewidth = 0.1) +
  scale_fill_viridis_c(name = "log(zinc)", alpha = 0.7) +
  labs(title = "Poligon Voronoi: nilai log(zinc) tiap poligon") +
  theme_minimal()

pred_closest <- function(train, test) {
  vv <- st_as_sf(voronoi(vect(train)))
  g  <- st_geometry(vv)
  st_crs(g) <- NA
  st_crs(g) <- st_crs(train)
  vv <- st_sf(pred = vv$logzinc, geometry = g)

  hit <- st_intersects(test, vv)
  sapply(hit, function(i) if (length(i)) vv$pred[i[1]] else NA_real_)
}

p_closest <- pred_closest(meuse, grid_sf)

2.2 Inverse Distance Weighting (IDW)

Semakin dekat jarak, semakin besar bobot yang diberikan pada nilai titik pengamatan. Menggunakan seluruh titik (nmax = nrow(meuse)) dan pangkat jarak idp.

idp_val <- 1

g_idw <- gstat(formula = logzinc ~ 1, data = meuse,
               nmax = nrow(meuse), set = list(idp = idp_val))

p_idw <- predict(g_idw, grid_sf, debug.level = 0)$var1.pred

2.3 Nearest Neighbors (NN)

Rata-rata nilai dari k = 5 titik terdekat. Semua tetangga berbobot sama, yaitu idp = 0.

k_nn <- 5

g_nn <- gstat(formula = logzinc ~ 1, data = meuse,
              nmax = k_nn, set = list(idp = 0))

p_nn <- predict(g_nn, grid_sf, debug.level = 0)$var1.pred

3 Metode Stokastik: Ordinary Kriging

Kriging memperhitungkan autokorelasi spasial melalui variogram.

3.1 Variogram

vc <- variogram(logzinc ~ 1, meuse, cloud = TRUE)
plot(vc, main = "Variogram cloud")

v_emp <- variogram(logzinc ~ 1, data = meuse)
plot(v_emp, main = "Variogram empiris")

Model variogram awal ditentukan dari bentuk variogram empiris: sill sekitar 0,5, model Spherical, range sekitar 900, dan nugget sekitar 0,1. Kemudian model di-fit.

vinitial <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)

fv <- fit.variogram(object = v_emp, model = vinitial)
fv
plot(v_emp, fv, cex = 1.5, main = "Variogram empiris dan model hasil fitting")

3.2 Prediksi pada Grid

g_krige <- gstat(formula = logzinc ~ 1, data = meuse, model = fv)
kpred   <- predict(g_krige, grid_sf, debug.level = 0)

p_krige <- kpred$var1.pred
ggplot() +
  geom_tile(data = cbind(grid_df, var = kpred$var1.var), aes(x, y, fill = var)) +
  geom_sf(data = meuse, size = 0.8) +
  scale_fill_viridis_c(name = "variance") +
  coord_sf(datum = NA) +
  labs(title = "Varians prediksi Kriging") +
  theme_minimal()

4 Ensemble

Ensemble menggabungkan beberapa metode. Di sini digunakan empat metode (Closest Observation, IDW, NN, dan Kriging) dengan bobot sama (1/4) untuk peta pada grid. Pada Cross Validation nanti, dibandingkan juga ensemble dengan bobot 1/RMSE.

P_grid <- cbind(Closest = p_closest, IDW = p_idw, NN = p_nn, Kriging = p_krige)
p_ens  <- rowMeans(P_grid)

5 Peta Hasil Interpolasi

hasil <- list(
  "Closest Observation" = p_closest,
  "IDW"                 = p_idw,
  "Nearest Neighbors"   = p_nn,
  "Ensemble"            = p_ens,
  "Kriging"             = p_krige
)

peta <- do.call(rbind, lapply(names(hasil), function(m) {
  data.frame(x = grid_df$x, y = grid_df$y, metode = m, pred = hasil[[m]])
}))
peta$metode <- factor(peta$metode, levels = names(hasil))

ggplot(peta, aes(x, y, fill = pred)) +
  geom_tile() +
  coord_equal() +
  facet_wrap(~ metode, ncol = 2) +
  scale_fill_viridis_c(name = "log(zinc)") +
  labs(title = "Perbandingan hasil interpolasi log(zinc)", x = NULL, y = NULL) +
  theme_minimal()

6 Evaluasi Model: K-fold Cross Validation

Data dibagi menjadi K = 5 lipatan dengan dismo::kfold(). Pada setiap iterasi satu lipatan menjadi data uji dan empat lipatan lainnya menjadi data latih. Variogram untuk Kriging di-fit ulang pada data latih di setiap lipatan.

Dua versi Ensemble dibandingkan:

  • Ensemble (bobot sama): rata-rata sederhana keempat metode.
  • Ensemble (bobot 1/RMSE): bobot berbanding terbalik dengan RMSE tiap metode pada lipatan tersebut (semakin kecil RMSE, semakin besar bobot), seperti pada materi praktikum.
set.seed(123)

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

K  <- 5
kf <- kfold(nrow(meuse), k = K)

metode <- c("Closest Observation", "IDW", "Nearest Neighbors", "Kriging",
            "Ensemble (bobot sama)", "Ensemble (bobot 1/RMSE)")
rmse_cv <- matrix(NA_real_, nrow = K, ncol = length(metode),
                  dimnames = list(paste("Fold", 1:K), metode))

for (k in 1:K) {
  test  <- meuse[kf == k, ]
  train <- meuse[kf != k, ]
  obs   <- test$logzinc

  # Closest observation
  p1 <- pred_closest(train, test)

  # IDW
  gi <- gstat(formula = logzinc ~ 1, data = train,
              nmax = nrow(train), set = list(idp = idp_val))
  p2 <- predict(gi, test, debug.level = 0)$var1.pred

  # Nearest neighbors
  gn <- gstat(formula = logzinc ~ 1, data = train,
              nmax = k_nn, set = list(idp = 0))
  p3 <- predict(gn, test, debug.level = 0)$var1.pred

  # Kriging (variogram di-fit ulang pada data latih)
  fv_k <- fit.variogram(variogram(logzinc ~ 1, train), vinitial)
  gk   <- gstat(formula = logzinc ~ 1, data = train, model = fv_k)
  p4   <- predict(gk, test, debug.level = 0)$var1.pred

  # RMSE tiap metode
  P <- cbind(p1, p2, p3, p4)
  r <- apply(P, 2, function(p) RMSE(obs, p))

  # Ensemble
  p_eq <- rowMeans(P)
  w    <- (1 / r) / sum(1 / r)
  p_w  <- as.vector(P %*% w)

  rmse_cv[k, ] <- c(r, RMSE(obs, p_eq), RMSE(obs, p_w))
}

6.1 RMSE pada Setiap Lipatan

kable(round(rmse_cv, 4), caption = "RMSE tiap lipatan (skala log zinc)")
RMSE tiap lipatan (skala log zinc)
Closest Observation IDW Nearest Neighbors Kriging Ensemble (bobot sama) Ensemble (bobot 1/RMSE)
Fold 1 0.6327 0.5496 0.4958 0.4281 0.4508 0.4458
Fold 2 0.5987 0.5828 0.4943 0.4183 0.4492 0.4439
Fold 3 0.6328 0.7843 0.6516 0.5055 0.5995 0.5864
Fold 4 0.4874 0.6333 0.5316 0.4198 0.4633 0.4542
Fold 5 0.5593 0.6293 0.4791 0.3840 0.4463 0.4347

6.2 Rata-rata RMSE

tabel_rmse <- data.frame(
  Metode = colnames(rmse_cv),
  RMSE   = colMeans(rmse_cv)
)
tabel_rmse <- tabel_rmse[order(tabel_rmse$RMSE), ]
tabel_rmse$Peringkat <- seq_len(nrow(tabel_rmse))

kable(tabel_rmse, digits = 4, row.names = FALSE,
      caption = "Rata-rata RMSE dari 5-fold Cross Validation")
Rata-rata RMSE dari 5-fold Cross Validation
Metode RMSE Peringkat
Kriging 0.4312 1
Ensemble (bobot 1/RMSE) 0.4730 2
Ensemble (bobot sama) 0.4818 3
Nearest Neighbors 0.5305 4
Closest Observation 0.5822 5
IDW 0.6358 6 
ggplot(tabel_rmse, aes(x = reorder(Metode, RMSE), y = RMSE, fill = Metode)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 3)), vjust = -0.4) +
  scale_fill_viridis_d() +
  labs(title = "Rata-rata RMSE per metode (5-fold CV)", x = NULL, y = "RMSE") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

7 Interpretasi

Pola datanya: Di peta, zinc yang tinggi mengumpul di bagian tepi wilayah (barat, utara, dan selatan), sedangkan bagian tengah sampai timur nilainya lebih rendah. Karena nilainya mengelompok dan tidak acak, lokasi yang berdekatan punya nilai yang mirip, jadi interpolasi cocok digunakan. Didukung juga dengan Variogramnya, dengan range sekitar 897 m dan nugget-nya kecil (0,0507 dibanding sill 0,6413), artinya ketergantungan spasialnya kuat.

Hasil perbandingan. Urutan dari RMSE terkecil: Kriging (0,4312), Ensemble bobot 1/RMSE (0,4730), Ensemble bobot sama (0,4818), NN (0,5305), Closest Observation (0,5822), dan IDW (0,6358).

  • Kriging paling bagus, dan terbaik di kelima fold. Hal ini karena Kriging memakai pola kemiripan nilai terhadap jarak dari variogram, bukan hanya jarak.
  • Ensemble lebih baik daripada metode deterministik, tetapi belum bisa mengalahkan Kriging karena ikut tertarik oleh metode yang lemah (IDW dan Closest Observation).
  • NN yang terbaik di antara metode deterministik, karena merata-ratakan 5 tetangga sehingga lebih stabil.
  • Closest Observation hanya melihat satu titik terdekat, jadi petanya kotak-kotak dan mudah terpengaruh satu titik yang nilainya aneh.
  • IDW paling jelek. Petanya hampir rata, karena dengan idp = 1 dan semua titik dipakai, titik yang jauh masih berpengaruh besar sehingga prediksi tertarik ke rata-rata.

8 Kesimpulan

Untuk data meuse, Kriging paling cocok dengan RMSE 0,4312. Untuk IDW, sebaiknya dicoba dengan parameter lain (misalnya idp lebih besar) agar hasilnya lebih baik.