Putri Dwi Suci Sri Mulyaningsih Sitanggang (M0401241033)

Library

# install.packages(c("sf", "sp", "gstat", "ggplot2", "knitr"))
library(sf)
library(sp)
library(gstat)
library(ggplot2)

Data

Variabel yang diinterpolasi: kandungan tembaga (copper) pada data meuse, ditransformasi log.

data(meuse); data(meuse.grid); data(meuse.riv)

meuse$log_cu <- log(meuse$copper)
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, remove = FALSE)
grid_sf  <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)
sungai   <- setNames(as.data.frame(meuse.riv), c("x", "y"))

n_obs <- nrow(meuse_sf)
summary(meuse$copper)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.00   23.00   31.00   40.32   49.50  128.00
batas_x <- range(meuse.grid$x) + c(-150, 150)
batas_y <- range(meuse.grid$y) + c(-150, 150)

ggplot() +
  geom_polygon(data = sungai, aes(x, y), fill = "#bcd8f2", colour = "#4a90c2", alpha = 0.7) +
  geom_point(data = meuse, aes(x, y, colour = copper, size = copper), alpha = 0.85) +
  scale_colour_viridis_c(option = "C", name = "Cu (ppm)") +
  scale_size_continuous(range = c(1, 5), guide = "none") +
  coord_equal(xlim = batas_x, ylim = batas_y) +
  labs(title = "Lokasi 155 titik sampel dan kandungan tembaga", x = "x (m)", y = "y (m)") +
  theme_minimal(base_size = 12)

Kandungan Cu menjulur ke kanan (kemencengan 1.39, turun menjadi 0.65 setelah di-log), sehingga seluruh analisis memakai log(Cu). Titik bernilai tinggi (kuning-oranye) berderet di sepanjang tepi sungai, sedangkan titik bernilai rendah (ungu) berada di bagian dalam dataran. Titik yang berdekatan bernilai mirip, jadi ada autokorelasi spasial dan interpolasi layak dilakukan.

Grid Prediksi

Lokasi prediksi adalah 3.103 sel grid meuse.grid.

ggplot() +
  geom_polygon(data = sungai, aes(x, y), fill = "#bcd8f2", colour = "#4a90c2", alpha = 0.7) +
  geom_tile(data = meuse.grid, aes(x, y), fill = "grey80") +
  geom_point(data = meuse, aes(x, y), size = 0.8) +
  coord_equal(xlim = batas_x, ylim = batas_y) +
  labs(title = "Grid prediksi dan titik sampel", x = "x (m)", y = "y (m)") +
  theme_minimal(base_size = 12)

RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))
batas_warna <- range(meuse$copper)

# Peta prediksi pada grid (input skala log, ditampilkan dalam ppm)
peta_grid <- function(pred_log, judul, opsi = "C", legenda = "Cu (ppm)", ke_ppm = TRUE, batas = batas_warna) {
  df <- data.frame(x = grid_sf$x, y = grid_sf$y, nilai = if (ke_ppm) exp(pred_log) else pred_log)
  p <- ggplot(df, aes(x, y, fill = nilai)) +
    geom_raster() +
    geom_polygon(data = sungai, aes(x, y), inherit.aes = FALSE, fill = NA, colour = "#4a90c2") +
    geom_point(data = meuse, aes(x, y), inherit.aes = FALSE, size = 0.5, colour = "grey15") +
    coord_equal(xlim = batas_x, ylim = batas_y) +
    labs(title = judul, x = "x (m)", y = "y (m)") +
    theme_minimal(base_size = 12)
  if (is.null(batas)) p + scale_fill_viridis_c(option = opsi, name = legenda)
  else p + scale_fill_viridis_c(option = opsi, name = legenda, limits = batas, oob = scales::squish)
}

Closest Observation (Voronoi)

vor <- st_collection_extract(st_voronoi(st_union(st_geometry(meuse_sf))), "POLYGON")
vor_sf <- st_sf(geometry = vor, crs = st_crs(meuse_sf))
idx <- st_intersects(vor_sf, meuse_sf)
vor_sf$cu <- meuse_sf$copper[sapply(idx, function(i) i[1])]

batas_studi <- st_convex_hull(st_union(st_geometry(grid_sf)))
st_agr(vor_sf) <- "constant"
vor_potong <- suppressWarnings(st_intersection(vor_sf, batas_studi))

ggplot() +
  geom_sf(data = vor_potong, aes(fill = cu), colour = "white", linewidth = 0.2) +
  geom_sf(data = meuse_sf, size = 0.6, colour = "grey15") +
  scale_fill_viridis_c(option = "C", name = "Cu (ppm)") +
  coord_sf(datum = st_crs(28992)) +
  labs(title = "Diagram Voronoi", x = "x (m)", y = "y (m)") +
  theme_minimal(base_size = 12)

# Nilai titik sampel terdekat (ekuivalen dengan poligon Voronoi)
p_closest <- function(train, test) train$log_cu[st_nearest_feature(test, train)]

pred_closest <- p_closest(meuse_sf, grid_sf)
peta_grid(pred_closest, "Closest Observation (Voronoi)")

Peta berblok dengan batas yang tegas dan hanya berisi 60 nilai berbeda, persis nilai asli titik sampel tanpa penghalusan. Satu titik bernilai ekstrem menguasai seluruh poligonnya, sehingga permukaannya diskontinu.

Nearest Neighbors (NN)

Nilai \(k\) dipilih dengan Leave-One-Out CV.

p_nn <- function(train, test, k) {
  idw(log_cu ~ 1, locations = train, newdata = test, idp = 0, nmax = k, debug.level = 0)$var1.pred
}
loo_rmse <- function(f) {
  pred <- sapply(seq_len(n_obs), function(i) f(meuse_sf[-i, ], meuse_sf[i, ]))
  RMSE(meuse_sf$log_cu, pred)
}

kandidat_k <- c(1, 2, 3, 4, 5, 7, 10, 15, 20)
rmse_k <- sapply(kandidat_k, function(k) loo_rmse(function(tr, te) p_nn(tr, te, k)))
k_opt  <- kandidat_k[which.min(rmse_k)]
knitr::kable(data.frame(k = kandidat_k, RMSE_LOOCV = rmse_k), digits = 4)
k RMSE_LOOCV
1 0.4441
2 0.3949
3 0.3387
4 0.3611
5 0.3714
7 0.3750
10 0.3838
15 0.4074
20 0.4238 
pred_nn <- p_nn(meuse_sf, grid_sf, k_opt)
peta_grid(pred_nn, paste0("Nearest Neighbors (k = ", k_opt, ")"))

\(k\) = 3 memberi RMSE LOOCV terkecil (0.3387), lebih baik dibanding \(k\) = 1 (0.4441) dan memburuk lagi saat \(k\) membesar karena tetangga yang jauh ikut dirata-ratakan. Peta lebih halus daripada Closest Observation; rentang prediksi menyempit menjadi 15.6 sampai 110.5 ppm (data asli 14 sampai 128 ppm) karena nilai ekstrem diredam oleh rata-rata.

Inverse Distance Weighting (IDW)

Nilai power \(p\) dipilih dengan Leave-One-Out CV.

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

kandidat_p <- c(0.5, 1, 1.5, 2, 3, 4, 5)
rmse_p <- sapply(kandidat_p, function(p) loo_rmse(function(tr, te) p_idw(tr, te, p)))
p_opt  <- kandidat_p[which.min(rmse_p)]
knitr::kable(data.frame(p = kandidat_p, RMSE_LOOCV = rmse_p), digits = 4)
p RMSE_LOOCV
0.5 0.4912
1.0 0.4597
1.5 0.4218
2.0 0.3905
3.0 0.3665
4.0 0.3720
5.0 0.3830 
pred_idw <- p_idw(meuse_sf, grid_sf, p_opt)
peta_grid(pred_idw, paste0("IDW (p = ", p_opt, ")"))

\(p\) = 3 memberi RMSE LOOCV terkecil (0.3665); pada \(p\) kecil kesalahan jauh lebih besar (0.4912 pada \(p\) = 0.5) karena titik jauh terlalu berpengaruh. Peta mulus dan kontinu, tetapi muncul pola lingkaran kecil di sekitar titik sampel bernilai tinggi, dan IDW tidak memberi ukuran ketidakpastian prediksi.

Ordinary Kriging

Variogram

v_cloud <- variogram(log_cu ~ 1, meuse_sf, cloud = TRUE)
ggplot(v_cloud, aes(dist, gamma)) +
  geom_point(alpha = 0.12, size = 0.8, colour = "#7570b3") +
  labs(title = "Variogram cloud log(Cu)", x = "Jarak (m)", y = "Semivariansi") +
  theme_minimal(base_size = 12)

v_emp <- variogram(log_cu ~ 1, meuse_sf)
ggplot(v_emp, aes(dist, gamma)) +
  geom_point(aes(size = np), colour = "#1b9e77") +
  scale_size_continuous(name = "Jumlah pasangan", range = c(2, 6)) +
  expand_limits(y = 0) +
  labs(title = "Variogram empiris log(Cu)", x = "Jarak (m)", y = "Semivariansi") +
  theme_minimal(base_size = 12)

Pada variogram cloud, pasangan titik berjarak dekat memiliki semivariansi yang lebih rendah dan sebarannya melebar pada jarak jauh. Variogram empiris naik dari 0.095 pada jarak terdekat lalu mendatar di kisaran 0.26 sampai 0.31. Pola ini menunjukkan adanya ketergantungan spasial: makin jauh dua titik, makin berbeda nilainya.

Model Semivariogram

Tiga model dicocokkan (Spherical, Exponential, Gaussian), dipilih yang SSErr-nya terkecil.

nilai_awal <- function(m, v) {
  vgm(psill = max(v$gamma) * 0.8, model = m, range = max(v$dist) / 3, nugget = min(v$gamma))
}
kandidat_model <- c("Sph", "Exp", "Gau")
fit_model <- setNames(lapply(kandidat_model, function(m) fit.variogram(v_emp, nilai_awal(m, v_emp))),
                      kandidat_model)

tabel_fit <- do.call(rbind, lapply(kandidat_model, function(m) {
  f <- fit_model[[m]]
  data.frame(Model = m, Nugget = f$psill[1], Partial_sill = f$psill[2],
             Sill_total = sum(f$psill), Range_m = f$range[2], SSErr = attr(f, "SSErr"))
}))
knitr::kable(tabel_fit, digits = 4)
Model Nugget Partial_sill Sill_total Range_m SSErr
Sph 0.0610 0.2275 0.2885 775.0711 0
Exp 0.0373 0.2748 0.3122 366.5063 0
Gau 0.0891 0.1933 0.2825 348.9513 0 
model_terpilih <- tabel_fit$Model[which.min(tabel_fit$SSErr)]
fit_terpilih   <- fit_model[[model_terpilih]]
par_t <- tabel_fit[tabel_fit$Model == model_terpilih, ]

garis_model <- do.call(rbind, lapply(kandidat_model, function(m) {
  g <- variogramLine(fit_model[[m]], maxdist = max(v_emp$dist)); g$Model <- m; g
}))
ggplot() +
  geom_point(data = v_emp, aes(dist, gamma), size = 3, colour = "grey25") +
  geom_line(data = garis_model, aes(dist, gamma, colour = Model), linewidth = 1) +
  scale_colour_brewer(palette = "Dark2") +
  expand_limits(y = 0) +
  labs(title = "Variogram empiris dan model terfit", x = "Jarak (m)", y = "Semivariansi") +
  theme_minimal(base_size = 12)

Model terpilih adalah Sph (SSErr terkecil = 2.27^{-6}) dengan nugget 0.061, partial sill 0.227 (sill total 0.288), dan range 775 m. Nugget hanya 21% dari sill total, artinya 79% variasi terstruktur secara spasial, dan titik yang terpisah lebih dari 775 m praktis tidak lagi berkorelasi.

Prediksi dan Variansi

# Variogram dicocokkan ulang pada data latih yang diberikan (dipakai juga saat CV)
p_krig <- function(train, test, model = model_terpilih) {
  v   <- variogram(log_cu ~ 1, train)
  fit <- suppressWarnings(try(fit.variogram(v, nilai_awal(model, v)), silent = TRUE))
  if (inherits(fit, "try-error") || isTRUE(attr(fit, "singular"))) {
    fit <- vgm(psill = var(train$log_cu) * 0.7, model = model,
               range = max(v$dist) / 3, nugget = var(train$log_cu) * 0.3)
  }
  krige(log_cu ~ 1, locations = train, newdata = test, model = fit, debug.level = 0)$var1.pred
}

krig_grid <- krige(log_cu ~ 1, locations = meuse_sf, newdata = grid_sf,
                   model = fit_terpilih, debug.level = 0)
pred_krig <- krig_grid$var1.pred

peta_grid(pred_krig, "Ordinary Kriging")

peta_grid(krig_grid$var1.var, "Variansi Ordinary Kriging (skala log)",
          opsi = "A", legenda = "Variansi", ke_ppm = FALSE, batas = NULL)

Permukaan kriging paling halus dan natural: tidak ada blok seperti Closest Observation maupun lingkaran kecil seperti IDW. Cu tinggi terkonsentrasi di dekat sungai dan rendah di bagian dalam dataran. Variansi prediksi berkisar 0.086 sampai 0.263; rendah di sekitar titik sampel dan tinggi di tepi area yang jarang sampel (korelasi dengan jarak ke titik sampel terdekat = 0.92), sehingga prediksi di area padat sampel lebih dapat dipercaya.

Ensemble

Dua versi ensemble:

  • Ensemble (Equal): NN + IDW + Kriging, bobot sama (1/3).
  • Ensemble (Weighted): Closest + NN + IDW + Kriging, bobot proporsional terhadap 1/MSE. Bobot dihitung dari CV internal pada data latih saja, bukan dari data uji.
semua_metode <- function(train, test) {
  cbind(Closest = p_closest(train, test),
        NN      = p_nn(train, test, k_opt),
        IDW     = p_idw(train, test, p_opt),
        Kriging = p_krig(train, test))
}
buat_fold <- function(n, K) sample(rep(1:K, length.out = n))

hitung_bobot <- function(train, K_dalam = 5) {
  fold_dalam <- buat_fold(nrow(train), K_dalam)
  P <- matrix(NA, nrow(train), 4)
  for (j in 1:K_dalam) {
    P[fold_dalam == j, ] <- semua_metode(train[fold_dalam != j, ], train[fold_dalam == j, ])
  }
  w <- 1 / apply(P, 2, function(p) mean((train$log_cu - p)^2))
  w / sum(w)
}

set.seed(2025)
bobot_penuh <- hitung_bobot(meuse_sf)
names(bobot_penuh) <- c("Closest", "NN", "IDW", "Kriging")
knitr::kable(data.frame(Metode = names(bobot_penuh), Bobot = bobot_penuh),
             digits = 3, row.names = FALSE)
Metode Bobot
Closest 0.188
NN 0.251
IDW 0.253
Kriging 0.308 
P_grid <- cbind(Closest = pred_closest, NN = pred_nn, IDW = pred_idw, Kriging = pred_krig)
pred_ens_eq <- rowMeans(P_grid[, c("NN", "IDW", "Kriging")])
pred_ens_w  <- as.numeric(P_grid %*% bobot_penuh)

peta_grid(pred_ens_eq, "Ensemble (Equal): NN + IDW + Kriging")

peta_grid(pred_ens_w,  "Ensemble (Weighted): Closest + NN + IDW + Kriging")

Kedua peta ensemble halus dan stabil, mewarisi karakter IDW, NN, dan Kriging. Bobot terbesar jatuh pada Kriging (0.308) dan terkecil pada Closest (0.188), artinya metode yang lebih akurat di validasi internal diberi kontribusi lebih besar. Selisih bobotnya tidak begitu ekstrem, sehingga kedua versi ensemble masih menghasilkan peta yang mirip.

Perbandingan Peta

nama_metode <- c("Closest Observation", paste0("Nearest Neighbors (k=", k_opt, ")"),
                 paste0("IDW (p=", p_opt, ")"), "Ordinary Kriging",
                 "Ensemble (Equal)", "Ensemble (Weighted)")
daftar_pred <- list(pred_closest, pred_nn, pred_idw, pred_krig, pred_ens_eq, pred_ens_w)

semua_peta <- do.call(rbind, lapply(seq_along(nama_metode), function(i) {
  data.frame(x = grid_sf$x, y = grid_sf$y, cu = exp(daftar_pred[[i]]),
             Metode = factor(nama_metode[i], levels = nama_metode))
}))

ggplot(semua_peta, aes(x, y, fill = cu)) +
  geom_raster() +
  facet_wrap(~Metode, ncol = 2) +
  scale_fill_viridis_c(option = "C", name = "Cu (ppm)", limits = batas_warna, oob = scales::squish) +
  coord_equal() +
  labs(title = "Perbandingan peta prediksi Cu", x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(axis.text = element_blank(), panel.grid = element_blank())

Semua metode menangkap pola besar yang sama (Cu tinggi di tepi sungai, rendah di bagian dalam), tetapi berbeda tekstur: Closest Observation berblok, NN berbercak halus, IDW mulus dengan lingkaran kecil di sekitar titik, sedangkan Kriging dan kedua Ensemble paling mulus.

Cross Validation

K-fold Cross Validation dengan \(K\) = 5. Semua metode dievaluasi pada fold yang sama yaitu variogram dan bobot ensemble dihitung ulang hanya dari data latih di setiap fold.

K <- 5
set.seed(2025)
fold_id <- buat_fold(n_obs, K)

nama_cv <- c("Closest", "NN", "IDW", "Kriging", "Ens_Equal", "Ens_Weighted")
rmse_log <- matrix(NA, K, 6, dimnames = list(paste("Fold", 1:K), nama_cv))
rmse_ppm <- rmse_log
bobot_fold <- matrix(NA, K, 4, dimnames = list(paste("Fold", 1:K), c("Closest", "NN", "IDW", "Kriging")))

for (i in 1:K) {
  latih <- meuse_sf[fold_id != i, ]
  uji   <- meuse_sf[fold_id == i, ]

  P <- semua_metode(latih, uji)
  w <- hitung_bobot(latih)
  bobot_fold[i, ] <- w
  P_cv <- cbind(P, Ens_Equal = rowMeans(P[, c("NN", "IDW", "Kriging")]),
                   Ens_Weighted = as.numeric(P %*% w))

  for (m in nama_cv) {
    rmse_log[i, m] <- RMSE(uji$log_cu, P_cv[, m])
    rmse_ppm[i, m] <- RMSE(uji$copper, exp(P_cv[, m]))
  }
}
knitr::kable(data.frame(Fold = 1:K, rmse_log), digits = 4, row.names = FALSE,
             caption = "RMSE tiap fold (skala log)")
RMSE tiap fold (skala log)
Fold Closest NN IDW Kriging Ens_Equal Ens_Weighted
1 0.4154 0.3481 0.3492 0.3056 0.3204 0.3226
2 0.3428 0.3909 0.3588 0.3440 0.3585 0.3423
3 0.4493 0.3215 0.3685 0.3213 0.3256 0.3413
4 0.4458 0.3245 0.3270 0.3230 0.3158 0.3231
5 0.5125 0.4838 0.4669 0.4041 0.4439 0.4407 
label <- c(Closest = "Closest Observation", NN = paste0("Nearest Neighbors (k=", k_opt, ")"),
           IDW = paste0("IDW (p=", p_opt, ")"), Kriging = "Ordinary Kriging",
           Ens_Equal = "Ensemble (Equal)", Ens_Weighted = "Ensemble (Weighted)")

rata_log <- colMeans(rmse_log)
sd_log   <- apply(rmse_log, 2, sd)
rata_ppm <- colMeans(rmse_ppm)
menang   <- table(factor(nama_cv[apply(rmse_log, 1, which.min)], levels = nama_cv))

ringkasan <- data.frame(Metode = label[nama_cv], RMSE_log = as.numeric(rata_log),
                        SD_antar_fold = as.numeric(sd_log), RMSE_ppm = as.numeric(rata_ppm),
                        Fold_terbaik = as.integer(menang))
ringkasan <- ringkasan[order(ringkasan$RMSE_log), ]
ringkasan <- cbind(Peringkat = seq_len(nrow(ringkasan)), ringkasan)
rownames(ringkasan) <- NULL
knitr::kable(ringkasan, digits = 4,
             col.names = c("Peringkat", "Metode", "RMSE rata-rata (log)", "SD antar fold",
                           "RMSE rata-rata (ppm)", "Jumlah fold terbaik"),
             caption = "Perbandingan RMSE rata-rata 5-fold Cross Validation")
Perbandingan RMSE rata-rata 5-fold Cross Validation
Peringkat Metode RMSE rata-rata (log) SD antar fold RMSE rata-rata (ppm) Jumlah fold terbaik
1 Ordinary Kriging 0.3396 0.0385 16.6472 3
2 Ensemble (Equal) 0.3528 0.0536 17.4384 1
3 Ensemble (Weighted) 0.3540 0.0494 17.3603 1
4 Nearest Neighbors (k=3) 0.3737 0.0675 18.2379 0
5 IDW (p=3) 0.3741 0.0541 18.4320 0
6 Closest Observation 0.4332 0.0616 19.9574 0 
cv_panjang <- data.frame(
  Fold   = rep(1:K, times = 6),
  Metode = factor(rep(label[nama_cv], each = K), levels = label[names(sort(rata_log))]),
  RMSE   = as.vector(rmse_log)
)
ggplot(cv_panjang, aes(Metode, RMSE, fill = Metode)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.12, size = 2, alpha = 0.8, show.legend = FALSE) +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 4, fill = "white") +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Sebaran RMSE 5-fold Cross Validation",
       subtitle = "Belah ketupat putih = rata-rata", x = NULL, y = "RMSE (skala log)") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

Interpretasi.

  1. Metode paling akurat adalah Ordinary Kriging dengan RMSE rata-rata 0.3396 (skala log) atau sekitar 16.6 ppm, unggul di 3 dari 5 fold. Urutan berikutnya: Ensemble (Equal) (0.3528), dan yang paling kurang akurat adalah Closest Observation (0.4332). Dibanding Closest Observation, RMSE metode terbaik lebih rendah sekitar 22%.
  2. Simpangan baku antar fold Ordinary Kriging adalah 0.0385, dibanding 0.0616 pada Closest Observation, jadi kualitas prediksinya juga lebih stabil.
  3. Ensemble (Equal) (RMSE 0.3528) lebih baik daripada NN dan IDW, sehingga penggabungan metode mengurangi kesalahan. Namun ensemble tidak mengalahkan Kriging karena Kriging sudah jauh lebih akurat daripada komponen lainnya, sehingga merata-ratakannya dengan metode yang lebih lemah justru menarik hasil ke arah kesalahan metode tersebut.
  4. Selisih RMSE ensemble weighted dan equal hanya 0.0012; bobot rata-rata antar fold: Closest = 0.181, NN = 0.253, IDW = 0.255, Kriging = 0.312.
  5. RMSE 0.34 pada skala log kira-kira setara kesalahan relatif sekitar 40% dari nilai Cu sebenarnya.

Kesimpulan

  1. Seluruh metode menghasilkan pola sebaran Cu yang samayakni tinggi di sepanjang Sungai Meuse dan menurun menjauhi sungai.
  2. Berdasarkan RMSE 5-fold Cross Validation, metode terbaik adalah Ordinary Kriging (RMSE rata-rata 0.3396 pada skala log), sedangkan yang paling kurang akurat adalah Closest Observation (0.4332).
  3. Kriging unggul karena autokorelasi spasial pada data ini jelas (variogram naik lalu mendatar) sehingga bobot dari semivariogram lebih tepat daripada bobot jarak sederhana, dan kriging juga menyediakan peta variansi prediksi.
  4. Ensemble menghasilkan prediksi yang stabil, tetapi tidak selalu lebih baik daripada metode tunggal terbaik.