tm_shape(d) +
  tm_dots(
    col = "lead",
    palette = "viridis",
    size = 0.3,
    title = "Lead"
  ) +
  tm_layout(
    title = "Persebaran Kandungan Lead",
    legend.outside = TRUE
  )

tm_shape(map) +
  tm_borders(col = "black") +
  tm_shape(coop) +
  tm_dots(
    col = "grey70",
    size = 0.02,
    alpha = 0.5
  ) +
  tm_shape(d) +
  tm_dots(
    col = "lead",
    palette = "viridis",
    size = 0.3,
    title = "Lead"
  ) +
  tm_layout(
    title = "Lokasi Data Pengamatan dan Titik Prediksi",
    legend.outside = TRUE
  )

summary(d$lead)

sd(d$lead)

var(d$lead)

ggplot(d, aes(x = lead)) +
  geom_histogram(
    bins = 20,
    fill = "steelblue",
    color = "black"
  ) +
  labs(
    title = "Distribusi Kandungan Timbal (Lead)",
    x = "Kandungan Lead",
    y = "Frekuensi"
  ) +
  theme_minimal()

# Menghitung experimental variogram
v_kriging <- variogram(lead ~ 1, d)

# Menampilkan experimental variogram
plot(
  v_kriging,
  main = "Experimental Variogram Lead"
)

# Menentukan nilai awal berdasarkan varians lead
nilai_var <- var(d$lead)

# Menyesuaikan model variogram Spherical
fv <- fit.variogram(
  object = v_kriging,
  model = vgm(
    psill = 0.8 * nilai_var,
    model = "Sph",
    range = 900,
    nugget = 0.2 * nilai_var
  )
)

# Menampilkan hasil fitting
fv

# Menampilkan plot experimental variogram dan model
plot(
  v_kriging,
  fv,
  main = "Fitting Variogram Lead"
)

# Membentuk model Ordinary Kriging
model_kriging <- gstat(
  formula = lead ~ 1,
  data = d,
  model = fv
)

# Melakukan prediksi pada titik prediksi
kpred <- predict(model_kriging, coop)

# Menyimpan hasil prediksi
kpred$pred <- kpred$var1.pred

# Membuat raster hasil prediksi
pred_kriging <- rasterize(
  vect(kpred),
  grid,
  field = "pred",
  fun = "mean"
)

c5 <- tm_shape(pred_kriging) +
  tm_raster(
    alpha = 0.6,
    palette = "viridis",
    title = "Prediksi Lead"
  ) +
  tm_layout(title = "Ordinary Kriging")

c5

# Menyimpan varians prediksi
kpred$varians <- kpred$var1.var

# Membuat raster varians prediksi
var_kriging <- rasterize(
  vect(kpred),
  grid,
  field = "varians",
  fun = "mean"
)

tm_shape(var_kriging) +
  tm_raster(
    palette = "viridis",
    title = "Varians"
  ) +
  tm_layout(
    title = "Varians Prediksi Kriging",
    legend.outside = TRUE
  )

# Menentukan seed
set.seed(123)

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

# Membagi data menjadi 5 fold
kf <- dismo::kfold(nrow(d), k = 5)

# Menyimpan nilai RMSE setiap metode
rmse1 <- rep(NA, 5) # Closest Observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest Neighbors
rmse4 <- rep(NA, 5) # Ensemble
rmse5 <- rep(NA, 5) # Kriging

for(k in 1:5){

  # Membagi data training dan testing
  test <- d[kf == k, ]
  train <- d[kf != k, ]

  # Memberikan ID pada data testing
  test$id <- seq_len(nrow(test))


  # 1. Closest Observation
  v_cv <- terra::voronoi(
    x = vect(train),
    bnd = vect(map)
  )

  v_cv <- st_as_sf(v_cv)

  # Menyimpan nilai lead sebagai prediksi
  v_cv$pred_co <- v_cv$lead

  # Mengambil prediksi pada lokasi data testing
  resp_cv <- st_intersection(
    v_cv["pred_co"],
    test["id"]
  )

  # Mencocokkan prediksi berdasarkan ID data testing
  p1 <- resp_cv$pred_co[
    match(test$id, resp_cv$id)
  ]

  # Memastikan seluruh data testing memperoleh prediksi
  stopifnot(
    length(p1) == nrow(test),
    !anyNA(p1)
  )

  rmse1[k] <- RMSE(test$lead, p1)


  # 2. IDW
  gs <- gstat(
    formula = lead ~ 1,
    locations = train,
    nmax = nrow(train),
    set = list(idp = 1)
  )

  p2 <- predict(gs, test)$var1.pred

  rmse2[k] <- RMSE(test$lead, p2)


  # 3. Nearest Neighbors
  nn <- gstat(
    formula = lead ~ 1,
    locations = train,
    nmax = 5,
    set = list(idp = 0)
  )

  p3 <- predict(nn, test)$var1.pred

  rmse3[k] <- RMSE(test$lead, p3)


  # 4. Ensemble (Equal Weights)
  weights <- c(1/3, 1/3, 1/3)

  p4 <- p1 * weights[1] +
        p2 * weights[2] +
        p3 * weights[3]

  rmse4[k] <- RMSE(test$lead, p4)


  # 5. Kriging
  vgram <- variogram(lead ~ 1, train)

  nilai_var_train <- var(train$lead)

  fv_cv <- fit.variogram(
    object = vgram,
    model = vgm(
      psill = 0.8 * nilai_var_train,
      model = "Sph",
      range = 900,
      nugget = 0.2 * nilai_var_train
    )
  )

  model_kriging_cv <- gstat(
    formula = lead ~ 1,
    data = train,
    model = fv_cv
  )

  p5 <- predict(model_kriging_cv, test)$var1.pred

  rmse5[k] <- RMSE(test$lead, p5)

}

hasil_rmse <- data.frame(
  Fold = 1:5,
  Closest_Observation = rmse1,
  IDW = rmse2,
  Nearest_Neighbors = rmse3,
  Ensemble = rmse4,
  Kriging = rmse5
)

hasil_rmse

rata_rmse <- data.frame(
  Metode = c(
    "Closest Observation",
    "IDW",
    "Nearest Neighbors",
    "Ensemble",
    "Kriging"
  ),
  RMSE = c(
    mean(rmse1),
    mean(rmse2),
    mean(rmse3),
    mean(rmse4),
    mean(rmse5)
  )
)

# Mengurutkan metode berdasarkan RMSE terkecil
rata_rmse <- rata_rmse[
  order(rata_rmse$RMSE),
]

rata_rmse