> library(BAMMtools) #read in tree >tree<-read.nexus("buzztreeBam.nex") #read in event data >edata <- getEventData(tree, eventdata = "Bamm6b_event_data.txt", burnin=0.1) #assess MCMC convergence #plot the log-likelihood trace of your MCMC output file > mcmcout <- read.csv("Bamm6b_mcmc_out.txt", header=T) > plot(mcmcout$logLik ~ mcmcout$generation) #discard burnin (first 10%) > burnstart <- floor(0.1 * nrow(mcmcout)) > postburn <- mcmcout[burnstart:nrow(mcmcout), ] #check effective sample sizes of the log-likelihood and the number of shift events present in each sample #do this using the coda library for R #these should be at least 200 > library(coda) > effectiveSize(postburn$N_shifts) > effectiveSize(postburn$logLik) #analysis of rate shifts #1-how many rate shifts #compute the posterior probabilities of models sampled post_probs <- table(postburn$N_shifts) / nrow(postburn) #see which models are part of the set that were sampled > names(post_probs) #any model that is not included in the names(post-probs) was so lousy it was never sampled, #if 0 is not in this set, this means you have overwhelming evidence for diversification rate heterogeneity #and that probability of no shifts if effectively 0 #alternatively, can summarize the posterior distribution of the number of shifts > shift_probs <- summary(edata) #Bayes factors for model comparison to get a pairwise matrix of Bayes factors #usually the best model is the model with the highest Bayes factor relative to the null model Mo > postfile <- "Bamm6b_mcmc_out.txt" > bfmat <- computeBayesFactors(postfile, expectedNumberOfShifts=1, burnin=0.1) > bfmat #look at first column and select number of shifts with highest bayes factor when compared to the model with the lowest number of supported shifts (often zero, but check to see if 0 is part of 95% credible set of rate shift configurations) #trying to see where to most likely positions on tree of rate shifts #bayesian credible sets of shift configurations #95% credible set is the set of distinct shift configurations that account for 95% of the probability of the data #estimate the credible set of rate shifts css <- credibleShiftSet(edata, expectedNumberOfShifts=1, threshold=5, set.limit = 0.95) #number of distinct shift configurations in the data css$number.distinct #view more info about the cridible set summary(css) #generate phylorate plots for each of the N shift conficuations with the highest post. probabilities plot.credibleshiftset(css) #text above plot=post.prob of each shift configuration #red circle=rate acceleration, blue= rate slowdown #find the single best shift configuration #if showing single set of rate shifts for pub, show the maximum a posteriori (MAP) probability # MAP is the one estimate of the overall best rate set of rate shifts given your data best <- getBestShiftConfiguration(edata, expectedNumberOfShifts=1) plot.bammdata(best, lwd = 2) addBAMMshifts(best, cex=2.5) #to plot to pdf instead best <- getBestShiftConfiguration(edata, expectedNumberOfShifts=1) pdf("Bamm6b.pdf", width=8, height=30) plot.bammdata(best, lwd = 2) addBAMMshifts(best, cex=2.5) dev.off() #correlation analysis dm<-read.table("bamm_buzzcodes.txt",header=T,sep=" ",stringsAsFactors = F) trait<-dm$buzz names(trait)<-dm$species trait<-trait[! is.na(trait)] #default rate='speciation', could do rate='net diversification'; method="m" = mann-whitney, options 'spearman'.. bee_dm<-traitDependentBAMM(ephy = edata,rate='speciation', traits = trait,reps = 10000,return.full = T,logrates = T, method= "m", two.tailed = T) #or: #bee_dm<-traitDependentBAMM(ephy = edata, rate='net diversification', traits = trait,reps = 10000,return.full = T,logrates = T, method= "m", two.tailed = T) bee_dm #check the difference in absolute spearman's correlation coefficients #bee_dm_diff<-abs(bee_dm$obs.corr)-abs(bee_dm$null)