find.pairs<-function(groups,morph.dist,pat.dist) { taxes.1<-vector(length=0) taxes.2<-vector(length=0) group.1<-vector(length=0) group.2<-vector(length=0) pair.morph<-vector(length=0) pair.pat<-vector(length=0) for(i in 1:nrow(morph.dist)) { for(j in i:ncol(morph.dist)) { if(i!=j) { tax.1<-rownames(morph.dist)[i] tax.2<-colnames(morph.dist)[j] taxes.1<-c(taxes.1,tax.1) taxes.2<-c(taxes.2,tax.2) group.1<-c(group.1,groups[which(groups[,1]==tax.1),2]) group.2<-c(group.2,groups[which(groups[,1]==tax.2),2]) pair.morph<-c(pair.morph,morph.dist[tax.1,tax.2]) pair.pat<-c(pair.pat,pat.dist[tax.1,tax.2]) } } } results<-data.frame(taxes.1,taxes.2,group.1,group.2,pair.morph,pair.pat) colnames(results)<-c("Taxon 1","Taxon 2","Group 1","Group 2","Morphological distance","Patristic distance") return(results) } plot.groups<-function(data,col.point,col.line,lwd.shadow=NULL,...) { groups<-unique(data$"Group 1") group.pairs<-vector(length=length(groups),mode="list") for(i in 1:nrow(data)) { if(data[i,"Group 1"] == data[i,"Group 2"]) { group.pairs[[which(groups==data[i,"Group 1"])]]<-c(group.pairs[[which(groups==data[i,"Group 1"])]],i) } } group.dists<-vector(length=length(groups),mode="list") for(i in 1:length(group.dists)) { group.dists[[i]]<-cbind(data[group.pairs[[i]],"Morphological distance"],data[group.pairs[[i]],"Patristic distance"]) } group.lin<-vector(length=length(groups),mode="list") for(i in 1:length(group.lin)) { group.lin[[i]]<-loess(V1~V2,as.data.frame(group.dists[[i]])) } plot(group.dists[[1]][,2],group.dists[[1]][,1],xlab="Patristic Distance", ylab="Morphological Dissimilarity",col=col.point[[1]], xlim=c(min(data$"Patristic distance"),max(data$"Patristic distance")), ylim=c(min(data$"Morphological distance"),max(data$"Morphological distance")),...) for(i in 2:length(groups)) { points(group.dists[[i]][,2],group.dists[[i]][,1],col=col.point[[i]],...) } if(is.null(lwd.shadow)==F) { for(i in 1:length(groups)) { lines(sort(group.dists[[i]][,2]),sort(predict(group.lin[[i]])),col="white",lwd=lwd.shadow) } } for(i in 1:length(groups)) { lines(sort(group.dists[[i]][,2]),sort(predict(group.lin[[i]])),col=col.line[[i]],...) } } char.sat<-function(tree,matrix,data,group,nsim,metric) { pat.dist<-as.matrix(distTips(tree,method = "patristic")) groups<-unique(data$"Group 1") group.pairs<-vector(length=length(groups),mode="list") for(i in 1:nrow(data)) { if(data[i,"Group 1"] == data[i,"Group 2"]) { group.pairs[[which(groups==data[i,"Group 1"])]]<-c(group.pairs[[which(groups==data[i,"Group 1"])]],i) } } group.dists<-vector(length=length(groups),mode="list") for(i in 1:length(group.dists)) { group.dists[[i]]<-cbind(data[group.pairs[[i]],"Morphological distance"],data[group.pairs[[i]],"Patristic distance"]) } sat.group<-vector(length=length(groups)) for(i in 1:length(sat.group)) { mm.group<-drm(V1~V2, data = as.data.frame(group.dists[[i]]), fct = MM.2()) sat.group[i]<-coef(mm.group)[1] } models<-vector(length=ncol(matrix[[2]][[3]]),mode="list") for(i in 1:length(models)) { char<-matrix[[2]][[3]][tree$tip.label,i] char[which(char=="")]<-NA test.tree<-drop.tip(tree,which(is.na(char))) char<-char[test.tree$tip.label] if(length(unique(char))>1) { models[[i]]<-fitMk(test.tree,char,model="ER") } else(models[[i]]<-NA) } characters<-vector(length=ncol(matrix[[2]][[3]]),mode="list") for(i in 1:length(models)) { if(length(models[[i]])>1) { sim.matrix<-matrix(nrow=nrow(models[[i]][[3]]),ncol=ncol(models[[i]][[3]]) ,data=models[[i]][[2]][models[[i]][[3]]]) rownames(sim.matrix)<-colnames(sim.matrix)<-1:nrow(sim.matrix) for (j in 1:nrow(sim.matrix)) { sim.matrix[j,j]<- -sum(sim.matrix[,j],na.rm=T) } characters[[i]]<-sim.Mk(tree,sim.matrix,nsim=nsim) } else(characters[[i]]<-NA) } sim.sat<-matrix(nrow=nsim,ncol=length(groups)) for(w in 1:nsim) { sim.matrix<-matrix sim.matrix[[2]][[3]]<-sim.matrix[[2]][[3]][tree$tip.label,] for(j in 1:ncol(sim.matrix[[2]][[3]])) { if(length(characters[[j]])>1) { sim.matrix[[2]][[3]][,j]<-characters[[j]][tree$tip.label,w] } else(sim.matrix[[2]][[3]][,j]<-NA) } sim.dist<-calculate_morphological_distances(sim.matrix,distance_metric=metric)$distance_matrix sim.dist<-sim.dist[rownames(pat.dist),colnames(pat.dist)] sim.pairs<-find.pairs(group,sim.dist,pat.dist) group.pairs<-vector(length=length(groups),mode="list") for(i in 1:nrow(sim.pairs)) { if(sim.pairs[i,"Group 1"] == sim.pairs[i,"Group 2"]) { group.pairs[[which(groups==sim.pairs[i,"Group 1"])]]<-c(group.pairs[[which(groups==sim.pairs[i,"Group 1"])]],i) } } group.dists<-vector(length=length(groups),mode="list") for(i in 1:length(group.dists)) { group.dists[[i]]<-cbind(sim.pairs[group.pairs[[i]],"Morphological distance"],sim.pairs[group.pairs[[i]],"Patristic distance"]) } for(i in 1:length(groups)) { mm.group<-drm(V1~V2, data = as.data.frame(group.dists[[i]]), fct = MM.2()) sim.sat[w,i]<-coef(mm.group)[1] } } names(sat.group)<-groups colnames(sim.sat)<-groups results<-list(sat.group,sim.sat) names(results)<-c("Observed","Simulated") return(results) }