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")
}

# Here the intersexes are excluded from the analysis

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)

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)
plot(tsdF.WA, resultmcmc = tsdF_mcmc.WA, replicate.CI = 2000, xlim=c(24, 36))

  
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)

plot(tsdF.EP, resultmcmc = tsdF_mcmc.EP, replicate.CI = 2000, xlim=c(24, 36))

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)
plot(tsdF.NEI, resultmcmc = tsdF_mcmc.NEI, replicate.CI = 2000, xlim=c(24, 36))

#########

py <- as.numeric(tsdF_mcmc.NEI$resultMCMC[[1]][, "SH"]+tsdF_mcmc.NEI$resultMCMC[[1]][, "SL"])
px <- as.numeric(tsdF_mcmc.NEI$resultMCMC[[1]][, "P"])
plot(x=px, y=py, pch=".", col="red", xlab="Pivotal temperature", ylab="TRT", las=1, bty="n", xlim=c(27, 32), cex=2)

py <- as.numeric(tsdF_mcmc.WA$resultMCMC[[1]][, "SH"]+tsdF_mcmc.WA$resultMCMC[[1]][, "SL"])
px <- as.numeric(tsdF_mcmc.WA$resultMCMC[[1]][, "P"])
points(x=px, y=py, pch=".", col="green", cex=2)

py <- as.numeric(tsdF_mcmc.EP$resultMCMC[[1]][, "SH"]+tsdF_mcmc.EP$resultMCMC[[1]][, "SL"])
px <- as.numeric(tsdF_mcmc.EP$resultMCMC[[1]][, "P"])
points(x=px, y=py, pch=".", col="blue", cex=2)

legend("topleft", legend=c("Northeastern Indian", "West Atlantic", "East Pacific"), pch=19, col=c("red", "green","blue"))

py <- as.numeric(tsdF_mcmc.NEI$resultMCMC[[1]][, "SH"])
px <- as.numeric(tsdF_mcmc.NEI$resultMCMC[[1]][, "SL"])
plot(x=px, y=py, pch=".", col="red", xlab="SL", ylab="SH", las=1, bty="n", xlim=c(0, 10), cex=2)

py <- as.numeric(tsdF_mcmc.WA$resultMCMC[[1]][, "SH"])
px <- as.numeric(tsdF_mcmc.WA$resultMCMC[[1]][, "SL"])
points(x=px, y=py, pch=".", col="green", cex=2)

py <- as.numeric(tsdF_mcmc.EP$resultMCMC[[1]][, "SH"])
px <- as.numeric(tsdF_mcmc.EP$resultMCMC[[1]][, "SL"])
points(x=px, y=py, pch=".", col="blue", cex=2)

legend("topright", legend=c("Northeastern Indian", "West Atlantic", "East Pacific"), pch=19, col=c("red", "green","blue"))

#########

# Aggregate the mcmc into a single file

dim(tsdF_mcmc.Global$WAIC)
dim(tsdF_mcmc.WA$WAIC)
dim(tsdF_mcmc.NEI$WAIC)
dim(tsdF_mcmc.EP$WAIC)

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]

#########
#Global model

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)

plot(tsdF.Global, resultmcmc = tsdF_mcmc.Global, replicate.CI = 2000, xlim=c(24, 36))


#########

library(loo)
tsdF_mcmc.Global_loo = loo(tsdF_mcmc.Global$WAIC)
tsdF_mcmc.RMUEffect_loo = loo(tsdF_mcmc.RMUEffect$WAIC)

loo_compare(list(global=tsdF_mcmc.Global_loo, 
                 RMU_Effect=tsdF_mcmc.RMUEffect_loo))

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