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*_{2},
but perhaps not so much in trait *y*_{1}. 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)
```

Hi Liam! How can I plot the legend for the x factor? I tried with add.simmap.legend(colors=x) but it doesnt work...

ReplyDeletethank you,

isabel

add.simmap.legend(colors=x) will automatically go to the interactive mode, so you have to click where on the plot you would like to add the legend. You can also use legend in base R which is super flexible.

DeleteI might want to thank you for the endeavors you have made in composing this article. I am trusting the same best work from you later on too.. cheap postcard printing

ReplyDelete