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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
