################################################################################ # # RevBayes Example: Bayesian inference of diversification rates under a # constant-rate pure-birth model # # # authors: Sebastian Hoehna and Tracy A. Heath # ################################################################################ ####################### # Reading in the Data # ####################### ###### This just defines a single model for all sites ####### ### Read in sequence data for both genes T <- readTrees("data/primates_tree.nex")[1] # Get some useful variables from the data. We need these later on. taxa <- T.taxa() # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() ###################### # pure-birth model # ###################### ### the birth rate is a stochastic random variable drawn from a uniform prior birth_rate ~ dnUniform(0,100.0) ### MCMC samples this variable using a scale proposal moves.append( mvScale(birth_rate,lambda=1.0,tune=true,weight=3.0) ) ### rho is the probability of sampling species at the present ### fix this to 233/367, since there are ~367 described species of primates ### and we have sampled 233 rho <- T.ntips()/367 ### the BDP is conditioned on the root time, we can get this value from our tree and set a constant node root_time <- T.rootAge() ### the time tree is a stochastic node modeled by the constant rate birth-death process (dnBDP) ### by setting mu to the constant value 0.0, we are specifying a pure-birth process timetree ~ dnBDP(lambda=birth_rate, mu=0.0, rho=rho, rootAge=root_time, samplingStrategy="uniform", condition="survival", taxa=taxa) ### clamp the model with the "observed" tree timetree.clamp(T) ############# # The Model # ############# ### workspace model wrapper ### mymodel = model(birth_rate) ### set up the monitors that will output parameter values to file and screen monitors.append( mnModel(filename="output/primates_Yule.log",printgen=10, separator = TAB) ) monitors.append( mnScreen(printgen=1000, birth_rate) ) ################ # The Analysis # ################ ### workspace mcmc ### mymcmc = mcmc(mymodel, monitors, moves, nruns=2, combine="mixed") ### run the MCMC ### mymcmc.run(generations=50000, tuningInterval=200) ## quit ## q()