fitCauchy: Model fitting for a Cauchy Process

In cauphy: Trait Evolution on Phylogenies Using the Cauchy Process

Model fitting for a Cauchy Process

Description

Fit the Cauchy process on a phylogeny, using numerical optimization.

Usage

fitCauchy(
  phy,
  trait,
  model = c("cauchy", "lambda"),
  method = c("reml", "random.root", "fixed.root"),
  starting.value = list(x0 = NULL, disp = NULL, lambda = NULL),
  lower.bound = list(disp = 0, lambda = 0),
  upper.bound = list(disp = Inf, lambda = NULL),
  root.edge = 100,
  hessian = FALSE,
  optim = c("local", "global"),
  method.init.disp = c("Qn", "Sn", "MAD", "IQR")
)

Arguments

phy	a phylogenetic tree of class phylo.
trait	named vector of traits at the tips.
model	a model for the trait evolution. One of "cauchy" or "lambda" (see Details).
method	the method used to fit the process. One of reml (the default), fixed.root or random.root. See Details.
starting.value	starting value for the parameters of the Cauchy. This should be a named list, with x0 and disp the root starting value and the dispersion parameter. The default initial values are computed from standard statistics used on (independent) Cauchy variables, see Details.
lower.bound	named list with lower bound values for the parameters. See Details for the default values.
upper.bound	named list with upper bound values for the parameters. See Details for the default values.
root.edge	multiplicative factor for the root dispersion, equal to the length of the root edge. Ignored if method!=random.root.
hessian	if TRUE, then the numerical hessian is computed, for confidence interval computations. See compute_vcov.
optim	if "local", only a local optimization around the initial parameter values is performed (the default). If "global", a global maximization is attempted using the "MLSL" approach (see nloptr).
method.init.disp	the initialization method for the dispersion. One of "Qn", "Sn", "MAD", "IQR". Default to the "Qn" statistics. See Details.

Details

For the default model="cauchy", the parameters of the Cauchy Process (CP) are disp, the dispersion of the process, and x0, 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).

Unless specified by the user, the initial values for the parameters are taken according to the following heuristics:

x0:

is the trimmed mean of the trait, keeping only 24% of the observations, as advocated in Rothenberg et al. 1964 (for method="fixed.root");

disp:

is initialized from the trait centered and normalized by tip heights, with one of the following statistics, taken from Rousseeuw & Croux 1993:

IQR:

half of the inter-quartile range (see IQR);

MAD:

median absolute deviation with constant equal to 1 (see mad);

Sn:

Sn statistics with constant 0.7071 (see Sn);

Qn:

(default) Qn statistics with constant 1.2071 (see Qn).

Unless specified by the user, x0 is taken to be unbounded, disp positive unbounded.

The method argument specifies the method used for the fit:

method="reml":

the dispersion parameter is fitted using the REML criterion, obtained by re-rooting the tree to one of the tips. See logDensityTipsCauchy for the default choice of the re-rooting tip;

method="random.root":

the root value is assumed to be a random Cauchy variable, centered at x0=0, and with a dispersion disp_root = disp * root.edge;

method="fixed.root":

the model is fitted conditionally on the root value x0, i.e. with a model where the root value is fixed and inferred from the data.

In the first two cases, the optimization is done on the dispersion only, while in the last case the optimization is on the root value and the dispersion.

The function uses nloptr for the numerical optimization of the (restricted) likelihood, computed with function logDensityTipsCauchy. It uses algorithms BOBYQA and MLSL_LDS for local and global optimization.

If model="lambda", the CP is fit on a tree with branch lengths re-scaled using the Pagel's lambda transform (see transf.branch.lengths), and the lambda value is estimated using numerical optimization. The default initial value for the lambda parameter is computed using adequate robust moments. The default maximum value is computed using phytools:::maxLambda, and is the ratio between the maximum height of a tip node over the maximum height of an internal node. This can be larger than 1. The default minimum value is 0.

Value

An object of S3 class cauphyfit, with fields:

x0	the fitted starting value (for method="fixed.root")
disp	the ml or reml estimate of the dispersion parameter
lambda	the ml or reml estimate of the lambda parameter (for model="lambda")
logLik	the maximum of the log (restricted) likelihood
p	the number of parameters of the model
aic	the AIC value of the model
trait	the named vector of traits at the tips used in the fit
y	the named vector of traits at the tips used in the fit
n	the number of tips in the tree
d	the number of dependent variables
call	the original call of the function
model	the phylogenetic model (one of "cauchy" or "lambda")
phy	the phylogenetic tree
method	the method used (one of "reml", "fixed.root", "random.root")
random.root	TRUE if method="random.root"
reml	TRUE if method="reml"
root_tip_reml	name of the tip used to reroot the tree (for method="reml")

References

Bastide, P. and Didier, G. 2023. The Cauchy Process on Phylogenies: a Tractable Model for Pulsed Evolution. Systematic Biology. doi:10.1093/sysbio/syad053.

Rothenberg T. J., Fisher F. M., Tilanus C. B. 1964. A Note on Estimation from a Cauchy Sample. Journal of the American Statistical Association. 59:460–463.

Rousseeuw P.J., Croux C. 1993. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association. 88:1273–1283.

See Also

confint.cauphyfit, profile.cauphyfit, ancestral, increment, logDensityTipsCauchy, cauphylm, fitContinuous

Examples

# Simulate tree and data
set.seed(1289)
phy <- ape::rphylo(20, 0.1, 0)
dat <- rTraitCauchy(n = 1, phy = phy, model = "cauchy",
                    parameters = list(root.value = 10, disp = 0.1))
# Fit the data
fit <- fitCauchy(phy, dat, model = "cauchy", method = "reml")
fit
# Approximate confidence intervals
confint(fit)
# Profile likelihood
pl <- profile(fit)
plot(pl)

cauphy documentation built on Oct. 1, 2024, 5:08 p.m.