Data

data(meuse)
data(meuse.grid)

meuse      <- st_as_sf(meuse,      coords = c("x", "y"), crs = 28992)
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Peubah yang diinterpolasi
meuse$vble <- log(meuse$zinc)

summary(meuse$zinc)

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

summary(meuse$vble)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.727   5.288   5.787   5.886   6.514   7.517

skew <- function(x) mean((x - mean(x))^3) / sd(x)^3
par(mfrow = c(1, 2))
hist(meuse$zinc, breaks = 20, col = "grey80", border = "white",
     main = "zinc", xlab = "zinc (ppm)")
hist(meuse$vble, breaks = 20, col = "steelblue", border = "white",
     main = "log(zinc)", xlab = "log(zinc)")

par(mfrow = c(1, 1))

Data terdiri atas 155 titik pengamatan, dan lokasi prediksi berupa 3103 titik grid (meuse.grid). Koefisien skewness zinc sebesar 1.46, turun menjadi 0.32 setelah transformasi log. Transformasi log membuat sebaran lebih simetris sehingga asumsi stasioneritas pada kriging lebih masuk akal dan variogram tidak didominasi nilai ekstrem. Seluruh RMSE di bawah berada pada skala log(zinc).

Sebaran lokasi pengamatan dan lokasi prediksi

ggplot() +
  geom_sf(data = meuse.grid, colour = "grey85", size = 0.3) +
  geom_sf(data = meuse, aes(colour = vble), size = 2) +
  scale_colour_viridis_c(name = "log(zinc)") +
  labs(title = "Titik pengamatan (berwarna) dan grid prediksi (abu-abu)") +
  theme_bw()

Titik pengamatan tidak tersebar merata, mengumpul di sepanjang tepi sungai, sedangkan grid prediksi menutup seluruh dataran banjir. Ini berarti sebagian lokasi prediksi berada jauh dari titik pengamatan terdekat, dan di sanalah perbedaan perilaku antarmetode akan paling terlihat. Nilai log(zinc) yang tinggi cenderung mengelompok, indikasi awal adanya autokorelasi spasial.

# Wadah hasil prediksi pada grid
xy   <- st_coordinates(meuse.grid)
pred <- data.frame(x = xy[, 1], y = xy[, 2])

# Fungsi bantu pemetaan
plot_pred <- function(df, kolom, judul, legenda = "log(zinc)") {
  ggplot(df, aes(x, y, fill = .data[[kolom]])) +
    geom_raster() +
    scale_fill_viridis_c(name = legenda) +
    coord_equal() +
    labs(title = judul, x = NULL, y = NULL) +
    theme_bw()
}

Closest Observation (Voronoi)

Setiap lokasi diberi nilai titik pengamatan terdekat. Wilayah yang dipengaruhi satu titik membentuk poligon Voronoi.

# Batas area
area <- st_convex_hull(st_union(c(st_geometry(meuse), st_geometry(meuse.grid))))
area <- st_buffer(area, 100)

v <- terra::voronoi(terra::vect(meuse), bnd = terra::vect(area))
v <- sf::st_as_sf(v)

ggplot() +
  geom_sf(data = v, aes(fill = vble), colour = "white", linewidth = 0.2) +
  geom_sf(data = meuse, size = 0.6) +
  scale_fill_viridis_c(name = "log(zinc)") +
  labs(title = "Diagram Voronoi") +
  theme_bw()

Setiap poligon memiliki satu nilai konstan dari titik pengamatan di dalamnya. Poligon di area padat titik berukuran kecil, sedangkan di area jarang titik berukuran besar, sehingga satu pengamatan “menguasai” wilayah luas.

idx <- st_nearest_feature(meuse.grid, meuse)
pred$closest <- meuse$vble[idx]

plot_pred(pred, "closest", "Closest Observation")

Peta prediksi berbentuk blok dengan perubahan nilai yang tiba-tiba di batas poligon. Metode ini sederhana dan mempertahankan nilai asli data, tetapi tidak menghasilkan permukaan yang mulus dan tidak memakai informasi dari titik lain selain yang terdekat, sehingga rentan terhadap noise pada satu titik.

Inverse Distance Weighting (IDW)

g_idw <- gstat(formula = vble ~ 1, locations = meuse,
               nmax = nrow(meuse), set = list(idp = 1))
pred$idw <- predict(g_idw, meuse.grid, debug.level = 0)$var1.pred

plot_pred(pred, "idw", "IDW (idp = 1, semua titik)")

Permukaan lebih mulus dibanding Closest Observation karena prediksi merupakan rata-rata berbobot dari banyak titik. Karena idp = 1 relatif kecil dan semua titik ikut dihitung, titik yang jauh tetap memberi pengaruh, sehingga prediksi cenderung tertarik ke nilai rata-rata global (over-smoothing), terutama di area yang jauh dari pengamatan. Menaikkan idp atau membatasi nmax akan membuat prediksi lebih lokal.

Nearest Neighbors (k = 5)

g_nn <- gstat(formula = vble ~ 1, locations = meuse,
              nmax = 5, set = list(idp = 0))
pred$nn <- predict(g_nn, meuse.grid, debug.level = 0)$var1.pred

plot_pred(pred, "nn", "Nearest Neighbors (k = 5)")

Hasilnya lebih halus dibanding Closest Observation karena merata-ratakan 5 titik, tetapi tetap tampak “bertingkat” saat himpunan 5 tetangga berganti. Pemilihan k menentukan keseimbangan bias-varians: k kecil membuat prediksi mengikuti data lokal (varians tinggi), k besar membuatnya lebih halus (bias lebih tinggi). Pada latihan ini k tidak di-tuning.

Variogram empiris

vc <- variogram(vble ~ 1, meuse, cloud = TRUE)
plot(vc)

v_emp <- variogram(vble ~ 1, data = meuse)
plot(v_emp)

Variogram cloud menampilkan semivarians seluruh pasangan titik dan sangat berisik, sehingga sulit dibaca langsung. Variogram empiris merangkumnya per kelas jarak. Semivarians rendah pada jarak dekat dan naik seiring jarak, yang menandakan adanya struktur ketergantungan spasial, titik berdekatan lebih mirip daripada titik berjauhan.

Model variogram

show.vgms(par.strip.text = list(cex = 0.75))

vinitial <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)
plot(v_emp, vinitial)

fv <- fit.variogram(object = v_emp, model = vinitial)
fv

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

plot(v_emp, fv, cex = 1.5)

nugget <- fv$psill[1]
sill   <- sum(fv$psill)
rng    <- fv$range[2]
rasio  <- nugget / sill

Model spherical hasil fitting memiliki nugget 0.051, sill 0.641, dan range 897 meter.

• Nugget 0.051 adalah variasi mikro-spasial dan galat pengukuran.
• Sill 0.641 adalah total variasi spasial.
• Range 897 m berarti pasangan titik yang berjarak lebih dari itu dianggap tidak lagi berkorelasi.
• Rasio nugget terhadap sill adalah 7.9%, termasuk ketergantungan spasial kuat (kriteria umum: < 25% kuat, 25-75% sedang, > 75% lemah).

Kurva model mengikuti pola titik variogram empiris dengan baik secara visual, tetapi kecocokan ini hanya berdasarkan inspeksi grafis dan nilai SSErr fitting, belum diuji lebih jauh (misalnya terhadap model eksponensial).

Prediksi Kriging

g_ok  <- gstat(formula = vble ~ 1, data = meuse, model = fv)
kpred <- predict(g_ok, meuse.grid, debug.level = 0)

pred$kriging     <- kpred$var1.pred
pred$kriging_var <- kpred$var1.var

plot_pred(pred, "kriging", "Ordinary Kriging prediksi log(zinc)")

Permukaan prediksi kriging paling mulus di antara metode sejauh ini, karena bobot tiap titik ditentukan dari struktur korelasi spasial (variogram), bukan hanya jarak. Titik-titik yang saling mengelompok juga tidak “dihitung ganda” sebagaimana pada IDW.

plot_pred(pred, "kriging_var", "Ordinary Kriging: varians prediksi", "varians")

Varians prediksi (0.085 sampai 0.5) kecil di sekitar titik pengamatan dan membesar di area yang jauh dari data, terutama di tepi grid. Inilah keunggulan utama kriging dibanding metode deterministik: kriging memberi ukuran ketidakpastian, bukan hanya prediksi titik.

Perbandingan seluruh peta

lab <- c(closest      = "Closest Observation",
         idw          = "IDW",
         nn           = "Nearest Neighbors (k = 5)",
         kriging      = "Ordinary Kriging",
         ensemble     = "Ensembel (bobot sama)")

pred_long <- pred %>%
  select(x, y, all_of(names(lab))) %>%
  pivot_longer(-c(x, y), names_to = "metode", values_to = "pred") %>%
  mutate(metode = factor(metode, levels = names(lab), labels = lab))

ggplot(pred_long, aes(x, y, fill = pred)) +
  geom_raster() +
  facet_wrap(~ metode, ncol = 2) +
  scale_fill_viridis_c(name = "log(zinc)") +
  coord_equal() +
  labs(x = NULL, y = NULL) +
  theme_bw()

Dengan skala warna yang sama, pola umum sebaran zinc (konsentrasi tinggi dekat sungai) tertangkap oleh semua metode. Perbedaannya ada pada kehalusan permukaan dan kekuatan kontras antarwilayah. Peta saja belum cukup untuk menentukan metode terbaik, sehingga diperlukan evaluasi kuantitatif.

Rancangan

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

Ada dua versi ensembel:

• Ensembel bobot sama: rata-rata sederhana keempat metode.
• Ensembel bobot 1/RMSE: bobot proporsional dengan kebalikan RMSE tiap metode, sehingga metode yang lebih akurat mendapat bobot lebih besar. Berbeda dengan contoh pada materi (yang menghitung bobot dari RMSE data uji pada lipatan yang sama, sehingga data uji ikut “bocor” ke proses pembobotan), bobot di sini dihitung hanya dari data latih lewat cross-validation bersarang (inner 5-fold).

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

# Fitting variogram pada data latih
fit_vgm <- function(train) {
  v_tr <- variogram(vble ~ 1, train)
  fit  <- tryCatch(
    suppressWarnings(fit.variogram(v_tr, vinitial)),
    error = function(e) NULL
  )
  if (is.null(fit)) vinitial else fit
}

# Prediksi keempat metode
predict_methods <- function(train, test) {
  p_closest <- train$vble[st_nearest_feature(test, train)]

  p_idw <- predict(
    gstat(formula = vble ~ 1, locations = train,
          nmax = nrow(train), set = list(idp = 1)),
    test, debug.level = 0)$var1.pred

  p_nn <- predict(
    gstat(formula = vble ~ 1, locations = train,
          nmax = 5, set = list(idp = 0)),
    test, debug.level = 0)$var1.pred

  p_ok <- predict(
    gstat(formula = vble ~ 1, data = train, model = fit_vgm(train)),
    test, debug.level = 0)$var1.pred

  data.frame(closest = p_closest, idw = p_idw, nn = p_nn, kriging = p_ok)
}

# Bobot 1/RMSE dari data latih saja
inv_rmse_weights <- function(train, k = 5) {
  kf_in <- dismo::kfold(nrow(train), k = k)
  P <- matrix(NA, nrow(train), 4)
  for (j in 1:k) {
    P[kf_in == j, ] <- as.matrix(
      predict_methods(train[kf_in != j, ], train[kf_in == j, ])
    )
  }
  w <- 1 / apply(P, 2, function(p) rmse(train$vble, p))
  w / sum(w)
}

Menjalankan cross-validation

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

res_rmse <- data.frame(fold = 1:K, closest = NA, idw = NA, nn = NA,
                       kriging = NA, ens_equal = NA, ens_weighted = NA)
w_store  <- matrix(NA, K, 4, dimnames = list(paste0("fold", 1:K), metode))

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

  P <- predict_methods(train, test)
  w <- inv_rmse_weights(train)
  w_store[k, ] <- w

  for (m in metode) res_rmse[k, m] <- rmse(test$vble, P[[m]])
  res_rmse$ens_equal[k]    <- rmse(test$vble, rowMeans(P))
  res_rmse$ens_weighted[k] <- rmse(test$vble, as.numeric(as.matrix(P) %*% w))
}

round(res_rmse, 4)

##   fold closest    idw     nn kriging ens_equal ens_weighted
## 1    1  0.6327 0.5496 0.4958  0.4281    0.4508       0.4476
## 2    2  0.5987 0.5828 0.4943  0.4183    0.4492       0.4441
## 3    3  0.6328 0.7843 0.6516  0.5055    0.5995       0.5891
## 4    4  0.4874 0.6333 0.5316  0.4198    0.4633       0.4566
## 5    5  0.5593 0.6293 0.4791  0.3840    0.4463       0.4364

Tabel di atas memuat RMSE tiap metode pada setiap lipatan, dengan ukuran data uji sekitar 31 titik per lipatan. RMSE antarlipatan bervariasi karena sebagian lipatan kebetulan memuat titik yang lebih sulit diprediksi. Karena itu, perbandingan metode sebaiknya dilihat dari rata-rata dan konsistensi antarlipatan.

Ringkasan RMSE

lab_cv <- c(closest      = "Closest Observation",
            idw          = "IDW",
            nn           = "Nearest Neighbors (k = 5)",
            kriging      = "Ordinary Kriging",
            ens_equal    = "Ensembel (bobot sama)",
            ens_weighted = "Ensembel (bobot 1/RMSE)")

rmse_mean <- colMeans(res_rmse[, names(lab_cv)])
rmse_sd   <- apply(res_rmse[, names(lab_cv)], 2, sd)

tab <- data.frame(
  kode          = names(lab_cv),
  Metode        = unname(lab_cv),
  RMSE_rata2    = round(rmse_mean, 4),
  SD_antar_fold = round(rmse_sd, 4),
  row.names     = NULL
)
tab <- tab[order(tab$RMSE_rata2), ]
tab$Peringkat <- seq_len(nrow(tab))
rownames(tab) <- NULL
tab[, c("Peringkat", "Metode", "RMSE_rata2", "SD_antar_fold")]

##   Peringkat                    Metode RMSE_rata2 SD_antar_fold
## 1         1          Ordinary Kriging     0.4312        0.0449
## 2         2   Ensembel (bobot 1/RMSE)     0.4748        0.0644
## 3         3     Ensembel (bobot sama)     0.4818        0.0661
## 4         4 Nearest Neighbors (k = 5)     0.5305        0.0704
## 5         5       Closest Observation     0.5822        0.0610
## 6         6                       IDW     0.6358        0.0899

res_long <- res_rmse %>%
  pivot_longer(-fold, names_to = "kode", values_to = "rmse") %>%
  mutate(Metode = factor(lab_cv[kode], levels = tab$Metode))

ggplot() +
  geom_col(data = tab %>% mutate(Metode = factor(Metode, levels = Metode)),
           aes(Metode, RMSE_rata2), fill = "steelblue", alpha = 0.6) +
  geom_point(data = res_long, aes(Metode, rmse),
             position = position_jitter(width = 0.08, height = 0), size = 1.8) +
  labs(x = NULL, y = "RMSE (skala log(zinc))",
       title = "RMSE rata-rata (batang) dan RMSE tiap fold (titik)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

round(w_store, 3)

##       closest   idw    nn kriging
## fold1   0.235 0.205 0.257   0.303
## fold2   0.217 0.208 0.270   0.305
## fold3   0.241 0.204 0.257   0.298
## fold4   0.216 0.203 0.260   0.320
## fold5   0.239 0.199 0.258   0.304

round(colMeans(w_store), 3)

## closest     idw      nn kriging 
##   0.230   0.204   0.260   0.306

Interpretasi.

• Metode dengan RMSE rata-rata terendah adalah Ordinary Kriging (0.4312), sedangkan yang tertinggi adalah IDW (0.6358). Selisihnya 32.2% terhadap yang tertinggi.
• Ordinary Kriging lebih akurat dibanding metode deterministik terbaik (RMSE rata-rata 0.4312 vs 0.5305), dan mengungguli ketiga metode deterministik pada 5 dari 5 lipatan.
• Ensembel bobot sama tidak lebih baik dibanding metode tunggal terbaik (0.4818 vs 0.4312). Ensembel bobot 1/RMSE menghasilkan 0.4748.
• Bobot rata-rata pada ensembel 1/RMSE: closest = 0.23, idw = 0.204, nn = 0.26, kriging = 0.306. Karena bobot ini proporsional dengan 1/RMSE, perbedaan antarmetode yang kecil menghasilkan bobot yang hampir sama, sehingga selisih kedua versi ensembel ikut kecil.
• Simpangan baku antarlipatan (kolom SD_antar_fold) perlu dibandingkan dengan selisih antarmetode: bila selisih rata-rata antara dua metode lebih kecil dari simpangan baku tersebut, urutan keduanya sebaiknya tidak dianggap pasti.