Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.