Deskripsi Data

quakes berisi 1000 kejadian gempa dengan magnitudo > 4.0 di dekat Fiji sejak 1964.

data(quakes)
head(quakes)

##      lat   long depth mag stations
## 1 -20.42 181.62   562 4.8       41
## 2 -20.62 181.03   650 4.2       15
## 3 -26.00 184.10    42 5.4       43
## 4 -17.97 181.66   626 4.1       19
## 5 -20.42 181.96   649 4.0       11
## 6 -19.68 184.31   195 4.0       12

summary(quakes)

##       lat              long           depth            mag      
##  Min.   :-38.59   Min.   :165.7   Min.   : 40.0   Min.   :4.00  
##  1st Qu.:-23.47   1st Qu.:179.6   1st Qu.: 99.0   1st Qu.:4.30  
##  Median :-20.30   Median :181.4   Median :247.0   Median :4.60  
##  Mean   :-20.64   Mean   :179.5   Mean   :311.4   Mean   :4.62  
##  3rd Qu.:-17.64   3rd Qu.:183.2   3rd Qu.:543.0   3rd Qu.:4.90  
##  Max.   :-10.72   Max.   :188.1   Max.   :680.0   Max.   :6.40  
##     stations     
##  Min.   : 10.00  
##  1st Qu.: 18.00  
##  Median : 27.00  
##  Mean   : 33.42  
##  3rd Qu.: 42.00  
##  Max.   :132.00

Penanganan Titik Duplikat

Pada data quakes terdapat banyak titik dengan koordinat yang sama. Titik kembar menyebabkan matriks kriging singular dan jarak 0 pada IDW/NN, sehingga nilainya dirata-ratakan per lokasi.

cat("Jumlah baris awal            :", nrow(quakes), "\n")

## Jumlah baris awal            : 1000

cat("Jumlah lokasi unik (lat,long):", nrow(distinct(quakes, lat, long)), "\n")

## Jumlah lokasi unik (lat,long): 998

gempa <- quakes %>%
  group_by(lat, long) %>%
  summarise(depth = mean(depth), mag = mean(mag), .groups = "drop")

cat("Jumlah titik setelah agregasi:", nrow(gempa), "\n")

## Jumlah titik setelah agregasi: 998

Proyeksi Koordinat

Koordinat lon/lat diproyeksikan ke sistem planar (azimuthal equidistant, satuan kilometer) agar jarak Euclidean bermakna untuk variogram dan pembobotan.

crs_km <- "+proj=aeqd +lat_0=-21 +lon_0=179 +x_0=0 +y_0=0 +datum=WGS84 +units=km"

gempa_sf <- st_as_sf(gempa, coords = c("long", "lat"), crs = 4326) %>%
  st_transform(crs_km)

# Cek apakah masih ada titik duplikat setelah proyeksi
stopifnot(!any(duplicated(st_coordinates(gempa_sf))))
n <- nrow(gempa_sf)
n

## [1] 998

Grid Prediksi

Grid reguler 50 km x 50 km, dibatasi pada convex hull titik sampel agar tidak terjadi ekstrapolasi jauh di luar sebaran data.

hull <- st_convex_hull(st_union(gempa_sf))

grid_sf <- st_make_grid(gempa_sf, cellsize = 50, what = "centers") %>%
  st_sf(geometry = .) %>%
  st_filter(hull)

nrow(grid_sf)

## [1] 1643

Eksplorasi Singkat

ggplot(gempa_sf) +
  geom_sf(aes(color = depth), size = 1.5) +
  scale_color_viridis_c(option = "C", direction = -1) +
  labs(title = "Sebaran Kedalaman Gempa (km)", color = "Depth (km)") +
  theme_minimal()

hist(gempa$depth, breaks = 30, col = "steelblue", border = "white",
     main = "Histogram Kedalaman Gempa", xlab = "Depth (km)")

Histogram kedalaman bersifat bimodal (gempa dangkal dan gempa dalam), hal ini wajar karena adanya zona subduksi di wilayah Fiji-Tonga.

Closest Observation

Prediksi suatu lokasi = nilai pada satu titik sampel yang paling dekat.

pred_closest <- function(train, test) {
  idw(depth ~ 1, locations = train, newdata = test, nmax = 1, debug.level = 0)$var1.pred
}

NN (Nearest Neighbour, k tetangga terdekat)

Prediksi = rata-rata nilai k tetangga terdekat. Nilai k ditentukan lewat Cross Validation.

pred_nn <- function(train, test, k = 5) {
  idw(depth ~ 1, locations = train, newdata = test,
      nmax = k, idp = 0, debug.level = 0)$var1.pred
}

IDW

pred_idw <- function(train, test, p = 2) {
  idw(depth ~ 1, locations = train, newdata = test, idp = p, debug.level = 0)$var1.pred
}

Kriging

v_emp <- variogram(depth ~ 1, gempa_sf)

v_fit <- fit.variogram(
  v_emp,
  model = vgm(psill = var(gempa$depth) * 0.7, model = "Sph",
              range = 500, nugget = var(gempa$depth) * 0.15)
)
v_fit

##   model    psill    range
## 1   Nug     0.00    0.000
## 2   Sph 69453.93 1068.548

plot(v_emp, v_fit, main = "Empirical Variogram dan Model Spherical")

pred_krige <- function(train, test) {
  v <- variogram(depth ~ 1, train)
  m <- tryCatch(
    fit.variogram(v, vgm(psill = var(train$depth) * 0.7, model = "Sph",
                         range = 500, nugget = var(train$depth) * 0.15)),
    error = function(e) NULL, warning = function(w) NULL
  )
  if (is.null(m) || any(m$psill < 0) || m$range[2] <= 0) m <- v_fit
  krige(depth ~ 1, locations = train, newdata = test, model = m,
        debug.level = 0)$var1.pred
}

Ensemble

Ensemble dibuat dari NN (k-tetangga), IDW, dan Kriging. Dua skema:

• Ensemble rata-rata (bobot sama)
• Ensemble berbobot (bobot proporsional dengan 1/RMSE tiap metode)

ensemble_mean <- function(pn, pi, pk) (pn + pi + pk) / 3

ensemble_weighted <- function(pn, pi, pk, rmse_n, rmse_i, rmse_k) {
  w <- 1 / c(rmse_n, rmse_i, rmse_k)
  w <- w / sum(w)
  w[1] * pn + w[2] * pi + w[3] * pk
}

Fungsi RMSE dan Pembuatan Fold

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))

set.seed(123)
K <- 10
fold_id <- sample(rep(1:K, length.out = n))
table(fold_id)

## fold_id
##   1   2   3   4   5   6   7   8   9  10 
## 100 100 100 100 100 100 100 100  99  99

Pemilihan Parameter (k untuk NN dan power untuk IDW)

cv_generic <- function(fun, ...) {
  pred <- rep(NA_real_, n)
  for (f in 1:K) {
    tr <- gempa_sf[fold_id != f, ]
    te <- gempa_sf[fold_id == f, ]
    pred[fold_id == f] <- fun(tr, te, ...)
  }
  pred
}

# Tuning k (NN)
k_grid <- 2:20
rmse_k <- sapply(k_grid, function(k) rmse(gempa_sf$depth, cv_generic(pred_nn, k = k)))

# Tuning power (IDW)
p_grid <- seq(0.5, 4, by = 0.5)
rmse_p <- sapply(p_grid, function(p) rmse(gempa_sf$depth, cv_generic(pred_idw, p = p)))

k_opt <- k_grid[which.min(rmse_k)]
p_opt <- p_grid[which.min(rmse_p)]

par(mfrow = c(1, 2))
plot(k_grid, rmse_k, type = "b", pch = 19, xlab = "k", ylab = "RMSE CV",
     main = paste("Tuning NN (k optimal =", k_opt, ")"))
abline(v = k_opt, lty = 2, col = "red")
plot(p_grid, rmse_p, type = "b", pch = 19, xlab = "Power (p)", ylab = "RMSE CV",
     main = paste("Tuning IDW (p optimal =", p_opt, ")"))
abline(v = p_opt, lty = 2, col = "red")

par(mfrow = c(1, 1))

cat("k optimal untuk NN :", k_opt, "\n")

## k optimal untuk NN : 6

cat("power optimal IDW  :", p_opt, "\n")

## power optimal IDW  : 3

Hasil Prediksi CV Semua Metode

cv_closest <- cv_generic(pred_closest)
cv_nn      <- cv_generic(pred_nn,  k = k_opt)
cv_idw     <- cv_generic(pred_idw, p = p_opt)
cv_krige   <- cv_generic(pred_krige)

cv_ens_mean <- ensemble_mean(cv_nn, cv_idw, cv_krige)

r_nn <- rmse(gempa_sf$depth, cv_nn)
r_id <- rmse(gempa_sf$depth, cv_idw)
r_kr <- rmse(gempa_sf$depth, cv_krige)

cv_ens_w <- ensemble_weighted(cv_nn, cv_idw, cv_krige, r_nn, r_id, r_kr)

Perbandingan RMSE

obs <- gempa_sf$depth

hasil <- tibble(
  Metode = c("Closest Observation", paste0("NN (k=", k_opt, ")"),
             paste0("IDW (p=", p_opt, ")"), "Kriging (Ordinary)",
             "Ensemble Rata-rata (NN+IDW+Kriging)",
             "Ensemble Berbobot (NN+IDW+Kriging)"),
  RMSE = c(rmse(obs, cv_closest), r_nn, r_id, r_kr,
           rmse(obs, cv_ens_mean), rmse(obs, cv_ens_w)),
  MAE  = c(mean(abs(obs - cv_closest)), mean(abs(obs - cv_nn)),
           mean(abs(obs - cv_idw)),     mean(abs(obs - cv_krige)),
           mean(abs(obs - cv_ens_mean)), mean(abs(obs - cv_ens_w)))
) %>% arrange(RMSE)

kable(hasil, digits = 2, caption = paste0("Perbandingan RMSE dan MAE (", K, "-fold CV)"))

ggplot(hasil, aes(x = reorder(Metode, -RMSE), y = RMSE, fill = RMSE)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 1)), hjust = -0.1) +
  coord_flip() +
  scale_fill_viridis_c(direction = -1) +
  expand_limits(y = max(hasil$RMSE) * 1.1) +
  labs(x = NULL, y = "RMSE (km)", title = "Perbandingan RMSE Cross Validation") +
  theme_minimal()

df_cv <- tibble(Observasi = obs,
                `Closest Obs` = cv_closest, NN = cv_nn, IDW = cv_idw,
                Kriging = cv_krige, Ensemble = cv_ens_mean) %>%
  pivot_longer(-Observasi, names_to = "Metode", values_to = "Prediksi")

ggplot(df_cv, aes(Observasi, Prediksi)) +
  geom_point(alpha = 0.4, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = 2) +
  facet_wrap(~Metode) +
  labs(title = "Observasi dgn Prediksi (Cross Validation)") +
  theme_minimal()

kcv <- krige.cv(depth ~ 1, gempa_sf, model = v_fit, nfold = K)
cat("RMSE krige.cv (variogram global):", round(rmse(kcv$observed, kcv$var1.pred), 2), "\n")

## RMSE krige.cv (variogram global): 65.96

Hasil Utama

Ensemble menghasilkan nilai RMSE paling kecil dibandingkan metode lainnya. Ensemble Berbobot memiliki RMSE sebesar 64,17 km, sedangkan Ensemble Rata-rata sebesar 64,20 km. Perbedaan keduanya hanya sekitar 0,03 km, sehingga penggunaan bobot tidak memberikan peningkatan hasil yang berarti dibandingkan rata-rata biasa. Kedua metode ini menghasilkan prediksi yang lebih baik dibandingkan metode tunggal. Hal ini menunjukkan bahwa penggabungan prediksi NN, IDW, dan Kriging dapat mengurangi kesalahan karena setiap metode memiliki kelebihan dan kekurangan yang berbeda.

Closest Observation memiliki RMSE terbesar, yaitu 77,21 km, serta MAE terbesar sebesar 43,54 km. Hal ini terjadi karena metode tersebut hanya menggunakan satu titik terdekat sebagai dasar prediksi. Akibatnya, hasil prediksi sangat dipengaruhi oleh nilai dari titik tersebut. Metode ini juga menghasilkan pola peta yang terlihat seperti kumpulan bidang dengan batas-batas yang jelas (poligon Voronoi).

Perbandingan Metode Tunggal

• NN (k = 6) menjadi metode tunggal dengan hasil terbaik berdasarkan RMSE, yaitu 65,88 km. Metode ini menggunakan rata-rata dari enam titik terdekat untuk menghasilkan prediksi. Hasilnya cukup baik karena nilai kedalaman gempa pada titik-titik yang berdekatan kemungkinan memiliki pola yang relatif mirip.
• Kriging memiliki RMSE sebesar 67,68 km, tetapi nilai MAE-nya (38,84 km) lebih kecil dibandingkan NN (39,96 km) dan IDW (41,52 km). Artinya, kesalahan prediksi Kriging pada umumnya lebih kecil, tetapi terdapat beberapa prediksi dengan kesalahan yang cukup besar sehingga meningkatkan nilai RMSE-nya. Kelebihan lain dari Kriging adalah kemampuannya menghasilkan peta variansi yang dapat digunakan untuk melihat tingkat ketidakpastian hasil prediksi.
• IDW (p = 3) memiliki RMSE sebesar 68,44 km. Nilai power yang cukup tinggi menunjukkan bahwa titik yang lokasinya lebih dekat memiliki pengaruh lebih besar terhadap hasil prediksi dibandingkan titik yang lebih jauh.
• Perbedaan RMSE antara NN, Kriging, dan IDW hanya sekitar 1 sampai 3 km. Dengan demikian, ketiga metode tersebut memiliki hasil yang relatif berdekatan dan tidak menunjukkan perbedaan akurasi yang terlalu besar.

Catatan tentang Variogram

Variogram yang dihasilkan menggunakan model Spherical dengan nugget sebesar 0, sill sekitar 69.454, dan range sekitar 1.069 km. Hasil ini menunjukkan bahwa hubungan antar nilai kedalaman gempa masih dapat terlihat pada jarak yang cukup jauh.

RMSE dari krige.cv dengan menggunakan variogram yang dibuat dari seluruh data adalah 65,96 km. Nilai ini lebih kecil dibandingkan RMSE Kriging sebesar 67,68 km yang menggunakan variogram hasil perhitungan ulang pada setiap fold. Perbedaan ini dapat terjadi karena variogram global dibuat dengan menggunakan seluruh data, termasuk data yang nantinya menjadi data pengujian. Kondisi tersebut membuat hasil evaluasi terlihat sedikit lebih baik dari kondisi sebenarnya. Oleh karena itu, nilai RMSE sebesar 67,68 km lebih tepat digunakan untuk membandingkan Kriging dengan metode lainnya.

Keterbatasan

1. Titik gempa banyak terkumpul di sepanjang zona subduksi Fiji-Tonga, sehingga terdapat beberapa wilayah yang memiliki sedikit atau bahkan tidak memiliki titik pengamatan. Akibatnya, hasil prediksi pada wilayah tersebut memiliki tingkat keandalan yang lebih rendah.
2. Data depth memiliki dua kelompok nilai yang cukup berbeda, yaitu gempa dangkal dan gempa dalam. Kondisi ini membuat pola data kurang sesuai dengan asumsi bahwa hubungan spasial memiliki karakteristik yang relatif sama di seluruh wilayah. Hal ini kemungkinan menjadi salah satu alasan mengapa Kriging belum mampu menghasilkan akurasi yang lebih baik dibandingkan metode lainnya.
3. Nilai k pada NN dan power pada IDW ditentukan menggunakan fold yang sama dengan proses evaluasi. Hal ini dapat membuat nilai RMSE yang diperoleh sedikit lebih rendah daripada hasil yang mungkin didapatkan jika menggunakan data pengujian yang benar-benar terpisah.
4. Cross Validation secara acak pada data yang titiknya saling berdekatan dapat menghasilkan nilai RMSE yang cenderung lebih rendah. Hal ini terjadi karena titik yang digunakan sebagai data pengujian sering kali memiliki titik tetangga yang sangat dekat di dalam data pelatihan. Oleh karena itu, Spatial Cross Validation dapat menjadi alternatif yang lebih sesuai untuk mengevaluasi data yang memiliki pola penyebaran spasial.

Kesimpulan Akhir

Berdasarkan hasil 10-fold Cross Validation, Ensemble yang menggabungkan NN, IDW, dan Kriging menghasilkan akurasi terbaik dengan nilai RMSE paling kecil. Sementara itu, Closest Observation menghasilkan akurasi paling rendah dengan nilai RMSE paling besar. Jika hanya membandingkan metode tunggal, NN memberikan hasil terbaik, meskipun perbedaannya dengan Kriging dan IDW tidak terlalu besar. Di sisi lain, Kriging dapat menjadi pilihan apabila diperlukan informasi tambahan mengenai tingkat ketidakpastian hasil prediksi melalui peta variansi. Namun, seluruh hasil tersebut tetap perlu diperhatikan dengan mempertimbangkan kondisi data yang titik pengamatannya mengelompok dan memiliki dua kelompok kedalaman gempa yang berbeda.