library(expoTree) library(TreePar) library(laser) #### Condensed version of TreePar Analysis ############################################ tree<-read.tree("YOUR.TREE.tre") l<-sort(branching.times(tree),decreasing=TRUE) l.noshift<-bd.shifts.optim(l, sampling=c(1), survival=1)[[2]] #above estimates loglikelihood, turnover (extinction), and diversification (speciation-extinction) #use miniall=(previous) to give it the sampling from your previous run and save time! l.oneshift <-bd.shifts.optim(l, sampling=c(1,1), miniall=l.noshift, survival=1, grid=1, start=0, end=60)[[2]] l.twoshift <-bd.shifts.optim(l, sampling=c(1,1,1), miniall=l.oneshift, survival=1, grid=1, start=0, end=60)[[2]] l.threeshift <-bd.shifts.optim(l, sampling=c(1,1,1,1), miniall=l.twoshift, survival=1, grid=1, start=0, end=60)[[2]] l.fourshift <-bd.shifts.optim(l, sampling=c(1,1,1,1,1), miniall=l.threeshift,survival=1, grid=1, start=0, end=60)[[2]] l.fiveshift <-bd.shifts.optim(l, sampling=c(1,1,1,1,1,1), miniall=l.fourshift,survival=1, grid=1, start=0, end=60)[[2]] l.sixshift<-bd.shifts.optim(l, sampling=c(1,1,1,1,1,1,1), miniall=l.fiveshift, survival=1, grid=1, start=0, end=60)[[2]] i<-1 comparison.0.1<-pchisq(2*(p.fourshift[[i]][1]-p.fourshift[[i+1]][1]),3) #compare noshift vs. oneshift via chi square (we reject once value < 0.95) comparison.1.2<-pchisq(2*(p.fourshift[[i+1]][1]-p.fourshift[[i+2]][1]),3) #oneshift vs. twoshift comparison.2.3<-pchisq(2*(p.fourshift[[i+2]][1]-p.fourshift[[i+3]][1]),3) #twoshift vs. threeshift comparison.3.4<-pchisq(2*(p.fourshift[[i+3]][1]-p.fourshift[[i+4]][1]),3) #threeshift vs. fourshift comparison.4.5<-pchisq(2*(p.sixshift[[i+4]][1]-p.sixshift[[i+5]][1]),3) #etc comparison.5.6<-pchisq(2*(p.sixshift[[i+5]][1]-p.sixshift[[i+6]][1]),3) #etc plot0<-bd.shifts.plot(list(l.noshift),0,60,0,1) # (num.shifts ,x.range, , y.range) plot1<-bd.shifts.plot(list(l.oneshift),1,60,0,1) # (num.shifts ,x.range, , y.range) plot2<-bd.shifts.plot(list(l.twoshift),2,60,0,1) # (num.shifts ,x.range, , y.range) plot3<-bd.shifts.plot(list(l.threeshift),3,60,0,1) # (num.shifts ,x.range, , y.range) plot4<-bd.shifts.plot(list(l.fourshift),4,60,0,1) # (num.shifts ,x.range, , y.range) plot5<-bd.shifts.plot(list(l.fiveshift),5,60,0,1) # (num.shifts ,x.range, , y.range) ###################################################################### # Create a loop to compare two (or more!) models across a set of trees ###################################################################### trees <- read.tree("YOUR TREES.trees") #multiPhylo object me.vs.shift <- NULL for (i in 1:100) { x <- sort(branching.times(trees[[i]]),decreasing=TRUE) #pulls and sorts the branching times in the tree me.res <- bd.shifts.optim(x, sampling=c(1,0.1), ME=F, all=T, start=27, grid=0.5, end=34)[[2]] #this is a model with mass extinction between 27-34 million years alt.me <- bd.shifts.optim(x, sampling=c(1,1), ME=F, start=27, grid=0.5, end=34)[[2]] #this is a model with a speciation rate shift between 27-34 million years #construct a vector of the likelihoods, named by the models candidateModels <- c("MassEx"=-(me.res[[2]][1]), "Shift30"=-(alt.me[[2]][1])) #make a grid with all possible combinations of the models LikelihoodGrid <- expand.grid(M0=names(candidateModels), M1=names(candidateModels)) #add a column that is the 2lnBF for each pair of models LikelihoodGrid$BF <- 2 * (candidateModels[LikelihoodGrid$M0] - candidateModels[LikelihoodGrid$M1]) #sort the comparisons by their 2lnBF in descending order LikelihoodGrid <- LikelihoodGrid[order(LikelihoodGrid$BF, decreasing=TRUE),] me.vs.shift <- rbind(me.vs.shift, LikelihoodGrid[1,]) } write.table(me.vs.shift, file="YOUR DIRECTORY/TreePar.Model.Comparison.table.xls", sep="\t", quote=F, row.names=F) # write it all to a file table(me.vs.shift[,1]) # you can summarize the frequency of the models with 'table'