glinv_gauss {glinvci}

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 class 'glinv_gauss'
lik(mod, par = NULL, ...)

## S3 method for class 'glinv_gauss'
grad(mod, par = NULL, lik = FALSE, num_threads = 2L, ...)

## S3 method for class 'glinv_gauss'
hess(
  mod,
  par = NULL,
  lik = FALSE,
  grad = FALSE,
  directions = NULL,
  num_threads = 2L,
  ...
)

## S3 method for class 'glinv_gauss'
print(x, ...)

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.

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.

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

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)