batas <- st_sf(geometry = st_convex_hull(
  st_union(c(st_geometry(d), st_geometry(grid_sf)))
))

ggplot() +
  geom_sf(data = batas, fill = NA, color = "red") +
  geom_sf(data = grid_sf, color = "grey70", size = 0.3) +
  geom_sf(data = d, size = 0.8) +
  labs(title = "Batas wilayah, grid prediksi, dan titik pengamatan") +
  theme_bw()

# Root Mean Square Error
RMSE <- function(observed, predicted) sqrt(mean((observed - predicted)^2))

# 1. Closest observation (poligon Voronoi) ---------------------------------
pred_closest <- function(train, newdata) {
  v <- st_as_sf(voronoi(x = vect(train), bnd = vect(batas)))
  v <- st_sf(pred = v$vble, geometry = st_geometry(v))

  lokasi <- st_sf(.id = seq_len(nrow(newdata)), geometry = st_geometry(newdata))
  j <- st_join(lokasi, v, join = st_intersects)
  j <- j[!duplicated(j$.id), ]          # titik tepat di sisi poligon dihitung sekali
  j <- j[order(j$.id), ]
  p <- j$pred

  # cadangan: bila ada titik di luar semua poligon, pakai titik terdekat
  if (anyNA(p)) {
    na <- which(is.na(p))
    p[na] <- train$vble[st_nearest_feature(lokasi[na, ], train)]
  }
  p
}

# 2. Nearest neighbors: rata-rata k tetangga terdekat (idp = 0 -> bobot sama)
pred_nn <- function(train, newdata, k = 5) {
  g <- gstat(formula = vble ~ 1, locations = train, nmax = k, set = list(idp = 0))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# 3. IDW: seluruh titik dipakai, bobot 1 / jarak^idp -------------------------
pred_idw <- function(train, newdata, idp = 1) {
  g <- gstat(formula = vble ~ 1, locations = train,
             nmax = nrow(train), set = list(idp = idp))
  predict(g, newdata, debug.level = 0)$var1.pred
}

# 5. Ordinary kriging: semivariogram di-fit ulang pada data `train` ---------
fit_variogram <- function(train) {
  v <- variogram(vble ~ 1, data = train)
  fit.variogram(v, model = vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
}

pred_ok <- function(train, newdata) {
  g <- gstat(formula = vble ~ 1, data = train, model = fit_variogram(train))
  predict(g, newdata, debug.level = 0)
}

vor <- st_as_sf(voronoi(x = vect(d), bnd = vect(batas)))

ggplot() +
  geom_sf(data = vor, aes(fill = vble), color = "white", linewidth = 0.1, alpha = 0.8) +
  geom_sf(data = d, size = 0.6) +
  scale_fill_viridis(name = "log(zinc)") +
  labs(title = "Diagram Voronoi (Closest Observation)") +
  theme_bw()

p_closest <- pred_closest(d, grid_sf)

p_nn <- pred_nn(d, grid_sf, k = 5)

p_idw <- pred_idw(d, grid_sf, idp = 1)

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

v_emp <- variogram(vble ~ 1, data = d)
plot(v_emp)

v_init <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)
fv <- fit.variogram(v_emp, model = v_init)
fv

plot(v_emp, fv, cex = 1.5)

g_ok <- gstat(formula = vble ~ 1, data = d, model = fv)
kpred <- predict(g_ok, grid_sf, debug.level = 0)
p_ok <- kpred$var1.pred

rmse_fold <- sapply(names(metode), function(m)
  tapply(seq_len(nrow(oof)), oof$fold, function(i) RMSE(oof$obs[i], oof[[m]][i])))
colnames(rmse_fold) <- metode
rownames(rmse_fold) <- paste("Fold", 1:K)

kable(rmse_fold, digits = 4, caption = "RMSE tiap fold (skala log(zinc))")

hasil <- data.frame(Metode = metode,
                    RMSE_rata2 = colMeans(rmse_fold),
                    SD_antar_fold = apply(rmse_fold, 2, sd),
                    row.names = NULL) %>%
  arrange(RMSE_rata2) %>%
  mutate(Peringkat = row_number(),
         Selisih_vs_terbaik = RMSE_rata2 - min(RMSE_rata2)) %>%
  select(Peringkat, Metode, RMSE_rata2, SD_antar_fold, Selisih_vs_terbaik)

kable(hasil, digits = 4,
      caption = "Rata-rata RMSE 5-fold CV (skala log(zinc)); makin kecil makin baik")

urutan <- rev(hasil$Metode)   # terbaik di atas setelah coord_flip
rmse_long <- data.frame(Metode = factor(rep(metode, each = K), levels = urutan),
                        Fold = rep(1:K, times = length(metode)),
                        RMSE = as.vector(rmse_fold))
hasil_plot <- transform(hasil, Metode = factor(Metode, levels = urutan))

ggplot(hasil_plot, aes(x = Metode, y = RMSE_rata2)) +
  geom_col(fill = "steelblue", alpha = 0.8) +
  geom_point(data = rmse_long, aes(y = RMSE), shape = 21, fill = "white", size = 2) +
  geom_text(aes(y = 0, label = sprintf("%.3f", RMSE_rata2)), hjust = -0.2, color = "white", fontface = "bold", size = 3.8) +
  coord_flip() +
  labs(title = "Perbandingan RMSE (5-fold CV)",
       subtitle = "Batang = rata-rata 5 fold; titik = RMSE tiap fold",
       x = NULL, y = "RMSE log(zinc)") +
  theme_bw()

rmse_asli <- sapply(names(metode), function(m)
  tapply(seq_len(nrow(oof)), oof$fold,
         function(i) RMSE(exp(oof$obs[i]), exp(oof[[m]][i]))))
colnames(rmse_asli) <- metode

kable(data.frame(Metode = metode, RMSE_zinc_ppm = colMeans(rmse_asli), row.names = NULL) %>%
        arrange(RMSE_zinc_ppm),
      digits = 2, caption = "Rata-rata RMSE 5-fold CV pada skala zinc (ppm)")

oof_long <- do.call(rbind, lapply(names(metode), function(m)
  data.frame(Metode = metode[[m]], obs = oof$obs, pred = oof[[m]])))
oof_long$Metode <- factor(oof_long$Metode, levels = metode)

ggplot(oof_long, aes(obs, pred)) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  geom_point(alpha = 0.5, size = 1) +
  facet_wrap(~ Metode, ncol = 3) +
  labs(title = "Nilai observasi vs prediksi CV", x = "Observasi log(zinc)", y = "Prediksi log(zinc)") +
  theme_bw()