dotTree fix for bug due to plotrix

I recently reported a bug in dotTree that is caused (but for reasons beyond my comprehension) by updating the package plotrix to it's most recently version.

The bug basically causes all the circles plotted at the tips of the tree to merge into a huge smear. Here's what I mean, using some simulated data:

library(phytools)
dotTree(utree,xu)

## Warning in cos(angles) * radius[circle] + x: longer object length is not a
## multiple of shorter object length

or, for multiple characters:

dotTree(utree,Xu)

## Warning in cos(angles) * radius[circle] + x: longer object length is not a
## multiple of shorter object length

Woah! Here is what a discrete character looks like:

dotTree(tree,y)

## Warning in cos(angles) * radius[circle] + x: longer object length is not a
## multiple of shorter object length

Obviously, this is no good - and the recommendation (such as it was) that I have been giving is to just not update plotrix, or to roll it back if you've already updated the package. What I did today, though, was “fix” dotTree instead (although I'd contend it was never broken) to work around this feature.

The way I did it was by using mapply to call the plotrix function draw.circle multiply across the elements of several vectors, instead of all at once (which previously worked just fine).

Here's a quick demo. Not that reading the source in & setting phylogram<-phytools:::phylogram should not be necessary if you merely update phytools from GitHub:

library(plotrix) ## get plotrix into name space
phylogram<-phytools:::phylogram
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/dotTree.R")
dotTree(utree,xu)

dotTree(utree,Xu)

dotTree(tree,x)

dotTree(tree,y)

On the very last plot, it looks like the legend is also messed up. This must be because add.simmap.legend also uses draw.circle from plotrix internally. I'll have to check this out.

The data for this demo were simulated as follows:

utree<-pbtree(n=26)
utree$tip.label<-LETTERS[26:1]
xu<-fastBM(utree)
Xu<-fastBM(utree,n=6)

tree<-rtree(n=26)
tree$tip.label<-LETTERS[26:1]
x<-fastBM(tree)

y<-sim.history(tree,matrix(c(-1,1,1,-1),2,2))$states

Some updates to phylo.to.map plotting method

I just finished adding some features to the phytools function phylo.to.map for projecting a phylogeny onto a geographic map, and, in particular, to the S3 plotting method.

Mostly, I have added features that I think better automate the spacing of tree, tree labels, links, and the map; but I also contributed some attributes that are aesthetic decisions regarding how the tree and plotted maps look. For instance, I changed the default point type from pch=19 to pch=21, and, by default, I added points to the tips of the tree - although these can be removed using the argument pts.

Let me know if you encounter any new bugs that might be due to this update!

Here is a demo using totally simulated, but somewhat realistic looking data:

library(phytools)
library(mapdata)

## load source code of function from GitHub
## also can be obtained by updating phytools from GitHub using devtools
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/phylo.to.map.R")

## first a world tree & dataset using default settings:
obj<-phylo.to.map(World.tree,World,plot=FALSE)

## objective: 130

## objective: 122
## objective: 122
## objective: 122

## objective: 118
## objective: 118
## objective: 118

## objective: 116
## objective: 116
## objective: 116

## objective: 104
## objective: 104

## objective: 100
## objective: 100
## objective: 100

## objective: 96
## objective: 96
## objective: 96
## objective: 96

## objective: 92

## objective: 88
## objective: 88
## objective: 88
## objective: 88
## objective: 88

obj

## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A geographic map with range:
##      -85.19N, 83.6N
##      -180W, 180W.
## 
## (3) A table containing 26 geographic coordinates.

plot(obj,ftype="i",psize=1.5)

## next, a simulated dataset for Argentina & the package mapdata
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
    regions="Argentina",plot=FALSE)

## objective: 30

## objective: 22

## objective: 16
## objective: 16

## objective: 14
## objective: 14
## objective: 14
## objective: 14
## objective: 14
## objective: 14

## objective: 12

obj

## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A geographic map with range:
##      -55.05N, -21.79N
##      -73.58W, -53.65W.
## 
## (3) A table containing 12 geographic coordinates.

plot(obj,direction="rightwards",ftype="i",colors=c("blue","white"),
    pts=FALSE,psize=1.2)

## finally, another simulated dataset, this time for the island
## of Puerto Rico
obj<-phylo.to.map(PuertoRico.tree,PuertoRico,database="worldHires",
    regions="Puerto Rico",plot=FALSE,rotate=FALSE)
obj

## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 20 tips and 19 internal nodes.
## 
## (2) A geographic map with range:
##      17.92N, 18.52N
##      -67.94W, -65.27W.
## 
## (3) A table containing 20 geographic coordinates.

plot(obj,ftype="off",colors="blue",pch=24,psize=1.5)

The data were simulated more or less as follows:

World.tree<-pbtree(n=26,scale=100)
World.tree$tip.label<-replicate(Ntip(World.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
lat<-fastBM(World.tree,sig2=10,bounds=c(-90,90))
long<-fastBM(World.tree,sig2=80,bounds=c(-180,180))
World<-cbind(lat,long)

Argentina<-locator(n=12) ## with world map plotted
Argentina<-cbind(Argentina$y,Argentina$x)
colnames(Argentina)<-c("lat","long")
Argentina.tree<-pbtree(n=nrow(Argentina),scale=100)
Argentina.tree$tip.label<-replicate(Ntip(Argentina.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
rownames(Argentina)<-Argentina.tree$tip.label

PuertoRico<-locator(n=20) ## with world map plotted
PuertoRico<-cbind(PuertoRico$y,PuertoRico$x)
colnames(PuertoRico)<-c("lat","long")
PuertoRico.tree<-pbtree(n=20,scale=1)
rownames(PuertoRico)<-PuertoRico.tree$tip.label

plot method for geiger's fitDiscrete function

As implicitly promised in a post yesterday, I have now added a plot method for fitted Mk models from the geiger function fitDiscrete. I felt like this was a good idea because I believe the geiger implementation to be more robust (that is, it uses a more rigorous optimization process) than does fitMk in my package. fitMk, by the way, recycles code extensively from ace in the ape package.

What I elected to do was to first pull out all the relevant attributes of the fitted model from the geiger object (an object of class "gfit" resulting from a call to "fitDiscrete"), then use these values to create a new object of class "fitMk" internally, then just call phytools S3 method on that object.

Here is what that wrapper code looks like. Note that an object of class "gfit" in geiger can be from a call to fitDiscrete or fitContinuous, so we have to check for that first.

## S3 plot method for objects resulting from fitDiscrete
plot.gfit<-function(x,...){
    if("mkn"%in%class(x$lik)==FALSE){
        stop("Sorry. No plot method presently available for objects of this type.")
    } else {
        obj<-list()
        QQ<-.Qmatrix.from.gfit(x)
        obj$states<-colnames(QQ)
        m<-length(obj$states)
        obj$index.matrix<-matrix(NA,m,m)
        k<-m*(m-1)
        obj$index.matrix[col(obj$index.matrix)!=row(obj$index.matrix)]<-1:k
        obj$rates<-QQ[sapply(1:k,function(x,y) which(x==y),obj$index.matrix)]
        class(obj)<-"fitMk"
        plot(obj,...)
    }
}

Now let's take our data & tree and fit a model:

library(phytools)
library(geiger)
tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  t67, t68, t76, t78, t79, t74, ...
## 
## Rooted; includes branch lengths.

x

##    t67    t68    t76    t78    t79    t74    t75    t46     t9    t26 
## purple  green    red    red   blue purple  green purple    red purple 
##    t62    t63    t69    t80    t81    t93    t94    t48    t49    t39 
## purple  green  green   blue   blue purple purple  green    red  green 
##    t51    t52    t35    t19     t4    t91    t92    t64    t65    t30 
## purple  green    red purple   blue   blue   blue   blue   blue    red 
##     t8     t2    t72    t86    t87    t66    t88    t89    t29    t10 
##    red   blue    red   blue purple purple   blue  green    red    red 
##    t11    t57    t58    t42    t43    t24    t25    t40    t41    t15 
## purple   blue   blue  green  green purple  green    red   blue   blue 
##    t37    t38    t28    t31    t32    t18    t16    t55    t56    t36 
##   blue  green purple  green   blue   blue   blue   blue   blue purple 
##     t1     t7    t70    t71    t12    t22    t23    t77    t84    t85 
## purple    red purple purple    red    red    red   blue   blue   blue 
##    t73    t95    t96    t50    t47    t27    t53    t54    t97    t98 
##   blue    red    red   blue  green   blue  green   blue  green  green 
##    t59    t90    t99   t100    t17     t3    t33    t34    t44    t45 
## purple    red    red    red  green   blue purple   blue    red purple 
##    t60    t61     t5     t6    t13    t14    t20    t21    t82    t83 
##  green    red   blue  green purple   blue    red   blue  green  green 
## Levels: green blue purple red

fit.phytools<-fitMk(tree,x,model="ARD")
plot(fit.phytools,main="fitMk model=\"ARD\"",show.zeros=FALSE,
    cex.traits=0.9)

fit.geiger<-fitDiscrete(tree,x,model="ARD")
plot(fit.geiger,main="fitDiscrete model=\"ARD\"",show.zeros=FALSE,
    cex.traits=0.9)

Both fitted models are basically identical. Although the configuration of the states in the plotting space are different, it is easy to see that the estimated rates are the same. In other cases, fitDiscrete and fitMk may differ, either because one or the other failed to converge on the MLE of Q, or because the two implementations make slightly different assumptions about the root.

FYI, the generating model in this case was as follows:

print(t(Q))

##        green blue purple red
## green     -2    1      1   0
## blue       1   -3      1   1
## purple     1    1     -3   1
## red        0    1      1  -2

so we are pretty close, particularly in that the two transition classes with their rates set to zero are the only ones with estimated rates of zero in the final fitted model.

Note that the function above uses the geiger internal function (borrowed by phytools) .Qmatrix.from.gfit, so to work that function will need to be loaded into the name space:

.Qmatrix.from.gfit<-phytools:::.Qmatrix.from.gfit

Finally, here's how the data were simulated:

tree<-pbtree(n=100)
Q<-matrix(c(-2,1,1,0,
    1,-3,1,1,
    1,1,-3,1,
    0,1,1,-2),4,4)
rownames(Q)<-colnames(Q)<-c("green","blue","purple","red")
x<-factor(sim.history(tree,Q)$states,
    levels=c("green","blue","purple","red"))