# Memuat pustaka yang diperlukan
library(sp)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)
library(gstat)
library(ggplot2)
library(mapview)
library(automap)

Data yang Digunakan

Data yang dipakai adalah meuse dari paket sp, berisi 155 titik sampel konsentrasi logam berat di tanah dekat Sungai Meuse, Belanda. Peubah yang dimodelkan adalah kadar zinc (ppm). Data pendukungnya meuse.grid, yaitu grid 40 m x 40 m yang menutupi wilayah studi dan dipakai sebagai lokasi prediksi.

# Memuat data meuse dari package sp
data(meuse, package = "sp")
data(meuse.grid, package = "sp")

# Mengubah data titik ke format SpatialPointsDataFrame
coordinates(meuse) <- ~x + y
proj4string(meuse) <- CRS("+init=epsg:28992")
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
# Mengubah data grid ke format SpatialPixelsDataFrame / SpatialGridDataFrame
coordinates(meuse.grid) <- ~x + y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) <- TRUE

# Konversi ke sf (Simple Features) untuk visualisasi
meuse_sf      <- st_as_sf(meuse)
meuse_grid_sf <- st_as_sf(meuse.grid)

# Data sampel (pengganti 'properties') dan peubah yang dimodelkan
d <- meuse_sf
d$vble <- d$zinc   # konsentrasi zinc (ppm)
# Ringkasan sebaran nilai zinc
summary(d$zinc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   113.0   198.0   326.0   469.7   674.5  1839.0

Nilai zinc di data ini berkisar dari 113 sampai 1839 ppm, dengan rata-rata 470 ppm dan median 326 ppm. Rata-rata yang lebih besar daripada median menandakan sebaran yang menjulur ke kanan: sebagian besar lokasi kadarnya rendah sampai sedang, dan hanya sedikit lokasi dengan kadar sangat tinggi. Hal ini nanti jadi alasan kriging dijalankan pada log(zinc).

# Visualisasi sebaran sampel data zinc
ggplot() +
  geom_sf(data = meuse_sf, aes(color = zinc, size = zinc)) +
  scale_color_viridis_c(option = "magma", name = "Zinc (ppm)") +
  theme_minimal() +
  labs(
    title = "Sebaran Lokasi Sampel Konsentrasi Zinc",
    subtitle = "Dataset Geostatistik Meuse",
    x = "Koordinat X", y = "Koordinat Y"
  )

Dari peta ini terlihat titik-titik dengan zinc tinggi (ukuran besar, warna terang) tidak tersebar acak, melainkan mengelompok, terutama di dekat sungai. Sebaliknya, titik yang lebih jauh dari sungai cenderung bernilai rendah. Artinya nilai di lokasi yang berdekatan mirip satu sama lain (ada autokorelasi spasial), dan itu syarat dasar supaya interpolasi masuk akal dilakukan.

# Batas wilayah studi (pengganti batas administratif Athena):
# gabungan sel grid meuse (resolusi 40 m) menjadi satu poligon
map <- meuse_grid_sf %>%
  st_buffer(dist = 20, endCapStyle = "SQUARE") %>%
  st_union() %>%
  st_sf()

# Batas untuk diagram Voronoi: convex hull yang mencakup seluruh titik sampel
# (supaya tidak ada sampel yang jatuh di luar batas)
bnd <- st_union(st_geometry(map), st_union(st_geometry(d))) %>%
  st_convex_hull() %>%
  st_buffer(100) %>%
  st_sf()

# Menampilkan peta zinc
tmap_mode("plot")
tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
  tm_dots("vble", palette = "viridis")

Poligon abu-abu adalah wilayah studi yang dibentuk dari gabungan sel grid, dan titik sampel berada di dalam atau di sekitar wilayah itu. Prediksi nanti hanya dibuat di dalam poligon ini, jadi semua yang dilakukan adalah interpolasi (menaksir di antara titik sampel), bukan ekstrapolasi.

# Membuat grid raster resolusi 40 m yang menutupi wilayah studi
grid <- rast(ext(vect(map)), resolution = 40, crs = "EPSG:28992")

# Lokasi prediksi = titik pusat sel grid meuse
coop <- st_sf(geometry = st_geometry(meuse_grid_sf))

# Menampilkan peta lokasi prediksi
qtm(coop)

Ada 3103 lokasi prediksi (titik pusat sel grid 40 m). Semua metode di bawah dijalankan pada lokasi yang sama, jadi perbedaan peta yang muncul murni berasal dari perbedaan metodenya.

Metode Deterministik

Closest Observation

Metode ini membagi wilayah dengan diagram Voronoi, satu poligon untuk satu titik sampel. Setiap lokasi di dalam poligon lebih dekat ke titik sampelnya daripada ke titik lain, sehingga nilai prediksinya sama dengan nilai sampel tersebut.

# Membuat diagram Voronoi dari lokasi sampel
v <- voronoi(x = vect(d), bnd = bnd)
plot(v)
# Menampilkan titik sampel
points(vect(d), cex = 0.5)

Diagram Voronoi membagi wilayah menjadi poligon-poligon, satu poligon untuk satu titik sampel. Setiap lokasi di dalam poligon itu lebih dekat ke titik sampelnya dibanding ke titik lain. Poligon kecil muncul di tempat sampel rapat, dan poligon besar di tempat sampel jarang.

# Prediction
v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)

Pada metode closest observation, prediksi di suatu lokasi sama dengan nilai sampel terdekatnya, sehingga satu poligon punya satu nilai. Hasilnya berbentuk blok-blok dengan perubahan nilai yang tajam di batas poligon, dan tidak ada nilai baru di luar nilai yang sudah teramati.

resp <- st_intersection(v, coop)
resp$pred <- resp$vble
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c1

Setelah dipotong dengan grid prediksi, peta raster ini masih terlihat berpetak-petak. Pola besarnya (zinc tinggi dekat sungai) sudah tertangkap, tetapi permukaannya kasar karena tidak ada penghalusan sama sekali.

Inverse Distance Weighting (IDW)

IDW menaksir nilai di lokasi baru sebagai rata-rata berbobot dari titik-titik sampel, dengan bobot yang menurun terhadap jarak. Makin dekat sebuah sampel, makin besar pengaruhnya. Besarnya penurunan diatur oleh pangkat jarak (idp), dan di sini dipakai idp = 1 dengan seluruh sampel dilibatkan.

# Menerapkan IDW menggunakan semua lokasi sampel
res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))

# Prediksi zinc untuk lokasi yang belum teramati
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

# Peta prediksi zinc dengan metode IDW
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c2

IDW menaksir nilai sebagai rata-rata berbobot dari seluruh sampel, dengan bobot sebesar 1/jarak (idp = 1), sehingga sampel yang dekat berpengaruh lebih besar. Karena nmax diisi sebanyak seluruh data dan pangkat jaraknya hanya 1, pengaruh sampel yang jauh tidak hilang begitu saja. Akibatnya peta lebih mulus dibanding Voronoi, tetapi nilai ekstrem cenderung “tertarik” ke arah rata-rata. Puncak zinc tinggi jadi terlihat lebih landai dan luas daripada kenyataan. Prediksi IDW juga tidak pernah keluar dari rentang nilai sampel.

Nearest Neighbors (NN)

NN menaksir nilai sebagai rata-rata sederhana dari k sampel terdekat, dengan bobot yang sama untuk semua tetangga. Pada kode, hal ini dicapai dengan nmax = 5 (k = 5) dan idp = 0 (jarak tidak ikut membobot).

res <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))

resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c3

Dengan nmax = 5 dan idp = 0, prediksi adalah rata-rata biasa dari 5 sampel terdekat tanpa memperhatikan seberapa dekat masing-masing. Hasilnya lebih halus dari Voronoi karena ada perataan, tetapi tetap muncul perubahan mendadak ketika kumpulan 5 tetangga berganti. Sifatnya sebagai penghalus lokal juga membuat nilai ekstrem teredam.

Ensemble

Ensemble menggabungkan prediksi dari beberapa metode dengan bobot tertentu supaya kelemahan satu metode dapat diimbangi metode lain. Metode yang digabung di sini adalah Closest Observation, IDW, dan NN. Untuk peta, bobotnya sama (1/3 masing-masing), sedangkan pada cross validation bobot ditentukan dari kebalikan RMSE.

# Closest observation (Voronoi = tetangga terdekat 1 titik)
# st_nearest_feature menjaga urutan baris sama dengan coop
p1 <- d$vble[st_nearest_feature(coop, d)]

# IDW
gs <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d),
            set = list(idp = 1))
p2 <- predict(gs, coop)$var1.pred
## [inverse distance weighted interpolation]
# Nearest neighbors
nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5,
            set = list(idp = 0))
p3 <- predict(nn, coop)$var1.pred
## [inverse distance weighted interpolation]
# Ensemble (bobot sama)
weights <- c(1/3, 1/3, 1/3)
p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
resp <- data.frame(x = st_coordinates(coop)[,1],
                   y = st_coordinates(coop)[,2],
                   pred = p4)

resp <- st_as_sf(resp, coords = c("x", "y"), crs = st_crs(map))
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c4

Ensemble menggabungkan ketiga prediksi di atas dengan bobot sama (masing-masing 1/3). Idenya, kelemahan satu metode bisa diimbangi metode lain: kekasaran Voronoi diperhalus oleh IDW dan nearest neighbors, sedangkan kecenderungan IDW meratakan nilai ekstrem sedikit dikoreksi oleh Voronoi yang memakai nilai asli sampel. Hasilnya adalah permukaan di antara ketiganya.

# Gabungkan semua plot dalam satu frame
tmap_arrange(c1, c2, c3, c4, ncol = 2)

Jika keempat peta dibandingkan, pola besarnya sama: zinc tinggi di sekitar sungai dan menurun ke arah yang lebih jauh. Perbedaannya ada pada tingkat kehalusan. Closest observation paling berpetak, IDW dan nearest neighbors lebih mulus, dan ensemble berada di tengah. Dari peta saja belum bisa dikatakan metode mana yang paling akurat, jadi perlu evaluasi kuantitatif.

Evaluasi Model

Kinerja tiap metode dinilai dengan K-fold cross validation (K = 5) dan RMSE. Data dibagi menjadi 5 bagian lewat dismo::kfold(). Bergiliran, satu bagian menjadi data uji dan empat bagian lainnya data latih, lalu prediksi dibandingkan dengan nilai sebenarnya. Cara ini menguji kemampuan model memprediksi data yang tidak ikut dipakai membangun model.

set.seed(123)

# Fungsi untuk menghitung RMSE
RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

# Pembagian data untuk K-fold
kf <- kfold(nrow(d), k = 5)

# Vektor penyimpan RMSE tiap metode
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Ensemble

for (k in 1:5) {
  # Split data menjadi test dan train
  test  <- d[kf == k, ]
  train <- d[kf != k, ]

  # Closest observation
  p1 <- train$vble[st_nearest_feature(test, train)]
  rmse1[k] <- RMSE(test$vble, p1)

  # IDW
  gs <- gstat(formula = vble ~ 1, locations = train,
              nmax = nrow(train), set = list(idp = 1))
  p2 <- predict(gs, test)$var1.pred
  rmse2[k] <- RMSE(test$vble, p2)

  # Nearest neighbors
  nn <- gstat(formula = vble ~ 1, locations = train,
              nmax = 5, set = list(idp = 0))
  p3 <- predict(nn, test)$var1.pred
  rmse3[k] <- RMSE(test$vble, p3)

  # Ensemble (bobot = invers RMSE, RMSE lebih kecil -> bobot lebih besar)
  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k])
  weights <- w / sum(w)
  p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
  rmse4[k] <- RMSE(test$vble, p4)
}
## [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]

Semakin kecil RMSE, semakin dekat prediksi ke nilai sebenarnya. Berbeda dengan peta ensemble sebelumnya, bobot ensemble di sini dibuat dari kebalikan RMSE tiap metode pada setiap lipatan, sehingga metode yang lebih akurat mendapat bobot lebih besar.

# RMSE untuk masing-masing dari 5 fold
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, ensemble = rmse4)
##   closest.obs      IDW nearest.neigh ensemble
## 1    309.9009 282.3744      265.3339 250.9364
## 2    385.8027 277.6231      275.7933 260.7226
## 3    323.5046 438.7651      364.2145 337.4435
## 4    237.6379 323.2493      255.5252 242.7105
## 5    299.4034 329.9387      294.7488 281.7473
# Rata-rata RMSE dari 5 fold
rata_rmse <- data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
                        nearest.neigh = mean(rmse3), ensemble = mean(rmse4))
rata_rmse
##   closest.obs      IDW nearest.neigh ensemble
## 1    311.2499 330.3901      291.1231 274.7121
terbaik <- names(rata_rmse)[which.min(unlist(rata_rmse))]

RMSE antar lipatan cukup bervariasi. Ini wajar karena nilai zinc ekstrem hanya ada di beberapa titik, sehingga lipatan yang kebetulan memuat titik-titik itu akan menghasilkan galat yang besar. Karena itu yang dibandingkan adalah rata-ratanya. Rata-rata RMSE terkecil dimiliki metode ensemble (274.7 ppm), sedangkan yang terbesar IDW (330.4 ppm). Artinya, secara rata-rata galat prediksi metode terbaik sekitar 275 ppm dari nilai zinc sebenarnya. Selisih RMSE antar metode juga perlu dilihat dalam konteks rata-rata zinc yang 470 ppm: jika selisihnya kecil, pemilihan metode tidak akan terlalu mengubah hasil.

Metode Stokastik

Kriging

Kriging adalah metode geostatistik yang memanfaatkan autokorelasi spasial lewat variogram, yaitu grafik yang menggambarkan bagaimana perbedaan nilai antar titik berubah menurut jarak. Model variogram (misalnya spherical, exponential, atau Gaussian) dicocokkan dengan variogram empiris, lalu dipakai untuk menentukan bobot sampel. Pada ordinary kriging, rata-rata populasi dianggap konstan tetapi tidak diketahui.

Tiga parameter pada model variogram:

  • Nugget: semivarians pada jarak mendekati nol, mewakili variasi skala mikro dan galat pengukuran.
  • Sill: nilai maksimum (dataran) semivarians, yaitu total varians spasial.
  • Range: jarak saat semivarians mencapai sill, batas jangkauan ketergantungan spasial.
# Peta interaktif sebaran zinc
mapview(meuse_sf, zcol = "zinc", map.types = "CartoDB.Voyager")

Peta interaktif ini mempertegas pengamatan sebelumnya, yaitu zinc tinggi berkumpul di sepanjang sungai. Kita bisa melihat posisi tiap titik di atas peta dasar untuk mengecek bahwa pola tersebut memang berkaitan dengan lokasi sungai.

# Variogram dengan cloud plot
vc <- variogram(log(zinc) ~ 1, meuse_sf, cloud = TRUE)
plot(vc)

Kriging memerlukan variogram, yaitu ukuran seberapa berbeda nilai dua titik sebagai fungsi jarak di antara keduanya. Peubah yang dipakai adalah log(zinc) untuk meredam kemiringan sebaran dan pengaruh nilai ekstrem. Variogram cloud menampilkan setiap pasangan titik, jumlahnya sangat banyak, sehingga plotnya terlihat menyebar dan sulit dibaca polanya. Gunanya lebih ke pengecekan awal, misalnya melihat apakah ada pasangan yang terlalu menyimpang.

# Variogram dengan plot standar
v <- variogram(log(zinc) ~ 1, data = meuse_sf)
plot(v)

Variogram standar merata-ratakan pasangan titik ke dalam kelompok jarak, sehingga polanya jelas. Semivarians kecil pada jarak dekat lalu naik seiring bertambahnya jarak sampai akhirnya mendatar. Artinya, titik yang berdekatan memang mirip, dan kemiripan itu hilang pada jarak tertentu. Bentuk seperti ini sesuai dengan model variogram yang memiliki sill dan range, misalnya model spherical.

vgm()
##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)
# Menampilkan daftar model variogram
show.vgms(par.strip.text = list(cex = 0.75))

Dua keluaran di atas hanya menunjukkan model variogram yang tersedia di gstat (Exp, Sph, Gau, dan lain-lain) beserta bentuk kurvanya. Dari bentuk variogram empiris sebelumnya, yang naik lalu mendatar secara jelas, model spherical dipilih sebagai kandidat.

# Model variogram awal
vinitial <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)

# Variogram empiris dengan model awal dan cutoff 1000
plot(v, vinitial, cutoff = 1000, cex = 1.5)

Model awal ini dibuat dari tebakan visual: nugget 0.1, partial sill 0.5, dan range 900 m. Kurvanya sudah cukup mengikuti titik-titik variogram empiris, tetapi belum tentu paling pas. Karena itu parameternya dioptimasi pada langkah berikutnya.

# Fitting variogram model "Sph"
fv <- fit.variogram(object = v,
                    model = vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
fv
##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044
plot(v, fv, cex = 1.5)

Hasil fitting memberikan nugget sebesar 0.051, partial sill 0.591, dan range 897 m. Nugget adalah variasi yang tidak dijelaskan oleh jarak (galat ukur atau variasi pada skala lebih kecil dari jarak antar sampel). Sill total (nugget + partial sill) adalah 0.641, dan range menunjukkan bahwa nilai dua titik masih saling berkaitan sampai kira-kira 897 m, lebih dari itu dianggap tidak berkorelasi. Rasio nugget terhadap sill sebesar 7.9%. Dengan patokan yang umum dipakai (di bawah 25% berarti ketergantungan spasial kuat), ketergantungan spasial pada data ini tergolong kuat, sehingga kriging layak dipakai. Kurva hasil fitting juga terlihat lebih menempel ke titik empiris dibanding model awal.

# Autofit variogram untuk data log-zinc
afvg <- autofitVariogram(log(zinc) ~ 1, meuse_sf)

# Plot hasil
plot(afvg)

autofitVariogram mencoba beberapa model secara otomatis dan memilih yang paling sesuai. Hasilnya dipakai sebagai pembanding fitting manual: kalau modelnya sama atau mirip dengan model spherical di atas, itu menguatkan pilihan sebelumnya. Perhatikan juga model yang dipilih dan parameternya pada judul plot.

# Objek gstat dengan model variogram hasil fitting (fv)
k <- gstat(formula = log(zinc) ~ 1, data = meuse_sf, model = fv)

# Prediksi pada grid
kpred <- predict(k, meuse_grid_sf)
## [using ordinary kriging]

Dengan formula ~ 1 (rata-rata dianggap konstan tapi tidak diketahui), yang dijalankan adalah ordinary kriging. Prediksi di tiap titik grid adalah kombinasi linear dari sampel di sekitarnya, dengan bobot yang ditentukan oleh model variogram fv. Selain prediksi, kriging juga menghasilkan varians prediksi.

# Prediksi log(zinc)
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = meuse_sf) +
  scale_color_viridis(name = "log(zinc)") + theme_bw()

Peta prediksi log(zinc) menunjukkan pola yang sejalan dengan metode deterministik: nilai tinggi di dekat sungai dan menurun ke arah yang lebih jauh. Permukaannya mulus dan tidak berpetak seperti Voronoi. Skala warnanya adalah log, jadi kalau ingin dibaca dalam ppm harus ditransformasi balik dengan exp().

# Varians prediksi
ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
  geom_sf(data = meuse_sf) +
  scale_color_viridis(name = "variance") + theme_bw()

Varians prediksi paling kecil di sekitar titik sampel dan membesar di lokasi yang jauh dari sampel atau berada di tepi wilayah dengan data yang jarang. Pola ini tidak bergantung pada nilai zinc, melainkan pada letak sampel dan model variogram. Ini keunggulan kriging dibanding metode deterministik, karena setiap prediksi disertai ukuran ketidakpastiannya, sehingga kita tahu bagian peta mana yang lebih bisa dipercaya.

Perbandingan kriging dengan metode deterministik

RMSE metode deterministik di atas dihitung pada skala ppm, sedangkan kriging dijalankan pada log(zinc), jadi angkanya belum bisa langsung dibandingkan. Supaya setara, kriging dievaluasi dengan 5 lipatan yang sama lalu hasil prediksinya ditransformasi balik ke ppm.

# Validasi silang kriging dengan lipatan yang sama (kf)
# Catatan: model variogram fv dipakai dari seluruh data
rmse5 <- rep(NA, 5)
for (j in 1:5) {
  test  <- d[kf == j, ]
  train <- d[kf != j, ]
  kk <- gstat(formula = log(zinc) ~ 1, data = train, model = fv)
  pk <- exp(predict(kk, test)$var1.pred)
  rmse5[j] <- RMSE(test$zinc, pk)
}
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
rata_rmse$kriging <- mean(rmse5)
rata_rmse
##   closest.obs      IDW nearest.neigh ensemble kriging
## 1    311.2499 330.3901      291.1231 274.7121 237.672
terbaik2 <- names(rata_rmse)[which.min(unlist(rata_rmse))]

Setelah kriging dimasukkan, metode dengan rata-rata RMSE terkecil adalah kriging (237.7 ppm), sedangkan kriging memiliki RMSE 237.7 ppm. Perlu diingat bahwa transformasi balik dari skala log bisa sedikit bias, dan model variogram dipakai dari seluruh data sehingga perbandingan ini sedikit menguntungkan kriging. Jadi perbandingannya sebaiknya dibaca sebagai gambaran kasar, bukan selisih yang pasti.

Kesimpulan

Konsentrasi zinc pada data meuse punya pola spasial yang jelas, yaitu tinggi di dekat sungai dan menurun menjauhinya, dan semua metode berhasil menangkap pola tersebut. Metode deterministik mudah dan cepat, tetapi hasilnya berbeda dalam tingkat kehalusan: Voronoi paling kasar, sedangkan IDW dan nearest neighbors lebih halus, dan ensemble menggabungkan ketiganya. Berdasarkan RMSE validasi silang, metode terbaik adalah kriging. Kriging punya nilai tambah berupa informasi varians prediksi yang tidak dimiliki metode deterministik, sehingga berguna kalau kita perlu tahu seberapa yakin hasil prediksi di setiap lokasi. Pemilihan akhirnya tergantung tujuan: kalau hanya butuh peta taksiran yang cepat, metode deterministik dengan RMSE terkecil sudah cukup, sedangkan kalau ketidakpastian juga penting, kriging lebih tepat.