In allanecology/BetaDivMultifun: Testing effects of multitrophic beta-diversity and land-use intensity (LUI) on ecosystem multifunctionality

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)

all trophic levels

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")

alpha-alpha plots, coloured by betadiversity.

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)

calc pairwise richness dist

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]
}

Associations among beta nestedness and richness

# 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)

allanecology/BetaDivMultifun documentation built on Nov. 9, 2023, 8:47 p.m.