Tugas Praktikum Regspas 10

Library

library(sp)
## Warning: package 'sp' was built under R version 4.5.3
library(sf)
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(terra)
## Warning: package 'terra' was built under R version 4.5.3
## terra 1.9.46
library(tmap)
## Warning: package 'tmap' was built under R version 4.5.3
library(viridis)
## Warning: package 'viridis' was built under R version 4.5.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.5.3
library(dismo)
## Warning: package 'dismo' was built under R version 4.5.3
## Loading required package: raster
## Warning: package 'raster' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## 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)
## Warning: package 'gstat' was built under R version 4.5.3
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(geodata)
## Warning: package 'geodata' was built under R version 4.5.3

Data

x <- read_excel("C:/Users/Abit/Downloads/[Metode Baru] Rata-rata Lama Sekolah, 2026.xlsx",
                skip = 3, col_names = c("wilayah", "vble"))
x$vble <- suppressWarnings(as.numeric(x$vble))
 
# Ambil baris provinsi (huruf kapital); buang "INDONESIA" dan provinsi Papua
# (batas GADM masih sebelum pemekaran Papua)
x <- x %>% filter(wilayah == toupper(wilayah), wilayah != "INDONESIA", !grepl("PAPUA", wilayah))
 
# Kunci penggabungan nama
key <- function(n) tolower(gsub("[^A-Za-z]", "", n))
x$key <- recode(key(x$wilayah), kepbangkabelitung = "bangkabelitung",
                dkijakarta = "jakartaraya", diyogyakarta = "yogyakarta")
crs <- "+proj=laea +lat_0=-2 +lon_0=118 +datum=WGS84 +units=m"
g <- gadm("IDN", level = 1, path = tempdir()) %>% st_as_sf() %>% st_transform(crs)
## Cached as: C:\Users\Abit\AppData\Local\Temp\RtmpYfDLnb/gadm/gadm41_IDN_1_pk.rds
g <- st_simplify(st_make_valid(g), dTolerance = 5000)
g$key <- key(g$NAME_1)
 
# Peta batas Indonesia dan nama yang belum cocok
map <- st_union(g) %>% st_sf()
setdiff(x$key, g$key)
## character(0)
# Titik pengamatan = centroid provinsi
d <- inner_join(g, x, by = "key") %>% st_centroid() %>% select(vble)
## Warning: st_centroid assumes attributes are constant over geometries
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.
tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
  tm_dots("vble", palette = "viridis")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

# Grid raster 20 km dan lokasi prediksi
grid <- rast(map, resolution = 20000)
xy <- xyFromCell(grid, 1:ncell(grid))
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))
coop <- st_filter(coop, map)

Closest Observation (Voronoi)

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

v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)
## 
## ── 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' 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`.

resp <- coop
resp$pred <- d$vble[st_nearest_feature(coop, d)]
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c1 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
## 
## ── 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`.
c1

# Inverse Distance Weighting (IDW)

res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$pred <- resp$var1.pred
 
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
## 
## ── 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`.
c2

# Nearest Neighbors (k = 5)

res <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$pred <- resp$var1.pred
 
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c3 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
## 
## ── 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`.
c3

# Ordinary Kriging

vg <- variogram(vble ~ 1, data = d)
fv <- fit.variogram(object = vg, model = vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1))
## Warning in fit.variogram(object = vg, model = vgm(psill = 0.5, model = "Sph", :
## singular model in variogram fit
fv
##   model psill range
## 1   Nug   0.1     0
## 2   Sph   0.5   900
plot(vg, fv, cex = 1.5)

k <- gstat(formula = vble ~ 1, data = d, model = fv)
resp <- predict(k, coop)
## [using ordinary kriging]
resp$pred <- resp$var1.pred
 
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
## 
## ── 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`.
c4
## The map variable `col` of the layer "raster" contains a unique value. Therefore
## a discrete scale is applied (tm_scale_discrete).

vg <- variogram(vble ~ 1, data = d)
plot(vg)

vinitial <- vgm(psill = 0.5, model = "Sph", range = 1000000, nugget = 0.3)
fv <- fit.variogram(object = vg, model = vinitial)
## Warning in fit.variogram(object = vg, model = vinitial): No convergence after
## 200 iterations: try different initial values?
fv
##   model     psill   range
## 1   Nug 0.4754689       0
## 2   Sph 0.7217403 3373121

Ensembel (Closest Obs + IDW + NN + Kriging, bobot sama)

p1 <- d$vble[st_nearest_feature(coop, d)]
p2 <- predict(gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1)), coop)$var1.pred
## [inverse distance weighted interpolation]
p3 <- predict(gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0)), coop)$var1.pred
## [inverse distance weighted interpolation]
p4 <- predict(gstat(formula = vble ~ 1, data = d, model = fv), coop)$var1.pred
## [using ordinary kriging]
weights <- c(1/4, 1/4, 1/4, 1/4)
p5 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3] + p4 * weights[4]
 
resp <- coop
resp$pred <- p5
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c5 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
## 
## ── 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`.
c5

# Semua Peta

tmap_arrange(c1, c2, c3, c4, c5, 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.
## The map variable `col` of the layer "raster" contains a unique value. Therefore
## a discrete scale is applied (tm_scale_discrete).

Evaluasi Model: 5-fold Cross Validation

set.seed(123)
 
RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}
 
kf <- kfold(nrow(d), k = 5)
 
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Kriging
rmse5 <- rep(NA, 5) # Ensemble
 
for(k in 1:5){
  test <- d[kf == k, ]
  train <- d[kf != k, ]
 
  # Closest observation
  p1 <- train$vble[st_nearest_feature(test, train)]
  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
  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)
 
  # Kriging (variogram dihitung ulang dari data train)
  vtr <- variogram(vble ~ 1, data = train)
  fvtr <- fit.variogram(vtr, model = vinitial)
  ok <- gstat(formula = vble ~ 1, data = train, model = fvtr)
  p4 <- predict(ok, test)$var1.pred
  rmse4[k] <- RMSE(test$vble, p4)
 
  # Ensemble (bobot = invers RMSE)
  w <- 1 / c(rmse1[k], rmse2[k], rmse3[k], rmse4[k])
  weights <- w / sum(w)
  p5 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3] + p4 * weights[4]
  rmse5[k] <- RMSE(test$vble, p5)
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## Warning in fit.variogram(vtr, model = vinitial): No convergence after 200
## iterations: try different initial values?
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## Warning in fit.variogram(vtr, model = vinitial): No convergence after 200
## iterations: try different initial values?
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## Warning in fit.variogram(vtr, model = vinitial): No convergence after 200
## iterations: try different initial values?
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## Warning in fit.variogram(vtr, model = vinitial): No convergence after 200
## iterations: try different initial values?
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## Warning in fit.variogram(vtr, model = vinitial): No convergence after 200
## iterations: try different initial values?
## [using ordinary kriging]
# RMSE tiap fold
data.frame(closest.obs = rmse1, IDW = rmse2, nearest.neigh = rmse3,
           kriging = rmse4, ensemble = rmse5)
##   closest.obs       IDW nearest.neigh   kriging  ensemble
## 1   1.2581070 0.7020675     0.7322213 0.7531451 0.7564833
## 2   1.2704667 1.0963358     1.1648237 1.1294399 1.1304561
## 3   0.5320088 0.5846024     0.7915424 0.6235385 0.5538616
## 4   1.2307141 0.9614094     1.0086288 0.9270655 0.9953090
## 5   1.1305308 1.1418525     1.2539577 1.1674295 1.1475052
# Rata-rata RMSE dari 5 fold
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2), nearest.neigh = mean(rmse3),
           kriging = mean(rmse4), ensemble = mean(rmse5))
##   closest.obs       IDW nearest.neigh   kriging  ensemble
## 1    1.084365 0.8972535     0.9902348 0.9201237 0.9167231
  1. IDW (0.897): paling kecil. IDW merata-ratakan semua provinsi dengan bobot jarak, jadi prediksinya “halus” mendekati rata-rata nasional. Lama sekolah antarprovinsi tidak terlalu jauh bedanya (sekitar 8 sampai 12 tahun), jadi prediksi yang halus ini justru aman.
  2. Ensembel (0.917): hampir sama dengan IDW. Menggabungkan beberapa metode membuat galatnya lebih stabil daripada metode tunggal yang lemah, tapi belum mengalahkan IDW.
  3. Kriging (0.920): selisihnya dengan IDW dan Ensembel sangat tipis. Kriging membutuhkan variogram yang bagus, sedangkan di sini hanya ada 32 titik provinsi, sehingga variogramnya kurang stabil dan keunggulan Kriging tidak muncul.
  4. Nearest Neighbors (1.006): memakai 5 tetangga terdekat. Provinsi yang bertetangga ternyata tidak selalu mirip nilainya (misalnya Kepulauan Riau 10.84 dan Riau 9.56, DIY 10.46 dan Jawa Tengah 8.16), jadi tetangga dekat tidak selalu jadi penduga yang baik.
  5. Closest Observation (1.084): paling kurang akurat, karena hanya meniru satu provinsi terdekat. Satu provinsi yang nilainya berbeda langsung membuat prediksi meleset, terutama di fold 1, 2, dan 4 (RMSE di atas 1.2).