Given a tree and a cladistic-type matrix uses either likelihood ratio tests or the Akaike Information Criterion to compare rate models across branches, clades, time bins, or character partitions.
test_rates(
time_tree,
cladistic_matrix,
time_bins,
branch_partitions = NULL,
character_partitions = NULL,
clade_partitions = NULL,
time_partitions = NULL,
change_times = "random",
test_type = "aic",
alpha = 0.01,
multiple_comparison_correction = "benjaminihochberg",
polymorphism_state = "missing",
uncertainty_state = "missing",
inapplicable_state = "missing",
time_binning_approach = "lloyd",
all_weights_integers = FALSE,
estimate_all_nodes = FALSE,
estimate_tip_values = FALSE,
inapplicables_as_missing = FALSE,
polymorphism_behaviour = "equalp",
uncertainty_behaviour = "equalp",
threshold = 0.01,
all_missing_allowed = FALSE
)The time binning used (NB: May be slightly altered from the input values).
Matrix of inferred character changes.
The global (mean) character rate in changes per million years. I.e, the average rate across all characters, branches and time bins, effectively the null hypothesis for any test performed.
Whether or not continuous characters were converted to discrete characters (important for handling the data in downstream analys(es)).
List of branch partition results (corresponding to branch_partitions. NULL if not requested.
List of character partition results (corresponding to character_partitions. NULL if not requested.
List of clade partition results (corresponding to clade_partitions. NULL if not requested.
List of time bin partition results (corresponding to time_partitions. NULL if not requested.
Matrix showing calculated rates for each branch. NULL if branch_partitions is not requested.
Matrix showing calculated rates for each character. NULL if character_partitions is not requested.
Matrix showing calculated rates for each clade. NULL if clade_partitions is not requested.
Matrix showing calculated rates for each time bin. NULL if time_partitions is not requested.
The time-scaled input tree used as input (provided as output for use with visualisation functions).
A tree (phylo object) with branch durations that represents the relationships of the taxa in cladistic_matrix.
A character-taxon matrix in the format imported by read_nexus_matrix.
An object of class timeBins.
A list of branch(es) (edge number) partitions to test as N-rate parameter model (where N is the total number of partitions). If NULL (the default) then no partition test(s) will be made.
A list of character partition(s) (character numbers) to test as N-rate parameter model (where N is the total number of partitions). If NULL (the default) then no partition test(s) will be made.
A list of clade partition(s) (node numbers) to test as N-rate parameter model (where N is the total number of partitions). If NULL (the default) then no partition test(s) will be made.
A list of time bin partition(s) (numbered 1 to N) to test as N-rate parameter model (where N is the total number of partitions). If NULL (the default) then no partition test(s) will be made.
The time at which to record the character changes. One of "midpoint" (changes occur at the midpoint of the branch), "spaced" (changes equally spaced along branch), or "random" (change times drawn at random from a uniform distribution; the default and recommended option). Note: this is only meaningful if testing for time bin partitions.
Whether to apply an Akaike Information Criterion ("aic"; the default) or likelihood ratio test ("lrt").
The alpha value to be used for the significance tests. The default is 0.01. Thsi is only relevant if using likelihood ratio tests.
Current options are: 1. "benjaminihochberg" (the Benjamini and Hochberg 1995 false discovery rate approach; default and recommended), or 2. "bonferroni" (the Bonferroni correction). This is only relevant if using likelihood ratio tests.
Current options are: 1. "missing" (converts polymorphic values to NA; the default), or 2. "random" (picks one of the possible polymorphic states at random).
Current options are: 1. "missing" (converts uncertain values to NA; the default), or 2. "random" (picks one of the possible uncertain states at random).
The only current option is "missing" (converts value to NA).
One of "close" or "lloyd" (the latter is the default and recommended option).
Logical for whether (TRUE) to reweight non-integer weights until all weights are integers or to leave them as they are (FALSE; the default).
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Option passed to estimate_ancestral_states.
Graeme T. Lloyd graemetlloyd@gmail.com and Steve C. Wang scwang@swarthmore.edu
Introduction
Morphological change can be captured by discrete characters and their evolution modelled as occurring along the branches of a phylogenetic tree. This function takes as primary input a character-taxon matrix of discrete characters (in the format imported by read_nexus_matrix) and a time-scaled phylogenetic tree (in the format of paleotree or strap) and begins by inferring ancestral states at the tree's internal nodes using the estimate_ancestral_states function. From here changes along individual branches can be estimated (only the minimum number of changes are inferred; see map_stochastic_changes for an alternative but unfinished approach) and hence rates can be calculated.
A discrete character rate can be expressed as the mean number of changes per million years (users may wish to normalise this by the number of characters for interpretation) and can be calculated for a branch (edge) of the tree, a clade (a mean rate for the edges descended from a single node), a character partition (the mean rate for a subset of the characters across all edges), or, most complex (see Lloyd 2016), the mean rate across the edges (or parts of edges) present in a time bin (defined by two values denoting the beginning and end of the time bin). In an ideal scenario these rates could be compared at face value, but that would require a large number of characters and very minimal (or zero) missing data. I.e., at an extreme of missing data if only one character can be observed along a branch it will either change (the maximum possible inferrable rate of evolution) or it will not (the minimum possible inferrable rate of evolution). In such cases it would be unwise to consider either outcome as being a significant departure from the mean rate.
Because of these complications Lloyd et al. (2012) introduced tests by which the significance of an edge (or other partitioning of the data, i.e., a clade, time bin etc.) could be considered to be significantly high or low in comparison to the mean rate for the whole tree (i.e., whether a two-rate model could be considered more likely than a one-rate model). This is achieved through a likelihood ratio test:
$$LR = value of likehood function under the null (one-rate) hypothesis / maximum possible value of likehood function under the alternative (two-rate) hypotheses$$
Typically we might expect the two hypotheses to be well defined a priori. E.g., an expectation that a specific branch of the tree might have a higher or lower rate than background due to some evolutionary shift. However, Lloyd et al. (2012) instead provided an exploratory approach whereby every possible one edge value was compared with the rate for the rest of the tree (and the equivalent with clades and time bins). This was the default in Claddis up to version 0.2, but this has now been replaced (since version 0.3) with a more customisable set of options that allows different types of hypotheses (e.g., partitioning the data by character), as well as more complex hypotheses (e.g., a three-rate model), to be tested. Since version 0.4 the option to replace likelihood ratio tests with the Akaike Information Criterion has also been added.
The four types of rate hypothesis
Following Cloutier (1991), Lloyd (2016) extended the two main types of rate hypotheses to four:
A branch rate (available here with the branch_partitions option).
A clade rate (available here with the clade_partitions option).
A time bin rate (available here with the time_partitions option).
A character partition rate (available here with the character_partitions option).
In Claddis (>=0.3) these partitions are defined as a list of lists of vectors where only the first N - 1 partitions need be defined. E.g., if comparing the first edge value (based on ape numbering, i.e., plot(tree); edgelabels()) to the rest of the tree then the user only needs to define the value "1" and the function will automatically add a second partition containing all other edges. This can be set with the option branch_partitions = list(list(1)). Similarly, to do what Lloyd et al. (2012) did and repeat the test for every edge in the tree (and assuming this variable is already named "tree") you could use, branch_partitions = lapply(X = as.list(x = 1:nrow(tree$edge)), as.list).
Because of the flexibility of this function the user can define any set of edges. For example, they could test whether terminal branches have a different rate from internal branches with branch_partitions = list(list(match(1:ape::Ntip(phy = tree), tree$edge[, 2]))). The clade_partitions is really just a special subset of this type of hypothesis, but with edges being defined as descending from a specific internal node in the tree. Once again, an exploratory approach like that of Lloyd et al. (2012) can be used with: clade_partitions = lapply(X = as.list(x = ape::Ntip(phy = tree) + (2:Nnode(tree))), as.list). Note that this excludes the root node as this would define a single partition and hence would represent the null hypothesis (a single rate model for the whole tree). (If using test_type = "aic" then the user typically will want a value for a single partition.) More generally clades must be defined by the node numbers they correspond to. In R an easy way to identify these is with: plot(tree); nodelabels().
Time bin partitions are defined in a similar way, but are numbered 1:N starting from the oldest time bin. So if wanting to do an exploratory test of single bin partitions (and only four time bins were specified) you could use: time_partitions = lapply(X = as.list(x = 1:4), as.list). Bins can be combined too, just as edges are above. For example, time bins 1 and 2 could form a single partition with: time_partitions = list(list(1:2)). Or if looking to test a model where each bin has its' own rate value you could use: time_partitions = list(as.list(x = 1:3)). Note, as before we do not need to specify the fourth bin as this will be automatically done by the function, however, time_partitions = list(as.list(x = 1:4)) will also work. Some caution needs to be applied with N-rate models (where N is three or larger) and test_type = "lrt" as a result favouring such models does not necessarily endorse N-separate rates. I.e., it could simply be that one bin has such a large excursion that overall the N-rate model fits better than the 1-rate model, but some 2-rate models might be better still. It is up to the user to check this themselves by exploring smaller combinations of bins and more genrally if exploring partitions of three or more use of the Akaike Information Criterion (test_type = "aic") is recommended.
Finally, character partitions allow the user to explore whether rates vary across different character types (numbers), e.g., skeletal characters versus soft tissue characters, or cranial characters versus postcranial characters. Here characters are simply numbered 1:N (across all blocks of a matrix), but single character partitions are less likely to be of interest. As an example of use lets say the first ten characters are what we are interested in as a partition (the second partition being the remaining characters), we could use: character_partitions = list(list(1:10)) to test for a two-rate model with test_type = "lrt".
Note that the list of lists structure is critical to defining partitions as it allows them to be of different sizes and numbers. For example, one partition of three and another of six, or one set of two partitions and another set of four partitions - structures not easily realizable using vectors or matrices. However, it may not be intuitive to some users so it is recommended that the user refers to the examples above as a guide.
Additionally, it should be noted that the user can test multiple types of hypotheses simultaneously with the function. For example, performing several branch tests whilst also performing clade tests. However, it is not necessary to perform all types simultaneously (as was the case up to version 0.2) and unused partition types can be set to NULL, the default in each case.
AIC vs LRT
Since Claddis version 0.4 the option to use the Akaike Information Criterion (AIC) instead of likelihood ratio tests (LRTs) has been added (although it was not properly functional until version 0.4.2). Practically speaking the AIC uses something similar to the denominator term from the LRT (see equation above) and adds a penalty for the number of parameters (partitions). However, it also fundamentally alters the comparative framework applied and hence needs more careful attention from the user to be applied correctly. Specifically, the LRT is by its nature comparative, always comparing an N-rate partition with a one-rate partition. By contrast the AIC does not directly apply any comparison and so the user must logically supply multiple partitionings of the data in order for the results to be meaningful. It might be assumed that a user will always want to apply a single partition that pools all the data for each type of test, whether this is all the edges (branches), time bins, or characters. This will thus be the obvious comparator for any multiple partition supplied, ensuring that any more complex partitioning is minimally superior to this. Additionally, it is also logical to consider each possible way of joining partitions simpler than the most complex partition being considered. E.g., if considering a four-partition model then the user should also consider all possible three-partition and two-partition combinations of that four-partition model. This can obviously lead to some complexity in supplying partitions on the user's part and so some automating of this process is planned in future (but is not fully available yet). For an example, see partition_time_bins.
Additionally, AIC values are not simple to compare as there is no direct equivalent of the alpha value from the LRT. Instead the user can modify the AIC values returned themselves to get delta-AIC or Akaike weight values (e.g., with geiger::aicw). (NB: I will not explain these here as there are better explanations online.) Furthermore, since version 0.4.2 sample-size corrected AIC (AICc) is also available in the output. Note that some caution should be used in applying this if the number of partitions is equal to the sample size or only one fewer. E.g., if you have ten time bins and ten partitions you can get a negative value due to the denominator term in the AICc calculation. Thus it is advised to use the raw AIC values as first approximation and be wary if the AICc flips the preferred model to a more complex model or models (i.e., those with more partitions) as this is the opposite of the intent of the AICc.
High versus low rates
Prior to Claddis version 0.3, rate results were classified as significantly high or significntly low as part of the output. This was done simply on whether the estimated p-value fell above or below the corrected alpha level (the significance threshold) and whether the first part of the two-rate partition had a higher or lower rate than the second part (the pooling of all other values). This simple interpretation is no longer valid here as the function can consider more than two partitions (high versus low is meaningless) and the allowing of AIC values means a significance test need not be performed either. Although the same interepretation can still be applied manually when only using two-partition tests and the LRT, other partition sizes and use of the AIC complicate this interpretation (in the same way an ANOVA is more complex than a two-sample t-test). This will also affect visualisation of the data (i.e. the simple pie chart coluring of non-significant, significntly high or significantly low rates seeen since Lloyd et al. 2012 will no longer apply). Instead the user should isolate the best model(s) and attempt to visualise these, perhaps using something like a heat map, with the mean rate for each partition being represented by an appropriate colour.
Other options
Since Claddis version 0.3 this function has allowed the user greater control with many more options than were offered previously and these should be considered carefully before running any tests.
Firstly, the user can pick an option for change_times which sets the times character changes are inferred to occur. This is only relevant when the user is performing time bin partition test(s) as this requires some inference to be made about when changes occur on branches that may span multiple time bins. The current options are: "midpoint" (all changes are inferred to occur midway along the branch, effectively mimicking the approach of Ruta et al. 2006), "spaced" (all changes are inferred to occur equally spaced along the branch, with changes occurring in character number order), or "random" (changes are assigned a random time by drawing from a uniform distribution between the beginning and end of each branch). The first of these is likely to lead to unrealistically "clumped" changes and by extension implies completely correlated character change that would violate the assumptions of the Poisson distribution that underlies the significance tests here (Lloyd et al. 2012). At the other extreme, the equally spaced option will likely unrealistically smooth out time series and potentially make it harder to reject the single-rate null (leading to Type II errors). For these reasons, the random option is recommended and is set as the default. However, because it is random this makes the function stochastic (the answer can vary each time it is run) and so the user should therefore run the function multiple times if using this option (i.e., by using a for loop) and aggregating the results at the end (e.g., as was done by previous authors; Lloyd et al. 2012; Close et al. 2015).
Secondly, the alpha value sets the significance threshold by which the likelihood ratio test's resulting p-value is compared (i.e., it is only reelevant when test_type = "lrt". Following Lloyd et al. (2012) this is set lower (0.01) than the standard 0.05 value by default as those authors found rates to be highly heterogenous in their data set (fossil lungfish). However, this should not be adopted as a "standard" value without question (just as 0.05 shouldn't). Note that the function also corrects for multiple comparisons (using the multiple_comparison_correction option) to avoid Type I errors (false positives). It does so (following Lloyd et al. 2012) using the Benjamini-Hochberg (Benjamini and Hochberg 1995) False Discovery Rate approach (see Lloyd et al. 2012 for a discussion of why), but the Bonferroni correction is also offered here (albeit not recommended).
Thirdly, polymorphisms and uncertainities create complications for assessing character changes along branches. These can occur at the tips (true polymorphisms or uncertainties in sampled taxa) and internal nodes (uncertainty over the estimated ancestral state). There are two options presented here, and applicable to both polymorphism_state and uncertainty_state (allowing these to be set separately). These are to convert such values to missing (NA) or to pick one of the possible states at random. Using missing values will increase overall uncertainty and potentially lead to Type II errors (false negatives), but represents a conservative solution. The random option is an attempt to avoid Type II errors, but can be considered unrealistic, especially if there are true polymorphisms. Additionally, the random option will again make the function stochastic meaning the user should run it multiple times amd aggregate the results. Note that if there are no polymorphisms or uncertainties in the character-taxon matrix the latter can still occur with ancestral state estimates, especially if the threshold value is set at a high value (see estimate_ancestral_states for details).
Fourthly, inapplicable characters can additionally complicate matters as they are not quite the same as missing data. I.e., they can mean that change in a particular character is not even possible along a branch. However, there is no easy way to deal with such characters at present so the user is not presented with a true option here - currently all inapplicable states are simply converted to missing values by the function. In future, though, other options may be available here. For now it is simply noted that users should be careful in making inferences if there are inapplicable characters in their data and should perhaps consider removing them with prune_cladistic_matrix to gauge their effect.
Fifthly, there are currenty two further options for assessing rates across time bins. As noted above a complication here is that character changes (the rate numerator) and character completeness (part of the rate denominator) are typically assessed on branches. However, branches will typically span time bin boundaries and hence many bins will contain only some portion of particular branches. The exact portion can be easily calculated for branch durations (the other part of the rate denominator) and the change_times option above is used to set the rate numerator, however, completeness remains complex to deal with. The first attempt to deal with this was made by Close et al. (2015) who simply used weighted mean completeness by calculating the proportion of a branch in each bin as the weight and multiplying this by each branch's completeness (the "close" option here). However, this may lead to unrealistic "smoothing" of the data and perhaps more importantly makes no account of which characters are known in a bin. Lloyd (2016) proposed an alternative "subtree" approach which assesses completeness by considering each character to be represented by a subtree where only branches that are complete are retained then branch durations in each bin are summed across subtrees such that the duration term automatically includes completeness (the "lloyd" option here). Here the latter is strongly recommended, for example, because this will lead to the same global rate across the whole tree as the branch, clade or character partitions, whereas the Close approach will not.
Sixthly, all character changes are weighted according to the weights provided in the input character-taxon matrix. In many cases these will simply all be one, although see the equalise weights option in read_nexus_matrix. However, when weights vary they can create some issues for the function. Specifically, changes are expected to be in the same (integer) units, but if weights vary then they have to be modelled accordingly. I.e., a character twice the weight of another may lead to a single change being counted as two changes. This is most problematic when the user has continuous characters which are automatically converted to gap-weighted (Thiele 1993) characters. However, this conversion creates drastically down-weighted characters and hence the user may wish to set the all_weights_integers option to TRUE. Note that reweighting will affect the results and hence shifting the weights of characters up or down will necessarily lead to shifts in the relative Type I and II errors. This is an unexplored aspect of such approaches, but is something the user should be aware of. More broadly it is recommended that continuous (or gap-weighted) characters be avoided when using this function.
Finally, the remaining options (estimate_all_nodes, estimate_tip_values, inapplicables_as_missing, polymorphism_behaviour, uncertainty_behaviour, and threshold) are all simply passed directly to estimate_ancestral_states for estimating the ancestral states and users should consult the help file for that function for further details.
Note that currently the function cannot deal with costmatrices and that the terminal versus internal option from Brusatte et al. (2014) is yet to be implemented.
Output
The output for each LRT test (i.e., the components of the branch_test_results, character_test_results, clade_test_results and time_test_results parts of the output) includes three main parts:
Rates.
p_value.
CorrectedAlpha.
Or for each AIC test there are:
Rates.
AIC.
AICc.
For each rate test the Rates part of the output is a vector of the absolute rate (number of changes per million years) for each partition in the test (in the order they were supplied to the function). So, for example, a branch rate for the sixth edge in a tree would be the rate for the sixth edge followed by the pooled rate for all other edges. The length of the vector is the length of the number of partitions.
The p_value is a single value indicating the probability that the likelihood ratio (see above and Lloyd et al. 2012) is one, i.e., that the likelihoods of the one-rate and N-rate models are the same.
The CorrectedAlpha is the alpha-value that should be used to determine the significance of the current partition test (i.e., The p_value, above). If the p_value exceeds the CorrectedAlpha then the null (single-rate) hypothesis should be accepted, if lower then the null should be rejected in favour of the N-rate hypthesis. Note that the CorrectedAlpha will not typically be the same for each partition and will also typically be different from the input alpha value due to the multiple_comparison_correction option used.
The AIC is the Akaike Information Criterion, and is relatively meaningless on its own and can only really be used to compare with the AIC values for other partitions of the data. The AICc is simply the sample-size corrected version of the AIC and is preferable when sample sizes are small.
Benjamini, Y. and Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 57, 289-300.
Brusatte, S. L., Lloyd, G. T., Wang, S. C. and Norell, M. A., 2014. Gradual assembly of avian body plan culminated in rapid rates of evolution across dinosaur-bird transition. Current Biology, 24, 2386-2392.
Close, R. A., Friedman, M., Lloyd, G. T. and Benson, R. B. J., 2015. Evidence for a mid-Jurassic adaptive radiation in mammals. Current Biology, 25, 2137-2142.
Cloutier, R., 1991. Patterns, trends, and rates of evolution within the Actinistia. Environmental Biology of Fishes, 32, 23-58.
Lloyd, G. T., 2016. Estimating morphological diversity and tempo with discrete character-taxon matrices: implementation, challenges, progress, and future directions. Biological Journal of the Linnean Society, 118, 131-151.
Lloyd, G. T., Wang, S. C. and Brusatte, S. L., 2012. Identifying heterogeneity in rates of morphological evolution: discrete character change in the evolution of lungfish (Sarcopterygii; Dipnoi). Evolution, 66, 330-348.
Ruta, M., Wagner, P. J. and Coates, M. I., 2006. Evolutionary patterns in early tetrapods. I. Rapid initial diversification followed by decrease in rates of character change. Proceedinsg of the Royal Society of London B, 273, 2107-2111.
Thiele, K.. 1993. The Holy Grail of the perfect character: the cladistic treatment of morphometric data. Cladistics, 9, 275-304.
# Set random seed:
set.seed(seed = 17)
# Generate a random tree for the Michaux data set:
time_tree <- ape::rtree(n = nrow(michaux_1989$matrix_1$matrix))
# Update taxon names to match those in the data matrix:
time_tree$tip.label <- rownames(x = michaux_1989$matrix_1$matrix)
# Set root time by making youngest taxon extant:
time_tree$root.time <- max(diag(x = ape::vcv(phy = time_tree)))
# Set four equal ength time bins spanning range of tree:
time_bins <- matrix(data = c(seq(from = time_tree$root.time, to = 0,
length.out = 5)[1:4], seq(from = time_tree$root.time, to = 0,
length.out = 5)[2:5]), ncol = 2, dimnames = list(LETTERS[1:4],
c("fad", "lad")))
# Set class as timeBins:
class(time_bins) <- "timeBins"
# Get discrete character rates:
x <- test_rates(
time_tree = time_tree,
cladistic_matrix = michaux_1989,
time_bins = time_bins,
branch_partitions = lapply(X = as.list(x = 1:nrow(time_tree$edge)), as.list),
character_partitions = lapply(X = as.list(x = 1:3), as.list),
clade_partitions = lapply(X = as.list(x = ape::Ntip(phy = time_tree) +
(2:ape::Nnode(phy = time_tree))), as.list),
time_partitions = lapply(X = as.list(x = 1:4), as.list),
change_times = "random",
alpha = 0.01,
polymorphism_state = "missing",
uncertainty_state = "missing",
inapplicable_state = "missing",
time_binning_approach = "lloyd"
)
Run the code above in your browser using DataLab