Skip to contents

Compute Weighted Lineage Rates

Source: R/lineage_rates.R
lineage_rates.Rd

Compute the lineage-rate summary statistic described by Berv et al. (2026) from mapped-state phylo/simmap tree objects (Paradis and Schliep 2019; Revell 2012, 2024), including the mapped trees returned by bifrost. The function calculates present-biased weighted summaries of the state values inherited at ancestral nodes along each root-to-tip path. For bifrost BMM fits, the state values are the scalar regime-rate parameters stored in bifrost_search$model_no_uncertainty$param; for generic input, set bifrost_search = NULL and provide tree with a named state_values vector. Results include shift counts, terminal-state values, and the summary statistic, returned on the log scale by default.

Usage

lineage_rates(
  bifrost_search = NULL,
  tree = NULL,
  state_values = NULL,
  decay_base = 2,
  half_life = 5,
  normalize_weights = TRUE,
  age_reference = c("tree", "tip"),
  log = TRUE,
  cores = 1L,
  progress = interactive()
)

Arguments

A bifrost_search object returned by bifrost, a plain list with the required bifrost components, or NULL. bifrost inputs must include a fitted model_no_uncertainty component that identifies a multi-regime Brownian-motion model (model = "BMM") and provides named numeric regime-rate parameters in bifrost_search$model_no_uncertainty$param. Use NULL for generic input mode, where tree and state_values are required. Supply either bifrost_search or both tree and state_values, but not both input modes at once. Plain-list bifrost inputs must include tree_no_uncertainty_untransformed and model_no_uncertainty.

tree

SIMMAP-style phylo/simmap tree (Paradis and Schliep 2019; Revell 2012, 2024) required when bifrost_search = NULL for generic input mode. Must be NULL when bifrost_search is supplied. The tree must include SIMMAP maps state data from which node and tip states can be recovered. For the age-decay interpretation, the tree should be time-scaled; non-ultrametric trees are accepted, with age calculations controlled by age_reference. Internal nodes must use the conventional phylo numbering in which the root is Ntip(tree) + 1.

state_values

Named numeric vector required in generic input mode and invalid for bifrost inputs. Names must include every SIMMAP state label in the tree; extra names are ignored. Values must be finite, and must be strictly positive when log = TRUE.

decay_base

Numeric exponential-decay base b in the weighting kernel. Must be at least 1. Larger values make weights decline more strongly as elements get farther from the present; decay_base = 1 gives equal age weights regardless of half_life.

half_life

Optional decay time scale T, in the same time units as the tree. With the default decay_base = 2, T is a true half-life: an element T time units farther from the present receives half the unnormalized weight. More generally, weights are multiplied by 1 / decay_base over each interval of length T. Set half_life = NULL to use the unscaled kernel b^(-age); this only changes the weights when decay_base > 1.

normalize_weights

Logical; divide raw age-decay weights by their lineage-specific sum before averaging. Berv et al. (2026) used TRUE, which keeps the statistic as a weighted mean and makes lineages comparable when root-to-tip paths include different numbers of ancestral nodes. If every included ancestral node has the same state, the normalized value equals that state's supplied value or fitted rate.

age_reference

Character; reference point for converting node heights into ages before the present. The default, "tree", measures ages relative to the maximum tree height and preserves the original global-present behavior. Use "tip" for heterochronous or otherwise non-ultrametric trees when recency should be measured before each terminal tip's own sampling point. These choices are identical for ultrametric trees and mainly affect results when normalize_weights = FALSE.

log

Logical; return log-scale lineage-rate columns by default. bifrost BMM rates must always be finite and strictly positive. Generic state_values must be strictly positive when log = TRUE; when log = FALSE, generic values can remain on their original scale.

cores

Integer number of future workers. Values greater than one use future::multisession through future.apply::future_lapply(), which is platform agnostic.

progress

Logical; show a local progressr text progress bar for the tip-wise computations. Defaults to interactive(), so progress is shown in interactive sessions and suppressed in non-interactive runs.

Value

A data.frame with one row per tip. The stable top-level columns are tip.label, shift_count, the ancestral-node summary statistic, tip_state, and tip_rate. With log = TRUE, the summary statistic is named log_lineage_rate, and log_tip_rate is included. With log = FALSE, the summary statistic is named lineage_rate.

Additional diagnostic tables are stored as attributes. branch_metrics contains the per-tip intermediate table for the top-level ancestral-node statistic. Its columns are Tip, Shift_Count, Total_Time, Shift_Rate_Per_Time, Shift_Rate_Per_Speciation, Weighted_Lineage_Value, Tip_State, and Tip_State_Value. shift_metrics contains the same concise shift-count and time diagnostics without Weighted_Lineage_Value; its columns are Tip, Shift_Count, Total_Time, Shift_Rate_Per_Time, Shift_Rate_Per_Speciation, Tip_State, and Tip_State_Value. Shift_Count and both shift-rate columns are unweighted counts or ratios of detected parent-child node-state transitions. The settings attribute records the resolved input mode, model family, tree name, weighting parameters, age reference, log setting, and execution settings used for the calculation.

Details

Overview

The summary statistic follows Berv et al. (2026) as a descriptive, present-biased summary of each tip's inferred mapped-state history. In the default BMM use case, it summarizes heterogeneous phenotypic tempo along a tip-to-root path while giving more influence to recently inherited regimes. The calculation uses the node and tip states recovered from the supplied SIMMAP-style mapped-state tree. The parameters decay_base, half_life, and normalize_weights define the weighting kernel used to explore how inferred lineage summaries change as deeper history is weighted more or less strongly.

In the default bifrost mode, lineage_rates() reads regime-specific Brownian variance rates from bifrost_search$model_no_uncertainty$param and summarizes how those rates are inherited along mapped root-to-tip histories. That interpretation requires a heterogeneous Brownian-motion model: a multi-regime BMM fit with at least two mapped states.

Method and formula

Using the notation of Berv et al. (2026), generalized to an arbitrary decay base b and state-associated value \(r_i\), the normalized summary statistic is

$$ \bar{r}_{\mathrm{WLR}} = \frac{1}{Z} \sum_{i=1}^{L} b^{-a_i/T} \cdot r_i $$

where the normalizing constant is

$$ Z = \sum_{j=1}^{L} b^{-a_j/T}. $$

Here, \(\bar{r}_{\mathrm{WLR}}\) is the weighted mean lineage value, \(r_i\) is the state-associated value for node or path element \(i\), \(a_i\) is that element's age before the selected reference point, \(b\) is decay_base, \(T\) is the decay time scale, and \(L\) is the number of nodes or path elements included in the statistic. For the branch-based summary statistic aligned with Berv et al. (2026), \(i\) indexes ancestral nodes along the lineage, including the root and excluding the terminal tip state. This is the expected structure of bifrost mapped trees, where regime changes are represented by changes in node states rather than integrated over within-branch SIMMAP segments.

In bifrost output, \(r_i\) is represented by the scalar BMM regime-rate parameters in bifrost_search$model_no_uncertainty$param; following Berv et al. (2026), these correspond to the arithmetic mean of the diagonal variance terms in each regime's estimated evolutionary variance-covariance matrix. This is the weighted phenotypic-rate interpretation used in Berv et al. (2026). In generic input mode, \(r_i\) is the corresponding value from state_values.

Weighting, normalization, and age reference

Berv et al. (2026) used the base-2 case, and the defaults (decay_base = 2, half_life = 5, and normalize_weights = TRUE) preserve that weighting scheme. With these defaults, a regime that is half_life time units farther from the present receives half the unnormalized weight. decay_base = 1 gives every element raw weight 1 regardless of half_life, so the summary statistic becomes an equal-weighted mean when normalize_weights = TRUE. Values below 1 are rejected because they would give greater weight to regimes farther from the present.

The kernel terms \(b^{-a_i/T}\) are raw age-decay weights. When normalize_weights = TRUE, the function divides those raw weights by \(Z\) within each lineage before averaging. This keeps the statistic as a weighted mean of state-associated values, rather than letting lineages with more ancestral nodes accumulate larger values simply because more terms were summed. If every included ancestral node is assigned to the same state, the normalized value equals that state's supplied value or fitted rate regardless of how many nodes occur along the path. When normalize_weights = FALSE, the \(1/Z\) normalization is skipped.

The mapped tree should be time-scaled for the age-decay weights to have a direct temporal interpretation. The age_reference argument controls how the ages \(a_i\) are measured. With the default age_reference = "tree", ages are measured relative to the maximum tree height, matching the original global-present interpretation. For heterochronous or otherwise non-ultrametric trees, age_reference = "tip" measures recency before each terminal tip's own sampling point. These choices are identical for ultrametric trees. When normalize_weights = TRUE, they give the same top-level normalized values because all raw weights within a lineage are rescaled by the same constant before normalization. The choice mainly affects results when normalize_weights = FALSE.

When half_life = NULL, the unscaled kernel \(b^{-a_i}\) is used instead of \(b^{-a_i/T}\). This removes the half-life divisor from the exponent and only changes the weights when decay_base > 1. Parameter combinations that make a mathematically positive decay weight unrepresentable in double precision are rejected rather than returning a non-finite summary.

Interpretation

Interpret this as a flexible descriptive statistic for an inferred mapped history, not as a model-derived estimator with a unique optimal weighting scheme. It is not an instantaneous tip-rate estimate or a strict recent-time-window average. The default statistic is node-based: it uses ancestral node states along each root-to-tip path and does not integrate over within-edge SIMMAP segment durations. For SIMMAP-style trees with mapped changes inside branches, those segment durations are not part of the top-level statistic, and within-edge changes are not counted as shifts unless they are represented by parent-child node-state differences.

Lineage and shift diagnostics

The top-level log_lineage_rate or lineage_rate column is computed from ancestral node states along each root-to-tip path. Because the statistic is node-based and normalized within each lineage, its temporal resolution depends on the ages and density of ancestral nodes. Lineages with many recent nodes can downweight regimes farther from the present more rapidly, whereas sparse lineages may retain more influence from regimes farther from the present.

The branch_metrics and shift_metrics attributes store concise unweighted shift diagnostics: the number of detected parent-child node-state changes, root-to-tip time, shifts per unit time, and shifts per speciation event. branch_metrics additionally stores the weighted ancestral-node intermediate used for the top-level lineage-rate column. tip_rate is the fitted rate or supplied value for the terminal SIMMAP state. Compare the terminal value to the lineage-rate columns when you want to assess how much the inherited state history differs from the terminal regime state.

Generic input mode

lineage_rates() can also operate directly on a SIMMAP-style tree whose mapped node states refer to externally estimated values supplied through state_values. In this generic input mode, no fitted model is inspected: the same node-based summary statistic is applied to the user-supplied state values. This can be useful for summarizing arbitrary state-associated quantities on a mapped tree when the node-state interpretation is appropriate. The output column names retain the rate terminology used by the bifrost workflow for compatibility, but the values represent the supplied state_values and their interpretation depends on what those values measure. For generic SIMMAP trees, use this helper when ancestral node states are the summary target; it is not a duration-weighted integrator over mapped within-edge segments.

References

Berv, J. S. et al. (2026). Rates of passerine body plan evolution in time and space. Nature Ecology & Evolution. doi:10.1038/s41559-026-03110-5.

Clavel, J., Escarguel, G., and Merceron, G. (2015). mvMORPH: an R package for fitting multivariate evolutionary models to morphometric data. Methods in Ecology and Evolution, 6, 1311-1319. doi:10.1111/2041-210X.12420.

Clavel, J., Aristide, L., and Morlon, H. (2019). A penalized likelihood framework for high-dimensional phylogenetic comparative methods and an application to new-world monkeys brain evolution. Systematic Biology, 68, 93-116. doi:10.1093/sysbio/syy045.

Paradis, E., and Schliep, K. (2019). ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics, 35, 526-528. doi:10.1093/bioinformatics/bty633.

Revell, L. J. (2012). phytools: an R package for phylogenetic comparative biology (and other things). Methods in Ecology and Evolution, 3, 217-223. doi:10.1111/j.2041-210X.2011.00169.x.

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

Examples

# Generic input mode: summarize arbitrary state-associated values on a mapped tree.
toy_tree <- ape::read.tree(text = "((a:1,b:1):1,c:2);")
toy_tree <- phytools::paintSubTree(
  tree = toy_tree,
  node = ape::Ntip(toy_tree) + 1L,
  state = "0"
)
toy_tree <- phytools::paintSubTree(
  tree = toy_tree,
  node = ape::Ntip(toy_tree) + 2L,
  state = "1"
)
state_values <- c("0" = 1.2, "1" = 3.4)
lineage_rates(
  tree = toy_tree,
  state_values = state_values,
  log = FALSE,
  progress = FALSE
)
#>   tip.label shift_count lineage_rate tip_state tip_rate
#> 1         a           1     2.376124         1      3.4
#> 2         b           1     2.376124         1      3.4
#> 3         c           0     1.200000         0      1.2

if (FALSE) { # \dontrun{
# Search-result mode: use the multi-regime BMM fit and SIMMAP tree.
lineage_summary <- lineage_rates(bifrost_search = search_result)
} # }