Additional method for fitPagel to fit Pagel's (1994) model of correlated binary character evolution

I just updated fitPagel, a phytools function that fit the Pagel (1994) method for testing whether theh evolutiokn of one binary character affects a second character (or vice versa).

The way this method works is it simply re-codes the two binary characters into a single four state character - e.g., 0|0, 0|1, etc. Then it fits two different models. In the first model, the 'independent' model, the rate of transition in the first character is constrained to be equal for a given transition type, irregardless of the state for the second trait. So, for instance, the rate of transition 0|0 -> 1|0 is the same as 0|1 -> 1|1 because both transitions involve the same change (0 -> 1) in character 1. The second model, the 'dependent' model, allows each type of transition to have a different rate. In other words 0|0 -> 1|0 and 0|1 -> 1|1 can have different rates (different probabilities of occuring on a given time interval), even though the exact same type of change is occuring for the same character in both instances. Obviously, transitions 0|0 -> 1|1 are constrained to have a rate of zero in both models.

The update now permits the phytools function fitMk to be used internally to fit the Mk model. Previously the two methods available were method = "ace" and method = "fitDiscrete", and these are both still available as options.

Let's compare them:

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)
library(phytools)

## Loading required package: ape
## Loading required package: maps

We can start by simulating data under the 'dependent' model

## simulate tree
tree<-pbtree(n=300,scale=1)
## simulate data for the two characters
Q<-matrix(c(0,0.4,0.4,0,2,0,0,2,2,0,0,2,0,0.4,0.4,0),4,4,byrow=TRUE)
rownames(Q)<-colnames(Q)<-c("aa","ab","ba","bb")
diag(Q)<--rowSums(Q)
tt<-sim.history(tree,Q)

## Note - the rate of substitution from i->j should be given by Q[j,i].
## Detecting that rows, not columns, of Q sum to zero :
##   Transposing Q for internal calculations.
## Done simulation(s).

## split the 'aa','ab','bb','ba' data into two characters
t1<-mergeMappedStates(tt,c("aa","ab"),"a")
t1<-mergeMappedStates(t1,c("ba","bb"),"b")
t2<-mergeMappedStates(tt,c("aa","ba"),"a")
t2<-mergeMappedStates(t2,c("ab","bb"),"b")
t1$states<-getStates(t1,"tips")
t2$states<-getStates(t2,"tips")
## visualize the correlated evolution of the two traits
par(mfrow=c(1,2))
plotSimmap(t1,setNames(c("red","blue"),letters[1:2]),lwd=1,ftype="off")
plotSimmap(t2,setNames(c("red","blue"),letters[1:2]),lwd=1,ftype="off",
    direction="leftwards")

## extract the tip data from our simulation
x<-getStates(t1,"tips")
y<-getStates(t2,"tips")

Now we're ready to fit our model using all three methods:

fitPagel(tree,x,y) ## default method is "fitMk"

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -0.8726028  0.6514011  0.2212017  0.0000000
## a|b  0.3220958 -0.5432975  0.0000000  0.2212017
## b|a  0.6189176  0.0000000 -1.2703187  0.6514011
## b|b  0.0000000  0.6189176  0.3220958 -0.9410134
## 
## Dependent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -0.5951317  0.3393103  0.2558215  0.0000000
## a|b  4.6486390 -4.7698708  0.0000000  0.1212318
## b|a  2.2086339  0.0000000 -4.5275279  2.3188940
## b|b  0.0000000  0.4113254  0.3105573 -0.7218827
## 
## Model fit:
##             log-likelihood
## independent      -172.8745
## dependent        -150.4001
## 
## Hypothesis test result:
##   likelihood-ratio:  44.94869 
##   p-value:  4.074835e-09 
## 
## Model fitting method used was fitMk

fitPagel(tree,x,y,method="ace")

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -0.8725883  0.6513844  0.2212039  0.0000000
## a|b  0.3220944 -0.5432984  0.0000000  0.2212039
## b|a  0.6189141  0.0000000 -1.2702985  0.6513844
## b|b  0.0000000  0.6189141  0.3220944 -0.9410086
## 
## Dependent model rate matrix:
##           a|a        a|b        b|a        b|b
## a|a -0.595134  0.3393112  0.2558228  0.0000000
## a|b  4.648654 -4.7698873  0.0000000  0.1212328
## b|a  2.208634  0.0000000 -4.5275253  2.3188912
## b|b  0.000000  0.4113245  0.3105577 -0.7218823
## 
## Model fit:
##             log-likelihood
## independent      -171.4882
## dependent        -149.0138
## 
## Hypothesis test result:
##   likelihood-ratio:  44.94869 
##   p-value:  4.074835e-09 
## 
## Model fitting method used was ace

fitPagel(tree,x,y,method="fitDiscrete")

## Loading required package: geiger
## Loading required package: parallel

## Warning in fitDiscrete(tree, xy, model = iQ): Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

## Warning in fitDiscrete(tree, xy, model = dQ): Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -0.8815200  0.6639046  0.2176153  0.0000000
## a|b  0.3054217 -0.5230370  0.0000000  0.2176153
## b|a  0.6199298  0.0000000 -1.2838344  0.6639046
## b|b  0.0000000  0.6199298  0.3054217 -0.9253515
## 
## Dependent model rate matrix:
##            a|a        a|b        b|a         b|b
## a|a -0.5943124  0.3350400  0.2592725  0.00000000
## a|b  4.0225133 -4.0854456  0.0000000  0.06293222
## b|a  2.0872453  0.0000000 -4.5241365  2.43689116
## b|b  0.0000000  0.3640895  0.2807990 -0.64488849
## 
## Model fit:
##             log-likelihood
## independent      -171.5709
## dependent        -149.5117
## 
## Hypothesis test result:
##   likelihood-ratio:  44.11832 
##   p-value:  6.06277e-09 
## 
## Model fitting method used was fitDiscrete

That's it.

Imperceptible update to rerootingMethod for ancestral state reconstruction

I just made an update to the phytools function rerootingMethod which does marginal ancestral state estimation using the 'rerooting method' of Yang (1996). This method is essentially redundant (now) with ace(...,type="discrete") under its default settings, although it does not permit assymetric models of trait evolution (i.e., models in which the backward & forward rates for a particular transition type are allowed to assume different values). The only advantage of this method is that the function does permit polytomies. Prior versions also permitted polytomies; however this was by resolving polytomies internally, estimating ancestral states, and then matching the nodes between the fully resolved tree and its original, multifurcating counterpart. Now that rerootingMethod uses my new Mk model fitting function internally (fitMk), this is no longer necessary.

We can see this as follows:

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)
library(phytools)

## Loading required package: ape
## Loading required package: maps

set.seed(1)
## simulate tree with polytomies & data
tree<-rtree(n=26,tip.label=LETTERS)
tree$edge.length[which(tree$edge[,2]==47)]<-0
tree$edge.length[which(tree$edge[,2]==38)]<-0
tree$edge.length[which(tree$edge[,2]==29)]<-0
tree<-di2multi(tree)
plotTree(tree)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x<-sim.history(tree,Q)$states

## Done simulation(s).

x

##   S   Y   J   R   U   E   N   C   V   G   A   F   M   W   O   Q   T   B 
## "a" "c" "a" "b" "c" "c" "b" "b" "c" "a" "b" "b" "b" "b" "a" "a" "a" "a" 
##   K   X   H   P   Z   I   D   L 
## "c" "c" "a" "c" "b" "a" "b" "b"

## fit model & estimate ancestral states using rerootingMethod
model<-matrix(c(0,1,0,1,0,1,0,1,0),3,3)
rownames(model)<-colnames(model)<-letters[1:3]
fit1<-rerootingMethod(tree,x,model=model)
plotTree(tree)
nodelabels(pie=fit1$marginal.anc,piecol=setNames(c("blue","red","green"),
    c("a","b","c")),cex=0.6)
tiplabels(pie=to.matrix(x[tree$tip.label],c("a","b","c")),
    piecol=setNames(c("blue","red","green"),c("a","b","c")),cex=0.3)

## compare to ace
fit2<-ace(x,tree,type="discrete",model=model) ## doesn't work

## Error in ace(x, tree, type = "discrete", model = model): "phy" is not rooted AND fully dichotomous.

fit3<-ace(x,multi2di(tree),type="discrete",model=model) ## doesn't work

## Error in ace(x, multi2di(tree), type = "discrete", model = model): some branches have length zero or negative

tt<-multi2di(tree)
tt$edge.length[tt$edge.length==0]<-1e-8
fit4<-ace(x,tt,type="discrete",model=model) ## works
## compare to rerootingMethod
M<-matchNodes(tree,tt) ## first match nodes between the trees
plotTree(tt)
nodelabels(pie=fit4$lik.anc,piecol=setNames(c("blue","red","green"),
    c("a","b","c")),cex=0.6)
tiplabels(pie=to.matrix(x[tt$tip.label],c("a","b","c")),
    piecol=setNames(c("blue","red","green"),c("a","b","c")),cex=0.3)

plot(fit1$marginal.anc,fit4$lik.anc[M[,2]-Ntip(tt),],
    xlab="rerootingMethod",ylab="ace")

That's it.

Important updates to internal function used by make.simmap

I just rewrote the function that is used internally to compute the likelihood of any transition matrix Q given the data & tree; and to estimate the ML value of Q. This function is highly similar (but slightly different in important ways), both functionally & structurally, to ace(...,type="discrete") in the ape package. It is also in some ways functionally redundant with ace and fitDiscrete in that they fit the same model - although they differ in that the prior distribution on the root node can be controlled in this function. In fact, this is the main reason for this 're-write' because in previous versions of phytools, I suspect that this prior was being used improperly in the calculations. Now, when this function is used internally by make.simmap (which is it's main purpose, although it can also be used alone) the prior, pi, can influence not only the states sampled at the root (and thus, indirectly, for other nodes in the tree), but it can also affect the estimated value of Q under likelihood.

Here's a quick demo of using the function alone to fit the Mk discrete character evolution model:

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)

Simulate tree & data:

library(phytools)
tree<-pbtree(n=100,scale=1)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
Q

##    a  b  c
## a -1  1  0
## b  1 -2  1
## c  0  1 -1

x<-sim.history(tree,Q,anc="a")$states

## Done simulation(s).

x

##  t16  t17  t63  t64  t39  t29  t54  t55   t5  t18  t50  t51  t33  t19  t20 
##  "a"  "a"  "a"  "a"  "a"  "a"  "c"  "b"  "c"  "b"  "c"  "c"  "c"  "b"  "c" 
##  t10  t99 t100  t88  t76  t79  t80  t93  t94  t89  t90  t44  t95  t96  t68 
##  "c"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "c" 
##  t69  t65  t81  t82  t77  t78  t60   t9  t40  t83  t84  t72  t73  t56  t57 
##  "c"  "c"  "c"  "c"  "c"  "c"  "c"  "a"  "c"  "b"  "b"  "b"  "b"  "c"  "c" 
##  t26  t21  t61  t62  t42  t43  t52  t53   t6   t7   t2  t58  t59  t11  t12 
##  "b"  "b"  "b"  "b"  "b"  "b"  "c"  "c"  "a"  "a"  "a"  "a"  "a"  "a"  "a" 
##   t1  t22  t23   t3   t4  t45  t46  t27  t85  t86  t35  t36  t37  t38  t28 
##  "a"  "c"  "c"  "b"  "a"  "a"  "a"  "a"  "a"  "a"  "b"  "b"  "b"  "a"  "a" 
##  t24  t25  t30  t74  t75  t14  t15  t13   t8  t87  t91  t92  t70  t71  t49 
##  "a"  "b"  "b"  "a"  "b"  "a"  "a"  "a"  "a"  "b"  "b"  "b"  "b"  "b"  "b" 
##  t47  t48  t41  t34  t66  t67  t31  t32  t97  t98 
##  "b"  "a"  "b"  "b"  "b"  "b"  "a"  "b"  "a"  "a"

Fit the model:

obj<-fitMk(tree,as.factor(x),model="SYM")
obj

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b        c
## a -0.930216  0.930216  0.00000
## b  0.930216 -2.094586  1.16437
## c  0.000000  1.164370 -1.16437
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -66.577856

It is also used internally by make.simmap. Here, with a prior that is not flat:

trees<-make.simmap(tree,x,pi=setNames(c(1,0,0),letters[1:3]),nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a          b        c
## a -0.8965016  0.8965016  0.00000
## b  0.8965016 -2.0898515  1.19335
## c  0.0000000  1.1933498 -1.19335

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

## a b c 
## 1 0 0

## Done.

obj<-summary(trees)
obj

## 100 trees with a mapped discrete character with states:
##  a, b, c 
## 
## trees have 34.17 changes between states on average
## 
## changes are of the following types:
##      a,b a,c  b,a  b,c c,a  c,b
## x->y  12   0 6.36 9.63   0 6.18
## 
## mean total time spent in each state is:
##               a         b         c    total
## raw  12.1857664 8.5412422 4.3792713 25.10628
## prop  0.4853673 0.3402034 0.1744293  1.00000

plot(obj,ftype="off",type="fan")

That's it.

Print method for fastAnc

I just added a print method for the ancestral state estimation function, fastAnc. Here is a quick demo of how it works:

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)
library(phytools)

## Loading required package: ape
## Loading required package: maps

Simulate a tree & data:

x<-fastBM(tree<-rtree(n=12))
tree

## 
## Phylogenetic tree with 12 tips and 11 internal nodes.
## 
## Tip labels:
##  t5, t6, t8, t7, t3, t4, ...
## 
## Rooted; includes branch lengths.

x

##          t5          t6          t8          t7          t3          t4 
##  0.08020791  0.90090545 -0.52120123 -0.70089226  0.60523627  0.32977174 
##         t10          t2         t12         t11          t1          t9 
## -1.35329157 -0.41231935 -0.12192103  0.56163234 -1.40509098 -1.52915450

obj<-fastAnc(tree,x)
obj

## Ancestral character estimates using fastAnc:
##        13        14        15        16        17        18        19 
## -0.613584 -0.333720  0.311369 -0.667771 -0.805572 -0.640828  0.133491 
##        20        21        22        23 
## -0.665823 -1.137190 -0.402258  0.504335

obj<-fastAnc(tree,x,CI=TRUE)
obj

## Ancestral character estimates using fastAnc:
##        13        14        15        16        17        18        19 
## -0.613584 -0.333720  0.311369 -0.667771 -0.805572 -0.640828  0.133491 
##        20        21        22        23 
## -0.665823 -1.137190 -0.402258  0.504335 
## 
## Lower & upper 95% CIs:
##        lower     upper
## 13 -1.872864  0.645696
## 14 -1.649728  0.982288
## 15 -0.648504  1.271242
## 16 -1.773040  0.437499
## 17 -1.911716  0.300572
## 18 -1.048660 -0.232996
## 19 -0.697674  0.964656
## 20 -1.724715  0.393068
## 21 -1.890327 -0.384054
## 22 -1.274650  0.470134
## 23  0.242778  0.765892

print(obj,printlen=6)

## Ancestral character estimates using fastAnc:
##         13       14       15        16        17        18     
##  -0.613584 -0.33372 0.311369 -0.667771 -0.805572 -0.640828 ....
## 
## Lower & upper 95% CIs:
##        lower     upper
## 13 -1.872864  0.645696
## 14 -1.649728  0.982288
## 15 -0.648504  1.271242
## 16  -1.77304  0.437499
## 17 -1.911716  0.300572
## 18  -1.04866 -0.232996
##         ....      ....

That's it.