Spatial Statistics
From Frequentis to Bayesian
I Gede Nyoman Mindra Jaya, Ph.D
04 December 2025
Kata Pengantar
1 Konsep Dasar Data Spasial
Mata kuliah Spatial Statistics dimulai dengan memahami konsep paling dasar yaitu lokasi, jarak, dan matriks bobot spasial. Tiga konsep ini merupakan fondasi penting sebelum masuk ke analisis spasial yang lebih lanjut seperti autokorelasi, regresi spasial, maupun geostatistik.
1.1 Lokasi dalam Data Spasial
1.1.1 Definisi
Lokasi adalah representasi posisi objek pada ruang geografis. Dalam data spasial, lokasi dapat direpresentasikan dalam berbagai bentuk: - Koordinat titik (point): misalnya lokasi rumah sakit, sekolah, kasus penyakit. - Garis (line): misalnya jalan, sungai, jaringan listrik. - Poligon (area): misalnya batas kecamatan, provinsi, atau wilayah administrasi.
Lokasi dapat ditentukan melalui sistem koordinat geografis (latitude–longitude) atau sistem koordinat proyeksi (misalnya UTM). Pemilihan sistem koordinat memengaruhi interpretasi jarak dan arah.
1.1.2 Representasi Data Lokasi
- Data vektor: menyimpan lokasi sebagai titik, garis, atau poligon.
- Data raster: membagi ruang menjadi grid/pixel; tiap sel mewakili area kecil dengan nilai tertentu.
1.1.3 Pentingnya Lokasi
Lokasi memungkinkan kita menghubungkan data dengan ruang. Misalnya: - Kasus demam berdarah di setiap kelurahan Bandung. - Pusat kesehatan terdekat bagi penduduk. - Perubahan tutupan lahan pada suatu wilayah.
1.2 Jarak dalam Analisis Spasial
1.2.1 Definisi Jarak
Jarak adalah ukuran kedekatan atau keterpisahan antara dua lokasi. Dalam data spasial, jarak sangat penting karena sering kali menentukan hubungan antar-unit.
1.2.2 Jenis Ukuran Jarak
Euclidean distance
Jarak garis lurus antara dua titik \((x_i, y_i)\) dan \((x_j, y_j)\):\[ d_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} \]
Manhattan distance
Jarak berbasis grid (seperti jalan kota):\[ d_{ij} = |x_i - x_j| + |y_i - y_j| \]
Great-circle distance (Haversine)
Digunakan pada permukaan bumi (koordinat lintang–bujur).
1.2.3 Pentingnya Jarak
- Menentukan kedekatan spasial: kasus penyakit cenderung berkelompok pada jarak tertentu.
- Digunakan dalam membangun matriks bobot spasial.
- Menjadi dasar dalam metode interpolasi dan klastering spasial.
1.3 Matriks Bobot Spasial
1.3.1 Definisi
Matriks bobot spasial (\(W\)) adalah representasi formal dari keterkaitan antar-unit spasial. Elemen \(w_{ij}\) menyatakan tingkat kedekatan atau hubungan antara unit \(i\) dan \(j\).
1.3.2 Bentuk Matriks Bobot
- Contiguity-based (berbatasan)
- Rook contiguity: dua wilayah bertetangga jika berbagi sisi.
- Queen contiguity: dua wilayah bertetangga jika berbagi sisi atau sudut.
- Distance-based (berbasis jarak)
Unit dianggap bertetangga jika jarak \(d_{ij} \leq d_0\) (threshold).
Bobot bisa diberikan sebagai fungsi jarak, misalnya:
\[ w_{ij} = \frac{1}{d_{ij}^2} \]
1.3.3 Normalisasi
Matriks \(W\) sering dinormalisasi agar setiap baris menjumlah ke 1:
\[ w_{ij}^{*} = \frac{w_{ij}}{\sum_j w_{ij}} \]
Hal ini memudahkan interpretasi dan analisis lanjutan.
1.3.4 Contoh Sederhana
Misalkan ada 4 kelurahan dengan pola tetangga sebagai berikut:
- A bertetangga dengan B dan C
- B bertetangga dengan A dan D
- C bertetangga dengan A
- D bertetangga dengan B
Maka matriks bobot \(W\) (rook contiguity) adalah:
\[ W = \begin{bmatrix} 0 & 1 & 1 & 0 \\\\ 1 & 0 & 0 & 1 \\\\ 1 & 0 & 0 & 0 \\\\ 0 & 1 & 0 & 0 \end{bmatrix} \]
1.4 Peran Lokasi, Jarak, dan Bobot Spasial dalam Analisis
- Autokorelasi spasial (Moran’s I, Geary’s C) membutuhkan definisi bobot spasial.
- Model regresi spasial (SLM, SEM) menggunakan \(W\) untuk menangkap dependensi.
- Geostatistik (Kriging) menggunakan jarak antar titik untuk membangun variogram.
1.5 Implementasi di R
library(spdep)
library(sf)
# Contoh data shapefile
nc <- st_read(system.file("shape/nc.shp", package="sf"))
# Membuat matriks bobot Queen contiguity
nb <- poly2nb(nc, queen=TRUE)
W <- nb2mat(nb, style="W") # Normalisasi row-standardized
# Menampilkan matriks bobot
print(W[1:5, 1:5])2 Visualisasi Distribusi Penyakit
Spatial mapping adalah proses visualisasi data yang memiliki dimensi
geografis.
Tujuan utama dari mapping adalah:
- Menyajikan data dalam bentuk peta agar pola spasial mudah
dikenali.
- Membantu dalam pengambilan keputusan berbasis wilayah.
- Menjadi dasar untuk analisis spasial lanjutan, seperti autokorelasi dan pemodelan spasial.
Dalam materi ini kita akan menggunakan data kasus Diare di Kota Bandung untuk contoh pemetaan.
2.1 Persiapan Library dan Data
library(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)
# Arahkan ke folder kerja Anda
# Data kasus
Diare <- read.csv("Diare.csv", sep=";")
# Peta administratif Indonesia level kecamatan (gadm36 level 3)
Indo_Kec <- readRDS('gadm36_IDN_3_sp.rds')
# Ambil hanya Kota Bandung
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
plot(Bandung)
# Tambahkan id untuk join
Bandung$id <- c(1:30)
# Ubah ke format sf
Bandung_sf <- st_as_sf(Bandung)
# Gabungkan dengan data kasus
Bandung_merged <- Bandung_sf %>%
left_join(Diare, by = "id")2.2 Pemetaan Kasus (Continuous Mapping)
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare), color=NA) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Peta Kasus Diare Kota Bandung",
fill = "Jumlah Kasus")
Interpretasi: Warna merah menandakan kasus lebih tinggi, kuning lebih rendah.
2.3 Pemetaan Diskrit (Discrete Mapping)
# Interval kategori
breaks <- c(-Inf, 550, 800, 1200, Inf)
labels <- c("Very Low", "Low", "High", "Very High")
# Buat variabel diskrit
Bandung_merged$Diare_Discrete <- cut(Bandung_merged$Diare,
breaks = breaks,
labels = labels,
right = TRUE)
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare_Discrete), color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow",
"Low" = "orange",
"High" = "red",
"Very High" = "red3")) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Kategori Kasus Diare Kota Bandung",
fill = "Kategori")
Interpretasi: Kategori “High” dan “Very High” dapat mengindikasikan wilayah dengan risiko lebih besar.
2.4 Autokorelasi Spasial (Moran’s I)
Autokorelasi spasial digunakan untuk melihat apakah data memiliki pola keterkaitan antar wilayah tetangga.
# Matriks bobot spasial
row.names(Bandung) <- as.character(1:30)
W <- poly2nb(Bandung, row.names(Bandung), queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE)
WL <- nb2listw(W)
# Plot jaringan tetangga
CoordK <- coordinates(Bandung)
plot(Bandung, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Bandung), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
# Global Moran's I
Global_Moran <- moran.test(Diare$Diare, WL)
Global_Moran##
## Moran I test under randomisation
##
## data: Diare$Diare
## weights: WL
##
## Moran I statistic standard deviate = 0.67226, p-value = 0.2507
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.04045578 -0.03448276 0.01242630Interpretasi: Nilai Moran’s I menunjukkan apakah ada pola clustering (positif) atau sebaran acak.
2.5 Local Moran’s I (LISA)
Local Moran’s I digunakan untuk mengidentifikasi cluster lokal (misalnya high-high, low-low).
Local_Moran <- localmoran(Diare$Diare, WL)
# Ambil informasi quadran
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Bandung_sf$mean_values <- mean_values_char
# Plot hasil LISA
ggplot() +
geom_sf(data = Bandung_sf, aes(fill = mean_values)) +
scale_fill_manual(values = c("Low-Low" = "blue",
"High-Low" = "green",
"Low-High" = "yellow",
"High-High" = "red")) +
labs(title = "Local Moran's I Kota Bandung",
fill = "Cluster Type") +
theme_minimal()
Interpretasi:
High-High (Merah): cluster wilayah dengan kasus tinggi dikelilingi wilayah tinggi.
Low-Low (Biru): cluster wilayah dengan kasus rendah dikelilingi wilayah rendah.
Low-High / High-Low: wilayah dengan nilai berbeda dari tetangganya (potensi outlier spasial).
2.6 Kesimpulan
- Mapping spasial membantu visualisasi distribusi penyakit.
- Continuous mapping menunjukkan gradien kasus.
- Discrete mapping memberi kategori untuk interpretasi lebih jelas.
- Moran’s I (global dan lokal) mendeteksi adanya pola spasial (clustering).
3 Spatial Autocorrelation
Spatial autocorrelation mengukur derajat kemiripan nilai antar lokasi. Dua ukuran global yang klasik adalah Moran’s \(I\) (Moran, 1950) dan Geary’s \(C\) (Geary, 1954). Analisis lokal menggunakan Local Moran’s \(I\) (Anselin, 1995) dan Getis–Ord untuk mengidentifikasi hot spots dan cold spots.
3.1 Notasi umum
Misalkan terdapat \(N\) area dengan variabel minat \(x_i\). Rata-rata \(\bar x\). Matriks bobot spasial \(W=(w_{ij})\) dengan \(w_{ii}=0\) kecuali dinyatakan lain. Jumlah total bobot \(S_0=\sum_i \sum_j w_{ij}\).
3.2 Rumus Indeks
3.2.1 Moran’s \(I\) (Global)
\[ I = \frac{N}{S_0} \cdot \frac{\sum_{i=1}^{N}\sum_{j=1}^{N} w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{N}(x_i-\bar{x})^2}. \]
Interpretasi:
- \(I > 0\): Nilai mirip
(tinggi–tinggi atau rendah–rendah) cenderung berdekatan \(\;\to\;\) positive spatial
autocorrelation.
- \(I < 0\): Pola
checkerboard (tinggi dikelilingi rendah, atau
sebaliknya).
- \(I \approx 0\): Pola acak.
Karakteristik:
- Mengukur kemiripan secara global.
- Lebih sensitif terhadap covariance antar lokasi.
- Dipandang sebagai analog spasial dari koefisien korelasi Pearson.
3.2.2 Local Moran’s \(I_i\) (LISA)
\[ I_i = \frac{x_i - \bar{x}}{m_2} \sum_{j=1}^{N} w_{ij}(x_j-\bar{x}), \qquad m_2 = \frac{1}{N}\sum_{k=1}^N (x_k-\bar{x})^2. \]
Interpretasi:
- Nilai \(I_i > 0\): Unit spasial
\(i\) memiliki nilai yang mirip dengan
tetangganya
→ membentuk cluster High–High (tinggi dikelilingi tinggi) atau Low–Low (rendah dikelilingi rendah).
- Nilai \(I_i < 0\): Unit spasial
\(i\) berbeda dengan tetangganya
→ pola High–Low (tinggi dikelilingi rendah) atau Low–High (rendah dikelilingi tinggi).
- Nilai \(I_i \approx 0\): Tidak ada pola lokal yang signifikan (acak).
Karakteristik:
- Merupakan ukuran lokal dari autokorelasi
spasial.
- Digunakan untuk mengidentifikasi lokasi cluster dan
outlier spasial.
- Sering dipakai dalam analisis eksplorasi data spasial (ESDA) untuk
membuat peta LISA cluster map.
- Melengkapi Moran’s \(I\) global dengan menunjukkan dimana autokorelasi terjadi.
3.2.3 Getis–Ord \(G_i\) dan \(G_i^\ast\) — raw statistic
Definisi asli dari statistik lokal Getis–Ord tidak terstandar adalah rasio massa tetangga dengan massa total (Ord & Getis, 1995). Untuk suatu band jarak \(d\) (atau bobot ketetanggaan biner) dengan \(w_{ij}(d)\in\{0,1\}\):
Tanpa memasukkan lokasi \(i\) (local \(G_i\)): \[ G_i(d) = \frac{\sum_{j \ne i} w_{ij}(d)\, x_j}{\sum_{j \ne i} x_j}. \]
Memasukkan lokasi \(i\) (local \(G_i^*\)): \[ G_i^*(d) = \frac{\sum_{j=1}^{N} w_{ij}(d)\, x_j}{\sum_{j=1}^{N} x_j}, \qquad \text{dengan } w_{ii}(d) = 1. \]
Nilai besar \(G_i^*\) mengindikasikan hot spot (nilai tinggi dikelilingi tinggi), sedangkan kecil mengindikasikan cold spot.
3.2.3.1 Standarisasi ke Z-skor (umum di perangkat lunak)
Untuk inferensi, \(G_i^*\) biasanya distandarisasi menjadi \(z(G_i^*)\). Bentuk z-skor (yang umum dihitung oleh paket perangkat lunak) adalah:
\[ z(G_i^*) = \frac{\sum_j w_{ij} x_j - \bar{x}\sum_j w_{ij}} { s \sqrt{\dfrac{ N \sum_j w_{ij}^2 - \left( \sum_j w_{ij} \right)^2 }{N-1}} }, \qquad s = \sqrt{\frac{\sum_k (x_k - \bar{x})^2}{N}}. \]
Catatan: Rumus di atas adalah z-skor dari \(G_i^*\), bukan definisi raw \(G_i^*\).
Interpretasi:
- Nilai tinggi \(G_i^*\) \(\;\to\;\) hot spot
(wilayah dengan nilai tinggi dikelilingi nilai tinggi).
- Nilai rendah \(G_i^*\) \(\;\to\;\) cold spot (wilayah dengan nilai rendah dikelilingi rendah).
Karakteristik:
- Digunakan untuk mendeteksi cluster lokal (hot spot
analysis).
- Tidak mendeteksi pola checkerboard (tinggi–rendah) seperti
Moran’s \(I\).
- Sangat populer untuk peta hot spot di epidemiologi & kriminologi.
3.2.4 Geary’s \(C\) (Global)
\[ C = \frac{(N-1)\sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - x_j)^2}{2 S_0 \sum_{i=1}^N (x_i-\bar{x})^2}. \]
Interpretasi:
- \(C < 1\): Autokorelasi positif
(tetangga mirip).
- \(C > 1\): Autokorelasi negatif
(tetangga berlawanan).
- \(C = 1\): Tidak ada autokorelasi (acak).
Karakteristik:
- Lebih sensitif terhadap perbedaan nilai antar tetangga
(contrast).
- Fokus pada variabilitas lokal, bukan kovarians.
- Bisa melengkapi Moran’s \(I\) (kadang hasilnya berbeda).
Ringkasan Perbedaan
|–||–|–|–| | Moran’s \(I\) | Global | Kovarians antar lokasi | \(I \approx 0\) | Kemiripan global (mirip Pearson) | | Geary’s \(C\) | Global | Variasi/perbedaan antar tetangga | \(C = 1\) | Lebih peka terhadap contrast lokal | | Getis–Ord \(G_i^*\) | Lokal | Deteksi hot/cold spots (klaster nilai) | \(Z\)-skor \(\approx 0\) | Identifikasi lokasi cluster tinggi/rendah |
3.3 Paket dan Data Peta Bandung
library(sf)
library(sp)
library(spdep)
library(dplyr)
library(ggplot2)
library(mapview)
library(spData)# ==== Opsi A: menggunakan file RDS Anda ====
Indo_Kec <- readRDS("gadm36_IDN_3_sp.rds")
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung", ]
Bandung$id <- seq_len(nrow(Bandung))
Bandung_sf <- sf::st_as_sf(Bandung)
# Untuk keperluan *knit* tanpa data eksternal, siapkan grid mock 30 poligon (DEMO).
if(!exists("Bandung_sf")) {
bb <- sf::st_as_sfc(sf::st_bbox(c(xmin=107.55, ymin=-6.98, xmax=107.72,
ymax=-6.85), crs=4326))
grid <- sf::st_make_grid(bb, n=c(6,5)) |> sf::st_as_sf() |> dplyr::mutate(id = dplyr::row_number())
Bandung_sf <- grid
}
Bandung_sf <- sf::st_make_valid(Bandung_sf)
Bandung_sf$id <- if(!"id" %in% names(Bandung_sf)) seq_len(nrow(Bandung_sf)) else Bandung_sf$idggplot(Bandung_sf) + geom_sf(fill="grey90", color="white") +
labs(title="Kota Bandung - Unit (contoh)") + theme_minimal()
Sketsa batas kecamatan (mock atau Bandung sebenarnya)
3.4 Simulasi Data Kasus Diare (Bandung)
# Neighbor queen
nb <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = TRUE)
# Listw untuk analitik global/lokal (row-standar & biner)
lwW <- spdep::nb2listw(nb, style="W") # row-standardized
lwB <- spdep::nb2listw(nb, style="B") # binary (untuk raw G)
N <- nrow(Bandung_sf)
z <- rnorm(N)
lz <- spdep::lag.listw(lwW, z)
beta0 <- 1.2; beta1 <- 0.6; beta2 <- 1.0
lambda <- exp(beta0 + beta1*scale(z) + beta2*scale(lz))
pop <- round(runif(N, 3000, 8000))
cases <- rpois(N, lambda=lambda * (pop/10000))
Diare <- dplyr::tibble(id = Bandung_sf$id, pop = pop, cases = cases,
rate10k = (cases/pop)*10000)
Bandung_merged <- dplyr::left_join(Bandung_sf, Diare, by="id")ggplot(Bandung_merged) +
geom_sf(aes(fill = rate10k), color="white", size=0.2) +
scale_fill_viridis_c(option="magma") +
labs(title="Rate Diare (per 10.000) — Simulasi", fill="Rate") +
theme_minimal()
Rate Diare (per 10.000) — Simulasi
3.5 Global Moran’s I
moran_res <- spdep::moran.test(Bandung_merged$rate10k, lwW,
randomisation = TRUE, alternative = "two.sided")
moran_res##
## Moran I test under randomisation
##
## data: Bandung_merged$rate10k
## weights: lwW
##
## Moran I statistic standard deviate = 4.8496, p-value = 1.237e-06
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.50481782 -0.03448276 0.012366373.6 Geary’s \(C\)
geary_res <- spdep::geary.test(Bandung_merged$rate10k, lwW,
randomisation = TRUE, alternative = "two.sided")
geary_res##
## Geary C test under randomisation
##
## data: Bandung_merged$rate10k
## weights: lwW
##
## Geary C statistic standard deviate = 4.3461, p-value = 1.386e-05
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.47515619 1.00000000 0.014583273.7 Local Moran’s I (LISA)
x <- scale(Bandung_merged$rate10k)[,1]
lagx <- spdep::lag.listw(lwW, x)
lisa <- spdep::localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
x >= 0 & lagx >= 0 ~ "High-High",
x < 0 & lagx < 0 ~ "Low-Low",
x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)
Bandung_LISA <- dplyr::bind_cols(Bandung_merged, lisa_df) |>
dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))ggplot(Bandung_LISA) +
geom_sf(aes(fill = quad), color="white", size=0.2) +
scale_fill_manual(values=c(
"High-High"="#d73027","Low-Low"="#4575b4",
"High-Low (Outlier)"="#fdae61","Low-High (Outlier)"="#74add1",
"Not significant"="grey85"
)) +
labs(title="Local Moran's I (LISA)", fill="Kategori") +
theme_minimal()
Local Moran’s I (LISA) — Cluster & Outlier (p≤0.05)
3.8 Getis–Ord:
x_raw <- Bandung_merged$rate10k
sum_x <- sum(x_raw)
# Matriks bobot biner (tidak distandarisasi)
Wb <- spdep::listw2mat(lwB) # diag=0 secara default
# G_i (tanpa i)
num_G <- as.numeric(Wb %*% x_raw) # Σ_{j≠i} w_ij x_j
den_G <- (sum_x - x_raw) # Σ_{j≠i} x_j
G_raw <- num_G / den_G
# G_i^* (dengan i) — set w_ii = 1
Wb_star <- Wb; diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_raw) # Σ_j w_ij^* x_j
den_Gs <- sum_x # Σ_j x_j
G_star_raw <- num_Gs / den_Gs
# Z-skor (yang dihitung spdep::localG)
Gz <- spdep::localG(x_raw, listw = lwW) # z(G_i^*) dengan listw yang diberikan
Bandung_G <- dplyr::mutate(Bandung_merged,
G_raw = G_raw,
G_star_raw = G_star_raw,
z_Gistar = as.numeric(Gz),
hotcold = dplyr::case_when(
z_Gistar >= 1.96 ~ "Hot spot (p≈0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p≈0.05)",
TRUE ~ "Not significant"
))
summary(dplyr::select(Bandung_G, G_raw, G_star_raw, z_Gistar))## G_raw G_star_raw z_Gistar geometry
## Min. :0.00000 Min. :0.01069 Min. :-1.7571 MULTIPOLYGON :30
## 1st Qu.:0.05115 1st Qu.:0.07107 1st Qu.:-1.1364 epsg:NA : 0
## Median :0.09574 Median :0.11959 Median :-0.9334 +proj=long...: 0
## Mean :0.13832 Mean :0.16419 Mean :-0.2353
## 3rd Qu.:0.16788 3rd Qu.:0.18923 3rd Qu.: 0.1827
## Max. :0.48942 Max. :0.53346 Max. : 4.1541ggplot(Bandung_G) +
geom_sf(aes(fill = G_star_raw), color="white", size=0.2) +
scale_fill_viridis_c() +
labs(title="Raw Getis–Ord G* (proporsi massa tetangga)", fill="G*_raw") +
theme_minimal()
ggplot(Bandung_G) +
geom_sf(aes(fill = hotcold), color="white", size=0.2) +
scale_fill_manual(values=c("Hot spot (p≈0.05)"="#b2182b",
"Cold spot (p≈0.05)"="#2166ac",
"Not significant"="grey85")) +
labs(title="Getis–Ord Gi* — Hot/Cold Spots (z-skor)", fill=NULL) +
theme_minimal()
3.9 Sensitivitas Bobot Spasial
# Rook: buat ulang neighbor dengan queen = FALSE
nb_rook <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = FALSE)
lw_rook <- spdep::nb2listw(nb_rook, style="W")
# kNN (k = 4)
coords <- sf::st_coordinates(sf::st_centroid(Bandung_sf))
nb_knn <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
lw_knn <- spdep::nb2listw(nb_knn, style="W")
# Bandingkan rata-rata z(Gi*) untuk tiga bobot
c(
mean_z_queen = mean(spdep::localG(Bandung_merged$rate10k, lwW)),
mean_z_rook = mean(spdep::localG(Bandung_merged$rate10k, lw_rook)),
mean_z_knn4 = mean(spdep::localG(Bandung_merged$rate10k, lw_knn))
)## mean_z_queen mean_z_rook mean_z_knn4
## -0.2353182 -0.2481379 -0.19051733.10 Interpretasi
- Raw \(G_i\)/\(G_i^*\) adalah ukuran proporsi massa tetangga dari variabel \(x\). Untuk pengambilan keputusan statistik, gunakan z-skor \(z(G_i^*)\).
- Hasil peka terhadap definisi bobot (queen/rook/kNN).
3.11 Lampiran: Ekspor
readr::write_csv(dplyr::as_tibble(Diare), "diare_simulasi_bandung.csv")
sf::st_write(Bandung_LISA, "bandung_lisa.geojson", delete_dsn=TRUE)
sf::st_write(Bandung_G, "bandung_getis.geojson", delete_dsn=TRUE)Referensi
- Moran, P. A. P. (1950). Notes on Continuous Stochastic
Phenomena. Biometrika.
- Geary, R. C. (1954). The Contiguity Ratio and Statistical
Mapping. The Incorporated Statistician.
- Anselin, L. (1995). Local Indicators of Spatial
Association—LISA. Geographical Analysis.
- Getis, A., & Ord, J. K. (1992, 1995). The Analysis of
Spatial Association by Use of Distance Statistics / Local Spatial
Statistics. Geographical Analysis.
- Moraga, P. (2023). Spatial Statistics for Data Science — Bab Spatial Autocorrelation. (https://www.paulamoraga.com/book-spatial/spatial-autocorrelation.html)
3.12 Tugas 1
Mahasiswa memahami konsep dasar autokorelasi spasial, mampu menghitung ukuran global dan lokal (Moran’s I, Geary’s C, Local Moran’s I, Getis–Ord), serta menginterpretasikan hasil dalam konteks nyata.
Bagian A. Konsep Teori
- Jelaskan dengan kata-kata Anda apa yang dimaksud dengan autokorelasi spasial positif dan autokorelasi spasial negatif.
- Berikan contoh sederhana dari fenomena yang bisa menunjukkan masing-masing pola tersebut.
- Tuliskan rumus matematis dari:
- Moran’s \(I\)
- Geary’s \(C\)
- Local Moran’s \(I_i\)
- Getis–Ord \(G_i\) dan \(G_i^*\)
- Jelaskan perbedaan utama antara ukuran global dan lokal.
Bagian B. Analisis Data (Simulasi)
- Gunakan data spasial kecamatan di Kota Bandung (atau gunakan grid simulasi jika data asli tidak tersedia).
- Buat data simulasi kasus penyakit menular (misalnya diare per 10.000 penduduk) untuk 30 kecamatan di Kota Bandung. Gunakan distribusi Poisson dengan rata-rata berbeda antar kecamatan.
- Buat peta choropleth dari data simulasi tersebut.
- Apa pola spasial yang terlihat secara visual?
Bagian C. Pengukuran Autokorelasi
- Hitung Moran’s I untuk data simulasi kasus diare di
Kota Bandung.
- Berapa nilai Moran’s I?
- Apakah signifikan secara statistik (uji permutasi)?
- Apa artinya bagi pola spasial penyakit?
- Hitung Geary’s C.
- Bagaimana perbandingannya dengan Moran’s I?
- Jelaskan perbedaan sensitivitas kedua ukuran ini.
- Hitung Local Moran’s I (LISA).
- Identifikasi kecamatan yang masuk kategori High-High, Low-Low, High-Low, dan Low-High.
- Buat peta cluster LISA.
- Apa interpretasi hasil ini untuk kasus penyakit menular?
- Hitung Getis–Ord \(G_i^*\).
- Tentukan kecamatan yang termasuk hot spot dan cold spot.
- Bandingkan hasilnya dengan peta LISA.
- Apakah ada perbedaan wilayah yang ditandai sebagai klaster signifikan?
Bagian D. Diskusi Kritis
- Bagaimana hasil analisis autokorelasi spasial bisa membantu dinas kesehatan dalam menyusun strategi pencegahan dan intervensi penyakit menular di Kota Bandung?
- Sebutkan keterbatasan dari analisis autokorelasi spasial, misalnya
terkait dengan:
- MAUP (Modifiable Areal Unit Problem)
- Ukuran bobot spasial (rook, queen, k-nearest neighbors)
- Masalah multiple testing pada analisis lokal
Penilaian
- Konsep Teori (20%) → kejelasan penjelasan konsep
dan rumus
- Analisis Data (30%) → kode R untuk simulasi data
dan peta choropleth
- Perhitungan Indeks (30%) → hasil Moran’s I, Geary’s
C, LISA, Getis–Ord beserta interpretasi
- Diskusi Kritis (20%) → refleksi dan keterbatasan
analisis
# Spatial Econometrics Spatial econometrics mempelajari data yang terindeks lokasi (wilayah, grid, titik koordinat) yang sering menunjukkan spatial dependence (unit saling memengaruhi) dan spatial heterogeneity (hubungan berbeda antar tempat).
Mengapa OLS biasa tidak cukup?
- Autokorelasi spasial pada residual melanggar asumsi independensi; standar error bias.
- Spillovers: perubahan pada suatu wilayah memengaruhi wilayah lain.
- Omitted variables yang berpola spasial menimbulkan error spatial correlation.
Notasi umum:
- \(y \in \mathbb{R}^{N}\), \(X \in \mathbb{R}^{N \times K}\), \(W \in \mathbb{R}^{N \times N}\) (baris-terstandar), \(I\) identitas.
- \(Wy\) atau \(WX\) adalah lag spasial terhadap \(y\) atau \(X\).
3.13 Data & Bobot Spasial (W)
Kita gunakan dua sumber data agar komprehensif:
- Simulasi grid (replikasi mudah).
- Data nyata Columbus (kriminalitas &
variabel sosio-ekonomi) dari paket
spData.
# Install jika perlu:
# install.packages(c("spdep", "spatialreg", "Matrix", "spData", "sp"))
library(spdep)
library(spatialreg)
library(Matrix)
library(spData) # untuk data columbus
library(sp)
set.seed(123)3.14 Simulasi grid & W (queen contiguity)
nrow_grid <- 15; ncol_grid <- 15
N <- nrow_grid * ncol_grid
nb <- cell2nb(nrow_grid, ncol_grid, type = "queen")
lw <- nb2listw(nb, style = "W", zero.policy = TRUE) # W row-standardized
W <- listw2mat(lw)
W_mat <- Matrix(W, sparse = TRUE)
# Kovariat & outcome untuk contoh
p <- 3
X <- cbind(1, matrix(rnorm(N*(p-1)), N, p-1))
colnames(X) <- c("(Intercept)", paste0("x", 1:(p-1)))
beta_true <- c(2, 1.0, -0.8)
rho_true <- 0.4; sigma_true <- 1.0
I_N <- Diagonal(N)
eps <- rnorm(N, sd = sigma_true)
y_mean <- as.numeric(X %*% beta_true)
y <- as.numeric(solve(I_N - rho_true * W_mat, y_mean + eps)) # SAR DGP
dat_grid <- data.frame(y = y, X[, -1, drop = FALSE])3.15 Data nyata: Columbus (kriminalitas)
data(columbus, package = "spData")
# 'columbus' umumnya hadir sebagai SpatialPolygonsDataFrame di beberapa versi spData/spdep.
# Jika tidak tersedia sebagai objek 'Spatial', kita buat poligon sederhana dari centroid via k-NN agar contoh tetap jalan.
# Jika columbus berupa 'Spatial*DataFrame'
if (inherits(columbus, "Spatial")) {
coords <- coordinates(columbus)
nb_c <- poly2nb(columbus) # ROOK dari poligon
} else {
# Fallback: gunakan koordinat X,Y lalu buat nb KNN
coords <- as.matrix(columbus[, c("X", "Y")])
knn <- knearneigh(coords, k = 5)
nb_c <- knn2nb(knn)
}
lw_c <- nb2listw(nb_c, style = "W")
d_col <- as.data.frame(columbus)
d_col$CRIME <- d_col$CRIME
d_col$INC <- d_col$INC
d_col$HOVAL <- d_col$HOVALCatatan: Perbedaan versi paket dapat mengubah kelas objek
columbus. Blok di atas membuat contoh tetap replikabel.
3.16 SLX — Spatial Lag of X
Ide dasar. Hanya spillover dari
kovariat tetangga dimodelkan, bukan dari \(y\). Tidak ada umpan-balik antar
respon.
Kapan digunakan. Saat kita percaya bahwa karakteristik
tetangga (mis. pendapatan, fasilitas publik) memengaruhi wilayah, tetapi
hasil/respon tidak saling menular.
Contoh aplikasi. Kemiskinan kota dipengaruhi oleh
pendidikan di kota tetangga; tingkat vaksinasi tetangga
memengaruhi risiko lokal meski prevalensi lokal tidak “menular”
antary.
Rumus. \[ y = X\beta + W X\theta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I). \]
R di data grid (simulasi).
slx_grid <- lmSLX(y ~ ., data = dat_grid, listw = lw, Durbin = TRUE)
summary(slx_grid)##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.351e+00 6.823e-02 4.911e+01 1.528e-120
## x1 9.767e-01 7.262e-02 1.345e+01 1.712e-30
## x2 -8.080e-01 6.825e-02 -1.184e+01 2.424e-25
## lag.x1 6.620e-01 2.003e-01 3.305e+00 1.108e-03
## lag.x2 -2.933e-02 1.845e-01 -1.590e-01 8.738e-01R di data Columbus (aplikasi nyata). Kita ingin melihat pengaruh INC (pendapatan) dan HOVAL (nilai rumah) lokal dan tetangga terhadap CRIME.
slx_col <- lmSLX(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, Durbin = TRUE)
summary(slx_col)##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.609e+01 6.168e+00 1.234e+01 7.058e-16
## INC -1.255e+00 3.857e-01 -3.253e+00 2.194e-03
## HOVAL -2.638e-01 1.031e-01 -2.559e+00 1.403e-02
## lag.INC -5.277e-01 6.667e-01 -7.916e-01 4.328e-01
## lag.HOVAL -1.481e-01 2.594e-01 -5.709e-01 5.710e-01Interpretasi. Koefisien pada \(X\) = efek langsung; pada \(WX\) = spillover. Tidak
perlu impacts() karena tidak ada lag \(y\).
3.17 SAR — Spatial Autoregressive (lag pada y)
Ide dasar. Respon suatu wilayah dipengaruhi respon
tetangga: umpan balik & difusi.
Kapan digunakan. Harga rumah, kriminalitas, adopsi
teknologi, pertumbuhan ekonomi regional.
Rumus.
\[ y = \rho Wy + X\beta + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, \sigma^2 I). \] Dengan \(A(\rho)=I-\rho W\), model setara: \(A(\rho)y = X\beta + \epsilon\).
Likelihoood ML (ringkas). \[
\ell(\beta,\rho,\sigma^2) = \log|A(\rho)| - \frac{N}{2}\log\sigma^2
- \frac{1}{2\sigma^2}\big(A(\rho)y - X\beta\big)^\top\big(A(\rho)y -
X\beta\big).
\] Likelihood terkonsentrasi: untuk setiap \(\rho\),
\(\hat{\beta}(\rho)=(X^\top X)^{-1}X^\top
A(\rho)y\),
\(\hat{\sigma}^2(\rho)=\frac{1}{N}\|A(\rho)y -
X\hat{\beta}(\rho)\|^2\).
Maksimalkan \(\ell_c(\rho)=\log|A(\rho)| -
\tfrac{N}{2}\log\hat{\sigma}^2(\rho)\).
R (grid & Columbus).
sar_grid <- lagsarlm(y ~ ., data = dat_grid, listw = lw, method = "eigen")
sar_col <- lagsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, method = "eigen")
summary(sar_grid)##
## Call:lagsarlm(formula = y ~ ., data = dat_grid, listw = lw, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.795576 -0.616588 0.048715 0.703882 2.411899
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.928892 0.257210 7.4993 6.417e-14
## x1 0.952970 0.068100 13.9936 < 2.2e-16
## x2 -0.805996 0.063997 -12.5943 < 2.2e-16
##
## Rho: 0.42618, LR test value: 27.234, p-value: 1.8027e-07
## Asymptotic standard error: 0.074348
## z-value: 5.7322, p-value: 9.9132e-09
## Wald statistic: 32.858, p-value: 9.9132e-09
##
## Log likelihood: -313.5954 for lag model
## ML residual variance (sigma squared): 0.92292, (sigma: 0.96069)
## Number of observations: 225
## Number of parameters estimated: 5
## AIC: 637.19, (AIC for lm: 662.42)
## LM test for residual autocorrelation
## test value: 0.80142, p-value: 0.37067summary(sar_col)##
## Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.87401 -6.05795 0.11128 5.94981 22.51885
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 40.745935 7.181891 5.6734 1.4e-08
## INC -0.910833 0.301872 -3.0173 0.002551
## HOVAL -0.260186 0.084861 -3.0660 0.002169
##
## Rho: 0.46489, LR test value: 14.121, p-value: 0.00017145
## Asymptotic standard error: 0.11549
## z-value: 4.0255, p-value: 5.6857e-05
## Wald statistic: 16.205, p-value: 5.6857e-05
##
## Log likelihood: -180.3169 for lag model
## ML residual variance (sigma squared): 88.291, (sigma: 9.3963)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 370.63, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 3.1695, p-value: 0.075025Dampak (direct/indirect/total). Pada SAR, koefisien
\(\beta\) bukan efek
marjinal jangka panjang. Gunakan impacts().
imp_grid <- impacts(sar_grid, listw = lw, R = 500)
imp_col <- impacts(sar_col, listw = lw_c, R = 500)
summary(imp_grid, zstats = TRUE)## Impact measures (lag, exact):
## Direct Indirect Total
## x1 0.9843937 0.6763428 1.660737
## x2 -0.8325739 -0.5720327 -1.404607
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 0.9874 0.07134 0.003190 0.003190
## x2 -0.8331 0.06504 0.002909 0.002342
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.8558 0.9385 0.9847 1.0353 1.1361
## x2 -0.9566 -0.8779 -0.8329 -0.7877 -0.7047
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 0.6996 0.2193 0.009805 0.009805
## x2 -0.5891 0.1835 0.008207 0.008207
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.3184 0.5526 0.6672 0.8417 1.1531
## x2 -1.0069 -0.7037 -0.5585 -0.4692 -0.2685
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 1.687 0.2577 0.011526 0.011526
## x2 -1.422 0.2152 0.009623 0.009623
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 1.215 1.522 1.663 1.857 2.223
## x2 -1.907 -1.551 -1.388 -1.283 -1.052
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## x1 0.07133920 0.2192545 0.2577311
## x2 0.06504308 0.1835187 0.2151677
##
## Simulated z-values:
## Direct Indirect Total
## x1 13.84085 3.190707 6.545473
## x2 -12.80801 -3.209856 -6.609454
##
## Simulated p-values:
## Direct Indirect Total
## x1 < 2.22e-16 0.0014192 5.9307e-11
## x2 < 2.22e-16 0.0013280 3.8574e-11summary(imp_col, zstats = TRUE)## Impact measures (lag, exact):
## Direct Indirect Total
## INC -0.9542578 -0.7478805 -1.7021382
## HOVAL -0.2725909 -0.2136376 -0.4862285
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -0.9554 0.31176 0.013942 0.013942
## HOVAL -0.2747 0.08672 0.003878 0.003878
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.6300 -1.1643 -0.9529 -0.7563 -0.4013
## HOVAL -0.4524 -0.3315 -0.2730 -0.2188 -0.1144
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -0.7911 0.3795 0.016972 0.016972
## HOVAL -0.2408 0.1417 0.006338 0.006909
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.6974 -0.9566 -0.7432 -0.5333 -0.27526
## HOVAL -0.6088 -0.3037 -0.2099 -0.1429 -0.05387
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -1.7465 0.5821 0.026031 0.026031
## HOVAL -0.5155 0.2061 0.009218 0.009893
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -3.084 -2.0691 -1.7338 -1.3510 -0.7374
## HOVAL -1.000 -0.6293 -0.4869 -0.3706 -0.1900
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## INC 0.3117554 0.3795156 0.5820792
## HOVAL 0.0867234 0.1417271 0.2061268
##
## Simulated z-values:
## Direct Indirect Total
## INC -3.064683 -2.084436 -3.000464
## HOVAL -3.167153 -1.699260 -2.500875
##
## Simulated p-values:
## Direct Indirect Total
## INC 0.0021790 0.037121 0.0026957
## HOVAL 0.0015394 0.089270 0.01238873.18 SEM — Spatial Error Model
Ide dasar. Tidak ada lag pada \(y\), tapi error saling
berkorelasi: variabel penting terabaikan berpola spasial.
Kapan digunakan. Jika Moran’s I pada residual OLS
signifikan, dan teori tidak mendukung endogenous interaction
pada \(y\).
Rumus.
\[ \begin{aligned} y &= X\beta + u,\\ u &= \lambda Wu + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, \sigma^2 I). \end{aligned} \] Dengan \(B(\lambda)=I-\lambda W\), log-likelihood: \[ \ell(\beta,\lambda,\sigma^2)=\log|B(\lambda)| - \tfrac{N}{2}\log\sigma^2 - \tfrac{1}{2\sigma^2}\,\|B(\lambda)(y - X\beta)\|^2. \]
R (grid & Columbus).
sem_grid <- errorsarlm(y ~ ., data = dat_grid, listw = lw, method = "eigen")
sem_col <- errorsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, method = "eigen")
summary(sem_grid)##
## Call:errorsarlm(formula = y ~ ., data = dat_grid, listw = lw, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.736741 -0.650379 0.030066 0.702579 2.270119
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.359068 0.125819 26.698 < 2.2e-16
## x1 0.908493 0.067459 13.467 < 2.2e-16
## x2 -0.795855 0.063042 -12.624 < 2.2e-16
##
## Lambda: 0.48915, LR test value: 23.206, p-value: 1.4554e-06
## Asymptotic standard error: 0.088634
## z-value: 5.5187, p-value: 3.4148e-08
## Wald statistic: 30.456, p-value: 3.4148e-08
##
## Log likelihood: -315.6094 for error model
## ML residual variance (sigma squared): 0.92952, (sigma: 0.96412)
## Number of observations: 225
## Number of parameters estimated: 5
## AIC: 641.22, (AIC for lm: 662.42)summary(sem_col)##
## Call:errorsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.19485 -5.49306 -0.10805 5.81121 21.92606
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 55.830269 5.949030 9.3848 < 2.2e-16
## INC -1.014431 0.306651 -3.3081 0.0009393
## HOVAL -0.246757 0.081733 -3.0191 0.0025356
##
## Lambda: 0.6788, LR test value: 15.814, p-value: 6.9889e-05
## Asymptotic standard error: 0.11407
## z-value: 5.9509, p-value: 2.6672e-09
## Wald statistic: 35.413, p-value: 2.6672e-09
##
## Log likelihood: -179.4703 for error model
## ML residual variance (sigma squared): 80.04, (sigma: 8.9465)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 368.94, (AIC for lm: 382.75)3.19 SDM — Spatial Durbin Model
Ide dasar. SAR + SLX: respon dipengaruhi \(Wy\) dan juga \(WX\).
Kapan digunakan. Ada interaksi endogen pada \(y\) dan spillover dari
kovariat; model kaya untuk menilai kebijakan yang menular antar
wilayah.
Rumus. \[
y = \rho Wy + X\beta + WX\theta + \epsilon.
\]
R (grid & Columbus) + impacts.
sdm_grid <- lagsarlm(y ~ ., data = dat_grid, listw = lw, Durbin = TRUE, method = "eigen")
sdm_col <- lagsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, Durbin = TRUE, method = "eigen")
summary(sdm_grid)##
## Call:lagsarlm(formula = y ~ ., data = dat_grid, listw = lw, Durbin = TRUE,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7233396 -0.6444892 0.0073549 0.6571302 2.3417941
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.841731 0.313645 5.8720 4.305e-09
## x1 0.942181 0.067925 13.8709 < 2.2e-16
## x2 -0.807055 0.063488 -12.7120 < 2.2e-16
## lag.x1 0.104119 0.217525 0.4787 0.63218
## lag.x2 0.313018 0.186639 1.6771 0.09352
##
## Rho: 0.45084, LR test value: 19.891, p-value: 8.1986e-06
## Asymptotic standard error: 0.091615
## z-value: 4.921, p-value: 8.6086e-07
## Wald statistic: 24.217, p-value: 8.6086e-07
##
## Log likelihood: -311.8131 for mixed model
## ML residual variance (sigma squared): 0.90482, (sigma: 0.95122)
## Number of observations: 225
## Number of parameters estimated: 7
## AIC: 637.63, (AIC for lm: 655.52)
## LM test for residual autocorrelation
## test value: 0.14516, p-value: 0.70321summary(sdm_col)##
## Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## Durbin = TRUE, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.20282 -6.19526 -0.28693 4.85019 22.19628
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.0866190 9.0368595 2.5547 0.0106274
## INC -1.0149657 0.3065914 -3.3105 0.0009314
## HOVAL -0.2542240 0.0820306 -3.0991 0.0019409
## lag.INC 0.8360698 0.5429754 1.5398 0.1236108
## lag.HOVAL 0.0051799 0.2116339 0.0245 0.9804731
##
## Rho: 0.6519, LR test value: 12.594, p-value: 0.00038703
## Asymptotic standard error: 0.12077
## z-value: 5.3979, p-value: 6.7437e-08
## Wald statistic: 29.137, p-value: 6.7437e-08
##
## Log likelihood: -179.0855 for mixed model
## ML residual variance (sigma squared): 79.639, (sigma: 8.9241)
## Number of observations: 49
## Number of parameters estimated: 7
## AIC: 372.17, (AIC for lm: 382.76)
## LM test for residual autocorrelation
## test value: 4.0639, p-value: 0.04381summary(impacts(sdm_grid, listw = lw, R = 500), zstats = TRUE)## Impact measures (mixed, exact):
## Direct Indirect Total
## x1 0.9864373 0.91884630 1.9052836
## x2 -0.8113090 -0.08831963 -0.8996287
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 0.9869 0.07527 0.003366 0.003366
## x2 -0.8167 0.06938 0.003103 0.003103
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.8538 0.9349 0.9853 1.0358 1.1363
## x2 -0.9505 -0.8625 -0.8190 -0.7696 -0.6794
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 0.93844 0.3603 0.01611 0.01611
## x2 -0.08369 0.3017 0.01349 0.01370
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.2867 0.6897 0.9105 1.1758 1.6537
## x2 -0.6199 -0.2913 -0.0753 0.1293 0.4956
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 1.9253 0.3958 0.01770 0.0177
## x2 -0.9004 0.3302 0.01477 0.0146
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 1.220 1.651 1.9088 2.1844 2.6696
## x2 -1.523 -1.116 -0.9059 -0.6726 -0.2693
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## x1 0.07527155 0.3602884 0.3958025
## x2 0.06938481 0.3016917 0.3302043
##
## Simulated z-values:
## Direct Indirect Total
## x1 13.11127 2.604704 4.864421
## x2 -11.77091 -0.277413 -2.726844
##
## Simulated p-values:
## Direct Indirect Total
## x1 < 2.22e-16 0.0091954 1.1479e-06
## x2 < 2.22e-16 0.7814630 0.0063943summary(impacts(sdm_col, listw = lw_c, R = 500), zstats = TRUE)## Impact measures (mixed, exact):
## Direct Indirect Total
## INC -0.9819288 0.4680021 -0.5139266
## HOVAL -0.2846355 -0.4308104 -0.7154459
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:498
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 498
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -0.9764 0.33158 0.01486 0.01392
## HOVAL -0.2928 0.09729 0.00436 0.00436
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.6035 -1.2049 -0.9875 -0.7494 -0.3363
## HOVAL -0.4857 -0.3604 -0.2915 -0.2305 -0.1097
##
## ========================================================
## Indirect:
##
## Iterations = 1:498
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 498
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC 0.6503 2.0586 0.09225 0.09225
## HOVAL -0.5220 0.8763 0.03927 0.03927
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -2.771 -0.3451 0.6053 1.44029 4.9712
## HOVAL -2.386 -0.8556 -0.4433 -0.06552 0.8535
##
## ========================================================
## Total:
##
## Iterations = 1:498
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 498
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -0.3261 2.208 0.09895 0.09895
## HOVAL -0.8148 0.934 0.04185 0.04185
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -3.963 -1.404 -0.3819 0.5853 4.3460
## HOVAL -2.801 -1.186 -0.7325 -0.3275 0.7446
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## INC 0.33158007 2.0585525 2.2082167
## HOVAL 0.09729073 0.8762607 0.9340098
##
## Simulated z-values:
## Direct Indirect Total
## INC -2.944688 0.3159212 -0.1476574
## HOVAL -3.009349 -0.5956923 -0.8723286
##
## Simulated p-values:
## Direct Indirect Total
## INC 0.0032328 0.75206 0.88261
## HOVAL 0.0026181 0.55138 0.383033.20 SDEM — Spatial Durbin Error Model
Ide dasar. Tidak ada lag \(y\), tetapi ada WX dan
error SEM.
Kapan digunakan. Respon tidak saling menular, namun
kovariat tetangga dan faktor tak terukur (berpola spasial) memengaruhi
respon.
Rumus. \[
\begin{aligned}
y &= X\beta + WX\theta + u,\\
u &= \lambda Wu + \epsilon.
\end{aligned}
\]
R (grid & Columbus).
sdem_grid <- errorsarlm(y ~ ., data = dat_grid, listw = lw, Durbin = TRUE, method = "eigen")
sdem_col <- errorsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, Durbin = TRUE, method = "eigen")
summary(sdem_grid)##
## Call:errorsarlm(formula = y ~ ., data = dat_grid, listw = lw, Durbin = TRUE,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.75174930 -0.63204781 -0.00014401 0.69128917 2.40517545
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.353648 0.115372 29.0682 < 2.2e-16
## x1 0.973797 0.070290 13.8539 < 2.2e-16
## x2 -0.809610 0.066648 -12.1475 < 2.2e-16
## lag.x1 0.608037 0.222470 2.7331 0.006274
## lag.x2 -0.058292 0.205883 -0.2831 0.777075
##
## Lambda: 0.44962, LR test value: 19.553, p-value: 9.7822e-06
## Asymptotic standard error: 0.09252
## z-value: 4.8596, p-value: 1.176e-06
## Wald statistic: 23.616, p-value: 1.176e-06
##
## Log likelihood: -311.9819 for error model
## ML residual variance (sigma squared): 0.90636, (sigma: 0.95203)
## Number of observations: 225
## Number of parameters estimated: 7
## AIC: 637.96, (AIC for lm: 655.52)summary(sdem_col)##
## Call:errorsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## Durbin = TRUE, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -31.15319 -5.83223 -0.95878 5.18274 22.76056
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 53.422602 12.940877 4.1282 3.656e-05
## INC -0.934886 0.308678 -3.0287 0.002456
## HOVAL -0.275318 0.085404 -3.2237 0.001265
## lag.INC 0.821946 0.797123 1.0311 0.302475
## lag.HOVAL -0.277941 0.280238 -0.9918 0.321294
##
## Lambda: 0.72231, LR test value: 12.918, p-value: 0.00032549
## Asymptotic standard error: 0.10205
## z-value: 7.0783, p-value: 1.4597e-12
## Wald statistic: 50.102, p-value: 1.4596e-12
##
## Log likelihood: -178.9235 for error model
## ML residual variance (sigma squared): 76.757, (sigma: 8.7611)
## Number of observations: 49
## Number of parameters estimated: 7
## AIC: 371.85, (AIC for lm: 382.76)3.21 SAC / SARAR — Spatial Autoregressive Combined
Ide dasar. Gabungan SAR + SEM: ada
Wy dan error SEM.
Kapan digunakan. Ada interaksi endogen dan juga omitted
variables berpola spasial.
Rumus. \[
\begin{aligned}
y &= \rho Wy + X\beta + u,\\
u &= \lambda Wu + \epsilon.
\end{aligned}
\]
R (grid & Columbus).
sac_grid <- sacsarlm(y ~ ., data = dat_grid, listw = lw, method = "eigen")
sac_col <- sacsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, method = "eigen")
summary(sac_grid)##
## Call:sacsarlm(formula = y ~ ., data = dat_grid, listw = lw, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.79687 -0.61385 0.03868 0.68742 2.39961
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.216086 0.464593 4.770 1.843e-06
## x1 0.947340 0.068164 13.898 < 2.2e-16
## x2 -0.814330 0.064009 -12.722 < 2.2e-16
##
## Rho: 0.34052
## Asymptotic standard error: 0.13667
## z-value: 2.4915, p-value: 0.012719
## Lambda: 0.17137
## Asymptotic standard error: 0.19116
## z-value: 0.89643, p-value: 0.37002
##
## LR test value: 28.029, p-value: 8.196e-07
##
## Log likelihood: -313.1979 for sac model
## ML residual variance (sigma squared): 0.92629, (sigma: 0.96244)
## Number of observations: 225
## Number of parameters estimated: 6
## AIC: 638.4, (AIC for lm: 662.42)summary(sac_col)##
## Call:sacsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.51350 -5.53228 -0.27094 5.52658 22.13005
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 48.866523 9.908562 4.9317 8.15e-07
## INC -1.032095 0.307752 -3.3537 0.0007975
## HOVAL -0.254910 0.083354 -3.0582 0.0022270
##
## Rho: 0.26441
## Asymptotic standard error: 0.23835
## z-value: 1.1093, p-value: 0.26729
## Lambda: 0.47403
## Asymptotic standard error: 0.25075
## z-value: 1.8904, p-value: 0.058701
##
## LR test value: 16.879, p-value: 0.0002161
##
## Log likelihood: -178.9375 for sac model
## ML residual variance (sigma squared): 82.309, (sigma: 9.0724)
## Number of observations: 49
## Number of parameters estimated: 6
## AIC: 369.87, (AIC for lm: 382.75)3.22 GNS — General Nested Spatial Model
Ide dasar. Model paling umum: gabungkan SAR + SLX +
SEM (sering disebut SAC + Durbin).
Kapan digunakan. Saat ingin kerangka paling kaya,
kemudian lakukan reduksi berdasarkan signifikansi parameter/teori.
Rumus. \[
\begin{aligned}
y &= \rho Wy + X\beta + WX\theta + u,\\
u &= \lambda Wu + \epsilon.
\end{aligned}
\]
R (grid & Columbus).
gns_grid <- sacsarlm(y ~ ., data = dat_grid, listw = lw, Durbin = TRUE, method = "eigen")
gns_col <- sacsarlm(CRIME ~ INC + HOVAL, data = d_col, listw = lw_c, Durbin = TRUE, method = "eigen")
summary(gns_grid)##
## Call:sacsarlm(formula = y ~ ., data = dat_grid, listw = lw, Durbin = TRUE,
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7527333 -0.6451218 0.0066546 0.6795943 2.3712426
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.217002 2.193647 1.0106 0.3122
## x1 0.950920 0.084273 11.2838 <2e-16
## x2 -0.807459 0.064255 -12.5665 <2e-16
## lag.x1 0.235516 0.798014 0.2951 0.7679
## lag.x2 0.218957 0.550015 0.3981 0.6906
##
## Rho: 0.33888
## Asymptotic standard error: 0.65396
## z-value: 0.51819, p-value: 0.60433
## Lambda: 0.14675
## Asymptotic standard error: 0.74906
## z-value: 0.19591, p-value: 0.84468
##
## LR test value: 30.899, p-value: 3.2095e-06
##
## Log likelihood: -311.7627 for sacmixed model
## ML residual variance (sigma squared): 0.91581, (sigma: 0.95698)
## Number of observations: 225
## Number of parameters estimated: 8
## AIC: 639.53, (AIC for lm: 662.42)summary(gns_col)##
## Call:sacsarlm(formula = CRIME ~ INC + HOVAL, data = d_col, listw = lw_c,
## Durbin = TRUE, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.19511 -5.64820 -0.33221 4.85169 22.57305
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 36.490647 118.632989 0.3076 0.7583922
## INC -1.017589 0.305562 -3.3302 0.0008678
## HOVAL -0.260051 0.099702 -2.6083 0.0090996
## lag.INC 0.817395 1.880142 0.4348 0.6637427
## lag.HOVAL -0.113744 0.735706 -0.1546 0.8771325
##
## Rho: 0.39954
## Asymptotic standard error: 2.0013
## z-value: 0.19964, p-value: 0.84176
## Lambda: 0.45966
## Asymptotic standard error: 1.8707
## z-value: 0.24572, p-value: 0.8059
##
## LR test value: 17.985, p-value: 0.0012423
##
## Log likelihood: -178.3846 for sacmixed model
## ML residual variance (sigma squared): 79.321, (sigma: 8.9062)
## Number of observations: 49
## Number of parameters estimated: 8
## AIC: 372.77, (AIC for lm: 382.75)summary(impacts(gns_grid, listw = lw, R = 300), zstats = TRUE)## Impact measures (sacmixed, exact):
## Direct Indirect Total
## x1 0.9829227 0.81165277 1.7945755
## x2 -0.8105957 -0.07955803 -0.8901537
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:253
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 253
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 0.9945 0.09743 0.006126 0.004856
## x2 -0.8117 0.07456 0.004688 0.004688
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.8481 0.9335 0.9857 1.0409 1.2098
## x2 -0.9668 -0.8504 -0.8153 -0.7585 -0.6874
##
## ========================================================
## Indirect:
##
## Iterations = 1:253
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 253
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 1.3689 2.667 0.1677 0.1677
## x2 -0.2642 1.986 0.1249 0.1249
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 0.2631 0.5489 0.73410 1.1124 9.3632
## x2 -2.4803 -0.2770 -0.05042 0.1076 0.5526
##
## ========================================================
## Total:
##
## Iterations = 1:253
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 253
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## x1 2.363 2.721 0.1710 0.1710
## x2 -1.076 2.015 0.1267 0.1267
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## x1 1.159 1.514 1.7255 2.1040 10.5083
## x2 -3.426 -1.090 -0.8699 -0.7092 -0.1838
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## x1 0.09743357 2.666743 2.720518
## x2 0.07456303 1.986151 2.015460
##
## Simulated z-values:
## Direct Indirect Total
## x1 10.20695 0.5133141 0.8687228
## x2 -10.88563 -0.1330090 -0.5337944
##
## Simulated p-values:
## Direct Indirect Total
## x1 < 2.22e-16 0.60773 0.38500
## x2 < 2.22e-16 0.89419 0.59348summary(impacts(gns_col, listw = lw_c, R = 300), zstats = TRUE)## Impact measures (sacmixed, exact):
## Direct Indirect Total
## INC -0.9838744 0.6504729 -0.3334015
## HOVAL -0.2779118 -0.3446018 -0.6225136
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:170
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 170
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC -1.0332 0.5625 0.04314 0.04314
## HOVAL -0.2753 0.1595 0.01223 0.01223
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.9037 -1.2999 -1.0117 -0.7742 -0.18730
## HOVAL -0.5208 -0.3433 -0.2869 -0.2088 -0.02505
##
## ========================================================
## Indirect:
##
## Iterations = 1:170
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 170
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC 1.0998 10.39 0.7971 0.7971
## HOVAL -0.2927 3.25 0.2493 0.2493
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -1.836 0.02435 0.4946 1.05626 8.6218
## HOVAL -3.107 -0.44547 -0.2301 -0.07977 0.8262
##
## ========================================================
## Total:
##
## Iterations = 1:170
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 170
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## INC 0.06665 10.81 0.8287 0.8287
## HOVAL -0.56801 3.37 0.2584 0.2584
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## INC -3.225 -0.9671 -0.6535 0.04396 7.8674
## HOVAL -3.455 -0.7109 -0.4742 -0.37465 0.6704
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## INC 0.5625249 10.393146 10.805067
## HOVAL 0.1595073 3.250207 3.369674
##
## Simulated z-values:
## Direct Indirect Total
## INC -1.836649 0.10582045 0.006168089
## HOVAL -1.725789 -0.09006702 -0.168565956
##
## Simulated p-values:
## Direct Indirect Total
## INC 0.066262 0.91572 0.99508
## HOVAL 0.084385 0.92823 0.866143.23 Diagnostik & Pemilihan Model
Langkah praktis:
- OLS awal → cek Moran’s I residual.
- LM tests (lag vs error) → indikasi pilih SAR atau SEM; robust LM jika dua-duanya signifikan.
- Pertimbangan teori: adakah mekanisme endogenous interaction? adakah variabel terabaikan?
- Bandingkan dengan SDM/SDEM/SAC/GNS bila bukti spillover kovariat atau error kuat.
- Interpretasi menggunakan impacts untuk model dengan lag \(y\).
# Grid (simulasi)
ols_grid <- lm(y ~ ., data = dat_grid)
moran_grid <- lm.morantest(ols_grid, listw = lw)
lmtests_grid <- lm.LMtests(ols_grid, listw = lw, test = c("LMlag","LMerr","RLMlag","RLMerr"))
moran_grid; lmtests_grid##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = y ~ ., data = dat_grid)
## weights: lw
##
## Moran I statistic standard deviate = 5.5575, p-value = 1.368e-08
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.190959210 -0.004498787 0.001236931##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ ., data = dat_grid)
## test weights: listw
##
## RSlag = 33.202, df = 1, p-value = 8.308e-09
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ ., data = dat_grid)
## test weights: listw
##
## RSerr = 28.591, df = 1, p-value = 8.94e-08
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ ., data = dat_grid)
## test weights: listw
##
## adjRSlag = 6.3766, df = 1, p-value = 0.01156
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = y ~ ., data = dat_grid)
## test weights: listw
##
## adjRSerr = 1.7658, df = 1, p-value = 0.1839# Columbus (nyata)
ols_col <- lm(CRIME ~ INC + HOVAL, data = d_col)
moran_col <- lm.morantest(ols_col, listw = lw_c)
lmtests_col <- lm.LMtests(ols_col, listw = lw_c, test = c("LMlag","LMerr","RLMlag","RLMerr"))
moran_col; lmtests_col##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = d_col)
## weights: lw_c
##
## Moran I statistic standard deviate = 4.7784, p-value = 8.836e-07
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.330982441 -0.033095098 0.005805358##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = d_col)
## test weights: listw
##
## RSlag = 15.933, df = 1, p-value = 6.561e-05
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = d_col)
## test weights: listw
##
## RSerr = 15.186, df = 1, p-value = 9.74e-05
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = d_col)
## test weights: listw
##
## adjRSlag = 3.8044, df = 1, p-value = 0.05112
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = CRIME ~ INC + HOVAL, data = d_col)
## test weights: listw
##
## adjRSerr = 3.0573, df = 1, p-value = 0.08037Checklist ringkas pemilihan model.
- SAR: Ada alasan kuat bahwa \(y\) wilayah memengaruhi \(y\) tetangga (umpan-balik).
- SEM: Tidak ada mekanisme endogen pada \(y\) tetapi ada error spasial (variabel terabaikan yang spatial).
- SDM: Seperti SAR, namun juga yakin ada spillover dari kovariat (\(WX\)).
- SDEM: Tidak ada lag \(y\), tapi ada \(WX\) dan error spasial.
- SAC/SARAR: SAR + SEM — interaksi endogen dan omitted variables spatial.
- SLX: Hanya ingin mengukur spillover kovariat.
- GNS: Model umum → reduksi berbasis uji/teori agar hemat parameter & interpretasi.
3.24 Keterbatasan & Tips
- Identifiability: Pada SDM/GNS, kolinearitas antara \(X\) dan \(WX\) bisa tinggi; gunakan regularisasi atau kurangi fitur.
- Boundary: Pastikan \(\rho\) dan \(\lambda\) berada pada domain stabilitas (spektrum \(W\)).
- Interpretasi: Gunakan
impacts()untuk lag-\(y\); jangan mengartikan koefisien \(\beta\) langsung sebagai efek marjinal jangka panjang. - W: Pilihan matriks bobot (queen/rook, KNN, jarak) sangat memengaruhi hasil; lakukan sensitivitas.
4 Geographically Weighted Regression (GWR)
Geographically Weighted Regression (GWR) adalah metode regresi yang memperhitungkan variasi spasial dalam data. GWR memungkinkan parameter model regresi berubah di setiap lokasi, memberikan fleksibilitas untuk menangkap pola lokal daripada satu model global yang seragam.
Kapan Menggunakan GWR?
GWR cocok digunakan ketika terdapat indikasi bahwa hubungan antara variabel dependen dan independen berubah-ubah berdasarkan lokasi geografis. GWR sering digunakan dalam penelitian yang melibatkan data spasial seperti ekologi, geografi, epidemiologi, dan ekonomi regional.
4.1 Model GWR
Model GWR dituliskan sebagai:
\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) x_{ik} + \varepsilon_i \]
Di mana:
\(y_i\) adalah nilai variabel dependen di lokasi \(i\).
\(\beta_k(u_i, v_i)\) adalah parameter regresi untuk variabel \(k\) pada lokasi geografis \((u_i, v_i)\).
\(x_{ik}\) adalah variabel independen ke-\(k\) pada lokasi \(i\).
\(\varepsilon_i\) adalah error residual pada lokasi \(i\).
Setiap lokasi \((u_i, v_i)\) memiliki nilai \(\beta_k\) yang berbeda, membuat GWR model lokal di setiap titik.
4.2 Pemilihan Bandwidth dalam GWR
Bandwidth dalam GWR menentukan seberapa banyak pengaruh lokasi tetangga di sekitar titik yang diamati. Pilihan bandwidth sangat mempengaruhi hasil estimasi parameter. Terdapat dua jenis bandwidth yang umum digunakan:
Fixed Bandwidth: Semua titik observasi menggunakan ukuran bandwidth yang sama, baik untuk lokasi yang padat maupun jarang.
Adaptive Bandwidth: Bandwidth disesuaikan dengan kepadatan titik observasi, sehingga di wilayah dengan observasi lebih sedikit, bandwidth diperluas untuk memasukkan lebih banyak tetangga.
4.2.1 3.1 Metode Pemilihan Bandwidth
Pemilihan bandwidth yang optimal dilakukan dengan metode berikut:
Cross-Validation (CV): Bandwidth dipilih yang meminimalkan nilai fungsi loss pada data yang tidak dilatih.
Akaike Information Criterion (AICc): Bandwidth dipilih berdasarkan nilai AIC yang disesuaikan, dengan memilih bandwidth yang meminimalkan AIC.
4.3 Metode Estimasi dalam GWR
Untuk mengestimasi parameter dalam model GWR, digunakan metode weighted least squares (WLS). Setiap titik observasi memiliki bobot yang berbeda berdasarkan jaraknya dari titik target. Secara matematis, estimasi \(\beta_k(u_i, v_i)\) diestimasi dengan:
\[ \hat{\beta}(u_i, v_i) = (X^T W_i X)^{-1} X^T W_i y \]
Di mana:
\(W_i\) adalah matriks bobot untuk lokasi ke-\(i\), dan elemen dari matriks ini ditentukan oleh fungsi kernel yang bergantung pada jarak antara titik \(i\) dan titik lainnya.
\(X\) adalah matriks desain (covariates), dan \(y\) adalah vektor variabel dependen.
4.4 Fungsi Kernel
Bobot dihitung menggunakan fungsi kernel yang mendekati nol ketika jarak semakin besar. Beberapa jenis kernel yang umum digunakan:
- Gaussian Kernel:
\[ w_{ij} = \exp \left( - \frac{d_{ij}^2}{2h^2} \right) \]
- Bisquare Kernel:
\[ w_{ij} = \left( 1 - \left( \frac{d_{ij}}{h} \right)^2 \right)^2 \text{ jika } d_{ij} < h, \text{ sebaliknya } w_{ij} = 0 \]
4.5 Pengujian dalam GWR
4.5.1 Pengujian Signifikansi Global
Untuk menguji apakah GWR lebih baik daripada model regresi global, dapat digunakan test statistik F. Pengujian ini membandingkan varians dari model global dengan varians yang dihasilkan oleh GWR.
4.5.2 Pengujian Signifikansi Lokal
Di setiap lokasi, koefisien \(\beta_k(u_i, v_i)\) dapat diuji secara lokal. t-statistik dapat dihitung untuk setiap koefisien pada setiap lokasi, guna melihat apakah parameter tersebut signifikan secara statistik pada level tertentu.
4.6 Teori Uji Asumsi dalam GWR
Sama seperti regresi linear biasa, GWR memiliki asumsi-asumsi yang harus diuji. Namun, beberapa asumsi ini memiliki nuansa berbeda karena adanya faktor spasial.
4.6.1 Uji Normalitas Residual
Normalitas residual dapat diuji dengan metode Shapiro-Wilk Test. Hipotesis nol dari uji ini adalah bahwa residual berdistribusi normal.
Rumus uji Shapiro-Wilk:
\[ W = \frac{\left( \sum_{i=1}^{n} a_i r_{(i)} \right)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]
Di mana:
\(r_{(i)}\) adalah residual yang diurutkan.
\(a_i\) adalah koefisien yang diperoleh dari matriks kovarian teoretis distribusi normal.
\(\bar{y}\) adalah rata-rata variabel respon.
Jika nilai \(p\) kecil (misalnya, kurang dari 0.05), maka dapat disimpulkan bahwa residual tidak berdistribusi normal.
4.6.2 Uji Homoskedastisitas
Untuk memeriksa homoskedastisitas (atau varians yang tetap), Breusch-Pagan Test dapat digunakan. Uji ini memeriksa apakah error residual bervariasi secara sistematis terhadap variabel prediktor.
Rumus uji Breusch-Pagan:
\[ BP = \frac{n}{2} R^2 \]
Di mana:
\(R^2\) adalah koefisien determinasi dari model regresi tambahan yang mengkuadratkan residual sebagai variabel dependen.
\(n\) adalah jumlah observasi.
Hipotesis nol menyatakan bahwa varians residual konstan (homoskedastisitas). Jika nilai \(p\) kecil, ada indikasi heteroskedastisitas.
4.6.3 Uji Autokorelasi Spasial
Uji autokorelasi spasial, seperti Moran’s I, digunakan untuk memeriksa apakah residual memiliki pola spasial yang signifikan.
Rumus Moran’s I:
\[ I = \frac{n}{W} \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij} (y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]
Di mana:
\(n\) adalah jumlah observasi.
\(w_{ij}\) adalah elemen dari matriks bobot spasial.
\(\bar{y}\) adalah rata-rata nilai residual.
Nilai \(I\) yang mendekati +1 menunjukkan autokorelasi positif, sementara nilai mendekati -1 menunjukkan autokorelasi negatif. Jika \(I\) signifikan, berarti ada pola spasial dalam residual.
4.6.4 Uji Multikolinearitas
Untuk mendeteksi multikolinearitas antara variabel independen, dapat digunakan Variance Inflation Factor (VIF).
Rumus VIF:
\[ VIF_k = \frac{1}{1 - R_k^2} \]
Di mana:
- \(R_k^2\) adalah koefisien determinasi dari regresi variabel \(X_k\) terhadap variabel independen lainnya.
Jika VIF lebih dari 10, hal ini menunjukkan adanya multikolinearitas yang serius, sehingga model perlu dikoreksi.
4.7 Kesimpulan
GWR adalah alat yang kuat untuk analisis data spasial karena mampu menangkap variasi lokal dalam hubungan antara variabel. Dengan pemilihan bandwidth yang tepat, estimasi yang akurat, dan pengujian asumsi yang lengkap, GWR dapat memberikan wawasan yang lebih dalam tentang bagaimana hubungan antara variabel berubah di ruang geografis.
4.8 APLIKASI - R
4.8.1 Derajat (Latitude longitude)
data_spasial <- data.frame(
y = runif(100),
x1 = runif(100),
x2 = runif(100),
longitude = runif(100, min = -100, max = -90),
latitude = runif(100, min = 30, max = 40)
)
coordinates(data_spasial) <- ~longitude+latitude4.8.2 Optimal Bandwidth
bandwidth_optimal <- gwr.sel(
y ~ x1 + x2,
data = data_spasial,
coords = cbind(data_spasial$longitude, data_spasial$latitude),
adapt = TRUE
)## Adaptive q: 0.381966 CV score: 8.118229
## Adaptive q: 0.618034 CV score: 8.128019
## Adaptive q: 0.236068 CV score: 8.143954
## Adaptive q: 0.4636324 CV score: 8.1206
## Adaptive q: 0.3262379 CV score: 8.126476
## Adaptive q: 0.4115305 CV score: 8.117914
## Adaptive q: 0.4037411 CV score: 8.117152
## Adaptive q: 0.3978179 CV score: 8.116531
## Adaptive q: 0.391763 CV score: 8.117982
## Adaptive q: 0.3955052 CV score: 8.117039
## Adaptive q: 0.3994483 CV score: 8.116204
## Adaptive q: 0.401088 CV score: 8.116415
## Adaptive q: 0.3996285 CV score: 8.116169
## Adaptive q: 0.3999755 CV score: 8.116103
## Adaptive q: 0.4004004 CV score: 8.116216
## Adaptive q: 0.4001378 CV score: 8.116139
## Adaptive q: 0.3999196 CV score: 8.116114
## Adaptive q: 0.4000162 CV score: 8.116103
## Adaptive q: 0.3999755 CV score: 8.1161034.8.3 Running GWR
model_gwr <- gwr(
y ~ x1 + x2,
data = data_spasial,
coords = cbind(data_spasial$longitude, data_spasial$latitude),
adapt = bandwidth_optimal,
hatmatrix = TRUE,
se.fit = TRUE
)4.8.4 Koefisien pada setiap lokasi
coef_gwr <- model_gwr$SDF@data[, c("x1", "x2")]4.8.5 summary
# Menampilkan ringkasan model
summary(model_gwr)## Length Class Mode
## SDF 100 SpatialPointsDataFrame S4
## lhat 10000 -none- numeric
## lm 11 -none- list
## results 14 -none- list
## bandwidth 100 -none- numeric
## adapt 1 -none- numeric
## hatmatrix 1 -none- logical
## gweight 1 -none- character
## gTSS 1 -none- numeric
## this.call 7 -none- call
## fp.given 1 -none- logical
## timings 12 -none- numeric4.8.6 UTM (Universal Transverse Mercator)
4.8.6.1 Mendefinisikan CRS latitude/longitude (WGS84)
proj_longlat <- CRS("+proj=longlat +datum=WGS84")4.8.6.2 Mendefinisikan CRS UTM, misalnya UTM zona 33N
proj_utm <- CRS("+proj=utm +zone=33 +datum=WGS84")4.8.6.3 Definisikan CRS untuk data longitude/latitude (WGS84)
proj4string(data_spasial) <- CRS("+proj=longlat +datum=WGS84")4.8.6.4 Konversi data spatial ke UTM
data_spasial_utm <- spTransform(data_spasial, proj_utm)4.8.6.5 Menjalankan GWR dengan data yang sudah terproyeksi ke UTM
bandwidth_optimal <- gwr.sel(y ~ x1 + x2, data = data_spasial_utm, adapt = TRUE)## Adaptive q: 0.381966 CV score: 8.131951
## Adaptive q: 0.618034 CV score: 8.129777
## Adaptive q: 0.763932 CV score: 8.124636
## Adaptive q: 0.854102 CV score: 8.112163
## Adaptive q: 0.9098301 CV score: 8.107845
## Adaptive q: 0.9748563 CV score: 8.101248
## Adaptive q: 0.9500185 CV score: 8.105789
## Adaptive q: 0.9844604 CV score: 8.100525
## Adaptive q: 0.9916947 CV score: 8.100276
## Adaptive q: 0.9951674 CV score: 8.100056
## Adaptive q: 0.9970133 CV score: 8.099937
## Adaptive q: 0.9981541 CV score: 8.099863
## Adaptive q: 0.9988592 CV score: 8.099817
## Adaptive q: 0.9992949 CV score: 8.099788
## Adaptive q: 0.9995642 CV score: 8.099771
## Adaptive q: 0.9997307 CV score: 8.09976
## Adaptive q: 0.9998336 CV score: 8.099753
## Adaptive q: 0.9998971 CV score: 8.099749
## Adaptive q: 0.9999378 CV score: 8.099746
## Adaptive q: 0.9999378 CV score: 8.099746model_gwr <- gwr(y ~ x1 + x2, data = data_spasial_utm, adapt = bandwidth_optimal)4.9 UJI ASUMSI
4.9.1 Heterogeneity
4.9.1.1 Model OLS untuk perbandingan
model_ols <- lm(y ~ x1 + x2, data = data_spasial)
summary(model_ols)##
## Call:
## lm(formula = y ~ x1 + x2, data = data_spasial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50760 -0.24753 0.04085 0.22413 0.60271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.39876 0.06920 5.762 9.8e-08 ***
## x1 -0.03077 0.08836 -0.348 0.7284
## x2 0.25947 0.09943 2.610 0.0105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2799 on 97 degrees of freedom
## Multiple R-squared: 0.06595, Adjusted R-squared: 0.0467
## F-statistic: 3.425 on 2 and 97 DF, p-value: 0.03655bptest(model_ols)##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 4.7724, df = 2, p-value = 0.091984.9.1.2 Lakukan perbandingan AIC (Semakin rendah AIC, semakin baik model)
4.9.1.3 Jika AIC model GWR lebih rendah dari model OLS, maka ada heterogenitas spasial.
AIC(model_ols)## [1] 34.10179model_gwr$results$AICc## NULL4.9.2 Multikolinearity
4.9.2.1 Gunakan VIF untuk mendeteksi multikolinearitas
vif(model_ols)## x1 x2
## 1.00338 1.003384.9.3 Uji Spatial Autokorelasi
4.9.3.1 Ekstraksi residual GWR
residual_gwr <- model_gwr$SDF$gwr.e4.9.3.2 Membuat matriks bobot spasial berdasarkan tetangga terdekat
coords <- cbind(data_spasial$longitude, data_spasial$latitude)
nb <- knn2nb(knearneigh(coords, k=4))
listw <- nb2listw(nb, style="W")4.9.3.3 Uji Moran’s I untuk residual spasial
moran.test(residual_gwr, listw)##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: listw
##
## Moran I statistic standard deviate = -0.76665, p-value = 0.7784
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.058962898 -0.010101010 0.0040620434.9.4 Uji Normality
shapiro.test(residual_gwr)##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.95395, p-value = 0.0015244.9.5 Membuat plot -spplot
4.9.6 Ekstrak koefisien lokal untuk variabel x1 dan x2
coef_x1 <- model_gwr$SDF$x1
coef_x2 <- model_gwr$SDF$x24.9.7 Plot koefisien regresi lokal untuk variabel x1
spplot(model_gwr$SDF, "x1", main = "Koefisien GWR untuk x1", col.regions = terrain.colors(100))
4.9.8 Plot koefisien regresi lokal untuk variabel x2
spplot(model_gwr$SDF, "x2", main = "Koefisien GWR untuk x2", col.regions = terrain.colors(100))
4.9.9 Membuat plot -gglot2
4.9.10 Konversi SpatialPointsDataFrame menjadi data frame biasa
gwr_df <- as.data.frame(model_gwr$SDF)
# Plot koefisien lokal untuk x1
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
geom_point(aes(color = x1), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x1", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x2
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
geom_point(aes(color = x2), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x2", color = "Koefisien") +
theme_minimal()
4.10 Membuat Surface Plot
4.11 IDW
4.11.1 Assuming data_spasial is in WGS84
proj4string(data_spasial)## [1] "+proj=longlat +datum=WGS84 +no_defs"data_spasial_wgs84 <- spTransform(data_spasial, CRS("+proj=longlat +datum=WGS84 +no_defs"))4.11.2 Create the grid with the same CRS
x.range <- range(data_spasial$longitude)
y.range <- range(data_spasial$latitude)4.11.3 Buat grid (titik-titik prediksi di mana kita ingin mengestimasi nilai koefisien)
grd <- expand.grid(
longitude = seq(from = x.range[1], to = x.range[2], length.out = 100),
latitude = seq(from = y.range[1], to = y.range[2], length.out = 100)
)
coordinates(grd) <- ~longitude + latitude
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84 +no_defs")4.11.4 Gabungkan koefisien dengan data spasial
data_spasial$coef_x1 <- coef_x14.11.5 Lakukan IDW untuk variabel koefisien x1
idw_output <- idw(coef_x1 ~ 1, data_spasial, newdata = grd, idp = 2.0) # idp = power parameter## [inverse distance weighted interpolation]4.12 Konversi output IDW menjadi data frame untuk ggplot
idw_df <- as.data.frame(idw_output)ggplot() +
geom_tile(data = idw_df, aes(x = longitude, y = latitude, fill = var1.pred)) +
scale_fill_viridis_c() +
labs(title = "Surface Plot of GWR Coefficients for x1 (IDW)",
fill = "Koefisien x1") +
theme_minimal()
5 Spatial Interpolation
5.1 Apa itu Spatial Interpolation?
Spatial interpolation adalah prosedur memperkirakan nilai variabel pada lokasi tidak teramati berdasarkan nilai di lokasi teramati di sekitarnya. Tujuannya membuat permukaan kontinu (surface) dari data titik, misalnya peta curah hujan, kualitas udara, kandungan logam tanah, atau ketinggian.
Intuisi: lokasi yang berdekatan cenderung memiliki nilai yang mirip (Tobler’s first law). Perbedaannya ada pada cara memberi bobot dan asumsi model yang digunakan untuk “menjahit” titik-titik menjadi peta kontinu.
5.2 Asumsi dan Pertimbangan Umum dalam Spatial Interpolation
Interpolasi spasial bukan sekadar masalah matematis — ia juga
memerlukan pemahaman mendalam tentang struktur spasial
data, asumsi model, dan implikasi
praktis dari setiap pendekatan.
Bagian ini menguraikan secara rinci asumsi-asumsi utama yang mendasari
sebagian besar metode interpolasi.
5.2.1 Autokorelasi Spasial
Autokorelasi spasial merupakan fondasi utama dari
semua metode interpolasi. Prinsip dasarnya dikenal sebagai Hukum
Pertama Geografi (Tobler, 1970):
> “Everything is related to everything else, but near things are more
related than distant things.”
5.2.1.1 Penjelasan Konseptual
Autokorelasi spasial menyatakan bahwa nilai suatu variabel di
lokasi tertentu cenderung mirip dengan nilai di lokasi lain yang
berdekatan secara geografis.
Jika \(Z(s)\) adalah nilai variabel di
lokasi \(s\), maka untuk dua lokasi
berdekatan \(s_i\) dan \(s_j\):
\[ \text{Cov}(Z(s_i), Z(s_j)) = C(h) \quad \text{dengan } h = \|s_i - s_j\| \]
Artinya, semakin kecil jarak \(h\), semakin besar korelasi antar nilai.
5.2.1.2 Implikasi Praktis
- Metode seperti IDW, Kriging, dan TPS mengandalkan
ide ini untuk menimbang pengaruh data sekitar.
- Jika data menunjukkan autokorelasi negatif (nilai yang berdekatan justru sangat berbeda), interpolasi konvensional menjadi tidak valid tanpa modifikasi model.
5.2.2 Stasioneritas dan Tren
5.2.2.1 Stasioneritas
Metode geostatistik klasik, terutama Kriging, mengasumsikan bahwa data berasal dari proses spasial stasioner, yakni:
- Mean konstan: \(E[Z(s)] =
\mu\), sama di seluruh wilayah.
- Kovarians hanya tergantung pada jarak: \(\text{Cov}(Z(s_i), Z(s_j)) = C(h)\), bukan pada lokasi absolut.
Namun dalam kenyataannya, banyak fenomena memiliki tren global (misal: ketinggian yang meningkat dari barat ke timur). Dalam kasus ini, asumsi stasioneritas ketat tidak terpenuhi.
5.2.2.2 Tren Global
Jika terdapat pola sistematik, misalnya perubahan rata-rata secara
spasial, pendekatan seperti Universal Kriging (UK) atau
Trend Surface Analysis lebih sesuai.
Keduanya mengizinkan mean yang bervariasi dengan lokasi:
\[ E[Z(s)] = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon(s) \]
Di mana \(\varepsilon(s)\) tetap dianggap stasioner secara lokal.
5.2.2.3 Implikasi Praktis
- Kriging biasa (OK) cocok untuk data dengan mean
konstan.
- Kriging universal (UK) cocok untuk data dengan
gradien atau tren global.
- Penting untuk memeriksa tren sebelum fitting variogram agar tidak salah tafsir korelasi spasial.
5.2.3 Skala dan Anisotropi
5.2.3.1 Skala (Scale Dependency)
Hubungan spasial bisa berbeda pada berbagai skala observasi. Misalnya, variabilitas suhu pada skala 1 km bisa berbeda drastis dengan skala 100 km. Oleh karena itu, pemilihan resolusi grid dan parameter variogram harus sesuai dengan skala analisis.
\[ C(h) = \sigma^2 \exp\left(-\frac{h}{a}\right) \]
Parameter \(a\) (range) mengatur seberapa jauh pengaruh spasial berlaku. Jika skala salah dipilih, interpolasi bisa terlalu halus atau terlalu kasar.
5.2.3.2 Anisotropi
Anisotropi terjadi ketika korelasi spasial bergantung pada arah. Misalnya, penyebaran polutan bisa lebih cepat ke arah angin dibanding arah tegak lurusnya.
Secara matematis, anisotropi dapat dimodelkan dengan transformasi jarak arah-spesifik:
\[ h' = \sqrt{ \left(\frac{x}{a_x}\right)^2 + \left(\frac{y}{a_y}\right)^2 } \]
di mana \(a_x\) dan \(a_y\) menunjukkan jangkauan korelasi yang berbeda menurut arah.
5.2.3.3 Implikasi Praktis
- Cek anisotropi sebelum fitting variogram (
vgm()di R bisa mendeteksi lewatalpha).
- Gunakan directional variogram untuk menganalisis
apakah pola berbeda menurut arah.
- Untuk fenomena seperti curah hujan atau penyebaran polusi udara, anisotropi seringkali signifikan.
5.2.4 Kerapatan dan Sebaran Sampel
Interpolasi membutuhkan titik pengamatan yang cukup rapat dan merata agar hasil representatif.
5.2.4.1 Masalah Umum
- Sparse sampling: terlalu sedikit titik menyebabkan
ketidakpastian besar di antara area kosong.
- Clustered sampling: terlalu banyak titik pada area kecil membuat model bias ke wilayah padat.
5.2.4.2 Prinsip Heuristik
- Setidaknya, setiap area berukuran “range” (dari
variogram) harus memiliki minimal satu titik data.
- Gunakan peta densitas titik untuk menilai distribusi observasi.
5.2.4.3 Implikasi Praktis
Distribusi titik menentukan pilihan metode:
- Data jarang → gunakan metode berbasis model (Kriging).
- Data rapat → metode deterministik (IDW, TPS) sering cukup baik.
5.2.5 Validasi Model
Tidak ada model interpolasi yang sempurna; karena itu validasi sangat penting.
5.2.5.1 Cross-validation
Gunakan Leave-One-Out (LOO) atau k-fold cross-validation untuk menguji seberapa baik model memprediksi data yang tidak digunakan dalam fitting.
\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\hat{Z}(s_i) - Z(s_i))^2} \]
Metrik lain yang umum:
- MAE (Mean Absolute Error): mengukur rata-rata
deviasi absolut.
- ME (Mean Error/Bias): menunjukkan kecenderungan over/under prediction.
5.2.5.2 Interpretasi Metrik
| Metrik | Interpretasi | Nilai Ideal |
|---|---|---|
| RMSE | Akurasi keseluruhan | Sekecil mungkin |
| MAE | Stabilitas prediksi | Sekecil mungkin |
| ME | Bias rata-rata | Mendekati 0 |
5.2.5.3 Validasi Visual
- Scatter plot antara nilai observasi dan
prediksi.
- Peta residual untuk melihat bias spasial.
- Variogram residual untuk memastikan tidak ada struktur spasial tersisa.
5.2.6 Kesimpulan
Asumsi-asumsi di atas tidak hanya formalitas, tetapi komponen
kunci keberhasilan interpolasi spasial.
Memahami autokorelasi, tren, anisotropi, dan distribusi data membantu
memastikan bahwa hasil interpolasi tidak hanya halus secara visual,
tetapi juga valid secara ilmiah.
“A map is only as good as the assumptions behind it.” — Adapted from Cressie (1993)
5.3 Metode Interpolasi yang Umum
Berikut kanonik metode yang sering dipakai (dari yang non-parametrik hingga berbasis model stokastik):
- Nearest Neighbor (NN): nilai lokasi terdekat; cepat tetapi kasar (peta bloky).
- Natural Neighbor / TIN: bobot berdasar Voronoi neighborhoods; halus dan konservatif (butuh lib tertentu).
- Inverse Distance Weighting (IDW): bobot \(w_i \propto d(s_0,s_i)^{-p}\) dengan orde \(p>0\); intuitif, tanpa model kovarians.
- Trend Surface / Polynomial Regression: permukaan deterministik berbasis fungsi koordinat (x, y).
- Splines (mis. Thin Plate Spline): permukaan halus
yang meminimalkan bending energy (butuh paket tambahan seperti
fields). - Kriging (SK/OK/UK): penduga linier tak bias dengan variansi minimum berbasis variogram/kovarians; varian: Simple, Ordinary, Universal.
- Co-Kriging: memanfaatkan variabel sekunder yang berkorelasi spasial.
- Radial Basis Functions (RBF): keluarga fungsi dasar radial (multikuadrik, Gaussian) untuk permukaan halus.
Pada modul ini, contoh kode difokuskan pada IDW dan Kriging dengan
gstat, karena keduanya populer dan tanpa dependensi tambahan.
5.4 Nearest Neighbor (NN)
Definisi
Metode Nearest Neighbor (NN) atau Thiessen interpolation adalah teknik paling sederhana dalam interpolasi spasial. Untuk suatu lokasi target \(s_0\), kita cari satu titik pengamatan terdekat \(s_{(1)}\), lalu nilai prediksinya langsung diambil dari titik itu:
\[ \hat{Z}(s_0) = Z\big(s_{(1)}\big). \]
Dengan kata lain, titik yang paling dekat “memenangkan” wilayahnya. Secara geometris, seluruh ruang terbagi menjadi poligon-poligon yang disebut Voronoi cells — setiap sel mencakup semua lokasi yang lebih dekat ke satu titik pengamatan dibandingkan ke titik lain.
Contoh Manual
Misalkan kita punya tiga titik pengamatan suhu udara (°C):
| Lokasi | Koordinat (x, y) | Nilai Z |
|---|---|---|
| A | (0, 0) | 28 |
| B | (2, 0) | 30 |
| C | (0, 3) | 25 |
Kita ingin memperkirakan suhu di titik baru \(s_0 = (1, 1)\).
Langkah Perhitungan
Hitung jarak Euclidean ke tiap titik:
\[ \begin{aligned} d(A, s_0) &= \sqrt{(1-0)^2 + (1-0)^2} = \sqrt{2} = 1.414, \\ d(B, s_0) &= \sqrt{(1-2)^2 + (1-0)^2} = \sqrt{2} = 1.414, \\ d(C, s_0) &= \sqrt{(1-0)^2 + (1-3)^2} = \sqrt{10} = 3.162. \end{aligned} \]
Titik A dan B sama-sama paling dekat.
Jika ada jarak sama (tie), ambil salah satu secara deterministik (misalnya yang pertama dalam urutan data), atau gunakan rata-rata dari kandidat terdekat jika ingin hasil lebih halus. Di sini kita ambil A.
Prediksi nilai di \(s_0\):
\[ \hat{Z}(s_0) = Z(A) = 28. \]
Implementasi Ringkas di R
Kode berikut mengilustrasikan perhitungan jarak, pemilihan tetangga terdekat, dan (opsional) sketsa mosaik Voronoi untuk tiga titik tersebut.
# Data titik pengamatan
pts <- data.frame(
id = c("A", "B", "C"),
x = c(0, 2, 0),
y = c(0, 0, 3),
Z = c(28, 30, 25)
)
# Titik target
s0 <- data.frame(x = 1, y = 1)
# Fungsi jarak Euclidean
edist <- function(x1, y1, x2, y2) sqrt((x1 - x2)^2 + (y1 - y2)^2)
# Hitung jarak ke setiap titik pengamatan
pts$dist_s0 <- with(pts, edist(x, y, s0$x, s0$y))
# Urutkan berdasarkan jarak
pts_ord <- pts[order(pts$dist_s0), ]
pts_ord## id x y Z dist_s0
## 1 A 0 0 28 1.414214
## 2 B 2 0 30 1.414214
## 3 C 0 3 25 2.236068# Ambil tetangga terdekat (tie-breaker: urutan data)
nearest <- pts_ord[1, ]
cat(sprintf("Nearest to s0=(%.1f,%.1f) is %s at distance %.3f with Z=%g\n",
s0$x, s0$y, nearest$id, nearest$dist_s0, nearest$Z))## Nearest to s0=(1.0,1.0) is A at distance 1.414 with Z=28Zhat <- nearest$Z
Zhat## [1] 28# Visualisasi Voronoi (opsional)
# Memerlukan paket deldir dan ggplot2
if (requireNamespace("deldir", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)) {
library(deldir)
library(ggplot2)
# Batas plot
xlim <- c(-0.5, 2.5)
ylim <- c(-0.5, 3.5)
# Delaunay/Voronoi tessellation
dd <- deldir::deldir(pts$x, pts$y, rw = c(xlim, ylim))
tiles <- deldir::tile.list(dd)
# Konversi tile ke data.frame untuk ggplot
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(
id = pts$id[i],
x = tiles[[i]]$x,
y = tiles[[i]]$y
)
}))
ggplot() +
geom_polygon(data = poly_df, aes(x = x, y = y, group = id),
fill = NA, color = "grey40") +
geom_point(data = pts, aes(x = x, y = y), size = 3) +
geom_text(data = pts, aes(x = x, y = y, label = id), vjust = -1) +
geom_point(data = s0, aes(x = x, y = y), shape = 4, size = 3, stroke = 1.2) +
coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
labs(title = "Sketsa Voronoi untuk Nearest Neighbor",
subtitle = "Titik s0 ditandai dengan silang (×)",
x = "x", y = "y") +
theme_minimal(base_size = 12)
} else {
message("Lewati plot: paket 'deldir' dan/atau 'ggplot2' belum terpasang.")
}
Sifat, Kelebihan, dan Kelemahan
- Kelebihan: sangat cepat, tanpa parameter, dan tidak membutuhkan variogram.
- Kelemahan: hasil tidak halus (blocky), tidak cocok untuk fenomena yang berubah gradual, dan sensitif terhadap distribusi titik.
Metode ini sering dipakai untuk visualisasi cepat (misalnya Thiessen polygons pada data curah hujan) sebelum beralih ke metode berbobot seperti IDW atau Kriging.
Contoh2
# Data meuse
data(meuse, package = "sp")
coordinates(meuse) <- ~x+y
meuse <- set_crs(meuse)
meuse$z <- log(meuse$zinc)
data(meuse.grid, package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS(proj4string(meuse))
# Sanity checks
stopifnot(!is.na(proj4string(meuse)), !is.na(proj4string(meuse.grid)))
# Cek ringkas
head(meuse@data[,c("zinc","copper")])## zinc copper
## 1 1022 85
## 2 1141 81
## 3 640 68
## 4 257 81
## 5 269 48
## 6 281 61# 1) Ambil koordinat observasi & grid (sebagai SpatialPoints)
Xobs <- coordinates(meuse)
grid_sp <- as(meuse.grid, "SpatialPoints")
# 2) Cari indeks tetangga terdekat untuk setiap titik grid
knn <- FNN::get.knnx(data = Xobs, query = coordinates(grid_sp), k = 1)
nn_ix <- knn$nn.index[,1]
# 3) Prediksi NN dan bangun SpatialPixelsDataFrame baru
pred_nn <- meuse$z[nn_ix]
stopifnot(length(pred_nn) == nrow(as.data.frame(meuse.grid)))
nn_grid <- SpatialPixelsDataFrame(points = grid_sp,
data = data.frame(z_nn = pred_nn),
proj4string = CRS(proj4string(meuse.grid)))
# 4) Plot aman — kolom 'z_nn' dijamin ada
spplot(nn_grid, "z_nn", main = "Nearest Neighbor (1-NN) — Prediksi log(Zinc)")
5.5 Natural Neighbor / TIN
Metode Natural Neighbor (atau Sibson interpolation) bekerja berdasarkan Voronoi diagram juga, namun dengan gagasan yang lebih halus. Jika Nearest Neighbor (NN) hanya mengambil satu titik terdekat, Natural Neighbor mempertimbangkan beberapa titik yang “terkait” dengan lokasi baru \(s_0\).
Langkah pikirannya:
- Kita punya sekumpulan titik data (mis. A, B, C, D).
- Dari titik-titik itu kita buat Voronoi diagram — setiap titik punya wilayahnya.
- Tambahkan lokasi baru \(s_0\). Diagram Voronoi berubah: sebagian area dari poligon-poligon tetangga “direbut” oleh \(s_0\).
- Luas potongan yang direbut dari masing-masing tetangga menjadi bobot interpolasi.
Secara konseptual:
\[ \hat{Z}(s_0) = \sum_i w_i \, Z(s_i), \quad \text{dengan } w_i \ge 0, \; \sum_i w_i = 1. \]
Bobot \(w_i\) sebanding dengan proporsi luas Voronoi yang “diambil alih” oleh \(s_0\). Karena bobot positif dan menjumlah ke 1, hasil interpolasi bersifat halus dan kontinu.
5.5.1 Hubungan dengan TIN (Triangulated Irregular Network)
Implementasi Sibson interpolation yang eksak cukup kompleks dan jarang tersedia langsung di R. Alternatif praktis adalah TIN linear interpolation: interpolasi linier di dalam Delaunay triangulation (dual dari Voronoi).
Intinya:
- Bentuk triangulasi Delaunay dari titik-titik data.
- Temukan segitiga yang mengandung \(s_0\).
- Hitung bobot barycentric (perbandingan luas sub-segitiga) terhadap tiga titik sudut.
- Nilai interpolasi adalah rata-rata tertimbang linier dengan bobot barycentric.
5.5.2 Contoh Manual — TIN Linear (Barycentric)
Misalkan tiga titik pengamatan membentuk satu segitiga:
| Lokasi | (x, y) | Nilai Z |
|---|---|---|
| A | (0, 0) | 10 |
| B | (4, 0) | 30 |
| C | (0, 4) | 20 |
Kita ingin menghitung nilai di \(s_0 = (1, 1)\).
1) Cek posisi \(s_0\)
\(s_0\) berada di dalam segitiga ABC karena koordinatnya positif dan \(1 + 1 < 4\).
2) Luas segitiga utama (ABC)
Gunakan rumus determinan (setengah dari determinan):
\[ L_{ABC} = \tfrac{1}{2} \big| x_A(y_B - y_C) + x_B(y_C - y_A) + x_C(y_A - y_B) \big|. \]
Substitusi koordinat menghasilkan \(L_{ABC} = 8\).
3) Luas tiga sub-segitiga yang melibatkan \(s_0\)
- \(L_{s_0BC} = \tfrac{1}{2} \big| 1(0-4) + 4(4-1) + 0(1-0) \big| = 4\)
- \(L_{As_0C} = \tfrac{1}{2} \big| 0(1-4) + 1(4-0) + 0(0-1) \big| = 2\)
- \(L_{ABs_0} = \tfrac{1}{2} \big| 0(0-1) + 4(1-0) + 1(0-0) \big| = 2\)
Total luas sub-segitiga = \(4 + 2 + 2 = 8\), konsisten dengan \(L_{ABC}\).
4) Bobot barycentric
\[ w_A = \frac{L_{s_0BC}}{L_{ABC}} = \frac{4}{8} = 0.5, \quad w_B = \frac{L_{As_0C}}{L_{ABC}} = \frac{2}{8} = 0.25, \quad w_C = \frac{L_{ABs_0}}{L_{ABC}} = \frac{2}{8} = 0.25. \]
5) Interpolasi nilai
\[ \hat{Z}(s_0) = 0.5\cdot 10 + 0.25\cdot 30 + 0.25\cdot 20 = 17.5. \]
Implementasi Ringkas di R
Contoh fungsi kecil untuk menghitung bobot barycentric dan nilai interpolasi pada \(s_0\):
tri_barycentric <- function(A, B, C, s0, Z) {
area <- function(P, Q, R) 0.5 * abs(P[1]*(Q[2]-R[2]) + Q[1]*(R[2]-P[2]) + R[1]*(P[2]-Q[2]))
L_ABC <- area(A, B, C)
L_s0BC <- area(s0, B, C)
L_As0C <- area(A, s0, C)
L_ABs0 <- area(A, B, s0)
wA <- L_s0BC / L_ABC
wB <- L_As0C / L_ABC
wC <- L_ABs0 / L_ABC
Zhat <- wA*Z[1] + wB*Z[2] + wC*Z[3]
list(weights = c(A = wA, B = wB, C = wC), Zhat = Zhat,
areas = c(ABC = L_ABC, s0BC = L_s0BC, As0C = L_As0C, ABs0 = L_ABs0))
}
A <- c(0, 0); B <- c(4, 0); C <- c(0, 4)
s0 <- c(1, 1)
Z <- c(10, 30, 20) # Z_A, Z_B, Z_C
res <- tri_barycentric(A, B, C, s0, Z)
res## $weights
## A B C
## 0.50 0.25 0.25
##
## $Zhat
## [1] 17.5
##
## $areas
## ABC s0BC As0C ABs0
## 8 4 2 2Visualisasi Voronoi Diagram
Visualisasi berikut memperlihatkan sel Voronoi untuk titik A, B, dan C, dengan \(s_0 = (1,1)\) ditandai silang (×).
if (requireNamespace("deldir", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)) {
library(deldir)
library(ggplot2)
pts <- data.frame(id = c("A","B","C"), x = c(0,4,0), y = c(0,0,4))
s0 <- data.frame(x = 1, y = 1)
dd <- deldir::deldir(pts$x, pts$y, rw = c(-0.5,4.5,-0.5,4.5))
tiles <- deldir::tile.list(dd)
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(id = pts$id[i], x = tiles[[i]]$x, y = tiles[[i]]$y)
}))
ggplot() +
geom_polygon(data = poly_df, aes(x=x, y=y, group=id), fill=NA, color="grey40") +
geom_point(data = pts, aes(x=x, y=y), size=3) +
geom_text(data = pts, aes(x=x, y=y, label=id), vjust=-1) +
geom_point(data = s0, aes(x=x, y=y), shape=4, size=3, stroke=1.2) +
coord_equal(xlim=c(-0.5,4.5), ylim=c(-0.5,4.5)) +
labs(title="Voronoi Diagram untuk A(0,0), B(4,0), C(0,4)",
subtitle="Titik s0=(1,1) ditandai dengan silang (×)",
x="x", y="y") +
theme_minimal(base_size = 12)
}
Catatan: Untuk implementasi TIN skala lebih besar beserta triangulasi Delaunay, Anda dapat menggunakan paket seperti
RTriangle,geometry(fungsidelaunayn), ataudeldir(untuk Voronoi/Delaunay 2D), lalu menerapkan interpolasi linier di dalam setiap segitiga.
Contoh 2
Xobs <- coordinates(meuse)
xg <- sort(unique(meuse.grid@coords[,1]))
yg <- sort(unique(meuse.grid@coords[,2]))
tin_out <- interp::interp(x = Xobs[,1], y = Xobs[,2], z = meuse$z,
xo = xg, yo = yg, duplicate = "mean", linear = TRUE)
zmat <- tin_out$z
grd <- expand.grid(x = tin_out$x, y = tin_out$y)
tin_sgdf <- SpatialPixelsDataFrame(points = SpatialPoints(grd),
data = data.frame(z_tin = as.vector(zmat)),
tolerance = 0.5)
proj4string(tin_sgdf) <- CRS(proj4string(meuse))
spplot(tin_sgdf, "z_tin", main = "TIN (Linear) — Perkiraan log(Zinc)")
xg <- sort(unique(meuse.grid@coords[,1]))
yg <- sort(unique(meuse.grid@coords[,2]))
tin_out <- interp::interp(x = Xobs[,1], y = Xobs[,2], z = meuse$z,
xo = xg, yo = yg, duplicate = "mean", linear = TRUE)
zmat <- tin_out$z
grd <- expand.grid(x = tin_out$x, y = tin_out$y)
tin_sgdf <- SpatialPixelsDataFrame(points = SpatialPoints(grd),
data = data.frame(z_tin = as.vector(zmat)),
tolerance = 0.5)
proj4string(tin_sgdf) <- CRS(proj4string(meuse))
spplot(tin_sgdf["z_tin"], main = "TIN (Linear Barycentric) — Perkiraan log(Zinc)")
5.6 Inverse Distance Weighting (IDW)
Inverse Distance Weighting (IDW) adalah metode interpolasi deterministik yang memperkirakan nilai di lokasi target \(s_0\) sebagai rata-rata tertimbang dari nilai-nilai di titik pengamatan \(s_i\), dengan bobot menurun terhadap jarak.
\[ w_i(s_0) \propto d(s_0, s_i)^{-p}, \qquad p > 0, \]
\[ \hat{Z}(s_0) = \frac{\sum_i w_i(s_0)\, Z(s_i)}{\sum_i w_i(s_0)}. \]
Umumnya kita gunakan jarak Euclidean \(d(\cdot,\cdot)\). Parameter power \(p\) mengontrol tingkat pengaruh jarak: makin besar \(p\), makin tajam penurunan bobot.
Sifat & Intuisi Parameter \(p\)
- \(p \to 0\): semua bobot hampir sama (mendekati rata-rata biasa).
- \(p=1\): penurunan bobot linear terhadap jarak (dalam skala log–log).
- \(p=2\): penurunan kuadratik; pilihan populer untuk fenomena yang meredup cepat.
- \(p \to \infty\): IDW mendekati Nearest Neighbor (hanya tetangga terdekat yang dominan).
IDW tidak memerlukan variogram dan tidak memodelkan struktur kovariansi. Ia lokal, cepat, dan deterministik; namun tidak optimal secara statistik (tidak BLUE seperti kriging ketika asumsi terpenuhi).
Kebijakan Tetangga (Neighborhood)
- Semua titik (global) – lebih halus, namun bisa bias bila distribusi titik tidak merata.
- \(k\) tetangga terdekat – stabil di tepi; pilih \(k\) (mis. 6–20) sesuai densitas sampel.
- Radius pencarian – hanya titik dalam radius \(r\) yang digunakan.
Catatan tepi (edge cases): - Jika ada titik tepat di \(s_0\) (\(d=0\)), definisikan \(\hat{Z}(s_0) = Z(s_i)\) untuk menghindari pembagian nol. - Titik duplikat/kolokasi sebaiknya diagregasi (mis. rata-rata) atau ditangani sebagai pengamatan terpisah dengan jarak \(>0\) sangat kecil. - Anisotropi bisa dicapai dengan metrik jarak berbobot/ter-skala (mis. skala sumbu x,y berbeda) atau rotasi sumbu.
Contoh Manual (Data Segitiga)
Tiga titik pengamatan:
| Lokasi | (x, y) | Z |
|---|---|---|
| A | (0, 0) | 10 |
| B | (4, 0) | 30 |
| C | (0, 4) | 20 |
Target: \(s_0 = (1, 1)\). Jarak Euclidean:
\[ d(A,s_0)=\sqrt{2}=1.4142,\quad d(B,s_0)=d(C,s_0)=\sqrt{10}=3.1623. \]
Perhitungan untuk \(p=2\) Bobot tak ternormalisasi: \(w_A=1/2=0.5,\; w_B=w_C=1/10=0.1\). Total = 0.7. Bobot normalisasi: \(w_A=0.7143,\; w_B=0.1429,\; w_C=0.1429\).
\[ \hat{Z}(s_0)=0.7143\cdot 10 + 0.1429\cdot 30 + 0.1429\cdot 20 = 7.143 + 4.286 + 2.857 \approx 14.29. \]
Perhitungan untuk \(p=1\)
Bobot tak ternormalisasi: \(w_A=1/1.4142=0.7071,\; w_B=w_C=1/3.1623=0.3162\). Total = 1.3395. Bobot normalisasi: \(w_A=0.528,\; w_B=0.236,\; w_C=0.236\).
\[ \hat{Z}(s_0)=0.528\cdot 10 + 0.236\cdot 30 + 0.236\cdot 20 = 5.28 + 7.08 + 4.72 = 17.08. \]
Terlihat: semakin besar \(p\), prediksi makin “menempel” pada titik terdekat (A), sehingga turun dari 17.08 (\(p=1\)) menjadi 14.29 (\(p=2\)).
Implementasi Ringkas di R
Fungsi IDW sederhana (dengan pilihan \(k\) tetangga terdekat dan penanganan \(d=0\)):
# x,y: koordinat pengamatan; z: nilai; x0,y0: target; p: power; k: k-NN (opsional)
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
d <- sqrt((x - x0)^2 + (y - y0)^2)
# Jika ada titik tepat di lokasi target
if (any(d == 0)) return(z[which.min(d)])
# Pilih k-NN jika diminta
if (!is.null(k) && k < length(d)) {
idx <- order(d)[1:k]
d <- d[idx]; z <- z[idx]
}
w <- 1 / (d^p)
sum(w * z) / sum(w)
}
# Contoh data
df <- data.frame(id = c("A","B","C"), x = c(0,4,0), y = c(0,0,4), z = c(10,30,20))
# Prediksi di s0 = (1,1)
idw_p2 <- idw_predict(df$x, df$y, df$z, 1, 1, p = 2)
idw_p1 <- idw_predict(df$x, df$y, df$z, 1, 1, p = 1)
c(p2 = idw_p2, p1 = idw_p1)## p2 p1
## 14.28571 17.08204Peta Permukaan IDW (Grid) — Visualisasi
Contoh kecil membuat permukaan IDW pada grid, lalu divisualisasikan.
library(ggplot2)
# Grid
xg <- seq(-0.5, 4.5, by = 0.1)
yg <- seq(-0.5, 4.5, by = 0.1)
G <- expand.grid(x = xg, y = yg)
# Prediksi untuk dua nilai p
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), G$x, G$y)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), G$x, G$y)
# Plot p=2 (lebih menempel ke titik terdekat)
plt2 <- ggplot(G, aes(x, y, fill = z_p2)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p2), color = "white", alpha = 0.6) +
geom_point(data = df, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text(data = df, aes(x, y, label = id), inherit.aes = FALSE, vjust = -1) +
coord_equal(xlim = c(-0.5, 4.5), ylim = c(-0.5, 4.5)) +
labs(title = "Permukaan IDW (p = 2)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
# Plot p=1 (lebih halus)
plt1 <- ggplot(G, aes(x, y, fill = z_p1)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p1), color = "white", alpha = 0.6) +
geom_point(data = df, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text(data = df, aes(x, y, label = id), inherit.aes = FALSE, vjust = -1) +
coord_equal(xlim = c(-0.5, 4.5), ylim = c(-0.5, 4.5)) +
labs(title = "Permukaan IDW (p = 1)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
plt1
plt2
Diskusi Singkat: IDW vs NN vs Kriging
- NN: ekstrem cepat, tetapi bloky/diskontinu.
- IDW: halus, parameter \(p\) mengontrol sharpness; tidak eksplisit memodelkan korelasi.
- Kriging: memerlukan model variogram; optimal (dalam asumsi), menyediakan ketidakpastian (varian kriging).
Praktik: mulai dengan IDW untuk baseline cepat; cek sensitivitas terhadap \(p\), \(k\), dan radius; jika pola spasial kuat dan data cukup, pertimbangkan kriging.
idw_p2 <- idw(z ~ 1, meuse, meuse.grid, idp = 2)## [inverse distance weighted interpolation]spplot(idw_p2["var1.pred"], main = "IDW (p=2) — Prediksi log(Zinc)")
5.7 Trend Surface / Polynomial Regression
Trend Surface Analysis (TSA) atau Polynomial Regression Surface adalah pendekatan global interpolation yang merepresentasikan variasi spasial suatu variabel (misalnya ketinggian, suhu, curah hujan) menggunakan fungsi polinomial terhadap koordinat spasial \(x\) dan \(y\):
\[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 x^2 + \beta_4 xy + \beta_5 y^2 + \cdots + \varepsilon. \]
Semakin tinggi orde polinomialnya, semakin kompleks bentuk permukaan yang bisa diwakili (misal gelombang, lengkungan, atau gradien kompleks).
Intuisi & Tujuan
- Menangkap tren global (gradien besar antar wilayah) dibandingkan fluktuasi lokal.
- Cocok untuk fenomena yang bervariasi secara halus dan kontinu (misalnya suhu rata-rata regional).
- Dapat dianggap sebagai model regresi OLS di ruang spasial, di mana \(x, y\) menjadi variabel prediktor.
Model Polinomial
1) Polinomial orde-1 (Linear) \[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon. \] Menunjukkan permukaan bidang miring (gradien konstan).
2) Polinomial orde-2 (Kuadratik) \[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 x^2 + \beta_4 xy + \beta_5 y^2 + \varepsilon. \] Menunjukkan permukaan paraboloid (melengkung); bisa menangkap cekungan, punggungan, atau lembah.
3) Polinomial orde-n (umum) \[ Z(x, y) = \sum_{i=0}^{n}\sum_{j=0}^{n-i} \beta_{ij} x^i y^j + \varepsilon. \] Makin tinggi orde, makin fleksibel, tetapi risiko overfitting juga meningkat.
Estimasi Parameter
Parameter \(\beta\) diestimasi menggunakan Ordinary Least Squares (OLS):
\[ \hat{\boldsymbol{\beta}} = (X^T X)^{-1} X^T Z, \]
- \(Z\): vektor nilai pengamatan (mis. curah hujan pada titik koordinat).
- \(X\): matriks desain yang berisi polinomial dari \(x\) dan \(y\).
- \(\varepsilon\): error (diasumsikan i.i.d., tanpa autokorelasi spasial).
Contoh Manual
Misalkan terdapat 4 titik pengamatan:
| Lokasi | x | y | Z |
|---|---|---|---|
| A | 0 | 0 | 10 |
| B | 1 | 0 | 15 |
| C | 0 | 1 | 12 |
| D | 1 | 1 | 20 |
Kita gunakan trend surface linear (orde-1): \[ Z = \beta_0 + \beta_1 x + \beta_2 y. \]
Bentuk matriksnya: \[ \begin{bmatrix} 10 \\ 15 \\ 12 \\ 20 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix}. \]
Kesimpulan: Trend surface membentuk fondasi bagi pendekatan universal kriging,
df_obs <- cbind(as.data.frame(meuse), x = coordinates(meuse)[,1], y = coordinates(meuse)[,2])
fit_trend <- lm(z ~ poly(x,2,raw=TRUE) + poly(y,2,raw=TRUE) + I(x*y), data = df_obs)
df_grd <- cbind(as.data.frame(meuse.grid), x = coordinates(meuse.grid)[,1], y = coordinates(meuse.grid)[,2])
pred_trend <- predict(fit_trend, newdata = df_grd)
trend_grid <- meuse.grid
trend_grid$z_trend <- pred_trend
spplot(trend_grid, "z_trend", main = "Trend Surface (Quadratic) — Prediksi log(Zinc)")
5.8 Thin Plate Splines (TPS)
Thin Plate Splines (TPS) adalah metode interpolasi smooth surface yang meminimalkan kelengkungan total permukaan — analogi matematisnya seperti membengkokkan lembaran logam tipis yang melewati semua titik data dengan gaya minimum.
TPS mencari fungsi \(f(x, y)\) yang:
- Melewati semua titik pengamatan \((x_i, y_i, z_i)\).
- Memiliki kelengkungan minimum, ditentukan oleh bending energy functional:
\[ J(f) = \iint \left[ \left(\frac{\partial^2 f}{\partial x^2}\right)^2 + 2\left(\frac{\partial^2 f}{\partial x \partial y}\right)^2 + \left(\frac{\partial^2 f}{\partial y^2}\right)^2 \right] dx\,dy. \]
Model Matematis
TPS dapat ditulis dalam bentuk:
\[ f(x, y) = a_0 + a_1 x + a_2 y + \sum_{i=1}^N w_i U(r_i), \]
dengan:
\[ U(r_i) = r_i^2 \log(r_i), \quad r_i = \sqrt{(x - x_i)^2 + (y - y_i)^2}. \]
Koefisien \(a_0, a_1, a_2, w_i\) diperoleh dengan menyelesaikan sistem persamaan linear yang meminimalkan bending energy sambil tetap memenuhi constraint interpolasi \(f(x_i, y_i) = z_i\).
Intuisi
- Komponen \(a_0 + a_1 x + a_2 y\) menangkap trend global linear.
- Komponen \(\sum w_i U(r_i)\) menangkap variasi lokal halus di sekitar titik data.
- Fungsi \(U(r) = r^2\log(r)\) menentukan bentuk basis radial yang memastikan permukaan halus dan kontinu hingga turunan kedua.
TPS termasuk dalam kelas Radial Basis Function (RBF) interpolators, dengan basis kernel \(\phi(r) = r^2\log(r)\).
Implementasi di R
R menyediakan fungsi Tps() dari paket
fields atau interp::interp() untuk TPS.
Contoh berikut menggunakan paket fields.
if (requireNamespace("fields", quietly = TRUE)) {
library(fields)
# Data contoh
set.seed(123)
x <- runif(10, 0, 1)
y <- runif(10, 0, 1)
z <- 5 + 2*x + 3*y + rnorm(10, 0, 0.2)
# Fit TPS model
tps_model <- fields::Tps(cbind(x, y), z)
summary(tps_model)
# Prediksi grid
grid <- expand.grid(x = seq(0, 1, length.out = 50), y = seq(0, 1, length.out = 50))
grid$z_hat <- predict(tps_model, grid)
# Visualisasi
library(ggplot2)
ggplot(grid, aes(x, y, fill = z_hat)) +
geom_raster(interpolate = TRUE) +
geom_point(data = data.frame(x = x, y = y),
mapping = aes(x = x, y = y),
inherit.aes = FALSE, size = 2, color = "black") +
geom_text(data = data.frame(x = x, y = y, z = z),
mapping = aes(x = x, y = y, label = round(z, 1)),
inherit.aes = FALSE, vjust = -1) +
scale_fill_viridis_c() +
coord_equal() +
labs(title = "Interpolasi Thin Plate Splines (TPS)", fill = "Ẑ") +
theme_minimal(base_size = 12)
}
Parameter Smoothing (λ)
TPS juga mendukung versi smoothing spline, dengan parameter \(\lambda\) yang menyeimbangkan fit dan smoothness:
\[ \text{Minimize } \sum_i [z_i - f(x_i, y_i)]^2 + \lambda J(f). \]
- \(\lambda = 0\): Exact interpolation (melewati semua titik).
- \(\lambda > 0\): permukaan lebih halus, tidak harus tepat melewati setiap titik.
Parameter \(\lambda\) sering dipilih
dengan Generalized Cross Validation (GCV) otomatis di
fungsi Tps().
Contoh Smoothing TPS
if (requireNamespace("fields", quietly = TRUE)) {
library(ggplot2)
# Pastikan grid, x, y, z tersedia; jika 'grid' sudah berisi kolom lain (mis. z_hat),
# gunakan hanya koordinat x,y saat memanggil predict.
if (!exists("x") || !exists("y") || !exists("z")) {
set.seed(123)
x <- runif(10, 0, 1); y <- runif(10, 0, 1); z <- 5 + 2*x + 3*y + rnorm(10, 0, 0.2)
}
if (!exists("grid")) {
grid <- expand.grid(x = seq(0, 1, length.out = 50), y = seq(0, 1, length.out = 50))
}
tps_smooth <- fields::Tps(cbind(x, y), z, lambda = 0.1)
coords <- grid[, c("x", "y")] # <- hanya kolom koordinat
grid$z_smooth <- predict(tps_smooth, coords)
ggplot(grid, aes(x, y, fill = z_smooth)) +
geom_raster(interpolate = TRUE) +
geom_point(data = data.frame(x = x, y = y),
mapping = aes(x = x, y = y),
inherit.aes = FALSE, size = 2, color = 'black') +
coord_equal() +
scale_fill_viridis_c() +
labs(title = "Thin Plate Splines dengan λ = 0.1 (Smoothed)", fill = "Ẑ") +
theme_minimal(base_size = 12)
}
Kelebihan & Keterbatasan
| Aspek | Penjelasan |
|---|---|
| Kelebihan | Permukaan sangat halus (hingga turunan kedua); tidak perlu parameter jarak atau tetangga; cocok untuk permukaan kontinyu. |
| Keterbatasan | Komputasi mahal untuk data besar (O(N³)); sulit mengatur anisotropy; hasil bisa oversmooth untuk gradien tajam. |
Hubungan dengan Metode Lain
- IDW: TPS lebih halus karena memperhitungkan turunan kedua.
- Polynomial trend: TPS dapat dianggap sebagai trend + komponen lokal fleksibel.
- Kriging: TPS serupa dengan kriging dengan thin-plate spline kernel (fungsi kovariansi yang setara).
# Aktifkan bila paket fields tersedia: install.packages("fields")
library(fields)
X <- coordinates(meuse)
tps_fit <- fields::Tps(X, meuse$z) # smoothing via GCV
pred_tps <- predict(tps_fit, coordinates(as(meuse.grid, "SpatialPoints")))
tps_grid <- meuse.grid
tps_grid$z_tps <- as.numeric(pred_tps)
spplot(tps_grid, "z_tps", main = "Thin Plate Splines (fields::Tps) — Prediksi log(Zinc)")
5.9 Kriging
Kriging adalah metode spatial interpolation berbasis teori proses stokastik yang mencari penduga linier tak bias dengan variansi minimum (BLUE) untuk nilai suatu peubah spasial \(Z(s)\) pada lokasi yang belum diobservasi \(s_0\). Intuisi singkatnya: gunakan tetangga yang informatif—dengan bobot optimal berbasis struktur autokorelasi—untuk menebak nilai di lokasi target.
Secara umum, penduga Kriging ditulis: \[ \hat{Z}(s_0) = \sum_{i=1}^n \lambda_i Z(s_i), \] dengan bobot \(\lambda_i\) dipilih agar tak bias dan meminimalkan varian galat prediksi.
Kekuatan Kriging ada pada variogram/kovarians: informasi jarak dan arah memengaruhi seberapa besar kontribusi tetangga terhadap prediksi.
5.10 Bidang Acak, Stasioneritas, Isotropi
5.10.1 Bidang Acak (Random Field)
Data spasial adalah realisasi dari suatu bidang acak \(Z(s)\), \(s \in \mathbb{R}^2\). Untuk dua lokasi \(s_i, s_j\), kovarians didefinisikan \[ \operatorname{Cov}[Z(s_i), Z(s_j)] = C(h), \quad h = s_i - s_j. \]
5.10.2 Stasioneritas Orde-1 dan Orde-2
- Orde-1: \(E[Z(s)] = \mu\) (konstan).
- Orde-2: \(\operatorname{Cov}[Z(s), Z(s+h)] = C(h)\) hanya fungsi lag \(h\) (bukan posisi absolut).
5.10.3 Isotropi vs Anisotropi
- Isotropi: struktur ketergantungan hanya bergantung \(\|h\|\) (panjang), tidak pada arah.
- Anisotropi: ketergantungan berubah menurut arah; dapat dimodelkan via transformasi koordinat atau variogram sektoral.
5.11 Variogram dan Kovarians
Variogram mengukur ketidaksamaan rata-rata kuadrat pada jarak \(h\): \[ \gamma(h) = \tfrac{1}{2}\operatorname{Var}[Z(s) - Z(s+h)]. \] Untuk proses stasioner orde-2, relasi dengan kovarians: \[ \gamma(h) = C(0) - C(h). \]
5.11.1 Komponen Variogram
- Nugget \(c_0\): disontinuitas di \(h\to 0\) (noise/pengukuran, variasi mikro).
- Sill \(c_0 + c\): batas nilai variogram pada jarak jauh.
- Range \(a\): jarak efektif korelasi memudar.
5.11.2 Model Variogram Teoretis
- Eksponensial: \(\gamma(h)=c_0+c\big[1-\exp(-h/a)\big]\).
- Sferis: \[\gamma(h)= \begin{cases} c_0 + c\!\left[\frac{3h}{2a} - \frac{1}{2}\!\left(\frac{h}{a}\right)^3\right], & h\le a,\\\\ c_0 + c, & h>a. \end{cases} \]
- Gaussian: \(\gamma(h)=c_0+c\big[1-\exp(-(h/a)^2)\big]\).
- (Opsional) Matérn: fleksibel; meliputi Gaussian/Eksponensial sebagai kasus batas.
5.12 Keluarga Kriging: SK, OK, UK
Tuliskan \(Z(s) = \mu(s) + \varepsilon(s)\) dengan \(E[\varepsilon(s)]=0\).
- Simple Kriging (SK): \(\mu(s)\equiv\mu\) diketahui (jarang realistis).
- Ordinary Kriging (OK): \(\mu(s)\equiv\mu\) tak diketahui tetapi konstan lokal.
- Universal Kriging (UK): \(\mu(s)=\mathbf{x}(s)^\top\beta\) (mis. polinomial koordinat).
5.13 Sistem Persamaan Kriging & Variansi Prediksi
Untuk OK, syarat tak bias: \(\sum_i \lambda_i = 1\), diberlakukan via pengali Lagrange \(\mu\): \[ \begin{bmatrix} \Sigma & \mathbf{1} \\\\ \mathbf{1}^\top & 0 \end{bmatrix} \begin{bmatrix} \lambda \\\\ \mu \end{bmatrix} = \begin{bmatrix} \sigma_0 \\\\ 1 \end{bmatrix}, \quad \sigma_K^2(s_0) = C(0) - 2\lambda^\top \sigma_0 + \lambda^\top \Sigma \lambda. \] Di mana \(\Sigma_{ij}=C(s_i-s_j)\), \(\sigma_{0,i}=C(s_i-s_0)\). Dalam bentuk variogram gunakan \(C(0)-C(h)=\gamma(h)\).
Bagian ini memberikan contoh perhitungan manual Ordinary
Kriging dengan dataset kecil (3 titik), lalu
memverifikasi hasilnya menggunakan paket
gstat di R.
Data & Target
Tiga titik pengamatan \(Z\) dan satu lokasi target \(s_0=(1,1)\):
| Titik | x | y | Z |
|---|---|---|---|
| A | 0 | 0 | 10 |
| B | 2 | 0 | 20 |
| C | 0 | 2 | 30 |
Target: \(s_0 = (1,1)\).
5.14 Model Variogram
Gunakan spherical variogram dengan parameter: - Nugget \(c_0=0\) - Sill parsial \(c=1\) - Range \(a=2\)
Bentuk fungsi variogram: \[ \gamma(h) = \begin{cases} c_0 + c\left[\tfrac{3h}{2a} - \tfrac{1}{2}\left(\tfrac{h}{a}\right)^3\right], & h \le a, \\ c_0 + c, & h > a. \end{cases} \]
Perhitungan Manual
1) Jarak antar titik & ke target
A <- c(0,0); B <- c(2,0); C <- c(0,2); s0 <- c(1,1)
edist <- function(p,q) sqrt(sum((p-q)^2))
h_AB <- edist(A,B)
h_AC <- edist(A,C)
h_BC <- edist(B,C)
h_As0 <- edist(A,s0)
h_Bs0 <- edist(B,s0)
h_Cs0 <- edist(C,s0)
c(h_AB=h_AB, h_AC=h_AC, h_BC=h_BC, h_As0=h_As0, h_Bs0=h_Bs0, h_Cs0=h_Cs0)## h_AB h_AC h_BC h_As0 h_Bs0 h_Cs0
## 2.000000 2.000000 2.828427 1.414214 1.414214 1.414214Hasil: \(h_{AB}=h_{AC}=2\); \(h_{BC}=2.828\); \(h_{A,s_0}=h_{B,s_0}=h_{C,s_0}=1.414\).
2) Fungsi variogram \(\gamma(h)\)
# Spherical variogram with c0=0, c=1, a=2
vgm_sph <- function(h, c0=0, c=1, a=2) {
out <- ifelse(h <= a, c0 + c*(1.5*(h/a) - 0.5*(h/a)^3), c0 + c)
as.numeric(out)
}
c(g_2 = vgm_sph(2), g_2_828 = vgm_sph(sqrt(8)), g_1_414 = vgm_sph(sqrt(2)))## g_2 g_2_828 g_1_414
## 1.0000000 1.0000000 0.8838835Untuk \(h=1.414\): \(\gamma \approx 1.038\); untuk \(h\ge a\), \(\gamma=1\).
3) Sistem Persamaan Kriging (Ordinary)
Sistem OK untuk 3 titik:
\[ \begin{bmatrix} 0 & \gamma_{AB} & \gamma_{AC} & 1 \\ \gamma_{BA} & 0 & \gamma_{BC} & 1 \\ \gamma_{CA} & \gamma_{CB} & 0 & 1 \\ 1 & 1 & 1 & 0 \end{bmatrix} \begin{bmatrix} \lambda_A \\ \lambda_B \\ \lambda_C \\ \mu \end{bmatrix} = \begin{bmatrix} \gamma_{A,s_0} \\ \gamma_{B,s_0} \\ \gamma_{C,s_0} \\ 1 \end{bmatrix}. \]
Dengan \(\gamma_{AB}=\gamma_{AC}=1\), \(\gamma_{BC}=1\), dan \(\gamma_{A,s_0}=\gamma_{B,s_0}=\gamma_{C,s_0}=1.038\).
G <- matrix(c(
0, 1, 1, 1,
1, 0, 1, 1,
1, 1, 0, 1,
1, 1, 1, 0
), nrow=4, byrow=TRUE)
rhs <- c(1.038, 1.038, 1.038, 1)
sol <- solve(G, rhs)
sol## [1] 0.3333333 0.3333333 0.3333333 0.3713333Komponen solusi: \(\lambda_A, \lambda_B, \lambda_C, \mu\). Periksa penjumlahan bobot = 1.
lambda <- sol[1:3]; mu <- sol[4]
sum_lambda <- sum(lambda)
Z <- c(A=10, B=20, C=30)
Zhat <- sum(lambda * Z)
list(weights=lambda, sum_lambda=sum_lambda, Zhat=Zhat, mu=mu)## $weights
## [1] 0.3333333 0.3333333 0.3333333
##
## $sum_lambda
## [1] 1
##
## $Zhat
## [1] 20
##
## $mu
## [1] 0.3713333Hasil manual: \(\lambda_A=\lambda_B=\lambda_C=1/3\), \(\hat Z(s_0)=20\).
5.15 Varian Kriging
\[ \sigma_k^2 = \sum_i \lambda_i\, \gamma(s_i, s_0) + \mu. \]
g_As0 <- vgm_sph(sqrt(2))
var_krig <- sum(lambda * c(g_As0, g_As0, g_As0)) + mu
var_krig## [1] 1.255217Varian kriging \(\approx 1.41\).
Verifikasi dengan gstat
library(sp)
library(gstat)
# Observasi sebagai SpatialPointsDataFrame
df <- data.frame(x=c(0,2,0), y=c(0,0,2), Z=c(10,20,30))
coordinates(df) <- ~x+y
# Lokasi prediksi (s0)
newd <- data.frame(x=1, y=1)
coordinates(newd) <- ~x+y
# Variogram model: Spherical, psill=1, range=2, nugget=0
vgm_model <- vgm(psill=1, model="Sph", range=2, nugget=0)
# Ordinary kriging
ok <- gstat(formula = Z~1, data=df, model=vgm_model)
kr <- predict(ok, newd, nsim=0)## [using ordinary kriging]kr## class : SpatialPointsDataFrame
## features : 1
## extent : 1, 1, 1, 1 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 2
## names : var1.pred, var1.var
## value : 20, 1.1011002862997Kolom var1.pred adalah \(\hat
Z(s_0)\) dan var1.var adalah varian kriging. Nilai
seharusnya sangat dekat dengan hasil manual (20 dan ~1.41).
Peta Permukaan Prediksi & Varian
# Grid untuk peta
xg <- seq(-0.5, 2.5, by=0.05)
yg <- seq(-0.5, 2.5, by=0.05)
newgrid <- expand.grid(x=xg, y=yg)
coordinates(newgrid) <- ~x+y
krg <- predict(ok, newgrid)## [using ordinary kriging]kr_df <- as.data.frame(krg)
library(ggplot2)
p1 <- ggplot(kr_df, aes(x, y, fill = var1.pred)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = var1.pred), color = "white", alpha = 0.6) +
geom_point(data = as.data.frame(df), aes(x, y), inherit.aes = FALSE, size = 2) +
coord_equal() +
labs(title = "Ordinary Kriging — Prediksi", fill = "Ẑ") +
theme_minimal(base_size = 12)
p2 <- ggplot(kr_df, aes(x, y, fill = var1.var)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = var1.var), color = "white", alpha = 0.6) +
geom_point(data = as.data.frame(df), aes(x, y), inherit.aes = FALSE, size = 2) +
coord_equal() +
labs(title = "Ordinary Kriging — Varian Prediksi", fill = "Var") +
theme_minimal(base_size = 12)
p1
p2
Catatan Praktis
- Pemilihan model variogram krusial: gunakan empirical variogram + fitting (OLS/WLS/REML) pada praktik nyata.
- Nugget menangkap variasi mikro/skala-pengukuran; tidak identik dengan error observasi semata.
- OK vs UK: Ordinary Kriging (mean konstan tak diketahui) vs Universal Kriging (mean fungsi kovariat/koordinat). Untuk tren kuat, pertimbangkan Regression–Kriging.
5.16 Implementasi Dasar
(CRS-Safe) — sp + gstat
# Fungsi bantu: set CRS aman ke EPSG:28992 (Amersfoort / RD New)
set_crs <- function(obj){
ok <- TRUE
tryCatch({
proj4string(obj) <- CRS("+init=epsg:28992")
}, error = function(e) ok <<- FALSE)
if(!ok){
proj4string(obj) <- CRS(SRS_string = "EPSG:28992")
}
obj
}
# Data meuse
data(meuse, package = "sp")
coordinates(meuse) <- ~x+y
meuse <- set_crs(meuse)
# Grid prediksi
data(meuse.grid, package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS(proj4string(meuse)) # samakan CRS
# Cek ringkas
head(meuse@data[,c("zinc","copper")])## zinc copper
## 1 1022 85
## 2 1141 81
## 3 640 68
## 4 257 81
## 5 269 48
## 6 281 615.16.1 Variogram Empirik & Fitting Model
v.emp <- variogram(log(zinc) ~ 1, meuse)
plot(v.emp, main = "Variogram Empirik log(Zinc)")
# Fit model Spherical (tebakan awal; silakan eksplorasi model lain)
v.fit <- fit.variogram(v.emp, vgm(psill = 1, model = "Sph", range = 900, nugget = 1))
v.fit## model psill range
## 1 Nug 0.05066243 0.0000
## 2 Sph 0.59060780 897.0209plot(v.emp, v.fit, main = "Fitting Variogram: Spherical")
5.16.2 Ordinary Kriging (OK)
ok <- krige(log(zinc) ~ 1, meuse, meuse.grid, model = v.fit)## [using ordinary kriging]spplot(ok["var1.pred"], main = "Prediksi log(Zinc) — Ordinary Kriging")
spplot(ok["var1.var"], main = "Kriging Variance — Ordinary Kriging")
5.16.3 Validasi Silang (10-fold)
set.seed(1)
cv <- krige.cv(log(zinc) ~ 1, meuse, model = v.fit, nfold = 10)
summary(cv$residual)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.927351 -0.209674 -0.007915 0.013349 0.212593 1.373407ggplot(cv@data, aes(x = observed, y = var1.pred)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
labs(x = "Observed log(Zinc)", y = "Predicted", title = "10-fold Cross-Validation (OK)")
5.17 Universal Kriging (UK)
Jika ada tren global, mis. \(\mu(s)=\beta_0+\beta_1 x + \beta_2 y\).
uk <- krige(log(zinc) ~ x + y, meuse, meuse.grid, model = v.fit)## [using universal kriging]spplot(uk["var1.pred"], main = "Universal Kriging (x + y) — Prediksi log(Zinc)")
Tips praktis: tambahkan polinomial atau kovariat lain (jarak ke sungai, elevasi) bila relevan.
5.18 Diagnostik & Praktik Baik
- Atur lebar bin variogram dengan bijak; terlalu halus → berisik, terlalu kasar → pola hilang.
- Cek anisotropi via variogram per sektor arah.
- Lihat peta residu (OK/UK) untuk tren sisa.
- Gunakan transformasi (mis. log) bila distribusi sangat skewed (seperti zinc).
5.19 Co-Kriging (Multivariat) — Linear Model of Coregionalization (LMC)
5.19.1 Konsep
Kita punya peubah target \(Z_1\) (zinc) dan peubah sekunder \(Z_2\) (copper) yang berkorelasi spasial. Bangun LMC yang memodelkan variogram masing-masing dan cross-variogram.
5.19.2 Implementasi
# Definisikan dua proses
gs <- gstat(NULL, id = "zinc", formula = log(zinc) ~ 1, data = meuse)
gs <- gstat(gs, id = "copper", formula = log(copper) ~ 1, data = meuse)
# Variogram & cross-variogram
v.cross <- variogram(gs)
plot(v.cross, main = "Variogram & Cross-Variogram (zinc, copper)")
# Fit LMC (seed Spherical)
vm <- vgm(psill = 1, model = "Sph", range = 900, nugget = 1)
lmc.fit <- fit.lmc(v.cross, gs, model = vm)
plot(v.cross, lmc.fit, main = "Fitting LMC (Spherical seed)")
# Prediksi co-kriging pada grid
ck.pred <- predict(lmc.fit, meuse.grid) # kolom: zinc.pred, copper.pred, dst.## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]spplot(ck.pred["zinc.pred"], main = "Prediksi Zinc — Co-Kriging")
Catatan: Co-Kriging efektif bila peubah sekunder informatif (korelasi spasial & cross variogram jelas), serta tersedia dengan cakupan spasial baik.
5.20 Kriging Ruang–Waktu
Untuk \(Z(s,t)\), gunakan \(C(h,u)\) dengan lag ruang \(h\) dan waktu \(u\). Model separabel: \(C(h,u)=C_s(h)C_t(u)\). Implementasi dapat
memakai gstat (variogram spatio-temporal) atau pendekatan
Bayesian (INLA) dengan efek acak ruang–waktu.
5.21 Kriging Bayesian via INLA–SPDE (Opsional)
5.21.1 Intuisi SPDE
SPDE mengaitkan GRF kontinu dan GMRF pada mesh triangulasi melalui \[ (\kappa^2 - \Delta)^{\alpha/2} x(s) = W(s). \] Keuntungan: komputasi efisien (matriks presisi sparse) dan inferensi probabilistik penuh.
5.21.2 Contoh Kode
(dibiarkan eval=FALSE)
suppressPackageStartupMessages(library(INLA))
# Mesh
mesh <- inla.mesh.2d(loc = coordinates(meuse), max.edge = c(200, 400))
plot(mesh); points(meuse, col = "red", pch = 19)
# SPDE Matérn dengan prior PC
spde <- inla.spde2.pcmatern(
mesh, alpha = 2,
prior.range = c(1000, 0.5), # P(range < 1000) = 0.5
prior.sigma = c(1, 0.01) # P(sigma > 1) = 0.01
)
# Index & proyeksi
iset <- inla.spde.make.index("spatial", n.spde = spde$n.spde)
A <- inla.spde.make.A(mesh, loc = coordinates(meuse))
# Stack
stk <- inla.stack(
data = list(y = log(meuse$zinc)),
A = list(A, 1),
effects = list(iset, data.frame(Intercept = 1)),
tag = "est"
)
# Fitting
formula <- y ~ 1 + f(spatial, model = spde)
res <- inla(
formula,
data = inla.stack.data(stk),
family = "gaussian",
control.predictor = list(A = inla.stack.A(stk), compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
summary(res)
# Posterior mean pada mesh
mean.field <- res$summary.random$spatial$mean
plot(mesh, main = "Posterior Mean Field (SPDE)"); points(meuse, pch = 20, cex = .5)5.22 Rangkuman & Latihan Singkat
Rangkuman: Variogram adalah jantung Kriging; OK/UK menangani mean lokal vs tren; Co-Kriging memanfaatkan peubah sekunder; SPDE memberi skala besar dan ketidakpastian penuh.
Latihan:
- Bandingkan RMSE CV antara OK vs UK (x+y). Manakah yang lebih baik
untuk data
meuse? - Ubah kernel variogram (Eksponensial/Gaussian/Matérn) dan evaluasi perubahannya.
- Tambahkan peubah sekunder lain (mis.
lead,cadmium) untuk Co-Kriging; bandingkan peta.
5.23 Ringkasan Kelebihan dan Keterbatasan Metode Interpolasi
Berikut perbandingan ringkas beberapa metode interpolasi spasial yang umum digunakan, termasuk Kriging sebagai pendekatan berbasis model stokastik yang paling kuat secara statistik.
5.23.1 Nearest Neighbor (NN)
Kelebihan:
- Super cepat dan mudah diimplementasikan.
- Cocok untuk baseline awal atau data dengan sebaran sangat
jarang.
Keterbatasan:
- Hasilnya “bloky” (diskontinu di batas Voronoi).
- Tidak menangkap transisi spasial yang halus.
5.23.2 TIN (Triangulated Irregular Network / Linear Interpolation)
Kelebihan:
- Menghasilkan permukaan halus di dalam setiap segitiga Delaunay.
- Konservatif: tidak memprediksi di luar convex hull data.
Keterbatasan:
- Sensitif terhadap kualitas dan distribusi triangulasi.
- Sulit diterapkan jika data sangat tidak merata.
5.23.3 Inverse Distance Weighting (IDW)
Kelebihan:
- Intuitif, parameter sederhana (pangkat jarak \(p\)).
- Memberikan hasil halus tanpa perlu model variogram.
Keterbatasan:
- Tidak menyediakan varian prediksi (ketidakpastian).
- Sensitif terhadap nilai ekstrem di sekitar lokasi target.
- Tidak mempertimbangkan struktur autokorelasi spasial.
5.23.4 Trend Surface (Polynomial Regression)
Kelebihan:
- Mudah dan cepat untuk menangkap pola global (misal gradien
timur–barat).
- Cocok untuk fenomena dengan variasi deterministik besar.
Keterbatasan:
- Mengabaikan autokorelasi residu.
- Tidak cocok untuk pola lokal yang kompleks.
5.23.5 Thin Plate Splines (TPS)
Kelebihan:
- Sangat halus dan fleksibel; parameter halus (\(\lambda\)) dapat dipilih otomatis lewat
GCV.
- Cocok untuk fenomena kontinu dan permukaan lembut.
Keterbatasan:
- Tidak menyediakan peta ketidakpastian.
- Dapat over-smooth jika \(\lambda\) tidak diatur dengan tepat.
- Komputasi cukup berat untuk dataset besar.
5.23.6 Kriging
Kelebihan:
- Penduga linier tak bias terbaik (BLUE).
- Mempertimbangkan struktur autokorelasi melalui
variogram/kovarians.
- Menghasilkan peta prediksi dan peta ketidakpastian
(varian) secara bersamaan.
- Dapat diperluas menjadi Co-Kriging, Universal Kriging, maupun Kriging
spasio-temporal.
Keterbatasan:
- Membutuhkan pemodelan variogram yang tepat dan validasi model.
- Asumsi stasioneritas sering sulit terpenuhi.
- Komputasi intensif untuk data besar (meskipun dapat diatasi dengan
pendekatan Bayesian/INLA-SPDE).
5.24 Kesimpulan Umum
|:-|:|:|:–| | NN | Non-parametrik | Cepat, sederhana | Bloky, tanpa transisi halus | | TIN | Geometrik | Halus dalam segitiga | Sensitif terhadap triangulasi | | IDW | Deterministik | Intuitif, parameter sederhana | Tidak ada varian prediksi | | Trend | Regresi | Menangkap pola global | Abaikan korelasi spasial | | TPS | Smoothing | Halus, fleksibel | Tidak sediakan varian prediksi | | Kriging | Stokastik | Optimal, sediakan ketidakpastian | Butuh variogram, berat komputasi |
6 Spatial Interpolation
6.1 Apa itu Spatial Interpolation?
Spatial interpolation adalah prosedur memperkirakan nilai variabel pada lokasi tidak teramati berdasarkan nilai di lokasi teramati di sekitarnya. Tujuannya membuat permukaan kontinu (surface) dari data titik, misalnya peta curah hujan, kualitas udara, kandungan logam tanah, atau ketinggian.
Intuisi: lokasi yang berdekatan cenderung memiliki nilai yang mirip (Tobler’s first law). Perbedaannya ada pada cara memberi bobot dan asumsi model yang digunakan untuk “menjahit” titik-titik menjadi peta kontinu.
6.2 Asumsi dan Pertimbangan Umum dalam Spatial Interpolation
Interpolasi spasial bukan sekadar masalah matematis — ia juga
memerlukan pemahaman mendalam tentang struktur spasial
data, asumsi model, dan implikasi
praktis dari setiap pendekatan.
Bagian ini menguraikan secara rinci asumsi-asumsi utama yang mendasari
sebagian besar metode interpolasi.
6.2.1 Autokorelasi Spasial
Autokorelasi spasial merupakan fondasi utama dari
semua metode interpolasi. Prinsip dasarnya dikenal sebagai Hukum
Pertama Geografi (Tobler, 1970):
> “Everything is related to everything else, but near things are more
related than distant things.”
6.2.1.1 Penjelasan Konseptual
Autokorelasi spasial menyatakan bahwa nilai suatu variabel di
lokasi tertentu cenderung mirip dengan nilai di lokasi lain yang
berdekatan secara geografis.
Jika \(Z(s)\) adalah nilai variabel di
lokasi \(s\), maka untuk dua lokasi
berdekatan \(s_i\) dan \(s_j\):
\[ \text{Cov}(Z(s_i), Z(s_j)) = C(h) \quad \text{dengan } h = \|s_i - s_j\| \]
Artinya, semakin kecil jarak \(h\), semakin besar korelasi antar nilai.
6.2.1.2 Implikasi Praktis
- Metode seperti IDW, Kriging, dan TPS mengandalkan
ide ini untuk menimbang pengaruh data sekitar.
- Jika data menunjukkan autokorelasi negatif (nilai yang berdekatan justru sangat berbeda), interpolasi konvensional menjadi tidak valid tanpa modifikasi model.
6.2.2 Stasioneritas dan Tren
6.2.2.1 Stasioneritas
Metode geostatistik klasik, terutama Kriging, mengasumsikan bahwa data berasal dari proses spasial stasioner, yakni:
- Mean konstan: \(E[Z(s)] =
\mu\), sama di seluruh wilayah.
- Kovarians hanya tergantung pada jarak: \(\text{Cov}(Z(s_i), Z(s_j)) = C(h)\), bukan pada lokasi absolut.
Namun dalam kenyataannya, banyak fenomena memiliki tren global (misal: ketinggian yang meningkat dari barat ke timur). Dalam kasus ini, asumsi stasioneritas ketat tidak terpenuhi.
6.2.2.2 Tren Global
Jika terdapat pola sistematik, misalnya perubahan rata-rata secara
spasial, pendekatan seperti Universal Kriging (UK) atau
Trend Surface Analysis lebih sesuai.
Keduanya mengizinkan mean yang bervariasi dengan lokasi:
\[ E[Z(s)] = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon(s) \]
Di mana \(\varepsilon(s)\) tetap dianggap stasioner secara lokal.
6.2.2.3 Implikasi Praktis
- Kriging biasa (OK) cocok untuk data dengan mean
konstan.
- Kriging universal (UK) cocok untuk data dengan
gradien atau tren global.
- Penting untuk memeriksa tren sebelum fitting variogram agar tidak salah tafsir korelasi spasial.
6.2.3 Skala dan Anisotropi
6.2.3.1 Skala (Scale Dependency)
Hubungan spasial bisa berbeda pada berbagai skala observasi. Misalnya, variabilitas suhu pada skala 1 km bisa berbeda drastis dengan skala 100 km. Oleh karena itu, pemilihan resolusi grid dan parameter variogram harus sesuai dengan skala analisis.
\[ C(h) = \sigma^2 \exp\left(-\frac{h}{a}\right) \]
Parameter \(a\) (range) mengatur seberapa jauh pengaruh spasial berlaku. Jika skala salah dipilih, interpolasi bisa terlalu halus atau terlalu kasar.
6.2.3.2 Anisotropi
Anisotropi terjadi ketika korelasi spasial bergantung pada arah. Misalnya, penyebaran polutan bisa lebih cepat ke arah angin dibanding arah tegak lurusnya.
Secara matematis, anisotropi dapat dimodelkan dengan transformasi jarak arah-spesifik:
\[ h' = \sqrt{ \left(\frac{x}{a_x}\right)^2 + \left(\frac{y}{a_y}\right)^2 } \]
di mana \(a_x\) dan \(a_y\) menunjukkan jangkauan korelasi yang berbeda menurut arah.
6.2.3.3 Implikasi Praktis
- Cek anisotropi sebelum fitting variogram (
vgm()di R bisa mendeteksi lewatalpha).
- Gunakan directional variogram untuk menganalisis
apakah pola berbeda menurut arah.
- Untuk fenomena seperti curah hujan atau penyebaran polusi udara, anisotropi seringkali signifikan.
6.2.4 Kerapatan dan Sebaran Sampel
Interpolasi membutuhkan titik pengamatan yang cukup rapat dan merata agar hasil representatif.
6.2.4.1 Masalah Umum
- Sparse sampling: terlalu sedikit titik menyebabkan
ketidakpastian besar di antara area kosong.
- Clustered sampling: terlalu banyak titik pada area kecil membuat model bias ke wilayah padat.
6.2.4.2 Prinsip Heuristik
- Setidaknya, setiap area berukuran “range” (dari
variogram) harus memiliki minimal satu titik data.
- Gunakan peta densitas titik untuk menilai distribusi observasi.
6.2.4.3 Implikasi Praktis
Distribusi titik menentukan pilihan metode:
- Data jarang → gunakan metode berbasis model (Kriging).
- Data rapat → metode deterministik (IDW, TPS) sering cukup baik.
6.2.5 Validasi Model
Tidak ada model interpolasi yang sempurna; karena itu validasi sangat penting.
6.2.5.1 Cross-validation
Gunakan Leave-One-Out (LOO) atau k-fold cross-validation untuk menguji seberapa baik model memprediksi data yang tidak digunakan dalam fitting.
\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\hat{Z}(s_i) - Z(s_i))^2} \]
Metrik lain yang umum:
- MAE (Mean Absolute Error): mengukur rata-rata
deviasi absolut.
- ME (Mean Error/Bias): menunjukkan kecenderungan over/under prediction.
6.2.5.2 Interpretasi Metrik
| Metrik | Interpretasi | Nilai Ideal |
|---|---|---|
| RMSE | Akurasi keseluruhan | Sekecil mungkin |
| MAE | Stabilitas prediksi | Sekecil mungkin |
| ME | Bias rata-rata | Mendekati 0 |
6.2.5.3 Validasi Visual
- Scatter plot antara nilai observasi dan
prediksi.
- Peta residual untuk melihat bias spasial.
- Variogram residual untuk memastikan tidak ada struktur spasial tersisa.
6.2.6 Kesimpulan
Asumsi-asumsi di atas tidak hanya formalitas, tetapi komponen
kunci keberhasilan interpolasi spasial.
Memahami autokorelasi, tren, anisotropi, dan distribusi data membantu
memastikan bahwa hasil interpolasi tidak hanya halus secara visual,
tetapi juga valid secara ilmiah.
“A map is only as good as the assumptions behind it.” — Adapted from Cressie (1993)
6.3 Metode Interpolasi yang Umum
Berikut kanonik metode yang sering dipakai (dari yang non-parametrik hingga berbasis model stokastik):
- Nearest Neighbor (NN): nilai lokasi terdekat; cepat tetapi kasar (peta bloky).
- Natural Neighbor / TIN: bobot berdasar Voronoi neighborhoods; halus dan konservatif (butuh lib tertentu).
- Inverse Distance Weighting (IDW): bobot \(w_i \propto d(s_0,s_i)^{-p}\) dengan orde \(p>0\); intuitif, tanpa model kovarians.
- Trend Surface / Polynomial Regression: permukaan deterministik berbasis fungsi koordinat (x, y).
- Splines (mis. Thin Plate Spline): permukaan halus
yang meminimalkan bending energy (butuh paket tambahan seperti
fields). - Kriging (SK/OK/UK): penduga linier tak bias dengan variansi minimum berbasis variogram/kovarians; varian: Simple, Ordinary, Universal.
- Co-Kriging: memanfaatkan variabel sekunder yang berkorelasi spasial.
- Radial Basis Functions (RBF): keluarga fungsi dasar radial (multikuadrik, Gaussian) untuk permukaan halus.
Pada modul ini, contoh kode difokuskan pada IDW dan Kriging dengan
gstat, karena keduanya populer dan tanpa dependensi tambahan.
6.4 Nearest Neighbor (NN)
Definisi
Metode Nearest Neighbor (NN) atau Thiessen interpolation adalah teknik paling sederhana dalam interpolasi spasial. Untuk suatu lokasi target \(s_0\), kita cari satu titik pengamatan terdekat \(s_{(1)}\), lalu nilai prediksinya langsung diambil dari titik itu:
\[ \hat{Z}(s_0) = Z\big(s_{(1)}\big). \]
Dengan kata lain, titik yang paling dekat “memenangkan” wilayahnya. Secara geometris, seluruh ruang terbagi menjadi poligon-poligon yang disebut Voronoi cells — setiap sel mencakup semua lokasi yang lebih dekat ke satu titik pengamatan dibandingkan ke titik lain.
Contoh Manual
Misalkan kita punya tiga titik pengamatan suhu udara (°C):
| Lokasi | Koordinat (x, y) | Nilai Z |
|---|---|---|
| A | (0, 0) | 28 |
| B | (2, 0) | 30 |
| C | (0, 3) | 25 |
Kita ingin memperkirakan suhu di titik baru \(s_0 = (1, 1)\).
Langkah Perhitungan
Hitung jarak Euclidean ke tiap titik:
\[ \begin{aligned} d(A, s_0) &= \sqrt{(1-0)^2 + (1-0)^2} = \sqrt{2} = 1.414, \\ d(B, s_0) &= \sqrt{(1-2)^2 + (1-0)^2} = \sqrt{2} = 1.414, \\ d(C, s_0) &= \sqrt{(1-0)^2 + (1-3)^2} = \sqrt{10} = 3.162. \end{aligned} \]
Titik A dan B sama-sama paling dekat.
Jika ada jarak sama (tie), ambil salah satu secara deterministik (misalnya yang pertama dalam urutan data), atau gunakan rata-rata dari kandidat terdekat jika ingin hasil lebih halus. Di sini kita ambil A.
Prediksi nilai di \(s_0\):
\[ \hat{Z}(s_0) = Z(A) = 28. \]
Implementasi Ringkas di R
Kode berikut mengilustrasikan perhitungan jarak, pemilihan tetangga terdekat, dan (opsional) sketsa mosaik Voronoi untuk tiga titik tersebut.
# Data titik pengamatan
pts <- data.frame(
id = c("A", "B", "C"),
x = c(0, 2, 0),
y = c(0, 0, 3),
Z = c(28, 30, 25)
)
# Titik target
s0 <- data.frame(x = 1, y = 1)
# Fungsi jarak Euclidean
edist <- function(x1, y1, x2, y2) sqrt((x1 - x2)^2 + (y1 - y2)^2)
# Hitung jarak ke setiap titik pengamatan
pts$dist_s0 <- with(pts, edist(x, y, s0$x, s0$y))
# Urutkan berdasarkan jarak
pts_ord <- pts[order(pts$dist_s0), ]
pts_ord## id x y Z dist_s0
## 1 A 0 0 28 1.414214
## 2 B 2 0 30 1.414214
## 3 C 0 3 25 2.236068# Ambil tetangga terdekat (tie-breaker: urutan data)
nearest <- pts_ord[1, ]
cat(sprintf("Nearest to s0=(%.1f,%.1f) is %s at distance %.3f with Z=%g\n",
s0$x, s0$y, nearest$id, nearest$dist_s0, nearest$Z))## Nearest to s0=(1.0,1.0) is A at distance 1.414 with Z=28Zhat <- nearest$Z
Zhat## [1] 28# Visualisasi Voronoi (opsional)
# Memerlukan paket deldir dan ggplot2
if (requireNamespace("deldir", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)) {
library(deldir)
library(ggplot2)
# Batas plot
xlim <- c(-0.5, 2.5)
ylim <- c(-0.5, 3.5)
# Delaunay/Voronoi tessellation
dd <- deldir::deldir(pts$x, pts$y, rw = c(xlim, ylim))
tiles <- deldir::tile.list(dd)
# Konversi tile ke data.frame untuk ggplot
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(
id = pts$id[i],
x = tiles[[i]]$x,
y = tiles[[i]]$y
)
}))
ggplot() +
geom_polygon(data = poly_df, aes(x = x, y = y, group = id),
fill = NA, color = "grey40") +
geom_point(data = pts, aes(x = x, y = y), size = 3) +
geom_text(data = pts, aes(x = x, y = y, label = id), vjust = -1) +
geom_point(data = s0, aes(x = x, y = y), shape = 4, size = 3, stroke = 1.2) +
coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
labs(title = "Sketsa Voronoi untuk Nearest Neighbor",
subtitle = "Titik s0 ditandai dengan silang (×)",
x = "x", y = "y") +
theme_minimal(base_size = 12)
} else {
message("Lewati plot: paket 'deldir' dan/atau 'ggplot2' belum terpasang.")
}
Sifat, Kelebihan, dan Kelemahan
- Kelebihan: sangat cepat, tanpa parameter, dan tidak membutuhkan variogram.
- Kelemahan: hasil tidak halus (blocky), tidak cocok untuk fenomena yang berubah gradual, dan sensitif terhadap distribusi titik.
Metode ini sering dipakai untuk visualisasi cepat (misalnya Thiessen polygons pada data curah hujan) sebelum beralih ke metode berbobot seperti IDW atau Kriging.
Contoh2
# Data meuse
data(meuse, package = "sp")
coordinates(meuse) <- ~x+y
meuse <- set_crs(meuse)
meuse$z <- log(meuse$zinc)
data(meuse.grid, package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS(proj4string(meuse))
# Sanity checks
stopifnot(!is.na(proj4string(meuse)), !is.na(proj4string(meuse.grid)))
# Cek ringkas
head(meuse@data[,c("zinc","copper")])## zinc copper
## 1 1022 85
## 2 1141 81
## 3 640 68
## 4 257 81
## 5 269 48
## 6 281 61# 1) Ambil koordinat observasi & grid (sebagai SpatialPoints)
Xobs <- coordinates(meuse)
grid_sp <- as(meuse.grid, "SpatialPoints")
# 2) Cari indeks tetangga terdekat untuk setiap titik grid
knn <- FNN::get.knnx(data = Xobs, query = coordinates(grid_sp), k = 1)
nn_ix <- knn$nn.index[,1]
# 3) Prediksi NN dan bangun SpatialPixelsDataFrame baru
pred_nn <- meuse$z[nn_ix]
stopifnot(length(pred_nn) == nrow(as.data.frame(meuse.grid)))
nn_grid <- SpatialPixelsDataFrame(points = grid_sp,
data = data.frame(z_nn = pred_nn),
proj4string = CRS(proj4string(meuse.grid)))
# 4) Plot aman — kolom 'z_nn' dijamin ada
spplot(nn_grid, "z_nn", main = "Nearest Neighbor (1-NN) — Prediksi log(Zinc)")
6.5 Natural Neighbor / TIN
Metode Natural Neighbor (atau Sibson interpolation) bekerja berdasarkan Voronoi diagram juga, namun dengan gagasan yang lebih halus. Jika Nearest Neighbor (NN) hanya mengambil satu titik terdekat, Natural Neighbor mempertimbangkan beberapa titik yang “terkait” dengan lokasi baru \(s_0\).
Langkah pikirannya:
- Kita punya sekumpulan titik data (mis. A, B, C, D).
- Dari titik-titik itu kita buat Voronoi diagram — setiap titik punya wilayahnya.
- Tambahkan lokasi baru \(s_0\). Diagram Voronoi berubah: sebagian area dari poligon-poligon tetangga “direbut” oleh \(s_0\).
- Luas potongan yang direbut dari masing-masing tetangga menjadi bobot interpolasi.
Secara konseptual:
\[ \hat{Z}(s_0) = \sum_i w_i \, Z(s_i), \quad \text{dengan } w_i \ge 0, \; \sum_i w_i = 1. \]
Bobot \(w_i\) sebanding dengan proporsi luas Voronoi yang “diambil alih” oleh \(s_0\). Karena bobot positif dan menjumlah ke 1, hasil interpolasi bersifat halus dan kontinu.
6.5.1 Hubungan dengan TIN (Triangulated Irregular Network)
Implementasi Sibson interpolation yang eksak cukup kompleks dan jarang tersedia langsung di R. Alternatif praktis adalah TIN linear interpolation: interpolasi linier di dalam Delaunay triangulation (dual dari Voronoi).
Intinya:
- Bentuk triangulasi Delaunay dari titik-titik data.
- Temukan segitiga yang mengandung \(s_0\).
- Hitung bobot barycentric (perbandingan luas sub-segitiga) terhadap tiga titik sudut.
- Nilai interpolasi adalah rata-rata tertimbang linier dengan bobot barycentric.
6.5.2 Contoh Manual — TIN Linear (Barycentric)
Misalkan tiga titik pengamatan membentuk satu segitiga:
| Lokasi | (x, y) | Nilai Z |
|---|---|---|
| A | (0, 0) | 10 |
| B | (4, 0) | 30 |
| C | (0, 4) | 20 |
Kita ingin menghitung nilai di \(s_0 = (1, 1)\).
1) Cek posisi \(s_0\)
\(s_0\) berada di dalam segitiga ABC karena koordinatnya positif dan \(1 + 1 < 4\).
2) Luas segitiga utama (ABC)
Gunakan rumus determinan (setengah dari determinan):
\[ L_{ABC} = \tfrac{1}{2} \big| x_A(y_B - y_C) + x_B(y_C - y_A) + x_C(y_A - y_B) \big|. \]
Substitusi koordinat menghasilkan \(L_{ABC} = 8\).
3) Luas tiga sub-segitiga yang melibatkan \(s_0\)
- \(L_{s_0BC} = \tfrac{1}{2} \big| 1(0-4) + 4(4-1) + 0(1-0) \big| = 4\)
- \(L_{As_0C} = \tfrac{1}{2} \big| 0(1-4) + 1(4-0) + 0(0-1) \big| = 2\)
- \(L_{ABs_0} = \tfrac{1}{2} \big| 0(0-1) + 4(1-0) + 1(0-0) \big| = 2\)
Total luas sub-segitiga = \(4 + 2 + 2 = 8\), konsisten dengan \(L_{ABC}\).
4) Bobot barycentric
\[ w_A = \frac{L_{s_0BC}}{L_{ABC}} = \frac{4}{8} = 0.5, \quad w_B = \frac{L_{As_0C}}{L_{ABC}} = \frac{2}{8} = 0.25, \quad w_C = \frac{L_{ABs_0}}{L_{ABC}} = \frac{2}{8} = 0.25. \]
5) Interpolasi nilai
\[ \hat{Z}(s_0) = 0.5\cdot 10 + 0.25\cdot 30 + 0.25\cdot 20 = 17.5. \]
Implementasi Ringkas di R
Contoh fungsi kecil untuk menghitung bobot barycentric dan nilai interpolasi pada \(s_0\):
tri_barycentric <- function(A, B, C, s0, Z) {
area <- function(P, Q, R) 0.5 * abs(P[1]*(Q[2]-R[2]) + Q[1]*(R[2]-P[2]) + R[1]*(P[2]-Q[2]))
L_ABC <- area(A, B, C)
L_s0BC <- area(s0, B, C)
L_As0C <- area(A, s0, C)
L_ABs0 <- area(A, B, s0)
wA <- L_s0BC / L_ABC
wB <- L_As0C / L_ABC
wC <- L_ABs0 / L_ABC
Zhat <- wA*Z[1] + wB*Z[2] + wC*Z[3]
list(weights = c(A = wA, B = wB, C = wC), Zhat = Zhat,
areas = c(ABC = L_ABC, s0BC = L_s0BC, As0C = L_As0C, ABs0 = L_ABs0))
}
A <- c(0, 0); B <- c(4, 0); C <- c(0, 4)
s0 <- c(1, 1)
Z <- c(10, 30, 20) # Z_A, Z_B, Z_C
res <- tri_barycentric(A, B, C, s0, Z)
res## $weights
## A B C
## 0.50 0.25 0.25
##
## $Zhat
## [1] 17.5
##
## $areas
## ABC s0BC As0C ABs0
## 8 4 2 2Visualisasi Voronoi Diagram
Visualisasi berikut memperlihatkan sel Voronoi untuk titik A, B, dan C, dengan \(s_0 = (1,1)\) ditandai silang (×).
if (requireNamespace("deldir", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)) {
library(deldir)
library(ggplot2)
pts <- data.frame(id = c("A","B","C"), x = c(0,4,0), y = c(0,0,4))
s0 <- data.frame(x = 1, y = 1)
dd <- deldir::deldir(pts$x, pts$y, rw = c(-0.5,4.5,-0.5,4.5))
tiles <- deldir::tile.list(dd)
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(id = pts$id[i], x = tiles[[i]]$x, y = tiles[[i]]$y)
}))
ggplot() +
geom_polygon(data = poly_df, aes(x=x, y=y, group=id), fill=NA, color="grey40") +
geom_point(data = pts, aes(x=x, y=y), size=3) +
geom_text(data = pts, aes(x=x, y=y, label=id), vjust=-1) +
geom_point(data = s0, aes(x=x, y=y), shape=4, size=3, stroke=1.2) +
coord_equal(xlim=c(-0.5,4.5), ylim=c(-0.5,4.5)) +
labs(title="Voronoi Diagram untuk A(0,0), B(4,0), C(0,4)",
subtitle="Titik s0=(1,1) ditandai dengan silang (×)",
x="x", y="y") +
theme_minimal(base_size = 12)
}
Catatan: Untuk implementasi TIN skala lebih besar beserta triangulasi Delaunay, Anda dapat menggunakan paket seperti
RTriangle,geometry(fungsidelaunayn), ataudeldir(untuk Voronoi/Delaunay 2D), lalu menerapkan interpolasi linier di dalam setiap segitiga.
Contoh 2
Xobs <- coordinates(meuse)
xg <- sort(unique(meuse.grid@coords[,1]))
yg <- sort(unique(meuse.grid@coords[,2]))
tin_out <- interp::interp(x = Xobs[,1], y = Xobs[,2], z = meuse$z,
xo = xg, yo = yg, duplicate = "mean", linear = TRUE)
zmat <- tin_out$z
grd <- expand.grid(x = tin_out$x, y = tin_out$y)
tin_sgdf <- SpatialPixelsDataFrame(points = SpatialPoints(grd),
data = data.frame(z_tin = as.vector(zmat)),
tolerance = 0.5)
proj4string(tin_sgdf) <- CRS(proj4string(meuse))
spplot(tin_sgdf, "z_tin", main = "TIN (Linear) — Perkiraan log(Zinc)")
xg <- sort(unique(meuse.grid@coords[,1]))
yg <- sort(unique(meuse.grid@coords[,2]))
tin_out <- interp::interp(x = Xobs[,1], y = Xobs[,2], z = meuse$z,
xo = xg, yo = yg, duplicate = "mean", linear = TRUE)
zmat <- tin_out$z
grd <- expand.grid(x = tin_out$x, y = tin_out$y)
tin_sgdf <- SpatialPixelsDataFrame(points = SpatialPoints(grd),
data = data.frame(z_tin = as.vector(zmat)),
tolerance = 0.5)
proj4string(tin_sgdf) <- CRS(proj4string(meuse))
spplot(tin_sgdf["z_tin"], main = "TIN (Linear Barycentric) — Perkiraan log(Zinc)")
6.6 Inverse Distance Weighting (IDW)
Inverse Distance Weighting (IDW) adalah metode interpolasi deterministik yang memperkirakan nilai di lokasi target \(s_0\) sebagai rata-rata tertimbang dari nilai-nilai di titik pengamatan \(s_i\), dengan bobot menurun terhadap jarak.
\[ w_i(s_0) \propto d(s_0, s_i)^{-p}, \qquad p > 0, \]
\[ \hat{Z}(s_0) = \frac{\sum_i w_i(s_0)\, Z(s_i)}{\sum_i w_i(s_0)}. \]
Umumnya kita gunakan jarak Euclidean \(d(\cdot,\cdot)\). Parameter power \(p\) mengontrol tingkat pengaruh jarak: makin besar \(p\), makin tajam penurunan bobot.
Sifat & Intuisi Parameter \(p\)
- \(p \to 0\): semua bobot hampir sama (mendekati rata-rata biasa).
- \(p=1\): penurunan bobot linear terhadap jarak (dalam skala log–log).
- \(p=2\): penurunan kuadratik; pilihan populer untuk fenomena yang meredup cepat.
- \(p \to \infty\): IDW mendekati Nearest Neighbor (hanya tetangga terdekat yang dominan).
IDW tidak memerlukan variogram dan tidak memodelkan struktur kovariansi. Ia lokal, cepat, dan deterministik; namun tidak optimal secara statistik (tidak BLUE seperti kriging ketika asumsi terpenuhi).
Kebijakan Tetangga (Neighborhood)
- Semua titik (global) – lebih halus, namun bisa bias bila distribusi titik tidak merata.
- \(k\) tetangga terdekat – stabil di tepi; pilih \(k\) (mis. 6–20) sesuai densitas sampel.
- Radius pencarian – hanya titik dalam radius \(r\) yang digunakan.
Catatan tepi (edge cases): - Jika ada titik tepat di \(s_0\) (\(d=0\)), definisikan \(\hat{Z}(s_0) = Z(s_i)\) untuk menghindari pembagian nol. - Titik duplikat/kolokasi sebaiknya diagregasi (mis. rata-rata) atau ditangani sebagai pengamatan terpisah dengan jarak \(>0\) sangat kecil. - Anisotropi bisa dicapai dengan metrik jarak berbobot/ter-skala (mis. skala sumbu x,y berbeda) atau rotasi sumbu.
Contoh Manual (Data Segitiga)
Tiga titik pengamatan:
| Lokasi | (x, y) | Z |
|---|---|---|
| A | (0, 0) | 10 |
| B | (4, 0) | 30 |
| C | (0, 4) | 20 |
Target: \(s_0 = (1, 1)\). Jarak Euclidean:
\[ d(A,s_0)=\sqrt{2}=1.4142,\quad d(B,s_0)=d(C,s_0)=\sqrt{10}=3.1623. \]
Perhitungan untuk \(p=2\) Bobot tak ternormalisasi: \(w_A=1/2=0.5,\; w_B=w_C=1/10=0.1\). Total = 0.7. Bobot normalisasi: \(w_A=0.7143,\; w_B=0.1429,\; w_C=0.1429\).
\[ \hat{Z}(s_0)=0.7143\cdot 10 + 0.1429\cdot 30 + 0.1429\cdot 20 = 7.143 + 4.286 + 2.857 \approx 14.29. \]
Perhitungan untuk \(p=1\)
Bobot tak ternormalisasi: \(w_A=1/1.4142=0.7071,\; w_B=w_C=1/3.1623=0.3162\). Total = 1.3395. Bobot normalisasi: \(w_A=0.528,\; w_B=0.236,\; w_C=0.236\).
\[ \hat{Z}(s_0)=0.528\cdot 10 + 0.236\cdot 30 + 0.236\cdot 20 = 5.28 + 7.08 + 4.72 = 17.08. \]
Terlihat: semakin besar \(p\), prediksi makin “menempel” pada titik terdekat (A), sehingga turun dari 17.08 (\(p=1\)) menjadi 14.29 (\(p=2\)).
Implementasi Ringkas di R
Fungsi IDW sederhana (dengan pilihan \(k\) tetangga terdekat dan penanganan \(d=0\)):
# x,y: koordinat pengamatan; z: nilai; x0,y0: target; p: power; k: k-NN (opsional)
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
d <- sqrt((x - x0)^2 + (y - y0)^2)
# Jika ada titik tepat di lokasi target
if (any(d == 0)) return(z[which.min(d)])
# Pilih k-NN jika diminta
if (!is.null(k) && k < length(d)) {
idx <- order(d)[1:k]
d <- d[idx]; z <- z[idx]
}
w <- 1 / (d^p)
sum(w * z) / sum(w)
}
# Contoh data
df <- data.frame(id = c("A","B","C"), x = c(0,4,0), y = c(0,0,4), z = c(10,30,20))
# Prediksi di s0 = (1,1)
idw_p2 <- idw_predict(df$x, df$y, df$z, 1, 1, p = 2)
idw_p1 <- idw_predict(df$x, df$y, df$z, 1, 1, p = 1)
c(p2 = idw_p2, p1 = idw_p1)## p2 p1
## 14.28571 17.08204Peta Permukaan IDW (Grid) — Visualisasi
Contoh kecil membuat permukaan IDW pada grid, lalu divisualisasikan.
library(ggplot2)
# Grid
xg <- seq(-0.5, 4.5, by = 0.1)
yg <- seq(-0.5, 4.5, by = 0.1)
G <- expand.grid(x = xg, y = yg)
# Prediksi untuk dua nilai p
G$z_p2 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 2), G$x, G$y)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df$x, df$y, df$z, x0, y0, p = 1), G$x, G$y)
# Plot p=2 (lebih menempel ke titik terdekat)
plt2 <- ggplot(G, aes(x, y, fill = z_p2)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p2), color = "white", alpha = 0.6) +
geom_point(data = df, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text(data = df, aes(x, y, label = id), inherit.aes = FALSE, vjust = -1) +
coord_equal(xlim = c(-0.5, 4.5), ylim = c(-0.5, 4.5)) +
labs(title = "Permukaan IDW (p = 2)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
# Plot p=1 (lebih halus)
plt1 <- ggplot(G, aes(x, y, fill = z_p1)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p1), color = "white", alpha = 0.6) +
geom_point(data = df, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text(data = df, aes(x, y, label = id), inherit.aes = FALSE, vjust = -1) +
coord_equal(xlim = c(-0.5, 4.5), ylim = c(-0.5, 4.5)) +
labs(title = "Permukaan IDW (p = 1)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
plt1
plt2
Diskusi Singkat: IDW vs NN vs Kriging
- NN: ekstrem cepat, tetapi bloky/diskontinu.
- IDW: halus, parameter \(p\) mengontrol sharpness; tidak eksplisit memodelkan korelasi.
- Kriging: memerlukan model variogram; optimal (dalam asumsi), menyediakan ketidakpastian (varian kriging).
Praktik: mulai dengan IDW untuk baseline cepat; cek sensitivitas terhadap \(p\), \(k\), dan radius; jika pola spasial kuat dan data cukup, pertimbangkan kriging.
idw_p2 <- idw(z ~ 1, meuse, meuse.grid, idp = 2)## [inverse distance weighted interpolation]spplot(idw_p2["var1.pred"], main = "IDW (p=2) — Prediksi log(Zinc)")
6.7 Trend Surface / Polynomial Regression
Trend Surface Analysis (TSA) atau Polynomial Regression Surface adalah pendekatan global interpolation yang merepresentasikan variasi spasial suatu variabel (misalnya ketinggian, suhu, curah hujan) menggunakan fungsi polinomial terhadap koordinat spasial \(x\) dan \(y\):
\[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 x^2 + \beta_4 xy + \beta_5 y^2 + \cdots + \varepsilon. \]
Semakin tinggi orde polinomialnya, semakin kompleks bentuk permukaan yang bisa diwakili (misal gelombang, lengkungan, atau gradien kompleks).
Intuisi & Tujuan
- Menangkap tren global (gradien besar antar wilayah) dibandingkan fluktuasi lokal.
- Cocok untuk fenomena yang bervariasi secara halus dan kontinu (misalnya suhu rata-rata regional).
- Dapat dianggap sebagai model regresi OLS di ruang spasial, di mana \(x, y\) menjadi variabel prediktor.
Model Polinomial
1) Polinomial orde-1 (Linear) \[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon. \] Menunjukkan permukaan bidang miring (gradien konstan).
2) Polinomial orde-2 (Kuadratik) \[ Z(x, y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 x^2 + \beta_4 xy + \beta_5 y^2 + \varepsilon. \] Menunjukkan permukaan paraboloid (melengkung); bisa menangkap cekungan, punggungan, atau lembah.
3) Polinomial orde-n (umum) \[ Z(x, y) = \sum_{i=0}^{n}\sum_{j=0}^{n-i} \beta_{ij} x^i y^j + \varepsilon. \] Makin tinggi orde, makin fleksibel, tetapi risiko overfitting juga meningkat.
Estimasi Parameter
Parameter \(\beta\) diestimasi menggunakan Ordinary Least Squares (OLS):
\[ \hat{\boldsymbol{\beta}} = (X^T X)^{-1} X^T Z, \]
- \(Z\): vektor nilai pengamatan (mis. curah hujan pada titik koordinat).
- \(X\): matriks desain yang berisi polinomial dari \(x\) dan \(y\).
- \(\varepsilon\): error (diasumsikan i.i.d., tanpa autokorelasi spasial).
Contoh Manual
Misalkan terdapat 4 titik pengamatan:
| Lokasi | x | y | Z |
|---|---|---|---|
| A | 0 | 0 | 10 |
| B | 1 | 0 | 15 |
| C | 0 | 1 | 12 |
| D | 1 | 1 | 20 |
Kita gunakan trend surface linear (orde-1): \[ Z = \beta_0 + \beta_1 x + \beta_2 y. \]
Bentuk matriksnya: \[ \begin{bmatrix} 10 \\ 15 \\ 12 \\ 20 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix}. \]
Kesimpulan: Trend surface membentuk fondasi bagi pendekatan universal kriging,
df_obs <- cbind(as.data.frame(meuse), x = coordinates(meuse)[,1], y = coordinates(meuse)[,2])
fit_trend <- lm(z ~ poly(x,2,raw=TRUE) + poly(y,2,raw=TRUE) + I(x*y), data = df_obs)
df_grd <- cbind(as.data.frame(meuse.grid), x = coordinates(meuse.grid)[,1], y = coordinates(meuse.grid)[,2])
pred_trend <- predict(fit_trend, newdata = df_grd)
trend_grid <- meuse.grid
trend_grid$z_trend <- pred_trend
spplot(trend_grid, "z_trend", main = "Trend Surface (Quadratic) — Prediksi log(Zinc)")
6.8 Thin Plate Splines (TPS)
Thin Plate Splines (TPS) adalah metode interpolasi smooth surface yang meminimalkan kelengkungan total permukaan — analogi matematisnya seperti membengkokkan lembaran logam tipis yang melewati semua titik data dengan gaya minimum.
TPS mencari fungsi \(f(x, y)\) yang:
- Melewati semua titik pengamatan \((x_i, y_i, z_i)\).
- Memiliki kelengkungan minimum, ditentukan oleh bending energy functional:
\[ J(f) = \iint \left[ \left(\frac{\partial^2 f}{\partial x^2}\right)^2 + 2\left(\frac{\partial^2 f}{\partial x \partial y}\right)^2 + \left(\frac{\partial^2 f}{\partial y^2}\right)^2 \right] dx\,dy. \]
Model Matematis
TPS dapat ditulis dalam bentuk:
\[ f(x, y) = a_0 + a_1 x + a_2 y + \sum_{i=1}^N w_i U(r_i), \]
dengan:
\[ U(r_i) = r_i^2 \log(r_i), \quad r_i = \sqrt{(x - x_i)^2 + (y - y_i)^2}. \]
Koefisien \(a_0, a_1, a_2, w_i\) diperoleh dengan menyelesaikan sistem persamaan linear yang meminimalkan bending energy sambil tetap memenuhi constraint interpolasi \(f(x_i, y_i) = z_i\).
Intuisi
- Komponen \(a_0 + a_1 x + a_2 y\) menangkap trend global linear.
- Komponen \(\sum w_i U(r_i)\) menangkap variasi lokal halus di sekitar titik data.
- Fungsi \(U(r) = r^2\log(r)\) menentukan bentuk basis radial yang memastikan permukaan halus dan kontinu hingga turunan kedua.
TPS termasuk dalam kelas Radial Basis Function (RBF) interpolators, dengan basis kernel \(\phi(r) = r^2\log(r)\).
Implementasi di R
R menyediakan fungsi Tps() dari paket
fields atau interp::interp() untuk TPS.
Contoh berikut menggunakan paket fields.
if (requireNamespace("fields", quietly = TRUE)) {
library(fields)
# Data contoh
set.seed(123)
x <- runif(10, 0, 1)
y <- runif(10, 0, 1)
z <- 5 + 2*x + 3*y + rnorm(10, 0, 0.2)
# Fit TPS model
tps_model <- fields::Tps(cbind(x, y), z)
summary(tps_model)
# Prediksi grid
grid <- expand.grid(x = seq(0, 1, length.out = 50), y = seq(0, 1, length.out = 50))
grid$z_hat <- predict(tps_model, grid)
# Visualisasi
library(ggplot2)
ggplot(grid, aes(x, y, fill = z_hat)) +
geom_raster(interpolate = TRUE) +
geom_point(data = data.frame(x = x, y = y),
mapping = aes(x = x, y = y),
inherit.aes = FALSE, size = 2, color = "black") +
geom_text(data = data.frame(x = x, y = y, z = z),
mapping = aes(x = x, y = y, label = round(z, 1)),
inherit.aes = FALSE, vjust = -1) +
scale_fill_viridis_c() +
coord_equal() +
labs(title = "Interpolasi Thin Plate Splines (TPS)", fill = "Ẑ") +
theme_minimal(base_size = 12)
}
Parameter Smoothing (λ)
TPS juga mendukung versi smoothing spline, dengan parameter \(\lambda\) yang menyeimbangkan fit dan smoothness:
\[ \text{Minimize } \sum_i [z_i - f(x_i, y_i)]^2 + \lambda J(f). \]
- \(\lambda = 0\): Exact interpolation (melewati semua titik).
- \(\lambda > 0\): permukaan lebih halus, tidak harus tepat melewati setiap titik.
Parameter \(\lambda\) sering dipilih
dengan Generalized Cross Validation (GCV) otomatis di
fungsi Tps().
Contoh Smoothing TPS
if (requireNamespace("fields", quietly = TRUE)) {
library(ggplot2)
# Pastikan grid, x, y, z tersedia; jika 'grid' sudah berisi kolom lain (mis. z_hat),
# gunakan hanya koordinat x,y saat memanggil predict.
if (!exists("x") || !exists("y") || !exists("z")) {
set.seed(123)
x <- runif(10, 0, 1); y <- runif(10, 0, 1); z <- 5 + 2*x + 3*y + rnorm(10, 0, 0.2)
}
if (!exists("grid")) {
grid <- expand.grid(x = seq(0, 1, length.out = 50), y = seq(0, 1, length.out = 50))
}
tps_smooth <- fields::Tps(cbind(x, y), z, lambda = 0.1)
coords <- grid[, c("x", "y")] # <- hanya kolom koordinat
grid$z_smooth <- predict(tps_smooth, coords)
ggplot(grid, aes(x, y, fill = z_smooth)) +
geom_raster(interpolate = TRUE) +
geom_point(data = data.frame(x = x, y = y),
mapping = aes(x = x, y = y),
inherit.aes = FALSE, size = 2, color = 'black') +
coord_equal() +
scale_fill_viridis_c() +
labs(title = "Thin Plate Splines dengan λ = 0.1 (Smoothed)", fill = "Ẑ") +
theme_minimal(base_size = 12)
}
Kelebihan & Keterbatasan
| Aspek | Penjelasan |
|---|---|
| Kelebihan | Permukaan sangat halus (hingga turunan kedua); tidak perlu parameter jarak atau tetangga; cocok untuk permukaan kontinyu. |
| Keterbatasan | Komputasi mahal untuk data besar (O(N³)); sulit mengatur anisotropy; hasil bisa oversmooth untuk gradien tajam. |
Hubungan dengan Metode Lain
- IDW: TPS lebih halus karena memperhitungkan turunan kedua.
- Polynomial trend: TPS dapat dianggap sebagai trend + komponen lokal fleksibel.
- Kriging: TPS serupa dengan kriging dengan thin-plate spline kernel (fungsi kovariansi yang setara).
# Aktifkan bila paket fields tersedia: install.packages("fields")
library(fields)
X <- coordinates(meuse)
tps_fit <- fields::Tps(X, meuse$z) # smoothing via GCV
pred_tps <- predict(tps_fit, coordinates(as(meuse.grid, "SpatialPoints")))
tps_grid <- meuse.grid
tps_grid$z_tps <- as.numeric(pred_tps)
spplot(tps_grid, "z_tps", main = "Thin Plate Splines (fields::Tps) — Prediksi log(Zinc)")
6.9 Kriging
Kriging adalah metode spatial interpolation berbasis teori proses stokastik yang mencari penduga linier tak bias dengan variansi minimum (BLUE) untuk nilai suatu peubah spasial \(Z(s)\) pada lokasi yang belum diobservasi \(s_0\). Intuisi singkatnya: gunakan tetangga yang informatif—dengan bobot optimal berbasis struktur autokorelasi—untuk menebak nilai di lokasi target.
Secara umum, penduga Kriging ditulis: \[ \hat{Z}(s_0) = \sum_{i=1}^n \lambda_i Z(s_i), \] dengan bobot \(\lambda_i\) dipilih agar tak bias dan meminimalkan varian galat prediksi.
Kekuatan Kriging ada pada variogram/kovarians: informasi jarak dan arah memengaruhi seberapa besar kontribusi tetangga terhadap prediksi.
6.9.1 Bidang Acak, Stasioneritas, Isotropi
6.9.2 Bidang Acak (Random Field)
Data spasial adalah realisasi dari suatu bidang acak \(Z(s)\), \(s \in \mathbb{R}^2\). Untuk dua lokasi \(s_i, s_j\), kovarians didefinisikan \[ \operatorname{Cov}[Z(s_i), Z(s_j)] = C(h), \quad h = s_i - s_j. \]
6.9.2.1 Stasioneritas Orde-1 dan Orde-2
- Orde-1: \(E[Z(s)] = \mu\) (konstan).
- Orde-2: \(\operatorname{Cov}[Z(s), Z(s+h)] = C(h)\) hanya fungsi lag \(h\) (bukan posisi absolut).
6.9.2.2 Isotropi vs Anisotropi
- Isotropi: struktur ketergantungan hanya bergantung \(\|h\|\) (panjang), tidak pada arah.
- Anisotropi: ketergantungan berubah menurut arah; dapat dimodelkan via transformasi koordinat atau variogram sektoral.
6.9.2.3 Variogram dan Kovarians
Variogram mengukur ketidaksamaan rata-rata kuadrat pada jarak \(h\): \[ \gamma(h) = \tfrac{1}{2}\operatorname{Var}[Z(s) - Z(s+h)]. \] Untuk proses stasioner orde-2, relasi dengan kovarians: \[ \gamma(h) = C(0) - C(h). \]
6.9.2.4 Komponen Variogram
- Nugget \(c_0\): disontinuitas di \(h\to 0\) (noise/pengukuran, variasi mikro).
- Sill \(c_0 + c\): batas nilai variogram pada jarak jauh.
- Range \(a\): jarak efektif korelasi memudar.
6.9.2.5 Model Variogram Teoretis
- Eksponensial: \(\gamma(h)=c_0+c\big[1-\exp(-h/a)\big]\).
- Sferis: \[\gamma(h)= \begin{cases} c_0 + c\!\left[\frac{3h}{2a} - \frac{1}{2}\!\left(\frac{h}{a}\right)^3\right], & h\le a,\\\\ c_0 + c, & h>a. \end{cases} \]
- Gaussian: \(\gamma(h)=c_0+c\big[1-\exp(-(h/a)^2)\big]\).
- (Opsional) Matérn: fleksibel; meliputi Gaussian/Eksponensial sebagai kasus batas.
6.9.2.6 Detail Variogram
Dokumen ini berfokus pada konsep variogram dan semivariogram dalam geostatistik. Tujuan utamanya:
- Menjelaskan definisi dan interpretasi semivariogram.
- Memperkenalkan parameter penting: nugget, sill, dan range.
- Menunjukkan contoh model variogram teoritis (khususnya model sferis).
- Menampilkan grafik semivariogram lengkap dengan keterangan nugget, sill, dan range.
- Menunjukkan contoh empirical variogram dari data simulasi dan fitting model teoritis sederhana.
6.9.2.7 Konsep Dasar Variogram dan Semivariogram
6.9.2.7.1 Proses spasial dan autokorelasi
Misalkan \(Z(s)\) adalah peubah acak yang terdefinisi pada lokasi spasial \(s\) di suatu domain dua dimensi (misalnya koordinat \(x, y\)). Dalam geostatistik, kita tertarik pada bagaimana nilai \(Z(s)\) pada suatu lokasi berkorelasi dengan nilai di lokasi lain yang berjarak \(h\).
Secara umum, fungsi kovarian didefinisikan sebagai:
\[ C(h) = \text{Cov}(Z(s), Z(s + h)). \]
Untuk proses spasial yang second-order stationary, kovarian hanya bergantung pada jarak \(h\), bukan pada posisi absolut \(s\).
6.9.2.7.2 Definisi semivariogram
Semivariogram (sering disebut variogram) didefinisikan sebagai:
\[ \gamma(h) = \frac{1}{2} \operatorname{Var}\big( Z(s) - Z(s + h) \big). \]
Untuk proses yang stationer, terdapat hubungan antara semivariogram dan kovarian:
\[ \gamma(h) = C(0) - C(h), \]
di mana \(C(0) = \operatorname{Var}(Z(s))\) adalah variansi pada lag nol.
Interpretasi:
- Jika dua lokasi dekat satu sama lain (\(h\) kecil), kita harapkan \(Z(s)\) dan \(Z(s+h)\) mirip, sehingga variance dari selisihnya kecil dan \(\gamma(h)\) juga kecil.
- Semakin jauh jaraknya, nilai menjadi kurang mirip, sehingga semivariansi meningkat sampai mencapai suatu batas.
6.9.2.7.3 Parameter Penting: Nugget, Sill, dan Range
6.9.2.8 Nugget
Nugget adalah nilai semivariogram saat \(h \to 0\). Dalam praktik, sering ada nugget effect berupa lompatan kecil di dekat nol.
Penyebab umum:
- Kesalahan pengukuran (measurement error).
- Variasi mikro pada skala lebih kecil daripada jarak pengukuran.
Secara grafik, nugget adalah intersep semivariogram di dekat \(h = 0\).
6.9.2.9 Sill
Sill adalah nilai semivariogram ketika kurva sudah mencapai plateau dan menjadi relatif stabil.
Secara teoritis:
\[ \text{sill} \approx C(0) = \operatorname{Var}(Z(s)). \]
Setelah mencapai sill, penambahan jarak tidak lagi membuat semivariansi meningkat secara signifikan. Hal ini berarti titik-titik yang berjarak lebih jauh dari range tidak lagi berkorelasi secara spasial.
6.9.2.10 Range
Range adalah jarak (lag) di mana semivariogram mencapai sill (atau mendekati sill). Interpretasi praktis:
Titik data yang berjarak kurang dari range masih memiliki korelasi spasial, sedangkan yang lebih jauh dari range dianggap tidak saling berkorelasi secara spasial.
6.9.2.10.1 Model Variogram Sferis (Spherical)
Salah satu model variogram teoritis yang paling sering digunakan adalah model sferis (spherical model).
Misalkan:
- Nugget = \(c_0\),
- Sill total (nugget + struktur) = \(c_0 + c_1\),
- Range efektif = \(a\).
Bentuk fungsi semivariogram sferis:
\[ \gamma(h) = \begin{cases} c_0 + c_1\left[1.5\left(\dfrac{h}{a}\right) - 0.5\left(\dfrac{h}{a}\right)^3\right], & 0 \le h \le a, \\[4pt] c_0 + c_1, & h > a. \end{cases} \]
Di sini:
- Nugget = \(c_0\),
- Sill = \(c_0 + c_1\),
- Range = \(a\).
6.9.2.10.2 Contoh Grafik Semivariogram Teoritis dengan Nugget, Sill, dan Range
Pada bagian ini, kita membuat plot semivariogram sferis dengan parameter tertentu dan menandai nugget, sill, dan range.
# Parameter contoh
nugget <- 3 # c0
c1 <- 22 # kontribusi struktur
sill <- nugget + c1 # total sill
range <- 15 # a
# Fungsi variogram sferis
gamma_spherical <- function(h, nugget, c1, range) {
gamma <- ifelse(
h <= range,
nugget + c1 * (1.5 * (h / range) - 0.5 * (h / range)^3),
nugget + c1
)
gamma
}
h_vals <- seq(0, 30, length.out = 300)
gamma_vals <- gamma_spherical(h_vals, nugget, c1, range)
# Plot dasar
plot(h_vals, gamma_vals, type = "l", lwd = 2,
xlab = "Lag distance (h)",
ylab = expression(gamma(h)),
main = "Semivariogram Sferis dengan Nugget, Sill, dan Range")
# Tambah garis sill
abline(h = sill, lty = 2)
text(x = max(h_vals) * 0.6, y = sill + 0.8, labels = "Sill", adj = 0)
# Tambah garis nugget
abline(h = nugget, lty = 3)
text(x = max(h_vals) * 0.05, y = nugget + 0.8, labels = "Nugget", adj = 0)
# Tambah garis range
abline(v = range, lty = 2)
text(x = range + 0.5, y = max(gamma_vals) * 0.5, labels = "Range", adj = 0)
grid()
Grafik di atas menggambarkan:
- Kurva semivariogram yang naik dari nugget hingga sill.
- Nugget sebagai lompatan awal di dekat \(h = 0\).
- Sill sebagai plateau ketika semivarians sudah stabil.
- Range sebagai jarak ketika kurva mencapai sill.
6.9.2.10.3 Contoh Semivariogram Empirik dari Data Simulasi
Pada bagian ini, kita mensimulasikan data spasial sederhana di grid
2D, lalu menghitung empirical variogram menggunakan
paket gstat dan membandingkannya dengan model sferis
teoritis.
library(sp)
library(gstat)
set.seed(123)
# Buat grid 10x10
x <- 1:10
y <- 1:10
grid <- expand.grid(x = x, y = y)
coordinates(grid) <- ~ x + y
# Definisikan model variogram 'true' (sferis)
vgm_true <- vgm(psill = c1, model = "Sph", range = range, nugget = nugget)
# Simulasi satu realisasi random field Gaussian
g <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE,
beta = 0, model = vgm_true, nmax = 50)
sim_data <- predict(g, newdata = grid, nsim = 1)## [using unconditional Gaussian simulation]names(sim_data)## [1] "sim1"# Hitung semivariogram empirik
emp_vario <- variogram(sim1 ~ 1, data = sim_data)
# Tampilkan beberapa lag pertama
head(emp_vario)## np dist gamma dir.hor dir.ver id
## 1 180 1.000000 4.557552 0 0 var1
## 2 162 1.414214 4.678094 0 0 var1
## 3 448 2.151758 6.065766 0 0 var1
## 4 268 2.918055 7.323752 0 0 var1
## 5 252 3.162278 7.247562 0 0 var1
## 6 224 3.605551 7.349904 0 0 var1Sekarang kita plot semivariogram empirik dan overlay model sferis teoritis sebagai referensi.
# Empirical variogram
plot(emp_vario$dist, emp_vario$gamma,
xlab = "Lag distance (h)",
ylab = expression(gamma(h)),
main = "Semivariogram Empirik vs Model Sferis",
pch = 16,
ylim = c(0, max(c(emp_vario$gamma, nugget + c1)) * 0.5)
)
# Gamma teoritis dihitung tepat di lag yang sama dengan empirical variogram
gamma_theo <- gamma_spherical(emp_vario$dist, nugget, c1, range)
# Urutkan berdasarkan jarak supaya garisnya halus
ord <- order(emp_vario$dist)
lines(emp_vario$dist[ord], gamma_theo[ord],
lwd = 2)
legend("bottomright",
legend = c("Empirik", "Model sferis teoritis"),
pch = c(16, NA),
lty = c(NA, 1),
lwd = c(NA, 2)
)
grid()
vgm_fit <- vgm(psill = c1, model = "Sph", range = range, nugget = nugget)
plot(emp_vario, model = vgm_fit)
Dari grafik di atas:
- Titik-titik hitam adalah semivariogram empirik yang dihitung dari data simulasi.
- Garis halus adalah model sferis teoritis dengan parameter yang digunakan saat simulasi.
- Dalam praktik, proses fitting variogram melakukan hal sebaliknya: kita punya semivariogram empirik, lalu menyesuaikan parameter (nugget, sill, range) supaya kurva teoritis mendekati pola empirik.
Ringkasan
- Semivariogram mengukur seberapa cepat kemiripan spasial menurun dengan jarak.
- Parameter kunci:
- Nugget: semivariansi pada lag mendekati nol, mencerminkan error pengukuran dan variasi mikro.
- Sill: nilai plateau semivariogram; mendekati variansi keseluruhan proses.
- Range: jarak di mana semivariogram mencapai sill; di luar range, korelasi spasial diabaikan.
- Model variogram teoritis (sferis, eksponensial, Gaussian, dsb.) digunakan dalam kriging dan metode prediksi spasial lainnya.
- Semivariogram empirik dihitung dari data, kemudian digunakan untuk menaksir parameter model teoritis.
Dokumen ini dapat dikembangkan lebih lanjut dengan:
- Menambah pembahasan variogram anisotropik,
- Membandingkan beberapa bentuk model (eksponensial, Gaussian, Matérn),
- Menghubungkan langsung dengan implementasi kriging (ordinary / universal kriging).
6.9.3 Keluarga Kriging: SK, OK, UK
Tuliskan \(Z(s) = \mu(s) + \varepsilon(s)\) dengan \(E[\varepsilon(s)]=0\).
- Simple Kriging (SK): \(\mu(s)\equiv\mu\) diketahui (jarang realistis).
- Ordinary Kriging (OK): \(\mu(s)\equiv\mu\) tak diketahui tetapi konstan lokal.
- Universal Kriging (UK): \(\mu(s)=\mathbf{x}(s)^\top\beta\) (mis. polinomial koordinat).
6.9.4 Sistem Persamaan Kriging & Variansi Prediksi
Untuk OK, syarat tak bias: \(\sum_i \lambda_i = 1\), diberlakukan via pengali Lagrange \(\mu\): \[ \begin{bmatrix} \Sigma & \mathbf{1} \\\\ \mathbf{1}^\top & 0 \end{bmatrix} \begin{bmatrix} \lambda \\\\ \mu \end{bmatrix} = \begin{bmatrix} \sigma_0 \\\\ 1 \end{bmatrix}, \quad \sigma_K^2(s_0) = C(0) - 2\lambda^\top \sigma_0 + \lambda^\top \Sigma \lambda. \] Di mana \(\Sigma_{ij}=C(s_i-s_j)\), \(\sigma_{0,i}=C(s_i-s_0)\). Dalam bentuk variogram gunakan \(C(0)-C(h)=\gamma(h)\).
Bagian ini memberikan contoh perhitungan manual Ordinary
Kriging dengan dataset kecil (3 titik), lalu
memverifikasi hasilnya menggunakan paket
gstat di R.
Data & Target
Tiga titik pengamatan \(Z\) dan satu lokasi target \(s_0=(1,1)\):
| Titik | x | y | Z |
|---|---|---|---|
| A | 0 | 0 | 10 |
| B | 2 | 0 | 20 |
| C | 0 | 2 | 30 |
Target: \(s_0 = (1,1)\).
6.9.5 Model Variogram
Gunakan spherical variogram dengan parameter: - Nugget \(c_0=0\) - Sill parsial \(c=1\) - Range \(a=2\)
Bentuk fungsi variogram: \[ \gamma(h) = \begin{cases} c_0 + c\left[\tfrac{3h}{2a} - \tfrac{1}{2}\left(\tfrac{h}{a}\right)^3\right], & h \le a, \\ c_0 + c, & h > a. \end{cases} \]
Perhitungan Manual
1) Jarak antar titik & ke target
A <- c(0,0); B <- c(2,0); C <- c(0,2); s0 <- c(1,1)
edist <- function(p,q) sqrt(sum((p-q)^2))
h_AB <- edist(A,B)
h_AC <- edist(A,C)
h_BC <- edist(B,C)
h_As0 <- edist(A,s0)
h_Bs0 <- edist(B,s0)
h_Cs0 <- edist(C,s0)
c(h_AB=h_AB, h_AC=h_AC, h_BC=h_BC, h_As0=h_As0, h_Bs0=h_Bs0, h_Cs0=h_Cs0)## h_AB h_AC h_BC h_As0 h_Bs0 h_Cs0
## 2.000000 2.000000 2.828427 1.414214 1.414214 1.414214Hasil: \(h_{AB}=h_{AC}=2\); \(h_{BC}=2.828\); \(h_{A,s_0}=h_{B,s_0}=h_{C,s_0}=1.414\).
2) Fungsi variogram \(\gamma(h)\)
# Spherical variogram with c0=0, c=1, a=2
vgm_sph <- function(h, c0=0, c=1, a=2) {
out <- ifelse(h <= a, c0 + c*(1.5*(h/a) - 0.5*(h/a)^3), c0 + c)
as.numeric(out)
}
c(g_2 = vgm_sph(2), g_2_828 = vgm_sph(sqrt(8)), g_1_414 = vgm_sph(sqrt(2)))## g_2 g_2_828 g_1_414
## 1.0000000 1.0000000 0.8838835Untuk \(h=1.414\): \(\gamma \approx 1.038\); untuk \(h\ge a\), \(\gamma=1\).
3) Sistem Persamaan Kriging (Ordinary)
Sistem OK untuk 3 titik:
\[ \begin{bmatrix} 0 & \gamma_{AB} & \gamma_{AC} & 1 \\ \gamma_{BA} & 0 & \gamma_{BC} & 1 \\ \gamma_{CA} & \gamma_{CB} & 0 & 1 \\ 1 & 1 & 1 & 0 \end{bmatrix} \begin{bmatrix} \lambda_A \\ \lambda_B \\ \lambda_C \\ \mu \end{bmatrix} = \begin{bmatrix} \gamma_{A,s_0} \\ \gamma_{B,s_0} \\ \gamma_{C,s_0} \\ 1 \end{bmatrix}. \]
Dengan \(\gamma_{AB}=\gamma_{AC}=1\), \(\gamma_{BC}=1\), dan \(\gamma_{A,s_0}=\gamma_{B,s_0}=\gamma_{C,s_0}=1.038\).
G <- matrix(c(
0, 1, 1, 1,
1, 0, 1, 1,
1, 1, 0, 1,
1, 1, 1, 0
), nrow=4, byrow=TRUE)
rhs <- c(1.038, 1.038, 1.038, 1)
sol <- solve(G, rhs)
sol## [1] 0.3333333 0.3333333 0.3333333 0.3713333Komponen solusi: \(\lambda_A, \lambda_B, \lambda_C, \mu\). Periksa penjumlahan bobot = 1.
lambda <- sol[1:3]; mu <- sol[4]
sum_lambda <- sum(lambda)
Z <- c(A=10, B=20, C=30)
Zhat <- sum(lambda * Z)
list(weights=lambda, sum_lambda=sum_lambda, Zhat=Zhat, mu=mu)## $weights
## [1] 0.3333333 0.3333333 0.3333333
##
## $sum_lambda
## [1] 1
##
## $Zhat
## [1] 20
##
## $mu
## [1] 0.3713333Hasil manual: \(\lambda_A=\lambda_B=\lambda_C=1/3\), \(\hat Z(s_0)=20\).
6.9.6 Varian Kriging
\[ \sigma_k^2 = \sum_i \lambda_i\, \gamma(s_i, s_0) + \mu. \]
g_As0 <- vgm_sph(sqrt(2))
var_krig <- sum(lambda * c(g_As0, g_As0, g_As0)) + mu
var_krig## [1] 1.255217Varian kriging \(\approx 1.41\).
Verifikasi dengan gstat
library(sp)
library(gstat)
# Observasi sebagai SpatialPointsDataFrame
df <- data.frame(x=c(0,2,0), y=c(0,0,2), Z=c(10,20,30))
coordinates(df) <- ~x+y
# Lokasi prediksi (s0)
newd <- data.frame(x=1, y=1)
coordinates(newd) <- ~x+y
# Variogram model: Spherical, psill=1, range=2, nugget=0
vgm_model <- vgm(psill=1, model="Sph", range=2, nugget=0)
# Ordinary kriging
ok <- gstat(formula = Z~1, data=df, model=vgm_model)
kr <- predict(ok, newd, nsim=0)## [using ordinary kriging]kr## class : SpatialPointsDataFrame
## features : 1
## extent : 1, 1, 1, 1 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 2
## names : var1.pred, var1.var
## value : 20, 1.1011002862997Kolom var1.pred adalah \(\hat
Z(s_0)\) dan var1.var adalah varian kriging. Nilai
seharusnya sangat dekat dengan hasil manual (20 dan ~1.41).
Peta Permukaan Prediksi & Varian
# Grid untuk peta
xg <- seq(-0.5, 2.5, by=0.05)
yg <- seq(-0.5, 2.5, by=0.05)
newgrid <- expand.grid(x=xg, y=yg)
coordinates(newgrid) <- ~x+y
krg <- predict(ok, newgrid)## [using ordinary kriging]kr_df <- as.data.frame(krg)
library(ggplot2)
p1 <- ggplot(kr_df, aes(x, y, fill = var1.pred)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = var1.pred), color = "white", alpha = 0.6) +
geom_point(data = as.data.frame(df), aes(x, y), inherit.aes = FALSE, size = 2) +
coord_equal() +
labs(title = "Ordinary Kriging — Prediksi", fill = "Ẑ") +
theme_minimal(base_size = 12)
p2 <- ggplot(kr_df, aes(x, y, fill = var1.var)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = var1.var), color = "white", alpha = 0.6) +
geom_point(data = as.data.frame(df), aes(x, y), inherit.aes = FALSE, size = 2) +
coord_equal() +
labs(title = "Ordinary Kriging — Varian Prediksi", fill = "Var") +
theme_minimal(base_size = 12)
p1
p2
Catatan Praktis
- Pemilihan model variogram krusial: gunakan empirical variogram + fitting (OLS/WLS/REML) pada praktik nyata.
- Nugget menangkap variasi mikro/skala-pengukuran; tidak identik dengan error observasi semata.
- OK vs UK: Ordinary Kriging (mean konstan tak diketahui) vs Universal Kriging (mean fungsi kovariat/koordinat). Untuk tren kuat, pertimbangkan Regression–Kriging.
6.9.7 Implementasi Dasar
(CRS-Safe) — sp + gstat
# Fungsi bantu: set CRS aman ke EPSG:28992 (Amersfoort / RD New)
set_crs <- function(obj){
ok <- TRUE
tryCatch({
proj4string(obj) <- CRS("+init=epsg:28992")
}, error = function(e) ok <<- FALSE)
if(!ok){
proj4string(obj) <- CRS(SRS_string = "EPSG:28992")
}
obj
}
# Data meuse
data(meuse, package = "sp")
coordinates(meuse) <- ~x+y
meuse <- set_crs(meuse)
# Grid prediksi
data(meuse.grid, package = "sp")
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS(proj4string(meuse)) # samakan CRS
# Cek ringkas
head(meuse@data[,c("zinc","copper")])## zinc copper
## 1 1022 85
## 2 1141 81
## 3 640 68
## 4 257 81
## 5 269 48
## 6 281 616.9.8 Variogram Empirik & Fitting Model
v.emp <- variogram(log(zinc) ~ 1, meuse)
plot(v.emp, main = "Variogram Empirik log(Zinc)")
# Fit model Spherical (tebakan awal; silakan eksplorasi model lain)
v.fit <- fit.variogram(v.emp, vgm(psill = 1, model = "Sph", range = 900, nugget = 1))
v.fit## model psill range
## 1 Nug 0.05066243 0.0000
## 2 Sph 0.59060780 897.0209plot(v.emp, v.fit, main = "Fitting Variogram: Spherical")
6.9.9 Ordinary Kriging (OK)
6.9.9.1 Detail Ordinary Kriging
Dokumen ini menjelaskan konsep dasar Ordinary Kriging
(OK), dengan contoh implementasi menggunakan data
meuse di R.
Secara singkat, ordinary kriging adalah metode interpolasi geostatistik yang:
- memandang fenomena spasial sebagai random field,
- menggunakan semivariogram untuk memodelkan korelasi spasial,
- menghasilkan prediksi linear tak-bias dengan varian minimum,
- menyediakan peta prediksi dan peta ketidakpastian (varian kriging).
6.9.9.2 Konsep Teoretis Ordinary Kriging
6.9.9.3 Random Field dan Stationarity
Misalkan \(Z(s)\) adalah nilai variabel di lokasi \(s = (x, y)\). Kita modelkan:
\[ Z(s) = \mu + e(s), \]
dengan:
- \(\mu\): rata-rata (mean) konstan pada neighborhood lokal (asumsi OK),
- \(e(s)\): komponen acak yang berkorelasi di ruang.
Ordinary kriging mengasumsikan bahwa:
- Rata-rata \(\mu\) konstan tapi tidak diketahui,
- Struktur korelasi diwakili oleh semivariogram \(\gamma(h)\).
6.9.9.4 Semivariogram
Semivariogram empiris pada jarak (lag) \(h\) didefinisikan sebagai:
\[ \hat\gamma(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} \big[ Z(s_i) - Z(s_i + h) \big]^2, \]
dengan \(N(h)\) adalah jumlah pasangan titik yang jaraknya sekitar \(h\).
Parameter penting:
- Nugget \(\gamma(0)\): semivarians saat jarak mendekati 0,
- Sill: level plateau semivariogram (kurang lebih sama dengan varians total),
- Range: jarak di mana semivariogram mencapai sill (di atas range, korelasi spasial \(\approx 0\)).
Berbagai bentuk model teoritis semivariogram yang digunakan antara lain:
- Spherical
- Exponential
- Gaussian
- Linear
Sebagai contoh, model spherical (tanpa nugget) dapat ditulis:
\[ \gamma(h) = \begin{cases} C_0 + C\left[ \frac{3h}{2a} - \frac{1}{2}\left(\frac{h}{a}\right)^3 \right], & 0 < h \le a \\ C_0 + C, & h > a \end{cases} \]
dengan \(C_0\) = nugget, \(C\) = partial sill, \(a\) = range.
6.9.9.5 Estimator Ordinary Kriging
Misalkan kita punya \(n\) observasi \(Z(s_1), \dots, Z(s_n)\). Ordinary kriging mencari estimator linear:
\[ \hat Z(s_0) = \sum_{i=1}^n \lambda_i Z(s_i), \]
dengan dua syarat:
Tak-bias: \[ \mathbb{E}[\hat Z(s_0) - Z(s_0)] = 0 \quad \Rightarrow \quad \sum_{i=1}^n \lambda_i = 1. \]
Varian galat minimum: \[ \text{Var}[\hat Z(s_0) - Z(s_0)] \to \min. \]
Dengan menggunakan Lagrange multiplier \(\mu\), kita peroleh sistem persamaan ordinary kriging:
\[ \begin{cases} \sum_{j=1}^n \lambda_j \gamma(s_i, s_j) + \mu = \gamma(s_i, s_0), & i=1,\dots,n, \\[4pt] \sum_{j=1}^n \lambda_j = 1. \end{cases} \]
Dalam bentuk matriks:
\[ \begin{bmatrix} \gamma_{11} & \dots & \gamma_{1n} & 1 \\ \vdots & \ddots & \vdots & \vdots \\ \gamma_{n1} & \dots & \gamma_{nn} & 1 \\ 1 & \dots & 1 & 0 \end{bmatrix} \begin{bmatrix} \lambda_1 \\ \vdots \\ \lambda_n \\ \mu \end{bmatrix} = \begin{bmatrix} \gamma(s_1, s_0) \\ \vdots \\ \gamma(s_n, s_0) \\ 1 \end{bmatrix}. \]
Setelah bobot \(\lambda_i\) diperoleh, prediksi di \(s_0\) adalah:
\[ \hat Z(s_0) = \sum_{i=1}^n \lambda_i Z(s_i), \]
sedangkan varian kriging (ketidakpastian) bisa dihitung sebagai:
\[ \sigma_K^2(s_0) = \sum_{i=1}^n \lambda_i \gamma(s_i, s_0) + \mu. \]
6.9.9.6 Manual
Dokumen ini menunjukkan perhitungan manual ordinary kriging untuk contoh kecil, langkah demi langkah, lalu diverifikasi menggunakan R.
Fokus kita:
- Menyiapkan data titik (koordinat dan nilai).
- Menghitung jarak antar titik dan ke titik prediksi.
- Menentukan model variogram teoritis (dengan parameter yang dianggap diketahui).
- Menghitung kovarian dari variogram.
- Menyusun dan menyelesaikan sistem persamaan ordinary kriging (dengan Lagrange multiplier).
- Menghitung nilai prediksi dan ragam kriging.
- Mengecek hasil menggunakan perhitungan matriks di R dan paket
gstat.
Catatan: Di dunia nyata, parameter variogram (nugget, sill, range) biasanya diestimasi dari data (fit variogram). Untuk tujuan latihan manual, di sini kita pilih parameter secara sengaja supaya perhitungannya relatif rapi dan simetris.
6.9.9.7 Data Contoh
Kita gunakan 4 titik data yang membentuk bujur sangkar, dengan titik prediksi di tengah.
- Titik data:
- P1: \((0, 0)\), nilai \(Z_1 = 10\)
- P2: \((0, 10)\), nilai \(Z_2 = 20\)
- P3: \((10, 0)\), nilai \(Z_3 = 30\)
- P4: \((10, 10)\), nilai \(Z_4 = 40\)
- Titik prediksi:
- P0: \((5, 5)\)
library(tibble)
data_pts <- tribble(
~id, ~x, ~y, ~z,
"P1", 0, 0, 10,
"P2", 0, 10, 20,
"P3", 10, 0, 30,
"P4", 10, 10, 40
)
pred_pt <- tibble(
id = "P0",
x = 5,
y = 5
)
data_pts## # A tibble: 4 × 4
## id x y z
## <chr> <dbl> <dbl> <dbl>
## 1 P1 0 0 10
## 2 P2 0 10 20
## 3 P3 10 0 30
## 4 P4 10 10 40pred_pt## # A tibble: 1 × 3
## id x y
## <chr> <dbl> <dbl>
## 1 P0 5 5Secara geometri:
- Keempat titik data berjarak 10 satuan secara horizontal/vertikal.
- Jarak diagonal antar titik data adalah \(\sqrt{10^2 + 10^2} = \sqrt{200} \approx 14.142\).
- Jarak dari titik prediksi \((5,5)\) ke keempat titik data sama, yaitu \[ d = \sqrt{5^2 + 5^2} = \sqrt{50} \approx 7.071. \]
6.9.9.8 Langkah 1 – Menghitung Jarak
6.9.9.8.1 Jarak Antar Titik Data
Kita hitung jarak Euclidean:
\[ d_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} \]
dist_matrix <- as.matrix(dist(data_pts[, c("x", "y")], method = "euclidean"))
rownames(dist_matrix) <- data_pts$id
colnames(dist_matrix) <- data_pts$id
dist_matrix## P1 P2 P3 P4
## P1 0.00000 10.00000 10.00000 14.14214
## P2 10.00000 0.00000 14.14214 10.00000
## P3 10.00000 14.14214 0.00000 10.00000
## P4 14.14214 10.00000 10.00000 0.00000Interpretasi (hasil di atas):
- \(d_{P1,P2} = 10\)
- \(d_{P1,P3} = 10\)
- \(d_{P1,P4} = 14.142...\)
- \(d_{P2,P3} = 14.142...\)
- \(d_{P2,P4} = 10\)
- \(d_{P3,P4} = 10\)
6.9.9.8.2 Jarak Titik Data ke Titik Prediksi
d_pred <- sqrt((data_pts$x - pred_pt$x)^2 + (data_pts$y - pred_pt$y)^2)
d_pred <- setNames(d_pred, data_pts$id)
d_pred## P1 P2 P3 P4
## 7.071068 7.071068 7.071068 7.071068Semua jarak data–prediksi sama: \[ d_{P1,P0} = d_{P2,P0} = d_{P3,P0} = d_{P4,P0} \approx 7.071. \]
6.9.9.9 Langkah 2 – Model Variogram Teoritis
Kita pilih model variogram sferis:
\[ \gamma(h) = \begin{cases} C_0\left[1.5\left(\dfrac{h}{a}\right) - 0.5\left(\dfrac{h}{a}\right)^3\right], & 0 \le h \le a, \\ C_0, & h > a. \end{cases} \]
Parameter yang kita pilih:
- Nugget = 0 (untuk menyederhanakan).
- Sill struktural \(C_0 = 25\).
- Range \(a = 15\).
Sehingga:
- Semivariogram: \(\gamma(h)\) seperti di atas.
- Kovarian: \[ C(h) = C(0) - \gamma(h), \quad C(0) = C_0 = 25. \]
Kita implementasikan fungsi \(\gamma(h)\) dan \(C(h)\) di R.
C0 <- 25 # sill
a <- 15 # range
gamma_spherical <- function(h, C0, a) {
h <- abs(h)
out <- ifelse(
h <= a,
C0 * (1.5 * (h / a) - 0.5 * (h / a)^3),
C0
)
out
}
cov_from_gamma <- function(h, C0, a) {
C0 - gamma_spherical(h, C0, a)
}
# Contoh nilai gamma dan kovarian untuk jarak penting
h_values <- c(0, 7.07106781187, 10, 14.1421356237)
gamma_values <- gamma_spherical(h_values, C0, a)
cov_values <- cov_from_gamma(h_values, C0, a)
cbind(h = h_values, gamma = gamma_values, cov = cov_values)## h gamma cov
## [1,] 0.000000 0.00000 25.000000
## [2,] 7.071068 16.36821 8.631787
## [3,] 10.000000 21.29630 3.703704
## [4,] 14.142136 24.87968 0.120317Dari tabel di atas (pembulatan):
- Untuk \(h = 0\):
- \(\gamma(0) = 0\)
- \(C(0) = 25\)
- Untuk \(h \approx 7.071\) (jarak
data–prediksi):
- \(\gamma(7.071) \approx 16.3682\)
- \(C(7.071) \approx 8.6318\)
- Untuk \(h = 10\) (jarak sisi bujur
sangkar):
- \(\gamma(10) \approx 21.2963\)
- \(C(10) \approx 3.7037\)
- Untuk \(h \approx 14.142\) (jarak
diagonal):
- \(\gamma(14.142) \approx 24.8797\)
- \(C(14.142) \approx 0.1203\)
6.9.9.10 Langkah 3 – Matriks Kovarian Antar Titik Data
Matriks kovarian \(4 \times 4\) dengan elemen:
\[ C_{ij} = C(h_{ij}) = C(\text{jarak antara titik } i \text{ dan } j). \]
# Fungsi untuk membangun matriks kovarian dari matriks jarak
cov_matrix <- cov_from_gamma(dist_matrix, C0 = C0, a = a)
cov_matrix## P1 P2 P3 P4
## P1 25.000000 3.703704 3.703704 0.120317
## P2 3.703704 25.000000 0.120317 3.703704
## P3 3.703704 0.120317 25.000000 3.703704
## P4 0.120317 3.703704 3.703704 25.000000Penjelasan:
- Diagonal: \(C(0) = 25\)
- Elemen dengan jarak 10: sekitar \(3.7037\)
- Elemen dengan jarak 14.142: sekitar \(0.1203\)
6.9.9.11 Langkah 4 – Vektor Kovarian ke Titik Prediksi
Vektor \(c\) berisi kovarian antara setiap titik data dengan titik prediksi:
\[ c_i = C(\text{jarak titik } i \text{ ke P0}). \]
Karena semua jaraknya sama (\(\approx 7.071\)), semua elemen vektor sama:
cov_pred <- cov_from_gamma(d_pred, C0 = C0, a = a)
cov_pred## P1 P2 P3 P4
## 8.631787 8.631787 8.631787 8.6317876.9.9.12 Langkah 5 – Sistem Persamaan Ordinary Kriging
Untuk ordinary kriging dengan \(n = 4\) titik, bobot \(w_1, w_2, w_3, w_4\) dan Lagrange multiplier \(\mu\) memenuhi:
\[ \begin{bmatrix} C_{11} & C_{12} & C_{13} & C_{14} & 1 \\ C_{21} & C_{22} & C_{23} & C_{24} & 1 \\ C_{31} & C_{32} & C_{33} & C_{34} & 1 \\ C_{41} & C_{42} & C_{43} & C_{44} & 1 \\ 1 & 1 & 1 & 1 & 0 \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \\ w_3 \\ w_4 \\ \mu \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \\ 1 \end{bmatrix} \]
Kita bentuk matriks dan selesaikan secara numerik di R (ini menghitung manual dalam arti memakai aljabar linear, bukan fungsi kriging “instan”).
# Matriks K (4x4 kovarian + kolom ones, baris ones di bawah)
K <- rbind(
cbind(cov_matrix, 1),
c(1, 1, 1, 1, 0)
)
# Vektor sisi kanan
rhs <- c(cov_pred, 1)
K## P1 P2 P3 P4
## P1 25.000000 3.703704 3.703704 0.120317 1
## P2 3.703704 25.000000 0.120317 3.703704 1
## P3 3.703704 0.120317 25.000000 3.703704 1
## P4 0.120317 3.703704 3.703704 25.000000 1
## 1.000000 1.000000 1.000000 1.000000 0rhs## P1 P2 P3 P4
## 8.631787 8.631787 8.631787 8.631787 1.000000Sekarang kita selesaikan sistem linear:
\[ K \cdot \begin{bmatrix}w \\ \mu\end{bmatrix} = \text{rhs} \]
sol <- solve(K, rhs)
sol## P1 P2 P3 P4
## 0.2500000 0.2500000 0.2500000 0.2500000 0.4998564Hasil sol berisi:
sol[1:4]= bobot \(w_1, w_2, w_3, w_4\)sol[5]= Lagrange multiplier \(\mu\)
6.9.9.12.1 Interpretasi Simetris (Analitik)
Karena konfigurasi titik simetris terhadap titik prediksi:
- Jarak dari P1, P2, P3, dan P4 ke P0 sama (\(\approx 7.071\)).
- Struktur kovarian antar titik juga simetris (misalnya titik-titik “posisi” yang sepadan).
Sangat kuat dugaan bahwa:
\[ w_1 = w_2 = w_3 = w_4 = w. \]
Dari constraint ordinary kriging:
\[ w_1 + w_2 + w_3 + w_4 = 1 \quad \Rightarrow \quad 4w = 1 \Rightarrow w = 0.25. \]
Kita cek apakah ini konsisten dengan persamaan kriging.
Untuk titik P1, misalnya:
\[ C_{11} w_1 + C_{12} w_2 + C_{13} w_3 + C_{14} w_4 + \mu = c_1. \]
Karena semua \(w_i = w\):
\[ (C_{11} + C_{12} + C_{13} + C_{14}) w + \mu = c_1. \]
Dengan:
- \(C_{11} = 25\)
- \(C_{12} = C_{13} = C(10) \approx 3.7037\)
- \(C_{14} = C(14.142) \approx 0.1203\)
- \(c_1 = C(7.071) \approx 8.6318\)
Total kovarian:
\[ C_{11} + C_{12} + C_{13} + C_{14} \approx 25 + 3.7037 + 3.7037 + 0.1203 \approx 32.5277. \]
Masukkan \(w = 0.25\):
\[ 32.5277 \times 0.25 + \mu = 8.6318 \]
\[ 8.1319 + \mu = 8.6318 \Rightarrow \mu \approx 0.4999 \approx 0.5. \]
Jadi solusi analitik:
- \(w_1 = w_2 = w_3 = w_4 = 0.25\)
- \(\mu \approx 0.5\)
Silakan bandingkan dengan hasil solve(K, rhs).
6.9.9.13 Langkah 6 – Prediksi Ordinary Kriging
Prediksi ordinary kriging:
\[ \hat{Z}(P0) = \sum_{i=1}^4 w_i Z_i. \]
Dengan \(w_i = 0.25\) dan nilai \(Z_i\):
- \(Z_1 = 10\)
- \(Z_2 = 20\)
- \(Z_3 = 30\)
- \(Z_4 = 40\)
\[ \hat{Z}(P0) = 0.25(10 + 20 + 30 + 40) = 0.25 \times 100 = 25. \]
Kita hitung di R untuk verifikasi.
w <- sol[1:4]
z <- data_pts$z
Z_hat <- sum(w * z)
Z_hat## [1] 25Secara manual dan numerik, hasil prediksi sekitar 25 (sangat rapi karena desain data simetris).
6.9.9.14 Langkah 7 – Ragam Kriging
Ragam kriging untuk ordinary kriging:
\[ \sigma_K^2 = C(0) - \sum_{i=1}^4 w_i c_i - \mu. \]
Dengan:
- \(C(0) = 25\)
- \(c_i = C(7.071) \approx 8.6318\) (semua sama)
- \(\sum w_i = 1\)
- \(\mu \approx 0.5\)
Maka:
\[ \sum_{i=1}^4 w_i c_i = \left(\sum_{i=1}^4 w_i\right) c_i = 1 \times 8.6318 = 8.6318. \]
\[ \sigma_K^2 = 25 - 8.6318 - 0.5 \approx 15.8682. \]
Kita hitung di R:
mu <- sol[5]
C0_val <- C0
c_dot_w <- sum(w * cov_pred)
sigma2_K <- C0_val - c_dot_w - mu
sigma2_K##
## 15.86836Jadi ragam kriging di titik P0 sekitar 15.87.
6.9.9.15 Langkah 8 –
Verifikasi dengan gstat::krige
Sebagai cek tambahan, kita gunakan paket gstat dan
memasukkan variogram model dengan parameter yang sama.
library(sp)
library(gstat)
# 1. Buat objek sp
coordinates(data_pts) <- ~ x + y
coordinates(pred_pt) <- ~ x + y
# 2. Definisikan model variogram:
# variogram model di gstat: psill = C0 (sill struktural), range = a, nugget = 0
vgm_model <- vgm(psill = C0, model = "Sph", range = a, nugget = 0)
vgm_model## model psill range
## 1 Nug 0 0
## 2 Sph 25 15Lakukan ordinary kriging:
ok_result <- krige(
formula = z ~ 1,
locations = data_pts,
newdata = pred_pt,
model = vgm_model
)## [using ordinary kriging]ok_result## class : SpatialPointsDataFrame
## features : 1
## extent : 5, 5, 5, 5 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 2
## names : var1.pred, var1.var
## value : 25, 15.8683561463474Biasanya output var1.pred dan var1.var akan
sedikit berbeda karena gstat menghitung variogram dengan
konvensi dan pembulatan tertentu. Namun nilainya harus sangat
dekat dengan hasil manual:
- Prediksi sekitar 25
- Ragam sekitar 15–16
6.9.9.16 Ringkasan Alur Perhitungan Manual
Data & titik prediksi
- Tentukan koordinat dan nilai \(Z_i\).
- Tentukan lokasi titik prediksi.
Hitung jarak
- Jarak antar titik data (\(h_{ij}\)).
- Jarak titik data ke titik prediksi (\(h_{i0}\)).
Pilih/estimasi model variogram
- Di sini: model sferis, nugget = 0, sill = 25, range = 15.
- Hitung \(\gamma(h)\) dan kovarian \(C(h) = C(0) - \gamma(h)\).
Bangun matriks kovarian & vektor kovarian
- Matriks \(C_{ij} = C(h_{ij})\).
- Vektor \(c_i = C(h_{i0})\).
Bangun sistem ordinary kriging
\[ \begin{bmatrix} C & 1 \\ 1^\top & 0 \end{bmatrix} \begin{bmatrix} w \\ \mu \end{bmatrix} = \begin{bmatrix} c \\ 1 \end{bmatrix} \]Selesaikan sistem
- Dapatkan bobot \(w_i\) dan \(\mu\) (dengan aljabar linear).
- Cek bahwa \(\sum w_i = 1\).
Hitung prediksi
\[ \hat{Z} = \sum_i w_i Z_i. \]Hitung ragam kriging
\[ \sigma_K^2 = C(0) - \sum_i w_i c_i - \mu. \]
Contoh ini sengaja dibuat simetris sehingga bobotnya sama semua (0.25). Untuk konfigurasi yang tidak simetris, langkah aljabar sama, hanya angka-angkanya saja yang lebih “liar”.
6.9.9.17 Langkah Analisis di R dengan Data Meuse
Pada bagian ini, kita melakukan ordinary kriging untuk kadar Zn pada
data meuse yang disediakan oleh paket
sp/gstat.
6.9.9.18 Persiapan Paket dan Data
# Install dulu jika belum pernah:
# install.packages(c("sf", "gstat", "ggplot2", "dplyr", "viridis"))
library(sf)
library(gstat)
library(ggplot2)
library(dplyr)
library(viridis)
# Memuat data meuse
data(meuse, package = "sp")
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 470# Konversi ke sf (CRS RD_New / EPSG:28992)
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
# Plot lokasi sampling Zn
ggplot(meuse_sf) +
geom_sf(aes(color = zinc)) +
scale_color_viridis() +
theme_minimal() +
labs(title = "Lokasi Sampling dan Kadar Zn (Data Meuse)",
color = "Zn (ppm)")
6.9.9.19 Variogram Empiris dan Model Teoritis
# Membuat objek gstat
gs <- gstat(id = "zinc", formula = zinc ~ 1, data = meuse_sf)
# Variogram empiris
vgm_emp <- variogram(gs)
vgm_emp## np dist gamma dir.hor dir.ver id
## 1 57 79.29244 37362.96 0 0 zinc
## 2 299 163.97367 72718.34 0 0 zinc
## 3 419 267.36483 82655.53 0 0 zinc
## 4 457 372.73542 111575.91 0 0 zinc
## 5 547 478.47670 123080.69 0 0 zinc
## 6 533 585.34058 147414.28 0 0 zinc
## 7 574 693.14526 142891.51 0 0 zinc
## 8 564 796.18365 153563.90 0 0 zinc
## 9 589 903.14650 160217.69 0 0 zinc
## 10 543 1011.29177 168867.36 0 0 zinc
## 11 500 1117.86235 186528.20 0 0 zinc
## 12 477 1221.32810 150320.76 0 0 zinc
## 13 452 1329.16407 180399.66 0 0 zinc
## 14 457 1437.25620 143139.84 0 0 zinc
## 15 415 1543.20248 144112.31 0 0 zincplot(vgm_emp, main = "Semivariogram Empiris Zn (Meuse)")
# Fit model variogram teoritis (misal spherical)
# Nilai awal (psill, range, nugget) bisa disesuaikan
vgm_start <- vgm(psill = 1e5, model = "Sph", range = 500, nugget = 0)
vgm_mod <- fit.variogram(vgm_emp, model = vgm_start)
vgm_mod## model psill range
## 1 Nug 24797.73 0.0000
## 2 Sph 134742.93 830.8785plot(vgm_emp, vgm_mod, main = "Semivariogram Empiris dan Model Spherical")
Tabel parameter vgm_mod akan berisi perkiraan nugget,
sill, dan range.
6.9.9.20 Pembuatan Grid Prediksi
# Membuat grid prediksi berdasarkan bounding box meuse
bbox <- st_bbox(meuse_sf)
res <- 50 # resolusi grid (meter)
x_seq <- seq(bbox["xmin"], bbox["xmax"], by = res)
y_seq <- seq(bbox["ymin"], bbox["ymax"], by = res)
grid_df <- expand.grid(x = x_seq, y = y_seq)
grid_sf <- st_as_sf(grid_df, coords = c("x", "y"), crs = st_crs(meuse_sf))
grid_sf[1:5, ]## Simple feature collection with 5 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 178605 ymin: 329714 xmax: 178805 ymax: 329714
## Projected CRS: Amersfoort / RD New
## geometry
## 1 POINT (178605 329714)
## 2 POINT (178655 329714)
## 3 POINT (178705 329714)
## 4 POINT (178755 329714)
## 5 POINT (178805 329714)6.9.9.21 Ordinary Kriging dengan R
# Ordinary kriging Zn
ok_result <- krige(
formula = zinc ~ 1,
locations = meuse_sf,
newdata = grid_sf,
model = vgm_mod
)## [using ordinary kriging]head(ok_result)## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 178605 ymin: 329714 xmax: 178855 ymax: 329714
## Projected CRS: Amersfoort / RD New
## var1.pred var1.var geometry
## 1 701.0351 114737.60 POINT (178605 329714)
## 2 708.6196 103485.53 POINT (178655 329714)
## 3 715.8440 92315.99 POINT (178705 329714)
## 4 718.6294 83048.55 POINT (178755 329714)
## 5 710.6198 77659.33 POINT (178805 329714)
## 6 690.8862 75130.79 POINT (178855 329714)Kolom penting:
var1.pred: prediksi Zn,var1.var: varian kriging (ketidakpastian).
6.9.9.22 Visualisasi Hasil Prediksi
ggplot(ok_result) +
geom_sf(aes(color = var1.pred)) +
scale_color_viridis() +
theme_minimal() +
labs(
title = "Prediksi Kadar Zn dengan Ordinary Kriging",
color = "Zn pred."
)
ggplot(ok_result) +
geom_sf(aes(color = var1.var)) +
scale_color_viridis() +
theme_minimal() +
labs(
title = "Varian Kriging (Ketidakpastian Prediksi Zn)",
color = "Var"
)
Dari peta varian kriging, biasanya kita melihat bahwa ketidakpastian lebih besar di area yang jauh dari titik pengamatan.
6.9.9.23 Cross-Validation (Validasi Model Kriging)
Kita dapat melakukan leave-one-out cross-validation untuk menilai seberapa baik model variogram dan kriging.
set.seed(123)
cv_ok <- krige.cv(
formula = zinc ~ 1,
locations = meuse_sf,
model = vgm_mod,
nfold = nrow(meuse_sf) # LOO
)
head(cv_ok)## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 181025 ymin: 333260 xmax: 181390 ymax: 333611
## Projected CRS: Amersfoort / RD New
## var1.pred var1.var observed residual zscore fold
## 1 890.6908 61747.46 1022 131.3091763 0.5284276463 1
## 2 868.8885 59922.24 1141 272.1115024 1.1116111476 2
## 3 661.2496 60456.28 640 -21.2495816 -0.0864230646 3
## 4 465.7753 72131.35 257 -208.7752608 -0.7773506972 4
## 5 268.8929 58749.38 269 0.1070912 0.0004418267 5
## 6 235.0607 73389.74 281 45.9392713 0.1695767933 6
## geometry
## 1 POINT (181072 333611)
## 2 POINT (181025 333558)
## 3 POINT (181165 333537)
## 4 POINT (181298 333484)
## 5 POINT (181307 333330)
## 6 POINT (181390 333260)# Menghitung RMSE
rmse_ok <- sqrt(mean((cv_ok$observed - cv_ok$var1.pred)^2))
rmse_ok## [1] 224.7866ggplot(as.data.frame(cv_ok), aes(x = observed, y = var1.pred)) +
geom_point(alpha = 0.7) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
theme_minimal() +
labs(
title = paste0("Cross-validation Ordinary Kriging (RMSE = ", round(rmse_ok, 2), ")"),
x = "Observasi Zn",
y = "Prediksi Zn (LOO)"
)
6.9.9.24 Ringkasan
- Ordinary kriging mengasumsikan mean lokal konstan dan tidak diketahui.
- Struktur korelasi spasial dimodelkan melalui semivariogram.
- Bobot kriging diperoleh dengan menyelesaikan sistem persamaan linier yang menjamin prediksi tak-bias dan ber-varian minimum.
- R menyediakan fungsi praktis untuk:
- menghitung variogram (
variogram), - mem-fit model (
fit.variogram), - melakukan kriging (
krige), - dan validasi (
krige.cv).
- menghitung variogram (
Dokumen ini dapat dikembangkan lebih lanjut dengan:
- membandingkan OK dengan IDW pada data yang sama,
- mempertimbangkan anisotropi,
- mencoba transformasi (log, Box–Cox) jika distribusi sangat skewed,
- serta menggabungkan kriging dalam MCDA (Multi-Criteria Decision Analysis).
library(sp)
library(gstat)
data(meuse)
coordinates(meuse) <- ~x+y
data(meuse.grid)
coordinates(meuse.grid) <- ~x+y # ← key step
gridded(meuse.grid) <- TRUE # ← then mark as grid
v <- variogram(log(zinc) ~ 1, meuse)
v.fit <- fit.variogram(v, model = vgm(1, "Sph", 300, 1))
ok <- krige(log(zinc) ~ 1, meuse, meuse.grid, model = v.fit)## [using ordinary kriging]spplot(ok["var1.pred"], main = "Prediksi log(Zinc) — Ordinary Kriging")
spplot(ok["var1.var"], main = "Kriging Variance — Ordinary Kriging")
6.9.9.25 Validasi Silang (10-fold)
set.seed(1)
cv <- krige.cv(log(zinc) ~ 1, meuse, model = v.fit, nfold = 10)
summary(cv$residual)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.927351 -0.209664 -0.007911 0.013350 0.212601 1.373383ggplot(cv@data, aes(x = observed, y = var1.pred)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
labs(x = "Observed log(Zinc)", y = "Predicted", title = "10-fold Cross-Validation (OK)")
6.9.10 Universal Kriging (UK)
Jika ada tren global, mis. \(\mu(s)=\beta_0+\beta_1 x + \beta_2 y\).
uk <- krige(log(zinc) ~ x + y, meuse, meuse.grid, model = v.fit)## [using universal kriging]spplot(uk["var1.pred"], main = "Universal Kriging (x + y) — Prediksi log(Zinc)")
Tips praktis: tambahkan polinomial atau kovariat lain (jarak ke sungai, elevasi) bila relevan.
6.9.10.1 Diagnostik & Praktik Baik
- Atur lebar bin variogram dengan bijak; terlalu halus → berisik, terlalu kasar → pola hilang.
- Cek anisotropi via variogram per sektor arah.
- Lihat peta residu (OK/UK) untuk tren sisa.
- Gunakan transformasi (mis. log) bila distribusi sangat skewed (seperti zinc).
6.9.11 Co-Kriging (Multivariat) — Linear Model of Coregionalization (LMC)
6.9.11.1 Konsep
Kita punya peubah target \(Z_1\) (zinc) dan peubah sekunder \(Z_2\) (copper) yang berkorelasi spasial. Bangun LMC yang memodelkan variogram masing-masing dan cross-variogram.
6.9.11.2 Implementasi
# Definisikan dua proses
gs <- gstat(NULL, id = "zinc", formula = log(zinc) ~ 1, data = meuse)
gs <- gstat(gs, id = "copper", formula = log(copper) ~ 1, data = meuse)
# Variogram & cross-variogram
v.cross <- variogram(gs)
plot(v.cross, main = "Variogram & Cross-Variogram (zinc, copper)")
# Fit LMC (seed Spherical)
vm <- vgm(psill = 1, model = "Sph", range = 900, nugget = 1)
lmc.fit <- fit.lmc(v.cross, gs, model = vm)
plot(v.cross, lmc.fit, main = "Fitting LMC (Spherical seed)")
# Prediksi co-kriging pada grid
ck.pred <- predict(lmc.fit, meuse.grid) # kolom: zinc.pred, copper.pred, dst.## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]spplot(ck.pred["zinc.pred"], main = "Prediksi Zinc — Co-Kriging")
Catatan: Co-Kriging efektif bila peubah sekunder informatif (korelasi spasial & cross variogram jelas), serta tersedia dengan cakupan spasial baik.
6.9.12 Kriging Ruang–Waktu
Untuk \(Z(s,t)\), gunakan \(C(h,u)\) dengan lag ruang \(h\) dan waktu \(u\). Model separabel: \(C(h,u)=C_s(h)C_t(u)\). Implementasi dapat
memakai gstat (variogram spatio-temporal) atau pendekatan
Bayesian (INLA) dengan efek acak ruang–waktu.
6.9.13 Kriging Bayesian via INLA–SPDE (Opsional)
6.9.13.1 Intuisi SPDE
SPDE mengaitkan GRF kontinu dan GMRF pada mesh triangulasi melalui \[ (\kappa^2 - \Delta)^{\alpha/2} x(s) = W(s). \] Keuntungan: komputasi efisien (matriks presisi sparse) dan inferensi probabilistik penuh.
6.9.13.2 Contoh Kode
(dibiarkan eval=FALSE)
suppressPackageStartupMessages(library(INLA))
# Mesh
mesh <- inla.mesh.2d(loc = coordinates(meuse), max.edge = c(200, 400))
plot(mesh); points(meuse, col = "red", pch = 19)
# SPDE Matérn dengan prior PC
spde <- inla.spde2.pcmatern(
mesh, alpha = 2,
prior.range = c(1000, 0.5), # P(range < 1000) = 0.5
prior.sigma = c(1, 0.01) # P(sigma > 1) = 0.01
)
# Index & proyeksi
iset <- inla.spde.make.index("spatial", n.spde = spde$n.spde)
A <- inla.spde.make.A(mesh, loc = coordinates(meuse))
# Stack
stk <- inla.stack(
data = list(y = log(meuse$zinc)),
A = list(A, 1),
effects = list(iset, data.frame(Intercept = 1)),
tag = "est"
)
# Fitting
formula <- y ~ 1 + f(spatial, model = spde)
res <- inla(
formula,
data = inla.stack.data(stk),
family = "gaussian",
control.predictor = list(A = inla.stack.A(stk), compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
summary(res)
# Posterior mean pada mesh
mean.field <- res$summary.random$spatial$mean
plot(mesh, main = "Posterior Mean Field (SPDE)"); points(meuse, pch = 20, cex = .5)6.9.14 Rangkuman & Latihan Singkat
Rangkuman: Variogram adalah jantung Kriging; OK/UK menangani mean lokal vs tren; Co-Kriging memanfaatkan peubah sekunder; SPDE memberi skala besar dan ketidakpastian penuh.
Latihan:
- Bandingkan RMSE CV antara OK vs UK (x+y). Manakah yang lebih baik
untuk data
meuse? - Ubah kernel variogram (Eksponensial/Gaussian/Matérn) dan evaluasi perubahannya.
- Tambahkan peubah sekunder lain (mis.
lead,cadmium) untuk Co-Kriging; bandingkan peta.
6.9.15 Ringkasan Kelebihan dan Keterbatasan Metode Interpolasi
Berikut perbandingan ringkas beberapa metode interpolasi spasial yang umum digunakan, termasuk Kriging sebagai pendekatan berbasis model stokastik yang paling kuat secara statistik.
6.9.16 Nearest Neighbor (NN)
Kelebihan:
- Super cepat dan mudah diimplementasikan.
- Cocok untuk baseline awal atau data dengan sebaran sangat
jarang.
Keterbatasan:
- Hasilnya “bloky” (diskontinu di batas Voronoi).
- Tidak menangkap transisi spasial yang halus.
6.9.17 TIN (Triangulated Irregular Network / Linear Interpolation)
Kelebihan:
- Menghasilkan permukaan halus di dalam setiap segitiga Delaunay.
- Konservatif: tidak memprediksi di luar convex hull data.
Keterbatasan:
- Sensitif terhadap kualitas dan distribusi triangulasi.
- Sulit diterapkan jika data sangat tidak merata.
6.9.18 Inverse Distance Weighting (IDW)
Kelebihan:
- Intuitif, parameter sederhana (pangkat jarak \(p\)).
- Memberikan hasil halus tanpa perlu model variogram.
Keterbatasan:
- Tidak menyediakan varian prediksi (ketidakpastian).
- Sensitif terhadap nilai ekstrem di sekitar lokasi target.
- Tidak mempertimbangkan struktur autokorelasi spasial.
6.9.19 Trend Surface (Polynomial Regression)
Kelebihan:
- Mudah dan cepat untuk menangkap pola global (misal gradien
timur–barat).
- Cocok untuk fenomena dengan variasi deterministik besar.
Keterbatasan:
- Mengabaikan autokorelasi residu.
- Tidak cocok untuk pola lokal yang kompleks.
6.9.20 Thin Plate Splines (TPS)
Kelebihan:
- Sangat halus dan fleksibel; parameter halus (\(\lambda\)) dapat dipilih otomatis lewat
GCV.
- Cocok untuk fenomena kontinu dan permukaan lembut.
Keterbatasan:
- Tidak menyediakan peta ketidakpastian.
- Dapat over-smooth jika \(\lambda\) tidak diatur dengan tepat.
- Komputasi cukup berat untuk dataset besar.
6.9.21 Kriging
Kelebihan:
- Penduga linier tak bias terbaik (BLUE).
- Mempertimbangkan struktur autokorelasi melalui
variogram/kovarians.
- Menghasilkan peta prediksi dan peta ketidakpastian
(varian) secara bersamaan.
- Dapat diperluas menjadi Co-Kriging, Universal Kriging, maupun Kriging
spasio-temporal.
Keterbatasan:
- Membutuhkan pemodelan variogram yang tepat dan validasi model.
- Asumsi stasioneritas sering sulit terpenuhi.
- Komputasi intensif untuk data besar (meskipun dapat diatasi dengan
pendekatan Bayesian/INLA-SPDE).
6.10 Kesimpulan Umum
|:-|:|:|:–| | NN | Non-parametrik | Cepat, sederhana | Bloky, tanpa transisi halus | | TIN | Geometrik | Halus dalam segitiga | Sensitif terhadap triangulasi | | IDW | Deterministik | Intuitif, parameter sederhana | Tidak ada varian prediksi | | Trend | Regresi | Menangkap pola global | Abaikan korelasi spasial | | TPS | Smoothing | Halus, fleksibel | Tidak sediakan varian prediksi | | Kriging | Stokastik | Optimal, sediakan ketidakpastian | Butuh variogram, berat komputasi |