# NOT RUN {
library("embryogrowth")
packageVersion("embryogrowth")
data(nest)
formated <- FormatNests(nest)
# The initial parameters value can be:
# "T12H", "DHA", "DHH", "Rho25"
# Or
# "T12L", "DT", "DHA", "DHH", "DHL", "Rho25"
x <- structure(c(115.758929130522, 428.649022170996, 503.687251738993,
12.2621455821612, 306.308841227278, 116.35048615105), .Names = c("DHA",
"DHH", "DHL", "DT", "T12L", "Rho25"))
# or
x <- structure(c(118.431040984352, 498.205702157603, 306.056280989839,
118.189669472381), .Names = c("DHA", "DHH", "T12H", "Rho25"))
# pfixed <- c(K=82.33) or rK=82.33/39.33
pfixed <- c(rK=2.093313)
################################################################################
#
# The values of rK=2.093313 and M0=1.7 were used in
# 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.
#
# Based on recent analysis on table of development for both Emys orbicularis and
# Caretta caretta, best value for pfixed should be 1.209 and M0 should be 0.34.
# 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.
#
# See the example in the stages datasets
#
################################################################################
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)
par(mar=c(4, 4, 1, 1))
plot(resultNest_4p_SSM4p$data[[1]][, 1]/60/24,resultNest_4p_SSM4p$data[[1]][, 2], bty="n", las=1,
xlab="Days of incubation", ylab="Temperatures in <U+00B0>C",
type="l", xlim=c(0,70),ylim=c(20, 35))
for (i in 2:resultNest_4p_SSM4p$data$IndiceT[3]) {
par(new=TRUE)
plot(resultNest_4p_SSM4p$data[[i]][, 1]/60/24,resultNest_4p_SSM4p$data[[i]][, 2],
bty="n", las=1, xlab="", ylab="", type="l", xlim=c(0,70),ylim=c(20, 35), axes = FALSE)
}
par(mar=c(4, 4, 1, 1))
pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM4p, accept=TRUE)
# Take care, it can be very long, sometimes several days
resultNest_mcmc_4p_SSM4p <- GRTRN_MHmcmc(result=resultNest_4p_SSM4p,
parametersMCMC=pMCMC, n.iter=10000, n.chains = 1, n.adapt = 0,
thin=1, trace=TRUE)
data(resultNest_mcmc_4p_SSM4p)
out <- as.mcmc(resultNest_mcmc_4p_SSM4p)
# This out obtained after as.mcmc can be used with coda package
# plot() can use the direct output of GRTRN_MHmcmc() function.
plot(resultNest_mcmc_4p_SSM4p, parameters=1, xlim=c(0,550))
plot(resultNest_mcmc_4p_SSM4p, parameters=3, xlim=c(290,320))
# But rather than to use the SD for each parameter independantly, it is
# more logical to estimate the distribution of the curves
new_result <- ChangeSSM(resultmcmc = resultNest_mcmc_4p_SSM4p, result = resultNest_4p_SSM4p,
temperatures = seq(from = 20, to = 35, by = 0.1),
initial.parameters = NULL)
par(mar=c(4, 4, 1, 5)+0.4)
plotR(result = resultNest_4p_SSM4p, parameters = new_result$par,
ylabH = "Temperatures\ndensity", ylimH=c(0, 0.3), atH=c(0, 0.1, 0.2),
ylim=c(0, 3), show.hist=TRUE)
# Beautiful density plots
plotR(result = resultNest_4p_SSM4p,
resultmcmc=resultNest_mcmc_4p_SSM4p,
curves = "MCMC quantiles", show.density=TRUE)
plotR(resultNest_6p_SSM6p, resultmcmc=resultNest_mcmc_6p_SSM6p,
ylim=c(0, 4), show.density=TRUE, show.hist=TRUE,
curves = "MCMC quantiles",
ylimH=c(0,0.5), atH=c(0, 0.1, 0.2))
# }
Run the code above in your browser using DataLab