A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

An integer vector of size Ntips, specifying the state of each tip in the tree as an integer from 1 to Nstates, where Nstates is the possible number of states (see below).

Either NULL, or an integer specifying the number of possible states of the trait. If NULL, then it will be computed based on the maximum value encountered in tip_states

Logical, specifying whether empirical frequencies should be normalized to represent probabilities. If FALSE, then the raw occurrence counts are returned.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

A 2D integer (if probabilities==FALSE) or numeric (if probabilities==TRUE) matrix, listing the frequency or probability of each state for each node. This matrix will have size Nnodes x Nstates, where Nstates was either explicitly provided as an argument or inferred from tip_states. The rows in this matrix will be in the order in which nodes are indexed in the tree, i.e. the [n,s]-th entry will be the frequency or probability of the s-th state for the n-th node. Note that the name was chosen for compatibility with other ASR functions.

Integer vector of length Nnodes, listing the ancestral states with highest probability. In the case of ties, the first state with maximum probability is chosen. This vector is computed based on ancestral_likelihoods.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A numeric vector of size Ntips, specifying the known state of each tip in the tree.

Logical, specifying whether to weight tips and nodes by the inverse length of their incoming edge, as in the original method by Felsenstein (1985). If FALSE, edge lengths are treated as if they were 1.

Logical, specifying whether to also calculate standard errors and confidence intervals for the reconstructed states under a Brownian motion model, as described by Garland et al (1999).

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

A numeric vector of size Nnodes, listing the reconstructed state of each node. The entries in this vector will be in the order in which nodes are indexed in the tree.

Numeric vector of size Nnodes, listing the phylogenetically estimated standard error for the state in each node, under a Brownian motion model. The standard errors have the same units as the trait and depend both on the tree topology as well as the tip states. Calculated as described by Garland et al. (1999, page 377). Only included if include_CI==TRUE.

Numeric vector of size Nnodes, listing the radius (half width) of the 95% confidence interval of the state in each node. Confidence intervals have same units as the trait and depend both on the tree topology as well as the tip states. For each node, the confidence interval is calculated according to the Student's t-distribution with Npics degrees of freedom, where Npics is the number of internally calculated independent contrasts descending from the node [Garland et al, 1999]. Only included if include_CI==TRUE.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

An integer vector of size Ntips, specifying the state of each tip in the tree as an integer from 1 to Nstates, where Nstates is the possible number of states (see below).

Either NULL, or an integer specifying the number of possible states of the trait. If NULL, then Nstates will be computed based on the maximum value encountered in tip_states

Either "all_equal","sequential", "proportional", "exponential", or a quadratic non-negatively valued matrix of size Nstates x Nstates, specifying the transition costs between all possible states (which can include 0 as well as Inf). The [r,c]-th entry of the matrix is the cost of transitioning from state r to state c. The option "all_equal" specifies that all transitions are permitted and are equally costly. "sequential" means that only transitions between adjacent states are permitted from a node to its children, and all permitted transitions are equally costly. "proportional" means that all transitions are permitted, but the cost increases proportional to the distance between states. "exponential" means that all transitions are permitted, but the cost increases exponentially with the distance between states. The options "sequential" and "proportional" only make sense if states exhibit an order relation (as reflected in their integer representation).

Non-negative real-valued number. Optional exponent for weighting transition costs by the inverse length of edge lengths. If 0, edge lengths do not influence the ancestral state reconstruction (this is the conventional max-parsimony). If >0, then at each edge the transition costs are multiplied by 1/L^e, where L is the edge length and e is the edge exponent. This parameter is mostly experimental; modify at your own discretion.

Logical, indicating whether to weight each optimal state of a node by the number of optimal maximum-parsimony scenarios in which the node is in that state. If FALSE, then all optimal states of a node are weighted equally (i.e. are assigned equal probabilities).

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Boolean, indicating whether ASR was successful. If FALSE, the remaining returned elements may be undefined.

A 2D numeric matrix, listing the probability of each node being in each state. This matrix will have size Nnodes x Nstates, where Nstates was either explicitly provided as an argument or inferred from tip_states. The rows in this matrix will be in the order in which nodes are indexed in the tree, i.e. the [n,s]-th entry will be the probability of the s-th state for the n-th node. These probabilities are calculated based on scenario_counts (see below), assuming that every maximum parsimony scenario is equally likely. Note that the name was chosen for compatibility with other ASR functions.

A 2D numeric matrix of size Nnodes x Nstates, listing for each node and each state the number of maximum parsimony scenarios in which the node was in the specific state. If only a single maximum parsimony scenario exists for the whole tree, then the sum of entries in each row will be one.

Integer vector of length Nnodes, listing the ancestral states with highest likelihoods. In the case of ties, the first state with maximum likelihood is chosen. This vector is computed based on ancestral_likelihoods.

Real number, specifying the total transition cost across the tree for the most parsimonious scenario. In the classical case where transition_costs="all_equal", the total_cost equals the total number of state changes in the tree under the most parsimonious scenario. Under some constrained transition models (e.g., "sequential"), total_cost may sometimes be Inf, which basically means that the data violates the model.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

An integer vector of size Ntips, specifying the state of each tip in the tree in terms of an integer from 1 to Nstates, where Ntips is the number of tips and Nstates is the number of possible states (see below). Can also be NULL. If tip_states==NULL, then tip_priors must not be NULL (see below).

Either NULL, or an integer specifying the number of possible states of the trait. If Nstates==NULL, then it will be computed based on the maximum value encountered in tip_states or based on the number of columns in tip_priors (whichever is non-NULL).

A 2D numeric matrix of size Ntips x Nstates, where Nstates is the possible number of states for the character modelled. Hence, tip_priors[i,s] is the likelihood of the observed state of tip i, if the tip's true state was in state s. For example, if you know for certain that a tip is in state k, then set tip_priors[i,s]=1 for s=k and tip_priors[i,s]=0 for all other s.

Rate model to be used for fitting the transition rate matrix. Can be "ER" (all rates equal), "SYM" (transition rate i–>j is equal to transition rate j–>i), "ARD" (all rates can be different), "SUEDE" (only stepwise transitions i–>i+1 and i–>i-1 allowed, all 'up' transitions are equal, all 'down' transitions are equal) or "SRD" (only stepwise transitions i–>i+1 and i–>i-1 allowed, and each rate can be different). Can also be an index matrix that maps entries of the transition matrix to the corresponding independent rate parameter to be fitted. Diagonal entries should map to 0, since diagonal entries are not treated as independent rate parameters but are calculated from the remaining entries in the transition matrix. All other entries that map to 0 represent a transition rate of zero. The format of this index matrix is similar to the format used by the ace function in the ape package. rate_model is only relevant if transition_matrix==NULL.

Either a numeric quadratic matrix of size Nstates x Nstates containing fixed transition rates, or NULL. The [r,c]-th entry in this matrix should store the transition rate from state r to state c. Each row in this matrix must have sum zero. If NULL, then the transition rates will be estimated using maximum likelihood, based on the rate_model specified.

Prior probability distribution of the root's states, used to calculate the model's overall likelihood from the root's marginal ancestral state likelihoods. Can be "flat" (all states equal), "empirical" (empirical probability distribution of states across the tree's tips), "stationary" (stationary probability distribution of the transition matrix), "likelihoods" (use the root's state likelihoods as prior) or "max_likelihood" (put all weight onto the state with maximum likelihood). If "stationary" and transition_matrix==NULL, then a transition matrix is first fitted using a flat root prior, and then used to calculate the stationary distribution. root_prior can also be a non-negative numeric vector of size Nstates and with total sum equal to 1.

Include the marginal ancestral likelihoods for each node (conditional scaled state likelihoods) in the return values. Note that this may increase the computation time and memory needed, so you may set this to FALSE if you don't need marginal ancestral states.

Reroot tree at each node when computing marginal ancestral likelihoods, according to Yang et al. (1995). This is the default and recommended behavior, but leads to increased computation time. If FALSE, ancestral likelihoods at each node are computed solely based on the subtree descending from that node, without rerooting. Caution: Rerooting is strictly speaking only valid if the Mk model is time-reversible (for example, if the transition matrix is symmetric). Do not use the rerooting method if rate_model="ARD".

Number of trials (starting points) for fitting the transition matrix. Only relevant if transition_matrix=NULL. A higher number may reduce the risk of landing in a local non-global optimum of the likelihood function, but will increase computation time during fitting.

Either "optim" or "nlminb", specifying which optimization algorithm to use for maximum-likelihood estimation of the transition matrix. Only relevant if transition_matrix==NULL.

Maximum number of iterations (per fitting trial) allowed for optimizing the likelihood function.

Relative tolerance (stop criterion) for optimizing the likelihood function.

Logical, specifying whether to pre-calculate and store exponentials of the transition matrix during calculation of ancestral likelihoods. This may reduce computation time because each exponential is only calculated once, but requires more memory since all exponentials are stored.

Only relevant if include_ancestral_likelihoods==TRUE, otherwise exponentials are never stored.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Number of parallel threads to use for running multiple fitting trials simultaneously. This only makes sense if your computer has multiple cores/CPUs and if Ntrials>1, and is only relevant if transition_matrix==NULL. This option is ignored on Windows, because Windows does not support forking.

Logical, indicating whether ASR was successful. If FALSE, all other return values may be NULL.

Integer, specifying the number of modeled trait states.

A numeric quadratic matrix of size Nstates x Nstates, containing the transition rates of the Markov model. The [r,c]-th entry is the transition rate from state r to state c. Will be the same as the input transition_matrix, if the latter was not NULL.

Log-likelihood of the observed tip states under the fitted (or provided) Mk model. If transition_matrix was NULL in the input, then this will be the log-likelihood maximized during fitting.

Numeric, the Akaike Information Criterion for the fitted Mk model (if applicable), defined as 2k-2\log(L), where k is the number of independent fitted parameters and L is the maximized likelihood. If the transition matrix was provided as input, then no fitting was performed and hence AIC will be NULL.

Optional, only returned if include_ancestral_likelihoods was TRUE. A 2D numeric matrix, listing the likelihood of each state at each node (marginal ancestral likelihoods). This matrix will have size Nnodes x Nstates, where Nstates was either explicitly provided as an argument, or inferred from tip_states or tip_priors (whichever was non-NULL). The rows in this matrix will be in the order in which nodes are indexed in the tree, i.e. the [n,s]-th entry will be the likelihood of the s-th state at the n-th node. For example, likelihoods[1,3] will store the likelihood of observing the tree's tip states (if reroot=TRUE) or the descending subtree's tip states (if reroot=FALSE), if the first node was in state 3. Note that likelihoods are rescaled (normalized) to sum to 1 for convenience and numerical stability. The marginal likelihoods at a node should not, however, be interpreted as a probability distribution among states.

Integer vector of length Nnodes, listing the ancestral states with highest likelihoods. In the case of ties, the first state with maximum likelihood is chosen. This vector is computed based on ancestral_likelihoods and is only included if include_ancestral_likelihoods=TRUE.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A numeric vector of size Ntips, specifying the known state of each tip in the tree.

Logical, specifying whether to weight transition costs by the inverted edge lengths. This corresponds to the "weighted squared-change parsimony" reconstruction by Maddison (1991) for a Brownian motion model of trait evolution.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

A numeric vector of size Nnodes, listing the reconstructed state of each node. The entries in this vector will be in the order in which nodes are indexed in the tree.

The total sum of squared changes, minimized by the (optionally weighted) squared-change parsimony algorithm. This is equation 7 in (Maddison, 1991).

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A numeric vector of size Ntips, specifying the known state of each tip in the tree.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Logical, indicating whether ASR was sucessful. If all input data are valid then this will always be TRUE, but it is provided for consistency with other ASR functions.

A numeric vector of size Nnodes, listing the reconstructed state (=average over descending tips) for each node. The entries in this vector will be in the order in which nodes are indexed in the tree.

A numeric vector of size Nnodes, listing the standard deviations corresponding to ancestral_stds.

A numeric vector of size Nnodes, listing the number of (descending) tips used to reconstruct the state of each node.

A rooted tree of class "phylo", where edge lengths represent time intervals (or similar).

Integer, number of equidistant age bins at which to calculate densities.

Numeric, minimum age to consider. If NULL, it will be set to the minimum possible.

Numeric, maximum age to consider. If NULL, it will be set to the maximum possible.

Logical, whether to normalize densities by the local number of lineages (in addition to dividing by the age interval). Hence, tip (or node) densities will represent number of tips (or nodes) per time unit per lineage.

Logical, specifying whether the input tree is guaranteed to be ultrametric, even in the presence of some numerical inaccuracies causing some tips not have exactly the same distance from the root.

Integer, indicating the number of discrete age bins.

Numeric vector of size Nbins, listing the centres of the age bins.

Numeric vector of size Nbins, listing the tip densities in the corresponding age bins.

Numeric vector of size Nbins, listing the node densities in the corresponding age bins.

A rooted tree of class "phylo".

Logical, indicating whether the root node should always be kept (i.e., even if it only has a single child).

Logical, indicating whether all edges should be assumed to have length 1. If TRUE, the outcome is the same as if the tree had no edges.

A new tree of class "phylo", containing only bifurcations (and multifurcations, if these existed in the input tree). The number of nodes in this tree, Nnodes_new, may be lower than of the input tree.

Integer vector of length Nnodes_new, mapping node indices in the new tree to node indices in the old tree.

Integer. Number of nodes (monofurcations) removed from the tree.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

Numeric, specifying the phylogenetic resolution at which to collapse the tree. This is the maximum distance a descending tip can have from a node, such that the node is collapsed into a new tip. If set to 0 (default), then only nodes whose descending tips are identical to the node will be collapsed.

Logical. Instead of considering edge lengths, consider edge counts as phylogenetic distance between nodes and tips. This is the same as if all edges had length equal to 1.

Logical, indicating whether collapsed nodes should be turned into tips at the same location (thus potentially shortening the tree). If FALSE, then the incoming edge of each collapsed node is extended by some length L, where L is the distance of the node to its farthest descending tip (thus maintaining the height of the tree).

Logical, indicating whether collapsed nodes should be renamed using a representative tip name (the farthest descending tip). See details below.

Character, specifying the criterion to use for collapsing (i.e. how to interpret resolution). 'max_tip_depth': Collapse nodes based on their maximum distance to any descending tip. 'sum_tip_paths': Collapse nodes based on the sum of descending edges (each edge counted once). 'max_tip_pair_dist': Collapse nodes based on the maximum distance between any pair of descending tips.

A new rooted tree of class "phylo", containing the collapsed tree.

Numeric, indicating the phylogenetic distance between the old and the new root. Will always be non-negative.

Integer vector, listing indices of collapsed nodes in the original tree (subset of 1,..,Nnodes).

Integer vector of the same length as collapsed_nodes, listing indices of the farthest tips for each collapsed node. Hence, farthest_tips[n] will be the index of a tip in the original tree that descended from node collapsed_nodes[n] and had the greatest distance from that node among all descending tips.

Integer vector of length equal to the number of tips+nodes in the collapsed tree, with values in 1,..,Ntips+Nnodes, mapping tip/node indices of the collapsed tree to tip/node indices in the original tree.

Integer vector of length equal to the number of edges in the collapsed tree, with values in 1,..,Nedges, mapping edge indices of the collapsed tree to edge indices in the original tree.

Integer vector of length equal to the number of tips+nodes in the original tree, mapping tip/node indices of the original tree to tip/node indices in the collapsed tree. Non-mapped tips/nodes (i.e., missing from the collapsed tree) will be represented by zeros.

A rooted tree object of class "phylo". Usually this tree will be time-calibrated (i.e. edge lengths represent time intervals).

A rooted tree object of class "phylo".

A table mapping a subset of target tips to a subset of reference tips, as described by Eastman et al (2013). Multiple target tips may map to the same reference tip, but not vice versa (i.e. every target tip can appear at most once in the mapping). In general, a tip mapped to in the reference tree is assumed to represent a monophyletic group of tips in the target tree, although this assumption may be violated in practice (Eastman et al. 2013).

The mapping must be in one of the following formats:

Option 1: A 2D integer array of size NM x 2 (with NM being the number of mapped target tips), listing target tip indices mapped to reference tip indices (mapping[m,1] (target tip) –> mapping[m,2] (reference tip)).

Option 2: A 2D character array of size NM x 2, listing target tip labels mapped to reference tip labels.

Option 3: A data frame of size NM x 1, whose row names are target tip labels and whose entries are either integers (reference tip indices) or characters (reference tip labels). This is the format used by geiger::congruify.phylo (v.206).

Option 4: A vector of size NM, whose names are target tip labels and whose entries are either integers (reference tip indices) or characters (reference tip labels).

Integer vector of length NC (where NC is the number of concordant node pairs found) and with values in 1,..,RNnodes, listing indices of reference nodes that could be matched with (i.e. were concordant to) a target node. Entries in Rnodes will correspond to entries in Tnodes and ages.

Integer vector of length NC and with values in 1,..,TNnodes, listing indices of target nodes that could be matched with (i.e. were concordant to) a reference node. Entries in Tnodes will correspond to entries in Rnodes and ages.

Numeric vector of length NC, listing divergence times (ages) of the reference nodes listed in Rnodes. These ages can be used as fixed anchors for time-calibrating the target tree using a separate program (such as PATHd8).

Numeric vector, listing discrete ages (time before present) on which the various model variables (e.g., \lambda, \mu etc) are specified. Listed ages must be strictly increasing, and must include the present-day (i.e. age 0).

Numeric vector, of the same size as age_grid, specifying the pulled speciation rate (PSR) (in units 1/time) at the ages listed in age_grid. The PSR is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case PSR is assumed to be time-independent.

Numeric vector, of the same size as age_grid, specifying the pulled diversification rate (PDR) (in units 1/time) at the ages listed in age_grid. The PDR is assumed to vary polynomially between grid points (see argument splines_degree). The PDR of a HBDS model is defined as PDR=\lambda-\mu-\psi+(1/\lambda)d\lambda/dt (where t denotes age). Can also be a single number, in which case PDR is assumed to be time-independent.

Numeric vector, of the same size as age_grid, specifying the product of speciation rate and sampling rate (\lambda\psi, in units 1/time^2) at the ages listed in age_grid. \lambda\psi is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \lambda\psi is assumed to be time-independent.

Numeric vector, of the same size as age_grid, specifying the speciation rate (\lambda, in units 1/time) at the ages listed in age_grid. The speciation rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \lambda is assumed to be time-independent. By providing \lambda, one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.

Numeric vector, of the same size as age_grid, specifying the extinction rate (\mu, in units 1/time) at the ages listed in age_grid. The extinction rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \mu is assumed to be time-independent. In an epidemiological context, \mu typically corresponds to the recovery rate plus the death rate of infected individuals. By providing \mu (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.

Numeric vector, of the same size as age_grid, specifying the (Poissonian) sampling rate (\psi, in units 1/time) at the ages listed in age_grid. The sampling rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \psi is assumed to be time-independent. By providing \psi, one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.

Numeric vector, of the same size as age_grid, specifying the effective reproduction ratio (R_e, unitless) at the ages listed in age_grid. The R_e is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case R_e is assumed to be time-independent. By providing R_e (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.

Numeric vector, of the same size as age_grid, specifying the removal rate (\mu+\psi, in units 1/time) at the ages listed in age_grid. IN an epidemiological context this is also known as "become uninfectious" rate. The removal rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case the removal rate is assumed to be time-independent. By providing \mu+\psi (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.

Numeric, specifying the speciation rate at the present-day (i.e., at age 0). Must be provided if and only if one of mu, Reff or removal_rate is provided.

Optional numeric vector, listing the ages of concentrated sampling attempts, in ascending order. Concentrated sampling is assumed to occur in addition to any continuous (Poissonian) sampling specified by psi.

Optional numeric vector of the same size as CSA_ages, listing pulled sampling probabilities at each concentrated sampling attempt (CSA). Note that in contrast to the sampling rates psi, the CSA_pulled_probs are interpreted as probabilities and must thus be between 0 and 1. CSA_pulled_probs must be provided if and only if CSA_ages is provided.

Optional numeric vector of the same size as CSA_ages, specifying the pulled sampling rate (PSR) during each concentrated sampling attempt. While in principle the PSR is already provided by the argument PSR, the PSR may be non-continuous at CSAs, which makes a representation as piecewise polynomial function difficult; hence, you must explicitly provide the correct PSR at each CSA. CSA_PSRs must be provided if and only if CSA_ages is provided.

Integer, either 0,1,2 or 3, specifying the polynomial degree of the provided time-dependent variables between grid points in age_grid. For example, if splines_degree==1, then the provided PDR, PSR and so on are interpreted as piecewise-linear curves; if splines_degree==2 they are interpreted as quadratic splines; if splines_degree==3 they are interpreted as cubic splines. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on.

Positive unitless number, specifying the default relative time step for internally used ordinary differential equation solvers. Typical values are 0.01-0.001.

Positive unitless number, specifying the relative difference between subsequent simulated and interpolated values, in internally used ODE solvers. Typical values are 1e-2 to 1e-5. A smaller ODE_relative_dy increases interpolation accuracy, but also increases memory requirements and adds runtime.

Logical, indicating whether the calculation was successful. If FALSE, then the returned list includes an additional 'error' element (character) providing a description of the error; all other return variables may be undefined.

Logical, indicating whether the returned model appears to be biologically valid (for example, does not have negative \lambda, \mu or \psi). In principle, a congruence class may include biologically invalid models, which might be returned depending on the input to congruent_hbds_model. Note that only biologically valid models can be subsequently simulated using simulate_deterministic_hbds.

Numeric vector of size NG, specifying the discrete ages (time before present) on which all returned time-curves are specified. Will always be equal to age_grid.

Numeric vector of size NG, listing the speciation rates \lambda of the returned model at the ages given in ages[].

Numeric vector of size NG, listing the extinction rates \mu of the returned model at the ages given in ages[].

Numeric vector of size NG, listing the (Poissonian) sampling rates \psi of the returned model at the ages given in ages[].

Numeric vector of size NG, listing the product \lambda\psi at the ages given in ages[].

Numeric vector of size NG, listing the effective reproduction ratio R_e of the returned model at the ages given in ages[].

Numeric vector of size NG, listing the removal rate (\mu+\psi, aka. "become uninfectious" rate) of the returned model at the ages given in ages[].

Numeric vector of size NG, listing the probability that a past lineage extant during ages[] will be missing from a tree generated by the model.

Numeric vector of the same size as CSA_ages, listing the sampling probabilities at each of the CSAs.

Numeric vector of the same size as CSA_ages, listing the probability that a past lineage extant during each of CSA_ages[] will be missing from a tree generated by the model.

A rooted tree of class "phylo".

Optional, character vector of length Ntips, listing taxonomic paths for the tips. If NULL, then tip labels in the tree are assumed to be tip taxonomies.

Character, the delimiter between taxonomic levels (e.g., ";" for SILVA and Greengenes taxonomies).

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A numeric vector of size Ntips indicating absence (value <=0) or presence (value >0) of a particular trait at each tip of the tree. Note that tip_states[i] (where i is an integer index) must correspond to the i-th tip in the tree.

Minimum fraction of tips in a clade exhibiting the trait, for the clade to be considered "positive" in the trait. In the original paper by Martiny et al (2013), this was 0.9.

Logical, specifying whether to include singletons in the statistics (tips positive in the trait, but not part of a larger positive clade). The phylogenetic depth of singletons is taken to be half the length of their incoming edge, as proposed by Martiny et al (2013). If FALSE, singletons are ignored.

Numeric, specifying the phylogenetic resolution at which to resolve singletons. Any clade found to be positive in a trait will be considered a singleton if the distance of the clade's root to all descending tips is below this threshold.

Whether to weight positive clades by their number of positive tips. If FALSE, each positive clades is weighted equally, as proposed by Martiny et al (2013).

Number of random permutations for estimating the statistical significance of the mean trait depth. If zero (default), the statistical significance is not calculated.

Mean phylogenetic depth of clades that are positive in the trait.

Variance of phylogenetic depths of clades that are positive in the trait.

Minimum phylogenetic depth of clades that are positive in the trait.

Maximum phylogenetic depth of clades that are positive in the trait.

Number of clades that are positive in the trait.

Statistical significance (P-value) of mean_depth, under a null model of random trait presences/absences (see details above). This is the probability that, under the null model, the mean_depth would be at least as high as observed in the data.

Mean random mean_depth, under a null model of random trait presences/absences (see details above).

Integer vector, listing indices of tips and nodes (from 1 to Ntips+Nnodes) that were found to be positive in the trait and counted towards the statistic.

Integer vector of size Ntips+Nnodes, listing the number of descending tips per clade (tip or node) that were positive in the trait.

Numeric vector of size Ntips+Nnodes, listing the mean phylogenetic depth of each clade (tip or node), i.e. the average distance to all its descending tips.

A rooted tree of class "phylo".

Numeric vector of size Ntips, specifying the latitudes (decimal degrees) of the tree's tips. By convention, positive latitudes correspond to the northern hemisphere. Note that tip_latitudes[i] must correspond to the i-th tip in the tree, i.e. as listed in tree$tip.label.

Numeric vector of size Ntips, specifying the longitudes (decimal degrees) of the tree's tips. By convention, positive longitudes correspond to the eastern hemisphere.

Character, one of "pearson", "spearman" or "kendall".

Numeric, maximum phylodistance between tips to consider, in the same units as the tree's edge lengths. If Inf, all tip pairs will be considered.

Integer, number of random permutations to consider for estimating the statistical significance (P value). If 0, the significance will not be computed. A larger number improves accuracy but at the cost of increased computing time.

Character, one of "two_sided", "right" or "left", specifying which part of the null model's distribution to use as P-value.

Optional numeric, radius to assume for the sphere. If 1, then all geodistances are measured in multiples of the sphere radius. This does not affect the correlation or P-value, but it affects the returned geodistances. Note that Earth's average radius is about 6371 km.

Numeric between -1 and 1, the correlation between phylodistances and geodistances.

Integer, the number of tip pairs considered.

Numeric between 0 and 1, estimated statistical significance of the correlation. Only returned if Npermutations>0.

Numeric between -1 and 1, the mean correlation obtained across all random permutations. Only returned if Npermutations>0.

Numeric vector of length Npairs, listing the pairwise phylodistances between tips, used to compute the correlation.

Numeric vector of length Npairs, listing the pairwise geodistances between tips, used to compute the correlation.

A rooted tree of class "phylo", where edge lengths represent time intervals (or similar).

Integer, number of equidistant time points at which to count lineages. Can also be NULL, in which case times must be provided.

Minimum time (distance from root) to consider. If NULL, this will be set to the minimum possible (i.e. 0). Only relevant if times==NULL.

Maximum time (distance from root) to consider. If NULL, this will be set to the maximum possible. Only relevant if times==NULL.

Integer vector, listing time points (in ascending order) at which to count lineages. Can also be NULL, in which case Ntimes must be provided.

Logical, specifying whether the slope and the relative slope of the returned clades-per-time-point curve should also be returned.

Logical, specifying whether the input tree is guaranteed to be ultrametric, even in the presence of some numerical inaccuracies causing some tips not have exactly the same distance from the root. If you know the tree is ultrametric, then this option helps the function choose a better time grid for the LTT.

Integer, specifying the "degree" of the LTT curve: LTT(t) will be the number of lineages in the tree at time t that have at least n descending tips in the tree. Typically order=1, which corresponds to the classical LTT curve.

Logical, specifying whether the automatically generated time grid should be regular (equal distances between grid points). This option only matters if times==NULL. If regular_grid==FALSE and times==NULL, the time grid will be irregular, with grid point density being roughly proportional to the square root of the number of lineages at any particular time (i.e., the grid becomes finer towards the tips).

Integer, indicating the number of returned time points. Equal to the provided Ntimes if applicable.

Numeric vector of size Ntimes, listing the time points at which the LTT was calculated. If times was provided as an argument to the function, then this will be the same as provided, otherwise times will be listed in ascending order.

Numeric vector of size Ntimes, listing the ages (time before the youngest tip) corresponding to the returned times[].

Integer vector of size Ntimes, listing the number of lineages represented in the tree at each time point that have at least degree descending tips, i.e. the LTT curve.

Numeric vector of size Ntimes, listing the slopes (finite-difference approximation of 1st derivative) of the LTT curve.

Numeric vector of size Ntimes, listing the relative slopes of the LTT curve, i.e. slopes divided by a sliding-window average of lineages.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

An integer vector or character vector of size Npairs, specifying the first of the two members of each pair of tips/nodes. If an integer vector, it must list indices of tips (from 1 to Ntips) and/or nodes (from Ntips+1 to Ntips+Nnodes). If a character vector, it must list tip and/or node names.

An integer vector or character vector of size Npairs, specifying the second of the two members of each pair of tips/nodes. If an integer vector, it must list indices of tips (from 1 to Ntips) and/or nodes (from Ntips+1 to Ntips+Nnodes). If a character vector, it must list tip and/or node names.

Integer vector of length Ntips+Nnodes, listing the discrete state of each tip and node in the tree. The order of entries must match the order of tips and nodes in the tree; this requirement is only verified if states has names and check_input==TRUE.

Logical, whether to perform basic validations of the input data. If you know for certain that your input is valid, you can set this to FALSE to reduce computation time.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

Integer, ranging between 1 and Nnodes. Index of the node to be used as dating anchor. If NULL, the tree's root is used as anchor.

Positive numeric. Age of the anchor node.

Logical, indicating whether the dating was successful. If FALSE, all other return values (except for error) may be undefined.

A new rooted tree of class "phylo", representing the dated tree.

Numeric vector of size Nnodes, listing the RED of each node in the input tree.

Character, listing any error message if success==FALSE.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A vector of size Ntips specifying the state at each tip. Note that tip_states[i] (where i is an integer index) must correspond to the i-th tip in the tree. This vector may be of any base data type, although character or integer are the most typical types.

Minimum fraction of tips in a clade that must have the dominant state, for the clade to be considered "uniform" in the trait.

Logical, specifying whether to consider singleton clades in the statistics (e.g., tips not part of a larger uniform clade). The phylogenetic depth of singletons is taken to be half the length of their incoming edge, as proposed by Martiny et al (2013). If FALSE, singletons are ignored. If you suspect a high risk of false positives in the detection of a trait, it may be worth setting count_singletons to FALSE to avoid skewing the distribution of conservation depths towards shallower depths due to false positives.

Numeric, specifying the phylogenetic resolution at which to resolve singletons. A clade will be considered a singleton if the distance of the clade's root to all descending tips is below this threshold.

Whether to weight uniform clades by their number of tips in the dominant state. If FALSE, each uniform clades is weighted equally.

Number of random permutations for estimating the statistical significance of the mean trait depth. If zero (default), the statistical significance is not calculated.

Vector of the same type as tip_states and of length Nstates, listing the unique possible states of the trait.

Numeric, specifying the mean phylogenetic depth of the trait, i.e., the mean depth of considered maximally uniform clades.

Numeric, specifying the variance of phylogenetic depths of considered maximally uniform clades.

Numeric, specifying the minimum phylogenetic depth of considered maximally uniform clades.

Numeric, specifying the maximum phylogenetic depth of considered maximally uniform clades.

Number of considered maximal uniform clades.

Numeric vector of size Nstates. Mean depth of considered maximally uniform clades, separately for each state and in the same order as unique_states. Hence, mean_depth_per_state[s] lists the mean depth of considered maximally uniform clades whose dominant state is unique_states[s].

Numeric vector of size Nstates. Variance of depths of considered maximally uniform clades, separately for each state and in the same order as unique_states

Numeric vector of size Nstates. Minimum phylogenetic depth of considered maximally uniform clades, separately for each state and in the same order as unique_states

Numeric vector of size Nstates. Maximum phylogenetic depth of considered maximally uniform clades, separately for each state and in the same order as unique_states

Integer vector of size Nstates. Number of considered maximally uniform clades, seperately for each state and in the same order as unique_states

Statistical significance (P-value) of mean_depth, under a null model of random tip states (see details above). This is the probability that, under the null model, the mean_depth would be at least as high as observed in the data.

Mean random mean_depth, under the null model of random tip states (see details above).

Numeric vector, listing x-values in ascending order.

Numeric vector of the same length as Xgrid, listing the values of Y on Xgrid.

Integer, either 0, 1, 2 or 3, specifying the polynomial degree of the spline curve Y between grid points. For example, 0 means Y is piecewise constant, 1 means Y is piecewise linear and so on.

Numeric vector, listing arbitrary X values on which to evaluate Y.

Character, specifying how to extrapolate Y beyond Xgrid if needed. Available options are "const" (i.e. use the value of Y on the nearest Xgrid point) or "splines" (i.e. use the polynomial coefficients from the nearest grid point).

Integer, specifying which derivative to return. To return the spline's value, set derivative=0. Currently only the options 0,1,2 are supported.

Character, the path to the input jplace file.

Logical, only keep the best placement of each query, i.e., the placement with maximum likelihood.

Positive integer, maximum number of sequence names to keep from each query. Only relevant if queries in the jplace file include multiple sequence names (these typically represent identical sequences). If greater than 1, and a query includes multiple sequence names, then each of these sequence names will be added as a tip to the tree.

Object of class "phylo", the extended tree constructed by adding the placements on the reference tree.

Integer vector, specifying which tips in the returned tree correspond to placements.

Object of class "phylo", the original reference tree loaded from the jplace file. This will be a subtree of tree.

Numeric, the diffusivity (aka. diffusion coefficient) of the SBM. The units of the diffusivity must be consistent with the units used for specifying the radius and time intervals (deltas); for example, if radius is in km and deltas are in years, then diffusivity must be specified in km^2/year.

Positive numeric, the radius of the sphere.

Numeric vector, listing time intervals for which to compute the expected geodesic distances.

A real-valued quadratic matrix of size N x N.

Vector of real-valued scalar scaling factors T, for each of which the exponential \exp(T\cdot A) should be calculated.

Maximum allowed absolute error for the returned approximations. A smaller allowed error implies a greater computational cost as more matrix polynomials need to be included (see below). The returned approximations may have a greater error if the parameter max_polynomials is set too low.

Minimum number of polynomials to include in the approximations (see equation below), even if max_absolute_error may be satisfied with fewer polynomials. If you don't know how to choose this, in most cases the default is fine. Note that regardless of min_polynomials and max_absolute_error, the number of polynomials used will not exceed max_polynomials.

Maximum allowed number of polynomials to include in the approximations (see equation below). Meant to provide a safety limit for the amount of memory and the computation time required. For example, a value of 1000 means that up to 1000 matrices (powers of A) of size N x N may be computed and stored temporarily in memory. Note that if max_polynomials is too low, the requested accuracy (via max_absolute_error) may not be achieved. That said, for large trees more memory may be required to store the actual result rather than the intermediate polynomials, i.e. for purposes of saving RAM it doesn't make much sense to choose max_polynomials much smaller than the length of scalings.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

Numeric, specifying the phylogenetic distance from the root to which tips are to be extended. If NULL or negative, then it is set to the maximum distance of any tip from the root.

A new rooted tree of class "phylo", representing the extended tree.

Numeric. The largest elongation added to a terminal edge.

A rooted tree of class "phylo".

Strictly positive integer, specifying the maximum number of splits to descend from the root. A larger value generally implies that more tips will be extracted, representing a larger number of splits.

Boolean, specifying whether to only return the subset of extracted tips, rather than the subtree spanned by those tips. If FALSE, a subtree is returned in addition to the tips (this comes at a computational cost).

Integer vector with values from 1 to Ntips, listing the indices of the extracted tips.

A new tree of class "phylo", the subtree spanned by the extracted tips. Only included if only_tips==FALSE.

A rooted tree of class "phylo".

Integer, specifying the number of constraints extracted.

2D character matrix of size Ntips x Nconstraints, with values '0', '1' or '-', specifying which side ("left" or "right") of a constraint (node) each tip is found on.

Integer vector of size Nconstraints, with values in 1,..,Nnodes, specifying the original node index used to define each constraint.

A rooted tree of class "phylo".

Either a character, specifying the name of the focal tip, or an integer between 1 and Ntips, specifying the focal tip's index.

Integer >=1, specifying how many splits backward (towards the root) to explore. A larger value of Nbackward will generally lead to a larger extracted subtree (i.e., including deeper splits).

Non-negative integer, specifying how many splits forward (towards the tips) to explore. A larger value of Nforward will generally lead to a larger extracted subtree (i.e., including more representative tips from each sister branch).

Optional integer or character list, specifying indices or names of tips to force-include in any case.

Logical, whether to actually extract the subtree, rather than just returning the inferred neighbor tips.

Integer vector with values in 1,..,Ntips, specifying which tips were found to be neighbors of the focal tip.

A new tree of class "phylo", containing a subset of the tips and nodes in the vicinity of the focal tip. Only returned if include_subtree=TRUE.

Integer vector of length Ntips_extracted (=number of tips in the extracted subtree) with values in 1,..,Ntips, mapping tip indices of the extracted subtree to tip indices in the original tree. In particular, tree$tip.label[new2old_tip] will be equal to subtree$tip.label. Only returned if include_subtree=TRUE.

A rooted tree of class "phylo".

Either a character, specifying the name of the focal tip, or an integer between 1 and Ntips, specifying the focal tip's index.

Non-negative numeric, specifying the patristic distance radius to consider around the focal tip.

Logical, whether to actually extract the subtree, rather than just returning the inferred within-radius tips.

Integer vector with values in 1,..,Ntips, specifying which tips were found to be within the specified radius of the focal tip.

A new tree of class "phylo", containing only the tips within the specified radius from the focal tip, and the nodes & edges connecting those tips to the root. Only returned if include_subtree=TRUE.

Integer vector of length Ntips_extracted (=number of tips in the extracted subtree, i.e., within the specified distance radius from the focal tip) with values in 1,..,Ntips, mapping tip indices of the extracted subtree to tip indices in the original tree. In particular, tree$tip.label[new2old_tip] will be equal to subtree$tip.label. Only returned if include_subtree=TRUE.

A rooted tree of class "phylo". While the tree must be rooted for technical reasons, the outcome does not actually depend on the rooting.

Logical, specifying whether to count phylogenetic distance in terms of edge counts instead of cumulative edge lengths. This is the same as setting all edge lengths to 1.

An integer between 1 and Ntips, specifying the first of the two most distant tips.

An integer between 1 and Ntips, specifying the second of the two most distant tips.

Numeric, specifying the phylogenetic (patristic) distance between the farthest_tip1 and farthest_tip2.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A logical indicating whether the farthest tip to a node or tip should be chosen from its descending tips only. If FALSE, then the whole set of possible target tips is considered.

Optional integer vector or character vector listing the subset of target tips to restrict the search to. If an integer vector, this should list tip indices (values in 1,..,Ntips). If a character vector, it should list tip names (in this case tree$tip.label must exist). If target_tips is NULL, then all tips of the tree are considered as target tips.

Logical, specifying whether to count phylogenetic distance in terms of edge counts instead of cumulative edge lengths. This is the same as setting all edge lengths to 1.

Logical, whether to perform basic validations of the input data. If you know for certain that your input is valid, you can set this to FALSE to reduce computation time.

An integer vector of size Ntips, listing the farthest target tip for each tip in the tree. Hence, farthest_tip_per_tip[i] is the index of the farthest tip (from the set of target tips), with respect to tip i (where i=1,..,Ntips). Some values may appear multiple times in this vector, if multiple tips share the same farthest target tip.

An integer vector of size Nnodes, listing the index of the farthest target tip for each node in the tree. Hence, farthest_tip_per_node[i] is the index of the farthest tip (from the set of target tips), with respect to node i (where i=1,..,Nnodes). Some values may appear multiple times in this vector, if multiple nodes share the same farthest target tip.

Integer vector of size Ntips. Phylogenetic ("patristic") distance of each tip in the tree to its farthest target tip. If only_descending_tips was set to TRUE, then farthest_distance_per_tip[i] will be set to infinity for any tip i that is not a target tip.

Integer vector of size Nnodes. Phylogenetic ("patristic") distance of each node in the tree to its farthest target tip. If only_descending_tips was set to TRUE, then farthest_distance_per_node[i] will be set to infinity for any node i that has no descending target tips.

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

A logical indicating whether the nearest tip to a node or tip should be chosen from its descending tips only. If FALSE, then the whole set of possible target tips is considered.

Optional integer vector or character vector listing the subset of target tips to restrict the search to. If an integer vector, this should list tip indices (values in 1,..,Ntips). If a character vector, it should list tip names (in this case tree$tip.label must exist). If target_tips is NULL, then all tips of the tree are considered as target tips.

Logical, specifying whether to count phylogenetic distance in terms of edge counts instead of cumulative edge lengths. This is the same as setting all edge lengths to 1.

Logical, whether to perform basic validations of the input data. If you know for certain that your input is valid, you can set this to FALSE to reduce computation time.

An integer vector of size Ntips, listing the nearest target tip for each tip in the tree. Hence, nearest_tip_per_tip[i] is the index of the nearest tip (from the set of target tips), with respect to tip i (where i=1,..,Ntips). Some values may appear multiple times in this vector, if multiple tips share the same nearest target tip.

An integer vector of size Nnodes, listing the index of the nearest target tip for each node in the tree. Hence, nearest_tip_per_node[i] is the index of the nearest tip (from the set of target tips), with respect to node i (where i=1,..,Nnodes). Some values may appear multiple times in this vector, if multiple nodes share the same nearest target tip.

Numeric vector of size Ntips. Phylogenetic ("patristic") distance of each tip in the tree to its nearest target tip. If only_descending_tips was set to TRUE, then nearest_distance_per_tip[i] will be set to infinity for any tip i that is not a target tip.

Numeric vector of size Nnodes. Phylogenetic ("patristic") distance of each node in the tree to its nearest target tip. If only_descending_tips was set to TRUE, then nearest_distance_per_node[i] will be set to infinity for any node i that has no descending target tips.

A tree of class "phylo". If the tree is not rooted, the function will return NA.

A tree object of class "phylo". Can be unrooted or rooted.

Character or integer vector, specifying the names or indices, respectively, of the target tips that should be turned monophyletic. If an integer vector, its elements must be between 1 and Ntips. If a character vector, its elements must be elements in tree$tip.label.

Logical, specifying whether the new root should become the MRCA of the target tips. If FALSE, the new root is chosen such that the MRCA of the target tips is the child of the new root.

Logical, specifying whether the input tree can be assumed to be rooted. If you are sure that the input tree is rooted, set this to TRUE for computational efficiency, otherwise to be on the safe side set this to FALSE.

Either a single rooted tree or a list of rooted trees, of class "phylo", corresponding to the first data set on which a BM model is to be fitted. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

Numeric state of each trait at each tip in each tree in the first data set. If trees1 is a single tree, then tip_states1 must either be a numeric vector of size Ntips or a 2D numeric matrix of size Ntips x Ntraits, listing the trait states for each tip in the tree. If trees1 is a list of Ntrees trees, then tip_states1 must be a list of length Ntrees, each element of which lists the trait states for the corresponding tree (as a vector or 2D matrix, similarly to the single-tree case).

Either a single rooted tree or a list of rooted trees, of class "phylo", corresponding to the second data set on which a BM model is to be fitted. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

Numeric state of each trait at each tip in each tree in the second data set, similarly to tip_states1.

Integer, specifying the number of parametric bootstraps to perform for calculating the confidence intervals of BM diffusivities fitted to each data set. If <=0, no bootstrapping is performed.

Integer, specifying the number of simulations to perform for assessing the statistical significance of the log-transformed difference between the diffusivities fitted to the two data sets, i.e. of |\log(D_1)-\log(D_2)|. Set to 0 to not calculate the statistical significance. See below for additional details.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Logical, specifying whether to print progress report messages to the screen.

Character, specifying a prefix to include in front of progress report messages on each line. Only relevant if verbose==TRUE.

Logical, indicating whether the fitting was successful for both data sets. If FALSE, then an additional return variable, error, will contain a description of the error; in that case all other return variables may be undefined.

A named list containing the fitting results for the first data set, in the same format as returned by fit_bm_model. In particular, the diffusivity fitted to the first data set will be stored in fit1$diffusivity.

A named list containing the fitting results for the second data set, in the same format as returned by fit_bm_model. In particular, the diffusivity fitted to the second data set will be stored in fit2$diffusivity.

The absolute difference between the log-transformed diffusivities, i.e. |\log(D_1)-\log(D_2)|. In the multivariate case, this will be a matrix of size Ntraits x Ntraits.

Numeric, statistical significance of the observed log-difference under the null hypothesis that the two data sets were generated by a common BM model. Only returned if Nsignificance>0.

A named list containing the fitting results for the two data sets combined, in the same format as returned by fit_bm_model. The common diffusivity, fit_common$diffusivity is used for the random simulations when assessing the statistical significance of the log-difference of the separately fitted diffusivities. Only returned if Nsignificance>0.

Either a single rooted tree or a list of rooted trees, of class "phylo", corresponding to the first data set on which an SBM model is to be fitted. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

Numeric vector listing the latitude (in decimal degrees) of each tip in each tree in the first data set. If trees1 is a single tree, then tip_latitudes1 must be a numeric vector of size Ntips, listing the latitudes for each tip in the tree. If trees1 is a list of Ntrees trees, then tip_latitudes1 must be a list of length Ntrees, each element of which lists the latitudes for the corresponding tree (as a vector, similarly to the single-tree case).

Similar to tip_latitudes1, but listing longitudes (in decimal degrees) of each tip in each tree in the first data set.

Either a single rooted tree or a list of rooted trees, of class "phylo", corresponding to the second data set on which an SBM model is to be fitted. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

Numeric vector listing the latitude (in decimal degrees) of each tip in each tree in the second data set, similarly to tip_latitudes1.

Numeric vector listing the longitude (in decimal degrees) of each tip in each tree in the second data set, similarly to tip_longitudes1.

Strictly positive numeric, specifying the radius of the sphere. For Earth, the mean radius is 6371 km.

Logical, specifying whether to estimate the diffusivity based on a planar approximation of the SBM model, i.e. by assuming that geographic distances between tips are as if tips are distributed on a 2D cartesian plane. This approximation is only accurate if geographical distances between tips are small compared to the sphere's radius.

Logical, specifying whether to only compare immediate sister tips, i.e., tips connected through a single parental node.

Logical, specifying whether to only compare tips at distinct geographic locations.

Numeric, specifying the minimum allowed time (distance from root) of the most recent common ancestor (MRCA) of sister tips considered in the fitting. In other words, an independent contrast is only considered if the two sister tips' MRCA has at least this distance from the root. Set min_MRCA_time<=0 to disable this filter.

Numeric, specifying the maximum allowed age (distance from youngest tip) of the MRCA of sister tips considered in the fitting. In other words, an independent contrast is only considered if the two sister tips' MRCA has at most this age (time to present). Set max_MRCA_age=Inf to disable this filter.

Numeric, maximum allowed geodistance for an independent contrast to be included in the SBM fitting. Set max_phylodistance=Inf to disable this filter.

Non-negative numeric, specifying the minimum possible diffusivity. If NULL, this is automatically chosen.

Non-negative numeric, specifying the maximum possible diffusivity. If NULL, this is automatically chosen.

Integer, specifying the number of parametric bootstraps to perform for calculating the confidence intervals of SBM diffusivities fitted to each data set. If <=0, no bootstrapping is performed.

Integer, specifying the number of simulations to perform for assessing the statistical significance of the linear difference and log-transformed difference between the diffusivities fitted to the two data sets, i.e. of |D_1-D_2| and of |\log(D_1)-\log(D_2)|. Set to 0 to not calculate statistical significances. See below for additional details.

SBM probability density functor object. Used internally and for debugging purposes. Unless you know what you're doing, you should keep this NULL.

Logical, specifying whether to print progress report messages to the screen.

Character, specifying a prefix to include in front of progress report messages on each line. Only relevant if verbose==TRUE.

Logical, indicating whether the fitting was successful for both data sets. If FALSE, then an additional return variable, error, will contain a description of the error; in that case all other return variables may be undefined.

A named list containing the fitting results for the first data set, in the same format as returned by fit_sbm_const. In particular, the diffusivity fitted to the first data set will be stored in fit1$diffusivity.

A named list containing the fitting results for the second data set, in the same format as returned by fit_sbm_const. In particular, the diffusivity fitted to the second data set will be stored in fit2$diffusivity.

The absolute difference between the two diffusivities, i.e. |D_1-D_2|.

The absolute difference between the two log-transformed diffusivities, i.e. |\log(D_1)-\log(D_2)|.

Numeric, statistical significance of the observed lin-difference under the null hypothesis that the two data sets were generated by a common SBM model. Only returned if Nsignificance>0.

Numeric, statistical significance of the observed log-difference under the null hypothesis that the two data sets were generated by a common SBM model. Only returned if Nsignificance>0.

A named list containing the fitting results for the two data sets combined, in the same format as returned by fit_sbm_const. The common diffusivity, fit_common$diffusivity is used for the random simulations when assessing the statistical significance of the lin-difference and log-difference of the separately fitted diffusivities. Only returned if Nsignificance>0.

Either a single rooted tree or a list of rooted trees, of class "phylo". The root of each tree is assumed to be the unique node with no incoming edge. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

Numeric state of each trait at each tip in each tree. If trees was a single tree, then tip_states must either be a numeric vector of size Ntips or a 2D numeric matrix of size Ntips x Ntraits, listing the trait states for each tip in the tree. If trees is a list of Ntrees trees, then tip_states must be a list of length Ntrees, each element of which lists the trait states for the corresponding tree (as a vector or 2D matrix, similarly to the single-tree case).

Logical, specifying whether diffusion should be assumed to be isotropic (i.e., independent of the direction). Hence, if isotropic=TRUE, then the diffusivity matrix is forced to be diagonal, with all entries being equal. If isotropic=FALSE, an arbitrary diffusivity matrix is fitted (i.e., the diffusivity matrix is only constrained to be symmetric and non-negative definite).

Integer, specifying the number of parametric bootstraps to perform for calculating the confidence intervals. If <=0, no bootstrapping is performed.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Logical, indicating whether the fitting was successful. If FALSE, then an additional return variable, error, will contain a description of the error; in that case all other return variables may be undefined.

Either a single non-negative number (if tip_states was a vector) or a 2D quadratic non-negative-definite matrix (if tip_states was a 2D matrix). The fitted diffusivity matrix of the multivariate Brownian motion model.

The log-likelihood of the fitted model, given the provided tip states data.

The AIC (Akaike Information Criterion) of the fitted model.

The BIC (Bayesian Information Criterion) of the fitted model.

Integer, number of independent contrasts used to estimate the diffusivity. This corresponds to the number of independent data points used.

Either a single numeric or a 2D numeric matrix of size Ntraits x Ntraits, listing the estimated standard errors of the estimated diffusivity, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Either a single numeric or a 2D numeric matrix of size Ntraits x Ntraits, listing the lower bounds of the 50% confidence interval for the estimated diffusivity (25-75% percentile), based on parametric bootstrapping. Only returned if Nbootstraps>0.

Either a single numeric or a 2D numeric matrix of size Ntraits x Ntraits, listing the upper bound of the 50% confidence interval for the estimated diffusivity, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Either a single numeric or a 2D numeric matrix of size Ntraits x Ntraits, listing the lower bound of the 95% confidence interval for the estimated diffusivity (2.5-97.5% percentile), based on parametric bootstrapping. Only returned if Nbootstraps>0.

Either a single numeric or a 2D numeric matrix of size Ntraits x Ntraits, listing the upper bound of the 95% confidence interval for the estimated diffusivity, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric between 0 and 1, estimated consistency of the data with the fitted model. If L denotes the loglikelihood of new data generated by the fitted model (under the same model) and M denotes the expectation of L, then consistency is the probability that |L-M| will be greater or equal to |X-M|, where X is the loglikelihood of the original data under the fitted model. Only returned if Nbootstraps>0. A low consistency (e.g., <0.05) indicates that the fitted model is a poor description of the data. See Lindholm et al. (2019) for background.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated reconstructed phylogeny of a set of extant sampled species.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which rho is defined. If age0>0, then rho0 refers to the sampling fraction at age age0, i.e. the fraction of lineages extant at age0 that are included in the tree. See below for more details.

Numeric vector, listing ages in ascending order, on which \lambda and \mu are allowed to vary independently. This grid must cover age0. If splines_degree>0 (see option below) then the age grid must also cover oldest_age. If NULL or of length <=1 (regardless of value), then \lambda and \mu are assumed to be time-independent.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted \lambda at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted \lambda at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric vector of length Ngrid, or a single numeric, specifying lower bounds for the fitted \mu at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted \mu at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric, specifying a lower bound for the fitted sampling fraction \rho (fraction of extant species included in the tree).

Numeric, specifying an upper bound for the fitted sampling fraction \rho.

Initial guess for \lambda at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \lambda at all ages) or a numeric vector of size Ngrid specifying a separate guess for \lambda at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for \mu at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \mu at all ages) or a numeric vector of size Ngrid specifying a separate guess for \mu at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Numeric, specifying an initial guess for the sampling fraction \rho at age0. Setting this to NULL or NA is the same as setting it to 1.

Optional fixed (i.e. non-fitted) \lambda values on one or more age-grid points. Either NULL (\lambda is not fixed anywhere), or a single numeric (\lambda fixed to the same value at all grid points) or a numeric vector of size Ngrid (\lambda fixed on one or more age-grid points, use NA for non-fixed values).

Optional fixed (i.e. non-fitted) \mu values on one or more age-grid points. Either NULL (\mu is not fixed anywhere), or a single numeric (\mu fixed to the same value at all grid points) or a numeric vector of size Ngrid (\mu fixed on one or more age-grid points, use NA for non-fixed values).

Numeric between 0 and 1, optionallly specifying a fixed value for the sampling fraction \rho. If NULL or NA, the sampling fraction \rho is estimated, however note that this may not always be meaningful (Stadler 2009, Stadler 2013).

Logical, specifying whether \lambda should be assumed constant across the grid, i.e. time-independent. Setting const_lambda=TRUE reduces the number of free (i.e., independently fitted) parameters. If \lambda is fixed on some grid points (i.e. via fixed_lambda), then only the non-fixed lambdas are assumed to be identical to one another.

Logical, specifying whether \mu should be assumed constant across the grid, i.e. time-independent. Setting const_mu=TRUE reduces the number of free (i.e., independently fitted) parameters. If \mu is fixed on some grid points (i.e. via fixed_mu), then only the non-fixed mus are assumed to be identical to one another.

Integer between 0 and 3 (inclusive), specifying the polynomial degree of \lambda and \mu between age-grid points. If 0, then \lambda and \mu are considered piecewise constant, if 1 then \lambda and \mu are considered piecewise linear, if 2 or 3 then \lambda and \mu are considered to be splines of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended, despite the fact that it has been historically popular. The case splines_degree=0 is also known as “skyline” model.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier. "none" is generally not recommended.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

The log-likelihood of the fitted model for the given timetree.

Numeric vector of size Ngrid, listing fitted or fixed speciation rates \lambda at each age-grid point. Between grid points \lambda should be interpreted as a piecewise polynomial function (natural spline) of degree splines_degree; to evaluate this function at arbitrary ages use the castor routine evaluate_spline.

Numeric vector of size Ngrid, listing fitted or fixed extinction rates \mu at each age-grid point. Between grid points \mu should be interpreted as a piecewise polynomial function (natural spline) of degree splines_degree; to evaluate this function at arbitrary ages use the castor routine evaluate_spline.

Numeric, specifying the fitted or fixed sampling fraction \rho.

Numeric vector of size Ngrid, specifying the initial guess for \lambda at each age-grid point.

Numeric vector of size Ngrid, specifying the initial guess for \mu at each age-grid point.

Numeric, specifying the initial guess for \rho.

The age-grid on which \lambda and \mu are defined. This will be the same as the provided age_grid, unless the latter was NULL or of length <=1.

Integer, number of free (i.e., independently) fitted parameters. If none of the \lambda, \mu and \rho were fixed, and const_lambda=FALSE and const_mu=FALSE, then NFP will be equal to 2*Ngrid+1.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Character, specifying what conditioning was root for the likelihood (e.g. "crown" or "stem").

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated reconstructed phylogeny of a set of extant sampled species.

Numeric vector, specifying fixed values for a some or all model parameters. For fitted (i.e., non-fixed) parameters, use NaN or NA. For example, the vector c(1.5,NA,40) specifies that the 1st and 3rd model parameters are fixed at the values 1.5 and 40, respectively, while the 2nd parameter is to be fitted. The length of this vector defines the total number of model parameters. If entries in this vector are named, the names are taken as parameter names. Names should be included if you'd like returned parameter vectors to have named entries, or if the functions lambda, mu or rho query parameter values by name (as opposed to numeric index).

Numeric vector of size NP, specifying a first guess for the value of each model parameter. For fixed parameters, guess values are ignored. Can be NULL only if all model parameters are fixed.

Optional numeric vector of size NP, specifying lower bounds for model parameters. If of size 1, the same lower bound is applied to all parameters. Use -Inf to omit a lower bound for a parameter. If NULL, no lower bounds are applied. For fixed parameters, lower bounds are ignored.

Optional numeric vector of size NP, specifying upper bounds for model parameters. If of size 1, the same upper bound is applied to all parameters. Use +Inf to omit an upper bound for a parameter. If NULL, no upper bounds are applied. For fixed parameters, upper bounds are ignored.

Optional numeric vector of size NP, specifying typical scales for model parameters. If of size 1, the same scale is assumed for all parameters. If NULL, scales are determined automatically. For fixed parameters, scales are ignored. It is strongly advised to provide reasonable scales, as this facilitates the numeric optimization algorithm.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which rho is defined. If age0>0, then rho0 refers to the sampling fraction at age age0, i.e. the fraction of lineages extant at age0 that are included in the tree. See below for more details.

Function specifying the speciation rate at any given age (time before present) and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with strictly positive entries. Can also be a single number (i.e., lambda is fixed).

Function specifying the extinction rate at any given age and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with non-negative entries. Can also be a single number (i.e., mu is fixed).

Function specifying the sampling fraction (fraction of extant species sampled at age0) for any given parameter values. This function must take exactly one argument, a numeric vector of size NP (parameter values), and return a numeric between 0 (exclusive) and 1 (inclusive). Can also be a single number (i.e., rho0 is fixed).

Numeric vector, specifying ages at which the lambda and mu functionals should be evaluated. This age grid must be fine enough to capture the possible variation in \lambda and \mu over time, within the permissible parameter range. If of size 1, then lambda & mu are assumed to be time-independent. Listed ages must be strictly increasing, and must cover at least the full considered age interval (from 0 to oldest_age). Can also be NULL or a vector of size 1, in which case the speciation rate and extinction rate is assumed to be time-independent.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier. "none" is generally not recommended.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of times to attempt finding a valid start point (per trial) before giving up on that trial. Randomly choosen extreme start parameters may occasionally result in Inf/undefined likelihoods, so this option allows the algorithm to keep looking for valid starting points.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Integer, specifying the number of parametric bootstraps to perform for estimating standard errors and confidence intervals of estimated model parameters. Set to 0 for no bootstrapping.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Character, specifying which optimization algorithm to use. Either "nlminb" or "subplex" are allowed.

Named list containing options for the nlminb or subplex optimization routine, depending on the choice of fit_algorithm. For example, for "nlminb" commonly modified options are iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package or of nloptr in the nloptr package.

Optional numeric matrix having NP columns and an arbitrary number of rows, listing combinations of parameter values of particular interest and for which the log-likelihoods should be returned. This may be used for diagnostic purposes, e.g., to examine the shape of the likelihood function.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Logical, specifying whether to print detailed information (such as model likelihoods) at every iteration of the fitting routine. For debugging purposes mainly.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

Numeric vector of size NP (number of model parameters), listing all fitted or fixed model parameters in their standard order (see details above). If param_names was provided, elements in fitted_params will be named.

Numeric vector of size NP, listing guessed or fixed values for all model parameters in their standard order. If param_names was provided, elements in param_guess will be named.

The log-likelihood of the fitted model for the given timetree.

Integer, number of fitted (i.e., non-fixed) model parameters.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Character, specifying what conditioning was root for the likelihood (e.g. "crown" or "stem").

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

Numeric vector of size Ntrials, listing the initial objective values (e.g., loglikelihoods) for each fitting trial, i.e. at the start parameter values.

Numeric vector of size Ntrials, listing the final maximized objective values (e.g., loglikelihoods) for each fitting trial.

Integer vector of size Ntrials, listing the number of start attempts for each fitting trial, until a starting point with valid likelihood was found.

Integer vector of size Ntrials, listing the number of iterations needed for each fitting trial.

Integer vector of size Ntrials, listing the number of likelihood evaluations needed for each fitting trial.

Numeric vector of size NP, estimated standard error of the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, median the estimated parameters across parametric bootstraps. Only returned if Nbootstraps>0.

Numeric vector of size NP, lower bound of the 50% confidence interval (25-75% percentile) for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, upper bound of the 50% confidence interval for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, lower bound of the 95% confidence interval (2.5-97.5% percentile) for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, upper bound of the 95% confidence interval for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric between 0 and 1, estimated consistency of the data with the fitted model (Lindholm et al. 2019). See the documentation of fit_hbds_model_on_grid for an explanation.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting. If age0>0, the tree essentially is trimmed at age0, omitting anything younger than age0, and the PDR and \rho\lambda_o are fitted to the trimmed tree while shifting time appropriately.

Numeric vector, listing alternative grid sizes to consider.

Logical, specifying whether to use uniform time grids (equal time intervals) or non-uniform time grids (more grid points towards the present, where more data are available).

Character, specifying which criterion to use for selecting the best grid. Options are "AIC" and "BIC".

Logical, whether to try all grid sizes before choosing the best one. If FALSE, the grid size is gradually increased until the selection criterio (e.g., AIC) starts becoming worse, at which point the search is halted. This avoids fitting models with excessive grid sizes when an optimum already seems to have been found at a smaller grid size.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted PDR at each point in the age grid. If a single numeric, the same lower bound applies at all ages. Note that in general the PDR may be negative as well as positive.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted PDR at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Strictly positive numeric, specifying the lower bound for the fitted \rho\lambda_o (sampling fraction times present-day extinction rate).

Strictly positive numeric, specifying the upper bound for the fitted \rho\lambda_o. Set to +Inf to omit this upper bound.

Initial guess for the PDR at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing a constant PDR at all ages), or a function handle (for generating guesses at each grid point; this function may also return NA at some time points for which a guess shall be computed automatically).

Numeric, specifying an initial guess for the product \rho\lambda_o. If NULL, a guess will be computed automatically.

Optional fixed (i.e. non-fitted) PDR values. Either NULL (none of the PDR values are fixed) or a function handle specifying the PDR for any arbitrary age (PDR will be fixed at any age for which this function returns a finite number). The function fixed_PDR() need not return finite values for all times, in fact doing so would mean that the PDR is not fitted anywhere.

Numeric, optionally specifying a fixed value for the product \rho\lambda_o. If NULL or NA, the product \rho\lambda_o is estimated.

Integer between 0 and 3 (inclusive), specifying the polynomial degree of the PDR between age-grid points. If 0, then the PDR is considered to be piecewise constant, if 1 then the PDR is considered piecewise linear, if 2 or 3 then the PDR is considered to be a spline of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying an optional number of bootstrap samplings to perform, for estimating standard errors and confidence intervals of maximum-likelihood fitted parameters. If 0, no bootstrapping is performed. Typical values are 10-100. At each bootstrap sampling, a random timetree is generated under the birth-death model according to the fitted PDR and \rho\lambda_o, the parameters are estimated anew based on the generated tree, and subsequently compared to the original fitted parameters. Each bootstrap sampling will use roughly the same information and similar computational resources as the original maximum-likelihood fit (e.g., same number of trials, same optimization parameters, same initial guess, etc). Bootstrapping is only performed for the best grid size.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether the function executed successfully. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

A named list containing the fitting results for the best grid size. This list has the same structure as the one returned by fit_hbd_pdr_on_grid.

Numeric vector, listing the grid sizes as provided during the function call.

Numeric vector of the same length as grid_sizes, listing the AIC for each considered grid size. Note that some entries may be NA, if the corresponding grid sizes were not considered (if exhaustive=FALSE).

Numeric vector of the same length as grid_sizes, listing the BIC for each considered grid size. Note that some entries may be NA, if the corresponding grid sizes were not considered (if exhaustive=FALSE).

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which rholambda0 is defined. If age0>0, then rholambda0 refers to the product of the sampling fraction at age age0 and the speciation rate at age age0. See below for more details.

Numeric vector, listing ages in ascending order at which the PDR is allowed to vary independently. This grid must cover at least the age range from age0 to oldest_age. If NULL or of length <=1 (regardless of value), then the PDR is assumed to be time-independent.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted PDR at each point in the age grid. If a single numeric, the same lower bound applies at all ages. Use -Inf to omit lower bounds.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted PDR at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Strictly positive numeric, specifying the lower bound for the fitted \rho\lambda_o (sampling fraction times present-day extinction rate).

Strictly positive numeric, specifying the upper bound for the fitted \rho\lambda_o. Set to +Inf to omit this upper bound.

Initial guess for the PDR at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same PDR at all ages) or a numeric vector of size Ngrid specifying a separate guess at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Numeric, specifying an initial guess for the product \rho\lambda_o. If NULL, a guess will be computed automatically.

Optional fixed (i.e. non-fitted) PDR values on one or more age-grid points. Either NULL (PDR is not fixed anywhere), or a single numeric (PDR fixed to the same value at all grid points) or a numeric vector of size Ngrid (PDR fixed at one or more age-grid points, use NA for non-fixed values).

Numeric, optionally specifying a fixed value for the product \rho\lambda_o. If NULL or NA, the product \rho\lambda_o is estimated.

Integer between 0 and 3 (inclusive), specifying the polynomial degree of the PDR between age-grid points. If 0, then the PDR is considered piecewise constant, if 1 then the PDR is considered piecewise linear, if 2 or 3 then the PDR is considered to be a spline of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying an optional number of bootstrap samplings to perform, for estimating standard errors and confidence intervals of maximum-likelihood fitted parameters. If 0, no bootstrapping is performed. Typical values are 10-100. At each bootstrap sampling, a random timetree is generated under the birth-death model according to the fitted PDR and \rho\lambda_o, the parameters are estimated anew based on the generated tree, and subsequently compared to the original fitted parameters. Each bootstrap sampling will use roughly the same information and similar computational resources as the original maximum-likelihood fit (e.g., same number of trials, same optimization parameters, same initial guess, etc).

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

The log-likelihood of the fitted model for the given timetree.

Numeric vector of size Ngrid, listing fitted or fixed pulled diversification rates (PDR) at each age-grid point. Between grid points the fitted PDR should be interpreted as a piecewise polynomial function (natural spline) of degree splines_degree; to evaluate this function at arbitrary ages use the castor routine evaluate_spline.

Numeric, specifying the fitted or fixed product \rho\lambda(0).

Numeric vector of size Ngrid, specifying the initial guess for the PDR at each age-grid point.

Numeric, specifying the initial guess for \rho\lambda(0).

The age-grid on which the PDR is defined. This will be the same as the provided age_grid, unless the latter was NULL or of length <=1.

Integer, number of fitted (i.e., non-fixed) parameters. If none of the PDRs or \rho\lambda0 were fixed, this will be equal to Ngrid+1.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

If Nbootstraps>0, this will be a named list containing the elements PDR (numeric matrix of size Nbootstraps x Ngrid, listing the fitted PDR at each grid point and for each bootstrap) and rholambda0 (a numeric vector of size Nbootstraps, listing the fitted \rho\lambda_o for each bootstrap).

If Nbootstraps>0, this will be a named list containing the elements PDR (numeric vector of size Ngrid, listing bootstrap-estimated standard errors for the fitted PDRs) and rholambda0 (a single numeric, bootstrap-estimated standard error for the fitted \rho\lambda_o).

If Nbootstraps>0, this will be a named list containing the elements PDR (numeric vector of size Ngrid, listing median fitted PDRs across bootstraps) and rholambda0 (a single numeric, median fitted \rho\lambda_o across bootstraps).

If Nbootstraps>0, this will be a named list containing the elements PDR (numeric vector of size Ngrid, listing bootstrap-estimated lower bounds of the 50-percent confidence intervals for the fitted PDRs) and rholambda0 (a single numeric, bootstrap-estimated lower bound of the 50-percent confidence intervals for the fitted \rho\lambda_o).

Similar to CI50lower, listing upper bounds of 50-percentile confidence intervals.

Similar to CI50lower, listing lower bounds of 95-percentile confidence intervals.

Similar to CI95lower, listing upper bounds of 95-percentile confidence intervals.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Numeric vector, specifying fixed values for a some or all model parameters. For fitted (i.e., non-fixed) parameters, use NaN or NA. For example, the vector c(1.5,NA,40) specifies that the 1st and 3rd model parameters are fixed at the values 1.5 and 40, respectively, while the 2nd parameter is to be fitted. The length of this vector defines the total number of model parameters. If entries in this vector are named, the names are taken as parameter names. Names should be included if you'd like returned parameter vectors to have named entries, or if the functions PDR or rho query parameter values by name (as opposed to numeric index).

Numeric vector of size NP, specifying a first guess for the value of each model parameter. For fixed parameters, guess values are ignored. Can be NULL only if all model parameters are fixed.

Optional numeric vector of size NP, specifying lower bounds for model parameters. If of size 1, the same lower bound is applied to all parameters. Use -Inf to omit a lower bound for a parameter. If NULL, no lower bounds are applied. For fixed parameters, lower bounds are ignored.

Optional numeric vector of size NP, specifying upper bounds for model parameters. If of size 1, the same upper bound is applied to all parameters. Use +Inf to omit an upper bound for a parameter. If NULL, no upper bounds are applied. For fixed parameters, upper bounds are ignored.

Optional numeric vector of size NP, specifying typical scales for model parameters. If of size 1, the same scale is assumed for all parameters. If NULL, scales are determined automatically. For fixed parameters, scales are ignored. It is strongly advised to provide reasonable scales, as this facilitates the numeric optimization algorithm.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which rholambda0 is defined. If age0>0, then rholambda0 refers to the product of the sampling fraction at age age0 and the speciation rate at age age0. See below for more details.

Function specifying the pulled diversification rate at any given age (time before present) and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument. Can also be a single number (i.e., PDR is fixed).

Function specifying the product \rho\lambda_o (sampling fraction times speciation rate at age0) for any given parameter values. This function must take exactly one argument, a numeric vector of size NP (parameter values), and return a strictly positive numeric. Can also be a single number (i.e., rholambda0 is fixed).

Numeric vector, specifying ages at which the PDR function should be evaluated. This age grid must be fine enough to capture the possible variation in the PDR over time, within the permissible parameter range. If of size 1, then the PDR is assumed to be time-independent. Listed ages must be strictly increasing, and must cover at least the full considered age interval (from age0 to oldest_age). Can also be NULL or a vector of size 1, in which case the PDR is assumed to be time-independent.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of times to attempt finding a valid start point (per trial) before giving up on that trial. Randomly choosen extreme start parameters may occasionally result in Inf/undefined likelihoods, so this option allows the algorithm to keep looking for valid starting points.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

Numeric vector of size NP (number of model parameters), listing all fitted or fixed model parameters in their standard order (see details above). If param_names was provided, elements in fitted_params will be named.

Numeric vector of size NP, listing guessed or fixed values for all model parameters in their standard order.

The log-likelihood of the fitted model for the given timetree.

Integer, number of fitted (i.e., non-fixed) model parameters.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

Numeric vector of size Ntrials, listing the initial objective values (e.g., loglikelihoods) for each fitting trial, i.e. at the start parameter values.

Numeric vector of size Ntrials, listing the final maximized objective values (e.g., loglikelihoods) for each fitting trial.

Integer vector of size Ntrials, listing the number of start attempts for each fitting trial, until a starting point with valid likelihood was found.

Integer vector of size Ntrials, listing the number of iterations needed for each fitting trial.

Integer vector of size Ntrials, listing the number of likelihood evaluations needed for each fitting trial.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting. If age0>0, the tree essentially is trimmed at age0, omitting anything younger than age0, and the PSR is fitted to the trimmed tree while shifting time appropriately.

Numeric vector, listing alternative grid sizes to consider.

Logical, specifying whether to use uniform time grids (equal time intervals) or non-uniform time grids (more grid points towards the present, where more data are available).

Character, specifying which criterion to use for selecting the best grid. Options are "AIC" and "BIC".

Logical, whether to try all grid sizes before choosing the best one. If FALSE, the grid size is gradually increased until the selection criterio (e.g., AIC) starts becoming worse, at which point the search is halted. This avoids fitting models with excessive grid sizes when an optimum already seems to have been found at a smaller grid size.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted PSR at each point in the age grid. If a single numeric, the same lower bound applies at all ages. Note that the PSR is never negative.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted PSR at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Initial guess for the PSR at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing a constant PSR at all ages), or a function handle (for generating guesses at each grid point; this function may also return NA at some time points for which a guess shall be computed automatically).

Optional fixed (i.e. non-fitted) PSR values. Either NULL (none of the PSR values are fixed) or a function handle specifying the PSR for any arbitrary age (PSR will be fixed at any age for which this function returns a finite number). The function fixed_PSR() need not return finite values for all times, in fact doing so would mean that the PSR is not fitted anywhere.

Integer between 0 and 3 (inclusive), specifying the polynomial degree of the PSR between age-grid points. If 0, then the PSR is considered piecewise constant, if 1 then the PSR is considered piecewise linear, if 2 or 3 then the PSR is considered to be a spline of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying an optional number of bootstrap samplings to perform, for estimating standard errors and confidence intervals of maximum-likelihood fitted parameters. If 0, no bootstrapping is performed. Typical values are 10-100. At each bootstrap sampling, a random timetree is generated under the birth-death model according to the fitted PSR, the parameters are estimated anew based on the generated tree, and subsequently compared to the original fitted parameters. Each bootstrap sampling will use roughly the same information and similar computational resources as the original maximum-likelihood fit (e.g., same number of trials, same optimization parameters, same initial guess, etc). Bootstrapping is only performed for the best grid size.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether the function executed successfully. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

A named list containing the fitting results for the best grid size. This list has the same structure as the one returned by fit_hbd_psr_on_grid.

Numeric vector, listing the grid sizes as provided during the function call.

Numeric vector of the same length as grid_sizes, listing the AIC for each considered grid size. Note that some entries may be NA, if the corresponding grid sizes were not considered (if exhaustive=FALSE).

Numeric vector of the same length as grid_sizes, listing the BIC for each considered grid size. Note that some entries may be NA, if the corresponding grid sizes were not considered (if exhaustive=FALSE).

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting. If age0>0, the tree essentially is trimmed at age0, omitting anything younger than age0, and the PSR is fitted to the trimmed tree while shifting time appropriately.

Numeric vector, listing ages in ascending order at which the PSR is allowed to vary independently. This grid must cover at least the age range from age0 to oldest_age. If NULL or of length <=1 (regardless of value), then the PSR is assumed to be time-independent.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted PSR at each point in the age grid. If a single numeric, the same lower bound applies at all ages. Note that the PSR is never negative.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted PSR at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Initial guess for the PSR at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same PSR at all ages) or a numeric vector of size Ngrid specifying a separate guess at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Optional fixed (i.e. non-fitted) PSR values on one or more age-grid points. Either NULL (PSR is not fixed anywhere), or a single numeric (PSR fixed to the same value at all grid points) or a numeric vector of size Ngrid (PSR fixed at one or more age-grid points, use NA for non-fixed values).

Integer between 0 and 3 (inclusive), specifying the polynomial degree of the PSR between age-grid points. If 0, then the PSR is considered piecewise constant, if 1 then the PSR is considered piecewise linear, if 2 or 3 then the PSR is considered to be a spline of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying an optional number of bootstrap samplings to perform, for estimating standard errors and confidence intervals of maximum-likelihood fitted parameters. If 0, no bootstrapping is performed. Typical values are 10-100. At each bootstrap sampling, a random timetree is generated under the birth-death model according to the fitted PSR, the parameters are estimated anew based on the generated tree, and subsequently compared to the original fitted parameters. Each bootstrap sampling will use roughly the same information and similar computational resources as the original maximum-likelihood fit (e.g., same number of trials, same optimization parameters, same initial guess, etc).

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Logical, specifying whether to print detailed information (such as model likelihoods) at every iteration of the fitting routine. For debugging purposes mainly.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

The log-likelihood of the fitted model for the given timetree.

Numeric vector of size Ngrid, listing fitted or fixed pulled speciation rates (PSR) at each age-grid point. Between grid points the fitted PSR should be interpreted as a piecewise polynomial function (natural spline) of degree splines_degree; to evaluate this function at arbitrary ages use the castor routine evaluate_spline.

Numeric vector of size Ngrid, specifying the initial guess for the PSR at each age-grid point.

The age-grid on which the PSR is defined. This will be the same as the provided age_grid, unless the latter was NULL or of length <=1.

Integer, number of fitted (i.e., non-fixed) parameters. If none of the PSRs were fixed, this will be equal to Ngrid.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters, and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

If Nbootstraps>0, this will be a numeric matrix of size Nbootstraps x Ngrid, listing the fitted PSR at each grid point and for each bootstrap.

If Nbootstraps>0, this will be a numeric vector of size NGrid, listing bootstrap-estimated standard errors for the fitted PSR at each grid point.

If Nbootstraps>0, this will be a numeric vector of size Ngrid, listing bootstrap-estimated lower bounds of the 50-percent confidence intervals for the fitted PSR at each grid point.

Similar to CI50lower, listing upper bounds of 50-percentile confidence intervals.

Similar to CI50lower, listing lower bounds of 95-percentile confidence intervals.

Similar to CI95lower, listing upper bounds of 95-percentile confidence intervals.

A rooted ultrametric timetree of class "phylo", representing the time-calibrated phylogeny of a set of extant sampled species.

Numeric vector, specifying fixed values for a some or all model parameters. For fitted (i.e., non-fixed) parameters, use NaN or NA. For example, the vector c(1.5,NA,40) specifies that the 1st and 3rd model parameters are fixed at the values 1.5 and 40, respectively, while the 2nd parameter is to be fitted. The length of this vector defines the total number of model parameters. If entries in this vector are named, the names are taken as parameter names. Names should be included if you'd like returned parameter vectors to have named entries, or if the function PSR queries parameter values by name (as opposed to numeric index).

Numeric vector of size NP, specifying a first guess for the value of each model parameter. For fixed parameters, guess values are ignored. Can be NULL only if all model parameters are fixed.

Optional numeric vector of size NP, specifying lower bounds for model parameters. If of size 1, the same lower bound is applied to all parameters. Use -Inf to omit a lower bound for a parameter. If NULL, no lower bounds are applied. For fixed parameters, lower bounds are ignored.

Optional numeric vector of size NP, specifying upper bounds for model parameters. If of size 1, the same upper bound is applied to all parameters. Use +Inf to omit an upper bound for a parameter. If NULL, no upper bounds are applied. For fixed parameters, upper bounds are ignored.

Optional numeric vector of size NP, specifying typical scales for model parameters. If of size 1, the same scale is assumed for all parameters. If NULL, scales are determined automatically. For fixed parameters, scales are ignored. It is strongly advised to provide reasonable scales, as this facilitates the numeric optimization algorithm.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is taken as the stem age, and the classical formula by Morlon et al. (2011) is used. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBD model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Non-negative numeric, specifying the youngest age (time before present) to consider for fitting, and with respect to which PSR is defined. See below for more details. Most users will typically keep age0=0.

Function specifying the pulled speciation rate at any given age (time before present) and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument.

Numeric vector, specifying ages at which the PSR function should be evaluated. This age grid must be fine enough to capture the possible variation in the PSR over time, within the permissible parameter range. If of size 1, then the PSR is assumed to be time-independent. Listed ages must be strictly increasing, and must cover at least the full considered age interval (from age0 to oldest_age). Can also be NULL or a vector of size 1, in which case the PSR is assumed to be time-independent.

Character, either "crown", "stem", "auto", "stemN" or "crownN" (where N is an integer >=2), specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root at that time. If "stem", the likelihood is conditioned on the survival of the stem lineage, with the process having started at oldest_age. Note that "crown" and "crownN"" really only make sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. If "stem2", the condition is that the process yielded at least two sampled tips, and similarly for "stem3" etc. If "crown3", the condition is that a splitting occurred at the root age, both child clades survived, and in total yielded at least 3 sampled tips (and similarly for "crown4" etc). If "auto", the condition is chosen according to the recommendations mentioned earlier.

Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time, when calculating the likelihood. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of times to attempt finding a valid start point (per trial) before giving up on that trial. Randomly choosen extreme start parameters may occasionally result in Inf/undefined likelihoods, so this option allows the algorithm to keep looking for valid starting points.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Logical, specifying whether to print detailed information (such as model likelihoods) at every iteration of the fitting routine. For debugging purposes mainly.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

Numeric vector of size NP (number of model parameters), listing all fitted or fixed model parameters in their standard order (see details above). If param_names was provided, elements in fitted_params will be named.

Numeric vector of size NP, listing guessed or fixed values for all model parameters in their standard order.

The log-likelihood of the fitted model for the given timetree.

Integer, number of fitted (i.e., non-fixed) model parameters.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

Numeric vector of size Ntrials, listing the initial objective values (e.g., loglikelihoods) for each fitting trial, i.e. at the start parameter values.

Numeric vector of size Ntrials, listing the final maximized objective values (e.g., loglikelihoods) for each fitting trial.

Integer vector of size Ntrials, listing the number of start attempts for each fitting trial, until a starting point with valid likelihood was found.

Integer vector of size Ntrials, listing the number of iterations needed for each fitting trial.

Integer vector of size Ntrials, listing the number of likelihood evaluations needed for each fitting trial.

A timetree of class "phylo", representing the time-calibrated reconstructed phylogeny of a set of extant and/or extinct species. Tips of the tree are interpreted as terminally sampled lineages, while monofurcating nodes are interpreted as non-terminally sampled lineages, i.e., lineages sampled at some past time point and with subsequently sampled descendants.

Positive numeric, specifying the age of the tree's root. Can be used to define a time offset, e.g. if the last tip was not actually sampled at the present. If NULL, this will be calculated from the tree and it will be assumed that the last tip was sampled at the present.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is interpreted as the stem age. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBDS model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Numeric vector, listing ages in ascending order, on which \lambda, \mu, \psi and \kappa are fitted and allowed to vary independently. This grid must cover at least the age range from the present (age 0) to oldest_age. If NULL or of length <=1 (regardless of value), then \lambda, \mu, \psi and \kappa are assumed to be time-independent.

Optional numeric vector, listing ages (in ascending order) at which concentrated sampling attempts (CSAs) occurred. If NULL, it is assumed that no concentrated sampling attempts took place and that all tips were sampled according to the continuous sampling rate psi.

Numeric vector of length Ngrid (=max(1,length(age_grid))), or a single numeric, specifying lower bounds for the fitted speciation rate \lambda at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted speciation rate \lambda at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric vector of length Ngrid, or a single numeric, specifying lower bounds for the fitted extinction rate \mu at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted extinction rate \mu at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric vector of length Ngrid, or a single numeric, specifying lower bounds for the fitted Poissonian sampling rate \psi at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted Poissonian sampling rate \psi at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric vector of length Ngrid, or a single numeric, specifying lower bounds for the fitted retention probability \kappa at each point in the age grid. If a single numeric, the same lower bound applies at all ages.

Numeric vector of length Ngrid, or a single numeric, specifying upper bounds for the fitted retention probability \kappa at each point in the age grid. If a single numeric, the same upper bound applies at all ages. Use +Inf to omit upper bounds.

Numeric vector of length NCSA (=length(CSA_ages)), or a single numeric, specifying lower bounds for the fitted sampling probabilities \rho_1,..,\rho_m at each concentrated sampling attempt. If a single numeric, the same lower bound applies at all CSAs. Note that, since \rho_1, \rho_2, ... are probabilities, min_CSA_probs should not be negative.

Numeric vector of length NCSA, or a single numeric, specifying upper bounds for the fitted sampling probabilities \rho_1, \rho_2, ... at each concentrated sampling attempt. If a single numeric, the same upper bound applies at all CSAs. Note that, since \rho_1, \rho_2, ... are probabilities, max_CSA_probs should not be greater than 1.

Numeric vector of length NCSA, or a single numeric, specifying lower bounds for the fitted retention probabilities \kappa_1, \kappa_2, ... at each concentrated sampling attempt. If a single numeric, the same lower bound applies at all CSAs. Note that, since \kappa_1, \kappa_2, ... are probabilities, min_CSA_kappas should not be negative.

Numeric vector of length NCSA, or a single numeric, specifying upper bounds for the fitted sampling probabilities \kappa_1, \kappa_2, ... at each concentrated sampling attempt. If a single numeric, the same upper bound applies at all CSAs. Note that, since \kappa_1, \kappa_2, .. are probabilities, max_CSA_kappas should not be greater than 1.

Initial guess for \lambda at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \lambda at all ages) or a numeric vector of size Ngrid specifying a separate guess for \lambda at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for \mu at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \mu at all ages) or a numeric vector of size Ngrid specifying a separate guess for \mu at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for \psi at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \psi at all ages) or a numeric vector of size Ngrid specifying a separate guess for \psi at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for \kappa at each age-grid point. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same \kappa at all ages) or a numeric vector of size Ngrid specifying a separate guess for \kappa at each age-grid point. To omit an initial guess for some but not all age-grid points, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for the \rho_1, \rho_2, ... at each concentrated sampling attempt. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same value at every CSA) or a numeric vector of size NCSA specifying a separate guess at each CSA. To omit an initial guess for some but not all CSAs, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Initial guess for the \kappa_1, \kappa_2, ... at each concentrated sampling attempt. Either NULL (an initial guess will be computed automatically), or a single numeric (guessing the same value at every CSA) or a numeric vector of size NCSA specifying a separate guess at each CSA. To omit an initial guess for some but not all CSAs, set their guess values to NA. Guess values are ignored for non-fitted (i.e., fixed) parameters.

Optional fixed (i.e. non-fitted) \lambda values on one or more age-grid points. Either NULL (\lambda is not fixed anywhere), or a single numeric (\lambda fixed to the same value at all grid points) or a numeric vector of size Ngrid (if fixed_age_grid=NULL; \lambda fixed on one or more age-grid points, use NA for non-fixed values) or a numeric vector of the same size as fixed_age_grid (if fixed_age_grid!=NULL, in which case all entries in fixed_lambda must be finite numbers).

Optional fixed (i.e. non-fitted) \mu values on one or more age-grid points. Either NULL (\mu is not fixed anywhere), or a single numeric (\mu fixed to the same value at all grid points) or a numeric vector of size Ngrid (if fixed_age_grid=NULL; \mu fixed on one or more age-grid points, use NA for non-fixed values) or a numeric vector of the same size as fixed_age_grid (if fixed_age_grid!=NULL, in which case all entries in fixed_mu must be finite numbers).

Optional fixed (i.e. non-fitted) \psi values on one or more age-grid points. Either NULL (\psi is not fixed anywhere), or a single numeric (\psi fixed to the same value at all grid points) or a numeric vector of size Ngrid (if fixed_age_grid=NULL; \psi fixed on one or more age-grid points, use NA for non-fixed values) or a numeric vector of the same size as fixed_age_grid (if fixed_age_grid!=NULL, in which case all entries in fixed_psi must be finite numbers).

Optional fixed (i.e. non-fitted) \kappa values on one or more age-grid points. Either NULL (\kappa is not fixed anywhere), or a single numeric (\kappa fixed to the same value at all grid points) or a numeric vector of size Ngrid (if fixed_age_grid=NULL; \kappa fixed on one or more age-grid points, use NA for non-fixed values) or a numeric vector of the same size as fixed_age_grid (if fixed_age_grid!=NULL, in which case all entries in fixed_kappa must be finite numbers).

Optional fixed (i.e. non-fitted) \rho_1, \rho_2, ... values on one or more age-grid points. Either NULL (none of the \rho_1, \rho_2,... are fixed), or a single numeric (\rho_1, \rho_2,... are fixed to the same value at all CSAs) or a numeric vector of size NCSA (one or more of the \rho_1, \rho_2, ... are fixed, use NA for non-fixed values).

Optional fixed (i.e. non-fitted) \kappa_1, \kappa_2, ... values on one or more age-grid points. Either NULL (none of the \kappa_1, \kappa_2,... are fixed), or a single numeric (\kappa_1, \kappa_2,... are fixed to the same value at all CSAs) or a numeric vector of size NCSA (one or more of the \kappa_1, \kappa_2, ... are fixed, use NA for non-fixed values).

Optional numeric vector, specifying an age grid on which fixed_lambda, fixed_mu, fixed_psi and fixed_kappa (whichever is provided) are defined instead of on the age_grid. If fixed_age_grid is provided, then each of fixed_lambda, fixed_mu, fixed_psi and fixed_kappa must be defined (i.e. have a finite non-negative value) on every point in fixed_age_grid. Entries in fixed_age_grid must be in ascending order and must cover at least the ages 0 to oldest_age.

This option may be useful if you want to fit some parameters on a coarse grid, but want to specify (fix) some other parameters on a much finer grid. Also note that if fixed_age_grid is used, all parameters lambda, mu, psi and kappa are internally re-interpolated onto fixed_age_grid when evaluating the likelihood; hence, in general fixed_age_grid should be much finer than age_grid. In most situations you would probably want to keep the default fixed_age_grid=NULL.

Logical, specifying whether \lambda should be assumed constant across the grid, i.e. time-independent. Setting const_lambda=TRUE reduces the number of free (i.e., independently fitted) parameters. If \lambda is fixed on some grid points (i.e. via fixed_lambda), then only the non-fixed lambdas are assumed to be identical to one another.

Logical, specifying whether \mu should be assumed constant across the grid, i.e. time-independent. Setting const_mu=TRUE reduces the number of free (i.e., independently fitted) parameters. If \mu is fixed on some grid points (i.e. via fixed_mu), then only the non-fixed mus are assumed to be identical to one another.

Logical, specifying whether \psi should be assumed constant across the grid, i.e. time-independent. Setting const_psi=TRUE reduces the number of free (i.e., independently fitted) parameters. If \psi is fixed on some grid points (i.e. via fixed_psi), then only the non-fixed psis are assumed to be identical to one another.

Logical, specifying whether \kappa should be assumed constant across the grid, i.e. time-independent. Setting const_kappa=TRUE reduces the number of free (i.e., independently fitted) parameters. If \kappa is fixed on some grid points (i.e. via fixed_kappa), then only the non-fixed kappas are assumed to be identical to one another.

Logical, specifying whether the \rho_1, \rho_2, ... should be the same across all CSAs. Setting const_CSA_probs=TRUE reduces the number of free (i.e., independently fitted) parameters. If some of the \rho_1, \rho_2, ... are fixed (i.e. via fixed_CSA_probs), then only the non-fixed CSA_probs are assumed to be identical to one another.

Logical, specifying whether the \kappa_1, \kappa_2, ... should be the same across all CSAs. Setting const_CSA_kappas=TRUE reduces the number of free (i.e., independently fitted) parameters. If some of the \kappa_1, \kappa_2, ... are fixed (i.e. via fixed_CSA_kappas), then only the non-fixed CSA_kappas are assumed to be identical to one another.

Integer between 0 and 3 (inclusive), specifying the polynomial degree of \lambda, \mu, \psi and \kappa between age-grid points. If 0, then \lambda, \mu, \psi and \kappa are considered piecewise constant, if 1 they are considered piecewise linear, if 2 or 3 they are considered to be splines of degree 2 or 3, respectively. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on. A degree of 0 is generally not recommended. The case splines_degree=0 is also known as “skyline” model.

Character, either "crown", "stem", "none" or "auto", specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root. If "stem", the likelihood is conditioned on the survival of the stem lineage. Note that "crown" really only makes sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. "none" is generally not recommended. If "auto", the condition is chosen according to the above recommendations.

Positive unitless number, specifying the default relative time step for the ordinary differential equation solvers. Typical values are 0.01-0.001.

Positive unitless number, specifying the relative difference between subsequent simulated and interpolated values, in internally used ODE solvers. Typical values are 1e-2 to 1e-5. A smaller ODE_relative_dy increases interpolation accuracy, but also increases memory requirements and adds runtime (scaling with the tree's age span, not with Ntips).

Non-negative numeric, in units of time, specfying the age radius around a concentrated sampling attempt, within which to assume that sampling events were due to that concentrated sampling attempt. If NULL, this is chosen automatically based on the anticipated scale of numerical rounding errors. Only relevant if concentrated sampling attempts are included.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of times to attempt finding a valid start point (per trial) before giving up on that trial. Randomly choosen extreme start parameters may occasionally result in Inf/undefined likelihoods, so this option allows the algorithm to keep looking for valid starting points.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Integer, specifying the number of parametric bootstraps to perform for estimating standard errors and confidence intervals of estimated parameters. Set to 0 for no bootstrapping.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Optional list, listing combinations of parameter values of particular interest and for which the log-likelihoods should be returned. Every element of this list should itself be a named list, containing the elements lambda, mu, psi and kappa (each being a numeric vector of size NG) as well as the elements CSA_probs and CSA_kappas (each being a numeric vector of size NCSA). This may be used e.g. for diagnostic purposes, e.g. to examine the shape of the likelihood function.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Logical, specifying whether to print detailed information (such as model likelihoods) at every iteration of the fitting routine. For debugging purposes mainly.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

The log-likelihood of the fitted model for the given timetree.

The log-likelihood of the guessed model for the given timetree.

Named list, specifying the fixed and fitted model parameters. This list will contain the elements lambda, mu, psi and kappa (each being a numeric vector of size NG, listing \lambda,\mu, \psi and \kappa at each age-grid point) as well as the elements CSA_probs and CSA_kappas (each being a numeric vector of size NCSA).

Named list, specifying the guessed model parameters. This list will contain the elements lambda, mu, psi and kappa (each being a numeric vector of size NG) as well as the elements CSA_probs and CSA_kappas (each being a numeric vector of size NCSA). Between grid points \lambda should be interpreted as a piecewise polynomial function (natural spline) of degree splines_degree; to evaluate this function at arbitrary ages use the castor routine evaluate_spline. The same also applies to \mu, \psi and \kappa.

Numeric vector of size NG, the age-grid on which \lambda, \mu, \psi and \kappa are defined. This will be the same as the provided age_grid, unless the latter was NULL or of length <=1.

Numeric vector of size NCSA, ting listhe ages at which concentrated sampling attempts occurred. This is the same as provided to the function.

Integer, number of free (i.e., independently) fitted parameters. If none of the \lambda, \mu and \rho were fixed, and const_lambda=FALSE and const_mu=FALSE, then NFP will be equal to 2*Ngrid+1.

Integer, the number of data points (sampling and branching events) used for fitting.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (number of branching times), and L is the maximized likelihood.

Character, specifying what conditioning was root for the likelihood (e.g. "crown" or "stem").

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

Named list specifying the standard errors of the parameters, based on parametric bootstrapping. This list will contain the elements lambda, mu, psi and kappa (each being a numeric vector of size NG) as well as the elements CSA_probs and CSA_kappas (each being a numeric vector of size NCSA). Only included if Nbootstraps>0. Note that the standard errors of non-fitted (i.e., fixed) parameters will be zero.

Named list specifying the lower end of the 50% confidence interval (i.e. the 25% quantile) for each parameter, based on parametric bootstrapping. This list will contain the elements lambda, mu, psi and kappa (each being a numeric vector of size NG) as well as the elements CSA_probs and CSA_kappas (each being a numeric vector of size NCSA). Only included if Nbootstraps>0.

Similar to CI50lower, but listing the upper end of the 50% confidence interval (i.e. the 75% quantile) for each parameter. For example, the confidence interval for \lambda at age age_grid[1] will be between CI50lower$lambda[1] and CI50upper$lambda[1]. Only included if Nbootstraps>0.

Similar to CI50lower, but listing the lower end of the 95% confidence interval (i.e. the 2.5% quantile) for each parameter. Only included if Nbootstraps>0.

Similar to CI50upper, but listing the upper end of the 95% confidence interval (i.e. the 97.5% quantile) for each parameter. Only included if Nbootstraps>0.

Numeric between 0 and 1, estimated consistency of the data with the fitted model. If L denotes the loglikelihood of new data generated by the fitted model (under the same model) and M denotes the expectation of L, then consistency is the probability that |L-M| will be greater or equal to |X-M|, where X is the loglikelihood of the original data under the fitted model. Only returned if Nbootstraps>0. A low consistency (e.g., <0.05) indicates that the fitted model is a poor description of the data. See Lindholm et al. (2019) for background.

A timetree of class "phylo", representing the time-calibrated reconstructed phylogeny of a set of extant and/or extinct species. Tips of the tree are interpreted as terminally sampled lineages, while monofurcating nodes are interpreted as non-terminally sampled lineages, i.e., lineages sampled at some past time point and with subsequently sampled descendants.

Numeric vector, specifying fixed values for a some or all model parameters. For fitted (i.e., non-fixed) parameters, use NaN or NA. For example, the vector c(1.5,NA,40) specifies that the 1st and 3rd model parameters are fixed at the values 1.5 and 40, respectively, while the 2nd parameter is to be fitted. The length of this vector defines the total number of model parameters. If entries in this vector are named, the names are taken as parameter names. Names should be included if the functions lambda, mu, psi, kappa, CSA_psi and CSA_kappa query parameter values by name (as opposed to numeric index).

Numeric vector of size NP, specifying a first guess for the value of each model parameter. For fixed parameters, guess values are ignored. Can be NULL only if all model parameters are fixed.

Optional numeric vector of size NP, specifying lower bounds for model parameters. If of size 1, the same lower bound is applied to all parameters. Use -Inf to omit a lower bound for a parameter. If NULL, no lower bounds are applied. For fixed parameters, lower bounds are ignored.

Optional numeric vector of size NP, specifying upper bounds for model parameters. If of size 1, the same upper bound is applied to all parameters. Use +Inf to omit an upper bound for a parameter. If NULL, no upper bounds are applied. For fixed parameters, upper bounds are ignored.

Optional numeric vector of size NP, specifying typical scales for model parameters. If of size 1, the same scale is assumed for all parameters. If NULL, scales are determined automatically. For fixed parameters, scales are ignored. It is strongly advised to provide reasonable scales, as this facilitates the numeric optimization algorithm.

Positive numeric, specifying the age of the tree's root. Can be used to define a time offset, e.g. if the last tip was not actually sampled at the present. If NULL, this will be calculated from the tree and it will be assumed that the last tip was sampled at the present.

Strictly positive numeric, specifying the oldest time before present (“age”) to consider when calculating the likelihood. If this is equal to or greater than the root age, then oldest_age is interpreted as the stem age. If oldest_age is less than the root age, the tree is split into multiple subtrees at that age by treating every edge crossing that age as the stem of a subtree, and each subtree is considered an independent realization of the HBDS model stemming at that age. This can be useful for avoiding points in the tree close to the root, where estimation uncertainty is generally higher. If oldest_age==NULL, it is automatically set to the root age.

Function specifying the speciation rate at any given age (time before present) and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with strictly positive entries. Can also be a single numeric (i.e., lambda is fixed).

Function specifying the extinction rate at any given age and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with non-negative entries. Can also be a single numeric (i.e., mu is fixed).

Function specifying the continuous (Poissonian) lineage sampling rate at any given age and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with non-negative entries. Can also be a single numeric (i.e., psi is fixed).

Function specifying the retention probability for continuously sampled lineages, at any given age and for any given parameter values. This function must take exactly two arguments, the 1st one being a numeric vector (one or more ages) and the 2nd one being a numeric vector of size NP (parameter values), and return a numeric vector of the same size as the 1st argument with non-negative entries. The retention probability is the probability of a sampled lineage remaining in the pool of extant lineages. Can also be a single numeric (i.e., kappa is fixed).

Numeric vector, specifying ages at which the lambda, mu, psi and kappa functionals should be evaluated. This age grid must be fine enough to capture the possible variation in \lambda, \mu, \psi and \kappa over time, within the permissible parameter range. Listed ages must be strictly increasing, and must cover at least the full considered age interval (from 0 to oldest_age). Can also be NULL or a vector of size 1, in which case \lambda, \mu, \psi and \kappa are assumed to be time-independent.

Optional numeric vector, listing ages (in ascending order) at which concentrated sampling attempts occurred. If NULL, it is assumed that no concentrated sampling attempts took place and that all tips were sampled according to the continuous sampling rate psi.

Function specifying the sampling probabilities during the various concentrated sampling attempts, depending on parameter values. Hence, for any choice of parameters, CSA_probs must return a numeric vector of the same size as CSA_ages. Can also be a single numeric (i.e., concentrated sampling probability is fixed).

Function specifying the retention probabilities during the various concentrated sampling attempts, depending on parameter values. Hence, for any choice of parameters, CSA_kappas must return a numeric vector of the same size as CSA_ages. Can also be a single numeric (i.e., retention probability during concentrated samplings is fixed).

Character, either "crown", "stem", "none" or "auto", specifying on what to condition the likelihood. If "crown", the likelihood is conditioned on the survival of the two daughter lineages branching off at the root. If "stem", the likelihood is conditioned on the survival of the stem lineage. Note that "crown" really only makes sense when oldest_age is equal to the root age, while "stem" is recommended if oldest_age differs from the root age. "none" is usually not recommended. If "auto", the condition is chosen according to the above recommendations.

Positive unitless number, specifying the default relative time step for the ordinary differential equation solvers. Typical values are 0.01-0.001.

Positive unitless number, specifying the relative difference between subsequent simulated and interpolated values, in internally used ODE solvers. Typical values are 1e-2 to 1e-5. A smaller ODE_relative_dy increases interpolation accuracy, but also increases memory requirements and adds runtime (scaling with the tree's age span, not with Ntips).

Non-negative numeric, in units of time, specfying the age radius around a concentrated sampling attempt, within which to assume that sampling events were due to that concentrated sampling attempt. If NULL, this is chosen automatically based on the anticipated scale of numerical rounding errors. Only relevant if concentrated sampling attempts are included.

Integer, specifying the number of independent fitting trials to perform, each starting from a random choice of model parameters. Increasing Ntrials reduces the risk of reaching a non-global local maximum in the fitting objective.

Integer, specifying the number of times to attempt finding a valid start point (per trial) before giving up on that trial. Randomly chosen extreme start parameters may occasionally result in Inf/undefined likelihoods, so this option allows the algorithm to keep looking for valid starting points.

Integer, specifying the number of parallel threads to use for performing multiple fitting trials simultaneously. This should generally not exceed the number of available CPUs on your machine. Parallel computing is not available on the Windows platform.

Optional numeric, specifying the maximum number of seconds to allow for each evaluation of the likelihood function. Use this to abort fitting trials leading to parameter regions where the likelihood takes a long time to evaluate (these are often unlikely parameter regions).

Integer, specifying the number of parametric bootstraps to perform for estimating standard errors and confidence intervals of estimated model parameters. Set to 0 for no bootstrapping.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Named list containing options for the nlminb optimization routine, such as iter.max, eval.max or rel.tol. For a complete list of options and default values see the documentation of nlminb in the stats package.

Optional numeric matrix having NP columns and an arbitrary number of rows, listing combinations of parameter values of particular interest and for which the log-likelihoods should be returned. This may be used for diagnostic purposes, e.g., to examine the shape of the likelihood function.

Logical, specifying whether to print progress reports and warnings to the screen. Note that errors always cause a return of the function (see return values success and error).

Logical, specifying whether to print detailed information (such as model likelihoods) at every iteration of the fitting routine. For debugging purposes mainly.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether model fitting succeeded. If FALSE, the returned list will include an additional “error” element (character) providing a description of the error; in that case all other return variables may be undefined.

The maximized fitting objective. Currently, only maximum-likelihood estimation is implemented, and hence this will always be the maximized log-likelihood.

The name of the objective that was maximized during fitting. Currently, only maximum-likelihood estimation is implemented, and hence this will always be “loglikelihood”.

The log-likelihood of the fitted model for the given timetree.

Numeric vector of size NP (number of model parameters), listing all fitted or fixed model parameters in their standard order (see details above). If param_names was provided, elements in fitted_params will be named.

Numeric vector of size NP, listing guessed or fixed values for all model parameters in their standard order. If param_names was provided, elements in param_guess will be named.

The loglikelihood of the data for the initial parameter guess (param_guess).

A numeric vector of the same size as nrow(focal_param_values), listing loglikelihoods for each of the focal parameter conbinations listed in focal_loglikelihoods.

Integer, number of fitted (i.e., non-fixed) model parameters.

Number of data points used for fitting, i.e., the number of sampling and branching events that occurred between ages 0 and oldest_age.

The Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of fitted parameters and L is the maximized likelihood.

The Bayesian information criterion for the fitted model, defined as \log(n)k-2\log(L), where k is the number of fitted parameters, n is the number of data points (Ndata), and L is the maximized likelihood.

Character, specifying what conditioning was root for the likelihood (e.g. "crown" or "stem").

Logical, specifying whether the maximum likelihood was reached after convergence of the optimization algorithm. Note that in some cases the maximum likelihood may have been achieved by an optimization path that did not yet converge (in which case it's advisable to increase iter.max and/or eval.max).

Integer, specifying the number of iterations performed during the optimization path that yielded the maximum likelihood.

Integer, specifying the number of likelihood evaluations performed during the optimization path that yielded the maximum likelihood.

Numeric vector of size Ntrials, listing the initial objective values (e.g., loglikelihoods) for each fitting trial, i.e. at the start parameter values.

Numeric vector of size Ntrials, listing the final maximized objective values (e.g., loglikelihoods) for each fitting trial.

Integer vector of size Ntrials, listing the number of start attempts for each fitting trial, until a starting point with valid likelihood was found.

Integer vector of size Ntrials, listing the number of iterations needed for each fitting trial.

Integer vector of size Ntrials, listing the number of likelihood evaluations needed for each fitting trial.

Numeric vector of size NP, estimated standard error of the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, median the estimated parameters across parametric bootstraps. Only returned if Nbootstraps>0.

Numeric vector of size NP, lower bound of the 50% confidence interval (25-75% percentile) for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, upper bound of the 50% confidence interval for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, lower bound of the 95% confidence interval (2.5-97.5% percentile) for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric vector of size NP, upper bound of the 95% confidence interval for the parameters, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric between 0 and 1, estimated consistency of the data with the fitted model (Lindholm et al. 2019). See the documentation of fit_hbds_model_on_grid for an explanation.

Either a single phylogenetic tree of class "phylo", or a list of phylogenetic trees. Edge lengths should correspond (or be analogous) to time. The trees don't need to be ultrametric.

Integer, specifying the number of possible discrete states that the trait can have.

Either an integer vector of size Ntips (only permitted if trees[] is a single tree) or a list containing Ntrees such integer vectors (if trees[] is a list of trees), listing the state of each tip in each tree. Note that tip_states cannot include NAs or NaNs; if the states of some tips are uncertain, you should use the option tip_priors instead. Can also be NULL, in which case tip_priors must be provided.

Either a numeric matrix of size Ntips x Nstates (only permitted if trees[] is a single tree), or a list containing Ntrees such matrixes (if trees[] is a list of trees), listing the likelihood of each state at each tip in each tree. Can also be NULL, in which case tip_states must be provided. Hence, tip_priors[t][i,s] is the likelihood of the observed state of tip i in tree t, if the tip's true state was in state s. For example, if you know for certain that a tip is in state k, then set tip_priors[t][i,s]=1 for s=k and tip_priors[t][i,s]=0 for all other s.

Rate model to be used for the transition rate matrix. Can be "ER" (all rates equal), "SYM" (transition rate i–>j is equal to transition rate j–>i), "ARD" (all rates can be different), "SUEDE" (only stepwise transitions i–>i+1 and i–>i-1 allowed, all 'up' transitions are equal, all 'down' transitions are equal) or "SRD" (only stepwise transitions i–>i+1 and i–>i-1 allowed, and each rate can be different). Can also be an index matrix that maps entries of the transition matrix to the corresponding independent rate parameter to be fitted. Diagonal entries should map to 0, since diagonal entries are not treated as independent rate parameters but are calculated from the remaining entries in the transition rate matrix. All other entries that map to 0 represent a transition rate of zero. The format of this index matrix is similar to the format used by the ace function in the ape package. rate_model is only relevant if transition_matrix==NULL.

Prior probability distribution of the root's states, used to calculate the model's overall likelihood from the root's marginal ancestral state likelihoods. Can be "flat" (all states equal), "empirical" (empirical probability distribution of states across the tree's tips), "stationary" (stationary probability distribution of the transition matrix), "likelihoods" (use the root's state likelihoods as prior), "max_likelihood" (put all weight onto the state with maximum likelihood) or “auto” (will be chosen automatically based on some internal logic). If "stationary" and transition_matrix==NULL, then a transition matrix is first fitted using a flat root prior, and then used to calculate the stationary distribution. root_prior can also be a non-negative numeric vector of size Nstates and with total sum equal to 1.

Optional numeric or numeric vector of size Ntrees, specifying the oldest age (time before present) for each tree to consider when fitting the Mk model. If NULL, the entire trees are considered from the present all the way to their root. If non-NULL, then each tree is “cut” at the corresponding oldest age, yielding multiple subtrees, each of which is assumed to be an independent realization of the Mk process. If oldest_ages is a single numeric, then all trees are cut at the same oldest age. This option may be useful if temporal variation is suspected in the Mk rates, and only data near the present are to be used for fitting to avoid violating the assumptions of a constant-rates Mk model.

Optional 2D numeric matrix, specifying a reasonable first guess for the transition rate matrix. May contain NA. May also be NULL, in which case a reasonable first guess is automatically generated.

Integer, number of trials (starting points) for fitting the transition rate matrix. A higher number may reduce the risk of landing in a local non-global optimum of the likelihood function, but will increase computation time during fitting.

Optional positive integer, number of randomly chosen starting points to consider for all fitting trials except the first one. Among all "scouted" starting points, the Ntrials-1 most promising ones will be considered. A greater number of scouts increases the chances of finding a global likelihood maximum. Each scout costs only one evaluation of the loglikelihood function. If NULL, Nscout is automatically chosen based on the number of fitted parameters and Ntrials. Only relevant if Ntrials>1, since the first trial always uses the default or provided parameter guess.

Optional positive numeric, specifying the maximum time (in seconds) allowed for a single evaluation of the likelihood function. If a specific Mk model takes longer than this threshold to evaluate, then its likelihood is set to -Inf. This option can be used to avoid badly parameterized models during fitting and can thus reduce fitting time. If NULL or <=0, this option is ignored.

Either "optim" or "nlminb", specifying which optimization algorithm to use for maximum-likelihood estimation of the transition matrix.

Maximum number of iterations (per fitting trial) allowed for optimizing the likelihood function.

Relative tolerance (stop criterion) for optimizing the likelihood function.

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Number of parallel threads to use for running multiple fitting trials simultaneously. This only makes sense if your computer has multiple cores/CPUs and if Ntrials>1. This option is ignored on Windows, because Windows does not support forking.

Integer, specifying the number of parametric bootstraps to perform for estimating standard errors and confidence intervals of estimated rate parameters. Set to 0 for no bootstrapping.

Integer, specifying the number of fitting trials to perform for each bootstrap sampling. If NULL, this is set equal to max(1,Ntrials). Decreasing Ntrials_per_bootstrap will reduce computation time, at the expense of potentially inflating the estimated confidence intervals; in some cases (e.g., for very large trees) this may be useful if fitting takes a long time and confidence intervals are very narrow anyway. Only relevant if Nbootstraps>0.

Logical, specifying whether to print progress reports and warnings to the screen.

Logical, specifying whether to print detailed diagnostic messages, mainly for debugging purposes.

Character, specifying the line prefix for printing progress reports to the screen.

Logical, indicating whether the fitting was successful. If FALSE, an additional element error (of type character) is included containing an explanation of the error; in that case the value of any of the other elements is undetermined.

Integer, the number of states assumed for the model.

A matrix of size Nstates x Nstates, the fitted transition rate matrix of the model. The [r,c]-th entry is the transition rate from state r to state c.

Numeric, the log-likelihood of the observed tip states under the fitted model.

Integer, the number of iterations required to reach the maximum log-likelihood. Depending on the optimization algorithm used (see optim_algorithm), this may be NA.

Integer, the number of evaluations of the likelihood function required to reach the maximum log-likelihood. Depending on the optimization algorithm used (see optim_algorithm), this may be NA.

Logical, indicating whether the fitting algorithm converged. Note that fit_Mk may return successfully even if convergence was not achieved; if this happens, the fitted transition matrix may not be reasonable. In that case it is recommended to change the optimization options, for example increasing optim_max_iterations.

Numeric, the initial guess used for the average transition rate, prior to fitting.

Numeric, the Akaike Information Criterion for the fitted model, defined as 2k-2\log(L), where k is the number of independent fitted parameters and L is the maximized likelihood.

Numeric matrix of size Nstates x Nstates, estimated standard error of the fitted transition rates, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric matrix of size Nstates x Nstates, lower bounds of the 50% confidence intervals (25-75% percentile) for the fitted transition rates, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric matrix of size Nstates x Nstates, upper bounds of the 50% confidence intervals for the fitted transition rates, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric matrix of size Nstates x Nstates, lower bounds of the 95% confidence intervals (2.5-97.5% percentile) for the fitted transition rates, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Numeric matrix of size Nstates x Nstates, upper bounds of the 95% confidence intervals for the fitted transition rates, based on parametric bootstrapping. Only returned if Nbootstraps>0.

Phylogenetic tree of class "phylo", representing the evolutionary relationships between sampled species/lineages. Unless Poissonian sampling is included in the model (option sampling_rates), this tree should be ultrametric.

Integer, specifying the number of possible discrete states a tip can have, influencing speciation/extinction rates. For example, if Nstates==2 then this corresponds to the common Binary State Speciation and Extinction (BiSSE) model (Maddison et al., 2007). In the case of a HiSSE/SecSSE model, Nstates refers to the total number of diversification rate categories. For example, in the case of the HiSSE model described by Beaulieu and O'meara (2016), Nstates=4.

Integer, optionally specifying a number of "proxy-states" that are observed instead of the underlying speciation/extinction-modulating states. To fit a HiSSE/SecSSE model, NPstates should be smaller than Nstates. Each state corresponds to a different proxy-state, as defined using the variable proxy_map (see below). For BiSSE/MuSSE with no hidden states, NPstates can be set to either NULL or equal to Nstates; in either case, NPstates will be considered equal to Nstates. For example, in the case of the HiSSE model described by Beaulieu and O'meara (2016), NPstates=2.

Integer vector of size Nstates and with values in 1,..NPstates, specifying the correspondence between states (i.e. diversification-rate categories) and proxy-states, in a HiSSE/SecSSE model. Specifically, proxy_map[s] indicates which proxy-state the state s is represented by. Each proxy-state can represent multiple states (i.e. proxies are ambiguous), but each state must be represented by exactly one proxy-state. For example, to setup the HiSSE model described by Beaulieu and O'meara (2016), use proxy_map=c(1,2,1,2). For non-HiSSE models, set this to NULL or to c(1:Nstates). See below for more details.

Optional character vector of size Nstates, specifying a name/description for each state. This does not influence any of the calculations. It is merely used to add human-readable row/column names (rather than integers) to the returned vectors/matrices. If NULL, no row/column names are added.

Integer vector of size Ntips, listing the proxy state at each tip, in the same order as tips are indexed in the tree. The vector may (but need not) include names; if it does, these are checked for consistency with the tree (if check_input==TRUE). Values must range from 1 to NPstates (which is assumed equal to Nstates in the case of BiSSE/MuSSE). States may also be NA, corresponding to unknown tip proxy states (no information available).

Numeric matrix of size Ntips x Nstates (or of size Ntips x NPstates), listing prior likelihoods of each state (or each proxy-state) at each tip. Can be provided as an alternative to tip_pstates. Thus, tip_priors[i,s] is the likelihood of observing the data (i.e., sampling tip i and observing the observed state) if the tip i was at state s (or proxy-state s). Hence, tip_priors should account for sampling fractions as well as reveal fractions. Either tip_pstates or tip_priors must be non-NULL, but not both.

Numeric vector of size NPstates, with values between 0 and 1, listing the sampling fractions of extant species depending on proxy-state. That is, sampling_fractions[p] is the probability that an extant species, having proxy state p, is included in the phylogeny at present-day. If all extant species are included in the tree with the same probability (i.e., independent of state), this can also be a single number. If NULL (default), all extant species are assumed to be included in the tree. Irrelevant if tip_priors is provided and valid for all tips.

Numeric vector of size NPstates, with values between 0 and 1, listing the probabilities of proxy-state identification depending on proxy-state. That is, reveal_fractions[p] is the probability that a species with proxy-state p will have a known ("revealed") state, conditional upon being included in the tree. This can be used to incorporate reveal biases for tips, depending on their proxy state. Can also be NULL or a single number (in which case reveal fractions are assumed to be independent of proxy-state). Note that only the relative values in reveal_fractions matter, for example c(1,2,1) has the same effect as c(0.5,1,0.5), because reveal_fractions is normalized internally anyway. Irrelevant if tip_priors is provided and valid for all tips.

Numeric vector of size NPstates, listing Poissonian per-lineage sampling rates over time. Hence, sampling_rates[p] is the rate at which lineages are sampled over time when they are in proxy state p. Can also be a single numeric, in which case sampling rates are the same for all proxy states. If NULL, Poissonian sampling is assumed to not occur. Note that earlier MuSSE/HiSSE models (e.g., by Beaulieu and O'Meara, 2016) do not include Poissonian sampling (i.e., all tips are assumed to have been sampled at present-day). Poissonian sampling through time is common in epidemiological models but uncommon in macroevolution models.

Either a character or a 2D integer matrix of size Nstates x Nstates, specifying the model for the transition rates between states. This option controls the parametric complexity of the state transition model, i.e. the number of independent rates and the correspondence between independent and dependent rates. If a character, then it must be one of "ER", "SYM", "ARD", "SUEDE" or "SRD", as used for Mk models (see the function asr_mk_model for details). For example, "ARD" (all rates different) specifies that all transition rates should be considered as independent parameters with potentially different values.

If an integer matrix, then it defines a custom parametric structure for the transition rates, by mapping entries of the transition matrix to a set of independent transition-rate parameters (numbered 1,2, and so on), similarly to the option rate_model in the function asr_mk_model, and as returned for example by the function get_transition_index_matrix. Entries must be between 1 and Nstates, however 0 may also be used to denote a fixed value of zero. For example, if transition_rate_model[1,2]=transition_rate_model[2,1], then the transition rates 1->2 and 2->1 are assumed to be equal. Entries on the diagonal are ignored, since the diagonal elements are always adjusted to ensure a valid Markov transition matrix. To construct a custom matrix with the proper structure, it may be convenient to first generate an "ARD" matrix using get_transition_index_matrix, and then modify individual entries to reduce the number of independent rates.

Either a character or an integer vector of length Nstates, specifying the model for the various birth (speciation) rates. This option controls the parametric complexity of the possible birth rates, i.e. the number of independent birth rates and the correspondence between independent and dependent birth rates. If a character, then it must be either "ER" (equal rates) or "ARD" (all rates different). If an integer vector, it must map each state to an indepedent birth-rate parameter (indexed 1,2,..). For example, the vector c(1,2,1) specifies that the birth-rates \lambda_1 and \lambda_3 must be the same, but \lambda_2 is independent.

Either a character or an integer vector of length Nstates, specifying the model for the various death (extinction) rates. Similar to birth_rate_model.

Either NULL or a 2D matrix of size Nstates x Nstates, specifying known (and thus fixed) transition rates between states. For example, setting some elements to 0 specifies that these transitions cannot occur directly. May also contain NA, indicating rates that are to be fitted. If NULL or empty, all rates are considered unknown and are therefore fitted. Note that, unless transition_rate_model=="ARD", values in transition_matrix are assumed to be consistent with the rate model, that is, rates specified to be equal under the transition rate model are expected to also have equal values in transition_matrix.

Either NULL, or a single number, or a numeric vector of length Nstates, specifying known (and thus fixed) birth rates for each state. May contain NA, indicating rates that are to be fitted. For example, the vector c(5,0,NA) specifies that \lambda_1=5, \lambda_2=0 and that \lambda_3 is to be fitted. If NULL or empty, all birth rates are considered unknown and are therefore fitted. If a single number, all birth rates are considered fixed at that given value.

Either NULL, or a single number, or a numeric vector of length Nstates, specifying known (and thus fixed) death rates for each state. Similar to birth_rates.

Either NULL, or a named list containing optional initial suggestions for various model parameters, i.e. start values for fitting. The list can contain any or all of the following elements:

• transition_matrix: A single number or a 2D numeric matrix of size Nstates x Nstates, specifying suggested start values for the transition rates. May contain NA, indicating rates that should be guessed automatically by the function. If a single number, then that value is used as a start value for all transition rates.

• birth_rates: A single number or a numeric vector of size Nstates, specifying suggested start values for the birth rates. May contain NA, indicating rates that should be guessed automatically by the function (by fitting a simple birth-death model, see fit_tree_model).

• death_rates: A single number or a numeric vector of size Nstates, specifying suggested start values for the death rates. May contain NA, indicating rates that should be guessed automatically by the function (by fitting a simple birth-death model, see fit_tree_model).

Start values are only relevant for fitted (i.e., non-fixed) parameters.

Either NULL or a named list containing optional lower bounds for various model parameters. The list can contain any or all of the elements transition_matrix, birth_rates and death_rates, structured similarly to first_guess. For example, list(transition_matrix=0.1, birth_rates=c(5,NA,NA)) specifies that all transition rates between states must be 0.1 or greater, that the birth rate \lambda_1 must be 5 or greater, and that all other model parameters have unspecified lower bound. For parameters with unspecified lower bounds, zero is used as a lower bound. Lower bounds only apply to fitted (i.e., non-fixed) parameters.

Either NULL or a named list containing optional upper bounds for various model parameters. The list can contain any or all of the elements transition_matrix, birth_rates and death_rates, structured similarly to upper. For example, list(transition_matrix=2, birth_rates=c(10,NA,NA)) specifies that all transition rates between states must be 2 or less, that the birth rate \lambda_1 must be 10 or less, and that all other model parameters have unspecified upper bound. For parameters with unspecified upper bounds, infinity is used as an upper bound. Upper bounds only apply to fitted (i.e., non-fixed) parameters.

Either a character or a numeric vector of size Nstates, specifying the prior probabilities of states for the root, i.e. the weights for obtaining a single model likelihood by averaging the root's state likelihoods. If a character, then it must be one of "flat", "empirical", "likelihoods", "max_likelihood" or "auto". "empirical" means the root's prior is set to the proportions of (estimated) extant species in each state (correcting for sampling fractions and reveal fractions, if applicable). "likelihoods" means that the computed state-likelihoods of the root are used, after normalizing to obtain a probability distribution; this is the approach used in the package hisse::hisse v1.8.9 under the option root.p=NULL, and the approach in the package diversitree::find.mle v0.9-10 under the option root=ROOT.OBS. If "max_likelihood", then the root's prior is set to a Dirac distribution, with full weight given to the maximum-likelihood state at the root (after applying the conditioning). If a numeric vector, root_prior specifies custom probabilities (weights) for each state. Note that if root_conditioning is "madfitz" or "herr_als" (see below), then the prior is set before the conditioning and not updated afterwards for consistency with other R packages.

Character, specifying an optional modification to be applied to the root's state likelihoods prior to averaging. Can be "none" (no modification), "madfitz", "herr_als", "crown" or "stem". "madfitz" and "herr_als" (after van Els, Etiene and Herrera-Alsina 2018) are the options implemented in the package hisse v1.8.9, conditioning the root's state-likelihoods based on the birth-rates and the computed extinction probability (after or before averaging, respectively). See van Els (2018) for a comparison between "madfitz" and "herr_als". The option "stem" conditions the state likelihoods on the probability that the stem lineage would survive until the present. The option "crown" conditions the state likelihoods on the probability that a split occurred at oldest_age and that the two child lineages survived until the present; this option is only recommended if oldest_age is equal to the root age.

Strictly positive numeric, specifying the oldest age (time before present) to consider for fitting. If this is smaller than the tree's root age, then the tree is split into multiple subtrees at oldest_age, and each subtree is considered as an independent realization of the same diversification/evolution process whose parameters are to be estimated. The root_conditioning and root_prior are applied separately to each subtree, prior to calculating the joint (product) likelihood of all subtrees. This option can be used to restrict the fitting to a small (recent) time interval, during which the MuSSE/BiSSE assumptions (e.g., time-independent speciation/extinction/transition rates) are more likely to hold. If oldest_age is NULL, it is automatically set to the root age. In principle oldest_age may also be older than the root age.

Non-negative integer, specifying the number of trials for fitting the model, using alternative (randomized) starting parameters at each trial. A larger Ntrials reduces the risk of landing on a local non-global optimum of the likelihood function, and thus increases the chances of finding the truly best fit. If 0, then no fitting is performed, and only the first-guess (i.e., provided or guessed start params) is evaluated and returned. Hence, setting Ntrials=0 can be used to obtain a reasonable set of start parameters for subsequent fitting or for Markov Chain Monte Carlo.