corHMM should be able to do this as long as the states are defined in the transition matrix (rate.mat.maker() to start, then modify). BioGeoBEARS could do this (since it has to for biogeography: not all combos of ancestral areas are observed at the tips), and one of its models might be equivalent to what you want. I imagine diversitree can, too.
Of course, this all might work best with a model with somewhat constrained rates. If rates are all free, and the model is unordered, the MLE for rate of moving into the unseen state will be zero, and rate for moving out of it will be Inf. Best, Brian _______________________________________ Brian O'Meara Assistant Professor Dept. of Ecology & Evolutionary Biology U. of Tennessee, Knoxville http://www.brianomeara.info Postdoc collaborators wanted: http://nimbios.org/postdocs/ Calendar: http://www.brianomeara.info/calendars/omeara On Wed, Jun 17, 2015 at 3:06 PM, David Bapst <dwba...@gmail.com> wrote: > 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 <liam.rev...@umb.edu> > 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: liam.rev...@umb.edu > > 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 - R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > Searchable archive at > http://www.mail-archive.com/r-sig-phylo@r-project.org/ > [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
