knitr::opts_chunk$set(echo = TRUE)
library(abind)
library(PCMBase)

Introduction

Model parameters are the key object of interest in every phylogenetic comparative method. PCMBase provides a powerful interface for specifying and manipulating model parameters. This interface is based on the S3 object system (see http://adv-r.had.co.nz/S3.html for an excellent introduction by Hadley Wickham).

In PCMBase, every model is an object of an S3 class, such as “OU”, inheriting from the base S3-class “PCM”. A PCM object represents a named list. Each element of that list can be one of the following:

a global parameter shared by all regimes in the model. For example, this would be the case for a non-heritable variance-covariance parameter \(\Sigma_{e}\) if it is assumed to be the same for every observed species (i.e. tip) in the tree.
a local stacked parameter that has a different value for each regime in the model. For example, in an OU model, the selection strength matrix \(H\) has a different value for each of the \(R\) regimes in the model. Therefore, an OU model contains a \(k\times k\times R\) array member called “H”. The element H[,,1] would correspond to regime 1, H[,,2] to regime 2, etc. As a second example, the long-term optimum parameter of an OU process is a \(k\)-dimensional vector \(\vec{\theta}_{r}\) for each regime in the model. In an OU model with \(R\) regimes, this would be represented by a \(k\times R\) array (matrix) called “Theta”, such that the vector Theta[, 1] corresponds to regime 1. Finally, a local scalar parameter (i.e. a number), would be represented by an \(R\)-vector.
a nested PCM object corresponding to a regime. This is the case for a mixed Gaussian phylogenetic model, where different types of Gaussian processes can be acting on different parts of the tree, represented by different regimes.

Generic functions for PCM parameters

The parametrization API provides the following S3 generic functions:

PCMParamLoadOrStore: Loading the parameters of a model from a numerical vector or storing them into a vector.
PCMParamCount: counting the parameters of a model;
PCMParamGetShortVector: Getting all values of the model parameters as a numeric vector;
PCMParamLowerLimit, PCMParamUpperLimit: Upper and lower limit for the parameters of a model;
PCMApplyTransformation: transforming a parameter.

All of the above generic functions can be called both on individual parameters, as well as on PCM objects. When called on a PCM object, the function will be called recursively on each parameter within the object. The behaviour of the above generic functions is specified by the parameter’s S3 class attribute described in the following section.

Parameter S3 class

Parameters can be numeric vectors, matrices or cubes (i.e. 3-dimensional arrays) with a “class” attribute specifying a set of S3 classes. The S3 classes of a parameter control the parameter’s:

Main type

The main type specifies the mathematical character of the parameter, i.e. as what does the parameter enter in the equations defining the model. The main type S3 class is mandatory and can be one of the following: “ScalarParameter”, “VectorParameter” and “MatrixParameter”. Here’s an example:

OU <- PCM("OU", k = 3, regimes = c("a", "b"))
class(OU$X0)

## [1] "VectorParameter" "_Global"         "numeric"

class(OU$H)

## [1] "MatrixParameter"

class(OU$Theta)

## [1] "VectorParameter" "matrix"

class(OU$Sigma_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"

class(OU$Sigmae_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"

Scope

Every parameter can be global, i.e. having the same value for each regime in a model, or local, i.e. having a different value for each regime. A parameter is assumed to be local scope, unless its class attribute includes the class "_Global". For instance, in the above code example, the root trait value X0 has a global scope. Another example would be an OU model with a global selection strength, drift and non-heritable variance matrices:

OU2 <- PCM(
  paste0(
    "OU_", 
    "_Global_X0_",
    "_Global_H_",
    "_Theta_",
    "_UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigma_x_",
    "_UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"), 
  k = 3, regimes = c("a", "b"))
class(OU2$X0)

## [1] "VectorParameter" "_Global"         "numeric"

class(OU2$H)

## [1] "MatrixParameter" "_Global"         "matrix"

class(OU2$Theta)

## [1] "VectorParameter" "matrix"

class(OU2$Sigma_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"     "_Global"                     
## [5] "matrix"

class(OU2$Sigmae_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"     "_Global"                     
## [5] "matrix"

Omission

In some cases, a parameter has to be omitted, i.e. not present in a model. For example, if a mixed Gaussian model has a global parameter \(\Sigma_{e}\), then having a local scope parameter \(\Sigma_{e}\) for each regime represents a non-identifiable case, i.e. infinitely many different combinations of values for the global and local scope \(\Sigma_{e}\) parameters would fit equally well to the data. To prevent this, it is possible to specify that the local parameter \(\Sigma_{e}\) should be omitted from each regime. This is done by including "_Omitted" in the class attribute of the parameter. Here is an example:

BMOU <- MixedGaussian(
  k = 3, 
  modelTypes = c(
    "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x",
    "OU__Omitted_X0__H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"),
  mapping = c(a = 1, b = 2))

class(BMOU$X0)

## [1] "VectorParameter" "_Global"         "numeric"

class(BMOU$Sigmae_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"     "_Global"                     
## [5] "matrix"

# There are no X0 and Sigmae_x parameters for the regimes because they are omitted:
names(BMOU$a)

## [1] "Sigma_x"

names(BMOU$b)

## [1] "H"       "Theta"   "Sigma_x"

Constancy

In some cases, a parameter will have a known and fixed value. Constant classes include: "_Fixed“,”_Ones“,”_Zeros“,”_Identity". For example, in a mixed Gaussian model, we can specify a fixed value for the global parameter Sigmae_x:

BMOU2 <- MixedGaussian(
  k = 3, 
  modelTypes = c(
    "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x",
    "OU__Omitted_X0__H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"),
  mapping = c(a = 1, b = 2), 
  Sigmae_x = structure(0, 
    class = c("MatrixParameter",
              "_UpperTriangularWithDiagonal", "_WithNonNegativeDiagonal", "_Fixed", "_Global"),
    description =
      "A fixed upper triangular factor of the non-phylogenetic variance-covariance"))

BMOU2$Sigmae_x[] <- rbind(
    c(0.5, 0.02, 0),
    c(0, 0.2, 0.01),
    c(0, 0, 0.1))

class(BMOU2$X0)

## [1] "VectorParameter" "_Global"         "numeric"

class(BMOU2$Sigmae_x)

## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"     "_Fixed"                      
## [5] "_Global"                      "matrix"

# There are no X0 and Sigmae_x parameters for the regimes because they are omitted:
names(BMOU2$a)

## [1] "Sigma_x"

names(BMOU2$b)

## [1] "H"       "Theta"   "Sigma_x"

# Notice that the BMOU model created in the previous example has more parameters than BMOU2:
PCMParamCount(BMOU)

## [1] 33

PCMParamCount(BMOU2)

## [1] 27

Transformation

Due to constraints arising from their mathematical formulation or to facilitate the biological interpretation, the parameters of Gaussian models are often subject to restrictions. For example, the \(k\times k\) matrix parameter \(\mathbf{\Sigma}\) of a BM or an OU process should always be a symmetric positive-definite matrix; the selection strength matrix \(\mathbf{H}\) of an OU process is usually restricted to be at least semi-positive-definite (but not necessarily symmetric), in order to limit the interpretation to cases of stabilizing selection rather than (usually unidentifiable) cases of repulsion. Such parameter restrictions can pose a challenge, during a parameter inference procedure, since many of the proposed parameter values may violate some of the imposed restrictions. For most of the \(\mathcal{G}_{LInv}\) models known to us, the common restrictions can be handled using an appropriate parametrization with a transformation applied to the parameters before calculating the functions \(\vec{\omega}\), \(\mathbf{\Phi}\) and \(\mathbf{V}\) and the model likelihood.

By default, the parameters named Sigma_x, Sigmae_x and Sigmaj_x are always transformed using the formula \(\mathbf{\Sigma} = \mathbf{\Sigma}_{x}\mathbf{\Sigma}_{x}^T\) (see also the runtim option “PCMBase.Transpose.Sigma_x” in ?PCMOptions for changing the above default behavior).

For other parameters, to specify that a parameter should be transformed, its class and the class of the encompassing PCM object(s) must include "_Transformable“. A further S3 class for the parameter specifies the type of transformation. Once a parameter is transformed via a call to the generic function PCMApplyTransformation, the”_Transformable" class is removed, i.e. the returned parameter is not transformable. For convenience, the name of the transformed parameter is kept the same as the name of the untransformed one. Currently the following transformations are defined:

_CholeskyFactor

If a matrix parameter \(\mathbf{S}\) has the classes "_CholeskyFactor" and "_Transformable“, it will be transformed as \(\mathbf{S}^{T}\mathbf{S}\). Note that this is analogical but not equivalent to the transformation \(\mathbf{S}\mathbf{S}^T\), which is the default for the parameters Sigma_x (see above). The transformed parameter will have classes”_SemiPositiveDefinite" and "_Transformed":

M <- structure(
  rbind(c(0.2, 0.5, 1.2),
        c(0, 0.1, 0.02),
        c(0, 0, 1.02)),
  class = c("MatrixParameter", "_CholeskyFactor", "_Transformable", "_Global"))

Mtransf <- PCMApplyTransformation(M)
Mtransf

##      [,1]  [,2]   [,3]
## [1,] 0.04 0.100 0.2400
## [2,] 0.10 0.260 0.6020
## [3,] 0.24 0.602 2.4808
## attr(,"class")
## [1] "MatrixParameter"       "_Global"               "_SemiPositiveDefinite"
## [4] "_Transformed"          "matrix"

_Schur

If both, the positive-definiteness and the symmetry of a matrix parameter are optional, e.g. the matrix \(\mathbf{H}\) of an OU process, such restrictions can be imposed through a Schur decomposition as shown previously by (Clavel, Escarguel, and Merceron 2015). Specifically, we define a \(k\times k\)-dimensional matrix \(\mathbf{H}_{S}\) as follows:

the upper triangle of \(\mathbf{H}_{S}\), excluding the diagonal, specifies \(k(k-1)/2\) rotation angles for Givens rotations (Golub and Van Loan 2013) to obtain a \(k\times k\)-dimensional orthoganal matrix \(\mathbf{Q}\);
the lower triangle of \(\mathbf{H}_{S}\) including the diagonal defines a \(k\times k\) triangular matrix \(\mathbf{T}\).

Then, \(\mathbf{H}\) is obtained from \(\mathbf{Q}\) and \(\mathbf{T}\) as follows (Clavel, Escarguel, and Merceron 2015): \[ \mathbf{H}=\mathbf{Q}\mathbf{T}^{T}\mathbf{Q}^{T} \] The matrix \(\mathbf{H}\) calculated in this way has all of its eigenvalues equal to the elements on the diagonal of \(\mathbf{H}_{S}\) (Clavel, Escarguel, and Merceron 2015). Thus, by restricting the diagonal of \(\mathbf{H}_{S}\) to non-negative or strictly positive values, we guarantee that \(\mathbf{H}\) will have all of its eigenvalues non-negative or strictly positive (hence, guaranteeing semi- or strict positive definiteness). Further, if \(\mathbf{H}_{S}\) is diagonal, then so is the matrix \(\mathbf{H}\); if \(\mathbf{H}_{S}\) is upper triangular, then \(\mathbf{T}\) is diagonal and the resulting matrix \(\mathbf{H}\) is symmetric. Finally, if \(\mathbf{H}_{S}\) is neither diagonal nor triangular, then the resulting matrix \(\mathbf{H}\) is asymmetric. Here is an example:

# Diagonal
Hs1 <- structure(
  rbind(c(0.2, 0, 0),
        c(0, 0, 0),
        c(0, 0, 1.02)),
  class = c("MatrixParameter", "_Schur", "_Transformable", "_Global"))

PCMApplyTransformation(Hs1)

##      [,1] [,2] [,3]
## [1,]  0.2    0 0.00
## [2,]  0.0    0 0.00
## [3,]  0.0    0 1.02
## attr(,"class")
## [1] "MatrixParameter" "_Global"         "_Transformed"    "matrix"

# Symmetric positive definite with eigenvalues 1.02, 0.1 and 0.02
Hs2 <- structure(
  rbind(c(1.02, 0.5, 1.2),
        c(0, 0.1, 0.02),
        c(0, 0, 0.02)),
  class = c("MatrixParameter", "_Schur", "_Transformable", "_Global"))

PCMApplyTransformation(Hs2)

##             [,1]        [,2]       [,3]
## [1,]  0.34835909 -0.08866649 -0.4527552
## [2,] -0.08866649  0.13163579  0.1593944
## [3,] -0.45275522  0.15939438  0.6600051
## attr(,"class")
## [1] "MatrixParameter" "_Global"         "_Transformed"    "matrix"

eigen(PCMApplyTransformation(Hs2))$values

## [1] 1.02 0.10 0.02

# Asymmetric positive definite with eigenvalues 1.02, 0.1 and 0.02
Hs3 <- structure(
  rbind(c(1.02, 0.5, 1.2),
        c(0.2, 0.1, 0.02),
        c(0.8, 0.1, 0.02)),
  class = c("MatrixParameter", "_Schur", "_Transformable", "_Global"))

PCMApplyTransformation(Hs3)

##            [,1]       [,2]       [,3]
## [1,]  0.7527670 -0.1047998 -0.1627623
## [2,] -0.1511891  0.1103937  0.1204947
## [3,] -0.9906014  0.1707989  0.2768393
## attr(,"class")
## [1] "MatrixParameter" "_Global"         "_Transformed"    "matrix"

eigen(PCMApplyTransformation(Hs3))$values

## [1] 1.02 0.10 0.02

Other restrictions

We may need to specify other restrictions for a parameter. For example, by definition the matrix parameters Sigma_x and Sigmae_x are required to be upper triangular with non-negative diagonal elements (Mitov et al. 2019). Restriction classes include: ’_AllEqual’, ’_NonNegative’, ’_WithNonNegativeDiagonal’, ’_LowerTriangular’, ’_LowerTriangularWithDiagonal’, ’_UpperTriangular’, ’_UpperTriangularWithDiagonal’, ’_Symmetric’ and ’_ScalarDiagonal’.

Generating model parametrizations

It can be helpful to examine the structure of a 3-regime bivariate BM model:

modelObject <- PCM("BM", k = 2L, regimes = c("a", "b", "c"))

# let's assign some values to the model parameters:
vec <- seq_len(PCMParamCount(modelObject))
PCMParamLoadOrStore(modelObject, vec, offset = 0, load=TRUE)

## [1] 20

str(modelObject)

## List of 3
##  $ X0      : 'VectorParameter' num [1:2] 1 2
##   ..- attr(*, "description")= chr "trait values at the root"
##  $ Sigma_x : 'MatrixParameter' num [1:2, 1:2, 1:3] 3 0 4 5 6 0 7 8 9 0 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "a" "b" "c"
##   ..- attr(*, "description")= chr "factor of the unit-time variance rate"
##  $ Sigmae_x: 'MatrixParameter' num [1:2, 1:2, 1:3] 12 0 13 14 15 0 16 17 18 0 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "a" "b" "c"
##   ..- attr(*, "description")= chr "factor of the non-heritable variance or the variance of the measurement error"
##  - attr(*, "class")= chr [1:3] "BM" "GaussianPCM" "PCM"
##  - attr(*, "k")= int 2
##  - attr(*, "regimes")= chr [1:3] "a" "b" "c"
##  - attr(*, "spec")=List of 3
##   ..$ X0      : 'VectorParameter' num 0
##   .. ..- attr(*, "description")= chr "trait values at the root"
##   ..$ Sigma_x : 'MatrixParameter' num 0
##   .. ..- attr(*, "description")= chr "factor of the unit-time variance rate"
##   ..$ Sigmae_x: 'MatrixParameter' num 0
##   .. ..- attr(*, "description")= chr "factor of the non-heritable variance or the variance of the measurement error"
##   ..- attr(*, "class")= chr [1:3] "BM" "GaussianPCM" "PCM"
##   ..- attr(*, "k")= int 2
##   ..- attr(*, "regimes")= chr [1:3] "a" "b" "c"
##  - attr(*, "p")= int 20
##  - attr(*, "modelTypes")= chr "BM"
##  - attr(*, "mapping")= int 1

Now let’s print the model parameters as a table:

options(digits = 1)
print(
  PCMTable(modelObject, addTransformed = FALSE, removeUntransformed = FALSE), 
  xtable = TRUE, type='html')

regime	\(X0\)	\(\Sigma_{x}\)	\(\Sigma_{e,x}\)
:global:	\(\begin{bmatrix}{} 1.0 \\ 2.0 \\ \end{bmatrix}\)
a		\(\begin{bmatrix}{} 3.0 & 4.0 \\ 0.0 & 5.0 \\ \end{bmatrix}\)	\(\begin{bmatrix}{} 12.0 & 13.0 \\ 0.0 & 14.0 \\ \end{bmatrix}\)
b		\(\begin{bmatrix}{} 6.0 & 7.0 \\ 0.0 & 8.0 \\ \end{bmatrix}\)	\(\begin{bmatrix}{} 15.0 & 16.0 \\ 0.0 & 17.0 \\ \end{bmatrix}\)
c		\(\begin{bmatrix}{} 9.0 & 10.0 \\ 0.0 & 11.0 \\ \end{bmatrix}\)	\(\begin{bmatrix}{} 18.0 & 19.0 \\ 0.0 & 20.0 \\ \end{bmatrix}\)

We notice that the parameter \(\Sigma_{e,x}\) is local for each regime in the model. In the context of a biological application, it can be more appropriate to assume that this parameter is shared by all regimes. Moreover, in many cases, e.g. when testing hypotheses about the presence of correlation between the traits, it can be necessary to assume that the matrices \(\Sigma_{x}\) and \(\Sigma_{e,x}\) are diagonal, meaning that the model assumes trait independence. Biologically, that would mean that any correlation between the traits originates from genetic (not environmental) factors. To specify such a scenario, first, we look at the list of available BM model types.

PCMModels(parentClass = "BM")

##  [1] "BM"                                                                                                                                              
##  [2] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                       
##  [3] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                    
##  [4] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                           
##  [5] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                 
##  [6] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"              
##  [7] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                     
##  [8] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                    
##  [9] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x" 
## [10] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"        
## [11] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                      
## [12] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                   
## [13] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                          
## [14] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                
## [15] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"             
## [16] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                    
## [17] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                   
## [18] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"
## [19] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"

We notice that the above list does not contain a BM model type with a global diagonal parameter Sigmae_x. This means that the S3 methods for such a model parametrization have not yet been generated and loaded. The function PCMListParameterizations returns a named list, where each named element is a list of the possible S3 class assignments for the corresponding parameter. So, let’s see what are the available S3 classes for the parameter Sigmae_x of a BM model:

PCMListParameterizations(structure(0.0, class="BM"))$Sigmae_x

## [[1]]
## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"    
## 
## [[2]]
## [1] "MatrixParameter"          "_Diagonal"               
## [3] "_WithNonNegativeDiagonal"
## 
## [[3]]
## [1] "MatrixParameter"          "_ScalarDiagonal"         
## [3] "_WithNonNegativeDiagonal"
## 
## [[4]]
## [1] "MatrixParameter"              "_UpperTriangularWithDiagonal"
## [3] "_WithNonNegativeDiagonal"     "_Global"                     
## 
## [[5]]
## [1] "MatrixParameter"          "_Diagonal"               
## [3] "_WithNonNegativeDiagonal" "_Global"                 
## 
## [[6]]
## [1] "MatrixParameter"          "_ScalarDiagonal"         
## [3] "_WithNonNegativeDiagonal" "_Global"                 
## 
## [[7]]
## [1] "MatrixParameter" "_Omitted"

We notice that, for Sigmae_x, the 5th element above contains "_Diagonal" and "_Global" in the list of S3-classes. We proceed to generate all BM parametrizations with this S3-class assignment for the parameter Sigmae_x.

# 1. Filter the list of parametrizations to avoid generating too many S3 methods.
# (note that we could do the same type of filtering for the other parameters).
listParameterizationsBM <- PCMListParameterizations(structure(0.0, class="BM"))
listParameterizationsBM$Sigmae_x <- listParameterizationsBM$Sigmae_x[5]
  
# 2. Generate a table of parametrizations for this list:
dtParameterizations <- PCMTableParameterizations(
  structure(0.0, class="BM"), listParameterizations = listParameterizationsBM)

print(dtParameterizations)

##                                    X0
##  1:           VectorParameter,_Global
##  2:    VectorParameter,_Fixed,_Global
##  3: VectorParameter,_AllEqual,_Global
##  4:          VectorParameter,_Omitted
##  5:           VectorParameter,_Global
##  6:    VectorParameter,_Fixed,_Global
##  7: VectorParameter,_AllEqual,_Global
##  8:          VectorParameter,_Omitted
##  9:           VectorParameter,_Global
## 10:    VectorParameter,_Fixed,_Global
## 11: VectorParameter,_AllEqual,_Global
## 12:          VectorParameter,_Omitted
##                                                                   Sigma_x
##  1: MatrixParameter,_UpperTriangularWithDiagonal,_WithNonNegativeDiagonal
##  2: MatrixParameter,_UpperTriangularWithDiagonal,_WithNonNegativeDiagonal
##  3: MatrixParameter,_UpperTriangularWithDiagonal,_WithNonNegativeDiagonal
##  4: MatrixParameter,_UpperTriangularWithDiagonal,_WithNonNegativeDiagonal
##  5:                    MatrixParameter,_Diagonal,_WithNonNegativeDiagonal
##  6:                    MatrixParameter,_Diagonal,_WithNonNegativeDiagonal
##  7:                    MatrixParameter,_Diagonal,_WithNonNegativeDiagonal
##  8:                    MatrixParameter,_Diagonal,_WithNonNegativeDiagonal
##  9:              MatrixParameter,_ScalarDiagonal,_WithNonNegativeDiagonal
## 10:              MatrixParameter,_ScalarDiagonal,_WithNonNegativeDiagonal
## 11:              MatrixParameter,_ScalarDiagonal,_WithNonNegativeDiagonal
## 12:              MatrixParameter,_ScalarDiagonal,_WithNonNegativeDiagonal
##                                                       Sigmae_x
##  1: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  2: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  3: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  4: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  5: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  6: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  7: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  8: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
##  9: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
## 10: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
## 11: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global
## 12: MatrixParameter,_Diagonal,_WithNonNegativeDiagonal,_Global

# 3. Generate the parametrizations (optionally, we could select a subset of the
# rows in the data.table)
PCMGenerateParameterizations(structure(0.0, class="BM"), 
                             tableParameterizations = dtParameterizations[])

Now the list of BM-models contains the generated model types with a diagonal Sigma_x and a global diagonal Sigmae_x:

PCMModels("BM")

##  [1] "BM"                                                                                                                                              
##  [2] "BM__AllEqual_Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                              
##  [3] "BM__AllEqual_Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                        
##  [4] "BM__AllEqual_Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"           
##  [5] "BM__Fixed_Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                                 
##  [6] "BM__Fixed_Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                           
##  [7] "BM__Fixed_Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"              
##  [8] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                                       
##  [9] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                       
## [10] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                    
## [11] "BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                           
## [12] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                                 
## [13] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                 
## [14] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"              
## [15] "BM__Global_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                     
## [16] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                    
## [17] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                    
## [18] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x" 
## [19] "BM__Global_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"        
## [20] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                                      
## [21] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                      
## [22] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                   
## [23] "BM__Omitted_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                          
## [24] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                                
## [25] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                                
## [26] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"             
## [27] "BM__Omitted_X0__ScalarDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"                    
## [28] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x"                   
## [29] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"                                                   
## [30] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Global_Sigmae_x"
## [31] "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x"       
## [32] "BM_drift"

BMModelGlobalSigmae <- paste0(
  "BM" , 
  "__Global_X0", 
  "__Diagonal_WithNonNegativeDiagonal_Sigma_x", 
  "__Diagonal_WithNonNegativeDiagonal_Global_Sigmae_x")

modelObject2 <- PCM(BMModelGlobalSigmae, k = 2, regimes = c("a", "b", "c"))
vec <- seq_len(PCMParamCount(modelObject2))
PCMParamLoadOrStore(modelObject2, vec, offset = 0, load=TRUE)

[1] 10

print(
  PCMTable(modelObject2, addTransformed = FALSE, removeUntransformed = FALSE), 
  xtable = TRUE, type='html')

regime	\(X0\)	\(\Sigma_{e,x}\)	\(\Sigma_{x}\)
:global:	\(\begin{bmatrix}{} 1.0 \\ 2.0 \\ \end{bmatrix}\)	\(\begin{bmatrix}{} 9.0 & 0.0 \\ 0.0 & 10.0 \\ \end{bmatrix}\)
a			\(\begin{bmatrix}{} 3.0 & 0.0 \\ 0.0 & 4.0 \\ \end{bmatrix}\)
b			\(\begin{bmatrix}{} 5.0 & 0.0 \\ 0.0 & 6.0 \\ \end{bmatrix}\)
c			\(\begin{bmatrix}{} 7.0 & 0.0 \\ 0.0 & 8.0 \\ \end{bmatrix}\)

References

Clavel, Julien, Gilles Escarguel, and Gildas Merceron. 2015. “mvmorph: an R package for fitting multivariate evolutionary models to morphometric data.” Methods in Ecology and Evolution 6 (11): 1311–9.

Golub, G H, and C F Van Loan. 2013. Matrix Computations. Baltimore: The Johns Hopkins University Press.

Mitov, Venelin, Krzysztof Bartoszek, Georgios Asimomitis, and Tanja Stadler. 2019. “Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts.” Theor. Popul. Biol., December. https://doi.org/10.1016/j.tpb.2019.11.005.

The PCMBase Parametrization API

Venelin Mitov

2021-06-07