diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index 54c33b0a0b17849abcac3eb0b7b2d4b4fb53ab12..d999dad447a3ec3f0119337d3d24f1defa32a113 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -10,10 +10,18 @@ source("99_HIDDEN_functions.R")
##### PART 1 ####
#### 1. traits data ####
-traits0 <- read_delim("_data/Mesobromion/traits3.txt", delim =";", col_names = T, locale = locale(encoding = 'latin1')) %>%
- dplyr::select(- c("PR_STAT_Indigen", "PR_STAT_Neophyt", "PR_STAT_Archaeophyt", "URBAN_urbanophob", "URBAN_maessig_urbanophob",
+traits0 <- read_delim("_data/Mesobromion/traits3_cov.txt", delim =";", ## manually corrected vowels with umlaut
+ col_names = T, locale = locale(encoding = 'UTF-8')) %>%
+ column_to_rownames("X1") %>%
+ dplyr::select(colnames(.)[which(colSums(., na.rm=T)!=0)]) %>%
+ dplyr::select(!starts_with("BLM")) %>%
+ dplyr::select(!starts_with("ZWT")) %>%
+ dplyr::select(!starts_with("Fabaceae")) %>%
+ dplyr::select(!starts_with("BEST_")) %>%
+ 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")) %>% #empty trait
+ "BL_ANAT_hydromorph", "WEID_V", "MAHD_V", "TRIT_V"))) %>% #empty trait
+ rownames_to_column("species") %>%
mutate(species0=species) %>%
rowwise() %>%
# quick and dirty clean up names
@@ -27,8 +35,8 @@ dim(traits0) #907 obs. of 75 variables:
#keep only traits with >=88 completeness
traits0 <- traits0 %>%
- dplyr::select_if(~mean(!is.na(.)) >= 0.88) # 907 x 67
-
+ dplyr::select_if(~mean(!is.na(.)) >= 0.88)
+dim(traits0)# 902 x 67
@@ -43,11 +51,13 @@ traits.asym.binary <- c('LEB_F_Makrophanerophyt','LEB_F_Nanophanerophyt',
'V_VER_Innovationsknospe.mit.Speicherwurzel','V_VER_Ausläuferknolle',
'V_VER_Brutsprösschen','V_VER_Fragmentation','V_VER_Turio','V_VER_Sprossknolle',
'V_VER_phyllogener_Spross','V_VER_Rhizompleiokorm','V_VER_Zwiebel',
- 'V_VER_Ausläuferrhizom','V_VER_Ausläuferzwiebel','V_VER_Bulbille',
+ 'V_VER_Ausläuferrhizom',#'V_VER_Ausläuferzwiebel',
+ 'V_VER_Bulbille',
'ROS_T_rosettenlose.Pflanzen','ROS_T_Halbrosettenpflanze','ROS_T_Ganzrosettenpflanzen',
'BL_AUSD_immergrün','BL_AUSD_sommergrün','BL_AUSD_vorsommergrün',
'BL_AUSD_überwinternd_grün','BL_ANAT_skleromorph','BL_ANAT_mesomorph',
- 'BL_ANAT_hygromorph','BL_ANAT_hydromorph','BL_ANAT_blattsukkulent','BL_ANAT_helomorph',
+ 'BL_ANAT_hygromorph',#'BL_ANAT_hydromorph',
+ 'BL_ANAT_blattsukkulent','BL_ANAT_helomorph',
'BL_FORM_gelappt_gefiedert_gefingert_mehrfach','BL_FORM_Vollblatt_Normalblatt',
'BL_FORM_Grasblatt_Langblatt','BL_FORM_nadelförmig','BL_FORM_röhrig',
'BL_FORM_schuppenförmig','BL_FORM_schwertförmig','REPR_T_Samen_Sporen',
@@ -56,7 +66,7 @@ traits.asym.binary <- c('LEB_F_Makrophanerophyt','LEB_F_Nanophanerophyt',
'STRAT_T_CSR','STRAT_T_R','STRAT_T_S','STRAT_T_SR')
traits0 <- traits0 %>%
- mutate_at(.vars=vars(any_of(traits.asym.binary)),
+ mutate_at(.vars=vars(all_of(traits.asym.binary)),
.funs=~(.>0)*1)
@@ -76,8 +86,8 @@ all.traits <- traits0 %>%
by="species")
traits <- all.traits %>%
filter(!is.na(LeafArea))
-dim(all.traits) #[1] 907 81
-dim(traits) #[1] 805 2
+dim(all.traits) #[1] 902 82
+dim(traits) #[1] 801 82
@@ -102,8 +112,8 @@ env.all <- env
### 3. Import species data ####
# columns in species correspond to those in env
# there is no PlotObservationID (yet)
-species0 <- read.table("_data/Mesobromion/GVRD_Mes2_veg1.csv", sep = ",", header=T)
-dim(species0) #6868 obs. of 907 variables:
+species0 <- read_csv("_data/Mesobromion/GVRD_Mes2_proz2.csv", locale = locale(encoding = 'latin1'))
+dim(species0) #6868 obs. of 903 variables:
rownames(species0) <- env0$RELEVE_NR
## select only plots already selected in env
@@ -117,7 +127,7 @@ species <- env %>%
dplyr::select(colnames(.)[which(colSums(.)!=0)])
#dplyr::select(traits$species0)
-dim(species) # [1] 5810 881
+dim(species) # [1] 5810 877
releve08trait <- species %>%
rownames_to_column("RELEVE_NR") %>%
@@ -147,10 +157,10 @@ releve08trait.samp <- sample(releve08trait, round(length(releve08trait)/10), rep
species <- species %>%
rownames_to_column("RELEVE_NR") %>%
filter(RELEVE_NR %in% releve08trait.samp) %>%
- #column_to_rownames("RELEVE_NR") %>%
#as.tbl() %>%
- dplyr::select(RELEVE_NR, one_of(traits$species0))
-
+ dplyr::select(RELEVE_NR, one_of(traits$species0)) %>%
+ column_to_rownames("RELEVE_NR") %>%
+ dplyr::select(colnames(.)[which(colSums(.)!=0)])
env <- env %>%
filter(RELEVE_NR %in% releve08trait.samp)
@@ -235,9 +245,9 @@ traits <- recode.traits(traits)
### ordered factors
-dim(species) #558 783
-dim(traits) #783 53
-dim(env) #558 8
+dim(species) #581 509
+dim(traits) #509 54
+dim(env) #581 8
@@ -282,74 +292,48 @@ env <- env %>%
species <- species %>%
rownames_to_column("RELEVE_NR")
-##export for Valerio
-write_delim(species, path="_data/Mesobromion/species.out.10perc.txt", delim="\t")
-write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
-write_delim(env, path="_data/Mesobromion/env.10perc.cov.txt", delim="\t")
+### 5. Transform and export data
+
+## presence absence
+species.pa <- species
+species.pa[,-1] <- (species[,-1]>0)*1
+## sqrt transform cov values
+species.cov <- species
+species.cov[,-1] <- sqrt(species.cov[,-1])
+#transform percentage cover to relative.cover
+species.cov <- species.cov %>%
+ mutate(sumVar = rowSums(.[-1])) %>%
+ mutate_at(.vars=vars(-RELEVE_NR),
+ .funs=~./sumVar) %>%
+ dplyr::select(-sumVar)
+dim(species.cov) #[1] 581 510
-## version without missing species
-empty <- which(colSums(species[,-1])==0)
-traits_nozero <- traits[-empty,]
-species_nozero <- species[,-(empty+1)]
+## export
+write_delim(species.pa, path="_data/Mesobromion/species.v2.10perc.pa.txt", delim="\t")
+write_delim(species.cov, path="_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
+write_delim(traits, path="_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
+write_delim(env, path="_data/Mesobromion/env.v2.10perc.cov.txt", delim="\t")
-write_delim(species_nozero , path="_data/Mesobromion/species.out.10perc_nozero.txt", delim="\t")
-write_delim(traits_nozero, path="_data/Mesobromion/traits.out.10perc_nozero.txt", delim="\t")
### create list of relevées
write_delim(species %>%
dplyr::select(RELEVE_NR),
- path="_derived/Mesobromion/ReleveList.txt", delim="\t")
+ path="_data/Mesobromion/ReleveList.txt", delim="\t")
-###### 5. Species data with cover values #####
-### version with cover values ### 4/08/2020
-env0 <- read_delim("_data/Mesobromion/GVRD_MES2_site.csv", delim = ",")
-species.proz <- read_csv("_data/Mesobromion/GVRD_Mes2_proz.csv", locale = locale(encoding = 'latin1'))
-species.proz$RELEVE_NR <- env0$RELEVE_NR
-species.proz <- species.proz %>%
- filter(RELEVE_NR %in% (species %>% pull(RELEVE_NR))) %>%
- ## delete species not appearing in any plot
- dplyr::select(colnames(.)[which(colSums(.)!=0)])
-dim(species.proz)
-write_delim(species.proz , path="_data/Mesobromion/species.out.10perc.proz.txt", delim="\t")
-#transform percentage cover to relative.cover
-species.proz <- species.proz %>%
- mutate(sumVar = rowSums(.[-1])) %>%
- mutate_at(.vars=vars(-RELEVE_NR),
- .funs=~./sumVar) %>%
- dplyr::select(-sumVar)
-dim(species.proz) #[1] 558 533
-write_delim(species.proz , path="_data/Mesobromion/species.out.10perc.cov.txt", delim="\t")
-
-## align traits to species in species.proz
-traits.proz <- recode.traits(all.traits)
-traits.proz <- data.frame(species=colnames(species.proz)[-1] ) %>%
- ### clean species names in both data.frames
- mutate(species0=as.character(species)) %>%
- rowwise() %>%
- # quick and dirty clean up names
- mutate(species=gsub(pattern="_agg_|_x_|_spec$|_agg$|_s_|_Sec_|__", replacement="_", x=species)) %>%
- mutate(species=gsub(pattern="_$", replacement = "", x = species)) %>%
- mutate(species=ifelse(is.na(word(species, 1, 2)), species, word(species, 1, 2))) %>%
- ungroup() %>%
- left_join(traits.proz %>%
- mutate(species=species0) %>%
- rowwise() %>%
- mutate(species=gsub(pattern="_agg_|_x_|_spec$|_agg$|_s_|_Sec_|__", replacement="_", x=species)) %>%
- mutate(species=gsub(pattern="_$", replacement = "", x = species)) %>%
- mutate(species=ifelse(is.na(word(species, 1, 2)), species, word(species, 1, 2))) %>%
- ungroup() %>%
- dplyr::select(-species0),
- by="species") %>%
- dplyr::select(-species) %>%
- dplyr::select(`species0`, everything())
##check for species without trait info
-traits.proz %>%
+traits %>%
filter_at(.vars=vars(-"species0"),
all_vars(is.na(.))) %>%
- dim() ## [1] 16 53 # species with no trait info
-write_delim(traits.proz, path="_data/Mesobromion/traits.out.10perc.cov.txt", delim="\t")
+ nrow() ## [1] no species with no trait info ##
+## simply because those 7 species without TRY data were excluded already
+
+traits %>%
+ filter_at(.vars=vars(-"species0"),
+ any_vars(is.na(.))) %>%
+ nrow() ## [1] 109 # species with no at least 1 NA in traits
+
#### CORRELATION BETWEEN FUZZY WEIGHTED AND BEALS MATRICES
#### WAS RUN IN THE CLUSTER WITH THE SCRIPT 01b_MesobromionCluster.R
diff --git a/01b_MesobromionCluster.R b/01b_MesobromionCluster.R
index 215c99d24923e9f43f785a9e870c17e99ecfb9d5..75528f350a98e87665fa3bf1622b2289d49ed933 100644
--- a/01b_MesobromionCluster.R
+++ b/01b_MesobromionCluster.R
@@ -56,7 +56,7 @@ get.best <- function(x, N){
Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY.bootstrap",
combinations=c("all", "sequential"), start.round=NA, relax.round=NA, max.inter.t,
- chunkn, chunk.i, nperm=199, ncores){
+ chunkn, chunk.i, nperm=199, ncores, exclude.na=F){
if(ncores>1){
require(parallel)
@@ -75,16 +75,23 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
column_to_rownames("RELEVE_NR")
traits <- read_delim(traits.path, delim="\t") %>%
as.data.frame() #%>%
-### TEMPORARY FOR TESTING!
- #dplyr::select(1:6)
-
-
+
traits <- traits %>%
- column_to_rownames("species0") %>%
rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.)) %>%
- mutate_if(~is.character(.), .funs=~as.factor(.))
-#temporary ### Use only a subset of traits
- #dplyr::select(LeafArea:Disp.unit.leng)
+ mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
+ column_to_rownames("species0")
+
+ if(exclude.na==T){
+ print("Excluding NAs and recalculating relative covers")
+ traits <- traits %>%
+ filter(complete.cases(.))
+ species <- species %>%
+ dplyr::select(rownames(traits)) %>%
+ mutate(sumVar = rowSums(.)) %>%
+ mutate_all(.funs=~./sumVar) %>%
+ dplyr::select(-sumVar)
+ }
+
if(combinations=="all") {
## create list of indices for each combination of traits up to a max number of interactions
diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R
index 77e3a7587d4aae8d7c90b090102fe1dd67c0725f..c5f7af2c46044e4afbf7e908d86d699aa10d154c 100644
--- a/02_Mesobromion_ExamineOutput.R
+++ b/02_Mesobromion_ExamineOutput.R
@@ -6,12 +6,11 @@ library(grid)
library(vegan)
source("99_HIDDEN_functions.R")
+source("01b_MesobromionCluster.R")
####1. Reimport data ################################
-## calculate corr between species composition matrix and traits
-species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
-species_nozero <- read_delim("_data/Mesobromion/species.out.10perc_nozero.txt", delim="\t")
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")
@@ -21,9 +20,24 @@ traits <- traits %>%
# filter(species0 %in% colnames(species)) %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0")
+
+traits.sign <- traits.sign %>%
+ 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))) %>%
+ rename(Trait.comb.new=Trait.comb) %>%
+ mutate(Trait.comb=1:7) %>%
+ 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
##check traits
@@ -57,66 +71,62 @@ trait.recap <- traits %>%
write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
- #dplyr::select(LeafArea:Disp.unit.leng)
- #dplyr::select(PlantHeight, LeafC.perdrymass, LeafN, StemDens, Stem.cond.dens, Seed.num.rep.unit, SLA)
-
+###### 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)
+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="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()
-###### 2. Import output from Cluster ####
-myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_round_[0-9]+.RData", full.names = T)
-#dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
-load("_derived/Mesobromion/HIDDEN_round_11.RData")
-
-
-#### 2.1 Import and adjust dataset ####
-
-#corXY = bind_rows(dataFiles) %>%
-# as.tbl()
-#rm( dataFiles)
-
-#trait.labs <- data.frame(Trait.comb=1:ncol(traits),
-# trait.name=colnames(traits)) %>% #[-which(colnames(traits) %in% c("species", "species0"))])
-# # shorten some labels
-# mutate(trait.name = fct_recode(trait.name,
-# V_VER_Wurzelknolle = "V_VER_Innovationsknopse.mit.Wurzelknolle",
-# V_VER_Speicherwurzel = "V_VER_Innovationsknospe.mit.Speicherwurzel"))
-#
-
-# ## 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.995) %>%
-# mutate(sign_minus=greater.than.perm<0.005) %>%
-# rowwise() %>%
-# mutate(ntraits= length(unlist(strsplit(Trait.comb, "_")))) %>%
-# ungroup() %>%
-# mutate_at(.vars=vars(starts_with("sign")),
-# .funs=~factor(.*1, levels=c(0,1), labels=c("FALSE", "TRUE")))
-
-corXY.ci <- corXY.ci %>%
+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:11)) %>%
- mutate_at(.vars=vars(trait1:trait11),
- .funs=~factor(., levels=trait.labs$Trait.comb, labels=trait.labs$Short_english_name)) %>%
+ 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)) %>%
arrange(ntraits, Coef.obs) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything())
-#### 2.2 Graph of r(XY) when taking traits one by one ####
+
+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)
+
+
+#### 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(ntraits==1) %>%
@@ -144,52 +154,53 @@ mydata.byone <- corXY.ci %>%
)
-ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI.png", dpi=400,
+ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI_pa.png", dpi=400,
width=4, height=6, plot = top.first)
-#### 2.3 (Huge) Plot with all trait combination ####
+#### 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:trait11), all_vars(. %in% traits.sign.alone | is.na(.)))%>%
+ filter_at(.vars=vars(trait1:trait7),
+ all_vars(. %in% traits.sign.alone | is.na(.)))%>%
mutate(seq=1:n())
-allpredictors <- top.first %+% (mydata) +
- scale_y_continuous(breaks=mydata$seq,
- labels=NULL) +
- # labels=mydata$Trait.comb) +
-
- theme_classic() +
- theme(axis.text.y = element_text(size=3))# +
- # facet_wrap(.~ntraits, scales = "free_y", nrow=3)
-
-ggsave(filename = "_pics/All_predictors_sign_individually_CI_faceted_0506.png", dpi=400,
- width=10, height=10, allpredictors)
-
-#### 2.4 Graph of r(XY) using best combination of traits at each level of interaction N
+# allpredictors <- top.first %+% (mydata) +
+# scale_y_continuous(breaks=mydata$seq,
+# labels=NULL) +
+# # labels=mydata$Trait.comb) +
+#
+# theme_classic() +
+# theme(axis.text.y = element_text(size=3))# +
+# # facet_wrap(.~ntraits, scales = "free_y", nrow=3)
+#
+# ggsave(filename = "_pics/All_predictors_sign_individually_CI_faceted_0506.png", dpi=400,
+# width=10, height=10, allpredictors)
+
+######## 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){
+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 <- trait.labs %>%
+ 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)
+top.one.by.one <- get.best(mydata, N=1, labs=trait.labs.sign)
## Routine to extract the best combination at each level of interaction (up to max traits)
-maxtraits <- 11
+maxtraits <- 7
for(nn in 1:maxtraits){
if(nn==1) {
- best.at.1 <- get.best(mydata, N=nn)
+ best.at.1 <- get.best(mydata, N=nn, labs=trait.labs.sign)
newdata <- mydata %>%
- filter_at(.vars=vars(trait1:trait11),
+ filter_at(.vars=vars(trait1:trait7),
.vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.)))
new.best.row <- newdata %>%
filter(Trait.comb==best.at.1$Trait.comb)
@@ -208,7 +219,7 @@ for(nn in 1:maxtraits){
pull(Trait.comb)
if(length(better$Trait.comb>0)){
- better$trait.name <- trait.labs %>%
+ better$trait.name <- trait.labs.sign %>%
filter(Trait.comb %in% strsplit(better$Trait.comb, split = "_")[[1]]) %>%
pull(trait.name)
@@ -231,15 +242,15 @@ for(nn in 1:maxtraits){
}
}
-best.5traits <- corXY.ci %>%
+best.2traits <- corXY.ci %>%
filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
- dplyr::select(trait1:trait11) %>%
+ dplyr::select(trait1:trait7) %>%
mutate_all(~as.character(.))
-best.5traits <- as.character(best.5traits[1,1:5])
+best.2traits <- as.character(best.2traits[1,1:2])
### Create dataset with best combinations + all the one-way combinations
mydata.best <- mydata %>%
- filter_at(.vars=vars(trait1:trait11),
+ filter_at(.vars=vars(trait1:trait7),
# .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)) %>%
@@ -259,7 +270,7 @@ mydata.best <- mydata %>%
write_csv(mydata.best %>%
dplyr::select(Trait.comb:sign_plus),
- path = "_derived/Mesobromion/BestSolutionTiers.csv")
+ path = "_derived/Mesobromion/S5_BestSolutionTiers.csv")
### Graph of all the best combinations with text legend
(top.all <- ggplot(data=(mydata.best %>%
@@ -284,7 +295,7 @@ tt2 <- ttheme_minimal(
rowhead = list(fg_params=list(cex = .7)))
(topall.leg <- cowplot::plot_grid(top.all,
- tableGrob(trait.labs %>%
+ tableGrob(trait.labs.sign %>%
mutate(Code=1:n()) %>%
dplyr::select(Code, Trait=Short_english_name) %>%
filter(Trait %in% traits.sign.alone),
@@ -294,27 +305,85 @@ 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.png", dpi=400,
+ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_pa.png", dpi=400,
width=6, height=3, topall.leg)
break()
+##### 2.2 Cover values ######
+
+#### ____TO DO ########
+### all traits taken alone
+myfilelist0 <- list.files(path="_derived/Mesobromion/Cover-pa/", pattern="HIDDENcov.nona_[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/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:3)) %>%
+ 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)
+
+### merge together
+corXY.ci <- corXY.comb.ci %>%
+ bind_rows(corXY.alone.ci %>%
+ filter(!trait1 %in% trait.labs.sign$Short_english_name)) %>%
+ arrange(Coef.obs)
+
+
+
##### 3. Calculate CWM ####
-#species <- species %>%
-# dplyr::select(colnames(.)[which(colSums(.)!=0)])
-
+#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
+
+
-CWM.wide <- species %>%
- #rownames_to_column("RELEVE_NR") %>%
- # reshape2::melt(.id="RELEVE_NR") %>%
- pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "pres") %>%
- #rename(species0=variable, pres=value) %>%
+CWM.wide <- species.cov %>%
+ pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "rel.cov") %>%
as.tbl() %>%
- filter(pres>0) %>%
+ filter(rel.cov>0) %>%
arrange(RELEVE_NR) %>%
## attach traits
left_join(traits %>%
@@ -323,13 +392,20 @@ CWM.wide <- species %>%
group_by(RELEVE_NR) %>%
select_if(.predicate=~is.numeric(.)) %>% # CWM calculated ONLY for numeric traits
summarize_at(.vars = vars(any_of(traits.sign.alone)), #GF_Macrophan:Rosette),
- .funs = list(~weighted.mean(., pres, na.rm=T))) %>%
+ .funs = list(~weighted.mean(., rel.cov, na.rm=T))) %>%
dplyr::select(RELEVE_NR, order(colnames(.))) %>%
column_to_rownames("RELEVE_NR")
#### 4. PCA Graphs ####
library(vegan)
library(ggrepel)
+
+##palette
+orange = "#d95f02"
+oilgreen = "#1b9e77"
+myblue = "#7570b3"
+
+
##from https://stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2
circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
r = diameter / 2
@@ -338,10 +414,121 @@ circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
yy <- center[2] + r * sin(tt)
return(data.frame(x = xx, y = yy))
}
+dat <- circleFun(diameter = 2, npoints = 100)
+### 4.0 PCA of X (Fuzzy-weighted) matrix + CWMs ####
+W.fuzzy <- matrix.x(comm=species.cov %>%
+ column_to_rownames("RELEVE_NR"),
+ traits=traits.cov %>%
+ 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!
+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
+env.cor <- cor(env %>%
+ dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC),
+ scores.pca[,1:4]) #double check
+
+myvectors <- as.data.frame(env.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(category="Env") %>%
+ bind_rows(as.data.frame(cwm.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(category="Trait")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.2traits, "bold", "italic")) %>%
+ mutate(category=as.factor(category)) %>%
+ mutate(mycol=ifelse(category=="Trait", oilgreen, orange))
+
+basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
+ theme_bw() +
+ scale_color_identity() +
+ scale_y_continuous(limits=c(-1, 1)) +
+ scale_x_continuous(limits=c(-1, 1)) + coord_equal() +
+ theme(panel.grid = element_blank())
+
+PCA.fuzz1_2 <- basemap0 +
+ geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:4]*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),
+ 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=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
+ 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) +
+ geom_segment(data=myvectors,
+ aes(x=0, xend=PC3, y=0, yend=PC4, 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,
+ 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=""))
+
+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())
+
+#### 4.0b Alternative showing species scores ####
+tmp <- as.data.frame(pca.fuzz$CA$v[,1:4]*7) %>%
+ mutate(species0=colnames(species.cov)[-1]) %>%
+ mutate(species=species0) %>%
+ separate(species0, sep="_", into=c("Gen", "Spe")) %>%
+ mutate(Gen=substr(Gen, 1, 3)) %>%
+ mutate(Spe=substr(Spe, 1, 3)) %>%
+ mutate(labels=paste(Gen, Spe, sep="_"))
+
+
+PCAfuzz1_2.sp <- basemap0 %+% tmp +
+ geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
+ geom_segment(data=myvectors,
+ aes(x=0, xend=PC1, y=0, yend=PC2, 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=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
+ 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) +
+ geom_segment(data=myvectors,
+ aes(x=0, xend=PC3, y=0, yend=PC4, 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,
+ 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=""))
+
+
+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)
+
+
+
+
+
+
#### 4.1 PCA of Y (Bealls) matrix + CWM ####
-W.beals <- as.data.frame(beals(species_nozero %>%
- column_to_rownames("RELEVE_NR"), include=T, type=2))
+W.beals <- as.data.frame(beals(species.cov %>%
+ column_to_rownames("RELEVE_NR") %>%
+ mutate_all(~(.>0)*1), #transform to p\a
+ include=T, type=2))
write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv")
pca.out <- rda(W.beals)
varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
@@ -364,24 +551,30 @@ 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])
sink("_derived/Mesobromion/EnvFit_CWMs_env.txt")
-cwms.envfit$vectors
-env.envfit$vectors
+cwms.cor
+env.cor
+fuzz.cor
sink()
-dat <- circleFun(diameter = 2, npoints = 100)
-
myvectors <- as.data.frame(env.cor) %>%
rownames_to_column("mylab") %>%
mutate(category="Env") %>%
bind_rows(as.data.frame(cwms.cor) %>%
rownames_to_column("mylab") %>%
mutate(category="Trait")) %>%
- mutate(fontface0=ifelse(mylab %in% best.5traits, "bold", "plain")) %>%
+ 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(category=as.factor(category)) %>%
- mutate(mycol=ifelse(category=="Trait", "Dark blue", "Dark red"))
+ mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
+ mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol))
+
+
basemap0 <- ggplot(data=as.data.frame(pca.out$CA$u[,1:4]*5)) +
theme_bw() +
@@ -398,10 +591,9 @@ PCA1_2 <- basemap0 +
aes(x=0, xend=PC1, y=0, yend=PC2, 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(data=myvectors %>%
- mutate_if(~is.numeric(.), ~(.)*1.2),
+ geom_label_repel(data=myvectors,
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
- position = position_dodge(1)) +#, segment.alpha=0.5, segment.colour=gray(0.8)) +
+ position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
@@ -434,20 +626,19 @@ tmp <- as.data.frame(pca.out$CA$v[,1:4]*5) %>%
mutate(labels=paste(Gen, Spe, sep="_"))
-PCA1_2.sp <- baseplot %+% tmp +
+PCA1_2.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
geom_segment(data=myvectors,
aes(x=0, xend=PC1, y=0, yend=PC2, 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(data=myvectors %>%
- mutate_if(~is.numeric(.), ~(.)*1.2),
+ geom_label_repel(data=myvectors,
aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2,
- position = position_dodge(1)) +#, segment.alpha=0.5, segment.colour=gray(0.8)) +
+ position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
-PCA3_4.sp <- baseplot %+% tmp +
+PCA3_4.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC3, y=PC4, label=labels), size=2, alpha=0.8) +
geom_segment(data=myvectors,
aes(x=0, xend=PC3, y=0, yend=PC4, col=mycol),
@@ -468,36 +659,36 @@ ggsave("_pics/Fig6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PC
-#### 4.2 PCA of CWMs #####
-CWM.pca <- vegan::rda(CWM.wide, scale=T)
-varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
-
-ggplot() +
- geom_point(data=as.data.frame(CWM.pca$CA$u[,1:2]),
- aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(CWM.pca$CA$v*2) %>%
- rownames_to_column("Trait") %>%
- mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
- aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_text(data=as.data.frame(CWM.pca$CA$v*2.1) %>%
- rownames_to_column("Trait") %>%
- mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
- aes(x=PC1, y=PC2, label=Trait, col=in.try), size=3 ) +
- scale_color_identity() +
- xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
- theme_bw() +
- theme(panel.grid = element_blank()) +
- title("PCA of CWMs")
+# #### 4.2 PCA of CWMs #####
+# CWM.pca <- vegan::rda(CWM.wide, scale=T)
+# varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
+#
+# ggplot() +
+# geom_point(data=as.data.frame(CWM.pca$CA$u[,1:2]),
+# aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) +
+# geom_segment(data=as.data.frame(CWM.pca$CA$v*2) %>%
+# rownames_to_column("Trait") %>%
+# mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
+# aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+# geom_text(data=as.data.frame(CWM.pca$CA$v*2.1) %>%
+# rownames_to_column("Trait") %>%
+# mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")),
+# aes(x=PC1, y=PC2, label=Trait, col=in.try), size=3 ) +
+# scale_color_identity() +
+# xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+# ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+# theme_bw() +
+# theme(panel.grid = element_blank()) +
+# title("PCA of CWMs")
#### 4.3 PCA of individual traits #####
# 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 %>%
+data2 <- traits.cov %>%
dplyr::select(any_of(traits.sign.alone)) %>%
rownames_to_column("species") %>%
- filter(species %in% colnames(species_nozero)) %>%
+ filter(species %in% colnames(species.cov)) %>%
column_to_rownames("species")
corr1 <- cor(data2, use="pairwise.complete.obs")
@@ -525,34 +716,34 @@ baseplot <- ggplot(data=as.data.frame(pca.scores*.2)) +
PCA.t1 <- baseplot +
geom_point(aes(x=pca1, y=pca2), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
- aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col="Dark green"),
+ aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col=oilgreen),
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=as.data.frame(pca2$loadings[,1:5]) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
- aes(x=Comp.1, y=Comp.2, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ mutate(fontface0=ifelse(Trait %in% best.2traits, "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")
PCA.t2 <- baseplot +
geom_point(aes(x=pca3, y=pca4), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
- aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col="Dark green"),
+ aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col=oilgreen),
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=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.5traits, "bold", "plain")),
- aes(x=Comp.3, y=Comp.4, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ mutate(fontface0=ifelse(Trait %in% best.2traits, "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="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
PC_traits <- cowplot::plot_grid(PCA.t1, PCA.t2, nrow=1)
-ggsave("_pics/FigS6_PCA_Traits_1-4_only11.png", width=10, height=5, dpi=300, PC_traits)
+ggsave("_pics/S6_PCA_Traits_1-4_only7.png", width=10, height=5, dpi=300, PC_traits)
##### 4.3b Alternative version of figS6, showing the species ####
tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
@@ -567,32 +758,32 @@ tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
PCA.t1.sp <- baseplot %+% tmp +
geom_text(aes(x=pca1, y=pca2, label=labels), size=2, alpha=0.7) +
geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
- aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col="Dark green"),
+ aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col=oilgreen),
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=as.data.frame(pca2$loadings[,1:5]) %>%
rownames_to_column("Trait") %>%
- mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
- aes(x=Comp.1, y=Comp.2, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ mutate(fontface0=ifelse(Trait %in% best.2traits, "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))
PCA.t2.sp <- baseplot %+% tmp +
geom_text(aes(x=pca3, y=pca4, label=labels), size=2, alpha=0.7) +
geom_segment(data=as.data.frame(pca2$loadings[,1:5]),
- aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col="Dark green"),
+ aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col=oilgreen),
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=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.5traits, "bold", "plain")),
- aes(x=Comp.3, y=Comp.4, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ mutate(fontface0=ifelse(Trait %in% best.2traits, "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="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
-ggsave("_pics/FigS6_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA.t1.sp)
-ggsave("_pics/FigS6_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.t2.sp)
+ggsave("_pics/S6a_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA.t1.sp)
+ggsave("_pics/S6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.t2.sp)
## Create list of species labels
write_delim(tmp %>%
@@ -603,43 +794,6 @@ write_delim(tmp %>%
-
-#### 4.4 PCoA of X (Fuzzy weighted) matrix ####
-library(vegan)
-fuzzy <- read.table("_data/Mesobromion/X_Dry-Grassland_581plots_488spp_R.txt", header=T)#, delim="\t")
-cwms <- read.table("_data/Mesobromion/X_scores_T_Dry-Grassland_581plots_5traits_Rao_R.txt", header=T) %>%
- dplyr::select(-Axis1, -Axis2) #, delim="\t")
-
-pcoa.out <- wcmdscale(dist(fuzzy), eig=T)
-varexpl <- round((pcoa.out$eig)/sum(pcoa.out$eig)*100,1)
-cwms.envfit <- envfit(pcoa.out, cwms)
-
-
-plot(pcoa.out, type = "p")
-plot(cwms.envfit)
-
-ggplot() +
- geom_point(data=as.data.frame(pcoa.out$points[,1:2]),
- aes(x=Dim1, y=Dim2), pch="+", size=2, alpha=0.8) +
- geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.0008),
- aes(x=0, xend=Dim1, y=0, yend=Dim2), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
- geom_text(data=as.data.frame(cwms.envfit$vectors$arrows*.00092) %>%
- rownames_to_column("Trait"),
- aes(x=Dim1, y=Dim2, label=Trait), col="Dark blue", size=3 ) +
- xlab(paste("PCoA1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PCoA2 (", varexpl[2], "%)", sep="")) +
- theme_bw() +
- theme(panel.grid = element_blank())
-
-ggsave("_data/Mesobromion/PCoA_fuzzy.png", width=6, height=4, dpi=300, last_plot())
-
-
-
-
-
-
-
-
#### 5 Map of Plot distribution #####
library(rgdal)
library(sp)
@@ -692,15 +846,15 @@ ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basema
#### 6 Correlation matrices ####
#### 6.1 Trait level #####
library(corrplot)
-traits11 <- traits %>%
+traits7 <- traits.cov %>%
rownames_to_column("species") %>%
- filter(species %in% colnames(species_nozero)) %>%
+ filter(species %in% colnames(species.cov)) %>%
column_to_rownames("species") %>%
dplyr::select(any_of(traits.sign.alone)) %>%
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(all_of(best.5traits), everything())
+ relocate(all_of(best.2traits), everything())
-res <- cor(traits11, use = "pairwise.complete.obs")
+res <- cor(traits7, use = "pairwise.complete.obs")
png(file="_pics/S7_Correlations_Trait.png", width=8, height=6.5, units = "in", res=300)
corrplot(res, type = "upper",
tl.col = "black", tl.srt = 45, number.cex=0.6, addCoef.col = "black", diag=F)
@@ -708,9 +862,9 @@ dev.off()
#### 6.2 CWM level ####
-res2 <- cor(CWM.wide0 %>%
+res2 <- cor(CWM.wide %>%
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(all_of(best.5traits), everything()))
+ relocate(all_of(best.2traits), everything()))
png(file="_pics/S8_Correlations_CWMs.png", width=8, height=6.5, units = "in", res=300)
corrplot(res2, type = "upper",
tl.col = "black", tl.srt = 45, number.cex=0.6, addCoef.col = "black", diag=F)
@@ -720,3 +874,6 @@ dev.off()
+
+
+
diff --git a/cli_01b.r b/cli_01b.r
index b7955685bc858ac7d09a047c68610479812086ab..cd615563217207f33a53f06aa1fed3fced60b6b8 100644
--- a/cli_01b.r
+++ b/cli_01b.r
@@ -71,6 +71,14 @@ options <- list (
default = NA,
help = "After what round should the selection of new best combos stop being based on c.i.?",
metavar = "4"
+ ),
+ make_option(
+ opt_str = c("-e", "--exclude.na"),
+ dest = "exclude.na",
+ type = "logical",
+ default = F,
+ help = "should species with ANY nas in traits be excluded? (and relative covers recalculated?)",
+ metavar = "4"
)
)
@@ -99,10 +107,11 @@ chunkn <- cli$options$chunkn
chunk.i <- cli$options$chunk.i
nperm <- cli$options$nperm
ncores <- cli$options$ncores
+exclude.na <- cli$options$exclude.na
# ------------------------------------------------------------------------------
# actual program
# ------------------------------------------------------------------------------
source("01b_MesobromionCluster.R")
-Mesobromion(species.path, traits.path, output, myfunction, combinations, start.round, relax.round, max.inter.t, chunkn, chunk.i, nperm, ncores)
+Mesobromion(species.path, traits.path, output, myfunction, combinations, start.round, relax.round, max.inter.t, chunkn, chunk.i, nperm, ncores, exclude.na)
diff --git a/session.R b/session.R
index 78b6f932caef24f863d324c3fbabd74494c9e1c3..c5e4d2517378f965ac0e5a3372e3a7844e719372 100644
--- a/session.R
+++ b/session.R
@@ -1,18 +1,19 @@
-species.path <- "_data/Mesobromion/species.out.10perc.cov.txt"
-traits.path <- "_data/Mesobromion/traits.out.10perc.cov.txt"
+species.path <- "_data/Mesobromion/species.v2.10perc.cov.txt"
+traits.path <- "_data/Mesobromion/traits.v2.10perc.txt"
output <- "_derived/Mesobromion/Cover/HIDDENtest"
myfunction <- "get.corXY.bootstrap"
max.inter.t <- 1
chunk.i <- 1
-nperm <- 49
+nperm <- 3
ncores <- 1
chunkn <- 1
combinations <- "all"
start.round <- 1
relax.round <- 1
+exclude.na <- T
source("01b_MesobromionCluster.R")
Mesobromion(species.path, traits.path, output, myfunction,
combinations, start.round, relax.round, max.inter.t,
- chunkn, chunk.i, nperm, ncores)
+ chunkn, chunk.i, nperm, ncores, exclude.na)
diff --git a/submit_01b.sh b/submit_01b.sh
index a62ea9ac2eaca5123d7432838ce2fae0875f6ee0..af5f79f62a8d13c5b911f8e4ecf43ff1db5c6f6c 100644
--- a/submit_01b.sh
+++ b/submit_01b.sh
@@ -28,6 +28,7 @@ Rscript \
--ncores ${NSLOTS:-1} \
--start.round 1 \
--relax.round 4 \
+ --exclude.na T \
/data/splot/HIDDEN/species.out.10perc.cov.txt \
/data/splot/HIDDEN/traits.out.10perc.cov.txt \
"$output" \
diff --git a/submit_01b_array.sh b/submit_01b_array.sh
index 175f5f186ce6113257c537a8094d633578428800..6565bce181840bf5b34cc0a7557a939e9496a988 100644
--- a/submit_01b_array.sh
+++ b/submit_01b_array.sh
@@ -30,9 +30,11 @@ Rscript \
--ncores ${NSLOTS:-1} \
--start.round 1 \
--relax.round 1 \
+ --exclude.na T \
/data/splot/HIDDEN/species.out.10perc.cov.txt \
/data/splot/HIDDEN/traits.out.10perc.cov.txt \
"$output" \
"get.corXY.bootstrap" \
"all" \
+