PCM.RdThis 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, ... )
| model | This argument can take one of the following forms:
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
|
| ... | additional parameters intended for use by sub-classes of the PCM class. |
an object of S3 class as defined by the argument model.
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.
# 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 #> #>#> 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#> [1] 8#> [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] 8print(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