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.

modelTypes

a character string vector specifying a set (family) of model-classes, to which the constructed model object belongs. These are used for model-selection.

k

integer denoting the number of traits (defaults to 1).

regimes

a character or integer vector denoting the regimes.

params

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.

vecParams

NULL (default) or a numeric vector the vector representation of the variable parameters in the model. See details.

offset

integer offset in vecParams; see Details.

spec

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.

See also

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] 2
PCMNumRegimes(modelBM)
#> [1] 1
PCMNumRegimes(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 #> #>
#> [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