MH_MCMC_FPK: MCMC estimation

Description Usage Arguments Details Value Author(s) Examples

View source: R/MH_MCMC_FPK.R

Description

The function estimates the parameters of the BBM+V model using an MCMC chain with the Metropolis Hastings algorithm and a Gibbs sampler.

Usage

1
2
3
4
5
6
7
MH_MCMC_FPK(tree, trait, bounds, Nsteps = 5e+05, record_every = 100, plot_every = 500, 
  Npts = 50, pars_init = c(0, 0, 0, 0, 25), prob_update = c(0.2, 0.2, 0.2, 0.2, 0.2), 
  verbose = TRUE, plot = TRUE, save_to = "MCMC_FPK_test.Rdata", save_every = 10000, 
  type_priors = c(rep("Normal", 4), "Uniform"), 
  shape_priors = list(c(0, 10), c(0, 10), c(0, 10), c(0, 10), NA), 
  proposal_type = "Uniform", proposal_sensitivity = c(0.1, 0.1, 0.1, 0.1, 1), 
  prior.only = F, burnin.plot = 0.1)

Arguments

tree

A phylogenetic tree in phylo format.

trait

A named vector of trait values for the tips of the tree. It should match tip labels in the phylogeny.

bounds

A vector with two elements giving the bounds on the trait interval.

Nsteps

The number of generations in the MCMC chain.

record_every

The frequency used for sampling the MCMC chain.

plot_every

The frequency at which the chain is plotted (if plot=TRUE).

Npts

The number of points on the grid between the bounds.

pars_init

A vector giving the initial parameters for starting the algorithm, which correspond to the following: c(log(sig^2/2),a,b,c,x0).

prob_update

A vector giving the relative frequencies of update of the different parameters of the model.

verbose

If TRUE, will print some generations of the chain to the screen.

plot

If TRUE, the chain is plotted from time to time.

save_to

The path to the file where the chain is saved (can be useful in case the chain crashes).

save_every

Sets how often the chain is saved.

type_priors

A character vector specifying the type of priors used. Either 'Uniform' or 'Normal'. See Details.

shape_priors

A list that gives the shape for each prior. See Details.

proposal_type

The type of proposal function, only 'Uniform' is available (the default).

proposal_sensitivity

A numeric vector specifying the width of the uniform proposal for each parameter. See Details.

prior.only

Default to FALSE for estimation of the posterior. If TRUE, the likelihood is not evaluated: this is mostly useful for internal test of the Gibbs sampler.

burnin.plot

The percentage of samples considered as burnin and thus not shown on the trace plot that the function produces.

Details

When specifying intial parameters yourself, be careful since x0 is actually the index of the point on the grid (between 1 and Npts_int), not the actual root value. Also the first parameter is the diffusion coefficient (log(sig^2/2)), not the evolutionary rate (sig^2). Finally, be careful that the bounds you propose must contain all trait values in you dataset.

Priors can be either 'Normal' (preferred) or 'Uniform' for log(sig^2/2), a, b and c. The only option for x0 is a discrete uniform prior, specified by 'Uniform'.

Each element of the shape_priors list should be a vector giving c(mean,sd) for normal priors and c(min,max) for continuous uniform priors. The shape is not specified for the root prior (it is set as 'NA' by default), since it is fixed to be discrete uniform on the grid.

Elements of the proposal_sensitivity vector can be any positive number for continuously varying parameters: c(log(sig^2/2),a,b,c). Default values should often be a good start. Only integer numbers are possible for x0 and give how many steps at a time can be travelled on the trait grid when updating these parameters. It is recommended to keep it to 1, as it is by default.

Value

A matrix of numeric values giving values of all parameters, the likelihood, prior and posterior at each generation sampled in the MCMC chain (one row per sample taken). The matrix has the following columns:

step

The number of the generation sampled.

sigsq

The evolutionary rate.

a

The coefficient of the x^4 term of the evolutionary potential.

b

The coefficient of the x^2 term of the evolutionary potential.

c

The coefficient of the x term of the evolutionary potential.

root

The value of the trait at the root of the tree.

lnprior

The logarithm of the prior.

lnlik

The logarithm of the likelihood.

quasi-lnpost

The logarithm of the (unnormalized) posterior.

Acceptance

Whether the proposed MCMC move was accepted (1) or not (0).

Par_updated

Which parameter was updated in this generation.

Author(s)

F. C. Boucher

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
## Not run: 
# Simulate data: tree + continuous trait
library(geiger)
tree=sim.bdtree(stop='taxa',n=10) # tree with few tips for quick tests
tree$edge.length=100*tree$edge.length/max(branching.times(tree)) # rescale the tree 
# Simulate trait evolving on a macroevolutionary landscape with two peaks of equal heights
x=seq(from=-1.5,to=1.5,length.out=100)
bounds=c(min(x),max(x)) # the bounds we use for simulating: for technical purposes only
V6=10*(x^4-0.5*(x^2)+0.*x) # this is the evolutionary potential: it has two wells
TRAIT= Sim_FPK(tree,x0=0,V=V6,sigma=10,bounds=c(-5, 5)) 
# Run a MCMC chain to fit the FPK model
MCMC=MH_MCMC_FPK(tree,trait=TRAIT,bounds=c(5,5),Nsteps=10000,record_every=100,
  plot_every=100,Npts=20,pars_init=c(0,-4,-4,0,1),prob_update=c(0.2,0.25,0.25,0.25,0.05),
  verbose=TRUE,plot=TRUE,save_to='MCMC_FPK_test.Rdata',save_every=100,
  type_priors=c(rep('Normal',4),'Uniform'),
  shape_priors=list(c(0,10),c(0,10),c(0,10),c(0,10),NA),proposal_type='Uniform',
  proposal_sensitivity=c(0.1,0.1,0.1,0.1,1),prior.only=F)

## End(Not run)

BBMV documentation built on April 30, 2018, 9:03 a.m.