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 <liam.rev...@umb.edu> 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: liam.rev...@umb.edu
> 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 <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/
