Package 'treeshapesim'

Assign node ages to a topology from internode intervals

Description

Internal nodes are ranked by a uniform random linear extension of the tree order (the ranking a Yule process implies given the topology): the next node to split is drawn among the available ones with probability proportional to the number of internal nodes in its subtree.

Usage

assign_ages(topo, g)

Arguments

topo	Binary, rooted phylo topology.
g	Internode intervals, as from intervals_lineage_rate() (length Ntip(topo) - 1, positive and finite).

Value

topo with branch lengths; ultrametric, with root depth sum(g).

Examples

set.seed(1)
assign_ages(rtopology_beta(10), intervals_lineage_rate(10))

---

Calibrate the shape parameters for a number of tips

Description

Estimates, for a given number of tips, the mean gamma statistic as a function of c and the mean normalized Colless index (see tree_shape()) as a function of beta. Because the topology and the internode intervals are drawn independently, each curve is estimated from its own process alone: gamma from intervals drawn with intervals_lineage_rate() and the Colless index from topologies drawn with rtopology_beta(), without building trees. The expected values are monotone (gamma increases with c and the Colless index decreases with beta), but the estimated means need not be, especially for few tips or few draws, so each curve is made monotone by isotonic regression (stats::isoreg()) before it is inverted with shape_params().

Usage

calibrate_shape(
  n,
  beta = c(-1.95, -1.8, -1.6, -1.4, -1.2, -1, -0.75, -0.5, 0, 0.5, 1, 2, 4, 8),
  c = seq(-2, 1, by = 0.25),
  reps = 200
)

Arguments

n	Number of tips, a single integer of at least 4 (with 3 tips there is a single topology, so the Colless index cannot be calibrated).
beta	Values of beta at which the Colless index is estimated (at least two).
c	Values of c at which gamma is estimated (at least two).
reps	Number of topologies or interval sets per value.

Value

An object of class shape_calibration: a list with n, reps and the data frames colless (columns beta, colless and colless_mean) and gamma (columns c, gamma and gamma_mean), where the ⁠_mean⁠ columns hold the mean statistic at each value and colless and gamma their monotone fit, which shape_params() inverts.

Examples

set.seed(1)
cal <- calibrate_shape(50, reps = 50)
cal$gamma

---

Estimate the shape parameters of a tree

Description

Estimates, by maximum likelihood, the beta of Aldous' beta-splitting model from the topology and the c of the lineage-dependent pure-birth process from the branching times, which places an empirical tree on the same scale as the trees of sim_tree_shape(). Optionally, a parametric bootstrap with the package's own generators gives percentile intervals.

Usage

fit_shape(
  tree,
  boot = 0,
  level = 0.95,
  beta_range = c(-1.999, 10),
  c_range = c(-5, 5),
  tol = .Machine$double.eps^0.5
)

Arguments

tree	A binary, rooted, ultrametric phylo object with at least four tips.
boot	Number of parametric bootstrap replicates (0 for none).
level	Coverage of the bootstrap percentile intervals.
beta_range, c_range	Search intervals for beta (lower bound above -2) and c.
tol	Tolerance of the ultrametricity check, as in ape::is.ultrametric() (scaled range of root-to-tip distances). Trees with rounded branch lengths may need a larger value, or can be made ultrametric first (for instance with phytools::force.ultrametric()).

Details

beta. Each internal node with a clade of $m$ tips split into $i$ and $m - i$ tips contributes the probability of that split under the beta-splitting model (see rtopology_beta()), and the log-likelihood is the sum over nodes (Aldous 1996; this is the estimator applied by Blum & Francois 2006 and implemented in maxlik.betasplit() of the archived package apTreeshape). It is maximized over beta_range. In the package's validation runs (data-raw/validate_estimators.R in the source repository), the estimate was nearly unbiased for beta of -1 or less, but had a positive bias and a large variance for balanced trees with few tips (beta of 0 or more with 50 to 100 tips). An estimate at a bound of beta_range only says that the topology is extreme for the model.

c. With $k$ lineages, each lineage splits at rate $r k^c$, which is the density-dependent exponential (DDX) model of the R package laser (Rabosky 2006), whose per-lineage rate is $r k^{-x}$, so that $c = -x$. As in laser, the splits at 2 to $n - 1$ lineages are events and the last interval, from the most recent split to the present, is not; $r$ is profiled out (its maximum likelihood value is $(n - 2) / \sum_k k^{1 + c} g_k$, where $g_k$ is the time with $k$ lineages). The estimate describes the shape of the lineage-through-time curve, not its cause: branching times alone do not separate diversity-dependent from time-dependent diversification (Pannetier et al. 2021).

Bootstrap. With boot > 0, topologies are drawn with rtopology_beta() at the estimated beta and intervals with intervals_lineage_rate() at the estimated c, and both parameters are re-estimated, following the recommendation of a parametric bootstrap for diversification models by Etienne et al. (2016). The spread of the bootstrap estimates around the point estimate includes the small-sample bias of the estimators, and, for c, a small additional bias because intervals_lineage_rate() draws the last interval from the same law, whereas the likelihood treats it as censored. Percentile intervals need a few hundred replicates.

Value

An object of class shape_fit: a list with n, beta, c, r (speciation rate per lineage at one lineage, in the time units of the tree), loglik_beta (log-likelihood of beta up to an additive constant, so not comparable with other software), loglik_c (profile log-likelihood of c, equal to the log-likelihood of laser's DDX model at its maximum), beta_at_bound and c_at_bound (estimate at a bound of the search interval), beta_range, c_range, shape (output of tree_shape()) and, with boot > 0, boot (data frame of bootstrap estimates) and ci (matrix of percentile intervals).

References

Aldous, D. (1996). Probability distributions on cladograms. In D. Aldous & R. Pemantle (Eds.), Random discrete structures (pp. 1–18). Springer. doi:10.1007/978-1-4612-0719-1_1

Blum, M. G. B., & Francois, O. (2006). Which random processes describe the tree of life? A large-scale study of phylogenetic tree imbalance. Systematic Biology, 55, 685–691. doi:10.1080/10635150600889625

Etienne, R. S., Pigot, A. L., & Phillimore, A. B. (2016). How reliably can we infer diversity-dependent diversification from phylogenies? Methods in Ecology and Evolution, 7, 1092–1099. doi:10.1111/2041-210X.12565

Pannetier, T., Martinez, C., Bunnefeld, L., & Etienne, R. S. (2021). Branching patterns in phylogenies cannot distinguish diversity-dependent diversification from time-dependent diversification. Evolution, 75, 25–38. doi:10.1111/evo.14124

Rabosky, D. L. (2006). LASER: a maximum likelihood toolkit for detecting temporal shifts in diversification rates from molecular phylogenies. Evolutionary Bioinformatics, 2. doi:10.1177/117693430600200024

Examples

set.seed(1)
tr <- sim_tree_shape(200, beta = -1, c = -0.8)
fit_shape(tr)
fit <- fit_shape(tr, boot = 50)
fit$ci

---

Internode intervals of a pure-birth process with lineage-dependent rates

Description

With k lineages each splits at rate k^c, so the interval with k lineages is Exp(k^(1 + c)). c = 0 is Yule; c < 0 slows diversification as lineages accumulate (gamma < 0); c > 0 accelerates it (gamma > 0). The final interval (n lineages, up to the present) is drawn from the same law.

Usage

intervals_lineage_rate(n, c = 0)

Arguments

n	Number of tips, a single integer of at least 2.
c	Lineage dependence of the split rate, a single finite number.

Value

Intervals g_2, ..., g_n (length n - 1), scaled to sum 1.

References

Pybus, O. G., & Harvey, P. H. (2000). Testing macro-evolutionary models using incomplete molecular phylogenies. Proceedings of the Royal Society B, 267, 2267–2272. doi:10.1098/rspb.2000.1278

Examples

set.seed(1)
intervals_lineage_rate(10, c = -1)

---

Random topology under Aldous' beta-splitting model

Description

A clade of m tips splits into i and m - i tips with probability proportional to Gamma(beta + i + 1) Gamma(beta + m - i + 1) / (Gamma(i + 1) Gamma(m - i + 1)). beta = 0 is Yule, beta = -1.5 is PDA, beta close to -2 is caterpillar-like and large beta is balanced.

Usage

rtopology_beta(n, beta = 0)

Arguments

n	Number of tips, a single integer of at least 2.
beta	Aldous' beta, a single finite number greater than -2. Lower values give more imbalanced trees.

Value

A phylo object without branch lengths, with tips t1, ..., tn.

References

Aldous, D. (1996). Probability distributions on cladograms. In D. Aldous & R. Pemantle (Eds.), Random discrete structures (pp. 1–18). Springer. doi:10.1007/978-1-4612-0719-1_1

Examples

set.seed(1)
rtopology_beta(10, beta = -1)

---

Grid of tree shapes

Description

Crosses levels of the gamma statistic with levels of the normalized Colless index and finds the beta and c that give each combination on average, as a design for simulation studies. The levels can be chosen, for instance, from quantiles of tree_shape() over empirical trees of similar size.

Usage

shape_grid(n, gamma, colless, calibration = NULL, reps = 200)

Arguments

n	Number of tips, a single integer of at least 2.
gamma, colless	Levels of each statistic.
calibration	Output of calibrate_shape() for the same n; if NULL, it is computed with reps draws per value.
reps	Passed to calibrate_shape() when calibration is NULL.

Value

A data frame with one row per combination of levels (gamma, colless, beta, c), with the calibration in attribute calibration. Each row can be passed to sim_tree_target() (targets) or to sim_tree_shape() (parameters).

Examples

set.seed(1)
grid <- shape_grid(50, gamma = c(-2, -1, 0), colless = c(0, 0.5, 1), reps = 50)
grid

---

Shape parameters that give target mean values of gamma and Colless

Description

Inverts the monotone curves of calibrate_shape() by linear interpolation. Where a curve is flat (tied values of its monotone fit), the parameter is the mean of the tied values.

Usage

shape_params(calibration, gamma, colless)

Arguments

calibration	Output of calibrate_shape().
gamma, colless	Target mean gamma statistic and normalized Colless index (see tree_shape()); vectors that are recycled to a common length, as in data.frame().

Value

A data frame with the targets and the corresponding beta and c.

Examples

set.seed(1)
cal <- calibrate_shape(50, reps = 50)
shape_params(cal, gamma = c(0, -2), colless = 0.5)

---

Ultrametric trees with given imbalance and branching-time parameters

Description

Draws the topology with rtopology_beta() and the internode intervals with intervals_lineage_rate(), and combines them with assign_ages(). The topology and the intervals are drawn independently, the gamma statistic depends only on the intervals and the Colless index only on the topology, so for a given n the realized gamma depends only on c and the realized Colless index only on beta. To draw trees near target values of both statistics, see sim_tree_target().

Usage

sim_tree_shape(n, beta = 0, c = 0, nsim = 1, root_age = 1)

Arguments

n	Number of tips, a single integer of at least 2.
beta	Aldous' beta, a single finite number greater than -2. Lower values give more imbalanced trees.
c	Lineage dependence of the split rate, a single finite number.
nsim	Number of trees.
root_age	Age of the root (positive).

Value

With nsim = 1, a binary, ultrametric phylo object with root depth root_age; otherwise a multiPhylo object with nsim such trees.

Reproducibility

The same seed gives the same trees on the same platform. Between platforms (for instance, arm64 and x86-64), differences in the last bit of floating-point functions can change a few draws of the topology, so trees may differ.

Examples

set.seed(1)
tr <- sim_tree_shape(50, beta = -1, c = -1)
ape::gammaStat(tr)
trees <- sim_tree_shape(50, beta = -1, c = -1, nsim = 10, root_age = 30)

---

Simulate trees near target values of gamma and Colless

Description

Draws trees whose gamma statistic and normalized Colless index (see tree_shape()) fall within a tolerance of target values. The parameters beta and c whose mean statistics equal the targets are found with calibrate_shape() and shape_params(). Topologies are then drawn with rtopology_beta() and kept if their Colless index is within tol of the target, interval sets are drawn with intervals_lineage_rate() and kept if their gamma is within tol of the target, and the kept topologies and interval sets are combined with assign_ages(). Because the topology and the intervals are independent, rejecting them separately gives the same distribution of trees as rejecting whole trees on both statistics, with a much higher acceptance rate. Random numbers are used in a different order than in sim_tree_shape(), so the same seed does not give the same trees as that function.

Usage

sim_tree_target(
  n,
  gamma,
  colless,
  nsim = 1,
  tol = 0.5,
  calibration = NULL,
  reps = 200,
  max_tries = 1000 * nsim,
  root_age = 1
)

Arguments

n	Number of tips, a single integer of at least 2.
gamma, colless	Target gamma statistic and normalized Colless index (single values).
nsim	Number of trees.
tol	Tolerance around the targets: one value for both statistics, or two values for gamma and Colless, in this order.
calibration	Output of calibrate_shape() for the same n; if NULL, it is computed with reps draws per value.
reps	Passed to calibrate_shape() when calibration is NULL.
max_tries	Maximum number of topologies, and of interval sets, drawn.
root_age	Age of the root.

Value

With nsim = 1, a phylo object; otherwise a multiPhylo object (also when fewer trees than nsim are found). Attributes: shape (realized gamma and Colless of each tree), params (beta and c used) and acceptance (proportion of topologies and of interval sets kept). If fewer than nsim trees are found within max_tries, the trees found are returned with a warning.

Examples

set.seed(1)
cal <- calibrate_shape(50, reps = 50)
trees <- sim_tree_target(50, gamma = -1, colless = 0.5, nsim = 5,
                         calibration = cal)
attr(trees, "shape")

---

Gamma statistic and Colless index of a tree

Description

Computes Pybus and Harvey's (2000) gamma statistic, with ape::gammaStat(), and the Colless index normalized against its expectation under the Yule model, $(I_C - n \log n - n(\gamma_E - 1 - \log 2)) / n$, where $n$ is the number of tips and $\gamma_E$ is Euler's constant (the "yule" normalization of treestats::colless()). Both are near zero for Yule trees: gamma is negative when branching times are concentrated near the root, and the normalized Colless index is positive for trees more imbalanced than Yule trees.

Usage

tree_shape(tree)

Arguments

tree	A binary, rooted phylo object with at least three tips. The gamma statistic assumes an ultrametric tree; this is not checked.

Value

A named vector with gamma and colless.

References

Pybus, O. G., & Harvey, P. H. (2000). Testing macro-evolutionary models using incomplete molecular phylogenies. Proceedings of the Royal Society B, 267, 2267–2272. doi:10.1098/rspb.2000.1278

Examples

set.seed(1)
tree_shape(sim_tree_shape(100, beta = -1, c = -1))

---