Hello Brian,
Could you please give a short example of how to define a transition matrix
with unobserved states? Assigning the states to the names of the rows and
columns doesn't seem to do it. Here's my attempt.
```
> library(corHMM)
>
> set.seed(802578)
>
> tr <- rtree(10)
>
> dat <- data.frame(cbind(tr$tip.label, c(0, 0, 1, 3, 3, 0, 4, 4, 0, 3))) #
No state 2
> # Recycling a bit of David's code. As far as I can see, the results are
equivalent to what rate.mat.maker can do.
> model<-matrix(NA, 5, 5)
> for(i in 1:5){for(j in 1:5){
+ if(abs(i-j)<2){
+ model[i,j]<-1
+ }
+ }}
>
> rownames(model)<-colnames(model)<-0:4
> diag(model) <- NA
> print(model)
0 1 2 3 4
0 NA 1 NA NA NA
1 1 NA 1 NA NA
2 NA 1 NA 1 NA
3 NA NA 1 NA 1
4 NA NA NA 1 NA
>
> res <- rayDISC(tr, dat, ntraits = 1, charnum = 1, rate.mat = model)
No model selected for 'node.states'. Will perform marginal ancestral state
estimation.
State distribution in data:
States: 0 1 3 4
Counts: 4 1 3 2
Initializing...
Finished. Beginning thorough search...
Finished. Inferring ancestral states using marginal reconstruction.
Error in `rownames<-`(`*tmp*`, value = c("0", "1", "3", "4")) :
length of 'dimnames' [1] not equal to array extent
In addition: There were 42 warnings (use warnings() to see them)
```
I feel I'm missing something really simple.
Cheers,
Eduardo
2015-06-17 21:25 GMT+02:00 Brian O'Meara <[email protected]>:
> 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 <[email protected]> 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 <[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]/
> >
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-phylo mailing list - [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at http://www.mail-archive.com/r-
> [email protected]/
>
[[alternative HTML version deleted]]
_______________________________________________
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]/
