Nothing
DivProfile <-
function(q.seq = seq(0, 2, .1), MC, Biased = TRUE, Correction = "Best", Tree = NULL, Normalize = TRUE, Z = NULL, NumberOfSimulations = 0, Alpha = 0.05, ShowProgressBar = TRUE, CheckArguments = TRUE)
{
if (CheckArguments)
CheckentropartArguments()
# Preprocess the tree
ppTree <- Preprocess.Tree(Tree)
if (Normalize) {
Height <- 1
} else {
Height <- ppTree$Height
}
# Calculate diversity profile. Parallelize.
Diversity.seq <- simplify2array(parallel::mclapply(q.seq, DivPart, MC=MC, Biased=Biased, Correction=Correction, Tree=ppTree, Normalize=Normalize, Z=Z, CheckArguments =FALSE))
# Rearrange complex structures
Dalpha <- unlist(Diversity.seq["CommunityAlphaDiversities", ])
arrDalpha <- simplify2array(tapply(Dalpha, names(Dalpha), c))
row.names(arrDalpha) <- q.seq
Ealpha <- unlist(Diversity.seq["CommunityAlphaEntropies", ])
arrEalpha <- simplify2array(tapply(Ealpha, names(Ealpha), c))
row.names(arrEalpha) <- q.seq
# Correction
Correctionlist <- Diversity.seq["Correction", ]
if (length(unique(unlist(Correctionlist))) == 1) {
# A single correction. Keep it.
Correctionlist <- unique(unlist(Correctionlist))
} else {
# Keep the whole list. Name q values.
names(Correctionlist) <- q.seq
}
# Prepare a list of results
DivProfile <- list(MetaCommunity = unlist(Diversity.seq["MetaCommunity", 1], use.names=FALSE),
Order = unlist(Diversity.seq["Order", ], use.names=FALSE),
Biased = unlist(Diversity.seq["Biased", 1], use.names=FALSE),
Correction = Correctionlist,
Normalized = unlist(Diversity.seq["Normalized", 1], use.names=FALSE),
CommunityAlphaDiversities = arrDalpha,
CommunityAlphaEntropies = arrEalpha,
TotalAlphaDiversity = unlist(Diversity.seq["TotalAlphaDiversity", ]),
TotalBetaDiversity = unlist(Diversity.seq["TotalBetaDiversity", ]),
GammaDiversity = unlist(Diversity.seq["GammaDiversity", ]),
TotalAlphaEntropy = unlist(Diversity.seq["TotalAlphaEntropy", ]),
TotalBetaEntropy = unlist(Diversity.seq["TotalBetaEntropy", ]),
GammaEntropy = unlist(Diversity.seq["GammaEntropy", ])
)
if(!is.null(Tree))
DivProfile$Tree <- ArgumentOriginalName(Tree)
if(is.null(Z)) {
DivProfile$Method <- "HCDT"
} else {
DivProfile$Method <- "Similarity-based"
DivProfile$Z <- ArgumentOriginalName(Z)
}
class(DivProfile) <- "DivProfile"
# Confidence interval
if (NumberOfSimulations > 0) {
Q <- length(q.seq)
ProgressBar <- utils::txtProgressBar(min=-3, max=Q, initial=-3)
## Obtain simulated Metacommunities
RedrawSpecies<- function(SpeciesAbundances){
# Very simplified (for speed) version of rCommunity with BootstrapMethod="Marcon"
MetaCommunity(stats::rmultinom(NumberOfSimulations, sum(SpeciesAbundances), SpeciesAbundances))
}
if(ShowProgressBar & interactive())
utils::setTxtProgressBar(ProgressBar, -2)
# Resample each community according to species abundances
ResampledCs <- apply(MC$Nsi, 2, function(Ns) RedrawSpecies(SpeciesAbundances=Ns))
if(interactive()) utils::setTxtProgressBar(ProgressBar, -1)
# Each MC of this list is a simulation set of each original community
# Build simulated MCs by picking simulated communities
SimMC <- lapply(seq_len(NumberOfSimulations), function(i) MetaCommunity(vapply(ResampledCs, function(mc) mc$Nsi[, i], FUN.VALUE=rep(0, MC$Nspecies)), Weights=MC$Wi))
if(ShowProgressBar & interactive())
utils::setTxtProgressBar(ProgressBar, -0)
## Calculate alpha, beta and gamma of each simulated MC, at each q
# Prepare a matrix to store envelopes.
Envelopes <- matrix(nrow=12, ncol=Q)
rownames(Envelopes) <- c("TotalAlphaEntropyLow", "TotalAlphaEntropyHigh",
"TotalBetaEntropyLow", "TotalBetaEntropyHigh",
"GammaEntropyLow", "GammaEntropyHigh",
"TotalAlphaDiversityLow", "TotalAlphaDiversityHigh",
"TotalBetaDiversityLow", "TotalBetaDiversityHigh",
"GammaDiversityLow", "GammaDiversityHigh")
for (qi in seq_len(Q)) {
# Calculate diversities. Parallelize.
Diversity.qi <- simplify2array(parallel::mclapply(SimMC, function(mc) DivPart(q=q.seq[qi], MC=mc, Biased=Biased, Correction=Correction, Tree=ppTree, Normalize=Normalize, Z=Z, CheckArguments=FALSE)))
# Put alpha and gamma simulated entropies into vectors
TotalAlphaEntropy <- unlist(Diversity.qi["TotalAlphaEntropy", ])
TotalBetaEntropy <- unlist(Diversity.qi["TotalBetaEntropy", ])
GammaEntropy <- unlist(Diversity.qi["GammaEntropy", ])
# Recenter entropy
TotalAlphaEntropy <- TotalAlphaEntropy + DivProfile$TotalAlphaEntropy[qi] - mean(TotalAlphaEntropy, na.rm=TRUE)
TotalBetaEntropy <- TotalBetaEntropy + DivProfile$TotalBetaEntropy[qi] - mean(TotalBetaEntropy, na.rm=TRUE)
GammaEntropy <- GammaEntropy + DivProfile$GammaEntropy[qi] - mean(GammaEntropy, na.rm=TRUE)
TotalAlphaDiversity <- expq(TotalAlphaEntropy, q.seq[qi])
GammaDiversity <- expq(GammaEntropy, q.seq[qi])
TotalBetaDiversity <- GammaDiversity / TotalAlphaDiversity
# Quantiles
Envelopes["TotalAlphaEntropyLow", qi] <- stats::quantile(TotalAlphaEntropy, probs=Alpha, na.rm=TRUE)
Envelopes["TotalAlphaEntropyHigh", qi] <- stats::quantile(TotalAlphaEntropy, probs=1-Alpha, na.rm=TRUE)
Envelopes["TotalBetaEntropyLow", qi] <- stats::quantile(TotalBetaEntropy, probs=Alpha, na.rm=TRUE)
Envelopes["TotalBetaEntropyHigh", qi] <- stats::quantile(TotalBetaEntropy, probs=1-Alpha, na.rm=TRUE)
Envelopes["GammaEntropyLow", qi] <- stats::quantile(GammaEntropy, probs=Alpha, na.rm=TRUE)
Envelopes["GammaEntropyHigh", qi] <- stats::quantile(GammaEntropy, probs=1-Alpha, na.rm=TRUE)
Envelopes["TotalAlphaDiversityLow", qi] <- stats::quantile(TotalAlphaDiversity, probs=Alpha, na.rm=TRUE)
Envelopes["TotalAlphaDiversityHigh", qi] <- stats::quantile(TotalAlphaDiversity, probs=1-Alpha, na.rm=TRUE)
Envelopes["TotalBetaDiversityLow", qi] <- stats::quantile(TotalBetaDiversity, probs=Alpha, na.rm=TRUE)
Envelopes["TotalBetaDiversityHigh", qi] <- stats::quantile(TotalBetaDiversity, probs=1-Alpha, na.rm=TRUE)
Envelopes["GammaDiversityLow", qi] <- stats::quantile(GammaDiversity, probs=Alpha, na.rm=TRUE)
Envelopes["GammaDiversityHigh", qi] <- stats::quantile(GammaDiversity, probs=1-Alpha, na.rm=TRUE)
# Progressbar
if(ShowProgressBar & interactive())
utils::setTxtProgressBar(ProgressBar, qi)
}
close(ProgressBar)
# Integrate the envelopes into the object
DivProfile$TotalAlphaEntropyLow <- Envelopes["TotalAlphaEntropyLow", ]
DivProfile$TotalAlphaEntropyHigh <- Envelopes["TotalAlphaEntropyHigh", ]
DivProfile$TotalBetaEntropyLow <- Envelopes["TotalBetaEntropyLow", ]
DivProfile$TotalBetaEntropyHigh <- Envelopes["TotalBetaEntropyHigh", ]
DivProfile$GammaEntropyLow <- Envelopes["GammaEntropyLow", ]
DivProfile$GammaEntropyHigh <- Envelopes["GammaEntropyHigh", ]
DivProfile$TotalAlphaDiversityLow <- Envelopes["TotalAlphaDiversityLow", ]
DivProfile$TotalAlphaDiversityHigh <- Envelopes["TotalAlphaDiversityHigh", ]
DivProfile$TotalBetaDiversityLow <- Envelopes["TotalBetaDiversityLow", ]
DivProfile$TotalBetaDiversityHigh <- Envelopes["TotalBetaDiversityHigh", ]
DivProfile$GammaDiversityLow <- Envelopes["GammaDiversityLow", ]
DivProfile$GammaDiversityHigh <- Envelopes["GammaDiversityHigh", ]
}
return (DivProfile)
}
is.DivProfile <-
function (x)
{
inherits(x, "DivProfile")
}
plot.DivProfile <-
function (x, ..., main = NULL, xlab = "Order of Diversity", ylab = NULL, Which = "All",
LineWidth = 2, ShadeColor = "grey75", BorderColor = "red")
{
# Save graphical parameters
if (Which == "All") {
op <- graphics::par(no.readonly = TRUE)
graphics::par(mfrow=c(2, 2))
}
if (Which == "All" | (Which == "Alpha" & is.null(main))) main <- "Total Alpha Diversity"
if (Which == "All" | (Which == "Alpha" & is.null(ylab))) ylab <- expression(paste(alpha, " diversity"))
if (Which == "All" | Which == "Alpha") {
graphics::plot(y=x$TotalAlphaDiversity, x=x$Order, type="n", main=main, xlab=xlab, ylab=ylab, ylim=c(min(x$TotalAlphaDiversity, x$TotalAlphaDiversityLow, na.rm=TRUE), max(x$TotalAlphaDiversity, x$TotalAlphaDiversityHigh, na.rm=TRUE)), ...)
if (!(is.null(x$TotalAlphaDiversityHigh) | is.null(x$TotalAlphaDiversityLow))) {
# Shaded polygon (adapted from Didzis Elferts,
# http://stackoverflow.com/questions/14069629/plotting-confidence-intervals)
graphics::polygon(c(x$Order, rev(x$Order)), c(pmax(x$TotalAlphaDiversityLow, graphics::par('usr')[3]), pmin(rev(x$TotalAlphaDiversityHigh), graphics::par('usr')[4])), col = ShadeColor, border = FALSE)
# Add red lines on borders of polygon
graphics::lines(x$Order, x$TotalAlphaDiversityHigh, col=BorderColor, lty=2)
graphics::lines(x$Order, x$TotalAlphaDiversityLow, col=BorderColor, lty=2)
}
graphics::lines(y=x$TotalAlphaDiversity, x=x$Order, main=main, xlab=xlab, ylab=ylab, lwd=LineWidth, ...)
}
if (Which == "All" | (Which == "Communities" & is.null(main))) main <- "Alpha Diversity of Communities"
if (Which == "All" | (Which == "Communities" & is.null(ylab))) ylab <- expression(paste(alpha, " diversity"))
if (Which == "All" | Which == "Communities") {
Palette <- grDevices::palette(grDevices::rainbow(ncol(x$CommunityAlphaDiversities)))
graphics::plot(x$CommunityAlphaDiversities[, 1] ~ x$Order, type="n", xlim=c(min(x$Order), max(x$Order)), ylim=c(min(x$CommunityAlphaDiversities), max(x$CommunityAlphaDiversities)), main=main, xlab=xlab, ylab=ylab, ...)
for (Community in seq_len(ncol(x$CommunityAlphaDiversities))) {
graphics::lines(x=x$Order, y=x$CommunityAlphaDiversities[, Community], lty=Community, col=Palette[Community])
}
if (Which == "Communities") {
graphics::legend("topright", colnames(x$CommunityAlphaDiversities), lty=seq_len(ncol(x$CommunityAlphaDiversities)), col=Palette, inset=0.01)
}
}
if (Which == "All" | (Which == "Beta" & is.null(main))) main <- "Beta Diversity"
if (Which == "All" | (Which == "Beta" & is.null(ylab))) ylab <- expression(paste(beta, " diversity"))
if (Which == "All" | Which == "Beta") {
graphics::plot(y=x$TotalBetaDiversity, x=x$Order, type="n", main=main, xlab=xlab, ylab=ylab, ylim=c(min(x$TotalBetaDiversity, x$TotalBetaDiversityLow, na.rm=TRUE), max(x$TotalBetaDiversity, x$TotalBetaDiversityHigh, na.rm=TRUE)), ...)
if (!(is.null(x$TotalBetaDiversityHigh) | is.null(x$TotalBetaDiversityLow))) {
# Shaded polygon (adapted from Didzis Elferts,
# http://stackoverflow.com/questions/14069629/plotting-confidence-intervals)
graphics::polygon(c(x$Order, rev(x$Order)), c(pmax(x$TotalBetaDiversityLow, graphics::par('usr')[3]), pmin(rev(x$TotalBetaDiversityHigh), graphics::par('usr')[4])), col = ShadeColor, border = FALSE)
# Add red lines on borders of polygon
graphics::lines(x$Order, x$TotalBetaDiversityHigh, col=BorderColor, lty=2)
graphics::lines(x$Order, x$TotalBetaDiversityLow, col=BorderColor, lty=2)
}
graphics::lines(y=x$TotalBetaDiversity, x=x$Order, main=main, xlab=xlab, ylab=ylab, lwd=LineWidth, ...)
}
if (Which == "All" | (Which == "Gamma" & is.null(main))) main <- "Gamma Diversity"
if (Which == "All" | (Which == "Gamma" & is.null(ylab))) ylab <- expression(paste(gamma, " diversity"))
if (Which == "All" | Which == "Gamma") {
graphics::plot(y=x$GammaDiversity, x=x$Order, type="n", main=main, xlab=xlab, ylab=ylab, ylim=c(min(x$GammaDiversity, x$GammaDiversityLow, na.rm=TRUE), max(x$GammaDiversity, x$GammaDiversityHigh, na.rm=TRUE)), ...)
if (!(is.null(x$GammaDiversityHigh) | is.null(x$GammaDiversityLow))) {
# Shaded polygon (adapted from Didzis Elferts,
# http://stackoverflow.com/questions/14069629/plotting-confidence-intervals)
graphics::polygon(c(x$Order, rev(x$Order)), c(pmax(x$GammaDiversityLow, graphics::par('usr')[3]), pmin(rev(x$GammaDiversityHigh), graphics::par('usr')[4])), col = ShadeColor, border = FALSE)
# Add red lines on borders of polygon
graphics::lines(x$Order, x$GammaDiversityHigh, col=BorderColor, lty=2)
graphics::lines(x$Order, x$GammaDiversityLow, col=BorderColor, lty=2)
}
graphics::lines(y=x$GammaDiversity, x=x$Order, main=main, xlab=xlab, ylab=ylab, lwd=LineWidth, ...)
}
# Restore parameters
if (Which == "All") {
graphics::par(op)
}
}
autoplot.DivProfile <-
function (object, ..., main = NULL, xlab = "Order of Diversity", ylab = NULL, Which = "All",
ShadeColor = "grey75", alpha = 0.3, BorderColor = "red", labels = NULL, font.label = list(size=11, face="plain"),
col = ggplot2::GeomLine$default_aes$colour,
lty = ggplot2::GeomLine$default_aes$linetype,
lwd = ggplot2::GeomLine$default_aes$size)
{
if (Which == "All" | (Which == "Alpha" & is.null(main))) main <- "Total Alpha Diversity"
if (Which == "All" | (Which == "Alpha" & is.null(ylab))) ylab <- expression(paste(alpha, " diversity"))
if (Which == "All" | Which == "Alpha") {
theData <- data.frame(object$Order, object$TotalAlphaDiversity)
if (!(is.null(object$TotalAlphaDiversityHigh) | is.null(object$TotalAlphaDiversityLow))) {
theData$TotalAlphaDiversityLow <- object$TotalAlphaDiversityLow
theData$TotalAlphaDiversityHigh <- object$TotalAlphaDiversityHigh
}
AlphaPlot <- ggplot2::ggplot(theData, ggplot2::aes(x=.data$object.Order, y=.data$object.TotalAlphaDiversity))
if (!(is.null(object$TotalAlphaDiversityHigh) | is.null(object$TotalAlphaDiversityLow))) {
AlphaPlot <- AlphaPlot +
ggplot2::geom_ribbon(ggplot2::aes(ymin=.data$TotalAlphaDiversityLow, ymax=.data$TotalAlphaDiversityHigh), fill=ShadeColor, alpha=alpha) +
# Add red lines on borders of polygon
ggplot2::geom_line(ggplot2::aes(y=.data$TotalAlphaDiversityLow), colour=BorderColor, linetype=2) +
ggplot2::geom_line(ggplot2::aes(y=.data$TotalAlphaDiversityHigh), colour=BorderColor, linetype=2)
}
AlphaPlot <- AlphaPlot +
ggplot2::geom_line(colour=col, linetype=lty, size=lwd) +
ggplot2::labs(title=main, x=xlab, y=ylab)
}
if (Which == "Alpha")
return(AlphaPlot)
if (Which == "All" | (Which == "Communities" & is.null(main))) main <- "Alpha Diversity of Communities"
if (Which == "All" | (Which == "Communities" & is.null(ylab))) ylab <- expression(paste(alpha, " diversity"))
if (Which == "All" | Which == "Communities") {
theData <- reshape2::melt(cbind(data.frame(object$Order), object$CommunityAlphaDiversities), id.vars="object.Order", variable.name = "Community")
CommunitiesPlot <- ggplot2::ggplot(theData, ggplot2::aes(x=.data$object.Order, y=.data$value, colour=.data$Community)) +
ggplot2::geom_line() +
ggplot2::labs(title=main, x=xlab, y=ylab)
}
if (Which == "Communities")
return(CommunitiesPlot)
if (Which == "All" | (Which == "Beta" & is.null(main))) main <- "Beta Diversity"
if (Which == "All" | (Which == "Beta" & is.null(ylab))) ylab <- expression(paste(beta, " diversity"))
if (Which == "All" | Which == "Beta") {
theData <- data.frame(object$Order, object$TotalBetaDiversity)
if (!(is.null(object$TotalBetaDiversityHigh) | is.null(object$TotalBetaDiversityLow))) {
theData$TotalBetaDiversityLow <- object$TotalBetaDiversityLow
theData$TotalBetaDiversityHigh <- object$TotalBetaDiversityHigh
}
BetaPlot <- ggplot2::ggplot(theData, ggplot2::aes(x=.data$object.Order, y=.data$object.TotalBetaDiversity))
if (!(is.null(object$TotalBetaDiversityHigh) | is.null(object$TotalBetaDiversityLow))) {
BetaPlot <- BetaPlot +
ggplot2::geom_ribbon(ggplot2::aes(ymin=.data$TotalBetaDiversityLow, ymax=.data$TotalBetaDiversityHigh), fill=ShadeColor, alpha=alpha) +
# Add red lines on borders of polygon
ggplot2::geom_line(ggplot2::aes(y=.data$TotalBetaDiversityLow), colour=BorderColor, linetype=2) +
ggplot2::geom_line(ggplot2::aes(y=.data$TotalBetaDiversityHigh), colour=BorderColor, linetype=2)
}
BetaPlot <- BetaPlot +
ggplot2::geom_line(colour=col, linetype=lty, size=lwd) +
ggplot2::labs(title=main, x=xlab, y=ylab)
}
if (Which == "Beta")
return(BetaPlot)
if (Which == "All" | (Which == "Gamma" & is.null(main))) main <- "Gamma Diversity"
if (Which == "All" | (Which == "Gamma" & is.null(ylab))) ylab <- expression(paste(gamma, " diversity"))
if (Which == "All" | Which == "Gamma") {
theData <- data.frame(object$Order, object$GammaDiversity)
if (!(is.null(object$GammaDiversityHigh) | is.null(object$GammaDiversityLow))) {
theData$GammaDiversityLow <- object$GammaDiversityLow
theData$GammaDiversityHigh <- object$GammaDiversityHigh
}
GammaPlot <- ggplot2::ggplot(theData, ggplot2::aes(x=.data$object.Order, y=.data$object.GammaDiversity))
if (!(is.null(object$GammaDiversityHigh) | is.null(object$GammaDiversityLow))) {
GammaPlot <- GammaPlot +
ggplot2::geom_ribbon(ggplot2::aes(ymin=.data$GammaDiversityLow, ymax=.data$GammaDiversityHigh), fill=ShadeColor, alpha=alpha) +
# Add red lines on borders of polygon
ggplot2::geom_line(ggplot2::aes(y=.data$GammaDiversityLow), colour=BorderColor, linetype=2) +
ggplot2::geom_line(ggplot2::aes(y=.data$GammaDiversityHigh), colour=BorderColor, linetype=2)
}
GammaPlot <- GammaPlot +
ggplot2::geom_line(colour=col, linetype=lty, size=lwd) +
ggplot2::labs(title=main, x=xlab, y=ylab)
}
if (Which == "Gamma")
return(GammaPlot)
# Which == "All": return a multiple plot
return(ggpubr::ggarrange(AlphaPlot, CommunitiesPlot, BetaPlot, GammaPlot, ncol = 2, nrow = 2, labels=labels, font.label=font.label))
}
summary.DivProfile <-
function(object, ...)
{
cat("Diversity profile 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("Diversity against its order:\n")
Values <- cbind(object$Order, object$TotalAlphaDiversity, object$TotalBetaDiversity, object$GammaDiversity)
colnames(Values) <- c("Order", "Alpha Diversity", "Beta Diversity", "Gamma Diversity")
print(Values)
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.