likelihood reconstruction of ancestral states

In a previous post, I talked about how summarizing the state at each internal node over many make.simmap mappings did not correspond exactly with the ace reconstructions.  I contacted Liam Revell about this, and he informed me that this may be because ace does not compute the scaled marginal likelihoods but the conditional likelihoods of the subtrees descending from each node.  He suggested I try rerooting the tree at each internal node and using the ace reconstruction at the root to find the marginal likelihoods.  Here is my attempt:

nodes = (exampleTree$Nnode+2):(exampleTree$Nnode*2+1) # This gives me the number associated with each internal node

reRootAnc = t(sapply(nodes, function(x) {
  tr = reroot(exampleTree, node = x, position = 0) # rerooting the tree at each internal node
  reconst = ace(x = discrete, phy = tr, type = "discrete", model = "SYM") # estimating the maximum likelihood ancestral state estimate
  reconst$lik.anc[1,] # taking only the value for the root
}))

This is the result I got for an example tree:

Seems a little strange to me...I wonder what might be going on.

I knew this shouldn't be so complicated...

I was trying to analyze my incomplete dataset, and I needed to remove data for species where I have measured fewer than four individuals.  Since my dataset was small enough, it was easy to just remove them by hand, doing something like this:

> df
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
5        b    7
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1
11       d    3
12       d    2

> dfremoved = df[-which(df$Species=="b"),]
> dfremoved = dfremoved[-which(dfremoved$Species=="d"),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1

But in order to do this systematically, I used a couple of steps.

> toofew = names(which(table(df$Species) < 4))
> toofew
[1] "b" "d"

First, I found the species names for species with fewer than four individuals.  With this, I can remove all rows where df$Species match any of these names.

> dfremoved = df[!(df$Species %in% toofew),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1

Avoiding Repetition

If there is anything I learned in STA 141, I learned the importance of the DRY principle:  don't repeat yourself.  Anything repetitive was heavily penalized in the grades, but it's also more prone to error and often takes much longer.  Even still, it's easy to be lazy and fall into the "copy/paste then change one word" routine, especially when just exploring a data set.  That's what I started out doing, but it turns out that doing it the 'right' way even easier!

For example, I want a plot of all morphological variables against size (standard length).  I want the points colored by species, and I want each point to be a unique value for the species.

plot(dat$standard.length, dat$head.length, col = as.factor(dat$Species), pch = as.character(dat$Number))

Now I can copy/paste this line and replace "head.length" with all of my other variables.  Simple enough.

But it turns out I have 24 variables.  So doing this will take much longer than this simple loop:

meas = names(dat)[9:length(dat)] # all morphological variables except standard length

sapply(meas, function(x) {
  png(file = paste(x, ".png", sep = ""))
  plot(dat$standard.length, dat[,x], col = as.factor(dat$Species), pch = as.character(dat$Number))
  dev.off()
})

I can do better by putting down axis labels and such, but now I have .png files of each of my morphological variables that I can browse through with my favorite image viewer.

Control Flow

I use if/else and for quite often, but rarely use while or repeat.  Even still, there are a couple of things with if/else that give me trouble if I haven't used it in a while.

If everything is on one line, all is good:

> x = 1

> y = if(x>0) 1 else 0

> y

[1] 1

But if I put the curly braces in the wrong place:

> if(x>0) {
+   x = -x
+   y = 1
+ }
> else {
Error: unexpected 'else' in "else"
>   x = x
>   y = 0
> }
Error: unexpected '}' in "}"

So I have to always make sure to do this:

> if(x>0) {

+   x = -x

+   y = 1

+ } else {

+   x = x

+   y = 0

+ }

> x

[1] -1

> y

[1] 1

No polytomies allowed?

I have recently been in a position to want an efficient way to obtain a binary Newick-format tree from a Newick-format tree that may or may not have polytomies.  Using a few functions from ape, this was fairly simple to obtain:


makebinary = function(newick) {
  tree = read.tree(text = newick)
  if(is.binary.tree(tree)) {
    return(newick)
  } else {
    return(write.tree(multi2di(tree)))
  }
}


This contains three very useful functions from ape:  read.tree, write.tree, is.binary.tree, and multi2di.  The first two are used to read/write Newick-format trees.  The functions read.nexus and write.nexus can be used for Nexus-format trees.  read.tree can be used for files as well.  is.binary.tree checks whether a tree has any polytomies, and multi2di converts a tree with multichotomies to a fully dichotomous tree with some branches of length 0.  There is also a function di2multi that will collapse any branches less than a tolerance level to a polytomy.

Animating SIMMAP trees

It took me quite a while to get this working, but I finally did:

The tree topology is from Rüber et al. (2004).  This is what I used to generate the image:

animateSimmap = function(phy, interval = .02, numb = length(phy), name = "animation.gif", ...) {
  if(class(phy) != "multiPhylo") stop("object must be multiPhylo")
  saveGIF(for(i in 1:numb) {
            dev.hold()
            plotSimmap(phy[[i]], ...)
            Sys.sleep(interval)
          }, movie.name = name)
}

So this function takes in a multiPhylo object that has stochastic character mappings, the time interval between trees (I'm not sure if this does anything because saveGIF might just have a set interval already), the number of mappings to animate, the name of the file to create, and any arguments you want to give to the plotSimmap function.  In terms of the mappings shown above, they're meaningless because the branches aren't proportional to anything meaningful, but it shows how this function would work.

plotSimmap is from the package phytools, and saveGIF is from the package animation.  To use saveGIF, you also need either ImageMagick (which is what I used), GraphicsMagick, or LyX.  

SIMMAP Trees

Thanks to Liam Revell, we can now produce simulated discrete character mappings within R.  The function make.simmap uses ace, from the package ape, to fit the model to use for the simulations.  Thus, the proportion at which a character state appears at a node among many iterations of simulations should be roughly equivalent to the likelihood of that state as estimated in ace.  I have heard this from many people, but I wanted to be able to summarize the actual simulated states.  In order to do this, I took a look at the components of a SIMMAP tree.  I'll illustrate this with a mock example:

We can take a look at what components this object has:

> names(exampleSimmap)
[1] "edge"        "edge.length" "tip.label"   "Nnode"       "maps"       
[6] "mapped.edge"

The first components are the same as any phylo object.  So maps and mapped.edge are what make a SIMMAP tree special.  Let's take a look (the middle elements removed to save space):

> exampleSimmap$maps
[[1]]
        1         0 
0.0958356 0.3197055 
[[2]]
        0         1 
0.2371619 0.1659281 
[[3]]
         1 
0.01380231 
...
[[17]]
        1 
0.1164989 
[[18]]
        1         0 
0.1161601 0.1769559 

This is precisely what we need!  Each element of exampleSimmap$maps represents a single branch, and the values represents the length of time that branch spends in each state, in this case 0 or 1.  That means we can simply take each branch's starting value (whatever is the name of the first element of that branch), and that is the value at the node where the branch starts.  Let's see if we can find this.

> exampleSimmap$edge
      [,1] [,2]
 [1,]   11   12
 [2,]   12   19
 [3,]   19    1
 [4,]   19    2
 [5,]   12   13
 [6,]   13   17
 [7,]   17    3
 [8,]   17    4
 [9,]   13   18
[10,]   18    5
[11,]   18    6
[12,]   11   14
[13,]   14   15
[14,]   15    7
[15,]   15   16
[16,]   16    8
[17,]   16    9
[18,]   14   10

The element named edge gives us the starting and ending node for each of the 18 edges in our tree.  That means we can use to figure out which node corresponds with which state.  Here I've written a function that takes in a SIMMAP tree and returns a named vector where the values are the node states and the names are the nodes.

mappedNode = function(phy) {
  # phy must be a SIMMAP tree
  nodes = phy$edge[,1] # this gives us the starting node for all edges
  map = sapply(phy$maps, function(x) attr(x, "names")[1]) # this gives us the starting value of each branch
  df = unique(data.frame(nodes = nodes, map = map)) # here we're removing the repeated values as interior nodes will have multiple branches
  mapping = df$map
  names(mapping) = df$nodes # naming the vector with node names
  mapping
}

Now we can run this function over all of the simulated mappings:

mappings = sapply(simmapTrees, function(x) mappedNode(x))
# change the 0 and 1 to numeric
nummaps = as.data.frame(sapply(1:length(simmapTrees), function(x) as.numeric(mappings[,x])))
# make sure the row names correspond to node names
rownames(nummaps) = rownames(mappings)
# get the number of simulated trees with  node state of 1
sums = sapply(rownames(nummaps), function(x) sum(nummaps[x,]))
# change that to a frequency
freq = sums/length(simmapTrees)
# plot
plot(exampleTree, label.offset = .05)
nodelabels(pie = freq, cex = .65, node = as.numeric(names(freq)))
tiplabels(pie = discTrait, cex = .65)
Looks fairly reasonable.  Now let's estimate the ancestral character using ace:

MLACE = ace(discTrait, exampleTree, type = "discrete", model = "SYM")

plot(exampleTree, label.offset = .05)
nodelabels(pie = 1 - MLACE$lik.anc, cex = .65)
tiplabels(pie = discTrait, cex = .65)

Pretty close, but not identical.  Maybe if we do more simulated mappings (I did 1000 here), they will start to look more similar?