#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
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)
The model name
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
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")
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")
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)
########### # 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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.