################################################################################ # # RevBayes Example: Inferring a Phylogeny of Fossil Bears Using the Mk Model # with a discretized Beta distribution for the stationary frequencies. # # This file: Runs the full MCMC ... # # authors: April M. Wright, Michael J. Landis, Sebastian Höhna # ################################################################################ ####################### # Reading in the Data # ####################### # Import the morphological character matrix # # this file contains only the taxa for which morphological characters are available # morpho <- readDiscreteCharacterData("data/bears.nex") ## helpers taxa <- morpho.names() num_taxa <- taxa.size() num_branches <- 2 * num_taxa - 2 # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() ############## # Tree model # ############## # Set up branch length hyperprior with a move br_len_lambda ~ dnExp(0.2) moves.append( mvScale(br_len_lambda, weight=2) ) # Define the tree parameter. # First, we generate the topology. # We assume a uniform prior on topology. phylogeny ~ dnUniformTopologyBranchLength(taxa, branchLengthDistribution=dnExponential(br_len_lambda)) # compute the tree length from the phylogeny tree_length := phylogeny.treeLength() moves.append( mvNNI(phylogeny, weight=num_branches/2.0) ) moves.append( mvSPR(phylogeny, weight=num_branches/10.0) ) moves.append( mvBranchLengthScale(phylogeny, weight=num_branches) ) ########################################### # Binary morphological substitution model # ########################################### # Specify the number of categories we would like to use to describe our data. For simplicity, # we will use 4. num_cats = 4 # Specify the two parameters to the Beta distribution, and the moves on these parameters. beta_scale ~ dnLognormal( 0.0, sd=2*0.587405 ) moves.append( mvScale(beta_scale, lambda=1, weight=5.0 ) ) # Create the Beta distribution, according to the two parameters and the number of categories. cats := fnDiscretizeBeta(beta_scale, beta_scale, num_cats) # Create the Q matrices. These data are binary, so we initialize the Jukes-Cantor matrix with # two states. The matrix is initialized from the median values of out four quantiles from # the Beta. for (i in 1:cats.size()) { Q[i] := fnF81(simplex(abs(1-cats[i]), cats[i])) } # Tell the model what the probability of a character going into any particular category. # This prior says that a character is equally likely to be put into any category. matrix_probs <- simplex( rep(1,num_cats) ) # Set up Gamma-distributed rate variation. alpha_morpho ~ dnUniform( 0.0, 1E6 ) rates_morpho := fnDiscretizeGamma( alpha_morpho, alpha_morpho, 4 ) # Moves on the parameters to the Gamma distribution. moves.append( mvScale(alpha_morpho, lambda=1, weight=2.0) ) # Combine all of our elements into a CTMC. Because we have not observed any invariant sites, # we specify the coding is 'variable'. Note that we are now feeding the site matrices to the CTMC. phyMorpho ~ dnPhyloCTMC(tree=phylogeny, siteRates=rates_morpho, Q=Q, type="Standard", coding="variable", siteMatrices=matrix_probs) phyMorpho.clamp(morpho) ######## # MCMC # ######## # initialize the model object # mymodel = model(phylogeny) # Create a vector of monitors # # 1. for the full model # monitors.append( mnModel(filename="output/mkv_discretized.log", printgen=10) ) # 2. the tree # monitors.append( mnFile(filename="output/mkv_discretized.trees", printgen=10, phylogeny) ) # 3. and a few select parameters to be printed to the screen # monitors.append( mnScreen(printgen=100) ) # Initialize the MCMC object # mymcmc = mcmc(mymodel, monitors, moves, nruns=2, combine="mixed") # Run the MCMC # mymcmc.run(generations=10000, tuningInterval=200) # check the performance of the MCMC/moves mymcmc.operatorSummary() # Read in the tree trace and construct the consensus tree tree # trace = readTreeTrace("output/mkv_discretized.trees", treetype="non-clock") trace.setBurnin(0.25) # Summarize tree trace and the consensus tree to file mapTree(trace, file="output/mkv_discretized.map.tre") consensusTree(trace, file="output/mkv_discretized.majrule.tre") # Quit RevBayes # q()