MH_MCMC_FPK_multiclades: MCMC estimation on multiple clades

In BBMV: Models for Continuous Traits Evolving in Macroevolutionary Landscapes of any Shape

Description

This function estimates parameter of the FPK model on multiple clades at once, making the assumption that they share the same macroevolutionary landscape but have different rates of evolution.

Usage

MH_MCMC_FPK_multiclades(trees, traits, bounds, Nsteps = 5e+05, record_every = 100, 
    plot_every = 500, Npts = 50, pars_init = NULL, prob_update = NULL, verbose = TRUE, 
    plot = TRUE, save_to = "MCMC_FPK_test.Rdata", save_every = 10000, type_priors = NULL, 
    shape_priors = NULL, proposal_type = "Normal", proposal_sensitivity = NULL, 
    prior.only = F, burnin.plot = 0.1)

Arguments

trees	A list of phylogenetic trees in 'phylo' format, one per clade.
traits	A list of trait vectors for species in each clade. Should be in the same order as trees.
bounds	The two bounds that constrain trait values when fitting the BBMV model. Specified by a numeric vector containing the minimum and maximum bound of the trait interval as the first and second element, respectively.
Nsteps	The number of steps in the MCMC chain.
record_every	How often to record a generation in the chain.
plot_every	How often to plot the trace of the chain.
Npts	The number of points used in the discretization procedure.
pars_init	The initial parameter values.
prob_update	A numeric vector with the relative probability of update of each parameter of the model.
verbose	If TRUE, prints generations to the screen.
plot	If TRUE, plots the trace of parameter values along iterations during the MCMC run.
save_to	The directory in which the chain should be saved.
save_every	How often to save the chain.
type_priors	The type of priors used, can be either normal (preferred) or uniform for log(sig2/2), a, b and c, ; and can only be discrete uniform for x0.
shape_priors	A list that gives the shape for each prior. (mean,sd) for normal priors and (min,max) for continuous uniform priors. The shape is not specified for the root prior, since it is fixed to be discrete uniform on the grid.
proposal_type	The type of proposal function, only uniform is available so far.
proposal_sensitivity	The width of the uniform proposal. The entire value for x0 gives how many steps at a time can be travelled on the trait grid (better to keep it to 1)
prior.only	If TRUE, only the prior is explored but the likelihood is ignored. Default to false for estimation of the posterior.
burnin.plot	The frequency of samples that should be discarded as burnin in trace plots.

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 n parameter, n being the number of clades studied, are diffusion coefficients (log(sig^2/2)), not evolutionary rates (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_clade_i	The evolutionary rate, one column per clade.
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.
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

## Not run: 
x=seq(from=-1.5,to=1.5,length.out=100)
bounds=c(min(x),max(x)) # the bounds we use for simulating
a=8 ; b=-4 ; c=1
V6=a*x^4+b*(x^2)+c*x
step_size=(max(bounds)-min(bounds))/(100-1)
V6_norm=exp(-V6)/sum(exp(-V6)*step_size)
par(mfrow=c(1,1))
plot(V6_norm,type='l')

# Now we simulate a tree and a continuous trait for 3 independent clades. 
tree=sim.bdtree(stop='taxa',n=25) # tree with few tips for quick tests
tree$edge.length=100*tree$edge.length/max(branching.times(tree))
TRAIT= Sim_FPK(tree,x0=0.5,V=V6,sigma=1,bounds=bounds)
tree1=tree ; TRAIT1=TRAIT

tree=sim.bdtree(stop='taxa',n=25)
tree$edge.length=100*tree$edge.length/max(branching.times(tree))
TRAIT= Sim_FPK(tree,x0=0.5,V=V6,sigma=0.5,bounds=bounds) 
tree2=tree ; TRAIT2=TRAIT

tree=sim.bdtree(stop='taxa',n=25)
tree$edge.length=100*tree$edge.length/max(branching.times(tree))
TRAIT= Sim_FPK(tree,x0=0.5,V=V6,sigma=0.1,bounds=bounds) 
tree3=tree ; TRAIT3=TRAIT
rm(tree) ; rm(TRAIT)

TREES=list(tree1,tree2,tree3)
TRAITS=list(TRAIT1,TRAIT2,TRAIT3)

# Fit the FPK model using ML: 
# In all clades the macroevolutionary landscape is the same 
#but they have different evolutionary rates
testbFPK4=lnl_FPK_multiclades_same_V_different_sig2(trees=TREES,
  traits=TRAITS,a=NULL,b=NULL,c=NULL,Npts=50)
fitbFPK4=find.mle_FPK_multiple_clades_same_V_different_sig2(model=testbFPK4,
  method='Nelder-Mead',init.optim=NULL)

# And now MCMC run
mcmc1=MH_MCMC_FPK_multiclades(trees=TREES,traits=TRAITS,
  bounds=fitmFPK4_SE$fits$fit_clade_1$par_fixed$bounds,Nsteps=10000,record_every=100,
  plot_every=200,Npts=25,pars_init=NULL,prob_update=NULL,verbose=TRUE,plot=TRUE,
  save_to='MCMC_FPK_test.Rdata',save_every=1000,type_priors=NULL,shape_priors=NULL,
  proposal_type='Normal',proposal_sensitivity=NULL,prior.only=F,burnin.plot=0.1)

get.landscape.FPK.MCMC(chain=mcmc1,bounds,Npts=100,burnin=0.1,probs.CI=c(0.25,0.75),
  COLOR_MEDIAN='red',COLOR_FILL='red',transparency=0.3,main='Macroevolutionary landscapes MCMC',
  ylab='N.exp(-V)',xlab='Trait',xlim=NULL,ylim=c(0,2))
  
## End(Not run)

BBMV documentation built on May 1, 2019, 10:26 p.m.