#================ 2.1 Input Data Stunting (sesuai file) ================#
library(readxl) 
library(dplyr)

## 
## 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

library(stringr)
library(ggplot2) 
library(forcats) 
library(sf)

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

Stunting <- read.csv("Stunting2.csv", 
                   header = TRUE, sep = ";", dec = ".")
head(Stunting)

##    No.  ID Provinsi           Kab.Kota Tahun Prev.Stunting
## 1 1029 448     ACEH      Kab. Simeulue  2021          18.9
## 2 1030   6     ACEH  Kab. Aceh Singkil  2021           0.0
## 3 1031   5     ACEH  Kab. Aceh Selatan  2021           0.0
## 4 1032   9     ACEH Kab. Aceh Tenggara  2021          18.9
## 5 1033  10     ACEH    Kab. Aceh Timur  2021          22.6
## 6 1034   8     ACEH   Kab. Aceh Tengah  2021           8.8
##   Tingkat.Pendidikan.Ibu  PDRB Tingkat.Kemiskinan      LAT      LON
## 1                     40  2460              18.98 2.617357 96.08060
## 2                  38.77  2714              20.36 2.349671 97.84902
## 3                   32.6  5970              13.18 3.164681 97.43539
## 4                  43.27  5422              13.41 3.371797 97.69460
## 5                  28.85 11698              14.45 4.628417 97.62920
## 6                  45.16  8032              15.26 4.530942 96.85939

#======== Choose Data: 2020 (selaraskan dengan contoh) ========#
Stunting20 <- Stunting[Stunting$Tahun == 2020, ]
head(Stunting20)

##  [1] No.                    ID                     Provinsi              
##  [4] Kab.Kota               Tahun                  Prev.Stunting         
##  [7] Tingkat.Pendidikan.Ibu PDRB                   Tingkat.Kemiskinan    
## [10] LAT                    LON                   
## <0 rows> (or 0-length row.names)

#======== Data Preprocessing (gaya contoh) ========#
# (angka kita sudah titik; tetap siapkan guard jika ada koma)
to_num <- function(x){
  if (is.numeric(x)) return(x)
  x |> as.character() |> str_replace(",", ".") |> as.numeric()
}

norm_wilayah <- function(x){
  x |>
    tolower() |>
    str_replace_all("\\s+", " ") |>
    str_trim() |>
    str_remove("^kab\\.?\\s+") |>
    str_remove("^kota\\s+adm(?:inistrasi)?\\s+") |>
    str_remove("^kota\\s+") |>
    str_replace_all("dki jakarta|daerah khusus ibu kota jakarta", "jakarta") |>
    str_replace_all("kepulauan seribu|adm\\. kepulauan seribu", "adm. kepulauan seribu")
}

Stunting20 <- Stunting20 |> transmute( provinsi = Provinsi, kabupaten_kota = kab/kota, tahun = Tahun, stunting = to_num(Prev.Stunting), pendidikan_ibu = to_num(tingkat.pendidikan.ibu), pdrb = to_num(pdrb), kemiskinan = to_num(tingkat.kemiskinan), lat = to_num(lAT), lon = to_num(LON), prov_clean = norm_wilayah(Provinsi), kab_clean = norm_wilayah(Kab/Kota) ) head(Stunting20)

#================ 2.1.1 Exploratory Data Analysis ================#
## 2.1.1.1 Statistics Descriptive (per provinsi)
st_stats_by_prov <- Stunting20 %>%
  group_by(Provinsi) %>%
  summarise(
    count      = n(),
    mean_st    = mean(Prev.Stunting, na.rm=TRUE),
    median_st  = median(Prev.Stunting, na.rm=TRUE),
    sd_st      = sd(Prev.Stunting, na.rm=TRUE),
    min_st     = min(Prev.Stunting, na.rm=TRUE),
    max_st     = max(Prev.Stunting, na.rm=TRUE)
  )

## Warning: There were 2 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `min_st = min(Prev.Stunting, na.rm = TRUE)`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

st_stats_by_prov

## # A tibble: 0 × 7
## # ℹ 7 variables: Provinsi <chr>, count <int>, mean_st <dbl>, median_st <dbl>,
## #   sd_st <dbl>, min_st <dbl>, max_st <dbl>

## Descriptive statistics: overall
mean(Stunting20$Prev.Stunting, na.rm=TRUE)

## [1] NaN

sd(Stunting20$Prev.Stunting,   na.rm=TRUE)

## [1] NA

min(Stunting20$Prev.Stunting,  na.rm=TRUE)

## Warning in min(Stunting20$Prev.Stunting, na.rm = TRUE): no non-missing
## arguments to min; returning Inf

## [1] Inf

max(Stunting20$Prev.Stunting,  na.rm=TRUE)

## Warning in max(Stunting20$Prev.Stunting, na.rm = TRUE): no non-missing
## arguments to max; returning -Inf

## [1] -Inf

#================ 2.1.1.2 Boxplot (Outlier Detection) ================#
ggplot(
  Stunting20 |> mutate(provinsi = fct_reorder(Provinsi, Prev.Stunting, median, .desc=TRUE)),
  aes(x = Provinsi, y = Prev.Stunting, fill = Provinsi)
) +
  geom_boxplot(color = "black") +
  labs(
    title = paste0("Boxplot of Stunting Prevalence by Province (", 2020, ")"),
    x = "Province", y = "Prevalence (%)"
  ) +
  guides(fill = "none") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1))

#================ 2.2 Input Indonesia Maps (RDS) ================#
indo_sf <- readRDS("map_sf_ringan")

#================ 2.3 Join Data dengan Semua Tahun ================#
indo_join_all <- indo_sf %>%
  left_join(
    Stunting %>% filter(Tahun %in% 2021:2024),
    by = c("KAB_KOTA" = "Kab.Kota")
  )

#================ 2.4 Visualisasi (Kategori 2021–2024) ================#
# 1) Skala kategori & palet
breaks <- c(0,10,20,30,40,Inf)
labels <- c("<10%","10–20%","20–30%","30–40%","≥40%")
pal_merah <- c("#FFFACD", "#FFD700", "#FFA500", "#FF4500", "#B22222")

# 2) Pastikan nilai numeric
indo_join_all <- indo_join_all %>%
  dplyr::mutate(
    .st_raw = as.character(`Prev.Stunting`),
    .st_raw = gsub(",", ".", .st_raw),
    stunting_num = as.numeric(.st_raw),
    stunting_cat = cut(stunting_num, breaks = breaks, labels = labels, include.lowest = TRUE)
  )

# 3) Plot facet (2021–2024)
ggplot2::ggplot(indo_join_all %>% filter(Tahun %in% 2021:2024)) +
  ggplot2::geom_sf(aes(fill = stunting_cat), color = "grey75", linewidth = 0.15) +
  ggplot2::scale_fill_manual(values = pal_merah, na.value = "grey85", drop = FALSE) +
  ggplot2::labs(
    title = "Sebaran Prevalensi Stunting Kabupaten/Kota di Indonesia (2021–2024)",
    fill  = "Kategori",
    caption = "NA = tidak ada data atau tidak match nama"
  ) +
  ggplot2::facet_wrap(~Tahun, ncol = 2) +
  ggplot2::theme_bw(base_size = 11) +
  ggplot2::theme(
    strip.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 13)
  )

# Buat kombinasi kabupaten × tahun
locations <- unique(Stunting[, c("LON","LAT","Kab.Kota","Provinsi")])
all_years <- data.frame(Tahun = 2021:2024)

all_combinations <- merge(locations, all_years)
final_stunting <- merge(all_combinations, Stunting, 
                        by = c("LON","LAT","Kab.Kota","Provinsi","Tahun"), all.x = TRUE)

# Pastikan Prev.Stunting numeric
final_stunting$Prev.Stunting <- as.numeric(final_stunting$Prev.Stunting)

# Peta facet per tahun
ggplot() +
  geom_sf(data = indo_sf, fill = "grey95", color = "white") +
  geom_point(data = final_stunting,
             aes(x = LON, y = LAT, color = Prev.Stunting), size = 1.8) +
  facet_wrap(~Tahun, ncol = 2) +
  scale_color_gradient(low = "yellow", high = "red", name = "Prevalensi (%)") +
  labs(title = "Distribusi Prevalensi Stunting per Tahun (2021–2024)",
       x = "Longitude", y = "Latitude") +
  theme_bw()

#================== 1. Siapkan interval kategori ==================#
breaks <- c(0, 10, 20, 30, 40, Inf)
labels <- c("<10%", "10–20%", "20–30%", "30–40%", "≥40%")

# Warna kuning → oranye → merah (seperti peta gempa)
pal_kuning_merah <- c("yellow", "#FECC5C", "#FD8D3C", "#F03B20", "#BD0026")

#================== 2. Mutasi data ke numerik & kategori ==================#
final_stunting <- final_stunting |>
  dplyr::mutate(
    .st_raw = as.character(`Prev.Stunting`),
    .st_raw = gsub(",", ".", .st_raw),
    Prev.Stunting = as.numeric(.st_raw),
    stunting_cat = cut(
      Prev.Stunting,
      breaks = breaks,
      labels = labels,
      include.lowest = TRUE
    )
  )

#================== 3. Plot dengan interval kategori ==================#
ggplot() +
  geom_sf(data = indo_sf, fill = "grey95", color = "white") +
  geom_point(
    data = final_stunting,
    aes(x = LON, y = LAT, color = stunting_cat),
    size = 1.8
  ) +
  facet_wrap(~Tahun, ncol = 2) +
  scale_color_manual(
    values = pal_kuning_merah,
    name = "Kategori Prevalensi",
    drop = FALSE,
    na.value = "grey80"
  ) +
  labs(
    title = "Distribusi Prevalensi Stunting per Tahun (2021–2024)",
    subtitle = "Warna menunjukkan interval kategori prevalensi stunting",
    x = "Longitude", y = "Latitude"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    strip.text = element_text(face = "bold"),
    legend.position = "right",
    panel.grid.major = element_line(color = "grey90", linetype = "dotted")
  )

indo_union <- indo_sf |> st_union() |> st_make_valid()
indo_utm <- st_transform(indo_union, 32748)  # zona 48S cukup representatif

# --- Buat grid 10 km ---
cellsize <- 10000
grid10 <- st_make_grid(indo_utm, cellsize = cellsize, what = "centers")

# --- Titik di dalam poligon ---
grid_inside <- grid10[lengths(st_within(grid10, indo_utm)) > 0]
grid_inside_sf <- st_sf(geometry = grid_inside)

# --- Plot dengan geom_sf (tidak ada garis silang!) ---
ggplot() +
  geom_sf(data = indo_utm, fill = NA, color = "red", linewidth = 0.4) +
  geom_sf(data = grid_inside_sf, color = "blue", size = 0.3) +
  labs(
    title = "Interval Grid 10 km x 10 km di dalam batas Indonesia",
    subtitle = "Grid titik dalam proyeksi UTM zone 48S",
    x = "Easting (m)", y = "Northing (m)"
  ) +
  theme_minimal()

make_grid_plot <- function(cellsize) {
  grid <- st_make_grid(indo_utm, cellsize = cellsize, what = "centers")
  inside <- grid[lengths(st_within(grid, indo_utm)) > 0]
  ggplot() +
    geom_sf(data = indo_utm, fill = NA, color = "red", linewidth = 0.3) +
    geom_sf(data = st_sf(geometry = inside), color = "blue", size = 0.3) +
    labs(title = paste0("Interval Grid: ", cellsize/1000, " km x ", cellsize/1000, " km")) +
    theme_minimal()
}

make_grid_plot(10000)

make_grid_plot(20000)

make_grid_plot(50000)

# ================== helper pembuat grid ==================
grid_points_inside <- function(boundary_utm, cell_km){
  cs <- cell_km * 1000
  pts <- st_make_grid(boundary_utm, cellsize = cs, what = "centers", square = TRUE)
  inside <- pts[lengths(st_within(pts, boundary_utm)) > 0]
  st_sf(grid_id = seq_along(inside), geometry = inside)
}

grid_polys_inside <- function(boundary_utm, cell_km){
  cs <- cell_km * 1000
  polys <- st_make_grid(boundary_utm, cellsize = cs, what = "polygons", square = TRUE)
  polys_sf <- st_sf(grid_id = seq_along(polys), geometry = polys)
  st_intersection(polys_sf, boundary_utm) |> 
    mutate(grid_id = as.integer(grid_id))
}

# ================== helper plotting ==================
plot_grid_pts_utm <- function(boundary_utm, grid_pts_sf, title_txt){
  ggplot() +
    geom_sf(data = boundary_utm, fill = NA, color = "red", linewidth = 0.4) +
    geom_sf(data = grid_pts_sf, color = "blue", size = 0.3) +
    labs(title = title_txt, x = "Easting (m)", y = "Northing (m)") +
    theme_minimal()
}

plot_grid_pts_lonlat <- function(boundary_utm, grid_pts_sf, title_txt){
  ggplot() +
    geom_sf(data = st_transform(boundary_utm, 4326), fill = NA, color = "red", linewidth = 0.4) +
    geom_sf(data = st_transform(grid_pts_sf, 4326), color = "blue", size = 0.3) +
    labs(title = title_txt, x = "Longitude", y = "Latitude") +
    theme_minimal()
}

# ================== contoh penggunaan (titik pusat sel) ==================
g10  <- grid_points_inside(indo_utm, 10)
g20  <- grid_points_inside(indo_utm, 20)
g50  <- grid_points_inside(indo_utm, 50)

# tampilkan dalam UTM (meter)
plot_grid_pts_utm(indo_utm, g10, "Interval Grid: 10 km x 10 km (UTM)")

plot_grid_pts_utm(indo_utm, g20, "Interval Grid: 20 km x 20 km (UTM)")

plot_grid_pts_utm(indo_utm, g50, "Interval Grid: 50 km x 50 km (UTM)")

# atau tampilkan dalam lon/lat (lebih enak buat laporan)
plot_grid_pts_lonlat(indo_utm, g10, "Interval Grid: 10 km x 10 km (Lon/Lat)")

plot_grid_pts_lonlat(indo_utm, g20, "Interval Grid: 20 km x 20 km (Lon/Lat)")

plot_grid_pts_lonlat(indo_utm, g50, "Interval Grid: 50 km x 50 km (Lon/Lat)")

###5. Prepare the data in the UTM system

st_crs(indo_sf)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

# ===============================================
# KONVERSI KOORDINAT & TAMBAH KOL. UTM_x, UTM_y, IT (TAHUN)
# ===============================================

library(sf)
library(dplyr)

# --- Data stunting nasional (514 kab/kota)
# Sudah punya kolom LAT, LON, Prev.Stunting, dan Tahun
# contoh: head(Stunting)

# 1️⃣ Pastikan nama kolom konsisten
Stunting <- Stunting |>
  rename(
    kabupaten = `Kab.Kota`,
    provinsi  = Provinsi,
    tahun     = Tahun,
    prev_stunting = `Prev.Stunting`,
    lon = LON,
    lat = LAT
  )

# 2️⃣ Ubah ke objek sf dengan CRS WGS84 (EPSG:4326)
Stunting_sf <- st_as_sf(
  Stunting,
  coords = c("lon", "lat"),
  crs = 4326,
  remove = FALSE
)

# 3️⃣ Transform ke UTM (zona 48S – meter)
Stunting_sf_utm <- st_transform(Stunting_sf, crs = 32748)

# 4️⃣ Ambil koordinat UTM dan tambahkan ke data
coords <- st_coordinates(Stunting_sf_utm)
Stunting_sf_utm$UTM_x <- coords[, 1]
Stunting_sf_utm$UTM_y <- coords[, 2]

# 5️⃣ Tambahkan kolom waktu tahunan (IT)
Stunting_sf_utm$IT <- Stunting_sf_utm$tahun

# 6️⃣ Cek hasil
head(Stunting_sf_utm[, c("provinsi", "kabupaten", "tahun", "prev_stunting", "UTM_x", "UTM_y", "IT")])

## Simple feature collection with 6 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -495516.6 ymin: 10261760 xmax: -297139.8 ymax: 10515860
## Projected CRS: WGS 84 / UTM zone 48S
##   provinsi          kabupaten tahun prev_stunting     UTM_x    UTM_y   IT
## 1     ACEH      Kab. Simeulue  2021          18.9 -495516.6 10292859 2021
## 2     ACEH  Kab. Aceh Singkil  2021           0.0 -297139.8 10261758 2021
## 3     ACEH  Kab. Aceh Selatan  2021           0.0 -342931.0 10352882 2021
## 4     ACEH Kab. Aceh Tenggara  2021          18.9 -313718.1 10375753 2021
## 5     ACEH    Kab. Aceh Timur  2021          22.6 -319778.6 10515863 2021
## 6     ACEH   Kab. Aceh Tengah  2021           8.8 -406070.2 10505927 2021
##                     geometry
## 1 POINT (-495516.6 10292859)
## 2 POINT (-297139.8 10261758)
## 3   POINT (-342931 10352882)
## 4 POINT (-313718.1 10375753)
## 5 POINT (-319778.6 10515863)
## 6 POINT (-406070.2 10505927)

# ========================================================
# MESH SPASIAL UNTUK KASUS STUNTING DI INDONESIA
# ========================================================
library(INLA)

## Loading required package: Matrix

## This is INLA_25.06.07 built 2025-06-11 19:15:03 UTC.
##  - See www.r-inla.org/contact-us for how to get help.
##  - List available models/likelihoods/etc with inla.list.models()
##  - Use inla.doc(<NAME>) to access documentation
##  - Consider upgrading R-INLA to testing[25.10.28] or stable[25.10.19].

library(sf)
library(ggplot2)

# --- Data input:
# Stunting_sf_utm : data stunting dengan kolom UTM_x, UTM_y (satuan meter)
# indo_sf         : shapefile batas Indonesia

# --- 1. Pastikan CRS sama (UTM 48S)
indo_utm <- st_transform(indo_sf, 32748)

# --- 2. Buat mesh triangulasi untuk INLA
stunting_mesh <- inla.mesh.2d(
  loc = cbind(Stunting_sf_utm$UTM_x, Stunting_sf_utm$UTM_y),
  offset = c(50000, 200000),     # margin tambahan luar dalam
  max.edge = c(100000, 1000000), # panjang sisi min–maks (meter)
  cutoff = 20000                 # jarak minimum antar titik (hindari duplikasi)
)

# --- 3. Plot versi dasar
plot(stunting_mesh, col = "gray80",
     main = "INLA Spatial Mesh – Prevalensi Stunting Indonesia (UTM 48S)")
points(Stunting_sf_utm$UTM_x, Stunting_sf_utm$UTM_y, col = "red", pch = 16, cex = 0.4)
plot(st_geometry(indo_utm), add = TRUE, border = "blue")

# --- 4. Plot versi ggplot (revisi tanpa error)
triangles <- stunting_mesh$graph$tv     # indeks segitiga dari mesh
mesh_coords <- as.data.frame(stunting_mesh$loc)
colnames(mesh_coords) <- c("x", "y")    # tambahkan nama kolom

# Buat data frame garis sisi segitiga (utk ggplot)
mesh_edges <- do.call(rbind, lapply(1:nrow(triangles), function(i) {
  tri <- triangles[i, ]
  data.frame(
    x = mesh_coords[tri[c(1, 2, 3, 1)], 1],
    y = mesh_coords[tri[c(1, 2, 3, 1)], 2],
    group = i
  )
}))

ggplot() +
  geom_path(data = mesh_edges, aes(x = x, y = y, group = group),
            color = "gray70", linewidth = 0.2) +
  geom_sf(data = indo_utm, fill = NA, color = "blue", linewidth = 0.3) +
  geom_point(data = as.data.frame(Stunting_sf_utm),
             aes(x = UTM_x, y = UTM_y),
             color = "red", size = 0.3) +
  labs(
    title = "Mesh Spasial INLA untuk Model Prevalensi Stunting Indonesia",
    subtitle = "Jaringan segitiga (mesh) untuk pendekatan SPDE",
    x = "Easting (m)", y = "Northing (m)"
  ) +
  theme_minimal()

# --- lokasi observasi (semua tahun)
coords_all <- as.matrix(data.frame(
  x = Stunting_sf_utm$UTM_x,
  y = Stunting_sf_utm$UTM_y
))

# --- Matriks A (spatial + temporal) ---

# Koordinat (UTM) untuk semua observasi
coords_all <- as.matrix(data.frame(
  x = Stunting_sf_utm$UTM_x,
  y = Stunting_sf_utm$UTM_y
))

# Buat index tahun berurutan 1..T (tidak pakai 2019 dst.)
years <- sort(unique(na.omit(Stunting_sf_utm$tahun)))
Stunting_sf_utm$year_id <- match(Stunting_sf_utm$tahun, years)  # 1..T

# Cek wajib:
stopifnot(nrow(coords_all) == nrow(Stunting_sf_utm))         # panjang sama
stopifnot(all(!is.na(Stunting_sf_utm$year_id)))              # tidak ada NA
stopifnot(all(Stunting_sf_utm$year_id >= 1))                 # mulai 1
stopifnot(max(Stunting_sf_utm$year_id) == length(years))     # max = n.group
n_year <- length(years)


A_est <- inla.spde.make.A(
  mesh = stunting_mesh,
  loc = coords_all,
  group = Stunting_sf_utm$year_id,   # tahun = group temporal
  n.group = length(unique(Stunting_sf_utm$year_id))
)
dim(A_est)

## [1]  2056 10852

=========================

6.3 Control parameter for INLA computation

=========================

control <- list( compute = list( dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE, config = TRUE, return.marginals.predictor = TRUE ) )

=========================

6.4 Define prior distribution (Penalized Complexity prior)

=========================

6.4.1 Grid 50 km × 50 km

(bisa diganti 10 km kalau mau mesh sangat halus)

Grid_utm50 <- as.data.frame(st_coordinates( st_make_grid(indo_utm, cellsize = 50000, what = “centers”) )) m1 <- nrow(Grid_utm50) Grid_utm50 <- as.matrix(Grid_utm50)

GroupTime <- rep(unique(Stunting_sf_utm$tahun), each = m1) AllGrid50 <- kronecker(matrix(1, length(unique(Stunting_sf_utm$tahun)), 1), Grid_utm50)

prior.median.sd = 0.1 prior.median.range = 50000 # 50 km dalam meter

stunting_spde_50_01 <- inla.spde2.pcmatern( mesh = stunting_mesh, prior.range = c(prior.median.range, 0.5), prior.sigma = c(prior.median.sd, 0.5) )

n_year <- length(unique(Stunting_sf_utm$tahun)) s_index_50_01 <- inla.spde.make.index( name = "spatial.field", n.spde = stunting_spde_50_01$n.spde, n.group = n_year )

— Stack observed (estimation) —

stack_est_50_01 <- inla.stack( data = list(y = Stunting_sf_utm$prev_stunting), A = list(1, A_est), effects = list( data.frame( Intercept = 1, Tahun = Stunting_sf_utm$tahun, UTM_x = Stunting_sf_utm$UTM_x / 1000, UTM_y = Stunting_sf_utm$UTM_y / 1000, PDRB = Stunting_sf_utm$PDRB, Kemiskinan = Stunting_sf_utm$Tingkat.Kemiskinan, Pendidikan = Stunting_sf_utm$Tingkat.Pendidikan.Ibu ), s_index_50_01 ), tag = “est” )

— Matriks A prediksi (grid tahunan) —

A_pred_50 <- inla.spde.make.A( mesh = stunting_mesh, loc = as.matrix(AllGrid50), group = GroupTime, n.group = n_year )

indo_union <- st_union(indo_utm) |> st_make_valid() grid50_ctr <- st_make_grid(indo_union, cellsize = 50000, what = “centers”, square = TRUE) grid50_in <- grid50_ctr[lengths(st_within(grid50_ctr, indo_union)) > 0] grid_xy <- st_coordinates(grid50_in) # kolom “X”, “Y” m1 <- nrow(grid_xy)

Buat semua kombinasi grid × tahun

AllGrid50 <- kronecker(matrix(1, n_year, 1), as.matrix(grid_xy)) GroupTime <- rep(1:n_year, each = m1)

Cek kesesuaian

stopifnot(nrow(AllGrid50) == length(GroupTime)) stopifnot(all(GroupTime %in% 1:n_year)) stopifnot(ncol(AllGrid50) == 2)

Bangun matriks A prediksi

A_pred_50 <- inla.spde.make.A( mesh = stunting_mesh, loc = AllGrid50, group = GroupTime, n.group = n_year )

Cov_pred <- data.frame( UTM_x = AllGrid50[, 1] / 1000, UTM_y = AllGrid50[, 2] / 1000, Tahun = GroupTime, PDRB = NA, Kemiskinan = NA, Pendidikan = NA )

— Stack prediksi —

stack_pred_50_01 <- inla.stack( data = list(y = NA), A = list(1, A_pred_50), effects = list( data.frame( Intercept = rep(1, nrow(AllGrid50)), Tahun = GroupTime, UTM_x = Cov_pred$UTM_x, UTM_y = Cov_pred$UTM_y, PDRB = Cov_pred$PDRB, Kemiskinan = Cov_pred$Kemiskinan, Pendidikan = Cov_pred$Pendidikan ), s_index_50_01 ), tag = “pred” )

stack_50_01 <- inla.stack(stack_est_50_01, stack_pred_50_01)

formula_50_01 <- y ~ -1 + Intercept + UTM_x + UTM_y + PDRB + Kemiskinan + Pendidikan + f(spatial.field, model = stunting_spde_50_01, group = spatial.field.group, control.group = list(model = “rw1”))

Jalankan (akan berat tergantung mesh dan grid)

output_50_01 <- inla( formula_50_01, data = inla.stack.data(stack_50_01, spde = stunting_spde_50_01), family = “gaussian”, control.compute = control$compute, control.predictor = list(A = inla.stack.A(stack_50_01), compute = TRUE) )

summary(output_50_01)

Ringkasan hasil

summary(output_50_01)$fixed summary(output_50_01)$hyperpar

DIC / WAIC

cat(“DIC :”, output_50_01$dic$dic, “”) cat(“WAIC:”, output_50_01$waic$waic, “”)

Ekstrak prediksi

id_pred <- inla.stack.index(stack_50_01, tag = “pred”)$data pred_mean <- output_50_01$summary.fitted.values[id_pred, “mean”] pred_ll <- output_50_01$summary.fitted.values[id_pred, "0.025quant"] pred_ul <- output_50_01$summary.fitted.values[id_pred, “0.975quant”]

Konversi matrix grid jadi data frame dan beri nama kolom

AllGrid50_df <- as.data.frame(AllGrid50) colnames(AllGrid50_df) <- c(“UTM_x”, “UTM_y”)

Ubah ke objek sf (koordinat UTM zona 48S)

grid50_sf <- st_as_sf(AllGrid50_df, coords = c(“UTM_x”, “UTM_y”), crs = 32748) ggplot() + geom_sf(data = indo_utm, fill = NA, color = “grey40”, linewidth = 0.3) + geom_sf(data = grid50_sf, aes(color = pred_mean), size = 1.2) + scale_color_gradientn( colors = c(“#FFFFB2”, “#FECC5C”, “#FD8D3C”, “#F03B20”, “#BD0026”), # kuning → oranye → merah name = “Pred (%)”, breaks = c(0, 10, 20, 30, 40, 50), # interval legend limits = c(0, 50), # rentang nilai labels = c(“0–10”, “10–20”, “20–30”, “30–40”, “40–50”, “>50”) ) + labs( title = “Prediksi Rata-rata Prevalensi Stunting (INLA–SPDE, Range 50 km, SD 0.1)”, subtitle = “Gradasi warna kuning–merah: semakin merah berarti prevalensi semakin tinggi”, x = “Longitude”, y = “Latitude” ) + theme_minimal(base_size = 13) + theme( legend.position = “right”, legend.title = element_text(size = 12, face = “bold”), legend.text = element_text(size = 11), plot.title = element_text(face = “bold”, size = 14) ) ```

# =========================

# 6.3 Control parameter for INLA computation

# =========================

control <- list( compute = list( dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE, config = TRUE, return.marginals.predictor = TRUE ) )

# =========================

# 6.4 Define prior distribution (Penalized Complexity prior)

# =========================

# 6.4.1 Grid 50 km × 50 km

# (bisa diganti 10 km kalau mau mesh sangat halus)

Grid_utm50 <- as.data.frame(st_coordinates( st_make_grid(indo_utm, cellsize = 50000, what = "centers") )) 
m1 <- nrow(Grid_utm50) 
Grid_utm50 <- as.matrix(Grid_utm50)

GroupTime <- rep(unique(Stunting_sf_utm$tahun), each = m1)
AllGrid50 <- kronecker(matrix(1, length(unique(Stunting_sf_utm$tahun)), 1), Grid_utm50)

prior.median.sd = 0.1 
prior.median.range = 50000 # 50 km dalam meter

stunting_spde_50_01 <- inla.spde2.pcmatern( mesh = stunting_mesh, prior.range = c(prior.median.range, 0.5), prior.sigma = c(prior.median.sd, 0.5) )

n_year <- length(unique(Stunting_sf_utm$tahun))
s_index_50_01 <- inla.spde.make.index(
  name = "spatial.field",
  n.spde = stunting_spde_50_01$n.spde, n.group = n_year )

# --- Stack observed (estimation) ---

stack_est_50_01 <- inla.stack( data = list(y = Stunting_sf_utm$prev_stunting),
                               A = list(1, A_est),
                               effects = list(
                                 data.frame(
                                   Intercept = 1,
                                   Tahun = Stunting_sf_utm$tahun, UTM_x = Stunting_sf_utm$UTM_x / 1000,
                                   UTM_y = Stunting_sf_utm$UTM_y / 1000, PDRB = Stunting_sf_utm$PDRB,
                                   Kemiskinan = Stunting_sf_utm$Tingkat.Kemiskinan, Pendidikan = Stunting_sf_utm$Tingkat.Pendidikan.Ibu ), s_index_50_01 ), tag = "est" )
# --- Matriks A prediksi (grid tahunan) ---



indo_union <- st_union(indo_utm) |> st_make_valid() 
grid50_ctr <- st_make_grid(indo_union, cellsize = 50000, what = "centers", square = TRUE) 
grid50_in <- grid50_ctr[lengths(st_within(grid50_ctr, indo_union)) > 0] 
grid_xy <- st_coordinates(grid50_in) # kolom "X", "Y" m1 <- nrow(grid_xy)

# Buat semua kombinasi grid × tahun

AllGrid50 <- kronecker(matrix(1, n_year, 1), as.matrix(grid_xy)) 
GroupTime <- rep(1:n_year, each = m1)

# Cek kesesuaian
#=================== 1. Tahun unik ===================#
years <- sort(unique(Stunting_sf_utm$tahun))
n_year <- length(years)

#=================== 2. Grid 50 km ===================#
indo_union <- st_union(indo_utm) |> st_make_valid()
grid_ctr   <- st_make_grid(indo_union, cellsize = 50000, what = "centers")
grid_in    <- grid_ctr[lengths(st_within(grid_ctr, indo_union)) > 0]

# Koordinat grid (UTM meter)
grid_xy <- st_coordinates(grid_in)   # matrix [n_grid x 2]
m_grid  <- nrow(grid_xy)

#=================== 3. Kombinasi grid x tahun ===================#
# Semua titik grid dikalikan dengan jumlah tahun (cross product)
AllGrid50 <- matrix(
  rep(as.numeric(t(grid_xy)), times = n_year),
  ncol = 2,
  byrow = TRUE
)
colnames(AllGrid50) <- c("X", "Y")

# GroupTime untuk tiap tahun
GroupTime <- rep(seq_len(n_year), each = m_grid)

#=================== 4. Uji kesesuaian ===================#
cat("Jumlah grid:", m_grid, "\n")

## Jumlah grid: 824

cat("Jumlah tahun:", n_year, "\n")

## Jumlah tahun: 4

cat("Total kombinasi:", nrow(AllGrid50), "\n")

## Total kombinasi: 3296

cat("Panjang GroupTime:", length(GroupTime), "\n")

## Panjang GroupTime: 3296

stopifnot(
  nrow(AllGrid50) == length(GroupTime),
  ncol(AllGrid50) == 2,
  all(GroupTime %in% 1:n_year)
)

# Bangun matriks A prediksi

A_pred_50 <- inla.spde.make.A( mesh = stunting_mesh, loc = AllGrid50, group = GroupTime, n.group = n_year )

Cov_pred <- data.frame( UTM_x = AllGrid50[, 1] / 1000, UTM_y = AllGrid50[, 2] / 1000, Tahun = GroupTime)


# --- Stack prediksi ---
stack_pred_50_01 <- inla.stack(
  data = list(y = NA),
  A = list(1, A_pred_50),
  effects = list(
    data.frame(
      Intercept = rep(1, nrow(AllGrid50)),
      Tahun = GroupTime,
      UTM_x = Cov_pred$UTM_x,
      UTM_y = Cov_pred$UTM_y
    ),
    s_index_50_01
  ),
  tag = "pred"
)

stack_50_01 <- inla.stack(stack_est_50_01, stack_pred_50_01)

formula_50_01 <- y ~ -1 + Intercept + UTM_x + UTM_y +
  f(spatial.field, model = stunting_spde_50_01, group = spatial.field.group,
    control.group = list(model = "rw1"))

# Jalankan (akan berat tergantung mesh dan grid)
output_50_01 <- inla(
  formula_50_01,
  data = inla.stack.data(stack_50_01, spde = stunting_spde_50_01),
  family = "gaussian",
  control.compute = control$compute,
  control.predictor = list(A = inla.stack.A(stack_50_01), compute = TRUE)
)

summary(output_50_01)

## Time used:
##     Pre = 1.13, Running = 104, Post = 2.52, Total = 107 
## Fixed effects:
##            mean     sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 0.102 31.045    -60.773    0.102     60.978 0.102   0
## UTM_x     0.002  0.044     -0.085    0.002      0.089 0.002   0
## UTM_y     0.001  0.009     -0.018    0.001      0.019 0.001   0
## 
## Random effects:
##   Name     Model
##     spatial.field SPDE2 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for the Gaussian observations 9.90e-02    0.007   8.60e-02 9.90e-02
## Range for spatial.field                 1.06e+05 7602.619   9.19e+04 1.06e+05
## Stdev for spatial.field                 4.21e+00    0.141   3.94e+00 4.21e+00
##                                         0.975quant     mode
## Precision for the Gaussian observations   1.14e-01 9.90e-02
## Range for spatial.field                   1.22e+05 1.05e+05
## Stdev for spatial.field                   4.50e+00 4.21e+00
## 
## Deviance Information Criterion (DIC) ...............: 11777.74
## Deviance Information Criterion (DIC, saturated) ....: 3275.52
## Effective number of parameters .....................: 1264.20
## 
## Watanabe-Akaike information criterion (WAIC) ...: 11975.74
## Effective number of parameters .................:  999.83
## 
## Marginal log-Likelihood:  15575.08 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

# Ringkasan hasil
summary(output_50_01)$fixed

##            mean     sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 0.102 31.045    -60.773    0.102     60.978 0.102   0
## UTM_x     0.002  0.044     -0.085    0.002      0.089 0.002   0
## UTM_y     0.001  0.009     -0.018    0.001      0.019 0.001   0

summary(output_50_01)$hyperpar

##                                               mean       sd 0.025quant
## Precision for the Gaussian observations      0.099    0.007      0.086
## Range for spatial.field                 105975.894 7602.619  91887.037
## Stdev for spatial.field                      4.214    0.141      3.943
##                                           0.5quant 0.975quant       mode
## Precision for the Gaussian observations      0.099      0.114      0.099
## Range for spatial.field                 105673.856 121806.106 105010.035
## Stdev for spatial.field                      4.212      4.499      4.207

# =========================

# 6.3 Control parameter for INLA computation

# =========================

control <- list( compute = list( dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE, config = TRUE, return.marginals.predictor = TRUE ) )

# =========================

# 6.4 Define prior distribution (Penalized Complexity prior)

# =========================

# 6.4.1 Grid 50 km × 50 km

# (bisa diganti 10 km kalau mau mesh sangat halus)

Grid_utm10 <- as.data.frame(st_coordinates( st_make_grid(indo_utm, cellsize = 10000, what = "centers") )) 
m1 <- nrow(Grid_utm10) 
Grid_utm10 <- as.matrix(Grid_utm10)

GroupTime <- rep(unique(Stunting_sf_utm$tahun), each = m1)
AllGrid10 <- kronecker(matrix(1, length(unique(Stunting_sf_utm$tahun)), 1), Grid_utm10)

prior.median.sd = 0.1 
prior.median.range = 10000 # 50 km dalam meter

stunting_spde_10_01 <- inla.spde2.pcmatern( mesh = stunting_mesh, prior.range = c(prior.median.range, 0.5), prior.sigma = c(prior.median.sd, 0.5) )

n_year <- length(unique(Stunting_sf_utm$tahun))
s_index_10_01 <- inla.spde.make.index(
  name = "spatial.field",
  n.spde = stunting_spde_10_01$n.spde, n.group = n_year )

# --- Stack observed (estimation) ---

stack_est_10_01 <- inla.stack( data = list(y = Stunting_sf_utm$prev_stunting),
                               A = list(1, A_est),
                               effects = list(
                                 data.frame(
                                   Intercept = 1,
                                   Tahun = Stunting_sf_utm$tahun, UTM_x = Stunting_sf_utm$UTM_x / 1000,
                                   UTM_y = Stunting_sf_utm$UTM_y / 1000, PDRB = Stunting_sf_utm$PDRB,
                                   Kemiskinan = Stunting_sf_utm$Tingkat.Kemiskinan, Pendidikan = Stunting_sf_utm$Tingkat.Pendidikan.Ibu ), s_index_10_01 ), tag = "est" )
# --- Matriks A prediksi (grid tahunan) ---

indo_union <- st_union(indo_utm) |> st_make_valid() 
grid10_ctr <- st_make_grid(indo_union, cellsize = 10000, what = "centers", square = TRUE) 
grid10_in <- grid10_ctr[lengths(st_within(grid10_ctr, indo_union)) > 0] 
grid_xy <- st_coordinates(grid10_in) # kolom "X", "Y" m1 <- nrow(grid_xy)

# Buat semua kombinasi grid × tahun

AllGrid50 <- kronecker(matrix(1, n_year, 1), as.matrix(grid_xy)) 
GroupTime <- rep(1:n_year, each = m1)

# Cek kesesuaian
#=================== 1. Tahun unik ===================#
years <- sort(unique(Stunting_sf_utm$tahun))
n_year <- length(years)

#=================== 2. Grid 50 km ===================#
indo_union <- st_union(indo_utm) |> st_make_valid()
grid_ctr   <- st_make_grid(indo_union, cellsize = 10000, what = "centers")
grid_in    <- grid_ctr[lengths(st_within(grid_ctr, indo_union)) > 0]

# Koordinat grid (UTM meter)
grid_xy <- st_coordinates(grid_in)   # matrix [n_grid x 2]
m_grid  <- nrow(grid_xy)

#=================== 3. Kombinasi grid x tahun ===================#
# Semua titik grid dikalikan dengan jumlah tahun (cross product)
AllGrid10 <- matrix(
  rep(as.numeric(t(grid_xy)), times = n_year),
  ncol = 2,
  byrow = TRUE
)
colnames(AllGrid10) <- c("X", "Y")

# GroupTime untuk tiap tahun
GroupTime <- rep(seq_len(n_year), each = m_grid)

#=================== 4. Uji kesesuaian ===================#
cat("Jumlah grid:", m_grid, "\n")

## Jumlah grid: 20843

cat("Jumlah tahun:", n_year, "\n")

## Jumlah tahun: 4

cat("Total kombinasi:", nrow(AllGrid50), "\n")

## Total kombinasi: 83372

cat("Panjang GroupTime:", length(GroupTime), "\n")

## Panjang GroupTime: 83372

stopifnot(
  nrow(AllGrid10) == length(GroupTime),
  ncol(AllGrid10) == 2,
  all(GroupTime %in% 1:n_year)
)

# Bangun matriks A prediksi

A_pred_10 <- inla.spde.make.A( mesh = stunting_mesh, loc = AllGrid50, group = GroupTime, n.group = n_year )

Cov_pred <- data.frame( UTM_x = AllGrid10[, 1] / 1000, UTM_y = AllGrid10[, 2] / 1000, Tahun = GroupTime)


# --- Stack prediksi ---
stack_pred_10_01 <- inla.stack(
  data = list(y = NA),
  A = list(1, A_pred_10),
  effects = list(
    data.frame(
      Intercept = rep(1, nrow(AllGrid10)),
      Tahun = GroupTime,
      UTM_x = Cov_pred$UTM_x,
      UTM_y = Cov_pred$UTM_y
    ),
    s_index_10_01
  ),
  tag = "pred"
)

stack_10_01 <- inla.stack(stack_est_10_01, stack_pred_10_01)

formula_10_01 <- y ~ -1 + Intercept + UTM_x + UTM_y +
  f(spatial.field, model = stunting_spde_10_01, group = spatial.field.group,
    control.group = list(model = "rw1"))

# Jalankan (akan berat tergantung mesh dan grid)
output_10_01 <- inla(
  formula_10_01,
  data = inla.stack.data(stack_10_01, spde = stunting_spde_10_01),
  family = "gaussian",
  control.compute = control$compute,
  control.predictor = list(A = inla.stack.A(stack_10_01), compute = TRUE)
)

summary(output_10_01)

## Time used:
##     Pre = 1.33, Running = 125, Post = 5.8, Total = 132 
## Fixed effects:
##            mean     sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 0.102 31.045    -60.775    0.102     60.979 0.102   0
## UTM_x     0.002  0.044     -0.085    0.002      0.089 0.002   0
## UTM_y     0.001  0.009     -0.018    0.001      0.019 0.001   0
## 
## Random effects:
##   Name     Model
##     spatial.field SPDE2 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for the Gaussian observations 9.90e-02    0.007   8.60e-02 9.90e-02
## Range for spatial.field                 1.06e+05 7597.878   9.17e+04 1.06e+05
## Stdev for spatial.field                 4.21e+00    0.141   3.94e+00 4.21e+00
##                                         0.975quant     mode
## Precision for the Gaussian observations   1.15e-01 9.90e-02
## Range for spatial.field                   1.22e+05 1.05e+05
## Stdev for spatial.field                   4.50e+00 4.21e+00
## 
## Deviance Information Criterion (DIC) ...............: 11777.75
## Deviance Information Criterion (DIC, saturated) ....: 3272.66
## Effective number of parameters .....................: 1261.89
## 
## Watanabe-Akaike information criterion (WAIC) ...: 11977.72
## Effective number of parameters .................:  999.54
## 
## Marginal log-Likelihood:  15573.73 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

# Ringkasan hasil
summary(output_10_01)$fixed

##            mean     sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 0.102 31.045    -60.775    0.102     60.979 0.102   0
## UTM_x     0.002  0.044     -0.085    0.002      0.089 0.002   0
## UTM_y     0.001  0.009     -0.018    0.001      0.019 0.001   0

summary(output_10_01)$hyperpar

##                                               mean       sd 0.025quant
## Precision for the Gaussian observations      0.099    0.007      0.086
## Range for spatial.field                 105794.186 7597.878  91676.954
## Stdev for spatial.field                      4.213    0.141      3.942
##                                           0.5quant 0.975quant       mode
## Precision for the Gaussian observations      0.099      0.115      0.099
## Range for spatial.field                 105505.014 121576.990 104894.027
## Stdev for spatial.field                      4.210      4.499      4.205

PERAMALAN DAN PEMETAAN PREVALENSI STUNTING DI INDONESIA MENGGUNAKAN

Kiki Amelia

2025-11-07

=========================

6.3 Control parameter for INLA computation

=========================

=========================

6.4 Define prior distribution (Penalized Complexity prior)

=========================

6.4.1 Grid 50 km × 50 km

(bisa diganti 10 km kalau mau mesh sangat halus)

— Stack observed (estimation) —

— Matriks A prediksi (grid tahunan) —

Buat semua kombinasi grid × tahun

Cek kesesuaian

Bangun matriks A prediksi

— Stack prediksi —

Jalankan (akan berat tergantung mesh dan grid)

Ringkasan hasil

DIC / WAIC

Ekstrak prediksi

Konversi matrix grid jadi data frame dan beri nama kolom

Ubah ke objek sf (koordinat UTM zona 48S)