#######################################################################################################################################
#
# R function to assess change in the Jaccard dissimilarity index based on the paper by Villéger & Brosse (Ecological Indicators)
#      for any question, contact Sébastien Villéger: sebastien.villeger@cnrs.fr
#
#
# INPUTS: two matrices (C x S) named "historical" and "current", 
#         with for the C communities of interest Presence/Absence (coded as 0/1) for the same pool of species S at two periods
#   - row and column names of "historical" and "current" have to be identical
#   - NA are replaced by "0"
#   - all communities have to contain at least one species for each of the two periods
#
#
# OUTPUTS: a list of 4 objects:
#
#   - "richness": a matrix (C x 2) with species richness in the C communites (rows) for the two periods considered (columns)
#
#   - "nmpairsC": a matrix (C*(C-1)/2 , 3) with codes for pairs of communities in the first column and corresponding communities' names in columns 2 & 3
#
#   - "dissJ": a matrix(C*(C-1)/2 , 3) with for each pair of communities: 
#               Jaccard's dissimilarity for historical ("Jacc_Hist") and current ("Jacc_Curr") situation
#               and the associated temporal change ("Jacc_delta"=Jacc_Curr-Jacc_Hist)
#             
#   - "diss_comp": a matrix(C*(C-1)/2 , 6) with for each pair of communities, 
#               the turnover ("Turn") and nestedness ("Nest") components for historical ("Hist") and current ("Curr") situations
#               and the associated changes (delta=Curr-Hist)
#
#   - "details": a list containing 2 matrices, which rows are the (C*(C-1)/2) pairwise combinations between communities 
#               - "abcefg": the number of species historically shared or not and the associated changes from historical to current situation 
#               - "intext": the number of species introduced or extirpated within each modality of introduction or extirpation
#                 (see figure 1 for notations)
#
# EXAMPLES are provided at the end of the script (including study cases shown on figure 3)
#
#######################################################################################################################################

chgedissJ<-function (historical , current )  {

##################################################
# data checking
mhist<-as.matrix(historical) ; mhist[which(is.na(mhist))]<-0
mcurr<-as.matrix(current) ; mcurr[which(is.na(mcurr))]<-0    
if (sum(row.names(mhist)!=row.names(mcurr))!=0) stop(" 'historical' and 'current' must have species composition for the same set of communities")
if (sum(colnames(mhist)!=colnames(mcurr))!=0) stop(" 'historical' and 'current' composition must be evaluated on the same set of species")
C<-nrow(mhist) ; nmC<-row.names(mhist)
S<-ncol(mhist) ; nmS<-colnames(mhist)

##################################################
# number of species
richness<-matrix(0,C,2,dimnames=list(nmC, c("Historical","Current")))

richness[,1]<-apply(mhist,1,sum)
richness[,2]<-apply(mcurr,1,sum)
if (length(which(richness[,1]==0))!=0) stop(paste("error :", nmC[which(richness[,1]==0)] , "had historically no species"),sep="")
if (length(which(richness[,2]==0))!=0) stop(paste("error :", nmC[which(richness[,2]==0)] , "has currently no species"),sep="")

##################################################
# names of communities pairs
nmpairsC<-matrix(0,C*(C-1)/2,3)
flag<-1
for (i in 1:(C-1))
for (j in (i+1):C)
{
nmpairsC[flag,]<-c(paste(nmC[i],nmC[j],sep="-"),nmC[i],nmC[j])
flag<-flag+1
} # end fo i,j
colnames(nmpairsC)<-c("commI-commII","commI","commII")

##################################################
# matrix with number of species unique or shared by each pair of communities and associated changes from historical to current situation
vardetails<-c("i","j","k","m","n","t","u","v","w","x","y","z")

abcefg<-matrix(0,nrow(nmpairsC),6, dimnames=list(nmpairsC[,1],c("a","b","c","e","f","g")) )

intext<-matrix(0,nrow(nmpairsC),length(vardetails), dimnames=list(nmpairsC[,1],vardetails) )

flag<-1

for (i in 1:(C-1))
for (j in (i+1):C)
{
abcefg[flag,"a"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==1) ) # initially present in comm i & comm j
abcefg[flag,"b"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==0) ) # initially present only in comm i
abcefg[flag,"c"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==1) ) # initially present only in comm j

intext[flag,"x"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==1) ) # introduced in comm i & comm j
intext[flag,"y"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==0) ) # introduced only in comm i
intext[flag,"z"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==1) ) # introduced only in comm j
intext[flag,"v"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==1) ) # translocated from comm i to community j
intext[flag,"w"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==1) ) # translocated from comm j to community i

intext[flag,"i"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==0) ) # extirpated from comm i & comm j
intext[flag,"j"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==1) ) # extirpated only from comm i, not from comm j
intext[flag,"k"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==0) ) # extirpated only from comm j, not from comm i

intext[flag,"m"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==0) )  # extirpated from comm i, never present in comm j
intext[flag,"n"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==0) )  # extirpated from comm j, never present in comm i

intext[flag,"t"]<-length( which(mhist[nmC[i],]==1 & mhist[nmC[j],]==0 & mcurr[nmC[i],]==0 & mcurr[nmC[j],]==1) )  # extirpated from comm i and introduced in comunity j
intext[flag,"u"]<-length( which(mhist[nmC[i],]==0 & mhist[nmC[j],]==1 & mcurr[nmC[i],]==1 & mcurr[nmC[j],]==0) )  # extirpated from comm j and introduced in comunity i

flag<-flag+1
} # end fo i,j

abcefg[,"e"]<-intext[,"x"]-intext[,"i"]-intext[,"j"]-intext[,"k"]+intext[,"v"]+intext[,"w"] # change in number of species shared by comm 1 and comm 2
abcefg[,"f"]<-intext[,"y"]+intext[,"j"]+intext[,"u"]-intext[,"m"]-intext[,"v"]-intext[,"t"] # change in number of species present only in comm 1
abcefg[,"g"]<-intext[,"z"]+intext[,"k"]+intext[,"t"]-intext[,"n"]-intext[,"w"]-intext[,"u"] # change in number of species present only in comm 2

##################################################
# computing Jaccard dissimlarity index

dissJ<-matrix(0,nrow(nmpairsC),3, dimnames=list(nmpairsC[,1],c("Jacc_Hist","Jacc_Curr","Jacc_delta")) )

# for each pair : total species richnesss and minimum number of species unique to one community
if (nrow(abcefg)==1) { 
        sumabc<-sum(abcefg[,c("a","b","c")])  ; sumabcefg<-sum(abcefg[,c("a","b","c","e","f","g")])       
        minbc<-min(abcefg[,c("b","c")])  ; minbfcg<-min(  (abcefg[,"b"]+abcefg[,"f"]), (abcefg[,"c"]+abcefg[,"g"]) )     } # end of if only one pair


if (nrow(abcefg)!=1) { 
        sumabc<-apply(abcefg[,c("a","b","c")],1,sum) ; sumabcefg<-apply(abcefg[,c("a","b","c","e","f","g")],1,sum) 
        minbc<-apply(abcefg[,c("b","c")],1,min)  ; minbfcg<-apply(  cbind(abcefg[,"b"]+abcefg[,"f"], abcefg[,"c"]+abcefg[,"g"]),1,min )     } # end of if more than 2 pairs


dissJ[,"Jacc_Hist"]<-(abcefg[,"b"]+abcefg[,"c"])/sumabc   # historical Jacccard dissimilarity
dissJ[,"Jacc_Curr"]<-(abcefg[,"b"]+abcefg[,"f"]+abcefg[,"c"]+abcefg[,"g"]) /sumabcefg   # current Jaccard dissimilarity
dissJ[,"Jacc_delta"]<-dissJ[,"Jacc_Curr"]-dissJ[,"Jacc_Hist"]   # difference in dissimilarity between historical and current situation

dissJ<-round(dissJ,3) # dissimilarity values rounded to 0.1%


##################################################
# components of dissimilarity: turnover and nestedness
diss_comp<-matrix(0,nrow(nmpairsC),6,dimnames=list(nmpairsC[,1], c("Turn_Hist","Turn_Curr","Turn_delta","Nest_Hist","Nest_Curr","Nest_delta") ) )

diss_comp[,"Turn_Hist"]<-2*minbc/(abcefg[,"a"]+2*minbc)  # taxonomic turnover for historical situation
diss_comp[,"Turn_Curr"]<-2*minbfcg/(abcefg[,"a"]+abcefg[,"e"]+2*minbfcg) # taxonomic turnover for current situation
diss_comp[,"Turn_delta"]<-diss_comp[,"Turn_Curr"]-diss_comp[,"Turn_Hist"] # change in taxonomic turnover

diss_comp[,"Nest_Hist"]<-abs(abcefg[,"b"]-abcefg[,"c"])*abcefg[,"a"]  / ( sumabc * (abcefg[,"a"]+2*minbc) ) # taxonomic nestedness for historical situation
diss_comp[,"Nest_Curr"]<-abs(abcefg[,"b"]+abcefg[,"f"]-(abcefg[,"c"]+abcefg[,"g"])) * (abcefg[,"a"]+abcefg[,"e"]) / ( sumabcefg * (abcefg[,"a"]+abcefg[,"e"]+2*minbfcg) )   # taxonomic nestedness for current situation
diss_comp[,"Nest_delta"]<-diss_comp[,"Nest_Curr"]-diss_comp[,"Nest_Hist"] # change in taxonomic nestedness

diss_comp<-round(diss_comp,3) # values rounded to 0.1%

##################################################
# listing results
res<-list( richness=richness,nmpairsC=nmpairsC,dissJ=dissJ,diss_comp=diss_comp,details=list(abcefg=abcefg,intext=intext) )
invisible(res)

} # end of function chgedissJ

#########################################################################################################################################
####################################################     END OF FUNCTION     ############################################################
#########################################################################################################################################



#########################################################################################################################################
# TWO EXAMPLES

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

if(show_example==T) {

####################################################
# basic example with 5 communities
ex_historical<-matrix(rbind(c(1,1,1,1,0,0,0,1,1,1,1,0,0,0,0),c(0,0,0,0,1,1,1,1,1,1,1,0,0,0,0),c(1,0,0,0,1,1,0,0,0,1,1,1,0,NA,1),
                      sample(c(0,1),15,T,prob=c(0.3,0.7)),sample(c(0,1),15,T,prob=c(0.5,0.5)) ),5,15,
                      dimnames=list(c("commA","commB","commC","commD","commE"),paste("sp",1:15,sep="")))
ex_current<-matrix(rbind(c(0,1,1,1,0,0,0,1,0,0,1,1,1,1,0),c(0,0,0,0,1,1,0,0,0,1,1,1,0,0,1),c(1,NA,0,0,1,1,0,0,0,1,1,1,1,1,1),
                      sample(c(0,1),15,T,prob=c(0.25,0.75)),sample(c(0,1),15,T,prob=c(0.3,0.7)) ),5,15,
                      dimnames=list(c("commA","commB","commC","commD","commE"),paste("sp",1:15,sep="")))
ex<-chgedissJ(ex_historical,ex_current)

print("Example with 5 communities: results for the 10 pairs of communities") 
print(ex)  

####################################################
# examples presented on figure 3

# same historical situation for the 8 study cases (a=6,b=5,c=4)
ex_hist<-matrix(0,2,21,T,dimnames=list(c("comm_I","comm_II"),paste("sp",1:21,sep="")))  
ex_hist[1,1:11]<-1 ;  ex_hist[2,6:15]<-1   ; ex_hist   
vard<-c("Jacc_delta","Turn_delta","Nest_delta")
 
# top left : e=1,f=0,g=1 
ex_curr<-ex_hist ; ex_curr[1,c(16)]<-1 ; ex_curr[2,c(16:17)]<-1 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 top left") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# top center : e=-2,f=-2,g=-1 
ex_curr<-ex_hist ; ex_curr[1,c(1:2,6:7)]<-0 ; ex_curr[2,c(6:7,12)]<-0 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 top center") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# top right : e=0,f=1,g=2 
ex_curr<-ex_hist ; ex_curr[1,c(19)]<-1 ; ex_curr[2,c(20:21)]<-1 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 top right") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# center left : e=2,f=-1,g=0 
ex_curr<-ex_hist ; ex_curr[1,c(16:17)]<-1 ; ex_curr[1,c(1)]<-0 ; ex_curr[2,c(16:17)]<-1 ;  ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 center left") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# center right : e=0,f=3,g=1 
ex_curr<-ex_hist ; ex_curr[1,c(16:18)]<-1 ; ex_curr[2,c(19)]<-1 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 center right") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# bottom left : e=1,f=1,g=-2 
ex_curr<-ex_hist ; ex_curr[1,c(16:17)]<-1 ; ex_curr[2,c(16)]<-1 ; ex_curr[2,c(12:13)]<-0 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 bottom left") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# bottom center : e=0,f=-2,g=2 
ex_curr<-ex_hist ; ex_curr[1,c(1:2)]<-0 ; ex_curr[2,c(16:17)]<-1 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 bottom center") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")

# bottom right : e=-2,f=1,g=-2 
ex_curr<-ex_hist ; ex_curr[1,c(16)]<-1 ; ex_curr[1,c(6:7)]<-0 ; ex_curr[2,c(6:7,12:13)]<-0 ; ex_chge<-chgedissJ(ex_hist,ex_curr)  
print("Fig. 3 bottom right") ; print(ex_chge$details$abcefg) ; cat("\n"); print(round(100*cbind(ex_chge$dissJ,ex_chge$diss_comp)[,vard],0)) ; cat("\n")


          
####################################################
} # end of if example
          
################################################################    END OF SCRIPT    ###################################################################