Calculate the variance covariance k x k matrix at time t, under a PCM model

Calculate the variance covariance k x k matrix at time t, under a PCM model

PCMVarAtTime(
  t,
  model,
  W0 = matrix(0, PCMNumTraits(model), PCMNumTraits(model)),
  SE = matrix(0, PCMNumTraits(model), PCMNumTraits(model)),
  regime = PCMRegimes(model)[1L],
  verbose = FALSE
)

Arguments

t	positive numeric denoting time
model	a PCM model object
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 k matrix specifying the upper triangular factor of the measurement error variance-covariance matrix. The product SE Default: SE = matrix(0.0, PCMNumTraits(model), PCMNumTraits(model)).
regime	an integer or a character denoting the regime in model for which to do the calculation; Defaults to PCMRegimes(model)[1L], meaning the first regime in the model.
verbose	a logical indicating if (debug) messages should be written on the console (Defaults to FALSE).

Value

A numeric k x k matrix

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.292899 1.143191 3.287773 -0.9373711 5.00441 1.808664 0.5926121 0.7527575
# 2. Then we load this random vector into the model.
PCMParamLoadOrStore(modelBM, randomParams, 0, PCMNumTraits(modelBM), PCMNumRegimes(modelBM), TRUE)
#> [1] 8

# PCMVarAtTime(1, modelBM)

# note that the variance at time 0 is not the 0 matrix because the model has a non-zero
# environmental deviation
PCMVarAtTime(0, modelBM)
#>           [,1]      [,2]
#> [1,] 3.6224531 0.4460931
#> [2,] 0.4460931 0.5666438