shape.statistic: Computes the log of the likelihood ratio (yule/pda)

In apTreeshape: Analyses of Phylogenetic Treeshape

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 log of 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 log 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 log 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 <[e-mail]>
Nicolas Bortolussi <[e-mail]>
Eric Durand <[e-mail]>
Olivier Fran?ois <[e-mail]>

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

data(universal.treeshape)
tree <- universal.treeshape
plot(tree)
summary(tree)

likelihood.test(tree,  model = "yule", alternative = "two.sided")
likelihood.test(tree,  model = "pda", alternative = "two.sided")

## 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(100,tip.number=50,model="yule"),FUN=shape.statistic,
      norm="yule"), freq=FALSE, main=main, xlab=xlab)
#Increase the numner of trees >100 for a better fit
## Does it fit the Gaussian distribution with mean=0 and sd=1 ?
x<-seq(-3,3,by=0.001)
lines(x,dnorm(x))

Example output

Loading required package: ape

Attaching package: 'apTreeshape'

The following object is masked from 'package:ape':

    is.binary.phylo

The "ward" method has been renamed to "ward.D"; note new "ward.D2"

Phylogenetic tree shape: object  tree  of class 'treeshape'

Number of tips: 26 

Colless' shape statistic: 84 
Expected value (Yule model):  55.69629  Standard Deviation:  21.15318 
Expected value (PDA model):  234.9822  Standard Deviation:  58.052 
Test of the Yule hypothesis: 
statistic = 4.330559 
p.value = 1.48731e-05 
alternative hypothesis: the tree does not fit the Yule model

Note: the p.value was computed according to a normal approximation
Test of the PDA hypothesis: 
statistic = -0.1646905 
p.value = 0.8691876 
alternative hypothesis: the tree does not fit the PDA model

Note: the p.value was computed according to a normal approximation

apTreeshape documentation built on Jan. 8, 2021, 2:07 a.m.