Package 'adespatial'

Construct asymmetric eigenvector maps (AEM)

Description

This function constructs eigenvectors of a site-by-link matrix. Weights can be applied to the links.

Usage

aem(
  aem.build.binary,
  binary.mat,
  weight,
  rm.link0 = FALSE,
  print.binary.mat = FALSE
)

Arguments

aem.build.binary	Object created by function aem.build.binary.
binary.mat	Site (n rows) by link (k columns) matrix. The 1s in the matrix represents the presence of a link influencing a site, directly or indirectly, otherwise the values are 0s.
weight	Vector of weights of length k, to be applied to the links.
rm.link0	Logical (TRUE, FALSE) determining if the links directly connecting a site to the origin (site 0) should be removed. Default value: FALSE. This parameter is only used when an object of class aem.build.binary is provided to the function.
print.binary.mat	Logical (TRUE, FALSE) determining if the site-by-link matrix used in the analysis should be printed.

Details

If only an object of class aem.build.binary is given to this function, The argument binary.mat is not considered. binary.mat is only considered when the argument aem.build.binary is missing.

If weights are applied to the links, the length of vector weight has to take into account wether the links connecting real sites to the origin (the fictitious site 0) have been kept or removed.

Value

value	A vector of singular values associated with the AEM.
vectors	A matrix of eigenvector. Each column is an AEM eigenfunction (or variable).
mod.binary.mat	A site-by-link matrix modified through the function.

Note

It sometimes happens that AEM eigenfunctions have equal singular values. In that case, different sets of AEM eigenfunctions may be produced on different plateforms.

Eigenvectors associated to an eigenvalue that is smaller than $10^{-12}$ are considered negligeable. They have been removed from the created AEM eigenfunctions.

Author(s)

F. Guillaume Blanchet

References

Blanchet F.G., P. Legendre and Borcard D. (2008) Modelling directional spatial processes in ecological data. Ecological Modelling, 215, 325-336.

See Also

aem.build.binary, svd

Examples

# Construction of object of class nb (spdep)
if(require("spdep", quietly = TRUE)){
nb <- cell2nb(5,5,"queen")

# Create fictitious geographical coordinates 
xy <- cbind(1:25,expand.grid(1:5,1:5))

# Build binary site-by-link matrix
bin.mat <- aem.build.binary(nb,xy)

# Construct AEM eigenfunctions from an object of class aem.build.binary
res <- aem(aem.build.binary=bin.mat,rm.link0=FALSE)
res$values

# Illustrate 4 AEM eigenfunctions using bubble plots
opal <- palette()
palette(c("black","white"))
oldpar <- par(mfrow=c(2,2))
symbols(x=xy[,2:3], circles=abs(res$vectors[,1]), inches=FALSE, asp=1,
 fg=ifelse(sign(-res$vectors[,1])+1>0,1,0), 
 bg=ifelse(sign(res$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
title("AEM 1")
symbols(x=xy[,2:3], circles=abs(res$vectors[,2]), inches=FALSE, 
asp=1, fg=ifelse(sign(-res$vectors[,2])+1>0,1,0),
 bg=ifelse(sign(res$vectors[,2])+1>0,1,0), xlab="x", ylab="y")
title("AEM 2")
symbols(x=xy[,2:3], circles=abs(res$vectors[,3]), inches=FALSE, 
asp=1, fg=ifelse(sign(-res$vectors[,3])+1>0,1,0), 
bg=ifelse(sign(res$vectors[,3])+1>0,1,0), xlab="x", ylab="y")
title("AEM 3")
symbols(x=xy[,2:3], circles=abs(res$vectors[,4]), inches=FALSE, asp=1,
 fg=ifelse(sign(-res$vectors[,4])+1>0,1,0), 
 bg=ifelse(sign(res$vectors[,4])+1>0,1,0), xlab="x", ylab="y")
title("AEM 4")

# Construct AEM eigenfunctions using only a site-by-link matrix
res2 <- aem(binary.mat=bin.mat[[1]])
res2$values

# Illustrate 4 AEM eigenfunctions using bubble plots
par(mfrow=c(2,2))
symbols(x=xy[,2:3], circles=abs(res2$vectors[,1]), inches=FALSE, 
asp=1, fg=ifelse(sign(-res2$vectors[,1])+1>0,1,0), 
bg=ifelse(sign(res2$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
title("AEM 1")
symbols(x=xy[,2:3], circles=abs(res2$vectors[,2]), inches=FALSE,
asp=1, fg=ifelse(sign(-res2$vectors[,2])+1>0,1,0), 
bg=ifelse(sign(res2$vectors[,2])+1>0,1,0), xlab="x", ylab="y")
title("AEM 2")
symbols(x=xy[,2:3], circles=abs(res2$vectors[,3]), inches=FALSE,
asp=1, fg=ifelse(sign(-res2$vectors[,3])+1>0,1,0), 
bg=ifelse(sign(res2$vectors[,3])+1>0,1,0), xlab="x", ylab="y")
title("AEM 3")
symbols(x=xy[,2:3], circles=abs(res2$vectors[,4]), inches=FALSE,asp=1,
 fg=ifelse(sign(-res2$vectors[,4])+1>0,1,0), 
 bg=ifelse(sign(res2$vectors[,4])+1>0,1,0), xlab="x", ylab="y")
title("AEM 4")

palette(opal)
par(oldpar)

# Construct AEM eigenfunctions with a function of the distance
# as weights to put on the links

# Construction of object of class nb (spdep)
nb<-cell2nb(5,5,"queen")

# Create fictitious geographical coordinates
xy <- cbind(1:25,expand.grid(1:5,1:5))

# Build binary site-by-link matrix
bin.mat <- aem.build.binary(nb,xy)

# Construct a matrix of distances
long.lien.mat<-as.matrix(dist(xy))

# Extract the edges, remove the ones directly linked to site 0
lien.b<-bin.mat$edges[-1:-5,]

# Construct a vector giving the length of each edge
long.lien<-vector(length=nrow(lien.b))

for(i in 1:nrow(lien.b)){
	long.lien[i]<-long.lien.mat[lien.b[i,1],lien.b[i,2]]
}

# Construct a vector of weights based on distance
weight.vec<-1-(long.lien/max(long.lien))^2

# Construct AEM eigenfunctions from an object of class aem.build.binary
res <- aem(aem.build.binary=bin.mat,weight=weight.vec,rm.link0=TRUE)
res

# Computing Moran's I for AEMs

# Building AEMs
xy <- cbind(1:25,expand.grid(1:5,1:5))
Wdist <- 1/as.matrix(dist(xy[,2:3]))

nb <- cell2nb(5,5,"queen")
bin.mat <- aem.build.binary(nb,xy)
linkBase <- bin.mat[[2]]
link <- linkBase[-which(linkBase[,1] == 0),]
weight <- numeric()

for(i in 1:nrow(link)){
   weight[i] <- Wdist[link[i,1],link[i,2]]
}

AEM <- aem(bin.mat, weight = weight, rm.link0 = TRUE)

# Constructing asymmetric matrix
matasym <- matrix(0,ncol=25, nrow=25)

for(i in 1:nrow(link)){
    matasym[link[i,1],link[i,2]]<- weight[i]
}

# Build a listw object from the asymmetric matrix
listwAsym <-  mat2listw(matasym, style = "B", zero.policy = TRUE)

# Calculate Moran's I for AEM
MoranIAEM <- moran.randtest(AEM$vectors, listwAsym)

}

---

Construct a site-by-edge binary matrix

Description

This function construct a site-by-edge binary matrix. It uses a set of sites coordinates and a connexion diagram (object of class nb from the spdep package). The 1s in the matrix represents the presence of a link influencing a site, directly or indirectly, otherwise the values are 0s. Graphically, the function is implemented such that the directional process is considered to be going from the bottom to the top of the screen in the graphical output of R. As such, the origin is underneath the set of points representing the sites. Prepare the table of site coordinates accordingly.

Usage

aem.build.binary(
  nb.object = NULL,
  coords,
  link = NULL,
  unit.angle = "degrees",
  rot.angle = 0,
  rm.same.y = TRUE,
  plot.connexions = TRUE
)

Arguments

nb.object	Object of class nb from library spdep.
coords	A three columns matrix or data frame. Colums 1: identifiers of the points (needs to be numeric). Columns 2 and 3: the X and Y coordinates of the points.
link	A two columns matrix. Each row define an edge. Column 1: The site from which the edge starts. Column 2: the site to which the edge ends. All values in link need to be integers.
unit.angle	Character. The measurement units in which the angle is defined: either "degrees" (default) or "radians".
rot.angle	Numeric. Angle of the vector describing the process influencing the sites. This argument generate a rotation of the site coordinates. The set of coordinates is rotated counterclockwise. Negative values will produce a clockwise rotation.
rm.same.y	Logical (TRUE, FALSE). Determines if the links perpendicular to the gradient should be removed. Default value: TRUE. If these links have already been removed this argument put to TRUE will make the function crash. See detail for more information.
plot.connexions	Logical (TRUE, FALSE). Determines if the sites and the associated connexion diagram should be plotted after rotation of the coordinates by gradient.angle.

Details

The lowest site in the gradient is the one that will connect to the fictitious site 0 to consider direction. Note that if there are multiple lowest sites, they will all be connected to the fictitious site 0 to consider direction.

The site-by-edge matrix created by this function and the list of edges include the links associated to a fictitious site upstream of all other, see Blanchet et al. (2008), for details. The decision regarding wether the origin and the edges associated with it should be kept or removed is left to the user. Removal of this site and of its associated edges can be done manually after the construction of the site-by-edge matrix and of the list edges. It can also be done when running the function aem.

If the connexion diagram was modified so that the links connecting sites that are exactly perpendicular to the gradient have been removed or if there is no sites exactly at the same level in the gradient, defining rm.same.y to TRUE will generate an error.

If all the sites have the same y coordinates after rotation, e.g. a horizontal transect perpendicular to the defined spatial asymmetry, this analysis should not be used.

The argument plot.connexions will plot the sites (coords) in black, after rotation, if any, and the connexion diagram (nb.object), in red. The site labels are also plotted on the graph. To show the direction of the spatial asymmetry considered by the function, a fictive site (in blue) was added upstream. This fictive site is linked (blue edges) to the site(s) that are the most upstream ones. Since this graph is generic, it might sometimes look odd, however, the information given will remain the accurate.

Value

se.mat	A binary (n x k) matrix of site (n rows) by link edges (k columns).
edges	A matrix describing the link edges. It has 2 columns (from, to) and as many rows as there are edges. The edges linked to the fictitious site of origin are found at the beginning of the list.

Author(s)

F. Guillaume Blanchet

References

Blanchet F.G., P. Legendre and Borcard D. (2008) Modelling directional spatial processes in ecological data. Ecological Modelling, 215, 325-336.

See Also

aem

Examples

### Create an object of class nb (spdep)
if(require("spdep", quietly = TRUE)){
nb<-cell2nb(5,5,"queen")

### Create fictitious geographical coordinates 
xy <- cbind(1:25,expand.grid(1:5,1:5))

### Build a binary site-by-link matrix; remove the site which have identical Y coordinate
### (by default argument: rm.same.y = TRUE)
bin.mat <- aem.build.binary(nb,xy)
str(bin.mat)

### Build a binary site-by-link matrix using the argument link: remove the site which
### have identical Y coordinate (by default argument: rm.same.y = TRUE)
edges<-expand.grid(1,2:25)
bin.mat <- aem.build.binary(coords=xy,link=edges)
str(bin.mat)

### Build a binary site-by-link matrix, making the process affect the points at 
### an angle of 45 degrees
bin.mat.45 <- aem.build.binary(nb,xy, rot.angle=45)
str(bin.mat.45)

### Build a binary site-by-link matrix, making the process affect the points at
### an angle of pi/3 radians
bin.mat.pi3 <- aem.build.binary(nb,xy,unit.angle="radians", rot.angle=pi/3)
str(bin.mat.pi3)
}

---

AEM for time series

Description

This function constructs AEM eigenfunctions for multi-scale analysis of a regular time series or spatial transect of univariate or multivariate data.

Usage

aem.time(n, w = NULL, moran = FALSE)

Arguments

n	Numeric. Number of points in the series.
w	A vector of weights to be applied to the edges (columns of matrix E). Equal weights are used if no vector w is provided. The length of vector w must be (n-1) where n is the number of points in the spatial or temporal series.
moran	Logical. If TRUE, Moran's I are computed for all AEM. If FALSE (default value), Moran's I are not computed.

Details

Time series represent a form of directional stochastic process. To emphasize the directional nature of the process influencing the data, AEM analysis, which was designed to take trends into account, should be applied to the non-detrended series. MEM analysis (see scores.listw) can be applied to data series that were detrended to remove the directional component as recommended by Blanchet et al. (2008, 2011) and Legendre & Legendre (2012, Subsection 14.1.2). Detrended palaeoecological sediment core data, for example, could be studied by MEM analysis.

No data file needs to be provided to this function. The AEM eigenvectors are constructed from a matrix E generated from the regular sequence of points along the series.

A vector of weights w can be provided, representing the ease of communication of matter, energy or information among the points. The most simple form would be the inverse of (d/dmax) where d is the distance between adjacent nodes and dmax is the maximum distance between adjacent nodes in the spatial or time series. More general forms of weights may represent the inverse of landscape resistance to the movement of organisms, propagules, genes, etc.

If the calculation of Moran's I is requested, the point coordinates are generated from the point positions along the series.

Value

E	Nodes-by-edges matrix E.
values	Eigenvalues of the principal component analysis of E.
aem	Matrix of AEM eigenfunctions normalized to unit length.
Moran	Moran's I statistics tested by a bilateral test with 999 permutations
listw	An object of class listw with the associated spatial weighting matrix

Author(s)

Pierre Legendre and F. Guillaume Blanchet

References

Blanchet F.G., P. Legendre and Borcard D. (2008) Modelling directional spatial processes in ecological data. Ecological Modelling, 215, 325-336.

Blanchet F.G., P. Legendre, R. Maranger, D. Monti, and P. Pepin. (2011) Modelling the effect of directional spatial ecological processes at different scales. Oecologia, 166, 357-368.

Legendre, P. and L. Legendre (2012) Numerical Ecology, 3rd English edition. Elsevier Science BV, Amsterdam.

Legendre, P. and O. Gauthier (2014) Statistical methods for temporal and space-time analysis of community composition data. Proceedings of the Royal Society B - Biological Sciences, 281, 20132728.

See Also

aem, scores.listw

Examples

# Time series containing 20 equispaced observations
out <- aem.time(20, moran = TRUE)

# Time series containing 20 observations with unequal spacing
# Generate (n-1) random interpoint distances
distances <- runif(19,1,5)

# Compute weights representing the ease of communication among points
w <- 1/(distances/max(distances))

# Compute the AEM eigenfunctions
out <- aem.time(20, w = w, moran = TRUE)

---

Weight edges when constructing AEM variables

Description

These functions construct a vector of weights that can be associated to the edges of the connexion diagram used as a basis to build AEM eigenfunctions. aem.weight.edges is general and can be used for 1 or 2 dimensional problems. aem.weight.time is meant to be used only for time series. It is a wrapper for aem.weight.edges.

Usage

aem.weight.edges(
  nb.object,
  coords,
  distmat = NULL,
  alpha = 2,
  beta = NULL,
  max.d = NULL,
  unit.angle = "degrees",
  rot.angle = 0,
  rm.same.y = TRUE,
  plot.connexions = TRUE
)

aem.weight.time(
  dates,
  distmat = NULL,
  alpha = 2,
  beta = NULL,
  max.d = NULL,
  unit.angle = "degrees",
  rot.angle = 0,
  rm.same.y = TRUE,
  plot.connexions = TRUE
)

Arguments

nb.object	Object with class 'nb', computed by the spdep package, containing a list of neighbours for each sampling unit (site or time).
coords	A three-column matrix or data frame. Column 1: identifiers of the points (must be numeric). Columns 2 and 3: the X and Y coordinates of the points.
distmat	Class 'matrix' or 'dist' object containing a dissimilarity or distance matrix. (See details).
alpha	Numeric. Exponent of the first weighting function. (See details).
beta	Numeric. Exponent of the second weighting function. (See details).
max.d	Numeric. Maximum distance for weighting. Default value if max.d=NULL: the maximum distance among a set of sites divided by 2 or the full span of a time series divided by 2 (not recommended in most problems, see details). A warning is given if max.d = NULL and the default value is used.
unit.angle	Character. The measurement units in which the angle is defined: either "degrees" (default) or "radians".
rot.angle	Numeric. Angle of the vector describing the process influencing the sites. This argument generates a rotation of the site coordinates. The set of coordinates is rotated counterclockwise. Negative values will produce a clockwise rotation.
rm.same.y	Logical (TRUE, FALSE). Determines if the links perpendicular to the gradient should be removed. Default value: TRUE. If these links have already been removed, this argument put to TRUE will make the function crash. See details for more information.
plot.connexions	Logical (TRUE, FALSE). Determines if the sites and the associated connexion diagram should be plotted after rotation of the coordinates by gradient.angle.
dates	A vector of dates, class 'numeric' or 'Date'.

Details

These functions should be used in close relationship with aem.build.binary, consequently many of the arguments in this function and in aem.build.binary are the same.

The argument distmat may contain general forms of dissimilarity, for example the difficulty of transfer of individuals, matter or energy among the sampling units through space or time.

In aem.weight.edges, two weighting functions, described in Legendre and Legendre (2012, eqs. 114.3 and 14.4) have been implemented, where $d_{ij}$ is the distance between sites $i$ and $j$:

Weighting function 1:	$1 - (d_{ij}/max(d))^\alpha$
Weighting function 2:	$1/d_{ij}^\beta$

Also note that if a value is provided for beta (that is, if it is not NULL), weighting function 2 is used regardless of whether alpha is defined or not.

In most applications, the default value of max.d is not optimal. A more meaningful solution in many applications is to compute a Moran's I correlogram (for univariate data) or a Mantel correlogram (for multivariate data), and provide the distance where the correlation becomes 0 as the value for max.d.

Value

A vector of weights associating a value to each edge of the graph.

Functions

• aem.weight.time():

Author(s)

Olivier Gauthier, Pierre Legendre and F. Guillaume Blanchet

References

Legendre, P. and L. Legendre (2012) Numerical Ecology, 3rd English edition. Elsevier Science BV, Amsterdam.

Legendre, P. and O. Gauthier (2014) Statistical methods for temporal and space-time analysis of community composition data. Proceedings of the Royal Society B - Biological Sciences, 281, 20132728.

See Also

aem.build.binary, sp.correlogram, mantel.correlog

Examples

### Time serie example
### Example - 12 dates (days from January 1st of year 1) 
### in a 6-year study starting September 5, 2000
if(require("spdep", quietly = TRUE)){
dates <- as.Date(c(129,269,500,631,864,976,1228,1352,1606,1730,1957,2087),origin="2000/1/1")
autocor.limit <- 522  # Limit of autcorrelation in the correlogram

### Using aem.weight.time()
(wtime <- aem.weight.time(dates, alpha=2, max.d=autocor.limit))
### Using aem.weight.edges()
n <- length(dates)
nb <- cell2nb(1, n)
xy.dates <- cbind(1:n, rep(1, n), dates)
(wtime <- aem.weight.edges(nb, xy.dates, alpha=2, max.d=autocor.limit))

n <- length(dates)
nb <- cell2nb(1, n)
xy.dates <- cbind(1:n, dates, rep(1, n)) ## Note the inversion of 'dates' and 'rep(1,n)'
wtime <- aem.weight.edges(nb, xy.dates, alpha=2, 
max.d=autocor.limit,rot.angle=90) # Note that 'rot.angle=90' was used

### Spatial example using default d.max (notice the warning)
###########################################################################
nb<-cell2nb(5,5,"queen")
xy <- cbind(1:25,expand.grid(1:5,1:5))
(wspace <- aem.weight.edges(nb,xy))
}

---

Bacterial production data set

Description

Longitude and latitude of 25 samples when bacterial production was sampled in Lake St. Pierre (Québec, Canada). These sampled were carried out in August 18, 2005.

Usage

data(bacProdxy)

Format

A data frame with the coordinates of the 25 sampled locations

References

Blanchet F.G., P. Legendre, R. Maranger, D. Monti, and P. Pepin. (2011) Modelling the effect of directional spatial ecological processes at different scales. Oecologia, 166, 357-368.

---

Beta diversity computed as Var(Y)

Description

Compute estimates of total beta diversity as the total variance in a community data matrix Y, as well as derived SCBD and LCBD statistics, for 19 dissimilarity coefficients or the raw data table. Computing beta diversity as Var(Y) for raw, untransformed community composition data is not recommended. Tests of significance of the LCBD indices are also produced.

Usage

beta.div(
  Y,
  method = "hellinger",
  sqrt.D = FALSE,
  samp = TRUE,
  nperm = 999,
  adj = TRUE,
  save.D = FALSE,
  clock = FALSE
)

Arguments

Y	Community composition data. The object class can be either data.frame or matrix.
method	One of the 19 dissimilarity coefficients available in the function: "hellinger", "chord", "log.chord", "chisquare", "profiles", "percentdiff", "ruzicka", "divergence", "canberra", "whittaker", "wishart", "kulczynski", "jaccard", "sorensen", "ochiai", "ab.jaccard", "ab.sorensen", "ab.ochiai", "ab.simpson", "euclidean". See Details. Names can be abbreviated to a non-ambiguous set of first letters. Default: method="hellinger".
sqrt.D	If sqrt.D=TRUE, the dissimilarities in matrix D are square-rooted before computation of SStotal, BDtotal and LCBD. This transformation may be useful for methods "manhattan", "whittaker", "divergence", "canberra", "percentdiff", "ruzicka", "wishart" since square-root transformation of the dissimilarities makes these D matrices Euclidean. Note 1 - Euclideanarity is useful for ordination by principal coordinate analysis; lack of this property does not adversely affect SStotal, BDtotal and LCBD. Note 2 - The logical value given to parameter sqrt.D has no incidence on calculations through methods "euclidean", "profiles", "hellinger", "log.chord", "chord", "chisquare" since no D matrix is computed in those cases. Note 3 - For methods "jaccard", "sorensen", "ochiai", that function produces the dissimilarity matrix in the form sqrt(D), which is Euclidean.
samp	If samp=TRUE, the abundance-based distances (ab.jaccard, ab.sorensen, ab.ochiai, ab.simpson) are computed for sample data. If samp=FALSE, they are computed for true population data.
nperm	Number of permutations for the tests of significance of LCBD indices.
adj	Compute adjusted p-values using the Holm method. Default: adj=TRUE
save.D	If save.D=TRUE, the distance matrix will appear in the output list.
clock	If clock=TRUE, the computation time is printed. Useful when nperm is large.

Details

Calculations may be carried out in two ways, depending on the selected method.

• For untransformed or transformed raw data, the total sum of squares (SStotal) is first computed, then the total beta diversity (BDtotal), which is SStotal divided by (n - 1), is calculated. This algorithm is used for methods "euclidean", "profiles", "hellinger", "chord", "log.chord", "chisquare". No transformation of the data is computed when the method is "euclidean". For methods "profiles", "hellinger", "chord", "log.chord", "chisquare", the algorithm begins with computation of the same-name transformation of the community data (Legendre and Gallagher 2001; Legendre and Legendre 2012, Section 7.7; Legendre and Borcard 2018); SStotal and BDtotal are then computed for the transformed data, followed by calculation of the SCBD and LCBD indices.

• Calculations of BDtotal can also be conducted from a dissimilarity matrix. SStotal is computed by summing the squared dissimilarities in the lower triangular dissimilarity matrix and dividing by n. Then, total beta diversity (BDtotal) is obtained by dividing SStotal by (n-1). With option sqrt.D = TRUE, the computation of SStotal is equivalent to summing the distances instead of the squared distances. Choices are: "whittaker", "divergence", "canberra", "percentdiff", "ruzicka", "wishart", "kulczynski", "ab.jaccard", "ab.sorensen", "ab.ochiai", "ab.simpson", "jaccard", "sorensen", "ochiai". Equations for these dissimilarities are presented in Table 1 of Legendre and De Cáceres (2013). The Ružička index is described in Legendre (2014); this coefficient is suitable for beta diversity studies. See Chao et al. (2006) for details about the abundance-based (ab) coefficients.

Community composition data can be log-transformed prior to analysis with the chord distance; see Legendre and Borcard (2018). The log(y+1) transformation (log1p function of base) reduces the asymmetry of the species distributions. The chord-log distance, readily available among the methods of the beta.div function, is the chord distance computed on log(y+1)-transformed data. This combined transformation is meaningful for community composition data because the log is one of the transformations in the Box-Cox series, corresponding to exponent 0; see Legendre and Legendre (2012, Section 1.5.6). Exponent 1 (no transformation of the data) followed by the chord transformation and calculation of the Euclidean distance would simply produce the chord distance. Exponent 0.5 (square root) followed by the chord transformation and the Euclidean distance would produce the Hellinger distance. The chord, Hellinger and log-chord distances represent a series where the data are increasingly transformed to reduce the asymmetry of the distributions. Note that it is meaningless to subject log-transformed community compostion data to the "profiles", "hellinger", or "chisquare" distances available in this function.

The Jaccard, Sørensen and Ochiai coefficients are the binary forms of 10 of the 12 dissimilarity coefficients (including the Ružička index) that are suitable for beta diversity assessment. The equivalences are described in Legendre and De Cáceres (2013, Table 1). These popular coefficients can be computed directly using function beta.div without going to the trouble of applying the quantitative forms of these coefficients to data reduced to presence-absence form. beta.div produces the dissimilarity matrix in the form sqrt(D), which is Euclidean. Hence for these three coefficients, function beta.div should be used with option sqrt.D=FALSE.

Species contributions to beta diversity (SCBD indices for the species) are computed for untransformed or transformed raw data, but they cannot be computed from dissimilarity matrices.

Local contributions to beta diversity (LCBD indices) represent the degree of uniqueness of the sites in terms of their species compositions. They can be computed in all cases: raw (not recommended) or transformed data, as well as dissimilarity matrices. See Legendre and De Cáceres (2013) for details. LCBD indices are tested for significance by random, independent permutations within the columns of Y. This permutation method tests H0 that the species are distributed at random, independently of one another, among the sites, while preserving the species abundance distributions in the observed data. See Legendre and De Cáceres (2013) for discussion.

This version of beta.div calls computer code written in C to speed up computation, especially for the permutation tests of the LCBD indices.

Value

A list containing the following results:

• beta: Total sum of squares and total beta diversity [= Var(Y)] of the data matrix. BDtotal statistics computed with the same D index are comparable among data sets having the same or different numbers of sampling units (n), provided that they are of the same size or represent the same sampling effort.

• SCBD: Vector of Species contributions to beta diversity (SCBD), if computed.

• LCBD: Vector of Local contributions to beta diversity (LCBD) for the sites.

• p.LCBD: P-values associated with the LCBD indices.

• p.adj: Corrected P-values for the LCBD indices, Holm correction.

• method: Method selected.

• note: Notes indicate whether the selected coefficient is Euclidean or not.

• D: The distance matrix if save.D=TRUE.

When all sites contain a different set of species with no species in common, the maximum value that BDtotal can take depends on the method used in the calculation.

• With methods "hellinger", "chord", "profiles", which have maximum values of sqrt(2), BDtotal produces an index in the range [0, 1] with a maximum value of 1.

• For dissimilarity indices with maximum values of 1, BDtotal has a maximum value of 0.5.

• Dissimilarity indices that do not have maximum values of 1 or sqrt(2) produce BDtotal values that do not have an upper bound; hence they cannot be compared across taxonomic groups or among study sites. This group includes the chi-square distance.

See Legendre & De Caceres (2013, p. 957-958), Table 2 and section Maximum value of BD.
For two sites only, the LCBD results are not interesting. With all coefficients, the two LCBD indices are equal to 0.5. The two associated p-values are 1 because LCBD is 0.5 for all columnwise permutations of the data.
The calculation is aborted when Y only contains two identical rows of data. In that case, SStotal and BDtotal are 0 and the LCBD indices cannot be computed (value NaN).

Author(s)

Pierre Legendre [email protected]

References

Chao, A., R. L. Chazdon, R. K. Colwell and T. J. Shen. 2006. Abundance-based similarity indices and their estimation when there are unseen species in samples. Biometrics 62: 361-371.

Legendre, P. 2014. Interpreting the replacement and richness difference components of beta diversity. Global Ecology and Biogeography 23: 1324-1334.

Legendre, P. and D. Borcard. 2018. Box-Cox-chord transformations for community composition data prior to beta diversity analysis. Ecography 41: 1820-1824.

Legendre, P. and M. De Cáceres. 2013. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters 16: 951-963.

Legendre, P. and E. D. Gallagher, E.D. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129: 271-280.

Legendre, P. and Legendre, L. 2012. Numerical Ecology. 3rd English edition. Elsevier Science BV, Amsterdam.

Examples

if(require("vegan", quietly = TRUE) & require("adegraphics", quietly = TRUE)){
data(mite)
res = beta.div(mite, "hellinger", nperm=999)

# Plot a map of the LCBD indices using the Cartesian coordinates
data(mite.xy)
s.value(mite.xy, res$LCBD, symbol = "circle", col = c("white", "brown"), main="Map of mite LCBD")

### Example using the mite abundance data and the percentage difference dissimilarity
res = beta.div(mite, "percentdiff", nperm=999, clock=TRUE)

# Plot a map of the LCBD indices
signif = which(res$p.LCBD <= 0.05)	# Which are the significant LCBD indices?
nonsignif = which(res$p.LCBD > 0.05)	# Which are the non-significant LCBD indices?
g1 <- s.value(mite.xy[signif,], res$LCBD[signif], ppoint.alpha = 0.5, plegend.drawKey = FALSE,
 symbol = "circle", col = c("white", "red"), main="Map of mite LCBD (red = significant indices)")
g2 <- s.value(mite.xy[nonsignif,], res$LCBD[nonsignif], ppoint.alpha = 0.5,
symbol = "circle", col = c("white", "blue"))
g2+g1
}

---

Decompose D in replacement and richness difference components

Description

Podani-family and Baselga-family decompositions of the Jaccard and Sørensen dissimilarity coefficients and their quantitative forms (Ruzicka and percentage difference) into replacement and richness difference components, for species presence-absence or abundance data, as described in Legendre (2014).

Usage

beta.div.comp(mat, coef = "J", quant = FALSE, save.abc = FALSE)

Arguments

mat	Community composition data (data.frame or matrix).
coef	Family of coefficients to be computed. "S" or "Sorensen": Podani family, Sørensen-based indices. "J" or "Jaccard": Podani family, Jaccard-based indices. "BS": Baselga family, Sørensen-based indices. "BJ": Baselga family, Jaccard-based indices. "N": Podani & Schmera (2011) relativized nestedness index. The quantitative form of the Sørensen dissimilarity is the percentage difference index. The quantitative form of the Jaccard dissimilarity is the Ruzicka index.
quant	If TRUE, compute the quantitative forms of replacement, nestedness and D. If FALSE, compute the presence-absence forms of the coefficients.
save.abc	If TRUE, save the matrices of parameters a, b and c used in presence-absence calculations.

Details

For species presence-absence data, the dissimilarity coefficients are Jaccard = (b+c)/(a+b+c) and Sørensen = (b+c)/(2*a+b+c) with the usual a,b,c notation. For species abundance data, the dissimilarity coefficients are the Ruzicka index = (B+C)/(A+B+C) and Odum's percentage difference = (B+C)/(2A+B+C) (aka Bray-Curtis in some packages), where

• A = sum of the intersections (or minima) of species abundances at two sites,

• B = sum of abundances at site 1 minus A,

• C = sum of abundances at site 2 minus A.

The binary (quant=FALSE) and quantitative (quant=TRUE) forms of the S and J indices return the same values when computed for presence-absence data.

Value

A list containing the following results:

• repl: Replacement matrix, class = dist.

• rich: Richness/abundance difference or nestedness matrix (class dist). With options "BJ", "BS" and "N", rich contains nestedness indices. With option "N", the repl[i,j] and rich[i,j] values do not add up to D[i,j].

• D: Dissimilarity matrix (classdist).

• part: Beta diversity partitioning vector:

  1. BDtotal (total beta diversity) = sum(D.ij)/(n*(n-1)) (Legendre & De Cáceres 2013). This is equal to sum(d.ij^2)/(n*(n-1)) where d.ij = sqrt(D.ij). The dissimilarities are square-rooted because the Jaccard, Sørensen, Ruzicka and percentage difference indices are not Euclidean.

  2. Repl = Total replacement diversity.

  3. RichDiff|Nes = Total richness difference diversity (or nestedness).

  4. Repl/BDtotal = Total replacement diversity/Total beta diversity.

  5. RichDiff/BDtotal = Total richness difference diversity (or nestedness)/Total beta diversity.

• note: Name of the dissimilarity coefficient.

The Jaccard and Sørensen dissimilarity coefficients and their quantitative forms, the Ruzicka and percentage difference indices, all have upper bounds (Dmax) of 1. Hence, when all sites contain a different set of species with no species in common, the maximum value that BDtotal can take is 0.5. See Legendre & De Caceres (2013, p. 958), section Maximum value of BD. This differs form the values produced by function beta.div(): with methods "hellinger", "chord" and "profiles", which have maximum values of sqrt(2), BDtotal has a maximum value of 1 for these dissimilarities.

Author(s)

Pierre Legendre [email protected]

References

Baselga, A. (2010) Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography, 19, 134–143.

Baselga, A. (2012) The relationship between species replacement, dissimilarity derived from nestedness, and nestedness. Global Ecology and Biogeography, 21, 1223-1232.

Baselga, A. (2013) Separating the two components of abundance-based dissimilarity: balanced changes in abundance vs. abundance gradients. Methods in Ecology and Evolution, 4, 552-557.

Carvalho, J.C., Cardoso, P., Borges, P.A.V., Schmera, D. & Podani, J. (2013) Measuring fractions of beta diversity and their relationships to nestedness: a theoretical and empirical comparison of novel approaches. Oikos, 122, 825-834.

Legendre, P. 2014. Interpreting the replacement and richness difference components of beta diversity. Global Ecology and Biogeography, 23, 1324-1334.

Legendre, P. and M. De Cáceres. 2013. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters 16: 951-963.

Podani, J., Ricotta, C. & Schmera, D. (2013) A general framework for analyzing beta diversity, nestedness and related community-level phenomena based on abundance data. Ecological Complexity, 15, 52-61.

Podani, J. & Schmera, D. 2011. A new conceptual and methodological framework for exploring and explaining pattern in presence-absence data. Oikos, 120, 1625-1638.

See Also

LCBD.comp

Examples

if(require(ade4, quietly = TRUE)){
data(doubs)
fish.sp = doubs$fish[-8,]   # Fish data; site 8 is removed because no fish were caught

# Compute and partition a matrix of Jaccard indices (presence-absence data)
out1 = beta.div.comp(fish.sp, coef="J", quant=FALSE)
out1$part

# Compute and partition a matrix of percentage difference indices
# (quantitative form of Sorensen index)
out2 = beta.div.comp(fish.sp, coef="S", quant=TRUE)
out2$part
# In paragraph Value, see the description of the 5 elements of vector part. 
# Is the fish beta diversity dominated by replacement or richness/abundance difference?
}

---

Function to choose a connection network

Description

The function chooseCN is a simple interface to build a connection network (CN) from xy coordinates. The user chooses from 6 types of graph and one additional weighting scheme. chooseCN calls functions from appropriate packages, handles non-unique coordinates and returns a connection network either with classe nb or listw. For graph types 1-4, duplicated locations are not accepted and will issue an error.

Usage

chooseCN(
  xy,
  ask = TRUE,
  type = NULL,
  result.type = "nb",
  d1 = NULL,
  d2 = NULL,
  k = NULL,
  a = NULL,
  dmin = NULL,
  plot.nb = TRUE,
  edit.nb = FALSE
)

Arguments

xy	an matrix or data.frame with two columns for x and y coordinates.
ask	a logical stating whether graph should be chosen interactively (TRUE,default) or not (FALSE). Set to FALSE if type is provided.
type	an integer giving the type of graph (see details).
result.type	a character giving the class of the returned object. Either "nb" (default) or "listw", both from spdep package. See details.
d1	the minimum distance between any two neighbours. Used if type=5.
d2	the maximum distance between any two neighbours. Used if type=5. Can also be a character: "dmin" for the minimum distance so that each site has at least one connection, or "dmax" to have all sites connected (despite the later has no sense).
k	the number of neighbours per point. Used if type=6.
a	the exponent of the inverse distance matrix. Used if type=7.
dmin	the minimum distance between any two distinct points. Used to avoid infinite spatial proximities (defined as the inversed spatial distances). Used if type=7.
plot.nb	a logical stating whether the resulting graph should be plotted (TRUE, default) or not (FALSE).
edit.nb	a logical stating whether the resulting graph should be edited manually for corrections (TRUE) or not (FALSE, default).

Details

There are 7 kinds of graphs proposed:
Delaunay triangulation (type 1)
Gabriel graph (type 2)
Relative neighbours (type 3)
Minimum spanning tree (type 4)
Neighbourhood by distance (type 5)
K nearests neighbours (type 6)
Inverse distances (type 7)

The last option (type=7) is not a true neighbouring graph: all sites are neighbours, but the spatial weights are directly proportional to the inversed spatial distances.
Also not that in this case, the output of the function is always a listw object, even if nb was requested.

The choice of the connection network has been discuted on the adegenet forum. Please search the archives from adegenet website (section 'contact') using 'graph' as keyword.

Value

Returns a connection network having the class nb or listw. The xy coordinates are passed as attribute to the created object.

Author(s)

Thibaut Jombart [email protected]

Examples

if(require("ade4", quietly = TRUE)){
data(mafragh)

oldpar <- par(mfrow=c(2,2))
cn1 <- chooseCN(mafragh$xy,ask=FALSE,type=1)
cn2 <- chooseCN(mafragh$xy,ask=FALSE,type=2)
cn3 <- chooseCN(mafragh$xy,ask=FALSE,type=3)
cn4 <- chooseCN(mafragh$xy,ask=FALSE,type=4)
par(oldpar)
}

---

Space- And Time-Constrained Clustering

Description

Function constr.hclust carries out space-constrained or time-constrained agglomerative clustering from a multivariate dissimilarity matrix.

Usage

constr.hclust(
  d,
  method = "ward.D2",
  links,
  coords,
  beta = -0.25,
  chron = FALSE,
  members = NULL
)

Arguments

d	A dist-class dissimilarity (distance) matrix.
method	The agglomeration method to be used (default: "ward.D2"; see details).
links	A list of edges (or links) connecting the points. May be omitted in some cases; see details and examples
coords	Coordinates of the observations (data rows) in the dissimilarity matrix d. The coordinates are used for plotting maps of the clustering results. This matrix may be omitted when the user does not wish to print maps of the clustering results or when no links file is provided. coords is a matrix or data frame with two columns, following the convention of the Cartesian plane: first column for abscissa, second column for ordinate. See examples.
beta	The beta parameter for beta-flexible clustering (default: beta = -0.25).
chron	Logical (TRUE or FALSE) indicating whether a chronological (i.e. time-constrained or spatial transect) clustering should be calculated (default: chron = FALSE).
members	NULL or a vector with length size of d (default: NULL; See details).

Details

The agglomeration method to be used should be (an unambiguous abbreviation of) one of "ward.D", "ward.D2", "single", "complete", "average" (UPGMA), "mcquitty" (WPGMA), "centroid" (UPGMC), "median" (WPGMC), or "flexible". Method "ward.D2" (default) implements the Ward (1963) clustering criterion, method "ward.D" does not (Murtagh and Legendre, 2014).

Agglomerative clustering can be carried out with a constraint of spatial or temporal contiguity. This means that only the objects that are linked in links are considered to be candidates for clustering: the next pair of objects to cluster will be the pair that has the lowest dissimilarity value among the pairs that are linked.

The same rule applies during the subsequent clustering steps, which involve groups of objects: the list of links is updated after each agglomeration step. All objects that are neighbours of one of the components that have fused are now neighbours of the newly formed cluster.

The edges (links) are specified using argument links, which can be an object of class nb (see, e.g., tri2nb), an object of class listw (see, e.g., nb2listw), a two-element list or an object coercible as a such (e.g., a two-column dataframe), or a two-column matrix with each row representing an edge and the columns representing the two ends of the edges. For lists with more than two elements, as well as dataframes or matrices with more than two-columns, only the first two elements or columns are used for the analysis. The edges are interpreted as being non directional; there is no need to specify an edge going from point a to point b and one going from point b to point a. While doing so is generally inconsequential for the analysis, it carries some penalty in terms of computation time. It is a good practice to place the nodes in increasing order of numbers from the top to the bottom and from the left to the right of the list but this is not mandatory. A word of caution: in cases where clusters with identical minimum distances occur, the order of the edges in the list may have an influence on the result. Alternative results would be statistically equivalent.

When argument link is omitted, regular (unconstrained) clustering is performed and a hclust-class object is returned unless argument chron = TRUE. When argument chron = TRUE, chronological clustering is performed, taking the order of observations as their positions in the sequence. Argument links is not used when chron = TRUE. Argument chron allows one to perform a chronological clustering in the case where observations are ordered chronologically. Here, the term "chronologically" should not be taken restrictively: the method remains applicable to other sequential data sets such as spatial series made of observations along a transect.

When the graph described by link is not entirely connected, a warning message is issued to warn the user about the presence and number of disjoint clusters and a procedure is suggested to identify the disjoint clusters. The disjoint clusters (or singletons) are merged in the order of their indices (i.e. the two clusters with smallest indices are merged first) and so on until all of disjoint clusters have been merged. The dissimilarity at which these clusters are merged is a missing value (NA) in vector height (i.e., unconnected clusters have undefined dissimilarities in constrained clustering).

If members != NULL, then d is taken to be a dissimilarity matrix between clusters instead of dissimilarities between individual objects. Then, members must be a vector giving the number of observations per cluster. In this way, the hierarchical clustering algorithm can be ‘started in the middle of the dendrogram’, e.g., in order to reconstruct the part of the tree above a cut. See examples in hclust for details on that functionality."

Memory storage and time to compute constrained clustering for N objects. The Lance and Williams algorithm for agglomerative clustering uses dissimilarity matrices. The amount of memory needed to store the dissimilarities among N observations as 64-bit double precision floating point variables (IEEE 754) is 8*N*(N-1)/2 bytes. For example, a dissimilarity matrix among 22 500 observations would require 2 024 910 000 bytes (1.89 GiB) of storage whereas one among 100 000 observations would take up 39 999 600 000 bytes (37.25 GiB). The implementation in this function needs to cache a copy of the dissimilarity matrix as its elements are modified following each merging of the closest clusters or singletons, thereby doubling the amounts of required memory shown above. Memory needed to store the other information associated with the clustering is much smaller. Users should make sure to have the necessary memory space (and system stability) before attempting to analyze large data sets. What is considered a large amount of memory has increased over time as computer hardware evolved with time. We let users apply contemporary common sense as to what sample sizes represent manageable clustering problems. Computation time grows with N at roughly the same speed as memory storage requirement to store the dissimilarity matrices increases. See the Benchmarking example below.

With large data sets, a manageable output describing the classification of the sites is obtained with function cutree(x, k) where k is the number of groups. A dendrogram would be unreadable.

Value

A constr.hclust-class object.

Author(s)

Pierre Legendre [email protected] (preliminary version coded in R) and Guillaume Guénard [email protected] (present version mostly coded in C)

References

Guénard, G. and P. Legendre. 2022. Hierarchical clustering with contiguity constraint in R. Journal of Statistical Software 103(7): 1-26 doi:10.18637/jss.v103.i07

Langfelder, P. and S. Horvath. 2012. Fast R functions for robust correlations and hierarchical clustering. Journal of Statistical Software 46(11): 1-17. doi:10.18637/jss.v046.i11

Legendre, P. and L. Legendre. 2012. Numerical ecology, 3rd English edition. Elsevier Science BV, Amsterdam. doi:10.1016/S0304-3800(00)00291-X

Murtagh, F. and P. Legendre. 2014. Ward’s hierarchical agglomerative clustering method: which algorithms implement Ward’s criterion? Journal of Classification 31: 274-295. doi:10.1007/s00357-014-9161-z

Ward, J. H. 1963. Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association 58: 236-244. doi:10.1080/01621459.1963.10500845

See Also

plot.constr.hclust, hclust, cutree, and ScotchWhiskey

Examples

## First example: Artificial map data from Legendre & Legendre
##                (2012, Fig. 13.26): n = 16

dat <- c(41,42,25,38,50,30,41,43,43,41,30,50,38,25,42,41)
coord.dat <- matrix(c(1,3,5,7,2,4,6,8,1,3,5,7,2,4,6,8,
                      4.4,4.4,4.4,4.4,3.3,3.3,3.3,3.3,
                      2.2,2.2,2.2,2.2,1.1,1.1,1.1,1.1),16,2)

## Obtaining a list of neighbours:
library(spdep)
listW <- nb2listw(tri2nb(coord.dat), style="B")
links.mat.dat <- listw2mat(listW)
neighbors <- listw2sn(listW)[,1:2]

## Calculating the (Euclidean) distance between points:
D.dat <- dist(dat)

## Display the points:
plot(coord.dat, type='n',asp=1)
title("Delaunay triangulation")
text(coord.dat, labels=as.character(as.matrix(dat)), pos=3)
for(i in 1:nrow(neighbors))
    lines(rbind(coord.dat[neighbors[i,1],],
          coord.dat[neighbors[i,2],]))

## Unconstrained clustring by hclust:
grpWD2_hclust <- hclust(D.dat, method="ward.D2")
plot(grpWD2_hclust, hang=-1)

## Clustering without a contiguity constraint;
## the result is represented as a dendrogram:
grpWD2_constr_hclust <- constr.hclust(D.dat, method="ward.D2")
plot(grpWD2_constr_hclust, hang=-1)

## Clustering with a contiguity constraint described by a list of
## links:
grpWD2cst_constr_hclust <-
    constr.hclust(
        D.dat, method="ward.D2",
        neighbors, coord.dat)

## To visualize using hclust's plotting method:
## stats:::plot.hclust(grpWD2cst_constr_hclust, hang=-1)

## Plot the results on a map with k=3 clusters:
plot(grpWD2cst_constr_hclust, k=3, links=TRUE, las=1, xlab="Eastings",
     ylab="Northings", cex=3, lwd=3)

## Generic functions from hclust can be used, for instance to obtain
## a list of members of each cluster:
cutree(grpWD2cst_constr_hclust, k=3)

## Now with k=5 clusters:
plot(grpWD2cst_constr_hclust, k=5, links=TRUE, las=1, xlab="Eastings",
     ylab="Northings", cex=3, lwd=3)
cutree(grpWD2cst_constr_hclust, k=5)

## End of the artificial map example


## Second example: Scotch Whiskey distilleries clustered using tasting
## scores (nose, body, palate, finish, and the four distances combined)
## constrained with respect to the distillery locations in Scotland.

## Documentation file about the Scotch Whiskey data: ?ScotchWhiskey

data(ScotchWhiskey)

## Cluster analyses for the nose, body, palate, and finish D
## matrices:

grpWD2cst_ScotchWhiskey <-
    lapply(
        ScotchWhiskey$dist,    ## A list of distance matrices
        constr.hclust,         ## The function called by function lapply
        links=ScotchWhiskey$neighbors@data,         ## The list of links
        coords=ScotchWhiskey$geo@coords/1000
    )

## The four D matrices (nose, body, palate, finish), represented as
## vectors in the ScotchWiskey data file, are combined as follows to
## produce a single distance matrix integrating all four types of
## tastes:

Dmat <- ScotchWhiskey$dist
ScotchWhiskey[["norm"]] <-
    sqrt(Dmat$nose^2 + Dmat$body^2 + Dmat$palate^2 + Dmat$finish^2)

## This example shows how to apply const.clust to a single D matrix when
## the data file contains several matrices.

grpWD2cst_ScotchWhiskey[["norm"]] <-
    constr.hclust(
        d=ScotchWhiskey[["norm"]],method="ward.D2",
        ScotchWhiskey$neighbors@data,
        coords=ScotchWhiskey$geo@coords/1000
    )

## A fonction to plot the Whiskey clustering results:

plotWhiskey <- function(wh, k) {
   par(fig=c(0,1,0,1))
   plot(grpWD2cst_ScotchWhiskey[[wh]], k=k, links=TRUE, las=1,
        xlab="Eastings (km)", ylab="Northings (km)", cex=0.1, lwd=3,
        main=sprintf("Feature: %s",wh))
   text(ScotchWhiskey$geo@coords/1000,labels=1:length(ScotchWhiskey$geo))
   legend(x=375, y=700, lty=1L, lwd=3, col=rainbow(1.2*k)[1L:k],
          legend=sprintf("Group %d",1:k), cex=1.25)
   SpeyZoom <- list(xlim=c(314.7,342.2), ylim=c(834.3,860.0))
   rect(xleft=SpeyZoom$xlim[1L], ybottom=SpeyZoom$ylim[1L],col="#E6E6E680",
        xright=SpeyZoom$xlim[2L], ytop=SpeyZoom$ylim[2L], lwd=2, lty=1L)
   par(fig=c(0.01,0.50,0.46,0.99), new=TRUE)
   plot(grpWD2cst_ScotchWhiskey[[wh]], xlim=SpeyZoom$xlim,
        ylim=SpeyZoom$ylim, k=k, links=TRUE, las=1, xlab="", ylab="",
        cex=0.1, lwd=3, axes=FALSE)
   text(ScotchWhiskey$geo@coords/1000,labels=1:length(ScotchWhiskey$geo))
   rect(xleft=SpeyZoom$xlim[1L], ybottom=SpeyZoom$ylim[1L],
        xright=SpeyZoom$xlim[2L], ytop=SpeyZoom$ylim[2L], lwd=2, lty=1L)
}

## Plot the clustering results on the map of Scotland for 5 groups.
## The inset map shows the Speyside distilleries in detail:
plotWhiskey("nose", 5L)
plotWhiskey("body", 5L)
plotWhiskey("palate", 5L)
plotWhiskey("finish", 5L)
plotWhiskey("norm", 5L)

## End of the Scotch Whiskey tasting data example

## Not run: 

## Third example: Fish community composition along the Doubs River,
## France. The sequence is analyzed as a case of chronological
## clustering, substituting space for time.

if(require("ade4", quietly = TRUE)){
data(doubs, package="ade4")
Doubs.D <- dist.ldc(doubs$fish, method="hellinger")
grpWD2cst_fish <- constr.hclust(Doubs.D, method="ward.D2", chron=TRUE,
                                coords=as.matrix(doubs$xy))
plot(grpWD2cst_fish, k=5, las=1, xlab="Eastings (km)",
     ylab="Northings (km)", cex=3, lwd=3)

## Repeat the plot with other values of k (number of groups)

## End of the Doubs River fish assemblages example

## Example with 6 connected points, shown in Fig. 2 of Guénard & Legendre paper 

var = c(1.5, 0.2, 5.1, 3.0, 2.1, 1.4)
ex.Y = data.frame(var)

## Site coordinates, matrix xy
x.coo = c(-1, -2, -0.5, 0.5, 2, 1)
y.coo = c(-2, -1, 0, 0, 1, 2)
ex.xy = data.frame(x.coo, y.coo)

## Matrix of connecting edges E
from = c(1,1,2,3,4,3,4)
to = c(2,3,3,4,5,6,6)
ex.E = data.frame(from, to)

## Carry out constrained clustering analysis
test.out <-
    constr.hclust(
        dist(ex.Y),       # Response dissimilarity matrix
        method="ward.D2", # Clustering method
        links=ex.E,       # File of link edges (constraint) E
        coords=ex.xy      # File of geographic coordinates
    )

par(mfrow=c(1,2))
## Plot the map of the results for k = 3
plot(test.out, k=3)
## Plot the dendrogram
stats:::plot.hclust(test.out, hang=-1)
}

## Same example modified: disjoint clusters
## Same ex.Y and ex.xy as in the previous example
var = c(1.5, 0.2, 5.1, 3.0, 2.1, 1.4)
ex.Y = data.frame(var)

## Site coordinates, matrix xy
x.coo = c(-1, -2, -0.5, 0.5, 2, 1)
y.coo = c(-2, -1, 0, 0, 1, 2)
ex.xy = data.frame(x.coo, y.coo)

## Matrix of connecting edges E2
from = c(1,1,2,4,4)
to = c(2,3,3,5,6)
ex.E2 = data.frame(from, to)

## Carry out constrained clustering analysis
test.out2 <-
    constr.hclust(
        dist(ex.Y),       # Response dissimilarity matrix
        method="ward.D2", # Clustering method
        links=ex.E2,      # File of link edges (constraint) E
        coords=ex.xy      # File of geographic coordinates
    )
cutree(test.out2, k=2)

par(mfrow=c(1,2))
## Plot the map of the results for k = 3
plot(test.out2, k=3)
## Plot the dendrogram showing the disconnected groups
stats:::plot.hclust(test.out2, hang=-1)
axis(2,at=0:ceiling(max(test.out2$height,na.rm=TRUE)))

## End of the disjoint clusters example


## End(Not run)
## End of examples

---

Class For Constrained Hiereachical Clustering

Description

Files belonging to this class hold information about the constrained agglomerative clustering and allows one to display results graphically.

Format

A file belonging to this class is a list with elements:

merge

A (n-1) by 2 matrix. Row i of file "merge" describes the merging of clusters at step i of the clustering. If an element j in the row is negative, it means that observation -j was merged at this stage. If j is positive, it means that the merge was with the cluster formed at the (earlier) stage j of the algorithm. Thus negative entries in file "merge" indicate agglomerations of singletons, and positive entries indicate agglomerations of non-singletons.

height

A set of (n-1) non-decreasing real values. The clustering height is the value of the criterion associated with the clustering method for the particular agglomeration.

order

A vector giving the permutation of the original observations suitable for plotting, in the sense that a cluster plot using this ordering and matrix merge will not have crossing branches.

labels

Labels for the clustered objects.

method

The agglomerative clustering method that has been used.

call

The call that produced the result.

dist.method

The distance that has been used to create dissimilarity matrix "d" (only returned if the dissimilarity matrix object has a "method" attribute attached to it).

links

A copy of the list of edges (if a matrix of edges was provided to the function).

coords

A copy of the coordinates (if coordinates were provided to the function).

Details

The class inherits from hclust-class and describes the tree produced by the constrained clustering procedure.

All class members except links and coords are identical to those in hclust-class. several methods designed to process these objects are expected to also work with constr.hclust-class objects.

See Also

hclust-class

---

Contingency periodogram

Description

Function to compute a contingency periodogram for a univariate series of qualitative data

Usage

Cperiodogram(x, T1 = 2, T2 = NULL, nperm = NULL, alpha = 0.05, graph = TRUE)

Arguments

x	a qualitative variable (factor)
T1	first period included in the calculations (default: T1 = 2)
T2	last period included in the calculations (default: T2 = n/2)
nperm	Number of permutations for the chi-square test. For chi-square tests using the chi- square distribution, use the default nperm=NULL
alpha	significance level for computation of the confidence limits
graph	a logical indicating if a graph is requested, by default TRUE.

Details

The contingency periodogram of Legendre et al. (1981) identifies periodic components in qualitative data vectors. The vector may contain classes of a qualitative variable or the classes obtained by hierarchical clustering or partitioning of a multivariate data table. The method is also described in Legendre & Legendre (2012). The optional graph produced by the function shows the following information:

• In red: the B statistics (information in common).

• In blue: Confidence limits for B without correction.

• In green: Bonferroni-corrected limits of the confidence intervals.

• In black: Confidence limits with progressive Bonferroni correction.

Value

A table with the statistics for the selected periods:

• Wilks' chi-square statistic (Wilks.chisq)

• information in common (B),

• degrees of freedom (df),

• p-value (prob)

Confidence interval limits:

• critical value of B without correction (B.crit),

• critical value of B with Bonferroni correction based on the number of periods studied in the periodogram (B.crit.Bonf),

• critical value of B with progressive Bonferroni correction (B.prog.Bonf).

Author(s)

Pierre Legendre [email protected]

References

Legendre, L., M. Fréchette & P. Legendre. 1981. The contingency periodogram: a method of identifying rhythms in series on nonmetric ecological data. Journal of Ecology 69: 965-979.

Legendre, P. and Legendre, L. 2012. Numerical Ecology. 3rd English ed. Elsevier, Amsterdam

Examples

# Data from the numerical example of Subsection 12.4.2 of Legendre and Legendre (2012).
test.vec <- c(1,1,2,3,3,2,1,2,3,2,1,1,2,3,3,1)
# Periodogram with tests using the chi-square distribution
res <- Cperiodogram(test.vec)
# Periodogram with permutation tests
res <- Cperiodogram(test.vec, nperm=2000, graph=FALSE)

---

Combine dbMEM matrices corresponding to groups of sites

Description

This function reads a file containing the Cartesian coordinates of sites forming different groups on the map, and constructs a combined staggered matrix of dbMEM spatial eigenvectors, ready for use in RDA. The method was first described and used in Declerck et al. (2011) and summarized in the Borcard et al. (2011) book, section 7.4.3.5. These publications provided preliminary versions of the present function. The present version is more completely documented. Furthermore, it uses the dbmem function of the adespatial package for computation of the eigenfunctions.

Usage

create.dbMEM.model(coord = NULL, D.mat = NULL, nsites)

Arguments

coord	Optional file containing the Cartesian coordinates of the sites.
D.mat	Optional distance matrix provided by user, class matrix or dist. If D.mat=NULL, the geographic distance matrix will be computed from the coordinates provided in file coord.
nsites	A vector containing the number of sites per group.

Details

The geographic positions of the sites are provided either in a file of geographic coordinates coord or as a geographic distance matrix D.mat.

The sites must, of course, be in the same order in file coord (or in file D.mat) and in the response data file used in the RDA. All sites of a group must be together in these two files, i.e. not interspersed. The numbers of sites in the groups are provided in vector nsites. See example.

File vector coord, if provided, must contain Cartesian coordinates of the sites, not coordinates in degrees. The Euclidean distance computed from the geographic coordinates is a meaningful representation of the geographic relationships only if the coordinates are Cartesian. Geodetic Cartesian coordinates can be derived from Lat-Lon data in degrees using the function geoXY of the SoDA package. Beware of UTM coordinates if the sites are not all located in the same UTM zone; UTM coordinates are Cartesian only within an UTM zone. See https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system.

Value

A matrix with n rows containing a set of k staggered matrices of dbMEM eigenfunctions in its diagonal portion; n is the total number of sites in the study and k is the number of groups. Each small matrix contains the dbMEM functions, modelling positive spatial correlation, describing the spatial relationships among the sites of a group. The remainder of the matrix is filled with zeros. Zero is the mean value of all eigenfunctions describing within-group relationships. This means that during the calculation of RDA, the sites of a focus group will have, with each other, relationships described by the dbMEM eigenfunctions of that group, whereas the sites outside that group will have weights of 0 in the regressions that concern these eigenfunctions.

Author(s)

Pierre Legendre [email protected], 2010. Adaptation to adespatial: Daniel Borcard and Pierre Legendre, 2016

References

Borcard, D., F. Gillet and P. Legendre. 2011. Numerical ecology with R. Use R! series, Springer Science, New York.

Declerck, S. A. J., J. S. Coronel, P. Legendre & L. Brendonck. 2011. Scale dependency of processes structuring metacommunities of cladocerans in temporary pools of High-Andes wetlands. Ecography 34: 296-305.

See Also

dbmem

Examples

{
 # Generate random coordinates for 35 sites forming 6 distinct groups on the map
 Easting <- runif(35)+c(rep(0,6),rep(1.5,7),rep(3,6), rep(0,5),rep(1.5,5),rep(3,6))
 Northing<- runif(35)+c(rep(2.8,6),rep(2.3,7),rep(2.8,6), rep(0,5),rep(0.5,5),rep(0,6))
 cartesian <- cbind(Easting,Northing)
 rownames(cartesian) <- paste("S",1:nrow(cartesian),sep='')
 nsites.per.group <- c(6,7,6,5,5,6)

 result <- create.dbMEM.model(coord=cartesian, nsites=nsites.per.group)

 # Draw a map to check the coding of the sites into the groups
 site.codes <- unlist(apply(cbind(1:6),1,n=nsites.per.group,function(a,n) rep(a,n[a])))

 col.vec <- c("green3","gray99","orange2","gold1","brown3","gray70")
 plot(cartesian, pch=22, col="black", bg=col.vec[site.codes], cex=2, ylim=c(0,4),asp=1)
 text(cartesian,labels=rownames(cartesian), cex=0.5, pos=3)

 # Examine the staggered matrix of dbMEM eigenfunctions
 # Not run:
 result
}

---

dbMEM spatial eigenfunctions

Description

Compute distance-based Moran's eigenvector maps (dbMEM, also called dbMEM spatial eigenfunctions) from a geographic distance matrix, in view of spatial eigenfunction analysis.

Usage

dbmem(
  xyORdist,
  thresh = NULL,
  MEM.autocor = c("positive", "non-null", "all", "negative"),
  store.listw = TRUE,
  silent = TRUE
)

Arguments

xyORdist	Either a matrix of spatial coordinates or a distance matrix (class dist).
thresh	A threshold value for truncation of the geographic distance matrix. If thresh=NULL, the length of the longest edge of the minimum spanning tree will be used as the threshold (as returned by the function give.thresh).
MEM.autocor	A string indicating if all MEMs must be returned or only those corresponding to non-null, positive or negative autocorrelation. The difference between options all and non-null is the following: when there are several null eigenvalues, option all removes only one of the eigenvectors with null eigenvalues and returns (n-1) eigenvectors, whereas non-null does not return any of the eigenvectors with null eigenvalues. Default: MEM.autocor="positive".
store.listw	A logical indicating if the spatial weighting matrix should be stored in the attribute listw of the returned object
silent	A logical indicating if some information should be printed during computation: truncation level and time to compute the dbmem

Details

dbMEM eigenfunctions were called PCNM in early papers (Borcard and Legendre 2002, Borcard et al. 2004). There is a small difference in the computation: to construct PCNMs, the distance matrix subjected to PCoA contained zeros on the diagonal. In dbMEM, the matrix contains 4*thresh values on the diagonal. The result is that the dbMEM eigenvalues are smaller than the PCNM eigenvalues by a constant (equal to (n.sites * (4*thresh)^2)/2). The dbMEM eigenvalues are proportional to Moran's I coefficient of spatial correlation (Dray et al. 2006; Legendre and Legendre 2012). The dbMEM eigenvectors only differ from the PCNM eigenvectors by a multiplicative constant; this has no impact on the use of MEMs as explanatory variables in linear models. In this implementation, dbMEM eigenvectors have a norm equal to 1 (using the uniform weigts 1/n.sites).

If a truncation value is not provided, the largest distance in a minimum spanning tree linking all sites on the map is computed (returned by the function give.thresh). That value is used as the truncation threshold value (thresh).

A square regular grid produces multiple eigenvalues (i.e. eigenvalues that are equal) and multiple eigenvalues have an infinity of eigenvector solutions. Hence, different eigenvectors may be produced by this function on computers with different operating systems or implementations of R. In addition, the eigenvectors found by the dbmem function from the site coordinates may differ from the eigenvectors computed from the geographic distance matrix among the sites. Nonetheless, the different complete sets of eigenvectors will have the exact same explanatory power (R-square) for a given response vector or matrix, despite the fact that they are not fully correlated on a one-to-one basis. This is, however, not the case for subsets of eigenvectors selected using stepwise procedures.

Value

An object of class orthobasisSp , subclass orthobasis. The dbMEM eigenfunctions (principal coordinates of the truncated distance matrix) are stored as a data.frame. It contains several attributes (see ?attributes) including:

• values: The dbMEM eigenvalues.

• listw: The associated spatial weighting matrix (if store.listw = TRUE).

Author(s)

Stéphane Dray [email protected], Pierre Legendre, Daniel Borcard and F. Guillaume Blanchet

References

Borcard, D. and P. Legendre. 2002. All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecological Modelling 153: 51-68.

Borcard, D., P. Legendre, C. Avois-Jacquet and H. Tuomisto. 2004. Dissecting the spatial structure of ecological data at multiple scales. Ecology 85: 1826-1832.

Dray, S., P. Legendre and P. R. Peres-Neto. 2006. Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecological Modelling 196: 483-493.

Legendre, P. and L. Legendre. 2012. Numerical ecology, 3rd English edition. Elsevier Science BV, Amsterdam.

See Also

give.thresh, mem

Examples

if(require("ade4", quietly = TRUE) & require("adegraphics", quietly = TRUE)){

data(oribatid)
mite <- oribatid$fau      # 70 peat cores, 35 species
mite.xy <- oribatid$xy    # Geographic coordinates of the 70 cores

# Example 1: Compute the MEMs corresponding to all non-null eigenvalues
# thresh=1.012 is the value used in Borcard and Legendre (2002)
mite.dbmem1 <- dbmem(mite.xy, thresh=1.012, MEM.autocor = "non-null", silent = FALSE)
mite.dbmem1

# Print the (n-1) non-null eigenvalues
attributes(mite.dbmem1)$values
# or:  attr(mite.dbmem1, "values")

# Plot the associated spatial weighting matrix
s.label(mite.xy, nb = attr(mite.dbmem1, "listw"))

# Plot maps of the first 3 dbMEM eigenfunctions
s.value(mite.xy, mite.dbmem1[,1:3])

# Compute and test associated Moran's I values
# Eigenvalues are proportional to Moran's I

test <- moran.randtest(mite.dbmem1, nrepet = 99)
plot(test$obs, attr(mite.dbmem1, "values"), xlab = "Moran's I", ylab = "Eigenvalues")

# Decreasing values of Moran's I for the successive MEM.
# The red line is the expected value of Moran's I under H0.

plot(test$obs, xlab="MEM rank", ylab="Moran's I")
abline(h=-1/(nrow(mite.xy) - 1), col="red")

# Example 2: Compute only the MEMs with positive eigenvalues (and positive Moran's I)
mite.dbmem2 <- dbmem(mite.xy, thresh=1.012)
# or:  mite.dbmem2 <- dbmem(dist(mite.xy), thresh=1.012, silent=FALSE)
mite.dbmem2

# Examine the eigenvalues
attributes(mite.dbmem2)$values
# or:  attr(mite.dbmem2, "values")

# Examine (any portion of) the dbmem spatial eigenvectors
tmp <- as.matrix(mite.dbmem2)
tmp[1:10,1:6]
}

---

Directional indices of community change

Description

Compute directional indices of community change along coenoclines or time.

Usage

directional.response(mat, method = "overlap", relativize = NULL)

Arguments

mat	A community composition data matrix with sites in rows and species in columns. The direction of the physical process is indicated by the order of the sampling units in mat. The class of mat can be either data.frame or matrix.
method	One of the 11 calculation methods available in the function: "overlap", "gain", "loss", "gaining.turnover", "neutral.turnover", "losing.turnover", "gaining.nestedness", "neutral.nestedness", "losing.nestedness", "gaining.strict.nestedness", "losing.strict.nestedness". The default value is method="overlap".
relativize	Compute relativized indices: relativize="J" for the Jaccard denominator (a+b+c) or relativize="S" for the Sorensen denominator (2*a+b+c). If relativize=NULL, the index is not divided by a denominator.

Details

The output matrix is non-symmetric, meaning that its upper triangle is not the mirror image of the lower triangle. Given the direction of the physical process through space or time indicated by the order of the sampling units, the output matrix mat.out reflects that direction in its non-symmetric presentation, From (rows of the matrix) To (columns of the matrix). Users of the function can extract one or the other of these triangular portions and analyse them separately. See Appendix xx for examples.

Value

A list containing the following results:

• mat.out: A square matrix with the chosen index computed among the sites. Depending on the method chosen, this matrix may be symmetric or non-symmetric. total.t: methods #4 to 6, a matrix with total turnover (b+c); else NA. total.n: For methods #7 to 9, a matrix with total nestedness a+abs(b-c) if a>0; else NA. total.strict.n: For methods #10 and 11, a matrix with total strict nestedness a+abs(b-c) if a>0 and b!=c; else NA. den: For calculation results with Jaccard or Sorensen denominators: a square matrix of denominators. If relativize=NULL, den=NA.

Author(s)

Dénes Schmera [email protected] and Pierre Legendre [email protected]

References

Schmera, D., P. Legendre, T. Eros, M. Toth, E. K. Magyari, B. Baur and J. Podani. 2022. New measures for quantifying directional changes in presence-absence community data. Ecological Indicators 136: 108618. https://doi.org/10.1016/j.ecolind.2022.108618

Verneaux, J. (1973) Cours d'eau de Franche-Comté (Massif du Jura). Recherches écologiques sur le réseau hydrographique du Doubs. Essai de biotypologie. Thèse d'État, Besançon. 1–257.

Examples

# Artificial Example
art <- c(1,1,1,0,0,0,
         0,0,0,1,1,0,
         0,0,0,0,0,1)
art.data <- matrix(art, nrow=3, ncol=6, byrow=TRUE)

art.out <- directional.response(art.data, method="overlap",relativize=NULL)

# Real data example: the Doubs River fish data (Verneaux 1973), available in ade4.
# 30 sites, 27 species. No fish had been caught at site 8; remove that site
if(require("ade4", quietly = TRUE)) {

data(doubs)
dim(doubs$fish)   
fish <- doubs$fish[-8,] 
dim(fish)
doubs.out <- directional.response(fish, method="gain", relativize="S")
}

---

Dissimilarity matrices for community composition data

Description

Compute dissimilarity indices for ecological data matrices. The dissimilarity indices computed by this function are those described in Legendre and De Cáceres (2013). In the name of the function, 'ldc' stands for the author's names. Twelve of these 21 indices are not readily available in other R package functions; four of them can, however, be computed in two computation steps in vegan.

Usage

dist.ldc(Y, method = "hellinger", binary = FALSE, samp = TRUE, silent = FALSE)

Arguments

Y	Community composition data. The object class can be either data.frame or matrix.
method	One of the 21 dissimilarity coefficients available in the function: "hellinger", "chord", "log.chord", "chisquare", "profiles", "percentdiff", "ruzicka", "divergence", "canberra", "whittaker", "wishart", "kulczynski", "jaccard", "sorensen", "ochiai", "ab.jaccard", "ab.sorensen", "ab.ochiai", "ab.simpson", "euclidean", "manhattan", "modmeanchardiff". See Details. Names can be abbreviated to a non-ambiguous set of first letters. Default: method="hellinger".
binary	If binary=TRUE, the data are transformed to presence-absence form before computation of the dissimilarities. Default value: binary=FALSE, except for the Jaccard, Sørensen and Ochiai indices where binary=TRUE.
samp	If samp=TRUE, the abundance-based distances (ab.jaccard, ab.sorensen, ab.ochiai, ab.simpson) are computed for sample data. If samp=FALSE, binary indices are computed for true population data.
silent	If silent=FALSE, informative messages sent to users will be printed to the R console. Use silent=TRUE is called on a numerical simulation loop, for example.

Details

The dissimilarities computed by this function are the following. Indices i and k designate two rows (sites) of matrix Y, j designates a column (species). D[ik] is the dissimilarity between rows i and k. p is the number of columns (species) in Y; pp is the number of species present in one or the other site, or in both. y[i+] is the sum of values in row i; same for y[k+]. y[+j] is the sum of values in column j. y[++] is the total sum of values in Y. The indices are computed by functions written in C for greater computation speed with large data matrices.

• Group 1 - D computed by transformation of Y followed by Euclidean distance

  • Hellinger D, D[ik] = sqrt(sum((sqrt(y[ij]/y[i+])-sqrt(y[kj]/y[k+]))^2))

  • chord D, D[ik] = sqrt(sum((y[ij]/sqrt(sum(y[ij]^2))-y[kj]/sqrt(sum(y[kj]^2)))^2))

  • log-chord D, D[ik] = chord D[ik] computed on log(y[ij]+1)-transformed data (Legendre and Borcard 2018)

  • chi-square D, D[ik] = sqrt(y[++] sum((1/j[+j])(y[ij]/y[i+]-y[kj]/y[k+])^2))

  • species profiles D, D[ik] = sqrt(sum((y[ij]/y[i+]-y[kj]/y[k+])^2))

• Group 2 - Other D functions appropriate for beta diversity studies where A = sum(min(y[ij],y[kj])), B = y[i+]-A, C = y[k+]-A

  • percentage difference D (aka Bray-Curtis), D[ik] = (sum(abs(y[ij]-y[k,j])))/(y[i+]+y[k+]) or else, D[ik] = (B+C)/(2A+B+C)

  • Ružička D, D[ik] = 1-(sum(min(y[ij],y[kj])/sum(max(y[ij],y[kj])) or else, D[ik] = (B+C)/(A+B+C)

  • coeff. of divergence D, D[ik] = sqrt((1/pp)sum(((y[ij]-y(kj])/(y[ij]+y(kj]))^2))

  • Canberra metric D, D[ik] = (1/pp)sum(abs(y[ij]-y(kj])/(y[ij]+y(kj]))

  • Whittaker D, D[ik] = 0.5*sum(abs(y[ij]/y[i+]-y(kj]/y[k+]))

  • Wishart D, D[ik] = 1-sum(y[ij]y[kj])/(sum(y[ij]^2)+sum(y[kj]^2)-sum(y[ij]y[kj]))

  • Kulczynski D, D[ik] = 1-0.5((sum(min(y[ij],y[kj])/y[i+]+sum(min(y[ij],y[kj])/y[k+]))

• Group 3 - Classical indices for binary data; they are appropriate for beta diversity studies. Value a is the number of species found in both i and k, b is the number of species in site i not found in k, and c is the number of species found in site k but not in i. The D matrices are square-root transformed, as in dist.binary of ade4; the user-oriented reason for this transformation is explained below.

  • Jaccard D, D[ik] = sqrt((b+c)/(a+b+c))

  • Sørensen D, D[ik] = sqrt((b+c)/(2a+b+c))

  • Ochiai D, D[ik] = sqrt(1 - a/sqrt((a+b)(a+c)))

• Group 4 - Abundance-based indices of Chao et al. (2006) for quantitative abundance data. These functions correct the index for species that have not been observed due to sampling errors. For the meaning of the U and V notations, see Chao et al. (2006, section 3). When samp=TRUE, the abundance-based distances (ab.jaccard, ab.sorensen, ab.ochiai, ab.simpson) are computed for sample data. If samp=FALSE, indices are computed for true population data. - Do not use indices of group 4 with samp=TRUE on presence-absence data; the indices are not meant to accommodate this type of data. If samp=FALSE is used with presence-absence data, the indices are the regular Jaccard, Sørensen, Ochiai and Simpson indices. On output, however, the D matrices are not square-rooted, contrary to the Jaccard, Sørensen and Ochiai indices in section 3 which are square-rooted.

  • abundance-based Jaccard D, D[ik] = 1-(UV/(U+V-UV))

  • abundance-based Sørensen D, D[ik] = 1-(2UV/(U+V))

  • abundance-based Ochiai D, D[ik] = 1-sqrt(UV)

  • abundance-based Simpson D, D[ik] = 1-(UV/(UV+min((U-UV),(V-UV))))

• Group 5 - General-purpose dissimilarities that do not have an upper bound (maximum D value). They are inappropriate for beta diversity studies.

  • Euclidean D, D[ik] = sqrt(sum(y[ij]-y[kj])^2)

  • Manhattan D, D[ik] = sum(abs(y[ij] - y[ik]))

  • modified mean character difference, D[ik] = (1/pp) sum(abs(y[ij] - y[ik]))

The properties of all dissimilarities available in this function (except Ružička D) were described and compared in Legendre & De Cáceres (2013), who showed that most of these dissimilarities are appropriate for beta diversity studies. Inappropriate are the Euclidean, Manhattan, modified mean character difference, species profile and chi-square distances. Most of these dissimilarities have a maximum value of either 1 or sqrt(2). Three dissimilarities (Euclidean, Manhattan, Modified mean character difference) do not have an upper bound and are thus inappropriate for beta diversity studies. The chi-square distance has an upper bound of sqrt(2*(sum(Y))).

The Euclidean, Hellinger, chord, chi-square and species profiles dissimilarities have the property of being Euclidean, meaning that they never produce negative eigenvalues in principal coordinate analysis. The Canberra, Whittaker, percentage difference, Wishart and Manhattan coefficients are Euclidean when they are square-root transformed (Legendre & De Cáceres 2013, Table 2). The distance forms (1-S) of the Jaccard, Sørensen and Ochiai similarity (S) coefficients are Euclidean after taking the square root of (1-S) (Legendre & Legendre 2012, Table 7.2). The D matrices resulting from these three coefficients are outputted in the form sqrt(1-S), as in function dist.binary of ade4, because that form is Euclidean and will thus produce no negative eigenvalues in principal coordinate analysis.

The Hellinger, chord, chi-square and species profile dissimilarities are computed using the two-step procedure developed by Legendre & Gallagher (2001). The data are first transformed using either the row marginals, or the row and column marginals in the case of the chi-square distance. The dissimilarities are then computed from the transformed data using the Euclidean distance formula. As a consequence, these four dissimilarities are necessarily Euclidean. D matrices for other binary coefficients can be computed in two ways: either by using function dist.binary of ade4, or by choosing option binary=TRUE, which transforms the abundance data to binary form, and using one of the quantitative indices of the present function. Table 1 of Legendre & De Cáceres (2013) shows the incidence-based (presence-absence-based) indices computed by the various indices using binary data.

The Euclidean distance computed on untransformed presence-absence or abundance data produces non-informative and incorrect ordinations, as shown in Legendre & Legendre (2012, p. 300) and in Legendre & De Cáceres (2013). However, the Euclidean distance computed on log-transformed abundance data produces meaningful ordinations in principal coordinate analysis (PCoA). Nonetheless, it is easier to compute a PCA of log-transformed abundance data instead of a PCoA; the resulting ordination with scaling 1 will be meaningful. Messages are printed to the R console indicating the Euclidean status of the computed dissimilarity matrices. Note that for the chi-square distance, the columns that sum to zero are eliminated before calculation of the distances, thus preventing divisions by zero in the calculation of the chi-square transformation.

Value

A dissimilarity matrix, with class dist.

Author(s)

Pierre Legendre [email protected] and Naima Madi

References

Chao, A., R. L. Chazdon, R. K. Colwell and T. J. Shen. 2006. Abundance-based similarity indices and their estimation when there are unseen species in samples. Biometrics 62: 361-371.

Legendre, P. and D. Borcard. 2018. Box-Cox-chord transformations for community composition data prior to beta diversity analysis. Ecography 41: 1820-1824.

Legendre, P. and M. De Cáceres. 2013. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters 16: 951-963.

Legendre, P. and E. D. Gallagher, E.D. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129: 271-280.

Legendre, P. and Legendre, L. 2012. Numerical Ecology. 3rd English edition. Elsevier Science BV, Amsterdam.

Examples

if(require("vegan", quietly = TRUE)) {
data(mite)
mat1  = as.matrix(mite[1:10, 1:15])   # No column has a sum of 0
mat2 = as.matrix(mite[61:70, 1:15])   # 7 of the 15 columns have a sum of 0

#Example 1: compute Hellinger distance for mat1
D.out = dist.ldc(mat1,"hellinger")

#Example 2: compute chi-square distance for mat2
D.out = dist.ldc(mat2,"chisquare")

#Example 3: compute percentage difference dissimilarity for mat2
D.out = dist.ldc(mat2,"percentdiff")

}

---

Perform a test of the shared space-environment fraction of a variation partitioning using torus-translation (TT) or Moran Spectral Randomisation (MSR)

Description

The function uses two different spatially-constrained null models to test the shared space-environment fraction (SSEF, or fraction [b]) of a variation partitioning of two explanatory components.

Usage

envspace.test(
  spe,
  env,
  coord,
  MEM.spe,
  listw.env,
  MEM.autocor = c("positive", "negative", "all"),
  regular = FALSE,
  nperm = 999,
  MSR.method = "singleton",
  alpha = 0.05
)

Arguments

spe	Vector, matrix, or dataframe of response variable(s) (e.g. species abundances)
env	Vector, matrix, or dataframe of environmental variables (rows = sites, columns = variables)
coord	Matrix or dataframe of spatial coordinates of the sampled sites
MEM.spe	Matrix or dataframe of spatial predictors (MEM variables) selected for spe
listw.env	An object of class listw (spatial weights) created by the functions of the spdep package or returned by listw.candidates
MEM.autocor	A string indicating the type of spatial structure of interest for env ("positive", "negative", or "all", for positive, negative, or both types of spatial autocorrelations, respectively); Default is "positive"
regular	Logical argument indicating whether a torus-translation test will be performed, in addition to the MSR. Set to TRUE only if the sampling design is regular (same number of sites on each line, same number of sites on each column). Set to FALSE otherwise; Default is FALSE
nperm	Number of permutations performed; Default is 999
MSR.method	Algorithm of msr to be used to perform the MSR. The three available procedures are "singleton" (default), "pair", and "triplet" (see msr for details)
alpha	Threshold value of null hypothesis rejection for the test of a spatial structure in the environment, and for the shared environment-space fraction of the variation partitioning; Default is 0.05

Details

The function tests the SSEF (also known as fraction [b]) of a variation partitioning of a response variable or matrix (y) between an environmental and a spatial component (env, and MEM.spe, respectively). The SSEF is the explained variation of y shared by env and MEM.spe. The adjusted R-squared (Peres-Neto et al. 2006; R2adj) of the SSEF is not an actual R2, as it is computed by subtracting the adjusted R2adj of other fractions and therefore has zero degree of freedom (Legendre and Legendre 2012). The SSEF can therefore not be computed in the classical way (residuals permutation; Anderson and Legendre 1999, Legendre and Legendre 2012).

The function envspace.test provides two ways of testing this fraction, that is, spatially-constrained null models based either on a torus-translation test (TT) (for regular sampling designs only), or on Moran spectral randomizations (MSR) (for any type of sampling design). The test of the SSEF should only be performed if both the global models of y against all the environmental variables and against all spatial variables are significant (see Bauman et al. 2018c). The function first checks whether the environment displays significant spatial structures, and then proceeds to the test of the SSEF if this condition is fulfilled (details in Bauman et al. 2018c).

spe can be a vector or a multicolumn matrix or dataframe (multivariate response data). If multivariate, it is greatly advised to transform spe prior to performing the variation partitioning and testing the SSEF (e.g., Hellinger transformation; see Legendre and Gallagher 2001).

MEM.spe is a set of spatial predictors (MEM variables). It is recommended to be a well-defined subset of MEM variables selected among the complete set generated from the spatial weighting matrix (SWM) (see review about spatial eigenvector selection in Bauman et al. 2018a). Optimising the selection of a subset of forward-selected MEM variables among a set of candidate SWMs has been shown to increase statistical power as well as R2-estimation accuracy (Bauman et al. 2018b). To do so, MEM.spe can be generated using listw.candidates followed by listw.select. If a SWM has already been selected in another way, then mem.select can be used to generate the MEM variables and to select an optimal subset among them, which can then be used as MEM.spe in envspace.test (see Details of function mem.select). listw.env corresponds to the SWM that will be used to test for a spatial structure in env, and to build the MEM variables for the MSR test. The choice of the SWM for env can also be optimised with listw.select. The SWMs selected for spe and env should be optimised separately to best model the spatial structure of both spe and env (see example).

To verify that env displays a significant spatial pattern, prior to performing the test of the SSEF, a residuals permutation test is performed on the global set of MEM variables (generated internally from listw.env) associated to the type of spatial structure of interest (see argument MEM.autocor). This test is performed with mem.select. The choice of MEM.autocor should be made according to the MEM.autocor argument used to build MEM.spe.

env is a dataset of environmental variables chosen by the user. We recommend dealing with collinearity issues prior to performing the variation partitioning and the test of the SSEF (see Dormann et al. 2013 for a review of methods to cope with collinearity).

regular is a logical argument indicating whether a TT test should be performed instead of the MSR to test the SSEF. Since the TT can only be performed on regular sampling designs, regular should only be set to TRUE if the sampling design is either a transect, or a grid displaying the same number of sites for all lines and columns (although the number of sites per column can differ from the number of sites per line).

listw.env is the SWM used by the MSR to generate spatially-constrained null environmental variables. It should ideally be a SWM optimised on the basis of env using the function listw.select, with the argument method = "global" (see Details of function mem.select for an explanation). This will allow detecting the spatial structures of env as accurately as possible, hence allowing MSR to generate null environmental variables as spatially faithful to the original ones. It is also on the basis of listw.env that MEM variables will be generated to test whether env is spatially structured (i.e. global test) prior to perform the test of the SSEF.

It is worth mentioning that, although a significant SSEF may provide evidence of an induced spatial dependence (Bauman et al. 2018c), a non-significant SSEF only indicates that no induced spatial dependence could be detected in relation with the chosen environmental variables. This does not exclude that this effect may exist with respect to some unmeasured variables.

Value

If the condition of env being spatially structured is fulfilled, the test is performed and the function returns an object of class randtest containing the results of the test.

Author(s)

David Bauman and Jason Vleminckx, [email protected], [email protected]

References

Anderson M. and Legendre P. (1999) An empirical comparison of permutation methods for tests of partial regression coefficients in a linear model. Journal of Statistical Computation and Simulation, 62(3), 271–303

Bauman D., Drouet T., Dray S. and Vleminckx J. (2018a) Disentangling good from bad practices in the selection of spatial or phylogenetic eigenvectors. Ecography, 41, 1–12

Bauman D., Fortin M-J, Drouet T. and Dray S. (2018b) Optimizing the choice of a spatial weighting matrix in eigenvector-based methods. Ecology

Bauman D., Vleminckx J., Hardy O., Drouet T. (2018c) Testing and interpreting the shared space-environment fraction in variation partitioning analyses of ecological data. Oikos

Blanchet G., Legendre P. and Borcard D. (2008) Forward selection of explanatory variables. Ecology, 89(9), 2623–2632

Legendre P., Gallagher E.D. (2001) Ecologically meaningful transformations for ordination of species data. Oecologia, 129(2), 271–280

Legendre P. and Legendre L. (2012) Numerical Ecology, Elsevier, Amsterdam

Peres-Neto P., Legendre P., Dray S., Borcard D. (2006) Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology, 87(10), 2614–2625

Peres-Neto P. and Legendre P. (2010) Estimating and controlling for spatial structure in the study of ecological communities. Global Ecology and Biogeography, 19, 174–184

See Also

varpart, listw.select, listw.candidates, mem.select

Examples

if(require(vegan)) { 
# Illustration of the test of the SSEF on the oribatid mite data
# (Borcard et al. 1992, 1994 for details on the dataset):
# Community data (response matrix):
data(mite)
# Hellinger-transformation of the community data (Legendre and Gallagher 2001):
Y <- decostand(mite, method = "hellinger")
# Environmental explanatory dataset:
data(mite.env)
# We only use two numerical explanatory variables:
env <- mite.env[, 1:2]
dim(Y)
dim(env)
# Coordinates of the 70 sites:
data(mite.xy)
coord <- mite.xy

### Building a list of candidate spatial weighting matrices (SWMs) for the 
### optimisation of the SWM selection, separately for 'Y' and 'env':
# We create five candidate SWMs: a connectivity matrix based on a Gabriel graphs, on
# a minimum spanning tree (i.e., two contrasted graph-based SWMs), either
# not weighted, or weighted by a linear function decreasing with the distance),
# and a distance-based SWM corresponding to the connectivity and weighting
# criteria of the original PCNM method:
candidates <- listw.candidates(coord, nb = c("gab", "mst", "pcnm"), weights = c("binary",
                                                                                "flin"))
### Optimisation of the selection of a SWM:
# SWM for 'Y' (based on the best forward-selected subset of MEM variables):
modsel.Y <- listw.select(Y, candidates, method = "FWD", MEM.autocor = "positive",
                         p.adjust = TRUE)
                         
names(candidates)[modsel.Y$best.id]                 # Best SWM selected
modsel.Y$candidates$Pvalue[modsel.Y$best.id]        # Adjusted p-value of the global model
modsel.Y$candidates$N.var[modsel.Y$best.id]         # Nb of forward-selected MEM variables
modsel.Y$candidates$R2Adj.select[modsel.Y$best.id]  # Adjusted R2 of the selected MEM var.

# SWM for 'env' (method = "global" for the optimisation, as all MEM variables are required
# to use MSR):
modsel.env <- listw.select(env, candidates, method = "global", MEM.autocor = "positive",
                           p.adjust = TRUE)

names(candidates)[modsel.env$best.id]                  # Best SWM selected
modsel.env$candidates$Pvalue[modsel.env$best.id]       # Adjusted p-value of the global model
modsel.env$candidates$N.var[modsel.env$best.id]        # Nb of forward-selected MEM variables
modsel.env$candidates$R2Adj.select[modsel.env$best.id] # Adjusted R2 of the selected MEM var.

### We perform the variation partitioning:
# Subset of selected MEM variables within the best SWM:
MEM.spe <- modsel.Y$best$MEM.select

VP <- varpart(Y, env, MEM.spe)
plot(VP)

# Test of the shared space-environment fraction (fraction [b]):
SSEF.test <- envspace.test(Y, env, coord, MEM.spe, 
                           listw.env = candidates[[modsel.env$best.id]], 
                           regular = FALSE, nperm = 999)
SSEF.test

# The SSEF is highly significant, indicating a potential induced spatial dependence.
}

---

Forward selection with multivariate Y using permutation under reducel model

Description

Performs a forward selection by permutation of residuals under reduced model. Y can be multivariate.

Usage

forward.sel(
  Y,
  X,
  K = nrow(X) - 1,
  R2thresh = 0.99,
  adjR2thresh = 0.99,
  nperm = 999,
  R2more = 0.001,
  alpha = 0.05,
  Xscale = TRUE,
  Ycenter = TRUE,
  Yscale = FALSE,
  verbose = TRUE
)

Arguments

Y	Response data matrix with n rows and m columns containing quantitative variables
X	Explanatory data matrix with n rows and p columns containing quantitative variables
K	Maximum number of variables to be selected. The default is one minus the number of rows
R2thresh	Stop the forward selection procedure if the R-square of the model exceeds the stated value. This parameter can vary from 0.001 to 1
adjR2thresh	Stop the forward selection procedure if the adjusted R-square of the model exceeds the stated value. This parameter can take any value (positive or negative) smaller than 1
nperm	The number of permutation to be used.The default setting is 999 permutation.
R2more	Stop the forward selection procedure if the difference in model R-square with the previous step is lower than R2more. The default setting is 0.001
alpha	Significance level. Stop the forward selection procedure if the p-value of a variable is higher than alpha. The default is 0.05 is TRUE
Xscale	Standardize the variables in table X to variance 1. The default setting is TRUE
Ycenter	Center the variables in table Y. The default setting is TRUE
Yscale	Standardize the variables in table Y to variance 1. The default setting is FALSE.
verbose	If 'TRUE' more diagnostics are printed. The default setting is TRUE

Details

The forward selection will stop when either K, R2tresh, adjR2tresh, alpha and R2more has its parameter reached.

Value

A dataframe with:

variables	The names of the variables
order	The order of the selection of the variables
R2	The R2 of the variable selected
R2Cum	The cumulative R2 of the variables selected
AdjR2Cum	The cumulative adjusted R2 of the variables selected
F	The F statistic
pval	The P-value statistic

Note

Not yet implemented for CCA (weighted regression) and with covariables.

Author(s)

Stephane Dray [email protected]

References

Canoco manual p.49

Examples

x <- matrix(rnorm(30),10,3)
y <- matrix(rnorm(50),10,5)
    
forward.sel(y,x,nperm=99, alpha = 0.5)

---

Parametric forward selection of explanatory variables in regression and RDA

Description

If Y is univariate, this function implements FS in regression. If Y is multivariate, this function implements FS using the F-test described by Miller and Farr (1971). This test requires that (i) the Y variables be standardized, and (ii) the error in the response variables be normally distributed (to be verified by the user).

Usage

forward.sel.par(
  Y,
  X,
  alpha = 0.05,
  K = nrow(X) - 1,
  R2thresh = 0.99,
  R2more = 0.001,
  adjR2thresh = 0.99,
  Yscale = FALSE,
  verbose = TRUE
)

Arguments

Y	Response data matrix with n rows and m columns containing quantitative variables
X	Explanatory data matrix with n rows and p columns containing quantitative variables
alpha	Significance level. Stop the forward selection procedure if the p-value of a variable is higher than alpha. The default is 0.05
K	Maximum number of variables to be selected. The default is one minus the number of rows
R2thresh	Stop the forward selection procedure if the R-square of the model exceeds the stated value. This parameter can vary from 0.001 to 1
R2more	Stop the forward selection procedure if the difference in model R-square with the previous step is lower than R2more. The default setting is 0.001
adjR2thresh	Stop the forward selection procedure if the adjusted R-square of the model exceeds the stated value. This parameter can take any value (positive or negative) smaller than 1
Yscale	Standardize the variables in table Y to variance 1. The default setting is FALSE. The setting is automatically changed to TRUE if Y contains more than one variable. This is a validity condition for the parametric test of significance (Miller and Farr 1971)
verbose	If 'TRUE' more diagnostics are printed. The default setting is TRUE

Details

The forward selection will stop when either K, R2tresh, adjR2tresh, alpha and R2more has its parameter reached.

Value

A dataframe with:

variables	The names of the variables
order	The order of the selection of the variables
R2	The R2 of the variable selected
R2Cum	The cumulative R2 of the variables selected
AdjR2Cum	The cumulative adjusted R2 of the variables selected
F	The F statistic
pval	The P-value statistic

Author(s)

Pierre Legendre [email protected] and Guillaume Blanchet

References

Miller, J. K. & S. D. Farr. 1971. Bimultivariate redundancy: a comprehensive measure of interbattery relationship. Multivariate Behavioral Research, 6, 313–324.

Examples

x <- matrix(rnorm(30),10,3)
y <- matrix(rnorm(50),10,5)
    
forward.sel.par(y,x, alpha = 0.5)

---

Compute the maximum distance of the minimum spanning tree based on a distance matrix

Description

It is used to select a truncation value for the dbMEM approach. It returns the minimum value that keep all samples connected.

Usage

give.thresh(matdist)

Arguments

matdist

A distance matrix (class dist or matrix)

Value

The maximum distance in the minimum spanning tree.

Author(s)

Stéphane Dray [email protected]

Examples

xy <- matrix(rnorm(60),30,2)
dxy <- dist(xy)
th <- give.thresh(dxy)

---

Global and local tests

Description

These two Monte Carlo tests are used to assess the existence of 'global' and 'local' spatial structures, corresponding respectively to positive and negative Moran's I .

Usage

global.rtest(X, listw, k = 1, nperm = 499)

Arguments

X	a data matrix, with variables in columns
listw	a list of weights of class listw. Can be obtained easily using the function chooseCN.
k	integer: the number of highest $R^2$ summed to form the test statistics
nperm	integer: the number of randomisations to be performed.

Details

They rely on the decomposition of a data matrix X into global and local components using multiple regression on Moran's Eigenvector Maps (MEMs). They require a data matrix (X) and a list of weights derived from a connection network. X is regressed onto global MEMs (U+) in the global test and on local ones (U-) in the local test. One mean $R^2$ is obtained for each MEM, the k highest being summed to form the test statistic.

The reference distribution of these statistics are obtained by randomly permuting the rows of X.

These tests were originally part of the adegenet package for R.

Value

An object of class randtest.

Author(s)

Thibaut Jombart [email protected]

References

Jombart, T., Devillard, S., Dufour, A.-B. and Pontier, D. 2008. Revealing cryptic spatial patterns in genetic variability by a new multivariate method. Heredity, 101, 92–103. doi: 10.1038/hdy.2008.34.

Examples

# wait for a generic dataset

---

Compute LCBD from any D matrix

Description

Compute LCBD indices (Legendre & De Cáceres 2013) from a symmetric dissimilarity matrix (D) or from a beta component matrix (Repl, RichDiff or AbDiff, or Nes) (Legendre 2014).

Usage

LCBD.comp(D, sqrt.D = TRUE, save.D = FALSE)

Arguments

D	A dissimilarity or beta diversity component matrix, class dist or matrix.
sqrt.D	Take the square root of the dissimilarities in matrix D before computing the LCBD indices.
save.D	If save.D is TRUE, the dissimilarity matrix will appear in the output list.

Details

Use sqrt.D = TRUE when computing LCBD indices for most of the replacement and richness/abundance difference indices computed by function beta.div.comp, as well as for the corresponding D matrices. See Table S1.4 in Appendix S1 of Legendre (2014) to identify the matrices that are Euclidean without taking the square root of the individual values. Only the RichDiffS (for presence-absence data) and AbDiff abundance data) of the Sørensen group in the Podani family have that property. In all other cases, use sqrt.D = TRUE.

When computing LCBD from a D matrix, use sqrt = TRUE if the D matrix is not Euclidean. The Euclidean property can be checked with function is.euclid of ade4.

BDtotal statistics are comparable among data sets having the same or different numbers of sampling units (n), provided that the sampling units are of the same size or represent the same sampling effort and that BDtotal is computed with the same D index.

Function LCBD.comp produces the same (SStotal, BDtotal, LCBD) results as function beta.div. Note, however, that the latter produces other interesting results (p.LCBD, SCBD). Function LCBD.comp should then only be used to compute LCBD indices from dissimilarity matrices that cannot be computed by function beta.div, e.g. genetic D matrices, or from replacement and richness difference matrices produced by function beta.div.comp. Significance of the LCBD indices cannot be tested when their calculation starts from a D matrix because the testing procedure involves permutation of the columns of raw data.

Value

A list containing the following results:

• beta: Total sum of squares and total beta diversity [= Var(Y)] of the data matrix.

• LCBD: Vector of Local contributions to beta diversity (LCBD) for the sites.

• D: The input dissimilarity matrix, class dist; only if save.D=TRUE

.

Author(s)

Pierre Legendre [email protected]

References

Legendre, P. 2014. Interpreting the replacement and richness difference components of beta diversity. Global Ecology and Biogeography 23: 1324-1334.

Legendre, P. & M. De Cáceres. 2013. Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters 16: 951-963.

Examples

### Example 1
### Compute the Hellinger distance, then the LCBD indices.
if(require("vegan", quietly = TRUE)){
data(mite)
mite.hel = decostand(mite, "hellinger")
mite.D = dist(mite.hel)
out.mite.D = LCBD.comp(mite.D, sqrt.D=FALSE)
}

### Example 2
if(require("ade4", quietly = TRUE) & require("adegraphics", quietly = TRUE)){
data(doubs)
fish.sp = doubs$fish[-8,]   # Fish data; site 8 is removed because no fish were caught

out.comp = beta.div.comp(fish.sp, coef="S", quant=TRUE)

out.fish.D = LCBD.comp(out.comp$D, sqrt.D=TRUE)   # out.comp.D is not Euclidean
out.fish.D$beta
out.fish.Repl = LCBD.comp(out.comp$repl, sqrt.D=TRUE)   # out.comp$repl is not Euclidean
out.fish.Repl$beta
out.fish.AbDiff = LCBD.comp(out.comp$rich, sqrt.D=FALSE)   # out.comp$rich is Euclidean
out.fish.AbDiff$beta

### Plot maps of the LCBD indices
fish.xy = doubs$xy[-8,]   # Geographic coordinates; site 8 removed because no fish were caught

# Map of LCBD indices for %difference dissimilarity
s.value(fish.xy, out.fish.D$LCBD, method="size", symbol = "circle",
col = c("white", "brown"), main = "Doubs fish LCBD, %difference D")

# Map of LCBD indices for replacement component of %difference dissimilarity
s.value(fish.xy, out.fish.Repl$LCBD, method="size", symbol = "circle",
col = c("white", "brown"), main = "Doubs fish replacement LCBD")

# Map of LCBD indices for abundance difference component of %difference dissimilarity
s.value(fish.xy, out.fish.AbDiff$LCBD, method="size", symbol = "circle", 
col = c("white", "brown"), main = "Doubs fish abundance diff. LCBD")
}


if(require("ade4", quietly = TRUE) & require("betapart", quietly = TRUE)){
### Example 3
### This example requires packages \code{"betapart"} and \code{"ade4"} for data. 
### For the Baselga-family indices, the same partitioning results are obtained using
### (1) beta.div.comp or (2) beta.pair.abund() of \code{"betapart"} and LCBD.comp()

data(doubs)   # Data available in \code{"ade4"}
fish.sp = doubs$fish[-8,]   
# Fish data; site 8 is removed because no fish were caught
# We use abundance data in this example, not presence-absence data

# Partition into Baselga-family replacement and nestedness components 
# using \code{"beta.div.comp"} with the percentage difference index (aka Bray-Curtis)
out.comp = beta.div.comp(fish.sp, coef="BS", quant=TRUE)
out.comp$part

# Compute the D and component matrices using \code{"beta.pair.abund"}
out3 = beta.pair.abund(fish.sp, index.family = "bray")
summary(out3)

is.euclid(out3$beta.bray)    # D matrix out3$beta.bray is not Euclidean
out3.D = LCBD.comp(out3$beta.bray, sqrt.D=TRUE)
out3.D$beta
# Compare BDtotal here to BDtotal in out.comp$part (above)

out3.Repl = LCBD.comp(out3$beta.bray.bal, sqrt.D=TRUE)
out3.Repl$beta
# Compare BDtotal here to RichDiff in out.comp$part (above)

out3.AbDiff = LCBD.comp(out3$beta.bray.gra, sqrt.D=TRUE)
out3.AbDiff$beta
# Compare BDtotal here to RichDiff/Nes in out.comp$part (above)
}

---

Function to create a list of spatial weighting matrices

Description

This function is a user-friendly way to create a list of one or several spatial weighting matrices (SWM) by selecting a set of predefined connectivity and weighting matrices (B and A matrices, respectively).

Usage

listw.candidates(
  coord,
  style = "B",
  nb = c("del", "gab", "rel", "mst", "pcnm", "dnear"),
  d1 = 0,
  d2,
  weights = c("binary", "flin", "fup", "fdown"),
  y_fdown = 5,
  y_fup = 0.5
)

Arguments

coord	Vector, matrix, or dataframe of point coordinates
style	Coding scheme style (see nb2listw of the spdep package). Can take values 'W', 'B', 'C', 'U', 'minmax', and 'S'; default is 'B'
nb	Defines how the B matrix (connectivity) is build: del Delaunay triangulation gab Gabriel's graph rel Relative neighbourhood graph mst Minimum spanning tree pcnm Distance-based SWM based on the principal coordinates of neighbour matrices (PCNM) criteria (see 'Details') dnear Distance-based
d1	Only considered if nb = "dnear". A single value defining the distance beyond which two sites are connected (i.e., minimum distance between two neighbor sites). The default value is 0 (no constraint on the min distance). d1 must be smaller than d2
d2	Only considered if nb = "dnear". It defines the connectivity distance threshold below which two sites are connected (i.e., maximum distance between two neighbors. It can either be a single value or a vector of values, in which case a different SWM will be generated for each threshold value. The default value is the minimum distance keeping all points connected (i.e., the largest edge of the minimum spanning tree)
weights	Defines how the A matrix (weighths) is build: binary without weights flin Linear weighting function fdown Concave-down weighting function(see Details below) fup Concave-up weighting function (see Details below)
y_fdown	Single value or vector of values of the y parameter in the concave-down weighting function; default is 5
y_fup	Single value or vector of values of the y parameter in the concave-up weighting function; default is 0.5

Details

The function allows constructing SWMs based on any combination of B and A matrices. The B matrices are either graph-based or distance-based. The function proposes the Delaunay triangulation, Gabriel graph, relative neighbourhood graph, and the minimum spanning tree criteria to build a graph-based B matrix. Distance-based SWMs can be built with the principal coordinates of neighbour matrices (PCNM; Borcard and Legendre 2002) criteria (see details below), or using another threshold distance to define the connected site pairs. The A matrix can be based on a binary, linear, concave-down, or concave-up function. The linear, concave-down, and concave-up weighting functions are defined by $1 - (D/dmax)$, $1 - (D/dmax)^y$, and $1 / D^y$, respectively, where D is the euclidean distance between the two sites considered, dmax is the maximum euclidean distance between two sites, and y is a user-defined parametre that can either be a single value or a vector of values. The choice nb = "pcnm" consists in constructing a distance-based SWM based on the largest edge of the minimum spanning tree as a connectivity distance threshold, and then by weighting the links by the function $1-(D/(4*t))^2$, where D is the euclidean distance between the sites, and t is the distance threshold below which two sites are considered connected (Dray et al. 2006). As optimizing the choice of a SWM has to be done with a p-value correction depending on the number of candidate SWMs tested (see function listw.select), Bauman et al. (2018) strongly encouraged plotting the concave-down and concave-up weighting functions with several parametre values in order to only choose the realistic ones to build the candidate W matrices (e.g., ranging between 0.1 and 1 for the concave-up function, as values over 1 would make no ecological sense). First visualizing the connectivity schemes with the listw.explore function may also help choosing the B matrices to select for the listw.candidates function.

Spatial eigenvectors can be generated from any candidate SWM obtained by listw.candidates using scores.listw, or can be generated and tested (recommended option for real data analysis) using mem.select. If several SWMs were created, the selection of an optimized SWM can be made using listw.select.

Value

A list of SWMs. Each element of the list was built by nb2listw (package spdep) and therefore is of class listw and nb. The name of each element of the list (SWM) is composed of the corresponding B and A matrices, followed (if any) by the y parameter value of the weighting function.

Author(s)

David Bauman ([email protected] or [email protected]) and Stéphane Dray

References

Bauman D., Fortin M-J., Drouet T. and Dray S. (2018) Optimizing the choice of a spatial weighting matrix in eigenvector-based methods. Ecology

Borcard D. and Legendre P. (2002) All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecological Modelling, 153, 51–68

Dray S., Legendre P. and Peres-Neto P. R. (2006) Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM). Ecological Modelling, 196, 483–493

See Also

listw.explore, scores.listw, mem.select, listw.select

Examples

### Create 100 random sampling locations in a squared grid of 120 x 120:
xy <- matrix(nrow = 100, ncol = 2)
xy[, 1] <- sample(c(1:120), 100, replace = FALSE)
xy[, 2] <- sample(c(1:120), 100, replace = FALSE)
### The function listw.candidates is used to build the spatial weighting matrices that
### we want to test and compare (with the listw.select function). We test a Gabriel's graph, 
### a minimum spanning tree, and a distance-based connectivity defined by a threshold
### distance corresponding to the smallest distance keeping all sites connected (i.e., 
### the defaut value of d2). These connectivity matrices are then either not weighted 
### (binary weighting), or weighted by the linearly decreasing function:
candidates <- listw.candidates(coord = xy, nb = c("gab", "mst", "dnear"), 
                               weights = c("binary", "flin"))
names(candidates)                              
plot(candidates[[1]], xy)
plot(candidates[[3]], xy)
### Construction of a different list of spatial weighting matrices. This time, the
### connexions are defined by a distance-based criterion based on the same threshold
### value, but the connections are weighted by the concave-down function with a y parameter
### varying between 2 and 5, and a concave-up function with a y parametre of 0.2.
candidates2 <- listw.candidates(coord = xy, nb = "dnear", weights = c("fdown", "fup"),
                                y_fdown = 1:5, y_fup = 0.2)
### Number of spatial weighting matrices generated:
length(candidates2) 
### A single SWM can also easily be generated with listw.candidates:
lw <- listw.candidates(xy, nb = "gab", weights = "bin")
plot(lw[[1]], xy)

### Generating MEM variables from an object of listw.candidates with scores.listw:
MEM <- scores.listw(lw[[1]])
### See functions mem.select and listw.select for examples of how to use an object
### created by listw.candidates with these functions.

---

Interactive tool to generate R code that creates a spatial weighting matrix

Description

Interactive tool to generate R code that creates a spatial weighting matrix

Usage

listw.explore()

Value

No return value

Author(s)

Stéphane Dray [email protected]

See Also

chooseCN

Examples

if(interactive()){
## a matrix or an object of class 'Spatial*' should be in the global environment
xy <- matrix(rnorm(50), 25)
listw.explore()
}

---

Function to optimize the selection of a spatial weighting matrix and select the best subset of eigenvectors (MEM, Moran's Eigenvector Maps)

Description

listw.select computes MEM variables (i.e., eigenvectors of a doubly centered spatial weighting matrix) for various definitions of spatial weighting matrices (SWM) and optimizes the selection of the SWM and of a subset of MEM variables. The optimization is done by maximizing the adjusted R-squared (R2) or by minimizing the residual spatial autocorrelation. The function controls the type I error rate by accounting for the number of tests performed. This function combine calls to the functions scores.listw and mem.select. The list of candidate SWMs can easily be generated using listw.candidates.

Usage

listw.select(
  x,
  candidates,
  MEM.autocor = c("positive", "negative", "all"),
  method = c("FWD", "MIR", "global"),
  MEM.all = FALSE,
  nperm = 999,
  nperm.global = 9999,
  alpha = 0.05,
  p.adjust = TRUE,
  verbose = FALSE
)

Arguments

x	Vector, matrix, or dataframe of the response variable(s)
candidates	A list of SWMs of the class listw; candidates can be created by listw.candidates
MEM.autocor	Sign of the spatial eigenvectors to generate; "positive", "negative", or "all", for positively, negatively autocorrelated eigenvectors, or both, respectively; default is "positive"
method	Criterion to select the best subset of MEM variables. Either forward (default option), "MIR" (for univariate x only), or "global" (see Details)
MEM.all	A logical indicating if the complete set of MEM variables for the best model should be returned
nperm	Number of permutations to perform the tests in the selection procedure; Default is 999
nperm.global	Number of permutations to perform the tests in the global test; Default is 9999
alpha	Significance threshold value for the tests; Default is 0.05
p.adjust	A logical indicating wheter the p-value of the global test performed on each SWM should be corrected for multiple tests (TRUE) or not (FALSE); default is TRUE
verbose	If 'TRUE' more diagnostics are printed. The default setting is FALSE

Details

While the selection of the SWM is the most critical step of the spatial eigenvector-based methods (Dray et al. 2006), Bauman et al. (2018) showed that optimizing the choice of the SWM led to inflated type I error rates if an explicit control of the number of SWMs tested was not applied. The function listw.select therefore applies a Sidak correction (Sidak 1967) for multiple tests to the p-value of the global test of each SWM (i.e., the model integrating the whole set of spatial predictors). The Sidak correction is computed as: $P_corrected = 1 - (1 - P)^n$, where $n$ is the number of tests performed, $P$ is the observed p-value, and $P_corrected$ is the new p-value after the correction. The p-value is first computed using nperm permutations and then corrected according to the total number of SWMs tested (if p.adjust = TRUE). Although the function can be run without this correction, using the default value is strongly recommended to avoid inflated type I error rates (Bauman et al. 2018).

As a consequence of the p-value correction, the significance threshold decreases as the number of SWMs increases, hence leading to a trade-off between the gain of accuracy and the power loss.

The optimization criterion of the SWM performed by listw.select is either based on the maximization of the significant adjusted R2 of all the generated spatial eigenvectors (also referred to as spatial predictors or MEM variables) (method = "global"), or is based on an optimized subset of eigenvectors (method = "FWD" and "MIR").

If the objective is only to optimize the selection of the SWM, without the intervention of the selection of a subset of predictors within each SWM (method = "global"), then the best SWM is the one maximizing the significant adjusted global R2, that is, the R2 of the model of x against the whole set of generated MEM variables which must be significant for the global test (method = "global").

The optimization of the SWM depends on the choosen method. See mem.select for a description of the situations in which method = "FWD", "MIR", and "global" should be preferred.

If a subset of MEM variables is needed, then the optimization of the subset of spatial predictors guides the optimization of the selection of SWM (method = "FWD" or "MIR"). If method = "FWD", listw.select performs the forward selection on the significant SWMs and selects among these the SWM for which the forward-selected subset of spatial eigenvectors yields the highest adjusted R2. If method = "MIR", listw.select performs the MIR selection on all the significant candidate SWMs, and selects the best SWM as the one with the smallest number of MIR-selected spatial eigenvectors. If two or more SWMs present the same smallest number of predictors, then the selection is made among them on the basis of the residual Moran's I. If MEM.autocor = "all", the optimization criteria described above are applied on the sum of the adjusted R2 or number of selected spatial eigenvectors, for method = "FWD" and "MIR", respectively. If no subset of MEM variable is required, then the optimization of the SWM is based on the maximization of the adjusted R2 of all the generated MEM variables (method = "global").

If MEM.autocor = "all", n-1 MEM variables are generated. In this case, if method = "global" or method = "FWD", the adjusted R2 is computed separately on the MEM associated to positive and negative eigenvalues (hereafter positive and negative MEM variables, respectively), and the SWM yielding the highest sum of the the two significant R2 values is selected. If method = "MIR", the MIR selection is performed separately on the positive and negative MEM variables, and the SWM is selected based on the sum of the number of positive and negative spatial predictors.

Value

listw.select returns a list that contains:

candidates

A data.frame that summarizes the results on all SWMs

best.id

The index and name of the best SWM

best

The results for the best SWM as returned by mem.select

Author(s)

Bauman David ([email protected] or [email protected]) and Stéphane Dray

References

Bauman D., Fortin M-J, Drouet T. and Dray S. (2018) Optimizing the choice of a spatial weighting matrix in eigenvector-based methods. Ecology

Blanchet G., Legendre P. and Borcard D. (2008) Forward selection of explanatory variables. Ecology, 89(9), 2623–2632

Dray S., Legendre P. and Peres-Neto P. R. (2006) Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM). Ecological Modelling, 196, 483–493

Sidak Z. (1967) Rectangular confidence regions for the means of multivariate normal distributions. Journal of the American Statistical Association, 62(318), 626–633

See Also

listw.candidates, mem.select, scores.listw

Examples

if(require(spdep)) {
### Create a grid of 15 x 15:
grid <- expand.grid(x = seq(1, 15, 1), y = seq(1, 15, 1))
### Generate a response variable Y structured at broad scale by linear combination of
### the first three MEM variables to which a normal noise is added:
nb <- cell2nb(nrow = 15, ncol = 15, "queen")
lw <- nb2listw(nb, style = "B")
MEM <- scores.listw(lw, MEM.autocor = "positive")
# Degree of spatial autocorrelation:
intensity <- 0.8
Y_space <- scale(MEM[, 1] + MEM[, 2] + MEM[, 3]) * intensity
Y_noise <- scale(rnorm(n = nrow(MEM), mean = 0, sd = 1)) * (1 - intensity)
Y <- Y_space + Y_noise
### Y is sampled in 100 randomly-chosen sites of the grid:
idx.sample <- sample(c(1:nrow(grid)), 100, replace = FALSE)
xy <- grid[idx.sample, ]
Y_sampled <- Y[idx.sample]
### The function listw.candidates is used to build the spatial weighting matrices that
### we want to test and compare (with the listw.select function). We test a Gabriel's graph,
### a minimum spanning tree, and a distance-based connectivity defined by a threshold
### distance corresponding to the smallest distance keeping all sites connected (i.e.,
### the defaut value of d2; see help of function listw.candidates).
### These connectivity matrices are then either not weighted (binary weighting), or
### weighted by the linearly decreasing function (see help of the function listw.candidates):
candidates <- listw.candidates(coord = xy, nb = c("gab", "mst"), weights = c("binary", "flin"))
### Number of candidate W matrices generated:
nbw <- length(candidates)
### Significance threshold value after p-value correction (Sidak correction):
1 - (1 - 0.05)^(1/nbw)
### Optimization of the selection of the SWM among the candidates generated above,
### using the corrected significance threshold calculated above for the global tests:
W_sel <- listw.select(Y_sampled, candidates, MEM.autocor = "positive", method = "FWD",
                    p.adjust = TRUE, nperm = 299)
### Some characteristics of the best spatial model:
# Best SWM:
W_sel$best.id
# Selected subset of spatial predictor within the best SWM:
W_sel$best$MEM.select
nrow(W_sel$best$summary)
# Corrected p-value of the global test of the best SWM:
W_sel$best$global.test$Pvalue
# Adjusted R2 of the subset of spatial predictors selected within the chosen SWM:
max(W_sel$best$summary$R2Adj)
# p-values of all the tested W matrices:
W_sel$candidates$Pvalue
# Adjusted R2 of the subset of spatial predictors selected for all the significant
# W matrices:
W_sel$candidates$R2Adj.select

# See Appendix S3 of Bauman et al. 2018 for more extensive examples and illustrations.
}

---

Mastigouche Lake network data set

Description

A list containing the XY coordinates of the lakes and a site-by-edge matrix describing how 42 lakes influence each other. The influence is defined by 66 edges.

Usage

data(mastigouche)

Format

A list that includes the centred coordinates of 42 lakes in the Mastigouche reserve and a site-by-edge matrix describing how the 42 lakes are influenced among each other. The influence is defined by 66 edges.

References

Magnan, P., Rodriguez, M.A., Legendre, P., Lacasse, S. (1994) Dietary variation in a freshwater fish species: relative contribution of biotic interactions, abiotic factors, and spatial structure. Canadian Journal of Fisheries and Aquatic Sciences 51, 2856-2865. Blanchet F.G., P. Legendre, and Borcard D. (2008) Modelling directional spatial processes in ecological data. Ecological Modelling, 215, 325-336.

---

Selection of the best subset of spatial eigenvectors (MEM, Moran's Eigenvector Maps)

Description

mem.select computes the spatial eigenvectors (MEM) of the spatial weighting matrix (SWM) provided (listw) and optimizes the selection of a subset of MEM variables relative to response variable(s) stored in x. The optimization is done either by maximizing the adjusted R-squared (R2) of all (method = "global") or a subset (method = "FWD") of MEM variables or by minimizing the residual spatial autocorrelation (method = "MIR") (see details in Bauman et al. 2018a).

Usage

mem.select(
  x,
  listw,
  MEM.autocor = c("positive", "negative", "all"),
  method = c("FWD", "MIR", "global"),
  MEM.all = FALSE,
  nperm = 999,
  nperm.global = 9999,
  alpha = 0.05,
  verbose = FALSE,
  ...
)

Arguments

x	A vector, matrix, or dataframe of response variable(s). The method = "MIR" is only implemented for a vector response. Note that x can also contain the residuals of a model when other (e.g., environmental) variables should be considered in the model.
listw	A spatial weighting matrix of class listw; can be created with functions of the package spdep, or with the user-friendly function listw.candidates. Note that, the function listw.candidates returns a list of listw and subselection by [[]] should be performed in this case (see Example)
MEM.autocor	Sign of the spatial eigenvectors to generate; "positive", "negative", or "all", for positively, negatively autocorrelated eigenvectors, or both, respectively; default is "positive"
method	Criterion to select the best subset of MEM variables. Either "FWD" (default option), "MIR" (for univariate x only), or "global" (see Details)
MEM.all	A logical indicating if the complete set of MEM variables should be returned
nperm	Number of permutations to perform the tests in the selection procedure; Default is 999
nperm.global	Number of permutations to perform the tests in the global test; Default is 9999
alpha	Significance threshold value for the tests; Default is 0.05
verbose	If 'TRUE' more diagnostics are printed. The default setting is FALSE
...	Other parameters (for internal use with listw.select)

Details

The function provides three different methods to select a subset of MEM variables. For all methods, a global test is firstly performed. If MEM.autocor = "all", two global tests are performed and p-values are corrected for multiple comparison (Sidak correction).

If the MEM variables are to be further used in a model including actual predictors (e.g. environmental), then a subset of spatial eigenvectors needs to be selected before proceeding to further analyses to avoid model overfitting and/or a loss of statistical power to detect the contribution of the environment to the variability of the response data (Griffith 2003, Dray et al. 2006, Blanchet et al. 2008, Peres-Neto and Legendre 2010, Diniz-Filho et al. 2012). Although several eigenvector selection approaches have been proposed to select a best subset of eigenvectors, Bauman et al. (2018b) showed that two main procedures should be preferred, depending on the underlying objective: the forward selection with double stopping criterion (Blanchet et al. 2008; method = "FWD") or the minimization of the residual spatial autocorrelation (Griffith and Peres-Neto 2006; MIR selection in Bauman et al. 2018a,b, method = "MIR"). The most powerful and accurate selection method, in terms of R2 estimation, is the forward selection. This method should be preferred when the objective is to capture as accurately as possible the spatial patterns of x. If the objective is to optimize the detection of the spatial patterns in the residuals of a model of the response variable(s) against a set of environmental predictors, for instance, then x should be the model residuals, and method = "FWD". This allows optimizing the detection of residual spatial patterns once the effect of the environmental predictors has been removed. If however the objective is only to remove the spatial autocorrelation from the residuals of a model of x against a set of actual predictors (e.g. environmental) with a small number of spatial predictors, then accuracy is not as important and one should focus mainly on the number of spatial predictors (Bauman et al. 2018b). In this case, method = "MIR" is more adapted, as it has the advantage to maintain the standard errors of the actual predictor coefficients as low as possible. Note that method = "MIR" can only be used for a univariate x, as the Moran's I is a univariate index. If x is multivariate, then the best criterion is the forward selection (see Bauman et al. 2018b). A third option is to not perform any selection of MEM variables (method = "global"). This option may be interesting when the complete set of MEM variables will be used, like in Moran spectral randomizations (Wagner and Dray 2015, Bauman et al. 2018c) or when using smoothed MEM (Munoz 2009).

For method = "MIR", the global test consists in computing the Moran's I of x (e.g. residuals of the model of the response variable against environmental variables) and tests it by permutation (results stored in global.test). If the Moran's I is significant, the function performs a selection procedure that searches among the set of generated spatial predictors the one that best minimizes the value of the Moran's I. A model of x against the selected eigenvector is built, and the significance of the Moran's I of the model residuals is tested again. The procedure goes on until the Moran's I of the model residuals is not significant anymore, hence the name of Minimization of moran's I in the Residuals (MIR).

For method = "global" and method = "FWD", the global test consists in computing the adjusted global R2, that is, the R2 of the model of x against the whole set of generated MEM variables and tests it by permutation (results stored in global.test).

For method = "global", if the adjusted global R2 is significant, the functions returns the whole set of generated MEM variables in MEM.select.

For method = "FWD", if the adjusted global R2 is significant, the function performs a forward selection with double stopping criterion that searches among the set of generated spatial predictors the one that best maximizes the R2 of the model. The procedure is repeated untill one of the two stopping criterion is reached (see Blanchet et al. 2008). Note that in a few cases, the forward selection does not select any variable even though the global model is significant. This can happen for example when a single variable has a strong relation with the response variable(s), because the integration of the variable alone yields an adjusted R2 slightly higher than the global adjusted R2. In this case, we recommend checking that this is indeed the reason why the first selected variable was rejected, and rerun the analysis with a second stopping criterion equal to the global adjusted R2 plus a small amount allowing avoiding this issue (e.g. 5 done through the argument adjR2thresh of function forward.sel, until the solution is implemented in mem.select.

For the method = "FWD" and method = "MIR", the MEM selected by the procedure are returned in MEM.select and a summary of the results is provided in summary. If no MEM are selected, then MEM.select and summary are not returned.

Value

The function returns a list with:

global.test

An object of class randtest containing the result of the global test associated to all MEM (adjusted R2 and p-value). If MEM.autocor = "all", a list with two elements (positive and negative) corresponding to the results of the global tests performed on positive and negative MEM respectively.

MEM.all

An object of class orthobasisSp containing the complete set of generated MEM variables (generated by scores.listw). Only returned if MEM.all = TRUE.

summary

A dataframe summarizing the results of the selection procedure

MEM.select

An object of class orthobasisSp containing the subset of significant MEM variables.

Author(s)

David Bauman ([email protected] or [email protected]) and Stéphane Dray

References

Bauman D., Fortin M-J, Drouet T. and Dray S. (2018a) Optimizing the choice of a spatial weighting matrix in eigenvector-based methods. Ecology, 99, 2159-2166

Bauman D., Drouet T., Dray S. and Vleminckx J. (2018b) Disentangling good from bad practices in the selection of spatial or phylogenetic eigenvectors. Ecography, 41, 1–12

Bauman D., Vleminckx J., Hardy O., Drouet T. (2018c) Testing and interpreting the shared space-environment fraction in variation partitioning analyses of ecological data. Oikos

Blanchet G., Legendre P. and Borcard D. (2008) Forward selection of explanatory variables. Ecology, 89(9), 2623–2632

Diniz-Filho J.A.F., Bini L.M., Rangel T.F., Morales-Castilla I. et al. (2012) On the selection of phylogenetic eigenvectors for ecological analyses. Ecography, 35, 239–249

Dray S., Legendre P. and Peres-Neto P. R. (2006) Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM). Ecological Modelling, 196, 483–493

Griffith D. (2003) Spatial autocorrelation and spatial filtering: gaining understanding through theory and scientific visualization. Springer, Berlin

Griffith D. and Peres-Neto P. (2006) Spatial modeling in Ecology: the flexibility of eigenfunction spatial analyses. Ecology, 87, 2603–2613

Munoz, F. 2009. Distance-based eigenvector maps (DBEM) to analyse metapopulation structure with irregular sampling. Ecological Modelling, 220, 2683–2689

Peres-Neto P. and Legendre P. (2010) Estimating and controlling for spatial structure in the study of ecological communities. Global Ecology and Biogeography, 19, 174–184

Wagner H., Dray S. (2015). Generating spatially constrained null models for irregularly spaced data using Moran spectral randomization methods. Methods in Ecology and Evolution, 6, 1169–1178

See Also

listw.candidates, listw.select, scores.listw

Examples

if(require(vegan)){ 
# Illustration of the MIR selection on the oribatid mite data
# (Borcard et al. 1992, 1994 for details on the dataset):
# *******************************************************
# Community data (response matrix):
data(mite)
# We will compute the example on a single species:
spe <- mite[, 2]
# Environmental explanatory dataset:
data(mite.env)
# We only use two numerical explanatory variables:
env <- mite.env[, 1:2]
dim(env)
# Coordinates of the 70 sites:
data(mite.xy)
coord <- mite.xy
# We build the model we are interested in:
mod <- lm(spe ~ ., data = env)


# In order to avoid possible type I error rate inflation issues, we check 
# whether the model residuals are independent, and if they are spatially
# autocorrelated, we select a small subset of MEM variables to add to the
# model as covariables with the MIR selection:

# 1) We build a spatial weighting matrix based on Gabriel graph with a
# weighting function decreasing linearly with the distance:
w <- listw.candidates(coord, nb = "gab", weights = "flin")


# 2) We test the spatial autocorrelation of the model residuals and, if
# necessary, select a subset of spatial predictors:
y <- residuals(mod)
MEM <- mem.select(x = y, listw = w[[1]], method = "MIR", MEM.autocor = "positive",
         nperm = 999, alpha = 0.05)
dim(MEM$MEM.select)
# The residuals of the model presented spatial autocorrelation. The selection
# of MEM variables is thus performed to remove residual autocorrelation.

# 3) We can reconstruct our model adding the selected MEM variable as covariables:
env2 <- cbind(env, MEM$MEM.select)
mod_complete <- lm(spe ~ ., data = env2)
summary(mod_complete)$coefficient[, 1]   # Coefficient estimates
summary(mod_complete)$coefficient[, 2]   # Standard errors
}

---

Multi-frequential periodogram analysis

Description

This function performs multi-frequential periodogram analysis for univariate temporal or spatial data series collected with equal intervals. Compared with the traditional periodogram used in spectral analysis, this method can detect overlapping signals with fractional frequencies. Fitting a joint polynomial-trigonometric model is achieved by Ordinary Least Squares (OLS) regression. The function also performs autocorrelation analysis of OLS residuals up to a number of lags determined by the user.

Usage

mfpa(
  y,
  MaxNFreq = 2,
  MinFreq = 3,
  MaxFreq = NA,
  ntrend = 0,
  nlags = 0,
  alpha = 0.05
)

## S3 method for class 'mfpa'
plot(x, xlab = "", ylab = "Values", ...)

## S3 method for class 'mfpa'
print(x, ...)

Arguments

y	Vector of n observations (vector of integer or real numbers, or one-column matrix).
MaxNFreq	Maximum number of frequencies to be estimated in the stepwise procedure (e.g. 2).
MinFreq	Minimum value for frequency estimates (e.g. 3.0).
MaxFreq	Maximum value for frequency estimates (e.g. 10). Must be larger than MinFreq and smaller than half of the number of observations in the series. If unspecified by the user (MaxFreq=NA), MaxFreq is set to n/4 by the function.
ntrend	Number (0 to 3) of orthogonal polynomial components estimating the broad-scale trend, to be included in the joint polynomial-trigonometric model. Use 0 to estimate no trend component, only an intercept.
nlags	Number of lags to be used for autocorrelation analysis of OLS residuals. Use 0 to bypass this analysis.
alpha	Significance threshold for including frequencies.
x	An object of class mfpa
xlab, ylab	Labels for x and y axes
...	Further arguments passed to or from other methods

Details

The fitting of a joint polynomial-trigonometric model is limited to ordinary least squares (OLS), with autocorrelation analysis of OLS residuals up to a certain lag. Orthogonal polynomials are used to model broad-scale trends, whereas cosines and sines model the periodic structures at intermediate scales. See Dutilleul (2011, section 6.5) and Legendre & Legendre (2012, section 12.4.4) for details. OLS regression could be replaced by an estimated generalized least squares (EGLS) procedure, as described in Dutilleul (2011).

In spectral analysis in general and in mfpa in particular, the cosines and sines are considered jointly in the search for the dominant frequency components since they are both required to fully account for a frequency component in a linear model. So, when either the cosine or the sine is significant, this is sufficient indication that a significant frequency component has been found. But see the first paragraph of the 'Recommendations to users' below.

The periodic phenomenon corresponding to each identified frequency is modelled by a cosine and a sine. The first pair ('cos 1', 'sin 1') corresponds to the first frequency, the second pair to the second frequency, and so on. An intercept is also computed, as well as a polynomial broad-scale trend if argument ntrend > 0. The coefficients shown for each periodic component ('cos' and 'sin') are the OLS regression coefficients. The tests of significance producing the p-values (called 'prob' in the output file) are 2-tailed parametric t-tests, as in standard OLS regression.

A global R-square statistic for the periodogram is computed as the variance of the fitted values divided by the variance of the data series. An R-squared corresponding to each frequency is also returned.

In the Dutilleul periodogram, the time unit is the length of the data series (in time units: seconds, hours, days, etc.). Hence, the frequency identified by a Dutilleul periodogram is the number of cycles of the periodic signal (how many full or partial cycles) along the time series. That number is an integer when the series contains an integer number of cycles; it may also be a real number when the number of cycles is fractional. The periodogram can identify several periodic phenomena with different frequencies. The estimated frequencies could be divided by an appropriate constant to produce numbers of cycles per second or day, or per meter or km, depending on the study.

To find the period (number of days, hours, etc.) of the process generating a periodic signal in the data, divide the length of the series (in days, hours, etc.) by the frequency identified by Dutilleul's periodogram.

Recommendations to users The mfpa code estimates the periodic frequencies to be included in the model through a combination of a stepwise procedure and non-linear optimisation. Following that, the contributions of the 'cos' and 'sin' components of all frequencies in the model are estimated by multiple linear regression in the presence of the intercept and trends (if any). Because the mfpa method estimates fractional frequencies, the cos-sin combinations are not orthogonal among the identified frequencies, and unnecessary frequencies may be selected as 'significant'.

1. It is important that users of this periodogram have hypotheses in mind about the frequencies of the processes that may be operating on the system under study and the number of periodic components they are expecting to find. If one asks for more components than the number of periodic phenomena at work on the system, the 'real' frequency usually has a strong or fairly strong R-squared and it is followed by other components with very small R-squared. Selection of frequencies of interest should thus be based more upon examination of the R-squares of the components rather than on the p-values. For short series in particular, the adjusted R-squared is an unbiased estimate of the variance of the data explained by the model. Even series of random numbers can produce 'significant' frequencies for periodic components; the associated (adjusted) R-squares will, however, be very small.

2. Function mfpa cannot detect frequencies < 1 (smaller than one cycle in the series) or larger than (n-1) where n is the number of observations in the series, the latter case corresponding to periods smaller than the interval between successive observations. When a periodic component with such a period is present in the data, Dutilleul's periodogram can detect harmonics of that frequency. Recommendation: when a frequency is detected that does not seem to correspond to a hypothesized process, one could check, using simulated data, if it could be produced by a process operating at a temporal scale (period) smaller than the interval between successive observations. An example is shown in Example 2.

3. When analysing a time series with unknown periodic structure, it is recommended to try first with more than one frequency, say 2 or 3, and also with a trend. Eliminate the non-significant components, step by step, in successive runs, starting with the trend(s), then eliminate the weakly significant periodic components, until there are only highly significant components left in the model.

Value

A list containing the following elements:

• frequencies: Vector of estimated frequencies of the model periodic components and associated R-squared. The frequencies are numbers of cycles in the whole (temporal or spatial) series under study.

• coefficients: Data frame containing OLS slope estimates, starting with the intercept, then the orthogonal polynomials modelling trend in increasing order, followed by the cosine and sine coefficients (alternating) in the order of the estimated frequencies. Columns: (1) coefficient: the OLS intercept or slope estimates; (2) prob: the associated probabilities.

• predicted: A vector (length n) of predicted response values (fitted values), including the trend if any. The data and predicted values can be plotted together using function plot.mfpa; type plot(name.of.output.object). The data values are represented by red circles and the fitted values by a black line.

• auto_coeff: If nlags > 0: data frame containing the following columns. (1) lag: lags at which autocorrelation analysis of the OLS residuals is performed; (2) auto_r: vector of sample autocorrelation coefficients calculated from OLS residuals for each lag; (3) prob: vector of probabilities associated with the tests of significance of the sample autocorrelation coefficients.

• y: the original data series (one-column matrix).

• X: the matrix of explanatory variables; it contains a column of "1" to estimate the intercept, a column for each of the trend components (if any), and two columns for each frequency component, each frequency being represented by a cosine and a sine.

• r.squared.global: The global R-squared of the model and the adjusted R-squared.

Author(s)

Guillaume Larocque <[email protected]> and Pierre Legendre.

References

Dutilleul, P. 1990. Apport en analyse spectrale d'un périodogramme modifié et modélisation des séries chronologiques avec répétitions en vue de leur comparaison en fréquence. Doctoral Dissertation, Université Catholique de Louvain, Louvain-la-Neuve, Belgium.

Dutilleul, P. 1998. Incorporating scale in ecological experiments: data analysis. Pp. 387-425 in: D. L. Peterson & V. T. Parker [eds.] Ecological scale - Theory and applications. Columbia University Press, New York.

Dutilleul, P. 2001. Multi-frequential periodogram analysis and the detection of periodic components in time series. Commun. Stat. - Theory Methods 30, 1063-1098.

Dutilleul, P. R. L. 2011. Spatio-temporal heterogeneity - Concepts and analyses. Cambridge University Press, Cambridge.

Dutilleul, P. and C. Till. 1992. Evidence of periodicities related to climate and planetary behaviors in ring-width chronologies of Atlas cedar (Cedrus atlantica) in Morocco. Can. J. For. Res. 22: 1469-1482.

Legendre, P. and P. Dutilleul. 1992. Introduction to the analysis of periodic phenomena. 11-25 in: M. A. Ali [ed.] Rhythms in fishes. NATO ASI Series, Vol. A-236. Plenum, New York.

Legendre, P. and L. Legendre. 2012. Numerical Ecology. 3rd English edition. Elsevier, Amsterdam.

Examples

### Example 1

# Simulate data with frequencies 2.3 and 6.1 and a random component, n = 100. 
# No trend, no autocorrelated residuals.

y <- as.matrix(0.4*(sin(2.3*2*pi*(1:100)/100)) +
0.4*(sin(6.1*2*pi*(1:100)/100)) + 0.2*rnorm(100))

res <- mfpa(y, MaxNFreq = 2, MinFreq = 2, ntrend = 0, nlags = 0)

# Compute the periods associated with the two periodic components. Each
# frequency in element $frequencies is a number of cycles in the whole series.
# The periods are expressed in numbers of time intervals of the data series. In
# this example, if the data are measured every min, the periods are in min.

periods <- 100/res$frequencies$frequency 

# Draw the data series and the fitted (or predicted) values

plot(res)

### Example 2

# Generate hourly periodic data with tide signal (tide period T = 12.42 h)
# during 1 year, hence 24*365 = 8760 hourly data. See
# https://en.wikipedia.org/wiki/Tide.

# In this simulation, constant (c = 0) puts the maximum value of the cosine at
# midnight on the first day of the series.

periodic.component <- function(x, T, c) cos((2*pi/T)*(x+c))

tide.h <- periodic.component(1:8760, 12.42, 0)

# The number of tides in the series is: 8760/12.42 = 705.314 tidal cycles
# during one year.

# Sample the hourly data series once a day at 12:00 noon every day. The
# periodic signal to be detected has a period smaller then the interval between
# consecutive observations and its frequency is larger than (n-1). The sequence
# of sampling hours for the tide.h data is:

h.noon <- seq(12, 8760, 24) 
tide.data <- tide.h[h.noon]
length(tide.data)   

# The series contains 365 sampling units

# Compute Dutilleul's multi-frequential periodogram

res.noon <- mfpa(tide.data, MaxNFreq = 1, MinFreq = 2, ntrend = 1, nlags = 2)

# Examine the frequency detected by the periodogram, element
# res.noon$frequencies. This is a harmonic of the tide signal in the original
# data series tide.h.

# Compute the period of the signal in the data series sampled hourly:

period <- 365/res.noon$frequencies$frequency

# Draw the data series and the adjusted values

plot(res.noon)

# Repeat this analysis after addition of random noise to the tide data

tide.noise <- tide.data + rnorm(365, 0, 0.25)

res.noise <- mfpa(tide.noise, MaxNFreq = 1, MinFreq = 2, ntrend = 1, nlags = 2)

plot(res.noise)

---

Function to compute extreme values of Moran's I

Description

This function computes the upper and lower bounds of Moran's I for a given spatial weighting matrix (stored in a listw object). These values are obtained by the eigendecomposition of the spatial weighting matrix.

Usage

moran.bounds(listw)

Arguments

listw

an object of class listw

Value

A vector containing the maximum and minimum of Moran's I for a given spatial weighting matrix value returned

Author(s)

Stéphane Dray [email protected]

References

de Jong, P., Sprenger, C., & van Veen, F. (1984). On extreme values of Moran's I and Geary's C. Geographical Analysis, 16(1), 17-24.

See Also

mem nb2listw

Examples

if(require("ade4", quietly = TRUE)){
 if(require("spdep", quietly = TRUE)){
     data(oribatid)
     nbtri <- tri2nb(as.matrix(oribatid$xy))
     lwB <- nb2listw(nbtri, style = "B")
     lwW <- nb2listw(nbtri, style = "W")
     scB <- mem(lwB)
     scW <- mem(lwW)
     moran.bounds(lwB)
     moran.mc(scB[,1], lwB, 9)
     moran.mc(scB[,69], lwB, 9)
     moran.bounds(lwW)
     moran.mc(scW[,1], lwW, 9)
     moran.mc(scW[,69], lwW, 9)
 }
}

---

Function to compute Moran's index of spatial autocorrelation

Description

This function computes Moran's I statistic and provide a testing procedure using random permutations. It is based on the moran.mc function of the spdep package. The moran.randtest is slightly different as it allows to consider several variables (x can have more than one columns) and its ouputs are objects of class randtest (one variable) or krandtest (several variables).

Usage

moran.randtest(x, listw, nrepet = 999, ...)

Arguments

x	a vector, matrix or data.frame with numeric data
listw	an object of class listw created for example by nb2listw
nrepet	an integer indicating the number of permutations used in the randomization procedure
...	other arguments to be passed to the as.randtest or as.krandtest functions.

Value

An object of class randtest (one variable) or krandtest (several variables)

Author(s)

Stéphane Dray [email protected]

References

Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37, 17-23.

See Also

moran.mc

Examples

if(require("ade4", quietly = TRUE)  & require("spdep", quiet = TRUE)){
data(mafragh)
tests <- moran.randtest(mafragh$env, nb2listw(mafragh$nb))
tests
plot(tests)

}

---