library(sf)
library(terra)
library(tmap)
library(viridis)
library(dplyr)
library(gstat)
library(automap)
library(ggplot2)
library(spdep)
library(knitr)

set.seed(123)

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

0.1 Memuat data

Data berisi kejadian gempa bumi (magnitudo ≥ 4) di sekitar Indonesia. Koordinat awal berupa lintang-bujur (EPSG:4326). Karena IDW dan kriging bekerja berdasarkan jarak, koordinat diproyeksikan ke sistem koordinat planar (Lambert Azimuthal Equal Area berpusat di Indonesia) dengan satuan kilometer, sehingga jarak dan range variogram nanti langsung dalam km.

dat <- read.csv("dataregspas.csv")

d <- st_as_sf(dat, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

crs_km <- "+proj=laea +lat_0=-3 +lon_0=118 +datum=WGS84 +units=km +no_defs"
d <- st_transform(d, crs_km)

# Ambil kolom yang diperlukan
d <- d[, c("time", "place", "depth", "mag", "magType", "latitude", "longitude")]

# Peubah yang diinterpolasi (mengikuti penamaan 'vble' pada modul)
d$vble <- d$depth

kable(head(st_drop_geometry(d)), digits = 2)
time place depth mag magType latitude longitude vble
2026-11-10T08:34:07.131Z 58 km NE of Fef, Indonesia 35.00 5.0 mb -0.49 132.86 35.00
2026-11-10T01:53:24.195Z Pulau Pulau Tanimbar, Indonesia 10.00 5.4 mb -7.87 130.29 10.00
2026-11-09T14:55:47.985Z 164 km E of Masohi, Indonesia 10.00 4.3 mb -3.21 130.45 10.00
2026-11-09T11:25:07.473Z 271 km N of Tobelo, Indonesia 98.29 4.7 mb 4.17 127.78 98.29
2026-11-09T09:53:38.135Z 183 km E of Masohi, Indonesia 10.00 4.8 mb -2.98 130.58 10.00
2026-11-09T09:48:27.794Z 160 km E of Masohi, Indonesia 10.00 5.4 mww -3.06 130.39 10.00 
summary(st_drop_geometry(d)[, c("depth", "mag")])
##      depth             mag       
##  Min.   : 10.00   Min.   :4.000  
##  1st Qu.: 10.00   1st Qu.:4.300  
##  Median : 46.50   Median :4.450  
##  Mean   : 93.23   Mean   :4.552  
##  3rd Qu.:129.45   3rd Qu.:4.800  
##  Max.   :566.43   Max.   :5.700

Data terdiri dari 116 kejadian gempa yang tercatat pada 2026-10-11 sampai 2026-11-10. Kedalaman (depth) memiliki rata-rata 93.2 km dan median 46.5 km, dengan rentang 10–566.4 km. Rata-rata yang jauh di atas median menunjukkan sebaran yang menjulur ke kanan: sebagian besar gempa dangkal, sebagian kecil sangat dalam. Sebanyak 37 kejadian (32%) tercatat tepat 10 km, yang merupakan kedalaman tetap yang biasa diberikan USGS ketika kedalaman tidak dapat ditentukan dengan baik. Ini perlu diingat karena menambah “derau” pada data.

par(mfrow = c(1, 2))
hist(d$depth, breaks = 30, col = "steelblue", border = "white",
     main = "Kedalaman gempa", xlab = "Kedalaman (km)")
hist(d$mag, breaks = 15, col = "tomato", border = "white",
     main = "Magnitudo", xlab = "Magnitudo")

par(mfrow = c(1, 1))

0.2 Pemilihan peubah yang diinterpolasi

Data memiliki dua peubah numerik yang kandidat untuk diinterpolasi: kedalaman (depth) dan magnitudo (mag). Interpolasi hanya bermakna jika peubahnya memiliki autokorelasi spasial (asumsi dasar pada modul: lokasi berdekatan bernilai mirip). Hal ini diperiksa dengan Moran’s I (tetangga = 5 titik terdekat).

nb <- knn2nb(knearneigh(st_coordinates(d), k = 5))
lw <- nb2listw(nb, style = "W")

mi_depth <- moran.test(d$depth, lw)
mi_mag   <- moran.test(d$mag, lw)

moran_tab <- data.frame(
  peubah  = c("depth", "mag"),
  Moran_I = c(mi_depth$estimate[[1]], mi_mag$estimate[[1]]),
  E_I     = c(mi_depth$estimate[[2]], mi_mag$estimate[[2]]),
  p_value = c(format.pval(mi_depth$p.value, digits = 3),
              format.pval(mi_mag$p.value, digits = 3))
)
kable(moran_tab, digits = 4)
peubah Moran_I E_I p_value
depth 0.2707 -0.0087 2.58e-08
mag 0.0456 -0.0087 0.151

Interpretasi. Moran’s I kedalaman sebesar 0.271 (p = 2.58e-08), sedangkan magnitudo hanya 0.046 (p = 0.151). Artinya kedalaman gempa memiliki autokorelasi spasial positif yang jelas, sedangkan magnitudo nyaris tidak berstruktur spasial (gempa besar dan kecil tersebar acak). Hal ini masuk akal secara geologi: kedalaman hiposenter mengikuti geometri zona subduksi (lempeng yang menunjam makin dalam menjauhi palung), sehingga gempa yang berdekatan cenderung memiliki kedalaman serupa. Interpolasi magnitudo hanya akan menghasilkan permukaan yang mendekati rata-rata, jadi peubah yang diinterpolasi adalah kedalaman (depth, dalam km).

crs_hull <- st_convex_hull(st_union(d)) %>% st_buffer(dist = 150)  # penyangga 150 km
map <- st_sf(geometry = crs_hull)

tmap_mode("plot")
tm_shape(map) + tm_polygons(alpha = 0.3) +
  tm_shape(d) + tm_dots("vble", palette = "viridis", size = 0.6, title = "Kedalaman (km)") +
  tm_layout(main.title = "Sebaran kedalaman gempa", legend.outside = TRUE)

Interpretasi. Gempa dangkal dan dalam membentuk pola yang berkelompok mengikuti busur kepulauan (sepanjang Sumatra–Jawa–Nusa Tenggara, Banda, Maluku–Sulawesi Utara, dan Papua). Titik-titiknya tidak tersebar merata (padat di zona subduksi, jarang di wilayah lain), sehingga akurasi interpolasi akan bervariasi menurut kerapatan data.

1 Area dan Grid Prediksi

Batas wilayah dibuat dari convex hull seluruh titik gempa yang diberi penyangga 150 km (data tidak menyertakan batas administratif seperti depmunic pada modul). Kemudian dibuat grid raster 100 × 100 sel dan titik pusat sel dipakai sebagai lokasi prediksi.

grid <- rast(vect(map), nrows = 100, ncols = 100)

xy <- xyFromCell(grid, 1:ncell(grid))
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))
coop <- st_filter(coop, map)

qtm(coop)

Terdapat 7094 lokasi prediksi di dalam area. Karena wilayahnya memanjang timur–barat (± 5237 km) dan lebih sempit utara–selatan (± 2090 km), sel grid berbentuk persegi panjang. Area ini mencakup perairan dan sebagian wilayah negara tetangga karena memang gempa-gempa tersebut berada di sana, dan prediksi di luar kumpulan titik (misalnya tepi hull) merupakan ekstrapolasi sehingga kurang dapat diandalkan.

2 Fungsi Bantu dan Penalaan Parameter

Agar perbandingan di bagian evaluasi adil dan tidak “mengintip” data uji, semua parameter metode (pangkat jarak IDW, jumlah tetangga NN, model variogram, dan bobot ensemble) ditentukan hanya dari data latih. Fungsi-fungsi berikut dipakai berulang pada setiap fold.

# predict() pada gstat mencetak teks "[inverse distance weighted interpolation]" dll.;
# dibungkus agar tidak muncul di dokumen
pg <- function(g, new) {
  invisible(capture.output(o <- predict(g, new)))
  o
}

# (1) Penalaan pangkat jarak IDW (idp) dan jumlah tetangga NN (k) dengan leave-one-out
#     k = 1 identik dengan Closest Observation, sehingga NN dipilih dari k >= 2
p_grid <- seq(0.5, 6, by = 0.5)
k_grid <- 1:10

tune_params <- function(train) {
  zz <- train$vble
  n  <- length(zz)
  Dm <- as.matrix(dist(st_coordinates(train)))

  rmse_p <- sapply(p_grid, function(p) {
    W <- 1 / Dm^p
    diag(W) <- 0
    RMSE(zz, as.vector(W %*% zz) / rowSums(W))
  })

  Dn <- Dm
  diag(Dn) <- Inf
  ord <- t(apply(Dn, 1, order))
  rmse_k <- sapply(k_grid, function(k) {
    RMSE(zz, sapply(1:n, function(i) mean(zz[ord[i, 1:k]])))
  })

  list(p = p_grid[which.min(rmse_p)],
       k = k_grid[-1][which.min(rmse_k[-1])],
       rmse_p = rmse_p, rmse_k = rmse_k)
}

# (2) Pencocokan variogram otomatis (automap), dengan cadangan fit.variogram bila gagal
fit_vgm <- function(train) {
  out <- tryCatch(autofitVariogram(vble ~ 1, train)$var_model,
                  error = function(e) NULL)
  if (is.null(out)) {
    v0  <- var(train$vble)
    out <- fit.variogram(variogram(vble ~ 1, train),
                         vgm(psill = v0, model = "Sph", range = 500, nugget = 0.1 * v0))
  }
  out
}

# (3) Prediksi keempat metode dasar: Closest, NN, IDW, Kriging
predict_semua <- function(train, new, prm, vg) {
  cbind(
    closest.obs = train$vble[st_nearest_feature(new, train)],
    NN = pg(gstat(formula = vble ~ 1, locations = train,
                  nmax = prm$k, set = list(idp = 0)), new)$var1.pred,
    IDW = pg(gstat(formula = vble ~ 1, locations = train,
                   nmax = nrow(train), set = list(idp = prm$p)), new)$var1.pred,
    kriging = pmax(pg(gstat(formula = vble ~ 1, locations = train, model = vg),
                      new)$var1.pred, 0)   # kedalaman tidak mungkin negatif
  )
}

# (4) Bobot ensemble = kebalikan RMSE leave-one-out pada data latih (RMSE kecil -> bobot besar)
bobot_inner <- function(train, prm, vg) {
  n  <- nrow(train)
  pr <- matrix(NA_real_, n, 4)
  for (i in seq_len(n)) pr[i, ] <- predict_semua(train[-i, ], train[i, ], prm, vg)
  r <- apply(pr, 2, function(x) RMSE(train$vble, x))
  names(r) <- c("closest.obs", "NN", "IDW", "kriging")
  w <- (1 / r) / sum(1 / r)
  list(w = w, rmse = r)
}

# (5) Pembuat peta raster dengan skala warna yang sama untuk semua metode
brks <- c(-Inf, 25, 50, 100, 150, 200, 300, 400, Inf)

peta_raster <- function(p, judul) {
  r <- coop
  r$pred <- p
  pred <- rasterize(vect(r), grid, field = "pred", fun = "mean")
  tm_shape(pred) +
    tm_raster(alpha = 0.6, palette = "viridis", breaks = brks, title = "Kedalaman (km)") +
    tm_layout(main.title = judul, main.title.size = 0.9, legend.outside = TRUE)
}

2.1 Hasil penalaan pada seluruh data

par_all <- tune_params(d)

par(mfrow = c(1, 2))
plot(p_grid, par_all$rmse_p, type = "b", pch = 16,
     xlab = "Pangkat jarak (idp)", ylab = "RMSE leave-one-out (km)", main = "IDW")
abline(v = par_all$p, lty = 2, col = "red")
plot(k_grid, par_all$rmse_k, type = "b", pch = 16,
     xlab = "Jumlah tetangga (k)", ylab = "RMSE leave-one-out (km)", main = "Nearest Neighbors")
abline(v = par_all$k, lty = 2, col = "red")

par(mfrow = c(1, 1))

Interpretasi. Untuk IDW, RMSE leave-one-out terkecil dicapai pada pangkat jarak idp = 6 (RMSE 69.1 km); idp = 1 yang dipakai pada modul memberi RMSE 105 km. Nilai idp yang besar berarti pengaruh titik terdekat sangat dominan, yang menandakan kedalaman gempa berubah cukup cepat di ruang (struktur lokal kuat). Untuk NN, jumlah tetangga terbaik (k ≥ 2) adalah k = 2 (RMSE 85.4 km), sedangkan k = 1 (identik dengan Closest Observation) memberi 71.9 km. Jika RMSE membesar seiring bertambahnya k, rata-rata dari tetangga yang lebih jauh mulai mencampur gempa dengan kedalaman berbeda.

3 Metode Deterministik

3.1 Closest Observation

Metode ini membagi ruang menjadi poligon Voronoi di sekitar tiap titik gempa. Setiap lokasi prediksi mengambil nilai titik terdekatnya.

v <- terra::voronoi(x = vect(d), bnd = vect(map))
plot(v)
points(vect(d), cex = 0.5)

v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6, title = "Kedalaman (km)") +
  tm_layout(legend.outside = TRUE)

Pada modul, prediksi diperoleh dari st_intersection(v, coop). Di sini digunakan st_nearest_feature() yang hasilnya setara (titik dalam suatu poligon Voronoi pasti paling dekat ke titik pembentuknya), tetapi urutan hasilnya dijamin sama dengan urutan coop. Hal ini penting saat prediksi dari beberapa metode digabung pada ensemble.

p1 <- d$vble[st_nearest_feature(coop, d)]

c1 <- peta_raster(p1, "Closest Observation")
c1

Interpretasi. Peta berbentuk blok-blok tidak halus, karena tiap poligon Voronoi bernilai konstan sehingga nilai berubah mendadak di batas poligon. Metode ini tidak melakukan perataan: setiap prediksi persis sama dengan salah satu nilai observasi (rentang prediksi 10–566 km). Kelebihannya sederhana dan mempertahankan perubahan tajam (misalnya batas gempa dangkal-dalam), kelemahannya sangat sensitif terhadap satu titik tunggal dan membuat poligon besar di area yang datanya jarang (tepi wilayah).

3.2 Nearest Neighbors

Prediksi adalah rata-rata kedalaman dari k gempa terdekat dengan bobot sama (idp = 0). Nilai k = 2 diperoleh dari penalaan di atas.

nn <- gstat(formula = vble ~ 1, locations = d, nmax = par_all$k, set = list(idp = 0))
p2 <- pg(nn, coop)$var1.pred

c2 <- peta_raster(p2, paste0("Nearest Neighbors (k = ", par_all$k, ")"))
c2

Interpretasi. Dibandingkan Closest Observation, peta NN lebih halus karena merupakan rata-rata 2 tetangga, dan nilai ekstrem terpangkas (rentang prediksi 10–561 km). Kelemahannya, rata-rata dapat mencampur gempa dangkal dan dalam yang kebetulan berdekatan sehingga menghasilkan kedalaman “menengah” yang tidak pernah teramati, khususnya di zona subduksi dengan perubahan kedalaman yang curam.

3.3 Inverse Distance Weighting (IDW)

Setiap observasi diberi bobot \(1/d^{p}\) dengan \(d\) jarak ke lokasi prediksi. Pangkat jarak \(p\) = 6 (hasil penalaan), dan seluruh titik dipakai (nmax = nrow(d)).

idw <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = par_all$p))
p3 <- pg(idw, coop)$var1.pred

c3 <- peta_raster(p3, paste0("IDW (idp = ", par_all$p, ")"))
c3

Interpretasi. Hasil IDW lebih mulus daripada Closest Observation namun tetap mengikuti pola lokal data (rentang prediksi 10–566 km). Prediksi tidak pernah melampaui rentang observasi, sehingga puncak/lembah di antara titik tidak dapat diprediksi, dan dengan pangkat besar muncul efek “bull’s-eye” (lingkaran kecil) di sekitar titik data. IDW juga tidak memberi ukuran ketidakpastian prediksi.

4 Metode Stokastik: Kriging

Kriging memperhitungkan autokorelasi spasial melalui variogram. Berikut tahapannya (mengikuti modul): variogram awan, variogram empiris, pencocokan model, dan prediksi Ordinary Kriging (rata-rata konstan namun tidak diketahui).

4.1 Variogram

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

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

Pada variogram empiris, semivarians mulai dari nilai rendah pada jarak dekat lalu meningkat seiring jarak, menandakan struktur ketergantungan spasial. Selanjutnya tiga model umum (spherical, exponential, Gaussian) dicocokkan dan dibandingkan berdasarkan jumlah kuadrat galat (SSErr; makin kecil makin baik).

v0 <- var(d$vble)
ini <- list(
  Sph = vgm(psill = v0, model = "Sph", range = 500, nugget = 0.1 * v0),
  Exp = vgm(psill = v0, model = "Exp", range = 500, nugget = 0.1 * v0),
  Gau = vgm(psill = v0, model = "Gau", range = 500, nugget = 0.1 * v0)
)
fits <- lapply(ini, function(m) fit.variogram(v_emp, model = m))
sse  <- sapply(fits, function(f) attr(f, "SSErr"))

kable(data.frame(model = names(sse), SSErr = sse), digits = 1, row.names = FALSE)
model SSErr
Sph 248896
Exp 2913408
Gau 3982431 
plot(v_emp, fits[[names(which.min(sse))]], cex = 1.5,
     main = paste("Model terbaik (manual):", names(which.min(sse))))

Pencocokan otomatis dengan automap::autofitVariogram() juga dijalankan, dan model hasilnya yang dipakai untuk kriging, karena fungsi yang sama dipakai ulang pada tiap fold cross-validation.

afv <- autofitVariogram(vble ~ 1, d)
plot(afv)

fv <- afv$var_model
fv
##   model    psill    range
## 1   Nug     0.00   0.0000
## 2   Gau 16645.98 115.0923
nug  <- sum(fv$psill[fv$model == "Nug"])
sill <- sum(fv$psill)
rng  <- max(fv$range)

Interpretasi. Model yang terpilih memiliki nugget 0, sill 1.6646^{4} (km²), dan range 115.1 km. Rasio nugget/sill = 0%, sehingga ketergantungan spasial tergolong kuat (< 25% kuat, 25–75% sedang, > 75% lemah). Range 115 km berarti gempa yang terpisah lebih jauh dari jarak tersebut dianggap tidak lagi berkorelasi dalam kedalamannya. Nugget yang tidak nol mencerminkan variasi skala sangat kecil, ketidakpastian pengukuran kedalaman, serta kedalaman tetap 10 km pada banyak kejadian. Jika semivarians masih naik hingga jarak terjauh, itu menandakan adanya tren skala luas (kedalaman membesar ke arah tertentu menjauhi palung) yang tidak dimodelkan oleh Ordinary Kriging.

4.2 Prediksi Kriging

k <- gstat(formula = vble ~ 1, locations = d, model = fv)
kpred <- pg(k, coop)

p4 <- pmax(kpred$var1.pred, 0)   # kedalaman tidak mungkin negatif

c4 <- peta_raster(p4, "Ordinary Kriging")
c4

ggplot() +
  geom_sf(data = kpred, aes(color = var1.pred), size = 1.2, shape = 15) +
  geom_sf(data = d, size = 0.8) +
  scale_color_viridis(name = "Kedalaman (km)") + theme_bw() +
  labs(title = "Prediksi kedalaman (Ordinary Kriging) dan lokasi gempa")

ggplot() +
  geom_sf(data = kpred, aes(color = var1.var), size = 1.2, shape = 15) +
  geom_sf(data = d, size = 0.8) +
  scale_color_viridis(name = "Varians") + theme_bw() +
  labs(title = "Varians prediksi kriging")

Interpretasi. Kriging menghasilkan permukaan yang mulus dengan rentang prediksi 0–2355 km. Berbeda dari metode lain, kriging juga memberi varians prediksi: varians rendah (0 km²) di area yang padat titik gempa dan makin besar hingga 1.7013^{4} km² di area jauh dari data (tepi wilayah dan perairan yang sepi gempa). Peta varians ini menunjukkan di mana prediksi dapat dipercaya dan di mana tidak. Prediksi negatif (jika ada) dipotong ke 0 karena kedalaman tidak mungkin negatif.

5 Ensemble

Ensemble menggabungkan prediksi beberapa metode. Di sini digabungkan empat metode (Closest Observation, NN, IDW, dan Kriging) dengan dua skema bobot:

  1. Bobot sama (1/4 tiap metode), seperti contoh modul.
  2. Bobot berdasarkan kinerja: bobot ∝ 1/RMSE leave-one-out sehingga metode yang lebih akurat berbobot lebih besar. Bobot dihitung hanya dari data yang tersedia (pada cross-validation hanya dari data latih), tidak dari data uji.
P_grid <- cbind(closest.obs = p1, NN = p2, IDW = p3, kriging = p4)

inner_all <- bobot_inner(d, par_all, fv)
w_all <- inner_all$w

p5 <- rowMeans(P_grid)               # ensemble bobot sama
p6 <- as.vector(P_grid %*% w_all)    # ensemble berbobot

kable(data.frame(metode = names(w_all),
                 RMSE_LOO = round(inner_all$rmse, 2),
                 bobot = round(w_all, 3)), row.names = FALSE)
metode RMSE_LOO bobot
closest.obs 71.92 0.287
NN 85.41 0.241
IDW 69.12 0.298
kriging 118.61 0.174 
c5 <- peta_raster(p5, "Ensemble (bobot sama)")
c6 <- peta_raster(p6, "Ensemble (bobot 1/RMSE)")
c5

c6

Interpretasi. Bobot terbesar diberikan kepada IDW (29.8%) dan terkecil kepada kriging (17.4%), sesuai RMSE leave-one-out masing-masing. Karena bobot hanya berbeda tipis jika RMSE antar-metode mirip, kedua skema ensemble cenderung memberi peta yang serupa. Secara konsep, ensemble efektif bila galat tiap metode tidak sepenuhnya berkorelasi sehingga saling mengimbangi; namun bila satu metode jauh lebih baik, ensemble dapat tidak lebih baik daripada metode tunggal tersebut. Hal ini dibuktikan lewat cross-validation di bawah.

6 Perbandingan Peta Seluruh Metode

tmap_arrange(c1, c2, c3, c4, c5, c6, ncol = 2)

Interpretasi. Seluruh metode menangkap pola besar yang sama: kedalaman rendah di sebagian besar area dan kedalaman besar pada jalur-jalur tertentu yang berkaitan dengan zona subduksi. Perbedaannya ada pada kehalusan dan tingkat ekstrem: Closest Observation paling “kasar” dan berblok; NN dan IDW lebih halus; Kriging paling mulus dan mempertahankan variasi regional; kedua ensemble berada di antara metode-metode penyusunnya. Perbedaan paling besar antar-metode terjadi di area yang datanya jarang, yaitu di tempat ketidakpastian tertinggi.

7 Evaluasi dengan K-fold Cross-Validation

Data dibagi menjadi K = 5 bagian (dismo::kfold, set.seed(123)). Pada setiap fold, 1 bagian menjadi data uji dan 4 bagian menjadi data latih. Seluruh parameter (idp, k, model variogram, bobot ensemble) ditentukan ulang dari data latih saja. Sebagai pembanding ditambahkan baseline, yaitu memprediksi dengan rata-rata data latih (tanpa informasi spasial).

Catatan: pada modul, bobot ensemble per fold dihitung dari RMSE pada data uji itu sendiri, sehingga ensemble “melihat” data uji dan RMSE-nya terlalu optimis. Di sini bobot dihitung dari data latih sehingga evaluasi lebih jujur.

set.seed(123)
kf <- dismo::kfold(nrow(d), k = 5)

metode <- c("closest.obs", "NN", "IDW", "kriging",
            "ens.sama", "ens.berbobot", "baseline.mean")

rmse_fold  <- matrix(NA_real_, 5, length(metode),
                     dimnames = list(paste("fold", 1:5), metode))
bobot_fold <- matrix(NA_real_, 5, 4,
                     dimnames = list(paste("fold", 1:5),
                                     c("closest.obs", "NN", "IDW", "kriging")))
param_fold <- data.frame(fold = 1:5, idp = NA_real_, k = NA_real_, model_variogram = NA_character_)
oof <- matrix(NA_real_, nrow(d), length(metode), dimnames = list(NULL, metode))

for (f in 1:5) {
  idx   <- which(kf == f)
  test  <- d[idx, ]
  train <- d[-idx, ]

  prm  <- tune_params(train)
  vg_f <- fit_vgm(train)

  P  <- predict_semua(train, test, prm, vg_f)
  bw <- bobot_inner(train, prm, vg_f)

  pe <- rowMeans(P)
  pw <- as.vector(P %*% bw$w)
  P  <- cbind(P, ens.sama = pe, ens.berbobot = pw, baseline.mean = mean(train$vble))

  oof[idx, ]       <- P
  rmse_fold[f, ]   <- apply(P, 2, function(x) RMSE(test$vble, x))
  bobot_fold[f, ]  <- bw$w
  param_fold[f, 2:4] <- list(prm$p, prm$k,
                             paste(as.character(vg_f$model[vg_f$model != "Nug"]), collapse = "+"))
}

7.1 Parameter dan bobot tiap fold

kable(param_fold, row.names = FALSE)
fold idp k model_variogram
1 5 2 Gau
2 3 2 Gau
3 6 2 Sph
4 4 4 Gau
5 4 2 Gau 
kable(round(bobot_fold, 3))
closest.obs NN IDW kriging
fold 1 0.255 0.263 0.266 0.216
fold 2 0.260 0.224 0.300 0.216
fold 3 0.261 0.214 0.261 0.264
fold 4 0.306 0.322 0.317 0.055
fold 5 0.261 0.203 0.291 0.245

7.2 RMSE tiap fold

rmse_mean <- colMeans(rmse_fold)
rmse_sd   <- apply(rmse_fold, 2, sd)

tab_rmse <- rbind(rmse_fold, `rata-rata` = rmse_mean)
kable(round(tab_rmse, 2))
closest.obs NN IDW kriging ens.sama ens.berbobot baseline.mean
fold 1 36.21 77.01 36.24 68.80 39.36 39.50 143.84
fold 2 73.78 104.03 81.37 108.29 79.40 78.89 85.94
fold 3 46.86 38.45 41.40 53.53 37.24 37.66 77.57
fold 4 118.52 155.15 117.30 171.57 119.99 121.23 176.70
fold 5 75.65 77.83 69.46 196.71 89.00 88.23 102.20
rata-rata 70.20 90.49 69.15 119.78 73.00 73.10 117.25

7.3 Ringkasan perbandingan

m_main <- c("closest.obs", "NN", "IDW", "kriging", "ens.sama", "ens.berbobot")

rmse_pool <- apply(oof, 2, function(x) RMSE(d$vble, x))
mae_pool  <- apply(oof, 2, function(x) mean(abs(d$vble - x)))

ringkas <- data.frame(
  metode        = metode,
  RMSE_rata2    = round(rmse_mean, 2),
  SD_antarfold  = round(rmse_sd, 2),
  RMSE_gabungan = round(rmse_pool, 2),
  MAE_gabungan  = round(mae_pool, 2),
  perbaikan_vs_baseline_persen = round(100 * (1 - rmse_mean / rmse_mean["baseline.mean"]), 1)
)
ringkas <- ringkas[order(ringkas$RMSE_rata2), ]
kable(ringkas, row.names = FALSE)
metode RMSE_rata2 SD_antarfold RMSE_gabungan MAE_gabungan perbaikan_vs_baseline_persen
IDW 69.15 32.87 74.92 41.46 41.0
closest.obs 70.20 31.93 75.59 41.50 40.1
ens.sama 73.00 35.06 79.18 45.79 37.7
ens.berbobot 73.10 35.23 79.34 44.98 37.7
NN 90.49 43.05 97.99 55.52 22.8
baseline.mean 117.25 41.90 122.77 84.42 0.0
kriging 119.78 62.69 131.78 71.90 -2.2 
ggplot(data.frame(metode = factor(metode, levels = metode[order(rmse_mean)]),
                  rmse = rmse_mean, sd = rmse_sd),
       aes(x = metode, y = rmse, fill = metode)) +
  geom_col(show.legend = FALSE) +
  geom_errorbar(aes(ymin = rmse - sd, ymax = rmse + sd), width = 0.2) +
  scale_fill_viridis_d() + theme_bw() +
  labs(x = NULL, y = "RMSE rata-rata 5-fold (km)",
       title = "Perbandingan RMSE Cross-Validation",
       subtitle = "Batang galat = simpangan baku antar-fold")

m_plot <- c("closest.obs", "NN", "IDW", "kriging", "ens.berbobot")
df_oof <- data.frame(
  metode = factor(rep(m_plot, each = nrow(d)), levels = m_plot),
  obs    = rep(d$vble, length(m_plot)),
  pred   = as.vector(oof[, m_plot])
)
ggplot(df_oof, aes(obs, pred)) +
  geom_point(alpha = 0.6) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  facet_wrap(~metode) + theme_bw() +
  labs(x = "Kedalaman observasi (km)", y = "Kedalaman prediksi CV (km)",
       title = "Observasi vs prediksi (out-of-fold)")

Interpretasi.

  • Metode dengan RMSE rata-rata terkecil adalah IDW (69.15 km) dan yang terbesar kriging (119.78 km). Baseline (rata-rata saja) memiliki RMSE 117.25 km, sehingga metode terbaik menurunkan galat sebesar 41%. Seluruh metode interpolasi tidak seluruhnya lebih baik daripada baseline, yang mengonfirmasi bahwa informasi lokasi bermanfaat untuk memprediksi kedalaman gempa.
  • Di antara metode tunggal, yang terbaik adalah IDW (69.15 km). Ensemble berbobot (73.1 km) tidak lebih baik daripada metode tunggal terbaik, artinya penggabungan tidak menambah informasi karena galat antar-metode saling berkorelasi tinggi dan metode yang lemah ikut menarik hasil. Ensemble bobot sama menghasilkan 73 km.
  • Selisih RMSE antar-metode terbaik relatif kecil dibandingkan simpangan baku antar-fold (lihat batang galat), sehingga peringkat metode yang sangat berdekatan sebaiknya tidak dianggap berbeda nyata. Dengan hanya 116 observasi, hasil dapat berubah bila pembagian fold diubah.
  • Pada grafik observasi vs prediksi, titik yang makin dekat garis 45° berarti prediksi makin akurat. Galat terbesar umumnya terjadi pada gempa sangat dalam (> 300 km), yang jarang dan terisolasi sehingga tetangganya tidak mewakili kedalamannya.

8 Kesimpulan

  1. Peubah dan data. Dari dua peubah numerik, hanya kedalaman gempa yang memiliki autokorelasi spasial jelas (Moran’s I = 0.271; magnitudo 0.046), sehingga kedalaman dipilih sebagai peubah yang diinterpolasi pada 116 kejadian gempa (2026-10-11 s.d. 2026-11-10).
  2. Metode. Kelima metode pada tugas (Closest Observation, NN, IDW, Ensemble, Kriging) telah diterapkan pada grid 100 × 100. Ensemble menggabungkan empat metode dengan dua skema bobot. Parameter IDW (idp = 6) dan NN (k = 2) dituning dengan leave-one-out; kriging memakai variogram otomatis dengan nugget 0, sill 1.6646^{4}, dan range 115 km.
  3. Hasil cross-validation 5-fold. Metode terbaik berdasarkan RMSE rata-rata adalah IDW dengan RMSE 69.15 km, dan terburuk adalah kriging (119.78 km). Metode tunggal terbaik adalah IDW. Ensemble berbobot (73.1 km) tidak mengungguli seluruh metode tunggal.
  4. Rekomendasi. Untuk data ini, IDW direkomendasikan sebagai metode interpolasi kedalaman gempa. Jika diperlukan ukuran ketidakpastian, kriging unggul karena menyediakan peta varians prediksi, sedangkan metode lain tidak.
  5. Keterbatasan.
    • Datanya hanya 116 titik selama sekitar sebulan dan tersebar tidak merata, sehingga prediksi di wilayah yang jarang gempa adalah ekstrapolasi.
    • 37 kejadian berkedalaman tetap 10 km (nilai bawaan USGS) menambah derau.
    • Kedalaman gempa sebenarnya merupakan struktur tiga dimensi (gempa dangkal pada kerak dapat berada tepat di atas gempa dalam pada lempeng yang menunjam), padahal interpolasi 2D hanya memberi satu nilai per lokasi.
    • Ordinary Kriging mengasumsikan rata-rata konstan dan kemiripan ke segala arah (isotropik); kedalaman gempa memiliki tren regional dan arah dominan mengikuti busur subduksi, sehingga kriging dengan tren (universal kriging) atau variogram anisotropik berpotensi lebih baik.
    • Area prediksi berupa convex hull (bukan batas administratif) dan mencakup perairan serta wilayah negara tetangga.