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

#TODO :  fill in code to work with significances, e.g. take out the non-significant maxsplines
#TODO : order of trophic groups : once all models are read, use effect size sequence and not alphabetically

Aim

This scripts visualises the GDM output. It is runned for each model individually.

Requirements : - GDM output as .rds, e.g. gdm_EFnestedness_LUI_input.Rds - find at planteco/.../BetaDivMultifun/analysis/output_datasets - cluster output as .rds, e.g. gdm_EFnestedness_LUI_permutation.Rds - TODO : not calculated yet (only for old)

Load data

The model name

select input

Depending on the model intented to run, chose the appropriate variables below. The model is selected from the table model_names.

Load table model_names.

# Load requirements from <..>nonpublic.R file
source("vignettes/analysis_nonpublic.R")

Chose model

# select a model to run, from the table `model_names` which is part of the R package data (see how to load in `analysis_nonpublic.R`)
if(exists("model_names_selection")){
  # model_names_selection is given from outside the script
  print("using selection of the model from outside the script")
} else {
  print("model selection is not given, default is EFturnover_0.7")
  model_names_selection <- model_names[which(model_names$modelname == "gdm_EFturnover_0.7_LUI"), ]
  print(c("working on model :", paste(model_names_selection, collapse = ", ")))
}

Load requirements

funs <- model_names_selection$funs
lui <- model_names_selection$lui
sections_to_be_loaded <- c("gdminput", "gdmoutput")
source("vignettes/analysis_nonpublic.R")
# will create "model_name" automatically
#    e.g. gdm_EFnestedness_0.7_LUI

permut <- F # T or F, depending wether permutations are available or not
#TODO : implement this

preparation

if(permut == T){
  stop("handling permutation values is currently not implemented")
  gdmperm <- readRDS(paste(paste(paste(pathtodata, "/analysis/output_datasets/pvalue_calculation/", sep = ""), model_name, sep = ""), "_perm.Rds", sep = ""))
  sign <- gdmperm[[3]][,1]
  sign <- data.table::data.table("names" = names(sign), "sign" = sign)
}


# get maxsplines
exSplines <- gdm::isplineExtract(gdmoutput)
maxsplines <- apply(exSplines$y, 2, max, na.rm = TRUE)

# define plot sequence
plotsequence_bio <- c("autotroph", "bacteria.RNA", "herbivore.arthropodsoillarvae", "secondary.consumer.arthropodsoillarvae", "herbivore.arthropod", "plantpathogen.fungi", "bacterivore.protist", "eukaryvore.protist", "omnivore.protist", "plantparasite.protist", "secondary.consumer", "decomposer.soilfungi", "pathotroph.soilfungi", "symbiont.soilfungi", "tertiary.consumer")
plotsequence_abio <- c(model_names_selection$lui, 
                       paste("delta", lui, sep = ""), "soil", "isolation", "geo")

Creating legend plots

Plots that are re-used in all the plots (Legend and model results)

###########
# create additional plots
###########
# model specifications and legends
labsplot <- ggplot(NULL) + #annotate("text", x = 0.8, y = 0.5, size=8, label = model_names_selection$modelname) +
  theme_void() + annotate("text", x = 0.4, y = 0.5, size = 3, 
                          label = paste(names(model_specs), model_specs, sep = ": ", collapse = "\n"))
biotic_legend <- get_nice_legend(type = "biotic") # plot with plot_grid(biotic_legend)
abiotic_legend <- get_nice_legend(type = "abiotic")

Line plots

naming convention : model_names_selection$modelname + "_lineplot_all_modelspec.pdf"

common y axis maximum is the maximum of the current model

###########
# create plot table
###########
    # one column with xaxis values and one column with yaxis values, and one column
    # containing names
x <- data.table(exSplines$x)
x[, ids := as.numeric(seq(1, nrow(x)))] # ids are required for later matching with y
y <- data.table(exSplines$y)
y[, ids := as.numeric(seq(1, nrow(y)))]
x <- melt(x, measure.vars = colnames(x), value.name = "xaxis", id.vars = "ids", variable.name = "names") # get to long format
y <- melt(y, measure.vars = colnames(y), value.name = "value", id.vars = "ids", variable.name = "names")
    # add x and y values
lineplot_data <- merge(x, y, by = c("ids", "names"))
lineplot_data[, ids := NULL]
lineplot_data <- lineplot_data[names != "ids",]
lineplot_data[, names := as.character(names)] # somehow required for merging
lineplot_data <- merge(lineplot_data, nicenames, by = "names", all.x = T)


###########
# create lineplots
###########
# above, below and abiotic subplots
p <- create_gdm_lineplot(lineplot_data[ground == "a"], legend = F, ymax = max(lineplot_data$value))
q <- create_gdm_lineplot(lineplot_data[ground == "b"], legend = F, ymax = max(lineplot_data$value))

# abiotic subplots
# LUI
luidat <- lineplot_data[names %in% model_names_selection$lui,]
deltaluidat <- lineplot_data[names %in% paste("delta", model_names_selection$lui, sep = ""),]
l <- create_gdm_lineplot(luidat, legend = F, ymax = max(lineplot_data$value))
dl <- create_gdm_lineplot(deltaluidat, legend = F, ymax = max(lineplot_data$value))

lgeo <- create_gdm_lineplot(lineplot_data[names %in% c("Geographic")], legend = F, ymax = max(lineplot_data$value))
lsoil <- create_gdm_lineplot(lineplot_data[names %in% c("edis_soil")], legend = F, ymax = max(lineplot_data$value))
lisol <- create_gdm_lineplot(lineplot_data[names %in% c("plot_isolation")], legend = F, ymax = max(lineplot_data$value))


###########
# add plots together
###########
# create classic lineplot
suppressWarnings(
  lin <- plot_grid(p, biotic_legend, q, labsplot, l, abiotic_legend,
                   labels = c('A', '', 'B', '', 'C', ''), nrow = 3,
                   rel_widths = c(0.6, 0.5))
  ) # suppresswarnings because of special characters
lin_abiotic <- plot_grid(dl, lgeo, lsoil, lisol,
                         labels = c("D", "", "", ""),
                         align = T, nrow = 1)
all_lines <- plot_grid(lin, lin_abiotic, nrow = 2, rel_heights = c(0.8, 0.2))


# Alternative plottings
# lin <- plot_grid(p, q, l, labels = c("A", "B", "C"), nrow = 3,
#                    align = T)
# suppressWarnings(
# lin <- plot_grid(p, labsplot, q, NULL, l,
#                  labels = c('A', '', 'B', '', 'C', ''), nrow = 3,
#                  rel_widths = c(0.7, 0.3))
# )
# p_annot <- plot_grid(biotic_legend, plot_grid(abiotic_legend, labsplot, NULL, nrow = 3))
# all_lines <- plot_grid(p, p_annot, q, l, dl, lgeo, lsoil, lisol,
#           labels = c("A", "B", "C", "D"),
#           nrow = 4)



###########
# save
###########
# `model_names_selection$modelname` + "_lineplot_all_modelspec.pdf"
# ggplot2::ggsave(plot = lin, 
#                 filename = paste(pathtoout, "/", model_names_selection$modelname, "_lineplot.pdf", sep = ""),
#                 device=cairo_pdf, width = 8.27, height = 11.96, units = "in", dpi = 900)

ggplot2::ggsave(plot = all_lines, 
                filename = paste(pathtoout, "/", model_names_selection$modelname, "_lineplot.pdf", sep = ""),
                device=cairo_pdf, width = 8.27, height = 11.96, units = "in", dpi = 900)

# ggplot2::ggsave(plot = lin_abiotic, filename = paste(pathtoout, "/", model_names_selection$modelname, "_lineplot_no_modelspec_abiotic.pdf", sep = ""),
#                 device=cairo_pdf, width = 11.96, height = 8.27, units = "in", dpi = 900)
# ggplot2::ggsave(plot = labsplot, filename = paste(pathtoout, "/", model_names_selection$modelname, "_lineplot_only_modelspec_abiotic.pdf", sep = ""),
#                 device=cairo_pdf, width = 11.96, height = 8.27, units = "in", dpi = 900)

Barplots

###########
# create table
###########
# maxsplines <- apply(exSplines$y, 2, max, na.rm = TRUE) # take maximum of each variable for barplot # is created above
restab <- create_restab0()
saveRDS(restab, file = paste(pathtoout, "/", model_names_selection$modelname, "_restab_data.RDS", sep = "")) #TODO no errorbars! is it really worth it to save this one?

###########
# produce barplots
###########
p <- create_bio_aboveground_barplot(type = "stacked") # type is  in c("stacked", "grouped")
b <- create_bio_belowground_barplot(type = "stacked")
q <- create_abio_barplot(type = "stacked") # here it does not matter if "stacked" or "grouped" is given.


###########
# produce overview bars
###########
# # # # #
# Produce data for plotting
df <- create_overviewbar_restab(restab, fun = "mean")
df2 <- create_overviewbar_restab(restab, fun = "sum")
df <- rbindlist(list(df, df2), use.names = T, fill = T)
rm(df2)
saveRDS(df, file = paste(pathtoout, "/", model_names_selection$modelname, "_overviewbar_data.RDS", sep = ""))
#TODO check if it's a problem in later code that instead of "averaged_scaled", the entry in the 
# column of the overviewbar data table is called "mean"
# checked with bash: never mentioned. keep this comment in case later problem occurs.

# # # # #
# Produce plot
ov <- create_overview_barplot(df)
ovL <- create_overview_barplot(df, legend = T)
# plot_grid(ovL) # to show the legend individually


###########
# add plots together
###########
multipanel_with_overview_bars <- plot_grid(p, b, q, ov, 
                                           labels = c('A', '', 'B', 'C'), 
                                             label_size = 12, nrow = 4,
                                             rel_heights = c(0.21, 0.32, 0.2, 0.15), 
                                           align = "v") # 2 summary bars


annotated_multipanel_with_overview_bars <- ggdraw(multipanel_with_overview_bars) +
  # draw_plot(labsplot, .65, .65, .35, .35) +
  draw_plot(ovL, 0.04, 0.02, .1, .2) + # (left, below, x, x)
  draw_plot(labsplot, 0.6, 0.6, 0.6, 0.6)

# extraplot <- plot_grid(labsplot, ov_legends)

###########
# save
###########
# -- probably delete below
#     # WITHOUT ANNOTATION
# # ggsave(filename=paste(model_names_selection$modelname, "_multipanel_barplot.pdf", sep = ""), plot = multipanel_without_overviews, device=cairo_pdf,
#        # width = 11.96, height = 8.27, units = "in", dpi = 900)
# # ggsave(filename = paste(model_names_selection$modelname, "_multipanel_barplot_overviewbars.pdf", sep = ""), plot = multipanel_with_overview_bars, device=cairo_pdf,
#        # width = 11.96, height = 8.27, units = "in", dpi = 900)
# 
#     # ANNOTATED PLOTS
# # ggsave(filename=paste(model_names_selection$modelname, "_multipanel_barplot_annotated.pdf", sep = ""), plot = annotated_multipanel_without_overviews, device=cairo_pdf,
#        # width = 11.96, height = 8.27, units = "in", dpi = 900)
# --

ggsave(filename = paste(pathtoout, "/", model_names_selection$modelname, "_barplot.pdf", sep = ""), plot = annotated_multipanel_with_overview_bars, device=cairo_pdf,
       width = 11.96, height = 8.27, units = "in", dpi = 900)

# ggsave(filename = paste(pathtoout, "/", model_names_selection$modelname, "_multipanel_barplot_overviewbars_extraplot_annotations.pdf", sep = ""), plot = extraplot, device=cairo_pdf,
#        width = 11.96, height = 8.27, units = "in", dpi = 900)

Note : the barplot of weighted average is done in the script GDM_multifun_thresholds.Rmd, together with pieces of code from this script.

allanecology/multiFunBetadivLuiGDM documentation built on Nov. 12, 2023, 6:16 a.m.