STRN_MHmcmc: Metropolis-Hastings algorithm for Sexualisation Thermal...

In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth

Metropolis-Hastings algorithm for Sexualisation Thermal Reaction Norm

Description

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 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.
If fill is NA, it will use the stored fill value in result.

Usage

STRN_MHmcmc(
  result = NULL,
  n.iter = 10000,
  parametersMCMC = NULL,
  n.chains = 1,
  n.adapt = 0,
  thin = 1,
  trace = NULL,
  traceML = FALSE,
  batchSize = sqrt(n.iter),
  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,
  fill = NA
)

Arguments

result	An object obtained after a STRN fit
n.iter	Number of iterations for each step
parametersMCMC	A set of parameters used as initial point for searching with information on priors
n.chains	Number of replicates
n.adapt	Number of iterations before to store outputs
thin	Number of iterations between each stored output
trace	TRUE or FALSE or period, shows progress
traceML	TRUE or FALSE to show ML
batchSize	Number of observations to include in each batch fo SE estimation
adaptive	Should an adaptive process for SDProp be used
adaptive.lag	Lag to analyze the SDProp value in an adaptive content
adaptive.fun	Function used to change the SDProp
parallel	Should parallel computing for info.nests() be used
intermediate	Period for saving intermediate result, NULL for no save
filename	If intermediate is not NULL, save intermediate result in this file
previous	Previous result to be continued. Can be the filename in which intermediate results are saved.
fill	Parameters to be sent to STRN().

Details

STRN_MHmcmc runs the Metropolis-Hastings algorithm for STRN (Bayesian MCMC)

Value

A list with resultMCMC being mcmc.list object, resultLnL being likelihoods and parametersMCMC being the parameters used

Author(s)

Marc Girondot [e-mail]

Examples

## Not run: 
library(embryogrowth)
MedIncubation_Cc <- subset(DatabaseTSD, Species=="Caretta caretta" & 
RMU=="Mediterranean" & Sexed!=0)
Med_Cc <- tsd(males=MedIncubation_Cc$Males, 
             females=MedIncubation_Cc$Females, 
             temperatures=MedIncubation_Cc$Incubation.temperature, 
             par=c(P=29.5, S=-0.1))
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_SSM$data))[-(1:2)])
sexed <- rep(10, length(males))
names(sexed) <- rev(rev(names(resultNest_4p_SSM$data))[-(1:2)])
Initial_STRN <- resultNest_4p_SSM$par[c("DHA", "DHH", "T12H")]
Initial_STRN <- structure(c(582.567096666926, 2194.0806711639, 3475.28414940385), 
                          .Names = c("DHA", "DHH", "T12H"))
fp <- c(Rho25=100)
fitSTRN <- STRN(Initial_STRN=Initial_STRN, 
                EmbryoGrowthTRN=resultNest_4p_SSM, 
                tsd=Med_Cc, 
                embryo.stages="Caretta caretta.SCL", 
                Sexed=sexed, Males=males, 
                fixed.parameters=fp, 
                sexratio="TSP.GrowthWeighted.STRNWeighted.sexratio")
                
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 = 10000, parametersMCMC = pMCMC,
                n.chains = 1, n.adapt = 0, thin = 1, trace = TRUE,
                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))


## End(Not run)

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