Data yang Dgunakan

Dalam bagian ini, kita akan menggunakan data harga rumah dari paket spData (Bivand et al., 2022) yang berisi informasi tentang penjualan rumah di Lucas County, Ohio, Amerika Serikat, pada periode 1993–1998. Data yang digunakan mencakup beberapa kolom seperti harga jual rumah (price) dan geometri yang memuat koordinat spasial dari masing-masing rumah.

Packages

library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(sf)
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(terra)
## terra 1.9.34
library(tmap)
library(viridis)
## Loading required package: viridisLite
library(dismo)
## Loading required package: raster
## Loading required package: sp
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gstat)
library(ggplot2)
library(sp)
library(mapview)
library(automap)

Dataset

# Memuat dataset house dari package spData
data(house)

# Mengambil sampel 1000 observasi secara acak
set.seed(123)
idx <- sample(1:nrow(house), 1000)

# Mengubah subset data menjadi objek sf
d <- st_as_sf(house[idx, ])

# Menambahkan kolom 'vble' sebagai variabel target interpolasi 
d$vble <- d$price

# Membuat batas wilayah berdasarkan sebaran titik sampel house
map <- st_union(d) %>% st_convex_hull() %>% st_sf()

# Menampilkan peta sebaran harga properti / rumah
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(map) + 
  tm_polygons(alpha = 0.1, border.col = "grey") + 
  tm_shape(d) + 
  tm_dots("vble", palette = "viridis", size = 0.04) +
  tm_layout(title = "Peta Sebaran Harga Sampel Rumah Lucas County")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

Membuat grid prediksi

# Membuat grid raster 100x100 sel untuk menutupi wilayah studi Lucas County
grid <- rast(map, nrows = 100, ncols = 100)

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

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

# Menyaring lokasi prediksi agar hanya berada di dalam batas wilayah studi (map)
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deterministik

Closest Observation

# Membuat diagram Voronoi dari lokasi sampel rumah 
v <- voronoi(x = vect(d), bnd = map)

# Menampilkan plot dasar diagram Voronoi
plot(v)

# Menambahkan titik lokasi rumah asli di atas diagram Voronoi
points(vect(d), cex = 0.5, col = "red")

Prediction & Visualisasi Peta Voronoi

# Mengonversi objek Voronoi (v) menjadi objek sf
v_sf <- st_as_sf(v)

# Menampilkan peta poligon Voronoi berdasarkan harga rumah (vble)
tm_shape(v_sf) + 
  tm_fill(col = "vble", palette = "viridis", alpha = 0.6, title = "Harga Rumah") +
  tm_borders(col = "grey80", alpha = 0.5) +
  tm_layout(title = "Prediksi Closest Observation (Voronoi)")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Mengonversi objek Voronoi ke format sf
v_sf <- st_as_sf(v)

# Menampilkan peta poligon Voronoi berdasarkan harga rumah 
tm_shape(v_sf) + 
  tm_fill(col = "vble", palette = "viridis", alpha = 0.6, title = "Prediksi Voronoi") +
  tm_shape(d) + tm_dots(size = 0.02, col = "black") +
  tm_layout(title = "Diagram Voronoi & Titik Sampel Rumah")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Melakukan prediksi ke titik grid  
resp_vor <- st_intersection(v_sf, coop)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
resp_vor$pred <- resp_vor$vble
pred_vor_rast <- rasterize(resp_vor, grid, field = "pred", fun = "mean")

# Menyimpan plot untuk perbandingan ensembel nantinya
c1 <- tm_shape(pred_vor_rast) + tm_raster(alpha = 0.6, palette = "viridis") + tm_layout(title = "Closest Observation")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).[v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
c1

Invers Distance Weight (IDW)

# Menerapkan model IDW menggunakan seluruh titik sampel rumah (d)
res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))

# Mendapatkan prediksi harga untuk lokasi grid yang belum teramati 
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
# Menyesuaikan kolom koordinat dan hasil prediksi untuk proses rasterisasi
resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

# Membuat peta prediksi harga dengan metode IDW
pred_idw <- rasterize(resp, grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred_idw) + 
  tm_raster(alpha = 0.6, palette = "viridis", title = "Harga ($)") +
  tm_layout(title = "IDW (idp = 1)")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
# Menampilkan peta c2
c2

## Nearest Neighbors

# Menerapkan Nearest Neighbors menggunakan 5 tetangga terdekat (nmax = 5) dan bobot sama (idp = 0)
res_nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))

# Memprediksi harga rumah di lokasi grid yang belum teramati
resp_nn <- predict(res_nn, coop)
## [inverse distance weighted interpolation]
# Menyesuaikan koordinat dan hasil prediksi untuk proses rasterisasi
resp_nn$x <- st_coordinates(resp_nn)[,1]
resp_nn$y <- st_coordinates(resp_nn)[,2]
resp_nn$pred <- resp_nn$var1.pred

# Mengubah titik hasil prediksi menjadi objek raster grid
pred_nn <- rasterize(resp_nn, grid, field = "pred", fun = "mean")

# Membuat peta prediksi harga dengan metode Nearest Neighbors dan menyimpannya ke c3
c3 <- tm_shape(pred_nn) + 
  tm_raster(alpha = 0.6, palette = "viridis", title = "Harga ($)") +
  tm_layout(title = "Nearest Neighbors (k=5)")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
# Menampilkan peta c3
c3

## Pendekatan Ensembel

# Closest Observation (Voronoi) 
v <- voronoi(x = vect(d), bnd = map)
v <- st_as_sf(v)
p1 <- st_intersection(v, coop)$vble
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# IDW 
gs <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))
p2 <- predict(gs, coop)$var1.pred
## [inverse distance weighted interpolation]
# Nearest Neighbors (k = 5) 
nn <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
p3 <- predict(nn, coop)$var1.pred
## [inverse distance weighted interpolation]
# Ensemble (Equal Weights: Bobot 1/3 untuk masing-masing metode)
weights <- c(1/3, 1/3, 1/3)
p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
# Menyusun kembali ke dalam data frame dan objek sf
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))

# Rasterisasi hasil ensembel dan membuat plot c4
pred_ens <- rasterize(resp, grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred_ens) + 
  tm_raster(alpha = 0.6, palette = "viridis", title = "Harga ($)") +
  tm_layout(title = "Ensemble Method")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
# Menampilkan peta ensembel (c4)
c4

# Menggabungkansemua plot dalam satu frame
tmap_arrange(c1, c2, c3, c4, ncol = 2)
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

# Evaluasi Model

set.seed(123)

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

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

# Vectors untuk menyimpan nilai RMSE masing-masing metode (+ Kriging)
rmse1 <- rep(NA, 5) # Closest observation (Voronoi)
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors (k=5)
rmse4 <- rep(NA, 5) # Ensemble (berbasis bobot invers RMSE)
rmse5 <- rep(NA, 5) # Ordinary Kriging

for(k in 1:5){
  # Split data menjadi testing dan training
  test  <- d[kf == k, ]
  train <- d[kf != k, ]
  
  # Konversi ke Spatial object untuk kebutuhan gstat variogram/kriging
  test_sp  <- as(test, "Spatial")
  train_sp <- as(train, "Spatial")

  # Closest Observation (Voronoi)
  v <- voronoi(x = vect(train), bnd = map)
  v <- st_as_sf(v)
  p1 <- st_intersection(v, test)$vble
  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)$var1.pred
  rmse2[k] <- RMSE(test$vble, p2)

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

  # Ensemble 
  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)

  # Ordinary Kriging
  try({
    v_var <- variogram(log(vble) ~ 1, data = train_sp)
    fit_v <- fit.variogram(v_var, model = vgm(model = "Sph"))
    k_mod <- gstat(formula = log(vble) ~ 1, data = train_sp, model = fit_v)
    k_pred <- predict(k_mod, test_sp)
    p5 <- exp(k_pred$var1.pred) # Back-transformation dari log ke harga asli ($)
    rmse5[k] <- RMSE(test$vble, p5)
  }, silent = TRUE)
}
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

Menampilkan Hasil RMSE

# Tabel RMSE untuk masing-masing fold (1 s.d. 5)
df_rmse_folds <- data.frame(
  Closest_Obs = rmse1,
  IDW = rmse2,
  Nearest_Neigh = rmse3,
  Ensemble = rmse4,
  Kriging = rmse5
)
print(" RMSE per Fold ")
## [1] " RMSE per Fold "
print(df_rmse_folds)
##   Closest_Obs      IDW Nearest_Neigh Ensemble  Kriging
## 1    49372.76 53183.18      46089.06 41646.79 38715.07
## 2    34011.64 34974.97      33786.93 28534.51 28591.93
## 3    34553.57 50871.90      41513.72 36590.89 36612.29
## 4    49888.53 60335.52      46986.69 44140.44 40441.00
## 5    31388.49 33249.47      30058.58 26717.89 29113.38
# Tabel Rata-rata RMSE 
df_rmse_avg <- data.frame(
  Closest_Obs = mean(rmse1, na.rm = TRUE),
  IDW = mean(rmse2, na.rm = TRUE),
  Nearest_Neigh = mean(rmse3, na.rm = TRUE),
  Ensemble = mean(rmse4, na.rm = TRUE),
  Kriging = mean(rmse5, na.rm = TRUE)
)
print(" Rata-rata RMSE (5-Fold Cross Validation) ")
## [1] " Rata-rata RMSE (5-Fold Cross Validation) "
print(df_rmse_avg)
##   Closest_Obs      IDW Nearest_Neigh Ensemble  Kriging
## 1       39843 46523.01         39687 35526.11 34694.73

Metode Stokastik

Kriging

# Mengonversi dataset rumah sampel d menjadi objek 'sf' (simple feature) dan grid ke format sf
d_sf <- st_as_sf(d)
coop_sf <- st_as_sf(coop)

# Menampilkan peta interaktif 
mapview(d_sf, zcol = "vble", map.types = "CartoDB.Voyager")
mapview(coop_sf, map.types = "CartoDB.Voyager")
## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)
# Mengonversi objek sf ke Spatial object untuk gstat
d_sp <- as(d, "Spatial")
coop_sp <- as(coop, "Spatial")
# Variogram dengan cloud plot
vc <- variogram(log(vble) ~ 1, d_sp, cloud = TRUE)
plot(vc)

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

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)
# Menampilkan daftar model variogram
show.vgms(par.strip.text = list(cex = 0.75))

# Menentukan model variogram awal dengan parameter psill, model, range, dan nugget
vinitial <- vgm(psill = var(log(d$vble)), model = "Sph", range = max(v$dist)/3, nugget = 0.1)

# Memvisualisasikan variogram dengan model awal (vinitial) dan batas cutoff max distance
plot(v, vinitial, cutoff = max(v$dist), cex = 1.5)

# Menyesuaikan variogram dengan model "Sph" menggunakan data v dan parameter
fv <- fit.variogram(object = v, model = vgm(psill = var(log(d$vble)), model = "Sph", range = max(v$dist)/3, nugget = 0.1))
fv
##   model      psill    range
## 1   Nug 0.05645424    0.000
## 2   Sph 0.51794475 7180.647
plot(v, fv, cex = 1.5)

# Autofit variogram untuk data log-vble
afvg <- autofitVariogram(log(vble) ~ 1, d_sp)

# Plot hasil autofit
plot(afvg)

# Membuat objek gstat dengan model variogram hasil fitting (fv) untuk data vble
k <- gstat(formula = log(vble) ~ 1, data = d_sp, model = fv)

# Melakukan prediksi variabel pada grid menggunakan model gstat yang sudah dibuat
kpred <- predict(k, coop_sp)
## [using ordinary kriging]
# Mengonversi kembali hasil prediksi ke sf untuk visualisasi ggplot
kpred_sf <- st_as_sf(kpred)
# Visualisasi prediksi log(vble) dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred_sf, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = "log(price)") + theme_bw()

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

Bandingkan RMSE

# Tabel RMSE untuk masing-masing fold (1 s.d. 5)
df_rmse_folds <- data.frame(
  Closest_Obs = rmse1,
  IDW = rmse2,
  Nearest_Neigh = rmse3,
  Ensemble = rmse4,
  Kriging = rmse5
)
print("RMSE per Fold")
## [1] "RMSE per Fold"
print(df_rmse_folds)
##   Closest_Obs      IDW Nearest_Neigh Ensemble  Kriging
## 1    49372.76 53183.18      46089.06 41646.79 38715.07
## 2    34011.64 34974.97      33786.93 28534.51 28591.93
## 3    34553.57 50871.90      41513.72 36590.89 36612.29
## 4    49888.53 60335.52      46986.69 44140.44 40441.00
## 5    31388.49 33249.47      30058.58 26717.89 29113.38
# Tabel Rata-rata RMSE
df_rmse_avg <- data.frame(
  Closest_Obs = mean(rmse1, na.rm = TRUE),
  IDW = mean(rmse2, na.rm = TRUE),
  Nearest_Neigh = mean(rmse3, na.rm = TRUE),
  Ensemble = mean(rmse4, na.rm = TRUE),
  Kriging = mean(rmse5, na.rm = TRUE)
)
print("Rata-rata RMSE (5-Fold Cross Validation)")
## [1] "Rata-rata RMSE (5-Fold Cross Validation)"
print(df_rmse_avg)
##   Closest_Obs      IDW Nearest_Neigh Ensemble  Kriging
## 1       39843 46523.01         39687 35526.11 34694.73