I just
pushed
a pretty substantive update to the phytools function `ratebytree`

for comparing rates (and
processes) of evolutionary change among trees. The method essentially extends the 'censored' model of
O'Meara et al. (2006), but applied to the problem of
comparing between different trees and modeling processes other than Brownian motion evolution of
continuous traits.

What I had not included, to this point, was any posthoc comparison of rates (or evolutionary parameter estimates - for processes other than BM) between trees. Thus - for more than two trees, we might be able to tell that evolutionary rates are different, but perhaps not which ones differ significantly one from the other.

I have now fixed this in two ways. First, I obtain the Hessian (curvature) matrix from the likelihood
surface around the MLE from which I could calculate variances and standard errors of the parameter
estimates for each model. These are now output as part of the `"ratebytree"`

object and
printed by the S3 `print`

method of the object class. Second, I developed a new method called
`posthoc`

which conducts a set of pairwise posthoc tests between each pair of trees and
reports the t-statistics and P-values for each comparison. The user can optionally supply a value to the
argument `p.adjust.method`

to adjust the reported P-values due to multiple testing.
(Alternatively, one could just use a stricter rejection threshold.)

Here I will demonstrate using some simulated trees & data.

```
library(phytools)
packageVersion("phytools")
```

```
## [1] '0.6.26'
```

```
print(trees,details=TRUE)
```

```
## 3 phylogenetic trees
## tree 1 : 40 tips
## tree 2 : 60 tips
## tree 3 : 40 tips
```

```
par(mfrow=c(1,3))
## x
nulo<-mapply(phenogram,tree=trees,x=x,
MoreArgs=list(spread.cost=c(1,0),ftype="i",
ylim=range(x),spread.range=range(x),
col=make.transparent("blue",0.5)))
```

```
## Optimizing the positions of the tip labels...
```

```
## Optimizing the positions of the tip labels...
```

```
## Optimizing the positions of the tip labels...
```

```
## y
nulo<-mapply(phenogram,tree=trees,x=y,
MoreArgs=list(spread.cost=c(1,0),ftype="i",
ylim=range(y),spread.range=range(y),
col=make.transparent("blue",0.5)))
```

```
## Optimizing the positions of the tip labels...
```

```
## Optimizing the positions of the tip labels...
```

```
## Optimizing the positions of the tip labels...
```

Now, at a glance, it may be kind of hard to tell which character (if either) has evolved with a difference in rate between trees - and which has not. Luckily, we have a method to test this! Let's see what it shows.

```
fitBM.x<-ratebytree(trees,x,model="BM")
fitBM.x
```

```
## ML common-rate model:
## s^2 a[1] a[2] a[3] k logL
## value 1.4619 0.9017 -0.6506 -0.2191 4 -219.2487
## SE 0.1747 1.0889 0.774 0.8131
##
## ML multi-rate model:
## s^2[1] s^2[2] s^2[3] a[1] a[2] a[3] k logL
## value 0.9336 2.1605 0.9423 0.9017 -0.6506 -0.2191 6 -213.2143
## SE 0.2088 0.3944 0.2107 0.8701 0.9409 0.6528
##
## Likelihood ratio: 12.0689
## P-value (based on X^2): 0.0024
##
## R thinks it has found the ML solution.
```

```
posthoc(fitBM.x,p.adjust.method="holm")
```

```
##
## Post-hoc test for "BM" model.
## (Comparison is of estimated values of sigma^2.)
##
## t df P
## tree1 vs. 2 -2.7492 92.1466 0.0216
## tree1 vs. 3 -0.0295 75.9934 0.9766
## tree2 vs. 3 2.7241 92.4300 0.0216
##
## P-values adjusted using method="holm".
```

```
fitBM.y<-ratebytree(trees,y,model="BM")
fitBM.y
```

```
## ML common-rate model:
## s^2 a[1] a[2] a[3] k logL
## value 1.0758 -1.6204 0.1808 0.5889 4 -197.7798
## SE 0.1286 0.934 0.6639 0.6975
##
## ML multi-rate model:
## s^2[1] s^2[2] s^2[3] a[1] a[2] a[3] k logL
## value 0.6775 1.2436 1.2222 -1.6204 0.1808 0.5889 6 -195.4351
## SE 0.1515 0.2271 0.2733 0.7412 0.7138 0.7435
##
## Likelihood ratio: 4.6894
## P-value (based on X^2): 0.0959
##
## R thinks it has found the ML solution.
```

```
posthoc(fitBM.y,p.adjust.method="holm")
```

```
##
## Post-hoc test for "BM" model.
## (Comparison is of estimated values of sigma^2.)
##
## t df P
## tree1 vs. 2 -2.0741 95.9917 0.1222
## tree1 vs. 3 -1.7433 59.3364 0.1729
## tree2 vs. 3 0.0602 71.1460 0.9522
##
## P-values adjusted using method="holm".
```

This result seems to show that trait *x* exhibits a difference in evolutionary rate between trees - while
trait *y* shows no evidence for a rate difference. In fact, this is precisely what we simulated:

```
trees<-c(pbtree(n=40),pbtree(n=60),pbtree(n=40))
class(trees)<-"multiPhylo"
x<-mapply(fastBM,tree=trees,sig2=c(1,2,1),SIMPLIFY=FALSE)
y<-sapply(trees,fastBM)
```

Neat.

## No comments:

## Post a Comment