Package 'Claddis'

Measuring Morphological Diversity and Evolutionary Tempo Measures morphological diversity from discrete character data and estimates evolutionary tempo on phylogenetic trees.

Description

Measures morphological diversity from discrete character data and estimates evolutionary tempo on phylogenetic trees. Imports morphological data from #NEXUS (Maddison et al. (1997) doi:10.1093/sysbio/46.4.590) format with read_nexus_matrix(), and writes to both #NEXUS and TNT format (Goloboff et al. (2008) doi:10.1111/j.1096-0031.2008.00217.x). Main functions are test_rates(), which implements AIC and likelihood ratio tests for discrete character rates introduced across Lloyd et al. (2012) doi:10.1111/j.1558-5646.2011.01460.x, Brusatte et al. (2014) doi:10.1016/j.cub.2014.08.034, Close et al. (2015) doi:10.1016/j.cub.2015.06.047, and Lloyd (2016) doi:10.1111/bij.12746, and calculate_morphological_distances(), which implements multiple discrete character distance metrics from Gower (1971) doi:10.2307/2528823, Wills (1998) doi:10.1006/bijl.1998.0255, Lloyd (2016) doi:10.1111/bij.12746, and Hopkins and St John (2018) doi:10.1098/rspb.2018.1784. This also includes the GED correction from Lehmann et al. (2019) doi:10.1111/pala.12430. Multiple functions implement morphospace plots: plot_chronophylomorphospace() implements Sakamoto and Ruta (2012) doi:10.1371/journal.pone.0039752, plot_morphospace() implements Wills et al. (1994) doi:10.1017/S009483730001263X, plot_changes_on_tree() implements Wang and Lloyd (2016) doi:10.1098/rspb.2016.0214, and plot_morphospace_stack() implements Foote (1993) doi:10.1017/S0094837300015864. Other functions include safe_taxonomic_reduction(), which implements Wilkinson (1995) doi:10.1093/sysbio/44.4.501, map_dollo_changes() implements the Dollo stochastic character mapping of Tarver et al. (2018) doi:10.1093/gbe/evy096, and estimate_ancestral_states() implements the ancestral state options of Lloyd (2018) doi:10.1111/pala.12380. calculate_tree_length() and reconstruct_ancestral_states() implements the generalised algorithms from Swofford and Maddison (1992; no doi).

Author(s)

Graeme T. Lloyd <[email protected]>

References

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.

Examples

# Get morphological distances for Michaux (1989) data set:
distances <- calculate_morphological_distances(cladistic_matrix = michaux_1989)

# Show distances:
distances

---

Adds polymorphisms to a costmatrix

Description

Given a costmatrix and choice of method, adds polymorphic transitions to a costmatrix.

Usage

add_polymorphisms_to_costmatrix(
  costmatrix,
  polymorphism_costs = "additive",
  polymorphism_geometry = "simplex",
  polymorphism_distance = "euclidean",
  message = TRUE
)

Arguments

costmatrix	An object of class "costMatrix".
polymorphism_costs	The method to use to assign costs to the transitions to and from polymorphic states. Must be one of "additive" (the default), "geometric", "maddison", or "stratigraphic". See details.
polymorphism_geometry	If using polymorphism_costs = "geometric", the type of geometry to use. Must be one of "hypercube", "hypersphere", or "simplex" (the default). See details.
polymorphism_distance	If using polymorphism_costs = "geometric", the type of distance to use. Must be one of "euclidean" (the default), "great_circle", or "manhattan". See details.
message	Logical indicating whether (TRUE, the default) or not (FALSE) to provide messages about the process of the function.

Details

Polymorphisms - the presence of two or more discrete character states in an operational taxonomic unit - represent perhaps the most complex problem for phylogenetic inference and ancestral state estimation. This is at least in part due to the varied nature their meaning can have (Swofford and Maddison 1992). For example, they may represent true variance within a species or composite variance in some higher taxon represented by a single OTU (something to be avoided if at all possible). Consequently, they cannot - and probably should not - be treated the same by all software in all situations.

One solution to the problem of polymorphisms is to pretend they do not exist. This is the approach used by common software such as TNT (Goloboff et al. 2008; Goloboff and Catalano 2016), which treat polymorphisms as uncertainties instead. I.e., they assume only one state is the "true" value and use parsimony to (implicitly) estimate what this state is. This approach can be adopted in Claddis too by simply selecting options that treat polymorphisms as uncertainties instead. However, misrepresenting true polymorphisms as uncertainties can leads to incorrect estimates of the amount of evolution (tree length under parsimony; Nixon and Davis 1991) and can disallow, potentially inappropriately, polymorphic values as ancestral state estimates (reconstructions). Claddis therefore offers options to treat polymorphisms as "real".

There is no single agreed method on dealing with "real" polymorphisms, although many have been proposed (Wiens 1999). Due to the constraints of how Claddis works the only options implemented here - aside from simple approaches like treating them as missing data - are those that work in the costmatrix setting. These are a mixture of novel approaches, extensions of other approaches, or taken directly from Maddison and Maddison (1987).

Generally speaking the costmatrix setting makes explicit the cost of individual transitions (e.g., from the single state 0 to the polymorphic state 1 and 2) and hence the only decision to make is how these costs should be assigned. Here three approaches are offered: 1) the generalised additive approach, 2) the geometric approach (for unordered characters), and 3) the Maddison approach (for ordered characters). These are described in more detail below.

The additive approach

One way of conceiving of polymorphic tip states is to conceptually "split" a terminal branch into multiple terminals, each with a single state from the polymorphism at its' tip. For example, a branch leading to the state "0&1" would become two branches, one leading to state 0 and one to state 1. Character costs can then be calculated for each branch separately and then, in order to "merge" these into a single branch, the sum of these costs can be taken. This approach was first suggested in Hoyal Cuthill and Lloyd (in prep) and is here termed the "additive" approach (polymorphism_costs = "additive"). A major advantage to this approach is that it's rules translate across the broadest range of character types. However, there is no logical interpretation of polymorphic ancestral states. (I.e., splitting an ancestral node into two or more nodes is not possible in the same way.) Consequently, transitions "from" polymorphic states are assigned a cost of infinity, precluding them from consideration as ancestral values.

An example of the additive approach, applied to an unordered character with three states (0, 1, and 2) would thus appear as the following costmatrix:

        ---------------------------------------------
        |  0  |  1  |  2  | 0&1 | 0&2 | 1&2 | 0&1&2 |
-----------------------------------------------------
|   0   |  0  |  1  |  1  |  1  |  1  |  2  |   2   |
-----------------------------------------------------
|   1   |  1  |  0  |  1  |  1  |  2  |  1  |   2   |
-----------------------------------------------------
|   2   |  1  |  1  |  0  |  2  |  1  |  1  |   2   |
-----------------------------------------------------
|  0&1  | Inf | Inf | Inf |  0  | Inf | Inf |  Inf  |
-----------------------------------------------------
|  0&2  | Inf | Inf | Inf | Inf |  0  | Inf |  Inf  |
-----------------------------------------------------
|  1&2  | Inf | Inf | Inf | Inf | Inf |  0  |  Inf  |
-----------------------------------------------------
| 0&1&2 | Inf | Inf | Inf | Inf | Inf | Inf |   0   |
-----------------------------------------------------

The geometric approach

An alternative approach based on an idea first suggested by Maddison and Maddison (1987) does allow polymorphisms to be treated as valid ancestral states. Here the conceptual idea is that transitions involve turning "on" or "off" particular states. Thus to go from state 1 to state 0&2 would involve turning "on" states 0 and 2 turning "off" state 1, giving a total cost of three. This approach is only suitable where costs between single states are equal and symmetric - i.e., the character is unordered.

An example of this approach, applied to an unordered character with three states (0, 1, and 2) would thus appear as the following costmatrix:

        ---------------------------------------------
        |  0  |  1  |  2  | 0&1 | 0&2 | 1&2 | 0&1&2 |
-----------------------------------------------------
|   0   |  0  |  2  |  2  |  1  |  1  |  3  |   2   |
-----------------------------------------------------
|   1   |  2  |  0  |  2  |  1  |  3  |  1  |   2   |
-----------------------------------------------------
|   2   |  2  |  2  |  0  |  3  |  1  |  1  |   2   |
-----------------------------------------------------
|  0&1  |  1  |  1  |  3  |  0  |  2  |  2  |   1   |
-----------------------------------------------------
|  0&2  |  1  |  3  |  1  |  2  |  0  |  2  |   1   |
-----------------------------------------------------
|  1&2  |  3  |  1  |  1  |  2  |  2  |  0  |   1   |
-----------------------------------------------------
| 0&1&2 |  2  |  2  |  2  |  1  |  1  |  1  |   0   |
-----------------------------------------------------

Note: to maintain a base cost of one the original single-state-to-single-state transitions are doubled from one to two. In order to avoid this affecting tree lengths the character weight is halved.

This Maddison and Maddison (1987) approach is extended here by drawing on an analogy, namely that this cost assignment is identical to considering a hypercube whose vertices lie at one of two (presence or absence) positions on orthogonal axes representing individual states. Thus the polymorphism 0&2 is at the coordinates (0 = present, 1 = absent, 2 = present). Here costs can be assigned by taking the minimum Manhattan distance between each pair of vertices. NB: here the origin (absence of all states) is unoccupied, but this does not affect the assignment of costs.

This analogy is why the approach is here termed the "geometric" approach (polymorphism_costs = "geometric"). Importantly, although the example above only represents the regular cube (three states = three dimensions) the approach translates to any number of states (i.e., any number of dimensions) without loss of generality, and hence this particular option is properly termed a hypercube (i.e., polymorphism_geometry = "hypercube" and polymorphism_distance = "manhattan").

In of itself this does not change what was proposed by Maddison and Maddison. However, the geometric setting allows us to consider both different shapes (topologies) and different distance measures, as long as the same generality to higher dimensions holds. Two additional shapes are offered here. These are the hypersphere (the N-dimensional circle) and the simplex (the N-dimensional equilateral triangle). In addition, two distances are also offered. These are Euclidean (the straight line distance between two points) and Great circle (the shortest distance across the "surface" of a hypersphere) distances.

These distances are intended to match respective straights (hypercube and Manhattan, hypersphere and Great Circle, simplex and Euclidean), but any combination is permitted. Furthermore, and as above, all costs are rescaled such that the base transition cost is always one with the character weight modified accordingly.

The Maddison approach

Although the above approach is based on Maddison and Maddison (1987) the Maddison name (polymorphism_costs = "maddison") is here reserved for an approach proposed by the same authors (Maddison and Maddison 2000) for use with ordered characters. The concept of this approach shares the idea of "switching" states on and off but here adds the idea of adjacency of switches. In other words, intermediate states (the essential element of an ordered character) are additional switches that must be turned "on" and "off" again in order to "reach" other states. Thus to go from state 0 to state 0 and 2 for the linear ordered character 0-1-2 we would have to turn "on" state to reach state "2", which is also turned on. Then we would turn off state "1", meaning total of three switches were flipped which would be the cost. Note: that here we are not turning "off" state 0 as it is present at both ends of the transition.

Extending our three state linear ordered example to all possible transitions we get the costmatrix:

        ---------------------------------------------
        |  0  |  1  |  2  | 0&1 | 0&2 | 1&2 | 0&1&2 |
-----------------------------------------------------
|   0   |  0  |  2  |  4  |  1  |  3  |  3  |   2   |
-----------------------------------------------------
|   1   |  2  |  0  |  2  |  1  |  3  |  1  |   2   |
-----------------------------------------------------
|   2   |  4  |  2  |  0  |  3  |  3  |  1  |   2   |
-----------------------------------------------------
|  0&1  |  1  |  1  |  3  |  0  |  2  |  2  |   1   |
-----------------------------------------------------
|  0&2  |  3  |  3  |  3  |  2  |  0  |  2  |   1   |
-----------------------------------------------------
|  1&2  |  3  |  1  |  1  |  2  |  2  |  0  |   1   |
-----------------------------------------------------
| 0&1&2 |  2  |  2  |  2  |  1  |  1  |  1  |   0   |
-----------------------------------------------------

NB: Again the single-state-to-single-state cost is modified (doubled) meaning the character weight is halved.

The stratigraphic approach

The stratigraphic case seems like it ought to follow the same rule as irreversible characters (see below), but here "polymorphisms" have a different logical meaning. Specifically they encode the presence of a taxon in multiple time units. Thus, transitions to a stratigraphic polymorphism are based on the oldest state. By contrast, polymorphisms can never be ancestral states as there would be no basis to choose between (e.g.) state 0 and state 0 and 1.

The option polymorphism_costs = "stratigraphic" is only available where costmatrix$type is "stratigraphic" and, where costmatrix$type is "stratigraphic", polymorphism_costs = "stratigraphic" is the only option.

Additional considerations

Aside from these three options there are some further considerations the user should make with respect to some more complex character types. These are discussed below under appropriate headers.

Dollo and irreversible polymorphisms

At first blush Dollo and irreversible (Camin-Sokal) characters appear to confound costmatrix construction. More specifically, for a Dollo character, the state 0 and 1 would suggest state 1 was acquired previously, but lost in some member(s) of the OTU. For an irreversible character the opposite would be true, and we would assume that state 1 was acquired independently by some member(s) of the OTU and state 0 was acquired previously. It is not possible to code either of these outcomes under the geometric or Maddison approaches. However, the additive approach perfectly captures this nuance. Thus tree length calculation and ancestral state estimation for Dollo and Camin-Sokal characters is possible under the additive approach and only the possibility of polymorphic ancestral states is excluded.

Transitions between uncertainties and polymorphisms

Because Claddis also permits users to include uncertainties in costmatrices transitions between these and polymorphisms must also be considered. Here this is automated by using the "minimum rule". In practice this means a transition from, for example, state 0 and 1 to state 1 or 2 will take the lowest cost of the two options (from state 0 and 1 to state 1 or from state 0 and 1 to state 2). In the case of the additive approach these costs will always be infinite. In all cases transitions from uncertainties are assigned infinite cost as these are not considered valid ancestral states.

Polymorphism limits and gap-weighted characters

Even where polymorphisms may be appropriate (i.e., for some ordered and unordered characters) they still represent a major increase to costmatrix size that can cause memory limits to be hit. Specifically, and as shown in the hypercube example above, there will be as many possible states as N^2 - 1, where N is the number of single (i.e., non-polymorphic) states and the minus one indicates the exclusion of the unrepresentable origin value. Thus the size of costmatrix required grows very quickly:

------------------------------
| N states | Costmatrix size |
------------------------------
|     2    |      3 x 3      |
|     3    |      7 x 7      |
|     4    |     15 x 15     |
|     5    |     31 x 31     |
|     6    |     63 x 63     |
|     7    |    127 x 127    |
|     8    |    255 x 255    |
|     9    |    511 x 511    |
|    10    |   1023 x 1023   |
|    11    |   2047 x 2047   |
|    12    |   4095 x 4095   |
|    13    |   8191 x 8191   |
|    14    |  16383 x 16383  |
------------------------------

Because of this, the function will become extremely slow for these higher values and here is hard-capped at no more than fourteen states. Consequently any gap-weighted characters are also most likely inappropriate for polymorphism use (as well as being too memory intensive).

Value

An object of class "costMatrix".

Author(s)

Graeme T. Lloyd [email protected]

References

Goloboff, P. A. and Catalano, S. A., 2016. TNT version 1.5, including a full implementation of phylogenetic morphometrics/ Cladistics, 32. 221-238

Goloboff, P., Farris, J. and Nixon, K., 2008. TNT, a free program for phylogenetic analysis. Cladistics, 24, 774-786.

Hoyal Cuthill and Lloyd, in prep.

Maddison, W. P. and Maddison, D. R., 1987. MacClade 2.1, computer program and manual. Cambridge, Massachusetts.

Nixon, K. C. and Davis, J. I., 1991. Polymorphic taxa, missing values and cladistic analysis. Cladistics, 7, 233-241.

Swofford, D. L. and Maddison, W. P., 1992. Parsimony, character-state reconstructions, and evolutionary inferences. In R. L. Mayden (ed.), Systematics, Historical Ecology, and North American Freshwater Fishes. Stanford University Press, Stanford. pp187-223.

Wiens, J. J., 1999. Polymorphism in systematics and comparative biology. Annual Review of Ecology and Systematics, 30, 327-362.

Examples

# Generate an example three-state unordered character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "unordered"
)

# Generate an example three-state ordered character costmatrix:
ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "ordered"
)

# Generate an example three-state ordered character costmatrix with uncertainties already included:
ordered_uncertainty_costmatrix <- list(
  size = 7,
  n_states = 3,
  single_states = c("0", "1", "2"),
  type = "ordered",
  costmatrix = matrix(
    data = c(
      0, 1, 2, 0, 0, 1, 0,
      1, 0, 1, 0, 1, 0, 0,
      2, 1, 0, 1, 0, 0, 0,
      Inf, Inf, Inf, 0, Inf, Inf, Inf,
      Inf, Inf, Inf, Inf, 0, Inf, Inf,
      Inf, Inf, Inf, Inf, Inf, 0, Inf,
      Inf, Inf, Inf, Inf, Inf, Inf, 0
    ),
    nrow = 7,
    byrow = TRUE,
    dimnames = list(
      c(as.character(x = 0:2), "0/1", "0/2", "1/2", "0/1/2"),
      c(as.character(x = 0:2), "0/1", "0/2", "1/2", "0/1/2")
    )
  ),
  symmetry = "Symmetric",
  includes_polymorphisms = FALSE,
  polymorphism_costs = "additive",
  polymorphism_geometry = "simplex",
  polymorphism_distance = "euclidean",
  includes_uncertainties = TRUE,
  pruned = FALSE,
  dollo_penalty = 999,
  base_age = 1,
  weight = 1
)
class(ordered_uncertainty_costmatrix) <- "costMatrix"

# Generate an example five-state stratigraphic character costmatrix:
stratigraphic_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 4,
  character_type = "stratigraphy",
  state_ages = c(0, 1.3, 5.3, 8.0, 11.1)
)

# Add polymorphisms to unordered costmatrix using additive method:
unordered_costmatrix_additive_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = unordered_costmatrix,
  polymorphism_costs = "additive"
)

# Show unordered costmatrix using additive method:
unordered_costmatrix_additive_polymorphisms$costmatrix

# Add polymorphisms to unordered costmatrix using geometric simplex method:
unordered_costmatrix_simplex_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = unordered_costmatrix,
  polymorphism_costs = "geometric",
  polymorphism_geometry = "simplex",
  polymorphism_distance = "euclidean"
)

# Show unordered costmatrix using geometric simplex method:
unordered_costmatrix_simplex_polymorphisms$costmatrix

# Add polymorphisms to unordered costmatrix using geometric hypercube method:
unordered_costmatrix_hypercube_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = unordered_costmatrix,
  polymorphism_costs = "geometric",
  polymorphism_geometry = "hypercube",
  polymorphism_distance = "manhattan"
)

# Show unordered costmatrix using geometric hypercube method:
unordered_costmatrix_hypercube_polymorphisms$costmatrix

# Add polymorphisms to ordered costmatrix using additive method:
ordered_costmatrix_additive_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = ordered_costmatrix,
  polymorphism_costs = "additive"
)

# Show ordered costmatrix using additive method:
ordered_costmatrix_additive_polymorphisms$costmatrix

# Add polymorphisms to ordered costmatrix using maddison method:
ordered_costmatrix_maddison_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = ordered_costmatrix,
  polymorphism_costs = "maddison"
)

# Show ordered costmatrix using Maddison method:
ordered_costmatrix_maddison_polymorphisms$costmatrix

# Add polymorphisms to ordered uncertainty costmatrix using additive method:
ordered_uncertainty_costmatrix_additive_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = ordered_uncertainty_costmatrix,
  polymorphism_costs = "additive"
)

# Show costmatrix, with polymorphism-to-uncertainty transitions as infinities:
ordered_uncertainty_costmatrix_additive_polymorphisms$costmatrix

# Add polymorphisms to ordered uncertainty costmatrix using Maddison method:
ordered_uncertainty_costmatrix_maddison_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = ordered_uncertainty_costmatrix,
  polymorphism_costs = "maddison"
)

# Show costmatrix, with polymorphism-to-uncertainty transitions
# interpolated using minimum cost rule:
ordered_uncertainty_costmatrix_maddison_polymorphisms$costmatrix

# Add polymorphisms to stratigraphic costmatrix using stratigraphic method:
stratigraphic_costmatrix_polymorphisms <- add_polymorphisms_to_costmatrix(
  costmatrix = stratigraphic_costmatrix,
  polymorphism_costs = "stratigraphic"
)

# Show stratigraphic costmatrix using stratigraphic method:
stratigraphic_costmatrix_polymorphisms$costmatrix

---

Adds uncertainties to a costmatrix

Description

Adds uncertainty transitions to a costmatrix.

Usage

add_uncertainties_to_costmatrix(costmatrix, message = TRUE)

Arguments

costmatrix	An object of class "costMatrix".
message	Logical indicating whether (TRUE, the default) or not (FALSE) to provide messages about the process of the function.

Details

Under a parsimony algorithm uncertainties indicate that the observed state cannot be restricted to a single value but some finite set of values. In #NEXUS format these are indicated with curly braces, , and within Claddis by a slash. For example, state 0 or 1 would be 0/1.

For most operations (calculating the length of a tree, estimating ancestral states) there is no need to formally add uncertainties to a costmatrix. Instead, the uncertainty can be established via setting up the tip states in the Swofford and Maddison (1992) approach. However, there are some usages within Claddis where it is useful to include them in the costmatrix itself, for example calculating minimum step counts for homoplasy metrics.

Here, uncertainties are added to a costmatrix in a similar way to polymorphisms. I.e., all possible uncertainty combinations are calculated and then the associated transition cost calculated and inserted. Note that here only one approach for calculating this cost is applied - the minimu rule of Hoyal Cuthill and Lloyd (in prep.).

Importantly costs of transitioning from uncertainties are assigned infinite cost. There aree two main reasons for this. Firstly, under parsimony and the minimum rule allowing uncertainties as ancestral states would lead to a tree length of zero and an ancestral state for every internal node that includes every possible state being favoured, regardless of the states at the tips. In other words, it would have no analytical value. Secondly, assigning uncertainty of ancestral states is best done by generating all (single or polymorphic) most parsimonious reconstructions and summarising these.

Note that in many practical cases an uncertainty involving all possible states - e.g., the uncertianty 0/1/2 for a three-state character - operates the same as a missing value (NA). However, it may still be useful in some circumstances. For example, to denote that the coding is definitively applicable (as opposed to inapplicable) or to properly record uncertainty rather than missing values for summary statistical purposes.

Finally, it should be noted that uncertainties are the only allowable breaking of the costmatrix rule that all off-diagonal transition costs must be positive. I.e., costs of zero from a single or polymorphic state to an uncertainty state are permitted.

Value

An object of class "costMatrix".

Author(s)

Graeme T. Lloyd [email protected]

References

Hoyal Cuthill and Lloyd, in prep.

Swofford, D. L. and Maddison, W. P., 1992. Parsimony, character-state reconstructions, and evolutionary inferences. In R. L. Mayden (ed.), Systematics, Historical Ecology, and North American Freshwater Fishes. Stanford University Press, Stanford. pp187-223.

Examples

# Generate an example three-state unordered character costmatrix:
# Generate an example three-state unordered character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "unordered"
)

# Generate an example three-state ordered character costmatrix:
ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "ordered"
)

# Generate an example three-state ordered character costmatrix with polymorphism included:
unordered_polymorphism_costmatrix <- list(
  size = 7,
  n_states = 3,
  single_states = c("0", "1", "2"),
  type = "unordered",
  costmatrix = matrix(
    data = c(
      0, 2, 2, 1, 1, 3, 2,
      2, 0, 2, 1, 3, 1, 2,
      2, 2, 0, 3, 1, 1, 2,
      1, 1, 3, 0, 2, 2, 1,
      1, 3, 1, 2, 0, 2, 1,
      3, 1, 1, 2, 2, 0, 1,
      2, 2, 2, 1, 1, 1, 0
    ),
    nrow = 7,
    byrow = TRUE,
    dimnames = list(
      c(as.character(x = 0:2), "0&1", "0&2", "1&2", "0&1&2"),
      c(as.character(x = 0:2), "0&1", "0&2", "1&2", "0&1&2")
    )
  ),
  symmetry = "Symmetric",
  includes_polymorphisms = TRUE,
  polymorphism_costs = "geometric",
  polymorphism_geometry = "hypercube",
  polymorphism_distance = "manhattan",
  includes_uncertainties = FALSE,
  pruned = FALSE,
  dollo_penalty = 999,
  base_age = 1,
  weight = 0.5
)
class(unordered_polymorphism_costmatrix) <- "costMatrix"

# Add uncertainties to unordered costmatrix:
unordered_costmatrix_uncertainties <- add_uncertainties_to_costmatrix(
  costmatrix = unordered_costmatrix
)

# Show unordered costmatrix with uncertainties included:
unordered_costmatrix_uncertainties$costmatrix

# Add uncertainties to ordered costmatrix:
ordered_costmatrix_uncertainties <- add_uncertainties_to_costmatrix(
  costmatrix = ordered_costmatrix
)

# Show ordered costmatrix with uncertainties included:
ordered_costmatrix_uncertainties$costmatrix

# Add uncertainties to unordered costmatrix with polymorphisms:
unordered_polymorphism_costmatrix_uncertainties <- add_uncertainties_to_costmatrix(
  costmatrix = unordered_polymorphism_costmatrix
)

# Show unordered polymorphism costmatrix with uncertainties included:
unordered_polymorphism_costmatrix_uncertainties$costmatrix

---

Aligns a phylogenetic matrix block

Description

Given a block of taxa and characters aligns text so each character block begins at same point.

Usage

align_matrix_block(matrix_block)

Arguments

matrix_block

The matrix block as raw input text.

Details

The function serves to help build NEXUS files by neatly aligning raw text blocks of taxa and characters. Or in simple terms it takes input that looks like this:

Allosaurus  012100?1011
Abelisaurus  0100???0000
Tyrannosaurus  01012012010
Yi  10101?0????

And turns it into something that looks like this:

Allosaurus     012100?1011
Abelisaurus    0100???0000
Tyrannosaurus  01012012010
Yi             10101?0????

I use this in building the NEXUS files on my site, graemetlloyd.com.

Value

Nothing is returned, instead the aligned block is sent to the clipboard ready for pasting into a text editor.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Build example block from above:
x <- paste(c(
  "Allosaurus  012100?1011",
  "Abelisaurus  0100???0000",
  "Tyrannosaurus  01012012010",
  "Yi  10101?0????"
), collapse = "\n")

# Look at block pre-alignment:
x

# Align block and place on clipboard:
## Not run: 
align_matrix_block(x)

## End(Not run)

# To test the response open a text editor and paste the
# contents of the clipboard.

---

Assign taxa to time bins

Description

Given a set of first and last appearances assigns a set of taxa to a series of time bins.

Usage

assign_taxa_to_bins(taxon_ages, time_bins)

Arguments

taxon_ages	A matrix of taxon ages, with columns for first ("fad") and last ("lad") appearances and rownames correspodning to taxon names.
time_bins	An object of class timeBins.

Details

The various disparity plotting functions (plot_chronophylomorphospace, plot_morphospace_stack, plot_morphospace, plot_multi_morphospace) are designed to allow assignment of taxa to named groups so that these groups may be assigned different colours when plotting. One way taxa may be grouped is temporally, by assignment to a series of time bins.

There are many ways this may be automated and this function provides a very simple one: if the first and last appearance dates of a taxon overlap with a time bin then it can be assigned to that time bin. (In practice, taxa often have multiple occurrences with "ranges" that really represent uncertainty around their true age.)

Note that it is recommended that time bins be named without special characters beyond letters and underscores.

Value

An object of class taxonGroups.

Author(s)

Graeme T. Lloyd [email protected]

See Also

plot_chronophylomorphospace, plot_morphospace_stack, plot_morphospace, plot_multi_morphospace, ordinate_cladistic_matrix

Examples

# Build example time bins:
time_bins <- matrix(data = c(443.8, 358.9, 358.9, 298.9, 298.9, 251.9,
  251.9, 201.3, 201.3, 145.0, 145.0, 65.5, 65.5, 23.03), ncol = 2,
  byrow = TRUE, dimnames = list(c("Silurodevonian", "Carboniferous",
  "Permian", "Triassic", "Jurassic", "Cretaceous", "Paleogene"),
  c("fad", "lad")))

# Set class as timeBins:
class(time_bins) <- "timeBins"

# Build example taxon ages:
taxon_ages <- matrix(data = c(385.3, 374.5, 407, 374.5, 251, 228, 385.3,
  251, 251, 251, 391.8, 251, 251, 228, 385.3, 391.8, 391.8, 385.3, 311.7,
  359.2, 359.2, 416, 407, 407, 407, 407, 385.3, 397.5, 385.3, 161.2, 385.3,
  345.3, 318.1, 385.3, 228, 385.3, 385.3, 385.3, 385.3, 385.3, 385.3, 385.3,
  385.3, 385.3, 391.8, 407, 391.8, 374.5, 407, 70.6, 311.7, 407, 145.5, 251,
  65.5, 251, 112, 374.5, 374.5, 374.5, 385.3, 311.7, 249.7, 359.2, 391.8,
  374.5, 385.3, 83.5, 418.7, 251, 385.3, 391.8, 374.5, 345.3, 385.3, 385.3,
  407, 411.2, 397.5, 345.3, 374.5, 407, 216.5, 326.4, 411.2, 411.2, 374.5,
  359.2, 391.8, 359.2, 245, 216.5, 374.5, 245, 245, 245, 385.3, 245, 245,
  199.6, 374.5, 385.3, 385.3, 374.5, 306.5, 345.3, 345.3, 411.2, 397.5,
  397.5, 397.5, 397.5, 374.5, 391.8, 374.5, 145.5, 374.5, 326.4, 311.7,
  374.5, 199.6, 374.5, 374.5, 374.5, 374.5, 374.5, 374.5, 374.5, 374.5,
  374.5, 385.3, 397.5, 385.3, 359.2, 397.5, 65.5, 306.5, 397.5, 99.6, 245,
  23.03, 245, 99.6, 359.2, 359.2, 359.2, 374.5, 306.5, 247.4, 318.1, 385.3,
  359.2, 374.5, 70.6, 416, 250.4, 374.5, 385.3, 359.2, 326.4, 374.5, 374.5,
  397.5, 407, 391.8, 326.4, 359.2, 397.5, 203.6, 318.1, 407, 407),
  ncol = 2, dimnames = list(c("Adololopas_moyasmithae", "Adelargo_schultzei",
  "Amadeodipterus_kencampbelli", "Andreyevichthys_epitomus",
  "Aphelodus_anapes", "Archaeoceratodus_avus", "Archaeonectes_pertusus",
  "Arganodus_atlantis", "Ariguna_formosa", "Asiatoceratodus_sharovi",
  "Barwickia_downunda", "Beltanodus_ambilobensis", "Ceratodus_formosa",
  "Ceratodus_latissimus", "Chirodipterus_australis",
  "Chirodipterus_onawwayensis", "Chirodipterus_rhenanus",
  "Chirodipterus_wildungensis", "Conchopoma_gadiforme", "Ctenodus_romeri",
  "Delatitia_breviceps", "Diabolepis_speratus", "Dipnorhynch_cathlesae",
  "Dipnorhynchus_sussmilchi", "Dipnorhynchus_kiandrensis",
  "Dipnorhynchus_kurikae", "Dipterus_cf_valenciennesi",
  "Dipterus_valenciennesi", "Eoctenodus_microsoma",
  "Ferganoceratodus_jurassicus", "Fleurantia_denticulata",
  "Ganopristodus_splendens", "Gnathorhiza_serrata", "Gogodipterus_paddyensis",
  "Gosfordia_truncata", "Griphognathus_minutidens", "Griphognathus_sculpta",
  "Griphognathus_whitei", "Grossipterus_crassus", "Holodipterus_elderae",
  "Holodipterus_gogoensis", "Robinsondipterus_longi",
  "Asthenorhynchus_meemannae", "Holodipterus_santacrucensis",
  "Howidipterus_donnae", "Ichnomylax_kurnai", "Iowadipterus_halli",
  "Jarvikia_arctica", "Jessenia_concentrica", "Lepidosiren_paradoxa",
  "Megapleuron_zangerli", "Melanognathus_canadensis",
  "Metaceratodus_wollastoni", "Microceratodus_angolensis",
  "Mioceratodus_gregoryi", "Namatozodia_pitikanta", "Neoceratodus_forsteri",
  "Nielsenia_nordica", "Oervigia_nordica", "Orlovichthys_limnatis",
  "Palaeodaphus_insignis", "Palaeophichthys_parvulus", "Paraceratodus_germaini",
  "Parasagenodus_sibiricus", "Pentlandia_macroptera",
  "Phaneropleuron_andersoni", "Pillararhynchus_longi", "Protopterus_annectens",
  "Psarolepis_romeri", "Ptychoceratodus_serratus", "Rhinodipterus_secans",
  "Rhinodipterus_ulrichi", "Rhynchodipterus_elginensis", "Sagenodus_inaequalis",
  "Scaumenacia_curta", "Soederberghia_groenlandica",
  "Sorbitorhynchus_deleaskitus", "Speonesydrion_iani", "Stomiahykus_thlaodus",
  "Straitonia_waterstoni", "Sunwapta_grandiceps", "Tarachomylax_oepiki",
  "Tellerodus_sturi", "Tranodis_castrensis", "Uranolophus_wyomingensis",
  "Westollrhynchus_lehmanni"), c("fad", "lad")))

# Assign taxa to time bins:
assign_taxa_to_bins(taxon_ages = taxon_ages, time_bins = time_bins)

---

Counts the changes in a series of time bins

Description

Given a vector of dates for a series of time bins and another for the times when a character change occurred will return the total number of changes in each bin.

Usage

bin_changes(change_times, time_bins)

Arguments

change_times	A vector of ages in millions of years at which character changes are hypothesised to have occurred.
time_bins	An object of class timeBins.

Details

Calculates the total number of evolutionary changes in a series of time bins. This is intended as an internal function for rate calculations, but could be used for other purposes (e.g., counting any point events in a series of time bins).

Value

A vector giving the number of changes for each time bin. Names indicate the maximum and minimum (bottom and top) values for each time bin.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a random dataset of 100 changes (between 100 and 0 Ma):
change_times <- stats::runif(n = 100, min = 0, max = 100)

# Create 10 equal-length time bins:
time_bins <- matrix(data = c(seq(from = 100, to = 10, length.out = 10),
  seq(from = 90, to = 0, length.out = 10)), ncol = 2,
  dimnames = list(LETTERS[1:10], c("fad", "lad")))

# Set class as timeBins:
class(time_bins) <- "timeBins"

# Get N changes for each bin:
bin_changes(change_times = change_times, time_bins = time_bins)

---

Phylogenetic character completeness in time-bins

Description

Given a cladistic matrix, time-scaled tree, and set of time bin boundaries will return the proportional character completeness in each bin.

Usage

bin_character_completeness(
  cladistic_matrix,
  time_tree,
  time_bins,
  plot = FALSE,
  confidence.interval = 0.95
)

Arguments

cladistic_matrix	A cladistic matrix in the form imported by read_nexus_matrix.
time_tree	A time-scaled phylogenetic tree containing all the taxa in cladistic_matrix.
time_bins	An object of class timeBins.
plot	An optional choice to plot the results (default is FALSE).
confidence.interval	The confidence interval to be used as a proportion (0 to 1). Default is 0.95 (i.e., 95%).

Details

Character completeness metrics have been used as an additional metric for comparing fossil record quality across time, space, and taxa. However, these only usually refer to point samples of fossils in bins, and not our ability to infer information along the branches of a phylogenetic tree.

This function returns the proportional phylogenetic character completeness for a set of time bins.

Value

A list summarising the mean, upper and lower confidence interval, and per character proportional character completeness in each time bin.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a random tree for the Day et al. 2016 data set:
day_2016tree <- ape::rtree(n = nrow(day_2016$matrix_1$matrix))
day_2016tree$tip.label <- rownames(x = day_2016$matrix_1$matrix)
day_2016tree$root.time <- max(diag(x = ape::vcv(phy = day_2016tree)))

# Build ten equal-length time bins spanning the tree:
time_bins <- matrix(data = c(seq(from = day_2016tree$root.time,
  to = day_2016tree$root.time - max(diag(x = ape::vcv(phy = day_2016tree))),
  length.out = 11)[1:10], seq(from = day_2016tree$root.time,
  to = day_2016tree$root.time - max(diag(x = ape::vcv(phy = day_2016tree))),
  length.out = 11)[2:11]), ncol = 2, dimnames = list(LETTERS[1:10], c("fad", "lad")))

# Set class as timeBins:
class(time_bins) <- "timeBins"

# Get proportional phylogenetic character completeness in ten equal-length
# time bins:
bin_character_completeness(
  cladistic_matrix = day_2016,
  time_tree = day_2016tree,
  time_bins = time_bins
)

# Same, but with a plot:
bin_character_completeness(
  cladistic_matrix = day_2016,
  time_tree = day_2016tree,
  time_bins = time_bins,
  plot = TRUE
)

---

Edge-lengths present in time-bins

Description

Given a time-scaled tree and set of time bin boundaries will sum the edge-lengths present in each bin.

Usage

bin_edge_lengths(time_tree, time_bins, pruned_tree = NULL)

Arguments

time_tree	A time-scaled tree in phylo format with a $root.time value.
time_bins	An object of class timeBins.
pruned_tree	A time-scaled tree in phylo format with a $root.time value that is a subset of time_tree.

Details

Calculates the total edge duration of a time-scaled tree present in a series of time bins. This is intended as an internal function for rate calculations, but may be of use to someone.

The option of using a pruned_tree allows the user to correctly classify internal and terminal branches in a subtree of the larger tree. So for example, if taxa A and B are sisters then after pruning B the subtree branch leading to A is composed of an internal and a terminal branch on the complete tree.

Value

binned_edge_lengths	A vector giving the summed values in millions of years for each time bin. Names indicate the maximum and minimum values for each time bin.
binned_terminal_edge_lengths	As above, but counting terminal edges only.
binned_internal_edge_lengths	As above, but counting internal edges only.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a random 10-taxon tree:
time_tree <- ape::rtree(n = 10)

# Add root age:
time_tree$root.time <- max(diag(ape::vcv(time_tree)))

# Create time bins:
time_bins <- matrix(data = c(seq(from = time_tree$root.time, to = 0,
  length.out = 11)[1:10], seq(from = time_tree$root.time, to = 0,
  length.out = 11)[2:11]), ncol = 2, dimnames = list(LETTERS[1:10],
  c("fad", "lad")))

# Set class:
class(time_bins) <- "timeBins"

# Get edge lengths for each bin:
bin_edge_lengths(time_tree = time_tree, time_bins = time_bins)

---

Creates a morphological data file from a matrix

Description

Creates a morphological data file from a character-taxon matrix.

Usage

build_cladistic_matrix(
  character_taxon_matrix,
  header = "",
  character_weights = NULL,
  ordering = NULL,
  symbols = NULL,
  equalise.weights = FALSE,
  ignore_duplicate_taxa = FALSE
)

Arguments

character_taxon_matrix	A Character-Taxon (columns-rows) matrix, with taxon names as rownames.
header	A scalar indicating any header text (defaults to an empty string: "").
character_weights	A vector specifying the weights used (if not specified defaults to 1).
ordering	A vector indicating whether characters are ordered ("ordered") or unordered ("unordered") (if no specified defaults to ordered).
symbols	The symbols to use if writing to a file (defaults to the numbers 0:9 then the letters A to V).
equalise.weights	Optional that overrides the weights specified above make all characters truly equally weighted.
ignore_duplicate_taxa	Logical indicating whether or not to ignore (allow; TRUE) duplicate taxa or not (FALSE; default).

Details

Claddis generally assumes that matrices will be imported into R from the #NEXUS format, but in some cases (e.g., when using simulated data) it might be desirable to build a matrix within R. This function allows the user to convert such a matrix into the format required by other Claddis functions as long as it only contains a single block.

NB: Currently the function cannot deal directly with costmatrices or continuous characters.

Value

topper

Contains any header text or costmatrices and pertains to the entire file.

matrix_N

One or more matrix blocks (numbered 1 to N) with associated information pertaining only to that matrix block. This includes the block name (if specificed, NA if not), the block datatype (one of "CONTINUOUS", "DNA", "NUCLEOTIDE", "PROTEIN", "RESTRICTION", "RNA", or "STANDARD"), the actual matrix (taxa as rows, names stored as rownames and characters as columns), the ordering type of each character ("ordered", "unordered" etc.), the character weights, the minimum and maximum values (used by Claddis' distance functions), and the original characters (symbols, missing, and gap values) used for writing out the data.

Author(s)

Graeme T. Lloyd [email protected]

See Also

compactify_cladistic_matrix, prune_cladistic_matrix, read_nexus_matrix, safe_taxonomic_reduction, write_nexus_matrix, write_tnt_matrix

Examples

# Create random 10-by-50 matrix:
character_taxon_matrix <- matrix(sample(c("0", "1", "0&1", NA, ""),
  500,
  replace = TRUE
),
nrow = 10, dimnames =
  list(apply(matrix(sample(LETTERS, 40,
    replace = TRUE
  ), nrow = 10), 1, paste,
  collapse = ""
  ), c())
)

# Reformat for use elsewhere in Claddis:
build_cladistic_matrix(character_taxon_matrix)

---

Calculate the maximum tree length, g, under parsimony

Description

Given a costmatrix and set of tip states returns the longest possible tree length under maximum parsimony.

Usage

calculate_g(
  costmatrix,
  tip_states,
  polymorphism_behaviour = "polymorphism",
  uncertainty_behaviour = "uncertainty"
)

Arguments

costmatrix	An object of class costMatrix.
tip_states	A character vector of tip states, with polymorphic states separated by &, uncertainties by /, missing values as NA, and inapplicables as empty strings "".
polymorphism_behaviour	One of either "missing", "uncertainty", "polymorphism", or "random". See details.
uncertainty_behaviour	One of either "missing", "uncertainty", "polymorphism", or "random". See details.

Details

The maximum cost a character could have on any tree under maximum parsimony, termed g, depends on both the individual state-to-state transition costs (captured by a costmatrix) and the sampled states (i.e., the tip_states input). In practice this is the maximum parsimony length on the star tree. This length cannot be exceeded by any other tree (Hoyal Cuthill and Lloyd, in review). Note: this is standard practice in phylogenetics software and is also how both PAUP* (Swofford 2003) and TNT (Goloboff et al. 2008; Goloboff and Catalano 2016) calculate maximum cost.

Special cases

A number of special cases apply to calculating g and are discussed further below.

Polymorphisms

Polymorphisms remain a complex problem in phylogenetics and here multiple options are provided to deal with them. These include: 1. "missing" - where they are simply replaced by a missing value (see below), 2. "uncertainty" - where they are treated as uncertainties instead (see below), 3. "polymorphism" - where they are treated as genuinely polymorphic, and 4. "random" - where one of the tip states is selected at random.

Options 1, 2, and 4 can be seen as undercounting the true amount of evolution that has occurred. However, how to correctly count this amount is unclear. If option 3 is chosen then polymorphic states must be present in costmatrix and users should refer to the add_polymorphisms_to_costmatrix function for details on available options.

Uncertainties

Uncertainties are much simpler to deal with than polymorphisms and a means to incorporate them into length counts was laid out in Swofford and Maddison (1992). Indeed, popular software such as PAUP* (Swofford 2003) and TNT (Goloboff et al. 2008; Goloboff and Catalano 2016) simply treat polymorphisms as uncertainties perhaps because of this. There is still a concern of undercounting evolutionary change for uncertainties in the maximum parsimony context as the cheapest possible state will in effect be used everytime, whereas future study that removes uncertainty may reveal s higher cost state to be the true value. As such the same options are offered for uncertainties as polymorphisms, including to treat them as polymorphisms although this should probably only be done where they were miscoded in the first place.

Again, if using uncertainties as uncertainties, these must be included in the costmatrix and this can be done by using the add_uncertainties_to_costmatrix function.

Missing values

In practice missing values (NA) may exist amongst tip_states. These are permitted, and in practice are mathematically and practically equivalent to a statement that a tip could be any state present in the costmatrix (i.e., a special case of an uncertainty where no state can be ruled out). However, it should be considered in interpretation that g will typically become smaller as the number of missing values increases.

Inapplicable values

Inapplicable values ("") may also exist amongst tip_states. These are conceptually different to missing values as there is no possibility that they can ever be (re)coded. Currently these are treated exactly the same as missing values, but again the user should apply caution in interpreting g in such cases as again it will be smaller than otherwise identical characters with fewer inapplicable values.

Value

A single value indicating the maximum length, g. Note: this is not modified by costmatrix$weight.

Author(s)

Graeme T. Lloyd [email protected] and Jen Hoyal Cuthill [email protected]

References

Goloboff, P. A. and Catalano, S. A., 2016. TNT version 1.5, including a full implementation of phylogenetic morphometrics/ Cladistics, 32. 221-238

Goloboff, P., Farris, J. and Nixon, K., 2008. TNT, a free program for phylogenetic analysis. Cladistics, 24, 774-786.

Hoyal Cuthill, J. F. and Lloyd, G. T., in press. Measuring homoplasy I: comprehensive measures of maximum and minimum cost under parsimony across discrete cost matrix character types. Cladistics, bold, .

Swofford, D. L., 2003. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts.

Swofford, D. L. and Maddison, W. P., 1992. Parsimony, character-state reconstructions, and evolutionary inferences. In R. L. Mayden (ed.), Systematics, Historical Ecology, and North American Freshwater Fishes. Stanford University Press, Stanford. pp187-223.

See Also

calculate_gmax

Examples

# Create a Type I character costmatrix:
constant_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 0,
  character_type = "unordered"
)

# Calculate g for the case of five state 0s:
calculate_g(
  costmatrix = constant_costmatrix,
  tip_states = c("0", "0", "0", "0", "0")
)

# Create a Type II character costmatrix:
binary_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "unordered"
)

# Calculate g for the case of two state 0s and three state 1s:
calculate_g(
  costmatrix = binary_symmetric_costmatrix,
  tip_states = c("0", "0", "1", "1", "1")
)

# Create a Type III character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "unordered"
)

# Calculate g for the case of two state 0s and three state 1s and two state 2s:
calculate_g(
  costmatrix = unordered_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2")
)

# Create a Type IV character costmatrix:
linear_ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "ordered"
)

# Calculate g for the case of two state 0s and three state 1s and two state 2s:
calculate_g(
  costmatrix = linear_ordered_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2")
)

# Create a Type V character costmatrix:
nonlinear_ordered_costmatrix <- convert_adjacency_matrix_to_costmatrix(
  adjacency_matrix = matrix(
    data = c(
      0, 1, 0, 0,
      1, 0, 1, 1,
      0, 1, 0, 0,
      0, 1, 0, 0
    ),
    nrow = 4,
    dimnames = list(0:3, 0:3)
  )
)

# Calculate g for the case of two state 0s, three state 1s, two state 2s and one state 3:
calculate_g(
  costmatrix = nonlinear_ordered_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2", "3")
)

# Create a Type VI character costmatrix:
binary_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "irreversible"
)

# Calculate g for the case of two state 0s and three state 1s:
calculate_g(
  costmatrix = binary_irreversible_costmatrix,
  tip_states = c("0", "0", "1", "1", "1")
)

# Create a Type VII character costmatrix:
multistate_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "irreversible"
)

# Calculate g for the case of two state 0s and three state 1s and two state 2s:
calculate_g(
  costmatrix = multistate_irreversible_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2")
)

# Create a Type VIII character costmatrix:
binary_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "dollo"
)

# Calculate g for the case of two state 0s and three state 1s:
calculate_g(
  costmatrix = binary_dollo_costmatrix,
  tip_states = c("0", "0", "1", "1", "1")
)

# Create a Type IX character costmatrix:
multistate_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "dollo"
)

# Calculate g for the case of two state 0s and three state 1s and two state 2s:
calculate_g(
  costmatrix = multistate_dollo_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2")
)

# Create a Type X character costmatrix:
multistate_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 5,
  character_type = "ordered"
)
multistate_symmetric_costmatrix$type <- "custom"
multistate_symmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 2, 3, 2, 3,
    1, 0, 3, 2, 1, 2,
    2, 3, 0, 3, 2, 1,
    3, 2, 3, 0, 1, 2,
    2, 1, 2, 1, 0, 1,
    3, 2, 1, 2, 1, 0
  ),
  nrow = multistate_symmetric_costmatrix$size,
  ncol = multistate_symmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_symmetric_costmatrix$single_states,
    multistate_symmetric_costmatrix$single_states
  )
)

# Calculate g for the case of two state 0s, three state 1s, two state 2s,
# one state 3, three state 4s and two state 5s:
calculate_g(
  costmatrix = multistate_symmetric_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2", "3", "4", "4", "4", "5", "5")
)

# Create a Type XI character costmatrix:
binary_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "ordered"
)
binary_asymmetric_costmatrix$type <- "custom"
binary_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1,
    10, 0
  ),
  nrow = binary_asymmetric_costmatrix$size,
  ncol = binary_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    binary_asymmetric_costmatrix$single_states,
    binary_asymmetric_costmatrix$single_states
  )
)
binary_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Calculate g for the case of two state 0s and three state 1s:
calculate_g(
  costmatrix = binary_asymmetric_costmatrix,
  tip_states = c("0", "0", "1", "1", "1")
)

# Create a Type XII character costmatrix:
multistate_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "ordered"
)
multistate_asymmetric_costmatrix$type <- "custom"
multistate_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 1,
    1, 0, 1,
    10, 10, 0
  ),
  nrow = multistate_asymmetric_costmatrix$size,
  ncol = multistate_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_asymmetric_costmatrix$single_states,
    multistate_asymmetric_costmatrix$single_states
  )
)
multistate_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Calculate g for the case of two state 0s and three state 1s and two state 2s:
calculate_g(
  costmatrix = multistate_asymmetric_costmatrix,
  tip_states = c("0", "0", "1", "1", "1", "2", "2")
)

---

Calculate the maximum possible tree length, gmax, under parsimony

Description

Given a costmatrix and number of tips, calculates the longest possible tree length under maximum parsimony.

Usage

calculate_gmax(costmatrix, n_taxa, allow_zeroes = FALSE)

Arguments

costmatrix	An object of class costMatrix.
n_taxa	The number of taxa (ree tips) to consider.
allow_zeroes	Logical indicating whether some states are allowed to be assigned zero tips. Defaults to FALSE (recommended). Note that if the number of states exceeds the number of tips then allow_zeroes = TRUE is forced.

Details

The concept of gmax was introduced by Hoyal Cuthill et al. (2010) and in the strictest sense represents the integer state frequency that maximizes g with the restriction that every state is sampled at least once. This function is intended to relax that slightly using the allow_zeroes option and indeed this is not possible in rare cases such as either high levels of missing data or gap-weighted (Thiele 1993) characters where it is possible for the number of states to exceed the (effective) number of tips. For a more detailed description of the problem and implemented solution(s) see Hoyal Cuthill and Lloyd (in press).

Value

A numeric indicating the maximum possible cost, gmax.

Author(s)

Graeme T. Lloyd [email protected] and Jen Hoyal Cuthill [email protected]

References

Hoyal Cuthill, J. F. and Lloyd, G. T., in press. Measuring homoplasy I: comprehensive measures of maximum and minimum cost under parsimony across discrete cost matrix character types. Cladistics, bold, .

Hoyal Cuthill, J. F., Braddy, S. J. and Donoghue, P. C. J., 2010. A formula for maximum possible steps in multistate characters: isolating matrix parameter effects on measures of evolutionary convergence. Cladistics, 26, 98-102.

Thiele, K.. 1993. The Holy Grail of the perfect character: the cladistic treatment of morphometric data. Cladistics, 9, 275-304.

See Also

calculate_g

Examples

# Create a Type I character costmatrix:
constant_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 0,
  character_type = "unordered"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = constant_costmatrix,
  n_taxa = 10
)

# Create a Type II character costmatrix:
binary_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "unordered"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = binary_symmetric_costmatrix,
  n_taxa = 10
)

# Create a Type III character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "unordered"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = unordered_costmatrix,
  n_taxa = 10
)

# Create a Type IV character costmatrix:
linear_ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "ordered"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = linear_ordered_costmatrix,
  n_taxa = 10
)

# Create a Type V character costmatrix:
nonlinear_ordered_costmatrix <- convert_adjacency_matrix_to_costmatrix(
  adjacency_matrix = matrix(
    data = c(
      0, 1, 0, 0,
      1, 0, 1, 1,
      0, 1, 0, 0,
      0, 1, 0, 0
    ),
    nrow = 4,
    dimnames = list(0:3, 0:3)
  )
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = nonlinear_ordered_costmatrix,
  n_taxa = 10
)


# Create a Type VI character costmatrix:
binary_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "irreversible"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = binary_irreversible_costmatrix,
  n_taxa = 10
)

# Create a Type VII character costmatrix:
multistate_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "irreversible"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = multistate_irreversible_costmatrix,
  n_taxa = 10
)

#' # Create a Type VIII character costmatrix:
binary_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "dollo"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = binary_dollo_costmatrix,
  n_taxa = 10
)

# Create a Type IX character costmatrix:
multistate_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "dollo"
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = multistate_dollo_costmatrix,
  n_taxa = 10
)

# Create a Type X character costmatrix:
multistate_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 5,
  character_type = "ordered"
)
multistate_symmetric_costmatrix$type <- "custom"
multistate_symmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 2, 3, 2, 3,
    1, 0, 3, 2, 1, 2,
    2, 3, 0, 3, 2, 1,
    3, 2, 3, 0, 1, 2,
    2, 1, 2, 1, 0, 1,
    3, 2, 1, 2, 1, 0
  ),
  nrow = multistate_symmetric_costmatrix$size,
  ncol = multistate_symmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_symmetric_costmatrix$single_states,
    multistate_symmetric_costmatrix$single_states
  )
)

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = multistate_symmetric_costmatrix,
  n_taxa = 10
)

# Create a Type XI character costmatrix:
binary_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "ordered"
)
binary_asymmetric_costmatrix$type <- "custom"
binary_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1,
    10, 0
  ),
  nrow = binary_asymmetric_costmatrix$size,
  ncol = binary_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    binary_asymmetric_costmatrix$single_states,
    binary_asymmetric_costmatrix$single_states
  )
)
binary_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = binary_asymmetric_costmatrix,
  n_taxa = 10
)

# Create a Type XII character costmatrix:
multistate_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "ordered"
)
multistate_asymmetric_costmatrix$type <- "custom"
multistate_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 1,
    1, 0, 1,
    10, 10, 0
  ),
  nrow = multistate_asymmetric_costmatrix$size,
  ncol = multistate_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_asymmetric_costmatrix$single_states,
    multistate_asymmetric_costmatrix$single_states
  )
)
multistate_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Calculate gmax for ten tips:
calculate_gmax(
  costmatrix = multistate_asymmetric_costmatrix,
  n_taxa = 10
)

---

Calculates a researcher's Kardashian Index

Description

Given counts of a researcher's Twitter followers and citations, returns their Kardashian Index.

Usage

calculate_kardashian_index(twitter_followers, total_citations)

Arguments

twitter_followers	The number of twitter followers the researcher has.
total_citations	The total number of citations across the researcher's publications (e.g., as garnered from a Google Scholar profile).

Details

This function implements the Kardashian Index of Hall (2014) and interested readers should consult that paper for more background.

Value

A scalar representing the ratio of expected Twitter followers (based on number of citations) to actual Twitter followers. Values greater than one indicate more Twitter followers than expected, those below one, fewer. According to Hall (2014), values above 5 are "Science Kardashians".

Author(s)

Graeme T. Lloyd [email protected]

References

Hall, N., 2014. The Kardashian index: a measure of discrepant social media profile for scientists. Genome Biology, 15, 424.

Examples

# Calculate the Kardashian Index of Sam Giles (@GilesPalaeoLab)
# as of 10/5/21:
calculate_kardashian_index(
  twitter_followers = 6534,
  total_citations = 550
)

# Calculate the Kardashian Index of Christopher Jackson (@seis_matters)
# as of 10/5/21:
calculate_kardashian_index(
  twitter_followers = 26000,
  total_citations = 6265
)

# Calculate the Kardashian Index of Graeme T. Lloyd (@GraemeTLloyd)
# as of 10/5/21:
calculate_kardashian_index(
  twitter_followers = 2133,
  total_citations = 2780
)

# Calculate the Kardashian Index of Katie Mack (@AstroKatie)
# as of 10/5/21:
calculate_kardashian_index(
  twitter_followers = 394900,
  total_citations = 1131
)

---

Get distance matrices from a cladistic matrix

Description

Takes a cladistic morphological dataset and converts it into a set of pairwise distances.

Usage

calculate_morphological_distances(
  cladistic_matrix,
  distance_metric = "mord",
  ged_type = "wills",
  distance_transformation = "arcsine_sqrt",
  polymorphism_behaviour = "min_difference",
  uncertainty_behaviour = "min_difference",
  inapplicable_behaviour = "missing",
  character_dependencies = NULL,
  alpha = 0.5
)

Arguments

cladistic_matrix	A character-taxon matrix in the format imported by read_nexus_matrix.
distance_metric	The distance metric to use. Must be one of "gc", "ged", "red", or "mord" (the default).
ged_type	The type of GED to use. Must be one of "legacy", "hybrid", or "wills" (the default). See details for an explanation.
distance_transformation	The type of distance transformation to perform. Options are "none", "sqrt", or "arcsine_sqrt" (the default). (Note: this is only really appropriate for the proportional distances, i.e., "gc" and "mord".)
polymorphism_behaviour	The distance behaviour for dealing with polymorphisms. Must be one of "mean_difference", "min_difference" (the default), or "random".
uncertainty_behaviour	The distance behaviour for dealing with uncertainties. Must be one of "mean_difference", "min_difference" (the default), or "random".
inapplicable_behaviour	The behaviour for dealing with inapplicables. Must be one of "missing" (default), or "hsj" (Hopkins and St John 2018; see details).
character_dependencies	Only relevant if using inapplicable_behaviour = "hsj". Must be a two-column matrix with colnames "dependent_character" and "independent_character" that specifies character hierarchies. See details.
alpha	The alpha value (sensu Hopkins and St John 2018). Only relevant if using inapplicable_behaviour = "hsj". See details.

Details

There are many options to consider when generating a distance matrix from morphological data, including the metric to use, how to treat inapplicable, polymorphic (e.g., 0&1), or uncertain (e.g., 0/1) states, and whether the output should be transformed (e.g., by taking the square root so that the distances are - or approximate - Euclidean distances). Some of these issues have been discussed previously in the literature (e.g., Lloyd 2016; Hopkins and St John 2018), but all likely require further study.

Claddis currently offers four different distance metrics: 1. Raw Euclidean Distance ("red") - this is only really applicable if there are no missing data, 2. The Gower Coefficient ("gc"; Gower 1971) - this rescales distances by the number of characters that can be coded for both taxa in each pairwise comparison thus correcting for missing data, 3. The Maximum Observable Rescaled Distance ("mord") - this was introduced by Lloyd (2016) as an extension of the "gc" designed to deal with the fact that multistate ordered characters can lead to "gc"s of greater than 1 and works by rescaling by the maximum possible distance that could be observed based on the number of characters codable in each pairwise comparison meaning all resulting distances are on a zero to one scale, and 4. The Generalised Euclidean Distance - this was introduced by Wills (1998) as a means of correcting for the fact that a "red" metric will become increasingly non-Euclidean as the amount of missing data increases and works by filling in missing distances (for characters that are coded as missing in at least one taxon in the pairwise comparison) by using the mean pairwise dissimilarity for that taxon pair as a substitute. In effect then, "red" makes no consideration of missing data, "gc" and "mord" normalise by the available data (and are identical if there are no ordered multistate characters), and "ged" fills in missing distances by extrapolating from the available data.

Note that Lloyd (2016) misidentified the substitute dissimilarity for the "ged" as the mean for the whole data set (Hopkins and St John 2018) and this was the way the GED implementation of Claddis operated up to version 0.2. This has now been amended (as of version 0.3) so that the function produces the "ged" in the form that Wills (1998) intended. However, this implementation can still be accessed as the "legacy" option for ged_type, with "wills" being the WIlls (1998) implementation. An advantage of this misinterpreted form of the GED is that it will always return a complete pairwise distance matrix, however it is not recommended (see Lloyd 2016). Instead a third option for ged_type - ("hybrid") - offers the same outcome but only uses the mean distance from the entire matrix in the case where there are no codable characters in common in a pairwise comparison. This new hybrid option has not been used in a published study.

Typically the resulting distance matrix will be used in an ordination procedure such as principal coordinates (effectively classical multidimensional scaling where k, the number of axes, is maximised at N - 1, where N is the number of rows (i.e., taxa) in the matrix). As such the distance should be - or approximate - Euclidean and hence a square root transformation is typically applied (distance_transformation with the "sqrt" option). However, if applying pre-ordination (i.e., ordination-free) disparity metrics (e.g., weighted mean pairwise distance) you may wish to avoid any transformation ("none" option). In particular the MORD will only fall on a zero to one scale if this is the case. However, if transforming the MORD for ordination this zero to one property may mean the arcsine square root ("arcsine_sqrt" option) is preferred. (Note that if using only unordered multistate or binary characters and the "gc" the zero to one scale will apply too.)

An unexplored option in distance matrix construction is how to deal with polymorphisms (Lloyd 2016). Up to version 0.2 of Claddis all polymorphisms were treated the same regardless of whether they were true polymorphisms (multiple states are observed in the taxon) or uncertainties (multiple, but not all states, are posited for the taxon). Since version 0.3, however, these two forms can be distinguished by using the different #NEXUS forms (Maddison et al. 1997), i.e., (01) for polymorphisms and {01} for uncertainties and within Claddis these are represented as 0&1 or 0/1, respectively. Thus, since 0.3 Claddis allows these two forms to be treated separately, and hence differently (with polymorphism_behaviour and uncertainty_behaviour). Again, up to version 0.2 of Claddis no options for polymorphism behaviour were offered, instead only a minimum distance was employed. I.e., the distance between a taxon coded 0&1 and a taxon coded 2 would be the smaller of the comparisons 0 with 2 or 1 with 2. Since version 0.3 this is encoded in the "min_difference" option. Currently two alternatives ("mean_difference" and "random") are offered. The first takes the mean of each possible difference and the second simply samples one of the states at random. Note this latter option makes the function stochastic and so it should be rerun multiple times (for example, with a for loop or apply function). In general this issue (and these options) are not explored in the literature and so no recommendation can be made beyond that users should think carefully about what this choice may mean for their individual data set(s) and question(s).

A final consideration is how to deal with inapplicable characters. Up to version 0.2 Claddis treated inapplicable and missing characters the same (as NA values, i.e., missing data). However, since Claddis version 0.3 these can be imported separately, i.e., by using the "MISSING" and "GAP" states in #NEXUS format (Maddison et al. 1997), with the latter typically representing the inapplicable character. These appear as NA and empty strings (""), respectively, in Claddis format. Hopkins and St John (2018) showed how inapplicable characters - typically assumed to represent secondary characters - could be treated in generating distance matrices. These are usually hierarchical in form. E.g., a primary character might record the presence or absence of feathers and a secondary character whether those feathers are symmetric or asymmetric. The latter will generate inapplicable states for taxa without feathers and without correcting for this ranked distances can be incorrect (Hopkins and St John 2018). Unfortunately, however, the #NEXUS format (Maddison et al. 1997) does not really allow explicit linkage between primary and secondary characters and so this information must be provided separately to use the Hopkins and St John (2018) approach. This is done here with the character_dependencies option. This must be in the form of a two-column matrix with column headers of "dependent_character" and "independent_character". The former being secondary characters and the latter the corresponding primary character. (Note that characters are to be numbered across the whole matrix from 1 to N and do not restart with each block of the matrix.) If using inapplicable_behaviour = "hsj" the user must also provide an alpha value between zero and one. When alpha = 0 the secondary characters contribute nothing to the distance and when alpha = 1 the primary character is not counted in the weight separately (see Hopkins and St John 2018). The default value (0.5) offers a compromise between these two extremes.

Here the implementation of this approach differs somewhat from the code available in the supplementary materials to Hopkins and St John (2018). Specifically, this approach is incorporated (and used) regardless of the overriding distance metric (i.e., the distance_metric option). Additionally, the Hopkins and St John function specifically allows an extra level of dependency (secondary and tertary characters) with these being applied recursively (tertiary first then secondary). Here, though, additional levels of dependency do not need to be defined by the user as this information is already encoded in the character_dependencies option. Furthermore, because of this any level of dependency is possible (if unlikely), e.g., quarternary etc.

Value

distance_metric	The distance metric used.
distance_matrix	The pairwise distance matrix generated.
comparable_character_matrix	The matrix of characters that can be compared for each pairwise distance.
comparable_weights_matrix	The matrix of weights for each pairwise distance.

Author(s)

Graeme T. Lloyd [email protected] and Thomas Guillerme [email protected]

References

Gower, J. C., 1971. A general coefficient of similarity and some of its properties. Biometrika, 27, 857-871.

Hopkins, M. J. and St John, K., 2018. A new family of dissimilarity metrics for discrete character matrices that include inapplicable characters and its importance for disparity studies. Proceedings of the Royal Society of London B, 285, 20181784.

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.

Maddison, D. R., Swofford, D. L. and Maddison, W. P., 1997. NEXUS: an extensible file format for systematic information. Systematic Biology, 46, 590-621.

Wills, M. A., 1998. Crustacean disparity through the Phanerozoic: comparing morphological and stratigraphic data. Biological Journal of the Linnean Society, 65, 455-500.

Examples

# Get morphological distances for the Day et al. (2016) data set:
distances <- calculate_morphological_distances(cladistic_matrix = day_2016)

# Show distance metric:
distances$distance_metric

# Show distance matrix:
distances$distance_matrix

# Show number of characters that can be scored for
# each pairwise comparison:
distances$comparable_character_matrix

# To repeat using the Hopkins and St John approach
# we first need to define the character dependency
# (here there is only one, character 8 is a
# secondary where 7 is the primary character):
character_dependencies <- matrix(c(8, 7),
  ncol = 2,
  byrow = TRUE, dimnames = list(
    c(),
    c(
      "dependent_character",
      "independent_character"
    )
  )
)

# Get morphological distances for the Day et
# al. (2016) data set using HSJ approach:
distances <- calculate_morphological_distances(
  cladistic_matrix = day_2016,
  inapplicable_behaviour = "hsj",
  character_dependencies = character_dependencies,
  alpha = 0.5
)

# Show distance metric:
distances$distance_metric

# Show distance matrix:
distances$distance_matrix

# Show number of characters that can be scored for
# each pairwise comparison:
distances$comparable_character_matrix

# Show weighting of each calculable pairwise distance:
distances$comparable_weights_matrix

---

Calculate mean pairwise distances

Description

Given distanceMatrices and taxonGroups objects calculates their mean pairwise distances.

Usage

calculate_MPD(distances, taxon_groups)

Arguments

distances	An object of class distanceMatrices.
taxon_groups	An object of class taxonGroups.

Details

Not all measures of disparity (morphological distance) require an ordination space. For example, the pariwise distances between taxa are themselves a disparity metric. This function takes the output from calculate_morphological_distances and a set of taxon groups and returns the mean pairwise distance for each of those groups.

Value

A labelled vector of weighted mean pairwise distances.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Get morphological distances for the Day et al. (2016) data set:
distances <- calculate_morphological_distances(
  cladistic_matrix = day_2016,
  distance_metric = "mord",
  distance_transformation = "none"
)

# Build simple taxonomic groups for Day et al. (2016) data set:
taxon_groups <- list(nonBurnetiamorpha = c("Biarmosuchus_tener", "Hipposaurus_boonstrai",
  "Bullacephalus_jacksoni", "Pachydectes_elsi", "Niuksenitia_sukhonensis", "Ictidorhinus_martinsi",
  "RC_20", "Herpetoskylax_hopsoni"), Burnetiamorpha = c("Lemurosaurus_pricei", "Lobalopex_mordax",
  "Lophorhinus_willodenensis", "Proburnetia_viatkensis", "Lende_chiweta",
  "Paraburnetia_sneeubergensis", "Burnetia_mirabilis", "BP_1_7098"))

# Set class as taxonGroups:
class(taxon_groups) <- "taxonGroups"

# Calculate mean pairwise distances:
calculate_MPD(distances, taxon_groups)

---

Calculates the parsimony length of a set of phylogenetic tree(s)

Description

Given a tree, or set of trees, and a cladistic matrix returns their parsimony length in number of steps.

Usage

calculate_tree_length(
  trees,
  cladistic_matrix,
  inapplicables_as_missing = FALSE,
  polymorphism_behaviour,
  uncertainty_behaviour,
  polymorphism_geometry,
  polymorphism_distance,
  state_ages,
  dollo_penalty
)

Arguments

trees	A tree (phylo object) or set of trees (multiPhylo object).
cladistic_matrix	A character-taxon matrix in the format imported by read_nexus_matrix. These should be discrete with rownames (ytaxon labels) matching the tip labels of trees.
inapplicables_as_missing	Logical that decides whether or not to treat inapplicables as missing (TRUE) or not (FALSE, the default and recommended option).
polymorphism_behaviour	One of either "missing", "uncertainty", "polymorphism", or "random". See details.
uncertainty_behaviour	One of either "missing", "uncertainty", "polymorphism", or "random". See details.
polymorphism_geometry	Argument passed to make_costmatrix.
polymorphism_distance	Argument passed to make_costmatrix.
state_ages	Argument passed to make_costmatrix.
dollo_penalty	Argument passed to make_costmatrix.

Details

Under the maximum parsimony criterion, a phylogenetic hypothesis is considered optimal if it requires the fewest number of evolutionary "steps" or - to generalise to non-discrete values - minimum total cost. In order to evalulate this criterion we must therefore be able to calculate a tree's "length" (total cost assuming the lowest cost for every character used). Given a set of phylogenetic hypothes(es) and a cladistic matrix this function calculates the minimum length for each tree.

Input data format

This function operates on a phylogenetic tree, or trees (in ape format), and a cladistic matrix (in cladisticMatrix format). However, the algorithm used is based on the generalised costmatrix approach of Swofford and Maddison (1992) and hence costmatrices need to be defined for each character (this is done internally by calling make_costmatrix), and some of the options are merely passed to this function.

Algorithm

Technically the Swofford and Maddison (1992) algorithm is designed for ancestral state reconstruction, but as its' first pass of the tree assigns lengths for each possible state at each node the minimum value of these options at the root is also the tree length for that character and hence by skipping the later steps this can be used as a tree length algorithm by simply summing the values across each character. The choice of the Swofford and Maddison algorithm, rather than the Wagner or Fitch algorithms (for ordered and unordered characters, respectively) is to generalize to the broadest range of character types, including asymmetric characters (Camin-Sokal, Dollo, stratigraphic), custom character types (specified using costmatrices or character state trees), as well as to any resolution of tree (i.e., including multifurcating trees - important for establishing maximum costs for homoplasy indices). The only restriction here is that the tree must be rooted such that time's arrow is explicitly present. This is essential, as the root defines the lengths across the whole tree, but also for asymmetric characters directionality must be explicit, as well as some downstream approaches (such as ACCTRAN and DELTRAN). The two obvious drawbacks to this algorithm are that it can be slower and that it is not appropriate for unrooted trees.

Costmatrices and costmatrix options

Costmatrices are described in detail in the make_costmatrix manual, as are the options that are passed from this function to that one. Thus, the user is directed there for a more in-depth discussion of options.

Inapplicable and missing characters

In practice these two character types are treated the same for length calculations - in effect these are "free" characters that do not constrain the tree length calculation in the same way that a coded character would (because a coded character's transition cost must be accounted for; Swofford and Maddison 1992). Note that there are reasons to take differences into account in phylogenetic inference itself (see papers by Brazeau et al. 2019 and Goloboff et al. in press). The option to treat them differently here is therefore only important in terms of downstream analyses, such as ancestral state reconstruction (see reconstruct_ancestral_states for details).

Polymorphisms and uncertainties

Polymorphisms (coded with empersands between states) and uncertainties (coded with slashes between states) can be interpreted in different ways, including those that affect estimates of tree length. Hence four options are provided to the user here:

1. Missing (polymorphism_behaviour = "missing" or uncertainty_behaviour = "missing"). Here polymorphisms are simply replaced by the missing character (NA). This removes polymorphisms and uncertainties from the calculation process completely (likely leading to undercounts), and hence is not generally recommended.

2. Uncertainty (polymorphism_behaviour = "uncertainty" or uncertainty_behaviour = "uncertainty"). This is the intended use of uncertain codings (e.g., 0/1) and constrains the tree length calculation to having to explain the least costly transition of those in the uncertainty. This is recommended for uncertain characters (although note that it biases the result towards the shortest possible length), but not truly polymorphic characters (as one or more state acquisitions are being missed, see Nixon and Davis 1991 and make_costmatrix for discussion of this). This is also - to the best of my knowledge - the approach used by most parsimony software, such as PAUP* (Swofford 2003) and TNT (Goloboff et al. 2008; Goloboff and Catalano 2016).

3. Polymorphism (polymorphism_behaviour = "polymorphism" or uncertainty_behaviour = "polymorphism"). If polymorphisms are real then some means of accounting for the changes that produce them seems appropriate, albeit difficult (see Nixon and Davis 1991 and Swofford and Maddison 1992 for discussions). If this option is applied it triggers the downstream options in make_costmatrix (by internally setting include_polymorphisms = TRUE), and the user should look there for more information. This is tentatively recommended for true polymorphisms (but note that it complicates interpretation), but not uncertainties.

4. Random (polymorphism_behaviour = "random" or uncertainty_behaviour = "random"). Another means of dealing with multiple-state characters is simply to sample a single state at random for each one, for example as Watanabe (2016) did with their PERDA algorithm. This simplifies the process, but also logically requires running the function multiple times to quantify uncertainty. This is not recommended for true polymorphisms (as interpretation is confounded), but may be appropriate for a less downwards biased tree count than "uncertainty".

These choices can also effect ancestral state estimation (see reconstruct_ancestral_states).

Polytomies

Polytomies are explicitly allowed by the function, but will always be treated as "hard" (i.e., literal multifurcations). Note that typically these will lead to higher tree lengths than fully bifurcating trees and indeed that the maximum cost is typically calculated from the star tree (single multifurcation).

Further constraints

In future the function will allow restrictions to be placed on the state at particular internal nodes. This can have multiple applications, including (for example) treating some taxa as ancestral such that their states are directly tied to specific nodes, e.g., in stratocladistics (Fisher 1994; Marcot and Fox 2008).

Character weights

Tree lengths output already include corrections for character weights as supplied in the cladistic_matrix input. So, for example, if a binary character costs two on the tree, but is weighted five then it will contribute a total cosr of 10 to the result.

Value

A list with multiple components, including:

input_trees	The tree(s) used as input.
input_matrix	The raw (unmodified) cladistic_matrix input.
input_options	The various input options used. Output for use by downstream functions, such as ancestral state estimation and stochastic character mapping.
costmatrices	The costmatrices (one for each character) used. These are typically generated automatically by the funcion, but are output here for later use in ancestral state estimation and stochastic character mapping functions.
character_matrix	The single character matrix object used. Essentially the input_matrix modified by the input_options.
character_lengths	A matrix of characters (rows) and trees (columns) with values indicating the costs. The column sums of this matrix are the tree_lengths values. This output can also be used for homoplasy metrics.
character_weights	A vector of the character weights used.
tree_lengths	The primary output - the length for each input tree in total cost.
node_values	The values (lengths for each state) for each node acrss trees and characters. This is used by reconstruct_ancestral_states for ancestral state reconstruction.

Author(s)

Graeme T. Lloyd [email protected]

References

Brazeau, M. D., Guillerme, T. and Smith, M. R., 2019. An algorithm for morphological phylogenetic analysis with inapplicable data. Systematic Biology, 68, 619-631.

Fisher, D. C., 1994. Stratocladistics: morphological and temporal patterns and their relation to phylogenetic process. In L. Grande and O. Rieppel (eds.), Interpreting the Hierarchy of Nature. Academic Press, San Diego. pp133–171.

Goloboff, P. A. and Catalano, S. A., 2016. TNT version 1.5, including a full implementation of phylogenetic morphometrics/ Cladistics, 32. 221-238

Goloboff, P., Farris, J. and Nixon, K., 2008. TNT, a free program for phylogenetic analysis. Cladistics, 24, 774-786.

Goloboff, P. A., De Laet, J., Rios-Tamayo, D. and Szumik, C. A., in press. A reconsideration of inapplicable characters, and an approximation with step‐matrix recoding. Cladistics.

Marcot, J. D. and Fox, D. L., 2008. StrataPhy: a new computer program for stratocladistic analysis. Palaeontologia Electronica, 11, 5A.

Nixon, K. C. and Davis, J. I., 1991. Polymorphic taxa, missing values and cladistic analysis. Cladistics, 7, 233-241.

Swofford, D. L., 2003. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts.

Swofford, D. L. and Maddison, W. P., 1992. Parsimony, character-state reconstructions, and evolutionary inferences. In R. L. Mayden (ed.) Systematics, Historical Ecology, and North American Freshwater Fishes. Stanford University Press, Stanford, p187-223.

Watanabe, A., 2016. The impact of poor sampling of polymorphism on cladistic analysis. Cladistics, 32, 317-334.

See Also

make_costmatrix, reconstruct_ancestral_states

Examples

# Use Gauthier 1986 as example matrix:
cladistic_matrix <- Claddis::gauthier_1986

# Use one of the MPTs from a TNT analysis as the tree:
tree <- ape::read.tree(
  text = paste(
    "(Outgroup,(Ornithischia,(Sauropodomorpha,(Ceratosauria,Procompsognathus,",
    "Liliensternus,(Carnosauria,(Ornithmimidae,Saurornitholestes,Hulsanpes,(Coelurus,",
    "Elmisauridae,(Compsognathus,(Ornitholestes,Microvenator,Caenagnathidae,",
    "(Deinonychosauria,Avialae))))))))));",
    sep = ""
  )
)

# Calculate tree length (and only use tree lengths from output):
calculate_tree_length(
  trees = tree,
  cladistic_matrix = cladistic_matrix,
  inapplicables_as_missing = TRUE,
  polymorphism_behaviour = "uncertainty",
  uncertainty_behaviour = "uncertainty",
  polymorphism_geometry = "simplex",
  polymorphism_distance = "euclidean",
  state_ages = c(200, 100),
  dollo_penalty = 999
)$tree_lengths

---

Calculate weighted mean pairwise distances

Description

Given distanceMatrices and taxonGroups objects calculates their weighted mean pairwise distances.

Usage

calculate_WMPD(distances, taxon_groups)

Arguments

distances	An object of class distanceMatrices.
taxon_groups	An object of class taxonGroups.

Details

Not all measures of disparity (morphological distance) require an ordination space. For example, the pariwise distances between taxa are themselves a disparity metric. However, due to variable amounts of missing data each pairwise distance should not necessarily be considered equal. Specifically, it could be argued that for a group of taxa the mean distance should be weighted by the number of characters that distance is based on, or more specifically the sum of the weights of those characters (e.g., Close et al. 2015).

This function takes the output from calculate_morphological_distances and a set of taxon groups and returns the weighted mean pairwise distance for those groups.

Value

A labelled vector of weighted mean pairwise distances.

Author(s)

Graeme T. Lloyd [email protected]

References

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.

Examples

# Get morphological distances for the Day et al. (2016) data set:
distances <- calculate_morphological_distances(
  cladistic_matrix = day_2016,
  distance_metric = "mord",
  distance_transformation = "none"
)

# Build simple taxonomic groups for Day et al. (2016) data set:
taxon_groups <- list(nonBurnetiamorpha = c("Biarmosuchus_tener", "Hipposaurus_boonstrai",
  "Bullacephalus_jacksoni", "Pachydectes_elsi", "Niuksenitia_sukhonensis", "Ictidorhinus_martinsi",
  "RC_20", "Herpetoskylax_hopsoni"), Burnetiamorpha = c("Lemurosaurus_pricei", "Lobalopex_mordax",
  "Lophorhinus_willodenensis", "Proburnetia_viatkensis", "Lende_chiweta",
  "Paraburnetia_sneeubergensis", "Burnetia_mirabilis", "BP_1_7098"))

# Set class as taxonGroups:
class(taxon_groups) <- "taxonGroups"

# Calculate mean pairiwise distances:
calculate_MPD(distances, taxon_groups)

# Now calculate weighted mean pairwise distances:
calculate_WMPD(distances, taxon_groups)

---

Check cladisticMatrix object for errors

Description

Internal function to check cladisticMatrix object for errors.

Usage

check_cladisticMatrix(cladistic_matrix)

Arguments

cladistic_matrix

An object of class cladisticMatrix.

Details

Internal Claddis function. Nothing to see here. Carry on.

Value

An error message or empty vector if no errors found.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Check that this is a valid cladisticMatrix object (will return error message as class
# is not set):
check_cladisticMatrix(cladistic_matrix = day_2016)

---

Check a costMatrix object for errors

Description

Internal function to check a costMatrix object for errors.

Usage

check_costMatrix(costmatrix)

Arguments

costmatrix

A costMatrix object.

Details

Costmatrix objects are more complex than what will typically be shown to the user. This function checks this hidden structure and reports any errors it finds.

These checks include rules 1-7 from Hoyal Cuthill and lloyd (i prep.).

Value

An error message or empty vector if no errors found.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Make an unordered costmatrix:
costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "unordered"
)

# Check that this is a valid costMatrix object (should return empty vector):
check_costMatrix(costmatrix = costmatrix)

---

Check a stateGraph object for errors

Description

Internal function to check a stateGraph object for errors.

Usage

check_stateGraph(stategraph)

Arguments

stategraph

A stateGraph object.

Details

Stategraph objects are more complex than what will typically be shown to the user. This function checks this hidden structure and reports any errors it finds.

These checks include rules 1-7 from Hoyal Cuthill and Lloyd (i prep.).

Value

An error message or empty vector if no errors found.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Make an unordered costmatrix:
costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "unordered"
)

# Convert costmatrix to stategraph:
stategraph <- convert_costmatrix_to_stategraph(costmatrix = costmatrix)

# Check that this is a valid stateGraph object (should return empty vector):
check_stateGraph(stategraph = stategraph)

---

Check taxonGroups object for errors

Description

Internal function to check taxonGroups object for errors.

Usage

check_taxonGroups(taxon_groups)

Arguments

taxon_groups

An object of class taxonGroups.

Details

Internal Claddis function. Nothing to see here. Carry on.

Value

An error message or empty vector if no errors found.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a taxon groups object:
taxon_groups <- list(
  Group_A = c("Species_1", "Species_2", "Species_3"),
  Group_B = c("Species_3", "Species_4"),
  Group_C = c("Species_5", "Species_6", "Species_7", "Species_8")
)

# Check that this is a valid taxonGroups object (will return error message as class
# is not set):
check_taxonGroups(taxon_groups = taxon_groups)

---

Check timeBins object for errors

Description

Internal function to check timeBins object for errors.

Usage

check_timeBins(time_bins)

Arguments

time_bins

A timeBins object.

Details

Internal Claddis function. Nothing to see here. Carry on.

Value

An error message or empty vector if no errors found.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a time bins object:
time_bins <- matrix(
  data = c(99.6, 93.5, 93.5, 89.3, 89.3, 85.8, 85.8, 83.5, 83.5, 70.6, 70.6, 65.5),
  ncol = 2,
  byrow = TRUE,
  dimnames = list(
    c("Cenomanian", "Turonian", "Coniacian", "Santonian", "Campanian", "Maastrichtian"),
    c("fad", "lad")
  )
)

# Check that this is a valid timeBins object (will return error message as class
# is not set):
check_timeBins(time_bins = time_bins)

---

Classify a costmatrix character

Description

Given a costmatrix, classifies it into one of twelve distinct types.

Usage

classify_costmatrix(costmatrix)

Arguments

costmatrix

An object of class costMatrix.

Details

Hoyal Cuthill and Lloyd (in press) classified discrete characters into twelve different types. This function applies their classification to a costmatrix object. The twelve different types are defined below alongside an empirical example for each, following Hoyal Cuthill and Lloyd (in review).

Type I - Constant

The simplest character type - only one state

Example: Character 18 of Sterli et al. (2013)

Costmatrix:

    -----
    | 0 |
---------
| 0 | 0 |
---------

State graph:

0

Type II - Binary symmetric

Exactly two states with equal and finite transition costs between them.

Example: Character 1 of Sterli et al. (2013)

Costmatrix:

    ---------
    | 0 | 1 |
-------------
| 0 | 0 | 1 |
-------------
| 1 | 1 | 0 |
-------------

State graph:

  1
0---1

Type III - Multistate unordered

Three or more states where every state-to-state transition has equal and finite cost.

Example: Character 53 of Sterli et al. (2013)

Costmatrix:

    -------------
    | 0 | 1 | 2 |
-----------------
| 0 | 0 | 1 | 1 |
-----------------
| 1 | 1 | 0 | 1 |
-----------------
| 2 | 1 | 1 | 0 |
-----------------

State graph:

    1
   / \
 1/   \1
 /     \
0-------2
    1

Type IV - linear ordered symmetric

Three or more states arranged as a single path with all connected states having equal and finite transition cost.

Example: Character 7 of Sterli et al. (2013)

Costmatrix:

    -------------
    | 0 | 1 | 2 |
-----------------
| 0 | 0 | 1 | 2 |
-----------------
| 1 | 1 | 0 | 1 |
-----------------
| 2 | 2 | 1 | 0 |
-----------------

State graph:

  1   1
0---1---2

Type V - Non-linear ordered symmetric

Four or more states arranged into a state graph with at least one branching point (vertex of degree three or greater) with all connected states having equal and finite transition costs.

Example: Character 71 of Hooker (2014)

Costmatrix:

    -------------
    | 0 | 1 | 2 | 3 |
---------------------
| 0 | 0 | 1 | 2 | 2 |
---------------------
| 1 | 1 | 0 | 1 | 1 |
---------------------
| 1 | 2 | 1 | 0 | 2 |
---------------------
| 1 | 2 | 1 | 2 | 0 |
---------------------

State graph:

          2
         /
       1/
   1   /
0-----1
       \
       1\
         \
          3

Type VI - Binary irreversible

Two states where transitions in one direction have infinite cost and the other finite cost.

Example: Character 119 of Gunnell et al. (2018)

Costmatrix:

    ---------
    | 0 | 1 |
-------------
| 0 | 0 | 1 |
-------------
| 1 | i | 0 |
-------------

(i = infinity)

State graph:

  1
0-->1

Type VII - Multistate irreversible

Three or more states where transitions in one direction have infinite cost and the other finite cost.

Example: Character 21 of Gunnell et al. (2018)

Costmatrix:

    -------------
    | 0 | 1 | 2 |
-----------------
| 0 | 0 | 1 | 2 |
-----------------
| 1 | i | 0 | 1 |
-----------------
| 2 | i | i | 0 |
-----------------

(i = infinity)

State graph:

  1   1
0-->1-->2

Type VIII - Binary Dollo

Two states where transitions in one direction can be made at most once, with no restriction in the other direction.

Example: Character 10 of Paterson et al. (2014)

Costmatrix:

    ----------
    | 0 | 1 |
-------------
| 0 | 0 | D |
-------------
| 1 | 1 | 0 |
-------------

(D = a Dollo penalty applied to avoid multiple transitions, see Swofford and Olsen 1990)

State graph:

      1
 --- <-- ---
| 0 |   | 1 |
 --- --> ---
      D

Type IX - Multistate Dollo

Three or more states where transitions in one direction can be made at most once, with no restriction in the other direction.

Example: Character 15 of Paterson et al. (2014)

Costmatrix:

    --------------
    | 0 | 1 |  2 |
------------------
| 0 | 0 | D | 2D |
------------------
| 1 | 1 | 0 |  D |
------------------
| 2 | 2 | 1 |  0 |
------------------

(D = a Dollo penalty applied to avoid multiple transitions, see Swofford and Olsen 1990)

State graph:

      1       1
 --- <-- --- <-- ---
| 0 |   | 1 |   | 2 |
 --- --> --- --> ---
      D       D

Type X - Multistate symmetric

Three or more states where connected states have symmetric finite costs, but where these are not all equal.

Example: Character 8 of Sumrall abd Brett (2002)

Costmatrix:

    -------------------------
    | 0 | 1 | 2 | 3 | 4 | 5 |
-----------------------------
| 0 | 0 | 1 | 2 | 3 | 2 | 3 |
-----------------------------
| 1 | 1 | 0 | 3 | 2 | 1 | 2 |
-----------------------------
| 2 | 2 | 3 | 0 | 3 | 2 | 1 |
-----------------------------
| 3 | 3 | 2 | 3 | 0 | 1 | 2 |
-----------------------------
| 4 | 2 | 1 | 2 | 1 | 0 | 1 |
-----------------------------
| 5 | 3 | 2 | 1 | 2 | 1 | 0 |
-----------------------------

State graph:

   1
 2---5   3
 |    \  |
2|    1\ |1
 |      \|
 0---1---4
   1   1

Type XI - Binary asymmetric

Two states where transition costs are finite but non-equal.

Example: Character 3 of Gheerbrant et al. (2014)

Costmatrix:

    ----------
    |  0 | 1 |
--------------
| 0 |  0 | 1 |
--------------
| 1 | 10 | 0 |
--------------

State graph:

      10
 --- <--- ---
| 0 |    | 1 |
 --- ---> ---
      1

Type XII - Multistate asymmetric

Three or more states where at least one transition is asymmetric, but the character does not meet the definition of other multistate asymmetric characters (i.e., Type VII and IX).

Example: Character 7 of Gheerbrant et al. (2014)

Costmatrix:

    ---------------
    |  0 |  1 | 2 |
-------------------
| 0 |  0 |  1 | 1 |
-------------------
| 1 |  1 |  0 | 1 |
-------------------
| 1 | 10 | 10 | 0 |
-------------------

State graph:

          10
    ______________
   /              \
  L   1       10   \
 --- <-- --- <--- ---
| 0 |   | 1 |    | 2 |
 --- --> --- ---> ---
  \   1        1   ^
   \______________/
           1

Other character types

The classification of Hoyal Cuthill and Lloyd (in review) was derived as a means to delineate a mathematical problem. As such it may be of use to consider where some special kinds of character would fall in this classification.

Stratigraphic characters

These would be classified as either Type VI, Type VII, Type XI, or Type XII characters (depending on state count and spacing of stratigraphic units).

Gap-weighted characters

A gap-weighted (Thiele 1993) character would be classified as either Type II or Type IV (depending on state count), but differ in the implicit assumption of whether all states are expected to be sampled. For a gap-weighted character there is no presumption that every state be sampled as the spacing between states is the primary information being considered.

Continuous characters

Continuous characters are necessarily not classifiable here as they are by definition non-discrete (except where translated into a gap-weighted character, above).

Value

A text string indicating the classification result (e.g., "Type VIII").

Author(s)

Graeme T. Lloyd [email protected]

References

Gunnell, G. F., Boyer, D. M., Friscia, A. R., Heritage, S., Manthi, F. K., Miller, E. R., Sallam, H. M., Simmons, N. B., Stevens, N. J. and Seiffert, E. R., 2018. Fossil lemurs from Egypt and Kenya suggest an African origin for Madagascar's aye-aye. Nature Communications, 9, 3193.

Hooker, J. J., 2014. New postcranial bones of the extinct mammalian family Nyctitheriidae (Paleogene, UK): primitive euarchontans with scansorial locomotion. Palaeontologia Electronica, 17.3.47A, 1-82.

Hoyal Cuthill, J. F. and Lloyd, G. T., in press. Measuring homoplasy I: comprehensive measures of maximum and minimum cost under parsimony across discrete cost matrix character types. Cladistics, bold, .

Paterson, A. M., Wallis, G. P., Kennedy, M. and Gray, R. D., 2014. Behavioural evolution in penguins does not reflect phylogeny. Cladistics, 30, 243-259.

Sterli, J., Pol, D. and Laurin, M., 2013. Incorporating phylogenetic uncertainty on phylogeny-based palaeontological dating and the timing of turtle diversification. Cladistics, 29, 233-246.

Sumrall, C. D. and Brett, C. E., 2002. A revision of Novacystis hawkesi Paul and Bolton 1991 (Middle Silurian: Glyptocystitida, Echinodermata) and the phylogeny of early callocystitids. Journal of Paleontology, 76, 733-740.

Swofford, D. L. and Olsen, G. J., 1990. Phylogeny reconstruction. In D. M. Hillis and C. Moritz (eds.), Molecular Systematics. Sinauer Associates, Sunderland. pp411-501.

Thiele, K.. 1993. The Holy Grail of the perfect character: the cladistic treatment of morphometric data. Cladistics, 9, 275-304.

Examples

# Create a Type I character costmatrix:
constant_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 0,
  character_type = "unordered"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = constant_costmatrix)

# Create a Type II character costmatrix:
binary_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 1,
  character_type = "unordered"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = binary_symmetric_costmatrix)

# Create a Type III character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "unordered"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = unordered_costmatrix)

# Create a Type IV character costmatrix:
linear_ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "ordered"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = linear_ordered_costmatrix)

# Create a Type V character costmatrix:
nonlinear_ordered_costmatrix <- convert_adjacency_matrix_to_costmatrix(
  adjacency_matrix = matrix(
    data = c(
      0, 1, 0, 0,
      1, 0, 1, 1,
      0, 1, 0, 0,
      0, 1, 0, 0
    ),
    nrow = 4,
    dimnames = list(0:3, 0:3)
  )
)

# Classify costmatrix:
classify_costmatrix(costmatrix = nonlinear_ordered_costmatrix)

# Create a Type VI character costmatrix:
binary_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "irreversible"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = binary_irreversible_costmatrix)

# Create a Type VII character costmatrix:
multistate_irreversible_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "irreversible"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = multistate_irreversible_costmatrix)

# Create a Type VIII character costmatrix:
binary_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "dollo"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = binary_dollo_costmatrix)

# Create a Type IX character costmatrix:
multistate_dollo_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "dollo"
)

# Classify costmatrix:
classify_costmatrix(costmatrix = multistate_dollo_costmatrix)

# Create a Type X character costmatrix:
multistate_symmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 5,
  character_type = "ordered"
)
multistate_symmetric_costmatrix$type <- "custom"
multistate_symmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 2, 3, 2, 3,
    1, 0, 3, 2, 1, 2,
    2, 3, 0, 3, 2, 1,
    3, 2, 3, 0, 1, 2,
    2, 1, 2, 1, 0, 1,
    3, 2, 1, 2, 1, 0
  ),
  nrow = multistate_symmetric_costmatrix$size,
  ncol = multistate_symmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_symmetric_costmatrix$single_states,
    multistate_symmetric_costmatrix$single_states
  )
)

# Classify costmatrix:
classify_costmatrix(costmatrix = multistate_symmetric_costmatrix)

# Create a Type XI character costmatrix:
binary_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 1,
  character_type = "ordered"
)
binary_asymmetric_costmatrix$type <- "custom"
binary_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1,
    10, 0
  ),
  nrow = binary_asymmetric_costmatrix$size,
  ncol = binary_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    binary_asymmetric_costmatrix$single_states,
    binary_asymmetric_costmatrix$single_states
  )
)
binary_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Classify costmatrix:
classify_costmatrix(costmatrix = binary_asymmetric_costmatrix)

# Create a Type XII character costmatrix:
multistate_asymmetric_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 2,
  character_type = "ordered"
)
multistate_asymmetric_costmatrix$type <- "custom"
multistate_asymmetric_costmatrix$costmatrix <- matrix(
  data = c(
    0, 1, 1,
    1, 0, 1,
    10, 10, 0
  ),
  nrow = multistate_asymmetric_costmatrix$size,
  ncol = multistate_asymmetric_costmatrix$size,
  byrow = TRUE,
  dimnames = list(
    multistate_asymmetric_costmatrix$single_states,
    multistate_asymmetric_costmatrix$single_states
  )
)
multistate_asymmetric_costmatrix$symmetry <- "Asymmetric"

# Classify costmatrix:
classify_costmatrix(costmatrix = multistate_asymmetric_costmatrix)

---

Collapses matrix to unique character state distributions

Description

Collapses a cladistic matrix to just unique character state distributions and taxon names.

Usage

compactify_cladistic_matrix(cladistic_matrix, message = TRUE)

Arguments

cladistic_matrix	The cladistic matrix in the format imported by read_nexus_matrix.
message	Logical indicating whether or not a message should be printed to the screen if the matrix cannot be compactified.

Details

Important: not recommended for general use.

This function is intended to make a matrix with redundant character state distributions smaller by collapsing these to single characters and upweighting them accordingly. It is intended purely for use with MRP matrices, but may have some very restricted uses elsewhere.

The function also deletes any characters weighted zero from the matrix and will merge duplicate taxon names into unique character strings.

Author(s)

Graeme T. Lloyd [email protected]

See Also

build_cladistic_matrix, prune_cladistic_matrix, read_nexus_matrix, safe_taxonomic_reduction, write_nexus_matrix, write_tnt_matrix

Examples

# Examine the matrix pre-compactification:
michaux_1989$matrix_1$matrix

# Examine the weights pre-compactification:
michaux_1989$matrix_1$character_weights

# Compactify the matrix:
michaux_1989compact <- compactify_cladistic_matrix(michaux_1989)

# Examine the matrix post-compactification:
michaux_1989compact$matrix_1$matrix

# Examine the weights post-compactification:
michaux_1989compact$matrix_1$character_weights

---

Converts an adjacency matrix to a costmatrix

Description

Takes an adjacency matrix as input and returns the corresponding costmatrix.

Usage

convert_adjacency_matrix_to_costmatrix(adjacency_matrix)

Arguments

adjacency_matrix

A labelled square matrix with zeroes denoting non-adjacencies and ones denoting adjacencies.

Details

This function is intended for internal use, but as it also generalizes to solving a general graph theory problem - generating a distance matrix corresponding to each shortest path between vertices of a connected graph represented as an adjacency matrix - it is made available explicitly here.

The process is best understood with an example. Imagine we have a graph like this:

  0-1-2
  |
  3

I.e., we have four labelled vertices, 0-3, and three edges (connections) between them:

  0-1
  1-2
  0-3

Note: here we assume symmetry, 0-1 = 1-0.

Graphs like this can be explicitly captured as adjacency matrices, where a one denotes two vertices are "adjacent" (connected by an edge) and a zero that they are not.

    _________________
    | 0 | 1 | 2 | 3 |
---------------------
| 0 | 0 | 1 | 0 | 1 |
---------------------
| 1 | 1 | 0 | 1 | 0 |
---------------------
| 2 | 0 | 1 | 0 | 0 |
---------------------
| 3 | 1 | 0 | 0 | 0 |
---------------------

But what such matrices do not tell us is how far every vertex-to-vertex path is in terms of edge counts. E.g., the path length from vertex 3 to vertex 2.

This function simply takes the adjacency matrix and returns the corresponding costmatrix, corresponding to every minimum vertex-to-vertex path length.

Value

An object of class costMatrix.

Author(s)

Graeme T. Lloyd [email protected]

See Also

convert_state_tree_to_adjacency_matrix, locate_bracket_positions

Examples

# Build the example adjacency matrix for the graph above:
adjacency_matrix <- matrix(
  data = c(
    0, 1, 0, 1,
    1, 0, 1, 0,
    0, 1, 0, 0,
    1, 0, 0, 0
  ),
  nrow = 4,
  ncol = 4,
  dimnames = list(0:3, 0:3)
)

# Convert this to a costmatrix:
convert_adjacency_matrix_to_costmatrix(
  adjacency_matrix = adjacency_matrix
)

---

Convert a costmatrix to a minimal state graph

Description

Given a costmatrix, returns the smallest possible state graph (fewest arcs).

Usage

convert_costmatrix_to_stategraph(costmatrix)

Arguments

costmatrix

An object of class costMatrix.

Details

A costmatrix summarises all possible state-to-state transition costs and hence each entry can also be considered as an arc of a directed state graph. However, many of these arcs could be removed and a complete description of the graph still be provided. For example, the diagonal (any transition from a state to itself - a loop) can be removed, as can any arc with infinite cost (as this arc would never be traversed in practice). Finally, some arcs are redundant as indirect paths already represent the same cost.

As an example, we can consider the linear ordered costmatrix:

    -------------
    | 0 | 1 | 2 |
-----------------
| 0 | 0 | 1 | 2 |
-----------------
| 1 | 1 | 0 | 1 |
-----------------
| 2 | 2 | 1 | 0 |
-----------------

A maximum directed graph representation would thus be:

----------------------
| from | to | weight |
----------------------
|   0  | 0  |   0    |
----------------------
|   0  | 1  |   1    |
----------------------
|   0  | 2  |   2    |
----------------------
|   1  | 0  |   1    |
----------------------
|   1  | 1  |   0    |
----------------------
|   1  | 2  |   1    |
----------------------
|   2  | 0  |   2    |
----------------------
|   2  | 1  |   1    |
----------------------
|   2  | 2  |   0    |
----------------------

But the following description is still complete, and minimal:

----------------------
| from | to | weight |
----------------------
|   0  | 1  |   1    |
----------------------
|   1  | 0  |   1    |
----------------------
|   1  | 2  |   1    |
----------------------
|   2  | 1  |   1    |
----------------------

This function effectively generates the latter (the minimal directed graph representation as a matrix of arcs).

Value

An object of class stateGraph.

Author(s)

Graeme T. Lloyd [email protected]

See Also

convert_adjacency_matrix_to_costmatrix

Examples

# Make a six-state unordered character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 5,
  character_type = "unordered"
)

# Find the minimal directed graph representation:
convert_costmatrix_to_stategraph(costmatrix = unordered_costmatrix)

# Make a six-state ordered character costmatrix:
ordered_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 5,
  character_type = "ordered"
)

# Find the minimal directed graph representation:
convert_costmatrix_to_stategraph(costmatrix = ordered_costmatrix)

# Make a six-state stratigraphic character costmatrix:
stratigraphic_costmatrix <- make_costmatrix(
  min_state = 0,
  max_state= 5,
  character_type = "stratigraphy",
  state_ages = c(103, 91.4, 78.2, 73.4, 66.0, 59.7)
)

# Find the minimal directed graph representation:
convert_costmatrix_to_stategraph(costmatrix = stratigraphic_costmatrix)

---

Converts a character state tree to an adjacency matrix

Description

Takes a character state tree as input and returns the corresponding adjacency matrix.

Usage

convert_state_tree_to_adjacency_matrix(state_tree)

Arguments

state_tree

A text string describing a character state tree.

Details

There are multiple ways to define discrete character states for use in phylogenetic analysis or ancestral state reconstruction in terms of transitions between states. The two most common being "ordered" ((Wagner optimization; Farris, 1970) or "unordered" (Fitch optimization; Fitch, 1971). These can be sonsidered as representing simple Markov models with symmetric transition probabilities, with the only difference being whether any direct state-tostae transition is possible (unordered) or only immediately "adjacent" states can be transitioned too.

In practice, however, there may be reason to believe in more complex transition models. For an example, see various characters in Hooker (2014). If these are still symmetric (i.e., going from state X to state Y is just as probable as going from state Y to state X) then they can be expressed using a "character state tree". For example, the model:

  0-1-2
  |
  3

This model has four states (labelled 0-3). However, it is clearly not an ordered model:

  0-1-2-3

Where only adjacent state transitions are possible, e.g., to get from state 0 to state 2 requires passing through state 1.

Nor is it an unordered model:

  0---1
  |\ /|
  | X |
  |/ \|
  2---3

Where any state-to-state transition is possible.

In other words we cannot simply label this model as "ordered" or "unordered". It can, however, be expressed directly as character state tree:

  (3,(2)1)0

Superficially these appear very similar to Newick strings (ways or representing bifurcating or multifurcating trees of species). However, the major difference is that with Newick strings we often only label the "tips" (terminal nodes, or vertices of degree one). Here there are no unlabelled "internal nodes", also known as "hypothetical ancestors". Instead, numbers after parentheses indicate the label of that node (vertex).

Thus we can read the above tree as node 2 is only directly connected to node 1, 3 and 1 are directly connected to node zero.

This is a very compact way of representing a Markov model, but unfortunately it is hard to read and difficult to do anything quantitative with. Thus this function takes such text strings and converts them to a more easily usable graph representation: the adjacency matrix.

Adjacency matrices represent graphs by showing links between vertices as values of one in symmetric square matrices. These always have a diagonal of zero (you cannot be adjacent to yourself) and any other vertices that are not directly linked are also coded as zero.

Note that here such matrices are also considered symmetric and hence the diagonal is also a line of reflection.

In practice more complex relationships may be desired, including differential weighting of specific transitions (without additional intermediate states; where going from state X to state Y has a cost other than zero or one), or asymmetric models (where going from state X to state Y has a different "cost" than going from state Y to state X). These are still possible in Claddis, and phylogenetic analysis in general, but require "costmatrices" to define them. (Character state trees are insufficent.)

Thus, costmatrices (or their probabilistic equivalent, Q-matrices) are the only generalisable form of defining any Markov model.

The output from this function can also be represented as a costmatrix with convert_adjacency_matrix_to_costmatrix.

Value

A symmetric adjacency matrix indicating which states are linked (value 1) by an edge.

Author(s)

Graeme T. Lloyd [email protected]

References

Farris, J. S., 1970. Methods for computing Wagner trees. Systematic Zoology, 19, 83-92.

Fitch, W. M., 1971. Towards defining the course of evolution: minimum change for a specified tree topology. Systematic Zoology, 20, 406-416.

Hooker, J. J., 2014. New postcranial bones of the extinct mammalian family Nyctitheriidae (Paleogene, UK): primitive euarchontans with scansorial locomotion. Palaeontologia Electronica, 17.3.47A, 1-82.

See Also

convert_adjacency_matrix_to_costmatrix, locate_bracket_positions

Examples

# Convert a simple state tree to an adjacency matrix:
convert_state_tree_to_adjacency_matrix(
 state_tree = "(3,(2)1)0"
)

# Convert a more complex state tree to an adjacency matrix:
convert_state_tree_to_adjacency_matrix(
  state_tree = "(((5)4)3,(2)1)0"
)

---

Convert a minimal state graph to a costmatrix

Description

Given a state graph returns the corresponding costmatrix.

Usage

convert_stategraph_to_costmatrix(stategraph)

Arguments

stategraph

An object of class stateGraph.

Details

A state graph describes the relationship between character states in terms of a graph of vertices (character states) and edges (weights). Edges can be "symmetric" (the graph is undirected) or "asymmetric" (the graph is directed, a digraph).

For example, a simple symmetric binary state graph might look ike this:

  1
A---B

Here the two states are A and B and the weight (cost of transitioning from A to B or B to A) is one.

In Claddis this graph can be represented using a stateGraph object including the arcs that describe the (di)graph:

from to weight
   0  1      1
   1  0      1

Each row represents an arc from one vertex to another, and the corresponding weight (transition cost). Note that for symmetric graphs the edge is stated using two arcs (from i to j and from j to i).

This function converts these state transitions to costmatrices, including interpolating any missing transitions by using the shortest indirect cost. In the example used here the costmatrix would look quite simple:

    ---------
    | 0 | 1 |
-------------
| 0 | 0 | 1 |
-------------
| 1 | 1 | 0 |
-------------

Value

An object of class costMatrix.

Author(s)

Graeme T. Lloyd [email protected]

See Also

convert_adjacency_matrix_to_costmatrix convert_costmatrix_to_stategraph

Examples

# Make state graph for a six-state linear ordered character:
stategraph <- list(
  n_vertices = 6,
  n_arcs = 10,
  n_states = 6,
  single_states = c("0", "1", "2", "3", "4", "5"),
  type = "ordered",
  arcs = data.frame(
    from = c("1", "0", "2", "1", "3", "2", "4", "3", "5", "4"),
    to = c("0", "1", "1", "2", "2", "3", "3", "4", "4", "5"),
    weight = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
  ),
  vertices = data.frame(
    label = c("0", "1", "2", "3", "4", "5"),
    in_degree = c(1, 2, 2, 2, 2, 1),
    out_degree = c(1, 2, 2, 2, 2, 1),
    eccentricity = c(5, 4, 3, 3, 4, 5),
    periphery = c(1, 0, 0, 0, 0, 1),
    centre = c(0, 0, 1, 1, 0, 0)
  ),
  radius = 3,
  diameter = 5,
  adjacency_matrix = matrix(
    data = c(
      0, 1, 0, 0, 0, 0,
      1, 0, 1, 0, 0, 0,
      0, 1, 0, 1, 0, 0,
      0, 0, 1, 0, 1, 0,
      0, 0, 0, 1, 0, 1,
      0, 0, 0, 0, 1, 0
    ),
    ncol = 6,
    byrow = TRUE,
    dimnames = list(
      c("0", "1", "2", "3", "4", "5"),
      c("0", "1", "2", "3", "4", "5")
    )
  ),
  directed = FALSE,
  includes_polymorphisms = FALSE,
  polymorphism_costs = "additive",
  polymorphism_geometry = "simplex",
  polymorphism_distance = "euclidean",
  includes_uncertainties = FALSE,
  pruned = FALSE,
  dollo_penalty = 999,
  base_age = 1,
  weight = 1
)

# Set calss as stateGraph:
class(x = stategraph) <- "stateGraph"

# View state graph:
stategraph

# Convert state graph to a costmatrix:
costmatrix <- convert_stategraph_to_costmatrix(stategraph = stategraph)

# Show costmatrix reflects linear ordered costs:
costmatrix$costmatrix

---

Counts the number of cherries in a tree

Description

Given a set of phylogenetic tree(s) returns the number of cherries in each one.

Usage

count_cherries(tree)

Arguments

tree	A tree (phylo or multiPhylo object).

Details

Cherries are components of a phylogenetic tree defined as internal nodes with exactly two terminal descendants.

This function simply counts the number present in a given tree.

Note that any fully dichotomous phylogenetic tree must have at least one cherry.

Value

Returns a vector of cherry counts for each tree retaining the order in which they were supplied.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create simple two-cherry tree:
tree <- ape::read.tree(text = "((A,B),(C,D));")

# Show count of cherries is two:
count_cherries(tree = tree)

# Create a star tree:
tree <- ape::read.tree(text = "(A,B,C,D);")

# Show count of cherries is zero:
count_cherries(tree = tree)

---

Returns node ages for a time-scaled tree

Description

Given a tree with branch-lengths scaled to time and a value for $root.time will return a vector of node ages.

Usage

date_nodes(time_tree)

Arguments

time_tree

A tree (phylo object) with branch lengths representing time and a value for $root.time.

Details

Returns a vector of node ages (terminal and internal) labelled by their node number.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create simple four-taxon tree with edge lengths all
# set to 1 Ma:
time_tree <- ape::read.tree(text = "(A:1,(B:1,(C:1,D:1):1):1);")

# Set root.time as 10 Ma:
time_tree$root.time <- 10

# Get node ages:
date_nodes(time_tree = time_tree)

---

Character-taxon matrix from Day et al. 2016

Description

The character-taxon matrix from Day et al. (2016).

Format

A character-taxon matrix in the format imported by read_nexus_matrix.

References

Day, M. O., Rubidge, B. S. and Abdala, F., 2016. A new mid-Permian burnetiamorph therapsid from the Main Karoo Basin of South Africa and a phylogenetic review of Burnetiamorpha. Acta Palaeontologica Polonica, 61, 701-719.

---

Drop tips from a time-scaled tree

Description

Drop tips from a time-scaled tree and update root.time accordingly

Usage

drop_time_tip(time_tree, tip_names, ...)

Arguments

time_tree	A time-scaled tree in phylo format where branch lengths are durations and where a $root.time value indicates the root age.
tip_names	A vector of tip names to be pruned from the tree.
...	Additional options to be passed to ape::drop.tip.

Details

(NB: This function is designed to only cope with trees containing at least three tips.)

Usually ape formatted trees are pruned with the drop.tip function in ape. However, trees time-scaled using either the paleotree or strap packages have an additional important component, the root age ($root.time) that may need updating when tips are removed. (See fix_root_time.) Thus this function is a modified version of drop.tip that also performs the fix_root_time step.

Note that dropPaleoTip in the paleotree package performs the exact same function, but is not called here to reduce the number of dependencies for Claddis.

Value

Returns a tree (phylo object) with pruned tips and corrected $root.time.

Author(s)

Graeme T. Lloyd [email protected]

See Also

fix_root_time

Examples

# Create a simple four-taxon tree with branch lengths:
tree <- ape::read.tree(text = "(A:1,(B:1,(C:1,D:1):1):1);")

# Set root age as 20 Ma:
tree$root.time <- 20

# Now prune taxon A:
pruned_tree <- ape::drop.tip(phy = tree, tip = "A")

# Show that drop.tip has not updated the tree's root time:
pruned_tree$root.time

# Use the function to fix the root time:
pruned_tree <- drop_time_tip(time_tree = tree, tip_names = "A")

# Show that the root time is now fixed (19 Ma):
pruned_tree$root.time

---

Ancestral Character State Estimation

Description

Given a tree and a cladistic matrix uses likelihood to estimate the ancestral states for every character.

Usage

estimate_ancestral_states(
  cladistic_matrix,
  time_tree,
  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
)

Arguments

cladistic_matrix	A character-taxon matrix in the format imported by read_nexus_matrix.
time_tree	A tree (phylo object) with branch lengths that represents the relationships of the taxa in cladistic_matrix.
estimate_all_nodes	Logical that allows the user to make estimates for all ancestral values. The default (FALSE) will only make estimates for nodes that link coded terminals (recommended).
estimate_tip_values	Logical that allows the user to make estimates for tip values. The default (FALSE) will only makes estimates for internal nodes (recommended).
inapplicables_as_missing	Logical that decides whether or not to treat inapplicables as missing (TRUE) or not (FALSE, the default and recommended option).
polymorphism_behaviour	One of either "equalp" or "treatasmissing".
uncertainty_behaviour	One of either "equalp" or "treatasmissing".
threshold	The threshold value to use when collapsing marginal likelihoods to discrete state(s).
all_missing_allowed	Logical to allow all missing character values (generally not recommended, hence default is FALSE).

Details

At its' core the function uses either the rerootingMethod (Yang et al. 1995) as implemented in the phytools package (for discrete characters) or the ace function in the ape package (for continuous characters) to make ancestral state estimates. For discrete characters these are collapsed to the most likely state (or states, given equal likelihoods or likelihood within a defined threshold value). In the latter case the resulting states are represented as an uncertainty (i.e., states separated by a slash, e.g., 0/1). This is the method developed for Brusatte et al. (2014).

The function can deal with ordered or unordered characters and does so by allowing only indirect transitions (from 0 to 2 must pass through 1) or direct transitions (from 0 straight to 2), respectively. However, more complex costmatrix transitions are not currently supported.

Ancestral state estimation is complicated where polymorphic or uncertain tip values exist. These are not currently well handled here, although see the fitpolyMk function in phytools for a way these could be dealt with in future. The only available options right now are to either treat multiple states as being equally probable of the "true" tip state (i.e., a uniform prior) or to avoid dealing with them completely by treating them as missing (NA) values.

It is also possible to try to use phylogenetic information to infer missing states, both for internal nodes (e.g., those leading to missing tip states) and for tips. This is captured by the estimate_all_nodes and estimate_tip_values options. These have been partially explored by Lloyd (2018), who cuationed against their use.

Value

The function will return the same cladistic_matrix, but with two key additions: 1. Internal nodes (numbered by ape formatting) will appear after taxa in each matrix block with estimated states coded for them, and 2. The time-scaled tree used will be added to cladistic_matrix as cladistic_matrix$topper$tree. Note that if using the estimate_tip_values = TRUE option then tip values may also be changed from those provided as input.

Author(s)

Graeme T. Lloyd [email protected] and Thomas Guillerme [email protected]

References

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.

Lloyd, G. T., 2018. Journeys through discrete-character morphospace: synthesizing phylogeny, tempo, and disparity. Palaeontology, 61, 637-645.

Yang, Z., Kumar, S. and Nei, M., 1995. A new method of inference of ancestral nucleotide and amino acid sequences. Genetics, 141, 1641-1650.

Examples

# Set random seed:
set.seed(4)

# Generate a random tree for the Day data set:
time_tree <- ape::rtree(n = nrow(day_2016$matrix_1$matrix))

# Update taxon names to match those in the data matrix:
time_tree$tip.label <- rownames(x = day_2016$matrix_1$matrix)

# Set root time by making youngest taxon extant:
time_tree$root.time <- max(diag(x = ape::vcv(phy = time_tree)))

# Use Day matrix as cladistic matrix:
cladistic_matrix <- day_2016

# Prune most characters out to make example run fast:
cladistic_matrix <- prune_cladistic_matrix(cladistic_matrix,
  characters2prune = c(2:3, 5:37)
)

# Estimate ancestral states:
estimate_ancestral_states(
  cladistic_matrix = cladistic_matrix,
  time_tree = time_tree
)

---

Estimate ancestral states for a continuous character under squared-change parsimony

Description

Given a phylogeny and set of continuous tip values, returns a single estimate of ancestral states that minimise the total squared-change cost on the tree.

Usage

estimate_squared_change_ancestors(tree, tip_values)

Arguments

tree	A tree (object of class phylo) with or without branch lengths and having tip labels that match the elements of tip_values.
tip_values	A labelled vector of continuous tip values, where the labels match the tip labels in tree. Tip values can be single numbers or a minimum-maximum range, see details.

Details

Multiple algorithms exist for estimating ancestral states for a univariate continuous character. This function specifically provides an ancestral state estimation that minimises the squared-change cost between each internal node and every other node that node is directly connected to. Note: in practice this approach can be identical to maximum likelihood ancestral state estimation under Brownian motion (e.g., Maddison 1991), leading Goloboff (2022) to argue that it isn't strictly a parsimony approach. However, this function extends squared-change parsimony to allow for ranges in tip values (see below) as well as implementing squared-change parsimony where branch lengths are variable and (hard) polytomies are permitted. As such it may still be of use to some users as distinct from other ancestral state estimation approaches for continuous characters.

Algorithm

Although squared-change parsimony ancestral state estimates can be directly calculated using the approach of Maddison (1991), the algorithm used here is an optimisation approach as this allows tip ranges to be accommodated (see below). In practice this leads to minimal speed reduction in most cases as an initial estimate is made using the fastAnc function from the phytools package (Revell 2024) that will frequently be identical to the final solution (except where ranges in tip values are used).

The optimisation used here is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm which was independently developed by Broyden (1970), Fletcher (1970), Goldfarb (1970), and Shanno (1970) and is implemented as the "L-BFGS-B" method in the R stats function optim(). The choice of this method is what allows ranges in tip values to be accommodated by treating these as an additional parameter to be estimated but with box constraints (i.e., minimum and maximum possible values - the ranges for the tip value) applied.

Ranges in tip values

Unlike other implementations of squared-change parsimony (to the best of my knowledge at least) this function removes the constraint that each tip must be represented by a single numeric value. Instead, any number of tip values can instead be represented by a range of values (i.e., a minimum and maximum value). These must be specified in the tip_values variable by separating the minimum and maximum values with an underscore character (_). E.g., "0.33_0.53" is the range 0.33 to 0.53. Note: an underscore is used rather than a dash as negative values are permitted and the dash symbol is reserved to indicate these.

In practice ranges are treated as an additional paramater to be estimated (alongside the internal node estimates) with the same criterion of minimising the total length, or cost (cum of squares), for the tree as a whole. As such this is truly a parsimony algorithm and negates the criticism of Goloboff (2022) that squared-change parsimony is not a true parsimony algorithm and also differs from (say) a maximum likelihood estimate that assumes a central tendency (i.e., midpoint) to the range of values at a tip.

The estimated tip values can also be examined a posteriori using the tip_values value from the function output.

Polytomies in the input tree

The function can also handle polytomies in the input tree, albeit it only does so by treating these as "hard" polytomies. I.e., it does not permute the minimum sum of squares for every possible resolution of a polytomy.

"Weighted" squared-change parsimony and branch lengths in the input tree

By default the function will apply so-called "unweighted" squared-change parsimony, meaning the user does not need to provide an input tree with branch lengths (only the topology). However, in practice "unweighted" parsimony really means treating every branch as having unit length (which is still a weighting). However, a user may wish to estimate ancestral states where variable branch lengths are accounted for (e.g., such that a cherry with variable terminal branch lengths means the subtending node is more strongly influenced by the tip value at the end of the shorter terminal).

In practice the function will implement weighted squared-change parsimony whenever the input tree already has edge-lengths applied. As such if the user does wish to implement "unweighted" squared-change parsimony they should be careful to supply an input tree where there are either no edge lengths or every edge is set to length one.

Value

The function returns a list with three elements:

ancestral_state_estimates	A single set of ancestral state estimates that represent a solution that minimises the total cost, or length, of the tree.
tip_values	The (estimated) tip values that minimise the sum of squares. This will be identical to the tip_values input if no ranges in tip values were used, but otherwise will show a single tip value that minimised the total cost, or length, of the tree.
sum_of_squares	A single numeric value indicating the total cost, or length, of the tree which is the sum of the squared length of each edge of the tree.

Author(s)

Graeme T. Lloyd [email protected]

References

Broyden, C. G., 1970. The convergence of a class of double-rank minimization algorithms. Journal of the Institute of Mathematics and Its Applications, 6, 76-90.

Fletcher, R., 1970. A new approach to variable metric algorithms. Computer Journal, 13, 317-322.

Goloboff, P. A., 2022. Phylogenetic Analysis of Morphological Data, Volume 2: Refining Phylogenetic Analyses. CRC Press, Boca Raton. 291 pp.

Goldfarb, D., 1970. A family of variable metric updates derived by variational means. Mathematics of Computation, 24, 23-26.

Maddison, W. P., 1991. Squared-change parsimony reconstructions of ancestral states for continuous-valued characters on a phylogenetic tree. Systematic Zoology, 40, 304-314.

Revell, L. J., 2024. phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and other things). PeerJ, 12, e16505.

Shanno, D. F., 1970. Conditioning of quasi-Newton methods for function minimization. Mathematics of Computation, 24, 647-656.

Examples

# Estimate squared-change parsimony for four tips with single values:
scp <- estimate_squared_change_ancestors(
  tree = ape::read.tree(text = "((A,B),(C,D));"),
  tip_values = tip_values <- c("A" = 3.6, "B" = 4.1, "C" = 7.3, "D" = 8.8)
)

# Plot results:
ape::plot.phylo(
  x = ape::read.tree(text = "((A,B),(C,D));"),
  main = paste("Sum of squares =", scp$sum_of_squares)
)
ape::tiplabels(text = scp$tip_values)
ape::nodelabels(text = round(x = scp$ancestral_state_estimates, digits = 2))

# Estimate squared-change parsimony for four tips with a ranged value:
scp <- estimate_squared_change_ancestors(
  tree = ape::read.tree(text = "((A,B),(C,D));"),
  tip_values = tip_values <- c("A" = "3.6", "B" = "4.1", "C" = "7.3", "D" = "8.8_9.3")
)

# Plot results:
ape::plot.phylo(
  x = ape::read.tree(text = "((A,B),(C,D));"),
  main = paste("Sum of squares =", scp$sum_of_squares)
)
ape::tiplabels(text = scp$tip_values)
ape::nodelabels(text = round(x = scp$ancestral_state_estimates, digits = 2))

# Estimate weighted squared-change parsimony for four tips with single values:
scp <- estimate_squared_change_ancestors(
  tree = ape::read.tree(text = "((A:1,B:2):2,(C:2,D:3):4);"),
  tip_values = tip_values <- c("A" = 3.6, "B" = 4.1, "C" = 7.3, "D" = 8.8)
)

# Plot results:
ape::plot.phylo(
  x = ape::read.tree(text = "((A:1,B:2):2,(C:2,D:3):4);"),
  main = paste("Sum of squares =", scp$sum_of_squares)
)
ape::tiplabels(text = scp$tip_values)
ape::nodelabels(text = round(x = scp$ancestral_state_estimates, digits = 2))

# Estimate squared-change parsimony for four tips with single values and a single polytomy:
scp <- estimate_squared_change_ancestors(
  tree = ape::read.tree(text = "(A,(B,C,D));"),
  tip_values = tip_values <- c("A" = 3.6, "B" = 4.1, "C" = 7.3, "D" = 8.8)
)

# Plot results:
ape::plot.phylo(
  x = ape::read.tree(text = "(A,(B,C,D));"),
  main = paste("Sum of squares =", scp$sum_of_squares)
)
ape::tiplabels(text = scp$tip_values)
ape::nodelabels(text = round(x = scp$ancestral_state_estimates, digits = 2))

---

Finds a minimum spanning tree of a costmatrix

Description

Given a costmatrix, returns a shortest tree connecting every state.

Usage

find_costmatrix_minimum_span(costmatrix)

Arguments

costmatrix

An object of class costMatrix.

Details

The minimum parsimony length a phylogenetic hypothesis can have depends on both the costmatrix of transition costs and the states actually sampled. If the costmatrix rows and columns already represent sampled states (the assumption here) then this minimum length is reduced to a graph theory problem - the minimum spanning tree or, if a directed graph, the minimum weight spanning arboresence. (NB: if there are unsampled states then the find_steiner_tree_of_stategraph function should be used instead.) This function returns one such shortest tree (although others may exist). The sum of the weights of the edges or arcs returned is the minimum cost.

As the algorithms used are graph theory based the function operates by simply calling convert_costmatrix_to_stategraph and find_stategraph_minimum_span. In practice, if the costmatrix represents a graph (transition costs are all symmetric) then Kruskal's algorithm is applied (Kruskal 1956). If costs are asymmetric, however, then the graph representation is a directed graph (or digraph) and so a version of Edmonds' algorithm is applied (Edmonds 1967).

Note that Dollo characters represent a special case solution as although a penalty weight is applied to the edges intended to only ever be traversed once this weight should not be used when calculating tree lengths. The function catches this and returns the edges with the weight that would actually be counted for a minimum weight spanning tree.

Value

A data.frame object describing a minimum spanning tree or minimum weight arboresence as a series of edges or arcs.

Author(s)

Graeme T. Lloyd [email protected]

References

Edmonds, J., 1967. Optimum branchings. Journal of Research of the National Bureau of Standards Section B, 71B, 233-240.

Kruskal, J. B., 1956. On the shortest spanning subtree of a graph and the traveling salesman problem. Proceedings of the American Mathematical Society, 7, 48-50.

See Also

find_shortest_costmatrix_path

Examples

# Make a four-state ordered character costmatrix:
ordered_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "ordered"
)

# Find length of shortest spanning tree of costmatrix:
find_costmatrix_minimum_span(costmatrix = ordered_costmatrix)

# Make a four-state unordered character costmatrix:
unordered_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "unordered"
)

# Find length of shortest spanning tree of costmatrix:
find_costmatrix_minimum_span(costmatrix = unordered_costmatrix)

# Make a four-state irreversible character costmatrix:
irreversible_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "irreversible"
)

# Find length of shortest spanning tree of costmatrix:
find_costmatrix_minimum_span(costmatrix = irreversible_costmatrix)

# Make a four-state Dollo character costmatrix:
dollo_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "dollo"
)

# Find length of shortest spanning tree of costmatrix:
find_costmatrix_minimum_span(costmatrix = dollo_costmatrix)

---

Gets descendant edges of an internal node

Description

Returns all descendant edges of an internal node for a phylo object.

Usage

find_descendant_edges(n, tree)

Arguments

n	An integer corresponding to the internal node for which the descendant edges are sought.
tree	A tree as a phylo object.

Details

Returns a vector of integers corresponding to row numbers in $edge or cells in $edge.length of the descendant edges of the internal node supplied.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create simple four-taxon tree:
tree <- ape::read.tree(text = "(A,(B,(C,D)));")

# Plot tree:
plot(tree)

# Show nodelabels:
nodelabels()

# Show edgelabels (note that edges 5 and 6
# are descendants of node 7):
edgelabels()

# Use find_descendant_edges to show that edges
# 5 and 6 are descendants of node 7:
find_descendant_edges(n = 7, tree = tree)

---

Find linked edges for a tree

Description

Given a tree finds edges that are linked to each other.

Usage

find_linked_edges(tree)

Arguments

tree	A tree (phylo object).

Details

Finds all edges that link (share a node) with each edge of a tree.

This is intended as an internal function, but may be of use to someone else.

Value

Returns a matrix where links are scored 1 and everything else 0. The diagonal is left as zero.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a simple four-taxon tree:
tree <- ape::read.tree(text = "(A,(B,(C,D)));")

# Find linked (1) edges matrix for tree:
find_linked_edges(tree)

---

Get edges of minimum spanning tree

Description

Returns edges of a minimum spanning tree given a distance matrix.

Usage

find_minimum_spanning_edges(distance_matrix)

Arguments

distance_matrix

A square matrix of distances between objects.

Details

This function is a wrapper for mst in the ape package, but returns a vector of edges rather than a square matrix of links.

Value

A vector of named edges (X->Y) with their distances. The sum of this vector is the length of the minimum spanning tree.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a simple square matrix of distances:
distance_matrix <- matrix(c(0, 1, 2, 3, 1, 0, 1, 2, 2, 1, 0, 1, 3, 2, 1, 0),
  nrow = 4,
  dimnames = list(LETTERS[1:4], LETTERS[1:4])
)

# Show matrix to confirm that the off diagonal has the shortest
# distances:
distance_matrix

# Use find_minimum_spanning_edges to get the edges for the minimum spanning
# tree:
find_minimum_spanning_edges(distance_matrix)

# Use sum of find_minimum_spanning_edges to get the length of the minimum
# spanning tree:
sum(find_minimum_spanning_edges(distance_matrix))

---

Find ancestor

Description

Finds the last common ancestor (node) of a set of two or more descendant tips.

Usage

find_mrca(descendant_names, tree)

Arguments

descendant_names	A vector of mode character representing the tip names for which an ancestor is sought.
tree	The tree as a phylo object.

Details

Intended for use as an internal function for trim_matrix, but potentially of more general use.

Value

ancestor_node

The ancestral node number.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a simple four-taxon tree:
tree <- ape::read.tree(text = "(A,(B,(C,D)));")

# Plot the tree:
ape::plot.phylo(tree)

# Add nodelabels and show that the most recent common
# ancestor of B, C, and D is node 6:
ape::nodelabels()

# Use find_mrca to show that the most recent common
# ancestor of B, C, and D is node 6:
find_mrca(
  descendant_names = c("B", "C", "D"),
  tree = tree
)

---

Finds the shortest path between two states in a costmatrix

Description

Given a start and end state, returns the shortest path through a costmatrix.

Usage

find_shortest_costmatrix_path(costmatrix, start, end)

Arguments

costmatrix	An object of class costMatrix.
start	The start state for the requested path.
end	The end state of the requested path.

Details

A common problem in graph theory is identifying the shortest path to take between two vertices of a connected graph. A costmatrix also describes a graph, with transition costs representing edge weights that can be asymmetric (e.g., going from 0 to 1 can have a different weight than going from 1 to 0). Finding the shortest path between states - i.e., the path that minimises total weight (cost) - from a costmatrix has two main applications in Claddis: 1. to check a costmatrix is internally consistent (no cost is misidentified due to a "cheaper" route being available - solving a problem identified in Maddison and Maddison 2003), and 2. to identify the minimum cost a character could have on a tree (an essential aspect of various homoplasy metrics, see Hoyal Cuthill et al. 2010).

The function returns a vector describing (one) shortest (i.e., lowest cost) path between start and end. If the direct path is shortest this will be simply start and end, but if an indirect route is cheaper then other node(s) will appear between these values.

In operation the function is inspired by Dijkstra's algorithm (Dijkstra 1959) but differs in some aspects to deal with the special case of a cladistic-style costmatrix. Essentially multiple paths are considered with the algorithm continuing until either the destination node (end) is reached or the accumulated cost (path length) exceeds the direct cost (meaning the path cannot be more optimal, lower cost, than the direct one).

Note: that because infinite costs are allowed in costmatrices to convey that a particular transition is not allowed these are always likely to be replaced (meanng the path does become possible) unless they apply to entire rows or columns of a costmatrix.

Note: negative costs are not allowed in costmatrices as they will confound the halting criteria of the algorithm. (They also do not make logical sense in a costmatrix anyway!)

Note: if multiple equally optimal solutions are possible, the function will only return one of them. I.e., just because a solution is not presented it cannot be assumed it is suboptimal.

Value

A vector of states describing (one) of the optimal path(s) in the order start to end.

Author(s)

Graeme T. Lloyd [email protected]

References

Dijkstra, E. W., 1959. A note on two problems in connexion with graphs. Numerische Mathematik, 1, 269–271.

Hoyal Cuthill, J. F., Braddy, S. J. and Donoghue, P. C. J., 2010. A formula for maximum possible steps in multistate characters: isolating matrix parameter effects on measures of evolutionary convergence. Cladistics, 26, 98-102.

Maddison, D. R. and Maddison, W. P., 2003. MacClade 4: Analysis of phylogeny and character evolution. Version 4.06. Sinauer Associates, Sunderland, Massachusetts.

See Also

convert_adjacency_matrix_to_costmatrix, make_costmatrix

Examples

# Make a four-state Dollo costmatrix:
costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 3,
  character_type = "dollo",
  dollo_penalty = 100
)

# Find the shortest path from state 0 to state 3:
find_shortest_costmatrix_path(
  costmatrix = costmatrix,
  start = "0",
  end = "3"
)

# Show that this is directional by asking for the reverse path:
find_shortest_costmatrix_path(
  costmatrix = costmatrix,
  start = "3",
  end = "0"
)

---

Finds a minimum spanning tree of a stategraph

Description

Given a stategraph, returns a shortest tree connecting every state.

Usage

find_stategraph_minimum_span(stategraph)

Arguments

stategraph

An object of class stateGraph.

Details

The minimum parsimony length a phylogenetic hypothesis can have depends on both the stategraph of transition costs and the states actually sampled. If the stategraph vertices already represent sampled states (the assumption here) then this minimum length is reduced to a graph theory problem - the minimum spanning tree or, if a directed graph, the minimum weight spanning arboresence. (NB: if there are unsampled states then the find_steiner_tree_of_stategraph function should be used instead.) This function returns one such shortest tree (although others may exist). The sum of the weights of the edges or arcs returned is the minimum cost.

As the algorithms used are graph theory based the function operates by simply calling convert_costmatrix_to_stategraph and find_stategraph_minimum_span. In practice, if the costmatrix represents a graph (transition costs are all symmetric) then Kruskal's algorithm is applied (Kruskal 1956). If costs are asymmetric, however, then the graph representation is a directed graph (or digraph) and so a version of Edmonds' algorithm is applied (Edmonds 1967).

Note that Dollo characters represent a special case solution as although a penalty weight is applied to the edges intended to only ever be traversed once this weight should not be used when calculating tree lengths. The function catches this and returns the edges with the weight that would actually be counted for a minimum weight spanning tree.

Value

A data.frame object describing a minimum spanning tree or minimum weight arboresence as a series of edges or arcs.

Author(s)

Graeme T. Lloyd [email protected]

References

Edmonds, J., 1967. Optimum branchings. Journal of Research of the National Bureau of Standards Section B, 71B, 233-240.

Kruskal, J. B., 1956. On the shortest spanning subtree of a graph and the traveling salesman problem. Proceedings of the American Mathematical Society, 7, 48-50.

See Also

find_shortest_costmatrix_path

Examples

# Make a four-state ordered character stategraph:
ordered_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "ordered"
)
ordered_stategraph <- convert_costmatrix_to_stategraph(costmatrix = ordered_costmatrix)

# Find length of shortest spanning tree of stategraph:
find_stategraph_minimum_span(stategraph = ordered_stategraph)

# Make a four-state unordered character stategraph:
unordered_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "unordered"
)
unordered_stategraph <- convert_costmatrix_to_stategraph(costmatrix = unordered_costmatrix)

# Find length of shortest spanning tree of stategraph:
find_stategraph_minimum_span(stategraph = unordered_stategraph)

# Make a four-state irreversible character stategraph:
irreversible_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "irreversible"
)
irreversible_stategraph <- convert_costmatrix_to_stategraph(costmatrix = irreversible_costmatrix)

# Find length of shortest spanning tree of stategraph:
find_stategraph_minimum_span(stategraph = irreversible_stategraph)

# Make a four-state Dollo character stategraph:
dollo_costmatrix <- make_costmatrix(
  min_state = "0",
  max_state = "3",
  character_type = "dollo"
)
dollo_stategraph <- convert_costmatrix_to_stategraph(costmatrix = dollo_costmatrix)

# Find length of shortest spanning tree of stategraph:
find_stategraph_minimum_span(stategraph = dollo_stategraph)

---

Find time bin midpoints

Description

Find the midpoint values for each bin from a timeBins object

Usage

find_time_bin_midpoints(time_bins)

Arguments

time_bins

A timeBins object.

Details

Frequently the midpoints of a series of time bins (defined by a beginning and ending) will be required, for example, when plotting binned data as a time series. Although the calculation involved is trivial (i.e., start date + end date / 2) this is a sufficiently common operation it is made into a formal function here.

Note that this function is designed to work specifically with objects of class "timeBins" - a format specific to Claddis that looks something like this:

                  fad  lad
    Cenomanian    99.6 93.5
    Turonian      93.5 89.3
    Coniacian     89.3 85.8
    Santonian     85.8 83.5
    Campanian     83.5 70.6
    Maastrichtian 70.6 65.5

I.e., a matrix with two columns (fad = first appearance date and lad = last appearance date) with rows corresponding to named time bins and indiviual values ages in millions of years ago (Ma). The object should also have class timeBins (see example below for hot to generate such an object). Note also that the convention here is to have time bins be ordered from oldest to youngest.

Value

A vector of time bin midpoint values.

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create a time bins object:
time_bins <- matrix(
  data = c(99.6, 93.5, 93.5, 89.3, 89.3, 85.8, 85.8, 83.5, 83.5, 70.6, 70.6, 65.5),
  ncol = 2,
  byrow = TRUE,
  dimnames = list(
    c("Cenomanian", "Turonian", "Coniacian", "Santonian", "Campanian", "Maastrichtian"),
    c("fad", "lad")
  )
)

# Set class as timeBins:
class(time_bins) <- "timeBins"

# Return midpoints for each time bin in sequence:
find_time_bin_midpoints(time_bins = time_bins)

---

Finds only the unique topologies amongst a set

Description

Given a set of trees with the same tip labels, returns just the unique topologies present.

Usage

find_unique_trees(trees)

Arguments

trees

An object of class multiPhylo.

Details

Where labelled topologies are generated randomly or modified by (e.g.) removing a tip, it may be useful to isolate just those that are truly unique. The ape package already has a function for this (unique.multiPhylo), but it can be slow when the number of trees is large. This function is thus intended as a faster version.

The function works by breaking down a tree into its' component bipartitions and treating the combination of these as the definition of the tree. It thus escapes problems due to the principle of free rotation. Specifically, these two trees are actually identical:

A  B  C   D  E
 \/    \   \/
  \     \  /
   \     \/
    \    /
     \  /
      \/

B  A  D  E   C
 \/    \/   /
  \     \  /
   \     \/
    \    /
     \  /
      \/

This becomes clearer if we decompose them into their bipartitions:

AB, DE, CDE, ABCDE

These correspond to the descendants of each internal node (branching point) and the last one is actually ignored (the root node) as it will be present in any tree.

Value

An object of class "multiPhylo".

Author(s)

Graeme T. Lloyd [email protected]

See Also

unique.multiPhylo

Examples

# Make a set of three identical trees (differing only in "rotation" of nodes):
trees <- ape::read.tree(text = c(
  "((A,B),(C,(D,E)));",
  "((C,(D,E)),(A,B));",
  "((B,A),(C,(E,D)));")
)

# Show that there is only one unique tree:
find_unique_trees(trees = trees)

---

Fixes a costmatrix that has inconsistent costs

Description

Given a costmatrix where transition costs are not self-consistent finds and returns a costmatrix that does.

Usage

fix_costmatrix(costmatrix, message = TRUE)

Arguments

costmatrix	A costMatrix object.
message	A logical indicating whether messages should be output (defaults to TRUE).

Details

A user may wish to consider a complicated set of transition costs between states when modelling discrete character evolution. This can be achieved with a custom costmatrix in Claddis (and elsewhere). However, some caution is urged when using such matrices to ensure that these costs are self-consistent (Maddison and Maddison 2003). More specifically, no direct state-to-state transition cost should be greater than is possible with an indirect path via one or more intermediate states.

This function offers a solution through an algorithm that will iteratively alter a costmatrix until all direct transition costs are self-consistent. It does so by finding the shortest state-to-state path for all possible transitions using the find_shortest_costmatrix_path function. Because the first solution may itself be inconsistent (as it relied on costs that have since updated) the algorithm is repeated until an equilibrium is reached. (This scenario is unlikely in most real world cases, but may be possible with very large matrices representing many states so was implemented here for safety.)

Note: infinite costs are allowed in costmatrices but unless they fill entire rows or columns (excluding the diagonal) they will not be self-consistent as there will always be a cheaper indirect cost.

Note: that both PAUP* (Swofford 2003) TNT (Goloboff et al. 2008; Goloboff and Catalano, 2016) offerthe same correction using the triangle inequality.

Note: other issues with a costmatrix may arise that are better revealed by using the check_costMatrix function, which returns informative error messages (with fix suggestions) where issues are found.

Value

A costMatrix object with self-consistent transition costs.

Author(s)

Graeme T. Lloyd [email protected]

References

Goloboff, P. A. and Catalano, S. A., 2016. TNT version 1.5, including a full implementation of phylogenetic morphometrics. Cladistics, 32, 221-238.

Goloboff, P., Farris, J. and Nixon, K., 2008. TNT, a free program for phylogenetic analysis. Cladistics, 24, 774-786.

Maddison, D. R. and Maddison, W. P., 2003. MacClade 4: Analysis of phylogeny and character evolution. Version 4.06. Sinauer Associates, Sunderland, Massachusetts.

Swofford, D. L., 2003. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts.

Examples

# Build a custom costmatrix with non-self consistent path lengths:
costmatrix <- make_costmatrix(
  min_state = 0,
  max_state = 2,
  character_type = "irreversible"
)
costmatrix$costmatrix[1:9] <- c(0, 2, 4, 1, 0, 3, 5, 3, 0)
costmatrix$symmetry <- "Asymmetric"
costmatrix$type <- "custom"

# Fix costmatrix:
fixed_costmatrix <- fix_costmatrix(costmatrix = costmatrix)

# Compare transition costs:
costmatrix$costmatrix
fixed_costmatrix$costmatrix

---

Fixes root.time after taxa have been pruned from a tree

Description

Fixes root.time after taxa have been pruned from a tree using ape::drop.tip

Usage

fix_root_time(original_tree, pruned_tree)

Arguments

original_tree	A tree in phylo format.
pruned_tree	A tree in phylo format that represents a pruned version of original_tree.

Details

(NB: This function is designed to only cope with trees containing at least three tips.)

When removing taxa from a time-scaled tree using drop.tip in ape $root.time is left unchanged. This can cause downstream problems if not fixed and that is what this function does.

Note that fix_root_time in the paleotree package performs the same function, but is not called here to reduce the number of libraries on which Claddis is dependent. Interested users should also refer to the dropPaleoTip function in paleotree.

Value

Returns a tree (phylo object) with a fixed $root.time.

Author(s)

Graeme T. Lloyd [email protected]

See Also

drop_time_tip

Examples

# Create a simple four-taxon tree with branch lengths:
tree <- ape::read.tree(text = "(A:1,(B:1,(C:1,D:1):1):1);")

# Set root age as 20 Ma:
tree$root.time <- 20

# Now prune taxon A:
pruned_tree <- ape::drop.tip(phy = tree, tip = "A")

# Show that drop.tip has not updated the tree's root time:
pruned_tree$root.time

# Use the function to fix the root time:
pruned_tree <- fix_root_time(original_tree = tree, pruned_tree = pruned_tree)

# Show that the root time is now fixed (19 Ma):
pruned_tree$root.time

---

Character-taxon matrix from Gauthier 1986

Description

The character-taxon matrix from Gauthier (1986).

Format

A character-taxon matrix in the format imported by read_nexus_matrix.

References

Gauthier, J. A., 1986. Saurischian monophyly and the origin of birds. In Padian, K. (ed.) The Origin of Birds and the Evolution of Flight. Towne and Bacon, San Francisco, CA, United States, 1-55.

---

Is a graph connected?

Description

Is a graph represented by an adjacenecy matrix connected?

Usage

is_graph_connected(adjacency_matrix)

Arguments

adjacency_matrix

An adjacency matrix where the diagonal is zeroes and the off-diagonal either ones (if the two vertices are directly connected) or zeroes (if not directly connected).

Details

Any undirected graph can be represented as an adjacency matrix and the properties of this matrix can be used to determine whether or not the graph is connected (i.e., a path between any two vertices exists) or not.

For example, the following graph:

6---4---5
    |   |\
    |   | 1
    |   |/
    3---2

Has the adjacency matrix:

    _________________________
    | 1 | 2 | 3 | 4 | 5 | 6 |
-----------------------------
| 1 | 0 | 1 | 0 | 0 | 1 | 0 |
-----------------------------
| 2 | 1 | 0 | 1 | 0 | 1 | 0 |
-----------------------------
| 3 | 0 | 1 | 0 | 1 | 0 | 0 |
-----------------------------
| 4 | 0 | 0 | 1 | 0 | 1 | 1 |
-----------------------------
| 5 | 1 | 1 | 0 | 1 | 0 | 0 |
-----------------------------
| 6 | 0 | 0 | 0 | 1 | 0 | 0 |
-----------------------------

This functions run through the following checks in order to confirm the connectivity of the graph:

1. As the graph has more than one vertex then further checks are required (a single vertex graph is considered connected).

2. As no vertice has degree zero (no links) then the graph may be connected (further checks are required).

3. As there are more than two vertices the graph may be disconnected (further checks are required).

4. As no missing paths are found (that would separate the graph into one or more disconnected subgraphs) then the graph must be connected.

This ordering means more complex queries are not triggered unless simpler tests do not provide a definitive answer.

Value

A logical (TRUE or FALSE).

Author(s)

Graeme T. Lloyd [email protected]

Examples

# Create the connected graph matrix:
x <- matrix(
  data = c(
    0, 1, 0, 0, 1, 0,
    1, 0, 1, 0, 1, 0,
    0, 1, 0, 1, 0, 0,
    0, 0, 1, 0, 1, 1,
    1, 1, 0, 1, 0, 0,
    0, 0, 0, 1, 0, 0
  ),
  ncol = 6,
  byrow = TRUE
)

# Check graph is connected:
is_graph_connected(adjacency_matrix = x)

---