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)
data(nest)

# Data from Girondot, M. & Kaska, Y. 2014. A model to predict the thermal reaction norm for the embryo growth rate from field data. J. Therm. Biol. (IF 2.361), 45, 96-102.
head(nest)


formated <- FormatNests(data=nest, previous=NULL, col.Time="Time", 
                        hatchling.metric.mean=39.33, 
                        hatchling.metric.sd=1.92)

#######################
# Gompertz model of SCL
#######################

# See nest video to understand the origin of these parameters
pfixed <- c(rK=1.210977)
M0 = 0.3517074 

# These values are just a first set of parameters
x4 <- c('DHA' = 109.31113503282113, 'DHH' = 617.80695919563857, 
        'T12H' = 306.38890489505093, 'Rho25' = 229.37265815800225)

plotR(parameters = x4, curve="ML")

resultNest_4p_SSM <- searchR(parameters=x4, fixed.parameters=pfixed, 
                             temperatures=formated, 
                             integral=integral.Gompertz, M0=M0)

pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM, accept=TRUE)

# Take care, it can be very long; several days
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(result=resultNest_4p_SSM , 
                                       adaptive = TRUE          , 
                                       WAIC=TRUE                , 
                                       parametersMCMC=pMCMC     , 
                                       n.iter=10000             , 
                                       n.chains = 1             ,  
                                       n.adapt = 100            , 
                                       thin=10                  , 
                                       trace=100               )

save(resultNest_mcmc_4p_SSM, file="resultNest_mcmc_4p_SSM.Rdata")
save(resultNest_4p_SSM, file="resultNest_4p_SSM.Rdata")

plotR(resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
      curve = "MCMC quantiles", ylim=c(0, 8), show.hist = TRUE, 
      ylimH = c(0, 1), atH=c(0, 0.1, 0.2))

plot(x=resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
     xlim=c(0,70), ylimT=c(22, 32), GTRN.CI="MCMC", 
     replicate.CI = 100, las=1, 
     ylimS=c(0,45), series=1, 
     show.TSP = FALSE, show.third = FALSE, show.stages = FALSE, 
     embryo.stages="Caretta caretta.SCL")

plot(x=resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
     xlim=c(0,70), ylimT=c(22, 32), GTRN.CI="MCMC", 
     replicate.CI = 100, las=1, 
     ylimS=c(0,45), series=1, 
     embryo.stages="Caretta caretta.SCL")

temperatures30 <- GenerateConstInc(durations = 100000, temperatures = 30, names = "X30")
temperatures30_f <- FormatNests(temperatures30)

plot(x=resultNest_4p_SSM, resultmcmc = resultNest_mcmc_4p_SSM, 
     temperatures = temperatures30_f, 
     metric.end.incubation = "hatchling.metric",
     hatchling.metric = c(Mean=39.33, SD=1.92), 
     xlim=c(0,50), ylimT=c(22, 32), GTRN.CI="MCMC", 
     stop.at.hatchling.metric = TRUE, 
     replicate.CI = 100, las=1, 
     ylimS=c(0,40), series=1, 
     show.TSP = FALSE, show.third = FALSE, show.stages = FALSE, 
     embryo.stages="Generic.ProportionDevelopment")

inests<- info.nests(x=resultNest_4p_SSM, 
                    resultmcmc = resultNest_mcmc_4p_SSM, 
                    GTRN.CI="MCMC", 
                    replicate.CI = 100, 
                    out="summary", 
                    embryo.stages="Caretta caretta.SCL")
inests <- inests$summary[, !is.na(inests$summary[1, ])]

plot_errbar(x=inests$TSP.TimeWeighted.temperature.quantile_0.5, 
            y=inests$MiddleThird.TimeWeighted.temperature.quantile_0.5, 
            x.plus = inests$TSP.TimeWeighted.temperature.quantile_0.975, 
            y.plus = inests$MiddleThird.TimeWeighted.temperature.quantile_0.975,
            x.minus = inests$TSP.TimeWeighted.temperature.quantile_0.025, 
            y.minus = inests$MiddleThird.TimeWeighted.temperature.quantile_0.025,
            xlab = "Average T in °C - TSP - Gompertz", 
            ylab = "Average T in °C - Middlethird - Gompertz", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for Middlethird = ", 
                   specify_decimal(x= mean(inests$MiddleThird.TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests$MiddleThird.TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(inests$TSP.TimeWeighted.temperature.quantile_0.5, 
            inests$TimeWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.TimeWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.TimeWeighted.temperature.quantile_0.975, 
            y.minus = inests$TimeWeighted.temperature.quantile_0.025,
            y.plus = inests$TimeWeighted.temperature.quantile_0.975,
            xlab = "Growth-W T in °C - TSP - Gompertz", 
            ylab = "Growth-W T in °C - Total - Gompertz", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for total = ", 
                   specify_decimal(x= mean(inests$TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests$TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(inests$TSP.GrowthWeighted.temperature.quantile_0.5, 
            inests$MiddleThird.GrowthWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests$MiddleThird.GrowthWeighted.temperature.quantile_0.025,
            y.plus = inests$MiddleThird.GrowthWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Gompertz", 
            ylab = "Growth-W T in °C - Middlethird - Gompertz", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for Middlethird = ", 
                   specify_decimal(x= mean(inests$MiddleThird.GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests$MiddleThird.GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(x=inests$TSP.GrowthWeighted.temperature.quantile_0.5, 
            y=inests$GrowthWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests$GrowthWeighted.temperature.quantile_0.025,
            y.plus = inests$GrowthWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Gompertz", 
            ylab = "Growth-W T in °C - Total - Gompertz", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for total = ", 
                   specify_decimal(x= mean(inests$GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests$GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(x=inests$TSP.GrowthWeighted.temperature.quantile_0.5, 
            y=inests$TSP.TimeWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests$TSP.TimeWeighted.temperature.quantile_0.025,
            y.plus = inests$TSP.TimeWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Gompertz", 
            ylab = "Average T in °C - Total - Gompertz", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for average = ", 
                   specify_decimal(x= mean(inests$TSP.TimeWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests$TSP.TimeWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

####################################
# linear model of progress of growth
####################################

pfixed <- NULL
M0 = 0

x4 <- c('DHA' = 64.868697530424186, 'DHH' = 673.18292743646771, 
        'T12H' = 400.90952554047749, 'Rho25' = 82.217237723502123)

resultNest_4p_SSM_Linear <- searchR(parameters=x4, fixed.parameters=pfixed, 
                                    temperatures=formated, integral=integral.linear, M0=M0, 
                                    hatchling.metric=c(Mean=39.33, SD=1.92)/39.33)

plotR(resultNest_4p_SSM_Linear, ylim=c(0, 2), scaleY= 100000, curve = "ML")

pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM_Linear, accept=TRUE)

# Take care, it can be very long; several days on old computer
resultNest_mcmc_4p_SSM_Linear <- GRTRN_MHmcmc(result=resultNest_4p_SSM , 
                                              adaptive = TRUE          , 
                                              WAIC=TRUE                , 
                                              parametersMCMC=pMCMC     , 
                                              n.iter=10000             , 
                                              n.chains = 1             ,  
                                              n.adapt = 100            , 
                                              thin=10                  , 
                                              trace=100               )

plotR(resultNest_4p_SSM_Linear, resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
      curve = "MCMC quantiles", ylim=c(0, 2), show.hist = TRUE, 
      ylimH = c(0, 1), atH=c(0, 0.1, 0.2))

plot(x=resultNest_4p_SSM_Linear, resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
     xlim=c(0,70), ylimT=c(22, 32), GTRN.CI="MCMC", 
     replicate.CI = 100, las=1, 
     ylimS=c(0,1), series=1, 
     show.TSP = TRUE, show.third = TRUE, show.stages = TRUE, 
     embryo.stages="Generic.ProportionDevelopment")


temperatures30 <- GenerateConstInc(durations = 100000, temperatures = 30, names = "X30")
temperatures30_f <- FormatNests(temperatures30)

plot(x=resultNest_4p_SSM_Linear, resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
     temperatures = temperatures30_f, 
     metric.end.incubation = "hatchling.metric",
     hatchling.metric = c(Mean=39.33, SD=1.92)/39.33, 
     xlim=c(0,50), ylimT=c(22, 32), GTRN.CI="MCMC", 
     stop.at.hatchling.metric = TRUE, 
     replicate.CI = 100, las=1, 
     ylimS=c(0,1), series=1, 
     show.TSP = TRUE, show.third = TRUE, show.stages = TRUE, 
     embryo.stages="Generic.ProportionDevelopment")


inests_linear <- info.nests(x=resultNest_4p_SSM_Linear, 
                            resultmcmc = resultNest_mcmc_4p_SSM_Linear, 
                            GTRN.CI="MCMC", 
                            replicate.CI = 100, 
                            out="summary", 
                            embryo.stages="Generic.ProportionDevelopment")
inests_linear <- inests_linear$summary[, !is.na(inests_linear$summary[1, ])]

plot_errbar(x=inests_linear$TSP.TimeWeighted.temperature.quantile_0.5, 
            y=inests_linear$MiddleThird.TimeWeighted.temperature.quantile_0.5, 
            x.plus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.975, 
            y.plus = inests_linear$MiddleThird.TimeWeighted.temperature.quantile_0.975,
            x.minus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.025, 
            y.minus = inests_linear$MiddleThird.TimeWeighted.temperature.quantile_0.025,
            xlab = "Average T in °C - TSP - Linear", 
            ylab = "Average T in °C - Middlethird - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for Middlethird = ", 
                   specify_decimal(x= mean(inests_linear$MiddleThird.TimeWeighted.temperature.mean-inests_linear$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$MiddleThird.TimeWeighted.temperature.mean-inests_linear$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(inests_linear$TSP.TimeWeighted.temperature.quantile_0.5, 
            inests_linear$TimeWeighted.temperature.quantile_0.5, 
            x.minus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.025, 
            x.plus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$TimeWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$TimeWeighted.temperature.quantile_0.975,
            xlab = "Growth-W T in °C - TSP - Linear", 
            ylab = "Growth-W T in °C - Total - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for total = ", 
                   specify_decimal(x= mean(inests_linear$TimeWeighted.temperature.mean-inests_linear$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$TimeWeighted.temperature.mean-inests_linear$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(inests_linear$TSP.GrowthWeighted.temperature.quantile_0.5, 
            inests_linear$MiddleThird.GrowthWeighted.temperature.quantile_0.5, 
            x.minus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$MiddleThird.GrowthWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$MiddleThird.GrowthWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Linear", 
            ylab = "Growth-W T in °C - Middlethird - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for Middlethird = ", 
                   specify_decimal(x= mean(inests_linear$MiddleThird.GrowthWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$MiddleThird.GrowthWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(x=inests_linear$TSP.GrowthWeighted.temperature.quantile_0.5, 
            y=inests_linear$GrowthWeighted.temperature.quantile_0.5, 
            x.minus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$GrowthWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$GrowthWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Linear", 
            ylab = "Growth-W T in °C - Total - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for total = ", 
                   specify_decimal(x= mean(inests_linear$GrowthWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$GrowthWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)

plot_errbar(x=inests_linear$TSP.GrowthWeighted.temperature.quantile_0.5, 
            y=inests_linear$TSP.TimeWeighted.temperature.quantile_0.5, 
            x.minus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Linear", 
            ylab = "Average T in °C - Total - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for average = ", 
                   specify_decimal(x= mean(inests_linear$TSP.TimeWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.TimeWeighted.temperature.mean-inests_linear$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)


###################################
# Comparison of TSP for both models
###################################

plot_errbar(inests$TSP.TimeWeighted.temperature.quantile_0.5, 
            inests_linear$TSP.TimeWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.TimeWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.TimeWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$TSP.TimeWeighted.temperature.quantile_0.975,
            
            xlab = "Average temperature within TSP for Gompertz model in °C", 
            ylab = "Average temperature within TSP for linear model in °C", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for linear = ", 
                   specify_decimal(x= mean(inests_linear$TSP.TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.TimeWeighted.temperature.mean-inests$TSP.TimeWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)


plot_errbar(inests$TSP.GrowthWeighted.temperature.quantile_0.5, 
            inests_linear$TSP.GrowthWeighted.temperature.quantile_0.5, 
            x.minus = inests$TSP.GrowthWeighted.temperature.quantile_0.025, 
            x.plus = inests$TSP.GrowthWeighted.temperature.quantile_0.975, 
            y.minus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.025,
            y.plus = inests_linear$TSP.GrowthWeighted.temperature.quantile_0.975,
            
            xlab = "Growth-W T in °C - TSP - Gompertz", 
            ylab = "Growth-W T in °C - TSP - Linear", las=1, bty="n", 
            xlim = c(27, 32), ylim = c(27, 32))
segments(x0=27, y0=27, x1=32, y1=32, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for linear = ", 
                   specify_decimal(x= mean(inests_linear$TSP.GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = "."), " °C (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.GrowthWeighted.temperature.mean-inests$TSP.GrowthWeighted.temperature.mean), 
                                   decimals = 2, decimal.point = ".")    
                   ," °C)"), 
     pos=4)


plot_errbar(inests$TSP.begin.quantile_0.5/24/60, 
            inests_linear$TSP.begin.quantile_0.5/24/60, 
            x.minus = inests$TSP.begin.quantile_0.025/24/60, 
            x.plus = inests$TSP.begin.quantile_0.975/24/60, 
            y.minus = inests_linear$TSP.begin.quantile_0.025/24/60,
            y.plus = inests_linear$TSP.begin.quantile_0.975/24/60,
            
            xlab = "Begin TSP for Gompertz model in days", 
            ylab = "Begin TSP for linear model in days", las=1, bty="n", 
            xlim = c(15, 33), ylim = c(15, 33)
)
segments(x0=15, y0=15, x1=33, y1=33, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for linear = ", 
                   specify_decimal(x= mean(inests_linear$TSP.begin.mean/24/60-inests$TSP.begin.mean/24/60), 
                                   decimals = 2, decimal.point = "."), " days (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.begin.mean/24/60-inests$TSP.begin.mean/24/60), 
                                   decimals = 2, decimal.point = ".")    
                   ," days)"), 
     pos=4)


plot_errbar(inests$TSP.end.quantile_0.5/24/60, 
            inests_linear$TSP.end.quantile_0.5/24/60, 
            x.minus = inests$TSP.end.quantile_0.025/24/60, 
            x.plus = inests$TSP.end.quantile_0.975/24/60, 
            y.minus = inests_linear$TSP.end.quantile_0.025/24/60,
            y.plus = inests_linear$TSP.end.quantile_0.975/24/60,
            
            xlab = "End TSP for Gompertz model in days", 
            ylab = "End TSP for linear model in days", las=1, bty="n", 
            xlim = c(30, 55), ylim = c(30, 55)
)
segments(x0=30, y0=30, x1=50, y1=55, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for linear = ", 
                   specify_decimal(x= mean(inests_linear$TSP.end.mean/24/60-inests$TSP.end.mean/24/60), 
                                   decimals = 2, decimal.point = "."), " days (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.end.mean/24/60-inests$TSP.end.mean/24/60), 
                                   decimals = 2, decimal.point = ".")    
                   ," days)"), 
     pos=4)

plot_errbar(inests$TSP.length.mean/24/60, 
            inests_linear$TSP.length.mean/24/60, 
            x.minus = inests$TSP.length.quantile_0.025/24/60, 
            x.plus = inests$TSP.length.quantile_0.975/24/60, 
            y.minus = inests_linear$TSP.length.quantile_0.025/24/60,
            y.plus = inests_linear$TSP.length.quantile_0.975/24/60,
            
            xlab = "TSP length for Gompertz model in days", 
            ylab = "TSP length for linear model in days", las=1, bty="n", 
            xlim = c(14, 24), ylim = c(14, 24)
)
segments(x0=14, y0=14, x1=24, y1=24, lty=2, col="red")

text(x=ScalePreviousPlot(x=0.1, y=0.7)$x, 
     y=ScalePreviousPlot(x=0.1, y=0.7)$y, 
     labels=paste0("Bias for linear = ", 
                   specify_decimal(x= mean(inests_linear$TSP.length.mean/24/60-inests$TSP.length.mean/24/60), 
                                   decimals = 2, decimal.point = "."), " days (SD ",  
                   specify_decimal(x= sd(inests_linear$TSP.length.mean/24/60-inests$TSP.length.mean/24/60), 
                                   decimals = 2, decimal.point = ".")    
                   ," days)"), 
     pos=4)