--- title: "Recent advances on the estimation of the thermal reaction norm for sex ratios in sea turtle nest incubations using olive ridley data as an example" author: "F. Alberto Abreu-Grobois, B. Alejandra Morales-Mérida, Catherine Hart, Jean-Michel Guillon, Matthew H. Godfrey, Erik Navarro, Marc Girondot" date: "`r format(Sys.time(), '%d %B %Y')`" output: word_document: reference_docx: template.docx chunk_output_type: console editor_options: chunk_output_type: console --- ```{r setup, include=FALSE, message=FALSE} knitr::opts_chunk$set(fig.path='Figs/Dc_', # dev='postscript', dev='png', concordance=TRUE, echo=TRUE, dpi=300) ``` **Worked out examples of the publication data using package _embryogrowth_** As a support to the publication, we use the data for _Lepidochelys olivacea_ mentioned to illustrate the latest methods for the estimation of sex ratios from egg incubations at constant temperatures. Interest on the thermal reaction norm for sex ratio at constant incubation temperatures for conservation is not trivial as the variation in TSD patterns can have profound consequence for conservation: Carter, A.L., Bodensteiner, B.L., Iverson, J.B., Milne‐Zelman, C.L., Mitchell, T.S., Refsnider, J.M., Warner, D.A., Janzen, F.J., White, C., 2019. Breadth of the thermal response captures individual and geographic variation in temperature‐dependent sex determination. Functional Ecology. Hulin, V., Delmas, V., Girondot, M., Godfrey, M.H., Guillon, J.-M., 2009. Temperature-dependent sex determination and global change: Are some species at greater risk? Oecologia 160, 493-506. *Note: this analysis incorporates new data from Mexico not included in the package distribution at the time of publication.* Version >7.6 of embryogrowth package must be installed, or from CRAN, or from http://max2.ese.u-psud.fr/epc/conservation/CRAN/ site. ```{r eval=FALSE} ToInstall <- FALSE if (!("embryogrowth" %in% installed.packages())) { ToInstall <- TRUE } else { x <- packageVersion("embryogrowth") if (x < "7.6") { ToInstall <- TRUE } } } if (ToInstall) { install.packages("http:///max2.ese.u-psud.fr/epc/conservation/CRAN/HelpersMG.tar.gz", repos=NULL, type="source") install.packages("http:///max2.ese.u-psud.fr/epc/conservation/CRAN/embryogrowth.tar.gz", repos=NULL, type="source") } ``` ```{r} library("embryogrowth") library("openxlsx") library("maps") library("mapdata") library("sp") library("car") DatabaseTSD$Version[1] packageVersion("embryogrowth") db <- subset(DatabaseTSD, subset=(Species=="Lepidochelys olivacea") & (!is.na(Sexed)) & (Sexed!=0), select=c("Area", "RMU", "Incubation.temperature", "Incubation.temperature.SD", "Incubation.temperature.Amplitude", "2ndThird.Incubation.temperature.Amplitude", "Sexed", "Males", "Females", "Intersexes", "Reference")) # write.xlsx(db, file="Table1.xlsx") unique(as.character(subset(DatabaseTSD, Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0) & (is.na(Incubation.temperature.Amplitude) | (Incubation.temperature.Amplitude<1)), select="RMU")[, 1])) ``` Only data with amplitude of temperatures during incubation lower than 2 °C are used. Indeed, incubating with large temperature variations is known to provoke higher feminization than expected (Georges, 1989). As a consequence, these data will not be used: Georges, A., 1989. Female turtles from hot nests: Is it duration of incubation or proportion of development at high temperatures that matters? Oecologia 81, 323-328. Sandoval Espinoza (2012) indicates the use of data from incubations at different varying temperatures (pages 53-55) but the sex ratio data are not available in the thesis. Sandoval Espinoza, S., 2012. Proporción sexual en crías de tortuga *Lepidochelys olivacea* en corrales de incubación del Pacífico mexicano, In Centro Interdisciplinario de Ciencias Marinas. p. 158. Instituto Politécnico Nacional, La Paz, BCS, México. **Data come from:** Castheloge, V.D., de D. dos Santos, M.R., de Castilhos, J.C., de J. Filho, P.R., de C. Gomes, L., Clemente-Carvalho, R.B.G., Ferreira, P.D., 2018. Pivotal temperature and hatchling sex ratio of olive ridley sea turtles *Lepidochelys olivacea* from the South Atlantic coast of Brazil. Herpetological Conservation and Biology 13, 488-496. McCoy, C.J., Vogt, R.C., Censky, E.J., 1983. Temperature-controlled sex determination in the sea turtle *Lepidochelys olivacea*. J. Herpetol. 17, 404-406. Wibbels, T., Rostal, D.C., Byles, R., 1998. High pivotal temperature in the sex determination of the olive ridley sea turtle, *Lepidochelys olivacea*, from Playa Nancite, Costa Rica. Copeia 1998, 1086-1088. Dimond, M.T., 1985. Some effects of temperature on turtle egg incubation, In Recent Advances in Developmental Biology of Animals. eds S.C. Goel, C.B.L. Srivastava, pp. 35-39. Indian Society of Developmental Biologists, Poona, India. López Correa, J.Y., 2010. Diferenciación gonádica en crias de *Lepidochelys olivacea* (Eschscholtz, 1829) (Testudinata: Cheloniidae), p. 108. Instituto Politécnico Nacional, Centro Interdisciplinaria de Ciencias Marinas, La Paz, BCS. Mohanty-Hejmadi, P., Behra, M., Dimond, M.T., 1985. Temperature dependent sex differentiation in the olive ridley *Lepidochelys olivacea* and its implications for conservation, In Symposium on Endangered Marine Animals and Marine Parks. pp. 1-5. Marine Biological Association of India, Cochin. Merchant-Larios, H., Ruiz-Ramirez, S., Moreno-Mendoza, N., Marmolejo-Valencia, A., 1997. Correlation among thermosensitive period, estradiol response, and gonad differentiation in the sea turtle *Lepidochelys olivacea*. General and Comparative Endocrinology 107, 373-385. Merchant-Larios, H., Villalpando-Fierro, I., Centeno-Urruiza, B., 1989. Gonadal morphogenesis under controlled temperature in the sea turtle *Lepidochelys olivacea*. Herpetol. Monographs 3, 43-61. Navarro Sánchez, E.J., 2015. Efecto de la temperatura de incubación y la diferenciación sexual sobre la morfología de crías de tortuga marina *Lepidochelys olivacea*, Ciencias del mar y limnología. Universidad Nacional Autónoma de México, México, D.F., p. 89. # Sex ratios from field studies These data cannot be used for TSD pattern because they come from field studies with large temperature variation: Maulany, R. I., D. T. Booth, and G. S. Baxter. 2012. Emergence Success and Sex Ratio of Natural and Relocated Nests of Olive Ridley Turtles from Alas Purwo National Park, East Java, Indonesia. Copeia :738-747. (season 2009 and 2010) Hernandez-Echeagaray, O. E., R. Hernandez-Cornejo, M. Harfush-Melendez, and A. Garcia-Gasca. 2012. Evaluation of Sex Ratios of the Olive Ridley Sea Turtle (*Lepidochelys olivacea*) on the Arribada Nesting Beach, La Escobilla, Mexico. Marine Turtle Newsletter 133:12-16. (Season 2009-1010) Garcia, A., G. Ceballos, and R. Adaya. 2003. Intensive beach management as an improved sea turtle conservation strategy in Mexico. Biological Conservation 111:253-261. (Playa Cuixmala, Jalisco, Mexico, season 1994) Sandoval Espinoza, S., 2008. Pronóstico de la temperatura de los nidos de tortuga golfina (Lepidochelys olivacea) en función de la temperatura ambiente, la profundidad y el calor metabólico, In Centro Interdisciplinario De Ciencias Marinas. p. 84. Instituto Politécnico Nacional, La Paz, Mexico. Sandoval Espinoza, S., 2012. Proporción sexual en crías de tortuga Lepidochelys olivacea en corrales de incubación del Pacífico mexicano, In Centro Interdisciplinario de Ciencias Marinas. p. 158. Instituto Politécnico Nacional, La Paz, BCS, México. # Load data ```{r} Lo_PacificE_Mexico <- subset(DatabaseTSD, RMU=="Pacific, E" & Country =="Mexico" & Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Sexed", "Incubation.temperature", "Reference", "Longitude", "Latitude", "Area")) Lo_PacificE_CostaRica <- subset(DatabaseTSD, RMU=="Pacific, E" & Country =="Costa Rica" & Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Sexed", "Incubation.temperature", "Reference", "Longitude", "Latitude", "Area")) Lo_PacificE <- subset(DatabaseTSD, RMU=="Pacific, E" & Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Sexed", "Incubation.temperature", "Reference", "Longitude", "Latitude", "Area")) Lo_AtlanticWest <- subset(DatabaseTSD, RMU=="Atlantic, West" & Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Sexed", "Incubation.temperature", "Reference", "Longitude", "Latitude", "Area")) Lo_IndianNE <- subset(DatabaseTSD, RMU=="Indian, NE" & Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Sexed", "Incubation.temperature")) Lo_Global <- subset(DatabaseTSD, Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0), select=c("Males", "Females", "Intersexes", "Sexed", "Incubation.temperature", "Reference", "Longitude", "Latitude", "Area")) ``` The total number of sexed embryos is `r sum(Lo_Global$Males, na.rm=TRUE)+sum(Lo_Global$Females, na.rm=TRUE)+sum(Lo_Global$Intersexes, na.rm=TRUE)` with `r sum(Lo_Global$Males, na.rm=TRUE)` males, `r sum(Lo_Global$Females, na.rm=TRUE)` females and `r sum(Lo_Global$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_PacificE)`. For Costa Rica, the total number of sexed embryos is `r sum(Lo_PacificE_CostaRica$Males, na.rm=TRUE)+sum(Lo_PacificE_CostaRica$Females, na.rm=TRUE)+sum(Lo_PacificE_CostaRica$Intersexes, na.rm=TRUE)` with `r sum(Lo_PacificE_CostaRica$Males, na.rm=TRUE)` males, `r sum(Lo_PacificE_CostaRica$Females, na.rm=TRUE)` females and `r sum(Lo_PacificE_CostaRica$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_PacificE_CostaRica)`. For Mexico, the total number of sexed embryos is `r sum(Lo_PacificE_Mexico$Males, na.rm=TRUE)+sum(Lo_PacificE_Mexico$Females, na.rm=TRUE)+sum(Lo_PacificE_Mexico$Intersexes, na.rm=TRUE)` with `r sum(Lo_PacificE_Mexico$Males, na.rm=TRUE)` males, `r sum(Lo_PacificE_Mexico$Females, na.rm=TRUE)` females and `r sum(Lo_PacificE_Mexico$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_PacificE_Mexico)`. At the scale of East Pacific RMU, the total number of sexed embryos is `r sum(Lo_PacificE$Males, na.rm=TRUE)+sum(Lo_PacificE$Females, na.rm=TRUE)+sum(Lo_PacificE$Intersexes, na.rm=TRUE)` with `r sum(Lo_PacificE$Males, na.rm=TRUE)` males, `r sum(Lo_PacificE$Females, na.rm=TRUE)` females and `r sum(Lo_PacificE$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_PacificE)`. From India, the total number of sexed embryos is `r sum(Lo_IndianNE$Males, na.rm=TRUE)+sum(Lo_IndianNE$Females, na.rm=TRUE)+sum(Lo_IndianNE$Intersexes, na.rm=TRUE)` with `r sum(Lo_IndianNE$Males, na.rm=TRUE)` males, `r sum(Lo_IndianNE$Females, na.rm=TRUE)` females and `r sum(Lo_IndianNE$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_IndianNE)`. From Brazil, the total number of sexed embryos is `r sum(Lo_AtlanticWest$Males, na.rm=TRUE)+sum(Lo_AtlanticWest$Females, na.rm=TRUE)+sum(Lo_AtlanticWest$Intersexes, na.rm=TRUE)` with `r sum(Lo_AtlanticWest$Males, na.rm=TRUE)` males, `r sum(Lo_AtlanticWest$Females, na.rm=TRUE)` females and `r sum(Lo_AtlanticWest$Intersexes, na.rm=TRUE)` intersexes. The total number of incubation temperatures is `r nrow(Lo_AtlanticWest)`. # Map with locations ```{r} # pdf(file = "Figure 1.pdf", width = 7, height = 4, pointsize = 12) par(mar=c(5, 5, 1, 1)) map('world', fill=TRUE, col="lightgrey", border=NA) points(unique(Lo_Global$Longitude), unique(Lo_Global$Latitude), col="white", pch=16, cex=3) points(unique(Lo_Global$Longitude), unique(Lo_Global$Latitude), col="black", pch=1, cex=3) # text(unique(Lo_Global$Longitude), unique(Lo_Global$Latitude), unique(Lo_Global$Area)) tl <- c(1, 2, 3, 4, "", "") tp <- c(1, 3, 2, 4, 5, 6) text(unique(Lo_Global$Longitude)[tp], unique(Lo_Global$Latitude)[tp], labels = tl) segments(x0=unique(Lo_Global$Longitude)[5], y0=unique(Lo_Global$Latitude)[5], x1=unique(Lo_Global$Longitude)[5]-20, y1=unique(Lo_Global$Latitude)[5]-3) segments(x0=unique(Lo_Global$Longitude)[6], y0=unique(Lo_Global$Latitude)[6], x1=unique(Lo_Global$Longitude)[5]-20, y1=unique(Lo_Global$Latitude)[6]+3) text(unique(Lo_Global$Longitude)[5]-25, unique(Lo_Global$Latitude)[5]-4, labels = "5") text(unique(Lo_Global$Longitude)[5]-25, unique(Lo_Global$Latitude)[5]+5, labels = "6") degAxis(side=1, cex.axis=1) # in package sp par(xpd=TRUE) degAxis(side= 2, las=1, cex.axis=1, pos=-170, at=seq(from=-90, to=90, length.out = 5)) # dev.off() ``` # Fit TSD model using maximum likelihood with all data The flexit model is selected but standard error of some parameters is very high. We will continue with logistic model. Hill and A-logistic model are constrained asymetrical model. If asymetry is tested, it is better to use flexit model. The choice of the model can have profound implication for the estimate of PT and TRT. It will be used as a prior for Bayesian analysis ```{r} tsdL_Lo_Global_logistic <- with (Lo_Global, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic")) tsdL_Lo_Global_logistic$deviance tsdL_Lo_Global_logistic$df tsdL_Lo_Global_logistic$pvalue par(mar=c(4, 4, 1, 1)+0.4) par(xpd=FALSE) hist(tsdL_Lo_Global_logistic$NullDeviance, main="", las=1, xlab="Deviance for null distribution", ylim=c(0, 250)) segments(x0=tsdL_Lo_Global_logistic$deviance, x1=tsdL_Lo_Global_logistic$deviance, y0=0, y1=ScalePreviousPlot()$ylim["end"], lty=2, lwd=2) text(x=tsdL_Lo_Global_logistic$deviance-1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels=round(1-tsdL_Lo_Global_logistic$NullDeviancePvalue, digits = 2), pos=2) text(x=tsdL_Lo_Global_logistic$deviance+1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels=round(tsdL_Lo_Global_logistic$NullDeviancePvalue, digits = 2), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=210, labels = expression(bold("Logistic")), pos=4, cex=1.1) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=180, labels = paste("Deviance=", round(x = tsdL_Lo_Global_logistic$deviance, digits = 2)), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=160, labels = paste("df=", round(x = tsdL_Lo_Global_logistic$df, digits = 0)), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=140, labels = paste("p-value=", round(x = tsdL_Lo_Global_logistic$pvalue, digits = 2)), pos=4) par(xpd=TRUE) text(x=tsdL_Lo_Global_logistic$deviance, y=ScalePreviousPlot()$ylim["end"]*1.1, labels ="Observed deviance") tsdL_Lo_Global_Hill <- with (Lo_Global, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="Hill", replicate.NullDeviance = 0)) tsdL_Lo_Global_Alogistic <- with (Lo_Global, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="A-logistic", parameters.initial = c(tsdL_Lo_Global_logistic$par, K = -3), replicate.NullDeviance = 0)) tsdL_Lo_Global_flexit <- with (Lo_Global, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="flexit", parameters.initial = c(tsdL_Lo_Global_logistic$par, K1 = -3, K2 = 200))) tsdL_Lo_Global_flexit$deviance tsdL_Lo_Global_flexit$df tsdL_Lo_Global_flexit$pvalue par(mar=c(4, 4, 1, 1)+0.4) par(xpd=FALSE) hist(tsdL_Lo_Global_flexit$NullDeviance, main="", las=1, xlab="Deviance for null distribution", ylim=c(0, 250), xlim=c(20, 45)) segments(x0=tsdL_Lo_Global_flexit$deviance, x1=tsdL_Lo_Global_flexit$deviance, y0=0, y1=ScalePreviousPlot()$ylim["end"], lty=2, lwd=2) text(x=tsdL_Lo_Global_flexit$deviance-1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels=round(1-tsdL_Lo_Global_flexit$NullDeviancePvalue, digits = 2), pos=2) text(x=tsdL_Lo_Global_flexit$deviance+5, y=ScalePreviousPlot()$ylim["end"]*0.95, labels=round(tsdL_Lo_Global_flexit$NullDeviancePvalue, digits = 2), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=210, labels = expression(bold("Flexit")), pos=4, cex=1.1) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=180, labels = paste("Deviance=", round(x = tsdL_Lo_Global_flexit$deviance, digits = 2)), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=160, labels = paste("df=", round(x = tsdL_Lo_Global_flexit$df, digits = 0)), pos=4) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.70, y=140, labels = paste("p-value=", round(x = tsdL_Lo_Global_flexit$pvalue, digits = 2)), pos=4) par(xpd=TRUE) text(x=tsdL_Lo_Global_flexit$deviance, y=ScalePreviousPlot()$ylim["end"]*1.1, labels ="Observed deviance") par(xpd=FALSE) (comp0 <- compare_AICc(logistic=tsdL_Lo_Global_logistic, # Hill=tsdL_Lo_Global_Hill, # Alogistic=tsdL_Lo_Global_Alogistic, flexit=tsdL_Lo_Global_flexit)) result <- tsdL_Lo_Global_logistic out1 <- data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), AICc=as.character(round(result$AICc, digits = 2)), "deltaAICc"=as.character(round(comp0[1, "DeltaAICc"], digits = 2)), Akaike.Weight=as.character(round(comp0[1, "Akaike_weight"], digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) result <- tsdL_Lo_Global_flexit out1 <- rbind(out1, data.frame(model="Flexit", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1=paste(round(result$par["K1"], digits = 2), "SE", round(result$SE["K1"], digits = 2)), K2=paste(round(result$par["K2"], digits = 2), "SE", round(result$SE["K2"], digits = 2)), LnL=as.character(round(result$value, digits = 2)), AICc=as.character(round(result$AICc, digits = 2)), "deltaAICc"=as.character(round(comp0[2, "DeltaAICc"], digits = 2)), Akaike.Weight=as.character(round(comp0[2, "Akaike_weight"], digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) ) # write.xlsx(x=out1, file="Table2.xlsx", asTable = FALSE) plot(tsdL_Lo_Global_logistic) text(x=25, y=0.1, labels = "Logistic model", pos=4) P_TRT(x=tsdL_Lo_Global_logistic, replicate.CI = 10000)$P_TRT_quantiles plot(tsdL_Lo_Global_Hill) text(x=25, y=0.1, labels = "Hill model", pos=4) P_TRT(x=tsdL_Lo_Global_Hill, replicate.CI = 10000)$P_TRT_quantiles plot(tsdL_Lo_Global_Alogistic) text(x=25, y=0.1, labels = "A-logistic model", pos=4) P_TRT(x=tsdL_Lo_Global_Alogistic, replicate.CI = 10000)$P_TRT_quantiles plot(tsdL_Lo_Global_flexit) text(x=25, y=0.1, labels = "Flexit model", pos=4) P_TRT(x=tsdL_Lo_Global_flexit, replicate.CI = 10000)$P_TRT_quantiles SEfromHessian(tsdL_Lo_Global_flexit$hessian) abs(SEfromHessian(tsdL_Lo_Global_flexit$hessian)/tsdL_Lo_Global_flexit$par) ``` ```{r eval=FALSE} cairo_pdf(filename = "Figure 2.pdf", width=7, height = 9, pointsize = 14) layout(1:2) plot(tsdL_Lo_Global_logistic, xlim=c(24, 34)) text(x=24, y=0.4, labels = expression(bold("Logistic")), pos=4, cex=1.1) text(x=24, y=0.2, labels = paste("AICc=", round(x = comp0["logistic", "AICc"], digits = 2)), pos=4) text(x=24, y=0.1, labels = paste("Akaike weight=", round(x = comp0["logistic", "Akaike_weight"], digits = 2)), pos=4) par(xpd=TRUE) text(x=23, y=1.2, labels = "A", cex=2) plot(tsdL_Lo_Global_flexit, xlim=c(24, 34)) text(x=24, y=0.4, labels = expression(bold("Flexit")), pos=4, cex=1.1) text(x=24, y=0.2, labels = paste("AICc=", round(x = comp0["flexit", "AICc"], digits = 2)), pos=4) text(x=24, y=0.1, labels = paste("Akaike weight=", round(x = comp0["flexit", "Akaike_weight"], digits = 2)), pos=4) par(xpd=TRUE) text(x=23, y=1.2, labels = "B", cex=2) dev.off() ``` # Fit TSD model for each region when it is possible with logistic model ```{r} tsdL_Lo_PacificE_CostaRica_logistic <- with (Lo_PacificE_CostaRica, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic", parameters.initial = tsdL_Lo_Global_logistic$par)) tsdL_Lo_PacificE_CostaRica_logistic$pvalue tsdL_Lo_PacificE_CostaRica_logistic$NullDeviancePvalue tsdL_Lo_PacificE_CostaRica_logistic$deviance tsdL_Lo_PacificE_Mexico_logistic <- with (Lo_PacificE_Mexico, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic", parameters.initial = tsdL_Lo_Global_logistic$par)) tsdL_Lo_PacificE_Mexico_logistic$pvalue tsdL_Lo_PacificE_Mexico_logistic$NullDeviancePvalue tsdL_Lo_PacificE_Mexico_logistic$deviance tsdL_Lo_PacificE_logistic <- with (Lo_PacificE, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic", parameters.initial = tsdL_Lo_Global_logistic$par)) tsdL_Lo_PacificE_logistic$pvalue tsdL_Lo_PacificE_logistic$NullDeviancePvalue tsdL_Lo_PacificE_logistic$deviance tsdL_Lo_AtlanticWest_logistic <- with (Lo_AtlanticWest, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic", parameters.initial = tsdL_Lo_Global_logistic$par)) tsdL_Lo_AtlanticWest_logistic$pvalue tsdL_Lo_AtlanticWest_logistic$NullDeviancePvalue tsdL_Lo_AtlanticWest_logistic$deviance tsdL_Lo_IndianNE_logistic <- with (Lo_IndianNE, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="logistic", replicate.NullDeviance = 10000, parameters.initial = tsdL_Lo_Global_logistic$par)) tsdL_Lo_IndianNE_logistic$pvalue tsdL_Lo_IndianNE_logistic$NullDeviancePvalue tsdL_Lo_IndianNE_logistic$deviance ``` # Flexit model ```{r} tsdL_Lo_PacificE_flexit <- with (Lo_PacificE, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="flexit", replicate.NullDeviance = 10000, parameters.initial = tsdL_Lo_Global_flexit$par)) P_TRT(tsdL_Lo_PacificE_flexit, replicate.CI = 10000)$P_TRT_quantiles tsdL_Lo_IndianNE_flexit <- with (Lo_IndianNE, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="flexit", replicate.NullDeviance = 10000, parameters.initial = tsdL_Lo_Global_flexit$par)) P_TRT(tsdL_Lo_IndianNE_flexit, replicate.CI = 10000)$P_TRT_quantiles tsdL_Lo_AtlanticWest_flexit <- with (Lo_AtlanticWest, tsd(males=Males, females=Females, temperatures=Incubation.temperature, equation="flexit", replicate.NullDeviance = 10000, parameters.initial = tsdL_Lo_Global_flexit$par)) P_TRT(tsdL_Lo_AtlanticWest_flexit, replicate.CI = 10000)$P_TRT_quantiles ``` # Plot the results ## Pacific East ```{r} plot(tsdL_Lo_PacificE_logistic) P_TRT(tsdL_Lo_PacificE_logistic, replicate.CI = 10000)$P_TRT_quantiles ``` ## Pacific East - Mexico ```{r} plot(tsdL_Lo_PacificE_Mexico_logistic) P_TRT(tsdL_Lo_PacificE_Mexico_logistic)$P_TRT_quantiles ``` ## Pacific East - Costa Rica ```{r} plot(tsdL_Lo_PacificE_CostaRica_logistic) P_TRT(tsdL_Lo_PacificE_CostaRica_logistic)$P_TRT_quantiles ``` ## Atlantic West ```{r} plot(tsdL_Lo_AtlanticWest_logistic) P_TRT(tsdL_Lo_AtlanticWest_logistic, replicate.CI = 10000)$P_TRT_quantiles ``` ## Indian Ocean ```{r} plot(x=tsdL_Lo_IndianNE_logistic, xlim=c(24, 34)) P_TRT(tsdL_Lo_IndianNE_logistic, replicate.CI = 10000)$P_TRT_quantiles ``` ## Plot of all series ```{r eval=FALSE} pdf(file = "Figure 4.pdf", width=10, height = 14, pointsize = 14) layout(mat = matrix(1:6, ncol=2, byrow = FALSE)) plot(tsdL_Lo_PacificE_logistic, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="A: East Pacific", cex=2, pos=4) plot(tsdL_Lo_PacificE_Mexico_logistic, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="B: Mexico", cex=2, pos=4) plot(tsdL_Lo_PacificE_CostaRica_logistic, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="C: Costa Rica", cex=2, pos=4) plot(tsdL_Lo_AtlanticWest_logistic, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="D: West Atlantic", cex=2, pos=4) plot(x=tsdL_Lo_IndianNE_logistic, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="E: Northeast Indian", cex=2, pos=4) dev.off() ``` # Compare a global model to a RMU specific model A single model for all locations is not selected. It proved that there is differences among locations. ```{r} compare_BIC(global=tsdL_Lo_Global_logistic, separate=list(tsdL_Lo_AtlanticWest_logistic, tsdL_Lo_PacificE_logistic, tsdL_Lo_IndianNE_logistic)) ``` A single model for Mexican and Costa Rican locations is selected. It proved that there is no difference among Pacific East locations. ```{r} compare_BIC(PacificE=tsdL_Lo_PacificE_logistic, separate=list(tsdL_Lo_PacificE_Mexico_logistic, tsdL_Lo_PacificE_CostaRica_logistic)) ``` # Bayesian model to compare Costa Rica and Mexico and other locations with logistic model For all locations, priors are obtained from the global model. ## Costa Rica, East Pacific ```{r} pMCMC_LoPacificE_CostaRica <- tsd_MHmcmc_p(tsdL_Lo_PacificE_CostaRica_logistic, accept=TRUE) pMCMC_LoPacificE_CostaRica[, "Prior1"] <- tsdL_Lo_Global_logistic$par result_mcmc_tsd_LoPacificE_CostaRica <- tsd_MHmcmc(result=tsdL_Lo_PacificE_CostaRica_logistic, parametersMCMC=pMCMC_LoPacificE_CostaRica, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_PacificE_CostaRica_logistic, resultmcmc = result_mcmc_tsd_LoPacificE_CostaRica) plot(result_mcmc_tsd_LoPacificE_CostaRica, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoPacificE_CostaRica, parameters = "P", xlim=c(25, 35)) P_TRT_result_mcmc_tsd_LoPacificE_CostaRica <- P_TRT(resultmcmc = result_mcmc_tsd_LoPacificE_CostaRica, equation = "logistic", replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE_CostaRica$P_TRT_quantiles ``` ## Mexico, East Pacific ```{r} pMCMC_LoPacificE_Mexico <- tsd_MHmcmc_p(tsdL_Lo_PacificE_Mexico_logistic, accept=TRUE) pMCMC_LoPacificE_Mexico[, "Prior1"] <- tsdL_Lo_Global_logistic$par result_mcmc_tsd_LoPacificE_Mexico <- tsd_MHmcmc(result=tsdL_Lo_PacificE_Mexico_logistic, parametersMCMC=pMCMC_LoPacificE_Mexico, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_PacificE_Mexico_logistic, resultmcmc = result_mcmc_tsd_LoPacificE_Mexico) plot(result_mcmc_tsd_LoPacificE_Mexico, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoPacificE_Mexico, parameters = "P", xlim=c(25, 35)) P_TRT_result_mcmc_tsd_LoPacificE_Mexico <- P_TRT(resultmcmc = result_mcmc_tsd_LoPacificE_Mexico, equation = "logistic", replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE_Mexico$P_TRT_quantiles ``` ## Northeast India ```{r} pMCMC_LoIndianNE <- tsd_MHmcmc_p(tsdL_Lo_IndianNE_logistic, accept=TRUE) pMCMC_LoIndianNE[, "Prior1"] <- tsdL_Lo_Global_logistic$par result_mcmc_tsd_LoIndianNE <- tsd_MHmcmc(result=tsdL_Lo_IndianNE_logistic, parametersMCMC=pMCMC_LoIndianNE, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_IndianNE_logistic, resultmcmc = result_mcmc_tsd_LoIndianNE) plot(result_mcmc_tsd_LoIndianNE, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoIndianNE, parameters = "P", xlim=c(25, 35)) P_TRT_result_mcmc_tsd_LoIndianNE <- P_TRT(resultmcmc = result_mcmc_tsd_LoIndianNE, equation = "logistic", replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoIndianNE$P_TRT_quantiles ``` ```{r} layout(matrix(1:3, ncol=3)) plot(result_mcmc_tsd_LoIndianNE, parameters = "P", xlim=c(25, 35), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="A", cex=2) plot(result_mcmc_tsd_LoIndianNE, parameters = "S", xlim=c(-3, 2), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="B", cex=2) rd <- sample(1:100000, 10000) plot(as.numeric(result_mcmc_tsd_LoIndianNE$resultMCMC[[1]][rd, "P"]), as.numeric(result_mcmc_tsd_LoIndianNE$resultMCMC[[1]][rd, "S"]), pch=".", xlab="P", ylab="S", bty="n", las=1) par(xpd=TRUE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.95, labels="C", cex=2) ``` ## West Atlantic ```{r} pMCMC_LoAtlanticWest <- tsd_MHmcmc_p(tsdL_Lo_AtlanticWest_logistic, accept=TRUE) pMCMC_LoAtlanticWest[, "Prior1"] <- tsdL_Lo_Global_logistic$par result_mcmc_tsd_LoAtlanticWest <- tsd_MHmcmc(result=tsdL_Lo_AtlanticWest_logistic, parametersMCMC=pMCMC_LoAtlanticWest, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_AtlanticWest_logistic, resultmcmc = result_mcmc_tsd_LoAtlanticWest) plot(result_mcmc_tsd_LoAtlanticWest, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoAtlanticWest, parameters = "P", xlim=c(25, 35)) P_TRT_result_mcmc_tsd_LoAtlanticWest <- P_TRT(resultmcmc = result_mcmc_tsd_LoAtlanticWest, equation = "logistic", replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT_quantiles ``` ```{r} layout(matrix(1:3, ncol=3)) plot(result_mcmc_tsd_LoAtlanticWest, parameters = "P", xlim=c(25, 35), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="A", cex=2) plot(result_mcmc_tsd_LoAtlanticWest, parameters = "S", xlim=c(-3, 2), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="B", cex=2) rd <- sample(1:100000, 10000) plot(as.numeric(result_mcmc_tsd_LoAtlanticWest$resultMCMC[[1]][rd, "P"]), as.numeric(result_mcmc_tsd_LoAtlanticWest$resultMCMC[[1]][rd, "S"]), pch=".", xlab="P", ylab="S", bty="n", las=1) par(xpd=TRUE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.95, labels="C", cex=2) ``` ## East Pacific ```{r} pMCMC_LoPacificE <- tsd_MHmcmc_p(tsdL_Lo_PacificE_logistic, accept=TRUE) pMCMC_LoPacificE[, "Prior1"] <- tsdL_Lo_Global_logistic$par result_mcmc_tsd_PacificE <- tsd_MHmcmc(result=tsdL_Lo_PacificE_logistic, parametersMCMC=pMCMC_LoPacificE, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_PacificE_logistic, resultmcmc = result_mcmc_tsd_PacificE) P_TRT_result_mcmc_tsd_LoPacificE <- P_TRT(resultmcmc = result_mcmc_tsd_PacificE, equation = "logistic", replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE$P_TRT_quantiles ``` ```{r} # pdf(file = "Figure 3.pdf", width=10, height = 4, pointsize = 12) layout(matrix(1:3, ncol=3)) plot(result_mcmc_tsd_PacificE, parameters = "P", xlim=c(25, 35), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="A", cex=2) plot(result_mcmc_tsd_PacificE, parameters = "S", xlim=c(-3, 2), ylim=c(0, 8), legend = FALSE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["end"]*0.95, labels="B", cex=2) rd <- sample(1:100000, 10000) plot(as.numeric(result_mcmc_tsd_PacificE$resultMCMC[[1]][rd, "P"]), as.numeric(result_mcmc_tsd_PacificE$resultMCMC[[1]][rd, "S"]), pch=".", xlab="P", ylab="S", bty="n", las=1) par(xpd=TRUE) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.1, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.95, labels="C", cex=2) # dev.off() ``` ## Plot of all series ```{r eval=FALSE} pdf(file = "Figure 5.pdf", width=10, height = 14, pointsize = 14) layout(mat = matrix(1:6, ncol=2, byrow = FALSE)) plot(tsdL_Lo_PacificE_logistic, resultmcmc = result_mcmc_tsd_PacificE, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="A: East Pacific", cex=2, pos=4) plot(tsdL_Lo_PacificE_Mexico_logistic, resultmcmc = result_mcmc_tsd_LoPacificE_Mexico, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="B: Mexico", cex=2, pos=4) plot(tsdL_Lo_PacificE_CostaRica_logistic, resultmcmc = result_mcmc_tsd_LoPacificE_CostaRica, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="C: Costa Rica", cex=2, pos=4) plot(tsdL_Lo_AtlanticWest_logistic, resultmcmc = result_mcmc_tsd_LoAtlanticWest, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="D: West Atlantic", cex=2, pos=4) plot(x=tsdL_Lo_IndianNE_logistic, resultmcmc = result_mcmc_tsd_LoIndianNE, xlim=c(24, 34)) text(x=ScalePreviousPlot()$xlim["begin"]+ScalePreviousPlot()$xlim["range"]*0.02, y=ScalePreviousPlot()$ylim["begin"]+ScalePreviousPlot()$ylim["range"]*0.1, labels="E: Northeast Indian", cex=2, pos=4) dev.off() ``` # TRT and P relationship for 4 locations from Bayesian analysis with logistic model Costa Rica and Mexico (East Pacific RMU) and India (Northeast Indian RMU) cannot be differentiated. Brazil (West Atlantic RMU) is different from all others. ```{r} # pdf(file = "Figure 6.pdf", width = 7, height = 7, pointsize = 14) layout(1) par(mar=c(4, 4, 1, 1)) rd <- runif(20000, min=1, max = 100001) plot(x = P_TRT_result_mcmc_tsd_LoPacificE_Mexico$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_Mexico$P_TRT[rd, "TRT"], pch=".", xlim=c(27, 33), ylim=c(0, 8), xlab="Pivotal temperature", ylab="Transitional range of temperatures 5%", bty="n", las=1, col=gray.colors(4)[4]) rd <- runif(20000, min=1, max = 100001) points(x = P_TRT_result_mcmc_tsd_LoIndianNE$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoIndianNE$P_TRT[rd, "TRT"], pch=".", col=gray.colors(4)[2]) rd <- runif(5000, min=1, max = 100001) points(x = P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "TRT"], pch=".", col=gray.colors(4)[1]) rd <- runif(5000, min=1, max = 100001) points(x = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica$P_TRT[rd, "TRT"], pch=".", col=gray.colors(4)[3]) rd <- runif(1000, min=1, max = 100001) points(x = P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "TRT"], pch=".", col=gray.colors(4)[1]) l <- 0.80 ellipse <- dataEllipse(x=P_TRT_result_mcmc_tsd_LoPacificE_Mexico$P_TRT[rd, "PT"], y=P_TRT_result_mcmc_tsd_LoPacificE_Mexico$P_TRT[rd, "TRT"], levels=l, draw=FALSE) polygon(ellipse[, 1], ellipse[, 2], border = "black", col=NULL, lwd=1) segments(x0=30.2, x1=31, y0=0.4, y1=0.1) text(x=31.1, y=0.08, labels = "Mexico", pos=4) ellipse <- dataEllipse(x=P_TRT_result_mcmc_tsd_LoIndianNE$P_TRT[rd, "PT"], y=P_TRT_result_mcmc_tsd_LoIndianNE$P_TRT[rd, "TRT"], levels=l, draw=FALSE) polygon(ellipse[, 1], ellipse[, 2], border ="black", col=NULL, lwd=1) segments(x0=28.8, x1=28, y0=2, y1=2) text(x=28, y=1.95, labels = "India", pos=2) ellipse <- dataEllipse(x=P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "PT"], y=P_TRT_result_mcmc_tsd_LoAtlanticWest$P_TRT[rd, "TRT"], levels=l, draw=FALSE) polygon(ellipse[, 1], ellipse[, 2], border = "black", col=NULL, lwd=1) segments(x0=31, x1=32, y0=2, y1=2) text(x=32.05, y=1.95, labels = "Brazil", pos=4) ellipse <- dataEllipse(x=P_TRT_result_mcmc_tsd_LoPacificE_CostaRica$P_TRT[rd, "PT"], y=P_TRT_result_mcmc_tsd_LoPacificE_CostaRica$P_TRT[rd, "TRT"], levels=l, draw=FALSE) polygon(ellipse[, 1], ellipse[, 2], border = "black", col=NULL, lwd=1) segments(x0=30.5, x1=32, y0=3, y1=3) text(x=32.05, y=2.95, labels = "Costa Rica", pos=4) legend("topleft", legend=c("Brazil, West Atlantic", "India, Northeast Indian", "Costa Rica, East Pacific", "Mexico, East Pacific"), pch=rep(19, 4), col=gray.colors(4), cex=0.8) # dev.off() ``` # Tentative to use MCMC with flexit model ## East Pacific ```{r} tsdL_Lo_PacificE_flexit <- tsdL_Lo_Global_flexit tsdL_Lo_PacificE_flexit$males <- Lo_PacificE$Males tsdL_Lo_PacificE_flexit$females <- Lo_PacificE$Females tsdL_Lo_PacificE_flexit$N <- Lo_PacificE$Males + Lo_PacificE$Females tsdL_Lo_PacificE_flexit$temperatures <- Lo_PacificE$Incubation.temperature pMCMC_LoPacificE <- tsd_MHmcmc_p(tsdL_Lo_PacificE_flexit, accept=TRUE) result_mcmc_tsd_LoPacificE_flexit <- tsd_MHmcmc(result=tsdL_Lo_PacificE_flexit, parametersMCMC=pMCMC_LoPacificE, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) layout(1) plot(tsdL_Lo_PacificE_flexit, resultmcmc = result_mcmc_tsd_LoPacificE_flexit) plot(result_mcmc_tsd_LoPacificE_flexit, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoPacificE_flexit, parameters = "P", xlim=c(25, 35)) plot(result_mcmc_tsd_LoPacificE_flexit, parameters = "K1", xlim=c(-10, 10)) plot(result_mcmc_tsd_LoPacificE_flexit, parameters = "K2", xlim=c(150, 250)) P_TRT_result_mcmc_tsd_LoPacificE_flexit <- P_TRT(resultmcmc = result_mcmc_tsd_LoPacificE_flexit, replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE_flexit$P_TRT_quantiles ``` ## Costa Rica, East Pacific ```{r} tsdL_Lo_PacificE_CostaRica_flexit <- tsdL_Lo_Global_flexit tsdL_Lo_PacificE_CostaRica_flexit$males <- Lo_PacificE_CostaRica$Males tsdL_Lo_PacificE_CostaRica_flexit$females <- Lo_PacificE_CostaRica$Females tsdL_Lo_PacificE_CostaRica_flexit$N <- Lo_PacificE_CostaRica$Males + Lo_PacificE_CostaRica$Females tsdL_Lo_PacificE_CostaRica_flexit$temperatures <- Lo_PacificE_CostaRica$Incubation.temperature pMCMC_LoPacificE_CostaRica <- tsd_MHmcmc_p(tsdL_Lo_PacificE_CostaRica_flexit, accept=TRUE) result_mcmc_tsd_LoPacificE_CostaRica_flexit <- tsd_MHmcmc(result=tsdL_Lo_PacificE_CostaRica_flexit, parametersMCMC=pMCMC_LoPacificE_CostaRica, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_PacificE_CostaRica_flexit, resultmcmc = result_mcmc_tsd_LoPacificE_CostaRica_flexit) plot(result_mcmc_tsd_LoPacificE_CostaRica_flexit, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoPacificE_CostaRica_flexit, parameters = "P", xlim=c(25, 35)) plot(result_mcmc_tsd_LoPacificE_CostaRica_flexit, parameters = "K1", xlim=c(-10, 10)) plot(result_mcmc_tsd_LoPacificE_CostaRica_flexit, parameters = "K2", xlim=c(150, 250)) P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit <- P_TRT(resultmcmc = result_mcmc_tsd_LoPacificE_CostaRica_flexit, replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit$P_TRT_quantiles ``` ## Mexico, East Pacific ```{r} tsdL_Lo_PacificE_Mexico_flexit <- tsdL_Lo_Global_flexit tsdL_Lo_PacificE_Mexico_flexit$males <- Lo_PacificE_Mexico$Males tsdL_Lo_PacificE_Mexico_flexit$females <- Lo_PacificE_Mexico$Females tsdL_Lo_PacificE_Mexico_flexit$N <- Lo_PacificE_Mexico$Males + Lo_PacificE_Mexico$Females tsdL_Lo_PacificE_Mexico_flexit$temperatures <- Lo_PacificE_Mexico$Incubation.temperature pMCMC_LoPacificE_Mexico <- tsd_MHmcmc_p(tsdL_Lo_PacificE_Mexico_flexit, accept=TRUE) result_mcmc_tsd_LoPacificE_Mexico_flexit <- tsd_MHmcmc(result=tsdL_Lo_PacificE_Mexico_flexit, parametersMCMC=pMCMC_LoPacificE_Mexico, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_PacificE_Mexico_flexit, resultmcmc = result_mcmc_tsd_LoPacificE_Mexico_flexit) plot(result_mcmc_tsd_LoPacificE_Mexico_flexit, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoPacificE_Mexico_flexit, parameters = "P", xlim=c(25, 35)) plot(result_mcmc_tsd_LoPacificE_Mexico_flexit, parameters = "K1", xlim=c(-100, 100)) plot(result_mcmc_tsd_LoPacificE_Mexico_flexit, parameters = "K2", xlim=c(150, 250)) P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit <- P_TRT(resultmcmc = result_mcmc_tsd_LoPacificE_Mexico_flexit, replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit$P_TRT_quantiles ``` ## Brazil, West Atlantic ```{r} tsdL_Lo_AtlanticWest_Brazil_flexit <- tsdL_Lo_Global_flexit tsdL_Lo_AtlanticWest_Brazil_flexit$males <- Lo_AtlanticWest$Males tsdL_Lo_AtlanticWest_Brazil_flexit$females <- Lo_AtlanticWest$Females tsdL_Lo_AtlanticWest_Brazil_flexit$N <- Lo_AtlanticWest$Males + Lo_AtlanticWest$Females tsdL_Lo_AtlanticWest_Brazil_flexit$temperatures <- Lo_AtlanticWest$Incubation.temperature pMCMC_LoAtlanticWest_Brazil <- tsd_MHmcmc_p(tsdL_Lo_AtlanticWest_Brazil_flexit, accept=TRUE) result_mcmc_tsd_LoAtlanticWest_flexit <- tsd_MHmcmc(result=tsdL_Lo_AtlanticWest_Brazil_flexit, parametersMCMC=pMCMC_LoAtlanticWest_Brazil, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_AtlanticWest_Brazil_flexit, resultmcmc = result_mcmc_tsd_LoAtlanticWest_flexit) plot(result_mcmc_tsd_LoAtlanticWest_flexit, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoAtlanticWest_flexit, parameters = "P", xlim=c(25, 35)) plot(result_mcmc_tsd_LoAtlanticWest_flexit, parameters = "K1", xlim=c(-100, 100)) plot(result_mcmc_tsd_LoAtlanticWest_flexit, parameters = "K2", xlim=c(150, 250)) P_TRT_result_mcmc_tsd_LoAtlanticWest_flexit <- P_TRT(resultmcmc = result_mcmc_tsd_LoAtlanticWest_flexit, replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoAtlanticWest_flexit$P_TRT_quantiles ``` ## India, Northeast India ```{r} tsdL_Lo_IndianNE_India_flexit <- tsdL_Lo_Global_flexit tsdL_Lo_IndianNE_India_flexit$males <- Lo_IndianNE$Males tsdL_Lo_IndianNE_India_flexit$females <- Lo_IndianNE$Females tsdL_Lo_IndianNE_India_flexit$N <- Lo_IndianNE$Males + Lo_IndianNE$Females tsdL_Lo_IndianNE_India_flexit$temperatures <- Lo_IndianNE$Incubation.temperature pMCMC_LoIndianNE_India <- tsd_MHmcmc_p(tsdL_Lo_AtlanticWest_Brazil_flexit, accept=TRUE) result_mcmc_tsd_LoIndianNE_India_flexit <- tsd_MHmcmc(result=tsdL_Lo_IndianNE_India_flexit, parametersMCMC=pMCMC_LoIndianNE_India, n.iter=100000, n.chains = 1, n.adapt = 0, thin=1, trace=FALSE, adaptive = TRUE) plot(tsdL_Lo_IndianNE_India_flexit, resultmcmc = result_mcmc_tsd_LoIndianNE_India_flexit) plot(result_mcmc_tsd_LoIndianNE_India_flexit, parameters = "S", xlim=c(-2, 2)) plot(result_mcmc_tsd_LoIndianNE_India_flexit, parameters = "P", xlim=c(25, 35)) plot(result_mcmc_tsd_LoIndianNE_India_flexit, parameters = "K1", xlim=c(-100, 100)) plot(result_mcmc_tsd_LoIndianNE_India_flexit, parameters = "K2", xlim=c(150, 250)) P_TRT_result_mcmc_tsd_LoIndianNE_India_flexit <- P_TRT(resultmcmc = result_mcmc_tsd_LoIndianNE_India_flexit, replicate.CI = 100000) P_TRT_result_mcmc_tsd_LoIndianNE_India_flexit$P_TRT_quantiles ``` # TRT and P relationship for 4 locations from Bayesian analysis with flexit model Costa Rica (RMU Pacific East), Mexico (RMU Pacific East), and Brazil (Atlantic West) cannot be differentiated. On the other hand, India (RMU Atlantic West) is different from all others. ```{r} par(mar=c(4, 4, 1, 1)) rd <- runif(20000, min=1, max = 100001) plot(x = P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit$P_TRT[rd, "TRT"], pch=".", xlim=c(27, 33), ylim=c(0, 8), xlab="Pivotal temperature", ylab="Transitional range of temperatures 5%", bty="n", las=1, col="black") points(x = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit$P_TRT[rd, "TRT"], pch=".", col="red") points(x = P_TRT_result_mcmc_tsd_LoIndianNE_India_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoIndianNE_India_flexit$P_TRT[rd, "TRT"], pch=".", col="green") points(x = P_TRT_result_mcmc_tsd_LoAtlanticWest_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoAtlanticWest_flexit$P_TRT[rd, "TRT"], pch=".", col="blue") points(x = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_CostaRica_flexit$P_TRT[rd, "TRT"], pch=".", col="red") points(x = P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit$P_TRT[rd, "PT"], y = P_TRT_result_mcmc_tsd_LoPacificE_Mexico_flexit$P_TRT[rd, "TRT"], pch=".", col="black") legend("topleft", legend=c("Mexico, East Pacific", "Costa Rica, East Pacific", "Brazil, West Atlantic", "India, Northeast Indian"), pch=rep(19, 4), col=c("black", "red", "blue", "green"), cex=0.8) ``` # Prepare the Table 3 ```{r eval=FALSE} result <- tsdL_Lo_IndianNE_logistic out1 <- data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) result <- tsdL_Lo_AtlanticWest_logistic out1 <- rbind(out1, data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) ) result <- tsdL_Lo_PacificE_logistic out1 <- rbind(out1, data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) ) result <- tsdL_Lo_PacificE_CostaRica_logistic out1 <- rbind(out1, data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) ) result <- tsdL_Lo_PacificE_Mexico_logistic out1 <- rbind(out1, data.frame(model="Logistic", P=paste(round(result$par["P"], digits = 2), "SE", round(result$SE["P"], digits = 2)), S=paste(round(result$par["S"], digits = 2), "SE", round(result$SE["S"], digits = 2)), K1="", K2="", LnL=as.character(round(result$value, digits = 2)), Deviance=as.character(round(result$deviance, digits = 2)), df=as.character(round(result$df, digits = 0)), p.value=as.character(round(result$pvalue, digits = 6)), p.Deviance.Null.model=as.character(round(result$NullDeviancePvalue, digits = 2)) ) ) # write.xlsx(x=out1, file="Table3.xlsx", asTable = FALSE) ``` ```{r eval=FALSE} save.image("ALLResults.Rdata") ```