As readers of this blog can probably imagine, I get *a lot* of emails inquiring on how to do this analysis or another. Sometimes these are using *phytools*, other times with *geiger* or other R phylogenetic comparative methods packages, and sometimes they even pertain to phylogeny inference (which I know a bit about, for sure – but have not actually ever really worked on). As often as I can I’ll do my best to respond, and sometimes I’ll even pass along a chapter or two from my recent book with Luke Harmon.

A common inquiry regards *ancestral state reconstruction*.

This is fair. I recently wrote a book chapter on this topic (which can be read in pre-print form here), so it might be reasonable to call me a bit of an “expert” (though a professed ‘ancestral state’ skeptic, at the same time).

In the book for ancestral state reconstruction of *continuous* traits we use `phytools::fastAnc`

. Even though there are other implementations out there, `fastAnc`

will still be a good option for many users & even provides the correct variances & 95% confidence intervals based on Rohlf (2001).

For discrete characters, on the other hand, we use the powerful R package *corHMM*, by Jeremy Beaulieu & al. (including James Boyko, whom I recently had the pleasure of co-teaching an R phylogenetic workshop with). *corHMM* works great; however, I have since *also* implemented marginal and joint ancestral state reconstruction for discrete characters in *phytools* via the new generic method `ancr`

.

There’s no *inherent* reason to prefer `phytools::ancr`

over *corHMM* – but if you read this blog, you’re probably already quite used to running code with *phytools*, so that’s one.

Unfortunately, we published the book before I’d implemented `ancr`

, so it’s not covered there. Since I can’t point to the book for a simple demo of how ancestral state estimation of discrete traits works using `phytools::ancr`

, I thought I might as well make one here!

(There is *one* published demo of `phytools::ancr`

in my monographic software note for *phytools* 2.0 in the open-access journal *PeerJ*.)

Here I’m going to provide not a single tutorial, but a rapid-fire set of demos: first of marginal ASR for a single M*k* transition model, then marginal ASR integrating over multiple models using model averaging, marginal ASR when some tips are not known with certainty, then marginal ASR extensions of the M*k* model, such as the hidden-rates model, and finally joint reconstruction.

Pretty much all of these methods (except for model-averaging, but this could easily be done ‘manually’) are also implemented in *corHMM*.

Before we begin, let’s load *phytools* and its dependencies, which we’ll be using throughout this tutorial.

```
## load packages
library(phytools)
```

### Marginal ancestral state reconstruction for a single Q matrix

For this first demo, I’ll use a dataset of squamate digit number based on Brandley et al. (2008). These data can be downloaded or read directly from the book website. Here I’ll do the latter, as follows.

```
squamate.tree<-read.nexus(
file="http://www.phytools.org/Rbook/6/squamate.tre")
squamate.tree
```

```
##
## Phylogenetic tree with 258 tips and 257 internal nodes.
##
## Tip labels:
## Abronia_graminea, Acontias_litoralis, Acontias_meleagris, Acontias_percivali, Acontophiops_lineatus, Acrochordus_granulatus, ...
##
## Rooted; includes branch lengths.
```

```
squamate.data<-read.csv(
file="http://www.phytools.org/Rbook/6/squamate-data.csv",
row.names=1)
head(squamate.data)
```

```
## rear.toes
## Acontias_meleagris 0
## Acontias_percivali 0
## Alopoglossus_atriventris 5
## Alopoglossus_copii 5
## Ameiva_ameiva 5
## Amphiglossus_igneocaudatus 5
```

There are some mismatches between our tree and data, so we can start by fixing these (which I’ll do using the still handy *geiger* function `name.check`

).

```
chk<-geiger::name.check(squamate.tree,
squamate.data)
summary(chk)
```

```
## 139 taxa are present in the tree but not the data:
## Abronia_graminea,
## Acontias_litoralis,
## Acontophiops_lineatus,
## Acrochordus_granulatus,
## Agamodon_anguliceps,
## Agkistrodon_contortrix,
## ....
## 1 taxon is present in the data but not the tree:
## Trachyboa_boulengeri
##
## To see complete list of mis-matched taxa, print object.
```

```
## prune tree
squamate.tree<-drop.tip(squamate.tree,
chk$tree_not_data)
## subsample data
squamate.data<-squamate.data[squamate.tree$tip.label,,
drop=FALSE]
## check again
geiger::name.check(squamate.tree,
squamate.data)
```

```
## [1] "OK"
```

Let’s pull out the data vector of our trait of interest. We’ll always need to do this; however, for ancestral state reconstruction of a discrete character our input data can be a named factor, character, or integer vector (or numeric *matrix*, more on that below).

```
squamate.toes<-setNames(
squamate.data$rear.toes,
rownames(squamate.data))
head(squamate.toes)
```

```
## Acontias_meleagris Acontias_percivali Alopoglossus_atriventris
## 0 0 5
## Alopoglossus_copii Ameiva_ameiva Amphiglossus_igneocaudatus
## 5 5 5
```

Next, I’m going to hypothesize a trait evolution model for these data & fit it. This model is one in which digits of the rear foot are invariably gained or loss in a sequential fashion: than is, to evolve from having 4 digits in the hindfoot to only two, you must first pass through the state of having 3 hindfoot digits. I won’t pre-suppose that these rates are equal among each other and certainly not in the two (loss vs. gain) directions. This seems pretty logical as an approximation.

To fit this particular model to our data, however, I first need to *build* it using a *design matrix* of the same dimensions as my trait space (that is \(k \times k\) for \(k\) levels of the trait), but containing 0s and positive integers, depending on the types of state changes that we want to include in our model and those that we don’t.

```
ordered_model<-matrix(c(
0,1,0,0,0,0,
2,0,3,0,0,0,
0,4,0,5,0,0,
0,0,6,0,7,0,
0,0,0,8,0,9,
0,0,0,0,10,0),6,6,
byrow=TRUE,
dimnames=list(0:5,0:5))
ordered_model
```

```
## 0 1 2 3 4 5
## 0 0 1 0 0 0 0
## 1 2 0 3 0 0 0
## 2 0 4 0 5 0 0
## 3 0 0 6 0 7 0
## 4 0 0 0 8 0 9
## 5 0 0 0 0 10 0
```

Since this model contains quite a few parameters (ten, in fact – though it could be worse), let’s run multiple optimization iterations to try to ensure that we converge to the true Maximum Likelihood solution.

I’ll do this in parallel using the *foreach* package as I demonstrated earlier this month on the blog. (I’m also going to use `try`

to make sure that if one of our optimization fails, we don’t lose all of them!)

```
library(foreach)
library(doParallel)
niter<-10 ## set iterations
## set ncores and open cluster
ncores<-min(niter,detectCores()-1)
mc<-makeCluster(ncores,type="PSOCK")
registerDoParallel(cl=mc)
all_fits<-foreach(i=1:niter)%dopar%{
obj<-list()
class(obj)<-"try-error"
while(inherits(obj,"try-error")){
obj<-try(phytools::fitMk(squamate.tree,
squamate.toes,model=ordered_model,pi="fitzjohn",
logscale=sample(c(FALSE,TRUE),1),
opt.method=sample(c("nlminb","optim"),1),
rand.start=TRUE))
}
obj
}
stopCluster(mc) ## stop cluster
lnL<-sapply(all_fits,logLik)
lnL
```

```
## [1] -114.9672 -115.2785 -114.9672 -115.2785 -114.9672 -115.2785 -115.2785 -114.9672 -115.2785
## [10] -114.9672
```

Here’s our best-fitting ordered model.

```
squamate_fit<-all_fits[[which.max(lnL)]]
plot(squamate_fit,width=TRUE,color=TRUE,mar=rep(0.1,4),
xlim=c(-1.8,1),ylim=c(-1.1,1.1),text=FALSE,
show.zeros=FALSE,max.lwd=6)
```

We shouldn’t worry too much about the *very high* transition rates between character levels 2 and 3. These just tell us that these two trait levels are probably virtually randomly distributed *with respect to one another* on the tree. That is to say, if a species has state 3 than it is just as likely that a very close relative will have state 2 as share state 3 and *vice versa*. Indeed, I predict that the same will be true of our estimated ancestral states.

With our fitted model in hand, we’re ready to run `ancr`

as follows.

```
squamate_ancr<-ancr(squamate_fit)
squamate_ancr
```

```
## Marginal ancestral state estimates:
## 0 1 2 3 4 5
## 120 0 0 0 0 0.125952 0.874048
## 121 0 0 0 0 0.027224 0.972776
## 122 0 0 0 0 0.027248 0.972752
## 123 0 0 0 0 0.041905 0.958095
## 124 0 0 0 0 0.074453 0.925547
## 125 0 0 0 0 0.063773 0.936227
## ...
##
## Log-likelihood = -114.967167
```

Very conveniently for us, the `"ancr"`

object class also comes with an accompanying `plot`

method. Here’s what that can look like.

```
cols<-hcl.colors(n=6)
plot(squamate_ancr,args.nodelabels=list(piecol=cols),
args.tiplabels=list(piecol=cols),direction="upwards")
```

Great. That’s not bad!

### Model-averaged ancestral states

For this next demo, I’ll use a phylogeny of primates and diel activity pattern from Kirk & Kay (2004). The data & tree are packaged with *phytools*, so there’s no need to download the files from the web (though they’re available from our book website).

```
## load tree and data
data("primate.tree")
primate.tree
```

```
##
## Phylogenetic tree with 90 tips and 89 internal nodes.
##
## Tip labels:
## Allenopithecus_nigroviridis, Cercopithecus_mitis, Cercopithecus_cephus, Cercopithecus_petaurista, Erythrocebus_patas, Chlorocebus_aethiops, ...
##
## Rooted; includes branch lengths.
```

```
data("primate.data")
head(primate.data)
```

```
## Group Skull_length Optic_foramen_area Orbit_area
## Allenopithecus_nigroviridis Anthropoid 98.5 7.0 298.7
## Alouatta_palliata Anthropoid 109.8 5.3 382.3
## Alouatta_seniculus Anthropoid 108.0 8.0 359.4
## Aotus_trivirgatus Anthropoid 60.5 3.1 297.4
## Arctocebus_aureus Prosimian 49.5 1.2 134.8
## Arctocebus_calabarensis Prosimian 53.8 1.6 156.6
## Activity_pattern Activity_pattern_code
## Allenopithecus_nigroviridis Diurnal 0
## Alouatta_palliata Diurnal 0
## Alouatta_seniculus Diurnal 0
## Aotus_trivirgatus Nocturnal 2
## Arctocebus_aureus Nocturnal 2
## Arctocebus_calabarensis Nocturnal 2
```

Next, just as before, let’s pull out our data as a single named vector as follows.

```
activity<-setNames(primate.data$Activity_pattern,
rownames(primate.data))
head(activity)
```

```
## Allenopithecus_nigroviridis Alouatta_palliata Alouatta_seniculus
## Diurnal Diurnal Diurnal
## Aotus_trivirgatus Arctocebus_aureus Arctocebus_calabarensis
## Nocturnal Nocturnal Nocturnal
## Levels: Cathemeral Diurnal Nocturnal
```

Our three character states are "cathemeral" (active randomly during the day and night", "diurnal," and "nocturnal."

Now we’re ready to fit our models. Even though I might have also contemplated other custom models, as we did for the squamate digits previously, for simplicity here I’ll fit just a model with equal back & forth transition rates between each pair of states (the ER model), a model with symmetric rates (the SYM model), and a model that allows all transition rates between different states to assume different values (the ARD model), as follows.

```
primate_er<-fitMk(primate.tree,activity,
model="ER",pi="fitzjohn")
primate_sym<-fitMk(primate.tree,activity,
model="SYM",pi="fitzjohn")
primate_ard<-fitMk(primate.tree,activity,
model="ARD",pi="fitzjohn")
```

After having fit each of these three models, let’s compare them using a generic `anova`

call. Take note that I’m saving the output of `anova`

to a new object called `primate_aov`

.

```
primate_aov<-anova(primate_er,
primate_sym,primate_ard)
```

```
## log(L) d.f. AIC weight
## primate_er -29.49184 1 60.98368 0.3348288
## primate_sym -27.69681 3 61.39363 0.2727746
## primate_ard -24.33319 6 60.66637 0.3923967
```

Our analysis reveals that the strength of evidence supporting each of these three different models is quite *similar* one to the other. As such, we might choose to use model averaging in our ancestral reconstruction. What this entails is just computing a set of weighted average ancestral states, in which the weights come from our Akaike weights. `ancr`

will do this *automatically*, in fact, if we pass it our `primate_aov`

object instead of any one of our original fitted models.

```
primate_ancr<-ancr(primate_aov)
primate_ancr
```

```
## Marginal ancestral state estimates:
## Cathemeral Diurnal Nocturnal
## 91 0.000018 0.011355 0.988626
## 92 0.002987 0.034054 0.962960
## 93 0.059849 0.881411 0.058740
## 94 0.010596 0.988241 0.001162
## 95 0.001634 0.998335 0.000030
## 96 0.000227 0.999767 0.000006
## ...
##
## Log-likelihood = -27.106777
```

We can finish this demo by visualizing the results we obtained.

```
cols<-setNames(c("goldenrod","lightblue","black"),
levels(activity))
node.cex<-apply(primate_ancr$ace,1,
function(x) if(any(x>0.95)) 0.3 else 0.6)
plot(primate_ancr,
args.plotTree=list(type="arc",arc_height=0.5,
fsize=0.5,offset=3),
args.nodelabels=list(cex=node.cex,piecol=cols),
args.tiplabels=list(cex=0.2,piecol=cols),
legend=FALSE)
legend(0.2*max(nodeHeights(primate.tree)),
0.3*max(nodeHeights(primate.tree)),
levels(activity),pch=16,col=cols,
horiz=FALSE,cex=0.8,bty="n",pt.cex=2,
y.intersp=1.2)
```

### Ancestral state reconstruction when tip values are uncertain

Something else that *phytools* can handle easily is ancestral state estimation when the tip values for some or many species are uncertain.

I’ve written about this extensively in the past, most comprehensively here.

In summary, we can either treat tip states as if the ambiguity was observed directly or as if we observe one condition for the character, but some probability exists that we might be incorrect. More details on the differences and logic underlying this distinction is in my previous post on this topic.

Let’s return to our primate example from before and imagine, for the purposes of this demo, that our diel activity behavioral data for primates is highly imperfect. Not only that, but it is *asymmetrically* imperfect: any species that is observed to be cathemeral (active randomly during the day and night) definitely *is* cathemeral; but species that are only observed as active during the day are probably diurnal, but may be cathemeral, with a similar logic applying to nocturnally active taxa.

Now, of course, plenty is known about primate diel activity pattern – so I’m in no way asserting this as a genuine attribute of our data; however, this property can & does often affect real interspecific datasets (like our data in urban tolerance in *Anolis* lizard from Winchell et al. 2020).

To start, I’ll just convert my factor vector `activity`

to a binary matrix using `phytools::to.matrix`

as follows.

```
Activity<-to.matrix(activity,levels(activity))
head(Activity)
```

```
## Cathemeral Diurnal Nocturnal
## Allenopithecus_nigroviridis 0 1 0
## Alouatta_palliata 0 1 0
## Alouatta_seniculus 0 1 0
## Aotus_trivirgatus 0 0 1
## Arctocebus_aureus 0 0 1
## Arctocebus_calabarensis 0 0 1
```

Next, for each species, I’ll (arbitrarily) determine that if the species is *coded* as `Diurnal`

or `Nocturnal`

it might actually be `Cathemeral`

with, say, a 10% probability.

```
for(i in 1:nrow(Activity)){
Activity[i,]<-if(Activity[i,"Diurnal"]==1) c(0.1,0.9,0) else
if(Activity[i,"Nocturnal"]==1) c(0.1,0,0.9) else
c(1,0,0)
}
head(Activity)
```

```
## Cathemeral Diurnal Nocturnal
## Allenopithecus_nigroviridis 0.1 0.9 0.0
## Alouatta_palliata 0.1 0.9 0.0
## Alouatta_seniculus 0.1 0.9 0.0
## Aotus_trivirgatus 0.1 0.0 0.9
## Arctocebus_aureus 0.1 0.0 0.9
## Arctocebus_calabarensis 0.1 0.0 0.9
```

This is what that looks like.

Now, just as before, we can run our models and estimate ancestral states.

```
primate_er2<-fitMk(primate.tree,Activity,
model="ER",pi="fitzjohn")
primate_sym2<-fitMk(primate.tree,Activity,
model="SYM",pi="fitzjohn")
primate_ard2<-fitMk(primate.tree,Activity,
model="ARD",pi="fitzjohn")
primate_aov2<-anova(primate_er2,primate_sym2,
primate_ard2)
```

```
## log(L) d.f. AIC weight
## primate_er2 -36.15873 1 74.31746 0.6472514
## primate_sym2 -35.32473 3 76.64947 0.2016893
## primate_ard2 -32.61379 6 77.22758 0.1510593
```

```
primate_ancr2<-ancr(primate_ard2)
primate_ancr2
```

```
## Marginal ancestral state estimates:
## Cathemeral Diurnal Nocturnal
## 91 0.000008 0.000537 0.999455
## 92 0.006294 0.014579 0.979127
## 93 0.156004 0.787791 0.056205
## 94 0.026345 0.972644 0.001011
## 95 0.003722 0.996256 0.000022
## 96 0.000495 0.999502 0.000004
## ...
##
## Log-likelihood = -32.613791
```

That’s all there is to it!

### Ancestral state reconstruction under the hidden-rates model

The *phytools* generic method `ancr`

is designed to work not only with `fitMk`

model objects, but also with many of the other discrete character evolution model object that *phytools* can produce. These include fitted models from `fitgammaMk`

(for \(\Gamma\) distributed rate heterogeneity among edges), `fitfnMk`

(for an ordered discrete character evolution model where the rate of evolution between adjacent states varies as a polynomial functional form), `fitpolyMk`

(for data with intraspecific polymorphism), and `fitHRM`

(for the hidden-rates model).

To illustrate ancestral state reconstruction under just one of these, the hidden-rates model, I’m going to use a dataset for pollen grain number in angiosperms from Williams et al. (2014). This is the same dataset that we used in our book and I wrote a previous post here comparing the hidden-rates and \(\Gamma\) models for these data. The Williams et al. (2014) can be downloaded from our book webpage, or read in directly from the web as I’ll do below.

First let’s get the data into R using `read.csv`

.

```
pollen_data<-read.csv(
file="http://www.phytools.org/Rbook/7/pollen-data.csv",
row.names=1)
head(pollen_data)
```

```
## V1
## Nuphar_lutea 3
## Nuphar_pumila 2
## Nuphar_advena 2
## Austrobaileya_scandens 2
## Schisandra_propinqua 2
## Kadsura_heteroclita 2
```

Next the tree, of course.

```
pollen_tree<-read.tree(
file="http://www.phytools.org/Rbook/7/pollen-tree.phy")
pollen_tree
```

```
##
## Phylogenetic tree with 511 tips and 510 internal nodes.
##
## Tip labels:
## Nuphar_lutea, Nuphar_pumila, Nuphar_advena, Austrobaileya_scandens, Schisandra_propinqua, Kadsura_heteroclita, ...
##
## Rooted; includes branch lengths.
```

As always, our next step must be to pull our discrete character out as a new vector with names corresponding to the tip labels of our tree.

```
pollen<-setNames(pollen_data$V1,
rownames(pollen_data))
```

To fit my hidden-rates model, I’m going to use a model with just two hidden levels. In the book I believe that we followed Williams et al. & used a model with three levels. This would take a lot longer to run, so for demonstrative purposes only (and with no assertion that this is a good model for our data) I’ll stick to the two-level model here.

```
pollen_hrm<-fitHRM(pollen_tree,pollen,ncat=2,
pi="fitzjohn",parallel=TRUE,niter=20)
```

```
##
## This is the design matrix of the fitted model.
## Does it make sense?
##
## 2 2* 3 3*
## 2 0 1 2 0
## 2* 3 0 0 4
## 3 5 0 0 6
## 3* 0 7 8 0
##
## Opened cluster with 15 cores.
## Running optimization iterations in parallel.
## Please wait....
```

```
pollen_hrm
```

```
## Object of class "fitHRM".
##
## Observed states: [ 2, 3 ]
## Number of rate categories per state: [ 2, 2 ]
##
## Fitted (or set) value of Q:
## 2 2* 3 3*
## 2 -0.090800 0.002171 0.088629 0.000000
## 2* 0.002895 -0.002895 0.000000 0.000000
## 3 0.961639 0.000000 -0.999650 0.038012
## 3* 0.000000 0.000000 0.003129 -0.003129
##
## Fitted (or set) value of pi:
## 2 2* 3 3*
## 0.108009 0.785332 0.101532 0.005126
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -201.944088
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
```

Finally, obtaining marginal ancestral states is done in the typical way, as follows.

```
pollen_ancr<-ancr(pollen_hrm)
pollen_ancr
```

```
## Marginal ancestral state estimates:
## 2 2* 3 3*
## 512 0.018264 0.965556 0.016139 0.000041
## 513 0.635040 0.249619 0.061181 0.054160
## 514 0.655555 0.273008 0.057236 0.014201
## 515 0.017118 0.981453 0.001399 0.000030
## 516 0.061269 0.932561 0.005148 0.001021
## 517 0.251698 0.727173 0.021026 0.000103
## ...
##
## Log-likelihood = -201.944088
```

```
cols<-c(hcl.colors(6)[c(1,2,4,5)],
c("2","2*","3","3*"))
plot(pollen_ancr,
args.plotTree=list(type="arc",arc_height=0.5),
args.nodelabels=list(piecol=cols,cex=0.25),
args.tiplabels=list(piecol=cols))
```

Terrific! That was easy.

### Joint ancestral state reconstruction

Marginal ancestral state reconstruction, which we’ve been doing so far, has the goal of computing the marginal scaled likelihoods that each node is in each of the observed (or hidden) states of a character, integrating of all states at all other nodes of the phylogeny. If we want to make explicit, probabilistic statements about *particular* nodes, then this is the ancestral state reconstruction method we should choose.

Joint ancestral state reconstruction, on the other hand, has a different aim. It involves identifying the set of node values at *all* nodes for which the probability of our observed data under the model (the likelihood) is the highest. Joint reconstruction is less useful for making explicit probabilistic statements about individual nodes, but is better if we want to make an observation about how our character is most likely to have evolved node by node in the tree.

To illustrate joint ancestral state reconstruction in *phytools*, I’ll use a dataset for bony fish obtained from Benun Sutton & Wilson (2019) but that modifies a phylogeny that was originally published by Betancur et al. (2017). These files are available from the book website, but are also packaged with *phytools* and thus most easily loaded directly as follows.

```
data(bonyfish.tree)
bonyfish.tree
```

```
##
## Phylogenetic tree with 90 tips and 89 internal nodes.
##
## Tip labels:
## Xenomystus_nigri, Chirocentrus_dorab, Talismania_bifurcata, Alepocephalus_tenebrosus, Misgurnus_bipartitus, Opsariichthys_bidens, ...
##
## Rooted; includes branch lengths.
```

```
data(bonyfish.data)
head(bonyfish.data)
```

```
## spawning_mode paternal_care
## Xenomystus_nigri pair male
## Chirocentrus_dorab group none
## Talismania_bifurcata group none
## Alepocephalus_tenebrosus group none
## Misgurnus_bipartitus pair none
## Opsariichthys_bidens pair none
```

For our analysis, we’ll model and reconstruct the history of the character of ‘spawning mode’, which has two different values in our dataset: `"pair"`

or `"group"`

.

```
spawning_mode<-setNames(
bonyfish.data$spawning_mode,
rownames(bonyfish.data))
head(spawning_mode)
```

```
## Xenomystus_nigri Chirocentrus_dorab Talismania_bifurcata
## pair group group
## Alepocephalus_tenebrosus Misgurnus_bipartitus Opsariichthys_bidens
## group pair pair
## Levels: group pair
```

```
bonyfish_er<-fitMk(bonyfish.tree,spawning_mode,
model="ER",pi="fitzjohn")
bonyfish_ard<-fitMk(bonyfish.tree,spawning_mode,
model="ARD",pi="fitzjohn")
anova(bonyfish_er,bonyfish_ard)
```

```
## log(L) d.f. AIC weight
## bonyfish_er -29.51819 1 61.03638 0.6726777
## bonyfish_ard -29.23851 2 62.47702 0.3273223
```

This tells us that, of these two models, the ER model is better-supported compared to its more parameter-rich counterpart. To undertake joint reconstruction, we just run the following code on our chosen fitted model.

```
bonyfish_joint<-ancr(bonyfish_er,type="joint")
bonyfish_joint
```

```
## Joint ancestral state estimates:
## state
## 91 pair
## 92 pair
## 93 pair
## 94 pair
## 95 group
## 96 pair
## ...
##
## Log-likelihood = -30.54223
```

Just for interest, let’s *also* perform marginal reconstruction under the same model, and compare our results.

```
bonyfish_marginal<-ancr(bonyfish_er,
type="marginal")
```

```
par(mfrow=c(2,1))
obj<-plot(bonyfish_joint,direction="upwards",
legend=FALSE)
legend("bottomright",levels(spawning_mode),
pch=16,col=obj$piecol,bty="n",cex=0.8,
pt.cex=1.5)
obj<-plot(bonyfish_marginal,direction="upwards",
legend=FALSE)
legend("bottomright",levels(spawning_mode),
pch=16,col=obj$piecol,bty="n",cex=0.8,
pt.cex=1.5)
```