Interpolasi Spasial (Data Harga Rumah Boston)

Data

Data yang digunakan adalah data harga rumah Boston dari paket spData (boston). Data berisi 506 census tract beserta koordinat (LON, LAT). Variabel yang diinterpolasi adalah CMEDV, yaitu nilai median rumah yang ditempati pemilik (dalam ribu dolar AS) dalam versi yang sudah dikoreksi. Variabel MEDV tidak dipakai karena nilainya terpotong di 50 pada sejumlah tract.

Langkah-langkah:

  1. Memuat data boston dan mengubahnya menjadi objek sf.
  2. Membuat variabel vble untuk menyimpan nilai CMEDV.
  3. Membuat batas wilayah (map) dari shapefile tract Boston yang disertakan dalam spData.
# Memuat pustaka yang diperlukan
library(dismo)
library(sf)
library(tmap)
library(viridis)
library(gstat)
library(ggplot2)
library(terra)  

sf_use_s2(FALSE)
set.seed(123)

# Memuat data harga rumah Boston
data(boston, package = "spData")
b <- boston.c

# Hapus koordinat kembar (bila ada)
b <- b[!duplicated(b[, c("LON", "LAT")]), ]

# Membuat objek sf, lalu proyeksikan ke UTM 19N (satuan meter)
d <- st_as_sf(b, coords = c("LON", "LAT"), crs = 4326)
d <- st_transform(d, 32619)

# Menambahkan kolom nilai median rumah (ribu dolar) sebagai variabel yang diinterpolasi
d$vble <- d$CMEDV
d <- d[, "vble"]

nrow(d)
## [1] 506
summary(d$vble)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   17.02   21.20   22.53   25.00   50.00

Berdasarkan hasil pemeriksaan, terdapat 506 observasi dengan nilai CMEDV berkisar antara 5 hingga 50 ribu dolar AS. Nilai rata-rata sebesar 22,53 ribu dolar dan median sebesar 21,20 ribu dolar menunjukkan bahwa harga rumah pada sebagian besar wilayah Boston berada di sekitar rentang tersebut.

map <- tryCatch(
  st_sf(geometry = st_concave_hull(st_union(d), ratio = 0.25)),
  error = function(e) NULL)

# Cadangan: gabungan lingkaran 2,5 km di sekitar tiap titik
if (is.null(map)) {
  map <- st_sf(geometry = st_union(st_buffer(d, 2500)))
}

tmap_mode("plot")
tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
  tm_dots("vble", palette = "viridis", size = 0.4)

Di sini, map adalah objek yang memuat batas wilayah Boston, yang nantinya akan digunakan untuk membatasi area perhitungan.

Grid raster dengan resolusi 500 meter dibuat untuk menentukan lokasi-lokasi yang akan menjadi titik prediksi. Titik pusat setiap sel grid kemudian digunakan sebagai lokasi interpolasi dan dibatasi sesuai wilayah Boston.

# Membuat grid raster resolusi 500 m
grid <- rast(vect(map), resolution = 500)

# Mengambil koordinat titik pusat dari setiap sel dalam grid
xy <- terra::xyFromCell(grid, 1:terra::ncell(grid))

# Membuat objek sf untuk lokasi prediksi
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))

# Menyaring lokasi prediksi yang berada dalam batas Boston
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deterministik

Closest Observation

Metode Closest Observation menggunakan nilai dari lokasi pengamatan terdekat untuk memperkirakan harga rumah pada lokasi yang belum teramati. Pembagian wilayah prediksi ditampilkan melalui diagram Voronoi, dengan setiap poligon mewakili pengaruh satu titik pengamatan terdekat.

# Geser titik sangat kecil (maks. 10 m) hanya agar Voronoi tidak gagal
dj <- st_jitter(d, amount = 10)

v <- terra::voronoi(x = vect(dj), bnd = vect(map))
terra::plot(v)
terra::points(vect(dj), cex = 0.5)

# Untuk chunk berikutnya: poligon + nilai tiap titik (tanpa dipotong)
v_sf <- st_sf(geometry = st_geometry(st_as_sf(v)))
st_crs(v_sf) <- st_crs(dj)
v_sf <- st_join(v_sf, dj["vble"])
# Prediction
tm_shape(v_sf) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)

pred_closest <- function(train, newdata) {
  train$vble[st_nearest_feature(newdata, train)]
}

resp <- coop
resp$pred <- pred_closest(d, coop)
pred <- terra::rasterize(vect(resp), grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c1

Peta hasil interpolasi memperlihatkan persebaran nilai prediksi berdasarkan titik pengamatan terdekat. Perubahan nilai antarlokasi terlihat mengikuti batas poligon Voronoi sehingga menghasilkan pola wilayah yang cenderung diskrit.

Inverse Distance Weighting

Interpolasi IDW dilakukan dengan memberikan bobot berdasarkan jarak antara lokasi prediksi dan seluruh titik pengamatan. Lokasi yang lebih dekat memperoleh bobot lebih besar sehingga memiliki pengaruh lebih kuat terhadap hasil prediksi.

# Menerapkan IDW menggunakan semua lokasi pengamatan
res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))

# Mendapatkan prediksi untuk lokasi yang belum teramati
resp <- predict(res, coop, debug.level = 0)

resp$x <- st_coordinates(resp)[, 1]
resp$y <- st_coordinates(resp)[, 2]
resp$pred <- resp$var1.pred

# Membuat peta prediksi dengan metode IDW
pred <- terra::rasterize(vect(resp), grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c2

Peta hasil interpolasi IDW menunjukkan persebaran harga rumah yang terbentuk dari kombinasi nilai titik-titik pengamatan. Pengaruh setiap titik berkurang seiring bertambahnya jarak sehingga menghasilkan permukaan prediksi yang lebih kontinu dibandingkan metode Closest Observation.

Nearest Neighbors

Metode Nearest Neighbors digunakan untuk memprediksi harga rumah berdasarkan rata-rata nilai dari lima titik pengamatan terdekat. Setiap titik tetangga diberikan bobot yang sama dalam menentukan nilai prediksi.

res <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))

resp <- predict(res, coop, debug.level = 0)

resp$x <- st_coordinates(resp)[, 1]
resp$y <- st_coordinates(resp)[, 2]
resp$pred <- resp$var1.pred

pred <- terra::rasterize(vect(resp), grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c3

Hasil interpolasi menunjukkan pola persebaran harga rumah berdasarkan rata-rata lima lokasi terdekat. Penggunaan beberapa tetangga membuat prediksi tidak hanya bergantung pada satu titik pengamatan, sehingga pola yang dihasilkan cenderung lebih merata.

Pendekatan Ensembel

Pendekatan Ensemble menggabungkan hasil prediksi Closest Observation, IDW, dan Nearest Neighbors dengan bobot yang sama, yaitu masing-masing sebesar sepertiga. Penggabungan ini dilakukan untuk memperoleh hasil prediksi yang memanfaatkan karakteristik dari ketiga metode tersebut.

# Closest observation
p1 <- pred_closest(d, coop)

# IDW
gs <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))
p2 <- predict(gs, coop, debug.level = 0)$var1.pred

# Nearest neighbors
nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
p3 <- predict(nn, coop, debug.level = 0)$var1.pred

# Ensemble (equal weights)
weights <- c(1/3, 1/3, 1/3)
p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]

resp <- data.frame(x = st_coordinates(coop)[, 1],
                   y = st_coordinates(coop)[, 2],
                   pred = p4)
resp <- st_as_sf(resp, coords = c("x", "y"), crs = st_crs(map))
pred <- terra::rasterize(vect(resp), grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c4

# Gabungkan semua plot dalam satu frame
tmap_arrange(c1, c2, c3, c4, ncol = 2)

Keempat peta memperlihatkan perbedaan pola hasil interpolasi. Closest Observation menghasilkan batas wilayah prediksi yang tegas, IDW memberikan pola yang lebih kontinu, sedangkan Nearest Neighbors menghasilkan prediksi berdasarkan rata-rata lima titik terdekat. Ensemble menggabungkan ketiga pendekatan tersebut sehingga menghasilkan pola prediksi yang merupakan kombinasi dari masing-masing metode.

Evaluasi Model

set.seed(123)

# Fungsi untuk menghitung RMSE
RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

# Pembagian data untuk K-fold
kf <- dismo::kfold(nrow(d), k = 5)

# Vectors to store the RMSE values obtained with each method
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Ensemble

for (k in 1:5) {
  # Split data into test and train
  test  <- d[kf == k, ]
  train <- d[kf != k, ]

  # Closest observation
  p1 <- pred_closest(train, test)
  rmse1[k] <- RMSE(test$vble, p1)

  # IDW
  gs <- gstat(formula = vble ~ 1, locations = train,
              nmax = nrow(train), set = list(idp = 1))
  p2 <- predict(gs, test, debug.level = 0)$var1.pred
  rmse2[k] <- RMSE(test$vble, p2)

  # Nearest neighbors
  nn <- gstat(formula = vble ~ 1, locations = train,
              nmax = 5, set = list(idp = 0))
  p3 <- predict(nn, test, debug.level = 0)$var1.pred
  rmse3[k] <- RMSE(test$vble, p3)

  # Ensemble (weights are inverse RMSE so lower RMSE higher weight)
  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k])
  weights <- w / sum(w)
  p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
  rmse4[k] <- RMSE(test$vble, p4)
}

# RMSE obtained for each of the 5 splits
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, ensemble = rmse4)
##   closest.obs      IDW nearest.neigh ensemble
## 1    7.202166 9.644621      6.974510 6.837791
## 2    5.402557 8.968279      5.823436 5.586777
## 3    5.210359 7.547950      5.282393 5.076404
## 4    6.936230 7.340753      5.369727 5.512547
## 5    5.575956 6.813363      4.477856 4.477634
# Average RMSE over the 5 splits
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), ensemble = mean(rmse4))
##   closest.obs      IDW nearest.neigh ensemble
## 1    6.065454 8.062993      5.585584 5.498231

Berdasarkan rata-rata RMSE dari lima fold, metode Ensemble menghasilkan nilai kesalahan terkecil sebesar 5,498 ribu dolar AS. Selanjutnya, Nearest Neighbors menghasilkan RMSE sebesar 5,586 ribu dolar, Closest Observation sebesar 6,065 ribu dolar, dan IDW sebesar 8,063 ribu dolar. Dengan demikian, Ensemble menjadi metode dengan kinerja prediksi terbaik di antara keempat metode tersebut.

Metode Stokastik

Kriging

Ordinary Kriging merupakan metode interpolasi yang mempertimbangkan hubungan spasial antarobservasi melalui variogram. Sebelum melakukan prediksi, dilakukan eksplorasi pola semivariogram dan penyesuaian model untuk memperoleh gambaran struktur ketergantungan spasial pada data.

library(sp)
library(mapview)

# Peta interaktif titik pengamatan (setara mapview(meuse, zcol = "zinc"))
mapview(d, zcol = "vble", map.types = "CartoDB.Voyager")
set.seed(123)
mapview(coop[sample(nrow(coop), min(3000, nrow(coop))), ],
        map.types = "CartoDB.Voyager")
# Variogram dengan cloud plot
vc <- variogram(vble ~ 1, d, cloud = TRUE)
plot(vc)

Variogram cloud memperlihatkan variasi semivarian pada berbagai jarak antarobservasi. Pola sebaran tersebut memberikan gambaran awal mengenai hubungan spasial antar lokasi pengamatan.

# Variogram dengan plot standar
v <- variogram(vble ~ 1, data = d)
plot(v)

Variogram menunjukkan perubahan semivarian terhadap jarak antarobservasi. Pola kenaikan semivarian mengindikasikan adanya ketergantungan spasial, sedangkan kecenderungan mendatar menunjukkan bahwa hubungan spasial semakin berkurang pada jarak tertentu.

# Menampilkan daftar model variogram
vgm()
##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)
show.vgms(par.strip.text = list(cex = 0.75))

Model variogram awal ditentukan menggunakan model spherical dengan nilai sill sebesar 70% varians data, nugget sebesar 30% varians, dan range sebesar sepertiga panjang diagonal wilayah pengamatan.

s2 <- var(d$vble)
bb <- st_bbox(d)
diag_len <- as.numeric(sqrt((bb["xmax"] - bb["xmin"])^2 + (bb["ymax"] - bb["ymin"])^2))

# Menentukan model variogram awal dengan parameter psill, model, range, dan nugget
vinitial <- vgm(psill = s2 * 0.7, model = "Sph", range = diag_len / 3, nugget = s2 * 0.3)

# Memvisualisasikan variogram dengan model awal
plot(v, vinitial, cex = 1.5)

# Menyesuaikan variogram dengan model "Sph"
fv <- fit.variogram(object = v, model = vinitial)
fv
##   model    psill   range
## 1   Nug  0.00000    0.00
## 2   Sph 93.82377 2157.72
plot(v, fv, cex = 1.5)

Hasil fitting menghasilkan model spherical dengan nilai partial sill sebesar 93,82377 dan range sebesar 2.157,72 meter. Nilai nugget sebesar 0 menunjukkan bahwa model hasil fitting tidak mengestimasi adanya komponen nugget.

library(automap)
# Autofit variogram (membutuhkan objek Spatial)
afvg <- tryCatch(autofitVariogram(vble ~ 1, as(d, "Spatial")),
                 error = function(e) { message("autofit dilewati: ", e$message); NULL })
if (!is.null(afvg)) plot(afvg)

# Membuat objek gstat dengan model variogram hasil fitting (fv)
k <- gstat(formula = vble ~ 1, data = d, model = fv, nmax = 50)

# Melakukan prediksi pada grid menggunakan model gstat yang sudah dibuat
kpred <- predict(k, coop, debug.level = 0)

# Visualisasi prediksi dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d, size = 0.4) +
  scale_color_viridis(name = "CMEDV (ribu $)") + theme_bw()

# Visualisasi prediksi varians dengan warna berdasarkan var1.var
ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
  geom_sf(data = d, size = 0.4) +
  scale_color_viridis(name = "variance") + theme_bw()

Peta hasil Ordinary Kriging memperlihatkan estimasi harga rumah pada lokasi yang tidak memiliki pengamatan langsung. Nilai prediksi ditampilkan menggunakan gradasi warna, sedangkan titik-titik pengamatan digunakan sebagai pembanding lokasi data aktual.

Perbandingan Seluruh Metode (Cross-Validation)

fit_vgm <- function(train) {
  vemp <- variogram(vble ~ 1, data = train)
  bb <- st_bbox(train)
  dl <- as.numeric(sqrt((bb["xmax"] - bb["xmin"])^2 + (bb["ymax"] - bb["ymin"])^2))
  s2 <- var(train$vble)
  best <- NULL; best_sse <- Inf
  for (m in c("Sph", "Exp", "Gau")) {
    fit <- tryCatch(
      suppressWarnings(fit.variogram(
        vemp, vgm(psill = s2 * 0.7, model = m, range = dl / 3, nugget = s2 * 0.3))),
      error = function(e) NULL)
    if (!is.null(fit) && !isTRUE(attr(fit, "singular"))) {
      sse <- attr(fit, "SSErr")
      if (!is.null(sse) && sse < best_sse) { best <- fit; best_sse <- sse }
    }
  }
  if (is.null(best)) best <- vgm(psill = s2, model = "Nug")
  best
}

rmse5 <- rep(NA, 5) # Kriging
for (k in 1:5) {
  test  <- d[kf == k, ]
  train <- d[kf != k, ]
  g <- gstat(formula = vble ~ 1, data = train, model = fit_vgm(train), nmax = 50)
  p5 <- predict(g, test, debug.level = 0)$var1.pred
  rmse5[k] <- RMSE(test$vble, p5)
}

hasil_cv <- data.frame(closest.obs = rmse1, IDW = rmse2, nearest.neigh = rmse3,
                       ensemble = rmse4, kriging = rmse5)
hasil_cv
##   closest.obs      IDW nearest.neigh ensemble   kriging
## 1    7.202166 9.644621      6.974510 6.837791 36.617380
## 2    5.402557 8.968279      5.823436 5.586777  5.282555
## 3    5.210359 7.547950      5.282393 5.076404 31.677523
## 4    6.936230 7.340753      5.369727 5.512547 15.504160
## 5    5.575956 6.813363      4.477856 4.477634 23.812472
# Rata-rata RMSE
rata <- colMeans(hasil_cv)
rata
##   closest.obs           IDW nearest.neigh      ensemble       kriging 
##      6.065454      8.062993      5.585584      5.498231     22.578818
hasil <- data.frame(metode = names(rata), rmse = as.numeric(rata))
ggplot(hasil, aes(x = reorder(metode, rmse), y = rmse)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(rmse, 3)), vjust = -0.4) +
  labs(x = "Metode", y = "Rata-rata RMSE (5-fold CV, ribu $)",
       title = "Perbandingan RMSE Metode Interpolasi") +
  theme_bw()

Berdasarkan hasil evaluasi menggunakan 5-fold cross-validation, metode Ensemble menghasilkan kinerja interpolasi terbaik dengan rata-rata RMSE terendah sebesar 5,498 ribu dolar AS, diikuti Nearest Neighbors (5,586), Closest Observation (6,065), IDW (8,063), dan Kriging (22,579). Nilai RMSE Kriging juga menunjukkan variasi yang cukup besar pada setiap fold, yaitu berkisar antara 5,283 hingga 36,617 ribu dolar AS. Hasil ini menunjukkan bahwa penggabungan prediksi dari Closest Observation, IDW, dan Nearest Neighbors melalui metode Ensemble mampu menghasilkan kesalahan prediksi paling rendah pada data yang digunakan, sedangkan Kriging memberikan kinerja paling rendah.