ericgoolsby / Rphylopars

Phylogenetic Comparative Tools for Missing Data and Within-Species Variation
28 stars 11 forks source link

equations for OU and EB models? #50

Closed r03ert0 closed 2 years ago

r03ert0 commented 2 years ago

Hello!

(sorry for the naive question...) What are the equations that Rphylopars implements for the OU and EB models? and how do the parameters of the fit correspond with the equations? For example, for the Ornstein-Uhlenbeck model, I was expecting something like:

$ \frac {d x_t} {dt} = \theta (\mu - x_t) dt + \sigma d W_t$.

In my p_OU fit result, I'd say p_OU$mu corresponds with the \mu value in the equation. Is p_OU$model$alpha = \theta, and \sigma maybe the sqrt of the diagonal of p_OU$model$stationary_cov?

thank you in advance 🙏 !

ericgoolsby commented 2 years ago

Hello! The p_OU$mu value returned by phylopars() is simply the estimated ancestral state of the root (i.e. you'll see the same result with p_OU$anc_recon[nspecies+1,]. Regarding how deviations from BM are handled, the branch lengths are transformed according to the evolutionary model parameter and then multiplied by "phylocov" (the "evolutionary rate" for BM, i.e. the interspecific trait covariance). Following tree transformation, the model is fit as if it were a BM model, and several different values for the evolutionary model parameter are evaluated until the likelihood is maximized. An important note for multivariate traits: the same evolutionary model parameter is assumed to be by all traits (with the exception of the "mvOU" model which allows a separate alpha for each trait and optionally for trait covariances as well).Also, if multiple observations per species are included, you also get a "phenocov" estimate (within-species trait covariance).

Rphylopars transforms the tree branch lengths in the same way as the phylolm package (described here: https://academic.oup.com/sysbio/article/63/3/397/1649891). For example, in the OU model, the branch length transformation is dictated by alpha. More specifically, the distance from the root to the common ancestor of any two tips i and j is either:

assuming the root trait value (mu) is fixed: sigma^2 / (2alpha) (1 - exp(-2alphat_ij)) exp(-alpha d_ij)

or if the root is integrated out as a random variable with the stationary distribution: sigma^2 / (2alpha) exp(-alpha * d_ij)

Here, d_ij equals the distance from taxon i to taxon j, and t_ij is the distance from the root to the common ancestor of i and j.

The EB branch lengths are given by:

(exp(rate_eb t_ij) - exp(rate_eb t_ij_anc)) / rate_eb

Julien Clavel's Methods in Ecology & Evolution paper for mvMORPH (https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12420) has more details about these models, and even more details about the OU and multivariate OU in Appendix 1.