R/plot.tsd.R

In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth

#' plot.tsd plot result of tsd() that best describe temperature-dependent sex determination
#' @title Plot results of tsd() that best describe temperature-dependent sex determination
#' @author Marc Girondot
#' @return Nothing
#' @param x A result file generated by tsd()
#' @param ... Parameters for plot()
#' @param show.observations Should the observations be shown
#' @param show.model Should the model be shown
#' @param show.PTRT Should the P and TRT information be shown
#' @param resultmcmc A result of tsd_MHmcmc()
#' @param chain What chain to be used is resultmcmc is provided
#' @param temperatures.plot Temperatures used for showing curves of sex ratio
#' @param durations.plot Durations used for showing curves of sex ratio
#' @param replicate.CI replicate.CI replicates from the hessian matrix to estimate CI 
#' @param range.CI The range of confidence interval for estimation, default=0.95
#' @param l Sex ratio limits to define TRT are l and 1-l
#' @param las.x las parameter for x axis
#' @param las.y las parameter for y axis
#' @param lab.PT Label to describe pivotal temperature
#' @param lab.TRT Label to describe transitional range of temperature
#' @param males.freq Should the graph uses males relative frequency [TRUE] or females [FALSE]
#' @param mar The par("mar") parameter
#' @param col.TRT The color of TRT
#' @param col.TRT.CI The color of CI of TRT based on range.CI
#' @param col.PT.CI The color of CI of PT based on range.CI
#' @param show.CI Do the CI for the curve should be shown
#' @param warn Do the warnings must be shown ? TRUE or FALSE
#' @description Plot the estimates that best describe temperature-dependent sex determination.
#' @references Girondot, M. 1999. Statistical description of temperature-dependent sex determination using maximum likelihood. Evolutionary Ecology Research, 1, 479-486.
#' @references Godfrey, M.H., Delmas, V., Girondot, M., 2003. Assessment of patterns of temperature-dependent sex determination using maximum likelihood model selection. Ecoscience 10, 265-272.
#' @references 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.
#' @references Girondot M., Submited. On the concept of embryological thermosensitive period for sex determination in reptiles.
#' @examples
#' \dontrun{
#' CC_AtlanticSW <- subset(DatabaseTSD, RMU=="Atlantic, SW" & 
#'                           Species=="Caretta caretta" & (!is.na(Sexed) & Sexed!=0))
#' tsdL <- with (CC_AtlanticSW, tsd(males=Males, females=Females, 
#'                                  temperatures=Incubation.temperature-Correction.factor, 
#'                                  equation="logistic"))
#' plot(tsdL)
#' }
#' @family Functions for temperature-dependent sex determination
#' @method plot tsd
#' @export



plot.tsd <- function(x, ..., 
                     show.observations=TRUE, 
                     show.model=TRUE, 
                     males.freq=TRUE, 
                     show.PTRT=TRUE, 
                     las.x=1, 
                     las.y=1, 
                     lab.PT=paste("Pivotal ", x$type), 
                     resultmcmc = NULL,
                     chain=1, 
                     l=0.05, 
                     replicate.CI=10000, 
                     range.CI=0.95, 
                     mar=c(4, 4, 4, 1)+0.4,
                     temperatures.plot=seq(from=25, to=35, by=0.1), 
                     durations.plot=seq(from=40, to=70, by=0.1), 
                     lab.TRT=paste0("Transitional range of ",  x$type, "s l=",x$l*100,"%"), 
                     col.TRT="gray", 
                     col.TRT.CI=rgb(0.8, 0.8, 0.8, 0.5), 
                     col.PT.CI=rgb(0.8, 0.8, 0.8, 0.5), 
                     show.CI=TRUE, 
                     warn = TRUE) {
  
  # show.observations=TRUE; males.freq=TRUE; show.PTRT = TRUE; las.x=1; las.y=1; lab.PT=paste("Pivotal ", x$type); resultmcmc = NULL; chain=1; l=0.05; replicate.CI=10000; temperatures.plot=seq(from=20, to=35, by=0.1); durations.plot=seq(from=40, to=70, by=0.1);TRT.limits=c(9, 90); precision=15; range.CI=0.95; mar=c(4, 4, 6, 1)+0.4; lab.TRT=paste0("Transitional range of ", x$type, "s l=",x$l*100,"%"); col.TRT="gray"; col.TRT.CI=rgb(0.8, 0.8, 0.8, 0.5); col.PT.CI=rgb(0.8, 0.8, 0.8, 0.5); show.CI=TRUE; warn=TRUE
  
  males <- x$males
  females <- x$females
  N <- x$N
  temperatures <- x$temperatures
  equation <- x$equation  
  
  if (x$type != "temperature") temperatures.plot <- durations.plot
  
  L <- list(...) # L <- list()
  par(mar=mar)
  
  xlll <- ifelse(x$type=="temperature", expression("Temperature in " * degree * "C"), 
                 "Duration in days")
  
  if (males.freq) {
    L1 <- modifyList(list(x=temperatures, y=males/N, bty="n", type="n", xlab=xlll, 
                          ylab="Male relative frequency"), L)
  } else {
    L1 <- modifyList(list(x=temperatures, y=females/N, bty="n", type="n", xlab=xlll, 
                          ylab="Female relative frequency"), L)	
  }
  L1 <- modifyList(L1, list(ylim=c(0,1), xaxt="n", las=las.y))
  
  if (is.null(L$xlim)) {
    L1 <- modifyList(L1, list(xlim=c(floor(min(temperatures)), floor(1+max(temperatures)))))
  }
  
  temperatures.plot <- seq(from=L1$xlim[1], to=L1$xlim[2], by=0.1) 
  
  o <- P_TRT(x=x, resultmcmc=resultmcmc, chain=chain, l=l, 
             replicate.CI=replicate.CI, temperatures=temperatures.plot, 
             probs = c((1-range.CI)/2, 0.5, 1-(1-range.CI)/2), warn=warn)
  
  
  L2 <- L1[!(names(L1) %in% c("errbar.tick", "errbar.lwd"
                              , "errbar.lty", "errbar.col"
                              , "errbar.y.polygon"
                              , "errbar.y.polygon.list"))]
  
  a <- do.call(plot, L2) 
  
  x2 <- (par("usr")[1]+par("usr")[2]*26)/27
  x1 <- x2*26-par("usr")[2]/0.04
  
  cex.x <- par("cex.axis")
  if (!is.null(L$cex.axis)) cex.x <- L$cex.axis
  
  ptax <- TRUE
  if (!is.null(L[["axes"]])) if (isFALSE(L[["axes"]])) ptax <- FALSE
  if (!is.null(L[["xaxt"]])) if (L[["xaxt"]] == "n") ptax <- FALSE
  
  
  if (ptax) axis(1, at=x1:x2, las=las.x, cex.axis=cex.x)
  
  # je trace la TRT centree sur P
  
  par(xpd=FALSE)
  
  if ((equation!="GSD") & (show.PTRT)) {
    if (any(colnames(o$P_TRT_quantiles)=="PT")) {
      par(xpd=FALSE)
      polygon(c(o$P_TRT_quantiles[2, "lower.limit.TRT"], o$P_TRT_quantiles[2, "lower.limit.TRT"], o$P_TRT_quantiles[2, "higher.limit.TRT"], o$P_TRT_quantiles[2, "higher.limit.TRT"]), c(0,1,1,0), border=NA, col=col.TRT)  
      # limites de la limite basse de la TRT
      polygon(c(o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la limite haute de la TRT
      polygon(c(o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la PT
      polygon(c(o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[3, "PT"], o$P_TRT_quantiles[3, "PT"]), c(0,1,1,0), border=NA, col=col.PT.CI)  
      par(xpd=TRUE)
      segments(o$P_TRT_quantiles[2, "PT"], 0, o$P_TRT_quantiles[2, "PT"], 1.05, lty=4)
      segments(o$P_TRT_quantiles[2, "lower.limit.TRT"], 0, o$P_TRT_quantiles[2, "lower.limit.TRT"], 1.15, lty=3)
      segments(o$P_TRT_quantiles[2, "higher.limit.TRT"], 0, o$P_TRT_quantiles[2, "higher.limit.TRT"], 1.15, lty=3)
      text(x=o$P_TRT_quantiles[2, "PT"], y=1.1, lab.PT)
      text(x=o$P_TRT_quantiles[2, "PT"], y=1.2, lab.TRT)
      par(xpd=FALSE)
    } else {
      par(xpd=FALSE)
      polygon(c(o$P_TRT_quantiles[2, "lower.limit.TRT_low"], o$P_TRT_quantiles[2, "lower.limit.TRT_low"], o$P_TRT_quantiles[2, "higher.limit.TRT_low"], o$P_TRT_quantiles[2, "higher.limit.TRT_low"]), c(0,1,1,0), border=NA, col=col.TRT)  
      # limites de la limite basse de la TRT
      polygon(c(o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la limite haute de la TRT
      polygon(c(o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la PT
      polygon(c(o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[3, "PT_low"], o$P_TRT_quantiles[3, "PT_low"]), c(0,1,1,0), border=NA, col=col.PT.CI)  
      par(xpd=TRUE)
      segments(o$P_TRT_quantiles[2, "PT_low"], 0, o$P_TRT_quantiles[2, "PT_low"], 1.05, lty=4)
      segments(o$P_TRT_quantiles[2, "lower.limit.TRT_low"], 0, o$P_TRT_quantiles[2, "lower.limit.TRT_low"], 1.15, lty=3)
      segments(o$P_TRT_quantiles[2, "higher.limit.TRT_low"], 0, o$P_TRT_quantiles[2, "higher.limit.TRT_low"], 1.15, lty=3)
      text(x=o$P_TRT_quantiles[2, "PT_low"], y=1.1, lab.PT)
      text(x=o$P_TRT_quantiles[2, "PT_low"], y=1.2, lab.TRT)
      par(xpd=FALSE)
      polygon(c(o$P_TRT_quantiles[2, "lower.limit.TRT_high"], o$P_TRT_quantiles[2, "lower.limit.TRT_high"], o$P_TRT_quantiles[2, "higher.limit.TRT_high"], o$P_TRT_quantiles[2, "higher.limit.TRT_high"]), c(0,1,1,0), border=NA, col=col.TRT)  
      # limites de la limite basse de la TRT
      polygon(c(o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la limite haute de la TRT
      polygon(c(o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"]), c(0,1,1,0), border=NA, col=col.TRT.CI)
      # limites de la PT
      polygon(c(o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[3, "PT_high"], o$P_TRT_quantiles[3, "PT_high"]), c(0,1,1,0), border=NA, col=col.PT.CI)  
      par(xpd=TRUE)
      segments(o$P_TRT_quantiles[2, "PT_high"], 0, o$P_TRT_quantiles[2, "PT_high"], 1.05, lty=4)
      segments(o$P_TRT_quantiles[2, "lower.limit.TRT_high"], 0, o$P_TRT_quantiles[2, "lower.limit.TRT_high"], 1.15, lty=3)
      segments(o$P_TRT_quantiles[2, "higher.limit.TRT_high"], 0, o$P_TRT_quantiles[2, "higher.limit.TRT_high"], 1.15, lty=3)
      text(x=o$P_TRT_quantiles[2, "PT_high"], y=1.1, lab.PT)
      text(x=o$P_TRT_quantiles[2, "PT_high"], y=1.2, lab.TRT)
      par(xpd=FALSE)
    }
  }
  
  if (show.observations) {
    if (males.freq) {  
      b <- getFromNamespace(".BinomialConfidence", ns="HelpersMG")(males,N)
      L1 <- modifyList(list(x=temperatures, y=males/N, bty="n", type="p", ylim=c(0,1), y.plus = b[,3], y.minus = b[,2]), L)
    } else {
      b <- getFromNamespace(".BinomialConfidence", ns="HelpersMG")(females,N)
      L1 <- modifyList(list(x=temperatures, y=females/N, bty="n", type="p", ylim=c(0,1), y.plus = b[,3], y.minus = b[,2]), L)
    }
    L1 <- modifyList(L1, list(ylim=c(0,1), xlab="", ylab="", 
                              main="", axes=FALSE, xlim=c(x1, x2)))
    
    par(xpd=FALSE)
    
    par(new=TRUE)
    
    a <- do.call(plot_errbar, L1) 
  }
  
  out_sr <- NULL
  
  if ((!is.null(o$sexratio_quantiles)) & (show.model)) {
    
    out_sr <- t(o$sexratio_quantiles)
    out_sr <- cbind(Temperatures=as.numeric(colnames(o$sexratio_quantiles)), 
                    out_sr)
    
    xi <- as.numeric(colnames(o$sexratio_quantiles))
    p <- o$sexratio_quantiles[2, ]
    
    par(new=TRUE)
    if (males.freq) {   
      L1 <- modifyList(list(x=xi, y=p, bty="n"), L)
    } else {
      L1 <- modifyList(list(x=xi, y=1-p, bty="n"), L)
    }
    L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", xlim=c(x1, x2)))
    
    L2 <- L1[!(names(L1) %in% c("errbar.tick", "errbar.lwd"
                                , "errbar.lty", "errbar.col"
                                , "errbar.y.polygon"
                                , "errbar.y.polygon.list"))]
    par(xpd=FALSE)
    a <- do.call(plot, L2)
    
    if (show.CI) {
      
      pm <- o$sexratio_quantiles[1, ]
      pp <- o$sexratio_quantiles[3, ]
      
      par(new=TRUE)
      if (males.freq) {
        L1 <- modifyList(list(x=xi, y=pm, bty="n"), L)
      } else {
        L1 <- modifyList(list(x=xi, y=1-pm, bty="n"), L)
      }
      L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", lty=2, xlim=c(x1, x2)))
      
      L2 <- L1[!(names(L1) %in% c("errbar.tick", "errbar.lwd"
                                  , "errbar.lty", "errbar.col"
                                  , "errbar.y.polygon"
                                  , "errbar.y.polygon.list"))]
      par(xpd=FALSE)
      a <- do.call(plot, L2) 
      
      par(new=TRUE)
      if (males.freq) {   
        L1 <- modifyList(list(x=xi, y=pp, bty="n"), L)
      } else {
        L1 <- modifyList(list(x=xi, y=1-pp, bty="n"), L)
      }
      L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", lty=2, xlim=c(x1, x2)))
      
      L2 <- L1[!(names(L1) %in% c("errbar.tick", "errbar.lwd"
                                  , "errbar.lty", "errbar.col"
                                  , "errbar.y.polygon"
                                  , "errbar.y.polygon.list"))]
      par(xpd=FALSE)
      a <- do.call(plot, L2) 
      
    }
  }
  return(invisible(out_sr))
}