##### 4 communities in a 2D functional space (convex hulls are rectangles)
traits.test <- cbind(c(1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5),
c(1, 2, 4, 1, 2, 3, 5, 1, 4, 3, 5))
dimnames(traits.test) <- list(paste("sp", 1:11, sep=""), c("Trait 1", "Trait 2"))
comm.test <- matrix(0, 4, 11, dimnames = list(c("A", "B", "C", "D"),
paste("sp", 1:11, sep="")))
comm.test["A", c(1, 2, 4, 5)] <- 1
comm.test["B", c(1, 3, 8, 9)] <- 1
comm.test["C", c(6, 7, 10, 11)] <- 1
comm.test["D", c(2, 4, 7, 9)] <- 1
plot(5, 5, xlim = c(0, 6), ylim = c(0, 6), type = "n", xlab = "Trait 1", ylab = "Trait 2")
points(traits.test[, 1], traits.test[, 2], pch = 21, cex = 1.5, bg = "black")
rect(1, 1, 4, 4, col = "#458B0050", border = "#458B00")
text(2.5, 2.5, "B" , col = "#458B00", cex = 1.5)
polygon(c(2, 1, 3, 4), c(1, 2, 5, 4), col = "#DA70D650", border = "#DA70D6")
text(2.5, 3, "D", col = "#DA70D6", cex = 1.5)
rect(1, 1, 2, 2, col = "#FF000050", border = "#FF0000")
text(1.5, 1.5, "A", col = "#FF0000", cex = 1.5)
rect(3, 3, 5, 5, col = "#1E90FF50", border = "#1E90FF")
text(4, 4.2, "C", col = "#1E90FF", cex = 1.5)
# for multiple dissimilarity, multi = TRUE
test.core <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = TRUE, return.details = FALSE)
test.core
# for pairwise dissimilarity, multi = FALSE
test.core <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = FALSE)
test.core
# to use systematilcally the "QJ" options
test.core <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = FALSE,
convhull.opt = list(conv1 = "QJ"))
# to use the "QJ" options only if the convhull function generates an error
# instead of returning NA
test.core <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = FALSE,
convhull.opt = list(conv2 = "QJ"))
# using functional.betapart.core to get details on intersections
# when only pairwise dissimilarity is computed
test.core.pair <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = TRUE)
test.core.pair
# for multiple dissimilarity
test.core.multi <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = TRUE, return.details = TRUE)
test.core.multi
# using core outputs to compute pairwise and multiple functional dissimilarities
functional.beta.pair(x = test.core.pair, index.family = "jaccard" )
functional.beta.multi(x = test.core.multi, index.family = "jaccard" )
##### using internal parallelisation to fasten pairiwse dissimilarity
# by default it use serial computation
test.core.pair <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = FALSE)
# by default it uses half of the cores and 1 task per run (this can be customised)
# test.core.pairp <- functional.betapart.core(x = comm.test, traits = traits.test,
# multi = FALSE, return.details = FALSE,
# fbc.step = FALSE, parallel = TRUE)
# you can set the number of core to use :
test.core.pairp <- functional.betapart.core(x = comm.test, traits = traits.test,
multi = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = TRUE,
opt.parallel = beta.para.control(nc = 2))
all.equal(test.core.pair, test.core.pairp)
# library(microbenchmark)
# microbenchmark(serial =
# functional.betapart.core(comm.test, traits.test, multi = FALSE,
# return.details = FALSE, fbc.step = FALSE,
# parallel = FALSE),
# nc2 =
# functional.betapart.core(comm.test, traits.test, multi = FALSE,
# return.details = FALSE, fbc.step = FALSE,
# parallel = TRUE,
# opt.parallel = beta.para.control(nc = 2)),
# nc4 =
# functional.betapart.core(comm.test, traits.test, multi = FALSE,
# return.details = FALSE, fbc.step = FALSE,
# parallel = TRUE,
# opt.parallel = beta.para.control(nc = 4))
# )
if (FALSE) {
# If the number of species is very different among communities
# load-balancing parallelisation could be very efficient
# especially when the number of community is high
test.core.pairp <- functional.betapart.core(comm.test, traits.test, multi = FALSE,
return.details = FALSE, fbc.step = FALSE,
parallel = TRUE,
opt.parallel =
beta.para.control(nc = 2, LB = TRUE))
# use you can use fork cluster (but not on Windows)
#test.core.pairp <- functional.betapart.core(comm.test, traits.test, multi = FALSE,
# return.details = FALSE, fbc.step = FALSE,
# parallel = TRUE,
# opt.parallel =
# beta.para.control(nc = 2, type = "FORK"))
# finally you can customise the number of task run at each time
test.core.pairp <- functional.betapart.core(comm.test, traits.test, multi = FALSE,
return.details = FALSE, fbc.step = FALSE,
parallel = TRUE,
opt.parallel =
beta.para.control(nc = 2, size = 6))
# using internal parallelisation is not always useful, especially on small data set
# load balancing is very helpful when species richness are highly variable
# Null model using 'external' parallel computing
# Example 1: pairwise functional beta diversity (functional.beta.pair)
# Note that this is an example with a small number of samples and null model
# permutations for illustration.
# Real null model analyses should have a much greater number of samples and permutations.
##### 4 communities in a 3D functional space
comm.test <- matrix(0, 4, 11, dimnames = list(c("A", "B", "C", "D"),
paste("sp", 1:11, sep = "")))
comm.test["A", c(1, 2, 4, 5)] <- 1
comm.test["B", c(1, 3, 8, 9)] <- 1
comm.test["C", c(6, 7, 10, 11)] <- 1
comm.test["D", c( 2, 4, 7, 9)] <- 1
set.seed(1)
traits.test <- matrix(rnorm(11*3, mean = 0, sd = 1), 11, 3)
dimnames(traits.test) <- list(paste("sp", 1:11, sep = "") ,
c("Trait 1", "Trait 2", "Trait 3"))
# Required packages
library(doSNOW)
library(picante)
library(fastmatch)
library(foreach)
# define number of cores
# Use parallel::detectCores() to determine number of cores available in your machine
nc <- 2
# 4 cores would be better (nc <- 4)
# create cluster
cl <- snow::makeCluster(nc)
# register parallel backend
doSNOW::registerDoSNOW(cl)
# define number of permutations for the null model (the usual is 1000)
# make sure that nperm/nc is a whole number so that all cores have the same number
# of permutations to work on
nperm <- 100
test.score <- functional.betapart.core(comm.test, traits.test, multi = FALSE,
warning.time = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = FALSE)
obs.pair.func.dis <- functional.beta.pair(x = test.score, index.family = "sorensen")
# transform functional.beta.pair results into a matrix
obs.pair.func.dis <- do.call(rbind, obs.pair.func.dis)
# set names for each pair of site
pair_names <- combn(rownames(comm.test), 2, FUN = paste, collapse = "_")
colnames(obs.pair.func.dis) <- pair_names
# export necessary variables and functions to the cluster of cores
snow::clusterExport(cl = cl, c("comm.test", "traits.test"),
envir = environment())
# creation of an iterator to run 1 comparaisons on each core at time
it <- itertools::isplitIndices(nperm, chunkSize = 1)
# parallel computation
null.pair.func.dis <-
foreach(n = it, .combine = c, .packages=c("picante","betapart","fastmatch")) %dopar% {
# it enables to adjust the number of permutations (nt) done on each run
nt <- length(n)
null.pair.temp <- vector(mode = "list", length = nt)
# for each core "n" perform "nt" permutations
for (j in 1:nt){
# randomize community with chosen null model
# for this particular example we used the "independent swap algorithm"
# but the user can choose other types of permutation
# or create it's own null model
null.comm.test <- randomizeMatrix(comm.test, null.model = "independentswap",
iterations=1000)
# execute functional.betapart.core function
null.test.score <-
try(functional.betapart.core(null.comm.test, traits = traits.test,
multi = FALSE, warning.time = FALSE,
return.details = FALSE, fbc.step = FALSE,
parallel = FALSE), silent = TRUE)
# using 'external' parallelisation it is necessary to set parralel to FALSE
# in this artificial example there are a few combinations of species that
# the convex hull cannot be calculated due to some odd geometric combination
# so we need to re-permute the community matrix
while(inherits(null.test.score, "try-error")){
null.comm.test <- randomizeMatrix(comm.test, null.model = "independentswap",
iterations = 1000)
null.test.score <-
try(functional.betapart.core(x = null.comm.test, traits = traits.test,
multi = FALSE, warning.time = FALSE,
return.details = FALSE, fbc.step = FALSE,
parallel = FALSE), silent = TRUE)
}
# compute the pairwise beta-diversity null values and input them in the
# temporary result matrix
res <- functional.beta.pair(x = null.test.score, index.family = "sorensen")
null.pair.temp[[j]] <- do.call(rbind, res)
}
#retrieve the results from each core
null.pair.temp
}
# stop cluster
snow::stopCluster(cl)
# compute the mean, standard deviation and p-values of dissimilarity metrics
# for each pair of site
mean.null.pair.func <- matrix(numeric(), ncol = ncol(obs.pair.func.dis),
nrow = nrow(obs.pair.func.dis))
sd.null.pair.func <- matrix(numeric(), ncol = ncol(obs.pair.func.dis),
nrow = nrow(obs.pair.func.dis))
p.pair.func.dis <- matrix(numeric(), ncol = ncol(obs.pair.func.dis),
nrow = nrow(obs.pair.func.dis))
# for each one of the 3 null dissimilarity metrics (SIN, SNE and SOR)
for (j in 1:nrow(obs.pair.func.dis)){
matnull <- sapply(null.pair.func.dis, function(x) x[j,])
mean.null.pair.func[j,] <- rowMeans(matnull)
sd.null.pair.func[j,] <- sqrt(rowSums((matnull - mean.null.pair.func[j,])^2)/(nperm-1))
p.pair.func.dis[j,] <- rowSums(matnull >= obs.pair.func.dis[j,])
p.pair.func.dis[j,] <- (pmin(p.pair.func.dis[j,],nperm-p.pair.func.dis[j,])+1)/(nperm+1)
# the +1 is to take into account that the observed value is one of the possibilities
}
# compute standardized effect sizes
ses.pair.func.dis <- (obs.pair.func.dis - mean.null.pair.func)/sd.null.pair.func
# Example 2: multiple functional beta diversity (functional.beta.multi)
# Note that this is an example with a small number of samples and null model
# permutations for illustration.
# Real null model analyses should have a much greater number of samples
# and permutations.
##### 4 communities in a 3D functional space
comm.test <- matrix(0, 4, 11,dimnames = list(c("A", "B", "C", "D"),
paste("sp", 1:11, sep = "")))
comm.test["A", c(1, 2, 4, 5)] <- 1
comm.test["B", c(1, 3, 8, 9)] <- 1
comm.test["C", c(6, 7, 10, 11)] <- 1
comm.test["D", c(2, 4, 7, 9)]<-1
set.seed(1)
traits.test <- matrix(rnorm(11*3, mean=0, sd=1), 11, 3)
dimnames(traits.test) <- list(paste("sp", 1:11, sep=""),
c("Trait 1", "Trait 2", "Trait 3"))
# Required packages
library(doSNOW)
library(picante)
library(fastmatch)
library(foreach)
# define number of cores
# Use parallel::detectCores() to determine number of cores available in your machine
nc <- 2
# create cluster
cl <- snow::makeCluster(nc)
# register parallel backend
doSNOW::registerDoSNOW(cl)
# define number of permutations for the null model (the usual is 1000)
# make sure that nperm/nc is a whole number so that all cores have the same number
# of permutations to work on
nperm <- 10
# compute observed values for multiple functional dissimilarities
test.score <- functional.betapart.core(comm.test, traits.test, multi = TRUE,
warning.time = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = FALSE)
obs.multi.func.dis <- do.call(cbind, functional.beta.multi(test.score,
index.family = "sorensen"))
# export necessary variables and functions to the cluster of cores
snow::clusterExport(cl = cl, c("comm.test", "traits.test"),
envir=environment())
it <- itertools::isplitIndices(nperm, chunkSize = 1)
null.multi.func.dis <-
foreach(n = it, .combine = rbind,
.packages = c("picante","betapart","fastmatch")) %dopar% {
# for each core, create temporary matrix to store 3 null multiple functional
# dissimilarity indices (SIN, SNE,SOR)
null.multi.temp <- matrix(numeric(), ncol = 3, nrow = length(n),
dimnames = list(NULL, c("funct.beta.SIM", "funct.beta.SNE",
"funct.beta.SOR")))
# number of tasks per core (i.e., permutations per core)
nt <- length(n)
# for each core "n" perform "nt" permutations
for (j in 1:nt) {
# randomize community matrix with chosen null model (for this example
# we chose the "independent swap" algorithm)
null.comm.test <- randomizeMatrix(comm.test, null.model="independentswap",
iterations=1000)
# execute functional.betapart.core function identifying each "n" core
# with the core.ident argument for external parallelization,
null.test.score <-
try(functional.betapart.core(null.comm.test, traits.test, multi = TRUE,
warning.time = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = FALSE),
silent = TRUE)
# in this artificial example there are a few combinations of species
# that the convex hull
# cannot be calculated due to some odd geometric combination so we
# need to re-permute the community matrix
while(inherits(null.test.score, "try-error")){
null.comm.test <- randomizeMatrix(comm.test, null.model="independentswap",
iterations=1000)
null.test.score <-
try(functional.betapart.core(null.comm.test, traits.test, multi = TRUE,
warning.time = FALSE, return.details = FALSE,
fbc.step = FALSE, parallel = FALSE),
silent = TRUE)
}
# input null values in the temporary result matrix
null.multi.temp[j,] <- unlist(functional.beta.multi(null.test.score,
index.family = "sorensen"))
}
# recover results from each core
null.multi.temp
}
# close cluster
snow::stopCluster(cl)
# result matrix
result <- matrix(numeric(), ncol = 3, nrow = 3,
dimnames = list(c("obs", "ses", "p"), colnames(obs.multi.func.dis)))
# input observed values for the multiple functional dissimilarity indices (SIN, SNE,SOR)
result[1,] = obs.multi.func.dis
# compute standardized effect sizes (ses) for the multiple functional
# dissimilarity indices (SIN, SNE,SOR)
result[2,] <- (obs.multi.func.dis-colMeans(null.multi.func.dis, na.rm=TRUE))/
apply(null.multi.func.dis,2, sd, na.rm=TRUE)
# compute p-values for the multiple functional dissimilarity indices (SIN, SNE,SOR)
for (i in 1:3) {
result[3, i] <- sum(obs.multi.func.dis[i]<=null.multi.func.dis[,i])
result[3, i] <- (pmin(result[3, i], nperm - result[3, i]) + 1)/(nperm+1)
}
# the +1 is to take into account that the observed value is one of the possibilities
result
###
}
Run the code above in your browser using DataLab