A heuristic-based guess of the optimal parameters of a model

This is an S3 generic function that returns a nxPCMParamCount(o) matrix of parameter vectors. Implementations try to deduce parameter values based on the passed tree and data, that should be suitable for starting points for likelihood optimization or MCMC inference.

GuessInitVecParams(
  o,
  regimes = as.character(PCMRegimes(o)),
  oParent = o,
  accessExpr = "",
  n = 1L,
  argsPCMParamLowerLimit = NULL,
  argsPCMParamUpperLimit = NULL,
  X = NULL,
  tree = NULL,
  SE = NULL,
  tableAnc = NULL,
  varyParams = FALSE,
  returnWithinBoundsOnly = TRUE,
  res = PCMParamRandomVecParams(o = oParent, k = PCMNumTraits(oParent), R =
    PCMNumRegimes(oParent), n = n, argsPCMParamLowerLimit = argsPCMParamLowerLimit,
    argsPCMParamUpperLimit = argsPCMParamUpperLimit),
  ...
)

Arguments

o	a PCM model object.
regimes	a vector of character strings denoting regimes in o, for which parameters must be guessed. This argument has special usage in S3 methods and should be left to its default setting (PCMRegimes(o)).
oParent	a PCM model object. This argument has special usage in S3 methods and should be left to its default setting (o).
accessExpr	a character string. This argument has special usage in S3 methods and should be left to its default setting ("").
n	integer denoting the number of parameter vectors to generate.
argsPCMParamLowerLimit, argsPCMParamUpperLimit	lists of parameters passed to PCMParamLowerLimit and PCMParamUpperLimit respectively. (see the argument returnWithinBoundsOnly).
X	a k x N numerical matrix with possible NA and NaN entries. For i=1,..., N, the column i of X contains the measured trait values for species i (the tip with integer identifier equal to i in tree). Missing values can be either not-available (NA) or not existing (NaN). These two values are treated differently when calculating likelihoods (see PCMPresentCoordinates).
tree	a phylo object with N tips.
SE	a k x N matrix specifying the standard error for each measurement in X. Alternatively, a k x k x N cube specifying an upper triangular k x k factor of the variance covariance matrix for the measurement error for each tip i=1, ..., N. When SE is a matrix, the k x k measurement error variance matrix for a tip i is calculated as VE[, , i] <- diag(SE[, i] * SE[, i], nrow = k). When SE is a cube, the way how the measurement variance matrix for a tip i is calculated depends on the runtime option PCMBase.Transpose.Sigma_x as follows: if getOption("PCMBase.Transpose.Sigma_x", FALSE) == FALSE (default): VE[, , i] <- SE[, , i] %*% t(SE[, , i]) if getOption("PCMBase.Transpose.Sigma_x", FALSE) == TRUE: VE[, , i] <- t(SE[, , i]) %*% SE[, , i] Note that the above behavior is consistent with the treatment of the model parameters Sigma_x, Sigmae_x and Sigmaj_x, which are also specified as upper triangular factors. Default: matrix(0.0, PCMNumTraits(model), PCMTreeNumTips(tree)).
tableAnc	an ancestor table for tree. Default NULL.
varyParams	logical indicating if fixed (FALSE) or varying (TRUE) parameter values should be returned for each of the n vectors. The default value is TRUE.
returnWithinBoundsOnly	logical indication if the returned parameter vectors should be within the lower and upper bound of the model.
res	an n x p matrix where p is the number of parameters in o. Internal use only.
...	additional arguments passed to implementations. Currently an N x N matrix argument named treeVCVMat can be passed which is equal to the phylogenetic covariance matrix for tree.

Value

an n x p matrix where p is the number of parameters in o.

Examples

library(PCMBase)
library(PCMFit)
modelBM.ab <- modelBM.ab.Guess <- PCM(
  PCMDefaultModelTypes()["B"], k = 2, regimes = c("a", "b"))
modelBM.ab$X0[] <- c(5, 2)
# in regime 'a' the traits evolve according to two independent BM
# processes (starting from the global vecto X0).
modelBM.ab$Sigma_x[,, "a"] <- rbind(c(1.6, 0.01),
                                    c(0, 2.4))
# in regime 'b' there is a correlation between the traits
modelBM.ab$Sigma_x[,, "b"] <- rbind(c(2.4, 0.08),
                                    c(0.0, 0.8))
modelOU.ab.Guess <- PCM(
  PCMDefaultModelTypes()["F"], k = 2, regimes = c("a", "b"))

modelMG.ab.Guess <- MixedGaussian(
k = 2, modelTypes = MGPMDefaultModelTypes(),
mapping = c(
  a = unname(MGPMDefaultModelTypes()["B"]),
  b = unname(MGPMDefaultModelTypes()["F"])))


param <- double(PCMParamCount(modelBM.ab))

# load the current model parameters into param
PCMParamLoadOrStore(modelBM.ab, param, offset=0, load=FALSE)
#> [1] 8
print(param)
#> [1] 5.00 2.00 1.60 0.01 2.40 2.40 0.08 0.80


# make results reproducible
set.seed(2, kind = "Mersenne-Twister", normal.kind = "Inversion")

# number of regimes
R <- 2

# number of extant tips
N <- 100

tree.a <- PCMTree(ape::rtree(n=N))
PCMTreeSetLabels(tree.a)
PCMTreeSetPartRegimes(
    tree.a, part.regime = c(`101` = "a"), setPartition = TRUE)

lstDesc <- PCMTreeListDescendants(tree.a)
splitNode <- names(lstDesc)[which(sapply(lstDesc, length) > N/2 &
   sapply(lstDesc, length) < 2*N/3)][1]

tree.ab <- PCMTreeInsertSingletons(
  tree.a, nodes = as.integer(splitNode),
  positions = PCMTreeGetBranchLength(tree.a, as.integer(splitNode))/2)
PCMTreeSetPartRegimes(
  tree.ab,
  part.regime = structure(
  c("a", "b"), names = as.character(c(N+1, splitNode))),
  setPartition = TRUE)

traits <- PCMSim(tree.ab, modelBM.ab, modelBM.ab$X0)

likFunBM <- PCMCreateLikelihood(traits, tree.ab, modelBM.ab.Guess)
likFunOU <- PCMCreateLikelihood(traits, tree.ab, modelOU.ab.Guess)
likFunMG <- PCMCreateLikelihood(traits, tree.ab, modelMG.ab.Guess)

vecGuessBM <-
  GuessInitVecParams(modelBM.ab.Guess, X = traits, tree = tree.ab)
vecGuessOU <-
  GuessInitVecParams(modelOU.ab.Guess, X = traits, tree = tree.ab)
vecGuessMG <-
  GuessInitVecParams(modelMG.ab.Guess, X = traits, tree = tree.ab)

# likelihood at true parameters used to generate the data
print(param)
#> [1] 5.00 2.00 1.60 0.01 2.40 2.40 0.08 0.80
likFunBM(param)
#> [1] -395.0189
#> attr(,"X0")
#> [1] 5 2

# likelihood at guessed parameter values for the BM model
print(vecGuessBM)
#>          [,1]     [,2]     [,3]       [,4]     [,5]     [,6]      [,7]
#> [1,] 5.195898 3.968937 1.501725 -0.1199992 1.878255 2.759574 0.1405427
#>           [,8]
#> [1,] 0.7858626
likFunBM(vecGuessBM)
#> [1] -398.6362
#> attr(,"X0")
#> [1] 5.195898 3.968937

# likelihood at guessed parameter values for the OU model
print(vecGuessOU)
#>          [,1]     [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]    [,11]
#> [1,] 5.195898 3.968937    0    0    0    0    0    0    0     0 3.636564
#>        [,12]    [,13]    [,14]    [,15]     [,16]    [,17]    [,18]      [,19]
#> [1,] 7.50923 4.084623 6.026438 1.484078 0.2442901 2.217484 2.067015 -0.0392634
#>          [,20]
#> [1,] 0.8843471
likFunOU(vecGuessOU)
#> [1] -397.1146
#> attr(,"X0")
#> [1] 5.195898 3.968937

# likelihood at guessed parameter values for the BM model
print(vecGuessMG)
#>          [,1]     [,2]     [,3]       [,4]     [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 5.195898 3.968937 1.227323 -0.3968675 1.710291    0    0    0    0     0
#>      [,11]    [,12]      [,13]     [,14]     [,15] [,16]      [,17]
#> [1,]     0 2.336112 -0.6563992 0.8052728 0.0135856     0 0.01747068
likFunMG(vecGuessMG)
#> [1] -410.007
#> attr(,"X0")
#> [1] 5.195898 3.968937

# likelihood at completely random parameter values
vecRand <- PCMParamRandomVecParams(modelBM.ab.Guess, k = 2, R = 2)
likFunBM(vecRand)
#> [1] -553.4146
#> attr(,"X0")
#> [1]  4.194961 -2.341388