Implemented alternative strategy - bootstrap + paired trait permutation

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=""))
 }

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)