Update to locate.yeti: REML and exact likelihood methods

I just posted code for a new version of the function locate.yeti. This function can be used to place recently extinct, missing, or cryptic taxa into an ultrametric phylogeny based on continuous character data.

The main updates are as follows:

(1) The function now includes a REML method which uses the contrasts algorithm of Felsenstein (1985). Note that this option is currently under development.

(2) The function also now includes exact likelihood calculation, rather then the approximate method that used orthogonalization of the original data based on the backbone tree. This obviously slows down computation, and I implemented it primarily to compare to the REML method (which does not need to use orthogonalization). It seems to not be prohibitively slow for relatively modest sized trees (say, 100 taxa or so).

Here's a quick demo.

First, simulate tree & data:

library(phytools)

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

library(phangorn)

## 
## Attaching package: 'phangorn'
## 
## The following object is masked from 'package:ape':
## 
##     as.prop.part

## simulate tree & data
N<-50 ## taxa in base tree
m<-10 ## number of continuous characters
## simulate tree
tt<-tree<-pbtree(n=N+1,tip.label=sample(c(paste("t",1:N,sep=""),"Yeti")))
## generate a covariance matrix for simulation
L<-matrix(rnorm(n=m*m),m,m)
L[upper.tri(L,diag=FALSE)]<-0
L<-L-diag(diag(L))+abs(diag(diag(L)))
V<-L%*%t(L)
X<-sim.corrs(tree,vcv=V)
tree<-drop.tip(tree,"Yeti")
## visualize trees
par(mfrow=c(1,2))
plotTree(tt,mar=c(0.1,0.1,4.1,0.1),fsize=0.8)
title("tree with Yeti")
plotTree(tree,mar=c(0.1,0.1,4.1,0.1),direction="leftwards",fsize=0.8)
title("tree without Yeti")

OK, now let's estimate using the ML method with & without rotation:

## without rotation
mltree<-locate.yeti(tree,X,plot=TRUE,search="exhaustive",rotate=FALSE)

## Optimizing the phylogenetic position of Yeti using ML. Please wait....

## Done.

mltree.rotate<-locate.yeti(tree,X,plot=TRUE,search="exhaustive")

## Optimizing the phylogenetic position of Yeti using ML. Please wait....

## Done.

mltree$logL

## [1] -455.4

mltree.rotate$logL

## [1] -456.3

RF.dist(mltree,mltree.rotate)

## [1] 0

Hopefully, we get more or less the same tree in both cases!

Now, let's try the REML method. We can compare our results to the true tree and to our ML tree from above:

remltree<-locate.yeti(tree,X,method="REML",plot=TRUE,search="exhaustive")

## ---------------------------------------------------------------
## | **Warning: method="REML" has not been thoroughly tested.    |
## |   Use with caution.**                                       |
## ---------------------------------------------------------------
## 
## Optimizing the phylogenetic position of Yeti using REML. Please wait....

## Done.

RF.dist(tt,remltree)

## [1] 4

RF.dist(mltree,remltree)

## [1] 0

In addition to the RF distance, let's compute the patristic distances on the ML and REML trees. Note that it doesn't mean much that most of these are precisely the same - these are distances that don't involve the unknown tip. We should instead compare the very few outliers that involve distances from our unknown tip to other taxa in the tree:

par(mfrow=c(1,2))
plot(cophenetic(tt)[tt$tip.label,tt$tip.label],
    cophenetic(mltree)[tt$tip.label,tt$tip.label],
    xlab="true distances",ylab="ML tree distances")
plot(cophenetic(tt)[tt$tip.label,tt$tip.label],
    cophenetic(remltree)[tt$tip.label,tt$tip.label],
    xlab="true distances",ylab="REML tree distances")

OK, well, that's all I have. This update, along with the new phytools function fitPagel, is in a new non-CRAN phytools version (phytools 0.4-40) that I recently posted. Let me know if you run into any difficulties with either function.

Simulating correlated evolution of discrete characters under Pagel's model

In a recent post, in which I gave an R function to run Pagel's (1994) test for correlation between two binary traits, I also posted (without much comment) some simple code to simulate the character histories of binary traits undergoing correlated evolution by Pagel's model.

Here is that code again with some explanation.

(1) First, just as a preliminary, let's load our packages & simulate a stochastic phylogeny:

library(phytools)

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

tree<-pbtree(n=300,scale=1)

(2) Next, we need to set up a matrix to simultaneously simulate our two traits. Remember, under Pagel's model, character correlation are just transition rates in one character that depend on the state of a second. To simulate under this model we actually have to simultaneously simulate states for the two characters at once. Here, as I did last time, I'll keep it simple. Let's say that the two characters have states a and b. I'll let the transition rates into states a,a or b,b be high (say 2.0), but I will make the transition rates to any state a,b or b,a low - say 0.4. This is what that looks like:

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)
Q

##      aa   ab   ba   bb
## aa -0.8  0.4  0.4  0.0
## ab  2.0 -4.0  0.0  2.0
## ba  2.0  0.0 -4.0  2.0
## bb  0.0  0.4  0.4 -0.8

(3) Now, we can simulate both characters simultaneously up the tree using the phytools function sim.history.

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).

(4) This is the history of our two character state data; however we still have to backtranslate this into the character histories for each of our two traits. To do this, we can use mergeMappedStates to merge a,a and a,b (for instance) into a for character 1 and so on. Here is what that looks like:

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")

(5) For fun, let's plot the two histories (just as we did last time) so that we can see that they are indeed highly correlated:

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")

(6) Finally, we can fit the Pagel (1994) model:

x<-getStates(t1,"tips")
y<-getStates(t2,"tips")
source("fitPagel.R")
.check.pkg<-phytools:::.check.pkg
fit<-fitPagel(tree,x,y)
fit

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -0.9753616  0.5606835  0.4146781  0.0000000
## a|b  2.7557800 -3.1704581  0.0000000  0.4146781
## b|a  2.7479439  0.0000000 -3.3086275  0.5606835
## b|b  0.0000000  2.7479439  2.7557800 -5.5037240
## 
## Dependent model rate matrix:
##            a|a        a|b        b|a       b|b
## a|a -0.7760794  0.4795664  0.2965130  0.000000
## a|b  5.1666938 -7.2422057  0.0000000  2.075512
## b|a  2.0876921  0.0000000 -4.0078768  1.920185
## b|b  0.0000000  3.8020580  0.7652618 -4.567320
## 
## Model fit:
##             log-likelihood
## independent      -186.3963
## dependent        -173.9505
## 
## Hypothesis test result:
##   likelihood-ratio:  24.89146 
##   p-value:  5.290212e-05

That's it. Hopefully this is helpful to someone!

R function for Pagel's 1994 correlation method

I just posted some code to fit the Pagel (1994) method for detecting “correlated” evolution of two binary traits. I say correlated because it might be more appropriate to call this a dependence - i.e., the rate of character 2 depends on 1 & vice versa.

This is actually quite simple to do because the test is merely a test of, for a two character combination written as [01] (meaning state 0 for the first character & state 1 for the second), the rates of transition from [00]->[01] equal [10]->[11]; [00]->[10] equal [01]->[11]; etc. For this we can just paste our two input characters together, fit both a model in which the dependence exists and one in which it does not, and then compare them. The dependence model has twice as many transition rates as the independence model, so this model will have (for two binary traits) 4 more parameters.

Note that I accomplished this by borrowing ace in the ape package or fitDiscrete in the geiger package to actually fit the models. All my code really does is pre- & post-process the data & results respectively!

Here's a demo:

## first load packages & source code
library(phytools)
library(geiger)
source("fitPagel.R")
.check.pkg<-phytools:::.check.pkg
## now let's simulate some uncorrelated data
tree<-pbtree(n=300,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
tt1<-sim.history(tree,Q)

## Done simulation(s).

tt2<-sim.history(tree,Q)

## Done simulation(s).

## these are uncorrelated, see:
par(mfrow=c(1,2))
plotSimmap(tt1,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1)
plotSimmap(tt2,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1,direction="leftwards")

x<-tt1$states
y<-tt2$states
fit.ape<-fitPagel(tree,x,y)
fit.ape

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.9712  1.0896  0.8816  0.0000
## a|b  0.9774 -1.8589  0.0000  0.8816
## b|a  0.7531  0.0000 -1.8426  1.0896
## b|b  0.0000  0.7531  0.9774 -1.7304
## 
## Dependent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.5606  0.8578  0.7028  0.0000
## a|b  0.9261 -1.8634  0.0000  0.9373
## b|a  0.3652  0.0000 -1.5764  1.2112
## b|b  0.0000  1.2754  1.0567 -2.3321
## 
## Model fit:
##             log-likelihood
## independent         -238.7
## dependent           -237.5
## 
## Hypothesis test result:
##   likelihood-ratio:  2.390211 
##   p-value:  0.6643971

fit.geiger<-fitPagel(tree,x,y,method="fitDiscrete")

## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41
## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

fit.geiger

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.8847  0.9942  0.8906  0.0000
## a|b  0.9699 -1.8605  0.0000  0.8906
## b|a  0.7112  0.0000 -1.7053  0.9942
## b|b  0.0000  0.7112  0.9699 -1.6811
## 
## Dependent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.5000  0.8192  0.6809  0.0000
## a|b  0.8987 -1.8596  0.0000  0.9609
## b|a  0.3704  0.0000 -1.4820  1.1116
## b|b  0.0000  1.1176  1.1075 -2.2251
## 
## Model fit:
##             log-likelihood
## independent         -239.3
## dependent           -238.3
## 
## Hypothesis test result:
##   likelihood-ratio:  2.096727 
##   p-value:  0.7179738

OK, now let's try correlated data. To do that I will have to build a big transition matrix myself for simulation!

Q<-matrix(c(0,0.5,0.5,0,2,0,0,2,2,0,0,2,0,0.5,0.5,0),4,4,byrow=TRUE)
rownames(Q)<-colnames(Q)<-c("aa","ab","ba","bb")
diag(Q)<--rowSums(Q)
tt<-sim.history(tree,t(Q))

## Note - the rate of substitution from i->j should be given by Q[j,i].
## Done simulation(s).

tt1<-mergeMappedStates(tt,c("aa","ab"),"a")
tt1<-mergeMappedStates(tt1,c("ba","bb"),"b")
tt2<-mergeMappedStates(tt,c("aa","ba"),"a")
tt2<-mergeMappedStates(tt2,c("ab","bb"),"b")
## these data are correlated, see:
par(mfrow=c(1,2))
plotSimmap(tt1,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1)
plotSimmap(tt2,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1,direction="leftwards")

x<-getStates(tt1,"tips")
y<-getStates(tt2,"tips")
fit.ape<-fitPagel(tree,x,y)
fit.ape

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.4325  0.9831  0.4493  0.0000
## a|b  0.9593 -1.4086  0.0000  0.4493
## b|a  1.0401  0.0000 -2.0233  0.9831
## b|b  0.0000  1.0401  0.9593 -1.9994
## 
## Dependent model rate matrix:
##        a|a     a|b     b|a    b|b
## a|a -1.043  0.8539  0.1893  0.000
## a|b  5.543 -7.9019  0.0000  2.359
## b|a  1.202  0.0000 -4.4114  3.209
## b|b  0.000  1.0849  0.4054 -1.490
## 
## Model fit:
##             log-likelihood
## independent         -241.1
## dependent           -207.8
## 
## Hypothesis test result:
##   likelihood-ratio:  66.63346 
##   p-value:  1.164753e-13

fit.geiger<-fitPagel(tree,x,y,method="fitDiscrete")

## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41
## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

fit.geiger

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a    a|b     b|a     b|b
## a|a -1.3892  1.006  0.3832  0.0000
## a|b  0.9236 -1.307  0.0000  0.3832
## b|a  1.0639  0.000 -2.0699  1.0061
## b|b  0.0000  1.064  0.9236 -1.9875
## 
## Dependent model rate matrix:
##         a|a     a|b        b|a    b|b
## a|a -0.8951  0.8951  1.799e-16  0.000
## a|b  5.4502 -7.7907  0.000e+00  2.340
## b|a  1.9288  0.0000 -3.958e+00  2.029
## b|b  0.0000  1.0086  4.105e-01 -1.419
## 
## Model fit:
##             log-likelihood
## independent         -241.9
## dependent           -208.8
## 
## Hypothesis test result:
##   likelihood-ratio:  66.16744 
##   p-value:  1.460391e-13

That's it, I guess. This is brand new so I have no idea how it works but feedback is welcome. I'm also well aware of the recent critique of this method by Maddision & Fitzjohn.

Testing a hypothesis that phylogenetic signal is different from 1.0

A phytools user asked today:

“A reviewer asked us to test if phylogenetic signal is equal to 1, but I have not found any package that does that. Do you have any suggestion about how to do this?

Since phylogenetic signal (measured using either Pagel's λ or Blomberg's K) has an expected value of 1.0 under Brownian motion; testing whether signal is different from 1.0 is essentially testing whether phylogenetic signal is smaller or greater than if the data arose by a Brownian process.

We can do this with a likelihood-ratio test (with some caveats) if we are using Pagel's λ to measure phylogenetic signal. If we are using Blomberg's K then we can use simulation (with fewer caveats).

Here's a quick demo.

First, let's simulate some data with & without signal. For the no-signal data I will just randomly permute the Brownian data:

library(phytools)
## tree
tree<-pbtree(n=26,tip.label=LETTERS)
## data with signal
x1<-fastBM(tree)
## data without signal
x2<-setNames(sample(x1),names(x1))

Now let's estimate Pagel's lambda for both datasets. In each case, we will compare the likelihood of the fitted model to the likelihood of λ = 1.

## Pagel's lambda
## first the data with signal
lam1<-phylosig(tree,x1,method="lambda")
lam1

## $lambda
## [1] 0.9413
## 
## $logL
## [1] -38.39

fitBrownian<-brownie.lite(paintSubTree(tree,tree$edge[1,1],state="1"),x1)
fitBrownian

## ML single-rate model:
##  s^2 se  a   k   logL
## value    1.1035  0.3061  -0.2291 2   -38.6817    
## 
## ML multi-rate model:
##  s^2(1)  se(1)   a   k   logL    
## value    1.1035      -0.2291 2   -38.6817
## 
## P-value (based on X^2): 1 
## 
## R thinks it has found the ML solution.

## likelihood ratio test
LR<-2*(lam1$logL-fitBrownian$logL1)
LR

## [1] 0.5774

P.lr<-pchisq(LR,df=1,lower.tail=FALSE)
P.lr

## [1] 0.4473

## now the data without signal
lam2<-phylosig(tree,x2,method="lambda")
lam2

## $lambda
## [1] 6.917e-05
## 
## $logL
## [1] -46.33

fitBrownian<-brownie.lite(paintSubTree(tree,tree$edge[1,1],state="1"),x2)
fitBrownian

## ML single-rate model:
##  s^2 se  a   k   logL
## value    4.5548  1.2633  0.0452  2   -57.112 
## 
## ML multi-rate model:
##  s^2(1)  se(1)   a   k   logL    
## value    4.5548  1   0.0452  2   -57.112
## 
## P-value (based on X^2): 1 
## 
## Optimization may not have converged.

## likelihood ratio test
LR<-2*(lam2$logL-fitBrownian$logL1)
LR

## [1] 21.56

P.lr<-pchisq(LR,df=1,lower.tail=FALSE)
P.lr

## [1] 3.428e-06

Now let's do roughly the same thing for Blomberg's K. In this case, however, we will use simulation by Brownian motion to get a null distribution for K.

## Blomberg's K
## first the data with signal
K1<-phylosig(tree,x1)
K1

## [1] 0.9221

## fit a Brownian model to the data to get parameters
## for simulation
fitBrownian<-brownie.lite(paintSubTree(tree,tree$edge[1,1],state="1"),x1)
X<-fastBM(tree,sig2=fitBrownian$sig2.single,a=fitBrownian$a.single,nsim=999)
nullK<-apply(X,2,phylosig,tree=tree)
mean(nullK)

## [1] 0.9857

hist(nullK,breaks=20,xlab="null distribution for K",
    main="Null distribution of K")
lines(c(K1,K1),c(0,par()$usr[4]),lty="dashed",col="red")
text(x=K1,y=0.985*par()$usr[4],"observed value of K",pos=4,offset=0.2)

The distribution is of K under the null hypothesis is highly assymetric, so it might make sense to conduct a test of K using abs(log(K)). (Just an idea.)

P.k<-mean(abs(log(c(K1,nullK)))>=abs(log(K1)))
P.k

## [1] 0.849

Now let's do the same with the no-signal data:

K2<-phylosig(tree,x2)
K2

## [1] 0.2261

## fit a Brownian model to the data to get parameters
## for simulation
fitBrownian<-brownie.lite(paintSubTree(tree,tree$edge[1,1],state="1"),x2)
X<-fastBM(tree,sig2=fitBrownian$sig2.single,a=fitBrownian$a.single,nsim=999)
nullK<-apply(X,2,phylosig,tree=tree)
mean(nullK)

## [1] 0.9844

hist(nullK,breaks=20,xlab="null distribution for K",
    main="Null distribution of K")
lines(c(K2,K2),c(0,par()$usr[4]),lty="dashed",col="red")
text(x=K2,y=0.985*par()$usr[4],"observed value of K",pos=4,offset=0.2)

P.k<-mean(abs(log(c(K2,nullK)))>=abs(log(K2)))
P.k

## [1] 0.001

That's it.