Hi David & Liam.You can fit a model with unobserved states with ace() because it considers the levels of 'x', not the observed states. You may set the levels of x with the function of the same name:
> x <- factor(sample(2, size = 30, replace = TRUE)) > levels(x) <- LETTERS[1:3] > table(x) x A B C 12 18 0You can then fit different models, for instance one model where the (direct) transition A <-> B is not possible would be specified with the following matrix:
> m A B C A 0 0 1 B 0 0 1 C 1 1 0 Best, Emmanuel Le 17/06/2015 23:32, Liam J. Revell a écrit :
Hi David. The diagonal of Q is defined such that row sum (or sometimes the column sum) is zero. This is the case for any transition rate matrix of a continuous time Markov chain & is required so that matrix exponentiation of Q*t give the matrix P[i,j] in which i,j gives the probability of going from i->j. In ace, rerootingMethod, etc. it is possible to specify the parameterization of the model (i.e., how many different transition rates to estimate, which ones to set equal to one another, and which ones to set to zero) - but not to constrain any particular rate to a specific value (other than zero). With make.simmap it is possible to fix Q at some arbitrary (valid) value. It will even give you the likelihood for that value, if I remember correctly. 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 4:29 PM, David Bapst wrote:Liam- Yes, that works! I apparently confused the sort of model matrix needed by rTraitDisc with the sort needed by ace. Here's a full working example, for future posterity: ################# library(ape) library(phytools) data(bird.orders) #make a matrix for a 4 step trait model<-matrix(c(0,1,0,0, 1,0,1,0, 0,1,0,1, 0,0,1,0),4,4) #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) ################# An interesting thing to note about the above example is that even with equal rates between sequential states, it infers a higher negative rate of staying in the missing states, so a better model might be one in we constrained the rate of staying in a state as the same across all states. However, that doesn't seem to be possible with ace/rerootingMethod as the diagonal of the model matrix is ignored, right? Brian, yes, corHMM seems like it should be able to do it; I overlooked that package. diversitree's ML reconstruction functions don't seem to accept a model matrix, at least based on current documentation. -Dave On Wed, Jun 17, 2015 at 1:30 PM, Liam J. Revell <[email protected]> wrote:Hi David. I think the problem may be that you're specifying the model matrix incorrectly. The model matrix is an integer matrix in which each number is a different rate type and the diagonal is zeroes. So, for instance, for a single-rate ordered model with four states you would do: model<-matrix(c(0,1,0,0, 1,0,1,0, 0,1,0,1, 0,0,1,0),4,4) (with row & column names of your states). A symmetric ordered model would be: model<-matrix(c(0,1,0,0, 1,0,2,0, 0,2,0,3, 0,0,3,0),4,4) ## I think Does this help? In my tests it seems to run fine for phytools. 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 3:06 PM, David Bapst 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_______________________________________________ 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]/
_______________________________________________ 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]/
