From fc416f5757642e479c554756fcddbefdc97663b2 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Thu, 20 Aug 2020 10:51:15 +0200
Subject: [PATCH] Results v2 finalized
---
00_Mesobromion_DataPreparation.R | 2 +-
02_Mesobromion_ExamineOutput.R | 366 +++++++++++++++++++++----------
2 files changed, 250 insertions(+), 118 deletions(-)
diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index 99650c5..9541231 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -330,7 +330,7 @@ write_delim(env, path="_data/Mesobromion/env.v2.10perc.cov.txt", delim="\t")
### create list of relevées
write_delim(species %>%
dplyr::select(RELEVE_NR),
- path="_data/Mesobromion/ReleveList.txt", delim="\t")
+ path="_output/S10_ReleveList.txt", delim="\t")
##check for species without trait info
diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R
index ab6dfe1..b99cbff 100644
--- a/02_Mesobromion_ExamineOutput.R
+++ b/02_Mesobromion_ExamineOutput.R
@@ -8,6 +8,21 @@ library(vegan)
source("99_HIDDEN_functions.R")
source("01b_MesobromionCluster.R")
+## ancillary function to get the best combination at a given Number of interactions
+get.best <- function(x, N, labs){
+ out <- x %>%
+ filter(ntraits==N) %>%
+ filter(sign_plus==TRUE) %>%
+ arrange(desc(Coef.obs)) %>%
+ slice(1) %>%
+ pull(Trait.comb)
+ lab <- labs %>%
+ filter(Trait.comb %in% strsplit(out, split = "_")[[1]]) %>%
+ pull(trait.name)
+ return(list(Trait.comb=out, trait.name=lab))
+}
+
+
####1. Reimport data ################################
traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
traits.sign.cov <- read_delim("_data/Mesobromion/traits.v2.10perc.cov.sign.txt", delim="\t")
@@ -26,30 +41,30 @@ traits <- traits %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0")
-traits.sign.cov <- traits.sign.cov %>%
- as.data.frame() %>%
- # filter(species0 %in% colnames(species)) %>%
- mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
- column_to_rownames("species0")
-
-traits.sign.pa <- traits.sign.pa %>%
- as.data.frame() %>%
- # filter(species0 %in% colnames(species)) %>%
- mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
- column_to_rownames("species0")
-
+# traits.sign.cov <- traits.sign.cov %>%
+# as.data.frame() %>%
+# # filter(species0 %in% colnames(species)) %>%
+# mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
+# column_to_rownames("species0")
+#
+# traits.sign.pa <- traits.sign.pa %>%
+# as.data.frame() %>%
+# # filter(species0 %in% colnames(species)) %>%
+# mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
+# column_to_rownames("species0")
+#
## adapt trait labs to sign traits only
-trait.labs.sign.cov <- trait.labs %>%
- filter(trait.name %in% colnames(traits.sign.cov)) %>%
- arrange(match(trait.name, colnames(traits.sign.cov))) %>%
- rename(Trait.comb.new=Trait.comb) %>%
- mutate(Trait.comb=1:8) %>%
- dplyr::select(Trait.comb, everything(), -Trait.comb.new)
+# trait.labs.sign.cov <- trait.labs %>%
+# filter(trait.name %in% colnames(traits.sign.cov)) %>%
+# arrange(match(trait.name, colnames(traits.sign.cov))) %>%
+# rename(Trait.comb.new=Trait.comb) %>%
+# mutate(Trait.comb=1:8) %>%
+# dplyr::select(Trait.comb, everything(), -Trait.comb.new)
#rename trait based on short labels
colnames(traits) <- trait.labs$Short_english_name
-colnames(traits.sign.cov) <- trait.labs.sign.cov$Short_english_name
+# colnames(traits.sign.cov) <- trait.labs.sign.cov$Short_english_name
##check traits
@@ -111,29 +126,29 @@ trait.coverage <- species.cov %>%
trait.recap <- trait.recap %>%
left_join(trait.coverage, by=c("Code"="trait"))
-write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
+write_csv(trait.recap, path="_output/S3_Traits.recap.csv")
###### 2. Import output from Cluster ####
##### 2.1 Cover values ######
-## traits onebyone
-myfilelist0 <- list.files(path="_derived/Mesobromion/Cover/onebyone/", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T)
-dataFiles0 = purrr::map(myfilelist0, function(x){get(load(x))})
-#load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
-corXY.alone = bind_rows(dataFiles0) %>%
- as.tbl() %>%
- distinct()
-corXY.alone.ci <- get.ci(corXY.alone)
-corXY.alone.ci <- corXY.alone.ci %>%
- mutate(trait1=Trait.comb) %>%
- mutate_at(.vars=vars(trait1),
- .funs=~factor(.,
- levels=trait.labs$Trait.comb,
- labels=trait.labs$Short_english_name)) %>%
- arrange(ntraits, desc(Coef.obs)) %>%
- dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
- mutate(Trait.comb=NA) %>%
- mutate(run="alone")
+# ## traits onebyone
+# myfilelist0 <- list.files(path="_derived/Mesobromion/Cover/onebyone/", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T)
+# dataFiles0 = purrr::map(myfilelist0, function(x){get(load(x))})
+# #load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
+# corXY.alone = bind_rows(dataFiles0) %>%
+# as.tbl() %>%
+# distinct()
+# corXY.alone.ci <- get.ci(corXY.alone)
+# corXY.alone.ci <- corXY.alone.ci %>%
+# mutate(trait1=Trait.comb) %>%
+# mutate_at(.vars=vars(trait1),
+# .funs=~factor(.,
+# levels=trait.labs$Trait.comb,
+# labels=trait.labs$Short_english_name)) %>%
+# arrange(ntraits, desc(Coef.obs)) %>%
+# dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
+# mutate(Trait.comb=NA) %>%
+# mutate(run="alone")
### sequential trait combo
@@ -156,37 +171,17 @@ corXY.all.ci <- corXY.all.ci %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
mutate(run="seq")
-### SEQUENTIAL of sign traits
-# myfilelist <- list.files(path="_derived/Mesobromion/Cover/_test1", pattern="HIDDENcov_round_[0-9]+.RData", full.names = T)
-# dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
-# #load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
-# corXY.comb = bind_rows(dataFiles) %>%
-# as.tbl()
-#
-# corXY.comb.ci <- get.ci(corXY.comb)
-# corXY.comb.ci <- corXY.comb.ci %>%
-# mutate(Trait.comb2=Trait.comb) %>%
-# #split strings of trait combinations and add labels
-# separate(Trait.comb2, into=paste0("trait", 1:8)) %>%
-# mutate_at(.vars=vars(trait1:trait3),
-# .funs=~factor(.,
-# levels=trait.labs$Trait.comb,
-# labels=trait.labs$Short_english_name)) %>%
-# arrange(ntraits, Coef.obs) %>%
-# dplyr::select(Trait.comb, Test, n, ntraits, everything())
-
-
-rm( dataFiles0, dataFiles1)
+rm( dataFiles1) #dataFiles0,
### merge together
-corXY.ci <- corXY.all.ci %>%
- bind_rows(corXY.alone.ci)
+corXY.ci <- corXY.all.ci # %>%
+ #bind_rows(corXY.alone.ci)
#### 2.1.2 One by one - Graph of r(XY) ####
### list traits being significant when taken one by one
traits.sign.alone <- corXY.ci %>%
- filter(run=="alone") %>%
+ #filter(run=="alone") %>%
filter(ntraits==1) %>%
filter(sign_plus==T) %>%
pull(trait1)
@@ -194,7 +189,7 @@ traits.sign.alone <- corXY.ci %>%
### Plot combinations one by one
mydata.byone <- corXY.ci %>%
- filter(run=="alone") %>%
+ #filter(run=="alone") %>%
filter(ntraits==1)%>%
arrange(Coef.obs) %>%
mutate(seq=1:n())
@@ -224,7 +219,7 @@ mydata <- corXY.ci %>%
all_vars(. %in% traits.sign.alone | is.na(.)))%>%
mutate(seq=1:n())
-#### 2.1.3 All traits - (Huge) Graph of r(XY) (w/ combinations) ####
+#### ~2.1.3 All traits - (Huge) Graph of r(XY) (w/ combinations) ####
### filter combinations where all traits belong to traits.sign.alone
# allpredictors <- top.first %+% (mydata) +
@@ -241,20 +236,6 @@ mydata <- corXY.ci %>%
######## 2.1.4 Best - Graph of r(XY) using best combination of traits at each level of interaction N ########
### extract best combinations of traits
-## ancillary function to get the best combination at a given Number of interactions
-get.best <- function(x, N, labs){
- out <- x %>%
- filter(ntraits==N) %>%
- filter(sign_plus==TRUE) %>%
- arrange(desc(Coef.obs)) %>%
- slice(1) %>%
- pull(Trait.comb)
- lab <- labs %>%
- filter(Trait.comb %in% strsplit(out, split = "_")[[1]]) %>%
- pull(trait.name)
- return(list(Trait.comb=out, trait.name=lab))
-}
-
top.one.by.one <- get.best(mydata, N=1, labs=trait.labs)
## Routine to extract the best combination at each level of interaction (up to max traits)
@@ -307,7 +288,7 @@ for(nn in 1:maxtraits){
best.4traits <- corXY.ci %>%
filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
- dplyr::select(trait1:trait7) %>%
+ dplyr::select(trait1:trait8) %>%
mutate_all(~as.character(.))
best.4traits <- as.character(best.4traits[1,1:4])
@@ -333,7 +314,7 @@ mydata.best <- mydata %>%
write_csv(mydata.best %>%
dplyr::select(Trait.comb:sign_plus),
- path = "_derived/Mesobromion/S5_BestSolutionTiers.csv")
+ path = "_output/S5_BestSolutionTiers.cov.csv")
### Graph of all the best combinations with text legend
(top.all <- ggplot(data=(mydata.best %>%
@@ -368,7 +349,7 @@ tt2 <- ttheme_minimal(
# dplyr::select(Code, Trait=trait.name) %>%
# slice(33:n()), theme=tt2),
nrow=1, rel_widths=c(0.70, 0.3)))
-ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_pa.png", dpi=400,
+ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_cov.png", dpi=400,
width=6, height=3, topall.leg)
@@ -379,54 +360,205 @@ break()
####### 2.2.1 Import and adjust dataset ####
### merge output from one-by-one and all combo of sign traits
-### all traits taken alone
-myfilelist0 <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa_[0-9]+_.RData", full.names = T)
-dataFiles0 = purrr::map(myfilelist0, function(x){get(load(x))})
+# ### all traits taken alone
+# myfilelist0 <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa_[0-9]+_.RData", full.names = T)
+# dataFiles0 = purrr::map(myfilelist0, function(x){get(load(x))})
+# #load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
+# corXY.alone = bind_rows(dataFiles0) %>%
+# as.tbl() %>%
+# distinct()
+# corXY.alone.ci <- get.ci(corXY.alone)
+# corXY.alone.ci <- corXY.alone.ci %>%
+# mutate(trait1=Trait.comb) %>%
+# mutate_at(.vars=vars(trait1),
+# .funs=~factor(.,
+# levels=trait.labs$Trait.comb,
+# labels=trait.labs$Short_english_name)) %>%
+# arrange(ntraits, Coef.obs) %>%
+# dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
+# mutate(Trait.comb=NA)
+
+
+### combo of sign traits
+myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa_round_[0-9]+.RData", full.names = T)
+dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
#load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
-corXY.alone = bind_rows(dataFiles0) %>%
+corXY.all = bind_rows(dataFiles) %>%
as.tbl() %>%
distinct()
-corXY.alone.ci <- get.ci(corXY.alone)
-corXY.alone.ci <- corXY.alone.ci %>%
- mutate(trait1=Trait.comb) %>%
- mutate_at(.vars=vars(trait1),
+
+corXY.ci.pa <- get.ci(corXY.all)
+corXY.ci.pa <- corXY.ci.pa %>%
+ mutate(Trait.comb2=Trait.comb) %>%
+ separate(Trait.comb2, into=paste0("trait", 1:14)) %>%
+ mutate_at(.vars=vars(trait1:trait14),
.funs=~factor(.,
levels=trait.labs$Trait.comb,
labels=trait.labs$Short_english_name)) %>%
- arrange(ntraits, Coef.obs) %>%
+ arrange(ntraits, desc(Coef.obs)) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
- mutate(Trait.comb=NA)
+ mutate(run="seq")
+rm( dataFiles1, corXY.all)
-### combo of sign traits
-myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpresabs_[0-9]+_.RData", full.names = T)
-dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
-#load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData")
-corXY.comb = bind_rows(dataFiles) %>%
- as.tbl()
-corXY.comb.ci <- get.ci(corXY.comb)
-corXY.comb.ci <- corXY.comb.ci %>%
- mutate(Trait.comb2=Trait.comb) %>%
- #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.sign$Trait.comb,
- labels=trait.labs.sign$Short_english_name)) %>%
+#### 2.2.2 One by one - Graph of r(XY) ####
+### list traits being significant when taken one by one
+traits.sign.alone.pa <- corXY.ci.pa %>%
+ filter(ntraits==1) %>%
+ filter(sign_plus==T) %>%
+ pull(trait1)
+#[1] Leaf_Scleroph Height GF_Nanophan BL_Dur SLA LeafP LeafCdrymass CSR
+#[9] LeafN SeedLength LDMC GF_Hemiphan
+
+### Plot combinations one by one
+mydata.byone <- corXY.ci.pa %>%
+ filter(ntraits==1)%>%
+ arrange(Coef.obs) %>%
+ mutate(seq=1:n())
+
+(top.first <- ggplot(data=mydata.byone %>%
+ mutate(sign_plus=fct_rev(as.factor(sign_plus)))) +
+ 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.byone$seq,
+ labels=mydata.byone$trait1, name=NULL) +
+ scale_x_continuous(name="RD correlation") +
+ labs(color = "Significant") +
+ theme_bw() +
+ theme(panel.grid.minor = element_blank(),
+ axis.text = element_text(size=7),
+ legend.position = c(0.8, 0.1))
+)
+
+ggsave(filename = "_pics/SXXX_TopFirstPredictors_CI_pa.png", dpi=400,
+ width=4, height=6, plot = top.first)
+
+
+######## 2.2.4 Best - Graph of r(XY) using best combination of traits at each level of interaction N ########
+mydata <- corXY.ci.pa %>%
+ filter_at(.vars=vars(trait1:trait14),
+ all_vars(. %in% traits.sign.alone.pa | is.na(.)))%>%
+ mutate(seq=1:n())
+
+### extract best combinations of traits
+top.one.by.one <- get.best(mydata, N=1, labs=trait.labs)
+
+## Routine to extract the best combination at each level of interaction (up to max traits)
+maxtraits <- 14
+for(nn in 1:maxtraits){
+ if(nn==1) {
+ best.at.1 <- get.best(mydata, N=nn, labs=trait.labs)
+ newdata <- mydata %>%
+ filter_at(.vars=vars(trait1:trait8),
+ .vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.)))
+ new.best.row <- newdata %>%
+ filter(Trait.comb==best.at.1$Trait.comb)
+ upper <- new.best.row$q975
+ lower <- new.best.row$q025
+ print(paste("new best at nn", nn, best.at.1$trait.name))
+ best.progr <- best.at.1$Trait.comb
+ }
+ if(nn>1){
+ better <- list()
+ better$Trait.comb <- newdata %>%
+ filter(ntraits==nn) %>%
+ filter(q025>upper) %>%
+ arrange(desc(Coef.obs)) %>%
+ slice(1) %>%
+ pull(Trait.comb)
+
+ if(length(better$Trait.comb>0)){
+ better$trait.name <- trait.labs %>%
+ filter(Trait.comb %in% strsplit(better$Trait.comb, split = "_")[[1]]) %>%
+ pull(trait.name)
+
+ newdata <- newdata %>%
+ rowwise() %>%
+ mutate(nmatching= sum(unlist(strsplit(Trait.comb, "_")) %in%
+ unlist(strsplit(better$Trait.comb, "_")),
+ na.rm=T)) %>%
+ ungroup() %>%
+ filter(nmatching==nn)
+
+ new.best.row <- newdata %>%
+ filter(Trait.comb==better$Trait.comb)
+ upper <- new.best.row$q975
+ lower <- new.best.row$q025
+ print(paste("new best at nn", nn, paste(better$trait.name, collapse=" ")))
+ best <- better
+ best.progr <- c(best.progr, better$Trait.comb)
+ }
+ }
+}
+
+best.4traits <- corXY.ci %>%
+ filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
+ dplyr::select(trait1:trait8) %>%
+ mutate_all(~as.character(.))
+best.4traits <- as.character(best.4traits[1,1:4])
+
+### Create dataset with best combinations + all the one-way combinations
+mydata.best <- mydata %>%
+ filter_at(.vars=vars(trait1:trait14),
+ # .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>%
+ .vars_predicate = all_vars(. %in% traits.sign.alone.pa | is.na(.))) %>%
+ #filter(ntraits <= length(best$trait.name)) %>%
+ filter(ntraits>1) %>%
+ filter(sign_plus==T) %>%
arrange(ntraits, Coef.obs) %>%
- dplyr::select(Trait.comb, Test, n, ntraits, everything())
+ group_by(ntraits) %>%
+ slice(n()) %>%
+ ungroup() %>%
+ bind_rows(corXY.ci.pa %>%
+ filter(run=="seq") %>%
+ filter(ntraits==1) %>%
+ filter(trait1 %in% traits.sign.alone.pa)) %>%
+ arrange(ntraits, Coef.obs) %>%
+ mutate(seq=1:n()) %>%
+ mutate(sign_plus=factor(Trait.comb %in% best.progr))
+
+write_csv(mydata.best %>%
+ dplyr::select(Trait.comb:sign_plus),
+ path = "_output/SXXX_BestSolutionTiers.pa.csv")
+
+### Graph of all the best combinations with text legend
+(top.all <- ggplot(data=(mydata.best %>%
+ mutate(size0=.6+(as.numeric(sign_plus)-1)*.6))) +
+ geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, col="a",
+ lwd=size0)) +
+ geom_point(aes(x=Coef.obs, y=seq), pch=15) +
+ scale_y_continuous(breaks=mydata.best$seq,
+ labels=mydata.best$Trait.comb, name=NULL) +
+ scale_x_continuous(name="RD correlation") +
+ scale_size_identity() +
+ theme_bw() +
+ theme(panel.grid.minor = element_blank(),
+ axis.text = element_text(size=7),
+ legend.position = "none"))
+
+
+tt2 <- ttheme_minimal(
+ core = list(fg_params=list(cex = .7),
+ padding=unit(c(1, 1), "mm")),
+ colhead = list(fg_params=list(cex = .7)),
+ rowhead = list(fg_params=list(cex = .7)))
+
+(topall.leg <- cowplot::plot_grid(top.all,
+ tableGrob(trait.labs %>%
+ mutate(Code=1:n()) %>%
+ dplyr::select(Code, Trait=Short_english_name) %>%
+ filter(Trait %in% traits.sign.alone.pa),
+ theme=tt2, rows = NULL),
+ nrow=1, rel_widths=c(0.70, 0.3)))
+ggsave(filename = "_pics/SXXX_Best_AllCombinations_CI_pa.png", dpi=400,
+ width=6, height=3, topall.leg)
+
-rm( dataFiles, dataFiles0)
-### merge together
-corXY.ci <- corXY.comb.ci %>%
- bind_rows(corXY.alone.ci %>%
- filter(!trait1 %in% trait.labs.sign$Short_english_name)) %>%
- arrange(Coef.obs)
-#### ____TO DO ########
@@ -659,7 +791,7 @@ cor(env %>%
fuzz.cor <- cor(pca.fuzz$CA$u[,1:3], scores.pca[,1:5])
-sink("_derived/Mesobromion/EnvFit_CWMs_env.txt")
+sink("_output/S9_EnvFit_CWMs_env.txt")
cwms.cor
env.cor
fuzz.cor
@@ -880,7 +1012,7 @@ ggsave("_pics/S6_PCA_Traits_1-4_only7.png", width=10, height=5, dpi=300, PC_trai
##### 4.3b Alternative version of figS6, showing the species ####
tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
- mutate(species0=traits$species0) %>%
+ mutate(species0=rownames(traits)) %>%
mutate(species=species0) %>%
separate(species0, sep="_", into=c("Gen", "Spe")) %>%
mutate(Gen=substr(Gen, 1, 3)) %>%
@@ -924,7 +1056,7 @@ ggsave("_pics/S6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.
## Create list of species labels
write_delim(tmp %>%
dplyr::select(species, label=labels),
- path = "_derived/Mesobromion/SpeciesList.txt")
+ path = "_output/S11_SpeciesList.txt")
--
GitLab