Hi Jingchun, Yes in this case I'd use the "true" model. The LRT may be not very powerful with your data.
Cheers, Emmanuel -----Original Message----- From: Jingchun Li <jingc...@umich.edu> Date: Fri, 5 Apr 2013 10:25:57 To: <emmanuel.para...@ird.fr> Cc: <r-sig-phylo@r-project.org> Subject: Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares Thank you very much Daniel and Emmanuel! Now things start to make sense. I guess my last question is, if there is biological evidence that the transition rates between the two traits are different, is it even meaningful to compare the two models? Should I simply avoid using the equal rates model? Cheers, Jingchun On Fri, Apr 5, 2013 at 9:12 AM, Emmanuel Paradis <emmanuel.para...@ird.fr>wrote: > ** > Jingchun, > > Thanks to Daniel for giving a detailed answer to your first question. > > You can use ace() to fit a model with irreversible change. You need to > build a custom model like this: > > > Q <- matrix(0, 2, 2) > > Q[2, 1] <- 1 > > Q > [,1] [,2] > [1,] 0 0 > [2,] 1 0 > > In this case there's one rate to be estimated and the other is fixed to > zero (the 0's on the diagonal are ignored). This can be generalized to more > complicated models (there's an example in my book). Then call ace with: > > ace(traits_data, tree, type = "d", model = Q) > > This time you cannot compare this model with the "ER" one with an LRT: > compare directly the log-liks since they have the same number of parameters > or use AIC(). BTW, there is an anova methods to compare nested models > fitted with ace: this avoids the pain to compute correctly df's. > > Cheers, > > Emmanuel > ------------------------------ > *From: * Jingchun Li <jingc...@umich.edu> > *Date: *Thu, 4 Apr 2013 22:25:09 -0400 > *To: *<emmanuel.para...@ird.fr>; <r-sig-phylo@r-project.org> > *Subject: *Re: [R-sig-phylo] ML ancestral state reconstruction using > different softwares > > Sorry, I forgot to include the list. there we go. > >> ------------------------------ >> *From: * Jingchun Li <jingc...@umich.edu> >> *Date: *Thu, 4 Apr 2013 21:58:56 -0400 >> *To: *<emmanuel.para...@ird.fr> >> *Subject: *Re: [R-sig-phylo] ML ancestral state reconstruction using >> different softwares >> >> Thank you very much, Dr.Paradis. You are right that the two parameter >> model does not fit the data significantly better. I guess I'm more puzzled >> about why the three softwares gave me different likelihoods. >> >> I did fit the symmetrical model using BayesTraits later and it gave me >> about the same transition rates as the asymmetrical model and a >> similar likelihood value (-14.08 for symmetrical and -14.00 for >> asymmetrical). However I do have biological evidence that transition from 0 >> to 1 might be extremely hard. So maybe I can also try set q01 to 0 (which >> is actually what ace gave me under the asymmetrical model)? >> >> Cheers, >> >> Jingchun >> >> >> On Thu, Apr 4, 2013 at 9:21 PM, Emmanuel Paradis <emmanuel.para...@ird.fr >> > wrote: >> >>> Hi, >>> >>> You cannot really say that "in both ace and Mesquite, it slightly favors >>> the asymmetrical model": the increase in log-likelihood is less than 1 >>> which is far from being significant with one additional parameter. So it >>> seems that all three pieces of software agree well on the estimate of rate >>> for the symmetrical model. You should fit also this model with BayeTraits >>> to do the full comparison. >>> >>> Cheers, >>> >>> Emmanuel >>> -----Original Message----- >>> From: Jingchun Li <jingc...@umich.edu> >>> Sender: r-sig-phylo-boun...@r-project.org >>> Date: Thu, 4 Apr 2013 10:15:05 >>> To: <r-sig-phylo@r-project.org> >>> Subject: [R-sig-phylo] ML ancestral state reconstruction using different >>> softwares >>> >>> Hi all, >>> >>> I am exploring different methods for ancestral state reconstruction for a >>> small phylogeny I have (26 OTUs). I have one binary trait coded as 0s and >>> 1s for all the taxa. And I am more interested in the transition rates >>> between the two states. >>> >>> If I understand correctly, a maximum likelihood ancestral state >>> reconstruction using a MK1 model or an asymmetrical rates model should >>> give >>> me the optimized estimated transition rates between the two states. The >>> formal will assume the two rates are the same, the later will >>> assume independent rates. >>> >>> So I firstly did this using the ape function in ace. >>> >>> ERreconstruction <- ace(traits_data, tree, type="discrete", model="ER") >>> ERreconstruction$rates >>> ERreconstruction$loglik >>> >>> ARDreconstruction <- ace(traits_data, tree, type="discrete", model="ARD") >>> ARDreconstruction$rates >>> ARDreconstruction$loglik >>> >>> Then I tried using Mesquite, using MK1 model and the "Asymmetrical >>> 2-param. >>> Markov-k Model " >>> >>> Then I tried using BayesTraits, with the MultiState -> Maximum Likelihood >>> option, with no restrictions on q10 and q01. >>> >>> What puzzled me is that I'm getting different results from the three >>> softwares. >>> -------------------------- >>> aceER: >>> transition rate: 0.0059, likelihood: -7.00 >>> aceARD: >>> transition rate: 0.019 and 0.00, likelihood: -6.11 >>> >>> Mesquite MK1: >>> transition rate: 0.0059, likelihood: -7.70 >>> Mesquite Asymmetrical: >>> transition rate: 0.015 and 0.005, likelihood: -7.30 >>> >>> BayesTraits no restriction: >>> transition rate: 0.0057 and 0.0057, likelihood: -14.09 >>> -------------------------------- >>> >>> I can see that in both ace and Mesquite, >>> it slightly favors the asymmetrical model with two different rates, and >>> the >>> rates and likelihoods are more of less comparable. But BayesTraits seems >>> to >>> think the two rates should be equal, and it has a different likelihood. >>> >>> Does this has something to do with different searching algorithms? Or am >>> I >>> missing something here? I am aware that the scaled likelihoods in ace are >>> the scaled conditional likelihoods, not the joint or marginal >>> reconstructions. But this should not affect the estimated transition >>> rates? >>> >>> >>> Thank you very much! >>> >>> [[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/ >>> >> >> > [[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/
