From 680c334b9c622126d8d0ce5d4f89f3a1b9d5e1f8 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Tue, 18 Aug 2020 15:33:14 +0200
Subject: [PATCH] Modified script to account for NAs in matrix.x and trait
permutation
---
00_Mesobromion_DataPreparation.R | 112 ++++++++++++++----------
01b_MesobromionCluster.R | 116 ++++++++++++-------------
99_HIDDEN_functions.R | 145 +++++++++++++++++--------------
session.R | 6 +-
4 files changed, 207 insertions(+), 172 deletions(-)
diff --git a/00_Mesobromion_DataPreparation.R b/00_Mesobromion_DataPreparation.R
index d999dad..682828e 100644
--- a/00_Mesobromion_DataPreparation.R
+++ b/00_Mesobromion_DataPreparation.R
@@ -5,6 +5,8 @@ library(gridExtra)
library(grid)
library(vegan)
+
+
source("99_HIDDEN_functions.R")
##### PART 1 ####
@@ -240,7 +242,7 @@ recode.traits <- function(x){
return(out)
}
traits <- recode.traits(traits)
-
+
### ordered factors
@@ -292,7 +294,7 @@ env <- env %>%
species <- species %>%
rownames_to_column("RELEVE_NR")
-### 5. Transform and export data
+###### 5. Transform and export data ######
## presence absence
species.pa <- species
@@ -332,7 +334,7 @@ traits %>%
traits %>%
filter_at(.vars=vars(-"species0"),
any_vars(is.na(.))) %>%
- nrow() ## [1] 109 # species with no at least 1 NA in traits
+ nrow() ## [1] 109 # species with at least 1 NA in traits
#### CORRELATION BETWEEN FUZZY WEIGHTED AND BEALS MATRICES
@@ -346,9 +348,9 @@ traits %>%
### PART 2 ####
source("01b_MesobromionCluster.R")
-#### 1. Traits individually significant for COVER data####
-traits <- read_delim("_data/Mesobromion/traits.out.10perc.cov.txt", delim="\t")
-myfilelist <- list.files(path="_derived/Mesobromion/Cover", pattern="HIDDENproz_[0-9]+_.RData", full.names = T)
+#### 1. Traits individually significant for COVER data#### na.exclude=T ########
+traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
+myfilelist <- list.files(path="_derived/Mesobromion/Cover", pattern="HIDDENcov_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
as_tibble()
@@ -356,39 +358,43 @@ rm( dataFiles)
trait.labs <- data.frame(Trait.comb=as.character(1:(ncol(traits)-1)),
trait.name=colnames(traits)[-1])
-
+
corXY.ci <- get.ci(corXY)
corXY.ci <- corXY.ci %>%
arrange(desc(sign_plus), desc(Coef.obs)) %>%
left_join(trait.labs, by="Trait.comb") %>%
dplyr::select(-Test)
-## NO significant TRAITS when using Cover values
-# traits.sign.alone <- corXY.ci %>%
-# filter(sign_plus) %>%
-# pull(trait.name)
-#
-# traits.sign <- traits %>%
-# dplyr::select(species0, any_of(traits.sign.alone))
-# write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.cov.sign.txt", delim="\t")
+traits.sign.alone <- corXY.ci %>%
+ filter(sign_plus) %>%
+ pull(trait.name)
-"# A tibble: 50 x 11
- Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+traits.sign <- traits %>%
+ dplyr::select(species0, any_of(traits.sign.alone))
+#write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.cov.sign.txt", delim="\t")
+
+### COV - NONAs
+" Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 17 0.335 0.205 0.275 0.436 0.903 999 FALSE FALSE 1 V_VER_Fragmentation
- 2 12 0.255 0.139 0.235 0.315 0.690 999 FALSE FALSE 1 V_VER_Rhizom
- 3 11 0.226 0.170 0.218 0.279 0.533 999 FALSE FALSE 1 V_VER_Ausläufer
- 4 5 0.199 0.340 0.194 0.247 0.502 999 FALSE FALSE 1 LEB_F_Hemiphanerophyt
- 5 2 0.198 0.232 0.191 0.247 0.447 999 FALSE FALSE 1 LEB_F_Nanophanerophyt
- 6 19 0.198 0.217 0.193 0.245 0.539 999 FALSE FALSE 1 V_VER_Sprossknolle
- 7 44 0.194 0.171 0.199 0.242 0.299 999 FALSE FALSE 1 BLU_KL
- 8 45 0.192 0.188 0.187 0.238 0.526 999 FALSE FALSE 1 REPR_T
- 9 10 0.191 0.187 0.187 0.240 0.366 999 FALSE FALSE 1 V_VER_Wurzelspross
-10 42 0.190 0.235 0.192 0.236 0.261 999 FALSE FALSE 1 Disp.unit.leng
-# … with 40 more rows"
+ 1 36 0.316 0.154 0.290 0.370 0.997 999 TRUE FALSE 1 PlantHeight
+ 2 50 0.286 0.154 0.257 0.339 0.977 999 TRUE FALSE 1 BL_ANAT
+ 3 30 0.259 0.114 0.232 0.308 0.974 999 TRUE FALSE 1 BL_DAU
+ 4 2 0.256 0.228 0.222 0.303 0.992 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 5 20 0.253 0.0703 0.205 0.324 0.992 999 TRUE FALSE 1 V_VER_Fragmentation
+ 6 49 0.252 0.116 0.226 0.303 0.967 999 TRUE FALSE 1 STRAT_T
+ 7 32 0.251 0.127 0.226 0.303 0.968 999 TRUE FALSE 1 SLA
+ 8 35 0.241 0.128 0.217 0.289 0.970 999 TRUE FALSE 1 LeafP "
+
+### COV - without deleting NAs
+"# A tibble: 53 x 11
+ Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+ <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
+ 1 36 0.304 0.117 0.280 0.356 0.995 999 TRUE FALSE 1 PlantHeight
+ 2 32 0.247 0.127 0.221 0.299 0.976 999 TRUE FALSE 1 SLA
+ 3 35 0.241 0.175 0.218 0.288 0.951 999 TRUE FALSE 1 LeafP "
#### 2. Traits individually significant for Presence|absence data####
-traits <- read_delim("_data/Mesobromion/traits.out.10perc.txt", delim="\t")
+traits <- read_delim("_data/Mesobromion/traits.v2.10perc.txt", delim="\t")
myfilelist <- list.files(path="_derived/Mesobromion/PresAbs", pattern="HIDDENpa_[0-9]+_.RData", full.names = T)
dataFiles = purrr::map(myfilelist, function(x){get(load(x))})
corXY = bind_rows(dataFiles) %>%
@@ -404,27 +410,41 @@ corXY.ci <- corXY.ci %>%
left_join(trait.labs, by="Trait.comb") %>%
dplyr::select(-Test)
-"# A tibble: 52 x 11
- Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
- <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
- 1 35 0.300 0.0714 0.270 0.350 0.993 999 TRUE FALSE 1 PlantHeight
- 2 34 0.279 0.201 0.250 0.326 0.993 999 TRUE FALSE 1 LeafP
- 3 31 0.253 0.165 0.231 0.303 0.987 999 TRUE FALSE 1 SLA
- 4 37 0.244 0.0847 0.218 0.293 0.986 999 TRUE FALSE 1 Seed.length
- 5 33 0.243 0.0984 0.218 0.290 0.994 999 TRUE FALSE 1 LeafN
- 6 32 0.240 0.0846 0.218 0.290 0.986 999 TRUE FALSE 1 LeafC.perdrymass
- 7 44 0.234 0.145 0.208 0.288 0.980 999 TRUE FALSE 1 Disp.unit.leng
- 8 2 0.234 0.263 0.216 0.283 0.593 1998 FALSE FALSE 1 LEB_F_Nanophanerophyt
- 9 5 0.233 0.267 0.214 0.287 0.571 999 FALSE FALSE 1 LEB_F_Hemiphanerophyt
-10 30 0.227 0.177 0.203 0.281 0.970 999 FALSE FALSE 1 LeafArea
-# … with 42 more rows"
-
-
traits.sign.alone <- corXY.ci %>%
filter(sign_plus) %>%
pull(trait.name)
traits.sign <- traits %>%
dplyr::select(species0, any_of(traits.sign.alone))
-write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t")
+#write_delim(traits.sign, path="_data/Mesobromion/traits.out.10perc.sign.txt", delim="\t")
+## Pres Abs No Nas
+"# A tibble: 13 x 11
+ Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+ <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
+ 1 36 0.313 0.0582 0.284 0.362 0.996 999 TRUE FALSE 1 PlantHeight
+ 2 50 0.309 0.0598 0.279 0.361 0.992 999 TRUE FALSE 1 BL_ANAT
+ 3 2 0.284 0.0746 0.248 0.331 0.994 999 TRUE FALSE 1 LEB_F_Nanophanerophyt
+ 4 32 0.260 0.140 0.234 0.305 0.981 999 TRUE FALSE 1 SLA
+ 5 30 0.256 0.215 0.223 0.311 0.980 999 TRUE FALSE 1 BL_DAU
+ 6 35 0.254 0.0577 0.226 0.308 0.971 999 TRUE FALSE 1 LeafP
+ 7 49 0.253 0.209 0.222 0.303 0.983 999 TRUE FALSE 1 STRAT_T
+ 8 33 0.251 0.272 0.223 0.307 0.979 999 TRUE FALSE 1 LeafC.perdrymass
+ 9 34 0.234 0.0945 0.209 0.283 0.968 999 TRUE FALSE 1 LeafN
+10 5 0.233 0.115 0.196 0.291 0.981 999 TRUE FALSE 1 LEB_F_Hemiphanerophyt
+11 38 0.222 0.152 0.198 0.274 0.966 999 TRUE FALSE 1 Seed.length
+12 9 0.217 0.110 0.179 0.272 0.968 999 TRUE FALSE 1 LEB_F_Chamaephyt
+13 31 0.211 0.0831 0.190 0.261 0.958 999 TRUE FALSE 1 LeafArea "
+
+
+## Pres Abs - Without excluding species with NA in traits
+"# A tibble: 7 x 11
+ Trait.comb Coef.obs Coef.perm q025 q975 greater.than.perm n sign_plus sign_minus ntraits trait.name
+ <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <lgl> <lgl> <int> <fct>
+1 36 0.307 0.109 0.279 0.358 0.999 999 TRUE FALSE 1 PlantHeight
+2 32 0.260 0.122 0.235 0.308 0.987 999 TRUE FALSE 1 SLA
+3 35 0.253 0.135 0.228 0.299 0.986 999 TRUE FALSE 1 LeafP
+4 33 0.240 0.184 0.210 0.296 0.975 999 TRUE FALSE 1 LeafC.perdrymass
+5 34 0.233 0.166 0.207 0.283 0.972 999 TRUE FALSE 1 LeafN
+6 38 0.228 0.0792 0.205 0.278 0.975 999 TRUE FALSE 1 Seed.length
+7 39 0.226 0.141 0.201 0.282 0.973 999 TRUE FALSE 1 LDMC "
diff --git a/01b_MesobromionCluster.R b/01b_MesobromionCluster.R
index ae5032b..8a805d8 100644
--- a/01b_MesobromionCluster.R
+++ b/01b_MesobromionCluster.R
@@ -32,8 +32,8 @@ get.ci <- function(corXY){
greater.than.perm=mean(Coef.obs>Coef.perm),
n=n())) %>%
#calculate significance using permuted correlations
- mutate(sign_plus=greater.than.perm>0.995) %>%
- mutate(sign_minus=greater.than.perm<0.005) %>%
+ mutate(sign_plus=greater.than.perm>0.95) %>%
+ mutate(sign_minus=greater.than.perm<0.05) %>%
ungroup() %>%
mutate(Trait.comb=as.character(Trait.comb)) %>%
rowwise() %>%
@@ -71,11 +71,11 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
myfunction <- get(myfunction)
## calculate corr between species composition matrix and traits
species <- read_delim(species.path, delim="\t") %>%
- as.data.frame() %>%
+ as.data.frame() %>%
column_to_rownames("RELEVE_NR")
traits <- read_delim(traits.path, delim="\t") %>%
- as.data.frame() #%>%
-
+ as.data.frame() #%>%
+
traits <- traits %>%
rename_all(.funs=~gsub(pattern=".mean$", replacement="", x=.)) %>%
mutate_if(~is.character(.), .funs=~as.factor(.)) %>%
@@ -87,14 +87,14 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
rownames_to_column("species0") %>%
filter(complete.cases(.)) %>%
column_to_rownames("species0")
-
+
species <- species %>%
dplyr::select(rownames(traits)) %>%
mutate(sumVar = rowSums(.)) %>%
mutate_all(.funs=~./sumVar) %>%
dplyr::select(-sumVar)
}
-
+
if(combinations=="all") {
## create list of indices for each combination of traits up to a max number of interactions
@@ -121,70 +121,70 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
tt <- unlist(allcomb.t[i])
corXY.output <- rbind(corXY.output,
myfunction(comm=species, trait=traits,
- trait.sel=tt, bootstrap=nperm))
+ trait.sel=tt, bootstrap=nperm))
}
save(corXY.output, file = paste(output,"_", chunk.i, "_.RData", sep=""))
}
if(combinations=="sequential"){
-
+
if(start.round>1){
print(paste("Load data from previous round=", start.round-1))
load(file = paste(output, "_round_", start.round-1, ".RData", sep=""))
- } else {corXY.ci <- NULL}
+ } else {corXY.ci <- NULL}
for(nround in start.round:max.inter.t){
corXY.output <- NULL
## select combination of traits based on best
if(nround==1) {
allcomb.t <- lapply(1:ncol(traits), function(x){return(x)})
nall <- length(allcomb.t)
- } else {
- if(nround >= relax.round){
- best <- corXY.ci %>%
- filter(sign_plus==T) %>%
- filter(ntraits==(nround-1)) %>%
- arrange(desc(Coef.obs)) %>%
- slice(1) %>%
- pull(Trait.comb)
-
- new.best.row <- corXY.ci %>%
- filter(Trait.comb==best)
- upper <- new.best.row$q975
- lower <- new.best.row$q025
- print(paste("Assumptions relaxed - new best at round=", nround, "is trait combination", best))
- } else {
- best.row <- corXY.ci %>%
- filter(Trait.comb==best)
- upper <- best.row$q975
- lower <- best.row$q025
- }
- best.split <- as.numeric(unlist(strsplit(best, "_")))
- max.inter.ti <- nround-length(best.split)
- #n.traits <- ncol(traits)
- list.traits <- as.numeric(traits.sign.alone)
- list.traits <- list.traits[-which(list.traits %in% best.split)] ## list of remaining traits without best
- if(length(list.traits) < max.inter.ti) {
- save(corXY.output, file = paste(output, ".RData", sep=""))
- if(ncores>1) {stopCluster(cl)}
- stop("Tested all combo of significant traits")
- }
- allcomb.t <- NULL
- #for(n.inter in 1:max.inter.t){
- allcomb.t <- c(allcomb.t, combn(list.traits, max.inter.ti, simplify=F))
- #}
- #add best to all combo
- allcomb.t <- lapply(allcomb.t, function(x){return(c(best.split, x))})
- nall <- length(allcomb.t)
+ } else {
+ if(nround >= relax.round){
+ best <- corXY.ci %>%
+ filter(sign_plus==T) %>%
+ filter(ntraits==(nround-1)) %>%
+ arrange(desc(Coef.obs)) %>%
+ slice(1) %>%
+ pull(Trait.comb)
+
+ new.best.row <- corXY.ci %>%
+ filter(Trait.comb==best)
+ upper <- new.best.row$q975
+ lower <- new.best.row$q025
+ print(paste("Assumptions relaxed - new best at round=", nround, "is trait combination", best))
+ } else {
+ best.row <- corXY.ci %>%
+ filter(Trait.comb==best)
+ upper <- best.row$q975
+ lower <- best.row$q025
}
+ best.split <- as.numeric(unlist(strsplit(best, "_")))
+ max.inter.ti <- nround-length(best.split)
+ #n.traits <- ncol(traits)
+ list.traits <- as.numeric(traits.sign.alone)
+ list.traits <- list.traits[-which(list.traits %in% best.split)] ## list of remaining traits without best
+ if(length(list.traits) < max.inter.ti) {
+ save(corXY.output, file = paste(output, ".RData", sep=""))
+ if(ncores>1) {stopCluster(cl)}
+ stop("Tested all combo of significant traits")
+ }
+ allcomb.t <- NULL
+ #for(n.inter in 1:max.inter.t){
+ allcomb.t <- c(allcomb.t, combn(list.traits, max.inter.ti, simplify=F))
+ #}
+ #add best to all combo
+ allcomb.t <- lapply(allcomb.t, function(x){return(c(best.split, x))})
+ nall <- length(allcomb.t)
+ }
## Divide in chunks if requested
-# if(chunkn>1){
-# print(paste("nround", nround, "--- divide in", chunkn, "chunks and run on chunk n=", chunk.i))
-# indices <- 1:length(allcomb.t)
-# chunks <- split(indices, sort(indices%%chunkn))
-# allcomb.t <- allcomb.t[chunks[[chunk.i]] ]
-# }
+ # if(chunkn>1){
+ # print(paste("nround", nround, "--- divide in", chunkn, "chunks and run on chunk n=", chunk.i))
+ # indices <- 1:length(allcomb.t)
+ # chunks <- split(indices, sort(indices%%chunkn))
+ # allcomb.t <- allcomb.t[chunks[[chunk.i]] ]
+ # }
#names.t <- unlist(lapply(allcomb.t, paste, collapse="_"))
## start main foreach loop
@@ -196,8 +196,8 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
for(i in 1:length(allcomb.i)){
tt <- unlist(allcomb.i[i])
corXY.output <- rbind(corXY.output,
- get.corXY.bootstrap(comm=species, trait=traits,
- trait.sel=tt, bootstrap=nperm))
+ get.corXY.bootstrap(comm=species, trait=traits,
+ trait.sel=tt, bootstrap=nperm))
}
return(corXY.output)
}
@@ -225,7 +225,7 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
}
if(nround>1){
- better <- corXY.ci %>%
+ better <- corXY.ci %>%
filter(ntraits==nround) %>%
filter(q025>upper) %>%
arrange(desc(Coef.obs)) %>%
@@ -247,8 +247,8 @@ Mesobromion <- function(species.path, traits.path, output, myfunction="get.corXY
print(paste("new best at round=", nround, "is trait combination", best))
} else {print(paste("no new best at round=", nround))}
}
- print(paste("save intermediate results at round", nround))
- save(corXY.output, best, traits.sign.alone, corXY.ci, file = paste(output, "_round_",nround, ".RData", sep=""))
+ print(paste("save intermediate results at round", nround))
+ save(corXY.output, best, traits.sign.alone, corXY.ci, file = paste(output, "_round_",nround, ".RData", sep=""))
}
}
if(ncores>1){stopCluster(cl)}
diff --git a/99_HIDDEN_functions.R b/99_HIDDEN_functions.R
index 23bab26..1b6bcfe 100644
--- a/99_HIDDEN_functions.R
+++ b/99_HIDDEN_functions.R
@@ -13,9 +13,11 @@ library(abind)
library(ade4)
library(energy)
-
#### Function 1b - CorXY bootstrap####
+
+
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="_")
@@ -38,26 +40,33 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){
syn.out.tmp <- matrix.x(comm=comm, traits=traits[,ii,drop=F], scale=T)$matrix.X #, asym.bin=binary.traits
W.beals <- as.data.frame(beals(comm, include=T, type=2))
# permute traits
- index.traits <- lapply(1:(bootstrap+1), function(x){sample(1:n.species, replace=F)})
+ ### Create vector of species to resample from, which exclude species with NAs
+ ids <- 1:nrow(traits)
+ id.nas <- ids[which(rowSums(is.na(traits))>0)]
+ toresample <- ids[!ids %in% id.nas]
+ ## create list of indices of permuted traits, but without shuffling traits with missing data
+ ## the length of the list is = to the number of bootstraps
+ index.traits <- lapply(1:(bootstrap+1), function(x){
+ out <- rep(NA, nrow(traits))
+ out[id.nas] <- id.nas
+ out[is.na(out)] <- sample(toresample, replace=F)
+ return(out)
+ })
+ #index.traits <- lapply(1:(bootstrap+1), function(x){sample(1:nrow(traits), replace=F)})
syn.out.perm.tmp <- matrix.x(comm=comm, traits=traits[index.traits[[bootstrap+1]],ii,drop=F],
scale=T)$matrix.X #, asym.bin=binary.traits
-
- #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
+ (RD.tmp <- dcor(W.beals, as.data.frame(syn.out.tmp))^2)
+ (RD.perm.tmp <- dcor(W.beals, as.data.frame(syn.out.perm.tmp))^2)
corXY <- NULL
corXY <- 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]],])
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[[b]],ii,drop=F],
scale=T)$matrix.X #, asym.bin=binary.traits
- #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="RD",
@@ -69,6 +78,68 @@ get.corXY.bootstrap <- function(comm, traits, trait.sel="all", bootstrap=199){
}
+## Cracked version of matrix.x from SYNCSA to account for asymmetric binary variables in gowdis
+## and to avoid error propagation after gowdis, when species have NAs in traits
+
+matrix.x <- function (comm, traits, scale = TRUE, ranks = TRUE, ord, notification = TRUE,
+ ...)
+{
+ comm <- as.matrix(comm)
+ vartype <- var.type(traits)
+ if (any(vartype == "n")) {
+ stop("\n trait must contain only numeric, binary or ordinal variables \n")
+ }
+ if (missing(ord)) {
+ for (i in 1:length(vartype)) {
+ if (ranks & vartype[i] == "o") {
+ traits[, i] <- rank(traits[, i], na.last = "keep")
+ }
+ traits[, i] <- as.numeric(traits[, i])
+ }
+ traits <- as.matrix(traits)
+ }
+ matrix.w <- sweep(comm, 1, rowSums(comm, na.rm = TRUE), "/")
+ w.NA <- apply(matrix.w, 2, is.na)
+ matrix.w[w.NA] <- 0
+ if (notification) {
+ if (any(w.NA)) {
+ warning("Warning: NA in community data", call. = FALSE)
+ }
+ }
+ x.NA <- apply(traits, 2, is.na)
+ if (notification) {
+ if (any(x.NA)) {
+ warning("Warning: NA in traits matrix", call. = FALSE)
+ }
+ }
+ if (scale) {
+ dist.traits <- FD::gowdis(traits, ...)
+ dist.traits <- as.matrix(dist.traits)
+ dist.traits[which(rowSums(is.na(dist.traits))==(ncol(dist.traits)-1)),] <- NA ### In case a species has NO trait info, replace the diag with NA to avoid propagation
+ similar.traits <- 1 - as.matrix(dist.traits)
+ matrix.traits <- 1/colSums(similar.traits, na.rm = TRUE)
+ matrix.u <- sweep(similar.traits, 1, matrix.traits, "*")
+ }
+ else {
+ dist.traits <- as.matrix(vegan::vegdist(traits, method = "euclidean",
+ diag = TRUE, upper = TRUE, na.rm = TRUE))
+ similar.traits <- 1 - (dist.traits/max(dist.traits, na.rm = TRUE))
+ matrix.traits <- 1/colSums(similar.traits, na.rm = TRUE)
+ matrix.u <- sweep(similar.traits, 1, matrix.traits, "*")
+ }
+ u.NA <- apply(matrix.u, 2, is.na)
+ if (notification) {
+ if (any(u.NA)) {
+ warning("Warning: NA in matrix U", call. = FALSE)
+ }
+ }
+ matrix.u[u.NA] <- 0
+ matrix.X <- matrix.w %*% matrix.u
+ matrix.X <- sweep(matrix.X, 1, rowSums(matrix.X), "/") ## RE-SCALE to avoid problems propagating from NAs
+ return(list(matrix.w = matrix.w, matrix.u = matrix.u, matrix.X = matrix.X))
+}
+
+
@@ -224,62 +295,6 @@ get.SES <- function(obs.df, perm.df, stat="RV") {
-## Cracked version of matrix.x from SYNCSA to account for asymmetric binary variables in gowdis
-
-matrix.x <- function (comm, traits, scale = TRUE, ranks = TRUE, ord, notification = TRUE,
- ...)
-{
- comm <- as.matrix(comm)
- vartype <- var.type(traits)
- if (any(vartype == "n")) {
- stop("\n trait must contain only numeric, binary or ordinal variables \n")
- }
- if (missing(ord)) {
- for (i in 1:length(vartype)) {
- if (ranks & vartype[i] == "o") {
- traits[, i] <- rank(traits[, i], na.last = "keep")
- }
- traits[, i] <- as.numeric(traits[, i])
- }
- traits <- as.matrix(traits)
- }
- matrix.w <- sweep(comm, 1, rowSums(comm, na.rm = TRUE), "/")
- w.NA <- apply(matrix.w, 2, is.na)
- matrix.w[w.NA] <- 0
- if (notification) {
- if (any(w.NA)) {
- warning("Warning: NA in community data", call. = FALSE)
- }
- }
- x.NA <- apply(traits, 2, is.na)
- if (notification) {
- if (any(x.NA)) {
- warning("Warning: NA in traits matrix", call. = FALSE)
- }
- }
- if (scale) {
- dist.traits <- FD::gowdis(traits, ...)
- similar.traits <- 1 - as.matrix(dist.traits)
- matrix.traits <- 1/colSums(similar.traits, na.rm = TRUE)
- matrix.u <- sweep(similar.traits, 1, matrix.traits, "*")
- }
- else {
- dist.traits <- as.matrix(vegan::vegdist(traits, method = "euclidean",
- diag = TRUE, upper = TRUE, na.rm = TRUE))
- similar.traits <- 1 - (dist.traits/max(dist.traits, na.rm = TRUE))
- matrix.traits <- 1/colSums(similar.traits, na.rm = TRUE)
- matrix.u <- sweep(similar.traits, 1, matrix.traits, "*")
- }
- u.NA <- apply(matrix.u, 2, is.na)
- if (notification) {
- if (any(u.NA)) {
- warning("Warning: NA in matrix U", call. = FALSE)
- }
- }
- matrix.u[u.NA] <- 0
- matrix.X <- matrix.w %*% matrix.u
- return(list(matrix.w = matrix.w, matrix.u = matrix.u, matrix.X = matrix.X))
-}
diff --git a/session.R b/session.R
index c5e4d25..08d6e08 100644
--- a/session.R
+++ b/session.R
@@ -4,10 +4,10 @@ traits.path <- "_data/Mesobromion/traits.v2.10perc.txt"
output <- "_derived/Mesobromion/Cover/HIDDENtest"
myfunction <- "get.corXY.bootstrap"
max.inter.t <- 1
-chunk.i <- 1
-nperm <- 3
+chunk.i <- 46
+nperm <- 5
ncores <- 1
-chunkn <- 1
+chunkn <- 53
combinations <- "all"
start.round <- 1
relax.round <- 1
--
GitLab