R package · phylogenetic diversification

Inference of diversification with endogenous drivers

Monte Carlo likelihood inference for birth–death models whose speciation and extinction rates depend on covariates the diversification process generates itself: how many lineages there are, how old each one is, how distinct it is, and how many splits its ancestry carries.

Install GitHub Usage

Given a dated phylogeny of the surviving species of a clade, emphasis estimates the speciation and extinction rates that generated it, and lets those rates depend on the state of the clade itself and on the position of each lineage within it.

A reconstructed phylogeny records only the lineages that survive to the present. The lineages that went extinct carry the information about the extinction rate, and under a lineage-level covariate they also change the covariate of their living relatives, so the likelihood of the observed tree is an integral over every history of hidden lineages consistent with it. emphasis evaluates that integral by data augmentation with importance sampling, maximises it by Monte Carlo EM, and chooses among covariates either by comparing fitted models or along one dgLARS path over all of them.

Data augmentation. A proposal completes the observed tree (solid) with hidden lineages (dashed) that die before the present. Each completion is weighted by the model's density over the proposal's, and the weighted average over many completions estimates the likelihood of the observed tree. The hidden part is not a number but a tree, and its size varies from draw to draw.

Install

devtools::install_github("franciscorichter/emphasis")

Imports Rcpp, RcppParallel, nloptr, ape, progress and DDD. The likelihood, the simulator and the thinning proposal are C++ with TBB parallelism. Suggested: TreeSim for the bound finder's tests, dglars as the reference the selection path is checked against, ggplot2 and patchwork for the diagnostic plots, mgcv for the optional survival surface. Checked on Linux, macOS and Windows by R CMD check at every push.

The model

Time runs forward from the crown (age $0$) to the present (age $T$). A clade starts with two lineages. Every lineage $s$ alive at $t$ speciates at rate $\lambda_s(t)$ and goes extinct at rate $\mu_s(t)$, each a linear predictor over covariates passed through a link $g$:

\[ \lambda_s(t) = g\big(\beta_0 + \beta_N N(t) + \beta_D D_s(t) + \beta_{ED}\,\mathrm{ED}_s(t) + \beta_K K_s(t)\big), \]

with the same form and coefficients $\gamma$ for $\mu_s(t)$. A speciation replaces $s$ by two lineages, one of which keeps $s$'s identity; an extinction removes $s$. The covariates are read on the complete tree at $t$, hidden lineages included, which is what makes them endogenous: the process writes the covariates it then reads.

Covariates

Let $\tau_s$ be the time lineage $s$ was born or last split and $E_s(t) = t - \tau_s$ its pendant age, the length of its current terminal edge. With $N$ lineages alive their pendant ages sum to the pendant phylogenetic diversity $P = \sum_s E_s$ and average $M = P/N$.

Symbol	Name	Definition at time $t$	Reads
$N(t)$	diversity	the number of lineages alive	the branching times
$D_s(t)$	centred pendant age	$E_s(t) - M(t)$; sums to zero over the lineages alive	the topology: older or younger than its contemporaries
$\mathrm{ED}_s(t)$	evolutionary distinctiveness	fair proportion $\sum_{b \succeq s} \ell_b / n_b$ over the branches from the crown to $s$, each branch's length divided by the lineages alive below it; sums over the clade to Faith's phylogenetic diversity	the whole ancestry, weighted by how much of it $s$ shares
$\mathrm{ED}^c_s(t)$	centred distinctiveness	$\mathrm{ED}_s - \overline{\mathrm{ED}}$, the within-clade contrast alone	as $\mathrm{ED}$, with the clade-level part removed
$K_s(t)$	ancestral splits	the number of speciation events on the path from the crown to $s$; one more on both daughters at a split, hidden splits counted	whether a speciation-rich ancestry speciates more

Pendant age and its centring. Each bar is a lineage's pendant age $E_s$. $N$ counts the lineages; $M$ (dashed) is their mean, used only to centre $D_s = E_s - M$. Since $\sum_s D_s = 0$ at every instant, $D$ carries no clade-level signal and is orthogonal to $N$. $\mathrm{ED}_s$ adds the shared ancestry above the pendant edge; $K_s$ counts the splits along it.

What the basis spans

Measured over 6,433 lineages on 30 complete diversity-dependent clades: a lineage's pendant age is its nearest-neighbour distance (rank correlation 0.97), its time since the crown is the mirror of $\mathrm{ED}$ (−0.87), and $D$ and $\mathrm{ED}$ share the pendant edge (0.54). Two candidate axes were not spanned by $\{N, D, \mathrm{ED}\}$: the part of $\mathrm{ED}$ above the pendant edge, and the ancestral-split count $K$, whose correlation with each of the three is at most 0.29 in magnitude. $K$ is now a covariate. Under constant rates $K$ among the lineages alive at $t$ is Poisson with mean $2\lambda t$, so on one clade a $K$ effect and a rate rising with time are not told apart by $K$ alone.

Models

A model activates a subset of the covariates, named by a shortcut or a formula. Every model carries the same covariates in $\lambda$ and in $\mu$, with separate coefficients.

Shortcut	Formula	Parameters	Hypothesis
`"cr"`	`~ 1`	$(\beta_0, \gamma_0)$	constant rates, the null
`"dd"`	`~ N`	$(\beta_0, \beta_N, \gamma_0, \gamma_N)$	diversity dependence
`"d"`	`~ D`	$(\beta_0, \beta_D, \gamma_0, \gamma_D)$	age-dependent speciation
`"nd"`	`~ N + D`	6	age dependence over and above diversity
`"ed"`	`~ ED`	$(\beta_0, \beta_{ED}, \gamma_0, \gamma_{ED})$	distinctiveness-dependent diversification
`"ned"`	`~ N + ED`	6	distinctiveness over and above diversity
`"edc"`, `"nedc"`	`~ EDc`, `~ N + EDc`	4, 6	the same on the centred covariate
`"k"`	`~ K`	$(\beta_0, \beta_K, \gamma_0, \gamma_K)$	inherited diversification: a speciation-rich ancestry speciates more, or less
`"nk"`	`~ N + K`	6	the same over and above diversity

The parameter vector lists the $\beta$ coefficients of the active covariates in the order above, then the $\gamma$ coefficients: $(\beta_0, \beta_N, \beta_D, \gamma_0, \gamma_N, \gamma_D)$ for "nd", $(\beta_0, \beta_N, \beta_K, \gamma_0, \gamma_N, \gamma_K)$ for "nk". "ep" is an older name for "d". The models with $D$, $\mathrm{ED}$ or $K$ need the tree's topology (a phylo, not a vector of branching times), because the covariate is read per lineage. See the biological meaning of D for what $\beta_D$ tests.

Link functions

The predictor is mapped to a non-negative rate by link = "linear" | "exponential" | "gaussian":

Link functions, rate against the covariate predictor $\eta_{\text{cov}}$. Linear is the identity truncated at zero, exponential is log-linear, and the Gaussian is bell-shaped with its peak at $\eta_{\text{cov}} = 1$.

Link	Rate	Where it applies
Linear	$\max(0,\ \beta_0 + \eta_{\text{cov}})$	every model; a lineage whose rate is truncated to zero cannot split, so an observed split by such a lineage gives a draw density zero unless a hidden split preceded it
Exponential	$\exp(\beta_0 + \eta_{\text{cov}})$	every model; the only link on which the selection path is defined, since the point process is then log-linear
Gaussian	$\beta_0 \, \exp\!\big(-(\eta_{\text{cov}} - 1)^2 / 2\big)$	`cr`, `dd`, `d`, `nd`; not offered with $\mathrm{ED}$ or $K$

The likelihood

Write $y$ for the observed tree and $z$ for everything the observation removed: the hidden lineages with their birth times, attachments and death times, and the alive lineages that were not sampled. The pair $(y, z)$ is the complete tree, whose density has a closed form, and the likelihood of the data is the marginal

\[ L(\theta) = \int f(y, z \mid \theta)\, dz, \]

the sum over every hidden history consistent with the tree in hand. Under constant rates the integral collapses to Nee, May and Harvey's product of closed-form terms, and under diversity dependence to a system of differential equations in the number of lineages alive (the DDD package). A lineage-level covariate needs more than a count: a hidden lineage lowers the $\mathrm{ED}$ of its living relatives while it lives, a hidden split raises its descendants' $K$ and restarts its parent's pendant age, and no finite family of equations carries that state. The integral is evaluated by Monte Carlo.

For one complete tree the log-density is that of a marked point process: a log-rate for every event, read on the lineage it belongs to, minus the compensator, the integrated total rate of everything that could have happened and did not, plus the sampling terms:

\[ \log f(y, z \mid \theta) = \sum_{\text{spec}} \log \lambda_{s_i}(t_i) + \sum_{\text{ext}} \log \mu_{s_i}(t_i) - \int_0^T \sum_{s \in A(t)} \big(\lambda_s(t) + \mu_s(t)\big)\, dt + n_{\text{obs}} \log\rho + n_{\text{uns}} \log(1 - \rho). \]

Between two consecutive events the alive set $A(t)$ does not change, so the compensator is a sum over segments. On a segment $N$ and every $K_s$ are fixed, $M$ is held at its value at the segment's end, and $D_s$ and $\mathrm{ED}_s$ grow at unit rate, so each lineage's predictor is a line in $t$ whose integral has a closed form under both links. lineage_table() exposes these numbers as a data frame with one row per lineage and segment; rebuilding the log-likelihood from that table reproduces the C++ value to $10^{-8}$, and the same table is the design matrix of the selection path.

Importance sampling

Draw $z_1, \dots, z_K$ from a proposal $q(z \mid y, \theta')$ whose support contains every $z$ the model can produce, and weight each draw:

\[ w_k = \frac{f(y, z_k \mid \theta)}{q(z_k \mid y, \theta')}, \qquad \hat L(\theta) = \frac{1}{K} \sum_k w_k, \qquad \mathrm{ESS} = \frac{(\sum_k w_k)^2}{\sum_k w_k^2}. \]

$\hat L$ is unbiased whatever the proposal; its variance is what the proposal decides, and the effective sample size $\mathrm{ESS}$ is the number of equally weighted draws that would give the same variance. Every E-step reports it.

The reported log-likelihood is biased down

A fit reports $\hat\ell = \log \hat L$, the log of a Monte Carlo mean, and by Jensen's inequality $\mathbb{E}\,\hat\ell \le \log L$. The amount grows as the weights spread; a second-order expansion gives $\text{gap} = \tfrac12\,(1/\mathrm{ESS} - 1/K)$, which every fit records. AIC inherits it doubled, and it does not cancel between models: the model that is harder to sample is penalised for being harder to sample. A draw with $\mathrm{ESS} < K/10$ is flagged as heavy, and its gap is treated as unquantified rather than as the number returned. The M-step is not affected: its objective averages $\log f$ directly.

The two proposals

Both proposals insert hidden births into the observed tree at an intensity of the form $N(t)\,\bar\lambda(t) \times \Pr\{\text{a lineage born at } t \text{ leaves no sampled descendant}\}$, draw each inserted lineage's lifetime, and attach it to a parent among the lineages alive. They differ in what they put in the two factors.

Thinning (sampling = "dynamic_fresh") works in forward time at the model's own segment rates and draws births as an inhomogeneous Poisson process with intensity

\[ nh(t) = N(t)\,\lambda_{\text{seg}}\,\big(1 - \rho\, e^{-\mu_{\text{seg}}(T - t)}\big), \]

the last factor being the probability that a lineage born at $t$, with an exponential lifetime, is dead by $T$ or alive and unsampled. For a model with $\mathrm{ED}$ or $K$ it reads the per-lineage rates instead and attaches each birth in proportion to the parent's rate. It is available for every model and link.

The surrogate process (sampling = "bdi") starts from the exact conditional distribution of the hidden part under lineage-independent rates: with $p(t)$ the probability that a lineage alive at $t$ leaves a sampled descendant, solving $p' = \lambda p^2 - (\lambda - \mu)\,p$ backward from $p(T) = \rho$, the hidden lineages form a birth–death process with immigration in which the observed lineages emit hidden births at rate $2k\lambda(1 - p)$ and each hidden lineage splits at $\lambda(1 - p)$ and dies at $\mu / (1 - p)$. Under constant rates this is the model's own conditional and the weights are equal. Under diversity dependence and the lineage-level models the rates are replaced by those of the clade's average lineage on a deterministic trajectory $\hat N(t)$, $\hat P(t)$, the fixed point of a backward–forward iteration with under-relaxation; the weights carry the difference. The attachment of a hidden birth does not have to be mean-field: for a $D$ model each candidate parent is weighted by the model's own rate at that instant, and for an $\mathrm{ED}$ model by a pendant-age proxy of its distinctiveness.

Model	Surrogate proposal	Default
`cr`	exact on every link, $\mathrm{ESS} = K$	surrogate
`dd`	mean field on the linear and exponential links; the effective sample grows with the tree, 42 to 169 of 200 draws from 25 to 400 tips. Refused on the gaussian link, where the rate is not monotone in $N$ and the iteration does not contract past its peak	surrogate
`d`, `nd`	mean field with exact attachment, linear and exponential links. On 12 trees at 100 tips it leaves an effective sample of 2–439 against 1.3–15 for thinning; in fits it lifts 3–6 to 13–33 in four cells of five, with recovery of $\beta_D$ at 0.43–1.06 of its value and selection at 0.90 for an effect of −0.35 but 0.20 for −0.15	thinning; the surrogate on request
`ed`, `ned`, `edc`, `nedc`	mean field through the clade mean of fair proportion, linear and exponential links; 41 times the effective sample of thinning on `ned` trees at complete sampling, and ahead at every $\rho$ down to 0.2. On `ed` it trails thinning below $\rho \approx 0.8$	surrogate
`k`, `nk`	refused: the mean field carries no count of ancestral splits	thinning

The gate is applied before a fit runs, with a message naming the reason, and the fit records the sampler that ran. Supported is not the same as better, and the table says where the measured ranking turns.

Monte Carlo EM

At iterate $\theta^{(m)}$ the E-step draws $z_1, \dots, z_K$ from the proposal, computes the weights, normalises them to $\bar w_k$, and the M-step maximises the weighted complete-data log-likelihood over a box $B$:

\[ \theta^{(m+1)} = \arg\max_{\theta \in B}\ \sum_k \bar w_k \log f(y, z_k \mid \theta). \]

The draws are the same at every $\theta$ the M-step tries, so the per-lineage tables are built once per E-step and cached. The maximisation uses Rowan's Subplex through NLopt, with the first simplex sized at a tenth of the box in every coordinate. The box comes from the observed tree alone, by auto_bounds(): a centre at the net rate the tip count and crown age imply, then a bisection along each axis toward the point where forward simulation stops producing clades within a tenth and ten times the observed size, each slope axis scaled by its covariate's size on the tree.

Two stopping rules. rel_change stops once the relative step, in units that do not depend on the tree's time scale, falls below $10^{-2}$ on three consecutive iterations. mc_error measures the step against its own Monte Carlo standard error from five batches of the draws; a step inside the noise raises the number of draws by half, and the run stops only once the draws reach their cap and the step has stayed inside the noise for three iterations, returning the mean of the last three iterates. Both are subject to an iteration cap and a clock, and a fit reports which of the four ended it.
The last E-step. After the loop one more E-step runs at the returned point, so that the reported log-likelihood, effective sample and gap describe the estimate rather than the iterate before it.
Notes a fit carries. A lineage-level covariate on fewer than 50 tips, the smallest tree on which such an effect has been measured as selectable; an effective sample below 10 at the last E-step, which the pre-registered rules exclude from every selection statement; and turnover above 0.5 at the intercepts, the boundary of the measured regime. Each note states the number it rests on.

At the fixed point the iterates fluctuate around the maximiser with the Monte Carlo error of one E-step, measured at 1–7% of the rates over the last 100 of 400 iterations and halving per fourfold draws. No convergence rule tighter than that floor can fire, and the draws must grow with the tree: 200 fixed draws drift the estimates up by 8–13% at 400–800 tips, where draws grown to the Monte Carlo error do not.

Incomplete sampling and conditioning

When the phylogeny holds a fraction $\rho$ of the living species, set rho in the control list or in simulate_tree(). Both proposals then insert unsampled extant lineages as well as extinct ones, the likelihood gains the binomial factor $n_{\text{obs}} \log \rho + n_{\text{uns}} \log(1 - \rho)$, and $\rho$ enters the surrogate's survival equation only through its terminal condition $p(T) = \rho$. With rho = 1 the tree is completely sampled.

The likelihood is not conditioned on the clade's survival by default. The conditioned form $\ell_{\text{cond}} = \ell - \log P_\theta(\text{survival})$ is available by passing cond, a survival surface that auto_bounds() trains by forward simulation over the box when train_surv_gam = TRUE.

Three routes to a model

Where the routes leave the iteration. The dashed frame is one Monte Carlo EM iteration. The path leaves from the draws, so it needs no M-step and runs at any $\theta$; the comparison of fits leaves from the fit, one per model, whether the models are named one by one or as the ladder.

Route	Call	What it answers
One model	`estimate_rates(tree, model = "nd", ...)`	the maximum-likelihood estimate of a named model, with its log-likelihood, effective sample, gap and notes
The ladder	`estimate_rates(tree, model = "all", ...)`	fits `cr`, `dd`, `nd`, `ned`, `nk`, each started from the estimate of the model it nests and boxed by `auto_bounds()`, and returns them with their `compare_models()` table; `control$models` restricts the set
The path	`covariate_path(tree, covariates = c("N", "D", "ED", "K"))`	the order in which the covariates enter one dgLARS path over all of them, and the set a criterion picks, for $\lambda$ and $\mu$ separately; estimation and selection at once, on the exponential link

Comparing fits. $\mathrm{AIC} = 2k - 2\hat\ell$ sits above the truth by twice each fit's gap, and the amounts do not cancel. For every model behind the leader compare_models() computes the swing by which correcting both gaps would close the AIC difference, and flags the ranking as not safe where the swing reaches the difference, or where the trailing fit's draw is heavy. Starting a fit from its nested estimate is what keeps the first E-step's augmented clade near the size the data imply: from the middle of a wide box an $\mathrm{ED}$ fit can augment a 54-tip tree to thousands of nodes and run out of time, and from the dd estimate the same fit takes seconds.

The path. Under the exponential link the complete-data log-likelihood is a log-linear point process on the per-lineage table: one row per lineage and segment, the segment's length as exposure, a count of one when that lineage's own event ends it, and each draw's rows entering with its normalised weight. The differential-geometric LARS of Augugliaro, Mineo and Wit (2013) traces the curve on which the Rao score statistics of the active covariates stay equal, from the null fit where the first covariate enters down to the full fit; a covariate enters when its score reaches the common level, and BIC or AIC picks a point on the curve. On a fully observed tree the path is exact; on augmented trees its resolution is the augmentation's effective sample. It is checked against the dglars package where that package applies.

What is measured, and where

Every statement of accuracy and cost comes from a pre-registered simulation experiment in the companion papers, with the experiment named. The numbers below are the ones a user needs to read a fit. Each selection rate is a floor, because a fit stopped on its budget carries a worse AIC than it would reach.

Question	Measured
Does the estimator reach the exact maximum-likelihood estimate where one exists?	For `cr` and `dd`, against `DDD`, on 30 cells times two proposals and in a size study to 800 tips. The draws must grow with the tree.
How small a tree carries a lineage-level effect?	Selectable by AIC from 50 tips at an effect of 0.35 of the speciation rate and turnover 0.25. At 100 tips an effect of 0.35 is selected on 0.75 of trees and one of 0.15 on 0.50. Below 50 tips, above turnover 0.5, and for $K$, nothing is measured yet.
What does the path choose when the hidden lineages are drawn at the generating value?	At 50, 100, 200 and 400 tips: the null on `cr` trees 0.85, 0.90, 1.00, 1.00; exactly $N$ on `dd` trees 0.40, 0.55, 0.70, 0.65, with $N$ in the set on 0.85–1.00 at every size; exactly $\{N, \mathrm{ED}\}$ on `ned` trees 0.30, 0.20, 0.05, 0.00, on an effective sample of 5 above 50 tips. Where $D$ or $\mathrm{ED}$ generated the tree the path admits both, since the two share the pendant edge, and above 100 tips it prefers $D$. A path costs a median 1–3 min per tree at 50–100 tips, 14–19 min at 200 and 51–78 min at 400.
Is $K$ admitted where it is absent, and recovered where it is present?	On the path, never on 20 `cr`, 20 `dd` and 20 `nd` trees, once on 20 `ned` trees; on `nk` trees it is in the chosen set on 0.70, at 0.97 of its value, with $\mathrm{ED}$ in the set on every one, since a speciation-rich ancestry is a low fair proportion. By fitting, on 100 `nk` trees at 100 tips, $\beta_K$ is recovered at 0.83–0.95 of its value with the right sign on 0.80–1.00, and no false positive at zero; AIC chooses `nk` over `dd` on 0.60 of trees at an effect of +0.35 and on 0.15 or fewer elsewhere, on fits that ended on their 30-minute clock on 93 of 100 trees.
Does a lineage-level covariate fitted alone recover its coefficient?	On 20 trees a cell at 100 tips, effects ±0.35: `ed` at 0.77 and 1.05 of its value, `edc` at 0.96 and 1.02, `d` at 0.67 and 1.01, `k` at 0.87 and 0.58, with the right sign on every tree. Chosen over `cr` on 0.70–1.00 of trees, except `k` at a positive effect, where the thinning proposal keeps 2.6 effective draws. False positives at zero effect: none for `edc` and `k`, 0.10 for `ed`, 0.15 for `d`, bought by an extinction slope the generator did not have. The centred covariate is the one of the four that meets every half of the rule.
What does an iteration cost?	One `ned` iteration with a strong $\mathrm{ED}$ effect costs 77 s at 100 tips and 5.1 h at 400 from a start far from the estimate; the difference is the starting point, which is why the ladder starts each model from the one it nests.

The ladder itself, model = "all" on trees from each of its five models at 100 and 200 tips, is being measured now; this page is updated when that result is in.

Simulation

simulate_tree() runs the forward process by Gillespie simulation from the crown, carrying every lineage's pendant age, distinctiveness and split count, for every model and link. It returns the reconstructed tree (tes), the complete tree with its extinct lineages (tas) and a status of "done", "extinct" or "too_large". Given an observed tree it draws from either proposal instead, returning the augmented trees with their $\log q$.

sim <- simulate_tree(
  pars  = c(0.8, -0.02, 0.15, 0.1, 0, 0),   # beta_0, beta_N, beta_D, gamma_0, gamma_N, gamma_D
  max_t = 10,
  model = "nd",
  link  = "linear"                         # or "exponential", "gaussian"
)
sim$status
plot(sim$tes)                               # the reconstructed tree

# draw hidden lineages for an observed tree, from either proposal
aug <- simulate_tree(sim$tes, pars = c(0.8, -0.02, 0.15, 0.1, 0, 0), model = "nd",
                     n_trees = 50, method = "bdi")      # or method = "thinning"
length(aug$trees); aug$log_q

The simulator page runs the same forward process in the browser, for every model and link and any sampling fraction, drawing the complete, reconstructed and sampled trees as the clade grows.

Usage

The three routes on one simulated tree. The exponential link is used because the path is defined on it; the ladder and single fits run on any link.

library(emphasis); library(ape)

# a clade under diversity dependence: lambda = exp(b0 + bN N), mu = exp(g0)
sim <- simulate_tree(pars = c(log(0.6), -0.01, log(0.05), 0), max_t = 12,
                     model = "dd", link = "exponential")
phy <- sim$tes

# 1. one model
box <- auto_bounds(phy, model = "dd", link = "exponential", train_surv_gam = FALSE)
fit <- estimate_rates(phy, model = "dd", link = "exponential",
                      control = list(lower_bound = box$lower_bound,
                                     upper_bound = box$upper_bound,
                                     sample_size = 200))
fit                        # estimates, log-likelihood, effective sample, gap, stop reason, notes

# 2. the ladder: cr, dd, nd, ned, nk from nested starts, ranked by AIC
fits <- estimate_rates(phy, model = "all", link = "exponential",
                       control = list(sample_size = 200))
fits$comparison            # AIC, ESS, gap, swing; attr "ranking_at_risk"
fits$fits$nd$pars

# 3. the path: which covariates the data ask for first
cp <- covariate_path(phy, covariates = c("N", "D", "ED", "K"), criterion = "BIC")
cp$speciation$entry_order; cp$speciation$active

# the numbers the likelihood is built on
head(lineage_table(phy))

Called on a phylo, covariate_path() treats the tree as fully observed: the speciation path is exact and the extinction path, with no event to read, is returned empty. The augmented route, with the hidden lineages entering at their importance weights, runs on the draws of an E-step and is what the experiments measure.

Choosing the proposal

# the surrogate on a D model (default is thinning there), or thinning anywhere
estimate_rates(phy, model = "nd", link = "exponential",
               control = list(sampling = "bdi", sample_size = 200,
                              lower_bound = lb, upper_bound = ub))
estimate_rates(phy, model = "dd", control = list(sampling = "dynamic_fresh", ...))

# incomplete sampling: a fifth of the living species are missing
estimate_rates(phy, model = "dd", control = list(rho = 0.8, ...))

Function reference

Function	Purpose
`simulate_tree()`	Forward simulation under every model and link; given a tree, hidden lineages from either proposal (`method = "bdi"` / `"thinning"`)
`estimate_rates()`	Monte Carlo EM for one model, or the nested ladder with `model = "all"`; the proposal, draws, stopping rule and relaxation in `control`; the trace of every iteration in `details$mcem`
`compare_models()`	AIC table over fits with effective sample, gap, swing and the ranking-at-risk flag
`covariate_path()`	The dgLARS path over the covariates, for speciation and extinction; entry order and the set BIC or AIC chooses
`lineage_table()`	One row per lineage and segment: the covariates and the event mark the likelihood reads
`auto_bounds()`	The parameter box from the observed tree, and the optional survival surface
`diagnose_mcem()`	The trace of a fit: iterates, log-likelihood, effective sample, weights

Also exported, from earlier versions: emphasis_pipeline(), train_GAM(), predict_survival(), diagnose_cem() and diagnose_gam(), with estimate_rates(method = "cem" | "gam"). They are a likelihood surface fitted on a grid and a cross-entropy search, not part of the method described on this page, and nothing above depends on them.

emphasis — Francisco Richter (richtf@usi.ch), with Thijs Janzen and Hanno Hildenbrandt. MIT license. Source and issues on GitHub.

Symbol	Name	Definition at time \(t\)	Reads
\(N(t)\)	diversity	the number of lineages alive	the branching times
\(D_s(t)\)	centred pendant age	\(E_s(t) - M(t)\); sums to zero over the lineages alive	the topology: older or younger than its contemporaries
\(\mathrm{ED}_s(t)\)	evolutionary distinctiveness	fair proportion \(\sum_{b \succeq s} \ell_b / n_b\) over the branches from the crown to \(s\), each branch's length divided by the lineages alive below it; sums over the clade to Faith's phylogenetic diversity	the whole ancestry, weighted by how much of it \(s\) shares
\(\mathrm{ED}^c_s(t)\)	centred distinctiveness	\(\mathrm{ED}_s - \overline{\mathrm{ED}}\), the within-clade contrast alone	as \(\mathrm{ED}\), with the clade-level part removed
\(K_s(t)\)	ancestral splits	the number of speciation events on the path from the crown to \(s\); one more on both daughters at a split, hidden splits counted	whether a speciation-rich ancestry speciates more

Shortcut	Formula	Parameters	Hypothesis
`"cr"`	`~ 1`	\((\beta_0, \gamma_0)\)	constant rates, the null
`"dd"`	`~ N`	\((\beta_0, \beta_N, \gamma_0, \gamma_N)\)	diversity dependence
`"d"`	`~ D`	\((\beta_0, \beta_D, \gamma_0, \gamma_D)\)	age-dependent speciation
`"nd"`	`~ N + D`	6	age dependence over and above diversity
`"ed"`	`~ ED`	\((\beta_0, \beta_{ED}, \gamma_0, \gamma_{ED})\)	distinctiveness-dependent diversification
`"ned"`	`~ N + ED`	6	distinctiveness over and above diversity
`"edc"`, `"nedc"`	`~ EDc`, `~ N + EDc`	4, 6	the same on the centred covariate
`"k"`	`~ K`	\((\beta_0, \beta_K, \gamma_0, \gamma_K)\)	inherited diversification: a speciation-rich ancestry speciates more, or less
`"nk"`	`~ N + K`	6	the same over and above diversity

Link	Rate	Where it applies
Linear	\(\max(0,\ \beta_0 + \eta_{\text{cov}})\)	every model; a lineage whose rate is truncated to zero cannot split, so an observed split by such a lineage gives a draw density zero unless a hidden split preceded it
Exponential	\(\exp(\beta_0 + \eta_{\text{cov}})\)	every model; the only link on which the selection path is defined, since the point process is then log-linear
Gaussian	\(\beta_0 \, \exp\!\big(-(\eta_{\text{cov}} - 1)^2 / 2\big)\)	`cr`, `dd`, `d`, `nd`; not offered with \(\mathrm{ED}\) or \(K\)

Model	Surrogate proposal	Default
`cr`	exact on every link, \(\mathrm{ESS} = K\)	surrogate
`dd`	mean field on the linear and exponential links; the effective sample grows with the tree, 42 to 169 of 200 draws from 25 to 400 tips. Refused on the gaussian link, where the rate is not monotone in \(N\) and the iteration does not contract past its peak	surrogate
`d`, `nd`	mean field with exact attachment, linear and exponential links. On 12 trees at 100 tips it leaves an effective sample of 2–439 against 1.3–15 for thinning; in fits it lifts 3–6 to 13–33 in four cells of five, with recovery of \(\beta_D\) at 0.43–1.06 of its value and selection at 0.90 for an effect of −0.35 but 0.20 for −0.15	thinning; the surrogate on request
`ed`, `ned`, `edc`, `nedc`	mean field through the clade mean of fair proportion, linear and exponential links; 41 times the effective sample of thinning on `ned` trees at complete sampling, and ahead at every \(\rho\) down to 0.2. On `ed` it trails thinning below \(\rho \approx 0.8\)	surrogate
`k`, `nk`	refused: the mean field carries no count of ancestral splits	thinning