1 Pendahuluan

1.1 Tentang Studi Kasus

Dokumen ini menerapkan empat varian Kriging yang dibahas pada Pertemuan 6 — Ordinary, Simple, Universal, dan Co-Kriging — lengkap dengan variogram, peta prediksi, peta ketidakpastian, dan perbandingan akurasi melalui cross-validation, pada studi kasus curah hujan tahunan di Provinsi Banten.

Dataset yang digunakan adalah data simulasi berisi 35 titik pemantauan yang tersebar acak di seluruh Provinsi Banten, dengan dua variabel:

curah_hujan — curah hujan tahunan (mm/tahun), variabel utama yang ingin diinterpolasi;
elevasi_m — elevasi tempat (m), variabel sekunder yang dipakai pada Co-Kriging.

Peta batas wilayah memakai shapefile resmi BPS (8 kabupaten/kota).

1.2 Rumusan Pertanyaan Analisis

Apakah curah hujan memiliki autokorelasi spasial (titik berdekatan cenderung mirip) sehingga layak diinterpolasi dengan Kriging?
Model variogram teoritis apa yang cocok, dan berapa nugget, sill, dan range-nya?
Bagaimana peta curah hujan Banten hasil Ordinary Kriging, beserta ketidakpastiannya?
Bagaimana hasil Simple, Universal, dan Co-Kriging, dan apakah berbeda secara berarti dibandingkan Ordinary Kriging?
Metode mana yang paling akurat menurut leave-one-out cross-validation?

2 Persiapan

2.1 Memuat Package

library(sf)        # data spasial: baca SHP, proyeksi, titik, grid
library(gstat)     # geostatistik: variogram dan kriging
library(dplyr)     # manipulasi data
library(ggplot2)   # visualisasi
library(readxl)    # membaca berkas Excel
library(viridis)   # palet warna ramah buta warna

2.2 Langkah 1 — Import Data dan Peta

setwd("C:/Users/mifta/OneDrive - untirta.ac.id/UNTIRTA/STATISTIKA/Materi MK Semester Ganjil/Statistika Lingkungan/6-Studi Kasus Kriging")
# --- Peta batas wilayah (8 kabupaten/kota) -----------------------------------
batas <- st_read("2-SHPBANTEN/Banten_ADMIN_BPS.shp", quiet = TRUE)
print(batas[, c("Kabupaten", "PROVINCE")])

## Simple feature collection with 8 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -7.016779 xmax: 106.7799 ymax: -5.807418
## Geodetic CRS:  WGS 84
##                Kabupaten PROVINCE                       geometry
## 1           Kota Cilegon   Banten MULTIPOLYGON (((105.928 -6....
## 2            Kota Serang   Banten MULTIPOLYGON (((106.1957 -6...
## 3         Kota Tangerang   Banten MULTIPOLYGON (((106.6872 -6...
## 4 Kota Tangerang Selatan   Banten MULTIPOLYGON (((106.6463 -6...
## 5                  Lebak   Banten MULTIPOLYGON (((106.2797 -6...
## 6             Pandeglang   Banten MULTIPOLYGON (((105.5498 -6...
## 7                 Serang   Banten MULTIPOLYGON (((105.8577 -6...
## 8              Tangerang   Banten MULTIPOLYGON (((106.6442 -6...

cat("CRS asli peta:", st_crs(batas)$input, "(koordinat dalam derajat)\n")

## CRS asli peta: WGS 84 (koordinat dalam derajat)

# --- Data 35 titik pengamatan -------------------------------------------------
berkas_data <- "6-Data Simulasi - Kriging Curah Hujan & Elevasi 35 Titik Banten.xlsx"
dat <- read_excel(berkas_data, sheet = "Data Kriging")
str(dat)

## tibble [35 × 5] (S3: tbl_df/tbl/data.frame)
##  $ id_titik   : chr [1:35] "T01" "T02" "T03" "T04" ...
##  $ lon        : num [1:35] 107 106 106 106 106 ...
##  $ lat        : num [1:35] -6.3 -6.8 -6.78 -6.3 -6.57 ...
##  $ curah_hujan: num [1:35] 2196 2211 2108 2236 2193 ...
##  $ elevasi_m  : num [1:35] 148 305 297 202 201 7 5 325 82 328 ...

head(dat)

Pemeriksaan kualitas data sebelum analisis spasial:

cat("Jumlah titik       :", nrow(dat), "\n")

## Jumlah titik       : 35

cat("Nilai hilang (NA)  :", sum(is.na(dat)), "\n")

## Nilai hilang (NA)  : 0

cat("Koordinat duplikat :", sum(duplicated(dat[, c("lon", "lat")])), "\n")

## Koordinat duplikat : 0

Interpretasi

Data terdiri dari 35 titik tanpa nilai hilang dan tanpa koordinat ganda — syarat dasar yang harus terpenuhi karena Kriging tidak dapat menangani dua titik pada lokasi yang persis sama (matriks kovariansi menjadi singular).

2.3 Langkah 2 — Proyeksi Koordinat

Variogram dihitung dari jarak antartitik, sedangkan lon/lat (derajat) bukan satuan jarak yang seragam. Karena itu data diproyeksikan ke UTM zona 48S (EPSG:32748) agar jarak dinyatakan dalam meter. Parameter range (mis. 45 km) baru bermakna bila satuan jarak konsisten.

UTM48S <- 32748
km <- 1000                                   # 1 km = 1000 m

batas_utm  <- st_transform(batas, UTM48S)
batas_prov <- st_union(batas_utm)            # batas luar provinsi (gabungan 8 kab/kota)

titik <- st_as_sf(dat, coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>%
  st_transform(UTM48S)

# Simpan koordinat X, Y (meter) sebagai kolom biasa -> dibutuhkan Universal Kriging
xy <- st_coordinates(titik)
titik$x <- xy[, 1]
titik$y <- xy[, 2]

cat("Jarak maksimum antar titik:", round(max(st_distance(titik)) / km, 1), "km\n")

## Jarak maksimum antar titik: 147.6 km

Catatan. Pada slide, data dibiarkan dalam derajat (EPSG:4326). Untuk koordinat geografis, gstat otomatis menghitung jarak great-circle dalam km, sehingga range = 45 pada slide sama maknanya dengan range = 45 * km pada dokumen ini. Kedua cara menghasilkan RMSE yang sama (dibuktikan pada Langkah 12).

3 Eksplorasi Data

3.1 Langkah 3 — Statistik Deskriptif dan Korelasi

ringkas <- data.frame(
  Variabel = c("Curah hujan (mm/tahun)", "Elevasi (m)"),
  n      = nrow(dat),
  Min    = c(min(dat$curah_hujan),    min(dat$elevasi_m)),
  Rerata = c(mean(dat$curah_hujan),   mean(dat$elevasi_m)),
  Median = c(median(dat$curah_hujan), median(dat$elevasi_m)),
  Maks   = c(max(dat$curah_hujan),    max(dat$elevasi_m)),
  SD     = c(sd(dat$curah_hujan),     sd(dat$elevasi_m))
)
ringkas[, -1] <- round(ringkas[, -1], 1)
ringkas

uji_r <- cor.test(dat$curah_hujan, dat$elevasi_m)
uji_r

## 
##  Pearson's product-moment correlation
## 
## data:  dat$curah_hujan and dat$elevasi_m
## t = 3.1428, df = 33, p-value = 0.003525
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1746474 0.7010733
## sample estimates:
##     cor 
## 0.47996

ggplot() +
  geom_sf(data = batas_utm, fill = "grey95", colour = "grey50") +
  geom_sf(data = titik, aes(colour = curah_hujan, size = elevasi_m)) +
  scale_colour_viridis_c(option = "turbo", name = "Curah hujan\n(mm/tahun)") +
  scale_size_continuous(name = "Elevasi (m)", range = c(1.5, 5)) +
  labs(title = "Sebaran 35 Titik Pengamatan - Provinsi Banten",
       subtitle = "Warna = curah hujan, ukuran titik = elevasi") +
  theme_minimal()

Interpretasi

Curah hujan tahunan berkisar 1942–2440 mm (rerata 2178.1 mm, SD 108.7 mm), sedangkan elevasi berkisar 5–485 m. Korelasi curah hujan–elevasi r = 0.48 (p = 0.0035): positif dan signifikan, namun hanya sedang — elevasi menjelaskan sekitar 23% variasi curah hujan. Inilah alasan Co-Kriging layak dicoba, tetapi sekaligus petunjuk bahwa perbaikan yang diperoleh tidak akan dramatis.

4 Variogram

4.1 Langkah 4 — Variogram Empiris

Variogram empiris mengukur seberapa berbeda dua titik sebagai fungsi jaraknya:

\[\hat{\gamma}(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}\left[z(s_i)-z(s_i+h)\right]^2\]

Bila ada autokorelasi spasial, \(\hat{\gamma}(h)\) kecil untuk jarak dekat lalu naik seiring jarak hingga mendatar (sill).

vgm_awan <- variogram(curah_hujan ~ 1, data = titik, cloud = TRUE)
plot(vgm_awan, main = "Variogram Awan (setiap pasangan titik)",
     xlab = "Jarak (m)", ylab = "Semivarians")

Variogram awan sangat tersebar sehingga sulit dibaca. Pasangan titik karenanya dirata-ratakan per kelas jarak (lebar 10 km, hingga 80 km):

vgm_emp <- variogram(curah_hujan ~ 1, data = titik,
                     width = 10 * km, cutoff = 80 * km)
vgm_emp

ggplot(vgm_emp, aes(dist / km, gamma)) +
  geom_point(size = 3, colour = "#00207E") +
  geom_text(aes(label = np), vjust = -1, size = 3) +
  labs(title = "Variogram Empiris Curah Hujan",
       subtitle = "Angka di atas titik = jumlah pasangan titik (np) pada kelas jarak tersebut",
       x = "Jarak (km)", y = "Semivarians") +
  ylim(0, max(vgm_emp$gamma) * 1.1) +
  theme_minimal()

Catatan teknis. width dan cutoff sengaja diatur. Bila dibiarkan default, gstat memakai 15 kelas jarak yang sempit dengan hanya 3–30 pasangan per kelas, sehingga kurva tampak sangat bergerigi. Aturan praktis: cutoff ≤ separuh jarak terjauh, dan tiap kelas sebaiknya memuat ≥ 30 pasangan.

Interpretasi

Semivarians naik dari 6981 pada jarak ≈ 7.1 km menjadi 14395 pada jarak ≈ 74.6 km. Pola “dekat = mirip, jauh = berbeda” ini adalah bukti autokorelasi spasial positif (Hukum Tobler), sehingga interpolasi Kriging layak dilakukan. Nilai pada jarak terdekat yang sudah mencapai ≈ 6981 — bukan mendekati nol — menandakan adanya efek nugget yang besar: variasi acak/galat ukur yang tidak dapat dijelaskan oleh jarak.

4.2 Langkah 5 — Model Variogram Teoritis

Variogram empiris hanyalah titik-titik; Kriging membutuhkan fungsi mulus yang valid. Parameter model:

Parameter	Arti
Nugget (\(c_0\))	loncatan pada \(h=0\): galat ukur + variasi skala sangat kecil
Partial sill (\(c_1\))	tinggi kenaikan di atas nugget; sill = \(c_0 + c_1\)
Range (\(a\))	jarak ketika korelasi spasial mulai hilang

Model Exponential yang dipakai (konvensi gstat):

\[\gamma(h)=c_0+c_1\left(1-e^{-h/a}\right)\]

4.2.1 Mengapa `fit.variogram()` otomatis perlu diwaspadai

awal_exp <- vgm(psill = 6500, model = "Exp", range = 45 * km, nugget = 7000)
fit_bebas <- fit.variogram(vgm_emp, model = awal_exp)   # semua parameter bebas

## Warning in fit.variogram(vgm_emp, model = awal_exp): No convergence after 200
## iterations: try different initial values?

fit_bebas

Peringatan No convergence dan range hasil fitting sebesar 1498 km — jauh melampaui ukuran Provinsi Banten (± 150 km) — menunjukkan optimasi gagal. Penyebabnya: variogram empiris belum mendatar dalam jarak 80 km, sehingga data tidak cukup “memberi tahu” di mana sill berada. Ini kejadian umum pada data berukuran kecil (n = 35).

4.2.2 Solusi praktis: tetapkan range secara visual

Pada praktiknya range ditentukan dari pembacaan variogram empiris (fit.ranges = FALSE), lalu nugget dan sill dioptimalkan. Tiga model kandidat dibandingkan:

fit_sph <- fit.variogram(vgm_emp, vgm(6500, "Sph", 60 * km, 7000), fit.ranges = FALSE)
fit_exp <- fit.variogram(vgm_emp, vgm(6500, "Exp", 45 * km, 7000), fit.ranges = FALSE)
fit_gau <- fit.variogram(vgm_emp, vgm(6500, "Gau", 40 * km, 7000), fit.ranges = FALSE)

# Model "kuliah": seluruh parameter ditetapkan dari pembacaan visual (sesuai slide)
model_fit <- vgm(psill = 6500, model = "Exp", range = 45 * km, nugget = 7000)

ssel <- function(f) attr(f, "SSErr")
kandidat <- list("Spherical (range 60 km, fit nugget & sill)"   = fit_sph,
                 "Exponential (range 45 km, fit nugget & sill)" = fit_exp,
                 "Gaussian (range 40 km, fit nugget & sill)"    = fit_gau,
                 "Exponential (ditetapkan visual: 7000/6500/45)" = model_fit)

# Pembanding yang adil: LOOCV Ordinary Kriging untuk tiap model
rmse_cv <- function(model) {
  cv <- krige.cv(curah_hujan ~ 1, titik, model = model, nfold = nrow(titik),
                 debug.level = 0)
  sqrt(mean(cv$residual^2))
}
tab_model <- data.frame(
  Model        = names(kandidat),
  Nugget       = sapply(kandidat, function(m) m$psill[1]),
  Partial_sill = sapply(kandidat, function(m) m$psill[2]),
  Range_km     = sapply(kandidat, function(m) m$range[2] / km),
  SSErr        = sapply(kandidat, function(m) if (is.null(attr(m, "SSErr"))) NA else ssel(m)),
  RMSE_LOOCV   = sapply(kandidat, rmse_cv),
  row.names = NULL
)
tab_model[, 2:3] <- round(tab_model[, 2:3])
tab_model[, 5] <- round(tab_model[, 5], 3)
tab_model[, 6] <- round(tab_model[, 6], 2)
tab_model

kurva <- bind_rows(lapply(names(kandidat)[c(1, 2, 3)], function(nm)
  variogramLine(kandidat[[nm]], maxdist = 80 * km) %>%
    mutate(Model = sub(" \\(.*", "", nm))))
kurva_manual <- variogramLine(model_fit, maxdist = 80 * km) %>%
  mutate(Model = "Exponential (ditetapkan)")

ggplot() +
  geom_line(data = kurva, aes(dist / km, gamma, colour = Model), linewidth = 1) +
  geom_line(data = kurva_manual, aes(dist / km, gamma, colour = Model),
            linewidth = 1.2, linetype = "dashed") +
  geom_point(data = vgm_emp, aes(dist / km, gamma), size = 3) +
  scale_colour_manual(values = c(Spherical = "#00207E", Exponential = "#FFA200",
                                 Gaussian = "#147D64",
                                 "Exponential (ditetapkan)" = "#C0392B")) +
  labs(title = "Variogram Empiris dan Model Teoritis",
       x = "Jarak (km)", y = "Semivarians") +
  theme_minimal()

nugget_hat <- model_fit$psill[1]
psill_hat  <- model_fit$psill[2]
range_hat  <- model_fit$range[2] / km
sill_hat   <- nugget_hat + psill_hat
cat(sprintf("Nugget = %.0f | Partial sill = %.0f | Sill = %.0f | Range = %.0f km\n",
            nugget_hat, psill_hat, sill_hat, range_hat))

## Nugget = 7000 | Partial sill = 6500 | Sill = 13500 | Range = 45 km

cat(sprintf("Rasio nugget/sill = %.1f%%  (<25%% kuat; 25-75%% sedang; >75%% lemah)\n",
            100 * nugget_hat / sill_hat))

## Rasio nugget/sill = 51.9%  (<25% kuat; 25-75% sedang; >75% lemah)

Interpretasi

Model yang dipakai pada seluruh metode berikutnya adalah Exponential dengan nugget = 7000, partial sill = 6500 (sill = 13500), dan range = 45 km, sesuai slide materi.
Rasio nugget/sill = 52% → ketergantungan spasial tergolong sedang (batas 25–75%). Sekitar separuh variasi curah hujan bersifat “acak” pada skala jarak antartitik; sisanya terstruktur secara spasial.
Pada jarak 65–75 km seluruh kurva model berada di bawah dua titik empiris terjauh (± 14184 dan 14395). Hal ini dapat terkait tren regional lemah yang belum dimodelkan (Langkah 9) — tren membuat semivarians terus naik pada jarak jauh — namun dugaan ini bukan kepastian.
Jangan memilih model hanya dari SSErr: selisih antarmodel kecil, hanya 8 titik empiris yang dicocokkan, dan SSErr terkecil tidak menjamin prediksi terbaik. Kriteria yang lebih bermakna adalah RMSE cross-validation: keempat model memberi RMSE antara 94.67 dan 96.59 mm — selisih yang kecil (< 3%). Artinya hasil Kriging tidak sensitif terhadap pilihan model pada data ini. Model yang ditetapkan secara visual memang kebetulan memberi RMSE terkecil, tetapi selisihnya terlalu kecil untuk dijadikan bukti bahwa model itu “paling benar”; ia cukup sebagai model kerja yang memadai.

5 Prediksi Spasial

5.1 Langkah 6 — Grid Prediksi

Grid titik-titik di dalam Banten tempat curah hujan akan diduga (80 × 80 sel, hanya yang berada di dalam batas provinsi).

grid <- st_make_grid(batas_prov, n = c(80, 80), what = "centers") %>%
  st_sf(geometry = .) %>%
  st_filter(batas_prov)                      # buang titik di luar provinsi

xyg <- st_coordinates(grid)
grid$x <- xyg[, 1]
grid$y <- xyg[, 2]
cat("Jumlah sel grid dalam batas Banten:", nrow(grid), "\n")

## Jumlah sel grid dalam batas Banten: 2411

Fungsi bantu untuk memetakan hasil (dipakai berulang):

peta_pred <- function(hasil, judul, kolom = "var1.pred",
                      legenda = "Curah hujan\n(mm/tahun)", opsi = "turbo") {
  ggplot() +
    geom_tile(data = cbind(st_drop_geometry(hasil), st_coordinates(hasil)),
              aes(X, Y, fill = .data[[kolom]])) +
    geom_sf(data = batas_utm, fill = NA, colour = "grey20", linewidth = 0.3) +
    geom_sf(data = titik, shape = 21, fill = "white", size = 1.6) +
    scale_fill_viridis_c(option = opsi, name = legenda) +
    labs(title = judul, x = NULL, y = NULL) +
    theme_minimal()
}

Satu titik uji di tengah Banten (106,15°BT; 6,55°LS; elevasi diasumsikan 250 m) dipakai untuk membandingkan keluaran tiap metode secara tepat:

titik_uji <- st_as_sf(data.frame(lon = 106.15, lat = -6.55),
                      coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(UTM48S)
xyu <- st_coordinates(titik_uji)
titik_uji$x <- xyu[, 1]; titik_uji$y <- xyu[, 2]; titik_uji$elevasi_m <- 250

5.2 Langkah 7 — Ordinary Kriging (OK)

Rata-rata tidak diketahui dan diduga secara lokal dari tetangga. Pendugaan berupa kombinasi linear bobot, dengan syarat jumlah bobot = 1 (dijamin lewat pengganda Lagrange \(m\)):

\[\hat{Z}(s_0)=\sum_{i=1}^{n}\lambda_i Z(s_i),\qquad \begin{bmatrix}\boldsymbol{\Gamma} & \mathbf{1}\\ \mathbf{1}^{\top} & 0\end{bmatrix} \begin{bmatrix}\boldsymbol{\lambda}\\ m\end{bmatrix}= \begin{bmatrix}\boldsymbol{\gamma}_0\\ 1\end{bmatrix},\qquad \sigma^2_{OK}(s_0)=\boldsymbol{\lambda}^{\top}\boldsymbol{\gamma}_0+m\]

hasil_ok <- krige(curah_hujan ~ 1, locations = titik, newdata = grid,
                  model = model_fit, debug.level = 0)
hasil_ok$sd <- sqrt(hasil_ok$var1.var)

summary(hasil_ok$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2068    2152    2193    2183    2218    2276

summary(hasil_ok$sd)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   92.47   97.02   98.68   98.98   99.94  117.96

peta_pred(hasil_ok, "Ordinary Kriging - Prediksi Curah Hujan Tahunan")

peta_pred(hasil_ok, "Ordinary Kriging - Simpangan Baku Kriging (Ketidakpastian)",
          kolom = "sd", legenda = "SD kriging\n(mm)", opsi = "magma")

uji_ok <- krige(curah_hujan ~ 1, titik, titik_uji, model = model_fit, debug.level = 0)
cat(sprintf("OK di titik uji: prediksi = %.1f mm | varians = %.1f | SD = %.1f mm\n",
            uji_ok$var1.pred, uji_ok$var1.var, sqrt(uji_ok$var1.var)))

## OK di titik uji: prediksi = 2211.6 mm | varians = 10035.3 | SD = 100.2 mm

Interpretasi

Peta prediksi. Curah hujan terduga berkisar 2068–2276 mm/tahun. Pola spasialnya: makin ke selatan (korelasi prediksi dengan lintang -0.7) dan makin ke barat (korelasi dengan bujur -0.49), curah hujan makin tinggi; wilayah utara-timur (Tangerang) paling kering.
Peta ketidakpastian. Simpangan baku kriging bernilai 92.5–118 mm. SD naik seiring jarak ke titik data terdekat (korelasi r = 0.85): rata-rata 96.3 mm pada sel yang berjarak < 5 km dari titik data, dan 106.4 mm pada sel yang berjarak > 15 km — informasi yang tidak dimiliki IDW. Namun perhatikan bahwa nilai terendahnya (92.5 mm) masih besar, karena nugget yang besar menetapkan “lantai” ketidakpastian yang tidak bisa dihilangkan oleh jarak yang dekat.
Titik uji. OK menduga 2211.6 mm dengan SD 100.2 mm; selang kepercayaan 95% kasarnya 2015–2408 mm (dengan asumsi galat berdistribusi normal).

5.3 Langkah 8 — Simple Kriging (SK)

Rata-rata konstan \(\mu\) DIKETAHUI. Pada gstat, cukup mengisi argumen beta. Di sini \(\mu\) = 2150 mm (misalnya dari normal iklim jangka panjang). Nilai ini sengaja tidak persis sama dengan rerata sampel (2178.1 mm).

\[\hat{Z}(s_0)=\mu+\sum_{i=1}^{n}\lambda_i\left[Z(s_i)-\mu\right],\qquad \mathbf{C}\boldsymbol{\lambda}=\mathbf{c}_0,\quad C(h)=\text{sill}-\gamma(h)\]

mu_diketahui <- 2150

hasil_sk <- krige(curah_hujan ~ 1, locations = titik, newdata = grid,
                  model = model_fit, beta = mu_diketahui, debug.level = 0)
hasil_sk$sd <- sqrt(hasil_sk$var1.var)
summary(hasil_sk$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2065    2150    2189    2179    2215    2271

peta_pred(hasil_sk, "Simple Kriging (mu = 2150 mm) - Prediksi Curah Hujan")

uji_sk <- krige(curah_hujan ~ 1, titik, titik_uji, model = model_fit,
                beta = mu_diketahui, debug.level = 0)
cat(sprintf("SK di titik uji: prediksi = %.1f mm | SD = %.1f mm\n",
            uji_sk$var1.pred, sqrt(uji_sk$var1.var)))

## SK di titik uji: prediksi = 2210.2 mm | SD = 100.2 mm

5.3.1 Sensitivitas terhadap asumsi \(\mu\)

Kelemahan utama SK: hasilnya bergantung pada \(\mu\) yang diasumsikan. Berikut dampak perubahan \(\mu\) terhadap prediksi di titik uji dan akurasi cross-validation:

sens <- data.frame(mu = c(1800, 1950, 2100, 2150, round(mean(dat$curah_hujan), 1),
                          2300, 2500))
sens$prediksi_titik_uji <- sapply(sens$mu, function(m)
  krige(curah_hujan ~ 1, titik, titik_uji, model = model_fit, beta = m,
        debug.level = 0)$var1.pred)
sens$RMSE_LOOCV <- sapply(sens$mu, function(m) {
  cv <- krige.cv(curah_hujan ~ 1, titik, model = model_fit, nfold = nrow(titik),
                 beta = m, debug.level = 0)
  sqrt(mean(cv$residual^2))
})
sens$prediksi_titik_uji <- round(sens$prediksi_titik_uji, 1)
sens$RMSE_LOOCV <- round(sens$RMSE_LOOCV, 2)
sens

# Pembanding: OK (tanpa asumsi mu)
cat("RMSE OK (tanpa asumsi mu):", round(rmse_cv(model_fit), 2), "\n")

## RMSE OK (tanpa asumsi mu): 94.67

Interpretasi

Pada \(\mu\) = 2150 mm, SK menduga 2210.2 mm di titik uji (OK: 2211.6 mm) — keduanya nyaris sama.
RMSE SK sangat bergantung pada \(\mu\). Terendah 94.23 mm pada \(\mu\) ≈ 2178.1 mm (rerata data), hanya sedikit di bawah OK (94.67 mm) — selisih yang tidak berarti. Bila \(\mu\) meleset ≥ 200 mm, RMSE membengkak menjadi 102.47 mm (\(\mu\) = 1950) dan 102.99 mm (\(\mu\) = 2500), serta 113.36 mm pada \(\mu\) = 1800 (20% lebih buruk daripada OK). Jadi SK hampir tidak menawarkan keuntungan bila \(\mu\) tepat, tetapi berisiko besar bila \(\mu\) keliru. SK layak dipakai hanya bila \(\mu\) benar-benar diketahui dari sumber independen yang andal; bila tidak, OK lebih aman.
Di mana pengaruh \(\mu\) terasa? Pengaruh \(\mu\) pada prediksi sebesar \((1-\sum\lambda_i)\) per satuan \(\mu\). Di titik uji — yang dikelilingi banyak data — tiap kenaikan \(\mu\) 100 mm hanya menggeser prediksi 3.8 mm. Di seluruh grid, kemiringan ini rata-rata 0.11 tetapi mencapai 0.67 pada sel yang jauh dari data (korelasi dengan jarak ke titik terdekat 0.61): di sana prediksi SK “tertarik” kuat ke \(\mu\).

5.4 Langkah 9 — Universal Kriging (UK)

Rata-rata tidak konstan: berubah mengikuti tren yang dimodelkan sebagai fungsi koordinat. Di sini tren linear \(\mu(s)=\beta_0+\beta_1 x+\beta_2 y\):

\[Z(s)=\sum_{k=0}^{p}\beta_k f_k(s)+\varepsilon(s),\qquad f_0=1,\; f_1=x,\; f_2=y\]

5.4.1 Apakah tren memang ada?

reg_tren <- lm(curah_hujan ~ I(x / km) + I(y / km), data = titik)
summary(reg_tren)

## 
## Call:
## lm(formula = curah_hujan ~ I(x/km) + I(y/km), data = titik)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -184.89  -61.48   14.60   54.04  210.91 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 13064.7591  5364.4978   2.435   0.0206 *
## I(x/km)        -0.7944     0.5047  -1.574   0.1254  
## I(y/km)        -1.1184     0.5898  -1.896   0.0670 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 97.94 on 32 degrees of freedom
## Multiple R-squared:  0.2358, Adjusted R-squared:  0.188 
## F-statistic: 4.936 on 2 and 32 DF,  p-value: 0.01354

hasil_uk <- krige(curah_hujan ~ x + y, locations = titik, newdata = grid,
                  model = model_fit, debug.level = 0)
hasil_uk$sd <- sqrt(hasil_uk$var1.var)
summary(hasil_uk$var1.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2060    2143    2199    2184    2226    2294

peta_pred(hasil_uk, "Universal Kriging (tren x + y) - Prediksi Curah Hujan")

uji_uk <- krige(curah_hujan ~ x + y, titik, titik_uji, model = model_fit,
                debug.level = 0)
cat(sprintf("UK di titik uji: prediksi = %.1f mm | SD = %.1f mm\n",
            uji_uk$var1.pred, sqrt(uji_uk$var1.var)))

## UK di titik uji: prediksi = 2213.1 mm | SD = 100.2 mm

Interpretasi

Regresi tren menghasilkan R² = 0.24 dan uji-F p = 0.014: secara keseluruhan tren linear terhadap koordinat signifikan, tetapi lemah (hanya 24% variasi). Koefisien \(x\) dan \(y\) bertanda negatif → curah hujan cenderung menurun ke arah timur dan utara (konsisten dengan peta OK), meskipun masing-masing koefisien belum signifikan sendiri-sendiri (p = 0.13 dan 0.07).
Karena trennya lemah, UK menghasilkan peta yang mirip OK; prediksi titik uji 2213.1 mm vs OK 2211.6 mm. Perbedaannya terkonsentrasi jauh dari data: selisih mutlak UK−OK rata-rata 7.7 mm pada sel < 5 km dari titik data, tetapi 26.4 mm pada sel > 15 km. Di 10% sel paling barat (ujung barat provinsi), UK menduga rata-rata 25.9 mm lebih tinggi daripada OK — ekstrapolasi tren linear, risiko khas UK di tepi sebaran data yang perlu diwaspadai.
Catatan metodologis. Idealnya variogram UK dihitung dari residu tren (variogram(curah_hujan ~ x + y, ...)), bukan dari data mentah. Di sini model variogram yang sama dipakai untuk OK/SK/UK agar perbandingannya adil.

5.4.2 Variasi: tren dari elevasi (Kriging with External Drift)

Tren tidak harus fungsi koordinat; kovariat seperti elevasi dapat dipakai (curah_hujan ~ elevasi_m). Syaratnya elevasi harus tersedia di setiap lokasi prediksi. Data simulasi hanya memiliki elevasi di 35 titik, sehingga untuk keperluan latihan elevasi grid diperoleh dari interpolasi IDW (di praktik nyata gunakan DEM/SRTM).

grid_ked <- grid
grid_ked$elevasi_m <- idw(elevasi_m ~ 1, locations = titik, newdata = grid,
                          idp = 2, debug.level = 0)$var1.pred

hasil_ked <- krige(curah_hujan ~ elevasi_m, locations = titik, newdata = grid_ked,
                   model = model_fit, debug.level = 0)
cv_ked <- krige.cv(curah_hujan ~ elevasi_m, titik, model = model_fit,
                   nfold = nrow(titik), debug.level = 0)
cat("RMSE LOOCV KED:", round(sqrt(mean(cv_ked$residual^2)), 2), "mm\n")

## RMSE LOOCV KED: 95.63 mm

peta_pred(hasil_ked, "Kriging dengan Tren Elevasi (KED) - Prediksi")

Interpretasi. RMSE KED = 95.63 mm, sedikit di atas OK (94.67 mm): memasukkan elevasi sebagai tren tidak memperbaiki akurasi pada data ini, walaupun korelasinya 0.48. Sebagian besar keterkaitan curah hujan–elevasi tampaknya bekerja pada skala struktur spasial (ditangkap Co-Kriging lewat variogram-silang), bukan sebagai tren global.

5.5 Langkah 10 — Co-Kriging (CoK)

Memanfaatkan variabel sekunder (elevasi) yang berkorelasi dengan curah hujan. Kebutuhan tambahan: variogram-silang dan Linear Model of Coregionalization (LMC), yaitu semua variogram (curah hujan, elevasi, dan silang) berbagi bentuk dan range yang sama, dengan matriks sill yang semi-definit positif:

\[\hat{Z}_1(s_0)=\sum_{i=1}^{n_1}\lambda_{1i}Z_1(s_i)+\sum_{j=1}^{n_2}\lambda_{2j}Z_2(s_j),\qquad \gamma_{12}(h)=\tfrac{1}{2N(h)}\sum\left[Z_1(s_i)-Z_1(s_i+h)\right]\left[Z_2(s_i)-Z_2(s_i+h)\right]\]

# Objek gstat multivariat: hujan (primer) dan elevasi (sekunder)
g <- gstat(NULL, id = "hujan", formula = curah_hujan ~ 1, data = titik)
g <- gstat(g,    id = "elev",  formula = elevasi_m   ~ 1, data = titik)

# Variogram langsung (hujan, elevasi) dan variogram-silang (hujan.elev)
vgm_gab <- variogram(g, width = 10 * km, cutoff = 80 * km)
plot(vgm_gab, pch = 19)

# Pasang LMC: struktur (Exp, range 45 km) sama untuk ketiganya; sill diestimasi
g_lmc <- fit.lmc(vgm_gab, g, model = model_fit, fit.ranges = FALSE)
g_lmc$model

## $hujan
##   model    psill range
## 1   Nug 5423.779     0
## 2   Exp 8528.387 45000
## 
## $hujan.elev
##   model     psill range
## 1   Nug -1354.122     0
## 2   Exp  8820.548 45000
## 
## $elev
##   model     psill range
## 1   Nug  2116.567     0
## 2   Exp 11730.547 45000

# Pemeriksaan LMC: matriks nugget dan sill harus semi-definit positif
mat <- function(k) {
  m <- g_lmc$model
  matrix(c(m$hujan$psill[k], m$hujan.elev$psill[k],
           m$hujan.elev$psill[k], m$elev$psill[k]), 2)
}
cat("Nilai eigen matriks nugget:", round(eigen(mat(1))$values, 1), "\n")

## Nilai eigen matriks nugget: 5907.5 1632.9

cat("Nilai eigen matriks sill  :", round(eigen(mat(2))$values, 1), "\n")

## Nilai eigen matriks sill  : 19094.1 1164.8

plot(vgm_gab, g_lmc$model, pch = 19)

hasil_cok <- predict(g_lmc, newdata = grid, debug.level = 0)
hasil_cok$sd <- sqrt(hasil_cok$hujan.var)
summary(hasil_cok$hujan.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2043    2136    2192    2184    2225    2322

peta_pred(hasil_cok, "Co-Kriging (curah hujan + elevasi) - Prediksi Curah Hujan",
          kolom = "hujan.pred")

uji_cok <- predict(g_lmc, newdata = titik_uji, debug.level = 0)
cat(sprintf("CoK di titik uji: prediksi = %.1f mm | SD = %.1f mm\n",
            uji_cok$hujan.pred, sqrt(uji_cok$hujan.var)))

## CoK di titik uji: prediksi = 2207.4 mm | SD = 93.1 mm

Interpretasi

LMC valid: semua nilai eigen matriks nugget dan sill positif (semi-definit positif), jadi sistem Co-Kriging dapat dipecahkan.
Korelasi pada skala struktur spasial (sill-silang / akar hasil kali sill) sebesar 0.88 — jauh lebih kuat daripada korelasi Pearson mentah (0.48). Selisihnya wajar: pada skala nugget (variasi acak) korelasinya justru negatif (-0.4; nugget-silang negatif tidak masalah selama matriksnya semi-definit positif), dan gabungan keduanya (korelasi total menurut LMC ≈ 0.54) sesuai dengan korelasi Pearson. Informasi pada skala struktur inilah yang dipinjam Co-Kriging dari elevasi.
Simpangan baku kriging rata-rata turun dari 99 mm (OK) menjadi 90 mm (CoK) — Co-Kriging lebih yakin karena memakai informasi tambahan. Prediksi titik uji: CoK 2207.4 mm vs OK 2211.6 mm.
Catatan: variabel sekunder di sini diukur di titik yang sama dengan curah hujan. Keunggulan Co-Kriging biasanya jauh lebih besar bila variabel sekunder tersedia lebih rapat (mis. citra satelit atau DEM).

5.6 Langkah 11 — Pembanding: IDW

hasil_idw <- idw(curah_hujan ~ 1, locations = titik, newdata = grid,
                 idp = 2, debug.level = 0)
peta_pred(hasil_idw, "IDW (p = 2) - Prediksi Curah Hujan")

6 Evaluasi

6.1 Langkah 12 — Cross-Validation (Leave-One-Out)

Satu titik dikeluarkan, diduga dari 34 titik lain, lalu dibandingkan dengan nilai sebenarnya; diulang untuk seluruh titik.

\[\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left[z(s_i)-\hat{z}_{-i}(s_i)\right]^2},\qquad \text{MAE}=\frac{1}{n}\sum_{i=1}^{n}\left|z(s_i)-\hat{z}_{-i}(s_i)\right|\]

ME (mean error) = rerata (pengamatan − prediksi); nilai positif berarti prediksi cenderung terlalu rendah. R² = kuadrat korelasi pengamatan–prediksi.

metrik <- function(obs, pred) {
  gal <- obs - pred
  c(RMSE = sqrt(mean(gal^2)), MAE = mean(abs(gal)), ME = mean(gal),
    R2 = cor(obs, pred)^2)
}
n <- nrow(titik)

# Ordinary, Simple, Universal: krige.cv()
cv_ok <- krige.cv(curah_hujan ~ 1,     titik, model = model_fit, nfold = n, debug.level = 0)
cv_sk <- krige.cv(curah_hujan ~ 1,     titik, model = model_fit, nfold = n,
                  beta = mu_diketahui, debug.level = 0)
cv_uk <- krige.cv(curah_hujan ~ x + y, titik, model = model_fit, nfold = n, debug.level = 0)

# Co-Kriging: gstat.cv() pada objek multivariat. Secara bawaan hanya pengamatan
# curah hujan di titik ke-i yang dikeluarkan; elevasi titik tersebut tetap dipakai
# (sesuai skenario nyata: elevasi dapat diukur/diketahui di titik target).
cv_cok <- suppressMessages(gstat.cv(g_lmc, nfold = n, debug.level = 0))

# IDW: leave-one-out manual
pred_idw <- sapply(seq_len(n), function(i)
  idw(curah_hujan ~ 1, titik[-i, ], titik[i, ], idp = 2, debug.level = 0)$var1.pred)

tab_cv <- rbind(
  "IDW (p = 2)"       = metrik(titik$curah_hujan, pred_idw),
  "Simple Kriging"    = metrik(titik$curah_hujan, cv_sk$var1.pred),
  "Ordinary Kriging"  = metrik(titik$curah_hujan, cv_ok$var1.pred),
  "Universal Kriging" = metrik(titik$curah_hujan, cv_uk$var1.pred),
  "Co-Kriging"        = metrik(titik$curah_hujan, cv_cok$hujan.pred)
)
tab_cv <- as.data.frame(round(tab_cv, 2))
tab_cv$Peringkat_RMSE <- rank(tab_cv$RMSE)
tab_cv[order(tab_cv$RMSE), ]

cat(sprintf("z-score OK : rerata = %.3f ; varians = %.3f  (ideal: 0 dan 1)\n",
            mean(cv_ok$zscore), var(cv_ok$zscore)))

## z-score OK : rerata = 0.005 ; varians = 0.928  (ideal: 0 dan 1)

ggplot(data.frame(obs = titik$curah_hujan, pred = cv_ok$var1.pred),
       aes(obs, pred)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_point(colour = "#00207E", size = 2.5) +
  labs(title = "Cross-Validation Ordinary Kriging: Pengamatan vs Prediksi",
       x = "Curah hujan pengamatan (mm)", y = "Curah hujan prediksi (mm)") +
  coord_equal() + theme_minimal()

6.1.1 Apakah selisih antarmetode bermakna?

Dengan hanya 35 titik, selisih RMSE beberapa mm bisa terjadi karena kebetulan. Galat LOOCV tiap titik dapat di-bootstrap (penarikan ulang 35 titik dengan pengembalian) untuk menaksir selang kepercayaan selisih RMSE Ordinary − Co-Kriging (nilai positif = Co-Kriging lebih baik). Catatan: galat LOOCV tidak sepenuhnya saling bebas, sehingga hasil ini hanya indikatif.

set.seed(2025)
res_ok  <- cv_ok$residual
res_cok <- cv_cok$residual
rmse_v  <- function(e) sqrt(mean(e^2))

selisih_boot <- replicate(5000, {
  i <- sample(n, replace = TRUE)
  rmse_v(res_ok[i]) - rmse_v(res_cok[i])
})
ci_boot <- quantile(selisih_boot, c(0.025, 0.975))
cat(sprintf("Selisih RMSE (OK - CoK) = %.2f mm | selang bootstrap 95%%: [%.2f ; %.2f] | P(CoK lebih baik) = %.2f\n",
            rmse_v(res_ok) - rmse_v(res_cok), ci_boot[1], ci_boot[2], mean(selisih_boot > 0)))

## Selisih RMSE (OK - CoK) = 3.44 mm | selang bootstrap 95%: [-7.39 ; 13.17] | P(CoK lebih baik) = 0.73

# Distribusi galat mutlak
galat <- data.frame(
  Metode = c("Ordinary", "Co-Kriging"),
  Median = c(median(abs(res_ok)), median(abs(res_cok))),
  P90    = c(quantile(abs(res_ok), 0.9), quantile(abs(res_cok), 0.9)),
  Maks   = c(max(abs(res_ok)), max(abs(res_cok)))
)
galat[, -1] <- round(galat[, -1], 1)
galat

cat("Titik dengan galat CoK lebih kecil daripada OK:", sum(abs(res_cok) < abs(res_ok)), "dari", n, "\n")

## Titik dengan galat CoK lebih kecil daripada OK: 16 dari 35

# Bukti: jalur slide (data tetap dalam derajat, jarak great-circle dalam km)
# memberi hasil yang sama dengan jalur UTM di atas.
titik_deg <- st_as_sf(dat, coords = c("lon", "lat"), crs = 4326)
cv_deg <- krige.cv(curah_hujan ~ 1, titik_deg, nfold = n, debug.level = 0,
                   model = vgm(psill = 6500, model = "Exp", range = 45, nugget = 7000))
cat("RMSE OK jalur derajat (range = 45 km):", round(sqrt(mean(cv_deg$residual^2)), 2), "\n")

## RMSE OK jalur derajat (range = 45 km): 94.67

cat("RMSE OK jalur UTM     (range = 45000 m):", round(tab_cv["Ordinary Kriging", "RMSE"], 2), "\n")

## RMSE OK jalur UTM     (range = 45000 m): 94.67

Interpretasi

Urutan RMSE (kecil = lebih akurat): Co-Kriging (91.23) paling kecil, diikuti Simple (94.61), Ordinary (94.67), Universal (95.81), dan IDW (96.57).
Simple ≈ Ordinary dalam RMSE. Selisihnya hanya 0.06 mm karena \(\mu\) = 2150 mm dekat dengan rerata data (2178 mm). Namun SK memiliki bias (ME = 5.48 mm vs OK 0.89 mm): \(\mu\) yang sedikit terlalu rendah menarik prediksi ke bawah. Sensitivitas pada Langkah 8 menunjukkan keunggulan RMSE ini rapuh: SK hanya setara OK bila \(\mu\) tepat, dan jauh lebih buruk bila meleset.
Universal sedikit di bawah Ordinary (95.81 vs 94.67) karena tren linear \(x+y\) lemah (R² 0.24); parameter tren tambahan tidak sebanding dengan manfaatnya pada n = 35.
Co-Kriging memiliki RMSE terkecil, 3.6% lebih rendah dari Ordinary — namun perbedaan ini tidak bermakna secara statistik: selang bootstrap 95% untuk selisih RMSE (-7.4 s.d. 13.2 mm) mencakup nol. Galat terbesar memang lebih kecil (maks 188.5 vs 239.8 mm), tetapi MAE-nya (76.27) justru lebih besar dari Ordinary (73.47), dan hanya 16 dari 35 titik yang galatnya lebih kecil dengan CoK. Jadi keunggulan CoK tidak konsisten; kesimpulan sebaiknya memakai beberapa metrik, bukan RMSE saja.
Kriging vs IDW. OK 2% lebih rendah RMSE-nya daripada IDW — selisih kecil. IDW juga memiliki bias yang lebih besar (ME = 10.45 mm). Keunggulan utama Kriging di sini bukan pada akurasi titik, melainkan pada kuantifikasi ketidakpastian (peta SD) yang tidak dimiliki IDW. Diagnostik z-score OK (rerata 0, varians 0.93) mendekati nilai ideal 0 dan 1, artinya ketidakpastian yang dilaporkan model cukup dapat dipercaya.
Prediksi “menyusut” ke rerata. Pada diagram pencar OK, kemiringan garis prediksi terhadap pengamatan hanya 0.24 (ideal 1): nilai ekstrem tidak tertangkap. Contohnya 5 titik dengan curah hujan ≥ 2300 mm (rerata pengamatan 2362 mm) hanya diduga rata-rata 2202 mm. Ini ciri khas Kriging (penghalusan) terlebih dengan nugget besar: peta Kriging cocok untuk pola umum, bukan untuk menaksir nilai ekstrem.
Keterbatasan. Dengan n = 35, selisih terjauh antarmetode hanya 5.9% (RMSE 91.23–96.57 mm) sehingga perbedaan peringkat jangan dilebih-lebihkan. Variogram dan LMC juga diestimasi dari data yang sama dengan yang divalidasi (praktik lazim pada LOOCV, tetapi sedikit optimistis).

6.2 Langkah 13 — Perbandingan Peta Empat Metode

gabung <- bind_rows(
  data.frame(st_coordinates(hasil_ok),  pred = hasil_ok$var1.pred,   metode = "Ordinary"),
  data.frame(st_coordinates(hasil_sk),  pred = hasil_sk$var1.pred,   metode = "Simple"),
  data.frame(st_coordinates(hasil_uk),  pred = hasil_uk$var1.pred,   metode = "Universal"),
  data.frame(st_coordinates(hasil_cok), pred = hasil_cok$hujan.pred, metode = "Co-Kriging")
)
gabung$metode <- factor(gabung$metode,
                        levels = c("Ordinary", "Simple", "Universal", "Co-Kriging"))
ggplot(gabung, aes(X, Y, fill = pred)) +
  geom_tile() +
  geom_sf(data = batas_utm, inherit.aes = FALSE, fill = NA, colour = "grey20",
          linewidth = 0.2) +
  facet_wrap(~ metode, ncol = 2) +
  scale_fill_viridis_c(option = "turbo", name = "mm/tahun") +
  labs(title = "Perbandingan Peta Prediksi Curah Hujan: Empat Metode Kriging",
       x = NULL, y = NULL) +
  theme_minimal()

Rata-rata prediksi per kabupaten/kota (zonal) — bentuk keluaran yang paling sering dibutuhkan untuk pelaporan:

grid_kab <- st_join(hasil_ok, batas_utm["Kabupaten"])
zonal <- data.frame(
  Kabupaten = grid_kab$Kabupaten,
  OK  = hasil_ok$var1.pred,
  SK  = hasil_sk$var1.pred,
  UK  = hasil_uk$var1.pred,
  CoK = hasil_cok$hujan.pred
) %>%
  group_by(Kabupaten) %>%
  summarise(across(c(OK, SK, UK, CoK), ~ round(mean(.x), 1)), n_sel = n()) %>%
  arrange(desc(OK))
zonal

Interpretasi

Keempat peta memiliki pola yang sama. Rata-rata per kabupaten/kota tertinggi di Lebak (2214.1 mm) dan terendah di Tangerang (2096.3 mm) menurut OK; wilayah tertinggi dan terendah sama pada keempat metode. Selisih antarmetode untuk satu wilayah rata-rata 25.4 mm (setara 23% SD data), dengan selisih terbesar 60.7 mm di Kota Cilegon — wilayah kecil (hanya 43 sel grid). Jadi kesimpulan spasial umum tidak bergantung pada varian Kriging yang dipilih, tetapi untuk wilayah kecil pilihan metode dapat menggeser angka rata-ratanya cukup terasa.

7 Contoh Soal Manual: Kriging 3 Titik dengan Matriks

Menghitung bobot Kriging “dengan tangan” di R agar mekanismenya dipahami. Tiga titik: A(0,0) = 1800 mm; B(4,0) = 2200 mm; C(2,3) = 2000 mm; target P(1,1). Model Exponential: \(c_0=0\), \(c_1=250000\), \(a=3\) → \(\gamma(h)=250000\left(1-e^{-h/3}\right)\).

pt  <- data.frame(id = c("A", "B", "C"), x = c(0, 4, 2), y = c(0, 0, 3),
                  z = c(1800, 2200, 2000))
P   <- c(x = 1, y = 1)
gam <- function(h) 250000 * (1 - exp(-h / 3))

D   <- as.matrix(dist(pt[, c("x", "y")]))                # jarak antartitik
d0  <- sqrt((pt$x - P["x"])^2 + (pt$y - P["y"])^2)       # jarak titik ke P
G   <- gam(D); diag(G) <- 0                              # matriks gamma
g0  <- gam(d0)                                           # vektor gamma ke P
n3  <- nrow(pt)

round(D, 3); round(d0, 3)

##       1     2     3
## 1 0.000 4.000 3.606
## 2 4.000 0.000 3.606
## 3 3.606 3.606 0.000

## [1] 1.414 3.162 2.236

7.0.1 Ordinary Kriging

A_ok <- rbind(cbind(G, 1), c(rep(1, n3), 0))             # sistem (n+1) x (n+1)
sol  <- solve(A_ok, c(g0, 1))
lam  <- sol[1:n3]; m <- sol[n3 + 1]

z_ok   <- sum(lam * pt$z)
var_ok <- sum(lam * g0) + m
cat("Bobot lambda :", round(lam, 4), "| jumlah =", round(sum(lam), 4), "\n")

## Bobot lambda : 0.5322 0.1579 0.3099 | jumlah = 1

cat("Lagrange m   :", round(m, 2), "\n")

## Lagrange m   : 10709.68

cat("Prediksi     :", round(z_ok, 1), "mm\n")

## Prediksi     : 1925.1 mm

cat("Varians      :", round(var_ok, 1), "| SD =", round(sqrt(var_ok), 1), "mm\n")

## Varians      : 127146.6 | SD = 356.6 mm

7.0.2 Simple Kriging (\(\mu\) = 2000)

C0   <- 250000
Cmat <- C0 - G                                            # C(h) = sill - gamma(h)
c0   <- C0 - g0
lam_sk <- solve(Cmat, c0)
mu     <- 2000

z_sk   <- mu + sum(lam_sk * (pt$z - mu))
var_sk <- C0 - sum(lam_sk * c0)
cat("Bobot lambda :", round(lam_sk, 4), "| jumlah =", round(sum(lam_sk), 4), "\n")

## Bobot lambda : 0.5045 0.1302 0.2837 | jumlah = 0.9185

cat("Prediksi     :", round(z_sk, 1), "mm\n")

## Prediksi     : 1925.1 mm

cat("Varians      :", round(var_sk, 1), "| SD =", round(sqrt(var_sk), 1), "mm\n")

## Varians      : 126273.3 | SD = 355.3 mm

7.0.3 Verifikasi dengan `gstat` dan pembanding Pertemuan 5

pt_sf  <- st_as_sf(pt, coords = c("x", "y"))
P_sf   <- st_as_sf(data.frame(x = 1, y = 1), coords = c("x", "y"))
mod_ct <- vgm(psill = 250000, "Exp", range = 3, nugget = 0)
cek_ok <- krige(z ~ 1, pt_sf, P_sf, model = mod_ct, debug.level = 0)
cek_sk <- krige(z ~ 1, pt_sf, P_sf, model = mod_ct, beta = 2000, debug.level = 0)

cat(sprintf("gstat OK: %.1f mm | varians %.1f   (manual: %.1f | %.1f)\n",
            cek_ok$var1.pred, cek_ok$var1.var, z_ok, var_ok))

## gstat OK: 1925.1 mm | varians 127146.6   (manual: 1925.1 | 127146.6)

cat(sprintf("gstat SK: %.1f mm | varians %.1f   (manual: %.1f | %.1f)\n",
            cek_sk$var1.pred, cek_sk$var1.var, z_sk, var_sk))

## gstat SK: 1925.1 mm | varians 126273.3   (manual: 1925.1 | 126273.3)

w <- 1 / d0^2
cat("IDW (p = 2)      :", round(sum(w * pt$z) / sum(w), 1), "mm\n")

## IDW (p = 2)      : 1900 mm

cat("Nearest Neighbor :", pt$z[which.min(d0)], "mm (titik terdekat:", pt$id[which.min(d0)], ")\n")

## Nearest Neighbor : 1800 mm (titik terdekat: A )

Interpretasi

Perhitungan matriks manual identik dengan keluaran gstat (prediksi dan varians sama), membuktikan bahwa krige() hanya menyelesaikan sistem persamaan linear di atas.
Bobot OK berjumlah tepat 1 (1); bobot SK berjumlah 0.9185 (< 1) — selisihnya 0.0815 “ditutup” oleh \(\mu\). Inilah perbedaan mendasar OK vs SK.
Titik A paling berpengaruh (bobot 0.53) karena paling dekat dengan P, sedangkan B berbobot terkecil (0.16) karena paling jauh. Bandingkan dengan bobot IDW (p = 2): A 0.62, B 0.12, C 0.25. Kriging memindahkan sebagian bobot dari A ke B dan C, sebab ia memperhitungkan posisi antartitik (bukan hanya jarak ke P): B dan C terletak di arah yang berbeda dari A sehingga membawa informasi yang tidak redundan.
Prediksi OK dan SK sama (1925.1 mm) dan sedikit di atas IDW (1900 mm); Nearest Neighbor hanya menyalin titik terdekat (1800 mm). Berbeda dari semua metode di Pertemuan 5, Kriging menyertakan SD kesalahan prediksi: 356.6 mm. SD yang besar wajar karena hanya ada 3 titik dengan nilai yang beragam.

8 Kesimpulan

Autokorelasi spasial ada (semivarians naik dari ≈ 6981 ke ≈ 14395), dengan efek nugget besar (52% dari sill) sehingga ketergantungan spasial tergolong sedang.
Model variogram Exponential (nugget 7000, partial sill 6500, range 45 km) memadai; fit.variogram() otomatis tidak konvergen pada data ini sehingga range harus ditetapkan dari pembacaan visual.
Ordinary Kriging menghasilkan peta dan peta ketidakpastian yang konsisten dengan pola geografis (lebih basah ke selatan-barat, lebih kering di utara-timur).
Menurut LOOCV, Co-Kriging memiliki RMSE terkecil (91.23 mm), tetapi selisih antarmetode kecil (91.23–96.57 mm), tidak bermakna secara statistik (selang bootstrap selisih OK−CoK: -7.4 s.d. 13.2 mm), dan MAE tidak searah dengan RMSE. Pilihan varian Kriging sebaiknya mengikuti karakteristik data: OK sebagai pilihan bawaan yang aman; SK hanya bila \(\mu\) sungguh diketahui; UK bila tren kuat; CoK bila variabel sekunder tersedia lebih rapat.
Keunggulan utama Kriging atas IDW adalah kuantifikasi ketidakpastian, bukan semata-mata akurasi titik.

Catatan. Seluruh data adalah simulasi untuk keperluan pembelajaran; angka di atas tidak mewakili kondisi curah hujan Banten yang sebenarnya.

9 Informasi Sesi

sessionInfo()

## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8  LC_CTYPE=English_Indonesia.utf8   
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=English_Indonesia.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] viridis_0.6.5     viridisLite_0.4.3 readxl_1.4.3      ggplot2_4.0.3    
## [5] dplyr_1.1.4       gstat_2.1-4       sf_1.0-19        
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.9         generics_0.1.3     class_7.3-22       KernSmooth_2.23-24
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.4.2         RColorBrewer_1.1-3 fastmap_1.2.0      cellranger_1.1.0  
## [13] jsonlite_1.9.0     e1071_1.7-16       DBI_1.2.3          gridExtra_2.3     
## [17] scales_1.4.0       jquerylib_0.1.4    abind_1.4-8        cli_3.6.4         
## [21] rlang_1.1.7        units_0.8-5        intervals_0.15.5   withr_3.0.3       
## [25] cachem_1.1.0       yaml_2.3.10        FNN_1.1.4.1        parallel_4.4.2    
## [29] tools_4.4.2        spacetime_1.3-3    vctrs_0.6.5        R6_2.6.1          
## [33] zoo_1.8-12         proxy_0.4-27       lifecycle_1.0.5    classInt_0.4-11   
## [37] pkgconfig_2.0.3    pillar_1.10.1      bslib_0.9.0        gtable_0.3.6      
## [41] glue_1.8.0         Rcpp_1.0.14        xfun_0.51          tibble_3.2.1      
## [45] tidyselect_1.2.1   rstudioapi_0.17.1  knitr_1.49         farver_2.1.2      
## [49] htmltools_0.5.8.1  labeling_0.4.3     rmarkdown_2.29     xts_0.14.1        
## [53] compiler_4.4.2     S7_0.2.1           stars_0.6-8        sp_2.2-0

Kriging untuk Interpolasi Spasial: Curah Hujan Tahunan di Provinsi Banten

Studi Kasus Pertemuan 6 — Mata Kuliah Statistika Lingkungan (STL623301)

Miftahus Sholihin, M.Si

27 October 2026

1 Pendahuluan

1.1 Tentang Studi Kasus

1.2 Rumusan Pertanyaan Analisis

2 Persiapan

2.1 Memuat Package

2.2 Langkah 1 — Import Data dan Peta

Interpretasi

2.3 Langkah 2 — Proyeksi Koordinat

3 Eksplorasi Data

3.1 Langkah 3 — Statistik Deskriptif dan Korelasi

Interpretasi

4 Variogram

4.1 Langkah 4 — Variogram Empiris

Interpretasi

4.2 Langkah 5 — Model Variogram Teoritis

4.2.1 Mengapa fit.variogram() otomatis perlu diwaspadai

4.2.2 Solusi praktis: tetapkan range secara visual

Interpretasi

5 Prediksi Spasial

5.1 Langkah 6 — Grid Prediksi

5.2 Langkah 7 — Ordinary Kriging (OK)

Interpretasi

5.3 Langkah 8 — Simple Kriging (SK)

5.3.1 Sensitivitas terhadap asumsi \(\mu\)

Interpretasi

5.4 Langkah 9 — Universal Kriging (UK)

5.4.1 Apakah tren memang ada?

Interpretasi

5.4.2 Variasi: tren dari elevasi (Kriging with External Drift)

5.5 Langkah 10 — Co-Kriging (CoK)

Interpretasi

5.6 Langkah 11 — Pembanding: IDW

6 Evaluasi

6.1 Langkah 12 — Cross-Validation (Leave-One-Out)

6.1.1 Apakah selisih antarmetode bermakna?

Interpretasi

6.2 Langkah 13 — Perbandingan Peta Empat Metode

Interpretasi

7 Contoh Soal Manual: Kriging 3 Titik dengan Matriks

7.0.1 Ordinary Kriging

7.0.2 Simple Kriging (\(\mu\) = 2000)

7.0.3 Verifikasi dengan gstat dan pembanding Pertemuan 5

Interpretasi

8 Kesimpulan

9 Informasi Sesi

4.2.1 Mengapa `fit.variogram()` otomatis perlu diwaspadai

7.0.3 Verifikasi dengan `gstat` dan pembanding Pertemuan 5