New bug for ordered, directional character evolution model in fitMk

I just discovered what seems to be a new bug in the phytools function fitMk for fitting the Mk model of discrete trait evolution. fitMk uses code modified from ape::ace, but has been adapted to do a variety of different things - such as accept ambiguous tip states (coded as prior probabilities in a matrix), and a flexible prior density at the global root of the tree. I believe I've tracked the bug down to the internally used function matexpo for matrix exponentiation from the ape package. When I substitute phytools::expm (which actually just wraps msm::MatrixExp the function works again.

I've only gotten the bug to appear for ordered, directional models. That is, for instance, a model in which we permit A→B→C, but not the reverse.

Here's an example using a dataset of squamates & digit number from a study by Matt Brandley et al. (2008). (Actually, the data are simplified because in the original data some species are variable in digit number for a particular limb, which is a nuance I've ignored.)

It might make sense to fit a model in which digit loss occurred 5→4→3→2→1→0, without jumps between non-adjacent digits numbers & without digit readquisition.* (I'm not proposing to imagine that neither is possible, let's just suppose we want to fit this seemingly reasonable model.)

We can do that as follows:

library(phytools)
library(geiger)
sqTree

## 
## Phylogenetic tree with 258 tips and 257 internal nodes.
## 
## Tip labels:
##  Abronia_graminea, Mesaspis_moreleti, Gerrhonotus_liocephalus, Barisia_imbricata, Elgaria_coerulea, Elgaria_kingii, ...
## 
## Rooted; includes branch lengths.

head(toes,n=30)

##     Agamodon_anguliceps        Amphisbaena_alba           Bipes_biporus 
##                       0                       0                       0 
##       Bipes_caniculatus       Bipes_tridactylus         Blanus_cinereus 
##                       0                       0                       0 
##   Chirindia_swimmertoni   Diplometopon_zarudnyi       Geocalamus_acutus 
##                       0                       0                       0 
##     Monopeltis_capensis      Rhineura_floridana   Trogonophis_wiegmanni 
##                       0                       0                       0 
##        Abronia_graminea         Anguis_fragilis   Anniella_geronimensis 
##                       5                       0                       0 
##        Anniella_pulchra       Barisia_imbricata   Celestus_enneagrammus 
##                       0                       5                       5 
##  Diploglossus_bilobatus      Diploglossus_pleei        Elgaria_coerulea 
##                       5                       5                       5 
##          Elgaria_kingii   Elgaria_multicarinata     Elgaria_panamintina 
##                       5                       5                       5 
##   Elgaria_paucicarinata Gerrhonotus_liocephalus       Mesaspis_moreleti 
##                       5                       5                       5 
##       Ophiodes_striatus       Ophisaurus_apodus   Ophisaurus_attenuatus 
##                       1                       1                       0 
## Levels: 0 1 2 3 4 5

I realize that this isn't the best way to visualize our data - but nevertheless:

plotTree.barplot(sqTree,
    setNames(as.numeric(toes)-1,names(toes)),
    args.plotTree=list(fsize=0.3),args.barplot=list(space=0,
    xlab="number of digits in the hindlimb",border="grey"))

Now we can set up & fit our model:

directional<-matrix(c(0,0,0,0,0,0,
    1,0,0,0,0,0,
    0,2,0,0,0,0,
    0,0,3,0,0,0,
    0,0,0,4,0,0,
    0,0,0,0,5,0),6,6,byrow=TRUE,
    dimnames=list(levels(toes),levels(toes)))
directional

##   0 1 2 3 4 5
## 0 0 0 0 0 0 0
## 1 1 0 0 0 0 0
## 2 0 2 0 0 0 0
## 3 0 0 3 0 0 0
## 4 0 0 0 4 0 0
## 5 0 0 0 0 5 0

fitDirectional<-fitDiscrete(sqTree,toes,model=directional)

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

## Warning in fitDiscrete(sqTree, toes, model = directional): Parameter estimates appear at bounds:
##  q12
##  q13
##  q14
##  q15
##  q16
##  q23
##  q24
##  q25
##  q26
##  q31
##  q34
##  q35
##  q36
##  q41
##  q42
##  q45
##  q46
##  q51
##  q52
##  q53
##  q56
##  q61
##  q62
##  q63
##  q64

fitDirectional

## GEIGER-fitted comparative model of discrete data
##  fitted Q matrix:
##                0           1          2          3            4
##     0 0.00000000  0.00000000  0.0000000  0.0000000  0.000000000
##     1 0.02314697 -0.02314697  0.0000000  0.0000000  0.000000000
##     2 0.00000000  0.17375162 -0.1737516  0.0000000  0.000000000
##     3 0.00000000  0.00000000  0.2123601 -0.2123601  0.000000000
##     4 0.00000000  0.00000000  0.0000000  0.1118374 -0.111837366
##     5 0.00000000  0.00000000  0.0000000  0.0000000  0.004129334
##                  5
##     0  0.000000000
##     1  0.000000000
##     2  0.000000000
##     3  0.000000000
##     4  0.000000000
##     5 -0.004129334
## 
##  model summary:
##  log-likelihood = -211.796485
##  AIC = 433.592970
##  AICc = 433.831065
##  free parameters = 5
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 0.06
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

We have also been able to fit this model using fitMk, but this appears to be broken:

fit.phytools<-fitMk(sqTree,toes,model=directional)

## Warning in log(comp[1:M + N]): NaNs produced

## Warning in log(comp[1:M + N]): NaNs produced

## Warning in log(comp[1:M + N]): NaNs produced

## Warning in log(comp[1:M + N]): NaNs produced

Let's load the fix and test it as follows:

source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/fitMk.R")
fit.phytools<-fitMk(sqTree,toes,model=directional)
fit.phytools

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##          0         1         2         3         4         5
## 0 0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
## 1 0.023147 -0.023147  0.000000  0.000000  0.000000  0.000000
## 2 0.000000  0.173751 -0.173751  0.000000  0.000000  0.000000
## 3 0.000000  0.000000  0.212357 -0.212357  0.000000  0.000000
## 4 0.000000  0.000000  0.000000  0.111838 -0.111838  0.000000
## 5 0.000000  0.000000  0.000000  0.000000  0.004129 -0.004129
## 
## Fitted (or set) value of pi:
##         0         1         2         3         4         5 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 
## 
## Log-likelihood: -213.588245 
## 
## Optimization method used was "nlminb"

Note that though very similar, fitMk and fitDiscrete are slightly different because they use different assumptions about the prior at the root of the tree. Nonetheless, it seems that fitMk is fixed and working again.

The bug will also affect other functions that use fitMk internally, such as make.simmap, but it can be fixed by updating phytools from GitHub.