R/plot.HatchingSuccess.R

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

#' plot.HatchingSuccess plot result of HatchingSuccess.fit() or HatchingSuccess.MHmcmc() that best describe hatching success
#' @title Plot results of HatchingSuccess.fit() that best describe hatching success
#' @author Marc Girondot
#' @return Nothing
#' @param x A result file generated by HatchingSuccess.fit()
#' @param xlim Range of temperatures
#' @param ylim Hatching success range for y-axis 
#' @param xlab x label
#' @param ylab y label
#' @param bty bty graphival parameter
#' @param las las graphical parameter
#' @param replicates Number of replicates to estimate confidence interval
#' @param resultmcmc Results obtained using HatchingSuccess.MHmcmc()
#' @param col.observations Color of observations
#' @param pch.observations Character used for observation (no observations if NULL)
#' @param cex.observations Size of characters for observations
#' @param show.CI.observations Should the confidence interval of the observations be shown ?
#' @param col.ML Color of the maximum likelihood model
#' @param lty.ML Line type of the maximum likelihood model (no line if NULL)
#' @param lwd.ML Line width of the maximum likelihood model
#' @param col.median Color of the median model
#' @param lty.median Line type of the median model (no line if NULL)
#' @param lwd.median Line width of the mean model
#' @param col.CI Color of the 95\% confidence interval lines
#' @param lty.CI Line type of the 95\% confidence interval lines (no line if NULL)
#' @param lwd.CI Line width of the 95\% confidence interval lines
#' @param polygon If TRUE, confidence interval is shown as a polygon
#' @param color.polygon The color used for polygon
#' @param what Indicate what to plot: "observations", "ML", "CI"
#' @param ... Parameters for plot()
#' @description Plot the estimates that best describe hatching success.\cr
#' If replicates is 0, it returns only the fitted model.\cr
#' If replicates is null and resultmcmc is not null, it will use all the mcmc data.\cr
#' if replicates is lower than the number of iterations in resultmcmc, it will use sequence of data regularly thined.
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' totalIncubation_Cc <- subset(DatabaseTSD, 
#'                              Species=="Caretta caretta" & 
#'                                Note != "Sinusoidal pattern" & 
#'                                !is.na(Total) & Total != 0)
#' 
#' par <- c(S.low=0.5, S.high=0.3, 
#'          P.low=25, deltaP=10, MaxHS=0.8)
#'          
#' HatchingSuccess.lnL(par=par, data=totalIncubation_Cc)
#' 
#' g <- HatchingSuccess.fit(par=par, data=totalIncubation_Cc)
#' 
#' HatchingSuccess.lnL(par=g$par, data=totalIncubation_Cc)
#' 
#' plot(g, replicates=0)
#' plot(g, replicates=10000)
#' 
#' pMCMC <- HatchingSuccess.MHmcmc_p(g, accept=TRUE)
#' mcmc <- HatchingSuccess.MHmcmc(result=g, parameters = pMCMC, 
#'                                 adaptive=TRUE, n.iter=100000, trace=1000)
#' plot(g, resultmcmc=mcmc)
#' plot(g, resultmcmc=mcmc, pch.observations=NULL, lty.mean=NULL)
#' }
#' @family Hatching success
#' @method plot HatchingSuccess
#' @export



plot.HatchingSuccess <- function(x, xlim=c(20, 40), ylim=c(0, 1), 
                                 xlab="Constant incubation temperatures", 
                                 ylab="Hatching success", bty="n", las=1, 
                                 col.observations="red", pch.observations=19, cex.observations=1,
                                 show.CI.observations=TRUE, 
                                 col.ML="black", lty.ML=1, lwd.ML=1,
                                 col.median="black", lty.median=2, lwd.median=1,
                                 col.CI="black", lty.CI=3, lwd.CI=1,
                                 replicates=NULL, resultmcmc=NULL, 
                                 polygon=TRUE, color.polygon=rgb(red = 0.8, green = 0, blue=0, alpha = 0.1), 
                                 what=c("observations", "ML", "CI"), ...) {
  
  # xlim=c(20, 40); ylim=c(0, 1)
  # xlab="Constant incubation temperatures"
  # ylab="Hatching success"; bty="n"; las=1 
  # col.observations="red"; pch.observations=19; cex.observations=1
  # show.CI.observations=TRUE
  # col.ML="black"; lty.ML=1; lwd.ML=1
  # col.median="black"; lty.median=2; lwd.median=1
  # col.CI="black"; lty.CI=3; lwd.CI=1
  # replicates=NULL; resultmcmc=NULL
  # polygon=TRUE; color.polygon=rgb(red = 0.8, green = 0, blue=0, alpha = 0.1) 
  # what=c("observations", "ML", "CI")
  # p3p <- list()
  
  p3p <- list(...)
  
  what <- tolower(what)
  what <- match.arg(what, choices = c("observations", "ml", "ci"), several.ok = TRUE)
  
  temperatures <- seq(from=xlim[1], to=xlim[2], length.out = 100)
  CIq <- predict(x, temperature=temperatures, replicates=0)
  
  if ((!is.null(lty.ML)) & any(what %in% "ml")) {
    type.ML <- "l"
  } else {
    type.ML <- "n"
  }
  
  
  p3p <- modifyList(p3p, list(x=temperatures, 
                              y=HatchingSuccess.model(par=c(x$par, x$fixed.parameters), temperature = temperatures), 
                              xlim=xlim, type=type.ML, xlab=xlab, ylab=ylab, bty=bty, las=las, 
                              ylim=ylim, col=col.ML, lty=lty.ML, lwd=lwd.ML))
  
  do.call(plot, args=p3p)
  
  if (!is.null(pch.observations) & any(what %in% "observations")) {
    
    points(x = x$data[, x$column.Incubation.temperature], 
           y = x$data[, x$column.Hatched]/(x$data[, x$column.Hatched]+x$data[, x$column.NotHatched]), 
           col=col.observations, pch=pch.observations, cex=cex.observations)
    if (show.CI.observations) {
      ciO <- getFromNamespace(x=".BinomialConfidence", ns="HelpersMG")(x=x$data[, x$column.Hatched], 
                                                                       n=x$data[, x$column.Hatched]+x$data[, x$column.NotHatched])
      for (i in 1:length(x$data[, x$column.Incubation.temperature])) {
        col.observations <- rep(col.observations, length(x$data[, x$column.Incubation.temperature]))[1:length(x$data[, x$column.Incubation.temperature])]
        segments(x0=x$data[i, x$column.Incubation.temperature], 
                 x1=x$data[i, x$column.Incubation.temperature], 
                 y0=ciO[i, 2], y1=ciO[i, 3], col=col.observations[i])
        tick <- ScalePreviousPlot()$xlim["range"]/60
        segments(x0=x$data[i, x$column.Incubation.temperature]-tick, 
                 x1=x$data[i, x$column.Incubation.temperature]+tick, 
                 y0=ciO[i, 2], y1=ciO[i, 2], col=col.observations[i])
        segments(x0=x$data[i, x$column.Incubation.temperature]-tick, 
                 x1=x$data[i, x$column.Incubation.temperature]+tick, 
                 y0=ciO[i, 3], y1=ciO[i, 3], col=col.observations[i])
      }
    }
  }
  
  
  show.CI <- TRUE
  if (any(what %in% "ci")) {
    if (is.null(resultmcmc)) {
      if (is.null(replicates)) {
        show.CI <- FALSE
      } else {
        if (replicates == 0) show.CI <- FALSE
      }
    }
  } else {
    show.CI <- FALSE
  }
  
  if (show.CI) {
    CIq <- predict(x, temperature=temperatures, replicates=replicates, resultmcmc=resultmcmc)
    
    if ((polygon) & (!is.null(lty.CI))) {
      polygon(x=c(temperatures, rev(temperatures)), y=c(CIq["2.5%", ], rev(CIq["97.5%", ])), 
              lty=lty.CI, lwd=lwd.CI, col=color.polygon)
    } else {
      if (!is.null(lty.CI)) {
        lines(x = temperatures, y=CIq["2.5%", ], lty=lty.CI, lwd=lwd.CI)
        lines(x = temperatures, y=CIq["97.5%", ], lty=lty.CI, lwd=lwd.CI)
      }
    }
    if (!is.null(lty.median)) {
      lines(x = temperatures, y=CIq["50%", ], lty=lty.median, lwd=lwd.median)
    }
  }
  
  # return(invisible())
}