In a recent
post,
in which I gave an R function to run Pagel's (1994) test for correlation between
two binary traits, I also posted (without much comment) some simple code to
*simulate* the character histories of binary traits *undergoing*
correlated evolution by Pagel's model.

Here is that code again with some explanation.

(1) First, just as a preliminary, let's load our packages & simulate a stochastic phylogeny:

```
library(phytools)
```

```
## Loading required package: ape
## Loading required package: maps
```

```
tree<-pbtree(n=300,scale=1)
```

(2) Next, we need to set up a matrix to simultaneously simulate our two
traits. Remember,
under
Pagel's model, character correlation are just transition rates in one
character that depend on the state of a second. To simulate under this model
we actually have to simultaneously simulate states for the two characters
at once. Here, as I did
last
time, I'll keep it simple. Let's say that the two characters have states
`a`

and `b`

. I'll let the transition rates into
states `a,a`

or `b,b`

be high (say `2.0`

),
but I will make the transition rates to any state `a,b`

or
`b,a`

*low* - say 0.4. This is what that looks like:

```
Q<-matrix(c(0,0.4,0.4,0,2,0,0,2,2,0,0,2,0,0.4,0.4,0),4,4,byrow=TRUE)
rownames(Q)<-colnames(Q)<-c("aa","ab","ba","bb")
diag(Q)<--rowSums(Q)
Q
```

```
## aa ab ba bb
## aa -0.8 0.4 0.4 0.0
## ab 2.0 -4.0 0.0 2.0
## ba 2.0 0.0 -4.0 2.0
## bb 0.0 0.4 0.4 -0.8
```

(3) Now, we can simulate both characters simultaneously up the tree using
the phytools function `sim.history`

.

```
tt<-sim.history(tree,Q)
```

```
## Note - the rate of substitution from i->j should be given by Q[j,i].
## Detecting that rows, not columns, of Q sum to zero :
## Transposing Q for internal calculations.
## Done simulation(s).
```

(4) This is the history of our two character state data; however we still
have to backtranslate this into the character histories for each of our
two traits. To do this, we can use `mergeMappedStates`

to merge
`a,a`

and `a,b`

(for instance) into `a`

for
character 1 and so on. Here is what that looks like:

```
t1<-mergeMappedStates(tt,c("aa","ab"),"a")
t1<-mergeMappedStates(t1,c("ba","bb"),"b")
t2<-mergeMappedStates(tt,c("aa","ba"),"a")
t2<-mergeMappedStates(t2,c("ab","bb"),"b")
t1$states<-getStates(t1,"tips")
t2$states<-getStates(t2,"tips")
```

(5) For fun, let's plot the two histories (just as we did last time) so that we can see that they are indeed highly correlated:

```
par(mfrow=c(1,2))
plotSimmap(t1,setNames(c("red","blue"),letters[1:2]),lwd=1,ftype="off")
plotSimmap(t2,setNames(c("red","blue"),letters[1:2]),lwd=1,ftype="off",
direction="leftwards")
```

(6) Finally, we can fit the Pagel (1994) model:

```
x<-getStates(t1,"tips")
y<-getStates(t2,"tips")
source("fitPagel.R")
.check.pkg<-phytools:::.check.pkg
fit<-fitPagel(tree,x,y)
fit
```

```
##
## Pagel's binary character correlation test:
##
## Indepedent model rate matrix:
## a|a a|b b|a b|b
## a|a -0.9753616 0.5606835 0.4146781 0.0000000
## a|b 2.7557800 -3.1704581 0.0000000 0.4146781
## b|a 2.7479439 0.0000000 -3.3086275 0.5606835
## b|b 0.0000000 2.7479439 2.7557800 -5.5037240
##
## Dependent model rate matrix:
## a|a a|b b|a b|b
## a|a -0.7760794 0.4795664 0.2965130 0.000000
## a|b 5.1666938 -7.2422057 0.0000000 2.075512
## b|a 2.0876921 0.0000000 -4.0078768 1.920185
## b|b 0.0000000 3.8020580 0.7652618 -4.567320
##
## Model fit:
## log-likelihood
## independent -186.3963
## dependent -173.9505
##
## Hypothesis test result:
## likelihood-ratio: 24.89146
## p-value: 5.290212e-05
```

That's it. Hopefully this is helpful to someone!

The Bestreplica watches. Here you can find almost swiss brand replica watches.Replica watches,one of the most famous brands,replica rolex watches ,Specialities watch for sale,Fast delivery and free shipping!

ReplyDelete