Nothing
DivEst <-
function(q = 0, MC, Biased = TRUE, Correction = "Best", Tree = NULL, Normalize = TRUE, Z = NULL, Simulations = 100, ShowProgressBar = TRUE, CheckArguments = TRUE)
{
if (CheckArguments)
CheckentropartArguments()
# Preprocess the tree
ppTree <- Preprocess.Tree(Tree)
if (Normalize) {
Height <- 1
} else {
Height <- ppTree$Height
}
# Estimation from data
RealEst <- DivPart(q, MC, Biased, Correction, ppTree, Normalize, Z, CheckArguments=FALSE)
# RedrawSpecies resamples a community according to species abundances.
RedrawSpecies<- function(SpeciesAbundances){
# Very simplified (for speed) version of rCommunity with BootstrapMethod="Marcon"
stats::rmultinom(1, sum(SpeciesAbundances), SpeciesAbundances)
}
# SimulateEntropy resamples all communities and calculates entropy
SimulateEntropy <- function(Progression) {
SimNsi <- apply(MC$Nsi, 2, RedrawSpecies)
rownames(SimNsi) <- rownames(MC$Nsi)
SimMC <- Preprocess.MC(SimNsi, MC$Wi)
NewSim <- DivPart(q, SimMC, Biased, Correction, Tree, Normalize, Z, CheckArguments=FALSE)
# update progress bar
if(ShowProgressBar & interactive())
utils::setTxtProgressBar(ProgressBar, Progression)
c(NewSim$TotalAlphaEntropy, NewSim$TotalBetaEntropy, NewSim$GammaEntropy)
}
# Simulate entropy
ProgressBar <- utils::txtProgressBar(min=0, max=Simulations)
RawSimulatedEntropy <- vapply(seq_len(Simulations), SimulateEntropy, FUN.VALUE=rep(0.0, 3))
close(ProgressBar)
# Recenter entropy
SimulatedEntropy <- RawSimulatedEntropy +
with(RealEst, c(TotalAlphaEntropy, TotalBetaEntropy, GammaEntropy)) -
apply(RawSimulatedEntropy, 1, mean)
# Transform entropy to diversity
if (q == 1) {
SimulatedDiversity <- exp(SimulatedEntropy)
} else {
SimulatedDiversity <- SimulatedEntropy
SimulatedDiversity[1,] <- expq(SimulatedEntropy[1,] / Height, q) * Height
SimulatedDiversity[3,] <- expq(SimulatedEntropy[3,] / Height, q) * Height
SimulatedDiversity[2,] <- SimulatedDiversity[3,] / SimulatedDiversity[1,]
# equals: expq(SimulatedEntropy[2,] / Height / (1 - (q-1)*SimulatedEntropy[1,]/Height), q) * Height
}
dimnames(SimulatedDiversity)[[1]] <- list("Alpha", "Beta", "Gamma")
DivEst <- RealEst
DivEst$SimulatedDiversity <- SimulatedDiversity
class(DivEst) <- c("DivEst", class(RealEst))
return (DivEst)
}
is.DivEst <-
function (x)
{
inherits(x, "DivEst")
}
plot.DivEst <-
function (x, ..., main = NULL, Which = "All",
Quantiles = c(0.025, 0.975), colValue = "red", lwdValue = 2, ltyValue = 2,
colQuantiles = "black", lwdQuantiles = 1, ltyQuantiles = 2)
{
# Save graphical parameters
op <- graphics::par(no.readonly = TRUE)
graphics::par(mfrow=c(2, 2))
if (Which == "All" | (Which == "Alpha" & is.null(main))) main <- "Alpha Diversity"
if (Which == "All" | Which == "Alpha") {
graphics::plot(as.SimTest(x$TotalAlphaDiversity, x$SimulatedDiversity["Alpha",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, lwdValue=lwdValue, ltyValue=ltyValue,
colQuantiles=colQuantiles, lwdQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
if (Which == "All" | (Which == "Beta" & is.null(main))) main <- "Beta Diversity"
if (Which == "All" | Which == "Beta") {
graphics::plot(as.SimTest(x$TotalBetaDiversity, x$SimulatedDiversity["Beta",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, lwdValue=lwdValue, ltyValue=ltyValue,
colQuantiles=colQuantiles, lwdQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
if (Which == "All" | (Which == "Gamma" & is.null(main))) main <- "Gamma Diversity"
if (Which == "All" | Which == "Gamma") {
graphics::plot(as.SimTest(x$GammaDiversity, x$SimulatedDiversity["Gamma",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, lwdValue=lwdValue, ltyValue=ltyValue,
colQuantiles=colQuantiles, lwdQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
# Legend and restore parameters
if (Which == "All") {
graphics::par(mar=c(0, 0, 0, 0))
graphics::plot(0:10, 0:10, type="n", xlab=NULL, frame.plot=FALSE, xaxt="n", yaxt="n", col.lab="white")
leg <- c("Null Distribution", "True Estimate", "95% confidence interval")
graphics::legend(2, 8, leg, col = c(1, 2, 1), lty = 1:3, merge = TRUE, cex=1)
graphics::par(op)
}
}
autoplot.DivEst <-
function (object, ..., main = NULL, Which = "All", labels = NULL, font.label = list(size=11, face="plain"),
Quantiles = c(0.025, 0.975), colValue = "red", colQuantiles = "black", ltyQuantiles = 2)
{
if (Which == "All" | (Which == "Alpha" & is.null(main))) main <- "Alpha Diversity"
if (Which == "All" | Which == "Alpha") {
AlphaPlot <- autoplot(as.SimTest(object$TotalAlphaDiversity, object$SimulatedDiversity["Alpha",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, colQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
if (Which == "Alpha")
return(AlphaPlot)
if (Which == "All" | (Which == "Beta" & is.null(main))) main <- "Beta Diversity"
if (Which == "All" | Which == "Beta") {
BetaPlot <- autoplot(as.SimTest(object$TotalBetaDiversity, object$SimulatedDiversity["Beta",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, colQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
if (Which == "Beta")
return(BetaPlot)
if (Which == "All" | (Which == "Gamma" & is.null(main))) main <- "Gamma Diversity"
if (Which == "All" | Which == "Gamma") {
GammaPlot <- autoplot(as.SimTest(object$GammaDiversity, object$SimulatedDiversity["Gamma",]), main=main,
Quantiles=Quantiles, ..., colValue=colValue, colQuantiles=colQuantiles, ltyQuantiles=ltyQuantiles)
}
if (Which == "Gamma")
return(GammaPlot)
# Which == "All": return a multiple plot
return(ggpubr::ggarrange(AlphaPlot + ggplot2::labs(x=NULL), BetaPlot + ggplot2::labs(x=NULL), GammaPlot,
ncol = 1, nrow = 3, labels=labels, font.label=font.label))
}
summary.DivEst <-
function(object, ...)
{
cat("Diversity partitioning of order", object$Order, "of MetaCommunity", object$MetaCommunity, fill=TRUE)
if (!object$Biased)
cat(" with correction:", object$Correction)
cat("\n")
if (!is.null(object$Tree)) {
cat("Phylogenetic or functional diversity was calculated according to the tree", object$Tree, "\n", fill=TRUE)
cat("Diversity is", ifelse(object$Normalized, "normalized", "not normalized"), "\n", fill=TRUE)
}
cat("Alpha diversity of communities:", "\n")
print(object$CommunityAlphaDiversities)
cat("Total alpha diversity of the communities:", "\n")
print(object$TotalAlphaDiversity)
cat("Beta diversity of the communities:", "\n")
print(object$TotalBetaDiversity)
cat("Gamma diversity of the metacommunity:", "\n")
print(object$GammaDiversity)
cat("Quantiles of simulations (alpha, beta and gamma diversity):\n")
quant <- c(0, 0.01, 0.025, 0.05, 0.1, seq(0.25, 0.75, 0.25), 0.9, 0.95, 0.975, 0.99, 1)
print(stats::quantile(object$SimulatedDiversity["Alpha", ], probs = quant))
print(stats::quantile(object$SimulatedDiversity["Beta", ], probs = quant))
print(stats::quantile(object$SimulatedDiversity["Gamma", ], probs = quant))
return(invisible(NULL))
}
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.