embryogrowth (version 8.0)

searchR: Fit the parameters that best represent nest incubation data.

Description

Fit the parameters that best represent data. test can be a data.frame with two columns Mean and SD and rownames with the nest name. If SD is na, then least squarre criteria is used for fitting. Function to fit thermal reaction norm can be expressed as :

• a 4-parameters Schoolfield, Sharpe, and Magnuson model (1981) with DHH, DHA, T12H, and Rho25;

• a 6-parameters Schoolfield, Sharpe, and Magnuson model (1981) with T12L, DT, DHH, DHL, DHA, and Rho25;

• Each of these two first models can be combined as low and high sets of parameters by adding the _L suffix to one set. Then you must add also transition_S and transition_P parameters and then the growth rate is 1/(1+exp((1/transition_S)*(P-transition_P))) with P being the proportion of development;

• The Rho25_b control the effect of hygrometry (or Rho25_b_L) (It is not fully functional still);

• a Weibull function with k (shape), lambda (scale) and theta parameters;

• a normal function with Peak, Scale, and sd parameters;

• an asymmetric normal fuction with Peak, Scale, sdH and sdL parameters;

• a symmetric trigonometric function with Length, Peak, and Max;

• an asymmetric trigonometric function with LengthB, LengthE, Peak, and Max.

• Dallwitz-Higgins model (1992) can be used using Dallwitz_b1, Dallwitz_b2, Dallwitz_b3, Dallwitz_b4 and Dallwitz_b5 parameters.

• If Dallwitz_b4 is not included, Dallwitz_b4 = 6 will be used.

• If Dallwitz_b5 is not included, Dallwitz_b5 = 0.4 will be used.

• It is possible also to add the parameter epsilon and then the model becomes X + epsilon with X being any of the above model;

• It is possible also to add the parameter epsilon_L and then the model becomes X_L + epsilon_L with X_L being any of the above model with suffix _L;

• If the name of the parameter is a number, then the mdel is a polynom anchored with the rate being the parameter value at this temperature (the name). see ChangeSSM() function.

Usage

searchR(
  parameters = stop("Initial set of parameters must be provided"),
  fixed.parameters = NULL,
  temperatures = stop("Formated temperature must be provided !"),
  derivate = dydt.Gompertz,
  test = c(Mean = 39.33, SD = 1.92),
  M0 = 1.7,
  saveAtMaxiter = FALSE,
  fileName = "intermediate",
  weight = NULL,
  control = list(trace = 1, REPORT = 100, maxit = 500)
)

Arguments

parameters

A set of parameters used as initial point for searching

fixed.parameters

A set of parameters that will not be changed

temperatures

Timeseries of temperatures after formated using FormatNests()

derivate

Function used to fit embryo growth: dydt.Gompertz, dydt.exponential or dydt.linear

test

A vector with Mean and SD of size of hatchlings, ex. test=c(Mean=39, SD=3). Can be a data.frame also. See description

M0

Measure of hatchling size or mass proxi at laying date

saveAtMaxiter

If True, each time number of interation reach maxiter, current data are saved in file with filename name

fileName

The intermediate results are saved in file with fileName.Rdata name

weight

A named vector of the weight for each nest for likelihood estimation

control

List for control parameters for optimx

Value

A result object

Details

searchR fits the parameters that best represent nest incubation data.

References

Angilletta, M.J., 2006. Estimating and comparing thermal performance curves. Journal of Thermal Biology 31, 541-545.

Dallwitz, M.J., Higgins, J.P., 1992. User<U+2019>s guide to DEVAR. A computer program for estimating development rate as a function of temperature. CSIRO Aust Div Entomol Rep 2, 1-23.

Georges, A., Beggs, K., Young, J.E., Doody, J.S., 2005. Modelling development of reptile embryos under fluctuating temperature regimes. Physiological and Biochemical Zoology 78, 18-30.

Girondot, M., Kaska, Y., 2014. A model to predict the thermal reaction norm for the embryo growth rate from field data. Journal of Thermal Biology 45, 96-102.

Schoolfield, R.M., Sharpe, P.J., Magnuson, C.E., 1981. Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory. Journal of Theoretical Biology 88, 719-731.

Examples

# NOT RUN {
library(embryogrowth)
data(nest)
formated <- FormatNests(nest)
# The initial parameters value can be:
# "T12H", "DHA",  "DHH", "Rho25"
# Or
# "T12L", "DT", "DHA",  "DHH", "DHL", "Rho25"
# K for Gompertz must be set as fixed parameter or being a constant K  
# or relative to the hatchling size rK
x <- structure(c(105.966881676793, 613.944134764125, 306.449533440186, 
                118.193882815108), .Names = c("DHA", "DHH", "T12H", "Rho25"))
# pfixed <- c(K=82.33) or rK=82.33/39.33
pfixed <- c(rK=2.093313)
resultNest_4p_SSM4p <- searchR(parameters=x, fixed.parameters=pfixed, 
	temperatures=formated, derivate=dydt.Gompertz, M0=1.7, 
	test=c(Mean=39.33, SD=1.92))
data(resultNest_4p_SSM4p)
plot(resultNest_4p_SSM4p, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,45), series=1, 
embryo.stages="Caretta caretta.SCL")
x <- structure(c(106.567809092008, 527.359011254683, 614.208632495199, 
2720.94506457237, 306.268259715624, 120.336791245212), .Names = c("DHA", 
"DHH", "DHL", "DT", "T12L", "Rho25"))
pfixed <- c(rK=2.093313)

# exemple of data.frame for test
ttest <- data.frame(Mean=rep(25.5, formated$IndiceT["NbTS"]), 
                     SD=rep(0.75, formated$IndiceT["NbTS"]), 
                     row.names=names(formated)[1:formated$IndiceT["NbTS"]])
resultNest_6p_SSM6p <- searchR(parameters=x, fixed.parameters=pfixed, 
	                        temperatures=formated, derivate=dydt.Gompertz, M0=1.7, 
	                        test=ttest)
	                        
data(resultNest_6p_SSM6p)
pMCMC <- TRN_MHmcmc_p(resultNest_6p_SSM6p, accept=TRUE)
# Take care, it can be very long, sometimes several days
resultNest_mcmc_6p_SSM6p <- GRTRN_MHmcmc(result=resultNest_6p_SSM6p,  
	parametersMCMC=pMCMC, n.iter=10000, n.chains = 1, n.adapt = 0,  
	thin=1, trace=TRUE)
# compare_AIC() is a function from the package "HelpersMG"
compare_AIC(test1=resultNest_4p_SSM4p, test2=resultNest_6p_SSM6p)
############ with new parametrization
data(resultNest_4p_SSM4p)
x0 <- resultNest_4p_SSM4p$par
t <- hist(resultNest_4p_SSM4p, plot=FALSE)
x <- c(3.4, 3.6, 5.4, 5.6, 7.6, 7.5, 3.2)
names(x) <- seq(from=range(t$temperatures)[1], to=range(t$temperatures)[2], 
     length.out=7)
newx <- ChangeSSM(temperatures = (200:350)/10, parameters = x0, 
       initial.parameters = x, 
       control=list(maxit=5000))
pfixed <- c(rK=2.093313)
resultNest_newp <- searchR(parameters=newx, fixed.parameters=pfixed,
 temperatures=formated, derivate=dydt.Gompertz, M0=1.7,
 test=c(Mean=39.33, SD=1.92))
plotR(resultNest_newp, ylim=c(0, 3), xlimR=c(23, 34), ylimH=c(0, 0.3), show.hist.TRUE)
compare_AIC(test4p=resultNest_4p_SSM4p, 
            test6p=resultNest_6p_SSM6p, 
            testAnchor=resultNest_newp)
            
############################################
# example with thermal reaction norm fitted from Weibull function
############################################

 x <- ChangeSSM(temperatures = (200:350)/10,
                parameters = resultNest_4p_SSM4p$par,
                initial.parameters = structure(c(73.4009010417375, 304.142079511996, 
                                                27.4671689276281), 
                                        .Names = c("k", "lambda", "scale")), 
                control=list(maxit=1000))
pfixed <- c(rK=2.093313)
resultNest_3p_Weibull <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=1.7, 
                         test=c(Mean=39.33, SD=1.92))
plotR(list(resultNest_4p_SSM4p, resultNest_3p_Weibull), ylim=c(0,3), col=c("Black", "red"))
compare_AIC(SSM4p=resultNest_4p_SSM4p, Weibull=resultNest_3p_Weibull)

###########################################
# example with thermal reaction norm fitted from asymmetric normal function
############################################

x <- ChangeSSM(temperatures = (200:350)/10,
               parameters = resultNest_4p_SSM4p$par,
               initial.parameters = structure(c(3, 7, 11, 32), 
                               .Names = c("Scale", "sdL", "sdH", "Peak")), 
               control=list(maxit=1000))
pfixed <- c(rK=2.093313)
resultNest_4p_normal <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=1.7, 
                         test=c(Mean=39.33, SD=1.92))
                         
###########################################
# example with thermal reaction norm fitted from trigonometric model
############################################

 x <- ChangeSSM(temperatures = (200:350)/10,
               parameters = resultNest_4p_SSM4p$par,
               initial.parameters = structure(c(3, 20, 40, 32), 
               .Names = c("Max", "LengthB", "LengthE", "Peak")), 
               control=list(maxit=1000))
pfixed <- c(rK=2.093313)
resultNest_4p_trigo <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=1.7, 
                         test=c(Mean=39.33, SD=1.92))
                         
###########################################
# example with thermal reaction norm fitted from Dallwitz model
# From Girondot, M., Monsinjon, J., Guillon, J.-M., Submitted. Delimitation of the 
# embryonic thermosensitive period for sex determination using an embryo growth model 
# reveals a potential bias for sex ratio prediction in turtles.
#  rK = 1.208968 
#  M0 = 0.3470893 
# See the example in stages datasets
############################################

x <- structure(c(4.88677476830268, 20.4051904475743, 31.5173105860335), 
.Names = c("Dallwitz_b1", "Dallwitz_b2", "Dallwitz_b3"))
pfixed <- c(rK=1.208968)
resultNest_3p_Dallwitz <- searchR(parameters=x, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=0.3470893, 
                         test=c(Mean=39.33, SD=1.92))
plotR(resultNest_3p_Dallwitz, ylim=c(0,8))

x <- structure(c(4.91191231405918, 12.7453211281394, 31.2670410811077, 
5.7449376569153, -0.825689964543813), .Names = c("Dallwitz_b1", 
"Dallwitz_b2", "Dallwitz_b3", "Dallwitz_b4", "Dallwitz_b5"))
pfixed <- c(rK=1.208968)
resultNest_5p_Dallwitz <- searchR(parameters=x, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=0.3470893, 
                         test=c(Mean=39.33, SD=1.92))
plotR(resultNest_5p_Dallwitz, ylim=c(0,8))

xp <- resultNest_6p_SSM6p$par
xp["Rho25"] <- 233
pfixed <- c(rK=1.208968)
resultNest_6p_SSM <- searchR(parameters=xp, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=0.3470893, 
                         test=c(Mean=39.33, SD=1.92))
plotR(resultNest_6p_SSM, ylim=c(0,8))

xp <- ChangeSSM(parameters = resultNest_3p_Dallwitz$par, 
                initial.parameters = resultNest_4p_SSM4p$par)
pfixed <- c(rK=1.208968)
resultNest_4p_SSM <- searchR(parameters=xp$par, fixed.parameters=pfixed, 
                         temperatures=formated, derivate=dydt.Gompertz, M0=0.3470893, 
                         test=c(Mean=39.33, SD=1.92))
plotR(resultNest_4p_SSM, ylim=c(0,8))

compare_AIC(Dallwitz3p=resultNest_3p_Dallwitz, Dallwitz5p=resultNest_5p_Dallwitz, 
             SSM4p=resultNest_4p_SSM, SSM6p=resultNest_6p_SSM)
                         
###########################################
# Example with thermal reaction norm of proportion of development
# fitted from Dallwitz model
# see Woolgar, L., Trocini, S., Mitchell, N., 2013. Key parameters describing 
# temperature-dependent sex determination in the southernmost population of loggerhead 
# sea turtles. Journal of Experimental Marine Biology and Ecology 449, 77-84.
############################################

x <- structure(c(1.48207559695689, 20.1100310234046, 31.5665036287242), 
 .Names = c("Dallwitz_b1", "Dallwitz_b2", "Dallwitz_b3"))
resultNest_PropDev_3p_Dallwitz <- searchR(parameters=x, fixed.parameters=NULL, 
                         temperatures=formated, derivate=dydt.linear, M0=0, 
                         test=c(Mean=1, SD=NA))
 plotR(resultNest_PropDev_3p_Dallwitz, ylim=c(0, 1.5), curves="ML")
 plot(x=resultNest_PropDev_3p_Dallwitz, ylimS=c(0,1), xlim=c(0,60), series=2, 
         embryo.stages="Generic.ProportionDevelopment")
         
x <- structure(c(1.48904182113431, 10.4170365155993, 31.2591665490154, 
6.32355497589913, -1.07425378667104), .Names = c("Dallwitz_b1", 
"Dallwitz_b2", "Dallwitz_b3", "Dallwitz_b4", "Dallwitz_b5"))
resultNest_PropDev_5p_Dallwitz <- searchR(parameters=x, fixed.parameters=NULL, 
                         temperatures=formated, derivate=dydt.linear, M0=0, 
                         test=c(Mean=1, SD=NA))
 plotR(resultNest_PropDev_5p_Dallwitz, ylim=c(0, 1.5))
 plot(x=resultNest_PropDev_5p_Dallwitz, ylimS=c(0,1), xlim=c(0,60), series=2, 
         embryo.stages="Generic.ProportionDevelopment")
         
 plotR(resultNest_PropDev_3p_Dallwitz, ylim=c(0, 1.5), curves="ML")
 plotR(resultNest_PropDev_5p_Dallwitz, ylim=c(0, 1.5), curves="ML", new=FALSE, col="red")
 compare_AICc(Dallwitz3p=resultNest_PropDev_3p_Dallwitz, 
              Dallwitz5p=resultNest_PropDev_5p_Dallwitz)

# }