# functional.betapart.core.pairwise Computes the basic quantities needed for computing the pairwise functional dissimilarity matrices. This function is similar to functional.betapart.core with multi=FALSE but it provides more options for computing convex hulls shaping each assemblage (through option passed to qhull algorithm in 'geometry::convhulln()' as well as their intersections (computed based on library 'geometry' whenever possible, else with 'rcdd' as in functional.betapart.core. ```r functional.betapart.core.pairwise(x, traits, return.details = TRUE, parallel = FALSE, opt.parallel = beta.para.control(), convhull.opt = qhull.opt(), progress = FALSE) ``` ## Arguments - `x`: A data frame, where rows are sites and columns are species. - `traits`: A data frame, where rows are species and columns are functional space dimensions (i.e. quantitative traits or synthetic axes after PCoA). Number of species in each site must be strictly higher than number of dimensions. - `return.details`: A logical value indicating if informations concerning the computation of the convexhull volumes should be returned (`TRUE`) or not (`FALSE`). See `Details`. - `parallel`: A logical value indicating if internal parallelization is used to compute pairwise dissimilarities (`TRUE`), see Examples.. - `opt.parallel`: A list of four values to modify default values used to define and run the parallel cluster. See `beta.para.control` and the `Details` section. - `convhull.opt`: A list of two named vectors of character `conv1` and `conv2` defining the options that will be passed to the `convhulln` function to compute the convexhull volumes ([http://www.qhull.org/html/qh-optq.htm](http://www.qhull.org/html/qh-optq.htm)). `conv1` sets the default option that will be passed tp `convhulln` (`'QJ'` by default instead of `'Qt'`), while `conv2` set the options if `convhulln` return an error with options set by `conv1`. See `Details` and the examples. - `progress`: A logical indicating if a progress bar should be displayed (`TRUE`) or not (by default) (`FALSE`). ## Returns The function returns an object of class `betapart` with the following elements: - **sumFRi**: The sum of the functional richness values of all sites - **FRt**: The total functional richness in the dataset `NA`. Kept for compatibility with `functional.betapart.core` - **a**: The multiple-site analog of the shared functional richness term, `NA`. Kept for compatibility with `functional.betapart.core` - **shared**: A matrix containing the functional richness shared between pairs of sites - **not.shared**: A matrix containing the functional richness not shared between pairs of sites: b, c - **sum.not.shared**: A matrix containing the total functional richness not shared between pairs of sites: b+c - **max.not.shared**: A matrix containing the total maximum functional richness not shared between pairs of sites: max(b,c) - **min.not.shared**: A matrix containing the total minimum functional richness not shared between pairs of sites: min(b,c) - **details**: `NA` if `return.details = FALSE`. Otherwise a list of two elements: - **`$FRi`**: A data frame with two columns: the `FRi` values and the qhull options used to compute them (`qhull.opt`). - **`$Intersection`**: A data frame with the pairs of communities (`Comms`), the function used to compute the volume of their intersections (`Inter`), and the qhull options used (`qhull.opt`). ## Details - **opt.parallel**: Among the four options (see `beta.para.control`), the number of cores (`nc`) and the number of convexhull volumes computed by each core at each iteration (`size`) are the most important with this function. As `inter_rcdd` is very fast, it is necessary to set a large value (>100) for `size`. Otherwise the parallelisation would not be so efficient. With a low number of communities, using internal parallelisation will slow down the function. - **convhull.opt**: Some specific distributions of points could generate errors when computing the convexhull volumes with the default qhull options ('Qt'), which is why `'QJ'` was preferred as the default value (`conv1`). Sometimes, it could be interesting to use alternative options such as 'Qt' or 'Qs' to achieve the computation. These alternative options could be used to compute all the convexhull volumes with `conv1`, or only when an error occurs using `conv2`. It is thus possible to define `conv1` as `'Qt'` to use the default 'qhull' options, but to prevent errors by setting `'QJ'` to the `conv2` argument. ## Examples ```r ##### 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" ) ## Not run: #### 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 ## End(Not run) # 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 ## Not run: # 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 ## End(Not run) # 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") ## Not run: ##### 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 ## End(Not run) ```

functional.betapart.core.pairwise function