tsd: Estimate the parameters that best describe...
In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth
tsd R Documentation
Estimate the parameters that best describe temperature-dependent sex determination
Description
Estimate the parameters that best describe the thermal reaction norm for sex ratio when temperature-dependent sex determination occurs.
It can be used also to evaluate the relationship between incubation duration and sex ratio.
The parameter l was defined in Girondot (1999). The TRT is defined from the difference between the two boundary temperatures giving sex ratios of l and 1 - l, respectively:
For logistic model (Girondot, 1999), it follows
TRT_{l}=abs\left ( S\: K_{l} \right )
where K_{l} is a constant equal to 2\: log\left ( \frac{l}{1-l} \right ).
In Girondot (1999), l was 0.05 and then the TRT was defined as being the range of temperatures producing from 5\
The default model is named logistic. This model (as well as the logit one) has the particularity to have a symmetric shape around P.
The logistic model is:
SR(T) = 1 / (1 + exp((1 / S) * (P - T))))
The logit model is:
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.
Hulin model is assymmetric but the control of asymmetry is difficult to manage.
If asymmetric model is selected, it is always better to use flexit model.
if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}
if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}
with:
S_1 = S/((4/K_1)*(2^{(-K_1)})^{(1/K_1+1)}*(2^{K_1}-1))
S_2 = S/((4/K_2)*(2^{(-K_2)})^{(1/K_2+1)}*(2^{K_2}-1))
The flexit* model is defined as (QBT is the Quasi-Binary Threshold):
QBT = 1/ (1 + exp(100 * (P - T)))
SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))
Usage
tsd(
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"
)
Arguments
df
A dataframe with at least two columns named males, females or N and temperatures, Incubation.temperature or durations column
males
A vector with male numbers
females
A vector with female numbers
N
A vector with total numbers
temperatures
The constant incubation temperatures used to fit sex ratio
durations
The duration of incubation or TSP used to fit sex ratio
l
Sex ratio limits to define TRT are l and 1-l (see Girondot, 1999)
parameters.initial
Initial values for P, S or K search as a vector, ex. c(P=29, S=-0.3)
males.freq
If TRUE data are shown as males frequency
fixed.parameters
Parameters that will not be changed
equation
Can be "logistic", "Hill", "A-logistic", "Hulin", "Double-A-logistic", "flexit", "flexit*", "GSD", "logit", "probit"
replicate.CI
Number of replicates to estimate confidence intervals
range.CI
The range of confidence interval for estimation, default=0.95
SE
If FALSE, does not estimate SE of parameters. Can be use when something wrong happens.
replicate.NullDeviance
Number of replicates to estimate null distribution of deviance
control
List of parameters used in optim.
print
Should the results be printed at screen? TRUE (default) or FALSE
method
method used for optim. Can be "BFGS", the most rapid or "Nelder-Mead" for special cases using n parameter.
Details
tsd estimates the parameters that best describe temperature-dependent sex determination
Value
A list the pivotal temperature, transitional range of temperatures and their SE
Author(s)
Marc Girondot marc.girondot@gmail.com
References
\insertRef1515embryogrowth
\insertRef3534embryogrowth
\insertRef11754embryogrowth
\insertRef5790embryogrowth
See Also
Other Functions for temperature-dependent sex determination:
DatabaseTSD,
DatabaseTSD.version(),
P_TRT(),
ROSIE,
ROSIE.version(),
TSP.list,
plot.tsd(),
predict.tsd(),
stages,
tsd_MHmcmc(),
tsd_MHmcmc_p()
Examples
## Not run:
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)
## End(Not run)
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.
| tsd | R Documentation |
Estimate the parameters that best describe temperature-dependent sex determination
Description
Estimate the parameters that best describe the thermal reaction norm for sex ratio when temperature-dependent sex determination occurs.
It can be used also to evaluate the relationship between incubation duration and sex ratio.
The parameter l was defined in Girondot (1999). The TRT is defined from the difference between the two boundary temperatures giving sex ratios of l and 1 - l, respectively:
For logistic model (Girondot, 1999), it follows
TRT_{l}=abs\left ( S\: K_{l} \right )
where K_{l} is a constant equal to 2\: log\left ( \frac{l}{1-l} \right ).
In Girondot (1999), l was 0.05 and then the TRT was defined as being the range of temperatures producing from 5\
The default model is named logistic. This model (as well as the logit one) has the particularity to have a symmetric shape around P.
The logistic model is:
SR(T) = 1 / (1 + exp((1 / S) * (P - T))))
The logit model is:
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.
Hulin model is assymmetric but the control of asymmetry is difficult to manage.
If asymmetric model is selected, it is always better to use flexit model.
if\ \ T < P\ \ then\ \ (1 + (2^{K_1} - 1) * exp(4 * S_1 * (P - T)))^{(-1/K_1)}
if\ \ T > P\ \ then\ \ 1-((1 + (2^{K_2} - 1) * exp(4 * S_2 * (T - P)))^{(-1/K_2)}
with:
S_1 = S/((4/K_1)*(2^{(-K_1)})^{(1/K_1+1)}*(2^{K_1}-1))
S_2 = S/((4/K_2)*(2^{(-K_2)})^{(1/K_2+1)}*(2^{K_2}-1))
The flexit* model is defined as (QBT is the Quasi-Binary Threshold):
QBT = 1/ (1 + exp(100 * (P - T)))
SR(T) = 1 / (1 + exp(4 * (SL * QBT + SH * (1 - QBT)) * (P - T))))
Usage
tsd(
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"
)
Arguments
df |
A dataframe with at least two columns named males, females or N and temperatures, Incubation.temperature or durations column |
males |
A vector with male numbers |
females |
A vector with female numbers |
N |
A vector with total numbers |
temperatures |
The constant incubation temperatures used to fit sex ratio |
durations |
The duration of incubation or TSP used to fit sex ratio |
l |
Sex ratio limits to define TRT are l and 1-l (see Girondot, 1999) |
parameters.initial |
Initial values for P, S or K search as a vector, ex. c(P=29, S=-0.3) |
males.freq |
If TRUE data are shown as males frequency |
fixed.parameters |
Parameters that will not be changed |
equation |
Can be "logistic", "Hill", "A-logistic", "Hulin", "Double-A-logistic", "flexit", "flexit*", "GSD", "logit", "probit" |
replicate.CI |
Number of replicates to estimate confidence intervals |
range.CI |
The range of confidence interval for estimation, default=0.95 |
SE |
If FALSE, does not estimate SE of parameters. Can be use when something wrong happens. |
replicate.NullDeviance |
Number of replicates to estimate null distribution of deviance |
control |
List of parameters used in optim. |
print |
Should the results be printed at screen? TRUE (default) or FALSE |
method |
method used for optim. Can be "BFGS", the most rapid or "Nelder-Mead" for special cases using n parameter. |
Details
tsd estimates the parameters that best describe temperature-dependent sex determination
Value
A list the pivotal temperature, transitional range of temperatures and their SE
Author(s)
Marc Girondot marc.girondot@gmail.com
References
\insertRef1515embryogrowth
\insertRef3534embryogrowth
\insertRef11754embryogrowth
\insertRef5790embryogrowth
See Also
Other Functions for temperature-dependent sex determination:
DatabaseTSD,
DatabaseTSD.version(),
P_TRT(),
ROSIE,
ROSIE.version(),
TSP.list,
plot.tsd(),
predict.tsd(),
stages,
tsd_MHmcmc(),
tsd_MHmcmc_p()
Examples
## Not run:
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)
## End(Not run)
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.