Data engineering and BI courses are free!
Run the Metropolis-Hastings algorithm for Sexualisation Thermal Reaction Norm.
The number of iterations is n.iter+n.adapt+1 because the initial likelihood is also displayed.
I recommend that thin=1 because the method to estimate SE uses resampling.
If initial point is maximum likelihood, n.adapt = 0 is a good solution.
To get the SE of the point estimates from result_mcmc <- STRN_MHmcmc(result=try), use:
result_mcmc$SD
coda package is necessary for this function.
The dataSTRN is a named list with the following objects:
EmbryoGrowthTRN= result of searchR
tsd= result of tsd
sexed= vector with number of sexed embryos
males= vector with number of males (could be also females=)
Temperatures= a text of the temperatures name used as CTE
The Temperatures text for CTE can be:
TimeWeighted.temperature.mean
TSP.TimeWeighted.temperature.mean
TSP.MassWeighted.temperature.mean
TSP.STRNWeighted.temperature.mean
TSP.MassWeighted.STRNWeighted.temperature.mean
MiddleThird.TimeWeighted.temperature.mean
They are explained in the info.nests function.
This function is not still fully described as it has not been published still.
The parameters intermediate and filename are used to save intermediate results every 'intermediate' iterations (for example 1000). Results are saved in a file of name filename.
The parameter previous is used to indicate the list that has been save using the parameters intermediate and filename. It permits to continue a mcmc search.
These options are used to prevent the consequences of computer crash or if the run is very very long and processes at time limited.
STRN_MHmcmc(result = NULL, n.iter = 10000, parametersMCMC = NULL,
n.chains = 1, n.adapt = 0, thin = 1, trace = NULL,
batchSize = sqrt(n.iter), dataSTRN = NULL, adaptive = FALSE,
adaptive.lag = 500, adaptive.fun = function(x) { ifelse(x > 0.234,
1.3, 0.7) }, parallel = TRUE, intermediate = NULL,
filename = "intermediate.Rdata", previous = NULL)An object obtained after a STRN fit
Number of iterations for each step
A set of parameters used as initial point for searching with information on priors
Number of replicates
Number of iterations before to store outputs
Number of iterations between each stored output
True or False, shows progress
Number of observations to include in each batch fo SE estimation
A named list data used to estimate likelihoods (see further in description)
Should an adaptive process for SDProp be used
Lag to analyze the SDProp value in an adaptive content
Function used to change the SDProp
Should parallel computing for info.nests() be used
Period for saving intermediate result, NULL for no save
If intermediate is not NULL, save intermediate result in this file
Previous result to be continued. Can be the filename in which intermediate results are saved.
A list with resultMCMC being mcmc.list object, resultLnL being likelihoods and parametersMCMC being the parameters used
STRN_MHmcmc runs the Metropolis-Hastings algorithm for STRN (Bayesian MCMC)
# NOT RUN {
library(embryogrowth)
MedIncubation_Cc <- subset(DatabaseTSD, Species=="Caretta caretta" &
RMU=="Mediterranean" & Sexed!=0)
Med_Cc <- with(MedIncubation_Cc, tsd(males=Males, females=Females,
temperatures=Incubation.temperature, par=c(P=29.5, S=-0.01)))
plot(Med_Cc, xlim=c(25, 35))
males <- c(7, 0, 0, 0, 0, 5, 6, 3, 5, 3, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0)
names(males) <- rev(rev(names(resultNest_4p_SSM4p$data))[-(1:2)])
sexed <- rep(10, length(males))
names(sexed) <- rev(rev(names(resultNest_4p_SSM4p$data))[-(1:2)])
Initial_STRN <- resultNest_4p_SSM4p$par[c("DHA", "DHH", "T12H")]
Initial_STRN <- structure(c(3460.21379985491, 588.062535503578, 2347.70617453574),
.Names = c("DHA", "DHH", "T12H"))
fp <- c(Rho25=100)
fitSTRN <- STRN(Initial_STRN, EmbryoGrowthTRN=resultNest_4p_SSM4p, tsd=Med_Cc,
Sexed=sexed, Males=males,
fixed.parameters=fp,
Temperatures="TSP.MassWeighted.STRNWeighted.temperature.mean")
pMCMC <- TRN_MHmcmc_p(fitSTRN, accept=TRUE)
pMCMC[, "Density"] <- "dunif"
pMCMC[, "Prior2"] <- pMCMC[, "Max"]<- 10000
pMCMC[, "Prior1"] <- pMCMC[, "Min"] <- 1
outMCMC <- STRN_MHmcmc(result = fitSTRN, n.iter = 50000, parametersMCMC = pMCMC,
n.chains = 1, n.adapt = 0, thin = 1, trace = TRUE,
dataSTRN = list(EmbryoGrowthTRN = resultNest_4p_SSM4p,
Temperatures = "TSP.STRNWeighted.temperature.mean",
fixed.parameters=fitSTRN$fixed.parameters,
tsd=Med_Cc,
Sexed=sexed, Males=males),
adaptive = TRUE, adaptive.lag = 500,
intermediate = 1000,
filename = "intermediate_mcmcSTRN.Rdata")
plot(outMCMC, parameters=1)
plot(outMCMC, parameters=2)
plot(outMCMC, parameters=3)
1-rejectionRate(as.mcmc(x = outMCMC))
# Take care,you computer will be 100% busy as all cores will be used intensively
outMCMC_parallel <- parallel::mclapply(1:detectCores(), function(x) {
STRN_MHmcmc(result = fitSTRN, n.iter = 50000/detectCores(),
parametersMCMC = pMCMC,
n.chains = 1, n.adapt = 0, thin = 1, trace = TRUE,
dataSTRN = list(EmbryoGrowthTRN = resultNest_4p_SSM4p,
Temperatures = "TSP.STRNWeighted.temperature.mean",
fixed.parameters=fitSTRN$fixed.parameters,
tsd=Med_Cc,
Sexed=sexed, Males=males),
parallel=FALSE,
adaptive = TRUE, adaptive.lag = 500,
intermediate = NULL)
})
outMCMC_parallel_merge <- outMCMC_parallel[[1]]
for (i in 2:length(outMCMC_parallel)) {
outMCMC_parallel_merge <- merge(outMCMC_parallel_merge, outMCMC_parallel[[i]])
}
plot(outMCMC_parallel_merge, parameters=1)
plot(outMCMC_parallel_merge, parameters=2)
plot(outMCMC_parallel_merge, parameters=3)
plotR(parameters = fitSTRN$par, fixed.parameters=fitSTRN$fixed.parameters,
curves = "MCMC quantiles", ylim=c(0, 5), resultmcmc = outMCMC_parallel_merge,
ylab="Relative contribution to sexualisation", xlim=c(28, 29))
# }
Run the code above in your browser using DataLab