# PREPARATION STEPS (NOT FOR PUBLICATION)

# Create data frame for export from ABCDR workflow
# df <- cbind(df2[, c('nett', 'AI', 'h')], y1, y2, y3, y4)
folder <- "C:/Users/thoma.TOM2025/Documents/Current Projects/Moore/ABCDR/figures/supplementary"
# write.csv(df, file = paste0(folder, '/Mal_J supplement.csv'))

# RUN STAN
library(rstan)
library(loo)
options(mc.cores = parallel::detectCores())

stan_results <- function(model, data_list = data_list) {
  fit <- rstan::stan(file = paste0(model, ".stan"),
                   model_name = model,
                   control = list(max_treedepth = 20),
                   data = data_list)
  # Use of loo package to compare fit of stan models
  log_lik_1 <- loo::extract_log_lik(fit, merge_chains = FALSE)
  r_eff <- loo::relative_eff(exp(log_lik_1), cores = getOption("mc.cores", 1))
  loo <- loo::loo(log_lik_1, r_eff = r_eff, cores = getOption("mc.cores", 1))
  return(list(fit = fit, loo = loo))
}

setwd(folder)

# DEFAULT STAN ANALYSIS
df <- read.csv('ABCDR data.csv')
y <- with(df, as.matrix(cbind(y1, y2, y3, y4)))
data_list <- list(N = nrow(df), nett = df$nett, chem = df$AI, y = y, h = df$h)
abcdr_results <- stan_results('ABCDR', data_list)
summary(abcdr_results$fit)$summary[1:8, c(4, 6, 8)]
abcdr_results$loo$estimates

# logistic model
library(dplyr)
library(tidyr)
dfb <- df %>% mutate(t1 = y1 + y2 + y3 + y4,
              t2 = y2 + y3 + y4,
              t3 = y3 + y4,
              t4 = y4) %>%
  pivot_longer(cols = starts_with("t"), names_to = "stept", values_to = "total") %>%
  mutate(step = recode(stept, `t1` = 1, `t2` = 2, 't3' = 3, 't4' = 4)) %>%
  filter(step < 4) %>%
  mutate(pos = case_when(step == 1 ~ y1, step == 2 ~ y2, step == 3 ~ y3))
data_list <- list(N = nrow(dfb),
                  step = dfb$step,
                  total = dfb$total,
                  pos = dfb$pos,
                  holed_area = dfb$h * 172950,
                  nett = dfb$nett,
                  chem = dfb$AI)
logistic_results <- stan_results('logisticABCDR', data_list)
summary(logistic_results$fit)$summary[1:13, c(4, 6, 8)]
logistic_results$loo$estimates