# read stratigraphic ranges taxa = readTaxonData(file = "data/dinosaur_ranges.tsv") # read fossil counts k <- readDataDelimitedFile(file = "data/dinosaur_fossil_counts.tsv", header = true, rownames = true) # interval boundaries timeline <- v(100, 145, 201) - 66 # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() alpha <- 10 # specify FBDR model parameters for(i in 1:(timeline.size()+1)) { mu[i] ~ dnExp(alpha) lambda[i] ~ dnExp(alpha) psi[i] ~ dnExp(alpha) div[i] := lambda[i] - mu[i] turnover[i] := mu[i]/lambda[i] moves.append( mvScale(mu[i], lambda = 0.01) ) moves.append( mvScale(mu[i], lambda = 0.1) ) moves.append( mvScale(mu[i], lambda = 1) ) moves.append( mvScale(lambda[i], lambda = 0.01) ) moves.append( mvScale(lambda[i], lambda = 0.1) ) moves.append( mvScale(lambda[i], lambda = 1) ) moves.append( mvScale(psi[i], lambda = 0.01) ) moves.append( mvScale(psi[i], lambda = 0.1) ) moves.append( mvScale(psi[i], lambda = 1) ) } rho <- 0 # model 3 bd ~ dnFBDRMatrix(taxa=taxa, lambda=lambda, mu=mu, psi=psi, rho=rho, timeline=timeline, k=k, binary=true) moves.append( mvMatrixElementScale(bd, lambda = 0.01, weight=taxa.size()) ) moves.append( mvMatrixElementScale(bd, lambda = 0.1, weight=taxa.size()) ) moves.append( mvMatrixElementScale(bd, lambda = 1, weight=taxa.size()) ) moves.append( mvMatrixElementSlide(bd, delta = 0.01, weight=taxa.size()) ) moves.append( mvMatrixElementSlide(bd, delta = 0.1, weight=taxa.size()) ) moves.append( mvMatrixElementSlide(bd, delta = 1, weight=taxa.size()) ) mymodel = model(alpha) # add monitors monitors.append( mnScreen(lambda, mu, psi, printgen=100) ) monitors.append( mnModel(filename="output/model3.log",printgen=10) ) # monitors to print RevGagets input monitors.append( mnFile(filename="output/model3_speciation_rates.log",lambda,printgen=10) ) monitors.append( mnFile(filename="output/model3_speciation_times.log",timeline,printgen=10) ) monitors.append( mnFile(filename="output/model3_extinction_rates.log",mu,printgen=10) ) monitors.append( mnFile(filename="output/model3_extinction_times.log",timeline,printgen=10) ) monitors.append( mnFile(filename="output/model3_sampling_rates.log",psi,printgen=10) ) monitors.append( mnFile(filename="output/model3_sampling_times.log",timeline,printgen=10) ) # run the analysis mymcmc = mcmc(mymodel, moves, monitors, moveschedule="random") mymcmc.run(30000) q()