Reconstructs ancestral states for internal nodes of a phylogenetic tree while accounting for heterogeneous sampling rates across states or locations.
Arguments
- phy
An object of class
phylospecifying the phylogenetic tree. The tree must be rooted, binary, and contain thetip.stateattribute. Use thecheck_tree_compatibility()function to verify thatphyis satisfies these conditions.- Q
An \(n\times n\) transition rate matrix, where \(n\) is the number of states. Row and column names must correspond to the states. Off-diagonal elements are the transition rates between states; diagonal elements are set such that rows sum to zero.
- pars
A
data.framewith the other parameters used in the ancestral state reconstruction algorithm. Must have the following columns:state: a vector of state names. These should match the column and row names ofQ.lambda: a vector of birth rates for the birth-death-sampling process.mu: a vector of death rates for the birth-death-sampling process.psi: a vector of sampling rates for the birth-death-sampling process.
It can also have
root_priorthat specifies the a priori distribution of the state of the root; if one is not provided a uniform distribution over the root state will be used.- sensitivity_test
A boolean indicating if a sensitivity test should be conducted (see Details). The default value is
FALSE.
Details
If sensitivity_test=TRUE, then a test of the sensitivity of saasi to the birth-death-sampling rate parameters is conducted.
The values of \(\lambda\), \(\mu\), and \(\psi\) are perturbed by sampling values at random from +/- 10% of the input value.
The ancestral state reconstruction is then run on the perturbed parameters; this is repeated 10 times.
To summarise sensitivity, the modal state for each internal node is first computed for each perturbed reconstruction.
Then the total number of differences in modal states are computed between the perturbed reconstructions and the original reconstruction.
The average number of these total differences across all 10 perturbed reconstructions is reported.
Examples
head(demo_pars) #contains state, lambda, mu, psi for each state
#> state rootprior lambda mu psi
#> 1 1 0.3333333 1.5 0.3 0.1
#> 2 2 0.3333333 1.5 0.3 0.5
#> 3 3 0.3333333 1.5 0.3 0.5
result <- saasi(phy = demo_tree_prepared,
Q = demo_Q,
pars = demo_pars)
#> Tree is compatible with SAASI
head(result)
#> 1 2 3
#> 14 0.693604070 0.2318290 0.0745668922
#> 15 0.263875644 0.6429590 0.0931653929
#> 16 0.001350377 0.9972886 0.0013609914
#> 17 0.000821105 0.9982972 0.0008816567
#> 18 0.004837719 0.9862182 0.0089440393
#> 19 0.004747732 0.9433805 0.0518717196