Skip to content
Snippets Groups Projects
Commit 680c334b authored by Francesco Sabatini's avatar Francesco Sabatini
Browse files

Modified script to account for NAs in matrix.x and trait permutation

parent 133db535
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
00_Mesobromion_DataPreparation.R
+ 66
− 46
View file @ 680c334b
......@@ -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 "
01b_MesobromionCluster.R
+ 58
− 58
View file @ 680c334b
......@@ -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)}
......
99_HIDDEN_functions.R
+ 80
− 65
View file @ 680c334b
......@@ -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))
}
......
session.R
+ 3
− 3
View file @ 680c334b
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment