R/Tsallis.R

In entropart: Entropy Partitioning to Measure Diversity

Tsallis <-
function(NorP, q = 1, ...) 
{
  UseMethod("Tsallis")
}


Tsallis.ProbaVector <-
function(NorP, q = 1, ..., CheckArguments = TRUE, Ps = NULL) 
{
  if (missing(NorP)){
    if (!missing(Ps)) {
      NorP <- Ps
    } else {
      stop("An argument NorP or Ps must be provided.")
    }
  }
  if (CheckArguments)
    CheckentropartArguments()
  
  Ps <- NorP
  dataTsallis <- -Ps^q * lnq(Ps, q)
  # Eliminate unobserved species
  dataTsallis[Ps == 0] <- 0
  
  entropy <- sum(dataTsallis)
  names(entropy) <- "None"
  return (entropy)
}


Tsallis.AbdVector <-
function(NorP, q = 1, Correction = "Best", Level = NULL, PCorrection="Chao2015", Unveiling="geom", RCorrection="Rarefy", ..., CheckArguments = TRUE, Ns = NULL) 
{
  if (missing(NorP)){
    if (!missing(Ns)) {
      NorP <- Ns
    } else {
      stop("An argument NorP or Ns must be provided.")
    }
  }
  if (is.null(Level)) {
    return (bcTsallis(Ns=NorP, q=q, Correction=Correction, CheckArguments=CheckArguments))
  } else {
    return (Tsallis.numeric(NorP, q=q, Correction=Correction, Level=Level, PCorrection=PCorrection, Unveiling=Unveiling, RCorrection=RCorrection, CheckArguments=CheckArguments))
  }
}


Tsallis.integer <-
function(NorP, q = 1, Correction = "Best", Level = NULL, PCorrection="Chao2015", Unveiling="geom", RCorrection="Rarefy", ..., CheckArguments = TRUE, Ns = NULL)
{
  if (missing(NorP)){
    if (!missing(Ns)) {
      NorP <- Ns
    } else {
      stop("An argument NorP or Ns must be provided.")
    }
  }
  if (is.null(Level)) {
    return(bcTsallis(Ns=NorP, q=q, Correction=Correction, CheckArguments=CheckArguments))
  } else {
    return (Tsallis.numeric(NorP, q=q, Correction=Correction, Level=Level, PCorrection=PCorrection, Unveiling=Unveiling, RCorrection=RCorrection, CheckArguments=CheckArguments))
  }
}


Tsallis.numeric <-
function(NorP, q = 1, Correction = "Best", Level = NULL, PCorrection = "Chao2015", Unveiling = "geom", RCorrection = "Rarefy", ..., CheckArguments = TRUE, Ps = NULL, Ns = NULL) 
{
  if (missing(NorP)){
    if (!missing(Ps)) {
      NorP <- Ps
    } else {
      if (!missing(Ns)) {
        NorP <- Ns
      } else {
        stop("An argument NorP or Ps or Ns must be provided.")
      }
    }
  }
  if (CheckArguments)
    CheckentropartArguments()
  
  if (abs(sum(NorP) - 1) < length(NorP)*.Machine$double.eps) {
    # Probabilities sum to 1, allowing rounding error
    return(Tsallis.ProbaVector(NorP, q=q, CheckArguments=FALSE))
  } else {
    # Abundances
    if (is.null(Level)) {
      return(bcTsallis(Ns=NorP, q=q, Correction=Correction, CheckArguments=FALSE))
    }
    # Eliminate 0
    NorP <- NorP[NorP > 0]
    N <- sum(NorP)
    if (Level == sum(NorP)) {
      # No interpolation/extrapolation needed: estimate with no correction
      return(Tsallis.ProbaVector(NorP/N, q=q, CheckArguments=FALSE))
    }
    # Interpolation or extrapolation
    if (q==0) {
      # Richness-1. Same result as general formula but faster
      return(Richness.numeric(NorP, Correction=Correction, Level=Level, PCorrection=PCorrection, Unveiling=Unveiling, RCorrection=RCorrection, CheckArguments=FALSE) - 1)
    } 
    if (q==1) {
      # Shannon. General formula is not defined at q=1
      return(Shannon.numeric(NorP, Correction=Correction, Level=Level, PCorrection=PCorrection, Unveiling=Unveiling, RCorrection=RCorrection, CheckArguments=CheckArguments))
    } 
    if (q==2) {
      # Simpson. No optional estimator: the uniased estimator is used.
      return(Simpson.numeric(NorP, Level=Level, CheckArguments=FALSE))
    }
    # non integer q
    # If Level is coverage, get size
    if (Level < 1) 
      Level <- Coverage2Size(NorP, SampleCoverage=Level, CheckArguments=FALSE)
    if (Level <= N) {
      # Obtain Abundance Frequence Count
      afc <- AbdFreqCount(NorP, Level=Level, CheckArguments=FALSE)
      # Calculate entropy (Chao et al., 2014, eq. 6)
      entropy <- (sum(((seq_len(Level))/Level)^q * afc[, 2]) - 1)/(1-q)
      names(entropy) <- attr(afc, "Estimator")
      return (entropy)
    } else {
      # Extrapolation.
      if (length(NorP) == 1) {
        # Single species: general formula won't work: log(1-PsU)
        entropy <- 0
        names(entropy) <- "Single Species"
      } else {
        # Unveil the full distribution that rarefies to the observed entropy
        PsU <- as.ProbaVector(NorP, Correction=PCorrection, Unveiling=Unveiling, RCorrection=RCorrection, q=q, CheckArguments=FALSE)
        # AbdFreqCount at Level (Chao et al., 2014, eq. 5)
        Slevel <- vapply(seq_len(Level), function(nu) sum(exp(lchoose(Level, nu) + nu*log(PsU) + (Level-nu)*log(1-PsU))), FUN.VALUE=0.0)
        # Estimate entropy (Chao et al., 2014, eq. 6)
        entropy <- (sum((seq_len(Level)/Level)^q * Slevel) - 1) / (1-q)
        names(entropy) <- RCorrection
      }
      return (entropy)
    }
  }
}


bcTsallis <-
function(Ns, q = 1, Correction = "Best", SampleCoverage = NULL, CheckArguments = TRUE) 
{
  if (CheckArguments)
    CheckentropartArguments()
  
  if (Correction == "Best") Correction <- "UnveilJ"
  
  # Eliminate 0
  Ns <- Ns[Ns > 0]
  N <- sum(Ns)
  # Exit if Ns contains no or a single species
  if (length(Ns) < 2) {
	  if (length(Ns) == 0) {
	    entropy <- NA
	    names(entropy) <- "No Species"
	    return (entropy)
	  } else {
	    entropy <- 0
	    names(entropy) <- "Single Species"
	    return (entropy)
	  }
  } 
  
  # Meta-Community estimation (Ns may not be integers, SampleCoverage is given)
  # SampleCoverage is between 0 and 1 (by CheckArguments), sum(Ns) must be an integer.
  # Correction may be ChaoShen or Marcon (max(ChoShen, Grassberger))
  if (!is.null(SampleCoverage) & is.IntValues(N) & (Correction == "ChaoShen" | Correction == "Marcon")) {
    CPs <- SampleCoverage*Ns/N
    ChaoShen <- -sum(CPs^q * lnq(CPs, q) /(1 - (1-CPs)^N))
    if (Correction == "Marcon") {
      # Calculate Grassberger's correction
      if (q == 1) {
        Grassberger <- sum(Ns/N*(log(N)-digamma(Ns)-(1-round(Ns)%%2*2)/(Ns+1)))
      } else {
        Grassberger <- (1-N^(-q)*sum(Enq(Ns, q)))/(q-1)
      }
    } else Grassberger <- 0
    # Take the max
    if (ChaoShen > Grassberger) {
      entropy <- ChaoShen
      names(entropy) <- "ChaoShen"
    } else {
      entropy <- Grassberger
      names(entropy) <- "Grassberger"
    }
    return (entropy)
  }
  
  # No correction
  if (Correction == "None") {
    return (Tsallis.ProbaVector(Ns/sum(Ns), q, CheckArguments = FALSE))
  } else {
    if (!is.IntValues(Ns)) {
      warning("Correction can't be applied to non-integer values.")
      # Correction <- "None"
      return (Tsallis.ProbaVector(Ns/sum(Ns), q, CheckArguments = FALSE))
    }
  }
  
  
  # Common code for ZhangGrabchak. Useless if EntropyEstimation is used.
  # if (Correction == "ZhangGrabchak" | Correction == "ChaoWangJost" | Correction == "ChaoJost") {
  #   Ps <- Ns/N
  #   V <- seq_len(N-1)
  #   # p_V_Ns is an array, containing (1 - (n_s-1)/(n-j)) for each species (lines) and all j from 1 to n-1
  #   p_V_Ns <- outer(Ns, V, function(Ns, j) 1- (Ns-1)/(N-j))
  #   # Useful values are products from j=1 to v, so prepare cumulative products
  #   p_V_Ns <- t(apply(p_V_Ns, 1, cumprod))
  #   # Sum of products weighted by w_v
  #   S_v <- function(s) {
  #     Usedv <- seq_len(N-Ns[s])
  #     return(sum(w_v[Usedv]*p_V_Ns[s, Usedv]))
  #   }
  # }
  
  # Shannon
  if (q == 1) {
    if (Correction == "Marcon") {
      ChaoShen <- bcShannon(Ns, Correction="ChaoShen", CheckArguments=FALSE)
      Grassberger <- bcShannon(Ns, Correction="Grassberger", CheckArguments=FALSE)
      if (ChaoShen > Grassberger) {
        entropy <- ChaoShen
        names(entropy) <- "ChaoShen"
      } else {
        entropy <- Grassberger
        names(entropy) <- "Grassberger"
      }
      return (entropy)
    } else {
      if (Correction == "ZhangGrabchak") {
        # Weights. Useless if EntropyEstimation is used.
        # w_v <- 1/V
        # entropy <- sum(Ps*vapply(seq_along(Ns), S_v, 0))
        # Use EntropyEstimation instead
        entropy <- EntropyEstimation::Entropy.z(Ns)
        names(entropy) <- Correction
        return (entropy)
      } else {
      return (bcShannon(Ns, Correction=Correction, CheckArguments=FALSE)) 
      }
    }
  }
  
  # Not Shannon
  if (Correction == "ZhangGrabchak" | Correction == "ChaoWangJost" | Correction == "ChaoJost") {
    # Weights. Useless here if EntropyEstimation is used, but weights are necessary for ChaoWangJost
    # i <- seq_len(N)
    # w_vi <- (i-q)/i
    # w_v <- cumprod(w_vi)
    # ZhangGrabchak <- sum(Ps*vapply(seq_along(Ns), S_v, 0))/(1-q)
    # Use EntropyEstimation instead
    if (q==0) ZhangGrabchak <- length(Ns)-1 else ZhangGrabchak <- EntropyEstimation::Tsallis.z(Ns, q)
    # ZhangGrabchak stops here, but ChaoWangJost adds an estimator of the bias
    if (Correction == "ZhangGrabchak") {
      names(ZhangGrabchak) <- Correction
      return(ZhangGrabchak)
    } 
    # Calculate abundance distribution to obtain A
    DistN <- tapply(Ns, Ns, length)
    Singletons <- DistN["1"]
    if (is.na(Singletons)) Singletons <- 0
    Doubletons <- DistN["2"]
    if (is.na(Doubletons)) Doubletons <- 0
    # Calculate A (Chao & Jost, 2015, eq. 6b)
    if (Doubletons) {
      A <- 2*Doubletons/((N-1)*Singletons+2*Doubletons)
    } else {
      if (Singletons) {
        A <- 2/((N-1)*(Singletons-1)+2)
      } else {
        A <- 1
      }
    }
    # Eq 7d in Chao & Jost (2015). Terms for r in seq_len(N-1) equal (-1)^r * w_v[r] * (A-1)^r. w_v is already available from ZhangGrabchak
    i <- seq_len(N)
    # Weights: here only if EntropyEstimation is used. Else, they have been calculated above.
    w_vi <- (i-q)/i
    w_v <- cumprod(w_vi)
    if (A == 1) {
      # The general formula of Eq 7d has a 0/0 part that must be forced to 0
      ChaoJostBias <- 0
    } else {
      Eq7dSum <- vapply(seq_len(N-1), function(r) w_v[r]*(1-A)^r, 0)
      # Calculate the estimator of the bias. Eq7dSum contains all terms of the sum except for r=0: the missing term equals 1.
      # The bias in Chao & Jost (2015) is that of the Hill number. It must be divided by 1-q to be applied to entropy.
      ChaoJostBias <- (Singletons/N*(1-A)^(1-N) * (A^(q-1) - sum(Eq7dSum) - 1))/(1-q)
    }
    entropy <- as.numeric(ZhangGrabchak + ChaoJostBias)
    names(entropy) <- Correction
    return (entropy)
  }
  if (Correction == "ChaoShen" | Correction == "GenCov" | Correction == "Marcon") {
    SampleCoverage <- Coverage(Ns, CheckArguments=FALSE)
  }
  if (Correction == "ChaoShen") {
    CPs <- SampleCoverage*Ns/N
  }
  if (Correction == "GenCov" | Correction == "Marcon") {
    CPs <- as.ProbaVector(Ns, Correction="Chao2015", CheckArguments = FALSE)
  } 
  if (Correction == "ChaoShen" | Correction == "GenCov" | Correction == "Marcon") {
    ChaoShen <- -sum(CPs^q * lnq(CPs, q) /(1 - (1-CPs)^N))
  }
  if (Correction == "ChaoShen" | Correction == "GenCov") {
    names(ChaoShen) <- Correction
    return(ChaoShen)
  }
  if (Correction == "Grassberger" | Correction == "Marcon") {
    Grassberger <- (1-N^(-q)*sum(Enq(Ns, q)))/(q-1)
  }
  if (Correction == "Grassberger") {
    names(Grassberger) <- Correction
    return(Grassberger)
  }
  if (Correction == "Marcon") {
    if (ChaoShen > Grassberger) {
      entropy <- ChaoShen
      names(entropy) <- "ChaoShen"
    } else {
      entropy <- Grassberger
      names(entropy) <- "Grassberger"
    }
    return (entropy)
  }
  if (Correction == "Holste") {
    entropy <- 1/(1-q)*(beta(length(Ns)+N, q)*sum(1/beta(Ns+1, q))-1)
    names(entropy) <- Correction
    return (entropy)
  } 
  if (Correction == "Bonachela") {
    entropy <- 1/(1-q)*(beta(2+N, q)*sum(1/beta(Ns+1, q))-1)
    names(entropy) <- Correction
    return (entropy)
  }
  if (Correction == "UnveilC") {
    TunedPs <- as.ProbaVector(Ns, Correction="Chao2015", Unveiling="geom", RCorrection = "Chao1", CheckArguments=FALSE)
  }
  if (Correction == "UnveiliC") {
    TunedPs <- as.ProbaVector(Ns, Correction="Chao2015", Unveiling="geom", RCorrection = "iChao1", CheckArguments=FALSE)
  }
  if (Correction == "UnveilJ") {
    TunedPs <- as.ProbaVector(Ns, Correction="Chao2015", Unveiling="geom", RCorrection = "Jackknife", CheckArguments=FALSE)
  }
  if (Correction == "UnveilC" | Correction == "UnveiliC" | Correction == "UnveilJ") {
    entropy <- Tsallis.ProbaVector(TunedPs, q, CheckArguments = FALSE)
    names(entropy) <- Correction
    return (entropy)
  }
  
  warning("Correction was not recognized")
  return (NA)
}