########################################################################################################################################
# betaQmult: R function to compute functional beta-diversity based on the multiplicative decomposition of the Rao's quadratic entropy  # 
# Four contrasted examples are provided at the end of the script
# Author: SÈbastien VillÈger, sebastien.villeger@cnrs.fr
#
# INPUTS:
#   - "abundances": matrix (C x S) of abundances of the S species in the C communities
#   - "functdist": matrix (S x S) or dist object with pairwise functional distances between the S species
#   - species names in "abundances" and in "functdist" must be identical
#      NA are automatically replaced by 0 in 'abundances'
#      NA are not allowed in 'functdist'
#
# OUTPUTS:  a list of 7 elements
#   - "nbc" : number of communities
#   - "nbsp" (vector, length C) : number of species in each community
#   - "weight" (vector, length C) : contribution of each community to the regional abundances
#   - "eqmQalpha" : equivalent number of species of the abundance-weighted mean of local Rao's quadratic entropy values(Qalpha)
#   - "eqQgamma" : equivalent number of species of the regional Rao's quadratic entropy 
#   - "Qbeta" : eqQgamma/eqmQalpha, raw value of functional beta-diversity
#   - "Qbetast" : standardized functional beta-diversity, =(Qbeta-1)/(C-1)
#####################################################################################################################################

betaQmult<-function(functdist,abundances) {

# checking NA absence in 'functdist'
if (length(which(is.na(functdist)==T))!=0)  stop("error : NA in 'functdist'")

# if necessary, conversion of 'functdist' into a matrix format
if (is.matrix(functdist)==F) functdist<-as.matrix(functdist)

# standardization of functional distances
stfunctdist<-functdist/max(functdist)

# checking 'abundances' and 'functdist' dimensions
if (dim(functdist)[1]!=dim(functdist)[2])  stop("error : check 'functdist' dimensions")
if (dim(abundances)[2]!=dim(functdist)[1]) stop("error : different number of species in 'abundances' and 'functdist' ")
if (sum(colnames(abundances)!=row.names(functdist))!=0) stop(" species names in 'functdist' and 'abundances' must be identical")

# replacement of NA by 0 in abundances
abundances[which(is.na(abundances))]<-0  

# number and names of communities
c<-dim(abundances)[1] ; namescomm<-row.names(abundances)
abundances<-as.matrix(abundances)

#  species richness and total abundances in each community
abloc<-apply(abundances,1,sum)
nbsp<-apply(abundances,1,function(x) {length(which(x>0))} )  
names(nbsp)<-row.names(abundances)

# species relative abundances in each community
locabrel<-abundances/abloc

# contribution of each community to the regional abundances 
weight<-abloc/sum(abloc) ; names(weight)<-row.names(abundances)

# species relative abundances at regional level
totabrel<-as.matrix(weight %*% locabrel)

# Rao's quadratic entropy of each community
Qalpha<-rep(0,c) ; names(Qalpha)<-namescomm 
for (i in 1:c)
{Qalpha[i]<-locabrel[i,] %*% stfunctdist %*% locabrel[i,] } # end of c

# alpha diversity = abundance-weighted mean of Qalpha and equivalent number of species
mQalpha<-( Qalpha%*%weight )[1]
eqmQalpha<-( 1/(1-mQalpha ))

# gamma diversity = Rao's quadratic entropy of the pool of community and equivalent number of species
Qgamma<-( totabrel %*% stfunctdist %*% t(totabrel) ) [1]
eqQgamma<-(1/(1-Qgamma))

# raw functional beta diversity
Qbeta<-eqQgamma/eqmQalpha

# standardized functional beta diversity given number of communities
Qbetast<-round( (Qbeta-1)/(c-1) ,6)

# list of results
resQ<-list(nbc=c, nbsp=nbsp, weight=weight, eqmQalpha=eqmQalpha, eqQgamma=eqQgamma, Qbeta=Qbeta, Qbetast=Qbetast )

invisible(resQ)

} # end of function betaQmult


#######################################################################################################################
#######################################################################################################################
#######################################################################################################################
# BASIC EXAMPLES WITH ONLY TWO COMMUNITIES

show_example<-F   # change value to 'TRUE' to see examples

if(show_example==T) {

# species traits and position in the 2D functional space  (a~c) != (e~f) != (d~g) != b  ['!=' means 'different']
tr<-matrix(c(0,1,4,1,0.25,1,2.5,4,2,0,2,0.25,2.5,3.75) ,7,2,T)   ; row.names(tr)<-c("a","b","c","d","e","f","g")
plot(tr,type="n",xlab="Trait 1", ylab="Trait2") ; text(tr[,1],tr[,2], row.names(tr))
# functional distance between species (Euclidean distance)
dfunct<-dist(tr)


# same species composition and minimal functional beta-diversity 
ab1<-matrix(1,2,7) ; row.names(ab1)<-c("I","II") ; colnames(ab1)<-row.names(tr)
ex1<-betaQmult(dfunct,ab1)
ex1


# same species composition but high functional beta diversity
ab2<-rbind(c(100,1,1,1,1,10,1) , c(1,100,1,10,1,1,1)) ; row.names(ab2)<-c("I","II") ; colnames(ab2)<-row.names(tr)
ex2<-betaQmult(dfunct,ab2) 
ex2


# no species in common but low functional beta diversity
ab3<-rbind(c(0,0,0,10,10,0,0) , c(0,0,0,0,0,10,10) ) ; row.names(ab3)<-c("I","II") ; colnames(ab3)<-row.names(tr)
ex3<-betaQmult(dfunct,ab3) 
ex3


# no species in common and high functional beta diversity 
# i.e. co-occuring species functionally similar but strong functional dissimilarity between species from different communities
ab4<-rbind(c(0,0,0,10,0,0,10) , c(0,0,0,0,10,10,0)) ; row.names(ab4)<-c("I","II") ; colnames(ab4)<-row.names(tr)
ex4<-betaQmult(dfunct,ab4) 
ex4


} # end of show example
#######################################################################################################################