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

Description

Fit the parameters that best represent data.
hatchling.metric 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 model 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 = c(rK = 1.208968),
  temperatures = stop("Formated temperature must be provided !"),
  integral = integral.Gompertz,
  derivate = NULL,
  hatchling.metric = c(Mean = 39.33, SD = 1.92),
  M0 = 0.3470893,
  saveAtMaxiter = FALSE,
  fileName = "intermediate",
  weight = NULL,
  control = list(trace = 1, REPORT = 100, maxit = 500)
)

Value

A result object

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()
integral: Function used to fit embryo growth: integral.Gompertz, integral.exponential or integral.linear
derivate: Function used to fit embryo growth: derivate.Gompertz, derivate.exponential or derivate.linear
hatchling.metric: A vector with Mean and SD of size of hatchlings, ex. hatchling.metric=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

Author

Marc Girondot marc.girondot@gmail.com

Details

searchR fits the parameters that best represent nest incubation data.

References

9039embryogrowth
10871embryogrowth
8566embryogrowth
10620embryogrowth
10039embryogrowth
3326embryogrowth
5322embryogrowth

Examples

Run this code

if (FALSE) {
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
############################################################################
# From Girondot M, Monsinjon J, Guillon J-M (2018) Delimitation of the 
# embryonic thermosensitive period for sex determination using an embryo 
# growth model reveals a potential bias for sex ratio prediction in turtles. 
# Journal of Thermal Biology 73: 32-40 
#  rK = 1.208968 
#  M0 = 0.3470893 
############################################################################
pfixed <- c(rK=1.208968)
M0 = 0.3470893 
############################################################################
# 4 parameters SSM
############################################################################
x <- c('DHA' = 109.31113503282113, 'DHH' = 617.80695919563857, 
    'T12H' = 306.38890489505093, 'Rho25' = 229.37265815800225)
    
resultNest_4p_SSM <- searchR(parameters=x, fixed.parameters=pfixed, 
	temperatures=formated, integral=integral.Gompertz, M0=M0, 
	hatchling.metric=c(Mean=39.33, SD=1.92))
	
data(resultNest_4p_SSM)
plot(resultNest_4p_SSM, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,45), series=1, 
embryo.stages="Caretta caretta.SCL")

############################################################################
# 6 parameters SSM
############################################################################
x <- structure(c(106.567809092008, 527.359011254683, 614.208632495199, 
2720.94506457237, 306.268259715624, 120.336791245212), .Names = c("DHA", 
"DHH", "DHL", "DT", "T12L", "Rho25"))

############################################################################
# example of data.frame for hatchling.metric
############################################################################
thatchling.metric <- data.frame(Mean=rep(39.33, formated$IndiceT["NbTS"]), 
                     SD=rep(1.92, formated$IndiceT["NbTS"]), 
                     row.names=names(formated)[1:formated$IndiceT["NbTS"]])
# It is sometimes difficult to find a good starting point for 
# SSM 6 parameters model. This function helps to find it based on a previoulsy 
# fitted model.

 x <- ChangeSSM(temperatures = (200:350)/10,
                parameters = resultNest_4p_SSM$par,
                initial.parameters = structure(c(106.567809092008, 
                527.359011254683, 614.208632495199, 
                4.94506457237, 306.268259715624, 120.336791245212), 
                .Names = c("DHA", "DHH", "DHL", "DT", "T12L", "Rho25")), 
                control=list(maxit=1000))
                     
resultNest_6p_SSM <- searchR(parameters=x$par, fixed.parameters=pfixed, 
	                        temperatures=formated, integral=integral.Gompertz, 
	                        M0=M0, 
	                        hatchling.metric=thatchling.metric)
	                        
plotR(resultNest_6p_SSM, ylim=c(0, 1))
	                        
data(resultNest_6p_SSM)
pMCMC <- TRN_MHmcmc_p(resultNest_6p_SSM, accept=TRUE)
# Take care, it can be very long, sometimes several days
resultNest_mcmc_6p_SSM <- GRTRN_MHmcmc(result=resultNest_6p_SSM,  
	                                      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_SSM, test2=resultNest_6p_SSM)
############################################################################

############################################################################
############ Example as linear progression of development
############ The development progress goes from 0 to 1
############################################################################

pfixed <- NULL
M0 = 0

############################################################################
# 4 parameters SSM
############################################################################
x <- c('DHA' = 64.868697530424186, 'DHH' = 673.18292743646771, 
       'T12H' = 400.90952554047749, 'Rho25' = 82.217237723502123)
resultNest_4p_SSM_Linear <- searchR(parameters=x, fixed.parameters=pfixed, 
	temperatures=formated, integral=integral.linear, M0=M0, 
	hatchling.metric=c(Mean=39.33, SD=1.92)/39.33)
plotR(resultNest_4p_SSM_Linear, ylim=c(0, 2), scaleY= 100000)
plot(resultNest_4p_SSM_Linear, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,1.1), 
series=1, embryo.stages="Generic.ProportionDevelopment")

tc <- GenerateConstInc(duration=300*24*60, temperatures = 28)
tc_f <- FormatNests(tc)
plot(resultNest_4p_SSM_Linear, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,1.1), 
series=1, embryo.stages="Generic.ProportionDevelopment", temperatures=tc_f)

############################################################################
############ with new parametrization based on anchor
############ This is a non-parametric version
############################################################################

data(resultNest_4p_SSM)
x0 <- resultNest_4p_SSM$par
t <- range(hist(resultNest_4p_SSM, plot=FALSE)$temperatures)

x <- getFromNamespace(".SSM", ns="embryogrowth")(T=seq(from=t[1], 
                                                       to=t[2], 
                                                       length.out=7), 
                                            parms=x0)[[1]]*1E5
names(x) <- as.character(seq(from=t[1], 
                          to=t[2], 
                          length.out=7))
     
M0 <- 0.3470893 
pfixed <- c(rK=1.208968)
resultNest_newp <- searchR(parameters=x, fixed.parameters=pfixed,
                           temperatures=formated, 
                           integral=integral.Gompertz, M0=M0, 
                           hatchling.metric=c(Mean=39.33, SD=1.92))
plotR(resultNest_newp, ylim=c(0, 2), 
      xlim=c(23, 34), ylimH=c(0, 3), show.hist=TRUE)
compare_AIC(test4p=resultNest_4p_SSM, 
            test6p=resultNest_6p_SSM, 
            testAnchor=resultNest_newp)
            
############################################
# example with thermal reaction norm fitted from Weibull function
############################################

 x <- ChangeSSM(temperatures = (200:350)/10,
                parameters = resultNest_4p_SSM$par,
                initial.parameters = structure(c(73.4009010417375, 304.142079511996, 
                                                27.4671689276281), 
                                        .Names = c("k", "lambda", "scale")), 
                control=list(maxit=1000))
M0 <- 0.3470893 
pfixed <- c(rK=1.208968)
resultNest_3p_Weibull <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=M0, 
                         hatchling.metric=c(Mean=39.33, SD=1.92))
plotR(resultNest_3p_Weibull, ylim=c(0,6), col="Black")
compare_AIC(SSM=resultNest_4p_SSM, Weibull=resultNest_3p_Weibull)

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

x <- ChangeSSM(temperatures = (200:350)/10,
               parameters = resultNest_4p_SSM$par,
               initial.parameters = structure(c(3, 7, 11, 32), 
                               .Names = c("Scale", "sdL", "sdH", "Peak")), 
               control=list(maxit=1000))
M0 <- 0.3470893 
pfixed <- c(rK=1.208968)
resultNest_4p_normal <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=M0, 
                         hatchling.metric=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_SSM$par,
               initial.parameters = structure(c(3, 20, 40, 32), 
               .Names = c("Max", "LengthB", "LengthE", "Peak")), 
               control=list(maxit=1000))
M0 <- 0.3470893 
pfixed <- c(rK=1.208968)
resultNest_4p_trigo <- searchR(parameters=x$par, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=M0, 
                         hatchling.metric=c(Mean=39.33, SD=1.92))
                         
###############################################################
# example with thermal reaction norm fitted from Dallwitz model
###############################################################
# See: Dallwitz, M.J., Higgins, J.P., 1992. User’s guide to DEVAR. A computer 
# program for estimating development rate as a function of temperature. CSIRO Aust 
# Div Entomol Rep 2, 1-23.

# Note that Dallwitz model has many problems and I recommend to not use it:
# - The 3-parameters is too highly constraint
# - The 5 parameters produced infinite outputs for some sets of parameters that
#   can be generated while using delta method.

x <- c('Dallwitz_b1' = 4.8854060791241816, 
       'Dallwitz_b2' = 20.398366565842029, 
       'Dallwitz_b3' = 31.510995256647092)
M0 <- 0.3470893 
pfixed <- c(rK=1.208968)
resultNest_3p_Dallwitz <- searchR(parameters=x, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=M0, 
                         hatchling.metric=c(Mean=39.33, SD=1.92))
plotR(resultNest_3p_Dallwitz, ylim=c(0,6))

x <- c('Dallwitz_b1' = 4.9104386262684656, 
       'Dallwitz_b2' = 7.515425231891359, 
       'Dallwitz_b3' = 31.221784599026638, 
       'Dallwitz_b4' = 7.0354467023505682, 
       'Dallwitz_b5' = -1.5955717975708577)
pfixed <- c(rK=1.208968)
resultNest_5p_Dallwitz <- searchR(parameters=x, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=0.3470893, 
                         hatchling.metric=c(Mean=39.33, SD=1.92))
plotR(resultNest_5p_Dallwitz, ylim=c(0,3), scaleY=10000)

xp <- resultNest_6p_SSM$par
xp["Rho25"] <- 233
pfixed <- c(rK=1.208968)
resultNest_6p_SSM <- searchR(parameters=xp, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=0.3470893, 
                         hatchling.metric=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_SSM$par)
pfixed <- c(rK=1.208968)
resultNest_4p_SSM <- searchR(parameters=xp$par, fixed.parameters=pfixed, 
                         temperatures=formated, integral=integral.Gompertz, M0=0.3470893, 
                         hatchling.metric=c(Mean=39.33, SD=1.92))
plotR(resultNest_4p_SSM, ylim=c(0,6))

compare_AIC(Dallwitz3p=resultNest_3p_Dallwitz, Dallwitz5p=resultNest_5p_Dallwitz, 
             SSM=resultNest_4p_SSM, SSM=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, integral=integral.linear, M0=0, 
                         hatchling.metric=c(Mean=1, SD=NA))
 plotR(resultNest_PropDev_3p_Dallwitz, ylim=c(0, 1.5), curve="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, integral=integral.linear, M0=0, 
                         hatchling.metric=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), curve="ML")
 plotR(resultNest_PropDev_5p_Dallwitz, ylim=c(0, 1.5), curve="ML", new=FALSE, col="red")
 compare_AICc(Dallwitz3p=resultNest_PropDev_3p_Dallwitz, 
              Dallwitz5p=resultNest_PropDev_5p_Dallwitz)
              
###########################################################################
# Dalwitz model with proportion of development and fitted SD for final size
###########################################################################

x <- c('Dallwitz_b1' = 1.4886497996404355, 
       'Dallwitz_b2' = 10.898310418085916, 
       'Dallwitz_b3' = 31.263224721068056, 
       'Dallwitz_b4' = 6.1624623077734535, 
       'Dallwitz_b5' = -1.0027132357973265, 
       'SD' = 0.041829475961912894)
resultNest_PropDev_5p_Dallwitz <- searchR(parameters=x, fixed.parameters=NULL, 
                         temperatures=formated, integral=integral.linear, M0=0, 
                         hatchling.metric=c(Mean=1))
 plotR(resultNest_PropDev_5p_Dallwitz, ylim=c(0, 1.5), curve="ML")
 # Note that the standard error of the curve cannot be estimated with delta method. 
 # MCMC should be used
 plot(x=resultNest_PropDev_5p_Dallwitz, ylimS=c(0,1), xlim=c(0,60), series=2, 
         embryo.stages="Generic.ProportionDevelopment")
         
##############################################################################         
# Parameters Threshold_Low and Threshold_High are used to truncate growth rate
##############################################################################         

plotR(result=resultNest_PropDev_5p_Dallwitz, 
     fixed.parameters=c(Threshold_Low=26, 
                        Threshold_High=33), 
     ylim=c(0, 1.5), curve="ML")
}

Run the code above in your browser using DataLab