Import Library

library(sp)
## Warning: package 'sp' was built under R version 4.5.3
library(sf)
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(viridis)
## Warning: package 'viridis' was built under R version 4.5.2
## Loading required package: viridisLite

Data

data(meuse)
data(meuse.grid)

meuse$logzinc <- log(meuse$zinc)

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

obs_df <- as.data.frame(meuse)

ggplot(obs_df, aes(x, y, color = logzinc)) +
  geom_point(size = 2) +
  scale_color_viridis_c() +
  coord_equal() + theme_minimal() +
  labs(title = "Lokasi Pengamatan Data Meuse", color = "log(zinc)")

Titik pengamatan data Meuse tersebar mengikuti aliran Sungai Meuse, dan sebarannya tidak merata. Titik padat di beberapa area dan jarang di area lain. Karena itu, metode yang bergantung pada jarak dan kerapatan titik akan dipengaruhi pola sebaran ini. Nilai log(zinc) yang tinggi cenderung berada di dekat sungai.

Fungsi Peta

peta <- function(pred, judul, kolom = "var1.pred", legend = "log(zinc)") {
  df <- as.data.frame(pred)
  df$nilai <- df[[kolom]]
  ggplot(df, aes(x, y, fill = nilai)) +
    geom_raster() +
    scale_fill_viridis_c() +
    coord_equal() + theme_minimal() +
    labs(title = judul, fill = legend, x = NULL, y = NULL)
}

1. Closest Observation (Voronoi)

pts <- st_as_sf(meuse)
st_crs(pts) <- 28992
hull <- st_convex_hull(st_union(pts))
vor  <- st_sf(geometry = st_collection_extract(st_voronoi(st_union(pts)), "POLYGON"))
st_crs(vor) <- 28992
vor  <- st_join(vor, pts["logzinc"])
vor  <- suppressWarnings(st_intersection(vor, hull))

ggplot(vor) +
  geom_sf(aes(fill = logzinc), color = "white", linewidth = 0.1) +
  geom_sf(data = pts, size = 0.5) +
  scale_fill_viridis_c() + theme_minimal() +
  labs(title = "Diagram Voronoi (Closest Observation)", fill = "log(zinc)")

m_closest <- gstat(formula = logzinc ~ 1, locations = meuse,
                   nmax = 1, set = list(idp = 0))
p_closest <- predict(m_closest, meuse.grid)
## [inverse distance weighted interpolation]
peta(p_closest, "Closest Observation")

Peta hasil interpolasi menunjukkan pola berbentuk blok-blok poligon dengan batas yang tegas. Setiap lokasi memperoleh nilai prediksi dari satu titik pengamatan terdekat, sehingga seluruh area dalam satu poligon memiliki nilai yang sama. Kondisi ini menyebabkan pola spasial terlihat kasar dan tidak mengalami penghalusan. Selain itu, nilai ekstrem pada suatu titik pengamatan dapat memengaruhi seluruh wilayah poligon di sekitarnya. Metode ini hanya mempertimbangkan kedekatan lokasi tanpa memperhitungkan kontribusi beberapa titik tetangga maupun struktur autokorelasi spasial.

2. Nearest Neighbors (k = 5)

m_nn <- gstat(formula = logzinc ~ 1, locations = meuse,
              nmax = 5, set = list(idp = 0))
p_nn <- predict(m_nn, meuse.grid)
## [inverse distance weighted interpolation]
peta(p_nn, "Nearest Neighbors (k = 5)")

Hasil interpolasi Nearest Neighbors menunjukkan pola yang lebih halus dibandingkan Closest Observation. Hal ini disebabkan oleh penggunaan rata-rata lima titik pengamatan terdekat dengan bobot yang sama dalam menentukan nilai prediksi. Penggunaan beberapa titik tetangga dapat mengurangi pengaruh nilai ekstrem dari satu titik pengamatan. Meskipun demikian, pola interpolasi masih menunjukkan batas-batas wilayah yang relatif terlihat karena metode ini tidak mempertimbangkan perbedaan jarak antar-titik dalam pemberian bobot.

3. IDW (idp = 2)

m_idw <- gstat(formula = logzinc ~ 1, locations = meuse, set = list(idp = 2))
p_idw <- predict(m_idw, meuse.grid)
## [inverse distance weighted interpolation]
peta(p_idw, "Inverse Distance Weighting (idp = 2)")

Peta hasil interpolasi IDW menunjukkan pola yang lebih halus dan kontinu dibandingkan metode sebelumnya. Hal ini terjadi karena setiap titik pengamatan diberikan bobot berdasarkan kebalikan dari jaraknya terhadap lokasi prediksi, sehingga titik yang lebih dekat memiliki pengaruh lebih besar. Dengan menggunakan parameter pangkat jarak sebesar 2 (idp = 2), pengaruh titik pengamatan akan semakin berkurang seiring bertambahnya jarak. Namun, metode ini masih berpotensi menghasilkan pola melingkar di sekitar titik-titik pengamatan yang memiliki nilai ekstrem. Selain itu, hasil interpolasi sangat dipengaruhi oleh pemilihan parameter pangkat jarak yang dalam analisis ini ditentukan sebesar 2.

4. Ensemble

p_ens <- p_closest
p_ens$var1.pred <- (p_closest$var1.pred + p_idw$var1.pred + p_nn$var1.pred) / 3
peta(p_ens, "Ensemble (Closest + IDW + NN)")

Metode Ensemble menggabungkan hasil prediksi dari tiga metode interpolasi, yaitu Closest Observation, Nearest Neighbors, dan IDW, dengan memberikan bobot yang sama pada masing-masing metode. Hasil interpolasi menunjukkan pola yang relatif lebih halus dibandingkan Closest Observation karena prediksi tidak hanya bergantung pada satu titik pengamatan. Penggabungan ketiga metode juga menghasilkan pola prediksi yang merepresentasikan karakteristik dari masing-masing metode penyusunnya. Dengan demikian, Ensemble dapat memberikan hasil interpolasi yang lebih seimbang, meskipun kualitas prediksinya tetap bergantung pada performa metode-metode yang digunakan sebagai penyusun.

5. Ordinary Kriging

v_emp <- variogram(logzinc ~ 1, meuse)
v_fit <- fit.variogram(v_emp, vgm(psill = 0.6, model = "Sph", range = 900, nugget = 0.05))
v_fit
##   model      psill    range
## 1   Nug 0.05066522   0.0000
## 2   Sph 0.59061054 897.0412
plot(v_emp, v_fit, main = "Variogram Empiris dan Model Spherical")

p_krig <- krige(logzinc ~ 1, meuse, meuse.grid, model = v_fit)
## [using ordinary kriging]
peta(p_krig, "Ordinary Kriging")

peta(p_krig, "Ordinary Kriging: Varians Prediksi", kolom = "var1.var", legend = "Varians")

A. Interpretasi Model Variogram

Berdasarkan model variogram spherical yang terbentuk, diperoleh nilai nugget sebesar 0,051, partial sill sebesar 0,591, dan range sebesar 897,041 meter. Nilai-nilai tersebut dapat diinterpretasikan sebagai berikut:

  1. Range sebesar 897,041 meter menunjukkan bahwa ketergantungan spasial antar-pengamatan semakin melemah seiring bertambahnya jarak dan mencapai kondisi sill pada jarak sekitar 897 meter. Dengan demikian, pengamatan yang memiliki jarak lebih dekat cenderung memiliki hubungan spasial yang lebih kuat dibandingkan pengamatan yang berjauhan.

  2. Total sill sebesar 0,641 merupakan penjumlahan nugget dan partial sill. Sementara itu, nilai nugget hanya sekitar 7,9% dari total sill. Proporsi nugget yang relatif kecil menunjukkan bahwa sebagian besar variasi data dijelaskan oleh struktur variasi spasial, sehingga terdapat indikasi autokorelasi spasial yang cukup kuat pada data zinc.

Hasil ini mendukung penggunaan Ordinary Kriging karena metode tersebut dapat memanfaatkan struktur ketergantungan spasial yang teridentifikasi melalui variogram.

B. Interpretasi Hasil Interpolasi Ordinary Kriging

Hasil interpolasi Ordinary Kriging menunjukkan pola permukaan yang relatif halus dan kontinu. Berbeda dengan metode interpolasi sebelumnya, Kriging mempertimbangkan struktur ketergantungan spasial antar-pengamatan melalui model variogram dalam menentukan bobot prediksi. Selain menghasilkan nilai prediksi, metode ini juga memberikan informasi mengenai varians prediksi yang dapat digunakan untuk menggambarkan tingkat ketidakpastian hasil interpolasi. Varians prediksi cenderung lebih kecil pada lokasi yang dikelilingi titik pengamatan dengan jarak relatif dekat dan lebih besar pada lokasi yang jauh dari titik pengamatan. Hal ini menunjukkan bahwa tingkat ketidakpastian prediksi dipengaruhi oleh konfigurasi spasial titik-titik pengamatan.

Gabungan Semua Peta

gabung <- rbind(
  data.frame(as.data.frame(p_closest)[, c("x","y","var1.pred")], Metode = "Closest Observation"),
  data.frame(as.data.frame(p_nn)[, c("x","y","var1.pred")],      Metode = "Nearest Neighbors"),
  data.frame(as.data.frame(p_idw)[, c("x","y","var1.pred")],     Metode = "IDW"),
  data.frame(as.data.frame(p_ens)[, c("x","y","var1.pred")],     Metode = "Ensemble"),
  data.frame(as.data.frame(p_krig)[, c("x","y","var1.pred")],    Metode = "Ordinary Kriging")
)
ggplot(gabung, aes(x, y, fill = var1.pred)) +
  geom_raster() + scale_fill_viridis_c() + coord_equal() +
  facet_wrap(~Metode, ncol = 3) + theme_minimal() +
  labs(fill = "log(zinc)", x = NULL, y = NULL,
       title = "Perbandingan Peta Prediksi log(zinc)")

Hasil perbandingan kelima peta interpolasi menunjukkan adanya perbedaan karakteristik permukaan prediksi dari setiap metode. Closest Observation menghasilkan pola paling kasar dengan batas wilayah yang tegas, sedangkan Nearest Neighbors dan IDW menghasilkan pola yang relatif lebih halus. Metode Ensemble memperlihatkan pola yang menggabungkan karakteristik ketiga metode penyusunnya. Sementara itu, Ordinary Kriging menghasilkan permukaan prediksi yang relatif halus dan kontinu dengan mempertimbangkan struktur autokorelasi spasial. Perbedaan tersebut menunjukkan bahwa pemilihan metode interpolasi dapat memengaruhi pola distribusi nilai prediksi pada lokasi yang tidak memiliki pengamatan.

Cross Validation 5-fold

RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))

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

metode <- c("Closest Observation", "Nearest Neighbors", "IDW", "Ensemble", "Ordinary Kriging")
rmse_log  <- matrix(NA, K, length(metode), dimnames = list(paste0("Fold ", 1:K), metode))
rmse_orig <- rmse_log

for (k in 1:K) {
  train <- meuse[folds != k, ]
  test  <- meuse[folds == k, ]
  obs   <- test$logzinc

  p1 <- predict(gstat(formula = logzinc ~ 1, locations = train, nmax = 1, set = list(idp = 0)), test)$var1.pred
  p2 <- predict(gstat(formula = logzinc ~ 1, locations = train, nmax = 5, set = list(idp = 0)), test)$var1.pred
  p3 <- predict(gstat(formula = logzinc ~ 1, locations = train, set = list(idp = 2)), test)$var1.pred
  p4 <- (p1 + p2 + p3) / 3

  vg <- variogram(logzinc ~ 1, train)
  vf <- fit.variogram(vg, vgm(psill = 0.6, model = "Sph", range = 900, nugget = 0.05))
  p5 <- krige(logzinc ~ 1, train, test, model = vf, debug.level = 0)$var1.pred

  preds <- list(p1, p2, p3, p4, p5)
  for (j in seq_along(preds)) {
    rmse_log[k, j]  <- RMSE(obs, preds[[j]])
    rmse_orig[k, j] <- RMSE(exp(obs), exp(preds[[j]]))
  }
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
round(rmse_log, 4)
##        Closest Observation Nearest Neighbors    IDW Ensemble Ordinary Kriging
## Fold 1              0.5710            0.5261 0.5616   0.5057           0.4214
## Fold 2              0.4728            0.3685 0.4216   0.3714           0.2919
## Fold 3              0.7532            0.5292 0.5505   0.5537           0.4411
## Fold 4              0.7641            0.5305 0.5362   0.5593           0.4685
## Fold 5              0.6685            0.5223 0.5362   0.5134           0.4157
hasil <- data.frame(
  Metode         = metode,
  RMSE_log       = round(colMeans(rmse_log), 4),
  RMSE_zinc_asli = round(colMeans(rmse_orig), 2),
  row.names = NULL
)
hasil <- hasil[order(hasil$RMSE_log), ]
hasil$Peringkat <- seq_len(nrow(hasil))
hasil
##                Metode RMSE_log RMSE_zinc_asli Peringkat
## 5    Ordinary Kriging   0.4077         224.40         1
## 2   Nearest Neighbors   0.4953         276.40         2
## 4            Ensemble   0.5007         272.49         3
## 3                 IDW   0.5212         295.91         4
## 1 Closest Observation   0.6459         346.82         5

Interpretasi:

  1. Ordinary Kriging menghasilkan nilai RMSE terendah, yaitu 0,4077 pada skala log dan 224,40 pada skala zinc asli. Hasil ini menunjukkan bahwa Ordinary Kriging memiliki tingkat kesalahan prediksi paling rendah dibandingkan keempat metode lainnya. Keunggulan tersebut juga terlihat secara konsisten pada seluruh fold, di mana Ordinary Kriging selalu menghasilkan nilai RMSE terkecil. Dibandingkan Closest Observation, nilai RMSE pada skala log lebih rendah sekitar 36,9%, sedangkan dibandingkan Nearest Neighbors lebih rendah sekitar 17,7%. Hasil ini mengindikasikan bahwa pemanfaatan struktur autokorelasi spasial melalui model variogram mampu meningkatkan akurasi prediksi pada data zinc Meuse.
  2. Closest Observation menghasilkan nilai RMSE tertinggi, yaitu 0,6459 pada skala log dan 346,82 pada skala zinc asli. Hal ini menunjukkan bahwa metode tersebut memiliki tingkat kesalahan prediksi paling besar dibandingkan metode lainnya. Selain itu, nilai RMSE pada kelima fold memiliki rentang yang cukup lebar, yaitu 0,4728 hingga 0,7641, yang menunjukkan bahwa performa metode ini relatif lebih bervariasi terhadap pembagian data. Kondisi tersebut dapat disebabkan oleh karakteristik metode yang hanya menggunakan satu titik pengamatan terdekat dalam menentukan nilai prediksi, sehingga kurang mampu menangkap variasi spasial di sekitar lokasi prediksi.
  3. Nearest Neighbors dan Ensemble menghasilkan nilai RMSE yang relatif berdekatan, masing-masing sebesar 0,4953 dan 0,5007 pada skala log. Hal ini menunjukkan bahwa kedua metode memiliki tingkat akurasi prediksi yang relatif serupa. Pada skala log, Nearest Neighbors menghasilkan RMSE yang sedikit lebih rendah dibandingkan Ensemble. Namun, pada skala zinc asli, Ensemble justru menghasilkan RMSE yang sedikit lebih rendah, yaitu 272,49 dibandingkan 276,40. Perbedaan hasil tersebut menunjukkan bahwa urutan performa kedua metode dapat berubah bergantung pada skala pengukuran yang digunakan. Meskipun demikian, selisih RMSE yang relatif kecil menunjukkan bahwa performa keduanya tidak jauh berbeda pada data ini.
  4. IDW menghasilkan nilai RMSE sebesar 0,5212 pada skala log dan 295,91 pada skala zinc asli, sehingga menempati peringkat keempat. Hasil ini menunjukkan bahwa performa IDW masih berada di bawah Nearest Neighbors dan Ensemble pada skala log. Meskipun IDW mempertimbangkan perbedaan jarak dalam pemberian bobot, penggunaan parameter pangkat jarak sebesar 2 belum mampu menghasilkan prediksi yang lebih akurat dibandingkan rata-rata lima tetangga terdekat pada data ini.
  5. Berdasarkan hasil 5-fold Cross Validation, Ordinary Kriging merupakan metode dengan performa prediksi terbaik karena menghasilkan nilai RMSE terendah pada kedua skala pengukuran. Sementara itu, Closest Observation memiliki performa terendah dengan nilai RMSE terbesar. Nearest Neighbors dan Ensemble menunjukkan performa yang relatif berdekatan, sedangkan IDW berada di antara kedua kelompok tersebut. Hasil ini menunjukkan bahwa metode yang mempertimbangkan struktur ketergantungan spasial mampu memberikan hasil prediksi yang lebih akurat pada data zinc Meuse dibandingkan metode yang hanya mengandalkan kedekatan titik pengamatan.

Boxplot RMSE

long_cv <- data.frame(
  Fold   = rep(1:K, times = length(metode)),
  Metode = rep(metode, each = K),
  RMSE   = as.vector(rmse_log)
)
ggplot(long_cv, aes(reorder(Metode, RMSE, median), RMSE)) +
  geom_boxplot(outlier.alpha = 0.5) +
  geom_jitter(width = 0.1, alpha = 0.6) +
  theme_minimal(base_size = 12) +
  labs(title = "Perbandingan RMSE antar Metode (5-fold CV)", x = NULL, y = "RMSE (log zinc)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Boxplot menunjukkan bahwa Ordinary Kriging memiliki sebaran nilai RMSE paling rendah dibandingkan metode lainnya, yang mengindikasikan tingkat kesalahan prediksi relatif kecil dan performa yang lebih konsisten pada kelima fold. Sebaliknya, Closest Observation memiliki sebaran RMSE paling tinggi dan paling lebar, menunjukkan bahwa metode tersebut memiliki tingkat kesalahan prediksi lebih besar serta performa yang lebih bervariasi antar-fold. Sementara itu, Nearest Neighbors, IDW, dan Ensemble memiliki sebaran RMSE yang relatif berdekatan. Hasil ini memperkuat temuan bahwa Ordinary Kriging memberikan performa interpolasi terbaik pada data Meuse berdasarkan hasil 5-fold Cross Validation.

Kesimpulan

Berdasarkan hasil interpolasi dan evaluasi menggunakan 5-fold Cross Validation, diperoleh bahwa Ordinary Kriging merupakan metode interpolasi terbaik untuk data zinc Meuse dengan nilai RMSE terendah, yaitu 0,4077 pada skala log dan 224,40 pada skala zinc asli. Keunggulan ini menunjukkan bahwa pemanfaatan struktur autokorelasi spasial melalui model variogram dapat meningkatkan akurasi prediksi pada data yang dianalisis.

Metode Nearest Neighbors dan Ensemble menunjukkan performa yang relatif berdekatan, diikuti oleh IDW, sedangkan Closest Observation menghasilkan nilai RMSE tertinggi sehingga menjadi metode dengan performa prediksi paling rendah. Hasil ini juga menunjukkan bahwa penggunaan beberapa titik pengamatan dan pembobotan berdasarkan jarak maupun struktur ketergantungan spasial dapat memberikan hasil interpolasi yang lebih baik dibandingkan hanya menggunakan satu titik pengamatan terdekat.

Selain menghasilkan nilai prediksi, Ordinary Kriging juga memberikan informasi varians prediksi yang dapat digunakan untuk menggambarkan tingkat ketidakpastian hasil interpolasi. Dengan demikian, Ordinary Kriging menjadi metode yang paling sesuai di antara metode yang diuji pada data zinc Meuse dalam analisis ini.