Expected variance-covariance matrix for each couple of tips (i,j)

Expected variance-covariance matrix for each couple of tips (i,j)

PCMVar(
  tree,
  model,
  W0 = matrix(0, PCMNumTraits(model), PCMNumTraits(model)),
  SE = matrix(0, PCMNumTraits(model), PCMTreeNumTips(tree)),
  metaI = PCMInfo(NULL, tree, model, verbose = verbose),
  internal = FALSE,
  diagOnly = FALSE,
  verbose = FALSE
)

Arguments

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).
W0	a numeric matrix denoting the initial k x k variance covariance matrix at the root (default is the k x k zero matrix).
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.
internal	a logical indicating if the per-node variance-covariances matrices for the internal nodes should be returned (see Value). Default FALSE.
diagOnly	a logical indicating if only the variance blocks for the nodes should be calculated. By default this is set to FALSE, meaning that the co-variances are calculated for all couples of nodes.
verbose	logical indicating if some debug-messages should printed.

Value

If internal is FALSE, a (k x N) x (k x N) matrix W, such that k x k block W[((i-1)*k)+(1:k), ((j-1)*k)+(1:k)] equals the expected covariance matrix between tips i and j. Otherwise, a list with an element 'W' as described above and a k x M matrix element 'Wii' containing the per-node variance covariance matrix for each node: The k x k block Wii[, (i-1)*k + (1:k)] represents the variance covariance matrix for node i.

Examples

# a Brownian motion model with one regime
modelBM <- PCM(model = "BM", k = 2)
# print the model
modelBM
#> Brownian motion model
#> S3 class: BM, GaussianPCM, PCM; k=2; p=8; regimes: 1. Parameters/sub-models:
#> X0 (VectorParameter, _Global, numeric; trait values at the root):
#> [1] 0 0
#> Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#> 
#> Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the non-heritable variance or the variance of the measurement error):
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#> 
#>  
# assign the model parameters at random: this will use uniform distribution
# with boundaries specified by PCMParamLowerLimit and PCMParamUpperLimit
# We do this in two steps:
# 1. First we generate a random vector. Note the length of the vector equals PCMParamCount(modelBM)
randomParams <- PCMParamRandomVecParams(modelBM, PCMNumTraits(modelBM), PCMNumRegimes(modelBM))
randomParams
#>          [,1]     [,2]     [,3]      [,4]     [,5]     [,6]      [,7]     [,8]
#> [1,] 2.098666 3.094479 3.531973 -4.594797 9.926841 6.334933 -5.735837 1.293723
# 2. Then we load this random vector into the model.
PCMParamLoadOrStore(modelBM, randomParams, 0, PCMNumTraits(modelBM), PCMNumRegimes(modelBM), TRUE)
#> [1] 8

# create a random tree of 10 tips
tree <- ape::rtree(10)
covMat <- PCMVar(tree, modelBM)