This is the entry-point function for creating model objects within the PCMBase framework representing a single model-type with one or several model-regimes of this type associated with the branches of a tree. For mixed Gaussian phylogenetic models, which enable multiple model-types, use the MixedGaussian function.

PCM(
model,
modelTypes = class(model)[1],
k = 1L,
regimes = 1L,
params = NULL,
vecParams = NULL,
offset = 0L,
spec = NULL,
...
)

## Arguments

model This argument can take one of the following forms: a character vector of the S3-classes of the model object to be created (one model object can have one or more S3-classes, with the class PCM at the origin of the hierarchy); an S3 object which's class inherits from the PCM S3 class. The Details section explains how these two types of input are processed. a character string vector specifying a set (family) of model-classes, to which the constructed model object belongs. These are used for model-selection. integer denoting the number of traits (defaults to 1). a character or integer vector denoting the regimes. NULL (default) or a list of parameter values (scalars, vectors, matrices, or arrays) or sub-models (S3 objects inheriting from the PCM class). See details. NULL (default) or a numeric vector the vector representation of the variable parameters in the model. See details. integer offset in vecParams; see Details. NULL or a list specifying the model parameters (see PCMSpecify). If NULL (default), the generic PCMSpecify is called on the created object of class model. additional parameters intended for use by sub-classes of the PCM class.

## Value

an object of S3 class as defined by the argument model.

## Details

This is an S3 generic. The PCMBase package defines three methods for it:

• PCM.PCM: A default constructor for any object with a class inheriting from "PCM".

• PCM.character: A default PCM constructor from a character string specifying the type of model.

• PCM.default: A default constructor called when no other constructor is found. When called this constructor raises an error message.

MixedGaussian

## 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
#>
#>
# a BM model with two regimes
modelBM.ab <- PCM("BM", k = 2, regimes = c("a", "b"))
modelBM.ab#> Brownian motion model
#> S3 class: BM, GaussianPCM, PCM; k=2; p=14; regimes: a, b. 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):
#> , , a
#>
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#>
#> , , b
#>
#>      [,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):
#> , , a
#>
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#>
#> , , b
#>
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#>
#>
# print a single parameter of the model (in this case, the root value)
modelBM.ab$X0#> [1] 0 0 #> attr(,"class") #> [1] "VectorParameter" "_Global" "numeric" #> attr(,"description") #> [1] "trait values at the root" # assign a value to this parameter (note that the brackets [] are necessary # to preserve the parameter attributes): modelBM.ab$X0[] <- c(5, 2)

PCMNumTraits(modelBM)#> [1] 2PCMNumRegimes(modelBM)#> [1] 1PCMNumRegimes(modelBM.ab)#> [1] 2
# number of numerical parameters in the model
PCMParamCount(modelBM)#> [1] 8
# Get a vector representation of all parameters in the model
PCMParamGetShortVector(modelBM)#> [1] 0 0 0 0 0 0 0 0
# Limits for the model parameters:
lowerLimit <- PCMParamLowerLimit(modelBM)
upperLimit <- PCMParamUpperLimit(modelBM)

# 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]
#> [1,] -8.384997 6.686661 6.007609 -6.855831 0.07399441 4.663935 -0.04445223
#>          [,8]
#> [1,] 2.897672# 2. Then we load this random vector into the model.
PCMParamLoadOrStore(modelBM, randomParams, 0, PCMNumTraits(modelBM), PCMNumRegimes(modelBM), TRUE)#> [1] 8
print(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] -8.384997  6.686661
#> Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
#> , , 1
#>
#>          [,1]        [,2]
#> [1,] 6.007609 -6.85583117
#> [2,] 0.000000  0.07399441
#>
#> Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the non-heritable variance or the variance of the measurement error):
#> , , 1
#>
#>          [,1]        [,2]
#> [1,] 4.663935 -0.04445223
#> [2,] 0.000000  2.89767245
#>
#>
PCMParamGetShortVector(modelBM)#> [1] -8.38499725  6.68666075  6.00760886 -6.85583117  0.07399441  4.66393497
#> [7] -0.04445223  2.89767245
# generate a random phylogenetic tree of 10 tips
tree <- ape::rtree(10)

#simulate the model on the tree
traitValues <- PCMSim(tree, modelBM, X0 = modelBM\$X0)

# calculate the likelihood for the model parameters, given the tree and the trait values
PCMLik(traitValues, tree, modelBM)#> [1] -63.06549
#> attr(,"X0")
#> [1] -8.384997  6.686661
# create a likelihood function for faster processing for this specific model.
# This function is convenient for calling in optim because it recieves and parameter
# vector instead of a model object.
likFun <- PCMCreateLikelihood(traitValues, tree, modelBM)
likFun(randomParams)#> [1] -63.06549
#> attr(,"X0")
#> [1] -8.384997  6.686661