likelihoodFunction: Likelihood function for the multivariate Brownian motion...

In Caetanods/ratematrix: Bayesian Estimation of the Evolutionary Rate Matrix

Likelihood function for the multivariate Brownian motion model

Description

Returns the log-likelihood for the multivariate Brownian motion model with 1 or more rate regimes mapped to the tree.

Usage

likelihoodFunction(data, phy, root, R)

Arguments

data	a matrix with the data. Species names need to be provided as rownames (rownames(data) == phy$tip.label).
phy	a phylogeny of the class "simmap" with the mapped regimes or "phylo" for a single rate model.
root	a numeric vector with the root value (phylogenetic mean).
R	a matrix or a list of matrices. If 'R' is a matrix then the likelihood for a single regime is calculated. If 'R' is a list of matrices, then each matrix will be fitted to a regime in 'phy' and the length of the list need to match the number of regimes fitted to the tree.

Details

If two or more rate regimes are mapped to the phylogenetic tree, then the function calculates the likelihood using the new prunning algorithm adapted to fit multiple rate regimes. The prunning algorithm is implemented in C++ using 'Rcpp' and 'RcppArmadillo'. Otherwise the function uses the three point algorithm (Ho and Ané, 2014) to make calculations for the single regime case.

Value

The log likelihood for the multivariate Brownian motion model.

Author(s)

Daniel S. Caetano and Luke J. Harmon

References

Ho, L. S. T. and Ané, C. (2014). "A linear-time algorithm for Gaussian and non-Gaussian trait evolution models". Systematic Biology *63*(3):397-408.

Examples

data( centrarchidae )
root <- colMeans( centrarchidae$data )
Rlist <- list( rbind(c(0.5, 0.1),c(0.1,0.5)), rbind(c(0.5, 0),c(0,0.5)) )
likelihoodFunction(data = centrarchidae$data, phy = centrarchidae$phy.map, root = root
                   , R = Rlist)
## Get the likelihood for a single regime model:
phy.single <- mergeSimmap(phy = centrarchidae$phy.map, drop.regimes = TRUE)
Rsingle <- rbind(c(0.5, 0.1),c(0.1,0.5))
likelihoodFunction(data = centrarchidae$data, phy = phy.single, root = root, R = Rsingle)

Caetanods/ratematrix documentation built on Jan. 4, 2023, 3:12 p.m.