Interpolasi Titik pada Data meuse

Data

library(sp)
library(gstat)
library(knitr)

data(meuse)
data(meuse.grid)

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

coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y
gridded(meuse.grid) <- TRUE

bubble(meuse, "zinc", col = "steelblue", main = "Sebaran Zinc (ppm) pada Titik Sampel")

Fungsi Metode Interpolasi

Setiap metode dibuat sebagai fungsi dengan input train (data latih) dan test (lokasi yang diprediksi). Dengan begitu fungsi yang sama dapat dipakai untuk membuat peta prediksi maupun untuk cross validation.

Closest Observation

Setiap lokasi diberi nilai dari satu titik pengamatan yang jaraknya paling dekat (setara poligon Voronoi/Thiessen).

pred_closest <- function(train, test) {
  d <- spDists(coordinates(test), coordinates(train))
  idx <- apply(d, 1, which.min)
  train$zinc[idx]
}

Nearest Neighbour (NN)

Nilai prediksi adalah rata-rata dari k = 5 titik terdekat (bobot sama, sehingga idp = 0).

K_NN <- 5

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

Inverse Distance Weighting (IDW)

Titik yang lebih dekat diberi bobot lebih besar, dengan bobot \(w_i \propto 1/d_i^{p}\). Di sini digunakan pangkat \(p = 2\).

P_IDW <- 2

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

Kriging

Ordinary Kriging membutuhkan model variogram. Berikut variogram empiris dan model yang difit pada seluruh data.

v_emp <- variogram(zinc ~ 1, meuse)

v_fit <- fit.variogram(
  v_emp,
  model = vgm(psill = var(meuse$zinc) * 0.7, model = "Sph",
              range = 900, nugget = var(meuse$zinc) * 0.3)
)
v_fit

##   model     psill    range
## 1   Nug  24799.13   0.0000
## 2   Sph 134743.94 830.9164

plot(v_emp, model = v_fit, main = "Variogram Empiris dan Model Spherical")

Pada cross validation, variogram dihitung ulang di setiap fold hanya dari data latih, agar tidak terjadi kebocoran informasi dari data uji.

pred_krige <- function(train, test) {
  v  <- variogram(zinc ~ 1, train)
  m  <- tryCatch(
    suppressWarnings(fit.variogram(
      v, vgm(psill = var(train$zinc) * 0.7, "Sph", range = 900,
             nugget = var(train$zinc) * 0.3))),
    error = function(e) NULL
  )
  if (is.null(m)) m <- v_fit   # cadangan bila fitting gagal
  krige(zinc ~ 1, locations = train, newdata = test,
        model = m, debug.level = 0)$var1.pred
}

Ensemble

Ensemble menggabungkan prediksi beberapa metode. Pada tugas ini digunakan tiga metode: NN, IDW, dan Kriging. Terdapat dua versi:

Ensemble rata-rata: rata-rata sederhana ketiga prediksi.
Ensemble berbobot: bobot sebanding dengan \(1/\mathrm{RMSE}^2\) tiap metode (bobot dihitung dari hasil cross validation, lihat bagian berikutnya).

pred_ens_mean <- function(train, test) {
  rowMeans(cbind(pred_nn(train, test),
                 pred_idw(train, test),
                 pred_krige(train, test)))
}

Peta Hasil Interpolasi

Seluruh data digunakan sebagai data latih untuk memprediksi meuse.grid.

grid_pred <- meuse.grid
grid_pred$Closest <- pred_closest(meuse, meuse.grid)
grid_pred$NN      <- pred_nn(meuse, meuse.grid)
grid_pred$IDW     <- pred_idw(meuse, meuse.grid)
grid_pred$Kriging <- pred_krige(meuse, meuse.grid)
grid_pred$Ensemble <- (grid_pred$NN + grid_pred$IDW + grid_pred$Kriging) / 3

rng <- range(meuse$zinc)

spplot(grid_pred,
       zcol = c("Closest", "NN", "IDW", "Kriging", "Ensemble"),
       col.regions = hcl.colors(50, "YlOrRd", rev = TRUE),
       at = seq(rng[1], rng[2], length.out = 51),
       main = "Interpolasi Zinc (ppm) Berdasarkan Berbagai Metode")

* Closest Observation menghasilkan permukaan berblok-blok (poligon Voronoi) dan tidak mulus. * NN (k = 5) lebih halus, tetapi masih terlihat berundak. * IDW mulus namun cenderung membentuk pola “mata sapi” di sekitar titik sampel. * Kriging menghasilkan permukaan yang halus dan memanfaatkan struktur autokorelasi spasial. * Ensemble merupakan kompromi dari ketiga metode.

Cross Validation

Digunakan 5-fold cross validation. Titik dibagi acak menjadi 5 kelompok. Secara bergantian satu kelompok menjadi data uji dan empat kelompok lainnya menjadi data latih. Fold yang sama dipakai untuk semua metode agar perbandingannya adil.

\[RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(\hat{z}_i - z_i)^2}\]

set.seed(123)
n_fold <- 5
folds  <- sample(rep(1:n_fold, length.out = nrow(meuse)))

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

cv_pred <- data.frame(obs = meuse$zinc,
                      Closest = NA_real_, NN = NA_real_,
                      IDW = NA_real_, Kriging = NA_real_)

for (f in 1:n_fold) {
  train <- meuse[folds != f, ]
  test  <- meuse[folds == f, ]

  cv_pred$Closest[folds == f] <- pred_closest(train, test)
  cv_pred$NN[folds == f]      <- pred_nn(train, test)
  cv_pred$IDW[folds == f]     <- pred_idw(train, test)
  cv_pred$Kriging[folds == f] <- pred_krige(train, test)
}

# Ensemble rata-rata (NN, IDW, Kriging)
cv_pred$Ensemble_Rata2 <- rowMeans(cv_pred[, c("NN", "IDW", "Kriging")])

head(cv_pred)

##    obs Closest    NN      IDW  Kriging Ensemble_Rata2
## 1 1022    1141 558.0 807.7687 884.0847       749.9511
## 2 1141    1022 684.2 754.4514 882.2215       773.6243
## 3  640    1022 636.8 689.1782 728.1113       684.6965
## 4  257     269 636.8 464.1078 496.1365       532.3481
## 5  269     281 237.8 347.7180 255.6693       280.3958
## 6  281     269 253.2 345.0371 278.9348       292.3906

Ensemble berbobot dibentuk dari RMSE ketiga metode penyusunnya:

rmse_base <- sapply(c("NN", "IDW", "Kriging"),
                    function(m) rmse(cv_pred$obs, cv_pred[[m]]))

w <- (1 / rmse_base^2) / sum(1 / rmse_base^2)
round(w, 3)

##      NN     IDW Kriging 
##   0.308   0.286   0.406

cv_pred$Ensemble_Berbobot <- as.numeric(
  as.matrix(cv_pred[, c("NN", "IDW", "Kriging")]) %*% w
)

Perbandingan RMSE

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

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

kable(hasil, digits = 2, row.names = FALSE,
      caption = "Perbandingan performa metode interpolasi (5-fold CV)")

Perbandingan performa metode interpolasi (5-fold CV)
Metode	RMSE	MAE	R2	Peringkat
Kriging	233.29	161.81	0.60	1
Ensemble_Berbobot	248.41	174.44	0.56	2
Ensemble_Rata2	251.43	177.06	0.55	3
NN	267.83	182.46	0.47	4
IDW	278.16	207.16	0.49	5
Closest	351.95	220.69	0.29	6

hs <- hasil[order(hasil$RMSE, decreasing = TRUE), ]
par(mar = c(4, 10, 3, 1))
bp <- barplot(hs$RMSE, names.arg = hs$Metode, horiz = TRUE, las = 1,
              col = ifelse(hs$RMSE == min(hs$RMSE), "darkgreen", "grey70"),
              xlab = "RMSE (ppm)", main = "RMSE Cross Validation per Metode",
              xlim = c(0, max(hs$RMSE) * 1.15))
text(hs$RMSE, bp, labels = round(hs$RMSE, 1), pos = 4, cex = 0.9)

Plot nilai observasi vs prediksi:

par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (m in metode) {
  plot(cv_pred$obs, cv_pred[[m]], pch = 16, col = rgb(0, 0, 1, 0.5),
       xlab = "Observasi", ylab = "Prediksi",
       main = sprintf("%s\nRMSE = %.1f", m, rmse(cv_pred$obs, cv_pred[[m]])))
  abline(0, 1, col = "red", lty = 2)
}

par(mfrow = c(1, 1))

Kesimpulan

best  <- hasil$Metode[1]
worst <- hasil$Metode[nrow(hasil)]
cat(sprintf(
"Berdasarkan 5-fold cross validation, metode dengan **RMSE terkecil** adalah **%s** (RMSE = %.2f ppm), sedangkan metode dengan RMSE terbesar adalah **%s** (RMSE = %.2f ppm).",
best, hasil$RMSE[1], worst, hasil$RMSE[nrow(hasil)]))

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