#####script to caculate several community diverisy metrics from phylogenetic, trait and community composition datasets.
####written by Marc W Cadotte
####July 2008

#Requires ape package
library(ape)

#phylogeny file (newick format), use read.nexus for nexus format file
my.tree<-read.tree("file name")#get tree

#this is file with 2 columns with NO headers. Column 1 is species names and col 2 is community names (e.g., 1, 2... or A, B ..., etc.)
sp.list<-read.table("file name",header=F)

#this is a trait matrix, where col 1 is species names and col 2...n are various continuous traits (or 1 & 0 for discrete states). This file needs headers. Make sure species names in all three files are identical, use the "match" command to find errors.
trait.data<-read.table(file="file name",header=T)

#make tree ultrametric
ultra.tree<-chronogram(my.tree)

species<-as.list(ultra.tree$tip.label)#list of species in tree

traits<-scale(trait.data[,2:ncol(trait.data)]) #standardize trait data to mean = 0 and var = 1

trait.clust<-hclust(dist(traits),method="average") #hierarchical clustering as per Petchey & Gaston

trait.clust$labels<-species[trait.clust$order] #label cluster objects as species names -plot by "plot(trait.clust)"

clust.phy <- as.phylo(trait.clust) #transform trait.clust into phylo object

comID<-unique(sp.list[,2]) #community names

#for output
PD<-numeric(length(comID)) #phylogenetic diversity
sp.rich<-numeric(length(comID)) #number of species
FD<-numeric(length(comID)) #Functional diversity (Petchey & Gaston 2002)
FAD<-numeric(length(comID)) # Functional attribute diversity (Walker et al. 1999)

SD.D<-data.frame()#dataframe for single variable standard deviations

for (i in 1:length(comID)) {#loop through community IDs & calculate PD, FD and FAD

	ID<-comID[i]#select community.ID
	myclade<-sp.list[sp.list[,2]==ID,]#get list of species in that  community
	dropme<-species[!species %in% myclade[,1]]#species in tree but not in community
	sub.tree<-drop.tip(ultra.tree,dropme) #create phylogeny for individual community
	PD[i]<-sum(sub.tree$edge.length) #calculate phylogenetic diversity
	sp.rich[i]<-length(sub.tree$tip.label) #number of species

	sub.clust<-drop.tip(clust.phy,dropme) #subset of trait dendrogram
	FD[i] <- sum(sub.clust$edge.length) #calculate functional diversity

	fad.spp<-match(myclade[,1],trait.data$species)
	sub.traits<-traits[fad.spp,]
	FAD[i] <-sum(dist(sub.traits))
	
	SD<-matrix(nrow=1, ncol=ncol(traits)) #temporarily hold sd values

	for (j in 1:ncol(traits))  {#loop to calculate sd values for each trait in community i
						
		SD[,j] <- sd(sub.traits[,j]) #calculate sd and add it to temp holder
					
	} #end j loop	
			
		
		SD.D<-data.frame(rbind(SD.D,as.data.frame(SD))) #add temp holder to perm dataframe


}#end i loop
	
colnames(SD.D)<-colnames(traits)
sp.rich<-log(sp.rich)


#OUTPUT RESULTS
results <- data.frame(comID,PD,sp.rich,FD,FAD)
results <- data.frame(cbind(results,SD.D))

write.csv(results, file="output file name")