knitr::opts_chunk$set(warning = FALSE, message = FALSE)Praktikum 4 : Regresi Spasial
Analisis Eksplorasi Data Geostatistik
Data Geostatistik
Data geostatistik adalah data pengukuran suatu fenomena yang bersifat kontinu secara spasial, tetapi diamati hanya pada sejumlah lokasi tertentu. Informasi dari lokasi pengamatan digunakan untuk mempelajari struktur spasial dan memprediksi nilai pada lokasi yang tidak diamati.
(Moraga, 2023)
Contohnya antara lain:
polusi udara yang diukur pada sejumlah stasiun pemantauan,
suhu yang diukur pada sejumlah lokasi,
konsentrasi logam dalam tanah,
dan fenomena lingkungan lain yang secara konseptual dapat berubah secara kontinu dalam ruang.
Contoh data geostatistik : meuse
library(sp)
library(spdep)data(meuse)
coordinates(meuse) <- ~x+y
plot(meuse, col = heat.colors(20)[cut(meuse$zinc,20)],
main="Konsentrasi Timbal (mg/kg tanah)")
Prevalensi Malaria
library(malariaAtlas)
library(ggplot2)
# Download data malaria (Zimbabwe)
zimbabwe <- getPR(country = "Zimbabwe", species = "pf")ggplot(zimbabwe, aes(x = longitude, y = latitude, color = pr)) +
borders("world", region = "Zimbabwe", fill="gray95") +
geom_point(size = 2) +
scale_color_viridis_c() +
theme_minimal() +
ggtitle("Prevalensi Malaria di Zimbabwe")
Interpolasi Spasial
Interpolasi spasial merupakan metode yang digunakan untuk memperkirakan nilai peubah spasial pada lokasi yang tidak diketahui berdasarkan data yang tersedia pada lokasi tetangga yang diketahui (Cressie, 1993). Interpolasi ini digunakan dalam berbagai bidang seperti geografi, meteorologi, pertanian, dan lingkungan untuk menghasilkan peta atau permukaan kontinu dari data yang awalnya hanya tersedia pada titik-titik diskrit.
Interpolasi spasial mengasumsikan bahwa nilai variabel yang diukur (misalnya, suhu, ketinggian, kelembaban tanah, dll.) di suatu lokasi akan serupa dengan nilai di lokasi-lokasi tetangga. Dengan kata lain, ada korelasi spasial, di mana lokasi yang lebih dekat cenderung memiliki nilai yang lebih mirip dibandingkan lokasi yang lebih jauh.
data(meuse)
class(meuse)[1] "data.frame"head(meuse) x y cadmium copper lead zinc elev dist om ffreq soil lime
1 181072 333611 11.7 85 299 1022 7.909 0.00135803 13.6 1 1 1
2 181025 333558 8.6 81 277 1141 6.983 0.01222430 14.0 1 1 1
3 181165 333537 6.5 68 199 640 7.800 0.10302900 13.0 1 1 1
4 181298 333484 2.6 81 116 257 7.655 0.19009400 8.0 1 2 0
5 181307 333330 2.8 48 117 269 7.480 0.27709000 8.7 1 2 0
6 181390 333260 3.0 61 137 281 7.791 0.36406700 7.8 1 2 0
landuse dist.m
1 Ah 50
2 Ah 30
3 Ah 150
4 Ga 270
5 Ah 380
6 Ga 470Perhatikan bahwa data meuse masih berupa data.frame. Oleh karenanya, kita akan mendefinisikan koordinat lokasi pada data meuse, agar data tersebut menjadi class data spatial. Terlihat pula pada gambar bahwa kandungan zinc cenderung lebih tinggi pada daerah di pinggir sungai (daerah paling kiri) yang ditunjukkan dengan lingkaran yang lebih besar.
coordinates(meuse) <- ~x+y # Menentukan kolom x dan y sebagai koordinat spasial
class(meuse)[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"plot(meuse,pch=1, cex=meuse$zinc/1000)
Selanjutnya, kita perlu mempersiapkan grid untuk interpolasi. Grid membantu dalam interpolasi karena kita dapat dengan mudah melakukan prediksi di setiap sel, memastikan cakupan yang merata di seluruh wilayah studi.
data(meuse.grid)
str(meuse.grid)'data.frame': 3103 obs. of 7 variables:
$ x : num 181180 181140 181180 181220 181100 ...
$ y : num 333740 333700 333700 333700 333660 ...
$ part.a: num 1 1 1 1 1 1 1 1 1 1 ...
$ part.b: num 0 0 0 0 0 0 0 0 0 0 ...
$ dist : num 0 0 0.0122 0.0435 0 ...
$ soil : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...Sama seperti sebelumnya, kita perlu mendefinisikan koordinat lokasi untuk data meuse.grid.
gridded(meuse.grid) <- ~x+y # Menjadikan meuse.grid sebagai grid spasial
plot(meuse.grid)
Kemudian setelah tersedia data titik dan data grid, kita dapat melakukan interpolasi pada data meuse.
Metode Interpolasi Spasial
Trend Surface Analysis (Analisis Permukaan Tren)
Trend Surface menggunakan fungsi polinomial untuk memodelkan pola spasial secara umum. Fungsi ini digunakan untuk memperkirakan nilai pada lokasi yang belum diketahui.
Cocok untuk menangkap pola umum atau tren yang mendasari dalam data.
library(gstat)
# Membuat model Trend Surface
trend_model <- gstat(formula = log(zinc) ~ poly(x, y, degree = 2), data = meuse)
# Prediksi dengan model Trend Surface
trend_pred <- predict(trend_model, meuse.grid)[ordinary or weighted least squares prediction]library(sp)
library(latticeExtra)
spplot(trend_pred, "var1.pred", main = "Trend Surface Interpolation of Zinc")+layer(sp.points(meuse, pch = 19, col = "black", cex = 0.5))
gstat(formula = zinc ~ poly(x, y, degree = 2), data = meuse): Membuat model Trend Surface menggunakan polinomial derajat 2 untuk variabelzinc.predict(trend_model, meuse.grid): Menggunakan model untuk memprediksi nilaizincpada grid spasial yang ada.spplot(): Fungsi untuk memplot hasil interpolasi.
# 1. Load library yang dibutuhkan
library(sp)
library(gstat)
library(latticeExtra)
# 2. Siapkan data meuse
data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x+y
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
# 3. Buat model TSA untuk Derajat 1, 2, dan 3
m1 <- lm(log(zinc) ~ x + y, data = meuse) # Derajat 1 (Linear)
m2 <- lm(log(zinc) ~ x + y + I(x^2) + I(y^2) + I(x*y), data = meuse) # Derajat 2 (Quadratic)
m3 <- lm(log(zinc) ~ poly(x, y, degree = 3, raw = TRUE), data = meuse) # Derajat 3 (Cubic)
# -------------------------------------------------------------------
# METODE 1: Uji ANOVA (Sequential F-Test)
# -------------------------------------------------------------------
# Membandingkan apakah peningkatan derajat memberikan signifikansi statistik
anova_result <- anova(m1, m2, m3)
print(anova_result)Analysis of Variance Table
Model 1: log(zinc) ~ x + y
Model 2: log(zinc) ~ x + y + I(x^2) + I(y^2) + I(x * y)
Model 3: log(zinc) ~ poly(x, y, degree = 3, raw = TRUE)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 152 58.802
2 149 39.464 3 19.338 24.337 7.118e-13 ***
3 149 39.464 0 0.000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# -------------------------------------------------------------------
# METODE 2: Akaike Information Criterion (AIC)
# -------------------------------------------------------------------
# Nilai AIC terendah menunjukkan model paling optimal dengan memperhatikan kompleksitas
aic_values <- AIC(m1, m2, m3)
print(aic_values) df AIC
m1 4 297.6365
m2 7 241.8250
m3 7 241.8250# -------------------------------------------------------------------
# METODE 3: Cross-Validation (LOOCV)
# -------------------------------------------------------------------
# Mengukur error prediksi (RMSE) pada data yang tidak dipakai dalam fitting
cv_m1 <- krige.cv(log(zinc) ~ 1, meuse, model = NULL, degree = 1)
cv_m2 <- krige.cv(log(zinc) ~ 1, meuse, model = NULL, degree = 2)
cv_m3 <- krige.cv(log(zinc) ~ 1, meuse, model = NULL, degree = 3)
rmse_cv <- data.frame(
Derajat = c(1, 2, 3),
RMSE = c(
sqrt(mean(cv_m1$residual^2)),
sqrt(mean(cv_m2$residual^2)),
sqrt(mean(cv_m3$residual^2))
)
)
print(rmse_cv) Derajat RMSE
1 1 0.6343133
2 2 0.5376588
3 3 0.6486922# -------------------------------------------------------------------
# PREDIKSI DAN VISUALISASI (Derajat Optimal: Derajat 2)
# -------------------------------------------------------------------
# interpollasi tren ke grid menggunakan krige() dengan derajat pilihan
trend_pred <- krige(log(zinc) ~ 1, meuse, meuse.grid, degree = 2)[ordinary or weighted least squares prediction]# Visualisasi spplot sesuai sintaks Anda
spplot(trend_pred, "var1.pred", main = "Trend Surface Interpolation of Zinc (Degree 2)") +
layer(sp.points(meuse, pch = 19, col = "black", cex = 0.5))
Thiessen Polygons
Metode ini membagi wilayah menjadi poligon di mana setiap poligon berisi satu titik pengamatan, dan setiap lokasi dalam poligon lebih dekat ke titik itu daripada ke titik lainnya.
library(gstat)
zn.tp <- krige(log(zinc) ~ 1, meuse, meuse.grid, nmax = 1)[inverse distance weighted interpolation]spplot(zn.tp, "var1.pred", main = "Thiessen Polygons Interpolation of Zinc") +
layer(sp.points(meuse, pch = 19, col = "black", cex = 0.5))
Inverse Distance Weight
IDW mengasumsikan bahwa nilai di suatu titik tak diketahui adalah rata-rata tertimbang dari nilai-nilai di titik-titik terdekat, dengan bobot yang berkurang seiring dengan meningkatnya jarak dari titik tersebut.
Semakin dekat titik pengamatan, semakin besar pengaruhnya terhadap nilai yang diestimasi.
# Membuat model IDW
idw_model <- gstat(formula = log(zinc) ~ 1, data = meuse, nmax = 7, set = list(idp = 2))
# Prediksi dengan model IDW
idw_pred <- predict(idw_model, meuse.grid)[inverse distance weighted interpolation]gstat(formula = zinc ~ 1, data = meuse, nmax = 7, set = list(idp = 2)): Membuat model IDW untuk variabelzinc, menggunakan parameteridp = 2untuk bobot jarak (semakin dekat titik, semakin besar bobotnya) dannmax = 7untuk jumlah titik yang digunakan dalam prediksi.Ketika
idplebih besar, pengaruh titik-titik yang lebih dekat akan lebih kuat, sementara titik yang lebih jauh akan memiliki pengaruh yang lebih kecil.Jika
nmaxdiatur ke angka tertentu (misalnya, 7), hanya 7 titik pengamatan terdekat yang akan digunakan untuk menghitung nilai yang diinterpolasi di setiap lokasi.
predict(idw_model, meuse.grid): Menggunakan model IDW untuk memprediksi nilaizincpada grid.
spplot(idw_pred, "var1.pred", main = "IDW Interpolation of Zinc") +
layer(sp.points(meuse, pch = 19, col = "black", cex = 0.5))
K-Point Means Interpolation
# K-Point Means Interpolation
kpoint_model <- krige(log(zinc) ~ 1, meuse, meuse.grid,
nmax = 5, set = list(idp = 0)) # rata-rata dari 5 tetangga terdekat[inverse distance weighted interpolation]spplot(kpoint_model, "var1.pred", main = "K-Point Means Interpolation") +
layer(sp.points(meuse, pch = 19, col = "black", cex = 0.5))
plots <- list(
"Trend Surface" = spplot(trend_pred, "var1.pred", main = "Trend Surface"),
"Thiessen" = spplot(zn.tp, "var1.pred", main = "Thiessen Polygons"),
"IDW" = spplot(idw_pred, "var1.pred", main = "IDW"),
"K-Point Means" = spplot(kpoint_model, "var1.pred", main = "K-Point Means")
)
gridExtra::grid.arrange(grobs = plots, ncol = 2,
top = "Perbandingan Metode Interpolasi Zinc (Meuse)")
| Metode | Karakteristik Hasil | Kelebihan | Kelemahan | ||||
| Trend Surface | Permukaan sangat halus, berbentuk gradien global. | Menangkap pola umum/tren global; mudah dipahami. | Mengabaikan detail lokal; bisa terlalu menyederhanakan data. | ||||
| Thiessen | Permukaan berupa poligon diskrit, tiap wilayah mewarisi nilai 1 titik terdekat. | Mudah dihitung; jelas batas wilayah pengaruh tiap titik. | Kasar; tidak cocok untuk variabel kontinu (hasil terlihat “blocky”). | ||||
| IDW | Permukaan halusdengan variasi lokal jelas; titik terdekat punya bobot besar. | Dapat menggambarkan variasi lokal; intuitif karena berbasis jarak. | Bisa overfitting; sensitif terhadap parameter idpdan distribusi titik. |
||||
| K-Point Means | Permukaan halus, lebih merata dibanding IDW; tiap prediksi rata-rata dari k titik. | Lebih stabil dari IDW; tidak terlalu dipengaruhi 1 titik ekstrim. | Variasi kecil hilang; hasil bisa terlalu “datar” jika k terlalu besar. |
Trend Surface → cocok kalau ingin tahu pola global
Thiessen → cocok kalau mau tahu “wilayah pengaruh” masing-masing titik sampel.
IDW → cocok untuk data kontinu dengan asumsi jarak berperan besar.
K-Point Means → cocok kalau ingin hasil halus tapi tidak terlalu dipengaruhi titik tunggal.
Referensi
Anisa, R. (2020). Interpolasi Spasial. https://rpubs.com/r_anisa/spatial-interpolation (accessed at 03 September 2024)
Moraga, P. (2024). Spatial Statistics for Data Science Theory and Practice with R. Boca Raton: CRC Press.
Pebesma, E & Bivand, R. (2023). Spatial Data Science with Applications in R. New York: CRC Press.