# Start to create helper variables for file handling, area_fn = "data/atlas_folded.txt" data_fn = "data/king_p.nex" tree_fn = "data/king.tre" out_fn = "output/p_dec_dist_noclado_folded" # read in the tree psi <- readTrees(tree_fn)[1] # load range observations (i.e. plumages) data = readDiscreteCharacterData(data_fn) # read in geographical info, epochs atlas = readAtlas(area_fn) n_areas = atlas.nAreas() # create counter for move vectors mvi = 0 # setup number of generations in MCMC analysis ngen = 100000 # save copy of the tree annoated with node indexes for comparing with output later write(psi, filename=tree_fn+".index.tre") # biogeographical clock for scaling rate of range evolution clock_bg ~ dnExponential(10) moves[++mvi] = mvScale(clock_bg, lambda=0.5, weight=5.0) # simplex of gain/loss rates (flat Dirichlet prior distribution) glr ~ dnDirichlet([1,1]) moves[++mvi] = mvSimplexElementScale(glr, alpha=30.0, weight=5.0) # deterministic nodes to give rates nicknames r_gain := glr[1] r_loss := glr[2] # average rate of area gain and loss per area q_area := fnFreeBinary(glr) # "dp" to represent beta parameter (importance of distance to dispersal) # beta ~ 0 (distance not important), beta !=0 (distance important) dp ~ dnExponential(10.0) moves[++mvi] = mvScale(x=dp, lambda=0.5, tune=true, weight=5.0) # uniform prior gives some weird results # dp ~ dnUniform(-1, 1) # moves[++mvi] = mvSlide(dp, weight=5.0) grm := fnBiogeoGRM(atlas=atlas, distancePower=dp, useDistance=true) # rate matrix Q q_range := fnBiogeoDE(gainLossRates=q_area, branchRates=clock_bg, geoRateMod=grm, numAreas=n_areas, forbidExtinction=false) # prior for cladogenetic event type probs clado_prob ~ dnDirichlet([1, 1e-3, 1e-3]) widespread_sympatry := clado_prob[1] subset_sympatry := clado_prob[2] allopatry := clado_prob[3] moves[++mvi] = mvSimplexElementScale(clado_prob, alpha=20.0, weight=5.0) # data evolving by a phylogenetic CTMC process (=continuous time markov chain, I think) m ~ dnPhyloDACTMC(tree=psi, Q=q_range, cladoProbs=clado_prob, type="Biogeo", forbidExtinction=false, useCladogenesis=false) # tell CTMC to observe data (primes model) m.clamp(data) m.lnProbability() # number of nodes (I expected ntips - 1, got 2X ntips.. ??) n_nodes <- psi.nnodes() # moderate character history updates for paths and nodes moves[++mvi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=0.2, type="Biogeo", graph="node", proposal="rejection", weight=n_nodes) moves[++mvi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=0.2, type="Biogeo", graph="branch", proposal="rejection", weight=n_nodes*2) # create model object my_model = model(m) # monitor simple parameters monitors[1] = mnScreen(clock_bg, r_gain, r_loss, dp, subset_sympatry, allopatry, widespread_sympatry, printgen=100) monitors[2] = mnFile(clock_bg, r_gain, r_loss, dp, subset_sympatry, allopatry, widespread_sympatry, filename=out_fn+".params.txt", printgen=10) # these monitors record sampled char history states for each node in tree # events - richer info along branch, anagenetic + cladogenetic monitors[3] = mnCharHistoryNewick(filename=out_fn+".events.txt", ctmc=m, tree=psi, printgen=100, style="events") # counts - report # gain/losses per branch in tab-delim file (good for Tracer) monitors[4] = mnCharHistoryNewick(filename=out_fn+".counts.txt", ctmc=m, tree=psi, printgen=100, style="counts") # create MCMC object my_mcmc = mcmc(my_model, monitors, moves, nruns=1) # and run it my_mcmc.run(generations=ngen) q()