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]/
