tag:blogger.com,1999:blog-8499895524521663926.post6794475090731947905..comments2019-07-21T05:31:35.548-04:00Comments on Phylogenetic Tools for Comparative Biology: Bayesian MCMC for the threshold modelLiam Revellhttp://www.blogger.com/profile/04314686830842384151noreply@blogger.comBlogger17125tag:blogger.com,1999:blog-8499895524521663926.post-87071290873970313502016-03-09T13:08:30.308-05:002016-03-09T13:08:30.308-05:00Basically, no. One thing that some collaborators o...Basically, no. One thing that some collaborators of mine elected to do was use ancThresh to sample liabilities under a multi-state univariate model, and then computed the correlation of those liabilities (for tips). This is a bit ad hoc, but may be reasonable to try....Liam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-12782661061793878092016-03-09T12:49:54.831-05:002016-03-09T12:49:54.831-05:00Hi Liam,
Has a >2-state discrete model been de...Hi Liam,<br /><br />Has a >2-state discrete model been developed? That is, will threshBayes handle >2 discrete states?<br /><br />Thanks!<br /><br />Frank Frank Burbrinkhttps://www.blogger.com/profile/08701731074140226071noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-11488088456595269592014-10-08T17:22:17.629-04:002014-10-08T17:22:17.629-04:00Hi David,
I'm trying to figure out how to use ...Hi David,<br />I'm trying to figure out how to use threshold models, so I just saw your question about ternary diagrams. Yes, they are used often in biomechanics, especially to examine the relative lengths of limb elements, and how this relates to different types of locomotion. See: Middleton and Gatesy 1998 and Ibrahim et al 2014. Just thought I'd post in case you were still curious, 2 years down the road.Talia Moorehttps://www.blogger.com/profile/02970041413911277466noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-27289030910998941282012-08-19T01:14:31.442-04:002012-08-19T01:14:31.442-04:00If Hadfield and Nakagawa use a probit function (a ...If Hadfield and Nakagawa use a probit function (a normal distribution function) to achieve a thrshold model in MCMCglimm, they are restricting their model too much. Then the environmental effect in each character would be independent. Better to use a step-function instead, and just have possibly-covarying environmental effects on the liability scale.<br /><br />The issue of how to go beyond a two-state model, or a model with three or more thresholds all on a single liability scale, is puzzling and I think we need more work on it. I don;t think that the multiple possible models are a recommendation to use an Mk model.Joe Felsensteinhttps://www.blogger.com/profile/06359126552631140000noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-92196760453741381392012-08-18T02:36:53.455-04:002012-08-18T02:36:53.455-04:00Oh wait, that would make state A3 more probable th...Oh wait, that would make state A3 more probable than states A2 and A1. Hrm, okay, but it definitely makes sense for pure unordered multistate!David Bapsthttps://www.blogger.com/profile/17606476387441191531noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-67916666977429194212012-08-18T02:34:39.465-04:002012-08-18T02:34:39.465-04:00Yeah, I don't know about the identifiability. ...Yeah, I don't know about the identifiability. I was just thinking out loud about how one could describe an unordered multistate under a threshold model.<br /><br />My ternary diagram allusion is incomplete, so let me abandon that. (It still applies, but its because the math on a ternary diagram is weird.) <br /><br />Instead, imagine a lineage at some time as having n continuously-evolving traits, where n is the number of unordered states. Each continuous trait corresponds to a different state, such that whichever trait is greatest is the expressed state.<br /><br />One could see how this could even account for further complications common to morphological data, like partially-ordered multistate characters... For example, a trait where there are three unordered states, but one state (state A) has three ordered sub-states (states A1, A2, A3). Then the only question is if state A is being expressed (is trait value A' greater than trait values B' and C') and is the value of trait A' within the proper thresholds for A1, A2, or A3? This would be a very simple model to (at least) simulate a complex morphological trait under.<br /><br />Hah, okay, this line of thinking turned out to be pretty tangential in retrospect.David Bapsthttps://www.blogger.com/profile/17606476387441191531noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-70323613764105606262012-08-17T21:37:30.507-04:002012-08-17T21:37:30.507-04:00Hi Vladimir.
I'm not sure that there isn'...Hi Vladimir.<br /><br />I'm not sure that there isn't information about ordering in phylogenetic data. . . . To be honest, however, I'm only starting to think about this problem - so, yes, I am worried, but it is an intriguing idea nonetheless.<br /><br />- LiamLiam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-22974116820443946162012-08-17T21:32:20.572-04:002012-08-17T21:32:20.572-04:00Hi David.
Hmmm, maybe. I'm not sure whether t...Hi David.<br /><br />Hmmm, maybe. I'm not sure whether this would fix the identifiability program of Vladimir. <br /><br />In addition, does putting the three-state trait on two axes create an implied ordering? (I.e., state A is below the threshold on axis 1; state B is above the threshold on axis 2; state C is past the threshold on axes 1 & 2? In this case, you can't get to state C from stat without going through state B. Am I thinking about your two axis model incorrectly?) <br /><br />- LiamLiam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-51923312775351315892012-08-17T20:58:04.472-04:002012-08-17T20:58:04.472-04:00Couldn't one do unordered multistate in a thre...Couldn't one do unordered multistate in a threshold model by just extending the 'threshold' value into multiple dimensions? <br /><br />I've got this mental image that an unordered three-state threshold model would kinda be like a little point moving around on a ternary diagram, which essentially a two-dimensional space... (Do ternary diagrams show up in biology at all? We use them in geology all the time.)David Bapsthttps://www.blogger.com/profile/17606476387441191531noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-18360972957556040782012-08-17T18:27:38.748-04:002012-08-17T18:27:38.748-04:00Liam, are worried about identifiability in the cas...Liam, are worried about identifiability in the case of unknown order? Usually, the more latent variables you introduce, the less identifiable the model becomes. If the order of observed categories is unknown, I don't see how it can be recovered...Vladimir Mininhttps://www.blogger.com/profile/00603257472903019482noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-24145860451698342182012-08-17T18:20:49.495-04:002012-08-17T18:20:49.495-04:00OK - I have no posted an example simulation & ...OK - I have no <a href="http://phytools.blogspot.com/2012/08/bayesian-mcmc-threshold-model.html" rel="nofollow">posted an example simulation & analysis</a>, as promised.Liam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-12487517397005682502012-08-17T17:44:47.390-04:002012-08-17T17:44:47.390-04:00@Matt - Very cool. You should try this with data a...@Matt - Very cool. You should try this with data actually simulated under the threshold model. I will post code to do this in a few minutes (but I'm sure you could easily figure out how to do this without my help).<br /><br />@Vladimir - Most likely you are right, although I'm forced to admit I don't totally understand it. You're right that the threshold model doesn't make sense for a multistate character with unordered states. What I have in mind is the situation where there is an ordering, but it is not known a priori. In this case, I imagine we could use MCMC to sample the orders as well (in addition to the relative distances between thresholds which matter when we have more than three states). (We might use rjMCMC to jump between ordering, but I'm not sure we'd need to.)Liam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-46209539211456827022012-08-17T17:33:27.711-04:002012-08-17T17:33:27.711-04:00Hi Liam,
I think all Hadfield and Nakagawa are tr...Hi Liam,<br /><br />I think all Hadfield and Nakagawa are trying to say is that once you convert a pedigree, in quantitative genetic models, or a phylogeny, in phylogenetic threshold models, into a covariance matrix that depends on hyperparameters (variance components in quantitative genetics and Brownian motion variance in the threshold model), the mathematical formulation of these two types of models is identical - just regular latent Gaussian models (aka probit regression models). So it is mathematical (algebraic) equivalence, with need to test it statistically, I think.<br /><br />It is my understanding that MCMCglmm already can handle arbitrary ordinal variables, not only binary. I don't immediately see how one could apply the threshold model to categorical (unordered) data.<br /><br />Vladimir<br /><br />P.S. Enjoying your blog and code! Keep up the good work.Vladimir Mininhttps://www.blogger.com/profile/00603257472903019482noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-4985478561761356562012-08-17T17:30:26.732-04:002012-08-17T17:30:26.732-04:00And even if this is equivalent (which again I am n...And even if this is equivalent (which again I am not convinced of), I don't believe it can be naturally extended to ordered states.mwpennellhttps://www.blogger.com/profile/02523892313236801738noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-85728608265932690502012-08-17T17:25:21.430-04:002012-08-17T17:25:21.430-04:00Vladimir and Liam,
After thinking about it, this ...Vladimir and Liam,<br /><br />After thinking about it, this may similar to fitting a phylogenetic glmm model with a logit link function though I am completely unsure about this. This code may be completely wrong but it runs.<br /><br />require(MCMCglmm)<br />require(geiger)<br /><br /><br /># i am just going input a simulated tree of 100 taxa<br /><br />tree <- drop.tip(birthdeath.tree(1, 0, taxa.stop=101), "101")<br /><br />## simulate two completely uncorrelated traits just for an example<br /><br /><br />q <- list(rbind(c(-.5, .5), c(.5, -.5)))<br /><br />trait1 <- sim.char(tree, q, model="discrete", n=1)<br /><br />trait2 <- sim.char(tree, q, model="discrete", n=1)<br /><br />data <- cbind(trait1, trait2)<br /><br />rownames(data) <- tree$tip.label<br /><br /><br /># The variable animal holds the tip labels of the tree<br /><br />animal <- tree$tip.label<br /><br /><br /># Create a data frame with all the variables (after converting them to factors if need be). <br /># One column needs to be the tip labels (this seems redundant but it is how i got it to work)<br /><br />traits <- as.data.frame(cbind(animal, data))<br />rownames(traits) <- animal<br /><br /><br /># Create a fixed effects model (the distinction between fixed and random effect in the context of comparative methods and in a Bayesian MCMC framework is not that clear to me)<br /><br />fixed <- trait1 ~ trait2<br /><br /><br /># Create a random effects model. This should at least include your phylogeny ("animal")<br />random <- ~ animal<br /><br /><br /># Since the dependent variable is categorical (same as binary for the purposes of this package)<br />family <- "categorical"<br /><br /><br /># Take the inverse of the phylogenetic variance-covariance matrix. This is actually technically not the inverse of the phy vcv matrix. It is actually a sparse matrix, as it includes empty values for all ancestoral nodes. For reasons detailed in Hadfield and Nakagawa 2010, this drastically increases the speed of the calculations (dont have to actually do the inversion) and memory requirements.<br />Ainv <- inverseA(pedigree=tree, nodes="ALL")$Ainv<br /><br /><br /># Assign priors for the R-matrix and G-matrix<br /># I am just going to use the very diffuse prior that is the default.<br /># Have not investigated the effect of the priors but this needs to be done.<br />prior = "NULL"<br /><br /><br /><br /># Run MCMC chain #<br /><br /># Instead of inserting the phylogeny as the argument "pedigree" I find it useful just to put the inverted vcv matrix in the argument ginverse, which i believe is equivalent<br /><br />res <- MCMCglmm(fixed=fixed, random=random, family=family, data=traits, ginverse=list(animal=Ainv), nitt=10000, burnin=1000)<br /><br /><br /><br /># look at the summary<br />summary(res)<br /><br /><br /><br />mwpennellhttps://www.blogger.com/profile/02523892313236801738noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-59108134400226646012012-08-17T16:12:39.074-04:002012-08-17T16:12:39.074-04:00Hi Vladimir.
Thanks for the comment.
This is cer...Hi Vladimir.<br /><br />Thanks for the comment.<br /><br />This is certainly possible, of course. The issue is specifically addressed (but without a satisfying resolution) in Felsenstein (2012; p. 146) in the passage beginning "Hadfield and Nakagawa (2010) have noted that all such models are equivalent to multivariate “mixed models” of quantitative genetics. For discrete traits. . . ."<br /><br />Do you know how you would design this test using the MCMCglmm package? Is this documented in the literature or on the web? A simple Google search of "threshold" with "MCMCglmm" doesn't turn up any useful results. How about natural extensions of this model to multiple ordered states and multiple states with an underlying (but unknown) order?<br /><br />Thanks for the suggestion! LiamLiam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-1831557466713760892012-08-17T15:55:45.657-04:002012-08-17T15:55:45.657-04:00Since you didn't mention this in the post, I w...Since you didn't mention this in the post, I want to point out that MCMCglmm package already has a Bayesian version of the threshold model implemented, I think.Vladimir Mininhttps://www.blogger.com/profile/00603257472903019482noreply@blogger.com