diff --git a/04_Additional_Figs.R b/04_Additional_Figs.R
new file mode 100644
index 0000000000000000000000000000000000000000..f20b40971b9aae7a7b98e5696e6bb4ecde70851c
--- /dev/null
+++ b/04_Additional_Figs.R
@@ -0,0 +1,173 @@
+library(tidyverse)
+
+#### Ancillary parsing functions ####
+# split final Dataset into a list, one element for each simulated dataset
+split.by.dataset <- function(toparse){
+ x <- read_lines(toparse)
+ tag <- "^SIMULATED DATA [0-9]+$"
+ last.tag <- "^FINAL RESULTS$"
+ splits <- grep(pattern=tag, x=x)
+ last.split <- grep(pattern=last.tag, x=x)-2
+ dataset.nr <- as.numeric(unlist(regmatches(x[splits], gregexpr("[[:digit:]]+", x[splits]))))
+
+ split.l <- list()
+ for(i in 1:length(splits)){
+ split.l[[i]] <- splits[i]: (c(splits, last.split)[i+1]-1)
+ }
+
+ x2 <- lapply(split.l, function(l){return(x[l])}) ### splits files into 100 datasets
+ names(x2) <- paste("Dataset", dataset.nr)
+ return(x2)
+}
+
+## split each element from the above list into introduction & 14 trait x environmental effect combinations
+split.trait.env.combos <- function( y ){
+ tag <- "^RESULTS WITH SIMULATED DATA [0-9]+$"
+ splits <- c(0, grep(pattern=tag, x=y))
+ dataset.nr <- as.numeric(unlist(regmatches(y[splits], gregexpr("[[:digit:]]+", y[splits]))))
+
+ split.l2 <- list()
+ for(i in 1:length(splits)){
+ split.l2[[i]] <- splits[i]: (c(splits, length(y))[i+1]-1)
+ }
+ trait.weights <- grep(pattern="Trait weights", x=y)
+ trait.weights <- regmatches(y[trait.weights], gregexpr("[[:digit:]]+", y[trait.weights]))
+ trait.weights <- sapply(trait.weights, function(t){paste(paste0("t",1:length(t),"=", t), collapse=" ")})
+
+ env.weights <- grep(pattern="Ecosystem variable weights", x=y)
+ env.weights <- regmatches(y[env.weights], gregexpr("[[:digit:]]+", y[env.weights]))
+ env.weights <- sapply(env.weights, function(e){paste(paste0("e",1:length(e),"=", e), collapse=" ")})
+
+ all.names <- c(paste("Simulated dataset", unique(dataset.nr)),
+ paste("Simulated dataset", dataset.nr, "-", trait.weights, "-",env.weights))
+ x3 <- lapply(split.l2, function(l){return(y[l])}) ### splits files into 14 results (trait x env combos)
+ names(x3) <- all.names
+ return(x3)
+ }
+
+## parse alpha beta and gamma diversity, as well as correlation coefficients from output of previous function
+extract.results <- function(z){
+ get.abg <- function(z, tag_string) {
+ W <- "Diversity descriptors of composition matrix W"
+ Y <- "Diversity descriptors of Beals smoothed composition matrix Y"
+ X <- "Diversity descriptors of fuzzy-weighted composition matrix X"
+ tag <- get(tag_string)
+ abg.lines <- z[grep(pattern=tag, x=z)+1:4]
+ abg.res <- NULL
+ if(length(abg.lines)>0){
+ abg.res <- regmatches(abg.lines, gregexpr("-*[[:digit:]]+\\.[[:digit:]]+", abg.lines))
+ abg.res <- sapply(abg.res, function(a){a[length(a)]})
+ names(abg.res) <- paste(c("alpha", "beta", "gamma", "propbeta"), tag_string, sep=".")}
+ return(abg.res)
+ }
+
+
+ get.cor <- function(z, tag_string) {
+ Y <- "Correlation coefficients between PCOA ordination community scores \\(based on Beals-smoothed composition matrix Y\\) and other community or ecosystem descriptors"
+ W <- "Correlation coefficients between PCOA ordination community scores \\(based on composition matrix W\\) and other community or ecosystem descriptors"
+ tag <- get(tag_string)
+ cor.lines <- z[grep(pattern=tag, x=z)+2:18]
+ cor.res <- NULL
+ if(length(cor.lines)>0){
+ cor.res <- unlist(regmatches(cor.lines, gregexpr("-*[[:digit:]]+\\.[[:digit:]]+", cor.lines)))
+ cor.with <- unlist(regmatches(cor.lines, gregexpr("CWM trait [[:digit:]]|Ecosystem variable [[:digit:]]", cor.lines)))
+ cor.with <- gsub(pattern="Ecosystem variable ", replacement="e", x=cor.with)
+ cor.with <- gsub(pattern="CWM trait ", replacement="t", x=cor.with)
+
+ names(cor.res) <- paste("cor", tag_string, cor.with, rep(paste0("a", 1:3), each=5), sep=".")
+ }
+ return(cor.res)
+ }
+
+
+ bind.res <- function(zz){
+ return(c(get.abg(zz, "W"),
+ get.abg(zz, "Y"),
+ get.abg(zz, "X"),
+ get.cor(zz, "Y"),
+ get.cor(zz, "W")))
+ }
+
+ tmp <- lapply(z, bind.res)
+ return( tmp %>%
+ bind_rows %>%
+ mutate(dataset=names(tmp)))
+ #mutate(dataset=rep(names(tmp), lapply(tmp, length))))
+}
+
+## Function to create Figure SXXV
+create.panel <- function(x){
+ gg.betaW <- ggplot(data=x %>%
+ dplyr::filter(metric=="propbeta") %>%
+ dplyr::filter(matrix=="W")) +
+ geom_vline(aes(xintercept=0), lty=2, alpha=0.6) +
+ geom_hline(aes(yintercept=0), lty=2, alpha=0.6) +
+ geom_density(aes(value)) +
+ xlab("Proportional Beta Diversity (W)") +
+ xlim(c(0,1.1)) +
+ ylim(c(0,60)) +
+ theme_classic()
+
+ gg.betaY <- gg.betaW %+%
+ (x %>%
+ dplyr::filter(metric=="propbeta") %>%
+ dplyr::filter(matrix=="Y")) +
+ xlab("Proportional Beta Diversity (Y)")
+
+ gg.corW <- gg.betaW %+%
+ (x %>%
+ dplyr::filter(metric=="cor") %>%
+ dplyr::filter(matrix=="W") %>%
+ dplyr::filter(with %in% c("e1", "e2", "e3")) %>%
+ group_by(dataset, metric, matrix) %>%
+ summarize(value=max(abs(value))))+
+ xlab("Cor(WE)") +
+ ylab(NULL)
+
+ gg.corY <- gg.betaW %+%
+ (x %>%
+ dplyr::filter(metric=="cor") %>%
+ dplyr::filter(matrix=="Y") %>%
+ dplyr::filter(with %in% c("e1", "e2", "e3")) %>%
+ group_by(dataset, metric, matrix) %>%
+ summarize(value=max(abs(value)))) +
+ xlab("Cor(YE)") +
+ ylab(NULL)
+
+ gg.panel <- cowplot::plot_grid(gg.betaW, gg.corW,
+ gg.betaY, gg.corY,
+ nrow=2, rel_widths = c(1,1.06))
+ return(gg.panel)
+}
+
+
+
+
+
+### Parse output file ####
+mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XY_SampleSize_Main=040_Inter=00_Corr=00_v21169"
+myfiles <- list.files(path=mypath, pattern = "FinalSimulatedData.txt", recursive = T, full.names = T)
+myfiles <- myfiles[grepl("_new", x=myfiles)]
+
+for(i in 1:length(myfiles)){
+ toparse <- myfiles[i]
+ sampleN <- regmatches(toparse, gregexpr("N=[[:digit:]]+", toparse))[[1]]
+ x2 <- split.by.dataset(toparse)
+ x3 <- lapply(x2, split.trait.env.combos)
+
+ div.summary <- lapply(x3, extract.results) %>%
+ bind_rows %>%
+ dplyr::relocate(dataset) %>%
+ pivot_longer(-dataset, names_to="metric", values_to="value") %>%
+ filter(!is.na(value)) %>%
+ separate(dataset, into=c("dataset", "trait", "env"), sep=" - ") %>%
+ dplyr::relocate(metric, .after=dataset) %>%
+ separate(metric, into=c("metric", "matrix", "with", "PCOA"), sep="\\.") %>%
+ mutate(value=as.numeric(value))
+
+ gg.out <- create.panel(div.summary)
+ ggsave(filename = paste0("_pics/R1/FigSXXZ_Panel_BetaCor_", sampleN, ".png"),
+ width=6, height=5, device="png", dpi = 300, plot = gg.out)
+}
+
+