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 diffent 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.

So let us create

modelObject <- structure(
  list(X0 = structure(c(0.0, 0.2),
                      class = c("VectorParameter", "_Global", "numeric")),
       Sigma = structure(abind(matrix(c(1,0,0,1), 2, 2),
                               matrix(c(2,0,0,2), 2, 2), along = 3),
                         class = c("MatrixParameter", "_UpperTriangularWithDiagonal", "_WithNonNegativeDiagonal", "matrix"))),
  class = c("PCM"),
  k = 2,
  regimes = 1:2)

PCMRegimes(modelObject)
## [1] 1 2
PCMParamCount(modelObject)
## [1] 8
## [1] 0.0 0.2 1.0 0.0 1.0 2.0 0.0 2.0
vec <- double(PCMParamCount(modelObject))
PCMParamLoadOrStore(modelObject, vec, offset = 0, load=FALSE)
## [1] 8
vec
## [1] 0.0 0.2 1.0 0.0 1.0 2.0 0.0 2.0
modelObjectLowerLimit <- PCMParamLowerLimit(modelObject)
PCMParamLoadOrStore(modelObjectLowerLimit, vec, 0, load=FALSE)
## [1] 8
vec
## [1] -10 -10   0 -10   0   0 -10   0
matParams <- PCMParamRandomVecParams(modelObject, n = 10)

matParamsX0 <- PCMParamRandomVecParams(modelObject$X0, 2, 2, n = 10)

matParamsSigma <- PCMParamRandomVecParams(modelObject$Sigma, 2, 2, n = 10)

Specifying a model type

library(data.table)
OUModelDummy <- list()
class(OUModelDummy) <- c("OU")

listParameterizationsOU <- PCMListParameterizations(OUModelDummy)
dtParameterizations <- PCMTableParameterizations(OUModelDummy)

PCMGenerateParameterizations(OUModelDummy, tableParameterizations = dtParameterizations[1:10])

attr(OUModelDummy, "k") <- 2
attr(OUModelDummy, "regimes") <- letters[1:3]
OUObject <- PCM("OU", k = 2, regimes = letters[1:3])


OUObject2 <- PCM("OU__Omitted_X0__H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigmae_x", k = 2, regimes = letters[1:3])
library(data.table)
OUModelDummy <- list()
class(OUModelDummy) <- c("OU")

listParameterizationsOU <- PCMListParameterizations(OUModelDummy)
dtParameterizations <- PCMTableParameterizations(OUModelDummy)

PCMGenerateParameterizations(OUModelDummy, tableParameterizations = dtParameterizations[1:10])

MGObject <- MixedGaussian(k = 2, modelTypes =  PCMModels(), mapping = c(1, 3, 5))

ul <- PCMParamUpperLimit(MGObject)
rv <- PCMParamRandomVecParams(MGObject)
PCMParamGetShortVector(ul)
##  [1] 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
## [24] 10 10 10 10 10 10 10 10 10 10 10 10 10