Generic "combine" (c) method for "simmap" and "multiSimmap" object classes

This afternoon I spent an inordinate amount of time trying to create generic c (“combine”) method for the phytools "simmap" and "multiSimmap" object types.

I’m not sure how I hadn’t done this before, but perhaps it’s because the c.phylo and c.multiPhylo methods from ape work perfectly well on these object types – just producing "multiPhylo" rather than "multiSimmap" objects in return!

(Normally, this will probably be OK, but there are some circumstances in which it will cause problems.)

At first, I imagined writing the simple method as follows.

c.simmap<-function(...){
  object<-c.phylo(...)
  class(object)<-union("simmap",class(object))
  object
}

Unfortunately, c.phylo is not exported by the namespace of ape and replacing it with ape:::c.phylo is not allowed by CRAN policies!

Instead, I thought I would copy & (lightly) adapt the code of ape:::c.phylo as follows.

## adapted from c.phylo in ape (Paradis & Schliep 2019)
c.simmap<-function(...,recursive=TRUE){
  obj<-list(...)
  classes<-lapply(obj,class)
  isphylo<-sapply(classes,function(x) "simmap" %in% x)
  if(all(isphylo)){
    class(obj)<-c("multiSimmap","multiPhylo")
    return(obj)
  }
  if(!recursive) return(obj)
  ismulti<-sapply(classes, function(x) "multiSimmap" %in% x)
  if(all(isphylo|ismulti)){
    result<-list()
    j<-1
    for(i in 1:length(isphylo)){
      if(isphylo[i]){ 
        result[[j]]<-obj[[i]]
        j<-j+1
      } else {
        n<-length(obj[[i]])
        result[0:(n-1)+j]<-.uncompressTipLabel(obj[[i]])
        j<-j+n
      }
    }
    class(result)<-c("multiSimmap","multiPhylo")
    obj<-result
  } else {
    msg<-paste("some objects not of class \"simmap\" or",
      "\"multiSimmap\": argument recursive=TRUE ignored")
    warning(msg)
  }
  return(obj)
}

We can do a similar thing for the "multiSimmap" object.

## adapted from c.multiPhylo in ape (Paradis & Schliep 2019)
c.multiSimmap<-function(...,recursive=TRUE){
  obj<-list(...)
  if(!recursive) return(obj)
  classes<-lapply(obj,class)
  isphylo<-sapply(classes,function(x) "simmap" %in% x)
  ismulti<-sapply(classes, function(x) "multiSimmap" %in% x)
  if(all(isphylo|ismulti)){
    result<-list()
    j<-1
    for(i in 1:length(isphylo)){
      if(isphylo[i]){ 
        result[[j]]<-obj[[i]]
        j<-j+1
      } else {
        n<-length(obj[[i]])
        result[0:(n-1)+j]<-.uncompressTipLabel(obj[[i]])
        j<-j+n
      }
    }
    class(result)<-c("multiSimmap","multiPhylo")
    obj<-result
  } else {
    msg<-paste("some objects not of class \"simmap\" or",
      "\"multiSimmap\": argument recursive=TRUE ignored")
    warning(msg)
  }
  return(obj)
}

Let’s try it out.

library(phytools)
data(sunfish.tree)
data(sunfish.data)
feeding_mode<-setNames(sunfish.data$feeding.mode,
  rownames(sunfish.data))
ard_fit<-fitMk(sunfish.tree,feeding_mode,
  model="ARD",pi="fitzjohn")

First, combine "simmap" objects.

first_simmap<-simmap(ard_fit,nsim=1)
first_simmap

## 
## Phylogenetic tree with 28 tips and 27 internal nodes.
## 
## Tip labels:
## 	Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
## 	non, pisc
## 
## Rooted; includes branch lengths.

second_simmap<-simmap(ard_fit,nsim=1)
second_simmap

## 
## Phylogenetic tree with 28 tips and 27 internal nodes.
## 
## Tip labels:
## 	Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
## 	non, pisc
## 
## Rooted; includes branch lengths.

two_simmaps<-c(first_simmap,second_simmap)
two_simmaps

## 2 phylogenetic trees with mapped discrete characters

Now, let’s combine two "multiSimmap" objects. I imagine this is more combine because these could be, for example, results from two separate stochastic mapping analyses.

first_hundred<-simmap(ard_fit)
first_hundred

## 100 phylogenetic trees with mapped discrete characters

second_hundred<-simmap(ard_fit)
second_hundred

## 100 phylogenetic trees with mapped discrete characters

twohundred_simmap<-c(first_hundred,second_hundred)
twohundred_simmap

## 200 phylogenetic trees with mapped discrete characters

Neat.

Just as with ape c.phylo and c.multiPhylo, we should be able to combine objects of the two different classes.

abunchof_simmaps<-c(first_hundred,first_simmap,
  second_hundred,second_simmap)
abunchof_simmaps

## 202 phylogenetic trees with mapped discrete characters

Just for Twitter, let’s see what 202 phylogenies look like plotted.

layout(matrix(c(203,1:202,204),12,17))
par(bg="black")
cols<-setNames(viridisLite::viridis(n=2),
  levels(feeding_mode))
cols<-setNames(palette()[c(4,2)],levels(feeding_mode))
plot(abunchof_simmaps,cols,ftype="off",lwd=2)

This will be in phytools soon.