shape.statistic {apTreeshape}R Documentation

Computes the log of the likelihood ratio (yule/pda) of the given tree

Description

shape.statistic computes the logarithm of the ratio of the likelihoods under the Yule model and the PDA model of the given tree.

Usage

shape.statistic(tree, norm="null")

Arguments

tree An object of class "treeshape".
norm A character string equals to "null" for no normalization, "yule" for the Yule model normalization or "pda" for the pda normalization.

Details

The likelihood ratio is proportional to

sum(log(N(v)-1)),

for all internal node v (where N(v) is the number of internal nodes descending from the node v ). The ratio of the likelihoods enables to build the most powerful test of the Yule model against the PDA one. (Neyman-Pearson lemma).

Under the PDA model, the ratio has approximate Gaussian distribution of mean ~ 2.03*n-3.545*sqrt(n-1) and variance ~ 2.45*(n-1)*log(n-1), where n is the number of tips of the tree. The Gaussian approximation is accurate for very large n (n greater than 10000(?)). The normalization of the ratio uses tabulated empirical values of variances estimated from Monte Carlo simulations. The normalization uses the formula:

variance ~ 1.570*n*log(n)-5.674*n+3.602*sqrt(n)+14.915

Under the Yule model, the ratio has approximate Gaussian distribution of mean = 1.204*n-log(n-1)-2 and variance = 0.168*n-0.710, where n is the number of tips of the tree. The Gaussian approximation is accurate for n greater than 20.

Value

An object of class numeric containing the shape statistic of the tree.

Author(s)

Michael Blum <michael.blum@imag.fr>
Nicolas Bortolussi <nicolas.bortolussi@imag.fr>
Eric Durand <eric.durand@imag.fr>
Olivier François <olivier.francois@imag.fr>

References

Fill, J. A. (1996), On the Distribution of Binary Search Trees under the Random Permutation Model. Random Structures and Algorithms, 8, 1 – 25, for more details about the normalization and proofs.

Examples


## Histogram of shape statistics for a list of Yule trees 
##      (may take some time to compute)
main="Histogram of shape statistics"; xlab="shape statistic"    
hist(sapply(rtreeshape(1000,tip.number=100,model="yule"),FUN=shape.statistic,
      norm="yule"), freq=FALSE, main=main, xlab=xlab)

## Does it fit the Gaussian distribution with mean=0 and sd=1 ?
x<-seq(-3,3,by=0.001)
lines(x,dnorm(x))

[Package apTreeshape version 1.0.0 Index]