simulasi ESF spasial temporal Pada Pendugaan Curah Hujan

library(spdep)

## Warning: package 'spdep' was built under R version 4.5.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.5.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Warning: package 'sf' was built under R version 4.5.3

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

# 1. Setup Spasial & Eigenvector
coords <- matrix(runif(20), ncol=2)
W <- as.matrix(dist(coords))
W_bin <- ifelse(W > 0 & W < 0.5, 1, 0)
ev <- eigen(W_bin)$vectors[, 1:2] # Ambil 2 filter utama

# 2. Setup Parameter
n <- 10; TT <- 10; N <- n * TT
E_long <- ev[rep(1:n, TT), ] # Replikasi filter untuk data panel
b0_hat <- NULL; b1_hat <- NULL

# 3. Looping Simulasi 100 Kali
for (i in 1:100) {
  time_vec <- rep(1:TT, each=n)
  X <- rnorm(N, 70, 5)
  eps <- rnorm(N, 0, 1)
  
  # Koefisien bervariasi (STVC)
  beta1_st <- 0.5 + (0.02 * time_vec) + (0.1 * E_long[,1])
  
  # Pembangkitan Y (Curah Hujan)
  Y <- 15 + (X * beta1_st) + (E_long[,1] * 2) + eps
  
  # Estimasi dengan OLS (sebagai pembanding dasar)
  model <- lm(Y ~ X + E_long)
  b0_hat <- c(b0_hat, coef(model)[1])
  b1_hat <- c(b1_hat, coef(model)[2])
}

# 4. Hasil Simulasi
hasil <- matrix(c(mean(b0_hat), sd(b0_hat), mean(b1_hat), sd(b1_hat)), nrow=2)
rownames(hasil) <- c("mean", "sd")
colnames(hasil) <- c("beta0", "beta1")
hasil

##          beta0     beta1
## mean 16.537351 0.5896852
## sd    6.775253 0.0964613

Hasil simulasi : Mean: Menunjukkan nilai rata-rata estimasi. Jika mendekati parameter awal (15 dan 0.5), maka model konsisten.SD: Menunjukkan tingkat variasi atau ketelitian model dalam menduga curah hujan di berbagai titik waktu dan ruang.

B. Visualisasi Peta Spasial (Spatial Plot) untuk melihat di mana curah hujan paling tinggi pada titik waktu tertentu.

library(ggplot2)
library(patchwork)

## Warning: package 'patchwork' was built under R version 4.5.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Menggabungkan koordinat, waktu, dan hasil prediksi (Y) ke dalam satu dataframe
df_vis <- data.frame(
  lon = rep(coords[,1], TT),
  lat = rep(coords[,2], TT),
  waktu = rep(1:TT, each = n),
  curah_hujan = Y, # Hasil Y dari simulasi terakhir
  prediksi = fitted(model)
)

plot_spasial <- ggplot(df_vis %>% filter(waktu == 1), aes(x = lon, y = lat)) +
  geom_point(aes(color = curah_hujan, size = curah_hujan)) +
  scale_color_viridis_c(option = "mako") +
  labs(title = "Sebaran Curah Hujan (Waktu = 1)",
       subtitle = "Model ESF-STVC",
       x = "Bujur", y = "Lintang", color = "CH (mm)") +
  theme_minimal()

print(plot_spasial)

C. Visualisasi Tren Temporal (Time Series Plot) Ini menunjukkan bagaimana model menangkap variasi curah hujan dari waktu ke waktu untuk satu lokasi tertentu.

plot_temporal <- ggplot(df_vis %>% filter(lon == coords[1,1]), aes(x = waktu)) +
  geom_line(aes(y = curah_hujan, color = "Observasi"), size = 1) +
  geom_line(aes(y = prediksi, color = "Prediksi ESF"), linetype = "dashed", size = 1) +
  labs(title = "Tren Curah Hujan di Lokasi A",
       x = "Periode Waktu", y = "Curah Hujan (mm)",
       color = "Keterangan") +
  theme_light()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(plot_temporal)

D. Visualisasi Error/Sisaan (Penting untuk Model Spasial) Untuk melihat apakah Spatial Filter (Eigenvector) berhasil menghilangkan autokorelasi, kita bisa melihat plot Residual vs Fitted.

plot_residual <- ggplot(data.frame(res = resid(model), fit = fitted(model)), aes(x = fit, y = res)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, color = "red") +
  labs(title = "Analisis Sisaan Model ESF",
       x = "Nilai Prediksi", y = "Residual") +
  theme_classic()
plot_residual