# # This file contains the bacula example of the graphical models paper (Example 2). # # Example file for a Bayesian analysis using MCMC. # The data set contains 274 species each with a single binary observation for the baculum (absent/present). # # author: Sebastian Hoehna # # read the data # the readCharacter function returns a vector of matrices. We just take the first one D <- readCharacterData("data/baculumData01.nex")[1] # get some useful variables from the data #taxaCount <- D.ntaxa() not used here nSites <- D.nchar()[1] #names <- D.names() not used here # create the parameter for the root frequencies # instead of a beta distribution we use the Dirichlet distribution because the data type needs to be a simplex # set a flat prior rf_prior <- [1,1] # create the random variable for the root frequencies rf ~ dirichlet(rf_prior) # create a move/proposal that changes the root-frequencies during the MCMC moves[1] <- mSimplexElementScale(rf, alpha=10.0, tune=true, weight=2.0) # now let us create the random variables for the continuous time Markov process # that changes the absent/present state along the tree # we use the F81 rate matrix again with a Dirichlet distribution as the prior on the base-frequencies bf_prior <- [1,1] # create the random variable for the base-frequencies bf ~ dirichlet(bf_prior) # the rate matrix Q := F81(bf) # create a move/proposal that changes the base-frequencies during the MCMC moves[2] <- mSimplexElementScale(bf, alpha=10.0, tune=true, weight=2.0) # We use a fixed tree. # tree (dos Reis et al.) from file tau <- readTrees("data/mammalia_dosReis.tree")[1] # Just use the default clock rate. # We could also ommit this parameter. clockRate <- 1.0 # the sequence evolution model seq ~ substModel(tree=tau, Q=Q, branchRates=clockRate, rootFrequencies=rf, nSites=nSites, type="Standard") # attach the data seq.clamp(D) # create the model from the DAG mymodel <- model(Q) monitors[1] <- modelmonitor(filename= "GraphicalModels_Example_2.log",printgen=10, separator = " ") monitors[2] <- screenmonitor(printgen=10, separator = " ", bf, rf) mymcmc <- mcmc(mymodel, monitors, moves) # If you choose more or different proposals, or different weights for the proposals, then the number of proposals changes per iteration. Currently there is only one proposal with weight 1.0. mymcmc.burnin(generations=2000,tuningInterval=100) mymcmc.run(generations=200000) result <- readTrace("GraphicalModels_Example_2.log") result[5]