get.post.traits: Sample posterior trait values at nodes given a fitted...

In bstaggmartin/evorates: Evolving Rate Models of Continuous Trait Evolution

Sample posterior trait values at nodes given a fitted evorates model

Description

This function samples (or extracts previously-sampled) trait values at tips or internal nodes (a form of stochastic ancestral state reconstruction) given an evorates_fit object.

Usage

get.post.traits(
  fit,
  select = seq_len(Ntip(fit) + Nnode(fit)),
  type = c("chains", "quantiles", "means", "diagnostics"),
  extra.select = NULL,
  simplify = TRUE,
  store.in.fit = FALSE,
  force.resample = FALSE
)

Arguments

fit	An object of class "evorates_fit".
select	A numeric vector specifying node indices in fit$call$tree for which to extract trait values. These can be negative to exclude nodes as well. As with all objects of class "phylo", the first n nodes correspond to tips (where n is the number of tips in the tree) and the rest correspond to internal nodes. All nodes (tip and internal alike) are extracted by default.
type	A string specifying whether to extract posterior samples ("chains"), quantiles ("quantiles"), means ("means"), or diagnostics ("diagnostics"). Can be any unambiguous abbreviation of these strings as well. Defaults to extracting posterior samples.
extra.select	A numeric, integer, or character vector specifying the specific samples/ quantiles/diagnostics to extract, depending on type. Defaults to NULL, which generally extracts all available samples/quantiles/diagnostics. See documentation on param_block operators for details (grapes-chains-grapes, grapes-quantiles-grapes, grapes-means-grapes, or grapes-diagnostics-grapes).
simplify	TRUE or FALSE: should the resulting param_block array be simplified? If TRUE (the default), dimensions of length 1 in the result are automatically collapsed, with corresponding information stored as attributes (this is the default behavior of param_block operators).
store.in.fit	(new in 0.1.1) TRUE or FALSE: should the resulting param_block array be added to the param_block arrays in fit? Often makes the fit object quite a bit larger, but makes later computation of trait values faster and more reproducible. Defaults to FALSE. Note that select, type, and extra.select are ignored if this is TRUE.
force.resample	(new in 0.1.1) TRUE or FALSE: should trait values be resampled even if already stored in fit? Defaults to FALSE.

Details

The tip and node indices of fit$call$tree can be viewed by running: plot(fit$call$tree); tiplabels(); nodelabels().

Note that this function does random sampling internally, so you will get different results every time you run it unless you use a consistent seed (e.g., using set.seed()) or store sampled trait values.

Value

If store.in.fit is FALSE: an array of class "param_block" with a param_type set to whatever type is. The dimension of these arrays will generally go in the order of iterations/quantiles/diagnostics, then parameters, then chains. Any dimensions of length 1 are collapsed and stored as attributes if simplify is TRUE. If store.in.fit is TRUE: an evorates_fit object including extra parameters corresponding to trait values. Parameters go by the name <trait.name>_i, where i is the label or index of the corresponding tip or internal node, respectively, in fit$call$tree. The trait name corresponds to the column name of fit$call$trait.data, which is "X1" if not provided in the original input to fit.evorates() or input.evorates().

See Also

Other parameter extraction functions: get.R(), get.bg.rate(), remove.trend()

Examples

#get whale/dolphin evorates fit
data("cet_fit")

#get posterior samples of all trait values
x <- get.post.traits(cet_fit)
#let's say we want to get some trait values towards the root
plot(cet_fit$call$tree); nodelabels()
#select particular edges
x <- get.post.traits(cet_fit, select = c(89, 90, 94, 103, 104, 105))
#maybe specific samples too?
x <- get.post.traits(cet_fit, select = c(89, 90, 94, 103, 104, 105),
                     extra.select=c(1, 23, 47))

#this may be a more common way to get node indices; say we want to look at the genus Mesoplodon
edges <- get.clade.edges(cet_fit$call$tree, "Mesoplodon")
nodes <- unique(as.vector(cet_fit$call$tree$edge[edges,]))
Meso.x <- get.post.traits(cet_fit, select = nodes)
#or at everything EXCEPT Mesoplodon
notMeso.x <- get.post.traits(cet_fit, select = -nodes)
#you could also rely on string matching if you're interested in the tips only
Meso.tips.x <- get.post.traits(cet_fit, select = "Mesoplodon")

#could also look at quantiles, means, or diagnostics
med.x <- get.post.traits(cet_fit, select = nodes,
                     type = "quantiles",
                     extra.select = 0.5)
mean.x <- get.post.traits(cet_fit, select = nodes,
                          type = "means")
init.x <- get.post.traits(cet_fit, select = nodes,
                          type = "diagnostics",
                          extra.select = c("inits", "ess"))
                 
#here's an example of what happens when you don't simplify the result
med.x <- get.R(cet_fit, select = nodes,
               type = "quantiles",
               extra.select = 0.5,
               simplify = FALSE)
                 
#note that, depending on your goals, it may be more efficient to do things like this instead
x <- get.post.traits(cet_fit)
med.Meso.x <- x %quantiles% list(nodes, 0.5)

bstaggmartin/evorates documentation built on May 31, 2024, 5:56 a.m.