I just
rewrote
the function that is used internally to compute the likelihood of any
transition matrix `Q`

given the data & tree; and to estimate
the ML value of `Q`

. This function is highly similar (but slightly
different in important ways), both functionally & structurally, to
`ace(...,type="discrete")`

in the ape package. It is also in some
ways functionally redundant with `ace`

and `fitDiscrete`

in that they fit the same model - although they differ in that the prior
distribution on the root node can be controlled in this function. In fact, this
is the main reason for this 're-write' because in previous versions of
phytools, I suspect that this prior was being used improperly in the
calculations. Now, when this function is used internally by
`make.simmap`

(which is it's main purpose, although it can also be
used alone) the prior, `pi`

, can influence not only the states
sampled at the root (and thus, indirectly, for other nodes in the tree),
but it can also affect the estimated value of `Q`

under likelihood.

Here's a quick demo of using the function alone to fit the M*k*
discrete character evolution model:

```
library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)
```

Simulate tree & data:

```
library(phytools)
tree<-pbtree(n=100,scale=1)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
Q
```

```
## a b c
## a -1 1 0
## b 1 -2 1
## c 0 1 -1
```

```
x<-sim.history(tree,Q,anc="a")$states
```

```
## Done simulation(s).
```

```
x
```

```
## t16 t17 t63 t64 t39 t29 t54 t55 t5 t18 t50 t51 t33 t19 t20
## "a" "a" "a" "a" "a" "a" "c" "b" "c" "b" "c" "c" "c" "b" "c"
## t10 t99 t100 t88 t76 t79 t80 t93 t94 t89 t90 t44 t95 t96 t68
## "c" "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "c"
## t69 t65 t81 t82 t77 t78 t60 t9 t40 t83 t84 t72 t73 t56 t57
## "c" "c" "c" "c" "c" "c" "c" "a" "c" "b" "b" "b" "b" "c" "c"
## t26 t21 t61 t62 t42 t43 t52 t53 t6 t7 t2 t58 t59 t11 t12
## "b" "b" "b" "b" "b" "b" "c" "c" "a" "a" "a" "a" "a" "a" "a"
## t1 t22 t23 t3 t4 t45 t46 t27 t85 t86 t35 t36 t37 t38 t28
## "a" "c" "c" "b" "a" "a" "a" "a" "a" "a" "b" "b" "b" "a" "a"
## t24 t25 t30 t74 t75 t14 t15 t13 t8 t87 t91 t92 t70 t71 t49
## "a" "b" "b" "a" "b" "a" "a" "a" "a" "b" "b" "b" "b" "b" "b"
## t47 t48 t41 t34 t66 t67 t31 t32 t97 t98
## "b" "a" "b" "b" "b" "b" "a" "b" "a" "a"
```

Fit the model:

```
obj<-fitMk(tree,as.factor(x),model="SYM")
obj
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b c
## a -0.930216 0.930216 0.00000
## b 0.930216 -2.094586 1.16437
## c 0.000000 1.164370 -1.16437
##
## Fitted (or set) value of pi:
## a b c
## 0.3333333 0.3333333 0.3333333
##
## Log-likelihood: -66.577856
```

It is also used internally by `make.simmap`

. Here, with a
prior that is not flat:

```
trees<-make.simmap(tree,x,pi=setNames(c(1,0,0),letters[1:3]),nsim=100)
```

```
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
```

```
## a b c
## a -0.8965016 0.8965016 0.00000
## b 0.8965016 -2.0898515 1.19335
## c 0.0000000 1.1933498 -1.19335
```

```
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
```

```
## a b c
## 1 0 0
```

```
## Done.
```

```
obj<-summary(trees)
obj
```

```
## 100 trees with a mapped discrete character with states:
## a, b, c
##
## trees have 34.17 changes between states on average
##
## changes are of the following types:
## a,b a,c b,a b,c c,a c,b
## x->y 12 0 6.36 9.63 0 6.18
##
## mean total time spent in each state is:
## a b c total
## raw 12.1857664 8.5412422 4.3792713 25.10628
## prop 0.4853673 0.3402034 0.1744293 1.00000
```

```
plot(obj,ftype="off",type="fan")
```

That's it.

## No comments:

## Post a Comment

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