diff --git a/02_Figures_Simulations.R b/02_Figures_Simulations.R
new file mode 100644
index 0000000000000000000000000000000000000000..15e9c74a30c64bd466eeead09d7f10e0e06644f4
--- /dev/null
+++ b/02_Figures_Simulations.R
@@ -0,0 +1,258 @@
+## this script extracts all the summary.txt files from all subfolder
+## and summarizes the output for each run, trait x environment combination, and statistics
+## It then plots the summarized output
+library(tidyverse)
+
+#### import function
+FormatData <- function(myfiles){
+ output <- NULL
+ for(ff in myfiles){
+ iter <- gsub(pattern="/Summary.txt$", replacement="", ff)
+ iter <- strsplit(iter, split = "_")[[1]]
+ iter <- as.integer(unlist(regmatches(iter, gregexpr("[[:digit:]]+", iter))))
+ tmp <- read_delim(paste(mypath, ff, sep="/"), delim="\t", col_names = F) %>%
+ dplyr::select(-X1, -X3, -X5, -X9, -X11, -X13) %>%
+ rename(simulated=X2, trait=X4, envir=X6, stat.type=X7, stat.obs=X8, pvalue=X10, SES=X12, exp.med=X14 ) %>%
+ mutate(stat.type=gsub(pattern = "^r\\(", replacement="", stat.type)) %>%
+ mutate(stat.type=gsub(pattern = "\\)\\=$", replacement="", stat.type)) %>%
+ mutate(stat.type=gsub(pattern = "_", replacement="\\.", stat.type)) %>%
+ mutate(trait=gsub(pattern="[[:space:]]+$", replacement="", trait)) %>%
+ mutate(envir=gsub(pattern="[[:space:]]+$", replacement="", envir)) %>%
+ mutate(main=iter[[1]]) %>%
+ mutate(inter=iter[[2]]) %>%
+ mutate(corr=iter[[3]]) %>%
+ dplyr::select(main:corr, everything())
+ output <- bind_rows(output, tmp)
+ }
+
+ outp.summary <- output %>%
+ dplyr::filter(!stat.type %in% c("TY", "XY", "XY.T", "XY.TR")) %>%
+ group_by(main, inter, corr, trait, envir, stat.type) %>%
+ summarize(stat.obs.med=median(stat.obs),
+ power=mean(pvalue<=0.05),
+ SES.med=median(SES),
+ exp.med.med=median(exp.med),
+ nsim=n()) %>%
+ bind_rows(output %>%
+ dplyr::filter(stat.type %in% c("TY", "XY", "XY.T", "XY.TR")) %>%
+ group_by(main, inter, corr, trait, stat.type) %>%
+ summarize(stat.obs.med=median(stat.obs),
+ power=mean(pvalue<=0.05),
+ SES.med=median(SES),
+ exp.med.med=median(exp.med),
+ nsim=n())) %>%
+ dplyr::select(main:stat.type, nsim, stat.obs.med:exp.med.med) %>%
+ arrange(stat.type, main, inter, corr, trait, envir)
+ return(outp.summary)
+}
+
+
+get.ntraits <- function(x){
+ tmp <- str_split(x, pattern = " ")[[1]]
+ return(length(tmp))
+}
+
+
+#### FIGURE 2 #####
+mypath <- "_data/Experiment_02Mar2020_FactorInteraction&TraitCorr"
+myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T)
+outp.summary <- FormatData(myfiles)
+
+
+#e41a1c ##red
+#377eb8 ##blue
+#ffffb2 ## yellow
+#4daf4a ##green
+#984ea3 ##violet
+#ff7f00 ##orange
+
+
+mypalette <- palette(c("#e41a1c", #1 - red)
+ "#ff7f00", #12 - orange
+ "#984ea3", #13 - violet
+ "##ffed6f", #2 - yellow
+ "#4daf4a", #23 - green
+ "#377eb8")) #3 - blue
+
+ggplot(data=outp.summary %>%
+ mutate(main=main/100) %>%
+ mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>%
+ mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY") %>%
+ dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>%
+ mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"),
+ labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn")))) + #%>%
+ #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ #mutate(inter=as.factor(inter))) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
+ #scale_colour_brewer(palette = "Dark2") +
+ scale_color_manual("Trait\ncomb.",
+ values=c("#e41a1c", #1 - red)
+ "#e6ab02", #2 - yellow
+ "#377eb8", #3 - blue
+ "#d95f02", #12 - orange
+ "#984ea3", #13 - violet
+ "#4daf4a" #23 - green
+ )
+ ) +
+ facet_grid(corr~inter) +
+ theme_bw() +
+ scale_x_continuous(name="Effect of factor e1 -> trait t1") +
+ scale_y_continuous(name="Power") +
+ theme(panel.grid = element_blank())
+
+ggsave(filename="_pics/Fig2_CorrInte_02March.png", width=6, height=5, device="png", dpi = 300, last_plot())
+
+
+
+
+#### FIGURE 3 #####
+mypath <- "_data/Experiment_04Mar2020_TraitNumber"
+myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T)
+outp.summary <- FormatData(myfiles) %>%
+ rename(ntraits=inter)
+
+## plotting power for XY with corr
+add.t.label <- function(x) {
+ x <- gsub(pattern=" ", replacement=" t", x=x, perl=T)
+ x <- paste0("t", x)
+ return(x)
+}
+
+outp.summary2 <- outp.summary %>%
+ ungroup() %>%
+ mutate(main=main/10) %>%
+ #mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>%
+ mutate(ntraits.lab=factor(ntraits, levels=1:5, labels=paste0("n. traits = ", 1:5))) %>%
+ rowwise() %>%
+ mutate(sel.ntraits=factor(get.ntraits(trait))) %>%
+ mutate(sel.ntraits.lab=factor(sel.ntraits, levels=levels(sel.ntraits), labels=paste0("Comb. tier = ", levels(sel.ntraits)))) %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY") %>%
+ mutate(trait=factor(trait)) %>%
+ mutate(trait=factor(trait, levels=levels(trait),
+ labels=add.t.label(levels(trait))))
+
+### rename null trait to "tn"
+outp.summary2 <- outp.summary2 %>%
+ mutate(tn.name=paste0("t", ntraits)) %>%
+ mutate(trait=str_replace(trait, pattern=tn.name, replacement="tn"))
+
+#reorder factors
+outp.summary2$trait <- factor(outp.summary2$trait, levels= c('t1', 't2','t3', 't4', 'tn',
+ 't1 t2', 't1 t3','t1 t4','t1 tn','t2 t3','t2 t4','t2 tn','t3 t4','t3 tn', 't4 tn',
+ 't1 t2 t3', 't1 t2 t4','t1 t2 tn','t1 t3 t4', 't1 t3 tn','t1 t4 tn','t2 t3 t4','t2 t3 tn', 't2 t4 tn','t3 t4 tn',
+ 't1 t2 t3 t4', 't1 t2 t3 tn','t1 t2 t4 tn','t1 t3 t4 tn','t2 t3 t4 tn',
+ 't1 t2 t3 t4 tn'))
+
+
+
+hugepalette0 <- c(RColorBrewer::brewer.pal(4, "Dark2"),
+ gray(0.2),
+ RColorBrewer::brewer.pal(10, "Paired"),
+ RColorBrewer::brewer.pal(10, "Set3"),
+ RColorBrewer::brewer.pal(5, "Pastel1"), "brown")
+#change tone of yellow of t1-t2-t4
+hugepalette0[17] <- "#ccebc5"
+
+hugepalette <- data.frame(trait=c('t1', 't2','t3', 't4', 'tn',
+ 't1 t2', 't1 t3','t1 t4','t1 tn','t2 t3','t2 t4','t2 tn','t3 t4','t3 tn', 't4 tn',
+ 't1 t2 t3', 't1 t2 t4','t1 t2 tn','t1 t3 t4', 't1 t3 tn','t1 t4 tn','t2 t3 t4','t2 t3 tn', 't2 t4 tn','t3 t4 tn',
+ 't1 t2 t3 t4', 't1 t2 t3 tn','t1 t2 t4 tn','t1 t3 t4 tn','t2 t3 t4 tn',
+ 't1 t2 t3 t4 tn'),
+ trait.col=factor(hugepalette0, levels=hugepalette0))
+
+outp.summary2 <- outp.summary2 %>%
+ left_join(hugepalette, by="trait")
+
+fig3 <- ggplot(data=outp.summary2) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait.col)) +
+ scale_x_continuous(name="Effect of factor e1 -> trait t1", n.breaks = 4) +
+ scale_y_continuous(name="Power") +
+ facet_grid(sel.ntraits.lab~ntraits.lab) +
+ scale_color_identity(guide = "legend",
+ labels= hugepalette$trait) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+
+
+
+ncols <- c(2,4,4,3,1)
+
+
+leg.list <- list()
+col.used <- 1
+for(tier in 1:5){
+ outp.summary.tier <- outp.summary2 %>% filter(sel.ntraits==tier)
+ ncombinations <- length(levels(factor(outp.summary.tier$trait)))
+ leg.list[[tier+1]] <- cowplot::get_legend(fig3 %+% outp.summary.tier +
+ guides(col=guide_legend(ncol=ncols[tier], byrow=TRUE))+
+ scale_color_identity(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)),
+ guide = "legend",
+ labels= hugepalette$trait[col.used:(col.used+ncombinations-1)])
+ )
+col.used <- col.used + ncombinations
+}
+leg.list[[tier+2]] <- NULL
+
+fig3.panel <- cowplot::plot_grid(fig3 + theme(legend.position = "none"),
+ cowplot::plot_grid(plotlist = leg.list, nrow = 7,
+ rel_heights = c(0.05, .2,.2,.2,.2,.2, 0.15), align="hv"),
+ nrow=1, rel_widths = c(0.6,.4))
+
+
+ggsave(filename="_pics/Fig3_TraitNumber.png",
+ width=9, height=7, device="png", dpi = 300, fig3.panel)
+
+
+
+#### FIGURE S4 #####
+mypath <- "_data/Experiment_08Jul2020_TraitSuppression&TraitCorr"
+myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T)
+outp.summary <- FormatData(myfiles)
+
+
+ggplot(data=outp.summary %>%
+ mutate(main=main/100) %>%
+ mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>%
+ mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY") %>%
+ dplyr::filter(trait %in% c("1", "2", "3", "4")) %>%
+ mutate(trait=factor(trait, levels=c("1", "2", "3", "4"),
+ labels=c("t1", "t2", "t3","tn")))) + #%>%
+ #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ #mutate(inter=as.factor(inter))) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
+ #scale_colour_brewer(palette = "Dark2") +
+ #scale_color_manual("Trait\ncomb.",
+ scale_color_manual("Trait",
+ values=c("#e41a1c", #1 - red)
+ "#e6ab02", #2 - yellow
+ "#4daf4a", #23 - green
+ "#377eb8", #3 - blue
+ "#d95f02", #12 - orange
+ "#984ea3" #13 - violet
+ )
+ ) +
+ facet_grid(corr~inter) +
+ theme_bw() +
+ scale_x_continuous(name="Effect of factor e1 -> trait t1") +
+ scale_y_continuous(name="Power") +
+ theme(panel.grid = element_blank())
+
+ggsave(filename="_pics/FigS4_Extra_CorrInte_08Jul20.png", width=6, height=5, device="png", dpi = 300, last_plot())
+
+
+
+
+
+
+
+
+
+
+
+
+