requirements : masterdiversity
in the file "master_diversity_alpha_beta.rds", see file analysis_nonpublic.R
- trlevels
- master_alpha
masterdiversity <- data.table(masterdiversity) require(cowplot) require(RColorBrewer) library(viridis)
number of species in plots
# load master_alpha from `analysis_nonpublic.R` master_alpha <- data.table::data.table(master_alpha) master_alpha <- master_alpha[Var1 %in% plotNAset, ] masterdiversity <- masterdiversity[Var1 %in% plotNAset,] masterdiversity <- masterdiversity[Var2 %in% plotNAset,]
Create table :Number of species per plot How many species are there on average per plot?
a <- data.frame(summary(master_alpha)) a <- a[-(1:6),-1] a <-data.table(a) a[, "type" := trimws(sub(":.*", "", a$Freq))] a[, "Value" := as.numeric(trimws(sub(".*:", "", a$Freq)))] a[, Freq := NULL] a <- dcast(a, Var2 ~ type, value.var = "Value") write.csv(a, file = "table_alphadiversities_per_plot.csv", quote = F, row.names = F)
# transform bacteria to 1/1000 [k#] gmaster_alpha <- master_alpha[, bacteria.RNA.alpha.k := bacteria.RNA.alpha / 1000] gmaster_alpha[, bacteria.RNA.alpha := NULL] gmaster_alpha <- melt(gmaster_alpha, id.vars = "Var1") gtrlevel_names <- unique(gmaster_alpha$variable) # ggplot(gmaster_alpha, aes(x = variable, y = value)) + # geom_point() i <- gtrlevel_names[2] ggplot(gmaster_alpha[variable == i], aes(x = variable, y = value)) + geom_point() + xlab("") + ylab("alpha diversity")
jitter <- position_jitter(width = 0.1, height = 0.1) for(t in trlevels){ t_colnames <- grep(t, colnames(masterdiversity), value = T) subset <- masterdiversity[, ..t_colnames] namesref <- c( "alpha" = paste(t, ".alpha", sep = ""), "alpha2" = paste(t, ".alpha.2", sep = ""), "beta.sim" = paste(t, ".beta.sim", sep = ""), "beta.sne" = paste(t, ".beta.sne", sep = "") ) plot_turnover <- ggplot(subset, aes(y = get(namesref[["alpha"]]), x = get(namesref[["alpha2"]]))) + geom_point(aes(colour = get(namesref[["beta.sim"]])), size = 4, position = jitter) + scale_colour_viridis_c() + ggtitle(label = t, subtitle = "coloured by turnover") + xlab("alpha diversity plot 1") + ylab("alpha diversity plot 2") + labs(colour = "beta") plot_nestedness <- ggplot(subset, aes(y = get(namesref[["alpha"]]), x = get(namesref[["alpha2"]]))) + geom_point(aes(colour = get(namesref[["beta.sne"]])), size = 4, position = jitter) + scale_colour_viridis_c() + ggtitle(label = t, subtitle = "coloured by nestedness") + xlab("alpha diversity plot 1") + ylab("alpha diversity plot 2") + labs(colour = "beta") plotname <- paste("plot_alpha_", t, sep = "") p <- plot_grid(plot_turnover, plot_nestedness,rel_widths = 1, rel_heights = 1, labels = c("to", "nes"), align = "h") assign(plotname, p) ggsave(filename=paste(plotname, "_coloured_by_beta.pdf", sep = ""), plot = p, device=cairo_pdf, width = 11.96, height = 8.27, units = "in", dpi = 900) } rm(p); rm(plot_nestedness); rm(plot_turnover); rm(subset); rm(namesref); rm(t); rm(t_colnames); rm(plotname)
richness_diff_beta <- data.table::copy(data.table(masterdiversity)) is.data.table(richness_diff_beta) for(t in trlevels){ t_colnames <- grep(t, colnames(richness_diff_beta), value = T) namesref <- c( "alpha" = paste(t, ".alpha", sep = ""), "alpha2" = paste(t, ".alpha.2", sep = ""), "beta.sim" = paste(t, ".beta.sim", sep = ""), "beta.sne" = paste(t, ".beta.sne", sep = ""), "alpha.diff" = paste(t, ".alpha.diff", sep = ""), "alpha.shared" = paste(t, ".alpha.shared", sep = ""), "gamma" = paste(t, ".gamma", sep = "") ) richness_diff_beta[, namesref["alpha.diff"] := abs(get(namesref[["alpha"]]) - get(namesref[["alpha2"]]))] richness_diff_beta[, namesref["alpha.shared"] := abs(get(namesref[["gamma"]]) - get(namesref[["alpha.diff"]]))] # deleteme <- c(namesref[["alpha"]], namesref[["alpha2"]]) # richness_diff_beta[, (deleteme) := NULL] }
# turnover and nestedness plot(richness_diff_beta$autotroph.beta.sim, richness_diff_beta$autotroph.beta.sne, main = "autotroph", xlab = "beta turnover", ylab = "beta nestedness") pdf(file = "alpha_plots_association_beta_nes_richness.pdf", paper = "a4r") par(mfrow = c(4, 3)) # autotrophs plot(richness_diff_beta$autotroph.beta.sne, richness_diff_beta$autotroph.alpha.diff, main = "autotroph", xlab = "beta nestedness", ylab = "diff in richness") plot(richness_diff_beta$autotroph.beta.sim, richness_diff_beta$autotroph.alpha.diff, main = "autotroph", xlab = "beta turnover", ylab = "diff in richness") plot(richness_diff_beta$autotroph.beta.sim + richness_diff_beta$autotroph.beta.sne, richness_diff_beta$autotroph.alpha.diff, main = "autotroph", xlab = "beta = beta sim + beta sne", ylab = "diff in richness") # secnodary consumers plot(richness_diff_beta$secondary.consumer.beta.sne, richness_diff_beta$secondary.consumer.alpha.diff, main = "secondary.consumer", xlab = "beta.sne", ylab = "diff in richness") plot(richness_diff_beta$secondary.consumer.beta.sim, richness_diff_beta$secondary.consumer.alpha.diff, main = "secondary.consumer", xlab = "beta.sim", ylab = "diff in richness") plot(richness_diff_beta$secondary.consumer.beta.sne + richness_diff_beta$secondary.consumer.beta.sim, richness_diff_beta$secondary.consumer.alpha.diff, main = "secondary.consumer", xlab = "beta sim + beta sne", ylab = "diff in richness") # omnivore protist plot(richness_diff_beta$omnivore.protist.beta.sne, richness_diff_beta$omnivore.protist.alpha.diff, main = "omnivore protist", xlab = "beta.sne", ylab = "diff in richness") plot(richness_diff_beta$omnivore.protist.beta.sim, richness_diff_beta$omnivore.protist.alpha.diff, main = "omnivore protist", xlab = "beta.sim", ylab = "diff in richness") plot(richness_diff_beta$omnivore.protist.beta.sim + richness_diff_beta$omnivore.protist.beta.sne, richness_diff_beta$omnivore.protist.alpha.diff, main = "omnivore protist", xlab = "beta turnover + beta nestedness", ylab = "diff in richness") # herbivore arthropod plot(richness_diff_beta$herbivore.arthropod.beta.sne, richness_diff_beta$herbivore.arthropod.alpha.diff, main = "herbivore.arthropod", xlab = "beta.sne", ylab = "diff in richness") plot(richness_diff_beta$herbivore.arthropod.beta.sim, richness_diff_beta$herbivore.arthropod.alpha.diff, main = "herbivore.arthropod", xlab = "beta.sim", ylab = "diff in richness") plot(richness_diff_beta$herbivore.arthropod.beta.sim + richness_diff_beta$herbivore.arthropod.beta.sne, richness_diff_beta$herbivore.arthropod.alpha.diff, main = "herbivore.arthropod", xlab = "beta turnover + beta nestedness", ylab = "diff in richness") dev.off()
# Correlation plot gcols <- grep("beta.s", colnames(richness_diff_beta), value = T) gcols <- c(gcols, grep("alpha.di", colnames(richness_diff_beta), value = T)) gcols <- c(gcols, grep("alpha.sha", colnames(richness_diff_beta), value = T)) M <- cor(richness_diff_beta[, ..gcols], use="complete.obs", method = "spearman") pdf("testalpha_plots_association_beta_nes_richness_corrplot.pdf", paper = "a4r") corrplot::corrplot(M, type = "upper", tl.col="black", tl.srt=90, diag = F, tl.cex = 0.4, order = "hclust") dev.off() rm(gcols)
In many cases, beta.sim is highly correlated with alpha diff.
linear model
# beta sim, beta sne, alpha diff, alpha shared lmrich <- data.table::copy(richness_diff_beta) lmrich[, "comparison" := paste(Var1, Var2, sep = ".")] lmcols <- c("comparison", gcols) lmrich <- lmrich[, ..lmcols] lmrich <- melt(lmrich, measure.vars = gcols) lmrich[grep("beta.sim", lmrich$variable), "type" := "beta.sim"] lmrich[grep("beta.sne", lmrich$variable), "type" := "beta.sne"] lmrich[grep("alpha.diff", lmrich$variable), "type" := "alpha.diff"] lmrich[grep("alpha.shared", lmrich$variable), "type" := "alpha.shared"] # add trophic levels for(t in trlevels){ lmrich[grep(t, lmrich$variable), "trlevel" := t] } lmrich <- dcast(lmrich, trlevel + comparison ~ type) summary(lm(beta.sne ~ alpha.diff + alpha.shared + beta.sim, data = lmrich)) summary(lm(beta.sne ~ alpha.diff + alpha.shared + beta.sim + trlevel, data = lmrich)) # beta.sne is mostly driven by beta.sim # but also by alpha.diff and alpha.shared # rather by alpha.diff than by alpha.shared # and this association strongly depends on the trophic level # --> secondary consumers are not having an effect --> they are driven by something else!
The difference is species which are in common.
Further considerations with gamma
plot(richness_diff_beta$autotroph.beta.sne, richness_diff_beta$autotroph.alpha.diff / richness_diff_beta$autotroph.gamma, main = "nestedness weighted by gamma") # because alpha.diff includes both nestedness and turnover. Take gamma away from alpha diff --> more correlation? plot(richness_diff_beta$autotroph.beta.sne + richness_diff_beta$autotroph.beta.sim, richness_diff_beta$autotroph.alpha.diff / richness_diff_beta$autotroph.gamma, main = "alpha weighted by gamma, full betadiversity") plot(richness_diff_beta$autotroph.beta.sne/ richness_diff_beta$autotroph.gamma, richness_diff_beta$autotroph.alpha.diff) plot(richness_diff_beta$autotroph.beta.sne, richness_diff_beta$autotroph.gamma)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.