98_SummarizeSimulations.R

## 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)

mypath <- "_data/Experiment_04Mar2020"

myfiles <- list.files(path=mypath, 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))))
  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(ntraits=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("XY", "XY.T", "XY.TR")) %>% 
  group_by(main, ntraits, 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("XY", "XY.T", "XY.TR")) %>% 
              group_by(main, ntraits, 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, ntraits, corr, trait, envir)


outp.summary 

  
get.ntraits <- function(x){ 
  tmp <- str_split(x, pattern = " ")[[1]]
  return(length(tmp))
}
  


## plotting power for XY with corr

ggplot(data=outp.summary %>% 
        ungroup() %>% 
        rowwise() %>% 
        mutate(sel.ntraits=factor(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))) + 
  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
  #scale_colour_brewer(palette = "Dark2") + 
  facet_grid(sel.ntraits~ntraits) + 
  theme_bw() + 
  theme(panel.grid = element_blank())

ggsave(filename="_data/Experiment_04Mar2020/corXY_obs_Exp04March2020.png", width=6, height=5, device="png", dpi = 300, last_plot())

## plotting non-parametric SES for XY with corr

ggplot(data=outp.summary %>% 
         ungroup() %>% 
         rowwise() %>% 
         mutate(sel.ntraits=factor(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))) + 
  geom_line(aes(x=main, y=SES.med, group=trait, col=trait)) + 
  #scale_colour_brewer(palette = "Dark2") + 
  #facet_grid(sel.ntraits~ntraits, scales = "free") + 
  theme_bw() + 
  theme(panel.grid = element_blank())

ggsave(filename="_data/Experiment_04Mar2020/corXY_SES_Exp04March2020.png", width=6, height=5, device="png", dpi = 300, last_plot())

## plotting XE

ggplot(data=outp.summary %>% 
         ungroup() %>% 
         dplyr::filter(stat.type=="XE") %>% 
         #dplyr::filter(envir=="1") %>% 
         filter(corr==0) %>% 
         filter(trait %in% c("1","2","3", "1 2")) %>% 
         mutate(ntraits=as.factor(ntraits))
       #mutate(group0=paste("t", trait, " - e", envir))
) + 
  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
  scale_colour_brewer(palette = "Dark2") + 
  facet_grid(envir~ntraits) + 
  theme_bw() + 
  theme(panel.grid = element_blank())

## plotting XY.T

ggplot(data=outp.summary %>% 
         ungroup() %>% 
         dplyr::filter(stat.type=="XY.T") %>% 
         dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
         mutate(ntraits=as.factor(ntraits))) + 
  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
  scale_colour_brewer(palette = "Dark2") + 
  facet_grid(corr~ntraits) + 
  theme_bw() + 
  theme(panel.grid = element_blank())

## plotting XY.TR

ggplot(data=outp.summary %>% 
         ungroup() %>% 
         dplyr::filter(stat.type=="XY.TR") %>% 
         dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
         mutate(ntraits=as.factor(ntraits))) + 
  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
  scale_colour_brewer(palette = "Dark2") + 
  facet_grid(corr~ntraits) + 
  theme_bw() + 
  theme(panel.grid = element_blank())