Hi,

I have reconstructed ancestral character states on a phylogeny using MuSSE in the diversitree package and plotted the character state probabilities as pie charts on the nodes. I would, however, like to colour the character states of my extant species, i.e. the tip labels, the same colours as my pie charts, such that all species in  state 1 are e.g. blue, species in state 2 red and  species in state 3 yellow, and have not been successful with my attempts. I am only able to colour them in repeating sets of 3e.g. sp1=blue, sp.2=red, sp.3=yellow, sp.4=blue, sp5=red, sp6=yellow etc. I am also wondering how to colour the branches or edges as the states transition from one to another over time (i.e. as in the "Analyzing diversification with diversitree" manual by Rich FitzJohn on page 23).

Code I've been working with is below:

library(diversitree) #loads library
tree<-read.tree("tree")#loads tree
tree<-chronopl(tree, lambda=1,CV=TRUE) #converts to ultrametric
states<-read.delim("states", header=TRUE)#load states
head(states)

#match states to tree
states<-structure(states$PC, names=states$Species)
names(states)<-tree$tip.label

#MuSSE
diversitree:::argnames.musse(NULL, 3) #number of states
lik<-make.musse(tree, states, 3)
argnames(lik)
#contstrain lambda
lik.base<-constrain(lik, lambda2~lambda1, lambda3~lambda1, mu2~mu1,mu3~mu1,q13~q12,q21~q12,q23~q12,q31~q12,q32~q12)
#find ML point for this model
p<-starting.point.musse(tree, 3)
fit.base<-find.mle(lik.base, p[argnames(lik.base)])

#unconstrained
lik.lambda<-constrain(lik,mu2~mu1,mu3~mu1,q13~q12,q21~q12,q23~q12,q31~q12,q32~q12)
fit.lambda<-find.mle(lik.lambda, p[argnames(lik.lambda)])
anova(fit.base, free.lambda=fit.lambda)

#find ancestral state probabilities
state.probs<-asr.marginal(lik.base, coef(fit.base))#ancestral state probabilities
state.probs
pie.probs<-t(state.probs)
pr<-apply(t(state.probs), 1, which.max)#max probability
tree$node.label<-pr #labels the nodes with the character states
write.tree(tree, file="T")#exports tree with nodes character states to directory #configuration in figtree - node label display label

#tree pie charts
pdf("TREE_PLOT.pdf", height=11, width=8.5)
plot(tree ,cex=.8)
nodelabels(pie=pie.probs,piecol=c("blue","red","yellow"), cex=.5)
dev.off()

	[[alternative HTML version deleted]]

