Since 2019
*phytools* has contained a function (`fitpolyMk`

) for fitting a polymorphic trait
evolution model. This is a relatively simple model in which the polymorphic condition
(e.g., *a* + *b*) is *intermediate* between the two, corresponding monomorphic states.
`fitpolyMk`

has no problem with higher degrees of polymorphism in which, for example,
the condition *a* + *b* + *c* might be intermediate between *a* + *b* and *b* + *c*.

I’ve presented on this model several times and I get lots of questions about it, so I thought it might be worthwhile to illustrate some ways in which the model & function can be used.

For this example I’m going to use simulated data for an imaginary clade of lizards in
which individual species might be ground-dwelling, bush-dwelling, tree-dwelling, and
perhaps various combinations of these three fundamental states. It might be the case
that evolutionary transitions occur ground ↔ bush ↔ tree (with all the
corresponding polymorphic conditions in between) – but we’re going to imagine that we
don’t know this *a priori*.

At the end I’ll reveal how the data were, in fact, simulated.

Note that to *fully* follow along with this post (as of the time of writing) you’ll
need to update *phytools* from GitHub. This
can be done with (for example) using the
*remotes* CRAN package.

```
packageVersion("phytools")
```

```
## [1] '1.4.1'
```

```
library(phytools)
```

So, here’s our data.

```
plotTree(lizard.tree,ftype="off",lwd=1,type="fan")
X<-strsplit(setNames(as.character(habs),names(habs)),"+",
fixed=TRUE)
pies<-matrix(0,Ntip(lizard.tree),3,
dimnames=list(lizard.tree$tip.label,
c("ground","bush","tree")))
for(i in 1:Ntip(lizard.tree))
pies[lizard.tree$tip.label[i],
X[[lizard.tree$tip.label[i]]]]<-
rep(1/length(X[[lizard.tree$tip.label[i]]]),
length(X[[lizard.tree$tip.label[i]]]))
cols<-c("gold1","olivedrab3","darkgreen")
par(fg="transparent")
tiplabels(pie=pies,piecol=cols,cex=0.3)
par(fg="black")
legend(x="topleft",legend=c("ground","bush","tree"),
pt.cex=2,pch=16,col=cols)
```

(Someone who’s not colorblind can tell me whether or not that’s a good palette.)

We can see from this plot that each tip of the tree is in the conditions *ground*, *bush*,
*tree*, or some combination of states.

To start, let’s examine the endeavor of **picking the best model**.

When we do this, we can think of two broad *flavors* of polymorphic trait evolution model:
one in which the monomorphic conditions have some kind of implied *order*; and a second
in which they do not.

A good example of an *ordered* polymorphic character might be a meristic trait: such at
the number of scales at the midline or the number of caudal vertebrae.

As a small sidebar, it’s probably worth noting that the *most general* unordered
polymorphic trait evolution model (that is, `model="ARD"`

) is guaranteed to have a higher
likelihood than all *ordered* models, because it has them as a special case in which the
transition rates between certain states or state combinations are zero.

Here, I’ll fit just *three* variants of the unordered polymorphic trait evolution: `"ER"`

(all transition rates equal);
`"transient"`

;
and `"ARD"`

.

```
unordered.er<-fitpolyMk(lizard.tree,habs,model="ER")
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## bush ground tree bush+ground bush+tree ground+tree bush+ground+tree
## bush 0 0 0 1 1 0 0
## ground 0 0 0 1 0 1 0
## tree 0 0 0 0 1 1 0
## bush+ground 1 1 0 0 0 0 1
## bush+tree 1 0 1 0 0 0 1
## ground+tree 0 1 1 0 0 0 1
## bush+ground+tree 0 0 0 1 1 1 0
```

```
unordered.transient<-fitpolyMk(lizard.tree,habs,
model="transient")
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## bush ground tree bush+ground bush+tree ground+tree bush+ground+tree
## bush 0 0 0 2 2 0 0
## ground 0 0 0 2 0 2 0
## tree 0 0 0 0 2 2 0
## bush+ground 1 1 0 0 0 0 2
## bush+tree 1 0 1 0 0 0 2
## ground+tree 0 1 1 0 0 0 2
## bush+ground+tree 0 0 0 1 1 1 0
```

```
unordered.ard<-fitpolyMk(lizard.tree,habs,model="ARD")
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## bush ground tree bush+ground bush+tree ground+tree bush+ground+tree
## bush 0 0 0 1 3 0 0
## ground 0 0 0 5 0 7 0
## tree 0 0 0 0 9 11 0
## bush+ground 2 6 0 0 0 0 13
## bush+tree 4 0 10 0 0 0 15
## ground+tree 0 8 12 0 0 0 17
## bush+ground+tree 0 0 0 14 16 18 0
```

```
layout(matrix(c(1,1,2,2,4,3,3,4),4,2))
plot(unordered.er)
mtext(paste("a) unordered ER, AIC:",
round(AIC(unordered.er),2)),adj=0)
plot(unordered.transient)
mtext(paste("b) unordered transient, AIC:",
round(AIC(unordered.transient),2)),adj=0)
plot(unordered.ard)
mtext(paste("c) unordered ARD, AIC:",
round(AIC(unordered.ard),2)),adj=0)
```

It’s evident which model is best-supported by our data among these three, but we can
also compare them using a new `anova`

method for the model class as follows.

```
anova(unordered.er,unordered.transient,unordered.ard)
```

```
## log(L) d.f. AIC weight
## unordered.er -265.1739 1 532.3479 1.106000e-09
## unordered.transient -243.6807 2 491.3613 8.787669e-01
## unordered.ard -229.6615 18 495.3229 1.212331e-01
```

This tells us fairly definitively that among the three models we’ve evaluated so far, the so-called “transient” model (featuring one rate for the acquisition of polymorphism and a second rate for its loss) is the best-supported by our data.

The other major flavor of polymorphic trait evolution that we’ll consider is the *ordered*
model. We can specify this model class by setting the `fitpolyMk`

argument `ordered=TRUE`

as follows.

```
ordered.er<-fitpolyMk(lizard.tree,habs,model="ER",ordered=TRUE)
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## bush bush+ground bush+ground+tree ground ground+tree tree
## bush 0 1 0 0 0 0
## bush+ground 1 0 1 1 0 0
## bush+ground+tree 0 1 0 0 1 0
## ground 0 1 0 0 1 0
## ground+tree 0 0 1 1 0 1
## tree 0 0 0 0 1 0
```

```
ordered.er
```

```
## Object of class "fitpolyMk".
##
## Evolution was modeled as 'ordered' (i.e., transitions are assumed
## to occur bush <-> bush+ground <-> bush+ground+tree and so on) using the "ER" model.
##
## Fitted (or set) value of Q:
## bush bush+ground bush+ground+tree ground ground+tree tree
## bush NaN NaN 0 0 0 0
## bush+ground NaN NaN NaN NaN 0 0
## bush+ground+tree 0 NaN NaN 0 NaN 0
## ground 0 NaN 0 NaN NaN 0
## ground+tree 0 0 NaN NaN NaN NaN
## tree 0 0 0 0 NaN NaN
##
## Fitted (or set) value of pi:
## bush bush+ground bush+ground+tree ground ground+tree
## 0.166667 0.166667 0.166667 0.166667 0.166667
## tree
## 0.166667
##
## Log-likelihood: -1e+50
##
## Optimization method used was "nlminb"
```

OK, something’s out-of-place here. When we take a closer at both our fitted model and our data we see the following mismatch.

```
ordered.er$states
```

```
## [1] "bush" "bush+ground" "bush+ground+tree" "ground"
## [5] "ground+tree" "tree"
```

```
levels(habs)
```

```
## [1] "bush" "bush+tree" "ground" "ground+bush"
## [5] "ground+bush+tree" "tree"
```

The default ordering of our monomorphic condition has created a set of polymorphic conditions
that is *different* from the set of conditions that we observe in our data!

Is there an ordered process that would produce the conditions that we observe? Yes, in fact, and it is the one we hypothesized at the beginning of this article: ground ↔ bush ↔ tree (with all the polymorphisms in between, of course).

This specific ordering *hypothesis* can be specified using the `fitpolyMk`

argument `order`

(as opposed to `ordered`

, which is a logical value), as follows.

```
levs<-c("ground","bush","tree")
ordered.er<-fitpolyMk(lizard.tree,habs,model="ER",
ordered=TRUE,order=levs)
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## ground bush+ground bush+ground+tree bush bush+tree tree
## ground 0 1 0 0 0 0
## bush+ground 1 0 1 1 0 0
## bush+ground+tree 0 1 0 0 1 0
## bush 0 1 0 0 1 0
## bush+tree 0 0 1 1 0 1
## tree 0 0 0 0 1 0
```

```
ordered.transient<-fitpolyMk(lizard.tree,habs,
model="transient",ordered=TRUE,order=levs)
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## ground bush+ground bush+ground+tree bush bush+tree tree
## ground 0 2 0 0 0 0
## bush+ground 1 0 2 1 0 0
## bush+ground+tree 0 1 0 0 1 0
## bush 0 2 0 0 2 0
## bush+tree 0 0 2 1 0 1
## tree 0 0 0 0 2 0
```

```
ordered.ard<-fitpolyMk(lizard.tree,habs,model="ARD",
ordered=TRUE,order=levs)
```

```
##
## This is the design matrix of the fitted model. Does it make sense?
##
## ground bush+ground bush+ground+tree bush bush+tree tree
## ground 0 1 0 0 0 0
## bush+ground 2 0 3 5 0 0
## bush+ground+tree 0 4 0 0 7 0
## bush 0 6 0 0 9 0
## bush+tree 0 0 8 10 0 11
## tree 0 0 0 0 12 0
```

Let’s graph these three fitted models as we did before.

```
layout(matrix(c(1,1,2,2,4,3,3,4),4,2))
plot(ordered.er)
mtext(paste("a) ordered ER, AIC:",
round(AIC(ordered.er),2)),adj=0)
plot(ordered.transient)
mtext(paste("b) ordered transient, AIC:",
round(AIC(ordered.transient),2)),adj=0)
plot(ordered.ard)
mtext(paste("c) ordered ARD, AIC:",
round(AIC(ordered.ard),2)),adj=0)
```

As before, we can also compare all of these models using an `anova`

call; however, in this
case we’ll also include our first set of three *unordered* models, as follows.

```
anova(ordered.er,unordered.er,
ordered.transient,unordered.transient,
ordered.ard,unordered.ard)
```

```
## log(L) d.f. AIC weight
## ordered.er -248.6455 1 499.2910 7.034847e-07
## unordered.er -265.1739 1 532.3479 4.667033e-14
## ordered.transient -233.4797 2 470.9594 9.985836e-01
## unordered.transient -243.6807 2 491.3613 3.708169e-05
## ordered.ard -230.0686 12 484.1373 1.373518e-03
## unordered.ard -229.6615 18 495.3229 5.115723e-06
```

With the entire set of models in hand, it’s now evident that the best-supported model in the six is clearly the ordered transient model.

For a next step, let’s try to *use* our fitted model to undertake **ancestral state
reconstruction**.

For this, I prefer the CRAN package *corHMM*;
however, to use *corHMM* we’re going to first have to do some bookkeeping with our data and
fitted model.

To see this, let’s first load the package.

```
library(corHMM)
```

Next, I’m going to pull my states off the fitted model object. It’s important that I pull
them from there because this ensures that (say) `a+b`

and `b+a`

have been synonymized (as
they should be).

```
xx<-apply(ordered.transient$data,1,
function(x,ss) ss[which(x==1)],
ss=colnames(ordered.transient$data))
```

I also discovered that *corHMM* does not like the `+`

character in state names. No problem.
We’ll use `gsub`

to subsitute `/`

.

```
xx<-gsub("+","/",xx,fixed=TRUE)
```

Now let’s build the data frame that `corHMM`

needs as input.

```
lizard.data<-data.frame(Genus_sp=names(xx),
habitat=xx)
head(lizard.data)
```

```
## Genus_sp habitat
## t220 t220 tree
## t221 t221 tree
## t199 t199 tree
## t36 t36 tree
## t37 t37 tree
## t81 t81 tree
```

So far so good. Next, I’ll pull out my model design matrix from the fitted object. Here, I
need to replace all the zeroes with `NA`

, rename the rows & columns substituting `/`

for
`+`

as I did in my state data, *and* reorder alphabetically.

```
rate.mat<-ordered.transient$index.matrix
rate.mat[rate.mat==0]<-NA
colnames(rate.mat)<-rownames(rate.mat)<-gsub("+","/",
colnames(rate.mat),fixed=TRUE)
ind<-order(colnames(rate.mat))
rate.mat<-rate.mat[ind,ind]
rate.mat
```

```
## bush bush/ground bush/ground/tree bush/tree ground tree
## bush NA 2 NA 2 NA NA
## bush/ground 1 NA 2 NA 1 NA
## bush/ground/tree NA 1 NA 1 NA NA
## bush/tree 1 NA 2 NA NA 1
## ground NA 2 NA NA NA NA
## tree NA NA NA 2 NA NA
```

Finally, I’m ready to get my ancestral states!

```
fit.marginal<-corHMM(lizard.tree,lizard.data,
rate.mat=rate.mat,
node.states="marginal",
rate.cat=1,p=ordered.transient$rates,
root.p=ordered.transient$pi)
```

```
## State distribution in data:
## States: 1 2 3 4 5 6
## Counts: 91 33 5 36 14 121
## Calculating likelihood from a set of fixed parameters
## Finished. Inferring ancestral states using marginal reconstruction.
```

```
fit.marginal
```

```
##
## Fit
## -lnL AIC AICc Rate.cat ntax
## -233.4797 470.9594 470.9998 1 300
##
## Legend
## 1 2 3 4
## "bush" "bush/ground" "bush/ground/tree" "bush/tree"
## 5 6
## "ground" "tree"
##
## Rates
## (1,R1) (2,R1) (3,R1) (4,R1) (5,R1) (6,R1)
## (1,R1) NA 0.6716599 NA 0.6716599 NA NA
## (2,R1) 2.056186 NA 0.6716599 NA 2.056186 NA
## (3,R1) NA 2.0561858 NA 2.0561858 NA NA
## (4,R1) 2.056186 NA 0.6716599 NA NA 2.056186
## (5,R1) NA 0.6716599 NA NA NA NA
## (6,R1) NA NA NA 0.6716599 NA NA
##
## Arrived at a reliable solution
```

I should see that the likelihood matches what I obtained from `fitpolyMk`

, which is a good
sign.

The marginal states are stored in `fit.marginal$states`

matrix. Let’s pull that out and
do some our previous bookkeeping again, but this time in reverse.

```
asr<-fit.marginal$states
colnames(asr)<-colnames(rate.mat)
colnames(asr)<-gsub("/","+",colnames(asr))
head(asr)
```

```
## bush bush+ground bush+ground+tree bush+tree ground tree
## [1,] 0.173153538 9.473746e-03 0.1274649528 0.67980046 1.212651e-05 0.01009518
## [2,] 0.258401494 9.076941e-03 0.0865099912 0.61484627 1.184951e-05 0.03115345
## [3,] 0.279675949 8.692735e-03 0.0607896743 0.59017966 1.823457e-05 0.06064374
## [4,] 0.045162574 2.810016e-04 0.0055312345 0.37299114 2.453109e-07 0.57603380
## [5,] 0.004614478 1.572653e-05 0.0013535985 0.22013706 1.738418e-09 0.77387914
## [6,] 0.000263540 1.605630e-06 0.0002711173 0.04837474 4.856456e-10 0.95108900
```

Last of all, plot our reconstructed states on the tree.

```
plotTree(lizard.tree,ftype="off",lwd=1,type="fan")
X<-strsplit(setNames(as.character(habs),names(habs)),"+",
fixed=TRUE)
pies<-matrix(0,Ntip(lizard.tree),3,
dimnames=list(lizard.tree$tip.label,
c("ground","bush","tree")))
for(i in 1:Ntip(lizard.tree))
pies[lizard.tree$tip.label[i],
X[[lizard.tree$tip.label[i]]]]<-
rep(1/length(X[[lizard.tree$tip.label[i]]]),
length(X[[lizard.tree$tip.label[i]]]))
cols<-setNames(c("gold1","olivedrab3","darkgreen"),
c("ground","bush","tree"))
par(fg="transparent")
tiplabels(pie=pies,piecol=cols,cex=0.3)
par(fg="black")
piecol<-vector()
for(i in 1:ncol(asr)){
nn<-strsplit(colnames(asr)[i],"+",fixed=TRUE)[[1]]
if(length(nn)==1) piecol[i]<-cols[nn]
else if(length(nn)==2) piecol[i]<-colorRampPalette(cols[nn])(3)[2]
else piecol[i]<-"black"
}
names(piecol)<-colnames(asr)
par(fg="transparent")
nodelabels(pie=asr,piecol=piecol,cex=0.5)
par(fg="black")
legend(x="topleft",legend=colnames(asr),
pt.cex=1.8,pch=16,cex=0.8,col=piecol,
bty="n")
```

Cool.

For our final analysis, I’ll show how to apply our fitted model to the method of
**stochastic character mapping**. Since both `fitpolyMk`

and `make.simmap`

are both
from *phytools* this is going to be pretty easy.

```
smap.trees<-make.simmap(lizard.tree,ordered.transient$data,
Q=as.Qmatrix(ordered.transient),pi=ordered.transient$pi,
nsim=100)
```

```
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## ground bush+ground bush+ground+tree bush bush+tree
## ground -0.6716599 0.6716599 0.0000000 0.000000 0.0000000
## bush+ground 2.0561858 -4.7840314 0.6716599 2.056186 0.0000000
## bush+ground+tree 0.0000000 2.0561858 -4.1123715 0.000000 2.0561858
## bush 0.0000000 0.6716599 0.0000000 -1.343320 0.6716599
## bush+tree 0.0000000 0.0000000 0.6716599 2.056186 -4.7840314
## tree 0.0000000 0.0000000 0.0000000 0.000000 0.6716599
## tree
## ground 0.0000000
## bush+ground 0.0000000
## bush+ground+tree 0.0000000
## bush 0.0000000
## bush+tree 2.0561858
## tree -0.6716599
## (specified by the user);
## and (mean) root node prior probabilities
## pi =
## ground bush+ground bush+ground+tree bush bush+tree
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## tree
## 0.1666667
```

```
## Done.
```

```
smap.trees
```

```
## 100 phylogenetic trees with mapped discrete characters
```

For fun, let’s look at our whole set of stochastic mapped trees.

```
par(mfrow=c(10,10))
plot(smap.trees,ftype="off",lwd=1,colors=piecol,type="fan")
```