Extension of ratebytree for discretely valued characters

As I mentioned the other day, during my sabbatical at the Universidad de los Andes here in Bogotá I submitted a paper with some colleagues describing the function ratebytree for comparing the evolutionary rates of a continuous character between trees.

The manuscript was favorably reviewed, but in their typically pesky* (*said with affection) manner, reviewers asked us to add different models of evolution - the feature that I added last week. They also asked us to extend the approach to discrete character evolution. This was something we pointed out could be done in the manuscript, but had not explicitly addressed. It is theoretically straightforward - although not completely trivial - particularly if you want your new function updates to be accompanied by legible print methods & so on.

I have, however, now added this feature to the function ratebytree. The function can take an argument, type (which should assume the value "continuous" or "discrete"); however, if not provided, the function will attempt to ascertain the character type.

Here is a quick demo using some simple simulated data:

library(phytools)
packageVersion("phytools")

## [1] '0.6.25'

First, our trees & data:

t1

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x1

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

t2

## 
## Phylogenetic tree with 40 tips and 39 internal nodes.
## 
## Tip labels:
##  t7, t34, t35, t3, t29, t30, ...
## 
## Rooted; includes branch lengths.

x2

##  t7 t34 t35  t3 t29 t30 t31 t14 t15  t2 t39 t40  t9  t8 t21 t22 t32 t33 
##   b   a   a   c   b   b   b   b   c   b   a   a   c   c   a   a   c   c 
## t23 t24 t19 t20 t11  t5  t6 t25 t26  t4 t36 t37 t38 t12 t13 t27 t28 t16 
##   c   c   b   b   c   c   a   c   c   a   b   b   b   c   c   c   b   c 
##  t1 t10 t17 t18 
##   a   a   a   a 
## Levels: a b c

## or
dotTree(t3,x3,ftype="off")

x<-list(x1,x2,x3)
x

## [[1]]
## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## c c c c a c b c b c a c b a b c b b a c a c c c c c 
## Levels: a b c
## 
## [[2]]
##  t7 t34 t35  t3 t29 t30 t31 t14 t15  t2 t39 t40  t9  t8 t21 t22 t32 t33 
##   b   a   a   c   b   b   b   b   c   b   a   a   c   c   a   a   c   c 
## t23 t24 t19 t20 t11  t5  t6 t25 t26  t4 t36 t37 t38 t12 t13 t27 t28 t16 
##   c   c   b   b   c   c   a   c   c   a   b   b   b   c   c   c   b   c 
##  t1 t10 t17 t18 
##   a   a   a   a 
## Levels: a b c
## 
## [[3]]
## t19 t20  t3 t11 t12 t36 t37 t21 t24 t25 t10 t13 t14  t2  t7  t8  t9 t16 
##   a   a   a   b   a   c   b   a   a   a   b   a   a   a   c   b   a   c 
## t28 t29 t32 t33 t17 t34 t35 t18  t1 t39 t40 t22 t23 t41 t42  t4  t5 t26 
##   b   b   b   b   b   b   b   b   b   c   c   c   c   b   b   b   a   b 
## t27  t6 t49 t50 t15 t43 t44 t45 t46 t30 t31 t47 t48 t38 
##   b   b   c   c   b   c   c   b   b   c   c   b   b   b 
## Levels: a b c

trees<-c(t1,t2,t3)
trees

## 3 phylogenetic trees

OK, now let's try it! First, I will assume a symmetric model of transitions between states:

fitSYM<-ratebytree(trees,x,type="discrete",model="SYM")
fitSYM

## ML common-rate model:
##          a<->b   a<->c   b<->c   k   logL
## value    1.3136  0.0815  3.0745  3   -100.7769
## 
## Model fitting method was "optim".
## 
## ML multi-rate model:
##          a<->b   a<->c   b<->c   k   logL
## tree1    3.5598  0   5.2209  9   -99.2075
## tree2    0.2588  0.5751  1.5436
## tree3    1.3966  0   3.2116
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 3.1388 
## P-value (based on X^2): 0.7912

This tells us that we have little evidence to reject a constant (set of) rate(s) across our trees. We can nonetheless try simpler or more complex models and see if the outcome is similar or different. Remember, I said the function would attempt to ascertain if the character is discretely or continuously-valued, so let's try that:

fitER<-ratebytree(trees,x)
fitER

## ML common-rate model:
##          q       k   logL
## value    1.1995  1   -104.2243
## 
## Model fitting method was "optimize".
## 
## ML multi-rate model:
##          q       k   logL
## tree1    2.1052  3   -103.4967
## tree2    0.8053
## tree3    1.4482
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 1.4552 
## P-value (based on X^2): 0.4831

fitARD<-ratebytree(trees,x,model="ARD")
fitARD

## ML common-rate model:
##          b->a    c->a    a->b    c->b    a->c    b->c    k   logL
## value    -0.0483 1.1847  2.3293  1.5367  -0.7321 3.8122  6   -99.4653
## 
## Model fitting method was "optim".
## 
## ML multi-rate model:
##          b->a    c->a    a->b    c->b    a->c    b->c    k   logL
## tree1    4.0703  0   4.1987  7.5158  0   16.8153 18  -94.8226
## tree2    0.1143  0.5863  0   1.2591  0.1455  2.3226
## tree3    0.6492  0.2584  1.9163  4.7043  0   2.7484
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 9.2853 
## P-value (based on X^2): 0.6784

By chance, these data were simulated as follows - and thus without a difference in rate or process between trees:

t1<-pbtree(n=26,tip.label=LETTERS,scale=1)
t2<-pbtree(n=40,scale=2)
t3<-pbtree(n=50,scale=1)
Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x1<-as.factor(sim.history(t1,Q)$states)
x2<-as.factor(sim.history(t2,Q)$states)
x3<-as.factor(sim.history(t3,Q)$states)

Why don't we simulate under conditions in which the process differs between trees and see if we obtain a different result:

t4<-pbtree(n=60,scale=1)
t5<-pbtree(n=70,scale=2)
t6<-pbtree(n=80,scale=1)
Q1<-matrix(c(-1,1,1,-1),2,2)
Q2<-matrix(c(-2,2,2,-2),2,2)
x4<-as.factor(sim.history(t4,Q1)$states)
x5<-as.factor(sim.history(t5,Q1)$states)
x6<-as.factor(sim.history(t6,Q2)$states)

Now we're ready to fit our models:

fitER<-ratebytree(c(t4,t5,t6),list(x4,x5,x6))
fitER

## ML common-rate model:
##          q       k   logL
## value    1.2509  1   -116.4188
## 
## Model fitting method was "optimize".
## 
## ML multi-rate model:
##          q       k   logL
## tree1    0.8768  3   -112.826
## tree2    0.7286
## tree3    2.7957
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 7.1856 
## P-value (based on X^2): 0.0275

This suggests that a model in which the different trees are permitted to have different rates of evolution fits better than one in which they are all constrained to have equal rates - just as we simulated.

Finally, to get this functionality you should update phytools from GitHub as follows:

library(devtools)
install_github("liamrevell/phytools")

More later.

Identifying all the tips consisting of all monophyletic clades in a tree with 11 taxa

A colleague just asked me the following:

“My question: if I have a particular number of taxa that I want to find a monophyletic group containing only that number of taxa out of a larger tree, is there such a function? For example, I have a tree of 100 taxa and I want a list (or lists) of exactly 11 taxa sharing an mrca, above that mrca should only contain 11 taxa and not more.”

This could be done multiple ways. Here's one using ape::extract.clade:

library(phytools)
set.seed(3)
tree<-pbtree(n=100)
nodes<-1:tree$Nnode+Ntip(tree)
Ndesc<-sapply(nodes,function(x,tree) Ntip(extract.clade(tree,x)),tree=tree)

plotTree(tree,ftype="off",lwd=1)
labelnodes(Ndesc,node=nodes,interactive=F,cex=0.5)

nodes11<-nodes[which(Ndesc==11)]
nodes11

## [1] 116 134 163

nodelabels(Ndesc[which(Ndesc==11)],node=nodes11,cex=0.8,frame="circle",
    bg="red",col="white")

tips<-sapply(nodes11,function(x,tree) extract.clade(tree,x)$tip.label,
    tree=tree)
colnames(tips)<-nodes11
tips

##       116   134   163  
##  [1,] "t47" "t37" "t72"
##  [2,] "t48" "t95" "t73"
##  [3,] "t36" "t96" "t32"
##  [4,] "t89" "t42" "t43"
##  [5,] "t90" "t74" "t44"
##  [6,] "t50" "t75" "t88"
##  [7,] "t52" "t82" "t97"
##  [8,] "t53" "t83" "t98"
##  [9,] "t57" "t69" "t7" 
## [10,] "t58" "t70" "t59"
## [11,] "t49" "t38" "t60"

ratebytree for comparing rates & regimes among trees with different evolutionary models

Along with several Universidad de los Andes colleagues from my sabbatical there (such as the famous Andrew Crawford and his postdoc Carlos Guarnizo), I recently submitted for publication a description of the phytools function ratebytree for comparing the evolutionary rate between trees - along with some of its statistical properties.

This method basically implements the censored test of O'Meara et al. (2006), but applied to the problem of comparing the rate of evolution between trees - rather than among clades or splits in a given tree.

One of the editor comments was that it would be straightforward to allow the comparison of different models or processes of evolution on the different trees. Although I prefer the simplicity of interpretation of a comparison of rates (a difference in which may be due to many different processes), the editor is right in that it is also possible to fit different models of trait evolution to the different trees & then ask if permitting the parameters of these models to take different values significantly better explains our data (or not).

So far, I have implemented the models "OU" (Ornstein-Uhlenbeck) and "EB" (early-burst). The OU model in comparative analyses is most often associated with Hansen (1997), although I believe Felsenstein suggested it earlier; and the EB model is linked to Blomberg et al. (2003).

Here's how it works, first with the "EB" model:

library(phytools)
packageVersion("phytools")

## [1] '0.6.24'

t1

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 40 tips and 39 internal nodes.
## 
## Tip labels:
##  t1, t4, t23, t24, t18, t16, ...
## 
## Rooted; includes branch lengths.

x1

##          A          B          C          D          E          F 
## -0.3441420 -0.3533002 -0.3115814 -0.3160447 -0.4460135 -0.4550335 
##          G          H          I          J          K          L 
## -0.4139192 -0.3953837 -0.4060349 -0.4077641 -0.8347388 -0.7452221 
##          M          N          O          P          Q          R 
## -0.7185576 -0.7635412 -0.7700484 -0.7695936 -0.7634808 -0.6948223 
##          S          T          U          V          W          X 
## -0.7412721 -0.8551189 -0.7677177 -0.7679566 -0.8786428 -0.9301932 
##          Y          Z 
## -0.9146077 -0.8141308

x2

##         t1         t4        t23        t24        t18        t16 
##  0.8566059  3.7251483  2.6191894  1.8225781  1.8057385  3.4079308 
##        t29        t30        t11        t12         t5        t10 
##  4.5958938  4.1774617  4.2133663  3.7181830  3.6334797  4.7123706 
##        t25        t26         t7        t17        t31        t32 
##  5.2011408  5.3905739  4.1595605  2.0951849  2.9624926  3.4058012 
##        t35        t36        t19        t15         t6        t27 
##  4.2541885  4.1980012  3.8449273  4.3596184  4.4675526  3.1422765 
##        t28        t20        t39        t40         t8        t13 
##  3.8031220  2.5404983  4.9569045  4.8507862  4.3520308  4.5914069 
##        t14         t2        t33        t34        t37        t38 
##  5.2167261  1.8276020  1.0305605  0.2938116  0.7720743  0.9339237 
##         t3        t21        t22         t9 
## -0.5178370  1.1703200  0.9172244  0.6591119

fitEB<-ratebytree(c(t1,t2),list(x1,x2),model="EB")
fitEB

## ML common-regime EB model:
##  s^2  a[1]   a[2]    r   k   logL
## value    12.0073 -0.7814 1.5506  -1.2953 4   -57.2681    
## 
## ML multi-regime EB model:
##   s^2[1] s^2[2]   a[1]   a[2]     r[1]   r[2]    k   logL
## value    0.3448  3.4066  -0.7853 1.6983  -1.5952 -0.6047 6   -8.2072 
## 
## Likelihood ratio: 98.1218 
## P-value (based on X^2): 0 
## 
## R thinks it has found the ML solution.

This tells us that under the "EB" model if we allow each tree to have different parameter values the fit is significantly better than if they are both constrained to have the same parameter values.

Next, we can try the "OU" model:

x3

##          A          B          C          D          E          F 
##  2.0994227  2.0572285  3.0075786  2.8504237  3.6951966  3.6323817 
##          G          H          I          J          K          L 
##  2.4724997  1.6972401  1.4710967  1.3452084  1.0863732  2.3565924 
##          M          N          O          P          Q          R 
##  3.4456653 -0.4734070  0.8027163  0.5108661  1.4534683  0.1115656 
##          S          T          U          V          W          X 
## -1.4205819  2.1651582  2.8256807  2.9683889  2.8473970 -1.1233798 
##          Y          Z 
## -1.5167766 -3.1592647

x4

##            t1            t4           t23           t24           t18 
## -0.4781630036 -0.4130760100 -0.7107502478  0.4403283556 -0.9813284805 
##           t16           t29           t30           t11           t12 
##  0.1100536930  0.1748046581 -0.0005799426 -0.4965174644  0.6461984642 
##            t5           t10           t25           t26            t7 
## -0.2498849416 -0.3849844114  0.1414963280 -0.2731847500 -0.3726606412 
##           t17           t31           t32           t35           t36 
##  0.1467272989 -0.4809644865 -0.0103007857  0.6234118866 -0.4649965648 
##           t19           t15            t6           t27           t28 
## -0.5483876423 -0.5286145444  0.6421465052  1.1954369835  0.4204277069 
##           t20           t39           t40            t8           t13 
##  1.2252202112 -0.2934648298 -0.4060123937  0.1507545508  0.5722774758 
##           t14            t2           t33           t34           t37 
##  0.3287662081  0.1404059566 -0.6361731193 -0.8295810175  0.6294745098 
##           t38            t3           t21           t22            t9 
##  0.7020644857 -0.6360564202  1.2895073650 -0.0730771239 -0.1221469493

fitOU<-ratebytree(c(t1,t2),list(x1,x2),model="OU")
fitOU

## ML common-regime OU model:
##  s^2  a[1]   a[2]    alpha   k   logL
## value    0.4868  -0.7499 1.9136  0   4   -65.6359    
## 
## ML multi-regime OU model:
##   s^2[1] s^2[2]   a[1]   a[2]     alp[1] alp[2]  k   logL
## value    0.0052  0.796   -0.7499 1.9136  -0.306  -0.1318 6   -16.5559    
## 
## Likelihood ratio: 98.1601 
## P-value (based on X^2): 0 
## 
## R thinks it has found the ML solution.

Once again, our fit allowing different regime parameterization in the two trees is significantly better than when they are forced to be the same.

(Note that here t1 has a negative value of α. Normally we wouldn't allow that, so I will have to fix that later.)

Finally, another example with "OU":

fitOU2<-ratebytree(c(t1,t2),list(x5,x6),model="OU")
fitOU2

## ML common-regime OU model:
##  s^2  a[1]   a[2]    alpha   k   logL
## value    1.0712  0.1545  0.092   1.4324  4   -54.3758    
## 
## ML multi-regime OU model:
##   s^2[1] s^2[2]   a[1]   a[2]     alp[1] alp[2]  k   logL
## value    1.658   0.7792  0.1563  0.0949  1.8076  1.2412  6   -53.373 
## 
## Likelihood ratio: 2.0056 
## P-value (based on X^2): 0.3669 
## 
## R thinks it has found the ML solution.

In this case, we cannot reject the possibility that both traits are evolving under the same OU process.

FYI, the trees & data for this example were simulated as follows:

t1<-pbtree(n=26,tip.label=LETTERS)
t2<-pbtree(n=40)
x1<-fastBM(phytools:::ebTree(t1,-2))
x2<-fastBM(t2)
x3<-fastBM(t1)
x4<-fastBM(t2,model="OU",alpha=2)
x5<-fastBM(t1,model="OU",alpha=1)
x6<-fastBM(t2,model="OU",alpha=1)

As an addendum, I thought I'd note the following things.

Firstly, the function can also be used in basically the same way as fitContinuous if we supply just a single tree & trait vector instead of a list of each. For instance:

obj<-list(t1)
class(obj)<-"multiPhylo"
fitBM1<-ratebytree(obj,list(x3))
fitContinuous(t1,x3)

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
##  sigsq = 1.059273
##  z0 = 0.540955
## 
##  model summary:
##  log-likelihood = -37.111074
##  AIC = 78.222148
##  AICc = 78.743887
##  free parameters = 2
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 1.00
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

fitEB1<-ratebytree(obj,list(x3),model="EB")
fitContinuous(t1,x3,model="EB")

## GEIGER-fitted comparative model of continuous data
##  fitted 'EB' model parameters:
##  a = -0.378183
##  sigsq = 3.254252
##  z0 = 0.403735
## 
##  model summary:
##  log-likelihood = -36.750942
##  AIC = 79.501884
##  AICc = 80.592793
##  free parameters = 3
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 0.77
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

and so on.

Secondly, the fitted multi-regime models should match the independent fits from fitContinuous in the geiger package - and the likelihood should be sum of the log-likelihoods from fitContinuous.

For instance:

fitBM2<-ratebytree(c(t2,t1),list(x2,x3))
fitC.bm1<-fitContinuous(t2,x2)
fitC.bm2<-fitContinuous(t1,x3)
fitC.bm1

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
##  sigsq = 0.799912
##  z0 = 1.913625
## 
##  model summary:
##  log-likelihood = -48.562175
##  AIC = 101.124350
##  AICc = 101.448675
##  free parameters = 2
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 1.00
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

fitC.bm2

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
##  sigsq = 1.059273
##  z0 = 0.540955
## 
##  model summary:
##  log-likelihood = -37.111074
##  AIC = 78.222148
##  AICc = 78.743887
##  free parameters = 2
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 1.00
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

logLik(fitC.bm1)+logLik(fitC.bm2)

## [1] -85.67325

fitBM2

## ML common-rate model:
##  s^2  a[1]   a[2]    k   logL
## value    0.9021  1.9136  0.541   3   -85.9892    
## 
## ML multi-rate model:
##   s^2[1] s^2[2]   a[1]   a[2]    k   logL
## value    0.7999  1.0593  1.9136  0.541   4   -85.6732    
## 
## Likelihood ratio: 0.6319 
## P-value (based on X^2): 0.4267 
## 
## R thinks it has found the ML solution.

ratebytree and fitContinuous are completely independent implementations of continuous character model-fitting in R, so the latter could provide a handy way to check the former, and vice versa.

That's all for now.

Pearson correlation with phylogenetic data

Today I was asked the following question:

“I am interested in testing the correlation between 2 traits taking into account for the phylogenetic structure in my dataset. PGLS is a good way to test it, but I was wondering if a simple correlation test (using a Pearson r) would be possible. Your phyl.vcv computes a phylogenetic trait variance-covariance matrix between two variables. Can I use this matrix to compute a phylogenetic Pearson r value? If so, can this r-value be used to compute the statistic to test the significance of the correlation (with n-2 df in a t distribution)?”

The question was also posted here.

The answer is basically, yes & yes - but it is not strictly necessary as PGLS will give us the same p-value. Let's see how:

library(phytools)
## our data
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

X

##              x           y           z
## A -0.723942909  0.96553388  0.50467165
## B -0.034600412  1.48143249  0.66533528
## C -0.001282309  0.62094719  0.60352347
## D -0.100971061  1.23102030  0.01866007
## E -0.363619319  0.06831715  0.91782821
## F -0.400834652 -0.07312087  0.77584466
## G -0.037785409  0.18406320  1.18509956
## H -1.284112703  0.26824899  0.15203366
## I -0.573447373  1.56312743  0.58464867
## J -0.275070942  0.05165640 -0.24120290
## K -1.050699705  0.31715393 -0.63071969
## L -0.689968116  0.22713074 -0.04087866
## M -0.572979102  0.16002629 -0.67389194
## N -0.353661739  0.31157683 -0.60997825
## O -1.720295244  0.57079006 -2.24647340
## P -2.830254149  0.78574612 -2.39094415
## Q -0.104129552  0.93488329  1.66983162
## R -0.324883255  1.00993141 -1.87523279
## S -1.330114250  0.53192888 -2.47716184
## T -0.394635253  0.84544343 -2.05151552
## U -2.788040191  0.75579827 -1.90934711
## V  0.076717481 -1.13500947 -0.01313445
## W  1.388265725 -1.15038331  0.50393587
## X  1.423667629 -1.41685533  1.29662644
## Y  1.459867791 -0.53755989  1.83210421
## Z  2.235835060 -1.28003540  0.84061852

First, let's compute the evolutionary (phylogenetic) correlation between x & y:

## our covariance matrix
obj<-phyl.vcv(X,vcv(tree),1)
obj$R

##            x          y         z
## x  0.6810667 -0.0592844 0.3212055
## y -0.0592844  0.4238996 0.0413417
## z  0.3212055  0.0413417 0.6319725

## correlation between x & y
r.xy<-cov2cor(obj$R)["x","y"]
## t-statistic & P-value
t.xy<-r.xy*sqrt((Ntip(tree)-2)/(1-r.xy^2))
P.xy<-2*pt(abs(t.xy),df=Ntip(tree)-2,lower.tail=F)
P.xy

## [1] 0.5915613

As noted before, this is the same P-value as we would obtain fitting a linear model using PGLS of y~x or x~y.

fit.yx<-gls(y~x,data=as.data.frame(X),correlation=corBrownian(1,tree))
anova(fit.yx)

## Denom. DF: 24 
##             numDF   F-value p-value
## (Intercept)     1 0.2645170  0.6117
## x               1 0.2957732  0.5916

fit.xy<-gls(x~y,data=as.data.frame(X),correlation=corBrownian(1,tree))
anova(fit.xy)

## Denom. DF: 24 
##             numDF   F-value p-value
## (Intercept)     1 0.1812473  0.6741
## y               1 0.2957732  0.5916

Now let's see an example in which a weak correlation exists between the two characters:

## correlation between x & z
r.xz<-cov2cor(obj$R)["x","z"]
## t-statistic & P-value
t.xz<-r.xz*sqrt((Ntip(tree)-2)/(1-r.xz^2))
P.xz<-2*pt(abs(t.xz),df=Ntip(tree)-2,lower.tail=F)
P.xz

## [1] 0.01112811

fit.zx<-gls(z~x,data=as.data.frame(X),correlation=corBrownian(1,tree))
anova(fit.zx)

## Denom. DF: 24 
##             numDF  F-value p-value
## (Intercept)     1 0.003181  0.9555
## x               1 7.566719  0.0111

## visualize it
phylomorphospace(tree,X[,c("x","z")],node.size=c(0,0))
points(X[,c("x","z")],pch=21,cex=1.5,bg="grey")
text(-3.2,2,paste("evolutionary correlation\nr = ",round(r.xz,4),
    ", P = ",round(P.xz,4),sep="" ),pos=4)

Neat.

The data for this example were simulated as follows:

tree<-pbtree(n=26,tip.label=LETTERS)
X<-sim.corrs(tree,
    vcv=matrix(c(1.0,0.0,0.3,
            0.0,1.0,0.0,
            0.3,0.0,1.0),3,3))
colnames(X)<-c("x","y","z")