R/MovingIncubation.R

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

#' MovingIncubation simulate incubation of a nest with the beginning varying day by day
#' @title Simulate incubation of a nest with the beginning of incubation varying
#' @author Marc Girondot
#' @return A dataframe with informations about thermosensitive period length and incubation length day by day of incubation
#' @param NestsResult A result file generated by searchR
#' @param resultmcmc A mcmc result. Will be used rather than SE if provided.
#' @param GTRN.CI How to estimate CI for embryo growth thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
#' @param SexualisationTRN A model for sexualisation thermal reaction norm during TSP obtained using STRN()
#' @param SexualisationTRN.mcmc MCMC object for STRN.
#' @param SexualisationTRN.CI How to estimate CI of sexualisation thermal reaction norm. Can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
#' @param tsd A object from tsd() that describe the thermal react norm of sex ratio at constant temperatures
#' @param tsd.mcmc A object from tsd_MHmcmc() 
#' @param tsd.CI How to estimate CI for sex ratio thermal reaction norm; Can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
#' @param temperatures.df A data.frame with 2 or 3 columns: Times, Temperatures and Temperatures.end.incubation (facultative)
#' @param temperature.heterogeneity SD of heterogeneity of temperatures. Can be 2 values, sd_low and sd_high and then HelpersMG::r2norm() is used.
#' @param average.incubation.duration The average time to complete incubation (not used if metabolic heating is setup)
#' @param metabolic.heating Degrees Celsius to be added at the end of incubation due to metabolic heating
#' @param max.time Maximum time of incubation
#' @param skip Number of data to skip between two runs
#' @param parameters A set of parameters if result is not provided.
#' @param fixed.parameters Another set of parameters if result is not provided.
#' @param SE Standard error for each parameter if not present in result is not provided
#' @param hessian A hessian matrix
#' @param integral Function used to fit embryo growth: integral.Gompertz, integral.exponential or integral.linear
#' @param derivate Function used to fit embryo growth: dydt.Gompertz, dydt.exponential or dydt.linear. It will replace the one in NestsResult.
#' @param hatchling.metric Mean and SD of size of hatchlings as a vector ie hatchling.metric=c(Mean=xx, SD=yy)
#' @param M0 Measure of hatchling size proxi at laying date
#' @param TSP.borders The limits of TSP
#' @param TSP.begin Where TSP begin during the stage of beginning? In relative proportion of the stage.
#' @param TSP.end Where TSP begin during the stage of ending? In relative proportion of the stage.
#' @param embryo.stages The embryo stages. At least TSP.borders stages must be provided to estimate TSP length
#' @param replicate.CI Number of randomizations to estimate CI
#' @param parallel Should parallel computing be used. TRUE or FALSE
#' @param progressbar Should a progress bar be shown ? TRUE or FALSE
#' @description Simulate incubation of a nest with the beginning varying day by day\cr
#' Temperatures must be in a data.frame with one column (Time) being the time and the second the temperatures (Temperature). A third columns can indicate the temperature at the end of incubation (Temperature.end.incubation). Do not use FormatNests() for this dataframe.
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' data(resultNest_4p_SSM)
#' ti <- seq(from=0, to=(60*24*100), by=60)
#' temperatures <- rnorm(length(ti), 29, 5)
#' temperatures <- temperatures+ti/(60*24*100)/2
#' layout(mat=1:3)
#' parpre <- par(mar=c(4, 4, 1, 1)+0.4)
#' plot(ti/(60*24), temperatures, type="l", xlab="Days", 
#'      ylab=expression("Nest temperature in "*degree*"C"), bty="n", las=1)
#' # The sexualisation thermal reaction norm is calculated for South Pacific RMU
#' out <- MovingIncubation(NestsResult=resultNest_4p_SSM, 
#'      temperatures.df=data.frame(Time=ti, Temperature=temperatures),
#'      metabolic.heating = 0, 
#'      SexualisationTRN = structure(c(71.922411148397, 613.773055147801, 
#'      318.059753164125, 120.327257089974), 
#'      .Names = c("DHA", "DHH", "T12H", "Rho25")))
#' with(out, plot(Time/(60*24), Incubation.length.mean/(60*24), 
#'      xlab="Days along the season", 
#'      ylab="Incubation duration", 
#'      type="l", bty="n", las=1, ylim=c(70, 80)))
#' with(out, plot(Time/(60*24), TSP.GrowthWeighted.STRNWeighted.temperature.mean, 
#'      xlab="Days along the season", 
#'      ylab=expression("CTE for sex ratio in "*degree*"C"), 
#'       type="l", bty="n", las=1, ylim=c(30, 31)))
#' par(mar=parpre)
#' layout(mat=c(1))
#' }
#' @export


MovingIncubation <-
  function(NestsResult=NULL                                       , 
           resultmcmc=NULL                                        ,
           GTRN.CI="Hessian"                                      ,
           tsd=NULL                                               , 
           tsd.CI=NULL                                            ,
           tsd.mcmc=NULL                                          , 
           SexualisationTRN=NULL                                  ,
           SexualisationTRN.CI="Hessian"                          ,
           SexualisationTRN.mcmc=NULL                             ,
           temperatures.df=stop("A data.frame must be provided")  ,
           temperature.heterogeneity=0                            , 
           metabolic.heating=0                                    , 
           average.incubation.duration=60*1440                    ,
           max.time=100*24*60                                     ,
           skip = 1                                               , 
           parameters=NULL                                        ,
           fixed.parameters=NULL                                  , 
           SE=NULL                                                , 
           hessian = NULL                                         , 
           integral=NULL                                          , 
           derivate=NULL                                          , 
           hatchling.metric=NULL                                  , 
           M0=NULL                                                , 
           embryo.stages="Caretta caretta.SCL"                    , 
           TSP.borders=c(21, 26)                                  , 
           TSP.begin=0                                            , 
           TSP.end=0.5                                            , 
           replicate.CI=1                                         , 
           parallel=TRUE                                          , 
           progressbar=TRUE                                       ) {
    
    # NestsResult=NULL; resultmcmc=NULL;temperatures.df=NULL; 
    # metabolic.heating=0;GTRN.CI="Hessian"; 
    # temperature.heterogeneity=0;
    # tsd=NULL; tsd.CI=NULL; tsd.mcmc=NULL
    # average.incubation.duration=60*1440; skip = 1; parameters=NULL; 
    # max.time=100*24*60;
    # fixed.parameters=NULL; SE=NULL; hessian=NULL; integral=NULL; derivate=NULL; hatchling.metric=NULL; 
    # M0=NULL; TSP.borders=c(21, 26); embryo.stages="Caretta caretta.SCL"; 
    # TSP.begin=0; TSP.end=0.5
    # SexualisationTRN=NULL; SexualisationTRN.CI="Hessian"; SexualisationTRN.mcmc=NULL
    # replicate.CI=1; parallel=TRUE
    
    # maintenant il n'est plus possible qu'il n'y ait pas de temperatures
    #  
    if (!is.null(GTRN.CI)) {
      GTRN.CI <- tolower(GTRN.CI)
      GTRN.CI <- match.arg(GTRN.CI, choices = c("null", "se", "mcmc", "pseudohessianfrommcmc", "hessian"))
    } else {
      GTRN.CI <- "null"
    }
    if (!is.null(SexualisationTRN.CI)) {
      SexualisationTRN.CI <- tolower(SexualisationTRN.CI)
      SexualisationTRN.CI <- match.arg(SexualisationTRN.CI, choices = c("null", "se", "pseudohessianfrommcmc", "mcmc", "hessian"))
    } else {
      SexualisationTRN.CI <- "null"
    }
    
    if (!is.null(tsd.CI)) {
      tsd.CI <- tolower(tsd.CI)
      tsd.CI <- match.arg(tsd.CI, choices = c("null", "se", "mcmc", "pseudohessianfrommcmc", "hessian"))
    } else {
      tsd.CI <- "null"
    }
    
    
    if (is.list(SexualisationTRN)) SexualisationTRN <- SexualisationTRN$par
    
    
    if (is.null(hatchling.metric)) {
      # si tous sont pareils, je reprends les memes
      # Correction d'un bug, rajout de [[1]] dans result$hatchling.metric["Mean"][[1]][1] 30/7/2012
      if (all(NestsResult$hatchling.metric["Mean"]==NestsResult$hatchling.metric["Mean"][[1]][1]) & all(NestsResult$hatchling.metric["SD"]==NestsResult$hatchling.metric["SD"][[1]][1])) {
        hatchling.metric <- c(Mean=NestsResult$hatchling.metric["Mean"][[1]][1], SD=NestsResult$hatchling.metric["SD"][[1]][1])
      } else {	
        stop("The size at hatching must be provided.")
      }
    }
    
    times <- temperatures.df$Time
    temperatures <- temperatures.df$Temperature
    if (metabolic.heating == 0) {
      if (is.null(temperatures.df$Temperature.end.incubation)) {
        temperatures.end.incubation <- temperatures
      } else {
        temperatures.end.incubation <- temperatures.df$Temperature.end.incubation
      }
    } else {
      temperatures.end.incubation <- temperatures + metabolic.heating
    }
    
    if (is.null(integral)) integral <- NestsResult$integral
    if (is.null(derivate)) derivate <- NestsResult$derivate
    if (is.null(M0))	M0 <- NestsResult$M0
    if (is.null(fixed.parameters)) fixed.parameters <- NestsResult$fixed.parameters
    if (is.null(parameters)) parameters <- NestsResult$par
    
    if (is.null(SexualisationTRN)) SexualisationTRN <- NestsResult$SexualisationTRN
    
    # je peux indiquer des SE en plus de ceux de result
    if (is.null(hessian)) {
      hessian <- NestsResult$hessian
    }
    if (is.null(SE)) SE <- NestsResult$SE
    

    if (is.null(resultmcmc) & (GTRN.CI=="mcmc")) {
      GTRN.CI <- "null"
    }
    if ((is.null(SE)) & (GTRN.CI == "se")) {
      GTRN.CI <- "null"
    }
    if ((is.null(hessian)) & (GTRN.CI == "hessian")) {
      GTRN.CI <- "null"
    }
    
    # if (GTRN.CI=="null") {
    #   GTRN.CI <- "null"
    # }
    
    
    # if(all(is.na(SE))) replicate.CI <- 1
    
    nbtp <- length(temperatures)
    
    result.out <- universalmclapply(seq(from=1, to=nbtp-2, by=skip), 
                                    FUN=function(temp) {
                                      # temp <-  1
                                      # temp <-  25
                                      # temp <-  49
                                      # temp <-  73
                                      # temp <-  97
                                      # temp <-  121
                                      # temp <-  145
                                      
                                      if (!progressbar) print(temp)
                                      dt <- floor(as.numeric(times[temp:nbtp]-times[temp]))
                                      tempencours <- temperatures[temp:nbtp]
                                      temperatures.end.incubation <- temperatures.end.incubation[temp:nbtp]
                                      
                                      if (!is.null(max.time)) {
                                        maxdt <- which(dt > max.time)
                                        if (!identical(maxdt, integer(0))) {
                                          maxdt <- min(maxdt)
                                          dt <- dt[1:(maxdt-1)]
                                          tempencours <- tempencours[1:(maxdt-1)]
                                          temperatures.end.incubation  <- temperatures.end.incubation[1:(maxdt-1)]
                                        }
                                      }
                                      
                                      if ((metabolic.heating == 0) & !is.null(temperatures.end.incubation)) {
                                        # et dans dt j ai le temps
                                        tempencours <- tempencours + (temperatures.end.incubation-tempencours)*(cumsum(dt)/average.incubation.duration)
                                      }
                                      
                                      
                                      formated <- FormatNests(data.frame(Time=dt, temp=tempencours))
                                      
                                      out.incubation <- info.nests(x=NULL                                                     , 
                                                                   NestsResult=NULL                                           , 
                                                                   parameters=parameters                                      , 
                                                                   fixed.parameters = fixed.parameters                        ,
                                                                   SE=SE                                                      , 
                                                                   temperatures = formated                                    , 
                                                                   integral=integral                                          , 
                                                                   derivate=derivate                                          , 
                                                                   tsd=tsd                                                    , 
                                                                   tsd.mcmc=tsd.mcmc                                          , 
                                                                   tsd.CI=tsd.CI                                              , 
                                                                   hatchling.metric=hatchling.metric                          , 
                                                                   metric.end.incubation = "hatchling.metric"                 , 
                                                                   stop.at.hatchling.metric = TRUE                            , 
                                                                   series="all"                                               ,
                                                                   M0=M0                                                      , 
                                                                   fill=NULL                                                  ,
                                                                   TSP.borders=TSP.borders                                    , 
                                                                   TSP.begin=TSP.begin                                        ,  
                                                                   TSP.end=TSP.end                                            , 
                                                                   GTRN.CI=GTRN.CI                                            , 
                                                                   hessian = hessian                                          , 
                                                                   SexualisationTRN.CI=SexualisationTRN.CI                    , 
                                                                   embryo.stages=embryo.stages                                , 
                                                                   metabolic.heating=metabolic.heating                        , 
                                                                   temperature.heterogeneity=temperature.heterogeneity        , 
                                                                   SexualisationTRN=SexualisationTRN                          , 
                                                                   SexualisationTRN.mcmc=SexualisationTRN.mcmc                , 
                                                                   progressbar = FALSE                                        ,  
                                                                   replicate.CI=replicate.CI                                  , 
                                                                   out="summary"                                              , 
                                                                   warnings=FALSE                                             , 
                                                                   resultmcmc=resultmcmc                                      , 
                                                                   parallel=FALSE                                             )$summary                                          
                                      
                                      # metric.end.incubation=c(Temp=hatchling.metric$Mean),
                                      return(cbind(Time=times[temp], out.incubation, row.names=NULL))
                                    }, 
                                    mc.cores = ifelse(parallel, detectCores(), 1), 
                                    clusterEvalQ=list(expression(library(embryogrowth))), 
                                    clusterExport=list(varlist=c("nbtp", "times", 
                                                                 "tsd", "tsd.mcmc", "tsd.CI", 
                                                                 "temperatures", "metabolic.heating", 
                                                                 "temperatures.end.incubation", "average.incubation.duration", 
                                                                 "resultmcmc", "replicate.CI", "parameters", "fixed.parameters", 
                                                                 "SE", "integral", "derivate", "hatchling.metric", "M0", 
                                                                 "TSP.borders", "embryo.stages", "progressbar", 
                                                                 "TSP.end", "TSP.begin", "SexualisationTRN.CI", 
                                                                 "GTRN.CI", "hessian", "SexualisationTRN.mcmc", 
                                                                 "temperature.heterogeneity", "SexualisationTRN", "skip"), 
                                                       envir = environment()), 
                                    progressbar = progressbar)
    
    
    df <- data.frame(matrix(unlist(result.out), ncol=177, byrow=TRUE))
    colnames(df) <- colnames(result.out[[1]])
    df[, "Time"] <- as.numeric(df[, "Time"])
    
    for (i in 3:ncol(df)) {
      df[, i] <- as.numeric(df[, i])
    }
    
    return(df)
    
  }