Liam,
Ah, I see, rerootingMethod(), unlike ace, will accept a matrix
representation of discrete trait values.
However, it doesn't appear that one can define a model matrix still?
Outside of molecular data, the missing state issue is often for
ordered data, hence the need to define a model matrix.
################################
library(ape)
library(phytools)
data(bird.orders)
#make a matrix for a 4 step trait
model<-matrix(0,4,4)
for(i in 1:4){for(j in 1:4){
if(abs(i-j)<2){
model[i,j]<-0.1
}
}}
rownames(model)<-colnames(model)<-0:3
#simulate data where you have only observed endmembers
trait <- factor(c(rep(0, 5), rep(3, 18)))
names(trait)<-bird.orders$tip.label
trait<-to.matrix(trait,seq=0:3)
rerootingMethod(tree=bird.orders, x=trait)
#works
rerootingMethod(tree=bird.orders, x=trait, model=model)
#fails
#################
Cheers,
-Dave
On Wed, Jun 17, 2015 at 12:06 PM, Liam J. Revell <[email protected]> wrote:
> Hi David.
>
> Actually, if I understand correctly, this does work in phytools. I suspect
> it can be set up in phangorn as well - Klaus can probably explain how.
>
> So, if I understand correctly, you want to get marginal ancestral states for
> a character that can assume, say, states "A", "C", "G", & "T", but for which
> you have only observed states "A", "C", & "G" (for instance) in your tip
> taxa?
>
> The way to do this in rerootingMethod is as follows:
>
> library(phytools)
> ## simulate some data to work with:
> tree<-pbtree(n=26,tip.label=letters,scale=1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> rownames(Q)<-colnames(Q)<-c("A","C","G")
> x<-sim.history(tree,Q)$states
> x ## our data
> ## now convert to a matrix representation:
> X<-to.matrix(x,seq=c("A","C","G","T"))
> ## reconstruct ancestral states
> fitER<-rerootingMethod(tree,X,model="ER")
> fitER
>
> Note that as we have not observed any evidence that changes to "T" occur for
> this character, if we use a symmetric model instead, we will find that the
> marginal likelihoods at all internal nodes for state "T" are zero. This is
> because the maximum likelihood q[i,"T"] & q["T",i] for all i will be zero
> (sensibly). For instance:
>
> fitSYM<-rerootingMethod(tree,X,model="SYM")
> fitSYM
>
> Let me know if this is what you had in mind.
>
> rerootingMethod just uses modified code from ace internally to do the
> calculations - so if ace does not permit this presently, it could easily be
> modified to allow it. phangorn too must allow this because for many
> nucleotide sites we want to reconstruct ancestral states - but we will have
> observed only a subset of the four nucleotides at that site.
>
> All the best, Liam
>
> Liam J. Revell, Assistant Professor of Biology
> University of Massachusetts Boston
> web: http://faculty.umb.edu/liam.revell/
> email: [email protected]
> blog: http://blog.phytools.org
>
>
> On 6/17/2015 1:46 PM, David Bapst wrote:
>>
>> Hello all,
>>
>> (I'm a troublemaker today.)
>>
>> Sometimes, in ordered discrete data, there are states we know might
>> exist as intermediary between observed states but aren't observed
>> themselves. I suppose this is probably common for meristic data. At
>> least to me, it seems like it should be possible to reconstruct node
>> states in a scenario with 'missing' states in a ML approach by
>> defining the corre
>>
>> Anyway, in a recent conversation with a friend, he noted that he
>> couldn't get the function for ML ancestral reconstruction of discrete
>> characters, specifically ace (in ape) or rerootingMethod (in phytools)
>> to accept data with missing states. I was a little shocked by this,
>> and I didn't believe him in my foolhardiness and investigated myself.
>>
>> #################################
>>
>> library(ape)
>> data(bird.orders)
>>
>> #make an ordered matrix for a 4 state trait
>> model<-matrix(0,4,4)
>> for(i in 1:4){for(j in 1:4){
>> if(abs(i-j)<2){
>> model[i,j]<-0.1
>> }
>> }}
>> rownames(model)<-colnames(model)<-0:3
>>
>> #simulate data where you have only observed endmembers
>> trait <- factor(c(rep(0, 5), rep(3, 18)))
>> names(trait)<-bird.orders$tip.label
>>
>> ace(trait, bird.orders, type = "discrete")
>> #works
>> ace(trait, bird.orders, type = "discrete", model=model)
>> #fail
>>
>> library(phytools)
>> rerootingMethod(tree=bird.orders, x=trait)
>> #works
>> rerootingMethod(tree=bird.orders, x=trait, model=model)
>> #fail
>>
>> #################################
>>
>> (And it doesn't look like ancestral.pml in phangorn accepts model
>> matrices.)
>>
>> So, are these functions rejecting a model with missing states due to a
>> real analytical issue that I'm blithely ignorant of, or is there some
>> other function that I've missed that will allow ML reconstruction with
>> missing states?
>>
>> Cheers,
>> -Dave
>>
>
--
David W. Bapst, PhD
Adjunct Asst. Professor, Geology and Geol. Eng.
South Dakota School of Mines and Technology
501 E. St. Joseph
Rapid City, SD 57701
http://webpages.sdsmt.edu/~dbapst/
http://cran.r-project.org/web/packages/paleotree/index.html
_______________________________________________
R-sig-phylo mailing list - [email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/[email protected]/
