######### ASSESSING CONSEQUENCES OF TRAIT CORRELATIONS FOR FUNCTIONAL DIVERSITY INDICES ######### 

#Author: Jon Lefcheck
#Last updated: 2014-07-23
#Contact: jslefche@vims.edu

#Load required libraries
library(MASS)
library(FD)
library(ggplot2)
library(plyr)

#Generate correlated traits for 2-10 traits from 10-100 species (increments of 10)
cormats.list=lapply(seq(10,100,10),function(S) { #Number of species
  lapply(2:9,function(traits) { #Number of traits
    lapply((1:99)/100,function(rho) { #Degree of correlation among traits
      Sigma=matrix(rho,nrow=traits,ncol=traits)
      diag(Sigma)=1
      mvrnorm(S,mu=rep(1,nrow(Sigma)),Sigma=Sigma,empirical=T)
    } ) 
  } )
} )

#Confirm the average pairwise correlation for each matrix is rho
all(laply(seq_along(cormats.list),function(S) {
  all(laply(seq_along(cormats.list[[S]]),function(traits) {
    all(laply(seq_along(cormats.list[[S]][[traits]]),function(rho) {
      round(mean(cor(cormats.list[[S]][[traits]][[rho]])[lower.tri(cor(cormats.list[[S]][[traits]][[rho]]))]),2)==rho/100
    } ) ) 
  } ) )
} ) )

#Calculate FD for each level of correlation, species richness, and trait number combination
fd=ldply(seq_along(cormats.list),.progress="text",function(S) { #Number of species
  ldply(seq_along(cormats.list[[S]]),function(traits) { #Number of traits
    ldply(seq_along(cormats.list[[S]][[traits]]),function(rho) { #Correlation
      dat=cormats.list[[S]][[traits]][[rho]]
      rownames(dat)=paste("sp",1:nrow(dat))
      tr=dbFD(dat,messages=F)
      rel.mat=matrix(rep(1,nrow(dat))/nrow(dat),nrow=1)
      RaoQ=apply(rel.mat,1,function(x) t(x) %*% as.matrix(gowdis(dat)) %*% x)
      data.frame(
        no.species=S*10,
        no.traits=traits+1,
        correlation=rho/100,
        index=c("FRic","FEve","FDiv","FDis","RaoQ.dbFD","RaoQ"),
        value=c(tr$FRic,tr$FEve,tr$FDiv,tr$FDis,tr$RaoQ,RaoQ) )
    } )
  } )
} )
  
#Scale FRic by maximum for each level of number of species
fd[fd$index=="FRic",]=ddply(fd[fd$index=="FRic",],c("no.species","no.traits"),function(x) cbind(x[,-ncol(x)],value=x$value/max(x$value)))

#Plot results
ggplot(subset(fd,index!="RaoQ.dbFD"),aes(x=correlation,y=value,group=no.traits,col=no.traits))+#,shape=as.factor(no.species))+
  #geom_point(size=4,alpha=0.1)+
  #geom_jitter(position=position_jitter(width=0.05),alpha=0.7)+
  geom_line(lwd=1,alpha=0.5)+
  #scale_shape_manual(values=c(0,1,2,3,4,5,6,7,8,9),guide=guide_legend(ncol=2),name="Richness")+
  scale_colour_gradient(high="black",low="grey50",name="Number of\ntraits")+
  scale_x_continuous(breaks=c(0,0.5,1),labels=c("0","0.5","1"))+
  facet_grid(index~no.species,scales="free")+
  labs(x="\nCorrelation",y="")+
  theme_bw(base_size=18)+
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),
        strip.text.x=element_text(size=12))

#Generate Figure 2 for the main text
levels(fd$index)=c("Functional dispersion","Functional divergence","Functional evenness","Functional richness","Rao's Q","Rao's Q (dbFD)")

ggplot(subset(fd,index!="Rao's Q (dbFD)" & no.species==100),aes(x=correlation,y=value,group=no.traits,col=no.traits))+
  geom_line(lwd=1.2,alpha=0.5)+
  scale_colour_gradient(high="black",low="grey50",name="Number\nof traits")+
  facet_wrap(~index,scales="free")+
  scale_x_continuous(breaks=c(0,0.25,0.5,0.75,1),labels=c("0","0.25","0.5","0.75","1"))+
  labs(x="\nCorrelation",y="")+
  theme_bw(base_size=18)+
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),
        strip.text.x=element_text(size=12),legend.position=c(0.75,0.25))