TUGAS PRAKTIKUM 9

Buat interpolasi menggunakan data jura dengan seluruh metode interpolasi, yaitu Closest Observation, Nearest Neighbors (NN), IDW, Ensembel (bebas memilih metode yang digabung, minimal 3), dan Kriging. Setelah itu bandingkan RMSE menggunakan Cross Validation. Variabel yang diinterpolasi adalah log(zinc) sehingga RMSE semua metode berada pada skala yang sama.

# Memuat pustaka yang diperlukan
library(sp)
library(sf)

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(terra)

## terra 1.9.34

library(gstat)
library(dismo)

## Loading required package: raster

library(interp)   

## 
## Attaching package: 'interp'

## The following object is masked from 'package:dismo':
## 
##     circles

## The following object is masked from 'package:raster':
## 
##     area

## The following objects are masked from 'package:terra':
## 
##     area, cells

library(ggplot2)
library(viridis)

## Loading required package: viridisLite

library(mapview)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spData)
library(tmap)

## Registered S3 method overwritten by 'stars':
##   method                  from
##   st_interpolate_aw.stars sf

# Memuat data jura
data(jura, package = "gstat")

# Data pengamatan sebagai objek sf (koordinat km, tanpa CRS)
d <- st_as_sf(jura.pred, coords = c("Xloc", "Yloc"))

# Menambahkan kolom peubah yang diinterpolasi (log Zn)
d$vble <- log(d$Zn)

# Grid prediksi bawaan data jura
grd <- st_as_sf(jura.grid, coords = c("Xloc", "Yloc"))

# Batas area studi: convex hull dari seluruh titik pengamatan dan grid
map <- st_sf(geometry = st_convex_hull(st_union(c(st_geometry(d), st_geometry(grd)))))

nrow(d)

## [1] 259

summary(d$Zn)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.20   55.00   73.56   75.08   89.92  219.32

all(d$Zn > 0)

## [1] TRUE

hist(d$Zn, main = "Distribusi Zn")

hist(d$vble, main = "Distribusi Log(Zn)")

# Menampilkan peta log(Zn)
tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.

tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
  tm_dots("vble", palette = "viridis")

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

# Membuat grid raster 100x100 sel untuk menutupi area studi
grid <- rast(vect(map), nrows = 100, ncols = 100)

# Mengambil koordinat titik pusat dari setiap sel dalam grid
xy <- xyFromCell(grid, 1:ncell(grid))

# Membuat objek sf untuk lokasi prediksi
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))

# Menyaring lokasi prediksi yang berada dalam batas area studi
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deterministik

# Membuat diagram Voronoi dari lokasi pengamatan
v <- terra::voronoi(x = vect(d), bnd = vect(map))
plot(v)
# Menampilkan titik pengamatan di atas diagram Voronoi
points(vect(d), cex = 0.5)

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

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

idx <- st_nearest_feature(coop, d)
resp <- coop
resp$pred <- d$vble[idx]
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.

c1

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

resp <- predict(res, coop)

## [inverse distance weighted interpolation]

resp$pred <- resp$var1.pred
p3_grid <- resp$pred

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

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.

c3

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

# Mendapatkan prediksi untuk lokasi yang belum teramati
resp <- predict(res, coop)

## [inverse distance weighted interpolation]

resp$pred <- resp$var1.pred
p2_grid <- resp$pred

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

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.

c2

Metode Stokastik

# Memuat paket yang dibutuhkan
library(sp)
library(gstat)
library(sf)
library(mapview)

# Memuat dataset jura yang tersedia di paket gstat
data(jura, package = "gstat")

# Mengonversi menjadi objek sf dengan koordinat 'Xloc' dan 'Yloc' (satuan km, tanpa CRS)
jura.pred <- st_as_sf(jura.pred, coords = c("Xloc", "Yloc"))
jura.grid <- st_as_sf(jura.grid, coords = c("Xloc", "Yloc"))

# Peta interaktif (tanpa peta dasar karena data tidak punya CRS)
mapview(jura.pred, zcol = "Zn")

mapview(jura.grid)

# Variogram dengan cloud plot
vc <- variogram(vble ~ 1, d, cloud = TRUE)
plot(vc)

# Variogram dengan plot standar
vg <- variogram(vble ~ 1, data = d)
plot(vg)

# Daftar model variogram yang tersedia
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)

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

# Model variogram awal: psill, model, range, dan nugget
# (koordinat dalam km, sehingga range awal = 1 km; psill & nugget dari varians data)
vinitial <- vgm(psill = 0.6 * var(d$vble), model = "Sph",
                range = 1, nugget = 0.1 * var(d$vble))

# Visualisasi variogram empiris dengan model awal
plot(vg, vinitial, cex = 1.5)

# Fitting model "Sph" ke variogram empiris
fv <- fit.variogram(object = vg, model = vinitial)
fv

##   model      psill     range
## 1   Nug 0.03028104 0.0000000
## 2   Sph 0.12833218 0.7616774

plot(vg, fv, cex = 1.5)

library(automap)
# Autofit variogram untuk data log(Zn)
afvg <- autofitVariogram(vble ~ 1, as(d, "Spatial"))

# Plot hasil
plot(afvg)

library(ggplot2)

# Objek gstat dengan model variogram hasil fitting (fv)
k <- gstat(formula = vble ~ 1, locations = d, model = fv)

# Prediksi pada grid
kpred <- predict(k, coop)

## [using ordinary kriging]

p4_grid <- kpred$var1.pred

# Visualisasi prediksi log(Zn)
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = "log(Zn)") + theme_bw()

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

# Peta raster Kriging (agar sebanding dengan metode lainnya)
resp <- kpred
resp$pred <- resp$var1.pred
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.

c4

Pendekatan Ensembel

Pendekatan ensembel menggabungkan hasil prediksi dari berbagai metode interpolasi (Closest Observation, IDW, dan Nearest Neighbors) dengan memberikan bobot pada setiap metode. Bobot dapat ditentukan berdasarkan kinerja masing-masing metode.

Langkah-langkah: 1. Menghitung prediksi harga dengan masing-masing metode. 2. Menggabungkan prediksi setiap metode menggunakan bobot tertentu.

library(gstat)

# 1. Closest observation (Voronoi / tetangga terdekat)
p1_grid <- d$vble[st_nearest_feature(coop, d)]

# 2. IDW (semua titik, power = 2)
gs <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d),
            set = list(idp = 2))
p2_grid <- predict(gs, coop)$var1.pred

## [inverse distance weighted interpolation]

# 3. Nearest neighbors (k = 5, bobot sama)
nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5,
            set = list(idp = 0))
p3_grid <- predict(nn, coop)$var1.pred

## [inverse distance weighted interpolation]

# 4. Ordinary Kriging (variogram di-fit dari data)
vg  <- variogram(vble ~ 1, d)
fv  <- fit.variogram(vg, vgm(psill = 0.6 * var(d$vble), model = "Sph",
                             range = 1, nugget = 0.1 * var(d$vble)))
k   <- gstat(formula = vble ~ 1, locations = d, model = fv)
p4_grid <- predict(k, coop)$var1.pred

## [using ordinary kriging]

# 5. Ensemble (bobot sama)
weights <- c(1/4, 1/4, 1/4, 1/4)
p5 <- p1_grid * weights[1] + p2_grid * weights[2] +
      p3_grid * weights[3] + p4_grid * weights[4]

# Peta ensemble
resp <- data.frame(x = st_coordinates(coop)[, 1],
                   y = st_coordinates(coop)[, 2], pred = p5)
resp <- st_as_sf(resp, coords = c("x", "y"), crs = st_crs(map))
pred <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c5 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.

c5

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

# Evaluasi Model

Untuk menilai kinerja setiap metode dilakukan K-fold cross-validation dan dihitung RMSE (Root Mean Square Error). Cross-validation digunakan untuk memastikan model tidak hanya cocok dengan data latih tetapi juga dapat memprediksi data baru dengan baik.

Langkah-langkah:

1. Membagi data menjadi \(K\) bagian dengan dismo::kfold().
2. Menghitung RMSE untuk setiap metode dengan membandingkan nilai prediksi dan nilai sebenarnya.

set.seed(123)

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

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

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

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

  # Closest observation (tetangga terdekat, setara poligon Voronoi)
  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 = 2))
  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)

  # Kriging (variogram dibangun ulang dari data latih)
  vg_tr <- variogram(vble ~ 1, train)
  fv_tr <- suppressWarnings(
    fit.variogram(vg_tr, vgm(psill = 0.6 * var(train$vble), model = "Sph",
                             range = 1, nugget = 0.1 * var(train$vble)))
  )
  kr <- gstat(formula = vble ~ 1, locations = train, model = fv_tr)
  p4 <- predict(kr, test)$var1.pred
  rmse4[k] <- RMSE(test$vble, p4)

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

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

# RMSE pada masing-masing dari 5 fold
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, kriging = rmse4, ensemble = rmse5)

##   closest.obs       IDW nearest.neigh   kriging  ensemble
## 1   0.2705985 0.2229132     0.2888462 0.2392041 0.2311712
## 2   0.2906030 0.2270119     0.2760069 0.2459429 0.2309471
## 3   0.2642238 0.2806751     0.3165698 0.2943364 0.2729755
## 4   0.3405167 0.2356816     0.3027757 0.2903673 0.2674239
## 5   0.3724947 0.3068936     0.3604904 0.3274605 0.3251891

# Rata-rata RMSE dari 5 fold
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), kriging = mean(rmse4),
           ensemble = mean(rmse5))

##   closest.obs       IDW nearest.neigh   kriging  ensemble
## 1   0.3076873 0.2546351     0.3089378 0.2794622 0.2655414

Cross Validation (5-Fold)

Pada setiap fold, semua model dibangun ulang dari data latih (termasuk variogram untuk Kriging) lalu diuji pada data uji. Dua jenis ensembel dihitung:

• Ensembel bobot sama dari Closest Obs, IDW, NN, dan Kriging.
• Ensembel bobot invers RMSE dari Closest Obs, IDW, NN, dan Kriging, dengan bobot dihitung dari RMSE data latih melalui CV internal sederhana (leave-one-out) agar tidak bocor informasi dari data uji.

set.seed(123)

RMSE <- function(observed, predicted) sqrt(mean((observed - predicted)^2))

# Fungsi: empat prediksi dasar untuk data uji `test` dari data latih `train`
pred_base <- function(train, test) {
  p1 <- train$vble[st_nearest_feature(test, train)]

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

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

  vg  <- variogram(vble ~ 1, train)
  vgf <- suppressWarnings(
    fit.variogram(vg, vgm(psill = 0.6 * var(train$vble), model = "Sph",
                          range = 1, nugget = 0.1 * var(train$vble)))
  )
  p4 <- predict(gstat(formula = vble ~ 1, locations = train, model = vgf),
                test)$var1.pred

  cbind(closest = p1, IDW = p2, NN = p3, kriging = p4)
}

K  <- 5
kf <- kfold(nrow(d), k = K)

hasil <- data.frame(fold = 1:K, closest.obs = NA, IDW = NA, NN = NA,
                    kriging = NA, ens.equal = NA, ens.invRMSE = NA)

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

  P <- pred_base(train, test)
  hasil[k, c("closest.obs", "IDW", "NN", "kriging")] <-
    apply(P, 2, function(p) RMSE(test$vble, p))

  # Ensembel bobot sama
  hasil$ens.equal[k] <- RMSE(test$vble, rowMeans(P))

  # Bobot invers RMSE dari CV internal (5-fold) pada data latih
  inner <- dismo::kfold(nrow(train), k = 5)
  Pin <- matrix(NA, nrow(train), 4)
  for (j in 1:5) {
    Pin[inner == j, ] <- pred_base(train[inner != j, ], train[inner == j, ])
  }
  w <- 1 / apply(Pin, 2, function(p) RMSE(train$vble, p))
  w <- w / sum(w)
  hasil$ens.invRMSE[k] <- RMSE(test$vble, as.vector(P %*% w))
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

# RMSE pada tiap fold
hasil

##   fold closest.obs       IDW        NN   kriging ens.equal ens.invRMSE
## 1    1   0.2705985 0.2229132 0.2888462 0.2392041 0.2327443   0.2319686
## 2    2   0.2906030 0.2270119 0.2760069 0.2459429 0.2319678   0.2315743
## 3    3   0.2642238 0.2806751 0.3165698 0.2943364 0.2743982   0.2750806
## 4    4   0.3405167 0.2356816 0.3027757 0.2903673 0.2712854   0.2705801
## 5    5   0.3724947 0.3068936 0.3604904 0.3274605 0.3265961   0.3255304

# Rata-rata RMSE
rmse_rata <- colMeans(hasil[, -1])
rmse_rata

## closest.obs         IDW          NN     kriging   ens.equal ens.invRMSE 
##   0.3076873   0.2546351   0.3089378   0.2794622   0.2673984   0.2669468

Perbandingan RMSE

df_rmse <- data.frame(metode = names(rmse_rata), RMSE = as.numeric(rmse_rata))
df_rmse <- df_rmse[order(df_rmse$RMSE), ]
df_rmse

##        metode      RMSE
## 2         IDW 0.2546351
## 6 ens.invRMSE 0.2669468
## 5   ens.equal 0.2673984
## 4     kriging 0.2794622
## 1 closest.obs 0.3076873
## 3          NN 0.3089378

ggplot(df_rmse, aes(x = reorder(metode, RMSE), y = RMSE)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(RMSE, 3)), vjust = -0.3) +
  labs(x = "Metode", y = "Rata-rata RMSE (5-fold CV)",
       title = "Perbandingan RMSE Metode Interpolasi (data jura, log(zinc))") +
  theme_bw()

#Kesimpulan Berdasarkan hasil yang didapatkan, metode IDW memiliki rata-rata RMSE terkecil (0.255), selanjutnya diikuti ensembel bobot invers RMSE (0.267), ensembel bobot sama (0.267), Kriging (0.279), Closest Observation (0.308), dan Nearest Neighbors (0.309), sehingga IDW merupakan metode interpolasi terbaik karena bobotnya bergantung pada jarak atau titik pengamatan yang lokasinya lebih dekat dengan titik yang akan diprediksi, sedangkan Closest Observation dan NN hanya memakai satu atau beberapa tetangga tanpa pembobotan jarak. Ensemble juga menunjukkan hasil yang cukup baik karena mampu menghasilkan RMSE lebih rendah dibandingkan Ordinary Kriging, Closest Observation, dan Nearest Neighbors. Dengan demikian, IDW menjadi metode interpolasi yang paling optimal di antara metode yang diuji untuk memprediksi kadar Zn pada data Jura berdasarkan nilai RMSE hasil Cross Validation.