1 Data

library(sp)
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(terra)
## terra 1.8.93
library(gstat)
library(automap)
library(tmap)

data(meuse)
data(meuse.grid)

meuse_sf      <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse_grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)
meuse_sf$vble <- log(meuse_sf$zinc)

map_bound <- meuse_grid_sf |>
  st_buffer(20, endCapStyle = "SQUARE") |>
  st_union() |>
  st_sf()

grid <- rast(meuse.grid[, c("x", "y", "dist")], type = "xyz", crs = "EPSG:28992")

RMSE <- function(observed, predicted) sqrt(mean((observed - predicted)^2))
plot_pred <- function(pred_values, judul) {
  pts <- vect(meuse_grid_sf)
  pts$pred <- pred_values
  r <- rasterize(pts, grid, field = "pred")
  r <- mask(r, vect(map_bound))
  tm_shape(r) +
    tm_raster(col.scale = tm_scale_continuous(values = "viridis"),
              col.legend = tm_legend(title = "log(Zn)"),
              col_alpha = 0.6) +
    tm_title(judul)
}
pred_voronoi <- function(train, newdata) {
  train$vble[st_nearest_feature(newdata, train)]
}

pred_idw <- function(train, newdata, idp = 1) {
  m <- gstat(formula = vble ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = idp))
  predict(m, newdata)$var1.pred
}

pred_nn <- function(train, newdata, nmax = 5) {
  m <- gstat(formula = vble ~ 1, locations = train,
             nmax = nmax, set = list(idp = 0))
  predict(m, newdata)$var1.pred
}

pred_ok <- function(train, newdata) {
  afvg <- autofitVariogram(vble ~ 1, as(train, "Spatial"))
  m <- gstat(formula = vble ~ 1, locations = train, model = afvg$var_model)
  predict(m, newdata)$var1.pred
}

set.seed(123)
kf <- dismo::kfold(nrow(meuse_sf), k = 5)

rmse1 <- rmse2 <- rmse3 <- rmse4 <- rmse5 <- rep(NA_real_, 5)

2 Closest Observation (Voronoi)

v <- terra::voronoi(vect(meuse_sf), bnd = vect(map_bound))
pred_c1 <- rasterize(v, grid, field = "vble")
pred_c1 <- mask(pred_c1, vect(map_bound))

c1 <- tm_shape(pred_c1) +
  tm_raster(col.scale = tm_scale_continuous(values = "viridis"),
            col.legend = tm_legend(title = "log(Zn)"),
            col_alpha = 0.6) +
  tm_title("Closest Observation (Voronoi)")
print(c1)

for (k in 1:5) {
  test  <- meuse_sf[kf == k, ]
  train <- meuse_sf[kf != k, ]
  rmse1[k] <- RMSE(test$vble, pred_voronoi(train, test))
}

3 Inverse Distance Weighting (IDW)

c2 <- plot_pred(pred_idw(meuse_sf, meuse_grid_sf), "IDW")
## [inverse distance weighted interpolation]
print(c2)

for (k in 1:5) {
  test  <- meuse_sf[kf == k, ]
  train <- meuse_sf[kf != k, ]
  rmse2[k] <- RMSE(test$vble, pred_idw(train, test))
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

4 Nearest Neighbor (NN)

c3 <- plot_pred(pred_nn(meuse_sf, meuse_grid_sf), "Nearest Neighbors (nmax = 5)")
## [inverse distance weighted interpolation]
print(c3)

# Cross validation
for (k in 1:5) {
  test  <- meuse_sf[kf == k, ]
  train <- meuse_sf[kf != k, ]
  rmse3[k] <- RMSE(test$vble, pred_nn(train, test))
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

5 Ensembel

Bobot tiap metode = kebalikan RMSE-nya (metode yang lebih akurat bobotnya lebih besar).

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

  p1 <- pred_voronoi(train, test)
  p2 <- pred_idw(train, test)
  p3 <- pred_nn(train, test)

  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k])
  weights <- w / sum(w)
  rmse4[k] <- RMSE(test$vble, p1 * weights[1] + p2 * weights[2] + p3 * weights[3])
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
# Peta (bobot dari rata-rata RMSE seluruh fold)
weights_all <- 1 / c(mean(rmse1), mean(rmse2), mean(rmse3))
weights_all <- weights_all / sum(weights_all)

p4_all <- pred_voronoi(meuse_sf, meuse_grid_sf) * weights_all[1] +
          pred_idw(meuse_sf, meuse_grid_sf)     * weights_all[2] +
          pred_nn(meuse_sf, meuse_grid_sf)      * weights_all[3]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
c4 <- plot_pred(p4_all, "Ensemble Interpolation")
print(c4)

6 Kriging

c5 <- plot_pred(pred_ok(meuse_sf, meuse_grid_sf), "Ordinary Kriging")
## [using ordinary kriging]
print(c5)

# Cross validation
for (k in 1:5) {
  test  <- meuse_sf[kf == k, ]
  train <- meuse_sf[kf != k, ]
  rmse5[k] <- tryCatch(RMSE(test$vble, pred_ok(train, test)),
                       error = function(e) NA_real_)
}
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]

7 Evaluasi

df_rmse_folds <- data.frame(
  Fold             = 1:5,
  Closest_Obs      = rmse1,
  IDW              = rmse2,
  Nearest_Neighbor = rmse3,
  Ensemble         = rmse4,
  Ordinary_Kriging = rmse5
)
print(df_rmse_folds)
##   Fold Closest_Obs       IDW Nearest_Neighbor  Ensemble Ordinary_Kriging
## 1    1   0.6327142 0.5495639        0.4957726 0.4686157        0.4263603
## 2    2   0.5986922 0.5828124        0.4943376 0.4723239        0.4136546
## 3    3   0.6328388 0.7842703        0.6516044 0.6317812        0.5116595
## 4    4   0.4874366 0.6332634        0.5315858 0.4851121        0.4224304
## 5    5   0.5593159 0.6292808        0.4791353 0.4771908        0.3777636
df_mean_rmse <- data.frame(
  Closest_Obs      = mean(rmse1, na.rm = TRUE),
  IDW              = mean(rmse2, na.rm = TRUE),
  Nearest_Neighbor = mean(rmse3, na.rm = TRUE),
  Ensemble         = mean(rmse4, na.rm = TRUE),
  Ordinary_Kriging = mean(rmse5, na.rm = TRUE)
)
print(df_mean_rmse)
##   Closest_Obs       IDW Nearest_Neighbor  Ensemble Ordinary_Kriging
## 1   0.5821996 0.6358382        0.5304872 0.5070047        0.4303737

8 Kesimpulan

Berdasarkan hasil validasi silang 5 fold, Ordinary Kriging menjadi metode interpolasi paling akurat untuk memprediksi pada data meuse, dengan rata-rata RMSE terkecil (0,430) dan nilai RMSE terendah di setiap fold (0,378-0,512). Diikuti oleh Ensemble (0,507), Nearest Neighbor (0,530), Closest Observation/Voronoi (0,582), dan IDW (0,636) sebagai yang paling kurang akurat. Keunggulan Kriging kemungkinan besar karena metode ini memodelkan struktur autokorelasi spasial lewat variogram, sedangkan metode lain hanya mengandalkan jarak atau kedekatan titik. Pendekatan Ensemble berhasil memperbaiki kinerja metode penyusunnya (Voronoi, IDW, dan NN) dengan menggabungkan prediksinya, tetapi belum cukup untuk mengalahkan Kriging. Dengan demikian, Ordinary Kriging direkomendasikan untuk interpolasi data ini, dan fold 3 yang RMSE-nya relatif tinggi di semua metode menunjukkan bahwa kualitas prediksi juga dipengaruhi oleh sebaran titik uji pada tiap fold.