GRTRN_MHmcmc: Metropolis-Hastings algorithm for Embryo Growth Rate Thermal...

Description Usage Arguments Details Value Author(s) Examples

View source: R/GRTRN_MHmcmc.R

Description

Run the Metropolis-Hastings algorithm for data.
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 <- GRTRN_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 named 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 with user limited time.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
GRTRN_MHmcmc(
  result = NULL,
  n.iter = 10000,
  parametersMCMC = NULL,
  n.chains = 1,
  n.adapt = 0,
  thin = 1,
  trace = NULL,
  traceML = FALSE,
  parallel = TRUE,
  adaptive = FALSE,
  adaptive.lag = 500,
  adaptive.fun = function(x) {     ifelse(x > 0.234, 1.3, 0.7) },
  intermediate = NULL,
  filename = "intermediate.Rdata",
  previous = NULL
)

Arguments

result

An object obtained after a SearchR 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

parallel

If true, try to use several cores using parallel computing

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

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.

Details

GRTRN_MHmcmc runs the Metropolis-Hastings algorithm for data (Bayesian MCMC)

Value

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

Author(s)

Marc Girondot

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
## Not run: 
library(embryogrowth)
data(nest)
formated <- FormatNests(nest)
# The initial parameters value can be:
# "T12H", "DHA",  "DHH", "Rho25"
# Or
# "T12L", "T12H", "DHA",  "DHH", "DHL", "Rho25"
############################################################################
pfixed <- c(rK=1.208968)
M0 = 0.3470893 
############################################################################
# 4 parameters
############################################################################
x c('DHA' = 109.31113503282113, 'DHH' = 617.80695919563857, 
    'T12H' = 306.38890489505093, 'Rho25' = 229.37265815800225)
    
resultNest_4p_SSM <- searchR(parameters=x, fixed.parameters=pfixed, 
	temperatures=formated, integral=integral.Gompertz, M0=M0, 
	hatchling.metric=c(Mean=39.33, SD=1.92))
data(resultNest_4p_SSM)
plot(resultNest_4p_SSM, xlim=c(0,70), ylimT=c(22, 32), ylimS=c(0,45), series=1, 
embryo.stages="Caretta caretta.SCL")
############################################################################
pMCMC <- TRN_MHmcmc_p(resultNest_4p_SSM, accept=TRUE)
# Take care, it can be very long; several days
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(result=resultNest_4p_SSM, 
   adaptive = TRUE,
		parametersMCMC=pMCMC, n.iter=10000, n.chains = 1,  
		n.adapt = 0, thin=1, trace=TRUE)
# The SDProp in pMCMC at the beginning of the Markov chain can 
# be considered as non-optimal. However, the values at the end 
# of the Markoc chain are better due to the use of the option
# adaptive = TRUE. Then a good strategy is to run again the MCMC
# with this final set:
pMCMC[, "SDProp"] <- resultNest_mcmc_4p_SSM$parametersMCMC$SDProp.end
# Also, take the set of parameters fitted as maximim likelihood
# as initial set of value
pMCMC[, "Init"] <- as.parameters(resultNest_mcmc_4p_SSM)
# Then I run again the MCMC; it will ensure to get the optimal distribution
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(result=resultNest_4p_SSM, 
   adaptive = TRUE,
		parametersMCMC=pMCMC, n.iter=10000, n.chains = 1,  
		n.adapt = 0, thin=1, trace=TRUE)
		
data(resultNest_mcmc_4p_SSM)
out <- as.mcmc(resultNest_mcmc_4p_SSM)
# This out can be used with coda package
# Test for stationarity and length of chain
require(coda)
heidel.diag(out)
raftery.diag(out)
# plot() can use the direct output of GRTRN_MHmcmc() function.
plot(resultNest_mcmc_4p_SSM, parameters=1, xlim=c(0,550))
plot(resultNest_mcmc_4p_SSM, parameters=3, xlim=c(290,320))
# summary() permits to get rapidly the standard errors for parameters
# They are store in the result also.
se <- result_mcmc_4p_SSM$SD
# the confidence interval is better estimated by:
apply(out[[1]], 2, quantile, probs=c(0.025, 0.975))
# The use of the intermediate method is as followed;
# Here the total mcmc iteration is 10000, but every 1000, intermediate
# results are saved in file intermediate1000.Rdata:
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(result=resultNest_4p_SSM, 
parametersMCMC=pMCMC, n.iter=10000, n.chains = 1,  
n.adapt = 0, thin=1, trace=TRUE, 
intermediate=1000, filename="intermediate1000.Rdata")
# If run has been stopped for any reason, it can be resumed with:
resultNest_mcmc_4p_SSM <- GRTRN_MHmcmc(previous="intermediate1000.Rdata")
# Example to use of the epsilon parameter to get confidence level
resultNest_4p_epsilon <- resultNest_4p
resultNest_4p_epsilon$fixed.parameters <- c(resultNest_4p_epsilon$par, 
                   resultNest_4p_epsilon$fixed.parameters)
resultNest_4p_epsilon$par <- c(epsilon = 0)
pMCMC <- TRN_MHmcmc_p(resultNest_4p_epsilon, accept = TRUE)
resultNest_mcmc_4p_epsilon <- GRTRN_MHmcmc(result = resultNest_4p_epsilon, 
           n.iter = 10000, parametersMCMC = pMCMC,
           n.chains = 1, n.adapt = 0, thin = 1, trace = TRUE, parallel = TRUE)
data(resultNest_mcmc_4p_epsilon)
plot(resultNest_mcmc_4p_epsilon, parameters="epsilon", xlim=c(-11, 11), las=1)
plotR(resultNest_4p_epsilon, SE=c(epsilon = unname(resultNest_mcmc_4p_epsilon$SD)), 
           ylim=c(0, 3), las=1)

## End(Not run)

embryogrowth documentation built on June 4, 2021, 9:05 a.m.