1 Pendahuluan

Interpolasi titik adalah proses menduga nilai suatu variabel pada lokasi yang tidak diamati berdasarkan nilai pada lokasi-lokasi yang diamati. Pada praktikum ini digunakan data meuse (kandungan logam berat pada tanah di sekitar Sungai Meuse, Belanda) dengan variabel yang diinterpolasi adalah kadar seng (zinc, ppm).

Metode interpolasi yang dibandingkan:

  1. Closest Observation (proximity / Thiessen polygon): nilai prediksi = nilai titik terdekat.
  2. Nearest Neighbour (NN): rata-rata dari nmax tetangga terdekat.
  3. Inverse Distance Weighting (IDW): rata-rata berbobot, bobot = 1 / jarak^p.
  4. Ensemble: gabungan dari beberapa metode (NN, IDW, dan Kriging).
  5. Kriging (Ordinary Kriging): prediksi geostatistik berdasarkan variogram.

Evaluasi dilakukan dengan k-fold Cross Validation dan ukuran akurasi RMSE (Root Mean Square Error).

2 Persiapan

2.1 Library

library(sp)       
library(sf)       
library(gstat)    
library(ggplot2)  
library(knitr)

2.2 Data meuse

data(meuse)
data(meuse.grid)

head(meuse)
str(meuse[, c("x", "y", "zinc")])
## 'data.frame':    155 obs. of  3 variables:
##  $ x   : num  181072 181025 181165 181298 181307 ...
##  $ y   : num  333611 333558 333537 333484 333330 ...
##  $ zinc: num  1022 1141 640 257 269 ...
summary(meuse$zinc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   113.0   198.0   326.0   469.7   674.5  1839.0
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992,
                     remove = FALSE, agr = "constant")
grid_sf  <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992,
                     remove = FALSE, agr = "constant")

nrow(meuse_sf)
## [1] 155
nrow(grid_sf)
## [1] 3103

2.3 Sebaran Titik dan Nilai zinc

ggplot() +
  geom_tile(data = meuse.grid, aes(x, y), fill = "grey90") +
  geom_point(data = meuse, aes(x, y, color = zinc, size = zinc)) +
  scale_color_viridis_c(option = "C") +
  coord_equal() +
  labs(title = "Sebaran titik pengamatan `zinc` (meuse)",
       x = "X (m)", y = "Y (m)", color = "zinc (ppm)", size = "zinc (ppm)") +
  theme_minimal()

ggplot(meuse, aes(zinc)) +
  geom_histogram(bins = 25, fill = "steelblue", color = "white") +
  labs(title = "Distribusi zinc", x = "zinc (ppm)", y = "Frekuensi") +
  theme_minimal()

3 Fungsi Prediksi untuk Tiap Metode

Setiap fungsi menerima data train dan data test (objek sf), lalu mengembalikan vektor prediksi zinc di lokasi test. Dengan format yang seragam ini, semua metode dapat dipakai baik untuk cross validation maupun untuk pemetaan grid.

3.1 1. Closest Observation

Setiap lokasi diberi nilai dari titik pengamatan yang jaraknya paling dekat (sama dengan Voronoi/Thiessen polygon).

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

3.2 2. Nearest Neighbour (NN)

Rata-rata dari nmax tetangga terdekat (tanpa pembobotan jarak, idp = 0).

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

3.3 3. Inverse Distance Weighting (IDW)

Bobot berbanding terbalik dengan jarak pangkat idp.

pred_idw <- function(train, test, idp = 2, nmax = Inf) {
  idw(zinc ~ 1, locations = train, newdata = test,
      idp = idp, nmax = nmax, debug.level = 0)$var1.pred
}

3.4 4. Kriging (Ordinary Kriging)

Langkah: (1) hitung variogram empiris, (2) fit model variogram, (3) lakukan kriging.

3.4.1 Variogram pada seluruh data

v_emp <- variogram(zinc ~ 1, meuse_sf)

init_model <- vgm(psill  = max(v_emp$gamma) * 0.8,
                  model  = "Sph",
                  range  = max(v_emp$dist) / 3,
                  nugget = min(v_emp$gamma))

v_fit <- fit.variogram(v_emp, init_model)
v_fit
plot(v_emp, v_fit, main = "Variogram empiris dan model Spherical (zinc)")

Model variogram (nugget, sill, range) selanjutnya dipakai untuk Ordinary Kriging.

3.4.2 Fungsi prediksi kriging

Pada saat cross validation, variogram difit ulang pada data train di setiap fold agar tidak terjadi kebocoran informasi dari data test.

fit_vgm <- function(data) {
  v    <- variogram(zinc ~ 1, data)
  init <- vgm(psill  = max(v$gamma) * 0.8,
              model  = "Sph",
              range  = max(v$dist) / 3,
              nugget = min(v$gamma))
  m <- tryCatch(suppressWarnings(fit.variogram(v, init)),
                error = function(e) init)
  if (isTRUE(attr(m, "singular"))) m <- init   # jaga-jaga jika fit gagal
  m
}

pred_ok <- function(train, test) {
  m <- fit_vgm(train)
  krige(zinc ~ 1, locations = train, newdata = test,
        model = m, debug.level = 0)$var1.pred
}

4 Cross Validation

Digunakan 5-fold cross validation: data dibagi acak menjadi 5 bagian; setiap bagian bergantian menjadi data test (diprediksi) sedangkan 4 bagian lainnya menjadi data train.

set.seed(123)
k     <- 5
folds <- sample(rep(1:k, length.out = nrow(meuse_sf)))
table(folds)
## folds
##  1  2  3  4  5 
## 31 31 31 31 31

4.1 Memilih parameter (opsional): nilai nmax pada NN dan idp pada IDW

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))

cv_generic <- function(fun, ...) {
  pred <- numeric(nrow(meuse_sf))
  for (i in 1:k) {
    tr <- meuse_sf[folds != i, ]
    te <- meuse_sf[folds == i, ]
    pred[folds == i] <- fun(tr, te, ...)
  }
  rmse(meuse_sf$zinc, pred)
}

tune_nn  <- data.frame(nmax = 2:10)
tune_nn$RMSE <- sapply(tune_nn$nmax, function(n) cv_generic(pred_nn, nmax = n))

tune_idw <- data.frame(idp = seq(0.5, 4, by = 0.5))
tune_idw$RMSE <- sapply(tune_idw$idp, function(p) cv_generic(pred_idw, idp = p))

kable(tune_nn,  digits = 2, caption = "RMSE NN untuk berbagai nmax")
RMSE NN untuk berbagai nmax
nmax RMSE
2 278.01
3 261.20
4 266.29
5 267.83
6 266.07
7 269.53
8 268.30
9 274.09
10 277.89 
kable(tune_idw, digits = 2, caption = "RMSE IDW untuk berbagai idp")
RMSE IDW untuk berbagai idp
idp RMSE
0.5 354.79
1.0 330.75
1.5 301.39
2.0 278.16
2.5 267.13
3.0 266.29
3.5 270.63
4.0 276.76 
nmax_opt <- tune_nn$nmax[which.min(tune_nn$RMSE)]
idp_opt  <- tune_idw$idp[which.min(tune_idw$RMSE)]
c(nmax_opt = nmax_opt, idp_opt = idp_opt)
## nmax_opt  idp_opt 
##        3        3

Parameter terbaik (nmax_opt=3 dan idp_opt=3.0) dipakai pada perbandingan akhir.

Catatan: parameter dipilih dari data yang sama dengan evaluasi sehingga RMSE NN dan IDW sedikit bersifat optimistik. Untuk evaluasi yang lebih ketat dapat digunakan nested cross validation.

4.2 Prediksi Cross Validation Semua Metode

cv_pred <- data.frame(obs = meuse_sf$zinc,
                      fold = folds,
                      closest = NA_real_, nn = NA_real_,
                      idw = NA_real_, kriging = NA_real_)

for (i in 1:k) {
  tr  <- meuse_sf[folds != i, ]
  te  <- meuse_sf[folds == i, ]
  sel <- folds == i

  cv_pred$closest[sel] <- pred_closest(tr, te)
  cv_pred$nn[sel]      <- pred_nn(tr, te, nmax = nmax_opt)
  cv_pred$idw[sel]     <- pred_idw(tr, te, idp = idp_opt)
  cv_pred$kriging[sel] <- pred_ok(tr, te)
}

head(cv_pred)

4.3 Ensemble

Ensemble menggabungkan prediksi beberapa metode. Di sini digunakan 3 metode: NN, IDW, dan Kriging, dengan dua skema:

  • Ensemble rata-rata (bobot sama: 1/3 tiap metode).
  • Ensemble berbobot RMSE: bobot sebanding dengan 1/RMSE (metode yang lebih akurat mendapat bobot lebih besar).
base_methods <- c("nn", "idw", "kriging")

# RMSE tiap metode dasar
rmse_base <- sapply(base_methods, function(m) rmse(cv_pred$obs, cv_pred[[m]]))

# Bobot = 1/RMSE, dinormalkan agar jumlahnya 1
w <- (1 / rmse_base) / sum(1 / rmse_base)
w
##        nn       idw   kriging 
## 0.3225264 0.3163590 0.3611145
cv_pred$ens_mean <- rowMeans(cv_pred[, base_methods])
cv_pred$ens_wt   <- as.vector(as.matrix(cv_pred[, base_methods]) %*% w)

Catatan: bobot dihitung dari RMSE cross validation yang sama, sehingga RMSE ens_wt sedikit optimistik dibanding ens_mean.

4.4 Perbandingan RMSE

metode <- c(closest = "Closest Observation",
            nn      = "Nearest Neighbour",
            idw     = "IDW",
            kriging = "Ordinary Kriging",
            ens_mean = "Ensemble (rata-rata NN+IDW+Kriging)",
            ens_wt   = "Ensemble (berbobot 1/RMSE)")

hasil <- data.frame(
  Metode = unname(metode),
  RMSE   = sapply(names(metode), function(m) rmse(cv_pred$obs, cv_pred[[m]])),
  MAE    = sapply(names(metode), function(m) mean(abs(cv_pred$obs - cv_pred[[m]])))
)
hasil <- hasil[order(hasil$RMSE), ]
hasil$Peringkat <- seq_len(nrow(hasil))
rownames(hasil) <- NULL

kable(hasil, digits = 2, caption = "Perbandingan RMSE (5-fold CV) seluruh metode interpolasi")
Perbandingan RMSE (5-fold CV) seluruh metode interpolasi
Metode RMSE MAE Peringkat
Ordinary Kriging 233.29 161.81 1
Ensemble (berbobot 1/RMSE) 245.53 170.81 2
Ensemble (rata-rata NN+IDW+Kriging) 246.30 171.32 3
Nearest Neighbour 261.20 177.29 4
IDW 266.29 186.02 5
Closest Observation 351.95 220.69 6 
ggplot(hasil, aes(x = reorder(Metode, -RMSE), y = RMSE, fill = RMSE)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 1)), hjust = -0.1, size = 3.5) +
  coord_flip() +
  scale_fill_viridis_c(option = "D", direction = -1) +
  expand_limits(y = max(hasil$RMSE) * 1.1) +
  labs(title = "RMSE hasil 5-fold Cross Validation",
       x = NULL, y = "RMSE (ppm)") +
  theme_minimal()

library(tidyr)  

long <- pivot_longer(cv_pred, cols = names(metode),
                     names_to = "metode", values_to = "pred")
long$metode <- factor(long$metode, levels = names(metode), labels = metode)

ggplot(long, aes(obs, pred)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  facet_wrap(~ metode) +
  labs(title = "Nilai observasi vs prediksi (Cross Validation)",
       x = "Observasi zinc (ppm)", y = "Prediksi zinc (ppm)") +
  theme_minimal()

5 Peta Hasil Interpolasi

Setelah dievaluasi, setiap metode diterapkan pada seluruh data untuk memprediksi zinc di seluruh meuse.grid.

pred_grid <- data.frame(
  closest = pred_closest(meuse_sf, grid_sf),
  nn      = pred_nn(meuse_sf, grid_sf, nmax = nmax_opt),
  idw     = pred_idw(meuse_sf, grid_sf, idp = idp_opt),
  kriging = krige(zinc ~ 1, meuse_sf, grid_sf, model = v_fit,
                  debug.level = 0)$var1.pred
)
pred_grid$ens_mean <- rowMeans(pred_grid[, base_methods])
pred_grid$ens_wt   <- as.vector(as.matrix(pred_grid[, base_methods]) %*% w)

grid_long <- cbind(meuse.grid[, c("x", "y")], pred_grid)
grid_long <- pivot_longer(grid_long, cols = names(metode),
                          names_to = "metode", values_to = "zinc")
grid_long$metode <- factor(grid_long$metode, levels = names(metode), labels = metode)
ggplot(grid_long, aes(x, y, fill = zinc)) +
  geom_raster() +
  geom_point(data = meuse, aes(x, y), inherit.aes = FALSE,
             size = 0.4, color = "black", alpha = 0.5) +
  facet_wrap(~ metode, ncol = 3) +
  coord_equal() +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Peta interpolasi zinc (meuse.grid)",
       x = "X (m)", y = "Y (m)", fill = "zinc (ppm)") +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.ticks = element_blank())

5.1 Peta Varians Kriging

Kelebihan Kriging dibanding metode lain adalah tersedianya ukuran ketidakpastian (varians prediksi).

ok_full <- krige(zinc ~ 1, meuse_sf, grid_sf, model = v_fit, debug.level = 0)

ggplot(data.frame(meuse.grid[, c("x", "y")], var = ok_full$var1.var),
       aes(x, y, fill = var)) +
  geom_raster() +
  geom_point(data = meuse, aes(x, y), inherit.aes = FALSE, size = 0.5) +
  coord_equal() +
  scale_fill_viridis_c(option = "B") +
  labs(title = "Varians prediksi Ordinary Kriging",
       x = "X (m)", y = "Y (m)", fill = "Varians") +
  theme_minimal()

6 Kesimpulan

Berdasarkan 5-fold cross validation, metode dengan RMSE terkecil adalah Ordinary Kriging (RMSE = 233.29 ppm), sedangkan metode dengan RMSE terbesar adalah Closest Observation (RMSE = 351.95 ppm).