Nothing
R/info.nests.R
In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth
Defines functions
info.nests
Documented in
info.nests
#' Calculate statistics about nests
#' @title Calculte statistics about nests
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return or the total likelihood or a list with $metric and $summary depending on out parameter
#' @param x A set of parameters to model the embryo growth thermal reaction norm or a NestsResult object.
#' @param parameters A set of parameters to model the embryo growth thermal reaction norm. It will replace the parameters included in NestsResult (same as x).
#' @param NestsResult A NestsResult object generated by searchR to model the embryo growth thermal reaction
#' @param resultmcmc A mcmc result for embryo growth thermal reaction norm
#' @param hessian An hessian matrix for embryo growth thermal reaction norm. It will replace the hessian matrix included in NestResult object.
#' @param GTRN.CI How to estimate CI for embryo growth thermal reaction norm; can be NULL, "SE", "MCMC", or "Hessian".
#' @param fixed.parameters A set of fixed parameters to model the embryo growth thermal reaction norm. It will replace the fixed parameters included in NestsResult.
#' @param SE Standard error for each parameter. It will replace the SE in NestsResult. Use SE=NA to remove SE from NestResult
#' @param temperatures Timeseries of temperatures formatted using formatNests(). It will replace the one in NestsResult.
#' @param integral Function used to fit embryo growth: integral.Gompertz, integral.exponential or integral.linear. It will replace the one in NestsResult.
#' @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. It will replace the one in NestsResult.
#' @param weight Weights of the different nests to estimate likelihood. It will replace the ones in NestsResult.
#' @param stop.at.hatchling.metric TRUE or FALSE. If TRUE, the model stops when proxy of size reached the mean hatchling.metric size.
#' @param M0 Measure of hatchling size proxi at laying date. It will replace the one in NestsResult.
#' @param SexualisationTRN A set of parameters used to model sexualisation thermal reaction norm during TSP or a result of STRN()
#' @param SexualisationTRN.mcmc A mcmc object obtained from STRN_MHmcmc() to generate variability for sexualisation thermal reaction norm during TSP
#' @param SexualisationTRN.CI How to estimate CI of sexualisation thermal reaction norm. Can be NULL, "SE", "MCMC", or "Hessian".
#' @param series The name or number of the series to be estimated.
#' @param TSP.borders The limits of TSP in stages. See embryo.stages parameter.
#' @param embryo.stages The embryo stages. At least TSP.borders stages must be provided to estimate TSP borders. See note.
#' @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 replicate.CI Number of replicates to estimate CI. See description
#' @param metabolic.heating Degrees Celsius to be added at the end of incubation due to metabolic heating.
#' @param temperature.heterogeneity SD of heterogeneity of temperatures. Can be 2 values, sd_low and sd_high and then HelpersMG::r2norm() is used.
#' @param fill Number of minutes between two records. Create new one if they do not exist. NULL does not change the time of temperature recordings.
#' @param probs Probabilities for metric quantiles.
#' @param out Can take the values of "likelihood", "summary", "metric" or "dynamic".
#' @param metric.end.incubation The metric at the end of incubation used to calibrate TSP size. Can be "hatchling.metric", or "observed".
#' @param progressbar If FALSE, the progress bar is not shown (useful for using with sweave or knitr)
#' @param warnings If FALSE, does not show warnings
#' @param parallel If TRUE use parallel version for nests estimation
#' @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", or "Hessian".
#' @param zero Value to replace 0 or 1.
#' @param verbose If TRUE, show more information.
#' @description This function calculates many statistics about nests.\cr
#' The embryo.stages is a named vector with relative size as compared to final size at the beginning of the stage. Names are the stages.\cr
#' For example for SCL in Caretta caretta:\cr
#' embryo.stages=structure(c(8.4, 9.4, 13.6, 13.8, 18.9, 23.5, 32.2, 35.2, 35.5, 38.5)/39.33), \cr
#' .Names = c("21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31"))\cr
#' indicates that the stages 21 begins at the relative size of 8.4/39.33 as compared to the final size.\cr
#' Series can be indicated as the name of the series, their numbers or series or succession of TRUE or FALSE. "all" indicates that all series must be analyzed.\cr
#' The likelihood object is just the total likelihood of the data in the model.\cr
#' If one parameter is named "pipping_emergence" it is used as the number of days between pipping and emergence to calculate the 1/3 and 2/3 of incubation.\cr
#' The summary object is a data.frame composed of these elements with the suffix .mean, .se or .quantile_x with x from the parameter probs.
#' \itemize{
#' \item \code{Temperature.max} Maximum temperature recorded during incubation
#' \item \code{TimeWeighted.temperature} Average temperature during all incubation
#' \item \code{GrowthWeighted.temperature} Average temperature weighted by the actual growth during all incubation
#' \item \code{TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during all incubation
#' \item \code{TSP.TimeWeighted.temperature} Average temperature during the TSP
#' \item \code{TSP.GrowthWeighted.temperature} Average temperature weighted by the actual growth during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during the TSP
#' \item \code{TSP.TimeWeighted.STRNWeighted.temperature} Average temperature weighted by the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.GrowthWeighted.STRNWeighted.temperature} Average temperature weighted by actual growth and the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature} Average temperature weighted by growth rate and the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.length} TSP duration
#' \item \code{TSP.begin} Beginning of the TSP
#' \item \code{TSP.end} End of the TSP
#' \item \code{TSP.PM.GrowthWeighted} Average of male probability for each temperature weighted by actual growth during the TSP
#' \item \code{TSP.PM.TimeWeighted.GrowthRateWeighted} Average of male probability for each temperature weighted by growth rate during the TSP
#' \item \code{TSP.PM.TimeWeighted} Average of male probability for each temperature during the TSP
#' \item \code{Incubation.length} Incubation length duration
#' \item \code{MiddleThird.length} Middle third incubation duration
#' \item \code{MiddleThird.begin} Beginning of the middle third incubation duration
#' \item \code{MiddleThird.end} End of the middle third incubation duration
#' \item \code{MiddleThird.TimeWeighted.temperature} Average temperature during the middle third incubation
#' \item \code{MiddleThird.GrowthWeighted.temperature} Average temperature weighted by the actual growth during the middle third incubation
#' \item \code{MiddleThird.TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during the middle third incubation
#' \item \code{TSP.TimeWeighted.sexratio} Sex ratio based on average temperature during the TSP
#' \item \code{TSP.GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by the actual growth during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by the growth rate during the TSP
#' \item \code{TSP.TimeWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.GrowthWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the actual growth and thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the growth rate and the thermal reaction norm of sexualization during the TSP
#' \item \code{MiddleThird.TimeWeighted.sexratio} Sex ratio based on average temperature during the middle third incubation
#' \item \code{MiddleThird.GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by actual growth during the middle third incubation
#' \item \code{MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by growth rate during the middle third incubation
#' \item \code{TimeWeighted.sexratio} Sex ratio based on average temperature during all incubation
#' \item \code{GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by actual growth during all incubation
#' \item \code{TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by growth rate during all incubation
#' }
#' If \code{out} is equal to \code{summary}, the return is a list with:
#' \itemize{
#' \item \code{summary} is a data.frame with statistics for each nest.
#' \item \code{dynamic.metric} object is a list composed of data.frames with the dynamics of growth for each nest. It showed only temperatures from original dataset.
#' \item \code{summary.dynamic.metric} is a data.frame with the following columns with the suffix .mean, .se or .quantile_x with x from the parameter probs.
#' }
#' If \code{out} is equal to \code{metric}, the return is a list with:
#' \itemize{
#' \item \code{dynamic.metric} object is a list composed of data.frames with the dynamics of growth for each nest
#' \item \code{indices.dynamic.metric} is a data.frame with the following columns.
#' }
#' The object \code{summary.dynamic.metric} or \code{indices.dynamic.metric} is a data.frame with the following columns:
#' \itemize{
#' \item \code{series} Name of the series
#' \item \code{metric.begin.tsp} Metric at the beginning of TSP
#' \item \code{metric.end.tsp} Metric at the end of TSP
#' \item \code{hatchling.metric.mean} Average expected size of hatchlings
#' \item \code{hatchling.metric.sd} standard deviation of expected size of hatchlings
#' \item \code{time.begin.tsp} Time at the beginning of TSP
#' \item \code{time.end.tsp} Time at the end of TSP
#' \item \code{time.begin.middlethird} Time at the beginning of the middle third incubation
#' \item \code{time.end.middlethird} Time at the end of the middle third incubation
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' }
#' If you indicate new set of temperatures, you must probably also indicate new hatchling.metric values.\cr
#' Note: four species have predefined embryo stages. embryo.stages parameter can take the values:\cr
#' \itemize{
#' \item \code{Caretta caretta.SCL}
#' \item \code{Chelonia mydas.SCL}
#' \item \code{Emys orbicularis.SCL}
#' \item \code{Emys orbicularis.mass}
#' \item \code{Podocnemis expansa.SCL}
#' \item \code{Lepidochelys olivacea.SCL}
#' \item \code{Generic.ProportionDevelopment}
#' }
#' But remember that mass is not the best proxy to describe the growth of an embryo because it can decrease if the substrate becomes dry.\cr
#' The progress bar is based on both replicates and timeseries progress. It necessitates the pbapply package.\cr
#' If replicate.CI is null or 0, only maximum likelihood is used and no confidence interval is calculated.\cr
#' If replicate.CI is 1, one random value for the parameters is used but no confidence interval is calculated.\cr
#' In other cases, replicate.CI random samples are used to estimate confidence interval.
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' data(resultNest_4p_SSM)
#' # Some basic calculations to show the advantage of parallel computing
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=FALSE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE, progressbar=TRUE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="likelihood",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE, progressbar=FALSE))
#'
#' # By default parallel computing is TRUE but progressbar is FALSE
#' # When out is "likelihood", it returns only the likelihood
#' # otherwise, it returns a list with 3 objects "summary",
#' # "dynamic.metric", and "summary.dynamic.metric".
#'
#' summary.nests <- info.nests(resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL",
#' replicate.CI=100,
#' resultmcmc=resultNest_mcmc_4p_SSM,
#' GTRN.CI="MCMC",
#' progressbar=TRUE)
#'
#' summary.nests <- info.nests(resultNest_4p_SSM,
#' embryo.stages="Caretta caretta.SCL",
#' out="summary", replicate.CI=100,
#' GTRN.CI="Hessian",
#' progressbar=TRUE)
#'
#' summary.nests <- info.nests(resultNest_4p_SSM,
#' series = 1,
#' embryo.stages="Caretta caretta.SCL",
#' out="summary", replicate.CI=100,
#' GTRN.CI="SE",
#' progressbar=TRUE)
#'
#' # Example of use of embryo.stages and TSP.borders:
#' summary.nests <- info.nests(resultNest_4p_SSM, out="summary",
#' embryo.stages=c("10"=0.33, "11"=0.33, "12"=0.66, "13"=0.66),
#' TSP.borders = c(10, 12),
#' replicate.CI=100,
#' progressbar=TRUE)
#'
#' #########################################
#' # Sex ratio using Massey et al. method PM
#' #########################################
#'
#' # Massey, M.D., Holt, S.M., Brooks, R.J., Rollinson, N., 2019. Measurement
#' # and modelling of primary sex ratios for species with temperature-dependent
#' # sex determination. J Exp Biol 222, 1-9.
#'
#' CC_Mediterranean <- subset(DatabaseTSD, RMU=="Mediterranean" &
#' Species=="Caretta caretta" & (!is.na(Sexed) & Sexed!=0))
#' tsdL <- with (CC_Mediterranean, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature,
#' equation="logistic", replicate.CI=NULL))
#'
#' PM <- info.nests(x=resultNest_4p_SSM,
#' GTRN.CI="Hessian", tsd.CI="Hessian",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=100,
#' out="summary", progressbar=TRUE, tsd=tsdL)
#'
#'
#' plot_errbar(x=PM$summary$TimeWeighted.temperature.mean,
#' y=PM$summary$TSP.PM.GrowthWeighted.mean,
#' y.minus=PM$summary$TSP.PM.GrowthWeighted.quantile_0.025,
#' y.plus=PM$summary$TSP.PM.GrowthWeighted.quantile_0.975,
#' xlab="CTE SCL growth",
#' ylab="PM Massey et al. 2016", xlim=c(26, 32), ylim=c(0, 1), las=1)
#'
#' # Relationship between growth and growth rate
#'
#' infoall.df <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL",
#' replicate.CI=100,
#' resultmcmc=resultNest_mcmc_4p_SSM,
#' GTRN.CI="MCMC",
#' progressbar=TRUE)
#'
#' layout(1)
#' plot(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_50%"],
#' type="l", las=1, bty="n",
#' xlab="Time in minute", ylab="Growth", ylim=c(0, 39), xlim=c(0, 100000))
#' lines(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_2.5%"], lty=2)
#' lines(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_97.5%"], lty=2)
#' }
#' @export
info.nests <- function(x=NULL ,
parameters=NULL ,
NestsResult=NULL ,
resultmcmc=NULL ,
hessian=NULL ,
GTRN.CI=NULL ,
fixed.parameters=NULL ,
SE=NULL ,
temperatures=NULL ,
integral=NULL ,
derivate=NULL ,
hatchling.metric=NULL ,
stop.at.hatchling.metric=FALSE ,
M0=NULL ,
series="all" ,
TSP.borders=NULL ,
embryo.stages=NULL ,
TSP.begin=0 ,
TSP.end=0.5 ,
replicate.CI=0 ,
weight=NULL ,
out="likelihood" ,
fill=NULL ,
probs= c(0.025, 0.5, 0.975) ,
SexualisationTRN=NULL ,
SexualisationTRN.mcmc=NULL ,
SexualisationTRN.CI=NULL ,
metric.end.incubation="observed" ,
metabolic.heating=0 ,
temperature.heterogeneity=0 ,
progressbar=FALSE ,
warnings=TRUE ,
parallel=TRUE ,
tsd=NULL ,
tsd.CI=NULL ,
tsd.mcmc=NULL ,
zero = 1E-9 ,
verbose = FALSE ) {
# parameters=NULL; temperatures=NULL; fixed.parameters=NULL
# integral=NULL; hatchling.metric=NULL; derivate=NULL
# weight=NULL; M0=NULL
# x=NULL; NestsResult=NULL; resultmcmc=NULL
# SE=NULL; hessian=NULL
# GTRN.CI=NULL; temperature.heterogeneity=0; metabolic.heating=0
# stop.at.hatchling.metric=FALSE; series="all"; TSP.borders=NULL
# embryo.stages=NULL
# TSP.begin=0; TSP.end=0.5; replicate.CI=100; out="Likelihood"
# probs= c(0.025, 0.5, 0.975)
# fill=NULL; SexualisationTRN.CI=NULL
# SexualisationTRN=NULL;SexualisationTRN.mcmc=NULL
# metric.end.incubation="observed"
# progressbar=FALSE;warnings=TRUE; parallel=TRUE
# tsd=NULL; tsd.CI=NULL; tsd.mcmc=NULL
# verbose = FALSE
# NestsResult=resultNest_Florida
# series=1:156
# embryo.stages="Caretta caretta.SCL"
# TSP.borders=NULL
# TSP.begin=0
# TSP.end=0.5
# SexualisationTRN=c('DHA' = 65.329132145227575, 'DHH' = 811.18443404950051, 'T12H' = 655.01294236106037)
# out="summary"
# replicate.CI = 0
# progressbar=FALSE
# parallel=TRUE
# x=resultNest_4p_SSM; out <- "summary"; embryo.stages="fitted"; replicate.CI=0; GTRN.CI="Hessian"; parallel=TRUE
# SexualisationTRN <- c('BeginTSP' = -1e+05, 'LengthTSP' = 0.96604492361218408)
# x=resultNest_4p_SSM; out <- "summary"; embryo.stages="Caretta caretta.SCL"; replicate.CI=100; GTRN.CI="Hessian"; parallel=TRUE
metric.end.incubation <- match.arg(arg=tolower(metric.end.incubation), choices=c("hatchling.metric", "observed"))
if ((metric.end.incubation == "observed") & stop.at.hatchling.metric) {
stop("The option metric.end.incubation == 'observed' is not compatible with stop.at.hatchling.metric being TRUE")
}
SSM <- getFromNamespace(".SSM", ns="embryogrowth")
papply <- ifelse(parallel, detectCores(), 1)
if (inherits(x, "NestsResult"))
NestsResult <- x
if (inherits(x, "numeric"))
parameters <- x
out <- tolower(out)
out <- match.arg(out, choices = c("likelihood", "metric", "summary", "dynamic"), several.ok = FALSE)
TSP.list <- embryogrowth::TSP.list
if (inherits(NestsResult, "NestsResult")) {
# temperatures est un objet Nests
if (is.null(hessian)) hessian <- NestsResult$hessian
if (is.null(temperatures)) {
temperatures <- NestsResult$data
if (is.null(hatchling.metric)) hatchling.metric <- NestsResult$hatchling.metric
}
if (is.null(integral)) integral <- NestsResult$integral
if (is.null(derivate)) derivate <- NestsResult$derivate
if (is.null(weight)) weight <- NestsResult$weight
if (is.null(M0)) M0 <- NestsResult$M0
if (is.null(fixed.parameters)) fixed.parameters <- NestsResult$fixed.parameters
if (is.null(SE)) SE <- NestsResult$SE
if (is.null(parameters)) parameters <- NestsResult$par
} else {
if (is.numeric(NestsResult) & is.null(parameters)) parameters <- NestsResult
}
if (!is.null(integral)) {
if (!is.null(derivate)) {
if (warnings) warning("Both derivate and integral cannot be supplied. Only integral parameter will be used")
derivate <- NULL
}
}
GTRN.CI <- ifelse(is.null(GTRN.CI), "null", GTRN.CI)
tsd.CI <- ifelse(is.null(tsd.CI), "null", tsd.CI)
SexualisationTRN.CI <- ifelse(is.null(SexualisationTRN.CI), "null", SexualisationTRN.CI)
if (is.null(replicate.CI)) replicate.CI <- 0
# Si je ne fais que la vraisemblance, je n'ai pas besoin de CI
if (out == "likelihood") {
if (replicate.CI != 0) {
if (warnings) warning("No replicate can be used for likelihood out.")
replicate.CI <- 0
}
}
if (replicate.CI == 0) {
GTRN.CI <- "null"
SexualisationTRN.CI <- "null"
tsd.CI <- "null"
}
if ((GTRN.CI == "null") & (tsd.CI == "null") & (SexualisationTRN.CI == "null")) {
# 15/9/2019: C'est une valeur mais sans resampling
replicate.CI <- 0
}
####################################### #
# vérifie la conformité de GTRN.CI #####
####################################### #
# GTRN.CI How to estimate CI for embryo growth thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc", or "Hessian".
GTRN.CI <- tolower(GTRN.CI)
GTRN.CI <- match.arg(GTRN.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
GTRN.CIr <- GTRN.CI
if ((GTRN.CI == "mcmc" | GTRN.CI == "pseudohessianfrommcmc") & is.null(resultmcmc)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as MCMC but no resultmcmc is available.")
}
if (!is.null(SE[1])) if (all(is.na(SE[]))) SE <- NULL
if ((GTRN.CI == "se") & is.null(SE)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as SE but no SE is available.")
}
if (GTRN.CI == "hessian" & is.null(hessian)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as Hessian but no Hessian matrix is available.")
}
GTRN.CI <- GTRN.CIr
################################################### #
# vérifie la conformité de SexualisationTRN.CI #####
################################################### #
# SexualisationTRN.CI How to estimate CI for sexualization thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
SexualisationTRN.CI <- tolower(SexualisationTRN.CI)
SexualisationTRN.CI <- match.arg(SexualisationTRN.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
SexualisationTRN.CIr <- SexualisationTRN.CI
if ((SexualisationTRN.CI == "mcmc" | SexualisationTRN.CI == "pseudohessianfrommcmc") & is.null(SexualisationTRN.mcmc)) {
SexualisationTRN.CIr <- "null"
if (warnings) warning("SexualisationTRN.CI is indicated as MCMC but no SexualisationTRN.mcmc is available.")
}
if (SexualisationTRN.CI == "hessian") {
if (is.list(SexualisationTRN)) {
if (is.null(SexualisationTRN$hessian)) {
SexualisationTRN.CIr <- "null"
}
} else {
SexualisationTRN.CIr <- "null"
}
if ((SexualisationTRN.CIr == "null") & warnings) warning("SexualisationTRN.CI is indicated as Hessian but no Hessian matrix is available.")
}
if (SexualisationTRN.CI == "se") {
if (is.list(SexualisationTRN)) {
if (is.null(SexualisationTRN$SE)) {
SexualisationTRN.CIr <- "null"
}
} else {
SexualisationTRN.CIr <- "null"
}
if ((SexualisationTRN.CIr == "null") & warnings) warning("SexualisationTRN.CI is indicated as SE but no SE is available.")
}
SexualisationTRN.CI <- SexualisationTRN.CIr
###################################### #
# vérifie la conformité de tsd.CI #####
###################################### #
# tsd.CI How to estimate CI for sex ratio thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
tsd.CI <- tolower(tsd.CI)
tsd.CI <- match.arg(tsd.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
tsd.CIr <- tsd.CI
if ((tsd.CI == "mcmc" | tsd.CI == "pseudohessianfrommcmc") & is.null(tsd.mcmc)) {
tsd.CIr <- "null"
if (warnings) warning("tsd.CI is indicated as MCMC but no tsd.mcmc is available.")
}
if (tsd.CI == "hessian") {
if (is.list(tsd)) {
if (is.null(tsd$hessian)) {
tsd.CIr <- "null"
}
} else {
tsd.CIr <- "null"
}
if ((tsd.CIr == "null") & warnings) warning("tsd.CI is indicated as Hessian but no Hessian matrix is available.")
}
if (tsd.CI == "se") {
if (is.list(tsd)) {
if (is.null(tsd$SE)) {
tsd.CIr <- "null"
}
} else {
tsd.CIr <- "null"
}
if ((tsd.CIr == "null") & warnings) warning("tsd.CI is indicated as SE but no SE is available.")
}
tsd.CI <- tsd.CIr
################################################### #
# Prepare logicTransition et pipping_emergence #####
################################################### #
# dans x j'ai les parametres a ajuster
# Il faut que je rajoute les fixes - 16/7/2012
totpar <- c(parameters, fixed.parameters)
if (is.na(totpar["pipping_emergence"])) fixed.parameters <- c(fixed.parameters, pipping_emergence=0)
# logictransition vaut FALSE s'il n'y en a pas
logicTransition <- !(is.na(totpar["transition_P"]) | is.na(totpar["transition_S"]))
##################################### #
# Créé la matrice des random GTRN #####
##################################### #
if (!is.null(SE)) SE <- ifelse(is.infinite(SE), 0, SE)
df_random <- RandomFromHessianOrMCMC(fitted.parameters = parameters ,
fixed.parameters = fixed.parameters ,
se = SE ,
mcmc = resultmcmc ,
Hessian = hessian ,
method = GTRN.CI ,
replicates = replicate.CI ,
chain=1 ,
regularThin=TRUE ,
MinMax=NULL ,
probs=c(0.025, 0.5, 0.975) ,
fn=NULL ,
silent=TRUE ,
ParTofn="par" )$random
##################################### #
# Créé la matrice des random STRN #####
##################################### #
if (is.list(SexualisationTRN)) {
fitted.parameters_SexualisationTRN <- SexualisationTRN$par
fixed.parameters_SexualisationTRN <- SexualisationTRN$fixed.parameters
SexualisationTRN_SE <- SexualisationTRN$SE
SexualisationTRN_Hessian <- SexualisationTRN$hessian
} else {
fitted.parameters_SexualisationTRN <- SexualisationTRN
fixed.parameters_SexualisationTRN <- NULL
SexualisationTRN_SE <- NULL
SexualisationTRN_Hessian <- NULL
}
if (is.null(fixed.parameters_SexualisationTRN) & is.null(fitted.parameters_SexualisationTRN)) {
df_random_SexualisationTRN <- NULL
SexualisationTRN.CI <- "null"
} else {
if (!is.null(SexualisationTRN_SE)) SexualisationTRN_SE <- ifelse(is.infinite(SexualisationTRN_SE), 0, SexualisationTRN_SE)
df_random_SexualisationTRN <- RandomFromHessianOrMCMC(fitted.parameters = fitted.parameters_SexualisationTRN,
fixed.parameters = fixed.parameters_SexualisationTRN,
mcmc = SexualisationTRN.mcmc,
se = SexualisationTRN_SE,
method = SexualisationTRN.CI,
Hessian = SexualisationTRN_Hessian,
replicates = replicate.CI)$random
}
##################################### #
# Créé la matrice des random tsd #####
##################################### #
tsd_fitted.parameters <- tsd$par
tsd_fixed.parameters <- tsd$fixed.parameters
if (is.null(tsd_fitted.parameters) & is.null(tsd_fixed.parameters)) {
df_random_tsd <- NULL
tsd.CI <- "null"
} else {
tsd_SE <- tsd$SE
tsd_SE <- ifelse(is.infinite(tsd_SE), 0, tsd_SE)
df_random_tsd <- RandomFromHessianOrMCMC(fitted.parameters = tsd_fitted.parameters,
fixed.parameters = tsd_fixed.parameters,
mcmc = tsd.mcmc,
se = tsd_SE,
Hessian = tsd$hessian,
method = tsd.CI,
replicates = replicate.CI)$random
}
##################################### #
# concatène les df_random #####
##################################### #
# if (!is.null(df_random_SexualisationTRN)) {
# df_random <- cbind(df_random, df_random_SexualisationTRN)
# }
# if (!is.null(df_random_tsd)) {
# df_random <- cbind(df_random, df_random_tsd)
# }
# df_random, df_random_tsd, df_random_SexualisationTRN
##################################### #
# Je sors les séries que je travaille
##################################### #
NBTs <- temperatures[["IndiceT"]]["NbTS"]
if (is.null(series[1]) | (series[1]=="all")) {
series <- rep(TRUE, NBTs)
}
# Je ne garde que les series que je dois travailler
temperatures_ec <- temperatures[1:NBTs][series]
if (any(is.na(names(temperatures_ec)))) {
stop("The names of 'series' does not exist in series of temperature within 'NestsResult'")
}
# Dans NBTs j'ai maintenant le nombre de series a travailler
NBTs <- length(temperatures_ec)
# dans name j'ai le nom des series que je vais travailler
name <- names(temperatures_ec)
names(name) <- name
# Je crée le df_random en rajoutant une colonne pour la série température
df_random <- as.data.frame(df_random)
df_random <- cbind(Series=rep(name[1], nrow(df_random)), df_random)
df_random_ec <- df_random
if (length(name) > 1) {
for (n in name[-1]) {
df_random_ec[, "Series"] <- rep(n, nrow(df_random_ec))
df_random <- rbind(df_random, df_random_ec)
}
}
df_random_tsd <- as.data.frame(df_random_tsd)
df_random_tsd <- cbind(Series=rep(name[1], nrow(df_random_tsd)), df_random_tsd)
df_random_tsd_ec <- df_random_tsd
if (length(name) > 1) {
for (n in name[-1]) {
df_random_tsd_ec[, "Series"] <- rep(n, nrow(df_random_tsd_ec))
df_random_tsd <- rbind(df_random_tsd, df_random_tsd_ec)
}
}
df_random_SexualisationTRN <- as.data.frame(df_random_SexualisationTRN)
df_random_SexualisationTRN <- cbind(Series=rep(name[1], nrow(df_random_SexualisationTRN)), df_random_SexualisationTRN)
df_random_SexualisationTRN_ec <- df_random_SexualisationTRN
if (length(name) > 1) {
for (n in name[-1]) {
df_random_SexualisationTRN_ec[, "Series"] <- rep(n, nrow(df_random_SexualisationTRN_ec))
df_random_SexualisationTRN <- rbind(df_random_SexualisationTRN, df_random_SexualisationTRN_ec)
}
}
##################################### #
# Les tailles finales attendues
##################################### #
# Je crée un dataframe dans hatchling.metric
if (inherits(hatchling.metric, "numeric")) {
hatchling.metric <- data.frame(Mean=rep(hatchling.metric["Mean"], NBTs),
SD=rep(hatchling.metric["SD"], NBTs), row.names=name)
}
if (is.null(hatchling.metric) | any(is.na(hatchling.metric[name, "Mean"]))) {
# 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])
hatchling.metric <- data.frame(Mean=rep(hatchling.metric["Mean"], NBTs),
SD=rep(hatchling.metric["SD"], NBTs), row.names=name)
} else {
stop("The hatchlings size (hatchling.metric parameter) must be provided for all series.")
}
}
########################################################################################## #
# J'ai donc 1 tableau: df_random avec une ligne par réplicat parameters/series.
# J'ai les tailles finales attendues dans un data.frame hatchling.metric avec deux colonnes #
# J'ai aussi les tailles finales dans SCL.final ? A vérifier #
# Dans name j'ai le nom des séries à travailler #
########################################################################################## #
# The theorem of Cauchy–Lipschitz (also known as the theorem of Picard–Lindelöf)
# gives a set of conditions under which an initial value problem (also named Cauchy problem) has a unique solution.
# Je dois calculer la température moyenne pondérée par le temps d'incubation
# Les températures sont dans temperatures. C'est une liste formattée
# Les séries à travailler sont dans name
# temperature[name]
# Pour chacune je calcule la température moyenne
# Mais cela n'a de sens que si je n'ai pas de metabolic heating ou temperature.heterogeneity et si j'ai bien une taille finale attendue donc !stop.at.hatchling.metric
################################### #
# Lecture des stades pour la TSP ####
################################### #
# Comme maintenant c'est une variable, je dois le rentrer dans la boucle
if (inherits(embryo.stages, "character")) {
if (embryo.stages == "fitted") {
embryo.stages = c("10"=invlogit(logit(0.33)), "11"=invlogit(logit(0.33)),
"12"=invlogit(logit(0.66)), "13"=invlogit(logit(0.66)))
TSP.borders = c(10, 12)
} else {
estages <- TSP.list[[gsub(" ", "_", embryo.stages)]]
if (is.null(estages)) {
stop("The TSP for ", embryo.stages, " does not exist")
} else {
embryo.stages <- estages[, "metric"]
names(embryo.stages) <- estages[, "stages"]
TSP.borders <- c(attributes(estages)$TSP.begin.stages, attributes(estages)$TSP.end.stages)
}
}
}
# dans name j'ai les noms des series
# prepare le data.frame pour les resultats ####
temperatures_ec <- universalmclapply(X=name, mc.cores = papply, progressbar=FALSE, FUN=function (serie_en_cours) {
df <- as.data.frame(temperatures_ec[[serie_en_cours]])
colnames(df) <- c("Time", "TempC", "TempK", "R", "SCL", "IndiceK")
if (df[1, "Time"] != 0) {
df <- rbind(df[1, ], df)
df[1, "Time"] <- 0
}
if (!is.null(fill)) {
newt <- seq(from=0, to=tail(df$Time, n=1L), by=fill)
aretirer <- -which(!is.na(match(newt, df$Time)))
if (!identical(aretirer, integer(0))) newt <- newt[aretirer]
if (!identical(newt, numeric(0))) {
df <- rbind(df,
data.frame(Time=newt, TempC=NA,
TempK=NA, R=NA, SCL=NA, IndiceK=NA))
df <- df[order(df$Time),]
}
}
df <- cbind(df, DeltaT=c(diff(df$Time), 0))
ind <- which(!is.na(df$TempC)) # get positions of non-missing values
if(is.na(df$TempC[1])) { # if it begins with a missing, add the
df$TempC[1] <- df$TempC[ind[1]]
df$TempK[1] <- df$TempK[ind[1]]
df$IndiceK[1] <- df$IndiceK[ind[1]]
ind <- c(1,ind)
}# first position to the indices
df$TempC <- rep(df$TempC[ind], times = diff( # repeat the values at these indices
c(ind, length(df$TempC) + 1) ))
df$TempK <- rep(df$TempK[ind], times = diff( # repeat the values at these indices
c(ind, length(df$TempK) + 1) ))
df$IndiceK <- rep(df$IndiceK[ind], times = diff( # repeat the values at these indices
c(ind, length(df$IndiceK) + 1) ))
df[, "SCL"] <- NA
df[1, "SCL"] <- M0
df[, "IndiceK"] <- 1
return(df)
}, clusterExport=list(varlist=c("temperatures_ec", "fill", "M0"),
envir = environment()),
clusterEvalQ=list(expr=expression(library(embryogrowth)))
)
# Je les classe
temperatures_ec <- temperatures_ec[name]
################################################################### #
# Je prépare le calcul avec tous les nids et tous les réplicats ####
################################################################### #
AnalyseTraces <- universalmclapply(progressbar=progressbar, mc.cores = papply, X=1:nrow(df_random),
FUN=function (ligne_df_random) {
# ligne_df_random <- 1
# print(ligne_df_random)
serie_en_cours <- df_random[ligne_df_random, "Series"]
df <- temperatures_ec[[serie_en_cours]]
setpar <- c(unlist(df_random[ligne_df_random, -1, drop=TRUE]), unlist(df_random_tsd[ligne_df_random, -1, drop=TRUE]))
setpar_STRN <- unlist(df_random_SexualisationTRN[ligne_df_random, -1, drop=TRUE])
# Si je n'ai pas setpar_STRN, je le mets en uniforme; dans dbeta_mu, j'ai le logit de mu (0=0.5)
if (is.null(setpar_STRN)) setpar_STRN <- c(dbeta_mu=0, dbeta_v=1/12)
if (all(substr(names(setpar_STRN), 1, 6) != "dbeta_")) {
setpar_STRN["dbeta_mu"] <- 0
setpar_STRN["dbeta_v"] <- 1/12
}
transition_P <- setpar["transition_P"]
transition_S <- setpar["transition_S"]
# "df_random_tsd", "df_random_SexualisationTRN"
if (!is.na(setpar["K"])) {
Kval <- unname(setpar["K"])
} else {
Kval <- NULL
}
if (length(temperature.heterogeneity)==1) {
tht <- 0
if (temperature.heterogeneity !=0) tht <- rnorm(1, 0, temperature.heterogeneity)
} else {
tht <- r2norm(1, 0, sd_low = temperature.heterogeneity[1], sd_high = temperature.heterogeneity[2])
}
df[, "TempC"] <- df[, "TempC"] + tht
df[, "TempK"] <- df[, "TempK"] + tht
meanSCL <- as.numeric(hatchling.metric[serie_en_cours, "Mean"])
sdSCL <- as.numeric(hatchling.metric[serie_en_cours, "SD"])
if (!is.na(setpar["rK"])) {
Kval <- unname(setpar["rK"]*meanSCL)
}
tmin <- df[, "Time"]
for (j in 2:nrow(df)) if (is.na(df[j, "TempK"])) df[j, "TempK"] <- df[j-1, "TempK"]
dT <- df[, "DeltaT"]
y0 <- ypre <- df[1, "SCL"]
if (!stop.at.hatchling.metric) {
if (!logicTransition) {
# if ((metabolic.heating == 0) &
# ((length(temperature.heterogeneity) == 1) &
# (temperature.heterogeneity[1] ==0)) & (out=="likelihood")) {
# indice.boucle <- nrow(df) - 1
# } else {
# indice.boucle <- 1
# }
indice.boucle <- 1
# Je prends de 1 jusqu'à l'avant-dernier
for (i in indice.boucle:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
a_inst <- SSM(CTE_inst, setpar, verbose=verbose)[[1]]
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
}
df[nrow(df), "TempK"] <- df[nrow(df), "TempK"] + metabolic.heating*(ypre/meanSCL)
} else {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
transition <- 1/(1+exp(transition_S*(ypre-transition_P)))
rT_inst < SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
}
df[nrow(df), "TempK"] <- df[nrow(df), "TempK"] + metabolic.heating*(ypre/meanSCL)
}
# je suis en stop.at.hatchling.metric
} else {
if (!logicTransition) {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
a_inst <- SSM(CTE_inst, setpar, verbose=verbose)[[1]]
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
if (ypre>meanSCL) {
break
}
}
} else {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
transition <- 1/(1+exp(transition_S*(ypre-transition_P)))
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
if (ypre>meanSCL) {
break
}
}
# fin du else logictransition
}
# fin du else !stop.at.hatchling.metric
}
df[, "TempC"] <- df[, "TempK"] - 273.15
if (stop.at.hatchling.metric) {
size.final <- ifelse(ypre < meanSCL, NA, meanSCL)
} else {
size.final <- ifelse(metric.end.incubation == "observed", ypre, meanSCL)
}
if (!is.na(size.final)) {
if (!is.null(embryo.stages)) {
if (!is.na(setpar_STRN["BeginTSP"])) {
xxbt <- invlogit(setpar_STRN["BeginTSP"])
if (xxbt == 0) xxbt <- zero
if (xxbt > 0.95) xxbt <- 0.95
embryo.stages[as.character(TSP.borders[1])] <- xxbt
embryo.stages[as.character(TSP.borders[1]+1)] <- xxbt
}
if (!is.na(setpar_STRN["EndTSP"])) {
xxbt <- invlogit(setpar_STRN["EndTSP"])
if (xxbt == 0) xxbt <- zero
if (xxbt >= 0.99) xxbt <- 0.99
embryo.stages[as.character(TSP.borders[2])] <- xxbt
embryo.stages[as.character(TSP.borders[2]+1)] <- xxbt
}
if (!is.na(setpar_STRN["LengthTSP"])) {
xxbt <- invlogit(setpar_STRN["BeginTSP"]+abs(setpar_STRN["LengthTSP"]))
if (xxbt == 0) xxbt <- zero
if (xxbt >= 0.99) xxbt <- 0.99
embryo.stages[as.character(TSP.borders[2])] <- xxbt
embryo.stages[as.character(TSP.borders[2]+1)] <- xxbt
}
# Je ne comprends pas le mean ici - 2021-05-12
size.begin.TSP <- mean(embryo.stages[as.character(TSP.borders[1])])
if (TSP.begin != 0) {
size.begin.TSP <- size.begin.TSP+(embryo.stages[as.character(TSP.borders[1]+1)]-size.begin.TSP)*TSP.begin
}
size.end.TSP <- mean(embryo.stages[as.character(TSP.borders[2])])
if (TSP.end != 0) {
size.end.TSP <- size.end.TSP+(embryo.stages[as.character(TSP.borders[2]+1)]-size.end.TSP)*TSP.end
}
size.begin.TSP_ec <- size.begin.TSP*size.final
size.end.TSP_ec <- size.end.TSP*size.final
} else {
size.end.TSP <- NA
size.begin.TSP <- NA
size.begin.TSP_ec <- NA
size.end.TSP_ec <- NA
}
if (out == "likelihood") {
if (!is.na(sdSCL)) {
likelihood <- (-dnorm(ypre, mean=meanSCL, sd=sdSCL, log=TRUE))
} else {
if (is.na(parameters["SD"])) {
likelihood <- ((ypre - meanSCL)^2)
} else {
# 28/3/2022
likelihood <- (-dnorm(ypre, mean=meanSCL, sd=parameters["SD"], log=TRUE))
}
}
return(list(likelihood=likelihood))
}
df_ec <- df
##### Estimate limit begin TSP #####
if ((size.final < size.begin.TSP_ec) | is.na(size.begin.TSP_ec) | (max(df_ec[,"SCL"], na.rm = TRUE) < size.begin.TSP_ec)) {
time.begin.TSP <- NA
indice.begin.tsp <- NA
} else {
if (size.begin.TSP_ec < df_ec[1,"SCL"]) {
size.begin.TSP_ec <- df_ec[1,"SCL"]+(df_ec[2,"SCL"] - df_ec[1,"SCL"])/10
}
indice.begin.tsp <- which(df_ec[,"SCL"]>(size.begin.TSP_ec))[1]-1
df_ec <- rbind(df_ec[1:indice.begin.tsp, ],
data.frame(Time=c(df_ec[indice.begin.tsp,"Time"]+df_ec[indice.begin.tsp,"DeltaT"]/3,
df_ec[indice.begin.tsp,"Time"]+df_ec[indice.begin.tsp,"DeltaT"]*2/3),
TempC=c(df_ec[indice.begin.tsp,"TempC"],
df_ec[indice.begin.tsp,"TempC"]),
TempK=c(df_ec[indice.begin.tsp,"TempK"],
df_ec[indice.begin.tsp,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.begin.tsp+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
df_ec[indice.begin.tsp, "DeltaT"] <- df_ec[indice.begin.tsp+1,"Time"]-df_ec[indice.begin.tsp,"Time"]
df_ec[indice.begin.tsp+1, "DeltaT"] <- df_ec[indice.begin.tsp+2,"Time"]-df_ec[indice.begin.tsp+1,"Time"]
df_ec[indice.begin.tsp+2, "DeltaT"] <- df_ec[indice.begin.tsp+3,"Time"]-df_ec[indice.begin.tsp+2,"Time"]
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
CTE_inst <- df_ec[indice.begin.tsp, "TempK"]
timesunique <- c(df_ec[indice.begin.tsp, "Time"], df_ec[indice.begin.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.begin.tsp], setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.tsp, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+1, "R"] <- a_inst
df_ec[indice.begin.tsp+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.begin.tsp+1, "Time"], df_ec[indice.begin.tsp+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
CTE_inst <- df_ec[indice.begin.tsp+1, "TempK"]
# rT <- SSM(CTE[indice.begin.tsp+1], setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre,
parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+2, "R"] <- a_inst
df_ec[indice.begin.tsp+2, "SCL"] <- ypre
# J'ai encadré le début de la TSP
t <- df_ec[indice.begin.tsp:(indice.begin.tsp+3), "Time"]
SCL <- df_ec[indice.begin.tsp:(indice.begin.tsp+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.begin.tsp <- predict(l, newdata = data.frame(SCL=size.begin.TSP_ec))
time.begin.tsp <- round(time.begin.tsp, 2)
if (time.begin.tsp <= 0) {
time.begin.tsp <- t[1]+zero
}
indice.begin.tsp <- which(df_ec[,"Time"]>(time.begin.tsp))[1]-1
df_ec <- rbind(df_ec[1:indice.begin.tsp, ],
data.frame(Time=time.begin.tsp,
TempC=df_ec[indice.begin.tsp,"TempC"],
TempK=df_ec[indice.begin.tsp,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.begin.tsp+1):nrow(df_ec), ])
# df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
df_ec[indice.begin.tsp, "DeltaT"] <- df_ec[indice.begin.tsp+1, "Time"]-df_ec[indice.begin.tsp, "Time"]
df_ec[indice.begin.tsp+1, "DeltaT"] <- df_ec[indice.begin.tsp+2, "Time"]-df_ec[indice.begin.tsp+1, "Time"]
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
timesunique <- c(df_ec[indice.begin.tsp, "Time"], df_ec[indice.begin.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.begin.tsp], setpar)
CTE_inst <- df_ec[indice.begin.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.tsp, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+1, "R"] <- a_inst
df_ec[indice.begin.tsp+1, "SCL"] <- ypre
indice.begin.tsp <- indice.begin.tsp+1
rownames(df_ec)[indice.begin.tsp] <- "Begin TSP"
time.begin.TSP <- df_ec[indice.begin.tsp,"Time"]
}
###### Estimate limit end TSP ####
if ((size.final < size.end.TSP_ec) | is.na(size.end.TSP_ec) | (max(df_ec[,"SCL"], na.rm = TRUE) < size.end.TSP_ec)) {
time.end.TSP <- NA
indice.end.tsp <- NA
} else {
if (size.end.TSP_ec < df_ec[1,"SCL"]) {
size.end.TSP_ec <- df_ec[1,"SCL"]+(df_ec[2,"SCL"] - df_ec[1,"SCL"])/10
}
ind_nonNA <- max(which(!is.na(df_ec[, "SCL"])))
if (size.end.TSP_ec >= df_ec[ind_nonNA,"SCL"]) {
size.end.TSP_ec <- df_ec[ind_nonNA-1,"SCL"] + (df_ec[ind_nonNA,"SCL"]-df_ec[ind_nonNA-1,"SCL"])/10
}
if (size.end.TSP_ec <= df_ec[indice.begin.tsp, "SCL"]) {
size.end.TSP_ec <- df_ec[indice.begin.tsp, "SCL"]+(df_ec[indice.begin.tsp+1, "SCL"]-df_ec[indice.begin.tsp, "SCL"])/10
}
indice.end.tsp <- which(df_ec[,"SCL"]>(size.end.TSP_ec))[1]-1
if (is.na(indice.end.tsp)) indice.end.tsp <- nrow(df_ec)
# if (is.na(indice.end.tsp)) {
# print("AAAA 1154")
# print(paste0("size.end.TSP_ec", as.character(size.end.TSP_ec)))
# save(df_ec, file="df_ec.Rdata")
# stop()
# }
df_ec <- rbind(df_ec[1:indice.end.tsp, ],
data.frame(Time=c(df_ec[indice.end.tsp,"Time"]+df_ec[indice.end.tsp,"DeltaT"]/3,
df_ec[indice.end.tsp,"Time"]+df_ec[indice.end.tsp,"DeltaT"]*2/3),
TempC=c(df_ec[indice.end.tsp,"TempC"],
df_ec[indice.end.tsp,"TempC"]),
TempK=c(df_ec[indice.end.tsp,"TempK"],
df_ec[indice.end.tsp,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.end.tsp+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# sum(dt[1:(indice.end.tsp+1)]*TK[1:(indice.end.tsp+1)])/sum(dt[1:(indice.end.tsp+1)])
# CTE <- cumsum(dt*TK)/cumsum(dt)
df_ec[indice.end.tsp, "DeltaT"] <- df_ec[indice.end.tsp+1,"Time"]-df_ec[indice.end.tsp,"Time"]
df_ec[indice.end.tsp+1, "DeltaT"] <- df_ec[indice.end.tsp+2,"Time"]-df_ec[indice.end.tsp+1,"Time"]
df_ec[indice.end.tsp+2, "DeltaT"] <- df_ec[indice.end.tsp+3,"Time"]-df_ec[indice.end.tsp+2,"Time"]
timesunique <- c(df_ec[indice.end.tsp, "Time"], df_ec[indice.end.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.end.tsp], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.end.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.tsp, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.tsp+1, "R"] <- a_inst
df_ec[indice.end.tsp+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.end.tsp+1, "Time"], df_ec[indice.end.tsp+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
# rT <- SSM(CTE[indice.end.tsp+1], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.end.tsp+1, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre,
parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.tsp+2, "R"] <- a_inst
df_ec[indice.end.tsp+2, "SCL"] <- ypre
t <- df_ec[indice.end.tsp:(indice.end.tsp+3), "Time"]
SCL <- df_ec[indice.end.tsp:(indice.end.tsp+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.end.tsp <- predict(l, newdata = data.frame(SCL=size.end.TSP_ec))
time.end.tsp <- round(time.end.tsp, 2)
# plot(t, SCL); points(time.begin.tsp, size.begin.TSP_ec, col="red"); segments(x0=0, y0=size.begin.TSP_ec, x1=100000, y1=size.begin.TSP_ec, col="green")
indice.end.tsp <- which(df_ec[,"Time"]>(time.end.tsp))[1]-1
# if (is.na(indice.end.tsp)) {
# indice.end.tsp <- nrow(df_ec)
#
# df_ec <- rbind(df_ec, df_ec[nrow(df_ec), ])
# } else {
#
df_ec <- rbind(df_ec[1:indice.end.tsp, ],
data.frame(Time=time.end.tsp,
TempC=df_ec[indice.end.tsp,"TempC"],
TempK=df_ec[indice.end.tsp,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.end.tsp+1):nrow(df_ec), ])
df_ec[indice.end.tsp, "DeltaT"] <- df_ec[indice.end.tsp+1, "Time"]-df_ec[indice.end.tsp, "Time"]
df_ec[indice.end.tsp+1, "DeltaT"] <- df_ec[indice.end.tsp+2, "Time"]-df_ec[indice.end.tsp+1, "Time"]
# }
timesunique <- c(df_ec[indice.end.tsp, "Time"], df_ec[indice.end.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.tsp, "SCL"] - transition_P)))
}
CTE_inst <- df_ec[indice.end.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.tsp, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
indice.end.tsp <- indice.end.tsp+1
df_ec[indice.end.tsp, "R"] <- a_inst
df_ec[indice.end.tsp, "SCL"] <- ypre
rownames(df_ec)[indice.end.tsp] <- "End TSP"
time.end.TSP <- df_ec[indice.end.tsp,"Time"]
}
#### End of incubation for stop.at.hatchling.metric ##########
# attSAT vaut NA si je ne suis pas en stop.at.hatchling.metric
# attSAT vaut TRUE si l'incubation a été jusquà meanSCL
# attSAT vaut FALSE si l'incubation n'a pas été jusqu'à meanSCL
if (stop.at.hatchling.metric) {
if (max(df_ec[,"SCL"], na.rm=TRUE) > meanSCL) {
attSAT <- TRUE
indice.fin.incubation <- which(df_ec[,"SCL"]> meanSCL)[1]-1
df_ec <- rbind(df_ec[1:indice.fin.incubation, ],
data.frame(Time=c(df_ec[indice.fin.incubation,"Time"]+df_ec[indice.fin.incubation,"DeltaT"]/3,
df_ec[indice.fin.incubation,"Time"]+df_ec[indice.fin.incubation,"DeltaT"]*2/3),
TempC=c(df_ec[indice.fin.incubation,"TempC"],
df_ec[indice.fin.incubation,"TempC"]),
TempK=c(df_ec[indice.fin.incubation,"TempK"],
df_ec[indice.fin.incubation,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.fin.incubation+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
#
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
#
# CTE <- cumsum(dt*TK)/cumsum(dt)
df_ec[indice.fin.incubation, "DeltaT"] <- df_ec[indice.fin.incubation+1,"Time"]-df_ec[indice.fin.incubation,"Time"]
df_ec[indice.fin.incubation+1, "DeltaT"] <- df_ec[indice.fin.incubation+2,"Time"]-df_ec[indice.fin.incubation+1,"Time"]
df_ec[indice.fin.incubation+2, "DeltaT"] <- df_ec[indice.fin.incubation+3,"Time"]-df_ec[indice.fin.incubation+2,"Time"]
timesunique <- c(df_ec[indice.fin.incubation, "Time"], df_ec[indice.fin.incubation+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.fin.incubation, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.fin.incubation, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.fin.incubation, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+1, "R"] <- a_inst
df_ec[indice.fin.incubation+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.fin.incubation+1, "Time"], df_ec[indice.fin.incubation+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation+1], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation+1, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+2, "R"] <- a_inst
df_ec[indice.fin.incubation+2, "SCL"] <- ypre
t <- df_ec[indice.fin.incubation:(indice.fin.incubation+3), "Time"]
SCL <- df_ec[indice.fin.incubation:(indice.fin.incubation+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.fin.incubation <- predict(l, newdata = data.frame(SCL=meanSCL))
time.fin.incubation <- round(time.fin.incubation, 2)
# plot(t, SCL); points(time.begin.tsp, size.begin.TSP_ec, col="red"); segments(x0=0, y0=size.begin.TSP_ec, x1=100000, y1=size.begin.TSP_ec, col="green")
indice.fin.incubation <- which(df_ec[,"Time"]>(time.fin.incubation))[1]-1
df_ec <- rbind(df_ec[1:indice.fin.incubation, ],
data.frame(Time=time.fin.incubation,
TempC=df_ec[indice.fin.incubation,"TempC"],
TempK=df_ec[indice.fin.incubation,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.fin.incubation+1):nrow(df_ec), ])
df_ec[indice.fin.incubation, "DeltaT"] <- df_ec[indice.fin.incubation+1, "Time"]-df_ec[indice.fin.incubation, "Time"]
df_ec[indice.fin.incubation+1, "DeltaT"] <- df_ec[indice.fin.incubation+2, "Time"]-df_ec[indice.fin.incubation+1, "Time"]
# df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
#
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
timesunique <- c(df_ec[indice.fin.incubation, "Time"], df_ec[indice.fin.incubation+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.fin.incubation, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.fin.incubation, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.fin.incubation, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+1, "R"] <- a_inst
df_ec[indice.fin.incubation+1, "SCL"] <- ypre
indice.fin.incubation <- indice.fin.incubation+1
rownames(df_ec)[indice.fin.incubation] <- "End Incubation"
time.end.incubation <- df_ec[indice.fin.incubation,"Time"]
# Maintenant je tronque df de faon ne garder que jusqu'a time.end.incubation
# garde
df_ec <- df_ec[1:indice.fin.incubation,]
# Je garde la fin de l'incubation
df_ec[indice.fin.incubation, "IndiceK"] <- 1
df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
} else {
attSAT <- FALSE
indice.fin.incubation <- max(which(!is.na(df_ec[, "SCL"])))
rownames(df_ec)[indice.fin.incubation] <- "End Incubation"
time.end.incubation <- df_ec[indice.fin.incubation,"Time"]
df_ec[indice.fin.incubation, "IndiceK"] <- 1
}
} else {attSAT <- NA}
################################## #
# indice.begin.middlethird ####
################################## #
# On peut la calculer seulement si attSAT == TRUE ou NA
if (attSAT | is.na(attSAT)) {
time.begin.middlethird <- unname((tail(df_ec$Time, n=1L)-setpar["pipping_emergence"])/3)
indice.begin.middlethird <- which(df_ec[,"Time"]>time.begin.middlethird)[1]-1
df_ec <- rbind(df_ec[1:indice.begin.middlethird, ],
data.frame(Time=time.begin.middlethird,
TempC=df_ec[indice.begin.middlethird, "TempK"]-273.15,
TempK=df_ec[indice.begin.middlethird, "TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT= NA),
df_ec[(indice.begin.middlethird+1):nrow(df_ec), ])
df_ec[indice.begin.middlethird,"DeltaT"] <- df_ec[indice.begin.middlethird + 1,"Time"] - df_ec[indice.begin.middlethird,"Time"]
df_ec[indice.begin.middlethird + 1,"DeltaT"] <- df_ec[indice.begin.middlethird + 2,"Time"] - df_ec[indice.begin.middlethird + 1,"Time"]
# sum(df_ec[1:10,"DeltaT"]); df_ec[10+1,"Time"]
# sum(df_ec[1:indice.begin.middlethird,"DeltaT"]); df_ec[indice.begin.middlethird+1,"Time"]
# timei1 <- df_ec[indice.begin.middlethird + 1,"Time"]
# dt <- df_ec[1:indice.begin.middlethird, "DeltaT"]
# TK <- df_ec[1:indice.begin.middlethird, "TempK"]
# CTE <- sum(dt*TK)/sum(dt)
CTE_inst <- df_ec[indice.begin.middlethird, "TempK"]
timesunique <- c(df_ec[indice.begin.middlethird,"Time"], df_ec[indice.begin.middlethird + 1,"Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.middlethird, "SCL"] - transition_P)))
}
# rT <- SSM(CTE, setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.middlethird, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.middlethird, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.middlethird+1, "R"] <- a_inst
df_ec[indice.begin.middlethird+1, "SCL"] <- ypre
rownames(df_ec)[indice.begin.middlethird+1] <- "Begin Middle-third"
size.begin.middlethird <- ypre
} else {
size.begin.middlethird <- NA
time.begin.middlethird <- NA
}
################################ #
# indice.end.middlethird ####
################################ #
# On peut la calculer seulement si attSAT == TRUE ou NA
if (attSAT | is.na(attSAT)) {
time.end.middlethird <- unname((tail(df_ec$Time, n=1L)-setpar["pipping_emergence"])*2/3)
indice.end.middlethird <- which(df_ec[,"Time"]>time.end.middlethird)[1]-1
# Je calcul la taille
df_ec <- rbind(df_ec[1:indice.end.middlethird, ],
data.frame(Time=time.end.middlethird,
TempC=df_ec[indice.end.middlethird, "TempK"]-273.15,
TempK=df_ec[indice.end.middlethird, "TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.end.middlethird+1):nrow(df_ec), ])
df_ec[indice.end.middlethird,"DeltaT"] <- df_ec[indice.end.middlethird + 1,"Time"] - df_ec[indice.end.middlethird,"Time"]
df_ec[indice.end.middlethird + 1,"DeltaT"] <- df_ec[indice.end.middlethird + 2,"Time"] - df_ec[indice.end.middlethird + 1,"Time"]
# timei1 <- df_ec[indice.end.middlethird+1,"Time"]
# dt <- df_ec[1:indice.end.middlethird, "DeltaT"]
# TK <- df_ec[1:indice.end.middlethird, "TempK"]
# CTE <- sum(dt*TK)/sum(dt)
CTE_inst <- df_ec[indice.end.middlethird, "TempK"]
timesunique <- c(df_ec[indice.end.middlethird,"Time"], df_ec[indice.end.middlethird+1,"Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.middlethird, "SCL"] - transition_P)))
}
# rT <- SSM(CTE, setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.middlethird, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.middlethird, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.middlethird+1, "R"] <- a_inst
df_ec[indice.end.middlethird+1, "SCL"] <- ypre
rownames(df_ec)[indice.end.middlethird+1] <- "End Middle-third"
size.end.middlethird <- ypre
} else {
size.end.middlethird <- NA
time.end.middlethird <- NA
}
df <- df_ec
# df[, "DeltaT"] <- c(diff(df[, "Time"]), 0)
indice.begin.middlethird <- which(rownames(df) == "Begin Middle-third")
if (identical(indice.begin.middlethird, integer(0))) indice.begin.middlethird <- NA
indice.end.middlethird <- which(rownames(df) == "End Middle-third")
if (identical(indice.end.middlethird, integer(0))) indice.end.middlethird <- NA
indice.fin.incubation <- nrow(df)
indice.begin.tsp <- which(rownames(df) == "Begin TSP")
indice.end.tsp <- which(rownames(df) == "End TSP")
if (identical(integer(0), indice.end.tsp)) indice.end.tsp <- NA
if (identical(integer(0), indice.begin.tsp)) indice.begin.tsp <- NA
## Je viens de calculer la croissance de l'embryon
if (out == "summary") {
# Maintenant je calcule les summary
TC <- df[1:(nrow(df)-1), "TempC"]
dt <- df[1:(nrow(df)-1), "DeltaT"]
dSCL <- diff(df[, "SCL"])
rT <- df[1:(nrow(df)-1), "R"]
Temperature.max <- max(TC)
TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
Incubation.length <- df[nrow(df), "Time"]
if (!is.na(indice.begin.tsp) & !is.na(indice.end.tsp)) {
TC <- df[indice.begin.tsp:(indice.end.tsp-1), "TempC"]
dt <- df[indice.begin.tsp:(indice.end.tsp-1), "DeltaT"]
dSCL <- diff(df[indice.begin.tsp:(indice.end.tsp), "SCL"])
rT <- df[indice.begin.tsp:(indice.end.tsp-1), "R"]
TSP.TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
TSP.GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
TSP.TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
TSP.length <- unname(time.end.tsp - time.begin.tsp)
TSP.begin <- unname(time.begin.tsp)
TSP.end <- unname(time.end.tsp)
if (!is.null(setpar_STRN)) {
if (indice.begin.tsp == indice.end.tsp-1) {
indice.end.tsp_ec <- indice.end.tsp+1
} else {
indice.end.tsp_ec <- indice.end.tsp
}
TK <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "TempK"]
TC <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "TempC"]
dt <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "DeltaT"]
dSCL <- diff(df[indice.begin.tsp:(indice.end.tsp_ec), "SCL"])
rT <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "R"]
# print("setpar_STRN")
# print(setpar_STRN)
if (any(is.finite(suppressWarnings(as.numeric(names(setpar_STRN)))))) {
sparSTRN <- setpar_STRN[!is.na(suppressWarnings(as.numeric(names(setpar_STRN))))]
# print("sparSTRN")
# print(sparSTRN)
STRN_model <- SSM(TK, sparSTRN, verbose=verbose)[[1]]
STRN_model <- ifelse(STRN_model <= 0, 0.001, STRN_model)
# print("SSM")
# print(STRN_model)
} else {
if (any(names(setpar_STRN) == "Rho25") | any(names(setpar_STRN) == "Peak")) {
STRN_model <- SSM(TK, setpar_STRN, verbose=verbose)[[1]]
} else {
#print("setpar_STRN")
# print(setpar_STRN)
STRN_model <- rep(1, length(TK))
# print("SSM")
# print(STRN_model)
}
}
SCL <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "SCL"]
SCL <- (SCL-SCL[1])/(SCL[length(SCL)]-SCL[1])
SCL[SCL==0] <- min(zero, min(SCL[SCL != 0])/10)
SCL[SCL==1] <- max(1-zero, max(SCL[SCL != 1]) + (1-max(SCL[SCL != 1]))/10)
if (!is.na(setpar_STRN["dbeta_mu"])) {
mu <- invlogit(setpar_STRN["dbeta_mu"])
v <- abs(setpar_STRN["dbeta_v"])
shape1 <- mu * (((mu * (1 - mu))/v) - 1)
shape2 <- shape1 * (1 - mu)/mu
model_beta <- dbeta(SCL, shape1=shape1, shape2=shape2)
} else {
shape1 <- abs(setpar_STRN["dbeta_shape1"])
shape2 <- abs(setpar_STRN["dbeta_shape2"])
model_beta <- dbeta(SCL, shape1=shape1, shape2=shape2)
}
TSP.TimeWeighted.STRNWeighted.temperature <- sum(TC * dt * STRN_model * model_beta)/sum(dt * STRN_model * model_beta)
TSP.GrowthWeighted.STRNWeighted.temperature <- sum(TC * dSCL * STRN_model * model_beta)/sum(dSCL * STRN_model * model_beta)
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- sum(TC * dt * rT * STRN_model * model_beta)/sum(dt * rT * STRN_model * model_beta)
} else {
TSP.TimeWeighted.STRNWeighted.temperature <- NA
TSP.GrowthWeighted.STRNWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- NA
}
} else {
TSP.PM.GrowthWeighted <- NA
TSP.PM.TimeWeighted.GrowthRateWeighted <- NA
TSP.PM.TimeWeighted <- NA
TSP.TimeWeighted.temperature <- NA
TSP.GrowthWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.temperature <- NA
TSP.TimeWeighted.STRNWeighted.temperature <- NA
TSP.GrowthWeighted.STRNWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- NA
TSP.length <- NA
}
if (!is.na(indice.begin.middlethird) & !is.na(indice.end.middlethird)) {
TC <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "TempC"])
dt <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "DeltaT"])
dSCL <- unname(diff(df[indice.begin.middlethird:(indice.end.middlethird), "SCL"]))
rT <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "R"])
MiddleThird.TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
MiddleThird.GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
MiddleThird.length <- unname(time.end.middlethird - time.begin.middlethird)
MiddleThird.begin <- unname(time.begin.middlethird)
MiddleThird.end <- unname(time.end.middlethird)
} else {
MiddleThird.length <- NA
MiddleThird.begin <- NA
MiddleThird.end <- NA
MiddleThird.TimeWeighted.temperature <- NA
MiddleThird.GrowthWeighted.temperature <- NA
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature <- NA
}
if (!is.null(tsd)) {
tsdec <- tsd
tsdec$par <- setpar[names(tsd$par)]
tsdec$fixed.parameters <- setpar[names(tsd$fixed.parameters)]
PM <- predict(tsdec, temperatures = TC, replicate.CI = NULL)
TSP.PM.GrowthWeighted <- sum(PM["50%", ]*dSCL)/sum(dSCL)
TSP.PM.TimeWeighted.GrowthRateWeighted <- sum(PM["50%", ]* rT * dt)/sum(rT * dt)
TSP.PM.TimeWeighted <- sum(PM["50%", ]* dt)/sum(dt)
TSP.TimeWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.GrowthWeighted.sexratio <- predict(tsd, temperatures = TSP.GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.TimeWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.GrowthWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TimeWeighted.sexratio <- predict(tsd, temperatures = TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
GrowthWeighted.sexratio <- predict(tsd, temperatures = GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
if (!is.na(TSP.TimeWeighted.STRNWeighted.temperature)) {
TSP.TimeWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.GrowthWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.GrowthWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
} else {
TSP.TimeWeighted.STRNWeighted.sexratio <- NA
TSP.GrowthWeighted.STRNWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- NA
}
} else {
TSP.PM.GrowthWeighted <- NA
TSP.PM.TimeWeighted.GrowthRateWeighted <- NA
TSP.PM.TimeWeighted <- NA
MiddleThird.TimeWeighted.sexratio <- NA
MiddleThird.GrowthWeighted.sexratio <- NA
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio <- NA
TimeWeighted.sexratio <- NA
GrowthWeighted.sexratio <- NA
TimeWeighted.GrowthRateWeighted.sexratio <- NA
TSP.TimeWeighted.sexratio <- NA
TSP.GrowthWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.sexratio <- NA
TSP.TimeWeighted.STRNWeighted.sexratio <- NA
TSP.GrowthWeighted.STRNWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- NA
}
summary <- data.frame(series=serie_en_cours,
Temperature.max=Temperature.max,
TimeWeighted.temperature=TimeWeighted.temperature,
GrowthWeighted.temperature=GrowthWeighted.temperature,
TimeWeighted.GrowthRateWeighted.temperature=TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.temperature=TSP.TimeWeighted.temperature,
TSP.GrowthWeighted.temperature=TSP.GrowthWeighted.temperature,
TSP.TimeWeighted.GrowthRateWeighted.temperature=TSP.TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.STRNWeighted.temperature=TSP.TimeWeighted.STRNWeighted.temperature,
TSP.GrowthWeighted.STRNWeighted.temperature=TSP.GrowthWeighted.STRNWeighted.temperature,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature=TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature,
TSP.length=TSP.length,
TSP.begin=TSP.begin,
TSP.end=TSP.end,
TSP.PM.GrowthWeighted=TSP.PM.GrowthWeighted,
TSP.PM.TimeWeighted.GrowthRateWeighted=TSP.PM.TimeWeighted.GrowthRateWeighted,
TSP.PM.TimeWeighted=TSP.PM.TimeWeighted,
Incubation.length=Incubation.length,
MiddleThird.length=MiddleThird.length,
MiddleThird.begin=MiddleThird.begin,
MiddleThird.end=MiddleThird.end,
MiddleThird.TimeWeighted.temperature=MiddleThird.TimeWeighted.temperature,
MiddleThird.GrowthWeighted.temperature=MiddleThird.GrowthWeighted.temperature,
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature=MiddleThird.TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.sexratio = TSP.TimeWeighted.sexratio,
TSP.GrowthWeighted.sexratio = TSP.GrowthWeighted.sexratio,
TSP.TimeWeighted.GrowthRateWeighted.sexratio = TSP.TimeWeighted.GrowthRateWeighted.sexratio,
TSP.TimeWeighted.STRNWeighted.sexratio = TSP.TimeWeighted.STRNWeighted.sexratio,
TSP.GrowthWeighted.STRNWeighted.sexratio = TSP.GrowthWeighted.STRNWeighted.sexratio,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio = TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio,
MiddleThird.TimeWeighted.sexratio = MiddleThird.TimeWeighted.sexratio,
MiddleThird.GrowthWeighted.sexratio = MiddleThird.GrowthWeighted.sexratio,
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio = MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio,
TimeWeighted.sexratio = TimeWeighted.sexratio,
GrowthWeighted.sexratio = GrowthWeighted.sexratio,
TimeWeighted.GrowthRateWeighted.sexratio = TimeWeighted.GrowthRateWeighted.sexratio)
# je ne garde que les valeurs d'origine
df <- df[df[, "IndiceK"] == 1, ]
indices.df <- data.frame(series=serie_en_cours,
metric.begin.tsp = unname(size.begin.TSP_ec),
metric.end.tsp = unname(size.end.TSP_ec),
hatchling.metric.mean = unname(meanSCL),
hatchling.metric.sd = unname(sdSCL),
time.end.tsp = unname(time.end.TSP),
time.begin.tsp = unname(time.begin.TSP),
metric.end.incubation = unname(df[nrow(df), "SCL"]),
time.begin.middlethird = unname(time.begin.middlethird),
time.end.middlethird = unname(time.end.middlethird),
stop.at.hatchling.metric = unname(attSAT))
} else {
# Sinon, je stocke aussi les indices
# rownames(df) <- 1:nrow(df)
summary <- NULL
indices.df <- data.frame(series=serie_en_cours,
metric.begin.tsp = unname(size.begin.TSP_ec),
metric.end.tsp = unname(size.end.TSP_ec),
hatchling.metric.mean = unname(meanSCL),
hatchling.metric.sd = unname(sdSCL),
time.end.tsp = unname(time.end.TSP),
time.begin.tsp = unname(time.begin.TSP),
metric.end.incubation = unname(df[nrow(df), "SCL"]),
time.begin.middlethird = unname(time.begin.middlethird),
time.end.middlethird = unname(time.end.middlethird),
stop.at.hatchling.metric = unname(attSAT),
indice.begin.middlethird = indice.begin.middlethird,
indice.end.middlethird = indice.end.middlethird,
indice.fin.incubation = indice.fin.incubation,
indice.begin.tsp = indice.begin.tsp,
indice.end.tsp = indice.end.tsp)
}
df[, "DeltaT"] <- c(diff(df[, "Time"]), 0)
} else {
# Je suis en stop.at.hatchling size mais trop court
indices.df <- NULL
summary <- data.frame(series=serie_en_cours,
Temperature.max=max(df[, "TempC"], na.rm = TRUE),
TimeWeighted.temperature=sum(df[, "TempC"]*df[, "DeltaT"])/sum(df[, "DeltaT"]),
GrowthWeighted.temperature=NA,
TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.temperature=NA,
TSP.GrowthWeighted.temperature=NA,
TSP.TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.STRNWeighted.temperature=NA,
TSP.GrowthWeighted.STRNWeighted.temperature=NA,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature=NA,
TSP.length=NA,
TSP.begin=NA,
TSP.end=NA,
TSP.PM.GrowthWeighted=NA,
TSP.PM.TimeWeighted.GrowthRateWeighted=NA,
TSP.PM.TimeWeighted=NA,
Incubation.length=NA,
MiddleThird.length=NA,
MiddleThird.begin=NA,
MiddleThird.end=NA,
MiddleThird.TimeWeighted.temperature=NA,
MiddleThird.GrowthWeighted.temperature=NA,
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.sexratio = NA,
TSP.GrowthWeighted.sexratio = NA,
TSP.TimeWeighted.GrowthRateWeighted.sexratio = NA,
TSP.TimeWeighted.STRNWeighted.sexratio = NA,
TSP.GrowthWeighted.STRNWeighted.sexratio = NA,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio = NA,
MiddleThird.TimeWeighted.sexratio = NA,
MiddleThird.GrowthWeighted.sexratio = NA,
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio = NA,
TimeWeighted.sexratio = NA,
GrowthWeighted.sexratio = NA,
TimeWeighted.GrowthRateWeighted.sexratio = NA)
df <- NULL
}
if (out == "summary") {
df <- NULL
indices.df <- NULL
}
if (out != "summary") {
summary <- NULL
}
returnlist <- list(x=list(metric=df, summary=summary, series=serie_en_cours, indices.df=indices.df))
names(returnlist) <- serie_en_cours
return(returnlist)
},
clusterEvalQ=list(expr=expression(library("embryogrowth"))),
clusterExport=list(varlist=c("temperatures_ec",
"temperature.heterogeneity",
"hatchling.metric",
"df_random",
"df_random_tsd", "df_random_SexualisationTRN",
"metabolic.heating",
"stop.at.hatchling.metric",
"embryo.stages",
"TSP.borders",
"TSP.begin",
"TSP.end",
"logicTransition",
"integral",
"derivate",
"out",
"metric.end.incubation", "tsd", "SSM"),
envir=environment())
)
if (out=="likelihood") {
return(sum(sapply(AnalyseTraces, FUN = function(x) x$likelihood)))
}
if (out=="metric") {
metric <- lapply(AnalyseTraces, FUN=function(x) x[[1]]$metric)
names(metric) <- sapply(AnalyseTraces, FUN=function(x) x[[1]]$series)
metric <- metric[name]
indices.metric <- lapply(AnalyseTraces, FUN=function(x) x[[1]]$indices.df)
names(indices.metric) <- sapply(AnalyseTraces, FUN=function(x) x[[1]]$series)
indices.metric <- indices.metric[name]
returnlist <- list(dynamic.metric=metric, indices.dynamic.metric=indices.metric)
return(returnlist)
}
# Je crée le summary
if (out == "summary") {
summary_df <- as.data.frame(t(sapply(AnalyseTraces, FUN = function(x) x[[1]]$summary)))
summary.mean <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, mean(unlist(x), na.rm=TRUE)))
colnames(summary.mean) <- c("Series", paste0(colnames(summary.mean)[-1], ".mean"))
summary.se <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, sd(unlist(x), na.rm=TRUE)))
colnames(summary.se) <- c("Series", paste0(colnames(summary.se)[-1], ".se"))
summary <- merge(x=summary.mean, y=summary.se, by.x="Series", by.y = "Series")
if (!is.null(probs)) {
for (pr in probs) {
summary.quantile <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, quantile(unlist(x), probs = pr, na.rm=TRUE)))
colnames(summary.quantile) <- c("Series", paste0(colnames(summary.quantile)[-1], paste0(".quantile_", as.character(pr))))
summary <- merge(x=summary, y=summary.quantile, by.x="Series", by.y = "Series")
}
}
rownames(summary) <- summary$Series
summary <- summary[name, ]
return(list(summary=summary))
}
# Là je suis en dynamic metric ou dynamic seul
# Je crée le dynamic.metric
if (!stop.at.hatchling.metric) {
df_list <- lapply(AnalyseTraces, FUN = function(x) x[[1]]$metric)
series_list <- sapply(AnalyseTraces, FUN = function(x) x[[1]]$series)
names(df_list) <- series_list
replicate.CI <- ifelse(replicate.CI == 0, 1, replicate.CI)
dynamic.metric <-list()
for (serie_en_cours in unique(series_list)) {
df_list_ec <- df_list[series_list == serie_en_cours]
# J'ai la série en cours dans df_list_ec
# Quel est l'intérêt de prendre df_list ici ? Il y a toutes les séries !
# Si j'ai plusieurs réplicats, ça me fait prendre toutes les séries
t_ec_count <- cbind(df_list_ec[[1]][, "Time", drop=FALSE], count=0)
for (s in seq_along(df_list_ec)) {
common_time <- t_ec_count$Time %in% df_list_ec[[s]][, "Time"]
t_ec_count[common_time, "count"] <- t_ec_count[common_time, "count"] + 1
}
t_ec <- unique(t_ec_count[t_ec_count[, "count"] == length(df_list_ec), 1])
SCL_ec <- sapply(df_list_ec, FUN = function(x) {
x[(x$Time %in% t_ec) & ((x$DeltaT != 0) | (is.na(x$R))), "SCL"]
})
SCL_ec <- as.data.frame(SCL_ec)
# Dans metric_ec j'ai les quantiles de tailles aux temps communs pour tous les réplicats
metric_ec <- t(apply(SCL_ec, MARGIN = 1, FUN = function(x) quantile(x, probs = probs)))
colnames(metric_ec) <- paste0("Metric_", colnames(metric_ec))
# Je sors les R pour tous les temps communs
R_ec <- sapply(df_list_ec, FUN = function(x) {
x[(x$Time %in% t_ec) & ((x$DeltaT != 0) | (is.na(x$R))), "R"]
})
R_ec <- as.data.frame(R_ec)
Rx_ec <- t(apply(R_ec, MARGIN = 1, FUN = function(x) quantile(x, probs = probs, na.rm = TRUE)))
colnames(Rx_ec) <- paste0("R_", colnames(Rx_ec))
# if (length(t_ec) != nrow(metric_ec)) print(serie_en_cours)
df_TC <- df_list_ec[[1]]
df_TC <- df_TC[(df_TC$Time %in% t_ec) & ((df_TC$DeltaT != 0) | (is.na(df_TC$R))), ]
metric_ec <- list(cbind(Time=t_ec, DeltaT=c(0, diff(t_ec)), TempC=df_TC[, "TempC"], metric_ec, Rx_ec))
names(metric_ec) <- serie_en_cours
dynamic.metric <- c(dynamic.metric, metric_ec)
}
summary_dynamic.metric <- as.data.frame(t(sapply(AnalyseTraces, FUN = function(x) x[[1]]$indices.df)))
summary_dynamic.metric.mean <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, mean(unlist(x), na.rm=TRUE)))
colnames(summary_dynamic.metric.mean) <- c("Series", paste0(colnames(summary_dynamic.metric.mean)[-1], ".mean"))
summary_dynamic.metric.se <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, sd(unlist(x), na.rm=TRUE)))
colnames(summary_dynamic.metric.se) <- c("Series", paste0(colnames(summary_dynamic.metric.se)[-1], ".se"))
sdm <- merge(x=summary_dynamic.metric.mean, y=summary_dynamic.metric.se, by.x="Series", by.y = "Series")
if (!is.null(probs)) {
for (pr in probs) {
summary_dynamic.metric.quantile <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, quantile(unlist(x), probs = pr, na.rm=TRUE)))
colnames(summary_dynamic.metric.quantile) <- c("Series", paste0(colnames(summary_dynamic.metric.quantile)[-1], paste0(".quantile_", as.character(pr))))
sdm <- merge(x=sdm, y=summary_dynamic.metric.quantile, by.x="Series", by.y = "Series")
}
}
} else {
sdm <- NULL
dynamic.metric <- NULL
}
dynamic.metric <- dynamic.metric[name]
rownames(sdm) <- sdm$Series
sdm <- sdm[name, ]
return(list(dynamic.metric=dynamic.metric, summary.dynamic.metric=sdm))
}
Try the embryogrowth package in your browser
Any scripts or data that you put into this service are public.
embryogrowth documentation built on Oct. 24, 2023, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions info.nests
Documented in info.nests
#' Calculate statistics about nests
#' @title Calculte statistics about nests
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return or the total likelihood or a list with $metric and $summary depending on out parameter
#' @param x A set of parameters to model the embryo growth thermal reaction norm or a NestsResult object.
#' @param parameters A set of parameters to model the embryo growth thermal reaction norm. It will replace the parameters included in NestsResult (same as x).
#' @param NestsResult A NestsResult object generated by searchR to model the embryo growth thermal reaction
#' @param resultmcmc A mcmc result for embryo growth thermal reaction norm
#' @param hessian An hessian matrix for embryo growth thermal reaction norm. It will replace the hessian matrix included in NestResult object.
#' @param GTRN.CI How to estimate CI for embryo growth thermal reaction norm; can be NULL, "SE", "MCMC", or "Hessian".
#' @param fixed.parameters A set of fixed parameters to model the embryo growth thermal reaction norm. It will replace the fixed parameters included in NestsResult.
#' @param SE Standard error for each parameter. It will replace the SE in NestsResult. Use SE=NA to remove SE from NestResult
#' @param temperatures Timeseries of temperatures formatted using formatNests(). It will replace the one in NestsResult.
#' @param integral Function used to fit embryo growth: integral.Gompertz, integral.exponential or integral.linear. It will replace the one in NestsResult.
#' @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. It will replace the one in NestsResult.
#' @param weight Weights of the different nests to estimate likelihood. It will replace the ones in NestsResult.
#' @param stop.at.hatchling.metric TRUE or FALSE. If TRUE, the model stops when proxy of size reached the mean hatchling.metric size.
#' @param M0 Measure of hatchling size proxi at laying date. It will replace the one in NestsResult.
#' @param SexualisationTRN A set of parameters used to model sexualisation thermal reaction norm during TSP or a result of STRN()
#' @param SexualisationTRN.mcmc A mcmc object obtained from STRN_MHmcmc() to generate variability for sexualisation thermal reaction norm during TSP
#' @param SexualisationTRN.CI How to estimate CI of sexualisation thermal reaction norm. Can be NULL, "SE", "MCMC", or "Hessian".
#' @param series The name or number of the series to be estimated.
#' @param TSP.borders The limits of TSP in stages. See embryo.stages parameter.
#' @param embryo.stages The embryo stages. At least TSP.borders stages must be provided to estimate TSP borders. See note.
#' @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 replicate.CI Number of replicates to estimate CI. See description
#' @param metabolic.heating Degrees Celsius to be added at the end of incubation due to metabolic heating.
#' @param temperature.heterogeneity SD of heterogeneity of temperatures. Can be 2 values, sd_low and sd_high and then HelpersMG::r2norm() is used.
#' @param fill Number of minutes between two records. Create new one if they do not exist. NULL does not change the time of temperature recordings.
#' @param probs Probabilities for metric quantiles.
#' @param out Can take the values of "likelihood", "summary", "metric" or "dynamic".
#' @param metric.end.incubation The metric at the end of incubation used to calibrate TSP size. Can be "hatchling.metric", or "observed".
#' @param progressbar If FALSE, the progress bar is not shown (useful for using with sweave or knitr)
#' @param warnings If FALSE, does not show warnings
#' @param parallel If TRUE use parallel version for nests estimation
#' @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", or "Hessian".
#' @param zero Value to replace 0 or 1.
#' @param verbose If TRUE, show more information.
#' @description This function calculates many statistics about nests.\cr
#' The embryo.stages is a named vector with relative size as compared to final size at the beginning of the stage. Names are the stages.\cr
#' For example for SCL in Caretta caretta:\cr
#' embryo.stages=structure(c(8.4, 9.4, 13.6, 13.8, 18.9, 23.5, 32.2, 35.2, 35.5, 38.5)/39.33), \cr
#' .Names = c("21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31"))\cr
#' indicates that the stages 21 begins at the relative size of 8.4/39.33 as compared to the final size.\cr
#' Series can be indicated as the name of the series, their numbers or series or succession of TRUE or FALSE. "all" indicates that all series must be analyzed.\cr
#' The likelihood object is just the total likelihood of the data in the model.\cr
#' If one parameter is named "pipping_emergence" it is used as the number of days between pipping and emergence to calculate the 1/3 and 2/3 of incubation.\cr
#' The summary object is a data.frame composed of these elements with the suffix .mean, .se or .quantile_x with x from the parameter probs.
#' \itemize{
#' \item \code{Temperature.max} Maximum temperature recorded during incubation
#' \item \code{TimeWeighted.temperature} Average temperature during all incubation
#' \item \code{GrowthWeighted.temperature} Average temperature weighted by the actual growth during all incubation
#' \item \code{TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during all incubation
#' \item \code{TSP.TimeWeighted.temperature} Average temperature during the TSP
#' \item \code{TSP.GrowthWeighted.temperature} Average temperature weighted by the actual growth during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during the TSP
#' \item \code{TSP.TimeWeighted.STRNWeighted.temperature} Average temperature weighted by the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.GrowthWeighted.STRNWeighted.temperature} Average temperature weighted by actual growth and the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature} Average temperature weighted by growth rate and the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.length} TSP duration
#' \item \code{TSP.begin} Beginning of the TSP
#' \item \code{TSP.end} End of the TSP
#' \item \code{TSP.PM.GrowthWeighted} Average of male probability for each temperature weighted by actual growth during the TSP
#' \item \code{TSP.PM.TimeWeighted.GrowthRateWeighted} Average of male probability for each temperature weighted by growth rate during the TSP
#' \item \code{TSP.PM.TimeWeighted} Average of male probability for each temperature during the TSP
#' \item \code{Incubation.length} Incubation length duration
#' \item \code{MiddleThird.length} Middle third incubation duration
#' \item \code{MiddleThird.begin} Beginning of the middle third incubation duration
#' \item \code{MiddleThird.end} End of the middle third incubation duration
#' \item \code{MiddleThird.TimeWeighted.temperature} Average temperature during the middle third incubation
#' \item \code{MiddleThird.GrowthWeighted.temperature} Average temperature weighted by the actual growth during the middle third incubation
#' \item \code{MiddleThird.TimeWeighted.GrowthRateWeighted.temperature} Average temperature weighted by the growth rate during the middle third incubation
#' \item \code{TSP.TimeWeighted.sexratio} Sex ratio based on average temperature during the TSP
#' \item \code{TSP.GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by the actual growth during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by the growth rate during the TSP
#' \item \code{TSP.TimeWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.GrowthWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the actual growth and thermal reaction norm of sexualization during the TSP
#' \item \code{TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio} Sex ratio based on average temperature weighted by the growth rate and the thermal reaction norm of sexualization during the TSP
#' \item \code{MiddleThird.TimeWeighted.sexratio} Sex ratio based on average temperature during the middle third incubation
#' \item \code{MiddleThird.GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by actual growth during the middle third incubation
#' \item \code{MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by growth rate during the middle third incubation
#' \item \code{TimeWeighted.sexratio} Sex ratio based on average temperature during all incubation
#' \item \code{GrowthWeighted.sexratio} Sex ratio based on average temperature weighted by actual growth during all incubation
#' \item \code{TimeWeighted.GrowthRateWeighted.sexratio} Sex ratio based on average temperature weighted by growth rate during all incubation
#' }
#' If \code{out} is equal to \code{summary}, the return is a list with:
#' \itemize{
#' \item \code{summary} is a data.frame with statistics for each nest.
#' \item \code{dynamic.metric} object is a list composed of data.frames with the dynamics of growth for each nest. It showed only temperatures from original dataset.
#' \item \code{summary.dynamic.metric} is a data.frame with the following columns with the suffix .mean, .se or .quantile_x with x from the parameter probs.
#' }
#' If \code{out} is equal to \code{metric}, the return is a list with:
#' \itemize{
#' \item \code{dynamic.metric} object is a list composed of data.frames with the dynamics of growth for each nest
#' \item \code{indices.dynamic.metric} is a data.frame with the following columns.
#' }
#' The object \code{summary.dynamic.metric} or \code{indices.dynamic.metric} is a data.frame with the following columns:
#' \itemize{
#' \item \code{series} Name of the series
#' \item \code{metric.begin.tsp} Metric at the beginning of TSP
#' \item \code{metric.end.tsp} Metric at the end of TSP
#' \item \code{hatchling.metric.mean} Average expected size of hatchlings
#' \item \code{hatchling.metric.sd} standard deviation of expected size of hatchlings
#' \item \code{time.begin.tsp} Time at the beginning of TSP
#' \item \code{time.end.tsp} Time at the end of TSP
#' \item \code{time.begin.middlethird} Time at the beginning of the middle third incubation
#' \item \code{time.end.middlethird} Time at the end of the middle third incubation
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' \item \code{stop.at.hatchling.metric} Take the value of NA if stop.at.hatchling.metric was FALSE. TRUE if at least one incubation series was longer than hatchling size and FALSE at contrary
#' }
#' If you indicate new set of temperatures, you must probably also indicate new hatchling.metric values.\cr
#' Note: four species have predefined embryo stages. embryo.stages parameter can take the values:\cr
#' \itemize{
#' \item \code{Caretta caretta.SCL}
#' \item \code{Chelonia mydas.SCL}
#' \item \code{Emys orbicularis.SCL}
#' \item \code{Emys orbicularis.mass}
#' \item \code{Podocnemis expansa.SCL}
#' \item \code{Lepidochelys olivacea.SCL}
#' \item \code{Generic.ProportionDevelopment}
#' }
#' But remember that mass is not the best proxy to describe the growth of an embryo because it can decrease if the substrate becomes dry.\cr
#' The progress bar is based on both replicates and timeseries progress. It necessitates the pbapply package.\cr
#' If replicate.CI is null or 0, only maximum likelihood is used and no confidence interval is calculated.\cr
#' If replicate.CI is 1, one random value for the parameters is used but no confidence interval is calculated.\cr
#' In other cases, replicate.CI random samples are used to estimate confidence interval.
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' data(resultNest_4p_SSM)
#' # Some basic calculations to show the advantage of parallel computing
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=FALSE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE, progressbar=TRUE))
#' system.time(summary.nests <- info.nests(x=resultNest_4p_SSM, out="likelihood",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=0, parallel=TRUE, progressbar=FALSE))
#'
#' # By default parallel computing is TRUE but progressbar is FALSE
#' # When out is "likelihood", it returns only the likelihood
#' # otherwise, it returns a list with 3 objects "summary",
#' # "dynamic.metric", and "summary.dynamic.metric".
#'
#' summary.nests <- info.nests(resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL",
#' replicate.CI=100,
#' resultmcmc=resultNest_mcmc_4p_SSM,
#' GTRN.CI="MCMC",
#' progressbar=TRUE)
#'
#' summary.nests <- info.nests(resultNest_4p_SSM,
#' embryo.stages="Caretta caretta.SCL",
#' out="summary", replicate.CI=100,
#' GTRN.CI="Hessian",
#' progressbar=TRUE)
#'
#' summary.nests <- info.nests(resultNest_4p_SSM,
#' series = 1,
#' embryo.stages="Caretta caretta.SCL",
#' out="summary", replicate.CI=100,
#' GTRN.CI="SE",
#' progressbar=TRUE)
#'
#' # Example of use of embryo.stages and TSP.borders:
#' summary.nests <- info.nests(resultNest_4p_SSM, out="summary",
#' embryo.stages=c("10"=0.33, "11"=0.33, "12"=0.66, "13"=0.66),
#' TSP.borders = c(10, 12),
#' replicate.CI=100,
#' progressbar=TRUE)
#'
#' #########################################
#' # Sex ratio using Massey et al. method PM
#' #########################################
#'
#' # Massey, M.D., Holt, S.M., Brooks, R.J., Rollinson, N., 2019. Measurement
#' # and modelling of primary sex ratios for species with temperature-dependent
#' # sex determination. J Exp Biol 222, 1-9.
#'
#' CC_Mediterranean <- subset(DatabaseTSD, RMU=="Mediterranean" &
#' Species=="Caretta caretta" & (!is.na(Sexed) & Sexed!=0))
#' tsdL <- with (CC_Mediterranean, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature,
#' equation="logistic", replicate.CI=NULL))
#'
#' PM <- info.nests(x=resultNest_4p_SSM,
#' GTRN.CI="Hessian", tsd.CI="Hessian",
#' embryo.stages="Caretta caretta.SCL", replicate.CI=100,
#' out="summary", progressbar=TRUE, tsd=tsdL)
#'
#'
#' plot_errbar(x=PM$summary$TimeWeighted.temperature.mean,
#' y=PM$summary$TSP.PM.GrowthWeighted.mean,
#' y.minus=PM$summary$TSP.PM.GrowthWeighted.quantile_0.025,
#' y.plus=PM$summary$TSP.PM.GrowthWeighted.quantile_0.975,
#' xlab="CTE SCL growth",
#' ylab="PM Massey et al. 2016", xlim=c(26, 32), ylim=c(0, 1), las=1)
#'
#' # Relationship between growth and growth rate
#'
#' infoall.df <- info.nests(x=resultNest_4p_SSM, out="summary",
#' embryo.stages="Caretta caretta.SCL",
#' replicate.CI=100,
#' resultmcmc=resultNest_mcmc_4p_SSM,
#' GTRN.CI="MCMC",
#' progressbar=TRUE)
#'
#' layout(1)
#' plot(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_50%"],
#' type="l", las=1, bty="n",
#' xlab="Time in minute", ylab="Growth", ylim=c(0, 39), xlim=c(0, 100000))
#' lines(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_2.5%"], lty=2)
#' lines(x=infoall.df$dynamic.metric[[1]][, "Time"],
#' y=infoall.df$dynamic.metric[[1]][, "Metric_97.5%"], lty=2)
#' }
#' @export
info.nests <- function(x=NULL ,
parameters=NULL ,
NestsResult=NULL ,
resultmcmc=NULL ,
hessian=NULL ,
GTRN.CI=NULL ,
fixed.parameters=NULL ,
SE=NULL ,
temperatures=NULL ,
integral=NULL ,
derivate=NULL ,
hatchling.metric=NULL ,
stop.at.hatchling.metric=FALSE ,
M0=NULL ,
series="all" ,
TSP.borders=NULL ,
embryo.stages=NULL ,
TSP.begin=0 ,
TSP.end=0.5 ,
replicate.CI=0 ,
weight=NULL ,
out="likelihood" ,
fill=NULL ,
probs= c(0.025, 0.5, 0.975) ,
SexualisationTRN=NULL ,
SexualisationTRN.mcmc=NULL ,
SexualisationTRN.CI=NULL ,
metric.end.incubation="observed" ,
metabolic.heating=0 ,
temperature.heterogeneity=0 ,
progressbar=FALSE ,
warnings=TRUE ,
parallel=TRUE ,
tsd=NULL ,
tsd.CI=NULL ,
tsd.mcmc=NULL ,
zero = 1E-9 ,
verbose = FALSE ) {
# parameters=NULL; temperatures=NULL; fixed.parameters=NULL
# integral=NULL; hatchling.metric=NULL; derivate=NULL
# weight=NULL; M0=NULL
# x=NULL; NestsResult=NULL; resultmcmc=NULL
# SE=NULL; hessian=NULL
# GTRN.CI=NULL; temperature.heterogeneity=0; metabolic.heating=0
# stop.at.hatchling.metric=FALSE; series="all"; TSP.borders=NULL
# embryo.stages=NULL
# TSP.begin=0; TSP.end=0.5; replicate.CI=100; out="Likelihood"
# probs= c(0.025, 0.5, 0.975)
# fill=NULL; SexualisationTRN.CI=NULL
# SexualisationTRN=NULL;SexualisationTRN.mcmc=NULL
# metric.end.incubation="observed"
# progressbar=FALSE;warnings=TRUE; parallel=TRUE
# tsd=NULL; tsd.CI=NULL; tsd.mcmc=NULL
# verbose = FALSE
# NestsResult=resultNest_Florida
# series=1:156
# embryo.stages="Caretta caretta.SCL"
# TSP.borders=NULL
# TSP.begin=0
# TSP.end=0.5
# SexualisationTRN=c('DHA' = 65.329132145227575, 'DHH' = 811.18443404950051, 'T12H' = 655.01294236106037)
# out="summary"
# replicate.CI = 0
# progressbar=FALSE
# parallel=TRUE
# x=resultNest_4p_SSM; out <- "summary"; embryo.stages="fitted"; replicate.CI=0; GTRN.CI="Hessian"; parallel=TRUE
# SexualisationTRN <- c('BeginTSP' = -1e+05, 'LengthTSP' = 0.96604492361218408)
# x=resultNest_4p_SSM; out <- "summary"; embryo.stages="Caretta caretta.SCL"; replicate.CI=100; GTRN.CI="Hessian"; parallel=TRUE
metric.end.incubation <- match.arg(arg=tolower(metric.end.incubation), choices=c("hatchling.metric", "observed"))
if ((metric.end.incubation == "observed") & stop.at.hatchling.metric) {
stop("The option metric.end.incubation == 'observed' is not compatible with stop.at.hatchling.metric being TRUE")
}
SSM <- getFromNamespace(".SSM", ns="embryogrowth")
papply <- ifelse(parallel, detectCores(), 1)
if (inherits(x, "NestsResult"))
NestsResult <- x
if (inherits(x, "numeric"))
parameters <- x
out <- tolower(out)
out <- match.arg(out, choices = c("likelihood", "metric", "summary", "dynamic"), several.ok = FALSE)
TSP.list <- embryogrowth::TSP.list
if (inherits(NestsResult, "NestsResult")) {
# temperatures est un objet Nests
if (is.null(hessian)) hessian <- NestsResult$hessian
if (is.null(temperatures)) {
temperatures <- NestsResult$data
if (is.null(hatchling.metric)) hatchling.metric <- NestsResult$hatchling.metric
}
if (is.null(integral)) integral <- NestsResult$integral
if (is.null(derivate)) derivate <- NestsResult$derivate
if (is.null(weight)) weight <- NestsResult$weight
if (is.null(M0)) M0 <- NestsResult$M0
if (is.null(fixed.parameters)) fixed.parameters <- NestsResult$fixed.parameters
if (is.null(SE)) SE <- NestsResult$SE
if (is.null(parameters)) parameters <- NestsResult$par
} else {
if (is.numeric(NestsResult) & is.null(parameters)) parameters <- NestsResult
}
if (!is.null(integral)) {
if (!is.null(derivate)) {
if (warnings) warning("Both derivate and integral cannot be supplied. Only integral parameter will be used")
derivate <- NULL
}
}
GTRN.CI <- ifelse(is.null(GTRN.CI), "null", GTRN.CI)
tsd.CI <- ifelse(is.null(tsd.CI), "null", tsd.CI)
SexualisationTRN.CI <- ifelse(is.null(SexualisationTRN.CI), "null", SexualisationTRN.CI)
if (is.null(replicate.CI)) replicate.CI <- 0
# Si je ne fais que la vraisemblance, je n'ai pas besoin de CI
if (out == "likelihood") {
if (replicate.CI != 0) {
if (warnings) warning("No replicate can be used for likelihood out.")
replicate.CI <- 0
}
}
if (replicate.CI == 0) {
GTRN.CI <- "null"
SexualisationTRN.CI <- "null"
tsd.CI <- "null"
}
if ((GTRN.CI == "null") & (tsd.CI == "null") & (SexualisationTRN.CI == "null")) {
# 15/9/2019: C'est une valeur mais sans resampling
replicate.CI <- 0
}
####################################### #
# vérifie la conformité de GTRN.CI #####
####################################### #
# GTRN.CI How to estimate CI for embryo growth thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc", or "Hessian".
GTRN.CI <- tolower(GTRN.CI)
GTRN.CI <- match.arg(GTRN.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
GTRN.CIr <- GTRN.CI
if ((GTRN.CI == "mcmc" | GTRN.CI == "pseudohessianfrommcmc") & is.null(resultmcmc)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as MCMC but no resultmcmc is available.")
}
if (!is.null(SE[1])) if (all(is.na(SE[]))) SE <- NULL
if ((GTRN.CI == "se") & is.null(SE)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as SE but no SE is available.")
}
if (GTRN.CI == "hessian" & is.null(hessian)) {
GTRN.CIr <- "null"
if (warnings) warning("GTRN.CI is indicated as Hessian but no Hessian matrix is available.")
}
GTRN.CI <- GTRN.CIr
################################################### #
# vérifie la conformité de SexualisationTRN.CI #####
################################################### #
# SexualisationTRN.CI How to estimate CI for sexualization thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
SexualisationTRN.CI <- tolower(SexualisationTRN.CI)
SexualisationTRN.CI <- match.arg(SexualisationTRN.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
SexualisationTRN.CIr <- SexualisationTRN.CI
if ((SexualisationTRN.CI == "mcmc" | SexualisationTRN.CI == "pseudohessianfrommcmc") & is.null(SexualisationTRN.mcmc)) {
SexualisationTRN.CIr <- "null"
if (warnings) warning("SexualisationTRN.CI is indicated as MCMC but no SexualisationTRN.mcmc is available.")
}
if (SexualisationTRN.CI == "hessian") {
if (is.list(SexualisationTRN)) {
if (is.null(SexualisationTRN$hessian)) {
SexualisationTRN.CIr <- "null"
}
} else {
SexualisationTRN.CIr <- "null"
}
if ((SexualisationTRN.CIr == "null") & warnings) warning("SexualisationTRN.CI is indicated as Hessian but no Hessian matrix is available.")
}
if (SexualisationTRN.CI == "se") {
if (is.list(SexualisationTRN)) {
if (is.null(SexualisationTRN$SE)) {
SexualisationTRN.CIr <- "null"
}
} else {
SexualisationTRN.CIr <- "null"
}
if ((SexualisationTRN.CIr == "null") & warnings) warning("SexualisationTRN.CI is indicated as SE but no SE is available.")
}
SexualisationTRN.CI <- SexualisationTRN.CIr
###################################### #
# vérifie la conformité de tsd.CI #####
###################################### #
# tsd.CI How to estimate CI for sex ratio thermal reaction norm; can be NULL, "SE", "MCMC", "pseudohessianfrommcmc" or "Hessian".
tsd.CI <- tolower(tsd.CI)
tsd.CI <- match.arg(tsd.CI, choices = c("null", "se", "mcmc", "hessian", "pseudohessianfrommcmc"))
tsd.CIr <- tsd.CI
if ((tsd.CI == "mcmc" | tsd.CI == "pseudohessianfrommcmc") & is.null(tsd.mcmc)) {
tsd.CIr <- "null"
if (warnings) warning("tsd.CI is indicated as MCMC but no tsd.mcmc is available.")
}
if (tsd.CI == "hessian") {
if (is.list(tsd)) {
if (is.null(tsd$hessian)) {
tsd.CIr <- "null"
}
} else {
tsd.CIr <- "null"
}
if ((tsd.CIr == "null") & warnings) warning("tsd.CI is indicated as Hessian but no Hessian matrix is available.")
}
if (tsd.CI == "se") {
if (is.list(tsd)) {
if (is.null(tsd$SE)) {
tsd.CIr <- "null"
}
} else {
tsd.CIr <- "null"
}
if ((tsd.CIr == "null") & warnings) warning("tsd.CI is indicated as SE but no SE is available.")
}
tsd.CI <- tsd.CIr
################################################### #
# Prepare logicTransition et pipping_emergence #####
################################################### #
# dans x j'ai les parametres a ajuster
# Il faut que je rajoute les fixes - 16/7/2012
totpar <- c(parameters, fixed.parameters)
if (is.na(totpar["pipping_emergence"])) fixed.parameters <- c(fixed.parameters, pipping_emergence=0)
# logictransition vaut FALSE s'il n'y en a pas
logicTransition <- !(is.na(totpar["transition_P"]) | is.na(totpar["transition_S"]))
##################################### #
# Créé la matrice des random GTRN #####
##################################### #
if (!is.null(SE)) SE <- ifelse(is.infinite(SE), 0, SE)
df_random <- RandomFromHessianOrMCMC(fitted.parameters = parameters ,
fixed.parameters = fixed.parameters ,
se = SE ,
mcmc = resultmcmc ,
Hessian = hessian ,
method = GTRN.CI ,
replicates = replicate.CI ,
chain=1 ,
regularThin=TRUE ,
MinMax=NULL ,
probs=c(0.025, 0.5, 0.975) ,
fn=NULL ,
silent=TRUE ,
ParTofn="par" )$random
##################################### #
# Créé la matrice des random STRN #####
##################################### #
if (is.list(SexualisationTRN)) {
fitted.parameters_SexualisationTRN <- SexualisationTRN$par
fixed.parameters_SexualisationTRN <- SexualisationTRN$fixed.parameters
SexualisationTRN_SE <- SexualisationTRN$SE
SexualisationTRN_Hessian <- SexualisationTRN$hessian
} else {
fitted.parameters_SexualisationTRN <- SexualisationTRN
fixed.parameters_SexualisationTRN <- NULL
SexualisationTRN_SE <- NULL
SexualisationTRN_Hessian <- NULL
}
if (is.null(fixed.parameters_SexualisationTRN) & is.null(fitted.parameters_SexualisationTRN)) {
df_random_SexualisationTRN <- NULL
SexualisationTRN.CI <- "null"
} else {
if (!is.null(SexualisationTRN_SE)) SexualisationTRN_SE <- ifelse(is.infinite(SexualisationTRN_SE), 0, SexualisationTRN_SE)
df_random_SexualisationTRN <- RandomFromHessianOrMCMC(fitted.parameters = fitted.parameters_SexualisationTRN,
fixed.parameters = fixed.parameters_SexualisationTRN,
mcmc = SexualisationTRN.mcmc,
se = SexualisationTRN_SE,
method = SexualisationTRN.CI,
Hessian = SexualisationTRN_Hessian,
replicates = replicate.CI)$random
}
##################################### #
# Créé la matrice des random tsd #####
##################################### #
tsd_fitted.parameters <- tsd$par
tsd_fixed.parameters <- tsd$fixed.parameters
if (is.null(tsd_fitted.parameters) & is.null(tsd_fixed.parameters)) {
df_random_tsd <- NULL
tsd.CI <- "null"
} else {
tsd_SE <- tsd$SE
tsd_SE <- ifelse(is.infinite(tsd_SE), 0, tsd_SE)
df_random_tsd <- RandomFromHessianOrMCMC(fitted.parameters = tsd_fitted.parameters,
fixed.parameters = tsd_fixed.parameters,
mcmc = tsd.mcmc,
se = tsd_SE,
Hessian = tsd$hessian,
method = tsd.CI,
replicates = replicate.CI)$random
}
##################################### #
# concatène les df_random #####
##################################### #
# if (!is.null(df_random_SexualisationTRN)) {
# df_random <- cbind(df_random, df_random_SexualisationTRN)
# }
# if (!is.null(df_random_tsd)) {
# df_random <- cbind(df_random, df_random_tsd)
# }
# df_random, df_random_tsd, df_random_SexualisationTRN
##################################### #
# Je sors les séries que je travaille
##################################### #
NBTs <- temperatures[["IndiceT"]]["NbTS"]
if (is.null(series[1]) | (series[1]=="all")) {
series <- rep(TRUE, NBTs)
}
# Je ne garde que les series que je dois travailler
temperatures_ec <- temperatures[1:NBTs][series]
if (any(is.na(names(temperatures_ec)))) {
stop("The names of 'series' does not exist in series of temperature within 'NestsResult'")
}
# Dans NBTs j'ai maintenant le nombre de series a travailler
NBTs <- length(temperatures_ec)
# dans name j'ai le nom des series que je vais travailler
name <- names(temperatures_ec)
names(name) <- name
# Je crée le df_random en rajoutant une colonne pour la série température
df_random <- as.data.frame(df_random)
df_random <- cbind(Series=rep(name[1], nrow(df_random)), df_random)
df_random_ec <- df_random
if (length(name) > 1) {
for (n in name[-1]) {
df_random_ec[, "Series"] <- rep(n, nrow(df_random_ec))
df_random <- rbind(df_random, df_random_ec)
}
}
df_random_tsd <- as.data.frame(df_random_tsd)
df_random_tsd <- cbind(Series=rep(name[1], nrow(df_random_tsd)), df_random_tsd)
df_random_tsd_ec <- df_random_tsd
if (length(name) > 1) {
for (n in name[-1]) {
df_random_tsd_ec[, "Series"] <- rep(n, nrow(df_random_tsd_ec))
df_random_tsd <- rbind(df_random_tsd, df_random_tsd_ec)
}
}
df_random_SexualisationTRN <- as.data.frame(df_random_SexualisationTRN)
df_random_SexualisationTRN <- cbind(Series=rep(name[1], nrow(df_random_SexualisationTRN)), df_random_SexualisationTRN)
df_random_SexualisationTRN_ec <- df_random_SexualisationTRN
if (length(name) > 1) {
for (n in name[-1]) {
df_random_SexualisationTRN_ec[, "Series"] <- rep(n, nrow(df_random_SexualisationTRN_ec))
df_random_SexualisationTRN <- rbind(df_random_SexualisationTRN, df_random_SexualisationTRN_ec)
}
}
##################################### #
# Les tailles finales attendues
##################################### #
# Je crée un dataframe dans hatchling.metric
if (inherits(hatchling.metric, "numeric")) {
hatchling.metric <- data.frame(Mean=rep(hatchling.metric["Mean"], NBTs),
SD=rep(hatchling.metric["SD"], NBTs), row.names=name)
}
if (is.null(hatchling.metric) | any(is.na(hatchling.metric[name, "Mean"]))) {
# 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])
hatchling.metric <- data.frame(Mean=rep(hatchling.metric["Mean"], NBTs),
SD=rep(hatchling.metric["SD"], NBTs), row.names=name)
} else {
stop("The hatchlings size (hatchling.metric parameter) must be provided for all series.")
}
}
########################################################################################## #
# J'ai donc 1 tableau: df_random avec une ligne par réplicat parameters/series.
# J'ai les tailles finales attendues dans un data.frame hatchling.metric avec deux colonnes #
# J'ai aussi les tailles finales dans SCL.final ? A vérifier #
# Dans name j'ai le nom des séries à travailler #
########################################################################################## #
# The theorem of Cauchy–Lipschitz (also known as the theorem of Picard–Lindelöf)
# gives a set of conditions under which an initial value problem (also named Cauchy problem) has a unique solution.
# Je dois calculer la température moyenne pondérée par le temps d'incubation
# Les températures sont dans temperatures. C'est une liste formattée
# Les séries à travailler sont dans name
# temperature[name]
# Pour chacune je calcule la température moyenne
# Mais cela n'a de sens que si je n'ai pas de metabolic heating ou temperature.heterogeneity et si j'ai bien une taille finale attendue donc !stop.at.hatchling.metric
################################### #
# Lecture des stades pour la TSP ####
################################### #
# Comme maintenant c'est une variable, je dois le rentrer dans la boucle
if (inherits(embryo.stages, "character")) {
if (embryo.stages == "fitted") {
embryo.stages = c("10"=invlogit(logit(0.33)), "11"=invlogit(logit(0.33)),
"12"=invlogit(logit(0.66)), "13"=invlogit(logit(0.66)))
TSP.borders = c(10, 12)
} else {
estages <- TSP.list[[gsub(" ", "_", embryo.stages)]]
if (is.null(estages)) {
stop("The TSP for ", embryo.stages, " does not exist")
} else {
embryo.stages <- estages[, "metric"]
names(embryo.stages) <- estages[, "stages"]
TSP.borders <- c(attributes(estages)$TSP.begin.stages, attributes(estages)$TSP.end.stages)
}
}
}
# dans name j'ai les noms des series
# prepare le data.frame pour les resultats ####
temperatures_ec <- universalmclapply(X=name, mc.cores = papply, progressbar=FALSE, FUN=function (serie_en_cours) {
df <- as.data.frame(temperatures_ec[[serie_en_cours]])
colnames(df) <- c("Time", "TempC", "TempK", "R", "SCL", "IndiceK")
if (df[1, "Time"] != 0) {
df <- rbind(df[1, ], df)
df[1, "Time"] <- 0
}
if (!is.null(fill)) {
newt <- seq(from=0, to=tail(df$Time, n=1L), by=fill)
aretirer <- -which(!is.na(match(newt, df$Time)))
if (!identical(aretirer, integer(0))) newt <- newt[aretirer]
if (!identical(newt, numeric(0))) {
df <- rbind(df,
data.frame(Time=newt, TempC=NA,
TempK=NA, R=NA, SCL=NA, IndiceK=NA))
df <- df[order(df$Time),]
}
}
df <- cbind(df, DeltaT=c(diff(df$Time), 0))
ind <- which(!is.na(df$TempC)) # get positions of non-missing values
if(is.na(df$TempC[1])) { # if it begins with a missing, add the
df$TempC[1] <- df$TempC[ind[1]]
df$TempK[1] <- df$TempK[ind[1]]
df$IndiceK[1] <- df$IndiceK[ind[1]]
ind <- c(1,ind)
}# first position to the indices
df$TempC <- rep(df$TempC[ind], times = diff( # repeat the values at these indices
c(ind, length(df$TempC) + 1) ))
df$TempK <- rep(df$TempK[ind], times = diff( # repeat the values at these indices
c(ind, length(df$TempK) + 1) ))
df$IndiceK <- rep(df$IndiceK[ind], times = diff( # repeat the values at these indices
c(ind, length(df$IndiceK) + 1) ))
df[, "SCL"] <- NA
df[1, "SCL"] <- M0
df[, "IndiceK"] <- 1
return(df)
}, clusterExport=list(varlist=c("temperatures_ec", "fill", "M0"),
envir = environment()),
clusterEvalQ=list(expr=expression(library(embryogrowth)))
)
# Je les classe
temperatures_ec <- temperatures_ec[name]
################################################################### #
# Je prépare le calcul avec tous les nids et tous les réplicats ####
################################################################### #
AnalyseTraces <- universalmclapply(progressbar=progressbar, mc.cores = papply, X=1:nrow(df_random),
FUN=function (ligne_df_random) {
# ligne_df_random <- 1
# print(ligne_df_random)
serie_en_cours <- df_random[ligne_df_random, "Series"]
df <- temperatures_ec[[serie_en_cours]]
setpar <- c(unlist(df_random[ligne_df_random, -1, drop=TRUE]), unlist(df_random_tsd[ligne_df_random, -1, drop=TRUE]))
setpar_STRN <- unlist(df_random_SexualisationTRN[ligne_df_random, -1, drop=TRUE])
# Si je n'ai pas setpar_STRN, je le mets en uniforme; dans dbeta_mu, j'ai le logit de mu (0=0.5)
if (is.null(setpar_STRN)) setpar_STRN <- c(dbeta_mu=0, dbeta_v=1/12)
if (all(substr(names(setpar_STRN), 1, 6) != "dbeta_")) {
setpar_STRN["dbeta_mu"] <- 0
setpar_STRN["dbeta_v"] <- 1/12
}
transition_P <- setpar["transition_P"]
transition_S <- setpar["transition_S"]
# "df_random_tsd", "df_random_SexualisationTRN"
if (!is.na(setpar["K"])) {
Kval <- unname(setpar["K"])
} else {
Kval <- NULL
}
if (length(temperature.heterogeneity)==1) {
tht <- 0
if (temperature.heterogeneity !=0) tht <- rnorm(1, 0, temperature.heterogeneity)
} else {
tht <- r2norm(1, 0, sd_low = temperature.heterogeneity[1], sd_high = temperature.heterogeneity[2])
}
df[, "TempC"] <- df[, "TempC"] + tht
df[, "TempK"] <- df[, "TempK"] + tht
meanSCL <- as.numeric(hatchling.metric[serie_en_cours, "Mean"])
sdSCL <- as.numeric(hatchling.metric[serie_en_cours, "SD"])
if (!is.na(setpar["rK"])) {
Kval <- unname(setpar["rK"]*meanSCL)
}
tmin <- df[, "Time"]
for (j in 2:nrow(df)) if (is.na(df[j, "TempK"])) df[j, "TempK"] <- df[j-1, "TempK"]
dT <- df[, "DeltaT"]
y0 <- ypre <- df[1, "SCL"]
if (!stop.at.hatchling.metric) {
if (!logicTransition) {
# if ((metabolic.heating == 0) &
# ((length(temperature.heterogeneity) == 1) &
# (temperature.heterogeneity[1] ==0)) & (out=="likelihood")) {
# indice.boucle <- nrow(df) - 1
# } else {
# indice.boucle <- 1
# }
indice.boucle <- 1
# Je prends de 1 jusqu'à l'avant-dernier
for (i in indice.boucle:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
a_inst <- SSM(CTE_inst, setpar, verbose=verbose)[[1]]
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
}
df[nrow(df), "TempK"] <- df[nrow(df), "TempK"] + metabolic.heating*(ypre/meanSCL)
} else {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
transition <- 1/(1+exp(transition_S*(ypre-transition_P)))
rT_inst < SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
}
df[nrow(df), "TempK"] <- df[nrow(df), "TempK"] + metabolic.heating*(ypre/meanSCL)
}
# je suis en stop.at.hatchling.metric
} else {
if (!logicTransition) {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
a_inst <- SSM(CTE_inst, setpar, verbose=verbose)[[1]]
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
if (ypre>meanSCL) {
break
}
}
} else {
for (i in 1:(nrow(df)-1)) {
timesunique <- c(tmin[i], tmin[i+1])
df[i, "TempK"] <- df[i, "TempK"]+metabolic.heating*(ypre/meanSCL)
CTE_inst <- df[i, "TempK"]
transition <- 1/(1+exp(transition_S*(ypre-transition_P)))
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df[i, "R"] <- a_inst
df[i+1, "SCL"] <- ypre
if (ypre>meanSCL) {
break
}
}
# fin du else logictransition
}
# fin du else !stop.at.hatchling.metric
}
df[, "TempC"] <- df[, "TempK"] - 273.15
if (stop.at.hatchling.metric) {
size.final <- ifelse(ypre < meanSCL, NA, meanSCL)
} else {
size.final <- ifelse(metric.end.incubation == "observed", ypre, meanSCL)
}
if (!is.na(size.final)) {
if (!is.null(embryo.stages)) {
if (!is.na(setpar_STRN["BeginTSP"])) {
xxbt <- invlogit(setpar_STRN["BeginTSP"])
if (xxbt == 0) xxbt <- zero
if (xxbt > 0.95) xxbt <- 0.95
embryo.stages[as.character(TSP.borders[1])] <- xxbt
embryo.stages[as.character(TSP.borders[1]+1)] <- xxbt
}
if (!is.na(setpar_STRN["EndTSP"])) {
xxbt <- invlogit(setpar_STRN["EndTSP"])
if (xxbt == 0) xxbt <- zero
if (xxbt >= 0.99) xxbt <- 0.99
embryo.stages[as.character(TSP.borders[2])] <- xxbt
embryo.stages[as.character(TSP.borders[2]+1)] <- xxbt
}
if (!is.na(setpar_STRN["LengthTSP"])) {
xxbt <- invlogit(setpar_STRN["BeginTSP"]+abs(setpar_STRN["LengthTSP"]))
if (xxbt == 0) xxbt <- zero
if (xxbt >= 0.99) xxbt <- 0.99
embryo.stages[as.character(TSP.borders[2])] <- xxbt
embryo.stages[as.character(TSP.borders[2]+1)] <- xxbt
}
# Je ne comprends pas le mean ici - 2021-05-12
size.begin.TSP <- mean(embryo.stages[as.character(TSP.borders[1])])
if (TSP.begin != 0) {
size.begin.TSP <- size.begin.TSP+(embryo.stages[as.character(TSP.borders[1]+1)]-size.begin.TSP)*TSP.begin
}
size.end.TSP <- mean(embryo.stages[as.character(TSP.borders[2])])
if (TSP.end != 0) {
size.end.TSP <- size.end.TSP+(embryo.stages[as.character(TSP.borders[2]+1)]-size.end.TSP)*TSP.end
}
size.begin.TSP_ec <- size.begin.TSP*size.final
size.end.TSP_ec <- size.end.TSP*size.final
} else {
size.end.TSP <- NA
size.begin.TSP <- NA
size.begin.TSP_ec <- NA
size.end.TSP_ec <- NA
}
if (out == "likelihood") {
if (!is.na(sdSCL)) {
likelihood <- (-dnorm(ypre, mean=meanSCL, sd=sdSCL, log=TRUE))
} else {
if (is.na(parameters["SD"])) {
likelihood <- ((ypre - meanSCL)^2)
} else {
# 28/3/2022
likelihood <- (-dnorm(ypre, mean=meanSCL, sd=parameters["SD"], log=TRUE))
}
}
return(list(likelihood=likelihood))
}
df_ec <- df
##### Estimate limit begin TSP #####
if ((size.final < size.begin.TSP_ec) | is.na(size.begin.TSP_ec) | (max(df_ec[,"SCL"], na.rm = TRUE) < size.begin.TSP_ec)) {
time.begin.TSP <- NA
indice.begin.tsp <- NA
} else {
if (size.begin.TSP_ec < df_ec[1,"SCL"]) {
size.begin.TSP_ec <- df_ec[1,"SCL"]+(df_ec[2,"SCL"] - df_ec[1,"SCL"])/10
}
indice.begin.tsp <- which(df_ec[,"SCL"]>(size.begin.TSP_ec))[1]-1
df_ec <- rbind(df_ec[1:indice.begin.tsp, ],
data.frame(Time=c(df_ec[indice.begin.tsp,"Time"]+df_ec[indice.begin.tsp,"DeltaT"]/3,
df_ec[indice.begin.tsp,"Time"]+df_ec[indice.begin.tsp,"DeltaT"]*2/3),
TempC=c(df_ec[indice.begin.tsp,"TempC"],
df_ec[indice.begin.tsp,"TempC"]),
TempK=c(df_ec[indice.begin.tsp,"TempK"],
df_ec[indice.begin.tsp,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.begin.tsp+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
df_ec[indice.begin.tsp, "DeltaT"] <- df_ec[indice.begin.tsp+1,"Time"]-df_ec[indice.begin.tsp,"Time"]
df_ec[indice.begin.tsp+1, "DeltaT"] <- df_ec[indice.begin.tsp+2,"Time"]-df_ec[indice.begin.tsp+1,"Time"]
df_ec[indice.begin.tsp+2, "DeltaT"] <- df_ec[indice.begin.tsp+3,"Time"]-df_ec[indice.begin.tsp+2,"Time"]
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
CTE_inst <- df_ec[indice.begin.tsp, "TempK"]
timesunique <- c(df_ec[indice.begin.tsp, "Time"], df_ec[indice.begin.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.begin.tsp], setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.tsp, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+1, "R"] <- a_inst
df_ec[indice.begin.tsp+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.begin.tsp+1, "Time"], df_ec[indice.begin.tsp+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
CTE_inst <- df_ec[indice.begin.tsp+1, "TempK"]
# rT <- SSM(CTE[indice.begin.tsp+1], setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre,
parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+2, "R"] <- a_inst
df_ec[indice.begin.tsp+2, "SCL"] <- ypre
# J'ai encadré le début de la TSP
t <- df_ec[indice.begin.tsp:(indice.begin.tsp+3), "Time"]
SCL <- df_ec[indice.begin.tsp:(indice.begin.tsp+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.begin.tsp <- predict(l, newdata = data.frame(SCL=size.begin.TSP_ec))
time.begin.tsp <- round(time.begin.tsp, 2)
if (time.begin.tsp <= 0) {
time.begin.tsp <- t[1]+zero
}
indice.begin.tsp <- which(df_ec[,"Time"]>(time.begin.tsp))[1]-1
df_ec <- rbind(df_ec[1:indice.begin.tsp, ],
data.frame(Time=time.begin.tsp,
TempC=df_ec[indice.begin.tsp,"TempC"],
TempK=df_ec[indice.begin.tsp,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.begin.tsp+1):nrow(df_ec), ])
# df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
df_ec[indice.begin.tsp, "DeltaT"] <- df_ec[indice.begin.tsp+1, "Time"]-df_ec[indice.begin.tsp, "Time"]
df_ec[indice.begin.tsp+1, "DeltaT"] <- df_ec[indice.begin.tsp+2, "Time"]-df_ec[indice.begin.tsp+1, "Time"]
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
timesunique <- c(df_ec[indice.begin.tsp, "Time"], df_ec[indice.begin.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.begin.tsp], setpar)
CTE_inst <- df_ec[indice.begin.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.tsp, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.tsp+1, "R"] <- a_inst
df_ec[indice.begin.tsp+1, "SCL"] <- ypre
indice.begin.tsp <- indice.begin.tsp+1
rownames(df_ec)[indice.begin.tsp] <- "Begin TSP"
time.begin.TSP <- df_ec[indice.begin.tsp,"Time"]
}
###### Estimate limit end TSP ####
if ((size.final < size.end.TSP_ec) | is.na(size.end.TSP_ec) | (max(df_ec[,"SCL"], na.rm = TRUE) < size.end.TSP_ec)) {
time.end.TSP <- NA
indice.end.tsp <- NA
} else {
if (size.end.TSP_ec < df_ec[1,"SCL"]) {
size.end.TSP_ec <- df_ec[1,"SCL"]+(df_ec[2,"SCL"] - df_ec[1,"SCL"])/10
}
ind_nonNA <- max(which(!is.na(df_ec[, "SCL"])))
if (size.end.TSP_ec >= df_ec[ind_nonNA,"SCL"]) {
size.end.TSP_ec <- df_ec[ind_nonNA-1,"SCL"] + (df_ec[ind_nonNA,"SCL"]-df_ec[ind_nonNA-1,"SCL"])/10
}
if (size.end.TSP_ec <= df_ec[indice.begin.tsp, "SCL"]) {
size.end.TSP_ec <- df_ec[indice.begin.tsp, "SCL"]+(df_ec[indice.begin.tsp+1, "SCL"]-df_ec[indice.begin.tsp, "SCL"])/10
}
indice.end.tsp <- which(df_ec[,"SCL"]>(size.end.TSP_ec))[1]-1
if (is.na(indice.end.tsp)) indice.end.tsp <- nrow(df_ec)
# if (is.na(indice.end.tsp)) {
# print("AAAA 1154")
# print(paste0("size.end.TSP_ec", as.character(size.end.TSP_ec)))
# save(df_ec, file="df_ec.Rdata")
# stop()
# }
df_ec <- rbind(df_ec[1:indice.end.tsp, ],
data.frame(Time=c(df_ec[indice.end.tsp,"Time"]+df_ec[indice.end.tsp,"DeltaT"]/3,
df_ec[indice.end.tsp,"Time"]+df_ec[indice.end.tsp,"DeltaT"]*2/3),
TempC=c(df_ec[indice.end.tsp,"TempC"],
df_ec[indice.end.tsp,"TempC"]),
TempK=c(df_ec[indice.end.tsp,"TempK"],
df_ec[indice.end.tsp,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.end.tsp+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# sum(dt[1:(indice.end.tsp+1)]*TK[1:(indice.end.tsp+1)])/sum(dt[1:(indice.end.tsp+1)])
# CTE <- cumsum(dt*TK)/cumsum(dt)
df_ec[indice.end.tsp, "DeltaT"] <- df_ec[indice.end.tsp+1,"Time"]-df_ec[indice.end.tsp,"Time"]
df_ec[indice.end.tsp+1, "DeltaT"] <- df_ec[indice.end.tsp+2,"Time"]-df_ec[indice.end.tsp+1,"Time"]
df_ec[indice.end.tsp+2, "DeltaT"] <- df_ec[indice.end.tsp+3,"Time"]-df_ec[indice.end.tsp+2,"Time"]
timesunique <- c(df_ec[indice.end.tsp, "Time"], df_ec[indice.end.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.tsp, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.end.tsp], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.end.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.tsp, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.tsp+1, "R"] <- a_inst
df_ec[indice.end.tsp+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.end.tsp+1, "Time"], df_ec[indice.end.tsp+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
# rT <- SSM(CTE[indice.end.tsp+1], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.end.tsp+1, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre,
parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.tsp+2, "R"] <- a_inst
df_ec[indice.end.tsp+2, "SCL"] <- ypre
t <- df_ec[indice.end.tsp:(indice.end.tsp+3), "Time"]
SCL <- df_ec[indice.end.tsp:(indice.end.tsp+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.end.tsp <- predict(l, newdata = data.frame(SCL=size.end.TSP_ec))
time.end.tsp <- round(time.end.tsp, 2)
# plot(t, SCL); points(time.begin.tsp, size.begin.TSP_ec, col="red"); segments(x0=0, y0=size.begin.TSP_ec, x1=100000, y1=size.begin.TSP_ec, col="green")
indice.end.tsp <- which(df_ec[,"Time"]>(time.end.tsp))[1]-1
# if (is.na(indice.end.tsp)) {
# indice.end.tsp <- nrow(df_ec)
#
# df_ec <- rbind(df_ec, df_ec[nrow(df_ec), ])
# } else {
#
df_ec <- rbind(df_ec[1:indice.end.tsp, ],
data.frame(Time=time.end.tsp,
TempC=df_ec[indice.end.tsp,"TempC"],
TempK=df_ec[indice.end.tsp,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.end.tsp+1):nrow(df_ec), ])
df_ec[indice.end.tsp, "DeltaT"] <- df_ec[indice.end.tsp+1, "Time"]-df_ec[indice.end.tsp, "Time"]
df_ec[indice.end.tsp+1, "DeltaT"] <- df_ec[indice.end.tsp+2, "Time"]-df_ec[indice.end.tsp+1, "Time"]
# }
timesunique <- c(df_ec[indice.end.tsp, "Time"], df_ec[indice.end.tsp+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.tsp, "SCL"] - transition_P)))
}
CTE_inst <- df_ec[indice.end.tsp, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.tsp, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.tsp, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
indice.end.tsp <- indice.end.tsp+1
df_ec[indice.end.tsp, "R"] <- a_inst
df_ec[indice.end.tsp, "SCL"] <- ypre
rownames(df_ec)[indice.end.tsp] <- "End TSP"
time.end.TSP <- df_ec[indice.end.tsp,"Time"]
}
#### End of incubation for stop.at.hatchling.metric ##########
# attSAT vaut NA si je ne suis pas en stop.at.hatchling.metric
# attSAT vaut TRUE si l'incubation a été jusquà meanSCL
# attSAT vaut FALSE si l'incubation n'a pas été jusqu'à meanSCL
if (stop.at.hatchling.metric) {
if (max(df_ec[,"SCL"], na.rm=TRUE) > meanSCL) {
attSAT <- TRUE
indice.fin.incubation <- which(df_ec[,"SCL"]> meanSCL)[1]-1
df_ec <- rbind(df_ec[1:indice.fin.incubation, ],
data.frame(Time=c(df_ec[indice.fin.incubation,"Time"]+df_ec[indice.fin.incubation,"DeltaT"]/3,
df_ec[indice.fin.incubation,"Time"]+df_ec[indice.fin.incubation,"DeltaT"]*2/3),
TempC=c(df_ec[indice.fin.incubation,"TempC"],
df_ec[indice.fin.incubation,"TempC"]),
TempK=c(df_ec[indice.fin.incubation,"TempK"],
df_ec[indice.fin.incubation,"TempK"]),
R=c(NA, NA),
SCL=c(NA, NA),
IndiceK=c(0, 0),
DeltaT=c(NA, NA)),
df_ec[(indice.fin.incubation+1):nrow(df_ec), ])
# df_ec[, 'DeltaT'] <- c(diff(df_ec[, "Time"]), 0)
#
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
#
# CTE <- cumsum(dt*TK)/cumsum(dt)
df_ec[indice.fin.incubation, "DeltaT"] <- df_ec[indice.fin.incubation+1,"Time"]-df_ec[indice.fin.incubation,"Time"]
df_ec[indice.fin.incubation+1, "DeltaT"] <- df_ec[indice.fin.incubation+2,"Time"]-df_ec[indice.fin.incubation+1,"Time"]
df_ec[indice.fin.incubation+2, "DeltaT"] <- df_ec[indice.fin.incubation+3,"Time"]-df_ec[indice.fin.incubation+2,"Time"]
timesunique <- c(df_ec[indice.fin.incubation, "Time"], df_ec[indice.fin.incubation+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.fin.incubation, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.fin.incubation, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.fin.incubation, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+1, "R"] <- a_inst
df_ec[indice.fin.incubation+1, "SCL"] <- ypre
timesunique <- c(df_ec[indice.fin.incubation+1, "Time"], df_ec[indice.fin.incubation+2, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(ypre - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation+1], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation+1, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=ypre, parms=param)[1, "metric"]
} else {
out1 <- lsoda(ypre, timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+2, "R"] <- a_inst
df_ec[indice.fin.incubation+2, "SCL"] <- ypre
t <- df_ec[indice.fin.incubation:(indice.fin.incubation+3), "Time"]
SCL <- df_ec[indice.fin.incubation:(indice.fin.incubation+3), "SCL"]
l <- lm(t ~ poly(SCL, degree = 2))
time.fin.incubation <- predict(l, newdata = data.frame(SCL=meanSCL))
time.fin.incubation <- round(time.fin.incubation, 2)
# plot(t, SCL); points(time.begin.tsp, size.begin.TSP_ec, col="red"); segments(x0=0, y0=size.begin.TSP_ec, x1=100000, y1=size.begin.TSP_ec, col="green")
indice.fin.incubation <- which(df_ec[,"Time"]>(time.fin.incubation))[1]-1
df_ec <- rbind(df_ec[1:indice.fin.incubation, ],
data.frame(Time=time.fin.incubation,
TempC=df_ec[indice.fin.incubation,"TempC"],
TempK=df_ec[indice.fin.incubation,"TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.fin.incubation+1):nrow(df_ec), ])
df_ec[indice.fin.incubation, "DeltaT"] <- df_ec[indice.fin.incubation+1, "Time"]-df_ec[indice.fin.incubation, "Time"]
df_ec[indice.fin.incubation+1, "DeltaT"] <- df_ec[indice.fin.incubation+2, "Time"]-df_ec[indice.fin.incubation+1, "Time"]
# df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
#
# dt <- df_ec[, "DeltaT"]
# TK <- df_ec[, "TempK"]
# CTE <- cumsum(dt*TK)/cumsum(dt)
timesunique <- c(df_ec[indice.fin.incubation, "Time"], df_ec[indice.fin.incubation+1, "Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.fin.incubation, "SCL"] - transition_P)))
}
# rT <- SSM(CTE[indice.fin.incubation], setpar)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
CTE_inst <- df_ec[indice.fin.incubation, "TempK"]
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.fin.incubation, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.fin.incubation, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.fin.incubation+1, "R"] <- a_inst
df_ec[indice.fin.incubation+1, "SCL"] <- ypre
indice.fin.incubation <- indice.fin.incubation+1
rownames(df_ec)[indice.fin.incubation] <- "End Incubation"
time.end.incubation <- df_ec[indice.fin.incubation,"Time"]
# Maintenant je tronque df de faon ne garder que jusqu'a time.end.incubation
# garde
df_ec <- df_ec[1:indice.fin.incubation,]
# Je garde la fin de l'incubation
df_ec[indice.fin.incubation, "IndiceK"] <- 1
df_ec[, "DeltaT"] <- c(diff(df_ec[, "Time"]), 0)
} else {
attSAT <- FALSE
indice.fin.incubation <- max(which(!is.na(df_ec[, "SCL"])))
rownames(df_ec)[indice.fin.incubation] <- "End Incubation"
time.end.incubation <- df_ec[indice.fin.incubation,"Time"]
df_ec[indice.fin.incubation, "IndiceK"] <- 1
}
} else {attSAT <- NA}
################################## #
# indice.begin.middlethird ####
################################## #
# On peut la calculer seulement si attSAT == TRUE ou NA
if (attSAT | is.na(attSAT)) {
time.begin.middlethird <- unname((tail(df_ec$Time, n=1L)-setpar["pipping_emergence"])/3)
indice.begin.middlethird <- which(df_ec[,"Time"]>time.begin.middlethird)[1]-1
df_ec <- rbind(df_ec[1:indice.begin.middlethird, ],
data.frame(Time=time.begin.middlethird,
TempC=df_ec[indice.begin.middlethird, "TempK"]-273.15,
TempK=df_ec[indice.begin.middlethird, "TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT= NA),
df_ec[(indice.begin.middlethird+1):nrow(df_ec), ])
df_ec[indice.begin.middlethird,"DeltaT"] <- df_ec[indice.begin.middlethird + 1,"Time"] - df_ec[indice.begin.middlethird,"Time"]
df_ec[indice.begin.middlethird + 1,"DeltaT"] <- df_ec[indice.begin.middlethird + 2,"Time"] - df_ec[indice.begin.middlethird + 1,"Time"]
# sum(df_ec[1:10,"DeltaT"]); df_ec[10+1,"Time"]
# sum(df_ec[1:indice.begin.middlethird,"DeltaT"]); df_ec[indice.begin.middlethird+1,"Time"]
# timei1 <- df_ec[indice.begin.middlethird + 1,"Time"]
# dt <- df_ec[1:indice.begin.middlethird, "DeltaT"]
# TK <- df_ec[1:indice.begin.middlethird, "TempK"]
# CTE <- sum(dt*TK)/sum(dt)
CTE_inst <- df_ec[indice.begin.middlethird, "TempK"]
timesunique <- c(df_ec[indice.begin.middlethird,"Time"], df_ec[indice.begin.middlethird + 1,"Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.begin.middlethird, "SCL"] - transition_P)))
}
# rT <- SSM(CTE, setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.begin.middlethird, "SCL"], parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.begin.middlethird, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.begin.middlethird+1, "R"] <- a_inst
df_ec[indice.begin.middlethird+1, "SCL"] <- ypre
rownames(df_ec)[indice.begin.middlethird+1] <- "Begin Middle-third"
size.begin.middlethird <- ypre
} else {
size.begin.middlethird <- NA
time.begin.middlethird <- NA
}
################################ #
# indice.end.middlethird ####
################################ #
# On peut la calculer seulement si attSAT == TRUE ou NA
if (attSAT | is.na(attSAT)) {
time.end.middlethird <- unname((tail(df_ec$Time, n=1L)-setpar["pipping_emergence"])*2/3)
indice.end.middlethird <- which(df_ec[,"Time"]>time.end.middlethird)[1]-1
# Je calcul la taille
df_ec <- rbind(df_ec[1:indice.end.middlethird, ],
data.frame(Time=time.end.middlethird,
TempC=df_ec[indice.end.middlethird, "TempK"]-273.15,
TempK=df_ec[indice.end.middlethird, "TempK"],
R=NA,
SCL=NA,
IndiceK=0,
DeltaT=NA),
df_ec[(indice.end.middlethird+1):nrow(df_ec), ])
df_ec[indice.end.middlethird,"DeltaT"] <- df_ec[indice.end.middlethird + 1,"Time"] - df_ec[indice.end.middlethird,"Time"]
df_ec[indice.end.middlethird + 1,"DeltaT"] <- df_ec[indice.end.middlethird + 2,"Time"] - df_ec[indice.end.middlethird + 1,"Time"]
# timei1 <- df_ec[indice.end.middlethird+1,"Time"]
# dt <- df_ec[1:indice.end.middlethird, "DeltaT"]
# TK <- df_ec[1:indice.end.middlethird, "TempK"]
# CTE <- sum(dt*TK)/sum(dt)
CTE_inst <- df_ec[indice.end.middlethird, "TempK"]
timesunique <- c(df_ec[indice.end.middlethird,"Time"], df_ec[indice.end.middlethird+1,"Time"])
if (!logicTransition) {
transition <- 1
} else {
transition <- 1/(1+exp(transition_S*(df_ec[indice.end.middlethird, "SCL"] - transition_P)))
}
# rT <- SSM(CTE, setpar)
rT_inst <- SSM(CTE_inst, setpar, verbose=verbose)
# a <- rT[[1]]*transition+rT[[2]]*(1-transition)
a_inst <- rT_inst[[1]]*transition+rT_inst[[2]]*(1-transition)
param <- c(alpha=unname(a_inst), K=Kval)
if (!is.null(integral)) {
ypre <- integral(t=timesunique[2]-timesunique[1], size=df_ec[indice.end.middlethird, "SCL"],
parms=param)[1, "metric"]
} else {
out1 <- lsoda(df_ec[indice.end.middlethird, "SCL"], timesunique, derivate, param)
ypre <- as.numeric(tail(out1[,2], n=1))
}
df_ec[indice.end.middlethird+1, "R"] <- a_inst
df_ec[indice.end.middlethird+1, "SCL"] <- ypre
rownames(df_ec)[indice.end.middlethird+1] <- "End Middle-third"
size.end.middlethird <- ypre
} else {
size.end.middlethird <- NA
time.end.middlethird <- NA
}
df <- df_ec
# df[, "DeltaT"] <- c(diff(df[, "Time"]), 0)
indice.begin.middlethird <- which(rownames(df) == "Begin Middle-third")
if (identical(indice.begin.middlethird, integer(0))) indice.begin.middlethird <- NA
indice.end.middlethird <- which(rownames(df) == "End Middle-third")
if (identical(indice.end.middlethird, integer(0))) indice.end.middlethird <- NA
indice.fin.incubation <- nrow(df)
indice.begin.tsp <- which(rownames(df) == "Begin TSP")
indice.end.tsp <- which(rownames(df) == "End TSP")
if (identical(integer(0), indice.end.tsp)) indice.end.tsp <- NA
if (identical(integer(0), indice.begin.tsp)) indice.begin.tsp <- NA
## Je viens de calculer la croissance de l'embryon
if (out == "summary") {
# Maintenant je calcule les summary
TC <- df[1:(nrow(df)-1), "TempC"]
dt <- df[1:(nrow(df)-1), "DeltaT"]
dSCL <- diff(df[, "SCL"])
rT <- df[1:(nrow(df)-1), "R"]
Temperature.max <- max(TC)
TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
Incubation.length <- df[nrow(df), "Time"]
if (!is.na(indice.begin.tsp) & !is.na(indice.end.tsp)) {
TC <- df[indice.begin.tsp:(indice.end.tsp-1), "TempC"]
dt <- df[indice.begin.tsp:(indice.end.tsp-1), "DeltaT"]
dSCL <- diff(df[indice.begin.tsp:(indice.end.tsp), "SCL"])
rT <- df[indice.begin.tsp:(indice.end.tsp-1), "R"]
TSP.TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
TSP.GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
TSP.TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
TSP.length <- unname(time.end.tsp - time.begin.tsp)
TSP.begin <- unname(time.begin.tsp)
TSP.end <- unname(time.end.tsp)
if (!is.null(setpar_STRN)) {
if (indice.begin.tsp == indice.end.tsp-1) {
indice.end.tsp_ec <- indice.end.tsp+1
} else {
indice.end.tsp_ec <- indice.end.tsp
}
TK <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "TempK"]
TC <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "TempC"]
dt <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "DeltaT"]
dSCL <- diff(df[indice.begin.tsp:(indice.end.tsp_ec), "SCL"])
rT <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "R"]
# print("setpar_STRN")
# print(setpar_STRN)
if (any(is.finite(suppressWarnings(as.numeric(names(setpar_STRN)))))) {
sparSTRN <- setpar_STRN[!is.na(suppressWarnings(as.numeric(names(setpar_STRN))))]
# print("sparSTRN")
# print(sparSTRN)
STRN_model <- SSM(TK, sparSTRN, verbose=verbose)[[1]]
STRN_model <- ifelse(STRN_model <= 0, 0.001, STRN_model)
# print("SSM")
# print(STRN_model)
} else {
if (any(names(setpar_STRN) == "Rho25") | any(names(setpar_STRN) == "Peak")) {
STRN_model <- SSM(TK, setpar_STRN, verbose=verbose)[[1]]
} else {
#print("setpar_STRN")
# print(setpar_STRN)
STRN_model <- rep(1, length(TK))
# print("SSM")
# print(STRN_model)
}
}
SCL <- df[indice.begin.tsp:(indice.end.tsp_ec-1), "SCL"]
SCL <- (SCL-SCL[1])/(SCL[length(SCL)]-SCL[1])
SCL[SCL==0] <- min(zero, min(SCL[SCL != 0])/10)
SCL[SCL==1] <- max(1-zero, max(SCL[SCL != 1]) + (1-max(SCL[SCL != 1]))/10)
if (!is.na(setpar_STRN["dbeta_mu"])) {
mu <- invlogit(setpar_STRN["dbeta_mu"])
v <- abs(setpar_STRN["dbeta_v"])
shape1 <- mu * (((mu * (1 - mu))/v) - 1)
shape2 <- shape1 * (1 - mu)/mu
model_beta <- dbeta(SCL, shape1=shape1, shape2=shape2)
} else {
shape1 <- abs(setpar_STRN["dbeta_shape1"])
shape2 <- abs(setpar_STRN["dbeta_shape2"])
model_beta <- dbeta(SCL, shape1=shape1, shape2=shape2)
}
TSP.TimeWeighted.STRNWeighted.temperature <- sum(TC * dt * STRN_model * model_beta)/sum(dt * STRN_model * model_beta)
TSP.GrowthWeighted.STRNWeighted.temperature <- sum(TC * dSCL * STRN_model * model_beta)/sum(dSCL * STRN_model * model_beta)
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- sum(TC * dt * rT * STRN_model * model_beta)/sum(dt * rT * STRN_model * model_beta)
} else {
TSP.TimeWeighted.STRNWeighted.temperature <- NA
TSP.GrowthWeighted.STRNWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- NA
}
} else {
TSP.PM.GrowthWeighted <- NA
TSP.PM.TimeWeighted.GrowthRateWeighted <- NA
TSP.PM.TimeWeighted <- NA
TSP.TimeWeighted.temperature <- NA
TSP.GrowthWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.temperature <- NA
TSP.TimeWeighted.STRNWeighted.temperature <- NA
TSP.GrowthWeighted.STRNWeighted.temperature <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature <- NA
TSP.length <- NA
}
if (!is.na(indice.begin.middlethird) & !is.na(indice.end.middlethird)) {
TC <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "TempC"])
dt <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "DeltaT"])
dSCL <- unname(diff(df[indice.begin.middlethird:(indice.end.middlethird), "SCL"]))
rT <- unname(df[indice.begin.middlethird:(indice.end.middlethird-1), "R"])
MiddleThird.TimeWeighted.temperature <- sum(TC*dt)/sum(dt)
MiddleThird.GrowthWeighted.temperature <- sum(TC * dSCL)/sum(dSCL)
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature <- sum(TC * rT * dt)/sum(rT * dt)
MiddleThird.length <- unname(time.end.middlethird - time.begin.middlethird)
MiddleThird.begin <- unname(time.begin.middlethird)
MiddleThird.end <- unname(time.end.middlethird)
} else {
MiddleThird.length <- NA
MiddleThird.begin <- NA
MiddleThird.end <- NA
MiddleThird.TimeWeighted.temperature <- NA
MiddleThird.GrowthWeighted.temperature <- NA
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature <- NA
}
if (!is.null(tsd)) {
tsdec <- tsd
tsdec$par <- setpar[names(tsd$par)]
tsdec$fixed.parameters <- setpar[names(tsd$fixed.parameters)]
PM <- predict(tsdec, temperatures = TC, replicate.CI = NULL)
TSP.PM.GrowthWeighted <- sum(PM["50%", ]*dSCL)/sum(dSCL)
TSP.PM.TimeWeighted.GrowthRateWeighted <- sum(PM["50%", ]* rT * dt)/sum(rT * dt)
TSP.PM.TimeWeighted <- sum(PM["50%", ]* dt)/sum(dt)
TSP.TimeWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.GrowthWeighted.sexratio <- predict(tsd, temperatures = TSP.GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.TimeWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.GrowthWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = MiddleThird.TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TimeWeighted.sexratio <- predict(tsd, temperatures = TimeWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
GrowthWeighted.sexratio <- predict(tsd, temperatures = GrowthWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TimeWeighted.GrowthRateWeighted.sexratio <- predict(tsd, temperatures = TimeWeighted.GrowthRateWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
if (!is.na(TSP.TimeWeighted.STRNWeighted.temperature)) {
TSP.TimeWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.GrowthWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.GrowthWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- predict(tsd, temperatures = TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature, replicate.CI = NULL)["50%", 1, drop=TRUE]
} else {
TSP.TimeWeighted.STRNWeighted.sexratio <- NA
TSP.GrowthWeighted.STRNWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- NA
}
} else {
TSP.PM.GrowthWeighted <- NA
TSP.PM.TimeWeighted.GrowthRateWeighted <- NA
TSP.PM.TimeWeighted <- NA
MiddleThird.TimeWeighted.sexratio <- NA
MiddleThird.GrowthWeighted.sexratio <- NA
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio <- NA
TimeWeighted.sexratio <- NA
GrowthWeighted.sexratio <- NA
TimeWeighted.GrowthRateWeighted.sexratio <- NA
TSP.TimeWeighted.sexratio <- NA
TSP.GrowthWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.sexratio <- NA
TSP.TimeWeighted.STRNWeighted.sexratio <- NA
TSP.GrowthWeighted.STRNWeighted.sexratio <- NA
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio <- NA
}
summary <- data.frame(series=serie_en_cours,
Temperature.max=Temperature.max,
TimeWeighted.temperature=TimeWeighted.temperature,
GrowthWeighted.temperature=GrowthWeighted.temperature,
TimeWeighted.GrowthRateWeighted.temperature=TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.temperature=TSP.TimeWeighted.temperature,
TSP.GrowthWeighted.temperature=TSP.GrowthWeighted.temperature,
TSP.TimeWeighted.GrowthRateWeighted.temperature=TSP.TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.STRNWeighted.temperature=TSP.TimeWeighted.STRNWeighted.temperature,
TSP.GrowthWeighted.STRNWeighted.temperature=TSP.GrowthWeighted.STRNWeighted.temperature,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature=TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature,
TSP.length=TSP.length,
TSP.begin=TSP.begin,
TSP.end=TSP.end,
TSP.PM.GrowthWeighted=TSP.PM.GrowthWeighted,
TSP.PM.TimeWeighted.GrowthRateWeighted=TSP.PM.TimeWeighted.GrowthRateWeighted,
TSP.PM.TimeWeighted=TSP.PM.TimeWeighted,
Incubation.length=Incubation.length,
MiddleThird.length=MiddleThird.length,
MiddleThird.begin=MiddleThird.begin,
MiddleThird.end=MiddleThird.end,
MiddleThird.TimeWeighted.temperature=MiddleThird.TimeWeighted.temperature,
MiddleThird.GrowthWeighted.temperature=MiddleThird.GrowthWeighted.temperature,
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature=MiddleThird.TimeWeighted.GrowthRateWeighted.temperature,
TSP.TimeWeighted.sexratio = TSP.TimeWeighted.sexratio,
TSP.GrowthWeighted.sexratio = TSP.GrowthWeighted.sexratio,
TSP.TimeWeighted.GrowthRateWeighted.sexratio = TSP.TimeWeighted.GrowthRateWeighted.sexratio,
TSP.TimeWeighted.STRNWeighted.sexratio = TSP.TimeWeighted.STRNWeighted.sexratio,
TSP.GrowthWeighted.STRNWeighted.sexratio = TSP.GrowthWeighted.STRNWeighted.sexratio,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio = TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio,
MiddleThird.TimeWeighted.sexratio = MiddleThird.TimeWeighted.sexratio,
MiddleThird.GrowthWeighted.sexratio = MiddleThird.GrowthWeighted.sexratio,
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio = MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio,
TimeWeighted.sexratio = TimeWeighted.sexratio,
GrowthWeighted.sexratio = GrowthWeighted.sexratio,
TimeWeighted.GrowthRateWeighted.sexratio = TimeWeighted.GrowthRateWeighted.sexratio)
# je ne garde que les valeurs d'origine
df <- df[df[, "IndiceK"] == 1, ]
indices.df <- data.frame(series=serie_en_cours,
metric.begin.tsp = unname(size.begin.TSP_ec),
metric.end.tsp = unname(size.end.TSP_ec),
hatchling.metric.mean = unname(meanSCL),
hatchling.metric.sd = unname(sdSCL),
time.end.tsp = unname(time.end.TSP),
time.begin.tsp = unname(time.begin.TSP),
metric.end.incubation = unname(df[nrow(df), "SCL"]),
time.begin.middlethird = unname(time.begin.middlethird),
time.end.middlethird = unname(time.end.middlethird),
stop.at.hatchling.metric = unname(attSAT))
} else {
# Sinon, je stocke aussi les indices
# rownames(df) <- 1:nrow(df)
summary <- NULL
indices.df <- data.frame(series=serie_en_cours,
metric.begin.tsp = unname(size.begin.TSP_ec),
metric.end.tsp = unname(size.end.TSP_ec),
hatchling.metric.mean = unname(meanSCL),
hatchling.metric.sd = unname(sdSCL),
time.end.tsp = unname(time.end.TSP),
time.begin.tsp = unname(time.begin.TSP),
metric.end.incubation = unname(df[nrow(df), "SCL"]),
time.begin.middlethird = unname(time.begin.middlethird),
time.end.middlethird = unname(time.end.middlethird),
stop.at.hatchling.metric = unname(attSAT),
indice.begin.middlethird = indice.begin.middlethird,
indice.end.middlethird = indice.end.middlethird,
indice.fin.incubation = indice.fin.incubation,
indice.begin.tsp = indice.begin.tsp,
indice.end.tsp = indice.end.tsp)
}
df[, "DeltaT"] <- c(diff(df[, "Time"]), 0)
} else {
# Je suis en stop.at.hatchling size mais trop court
indices.df <- NULL
summary <- data.frame(series=serie_en_cours,
Temperature.max=max(df[, "TempC"], na.rm = TRUE),
TimeWeighted.temperature=sum(df[, "TempC"]*df[, "DeltaT"])/sum(df[, "DeltaT"]),
GrowthWeighted.temperature=NA,
TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.temperature=NA,
TSP.GrowthWeighted.temperature=NA,
TSP.TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.STRNWeighted.temperature=NA,
TSP.GrowthWeighted.STRNWeighted.temperature=NA,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.temperature=NA,
TSP.length=NA,
TSP.begin=NA,
TSP.end=NA,
TSP.PM.GrowthWeighted=NA,
TSP.PM.TimeWeighted.GrowthRateWeighted=NA,
TSP.PM.TimeWeighted=NA,
Incubation.length=NA,
MiddleThird.length=NA,
MiddleThird.begin=NA,
MiddleThird.end=NA,
MiddleThird.TimeWeighted.temperature=NA,
MiddleThird.GrowthWeighted.temperature=NA,
MiddleThird.TimeWeighted.GrowthRateWeighted.temperature=NA,
TSP.TimeWeighted.sexratio = NA,
TSP.GrowthWeighted.sexratio = NA,
TSP.TimeWeighted.GrowthRateWeighted.sexratio = NA,
TSP.TimeWeighted.STRNWeighted.sexratio = NA,
TSP.GrowthWeighted.STRNWeighted.sexratio = NA,
TSP.TimeWeighted.GrowthRateWeighted.STRNWeighted.sexratio = NA,
MiddleThird.TimeWeighted.sexratio = NA,
MiddleThird.GrowthWeighted.sexratio = NA,
MiddleThird.TimeWeighted.GrowthRateWeighted.sexratio = NA,
TimeWeighted.sexratio = NA,
GrowthWeighted.sexratio = NA,
TimeWeighted.GrowthRateWeighted.sexratio = NA)
df <- NULL
}
if (out == "summary") {
df <- NULL
indices.df <- NULL
}
if (out != "summary") {
summary <- NULL
}
returnlist <- list(x=list(metric=df, summary=summary, series=serie_en_cours, indices.df=indices.df))
names(returnlist) <- serie_en_cours
return(returnlist)
},
clusterEvalQ=list(expr=expression(library("embryogrowth"))),
clusterExport=list(varlist=c("temperatures_ec",
"temperature.heterogeneity",
"hatchling.metric",
"df_random",
"df_random_tsd", "df_random_SexualisationTRN",
"metabolic.heating",
"stop.at.hatchling.metric",
"embryo.stages",
"TSP.borders",
"TSP.begin",
"TSP.end",
"logicTransition",
"integral",
"derivate",
"out",
"metric.end.incubation", "tsd", "SSM"),
envir=environment())
)
if (out=="likelihood") {
return(sum(sapply(AnalyseTraces, FUN = function(x) x$likelihood)))
}
if (out=="metric") {
metric <- lapply(AnalyseTraces, FUN=function(x) x[[1]]$metric)
names(metric) <- sapply(AnalyseTraces, FUN=function(x) x[[1]]$series)
metric <- metric[name]
indices.metric <- lapply(AnalyseTraces, FUN=function(x) x[[1]]$indices.df)
names(indices.metric) <- sapply(AnalyseTraces, FUN=function(x) x[[1]]$series)
indices.metric <- indices.metric[name]
returnlist <- list(dynamic.metric=metric, indices.dynamic.metric=indices.metric)
return(returnlist)
}
# Je crée le summary
if (out == "summary") {
summary_df <- as.data.frame(t(sapply(AnalyseTraces, FUN = function(x) x[[1]]$summary)))
summary.mean <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, mean(unlist(x), na.rm=TRUE)))
colnames(summary.mean) <- c("Series", paste0(colnames(summary.mean)[-1], ".mean"))
summary.se <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, sd(unlist(x), na.rm=TRUE)))
colnames(summary.se) <- c("Series", paste0(colnames(summary.se)[-1], ".se"))
summary <- merge(x=summary.mean, y=summary.se, by.x="Series", by.y = "Series")
if (!is.null(probs)) {
for (pr in probs) {
summary.quantile <- aggregate(summary_df[, -1], by=list(unlist(summary_df$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, quantile(unlist(x), probs = pr, na.rm=TRUE)))
colnames(summary.quantile) <- c("Series", paste0(colnames(summary.quantile)[-1], paste0(".quantile_", as.character(pr))))
summary <- merge(x=summary, y=summary.quantile, by.x="Series", by.y = "Series")
}
}
rownames(summary) <- summary$Series
summary <- summary[name, ]
return(list(summary=summary))
}
# Là je suis en dynamic metric ou dynamic seul
# Je crée le dynamic.metric
if (!stop.at.hatchling.metric) {
df_list <- lapply(AnalyseTraces, FUN = function(x) x[[1]]$metric)
series_list <- sapply(AnalyseTraces, FUN = function(x) x[[1]]$series)
names(df_list) <- series_list
replicate.CI <- ifelse(replicate.CI == 0, 1, replicate.CI)
dynamic.metric <-list()
for (serie_en_cours in unique(series_list)) {
df_list_ec <- df_list[series_list == serie_en_cours]
# J'ai la série en cours dans df_list_ec
# Quel est l'intérêt de prendre df_list ici ? Il y a toutes les séries !
# Si j'ai plusieurs réplicats, ça me fait prendre toutes les séries
t_ec_count <- cbind(df_list_ec[[1]][, "Time", drop=FALSE], count=0)
for (s in seq_along(df_list_ec)) {
common_time <- t_ec_count$Time %in% df_list_ec[[s]][, "Time"]
t_ec_count[common_time, "count"] <- t_ec_count[common_time, "count"] + 1
}
t_ec <- unique(t_ec_count[t_ec_count[, "count"] == length(df_list_ec), 1])
SCL_ec <- sapply(df_list_ec, FUN = function(x) {
x[(x$Time %in% t_ec) & ((x$DeltaT != 0) | (is.na(x$R))), "SCL"]
})
SCL_ec <- as.data.frame(SCL_ec)
# Dans metric_ec j'ai les quantiles de tailles aux temps communs pour tous les réplicats
metric_ec <- t(apply(SCL_ec, MARGIN = 1, FUN = function(x) quantile(x, probs = probs)))
colnames(metric_ec) <- paste0("Metric_", colnames(metric_ec))
# Je sors les R pour tous les temps communs
R_ec <- sapply(df_list_ec, FUN = function(x) {
x[(x$Time %in% t_ec) & ((x$DeltaT != 0) | (is.na(x$R))), "R"]
})
R_ec <- as.data.frame(R_ec)
Rx_ec <- t(apply(R_ec, MARGIN = 1, FUN = function(x) quantile(x, probs = probs, na.rm = TRUE)))
colnames(Rx_ec) <- paste0("R_", colnames(Rx_ec))
# if (length(t_ec) != nrow(metric_ec)) print(serie_en_cours)
df_TC <- df_list_ec[[1]]
df_TC <- df_TC[(df_TC$Time %in% t_ec) & ((df_TC$DeltaT != 0) | (is.na(df_TC$R))), ]
metric_ec <- list(cbind(Time=t_ec, DeltaT=c(0, diff(t_ec)), TempC=df_TC[, "TempC"], metric_ec, Rx_ec))
names(metric_ec) <- serie_en_cours
dynamic.metric <- c(dynamic.metric, metric_ec)
}
summary_dynamic.metric <- as.data.frame(t(sapply(AnalyseTraces, FUN = function(x) x[[1]]$indices.df)))
summary_dynamic.metric.mean <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, mean(unlist(x), na.rm=TRUE)))
colnames(summary_dynamic.metric.mean) <- c("Series", paste0(colnames(summary_dynamic.metric.mean)[-1], ".mean"))
summary_dynamic.metric.se <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, sd(unlist(x), na.rm=TRUE)))
colnames(summary_dynamic.metric.se) <- c("Series", paste0(colnames(summary_dynamic.metric.se)[-1], ".se"))
sdm <- merge(x=summary_dynamic.metric.mean, y=summary_dynamic.metric.se, by.x="Series", by.y = "Series")
if (!is.null(probs)) {
for (pr in probs) {
summary_dynamic.metric.quantile <- aggregate(summary_dynamic.metric[, -1], by=list(unlist(summary_dynamic.metric$series)), FUN=function(x) ifelse(all(is.na(unlist(x))), NA, quantile(unlist(x), probs = pr, na.rm=TRUE)))
colnames(summary_dynamic.metric.quantile) <- c("Series", paste0(colnames(summary_dynamic.metric.quantile)[-1], paste0(".quantile_", as.character(pr))))
sdm <- merge(x=sdm, y=summary_dynamic.metric.quantile, by.x="Series", by.y = "Series")
}
}
} else {
sdm <- NULL
dynamic.metric <- NULL
}
dynamic.metric <- dynamic.metric[name]
rownames(sdm) <- sdm$Series
sdm <- sdm[name, ]
return(list(dynamic.metric=dynamic.metric, summary.dynamic.metric=sdm))
}
Try the embryogrowth package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.