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))
head(X,20)
## 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))
head(X1,20)
## 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
## Piliocolobus_badius 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)
add.arrow(primate.tree,tips,arrl=3,offset=2,lwd=2,
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))
head(X2,20)
## 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
## Piliocolobus_badius 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)
add.arrow(primate.tree,tips,arrl=3,offset=2,lwd=2,
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)
head(X3,40)
## 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
## Daubentonia_madagascariensis 0.1 0.0 0.9
## 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
## Piliocolobus_badius 7.6e-05 0.999924 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))
head(X4,20)
## 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))
head(X5,20)
## 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!
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.