Expected mean vector at each tip conditioned on a trait-value vector at the root

Expected mean vector at each tip conditioned on a trait-value vector at the root

PCMMean(
  tree,
  model,
  X0 = model$X0,
  metaI = PCMInfo(NULL, tree, model, verbose = verbose),
  internal = 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).
X0

a k-vector denoting the root trait
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 ig the per-node mean vectors should be returned (see Value). Default FALSE.
verbose

logical indicating if some debug-messages should printed.

Value

If internal is FALSE (default), then a k x N matrix Mu, such that Mu[, i] equals the expected mean k-vector at tip i, conditioned on X0 and the tree. Otherwise, a k x M matrix Mu containing the mean vector for each node.

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,] 8.481489 1.975219 9.761707 4.63585 3.567269 4.314737 -7.035769 0.1307758
# 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)
PCMMean(tree, modelBM)
#>          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#> [1,] 8.481489 8.481489 8.481489 8.481489 8.481489 8.481489 8.481489 8.481489
#> [2,] 1.975219 1.975219 1.975219 1.975219 1.975219 1.975219 1.975219 1.975219
#>          [,9]    [,10]
#> [1,] 8.481489 8.481489
#> [2,] 1.975219 1.975219