1. Persiapan Paket dan Data

# Paket yang dibutuhkan
library(sf)            # untuk data spasial vector (.shp)
library(terra)         # untuk data raster (.tif)
library(dplyr)         # untuk manipulasi data
library(readr)         # untuk membaca file CSV
library(tidyr)         # untuk fungsi separate
library(spOccupancy)   # untuk modeling stPGOcc
library(stringr)       # untuk manipulasi string
# Load shapefile: titik perjumpaan Kelempiau - Hylobates albibarbis, jalur patroli, dan grid kawasan
Kelempiau <- st_read("input/Kelempiau_2020-2024.shp")                 # titik perjumpaan Kelempiau
## Reading layer `Kelempiau_2020-2024' from data source 
##   `D:\ARIF\WORK\YPI\MEL\R\tesis_hendra\input\Kelempiau_2020-2024.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 129 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 490021.3 ymin: 9939979 xmax: 518047 ymax: 9960096
## Projected CRS: WGS 84 / UTM zone 49S
patrol_track <- st_read("input/track_2020-2024.shp")         # track patroli
## Reading layer `track_2020-2024' from data source 
##   `D:\ARIF\WORK\YPI\MEL\R\tesis_hendra\input\track_2020-2024.shp' 
##   using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 159 features and 3 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: 482943.5 ymin: 9933999 xmax: 523443.8 ymax: 9965050
## z_range:       zmin: 1.595325e+12 zmax: 1.726392e+12
## Projected CRS: WGS 84 / UTM zone 49S
protected_area_grid <- st_read("input/grid_hlgn.shp")        # grid kawasan 2x2 km
## Reading layer `grid_hlgn' from data source 
##   `D:\ARIF\WORK\YPI\MEL\R\tesis_hendra\input\grid_hlgn.shp' using driver `ESRI Shapefile'
## Simple feature collection with 693 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 478628.2 ymin: 9893919 xmax: 555381.1 ymax: 9966136
## Projected CRS: WGS 84 / UTM zone 49S
# Load dan susun raster variabel lingkungan
var_envi <- c("output/Elevasi.tif", "output/Jarak Dari Jalan.tif", "output/Jarak Dari Permukiman.tif",
              "output/Jarak Dari Sungai.tif", "output/LAI.tif", "output/NDVI.tif", 
              "output/Slope.tif", "output/TUPLAH.tif")
rasters_var_env <- rast(var_envi)
names(rasters_var_env) <- c("elev", "dist_jalan", "dist_permukiman", 
                            "dist_sungai", "lai", "ndvi", "slope", "tuplah")

2. Proses Data Spasial

# Tambah kolom tahun-bulan (format "YYYY-MM") untuk agregasi waktu
Kelempiau <- Kelempiau %>% mutate(year_month = format(as.Date(date), "%Y-%m"))
patrol_track <- patrol_track %>% mutate(year_month = format(as.Date(date), "%Y-%m"))

# Gabungkan titik Kelempiau ke grid → identifikasi grid dengan perjumpaan
Kelempiau_grid <- st_join(Kelempiau, protected_area_grid, join = st_within) %>%
  filter(!is.na(id_grid)) %>%
  st_drop_geometry() %>%
  select(id_grid, year_month) %>%
  distinct() %>%
  mutate(detection = 1)  # 1 = Kelempiau terdeteksi

# Potong jalur patroli berdasarkan grid → identifikasi grid yang dipatroli
patrol_grid <- st_intersection(patrol_track, protected_area_grid) %>%
  st_drop_geometry() %>%
  select(id_grid, year_month) %>%
  distinct() %>%
  mutate(effort = 1)  # 1 = ada effort patroli

# Gabungkan data effort patroli dan deteksi Kelempiau
detection_effort <- full_join(patrol_grid, Kelempiau_grid, by = c("id_grid", "year_month")) %>%
  mutate(detection = ifelse(is.na(detection) & effort == 1, 0, detection)) %>% # tidak deteksi → 0
  filter(!is.na(effort)) %>%  # buang data tanpa effort
  arrange(id_grid, year_month)

# Hapus komponen spasial sepenuhnya
detection_effort_clean <- detection_effort %>% st_drop_geometry()
write.csv(detection_effort_clean, "output/patrol_Kelempiau_detection_effort.csv", row.names = FALSE)

3. Bentuk Array Deteksi

grid_ids <- sort(unique(detection_effort$id_grid))
periods <- sort(unique(detection_effort$year_month))

detection_effort <- detection_effort %>%
  mutate(site_idx = match(id_grid, grid_ids),
         time_idx = match(year_month, periods))

n_sites <- length(grid_ids)
n_time <- length(periods)
y_array <- array(NA, dim = c(n_sites, n_time, 1))

for (i in 1:nrow(detection_effort)) {
  y_array[detection_effort$site_idx[i],
          detection_effort$time_idx[i], 1] <- detection_effort$detection[i]
}

4. Ekstraksi Kovariat Lingkungan

# Konversi grid ke format terra
protected_area_grid_terra <- vect(protected_area_grid)

# Crop & mask raster agar hanya area grid
rasters_var_env_crop <- crop(rasters_var_env, protected_area_grid_terra)
rasters_var_env_mask <- mask(rasters_var_env_crop, protected_area_grid_terra)
## |---------|---------|---------|---------|=========================================                                          
# Ekstrak nilai rata-rata raster per grid
val_env <- terra::extract(rasters_var_env_mask, protected_area_grid_terra, fun = mean, na.rm = TRUE)
grid_env_covariates <- cbind(protected_area_grid_terra, val_env)

# Filter NA dan simpan id_grid yang valid
grid_env_df <- grid_env_covariates %>%
  as.data.frame() %>%
  filter(!if_any(c(elev, dist_jalan, dist_permukiman, dist_sungai, lai, ndvi, slope, tuplah), is.na))

valid_grids <- grid_env_df$id_grid

5. Persiapan Data stPGOcc

data_Kelempiau <- read_csv("output/patrol_Kelempiau_detection_effort.csv") # atau detection_effort_clean
data_Kelempiau_clean <- data_Kelempiau %>% filter(id_grid %in% valid_grids)

data_Kelempiau2 <- data_Kelempiau_clean %>%
  arrange(id_grid, year_month) %>%
  group_by(id_grid, year_month) %>%
  mutate(visit = row_number()) %>%
  ungroup() %>%
  separate(year_month, into = c("tahun", "bulan"), sep = "-") %>%
  mutate(waktu = paste0("y", tahun, "_m", sprintf("%02d", as.integer(bulan))))

grids <- sort(unique(data_Kelempiau2$id_grid))
waktus <- sort(unique(data_Kelempiau2$waktu))
max_visit <- max(data_Kelempiau2$visit)

# Inisialisasi array deteksi dan effort
deteksi_array <- array(NA, dim = c(length(grids), length(waktus), max_visit))
effort_array <- array(NA, dim = c(length(grids), length(waktus), max_visit))

for (i in seq_along(grids)) {
  for (j in seq_along(waktus)) {
    subset_df <- filter(data_Kelempiau2, id_grid == grids[i], waktu == waktus[j])
    if (nrow(subset_df) > 0) {
      deteksi_array[i, j, 1:nrow(subset_df)] <- subset_df$detection
      effort_array[i, j, 1:nrow(subset_df)] <- subset_df$effort
    }
  }
}

# Ganti NA jadi 0 (tidak ada deteksi/effort)
deteksi_array[is.na(deteksi_array)] <- 0
effort_array[is.na(effort_array)] <- 0

6. Siapkan Kovariat dan Koordinat

site_covs <- grid_env_df %>%
  filter(id_grid %in% grid_ids) %>%
  arrange(match(id_grid, grid_ids)) %>%
  select(elev, dist_jalan, dist_permukiman, dist_sungai, lai, ndvi, slope, tuplah)

site_covs_scaled <- scale(site_covs)

# Buat koordinat centroid dari grid
grid_coords <- st_coordinates(st_centroid(st_as_sf(protected_area_grid_terra))) %>%
  as.data.frame() %>%
  rename(x = X, y = Y) %>%
  slice(match(grids, protected_area_grid_terra$id_grid))

7. Fit Model Okupansi

set.seed(123)

# Tentukan formula
occ_formula <- ~ elev + dist_jalan + dist_permukiman + dist_sungai + lai + ndvi + slope + tuplah
det_formula <- ~ effort

model_Kelempiau <- stPGOcc(
  occ.formula = occ_formula,
  det.formula = det_formula,
  data = list(
    y = deteksi_array,
    occ.covs = as.data.frame(site_covs_scaled),
    det.covs = list(effort = effort_array),
    coords = grid_coords
  ),
  n.batch = 1000,
  batch.length = 10,
  n.burn = 2000,
  n.thin = 5,
  n.omp.threads = 4,
  verbose = TRUE,
  n.chains = 3  # GUNAKAN ≥3 CHAIN UNTUK KONVERGENSI
)
## ----------------------------------------
##  Preparing the data
## ----------------------------------------
## ----------------------------------------
##  Building the neighbor list
## ----------------------------------------
## ----------------------------------------
## Building the neighbors of neighbors list
## ----------------------------------------
## ----------------------------------------
##  Model description
## ----------------------------------------
## Spatial NNGP Multi-season Occupancy Model with Polya-Gamma latent
## variable fit with 130 sites and 49 primary time periods.
## 
## Samples per chain: 10000 (1000 batches of length 10)
## Burn-in: 2000 
## Thinning Rate: 5 
## Number of Chains: 3 
## Total Posterior Samples: 4800 
## 
## Using the exponential spatial correlation model.
## 
## Using 15 nearest neighbors.
## 
## Source compiled with OpenMP support and model fit using 4 thread(s).
## 
## Adaptive Metropolis with target acceptance rate: 43.0
## ----------------------------------------
##  Chain 1
## ----------------------------------------
## Sampling ... 
## Batch: 100 of 1000, 10.00%
##  Parameter   Acceptance  Tuning
##  phi     60.0        0.62500
## -------------------------------------------------
## Batch: 200 of 1000, 20.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.54335
## -------------------------------------------------
## Batch: 300 of 1000, 30.00%
##  Parameter   Acceptance  Tuning
##  phi     70.0        0.71892
## -------------------------------------------------
## Batch: 400 of 1000, 40.00%
##  Parameter   Acceptance  Tuning
##  phi     30.0        0.63763
## -------------------------------------------------
## Batch: 500 of 1000, 50.00%
##  Parameter   Acceptance  Tuning
##  phi     80.0        0.81058
## -------------------------------------------------
## Batch: 600 of 1000, 60.00%
##  Parameter   Acceptance  Tuning
##  phi     70.0        0.73345
## -------------------------------------------------
## Batch: 700 of 1000, 70.00%
##  Parameter   Acceptance  Tuning
##  phi     80.0        0.74826
## -------------------------------------------------
## Batch: 800 of 1000, 80.00%
##  Parameter   Acceptance  Tuning
##  phi     40.0        0.62500
## -------------------------------------------------
## Batch: 900 of 1000, 90.00%
##  Parameter   Acceptance  Tuning
##  phi     90.0        0.69073
## -------------------------------------------------
## Batch: 1000 of 1000, 100.00%
## ----------------------------------------
##  Chain 2
## ----------------------------------------
## Sampling ... 
## Batch: 100 of 1000, 10.00%
##  Parameter   Acceptance  Tuning
##  phi     40.0        0.66365
## -------------------------------------------------
## Batch: 200 of 1000, 20.00%
##  Parameter   Acceptance  Tuning
##  phi     70.0        0.57695
## -------------------------------------------------
## Batch: 300 of 1000, 30.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.65051
## -------------------------------------------------
## Batch: 400 of 1000, 40.00%
##  Parameter   Acceptance  Tuning
##  phi     60.0        0.58860
## -------------------------------------------------
## Batch: 500 of 1000, 50.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.65051
## -------------------------------------------------
## Batch: 600 of 1000, 60.00%
##  Parameter   Acceptance  Tuning
##  phi     10.0        0.51171
## -------------------------------------------------
## Batch: 700 of 1000, 70.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.53259
## -------------------------------------------------
## Batch: 800 of 1000, 80.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.50158
## -------------------------------------------------
## Batch: 900 of 1000, 90.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.60050
## -------------------------------------------------
## Batch: 1000 of 1000, 100.00%
## ----------------------------------------
##  Chain 3
## ----------------------------------------
## Sampling ... 
## Batch: 100 of 1000, 10.00%
##  Parameter   Acceptance  Tuning
##  phi     60.0        0.54335
## -------------------------------------------------
## Batch: 200 of 1000, 20.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.55433
## -------------------------------------------------
## Batch: 300 of 1000, 30.00%
##  Parameter   Acceptance  Tuning
##  phi     40.0        0.53259
## -------------------------------------------------
## Batch: 400 of 1000, 40.00%
##  Parameter   Acceptance  Tuning
##  phi     60.0        0.58860
## -------------------------------------------------
## Batch: 500 of 1000, 50.00%
##  Parameter   Acceptance  Tuning
##  phi     40.0        0.55433
## -------------------------------------------------
## Batch: 600 of 1000, 60.00%
##  Parameter   Acceptance  Tuning
##  phi     60.0        0.57695
## -------------------------------------------------
## Batch: 700 of 1000, 70.00%
##  Parameter   Acceptance  Tuning
##  phi     50.0        0.54335
## -------------------------------------------------
## Batch: 800 of 1000, 80.00%
##  Parameter   Acceptance  Tuning
##  phi     70.0        0.60050
## -------------------------------------------------
## Batch: 900 of 1000, 90.00%
##  Parameter   Acceptance  Tuning
##  phi     30.0        0.61263
## -------------------------------------------------
## Batch: 1000 of 1000, 100.00%

8. Simpan & Ringkasan Model

# Simpan model untuk keperluan validasi / penggunaan ulang
saveRDS(model_Kelempiau, file = "output/model_Kelempiau_stpgocc.rds")

# Load kembali bila perlu:
# model_Kelempiau <- readRDS("output/model_Kelempiau_stpgocc.rds")

9. Visualisasi

# MULAI VISUALISASI
library(ggplot2)
library(cowplot)
library(coda)    # Untuk traceplot MCMC

# Cek dan investigasi data
# Struktur model
str(model_Kelempiau)
## List of 28
##  $ rhat          :List of 3
##   ..$ beta : num [1:9] 2.08 1.06 1.18 1.28 1.04 ...
##   ..$ alpha: num [1:2] 1.32 1.08
##   ..$ theta: num [1:2] 1.54 1.38
##  $ beta.samples  : 'mcmc' num [1:4800, 1:9] -3.23 -3.35 -3.11 -3.19 -3.17 ...
##   ..- attr(*, "mcpar")= num [1:3] 1 4800 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:9] "(Intercept)" "elev" "dist_jalan" "dist_permukiman" ...
##  $ alpha.samples : 'mcmc' num [1:4800, 1:2] -4.82 -4.78 -5.41 -5.23 -5.41 ...
##   ..- attr(*, "mcpar")= num [1:3] 1 4800 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "(Intercept)" "effort"
##  $ theta.samples : 'mcmc' num [1:4800, 1:2] 2.05 1.76 1.94 2.44 2.8 ...
##   ..- attr(*, "mcpar")= num [1:3] 1 4800 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "sigma.sq" "phi"
##  $ coords        : num [1:130, 1:2] 515628 517628 519628 521628 515628 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "x" "y"
##  $ X             : num [1:130, 1:49, 1:9] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : chr [1:9] "(Intercept)" "elev" "dist_jalan" "dist_permukiman" ...
##  $ X.re          : int[1:130, 1:49, 0 ] 
##  $ w.samples     : 'mcmc' num [1:4800, 1:130] -1.578 -0.805 -1.514 -1.67 -0.109 ...
##   ..- attr(*, "mcpar")= num [1:3] 1 4800 1
##  $ z.samples     : num [1:4800, 1:130, 1:49] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : NULL
##  $ psi.samples   : num [1:4800, 1:130, 1:49] 0.00288 0.00363 0.00296 0.00228 0.00559 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : NULL
##  $ like.samples  : num [1:4800, 1:130, 1:49] 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : NULL
##  $ y             : num [1:130, 1:49, 1] 0 0 0 0 0 0 0 0 0 0 ...
##  $ X.p           : num [1:6370, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "(Intercept)" "effort"
##  $ X.p.re        : int[1:6370, 0 ] 
##  $ ESS           :List of 3
##   ..$ beta : Named num [1:9] 18.6 143.7 56 74.7 89 ...
##   .. ..- attr(*, "names")= chr [1:9] "(Intercept)" "elev" "dist_jalan" "dist_permukiman" ...
##   ..$ alpha: Named num [1:2] 257 359
##   .. ..- attr(*, "names")= chr [1:2] "(Intercept)" "effort"
##   ..$ theta: Named num [1:2] 24.3 114.4
##   .. ..- attr(*, "names")= chr [1:2] "sigma.sq" "phi"
##  $ call          : language stPGOcc(occ.formula = occ_formula, det.formula = det_formula, data = list(y = deteksi_array,      occ.covs = as.d| __truncated__ ...
##  $ n.samples     : int 10000
##  $ n.neighbors   : int 15
##  $ cov.model.indx: int 0
##  $ type          : chr "NNGP"
##  $ n.post        : int 1600
##  $ n.thin        : int 5
##  $ n.burn        : int 2000
##  $ n.chains      : int 3
##  $ ar1           : logi FALSE
##  $ pRE           : logi FALSE
##  $ psiRE         : logi FALSE
##  $ run.time      : 'proc_time' Named num [1:5] 245.4 33.7 234 NA NA
##   ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ...
##  - attr(*, "class")= chr "stPGOcc"
# Ringkasan parameter
summary(model_Kelempiau)
## 
## Call:
## stPGOcc(occ.formula = occ_formula, det.formula = det_formula, 
##     data = list(y = deteksi_array, occ.covs = as.data.frame(site_covs_scaled), 
##         det.covs = list(effort = effort_array), coords = grid_coords), 
##     n.batch = 1000, batch.length = 10, n.omp.threads = 4, verbose = TRUE, 
##     n.burn = 2000, n.thin = 5, n.chains = 3)
## 
## Samples per Chain: 10000
## Burn-in: 2000
## Thinning Rate: 5
## Number of Chains: 3
## Total Posterior Samples: 4800
## Run Time (min): 3.9007
## 
## Occurrence (logit scale): 
##                    Mean     SD    2.5%     50%   97.5%   Rhat ESS
## (Intercept)     -2.3567 0.8742 -3.9159 -2.4418 -0.6376 2.0770  19
## elev            -0.0104 0.5662 -1.0882 -0.0387  1.2462 1.0557 144
## dist_jalan      -0.3390 0.7343 -2.0359 -0.2902  0.9541 1.1812  56
## dist_permukiman  1.4268 0.6280  0.1920  1.3766  2.7138 1.2829  75
## dist_sungai      0.0771 0.4411 -0.7537  0.0669  0.9602 1.0442  89
## lai              0.0312 0.4732 -0.8343  0.0156  1.0212 1.0844  98
## ndvi            -0.0102 0.6184 -1.3832  0.0273  1.1931 1.0366  68
## slope            0.0469 0.3788 -0.6149  0.0200  0.8642 1.0292 119
## tuplah           0.7548 0.5417 -0.1830  0.7230  1.9001 1.0655 147
## 
## Detection (logit scale): 
##                Mean     SD    2.5%     50%   97.5%   Rhat ESS
## (Intercept) -5.7040 0.5673 -6.9080 -5.6654 -4.6855 1.3190 257
## effort       5.5163 0.5965  4.4552  5.4857  6.8070 1.0826 359
## 
## Spatial Covariance: 
##            Mean     SD   2.5%    50%   97.5%   Rhat ESS
## sigma.sq 5.4799 4.0641 1.3520 4.1959 16.8274 1.5418  24
## phi      0.0002 0.0001 0.0001 0.0002  0.0004 1.3798 114
# Nama grid dan dimensi
head(model_Kelempiau$coords)
##             x       y
## [1,] 515628.2 9934919
## [2,] 517628.2 9934919
## [3,] 519628.2 9934919
## [4,] 521628.2 9934919
## [5,] 515628.2 9936919
## [6,] 519628.2 9936919
# Jika ada shapefile/spatial grid:
grid_viz <- st_read("input/grid_hlgn.shp")
## Reading layer `grid_hlgn' from data source 
##   `D:\ARIF\WORK\YPI\MEL\R\tesis_hendra\input\grid_hlgn.shp' using driver `ESRI Shapefile'
## Simple feature collection with 693 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 478628.2 ymin: 9893919 xmax: 555381.1 ymax: 9966136
## Projected CRS: WGS 84 / UTM zone 49S
# untuk cek berapa dimensi
dim(model_Kelempiau$psi.samples)
## [1] 4800  130   49
# 1. Plot α (Alpha Plot) – Koefisien Deteksi -----------------------------------
beta_mcmc <- model_Kelempiau$beta.samples  # mcmc object, dim: 4800 x 9

beta_summary <- summary(beta_mcmc)
beta_sum_df <- data.frame(
  Parameter = colnames(beta_mcmc),
  Mean = beta_summary$statistics[,"Mean"],
  SD = beta_summary$statistics[,"SD"],
  Lower = beta_summary$quantiles[,"2.5%"],
  Upper = beta_summary$quantiles[,"97.5%"]
)

alpha_mcmc <- model_Kelempiau$alpha.samples # mcmc object 4800 x 2
alpha_summary <- summary(alpha_mcmc)
alpha_sum_df <- data.frame(
  Parameter = colnames(alpha_mcmc),
  Mean = alpha_summary$statistics[,"Mean"],
  SD = alpha_summary$statistics[,"SD"],
  Lower = alpha_summary$quantiles[,"2.5%"],
  Upper = alpha_summary$quantiles[,"97.5%"]
)

p_alpha <- ggplot(alpha_sum_df, aes(x=Parameter, y=Mean)) +
  geom_point(size = 3, color = "gray30") +  # titik lebih ringan) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, linewidth = 1, 
                color = "gray30") +
  theme_minimal(base_size = 14) +
  ggtitle("Alpha Plot – Koefisien Deteksi") +
  ylab("Nilai Koefisien (logit scale)") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # judul di tengah
    axis.text.x = element_text(angle = 45, hjust = 1, size = 13),
    axis.text.y = element_text(size = 13),
    panel.grid.major = element_line(color = "lightblue", size = 0.5),
    panel.grid.minor = element_line(color = "lightblue", size = 0.25)
  )
p_alpha # <- ditampilkan sebagai output

# 2. Visualisasi Koefisien Kovariat Lingkungan dari Proses Okupansi

# Tambahkan kolom penanda jenis proses (Okupansi atau Deteksi) — meskipun hanya akan digunakan 'Occ'
beta_sum_df$Type <- "Occ"
alpha_sum_df$Type <- "Det"

# Filter hanya parameter okupansi yang merupakan kovariat lingkungan (tanpa intercept)
# dan ubah nama parameter menjadi label yang lebih informatif
coef_df <- beta_sum_df %>%
  filter(Parameter != "(Intercept)") %>%
  mutate(
    CleanParam = recode(Parameter,
                        "elev" = "Elevasi",
                        "dist_jalan" = "J_Jalan",
                        "dist_permukiman" = "J_Permukiman",
                        "dist_sungai" = "J_Sungai",
                        "lai" = "LAI",
                        "ndvi" = "NDVI",
                        "slope" = "Lereng",
                        "tuplah" = "Tutupan Lahan"),
    CleanParam = factor(CleanParam, levels = unique(CleanParam))  # Atur urutan tampil di plot
  )

# Buat plot dot-error bar untuk koefisien kovariat lingkungan (tanpa intercept)
p_coefs <- ggplot(coef_df, aes(x = CleanParam, y = Mean)) +
  geom_point(size = 3, color = "gray30") +  # Titik koefisien
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, linewidth = 1, 
                color = "gray30") +  # Error bar
  theme_minimal(base_size = 14) +
  labs(
    title = "Koefisien Kovariat Lingkungan pada Proses Okupansi",
    x = "Kovariat Lingkungan",
    y = "Nilai Koefisien (logit scale)"
  ) +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),  # Judul di tengah
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),  # Label x miring
    axis.text.y = element_text(size = 13),
    panel.grid.major = element_line(color = "lightblue", size = 0.5),
    panel.grid.minor = element_line(color = "lightblue", size = 0.25)
  )

# Tampilkan plot
p_coefs

# 3. Traceplot untuk parameter spasial & intercept --------------------------------

# Extract mcmc samples for intercept dan spatial params (theta)
beta_mcmc_list <- as.mcmc.list(lapply(1:model_Kelempiau$n.chains, function(i){
  window(model_Kelempiau$beta.samples, start=1, end=model_Kelempiau$n.samples/model_Kelempiau$n.chains, thin=model_Kelempiau$n.thin)
}))

theta_mcmc_list <- as.mcmc.list(lapply(1:model_Kelempiau$n.chains, function(i){
  window(model_Kelempiau$theta.samples, start=1, end=model_Kelempiau$n.samples/model_Kelempiau$n.chains, thin=model_Kelempiau$n.thin)
}))

# Traceplot Theta (spatial) - default judul berdasarkan nama parameter
par(mfrow = c(1, 2))  # atur layout sesuai jumlah parameter
traceplot(theta_mcmc_list)

# Traceplot Semua Koefisien Beta (Okupansi)
par(mfrow = c(3, 3),             # 3x3 untuk 9 parameter
    mar = c(5, 5, 4, 2),         
    cex.main = 1.6,              
    cex.lab = 1.4,               
    cex.axis = 1.2)              

param_names <- colnames(model_Kelempiau$beta.samples)

for (param in param_names) {
  traceplot(model_Kelempiau$beta.samples[, param], 
            main = paste("Traceplot:", param),
            ylab = "Nilai Sampel",
            xlab = "Iterasi")
}

# Simpan traceplot theta secara terpisah (menjadi 2 file .png)
# Buat folder jika belum ada
dir.create("R/plot/traceplots_theta", recursive = TRUE, showWarnings = FALSE)

# Ambil nama parameter dari theta
theta_param_names <- colnames(as.matrix(theta_mcmc_list[[1]]))

# Simpan traceplot setiap parameter theta
for (param in theta_param_names) {
  png(filename = paste0("R/plot/traceplots_theta/Traceplot_Theta_", param, ".png"),
      width = 1600, height = 1200, res = 150)
  
  traceplot(theta_mcmc_list[, param], 
            main = paste("Traceplot Theta:", param),
            ylab = "Nilai Sampel",
            xlab = "Iterasi")
  
  dev.off()
}


# Simpan traceplot beta secara terpisah (menjadi 9 file .png)
# Buat folder jika belum ada
dir.create("R/plot/traceplots_beta", recursive = TRUE, showWarnings = FALSE)

# Ambil nama parameter beta
param_names <- colnames(model_Kelempiau$beta.samples)

# Simpan traceplot setiap parameter beta
for (param in param_names) {
  png(filename = paste0("R/plot/traceplots_beta/Traceplot_Beta_", param, ".png"),
      width = 1600, height = 1200, res = 150)
  
  traceplot(model_Kelempiau$beta.samples[, param], 
            main = paste("Traceplot Beta:", param),
            ylab = "Nilai Sampel",
            xlab = "Iterasi")
  
  dev.off()
}
par(mfrow = c(1, 1))  # reset layout
# 4. Peta probabilitas okupansi (ψ_mean) atas grid ------------------------------

# Hitung mean probabilitas okupansi tiap grid pada periode tertentu (misal waktu terakhir 49)
psi_samples <- model_Kelempiau$psi.samples  # 4800 x 130 x 49
psi_mean_last <- apply(psi_samples[,,49], 2, mean)   # rata-rata sampel untuk tiap grid pada waktu ke-49
psi_sd_last <- apply(psi_samples[,,49], 2, sd)       # standar deviasi (ketidakpastian)

# Siapkan data spatial grid
grid_viz <- st_read("input/grid_hlgn.shp")
## Reading layer `grid_hlgn' from data source 
##   `D:\ARIF\WORK\YPI\MEL\R\tesis_hendra\input\grid_hlgn.shp' using driver `ESRI Shapefile'
## Simple feature collection with 693 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 478628.2 ymin: 9893919 xmax: 555381.1 ymax: 9966136
## Projected CRS: WGS 84 / UTM zone 49S
head(model_Kelempiau$coords)
##             x       y
## [1,] 515628.2 9934919
## [2,] 517628.2 9934919
## [3,] 519628.2 9934919
## [4,] 521628.2 9934919
## [5,] 515628.2 9936919
## [6,] 519628.2 9936919
grid_centroids <- st_centroid(grid_viz)
coords_shp <- st_coordinates(grid_centroids)
head(coords_shp)
##             X       Y
## [1,] 509930.5 9895444
## [2,] 511594.4 9895414
## [3,] 513587.9 9895440
## [4,] 515384.9 9895677
## [5,] 517844.6 9895658
## [6,] 519754.4 9895272
# Buat sf object titik model
model_points <- st_as_sf(as.data.frame(model_Kelempiau$coords), coords = c("x", "y"), crs = st_crs(grid_viz))

# Hitung jarak tiap grid centroid ke titik model
grid_centroids <- st_centroid(grid_viz)

# Cari grid yang jaraknya paling dekat ke minimal satu titik model
# Misal threshold jarak 5000 meter (5 km)
dist_matrix <- st_distance(grid_centroids, model_points)  # matrix jarak

min_dist <- apply(dist_matrix, 1, min)  # jarak terdekat tiap grid

# Subset grid yang jaraknya kurang dari threshold (5 km)
grid_subset <- grid_viz[min_dist < 5000, ]

# Centroid dari subset grid
grid_subset_centroids <- st_centroid(grid_subset)

# Cari nearest grid untuk tiap titik model
nearest_idx <- st_nearest_feature(model_points, grid_subset_centroids)

# Buat data frame model dengan psi_mean dan psi_sd
df_model <- as.data.frame(model_Kelempiau$coords) %>%
  mutate(psi_mean = psi_mean_last,
         psi_sd = psi_sd_last)

# Tambahkan hasil model ke subset grid berdasarkan nearest
grid_subset$psi_mean <- NA
grid_subset$psi_sd <- NA

grid_subset$psi_mean[nearest_idx] <- df_model$psi_mean
grid_subset$psi_sd[nearest_idx] <- df_model$psi_sd
# Sekarang grid_subset sudah mengandung hasil model

# Pastikan grid_subset sudah punya kolom psi_mean (dari langkah sebelumnya)

ggplot(grid_subset) +
  geom_sf(aes(fill = psi_mean), color = NA) +   # polygon diisi warna probabilitas
  scale_fill_viridis_c(option = "plasma", name = "Probabilitas\nOkupansi") +
  theme_minimal() +
  labs(title = "Peta Probabilitas Okupansi Kelempiau",
       subtitle = "Periode terakhir (waktu ke-49)",
       caption = "Data model dan grid subset") +
  theme(
    legend.position = "right",
    plot.title = element_text(size = 16, face = "bold"),
    plot.subtitle = element_text(size = 12)
  )

10. EXTRAPOLASI DARI MODEL BERBASIS GRID 2X2 KM > KOVARIAT LINGKUKANGN 30M

##### Langkah-langkah Skrip Ekstrapolasi:
# 1. Ambil mean beta dari hasil model
beta_mean <- colMeans(model_Kelempiau$beta.samples)

# 2. Raster kovariat (sudah di-scale dengan mean & sd dari training data!)
# Catatan: Jika belum diskalakan, gunakan nilai mean dan sd dari model training:
means <- attr(site_covs_scaled, "scaled:center")
sds <- attr(site_covs_scaled, "scaled:scale")

# Load ulang raster kovariat (pastikan urutannya sama dengan model)
rasters_var_env <- rast(c("output/Elevasi.tif", "output/Jarak Dari Jalan.tif",
                          "output/Jarak Dari Permukiman.tif", "output/Jarak Dari Sungai.tif",
                          "output/LAI.tif", "output/NDVI.tif", "output/Slope.tif",
                          "output/TUPLAH.tif"))

names(rasters_var_env) <- c("elev", "dist_jalan", "dist_permukiman", 
                            "dist_sungai", "lai", "ndvi", "slope", "tuplah")

# 3. Skala ulang semua raster menggunakan mean & sd dari model
scale_raster <- function(r, mean_val, sd_val) {
  (r - mean_val) / sd_val
}

rasters_scaled <- c()
for (i in 1:nlyr(rasters_var_env)) {
  layer_name <- names(rasters_var_env)[i]
  scaled_layer <- scale_raster(rasters_var_env[[i]], means[layer_name], sds[layer_name])
  names(scaled_layer) <- layer_name
  rasters_scaled <- c(rasters_scaled, scaled_layer)
}

# 4. Hitung prediksi ψ = inv.logit(Xβ)
# Pertama stack semua raster jadi 1
X_stack <- rasters_scaled

# Terapkan logit(ψ) = β0 + β1*x1 + β2*x2 + ...
logit_psi_raster <- beta_mean[1]  # intercept
for (i in 1:(length(beta_mean) - 1)) {
  logit_psi_raster <- logit_psi_raster + X_stack[[i]] * beta_mean[i + 1]
}

# 5. Transformasi ke skala probabilitas (logit−1)
inv_logit <- function(x) { 1 / (1 + exp(-x)) }
psi_raster <- app(logit_psi_raster, inv_logit)

# 6. Simpan hasil sebagai raster prediksi ψ
writeRaster(psi_raster, "output/prediksi_okupansi_Kelempiau_30m.tif", overwrite = TRUE)

# (Opsional) Plot
plot(psi_raster, main = "Prediksi Probabilitas Okupansi Kelempiau")

11. MODEL AVERAGING (Weighted Averaging / Ensemble)

# Tujuan: Menggabungkan prediksi dari dua model habitat suitability
# (Model Okupansi dan MaxEnt) menggunakan pendekatan Weighted Averaging
# berdasarkan nilai AUC masing-masing model.

# 1. BACA DATA RASTER DARI MODEL MAXENT
# Buka file data hasil model MaxEnt (.ascii)
maxent_raster <- rast("output/Percobaan_2_kelempiau_-_Hylobates_albibarbis/kelempiau_-_Hylobates_albibarbis_avg.asc")
r_min <- rast("output/Percobaan_2_kelempiau_-_Hylobates_albibarbis/kelempiau_-_Hylobates_albibarbis_min.asc")
r_stddev <- rast("output/Percobaan_2_kelempiau_-_Hylobates_albibarbis/kelempiau_-_Hylobates_albibarbis_stddev.asc")

par(mfrow = c(1, 1)) 
# Tampilkan peta dan statistik dasar
plot(maxent_raster, main = "Pictures of the MaxEnt model (avg)")

plot(r_min, main = "Pictures of the MaxEnt model (min)")

plot(r_stddev, main = "Pictures of the MaxEnt model (stddev)")

# 2. BACA DATA RASTER DARI MODEL OKUPANSI (stPGOcc) 
# Raster output dari prediksi rata-rata okupansi dari model `stPGOcc`
psi_raster <- rast("output/prediksi_okupansi_Kelempiau_30m.tif")


# 3. SESUAIKAN SISTEM KOORDINAT (CRS) JIKA BERBEDA 
# Pastikan kedua raster memiliki sistem proyeksi yang sama
# if (!compareCRS(psi_raster, maxent_raster)) { psi_raster <- projectRaster(psi_raster, crs = crs(maxent_raster))}

crs(maxent_raster)
## [1] ""
crs(psi_raster)
## [1] "PROJCRS[\"WGS 84 / UTM zone 49S\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 49S\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",111,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",10000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 108°E and 114°E, southern hemisphere between 80°S and equator, onshore and offshore. Australia. Indonesia.\"],\n        BBOX[-80,108,0,114]],\n    ID[\"EPSG\",32749]]"
# Assign CRS dari psi_raster ke maxent_raster
crs(maxent_raster) <- crs(psi_raster)

# Cek ulang
crs(maxent_raster)
## [1] "PROJCRS[\"WGS 84 / UTM zone 49S\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 49S\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",111,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",10000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 108°E and 114°E, southern hemisphere between 80°S and equator, onshore and offshore. Australia. Indonesia.\"],\n        BBOX[-80,108,0,114]],\n    ID[\"EPSG\",32749]]"
# Setelah ini, coba lagi projectRaster kalau perlu
if (!same.crs(maxent_raster, psi_raster)) {
  psi_raster <- project(psi_raster, maxent_raster)
}

# 4. RESAMPLE RASTER OKUPANSI UNTUK MENYESUAIKAN MAXENT
# Samakan resolusi, dimensi, dan extent raster okupansi agar cocok dengan MaxEnt
# Gunakan interpolasi bilinear karena nilai bersifat kontinu (probabilitas)
psi_raster_resampled <- resample(psi_raster, maxent_raster, method = "bilinear")


# 5. EKSTRAK NILAI DAN SKALAKAN KE 0–1 (NORMALISASI)
# Ekstrak nilai raster dari kedua model
maxent_vals <- values(maxent_raster)
psi_vals <- values(psi_raster_resampled)

# Identifikasi sel raster yang memiliki nilai valid pada kedua raster
valid_idx <- which(!is.na(maxent_vals) & !is.na(psi_vals))

# Buat vektor kosong dengan panjang sama, isi dengan NA
maxent_scaled <- rep(NA, length(maxent_vals))
psi_scaled <- rep(NA, length(psi_vals))

# Normalisasi nilai MaxEnt ke rentang 0–1 (hanya untuk nilai valid)
maxent_scaled[valid_idx] <- (maxent_vals[valid_idx] - min(maxent_vals[valid_idx])) /
  (max(maxent_vals[valid_idx]) - min(maxent_vals[valid_idx]))

# Normalisasi nilai okupansi ke rentang 0–1 (hanya untuk nilai valid)
psi_scaled[valid_idx] <- (psi_vals[valid_idx] - min(psi_vals[valid_idx])) /
  (max(psi_vals[valid_idx]) - min(psi_vals[valid_idx]))

# 6. HITUNG PREDIKSI ENSEMBLE BERDASARKAN BOBOT AUC

# Bobot berdasarkan Area Under Curve (AUC) dari validasi masing-masing model
auc_maxent <- 0.929
auc_occupancy <- 0.8210285

# Hitung bobot normalisasi
w_maxent <- auc_maxent / (auc_maxent + auc_occupancy)
w_occupancy <- auc_occupancy / (auc_maxent + auc_occupancy)

# Hitung prediksi ensemble sebagai rata-rata tertimbang
ensemble_vals <- rep(NA, length(maxent_vals))
ensemble_vals[valid_idx] <- (w_maxent * maxent_scaled[valid_idx]) +
  (w_occupancy * psi_scaled[valid_idx])

# 7. SIMPAN DAN VISUALISASI RASTER ENSEMBLE
# Buat objek raster baru dengan struktur yang sama seperti MaxEnt
ensemble_raster <- rast(maxent_raster)
values(ensemble_raster) <- ensemble_vals

# Simpan raster ensemble ke file (GeoTIFF)
writeRaster(ensemble_raster, "output/prediksi_ensemble_Kelempiau_30m.tif", overwrite = TRUE)

# Tampilkan peta ensemble hasil penggabungan dua model
plot(ensemble_raster, main = "Prediksi Ensemble Okupansi Kelempiau")

ensemble_raster
## class       : SpatRaster 
## size        : 2407, 2559, 1  (nrow, ncol, nlyr)
## resolution  : 30.00311, 30.00311  (x, y)
## extent      : 478628.2, 555406.2, 9893919, 9966136  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 49S (EPSG:32749) 
## source(s)   : memory
## varname     : kelempiau_-_Hylobates_albibarbis_avg 
## name        : kelempiau_-_Hylobates_albibarbis_avg 
## min value   :                         0.0001359526 
## max value   :                         0.7516768775
summary(ensemble_raster)
##  kelempiau_._Hylobates_albibarbis_avg
##  Min.   :0.000                       
##  1st Qu.:0.128                       
##  Median :0.410                       
##  Mean   :0.315                       
##  3rd Qu.:0.469                       
##  Max.   :0.721                       
##  NA's   :59159