The Box-Einstein surface of mathematical models

As a mathematical modeller in evolutionary biology, my seminar bingo card has four prime boxes. Watching a talk about evolution, I count down the minutes to the first appearance of Dobzhansky’s “nothing in biology” quote (or some variant thereof) or a picture of Darwin’s “I think” sketch. For mathematical modelling, it’ll be either Albert Einstein or George Box:

“All models are wrong but some are useful” – George Box

“Everything should be made as simple as possible, but not simpler” – probably not Albert Einstein

Of course, such quotes are popular for good reason, and I’m not criticising those who use them to good effect, but all the same it can be fun to try to find a new way of presenting familiar material. That’s why in spring 2015 I came up with and tweeted a visual summary of the latter two aphorisms, which I named the Box-Einstein surface of mathematical models:

plot

The grey region in the plot ensures that all possible models have some degree of “wrongness”, but the contours in the remaining region tell us that some models are useful all the same. To find the most useful description of a particular phenomenon, we must reduce complexity without overly increasing wrongness.

A key thing to understand about this diagram is that although the boundary of the grey region is invariant, the surface is changeable. If our empirical knowledge of the system becomes richer, or if we change the scope of our enquiry, the most useful model may be more or less complex than before.

Einstein’s quote can be seen as simply paraphrasing Occam’s razor, but I think it has additional meaning with regard to what Artem Kaznatcheev calls heuristic and abstract mathematical models, such as are generally used in biology. In statistics, a simple model has few degrees of freedom, which is desirable to reduce overfitting. However, statisticians should also beware what JP Simmons and colleagues termed “researcher degrees of freedom”:

“In the course of collecting and analyzing data, researchers have many decisions to make: Should more data be collected? Should some observations be excluded? Which conditions should be combined and which ones compared? Which control variables should be considered? Should specific measures be combined or transformed or both?

“It is rare, and sometimes impractical, for researchers to make all these decisions beforehand. Rather, it is common (and accepted practice) for researchers to explore various analytic alternatives, to search for a combination that yields “statistical significance,” and to then report only what “worked.” The problem, of course, is that the likelihood of at least one (of many) analyses producing a falsely positive finding at the 5% level is necessarily greater than 5%.”

Likewise, when a researcher makes a mathematical model of a dynamical system – be it a set of differential equations or a stochastic agent-based model – he or she makes numerous decisions, usually with more or less full knowledge of the empirical data against which the model will be judged.

But there’s an important difference between the process of collecting data and that of creating a mathematical model. Ideally, the experimentalist can minimise researcher degrees of freedom by following a suitable experimental design and running controls that enable him or her to test a hypothesis against a null according to a predetermined statistical model. For most mathematical models there is no such template, and a process of trial and improvement is unavoidable, forgivable, and even desirable (inasmuch as it strengthens understanding of why the model works). The role of mathematical modeller is somewhere between experimentalist and pure mathematician. By making our models as simple as possible, we shift ourselves further toward the latter role, and our experimentation becomes less about exploiting our freedom and more about honing our argument.

For further reading, check out Artem Kaznatcheev’s insightful post about what “wrong” might mean, and why Box’s quote doesn’t necessarily apply to all types of model.

Visualizing evolutionary dynamics with ggmuller

A how-to guide for making beautiful plots of your data.

Readers of experimental evolution studies by Richard Lenski, Jeffrey Barrick and others, or of cancer evolution reviews, will have seen stacked area plots that combine information about frequency dynamics and phylogeny. The example below is from a recent preprint by Rohan Maddamsetti, Lenski and Barrick. The horizontal axis is time (measured in generations), and each coloured shape represents a genotype, with height corresponding to genotype frequency (or relative abundance) at the corresponding time point. Descendant genotypes are shown emerging from inside their parents.

Lenski

Such diagrams are sometimes termed Muller plots in honour of Hermann Joseph Muller (of ratchet fame), who used them in 1932 to illustrate an evolutionary advantage of sex. I wanted to draw Muller plots of some evolutionary simulations. Unable to find generic software, I spent a few days coding in R and tweeted an example result:

Tweet

Encouraged by the response, I expanded my code into the fully documented R package ggmuller (thanks to Sean Leonard for suggesting the name). Ggmuller is my attempt to set an easy-to-use standard for representing evolutionary dynamics. You can install the package from github with a few lines of R code:

install.packages("devtools") # enables installation from github
library(devtools)
install_github("robjohnnoble/ggmuller")
library(ggmuller)

Some of the code that follows depends on the latest version of ggmuller; if you installed the package before 20th August then please reinstall it for full functionality.

How does it work?

The core function is get_Muller_df, which constructs a specially ordered data frame from which we can draw a stacked area plot. get_Muller_df takes two data frames as input: first an adjacency matrix that defines a phylogeny (with columns named “Identity” and “Parent”), and second the populations of each genotype over time (with columns named “Identity”, “Generation” and “Population”). The genotype names in the “Identity” and “Parent” columns can be any words or numbers. You can use a phylo object instead of an adjacency matrix, provided the genotype names in the population data are integers that correspond to the node numbers in the phylo object.

get_Muller_df records a path through the phylogenetic tree. Each genotype appears exactly twice in the path: once before and once after all of the genotype’s descendants. The function then associates each instance of each genotype with half of the genotype’s frequency over time.

To draw plots, ggmuller exploits Hadley Wickham’s ggplot package (which accounts for the first part of its name). Because ggplot’s stacked area plot respects row order, and because the two instances of each genotypes are coloured identically, the resulting plot shows each genotype emerging from the centre of its parent (following Jeffrey Barrick’s example, multiple descendants are stacked from top to bottom in chronological order of appearance). Whereas I could have made descendants emerge at different heights – as in the Lenski lab figure above – it seems simpler and less ambiguous to have all descendants emerge exactly in the middle of their parents. I propose this as a standard to facilitate comparing plots of this type.

Basic usage

If your adjacency matrix and your population data frame are properly formatted, you can use ggmuller to visualize your results with just two lines of R code:

Muller_df <- get_Muller_df(example_edges, example_pop_df)
Muller_plot(Muller_df)

Simply replace example_edges and example_pop_df with the names of the data frames that contain your adjacency matrix and population data, respectively.

That’s all you need for basic usage. The rest of this post explains why these plots are such a powerful way of visualizing evolutionary dynamics, and introduces some additional features and options.

A couple of minimal examples

To get a better idea of why you might want to use ggmuller, let’s recreate Figure 1 of a recent commentary by Mark Robertson-Tessi & Sandy Anderson. The figure (drawn by Chandler Gatenbee) comprises two diagrams contrasting the “traditional” selective sweep model of cancer progression and the“Big Bang” neutral evolution model proposed by Andrea Sottoriva and colleagues.

Most of the code that follows is to create an appropriate data set. First we make an adjacency matrix that defines the phylogeny. In this minimal example, the gentoype names are integers, and each genotype has exactly one descendant (1 begets 2, which begets 3, which begets 4):

edges1 <- data.frame(Parent = 1:3, Identity = 2:4)

Next we construct a data frame of genotype populations over time. In the neutral evolution case, all genotypes (except the background genotype) have the same growth rate:

# a function for generating exponential growth curves:
pop_seq <- function(gens, lambda, start_gen) c(rep(0, start_gen),
 exp(lambda * gens[0:(length(gens) - start_gen)]))
gens <- 0:150 # generations
lambda <- 0.1 # baseline exponential growth rate
pop1 <- data.frame(Generation = rep(gens, 4),
 Identity = rep(1:4, each = length(gens)), 
 Population = c(1E2 * pop_seq(gens, lambda, 0), 
 pop_seq(gens, 2*lambda, 0), 
 pop_seq(gens, 2*lambda, 3), 
 pop_seq(gens, 2*lambda, 8)))

We use get_Muller_df to combine the phylogeny and population dynamics into a data frame we can plot from:

Muller_df1 <- get_Muller_df(edges1, pop1)

We create the plot using Muller_plot (which is a wrapper for ggplot’s stacked area plot) using a custom colour palette (optional):

my_palette <- c("white", "skyblue2", "red", "green3")
plot1 <- Muller_plot(Muller_df1, palette = my_palette)

We proceed similarly for the second panel, in which genotypes have different growth rates (for example, due to the accumulation of advantageous mutations):

edges2 <- data.frame(Parent = 1:3, 
 Identity = 2:4)
pop2 <- data.frame(Generation = rep(gens, 4),
 Identity = rep(1:4, each = length(gens)),
 Population = c(1E2 * pop_seq(gens, lambda, 0), 
 pop_seq(gens, 2*lambda, 0), 
 pop_seq(gens, 4*lambda, 40), 
 pop_seq(gens, 8*lambda, 80)))
Muller_df2 <- get_Muller_df(edges2, pop2)
plot2 <- Muller_plot(Muller_df2, palette = my_palette)

Finally, we draw the plots together using the gridExtra package:

library(gridExtra)
grid.arrange(plot1, plot2)

Simple

And here we see the power of the Muller vizualisation: even though these systems have identical phylogenies, we can see at a glance that they have very different evolutionary dynamics.

A more sophisticated phylogeny

Now let’s look at a more interesting, branched phylogeny:

edges3 <- data.frame(Parent = paste0("clone_", 
 LETTERS[c(rep(1:3, each = 2), 2, 5)]), 
 Identity = paste0("clone_", LETTERS[2:9]))

The population data are generated using diverse fitness values:

fitnesses <- c(1, 2, 2.2, 2.5, 3, 3.2, 3.5, 3.5, 3.8)
pop3 <- data.frame(Generation = rep(gens, 9),
 Identity = paste0("clone_", LETTERS[rep(1:9, each = length(gens))]),
 Population = c(1E2 * pop_seq(gens, fitnesses[1]*lambda, 0), 
 pop_seq(gens, fitnesses[2]*lambda, 0), 
 pop_seq(gens, fitnesses[3]*lambda, 10), 
 pop_seq(gens, fitnesses[4]*lambda, 20),
 pop_seq(gens, fitnesses[5]*lambda, 30),
 pop_seq(gens, fitnesses[6]*lambda, 40),
 pop_seq(gens, fitnesses[7]*lambda, 50),
 pop_seq(gens, fitnesses[8]*lambda, 50),
 pop_seq(gens, fitnesses[9]*lambda, 60)),
 Fitness = rep(fitnesses, each = length(gens)))

Of course there are more efficient ways of coding this in R, but perhaps the long-winded version is clearer for readers unfamiliar with lapply and its siblings. Note the inclusion of an optional “Fitness” column containing the fitness values. We combine the data using get_Muller_df as before:

Muller_df3 <- get_Muller_df(edges3, pop3)

Finally, we create two versions of the plot, first coloured by genotype identity (the default) and then (using the RColorBrewer package) coloured by fitness, making use of that optional column:

plot3 <- Muller_plot(Muller_df3, add_legend = TRUE)
library(RColorBrewer)
num_cols <- length(unique(Muller_df3$Fitness)) + 1
Muller_df3$Fitness <- as.factor(Muller_df3$Fitness)
plot3a <- Muller_plot(Muller_df3, colour_by = "Fitness", 
 palette = rev(colorRampPalette(brewer.pal(9, "YlOrRd"))(num_cols)), 
 add_legend = TRUE)
grid.arrange(plot3, plot3a)

Branching

 Optional extras

Since simulations can result in very large sets of genotypes, many of which never become abundant, get_Muller_df has a threshold option to exclude rarities. By default, the data frame includes only rows with nonzero population; this can be overridden by setting add_zeroes = TRUE, in which case all genotypes are represented at every generation (which makes it easier to add new columns, for example).

Converting between adjacency matrix and phylo representations

We can visualize the tree from our branched phylogeny example using Emmanuel Paradis’ ape package, but first we need to use the ggmuller function adj_matrix_to_tree to convert our adjacency matrix (with arbitrary node names) to a phylo object that ape can understand. Phylo is the standard way of representing trees in R, and it requires nodes to be numbered in a particular order.

tree <- adj_matrix_to_tree(edges3)

After optionally labeling the nodes, we use ape’s plot.phylo function to draw the tree:

library(ape)
tree$tip.label <- 1:length(tree$tip.label) # optional
tree$node.label <- (length(tree$tip.label) + 1):10 # optional
plot(tree, show.node.label = TRUE, show.tip.label = TRUE, tip.color = "red")

Tree

If you already have a tree in phylo format, you can use that instead of an adjacency matrix as input to get_Muller_df, as long as the genotype names in the population data are integers that correspond to the node numbers in the phylo object.

Showing changes in population size

Although ggmuller deals in frequencies, you can hack it to show changes in the total population size by adding a dummy common ancestor to your phylogeny, with population equal to the maximum total population size minus the total population size at each time point, and then putting a new scale on the y-axis. By default, the dummy common ancestor (which corresponds to unfilled space) will be shown in black, whereas actual genotypes will have paler colours. I may automate this process in future versions.

Next steps

I hope ggmuller will prove useful to researchers in experimental evolution, cancer research, population genetics, and beyond. I plan to submit the package to CRAN and will continue to add functionality. Meanwhile, if you use ggmuller in a publication, please cite the github version and notify me via email or twitter so I can keep track of uptake. Bug reports and feature requests would be most welcome.

[Edited on 5th September 2016 to fix a typo in the code defining “edges3”.]