Expected variance-covariance matrix for each couple of tips (i,j)
PCMVar(tree, model, W0 = matrix(0, PCMNumTraits(model), PCMNumTraits(model)), metaI = PCMInfo(NULL, tree, model, verbose), internal = FALSE, verbose = FALSE)
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). |
metaI | a list returned from a call to |
internal | a logical indicating if the per-node variance-covariances matrices for the internal nodes should be returned (see Value). Default FALSE. |
verbose | logical indicating if some debug-messages should printed. |
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.
# 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; Choleski factor of the unit-time variance rate): #> , , 1 #> #> [,1] [,2] #> [1,] 0 0 #> [2,] 0 0 #> #> Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; Choleski 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] 1.627321 5.460170 9.951230 4.219425 2.149426 2.917576 4.435195 8.666157# 2. Then we load this random vector into the model. PCMParamLoadOrStore(modelBM, randomParams, 0, PCMNumTraits(modelBM), PCMNumRegimes(modelBM), TRUE)#> [1] 8