FunD_Rcode.txt

## R code for obtaining functional diversity (FD) based on Chao et al. (2018) 
# Please cite Chao et al. (2018) An attribute-diversity approach to functional diversity, functional beta diversity, and related (dis)similarity measures. 
# Under revision, Ecological Monographs. 
# NOTE: The packages "cluster", "dplyr" and "ggplot2" must be installed and loaded before running the scripts. 
# These three packages can be downloaded from CRAN.
#
# input data:
# fs denotes Fushan (second census), Kt denotes Kenting (the third census) 
fskt <- read.table("Plant_Abundance.txt")
fskt = list(fs=fskt$Fushan, kt=fskt$Kenting)
dij_fskt = read.table("Plant_Distance_Matrix.txt")
dij_fskt = as.matrix(dij_fskt)

#------------------------------------------------------------------------------
# Functional Diversity of a single site for specified values of tau and q
#------------------------------------------------------------------------------
#' FD_MLE(data, dij, tau, q) is a function of obtaining FD index of order q.
#' @param data a vector of species sample frequencies.
#' @param dij a matrix of species-pairwise distances.
#' @param tau a numeric for a specified level of threshold distinctiveness.
#' @param q a numeric for a specified diversity order q. 
#' @return a numeric value of FD.
FD_MLE <- function(data, dij, tau, q){
  dij <- as.matrix(dij)
  dij[which(dij>tau,arr.ind = T)] <- tau
  a <- as.vector((1 - dij/tau) %*% data )  
  data <- data[a!=0]
  a <- a[a!=0]
  v <- data/a
  nplus <- sum(data)
  if(q==1){
    exp(sum(-v*a/nplus*log(a/nplus)))
  }else{
    (sum(v*(a/nplus)^q))^(1 / (1-q))
  }
}
FD_MLE(fskt$fs, dij_fskt, max(dij_fskt), 2)

#----------------------------------------------------------------
# Functional Diversity of N sites for various values of tau and q
#----------------------------------------------------------------
#' FunD(data, dij, tau, q, boot, datatype) is a function of obtaining FD tau-table or q-table.
#' @param data a list with N sites; each element of list is species abundances or species-by-sampling-unit incidence matrix.
#' @param dij a matrix of species-pairwise distances.
#' @param tau a numeric or a vector of levels of threshold distinctiveness.
#' @param q a numeric or a vector of diversity orders; the suggested range for q is [0, 3].
#' @param boot a numeric of number of bootstrap replications.
#' @param datatype data type of input data: individual-based abundance data (datatype = "abundance") or species-by-sampling-unit incidence matrix (datatype = "incidence_raw").
#' @return two matrices of FD; the first matrix is the q-profile information, the second matrix is the tau profile information.
FunD <- function(data, dij, tau, q, boot, datatype){
  EstiBootComm.Func = function(data, distance, datatype){
    if (datatype=="abundance") {
      n = sum(data)
    } else if (datatype=="incidence_raw") {
      n <- data[1]
      data <- data[-1]
      u=sum(data)
    }
    distance = as.matrix(distance)
    dij = distance[data!=0, data!=0]
    
    X = data[data>0]
    f1 <- sum(X == 1) ; f2 <- sum(X == 2) 	
    f0.hat <- ceiling(ifelse(f2>0, ((n-1)/n)*f1^2/2/f2, ((n-1)/n)*f1*(f1-1)/2))
    if (datatype=="abundance") {
      C1 = ifelse(f2>0, 1-f1*(n-1)*f1/n/((n-1)*f1+2*f2), 1-f1*(n-1)*(f1-1)/n/((n-1)*(f1-1)+2))
      W <- (1 - C1)/sum(X/n*(1-X/n)^n) 
      Prob.hat.Unse <- rep((1-C1)/f0.hat, f0.hat)
    } else if (datatype=="incidence_raw") {
      C1 = ifelse(f2>0, 1-f1/u*(n-1)*f1/((n-1)*f1+2*f2), 1-f1*(n-1)*(f1-1)/u/((n-1)*(f1-1)+2))
      W <- (1 - C1)/sum(X/u*(1-X/n)^n) 
      Prob.hat.Unse <- rep(u/n*(1-C1)/f0.hat, f0.hat)	
    }
    
    Prob.hat <- X/n*(1-W*(1-X/n)^n)
    Prob <- c(Prob.hat, Prob.hat.Unse)
    
    F.1 <- sum(dij[, X==1]) ; F.2 <- sum(dij[, X==2])
    F11 <- sum(dij[X==1, X==1]) ; F22 <- sum(dij[X==2, X==2])
    #
    if (datatype=="abundance") {
      F.0hat <- ifelse(F.2 > 0, ((n-1)/n) * (F.1^2/(2 * F.2)), ((n-1)/n)*(F.1*(F.1-0.01)/(2)))
      F00hat <- ifelse(F22 > 0, ((n-2)* (n-3)* (F11^2)/(4* n* (n-1)* F22)), ((n-2)* (n-3)* (F11*(F11-0.01))/(4 *n * (n-1))) )
    } else if (datatype=="incidence_raw") {
      F.0hat <- ifelse(F.2 > 0, ((n-1)/n) * (F.1^2/(2 * F.2)), ((n-1)/n)*(F.1*(F.1-0.01)/(2)))
      F00hat <- ifelse(F22 > 0, ((n-1)^2 * (F11^2)/(4* n* n* F22)), ((n-1)* (n-1)* (F11*(F11-0.01))/(4 *n * n)) )
    }
    if (f0.hat==0) {
      d=dij
    } else if (f0.hat==1) {
      random_dij = as.vector(rmultinom(1, 1000, rep(1/(length(X)*f0.hat), length(X)*f0.hat) ) )/1000
      d.0bar <- matrix(random_dij*F.0hat, length(X), f0.hat, byrow = T)
      d00 = matrix(0, f0.hat, f0.hat)
      d <- cbind(dij, d.0bar )
      aa <- cbind(t(d.0bar), d00 )
      d <- rbind(d, aa)
      diag(d) = 0
    } else {
      random_dij = as.vector(rmultinom(1, 1000, rep(1/(length(X)*f0.hat), length(X)*f0.hat) ) )/1000
      d.0bar <- matrix(random_dij*F.0hat, length(X), f0.hat, byrow = T)
      
      fo.num = (f0.hat * (f0.hat-1) )/2
      random_d00 = as.vector(rmultinom(1, 1000, rep(1/fo.num, fo.num) ) )/1000
      d00 = matrix(0, f0.hat, f0.hat)
      d00[upper.tri(d00)] = (F00hat/2)*random_d00
      d00 <- pmax(d00, t(d00))###signmatrix
      d <- cbind(dij, d.0bar )
      aa <- cbind(t(d.0bar), d00 )
      d <- rbind(d, aa)
      diag(d) = 0
    }
    return(list("pi" = Prob,"dij" = d))
  }
  dij <-  as.matrix(dij)
  out <- as.vector(dij)
  out <- out[out!=0]
  dmin <- min(out)
  dmax <- max(out)
  if (datatype=="incidence_raw") {
    data <- lapply(data, function(i) {
      c(ncol(i), rowSums(i))
    })
  }
  if (datatype=="abundance") {
    if(length(data)!=1){
      tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x/sum(x))), 1, mean)
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
    }else{
      tmp <- data[[1]]/sum(data[[1]])
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij)   
    }
  } else {
    if(length(data)!=1){
      tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x[-1]/sum(x[-1]))), 1, mean)
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
    }else{
      tmp <- data[[1]][-1]/sum(data[[1]][-1])
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij)   
    }
  }
  FD.CI = function(data, dij, tau, q, datatype){
    if (datatype == "abundance") {
      qFun = FD_MLE(data, dij, tau, q)
    } else {
      qFun = FD_MLE(data[-1], dij, tau, q)
    }
    if(boot!=0){
      BT = EstiBootComm.Func(data, dij, datatype)
      p_hat = BT[[1]]
      dij_boot = BT[[2]]
      dij_boot <-  as.matrix(dij_boot)
      dij_boot <- replace(dij_boot, dij_boot==0, 10^(-10))
      for (i in seq_len(nrow(dij_boot))) {
        dij_boot[i, i] <- 0
      }
      if (datatype=="abundance") {
        n=sum(data)
        Boot.X = rmultinom(boot, n, p_hat)
      } else {
        n=data[1]
        Boot.X = t(sapply(p_hat,function(i) rbinom(boot, n, i)))
      }
      qFun_sd = sd(sapply(seq_len(ncol(Boot.X)), function(i) {
        FD_MLE(Boot.X[, i], dij_boot, tau, q)
      }))
    }else{
      qFun_sd = 0
    }
    LCL = max(0, qFun - qnorm(0.975) * qFun_sd)
    UCL = qFun + qnorm(0.975) * qFun_sd
    a = round(c(qFun, qFun_sd, LCL, UCL), 4)
    a
  }
  
  Funq <- function(data, datatype){
    dminFDforq <- t(sapply(q, FUN = function(q) FD.CI(data, dij, dmin, q, datatype) ))
    dmaxFDforq <- t(sapply(q, FUN = function(q) FD.CI(data, dij, dmax, q, datatype) ))
    dmeanFDforq <-t(sapply(q, FUN = function(q) FD.CI(data, dij, dmean, q, datatype) ))
    out <- data.frame(rep(q,3), rbind(dminFDforq,dmaxFDforq,dmeanFDforq),rep(c("dmin","dmax","dmean"),each=length(q)))
  }
  Funtau <- function(data, datatype){
    q0FDfortau <- t(sapply(tau, FUN = function(tau) FD.CI(data, dij, tau, 0, datatype) ))
    q1FDfortau <- t(sapply(tau, FUN = function(tau) FD.CI(data, dij, tau, 1, datatype) ))
    q2FDfortau <- t(sapply(tau, FUN = function(tau) FD.CI(data, dij, tau, 2, datatype) ))
    out <- data.frame(rep(tau,3), rbind(q0FDfortau,q1FDfortau,q2FDfortau),rep(c("0","1","2"),each=length(tau)))
  }
  
  if(length(data)!=1){
    name = names(data)
    Outputforq <- data.frame(do.call(rbind,lapply(data, Funq, datatype=datatype)), rep(name, each=3*length(q)), row.names = NULL)
    Outputfortau <- data.frame(do.call(rbind,lapply(data, Funtau, datatype=datatype)), rep(name, each=3*length(tau)), row.names = NULL)
  }else{
    name = names(data)
    Outputforq <- data.frame(Funq(data[[1]], datatype), name, row.names = NULL)
    Outputfortau <- data.frame(Funtau(data[[1]], datatype), name, row.names = NULL)
  }
  colnames(Outputforq) <- c("q","estimate", "s.e.", "LCL", "UCL", "tau","site")
  colnames(Outputfortau) <- c("tau","estimate", "s.e.", "LCL", "UCL", "q","site")
  
  Output <- list(forq = Outputforq, fortau = Outputfortau)
  return(Output)
}
FunD(fskt, dij_fskt, seq(min(dij_fskt[dij_fskt>0]), max(dij_fskt), length.out = 5), c(0, 1, 2), 50, "abundance")
#-----------------------------------
# Plot FD q-profile and tau-profile
#-----------------------------------
#' PlotFD(out, data, dij, datatype) is a function to plot the output values of FunD
#' @param out a list with 2 elements. First is the q-profile matrix, the second is the tau-profile matrix.
#' @param data a list with N sites. Each element of list is species abundances or species-by-sampling unit incidence matrix.
#' @param dij a matrix of species-pairwise distances.
#' @param datatype data type of input data: individual-based abundance data (datatype = "abundance") or species-by-sampling-units incidence matrix (datatype = "incidence_raw").
#' @return two plots. The first one is the q-profile, the second is the tau-profile.
PlotFD <- function(out, data, dij, datatype){
  dij <-  as.matrix(dij)
  dij_v <- as.vector(dij)
  dij_v <- dij_v[dij_v!=0]
  dmin <- min(dij_v)
  dmax <- max(dij_v)
  if (datatype=="incidence_raw") {
    data <- lapply(data, function(i) {
      c(ncol(i), rowSums(i))
    })
  }
  if (datatype=="abundance") {
    if(length(data)!=1){
      tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x/sum(x))), 1, mean)
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
    }else{
      tmp <- data[[1]]/sum(data[[1]])
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij)   
    }
  } else {
    if(length(data)!=1){
      tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x[-1]/sum(x[-1]))), 1, mean)
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
    }else{
      tmp <- data[[1]][-1]/sum(data[[1]][-1])
      dmean <-  sum ( (tmp %*% t(tmp) ) * dij)   
    }
  }
  forq <- out[[1]]
  forq$tau = factor(forq$tau, levels = c("dmin","dmean","dmax"))
  fortau <- out[[2]]
  site <- unique(forq$site)
  if(length(site)==1){
    p1 <- ggplot(data=forq) + 
      geom_line(aes(x=q, y=estimate, color=tau, group=tau), linewidth=1.2)+
      geom_ribbon(aes(x=q,ymin=LCL, ymax=UCL, fill=tau, group=tau), alpha=0.3)+
      theme(legend.position = "bottom")+
      scale_color_discrete(name="method")+
      xlab("Order q")+ylab("Functional Diversity")+
      guides(fill="none", linetype="none")
    p2=ggplot(fortau) + 
      geom_line(aes(x=tau, y=estimate, color=q,group=q), linewidth=1.5)  + 
      geom_ribbon(aes(x=tau, ymin=LCL, ymax=UCL, fill=q), alpha=0.3)+
      annotate('text',x=dmin, y=0.1,label="d[min]" ,parse = TRUE,size=10) +
      annotate('text',x=dmean, y=0.1, label="d[mean]",parse = TRUE,size=10) +
      annotate('text',x=dmax, y=0.1, label="d[max]",parse = TRUE,size=10) +
      geom_vline(xintercept = c(dmin, dmean, dmax), linetype = "longdash",linewidth=0.5) +
      theme(legend.position = "bottom")+
      xlab("Tau (Level of threshold distinctiveness)")+ylab("Functional Diversity")+
      scale_color_discrete(name="method")+
      guides(fill="none")
  }else{
    p1 <- ggplot(forq, aes(x=q, y=estimate, color=site,linetype=site)) + geom_line(linewidth=1.5)  + 
      geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=site), alpha=0.3, linetype=0) + facet_wrap(~tau, scales = "free")+
      theme(legend.key.width=unit(2,"cm"))
    p1 <- ggplot(data=forq) + 
      geom_line(aes(x=q, y=estimate, color=site, group=site, linetype=site), linewidth=1.2)+
      geom_ribbon(aes(x=q,ymin=LCL, ymax=UCL, fill=site, group=site), alpha=0.3)+
      theme(legend.position = "bottom")+
      scale_color_discrete(name="method")+
      #coord_cartesian(xlim = c(0, 1))+
      facet_wrap(~tau, scales = "free")+
      #xlim("0", "1")+
      guides(fill="none", linetype="none")
    p2 <- ggplot(fortau, aes(x=tau, y=estimate, color=site,linetype=site)) + geom_line(linewidth=1.5)  + 
      geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=site), alpha=0.3, linetype=0)+
      annotate('text',x=dmin, y=0.1,label="d[min]" ,parse = TRUE,size=10) +
      annotate('text',x=dmean, y=0.1, label="d[mean]",parse = TRUE,size=10) +
      annotate('text',x=dmax, y=0.1, label="d[max]",parse = TRUE,size=10) +
      geom_vline(xintercept = c(dmin,dmean,dmax), linetype = "longdash",linewidth=0.5) +
      facet_wrap(~q, scales = "free")+
      theme(legend.key.width=unit(2,"cm"))
    p2=ggplot(fortau) + 
      geom_line(aes(x=tau, y=estimate, color=site,group=site, linetype=site), linewidth=1.5)  + 
      geom_ribbon(aes(x=tau, ymin=LCL, ymax=UCL, fill=site), alpha=0.3)+
      annotate('text',x=dmin, y=0.1,label="d[min]" ,parse = TRUE,size=10) +
      annotate('text',x=dmean, y=0.1, label="d[mean]",parse = TRUE,size=10) +
      annotate('text',x=dmax, y=0.1, label="d[max]",parse = TRUE,size=10) +
      geom_vline(xintercept = c(dmin,
                                dmean,
                                dmax), linetype = "longdash",linewidth=0.5) +
      facet_wrap(~q, scales = "free")+
      theme(legend.position = "bottom")+
      labs(x="tau")+
      scale_color_discrete(name="method")+
      guides(fill="none", linetype="none")
    p1 <-  p1 +xlab("Order q")+ylab("Functional Diversity") 
    p2 <-  p2 +xlab("Tau (Level of threshold distinctiveness)")+ylab("Functional Diversity")
  }
  output <- list(pic1= p1, pic2 = p2)
  return(output)
}
PlotFD(FunD(fskt, dij_fskt, seq(min(dij_fskt[dij_fskt>0]), max(dij_fskt), length.out = 5), c(0, 1, 2), 50, "abundance"), fskt, dij_fskt, "abundance")


#-----------------------------------------------------------------
# Functional Beta/Dissimilarity for specified values of tau and q
#------------------------------------------------------------------
#' FD_Beta(data, dij, tau, q, type) is a function of obtaining functional dissimilarity measure of order q.
#' @param data a SxN matrix of species sample frequencies with S species, N sites.
#' @param dij a matrix of species pairwise distances.
#' @param tau a numeric of a specified level of threshold distinctiveness.
#' @param q a numeric of a specified diversity order.
#' @param type a character to choose index: "C" for 1-CqN ; "U" for 1-UqN.
#' @return a numerical of functional dissimilarity.
FD_Beta <- function(data, dij, tau, q, type){
  N <- ncol(data)
  if (class(dij)[1]=="matrix") {
    dij <- lapply(seq_len(N), function(i) dij)
  }
  dij <- lapply(seq_along(dij), function(i) {
    d <- as.matrix(dij[[i]])
    d[which(d>tau,arr.ind = T)] <- tau
    d
  })
  aik <- sapply(seq_along(dij), function(i) {
    (1-dij[[i]]/tau) %*% data[, i] 
  })
  aiplus <- apply(aik, 1, sum)
  vi <- apply(data, 1, sum)/aiplus  
  alpha_v <- rep(vi, N)
  nplus <- sum(data)
  aik <- as.vector(aik)
  alpha_v <- alpha_v[aik!=0]
  aik <- aik[aik!=0]
  if(q==1){
    gamma=exp(sum(-vi*aiplus/nplus*log(aiplus/nplus)))
    alpha=1/N*exp(sum(-alpha_v*aik/nplus*log(aik/nplus)))
    beta = max(1, gamma/alpha)
    out <- log(beta)/log(N)
  }else{
    gamma=(sum(vi*(aiplus/nplus)^q))^(1 / (1-q))
    alpha=1/N*(sum(alpha_v*(aik/nplus)^q))^(1 / (1-q))
    beta = max(1, gamma/alpha)
    if(type == "C") out <- (1-(beta)^(1-q))/(1-(N)^(1-q))
    if(type == "U") out <- (1-(beta)^(q-1))/(1-(N)^(q-1))
  }
  out
}
FD_Beta(do.call(cbind, fskt), dij_fskt, max(dij_fskt), 2, "C")

#----------------------------------------------------------
# Functional Beta/Dissimilarity for various values of tau and q
#-----------------------------------------------------------
#' MultiFD.Beta(data, dij, tau, q, boot, datatype, relative = T) is a function of obtaining functional dissimilarity tau-table or q-table.
#' @param data a list with N sites. Each element of list is species abundances or species-by-sampling-unit incidence matrix.
#' @param dij a matrix of species pairwise distances.
#' @param tau a numeric or a vector of levels of threshold distinctiveness. 
#' @param q a numeric or a vector of diversity order; The suggested range for q is [0, 3].
#' @param boot a numeric of number of bootstrap replications.
#' @param datatype data type of input data: individual-based abundance data (datatype = "abundance") or species-by-sampling-unit incidence matrix (datatype = "incidence_raw").
#' @param relative data are species relative abundances or absolute abundances, default is TRUE
#' @return four matrices of functional dissimilarity measures: (1) 1-CqN q-profile table, (2) 1-UqN q-profile table, (3) 1-CqN tau-profile table, and (4) 1-UqN tau-profile table
MultiFD.Beta <- function(data, dij, tau, q, boot, datatype, relative = TRUE){
  Boots.beta = function(data, dij, datatype){
    #data is pooled matrix
    if(datatype == "abundance"){
      N = ncol(data); n = colSums(data)
      pool=rowSums(data) ; OBS=length(pool[pool>0])
      data=data[pool>0,]
      obs=sapply(1:N,function(k) length(data[,k][data[,k]>0]))
      F1=sum(pool==1);F2=sum(pool==2)
      F0=ceiling((sum(n)-1)/sum(n)*ifelse(F2==0,F1*(F1-1)/2,F1^2/(2*F2)))
      f1=sapply(1:N,function(k) sum(data[,k]==1))
      f2=sapply(1:N,function(k) sum(data[,k]==2))
      C=sapply(1:N, function(k) {
        ifelse(f2[k]>0, 1-f1[k]*(n[k]-1)*f1[k]/n[k]/((n[k]-1)*f1[k]+2*f2[k]), 
               1-f1[k]*(n[k]-1)*(f1[k]-1)/n[k]/((n[k]-1)*(f1[k]-1)+2))
      })
      f0=ceiling(sapply(1:N,function(k) ((n[k]-1)/n[k])*ifelse(f2[k]==0,f1[k]*(f1[k]-1)/2,f1[k]^2/(2*f2[k]))))
      r.data=sapply(1:N,function(k) data[,k]/n[k])
      W=sapply(1:N,function(k) (1-C[k])/sum(r.data[,k]*(1-r.data[,k])^n[k]))
      if(F0>0){boots.pop=rbind(r.data,matrix(0,ncol=N,nrow=F0))
      }else{boots.pop=r.data}
      for(i in 1:N)
      {
        if(f0[i]>0)
        {
          f0[i]=ifelse(f0[i]+obs[i]>OBS+F0, OBS+F0-obs[i],f0[i])
          boots.pop[,i][1:OBS]=boots.pop[,i][1:OBS]*(1-W[i]*(1-boots.pop[,i][1:OBS])^n[i])   #
          I=which(boots.pop[,i]==0);II=sample(I,f0[i])
          boots.pop[II,i]=rep((1-C[i])/f0[i],f0[i])
        }
      }
    }
    if(datatype == "incidence_raw"){
      dat = data ; n = data[1, ] %>% as.numeric()
      data = data[-1, ]
      Tt=sum(n) ; N=ncol(data)
      pool=rowSums(data);OBS=length(pool[pool>0]);
      data=data[pool>0,];
      obs=sapply(1:N,function(k) length(data[,k][data[,k]>0]));
      Q1=sum(pool==1);Q2=sum(pool==2);
      Q0=ceiling((Tt-1)/Tt*ifelse(Q2==0,Q1*(Q1-1)/2,Q1^2/(2*Q2)))
      q1=sapply(1:N,function(k) sum(data[,k]==1));
      q2=sapply(1:N,function(k) sum(data[,k]==2));
      C=sapply(1:N, function(k) {
        ifelse(q2[k]>0, 1-q1[k]*(n[k]-1)*q1[k]/Tt/((n[k]-1)*q1[k]+2*q2[k]), 
               1-q1[k]*(n[k]-1)*(q1[k]-1)/Tt/((n[k]-1)*(q1[k]-1)+2))
      })
      q0=ceiling(sapply(1:N,function(k) ((n[k]-1)/n[k])*ifelse(q2[k]==0,q1[k]*(q1[k]-1)/2,q1[k]^2/(2*q2[k]))));
      r.data=sapply(1:N,function(k) data[,k]/n[k]);
      W=sapply(1:N,function(k) Tt/n[k]*(1-C[k])/sum(r.data[,k]*(1-r.data[,k])^n[k]))
      if(Q0>0){ boots.pop=rbind(r.data,matrix(0,ncol=N,nrow=Q0))
      }else{boots.pop=r.data}
      for(i in 1:N)
      {
        if(q0[i]>0)
        {
          q0[i]=ifelse(q0[i]+obs[i]>OBS+Q0, OBS+Q0-obs[i],q0[i])
          boots.pop[,i][1:OBS]=boots.pop[,i][1:OBS]*(1-W[i]*(1-boots.pop[,i][1:OBS])^n[i])
          I=which(boots.pop[,i]==0);II=sample(I,q0[i])
          boots.pop[II,i]=rep((q1[i]/n[i])/q0[i],q0[i])
        }
      }
    }
    d <- lapply(seq_along(dij), function(i) {
      X = data[, i];# names(X) <- rownames(data)
      dij=dij[[i]]; #dij=dij[names(X), names(X)]
      F.1 <- sum(dij[, X==1]) ; F.2 <- sum(dij[, X==2])
      F11 <- sum(dij[X==1, X==1]) ; F22 <- sum(dij[X==2, X==2])
      if (datatype=="abundance") {
        F.0hat <- ifelse(F.2 > 0, ((n[i]-1)/n[i]) * (F.1^2/(2 * F.2)), ((n[i]-1)/n[i])*(F.1*(F.1-0.01)/(2)))
        F00hat <- ifelse(F22 > 0, ((n[i]-2)* (n[i]-3)* (F11^2)/(4* n[i]* (n[i]-1)* F22)), ((n[i]-2)* (n[i]-3)* (F11*(F11-0.01))/(4 *n[i] * (n[i]-1))) )
        f0=F0
      } else if (datatype=="incidence_raw") {
        F.0hat <- ifelse(F.2 > 0, ((n[i]-1)/n[i]) * (F.1^2/(2 * F.2)), ((n[i]-1)/n[i])*(F.1*(F.1-0.01)/(2)))
        F00hat <- ifelse(F22 > 0, ((n[i]-1)^2 * (F11^2)/(4* n[i]* n[i]* F22)), ((n[i]-1)* (n[i]-1)* (F11*(F11-0.01))/(4 *n[i] * n[i])) )
        f0=Q0
      }
      
      if (f0==0) {
        d <- dij
      } else if (f0==1) {
        random_dij = as.vector(rmultinom(1, 1000, rep(1/(length(X)*f0), length(X)*f0) ) )/1000
        d.0bar <- matrix(random_dij*F.0hat, length(X), f0, byrow = T)
        d00 = matrix(0, f0, f0)
        d <- cbind(dij, d.0bar ) %>% as.matrix()
        aa <- cbind(t(d.0bar), d00 ) %>% as.matrix()
        d <- rbind(d, aa)
        diag(d) = 0
      } else {
        random_dij = as.vector(rmultinom(1, 1000, rep(1/(length(X)*f0), length(X)*f0) ) )/1000
        d.0bar <- matrix(random_dij*F.0hat, length(X), f0, byrow = T)
        f0.num = (f0 * (f0-1) )/2
        random_d00 = as.vector(rmultinom(1, 1000, rep(1/f0.num, f0.num) ) )/1000
        d00 = matrix(0, f0, f0)
        d00[upper.tri(d00)] = (F00hat/2)*random_d00
        d00 <- pmax(d00, t(d00))###signmatrix
        d <- cbind(dij, d.0bar ) %>% as.matrix()
        aa <- cbind(t(d.0bar), d00 ) %>% as.matrix()
        d <- rbind(d, aa)
        diag(d) = 0
      }
      d
    })
    return(list(pop=boots.pop, dij=d))
  }
  dij <-  as.matrix(dij)
  dij_v <- as.vector(dij)
  dij_v <- dij_v[dij_v!=0]
  dmin <- min(dij_v)
  dmax <- max(dij_v)
  if (datatype=="incidence_raw") {
    data <- lapply(data, function(i) {
      c(ncol(i), rowSums(i))
    })
    data2=data
    if (relative==T) {
      data2 <- lapply(data, function(i) {
        c(i[1], i[-1]/sum(i[-1]))
      })
    }
    tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x[-1]/sum(x[-1]))), 1, mean)
    dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
  } else if (datatype=="abundance") {
    data2=data
    if (relative==T) {
      data2 <- lapply(data, function(i) i/sum(i))
    }
    tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x/sum(x))), 1, mean)
    dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
  } 
  data <- do.call(cbind, data)
  data2 <- do.call(cbind, data2)
  dij <- lapply(1:ncol(data), function(i) dij)
  FD.CI.beta = function(data, dij, tau, q, datatype, type, boot, data2){
    if (datatype == "abundance") {
      qFun = FD_Beta(data2, dij, tau, q, type)
    } else {
      da_inc <- data2[-1, ]
      qFun = FD_Beta(da_inc, dij, tau, q, type)
    }
    if(boot!=0){
      BT = Boots.beta(data, dij, datatype)
      p_hat = BT$pop
      dij_boot = BT$dij
      dij_boot <- lapply(dij_boot, function(i) {
        d <- as.matrix(i)
        d <- replace(d, d==0, 10^(-10))
        diag(d) <- 0
        d
      })
      
      if (datatype=="abundance") {
        n=colSums(data)
        Boot.X <- sapply(seq_along(n), function(i) rmultinom(boot, n[i], p_hat[, i]), simplify = "array")
      } else {
        n=data[1, ]
        Boot.X <- sapply(seq_along(n), function(i) t(sapply(p_hat[, i], function(j) rbinom(boot, n[i], j))), simplify = "array")
      }
      boot.sample <- array(0, dim=c(dim(Boot.X)[1], dim(Boot.X)[3], dim(Boot.X)[2]))
      for (i in seq_len(dim(Boot.X)[2])) {
        for (j in seq_len(dim(Boot.X)[3])) {
          boot.sample[, j, i] <- Boot.X[, i, j]
        }
      }
      qFun_sd = sd(sapply(seq_len(dim(boot.sample)[3]), function(i) {
        da.sample <- boot.sample[,,i][rowSums(boot.sample[,,i])>0, ]
        dij_sample <- lapply(dij_boot, function(j) j[rowSums(boot.sample[,,i])>0, rowSums(boot.sample[,,i])>0])
        FD_Beta(da.sample, dij_sample, tau, q, type)
      }))
    }else{
      qFun_sd = 0
    }
    LCL = max(0, qFun - qnorm(0.975) * qFun_sd)
    UCL = qFun + qnorm(0.975) * qFun_sd
    a = round(c(qFun, qFun_sd, LCL, UCL), 4)
    a
  }
  
  Funq <- function(data, type, datatype, data2){
    dminFDforq <- t(sapply(q, FUN = function(q) FD.CI.beta(data, dij, dmin, q, datatype, type, boot, data2) ))
    dmaxFDforq <- t(sapply(q, FUN = function(q) FD.CI.beta(data, dij, dmax, q, datatype, type, boot, data2) ))
    dmeanFDforq <- t(sapply(q, FUN = function(q) FD.CI.beta(data, dij, dmean, q, datatype, type, boot, data2) ))
    out <- data.frame(rep(q,3), rbind(dminFDforq,dmeanFDforq,dmaxFDforq),rep(c("dmin","dmean","dmax"),each=length(q)))
    out
  }
  Funtau <- function(data, type, datatype, data2){
    q0FDfortau <- t(sapply(tau, FUN = function(i) FD.CI.beta(data, dij, i, 0, datatype, type, boot, data2) ))
    q1FDfortau <- t(sapply(tau, FUN = function(i) FD.CI.beta(data, dij, i, 1, datatype, type, boot, data2) ))
    q2FDfortau <- t(sapply(tau, FUN = function(i) FD.CI.beta(data, dij, i, 2, datatype, type, boot, data2) ))
    out <- data.frame(rep(tau,3), rbind(q0FDfortau,q1FDfortau,q2FDfortau),rep(c("0","1","2"),each=length(tau)))
    out
  }
  
  CqNforq <- data.frame(Funq(data,"C", datatype, data2), row.names = NULL)
  CqNfortau <- data.frame(Funtau(data,"C", datatype, data2), row.names = NULL)
  UqNforq <- data.frame(Funq(data,"U", datatype, data2), row.names = NULL)
  UqNfortau <- data.frame(Funtau(data,"U", datatype, data2), row.names = NULL)
  colnames(CqNforq) <- colnames(UqNforq) <- c("q","estimate", "s.e.", "LCL", "UCL", "tau")
  colnames(CqNfortau) <- colnames(UqNfortau) <- c("tau","estimate", "s.e.", "LCL", "UCL", "q")
  Output <- list(CqNforq=CqNforq, UqNforq=UqNforq, CqNfortau=CqNfortau, UqNfortau=UqNfortau)
  return(Output)
}
MultiFD.Beta(fskt, dij_fskt, seq(min(dij_fskt[dij_fskt>0]), max(dij_fskt), length.out = 5), c(0, 1, 2), 50, "abundance", T)
#--------------------------------------------------------
# Plot functional dissimilarity q-profile and tau-profile
#--------------------------------------------------------
#' PlotBeta(out) is a function to plot the output values of MultiFD.Beta
#' @param out a list with four elements: (1) 1-CqN q-profile table, (2) 1-UqN q-profile table, (3) 1-CqN tau-profile table, and (4) 1-UqN tau-profile table
#' @param data a list with N sites; each element of list is species abundances or species-by-sampling-unit incidence matrix.
#' @param dij a matrix of species-pairwise distances.
#' @param datatype data type of input data: individual-based abundance data (datatype = "abundance") or species-by-sampling-unit incidence matrix (datatype = "incidence_raw").
#' @return four plots. (1) 1-CqN q-profile, (2) 1-UqN q-profile, (3) 1-CqN tau-profile, and (4) 1-UqN tau-profile
PlotBeta <- function(out, data, dij, datatype, relative = TRUE){
  dij <-  as.matrix(dij)
  dij_v <- as.vector(dij)
  dij_v <- dij_v[dij_v!=0]
  dmin <- min(dij_v)
  dmax <- max(dij_v)
  if (datatype=="incidence_raw") {
    data <- lapply(data, function(i) {
      c(ncol(i), rowSums(i))
    })
    if (relative==T) {
      data <- lapply(data, function(i) {
        c(i[1], i[-1]/sum(i[-1]))
      })
    }
    tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x[-1]/sum(x[-1]))), 1, mean)
    dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
  } else if (datatype=="abundance") {
    if (relative==T) {
      data <- lapply(data, function(i) i/sum(i))
    }
    tmp <- apply(do.call(cbind,lapply(data, FUN = function(x) x/sum(x))), 1, mean)
    dmean <-  sum ( (tmp %*% t(tmp) ) * dij) 
  } 
  Cforq <- out[[1]]
  Uforq <- out[[2]]
  Cfortau <- out[[3]]
  Ufortau <- out[[4]]
  Cforq <-   Cforq %>% transmute(q=q, estimate=estimate, 
                                 s.e.=s.e., LCL=ifelse(LCL>=1, 1, ifelse(LCL<=0, 0, LCL)), 
                                 UCL=ifelse(UCL>=1, 1, ifelse(UCL<=0, 0, UCL)), tau=tau)
  Uforq <- Uforq %>% transmute(q=q, estimate=estimate, 
                               s.e.=s.e., LCL=ifelse(LCL>=1, 1, ifelse(LCL<=0, 0, LCL)), 
                               UCL=ifelse(UCL>=1, 1, ifelse(UCL<=0, 0, UCL)), tau=tau)
  Cfortau <- Cfortau %>% transmute(tau=tau, estimate=estimate, 
                                   s.e.=s.e., LCL=ifelse(LCL>=1, 1, ifelse(LCL<=0, 0, LCL)), 
                                   UCL=ifelse(UCL>=1, 1, ifelse(UCL<=0, 0, UCL)), q=q)
  Ufortau <- Ufortau %>% transmute(tau=tau, estimate=estimate, 
                                   s.e.=s.e., LCL=ifelse(LCL>=1, 1, ifelse(LCL<=0, 0, LCL)), 
                                   UCL=ifelse(UCL>=1, 1, ifelse(UCL<=0, 0, UCL)), q=q)
  p1Cforq <- ggplot(Cforq, aes(x=q, y=estimate, color=tau, linetype=tau)) + 
    geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=tau, colour=NULL), alpha=0.3, linetype=0)+
    geom_line(linewidth=1.5) + 
    theme(legend.position = "bottom")+
    ylim(c(0,1))
  p1Uforq <- ggplot(Uforq, aes(x=q, y=estimate, color=tau, linetype=tau)) + 
    geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=tau, colour=NULL), alpha=0.3, linetype=0)+
    geom_line(linewidth=1.5) + 
    theme(legend.position = "bottom")+
    ylim(c(0,1))
  p2Cfortau <- ggplot(Cfortau, aes(x=tau, y=estimate, color=q, linetype=q)) + 
    geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=q, colour=NULL), alpha=0.3, linetype=0)+
    geom_line(linewidth=1.5) + 
    theme(legend.position = "bottom")+
    annotate('text',x=dmin, y=0.1,label="d[min]" ,parse = TRUE,size=10) +
    annotate('text',x=dmean, y=0.1, label="d[mean]",parse = TRUE,size=10) +
    annotate('text',x=dmax, y=0.1, label="d[max]",parse = TRUE,size=10) +
    geom_vline(xintercept = c(dmin, dmean, dmax), linetype = "longdash",linewidth=0.5) +
    ylim(c(0,1))
  p2Ufortau <- ggplot(Ufortau, aes(x=tau, y=estimate, color=q, linetype=q)) + 
    geom_ribbon(aes(ymin=LCL, ymax=UCL, fill=q, colour=NULL), alpha=0.3, linetype=0)+
    geom_line(linewidth=1.5) + 
    theme(legend.position = "bottom")+
    annotate('text',x=dmin, y=0.1,label="d[min]" ,parse = TRUE,size=10) +
    annotate('text',x=dmean, y=0.1, label="d[mean]",parse = TRUE,size=10) +
    annotate('text',x=dmax, y=0.1, label="d[max]",parse = TRUE,size=10) +
    geom_vline(xintercept = c(dmin, dmean, dmax), linetype = "longdash",linewidth=0.5) +
    ylim(c(0,1))
  p1Cforq <-  p1Cforq +xlab("Order q")+ylab("Functional Differentiation") +theme(text=element_text(size=18),legend.position="bottom")+
    theme(legend.key.width=unit(2,"cm"))
  p1Uforq <-  p1Uforq +xlab("Order q")+ylab("Functional Differentiation") +theme(text=element_text(size=18),legend.position="bottom")+
    theme(legend.key.width=unit(2,"cm"))
  p2Cfortau <-  p2Cfortau +xlab("Tau (Level of threshold distinctiveness)")+ylab("Functional Differentiation")+theme(text=element_text(size=18),legend.position="bottom")+
    theme(legend.key.width=unit(2,"cm"))
  p2Ufortau <-  p2Ufortau +xlab("Tau (Level of threshold distinctiveness)")+ylab("Functional Differentiation")+theme(text=element_text(size=18),legend.position="bottom")+
    theme(legend.key.width=unit(2,"cm"))
  output <- list(p1Cforq= p1Cforq, p1Uforq = p1Uforq, p2Cfortau=p2Cfortau, p2Ufortau=p2Ufortau )
  return(output)
}
PlotBeta(MultiFD.Beta(fskt, dij_fskt, seq(min(dij_fskt[dij_fskt>0]), max(dij_fskt), length.out = 17), seq(0, 2, length.out =17), 100, "abundance", T), 
         fskt, dij_fskt, "abundance", T)