R: functional.betapart.core.pairwise

functional.betapart.core.pairwise {betapart}

R Documentation

functional.betapart.core.pairwise

Description

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.

Usage

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). `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`).

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 distribution of points could generate errors when computing the convexhull volumes with the default qhull options ('Qt'), that is why 'QJ' was preferred as default values (conv1). Sometimes, it could be interesting to use alternative options such as 'Qt' or 'Qs' or to use alternative options 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 error by setting 'QJ' to the conv2 argument.

Value

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)

Examples

##### 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)

[Package betapart version 1.6 Index]