Compute the log density of the vector of trait at the tips of the phylogenetic tree, assuming a Cauchy process.
logDensityTipsCauchy(
tree,
tipTrait,
root.value = NULL,
disp,
method = c("reml", "random.root", "fixed.root"),
rootTip = NULL,
do_checks = TRUE
)
a phylogenetic tree of class phylo
.
a names vector of tip trait values, with names matching the tree labels.
the root starting value of the process.
the dispersion value.
the method used to compute the likelihood.
One of reml
(the default), fixed.root
or random.root
.
See Details.
the tip used to re-root the tree, when the REML method is used.
If NULL
, the tip with the smallest distance to all the others is used (see Details).
Ignored in method != "reml"
.
if FALSE
, the entry parameters are not checked for consistency.
This can be useful when doing multiple calls to the function, as in numerical optimization.
Default to TRUE
.
the log density value.
The parameters of the Cauchy Process (CP)
are disp
, the dispersion of the process,
and root.value
, the starting value of the process at the root (for method="fixed.root"
).
The model assumes that each increment of the trait \(X\) on a branch going from node \(k\) to \(l\) follows a Cauchy distribution, with a dispersion proportional to the length \(t_l\) of the branch: $$X_l - X_k \sim \mathcal{C}(0, \mbox{disp} \times t_l).$$
The method
argument specifies the type of likelihood that is computed:
method="reml"
:the dispersion parameter is fitted using the REML criterion,
obtained by re-rooting the tree to one of the tips.
The default tip used to reroot the tree is:
rootTip = which.min(colSums(cophenetic.phylo(tree)))
.
Any tip can be used, but this default empirically proved to be the most robust numerically;
method="random.root"
:the root value is assumed to be a random Cauchy variable, centered at root.value=0
,
and with a dispersion disp_root = disp * root.edge
;
method="fixed.root"
:the model is fitted conditionally on the root value root.value
,
i.e. with a model where the root value is fixed and inferred from the data.
phy <- ape::rphylo(5, 0.1, 0)
dat <- rTraitCauchy(n = 1, phy = phy, model = "cauchy", parameters = list(root.value = 0, disp = 1))
logDensityTipsCauchy(phy, dat, 0, 1, method = "fixed.root")
#> [1] -17.72164