1 Pendahuluan & Eksplorasi Data `meuse`

Dataset meuse (dari paket sp) berisi data konsentrasi empat logam berat (kadmium, tembaga, timbal, dan seng/zinc) beserta variabel lingkungan di dataran banjir Sungai Meuse, Belanda (\(n = 155\) titik sampel). Selain itu, tersedia meuse.grid berupa grid reguler beresolusi \(40\,\text{m} \times 40\,\text{m}\) sebanyak \(3.103\) sel sebagai domain target prediksi spasial.

Pada analisis ini, variabel target yang diinterpolasi adalah konsentrasi Zinc (zinc, satuan ppm). Karena sebaran zinc sangat menceng kanan (right-skewed), transformasi logaritma natural log(zinc) digunakan sebagai target pemodelan standar dalam geostatistika agar lebih mendekati asumsi normalitas dan homogenitas varians. Seluruh metode dievaluasi pada skala yang konsisten (log(zinc)), dan juga dilaporkan hasil transformasi baliknya ke skala asli (zinc ppm).

library(sp)
library(gstat)
library(mgcv)
library(deldir)
library(ggplot2)
library(dplyr)
library(tidyr)
library(viridis)
library(patchwork)
library(knitr)

data(meuse)
data(meuse.grid)

# Variabel target log(zinc)
meuse$log_zinc <- log(meuse$zinc)

# Objek SpatialPointsDataFrame untuk gstat
meuse_sp <- meuse
coordinates(meuse_sp) <- ~x + y

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

cat("Jumlah titik observasi meuse :", nrow(meuse), "\n")

## Jumlah titik observasi meuse : 155

cat("Jumlah sel grid meuse.grid   :", nrow(meuse.grid), "\n")

## Jumlah sel grid meuse.grid   : 3103

ringkasan <- data.frame(
  Statistik = c("Minimum", "Kuartil 1 (Q1)", "Median", "Mean", "Kuartil 3 (Q3)", "Maksimum", "Std. Deviasi"),
  Zinc_ppm = c(
    min(meuse$zinc), quantile(meuse$zinc, 0.25), median(meuse$zinc),
    mean(meuse$zinc), quantile(meuse$zinc, 0.75), max(meuse$zinc), sd(meuse$zinc)
  ),
  Log_Zinc = c(
    min(meuse$log_zinc), quantile(meuse$log_zinc, 0.25), median(meuse$log_zinc),
    mean(meuse$log_zinc), quantile(meuse$log_zinc, 0.75), max(meuse$log_zinc), sd(meuse$log_zinc)
  )
)

kable(ringkasan, digits = 3, caption = "Ringkasan Statistik Variabel Target Zinc dan log(Zinc)")

Ringkasan Statistik Variabel Target Zinc dan log(Zinc)
Statistik	Zinc_ppm	Log_Zinc
Minimum	113.000	4.727
Kuartil 1 (Q1)	198.000	5.288
Median	326.000	5.787
Mean	469.716	5.886
Kuartil 3 (Q3)	674.500	6.514
Maksimum	1839.000	7.517
Std. Deviasi	367.074	0.722

p_hist1 <- ggplot(meuse, aes(x = zinc)) +
  geom_histogram(bins = 22, fill = "#2c7fb8", color = "white", alpha = 0.85) +
  labs(title = "Distribusi Zinc (ppm)", x = "zinc (ppm)", y = "Frekuensi") +
  theme_minimal()

p_hist2 <- ggplot(meuse, aes(x = log_zinc)) +
  geom_histogram(bins = 22, fill = "#41ab5d", color = "white", alpha = 0.85) +
  labs(title = "Distribusi log(Zinc)", x = "log(zinc)", y = "Frekuensi") +
  theme_minimal()

p_map_pts <- ggplot() +
  geom_raster(data = meuse.grid, aes(x = x, y = y), fill = "grey90") +
  geom_point(data = meuse, aes(x = x, y = y, color = log_zinc, size = zinc), alpha = 0.85) +
  scale_color_viridis_c(option = "plasma", name = "log(zinc)") +
  scale_size_continuous(range = c(1.5, 5), name = "zinc (ppm)") +
  coord_equal() +
  labs(title = "Sebaran Spasial Titik Sampel Meuse", x = "Koordinat X (m)", y = "Koordinat Y (m)") +
  theme_minimal() +
  theme(axis.text = element_text(size = 7))

(p_hist1 / p_hist2) | p_map_pts

Terlihat pola spasial yang sangat jelas: konsentrasi zinc tinggi terkonsentrasi di sepanjang tepi Sungai Meuse (bagian barat-barat laut area studi) dan menurun seiring menjauhnya lokasi dari sungai.

2 Implementasi Metode Interpolasi Spasial

Pada bagian ini, seluruh metode interpolasi dijalankan pada data penuh untuk memetakan prediksi ke domain meuse.grid:

Closest Observation (1-NN / Poligon Voronoi)
Nearest Neighbour (NN / \(k\)-NN, \(k = 5\))
Inverse Distance Weighting (IDW, \(p = 2, \text{nmax} = 10\))
Ensemble (menggabungkan 3 metode: NN, IDW, dan Thin Plate Spline / GAM)
Ordinary Kriging

2.1 Closest Observation (1-Nearest Neighbor / Voronoi)

Metode Closest Observation menetapkan nilai prediksi di lokasi target \(s_0\) sama persis dengan nilai 1 titik observasi terdekat (\(k = 1\)):

\[\hat{z}(s_0) = z(s_i), \quad i = \arg\min_{j \in \{1,\dots,n\}} \|s_0 - s_j\|\]

Secara geometris, pendekatan ini membentuk partisi Poligon Voronoi / Thiessen yang bersifat konstan per poligon (piecewise constant) dan tidak kontinu di sepanjang batas poligon.

model_closest <- gstat(formula = log_zinc ~ 1, locations = meuse_sp, nmax = 1)
pred_closest_grid <- predict(model_closest, grid_sp, debug.level = 0)

# Hitung segmen batas poligon Voronoi menggunakan deldir untuk visualisasi
vor_res <- deldir(meuse$x, meuse$y, rw = c(
  min(meuse.grid$x) - 20, max(meuse.grid$x) + 20,
  min(meuse.grid$y) - 20, max(meuse.grid$y) + 20
))

df_grid <- as.data.frame(meuse.grid)
df_grid$closest <- pred_closest_grid$var1.pred

ggplot() +
  geom_raster(data = df_grid, aes(x = x, y = y, fill = closest)) +
  geom_segment(
    data = vor_res$dirsgs,
    aes(x = x1, y = y1, xend = x2, yend = y2),
    color = "white", linewidth = 0.25, alpha = 0.65
  ) +
  geom_point(data = meuse, aes(x = x, y = y), size = 0.7, color = "black", alpha = 0.5) +
  scale_fill_viridis_c(option = "viridis", name = "Prediksi\nlog(zinc)") +
  coord_equal() +
  labs(
    title = "Interpolasi Closest Observation (1-NN) & Batas Poligon Voronoi",
    x = "Koordinat X (m)", y = "Koordinat Y (m)"
  ) +
  theme_minimal()

2.2 Nearest Neighbour (NN / \(k\)-Nearest Neighbors)

Metode Nearest Neighbour (NN) memperhalus Closest Observation dengan mengambil rata-rata tak terboboti (equal-weight average) dari \(k\) titik observasi terdekat:

\[\hat{z}(s_0) = \frac{1}{k} \sum_{i \in N_k(s_0)} z(s_i)\]

Di paket gstat, hal ini diimplementasikan dengan nmax = k dan parameter pangkat jarak idp = 0. Di sini digunakan \(k = 5\) (nilai standar yang menjaga keseimbangan antara kehalusan lokal dan detail spasial). Kita juga dapat memeriksa sensitivitas nilai \(k\) pada bagian evaluasi.

model_nn5 <- gstat(formula = log_zinc ~ 1, locations = meuse_sp, nmax = 5, set = list(idp = 0))
pred_nn5_grid <- predict(model_nn5, grid_sp, debug.level = 0)
df_grid$nn5 <- pred_nn5_grid$var1.pred

ggplot(df_grid, aes(x = x, y = y, fill = nn5)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "Prediksi\nlog(zinc)") +
  coord_equal() +
  labs(
    title = "Interpolasi Nearest Neighbour (k = 5, Bobot Seragam)",
    x = "Koordinat X (m)", y = "Koordinat Y (m)"
  ) +
  theme_minimal()

2.3 Inverse Distance Weighting (IDW)

Metode Inverse Distance Weighting (IDW) memberi bobot yang semakin besar pada titik observasi yang lebih dekat terhadap lokasi target \(s_0\):

\[\hat{z}(s_0) = \frac{\sum_{i \in N_m(s_0)} w_i(s_0) \, z(s_i)}{\sum_{i \in N_m(s_0)} w_i(s_0)}, \quad w_i(s_0) = \frac{1}{\|s_0 - s_i\|^p}\]

Parameter pangkat jarak standar \(p = 2\) (inverse distance squared) dengan lingkungan tetangga lokal \(\text{nmax} = 10\) digunakan agar titik-titik jauh di seberang dataran banjir tidak meredam pola gradien lokal di dekat sungai.

model_idw <- gstat(formula = log_zinc ~ 1, locations = meuse_sp, nmax = 10, set = list(idp = 2))
pred_idw_grid <- predict(model_idw, grid_sp, debug.level = 0)
df_grid$idw <- pred_idw_grid$var1.pred

ggplot(df_grid, aes(x = x, y = y, fill = idw)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "Prediksi\nlog(zinc)") +
  coord_equal() +
  labs(
    title = "Interpolasi Inverse Distance Weighting (IDW, p = 2, nmax = 10)",
    x = "Koordinat X (m)", y = "Koordinat Y (m)"
  ) +
  theme_minimal()

2.4 Metode Ensemble (Kombinasi 3 Metode Interpolasi)

Untuk metode Ensemble, disyaratkan menggabungkan minimal tiga metode interpolasi. Di sini kita menggabungkan tiga metode interpolasi dengan karakteristik mekanisme yang berbeda:

Nearest Neighbour (NN, \(k = 5\)): berbasis rata-rata tetangga lokal tak terboboti.
Inverse Distance Weighting (IDW, \(p = 2, \text{nmax} = 10\)): berbasis pembobotan jarak lokal kontinu.
Thin Plate Spline via Generalized Additive Model (GAM / TPS, mgcv::gam(log_zinc ~ s(x, y, k = 30))): berbasis pemulusan permukaan global (thin-plate regression spline) yang menangkap tren spasial kontinu secara non-parametrik.

Kita membangun dua varian ensemble:

Ensemble Rata-rata Sederhana (Equal-Weight Ensemble): \[\hat{z}_{\text{ens-eq}}(s_0) = \frac{1}{3}\Big(\hat{z}_{\text{NN}}(s_0) + \hat{z}_{\text{IDW}}(s_0) + \hat{z}_{\text{GAM}}(s_0)\Big)\]
Ensemble Terboboti (Inverse-RMSE\(^2\) Weighted Ensemble): \[\hat{z}_{\text{ens-w}}(s_0) = w_1 \hat{z}_{\text{NN}}(s_0) + w_2 \hat{z}_{\text{IDW}}(s_0) + w_3 \hat{z}_{\text{GAM}}(s_0), \quad w_m \propto \frac{1}{\text{RMSE}_m^2}\] di mana bobot \(w_m\) diestimasi dari inner cross-validation pada data latih (sehingga bebas dari data leakage saat evaluasi CV).

# Fit anggota ketiga: Thin Plate Spline (GAM) pada data penuh
fit_gam_full <- gam(log_zinc ~ s(x, y, k = 30), data = meuse)
df_grid$gam_tps <- as.numeric(predict(fit_gam_full, newdata = meuse.grid))

# Estimasi bobot inverse-RMSE^2 menggunakan 5-fold CV pada data penuh untuk peta akhir
set.seed(123)
rmse_fn <- function(obs, pred) sqrt(mean((obs - pred)^2))

compute_inv_weights <- function(df_train, inner_k = 5) {
  n_tr <- nrow(df_train)
  in_folds <- sample(rep(1:inner_k, length.out = n_tr))
  p_nn <- p_idw <- p_gam <- numeric(n_tr)
  for (f in 1:inner_k) {
    itr <- df_train[in_folds != f, ]
    ite <- df_train[in_folds == f, ]
    itr_sp <- itr; coordinates(itr_sp) <- ~x + y
    ite_sp <- ite; coordinates(ite_sp) <- ~x + y

    p_nn[in_folds == f]  <- predict(gstat(formula = log_zinc ~ 1, locations = itr_sp, nmax = 5, set = list(idp = 0)), ite_sp, debug.level = 0)$var1.pred
    p_idw[in_folds == f] <- predict(gstat(formula = log_zinc ~ 1, locations = itr_sp, nmax = 10, set = list(idp = 2)), ite_sp, debug.level = 0)$var1.pred
    p_gam[in_folds == f] <- as.numeric(predict(gam(log_zinc ~ s(x, y, k = 30), data = itr), newdata = ite))
  }
  rmses <- c(
    NN = rmse_fn(df_train$log_zinc, p_nn),
    IDW = rmse_fn(df_train$log_zinc, p_idw),
    GAM_TPS = rmse_fn(df_train$log_zinc, p_gam)
  )
  weights <- (1 / rmses^2) / sum(1 / rmses^2)
  list(weights = weights, rmses = rmses)
}

ens_info_full <- compute_inv_weights(meuse, inner_k = 5)
w_full <- ens_info_full$weights

kable(
  data.frame(
    Anggota_Ensemble = c("Nearest Neighbour (k=5)", "IDW (p=2, nmax=10)", "Thin Plate Spline (GAM)"),
    Inner_CV_RMSE = ens_info_full$rmses,
    Bobot_Ensemble = w_full
  ),
  digits = 4,
  row.names = FALSE,
  caption = "Bobot Optimal Ensemble (Berbasis Inverse-RMSE^2 pada Inner 5-Fold CV)"
)

Bobot Optimal Ensemble (Berbasis Inverse-RMSE^2 pada Inner 5-Fold CV)
Anggota_Ensemble	Inner_CV_RMSE	Bobot_Ensemble
Nearest Neighbour (k=5)	0.4994	0.2732
IDW (p=2, nmax=10)	0.4724	0.3053
Thin Plate Spline (GAM)	0.4020	0.4215

df_grid$ens_eq <- (df_grid$nn5 + df_grid$idw + df_grid$gam_tps) / 3
df_grid$ens_w  <- w_full["NN"] * df_grid$nn5 + w_full["IDW"] * df_grid$idw + w_full["GAM_TPS"] * df_grid$gam_tps

p_ens1 <- ggplot(df_grid, aes(x = x, y = y, fill = ens_eq)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "log(zinc)") +
  coord_equal() +
  labs(title = "Ensemble Rata-Rata (Equal)", x = "X (m)", y = "Y (m)") +
  theme_minimal()

p_ens2 <- ggplot(df_grid, aes(x = x, y = y, fill = ens_w)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "log(zinc)") +
  coord_equal() +
  labs(title = "Ensemble Terboboti (Inv-RMSE^2)", x = "X (m)", y = "Y (m)") +
  theme_minimal()

p_ens1 + p_ens2

2.5 Metode Kriging (Ordinary Kriging)

Ordinary Kriging merupakan estimator linier tak bias terbaik (Best Linear Unbiased Estimator / BLUE) yang memanfaatkan fungsi korelasi spasial empiris melalui semivariogram:

\[\hat{\gamma}(h) = \frac{1}{2|N(h)|} \sum_{(i,j) \in N(h)} \big[z(s_i) - z(s_j)\big]^2\]

Tahapan yang dilakukan:

Menghitung semivariogram empiris untuk log(zinc).
Melakukan fitting model semivariogram teoretis (Spherical, Exponential, atau Matérn) dan memilih model dengan kecocokan terbaik (minimum sum of squared errors).
Melakukan prediksi Ordinary Kriging pada meuse.grid sekaligus memperoleh varians / standar error Kriging sebagai ukuran ketidakpastian spasial.

v_emp <- variogram(log_zinc ~ 1, meuse_sp)
v_fit <- fit.variogram(v_emp, model = vgm(c("Sph", "Exp", "Mat")))

print(v_fit)

##   model      psill    range
## 1   Nug 0.05065971   0.0000
## 2   Sph 0.59060511 897.0011

# Kurva variogram teoretis untuk plot ggplot2
v_line <- variogramLine(v_fit, maxdist = max(v_emp$dist))

ggplot() +
  geom_point(data = v_emp, aes(x = dist, y = gamma, size = np), color = "#2c7fb8", alpha = 0.85) +
  geom_line(data = v_line, aes(x = dist, y = gamma), color = "#d95f02", linewidth = 1.1) +
  scale_size_continuous(name = "Jumlah Pasangan\n(np)", range = c(2, 5)) +
  labs(
    title = paste("Semivariogram Empiris & Model Teoretis Terpilih:", paste(v_fit$model, collapse = " + ")),
    subtitle = sprintf("Nugget = %.4f | Partial Sill = %.4f | Range = %.1f m",
                       v_fit$psill[1], v_fit$psill[2], v_fit$range[2]),
    x = "Jarak h (meter)",
    y = "Semivariansi gamma(h)"
  ) +
  theme_minimal()

pred_ok_grid <- krige(log_zinc ~ 1, meuse_sp, grid_sp, model = v_fit, debug.level = 0)
df_grid$ok_pred <- pred_ok_grid$var1.pred
df_grid$ok_se   <- sqrt(pred_ok_grid$var1.var)

p_ok1 <- ggplot(df_grid, aes(x = x, y = y, fill = ok_pred)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "Prediksi\nlog(zinc)") +
  coord_equal() +
  labs(title = "Prediksi Ordinary Kriging", x = "X (m)", y = "Y (m)") +
  theme_minimal()

p_ok2 <- ggplot(df_grid, aes(x = x, y = y, fill = ok_se)) +
  geom_raster() +
  geom_point(data = meuse, aes(x = x, y = y), inherit.aes = FALSE, size = 0.7, color = "white", alpha = 0.6) +
  scale_fill_viridis_c(option = "magma", name = "Std. Error\nlog(zinc)") +
  coord_equal() +
  labs(title = "Standar Error Ordinary Kriging", x = "X (m)", y = "Y (m)") +
  theme_minimal()

p_ok1 + p_ok2

Terlihat bahwa standar error Kriging paling kecil di sekitar titik-titik sampel dan meningkat di daerah yang jauh dari observasi.

3 Evaluasi & Perbandingan RMSE dengan Cross-Validation

Untuk membandingkan akurasi seluruh metode secara adil dan bebas kebocoran data (zero data leakage), digunakan partisi cross-validation yang identik untuk semua metode, dengan dua skema validasi:

10-Fold Cross-Validation (set.seed(123)): data dibagi menjadi 10 partisi; pada setiap iterasi, 9 fold digunakan untuk melatih model (termasuk refitting model GAM, estimasi bobot inner-CV ensemble terboboti, dan refitting semivariogram Kriging) lalu diprediksi pada 1 fold uji.
Leave-One-Out Cross-Validation (LOOCV, \(k = n = 155\)): skema deterministik tanpa variasi pengacakan fold, di mana setiap 1 titik observasi diprediksi dari \(154\) titik lainnya.

Metrik utama perbandingan adalah Root Mean Squared Error (RMSE) pada skala log(zinc):

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

Dilengkapi pula metrik MAE, \(R^2\), serta RMSE pada skala asli zinc (ppm) setelah transformasi balik \(\hat{z}_i = \exp(\hat{y}_i)\).

3.1 10-Fold Cross-Validation

set.seed(123)
n <- nrow(meuse)
k_folds <- 10
folds <- sample(rep(1:k_folds, length.out = n))

mae_fn <- function(obs, pred) mean(abs(obs - pred))
r2_fn  <- function(obs, pred) 1 - sum((obs - pred)^2) / sum((obs - mean(obs))^2)

cv10_preds <- data.frame(
  obs = meuse$log_zinc,
  obs_ppm = meuse$zinc,
  closest = numeric(n),
  nn5 = numeric(n),
  idw = numeric(n),
  gam_tps = numeric(n),
  ens_eq = numeric(n),
  ens_w = numeric(n),
  ok = numeric(n)
)

for (f in 1:k_folds) {
  tr_idx <- which(folds != f)
  te_idx <- which(folds == f)
  tr_df <- meuse[tr_idx, ]
  te_df <- meuse[te_idx, ]
  tr_sp <- tr_df; coordinates(tr_sp) <- ~x + y
  te_sp <- te_df; coordinates(te_sp) <- ~x + y

  # 1. Closest Observation (1-NN)
  p_c <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 1), te_sp, debug.level = 0)$var1.pred

  # 2. Nearest Neighbour (k = 5)
  p_nn <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 5, set = list(idp = 0)), te_sp, debug.level = 0)$var1.pred

  # 3. IDW (p = 2, nmax = 10)
  p_idw <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 10, set = list(idp = 2)), te_sp, debug.level = 0)$var1.pred

  # 4. Thin Plate Spline / GAM (anggota ke-3 ensemble)
  fit_gam_f <- gam(log_zinc ~ s(x, y, k = 30), data = tr_df)
  p_gam <- as.numeric(predict(fit_gam_f, newdata = te_df))

  # 5. Ensemble Equal & Weighted (bobot diestimasi ketat di dalam training fold)
  p_ens_eq <- (p_nn + p_idw + p_gam) / 3
  w_f <- compute_inv_weights(tr_df, inner_k = 5)$weights
  p_ens_w <- w_f["NN"] * p_nn + w_f["IDW"] * p_idw + w_f["GAM_TPS"] * p_gam

  # 6. Ordinary Kriging (refit semivariogram di setiap training fold)
  v_f <- variogram(log_zinc ~ 1, tr_sp)
  m_f <- fit.variogram(v_f, vgm(c("Sph", "Exp", "Mat")))
  p_ok <- krige(log_zinc ~ 1, tr_sp, te_sp, model = m_f, debug.level = 0)$var1.pred

  cv10_preds$closest[te_idx] <- p_c
  cv10_preds$nn5[te_idx]     <- p_nn
  cv10_preds$idw[te_idx]     <- p_idw
  cv10_preds$gam_tps[te_idx] <- p_gam
  cv10_preds$ens_eq[te_idx]  <- p_ens_eq
  cv10_preds$ens_w[te_idx]   <- p_ens_w
  cv10_preds$ok[te_idx]      <- p_ok
}

method_labels <- c(
  closest = "Closest Observation (1-NN)",
  nn5     = "Nearest Neighbour (k = 5)",
  idw     = "IDW (p = 2, nmax = 10)",
  gam_tps = "Thin Plate Spline / GAM (Anggota Ensemble)",
  ens_eq  = "Ensemble Rata-Rata (NN + IDW + GAM)",
  ens_w   = "Ensemble Terboboti (NN + IDW + GAM)",
  ok      = "Ordinary Kriging"
)

eval_10fold <- lapply(names(method_labels), function(col) {
  pred_log <- cv10_preds[[col]]
  pred_ppm <- exp(pred_log)
  data.frame(
    Metode = method_labels[[col]],
    RMSE_log = rmse_fn(cv10_preds$obs, pred_log),
    MAE_log  = mae_fn(cv10_preds$obs, pred_log),
    R2_log   = r2_fn(cv10_preds$obs, pred_log),
    RMSE_ppm = rmse_fn(cv10_preds$obs_ppm, pred_ppm)
  )
}) %>%
  bind_rows() %>%
  arrange(RMSE_log)

kable(eval_10fold, digits = 4, caption = "Hasil Perbandingan 10-Fold Cross-Validation (Diurutkan berdasarkan RMSE Terkecil)")

Hasil Perbandingan 10-Fold Cross-Validation (Diurutkan berdasarkan RMSE Terkecil)
Metode	RMSE_log	MAE_log	R2_log	RMSE_ppm
Thin Plate Spline / GAM (Anggota Ensemble)	0.4067	0.3057	0.6805	233.1199
Ordinary Kriging	0.4106	0.3105	0.6744	234.9977
Ensemble Terboboti (NN + IDW + GAM)	0.4327	0.3262	0.6383	248.6839
Ensemble Rata-Rata (NN + IDW + GAM)	0.4405	0.3320	0.6252	253.4406
IDW (p = 2, nmax = 10)	0.4660	0.3529	0.5806	269.1210
Nearest Neighbour (k = 5)	0.5057	0.3855	0.5060	284.8565
Closest Observation (1-NN)	0.5817	0.4101	0.3464	329.8958

3.2 Leave-One-Out Cross-Validation (LOOCV)

lo_preds <- data.frame(
  obs = meuse$log_zinc,
  obs_ppm = meuse$zinc,
  closest = numeric(n),
  nn5 = numeric(n),
  idw = numeric(n),
  gam_tps = numeric(n),
  ens_eq = numeric(n),
  ok = numeric(n)
)

for (i in 1:n) {
  tr_df <- meuse[-i, ]
  te_df <- meuse[i, , drop = FALSE]
  tr_sp <- tr_df; coordinates(tr_sp) <- ~x + y
  te_sp <- te_df; coordinates(te_sp) <- ~x + y

  lo_preds$closest[i] <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 1), te_sp, debug.level = 0)$var1.pred
  lo_preds$nn5[i]     <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 5, set = list(idp = 0)), te_sp, debug.level = 0)$var1.pred
  lo_preds$idw[i]     <- predict(gstat(formula = log_zinc ~ 1, locations = tr_sp, nmax = 10, set = list(idp = 2)), te_sp, debug.level = 0)$var1.pred
  lo_preds$gam_tps[i] <- as.numeric(predict(gam(log_zinc ~ s(x, y, k = 30), data = tr_df), newdata = te_df))
  lo_preds$ens_eq[i]  <- (lo_preds$nn5[i] + lo_preds$idw[i] + lo_preds$gam_tps[i]) / 3
  lo_preds$ok[i]      <- krige(log_zinc ~ 1, tr_sp, te_sp, model = v_fit, debug.level = 0)$var1.pred
}

lo_labels <- method_labels[c("closest", "nn5", "idw", "gam_tps", "ens_eq", "ok")]

eval_loocv <- lapply(names(lo_labels), function(col) {
  pred_log <- lo_preds[[col]]
  pred_ppm <- exp(pred_log)
  data.frame(
    Metode = lo_labels[[col]],
    RMSE_LOOCV_log = rmse_fn(lo_preds$obs, pred_log),
    MAE_LOOCV_log  = mae_fn(lo_preds$obs, pred_log),
    R2_LOOCV_log   = r2_fn(lo_preds$obs, pred_log),
    RMSE_LOOCV_ppm = rmse_fn(lo_preds$obs_ppm, pred_ppm)
  )
}) %>%
  bind_rows() %>%
  arrange(RMSE_LOOCV_log)

kable(eval_loocv, digits = 4, caption = "Hasil Perbandingan Leave-One-Out Cross-Validation (LOOCV)")

Hasil Perbandingan Leave-One-Out Cross-Validation (LOOCV)
Metode	RMSE_LOOCV_log	MAE_LOOCV_log	R2_LOOCV_log	RMSE_LOOCV_ppm
Ordinary Kriging	0.3918	0.2921	0.7035	224.8700
Thin Plate Spline / GAM (Anggota Ensemble)	0.3989	0.2966	0.6927	227.7349
Ensemble Rata-Rata (NN + IDW + GAM)	0.4252	0.3223	0.6509	249.2406
IDW (p = 2, nmax = 10)	0.4501	0.3414	0.6087	264.7435
Nearest Neighbour (k = 5)	0.4738	0.3624	0.5664	275.9611
Closest Observation (1-NN)	0.5655	0.3964	0.3824	309.9685

3.3 Visualisasi Perbandingan RMSE & Diagnostik Error

comp_plot_df <- eval_10fold %>%
  dplyr::select(Metode, `10-Fold CV` = RMSE_log) %>%
  left_join(eval_loocv %>% dplyr::select(Metode, `LOOCV` = RMSE_LOOCV_log), by = "Metode") %>%
  pivot_longer(cols = c(`10-Fold CV`, `LOOCV`), names_to = "Skema_CV", values_to = "RMSE") %>%
  filter(!is.na(RMSE))

# Urutkan metode berdasarkan RMSE 10-Fold CV
ordered_methods <- eval_10fold$Metode
comp_plot_df$Metode <- factor(comp_plot_df$Metode, levels = rev(ordered_methods))

ggplot(comp_plot_df, aes(x = Metode, y = RMSE, fill = Skema_CV)) +
  geom_col(position = position_dodge(width = 0.72), width = 0.65) +
  geom_text(
    aes(label = sprintf("%.4f", RMSE)),
    position = position_dodge(width = 0.72),
    hjust = -0.1, size = 3.2
  ) +
  coord_flip(ylim = c(0, max(comp_plot_df$RMSE, na.rm = TRUE) * 1.15)) +
  scale_fill_manual(values = c("10-Fold CV" = "#2c7fb8", "LOOCV" = "#41ab5d"), name = "Skema CV") +
  labs(
    title = "Perbandingan RMSE Cross-Validation Antar Metode Interpolasi",
    subtitle = "Evaluasi pada target log(zinc) — nilai lebih rendah menunjukkan akurasi lebih baik",
    x = NULL, y = "RMSE log(zinc)"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

scatter_df <- cv10_preds %>%
  dplyr::select(obs, closest, nn5, idw, ens_w, ok) %>%
  pivot_longer(cols = -obs, names_to = "kode", values_to = "pred") %>%
  mutate(Metode = factor(kode,
    levels = c("closest", "nn5", "idw", "ens_w", "ok"),
    labels = c("Closest Observation (1-NN)", "Nearest Neighbour (k=5)",
               "IDW (p=2, nmax=10)", "Ensemble Terboboti", "Ordinary Kriging")
  ))

ggplot(scatter_df, aes(x = obs, y = pred)) +
  geom_point(alpha = 0.65, color = "#2c7fb8", size = 1.6) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#d95f02", linewidth = 0.8) +
  facet_wrap(~Metode, ncol = 3) +
  coord_equal() +
  labs(
    title = "Nilai Aktual vs Prediksi (10-Fold Cross-Validation)",
    x = "Observasi log(zinc)", y = "Prediksi log(zinc)"
  ) +
  theme_minimal()

4 Perbandingan Peta Interpolasi Keseluruhan

Sebagai rangkuman visual, berikut peta prediksi seluruh metode interpolasi utama pada domain meuse.grid dengan skala warna yang diseragamkan agar perbedaan karakter kehalusan (smoothness) dan detail lokal dapat diamati langsung.

# Pastikan kolom prediksi Ordinary Kriging tersedia meskipun chunk dijalankan tidak berurutan di RStudio
if (!"ok_pred" %in% names(df_grid)) {
  if (!exists("v_fit")) {
    v_emp <- variogram(log_zinc ~ 1, meuse_sp)
    v_fit <- fit.variogram(v_emp, model = vgm(c("Sph", "Exp", "Mat")))
  }
  pred_ok_grid <- krige(log_zinc ~ 1, meuse_sp, grid_sp, model = v_fit, debug.level = 0)
  df_grid$ok_pred <- pred_ok_grid$var1.pred
  df_grid$ok_se   <- sqrt(pred_ok_grid$var1.var)
}

maps_long <- df_grid %>%
  dplyr::select(x, y, closest, nn5, idw, ens_w, ok_pred) %>%
  pivot_longer(cols = c(closest, nn5, idw, ens_w, ok_pred),
               names_to = "Metode", values_to = "pred") %>%
  mutate(Metode = factor(Metode,
    levels = c("closest", "nn5", "idw", "ens_w", "ok_pred"),
    labels = c(
      "1. Closest Observation (1-NN)",
      "2. Nearest Neighbour (k=5)",
      "3. IDW (p=2, nmax=10)",
      "4. Ensemble Terboboti (NN+IDW+GAM)",
      "5. Ordinary Kriging"
    )
  ))

ggplot(maps_long, aes(x = x, y = y, fill = pred)) +
  geom_raster() +
  scale_fill_viridis_c(option = "viridis", name = "Prediksi\nlog(zinc)") +
  facet_wrap(~Metode, ncol = 3) +
  coord_equal() +
  labs(
    title = "Perbandingan Peta Hasil Interpolasi Spasial log(zinc) pada Meuse Grid",
    x = "Koordinat X (m)", y = "Koordinat Y (m)"
  ) +
  theme_minimal() +
  theme(axis.text = element_text(size = 7))

5 Kesimpulan

Berdasarkan hasil interpolasi spasial log(zinc) pada dataset meuse dan evaluasi cross-validation (10-fold CV maupun LOOCV):

Closest Observation (1-NN) menghasilkan error terbesar (RMSE 10-fold CV = 0.5817; LOOCV = 0.5655) karena hanya mengandalkan satu titik terdekat sehingga sangat sensitif terhadap noise lokal dan menghasilkan permukaan diskontinu berbentuk poligon Voronoi.
Nearest Neighbour (\(k = 5\)) menurunkan error secara signifikan (RMSE 10-fold CV = 0.5057) karena merata-ratakan variasi acak lokal, namun masih menunjukkan pola transisi bertahap (step-like).
IDW (\(p = 2, \text{nmax} = 10\)) menghasilkan permukaan lebih halus dan akurasi lebih baik daripada NN dan Closest Observation (RMSE 10-fold CV = 0.4660; LOOCV = 0.4501) karena memperhitungkan jarak secara kontinu.
Ensemble (NN + IDW + Thin Plate Spline/GAM) mampu menggabungkan kekuatan interpolator lokal (NN, IDW) dengan pemulusan permukaan global (GAM / TPS), menurunkan RMSE dibandingkan metode lokal individual (RMSE 10-fold CV Ensemble Terboboti = 0.4327).
Ordinary Kriging memberikan performa yang sangat unggul dan stabil (RMSE 10-fold CV = 0.4106; LOOCV = 0.3918) — terbaik pada LOOCV dan bersaing ketat dengan anggota spline GAM pada 10-fold CV — sekaligus menyediakan ukuran ketidakpastian spasial (standar error Kriging) melalui pemodelan semivariogram.

Perbandingan Metode Interpolasi Spasial pada Data Meuse

Closest Observation, Nearest Neighbour (NN), IDW, Ensemble, dan Ordinary Kriging dengan Evaluasi Cross-Validation

Analisis Geostatistika & Interpolasi Spasial

10 November 2026

1 Pendahuluan & Eksplorasi Data `meuse`

2 Implementasi Metode Interpolasi Spasial

2.1 Closest Observation (1-Nearest Neighbor / Voronoi)

2.2 Nearest Neighbour (NN / \(k\)-Nearest Neighbors)

2.3 Inverse Distance Weighting (IDW)

2.4 Metode Ensemble (Kombinasi 3 Metode Interpolasi)

2.5 Metode Kriging (Ordinary Kriging)

3 Evaluasi & Perbandingan RMSE dengan Cross-Validation

3.1 10-Fold Cross-Validation

3.2 Leave-One-Out Cross-Validation (LOOCV)

3.3 Visualisasi Perbandingan RMSE & Diagnostik Error

4 Perbandingan Peta Interpolasi Keseluruhan

5 Kesimpulan

Perbandingan Metode Interpolasi Spasial pada Data Meuse

Closest Observation, Nearest Neighbour (NN), IDW, Ensemble, dan Ordinary Kriging dengan Evaluasi Cross-Validation

Analisis Geostatistika & Interpolasi Spasial

10 November 2026

1 Pendahuluan & Eksplorasi Data meuse

2 Implementasi Metode Interpolasi Spasial

2.1 Closest Observation (1-Nearest Neighbor / Voronoi)

2.2 Nearest Neighbour (NN / \(k\)-Nearest Neighbors)

2.3 Inverse Distance Weighting (IDW)

2.4 Metode Ensemble (Kombinasi 3 Metode Interpolasi)

2.5 Metode Kriging (Ordinary Kriging)

3 Evaluasi & Perbandingan RMSE dengan Cross-Validation

3.1 10-Fold Cross-Validation

3.2 Leave-One-Out Cross-Validation (LOOCV)

3.3 Visualisasi Perbandingan RMSE & Diagnostik Error

4 Perbandingan Peta Interpolasi Keseluruhan

5 Kesimpulan

1 Pendahuluan & Eksplorasi Data `meuse`