Minor fix to pictures and ordinations

02_Mesobromion_ExamineOutput.R

@@ -72,6 +72,9 @@ is.binary <- function(x){(all(na.omit(x) %in% 0:1))}
 
 ####1.1 Table S3 - Trait recap #####
 trait.recap <- traits %>% 
+  ##back transform log-transformed traits from TRY
+  mutate_at(.vars=vars(LA:DispUnitLength),
+            .funs=~exp(.)) %>% 
   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", 
@@ -82,7 +85,7 @@ trait.recap <- traits %>%
                                ifelse(is.binary(.),
                                       paste(names(table(.)), "=", table(.), collapse="; "),
                                       ifelse(is.numeric(.), 
-                                             paste(range(., na.rm=T), collapse=" - "),
+                                             paste(round(quantile(.,probs = c(0.05,.95), na.rm=T),1), collapse=" - "),
                                              ifelse(is.factor(.), 
                                                     paste(levels(.), collapse=", "),
                                                     NA)))))) %>% 
@@ -97,6 +100,21 @@ trait.recap <- traits %>%
             by=c("Variable" = "trait.name")) %>% 
   dplyr::select(Code=Variable, Variable=Long_English_name, everything())
 
+trait.recap <- trait.recap %>% 
+  mutate(Variable=ifelse(Code %in% c("FP_Beg", "FP_End", "FP_Dur"), paste0(Variable, " (month)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="LA", paste0(Variable, " (mm²)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="SLA", paste0(Variable, " (m2kg-1)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code %in% c("LeafNPratio", "LDMC"), paste0(Variable, " (g g-1)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code %in% c("LeafN","LeafP", "LeafCdrymass"), paste0(Variable, " (mg g-1)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="Height", paste0(Variable, " (m)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="SeedMass", paste0(Variable, " (mg)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code %in% c("SeedLength","DispUnitLength"), paste0(Variable, " (mm)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="LeafNarea", paste0(Variable, " (g m-2)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="Leaf.delta.15N", paste0(Variable, " (ppm)"), Variable)) %>% 
+  mutate(Variable=ifelse(Code=="LeafFreshMass", paste0(Variable, " (g)"), Variable))
+  
+
+  
 ## calculate trait coverage across plots, both on p\a and cover
 trait.coverage <- species.cov %>% 
   pivot_longer(values_to="cover", names_to = "species", -RELEVE_NR) %>% 
@@ -128,29 +146,29 @@ trait.recap <- trait.recap %>%
 
 write_csv(trait.recap, path="_output/S3_Traits.recap.csv")  
 
-
+###### ___________________ ######
 ###### 2. Import output from Cluster #### 
 ##### 2.1 Cover values ######
 # ## traits onebyone
- myfilelist0 <- list.files(path="_derived/Mesobromion/Cover/onebyone/", pattern="HIDDENcov-nona2_[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_tibble() %>% 
-   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$trait.name)) %>% 
-   arrange(ntraits, desc(Coef.obs)) %>% 
-   dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>% 
-   #mutate(Trait.comb=NA) %>% 
-   mutate(run="alone") %>% 
-   arrange(desc(sign_plus, Coef.obs))
-  
+# myfilelist0 <- list.files(path="_derived/Mesobromion/Cover/onebyone/", pattern="HIDDENcov-nona2_[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_tibble() %>% 
+#   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$trait.name)) %>% 
+#   arrange(ntraits, desc(Coef.obs)) %>% 
+#   dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>% 
+#   #mutate(Trait.comb=NA) %>% 
+#   mutate(run="alone") %>% 
+#   arrange(desc(sign_plus, Coef.obs))
+#  
 
 ### sequential trait combo
 myfilelist1 <- list.files(path="_derived/Mesobromion/Cover/", pattern="HIDDENcov_round_[0-9]+.RData", full.names = T)
@@ -181,7 +199,7 @@ 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.cov <- corXY.all.ci %>% 
+traits.sign.alone.cov <- corXY.ci %>% 
   #filter(run=="alone") %>% 
   filter(ntraits==1) %>% 
   filter(sign_plus==T) %>% 
@@ -194,12 +212,13 @@ mydata.byone <- corXY.ci %>%
   #filter(run=="alone") %>% 
   filter(ntraits==1)%>% 
   arrange(Coef.obs) %>% 
-  mutate(seq=1:n())
+  mutate(seq=1:n()) %>% 
+  left_join(trait.labs, by=c("trait1"="trait.name"))
 
 write_csv(mydata.byone %>% 
             dplyr::select(Trait.comb:sign_plus) %>% 
             arrange(desc(sign_plus), desc(Coef.obs)), 
-          path = "_output/S992_SolutionsOneByOIne.cov.csv")
+          path = "_output/S992_SolutionsOneByOne.pa.csv")
 
 (top.first <- ggplot(data=mydata.byone %>% 
                        mutate(sign_plus=fct_rev(as.factor(sign_plus)))) + 
@@ -207,7 +226,7 @@ write_csv(mydata.byone %>%
                    col=sign_plus)) + 
   geom_point(aes(x=Coef.obs, y=seq), pch=15) + 
   scale_y_continuous(breaks=mydata.byone$seq, 
-                     labels=mydata.byone$trait1, name=NULL)  +  
+                     labels=mydata.byone$Long_English_name , name=NULL)  +  
   scale_x_continuous(name="RD correlation") + 
   labs(color = "Significant") + 
   theme_bw() + 
@@ -218,7 +237,7 @@ write_csv(mydata.byone %>%
 
 
 ggsave(filename = "_pics/Fig4_TopFirstPredictors_CI_cov.png", dpi=400, 
-       width=4, height=6, plot = top.first)
+       width=5, height=6, plot = top.first)
 
 mydata <- corXY.ci %>% 
   filter(run=="seq") %>% 
@@ -348,30 +367,40 @@ tt2 <- ttheme_minimal(
   colhead = list(fg_params=list(cex = .7)),
   rowhead = list(fg_params=list(cex = .7)))
 
-(topall.leg <- cowplot::plot_grid(top.all,
-                                 tableGrob(trait.labs %>% 
+
+ttlabs <- trait.labs %>% 
   mutate(Code=1:n()) %>% 
-                                             dplyr::select(Code, Trait=trait.name) %>% 
-                                             filter(Trait %in% traits.sign.alone.cov),
-                                             theme=tt2, rows = NULL), 
-                                 #tableGrob(trait.labs %>% 
-                                #             mutate(Code=1:n()) %>% 
-                                #             dplyr::select(Code, Trait=trait.name) %>% 
-                                #             slice(33:n()), theme=tt2), 
-                         nrow=1, rel_widths=c(0.70, 0.3)))
+  filter(trait.name %in% traits.sign.alone.cov) 
+
+tobold <- which(ttlabs$trait.name %in% best.traits.cov)  
+tg <-  tableGrob(ttlabs %>% 
+                   dplyr::select(Code, Trait=Long_English_name) %>% 
+                   mutate(Trait=replace(x = Trait, 
+                                        list = Trait=="Vegetative Propagation - Fragmentation", 
+                                        values = "Veg. Propag. - Fragmentation")), 
+                 theme=tt2, rows = NULL)
+
+## this needs to be triple checked
+for (i in (11 + tobold)) {
+  tg$grobs[[i]] <- editGrob(tg$grobs[[i]], gp=gpar(fontface="bold"))
+}
+
+
+(topall.leg <- cowplot::plot_grid(top.all, tg,
+                         nrow=1, rel_widths=c(0.60, 0.4)))
 ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_cov.png", dpi=400, 
-       width=6, height=3, topall.leg)
+       width=6, height=2, topall.leg)
 
 
 break()
-
+###### ___________________ ######
 
 ###### 2.2 Presence/Absence data #### 
 ####### 2.2.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)
+# myfilelist0 <- list.files(path="_derived/Mesobromion/PresAbs/onebyone", pattern="HIDDENpa-nona2_[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) %>% 
@@ -383,10 +412,11 @@ break()
 #   mutate_at(.vars=vars(trait1),
 #             .funs=~factor(., 
 #                           levels=trait.labs$Trait.comb, 
-#                           labels=trait.labs$Short_english_name)) %>% 
+#                           labels=trait.labs$trait.name)) %>% 
 #   arrange(ntraits, Coef.obs) %>% 
 #   dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>% 
-#   mutate(Trait.comb=NA)
+#   arrange(desc(sign_plus), desc(Coef.obs)) 
+#   #mutate(Trait.comb=NA)
 
 
 ### combo of sign traits 
@@ -421,7 +451,7 @@ traits.sign.alone.pa <- corXY.ci.pa %>%
 ### 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         
+### NEW ### [1] Leaf_Scleroph Height        GF_Nanophan   SLA           LeafCdrymass  FP_Dur        GF_Hemiphan   LeafN         
 # LeafP         SeedLength   GF_Macrophan 
 
 ### Plot combinations one by one
@@ -577,7 +607,7 @@ ggsave(filename = "_pics/SXXX_Best_AllCombinations_CI_pa.png", dpi=400,
 
 
 
-
+###### ___________________ ######
 ##### 3. Calculate CWM ####
 #species <- read_delim("_data/Mesobromion/species.out.10perc.txt", delim="\t")
 species.cov <- read_delim("_data/Mesobromion/species.v2.10perc.cov.txt", delim="\t")
@@ -646,6 +676,7 @@ circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
   return(data.frame(x = xx, y = yy))
 }
 dat <- circleFun(diameter = 2, npoints = 100)
+###### ___________________ ######
 ### 4.0 PCA of X (Fuzzy-weighted) matrix + CWMs ####
 #best.traits.cov <- c("Height", "BL_Dur", "VP_Fragm", "CSR")
 source("01b_MesobromionCluster.R")
@@ -654,7 +685,7 @@ W.fuzzy <- matrix.x(comm=species.cov %>%
                     traits=traits %>% 
                       filter(species0 %in% colnames(species.cov)) %>% 
                       column_to_rownames("species0") %>% 
-                      dplyr::select(all_of(traits.sign.alone.cov)), 
+                      dplyr::select(all_of(best.traits.cov)), 
                     scale = T)$matrix.X
 pca.fuzz <- rda(W.fuzzy, scale = F) 
 varexpl <- round((pca.fuzz$CA$eig)/sum(pca.fuzz$CA$eig)*100,1)
@@ -690,9 +721,9 @@ basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
                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,# %>% 
                #mutate_if(~is.numeric(.), ~(.)*1.2),
              aes(x=PC1, y=PC2, label=mylab, col=mycol, fontface=fontface0), size=2, 
@@ -701,26 +732,26 @@ basemap0 <- ggplot(data=as.data.frame(scores.pca*5)) +
   ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
 )
 
-PCA.fuzz2_3 <- basemap0 + 
+PCA.fuzz1_3 <- basemap0 + 
   geom_point(data=as.data.frame(pca.fuzz$CA$u[,1:3]*5), 
-             aes(x=PC2, y=PC3), pch="+", size=2, alpha=0.8) + 
+             aes(x=PC1, y=PC3), pch="+", size=2, alpha=0.8) + 
   geom_segment(data=myvectors%>% 
                  filter(categorical==0), 
-               aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol), 
+               aes(x=0, xend=PC1, 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=PC1, y=PC3, col=mycol), pch=16, size=1, alpha=0.8) + 
   geom_label_repel(data=myvectors,# %>% 
                #mutate_if(~is.numeric(.), ~(.)*1.2),
-             aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2, 
+             aes(x=PC1, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2, 
              position = position_dodge(1), segment.alpha=0.8, segment.colour=gray(0.7)) +
-  xlab(paste("PC2 (", varexpl[2], "%)", sep="")) + 
+  xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + 
   ylab(paste("PC3 (", varexpl[3], "%)", sep=""))
 
-PC_fuzzy <- cowplot::plot_grid(PCA.fuzz1_2,PCA.fuzz2_3, nrow=1)
-ggsave("_pics/FigS991_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot())
+PC_fuzzy <- cowplot::plot_grid(PCA.fuzz1_2,PCA.fuzz1_3, nrow=1)
+ggsave("_pics/S11_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot())
 
 #### 4.0.1 Alternative showing species scores ####
 tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>% 
@@ -760,35 +791,35 @@ PCAfuzz1_2.sp <- basemap0 %+% tmp +
   xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + 
   ylab(paste("PC2 (", varexpl[2], "%)", sep=""))
 
-PCAfuzz2_3.sp <- basemap0 %+% tmp + 
-  geom_text(aes(x=PC2, y=PC3, label=labels), size=2, alpha=0.8) + 
+PCAfuzz1_3.sp <- basemap0 %+% tmp + 
+  geom_text(aes(x=PC1, y=PC3, label=labels), size=2, alpha=0.8) + 
   geom_segment(data=myvectors %>% 
                  filter(categorical==0), 
-               aes(x=0, xend=PC2, y=0, yend=PC3, col=mycol), 
+               aes(x=0, xend=PC1, 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_label_repel(data=myvectors, 
-                   aes(x=PC2, y=PC3, label=mylab, col=mycol, fontface=fontface0), size=2, 
+                   aes(x=PC1, 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) +
-  xlab(paste("PC2 (", varexpl[2], "%)", sep="")) + 
+  xlab(paste("PC1 (", varexpl[2], "%)", sep="")) + 
   ylab(paste("PC3 (", varexpl[3], "%)", sep="")) 
 
 
-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)
-
+ggsave("_pics/S11a_PCA_Fuzzy_1-2_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_2.sp)
+ggsave("_pics/S11b_PCA_Fuzzy_1-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_3.sp)
 
 
 
 
+###### _ ######
 #### 4.1 PCA of Y (Bealls) matrix + CWM ####
 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")
+#write.table(W.beals, sep="\t", file="_derived/Mesobromion/MatrixY_Beals.csv")
 pca.out <- rda(W.beals, scale=F)
 varexpl <- round((pca.out$CA$eig)/sum(pca.out$CA$eig)*100,1)
 scores.pca <- pca.out$CA$u[,1:4]
@@ -813,7 +844,7 @@ myvectors <- as.data.frame(env.cor) %>%
               mutate(category="Trait")) %>% 
   bind_rows(as.data.frame(fuzz.cor) %>% 
               rownames_to_column("mylab") %>%
-              mutate(mylab=paste0("FWM-", mylab)) %>% 
+              mutate(mylab=paste0("FW-", mylab)) %>% 
               mutate(category="Fuzzy-Weighted")) %>% 
   mutate(fontface0=ifelse(mylab %in% best.traits.cov, "bold", "italic")) %>% 
   mutate(category=as.factor(category)) %>% 
@@ -922,7 +953,7 @@ ggsave("_pics/Fig6a_PCA_Beals_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA
 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
@@ -999,7 +1030,6 @@ 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) %>% 
@@ -1049,7 +1079,7 @@ ggsave("_pics/Fig6v2b_RDA_Beals_1-3_wSpecies.png", width=8, height=8, dpi=300, R
 
 
 
-
+###### _ ######
 #### 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))
@@ -1203,7 +1233,7 @@ write_delim(tmp %>%
 
 
 
-
+###### ___________________ ######
 #### 5 Map of Plot distribution #####
 library(rgdal)
 library(sp)
@@ -1253,7 +1283,7 @@ ggsave(filename="_pics/S2_PlotDistribution.png", height
 
 
 
-
+###### ___________________ ######
 #### 6 Correlation matrices ####
 #### 6.1 Trait level #####
 library(corrplot)