From a554d5b7b78570ba1cfcbdeeca1795bfbd951432 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Wed, 17 Jun 2020 11:04:52 +0200
Subject: [PATCH] Updated Fig S6, Created correlation graphs
---
01_Mesobromion.R | 85 +++++++++++++++++++++++++++++++++---------------
1 file changed, 58 insertions(+), 27 deletions(-)
diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index c838847..1a4857f 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -317,12 +317,14 @@ colnames(traits) <- trait.labs$Short_english_name
##check traits
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)) %>%
summarize_all(.funs = list(xxxType.of.variable=~ifelse(is.binary(.), "binary",
ifelse(is.ordered(.), "ordered",
- ifelse(is.numeric(.), "numeric",
- ifelse(is.factor(.), "factor", NA)))),
+ ifelse(is.numeric(.), "quantitative",
+ ifelse(is.factor(.), "nominal", NA)))),
xxxLevels=~(
ifelse(is.binary(.),
paste(names(table(.)), "=", table(.), collapse="; "),
@@ -338,7 +340,6 @@ trait.recap <- traits %>%
mutate(Variable=factor(Variable, levels=colnames(traits))) %>%
arrange(Variable)
-## Shoudln't BL_ANAT better treated as ordered.factor?
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)
@@ -407,7 +408,7 @@ traits.sign.alone <- corXY.ci %>%
pull(trait1)
### Plot combinations one by one
-mydata.byone <- mydata %>%
+mydata.byone <- corXY.ci %>%
filter(ntraits==1)%>%
mutate(seq=1:n())
@@ -433,7 +434,7 @@ ggsave(filename = "_pics/TopFirstPredictors_CI.png", dpi=400,
#### 2.3 (Huge) Plot with all trait combination ####
### filter combinations where all traits belong to traits.sign.alone
mydata <- corXY.ci %>%
- filter_at(.vars=vars(trait1:trait22), all_vars(. %in% traits.sign.alone | is.na(.)))%>%
+ filter_at(.vars=vars(trait1:trait11), all_vars(. %in% traits.sign.alone | is.na(.)))%>%
mutate(seq=1:n())
allpredictors <- top.first %+% (mydata) +
@@ -514,6 +515,12 @@ for(nn in 1:maxtraits){
}
}
+best.5traits <- corXY.ci %>%
+ filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>%
+ dplyr::select(trait1:trait11) %>%
+ mutate_all(~as.character(.))
+best.5traits <- as.character(best.5traits[1,1:5])
+
### Create dataset with best combinations + all the one-way combinations
mydata.best <- mydata %>%
filter_at(.vars=vars(trait1:trait11),
@@ -600,15 +607,15 @@ CWM.wide <- species %>%
column_to_rownames("RELEVE_NR")
#### 4. PCA Graphs ####
-#### 4.1 PCA of Y (Bealls) matrix ####
+#### 4.1 PCA of Y (Bealls) matrix + CWM ####
library(vegan)
W.beals <- as.data.frame(beals(species, 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)
-cwms.envfit <- envfit(pca.out, CWM.wide, na.rm = T, choices = 1:3)
+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:3)
+ dplyr::select(Temp, Prec, pH=PHIPHOX, C.org=ORCDRC), choices=1:5)
PCA1_2 <- ggplot() +
geom_point(data=as.data.frame(pca.out$CA$u[,1:2]),
@@ -621,9 +628,10 @@ PCA1_2 <- ggplot() +
rownames_to_column("Env"),
aes(x=PC1, y=PC2, label=Env), col="Dark red", size=2,
position = position_dodge(1) ) +
- geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Trait"),
- aes(x=PC1, y=PC2, label=Trait), col="Dark blue", size=2,
+ geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
+ aes(x=PC1, y=PC2, label=Trait, fontface=fontface0), col="Dark blue", size=2,
position = position_dodge(1) ) +
xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
@@ -632,37 +640,37 @@ PCA1_2 <- ggplot() +
scale_x_continuous(limits=c(-0.25, 0.25)) + coord_equal() +
theme(panel.grid = element_blank())
-(PCA1_3 <- ggplot() +
- geom_point(data=as.data.frame(pca.out$CA$u[,1:3]),
- aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) +
+(PCA3_4 <- ggplot() +
+ geom_point(data=as.data.frame(pca.out$CA$u[,1:4]),
+ aes(x=PC3, y=PC4), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(cwms.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC1, y=0, yend=PC3), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ aes(x=0, xend=PC3, y=0, yend=PC4), col="Dark blue", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_segment(data=as.data.frame(env.envfit$vectors$arrows*.2),
- aes(x=0, xend=PC1, y=0, yend=PC3), col="Dark red", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ aes(x=0, xend=PC3, y=0, yend=PC4), col="Dark red", arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(env.envfit$vectors$arrows*.23) %>%
rownames_to_column("Env"),
- aes(x=PC1, y=PC3, label=Env), col="Dark red", size=2,
+ aes(x=PC3, y=PC4, label=Env), col="Dark red", size=2,
position = position_dodge(1)) +
- geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
- rownames_to_column("Trait"),
- aes(x=PC1, y=PC3, label=Trait), col="Dark blue", size=2,
- position = position_dodge(1) ) +
- xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PC3 (", varexpl[3], "%)", sep="")) +
+ geom_label(data=as.data.frame(cwms.envfit$vectors$arrows*.23) %>%
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, 2, 1)),
+ aes(x=PC3, y=PC4, label=Trait, fontface=fontface0 ), col="Dark blue", size=2,
+ position = position_dodge(1)
+ ) +
+ xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
+ ylab(paste("PC4 (", varexpl[4], "%)", sep="")) +
theme_bw() +
scale_y_continuous(limits=c(-0.25, 0.25)) +
scale_x_continuous(limits=c(-0.25, 0.25)) +
coord_equal() +
theme(panel.grid = element_blank())
)
-PC_beals <- cowplot::plot_grid(PCA1_2,PCA1_3, nrow=1)
+PC_beals <- cowplot::plot_grid(PCA1_2,PCA3_4, nrow=1)
-ggsave("_pics/PC_Beals.png", width=10, height=5, dpi=300, last_plot())
+ggsave("_pics/PC_Beals_1-4.png", width=10, height=5, dpi=300, last_plot())
-#### Old down from here
-
#### 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)
@@ -806,6 +814,29 @@ ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basema
+#### 6 Correlation matrices ####
+#### 6.1 Trait level #####
+library(corrplot)
+traits11 <- traits %>%
+ dplyr::select(any_of(traits.sign.alone)) %>%
+ select(sort(tidyselect::peek_vars()))
+
+res <- cor(traits11, use = "pairwise.complete.obs")
+png(file="_pics/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)
+dev.off()
+
+
+#### 6.2 CWM level ####
+res <- cor(CWM.wide %>%
+ select(sort(tidyselect::peek_vars())))
+png(file="_pics/Correlations_CWMs.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)
+dev.off()
+
+
--
GitLab