diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index 391b9a4d33ff2847851acd95f6d1fd5fc9783c67..69b1d1d7f0f667415ae20cc685c2d53ed8a94b91 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -21,6 +21,7 @@ traits0 <- read_delim("_data/Mesobromion/traits4.csv", delim =",", ## manually
dplyr::select(!starts_with("ZWT")) %>%
dplyr::select(!starts_with("Fabaceae")) %>%
dplyr::select(!starts_with("BEST_")) %>%
+ dplyr::select(!ends_with("1")) %>%
dplyr::select(- any_of(c("PR_STAT_Indigen", "PR_STAT_Neophyt", "PR_STAT_Archaeophyt", "URBAN_urbanophob", "URBAN_maessig_urbanophob",
"URBAN_urbanoneutral", "URBAN_maessig_urbanophil", "URBAN_urbanophil", "WUH_von", "WUH_bis", "ARL_c_I_von", "ARL_c_I_bis",
"BL_ANAT_hydromorph", "WEID_V", "MAHD_V", "TRIT_V"))) %>% #empty trait
@@ -33,7 +34,7 @@ traits0 <- read_delim("_data/Mesobromion/traits4.csv", delim =",", ## manually
mutate(species=gsub(pattern=" $", replacement = "", x = species)) %>%
mutate(species=ifelse(is.na(word(species, 1, 2)), species, word(species, 1, 2))) %>%
ungroup()
-dim(traits0) #902 obs. of 67 variables:
+dim(traits0) #902 obs. of 65 variables:
## remove species with NAs
species.to.remove <- traits0 %>%
@@ -96,8 +97,8 @@ all.traits <- traits0 %>%
by="species")
traits <- all.traits %>%
filter(!is.na(LeafArea))
-dim(all.traits) #[1] 898 82
-dim(traits) #[1] 799 82
+dim(all.traits) #[1] 898 80
+dim(traits) #[1] 799 80
@@ -252,12 +253,13 @@ recode.traits <- function(x){
}
traits <- recode.traits(traits)
-
-### ordered factors
+##exclude traits being all == 0
+traits <- traits %>%
+ dplyr::select(-colnames(.)[which(colSums(.!=0)==0)])
dim(species) #581 509
-dim(traits) #509 54
+dim(traits) #509 49
dim(env) #581 8
@@ -359,7 +361,14 @@ traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
myfilelist <- list.files(path="_derived/Mesobromion/Cover", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
- as_tibble()
+ as_tibble() %>%
+ distinct()
+##find missing
+#missing <- gsub(pattern="_derived/Mesobromion/Cover/HIDDENcov-nona2_", replacement = "", x = myfilelist)
+#missing <- gsub(pattern="_.RData", replacement = "", x = missing)
+#missing <- sort(as.numeric(missing))
+#which(!1:127 %in% missing)
+
rm( dataFiles)
trait.labs <- data.frame(Trait.comb=as.character(1:(ncol(traits)-1)),
@@ -377,7 +386,7 @@ traits.sign.alone <- corXY.ci %>%
traits.sign <- traits %>%
dplyr::select(species0, any_of(traits.sign.alone))
-#write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.cov.sign.txt", delim="\t")
+write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.cov.sign.txt", delim="\t")
### COV - NONAs - all species with at least 1 NAs in traits excluded BEFORE The analysis
" Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
@@ -392,16 +401,17 @@ traits.sign <- traits %>%
8 35 0.241 0.128 0.217 0.289 0.970 999 TRUE FALSE 1 LeafP "
### COV - NONA2 - modified Matrix.x to deal with NAs INSIDE the analysis
-" Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
- <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 36 0.304 0.134 0.281 0.356 0.990 199 TRUE FALSE 1 PlantHeight
- 2 50 0.290 0.126 0.261 0.339 0.985 199 TRUE FALSE 1 BL_ANAT
- 3 2 0.252 0.0621 0.214 0.303 0.985 199 TRUE FALSE 1 LEB_F_Nanophaneroph…
- 4 20 0.251 0.145 0.206 0.320 0.975 199 TRUE FALSE 1 V_VER_Fragmentation
- 5 30 0.250 0.137 0.228 0.298 0.970 199 TRUE FALSE 1 BL_DAU
- 6 32 0.247 0.0932 0.225 0.292 0.965 199 TRUE FALSE 1 SLA
- 7 49 0.236 0.0784 0.212 0.284 0.975 199 TRUE FALSE 1 STRAT_T
- 8 35 0.241 0.127 0.218 0.288 0.940 199 FALSE FALSE 1 LeafP "
+"# A tibble: 53 x 11
+ Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+ <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
+ 1 36 0.304 0.151 0.280 0.355 0.998 999 TRUE FALSE 1 PlantHeight
+ 2 50 0.292 0.0641 0.265 0.346 0.988 999 TRUE FALSE 1 BL_ANAT
+ 3 2 0.262 0.0633 0.223 0.312 0.989 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 4 30 0.259 0.145 0.230 0.307 0.976 999 TRUE FALSE 1 BL_DAU
+ 5 20 0.251 0.109 0.205 0.324 0.984 999 TRUE FALSE 1 V_VER_Fragmentation
+ 6 32 0.247 0.238 0.224 0.297 0.964 999 TRUE FALSE 1 SLA
+ 7 35 0.241 0.200 0.221 0.287 0.965 999 TRUE FALSE 1 LeafP
+ 8 49 0.238 0.0786 0.217 0.286 0.965 999 TRUE FALSE 1 STRAT_T "
### COV - without deleting NAs
"# A tibble: 53 x 11
@@ -413,7 +423,7 @@ traits.sign <- traits %>%
#### 2. Traits individually significant for Presence|absence data####
traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
-myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa_[0-9]+_.RData", full.names = T)
+myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa-nona2_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
as_tibble()
@@ -434,7 +444,7 @@ traits.sign.alone <- corXY.ci %>%
traits.sign <- traits %>%
dplyr::select(species0, any_of(traits.sign.alone))
-#write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t")
+write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.pa.sign.txt", delim="\t")
## Pres Abs No Nas
"# A tibble: 13 x 11
@@ -454,6 +464,25 @@ traits.sign <- traits %>%
12 9 0.217 0.110 0.179 0.272 0.968 999 TRUE FALSE 1 LEB_F_Chamaephyt
13 31 0.211 0.0831 0.190 0.261 0.958 999 TRUE FALSE 1 LeafArea "
+## Pres Abs NoNa2 -modified Matrix.x to deal with NAs INSIDE the analysis
+## Data Filled up by HB used in this analysis
+" Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+ <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
+ 1 50 0.311 0.0624 0.284 0.362 0.996 999 TRUE FALSE 1 BL_ANAT
+ 2 36 0.307 0.190 0.279 0.359 0.995 999 TRUE FALSE 1 PlantHeight
+ 3 2 0.295 0.105 0.254 0.346 0.997 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 4 30 0.260 0.124 0.229 0.312 0.991 999 TRUE FALSE 1 BL_DAU
+ 5 32 0.260 0.243 0.237 0.305 0.987 999 TRUE FALSE 1 SLA
+ 6 35 0.253 0.131 0.230 0.299 0.983 999 TRUE FALSE 1 LeafP
+ 7 33 0.240 0.0910 0.211 0.296 0.979 999 TRUE FALSE 1 LeafC.perdrymass
+ 8 49 0.240 0.0954 0.212 0.287 0.975 999 TRUE FALSE 1 STRAT_T
+ 9 34 0.233 0.294 0.210 0.281 0.963 999 TRUE FALSE 1 LeafN
+10 38 0.228 0.0907 0.203 0.278 0.958 999 TRUE FALSE 1 Seed.length
+11 39 0.226 0.106 0.198 0.282 0.963 999 TRUE FALSE 1 LDMC
+12 5 0.226 0.103 0.192 0.281 0.975 999 TRUE FALSE 1 LEB_F_Hemiphanerophyt
+13 1 0.202 0.0560 0.170 0.248 0.953 999 TRUE FALSE 1 LEB_F_Makrophanerophyt
+14 9 0.196 0.127 0.162 0.249 0.954 999 TRUE FALSE 1 LEB_F_Chamaephyt "
+
## Pres Abs - Without excluding species with NA in traits
"# A tibble: 7 x 11
diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R
index c5f7af2c46044e4afbf7e908d86d699aa10d154c..83d0d276ad9a961902b46b7150e9df95f576fd91 100644
--- a/02_Mesobromion_ExamineOutput.R
+++ b/02_Mesobromion_ExamineOutput.R
@@ -9,11 +9,15 @@ source("99_HIDDEN_functions.R")
source("01b_MesobromionCluster.R")
####1. Reimport data ################################
-traits <- read_delim("_data/Mesobromion/traits.out.10perc.txt", delim="\t")
-traits.sign <- read_delim("_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t")
-trait.labs <- read_delim("_data/Mesobromion/TraitLabels_Long.csv", delim=",") %>%
- rownames_to_column("Trait.comb")
-env <- read_delim("_data/Mesobromion/env.10perc.txt", delim="\t")
+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")
+traits.sign.pa <- read_delim("_data/Mesobromion/traits.v2.10perc.pa.sign.txt", delim="\t")
+trait.labs <- data.frame(trait.name=colnames(traits)[-1]) %>%
+ left_join(read_delim("_data/Mesobromion/TraitLabels_Long.csv", delim=","),
+ by="trait.name") %>%
+ rownames_to_column("Trait.comb") %>%
+ mutate_at(.vars=vars(Short_english_name:Long_English_name), ~ifelse(is.na(.), trait.name, {.}))
+env <- read_delim("_data/Mesobromion/env.v2.10perc.cov.txt", delim="\t")
traits <- traits %>%
as.data.frame() %>%
@@ -21,23 +25,30 @@ traits <- traits %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0")
-traits.sign <- traits.sign %>%
+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 <- trait.labs %>%
- filter(trait.name %in% colnames(traits.sign)) %>%
- arrange(match(trait.name, colnames(traits.sign))) %>%
+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:7) %>%
+ 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) <- trait.labs.sign$Short_english_name
+colnames(traits.sign.cov) <- trait.labs.sign.cov$Short_english_name
##check traits
@@ -45,7 +56,8 @@ is.binary <- function(x){(all(na.omit(x) %in% 0:1))}
####1.1 Table S3 - Trait recap #####
trait.recap <- traits %>%
- mutate_if(.predicate = ~is.numeric(.), ~round(.,3)) %>%
+ dplyr::select(-ends_with("1")) %>% ## exclude two traits that shouldn't be there V_VER_present1 & V_VER_absent1
+ mutate_if(.predicate = ~is.numeric(.), .funs=~round(.,3)) %>%
summarize_all(.funs = list(xxxType.of.variable=~ifelse(is.binary(.), "binary",
ifelse(is.ordered(.), "ordered",
ifelse(is.numeric(.), "quantitative",
@@ -74,12 +86,9 @@ write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
###### 2. Import output from Cluster ####
-###### 2.1 Presence/Absence data ####
-####### 2.1.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)
+##### 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) %>%
@@ -92,50 +101,74 @@ corXY.alone.ci <- corXY.alone.ci %>%
.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(Trait.comb=NA) %>%
+ mutate(run="alone")
+
-### 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))})
+### sequential trait combo
+myfilelist1 <- list.files(path="_derived/Mesobromion/Cover/", pattern="HIDDENcov_round_[0-9]+.RData", full.names = T)
+dataFiles1 = purrr::map(myfilelist1, 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 %>%
+corXY.all = bind_rows(dataFiles1) %>%
+ as.tbl() %>%
+ distinct()
+corXY.all.ci <- get.ci(corXY.all)
+corXY.all.ci <- corXY.all.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),
+ separate(Trait.comb2, into=paste0("trait", 1:8)) %>%
+ mutate_at(.vars=vars(trait1:trait8),
.funs=~factor(.,
- levels=trait.labs.sign$Trait.comb,
- labels=trait.labs.sign$Short_english_name)) %>%
- arrange(ntraits, Coef.obs) %>%
- dplyr::select(Trait.comb, Test, n, ntraits, everything())
+ levels=trait.labs$Trait.comb,
+ labels=trait.labs$Short_english_name)) %>%
+ arrange(ntraits, desc(Coef.obs)) %>%
+ #filter(ntraits>1) %>%
+ 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( dataFiles, dataFiles0)
+rm( dataFiles0, dataFiles1)
### merge together
-corXY.ci <- corXY.comb.ci %>%
- bind_rows(corXY.alone.ci %>%
- filter(!trait1 %in% trait.labs.sign$Short_english_name)) %>%
- arrange(Coef.obs)
-
+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(ntraits==1) %>%
filter(sign_plus==T) %>%
pull(trait1)
+#"Height Leaf_Anat GF_Nanophan BL_Dur VP_Fragm SLA LeafP CSR "
### Plot combinations one by one
mydata.byone <- corXY.ci %>%
+ filter(trait1 != c("V_VER_absent1", "V_VER_present1")) %>% ##delete traits that shouldn't be there
+ filter(run=="alone") %>%
filter(ntraits==1)%>%
+ arrange(Coef.obs) %>%
mutate(seq=1:n())
(top.first <- ggplot(data=mydata.byone %>%
@@ -154,16 +187,18 @@ mydata.byone <- corXY.ci %>%
)
-ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI_pa.png", dpi=400,
+ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI_cov.png", dpi=400,
width=4, height=6, plot = top.first)
-#### 2.1.3 All traits - (Huge) Graph of r(XY) (w/ combinations) ####
-### filter combinations where all traits belong to traits.sign.alone
mydata <- corXY.ci %>%
- filter_at(.vars=vars(trait1:trait7),
+ filter(run=="seq") %>%
+ filter_at(.vars=vars(trait1:trait8),
all_vars(. %in% traits.sign.alone | is.na(.)))%>%
mutate(seq=1:n())
+#### 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) +
# scale_y_continuous(breaks=mydata$seq,
# labels=NULL) +
@@ -192,15 +227,15 @@ get.best <- function(x, N, labs){
return(list(Trait.comb=out, trait.name=lab))
}
-top.one.by.one <- get.best(mydata, N=1, labs=trait.labs.sign)
+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 <- 7
+maxtraits <- 8
for(nn in 1:maxtraits){
if(nn==1) {
- best.at.1 <- get.best(mydata, N=nn, labs=trait.labs.sign)
+ best.at.1 <- get.best(mydata, N=nn, labs=trait.labs)
newdata <- mydata %>%
- filter_at(.vars=vars(trait1:trait7),
+ 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)
@@ -219,7 +254,7 @@ for(nn in 1:maxtraits){
pull(Trait.comb)
if(length(better$Trait.comb>0)){
- better$trait.name <- trait.labs.sign %>%
+ better$trait.name <- trait.labs %>%
filter(Trait.comb %in% strsplit(better$Trait.comb, split = "_")[[1]]) %>%
pull(trait.name)
@@ -242,28 +277,28 @@ for(nn in 1:maxtraits){
}
}
-best.2traits <- corXY.ci %>%
+best.4traits <- corXY.ci %>%
filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
dplyr::select(trait1:trait7) %>%
mutate_all(~as.character(.))
-best.2traits <- as.character(best.2traits[1,1:2])
+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:trait7),
+ filter_at(.vars=vars(trait1:trait8),
# .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>%
.vars_predicate = all_vars(. %in% traits.sign.alone | is.na(.))) %>%
#filter(ntraits <= length(best$trait.name)) %>%
+ filter(ntraits>1) %>%
filter(sign_plus==T) %>%
arrange(ntraits, Coef.obs) %>%
group_by(ntraits) %>%
slice(n()) %>%
ungroup() %>%
- bind_rows(mydata %>%
- filter(ntraits==1) %>%
- filter(trait1 %in% traits.sign.alone) %>%
- arrange(Coef.obs) %>%
- slice(1:n()-1)) %>%
+ bind_rows(corXY.ci %>%
+ filter(run=="seq") %>%
+ filter(ntraits==1) %>%
+ filter(trait1 %in% traits.sign.alone)) %>%
arrange(ntraits, Coef.obs) %>%
mutate(seq=1:n()) %>%
mutate(sign_plus=factor(Trait.comb %in% best.progr))
@@ -295,7 +330,7 @@ tt2 <- ttheme_minimal(
rowhead = list(fg_params=list(cex = .7)))
(topall.leg <- cowplot::plot_grid(top.all,
- tableGrob(trait.labs.sign %>%
+ tableGrob(trait.labs %>%
mutate(Code=1:n()) %>%
dplyr::select(Code, Trait=Short_english_name) %>%
filter(Trait %in% traits.sign.alone),
@@ -312,11 +347,12 @@ ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_pa.png", dpi=400,
break()
-##### 2.2 Cover values ######
+###### 2.2 Presence/Absence data ####
+####### 2.2.1 Import and adjust dataset ####
+### merge output from one-by-one and all combo of sign traits
-#### ____TO DO ########
### all traits taken alone
-myfilelist0 <- list.files(path="_derived/Mesobromion/Cover-pa/", pattern="HIDDENcov.nona_[0-9]+_.RData", full.names = T)
+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) %>%
@@ -330,11 +366,12 @@ corXY.alone.ci <- corXY.alone.ci %>%
levels=trait.labs$Trait.comb,
labels=trait.labs$Short_english_name)) %>%
arrange(ntraits, Coef.obs) %>%
- dplyr::select(Trait.comb, Test, n, ntraits, everything())# %>%
+ dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
mutate(Trait.comb=NA)
+
### combo of sign traits
-myfilelist <- list.files(path="_derived/Mesobromion/Cover/_test1", pattern="HIDDENcov_round_[0-9]+.RData", full.names = T)
+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) %>%
@@ -344,11 +381,11 @@ 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:3)) %>%
- mutate_at(.vars=vars(trait1:trait3),
+ separate(Trait.comb2, into=paste0("trait", 1:7)) %>%
+ mutate_at(.vars=vars(trait1:trait7),
.funs=~factor(.,
- levels=trait.labs$Trait.comb,
- labels=trait.labs$Short_english_name)) %>%
+ levels=trait.labs.sign$Trait.comb,
+ labels=trait.labs.sign$Short_english_name)) %>%
arrange(ntraits, Coef.obs) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything())
@@ -361,24 +398,14 @@ corXY.ci <- corXY.comb.ci %>%
filter(!trait1 %in% trait.labs.sign$Short_english_name)) %>%
arrange(Coef.obs)
+#### ____TO DO ########
##### 3. Calculate CWM ####
#species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
-species.cov <- read_delim("_data/Mesobromion/species.out.10perc.cov.txt", delim="\t")
-traits.cov <- read_delim("_data/Mesobromion/traits.out.10perc.cov.txt", delim="\t")
-
-traits.cov <- traits.cov %>%
- as.data.frame() %>%
- mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
- column_to_rownames("species0")
-
-#rename trait based on short labels
-colnames(traits.cov) <- trait.labs$Short_english_name
-
-
+species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
CWM.wide <- species.cov %>%
pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "rel.cov") %>%
@@ -416,21 +443,22 @@ circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
}
dat <- circleFun(diameter = 2, npoints = 100)
### 4.0 PCA of X (Fuzzy-weighted) matrix + CWMs ####
+source("01b_MesobromionCluster.R")
W.fuzzy <- matrix.x(comm=species.cov %>%
column_to_rownames("RELEVE_NR"),
- traits=traits.cov %>%
+ traits=traits %>%
rownames_to_column("species0") %>%
filter(species0 %in% colnames(species.cov)) %>%
column_to_rownames("species0") %>%
dplyr::select(all_of(traits.sign.alone)),
scale = T)$matrix.X
-pca.fuzz <- rda(W.fuzzy,scale = T) ##### DOUBLE CHECK HERE: IT SHOULDN'T be scaled!
+pca.fuzz <- rda(W.fuzzy)
varexpl <- round((pca.fuzz$CA$eig)/sum(pca.fuzz$CA$eig)*100,1)
scores.pca <- pca.fuzz$CA$u
-cwm.cor <- cor(CWM.wide, scores.pca[,1:4], use = "pairwise.complete.obs") #double check
+cwm.cor <- cor(CWM.wide, scores.pca[,1:3], use = "pairwise.complete.obs") #double check
env.cor <- cor(env %>%
dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC),
- scores.pca[,1:4]) #double check
+ scores.pca[,1:3]) #double check
myvectors <- as.data.frame(env.cor) %>%
rownames_to_column("mylab") %>%
@@ -438,7 +466,7 @@ myvectors <- as.data.frame(env.cor) %>%
bind_rows(as.data.frame(cwm.cor) %>%
rownames_to_column("mylab") %>%
mutate(category="Trait")) %>%
- mutate(fontface0=ifelse(mylab %in% best.2traits, "bold", "italic")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
mutate(mycol=ifelse(category=="Trait", oilgreen, orange))
@@ -450,7 +478,7 @@ basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
theme(panel.grid = element_blank())
PCA.fuzz1_2 <- basemap0 +
- geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:4]*5),
+ geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:3]*5),
aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
geom_segment(data=myvectors,
aes(x=0, xend=PC1, y=0, yend=PC2, col=mycol),
@@ -463,25 +491,25 @@ PCA.fuzz1_2 <- basemap0 +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
-PCA.fuzz3_4 <- basemap0 +
- geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:4]*5),
- aes(x=PC3, y=PC4), pch="+", size=2, alpha=0.8) +
+PCA.fuzz2_3 <- basemap0 +
+ geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:3]*5),
+ aes(x=PC2, y=PC3), pch="+", size=2, alpha=0.8) +
geom_segment(data=myvectors,
- aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
+ aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
geom_label_repel(data=myvectors,# %>%
#mutate_if(~is.numeric(.), ~(.)*1.2),
- aes(x=PC3, y=PC4, label=mylab, col=mycol, fontface=fontface0), size=2,
+ aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
- xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
- ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
+ xlab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ ylab(paste("PC3 (", varexpl[3], "%)", sep=""))
-PC_fuzzy <- cowplot::plot_grid(PCA.fuzz1_2,PCA.fuzz3_4, nrow=1)
-ggsave("_pics/FigSXXX_PC_Fuzzy_1-4.png", width=10, height=5, dpi=300, last_plot())
+PC_fuzzy <- cowplot::plot_grid(PCA.fuzz1_2,PCA.fuzz2_3, nrow=1)
+ggsave("_pics/FigSXXX_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot())
#### 4.0b Alternative showing species scores ####
-tmp <- as.data.frame(pca.fuzz$CA$v[,1:4]*7) %>%
+tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>%
mutate(species0=colnames(species.cov)[-1]) %>%
mutate(species=species0) %>%
separate(species0, sep="_", into=c("Gen", "Spe")) %>%
@@ -502,21 +530,21 @@ PCAfuzz1_2.sp <- basemap0 %+% tmp +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
-PCAfuzz3_4.sp <- basemap0 %+% tmp +
- geom_text(aes(x=PC3, y=PC4, label=labels), size=2, alpha=0.8) +
+PCAfuzz2_3.sp <- basemap0 %+% tmp +
+ geom_text(aes(x=PC2, y=PC3, label=labels), size=2, alpha=0.8) +
geom_segment(data=myvectors,
- aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
+ aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
geom_label_repel(data=myvectors,
- aes(x=PC3, y=PC4, label=mylab, col=mycol, fontface=fontface0), size=2,
+ aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2,
position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7), segment.size = 0.5) +
- xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
- ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
+ xlab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ ylab(paste("PC3 (", varexpl[3], "%)", sep=""))
ggsave("_pics/FigSXXXa_PCA_Fuzzy_1-2_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_2.sp)
-ggsave("_pics/FigSXXXb_PCA_Fuzzy_3-4_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz3_4.sp)
+ggsave("_pics/FigSXXXb_PCA_Fuzzy_2-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz2_3.sp)
@@ -551,7 +579,7 @@ env.cor <- arrow_heads * sqrt (r2)
cor(env %>%
dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), scores.pca[,1:5]) #double check
-fuzz.cor <- cor(pca.fuzz$CA$u[,1:4], scores.pca[,1:5])
+fuzz.cor <- cor(pca.fuzz$CA$u[,1:3], scores.pca[,1:5])
sink("_derived/Mesobromion/EnvFit_CWMs_env.txt")
cwms.cor
@@ -569,7 +597,7 @@ myvectors <- as.data.frame(env.cor) %>%
bind_rows(as.data.frame(fuzz.cor) %>%
rownames_to_column("mylab") %>%
mutate(category="Fuzzy-Weighted")) %>%
- mutate(fontface0=ifelse(mylab %in% best.2traits, "bold", "italic")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol))
@@ -685,7 +713,8 @@ ggsave("_pics/Fig6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PC
# assign weights to traits [Growth forms and root structures are coded as binary, but should be downweighted]
#myweights <- c(rep(1/9, 9), rep(1/17, 17), rep(1, 26))
#traits.gowdis <- FD::gowdis(traits, w=myweights, ord="podani")
-data2 <- traits.cov %>%
+data2 <- traits %>%
+ dplyr::select_if(.predicate = ~is.numeric(.)) %>% ### CAUTION -- ONLY NUMERIC TRAITS SHOWN
dplyr::select(any_of(traits.sign.alone)) %>%
rownames_to_column("species") %>%
filter(species %in% colnames(species.cov)) %>%
@@ -721,7 +750,7 @@ PCA.t1 <- baseplot +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.2traits, "bold", "italic")),
+ mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
aes(x=Comp.1, y=Comp.2, label=Trait, col=oilgreen, fontface=fontface0), size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
#ggtitle("PCoA of species-level traits")
@@ -735,7 +764,7 @@ PCA.t2 <- baseplot +
geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
#geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.2traits, "bold", "italic")),
+ mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
aes(x=Comp.3, y=Comp.4, label=Trait, col=oilgreen, fontface=fontface0), size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
@@ -763,7 +792,7 @@ PCA.t1.sp <- baseplot %+% tmp +
geom_path(data=dat,aes(x,y), col=gray(0.8), lwd=0.5) + #add correlation circle
geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.2traits, "bold", "italic")),
+ mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
aes(x=Comp.1, y=Comp.2, label=Trait, col=oilgreen, fontface=fontface0), size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
@@ -776,7 +805,7 @@ PCA.t2.sp <- baseplot %+% tmp +
geom_label_repel(data=as.data.frame(pca2$loadings[,1:5]) %>%
#geom_label(data=as.data.frame(pca2$loadings[,1:5]*8) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.2traits, "bold", "italic")),
+ mutate(fontface0=ifelse(Trait %in% best.4traits, "bold", "italic")),
aes(x=Comp.3, y=Comp.4, label=Trait, col=oilgreen, fontface=fontface0), size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
@@ -838,7 +867,8 @@ env.all.sf <- SpatialPointsDataFrame(coords=env.all %>% dplyr::select(LON, LAT),
theme_bw()
)
-ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basemap2)
+ggsave(filename="_pics/PlotDistribution.png", height
+ =8, width=6, dpi=300, basemap2)
diff --git a/session.R b/session.R
index 08d6e08f4dd5e416abc27471f51dcda05a9f3991..d524c2515d23d3877881c1e76495d859163c6bd5 100644
--- a/session.R
+++ b/session.R
@@ -11,7 +11,7 @@ chunkn <- 53
combinations <- "all"
start.round <- 1
relax.round <- 1
-exclude.na <- T
+exclude.na <- F
source("01b_MesobromionCluster.R")
Mesobromion(species.path, traits.path, output, myfunction,