As implicitly
promised
in a post yesterday, I have now added a `plot`

method for fitted M*k*
models from the geiger function `fitDiscrete`

. I felt like this was a good
idea because I believe the geiger implementation to be more robust (that is, it uses
a more rigorous optimization process) than does `fitMk`

in my package.
`fitMk`

, by the way, recycles code extensively from `ace`

in the
ape package.

What I elected to do was to first pull out all the relevant attributes of the fitted model
from the geiger object (an object of class `"gfit"`

resulting from a call to
`"fitDiscrete"`

), then use these values to create a new object of class
`"fitMk"`

internally, then just call phytools S3 method on that object.

Here is what that wrapper code looks like. Note that an object of class `"gfit"`

in geiger can be from a call to `fitDiscrete`

*or* `fitContinuous`

,
so we have to check for that first.

```
## S3 plot method for objects resulting from fitDiscrete
plot.gfit<-function(x,...){
if("mkn"%in%class(x$lik)==FALSE){
stop("Sorry. No plot method presently available for objects of this type.")
} else {
obj<-list()
QQ<-.Qmatrix.from.gfit(x)
obj$states<-colnames(QQ)
m<-length(obj$states)
obj$index.matrix<-matrix(NA,m,m)
k<-m*(m-1)
obj$index.matrix[col(obj$index.matrix)!=row(obj$index.matrix)]<-1:k
obj$rates<-QQ[sapply(1:k,function(x,y) which(x==y),obj$index.matrix)]
class(obj)<-"fitMk"
plot(obj,...)
}
}
```

Now let's take our data & tree and fit a model:

```
library(phytools)
library(geiger)
tree
```

```
##
## Phylogenetic tree with 100 tips and 99 internal nodes.
##
## Tip labels:
## t67, t68, t76, t78, t79, t74, ...
##
## Rooted; includes branch lengths.
```

```
x
```

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

```
fit.phytools<-fitMk(tree,x,model="ARD")
plot(fit.phytools,main="fitMk model=\"ARD\"",show.zeros=FALSE,
cex.traits=0.9)
```

```
fit.geiger<-fitDiscrete(tree,x,model="ARD")
plot(fit.geiger,main="fitDiscrete model=\"ARD\"",show.zeros=FALSE,
cex.traits=0.9)
```

Both fitted models are basically identical. Although the configuration of the states in the plotting space are different, it is easy to see that the estimated rates are the same. In other cases, `fitDiscrete`

and `fitMk`

may differ, either because one or the other failed to converge on the MLE of `Q`

, or because the two implementations make slightly different assumptions about the root.

FYI, the generating model in this case was as follows:

```
print(t(Q))
```

```
## green blue purple red
## green -2 1 1 0
## blue 1 -3 1 1
## purple 1 1 -3 1
## red 0 1 1 -2
```

so we are pretty close, particularly in that the two transition classes with their rates set to zero are the only ones with estimated rates of zero in the final fitted model.

Note that the function above uses the geiger internal function (borrowed by phytools)
`.Qmatrix.from.gfit`

, so to work that function will need to be loaded into the
name space:

```
.Qmatrix.from.gfit<-phytools:::.Qmatrix.from.gfit
```

Finally, here's how the data were simulated:

```
tree<-pbtree(n=100)
Q<-matrix(c(-2,1,1,0,
1,-3,1,1,
1,1,-3,1,
0,1,1,-2),4,4)
rownames(Q)<-colnames(Q)<-c("green","blue","purple","red")
x<-factor(sim.history(tree,Q)$states,
levels=c("green","blue","purple","red"))
```

## No comments:

## Post a Comment