TUGAS PRAKTIKUM 2

1. Library

# install.packages(c("sp", "sf", "gstat", "ggplot2", "dplyr", "tidyr"))
library(sp)
library(sf)
library(gstat)
library(ggplot2)
library(dplyr)
library(tidyr)

2. Memuat Data

Variabel yang diinterpolasi: zinc (konsentrasi seng pada tanah).

data(meuse)       # 155 titik observasi
data(meuse.grid)  # grid lokasi yang akan diprediksi

# Konversi ke objek sf (CRS: Amersfoort / RD New = EPSG 28992)
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, remove = FALSE)
grid_sf  <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992, remove = FALSE)

head(meuse_sf)

## Simple feature collection with 6 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 181025 ymin: 333260 xmax: 181390 ymax: 333611
## Projected CRS: Amersfoort / RD New
##        x      y cadmium copper lead zinc  elev       dist   om ffreq soil lime
## 1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1
## 2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1
## 3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1
## 4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0
## 5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0
## 6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0
##   landuse dist.m              geometry
## 1      Ah     50 POINT (181072 333611)
## 2      Ah     30 POINT (181025 333558)
## 3      Ah    150 POINT (181165 333537)
## 4      Ga    270 POINT (181298 333484)
## 5      Ah    380 POINT (181307 333330)
## 6      Ga    470 POINT (181390 333260)

summary(meuse$zinc)

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

3. Eksplorasi Data

ggplot(meuse_sf) +
  geom_sf(aes(color = zinc), size = 2.5) +
  scale_color_viridis_c(option = "C") +
  labs(title = "Sebaran Observasi Zinc (Meuse)", color = "Zinc (ppm)") +
  theme_minimal()

hist(meuse$zinc, breaks = 20, col = "steelblue",
     main = "Histogram Zinc", xlab = "Zinc (ppm)")

Dari grafik di atas, dapat diketahui bahwa data berbentuk skewed ke kanan. - Sebagian besar lokasi punya kadar zinc relatif rendah, kira-kira di bawah 500 ppm. - Sebagian kecil lokasi punya kadar sangat tinggi (di atas 1.200 ppm), yang kemungkinan besar adalah titik-titik yang dekat dengan sungai atau sumber pencemar. Frekuensinya hanya 1-3 titik per kelas, jadi bisa dianggap pencilan atau nilai ekstrem. - Rata-rata (mean) akan lebih besar dari median karena tertarik oleh nilai-nilai ekstrem di ekor kanan.

4. Fungsi Prediksi untuk Setiap Metode

Satu fungsi dipakai untuk semua metode agar bisa dipakai untuk peta dan cross validation.

Closest Observation: nilai prediksi = nilai titik terdekat (1 tetangga).
NN (Nearest Neighbour): rata-rata k tetangga terdekat (bobot sama, idp = 0).
IDW: Inverse Distance Weighting dengan power idp = 2.
Kriging: Ordinary Kriging dengan model variogram hasil fitting.

# Variogram empiris
v_emp <- variogram(zinc ~ 1, meuse_sf)

# Fitting model variogram (spherical)
sill_awal <- max(v_emp$gamma)
v_fit <- fit.variogram(
  v_emp,
  model = vgm(psill = 0.7 * sill_awal, model = "Sph",
              range = 900, nugget = 0.3 * sill_awal)
)
v_fit

##   model     psill    range
## 1   Nug  24816.48   0.0000
## 2   Sph 134756.29 831.3829

plot(v_emp, model = v_fit, main = "Variogram Zinc")

Variogram zinc menunjukkan adanya autokorelasi spasial yang jelas: semivariansi naik seiring bertambahnya jarak (titik yang berdekatan memiliki nilai zinc yang mirip, sedangkan titik yang berjauhan semakin berbeda), lalu mendatar pada sekitar 160.000 (sill) di jarak kira-kira 800 m (range), sehingga di luar jarak tersebut observasi dianggap tidak lagi saling berkorelasi. Model spherical yang di-fit mengikuti titik-titik empiris dengan cukup baik, dengan nilai nugget sekitar 25.000 yang menandakan adanya variasi pada skala sangat kecil atau galat pengukuran. Efeknya pada interpolasi: struktur spasial yang kuat ini membuat Kriging layak digunakan karena bobot tetangganya ditentukan oleh pola variogram, bukan sekadar jarak seperti pada IDW, sehingga Kriging cenderung menghasilkan RMSE yang lebih kecil. Hanya titik dalam radius sekitar 800 m yang efektif berkontribusi pada prediksi, dan nugget yang tidak nol membuat hasil Kriging sedikit lebih halus (smooth) dibanding metode yang sekadar meniru nilai terdekat seperti Closest Observation.

k_nn  <- 5   # jumlah tetangga untuk NN
p_idw <- 2   # power IDW

pred_dasar <- function(train, test, metode) {
  switch(metode,
    closest = idw(zinc ~ 1, train, test, nmax = 1,
                  idp = 1, debug.level = 0)$var1.pred,
    nn      = idw(zinc ~ 1, train, test, nmax = k_nn,
                  idp = 0, debug.level = 0)$var1.pred,
    idw     = idw(zinc ~ 1, train, test,
                  idp = p_idw, debug.level = 0)$var1.pred,
    kriging = krige(zinc ~ 1, train, test,
                    model = v_fit, debug.level = 0)$var1.pred
  )
}

# Ensemble: rata-rata dari 3 metode (NN, IDW, Kriging)
metode_ensemble <- c("nn", "idw", "kriging")

pred_ensemble <- function(train, test) {
  hasil <- sapply(metode_ensemble, function(m) pred_dasar(train, test, m))
  if (is.null(dim(hasil))) hasil <- matrix(hasil, nrow = 1)
  rowMeans(hasil)
}

5. Interpolasi ke Seluruh Grid

hasil_grid <- data.frame(
  x = meuse.grid$x,
  y = meuse.grid$y,
  Closest_Observation = pred_dasar(meuse_sf, grid_sf, "closest"),
  NN                  = pred_dasar(meuse_sf, grid_sf, "nn"),
  IDW                 = pred_dasar(meuse_sf, grid_sf, "idw"),
  Kriging             = pred_dasar(meuse_sf, grid_sf, "kriging")
)

# Ensemble = rata-rata NN, IDW, dan Kriging
hasil_grid$Ensemble <- rowMeans(hasil_grid[, c("NN", "IDW", "Kriging")])

head(hasil_grid)

##        x      y Closest_Observation    NN      IDW  Kriging Ensemble
## 1 181180 333740                1022 681.2 633.6864 762.1218 692.3361
## 2 181140 333700                1022 681.2 712.5450 823.7930 739.1793
## 3 181180 333700                1022 681.2 654.1617 758.4210 697.9276
## 4 181220 333700                 640 681.2 604.4422 695.7806 660.4743
## 5 181100 333660                1022 681.2 857.2558 899.3148 812.5902
## 6 181140 333660                1022 681.2 755.5061 822.1090 752.9384

6. Visualisasi Peta Hasil Interpolasi

peta_long <- hasil_grid %>%
  pivot_longer(cols = -c(x, y), names_to = "Metode", values_to = "Zinc") %>%
  mutate(Metode = factor(Metode,
         levels = c("Closest_Observation", "NN", "IDW", "Ensemble", "Kriging")))

ggplot(peta_long, aes(x = x, y = y, fill = Zinc)) +
  geom_tile() +
  coord_equal() +
  facet_wrap(~ Metode, ncol = 3) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Hasil Interpolasi Zinc Berdasarkan Metode",
       fill = "Zinc (ppm)", x = "X", y = "Y") +
  theme_minimal()

7. Cross Validation (Leave-One-Out)

Setiap titik dikeluarkan satu per satu, diprediksi menggunakan titik lainnya, lalu dibandingkan dengan nilai aslinya.

n <- nrow(meuse_sf)
metode_dasar <- c("closest", "nn", "idw", "kriging")

cv_pred <- matrix(NA, nrow = n, ncol = length(metode_dasar),
                  dimnames = list(NULL, metode_dasar))

for (i in 1:n) {
  train <- meuse_sf[-i, ]
  test  <- meuse_sf[i, ]
  for (m in metode_dasar) {
    cv_pred[i, m] <- pred_dasar(train, test, m)
  }
}

cv_df <- as.data.frame(cv_pred)
cv_df$ensemble <- rowMeans(cv_df[, metode_ensemble])  # NN, IDW, Kriging
cv_df$observasi <- meuse_sf$zinc

head(cv_df)

##   closest    nn      idw  kriging ensemble observasi
## 1    1141 558.0 793.8598 890.6618 747.5072      1022
## 2    1022 702.0 727.3774 868.8288 766.0687      1141
## 3    1022 634.4 605.2803 661.2713 633.6505       640
## 4     640 511.6 477.0037 465.8510 484.8182       257
## 5     281 269.2 368.8423 268.9481 302.3301       269
## 6     269 248.8 347.1991 235.1268 277.0420       281

8. Perhitungan dan Perbandingan RMSE

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

tabel_rmse <- data.frame(
  Metode = c("Closest Observation", "NN", "IDW", "Ensemble", "Kriging"),
  RMSE   = c(rmse(cv_df$observasi, cv_df$closest),
             rmse(cv_df$observasi, cv_df$nn),
             rmse(cv_df$observasi, cv_df$idw),
             rmse(cv_df$observasi, cv_df$ensemble),
             rmse(cv_df$observasi, cv_df$kriging))
) %>%
  arrange(RMSE)

tabel_rmse

##                Metode     RMSE
## 1             Kriging 224.7995
## 2            Ensemble 247.1178
## 3                  NN 263.9874
## 4                 IDW 278.2734
## 5 Closest Observation 309.9685

Berdasarkan niali RMSE yang didapatkan, Kriging merupakan metode terbaik (RMSE = 224,80), diikuti Ensemble (247,12), NN (263,99), IDW (278,27), dan Closest Observation (309,97). Hal ini menunjukkan bahwa pemanfaatan struktur autokorelasi spasial melalui variogram memberikan prediksi yang lebih akurat dibandingkan metode yang hanya mempertimbangkan jarak.

ggplot(tabel_rmse, aes(x = reorder(Metode, RMSE), y = RMSE, fill = Metode)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 2)), vjust = -0.4) +
  labs(title = "Perbandingan RMSE (Leave-One-Out Cross Validation)",
       x = "Metode", y = "RMSE") +
  theme_minimal()

cv_df %>%
  pivot_longer(cols = c(closest, nn, idw, ensemble, kriging),
               names_to = "Metode", values_to = "Prediksi") %>%
  ggplot(aes(x = observasi, y = Prediksi)) +
  geom_point(alpha = 0.6) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  facet_wrap(~ Metode) +
  labs(title = "Observasi vs Prediksi (LOOCV)",
       x = "Observasi", y = "Prediksi") +
  theme_minimal()

9. Kesimpulan

terbaik <- tabel_rmse$Metode[1]
cat("Metode dengan RMSE terkecil:", terbaik,
    "dengan RMSE =", round(tabel_rmse$RMSE[1], 2), "\n")

## Metode dengan RMSE terkecil: Kriging dengan RMSE = 224.8

Berdasarkan hasil Leave-One-Out Cross Validation, Kriging merupakan metode interpolasi terbaik untuk data zinc pada meuse dengan RMSE terkecil (224,80), diikuti Ensemble (247,12), NN (263,99), IDW (278,27), dan Closest Observation (309,97). Keunggulan Kriging disebabkan oleh adanya autokorelasi spasial yang kuat pada variogram, yang dimanfaatkan untuk menentukan bobot tetangga secara lebih tepat dibanding metode yang hanya mempertimbangkan jarak. Ensemble mampu mengungguli NN dan IDW, tetapi tidak mengungguli Kriging, sehingga penggabungan metode tidak memberi perbaikan ketika salah satu komponennya jauh lebih akurat. Dengan demikian, Kriging direkomendasikan untuk menginterpolasi zinc pada data ini.