R/HeterogeneityNests.R

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

#' HeterogeneityNests models heterogeneity of temperatures.
#' @title Model heterogeneity of temperatures.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Nothing.
#' @param nests An object of class Nests, Nests2, NestsResult, or mcmcComposite.
#' @param control.legend.total A list of options for legend.
#' @param control.legend.metabolicheating A list of options for legend.
#' @param show.full.incubation Show the plot with full incubation?
#' @param show.first.half.incubation Show the plot with first half incubation?
#' @param col.absolute.difference.full The color for absolute difference of temperatures for full incubation.
#' @param pch.absolute.difference.full The pch for absolute difference of temperatures for full incubation.
#' @param col.absolute.difference.first The color for absolute difference of temperatures for first-half incubation.
#' @param pch.absolute.difference.first The pch for absolute difference of temperatures for first-half incubation.
#' @param show.se The standard error of the modelled heterogeneity.
#' @param show.sd The standard deviation of the observed heterogeneity.
#' @param probs A vector of two values to set the range of quantiles to define heterogeneity. Use 0 and 1 for min and max.
#' @param n.iter Number of iteration for MCMC
#' @param n.adapt Number of adaptation for MCMC
#' @param thin Number of thin for MCMC
#' @param adaptive Do adaptive MCMC should be used
#' @param return.mcmc Should the mcmc object be returned?
#' @param trace Show intermediate results.
#' @param rules The rules for max and min for each parameters.
#' @param fitted.parameters The set of fitted parameters with plausible values based on rules.
#' @param fixed.parameters The values of parameters that will not be changed during optimisation.
#' @param ... Parameters used for plot.
#' @description Generate a model of heterogeneity of temperatures.\cr
#' @examples 
#' \dontrun{
#' library(embryogrowth)
#' data(nest)
#' Laying.Time <- matrix(c("DY.1", "15/05/2010", 
#'                  "DY.17", "24/05/2010", 
#'                  "DY.16", "24/05/2010", 
#'                  "DY.18", "25/05/2010", 
#'                  "DY.20", "25/05/2010", 
#'                  "DY.21", "26/05/2010", 
#'                  "DY.22", "26/05/2010", 
#'                  "DY.23", "26/05/2010", 
#'                  "DY.24", "27/05/2010", 
#'                  "DY.25", "27/05/2010", 
#'                  "DY.28", "28/05/2010", 
#'                  "DY.26", "28/05/2010", 
#'                  "DY.27", "28/05/2010", 
#'                  "DY.146", "20/06/2010", 
#'                  "DY.147", "20/06/2010", 
#'                  "DY.172", "24/06/2010", 
#'                  "DY.175", "24/06/2010", 
#'                  "DY.170", "24/06/2010", 
#'                  "DY.260", "06/07/2010", 
#'                  "DY.282", "12/07/2010", 
#'                  "DY.310", "18/07/2010", 
#'                  "DY.309", "18/07/2010", 
#'                  "DY.328", "25/07/2010", 
#'                  "DY.331", "26/07/2010"), byrow=TRUE, ncol=2)
#' tz <- OlsonNames()[grepl("Asia/Istanbul", OlsonNames())]
#' Laying.Time_f <- as.POSIXlt.character(Laying.Time[, 2], format = "%d/%m/%Y", tz=tz)
#' names(Laying.Time_f) <- Laying.Time[, 1]
#' nests <- FormatNests(data=nest, previous=NULL, col.Time="Time", 
#'                       LayingTime=Laying.Time_f, simplify=FALSE)
#' HeterogeneityNests(nests, ylim=c(0, 30))
#' }
#' @export


HeterogeneityNests <-
  function(nests=stop("An object of class Nests2, NestsResult, or mcmcComposite.")     , 
           probs = c(0.025, 0.975)                                                     ,
           control.legend.total=list()                                                 , 
           control.legend.metabolicheating=list()                                      , 
           show.full.incubation = TRUE                                                 , 
           show.first.half.incubation = TRUE                                           , 
           col.absolute.difference.full = rgb(blue=0.1, green=0.1, red=1, alpha=0.01)  ,
           pch.absolute.difference.full = 19                                           , 
           col.absolute.difference.first = rgb(blue=1, green=0.1, red=0.1, alpha=0.01) ,
           pch.absolute.difference.first = 19                                          , 
           show.se = TRUE                                                              , 
           show.sd = TRUE                                                              , 
           n.iter = 10000                                                              , 
           n.adapt = 1000                                                              , 
           thin = 5                                                                    , 
           adaptive =TRUE                                                              , 
           return.mcmc=FALSE                                                           , 
           trace = FALSE                                                               , 
           rules = rbind(data.frame(Name="min", Min=0.01, Max=2),
                         data.frame(Name="max", Min=0.02, Max=20),
                         data.frame(Name="S", Min=0, Max=10),
                         data.frame(Name="P", Min=0, Max=50),
                         data.frame(Name="asd", Min=0.01, Max=2),
                         data.frame(Name="bsd", Min=0.01, Max=1))                          ,
           fitted.parameters = c('min' = 1,
                      'max' = 3.95,
                      'S' = 0.62,
                      'P' = 5.70,
                      'asd' = 0.014,
                      'bsd' = 1.19)                                                     ,
           fixed.parameters = NULL                                                      ,
           ...                                                                         ) {
    
    p3p <- tryCatch(list(...), error=function(e) list())
    if (!is.null(p3p)) 
      if (!is.null(p3p$xlim)) 
        p3p$xlim[1] <- 2
    
    Min <- rules[, "Min"]; Max <- rules[, "Max"]; names(Max) <- names(Min) <- rules[, "Name"]
    
    Min <- Min[names(fitted.parameters)]
    Max <- Max[names(fitted.parameters)]
    
    if (inherits(nests, "mcmcComposite")) {
      nests <- nests$parametersMCMC$control$temperatures
    }
    
    if (inherits(nests, "Nests")) {
      nests <- UpdateNests(nests)
    }
    
    if (inherits(nests, "NestsResult")) {
      nests <- nests$data
    }
    
    # Model
    modelH <- function(par, data, fixed.parameters=NULL) {
      par <- c(par, fixed.parameters)
      theoric <- par["min"]+(par["max"]-par["min"])*(1/(1+exp((1/par["S"])*(par["P"]-data[ ,"nb"]))))
      sd <- abs(par["asd"])*theoric+abs(par["bsd"])
      sd <- ifelse(sd == 0, 1E-3, sd)
      # print(par)
      # print(theoric)
      # print(sd)
      # if (any(abs(par["asd"])*theoric+abs(par["bsd"]) <= 0)) return(1E9)
      l <- sum(-dnorm(x=data[ ,"H"], mean=theoric, sd=sd, log = TRUE))
      # print(l)
      return(l)
    }
    
    if (!identical(probs, c(0.025, 0.975))) stop("Only 0.025 and 0.975 quantiles are allowed")
    qnames <- names(quantile(c(0, 0, 0), probs = probs))
    
    name_nests <- nests$Names
    
    units <- lapply(nests$Nests, FUN=function(x) unname(x$UnitTime))
    units <- do.call("c", units)
    
    if (any(units != "minutes")) stop("Only time units in minutes is managed at this time.")
    
    LayingTime <- lapply(nests$Nests, FUN=function(x) unname(as.POSIXct(x$LayingTime)))
    LayingTime <- do.call("c", LayingTime)
    
    if (is.null(LayingTime)) {
      stop("Heterogeneity.Nests requires LayingTime being defined in FormatNests()")
    }
    
    BeginTime <- as.POSIXlt(min(LayingTime))
    BeginTime$mday <- 1
    BeginTime$mon <- 0
    BeginTime$hour <- 0
    BeginTime$min <- 0
    BeginTime$sec <- 0
    
    EndTime <- lapply(nests$Nests, FUN=function(x) unname(as.POSIXct(x$LayingTime) + max(x$data[, "Time"])*ifelse(x$UnitTime == "minutes", 60, NA)))
    EndTime <- do.call("c", EndTime)
    
    EndTime <- as.POSIXlt(max(EndTime))
    EndTime$mday <- 1
    EndTime$mon <- 0
    EndTime$hour <- 0
    EndTime$min <- 0
    EndTime$sec <- 0
    
    
    dates <- lapply(name_nests, FUN = function(x) as.POSIXct(nests$Nests[[x]]$data[, "Time"] * ifelse(nests$Nests[[x]]$UnitTime == "minutes", 60, NA) + LayingTime[x]))
    dates <- do.call("c", dates)
    
    dates <- sort(unique(dates))
    tempsmatrix <- matrix(data=NA, nrow=length(name_nests), ncol = length(dates))
    
    BeginTime <- as.POSIXct(BeginTime)
    EndTime <- as.POSIXct(EndTime)
    
    psdprop <- c('min' = 0.052968095149530374, 
                 'max' = 0.25326307291516775, 
                 'S' = 0.47408803190847471, 
                 'P' = 0.58011009250310686, 
                 'asd' = 0.033787381019768449, 
                 'bsd' = 0.034968838762981208)[names(fitted.parameters)]
    
    
    if (show.full.incubation) {
      
      for (j in name_nests) {
        Eret <- nests$Nests[[j]]$data
        d <- nests$Nests[[j]]$LayingTime + Eret[, "Time"]*ifelse(nests$Nests[[j]]$UnitTime == "minutes", 60, NA)
        tempsmatrix[which(j == name_nests), match(d, dates)] <- Eret[, "Temperatures C"]
      }
      
        qt_temp <- apply(tempsmatrix, MARGIN = 2, FUN = function(x) quantile(x, probs=c(probs[1], 0.5, probs[2]), na.rm=TRUE))
        
        qt_median_ma <- qt_temp[2, , drop=TRUE]
        qt_low_ma <- qt_temp[1, , drop=TRUE]
        qt_high_ma <- qt_temp[3, , drop=TRUE]
 
      nb_temp <- apply(tempsmatrix, MARGIN = 2, FUN = function(x) length(x[!is.na(x)]))
      
      dta <- data.frame(nb=nb_temp, H=qt_high_ma-qt_low_ma)
      dta <- dta[dta$nb > 1, ]
      # plot(dta$nb, dta$H)
      
      mH <- optim(par=fitted.parameters,
                  fn=modelH, data=dta, hessian = TRUE,
                  method = "L-BFGS-B", 
                  lower=Min,
                  upper=Max,
                  fixed.parameters=fixed.parameters,
                  control = list(maxit=5000, trace=0))
      par <- mH$par
      
      if (trace) print(par)
      
      # MCMC estimation of parameters
      priors <- setPriors(
        par = par,
        se = suppressWarnings(SEfromHessian(mH$hessian)),
        density = "dunif",
        rules = rules[rules$Name %in% names(fitted.parameters), ],
        silent = TRUE
      )
      
      priors$SDProp <- psdprop
      
      mcmc <- MHalgoGen(likelihood = modelH, parameters=priors, parameters_name = "par",
                        n.adapt = n.adapt, thin = thin, n.iter = n.iter,
                        data=dta, fixed.parameter=fixed.parameters,  trace = FALSE, adaptive = adaptive)
      
      if (adaptive) psdprop <- mcmc$parametersMCMC$SDProp.end
      # se <- SEfromHessian(mH$hessian)
      
      maxx <- max(c(max(nb_temp), p3p$xlim[2]))
      x <- seq(from=2, to=maxx, by=0.1)
      
      outpar <- RandomFromHessianOrMCMC(
        mcmc = mcmc,
        chain = 1,
        regularThin = TRUE,
        MinMax = NULL,
        fitted.parameters = NULL,
        fixed.parameters = NULL,
        method = "MCMC",
        replicates = floor(n.iter/thin),
        probs = c(0.025, 0.5, 0.975),
        fn = function(par) {y <- (par["min"]+(par["max"]-par["min"])*(1/(1+exp((1/par["S"])*(par["P"]-x)))));
        return(rnorm(length(y), y,  abs(par["asd"]) *y + abs(par["bsd"]) ))},
        silent = TRUE,
        ParTofn = "par"
      )
      
      outpar_se <- RandomFromHessianOrMCMC(
        mcmc = mcmc,
        chain = 1,
        regularThin = TRUE,
        MinMax = NULL,
        fitted.parameters = NULL,
        fixed.parameters = NULL,
        method = "MCMC",
        replicates = floor(n.iter/thin),
        probs = c(0.025, 0.5, 0.975),
        fn = function(par) {y <- (par["min"]+(par["max"]-par["min"])*(1/(1+exp((1/par["S"])*(par["P"]-x)))));
        return(y)},
        silent = TRUE,
        ParTofn = "par"
      )
      
      if (trace) print(as.parameters(mcmc))
      
      mcmc_tot <- mcmc
      
      par(mar=c(4, 4, 1, 1))
      
      param.plot <- modifyList(modifyList(list(x=dta$nb, y=dta$H, xlim=c(2, max(dta$nb)), ylim=c(0, 8),
                                               pch=19, xlab="Number of overlapping timeseries", ylab = expression("Range of temperatures in "*~degree*C*""),
                                               bty="n", las=1, col=rgb(blue=0.1, green=0.1, red=1, alpha=0.01), xaxt="n"), 
                                          p3p), list(type = "n"))
      
      do.call("plot", param.plot)
      axis(1, at=2:maxx)
      
      points(x=dta$nb, y=dta$H, col=col.absolute.difference.full, pch=pch.absolute.difference.full)
      
      # lines(x, outpar$quantiles["50%", ])
      if (show.sd)
        polygon(x=c(x, rev(x)), y=c(outpar$quantiles[qnames[1], ], rev(outpar$quantiles[qnames[2], ])),
                border = NA, col=rgb(red = 1, green = 0, blue = 0, alpha = 0.2))
      if (show.se)
        polygon(x=c(x, rev(x)), y=c(outpar_se$quantiles[qnames[1], ], rev(outpar_se$quantiles[qnames[2], ])),
                border = NA, col=rgb(red = 1, green = 0, blue = 0, alpha = 0.6))
      
      segments(x0=1, x1=maxx, y0=summary(mcmc)$quantiles["max", "50%"], y1=summary(mcmc)$quantiles["max", "50%"], col="red", lty=4, lwd=2)
      segments(x0=1, x1=maxx, y0=summary(mcmc)$quantiles["max", qnames[1]], y1=summary(mcmc)$quantiles["max", qnames[1]], col="red", lty=3, lwd=2)
      segments(x0=1, x1=maxx, y0=summary(mcmc)$quantiles["max", qnames[2]], y1=summary(mcmc)$quantiles["max", qnames[2]], col="red", lty=3, lwd=2)
      cat("Taking into account all data in time series\n")
      cat("-------------------------------------------\n")
      cat(paste("The model converges to an heterogeneity of temperatures of",
                specify_decimal(unname(summary(mcmc)$quantiles["max", "50%"]), decimals=2),
                "\u00B0C (95% Credible Interval between",
                specify_decimal(unname(summary(mcmc)$quantiles["max", qnames[1]]), decimals=2),
                "and", specify_decimal(unname(summary(mcmc)$quantiles["max", qnames[2]]), decimals=2),
                ") at the scale of the beach.\n"))
      cat(paste("The 95% maximum predicted heterogeneity is",
                specify_decimal(unname(summary(mcmc)$quantiles["bsd", "50%"])+unname(summary(mcmc)$quantiles["asd", "50%"])*unname(summary(mcmc)$quantiles["max", "50%"]) * 1.96 + unname(summary(mcmc)$quantiles["max", "50%"]), decimals=2),
                "\u00B0C.\n"))
      
      pari <- as.parameters(mcmc, index="median")
      
    } else {
      pari <- fitted.parameters
    }
    
    if (show.first.half.incubation) {
      
      # With only first half of timeseries
      tempsmatrix[] <- NA
      for (j in name_nests) {
        Eret <- nests$Nests[[j]]$data
        maxt <- max(Eret[, "Time"])/2
        Eret <- Eret[Eret[, "Time"] < maxt, ]
        d <- nests$Nests[[j]]$LayingTime + Eret[, "Time"]*60
        
        tempsmatrix[which(j == name_nests), match(d, dates)] <- Eret[, "Temperatures C"]
      }
      
      qt_temp <- apply(tempsmatrix, MARGIN = 2, FUN = function(x) quantile(x, probs=c(probs[1], 0.5, probs[2]), na.rm=TRUE))

      qt_median_ma <- qt_temp[2, , drop=TRUE]
      qt_low_ma <- qt_temp[1, , drop=TRUE]
      qt_high_ma <- qt_temp[3, , drop=TRUE]
      
      
      nb_temp <- apply(tempsmatrix, MARGIN = 2, FUN = function(x) length(x[!is.na(x)]))
      
      dta <- data.frame(nb=nb_temp, H=qt_high_ma-qt_low_ma)
      dta <- dta[dta$nb > 1, ]
      
      
      mH <- optim(par=pari,
                  fn=modelH, data=dta, hessian = TRUE,
                  fixed.parameters=fixed.parameters,
                  method = "L-BFGS-B", 
                  lower=Min,
                  upper=Max,
                  control = list(maxit=5000))
      par <- mH$par
      
      if (trace) print(par)
      
      # MCMC estimation of parameters
      priors <- setPriors(
        par = par,
        se = suppressWarnings(SEfromHessian(mH$hessian)),
        density = "dunif",
        rules = rules[rules$Name %in% names(fitted.parameters), ],
        silent = TRUE
      )
      
      
      
      priors$SDProp <- psdprop
      
      mcmc <- MHalgoGen(likelihood = modelH, parameters=priors, parameters_name = "par",
                        n.adapt = n.adapt, thin = thin, n.iter = n.iter,
                        fixed.parameters=fixed.parameters,
                        data=dta, trace = FALSE, adaptive = adaptive)
      
      # se <- SEfromHessian(mH$hessian)
      
      if (trace) print(as.parameters(mcmc))
      
      maxx <- max(c(max(nb_temp), p3p$xlim[2]))
      x <- seq(from=2, to=maxx, by=0.1)
      
      outpar <- RandomFromHessianOrMCMC(
        mcmc = mcmc,
        chain = 1,
        regularThin = TRUE,
        MinMax = NULL,
        fitted.parameters = NULL,
        fixed.parameters = fixed.parameters,
        method = "MCMC",
        replicates = floor(n.iter/thin),
        probs = c(0.025, 0.5, 0.975),
        fn = function(par) {y <- (par["min"]+(par["max"]-par["min"])*(1/(1+exp((1/par["S"])*(par["P"]-x)))));
        return(rnorm(length(y), y, abs(par["asd"]) * y + abs(par["bsd"])))},
        silent = TRUE,
        ParTofn = "par"
      )
      
      outpar_se <- RandomFromHessianOrMCMC(
        mcmc = mcmc,
        chain = 1,
        regularThin = TRUE,
        MinMax = NULL,
        fitted.parameters = NULL,
        fixed.parameters = fixed.parameters,
        method = "MCMC",
        replicates = floor(n.iter/thin),
        probs = c(0.025, 0.5, 0.975),
        fn = function(par) {y <- (par["min"]+(par["max"]-par["min"])*(1/(1+exp((1/par["S"])*(par["P"]-x)))));
        return(y)},
        silent = TRUE,
        ParTofn = "par"
      )
      
      if (!show.full.incubation) {
        par(mar=c(4, 4, 1, 1))
        
        param.plot <- modifyList(modifyList(list(x=dta$nb, y=dta$H, xlim=c(2, max(dta$nb)), ylim=c(0, 8),
                                                 pch=19, xlab="Number of overlapping timeseries", ylab = expression("Range of temperatures in "*~degree*C*""),
                                                 bty="n", las=1, col=rgb(blue=0.1, green=0.1, red=1, alpha=0.01), xaxt="n"), 
                                            p3p), list(type = "n"))
        
        do.call("plot", param.plot)
        axis(1, at=2:maxx)
        
      }
      
      points(x=dta$nb, y=dta$H, 
             col=col.absolute.difference.first, pch=pch.absolute.difference.first)
      
      if (show.sd)
        polygon(x=c(x, rev(x)), y=c(outpar$quantiles[qnames[1], ], rev(outpar$quantiles[qnames[2], ])),
                border = NA, col=rgb(red = 0, green = 0, blue = 1, alpha = 0.2))
      if (show.se)
        polygon(x=c(x, rev(x)), y=c(outpar_se$quantiles[qnames[1], ], rev(outpar_se$quantiles[qnames[2], ])),
                border = NA, col=rgb(red = 0, green = 0, blue = 1, alpha = 0.6))
      
      pari <- as.parameters(mcmc, index="median")
      
      segments(x0=2, x1=maxx, y0=unname(summary(mcmc)$quantiles["max", "50%"]), y1=unname(summary(mcmc)$quantiles["max", "50%"]), col="blue", lty=4, lwd=2)
      
      segments(x0=2, x1=maxx, y0=summary(mcmc)$quantiles["max", qnames[1]], y1=summary(mcmc)$quantiles["max", qnames[1]], col="blue", lty=3, lwd=2)
      segments(x0=2, x1=maxx, y0=summary(mcmc)$quantiles["max", qnames[2]], y1=summary(mcmc)$quantiles["max", qnames[2]], col="blue", lty=3, lwd=2)
      
      cat("Taking into account only first half time series\n")
      cat("-----------------------------------------------\n")
      cat(paste("The model converges to an heterogeneity of temperatures of",
                specify_decimal(unname(summary(mcmc)$quantiles["max", "50%"]), decimals=2),
                "\u00B0C (95% Credible Interval between",
                specify_decimal(unname(summary(mcmc)$quantiles["max", qnames[1]]), decimals=2),
                "and", specify_decimal(unname(summary(mcmc)$quantiles["max", qnames[2]]), decimals=2),
                ") at the scale of the beach.\n"))
      cat(paste("The 95% maximum predicted heterogeneity is",
                specify_decimal(unname(summary(mcmc)$quantiles["bsd", "50%"])+unname(summary(mcmc)$quantiles["asd", "50%"])*unname(summary(mcmc)$quantiles["max", "50%"]) * 1.96 + unname(summary(mcmc)$quantiles["max", "50%"]), decimals=2),
                "\u00B0C.\n"))
    }
    #
    if (show.se & show.sd) {
      control.legend.total <- modifyList(list(x="bottomright", legend = c("Observations", "95% SE", "95% Bayesian CI model", "Heterogeneity of the beach"),
                                              col=c("lightgrey", "grey", "black", "black"), 
                                              pch=c(pch.absolute.difference.first, 15, NA, NA), lty=c(NA, NA, 1, 3), lwd=c(NA, NA, 3, 1)), 
                                         control.legend.total)
    } else {
      if (show.se) {
        control.legend.total <- modifyList(list(x="bottomright", legend = c("Observations", "95% Bayesian CI model", "Heterogeneity of the beach"),
                                                col=c("lightgrey", "black", "black"), 
                                                pch=c(pch.absolute.difference.first, NA, NA), lty=c(NA, 1, 3), lwd=c(NA, 3, 1)), 
                                           control.legend.total)
      } else {
        if (show.sd) {
          control.legend.total <- modifyList(list(x="bottomright", legend = c("Observations", "95% SD", "Heterogeneity of the beach"),
                                                  col=c("lightgrey", "grey", "black"), 
                                                  pch=c(pch.absolute.difference.first, 15, NA), lty=c(NA, NA, 3), lwd=c(NA, NA, 1)), 
                                             control.legend.total)
        } else {
          control.legend.total <- modifyList(list(x="bottomright", legend = c("Observations", "Heterogeneity of the beach"),
                                                  col=c("lightgrey", "black"), 
                                                  pch=c(pch.absolute.difference.first, NA), lty=c(NA, 3), lwd=c(NA, 1)), 
                                             control.legend.total)
        }
      }
    }
    do.call("legend", control.legend.total)
    
    control.legend.metabolicheating <- modifyList(list(x="topleft", legend = c("With Metabolic heating", "Without metabolic heating"), 
                                                       lty=c(1, 1), col=c("red", "blue")), 
                                                  control.legend.metabolicheating)
    
    do.call("legend", control.legend.metabolicheating)
    
    out <- NULL
    if (show.first.half.incubation & show.full.incubation)
      out <- c(lower.max.half=unname(summary(mcmc)$quantiles["max", qnames[1]]), 
               median.max.half=unname(summary(mcmc)$quantiles["max", "50%"]), 
               upper.max.half=unname(summary(mcmc)$quantiles["max", qnames[2]]), 
               lower.max=unname(summary(mcmc_tot)$quantiles["max", qnames[1]]), 
               median.max=unname(summary(mcmc_tot)$quantiles["max", "50%"]), 
               upper.max=unname(summary(mcmc_tot)$quantiles["max", qnames[2]]))
    
    if (show.first.half.incubation & !show.full.incubation)
      out <- c(lower.max.half=unname(summary(mcmc)$quantiles["max", qnames[1]]), 
               median.max.half=unname(summary(mcmc)$quantiles["max", "50%"]), 
               upper.max.half=unname(summary(mcmc)$quantiles["max", qnames[2]]))
    
    if (!show.first.half.incubation & show.full.incubation)
      out <- c(lower.max=unname(summary(mcmc_tot)$quantiles["max", qnames[1]]), 
               median.max=unname(summary(mcmc_tot)$quantiles["max", "50%"]), 
               upper.max=unname(summary(mcmc_tot)$quantiles["max", qnames[2]]))
    
    if (return.mcmc) {
      if (show.first.half.incubation & show.full.incubation)
        out <- list(summary=out, 
                    mcmc=mcmc, 
                    mcmc_tot=mcmc_tot)
      if (show.first.half.incubation & !show.full.incubation)
        out <- list(summary=out, 
                    mcmc=mcmc)
      if (!show.first.half.incubation & show.full.incubation)
        out <- list(summary=out, 
                    mcmc_tot=mcmc_tot)
    } else {
      out <- list(summary=out)
    }
    
    return(invisible(out))
  }