functional.betapart.core.pairwise: functional.betapart.core.pairwise

Description Usage Arguments Details Value Examples

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

1
2
3
4
5
6
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

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
##### 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)

betapart documentation built on April 11, 2021, 9:06 a.m.