Nothing
R/tsd.R
In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth
Defines functions
tsd
Documented in
tsd
#' tsd estimates the parameters that best describe temperature-dependent sex determination
#' @title Estimate the parameters that best describe temperature-dependent sex determination
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return A list the pivotal temperature, transitional range of temperatures and their SE
#' @param males A vector with male numbers
#' @param females A vector with female numbers
#' @param N A vector with total numbers
#' @param temperatures The constant incubation temperatures used to fit sex ratio
#' @param durations The duration of incubation or TSP used to fit sex ratio
#' @param df A dataframe with at least two columns named males, females or N and temperatures, Incubation.temperature or durations column
#' @param l Sex ratio limits to define TRT are l and 1-l (see Girondot, 1999)
#' @param parameters.initial Initial values for P, S or K search as a vector, ex. c(P=29, S=-0.3)
#' @param fixed.parameters Parameters that will not be changed
#' @param males.freq If TRUE data are shown as males frequency
#' @param equation Can be "logistic", "Hill", "A-logistic", "Hulin", "Double-A-logistic", "flexit", "flexit*", "GSD", "logit", "probit"
#' @param replicate.CI Number of replicates to estimate confidence intervals
#' @param replicate.NullDeviance Number of replicates to estimate null distribution of deviance
#' @param SE If FALSE, does not estimate SE of parameters. Can be use when something wrong happens.
#' @param range.CI The range of confidence interval for estimation, default=0.95
#' @param print Should the results be printed at screen? TRUE (default) or FALSE
#' @param control List of parameters used in optim.
#' @param method method used for optim. Can be "BFGS", the most rapid or "Nelder-Mead" for special cases using n parameter.
#' @description Estimate the parameters that best describe the thermal reaction norm for sex ratio when temperature-dependent sex determination occurs.\cr
#' It can be used also to evaluate the relationship between incubation duration and sex ratio.\cr
#' The parameter l was defined in Girondot (1999). The TRT is defined from the difference between the two boundary temperatures giving sex ratios of \eqn{l} and \eqn{1 - l}, respectively:\cr
#' For logistic model (Girondot, 1999), it follows \deqn{TRT_{l}=abs\left ( S\: K_{l} \right )}{TRTl = abs( S.Kl )}
#' where \eqn{K_{l}}{Kl} is a constant equal to \eqn{2\: log\left ( \frac{l}{1-l} \right )}{2 log(l / ( 1 - l))}.\cr
#' In Girondot (1999), l was 0.05 and then the TRT was defined as being the range of temperatures producing from 5\% to 95\% of each sex.\cr
#' The default model is named logistic. This model (as well as the logit one) has the particularity to have a symmetric shape around P.\cr
#' The \emph{logistic} model is:
#' \deqn{SR(T) = 1 / (1 + exp((1 / S) * (P - T))))}{SR(T) = 1 / (1 + exp((1 / S) * (P - T))))}
#' The \emph{logit} model is:
#' \deqn{SR(T) = 1 / (1 + exp(4 * S * (P - T))))}{SR(T) = 1 / (1 + exp(4 * S * (P - T))))}
#' The other models have been built to alleviate this constraint. Hill and A-logistic models can be asymmetric, but it is impossible to control independently the low and high transitions.\cr
#' Hulin model is assymmetric but the control of asymmetry is difficult to manage.\cr
#' If asymmetric model is selected, it is always better to use \emph{flexit} model.
#' \deqn{if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}}{if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}}
#' \deqn{if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}}{if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}}
#' with:\cr
#' \deqn{S_1 = S/((4/K_1)*(2^{(-K_1)})^{(1/K_1+1)}*(2^{K_1}-1))}{S_1 = S/((4/K_1)*(2^{(-K_1)})^(1/K_1+1)*(2^{K_1}-1))}
#' \deqn{S_2 = S/((4/K_2)*(2^{(-K_2)})^{(1/K_2+1)}*(2^{K_2}-1))}{S_2 = S/((4/K_2)*(2^{(-K_2)})^(1/K_2+1)*(2^{K_2}-1))}
#' The \emph{flexit*} model is defined as (QBT is the Quasi-Binary Threshold):
#' \deqn{QBT = 1/ (1 + exp(100 * (P - T)))}{QBT = 1 / (1 + exp(100 * (P - T)))}
#' \deqn{SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))}{SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))}
#' @references
#' \insertRef{1515}{embryogrowth}\cr
#' \insertRef{3534}{embryogrowth}\cr
#' \insertRef{11754}{embryogrowth}\cr
#' \insertRef{5790}{embryogrowth}\cr
#' @family Functions for temperature-dependent sex determination
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Caretta caretta" & (!is.na(Sexed) & Sexed!=0))
#' tsdL <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="logistic", replicate.CI=NULL))
#' tsdH <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Hill", replicate.CI=NULL))
#' tsdR <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="A-logistic", replicate.CI=NULL))
#' tsdF <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Flexit", replicate.CI=NULL))
#' tsdF1 <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Flexit*", replicate.CI=NULL))
#' tsdDR <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Double-A-logistic", replicate.CI=NULL))
#' gsd <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="GSD", replicate.CI=NULL))
#' compare_AIC(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' flexit=tsdF,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd)
#' compare_AICc(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd, factor.value = -1)
#' compare_BIC(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd, factor.value = -1)
#' ##############
#' eo <- subset(DatabaseTSD, Species=="Emys orbicularis", c("Males", "Females",
#' "Incubation.temperature"))
#'
#' eo_Hill <- with(eo, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Hill"))
#' eo_Hill <- tsd(df=eo, equation="Hill", replicate.CI=NULL)
#' eo_logistic <- tsd(eo, replicate.CI=NULL)
#' eo_Alogistic <- with(eo, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="a-logistic", replicate.CI=NULL))
#' ### The Hulin model is a modification of A-logistic (See Hulin et al. 2009)
#'
#' ########## Caution
#' ### It should not be used anymore as it can produce unexpected results
#' par <- eo_Alogistic$par
#' names(par)[which(names(par)=="K")] <- "K2"
#' par <- c(par, K1=0)
#' eo_Hulin <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="Hulin", replicate.CI=NULL))
#'
#' ### The Double-A-logistic model is a A-logistic model with K1 and K2 using respectively
#' ### below and above P
#'
#' ########## Caution
#' ### The curve is not smooth at pivotal temperature
#'
#' par <- eo_Alogistic$par
#' names(par)[which(names(par)=="K")] <- "K2"
#' par <- c(par, K1=as.numeric(par["K2"])*0.8)
#' par["K2"] <- par["K2"]*0.8
#' eo_Double_Alogistic <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="Double-a-logistic", replicate.CI=NULL))
#'
#' ### The flexit model is modeled with K1 and K2 using respectively
#' ### below and above P and smooth transition at P; S is the slope at P
#'
#' par <- c(eo_logistic$par["P"], 1/4*eo_logistic$par["S"], K1=1, K2=1)
#' eo_flexit <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="flexit", replicate.CI=NULL))
#'
#' compare_AIC(Logistic=eo_logistic, Hill=eo_Hill, Alogistic=eo_Alogistic,
#' Hulin=eo_Hulin, Double_Alogistic=eo_Double_Alogistic,
#' flexit=eo_flexit)
#' ## Note that SE for lower limit of TRT is wrong
#' plot(eo_flexit)
#' ## To get correct confidence interval, check \code{tsd_MHmcmc()}.
#'
#' ### Note the asymmetry of the Double-A-logistic and flexit models
#' predict(eo_Double_Alogistic,
#' temperatures=c(eo_Double_Alogistic$par["P"]-0.2, eo_Double_Alogistic$par["P"]+0.2))
#' predict(eo_Double_Alogistic)
#'
#' (p <- predict(eo_flexit,
#' temperatures=c(eo_flexit$par["P"]-0.3, eo_flexit$par["P"]+0.3)))
#' p["50%", 1]-0.5; 0.5-p["50%", 2]
#' predict(eo_flexit)
#'
#' ### It can be used also for incubation duration
#' CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Caretta caretta" & Sexed!=0)
#' tsdL_IP <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' durations=IP.mean,
#' equation="logistic", replicate.CI=NULL))
#' plot(tsdL_IP, xlab="Incubation durations in days")
#' # Example with Chelonia mydas
#' cm <- subset(DatabaseTSD, Species=="Chelonia mydas" & !is.na(Sexed), c("Males", "Females",
#' "Incubation.temperature", "RMU.2010"))
#' tsd(subset(cm, subset=RMU.2010=="Pacific, SW"))
#' tsd(subset(cm, subset=RMU.2010=="Pacific, Northwest"))
#' tsd(subset(cm, subset=RMU.2010=="Atlantic, S Caribbean"))
#'
#' ### Eretmochelys imbricata
#' Ei_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" &
#' Species=="Eretmochelys imbricata")
#' Ei_AtlanticW <- subset(DatabaseTSD, RMU.2010=="Atlantic, W (Caribbean and E USA)" &
#' Species=="Eretmochelys imbricata")
#' Ei_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Eretmochelys imbricata")
#' Ei_PacSW <- tsd(Ei_PacificSW)
#' Ei_AtlW <- tsd(Ei_AtlanticW)
#' Ei_AtlSW <- tsd(Ei_AtlanticSW)
#'
#' plot(Ei_PacSW, xlim=c(27, 33), show.PTRT = FALSE, main=expression(italic("Eretmochelys imbricata")))
#' par(new=TRUE)
#' plot(Ei_AtlW, xlim=c(27, 33), col="red", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="red")
#' par(new=TRUE)
#' plot(Ei_AtlSW, xlim=c(27, 33), col="blue", xlab="", ylab="", axes=FALSE,
#' xaxt="n", show.PTRT = FALSE, errbar.col="blue")
#' legend("topright", legend=c("Pacific, SW", "Atlantic, W", "Atlantic, SW"), lty=1,
#' col=c("black", "red", "blue"))
#'
#' ### Chelonia mydas
#' Cm_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_PacificNW <- subset(DatabaseTSD, RMU.2010=="Pacific, NW" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_AtlanticSC <- subset(DatabaseTSD, RMU.2010=="Atlantic, S Caribbean" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_IndianSE <- subset(DatabaseTSD, RMU.2010=="Indian, SE" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_PacSW <- tsd(Cm_PacificSW)
#' Cm_PacNW <- tsd(Cm_PacificNW)
#' Cm_IndSE <- tsd(Cm_IndianSE)
#' Cm_AtlSC <- tsd(Cm_AtlanticSC)
#'
#' plot(Cm_PacSW, xlim=c(24, 34), show.PTRT = FALSE, main=expression(italic("Chelonia mydas")))
#' par(new=TRUE)
#' plot(Cm_PacNW, xlim=c(24, 34), col="red", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="red")
#' par(new=TRUE)
#' plot(Cm_IndSE, xlim=c(24, 34), col="blue", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="blue")
#' par(new=TRUE)
#' plot(Cm_AtlSC, xlim=c(24, 34), col="green", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="green")
#'
#' # To fit a TSDII or FMF TSD pattern, you must indicate P_low, S_low, P_high, and S_high
#' # for logistic model and P_low, S_low, K1_low, K2_low, P_high, S_high, K1_high, and K2_high for
#' # flexit model
#' # The model must be 0-1 for low and 1-0 for high with P_low < P_high
#'
#' Chelydra_serpentina <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Chelydra serpentina")
#'
#' model_TSDII <- tsd(Chelydra_serpentina, males.freq=FALSE,
#' parameters.initial=c(P_low=21, S_low=0.3, P_high=28, S_high=-0.4),
#' equation="logistic")
#' plot(model_TSDII, lab.TRT = "TRT l = 5 %")
#' priors <- tsd_MHmcmc_p(result=model_TSDII, accept=TRUE)
#' out_mcmc <- tsd_MHmcmc(result=model_TSDII, n.iter=10000, parametersMCMC=priors)
#' plot(model_TSDII, resultmcmc=out_mcmc, lab.TRT = "TRT l = 5 %")
#' predict(model_TSDII, temperatures=25:35)
#'
#' # Podocnemis expansa
#' Podocnemis_expansa <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Podocnemis expansa")
#' Podocnemis_expansa_Valenzuela_2001 <- subset(Podocnemis_expansa,
#' Reference=="Valenzuela, 2001")
#' PeL2001 <- tsd(df=Podocnemis_expansa_Valenzuela_2001)
#' # The pivotal temperature is 32.133 °C (CI 95% 31.495;32.766)
#' # In Valenzuela, 2001: "Using data from the present study alone,
#' # the critical temperature was 32.2 °C by both methods and the 95%
#' # confidence limits were 31.4 °C and 32.9 °C."
#' # Data are close but not identical to what was published.
#'
#' # The pivotal temperature calculated by maximum likelihood and by inverse
#' # prediction from logistic regression, was 32.6°C using raw data from
#' # 1991 (N. Valenzuela, unpublished data) and from this study. The lower
#' # and upper 95% confidence limits of the pivotal temperature were 32.2°C
#' # and 33.2°C,
#'
#' Podocnemis_expansa_Valenzuela_1997 <- subset(Podocnemis_expansa,
#' subset=(((Reference=="Lance et al., 1992; Valenzuela et al., 1997") |
#' (Reference=="Valenzuela, 2001")) &
#' (!is.na(Sexed)) & (Sexed != 0)))
#'
#' PeL1997 <- tsd(df=Podocnemis_expansa_Valenzuela_1997)
#'
#' # Gekko japonicus
#'
#' Gekko_japonicus <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Gekko japonicus")
#' model_TSDII_gj <- tsd(Gekko_japonicus, males.freq=TRUE,
#' parameters.initial=c(P_low=26, S_low=1.5,
#' P_high=31, S_high=-1.5),
#' equation="logistic")
#' plot(model_TSDII_gj, lab.TRT = "TRT l = 5 %")
#' print(model_TSDII_gj)
#' prior <- tsd_MHmcmc_p(result = model_TSDII_gj, accept = TRUE)
#' prior <- structure(list(
#' Density = c("dnorm", "dnorm", "dnorm", "dnorm"),
#' Prior1 = c(26, 0.3, 31, -0.4),
#' Prior2 = c(2, 1, 2, 1),
#' SDProp = c(2, 0.5, 2, 0.5),
#' Min = c(25, -2, 25, -2),
#' Max = c(35, 2, 35, 2),
#' Init = c(26, 0.3, 31, -0.4)),
#' row.names = c("P_low", "S_low", "P_high", "S_high"),
#' class = "data.frame")
#'
#' result_mcmc_tsd_gj <- tsd_MHmcmc(result=model_TSDII_gj,
#' parametersMCMC=prior, n.iter=10000, n.chains = 1,
#' n.adapt = 0, thin=1, trace=FALSE, adaptive=TRUE)
#' summary(result_mcmc_tsd_gj)
#' plot(result_mcmc_tsd_gj, parameters="P_low", scale.prior=TRUE, xlim=c(20, 30), las=1)
#' plot(result_mcmc_tsd_gj, parameters="P_high", scale.prior=TRUE, xlim=c(25, 35), las=1)
#' plot(model_TSDII_gj, resultmcmc = result_mcmc_tsd_gj)
#'
#' # Trachemys scripta elegans
#' # Take care, the pattern reflects large population variation
#'
#' Tse <- subset(DatabaseTSD, Species=="Trachemys scripta" & Subspecies == "elegans" & !is.na(Sexed))
#' Tse_logistic <- tsd(Tse)
#' plot(Tse_flexit)
#' compare_AICc(logistic=Tse_logistic, flexit=Tse_flexit)
#' plot(Tse_flexit)
#'
#' # Exemple when only proportion is known; experimental
#' Ei_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" &
#' Species=="Eretmochelys imbricata")
#' males <- Ei_PacificSW$Males/(Ei_PacificSW$Males+Ei_PacificSW$Females)*100
#' females <- 100-(Ei_PacificSW$Males/(Ei_PacificSW$Males+Ei_PacificSW$Females)*100)
#' temperatures <- Ei_PacificSW$Incubation.temperature
#' Ei_PacSW <- tsd(Ei_PacificSW)
#' par <- c(Ei_PacSW$par, n=10)
#' embryogrowth:::.tsd_fit(par=par, males=males, N=males+females, temperatures=temperatures,
#' equation="logistic")
#' Ei_PacSW_NormalApproximation <- tsd(males=males, females=females,
#' temperatures=temperatures,
#' parameters.initial=par)
#' Ei_PacSW_NormalApproximation$par
#' Ei_PacSW$par
#' # The data looks like only n=0.01 observations were done
#' # This is the reason of the large observed heterogeneity
#' plot(Ei_PacSW_NormalApproximation)
#' }
#' @export
tsd <- function(df=NULL ,
males=NULL ,
females=NULL ,
N=NULL ,
temperatures=NULL ,
durations=NULL ,
l=0.05 ,
parameters.initial=c(P=30, S=-2, K=0, K1=1, K2=0,
SL = -1, SH = -1) ,
males.freq=TRUE ,
fixed.parameters=NULL ,
equation="logistic" ,
replicate.CI=10000 ,
range.CI=0.95 ,
SE=TRUE ,
replicate.NullDeviance=1000 ,
control=list(maxit=1000) ,
print=TRUE ,
method="BFGS" ) {
# df=NULL; males=NULL; females=NULL; N=NULL; temperatures=NULL; durations=NULL; l=0.05; parameters.initial=c(P=NA, S=-0.5, K=0, K1=1, K2=0, SL = -1, SH = -1); males.freq=TRUE; fixed.parameters=NULL; equation="logistic"; replicate.CI=10000; range.CI=0.95; SE=TRUE;print=TRUE; control=list(maxit=1000); replicate.NullDeviance=1000; method="BFGS"
# CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" & Species=="Caretta caretta" & Sexed!=0)
# males=CC_AtlanticSW$Males; females=CC_AtlanticSW$Females; temperatures=CC_AtlanticSW$Incubation.temperature; equation="logistic"
# parameters.initial = c(P=29.19); fixed.parameters = c(S=-0.21)
equation <- tolower(equation)
if (is.null(replicate.CI)) replicate.CI <- 0
equation <- match.arg(equation,
choices = c("logistic", "hill", "a-logistic", "hulin",
"double-a-logistic", "flexit", "flexit*", "gsd",
"probit", "logit"),
several.ok = FALSE)
# method <- tolower(method)
method <- match.arg(method, choices=c("BFGS", "Nelder-Mead"),
several.ok = FALSE)
if (!is.null(df)) {
if ((!inherits(df, "data.frame")) & (!inherits(df, "matrix"))) {
stop("df parameter must be a data.frame or a matrix")
}
namesdf <- tolower(colnames(df))
males <- NULL
females <- NULL
N <- NULL
if (any(namesdf=="males")) males <- df[, which(namesdf=="males")]
if (any(namesdf=="male")) males <- df[, which(namesdf=="male")]
if (any(namesdf=="females")) females <- df[, which(namesdf=="females")]
if (any(namesdf=="female")) females <- df[, which(namesdf=="female")]
if (any(namesdf=="n")) N <- df[, which(namesdf=="n")]
if (any(namesdf=="sexed")) N <- df[, which(namesdf=="sexed")]
if (any(namesdf=="temperatures")) temperatures <- df[, which(namesdf=="temperatures")]
if (any(namesdf=="temperature")) temperatures <- df[, which(namesdf=="temperature")]
if (any(namesdf=="incubation.temperature.set")) temperatures <- df[, which(namesdf=="incubation.temperature.set")]
if (any(namesdf=="durations")) durations <- df[, which(namesdf=="durations")]
if (any(namesdf=="duration")) durations <- df[, which(namesdf=="duration")]
if (any(namesdf=="IP.mean")) durations <- df[, which(namesdf=="IP.mean")]
}
if (is.null(durations) & is.null(temperatures)) {
stop("Or temperatures or durations must be provided")
}
if (!is.null(durations) & !is.null(temperatures)) {
stop("Temperatures and durations cannot be both provided")
}
# si je retire des lignes, c'est decale
# if (!is.null(males)) males <- na.omit(males)
# if (!is.null(females)) females <- na.omit(females)
# if (!is.null(N)) N <- na.omit(N)
if (is.null(temperatures)) {
IP <- TRUE
temperatures <- durations
} else {
IP <- FALSE
}
cpt1 <- ifelse(is.null(males), 0, 1)+ifelse(is.null(females), 0, 1)+ifelse(is.null(N), 0, 1)
if (cpt1<2) {
stop("At least two informations from males, females or N must be provided")
}
if (is.null(males)) males <- N-females
if (is.null(females)) females <- N-males
if (is.null(N)) N <- males+females
if (length(temperatures) != length(N)) {
stop("A temperature or duration must be provided for each experiment")
}
if (all(males * females == 0)) {
warning("At least one incubation temperature should produce mixed sex ratio. Use MCMC rather.")
}
if (equation == "gsd") {
value <- -sum(dbinom(x=males, size=males+females, prob=0.5, log=TRUE))
result <- list(par = NULL, SE=NULL, hessian=NULL,
TRT=NULL,
SE_TRT=NULL,
message=NULL,
convergence=0,
counts=NULL,
value = value,
AIC= 2*value + 2*0,
AICc= 2*value + (2*0*(0+1))/(sum(N)-0-1),
BIC= 2*value + 0
)
# limit.low.TRT <- min(temperatures)
# limit.high.TRT <- max(temperatures)
} else {
ppi <- c(parameters.initial, fixed.parameters)
if (is.na(ppi["P_low"])) {
if (is.na(ppi["P"])) {
if ((equation!="probit") & (equation!="logit")) {
ppi["P"] <- temperatures[which.min(abs((males/(males+females))-0.5))]
} else {
ppi["P"] <- 0
}
}
if (is.na(ppi["S"])) {
if ((equation!="probit")) {
pente <- lm(males/(males+females) ~ temperatures)
ppi["S"] <- pente$coefficients["temperatures"]
if ((equation != "flexit") & (equation != "flexit*")) ppi["S"] <- 1/(4*ppi["S"])
} else {
ppi["S"] <- 0
}
}
}
# if (IP) ppi["S"] <- abs(ppi["S"])
if (equation != "a-logistic") ppi <- ppi[which(names(ppi) != "K")]
if ((equation != "hulin") & (equation != "double-a-logistic") & (equation != "flexit")) {
ppi <- ppi[which(names(ppi) != "K1")]
ppi <- ppi[which(names(ppi) != "K2")]
}
if (equation != "flexit*") {
ppi <- ppi[which(names(ppi) != "SL")]
ppi <- ppi[which(names(ppi) != "SH")]
} else {
ppi <- ppi[which(names(ppi) != "S")]
}
if (!is.null(fixed.parameters))
ppi[names(fixed.parameters)] <- unname(fixed.parameters)
repeat {
# result <- optim(par, embryogrowth:::.tsd_fit, fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures, equation=equation, method="BFGS", hessian=TRUE, control = list(maxit=1000))
result <- try(optim(ppi, getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures,
equation=equation, method=method, hessian=FALSE, control = control), silent=TRUE)
if (inherits(result, "try-error")) {
stop("Try using method='Nelder-Mead' or use MCMC but it will be complicated ! Contact me.")
}
if (result$convergence != 1) break
ppi <- result$par
if (print) print("Convergence is not acheived. Optimization continues !")
}
result_hessian <- try(optim(result$par, getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures,
equation=equation, method=method, hessian=TRUE, control = control), silent=TRUE)
ppi <- c(result$par, fixed.parameters)
if (inherits(result_hessian, "try-error")) {
warning("Likelihood is flat close to ML point and Hessian cannot be estimated")
result$hessian <- NULL
result$SE <- NULL
replicate.CI = 0
} else {
result$hessian <- result_hessian$hessian
if (SE) {
rh <- SEfromHessian(result$hessian, hessian=TRUE)
result$SE <- rh$SE
}
}
# result$hessian <- result$hessian
result$AIC <- 2*result$value+2*length(result$par)
# Correction le 21/10/2020
result$AICc <- result$AIC +(2*length(result$par)*(length(result$par)+1))/(sum(N)-length(result$par)-1)
# Correction le 21/10/2020
result$BIC <- 2*result$value+ length(result$par)*log(sum(N))
}
result$range.CI <- range.CI
result$l <- l
result$equation <- equation
result$males <- males
result$females <- females
result$N <- N
result$temperatures <- temperatures
result$males.freq <- males.freq
result$fixed.parameters <- fixed.parameters
result$type <- ifelse(IP, "duration", "temperature")
if (all(names(ppi) != "n")) {
result$deviance <- -2*(-result$value - sum(dbinom(x = males, size=N,
prob = males/N, log = TRUE)))
} else {
pt <- males/N
pt <- ifelse(pt<=1E-9, 1E-9, pt)
pt <- ifelse(pt>=1-1E-9, 1-1E-9, pt)
sd <- sqrt((pt)*(1-(pt))/ppi["n"])
pr <- dnorm(males/N, mean=pt, sd=sd, log=TRUE)
result$deviance <- -2*(-result$value - sum(pr))
}
# degrees of freedom calculated from the difference of the number of parameters in the saturated and the fitted model.
result$df <- length(males) - length(result$par)
result$pvalue <- pchisq(q=result$deviance, df=result$df, lower.tail = FALSE)
if (equation != "gsd") {
result <- addS3Class(result, "tsd")
# result <<- result
if (!is.null(result$hessian)) {
vcov <- solve(result$hessian)
if (any(eigen(vcov)$values < 0)) {
warning("Non-positive definite variance-covariance matrix; use MCMC to get credible interval.")
result$hessian <- NULL
result$SE <- NULL
replicate.CI = 0
}
}
o <- P_TRT(x=result, l=l,
replicate.CI = replicate.CI,
probs = c((1-range.CI)/2, 0.5, 1-(1-range.CI)/2))
result$P_TRT <- o$P_TRT_quantiles
}
probtheor <- getFromNamespace(x=".modelTSD", ns="embryogrowth")(c(result$par, result$fixed.parameters),
result$temperatures,
equation=equation)
result$predictMaleFrequency <- probtheor
if (!is.null(replicate.NullDeviance) & (replicate.NullDeviance != 0)) {
pt <- NULL
# Bug correct 2020-01-11: change 1000 to replicate.NullDeviance
for (i in 1:replicate.NullDeviance) {
m <- rbinom(n=length(result$males), size=result$N, prob = probtheor)
result_dev <- optim(par=result$par, fn=getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=result$fixed.parameters, males=m, N=result$N,
temperatures=result$temperatures,
equation=result$equation, method=method,
hessian=FALSE, control = control)
deviance_dev <- -2*(-result$value - sum(dbinom(x = m, size=N,
prob = m/N, log = TRUE)))
pt <- c(pt, deviance_dev)
}
result$NullDeviance <- pt
result$NullDeviancePvalue <- sum(result$deviance < pt)/length(pt)
}
if (equation!="gsd" & print) {
# SI je n'ai qu'un PT c'est que c'est un profile TSDIA ou IB
if (any(colnames(o$P_TRT_quantiles)=="PT")) {
if (!is.null(replicate.CI) & (replicate.CI != 0)) {
print(paste0("The pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[3, "PT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[3, "PT"]))))
print(paste0("The transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT"], o$P_TRT_quantiles[3, "TRT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT"], o$P_TRT_quantiles[3, "TRT"]))))
print(paste0("The lower limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"]))))
print(paste0("The upper limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"]))))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
} else {
print(paste0("The pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT"])))
print(paste0("The transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT"])))
print(paste0("The lower limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT"])))
print(paste0("The upper limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT"])))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
}
} else {
if (!is.null(replicate.CI) & (replicate.CI != 0)) {
print(paste0("The lower pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[3, "PT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[3, "PT_low"]))))
print(paste0("The lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT_low"], o$P_TRT_quantiles[3, "TRT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT_low"], o$P_TRT_quantiles[3, "TRT_low"]))))
print(paste0("The lower limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"]))))
print(paste0("The upper limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"]))))
print(paste0("The upper pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[3, "PT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[3, "PT_high"]))))
print(paste0("The upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT_high"], o$P_TRT_quantiles[3, "TRT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT_high"], o$P_TRT_quantiles[3, "TRT_high"]))))
print(paste0("The lower limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"]))))
print(paste0("The upper limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"]))))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
} else {
print(paste0("The lower pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_low"])))
print(paste0("The lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_low"])))
print(paste0("The lower limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_low"])))
print(paste0("The upper limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_low"])))
print(paste0("The upper pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_high"])))
print(paste0("The upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_high"])))
print(paste0("The lower limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_high"])))
print(paste0("The upper limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_high"])))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
}
}
}
result <- addS3Class(result, "tsd")
return(invisible(result))
}
Try the embryogrowth package in your browser
Any scripts or data that you put into this service are public.
embryogrowth documentation built on Sept. 11, 2024, 8:16 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 tsd
Documented in tsd
#' tsd estimates the parameters that best describe temperature-dependent sex determination
#' @title Estimate the parameters that best describe temperature-dependent sex determination
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return A list the pivotal temperature, transitional range of temperatures and their SE
#' @param males A vector with male numbers
#' @param females A vector with female numbers
#' @param N A vector with total numbers
#' @param temperatures The constant incubation temperatures used to fit sex ratio
#' @param durations The duration of incubation or TSP used to fit sex ratio
#' @param df A dataframe with at least two columns named males, females or N and temperatures, Incubation.temperature or durations column
#' @param l Sex ratio limits to define TRT are l and 1-l (see Girondot, 1999)
#' @param parameters.initial Initial values for P, S or K search as a vector, ex. c(P=29, S=-0.3)
#' @param fixed.parameters Parameters that will not be changed
#' @param males.freq If TRUE data are shown as males frequency
#' @param equation Can be "logistic", "Hill", "A-logistic", "Hulin", "Double-A-logistic", "flexit", "flexit*", "GSD", "logit", "probit"
#' @param replicate.CI Number of replicates to estimate confidence intervals
#' @param replicate.NullDeviance Number of replicates to estimate null distribution of deviance
#' @param SE If FALSE, does not estimate SE of parameters. Can be use when something wrong happens.
#' @param range.CI The range of confidence interval for estimation, default=0.95
#' @param print Should the results be printed at screen? TRUE (default) or FALSE
#' @param control List of parameters used in optim.
#' @param method method used for optim. Can be "BFGS", the most rapid or "Nelder-Mead" for special cases using n parameter.
#' @description Estimate the parameters that best describe the thermal reaction norm for sex ratio when temperature-dependent sex determination occurs.\cr
#' It can be used also to evaluate the relationship between incubation duration and sex ratio.\cr
#' The parameter l was defined in Girondot (1999). The TRT is defined from the difference between the two boundary temperatures giving sex ratios of \eqn{l} and \eqn{1 - l}, respectively:\cr
#' For logistic model (Girondot, 1999), it follows \deqn{TRT_{l}=abs\left ( S\: K_{l} \right )}{TRTl = abs( S.Kl )}
#' where \eqn{K_{l}}{Kl} is a constant equal to \eqn{2\: log\left ( \frac{l}{1-l} \right )}{2 log(l / ( 1 - l))}.\cr
#' In Girondot (1999), l was 0.05 and then the TRT was defined as being the range of temperatures producing from 5\% to 95\% of each sex.\cr
#' The default model is named logistic. This model (as well as the logit one) has the particularity to have a symmetric shape around P.\cr
#' The \emph{logistic} model is:
#' \deqn{SR(T) = 1 / (1 + exp((1 / S) * (P - T))))}{SR(T) = 1 / (1 + exp((1 / S) * (P - T))))}
#' The \emph{logit} model is:
#' \deqn{SR(T) = 1 / (1 + exp(4 * S * (P - T))))}{SR(T) = 1 / (1 + exp(4 * S * (P - T))))}
#' The other models have been built to alleviate this constraint. Hill and A-logistic models can be asymmetric, but it is impossible to control independently the low and high transitions.\cr
#' Hulin model is assymmetric but the control of asymmetry is difficult to manage.\cr
#' If asymmetric model is selected, it is always better to use \emph{flexit} model.
#' \deqn{if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}}{if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}}
#' \deqn{if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}}{if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}}
#' with:\cr
#' \deqn{S_1 = S/((4/K_1)*(2^{(-K_1)})^{(1/K_1+1)}*(2^{K_1}-1))}{S_1 = S/((4/K_1)*(2^{(-K_1)})^(1/K_1+1)*(2^{K_1}-1))}
#' \deqn{S_2 = S/((4/K_2)*(2^{(-K_2)})^{(1/K_2+1)}*(2^{K_2}-1))}{S_2 = S/((4/K_2)*(2^{(-K_2)})^(1/K_2+1)*(2^{K_2}-1))}
#' The \emph{flexit*} model is defined as (QBT is the Quasi-Binary Threshold):
#' \deqn{QBT = 1/ (1 + exp(100 * (P - T)))}{QBT = 1 / (1 + exp(100 * (P - T)))}
#' \deqn{SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))}{SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))}
#' @references
#' \insertRef{1515}{embryogrowth}\cr
#' \insertRef{3534}{embryogrowth}\cr
#' \insertRef{11754}{embryogrowth}\cr
#' \insertRef{5790}{embryogrowth}\cr
#' @family Functions for temperature-dependent sex determination
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Caretta caretta" & (!is.na(Sexed) & Sexed!=0))
#' tsdL <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="logistic", replicate.CI=NULL))
#' tsdH <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Hill", replicate.CI=NULL))
#' tsdR <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="A-logistic", replicate.CI=NULL))
#' tsdF <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Flexit", replicate.CI=NULL))
#' tsdF1 <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Flexit*", replicate.CI=NULL))
#' tsdDR <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Double-A-logistic", replicate.CI=NULL))
#' gsd <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="GSD", replicate.CI=NULL))
#' compare_AIC(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' flexit=tsdF,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd)
#' compare_AICc(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd, factor.value = -1)
#' compare_BIC(Logistic_Model=tsdL, Hill_model=tsdH, Alogistic_model=tsdR,
#' DoubleAlogistic_model=tsdDR, GSD_model=gsd, factor.value = -1)
#' ##############
#' eo <- subset(DatabaseTSD, Species=="Emys orbicularis", c("Males", "Females",
#' "Incubation.temperature"))
#'
#' eo_Hill <- with(eo, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="Hill"))
#' eo_Hill <- tsd(df=eo, equation="Hill", replicate.CI=NULL)
#' eo_logistic <- tsd(eo, replicate.CI=NULL)
#' eo_Alogistic <- with(eo, tsd(males=Males, females=Females,
#' temperatures=Incubation.temperature.set,
#' equation="a-logistic", replicate.CI=NULL))
#' ### The Hulin model is a modification of A-logistic (See Hulin et al. 2009)
#'
#' ########## Caution
#' ### It should not be used anymore as it can produce unexpected results
#' par <- eo_Alogistic$par
#' names(par)[which(names(par)=="K")] <- "K2"
#' par <- c(par, K1=0)
#' eo_Hulin <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="Hulin", replicate.CI=NULL))
#'
#' ### The Double-A-logistic model is a A-logistic model with K1 and K2 using respectively
#' ### below and above P
#'
#' ########## Caution
#' ### The curve is not smooth at pivotal temperature
#'
#' par <- eo_Alogistic$par
#' names(par)[which(names(par)=="K")] <- "K2"
#' par <- c(par, K1=as.numeric(par["K2"])*0.8)
#' par["K2"] <- par["K2"]*0.8
#' eo_Double_Alogistic <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="Double-a-logistic", replicate.CI=NULL))
#'
#' ### The flexit model is modeled with K1 and K2 using respectively
#' ### below and above P and smooth transition at P; S is the slope at P
#'
#' par <- c(eo_logistic$par["P"], 1/4*eo_logistic$par["S"], K1=1, K2=1)
#' eo_flexit <- with(eo, tsd(males=Males, females=Females,
#' parameters.initial=par,
#' temperatures=Incubation.temperature.set,
#' equation="flexit", replicate.CI=NULL))
#'
#' compare_AIC(Logistic=eo_logistic, Hill=eo_Hill, Alogistic=eo_Alogistic,
#' Hulin=eo_Hulin, Double_Alogistic=eo_Double_Alogistic,
#' flexit=eo_flexit)
#' ## Note that SE for lower limit of TRT is wrong
#' plot(eo_flexit)
#' ## To get correct confidence interval, check \code{tsd_MHmcmc()}.
#'
#' ### Note the asymmetry of the Double-A-logistic and flexit models
#' predict(eo_Double_Alogistic,
#' temperatures=c(eo_Double_Alogistic$par["P"]-0.2, eo_Double_Alogistic$par["P"]+0.2))
#' predict(eo_Double_Alogistic)
#'
#' (p <- predict(eo_flexit,
#' temperatures=c(eo_flexit$par["P"]-0.3, eo_flexit$par["P"]+0.3)))
#' p["50%", 1]-0.5; 0.5-p["50%", 2]
#' predict(eo_flexit)
#'
#' ### It can be used also for incubation duration
#' CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Caretta caretta" & Sexed!=0)
#' tsdL_IP <- with (CC_AtlanticSW, tsd(males=Males, females=Females,
#' durations=IP.mean,
#' equation="logistic", replicate.CI=NULL))
#' plot(tsdL_IP, xlab="Incubation durations in days")
#' # Example with Chelonia mydas
#' cm <- subset(DatabaseTSD, Species=="Chelonia mydas" & !is.na(Sexed), c("Males", "Females",
#' "Incubation.temperature", "RMU.2010"))
#' tsd(subset(cm, subset=RMU.2010=="Pacific, SW"))
#' tsd(subset(cm, subset=RMU.2010=="Pacific, Northwest"))
#' tsd(subset(cm, subset=RMU.2010=="Atlantic, S Caribbean"))
#'
#' ### Eretmochelys imbricata
#' Ei_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" &
#' Species=="Eretmochelys imbricata")
#' Ei_AtlanticW <- subset(DatabaseTSD, RMU.2010=="Atlantic, W (Caribbean and E USA)" &
#' Species=="Eretmochelys imbricata")
#' Ei_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" &
#' Species=="Eretmochelys imbricata")
#' Ei_PacSW <- tsd(Ei_PacificSW)
#' Ei_AtlW <- tsd(Ei_AtlanticW)
#' Ei_AtlSW <- tsd(Ei_AtlanticSW)
#'
#' plot(Ei_PacSW, xlim=c(27, 33), show.PTRT = FALSE, main=expression(italic("Eretmochelys imbricata")))
#' par(new=TRUE)
#' plot(Ei_AtlW, xlim=c(27, 33), col="red", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="red")
#' par(new=TRUE)
#' plot(Ei_AtlSW, xlim=c(27, 33), col="blue", xlab="", ylab="", axes=FALSE,
#' xaxt="n", show.PTRT = FALSE, errbar.col="blue")
#' legend("topright", legend=c("Pacific, SW", "Atlantic, W", "Atlantic, SW"), lty=1,
#' col=c("black", "red", "blue"))
#'
#' ### Chelonia mydas
#' Cm_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_PacificNW <- subset(DatabaseTSD, RMU.2010=="Pacific, NW" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_AtlanticSC <- subset(DatabaseTSD, RMU.2010=="Atlantic, S Caribbean" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_IndianSE <- subset(DatabaseTSD, RMU.2010=="Indian, SE" & !is.na(Sexed) &
#' Species=="Chelonia mydas")
#' Cm_PacSW <- tsd(Cm_PacificSW)
#' Cm_PacNW <- tsd(Cm_PacificNW)
#' Cm_IndSE <- tsd(Cm_IndianSE)
#' Cm_AtlSC <- tsd(Cm_AtlanticSC)
#'
#' plot(Cm_PacSW, xlim=c(24, 34), show.PTRT = FALSE, main=expression(italic("Chelonia mydas")))
#' par(new=TRUE)
#' plot(Cm_PacNW, xlim=c(24, 34), col="red", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="red")
#' par(new=TRUE)
#' plot(Cm_IndSE, xlim=c(24, 34), col="blue", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="blue")
#' par(new=TRUE)
#' plot(Cm_AtlSC, xlim=c(24, 34), col="green", xlab="", ylab="",
#' axes=FALSE, xaxt="n", show.PTRT = FALSE, errbar.col="green")
#'
#' # To fit a TSDII or FMF TSD pattern, you must indicate P_low, S_low, P_high, and S_high
#' # for logistic model and P_low, S_low, K1_low, K2_low, P_high, S_high, K1_high, and K2_high for
#' # flexit model
#' # The model must be 0-1 for low and 1-0 for high with P_low < P_high
#'
#' Chelydra_serpentina <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Chelydra serpentina")
#'
#' model_TSDII <- tsd(Chelydra_serpentina, males.freq=FALSE,
#' parameters.initial=c(P_low=21, S_low=0.3, P_high=28, S_high=-0.4),
#' equation="logistic")
#' plot(model_TSDII, lab.TRT = "TRT l = 5 %")
#' priors <- tsd_MHmcmc_p(result=model_TSDII, accept=TRUE)
#' out_mcmc <- tsd_MHmcmc(result=model_TSDII, n.iter=10000, parametersMCMC=priors)
#' plot(model_TSDII, resultmcmc=out_mcmc, lab.TRT = "TRT l = 5 %")
#' predict(model_TSDII, temperatures=25:35)
#'
#' # Podocnemis expansa
#' Podocnemis_expansa <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Podocnemis expansa")
#' Podocnemis_expansa_Valenzuela_2001 <- subset(Podocnemis_expansa,
#' Reference=="Valenzuela, 2001")
#' PeL2001 <- tsd(df=Podocnemis_expansa_Valenzuela_2001)
#' # The pivotal temperature is 32.133 °C (CI 95% 31.495;32.766)
#' # In Valenzuela, 2001: "Using data from the present study alone,
#' # the critical temperature was 32.2 °C by both methods and the 95%
#' # confidence limits were 31.4 °C and 32.9 °C."
#' # Data are close but not identical to what was published.
#'
#' # The pivotal temperature calculated by maximum likelihood and by inverse
#' # prediction from logistic regression, was 32.6°C using raw data from
#' # 1991 (N. Valenzuela, unpublished data) and from this study. The lower
#' # and upper 95% confidence limits of the pivotal temperature were 32.2°C
#' # and 33.2°C,
#'
#' Podocnemis_expansa_Valenzuela_1997 <- subset(Podocnemis_expansa,
#' subset=(((Reference=="Lance et al., 1992; Valenzuela et al., 1997") |
#' (Reference=="Valenzuela, 2001")) &
#' (!is.na(Sexed)) & (Sexed != 0)))
#'
#' PeL1997 <- tsd(df=Podocnemis_expansa_Valenzuela_1997)
#'
#' # Gekko japonicus
#'
#' Gekko_japonicus <- subset(DatabaseTSD, !is.na(Sexed) & (Sexed != 0) &
#' Species=="Gekko japonicus")
#' model_TSDII_gj <- tsd(Gekko_japonicus, males.freq=TRUE,
#' parameters.initial=c(P_low=26, S_low=1.5,
#' P_high=31, S_high=-1.5),
#' equation="logistic")
#' plot(model_TSDII_gj, lab.TRT = "TRT l = 5 %")
#' print(model_TSDII_gj)
#' prior <- tsd_MHmcmc_p(result = model_TSDII_gj, accept = TRUE)
#' prior <- structure(list(
#' Density = c("dnorm", "dnorm", "dnorm", "dnorm"),
#' Prior1 = c(26, 0.3, 31, -0.4),
#' Prior2 = c(2, 1, 2, 1),
#' SDProp = c(2, 0.5, 2, 0.5),
#' Min = c(25, -2, 25, -2),
#' Max = c(35, 2, 35, 2),
#' Init = c(26, 0.3, 31, -0.4)),
#' row.names = c("P_low", "S_low", "P_high", "S_high"),
#' class = "data.frame")
#'
#' result_mcmc_tsd_gj <- tsd_MHmcmc(result=model_TSDII_gj,
#' parametersMCMC=prior, n.iter=10000, n.chains = 1,
#' n.adapt = 0, thin=1, trace=FALSE, adaptive=TRUE)
#' summary(result_mcmc_tsd_gj)
#' plot(result_mcmc_tsd_gj, parameters="P_low", scale.prior=TRUE, xlim=c(20, 30), las=1)
#' plot(result_mcmc_tsd_gj, parameters="P_high", scale.prior=TRUE, xlim=c(25, 35), las=1)
#' plot(model_TSDII_gj, resultmcmc = result_mcmc_tsd_gj)
#'
#' # Trachemys scripta elegans
#' # Take care, the pattern reflects large population variation
#'
#' Tse <- subset(DatabaseTSD, Species=="Trachemys scripta" & Subspecies == "elegans" & !is.na(Sexed))
#' Tse_logistic <- tsd(Tse)
#' plot(Tse_flexit)
#' compare_AICc(logistic=Tse_logistic, flexit=Tse_flexit)
#' plot(Tse_flexit)
#'
#' # Exemple when only proportion is known; experimental
#' Ei_PacificSW <- subset(DatabaseTSD, RMU.2010=="Pacific, SW" &
#' Species=="Eretmochelys imbricata")
#' males <- Ei_PacificSW$Males/(Ei_PacificSW$Males+Ei_PacificSW$Females)*100
#' females <- 100-(Ei_PacificSW$Males/(Ei_PacificSW$Males+Ei_PacificSW$Females)*100)
#' temperatures <- Ei_PacificSW$Incubation.temperature
#' Ei_PacSW <- tsd(Ei_PacificSW)
#' par <- c(Ei_PacSW$par, n=10)
#' embryogrowth:::.tsd_fit(par=par, males=males, N=males+females, temperatures=temperatures,
#' equation="logistic")
#' Ei_PacSW_NormalApproximation <- tsd(males=males, females=females,
#' temperatures=temperatures,
#' parameters.initial=par)
#' Ei_PacSW_NormalApproximation$par
#' Ei_PacSW$par
#' # The data looks like only n=0.01 observations were done
#' # This is the reason of the large observed heterogeneity
#' plot(Ei_PacSW_NormalApproximation)
#' }
#' @export
tsd <- function(df=NULL ,
males=NULL ,
females=NULL ,
N=NULL ,
temperatures=NULL ,
durations=NULL ,
l=0.05 ,
parameters.initial=c(P=30, S=-2, K=0, K1=1, K2=0,
SL = -1, SH = -1) ,
males.freq=TRUE ,
fixed.parameters=NULL ,
equation="logistic" ,
replicate.CI=10000 ,
range.CI=0.95 ,
SE=TRUE ,
replicate.NullDeviance=1000 ,
control=list(maxit=1000) ,
print=TRUE ,
method="BFGS" ) {
# df=NULL; males=NULL; females=NULL; N=NULL; temperatures=NULL; durations=NULL; l=0.05; parameters.initial=c(P=NA, S=-0.5, K=0, K1=1, K2=0, SL = -1, SH = -1); males.freq=TRUE; fixed.parameters=NULL; equation="logistic"; replicate.CI=10000; range.CI=0.95; SE=TRUE;print=TRUE; control=list(maxit=1000); replicate.NullDeviance=1000; method="BFGS"
# CC_AtlanticSW <- subset(DatabaseTSD, RMU.2010=="Atlantic, SW" & Species=="Caretta caretta" & Sexed!=0)
# males=CC_AtlanticSW$Males; females=CC_AtlanticSW$Females; temperatures=CC_AtlanticSW$Incubation.temperature; equation="logistic"
# parameters.initial = c(P=29.19); fixed.parameters = c(S=-0.21)
equation <- tolower(equation)
if (is.null(replicate.CI)) replicate.CI <- 0
equation <- match.arg(equation,
choices = c("logistic", "hill", "a-logistic", "hulin",
"double-a-logistic", "flexit", "flexit*", "gsd",
"probit", "logit"),
several.ok = FALSE)
# method <- tolower(method)
method <- match.arg(method, choices=c("BFGS", "Nelder-Mead"),
several.ok = FALSE)
if (!is.null(df)) {
if ((!inherits(df, "data.frame")) & (!inherits(df, "matrix"))) {
stop("df parameter must be a data.frame or a matrix")
}
namesdf <- tolower(colnames(df))
males <- NULL
females <- NULL
N <- NULL
if (any(namesdf=="males")) males <- df[, which(namesdf=="males")]
if (any(namesdf=="male")) males <- df[, which(namesdf=="male")]
if (any(namesdf=="females")) females <- df[, which(namesdf=="females")]
if (any(namesdf=="female")) females <- df[, which(namesdf=="female")]
if (any(namesdf=="n")) N <- df[, which(namesdf=="n")]
if (any(namesdf=="sexed")) N <- df[, which(namesdf=="sexed")]
if (any(namesdf=="temperatures")) temperatures <- df[, which(namesdf=="temperatures")]
if (any(namesdf=="temperature")) temperatures <- df[, which(namesdf=="temperature")]
if (any(namesdf=="incubation.temperature.set")) temperatures <- df[, which(namesdf=="incubation.temperature.set")]
if (any(namesdf=="durations")) durations <- df[, which(namesdf=="durations")]
if (any(namesdf=="duration")) durations <- df[, which(namesdf=="duration")]
if (any(namesdf=="IP.mean")) durations <- df[, which(namesdf=="IP.mean")]
}
if (is.null(durations) & is.null(temperatures)) {
stop("Or temperatures or durations must be provided")
}
if (!is.null(durations) & !is.null(temperatures)) {
stop("Temperatures and durations cannot be both provided")
}
# si je retire des lignes, c'est decale
# if (!is.null(males)) males <- na.omit(males)
# if (!is.null(females)) females <- na.omit(females)
# if (!is.null(N)) N <- na.omit(N)
if (is.null(temperatures)) {
IP <- TRUE
temperatures <- durations
} else {
IP <- FALSE
}
cpt1 <- ifelse(is.null(males), 0, 1)+ifelse(is.null(females), 0, 1)+ifelse(is.null(N), 0, 1)
if (cpt1<2) {
stop("At least two informations from males, females or N must be provided")
}
if (is.null(males)) males <- N-females
if (is.null(females)) females <- N-males
if (is.null(N)) N <- males+females
if (length(temperatures) != length(N)) {
stop("A temperature or duration must be provided for each experiment")
}
if (all(males * females == 0)) {
warning("At least one incubation temperature should produce mixed sex ratio. Use MCMC rather.")
}
if (equation == "gsd") {
value <- -sum(dbinom(x=males, size=males+females, prob=0.5, log=TRUE))
result <- list(par = NULL, SE=NULL, hessian=NULL,
TRT=NULL,
SE_TRT=NULL,
message=NULL,
convergence=0,
counts=NULL,
value = value,
AIC= 2*value + 2*0,
AICc= 2*value + (2*0*(0+1))/(sum(N)-0-1),
BIC= 2*value + 0
)
# limit.low.TRT <- min(temperatures)
# limit.high.TRT <- max(temperatures)
} else {
ppi <- c(parameters.initial, fixed.parameters)
if (is.na(ppi["P_low"])) {
if (is.na(ppi["P"])) {
if ((equation!="probit") & (equation!="logit")) {
ppi["P"] <- temperatures[which.min(abs((males/(males+females))-0.5))]
} else {
ppi["P"] <- 0
}
}
if (is.na(ppi["S"])) {
if ((equation!="probit")) {
pente <- lm(males/(males+females) ~ temperatures)
ppi["S"] <- pente$coefficients["temperatures"]
if ((equation != "flexit") & (equation != "flexit*")) ppi["S"] <- 1/(4*ppi["S"])
} else {
ppi["S"] <- 0
}
}
}
# if (IP) ppi["S"] <- abs(ppi["S"])
if (equation != "a-logistic") ppi <- ppi[which(names(ppi) != "K")]
if ((equation != "hulin") & (equation != "double-a-logistic") & (equation != "flexit")) {
ppi <- ppi[which(names(ppi) != "K1")]
ppi <- ppi[which(names(ppi) != "K2")]
}
if (equation != "flexit*") {
ppi <- ppi[which(names(ppi) != "SL")]
ppi <- ppi[which(names(ppi) != "SH")]
} else {
ppi <- ppi[which(names(ppi) != "S")]
}
if (!is.null(fixed.parameters))
ppi[names(fixed.parameters)] <- unname(fixed.parameters)
repeat {
# result <- optim(par, embryogrowth:::.tsd_fit, fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures, equation=equation, method="BFGS", hessian=TRUE, control = list(maxit=1000))
result <- try(optim(ppi, getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures,
equation=equation, method=method, hessian=FALSE, control = control), silent=TRUE)
if (inherits(result, "try-error")) {
stop("Try using method='Nelder-Mead' or use MCMC but it will be complicated ! Contact me.")
}
if (result$convergence != 1) break
ppi <- result$par
if (print) print("Convergence is not acheived. Optimization continues !")
}
result_hessian <- try(optim(result$par, getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=fixed.parameters, males=males, N=N, temperatures=temperatures,
equation=equation, method=method, hessian=TRUE, control = control), silent=TRUE)
ppi <- c(result$par, fixed.parameters)
if (inherits(result_hessian, "try-error")) {
warning("Likelihood is flat close to ML point and Hessian cannot be estimated")
result$hessian <- NULL
result$SE <- NULL
replicate.CI = 0
} else {
result$hessian <- result_hessian$hessian
if (SE) {
rh <- SEfromHessian(result$hessian, hessian=TRUE)
result$SE <- rh$SE
}
}
# result$hessian <- result$hessian
result$AIC <- 2*result$value+2*length(result$par)
# Correction le 21/10/2020
result$AICc <- result$AIC +(2*length(result$par)*(length(result$par)+1))/(sum(N)-length(result$par)-1)
# Correction le 21/10/2020
result$BIC <- 2*result$value+ length(result$par)*log(sum(N))
}
result$range.CI <- range.CI
result$l <- l
result$equation <- equation
result$males <- males
result$females <- females
result$N <- N
result$temperatures <- temperatures
result$males.freq <- males.freq
result$fixed.parameters <- fixed.parameters
result$type <- ifelse(IP, "duration", "temperature")
if (all(names(ppi) != "n")) {
result$deviance <- -2*(-result$value - sum(dbinom(x = males, size=N,
prob = males/N, log = TRUE)))
} else {
pt <- males/N
pt <- ifelse(pt<=1E-9, 1E-9, pt)
pt <- ifelse(pt>=1-1E-9, 1-1E-9, pt)
sd <- sqrt((pt)*(1-(pt))/ppi["n"])
pr <- dnorm(males/N, mean=pt, sd=sd, log=TRUE)
result$deviance <- -2*(-result$value - sum(pr))
}
# degrees of freedom calculated from the difference of the number of parameters in the saturated and the fitted model.
result$df <- length(males) - length(result$par)
result$pvalue <- pchisq(q=result$deviance, df=result$df, lower.tail = FALSE)
if (equation != "gsd") {
result <- addS3Class(result, "tsd")
# result <<- result
if (!is.null(result$hessian)) {
vcov <- solve(result$hessian)
if (any(eigen(vcov)$values < 0)) {
warning("Non-positive definite variance-covariance matrix; use MCMC to get credible interval.")
result$hessian <- NULL
result$SE <- NULL
replicate.CI = 0
}
}
o <- P_TRT(x=result, l=l,
replicate.CI = replicate.CI,
probs = c((1-range.CI)/2, 0.5, 1-(1-range.CI)/2))
result$P_TRT <- o$P_TRT_quantiles
}
probtheor <- getFromNamespace(x=".modelTSD", ns="embryogrowth")(c(result$par, result$fixed.parameters),
result$temperatures,
equation=equation)
result$predictMaleFrequency <- probtheor
if (!is.null(replicate.NullDeviance) & (replicate.NullDeviance != 0)) {
pt <- NULL
# Bug correct 2020-01-11: change 1000 to replicate.NullDeviance
for (i in 1:replicate.NullDeviance) {
m <- rbinom(n=length(result$males), size=result$N, prob = probtheor)
result_dev <- optim(par=result$par, fn=getFromNamespace(".tsd_fit", ns="embryogrowth"),
fixed.parameters=result$fixed.parameters, males=m, N=result$N,
temperatures=result$temperatures,
equation=result$equation, method=method,
hessian=FALSE, control = control)
deviance_dev <- -2*(-result$value - sum(dbinom(x = m, size=N,
prob = m/N, log = TRUE)))
pt <- c(pt, deviance_dev)
}
result$NullDeviance <- pt
result$NullDeviancePvalue <- sum(result$deviance < pt)/length(pt)
}
if (equation!="gsd" & print) {
# SI je n'ai qu'un PT c'est que c'est un profile TSDIA ou IB
if (any(colnames(o$P_TRT_quantiles)=="PT")) {
if (!is.null(replicate.CI) & (replicate.CI != 0)) {
print(paste0("The pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[3, "PT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT"], o$P_TRT_quantiles[3, "PT"]))))
print(paste0("The transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT"], o$P_TRT_quantiles[3, "TRT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT"], o$P_TRT_quantiles[3, "TRT"]))))
print(paste0("The lower limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT"], o$P_TRT_quantiles[3, "lower.limit.TRT"]))))
print(paste0("The upper limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT"], o$P_TRT_quantiles[3, "higher.limit.TRT"]))))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
} else {
print(paste0("The pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT"])))
print(paste0("The transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT"])))
print(paste0("The lower limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT"])))
print(paste0("The upper limit of transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT"])))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
}
} else {
if (!is.null(replicate.CI) & (replicate.CI != 0)) {
print(paste0("The lower pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[3, "PT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT_low"], o$P_TRT_quantiles[3, "PT_low"]))))
print(paste0("The lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT_low"], o$P_TRT_quantiles[3, "TRT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT_low"], o$P_TRT_quantiles[3, "TRT_low"]))))
print(paste0("The lower limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT_low"], o$P_TRT_quantiles[3, "lower.limit.TRT_low"]))))
print(paste0("The upper limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_low"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT_low"], o$P_TRT_quantiles[3, "higher.limit.TRT_low"]))))
print(paste0("The upper pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f", min(o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[3, "PT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "PT_high"], o$P_TRT_quantiles[3, "PT_high"]))))
print(paste0("The upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "TRT_high"], o$P_TRT_quantiles[3, "TRT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "TRT_high"], o$P_TRT_quantiles[3, "TRT_high"]))))
print(paste0("The lower limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"])), ";", sprintf("%.3f", max(o$P_TRT_quantiles[1, "lower.limit.TRT_high"], o$P_TRT_quantiles[3, "lower.limit.TRT_high"]))))
print(paste0("The upper limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_high"]), " CI", as.character(range.CI*100), "% ", sprintf("%.3f",min(o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"])), ";", sprintf("%.3f",max(o$P_TRT_quantiles[1, "higher.limit.TRT_high"], o$P_TRT_quantiles[3, "higher.limit.TRT_high"]))))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
} else {
print(paste0("The lower pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_low"])))
print(paste0("The lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_low"])))
print(paste0("The lower limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_low"])))
print(paste0("The upper limit of lower transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_low"])))
print(paste0("The upper pivotal ", result$type, " is ", sprintf("%.3f",o$P_TRT_quantiles[2, "PT_high"])))
print(paste0("The upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "TRT_high"])))
print(paste0("The lower limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "lower.limit.TRT_high"])))
print(paste0("The upper limit of upper transitional range of ", result$type, "s is ", sprintf("%.3f",o$P_TRT_quantiles[2, "higher.limit.TRT_high"])))
if (!is.na(result$par["S"])) print(paste0("The S parameter value is ", sprintf("%.3f",result$par["S"])))
}
}
}
result <- addS3Class(result, "tsd")
return(invisible(result))
}
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.