Updated Fig S6, Created correlation graphs

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()
+
+