embryogrowth (version 5.0)

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

Description

Fit the parameters that best represent data.

Usage

searchR(parameters = stop("Initial set of parameters must be provided"),
  fixed.parameters = NULL,
  temperatures = stop("Temperature data must be provided !"),
  derivate = dydt.Gompertz, test = c(Mean = 39.33, SD = 1.92), M0 = 1.7,
  method = "BFGS", maxiter = 200, saveAtMaxiter = FALSE,
  fileName = "intermediate", weight = NULL, hessian = TRUE,
  parallel = (.Platform$OS.type == "unix"))

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

Mean and SD of size of hatchlings

M0

Measure of hatchling size or mass proxi at laying date

method

Method uses for searching. Can be "x" then it uses "Nelder-Mead" and "BFGS", or any method from optim function

maxiter

After maxiter iteration, the value of parameters is displayed but it continues if convergence is not acheived

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.

parallel

If true, try to use several cores using parallel computing. Must be FALSE in Windows.

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(118.768297442004, 475.750095909406, 306.243694918151,
116.055824800264), .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), method = "BFGS", maxiter = 200)
data(resultNest_4p)
plot(resultNest_4p, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,45), series=1)
x <- structure(c(115.758929130522, 428.649022170996, 503.687251738993,
12.2621455821612, 306.308841227278, 116.35048615105), .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), method = "BFGS", maxiter = 200)
data(resultNest_6p)
pMCMC <- embryogrowth_MHmcmc_p(resultNest_6p, accept=TRUE)
# Take care, it can be very long, sometimes several days
result_mcmc_6p <- embryogrowth_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 "phenology"
outputAIC<-compare_AIC(list(test1=resultNest_4p, test2=resultNest_6p))
############ with new parametrization
data(resultNest_4p)
x0 <- resultNest_4p$par
t <- hist(resultNest_4p, plot=FALSE)
temperatures <- seq(from=floor(range(t$temperatures)[1]+273.15-1),
to=floor(range(t$temperatures)[2]+273.15+1), length.out=7)
newx <- embryogrowth:::.SSM(temperatures, x0)[[1]]*1E5
names(newx) <- temperatures
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), method = "BFGS", maxiter = 200)
plotR_hist(resultNest_newp, ylim=c(0,0.3), xlimR=c(23, 34), ylimH=c(0, 0.3))
outputAIC<-compare_AIC(list(test4p=resultNest_4p, test6p=resultNest_6p,
 testAnchor=resultNest_newp))