Library

library(sp)

## Warning: package 'sp' was built under R version 4.5.3

library(sf)

## Warning: package 'sf' was built under R version 4.5.3

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(terra)

## Warning: package 'terra' was built under R version 4.5.3

## terra 1.9.34

library(gstat)

## Warning: package 'gstat' was built under R version 4.5.3

library(dismo)

## Warning: package 'dismo' was built under R version 4.5.3

## Loading required package: raster

## Warning: package 'raster' was built under R version 4.5.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.5.2

library(viridis)

## Warning: package 'viridis' was built under R version 4.5.3

## Loading required package: viridisLite

## Warning: package 'viridisLite' was built under R version 4.5.2

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.5.2

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:raster':
## 
##     area

## The following object is masked from 'package:terra':
## 
##     area

Data

data(meuse)
data(meuse.grid)
 
meuse$logzinc <- log(meuse$zinc)
 
d    <- st_as_sf(meuse,      coords = c("x", "y"), crs = 28992)
coop <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

semua_titik <- as.data.frame(rbind(st_coordinates(d), st_coordinates(coop)))
map <- st_as_sf(semua_titik, coords = c("X", "Y"), crs = 28992) %>%
  st_combine() %>%
  st_convex_hull() %>%
  st_sf()
 
ggplot() +
  geom_sf(data = coop, color = "grey85", size = 0.4) +
  geom_sf(data = d, aes(color = logzinc), size = 1.8) +
  scale_color_viridis(name = "log(zinc)") +
  theme_bw()

Fungsi prediksi per metode

rmse <- function(observed, predicted) sqrt(mean((observed - predicted)^2))
 
K_NEIGHBOR <- 5   # jumlah tetangga untuk NN
IDP_IDW    <- 1   # pangkat jarak IDW (materi memakai 1)
 

pred_closest <- function(train, newdata) {
  train$logzinc[st_nearest_feature(newdata, train)]
}
 
pred_idw <- function(train, newdata) {
  g <- gstat(formula = logzinc ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = IDP_IDW))
  predict(g, newdata, debug.level = 0)$var1.pred
}
 
pred_nn <- function(train, newdata) {
  g <- gstat(formula = logzinc ~ 1, locations = train,
             nmax = K_NEIGHBOR, set = list(idp = 0))
  predict(g, newdata, debug.level = 0)$var1.pred
}
 
fit_variogram_sph <- function(train) {
  v <- variogram(logzinc ~ 1, data = train)
  fit.variogram(v, model = vgm(psill = 0.5, model = "Sph",
                               range = 900, nugget = 0.1))
}
 
pred_kriging <- function(train, newdata) {
  fv <- fit_variogram_sph(train)
  krige(logzinc ~ 1, locations = train, newdata = newdata,
        model = fv, debug.level = 0)$var1.pred
}

Closest Observation (Voronoi)

vor <- voronoi(x = vect(d), bnd = vect(map)) %>% st_as_sf()
 
ggplot(vor) +
  geom_sf(aes(fill = logzinc), color = "white", linewidth = 0.1) +
  geom_sf(data = d, size = 0.5) +
  scale_fill_viridis(name = "log(zinc)") +
  theme_bw()

coop$closest <- pred_closest(d, coop)

IDW & Nearest Neighbors

coop$idw <- pred_idw(d, coop)
coop$nn  <- pred_nn(d, coop)

Ensemble

inverse_rmse_weights <- function(train, k = 5) {
  inner <- kfold(nrow(train), k = k)
  p <- matrix(NA_real_, nrow(train), 3)
  for (j in 1:k) {
    tr <- train[inner != j, ]
    te <- train[inner == j, ]
    p[inner == j, ] <- cbind(pred_closest(tr, te),
                             pred_idw(tr, te),
                             pred_nn(tr, te))
  }
  w <- 1 / apply(p, 2, function(m) rmse(train$logzinc, m))
  w / sum(w)
}
 
pred_ensemble <- function(train, newdata, w = inverse_rmse_weights(train)) {
  w[1] * pred_closest(train, newdata) +
    w[2] * pred_idw(train, newdata) +
    w[3] * pred_nn(train, newdata)
}
 
w_full <- inverse_rmse_weights(d)
names(w_full) <- c("closest", "idw", "nn")
w_full   # bobot ensemble pada seluruh data

##   closest       idw        nn 
## 0.3307573 0.2819060 0.3873367

coop$ensemble <- w_full["closest"] * coop$closest +
  w_full["idw"] * coop$idw +
  w_full["nn"] * coop$nn

Ordinary Kriging

v_emp <- variogram(logzinc ~ 1, data = d)
fv    <- fit_variogram_sph(d)
fv

##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044

plot(v_emp, fv, cex = 1.5)

krig <- gstat(formula = logzinc ~ 1, data = d, model = fv)
kp   <- predict(krig, coop)

## [using ordinary kriging]

coop$kriging <- kp$var1.pred
coop$krig_var <- kp$var1.var

Peta perbandingan semua metode

peta <- function(kolom, judul) {
  ggplot(coop) +
    geom_sf(aes(color = .data[[kolom]]), size = 0.7, shape = 15) +
    scale_color_viridis(name = "log(zinc)", limits = range(d$logzinc)) +
    labs(title = judul) +
    theme_bw()
}
 
(peta("closest", "Closest Observation") | peta("nn", "Nearest Neighbors")) /
(peta("idw", "IDW") | peta("ensemble", "Ensemble")) /
(peta("kriging", "Ordinary Kriging") | plot_spacer()) +
  plot_layout(guides = "collect")

ggplot(coop) +
  geom_sf(aes(color = krig_var), size = 0.7, shape = 15) +
  geom_sf(data = d, size = 0.6) +
  scale_color_viridis(name = "Varians kriging", option = "magma") +
  theme_bw()

5-fold Cross Validation

set.seed(123)
K  <- 5
kf <- kfold(nrow(d), k = K)    # fold yang sama dipakai semua metode
 
metode <- c("Closest", "IDW", "NN", "Ensemble", "Kriging")
rmse_cv <- matrix(NA_real_, nrow = K, ncol = length(metode),
                  dimnames = list(paste0("fold", 1:K), metode))
 
for (k in 1:K) {
  train <- d[kf != k, ]
  test  <- d[kf == k, ]
 
  p_closest <- pred_closest(train, test)
  p_idw     <- pred_idw(train, test)
  p_nn      <- pred_nn(train, test)
  p_ens     <- pred_ensemble(train, test)
  p_krig    <- pred_kriging(train, test)
 
  rmse_cv[k, ] <- c(rmse(test$logzinc, p_closest),
                    rmse(test$logzinc, p_idw),
                    rmse(test$logzinc, p_nn),
                    rmse(test$logzinc, p_ens),
                    rmse(test$logzinc, p_krig))
}
 
round(rmse_cv, 4)

##       Closest    IDW     NN Ensemble Kriging
## fold1  0.6327 0.5496 0.4958   0.4702  0.4281
## fold2  0.5987 0.5828 0.4943   0.4718  0.4183
## fold3  0.6328 0.7843 0.6516   0.6330  0.5055
## fold4  0.4874 0.6333 0.5316   0.4907  0.4198
## fold5  0.5593 0.6293 0.4791   0.4768  0.3840

hasil <- data.frame(metode    = metode,
                    rmse_rata = colMeans(rmse_cv),
                    rmse_sd   = apply(rmse_cv, 2, sd)) %>%
  arrange(rmse_rata)
hasil

##            metode rmse_rata    rmse_sd
## Kriging   Kriging 0.4311590 0.04487027
## Ensemble Ensemble 0.5085005 0.07005186
## NN             NN 0.5304872 0.07039087
## Closest   Closest 0.5821996 0.06101586
## IDW           IDW 0.6358382 0.08991205

ggplot(hasil, aes(x = reorder(metode, rmse_rata), y = rmse_rata)) +
  geom_col(fill = "grey35", width = 0.6) +
  geom_errorbar(aes(ymin = rmse_rata - rmse_sd, ymax = rmse_rata + rmse_sd),
                width = 0.15) +
  geom_text(aes(label = round(rmse_rata, 3)), vjust = -0.6, hjust = 1.4) +
  labs(x = NULL, y = "RMSE rata-rata 5-fold CV (skala log zinc)") +
  theme_bw()

Berdasarkan hasil evaluasi 5-fold Cross Validation pada data meuse (skala log zinc), diperoleh nilai rata-rata Root Mean Square Error (RMSE) dan simpangan baku untuk masing-masing metode interpolasi sebagai berikut:

• Kriging (RMSE Rata-rata: 0.4312, SD: 0.0449)
  Interpretasi: Metode Kriging menghasilkan nilai RMSE rata-rata paling kecil dibanding seluruh metode lainnya. Hal ini menunjukkan bahwa Kriging memberikan tingkat akurasi prediksi yang paling tinggi dan stabil (memiliki deviasi standar terendah). Keunggulan ini terjadi karena Kriging memperhitungkan struktur spasial secara geostatistik melalui analisis variogram.

• Ensemble (RMSE Rata-rata: 0.5085, SD: 0.0701)
  Interpretasi: Metode Ensemble menduduki peringkat kedua terbaik. Penggabungan beberapa model terbukti mampu memperbaiki performa dibandingkan sebagian besar metode tunggal (seperti NN, Closest, dan IDW), memberikan alternatif prediktor yang cukup handal dan seimbang.

• Nearest Neighbor / NN (RMSE Rata-rata: 0.5305, SD: 0.0704)
  Interpretasi: Metode NN berada di posisi menengah, menunjukkan performa yang lebih baik daripada IDW dan Closest Observation, namun masih kalah akurat dibandingkan pendekatan geostatistik (Kriging) dan Ensemble.

• Closest Observation (RMSE Rata-rata: 0.5822, SD: 0.0610)
  Interpretasi: Metode ini hanya mengambil pengamatan terdekat tanpa mempertimbangkan tren jarak bobot yang lebih luas atau struktur spasial, sehingga menghasilkan error yang lebih tinggi.

• IDW / Inverse Distance Weighting (RMSE Rata-rata: 0.6358, SD: 0.0891)
  Interpretasi: IDW mencatatkan nilai RMSE rata-rata tertinggi dan variasi error (SD) yang paling besar di antara kelima metode. Hal ini mengindikasikan bahwa untuk sebaran data meuse ini, IDW kurang optimal dalam menangkap variasi spasial kadar zinc dibandingkan metode lainnya.

Kesimpulan: Secara keseluruhan, Kriging adalah metode interpolasi terbaik untuk dataset meuse pada pengujian ini karena menghasilkan error terkecil dan paling konsisten di setiap fold pengujian.