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.,
0 ↔ 0 + 1 ↔ 1 ↔ 1 + 2 and so on) or unordered (e.g., 0 ↔
0 + 2 ↔ 0 + 1 + 2 ↔ 1 + 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:
0 → 0 + 1 → 0 + 1 + 2 → 1 + 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])
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)
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)
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)
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)
Is there a way to plot the node pies on the tree for fitpolymk, as one can do for fitmk?
ReplyDeleteYes. 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