diff --git a/02_Figures.R b/02_Figures.R
new file mode 100644
index 0000000000000000000000000000000000000000..cee5870a65f46da0ac107dde4b84a68d29f240b9
--- /dev/null
+++ b/02_Figures.R
@@ -0,0 +1,295 @@
+## 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", "t3", "t1 t2", "t1 t3", "t2 t3")))) + #%>%
+ #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 = 6) +
+ 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)
+
+
+
+
+
+
+#### test to duplicate rows in outp.summary and draw multiple lines next to each other
+###create unique palette
+gg_color_hue <- function(n) {
+ hues = seq(15, 375, length = n + 1)
+ hcl(h = hues, l = 65, c = 100)[1:n]
+}
+n = 5 ## 31 combination of traits
+mycols = gg_color_hue(n)
+
+outp.summary2.multicolor <- outp.summary2 %>%
+ ungroup() %>%
+ uncount(as.numeric(as.character(sel.ntraits)), .id="replicate") %>%
+ group_by(ntraits, trait) %>%
+ ungroup() %>%
+ rowwise() %>%
+ mutate(which.traits=str_split(trait, pattern=" ")) %>%
+ mutate(trait.color=which.traits[replicate]) %>%
+ # left_join(data.frame(trait.color=c("t1", "t2", "t3", "t4", "t5"),
+ # color=mycols[1:5]),
+ # by="trait.color") %>%
+ #mutate(power=power + replicate/100) %>%
+ mutate(trait.repl = str_glue(as.character(trait)," r", replicate)) %>%
+ mutate(trait.repl = factor(trait.repl))
+
+(fig3.multicolor <- ggplot(data=outp.summary2.multicolor) +
+ geom_line(aes(x=main, y=power, group=trait.repl, col=trait.color, lty=as.factor(replicate))) +
+ scale_x_continuous(name="Effect of factor e1 -> trait t1") +
+ scale_y_continuous(name="Power") +
+ #scale_colour_brewer(palette = "Dark2") +
+ facet_grid(sel.ntraits.lab~ntraits) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+)
+
+
+
+
+
+#### Alternative plotting to create grobs with individual legend
+mydata <- outp.summary %>%
+ ungroup() %>%
+ rowwise() %>%
+ mutate(sel.ntraits=(get.ntraits(trait))) %>%
+ ungroup() %>%
+ dplyr::filter(stat.type=="XY") %>%
+ #filter(ntraits==3) %>%
+ #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+ #dplyr::filter(trait %in% c("1", "2", "3")) %>%
+ #dplyr::filter(trait %in% c("1", "2", "3", "1 2", "1 2 3")) %>%
+ mutate(ntraits=as.factor(ntraits)) %>%
+ mutate(sel.ntraits=as.factor(sel.ntraits))
+
+gglist <- list()
+tick <- 1
+for(sel in levels(mydata$sel.ntraits)){
+ for(nn in levels(mydata$ntraits)){
+ gglist[[tick]] <- ggplot(data=mydata %>%
+ filter(ntraits==nn) %>%
+ filter(sel.ntraits==sel)) +
+ geom_line(aes(x=main, y=power, group=trait, col=trait)) +
+ #scale_colour_brewer(palette = "Dark2") +
+ theme_bw() +
+ guides(fill=guide_legend(ncol=2)) +
+ theme(panel.grid = element_blank(),
+ legend.position = c(0.9,0.5),
+ legend.title = element_blank(),
+ legend.text = element_text(size=4),
+ legend.background = element_blank())
+ tick <- tick + 1
+ }
+}
+gglist[[1]]
+
+cowplot::plot_grid(plotlist=gglist,
+ nrow=5, ncol=3)
+
+ggsave(filename="_data/corXY_obs_Exp04March2020_TraitNumber.png", width=6, height=5, device="png", dpi = 300, last_plot())
+
+
+
+
+
+