R Hackathon 1/Ancestral State Reconstruction
The primary ancestral state reconstruction algorithms in R are accessed through the function "ace" in package "ape". To work through the following example using the Geospiza dataset, make sure that you have installed and loaded ape into your R session and loaded the Geospiza phylogeny and tip data into memory.
library(ape) geotree <- read.nexus("geospiza.nex") geodata <- read.table("geospiza.txt")
Tip data is not available for the outgroup "olivacea", so drop that taxon from the analysis.
geotree <- drop.tip(geotree, "olivacea")
IMPORTANT. The row.names of your dataframe must match the tip.labels of your phylogeny. In the example above, the row.names for geodata will match the tip.labels for geotree after "olivacea" has been culled. However, the individual columns of geodata (e.g. geodata$wingL) do not automatically have the row.names of the whole data table associated with them! If you call a column of the data.table with ace without first dealing with this issue, the analysis will run, but the tip data will be disassociated from the proper tips! There are two workarounds.
One option is to sort the datatable so that the taxa appear in the same order in the table as in the phylogeny. For example:
geodata <- geodata[geotree$tip.label, ]
The other, and better option is to extract the column of data of interest from the overall datatable, creating a new vector, and to transfer the row.names of the datatable to the NAMES of the vector. This is the option least likely to inadverently disassociate the tip data from the tips.
wingL <- geodata$wingL names(wingL) <- row.names(geodata)
The worked examples below assume that you have used the second solution (vector extraction and name assignment).
Reconstructing Ancestral States for Continuous Variables
There are three general options for continuous variables: reconstructions based on maximum likelihood (ML), reconstructions based upon phylogenetic independent contrasts (pic) and reconstructions based on generalized least squares (GLS). Be aware that the generalized least squares algorithms seem to give spurious results near the root of the phylogeny at present.
This syntax will reconstruct the ancestral states for the variable "wingL" (extracted from geodata) using the Brownian motion-based maximum likelihood (ML) estimator of Schluter et. al. (1997). This is the default method.
MLreconstruction <- ace(wingL, geotree, type="continuous", method="ML")
Why am I getting all these warning messages?
ML reconstruction using ace tends to generate a large number of not-a-number error messages. These result when the program calculates the likelihood of particularly poor fits. You can safely ignore these messages.
PHYLOGENETIC INDEPENDENT CONSTRASTS
This syntax will reconstruct the ancestral states for the variable "wingL" (extracted from geodata) using Felsenstein's (1985) phylogenetic independent contrasts (pic). This is also a Brownian-motion based estimator, but it only takes descendants of each node into account in reconstructing the state at that node. More basal nodes are ignored.
picreconstruction <- ace(wingL, geotree, type="continuous", method="pic")
GENERALIZED LEAST SQUARES
The generalized least squares algorithm in ace seems to given spurious results near the root of the phylogeny, at least at present. Specifically, the root node is always reconstructed at possessing state zero, and internal nodes near the root may receive reconstructed values out of the range of values observed among the tips.
Reconstructions based upon generalized least squares require specification of a correlation structure for the generalized linear model. There are three basic options for the specification of the correlation structure: 1) corBrownian, which uses a simple Brownian motion model, corMartins, which uses an Ornstein-Uhlenbeck (constrained random-walk) model, and corGrafen, which is a modified Brownian motion model (Grafen 1989) (ADD MORE INFO ON WHAT EXACTLY THE GRAFEN MODEL DOES).
This is the syntax for the simple Brownian model
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corBrownian(1, geotree))
The syntax for the Ornstein-Uhlenbeck model requires specifying an alpha parameter (here, 0.5).
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corMartins(0.5, geotree))
The syntax for the Grafen model requires specifying a rho parameter (here, 1) which exponentiates the recalculated branch lengths.
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corGrafen(1, geotree))
WHAT ABOUT SQUARED CHANGE PARSIMONY?
Squared change parsimony is mathematically equivalent to a special case of Schluter et. al.'s maximum likelihood method in which branch length information is ignored. Thus, to perform a squared change parsimony reconstruction, set all branch lengths equal to 1 and calculate a maximum likelihood solution.
geotreeones <- compute.brlen(geotree, 1) SQPreconstruction <- ace(wingL, geotreeones, type="continuous", method="ML")
HOW DO I USE THE OUTPUT FROM ACE?
Each of the objects (for example: MLreconstruction) that we reconstructed above is a list of several elements. These are loglik (the log-likelihood of the most likely reconstruction), ace (the vector of reconstructed node values), sigma2 (a two element vector of the rate estimate and its standard error) and CI95, the 95% confidence intervals around the vector of reconstructed node values. Thus
returns the vector of node values for wing length (wingL) that were reconstructed using maximum likelihood, and
wingLfinal <- c(wingL, MLreconstruction$ace)
concatenates the original tip data with the reconstructed node data into a single vector.
HOW DO I PLOT THE OUTPUT FROM ACE? (drawn from a course handout by Gene Hunt)
Plot.phylo can scale symbols at tips and nodes of your tree to your character data. Try this to visualize the evolution of wing length (wingL) across the Geospiza phylogeny.
plot.phylo(geotree) tiplabels(pch = 21, cex=wingL) nodelabels(pch = 21, cex=GLSreconstruction1$ace) ## (pch = 21 is just telling plot.phylo which symbol to use)
Note that the blue internal nodes are all similar in size to the yellow tip labels. This is because there is little variation in wing length across the dataset.
HOW DO I FIT A MODEL OF CHARACTER CHANGE TO CONTINUOUS DATA?
While ace returns a rate of Brownian evolution (e.g. MLreconstruction$sigms2), that rate and its associated likelihood vale (e.g. MLreconstruction$loglik) is conditioned on the specific ancestral states reconstructed by ace. More generally applicable model fits using Brownian motion and and the Ornstein-Uhlenbeck model are available in the packages ouch and geiger. Please see here for more information.
Reconstructing Ancestral States for Discrete Variables
The ace command in ape can reconstruct ancestral states for discrete characters using maximum likelihood.
First, ensure that you have ape loaded.
Load or reload the Geospiza phylogeny, which will restore the taxon "olivacea" that you may have removed in the continuous example above.
geotree <- read.nexus("geospiza.nex")
These commands will set up a new discrete character vector "char1" and assign the names of the terminal taxa in the Geospiza phylogeny to the elements in the new vector. This is same vector used in the example on modeling discrete character evolution in R.
You will need to choose among several different options for the matrix that specifies the transition probabilitites between the states of your discrete character. Ace offers shortcuts for a one-parameter equal rates model (ER), a symmetric model (SYM) in which forwards and reverse transitions between states are constrained to be equal, and an all rates different matrix (ARD) where all possible transitions between states receive distinct parameters. You can also specify your own matrix.
Let's reconstruct ancestral states for the Geospiza data using the three standard methods (ER, SYM and ARD).
ERreconstruction <- ace(char1, geotree, type="discrete", model="ER") SYMreconstruction <- ace(char1, geotree, type="discrete", model="SYM") ARDreconstruction <- ace(char1, geotree, type="discrete", model="ARD")
These commands issue lots of output, of which you are likely to be most interested in the log likelihoods and the inferred transition rates. You can access the log likelihoods for each model as follows
ERreconstruction$loglik SYMreconstruction$loglik ARDreconstruction$loglik
You should receive values of -9.00, -9.00 and -7.28, respectively.
The transition rates are accessed with
ERreconstruction$rates SYMreconstruction$rates ARDreconstruction$rates
ER and SYM return a single rate of 51.6, while ARD returns a forward and reverse rate (6.00 and 17.9, respectively).
You may have noticed that the output for the ER model matches that of the SYM model. This is not surprising. For a binary character these models are identical, though they are distinct for characters with three or more states. For a three-state character, ER is a one parameter model, SYM a three parameter model, and ARD a six parameter model.
Which model should you prefer?
In this example, the ARD model gives the highest likelihood. However, it also includes more parameters that the ER and SYM models, and it is well known that adding parameters to models generally increase likelihoods. To determine whether the use of the more heavily parameterized mode is appropriate, you should execute a likelihood ratio test.
Twice the difference in log likelihood between two models (known as the G statistic) is distributed as chi square, with degrees of freedom equal to the number of parameters that differ between the models. Thus a likelihood test of these models asks whether the difference in likelihoods is large enough to lie in the rightmost tail of the chisquare distribution, such as the largest 5% of values.
The command pchisq(value, df) gives the percentage of the cumulative distribution function for chisquare that lies to the left of the given value for a desired degree of freedom. Thus, the following command gives the percentage of the chisquare distribution that lies to the right of the observed likelihood difference for the ER and ARD models, given the Geospiza data.
1-pchisq(2*abs(ERreconstruction$loglik - ARDreconstruction$loglik), 1)
Evaluation of the equation above yields the p-value for the likelihood test, in this case 0.064. This is not a significant result at a 5% error rate (though it is very close to significance), so we do not accept the more heavily parameterized ARD model, and instead choose the single-parameter ER model.
For a single degree of freedom, the significance cutoff for the likelihood values lies at about 1.93 log-likelihood units.
Just the absolute difference in log likelihoods.
(abs(ERreconstruction$loglik) - abs(ARDreconstruction$loglik)
How do I define a custom transition model.
Defining a custom four-parameter matrix for a three state character.
transitions <- matrix (c(0, 1, 2, 1, 0, 3, 2, 4, 0), nrow=3)
Matrix should look like.
[,1] [,2] [,3] [1,] 0 1 2 [2,] 1 0 4 [3,] 2 3 0