Preliminary results

01_Mesobromion.R

@@ -183,13 +183,90 @@ for(ff in myfilelist){
 corXY %>% 
   arrange(Test, desc(Coef))
 
-aa <- data.frame(Trait.comb=paste0("t", 1:95), trait.name=colnames(traits)[-which(colnames(traits) %in% c("species", "species0"))])
+aa <- data.frame(Trait.comb=paste0("t", 1:80), 
+                 trait.name=colnames(traits))#[-which(colnames(traits) %in% c("species", "species0"))])
 bb <- get.SES(obs.df = corXY, perm.df = corXY.perm, stat="RV") %>% 
-  left_join(aa, by="Trait.comb") %>% 
+  separate(Trait.comb, c("trait1", "trait2")) %>% 
+  mutate(trait2=paste0("t", trait2)) %>% 
+  left_join(aa %>% 
+              rename(trait1=Trait.comb), by="trait1") %>% 
+  dplyr::select(-trait1) %>% 
+  rename(trait1=trait.name) %>% 
+  left_join(aa %>% 
+              rename(trait2=Trait.comb), by="trait2") %>% 
+  dplyr::select(-trait2) %>% 
+  rename(trait2=trait.name) %>% 
+  dplyr::select(trait1, trait2, q025:conf.p) %>% 
   arrange(desc(SES.np))
 print(bb, n=20)
 
 
+#### Calculate and plot PCA of CWMs
+CWM.wide <- species %>% 
+  rownames_to_column("RELEVE_NR") %>% 
+  reshape2::melt(.id="RELEVE_NR") %>% 
+  rename(species0=variable, pres=value) %>% 
+  as.tbl() %>% 
+  filter(pres>0) %>% 
+  arrange(RELEVE_NR)  %>% 
+  ## attach traits 
+  left_join(traits %>% 
+              rownames_to_column("species0"), 
+            by="species0") %>% 
+  group_by(RELEVE_NR) %>% 
+  summarize_at(.vars = vars(LEB_F_Makrophanerophyt:Disp.unit.leng), 
+            .funs = list(~weighted.mean(., pres, na.rm=T))) %>%
+  dplyr::select(RELEVE_NR, order(colnames(.))) %>% 
+  column_to_rownames("RELEVE_NR")
+
+#PCA of CWMs
+CWM.pca <- vegan::rda(CWM.wide, scale=T)
+varexpl <- round(CWM.pca$CA$eig/sum(CWM.pca$CA$eig)*100,1)
+
+ggplot() + 
+  geom_point(data=as.data.frame(CWM.pca$CA$u[,1:2]), 
+             aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) + 
+  geom_segment(data=as.data.frame(CWM.pca$CA$v*2) %>% 
+                 rownames_to_column("Trait") %>% 
+                 mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), 
+               aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+  geom_text(data=as.data.frame(CWM.pca$CA$v*2.1) %>% 
+              rownames_to_column("Trait") %>% 
+              mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), 
+            aes(x=PC1, y=PC2, label=Trait, col=in.try),  size=3 ) +
+  scale_color_identity() + 
+  xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + 
+  ylab(paste("PC2 (", varexpl[2], "%)", sep="")) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank()) + 
+  title("PCA of CWMs")
+
+#PCA of individual traits
+trait.pca <- vegan::rda(as.matrix(traits %>% 
+                                    filter(!is.na(rowSums(.)))), na.action=na.omit)
+varexpl.t <- round(trait.pca$CA$eig/sum(trait.pca$CA$eig)*100,1)
+
+ggplot() + 
+  geom_point(data=as.data.frame(trait.pca$CA$u[,1:2]), 
+             aes(x=PC1, y=PC2), pch="+", size=2, alpha=0.8) + 
+  geom_segment(data=as.data.frame(trait.pca$CA$v*2) %>% 
+                 rownames_to_column("Trait") %>% 
+                 mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), 
+               aes(x=0, xend=PC1, y=0, yend=PC2, col=in.try), arrow = arrow(length = unit(0.08, "inches")), alpha=0.8) +
+  geom_text(data=as.data.frame(trait.pca$CA$v*2.1) %>% 
+              rownames_to_column("Trait") %>% 
+              mutate(in.try=ifelse(Trait %in% colnames(alltry), "Dark red", "Dark blue")), 
+            aes(x=PC1, y=PC2, label=Trait, col=in.try),  size=3 ) +
+  scale_color_identity() + 
+  xlab(paste("PC1 (", varexpl[1], "%)", sep="")) + 
+  ylab(paste("PC2 (", varexpl[2], "%)", sep="")) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank()) + 
+  title("PCA of species-level traits")
+
+
+
+
 #### Map of plots
 library(rgdal)
 library(sp)

98_SummarizeSimulations.R

@@ -3,7 +3,7 @@
 ## It then plots the summarized output
 library(tidyverse)
 
-mypath <- "_data/Experiment_27Feb2020"
+mypath <- "_data/Experiment_29Feb2020"
 
 myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T)
 output <- NULL
@@ -27,7 +27,7 @@ for(ff in myfiles){
 }
 
 outp.summary <- output %>% 
-  filter(!stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
+  dplyr::filter(!stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
   group_by(main, inter, feedb, trait, envir, stat.type) %>% 
   summarize(stat.obs.med=median(stat.obs),
             power=mean(pvalue<=0.05),
@@ -35,7 +35,7 @@ outp.summary <- output %>%
             exp.med.med=median(exp.med), 
             nsim=n()) %>% 
   bind_rows(output %>% 
-              filter(stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
+              dplyr::filter(stat.type %in% c("XY", "XY.T", "XY.TR")) %>% 
               group_by(main, inter, feedb, trait, stat.type) %>% 
               summarize(stat.obs.med=median(stat.obs),
                         power=mean(pvalue<=0.05),
@@ -46,33 +46,78 @@ outp.summary <- output %>%
   arrange(stat.type, main, inter, feedb, trait, envir)
 
 
-## plotting XY
+## plotting power for XY with feedback
 
 ggplot(data=outp.summary %>% 
          ungroup() %>% 
-         filter(stat.type=="XY") %>% 
-         filter(trait %in% c("1", "2", "1 2", "3")) %>% 
+         dplyr::filter(stat.type=="XY") %>% 
+         filter(feedb==0) %>% 
+         #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>% 
+         dplyr::filter(trait %in% c("1", "2", "3")) %>% 
          mutate(inter=as.factor(inter))) + 
   geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
   scale_colour_brewer(palette = "Dark2") + 
-  facet_grid(feedb~inter) + 
+  facet_grid(.~inter) + 
   theme_bw() + 
   theme(panel.grid = element_blank())
 
+## plotting non-parametric SES for XY with feedback
+
+ggplot(data=outp.summary %>% 
+         ungroup() %>% 
+         dplyr::filter(stat.type=="XY") %>% 
+         filter(feedb==0) %>% 
+         #dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>% 
+         dplyr::filter(trait %in% c("1", "2", "3")) %>% 
+         mutate(inter=as.factor(inter))) + 
+  geom_line(aes(x=main, y=SES.med, group=trait, col=trait)) + 
+  scale_colour_brewer(palette = "Dark2") + 
+  facet_grid(.~inter) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank())
+
+
 
 ## plotting XE
 
 ggplot(data=outp.summary %>% 
          ungroup() %>% 
-         filter(stat.type=="XE") %>% 
-         filter(envir=="1 ") %>% 
+         dplyr::filter(stat.type=="XE") %>% 
+         #dplyr::filter(envir=="1") %>% 
+         filter(feedb==0) %>% 
+         filter(trait %in% c("1","2","3", "1 2")) %>% 
          mutate(inter=as.factor(inter))
        #mutate(group0=paste("t", trait, " - e", envir))
 ) + 
-  geom_line(aes(x=main, y=power, group=trait, col=trait), lwd=1.5, alpha=0.3) + 
+  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
+  scale_colour_brewer(palette = "Dark2") + 
+  facet_grid(envir~inter) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank())
+
+## plotting XY.T
+
+ggplot(data=outp.summary %>% 
+         ungroup() %>% 
+         dplyr::filter(stat.type=="XY.T") %>% 
+         dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+         mutate(inter=as.factor(inter))) + 
+  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
   scale_colour_brewer(palette = "Dark2") + 
   facet_grid(feedb~inter) + 
   theme_bw() + 
   theme(panel.grid = element_blank())
 
+## plotting XY.TR
+
+ggplot(data=outp.summary %>% 
+         ungroup() %>% 
+         dplyr::filter(stat.type=="XY.TR") %>% 
+         dplyr::filter(trait %in% c("1", "2", "1 2", "3")) %>%
+         mutate(inter=as.factor(inter))) + 
+  geom_line(aes(x=main, y=power, group=trait, col=trait)) + 
+  scale_colour_brewer(palette = "Dark2") + 
+  facet_grid(feedb~inter) + 
+  theme_bw() + 
+  theme(panel.grid = element_blank())