################################################################################ # # RevBayes Example: Marginal likelihood estimation under a Jukes-Cantor # substitution model on a single gene. # # authors: Sebastian Hoehna, Tracy A. Heath, Michael Landis and Brian R. Moore # ################################################################################ ####################### # Reading in the Data # ####################### ###### This just defines a single model for all sites ####### ### Read in sequence data for one gene data <- readDiscreteCharacterData("data/primates_and_galeopterus_cytb.nex") # Get some useful variables from the data. We need these later on. n_species <- data.ntaxa() n_branches <- 2 * n_species - 3 taxa <- data.taxa() # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() ###################### # Substitution Model # ###################### #### specify the Jukes-Cantor substitution model applied uniformly to all sites ### Q <- fnJC(4) ############## # Tree model # ############## out_group = clade("Galeopterus_variegatus") # Prior distribution on the tree topology topology ~ dnUniformTopology(taxa, outgroup=out_group) moves.append( mvNNI(topology, weight=5.0) ) moves.append( mvSPR(topology, weight=1.0) ) # Branch length prior for (i in 1:n_branches) { bl[i] ~ dnExponential(10.0) moves.append( mvScale(bl[i]) ) } TL := sum(bl) psi := treeAssembly(topology, bl) ################### # PhyloCTMC Model # ################### # the sequence evolution model seq ~ dnPhyloCTMC(tree=psi, Q=Q, type="DNA") # attach the data seq.clamp(data) ############# # The Model # ############# # We define our model. # We can use any node of our model as a handle, here we chose to use the rate matrix. mymodel = model(Q) monitors.append( mnModel(filename="output/primates_cytb_JC.log",printgen=10, separator = TAB) ) monitors.append( mnFile(filename="output/primates_cytb_JC.trees",printgen=10, separator = TAB, psi) ) ### Compute power posterior distributions pow_p = powerPosterior(mymodel, moves, monitors, "output/pow_p_JC.out", cats=50, sampleFreq=10) pow_p.burnin(generations=10000,tuningInterval=250) pow_p.run(generations=1000) ### Use stepping-stone sampling to calculate marginal likelihoods ss = steppingStoneSampler(file="output/pow_p_JC.out", powerColumnName="power", likelihoodColumnName="likelihood") ss.marginal() ### Use path-sampling to calculate marginal likelihoods ps = pathSampler(file="output/pow_p_JC.out", powerColumnName="power", likelihoodColumnName="likelihood") ps.marginal() # you may want to quit RevBayes now q()