##### 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 pairwise dissimilarity
test.core <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test,
return.details = FALSE)
test.core
# equivalent to
test <- functional.betapart.core(x = comm.test, traits = traits.test,
return.details = FALSE,
multi = FALSE)
all.equal(test.core, test)
# using core outputs to compute pairwise and multiple functional dissimilarities
functional.beta.pair(x = test.core, index.family = "jaccard" )
if (FALSE) {
#### using convhulln options
# a data set that generates NA (due to errors) with functional.betapart.core
data(betatest)
# a list of 2 data.frame : comm.test & traits.test
comm.test <- betatest$comm.test
traits.test <- betatest$traits.test
test <- functional.betapart.core(x = comm.test, traits = traits.test,
return.details = FALSE,
multi = FALSE)
any(is.na(test$shared))
# no NA because the default option was set to QJ
# if we use the default option of qhull (Qt) :
test <- functional.betapart.core(x = comm.test, traits = traits.test,
return.details = FALSE,
multi = FALSE,
convhull.opt = list(conv1 = "Qt"))
any(is.na(test$shared))
# some NA arise
}
# with functional.betapart.core.pairwise
test.core <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test,
return.details = FALSE,
convhull.opt = list(conv1 = "Qt"))
any(is.na(test.core$shared))
# here no NA were generated because the volumes of the intersections
# were computed only with inter_geom
# to know which functions were used, set return.details to TRUE
test.core <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test,
return.details = TRUE,
convhull.opt = list(conv1 = "Qt"))
test.core$details$Intersection
#### convhull options
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
# simulating a case with species very close to each other
# here species 2, 4, 3, 8 close to species 1
# so it is like commA and commB have only 2 species instead of 4 in the 2D space
traits.test0<-traits.test
traits.test0[c(2,5),]<-traits.test0[1,]+10^-6
traits.test0[c(3,8),]<-traits.test0[1,]+10^-6
traits.test0
if (FALSE) {
# trying .core function with default qhull option
core.test0 <- functional.betapart.core(x = comm.test, traits = traits.test0, multi=FALSE,
convhull.opt = list(conv1 = "Qt"))
core.test0 # crashing because of coplanarity
# trying new .core.pairwise with default qhull option for convex hull
core.pair_test0 <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test0,
convhull.opt = list(conv1 = "Qt"))
# with default qhull options (Qt) it coud be impossible to compute the functional volumes
# this is why 'Qj' is now used as the default option for the convexhull volumes
# but other options could be passed to convexhull
# trying new core.pairwise with Qs for convex hull
core.pair_test0_Qs <- functional.betapart.core.pairwise(x = comm.test,
traits = traits.test0,
convhull.opt = list(conv1= "Qs"))
# not working
}
# trying new .pairwise with QJ (default option) for convex hull
core.pair_test0_Qj <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test0,
convhull.opt = list(conv1 = "QJ"),
return.details = TRUE)
# OK and QJ applied to each volume computation
core.pair_test0_Qj
# equivalent to
core.pair_test0_Qj <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test0,
return.details = TRUE)
# but QJ could be applied only when error are generated with other options :
core.pair_test0_Qja <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test0,
convhull.opt = list(conv1 = "Qt",
conv2 = "QJ"),
return.details = TRUE)
# OK and QJ applied only for one volume computation (community 'B')
core.pair_test0_Qja
# numerous intersection had to be computed with inter_rcdd
# the results are comparable
all.equal(core.pair_test0_Qj[-9], core.pair_test0_Qja[-9]) # -9 to remove details
# the pairwise functional functional betadiversity
functional.beta.pair(core.pair_test0_Qj, index.family = "jaccard")
if (FALSE) {
##### using internal parallelisation to fasten pairiwse dissimilarity
# by default (parallel = FALSE) the code is run in serial
test.core.pair <- functional.betapart.core.pairwise(x = comm.test, traits = traits.test,
parallel = FALSE)
# when using internal parallelisation and default options
# 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.pairwise(x = comm.test, traits = traits.test,
parallel = TRUE,
opt.parallel = beta.para.control(nc = 2))
all.equal(test.core.pair, test.core.pairp)
# as inter_geom is much more faster than inter_rccd it is useful to increase the number
# of calculus run per iteration to limit the time consumed by the cluster itself
# you can play on size (this would not have sense here as there is only few computation)
test.core.pairp <-
functional.betapart.core.pairwise(x = comm.test, traits = traits.test,
parallel = TRUE,
opt.parallel =
beta.para.control(nc = 2,
size = 100))
all.equal(test.core.pair, test.core.pairp)
library(microbenchmark)
microbenchmark(
serial = functional.betapart.core.pairwise(comm.test, traits.test),
nc2 = functional.betapart.core.pairwise(comm.test, traits.test,
parallel = TRUE,
opt.parallel = beta.para.control(nc = 2)),
nc4 = functional.betapart.core.pairwise(comm.test, traits.test, multi = FALSE,
parallel = TRUE,
opt.parallel = beta.para.control(nc = 4))
)
# If the number of species is very different among communities
# load-balancing parallelisation could be more efficient
# especially when the number of community is high
test.core.pairp <-
functional.betapart.core.pairwise(comm.test, traits.test,
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.pairwise(comm.test, traits.test,
parallel = TRUE,
opt.parallel =
beta.para.control(nc = 2, type = "FORK"))
# a progress bar can be displayed to asses the evolution of the computations
test.core.pairp <-
functional.betapart.core.pairwise(comm.test, traits.test,
parallel = TRUE,
opt.parallel =
beta.para.control(nc = 2, LB = TRUE),
progress = TRUE)
# 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(foreach)
library(itertools)
# 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.pairwise(comm.test, traits.test)
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
# 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)
# 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 = 10)
null.pair.func.dis <-
foreach(n = it, .combine = c,
.packages=c("picante","betapart","fastmatch", "rcdd", "geometry")) %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 <-
functional.betapart.core.pairwise(null.comm.test, traits = traits.test,
parallel = FALSE)
# 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 specify to use the 'QJ' options in case of errors
# 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
}
Run the code above in your browser using DataLab