More user control in co-phylogenetic plotting with phytools

I just added a little more user control of the S3 plot method for the object class "cophylo" for co-phylogenetic plotting.

Firstly, I allow the user to access the internally-used argument part which can be adjusted to permit more space between the two plotted trees. For instance:

library(phytools)
t1

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  R, Q, A, P, F, B, ...
## 
## Rooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 20 tips and 19 internal nodes.
## 
## Tip labels:
##  t4, t6, t13, t16, t15, t14, ...
## 
## Rooted; includes branch lengths.

assoc

##    tree 1 tree 2
## 1       S    t14
## 2       V    t18
## 3       L    t19
## 4       X    t12
## 5       D    t20
## 6       E    t16
## 7       I    t15
## 8       Z    t11
## 9       N     t7
## 10      R     t2
## 11      F     t1
## 12      G     t4
## 13      B     t3
## 14      Y    t13
## 15      U     t8
## 16      Q     t9
## 17      F     t5
## 18      J    t10
## 19      N     t6
## 20      C    t17

obj<-cophylo(t1,t2,assoc=assoc)

## Rotating nodes to optimize matching...
## Done.

plot(obj) ## the default

plot(obj,part=0.3)

Secondly, I added control over the lines linking the tips of the tree to the tip labels. This is because sometimes when rendered as a PDF this lines can seem just too fine! So, for instance:

plot(obj,tip.lwd=2,link.type="curved",link.lwd=3,
    link.col=make.transparent("grey",0.4),
    part=0.35,lwd=2)

## or
plot(obj,tip.lty="solid",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.2),
    part=0.3,lwd=3,pts=F)

Kind of neat. Of course, in this case the two trees have no genuine association - but nonetheless.

These updates can be obtained by installing phytools directly from GitHub using the package devtools.

Running make.simmap in parallel in R for Windows using snow

Shortly after I posted earlier today about parallelizing make.simmap using mcapply, which does not work in R for Windows, a colleague commented that this could indeed be done in Windows, just with slightly greater difficulty.

In fact, this turns out to be precisely true using the CRAN package snow.

Below I demonstrate that this works & that it results (on my machine, which has two cores but 4 hyper-threaded logical processors) in about a 50% speed-up relative to running the same process on a single core.

First, here's our data again:

library(phytools)
cols<-setNames(colorRampPalette(c("blue","red"))(3),
    c("a","b","c"))
dotTree(tree,x,fsize=0.7,ftype="i",colors=cols)

Next we'll fit our model so we don't have to recompute the transition matrix Q for each process.

fit<-fitMk(tree,x,model="ARD")
fittedQ<-matrix(NA,length(fit$states),length(fit$states))
fittedQ[]<-c(0,fit$rates)[fit$index.matrix+1]
diag(fittedQ)<-0
diag(fittedQ)<--rowSums(fittedQ)
colnames(fittedQ)<-rownames(fittedQ)<-fit$states

Now let's run our analysis, first on a single core:

system.time(trees.sc<-make.simmap(tree=tree,x=x,Q=fittedQ,nsim=1000))

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a           b          c
## a -0.6678208  0.66782077  0.0000000
## b  0.3374750 -0.71533391  0.3778589
## c  0.9894179  0.07069738 -1.0601152
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

##    user  system elapsed 
##   47.60    0.06   47.90

trees.sc

## 1000 phylogenetic trees with mapped discrete characters

now on multiple cores using snow:

library(snow)
cl<-makeSOCKcluster(rep("localhost",4))
system.time(trees.mc<-clusterApply(cl,x=replicate(4,x,simplify=FALSE),
    fun=make.simmap,tree=tree,Q=fittedQ,nsim=250))

##    user  system elapsed 
##    0.11    0.07   24.94

trees.mc<-do.call(c,trees.mc)
if(!("multiSimmap"%in%class(trees.mc))) 
    class(trees.mc)<-c("multiSimmap",class(trees.mc))
stopCluster(cl)
trees.mc

## 1000 phylogenetic trees with mapped discrete characters

I did one peculiar thing which is sending a list of length equivalent to the number of cores consisting of replicates of our character vector x. This is because both clusterApply and make.simmap take the argument x. I could have instead written a simple wrapper function with different argument names if I wanted.

It's easy to show that our results are (more or less - remember this is stochastic mapping) identical in the two cases:

plot(summary(trees.sc)$ace,summary(trees.mc)$ace,cex=1.5,pch=21,
    bg="grey",xlab="posterior probabilities from single core",
    ylab="posterior probabilities from multi-core")

That's pretty neat right?

Thanks to Will Gearty for the tip.

Running make.simmap in parallel

Today a phytools user asked on my GitHub page if there was any way to permit make.simmap run in parallel.

Note that the issue motivating this report likely stems from a replacement of the method to do matrix exponentiation in the latest CRAN version of phytools. This has already been fixed - although, of course, if we are using make.simmap with large trees, complex models, or generating many samples, it will still be slow!

make.simmap does not have a parallelization option; however, for the basic method we are just estimating the transition matrix, Q, via Maximum Likelihood, and then sampling histories for the discrete character conditioned on this fitted model. Consequently, we ought to be able to first compute Q once and then run our sampler on multiple cores. In theory, this can be done using parallel::mclapply.

Unfortunately, parallelization is not supported in R for Windows. Nonetheless, it is possible to run mclapply which pseudo-parallelizes our analysis. This is so that functions & packages using parallelization don't break on Windows - but it will not run any faster!

Here is what that might look like. I have assumed we have 4 cores to play with.

library(phytools)
## our data:
cols<-setNames(colorRampPalette(c("blue","red"))(3),
    c("a","b","c"))
dotTree(tree,x,fsize=0.7,ftype="i",colors=cols)

## first fit our model
fit<-fitMk(tree,x,model="ARD")
fit

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.667821  0.667821  0.000000
## b  0.337475 -0.715334  0.377859
## c  0.989418  0.070697 -1.060115
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -51.296883 
## 
## Optimization method used was "nlminb"

## extracted the fitted transition matrix:
fittedQ<-matrix(NA,length(fit$states),length(fit$states))
fittedQ[]<-c(0,fit$rates)[fit$index.matrix+1]
diag(fittedQ)<-0
diag(fittedQ)<--rowSums(fittedQ)
colnames(fittedQ)<-rownames(fittedQ)<-fit$states
## ready to run our analysis
library(parallel)
trees<-mclapply(1:4,function(n,tree,x,fixedQ) make.simmap(tree,x,Q=fixedQ,nsim=100),
    tree=tree,x=x,fixedQ=fittedQ,mc.cores=if(.Platform$OS.type=="windows") 1L else 4L)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a           b          c
## a -0.6678208  0.66782077  0.0000000
## b  0.3374750 -0.71533391  0.3778589
## c  0.9894179  0.07069738 -1.0601152
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a           b          c
## a -0.6678208  0.66782077  0.0000000
## b  0.3374750 -0.71533391  0.3778589
## c  0.9894179  0.07069738 -1.0601152
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a           b          c
## a -0.6678208  0.66782077  0.0000000
## b  0.3374750 -0.71533391  0.3778589
## c  0.9894179  0.07069738 -1.0601152
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a           b          c
## a -0.6678208  0.66782077  0.0000000
## b  0.3374750 -0.71533391  0.3778589
## c  0.9894179  0.07069738 -1.0601152
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

trees

## [[1]]
## 100 phylogenetic trees with mapped discrete characters
## 
## [[2]]
## 100 phylogenetic trees with mapped discrete characters
## 
## [[3]]
## 100 phylogenetic trees with mapped discrete characters
## 
## [[4]]
## 100 phylogenetic trees with mapped discrete characters

At the end, though, instead of one "multiSimmap" object, we have a list of 4! Fortunately, we can combine these in a fairly easy way as follows:

trees<-do.call(c,trees)
if(!("multiSimmap"%in%class(trees))) class(trees)<-c("multiSimmap",class(trees))
trees

## 400 phylogenetic trees with mapped discrete characters

Neat.

Now, let's compute a summary & visualize the results:

obj<-summary(trees)
obj

## 400 trees with a mapped discrete character with states:
##  a, b, c 
## 
## trees have 43.7325 changes between states on average
## 
## changes are of the following types:
##        a,b a,c   b,a   b,c    c,a  c,b
## x->y 18.29   0 7.945 8.745 8.1425 0.61
## 
## mean total time spent in each state is:
##               a          b         c    total
## raw  27.2549076 23.1011251 8.2695555 58.62559
## prop  0.4648978  0.3940451 0.1410571  1.00000

plot(obj,colors=cols,ftype="i",fsize=0.7)
add.simmap.legend(colors=cols,prompt=FALSE,vertical=FALSE,
    x=0,y=2)

The tree & data for this example were simulated as follows:

tree<-pbtree(n=60)
Q<-matrix(c(-1,1,0,
    1,-2,1,
    0,1,-1),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x<-as.factor(sim.history(tree,Q)$states)

I would love to hear a report from someone working on a Unix-alike system (e.g., Mac OS) who can report if this indeed works on their machine!