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

knitr::opts_chunk$set(echo = TRUE)

Aim

Create plots which are bulit from multiple models. This script is only runned once.

Once-in lifetime plots (only run once)

Requirements

if(!exists("model_names")){
  sections_to_be_loaded <- c()
  source("vignettes/analysis_nonpublic.R")
}

Reads input from the GDM modelling as well as potentially from the modelwise plots.

TODO

• split into main and appendix figures
• [x] Move the creation of the average-threshold computation to a separate file and put that to the group of scripts that only need to be runned once. DONE : threshold script
• Check if the output format of the average over threshold models can be used as input here.

Legends

###########
# legends
###########
biotic_legend <- get_nice_legend(type = "biotic") # plot with plot_grid(biotic_legend)
abiotic_legend <- get_nice_legend(type = "abiotic")

all_legends <- plot_grid(biotic_legend, abiotic_legend, ovL, nrow = 4, rel_heights = c(0.3, 0.3, 0.1, 0.1))
all_legends <- plot_grid(NULL, NULL, NULL, NULL,  all_legends, NULL, NULL, NULL, nrow = 3, rel_heights = c(0.1, 0.8, 0.1))
ggsave(filename = paste(pathtoout, "/", "all_legends_incl_overview.pdf", sep = ""), plot = all_legends, device=cairo_pdf,
       width = 8.27, height = 11.96, units = "in", dpi = 900)

single Ecosystem Functions

heatmap single EF

1. Get dataset of effects sizes at maximum distance (variable importane) for all single EF models.
2. Plot a heatmap from this dataset.

###########
# load data
###########
singleEFdist <- readRDS(paste(pathtodata, "/data_assembly/output_data/singleEFdist.rds", sep = "")) # already scaled
i <- 1
ni <- colnames(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")][i]
i <- paste(pathtodata, "/analysis/output_datasets/gdm_", ni, "_LUI_input.Rds", sep = "")
gdmin <- readRDS(i)
permut <- "F"

gdm_output <- gdm::gdm(gdmin, geo = F, splines = NULL, knots = NULL)
exSplines <- gdm::isplineExtract(gdm_output)

# take maximum of each variable for barplot
maxsplines <- apply(exSplines$y, 2, max, na.rm = TRUE)
# get significances
if(permut == T){sign <- gdmperm[[3]][,1]}

# create restults table
restab <- create_restab()
data.table::setnames(restab, old = "maxsplines", new = ni)
restab2 <- data.table::copy(restab)

i <- 2
for(i in 2:(length(colnames(singleEFdist))-2)){
  print(i)
  ni <- colnames(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")][i]
  i <- paste(pathtodata, "/analysis/output_datasets/gdm_", ni, "_LUI_input.Rds", sep = "")
  gdmin <- readRDS(i)
  permut <- "F"
  gdm_output <- gdm::gdm(gdmin, geo = F, splines = NULL, knots = NULL)
  exSplines <- gdm::isplineExtract(gdm_output)
  # take maximum of each variable for barplot
  maxsplines <- apply(exSplines$y, 2, max, na.rm = TRUE)
  # get significances
  if(permut == T){sign <- gdmperm[[3]][,1]}

  restab <- create_restab()
  restab <- data.table::data.table(restab)
  data.table::setnames(restab, old = "maxsplines", new = ni)
  include <- c(ni, "names")
  restab2 <- merge(restab2, restab[, ..include], by = "names", all.x = T)
}
restab2 <- data.table::data.table(restab2)

# report the number of functions affected by each group / affected by above- belowground / abiotic
singleEFnames <- names(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")]
x <- apply(restab2[, ..singleEFnames], 1, function(x) sum(x > 0)) # count how many models have effect size > 0
restab2[, n_affected_fun := x]
rm(x)
#TODO save this dataset?? (July 2023)
#TODO HERE

Create heatmap

###########
# Create heatmap
###########
include <- c("nicenames", colnames(singleEFdist)[!colnames(singleEFdist) %in% c("Var1", "Var2")], "n_affected_fun")

## chose colour palette ----
# pal <- rev(RColorBrewer::brewer.pal(9, "RdYlBu"))
pal <- RColorBrewer::brewer.pal(9, "YlOrRd")
# pal <- c("white", RColorBrewer::brewer.pal(9, "YlOrRd"))[-2]
# colgrad <- melt(restab2[, ..include], id.vars = c("nicenames", "n_affected_fun"))
# colgrad <- colgrad$value  

## biotic aboveground sub-plot ----
df <- restab2[ground == "a", ]
df <- data.table::melt(df[, ..include], id.vars = c("nicenames", "n_affected_fun"))
df[value == 0, value := NA] # 0 values are set to NA to make them appear grey in the plot

a <- ggplot(data = df, aes(x = variable, y = nicenames, fill = value)) +
  geom_tile() +
  scale_fill_gradientn(colours = pal, na.value = "grey80") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, 
    size = 9, hjust = 1),
    plot.margin = margin(l = 50, r = 20),
    axis.line.x=element_blank(),
    axis.text.y = element_text(size=9, angle = 0),
    axis.ticks.x = element_blank(),
    panel.grid.major.x = element_line(color = "grey"),
    axis.title = element_blank(),
    legend.position = "left", # note : put legend left, remove below and put again after.
    legend.title = element_blank(),
    legend.text = element_text(size = 9)) +
  scale_x_discrete(position = "top") +
  coord_cartesian(clip = 'off') +  
  geom_text(data = unique(df[, .(nicenames, n_affected_fun)]), aes(label = n_affected_fun, x = length(unique(df$variable)) + 1, y = nicenames), color="black", size=3, inherit.aes = FALSE)
# get legend
legend_a <- get_legend(a)
# plot(legend_a)
# remove legend from plot
a <- a + theme(legend.position = "none")

## biotic belowground sub-plot ----
df <- restab2[ground == "b", ]
df <- data.table::melt(df[, ..include], id.vars = c("nicenames", "n_affected_fun"))
df[value == 0, value := NA]
b <- ggplot(data = df, aes(x = variable, y = nicenames, fill = value)) +
  geom_tile() +
  scale_fill_gradientn(colours = pal, na.value = "grey80") +
  theme(legend.position = "none", 
        axis.title = element_blank(),
        axis.text.y = element_text(size=9, angle = 0),
        plot.margin = margin(l = 50, r = 20),
        axis.text.x = element_blank(),
        axis.line.x=element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid.major.x = element_line(color = "grey")) +
    coord_cartesian(clip = 'off') +  
  geom_text(data = unique(df[, .(nicenames, n_affected_fun)]), aes(label = n_affected_fun, x = length(unique(df$variable)) + 1, y = nicenames), color="black", size=3, inherit.aes = FALSE)

## abiotic sub-plot ----
df <- restab2[ground == "x", ]
df <- data.table::melt(df[, ..include], id.vars = c("nicenames", "n_affected_fun"))
df[value == 0, value := NA]
c <- ggplot(data = df, aes(x = variable, y = nicenames, fill = value)) +
  geom_tile() +
  scale_fill_gradientn(colours = pal, na.value = "grey80") +
  theme(legend.position = "none", 
        axis.title = element_blank(),
        axis.text.y = element_text(size=9, angle = 0),
        plot.margin = margin(l = 50, r = 20),
        axis.text.x = element_blank(),
        axis.line.x=element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid.major.x = element_line(color = "grey")) +
      coord_cartesian(clip = 'off') +  
  geom_text(data = unique(df[, .(nicenames, n_affected_fun)]), aes(label = n_affected_fun, x = length(unique(df$variable)) + 1, y = nicenames), color="black", size=3, inherit.aes = FALSE)

## combining plot ----
# plot_grid(a, b, c, labels = c('A', '', 'B'), label_size = 12, nrow = 3, rel_heights = c(0.31, 0.51, 0.18), align = "v")
heatmap <- plot_grid(a, b, c, labels = c('A', '', 'B'), label_size = 12, nrow = 3, rel_heights = c(0.46, 0.42, 0.09), align = "v")

heatmap <- plot_grid(heatmap, legend_a, nrow = 1, rel_widths = c(0.9, 0.1))
# plot_grid(a, b, c, labels = c('A', '', 'B'), label_size = 12, nrow = 3, rel_heights = c(0.58, 0.44, 0.10), align = "v") # grayscale with legend

# plot_grid(p, b, q, labels = c('aboveground', 'belowground', 'abiotic'), label_size = 12, nrow = 3, rel_heights = c(0.31, 0.51, 0.18), align = "v")

ggsave(filename = paste(pathtoout, "/plot_GDM_singleEFdist_heatmap.pdf", sep = ""), plot = heatmap, device=cairo_pdf,
       width = 11.96, height = 8.27, units = "in", dpi = 900)

overview bars of single EFs

Summarise the effects on single EFs per group (biotic above/ belowground, abiotic).

Visualisation as stacked barplots.

# requires : restab2, created in chunk 1 under "# single EFs heatmap"
singfun_ovbar <- create_single_funs_overviewbars(restab2 = restab2, rel_colnames = singleEFnames)

# bring to ggplot format
ov_ab_singleEFmods <- melt(singfun_ovbar$above_below, id.vars = c("ground", "color"))
ov_tn_singleEFmods <- melt(singfun_ovbar$turnover_nestedess, id.vars = c("component", "color"))
# create plots
p_ov_ab <- create_single_funs_overviewbar_plot(singleF_restab = ov_ab_singleEFmods, pos = "stack") # note dodged barplots would be possible by setting pos = "dodge"
p_ov_tn <- create_single_funs_overviewbar_plot(singleF_restab = ov_tn_singleEFmods, pos = "stack")

p_ov_abtn <- plot_grid(p_ov_ab, p_ov_tn, labels = c("A", "B"))

# saved as : <date>_plot_GDM_singleEFdist_overview_bars
ggsave(filename = paste(pathtoout, "/plot_GDM_singleEFdist_overview_bars.pdf", sep = ""), 
       plot = p_ov_abtn, 
       device=cairo_pdf,
       width = 11.96, height = 8.27, units = "in", dpi = 900)

Summary numbers about number of average biotic/ abiotic/ LUI drivers driving single functions

# use vector "singleEFnames" (created above) to get names of single functions to address each column separately

# add column to restab2 with needed categories
restab2[, bio_abio_lui := lui_ground]
restab2[bio_abio_lui == "a", bio_abio_lui := "b"]

calc_avg_number_of_biotic_abiotic_drivers_of_single_functions(restab2 = restab2,
                                                              singleEFnames = singleEFnames)
# total available categories
table(restab2$bio_abio_lui)

Result : on average, 13.235294 biotic drivers, 1.666667 lui and 1.214286 abiotic drivers per single function.

Of total 30 biotic drivers, 2 lui and 2 abiotic drivers in the model.

Heatmap with overviewbars

Creation of a figure of single EF heatmaps, combined with the overviewbars.

# remove geom_text layer in order to get nice alignment
a_test <- a
a_test$layers[[2]] <- NULL

b_test <- b
b_test$layers[[2]] <- NULL

c_test <- c
c_test$layers[[2]] <- NULL

p_ov_ab_test <- p_ov_ab + theme(axis.title.x = element_blank(),
                              axis.line.x = element_blank(),
                              axis.text.x = element_blank(),
                              axis.ticks.x = element_blank(),
                              axis.text.y = element_text(size=9)) +
                        theme(plot.margin = unit(rep(0, 4), "pt"))
# p_ov_tn_test <- p_ov_tn + theme(axis.text.x = element_text(size=9))
p_ov_tn_test <- p_ov_tn + theme(axis.title.x = element_blank(),
                              axis.line.x = element_blank(),
                              axis.text.x = element_blank(),
                              axis.ticks.x = element_blank(),
                              axis.text.y = element_text(size=9)) +
                        theme(plot.margin = unit(rep(0, 4), "pt"))

singlefun_fullplot <- plot_grid(a_test,
                      b_test,
                      c_test,
                      p_ov_ab_test,
                      p_ov_tn_test,
                      ncol = 1,
                      align = "v",
                      rel_heights = c(0.2, 0.18, 0.05, 0.05, 0.05))
singlefun_fullplot

ggsave(filename = paste(pathtoout, "/plot_GDM_singleEFdist_heatmap_overview_bars.pdf", sep = ""), 
       plot = singlefun_fullplot,
       device=cairo_pdf,
       width = 11.96, height = 8.27, units = "in", dpi = 900)

clean up

# clean workspace
rm(a_test, b_test, c_test, p_ov_ab_test, p_ov_tn_test)
rm(a, b, c, df, exSplines, gdm_output, gdmin, heatmap, legend_a, ov_ab_singleEFmods, ov_tn_singleEFmods,
   p_ov_ab, p_ov_abtn, p_ov_tn, restab, restab2, singfun_ovbar, singlefun_fullplot, i, include, 
   ni, pal, singleEFdist, singleEFnames)

Threshold selection sd and devexpl

Aim : Show that a selection of a given threshold is impossible.

Requirements

sections_to_be_loaded <- c("thresholds", "results", "gdminput", "gdmoutput")
source("vignettes/analysis_nonpublic.R")

Variability in Response vs. deviance explained

Aim : over all thresholds, show the varaibility of the response variable as well as the percent deviance explained of the fitted model.

Included are EFnes and EFturnover model, thresholds 0.1 - 0.9.

We prefer a selection based on variability of the response variable, because this is not a post-hoc decision.

# # # # # # # # #
# CALC VARIABILITY OF RESPONSE VARIABLE
#
# require : input data : EFmaster_all_thresholds.Rds : load from script "analysis_nonpublic.R"
# is the response variable of all the models including a threshold
names <- grep(c("nestedness|turnover"), 
                model_names[model_class %in% c("multifun") & lui == "LUI", funs], 
                value = T)
# EFmaster[, ..names]

EFmaster_sd <- apply(EFmaster[, ..names], 2, function(x) sd(x, na.rm = F))
EFmaster_var <- apply(EFmaster[, ..names], 2, function(x) var(x, na.rm = F))

RespVar <- data.table(funs = names(EFmaster_var), sd = EFmaster_sd, var = EFmaster_var)

# add column with type of model
RespVar[, type := "EFdistance"]
RespVar[grep("nestedness", RespVar$funs), type := "nestedness"]
RespVar[grep("turnover", RespVar$funs), type := "turnover"]
# EFmaster[grep("beta", EFmaster$funs), type := "beta"] # note that this is not currently included

# # # # # # # # #
# MERGE VARIABILITY AND DEVIANCE EXPLAINED
#
# Read in output from the script "summarise_GDM_results.Rmd" (is done automatically)
devexpl <- merge(devexpl, model_names[, modelname, funs], by = "modelname")
#
model_results <- merge(RespVar, devexpl, by = "funs")
model_results[, modelname := NULL]
setorder(model_results, cols = "funs")
model_results[, type := factor(type, ordered = T, levels = c("EFdistance", "singleEF", "abund", "beta", "turnover", "nestedness"))] # assign order for plotting
model_results <- model_results[!is.na(GDM_dev_expl)]
model_results$funs <- sub("_", " ", model_results$funs)

# # # # # # # # #
# PLOT
#
# # # # # #
textsize_tr <- 13
# PLOT RESPONSE VARIANCE
#
sd_plot <- ggplot(model_results, aes(x = funs, y = sd, fill = type)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = textsize_tr),
        axis.text.y = element_text(size = textsize_tr),
        axis.title.x = element_blank(),
        text = element_text(size = textsize_tr),
        legend.title = element_blank()) +
  ylab("standard deviation")
# ggsave(filename = paste(pathtoout, "/plot_GDM_EFbeta_tresholds_sd.pdf", sep = ""), 
#        plot = sd_plot, device=cairo_pdf,
#        width = 11.96, height = 8.27, units = "in", dpi = 900)

# # # # # #
# PLOT DEVIANCE EXPLAINED
#
devexpl_plot <- ggplot(model_results, aes(x = funs, y = GDM_dev_expl, fill = type)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = textsize_tr),
        axis.text.y = element_text(size = textsize_tr),
        text = element_text(size = textsize_tr),
        axis.title.x = element_blank(),
        legend.title = element_blank()) +
  ylab("% deviance explained")
# ggsave(filename = paste(pathtoout, "/plot_GDM_EFbeta_tresholds_var_expl.pdf", sep = ""), 
#        plot = devexpl_plot, device=cairo_pdf,
#        width = 11.96, height = 8.27, units = "in", dpi = 900)

# # # # # #
# PLOT BOTH
#
# STANDARD DEVIATION AND PERCENT DEVIANCE EXPLAINED

threshold_selection <- plot_grid(sd_plot + theme(axis.text.x = element_blank()), 
                                 devexpl_plot,
                                 nrow = 2, align = "v",
                                 rel_heights = c(0.4, 0.6))

ggsave(filename = paste(pathtoout, "/plot_GDM_EFbeta_tresholds_var_expl_and_sd.pdf", sep = ""), 
       plot = threshold_selection, device=cairo_pdf,
       width = 11.96, height = 8.27, units = "in", dpi = 900)
ggsave(filename = paste(pathtoout, "/plot_GDM_EFbeta_tresholds_var_expl_and_sd.png", sep = ""), 
       plot = threshold_selection, device = "png",
       width = 11.96, height = 8.27, units = "in", dpi = 900)

Based on this plot, no ideal model can be selected. Therefore, we use all thresholds.

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