cospeciation method applied to empirical dataset

Recently I gave a simple cospeciation test which, given two input trees & a matrix of associations between their respective tips, compares topological tree distance with tree distances between random trees - either by permuting the labels of each tree, or by simulating trees under a Yule process (the default). The null hypothesis of the test is thus no similarity between the trees; and if the null hypothesis is rejected, this suggests that the association of the trees is greater than expected by chance. Note that this is not the same as Page's (1994) method of counting cospeciation events. I may look into that subsequently. I also didn't invent this method (it seems to have been used in a number of empirically studies) - though nor have I yet identified the original reference.

I thought it might be interesting to apply this method to the data of Lopez-Vaamonde et al. (2001) on the co-evolution of pollinators & parasites in fig wasps. (I used these data in an earlier example of co-phylogenetic plotting.)

library(phytools)
Pleistodontes<-read.tree(text="(((1,(2,3)100)50,(((((((4,(5,(6,7)100)69)100,(8,(9,10))81)61,11),(12,13)100),(14,15)100),(16,17)96),18)),19);")
Pleistodontes$tip.label<-c("P._froggatti",
    "P._blandus",
    "P._near_blandus",
    "P._proximus",
    "P._athysanus",
    "P._cuneatus",
    "P._astrabocheilus",
    "P._macrocainus",
    "P._rigisamos",
    "P._imperialis",
    "P._nigriventris",
    "P._xanthocephalus",
    "P._greenwoodi",
    "P._rieki",
    "P._plebejus",
    "P._schizodontes",
    "P._nitens",
    "P._spec._nov._1",
    "P._regalis")
## rotate all nodes to match orientation in article
Pleistodontes<-rotateNodes(Pleistodontes,"all")
Sycoscapter<-read.tree(text="(1,((((2,3)100,((4,5),(6,7)78)100),((8,9)50,10)99),(((11,12),13)91,(14,15)77)));")
Sycoscapter$tip.label<-c("S._australis_{macrophylla}",
    "S._3_{playpoda}",
    "S._5_{cerasicarpa}",
    "S._6_{lilliputiana}",
    "S._7_{subpuberula}",
    "S._1_{brachypoda}",
    "S._15_{rubiginosa}",
    "S._9_{triradiata}",
    "S._10_{crassipes}",
    "S._8_{pleurocarpa}",
    "S._11_{glandifera}",
    "S._12_{xylosycia}",
    "S._14_{hesperidiiformis}",
    "S._2_{aff._obliqua}",
    "S._4_{obliqua}")
Sycoscapter<-rotateNodes(Sycoscapter,"all")
assoc<-cbind(
    Pleistodontes$tip.label[c(7,8,5,6,18,1,3,4,10,15,13,16,12,14,19)],
    Sycoscapter$tip.label)
assoc

##       [,1]                [,2]                        
##  [1,] "P._greenwoodi"     "S._4_{obliqua}"            
##  [2,] "P._xanthocephalus" "S._2_{aff._obliqua}"       
##  [3,] "P._plebejus"       "S._14_{hesperidiiformis}"  
##  [4,] "P._rieki"          "S._12_{xylosycia}"         
##  [5,] "P._blandus"        "S._11_{glandifera}"        
##  [6,] "P._regalis"        "S._8_{pleurocarpa}"        
##  [7,] "P._nitens"         "S._10_{crassipes}"         
##  [8,] "P._schizodontes"   "S._9_{triradiata}"         
##  [9,] "P._imperialis"     "S._15_{rubiginosa}"        
## [10,] "P._athysanus"      "S._1_{brachypoda}"         
## [11,] "P._astrabocheilus" "S._7_{subpuberula}"        
## [12,] "P._proximus"       "S._6_{lilliputiana}"       
## [13,] "P._macrocainus"    "S._5_{cerasicarpa}"        
## [14,] "P._cuneatus"       "S._3_{playpoda}"           
## [15,] "P._froggatti"      "S._australis_{macrophylla}"

Here's our co-phylogenetic object:

obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc)

## Rotating nodes to optimize matching...
## Done.

obj

## Object of class "cophylo" containing:
## 
## (1) 2 (possibly rotated) phylogenetic trees in an object of class "multiPhylo".
## 
## (2) A table of associations between the tips of both trees.

and our cospeciation test:

test<-cospeciation(Pleistodontes,Sycoscapter,assoc=assoc,
    nsim=1000)
test

## 
## Co-speciation test based on RF distance.
## P-value obtained via simulation .
## 
##    RF distance: 18 
##    Mean(SD) from null: 22.7(1.5)
##    P-value: 0.026973

Finally, here are the results plotted together:

par(mfrow=c(2,1))
plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8,
    mar=c(0.1,0.1,3.1,0.1))
title(main="a) co-phylogenetic plot",cex.main=1.1)

par(mar=c(4.1,4.1,3.1,1.1))
plot(test) 
title(main=
    "b) topological dissimilarity between parasite & host phylogenies",
    cex.main=1.1)