From 77cc017eb6797576889fd84ccf8d30286f3e2631 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Wed, 22 Jul 2020 15:33:25 +0200
Subject: [PATCH] New figure version in 01_Mesobromion.R
---
01_Mesobromion.R | 172 +++++++++++++++++++++++++++++++++++++++++------
session.R | 6 +-
2 files changed, 155 insertions(+), 23 deletions(-)
diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index 1a4857f..4ac0a19 100644
--- a/01_Mesobromion.R
+++ b/01_Mesobromion.R
@@ -301,7 +301,7 @@ write_delim(traits_nozero, path="_data/Mesobromion/traits.out.10perc_nozero.txt"
## calculate corr between species composition matrix and traits
species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
traits <- read_delim("_data/Mesobromion/traits.out.10perc.txt", delim="\t")
-trait.labs <- read_delim("_data/Mesobromion/TraitLabels.csv", delim=",") %>%
+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")
@@ -338,7 +338,12 @@ trait.recap <- traits %>%
spread(key=feature, value = value) %>%
rename(`Range/Levels`=Levels) %>%
mutate(Variable=factor(Variable, levels=colnames(traits))) %>%
- arrange(Variable)
+ arrange(Variable) %>%
+ left_join(trait.labs %>%
+ dplyr::select(Short_english_name, Long_English_name),
+ by=c("Variable" = "Short_english_name")) %>%
+ dplyr::select(Code=Variable, Variable=Long_English_name, everything())
+
write_csv(trait.recap, path="_derived/Mesobromion/Traits.recap.csv")
#dplyr::select(LeafArea:Disp.unit.leng)
@@ -428,7 +433,7 @@ mydata.byone <- corXY.ci %>%
)
-ggsave(filename = "_pics/TopFirstPredictors_CI.png", dpi=400,
+ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI.png", dpi=400,
width=4, height=6, plot = top.first)
#### 2.3 (Huge) Plot with all trait combination ####
@@ -578,7 +583,7 @@ tt2 <- ttheme_minimal(
# dplyr::select(Code, Trait=trait.name) %>%
# slice(33:n()), theme=tt2),
nrow=1, rel_widths=c(0.70, 0.3)))
-ggsave(filename = "_pics/Best_AllCombinations_CI.png", dpi=400,
+ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI.png", dpi=400,
width=6, height=3, topall.leg)
@@ -588,6 +593,10 @@ break()
##### 3. Calculate CWM ####
+#species <- species %>%
+# dplyr::select(colnames(.)[which(colSums(.)!=0)])
+
+
CWM.wide <- species %>%
rownames_to_column("RELEVE_NR") %>%
reshape2::melt(.id="RELEVE_NR") %>%
@@ -671,6 +680,26 @@ ggsave("_pics/PC_Beals_1-4.png", width=10, height=5, dpi=300, last_plot())
+
+### Transform to correlations and sink envfits
+### see https://www.davidzeleny.net/anadat-r/doku.php/en:suppl_vars_examples for procedure
+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)
+
+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)
+
+
+sink("_derived/Mesobromion/EnvFit_toCor_CWMs_env.txt")
+rbind(cwms.cor, env.cor)
+sink()
+
+
+
+
+
#### 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)
@@ -695,15 +724,80 @@ ggplot() +
#### 4.3 PCoA 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")
+#myweights <- c(rep(1/9, 9), rep(1/17, 17), rep(1, 26))
+#traits.gowdis <- FD::gowdis(traits, w=myweights, ord="podani")
+data2 <- traits %>%
+ dplyr::select(any_of(traits.sign.alone))
+corr1 <- cor(data2, use="pairwise.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)
+#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
+varexpl <- round((e1$values)/sum(e1$values)*100,1)
+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
+
+PCA.t1 <- ggplot() +
+ geom_point(data=as.data.frame(pca.scores[,1:2]),
+ aes(x=pca1, y=pca2), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]*8),
+ aes(x=0, xend=Comp.1, y=0, yend=Comp.2, col="Dark green"),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_label(data=as.data.frame(jitter(pca2$loadings[,1:5]*8, factor = 300)) %>%
+ #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.1, y=Comp.2, label=Trait, col="Dark green", fontface=fontface0), size=2,
+ position = position_dodge(2)) +
+ scale_color_identity() +
+ coord_equal(xlim = c(-5, 5), ylim=c(-5, 5)) +
+ xlab(paste("PC1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PC2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank())# +
+#ggtitle("PCoA of species-level traits")
+
+PCA.t2 <- ggplot() +
+ geom_point(data=as.data.frame(pca.scores[,1:4]),
+ aes(x=pca3, y=pca4), pch="+", size=2, alpha=0.8) +
+ geom_segment(data=as.data.frame(pca2$loadings[,1:5]*6),
+ aes(x=0, xend=Comp.3, y=0, yend=Comp.4, col="Dark green"),
+ arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+ geom_label(data=as.data.frame(jitter(pca2$loadings[,1:5]*6, factor = 300)) %>%
+ #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,
+ position = position_dodge(2)) +
+ scale_color_identity() +
+ coord_equal(xlim = c(-5, 5), ylim=c(-5, 5)) +
+ xlab(paste("PC3 (", varexpl[3], "%)", sep="")) +
+ ylab(paste("PC4 (", varexpl[4], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank())
+
+PC_traits <- cowplot::plot_grid(PCA.t1, PCA.t2, nrow=1)
+
+ggsave("_pics/PCA_Traits_1-4_only11.png", width=10, height=5, dpi=300, PC_traits)
+
+
+
+
+traits.gowdis <- FD::gowdis(traits %>%
+ dplyr::select(any_of(traits.sign.alone)), ord="podani")
+
pcoa.species <- wcmdscale(traits.gowdis, eig=T)
varexpl <- round((pcoa.species$eig)/sum(pcoa.species$eig)*100,1)
traits.envfit <- envfit(pcoa.species, traits %>%
dplyr::select(any_of(traits.sign.alone)),
- na.rm = T, choices = 1:3)
+ na.rm = T, choices = 1:4)
+
+
-set.seed(14)
+set.seed(15)
PCoA.t1 <- ggplot() +
geom_point(data=as.data.frame(pcoa.species$points[,1:2]),
aes(x=Dim1, y=Dim2), pch="+", size=2, alpha=0.8) +
@@ -711,8 +805,9 @@ PCoA.t1 <- ggplot() +
aes(x=0, xend=Dim1, y=0, yend=Dim2, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(jitter(traits.envfit$vectors$arrows*0.22, factor = 300)) %>%
- rownames_to_column("Trait"),
- aes(x=Dim1, y=Dim2, label=Trait, col="Dark green"), size=2,
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
+ aes(x=Dim1, y=Dim2, label=Trait, col="Dark green", fontface=fontface0), size=2,
position = position_dodge(2)) +
scale_color_identity() +
coord_equal(xlim = c(-.25, .25), ylim=c(-.25, .25)) +
@@ -723,25 +818,48 @@ PCoA.t1 <- ggplot() +
#ggtitle("PCoA of species-level traits")
PCoA.t2 <- ggplot() +
- geom_point(data=as.data.frame(pcoa.species$points[,1:3]),
- aes(x=Dim1, y=Dim3), pch="+", size=2, alpha=0.8) +
+ geom_point(data=as.data.frame(pcoa.species$points[,1:4]),
+ aes(x=Dim3, y=Dim4), pch="+", size=2, alpha=0.8) +
geom_segment(data=as.data.frame(traits.envfit$vectors$arrows*0.2),
- aes(x=0, xend=Dim1, y=0, yend=Dim3, col="Dark green"),
+ aes(x=0, xend=Dim3, y=0, yend=Dim4, col="Dark green"),
arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
geom_label(data=as.data.frame(traits.envfit$vectors$arrows*0.22) %>%
- rownames_to_column("Trait"),
- aes(x=Dim1, y=Dim3, label=Trait, col="Dark green"), size=2,
+ rownames_to_column("Trait") %>%
+ mutate(fontface0=ifelse(Trait %in% best.5traits, "bold", "plain")),
+ aes(x=Dim3, y=Dim4, label=Trait, col="Dark green", fontface=fontface0), size=2,
position = position_dodge(2)) +
scale_color_identity() +
coord_equal(xlim = c(-.25, .25), ylim=c(-.25, .25)) +
- xlab(paste("PCo1 (", varexpl[1], "%)", sep="")) +
- ylab(paste("PCo3 (", varexpl[3], "%)", sep="")) +
+ xlab(paste("PCo3 (", varexpl[3], "%)", sep="")) +
+ ylab(paste("PCo4 (", varexpl[4], "%)", sep="")) +
theme_bw() +
theme(panel.grid = element_blank())
PC_traits <- cowplot::plot_grid(PCoA.t1, PCoA.t2, nrow=1)
-ggsave("_pics/PCoA_Traits.png", width=10, height=5, dpi=300, PC_traits)
+ggsave("_pics/PCoA_Traits_1-4_only11.png", width=10, height=5, dpi=300, PC_traits)
+
+
+tmp <- as.data.frame(pcoa.species$points[,1:2]) %>%
+ rownames_to_column("labels") %>%
+ separate(labels, sep="_", into=c("Gen", "Spe")) %>%
+ mutate(Gen=substr(Gen, 1, 3)) %>%
+ mutate(Spe=substr(Spe, 1, 3)) %>%
+ mutate(labels=paste(Gen, Spe, sep="_"))
+
+
+ggplot() +
+ geom_text(data=tmp,
+ aes(x=Dim1, y=Dim2, label=labels), size=2, alpha=0.7) +
+ scale_color_identity() +
+ coord_equal(xlim = c(-.25, .25), ylim=c(-.25, .25)) +
+ xlab(paste("PCo1 (", varexpl[1], "%)", sep="")) +
+ ylab(paste("PCo2 (", varexpl[2], "%)", sep="")) +
+ theme_bw() +
+ theme(panel.grid = element_blank())# +
+#ggtitle("PCoA of species-level traits")
+ggsave("_pics/ForUTE_PCoA_Traits_12.png", width=10, height=9, dpi=300, last_plot())
+
#### 4.4 PCoA of X (Fuzzy weighted) matrix ####
@@ -785,6 +903,18 @@ library(rgdal)
library(sp)
library(sf)
env <- read_delim(file = "_data/Mesobromion/env.10perc.txt", delim = "\t")
+
+env0 <- read_delim("_data/Mesobromion/GVRD_MES2_site.csv", delim = ",")
+header <- "/data/sPlot/users/Francesco/Project_11/Germany/_data/tvhabita.dbf"
+env.all <- env0 %>%
+ left_join(foreign::read.dbf(header) %>%
+ as.data.frame() %>%
+ dplyr::select(RELEVE_NR, LAT, LON),
+ by="RELEVE_NR") %>%
+ filter(!is.na(LAT)) %>%
+ filter(!(LAT==0 | LON==0))
+
+
DE_NUTS1 <- readOGR(dsn="/data/sPlot/users/Francesco/Project_11/Germany/_data/Spatial", layer="Germany_LVL1_WGS84")
DE_NUTS_sf <- DE_NUTS1 %>%
st_as_sf()
@@ -819,7 +949,8 @@ ggsave(filename="_pics/PlotDistribution.png", height=8, width=6, dpi=300, basema
library(corrplot)
traits11 <- traits %>%
dplyr::select(any_of(traits.sign.alone)) %>%
- select(sort(tidyselect::peek_vars()))
+ select(sort(tidyselect::peek_vars())) %>%
+ relocate(all_of(best.5traits), everything())
res <- cor(traits11, use = "pairwise.complete.obs")
png(file="_pics/Correlations_Trait.png", width=8, height=6.5, units = "in", res=300)
@@ -830,7 +961,8 @@ dev.off()
#### 6.2 CWM level ####
res <- cor(CWM.wide %>%
- select(sort(tidyselect::peek_vars())))
+ select(sort(tidyselect::peek_vars())) %>%
+ relocate(all_of(best.5traits), everything()))
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)
diff --git a/session.R b/session.R
index f7a2323..0c16f8d 100644
--- a/session.R
+++ b/session.R
@@ -3,13 +3,13 @@ species.path <- "_data/Mesobromion/species.out.10perc.txt"
traits.path <- "_data/Mesobromion/traits.out.10perc.txt"
output <- "_derived/Mesobromion/HIDDEN"
myfunction <- "get.corXY.bootstrap"
-max.inter.t <- 7
+max.inter.t <- 2
chunk.i <- NA
-nperm <- 99
+nperm <- 5
ncores <- 8
chunkn <- 3*ncores
combinations <- "sequential"
-start.round <- 2
+start.round <- 1
relax.round <- 2
source("01b_MesobromionCluster.R")
--
GitLab