1 Pendahuluan

Tugas ini membahas pendugaan titik (interpolasi spasial) menggunakan data meuse dengan tujuan:

Memprediksi konsentrasi seng (zinc) di seluruh wilayah studi dengan lima pendekatan: Closest Observation, Nearest Neighbors, Inverse Distance Weighting (IDW), Ensemble, dan Kriging.
Membandingkan akurasi kelima pendekatan dengan K-fold cross-validation memakai RMSE (Root Mean Squared Error).

Dasar metodenya adalah autokorelasi spasial, yaitu lokasi yang berdekatan cenderung memiliki nilai yang mirip. Karena itu, nilai di lokasi yang tidak diukur bisa ditaksir dari nilai di sekitarnya.

2 Paket dan Data

2.1 Library

library(sp)        # data meuse dan meuse.grid
library(sf)        # objek spasial sf
library(gstat)     # IDW, nearest neighbors, variogram, kriging
library(ggplot2)   # visualisasi
library(tidyr)     # reshape data untuk plot
library(dplyr)     # manipulasi data
library(knitr)     # tabel

Fungsi utama tiap paket: sp menyediakan data meuse dan meuse.grid; sf menangani objek spasial; gstat berisi semua metode interpolasi dan variogram; sisanya untuk visualisasi dan pengolahan tabel.

2.2 Memuat data meuse

Data meuse berisi 155 sampel tanah lapisan atas dari dataran banjir Sungai Meuse (Belanda). Peubah utama adalah zinc (ppm). Koordinat x dan y memakai sistem proyeksi RD New (EPSG: 28992), dalam satuan meter.

data(meuse)
data(meuse.grid)

str(meuse[, c("x", "y", "zinc")])

## 'data.frame':    155 obs. of  3 variables:
##  $ x   : num  181072 181025 181165 181298 181307 ...
##  $ y   : num  333611 333558 333537 333484 333330 ...
##  $ zinc: num  1022 1141 640 257 269 ...

summary(meuse$zinc)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   113.0   198.0   326.0   469.7   674.5  1839.0

str() memastikan ada 155 baris dan koordinat berupa angka, sehingga siap dijadikan objek spasial. summary() memberi ringkasan nilai zinc. Perhatikan bahwa nilai maksimum jauh lebih besar dari median, yang menandakan distribusi yang miring ke kanan.

2.3 Distribusi peubah zinc

par(mfrow = c(1, 2))
hist(meuse$zinc, main = "zinc (skala asli)", xlab = "zinc (ppm)",
     col = "#99F6E4", border = "#0D9488")
hist(log(meuse$zinc), main = "log(zinc)", xlab = "log(zinc)",
     col = "#5EEAD4", border = "#0D9488")

par(mfrow = c(1, 1))

Histogram skala asli menunjukkan ekor panjang ke kanan. Setelah ditransformasi dengan logaritma natural, bentuknya lebih simetris dan mendekati normal. Transformasi ini dipakai karena:

Kriging dan variogram lebih stabil jika distribusi peubah tidak terlalu miring. Nilai ekstrem yang besar bisa membuat variogram empiris sangat tidak stabil.
Log membuat hubungan antar titik lebih proporsional, sehingga selisih kecil di nilai rendah dan tinggi diperlakukan lebih setara.

Konsekuensinya, semua prediksi dihasilkan dalam skala log dan harus ditransformasi balik dengan exp() jika ingin dibaca dalam ppm. Transformasi ini tidak menjamin asumsi stasioneritas sepenuhnya terpenuhi. Asumsi itu tetap perlu diperiksa lewat variogram.

2.4 Sebaran titik sampel

ggplot(meuse, aes(x = x, y = y, color = zinc)) +
  geom_point(size = 2.5) +
  scale_color_gradient(low = "#CCFBF1", high = "#115E59", name = "zinc (ppm)") +
  coord_equal() +
  labs(title = "Sebaran lokasi pengukuran zinc", x = "Easting", y = "Northing") +
  theme_minimal(base_size = 12)

Titik-titik sampel tidak tersebar merata. Kepadatan tinggi di bagian tertentu, dan titik-titik dengan zinc tinggi cenderung berkumpul di area yang sama. Pola berkelompok ini menunjukkan adanya struktur spasial, yaitu nilai yang mirip saling berdekatan. Ini alasan utama metode spasial dipakai, bukan sekadar merata-ratakan semua nilai. Sebaran yang tidak merata juga penting untuk diingat, karena di area yang jarang sampel, prediksi biasanya kurang akurat.

3 Persiapan Data

meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, remove = FALSE)
meuse_sf$logzinc <- log(meuse_sf$zinc)

grid_sf <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)

cat("Jumlah titik data :", nrow(meuse_sf), "\n")

## Jumlah titik data : 155

cat("Jumlah titik grid :", nrow(grid_sf), "\n")

## Jumlah titik grid : 3103

Jumlah titik data harus 155. Jumlah titik grid adalah jumlah sel prediksi. Keduanya perlu diperiksa karena menentukan ukuran peta dan jumlah pasangan titik untuk cross-validation.

4 Metode Interpolasi

Setiap metode dijelaskan dulu secara teori, lalu diterapkan untuk memprediksi logzinc di seluruh grid. Peta ini hanya untuk visual. Evaluasi akurasi dilakukan terpisah di bagian Cross-Validation, karena prediksi di grid tidak bisa dibandingkan dengan nilai sebenarnya.

4.1 Closest Observation (Poligon Thiessen)

Metode ini membagi ruang menjadi poligon Voronoi (Thiessen). Setiap poligon berisi satu titik pengamatan. Setiap lokasi di dalam poligon mendapat nilai titik pengamatan tersebut. Lokasi di dalam poligon suatu titik adalah lokasi yang paling dekat dengan titik itu. Jadi prediksi Thiessen sama dengan pengamatan terdekat tunggal (1-NN).

\[ \hat{Z}(x_0) = Z(x_i) \quad \text{dengan } d(x_0, x_i) = \min_{j} d(x_0, x_j) \]

Implementasi memakai st_nearest_feature(), yang mengembalikan indeks titik terdekat. Cara ini setara dengan mencari poligon yang memuat titik prediksi, tetapi lebih sederhana dan tidak perlu membangun geometri poligon.

idx_co <- st_nearest_feature(grid_sf, meuse_sf)
grid_sf$CO <- meuse_sf$logzinc[idx_co]

summary(grid_sf$CO)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.727   5.118   5.412   5.697   6.315   7.517

idx_co berisi, untuk setiap sel grid, nomor titik sampel terdekatnya. grid_sf$CO lalu mengambil nilai logzinc dari titik tersebut. Akibatnya, nilai prediksi hanya ada sejumlah nilai dari 155 titik. Peta yang dihasilkan akan terlihat seperti mosaik, dengan blok-blok nilai konstan. Di batas dua poligon, nilai berubah secara mendadak. Inilah kelemahan utama metode ini. Ringkasan summary() menunjukkan rentang nilai prediksi, yang seharusnya tidak melampaui rentang logzinc data asli, karena metode ini tidak mengekstrapolasi.

4.2 Nearest Neighbors (rata-rata k tetangga)

Prediksi adalah rata-rata aritmetik dari $k$ titik terdekat dengan bobot sama.

\[ \hat{Z}(x_0) = \frac{1}{k} \sum_{i=1}^{k} Z(x_i) \]

g_nn_full <- gstat(formula = logzinc ~ 1, locations = meuse_sf,
                   nmax = 5, set = list(idp = 0))
pred_nn <- predict(g_nn_full, newdata = grid_sf, debug.level = 0)
grid_sf$NN <- pred_nn$var1.pred

summary(grid_sf$NN)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.821   5.296   5.597   5.690   6.069   7.173

gstat(..., nmax = 5, idp = 0) mencari 5 titik terdekat untuk setiap sel, lalu menghitung rata-ratanya. idp = 0 berarti bobot tidak bergantung jarak. Hasilnya lebih halus dibanding Thiessen, karena nilai prediksi berubah perlahan saat lokasi bergeser, dan tidak melompat di batas poligon. Namun, metode ini tetap tidak memperhitungkan jarak secara eksplisit. Titik ke-5 yang jauh diberi bobot sama dengan titik pertama yang dekat.

4.3 Inverse Distance Weighting (IDW)

Bobot berbanding terbalik dengan jarak dipangkatkan $p$. Di sini dipakai $p = 2$, pilihan yang paling umum.

\[ \hat{Z}(x_0) = \frac{\sum_{i=1}^{n} Z(x_i)\, d(x_0, x_i)^{-p}}{\sum_{i=1}^{n} d(x_0, x_i)^{-p}} \]

Semua titik dipakai (nmax = nrow).

g_idw_full <- gstat(formula = logzinc ~ 1, locations = meuse_sf,
                    nmax = nrow(meuse_sf), set = list(idp = 2))
pred_idw <- predict(g_idw_full, newdata = grid_sf, debug.level = 0)
grid_sf$IDW <- pred_idw$var1.pred

summary(grid_sf$IDW)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.791   5.484   5.694   5.777   6.041   7.482

Dengan $p = 2$, bobot titik turun sangat cepat seiring jarak. Titik sampel yang berjarak dua kali lebih jauh hanya memberi bobot seperempat. Akibatnya, nilai prediksi di sekitar titik sampel mendekati nilai titik itu, dan di tengah dua titik nilainya berada di antara keduanya. Peta IDW biasanya terlihat halus, tetapi di area yang jauh dari sampel, prediksi cenderung ditarik ke rata-rata lokal yang kurang informatif. Pilihan $p$ memengaruhi hasil: $p$ yang lebih besar membuat prediksi lebih “lokal”. Nilai $p$ yang sama sekali tidak diuji di sini. Itu bisa menjadi bahan diskusi tambahan.

4.4 Kriging (Ordinary Kriging)

Kriging memakai variogram untuk menentukan bobot. Langkahnya:

Hitung variogram empiris $\gamma(h) = \frac{1}{2N(h)} \sum [Z(x_i) - Z(x_i + h)]^2$.
Cocokkan model teoritis (Spherical, Exponential, Gaussian, dan lainnya) untuk memperoleh nugget $C_0$, sill $C_0 + C$, dan range $A$.
Selesaikan sistem Ordinary Kriging $\begin{bmatrix} C & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \mathbf{c}_0 \\ 1 \end{bmatrix}$ dengan kendala $\sum \lambda_i = 1$.

4.4.1 Variogram empiris

v_emp <- variogram(logzinc ~ 1, data = meuse_sf)
plot(v_emp, main = "Variogram empiris log(zinc)", col = "#0D9488", pch = 19)

Sumbu horizontal adalah jarak antar pasangan titik, dan sumbu vertikal adalah semivarian (setengah rata-rata selisih kuadrat). Pola yang diharapkan adalah semivarian kecil di jarak pendek, lalu naik dan mendatar. Jika semivarian naik sejak jarak paling kecil, itu menandakan variasi mikro atau galat pengukuran. Jika naik terus tanpa mendatar, itu tanda tren yang belum dimodelkan. Pada data ini, pola naik lalu mendatar terlihat jelas, yang menunjukkan adanya struktur spasial dan layak dilanjutkan ke kriging. Titik-titik variogram di jarak besar biasanya kurang andal karena jumlah pasangannya sedikit.

4.4.2 Pencocokan model variogram

v_init <- vgm(psill = 0.6, model = "Sph", range = 900, nugget = 0.05)
fv_full <- fit.variogram(v_emp, model = v_init)
fv_full

##   model      psill    range
## 1   Nug 0.05066522   0.0000
## 2   Sph 0.59061054 897.0412

plot(v_emp, fv_full, main = "Variogram empiris dan model Spherical",
     col = "#0D9488", pch = 19)

fit.variogram() mencari parameter model Spherical yang paling cocok dengan variogram empiris. Output tabel menampilkan:

psill (partial sill, $C$): selisih antara sill dan nugget. Ini bagian varians yang berasal dari struktur spasial.
range (A): jarak sampai semivarian mencapai sill. Di luar jarak ini, dua titik dianggap tidak saling memengaruhi.
Nugget (nug): semivarian saat jarak mendekati nol. Nilai nugget yang besar relatif terhadap sill berarti sebagian besar variasi bersifat acak atau terjadi pada skala lebih kecil dari jarak sampel terdekat.

Nilai awal v_init hanya titik awal optimasi. Jika hasil fitting terlihat tidak masuk akal, misalnya range sangat kecil atau psill negatif, nilai awal perlu diubah. Gambar memperlihatkan apakah kurva model mengikuti titik empiris. Perhatikan terutama jarak pendek, karena bagian itu paling menentukan bobot kriging untuk tetangga terdekat.

4.4.3 Prediksi Ordinary Kriging

krig_full <- krige(logzinc ~ 1, locations = meuse_sf, newdata = grid_sf,
                   model = fv_full, debug.level = 0)
grid_sf$KRIG <- krig_full$var1.pred
grid_sf$KRIG_var <- krig_full$var1.var

summary(grid_sf$KRIG)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.777   5.238   5.573   5.707   6.172   7.440

krige() menghitung prediksi var1.pred dan varians kriging var1.var untuk setiap sel. Varians ini adalah ukuran ketidakpastian prediksi. Varians kecil di dekat titik sampel dan membesar di area yang jauh dari sampel. Ini sifat khas kriging, dan tidak dimiliki IDW atau NN. Kriging juga memberi bobot yang mempertimbangkan kelompok titik. Misalnya, dua titik yang berdekatan satu sama lain tidak akan menggandakan pengaruhnya secara berlebihan. Keunggulan ini membuat kriging biasanya lebih akurat jika variogram terpasang dengan baik.

4.5 Ensemble

Ensemble menggabungkan prediksi beberapa metode dengan bobot $w_i$, di mana $\sum w_i = 1$.

\[ \hat{Z}_{ens}(s_0) = \sum_{i=1}^{M} w_i \, \hat{Z}_i(s_0) \]

Dipakai empat metode dasar (Closest Observation, NN, IDW, dan Kriging) dengan dua skema bobot:

ENS_EQ: bobot sama, $w_i = 1/4$. Bobot tidak dipilih dari data, sehingga tidak ada kebocoran informasi dari evaluasi.
ENS_W: bobot berbanding terbalik dengan RMSE cross-validation, $w_i \propto 1/\text{RMSE}_i$. Bobot ini baru bisa dihitung setelah cross-validation, sehingga dijelaskan di bagian berikutnya.

grid_sf$ENS_EQ <- rowMeans(st_drop_geometry(grid_sf)[, c("CO", "NN", "IDW", "KRIG")])
summary(grid_sf$ENS_EQ)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.789   5.317   5.576   5.718   6.166   7.395

rowMeans() menghitung rata-rata keempat prediksi untuk setiap sel grid. Karena bobotnya sama, ENS_EQ menyeimbangkan kelebihan dan kekurangan tiap metode. Thiessen yang kasar tetapi lokal, IDW yang halus, NN yang rata-rata lokal, dan kriging yang memakai struktur spasial akan saling mengimbangi. Ensemble dengan bobot sama tidak mengungguli kriging pada data ini. Nilai RMSE-nya (0,476) berada di bawah NN dan IDW, tetapi masih di atas kriging (0,412). Ini menunjukkan bahwa menggabungkan metode yang lebih lemah tidak selalu menghasilkan prediksi yang lebih baik. Namun, ENS_EQ tidak mempertimbangkan metode mana yang sebenarnya lebih baik. Itu yang dicoba diatasi oleh ENS_W.

5 Cross-Validation

Cross-validation menilai akurasi prediksi pada titik yang tidak dipakai untuk melatih model. Dengan begitu, hasilnya menggambarkan kinerja pada lokasi baru, bukan hanya seberapa baik model cocok dengan data latih.

5.1 Pembagian fold

set.seed(123)
K <- 5
n <- nrow(meuse_sf)
fold <- sample(rep(1:K, length.out = n))

table(fold)

## fold
##  1  2  3  4  5 
## 31 31 31 31 31

Data dibagi menjadi 5 bagian. sample(rep(1:K, ...)) mengacak keanggotaan tiap titik secara merata, sehingga jumlah titik per fold kira-kira sama (sekitar 31 titik per fold, karena 155/5 = 31). set.seed(123) memastikan pembagian yang sama setiap kali file dijalankan, sehingga hasil bisa direproduksi. Tabel di atas harus menunjukkan jumlah yang hampir sama di setiap fold. Jika tidak, ada yang salah dalam pembagian.

5.2 Prediksi out-of-fold

Untuk setiap fold $k$:

Data latih adalah semua fold selain $k$.
Data uji adalah fold $k$.
Model, termasuk variogram kriging, dilatih ulang hanya dengan data latih. Ini penting agar evaluasi tidak bias.

base_methods <- c("CO", "NN", "IDW", "KRIG")

cv_pred <- data.frame(
  id   = seq_len(n),
  obs  = meuse_sf$logzinc,
  fold = fold,
  CO   = NA_real_,
  NN   = NA_real_,
  IDW  = NA_real_,
  KRIG = NA_real_
)

for (k in 1:K) {
  tr_idx <- which(fold != k)
  te_idx <- which(fold == k)

  tr <- meuse_sf[tr_idx, ]
  te <- meuse_sf[te_idx, ]

  # 1. Closest Observation (Thiessen)
  nearest <- st_nearest_feature(te, tr)
  cv_pred$CO[te_idx] <- tr$logzinc[nearest]

  # 2. Nearest Neighbors (k = 5)
  g_nn <- gstat(formula = logzinc ~ 1, locations = tr,
                nmax = 5, set = list(idp = 0))
  cv_pred$NN[te_idx] <- predict(g_nn, newdata = te, debug.level = 0)$var1.pred

  # 3. IDW (p = 2)
  g_idw <- gstat(formula = logzinc ~ 1, locations = tr,
                 nmax = nrow(tr), set = list(idp = 2))
  cv_pred$IDW[te_idx] <- predict(g_idw, newdata = te, debug.level = 0)$var1.pred

  # 4. Kriging (variogram di-fit ulang pada data latih)
  v_tr  <- variogram(logzinc ~ 1, data = tr)
  fv_tr <- fit.variogram(v_tr, model = v_init)
  cv_pred$KRIG[te_idx] <- krige(logzinc ~ 1, locations = tr, newdata = te,
                                model = fv_tr, debug.level = 0)$var1.pred
}

head(cv_pred)

##   id      obs fold       CO       NN      IDW     KRIG
## 1  1 6.929517    4 7.039660 6.180599 6.544482 6.764230
## 2  2 7.039660    5 6.929517 6.389167 6.491090 6.808790
## 3  3 6.461468    3 6.929517 6.283336 6.340148 6.446181
## 4  4 5.549076    3 5.594711 6.283336 5.943268 6.020710
## 5  5 5.594711    5 5.638355 5.453975 5.721295 5.565778
## 6  6 5.638355    3 5.594711 5.504719 5.692852 5.499814

Loop ini menghasilkan satu prediksi untuk setiap titik, dari model yang tidak pernah melihat titik itu. Setiap titik hanya diprediksi satu kali, sehingga kolom CO, NN, IDW, dan KRIG berisi prediksi out-of-fold untuk seluruh 155 titik. Ada dua hal yang perlu diperhatikan:

Variogram kriging di-fit ulang di setiap fold. Jika variogram dipasang sekali pada seluruh data, informasi dari titik uji sudah ikut masuk ke model, dan RMSE kriging akan terlalu optimis.
head() hanya menunjukkan lima baris pertama. Pastikan tidak ada NA di kolom prediksi, karena NA berarti ada titik yang tidak diprediksi.

5.3 Ensemble pada hasil cross-validation

# ENS_EQ: bobot sama
cv_pred$ENS_EQ <- rowMeans(cv_pred[, base_methods])

# RMSE pooled tiap metode dasar, dipakai untuk bobot ENS_W
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))

rmse_base <- sapply(base_methods, function(m) rmse(cv_pred$obs, cv_pred[[m]]))
w_inv <- (1 / rmse_base) / sum(1 / rmse_base)

# ENS_W: bobot berbanding terbalik dengan RMSE
cv_pred$ENS_W <- as.matrix(cv_pred[, base_methods]) %*% w_inv

round(data.frame(RMSE = rmse_base, Bobot_ENS_W = w_inv), 4)

##        RMSE Bobot_ENS_W
## CO   0.6554      0.1940
## NN   0.4994      0.2547
## IDW  0.5237      0.2428
## KRIG 0.4122      0.3085

Fungsi rmse() menghitung akar rata-rata selisih kuadrat, sesuai rumus pada bagian evaluasi. Tabel memperlihatkan RMSE tiap metode dasar dan bobot ENS_W. Metode dengan RMSE lebih kecil mendapat bobot lebih besar, dan jumlah seluruh bobot harus tepat 1.

Bobot ENS_W diturunkan dari RMSE pada data yang juga dipakai untuk evaluasi. Akibatnya, RMSE ENS_W di tabel perbandingan cenderung terlalu optimis. Untuk penilaian yang lebih ketat, perlu nested cross-validation, yaitu bobot dihitung dari data latih di dalam setiap fold.

grid_sf$ENS_W <- as.vector(
  as.matrix(st_drop_geometry(grid_sf)[, base_methods]) %*% w_inv
)

Bobot yang dipakai sama dengan bobot cross-validation di atas, sehingga peta ENS_W konsisten dengan tabel RMSE.

6 Perbandingan RMSE

6.1 Fungsi evaluasi dan label metode

method_labels <- c(
  CO     = "Closest Observation",
  NN     = "Nearest Neighbors (k=5)",
  IDW    = "IDW (p=2)",
  KRIG   = "Ordinary Kriging",
  ENS_EQ = "Ensemble (bobot sama)",
  ENS_W  = "Ensemble (bobot 1/RMSE)"
)
all_methods <- names(method_labels)

RMSE dihitung dengan rumus

\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} \left(y_i - \hat{y}_i\right)^2} \]

Tiga ukuran dihasilkan:

RMSE pooled (skala log): semua prediksi out-of-fold digabung, lalu dihitung RMSE-nya.
RMSE rata-rata fold (skala log): RMSE tiap fold dihitung, lalu dirata-ratakan.
RMSE skala asli (ppm): prediksi ditransformasi balik dengan $\exp(\cdot)$, lalu dibandingkan dengan zinc asli.

Ketiga ukuran ini sebaiknya dibaca bersama. RMSE pooled adalah ukuran utama karena memakai seluruh 155 prediksi. RMSE rata-rata fold menunjukkan apakah hasil konsisten antar bagian data. Jika kedua ukuran memberi urutan metode yang berbeda, hal tersebut menandakan perbedaan antar metode tidak stabil.

Untuk RMSE dalam ppm, perlu hati-hati. $\exp(\hat{y})$ memberi median prediksi, bukan rata-rata, sehingga ada bias kecil. Karena itu ukuran ini sebaiknya dipakai sebagai pelengkap, bukan dasar utama perbandingan. RMSE log lebih jujur untuk membandingkan metode.

6.2 RMSE per fold

rmse_fold <- sapply(all_methods, function(m) {
  sapply(1:K, function(k) {
    i <- which(cv_pred$fold == k)
    rmse(cv_pred$obs[i], cv_pred[[m]][i])
  })
})
rownames(rmse_fold) <- paste("Fold", 1:K)
colnames(rmse_fold) <- method_labels[all_methods]

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

RMSE tiap fold (skala log)
	Closest Observation	Nearest Neighbors (k=5)	IDW (p=2)	Ordinary Kriging	Ensemble (bobot sama)	Ensemble (bobot 1/RMSE)
Fold 1	0.5710	0.5261	0.5616	0.4214	0.4800	0.4756
Fold 2	0.4728	0.3685	0.4216	0.2919	0.3447	0.3380
Fold 3	0.7532	0.5292	0.5505	0.4411	0.5217	0.5087
Fold 4	0.7641	0.5305	0.5362	0.4685	0.5313	0.5186
Fold 5	0.6685	0.5223	0.5362	0.4157	0.4778	0.4668

Setiap baris adalah satu fold, dan setiap kolom adalah satu metode. Metode yang selalu lebih kecil di hampir semua fold lebih stabil dan layak dipercaya. Metode yang menang di satu fold tetapi kalah di fold lain hasilnya sensitif terhadap titik mana yang masuk data uji. Fold dengan RMSE tinggi untuk semua metode biasanya berisi titik-titik ekstrem, misalnya nilai zinc yang sangat tinggi di area kecil yang sulit diprediksi.

Tabel ini membantu memastikan bahwa keunggulan satu metode bukan kebetulan pada satu pembagian data.

6.3 Ringkasan perbandingan

rmse_summary <- data.frame(
  Metode = unname(method_labels[all_methods]),
  RMSE_pooled_log = sapply(all_methods, function(m) rmse(cv_pred$obs, cv_pred[[m]])),
  RMSE_rata_fold_log = colMeans(rmse_fold),
  RMSE_skala_ppm = sapply(all_methods, function(m) {
    rmse(meuse$zinc, exp(cv_pred[[m]]))
  })
)
rmse_summary <- rmse_summary[order(rmse_summary$RMSE_pooled_log), ]
rownames(rmse_summary) <- NULL

kable(
  rmse_summary %>%
    mutate(across(where(is.numeric), ~ round(.x, 4))),
  caption = "Perbandingan RMSE metode interpolasi (cross-validation 5-fold)",
  col.names = c("Metode", "RMSE pooled (log)",
                "RMSE rata-rata fold (log)", "RMSE (ppm)")
)

Perbandingan RMSE metode interpolasi (cross-validation 5-fold)
Metode	RMSE pooled (log)	RMSE rata-rata fold (log)	RMSE (ppm)
Ordinary Kriging	0.4122	0.4077	229.9821
Ensemble (bobot 1/RMSE)	0.4661	0.4615	258.9753
Ensemble (bobot sama)	0.4758	0.4711	261.6825
Nearest Neighbors (k=5)	0.4994	0.4953	279.0649
IDW (p=2)	0.5237	0.5212	299.4244
Closest Observation	0.6554	0.6459	351.9517

Tabel diurutkan dari RMSE pooled terkecil ke terbesar, sehingga baris pertama adalah metode terbaik menurut ukuran utama. Metode terbaik pada data ini adalah Ordinary Kriging dengan RMSE pooled log sebesar 0.4122. Metode dengan RMSE terbesar adalah Closest Observation (0.6554).

Selisih antar metode perlu dinilai dalam konteks. Jika selisih RMSE kecil, misalnya di bawah sekitar 0,01 pada skala log, perbedaan itu belum tentu bermakna. Dalam kasus seperti itu, pilihan metode bisa didasarkan pada kemudahan atau kesesuaian dengan tujuan analisis. Ringkasan ini juga tidak menunjukkan ketidakpastian. Untuk itu, perlu melihat variasi antar fold di tabel sebelumnya.

6.4 Visualisasi RMSE

par(mar = c(9, 4, 3, 1))
bp <- barplot(rmse_summary$RMSE_pooled_log,
              names.arg = rmse_summary$Metode,
              las = 2, col = "#5EEAD4", border = "#0D9488",
              ylab = "RMSE (log)", main = "RMSE pooled cross-validation 5-fold",
              ylim = c(0, max(rmse_summary$RMSE_pooled_log) * 1.15))
text(bp, rmse_summary$RMSE_pooled_log,
     labels = round(rmse_summary$RMSE_pooled_log, 3), pos = 3, cex = 0.8)

par(mar = c(5, 4, 4, 2) + 0.1)

Grafik batang membuat perbedaan RMSE lebih mudah dilihat. Batang yang lebih pendek berarti galat lebih kecil. Perhatikan bahwa sumbu vertikal tidak dimulai dari nol secara ekstrem, sehingga selisih kecil bisa tampak besar. Angka di atas batang memberi nilai pastinya, sehingga grafik dan tabel sebaiknya dibaca bersama.

6.5 Observasi vs prediksi

cv_long <- cv_pred %>%
  select(obs, all_of(all_methods)) %>%
  pivot_longer(-obs, names_to = "Metode", values_to = "pred") %>%
  mutate(Metode = factor(method_labels[Metode], levels = method_labels))

ggplot(cv_long, aes(x = pred, y = obs)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#115E59") +
  geom_point(color = "#0D9488", alpha = 0.6, size = 1.6) +
  facet_wrap(~ Metode, ncol = 3) +
  labs(title = "Nilai observasi vs prediksi out-of-fold (skala log)",
       x = "Prediksi log(zinc)", y = "Observasi log(zinc)") +
  theme_minimal(base_size = 11) +
  theme(strip.text = element_text(color = "#115E59", face = "bold"))

Garis putus-putus adalah garis prediksi sempurna, yaitu observasi sama dengan prediksi. Makin dekat titik ke garis, makin akurat metode. Yang perlu diperhatikan:

Sebaran di sekitar garis: sebaran yang sempit berarti galat kecil dan konsisten.
Pola sistematis: jika titik cenderung berada di atas garis di prediksi rendah dan di bawah garis di prediksi tinggi, metode cenderung terlalu “menarik” nilai ke rata-rata.
Titik ekstrem: titik yang sangat jauh dari garis adalah sumber utama RMSE besar.

Metode yang sebarannya paling dekat dengan garis dan paling sedikit pola sistematis biasanya adalah yang terbaik.

7 Peta Prediksi

Peta dibuat dengan model yang dilatih pada seluruh data. Setiap sel grid diberi warna sesuai nilai prediksi log(zinc).

grid_map <- grid_sf %>%
  st_drop_geometry() %>%
  select(x, y, all_of(all_methods)) %>%
  pivot_longer(all_of(all_methods), names_to = "Metode", values_to = "pred") %>%
  mutate(Metode = factor(method_labels[Metode], levels = method_labels))

ggplot(grid_map, aes(x = x, y = y, fill = pred)) +
  geom_tile(width = 40, height = 40) +
  geom_point(data = meuse, aes(x = x, y = y), inherit.aes = FALSE,
             size = 0.6, color = "black", alpha = 0.6) +
  facet_wrap(~ Metode, ncol = 2) +
  scale_fill_viridis_c(name = "log(zinc)", option = "viridis") +
  coord_equal() +
  labs(title = "Peta prediksi log(zinc) pada grid meuse",
       x = "Easting", y = "Northing") +
  theme_minimal(base_size = 11) +
  theme(strip.text = element_text(color = "#115E59", face = "bold"),
        legend.position = "right")

Setiap panel adalah satu metode dengan skala warna yang sama sehingga bisa dibandingkan langsung. Titik hitam adalah lokasi sampel. Hal yang perlu dibandingkan:

Closest Observation: tampak seperti mosaik dengan batas poligon yang tegas.
NN dan IDW: lebih halus. IDW biasanya menghasilkan “bintik” di sekitar titik sampel.
Kriging: memperlihatkan pola yang lebih menyerupai kontur, dan pola ini konsisten dengan struktur spasial di variogram.
Ensemble: merupakan campuran metode, sehingga bentuknya di antara komponennya.

Peta menunjukkan pola, tetapi tidak menunjukkan akurasi. Peta yang paling halus atau paling indah belum tentu yang paling akurat. Akurasi dinilai dari RMSE. Peta hanya membantu menjelaskan mengapa suatu metode bekerja seperti itu.

8 Kesimpulan

best  <- rmse_summary[1, ]
worst <- rmse_summary[nrow(rmse_summary), ]
cat(sprintf(
  "- Berdasarkan RMSE pooled (skala log), metode terbaik adalah **%s** (%.4f), dan terburuk adalah **%s** (%.4f).\n",
  best$Metode, best$RMSE_pooled_log, worst$Metode, worst$RMSE_pooled_log))

Berdasarkan RMSE pooled (skala log), metode terbaik adalah Ordinary Kriging (0.4122), dan terburuk adalah Closest Observation (0.6554).

cat(sprintf(
  "- Selisih terbesar antara terbaik dan terburuk adalah %.4f pada skala log, yaitu sekitar %.1f%% dari RMSE terbaik.\n",
  worst$RMSE_pooled_log - best$RMSE_pooled_log,
  100 * (worst$RMSE_pooled_log - best$RMSE_pooled_log) / best$RMSE_pooled_log))

Selisih terbesar antara terbaik dan terburuk adalah 0.2432 pada skala log, yaitu sekitar 59.0% dari RMSE terbaik.

cat("- Kriging memakai struktur spasial dari variogram, sehingga bobotnya mempertimbangkan korelasi antar titik, bukan hanya jarak.\n")

Kriging memakai struktur spasial dari variogram, sehingga bobotnya mempertimbangkan korelasi antar titik, bukan hanya jarak.

cat("- Thiessen (Closest Observation) paling sederhana, tetapi prediksinya melompat di batas poligon. Kelemahan ini biasanya membuat galatnya lebih besar.\n")

Thiessen (Closest Observation) paling sederhana, tetapi prediksinya melompat di batas poligon. Kelemahan ini biasanya membuat galatnya lebih besar.

cat("- Ensemble berbobot 1/RMSE sedikit optimistis karena bobotnya diturunkan dari data evaluasi yang sama. Hasilnya perlu dibaca dengan hati-hati.\n")

Ensemble berbobot 1/RMSE sedikit optimistis karena bobotnya diturunkan dari data evaluasi yang sama. Hasilnya perlu dibaca dengan hati-hati.

Pemilihan metode interpolasi sebaiknya didasarkan pada RMSE cross-validation, bukan pada tampilan peta. Peta yang halus belum tentu akurat, sehingga evaluasi numerik tetap diperlukan. Keterbatasan laporan ini adalah belum ada nested cross-validation untuk ENS_W, dan belum ada uji sensitivitas untuk parameter $p$ pada IDW serta nilai $k$ pada NN. Meskipun ENS_W sedikit optimistis, nilai RMSE-nya tetap lebih tinggi daripada kriging, sehingga kesimpulan bahwa kriging adalah metode terbaik pada data ini tidak berubah.

Tugas Praktikum 10 Regresi Spasial - Interpolasi Data Meuse

Arifatu Dzikra