1. Data

library(rlang)
library(sp)
library(gstat)
data(meuse)

head(meuse)

##        x      y cadmium copper lead zinc  elev       dist   om ffreq soil lime
## 1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1
## 2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1
## 3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1
## 4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0
## 5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0
## 6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0
##   landuse dist.m
## 1      Ah     50
## 2      Ah     30
## 3      Ah    150
## 4      Ga    270
## 5      Ah    380
## 6      Ga    470

library(spData)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)
library(ggplot2)

data(meuse)        
data(meuse.grid)   

meuse_sf <- st_as_sf(meuse,      coords = c("x", "y"), crs = 28992, remove = FALSE)
grid_sf  <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)
summary(meuse_sf$zinc)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   113.0   198.0   326.0   469.7   674.5  1839.0

ggplot(meuse_sf) +
  geom_sf(aes(color = zinc), size = 2) +
  scale_color_viridis_c(option = "C") +
  labs(title = "Sebaran Konsentrasi Zinc (meuse)", color = "Zinc (ppm)") +
  theme_minimal()

2. Fungsi RMSE

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

3. Fungsi tiap metode interpolasi

pred_closest <- function(train, test) {
  idx <- st_nearest_feature(test, train)
  train$zinc[idx]
}

pred_nn <- function(train, test, k = 5) {
  g <- gstat(formula = zinc ~ 1, data = train, nmax = k, set = list(idp = 0))
  predict(g, newdata = test, debug.level = 0)$var1.pred
}

pred_idw <- function(train, test, idp = 2) {
  g <- gstat(formula = zinc ~ 1, data = train, set = list(idp = idp))
  predict(g, newdata = test, debug.level = 0)$var1.pred
}

fit_variogram <- function(train) {
  v <- variogram(zinc ~ 1, train)

  buat_awal <- function(tipe) {
    vgm(psill  = var(train$zinc) * 0.8,
        model  = tipe,
        range  = max(v$dist) / 3,
        nugget = var(train$zinc) * 0.2)
  }

  hasil <- lapply(c(Sph = "Sph", Exp = "Exp", Gau = "Gau"), function(tipe) {
    fit.variogram(v, model = buat_awal(tipe))
  })

  sse <- sapply(hasil, function(m) attr(m, "SSErr"))

  list(emp = v, model = hasil[[which.min(sse)]], sse = sse)
}

pred_krige <- function(train, test) {
  vm <- fit_variogram(train)$model
  k <- krige(zinc ~ 1, locations = train, newdata = test, model = vm,
             debug.level = 0)
  k$var1.pred
}

4. Variogram

vg <- fit_variogram(meuse_sf)
print(vg$model)

##   model      psill    range
## 1   Nug   9486.627   0.0000
## 2   Exp 163285.480 381.7101

plot(vg$emp, vg$model, main = "Variogram Eksperimental & Model (zinc)")

print(round(vg$sse, 2))

##     Sph     Exp     Gau 
## 2223258 1791466 3729204

Variogram Zinc menunjukkan adanya struktur ketergantungan spasial, yang ditandai dengan meningkatnya semivarian seiring bertambahnya jarak hingga mencapai kondisi relatif stabil. Hal ini mengindikasikan bahwa lokasi yang berdekatan cenderung memiliki konsentrasi Zinc yang lebih mirip dibandingkan lokasi yang berjauhan. Berdasarkan nilai SSErr terkecil, model variogram yang terpilih adalah Exp.

5. K-Fold Cross Validation

set.seed(123)
K    <- 5
n    <- nrow(meuse_sf)
fold <- sample(rep(1:K, length.out = n))
 
cv_pred <- data.frame(
  obs     = meuse_sf$zinc,
  fold    = fold,
  closest = NA_real_,
  nn      = NA_real_,
  idw     = NA_real_,
  kriging = NA_real_
)
 
for (f in 1:K) {
  train <- meuse_sf[fold != f, ]
  test  <- meuse_sf[fold == f, ]
  i     <- which(fold == f)
 
  cv_pred$closest[i] <- pred_closest(train, test)
  cv_pred$nn[i]      <- pred_nn(train, test)
  cv_pred$idw[i]     <- pred_idw(train, test)
  cv_pred$kriging[i] <- pred_krige(train, test)
}

6. Ensemble

cv_pred$ens_mean <- rowMeans(cv_pred[, c("idw", "nn", "kriging")])
rm_comp <- sapply(c("idw", "nn", "kriging"), function(m) rmse(cv_pred$obs, cv_pred[[m]]))
w <- (1 / rm_comp) / sum(1 / rm_comp)
print(round(w, 3))

##     idw      nn kriging 
##   0.311   0.323   0.367

cv_pred$ens_weighted <- as.numeric(as.matrix(cv_pred[, c("idw", "nn", "kriging")]) %*% w)

Ensemble dibangun dari tiga metode, yaitu IDW, Nearest Neighbour, dan Ordinary Kriging. Dua versi dibandingkan:

1. rata-rata sederhana

2. rata-rata terbobot dengan bobot proporsional terhadap 1/RMSE tiap metode

metode dengan kesalahan lebih kecil mendapat bobot lebih besar. Bobot hasil perhitungan ditampilkan pada keluaran di atas.

7. Perbandingan RMSE

null_pred <- sapply(1:n, function(i) mean(meuse_sf$zinc[fold != fold[i]]))  # model nol (rata-rata)
hasil_rmse <- data.frame(
  Metode = c("Model Nol (rata-rata)", "Closest Observation", "Nearest Neighbour (k=5)",
             "IDW (p=2)", "Kriging (Ordinary)",
             "Ensemble (rata-rata)", "Ensemble (terbobot 1/RMSE)"),
  RMSE   = c(rmse(cv_pred$obs, null_pred),
             rmse(cv_pred$obs, cv_pred$closest),
             rmse(cv_pred$obs, cv_pred$nn),
             rmse(cv_pred$obs, cv_pred$idw),
             rmse(cv_pred$obs, cv_pred$kriging),
             rmse(cv_pred$obs, cv_pred$ens_mean),
             rmse(cv_pred$obs, cv_pred$ens_weighted))
) %>%
  mutate(RMSE = round(RMSE, 2),
         Peningkatan_vs_Nol = round((1 - RMSE / RMSE[1]) * 100, 1)) %>%
  arrange(RMSE)

print(hasil_rmse)

##                       Metode   RMSE Peningkatan_vs_Nol
## 1         Kriging (Ordinary) 235.75               35.9
## 2 Ensemble (terbobot 1/RMSE) 250.58               31.9
## 3       Ensemble (rata-rata) 251.93               31.5
## 4    Nearest Neighbour (k=5) 267.83               27.2
## 5                  IDW (p=2) 278.16               24.4
## 6        Closest Observation 351.95                4.4
## 7      Model Nol (rata-rata) 368.01                0.0

ggplot(hasil_rmse, aes(x = reorder(Metode, -RMSE), y = RMSE, fill = RMSE)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = RMSE), hjust = -0.1) +
  coord_flip() +
  scale_fill_viridis_c(direction = -1) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +   
  labs(title = "Perbandingan RMSE (5-Fold Cross Validation)",
       x = NULL, y = "RMSE (zinc, ppm)") +
  theme_minimal()

cat("\nMetode terbaik:", hasil_rmse$Metode[1], "dengan RMSE =", hasil_rmse$RMSE[1], "\n")

## 
## Metode terbaik: Kriging (Ordinary) dengan RMSE = 235.75

r_ <- setNames(hasil_rmse$RMSE, hasil_rmse$Metode)
peningkatan <- setNames(hasil_rmse$Peningkatan_vs_Nol, hasil_rmse$Metode)

Berdasarkan hasil 5-Fold Cross Validation, Kriging (Ordinary) merupakan metode interpolasi terbaik dengan RMSE sebesar 235.75 ppm, atau menurunkan kesalahan sekitar 35.9% dibandingkan model nol (RMSE 368.01 ppm). Pemodelan variogram yang memanfaatkan struktur ketergantungan spasial antar lokasi menghasilkan estimasi Zinc yang lebih akurat dibandingkan metode yang hanya bergantung pada jarak.

Closest Observation hanya memberikan sedikit perbaikan dari model nol (RMSE 351.95 ppm) karena hanya meminjam satu titik terdekat sehingga sangat sensitif terhadap variasi lokal. Nearest Neighbour (267.83 ppm) dan IDW (278.16 ppm) lebih baik karena merata-ratakan beberapa titik. Ensemble (terbobot 250.58 ppm; rata-rata 251.93 ppm) lebih baik dari NN dan IDW, namun belum melampaui Kriging karena ikut memuat metode yang lebih lemah. Selisih antara ensemble terbobot dan rata-rata sederhana sangat kecil, sehingga pembobotan tidak banyak menambah manfaat.

8. Peta hasil interpolasi

peta <- data.frame(st_drop_geometry(grid_sf)[, c("x", "y")])
peta$Closest  <- pred_closest(meuse_sf, grid_sf)
peta$NN       <- pred_nn(meuse_sf, grid_sf)
peta$IDW      <- pred_idw(meuse_sf, grid_sf)
peta$Kriging  <- pred_krige(meuse_sf, grid_sf)
peta$Ensemble <- w["idw"] * peta$IDW + w["nn"] * peta$NN + w["kriging"] * peta$Kriging

metode <- c("Closest", "NN", "IDW", "Ensemble", "Kriging")

peta_long <- do.call(rbind, lapply(metode, function(m) {
  data.frame(x = peta$x, y = peta$y, Metode = m, Zinc = peta[[m]])
}))
peta_long$Metode <- factor(peta_long$Metode, levels = metode)

ggplot(peta_long, aes(x = x, y = y, fill = Zinc)) +
  geom_tile() +
  facet_wrap(~ Metode, ncol = 3) +
  scale_fill_viridis_c(option = "C") +
  coord_equal() +
  labs(title = "Hasil Interpolasi Zinc (data meuse)", fill = "Zinc (ppm)") +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.title = element_blank())

Secara keseluruhan, kelima metode mampu menggambarkan pola persebaran konsentrasi Zinc pada wilayah pengamatan. Namun, terdapat perbedaan karakteristik hasil interpolasi berdasarkan pendekatan masing-masing metode. Jika dikaitkan dengan hasil 5-Fold Cross Validation sebelumnya, Kriging (Ordinary) merupakan metode terbaik dengan RMSE sebesar 235.75 ppm. Dengan demikian, meskipun secara visual beberapa metode menghasilkan pola yang relatif serupa, Ordinary Kriging memberikan tingkat kesalahan prediksi paling rendah sehingga lebih direkomendasikan untuk interpolasi konsentrasi Zinc pada data Meuse.