ratematrix (version 1.2)

likelihoodFunction: 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.

Value

The log likelihood for the multivariate Brownian motion model.

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<U+00E9>, 2014) to make calculations for the single regime case.

References

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

Examples

# NOT RUN {
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)
# }