Preparation

library(sf) #Menangani data spasial
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(gstat) #Untuk NN, IDW, dan Kriging
## Warning: package 'gstat' was built under R version 4.6.1
library(ggplot2) #Untuk visualisasi data
library(dplyr) #Untuk manipulasi data
## Warning: package 'dplyr' was built under R version 4.6.1
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr) #Untuk struktur data
library(sp)
## Warning: package 'sp' was built under R version 4.6.1
data(meuse) 
data(meuse.grid)
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)       
grid_sf  <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

Mengubah data menjadi data spasial

Data Exploration

Kesimpulan

summary(meuse_sf$zinc)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   113.0   198.0   326.0   469.7   674.5  1839.0

Histogram

hist(meuse_sf$zinc, breaks = 20, main = "Zinc distribution", xlab = "Zinc (mg/kg)")

Scatter Plot

ggplot(meuse_sf) +
  geom_sf(aes(colour = zinc), size = 2) +
  scale_colour_viridis_c() +
  theme_minimal()

Cross-Validation

Closest Observation

cv_closest <- krige.cv(zinc ~ 1, meuse_sf, nmax = 1)
cv_closest

Proses Cross-Validation untuk metode Closest Observation (Nearest Neighbor) dilakukan menggunakan fungsi krige.cv() dengan opsi nmax = 1. Output menghasilkan objek sf berjumlah 155 titik yang memuat kolom var1.pred (nilai prediksi) dan observed (nilai aktual). Nilai var1.var dan zscore bernilai NA karena metode ini merupakan metode deterministic berbasis jarak terdekat yang tidak memperhitungkan struktur varians geostatistik. Kolom residual yang dihasilkan akan digunakan lebih lanjut untuk menghitung evaluasi kebaikan model, seperti Root Mean Square Error (RMSE).

NN

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))

k_vals <- seq(2, 14, by = 2)
rmse_k <- sapply(k_vals, function(k) {
  cv <- krige.cv(zinc ~ 1, meuse_sf, nmax = k, set = list(idp = 0))
  rmse(cv$observed, cv$var1.pred)
})

best_k <- k_vals[which.min(rmse_k)]
best_k
## [1] 4
plot(k_vals, rmse_k, type = "b", xlab = "k", ylab = "RMSE")

cv_nn <- krige.cv(zinc ~ 1, meuse_sf, nmax = best_k, set = list(idp = 0))
cv_nn

Untuk menentukan jumlah tetangga terdekat (\(k\)) optimal pada interpolasi k-Nearest Neighbor (k-NN), dilakukan pengujian cross-validation teriterasi (LOOCV) pada rentang \(k = 2\) hingga \(14\) dengan mematikan bobot jarak (idp = 0). Nilai evaluasi kebaikan model diukur menggunakan Root Mean Square Error (RMSE).Berdasarkan plot evaluasi parameter \(k\), diperoleh bahwa nilai \(k = 4\) menghasilkan error terendah (RMSE minimal sekitar 245). Penambahan jumlah tetangga melebihi \(k = 4\) menyebabkan kenaikan nilai RMSE secara monoton akibat efek over-smoothing dari lokasi pengamatan yang lebih jauh. Oleh karena itu, \(k = 4\) dipilih sebagai parameter terbaik (best_k) untuk pemodelan Nearest Neighbor.

IDW

idp_vals <- seq(0.5, 5, by = 0.5)
rmse_idp <- sapply(idp_vals, function(p) {
  cv <- krige.cv(zinc ~ 1, meuse_sf, set = list(idp = p))
  rmse(cv$observed, cv$var1.pred)
})

best_idp <- idp_vals[which.min(rmse_idp)]
best_idp
## [1] 3
cv_idw <- krige.cv(zinc ~ 1, meuse_sf, set = list(idp = best_idp))
cv_idw

Proses optimasi parameter Inverse Distance Weighting (IDW) dilakukan dengan menguji rentang nilai Inverse Distance Power (\(p\)) dari 0.5 hingga 5.0 menggunakan metode Leave-One-Out Cross-Validation (LOOCV). Nilai \(p\) optimal (best_idp) dipilih berdasarkan nilai RMSE terendah.Hasil cross-validation dengan parameter \(p\) optimal disimpan dalam objek cv_idw yang memuat 155 titik lokasi pengamatan. Nilai prediksi (var1.pred) dihitung berdasarkan kontribusi bobot jarak berkebalikan terhadap titik tetangga. Sama halnya dengan metode Nearest Neighbor, kolom var1.var dan zscore bernilai NA karena IDW tergolong sebagai metode deterministik. Nilai eror pada kolom residual selanjutnya dapat digunakan untuk evaluasi RMSE akhir model IDW.

Kriging

Variogram Empiris

v <- variogram(zinc ~ 1, meuse_sf)
plot(v)

Plot semivariogram eksperimental variabel zinc yang dihasilkan dari fungsi variogram(zinc ~ 1, meuse_sf) menggambarkan struktur autokorelasi spasial data. Terlihat bahwa semivarians meningkat tajam seiring bertambahnya jarak hingga mencapai titik pelandai (sill) di sekitar jarak \(900\) hingga \(1.100\) meter (range). Hal ini mengonfirmasi adanya ketergantungan spasial yang kuat pada skala lokal, di mana titik-titik sampel yang berdekatan memiliki karakteristik konsentrasi seng yang mirip. Plot empiris ini selanjutnya menjadi dasar untuk melakukan fitting model teoritis (seperti Spherical, Exponential, atau Gaussian) sebelum digunakan pada proses interpolasi Ordinary Kriging.

Model Teoritis

init <- vgm(psill = max(v$gamma) * 0.8,
            model = "Sph",
            range = max(v$dist) / 2,
            nugget = min(v$gamma))

v_fit <- fit.variogram(v, model = init)
v_fit
plot(v, v_fit)

Pemodelan variogram teoritis dilakukan menggunakan model Spherical (model = “Sph”) yang di-fit menggunakan metode Weighted Least Squares via fungsi fit.variogram(). Hasil fitting menghasilkan parameter struktur geostatistik berupa nilai Nugget (\(C_0\)) sebesar 24.800,95, Partial Sill (\(C\)) sebesar 134.745,26 (sehingga Sill Total \(C_0 + C = 159.546,21\)), serta jarak jangkauan spasial (range) sebesar 830,97 meter. Visualisasi plot plot(v, v_fit) menunjukkan bahwa kurva Spherical mampu merepresentasikan tren autokorelasi spasial data empiris secara akurat. Model variogram terpasang (v_fit) ini siap digunakan sebagai pembobot rasional dalam interpolasi Ordinary Kriging.

LOOCV

cv_krige <- krige.cv(zinc ~ 1, meuse_sf, model = v_fit)
cv_krige

Evaluasi kebaikan model Ordinary Kriging dilakukan melalui metode Leave-One-Out Cross-Validation (LOOCV) menggunakan fungsi krige.cv() dengan menyertakan model variogram Spherical yang telah di-fit (v_fit). Hasil pengujian disimpan dalam objek cv_krige yang memuat 155 lokasi sampel. Berbeda dengan metode deterministik seperti Nearest Neighbor dan IDW, Ordinary Kriging menghasilkan estimasi varians eror pada kolom var1.var serta nilai residual terstandardisasi pada kolom zscore. Keberadaan kolom zscore memfasilitasi evaluasi konsistensi ketidakpastian model geostatistik, sementara nilai pada kolom residual digunakan untuk menghitung nilai RMSE akhir pemodelan Ordinary Kriging.

Ensemble

ens_pred <- (cv_nn$var1.pred + cv_idw$var1.pred + cv_krige$var1.pred) / 3

r <- c(nn   = rmse(cv_nn$observed,    cv_nn$var1.pred),
       idw  = rmse(cv_idw$observed,   cv_idw$var1.pred),
       krig = rmse(cv_krige$observed, cv_krige$var1.pred))
w <- (1 / r) / sum(1 / r)

ens_pred_w <- w["nn"] * cv_nn$var1.pred +
              w["idw"] * cv_idw$var1.pred +
              w["krig"] * cv_krige$var1.pred
ens_pred_w
##   [1]  840.6244  795.3871  683.2033  442.9446  284.7712  249.4139  408.1628
##   [8]  599.8271  359.2549  234.9506  229.1274  282.8288  695.4083  537.2746
##  [15]  369.1887  697.3464  622.9691  564.0524  601.1452  663.1946  593.7352
##  [22]  555.1598  560.2541  202.5803  194.8217  205.1220  210.9249  230.8952
##  [29]  246.2986  248.4042  237.9838  233.4991  238.1647  228.3315  236.1808
##  [36]  331.4342  473.6181  878.4128  945.9911  775.6291  539.4412  274.4714
##  [43]  428.6031  325.5810  701.7043  575.1201  382.2879  256.9496  235.6151
##  [50]  220.1674  260.0639 1070.0175 1300.2677 1207.1811  950.0424  977.9985
##  [57]  576.1930  906.0160 1087.4339  644.7446  616.4526  522.8085  419.4545
##  [64]  425.3249  513.2865  671.9888  396.1112  291.3091  220.4536  833.0452
##  [71]  679.9365  720.9075  635.8796  546.4642  546.7292  386.4725  366.9983
##  [78]  661.4699  612.7993  770.8988  838.1638  488.3777  826.2369  292.6127
##  [85]  523.2390  259.6267  958.2194  804.8499  487.1911  412.8149  395.0276
##  [92]  462.1989  516.8256  228.6665  180.6401  216.1269  257.8184  440.7154
##  [99]  471.6309  331.2419  212.7925  417.3820  419.6762  143.1344  225.8137
## [106]  285.0476  154.7797  179.7715  318.9332  337.4048  271.5459  331.1620
## [113]  186.6584  220.0764  512.3580  486.0677  434.7992  966.7186  262.5600
## [120]  243.5176  436.4462  723.9004  885.2891  899.3552  962.6293  256.3503
## [127]  193.7233  239.9501  586.7938  636.1892  138.5509  315.8091  270.2845
## [134]  184.4796  173.2653  157.4303  173.4804  515.6618  488.1308  393.7814
## [141]  276.7348  211.9404  315.7309  379.5792  355.1169  440.4014  511.2454
## [148]  449.5808  435.6203  335.2028  278.6327  280.7974  386.9245  171.3966
## [155]  640.4728

Perbandingan RMSE

obs <- cv_krige$observed

results <- data.frame(
  Method = c("Closest Observation", paste0("NN (k=", best_k, ")"),
             paste0("IDW (p=", best_idp, ")"), "Ensemble (NN+IDW+Kriging)",
             "Ensemble (weighted)", "Ordinary Kriging"),
  RMSE = c(rmse(obs, cv_closest$var1.pred),
           rmse(obs, cv_nn$var1.pred),
           rmse(obs, cv_idw$var1.pred),
           rmse(obs, ens_pred),
           rmse(obs, ens_pred_w),
           rmse(obs, cv_krige$var1.pred))
) %>% arrange(RMSE)

results
ggplot(results, aes(x = reorder(Method, RMSE), y = RMSE)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(RMSE, 1)), hjust = -0.1) +
  coord_flip() + labs(x = NULL, y = "LOOCV RMSE") +
  theme_minimal()

Tahap akhir analisis dilakukan dengan membandingkan performa seluruh metode interpolasi spasial berbasis nilai Root Mean Square Error (RMSE) hasil Leave-One-Out Cross-Validation (LOOCV). Hasil pengujian diurutkan dari model dengan presisi tertinggi hingga terendah dan divisualisasikan menggunakan diagram batang ggplot2. Berdasarkan grafik evaluasi, Ordinary Kriging terpilih sebagai metode terbaik dengan nilai RMSE terendah yaitu 224,8. Pendekatan geostatistik ini terbukti paling superior dibanding metode deterministik karena mampu memodelkan autokorelasi spasial secara probabilistik melalui pemodelan variogram. Penggunaan pendekatan Ensemble (baik terbobot sebesar 235,0 maupun rata-rata sederhana sebesar 235,7) berhasil menjadi alternatif terkuat kedua yang mampu menstabilkan variansi eror prediksi. Sebaliknya, metode Closest Observation (\(nmax=1\)) menghasilkan tingkat eror tertinggi dengan RMSE sebesar 310,0. Hal ini mengonfirmasi bahwa mengandalkan satu titik pengamatan tunggal tanpa pembobotan jarak atau smoothing lokal sangat tidak direkomendasikan untuk pemodelan spasial konsentrasi seng pada dataset meuse.