MH_MCMC_FPK {BBMV} | R Documentation |

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

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)

`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. |

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.

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. |

F. C. Boucher

## 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)

[Package *BBMV* version 2.1 Index]