Monday, July 6, 2020

A new polymorphic ordered discrete character evolution model for phytools

Last year I developed a new method in phytools (fitpolyMk) to fit the Mk model to discrete character data with intraspecific polymorphism.

The idea is pretty simple - we just build a model in which the polymorphic state (say, 0 + 1) is an intermediate between the two corresponding monomorphic conditions (0 and 1).

Things get a tiny bit more complicated if we have more than two states for our discrete character, because in that case the evolution of our character might be ordered (i.e., 00 + 111 + 2 and so on) or unordered (e.g., 00 + 20 + 1 + 21 + 2 and so on).

Without thinking too much about it, however, I supposed that ordered evolution automatically implied that a lineage could only have a maximum of two states for the character. I later realized that this is only one type of ordered character evolution: 00 + 10 + 1 + 21 + 2 is also ordered. In fact, this latter form of ordered discrete character evolution may be common for some traits - such as meristic characters.

As of today, I have added this model to fitpolyMk. It can be specified by updating the optional argument max.states which presently defaults to 2, but can be given any value up to and including the number of distinct states of the character.

This update can be obtained from GitHub with the devtools package installed:

library(devtools)
install_github("liamrevell/phytools")
library(phytools)
packageVersion("phytools")
## [1] '0.7.54'

Let's try it out.

First, here are some data.

## plot our data:
plotTree(tree,ftype="off",lwd=1,type="fan")
X<-strsplit(setNames(as.character(x),names(x)),"+",
    fixed=TRUE)
pies<-matrix(0,Ntip(tree),5,dimnames=list(tree$tip.label,
    0:4))
for(i in 1:Ntip(tree)) 
    pies[tree$tip.label[i],X[[tree$tip.label[i]]]]<-
        rep(1/length(X[[tree$tip.label[i]]]),
        length(X[[tree$tip.label[i]]]))
tiplabels(pie=pies,piecol=palette()[1:5],
    cex=0.35)
legend(x="topleft",legend=0:4,pt.cex=2,pch=21,
    pt.bg=palette()[1:5])

plot of chunk unnamed-chunk-3

Let's fit a polymorphic character model to these data. Here, I'll use the 'transient' model in which polymorphism is acquired with one rate & then lost with another:

fit.ordered<-fitpolyMk(tree,x,ordered=TRUE,
    max.states=5,pi="fitzjohn",model="transient")
## 
## This is the design matrix of the fitted model. Does it make sense?
## 
##           0 0+1 0+1+2 0+1+2+3 0+1+2+3+4 1 1+2 1+2+3 1+2+3+4 2 2+3 2+3+4 3 3+4
## 0         0   2     0       0         0 0   0     0       0 0   0     0 0   0
## 0+1       1   0     2       0         0 1   0     0       0 0   0     0 0   0
## 0+1+2     0   1     0       2         0 0   1     0       0 0   0     0 0   0
## 0+1+2+3   0   0     1       0         2 0   0     1       0 0   0     0 0   0
## 0+1+2+3+4 0   0     0       1         0 0   0     0       1 0   0     0 0   0
## 1         0   2     0       0         0 0   2     0       0 0   0     0 0   0
## 1+2       0   0     2       0         0 1   0     2       0 1   0     0 0   0
## 1+2+3     0   0     0       2         0 0   1     0       2 0   1     0 0   0
## 1+2+3+4   0   0     0       0         2 0   0     1       0 0   0     1 0   0
## 2         0   0     0       0         0 0   2     0       0 0   2     0 0   0
## 2+3       0   0     0       0         0 0   0     2       0 1   0     2 1   0
## 2+3+4     0   0     0       0         0 0   0     0       2 0   1     0 0   1
## 3         0   0     0       0         0 0   0     0       0 0   2     0 0   2
## 3+4       0   0     0       0         0 0   0     0       0 0   0     2 1   0
## 4         0   0     0       0         0 0   0     0       0 0   0     0 0   2
##           4
## 0         0
## 0+1       0
## 0+1+2     0
## 0+1+2+3   0
## 0+1+2+3+4 0
## 1         0
## 1+2       0
## 1+2+3     0
## 1+2+3+4   0
## 2         0
## 2+3       0
## 2+3+4     0
## 3         0
## 3+4       1
## 4         0
print(fit.ordered,digits=4)
## Object of class "fitpolyMk".
## 
## Evolution was modeled as 'ordered' (i.e., transitions are assumed
## to occur 0 <-> 0+1 <-> 0+1+2 and so on) using the "transient" model.
## 
## Fitted (or set) value of Q:
##                 0     0+1   0+1+2 0+1+2+3 0+1+2+3+4       1     1+2   1+2+3
## 0         -1.4243  1.4243  0.0000  0.0000    0.0000  0.0000  0.0000  0.0000
## 0+1        3.3043 -8.0328  1.4243  0.0000    0.0000  3.3043  0.0000  0.0000
## 0+1+2      0.0000  3.3043 -8.0328  1.4243    0.0000  0.0000  3.3043  0.0000
## 0+1+2+3    0.0000  0.0000  3.3043 -8.0328    1.4243  0.0000  0.0000  3.3043
## 0+1+2+3+4  0.0000  0.0000  0.0000  3.3043   -6.6085  0.0000  0.0000  0.0000
## 1          0.0000  1.4243  0.0000  0.0000    0.0000 -2.8486  1.4243  0.0000
## 1+2        0.0000  0.0000  1.4243  0.0000    0.0000  3.3043 -9.4571  1.4243
## 1+2+3      0.0000  0.0000  0.0000  1.4243    0.0000  0.0000  3.3043 -9.4571
## 1+2+3+4    0.0000  0.0000  0.0000  0.0000    1.4243  0.0000  0.0000  3.3043
## 2          0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  1.4243  0.0000
## 2+3        0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  0.0000  1.4243
## 2+3+4      0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  0.0000  0.0000
## 3          0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  0.0000  0.0000
## 3+4        0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  0.0000  0.0000
## 4          0.0000  0.0000  0.0000  0.0000    0.0000  0.0000  0.0000  0.0000
##           1+2+3+4       2     2+3   2+3+4       3     3+4       4
## 0          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3+4  3.3043  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2        0.0000  3.3043  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2+3      1.4243  0.0000  3.3043  0.0000  0.0000  0.0000  0.0000
## 1+2+3+4   -8.0328  0.0000  0.0000  3.3043  0.0000  0.0000  0.0000
## 2          0.0000 -2.8486  1.4243  0.0000  0.0000  0.0000  0.0000
## 2+3        0.0000  3.3043 -9.4571  1.4243  3.3043  0.0000  0.0000
## 2+3+4      1.4243  0.0000  3.3043 -8.0328  0.0000  3.3043  0.0000
## 3          0.0000  0.0000  1.4243  0.0000 -2.8486  1.4243  0.0000
## 3+4        0.0000  0.0000  0.0000  1.4243  3.3043 -8.0328  3.3043
## 4          0.0000  0.0000  0.0000  0.0000  0.0000  1.4243 -1.4243
## 
## Fitted (or set) value of pi:
##         0       0+1     0+1+2   0+1+2+3 0+1+2+3+4         1       1+2 
##    0.0003    0.0151    0.1388    0.1381    0.0453    0.0589    0.2475 
##     1+2+3   1+2+3+4         2       2+3     2+3+4         3       3+4 
##    0.1440    0.0344    0.1292    0.0393    0.0065    0.0025    0.0003 
##         4 
##    0.0000 
## 
## Log-likelihood: -305.0916 
## 
## Optimization method used was "nlminb"
par(mar=rep(0.1,4))
plot(fit.ordered)

plot of chunk unnamed-chunk-4

Since we never observe 0+1+2+3+4 in these data, we might decide to assume that the maximum number of states is 4 instead of 5. Let's set it to have this value & fit our model again:

fit.ordered.4state<-fitpolyMk(tree,x,ordered=TRUE,
    max.states=4,pi="fitzjohn",model="transient")
## 
## This is the design matrix of the fitted model. Does it make sense?
## 
##         0 0+1 0+1+2 0+1+2+3 1 1+2 1+2+3 1+2+3+4 2 2+3 2+3+4 3 3+4 4
## 0       0   2     0       0 0   0     0       0 0   0     0 0   0 0
## 0+1     1   0     2       0 1   0     0       0 0   0     0 0   0 0
## 0+1+2   0   1     0       2 0   1     0       0 0   0     0 0   0 0
## 0+1+2+3 0   0     1       0 0   0     1       0 0   0     0 0   0 0
## 1       0   2     0       0 0   2     0       0 0   0     0 0   0 0
## 1+2     0   0     2       0 1   0     2       0 1   0     0 0   0 0
## 1+2+3   0   0     0       2 0   1     0       2 0   1     0 0   0 0
## 1+2+3+4 0   0     0       0 0   0     1       0 0   0     1 0   0 0
## 2       0   0     0       0 0   2     0       0 0   2     0 0   0 0
## 2+3     0   0     0       0 0   0     2       0 1   0     2 1   0 0
## 2+3+4   0   0     0       0 0   0     0       2 0   1     0 0   1 0
## 3       0   0     0       0 0   0     0       0 0   2     0 0   2 0
## 3+4     0   0     0       0 0   0     0       0 0   0     2 1   0 1
## 4       0   0     0       0 0   0     0       0 0   0     0 0   2 0
print(fit.ordered.4state,digits=4)
## Object of class "fitpolyMk".
## 
## Evolution was modeled as 'ordered' (i.e., transitions are assumed
## to occur 0 <-> 0+1 <-> 0+1+2 and so on) using the "transient" model.
## 
## Fitted (or set) value of Q:
##               0     0+1   0+1+2 0+1+2+3       1     1+2   1+2+3 1+2+3+4
## 0       -1.4206  1.4206  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1      3.3658 -8.1523  1.4206  0.0000  3.3658  0.0000  0.0000  0.0000
## 0+1+2    0.0000  3.3658 -8.1523  1.4206  0.0000  3.3658  0.0000  0.0000
## 0+1+2+3  0.0000  0.0000  3.3658 -6.7317  0.0000  0.0000  3.3658  0.0000
## 1        0.0000  1.4206  0.0000  0.0000 -2.8411  1.4206  0.0000  0.0000
## 1+2      0.0000  0.0000  1.4206  0.0000  3.3658 -9.5728  1.4206  0.0000
## 1+2+3    0.0000  0.0000  0.0000  1.4206  0.0000  3.3658 -9.5728  1.4206
## 1+2+3+4  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  3.3658 -6.7317
## 2        0.0000  0.0000  0.0000  0.0000  0.0000  1.4206  0.0000  0.0000
## 2+3      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  1.4206  0.0000
## 2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  1.4206
## 3        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 4        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
##               2     2+3   2+3+4       3     3+4       4
## 0        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2      3.3658  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2+3    0.0000  3.3658  0.0000  0.0000  0.0000  0.0000
## 1+2+3+4  0.0000  0.0000  3.3658  0.0000  0.0000  0.0000
## 2       -2.8411  1.4206  0.0000  0.0000  0.0000  0.0000
## 2+3      3.3658 -9.5728  1.4206  3.3658  0.0000  0.0000
## 2+3+4    0.0000  3.3658 -8.1523  0.0000  3.3658  0.0000
## 3        0.0000  1.4206  0.0000 -2.8411  1.4206  0.0000
## 3+4      0.0000  0.0000  1.4206  3.3658 -8.1523  3.3658
## 4        0.0000  0.0000  0.0000  0.0000  1.4206 -1.4206
## 
## Fitted (or set) value of pi:
##       0     0+1   0+1+2 0+1+2+3       1     1+2   1+2+3 1+2+3+4       2 
##  0.0002  0.0139  0.1457  0.2205  0.0493  0.2189  0.1568  0.0395  0.1088 
##     2+3   2+3+4       3     3+4       4 
##  0.0371  0.0066  0.0023  0.0003  0.0000 
## 
## Log-likelihood: -304.1081 
## 
## Optimization method used was "nlminb"
plot(fit.ordered.4state)

plot of chunk unnamed-chunk-5

Finally, we could compare this model to an unordered transient model:

fit.unordered<-fitpolyMk(tree,x,pi="fitzjohn",
    model="transient")
## 
## This is the design matrix of the fitted model. Does it make sense?
## 
##           0 1 2 3 4 0+1 0+2 0+3 0+4 1+2 1+3 1+4 2+3 2+4 3+4 0+1+2 0+1+3 0+1+4
## 0         0 0 0 0 0   2   2   2   2   0   0   0   0   0   0     0     0     0
## 1         0 0 0 0 0   2   0   0   0   2   2   2   0   0   0     0     0     0
## 2         0 0 0 0 0   0   2   0   0   2   0   0   2   2   0     0     0     0
## 3         0 0 0 0 0   0   0   2   0   0   2   0   2   0   2     0     0     0
## 4         0 0 0 0 0   0   0   0   2   0   0   2   0   2   2     0     0     0
## 0+1       1 1 0 0 0   0   0   0   0   0   0   0   0   0   0     2     2     2
## 0+2       1 0 1 0 0   0   0   0   0   0   0   0   0   0   0     2     0     0
## 0+3       1 0 0 1 0   0   0   0   0   0   0   0   0   0   0     0     2     0
## 0+4       1 0 0 0 1   0   0   0   0   0   0   0   0   0   0     0     0     2
## 1+2       0 1 1 0 0   0   0   0   0   0   0   0   0   0   0     2     0     0
## 1+3       0 1 0 1 0   0   0   0   0   0   0   0   0   0   0     0     2     0
## 1+4       0 1 0 0 1   0   0   0   0   0   0   0   0   0   0     0     0     2
## 2+3       0 0 1 1 0   0   0   0   0   0   0   0   0   0   0     0     0     0
## 2+4       0 0 1 0 1   0   0   0   0   0   0   0   0   0   0     0     0     0
## 3+4       0 0 0 1 1   0   0   0   0   0   0   0   0   0   0     0     0     0
## 0+1+2     0 0 0 0 0   1   1   0   0   1   0   0   0   0   0     0     0     0
## 0+1+3     0 0 0 0 0   1   0   1   0   0   1   0   0   0   0     0     0     0
## 0+1+4     0 0 0 0 0   1   0   0   1   0   0   1   0   0   0     0     0     0
## 0+2+3     0 0 0 0 0   0   1   1   0   0   0   0   1   0   0     0     0     0
## 0+2+4     0 0 0 0 0   0   1   0   1   0   0   0   0   1   0     0     0     0
## 0+3+4     0 0 0 0 0   0   0   1   1   0   0   0   0   0   1     0     0     0
## 1+2+3     0 0 0 0 0   0   0   0   0   1   1   0   1   0   0     0     0     0
## 1+2+4     0 0 0 0 0   0   0   0   0   1   0   1   0   1   0     0     0     0
## 1+3+4     0 0 0 0 0   0   0   0   0   0   1   1   0   0   1     0     0     0
## 2+3+4     0 0 0 0 0   0   0   0   0   0   0   0   1   1   1     0     0     0
## 0+1+2+3   0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     1     1     0
## 0+1+2+4   0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     1     0     1
## 0+1+3+4   0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     0     1     1
## 0+2+3+4   0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     0     0     0
## 1+2+3+4   0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     0     0     0
## 0+1+2+3+4 0 0 0 0 0   0   0   0   0   0   0   0   0   0   0     0     0     0
##           0+2+3 0+2+4 0+3+4 1+2+3 1+2+4 1+3+4 2+3+4 0+1+2+3 0+1+2+4 0+1+3+4
## 0             0     0     0     0     0     0     0       0       0       0
## 1             0     0     0     0     0     0     0       0       0       0
## 2             0     0     0     0     0     0     0       0       0       0
## 3             0     0     0     0     0     0     0       0       0       0
## 4             0     0     0     0     0     0     0       0       0       0
## 0+1           0     0     0     0     0     0     0       0       0       0
## 0+2           2     2     0     0     0     0     0       0       0       0
## 0+3           2     0     2     0     0     0     0       0       0       0
## 0+4           0     2     2     0     0     0     0       0       0       0
## 1+2           0     0     0     2     2     0     0       0       0       0
## 1+3           0     0     0     2     0     2     0       0       0       0
## 1+4           0     0     0     0     2     2     0       0       0       0
## 2+3           2     0     0     2     0     0     2       0       0       0
## 2+4           0     2     0     0     2     0     2       0       0       0
## 3+4           0     0     2     0     0     2     2       0       0       0
## 0+1+2         0     0     0     0     0     0     0       2       2       0
## 0+1+3         0     0     0     0     0     0     0       2       0       2
## 0+1+4         0     0     0     0     0     0     0       0       2       2
## 0+2+3         0     0     0     0     0     0     0       2       0       0
## 0+2+4         0     0     0     0     0     0     0       0       2       0
## 0+3+4         0     0     0     0     0     0     0       0       0       2
## 1+2+3         0     0     0     0     0     0     0       2       0       0
## 1+2+4         0     0     0     0     0     0     0       0       2       0
## 1+3+4         0     0     0     0     0     0     0       0       0       2
## 2+3+4         0     0     0     0     0     0     0       0       0       0
## 0+1+2+3       1     0     0     1     0     0     0       0       0       0
## 0+1+2+4       0     1     0     0     1     0     0       0       0       0
## 0+1+3+4       0     0     1     0     0     1     0       0       0       0
## 0+2+3+4       1     1     1     0     0     0     1       0       0       0
## 1+2+3+4       0     0     0     1     1     1     1       0       0       0
## 0+1+2+3+4     0     0     0     0     0     0     0       1       1       1
##           0+2+3+4 1+2+3+4 0+1+2+3+4
## 0               0       0         0
## 1               0       0         0
## 2               0       0         0
## 3               0       0         0
## 4               0       0         0
## 0+1             0       0         0
## 0+2             0       0         0
## 0+3             0       0         0
## 0+4             0       0         0
## 1+2             0       0         0
## 1+3             0       0         0
## 1+4             0       0         0
## 2+3             0       0         0
## 2+4             0       0         0
## 3+4             0       0         0
## 0+1+2           0       0         0
## 0+1+3           0       0         0
## 0+1+4           0       0         0
## 0+2+3           2       0         0
## 0+2+4           2       0         0
## 0+3+4           2       0         0
## 1+2+3           0       2         0
## 1+2+4           0       2         0
## 1+3+4           0       2         0
## 2+3+4           2       2         0
## 0+1+2+3         0       0         2
## 0+1+2+4         0       0         2
## 0+1+3+4         0       0         2
## 0+2+3+4         0       0         2
## 1+2+3+4         0       0         2
## 0+1+2+3+4       1       1         0
print(fit.unordered,digits=4)
## Object of class "fitpolyMk".
## 
## Evolution was modeled as 'unordered' using the "transient" model.
## 
## Fitted (or set) value of Q:
##                 0       1       2       3       4     0+1     0+2     0+3
## 0         -1.6103  0.0000  0.0000  0.0000  0.0000  0.4026  0.4026  0.4026
## 1          0.0000 -1.6103  0.0000  0.0000  0.0000  0.4026  0.0000  0.0000
## 2          0.0000  0.0000 -1.6103  0.0000  0.0000  0.0000  0.4026  0.0000
## 3          0.0000  0.0000  0.0000 -1.6103  0.0000  0.0000  0.0000  0.4026
## 4          0.0000  0.0000  0.0000  0.0000 -1.6103  0.0000  0.0000  0.0000
## 0+1        2.0372  2.0372  0.0000  0.0000  0.0000 -5.2821  0.0000  0.0000
## 0+2        2.0372  0.0000  2.0372  0.0000  0.0000  0.0000 -5.2821  0.0000
## 0+3        2.0372  0.0000  0.0000  2.0372  0.0000  0.0000  0.0000 -5.2821
## 0+4        2.0372  0.0000  0.0000  0.0000  2.0372  0.0000  0.0000  0.0000
## 1+2        0.0000  2.0372  2.0372  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+3        0.0000  2.0372  0.0000  2.0372  0.0000  0.0000  0.0000  0.0000
## 1+4        0.0000  2.0372  0.0000  0.0000  2.0372  0.0000  0.0000  0.0000
## 2+3        0.0000  0.0000  2.0372  2.0372  0.0000  0.0000  0.0000  0.0000
## 2+4        0.0000  0.0000  2.0372  0.0000  2.0372  0.0000  0.0000  0.0000
## 3+4        0.0000  0.0000  0.0000  2.0372  2.0372  0.0000  0.0000  0.0000
## 0+1+2      0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  2.0372  0.0000
## 0+1+3      0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  0.0000  2.0372
## 0+1+4      0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  0.0000  0.0000
## 0+2+3      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  2.0372
## 0+2+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  0.0000
## 0+3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372
## 1+2+3      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 2+3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3+4  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
##               0+4     1+2     1+3     1+4     2+3     2+4     3+4   0+1+2
## 0          0.4026  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1          0.0000  0.4026  0.4026  0.4026  0.0000  0.0000  0.0000  0.0000
## 2          0.0000  0.4026  0.0000  0.0000  0.4026  0.4026  0.0000  0.0000
## 3          0.0000  0.0000  0.4026  0.0000  0.4026  0.0000  0.4026  0.0000
## 4          0.4026  0.0000  0.0000  0.4026  0.0000  0.4026  0.4026  0.0000
## 0+1        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.4026
## 0+2        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.4026
## 0+3        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+4       -5.2821  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2        0.0000 -5.2821  0.0000  0.0000  0.0000  0.0000  0.0000  0.4026
## 1+3        0.0000  0.0000 -5.2821  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+4        0.0000  0.0000  0.0000 -5.2821  0.0000  0.0000  0.0000  0.0000
## 2+3        0.0000  0.0000  0.0000  0.0000 -5.2821  0.0000  0.0000  0.0000
## 2+4        0.0000  0.0000  0.0000  0.0000  0.0000 -5.2821  0.0000  0.0000
## 3+4        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -5.2821  0.0000
## 0+1+2      0.0000  2.0372  0.0000  0.0000  0.0000  0.0000  0.0000 -6.9167
## 0+1+3      0.0000  0.0000  2.0372  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+4      2.0372  0.0000  0.0000  2.0372  0.0000  0.0000  0.0000  0.0000
## 0+2+3      0.0000  0.0000  0.0000  0.0000  2.0372  0.0000  0.0000  0.0000
## 0+2+4      2.0372  0.0000  0.0000  0.0000  0.0000  2.0372  0.0000  0.0000
## 0+3+4      2.0372  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  0.0000
## 1+2+3      0.0000  2.0372  2.0372  0.0000  2.0372  0.0000  0.0000  0.0000
## 1+2+4      0.0000  2.0372  0.0000  2.0372  0.0000  2.0372  0.0000  0.0000
## 1+3+4      0.0000  0.0000  2.0372  2.0372  0.0000  0.0000  2.0372  0.0000
## 2+3+4      0.0000  0.0000  0.0000  0.0000  2.0372  2.0372  2.0372  0.0000
## 0+1+2+3    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372
## 0+1+2+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  2.0372
## 0+1+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1+2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3+4  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
##             0+1+3   0+1+4   0+2+3   0+2+4   0+3+4   1+2+3   1+2+4   1+3+4
## 0          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 1          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 2          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 3          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 4          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1        0.4026  0.4026  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+2        0.0000  0.0000  0.4026  0.4026  0.0000  0.0000  0.0000  0.0000
## 0+3        0.4026  0.0000  0.4026  0.0000  0.4026  0.0000  0.0000  0.0000
## 0+4        0.0000  0.4026  0.0000  0.4026  0.4026  0.0000  0.0000  0.0000
## 1+2        0.0000  0.0000  0.0000  0.0000  0.0000  0.4026  0.4026  0.0000
## 1+3        0.4026  0.0000  0.0000  0.0000  0.0000  0.4026  0.0000  0.4026
## 1+4        0.0000  0.4026  0.0000  0.0000  0.0000  0.0000  0.4026  0.4026
## 2+3        0.0000  0.0000  0.4026  0.0000  0.0000  0.4026  0.0000  0.0000
## 2+4        0.0000  0.0000  0.0000  0.4026  0.0000  0.0000  0.4026  0.0000
## 3+4        0.0000  0.0000  0.0000  0.0000  0.4026  0.0000  0.0000  0.4026
## 0+1+2      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+3     -6.9167  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+4      0.0000 -6.9167  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+2+3      0.0000  0.0000 -6.9167  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+2+4      0.0000  0.0000  0.0000 -6.9167  0.0000  0.0000  0.0000  0.0000
## 0+3+4      0.0000  0.0000  0.0000  0.0000 -6.9167  0.0000  0.0000  0.0000
## 1+2+3      0.0000  0.0000  0.0000  0.0000  0.0000 -6.9167  0.0000  0.0000
## 1+2+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -6.9167  0.0000
## 1+3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -6.9167
## 2+3+4      0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
## 0+1+2+3    2.0372  0.0000  2.0372  0.0000  0.0000  2.0372  0.0000  0.0000
## 0+1+2+4    0.0000  2.0372  0.0000  2.0372  0.0000  0.0000  2.0372  0.0000
## 0+1+3+4    2.0372  2.0372  0.0000  0.0000  2.0372  0.0000  0.0000  2.0372
## 0+2+3+4    0.0000  0.0000  2.0372  2.0372  2.0372  0.0000  0.0000  0.0000
## 1+2+3+4    0.0000  0.0000  0.0000  0.0000  0.0000  2.0372  2.0372  2.0372
## 0+1+2+3+4  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
##             2+3+4 0+1+2+3 0+1+2+4 0+1+3+4 0+2+3+4 1+2+3+4 0+1+2+3+4
## 0          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 1          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 2          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 3          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 4          0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 0+1        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 0+2        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 0+3        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 0+4        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 1+2        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 1+3        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 1+4        0.0000  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 2+3        0.4026  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 2+4        0.4026  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 3+4        0.4026  0.0000  0.0000  0.0000  0.0000  0.0000    0.0000
## 0+1+2      0.0000  0.4026  0.4026  0.0000  0.0000  0.0000    0.0000
## 0+1+3      0.0000  0.4026  0.0000  0.4026  0.0000  0.0000    0.0000
## 0+1+4      0.0000  0.0000  0.4026  0.4026  0.0000  0.0000    0.0000
## 0+2+3      0.0000  0.4026  0.0000  0.0000  0.4026  0.0000    0.0000
## 0+2+4      0.0000  0.0000  0.4026  0.0000  0.4026  0.0000    0.0000
## 0+3+4      0.0000  0.0000  0.0000  0.4026  0.4026  0.0000    0.0000
## 1+2+3      0.0000  0.4026  0.0000  0.0000  0.0000  0.4026    0.0000
## 1+2+4      0.0000  0.0000  0.4026  0.0000  0.0000  0.4026    0.0000
## 1+3+4      0.0000  0.0000  0.0000  0.4026  0.0000  0.4026    0.0000
## 2+3+4     -6.9167  0.0000  0.0000  0.0000  0.4026  0.4026    0.0000
## 0+1+2+3    0.0000 -8.5514  0.0000  0.0000  0.0000  0.0000    0.4026
## 0+1+2+4    0.0000  0.0000 -8.5514  0.0000  0.0000  0.0000    0.4026
## 0+1+3+4    0.0000  0.0000  0.0000 -8.5514  0.0000  0.0000    0.4026
## 0+2+3+4    2.0372  0.0000  0.0000  0.0000 -8.5514  0.0000    0.4026
## 1+2+3+4    2.0372  0.0000  0.0000  0.0000  0.0000 -8.5514    0.4026
## 0+1+2+3+4  0.0000  2.0372  2.0372  2.0372  2.0372  2.0372  -10.1860
## 
## Fitted (or set) value of pi:
##         0         1         2         3         4       0+1       0+2 
##    0.0001    0.0005    0.0090    0.0001    0.0000    0.0020    0.0230 
##       0+3       0+4       1+2       1+3       1+4       2+3       2+4 
##    0.0005    0.0000    0.0680    0.0014    0.0001    0.0068    0.0009 
##       3+4     0+1+2     0+1+3     0+1+4     0+2+3     0+2+4     0+3+4 
##    0.0000    0.2596    0.0057    0.0004    0.0281    0.0031    0.0001 
##     1+2+3     1+2+4     1+3+4     2+3+4   0+1+2+3   0+1+2+4   0+1+3+4 
##    0.0821    0.0092    0.0002    0.0009    0.3796    0.0400    0.0009 
##   0+2+3+4   1+2+3+4 0+1+2+3+4 
##    0.0043    0.0124    0.0610 
## 
## Log-likelihood: -355.6462 
## 
## Optimization method used was "nlminb"
plot(fit.unordered)

plot of chunk unnamed-chunk-6

data.frame(
    transition_model=c("transient","transient","transient"),
    ordered=c("yes","yes","no"),
    max_states=c(5,4,NA),
    logLik=c(logLik(fit.ordered),logLik(fit.ordered.4state),
    logLik(fit.unordered)),
    k=c(attr(logLik(fit.ordered),"df"),
    attr(logLik(fit.ordered.4state),"df"),
    attr(logLik(fit.unordered),"df")),
    AIC=c(AIC(fit.ordered),AIC(fit.ordered.4state),
    AIC(fit.unordered)))
##   transition_model ordered max_states    logLik k      AIC
## 1        transient     yes          5 -305.0916 2 614.1832
## 2        transient     yes          4 -304.1081 2 612.2161
## 3        transient      no         NA -355.6462 2 715.2925

Clearly, the ordered transient model is best (although the 5 & 4 state models had almost identical fits). That's good, at least, because this is the model I used to simulate the data:

set.seed(77)
tree<-pbtree(n=200,scale=1)
## 5-state 'transient' ordered
Q<-matrix(0,15,15)
rownames(Q)<-colnames(Q)<-c(
    "0","0+1","0+1+2","0+1+2+3","0+1+2+3+4",
    "1","1+2","1+2+3","1+2+3+4",
    "2","2+3","2+3+4",
    "3","3+4",
    "4")
Q[1,2]<-Q[2,3]<-Q[3,4]<-Q[4,5]<-Q[6,2]<-Q[6,7]<-
    Q[7,3]<-Q[7,8]<-Q[8,4]<-Q[8,9]<-Q[9,5]<-
    Q[10,7]<-Q[10,11]<-Q[11,8]<-Q[11,12]<-
    Q[12,9]<-Q[13,11]<-Q[13,14]<-Q[14,12]<-
    Q[14,15]<-Q[15,14]<-1
Q[2,1]<-Q[2,6]<-Q[3,2]<-Q[3,7]<-Q[4,3]<-Q[4,8]<-
    Q[5,4]<-Q[5,9]<-Q[7,6]<-Q[7,10]<-Q[8,7]<-
    Q[8,11]<-Q[9,8]<-Q[9,12]<-Q[11,10]<-
    Q[11,13]<-Q[12,11]<-Q[12,14]<-Q[14,13]<-
    Q[14,15]<-5
diag(Q)<--rowSums(Q)
x<-sim.Mk(tree,Q)

2 comments:

  1. Is there a way to plot the node pies on the tree for fitpolymk, as one can do for fitmk?

    ReplyDelete
    Replies
    1. Yes. The fitpolyMk function returns an object that contains a design matrix that can be used in other functions (such as ace, make.simmap, or rerootingMethod) that estimate ancestral states. Try it!

      Delete

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