Dear all,
I need some help with the tree plotting functions of surface and bayou.
I have tree outputs for my tested character in surface and also in bayou.
If I understood it right the surface output of the function surfaceTreeplot () shows me where changes happened. And the tree looks rather nice even if we have a lot of changes here.
The bayou tree should show me two things using the function plotSimmap.mcmc (), to plot the object of the former used function load.bayou()
1) the color of the edgelabels should show me whether we have an increase or a decrease in our trait and the size of the circle should show me the probability of it.
Firstly, my circles are all red :( which shouldn't be the case
and secondly the whole plot looks like a mess.
But my real problem is that I want to macth both trees and I just don't know how.
I read these really nice papers (e.g. Cuff et al. 2015; Big cat, small cat:reconstructing body size evolution in living and extinct Felidae), where they matched the trees and showed the surface and the bayou output in the same tree.
Maybe you guys can help me with that.
###create first chain
par(mfrow=c(2,3))
fit1 <- bayou.mcmc(tree, dat, SE=SE, model="OU", prior, ngen=1000000, new.dir=getwd(), plot.freq=10000, ticker.freq=1000)
fit1
surfaceTreePlot2 <-
function(tree,hansenfit,cols=NULL,convcol=TRUE,labelshifts=FALSE,...){
fit<-hansenfit$fit[[1]]
otree<-as(fit,"data.frame")
otree<-data.frame(otree,shifts=rep(NA,length(otree$nodes)))
otree$shifts[match(names(hansenfit$savedshifts),otree$nodes)]<-1:length(hansenfit$savedshifts)
ntip<-(dim(otree)[1]+1)/2;nnode<-ntip-1
otree2<-otree[match(c(tree$tip.label,tree$node.label),otree$labels),]
otree2<-otree2[tree$edge[,2],]
if(length(cols)==1)cols<-rep(cols,length(unique(hansenfit$savedshifts)))
if(is.null(cols)){
xx<-summary(factor(hansenfit$savedshifts))
if(convcol){
cols<-character(length(xx))
cols[xx>1]<-rainbow(sum(xx>1))
if(any(xx==1))
cols[xx==1]<-c("black",grey(seq(0.7,0.3,length.out=sum(xx==1)-1)))
}else{
cols<-c("black",rainbow(length(xx)-1))
}
}
edgecols<-cols[as.numeric(factor(otree2[,5]))]
plot(tree,cex=0.4,type="fan",tips=NULL,edge.color=edgecols,...)
if(labelshifts){
nodelabels(node=tree$edge[,2][which(!is.na(otree2$shifts))],bg="white", text=otree2$shifts[!is.na(otree2$shifts)],cex=0.3,frame="circle")
}
}
Thank you in advance for any help
####bayou
### create prior
par(mar=c(1,1,1,1))
prior <- make.prior(tree, dists=list(dalpha="dhalfcauchy", dsig2="dhalfcauchy",dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(scale=1), dsig2=list(scale=1), dk=list(lambda=1, kmax=12), dsb=list(bmax=1,prob=1), dtheta=list(mean=mean(dat), sd=2)))
prior <- make.prior(tree, dists=list(dalpha="dhalfcauchy", dsig2="dhalfcauchy",dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(scale=1), dsig2=list(scale=1), dk=list(lambda=1, kmax=12), dsb=list(bmax=1,prob=1), dtheta=list(mean=mean(dat), sd=2)))
###create first chain
par(mfrow=c(2,3))
fit1 <- bayou.mcmc(tree, dat, SE=SE, model="OU", prior, ngen=1000000, new.dir=getwd(), plot.freq=10000, ticker.freq=1000)
fit1
chain <- load.bayou(fit1, save.Rdata=FALSE, cleanup=TRUE)
chain <- set.burnin(chain, 0.3)
chain <- set.burnin(chain, 0.3)
out <- summary(chain)
plot(chain)
### plot tree & phenogram
par(mfrow=c(1,1))
plotSimmap.mcmc(chain, burnin=0.3, lwd=2, edge.type="theta", pal=colorRampPalette(c("black", "gray90","red","blue")), show.tip.label=TRUE,cex=0.4,type="fan")
plotSimmap.mcmc(chain, burnin=0.3, lwd=2, edge.type="theta", pal=colorRampPalette(c("black", "gray90","red","blue")), show.tip.label=TRUE,cex=0.4,type="fan")
phenogram.density(tree, dat, chain=chain, burnin=0.3, pp.cutoff=0.3)
####surface
surfaceTreePlot2 <-
function(tree,hansenfit,cols=NULL,convcol=TRUE,labelshifts=FALSE,...){
fit<-hansenfit$fit[[1]]
otree<-as(fit,"data.frame")
otree<-data.frame(otree,shifts=rep(NA,length(otree$nodes)))
otree$shifts[match(names(hansenfit$savedshifts),otree$nodes)]<-1:length(hansenfit$savedshifts)
ntip<-(dim(otree)[1]+1)/2;nnode<-ntip-1
otree2<-otree[match(c(tree$tip.label,tree$node.label),otree$labels),]
otree2<-otree2[tree$edge[,2],]
if(length(cols)==1)cols<-rep(cols,length(unique(hansenfit$savedshifts)))
if(is.null(cols)){
xx<-summary(factor(hansenfit$savedshifts))
if(convcol){
cols<-character(length(xx))
cols[xx>1]<-rainbow(sum(xx>1))
if(any(xx==1))
cols[xx==1]<-c("black",grey(seq(0.7,0.3,length.out=sum(xx==1)-1)))
}else{
cols<-c("black",rainbow(length(xx)-1))
}
}
edgecols<-cols[as.numeric(factor(otree2[,5]))]
plot(tree,cex=0.4,type="fan",tips=NULL,edge.color=edgecols,...)
if(labelshifts){
nodelabels(node=tree$edge[,2][which(!is.na(otree2$shifts))],bg="white", text=otree2$shifts[!is.na(otree2$shifts)],cex=0.3,frame="circle")
}
}
surfaceTreePlot2(tree,newout$bwd[[newkk]],labelshifts = F, convcol = TRUE)
_______________________________________________ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
