# We are going to have a little fun comparing two new approaches to identifying rate shifts on a topology. These are the approaches of Revell et al. 2012 and Eastman et al. 2011. This is just a quick introduction to these new methods to show you what is possible, as they both implement MCMC we are not going to be able to get post burn-in results during the time available.

# Both these methods implement MCMC so you will need to apply what you have learnt throughout the workshop to implement them correctly - in particular refer to Brian Moore's tutorial on MCMC diagnosis. So if you want to implement these methods in your thesis or papers please spend *a lot* of time familiarizing yourself with the method you choose; its benefits and pitfalls as well as how to implement it correctly i.e. checking that the MCMC has converged, is appropriately mixed and not autocorrelated etc..


## Eastman et al. auteur package
# auteur implements a reversible-jump MCMC approach to find one or more shifts in the Brownian motion rate throughout the topology. For precise details about the RJMCMC approach they implement please see the paper: Eastman et al. 2011 Evolution 65(12), 3578-3589. This method is ideal for data exploration, to help you identify what may be interesting patterns pursue, it is unlikely to be sufficient on its own.   

trait<-dat[,1]
names(trait)<-row.names(dat)

# running two short reversible-jump Markov chains

rjmcmc.bm(gruntrs[[1]], trait, summary=TRUE, reml=TRUE, ngen=100000, sample.freq=500, constrainK=1, fileBase="trait1tr1_test1" ) # We are constraining it to only evaluate models with an additional independent rate using constrainK=1 so that it is comparable to the Revell et al. approach. 
rjmcmc.bm(gruntrs[[1]], trait, summary=TRUE,reml=TRUE, ngen=100000,  sample.freq=500, fileBase="trait1tr1_test2")

# This is where you should do your MCMC diagnosis - we are going to pretend you have done it for now!

pool.rjmcmcsamples(base.dirs=c("BM.trait1tr1_test1.parameters", "BM.trait1tr1_test2.parameters"), lab="trait1tr1")

shifts.plot(phy=gruntrs[[1]], base.dir="trait1tr1.combined.rjmcmc", burnin=0.25, legend=TRUE, edge.width=2)

## Revell et al. 2012 evol.rate.mcmc and companion functions in the phytools package
# This method uses a Bayesian MCMC approach to estimate a single rate shift upon a tree. 

gruntr<-read.tree(text=write.tree(gruntrs[[1]]))

gruntres<-evol.rate.mcmc(gruntr, trait, ngen=100000, control=list(sample=200))

#Defaults symmetric exponential proposal distribution and Gaussian proposal distributions centered on 0 for the two rates and the ancestral value (at root). Here we have set the 100,000 generations and sampling every 200. 

# This is where you should do all your MCMC diagnosis - for now we are going to pretend you have done it! For now we will remove the first 25% as burn-in which is done by only using the last 75% of the output gruntres$mcmc[125:501, ]

# Now calculate the split with the smallest summed distances to all other splits.

gruntsminsplit<-minSplit(gruntr, gruntres$mcmc[125:501, c("node", "bp")], method="sum")

# Extract the posterior sample of the average split

gruntpost<-posterior.evolrate(gruntr, ave.shift=gruntsminsplit,gruntres$mcmc[125:501,], gruntres$tips[125:501])

pies<-hist(gruntres$mcmc[125:501,][,"node"], breaks=0.5:length(gruntr$tip)+gruntr$Nnode+0.5)

plot(gruntr)

# We can plot the posterior probabilities on the tree using pie charts on the edges (rate shift has to appear somewhere so the probabilities across the tree have to sum to 1)
edgelabels(edge=match(which(pies$density>0.01), gruntr$edge[,2]), pie=pies$density[which(pies$density>0.01)], piecol=c("red", "black"), cex=0.5)

# Really useful to look at the code used to run the analyses in paper presented in the appendix SI1. Also see Liam Revell's blog (http://phytools.blogspot.com)

# The most obvious difference between the methods in aueter and phytools is that auteur allows multiple rate shifts and even when we constrain it to two rate classes to be comparable to the Revell et al. method it can assign branches across the tree to the new rate (doesn't have to be a monophyletic clade).  

# These methods are new and relatively untested, use at your own risk and ensure you have run the appropriate MCMC diagnostics.