From dc436044826cc89471aacb11efea76085ed60571 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Fri, 20 Mar 2020 14:56:11 +0100
Subject: [PATCH] Implemented alternative strategy - bootstrap + paired trait
permutation
---
01b_MesobromionCluster.R | 50 ++++++++++++++++++++++------------------
session.R | 4 ++--
2 files changed, 29 insertions(+), 25 deletions(-)
diff --git a/01b_MesobromionCluster.R b/01b_MesobromionCluster.R
index 375929f..3b60cdb 100644
--- a/01b_MesobromionCluster.R
+++ b/01b_MesobromionCluster.R
@@ -170,7 +170,9 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
traits <- traits %>%
column_to_rownames("species0") %>%
- rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.))
+ rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.)) %>%
+#temporary ### Use only a subset of traits
+ dplyr::select(LeafArea:Disp.unit.leng)
## create list of indices for each combination of traits up to a max number of interactions
n.traits <- ncol(traits)
@@ -190,34 +192,36 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
print(paste("Running on", length(allcomb.t),"out of", nall, "combinations"))
names.t <- unlist(lapply(allcomb.t, paste, collapse="_"))
- print("Run on observed data")
- cor.obs <- NULL
+ print("Start main loop")
+ cor.obs <- NA
+ cor.bootstr <- matrix(NA, nrow=length(allcomb.t), ncol=nperm,
+ dimnames = list(paste("t", names.t, sep=""), paste("p",1:nperm, sep="")))
+ cor.perm <- matrix(NA, nrow=length(allcomb.t), ncol=nperm,
+ dimnames = list(paste("t", names.t, sep=""), paste("p",1:nperm, sep="")))
+
+ speciesb <- species
+ traitsb <- traits
for(i in 1:length(allcomb.t)) {
tt <- unlist(allcomb.t[i])
- cor.obs <- rbind(cor.obs,
- myfunction(comm=species, trait=traits, trait.sel=tt, stat="RV"))
- print(i)
- }
-
- print("Permute traits data")
- set.seed(1984)
- traits.perm <- lapply(1:nperm, FUN=function(x){
- tmp <- traits[sample(1:nrow(traits)),]
- rownames(tmp) <- rownames(traits)
- return(tmp)}
- )
-
- print("Run on permuted data")
- cor.perm <- matrix(NA, nrow=length(allcomb.t), ncol=nperm, dimnames = list(paste("t", names.t, sep=""),
- paste("p",1:nperm, sep="")))
- for(j in 1:length(allcomb.t)){
- tt <- unlist(allcomb.t[j])
+ #bootstrap species matrix
+# set.seed(1984)
for(b in 1:nperm){
- cor.perm[j,b] <- myfunction(comm=species, trait=traits.perm[[b]], trait.sel=tt, stat="RV")$Coef
+ if(b>1){
+ speciesb <- species[sample(1:nrow(species), replace=T),] #resample plots
+ speciesb <- speciesb[,-which(colSums(speciesb)==0)] # delete empty species
+ traitsb <- traits[which(rownames(traits) %in% colnames(speciesb)),] #subset traits
+ }
+ traitsb.perm <- traitsb[sample(1:nrow(traitsb)),] #permute traits
+ tmp <- myfunction(comm=speciesb, trait=traitsb, trait.sel=tt, stat="RV")
+ if(b==1) {cor.obs <- rbind(cor.obs, tmp)}
+ cor.bootstr[i,b] <- tmp$Coef
+ cor.perm[i,b] <- myfunction(comm=speciesb, trait=traitsb.perm, trait.sel=tt, stat="RV")$Coef
print(paste("trait", tt, "perm", b))
}
+ print(i)
}
+
#save(corXY, file="_data/Mesobromion/corXY/corXY.RData")
- save(cor.obs, cor.perm, file = paste(output, "_", chunk.i, ".RData", sep=""))
+ save(cor.obs, cor.bootstr, cor.perm, file = paste(output, "_", chunk.i, ".RData", sep=""))
}
diff --git a/session.R b/session.R
index 16ae1cb..a12c291 100644
--- a/session.R
+++ b/session.R
@@ -4,9 +4,9 @@ species.path <- "_data/Mesobromion/species.out.10perc.txt"
traits.path <- "_data/Mesobromion/traits.out.10perc.txt"
output <- "test"
myfunction <- "get.corXY"
-max.inter.t <- 1
+max.inter.t <- 2
chunkn <- 40
chunk.i <- 1
-nperm <- 99
+nperm <- 3
Mesobromion(species.path, traits.path, output, myfunction, max.inter.t, chunkn, chunk.i, nperm)
--
GitLab