Wright’s population genetic threshold model for discrete traits is an attractive trait model that may be biologically more realistic than alternatives like the Mk model. It assumes that the discrete traits we observe are the expression of an unobserved continuous trait (liability). Depending on what side of a threshold the liability value for an individual falls determines what state the individual expresses. Felsenstein has extended this model into phylogenetics where we imagine that these liabilities evolve via a Brownian motion model along a phylogeny (paper). Liam Revell has developed a nice implementations that allows us to estimate these underlying liability scores and analyze discrete traits in the same ways we would analyze continuous traits. One current limitation of these models though is that they are limited to binary traits (although multi-state ordered traits would be a simple extension from current methods).
For a project I am working on I wanted to extend this to multi-state unordered traits. As a start I thought that I might first tackle simulating under this model for a simple 3-state unordered character. In this case we can imagine that the liability score is represented by a value evolving in two dimensions depicted here:
Now we can think about this process playing out with a set of thresholds that divide this 2-dimensional space into 3 areas. Like this:
Next we simply move this process onto a phylogeny where this evolution is occurring along the branches of phylogeny. Here is the code that I use to do this:
SimThresh3 <- function(tree,
liabilities = F,
sig2 = 1,
root = c(0, 0)){
returnState <- function(x,y){
state <- vector()
# angle A is at origin
# angle B is at 1,0
# angle C is at x,y
# so we can calculate the sides like this
a <- sqrt((x-1)^2 + (y)^2)
c <- sqrt(x^2+y^2)
b <- sqrt(1)
# and then we can calculate the angle at the origin like this:
A <- (acos((b^2 + c^2 - a^2)/(2*b*c))) * (180/pi)
state <- vector(mode="numeric", length=length(y))
state[(y >= 0 & A < 90) | (y < 0 & A < 30)] <- 1
state[(y >= 0 & A > 90) | (y < 0 & A > 150)] <- 2
state[y < 0 & (A > 30 & A < 150)] <- 3
names(state) <- names(x)
return(state)
}
x<-fastBM(tree, sig2 = sig2, a = root[1])
y<-fastBM(tree, sig2 = sig2, a = root[2])
tip.state <- returnState(x,y)
if(liabilities == F) return(tip.state)
if(liabilities == T){
results <- vector("list", 3)
results[[1]] <- tip.state
results[[2]] <- x
results[[3]] <- y
names(results) <- c("observed", "liab1", "liab2")
return(results)
}
}
We can now run this code on a simulated tree and see what it looks like:
library(phytools)
library(evobiR)
tree<-pbtree n="500," scale="1)</font">
tip.state.full <- liabilities="T)</font" simthresh3="" tree="">
ShowTree(tree, tip.state.full[[1]],
cols = c("red", "blue", "green"),
tip.cex=.5, type="fan")
Alternatively we could visualize where our tips ended up in this two dimensional space:
One of the particularly attractive characteristics of the threshold model is that it implicitly allows for variation in the rate of transitions across the tree. For instance if a lineage is very close to a threshold value then transitions between states should be common. In contrast if an individual is very far from any threshold then we would predict that transitions should be very rare. Below points A and B respectively illustrate these two possibilities.
We can do a spot check of this just by looking at our tree and looking up the liability values of two tips that seem to illustrate these extremes. Here I have zoomed in on the lower right quadrant of our tree and we can see two region that we would predict would have liabilities roughly in the areas indicated above. B being largely fixed for blue far from the thresholds with red and green, and A somewhere close the the threshold between red and green.
When we then plot these two points on tip of the whole distribution we see just what we expected.
This functions (SimThresh3, and ShowTree) are both in my R package evobiR but I haven't done any further testing than what you see here so if you use them make sure they are really doing what you want them to as they will likely change significantly as I have a chance to work on them.
Contents of this website are good and appreciative. Congratulation
ReplyDelete3 Track Sliding Windows Bhilai