In allanecology/BetaDivMultifun: Testing effects of multitrophic beta-diversity and land-use intensity (LUI) on ecosystem multifunctionality
This script :
- calculates the multifunctionalities from the imputed functions dataset
- select functions for the analysis
- combine years
- calculate functional dissimilarity
cols <- RColorBrewer::brewer.pal(8, "Dark2")
Requirements :
- imputed_grlfuns : the imputed functions data
- usedforBetadivMultifun, info_data
sections_to_be_loaded <- c("imputed_functions")
source("vignettes/analysis_nonpublic.R")
Calculate mini-multifunctionalities
The mini-multifunctionalities are calculated from the imputed values.
soilCflxs
Enzymes related to soil C. All measured in the year 2011.
sce <- imputed_grlfuns[, c("beta_Glucosidase", "N_Acetyl_beta_Glucosaminidase", "Xylosidase")]
# take z-scores of the functions
imputed_grlfuns[, "soilCflxs" := multidiv(sce, sc="sd", cent=TRUE)[,1]]
imputed_grlfuns[, c("beta_Glucosidase", "N_Acetyl_beta_Glucosaminidase", "Xylosidase") := NULL]
rm(sce)
about soil processes related to N
based on clustering and biology:
- ammonium is oxidised by bacteria and archaea (AOB, AOA)
- nitrate oxidation by bacteria (NS and NB)
these two steps are separated steps. The correlation of AO and NO is based on co-occurrence.
- ammonium oxidising group : AOA, AOB, urease (because this enzyme hydrolyses urea to ammonia and CO2)
- nitrate oxidising group : NS, NB, nifH, DEA.inverted
- potential nitrification is a separate process
soilNitrateflxs
It could be that nitrospira and nitrobacter have slightly different habitat preferences --> sum their abundance as 'nitrite-oxidising functional gene abundance'.
Inversion of DEA such that functions "point to same functional direction".
nce <- imputed_grlfuns[,c("Plotn", "DEA.inverted","nifH", "nxrA_NS", "16S_NB")]
# sum abundances of nitrite oxidising functional genes
nce[, "nitOx_fga" := nxrA_NS + `16S_NB`]
nce[, c("nxrA_NS", "16S_NB") := NULL]
# mini-multifunctionality
imputed_grlfuns[, "soilNitrateflxs" := multidiv(nce[, !"Plotn", with=F], sc="sd", cent=T)[,1]]
rm(nce)
soilAmmoniaflxs
nce <- imputed_grlfuns[, c("Plotn", "Urease","amoA_AOB.2011","amoA_AOA.2011", "amoA_AOB.2016", "amoA_AOA.2016")]
# sum per year and take the mean of 2 years after
nce[, "amOX_fga2011" := amoA_AOB.2011 + amoA_AOA.2011]
nce[, "amOX_fga2016" := amoA_AOB.2016 + amoA_AOA.2016]
nce[, "amOX_fga" := apply(nce[,c("amOX_fga2011", "amOX_fga2016")],1, function(x) mean(x, na.rm = T))]
nce[, c("amOX_fga2011", "amOX_fga2016", "amoA_AOB.2011", "amoA_AOB.2016", "amoA_AOA.2011", "amoA_AOA.2016") := NULL]
imputed_grlfuns[, "soilAmmoniaflxs" := multidiv(nce[, !"Plotn", with=F], sc="sd", cent=T)[,1]]
rm(nce)
Select functions
selecting the functions which are used in analysis.
# select variables which are included in betadiv multifun analysis
usedforBetadivMultifun <- usedforBetadivMultifun[!usedforBetadivMultifun == "Plot"]
imputed_grlfuns <- imputed_grlfuns[, ..usedforBetadivMultifun]
Check Correlations
Of the non-factor columns only
# the dataset is complete now, use "complete.obs" for the cor() function now
M <- cor(imputed_grlfuns[, !colnames(imputed_grlfuns) %in% "Plotn", with=F], use="complete.obs", method = "spearman")
# imputed_grlfun_corrplot1
corrplot::corrplot(M,type="lower", addCoef.col = "black", method="color", diag=F, tl.srt=1, tl.col="black", mar=c(0,0,0,0), number.cex=0.4, order = "hclust", tl.cex = 0.6)
# saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman.pdf
# imputed_grlfun_corrplot2
corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F,
tl.cex = 0.3, order = "hclust")
# only show medium high correlations
tres <- 0.5
M[which(M < tres & M > -tres)] <- 0
corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F,
tl.cex = 0.6, order = "hclust")
# saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman_05.pdf
# only show high correlations
tres <- 0.7
M[which(M < tres & M > -tres)] <- 0
corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F,
tl.cex = 0.6, order = "hclust")
# saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman_07.pdf
rm(M)
nimputed_grlfuns <- data.table::copy(imputed_grlfuns) # save an unscaled version for EFturnover / nestedness
Dissimilarity of single functions
This part of the script prepares the single function dissimilarities for the GDM models. Each function is fitted individually.
Dissimilarity between plots is the euclidean distance - which is just a substraction in the 1D case.
For each function : build dissimilarity matrix (euclidean distance) and then center & normalise
# get a reference Table with the Plot combinations
singleEFdist <- data.table::data.table(t(combn(imputed_grlfuns$Plotn, 2)))
data.table::setnames(singleEFdist, old = c("V1", "V2"), new = c("Var1", "Var2"))
for(f in usedforBetadivMultifun[!usedforBetadivMultifun == "Plotn"]){
# diffname <- paste("sEFdist", f, sep = "_")
# print(diffname)
small <- data.frame(imputed_grlfuns[, .(Plotn, get(f))])
rownames(small) <- small$Plotn
small$Plotn <- NULL
diffsmall <- vegan::vegdist(small, method = "euclid")
diffsmall <- as.matrix(diffsmall)
diffsmall[!lower.tri(diffsmall)] <- NA
diffsmall <- reshape2::melt(diffsmall, value.name = "X")
diffsmall <- diffsmall[!is.na(diffsmall[, 3]),]
diffsmall <- data.table::data.table(diffsmall)
setnames(diffsmall, old = "X", new = f)
# assign(diffname, diffsmall)
singleEFdist <- merge(singleEFdist, diffsmall, by = c("Var1", "Var2"), all.y = T)
# saveRDS(diffsmall, file = paste(diffname, ".rds", sep = "")) # note : is not saved, only produce on demand
}
# standardisation
include <- colnames(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")]
singleEFdist[, (include) := lapply(.SD, scale01), .SDcols = include]
saveRDS(singleEFdist, "singleEFdist.rds")
Dissimilarity of functions
2 ways of dissimilarity calculation :
- EFturnover and EFnestedness : calculate presence-absence of functions by 0.5 treshold, and then calculate betadiversity turnover and nestedness components of them.
- EFdistance : calculate PCA, identify most important axis, take euclidean distance of them to get a single value for each plot.
- weighting of the PCA axis: Weighting by var explained gives the correlated functions too much weight. Equal weighting is chosen. (soil axis should not get too much weight)
EFdistance calculation DEPRECATED
From PCA of all functions, take most important axes and calculate euclidean distance of them (using equal weights).
We expect that one of the important axis - possibly the first to represent soil. As we have many soil functions compared to aboveground functions.
Selection of relevant PC axes :
- check how many correlated functions there are
- 9 correlations > 0.5 (see corrplots above)
- (number of selected axes) = (number of PCA) - (number of correlated functions)
- check if this number corresponds to around 90% variance explained.
# STANDARDISATION
nimputed_grlfuns <- data.table::copy(imputed_grlfuns) # save an unscaled version for EFturnover / nestedness
include <- names(imputed_grlfuns)[!names(imputed_grlfuns) %in% "Plotn"]
imputed_grlfuns[, (include) := lapply(.SD, scale01),.SDcols=include]
# PCA
pc_grlfuns <- imputed_grlfuns[, ..include]
pc_grlfuns <- as.data.frame(pc_grlfuns)
rownames(pc_grlfuns) <- imputed_grlfuns$Plotn
pca_grlfuns <- stats::prcomp(pc_grlfuns, scale=F)
summ <- summary(pca_grlfuns)
summ$importance[3, ] <= 0.9 # take axes 1 - 11
capture.output(summary(pca_grlfuns), file = "calc_multifun_from_imputed_pca_summary.txt")
par(mfrow = c(2, 2))
plot(pca_grlfuns)
biplot(pca_grlfuns)
# plot cumulative prop of variance explained (same as output by pca summary)
plot(cumsum(pca_grlfuns$sdev^2 / sum(pca_grlfuns$sdev^2)), pch = 16,
xlab = "PCA axes", ylab = "Cumulative Proportion of Variance", xaxt = "n", yaxt = "n")
axis(1, at = seq(1, 17, 2))
# axis(2, at = c(seq(0, 0.8, 0.2), 0.9161965, 1), labels = T)
axis(2, at = seq(0, 1, 0.2))
axis(2, at = 0.9161965, labels = "0.92", las = 2, cex = 0.04)
abline(v = 11, lty = 2)
abline(h = 0.9161965, lty = 2)
# save as : "22-05-16_calc_multifun_from_imputed_pca_output_cumvar.pdf"
dev.off()
Plotting first 5 axis
temp <- cbind(pc_grlfuns, "region" = sub("[0-9][0-9]", "", rownames(pc_grlfuns)))
pcaplots <- list(autoplot(pca_grlfuns, data = temp, colour = "region"),
autoplot(pca_grlfuns, data = temp, colour = "region", x = 3, y = 4),
autoplot(pca_grlfuns, data = temp, colour = "region", x = 5, y = 6),
autoplot(pca_grlfuns, data = temp, colour = "region", x = 7, y = 8))
plot_grid(plotlist = pcaplots, nrow = 2)
rm(temp)
# saved as : "21-03-09_calc_multifun_from_imputed_pca_autoplot_axes_1_to_7.pdf"
p <- plot(pca_grlfuns)
labelsize <- 3
p1 <- ggplot2::autoplot(pca_grlfuns, loadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize)
p2 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize,
x = 3, y = 2)
p3 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize,
x = 4, y = 3)
p4 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize,
x = 4, y = 5)
cowplot::plot_grid(plotlist = list(p1, p2, p3, p4), labels = c("A", "B", "C", "D"), align = T)
# cowplot::plot_grid(plotlist = list(p1, p2), labels = c("A", "B"), align = T)
# saved as : "calc_multifun_from_imputed_pca_axes_1_to_5.pdf"
p5 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T,
x = 5, y = 6)
p6 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T,
x = 6, y = 7)
p7 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T,
x = 7, y = 8)
p8 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black",
loadings.label.colour = "black", loadings.label.repel = T,
x = 8, y = 9)
cowplot::plot_grid(plotlist = list(p5, p6, p7, p8), labels = c("E", "F", "G", "H"), align = T)
# saved as : "23-11-09_calc_multifun_from_imputed_pca_axes_5_to_9.pdf"
rm(p); rm(p1); rm(p2); rm(p3); rm(p4); rm(p5)
axis selection :
Use axes until 90% variance explained : take axes 1 - 11
pca_grlfuns <- pca_grlfuns$x[, 1:11]
# euclidean distance
EFmaster <- vegan::vegdist(pca_grlfuns, method = "euclid")
EFmaster <- as.matrix(EFmaster)
EFmaster[!lower.tri(EFmaster)] <- NA
EFmaster <- reshape2::melt(EFmaster, value.name = "EFdistance")
EFmaster <- EFmaster[!is.na(EFmaster[, 3]),]
EFmaster <- data.table::data.table(EFmaster)
# standardise
EFmaster[, EFdistance := scale01(EFdistance)]
saveRDS(EFmaster, "EFdistance.rds")
rm(pca_grlfuns); rm(pc_grlfuns)
EFturnover and EFnestedness
Based on a 50% treshold, classify for each function and each plot presence or absence of the given function. Based on the scaled functions.
# percentage of measured functions that exceed a given treshold of their maximum observed level across all study sites.
# maximum observed level: average of the top five sites.
threshold <- 0.1
# tresholds used in Solivieres paper : 0.25,0.5,0.75,0.9
Plot for visualisation
# plot all 18 functions on one page
# x <- nimputed_grlfuns$Root.biomass
require(cowplot)
for(i in colnames(nimputed_grlfuns)[which(colnames(nimputed_grlfuns) != "Plotn")]){
x <- nimputed_grlfuns[, get(i)]
names(x) <- nimputed_grlfuns$Plotn
# calc mean of top 5 values
upper <- mean(sort(x, decreasing = T, na.last=T)[1:5])
upperq <- quantile(x, probs = threshold)
df <- cbind(len = x, deleteme = x)
p <- ggplot2::ggplot(df, ggplot2::aes(y = len, x= rownames(df))) +
ggplot2::geom_bar(stat="identity") +
ggplot2::geom_abline(intercept = upper, slope = 0, color = "lightgray") +
ggplot2::geom_abline(intercept = upper*threshold, slope = 0, color =cols[1]) + # turkis
ggplot2::geom_abline(intercept = upperq, slope = 0, color = cols[2]) + # orange
theme(axis.text.x=element_blank()) + xlab("Plots") + ylab(i)
plotname <- paste("plot", i, sep = "")
assign(plotname, p)
rm(p)
}
rm(i)
plotnames <- paste("plot", colnames(nimputed_grlfuns)[which(colnames(nimputed_grlfuns) != "Plotn")], sep = "")
plot_grid(plotlist = list(plotBiomass, plotdung.removal, plotGroundwater.recharge, plotherbivory.20172018, plotLitter.decomposition, plotN_leaching_risk2015, plotpathogen.infection, plotPhosphatase, plotPotential.nitrification),
align = T)
plot_grid(plotlist = list( plotcaterpillars.predation, plotRoot.biomass, plotRoot.decomposition, plotseed.depletion, plotsoilCflxs, plotP_leaching_risk_comb, plotsoilNitrateflxs, plotsoilAmmoniaflxs),
align = T)
# saved name : 21-03-09_calc_multifun_from_imputed_EFbeta_threshold_method_visual_tr01_part1 , change to tr = 0.7 etc, change part 1 and part2
# barplot(x, main="Root Biomass", xlab="Plots")
# abline(a=upper,b=0,col="lightgray")
# abline(a = upperq, b = 0, col = "darkgreen")
# abline(a=upper*threshold, b=0, col="orange")
rm(df); rm(plotBiomass, plotdung.removal, plotGroundwater.recharge, plotherbivory.20172018, plotLitter.decomposition, plotN_leaching_risk2015, plotpathogen.infection, plotPhosphatase, plotPotential.nitrification, plotcaterpillars.predation, plotRoot.biomass, plotRoot.decomposition, plotseed.depletion, plotsoilCflxs, plotP_leaching_risk_comb, plotsoilNitrateflxs, plotsoilAmmoniaflxs); gc()
With threshold max 5 method
include <- names(nimputed_grlfuns)[!names(nimputed_grlfuns) %in% "Plotn"]
pa_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for max 5 functions
threshold <- 0.9
Note : to make EF master for all thresholds : change manually above and run again and again the chunk below AND above.
# CONVERT TO PRESENCE ABSENCE
pa_grlfuns[, (include) := lapply(.SD, function(c) calc_presenceabsence(c, threshold = threshold, type = "max5")), .SDcols = include]
pa_grlfuns[, (include) := (.SD * 1), .SDcols = include]
setcolorder(pa_grlfuns, neworder = c("Plotn", include))
# saveRDS(pa_grlfuns, file = paste(pathtoout, paste(paste("/single_functions_passing_threshold_", threshold, sep = ""), "Rds", sep = "."), sep = ""))
# CALC BETADIVERSITY
pa_grlfuns <- data.frame(pa_grlfuns[, !"Plotn", with=F], row.names = pa_grlfuns$Plotn)
pa_grlfuns <- BetaDivMultifun::beta.pair_zerospecies(pa_grlfuns, index.family = "sorensen")
EFbeta <- pa_grlfuns$beta.sor
t <- as.matrix(EFbeta)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = paste("EFbeta", threshold, sep = "_"))
EFbeta <- t[!is.na(t[, 3]), ]
if(!exists("EFmaster")){ # if EFmaster does not exist yet, create
EFmaster <- EFbeta
} else {
EFmaster <- merge(EFmaster, EFbeta, by = c("Var1", "Var2"), all = T)
}
EFturnover <- pa_grlfuns$beta.sim
t <- as.matrix(EFturnover)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = paste("EFturnover", threshold, sep = "_"))
EFturnover <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, EFturnover, by = c("Var1", "Var2"), all = T)
EFnestedness <- pa_grlfuns$beta.sne
t <- as.matrix(EFnestedness)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = paste("EFnestedness", threshold, sep = "_"))
EFnestedness <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, EFnestedness, by = c("Var1", "Var2"), all = T)
standardise and save
# saveRDS(EFturnover, "EFturnover.rds")
# saveRDS(EFnestedness, "EFnestedness.rds")
saveRDS(EFmaster, "EFmaster_all_thresholds.rds") # if all threshods were calculated and added to EFmaster
saveRDS(EFmaster, "EFmaster.rds")
Plot components of EFbeta
# # check EFmaster
# EFcheck <- data.table::copy(data.table::data.table(EFmaster))
# EFcheck[, check := EFturnover + EFnestedness - EFbeta]
# plot(EFcheck$check, ylim = c(-1, 1))
# rm(EFcheck)
# # TRUE : EFturnover + EFnestedness - EFbeta is always 0 (or very near 0)
plotEFbeta <- reshape2::melt(EFmaster[, c("Var1", "Var2", "EFbeta_0.5", "EFturnover_0.5", "EFnestedness_0.5")], id = c("Var1", "Var2"))
p <- ggplot2::ggplot(data = plotEFbeta, aes(x = variable, y = value)) +
geom_violin() + geom_jitter(shape=".", position=position_jitter(0.2), color = "gray") +
labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "threshold = 0.5", xlab = "") +
scale_x_discrete(labels=c("beta", "turnover", "nestedness")) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())
# Plot all thresholds
plotEFbeta <- data.table(reshape2::melt(EFmaster, id = c("Var1", "Var2")))
plotEFbeta[, type := sub("_[0-9].[0-9]", "", variable)]
p1 <- ggplot2::ggplot(data = plotEFbeta[variable %in% unique(plotEFbeta$variable)[1:14]],
aes(x = variable, y = value, fill = type)) +
geom_violin() +
geom_jitter(shape=".", position=position_jitter(0.2), color = "gray", alpha = 0.2) +
theme(axis.text.x = element_text(angle = 90)) +
labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "all thresholds", xlab = "") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
# scale_fill_brewer(palette="Dark2")
scale_fill_manual(values = c("#D95F02", "#1B9E77", "#7570B3"))
p2 <- ggplot2::ggplot(data = plotEFbeta[variable %in% unique(plotEFbeta$variable)[14:28]],
aes(x = variable, y = value, fill = type)) +
geom_violin() +
geom_jitter(shape=".", position=position_jitter(0.2), color = "gray", alpha = 0.2) +
theme(axis.text.x = element_text(angle = 90)) +
labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "all thresholds", xlab = "") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
# scale_fill_brewer(palette="Dark2")
scale_fill_manual(values = c("#D95F02", "#1B9E77", "#7570B3"))
plot_grid(p1, p2, nrow = 2)
# save as : 21-10-04_calc_multifun_from_imputed_all_EFdistances_violinplot
Plot various thresholds
# only run if EFmaster contains various thresholds
temp <- c(grep("turnover|nestedness", names(EFmaster), value = T))
# temp <- grep("turnover", names(EFmaster), value = T)
# pairs(EFmaster[, ..temp], upper.panel = NULL )
EFmaster <- data.table(EFmaster)
m <- cor(EFmaster[, ..temp], method = "spearman")
pdf("calc_multifun_from_imputed_all_EFdistances_corrplot.pdf", paper = "a4r")
corrplot::corrplot(m, type = "lower", addCoef.col = "black", method = "color", diag = F, tl.srt = 50, tl.col = "black", order = "hclus", tl.cex = 0.5, cl.cex = 0.4, number.cex = 0.4)
dev.off()
# EFturnover and nestedness are additive
plot(EFmaster$EFbeta_0.5, EFmaster$EFturnover_0.5 + EFmaster$EFnestedness_0.5)
With quantile method DEPRECATED
threshold <- 0.5
include <- names(nimputed_grlfuns)[!names(imputed_grlfuns) %in% "Plotn"]
paq_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for quantile functions
# CONVERT TO PRESENCE ABSENCE
# quantile method
paq_grlfuns[, (include) := lapply(.SD, function(c) calc_presenceabsence(c, threshold = threshold, type = "quantile")), .SDcols = include]
paq_grlfuns[, (include) := (.SD * 1), .SDcols = include]
# CALC BETADIVERSITY
paq_grlfuns <- data.frame(paq_grlfuns[, !"Plotn", with = F], row.names = paq_grlfuns$Plotn)
paq_grlfuns <- betapart::beta.pair(paq_grlfuns, index.family = "sorensen")
EFturnoverQ <- paq_grlfuns$beta.sim
t <- as.matrix(EFturnoverQ)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFturnover_median")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
EFnestednessQ <- paq_grlfuns$beta.sne
t <- as.matrix(EFnestednessQ)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFnestedness_median")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
EFbetaQ <- paq_grlfuns$beta.sor
t <- as.matrix(EFbetaQ)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFbeta_median")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
saveRDS(EFmaster, "EFquantile_median.rds")
abundance betadiversity DEPRECATED
Abundance betadiversity of z-scores of functions.
abund_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for quantile functions
# CALC BETADIVERSITY
abund_grlfuns <- data.frame(abund_grlfuns[, !"Plotn", with = F], row.names = abund_grlfuns$Plotn)
# scale to be between 0 and 1000 (promille) and round to 0 decimals
# betapart does not allow negative values, and decimals do not make sense.
abund_grlfuns <- apply(abund_grlfuns, 2, function(x) round(scale01(x)*1000, digits = 0))
abund_grlfuns <- betapart::beta.pair.abund(abund_grlfuns, index.family = "bray")
temp <- abund_grlfuns$beta.bray
t <- as.matrix(temp)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFabund_bray")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
temp <- abund_grlfuns$beta.bray.bal
t <- as.matrix(temp)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFabund_bray_bal")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
temp <- abund_grlfuns$beta.bray.gra
t <- as.matrix(temp)
t[!lower.tri(t)] <- NA
t <- reshape2::melt(t, value.name = "EFabund_bray_gra")
t <- t[!is.na(t[, 3]), ]
EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2"))
saveRDS(EFmaster[, .(Var1, Var2, EFabund_bray, EFabund_bray_bal, EFabund_bray_gra)], "EFbeta_abund.rds")
rm(nimputed_grlfuns); rm(imputed_grlfuns); rm(pa_grlfuns); rm(raw_grlfuns)
rm(include)
rm(p); rm(cols); rm(tres); rm(threshold); rm(upper); rm(upperq); rm(x)
allanecology/BetaDivMultifun documentation built on Nov. 9, 2023, 8:47 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
This script : - calculates the multifunctionalities from the imputed functions dataset - select functions for the analysis - combine years - calculate functional dissimilarity
cols <- RColorBrewer::brewer.pal(8, "Dark2")
Requirements : - imputed_grlfuns : the imputed functions data - usedforBetadivMultifun, info_data
sections_to_be_loaded <- c("imputed_functions") source("vignettes/analysis_nonpublic.R")
Calculate mini-multifunctionalities
The mini-multifunctionalities are calculated from the imputed values.
soilCflxs
Enzymes related to soil C. All measured in the year 2011.
sce <- imputed_grlfuns[, c("beta_Glucosidase", "N_Acetyl_beta_Glucosaminidase", "Xylosidase")] # take z-scores of the functions imputed_grlfuns[, "soilCflxs" := multidiv(sce, sc="sd", cent=TRUE)[,1]] imputed_grlfuns[, c("beta_Glucosidase", "N_Acetyl_beta_Glucosaminidase", "Xylosidase") := NULL] rm(sce)
about soil processes related to N
based on clustering and biology: - ammonium is oxidised by bacteria and archaea (AOB, AOA) - nitrate oxidation by bacteria (NS and NB)
these two steps are separated steps. The correlation of AO and NO is based on co-occurrence.
- ammonium oxidising group : AOA, AOB, urease (because this enzyme hydrolyses urea to ammonia and CO2)
- nitrate oxidising group : NS, NB, nifH, DEA.inverted
- potential nitrification is a separate process
soilNitrateflxs
It could be that nitrospira and nitrobacter have slightly different habitat preferences --> sum their abundance as 'nitrite-oxidising functional gene abundance'.
Inversion of DEA such that functions "point to same functional direction".
nce <- imputed_grlfuns[,c("Plotn", "DEA.inverted","nifH", "nxrA_NS", "16S_NB")] # sum abundances of nitrite oxidising functional genes nce[, "nitOx_fga" := nxrA_NS + `16S_NB`] nce[, c("nxrA_NS", "16S_NB") := NULL] # mini-multifunctionality imputed_grlfuns[, "soilNitrateflxs" := multidiv(nce[, !"Plotn", with=F], sc="sd", cent=T)[,1]] rm(nce)
soilAmmoniaflxs
nce <- imputed_grlfuns[, c("Plotn", "Urease","amoA_AOB.2011","amoA_AOA.2011", "amoA_AOB.2016", "amoA_AOA.2016")] # sum per year and take the mean of 2 years after nce[, "amOX_fga2011" := amoA_AOB.2011 + amoA_AOA.2011] nce[, "amOX_fga2016" := amoA_AOB.2016 + amoA_AOA.2016] nce[, "amOX_fga" := apply(nce[,c("amOX_fga2011", "amOX_fga2016")],1, function(x) mean(x, na.rm = T))] nce[, c("amOX_fga2011", "amOX_fga2016", "amoA_AOB.2011", "amoA_AOB.2016", "amoA_AOA.2011", "amoA_AOA.2016") := NULL] imputed_grlfuns[, "soilAmmoniaflxs" := multidiv(nce[, !"Plotn", with=F], sc="sd", cent=T)[,1]] rm(nce)
Select functions
selecting the functions which are used in analysis.
# select variables which are included in betadiv multifun analysis usedforBetadivMultifun <- usedforBetadivMultifun[!usedforBetadivMultifun == "Plot"] imputed_grlfuns <- imputed_grlfuns[, ..usedforBetadivMultifun]
Check Correlations
Of the non-factor columns only
# the dataset is complete now, use "complete.obs" for the cor() function now M <- cor(imputed_grlfuns[, !colnames(imputed_grlfuns) %in% "Plotn", with=F], use="complete.obs", method = "spearman") # imputed_grlfun_corrplot1 corrplot::corrplot(M,type="lower", addCoef.col = "black", method="color", diag=F, tl.srt=1, tl.col="black", mar=c(0,0,0,0), number.cex=0.4, order = "hclust", tl.cex = 0.6) # saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman.pdf # imputed_grlfun_corrplot2 corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F, tl.cex = 0.3, order = "hclust") # only show medium high correlations tres <- 0.5 M[which(M < tres & M > -tres)] <- 0 corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F, tl.cex = 0.6, order = "hclust") # saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman_05.pdf # only show high correlations tres <- 0.7 M[which(M < tres & M > -tres)] <- 0 corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=40, diag = F, tl.cex = 0.6, order = "hclust") # saved as : 21-09-06_calc_multifun_from_imputed_corrplot_spearman_07.pdf rm(M)
nimputed_grlfuns <- data.table::copy(imputed_grlfuns) # save an unscaled version for EFturnover / nestedness
Dissimilarity of single functions
This part of the script prepares the single function dissimilarities for the GDM models. Each function is fitted individually.
Dissimilarity between plots is the euclidean distance - which is just a substraction in the 1D case.
For each function : build dissimilarity matrix (euclidean distance) and then center & normalise
# get a reference Table with the Plot combinations singleEFdist <- data.table::data.table(t(combn(imputed_grlfuns$Plotn, 2))) data.table::setnames(singleEFdist, old = c("V1", "V2"), new = c("Var1", "Var2")) for(f in usedforBetadivMultifun[!usedforBetadivMultifun == "Plotn"]){ # diffname <- paste("sEFdist", f, sep = "_") # print(diffname) small <- data.frame(imputed_grlfuns[, .(Plotn, get(f))]) rownames(small) <- small$Plotn small$Plotn <- NULL diffsmall <- vegan::vegdist(small, method = "euclid") diffsmall <- as.matrix(diffsmall) diffsmall[!lower.tri(diffsmall)] <- NA diffsmall <- reshape2::melt(diffsmall, value.name = "X") diffsmall <- diffsmall[!is.na(diffsmall[, 3]),] diffsmall <- data.table::data.table(diffsmall) setnames(diffsmall, old = "X", new = f) # assign(diffname, diffsmall) singleEFdist <- merge(singleEFdist, diffsmall, by = c("Var1", "Var2"), all.y = T) # saveRDS(diffsmall, file = paste(diffname, ".rds", sep = "")) # note : is not saved, only produce on demand } # standardisation include <- colnames(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")] singleEFdist[, (include) := lapply(.SD, scale01), .SDcols = include] saveRDS(singleEFdist, "singleEFdist.rds")
Dissimilarity of functions
2 ways of dissimilarity calculation : - EFturnover and EFnestedness : calculate presence-absence of functions by 0.5 treshold, and then calculate betadiversity turnover and nestedness components of them. - EFdistance : calculate PCA, identify most important axis, take euclidean distance of them to get a single value for each plot. - weighting of the PCA axis: Weighting by var explained gives the correlated functions too much weight. Equal weighting is chosen. (soil axis should not get too much weight)
EFdistance calculation DEPRECATED
From PCA of all functions, take most important axes and calculate euclidean distance of them (using equal weights).
We expect that one of the important axis - possibly the first to represent soil. As we have many soil functions compared to aboveground functions.
Selection of relevant PC axes : - check how many correlated functions there are - 9 correlations > 0.5 (see corrplots above) - (number of selected axes) = (number of PCA) - (number of correlated functions) - check if this number corresponds to around 90% variance explained.
# STANDARDISATION nimputed_grlfuns <- data.table::copy(imputed_grlfuns) # save an unscaled version for EFturnover / nestedness include <- names(imputed_grlfuns)[!names(imputed_grlfuns) %in% "Plotn"] imputed_grlfuns[, (include) := lapply(.SD, scale01),.SDcols=include] # PCA pc_grlfuns <- imputed_grlfuns[, ..include] pc_grlfuns <- as.data.frame(pc_grlfuns) rownames(pc_grlfuns) <- imputed_grlfuns$Plotn pca_grlfuns <- stats::prcomp(pc_grlfuns, scale=F) summ <- summary(pca_grlfuns) summ$importance[3, ] <= 0.9 # take axes 1 - 11 capture.output(summary(pca_grlfuns), file = "calc_multifun_from_imputed_pca_summary.txt") par(mfrow = c(2, 2)) plot(pca_grlfuns) biplot(pca_grlfuns) # plot cumulative prop of variance explained (same as output by pca summary) plot(cumsum(pca_grlfuns$sdev^2 / sum(pca_grlfuns$sdev^2)), pch = 16, xlab = "PCA axes", ylab = "Cumulative Proportion of Variance", xaxt = "n", yaxt = "n") axis(1, at = seq(1, 17, 2)) # axis(2, at = c(seq(0, 0.8, 0.2), 0.9161965, 1), labels = T) axis(2, at = seq(0, 1, 0.2)) axis(2, at = 0.9161965, labels = "0.92", las = 2, cex = 0.04) abline(v = 11, lty = 2) abline(h = 0.9161965, lty = 2) # save as : "22-05-16_calc_multifun_from_imputed_pca_output_cumvar.pdf" dev.off()
Plotting first 5 axis
temp <- cbind(pc_grlfuns, "region" = sub("[0-9][0-9]", "", rownames(pc_grlfuns))) pcaplots <- list(autoplot(pca_grlfuns, data = temp, colour = "region"), autoplot(pca_grlfuns, data = temp, colour = "region", x = 3, y = 4), autoplot(pca_grlfuns, data = temp, colour = "region", x = 5, y = 6), autoplot(pca_grlfuns, data = temp, colour = "region", x = 7, y = 8)) plot_grid(plotlist = pcaplots, nrow = 2) rm(temp) # saved as : "21-03-09_calc_multifun_from_imputed_pca_autoplot_axes_1_to_7.pdf" p <- plot(pca_grlfuns) labelsize <- 3 p1 <- ggplot2::autoplot(pca_grlfuns, loadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize) p2 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize, x = 3, y = 2) p3 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize, x = 4, y = 3) p4 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, loadings.label.size = labelsize, x = 4, y = 5) cowplot::plot_grid(plotlist = list(p1, p2, p3, p4), labels = c("A", "B", "C", "D"), align = T) # cowplot::plot_grid(plotlist = list(p1, p2), labels = c("A", "B"), align = T) # saved as : "calc_multifun_from_imputed_pca_axes_1_to_5.pdf" p5 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, x = 5, y = 6) p6 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, x = 6, y = 7) p7 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, x = 7, y = 8) p8 <- ggplot2::autoplot(pca_grlfuns, loadingsloadings = T, loadings.label= T, colour = "gray", loadings.colour = "black", loadings.label.colour = "black", loadings.label.repel = T, x = 8, y = 9) cowplot::plot_grid(plotlist = list(p5, p6, p7, p8), labels = c("E", "F", "G", "H"), align = T) # saved as : "23-11-09_calc_multifun_from_imputed_pca_axes_5_to_9.pdf"
rm(p); rm(p1); rm(p2); rm(p3); rm(p4); rm(p5)
axis selection : Use axes until 90% variance explained : take axes 1 - 11
pca_grlfuns <- pca_grlfuns$x[, 1:11] # euclidean distance EFmaster <- vegan::vegdist(pca_grlfuns, method = "euclid") EFmaster <- as.matrix(EFmaster) EFmaster[!lower.tri(EFmaster)] <- NA EFmaster <- reshape2::melt(EFmaster, value.name = "EFdistance") EFmaster <- EFmaster[!is.na(EFmaster[, 3]),] EFmaster <- data.table::data.table(EFmaster) # standardise EFmaster[, EFdistance := scale01(EFdistance)] saveRDS(EFmaster, "EFdistance.rds") rm(pca_grlfuns); rm(pc_grlfuns)
EFturnover and EFnestedness
Based on a 50% treshold, classify for each function and each plot presence or absence of the given function. Based on the scaled functions.
# percentage of measured functions that exceed a given treshold of their maximum observed level across all study sites. # maximum observed level: average of the top five sites. threshold <- 0.1 # tresholds used in Solivieres paper : 0.25,0.5,0.75,0.9
Plot for visualisation
# plot all 18 functions on one page # x <- nimputed_grlfuns$Root.biomass require(cowplot) for(i in colnames(nimputed_grlfuns)[which(colnames(nimputed_grlfuns) != "Plotn")]){ x <- nimputed_grlfuns[, get(i)] names(x) <- nimputed_grlfuns$Plotn # calc mean of top 5 values upper <- mean(sort(x, decreasing = T, na.last=T)[1:5]) upperq <- quantile(x, probs = threshold) df <- cbind(len = x, deleteme = x) p <- ggplot2::ggplot(df, ggplot2::aes(y = len, x= rownames(df))) + ggplot2::geom_bar(stat="identity") + ggplot2::geom_abline(intercept = upper, slope = 0, color = "lightgray") + ggplot2::geom_abline(intercept = upper*threshold, slope = 0, color =cols[1]) + # turkis ggplot2::geom_abline(intercept = upperq, slope = 0, color = cols[2]) + # orange theme(axis.text.x=element_blank()) + xlab("Plots") + ylab(i) plotname <- paste("plot", i, sep = "") assign(plotname, p) rm(p) } rm(i) plotnames <- paste("plot", colnames(nimputed_grlfuns)[which(colnames(nimputed_grlfuns) != "Plotn")], sep = "") plot_grid(plotlist = list(plotBiomass, plotdung.removal, plotGroundwater.recharge, plotherbivory.20172018, plotLitter.decomposition, plotN_leaching_risk2015, plotpathogen.infection, plotPhosphatase, plotPotential.nitrification), align = T) plot_grid(plotlist = list( plotcaterpillars.predation, plotRoot.biomass, plotRoot.decomposition, plotseed.depletion, plotsoilCflxs, plotP_leaching_risk_comb, plotsoilNitrateflxs, plotsoilAmmoniaflxs), align = T) # saved name : 21-03-09_calc_multifun_from_imputed_EFbeta_threshold_method_visual_tr01_part1 , change to tr = 0.7 etc, change part 1 and part2 # barplot(x, main="Root Biomass", xlab="Plots") # abline(a=upper,b=0,col="lightgray") # abline(a = upperq, b = 0, col = "darkgreen") # abline(a=upper*threshold, b=0, col="orange") rm(df); rm(plotBiomass, plotdung.removal, plotGroundwater.recharge, plotherbivory.20172018, plotLitter.decomposition, plotN_leaching_risk2015, plotpathogen.infection, plotPhosphatase, plotPotential.nitrification, plotcaterpillars.predation, plotRoot.biomass, plotRoot.decomposition, plotseed.depletion, plotsoilCflxs, plotP_leaching_risk_comb, plotsoilNitrateflxs, plotsoilAmmoniaflxs); gc()
With threshold max 5 method
include <- names(nimputed_grlfuns)[!names(nimputed_grlfuns) %in% "Plotn"] pa_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for max 5 functions threshold <- 0.9
Note : to make EF master for all thresholds : change manually above and run again and again the chunk below AND above.
# CONVERT TO PRESENCE ABSENCE pa_grlfuns[, (include) := lapply(.SD, function(c) calc_presenceabsence(c, threshold = threshold, type = "max5")), .SDcols = include] pa_grlfuns[, (include) := (.SD * 1), .SDcols = include] setcolorder(pa_grlfuns, neworder = c("Plotn", include)) # saveRDS(pa_grlfuns, file = paste(pathtoout, paste(paste("/single_functions_passing_threshold_", threshold, sep = ""), "Rds", sep = "."), sep = "")) # CALC BETADIVERSITY pa_grlfuns <- data.frame(pa_grlfuns[, !"Plotn", with=F], row.names = pa_grlfuns$Plotn) pa_grlfuns <- BetaDivMultifun::beta.pair_zerospecies(pa_grlfuns, index.family = "sorensen") EFbeta <- pa_grlfuns$beta.sor t <- as.matrix(EFbeta) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = paste("EFbeta", threshold, sep = "_")) EFbeta <- t[!is.na(t[, 3]), ] if(!exists("EFmaster")){ # if EFmaster does not exist yet, create EFmaster <- EFbeta } else { EFmaster <- merge(EFmaster, EFbeta, by = c("Var1", "Var2"), all = T) } EFturnover <- pa_grlfuns$beta.sim t <- as.matrix(EFturnover) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = paste("EFturnover", threshold, sep = "_")) EFturnover <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, EFturnover, by = c("Var1", "Var2"), all = T) EFnestedness <- pa_grlfuns$beta.sne t <- as.matrix(EFnestedness) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = paste("EFnestedness", threshold, sep = "_")) EFnestedness <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, EFnestedness, by = c("Var1", "Var2"), all = T)
standardise and save
# saveRDS(EFturnover, "EFturnover.rds") # saveRDS(EFnestedness, "EFnestedness.rds") saveRDS(EFmaster, "EFmaster_all_thresholds.rds") # if all threshods were calculated and added to EFmaster saveRDS(EFmaster, "EFmaster.rds")
Plot components of EFbeta
# # check EFmaster # EFcheck <- data.table::copy(data.table::data.table(EFmaster)) # EFcheck[, check := EFturnover + EFnestedness - EFbeta] # plot(EFcheck$check, ylim = c(-1, 1)) # rm(EFcheck) # # TRUE : EFturnover + EFnestedness - EFbeta is always 0 (or very near 0) plotEFbeta <- reshape2::melt(EFmaster[, c("Var1", "Var2", "EFbeta_0.5", "EFturnover_0.5", "EFnestedness_0.5")], id = c("Var1", "Var2")) p <- ggplot2::ggplot(data = plotEFbeta, aes(x = variable, y = value)) + geom_violin() + geom_jitter(shape=".", position=position_jitter(0.2), color = "gray") + labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "threshold = 0.5", xlab = "") + scale_x_discrete(labels=c("beta", "turnover", "nestedness")) + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) # Plot all thresholds plotEFbeta <- data.table(reshape2::melt(EFmaster, id = c("Var1", "Var2"))) plotEFbeta[, type := sub("_[0-9].[0-9]", "", variable)] p1 <- ggplot2::ggplot(data = plotEFbeta[variable %in% unique(plotEFbeta$variable)[1:14]], aes(x = variable, y = value, fill = type)) + geom_violin() + geom_jitter(shape=".", position=position_jitter(0.2), color = "gray", alpha = 0.2) + theme(axis.text.x = element_text(angle = 90)) + labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "all thresholds", xlab = "") + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) + # scale_fill_brewer(palette="Dark2") scale_fill_manual(values = c("#D95F02", "#1B9E77", "#7570B3")) p2 <- ggplot2::ggplot(data = plotEFbeta[variable %in% unique(plotEFbeta$variable)[14:28]], aes(x = variable, y = value, fill = type)) + geom_violin() + geom_jitter(shape=".", position=position_jitter(0.2), color = "gray", alpha = 0.2) + theme(axis.text.x = element_text(angle = 90)) + labs(title="Threshold-based dissimiliarity of Ecosystem Functions", subtitle = "all thresholds", xlab = "") + theme(axis.title.x = element_blank(), axis.title.y = element_blank()) + # scale_fill_brewer(palette="Dark2") scale_fill_manual(values = c("#D95F02", "#1B9E77", "#7570B3")) plot_grid(p1, p2, nrow = 2) # save as : 21-10-04_calc_multifun_from_imputed_all_EFdistances_violinplot
Plot various thresholds
# only run if EFmaster contains various thresholds temp <- c(grep("turnover|nestedness", names(EFmaster), value = T)) # temp <- grep("turnover", names(EFmaster), value = T) # pairs(EFmaster[, ..temp], upper.panel = NULL ) EFmaster <- data.table(EFmaster) m <- cor(EFmaster[, ..temp], method = "spearman") pdf("calc_multifun_from_imputed_all_EFdistances_corrplot.pdf", paper = "a4r") corrplot::corrplot(m, type = "lower", addCoef.col = "black", method = "color", diag = F, tl.srt = 50, tl.col = "black", order = "hclus", tl.cex = 0.5, cl.cex = 0.4, number.cex = 0.4) dev.off()
# EFturnover and nestedness are additive plot(EFmaster$EFbeta_0.5, EFmaster$EFturnover_0.5 + EFmaster$EFnestedness_0.5)
With quantile method DEPRECATED
threshold <- 0.5 include <- names(nimputed_grlfuns)[!names(imputed_grlfuns) %in% "Plotn"] paq_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for quantile functions # CONVERT TO PRESENCE ABSENCE # quantile method paq_grlfuns[, (include) := lapply(.SD, function(c) calc_presenceabsence(c, threshold = threshold, type = "quantile")), .SDcols = include] paq_grlfuns[, (include) := (.SD * 1), .SDcols = include] # CALC BETADIVERSITY paq_grlfuns <- data.frame(paq_grlfuns[, !"Plotn", with = F], row.names = paq_grlfuns$Plotn) paq_grlfuns <- betapart::beta.pair(paq_grlfuns, index.family = "sorensen") EFturnoverQ <- paq_grlfuns$beta.sim t <- as.matrix(EFturnoverQ) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFturnover_median") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) EFnestednessQ <- paq_grlfuns$beta.sne t <- as.matrix(EFnestednessQ) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFnestedness_median") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) EFbetaQ <- paq_grlfuns$beta.sor t <- as.matrix(EFbetaQ) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFbeta_median") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) saveRDS(EFmaster, "EFquantile_median.rds")
abundance betadiversity DEPRECATED
Abundance betadiversity of z-scores of functions.
abund_grlfuns <- data.table::copy(nimputed_grlfuns) # dataset for quantile functions # CALC BETADIVERSITY abund_grlfuns <- data.frame(abund_grlfuns[, !"Plotn", with = F], row.names = abund_grlfuns$Plotn) # scale to be between 0 and 1000 (promille) and round to 0 decimals # betapart does not allow negative values, and decimals do not make sense. abund_grlfuns <- apply(abund_grlfuns, 2, function(x) round(scale01(x)*1000, digits = 0)) abund_grlfuns <- betapart::beta.pair.abund(abund_grlfuns, index.family = "bray") temp <- abund_grlfuns$beta.bray t <- as.matrix(temp) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFabund_bray") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) temp <- abund_grlfuns$beta.bray.bal t <- as.matrix(temp) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFabund_bray_bal") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) temp <- abund_grlfuns$beta.bray.gra t <- as.matrix(temp) t[!lower.tri(t)] <- NA t <- reshape2::melt(t, value.name = "EFabund_bray_gra") t <- t[!is.na(t[, 3]), ] EFmaster <- merge(EFmaster, t, by = c("Var1", "Var2")) saveRDS(EFmaster[, .(Var1, Var2, EFabund_bray, EFabund_bray_bal, EFabund_bray_gra)], "EFbeta_abund.rds")
rm(nimputed_grlfuns); rm(imputed_grlfuns); rm(pa_grlfuns); rm(raw_grlfuns) rm(include) rm(p); rm(cols); rm(tres); rm(threshold); rm(upper); rm(upperq); rm(x)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.