## Saturday, May 6, 2023

### Fitting discrete character evolution models & reconstructing marginal ancestral states with uncertainty: Two approaches

Many people ask me about handling uncertain tip states in both Mk (and extended Mk) model-fitting, and marginal ancestral state reconstruction, using Maximum Likelihood.

There are effectively two ways that we can handle uncertainty in tip states.

The first, which is what is typically done in phylogeny inference by ML, is to treat the uncertainty (and this is my characterization, which I hope to justify below) as if it (that is, the uncertainty) was observed directly. For example, we might have “observed” that a nucleotide at a particular locus is a pyrimidine, but we don’t know whether it’s a T or a C. In that case, we compute the probability of our data given any particular value for our transition matrix (Q) as if the nucleotide was a T; compute the probability as if it was a C; and then add them. If we find the value of Q that maximizes this probability, we’ve got the MLE.

The second is to act as if the uncertainty (again, my characterization) was not observed directly. Instead, we imagine having observed one or the other condition of the character, but that this observation was made with some known or estimated degree of error. That is to say, perhaps we have a sequencing technology in which every time a thymine (T) is scored at a site there is some probability (Îµ) that the true nucleotide is actually a cytosine (C) – meaning that the chances it’s a T is in fact 1 - Îµ. In this case, for any value of Q, we’d calculate the total probability of our data by multiplying the probability given that the state is T (the state that we observed) by 1 - Îµ (the chances that T is correct, given that it was observed); and then adding this quantity to the probability of data given that the nucleotide is in condition C, multiplied by Îµ (the probability that it’s true state is C, even though T was observed). Does this make sense?

The most concise description of these two scenarios that I’m aware of is given in Chapter 4 of Yang (2014). The former approach is what we use in phylogeny inference for ambiguity (e.g., Y for a T or C; N for missing data at a site; etc.); whereas the latter is used for sequencing error (e.g., we observed T, but the true nucleotide might be something other than T, with some probability that can be defined by the investigator).

To see how we might go about undertaking either of the approaches using phytools in R, let’s first see an example of a typical workflow that we might use to reconstruct ancestral states at internal nodes without ambiguity.

``````library(phytools)
``````

For our example, we’ll use a primate phylogeny & dataset from Kirk & Kay (2004).

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

Here, our input data, `activity`, is a factor vector without ambiguity: it could also be a character vector or integer numeric vector.

``````head(activity,20)
``````
``````## Allenopithecus_nigroviridis           Alouatta_palliata          Alouatta_seniculus           Aotus_trivirgatus
##                     Diurnal                     Diurnal                     Diurnal                   Nocturnal
##           Arctocebus_aureus     Arctocebus_calabarensis            Ateles_fusciceps            Ateles_geoffroyi
##                   Nocturnal                   Nocturnal                     Diurnal                     Diurnal
##             Ateles_paniscus               Avahi_laniger      Cacajao_melanocephalus           Callicebus_moloch
##                     Diurnal                   Nocturnal                     Diurnal                     Diurnal
##        Callicebus_torquatus           Callimico_goeldii        Callithrix_argentata          Callithrix_jacchus
##                     Diurnal                     Diurnal                     Diurnal                     Diurnal
##          Callithrix_pygmaea             Cebus_albifrons             Cebus_capucinus        Cercopithecus_cephus
##                     Diurnal                     Diurnal                     Diurnal                     Diurnal
## Levels: Cathemeral Diurnal Nocturnal
``````

For fun, and because in this instance we’re not really interested in saving our intermediate results, let’s use the R pipe operator to fit a discrete character (extended Mk) trait evolution model to these data, and then compute and visualize marginal ancestral states.

``````library(magrittr)
fitMk(primate.tree,activity,model="ARD",
pi="fitzjohn",logscale=TRUE,lik.func="pruning") %>%
ancr %>%
plot(args.plotTree=list(direction="upwards"))
``````

OK, now let’s imagine that instead of knowing each tip state with certainty, we have some degree of ambiguity about the diel activity pattern state for one terminal taxon or another. To handle this ambiguity in R, we first need to convert our data into a matrix instead of a character or factor vector. Absent ambiguity, this matrix looks as follows.

``````X<-to.matrix(activity,levels(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
## Ateles_fusciceps                     0       1         0
## Ateles_geoffroyi                     0       1         0
## Ateles_paniscus                      0       1         0
## Avahi_laniger                        0       0         1
## Cacajao_melanocephalus               0       1         0
## Callicebus_moloch                    0       1         0
## Callicebus_torquatus                 0       1         0
## Callimico_goeldii                    0       1         0
## Callithrix_argentata                 0       1         0
## Callithrix_jacchus                   0       1         0
## Callithrix_pygmaea                   0       1         0
## Cebus_albifrons                      0       1         0
## Cebus_capucinus                      0       1         0
## Cercopithecus_cephus                 0       1         0
``````

Now, to simulate ambiguity (but since our data in this example don’t actually contain any real ambiguity, so far as I know), I’ll just pick a set of random taxon labels and then (following our first approach), I’ll assign them values of 1 for each of the three states of our character.

``````labs<-sample(primate.tree\$tip.label,10)
labs
``````
``````##  [1] "Pongo_pygmaeus"        "Eulemur_rubriventer"   "Lepilemur_mustelinus"  "Hapalemur_griseus"
##  [5] "Nasalis_larvatus"      "Propithecus_verreauxi" "Cheirogaleus_major"    "Macaca_mulatta"
##  [9] "Hylobates_lar"         "Callimico_goeldii"
``````
``````X1<-to.matrix(activity,levels(activity))
X1[labs,]<-rep(1,ncol(X1))
``````
``````##                             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
## Ateles_fusciceps                     0       1         0
## Ateles_geoffroyi                     0       1         0
## Ateles_paniscus                      0       1         0
## Avahi_laniger                        0       0         1
## Cacajao_melanocephalus               0       1         0
## Callicebus_moloch                    0       1         0
## Callicebus_torquatus                 0       1         0
## Callimico_goeldii                    1       1         1
## Callithrix_argentata                 0       1         0
## Callithrix_jacchus                   0       1         0
## Callithrix_pygmaea                   0       1         0
## Cebus_albifrons                      0       1         0
## Cebus_capucinus                      0       1         0
## Cercopithecus_cephus                 0       1         0
``````

We can then proceed to fit our model, just substituting `X1` (in this instance) for our factor vector of before.

``````fit1<-fitMk(primate.tree,X1,model="ARD",pi="fitzjohn",
lik.func="pruning",logscale=TRUE)
fit1
``````
``````## Object of class "fitMk".
##
## Fitted (or set) value of Q:
##            Cathemeral   Diurnal Nocturnal
## Cathemeral  -0.038595  0.038595  0.000000
## Diurnal      0.000000 -0.002050  0.002050
## Nocturnal    0.002557  0.003940 -0.006497
##
## Fitted (or set) value of pi:
## Cathemeral    Diurnal  Nocturnal
##   0.002321   0.022704   0.974975
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -22.3311
##
## Optimization method used was "nlminb"
``````

Finally, let’s calculate and plot ancestral states for our trait. In this case, we can even include the marginal probabilities for our ambiguous tip nodes – which we can obtain just by including them in our pre-order, marginal ancestral state reconstruction tree traversal. (This is done if we set `tips=TRUE` in `phytools::ancr`.)

``````ace1<-ancr(fit1,tips=TRUE)
print(ace1,printlen=20)
``````
``````## Marginal ancestral state estimates:
##                             Cathemeral  Diurnal Nocturnal
## Allenopithecus_nigroviridis   0.000000 1.000000  0.000000
## Cercopithecus_mitis           0.000000 1.000000  0.000000
## Cercopithecus_cephus          0.000000 1.000000  0.000000
## Cercopithecus_petaurista      0.000000 1.000000  0.000000
## Erythrocebus_patas            0.000000 1.000000  0.000000
## Chlorocebus_aethiops          0.000000 1.000000  0.000000
## Miopithecus_talapoin          0.000000 1.000000  0.000000
## Macaca_nemestrina             0.000000 1.000000  0.000000
## Macaca_fascicularis           0.000000 1.000000  0.000000
## Macaca_mulatta                0.000091 0.989701  0.010208
## Colobus_guereza               0.000000 1.000000  0.000000
## Colobus_polykomos             0.000000 1.000000  0.000000
## Nasalis_larvatus              0.000218 0.980750  0.019032
## Pygathrix_nemaeus             0.000000 1.000000  0.000000
## Rhinopithecus_roxellana       0.000000 1.000000  0.000000
## Semnopithecus_entellus        0.000000 1.000000  0.000000
## Trachypithecus_cristatus      0.000000 1.000000  0.000000
## Presbytis_comata              0.000000 1.000000  0.000000
## Presbytis_melalophos          0.000000 1.000000  0.000000
## ...
##
## Log-likelihood = -22.3311
``````
``````plot(ace1,args.plotTree=list(direction="upwards"))
tips<-sapply(labs,function(x,y) which(y==x),y=primate.tree\$tip.label)
col=palette()[4])
``````

The blue arrows mark all the tips that we “ambiguated” in our previous step.

The other way to treat uncertainty in the tip states, which (recall) is more akin to sequencing error (following Yang, 2014), would be to code one state as 1 - (k - 1) × Îµ, and the others as Îµ. If Îµ is said to be 1/k (for k states) then this is equivalent to assigning each level of the trait a value of 1/k.

Let’s try it.

``````X2<-to.matrix(activity,levels(activity))
X2[labs,]<-rep(1/ncol(X2),ncol(X2))
``````
``````##                             Cathemeral   Diurnal Nocturnal
## Allenopithecus_nigroviridis  0.0000000 1.0000000 0.0000000
## Alouatta_palliata            0.0000000 1.0000000 0.0000000
## Alouatta_seniculus           0.0000000 1.0000000 0.0000000
## Aotus_trivirgatus            0.0000000 0.0000000 1.0000000
## Arctocebus_aureus            0.0000000 0.0000000 1.0000000
## Arctocebus_calabarensis      0.0000000 0.0000000 1.0000000
## Ateles_fusciceps             0.0000000 1.0000000 0.0000000
## Ateles_geoffroyi             0.0000000 1.0000000 0.0000000
## Ateles_paniscus              0.0000000 1.0000000 0.0000000
## Avahi_laniger                0.0000000 0.0000000 1.0000000
## Cacajao_melanocephalus       0.0000000 1.0000000 0.0000000
## Callicebus_moloch            0.0000000 1.0000000 0.0000000
## Callicebus_torquatus         0.0000000 1.0000000 0.0000000
## Callimico_goeldii            0.3333333 0.3333333 0.3333333
## Callithrix_argentata         0.0000000 1.0000000 0.0000000
## Callithrix_jacchus           0.0000000 1.0000000 0.0000000
## Callithrix_pygmaea           0.0000000 1.0000000 0.0000000
## Cebus_albifrons              0.0000000 1.0000000 0.0000000
## Cebus_capucinus              0.0000000 1.0000000 0.0000000
## Cercopithecus_cephus         0.0000000 1.0000000 0.0000000
``````
``````fit2<-fitMk(primate.tree,X2,model="ARD",pi="fitzjohn",
logscale=TRUE,lik.func="pruning")
fit2
``````
``````## Object of class "fitMk".
##
## Fitted (or set) value of Q:
##            Cathemeral   Diurnal Nocturnal
## Cathemeral  -0.038595  0.038595  0.000000
## Diurnal      0.000000 -0.002050  0.002050
## Nocturnal    0.002557  0.003940 -0.006497
##
## Fitted (or set) value of pi:
## Cathemeral    Diurnal  Nocturnal
##   0.002321   0.022704   0.974975
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -33.317223
##
## Optimization method used was "nlminb"
``````

Something that probably jumps out at us straight away is that our likelihood has changed – but the parameter estimates of our fitted model are identical to when we set all ambiguous tip states to be vectors of 1s!

Let’s see if the same is true of marginal ancestral states.

``````ace2<-ancr(fit2,tips=TRUE)
print(ace2,printlen=20)
``````
``````## Marginal ancestral state estimates:
##                             Cathemeral  Diurnal Nocturnal
## Allenopithecus_nigroviridis   0.000000 1.000000  0.000000
## Cercopithecus_mitis           0.000000 1.000000  0.000000
## Cercopithecus_cephus          0.000000 1.000000  0.000000
## Cercopithecus_petaurista      0.000000 1.000000  0.000000
## Erythrocebus_patas            0.000000 1.000000  0.000000
## Chlorocebus_aethiops          0.000000 1.000000  0.000000
## Miopithecus_talapoin          0.000000 1.000000  0.000000
## Macaca_nemestrina             0.000000 1.000000  0.000000
## Macaca_fascicularis           0.000000 1.000000  0.000000
## Macaca_mulatta                0.000091 0.989701  0.010208
## Colobus_guereza               0.000000 1.000000  0.000000
## Colobus_polykomos             0.000000 1.000000  0.000000
## Nasalis_larvatus              0.000218 0.980750  0.019032
## Pygathrix_nemaeus             0.000000 1.000000  0.000000
## Rhinopithecus_roxellana       0.000000 1.000000  0.000000
## Semnopithecus_entellus        0.000000 1.000000  0.000000
## Trachypithecus_cristatus      0.000000 1.000000  0.000000
## Presbytis_comata              0.000000 1.000000  0.000000
## Presbytis_melalophos          0.000000 1.000000  0.000000
## ...
##
## Log-likelihood = -33.317223
``````
``````plot(ace2,args.plotTree=list(direction="upwards"))
tips<-sapply(labs,function(x,y) which(y==x),y=primate.tree\$tip.label)
col=palette()[4])
``````

What we see is that even though we obtain different likelihoods from these two different approaches to uncertainty, the fitted models and marginal states end up the same either way.

Which is the correct way to treat uncertain tip values? Well, someone smarter than me will have to answer that – but an advantage of treating uncertainty in the way that “sequencing error” is dealt with in phylogeny inference is that we can easily encode partial ambiguity in the character state.

What do I mean by that?

Well, just taking this simple example, let’s think about what each of our primate diel activity pattern states mean.

We have the state “cathemeral” which means active randomly during the day or night. And then we have the states diurnal & nocturnal: active during the day or during the night, respectively.

But a species that’s observed being active during the night could be cathemeral (but just, by chance, never observed as active during the day); and, conversely, a species that has been observed in the day could be cathemeral (but, by chance, never observed as active during the night). (Arguably, the converse is not true: a species observed as active during both the day & night must be cathemeral!)

Let’s add this uncertainty to our input matrix as follows:

``````X3<-to.matrix(activity,levels(activity))
X3[,1]<-sapply(X3[,1], function(x) if(x==0) 0.1 else 1)
X3[,2]<-sapply(X3[,2], function(x) if(x==1) 0.9 else 0)
X3[,3]<-sapply(X3[,3], function(x) if(x==1) 0.9 else 0)
``````
``````##                              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
## Ateles_fusciceps                    0.1     0.9       0.0
## Ateles_geoffroyi                    0.1     0.9       0.0
## Ateles_paniscus                     0.1     0.9       0.0
## Avahi_laniger                       0.1     0.0       0.9
## Cacajao_melanocephalus              0.1     0.9       0.0
## Callicebus_moloch                   0.1     0.9       0.0
## Callicebus_torquatus                0.1     0.9       0.0
## Callimico_goeldii                   0.1     0.9       0.0
## Callithrix_argentata                0.1     0.9       0.0
## Callithrix_jacchus                  0.1     0.9       0.0
## Callithrix_pygmaea                  0.1     0.9       0.0
## Cebus_albifrons                     0.1     0.9       0.0
## Cebus_capucinus                     0.1     0.9       0.0
## Cercopithecus_cephus                0.1     0.9       0.0
## Cercopithecus_mitis                 0.1     0.9       0.0
## Cercopithecus_petaurista            0.1     0.9       0.0
## Cheirogaleus_major                  0.1     0.0       0.9
## Cheirogaleus_medius                 0.1     0.0       0.9
## Chiropotes_satanas                  0.1     0.9       0.0
## Chlorocebus_aethiops                0.1     0.9       0.0
## Colobus_guereza                     0.1     0.9       0.0
## Colobus_polykomos                   0.1     0.9       0.0
## Erythrocebus_patas                  0.1     0.9       0.0
## Eulemur_fulvus                      1.0     0.0       0.0
## Eulemur_macaco                      1.0     0.0       0.0
## Eulemur_mongoz                      1.0     0.0       0.0
## Eulemur_rubriventer                 1.0     0.0       0.0
## Euoticus_elegantulus                0.1     0.0       0.9
## Galago_alleni                       0.1     0.0       0.9
## Galago_matschiei                    0.1     0.0       0.9
## Galago_moholi                       0.1     0.0       0.9
## Galago_senegalensis                 0.1     0.0       0.9
## Galagoides_demidoff                 0.1     0.0       0.9
``````

(The specific values I used here are arbitrary – except inasmuch as our rows should sum to 1, and the values should correspond to our prior probabilities that each tip is in each of our three states!)

``````fit3<-fitMk(primate.tree,X3,model="ARD",pi="fitzjohn",
logscale=TRUE,lik.func="pruning")
fit3
``````
``````## Object of class "fitMk".
##
## Fitted (or set) value of Q:
##            Cathemeral   Diurnal Nocturnal
## Cathemeral  -0.042896  0.042896  0.000000
## Diurnal      0.000000 -0.001968  0.001968
## Nocturnal    0.003170  0.003010 -0.006180
##
## Fitted (or set) value of pi:
## Cathemeral    Diurnal  Nocturnal
##   0.002789   0.022582   0.974629
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -32.613791
##
## Optimization method used was "nlminb"
``````
``````ace3<-ancr(fit3,tips=TRUE)
print(ace3,printlen=20)
``````
``````## Marginal ancestral state estimates:
##                             Cathemeral  Diurnal Nocturnal
## Allenopithecus_nigroviridis    5.4e-05 0.999946         0
## Cercopithecus_mitis            1.5e-05 0.999985         0
## Cercopithecus_cephus           4.0e-06 0.999996         0
## Cercopithecus_petaurista       4.0e-06 0.999996         0
## Erythrocebus_patas             2.6e-05 0.999974         0
## Chlorocebus_aethiops           2.6e-05 0.999974         0
## Miopithecus_talapoin           4.9e-05 0.999951         0
## Macaca_nemestrina              3.1e-05 0.999969         0
## Macaca_fascicularis            1.9e-05 0.999981         0
## Macaca_mulatta                 1.9e-05 0.999981         0
## Colobus_guereza                1.5e-05 0.999985         0
## Colobus_polykomos              1.5e-05 0.999985         0
## Nasalis_larvatus               3.2e-05 0.999968         0
## Pygathrix_nemaeus              3.2e-05 0.999968         0
## Rhinopithecus_roxellana        3.6e-05 0.999964         0
## Semnopithecus_entellus         4.5e-05 0.999955         0
## Trachypithecus_cristatus       4.5e-05 0.999955         0
## Presbytis_comata               1.9e-05 0.999981         0
## Presbytis_melalophos           1.9e-05 0.999981         0
## ...
##
## Log-likelihood = -32.613791
``````
``````plot(ace3,args.plotTree=list(direction="upwards"))
``````

Sensibly, since these are pretty-much the same* data, we obtain predictably similar results as before.

One would think that we need to decide a priori how uncertainty is being treated across all the tips of the tree and be consistent; however, surprisingly, one would be wrong.

To see that, let’s combine our uncertainty in the code of diurnal & nocturnal species (above) with the total ambiguity of a subset of tips that we imagined earlier, but using our two different methods (all 1s or 1/k for each level) of encoding that total ambiguity.

First, the former.

``````X4<-X3
X4[labs,]<-rep(1,ncol(X4))
``````
``````##                             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
## Ateles_fusciceps                   0.1     0.9       0.0
## Ateles_geoffroyi                   0.1     0.9       0.0
## Ateles_paniscus                    0.1     0.9       0.0
## Avahi_laniger                      0.1     0.0       0.9
## Cacajao_melanocephalus             0.1     0.9       0.0
## Callicebus_moloch                  0.1     0.9       0.0
## Callicebus_torquatus               0.1     0.9       0.0
## Callimico_goeldii                  1.0     1.0       1.0
## Callithrix_argentata               0.1     0.9       0.0
## Callithrix_jacchus                 0.1     0.9       0.0
## Callithrix_pygmaea                 0.1     0.9       0.0
## Cebus_albifrons                    0.1     0.9       0.0
## Cebus_capucinus                    0.1     0.9       0.0
## Cercopithecus_cephus               0.1     0.9       0.0
``````
``````fit4<-fitMk(primate.tree,X4,model="ARD",pi="fitzjohn",
logscale=TRUE,lik.func="pruning")
fit4
``````
``````## Object of class "fitMk".
##
## Fitted (or set) value of Q:
##            Cathemeral   Diurnal Nocturnal
## Cathemeral  -0.031320  0.031320  0.000000
## Diurnal      0.000000 -0.002018  0.002018
## Nocturnal    0.002865  0.003715 -0.006581
##
## Fitted (or set) value of pi:
## Cathemeral    Diurnal  Nocturnal
##   0.001659   0.023625   0.974716
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -29.781298
##
## Optimization method used was "nlminb"
``````
``````ancr(fit4)
``````
``````## Marginal ancestral state estimates:
##    Cathemeral  Diurnal Nocturnal
## 91   0.000003 0.000587  0.999410
## 92   0.004205 0.017513  0.978282
## 93   0.099599 0.838156  0.062245
## 94   0.015881 0.982774  0.001345
## 95   0.001830 0.998136  0.000033
## 96   0.000217 0.999777  0.000006
## ...
##
## Log-likelihood = -29.781298
``````

Now, the latter.

``````X5<-X3
X5[labs,]<-rep(1/ncol(X5),ncol(X5))
``````
``````##                             Cathemeral   Diurnal Nocturnal
## Allenopithecus_nigroviridis  0.1000000 0.9000000 0.0000000
## Alouatta_palliata            0.1000000 0.9000000 0.0000000
## Alouatta_seniculus           0.1000000 0.9000000 0.0000000
## Aotus_trivirgatus            0.1000000 0.0000000 0.9000000
## Arctocebus_aureus            0.1000000 0.0000000 0.9000000
## Arctocebus_calabarensis      0.1000000 0.0000000 0.9000000
## Ateles_fusciceps             0.1000000 0.9000000 0.0000000
## Ateles_geoffroyi             0.1000000 0.9000000 0.0000000
## Ateles_paniscus              0.1000000 0.9000000 0.0000000
## Avahi_laniger                0.1000000 0.0000000 0.9000000
## Cacajao_melanocephalus       0.1000000 0.9000000 0.0000000
## Callicebus_moloch            0.1000000 0.9000000 0.0000000
## Callicebus_torquatus         0.1000000 0.9000000 0.0000000
## Callimico_goeldii            0.3333333 0.3333333 0.3333333
## Callithrix_argentata         0.1000000 0.9000000 0.0000000
## Callithrix_jacchus           0.1000000 0.9000000 0.0000000
## Callithrix_pygmaea           0.1000000 0.9000000 0.0000000
## Cebus_albifrons              0.1000000 0.9000000 0.0000000
## Cebus_capucinus              0.1000000 0.9000000 0.0000000
## Cercopithecus_cephus         0.1000000 0.9000000 0.0000000
``````
``````fit5<-fitMk(primate.tree,X5,model="ARD",pi="fitzjohn",
logscale=TRUE,lik.func="pruning")
fit5
``````
``````## Object of class "fitMk".
##
## Fitted (or set) value of Q:
##            Cathemeral   Diurnal Nocturnal
## Cathemeral  -0.031320  0.031320  0.000000
## Diurnal      0.000000 -0.002018  0.002018
## Nocturnal    0.002865  0.003715 -0.006581
##
## Fitted (or set) value of pi:
## Cathemeral    Diurnal  Nocturnal
##   0.001659   0.023625   0.974716
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -40.767421
##
## Optimization method used was "nlminb"
``````
``````ancr(fit5)
``````
``````## Marginal ancestral state estimates:
##    Cathemeral  Diurnal Nocturnal
## 91   0.000003 0.000587  0.999410
## 92   0.004205 0.017513  0.978282
## 93   0.099598 0.838157  0.062245
## 94   0.015881 0.982775  0.001345
## 95   0.001830 0.998136  0.000033
## 96   0.000217 0.999777  0.000006
## ...
##
## Log-likelihood = -40.767421
``````

Once again, we see a large effect on the likelihood – but none whatsoever on the parameter estimates or ancestral states. Has your mind been blown yet?

In the end, I’ll let the reader judge for her or himself; however, it makes sense to me to treat uncertainty in tip states in the second of our two approaches, in which each terminus is assigned a (prior) probability of being in each of the tip states (i.e., like sequencing error). This in turn makes the total probability of our data equal to the probability of each possible unobserved state multiplied by the (prior) probability that it is (in fact) in each of the states. Doing so also has the very convenient property of allowing us to include unequal ambiguity between each of the possible states for a trait at all tips.

That’s it, folks!