glinv_gauss: Construct an object representing a GLInv model with respect to the underlying Gaussian process parameters.

Description

The glinv_gauss function constructs an object of class glinv_gauss, which represents a lower-level GLInv model with respect to the underlying Gaussian process space. The likelihood Hessian of, for example, Brownian motion and Ornstein-Uhlenbeck models can be computed by applying the calculus chain rule to the output of Jacobians and Hessians from this class.

The lik.glinv_gauss function computes the likelihood of a full glinv_gauss model.

The grad.glinv_gauss function computes the log-likelihood gradients of a glinv_gauss models. If par is NULL, it returns a function that, when called, returns the same thing as if grad.glinv_gauss were called with par argument.

The hess.glinv_gauss function computes the log-likelihood Hessian of a glinv_gauss models.

Usage

glinv_gauss(tree, x0, dimtab = NULL, X = NULL)
# S3 method for glinv_gauss
lik(mod, par = NULL, ...)
# S3 method for glinv_gauss
grad(mod, par = NULL, lik = FALSE, num_threads = 2L, ...)
# S3 method for glinv_gauss
hess(
  mod,
  par = NULL,
  lik = FALSE,
  grad = FALSE,
  directions = NULL,
  num_threads = 2L,
  ...
)
# S3 method for glinv_gauss
print(x, ...)

Value

An object of S3 class glinv_gauss with the following components

ctree: A pointer to an internal C structure.
apetree: Same as the tree argument but with some pre-processing in its edge table
origtree: The tree argument.
x0: The trait vector at the root of the tree.
dimtab: Identical to the dimtab argument.
gaussdim: The number of dimension of the parameter space of this model.

Arguments

tree: A tree of class ape::phylo.
x0: A vector representing the root's trait vector.
dimtab: An integer, a vector of integers, or NULL, specifying the number of dimensions of each nodes of the tree. If it is a vector, dimtab[n] is the trait vector dimension of node n. If it is only a single integer than all nodes are assumed to have the same amount of dimensions. If it is NULL then all nodes are asummed to have the same amount of dimensions as x0.
X: Trait values, either a matrix in which X[p,n] stores the $p$-th dimension of the multivariate trait of the $n$-th tip of the phylogeny, or a list in which X[[n]] is a numeric vector representing the multivariate trait of the $n$-th tip. The latter form is required if not all the tips has the same number of dimensions.
mod: A model object of class glinv_gauss.
par: A vector, containing the parameters at which the likelihood should be computed.
...: Not used.
lik: If TRUE, grad.glinv_gauss and hess.glinv_gauss returns also the log-likelihood.
num_threads: Number of threads to use.
grad: If TRUE, hess.glinv_gauss returns also the gradient.
directions: Either NULL or a matrix with mod$nparams many rows and arbitrarily many columns. If NULL, `hess.glinv_gauss` returns the Hessian matrix itself, which is typically a huge matrix; otherwise, the funciton returns a square matrix $M$ such that $M_ij$ contains $d_i^T H d_j$, where $d_i$ is the $i$-th column of directions and $H$ is the huge Hessian matrix, without storing $H$ itself in memory.
x: An object of class glinv_gauss.

Details

The glinv_gauss class does not include any information for dealing with evolutionary regimes, lost traits, and missing data, nor does it facilitate reparametrisation. These are all functionality of the glinv class instead. The member variables of the objects of the glinv_gauss class only are for the users' convenience to read the information about the model, and the user should not modify its member variables directly.

For each non-root node $i$ in the phylogeny, the multivariate trait vector $x_i$ follows a Gaussian distribution with mean $\Phi_i x_j + w_i$ and variance $V_i$ when conditioned on the mother's trait vector $x_j$. The `parameters' of this model is, therefore, the joint of all $(\Phi_i, w_i V'_i)$ for all nodes $i$. The root does not have any associated parameters.

The parameter vector par should be the concatenation of all $(\Phi_i, w_i, V'_i)$ in accending order sorted by $i$, the node number (which is the same node numbers as in tree$edge). The matrix $\Phi_i$ is flattened in column-major order and $V'_i$ is the lower-triangular part of V_i, column-major-flattened. Since the root does not have parameters, its entry is simply skipped. For example, if a binary tree has 10 non-root nodes in total and each of them are 3 dimensional, then each $(\Phi_i, w_i, V'_i)$ should have $9+3+6=18$ elements; thus after concatenation par should be a 180 elements.

References

Mitov V, Bartoszek K, Asimomitis G, Stadler T (2019). “Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts.” Theor. Popul. Biol.. https://doi.org/10.1016/j.tpb.2019.11.005.

Examples

Run this code

tr = ape::rtree(3)
model = glinv_gauss(tr, x0=c(0,0), X=matrix(rnorm(6),2,3))
par = unlist(
 list(
   list('Phi' = c(1,0,0,1), # Parameters for node #1, a tip
        'w'   = c(-1,1),
        'V'   = c(1,0,1)),  # Lower triangular part of a 2D identity matrix
   list('Phi' = c(2,0,0,2), # For node #2, a tip
        'w'   = c(-2,2),
        'V'   = c(2,0,2)),
   list('Phi' = c(3,0,0,3), # For node #3, a tip
        'w'   = c(-3,3),
        'V'   = c(3,0,3)),
   list('Phi' = c(4,0,0,4), # For node #5. Node #4 skipped as it is the root
        'w'   = c(-4,4),
        'V'   = c(4,0,4))
   ))
print(par)
lik(model, par)
grad(model, par)
hess(model, par)

Run the code above in your browser using DataLab