From 0bc2be36b6f3d5288209439466241f1e1784aab3 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Wed, 26 Aug 2020 16:38:38 +0200
Subject: [PATCH] Aligned to new results
---
00_Mesobromion_DataPreparation.R | 85 +++---
02_Mesobromion_ExamineOutput.R | 445 ++++++++++++++++++++-----------
2 files changed, 325 insertions(+), 205 deletions(-)
diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index 301734e..658f902 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -412,19 +412,22 @@ traits.sign <- traits %>%
dplyr::select(species0, any_of(traits.sign.alone))
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
- <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 36 0.316 0.154 0.290 0.370 0.997 999 TRUE FALSE 1 PlantHeight
- 2 50 0.286 0.154 0.257 0.339 0.977 999 TRUE FALSE 1 BL_ANAT
- 3 30 0.259 0.114 0.232 0.308 0.974 999 TRUE FALSE 1 BL_DAU
- 4 2 0.256 0.228 0.222 0.303 0.992 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
- 5 20 0.253 0.0703 0.205 0.324 0.992 999 TRUE FALSE 1 V_VER_Fragmentation
- 6 49 0.252 0.116 0.226 0.303 0.967 999 TRUE FALSE 1 STRAT_T
- 7 32 0.251 0.127 0.226 0.303 0.968 999 TRUE FALSE 1 SLA
- 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
+### COV - NONAs - NEW DATASET
+"# A tibble: 49 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 46 0.317 0.106 0.289 0.365 0.997 999 TRUE FALSE 1 Leaf_Scleroph
+ 2 33 0.267 0.168 0.244 0.317 0.990 999 TRUE FALSE 1 Height
+ 3 29 0.251 0.206 0.231 0.300 0.976 999 TRUE FALSE 1 SLA
+ 4 2 0.250 0.111 0.218 0.296 0.988 999 TRUE FALSE 1 GF_Nanophan
+ 5 27 0.237 0.167 0.213 0.284 0.975 999 TRUE FALSE 1 FP_Dur
+ 6 30 0.212 0.125 0.192 0.260 0.953 999 TRUE FALSE 1 LeafCdrymass
+ 7 5 0.210 0.0489 0.176 0.267 0.977 999 TRUE FALSE 1 GF_Hemiphan
+ 8 17 0.207 0.123 0.166 0.273 0.965 999 TRUE FALSE 1 VP_Fragm
+ 9 32 0.214 0.0827 0.195 0.262 0.939 999 FALSE FALSE 1 LeafP
+10 31 0.195 0.133 0.176 0.241 0.919 999 FALSE FALSE 1 LeafN"
+
+### COV - NONAs - old dataset - First selection of plots
"# A tibble: 48 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>
@@ -437,13 +440,6 @@ write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.cov.sign.txt",
7 30 0.241 0.207 0.213 0.286 0.951 999 TRUE FALSE 1 LeafP
8 44 0.238 0.187 0.213 0.285 0.971 999 TRUE FALSE 1 STRAT_T "
-### COV - without deleting NAs
-"# 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.117 0.280 0.356 0.995 999 TRUE FALSE 1 PlantHeight
- 2 32 0.247 0.127 0.221 0.299 0.976 999 TRUE FALSE 1 SLA
- 3 35 0.241 0.175 0.218 0.288 0.951 999 TRUE FALSE 1 LeafP "
#### 2. Traits individually significant for Presence|absence data####
traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
@@ -470,26 +466,24 @@ traits.sign <- traits %>%
dplyr::select(species0, any_of(traits.sign.alone))
write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.pa.sign.txt", delim="\t")
-## Pres Abs No Nas
-"# A tibble: 13 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.313 0.0582 0.284 0.362 0.996 999 TRUE FALSE 1 PlantHeight
- 2 50 0.309 0.0598 0.279 0.361 0.992 999 TRUE FALSE 1 BL_ANAT
- 3 2 0.284 0.0746 0.248 0.331 0.994 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
- 4 32 0.260 0.140 0.234 0.305 0.981 999 TRUE FALSE 1 SLA
- 5 30 0.256 0.215 0.223 0.311 0.980 999 TRUE FALSE 1 BL_DAU
- 6 35 0.254 0.0577 0.226 0.308 0.971 999 TRUE FALSE 1 LeafP
- 7 49 0.253 0.209 0.222 0.303 0.983 999 TRUE FALSE 1 STRAT_T
- 8 33 0.251 0.272 0.223 0.307 0.979 999 TRUE FALSE 1 LeafC.perdrymass
- 9 34 0.234 0.0945 0.209 0.283 0.968 999 TRUE FALSE 1 LeafN
-10 5 0.233 0.115 0.196 0.291 0.981 999 TRUE FALSE 1 LEB_F_Hemiphanerophyt
-11 38 0.222 0.152 0.198 0.274 0.966 999 TRUE FALSE 1 Seed.length
-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
+## NONAs - New data
+"# A tibble: 49 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 46 0.326 0.190 0.297 0.376 0.999 999 TRUE FALSE 1 Leaf_Scleroph
+ 2 33 0.285 0.0765 0.256 0.334 0.997 999 TRUE FALSE 1 Height
+ 3 2 0.281 0.0683 0.250 0.329 0.998 999 TRUE FALSE 1 GF_Nanophan
+ 4 29 0.265 0.0748 0.243 0.311 0.990 999 TRUE FALSE 1 SLA
+ 5 30 0.254 0.105 0.227 0.306 0.988 999 TRUE FALSE 1 LeafCdrymass
+ 6 27 0.245 0.119 0.220 0.296 0.982 999 TRUE FALSE 1 FP_Dur
+ 7 5 0.245 0.131 0.213 0.303 0.986 999 TRUE FALSE 1 GF_Hemiphan
+ 8 31 0.235 0.0744 0.214 0.285 0.974 999 TRUE FALSE 1 LeafN
+ 9 32 0.231 0.189 0.204 0.283 0.976 999 TRUE FALSE 1 LeafP
+10 35 0.207 0.0881 0.185 0.255 0.951 999 TRUE FALSE 1 SeedLength
+11 1 0.191 0.0483 0.161 0.239 0.961 999 TRUE FALSE 1 GF_Macrophan
+12 34 0.197 0.0669 0.179 0.244 0.948 999 FALSE FALSE 1 SeedMass "
+
+## NONAs - Old data
"# A tibble: 48 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>
@@ -508,14 +502,3 @@ write_delim(traits.sign, path="_data/Mesobromion/traits.v2.10perc.pa.sign.txt",
13 40 0.209 0.128 0.185 0.259 0.957 999 TRUE FALSE 1 Disp.unit.leng "
-## Pres Abs - Without excluding species with NA in traits
-"# A tibble: 7 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.307 0.109 0.279 0.358 0.999 999 TRUE FALSE 1 PlantHeight
-2 32 0.260 0.122 0.235 0.308 0.987 999 TRUE FALSE 1 SLA
-3 35 0.253 0.135 0.228 0.299 0.986 999 TRUE FALSE 1 LeafP
-4 33 0.240 0.184 0.210 0.296 0.975 999 TRUE FALSE 1 LeafC.perdrymass
-5 34 0.233 0.166 0.207 0.283 0.972 999 TRUE FALSE 1 LeafN
-6 38 0.228 0.0792 0.205 0.278 0.975 999 TRUE FALSE 1 Seed.length
-7 39 0.226 0.141 0.201 0.282 0.973 999 TRUE FALSE 1 LDMC "
diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R
index da06627..da518f6 100644
--- a/02_Mesobromion_ExamineOutput.R
+++ b/02_Mesobromion_ExamineOutput.R
@@ -25,18 +25,21 @@ get.best <- function(x, N, labs){
####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")
-traits.sign.pa <- read_delim("_data/Mesobromion/traits.v2.10perc.pa.sign.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") %>%
+ left_join(read_delim("_data/Mesobromion/TraitLabels_Long.csv", delim=",") %>%
+ dplyr::select(-trait.name),
+ by=c("trait.name"="Short_english_name")) %>%
rownames_to_column("Trait.comb") %>%
- mutate_at(.vars=vars(Short_english_name:Long_English_name), ~ifelse(is.na(.), trait.name, {.}))
+ mutate_at(.vars=vars(Long_English_name), ~ifelse(is.na(.), trait.name, {.}))
env <- read_delim("_data/Mesobromion/env.v2.10perc.cov.txt", delim="\t")
species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
traits <- traits %>%
as.data.frame() %>%
+ mutate_at(.vars=vars(Leaf_Scleroph, LifeSpan, Rosette),
+ .funs=~as.ordered(.)) %>%
# filter(species0 %in% colnames(species)) %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
column_to_rownames("species0")
@@ -62,9 +65,6 @@ traits <- traits %>%
# 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
##check traits
@@ -72,8 +72,6 @@ is.binary <- function(x){(all(na.omit(x) %in% 0:1))}
####1.1 Table S3 - Trait recap #####
trait.recap <- traits %>%
- mutate_at(.vars=vars(Leaf_Scleroph, LifeSpan, Rosette),
- .funs=~as.ordered(.)) %>%
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",
@@ -95,8 +93,8 @@ trait.recap <- traits %>%
mutate(Variable=factor(Variable, levels=colnames(traits))) %>%
arrange(Variable) %>%
left_join(trait.labs %>%
- dplyr::select(Short_english_name, Long_English_name),
- by=c("Variable" = "Short_english_name")) %>%
+ dplyr::select(trait.name, Long_English_name),
+ by=c("Variable" = "trait.name")) %>%
dplyr::select(Code=Variable, Variable=Long_English_name, everything())
## calculate trait coverage across plots, both on p\a and cover
@@ -163,11 +161,11 @@ corXY.all = bind_rows(dataFiles1) %>%
corXY.all.ci <- get.ci(corXY.all)
corXY.all.ci <- corXY.all.ci %>%
mutate(Trait.comb2=Trait.comb) %>%
- separate(Trait.comb2, into=paste0("trait", 1:8)) %>%
- mutate_at(.vars=vars(trait1:trait8),
+ 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)) %>%
+ labels=trait.labs$trait.name)) %>%
arrange(ntraits, desc(Coef.obs)) %>%
#filter(ntraits>1) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
@@ -182,12 +180,13 @@ corXY.ci <- corXY.all.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 %>%
+traits.sign.alone.cov <- 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 "
+### OLD ### #"Height Leaf_Anat GF_Nanophan BL_Dur VP_Fragm SLA LeafP CSR "
+### NEW ### #Leaf_Scleroph Height SLA GF_Nanophan FP_Dur GF_Hemiphan VP_Fragm
### Plot combinations one by one
mydata.byone <- corXY.ci %>%
@@ -217,8 +216,8 @@ ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI_cov.png", dpi=400,
mydata <- corXY.ci %>%
filter(run=="seq") %>%
- filter_at(.vars=vars(trait1:trait8),
- all_vars(. %in% traits.sign.alone | is.na(.)))%>%
+ filter_at(.vars=vars(trait1:trait7),
+ all_vars(. %in% traits.sign.alone.cov | is.na(.)))%>%
mutate(seq=1:n())
#### ~2.1.3 All traits - (Huge) Graph of r(XY) (w/ combinations) ####
@@ -241,12 +240,12 @@ mydata <- corXY.ci %>%
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 <- 8
+maxtraits <- 7
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),
+ 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)
@@ -288,17 +287,20 @@ for(nn in 1:maxtraits){
}
}
-best.4traits <- corXY.ci %>%
+best.traits.cov <- 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])
+ dplyr::select(trait1:trait7) %>%
+ mutate_all(~as.character(.)) %>%
+ dplyr::select(colnames(.)[which(colSums(is.na(.))==0)])
+
+best.traits.cov <- as.character(best.traits.cov[1,])
+#"Leaf_Scleroph" "FP_Dur" "VP_Fragm" "Height" "SLA"
### Create dataset with best combinations + all the one-way combinations
mydata.best <- mydata %>%
- filter_at(.vars=vars(trait1:trait8),
+ 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(.))) %>%
+ .vars_predicate = all_vars(. %in% traits.sign.alone.cov | is.na(.))) %>%
#filter(ntraits <= length(best$trait.name)) %>%
filter(ntraits>1) %>%
filter(sign_plus==T) %>%
@@ -309,7 +311,7 @@ mydata.best <- mydata %>%
bind_rows(corXY.ci %>%
filter(run=="seq") %>%
filter(ntraits==1) %>%
- filter(trait1 %in% traits.sign.alone)) %>%
+ filter(trait1 %in% traits.sign.alone.cov)) %>%
arrange(ntraits, Coef.obs) %>%
mutate(seq=1:n()) %>%
mutate(sign_plus=factor(Trait.comb %in% best.progr))
@@ -343,8 +345,8 @@ tt2 <- ttheme_minimal(
(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),
+ dplyr::select(Code, Trait=trait.name) %>%
+ filter(Trait %in% traits.sign.alone.cov),
theme=tt2, rows = NULL),
#tableGrob(trait.labs %>%
# mutate(Code=1:n()) %>%
@@ -392,11 +394,11 @@ corXY.all = bind_rows(dataFiles) %>%
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),
+ 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)) %>%
+ labels=trait.labs$trait.name)) %>%
arrange(ntraits, desc(Coef.obs)) %>%
dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>%
mutate(run="seq")
@@ -410,9 +412,12 @@ 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
+### OLD ###[1] Leaf_Scleroph Height GF_Nanophan BL_Dur SLA LeafP LeafCdrymass CSR
#[9] LeafN SeedLength LDMC GF_Hemiphan
+# [1] Leaf_Scleroph Height GF_Nanophan SLA LeafCdrymass FP_Dur GF_Hemiphan LeafN
+# LeafP SeedLength GF_Macrophan
+
### Plot combinations one by one
mydata.byone <- corXY.ci.pa %>%
filter(ntraits==1)%>%
@@ -440,7 +445,7 @@ ggsave(filename = "_pics/SXXX_TopFirstPredictors_CI_pa.png", dpi=400,
######## 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),
+ filter_at(.vars=vars(trait1:trait11),
all_vars(. %in% traits.sign.alone.pa | is.na(.)))%>%
mutate(seq=1:n())
@@ -448,12 +453,12 @@ mydata <- corXY.ci.pa %>%
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
+maxtraits <- 11
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),
+ filter_at(.vars=vars(trait1:trait11),
.vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.)))
new.best.row <- newdata %>%
filter(Trait.comb==best.at.1$Trait.comb)
@@ -495,15 +500,17 @@ for(nn in 1:maxtraits){
}
}
-best.4traits.pa <- corXY.ci.pa %>%
+best.traits.pa <- corXY.ci.pa %>%
filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
dplyr::select(trait1:trait4) %>%
- mutate_all(~as.character(.))
-best.4traits.pa <- as.character(best.4traits.pa[1,1:4])
+ mutate_all(~as.character(.))%>%
+ dplyr::select(colnames(.)[which(colSums(is.na(.))==0)])
+best.traits.pa <- as.character(best.traits.pa[1,])
+#[1] "Leaf_Scleroph" "FP_Dur" "LeafCdrymass" "Height"
### Create dataset with best combinations + all the one-way combinations
mydata.best <- mydata %>%
- filter_at(.vars=vars(trait1:trait14),
+ filter_at(.vars=vars(trait1:trait11),
# .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)) %>%
@@ -550,7 +557,7 @@ tt2 <- ttheme_minimal(
(topall.leg <- cowplot::plot_grid(top.all,
tableGrob(trait.labs %>%
mutate(Code=1:n()) %>%
- dplyr::select(Code, Trait=Short_english_name) %>%
+ dplyr::select(Code, Trait=trait.name) %>%
filter(Trait %in% traits.sign.alone.pa),
theme=tt2, rows = NULL),
nrow=1, rel_widths=c(0.70, 0.3)))
@@ -570,14 +577,14 @@ ggsave(filename = "_pics/SXXX_Best_AllCombinations_CI_pa.png", dpi=400,
species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
#traits <- traits.backup
+traits <- traits %>%
+ rownames_to_column("species0") %>%
+ mutate_if(.predicate=~is.ordered(.),
+ .funs=~as.numeric(as.character(.)))
categorical.traits <- colnames(traits)[which(sapply(traits, "is.factor"))]
-#traits$Pollination <- factor(traits$Pollination,
-# levels=c("NEKTAR_HONIG_INSEKTEN","POLLEN","WIND" ),
-# labels=c("insects", "pollen", "wind"))
##spread categorical traits into multiple columns
-traits.dummy <- traits %>%
- rownames_to_column("species0")
+traits.dummy <- traits
while(any(sapply(traits.dummy, class)=="factor")){
id <- which(sapply(traits.dummy, class) == "factor")[1]
id.name <- names(id)
@@ -596,8 +603,8 @@ while(any(sapply(traits.dummy, class)=="factor")){
#print(id.name)
}
+
## match trait selection
-#e.g., traits.sign.alone <- c("SLA", "SeedMass", "CSR", "LeafPersistence")
CWM.wide <- species.cov %>%
pivot_longer(-RELEVE_NR, names_to = "species0", values_to = "rel.cov") %>%
@@ -609,7 +616,7 @@ CWM.wide <- species.cov %>%
by="species0") %>%
group_by(RELEVE_NR) %>%
select_if(.predicate=~is.numeric(.)) %>%
- summarize_at(.vars = vars(any_of(starts_with(as.character(traits.sign.alone)))),
+ summarize_at(.vars = vars(any_of(starts_with(as.character(traits.sign.alone.cov)))),
.funs = list(~weighted.mean(., rel.cov, na.rm=T))) %>%
dplyr::select(RELEVE_NR, order(colnames(.))) %>%
column_to_rownames("RELEVE_NR")
@@ -634,23 +641,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 ####
-#best.4traits <- c("Height", "BL_Dur", "VP_Fragm", "CSR")
+#best.traits.cov <- c("Height", "BL_Dur", "VP_Fragm", "CSR")
source("01b_MesobromionCluster.R")
W.fuzzy <- matrix.x(comm=species.cov %>%
column_to_rownames("RELEVE_NR"),
traits=traits %>%
- rownames_to_column("species0") %>%
filter(species0 %in% colnames(species.cov)) %>%
column_to_rownames("species0") %>%
- dplyr::select(all_of(traits.sign.alone)),
+ dplyr::select(all_of(traits.sign.alone.cov)),
scale = T)$matrix.X
-pca.fuzz <- rda(W.fuzzy)
+pca.fuzz <- rda(W.fuzzy, scale = F)
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:3], use = "pairwise.complete.obs") #double check
env.cor <- cor(env %>%
dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC),
- scores.pca[,1:3]) #double check
+ scores.pca[,1:3], use = "pairwise.complete.obs") #double check
myvectors <- as.data.frame(env.cor) %>%
rownames_to_column("mylab") %>%
@@ -658,12 +664,10 @@ 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.4traits, "bold", "italic")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.traits.cov, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
- mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
- rowwise() %>%
- mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab))
+ mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0))
basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
theme_bw() +
@@ -710,9 +714,9 @@ PCA.fuzz2_3 <- basemap0 +
ylab(paste("PC3 (", varexpl[3], "%)", sep=""))
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())
+ggsave("_pics/FigS991_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot())
-#### 4.0b Alternative showing species scores ####
+#### 4.0.1 Alternative showing species scores ####
tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>%
mutate(species0=colnames(species.cov)[-1]) %>%
mutate(species=species0) %>%
@@ -721,6 +725,18 @@ tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>%
mutate(Spe=substr(Spe, 1, 3)) %>%
mutate(labels=paste(Gen, Spe, sep="_"))
+fix.duplicate.labels <- function(x){
+ duplicated.labels.to.correct <- x %>%
+ dplyr::select(labels) %>%
+ mutate(ID=row_number()) %>%
+ group_by(labels) %>%
+ filter(n()>1) %>%
+ mutate(labels=paste0(labels, c("",2)))
+ x[duplicated.labels.to.correct$ID, "labels"] <- duplicated.labels.to.correct$labels
+ return(x)
+}
+
+tmp <- fix.duplicate.labels(tmp)
PCAfuzz1_2.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
@@ -729,9 +745,9 @@ PCAfuzz1_2.sp <- basemap0 %+% tmp +
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_point(data=myvectors %>%
- filter(categorical==1),
- aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
+# geom_point(data=myvectors %>%
+# filter(categorical==1),
+# aes(x=PC1, y=PC2, col=mycol), pch=16, size=1, alpha=0.8) +
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)) +
@@ -745,9 +761,9 @@ PCAfuzz2_3.sp <- basemap0 %+% tmp +
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_point(data=myvectors %>%
- filter(categorical==1),
- aes(x=PC2, y=PC3, col=mycol), pch=16, size=1, alpha=0.8) +
+# geom_point(data=myvectors %>%
+# filter(categorical==1),
+# aes(x=PC2, y=PC3, col=mycol), pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvectors,
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) +
@@ -755,8 +771,8 @@ PCAfuzz2_3.sp <- basemap0 %+% tmp +
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_2-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz2_3.sp)
+ggsave("_pics/FigS991a_PCA_Fuzzy_1-2_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_2.sp)
+ggsave("_pics/FigS991b_PCA_Fuzzy_2-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz2_3.sp)
@@ -767,28 +783,14 @@ W.beals <- as.data.frame(beals(species.cov %>%
column_to_rownames("RELEVE_NR"),
include=T, type=2))
write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv")
-pca.out <- rda(W.beals)
+pca.out <- rda(W.beals, scale=F)
varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
-cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:5)
-env.envfit <- envfit(pca.out, env %>%
- dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), choices=1:4, na.rm = T)
-
-### Transform to correlations and sink envfits
-### see https://www.davidzeleny.net/anadat-r/doku.php/en:suppl_vars_examples for procedure
-scores.pca <- pca.out$CA$u
-arrow_heads <- cwms.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
-r2 <- cwms.envfit$vectors$r # extracts vector of r2 for each env. variable
-cwms.cor <- arrow_heads * sqrt (r2)
-cor(CWM.wide, scores.pca[,1:5]) #double check
-
-
-arrow_heads <- env.envfit$vectors$arrows # extracts matrix of coordinates of arrow heads from ef
-r2 <- env.envfit$vectors$r # extracts vector of r2 for each env. variable
-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:3], scores.pca[,1:5])
+scores.pca <- pca.out$CA$u[,1:4]
+cwms.cor <- cor(CWM.wide, scores.pca)
+env.cor <- cor(env %>%
+ dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC),
+ scores.pca, use = "pairwise.complete.obs") #double check
+fuzz.cor <- cor(pca.fuzz$CA$u[,1:3], scores.pca)
sink("_output/S9_EnvFit_CWMs_env.txt")
cwms.cor
@@ -803,10 +805,11 @@ myvectors <- as.data.frame(env.cor) %>%
bind_rows(as.data.frame(cwms.cor) %>%
rownames_to_column("mylab") %>%
mutate(category="Trait")) %>%
- #bind_rows(as.data.frame(fuzz.cor) %>%
- # rownames_to_column("mylab") %>%
- # mutate(category="Fuzzy-Weighted")) %>%
- mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic")) %>%
+ bind_rows(as.data.frame(fuzz.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(mylab=paste0("FWM-", mylab)) %>%
+ mutate(category="Fuzzy-Weighted")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.traits.cov, "bold", "italic")) %>%
mutate(category=as.factor(category)) %>%
mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol)) %>%
@@ -864,7 +867,8 @@ PC_beals <- cowplot::plot_grid(PCA1_2,PCA3_4, nrow=1)
ggsave("_pics/Fig6_PC_Beals_1-4.png", width=10, height=5, dpi=300, last_plot())
-#### 4.1b Alternative showing species scores ####
+
+#### 4.1.1 Alternative showing species scores ####
tmp <- as.data.frame(pca.out$CA$v[,1:4]*5) %>%
mutate(species0=rownames(pca.out$CA$v)) %>%
mutate(species=species0) %>%
@@ -873,6 +877,7 @@ tmp <- as.data.frame(pca.out$CA$v[,1:4]*5) %>%
mutate(Spe=substr(Spe, 1, 3)) %>%
mutate(labels=paste(Gen, Spe, sep="_"))
+tmp <- fix.duplicate.labels(tmp)
PCA1_2.sp <- basemap0 %+% tmp +
geom_text(aes(x=PC1, y=PC2, label=labels), size=2, alpha=0.8) +
@@ -907,70 +912,171 @@ PCA3_4.sp <- basemap0 %+% tmp +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
-ggsave("_pics/Fig6a_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA1_2.sp)
-ggsave("_pics/Fig6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA3_4.sp)
+ggsave("_pics/Fig6a_PCA_Beals_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA1_2.sp)
+ggsave("_pics/Fig6b_PCA_Beals_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA3_4.sp)
+
+
+
+#### 4.2 RDA of Beals ~ FWMs ####
+RDA.beals <- rda(W.beals ~ scores(pca.fuzz, choices=1:3)$sites, scale=F)
+# var explained by CONSTRAINED axes
+varexpl <- round((RDA.beals$CCA$eig)/(sum(RDA.beals$CA$eig) + sum(RDA.beals$CCA$eig))*100,1)
+
+scores.rda <- scores(RDA.beals, choices = 1:3)$sites #
+#scores.rda <- RDA.beals$CA$u[,1:3]
+(cwms.cor <- cor(CWM.wide, RDA.beals$CCA$u[,1:3]))
+env.cor <- cor(env %>%
+ dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC),
+ RDA.beals$CCA$u[,1:3], use = "pairwise.complete.obs") #double check
+(fuzz.cor <- cor(pca.fuzz$CA$u, RDA.beals$CCA$u[,1:3])) #RDA.beals$CCA$biplot #
+
+myvectors.rda <- 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")) %>%
+ bind_rows(as.data.frame(fuzz.cor) %>%
+ rownames_to_column("mylab") %>%
+ mutate(mylab=paste0("FWM-", mylab)) %>%
+ mutate(category="Fuzzy-Weighted")) %>%
+ mutate(fontface0=ifelse(mylab %in% best.traits.cov, "bold", "italic")) %>%
+ mutate(category=as.factor(category)) %>%
+ mutate(mycol=ifelse(category=="Trait", oilgreen, orange)) %>%
+ mutate(mycol=ifelse(category=="Fuzzy-Weighted", myblue, mycol)) %>%
+ mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
+ rowwise() %>%
+ mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab))
+
+
+
+basemap0 <- ggplot(data=as.data.frame(scores.rda)) +
+ 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())
+
+
+RDA1_2 <- basemap0 +
+ geom_point(data=as.data.frame(scores.rda),
+ aes(x=RDA1, y=RDA2), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=myvectors.rda %>%
+ filter(categorical==0),
+ aes(x=0, xend=RDA1, y=0, yend=RDA2, 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.rda,
+ aes(x=RDA1, y=RDA2, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
+ xlab(paste("RDA1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("RDA2 (", varexpl[2], "%)", sep=""))
+
+ggsave("_pics/Fig6v2_RDA_Beals_1-2.png", width=10, height=5, dpi=300, RDA1_2)
+
+
+RDA1_3 <- basemap0 +
+ geom_point(data=as.data.frame(scores.rda),
+ aes(x=RDA1, y=RDA3), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=myvectors.rda %>%
+ filter(categorical==0),
+ aes(x=0, xend=RDA1, y=0, yend=RDA3, 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.rda,
+ aes(x=RDA1, y=RDA3, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
+ xlab(paste("RDA1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("RDA3 (", varexpl[3], "%)", sep=""))
+
+panel.RDA_beals <- cowplot::plot_grid(RDA1_2,RDA1_3, nrow=1)
+
+ggsave("_pics/Fig6v2_RDA_Beals_1-2-3.png", width=10, height=5, dpi=300, panel.RDA_beals)
+
+
+#### 4.2.1 Alternative showing species scores ####
+tmp <- as.data.frame(scores(RDA.beals, choices = 1:3)$species*4) %>%
+#tmp <- as.data.frame(RDA.beals$CCA$v*4) %>%
+ mutate(species0=rownames(RDA.beals$CCA$v)) %>%
+ 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="_"))
+
+tmp <- fix.duplicate.labels(tmp)
+
+
+RDA1_2.sp <- basemap0 %+% tmp +
+ geom_text(aes(x=RDA1, y=RDA2, label=labels), size=2, alpha=0.8) +
+ geom_segment(data=myvectors.rda%>%
+ filter(categorical==0),
+ aes(x=0, xend=RDA1, y=0, yend=RDA2, 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.rda,
+ aes(x=RDA1, y=RDA2, label=mylab, col=mycol, fontface=fontface0), size=2,
+ position = position_dodge(1), segment.alpha=0.5, segment.colour=gray(0.8)) +
+ xlab(paste("RDA1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("RDA2 (", varexpl[2], "%)", sep=""))
+
+RDA1_3.sp <- basemap0 %+% tmp +
+ geom_text(aes(x=RDA1, y=RDA3, label=labels), size=2, alpha=0.8) +
+ geom_segment(data=myvectors.rda%>%
+ filter(categorical==0),
+ aes(x=0, xend=RDA1, y=0, yend=RDA3, 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.rda,
+ aes(x=RDA1, y=RDA3, 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("RDA1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("RDA3 (", varexpl[3], "%)", sep=""))
+
+
+ggsave("_pics/Fig6v2a_RDA_Beals_1-2_wSpecies.png", width=8, height=8, dpi=300, RDA1_2.sp)
+ggsave("_pics/Fig6v2b_RDA_Beals_1-3_wSpecies.png", width=8, height=8, dpi=300, RDA1_3.sp)
+
-# #### 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.dummy %>%
+ column_to_rownames("species0") %>%
dplyr::select_if(.predicate = ~is.numeric(.)) %>% ### CAUTION -- ONLY NUMERIC TRAITS SHOWN
- dplyr::select(any_of(starts_with(as.character(traits.sign.alone))))
+ dplyr::select(any_of(starts_with(as.character(traits.sign.alone.cov))))
-corr1 <- cor(data2, use="complete.obs")
-pca2 <- princomp(covmat=corr1)
-#plot(pca2$loadings[,2]~pca2$loadings[,1])
-#text(pca2$loadings[,1], pca2$loadings[,2], dimnames(pca2$loadings)[[1]])
-#arrows(rep(0,8),rep(0,8),pca2$loadings[,1],pca2$loadings[,2], length=0.1)
+corr1 <- cor(data2, use="pairwise.complete.obs")
+#pca2 <- princomp(covmat=corr1)
zdat <- scale(data2) #this is just to standardize the original data, M = 0, SD =1
-e1 <- eigen(corr1) #solving for the eigenvalues and eigenvectors from the correlation matrix
+e1 <- eigen(corr1) #solving for the eigenvalues and eigenvectors from the correlation matrix
varexpl <- round((e1$values)/sum(e1$values)*100,1)
-pca.scores<- zdat %*% e1$vectors #scaled values x vectors
+zdat[which(is.na(zdat), arr.ind=T)] <- 0
+pca.scores <- zdat %*% e1$vectors #scaled values x vectors
pca.scores <- pca.scores[,1:5]
colnames(pca.scores)<-c('pca1','pca2','pca3','pca4','pca5') #just quickly naming the columns
dat <- circleFun(diameter = 2, npoints = 100)
-myvector3 <- as.data.frame(pca2$loadings[,1:5]) %>%
- rownames_to_column("mylab") %>%
+myvector3 <- as.data.frame(e1$vectors[,1:4]) %>%
+ rename(Comp.1=V1, Comp.2=V2, Comp.3=V3, Comp.4=V4) %>%
+ mutate(mylab=rownames(corr1)) %>%
+ #rownames_to_column("mylab") %>%
mutate(categorical=ifelse(grepl(pattern=paste(categorical.traits, collapse="|"), mylab), 1, 0)) %>%
- rowwise() %>%
- mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab)) %>%
- mutate(fontface0=ifelse(mylab %in% best.4traits, "bold", "italic"))
+ #rowwise() %>%
+ #mutate(mylab=gsub(pattern="LeafPersistence", replacement = "LeafPers", x = mylab)) %>%
+ mutate(fontface0=ifelse(mylab %in% best.traits.cov, "bold", "italic"))
baseplot <- ggplot(data=as.data.frame(pca.scores*.2)) +
scale_color_identity() +
coord_equal(xlim = c(-1, 1), ylim=c(-1, 1)) +
- xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
theme_bw() +
theme(panel.grid = element_blank())
@@ -981,12 +1087,14 @@ PCA.t1 <- baseplot +
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_point(data=myvector3 %>%
- filter(categorical==1),
- aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ #geom_point(data=myvector3 %>%
+ # filter(categorical==1),
+ # aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvector3,
aes(x=Comp.1, y=Comp.2, label=mylab, fontface=fontface0), col=oilgreen, size=2,
- position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
#ggtitle("PCoA of species-level traits")
PCA.t2 <- baseplot +
@@ -996,9 +1104,9 @@ PCA.t2 <- baseplot +
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_point(data=myvector3 %>%
- filter(categorical==1),
- aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ #geom_point(data=myvector3 %>%
+ # filter(categorical==1),
+ # aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvector3,
aes(x=Comp.3, y=Comp.4, label=mylab, fontface=fontface0), col=oilgreen, size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
@@ -1010,14 +1118,15 @@ PC_traits <- cowplot::plot_grid(PCA.t1, PCA.t2, nrow=1)
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) %>%
- mutate(species0=rownames(traits)) %>%
+tmp <- as.data.frame(pca.scores[,1:4]*.2) %>%
+ rownames_to_column("species0") %>%
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="_"))
+tmp <- fix.duplicate.labels(tmp)
PCA.t1.sp <- baseplot %+% tmp +
geom_text(aes(x=pca1, y=pca2, label=labels), size=2, alpha=0.7) +
@@ -1026,12 +1135,14 @@ PCA.t1.sp <- baseplot %+% tmp +
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_point(data=myvector3 %>%
- filter(categorical==1),
- aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ #geom_point(data=myvector3 %>%
+ # filter(categorical==1),
+ # aes(x=Comp.1, y=Comp.2),col=oilgreen, pch=16, size=1, alpha=0.8) +
geom_label_repel(data=myvector3,
aes(x=Comp.1, y=Comp.2, label=mylab, fontface=fontface0), col=oilgreen, size=2,
- position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8))
+ position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
PCA.t2.sp <- baseplot %+% tmp +
geom_text(aes(x=pca3, y=pca4, label=labels), size=2, alpha=0.7) +
@@ -1043,20 +1154,46 @@ PCA.t2.sp <- baseplot %+% tmp +
geom_label_repel(data=myvector3,
aes(x=Comp.3, y=Comp.4, label=mylab, fontface=fontface0),col=oilgreen, size=2,
position = position_dodge(2),segment.alpha=0.5, segment.colour=gray(0.8)) +
- geom_point(data=myvector3 %>%
- filter(categorical==1),
- aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
+ #geom_point(data=myvector3 %>%
+ # filter(categorical==1),
+ # aes(x=Comp.3, y=Comp.4),col=oilgreen, pch=16, size=1, alpha=0.8) +
xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
ylab(paste("PC4 (", varexpl[4], "%)", sep=""))
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)
+
+#traits.dummy %>% filter(species0 %in% (tmp %>% filter(labels %in% c("Fes_pal", "Ses_alb", "Car_hum")) %>% pull(species))) %>% dplyr::select(any_of(starts_with(as.character(traits.sign.alone.cov))))
+
+
+
## Create list of species labels
write_delim(tmp %>%
dplyr::select(species, label=labels),
path = "_output/S11_SpeciesList.txt")
+# #### depr 4.4 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")
@@ -1115,16 +1252,16 @@ ggsave(filename="_pics/S2_PlotDistribution.png", height
#### 6.1 Trait level #####
library(corrplot)
traits7 <- traits %>%
- dplyr::select(any_of(traits.sign.alone)) %>%
+ dplyr::select(any_of(traits.sign.alone.cov)) %>%
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(all_of(best.4traits), everything()) %>%
+ relocate(all_of(best.traits.cov), everything()) %>%
select_if(.predicate=~is.numeric(.)) ## CAUTION only numerical variables shown
## alternative using also dummy variables
#traits7 <- traits.dummy %>%
-# dplyr::select(any_of(starts_with(as.character(traits.sign.alone)))) %>%
+# dplyr::select(any_of(starts_with(as.character(traits.sign.alone.cov)))) %>%
# dplyr::select(sort(tidyselect::peek_vars())) %>%
-# relocate(all_of(starts_with(as.character(best.4traits))), everything())
+# relocate(all_of(starts_with(as.character(best.traits.cov))), everything())
res <- cor(traits7, use = "pairwise.complete.obs")
png(file="_pics/S7_Correlations_Trait.png", width=8, height=6.5, units = "in", res=300)
@@ -1135,9 +1272,9 @@ dev.off()
#### 6.2 CWM level ####
res2 <- cor(CWM.wide %>%
- dplyr::select(any_of(traits.sign.alone)) %>% ## caution selecting only numerical variables
+ dplyr::select(any_of(traits.sign.alone.cov)) %>% ## caution selecting only numerical variables
dplyr::select(sort(tidyselect::peek_vars())) %>%
- relocate(any_of(best.4traits), everything()))
+ relocate(any_of(best.traits.cov), 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)
--
GitLab