From 9af7c6e29a9903feeff76fbae50a3dbd93d12b51 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Wed, 15 Apr 2020 16:00:53 +0200
Subject: [PATCH] Switched to RD
---
01_Mesobromion.R | 42 ++++++++++++++++++++++++++++++----------
01b_MesobromionCluster.R | 3 ++-
99_HIDDEN_functions.R | 21 +++++++++++++-------
session.R | 4 ++--
4 files changed, 50 insertions(+), 20 deletions(-)
diff --git a/01_Mesobromion.R b/01_Mesobromion.R
index 09c4102..2ad93f9 100644
--- a/01_Mesobromion.R
+++ b/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)
diff --git a/01b_MesobromionCluster.R b/01b_MesobromionCluster.R
index f3e1bf7..b7c56c4 100644
--- a/01b_MesobromionCluster.R
+++ b/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=""))
}
}
diff --git a/99_HIDDEN_functions.R b/99_HIDDEN_functions.R
index 00659b1..7d5580d 100644
--- a/99_HIDDEN_functions.R
+++ b/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)
diff --git a/session.R b/session.R
index a12c291..cf3d700 100644
--- a/session.R
+++ b/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
--
GitLab