From 1dc7e5d5b769e1cf3b4a0fd3994c54059c292f86 Mon Sep 17 00:00:00 2001 From: Francesco Sabatini <francesco.sabatini@idiv.de> Date: Wed, 15 Apr 2020 10:33:46 +0200 Subject: [PATCH] Aligned 01 --- 01_Mesobromion.R | 148 ++++++++++++++--------------------------------- 1 file changed, 45 insertions(+), 103 deletions(-) diff --git a/01_Mesobromion.R b/01_Mesobromion.R index 9d6bf71..09c4102 100644 --- a/01_Mesobromion.R +++ b/01_Mesobromion.R @@ -169,127 +169,69 @@ traits <- traits %>% #### ## Import output #### - myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN(\\.[A-Za-z]+|_[A-Za-z]+)_[0-9]+.RData", full.names = T) -corXY <- NULL -#corXY.perm <- NULL -#corXY.bootstr <- NULL -for(ff in myfilelist){ - index <- as.numeric(regmatches(ff, gregexpr("[[:digit:]]+", ff))[[1]]) - load(ff) - corXY <- bind_rows(corXY, cor.obs) - #corXY.perm <- rbind(corXY.perm, cor.perm) - #corXY.bootstr <- rbind(corXY.bootstr, cor.bootstr) - print(index) -} - - -dataFiles = map(myfilelist, load) -dat = bind_rows(dataFiles) - - +dataFiles = purrr::map(myfilelist, function(x){get(load(x))}) +corXY = bind_rows(dataFiles) %>% + as.tbl() +rm(corXY.output, dataFiles) -diffs <- corXY.bootstr - corXY.perm -sign <- apply(diffs, MARGIN=1, function(x){sum(x>0)>190}) -sign <- data.frame(Trait.comb=names(sign), sign=sign) - -corXY %>% - arrange(Test, desc(Coef)) - -aa <- data.frame(Trait.comb=paste0("t", 1:ncol(traits)), +trait.labs <- data.frame(Trait.comb=1:ncol(traits), trait.name=colnames(traits))#[-which(colnames(traits) %in% c("species", "species0"))]) -corXY.ci <- data.frame(Trait.comb=rownames(corXY.bootstr), - obs=corXY.bootstr[,1], - q025=apply(corXY.bootstr[,-1], MARGIN=1, quantile, probs=0.025), - q975=apply(corXY.bootstr[,-1], MARGIN=1, quantile, probs=0.975)) %>% - left_join(sign, by="Trait.comb") %>% - as.tbl() %>% +## calculate confidence intervals for bootstrapped correlations +corXY.ci <- corXY %>% + filter(bootstr.n==0) %>% + dplyr::select(-bootstr.n) %>% + group_by(Trait.comb) %>% + slice(1) %>% + left_join(corXY %>% + filter(bootstr.n!=0) %>% + group_by(Trait.comb) %>% + summarize(q025=quantile(Coef.obs, 0.025), + q975=quantile(Coef.obs, 0.975), + greater.than.perm=mean(Coef.obs>Coef.perm), + n=n())) %>% + #calculate significance using permuted correlations + mutate(sign_plus=greater.than.perm>0.975) %>% + mutate(sign_minus=greater.than.perm<0.025) %>% mutate(Trait.comb2=Trait.comb) %>% - separate(Trait.comb2, c("trait1", "trait2")) %>% #, "trait3")) %>% - mutate(trait2=paste0("t", trait2)) %>% - #mutate(trait3=paste0("t", trait3)) %>% - left_join(aa %>% - rename(trait1=Trait.comb), by="trait1") %>% - dplyr::select(-trait1) %>% - rename(trait1=trait.name) %>% - left_join(aa %>% - rename(trait2=Trait.comb), by="trait2") %>% - dplyr::select(-trait2) %>% - rename(trait2=trait.name) %>% - #left_join(aa %>% - # rename(trait3=Trait.comb), by="trait3") %>% - #dplyr::select(-trait3) %>% - #rename(trait3=trait.name) %>% - mutate(ntraits= (!is.na(trait1)) + (!is.na(trait2))) %>% - arrange(ntraits, obs) + #split strings of trait combinations and add labels + separate(Trait.comb2, into=paste0("trait", 1:7)) %>% + mutate_at(.vars=vars(trait1:trait7), + .funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$trait.name)) %>% + rowwise() %>% + mutate(ntraits= length(unlist(strsplit(Trait.comb, "_")))) %>% + ungroup() %>% + arrange(ntraits, Coef.obs) %>% + dplyr::select(Trait.comb, Test, n, ntraits, everything()) -### spider plot -# faceted spider plot -# mydata <- corXY.ci %>% -# filter(grepl(pattern = "t2", Trait.comb)) %>% -# #filter(grepl(pattern = "_6", Trait.comb)) %>% -# #filter(sign==T) %>% -# filter(ntraits<=2) %>% -# mutate(Trait.comb=factor(Trait.comb)) %>% -# mutate(seq=1:nrow(.)) - - - -corXY.SES <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>% - mutate(Trait.comb0=Trait.comb) %>% - separate(Trait.comb0, c("trait1", "trait2")) %>% #, "trait3")) %>% - mutate(trait2=paste0("t", trait2)) %>% - #mutate(trait3=paste0("t", trait3)) %>% - left_join(aa %>% - rename(trait1=Trait.comb), by="trait1") %>% - dplyr::select(-trait1) %>% - rename(trait1=trait.name) %>% - left_join(aa %>% - rename(trait2=Trait.comb), by="trait2") %>% - dplyr::select(-trait2) %>% - rename(trait2=trait.name) %>% -# left_join(aa %>% -# rename(trait3=Trait.comb), by="trait3") %>% -# dplyr::select(-trait3) %>% -# rename(trait3=trait.name) %>% - dplyr::select(Trait.comb, trait1, trait2, q025:conf.p) %>% #trait3, - arrange(desc(SES.np)) -print(corXY.SES, n=20) - - -mydata.SES <- corXY.ci %>% - left_join(corXY.SES %>% - mutate(sign.SES=ifelse(smaller.than.obs>0.975, T, F)) %>% # smaller.than.obs>0.995 | - dplyr::select(Trait.comb, sign.SES), - by="Trait.comb") - - -traits.sign.alone <- mydata.SES %>% +traits.sign.alone <- corXY.ci %>% filter(ntraits==1) %>% - filter(sign.SES==T) %>% + filter(sign_plus==T) %>% pull(trait1) -mydata <- mydata.SES %>% +mydata <- corXY.ci %>% filter(ntraits==1)%>% - mutate(seq=1:nrow(.)) + mutate(seq=1:n()) ggplot(data=mydata) + - geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, col=sign.SES)) + - geom_point(aes(x=obs, y=seq), pch=15) + + geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, + col=sign_plus)) + + geom_point(aes(x=Coef.obs, y=seq), pch=15) + scale_y_continuous(breaks=mydata$seq, labels=mydata$Trait.comb) -mydata <- mydata.SES %>% +mydata <- mydata %>% #filter(ntraits==2)%>% - filter(sign.SES==T) %>% + filter(sign_plus==T) %>% filter(trait1 %in% traits.sign.alone & (is.na(trait2) | trait2 %in% traits.sign.alone)) %>% mutate(seq=1:nrow(.)) # %>% # arrange(desc(obs)) %>% # slice(1:20) ggplot(data=mydata) + - geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, col=sign.SES)) + - geom_point(aes(x=obs, y=seq), pch=15) + + geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, + col=sign_plus)) + + geom_point(aes(x=Coef.obs, y=seq), pch=15) + scale_y_continuous(breaks=mydata$seq, labels=mydata$Trait.comb, name="RV correlation") @@ -297,9 +239,9 @@ ggplot(data=mydata) + -traits.sign <- traits %>% - dplyr::select(species0, traits.sign.alone) -write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t") +#traits.sign <- traits %>% +# dplyr::select(species0, traits.sign.alone) +#write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t") break() -- GitLab