1 Mengenal interpolasi titik

Interpolasi titik adalah proses menduga nilai suatu variabel pada lokasi yang belum diamati menggunakan pengamatan di lokasi lain. Contohnya, suhu di lokasi tanpa alat ukur diperkirakan menggunakan suhu beberapa stasiun di sekitarnya.

Capaian praktikum

  1. Menjelaskan prinsip setiap metode.
  2. Menghitung prediksi satu lokasi dan membuat peta prediksi.
  3. Membandingkan kinerja pada data yang tidak digunakan untuk membentuk prediksi.
  4. Menggabungkan beberapa metode tanpa kebocoran data.
  5. Menginterpretasikan angka, peta, dan ukuran kesalahan dengan satuan yang tepat.

1.1 Konsep yang perlu dipahami

  • Lokasi yang dekat sering memiliki karakteristik mirip. Namun, kedekatan tidak menjamin kesamaan; pertimbangkan sifat fenomena yang dipelajari.
  • Tidak semua metode menghasilkan permukaan yang halus atau kontinu. NN, Thiessen, dan local averaging dapat berubah mendadak.
  • Interpolasi menduga nilai di antara lokasi pengamatan. Ekstrapolasi menduga di luar cakupan lokasi tersebut dan perlu ditafsirkan lebih hati-hati.
  • Koordinat latihan menggunakan sistem Kartesian lokal dalam kilometer, bukan lintang dan bujur. Untuk data nyata, gunakan proyeksi yang sesuai sebelum menghitung jarak Euclidean; selisih derajat bukan jarak kilometer.
  • Data dalam modul adalah simulasi untuk pembelajaran, bukan hasil pengukuran wilayah tertentu.

2 Persiapan praktikum

2.1 Sintaks

# Semua kode muncul pada hasil Knit.
knitr::opts_chunk$set(
  echo = TRUE, eval = TRUE, message = FALSE, warning = TRUE,
  fig.width = 8, fig.height = 5.5, fig.align = "center",
  dpi = 120, cache = FALSE, error = FALSE
)
options(width = 85, digits = 5)

# Instal hanya paket yang belum tersedia.
paket <- c("deldir", "knitr", "rmarkdown")
belum_ada <- paket[!vapply(paket, requireNamespace, logical(1), quietly = TRUE)]
if (length(belum_ada) > 0) {
  install.packages(belum_ada, repos = "https://cloud.r-project.org")
}
stopifnot(all(vapply(paket, requireNamespace, logical(1), quietly = TRUE)))

2.1.1 Data latihan dan target

Wilayah latihan berbentuk persegi 10 × 10 km. Empat stasiun berada di sudut wilayah dan 32 stasiun berada di bagian dalam.

set.seed(1352)
stasiun <- data.frame(
  id = sprintf("S%02d", 1:36),
  x = c(0, 10, 10, 0, runif(32, 0.5, 9.5)),
  y = c(0, 0, 10, 10, runif(32, 0.5, 9.5))
)
# Tren spasial, variasi lokal, dan galat acak
stasiun$suhu <- with(stasiun,
  25 + 0.45*x - 0.25*y +
    2*exp(-((x - 7)^2 + (y - 3)^2)/8) + rnorm(36, 0, 0.6)
)
target <- data.frame(x = 5, y = 5)
stopifnot(
  all(complete.cases(stasiun)),
  all(is.finite(as.matrix(stasiun[c("x", "y", "suhu")]))),
  !anyDuplicated(stasiun[c("x", "y")])
)
knitr::kable(head(stasiun, 8), digits = 2,
             caption = "Delapan baris pertama data suhu simulasi")
Delapan baris pertama data suhu simulasi
id x y suhu
S01 0.00 0.00 23.76
S02 10.00 0.00 29.63
S03 10.00 10.00 27.73
S04 0.00 10.00 21.68
S05 5.78 7.23 27.87
S06 9.16 6.92 27.83
S07 1.82 5.57 24.20
S08 2.60 4.48 24.73 
summary(stasiun[c("x", "y", "suhu")])
##        x               y              suhu     
##  Min.   : 0.00   Min.   : 0.00   Min.   :21.7  
##  1st Qu.: 2.32   1st Qu.: 2.91   1st Qu.:25.0  
##  Median : 4.88   Median : 5.21   Median :26.2  
##  Mean   : 4.84   Mean   : 5.33   Mean   :26.2  
##  3rd Qu.: 7.13   3rd Qu.: 7.68   3rd Qu.:27.6  
##  Max.   :10.00   Max.   :10.00   Max.   :29.6
plot(stasiun$x, stasiun$y, pch = 21, bg = "#1495a7",
     xlim = c(0, 10), ylim = c(0, 10), asp = 1,
     xlab = "Koordinat x (km)", ylab = "Koordinat y (km)",
     main = "Lokasi stasiun dan target")
text(stasiun$x, stasiun$y, labels = stasiun$id, pos = 3, cex = 0.6)
points(target$x, target$y, pch = 8, col = "#d44a30", cex = 1.8, lwd = 2)
legend("topleft", c("Stasiun", "Target (5,5)"),
       pch = c(21, 8), pt.bg = "#1495a7",
       col = c("black", "#d44a30"), bty = "n")

2.1.2 Grid dan fungsi bantuan peta

Grid adalah kumpulan lokasi yang akan diprediksi. Grid lebih rapat membuat gambar lebih rinci, tetapi tidak otomatis meningkatkan akurasi model.

gx <- seq(0, 10, length.out = 51)
gy <- seq(0, 10, length.out = 51)
grid <- expand.grid(x = gx, y = gy)

jarak_ke <- function(data, titik) {
  sqrt((data$x - titik$x)^2 + (data$y - titik$y)^2)
}

# Semua peta memakai skala warna yang sama.
palet <- grDevices::hcl.colors(60, "YlOrRd", rev = TRUE)
batas_warna <- range(stasiun$suhu)
peta <- function(prediksi, judul) {
  stopifnot(length(prediksi) == nrow(grid))
  image(gx, gy, matrix(prediksi, nrow = length(gx)),
        col = palet, zlim = batas_warna, asp = 1,
        xlab = "Koordinat x (km)", ylab = "Koordinat y (km)", main = judul)
  contour(gx, gy, matrix(prediksi, nrow = length(gx)),
          add = TRUE, drawlabels = FALSE, col = "#55555555")
  points(stasiun$x, stasiun$y, pch = 21, bg = "white", cex = 0.65)
  points(target$x, target$y, pch = 8, col = "black", cex = 1.3)
  box()
}
# Legenda warna bersama
image(seq(0, 1, length.out = 60), 1,
      matrix(seq(batas_warna[1], batas_warna[2], length.out = 60), 60, 1),
      col = palet, axes = FALSE, xlab = "Suhu (derajat Celsius)", ylab = "")
axis(1, at = seq(0, 1, length.out = 5),
     labels = round(seq(batas_warna[1], batas_warna[2], length.out = 5), 1))

2.2 Interpretasi

Kolom x dan y menyatakan lokasi dalam km, sedangkan suhu adalah variabel yang diprediksi dalam °C. Tanda bintang menunjukkan target yang tidak memiliki suhu pengamatan. Warna peta menunjukkan suhu dugaan, bukan ketidakpastian.

Peta pada bagian demonstrasi menggunakan semua stasiun. Penilaian akurasi dilakukan secara terpisah menggunakan pembagian data latih dan data uji.

3 Inverse Distance Weighting (IDW)

3.1 Pengantar

IDW memberi bobot lebih besar kepada stasiun yang lebih dekat:

\[ \widehat{Z}(s_0)= \frac{\sum_{i=1}^{n}d(s_0,s_i)^{-p}Z(s_i)} {\sum_{i=1}^{n}d(s_0,s_i)^{-p}}. \]

Parameter \(p>0\) mengatur penurunan pengaruh terhadap jarak. Ketika \(p\) semakin besar, titik paling dekat semakin dominan. Nilai \(p=2\) dipakai untuk contoh, bukan nilai yang selalu terbaik.

3.2 Sintaks

3.2.1 Contoh perhitungan IDW

Misalkan tersedia tiga stasiun A, B, dan C dengan suhu masing-masing 25 °C, 30 °C, dan 28 °C. Jarak setiap stasiun ke lokasi target berturut-turut adalah 2 km, 3 km, dan 5 km. Dengan parameter \(p=2\), suhu di lokasi target dihitung sebagai berikut.

contoh_idw <- data.frame(
  titik = c("A", "B", "C"),
  suhu = c(25, 30, 28), jarak = c(2, 3, 5)
)
contoh_idw$bobot_awal <- 1 / contoh_idw$jarak^2
contoh_idw$bobot <- contoh_idw$bobot_awal / sum(contoh_idw$bobot_awal)
contoh_idw$kontribusi <- contoh_idw$bobot * contoh_idw$suhu
knitr::kable(contoh_idw, digits = 4)
titik suhu jarak bobot_awal bobot kontribusi
A 25 2 0.2500 0.6233 15.5817
B 30 3 0.1111 0.2770 8.3102
C 28 5 0.0400 0.0997 2.7922 
idw_manual <- sum(contoh_idw$kontribusi)
idw_manual
## [1] 26.684

3.2.2 Fungsi dan peta IDW

Jika target tepat sama dengan lokasi stasiun, fungsi mengembalikan nilai pengamatan untuk menghindari pembagian dengan nol.

pred_idw <- function(data, baru, p = 2) {
  stopifnot(p > 0)
  vapply(seq_len(nrow(baru)), function(j) {
    d <- jarak_ke(data, baru[j, ])
    if (any(d < 1e-10)) return(mean(data$suhu[d < 1e-10]))
    # Setara dengan 1/d^p setelah normalisasi.
    w <- (min(d) / d)^p
    sum(w * data$suhu) / sum(w)
  }, numeric(1))
}
idw_target <- pred_idw(stasiun, target, p = 2)
idw_grid <- pred_idw(stasiun, grid, p = 2)
idw_target
## [1] 26.801
peta(idw_grid, "Interpolasi IDW (p = 2)")

knitr::kable(data.frame(
  p = c(1, 2, 3),
  prediksi = vapply(c(1, 2, 3),
    function(p) pred_idw(stasiun, target, p), numeric(1))
), digits = 3)
p prediksi
1 26.516
2 26.801
3 26.979

3.3 Interpretasi

Pada contoh tiga stasiun, metode IDW dengan parameter \(p=2\) menghasilkan dugaan suhu sebesar 26,684 °C. Stasiun A memberikan pengaruh terbesar karena jaraknya paling dekat dengan lokasi target. Dengan demikian, kontribusi setiap stasiun ditentukan oleh jaraknya, bukan hanya oleh besar suhunya.

Pada data simulasi, suhu di lokasi (5,5) diprediksi sebesar 26,8 °C menggunakan metode IDW dengan \(p=2\). Angka ini merupakan perkiraan suhu pada lokasi yang belum diamati, berdasarkan suhu stasiun-stasiun di sekitarnya.

4 Tetangga terdekat (Nearest Neighbor)

4.1 Pengantar

Metode ini mengambil nilai dari satu stasiun terdekat, tanpa rata-rata beberapa stasiun:

\[ \widehat Z(s_0)=Z(s_{i^*}),\qquad i^*=\arg\min_i d(s_0,s_i). \]

4.2 Sintaks

# Ilustrasi slide 10–11
contoh_nn <- data.frame(
  id = c("A", "B", "C"),
  x = c(1, 4, 6), y = c(2, 4, 1), suhu = c(25, 30, 28)
)
contoh_nn$jarak <- jarak_ke(contoh_nn, data.frame(x = 3, y = 3))
knitr::kable(contoh_nn, digits = 3)
id x y suhu jarak
A 1 2 25 2.236
B 4 4 30 1.414
C 6 1 28 3.606 
contoh_nn$suhu[which.min(contoh_nn$jarak)]
## [1] 30
pred_nn <- function(data, baru) {
  vapply(seq_len(nrow(baru)), function(j) {
    d <- jarak_ke(data, baru[j, ])
    data$suhu[which.min(d)]
  }, numeric(1))
}
# Jika jarak seri, which.min memilih baris pertama.
tetangga <- which.min(jarak_ke(stasiun, target))
knitr::kable(stasiun[tetangga, ], digits = 3)
id x y suhu
10 S10 4.833 5.636 27.133 
nn_target <- pred_nn(stasiun, target)
nn_grid <- pred_nn(stasiun, grid)
nn_target
## [1] 27.133
peta(nn_grid, "Interpolasi tetangga terdekat")

4.3 Interpretasi

Pada contoh tiga stasiun, stasiun B merupakan lokasi terdekat dari target (3,3), dengan jarak sekitar 1,414 satuan. Oleh karena itu, suhu target diprediksi sebesar 30 °C, sama dengan suhu di stasiun B.

Pada data simulasi, stasiun terdekat dari target (5,5) adalah S10, sehingga suhu target diprediksi sebesar 27,13 °C. Metode Nearest Neighbor mengambil langsung nilai dari satu stasiun terdekat, tanpa menghitung rata-rata atau mempertimbangkan suhu stasiun lainnya.

Pada peta, setiap daerah berwarna seragam menunjukkan lokasi-lokasi yang memiliki stasiun terdekat yang sama. Perubahan warna dapat terjadi secara mendadak ketika berpindah ke daerah pengaruh stasiun lain.

5 Triangulasi Delaunay dan interpolasi linier

5.1 Pengantar

Triangulasi Delaunay membentuk jaringan segitiga dari lokasi pengamatan. Setelah ditemukan segitiga yang memuat target, prediksi dihitung dari tiga simpul dengan bobot barycentric:

\[ \widehat Z(s_0)=w_A Z_A+w_B Z_B+w_C Z_C,\qquad w_A+w_B+w_C=1. \]

Bobot dapat diperoleh dari perbandingan luas subsegitiga dengan luas segitiga utama. Untuk target di dalam segitiga, bobot tidak negatif. Modul tidak melakukan ekstrapolasi. Target di luar convex hull, yaitu batas cembung terluar lokasi stasiun, menghasilkan NA.

5.2 Sintaks

5.2.1 Contoh Perhitungan Interpolasi Linier dalam Segitiga

abc <- data.frame(
  id = c("A", "B", "C"),
  x = c(1, 4, 2), y = c(1, 1, 4), suhu = c(20, 30, 25)
)
# Persamaan: sum(w)=1, sum(w*x)=x_target, sum(w*y)=y_target.
M <- rbind(1, abc$x, abc$y)
w_abc <- as.numeric(solve(M, c(1, 2, 2)))
tri_manual <- sum(w_abc * abc$suhu)
knitr::kable(data.frame(abc, bobot = w_abc), digits = 5)
id x y suhu bobot
A 1 1 20 0.44444
B 4 1 30 0.22222
C 2 4 25 0.33333 
tri_manual
## [1] 23.889

5.2.2 Jaringan dan fungsi prediksi

Paket deldir membentuk segitiga. Fungsi berikut menyelesaikan persamaan bobot dan mencari segitiga yang memuat target.

pred_tri <- function(data, baru) {
  jaringan <- deldir::deldir(data$x, data$y, round = FALSE)
  segitiga <- deldir::triang.list(jaringan)
  hasil <- rep(NA_real_, nrow(baru))
  for (j in seq_len(nrow(baru))) {
    for (s in segitiga) {
      M <- rbind(1, s$x, s$y)
      w <- as.numeric(solve(M, c(1, baru$x[j], baru$y[j])))
      if (all(w >= -1e-8) && all(w <= 1 + 1e-8)) {
        hasil[j] <- sum(w * data$suhu[s$ptNum])
        break
      }
    }
  }
  hasil
}

jaringan <- deldir::deldir(stasiun$x, stasiun$y, round = FALSE)
plot(stasiun$x, stasiun$y, asp = 1, pch = 19,
     xlab = "x (km)", ylab = "y (km)", main = "Jaringan segitiga Delaunay")
for (s in deldir::triang.list(jaringan)) {
  polygon(s$x, s$y, border = "#168a9b")
}
points(stasiun$x, stasiun$y, pch = 19, cex = 0.6)
points(target$x, target$y, pch = 8, col = "#d44a30", cex = 1.5)

tri_target <- pred_tri(stasiun, target)
tri_grid <- pred_tri(stasiun, grid)
tri_target
## [1] 27.243
peta(tri_grid, "Interpolasi linier berbasis triangulasi")

# Di luar convex hull:
pred_tri(stasiun, data.frame(x = 12, y = 12))
## [1] NA

5.3 Interpretasi

Berdasarkan interpolasi linier berbasis triangulasi Delaunay, suhu pada lokasi target (5,5) diprediksi sebesar 27,24 °C. Prediksi diperoleh dari suhu di tiga stasiun yang menjadi simpul segitiga tempat target berada, dengan bobot sesuai posisi target di dalam segitiga.

Pada peta, suhu berubah secara linier di dalam setiap segitiga. Nilai suhu tetap menyambung pada batas antarsegi­tiga, tetapi arah atau kemiringan perubahannya dapat berbeda. Hasil NA pada lokasi di luar cakupan titik pengamatan berarti metode tidak menghasilkan prediksi di lokasi tersebut, bukan berarti suhunya nol.

6 Poligon Thiessen (Voronoi)

6.1 Pengantar

Poligon Thiessen membagi wilayah menjadi daerah pengaruh stasiun. Lokasi di dalam sebuah poligon lebih dekat ke stasiun pemiliknya daripada ke stasiun lain. Untuk prediksi titik, hasilnya sama dengan nearest neighbor.

Thiessen juga digunakan untuk rata-rata wilayah:

\[ \overline Z_A=\frac{\sum_i A_i Z_i}{\sum_i A_i}, \]

dengan \(A_i\) luas bagian poligon dalam wilayah kajian. Rata-rata wilayah berbeda dari prediksi satu lokasi.

6.2 Sintaks

6.2.1 Contoh

contoh_thiessen <- data.frame(
  stasiun = c("A", "B", "C"),
  hujan_mm = c(100, 120, 80), luas_km2 = c(40, 35, 25)
)
contoh_thiessen$bobot <- contoh_thiessen$luas_km2 /
  sum(contoh_thiessen$luas_km2)
knitr::kable(contoh_thiessen, digits = 3)
stasiun hujan_mm luas_km2 bobot
A 100 40 0.40
B 120 35 0.35
C 80 25 0.25 
weighted.mean(contoh_thiessen$hujan_mm, contoh_thiessen$luas_km2)
## [1] 102

6.2.2 Poligon simulasi dan rata-rata wilayah

Batas persegi ditentukan melalui rw. Pada wilayah nyata yang tidak persegi, poligon harus dipotong sesuai batas wilayah sebelum luas dihitung.

vor <- deldir::deldir(stasiun$x, stasiun$y,
                     rw = c(0, 10, 0, 10), round = FALSE)
ubin <- deldir::tile.list(vor)
plot(NA, xlim = c(0, 10), ylim = c(0, 10), asp = 1,
     xlab = "x (km)", ylab = "y (km)", main = "Poligon Thiessen")
for (u in ubin) {
  suhu_i <- stasiun$suhu[u$ptNum]
  indeks <- findInterval(suhu_i,
    seq(batas_warna[1], batas_warna[2], length.out = 61),
    all.inside = TRUE)
  polygon(u$x, u$y, col = palet[indeks], border = "white")
}
points(stasiun$x, stasiun$y, pch = 19, cex = 0.6)
points(target$x, target$y, pch = 8, cex = 1.4)
box()

luas <- data.frame(
  indeks = vapply(ubin, function(u) u$ptNum, integer(1)),
  luas_km2 = vapply(ubin, function(u) u$area, numeric(1))
)
luas$id <- stasiun$id[luas$indeks]
luas$suhu <- stasiun$suhu[luas$indeks]
luas$bobot <- luas$luas_km2 / sum(luas$luas_km2)
knitr::kable(head(luas, 8), digits = 3)
indeks luas_km2 id suhu bobot
pt.1 1 1.223 S01 23.756 0.012
pt.2 2 4.565 S02 29.632 0.046
pt.3 3 0.951 S03 27.733 0.010
pt.4 4 2.122 S04 21.683 0.021
pt.5 5 1.916 S05 27.871 0.019
pt.6 6 5.447 S06 27.830 0.054
pt.7 7 5.017 S07 24.199 0.050
pt.8 8 0.514 S08 24.729 0.005 
sum(luas$luas_km2)  # Sekitar 100 km2
## [1] 100
suhu_wilayah <- weighted.mean(luas$suhu, luas$luas_km2)
suhu_wilayah
## [1] 26.402
thiessen_target <- pred_nn(stasiun, target)
thiessen_target
## [1] 27.133
stopifnot(identical(thiessen_target, nn_target))

6.3 Interpretasi

Total luas poligon Thiessen adalah 100 km², sesuai dengan luas wilayah kajian. Rata-rata suhu wilayah berdasarkan pembobotan luas poligon diperoleh sebesar 26,402 °C, atau sekitar 26,40 °C. Stasiun dengan daerah pengaruh lebih luas memperoleh bobot lebih besar dalam perhitungan rata-rata tersebut.

Sementara itu, suhu pada target (5,5) diprediksi sebesar 27,133 °C, atau sekitar 27,13 °C, mengikuti suhu stasiun terdekat. Hasil ini sama dengan metode Nearest Neighbor karena keduanya menetapkan nilai target berdasarkan stasiun terdekat.

Kedua nilai tersebut memiliki makna berbeda: 26,40 °C menggambarkan rata-rata suhu seluruh wilayah, sedangkan 27,13 °C merupakan prediksi suhu pada satu lokasi target.

7 Local averaging

7.1 Pengantar

Local averaging menghitung rata-rata dari \(k\) stasiun terdekat:

\[ \widehat Z(s_0)=\frac{1}{k}\sum_{i\in N_k(s_0)} Z(s_i). \]

Jarak dipakai untuk memilih tetangga, tetapi bobot setiap tetangga terpilih sama, yaitu \(1/k\). Pada IDW, bobot berubah mengikuti jarak. Dengan \(k=1\), local averaging sama dengan NN.

7.2 Sintaks

# Contoh slide 23:
mean(c(80, 100, 70, 90))  # 85 mm
## [1] 85
pred_local <- function(data, baru, k = 4) {
  stopifnot(k >= 1, k <= nrow(data), k == as.integer(k))
  vapply(seq_len(nrow(baru)), function(j) {
    urutan <- order(jarak_ke(data, baru[j, ]))
    mean(data$suhu[urutan[seq_len(k)]])
  }, numeric(1))
}
terpilih <- order(jarak_ke(stasiun, target))[1:4]
knitr::kable(data.frame(
  stasiun[terpilih, ],
  jarak_km = jarak_ke(stasiun, target)[terpilih], bobot = 1/4
), digits = 3)
id x y suhu jarak_km bobot
10 S10 4.833 5.636 27.133 0.658 0.25
25 S25 5.783 4.430 26.924 0.969 0.25
23 S23 3.356 5.432 25.083 1.700 0.25
26 S26 6.645 5.643 26.776 1.767 0.25 
local_target <- pred_local(stasiun, target, k = 4)
local_grid <- pred_local(stasiun, grid, k = 4)
local_target
## [1] 26.479
peta(local_grid, "Local averaging (k = 4)")

knitr::kable(data.frame(
  k = c(1, 4, 8),
  prediksi = vapply(c(1, 4, 8),
    function(k) pred_local(stasiun, target, k), numeric(1))
), digits = 3)
k prediksi
1 27.133
4 26.479
8 26.762

7.3 Interpretasi

Contoh menghasilkan 85 mm. Pada simulasi, empat stasiun terdekat menghasilkan prediksi 26.48 °C.

Nilai \(k\) lebih besar cenderung meratakan variasi lokal. Permukaan dapat berubah mendadak ketika kelompok tetangga terpilih berubah. Nilai \(k=4\) ditetapkan untuk latihan; pemilihan nilai terbaik memerlukan validasi pada data latih.

8 Validasi silang dan evaluasi

8.1 Pengantar

Peta yang menarik belum membuktikan akurasi. Prediksi perlu dibandingkan dengan pengamatan yang tidak digunakan untuk membentuk prediksi tersebut.

  1. Sisihkan sekitar 20% stasiun sebagai data uji akhir.
  2. Pada data latih, jalankan 5-fold cross-validation untuk memperoleh prediksi out-of-fold (OOF). Setiap stasiun diprediksi tanpa menggunakan dirinya dan anggota fold validasinya sebagai data latih.
  3. Gunakan OOF untuk menilai metode dan membentuk ensemble.
  4. Prediksi uji akhir dari seluruh data latih. Respons uji tidak dipakai untuk memilih bobot atau melatih meta-model.

\[ RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\widehat y_i)^2}, \qquad MAE=\frac{1}{n}\sum_{i=1}^{n}|y_i-\widehat y_i|. \]

Keduanya bersatuan °C. Nilai lebih kecil menunjukkan galat lebih rendah pada skema yang sama. RMSE lebih peka terhadap galat besar.

Validasi acak di sini merupakan pengenalan untuk mengisi lokasi yang tersebar di wilayah serupa. Untuk menilai prediksi wilayah baru yang berjauhan, gunakan validasi berbasis blok spasial. Kedekatan lokasi latih dan uji dapat membuat evaluasi acak terlalu optimistis.

8.2 Sintaks

8.2.1 Pembagian data dan fungsi metrik

Thiessen tidak menjadi kolom terpisah karena prediksi titiknya sama dengan NN. Menggabungkan keduanya akan menghitung informasi yang sama dua kali.

set.seed(2026)
indeks_uji <- sample(seq_len(nrow(stasiun)), size = 7)
latih <- stasiun[-indeks_uji, ]
uji <- stasiun[indeks_uji, ]
data.frame(n_latih = nrow(latih), n_uji = nrow(uji))
##   n_latih n_uji
## 1      29     7
pred_semua <- function(data, baru) {
  data.frame(
    IDW = pred_idw(data, baru, p = 2),
    NN = pred_nn(data, baru),
    Triangulasi = pred_tri(data, baru),
    Local = pred_local(data, baru, k = 4)
  )
}
nilai_metrik <- function(y, pred) {
  ok <- is.finite(y) & is.finite(pred)
  n_ok <- sum(ok)
  if (n_ok == 0) return(c(n = 0, RMSE = NA_real_, MAE = NA_real_))
  c(n = n_ok,
    RMSE = sqrt(mean((y[ok] - pred[ok])^2)),
    MAE = mean(abs(y[ok] - pred[ok])))
}
ringkas_metrik <- function(y, pred) {
  out <- do.call(rbind, lapply(pred, function(z) nilai_metrik(y, z)))
  data.frame(metode = rownames(out), out, row.names = NULL)
}

8.2.2 Lima fold pada data latih

set.seed(27)
K <- 5
fold <- sample(rep(seq_len(K), length.out = nrow(latih)))
table(fold)
## fold
## 1 2 3 4 5 
## 6 6 6 6 5
oof <- as.data.frame(matrix(NA_real_, nrow(latih), 4))
names(oof) <- c("IDW", "NN", "Triangulasi", "Local")
for (f in seq_len(K)) {
  valid <- which(fold == f)
  oof[valid, ] <- pred_semua(latih[-valid, ], latih[valid, ])
}
hasil_cv <- ringkas_metrik(latih$suhu, oof)
hasil_cv$cakupan_persen <- 100 * hasil_cv$n / nrow(latih)
knitr::kable(hasil_cv, digits = 3,
  caption = "Galat OOF gabungan dan cakupan lokasi prediksi")
Galat OOF gabungan dan cakupan lokasi prediksi
metode n RMSE MAE cakupan_persen
IDW 29 1.424 1.041 100.000
NN 29 1.053 0.826 100.000
Triangulasi 18 0.634 0.508 62.069
Local 29 1.153 0.928 100.000 
# Perbandingan pada lokasi yang sama untuk seluruh metode
sama_cv <- complete.cases(oof)
hasil_cv_sama <- ringkas_metrik(latih$suhu[sama_cv],
                               oof[sama_cv, , drop = FALSE])
knitr::kable(hasil_cv_sama, digits = 3,
  caption = "Perbandingan OOF pada irisan lokasi yang sama")
Perbandingan OOF pada irisan lokasi yang sama
metode n RMSE MAE
IDW 18 0.668 0.588
NN 18 0.885 0.669
Triangulasi 18 0.634 0.508
Local 18 0.840 0.688 
# RMSE per fold dan rata-ratanya, sesuai pembahasan PDF
rmse_fold <- do.call(rbind, lapply(seq_len(K), function(f) {
  idx <- fold == f
  a <- ringkas_metrik(latih$suhu[idx], oof[idx, , drop = FALSE])
  data.frame(fold = f, a)
}))
knitr::kable(rmse_fold, digits = 3)
fold metode n RMSE MAE
1 IDW 6 0.538 0.458
1 NN 6 0.898 0.681
1 Triangulasi 6 0.464 0.404
1 Local 6 0.457 0.374
2 IDW 6 0.814 0.727
2 NN 6 0.580 0.455
2 Triangulasi 5 0.602 0.385
2 Local 6 0.685 0.629
3 IDW 6 1.681 1.408
3 NN 6 1.046 0.846
3 Triangulasi 3 0.627 0.555
3 Local 6 1.343 1.235
4 IDW 6 2.262 1.762
4 NN 6 1.367 1.142
4 Triangulasi 2 0.537 0.479
4 Local 6 1.720 1.479
5 IDW 5 1.040 0.813
5 NN 5 1.228 1.040
5 Triangulasi 2 1.088 1.083
5 Local 5 1.085 0.921 
aggregate(RMSE ~ metode, data = rmse_fold, FUN = mean)
##        metode    RMSE
## 1         IDW 1.26687
## 2       Local 1.05792
## 3          NN 1.02390
## 4 Triangulasi 0.66364

8.2.3 Prediksi uji akhir

pred_uji <- pred_semua(latih, uji)
hasil_uji <- ringkas_metrik(uji$suhu, pred_uji)
hasil_uji$cakupan_persen <- 100 * hasil_uji$n / nrow(uji)
knitr::kable(hasil_uji, digits = 3,
  caption = "Kinerja uji akhir dan cakupan prediksi")
Kinerja uji akhir dan cakupan prediksi
metode n RMSE MAE cakupan_persen
IDW 7 0.920 0.765 100
NN 7 1.402 1.155 100
Triangulasi 7 1.091 0.858 100
Local 7 0.858 0.746 100 
sama_uji <- complete.cases(pred_uji)
hasil_uji_sama <- ringkas_metrik(uji$suhu[sama_uji],
                               pred_uji[sama_uji, , drop = FALSE])
knitr::kable(hasil_uji_sama, digits = 3,
  caption = "Kinerja uji pada lokasi yang sama untuk seluruh metode")
Kinerja uji pada lokasi yang sama untuk seluruh metode
metode n RMSE MAE
IDW 7 0.920 0.765
NN 7 1.402 1.155
Triangulasi 7 1.091 0.858
Local 7 0.858 0.746

8.3 Interpretasi

Kolom n adalah jumlah lokasi yang berhasil diprediksi. Triangulasi dapat kehilangan cakupan jika titik validasi atau uji berada di luar convex hull data latih. Jangan menyimpulkan metode lebih baik hanya karena RMSE-nya dihitung pada lokasi yang lebih mudah diprediksi.

Tabel irisan lokasi membandingkan metode pada pengamatan yang sama. Jika irisan kecil, kesimpulannya terbatas. Pertimbangkan galat dan cakupan sekaligus.

RMSE OOF gabungan dihitung dari seluruh galat OOF yang tersedia. Rata-rata RMSE fold merata-ratakan RMSE setiap fold. Keduanya tidak harus sama karena ukuran fold dan operasi akar kuadrat berbeda. Jika jumlah fold sama dengan jumlah observasi, prosedurnya disebut leave-one-out cross-validation (LOOCV).

Parameter \(p=2\) dan \(k=4\) ditetapkan sebelum evaluasi. Jika melakukan tuning, pilih parameter menggunakan validasi di dalam data latih; jangan memakai galat uji akhir untuk mengubah parameter.

9 Ensemble interpolasi

9.1 Pengantar

Ensemble menggabungkan prediksi beberapa metode. Modul menggunakan IDW, NN, dan local averaging. Thiessen diwakili NN karena identik untuk titik. Triangulasi dipertahankan sebagai pembanding tersendiri karena tidak selalu menghasilkan prediksi OOF di luar convex hull.

Bentuk Cara penggabungan Sumber aturan penggabungan
Rata-rata sederhana Bobot setiap metode \(1/3\) Ditetapkan langsung
Rata-rata berbobot Bobot proporsional terhadap \(1/RMSE_{OOF}\) OOF data latih
Stacking (pengayaan) Regresi linier sebagai meta-model Prediktor OOF data latih

Untuk rata-rata berbobot,

\[ \widehat Z_{\mathrm{ens}}(s_0)=\sum_{m=1}^{M} w_m\widehat Z_m(s_0), \quad w_m=\frac{1/RMSE_m}{\sum_{\ell=1}^{M}1/RMSE_\ell}. \]

Aturan inverse-RMSE merupakan pilihan sederhana, bukan jaminan bobot optimal. Stacking mempelajari intersep dan koefisien melalui regresi; koefisiennya tidak harus positif atau berjumlah satu. Ensemble tidak dijamin mengungguli metode tunggal.

9.2 Sintaks

9.2.1 Rata-rata sederhana dan berbobot

Bobot ditentukan dari OOF data latih, bukan data uji.

basis <- c("IDW", "NN", "Local")
rmse_basis <- vapply(oof[basis], function(z) {
  sqrt(mean((latih$suhu - z)^2))
}, numeric(1))
stopifnot(all(is.finite(rmse_basis)))
bobot <- 1 / pmax(rmse_basis, 1e-8)
bobot <- bobot / sum(bobot)
knitr::kable(data.frame(metode = basis, RMSE_OOF = rmse_basis,
                       bobot = as.numeric(bobot)), digits = 4)
metode RMSE_OOF bobot
IDW IDW 1.4236 0.2788
NN NN 1.0534 0.3768
Local Local 1.1525 0.3444 
sum(bobot)
## [1] 1
uji_rata <- rowMeans(pred_uji[basis])
uji_bobot <- as.numeric(as.matrix(pred_uji[basis]) %*% bobot)

9.2.2 Stacking dengan meta-model regresi

# Tiap baris prediktor berikut merupakan prediksi OOF.
data_meta <- data.frame(suhu = latih$suhu, oof[basis])
meta <- lm(suhu ~ IDW + NN + Local, data = data_meta)
summary(meta)
## 
## Call:
## lm(formula = suhu ~ IDW + NN + Local, data = data_meta)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.820 -0.425  0.079  0.693  1.399 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -1.312      5.959   -0.22    0.828  
## IDW           -0.665      0.416   -1.60    0.123  
## NN             0.812      0.329    2.47    0.021 *
## Local          0.898      0.425    2.11    0.045 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.01 on 25 degrees of freedom
## Multiple R-squared:  0.756,  Adjusted R-squared:  0.727 
## F-statistic: 25.8 on 3 and 25 DF,  p-value: 7.92e-08
# Base-model memprediksi uji dari seluruh data latih.
uji_stack <- as.numeric(predict(meta, newdata = pred_uji[basis]))

evaluasi_ensemble <- data.frame(
  IDW = pred_uji$IDW, NN = pred_uji$NN, Local = pred_uji$Local,
  Rata_rata = uji_rata, Berbobot = uji_bobot, Stacking = uji_stack
)
hasil_akhir <- ringkas_metrik(uji$suhu, evaluasi_ensemble)
hasil_akhir <- hasil_akhir[order(hasil_akhir$RMSE), ]
knitr::kable(hasil_akhir, digits = 3,
  caption = "Metode dasar dan ensemble pada tujuh lokasi uji yang sama")
Metode dasar dan ensemble pada tujuh lokasi uji yang sama
metode n RMSE MAE
3 Local 7 0.858 0.746
1 IDW 7 0.920 0.765
4 Rata_rata 7 1.038 0.888
5 Berbobot 7 1.058 0.905
6 Stacking 7 1.299 1.053
2 NN 7 1.402 1.155 
plot(uji$suhu, uji_bobot, pch = 21, bg = "#1495a7",
     xlab = "Suhu pengamatan (C)", ylab = "Prediksi ensemble berbobot (C)",
     xlim = range(c(uji$suhu, uji_bobot)),
     ylim = range(c(uji$suhu, uji_bobot)), asp = 1,
     main = "Pengamatan vs prediksi uji")
abline(0, 1, col = "#d44a30", lty = 2, lwd = 2)

9.2.3 Peta dari model yang dievaluasi

Peta ini memakai data latih dan bobot yang sudah ditetapkan. Data uji ditandai agar lokasi evaluasi terlihat.

basis_grid <- data.frame(
  IDW = pred_idw(latih, grid),
  NN = pred_nn(latih, grid),
  Local = pred_local(latih, grid, k = 4)
)
grid_bobot <- as.numeric(as.matrix(basis_grid) %*% bobot)
peta(grid_bobot, "Ensemble berbobot dari data latih")
points(uji$x, uji$y, pch = 4, col = "#194f98", lwd = 2, cex = 1.3)
legend("topleft", c("Stasiun", "Data uji", "Target"),
       pch = c(21, 4, 8), pt.bg = "white",
       col = c("black", "#194f98", "black"), bty = "n")

9.3 Interpretasi

Metode dengan RMSE OOF lebih kecil mendapat bobot lebih besar menurut aturan yang dipakai. Bobot adalah kontribusi dalam rumus ensemble, bukan probabilitas bahwa suatu metode benar.

Pada tujuh lokasi uji simulasi, metode dengan RMSE terkecil dalam tabel adalah Local, dengan RMSE 0.858 °C. Hasil ini berlaku pada simulasi dan pembagian tersebut, bukan semua data spasial.

Titik dekat garis merah menunjukkan prediksi mendekati pengamatan. Titik di atas garis menunjukkan prediksi terlalu tinggi; di bawah garis menunjukkan prediksi terlalu rendah.

Koefisien, p-value, dan \(R^2\) pada summary(meta) menjelaskan penyesuaian meta-model pada data latih OOF, bukan pengganti evaluasi uji akhir. Galat ensemble pada baris OOF yang juga dipakai untuk menentukan bobot atau melatih meta-model tidak digunakan sebagai estimasi kinerja akhir.

10 Perbandingan pada satu target

10.1 Pengantar

Tabel berikut kembali memakai seluruh stasiun, seperti demonstrasi awal. Ensemble tidak dicampurkan karena contoh evaluasinya memakai data latih.

10.2 Sintaks

tabel_target <- data.frame(
  metode = c("IDW (p=2)", "Nearest neighbor", "Triangulasi linier",
             "Thiessen (prediksi titik)", "Local averaging (k=4)"),
  suhu_prediksi = c(idw_target, nn_target, tri_target,
                   thiessen_target, local_target)
)
knitr::kable(tabel_target, digits = 3,
  caption = "Prediksi target (5,5) dari semua stasiun")
Prediksi target (5,5) dari semua stasiun
metode suhu_prediksi
IDW (p=2) 26.801
Nearest neighbor 27.133
Triangulasi linier 27.243
Thiessen (prediksi titik) 27.133
Local averaging (k=4) 26.479

10.3 Interpretasi

Perbedaan hasil berasal dari aturan pemilihan stasiun dan pembobotan. NN dan Thiessen menghasilkan angka sama. Target tidak memiliki pengamatan, sehingga tabel ini tidak dapat menentukan metode paling akurat.

11 Latihan mahasiswa

  1. Ganti target menjadi (4,6). Identifikasi tetangga terdekat dan bandingkan hasil semua metode.
  2. Gunakan IDW dengan \(p=1,2,3\). Jelaskan perubahan bobot dan pola peta.
  3. Gunakan local averaging dengan \(k=2,4,8\). Bandingkan pengaruh jumlah tetangga.
  4. Ulangi hitungan triangulasi PDF tanpa pembulatan antara. Mengapa pembulatan memengaruhi hasil?
  5. Mengapa prediksi titik Thiessen sama dengan NN, tetapi rata-rata wilayah Thiessen berbeda dari rata-rata sederhana seluruh stasiun?
  6. Bandingkan hasil validasi pada seluruh lokasi tersedia dan irisan lokasi. Mengapa cakupan perlu dilaporkan?
  7. Apakah ensemble selalu memperbaiki prediksi dibandingkan metode tunggal? Gunakan tabel uji sebagai bukti.
  8. Mengapa meta-model stacking tidak boleh dilatih dari prediksi base-model pada baris yang juga digunakan melatih base-model tersebut?
  9. Pengayaan: rancang pembagian blok spasial. Jelaskan mengapa hasilnya dapat berbeda dari validasi acak.

12 Ringkasan dan rujukan

Metode Informasi yang digunakan Ciri keluaran
IDW Nilai dan bobot jarak Titik dekat lebih berpengaruh
NN Satu stasiun terdekat Konstan dalam daerah pengaruh
Triangulasi linier Tiga simpul segitiga Linier per segitiga; NA di luar convex hull
Thiessen Daerah pengaruh terdekat Setara NN untuk titik; luas untuk rata-rata wilayah
Local averaging \(k\) tetangga terdekat Rata-rata berbobot sama
Ensemble Prediksi beberapa metode Bergantung pada komponen dan penggabungannya

Penyesuaian terhadap PDF: hitungan IDW dilengkapi; hasil triangulasi dihitung tanpa pembulatan antara (slide 15); nilai B pada contoh Thiessen diluruskan (slide 19–20); RMSE ditulis eksplisit menggunakan akar kuadrat. Validasi ditempatkan sebelum ensemble agar sumber bobot OOF dapat diikuti.

Rujukan

  • Djuraidah, A. Tanpa tahun. K9A. Interpolasi Titik. Materi STA1352 Regresi Spasial, IPB University. PDF yang disediakan untuk penyusunan modul.
  • Dokumentasi R paket deldir: fungsi deldir, triang.list, dan tile.list. Dokumentasi lokal dapat dibuka dengan kode di bawah.
help("deldir", package = "deldir")
help("triang.list", package = "deldir")
help("tile.list", package = "deldir")
help("html_document", package = "rmarkdown")
# Versi R dan paket untuk menelusuri hasil praktikum
sessionInfo()
## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## 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     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     R6_2.6.1          fastmap_1.2.0     xfun_0.52        
##  [5] cachem_1.1.0      knitr_1.50        deldir_2.0-4      htmltools_0.5.8.1
##  [9] rmarkdown_2.29    lifecycle_1.0.5   cli_3.6.5         sass_0.4.10      
## [13] jquerylib_0.1.4   compiler_4.5.3    rstudioapi_0.19.0 tools_4.5.3      
## [17] evaluate_1.0.5    bslib_0.9.0       yaml_2.3.10       rlang_1.3.0      
## [21] jsonlite_2.0.0