The PCMBase Parametrization API

Venelin Mitov

2019-03-15

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:

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
PCMParamGetShortVector(modelObject)
## [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)
## Warning: package 'data.table' was built under R version 3.5.2
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