if (FALSE) {
  # Run this part of the code only once
  install.packages("embryogrowth")
  
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/HelpersMG.tar.gz", repos=NULL, type="source")
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/embryogrowth.tar.gz", repos=NULL, type="source")
}

library(embryogrowth)
Lepidochelys.olivacea = subset(DatabaseTSD, Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0))
colnames(Lepidochelys.olivacea)
unique(Lepidochelys.olivacea$RMU.2023)
nrow(Lepidochelys.olivacea)

# Fixed effect

Lepidochelys.olivacea.WA = subset(Lepidochelys.olivacea, subset=RMU.2023 == "West Atlantic")
nrow(Lepidochelys.olivacea.WA)

tsdF.WA = tsd(df=Lepidochelys.olivacea.WA, equation = "flexit**", 
              parameters.initial = c(P=30, SH=1, SL=1))

priorsN = tsd_MHmcmc_p(result = tsdF.WA, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF.WA, default = "dunif", accept = TRUE)
priors.WA = rbind(priorsN[1, , drop=FALSE], priorsU[2:3, ])
priors.WA[1, "Prior1"] = 30
priors.WA[2:3, c("Prior1", "Min")] = 0
priors.WA[2:3, c("Prior2", "Max")] = 10

tsdF_mcmc.WA = tsd_MHmcmc(result = tsdF.WA, parametersMCMC = priors.WA, 
                          n.iter = 20000, thin = 10, n.adapt=1000)


Lepidochelys.olivacea.EP = subset(Lepidochelys.olivacea, subset=RMU.2023 == "East Pacific")
nrow(Lepidochelys.olivacea.EP)

tsdF.EP = tsd(df=Lepidochelys.olivacea.EP, equation = "flexit**", 
              parameters.initial = c(P=30, SH=1, SL=1))
priorsN = tsd_MHmcmc_p(result = tsdF.EP, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF.EP, default = "dunif", accept = TRUE)
priors.EP = rbind(priorsN[1, , drop=FALSE], priorsU[2:3, ])
priors.EP[1, "Prior1"] = 30
priors.EP[2:3, c("Prior1", "Min")] = 0
priors.EP[2:3, c("Prior2", "Max")] = 10

tsdF_mcmc.EP = tsd_MHmcmc(result = tsdF.EP, parametersMCMC = priors.EP, 
                          n.iter = 20000, thin = 10, n.adapt=1000)


Lepidochelys.olivacea.NEI = subset(Lepidochelys.olivacea, subset=RMU.2023 == "Northeast Indian")
nrow(Lepidochelys.olivacea.NEI)

tsdF.NEI = tsd(df=Lepidochelys.olivacea.NEI, equation = "flexit**", 
               parameters.initial = c(P=30, SH=1, SL=1))

priorsN = tsd_MHmcmc_p(result = tsdF.NEI, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF.NEI, default = "dunif", accept = TRUE)
priors.NEI = rbind(priorsN[1, , drop=FALSE], priorsU[2:3, ])
priors.NEI[1, "Prior1"] = 30
priors.NEI[2:3, c("Prior1", "Min")] = 0
priors.NEI[2:3, c("Prior2", "Max")] = 10

tsdF_mcmc.NEI = tsd_MHmcmc(result = tsdF.NEI, parametersMCMC = priors.NEI, 
                           n.iter = 20000, thin = 10, n.adapt=1000)

## Generate object mcmc with RMU as fixed effect

# Aggregate the mcmc into a single file
tsdF_mcmc.RMUEffect <- tsdF_mcmc.Global
tsdF_mcmc.RMUEffect$WAIC[, , 1:12] <- tsdF_mcmc.WA$WAIC[, , 1:12]
tsdF_mcmc.RMUEffect$WAIC[, , 12+1:6] <- tsdF_mcmc.NEI$WAIC[, , 1:6]
tsdF_mcmc.RMUEffect$WAIC[, , 18+1:22] <- tsdF_mcmc.EP$WAIC[, , 1:22]

# No RMU Effect

tsdF.Global <- tsd(df=Lepidochelys.olivacea, equation = "flexit**", 
                  parameters.initial = c(P=30, SH=1, SL=1))
plot(tsdF.Global)

priorsN = tsd_MHmcmc_p(result = tsdF.Global, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF.Global, default = "dunif", accept = TRUE)
priors.Global = rbind(priorsN[1, , drop=FALSE], priorsU[2:3, ])
priors.Global[1, "Prior1"] = 30
priors.Global[2:3, c("Prior1", "Min")] = 0
priors.Global[2:3, c("Prior2", "Max")] = 10
priors.Global

tsdF_mcmc.Global = tsd_MHmcmc(result = tsdF.Global, parametersMCMC = priors.Global, 
                              n.iter = 20000, thin = 10, n.adapt=1000)

# RMU as random effect

library(brms)

prior_ini = set_prior(prior = "normal(30,5)", class = "b", coef="", 
                       nlpar = "P", lb = 27, ub = 35, check = TRUE) + 
             set_prior(prior = "uniform(0, 12)", class = "b", coef="", 
                       nlpar = "SL", lb = 0, ub = 12, check = TRUE) + 
             set_prior(prior = "uniform(0, 12)", class = "b", coef="", 
                       nlpar = "SH", lb = 0, ub = 12, check = TRUE)

options(mc.cores = 1)

control = list(
  adapt_engaged = TRUE,
  adapt_delta = 0.99, #increased from default of 0.8
  stepsize = 0.05, # 0.05 default
  max_treedepth = 20
)

flexit_RMU = brm(bf(Males | trials(Sexed) ~  
                                 1/(
                                   1+exp(
                                     ((-2.94443897916644*(P-Incubation.temperature.set))
                                      /
                                        (
                                          (SL/ (1+ exp(100*(Incubation.temperature.set - P))))+
                                            (SH/ (1+ exp(100*(P - Incubation.temperature.set))))
                                        )
                                     )
                                   )
                                 ), 
                               P ~  1  + (1 | RMU.2023), 
                               SL ~ 1 + (1 | RMU.2023), 
                               SH ~ 1 + (1 | RMU.2023),
                               nl = TRUE),
                            family = binomial(link="identity"), 
                            data = Lepidochelys.olivacea, 
                            save_pars = save_pars(all = TRUE), 
                            control = control, 
                            sample_prior = TRUE, 
                            prior = prior_ini,
                            warmup = 5000, 
                            thin = 10, 
                            iter = 100000, 
                            chains = 1, 
                            cores=1,
                            seed = 123, 
                            init=rep(list(list(b_P_Intercept=31, b_SL_Intercept=1.4, b_SH_Intercept=1.4)), 1))

plot(density(prior_draws(flexit_RMU, variable = "b_P_Intercept")[, 1]))
plot(density(prior_draws(flexit_RMU, variable = "b_SL_Intercept")[, 1]))
plot(density(prior_draws(flexit_RMU, variable = "b_SH_Intercept")[, 1]))

plot(mcmc_plot(flexit_RMU, type = "areas", prob = 0.95, variable = "b_P_Intercept"))
plot(mcmc_plot(flexit_RMU, type = "areas", prob = 0.95, variable = c("b_SL_Intercept", "b_SH_Intercept")))

plot(flexit_RMU, nvariables = 3, combo = c("hist", "trace"))

conditions = make_conditions(flexit_RMU, vars="RMU.2023")
k = conditional_effects(flexit_RMU, categorical = FALSE, 
                        effects='Incubation.temperature.set', 
                        condition=conditions, robust = TRUE, mean=TRUE, 
                        prob = 0.95,
                        re_formula=NULL, spaghetti = FALSE, ndraws = 1500)
plot(k)

# Estimation of WAIC

flexit_RMU_RE = loo(flexit_RMU)
tsdF_mcmc.Global_loo <- loo(tsdF_mcmc.Global$WAIC)
tsdF_mcmc.RMUEffect_loo <- loo(tsdF_mcmc.RMUEffect$WAIC)

# Comparison

loo_compare(x=list(Global=tsdF_mcmc.Global_loo, 
                   RMU.effect=tsdF_mcmc.RMUEffect_loo, 
                   RMU.Effect.RE=flexit_RMU_RE))

ic = c(global=tsdF_mcmc.Global_loo$estimates["looic", "Estimate"], 
        RMU.Effect=tsdF_mcmc.RMUEffect_loo$estimates["looic", "Estimate"], 
        RMU.Effect.RE=flexit_RMU_RE$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))

loo_compare(x=list(RMU.effect.FE=tsdF_mcmc.RMUEffect_loo, 
                   RMU.Effect.RE=flexit_RMU_RE))
ic = c(RMU.Effect=tsdF_mcmc.RMUEffect_loo$estimates["looic", "Estimate"], 
        RMU.Effect.RE=flexit_RMU_RE$estimates["looic", "Estimate"])
exp(ic) / sum(exp(ic))