As I mentioned
the
other day, during my sabbatical at the *Universidad de los Andes* here in
Bogotá I submitted a paper with some colleagues describing the function
`ratebytree`

for comparing the evolutionary rates of a continuous character between
trees.

The manuscript was favorably reviewed, but in their typically pesky* (*said with affection)
manner, reviewers asked us to add different models of evolution - the feature that I added
last
week. They *also* asked us to extend the approach to discrete character evolution.
This was something we pointed out could be done in the manuscript, but had not explicitly
addressed. It is theoretically straightforward - although not
completely
trivial - particularly if you want your new function updates to be accompanied by
legible `print`

methods & so on.

I have, however, now added this feature to the function `ratebytree`

. The function
can take an argument, `type`

(which should assume the value `"continuous"`

or `"discrete"`

); however, if not provided, the function will attempt to ascertain
the character type.

Here is a quick demo using some simple simulated data:

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

```
## [1] '0.6.25'
```

First, our trees & data:

```
t1
```

```
##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
##
## Rooted; includes branch lengths.
```

```
x1
```

```
## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
## c c c c a c b c b c a c b a b c b b a c a c c c c c
## Levels: a b c
```

```
t2
```

```
##
## Phylogenetic tree with 40 tips and 39 internal nodes.
##
## Tip labels:
## t7, t34, t35, t3, t29, t30, ...
##
## Rooted; includes branch lengths.
```

```
x2
```

```
## t7 t34 t35 t3 t29 t30 t31 t14 t15 t2 t39 t40 t9 t8 t21 t22 t32 t33
## b a a c b b b b c b a a c c a a c c
## t23 t24 t19 t20 t11 t5 t6 t25 t26 t4 t36 t37 t38 t12 t13 t27 t28 t16
## c c b b c c a c c a b b b c c c b c
## t1 t10 t17 t18
## a a a a
## Levels: a b c
```

```
## or
dotTree(t3,x3,ftype="off")
```

```
x<-list(x1,x2,x3)
x
```

```
## [[1]]
## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
## c c c c a c b c b c a c b a b c b b a c a c c c c c
## Levels: a b c
##
## [[2]]
## t7 t34 t35 t3 t29 t30 t31 t14 t15 t2 t39 t40 t9 t8 t21 t22 t32 t33
## b a a c b b b b c b a a c c a a c c
## t23 t24 t19 t20 t11 t5 t6 t25 t26 t4 t36 t37 t38 t12 t13 t27 t28 t16
## c c b b c c a c c a b b b c c c b c
## t1 t10 t17 t18
## a a a a
## Levels: a b c
##
## [[3]]
## t19 t20 t3 t11 t12 t36 t37 t21 t24 t25 t10 t13 t14 t2 t7 t8 t9 t16
## a a a b a c b a a a b a a a c b a c
## t28 t29 t32 t33 t17 t34 t35 t18 t1 t39 t40 t22 t23 t41 t42 t4 t5 t26
## b b b b b b b b b c c c c b b b a b
## t27 t6 t49 t50 t15 t43 t44 t45 t46 t30 t31 t47 t48 t38
## b b c c b c c b b c c b b b
## Levels: a b c
```

```
trees<-c(t1,t2,t3)
trees
```

```
## 3 phylogenetic trees
```

OK, now let's try it! First, I will assume a symmetric model of transitions between states:

```
fitSYM<-ratebytree(trees,x,type="discrete",model="SYM")
fitSYM
```

```
## ML common-rate model:
## a<->b a<->c b<->c k logL
## value 1.3136 0.0815 3.0745 3 -100.7769
##
## Model fitting method was "optim".
##
## ML multi-rate model:
## a<->b a<->c b<->c k logL
## tree1 3.5598 0 5.2209 9 -99.2075
## tree2 0.2588 0.5751 1.5436
## tree3 1.3966 0 3.2116
##
## Model fitting method was "nlminb".
##
## Likelihood ratio: 3.1388
## P-value (based on X^2): 0.7912
```

This tells us that we have little evidence to reject a constant (set of) rate(s) across our trees. We can nonetheless try simpler or more complex models and see if the outcome is similar or different. Remember, I said the function would attempt to ascertain if the character is discretely or continuously-valued, so let's try that:

```
fitER<-ratebytree(trees,x)
fitER
```

```
## ML common-rate model:
## q k logL
## value 1.1995 1 -104.2243
##
## Model fitting method was "optimize".
##
## ML multi-rate model:
## q k logL
## tree1 2.1052 3 -103.4967
## tree2 0.8053
## tree3 1.4482
##
## Model fitting method was "nlminb".
##
## Likelihood ratio: 1.4552
## P-value (based on X^2): 0.4831
```

```
fitARD<-ratebytree(trees,x,model="ARD")
fitARD
```

```
## ML common-rate model:
## b->a c->a a->b c->b a->c b->c k logL
## value -0.0483 1.1847 2.3293 1.5367 -0.7321 3.8122 6 -99.4653
##
## Model fitting method was "optim".
##
## ML multi-rate model:
## b->a c->a a->b c->b a->c b->c k logL
## tree1 4.0703 0 4.1987 7.5158 0 16.8153 18 -94.8226
## tree2 0.1143 0.5863 0 1.2591 0.1455 2.3226
## tree3 0.6492 0.2584 1.9163 4.7043 0 2.7484
##
## Model fitting method was "nlminb".
##
## Likelihood ratio: 9.2853
## P-value (based on X^2): 0.6784
```

By chance, these data were *simulated* as follows - and thus without a difference in rate
or process between trees:

```
t1<-pbtree(n=26,tip.label=LETTERS,scale=1)
t2<-pbtree(n=40,scale=2)
t3<-pbtree(n=50,scale=1)
Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x1<-as.factor(sim.history(t1,Q)$states)
x2<-as.factor(sim.history(t2,Q)$states)
x3<-as.factor(sim.history(t3,Q)$states)
```

Why don't we simulate under conditions in which the process *differs* between trees and
see if we obtain a different result:

```
t4<-pbtree(n=60,scale=1)
t5<-pbtree(n=70,scale=2)
t6<-pbtree(n=80,scale=1)
Q1<-matrix(c(-1,1,1,-1),2,2)
Q2<-matrix(c(-2,2,2,-2),2,2)
x4<-as.factor(sim.history(t4,Q1)$states)
x5<-as.factor(sim.history(t5,Q1)$states)
x6<-as.factor(sim.history(t6,Q2)$states)
```

Now we're ready to fit our models:

```
fitER<-ratebytree(c(t4,t5,t6),list(x4,x5,x6))
fitER
```

```
## ML common-rate model:
## q k logL
## value 1.2509 1 -116.4188
##
## Model fitting method was "optimize".
##
## ML multi-rate model:
## q k logL
## tree1 0.8768 3 -112.826
## tree2 0.7286
## tree3 2.7957
##
## Model fitting method was "nlminb".
##
## Likelihood ratio: 7.1856
## P-value (based on X^2): 0.0275
```

This suggests that a model in which the different trees are permitted to have different rates of evolution fits better than one in which they are all constrained to have equal rates - just as we simulated.

Finally, to get this functionality you should update phytools from GitHub as follows:

```
library(devtools)
install_github("liamrevell/phytools")
```

More later.

## No comments:

## Post a Comment