R/dbFD_joggle.R

dbFD_joggle <- function (x, a, w, w.abun = TRUE, stand.x = TRUE, ord = c("podani",
                                                          "metric"), asym.bin = NULL, corr = c("sqrt", "cailliez",
                                                                                               "lingoes", "none"), calc.FRic = TRUE, m = "max", stand.FRic = FALSE,
          scale.RaoQ = FALSE, calc.FGR = FALSE, clust.type = "ward",
          km.inf.gr = 2, km.sup.gr = nrow(x) - 1, km.iter = 100, km.crit = c("calinski",
                                                                             "ssi"), calc.CWM = TRUE, CWM.type = c("dom", "all"),
          calc.FDiv = TRUE, dist.bin = 2, print.pco = FALSE, messages = TRUE)
{
  tol <- .Machine$double.eps
  corr <- match.arg(corr)
  ord <- match.arg(ord)
  CWM.type <- match.arg(CWM.type)
  km.crit <- match.arg(km.crit)
  if (!is.logical(messages))
    stop("'messages' must be TRUE or FALSE.", "\n")
  if (!is.logical(stand.FRic))
    stop("'stand.FRic' must be TRUE or FALSE.", "\n")
  if (!is.logical(stand.x))
    stop("'stand.x' must be TRUE or FALSE.", "\n")
  if (!is.logical(w.abun))
    stop("'w.abun' must be TRUE or FALSE.", "\n")
  if (!is.logical(calc.FRic))
    stop("'calc.FRic' must be TRUE or FALSE.", "\n")
  if (!is.logical(calc.FDiv))
    stop("'calc.FDiv' must be TRUE or FALSE.", "\n")
  if (!is.logical(calc.FGR))
    stop("'calc.FGR' musts be TRUE or FALSE.", "\n")
  if (!is.logical(calc.CWM))
    stop("'calc.CWM' must be TRUE or FALSE.", "\n")
  if (!is.logical(scale.RaoQ))
    stop("'scale.RaoQ' must be TRUE or FALSE.", "\n")
  if (!is.logical(print.pco))
    stop("'print.pco' must be TRUE or FALSE.", "\n")
  if (is.matrix(x) | is.data.frame(x)) {
    is.dist.x <- FALSE
    s.x <- dim(x)[1]
    t.x <- dim(x)[2]
    if (is.null(row.names(x)))
      stop("'x' must have row names.", "\n")
    else x.rn <- row.names(x)
  }
  if (is.vector(x) | is.factor(x)) {
    is.dist.x <- FALSE
    s.x <- length(x)
    t.x <- 1
    if (is.null(names(x)))
      stop("'x' must have names.", "\n")
    else x.rn <- names(x)
  }
  if (class(x)[1] == "dist" | class(x)[1] == "dissimilarity") {
    is.dist.x <- TRUE
    s.x <- attr(x, "Size")
    t.x <- 1
    if (is.null(attr(x, "Labels")))
      stop("'x' must have labels.", "\n")
    else x.rn <- attr(x, "Labels")
  }
  if (missing(a)) {
    ab.names <- list("Community1", x.rn)
    a <- matrix(1, 1, s.x, dimnames = ab.names)
  }
  else {
    if (is.matrix(a) | is.data.frame(a)) {
      s.a <- dim(a)[2]
      ab.t <- t(a)
      if (is.null(row.names(ab.t)))
        stop("'a' must have column names.", "\n")
      else ab.t.row <- row.names(ab.t)
      a <- as.matrix(a)
    }
    if (is.vector(a)) {
      s.a <- length(a)
      if (is.null(names(a)))
        stop("'a' must have names.", "\n")
      else ab.t.row <- names(a)
      ab.names <- list("Community1", ab.t.row)
      a <- matrix(a, 1, s.a, dimnames = ab.names)
    }
    if (s.x != s.a)
      stop("Different number of species in 'x' and 'a'.",
           "\n")
    if (any(ab.t.row != x.rn))
      stop("Species labels in 'x' and 'a' need to be identical and ordered alphabetically (or simply in the same order).",
           "\n")
  }
  a <- as.matrix(a)
  a[which(is.na(a))] <- 0
  abun.sum <- apply(a, 1, sum)
  if (any(abun.sum == 0))
    stop("At least one community has zero-sum abundances (no species).",
         "\n")
  abun.sum2 <- apply(a, 2, sum)
  if (any(abun.sum2 == 0))
    stop("At least one species does not occur in any community (zero total abundance across all communities).",
         "\n")
  if (!missing(w) & is.dist.x)
    stop("When 'x' is a distance matrix, 'w' should be left missing.",
         "\n")
  if (!missing(w) & !is.dist.x) {
    if (!is.numeric(w) | length(w) != t.x)
      stop("'w' should be a numeric vector of length = number of traits.",
           "\n")
    else w <- w/sum(w)
  }
  if (missing(w))
    w <- rep(1, t.x)/sum(rep(1, t.x))
  if (is.matrix(x) | is.data.frame(x)) {
    x <- data.frame(x)
    if (t.x >= 2) {
      x.class <- sapply(x, data.class)
      if (any(x.class == "character"))
        x[, x.class == "character"] <- as.factor(x[,
                                                   x.class == "character"])
      else x <- x
      if (all(x.class == "numeric") & all(!is.na(x))) {
        if (length(unique(w)) == 1) {
          x.s <- apply(x, 2, scale, center = TRUE, scale = stand.x)
          x.dist <- dist(x.s)
        }
        else {
          x.dist <- FD::gowdis(x, w = w, ord = ord, asym.bin = asym.bin)
        }
      }
      else {
        x.dist <- FD::gowdis(x, w = w, ord = ord, asym.bin = asym.bin)
      }
    }
    if (t.x == 1) {
      if (is.numeric(x[, 1])) {
        if (all(!is.na(x))) {
          x.s <- apply(x, 2, scale, center = TRUE, scale = stand.x)
          x.dist <- dist(x.s)
        }
        if (any(is.na(x))) {
          pos.NA <- which(is.na(x), arr.ind = TRUE)
          x <- na.omit(x)
          x.s <- apply(x, 2, scale, center = TRUE, scale = stand.x)
          x.dist <- dist(x.s)
          row.excl.ab <- pos.NA[, 1]
          a <- a[, -row.excl.ab]
          if (messages)
            cat("Warning: Species with missing trait values have been excluded.",
                "\n")
        }
      }
      if (is.factor(x[, 1]) | is.character(x[, 1])) {
        if (is.ordered(x[, 1]))
          x <- x
        else x[, 1] <- as.factor(x[, 1])
        if (any(is.na(x))) {
          pos.NA <- which(is.na(x), arr.ind = TRUE)
          x <- na.omit(x)
          row.excl.ab <- pos.NA[, 1]
          a <- a[, -row.excl.ab]
          x.rn <- x.rn[-pos.NA]
          if (messages)
            cat("Warning: Species with missing trait values have been excluded.",
                "\n")
        }
        if (is.ordered(x[, 1])) {
          x.s <- data.frame(rank(x[, 1]))
          names(x.s) <- x.rn
          x.dist <- dist(x.s)
        }
        else {
          x.f <- as.factor(x[, 1])
          x.dummy <- diag(nlevels(x.f))[x.f, ]
          x.dummy.df <- data.frame(x.dummy, row.names = x.rn)
          sequence <- 1:10
          if (all(dist.bin != sequence[any(sequence)]))
            stop("'dist.bin' must be an integer between 1 and 10.",
                 "\n")
          x.dist <- dist.binary(x.dummy.df, method = dist.bin)
        }
      }
    }
  }
  if (is.vector(x) & is.numeric(x)) {
    if (any(is.na(x))) {
      pos.NA <- which(is.na(x))
      x <- na.omit(x)
      a <- a[, -pos.NA]
      x.rn <- x.rn[-pos.NA]
      if (messages)
        cat("Warning: Species with missing trait values have been excluded.",
            "\n")
    }
    else x <- x
    x.s <- scale(x, center = T, scale = stand.x)
    x.dist <- dist(x.s)
    x <- data.frame(x)
    dimnames(x) <- list(x.rn, "Trait")
  }
  if (is.vector(x) & is.character(x)) {
    x <- as.factor(x)
    if (any(is.na(x))) {
      pos.NA <- which(is.na(x))
      x <- na.omit(x)
      a <- a[, -pos.NA]
      x.rn <- x.rn[-pos.NA]
      if (messages)
        cat("Warning: Species with missing trait values have been excluded.",
            "\n")
    }
    else x <- x
    dimnames(x) <- list(x.rn, "Trait")
    x.dummy <- diag(nlevels(x))[x, ]
    x.dummy.df <- data.frame(x.dummy, row.names = x.rn)
    sequence <- 1:10
    if (all(dist.bin != sequence[any(sequence)]))
      stop("'dist.bin' must be an integer between 1 and 10.",
           "\n")
    x <- data.frame(x)
    x.dist <- dist.binary(x.dummy.df, method = dist.bin)
  }
  if (is.ordered(x)) {
    if (any(is.na(x))) {
      pos.NA <- which(is.na(x))
      x <- na.omit(x)
      a <- a[, -pos.NA]
      x.rn <- x.rn[-pos.NA]
      cat("Warning: Species with missing trait values have been excluded.",
          "\n")
    }
    else x <- x
    x <- data.frame(x)
    dimnames(x) <- list(x.rn, "Trait")
    x.dist <- FD::gowdis(x, w = w, ord = ord, asym.bin = asym.bin)
  }
  if (is.factor(x) & !is.ordered(x)) {
    if (any(is.na(x))) {
      pos.NA <- which(is.na(x))
      x <- na.omit(x)
      a <- a[, -pos.NA]
      x.rn <- x.rn[-pos.NA]
      if (messages)
        cat("Warning: Species with missing trait values have been excluded.",
            "\n")
    }
    else x <- x
    x.dummy <- diag(nlevels(x))[x, ]
    x.dummy.df <- data.frame(x.dummy, row.names = x.rn)
    sequence <- 1:10
    if (all(dist.bin != sequence[any(sequence)]))
      stop("'dist.bin' must be an integer between 1 and 10.",
           "\n")
    x.dist <- dist.binary(x.dummy.df, method = dist.bin)
    x <- data.frame(x)
    dimnames(x) <- list(x.rn, "Trait")
  }
  if (class(x)[1] == "dist" | class(x)[1] == "dissimilarity") {
    if (any(is.na(x)))
      stop("When 'x' is a distance matrix, it cannot have missing values (NA).",
           "\n")
    x.dist <- x
  }
  if (any(is.na(x.dist)))
    stop("NA's in the distance matrix.", "\n")
  if (!is.dist.x) {
    no.traits <- apply(x, 1, function(v) length(v[!is.na(v)]))
    if (any(no.traits == 0))
      stop("At least one species has no trait data.",
           "\n")
  }
  c <- dim(a)[1]
  if (!w.abun)
    for (h in 1:c) {
      abpos <- which(a[h, ] > 0)
      a[h, abpos] <- 1
    }
  attr(x.dist, "Labels") <- x.rn
  if (is.euclid(x.dist))
    x.dist2 <- x.dist
  if (!is.euclid(x.dist)) {
    if (corr == "lingoes") {
      x.dist2 <- lingoes(x.dist)
      if (messages)
        cat("Species x species distance matrix was not Euclidean. Lingoes correction was applied.",
            "\n")
    }
    if (corr == "cailliez") {
      x.dist2 <- cailliez(x.dist)
      if (messages)
        cat("Species x species distance matrix was not Euclidean. Cailliez correction was applied.",
            "\n")
    }
    if (corr == "sqrt") {
      x.dist2 <- sqrt(x.dist)
      if (!is.euclid(x.dist2))
        stop("Species x species distance matrix was still is not Euclidean after 'sqrt' correction. Use another correction method.",
             "\n")
      if (is.euclid(x.dist2))
        if (messages)
          cat("Species x species distance matrix was not Euclidean. 'sqrt' correction was applied.",
              "\n")
    }
    if (corr == "none") {
      x.dist2 <- quasieuclid(x.dist)
      if (messages)
        cat("Species x species distance was not Euclidean, but no correction was applied. Only the PCoA axes with positive eigenvalues were kept.",
            "\n")
    }
  }
  x.pco <- dudi.pco(x.dist2, scannf = FALSE, full = TRUE)
  traits <- round(x.pco$li, .Machine$double.exponent)
  nb.sp <- numeric(c)
  for (i in 1:c) {
    sp.pres <- which(a[i, ] > 0)
    traits.sp.pres <- traits[sp.pres, , drop = F]
    traits.sp.pres[traits.sp.pres != 0 & abs(traits.sp.pres) <
                     tol] <- 0
    nb.sp[i] <- nrow(unique(traits.sp.pres))
  }
  names(nb.sp) <- row.names(a)
  min.nb.sp <- min(nb.sp)
  if (min.nb.sp < 3)
    if (messages)
      cat("FEVe: Could not be calculated for communities with <3 functionally singular species.",
          "\n")
  if (min.nb.sp < 2)
    if (messages)
      cat("FDis: Equals 0 in communities with only one functionally singular species.",
          "\n")
  if (calc.FRic) {
    x.class2 <- sapply(x, data.class)
    if (all(x.class2 == "factor" | x.class2 == "ordered")) {
      if (length(x.class2) == 1 & x.class2[1] == "ordered") {
        traits.FRic1 <- rank(x[, 1])
        names(traits.FRic1) <- x.rn
        traits.FRic <- data.frame(traits.FRic1)
        qual.FRic = 1
        if (messages)
          cat("FRic: Only one ordinal trait present in 'x'. FRic was measured as the range of the ranks, NOT as the convex hull volume.",
              "\n")
        if (calc.FDiv) {
          calc.FDiv <- FALSE
          if (messages)
            cat("FDiv: Cannot be computed when 'x' is a single ordinal trait.",
                "\n")
        }
        if (stand.FRic) {
          traits.range <- range(traits.FRic[, 1])
          FRic.all <- traits.range[2] - traits.range[1]
        }
      }
      else {
        traits.FRic <- x
        qual.FRic = 1
        if (messages)
          cat("FRic: Only categorical and/or ordinal trait(s) present in 'x'. FRic was measured as the number of unique trait combinations, NOT as the convex hull volume.",
              "\n")
        if (stand.FRic)
          FRic.all <- nrow((unique(traits.FRic)))
        if (calc.FDiv) {
          calc.FDiv <- FALSE
          if (messages)
            cat("FDiv: Cannot be computed when only categorical and/or ordinal trait(s) present in 'x'.",
                "\n")
        }
      }
    }
    else {
      if (x.pco$nf == 1) {
        traits.FRic <- x.pco$li
        qual.FRic = 1
        if (messages)
          cat("FRic: Only one continuous trait or dimension in 'x'. FRic was measured as the range, NOT as the convex hull volume.",
              "\n")
        if (calc.FDiv) {
          calc.FDiv <- FALSE
          if (messages)
            cat("FDiv: Cannot not be computed when 'x' contains one single continuous trait or dimension.",
                "\n")
        }
        if (stand.FRic) {
          traits.range <- range(traits.FRic[, 1])
          FRic.all <- traits.range[2] - traits.range[1]
        }
      }
      if (x.pco$nf > 1) {
        warning <- FALSE
        m.max <- min.nb.sp - 1
        if (m == "min") {
          warning <- TRUE
          if (min.nb.sp < 4) {
            nb.sp2 <- nb.sp[nb.sp > 3]
            m.min <- floor(log2(min(nb.sp2)))
            if (messages)
              cat("FRic: To respect s >= 2^t, FRic could not be calculated for communities with <4 functionally singular species.",
                  "\n")
          }
          else m.min <- floor(log2(min.nb.sp))
        }
        else {
          if (min.nb.sp < 3) {
            nb.sp2 <- nb.sp[nb.sp > 2]
            m.max <- min(nb.sp2) - 1
            if (messages)
              cat("FRic: To respect s > t, FRic could not be calculated for communities with <3 functionally singular species.",
                  "\n")
          }
          else m.max <- m.max
        }
        if (is.numeric(m) & m <= 1)
          stop("When 'm' is an integer, it must be >1.",
               "\n")
        if (is.numeric(m) & m > m.max)
          m <- m.max
        if (m == "min")
          m <- m.min
        if (m == "max")
          m <- m.max
        if (!is.numeric(m) & m != "min" & m != "max")
          stop("'m' must be an integer >1, 'min', or 'max'.",
               "\n")
        if (m < x.pco$nf) {
          traits.FRic <- x.pco$li[, 1:m]
          if (x.pco$nf - m == 1)
            if (messages)
              cat("FRic: Dimensionality reduction was required. The last PCoA axis (out of",
                  x.pco$nf, "in total) was removed.",
                  "\n")
          if (x.pco$nf - m > 1)
            if (messages)
              cat("FRic: Dimensionality reduction was required. The last",
                  x.pco$nf - m, "PCoA axes (out of", x.pco$nf,
                  "in total) were removed.", "\n")
          if (is.euclid(x.dist)) {
            qual.FRic <- sum(x.pco$eig[1:m])/sum(x.pco$eig)
            if (messages)
              cat("FRic: Quality of the reduced-space representation =",
                  qual.FRic, "\n")
          }
          if (!is.euclid(x.dist) & corr != "none") {
            qual.FRic <- sum(x.pco$eig[1:m])/sum(x.pco$eig)
            if (messages)
              cat("FRic: Quality of the reduced-space representation (based on corrected distance matrix) =",
                  qual.FRic, "\n")
          }
          if (!is.euclid(x.dist) & corr == "none") {
            delta <- -0.5 * bicenter.wt(x.dist * x.dist)
            lambda <- eigen(delta, symmetric = TRUE,
                            only.values = TRUE)$values
            sum.m <- sum(lambda[1:m])
            sum.n <- sum(lambda)
            lambda.neg <- c(lambda[lambda < 0])
            max.neg <- abs(min(lambda.neg))
            qual.FRic <- (sum.m + (length(lambda[1:m]) *
                                     max.neg))/(sum.n + ((length(lambda) -
                                                            1) * max.neg))
            if (messages)
              cat("FRic: Quality of the reduced-space representation (taking into account the negative eigenvalues) =",
                  qual.FRic, "\n")
          }
        }
        if (m >= x.pco$nf) {
          qual.FRic = 1
          traits.FRic <- x.pco$li
          if (x.pco$nf == 2)
            if (messages)
              cat("FRic: No dimensionality reduction was required. The 2 PCoA axes were kept as 'traits'.",
                  "\n")
          if (x.pco$nf > 2)
            if (messages)
              cat("FRic: No dimensionality reduction was required. All",
                  x.pco$nf, "PCoA axes were kept as 'traits'.",
                  "\n")
        }
        if (stand.FRic) {
          hull.all <- convhulln(traits.FRic, c("FA", "QJ"))
          FRic.all <- hull.all$vol
        }
      }
    }
  }
  if (!calc.FRic & calc.FDiv)
    cat("FDiv: Cannot be computed when 'calc.FRic' is FALSE.",
        "\n")
  if (calc.FRic & calc.FDiv)
    if (min.nb.sp < 3)
      if (messages)
        cat("FDiv: Could not be calculated for communities with <3 functionally singular species.",
            "\n")
  if (calc.FGR) {
    if (clust.type == "kmeans") {
      tr.clust <- cascadeKM(traits, km.inf.gr, km.sup.gr,
                            km.iter, km.crit)
      cat("FGR: Summary of kmeans clustering\n")
      cat("\nPartition\n")
      print(tr.clust$partition)
      cat("\nResults\n")
      print(tr.clust$results)
      cat("\nSize\n")
      print(tr.clust$size)
      plot(tr.clust)
      part.names <- colnames(tr.clust$partition)
      part.names <- as.numeric(substr(part.names, 1, 1))
      cat("\nFGR: How many groups?", "\n")
      cut.g <- toupper(scan(file = "", what = "character",
                            nlines = 1, quiet = T))
      cut.gr <- as.integer(cut.g)
      if (cut.gr < km.inf.gr | cut.gr > km.sup.gr)
        stop("You must type an integer between 'km.ing.gr' and 'km.sup.gr'.",
             "\n")
      spfgr.all <- tr.clust$partition[, part.names ==
                                        cut.gr]
      names(spfgr.all) <- x.rn
    }
    else {
      tr.clust <- hclust(x.dist, method = clust.type)
      plot(tr.clust, main = "Cluster dengrogram of species based on functional traits")
      cat("FGR: Do you want to cut the dendrogram from height or from the number of groups? Type 'h' for height, 'g' for groups.",
          "\n")
      cut <- toupper(scan(file = "", what = "character",
                          nlines = 1, quiet = T))
      if (cut == "H") {
        cat("FGR: At what height do you want the dendrogram to be cut?",
            "\n")
        cut.d <- toupper(scan(file = "", what = "character",
                              nlines = 1, quiet = T))
        cut.dist <- as.numeric(cut.d)
        spfgr.all <- cutree(tr.clust, h = cut.dist)
      }
      if (cut == "G") {
        cat("FGR: How many groups?", "\n")
        cut.g <- toupper(scan(file = "", what = "character",
                              nlines = 1, quiet = T))
        cut.gr <- as.integer(cut.g)
        spfgr.all <- cutree(tr.clust, k = cut.gr)
      }
      if (cut != "H" & cut != "G")
        stop("You must type 'h' or 'g'", "\n")
    }
    a.t <- t(a)
    by.gr <- list(spfgr.all)
    gr.abun <- aggregate(a.t, by.gr, sum)
    lab <- paste("group", gr.abun[, 1], sep = "")
    gr.abun <- data.frame(t(gr.abun[, -1]))
    colnames(gr.abun) <- lab
    rownames(gr.abun) <- rownames(a)
  }
  if (is.matrix(x) | is.data.frame(x) & calc.CWM) {
    CWM <- functcomp(x, a, CWM.type = CWM.type)
  }
  if (calc.CWM & class(x)[1] == "dist" | class(x)[1] == "dissimilarity")
    if (messages)
      cat("CWM: When 'x' is a distance matrix, CWM cannot be calculated.",
          "\n")
  divc <- function(df, dis = NULL, scale = FALSE) {
    if (!inherits(df, "data.frame"))
      stop("Non convenient df")
    if (any(df < 0))
      stop("Negative value in df")
    if (!is.null(dis)) {
      if (!inherits(dis, "dist"))
        stop("Object of class 'dist' expected for distance")
      dis <- as.matrix(dis)
      if (nrow(df) != nrow(dis))
        stop("Non convenient df")
      dis <- as.dist(dis)
    }
    if (is.null(dis))
      dis <- as.dist((matrix(1, nrow(df), nrow(df)) -
                        diag(rep(1, nrow(df)))) * sqrt(2))
    div <- as.data.frame(rep(0, ncol(df)))
    names(div) <- "diversity"
    rownames(div) <- names(df)
    for (i in 1:ncol(df)) {
      if (sum(df[, i]) < 1e-16)
        div[i, ] <- 0
      else div[i, ] <- (t(df[, i]) %*% (as.matrix(dis)^2) %*%
                          df[, i])/2/(sum(df[, i])^2)
    }
    if (scale == TRUE) {
      divmax <- divcmax(dis)$value
      div <- div/divmax
    }
    return(div)
  }
  RaoQ <- divc(data.frame(t(a)), x.dist, scale = scale.RaoQ)
  RaoQ <- RaoQ[, 1]
  names(RaoQ) <- rownames(a)
  disp <- fdisp(x.dist, a)
  FDis <- disp$FDis
  nbsp <- rep(NA, c)
  names(nbsp) <- row.names(a)
  FRic <- rep(NA, c)
  names(FRic) <- row.names(a)
  FEve <- rep(NA, c)
  names(FEve) <- row.names(a)
  FGR <- rep(NA, c)
  names(FGR) <- row.names(a)
  FDiv <- rep(NA, c)
  names(FDiv) <- row.names(a)
  for (i in 1:c) {
    sppres <- which(a[i, ] > 0)
    S <- length(sppres)
    nbsp[i] <- S
    tr <- data.frame(traits[sppres, ])
    if (calc.FRic)
      tr.FRic <- data.frame(traits.FRic[sppres, ])
    ab <- as.matrix(a[i, sppres])
    abundrel <- ab/sum(ab)
    if (calc.FRic) {
      if (all(x.class2 == "factor" | x.class2 == "ordered")) {
        if (length(x.class2) == 1 & x.class2[1] == "ordered") {
          tr.range <- range(tr.FRic[, 1])
          t.range <- tr.range[2] - tr.range[1]
          if (!stand.FRic)
            FRic[i] <- t.range
          if (stand.FRic)
            FRic[i] <- t.range/FRic.all
        }
        else {
          if (!stand.FRic)
            FRic[i] <- nrow((unique(tr.FRic)))
          if (stand.FRic)
            FRic[i] <- nrow((unique(tr.FRic)))/FRic.all
        }
      }
      else {
        if (dim(tr.FRic)[2] > 1 & nb.sp[i] >= 3) {
          if (warning)
            thresh <- 4
          if (!warning)
            thresh <- 3
          if (nb.sp[i] >= thresh) {
            convhull <- convhulln(tr.FRic, c("FA", "QJ"))
            if (!stand.FRic)
              FRic[i] <- convhull$vol
            if (stand.FRic)
              FRic[i] <- convhull$vol/FRic.all
          }
          else {
          }
        }
        if (dim(tr.FRic)[2] == 1) {
          tr.range <- range(tr.FRic[, 1])
          t.range <- tr.range[2] - tr.range[1]
          if (!stand.FRic)
            FRic[i] <- t.range
          if (stand.FRic)
            FRic[i] <- t.range/FRic.all
        }
      }
    }
    if (nb.sp[i] >= 3) {
      tr.dist <- dist(tr)
      linkmst <- mst(tr.dist)
      mstvect <- as.dist(linkmst)
      abund2 <- matrix(0, nrow = S, ncol = S)
      for (q in 1:S) for (r in 1:S) abund2[q, r] <- abundrel[q] +
        abundrel[r]
      abund2vect <- as.dist(abund2)
      EW <- rep(0, S - 1)
      flag <- 1
      for (m in 1:((S - 1) * S/2)) {
        if (mstvect[m] != 0) {
          EW[flag] <- tr.dist[m]/(abund2vect[m])
          flag <- flag + 1
        }
      }
      minPEW <- rep(0, S - 1)
      OdSmO <- 1/(S - 1)
      for (l in 1:(S - 1)) minPEW[l] <- min((EW[l]/sum(EW)),
                                            OdSmO)
      FEve[i] <- ((sum(minPEW)) - OdSmO)/(1 - OdSmO)
    }
    if (calc.FDiv & calc.FRic) {
      if (any(x.class2 == "numeric") & dim(tr.FRic)[2] >
          1 & nb.sp[i] >= 3) {
        vert0 <- convhulln(tr.FRic, c("Fx TO 'vert.txt'", "QJ"))
        vert1 <- scan("vert.txt", quiet = T)
        vert2 <- vert1 + 1
        vertices <- vert2[-1]
        trvertices <- tr.FRic[vertices, ]
        baryv <- apply(trvertices, 2, mean)
        distbaryv <- rep(0, S)
        for (j in 1:S) distbaryv[j] <- (sum((tr.FRic[j,
                                                     ] - baryv)^2))^0.5
        meandB <- mean(distbaryv)
        devdB <- distbaryv - meandB
        abdev2 <- abundrel * devdB
        ababsdev2 <- abundrel * abs(devdB)
        FDiv[i] <- (sum(abdev2) + meandB)/(sum(ababsdev2) +
                                             meandB)
      }
    }
    if (calc.FGR)
      FGR[i] <- length(unique(spfgr.all[sppres]))
  }
  res <- list()
  res$nbsp <- nbsp
  res$sing.sp <- nb.sp
  if (calc.FRic)
    res$FRic <- FRic
  if (calc.FRic)
    res$qual.FRic <- qual.FRic
  res$FEve <- FEve
  if (calc.FDiv)
    res$FDiv <- FDiv
  res$FDis <- FDis
  res$RaoQ <- RaoQ
  if (calc.FGR) {
    res$FGR <- FGR
    res$spfgr <- spfgr.all
    res$gr.abun <- gr.abun
  }
  if (is.matrix(x) | is.data.frame(x) & calc.CWM)
    res$CWM <- CWM
  if (print.pco) {
    res$x.values <- x.pco$eig
    res$x.axes <- x.pco$li
  }
  invisible(res)
}
karinorman/biodivTS documentation built on Dec. 23, 2024, 10 p.m.