print method for phylANOVA (the phylogenetic ANOVA)

A simple method that has long been fairly popular in the phytools package is the function phylANOVA that implements Ted Garland et al.'s (1993) simulation-based phylogenetic ANOVA method.

I just pushed a small update that adds a nice print method to the function. Here is an example.

library(phytools)

First, our data:

phenogram(tree,y1,ftype="off",col=make.transparent("blue",0.5))
tiplabels(pie=to.matrix(x,letters[1:3]),
    piecol=colorRampPalette(c("blue", "yellow"))(3),
    cex=0.4)

phenogram(tree,y2,ftype="off",col=make.transparent("blue",0.5))
tiplabels(pie=to.matrix(x,letters[1:3]),
    piecol=colorRampPalette(c("blue", "yellow"))(3),
    cex=0.4)

It looks like there might be an interesting pattern in y₂, but perhaps not so much in trait y₁. In fact, that is precisely what phylANOVA will help us quantify:

phylANOVA(tree,x,y1)

## ANOVA table: Phylogenetic ANOVA
## 
## Response: y
##              Sum Sq  Mean Sq  F value Pr(>F)
## x          9.686996 4.843498 2.035856   0.19
## Residual 230.772367 2.379097                
## 
## P-value based on simulation.
## ---------
## 
## Pairwise posthoc test using method = "holm"
## 
## Pairwise t-values:
##           a         b        c
## a  0.000000  1.007945 2.015685
## b -1.007945  0.000000 0.869513
## c -2.015685 -0.869513 0.000000
## 
## Pairwise corrected P-values:
##       a     b     c
## a 1.000 0.798 0.186
## b 0.798 1.000 0.798
## c 0.186 0.798 1.000
## ---------

phylANOVA(tree,x,y2)

## ANOVA table: Phylogenetic ANOVA
## 
## Response: y
##            Sum Sq   Mean Sq   F value Pr(>F)
## x        79.35828 39.679140 42.007471  0.001
## Residual 91.62362  0.944573                 
## 
## P-value based on simulation.
## ---------
## 
## Pairwise posthoc test using method = "holm"
## 
## Pairwise t-values:
##          a         b         c
## a 0.000000 -3.225067 -8.964991
## b 3.225067  0.000000 -5.243302
## c 8.964991  5.243302  0.000000
## 
## Pairwise corrected P-values:
##       a     b     c
## a 1.000 0.005 0.003
## b 0.005 1.000 0.003
## c 0.003 0.003 1.000
## ---------

Neat. The data were simulated as follows:

tree<-pbtree(n=100)
Q<-matrix(c(-2,1,1,
    1,-2,1,
    1,1,-2),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x<-as.factor(sim.history(tree,Q)$states)
y1<-fastBM(tree)
y2<-fastBM(tree,sig2=0.5)+as.numeric(x)