Nothing
#' 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.