diff --git a/98_SummarizeSimulations.R b/98_SummarizeSimulations.R
new file mode 100644
index 0000000000000000000000000000000000000000..68b7dc703557d551a5cfacaafb957d0f504b3c7e
--- /dev/null
+++ b/98_SummarizeSimulations.R
@@ -0,0 +1,77 @@
+## 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)
+
+
+myfiles <- list.files(pattern = "Summary.txt", recursive = T)
+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))))/10
+  tmp <- read_delim(ff, 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(feedb=iter[[3]]) %>% 
+    dplyr::select(main:feedb, everything())
+  output <- bind_rows(output, tmp)
+}
+
+outp.summary <- output %>% 
+  filter(!stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
+  group_by(main, inter, feedb, 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 %>% 
+              filter(stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
+              group_by(main, inter, feedb, 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, feedb, trait, envir)
+
+
+## plotting XY
+
+ggplot(data=outp.summary %>% 
+         ungroup() %>% 
+         filter(stat.type=="XY") %>% 
+         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") + 
+  facet_grid(feedb~inter) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank())
+
+
+## plotting XE
+
+ggplot(data=outp.summary %>% 
+         ungroup() %>% 
+         filter(stat.type=="XE") %>% 
+         filter(envir=="1 ") %>% 
+         mutate(inter=as.factor(inter))
+       #mutate(group0=paste("t", trait, " - e", envir))
+) + 
+  geom_line(aes(x=main, y=power, group=trait, col=trait), lwd=1.5, alpha=0.3) + 
+  scale_colour_brewer(palette = "Dark2") + 
+  facet_grid(feedb~inter) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank())
+
+