Among the list of *phytools* functions that may sometimes be overlooked,
I'd like to mention the simple but neat `density`

S3 method for
the *phytools* object class `"multiSimmap"`

.

The function computes (and an associated `plot`

method graphs)
a posterior density on the number of changes between states for the
character across our sample of stochastic character maps. The following
is a short demo of how that looks.

First, here are some data & a tree:

```
library(phytools)
tree
```

```
##
## Phylogenetic tree with 150 tips and 149 internal nodes.
##
## Tip labels:
## t23, t50, t51, t72, t73, t102, ...
##
## Rooted; includes branch lengths.
```

```
x
```

```
## t23 t50 t51 t72 t73 t102 t103 t93 t64 t94 t95 t11 t19 t89 t90
## a a a a a a a a a a a b a a a
## t3 t54 t55 t56 t81 t82 t27 t57 t91 t92 t38 t39 t37 t47 t145
## a a a a b b b b a b b b b b b
## t146 t2 t24 t58 t59 t60 t61 t123 t124 t6 t42 t43 t100 t101 t16
## b b a a a a a b b a a b a b a
## t8 t9 t7 t74 t75 t69 t143 t144 t41 t132 t133 t1 t129 t130 t35
## a b a a a a b b a b b b b b b
## t36 t31 t65 t119 t120 t118 t87 t52 t53 t108 t109 t141 t142 t20 t21
## b b b a b b b b b a a b b a b
## t25 t26 t45 t46 t104 t105 t98 t99 t30 t83 t84 t138 t139 t106 t107
## a a a a a a a a a a a a a a a
## t85 t86 t71 t12 t67 t68 t62 t63 t28 t29 t33 t34 t136 t137 t125
## a a b a a a a a b b a a b b a
## t126 t10 t4 t5 t17 t18 t114 t115 t22 t127 t128 t96 t97 t70 t147
## a a b a a a a a a a a a a a a
## t148 t131 t48 t32 t79 t80 t66 t116 t117 t88 t112 t113 t49 t121 t122
## a a a a a a b b b b b b b b b
## t44 t76 t77 t110 t111 t13 t14 t140 t149 t150 t78 t134 t135 t40 t15
## b b b b b b b a a a a a a a b
## Levels: a b
```

Now, let's undertake our stochastic mapping:

```
mtrees<-make.simmap(tree,x,nsim=500,model="ARD")
```

```
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## a b
## a -0.8282081 0.8282081
## b 0.9020183 -0.9020183
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## a b
## 0.5 0.5
```

```
## Done.
```

```
mtrees
```

```
## 500 phylogenetic trees with mapped discrete characters
```

Finally, we'll use our `density`

method:

```
d<-density(mtrees)
d
```

```
##
## Distribution of changes from stochastic mapping:
## a->b b->a
## Min. :10 Min. :3
## Median :20 Median :10
## Mean :19.38 Mean :10.64
## Max. :28 Max. :23
##
## 95% HPD interval(a->b): [13, 24]
## 95% HPD interval(b->a): [4, 17]
```

```
par(family="mono")
plot(d)
```

Just for interest, the tree & data were simulated as follows:

```
library(phytools)
Q<-matrix(c(-0.8,0.8,0.4,-0.4),2,2,byrow=TRUE)
rownames(Q)<-colnames(Q)<-letters[1:2]
tree<-pbtree(n=150,scale=1)
x<-sim.Mk(tree,Q)
```

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.