New figure version in 01_Mesobromion.R

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)

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")