Switched to RD

01_Mesobromion.R

@@ -160,7 +160,7 @@ write_delim(traits, path="_data/Mesobromion/traits.out.10perc.txt", delim="\t")
 ############################
 ## 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")
+traits <- read_delim("_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t")
 
 traits <- traits %>% 
   column_to_rownames("species0") # %>% 
@@ -169,7 +169,7 @@ traits <- traits %>%
 
 
 #### ## Import output #### 
-myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN(\\.[A-Za-z]+|_[A-Za-z]+)_[0-9]+.RData", full.names = T)
+myfilelist <- list.files(path="_derived/Mesobromion/", pattern="HIDDEN_signpair_[0-9]+.RData", full.names = T)
 dataFiles = purrr::map(myfilelist, function(x){get(load(x))}) 
 corXY = bind_rows(dataFiles) %>% 
   as.tbl()
@@ -205,6 +205,15 @@ corXY.ci <- corXY %>%
   arrange(ntraits, Coef.obs) %>% 
   dplyr::select(Trait.comb, Test, n, ntraits, everything())
   
+corXY.ci.sign <- corXY.ci 
+corXY.ci.pair <- corXY.ci
+
+corXY.ci <- bind_rows(corXY.ci.sign, corXY.ci.pair)
+
+corXY.ci %>% 
+  arrange(ntraits, Trait.comb, Coef.obs) %>% 
+  dplyr::select(Trait.comb, Test, n, ntraits, everything())
+
 traits.sign.alone <- corXY.ci %>% 
   filter(ntraits==1) %>% 
   filter(sign_plus==T) %>% 
@@ -218,24 +227,37 @@ ggplot(data=mydata) +
                    col=sign_plus)) + 
   geom_point(aes(x=Coef.obs, y=seq), pch=15) + 
   scale_y_continuous(breaks=mydata$seq, 
-                     labels=mydata$Trait.comb) 
+                     labels=mydata$trait1)  +  
+  scale_x_continuous(name="RV correlation")
 
 
-mydata <- mydata %>% 
+mydata <- corXY.ci %>% 
   #filter(ntraits==2)%>% 
-  filter(sign_plus==T) %>% 
-  filter(trait1 %in% traits.sign.alone & (is.na(trait2) | trait2 %in% traits.sign.alone)) %>% 
+  #filter(sign_plus==T) %>% 
+  #filter(trait1 %in% traits.sign.alone & (is.na(trait2) | trait2 %in% traits.sign.alone)) %>% 
   mutate(seq=1:nrow(.)) # %>% 
 #  arrange(desc(obs)) %>% 
 #  slice(1:20)
-ggplot(data=mydata) + 
+top7 <- ggplot(data=mydata %>% 
+                 rename(Significant=sign_plus)) + 
   geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, 
-                   col=sign_plus)) + 
+                   col=Significant)) + 
   geom_point(aes(x=Coef.obs, y=seq), pch=15) + 
   scale_y_continuous(breaks=mydata$seq, 
                      labels=mydata$Trait.comb, 
-                     name="RV correlation") 
-
+                     name="Trait combination") +  
+  scale_x_continuous(name="RV correlation") + 
+  theme(legend.position="bottom", legend.box = "horizontal")
+
+library(gridExtra)
+library(grid)
+tt2 <- ttheme_minimal()
+top7.leg <- cowplot::plot_grid(top7,
+                         tableGrob(trait.labs %>% dplyr::select(trait.name), theme=tt2), 
+                         nrow=1, rel_widths=c(0.65, 0.35))
+
+ggsave(filename = "_pics/Top7Predictors_CI.png", dpi=400, 
+       width=8, height=11, top7.leg)
 
 
 

01b_MesobromionCluster.R

@@ -66,7 +66,8 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
     #bootstrap species matrix 
 #    set.seed(1984)
     corXY.output <- rbind(corXY.output, 
-                          get.corXY.bootstrap(comm=species, trait=traits, trait.sel=tt, stat="RV", bootstrap=nperm))
+                          get.corXY.bootstrap(comm=species, trait=traits, 
+                                              trait.sel=tt, bootstrap=nperm))
     save(corXY.output, file = paste(output, "_", chunk.i, ".RData", sep=""))
   }
 }

99_HIDDEN_functions.R

@@ -11,6 +11,7 @@ library(SYNCSA)
 library(vegan)
 library(abind)
 library(ade4)
+library(energy)
 
 
 #### Function 1 - CorXY ####
@@ -48,7 +49,7 @@ get.corXY <- function(comm, traits, trait.sel="all", stat=c("mantel", "RV", "pro
 
 
 #### Function 1b - CorXY bootstrap####
-get.corXY.bootstrap <- function(comm, traits, trait.sel="all", stat="RV", bootstrap=199){
+get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){
   if(identical(trait.sel, "all")) {trait.sel <- 1:ncol(traits)}
   ii <- trait.sel
   lab.tmp <- paste(ii, collapse="_")
@@ -68,18 +69,24 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", stat="RV", bootst
   syn.out.perm.tmp <- matrix.x(comm=comm, traits=traits[index.traits[[bootstrap+1]],ii,drop=F], scale=T)$matrix.X
   
   corXY <- NULL
-  RV.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.tmp))
-  RV.perm.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.perm.tmp))
+  #RD.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.tmp), nrepet = 0)
+  RD.tmp <- dcor(W.beals, as.data.frame(syn.out.tmp))^2
+  #RV.perm.tmp <- RV.rtest(W.beals, as.data.frame(syn.out.perm.tmp), nrepet = 0)
+  RD.perm.tmp <- dcor(W.beals, as.data.frame(syn.out.perm.tmp))^2
   corXY <- rbind(corXY,
-                 data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RV", Coef.obs=RV.tmp$obs, Coef.perm=RV.perm.tmp$obs))
+                 data.frame(Trait.comb=lab.tmp, bootstr.n=0, Test="RD", 
+                            Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp))
   
   index.bootstr <- lapply(1:bootstrap, function(x){sample(1:n.sites, replace=T)})
   for(b in 1:bootstrap){
-    RV.tmp <- RV.rtest(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.tmp)[index.bootstr[[b]],])
+    #RV.tmp <- RV.rtest(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.tmp)[index.bootstr[[b]],])
+    RD.tmp <- dcor(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.tmp)[index.bootstr[[b]],])^2
     syn.out.perm.tmp <- matrix.x(comm=comm, traits=traits[index.traits[[bootstrap]],ii,drop=F], scale=T)$matrix.X
-    RV.perm.tmp <- RV.rtest(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.perm.tmp)[index.bootstr[[b]],])
+    #RV.perm.tmp <- RV.rtest(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.perm.tmp)[index.bootstr[[b]],])
+    RD.perm.tmp <- dcor(W.beals[index.bootstr[[b]],], as.data.frame(syn.out.perm.tmp)[index.bootstr[[b]],])^2
     corXY <- rbind(corXY,
-                   data.frame(Trait.comb=lab.tmp, bootstr.n=b, Test="RV", Coef.obs=RV.tmp$obs, Coef.perm=RV.perm.tmp$obs))
+                   data.frame(Trait.comb=lab.tmp, bootstr.n=b, Test="RD", 
+                              Coef.obs=RD.tmp, Coef.perm=RD.perm.tmp))
     print(paste("trait", paste0(ii, collapse="_"), "perm", b))
   }
   return(corXY)

session.R

@@ -3,8 +3,8 @@ source("01b_MesobromionCluster.R")
 species.path <- "_data/Mesobromion/species.out.10perc.txt"
 traits.path  <- "_data/Mesobromion/traits.out.10perc.txt"
 output  <- "test"
-myfunction <- "get.corXY"
-max.inter.t <- 2
+myfunction <- "get.corXY.bootstrap"
+max.inter.t <- 1
 chunkn <- 40
 chunk.i <- 1
 nperm <- 3