fitCauchy {cauphy}  R Documentation 
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 
trait 
named vector of traits at the tips. 
model 
a model for the trait evolution. One of 
method 
the method used to fit the process.
One of 
starting.value 
starting value for the parameters of the Cauchy.
This should be a named list, with 
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 
hessian 
if 
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 
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 interquartile 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 rerooting the tree to one of the tips. See
logDensityTipsCauchy
for the default choice of the rerooting tip; method="random.root"
:
the root value is assumed to be a random Cauchy variable, centered at
x0=0
, and with a dispersiondisp_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 rescaled
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 
disp 
the ml or reml estimate of the dispersion parameter 
lambda 
the ml or reml estimate of the lambda parameter (for 
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 
phy 
the phylogenetic tree 
method 
the method used (one of 
random.root 

reml 

root_tip_reml 
name of the tip used to reroot the tree (for 
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)