embryogrowth (version 6.1.1)

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

Description

Fit the parameters that best represent data. test can be a list with two elements Mean and SD and each element is a named vector with the nest name.

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,
  hessian = TRUE, 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)

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

hessian

If TRUE, the hessian matrix is estimated and the SE of parameters estimated.

control

List for control parameters for optimx

Value

• A result object

Details

searchR fits the parameters that best represent nest incubation data.

Examples

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(106.59891311201, 614.181133951497, 306.267053513175,
120.327257089974), .Names = c("DHA", "DHH", "T12H", "Rho25"))
# pfixed <- c(K=82.33) or rK=82.33/39.33
pfixed <- c(rK=2.093313)
resultNest_4p <- 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)
plot(resultNest_4p, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,45), series=1)
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)
resultNest_6p <- searchR(parameters=x, fixed.parameters=pfixed,
	temperatures=formated, derivate=dydt.Gompertz, M0=1.7,
	test=c(Mean=39.33, SD=1.92))
data(resultNest_6p)
pMCMC <- TRN_MHmcmc_p(resultNest_6p, accept=TRUE)
# Take care, it can be very long, sometimes several days
result_mcmc_6p <- GRTRN_MHmcmc(result=resultNest_6p,
	parametersMCMC=pMCMC, n.iter=10000, n.chains = 1, n.adapt = 0,
	thin=1, trace=TRUE)
data(result_mcmc_6p)
# compare_AIC() is a function from the package "HelpersMG"
compare_AIC(test1=resultNest_4p, test2=resultNest_6p)
############ with new parametrization
data(resultNest_4p)
x0 <- resultNest_4p$par
t <- hist(resultNest_4p, 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_hist(resultNest_newp, ylim=c(0,0.3), xlimR=c(23, 34), ylimH=c(0, 0.3))
compare_AIC(test4p=resultNest_4p,
            test6p=resultNest_6p,
            testAnchor=resultNest_newp)