Likelihood of a multivariate Gaussian phylogenetic comparative model with non-interacting lineages

The likelihood of a PCM represents the probability density function of observed trait values (data) at the tips of a tree given the tree and the model parameters. Seen as a function of the model parameters, the likelihood is used to fit the model to the observed trait data and the phylogenetic tree (which is typically inferred from another sort of data, such as an alignment of genetic sequences for the species at the tips of the tree). The PCMLik function provides a common interface for calculating the (log-)likelihood of different PCMs. Below we denote by N the number of tips, by M the total number of nodes in the tree including tips, internal and root node, and by k - the number of traits.

PCMLik(
  X,
  tree,
  model,
  SE = matrix(0, PCMNumTraits(model), PCMTreeNumTips(tree)),
  metaI = PCMInfo(X = X, tree = tree, model = model, SE = SE, verbose = verbose),
  log = TRUE,
  verbose = FALSE
)

Arguments

X	a k x N numerical matrix with possible NA and NaN entries. For i=1,..., N, the column i of X contains the measured trait values for species i (the tip with integer identifier equal to i in tree). Missing values can be either not-available (NA) or not existing (NaN). These two values are treated differently when calculating likelihoods (see PCMPresentCoordinates).
tree	a phylo object with N tips.
model	an S3 object specifying both, the model type (class, e.g. "OU") as well as the concrete model parameter values at which the likelihood is to be calculated (see also Details).
SE	a k x N matrix specifying the standard error for each measurement in X. Alternatively, a k x k x N cube specifying an upper triangular k x k factor of the variance covariance matrix for the measurement error for each tip i=1, ..., N. When SE is a matrix, the k x k measurement error variance matrix for a tip i is calculated as VE[, , i] <- diag(SE[, i] * SE[, i], nrow = k). When SE is a cube, the way how the measurement variance matrix for a tip i is calculated depends on the runtime option PCMBase.Transpose.Sigma_x as follows: if getOption("PCMBase.Transpose.Sigma_x", FALSE) == FALSE (default): VE[, , i] <- SE[, , i] %*% t(SE[, , i]) if getOption("PCMBase.Transpose.Sigma_x", FALSE) == TRUE: VE[, , i] <- t(SE[, , i]) %*% SE[, , i] Note that the above behavior is consistent with the treatment of the model parameters Sigma_x, Sigmae_x and Sigmaj_x, which are also specified as upper triangular factors. Default: matrix(0.0, PCMNumTraits(model), PCMTreeNumTips(tree)).
metaI	a list returned from a call to PCMInfo(X, tree, model, SE), containing meta-data such as N, M and k. Alternatively, this can be a character string naming a function or a function object that returns such a list, e.g. the functionPCMInfo or the function PCMInfoCpp from the PCMBaseCpp package.
log	logical indicating whether a log-likelehood should be calculated. Default is TRUE.
verbose	logical indicating if some debug-messages should printed.

Value

a numerical value with named attributes as follows:

X0

A numerical vector of length k specifying the value at the root for which the likelihood value was calculated. If the model contains a member called X0, this vector is used; otherwise the value of X0 maximizing the likelihood for the given model parameters is calculated by maximizing the quadratic polynomial 'X0 * L_root * X0 + m_root * X0 + r_root';

error

A character string with information if a numerical or other logical error occurred during likelihood calculation.

If an error occured during likelihood calculation, the default behavior is to return NA with a non-NULL error attribute. This behavior can be changed in using global options:

"PCMBase.Value.NA"

Allows to specify a different NA value such as -Inf or -1e20 which can be used in combination with log = TRUE when using optim to maximize the log-likelihood;

"PCMBase.Errors.As.Warnings"

Setting this option to FALSE will cause any error to result in calling the stop R-base function. If not caught in a tryCatch, this will cause the inference procedure to abort at the occurence of a numerical error. By default, this option is set to TRUE, which means that getOption("PCMBase.Value.NA", as.double(NA)) is returned with an error attribute and a warning is issued.

Details

For efficiency, the argument metaI can be provided explicitly, because this is not supposed to change during a model inference procedure such as likelihood maximization.

See also

PCMInfo PCMAbCdEf PCMLmr PCMSim PCMCond

Examples

N <- 10
tr <- PCMTree(ape::rtree(N))

model <- PCMBaseTestObjects$model_MixedGaussian_ab

PCMTreeSetPartRegimes(tr, c(`11` = 'a'), setPartition = TRUE)

set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
X <- PCMSim(tr, model, X0 = rep(0, 3))

PCMLik(X, tr, model)
#> [1] -67.52985
#> attr(,"X0")
#> [1] 5 2 1