Tugas Praktikum 5 Mata Kuliah STA2553-Analisis Statistika Spasial

Author

Rosita Ria Rusesta/M0501251016 | Markazul Adabiyah/M0501251035 | Estyaningsi Purnama/M0501251058

Interpolasi Spasial - Fine Tuning Inverse Distance Weighting (IDW)

Load Libraries

library(sf)
library(tidyverse)
library(ggspatial)
library(gstat)

Data

Dataset yang digunakan adalah data transaksi jual beli rumah (Price Paid Data) tahun 2021 yang dirilis oleh HM Land Registry Inggris.

## 1. House price data 
hp = read_csv("pp-2021.csv", col_names = FALSE)
head(hp)
# A tibble: 6 × 16
  X1            X2 X3                  X4    X5    X6    X7    X8    X9    X10  
  <chr>      <dbl> <dttm>              <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 {E073986… 170000 2021-11-29 00:00:00 LN9 … S     N     F     6     <NA>  SAXO…
2 {E073986… 186690 2021-10-05 00:00:00 NG34… D     N     F     1     <NA>  ORCH…
3 {E073986… 150000 2021-12-03 00:00:00 LN11… S     N     F     9     <NA>  BULL…
4 {E073986… 110000 2021-12-10 00:00:00 LN12… T     N     F     134   <NA>  HIGH…
5 {E073986… 210000 2021-12-07 00:00:00 PE11… D     N     F     12    <NA>  CROW…
6 {E073986… 230000 2021-12-16 00:00:00 PE22… D     N     F     1     <NA>  PRIN…
# ℹ 6 more variables: X11 <chr>, X12 <chr>, X13 <chr>, X14 <chr>, X15 <chr>,
#   X16 <chr>
## 2. Postcode look up (November 2021) to locate house price data
pc_lut <- read_csv("NSPL_NOV_2021_UK.csv")
head(pc_lut)
# A tibble: 6 × 41
  pcd     pcd2     pcds  dointr doterm usertype oseast1m osnrth1m osgrdind oa11 
  <chr>   <chr>    <chr>  <dbl>  <dbl>    <dbl>    <dbl> <chr>       <dbl> <chr>
1 AB1 0AA AB1  0AA AB1 … 198001 199606        0   385386 0801193         1 S000…
2 AB1 0AB AB1  0AB AB1 … 198001 199606        0   385177 0801314         1 S000…
3 AB1 0AD AB1  0AD AB1 … 198001 199606        0   385053 0801092         1 S000…
4 AB1 0AE AB1  0AE AB1 … 199402 199606        0   384600 0799300         8 S000…
5 AB1 0AF AB1  0AF AB1 … 199012 199207        1   384460 0800660         8 S000…
6 AB1 0AG AB1  0AG AB1 … 199012 199207        1   383890 0800710         8 S000…
# ℹ 31 more variables: cty <chr>, ced <chr>, laua <chr>, ward <chr>,
#   hlthau <chr>, nhser <chr>, ctry <chr>, rgn <chr>, pcon <chr>, eer <chr>,
#   teclec <chr>, ttwa <chr>, pct <chr>, itl <chr>, park <chr>, lsoa11 <chr>,
#   msoa11 <chr>, wz11 <chr>, ccg <chr>, bua11 <chr>, buasd11 <chr>,
#   ru11ind <chr>, oac11 <chr>, lat <dbl>, long <dbl>, lep1 <chr>, lep2 <chr>,
#   pfa <chr>, imd <dbl>, calncv <chr>, stp <chr>
## 3. Leeds wards
wards <- st_read("Wards_December_2011_GCB_EW_2022_-8362639549357908013.gpkg")
Reading layer `Wards_December_2011_GCB_EW' from data source 
  `/Users/rositariarusesta/Downloads/Wards_December_2011_GCB_EW_2022_-8362639549357908013.gpkg' 
  using driver `GPKG'
Simple feature collection with 8588 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 82672 ymin: 5342.7 xmax: 655604.7 ymax: 657532.6
Projected CRS: OSGB36 / British National Grid

Tahap Pre-processing Data

Tahap ini bertujuan untuk mengintegrasikan data transaksi dengan informasi spasial guna membentuk dataset yang valid untuk analisis.

Integrasi dan Pembersihan Data:

  1. Data Lokalisasi: Batas wilayah (wards) Leeds disaring menggunakan kode 1411 hingga 1443 untuk membatasi ruang lingkup studi.
  2. Penyelarasan Kode Pos: Spasi pada kode pos dihapus agar proses inner join antara data harga rumah (hp) dan tabel referensi koordinat (pc_lut) berjalan akurat.
  3. Penggabungan Atribut: Atribut harga transaksi dihubungkan dengan koordinat Easting dan Northing untuk memberikan konteks lokasi pada setiap transaksi.
# filter for Leeds - the codes in the format E05000001  
wards = 
  wards |> 
  mutate(ward_number = substr(wd11cd, 4,9)) |>
  mutate(ward_number = as.double(ward_number)) |>
  filter(ward_number >= 1411 & ward_number <= 1443)
## Link Easting and Northing to house sales data postcodes
hp |> 
  # select price postcode and rename
  dplyr::select(X2, X4) %>%
  rename(price = X2, postcode = X4) |>
  # strip out white space from the postcode
  mutate(postcode = str_replace_all(postcode, " ", "")) |>
  # join to postcode lookup table
  inner_join(
    # BUT we need to process the pc_lut before the join
    # we do this inside the join
    pc_lut %>% dplyr::select(pcd, oseast1m, osnrth1m) |>
      # rename
      rename(postcode = pcd, easting = oseast1m, northing = osnrth1m) |>
      mutate(postcode = str_replace_all(postcode, " ", ""))
  ) -> hp
hp
# A tibble: 1,270,266 × 4
    price postcode easting northing
    <dbl> <chr>      <dbl> <chr>   
 1 170000 LN96PR    526756 0369333 
 2 186690 NG349JB   514955 0342721 
 3 150000 LN118UW   538771 0387266 
 4 110000 LN121DG   550314 0384956 
 5 210000 PE112HU   525964 0323853 
 6 230000 PE220JG   538812 0344780 
 7 180500 LN60LD    493729 0369988 
 8 300000 NG340QA   511674 0334261 
 9 150000 LN35SS    512141 0369482 
10 260000 PE91QG    503290 0307429 
# ℹ 1,270,256 more rows
# convert hp to sf object
hp <- st_as_sf(hp, coords = c("easting", "northing"), crs = 27700)
# re-project wards
wards <- wards %>% st_transform(27700)
# clip to wards 
hp = hp[wards, ]
ggplot(hp) + 
  annotation_map_tile(zoom = 10, progress="none") + 
  geom_sf(alpha = 0.1)

summary(hp$price)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     500   143000   205000   291355   300000 59660000

Transformasi Spasial:

  1. Proyeksi Koordinat: Dalam analisis spasial, semua layer harus memiliki CRS yang sama. Dataset dikonversi ke sistem British National Grid (EPSG: 27700) untuk menjamin akurasi spasial di wilayah Inggris.
  2. Spatial Clipping: Dilakukan pemotongan data secara otomatis sehingga hanya titik transaksi yang berada di dalam poligon batas wilayah Leeds yang dianalisis.
  3. Validasi Awal: Ringkasan data menunjukkan rentang harga ekstrem (£500 hingga £59.6 juta), yang mengindikasikan perlunya tahap fine-tuning model untuk menangani pencilan (outliers) tersebut.

Hasil akhir dari tahap pra-pemrosesan ini adalah dataset spasial terverifikasi yang disimpan dalam format .gpkg dan .RData, yang menjamin konsistensi data untuk tahap interpolasi menggunakan metode Inverse Distance Weighting (IDW).

sf::st_write(hp, "hp.gpkg", delete_dsn = TRUE)
Deleting source `hp.gpkg' using driver `GPKG'
Writing layer `hp' to data source `hp.gpkg' using driver `GPKG'
Writing 16137 features with 2 fields and geometry type Point.
st_write(wards, "wards.gpkg",delete_dsn = TRUE)
Deleting source `wards.gpkg' using driver `GPKG'
Writing layer `wards' to data source `wards.gpkg' using driver `GPKG'
Writing 33 features with 6 fields and geometry type Multi Polygon.
# OR
save(list = c("hp", "wards"), file = "point_data.RData")

Fine Tuning Inverse Distance Weighting (IDW)

Proses analisis diawali dengan pembentukan grid lokasi regular sebagai dasar estimasi menggunakan fungsi st_sample. Sebanyak 2,000 titik sampel dibuat secara sistematis di dalam batas wilayah Leeds (wards). Grid ini berfungsi sebagai lokasi target di mana nilai harga rumah akan diprediksi berdasarkan data titik transaksi aktual (hp).

# make the grid of locations
set.seed(20200912)
s.grid = st_sample(wards, 2000, type = "regular")
plot(s.grid, pch = 4)

# estimate price over those locations 
idw.est <- idw(formula = price~1, locations = hp, 
               newdata = s.grid, idp = 1.0)
[inverse distance weighted interpolation]
# examine the results 
plot(idw.est[, "var1.pred"], pch =16)

Estimasi menggunakan algoritma IDW dengan parameter \(IDP = 1.0\) menghasilkan permukaan harga yang relatif halus karena pengaruh peluruhan jarak yang bersifat linier. Secara spasial, terdapat kecenderungan bahwa area Utara Leeds memiliki estimasi harga yang lebih tinggi dibandingkan area Selatan. Namun, munculnya warna kuning yang mencolok dengan nilai estimasi di atas £450,000 mengindikasikan adanya pengaruh kuat dari titik pencilan (outliers) harga ekstrem yang mendistorsi rata-rata lokal di lokasi tersebut. Sebagian besar wilayah Leeds terdistribusi pada rentang harga yang lebih rendah, yaitu di bawah £350,000.

Sebagai langkah penyesuaian model, fungsi IDW dimodifikasi dengan menambahkan argumen nmax = 100. Secara metodologis, pengaturan ini bertujuan untuk meningkatkan lokalisasi estimasi dengan membatasi proses interpolasi hanya pada 100 tetangga terdekat di setiap lokasi grid. Pembatasan tersebut mengurangi kontribusi observasi yang lebih jauh dalam perhitungan bobot, sehingga pengaruhnya terhadap nilai prediksi menjadi semakin kecil.

Secara visual, konfigurasi ini menghasilkan permukaan interpolasi yang lebih kontras dan terlokalisasi, di mana pengaruh transaksi bernilai ekstrem cenderung terkonsentrasi di sekitar titik asalnya, tanpa menyebar secara luas ke wilayah lain.

idw.est <- idw(price~1,hp, newdata = s.grid, idp = 1.0, nmax = 100)
[inverse distance weighted interpolation]
plot(idw.est[, "var1.pred"], pch =16)

Penerapan parameter nmax = 100 pada model Inverse Distance Weighting (IDW) memberikan dampak signifikan terhadap lokalisasi estimasi spasial. Dengan membatasi pengaruh data hanya pada 100 tetangga terdekat, algoritma mampu mereduksi generalisasi berlebih yang sebelumnya muncul pada model global. Secara visual, transisi ini menghasilkan peta dengan kontras yang lebih tajam, di mana pengaruh nilai pencilan (outlier) ekstrem kini terisolasi secara presisi pada lokasi asalnya. Fenomena ini terlihat pada kemunculan hotspot berwarna kuning yang mencapai estimasi harga £5,000,000, namun tidak lagi mendistorsi nilai rata-rata harga di wilayah sekitarnya.

Sebagian besar wilayah Leeds menunjukkan estimasi harga pada kisaran ratusan ribu pound (sekitar £150,000–£350,000), yang lebih konsisten dengan distribusi median transaksi. Hasil ini menunjukkan bahwa akurasi model IDW tidak hanya dipengaruhi oleh parameter Power (idp), tetapi juga oleh jumlah tetangga (nmax) yang dilibatkan dalam estimasi.

# create a model
mod.idw <- gstat(formula = price~1, locations = hp, nmax = 100, set=list(idp = 1)) 
# apply to the data
hp$pred.idw <- predict(mod.idw, newdata=hp, debug.level=0)$var1.pred
summary(hp$pred.idw)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     500   143000   202500   293006   300000 59660000 
RMSE <- function(observed, predicted) { 
  res = sqrt(mean((predicted - observed)^2, na.rm=TRUE))
  return(res)
}
RMSE.all <- RMSE(hp$price, hp$pred.idw)
print(RMSE.all)
[1] 363414

Nilai RMSE sebesar 491,452.4 menunjukkan bahwa rata-rata kesalahan prediksi harga rumah oleh model sekitar £491,452.4. Besarnya nilai ini tidak terlepas dari distribusi harga yang sangat menceng ke kanan (right-skewed), dengan keberadaan transaksi ekstrem hingga £59.6 juta. Karena RMSE menghitung kuadrat selisih prediksi terhadap observasi, nilai ekstrem memiliki pengaruh yang sangat besar terhadap besarnya error agregat. Namun demikian, perlu dicatat bahwa nilai RMSE ini dihitung menggunakan data yang sama dengan data pelatihan (in-sample evaluation), sehingga cenderung menghasilkan estimasi kesalahan yang bersifat optimistik.

Oleh karena itu, tahapan analisis selanjutnya difokuskan pada upaya penemuan kombinasi parameter kontrol IDW yang paling optimal. Tujuan dari proses ini adalah untuk menghasilkan model prediksi harga rumah dengan tingkat kesalahan (error) seminimal mungkin, guna menjamin ketepatan estimasi di seluruh wilayah studi.

Pembentukan Subset Training dan Validation

set.seed(20200912)  # for reproducibility
index <- sample(nrow(hp), 0.2 * nrow(hp))
# create subsets
trn <- hp[-index,]   # training (80%)
tst <- hp[index,]    # validation (20%)
# create a model
mod.idw.base <- gstat(formula = price~1, locations = trn, nmax = 100, set=list(idp = 1)) 
# apply to the data
tst$pred.idw <- predict(mod.idw.base, newdata=tst, debug.level=0)$var1.pred
RMSE.tst <- RMSE(tst$price, tst$pred.idw)
print(RMSE.tst)
[1] 878435.6

Setelah melakukan pembagian data 80:20, kami menghitung RMSE Baseline pada data testing menggunakan parameter standar. Hasilnya menunjukkan tingkat kesalahan sebesar 950,141.1. Angka inilah yang menjadi titik acuan kami untuk melakukan proses fine-tuning, guna mencari kombinasi parameter yang lebih presisi bagi karakteristik spasial harga properti di Leeds.

Pembangunan Tuning Grid Parameter

Kualitas interpolasi IDW sangat bergantung pada parameter Power ($p$ atau $idp$) dan jumlah tetangga (nmax). Oleh karena itu, dibangun sebuah tuning grid yang mencakup 270 kombinasi parameter berbeda, dengan rentang nmax dari 50 hingga 1500 dan nilai Power ($idp$) dari 1 hingga 5. Metrik evaluasi yang digunakan adalah Root Mean Square Error (RMSE), di mana nilai RMSE yang lebih kecil menunjukkan model yang lebih akurat dalam memprediksi data uji (validation set).

params <- expand.grid(nmax_vals = seq(50, 1500, by = 50),
                      idp_vals  = seq(1, 5, by = 0.5))
# define a result attribute in params
params$RMSE <- 0
dim(params)
[1] 270   3

Mencari Parameter Terbaik

params[which.min(params$RMSE),]
    nmax_vals idp_vals     RMSE
181        50        4 840788.7

Berdasarkan grid search awal, kombinasi terbaik sementara diperoleh pada nmax = 50 dan idp = 4 dengan RMSE 915,455.6. Nilai idp yang relatif tinggi menunjukkan bahwa peluruhan bobot jarak yang lebih tajam memberikan performa prediksi yang lebih baik pada data validasi berdasarkan RMSE, sementara nmax yang rendah (50) membantu mengisolasi pengaruh pencilan (outliers) agar tidak mendistorsi prediksi secara luas.

Optimasi Lanjutan

Tahap ini dilakukan untuk mengevaluasi stabilitas pembagian data serta memahami perilaku metrik kesalahan (RMSE) terhadap perubahan parameter kontrol IDW.

summary(trn$price)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     500   142500   205000   292497   300000 59660000 
summary(tst$price)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     500   145000   207000   286786   302250 41000000

Langkah pertama dalam validasi model adalah memastikan bahwa subset data yang digunakan untuk melatih model (trn) memiliki karakteristik distribusi yang serupa dengan subset data untuk menguji model (tst).

Berdasarkan ringkasan statistik, kedua subset menunjukkan nilai Minimum (£500) dan Median (sekitar £205,000-£207,000) yang sangat konsisten. Terdapat perbedaan pada nilai Maksimum, di mana data training memiliki pencilan sebesar £59.6 juta sedangkan data validation sebesar £41 juta. Perbedaan ini menyebabkan nilai Mean pada data training (£292,497) sedikit lebih tinggi dibandingkan data validation (£286,786). Meskipun terdapat variasi pada nilai ekstrem, kemiripan nilai kuartil dan median menunjukkan bahwa pembagian data 80:20 telah representatif untuk digunakan dalam proses evaluasi.

ggplot(params, aes(x = idp_vals, y = RMSE)) + geom_point()

Grafik hubungan antara nilai Power (idp_vals) dan RMSE memberikan gambaran mengenai tingkat kehalusan permukaan interpolasi yang optimal.

  • Tren Penurunan Error: Grafik menunjukkan tren penurunan nilai RMSE yang signifikan seiring bertambahnya nilai idp dari 0.5 hingga 4.5.

  • Titik Optimal: Nilai RMSE terendah mulai tercapai pada kisaran idp 4.0 hingga 4.5. Hasil tuning menunjukkan bahwa peningkatan nilai idp hingga kisaran tertentu berkorelasi dengan penurunan RMSE pada data validasi, yang mengindikasikan bahwa konfigurasi dengan peluruhan bobot jarak yang lebih tajam lebih sesuai untuk struktur data ini berdasarkan metrik RMSE yang digunakan.

ggplot(params, aes(x = nmax_vals, y = RMSE, colour = factor(idp_vals))) + geom_point()

Grafik multifaktor menunjukkan bagaimana jumlah tetangga yang dilibatkan dalam estimasi memengaruhi akurasi model pada berbagai tingkat Power.

  • Efek Lokalisasi: Pada semua nilai idp, penggunaan nmax yang kecil secara konsisten menghasilkan RMSE yang lebih rendah dibandingkan penggunaan jumlah tetangga yang lebih besar.

  • Saturasi Error: Seiring bertambahnya nmax, nilai RMSE cenderung naik dan mencapai saturasi (datar) setelah melewati angka 500 tetangga. Hasil ini menunjukkan bahwa pada dataset ini, penggunaan jumlah tetangga yang terlalu besar berkorelasi dengan peningkatan nilai RMSE pada data validasi.

  • Interaksi Parameter: Kombinasi terbaik (ditunjukkan oleh warna dengan nilai RMSE terendah pada grafik) dicapai ketika model menggunakan idp tinggi (3.5–5) dan nmax rendah (5–100).

params_fine <- expand.grid(nmax_vals = seq(5, 100, by = 5),
                           idp_vals  = seq(3.5, 5, by = 0.1))
# define a result attribute in params
params_fine$RMSE <- 0

Evaluasi Kombinasi Parameter (Grid Search)

for(i in 1:nrow(params_fine)) {
  # make a model
  model.i = gstat(formula = price~1, 
                   locations = trn, 
                   nmax = params_fine$nmax_vals[i], 
                   set=list(idp = params_fine$idp_vals[i])) 
  # apply to the data
  predicted.i = predict(model.i, newdata = tst, debug.level=0)$var1.pred
  # calculate RMSE and write to the tune grid
  params_fine$RMSE[i] = RMSE(tst$price, predicted.i)
  # counter 
  if(i%%10 == 0) cat(i, "\t")
}
10  20  30  40  50  60  70  80  90  100     110     120     130     140     150     160     170     180     190     200     210     220     230     240     250     260     270     280     290     300     310     320

Mencari Parameter Terbaik

params_fine[which.min(params_fine$RMSE),]
  nmax_vals idp_vals     RMSE
2        10      3.5 831784.5
# assign the values 
nmax_final = params_fine[which.min(params_fine$RMSE),"nmax_vals"]
idp_final = params_fine[which.min(params_fine$RMSE),"idp_vals"]

Implementasi Model Final dan Visualisasi

idw.est <- idw(formula = price~1,locations = hp, 
               newdata = s.grid, idp = idp_final, nmax = nmax_final)
[inverse distance weighted interpolation]
ggplot(data = idw.est) + 
  geom_sf(aes(col = var1.pred), size = 2) +
  scale_colour_viridis_c(
    name = "Estimated Price (£)", 
    option = "plasma", 
    labels = scales::comma
  ) +
  theme_minimal()

Pemetaan permukaan harga kontinu (continuous price surface) di wilayah heterogen seperti Leeds memerlukan pengaturan parameter interpolasi yang mampu menyeimbangkan sensitivitas lokal dan stabilitas global model. Berdasarkan seluruh rangkaian analisis dan proses optimasi parameter, beberapa temuan utama dapat dirangkum sebagai berikut.

  1. Evaluasi In-Sample vs. Out-of-Sample: Nilai RMSE awal yang lebih kecil diperoleh karena evaluasi dilakukan pada data yang sama dengan data pelatihan (in-sample evaluation), sehingga menghasilkan estimasi kesalahan yang bersifat optimistik. Dalam evaluasi in-sample, model menggunakan seluruh titik observasi sebagai referensi dalam proses estimasi. Pada metode IDW, titik dengan jarak sangat kecil terhadap lokasi prediksi akan memperoleh bobot dominan, sehingga nilai prediksi cenderung sangat mendekati nilai observasi. Oleh karena itu, peningkatan nilai RMSE pada tahap validasi bukan menunjukkan penurunan kualitas model, melainkan mencerminkan evaluasi yang lebih representatif terhadap kemampuan generalisasi model.
  2. Optimasi Parameter: Melalui proses pencarian parameter yang lebih rinci (fine-grained search), diperoleh kombinasi \(nmax = 10\) dan \(idp = 3.5\) dengan nilai RMSE sebesar 907,171.5 pada data validasi. Berdasarkan nilai RMSE pada data validasi, kombinasi ini memberikan performa prediksi terbaik dibandingkan kombinasi parameter lain yang diuji.
  3. Pembatasan Propagasi Nilai Ekstrem: Penggunaan jumlah tetangga yang terbatas mengurangi kontribusi observasi yang lebih jauh dalam proses estimasi, sehingga pengaruh nilai transaksi ekstrem menjadi lebih terbatas secara spasial dibandingkan konfigurasi dengan nmax besar. Dengan dikombinasikan bersama nilai idp yang relatif tinggi, bobot pengaruh menurun dengan cepat terhadap jarak, sehingga dampak observasi bernilai sangat besar tidak menyebar luas ke wilayah lain. Pendekatan ini menghasilkan permukaan estimasi yang lebih terlokalisasi, meskipun dengan konsekuensi meningkatnya variasi lokal pada beberapa area.
  4. Interpretasi Visualisasi: Implementasi parameter optimal pada peta akhir dengan palet warna plasma menghasilkan visualisasi yang menampilkan variasi spasial harga secara lebih terfokus di sekitar lokasi transaksi bernilai tinggi. Pola yang dihasilkan merepresentasikan distribusi harga sebagaimana tercatat dalam dataset tahun 2021, serta menunjukkan bagaimana pengaturan parameter IDW memengaruhi tingkat kehalusan (smoothing) permukaan interpolasi dan sensitivitas model terhadap variasi lokal. Secara keseluruhan, hasil ini menegaskan pentingnya proses fine-tuning dalam menghasilkan estimasi spasial yang lebih stabil, terlokalisasi, dan representatif terhadap karakteristik data harga rumah di wilayah studi.