HatchingSuccess.fit: Fit a hatching success model to data using maximum likelihood
In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth

HatchingSuccess.fit

R Documentation

Fit a hatching success model to data using maximum likelihood

Description

Set of functions to study the hatching success.
The first version of the model was published in:
Laloë, J.-O., Monsinjon, J., Gaspar, C., Touron, M., Genet, Q., Stubbs, J., Girondot, M. & Hays, G.C. (2020) Production of male hatchlings at a remote South Pacific green sea turtle rookery: conservation implications in a female-dominated world. Marine Biology, 167, 70.
The version available here is enhanced by using a double flexit model rather than a double logistic model. The flexit model is described here:
Abreu-Grobois, F.A., Morales-Mérida, B.A., Hart, C.E., Guillon, J.-M., Godfrey, M.H., Navarro, E. & Girondot, M. (2020) Recent advances on the estimation of the thermal reaction norm for sex ratios. PeerJ, 8, e8451.

Usage

HatchingSuccess.fit(
  par = NULL,
  data = stop("data must be provided"),
  fixed.parameters = NULL,
  column.Incubation.temperature = "Incubation.temperature",
  column.Hatched = "Hatched",
  column.NotHatched = "NotHatched",
  hessian = TRUE
)

Arguments

`par`	A set of parameters.
`data`	A dataset in a data.frame with a least three columns: Incubation.temperature, Hatched and NotHatched
`fixed.parameters`	A set of parameters that must not be fitted.
`column.Incubation.temperature`	Name of the column with incubation temperatures
`column.Hatched`	Name of the column with hatched number
`column.NotHatched`	Name of the column with not hatched number
`hessian`	Should Hessian matrix be estimated?

Details

HatchingSuccess.fit fits a hatching success model to data

Value

Return a object of class HatchingSuccess

Author(s)

Marc Girondot

Examples

## Not run: 
library(embryogrowth)
totalIncubation_Cc <- subset(DatabaseTSD, 
                             Species=="Caretta caretta" & 
                               Note != "Sinusoidal pattern" & 
                               !is.na(Total) & Total != 0)

par <- c(S.low=0.5, S.high=0.3, 
         P.low=25, deltaP=10, MaxHS=0.8)
         
g.logistic <- HatchingSuccess.fit(par=par, data=totalIncubation_Cc)
         
HatchingSuccess.lnL(par=g.logistic$par, data=totalIncubation_Cc)

plot(g.logistic)

par <- c(S.low=0.5, S.high=0.3, 
         P.low=25, deltaP=10, 
         K1.low=1, K2.low=-1, K1.high=1, K2.high=-1, 
         MaxHS=0.8)

g.flexit <- HatchingSuccess.fit(par=par, data=totalIncubation_Cc)

HatchingSuccess.lnL(par=g.flexit$par, data=totalIncubation_Cc)

compare_AICc(logistic=g.logistic, flexit=g.flexit)
plot(x=g.logistic, what = c("observations", "ML", "CI"), replicates=10000)

pMCMC <- HatchingSuccess.MHmcmc_p(result = g.logistic, accept = TRUE)
MCMC <- HatchingSuccess.MHmcmc(result = g.logistic, parametersMCMC = pMCMC,
                            n.iter = 100000, 
                           adaptive = TRUE)

plot(MCMC, parameters = "S.low")
plot(MCMC, parameters = "S.high")
plot(MCMC, parameters = "P.low")
plot(MCMC, parameters = "deltaP")
plot(MCMC, parameters = "MaxHS")

plot(x=g.logistic, what = c("observations", "ML", "CI"), 
        replicates=10000, resultmcmc=MCMC)

####### Exemple with Chelonia mydas
totalIncubation_Cm <- subset(DatabaseTSD, 
Species=="Chelonia mydas" & 
  Note != "Sinusoidal pattern" & 
  !is.na(Total) & Total != 0 & 
  !is.na(NotHatched) & !is.na(Hatched))
  
  totalIncubation_Cm$NotHatched <- totalIncubation_Cm$NotHatched + 
  ifelse(!is.na(totalIncubation_Cm$Undeveloped), totalIncubation_Cm$Undeveloped, 0)
  
  
  plot(x=totalIncubation_Cm$Incubation.temperature, 
  y=totalIncubation_Cm$Hatched/totalIncubation_Cm$Total, bty="n", las=1, 
  xlab="Constant incubation temperature", ylab="Proportion of hatching")

par <- c(S.low=0.5, S.high=0.3, 
         P.low=25, deltaP=10, MaxHS=0.8)

g.logistic <- HatchingSuccess.fit(par=par, data=totalIncubation_Cm)
plot(g.logistic)

pMCMC <- HatchingSuccess.MHmcmc_p(g.logistic, accept=TRUE)
mcmc <- HatchingSuccess.MHmcmc(result=g.logistic, parameters = pMCMC, 
                  adaptive=TRUE, n.iter=100000, trace=1000)
par <- as.parameters(mcmc)
par <- as.parameters(mcmc, index="median")

plot(mcmc, parameters=c("P.low"))
plot(mcmc, parameters=c("deltaP"))
plot(mcmc, parameters=c("S.low"))
plot(mcmc, parameters=c("S.high"))
plot(mcmc, parameters=c("MaxHS"))

plot(g.logistic, resultmcmc=mcmc, what = c("observations", "CI"))

## End(Not run)

embryogrowth documentation built on Oct. 24, 2023, 5:07 p.m.