Package 'vegan'

Community Ecology Package: Ordination, Diversity and Dissimilarities

Description

The vegan package provides tools for descriptive community ecology. It has most basic functions of diversity analysis, community ordination and dissimilarity analysis. Most of its multivariate tools can be used for other data types as well.

Details

The functions in the vegan package contain tools for diversity analysis, ordination methods and tools for the analysis of dissimilarities. Together with the labdsv package, the vegan package provides most standard tools of descriptive community analysis. Package ade4 provides an alternative comprehensive package, and several other packages complement vegan and provide tools for deeper analysis in specific fields. Package BiodiversityR provides a GUI for a large subset of vegan functionality.

The vegan package is developed at GitHub (https://github.com/vegandevs/vegan/). GitHub provides up-to-date information and forums for bug reports.

Most important changes in vegan documents can be read with news(package="vegan") and vignettes can be browsed with browseVignettes("vegan"). The vignettes include a vegan FAQ, discussion on design decisions, short introduction to ordination and discussion on diversity methods.

To see the preferable citation of the package, type citation("vegan").

Author(s)

The vegan development team is Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Helene Wagner. Many other people have contributed to individual functions: see credits in function help pages.

Examples

### Example 1: Unconstrained ordination
## NMDS
data(varespec)
data(varechem)
ord <- metaMDS(varespec)
plot(ord, type = "t")
## Fit environmental variables
ef <- envfit(ord, varechem)
ef
plot(ef, p.max = 0.05)
### Example 2: Constrained ordination (RDA)
## The example uses formula interface to define the model
data(dune)
data(dune.env)
## No constraints: PCA
mod0 <- rda(dune ~ 1, dune.env)
mod0
plot(mod0)
## All environmental variables: Full model
mod1 <- rda(dune ~ ., dune.env)
mod1
plot(mod1)
## Automatic selection of variables by permutation P-values
mod <- ordistep(mod0, scope=formula(mod1))
mod
plot(mod)
## Permutation test for all variables
anova(mod)
## Permutation test of "type III" effects, or significance when a term
## is added to the model after all other terms
anova(mod, by = "margin")
## Plot only sample plots, use different symbols and draw SD ellipses 
## for Managemenet classes
plot(mod, display = "sites", type = "n")
with(dune.env, points(mod, disp = "si", pch = as.numeric(Management)))
with(dune.env, legend("topleft", levels(Management), pch = 1:4,
  title = "Management"))
with(dune.env, ordiellipse(mod, Management, label = TRUE))
## add fitted surface of diversity to the model
ordisurf(mod, diversity(dune), add = TRUE)
### Example 3: analysis of dissimilarites a.k.a. non-parametric
### permutational anova
adonis2(dune ~ ., dune.env)
adonis2(dune ~ Management + Moisture, dune.env)

---

Add or Drop Single Terms to a Constrained Ordination Model

Description

Compute all single terms that can be added to or dropped from a constrained ordination model.

Usage

## S3 method for class 'cca'
add1(object, scope, test = c("none", "permutation"),
    permutations = how(nperm=199), ...)
## S3 method for class 'cca'
drop1(object, scope, test = c("none", "permutation"), 
    permutations = how(nperm=199), ...)

Arguments

object	A constrained ordination object from cca, rda, dbrda or capscale.
scope	A formula giving the terms to be considered for adding or dropping; see add1 for details.
test	Should a permutation test be added using anova.cca.
permutations	a list of control values for the permutations as returned by the function how, or the number of permutations required, or a permutation matrix where each row gives the permuted indices.
...	Other arguments passed to add1.default, drop1.default, and anova.cca.

Details

With argument test = "none" the functions will only call add1.default or drop1.default. With argument test = "permutation" the functions will add test results from anova.cca. Function drop1.cca will call anova.cca with argument by = "margin". Function add1.cca will implement a test for single term additions that is not directly available in anova.cca.

Functions are used implicitly in step, ordiR2step and ordistep. The deviance.cca and deviance.rda used in step have no firm basis, and setting argument test = "permutation" may help in getting useful insight into validity of model building. Function ordistep calls alternately drop1.cca and add1.cca with argument test = "permutation" and selects variables by their permutation $P$-values. Meticulous use of add1.cca and drop1.cca will allow more judicious model building.

The default number of permutations is set to a low value, because permutation tests can take a long time. It should be sufficient to give a impression on the significances of the terms, but higher values of permutations should be used if $P$ values really are important.

Value

Returns a similar object as add1 and drop1.

Author(s)

Jari Oksanen

See Also

add1, drop1 and anova.cca for basic methods. You probably need these functions with step and ordistep. Functions deviance.cca and extractAIC.cca are used to produce the other arguments than test results in the output.

Examples

data(dune)
data(dune.env)
## Automatic model building based on AIC but with permutation tests
step(cca(dune ~  1, dune.env), reformulate(names(dune.env)), test="perm")
## see ?ordistep to do the same, but based on permutation P-values
## Not run: 
ordistep(cca(dune ~  1, dune.env), reformulate(names(dune.env)))

## End(Not run)
## Manual model building
## -- define the maximal model for scope
mbig <- rda(dune ~  ., dune.env)
## -- define an empty model to start with
m0 <- rda(dune ~ 1, dune.env)
## -- manual selection and updating
add1(m0, scope=formula(mbig), test="perm")
m0 <- update(m0, . ~ . + Management)
add1(m0, scope=formula(mbig), test="perm")
m0 <- update(m0, . ~ . + Moisture)
## -- included variables still significant?
drop1(m0, test="perm")
add1(m0, scope=formula(mbig), test="perm")

---

Additive Diversity Partitioning and Hierarchical Null Model Testing

Description

In additive diversity partitioning, mean values of alpha diversity at lower levels of a sampling hierarchy are compared to the total diversity in the entire data set (gamma diversity). In hierarchical null model testing, a statistic returned by a function is evaluated according to a nested hierarchical sampling design (hiersimu).

Usage

adipart(...)
## Default S3 method:
adipart(y, x, index=c("richness", "shannon", "simpson"),
    weights=c("unif", "prop"), relative = FALSE, nsimul=99,
    method = "r2dtable", ...)
## S3 method for class 'formula'
adipart(formula, data, index=c("richness", "shannon", "simpson"),
    weights=c("unif", "prop"), relative = FALSE, nsimul=99,
    method = "r2dtable", ...)

hiersimu(...)
## Default S3 method:
hiersimu(y, x, FUN, location = c("mean", "median"),
    relative = FALSE, drop.highest = FALSE, nsimul=99,
    method = "r2dtable", ...)
## S3 method for class 'formula'
hiersimu(formula, data, FUN, location = c("mean", "median"),
    relative = FALSE, drop.highest = FALSE, nsimul=99,
    method = "r2dtable", ...)

Arguments

y	A community matrix.
x	A matrix with same number of rows as in y, columns coding the levels of sampling hierarchy. The number of groups within the hierarchy must decrease from left to right. If x is missing, function performs an overall decomposition into alpha, beta and gamma diversities.
formula	A two sided model formula in the form y ~ x, where y is the community data matrix with samples as rows and species as column. Right hand side (x) must be grouping variables referring to levels of sampling hierarchy, terms from right to left will be treated as nested (first column is the lowest, last is the highest level). The formula will add a unique indentifier to rows and constant for the rows to always produce estimates of row-level alpha and overall gamma diversities. You must use non-formula interface to avoid this behaviour. Interaction terms are not allowed.
data	A data frame where to look for variables defined in the right hand side of formula. If missing, variables are looked in the global environment.
index	Character, the diversity index to be calculated (see Details).
weights	Character, "unif" for uniform weights, "prop" for weighting proportional to sample abundances to use in weighted averaging of individual alpha values within strata of a given level of the sampling hierarchy.
relative	Logical, if TRUE then alpha and beta diversity values are given relative to the value of gamma for function adipart.
nsimul	Number of permutations to use. If nsimul = 0, only the FUN argument is evaluated. It is thus possible to reuse the statistic values without a null model.
method	Null model method: either a name (character string) of a method defined in make.commsim or a commsim function. The default "r2dtable" keeps row sums and column sums fixed. See oecosimu for Details and Examples.
FUN	A function to be used by hiersimu. This must be fully specified, because currently other arguments cannot be passed to this function via ....
location	Character, identifies which function (mean or median) is to be used to calculate location of the samples.
drop.highest	Logical, to drop the highest level or not. When FUN evaluates only arrays with at least 2 dimensions, highest level should be dropped, or not selected at all.
...	Other arguments passed to functions, e.g. base of logarithm for Shannon diversity, or method, thin or burnin arguments for oecosimu.

Details

Additive diversity partitioning means that mean alpha and beta diversities add up to gamma diversity, thus beta diversity is measured in the same dimensions as alpha and gamma (Lande 1996). This additive procedure is then extended across multiple scales in a hierarchical sampling design with $i = 1, 2, 3, \ldots, m$ levels of sampling (Crist et al. 2003). Samples in lower hierarchical levels are nested within higher level units, thus from $i=1$ to $i=m$ grain size is increasing under constant survey extent. At each level $i$, $\alpha_i$ denotes average diversity found within samples.

At the highest sampling level, the diversity components are calculated as

$\beta_m = \gamma - \alpha_m$

For each lower sampling level as

$\beta_i = \alpha_{i+1} - \alpha_i$

Then, the additive partition of diversity is

$\gamma = \alpha_1 + \sum_{i=1}^m \beta_i$

Average alpha components can be weighted uniformly (weight="unif") to calculate it as simple average, or proportionally to sample abundances (weight="prop") to calculate it as weighted average as follows

$\alpha_i = \sum_{j=1}^{n_i} D_{ij} w_{ij}$

where $D_{ij}$ is the diversity index and $w_{ij}$ is the weight calculated for the $j$th sample at the $i$th sampling level.

The implementation of additive diversity partitioning in adipart follows Crist et al. 2003. It is based on species richness ($S$, not $S-1$), Shannon's and Simpson's diversity indices stated as the index argument.

The expected diversity components are calculated nsimul times by individual based randomisation of the community data matrix. This is done by the "r2dtable" method in oecosimu by default.

hiersimu works almost in the same way as adipart, but without comparing the actual statistic values returned by FUN to the highest possible value (cf. gamma diversity). This is so, because in most of the cases, it is difficult to ensure additive properties of the mean statistic values along the hierarchy.

Value

An object of class "adipart" or "hiersimu" with same structure as oecosimu objects.

Author(s)

Péter Sólymos, [email protected]

References

Crist, T.O., Veech, J.A., Gering, J.C. and Summerville, K.S. (2003). Partitioning species diversity across landscapes and regions: a hierarchical analysis of $\alpha$, $\beta$, and $\gamma$-diversity. Am. Nat., 162, 734–743.

Lande, R. (1996). Statistics and partitioning of species diversity, and similarity among multiple communities. Oikos, 76, 5–13.

See Also

See oecosimu for permutation settings and calculating $p$-values. multipart for multiplicative diversity partitioning.

Examples

## NOTE: 'nsimul' argument usually needs to be >= 99
## here much lower value is used for demonstration

data(mite)
data(mite.xy)
data(mite.env)
## Function to get equal area partitions of the mite data
cutter <- function (x, cut = seq(0, 10, by = 2.5)) {
    out <- rep(1, length(x))
    for (i in 2:(length(cut) - 1))
        out[which(x > cut[i] & x <= cut[(i + 1)])] <- i
    return(out)}
## The hierarchy of sample aggregation
levsm <- with(mite.xy, data.frame(
    l1=1:nrow(mite),
    l2=cutter(y, cut = seq(0, 10, by = 2.5)),
    l3=cutter(y, cut = seq(0, 10, by = 5)),
    l4=rep(1, nrow(mite))))
## Let's see in a map
par(mfrow=c(1,3))
plot(mite.xy, main="l1", col=as.numeric(levsm$l1)+1, asp = 1)
plot(mite.xy, main="l2", col=as.numeric(levsm$l2)+1, asp = 1)
plot(mite.xy, main="l3", col=as.numeric(levsm$l3)+1, asp = 1)
par(mfrow=c(1,1))
## Additive diversity partitioning
adipart(mite, index="richness", nsimul=19)
## the next two define identical models
adipart(mite, levsm, index="richness", nsimul=19)
adipart(mite ~ l2 + l3, levsm, index="richness", nsimul=19)
## Hierarchical null model testing
## diversity analysis (similar to adipart)
hiersimu(mite, FUN=diversity, relative=TRUE, nsimul=19)
hiersimu(mite ~ l2 + l3, levsm, FUN=diversity, relative=TRUE, nsimul=19)
## Hierarchical testing with the Morisita index
morfun <- function(x) dispindmorisita(x)$imst
hiersimu(mite ~., levsm, morfun, drop.highest=TRUE, nsimul=19)

---

Permutational Multivariate Analysis of Variance Using Distance Matrices

Description

Analysis of variance using distance matrices — for partitioning distance matrices among sources of variation and fitting linear models (e.g., factors, polynomial regression) to distance matrices; uses a permutation test with pseudo-$F$ ratios.

Usage

adonis2(formula, data, permutations = 999, method = "bray",
    sqrt.dist = FALSE, add = FALSE, by = NULL,
    parallel = getOption("mc.cores"), na.action = na.fail,
    strata = NULL, ...)

Arguments

formula	Model formula. The left-hand side (LHS) of the formula must be either a community data matrix or a dissimilarity matrix, e.g., from vegdist or dist. If the LHS is a data matrix, function vegdist will be used to find the dissimilarities. The right-hand side (RHS) of the formula defines the independent variables. These can be continuous variables or factors, they can be transformed within the formula, and they can have interactions as in a typical formula.
data	the data frame for the independent variables, with rows in the same order as the community data matrix or dissimilarity matrix named on the LHS of formula.
permutations	a list of control values for the permutations as returned by the function how, or the number of permutations required, or a permutation matrix where each row gives the permuted indices.
method	the name of any method used in vegdist to calculate pairwise distances if the left hand side of the formula was a data frame or a matrix.
sqrt.dist	Take square root of dissimilarities. This often euclidifies dissimilarities.
add	Add a constant to the non-diagonal dissimilarities such that all eigenvalues are non-negative in the underlying Principal Co-ordinates Analysis (see wcmdscale for details). Choice "lingoes" (or TRUE) use the recommended method of Legendre & Anderson (1999: “method 1”) and "cailliez" uses their “method 2”.
by	by = NULL will assess the overall significance of all terms together, by = "terms" will assess significance for each term (sequentially from first to last), setting by = "margin" will assess the marginal effects of the terms (each marginal term analysed in a model with all other variables), by = "onedf" will analyse one-degree-of-freedom contrasts sequentially. The argument is passed on to anova.cca.
parallel	Number of parallel processes or a predefined socket cluster. With parallel = 1 uses ordinary, non-parallel processing. The parallel processing is done with parallel package.
na.action	Handling of missing values on the right-hand-side of the formula (see na.fail for explanation and alternatives). Missing values are not allowed on the left-hand-side. NB, argument subset is not implemented.
strata	Groups within which to constrain permutations. The traditional non-movable strata are set as Blocks in the permute package, but some more flexible alternatives may be more appropriate.
...	Other arguments passed to vegdist.

Details

adonis2 is a function for the analysis and partitioning sums of squares using dissimilarities. The function is based on the principles of McArdle & Anderson (2001) and can perform sequential, marginal and overall tests. The function also allows using additive constants or squareroot of dissimilarities to avoid negative eigenvalues, but can also handle semimetric indices (such as Bray-Curtis) that produce negative eigenvalues. The adonis2 tests are identical to anova.cca of dbrda. With Euclidean distances, the tests are also identical to anova.cca of rda.

The function partitions sums of squares of a multivariate data set, and they are directly analogous to MANOVA (multivariate analysis of variance). McArdle and Anderson (2001) and Anderson (2001) refer to the method as “permutational MANOVA” (formerly “nonparametric MANOVA”). Further, as the inputs are linear predictors, and a response matrix of an arbitrary number of columns, they are a robust alternative to both parametric MANOVA and to ordination methods for describing how variation is attributed to different experimental treatments or uncontrolled covariates. The method is also analogous to distance-based redundancy analysis and algorithmically similar to dbrda (Legendre and Anderson 1999), and provides an alternative to AMOVA (nested analysis of molecular variance, Excoffier, Smouse, and Quattro, 1992; amova in the ade4 package) for both crossed and nested factors.

Value

The function returns an anova.cca result object with a new column for partial $R^2$: This is the proportion of sum of squares from the total, and in marginal models (by = "margin") the $R^2$ terms do not add up to 1.

Note

Anderson (2001, Fig. 4) warns that the method may confound location and dispersion effects: significant differences may be caused by different within-group variation (dispersion) instead of different mean values of the groups (see Warton et al. 2012 for a general analysis). However, it seems that adonis2 is less sensitive to dispersion effects than some of its alternatives (anosim, mrpp). Function betadisper is a sister function to adonis2 to study the differences in dispersion within the same geometric framework.

Author(s)

Martin Henry H. Stevens and Jari Oksanen.

References

Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46.

Excoffier, L., P.E. Smouse, and J.M. Quattro. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics, 131:479–491.

Legendre, P. and M.J. Anderson. 1999. Distance-based redundancy analysis: Testing multispecies responses in multifactorial ecological experiments. Ecological Monographs, 69:1–24.

McArdle, B.H. and M.J. Anderson. 2001. Fitting multivariate models to community data: A comment on distance-based redundancy analysis. Ecology, 82: 290–297.

Warton, D.I., Wright, T.W., Wang, Y. 2012. Distance-based multivariate analyses confound location and dispersion effects. Methods in Ecology and Evolution, 3, 89–101.

See Also

mrpp, anosim, mantel, varpart.

Examples

data(dune)
data(dune.env)
## default is overall (omnibus) test
adonis2(dune ~ Management*A1, data = dune.env)
## sequential tests
adonis2(dune ~ Management*A1, data = dune.env, by = "terms")

### Example of use with strata, for nested (e.g., block) designs.
dat <- expand.grid(rep=gl(2,1), NO3=factor(c(0,10)),field=gl(3,1) )
dat
Agropyron <- with(dat, as.numeric(field) + as.numeric(NO3)+2) +rnorm(12)/2
Schizachyrium <- with(dat, as.numeric(field) - as.numeric(NO3)+2) +rnorm(12)/2
total <- Agropyron + Schizachyrium
dotplot(total ~ NO3, dat, jitter.x=TRUE, groups=field,
        type=c('p','a'), xlab="NO3", auto.key=list(columns=3, lines=TRUE) )

Y <- data.frame(Agropyron, Schizachyrium)
mod <- metaMDS(Y, trace = FALSE)
plot(mod)
### Ellipsoid hulls show treatment
with(dat, ordiellipse(mod, NO3, kind = "ehull", label = TRUE))
### Spider shows fields
with(dat, ordispider(mod, field, lty=3, col="red", label = TRUE))

### Incorrect (no strata)
adonis2(Y ~ NO3, data = dat, permutations = 199)
## Correct with strata
with(dat, adonis2(Y ~ NO3, data = dat, permutations = 199, strata = field))

---

Analysis of Similarities

Description

Analysis of similarities (ANOSIM) provides a way to test statistically whether there is a significant difference between two or more groups of sampling units.

Usage

anosim(x, grouping, permutations = 999, distance = "bray", strata = NULL,
    parallel = getOption("mc.cores"))

Arguments

x	Data matrix or data frame in which rows are samples and columns are response variable(s), or a dissimilarity object or a symmetric square matrix of dissimilarities.
grouping	Factor for grouping observations.
permutations	a list of control values for the permutations as returned by the function how, or the number of permutations required, or a permutation matrix where each row gives the permuted indices.
distance	Choice of distance metric that measures the dissimilarity between two observations. See vegdist for options. This will be used if x was not a dissimilarity structure or a symmetric square matrix.
strata	An integer vector or factor specifying the strata for permutation. If supplied, observations are permuted only within the specified strata.
parallel	Number of parallel processes or a predefined socket cluster. With parallel = 1 uses ordinary, non-parallel processing. The parallel processing is done with parallel package.

Details

Analysis of similarities (ANOSIM) provides a way to test statistically whether there is a significant difference between two or more groups of sampling units. Function anosim operates directly on a dissimilarity matrix. A suitable dissimilarity matrix is produced by functions dist or vegdist. The method is philosophically allied with NMDS ordination (monoMDS), in that it uses only the rank order of dissimilarity values.

If two groups of sampling units are really different in their species composition, then compositional dissimilarities between the groups ought to be greater than those within the groups. The anosim statistic $R$ is based on the difference of mean ranks between groups ($r_B$) and within groups ($r_W$):

$R = (r_B - r_W)/(N (N-1) / 4)$

The divisor is chosen so that $R$ will be in the interval $-1 \dots +1$, value $0$ indicating completely random grouping.

The statistical significance of observed $R$ is assessed by permuting the grouping vector to obtain the empirical distribution of $R$ under null-model. See permutations for additional details on permutation tests in Vegan. The distribution of simulated values can be inspected with the permustats function.

The function has summary and plot methods. These both show valuable information to assess the validity of the method: The function assumes that all ranked dissimilarities within groups have about equal median and range. The plot method uses boxplot with options notch=TRUE and varwidth=TRUE.

Value

The function returns a list of class "anosim" with following items:

call	Function call.
statistic	The value of ANOSIM statistic $R$
signif	Significance from permutation.
perm	Permutation values of $R$. The distribution of permutation values can be inspected with function permustats.
class.vec	Factor with value Between for dissimilarities between classes and class name for corresponding dissimilarity within class.
dis.rank	Rank of dissimilarity entry.
dissimilarity	The name of the dissimilarity index: the "method" entry of the dist object.
control	A list of control values for the permutations as returned by the function how.

Note

The anosim function can confound the differences between groups and dispersion within groups and the results can be difficult to interpret (cf. Warton et al. 2012). The function returns a lot of information to ease studying its performance. Most anosim models could be analysed with adonis2 which seems to be a more robust alternative.

Author(s)

Jari Oksanen, with a help from Peter R. Minchin.

References

Clarke, K. R. (1993). Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18, 117–143.

Warton, D.I., Wright, T.W., Wang, Y. 2012. Distance-based multivariate analyses confound location and dispersion effects. Methods in Ecology and Evolution, 3, 89–101

See Also

mrpp for a similar function using original dissimilarities instead of their ranks. dist and vegdist for obtaining dissimilarities, and rank for ranking real values. For comparing dissimilarities against continuous variables, see mantel. Function adonis2 is a more robust alternative that should preferred.

Examples

data(dune)
data(dune.env)
dune.dist <- vegdist(dune)
dune.ano <- with(dune.env, anosim(dune.dist, Management))
summary(dune.ano)
plot(dune.ano)

---

Permutation Test for Constrained Correspondence Analysis, Redundancy Analysis and Constrained Analysis of Principal Coordinates

Description

The function performs an ANOVA like permutation test for Constrained Correspondence Analysis (cca), Redundancy Analysis (rda) or distance-based Redundancy Analysis (dbRDA, dbrda) to assess the significance of constraints.

Usage

## S3 method for class 'cca'
anova(object, ..., permutations = how(nperm=999),
     by = NULL, model = c("reduced", "direct", "full"),
     parallel = getOption("mc.cores"), strata = NULL,
     cutoff = 1, scope = NULL)
## S3 method for class 'cca'
permutest(x, permutations = how(nperm = 99),
     model = c("reduced", "direct", "full"), by = NULL, first = FALSE,
     strata = NULL, parallel = getOption("mc.cores"), ...)

Arguments

object	One or several result objects from cca, rda, dbrda or capscale. If there are several result objects, they are compared against each other in the order they were supplied. For a single object, a test specified in by or an overall test is given.
x	A single ordination result object.
permutations	a list of control values for the permutations as returned by the function how, or the number of permutations required, or a permutation matrix where each row gives the permuted indices.
by	Setting by = "axis" will assess significance for each constrained axis, and setting by = "terms" will assess significance for each term (sequentially from first to last), and setting by = "margin" will assess the marginal effects of the terms (each marginal term analysed in a model with all other variables), and by = "onedf" will assess sequentially one-degree-of-freedom contrasts of split factors.
model	Permutation model: model="direct" permutes community data, model="reduced" permutes residuals of the community data after Conditions (partial model), model = "full" permutes residuals after Conditions and Constraints.
parallel	Use parallel processing with the given number of cores.
strata	An integer vector or factor specifying the strata for permutation. If supplied, observations are permuted only within the specified strata. It is an error to use this when permutations is a matrix, or a how defines blocks. This is a legacy argument that will be deprecated in the future: use permutations = how(..., blocks) instead.
cutoff	Only effective with by="axis" where stops permutations after an axis equals or exceeds the cutoff $p$-value.
scope	Only effective with by="margin" where it can be used to select the marginal terms for testing. The default is to test all marginal terms in drop.scope.
first	Analyse only significance of the first axis.
...	Parameters passed to other functions. anova.cca passes all arguments to permutest.cca. In anova with by = "axis" you can use argument cutoff (defaults 1) which stops permutations after exceeding the given level.

Details

Functions anova.cca and permutest.cca implement ANOVA like permutation tests for the joint effect of constraints in cca, rda, dbrda or capscale. Function anova is intended as a more user-friendly alternative to permutest (that is the real workhorse).

Function anova can analyse a sequence of constrained ordination models. The analysis is based on the differences in residual deviance in permutations of nested models.

The default test is for the sum of all constrained eigenvalues. Setting first = TRUE will perform a test for the first constrained eigenvalue. Argument first can be set either in anova.cca or in permutest.cca. It is also possible to perform significance tests for each axis or for each term (constraining variable) using argument by in anova.cca. Setting by = "axis" will perform separate significance tests for each constrained axis. All previous constrained axes will be used as conditions (“partialled out”) and a test for the first constrained eigenvalues is performed (Legendre et al. 2011). You can stop permutation tests after exceeding a given significance level with argument cutoff to speed up calculations in large models. Setting by = "terms" will perform separate significance test for each term (constraining variable). The terms are assessed sequentially from first to last, and the order of the terms will influence their significances. Setting by = "onedf" will perform a similar sequential test for one-degree-of-freedom effects, where multi-level factors are split in their contrasts. Setting by = "margin" will perform separate significance test for each marginal term in a model with all other terms. The marginal test also accepts a scope argument for the drop.scope which can be a character vector of term labels that are analysed, or a fitted model of lower scope. The marginal effects are also known as “Type III” effects, but the current function only evaluates marginal terms. It will, for instance, ignore main effects that are included in interaction terms. In calculating pseudo-$F$, all terms are compared to the same residual of the full model.

Community data are permuted with choice model="direct", and residuals after partial CCA/ RDA/ dbRDA with choice model="reduced" (default). If there is no partial CCA/ RDA/ dbRDA stage, model="reduced" simply permutes the data and is equivalent to model="direct". The test statistic is “pseudo-$F$”, which is the ratio of constrained and unconstrained total Inertia (Chi-squares, variances or something similar), each divided by their respective degrees of freedom. If there are no conditions (“partial” terms), the sum of all eigenvalues remains constant, so that pseudo-$F$ and eigenvalues would give equal results. In partial CCA/ RDA/ dbRDA, the effect of conditioning variables (“covariables”) is removed before permutation, and the total Chi-square is not fixed, and test based on pseudo-$F$ would differ from the test based on plain eigenvalues.

Value

The function anova.cca calls permutest.cca and fills an anova table. Additional attributes are Random.seed (the random seeds used), control (the permutation design, see how) and F.perm (the permuted test statistics).

Author(s)

Jari Oksanen

References

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

Legendre, P., Oksanen, J. and ter Braak, C.J.F. (2011). Testing the significance of canonical axes in redundancy analysis. Methods in Ecology and Evolution 2, 269–277.

See Also

anova.cca, cca, rda, dbrda to get something to analyse. Function drop1.cca calls anova.cca with by = "margin", and add1.cca an analysis for single terms additions, which can be used in automatic or semiautomatic model building (see deviance.cca).

Examples

data(dune, dune.env)
mod <- cca(dune ~ Moisture + Management, dune.env)
## overall test
anova(mod)
## tests for individual terms
anova(mod, by="term")
anova(mod, by="margin")
## sequential test for contrasts
anova(mod, by = "onedf")
## test for adding all environmental variables
anova(mod, cca(dune ~ ., dune.env))

---

Averaged Subsampled Dissimilarity Matrices

Description

The function computes the dissimilarity matrix of a dataset multiple times using vegdist while randomly subsampling the dataset each time. All of the subsampled iterations are then averaged (mean) to provide a distance matrix that represents the average of multiple subsampling iterations. This emulates the behavior of the distance matrix calculator within the Mothur microbial ecology toolkit.

Usage

avgdist(x, sample, distfun = vegdist, meanfun = mean,
    transf = NULL, iterations = 100, dmethod = "bray",
    diag = TRUE, upper = TRUE, ...)

Arguments

x	Community data matrix.
sample	The subsampling depth to be used in each iteration. Samples that do not meet this threshold will be removed from the analysis, and their identity returned to the user in stdout.
distfun	The dissimilarity matrix function to be used. Default is the vegan vegdist
meanfun	The calculation to use for the average (mean or median).
transf	Option for transforming the count data before calculating the distance matrix. Any base transformation option can be used (e.g. sqrt)
iterations	The number of random iterations to perform before averaging. Default is 100 iterations.
dmethod	Dissimilarity index to be used with the specified dissimilarity matrix function. Default is Bray-Curtis
diag, upper	Return dissimilarities with diagonal and upper triangle. NB. the default differs from vegdist and returns symmetric "dist" structure instead of lower diagonal. However, the object cannot be accessed with matrix indices unless cast to matrix with as.matrix.
...	Any additional arguments to add to the distance function or mean/median function specified.

Note

The function builds on the function rrarefy and and additional distance matrix function (e.g. vegdist) to add more meaningful representations of distances among randomly subsampled datasets by presenting the average of multiple random iterations. This function runs using the vegdist. This functionality has been utilized in the Mothur standalone microbial ecology toolkit here.

Author(s)

Geoffrey Hannigan, with some minor tweaks by Gavin L. Simpson.

See Also

This function utilizes the vegdist and rrarefy functions.

Examples

# Import an example count dataset
data(BCI)
# Test the base functionality
mean.avg.dist <- avgdist(BCI, sample = 50, iterations = 10)
# Test the transformation function
mean.avg.dist.t <- avgdist(BCI, sample = 50, iterations = 10, transf = sqrt)
# Test the median functionality
median.avg.dist <- avgdist(BCI, sample = 50, iterations = 10, meanfun = median)
# Print the resulting tables
head(as.matrix(mean.avg.dist))
head(as.matrix(mean.avg.dist.t))
head(as.matrix(median.avg.dist))
# Run example to illustrate low variance of mean, median, and stdev results
# Mean and median std dev are around 0.05
sdd <- avgdist(BCI, sample = 50, iterations = 100, meanfun = sd)
summary(mean.avg.dist)
summary(median.avg.dist)
summary(sdd)
# Test for when subsampling depth excludes some samples
# Return samples that are removed for not meeting depth filter
depth.avg.dist <- avgdist(BCI, sample = 450, iterations = 10)
# Print the result
depth.avg.dist

---

Barro Colorado Island Tree Counts

Description

Tree counts in 1-hectare plots in the Barro Colorado Island and associated site information.

Usage

data(BCI)
data(BCI.env)

Format

A data frame with 50 plots (rows) of 1 hectare with counts of trees on each plot with total of 225 species (columns). Full Latin names are used for tree species. The names were updated against http://www.theplantlist.org and Kress et al. (2009) which allows matching 207 of species against doi:10.5061/dryad.63q27 (Zanne et al., 2014). The original species names are available as attribute original.names of BCI. See Examples for changed names.

For BCI.env, a data frame with 50 plots (rows) and nine site variables derived from Pyke et al. (2001) and Harms et al. (2001):

UTM.EW:

UTM coordinates (zone 17N) East-West.

UTM.NS:

UTM coordinates (zone 17N) North-South.

Precipitation:

Precipitation in mm per year.

Elevation:

Elevation in m above sea level.

Age.cat:

Forest age category.

Geology:

The Underlying geological formation.

Habitat:

Dominant habitat type based on the map of habitat types in 25 grid cells in each plot (Harms et al. 2001, excluding streamside habitat). The habitat types are Young forests (ca. 100 years), old forests on > 7 degree slopes (OldSlope), old forests under 152 m elevation (OldLow) and at higher elevation (OldHigh) and Swamp forests.

River:

"Yes" if there is streamside habitat in the plot.

EnvHet:

Environmental Heterogeneity assessed as the Simpson diversity of frequencies of Habitat types in 25 grid cells in the plot.

Details

Data give the numbers of trees at least 10 cm in diameter at breast height (DBH) in each one hectare quadrat in the 1982 BCI plot. Within each plot, all individuals were tallied and are recorded in this table. The full survey included smaller trees with DBH 1 cm or larger, but the BCI dataset is a subset of larger trees as compiled by Condit et al. (2002). The full data with thinner trees has densities above 4000 stems per hectare, or about ten times more stems than these data. The dataset BCI was provided (in 2003) to illustrate analysis methods in vegan. For scientific research on ecological issues we strongly recommend to access complete and more modern data (Condit et al. 2019) with updated taxonomy (Condit et al. 2020).

The data frame contains only the Barro Colorado Island subset of the full data table of Condit et al. (2002).

The quadrats are located in a regular grid. See BCI.env for the coordinates.

A full description of the site information in BCI.env is given in Pyke et al. (2001) and Harms et al. (2001). N.B. Pyke et al. (2001) and Harms et al. (2001) give conflicting information about forest age categories and elevation.

Source

https://www.science.org/doi/10.1126/science.1066854 for community data and References for environmental data. For updated complete data (incl. thinner trees down to 1 cm), see Condit et al. (2019).

References

Condit, R, Pitman, N, Leigh, E.G., Chave, J., Terborgh, J., Foster, R.B., Nuñez, P., Aguilar, S., Valencia, R., Villa, G., Muller-Landau, H.C., Losos, E. & Hubbell, S.P. (2002). Beta-diversity in tropical forest trees. Science 295, 666–669.

Condit R., Pérez, R., Aguilar, S., Lao, S., Foster, R. & Hubbell, S. (2019). Complete data from the Barro Colorado 50-ha plot: 423617 trees, 35 years [Dataset]. Dryad. doi:10.15146/5xcp-0d46

Condit, R., Aguilar, S., Lao, S., Foster, R., Hubbell, S. (2020). BCI 50-ha Plot Taxonomy [Dataset]. Dryad. doi:10.15146/R3FH61

Harms K.E., Condit R., Hubbell S.P. & Foster R.B. (2001) Habitat associations of trees and shrubs in a 50-ha neotropical forest plot. J. Ecol. 89, 947–959.

Kress W.J., Erickson D.L, Jones F.A., Swenson N.G, Perez R., Sanjur O. & Bermingham E. (2009) Plant DNA barcodes and a community phylogeny of a tropical forest dynamics plot in Panama. PNAS 106, 18621–18626.

Pyke, C. R., Condit, R., Aguilar, S., & Lao, S. (2001). Floristic composition across a climatic gradient in a neotropical lowland forest. Journal of Vegetation Science 12, 553–566. doi:10.2307/3237007

Zanne A.E., Tank D.C., Cornwell, W.K., Eastman J.M., Smith, S.A., FitzJohn, R.G., McGlinn, D.J., O’Meara, B.C., Moles, A.T., Reich, P.B., Royer, D.L., Soltis, D.E., Stevens, P.F., Westoby, M., Wright, I.J., Aarssen, L., Bertin, R.I., Calaminus, A., Govaerts, R., Hemmings, F., Leishman, M.R., Oleksyn, J., Soltis, P.S., Swenson, N.G., Warman, L. & Beaulieu, J.M. (2014) Three keys to the radiation of angiosperms into freezing environments. Nature 506, 89–92. doi:10.1038/nature12872 (published online Dec 22, 2013).

See Also

Extra-CRAN package natto (https://github.com/jarioksa/natto) has data set BCI.env2 with original grid data of Harms et al. (2001) habitat classification, and data set BCI.taxon of APG III classification of tree species.

Examples

data(BCI, BCI.env)
head(BCI.env)
## see changed species names
oldnames <- attr(BCI, "original.names")
taxa <- cbind("Old Names" = oldnames, "Current Names" = names(BCI))
noquote(taxa[taxa[,1] != taxa[,2], ])

---

Beals Smoothing and Degree of Absence

Description

Beals smoothing replaces each entry in the community data with a probability of a target species occurring in that particular site, based on the joint occurrences of the target species with the species that actually occur in the site. Swan's (1970) degree of absence applies Beals smoothing to zero items so long that all zeros are replaced with smoothed values.

Usage

beals(x, species = NA, reference = x, type = 0, include = TRUE)
swan(x, maxit = Inf, type = 0)

Arguments

x	Community data frame or matrix.
species	Column index used to compute Beals function for a single species. The default (NA) indicates that the function will be computed for all species.
reference	Community data frame or matrix to be used to compute joint occurrences. By default, x is used as reference to compute the joint occurrences.
type	Numeric. Specifies if and how abundance values have to be used in function beals. See details for more explanation.
include	This logical flag indicates whether the target species has to be included when computing the mean of the conditioned probabilities. The original Beals (1984) definition is equivalent to include=TRUE, while the formulation of Münzbergová and Herben is equal to include=FALSE.
maxit	Maximum number of iterations. The default Inf means that iterations are continued until there are no zeros or the number of zeros does not change. Probably only maxit = 1 makes sense in addition to the default.

Details

Beals smoothing is the estimated probability $p_{ij}$ that species $j$ occurs at site $i$. It is defined as $p_{ij} = \frac{1}{S_i} \sum_k \frac{N_{jk} I_{ik}}{N_k}$, where $S_i$ is the number of species at site $i$, $N_{jk}$ is the number of joint occurrences of species $j$ and $k$, $N_k$ is the number of occurrences of species $k$, and $I$ is the incidence (0 or 1) of species (this last term is usually omitted from the equation, but it is necessary). As $N_{jk}$ can be interpreted as a mean of conditional probability, the beals function can be interpreted as a mean of conditioned probabilities (De Cáceres & Legendre 2008). The present function is generalized to abundance values (De Cáceres & Legendre 2008).

The type argument specifies if and how abundance values have to be used. type = 0 presence/absence mode. type = 1 abundances in reference (or x) are used to compute conditioned probabilities. type = 2 abundances in x are used to compute weighted averages of conditioned probabilities. type = 3 abundances are used to compute both conditioned probabilities and weighted averages.

Beals smoothing was originally suggested as a method of data transformation to remove excessive zeros (Beals 1984, McCune 1994). However, it is not a suitable method for this purpose since it does not maintain the information on species presences: a species may have a higher probability of occurrence at a site where it does not occur than at sites where it occurs. Moreover, it regularizes data too strongly. The method may be useful in identifying species that belong to the species pool (Ewald 2002) or to identify suitable unoccupied patches in metapopulation analysis (Münzbergová & Herben 2004). In this case, the function should be called with include=FALSE for cross-validation smoothing for species; argument species can be used if only one species is studied.

Swan (1970) suggested replacing zero values with degrees of absence of a species in a community data matrix. Swan expressed the method in terms of a similarity matrix, but it is equivalent to applying Beals smoothing to zero values, at each step shifting the smallest initially non-zero item to value one, and repeating this so many times that there are no zeros left in the data. This is actually very similar to extended dissimilarities (implemented in function stepacross), but very rarely used.

Value

The function returns a transformed data matrix or a vector if Beals smoothing is requested for a single species.

Author(s)

Miquel De Cáceres and Jari Oksanen

References

Beals, E.W. 1984. Bray-Curtis ordination: an effective strategy for analysis of multivariate ecological data. Pp. 1–55 in: MacFadyen, A. & E.D. Ford [eds.] Advances in Ecological Research, 14. Academic Press, London.

De Cáceres, M. & Legendre, P. 2008. Beals smoothing revisited. Oecologia 156: 657–669.

Ewald, J. 2002. A probabilistic approach to estimating species pools from large compositional matrices. J. Veg. Sci. 13: 191–198.

McCune, B. 1994. Improving community ordination with the Beals smoothing function. Ecoscience 1: 82–86.

Münzbergová, Z. & Herben, T. 2004. Identification of suitable unoccupied habitats in metapopulation studies using co-occurrence of species. Oikos 105: 408–414.

Swan, J.M.A. 1970. An examination of some ordination problems by use of simulated vegetational data. Ecology 51: 89–102.

See Also

decostand for proper standardization methods, specpool for an attempt to assess the size of species pool. Function indpower assesses the power of each species to estimate the probabilities predicted by beals.

Examples

data(dune)
## Default
x <- beals(dune)
## Remove target species
x <- beals(dune, include = FALSE)
## Smoothed values against presence or absence of species
pa <- decostand(dune, "pa")
boxplot(as.vector(x) ~ unlist(pa), xlab="Presence", ylab="Beals")
## Remove the bias of tarbet species: Yields lower values.
beals(dune, type =3, include = FALSE)
## Uses abundance information.
## Vector with beals smoothing values corresponding to the first species
## in dune.
beals(dune, species=1, include=TRUE)

---

Multivariate homogeneity of groups dispersions (variances)

Description

Implements Marti Anderson's PERMDISP2 procedure for the analysis of multivariate homogeneity of group dispersions (variances). betadisper is a multivariate analogue of Levene's test for homogeneity of variances. Non-euclidean distances between objects and group centres (centroids or medians) are handled by reducing the original distances to principal coordinates. This procedure has latterly been used as a means of assessing beta diversity. There are anova, scores, plot and boxplot methods.

TukeyHSD.betadisper creates a set of confidence intervals on the differences between the mean distance-to-centroid of the levels of the grouping factor with the specified family-wise probability of coverage. The intervals are based on the Studentized range statistic, Tukey's 'Honest Significant Difference' method.

Usage

betadisper(d, group, type = c("median","centroid"), bias.adjust = FALSE,
       sqrt.dist = FALSE, add = FALSE)

## S3 method for class 'betadisper'
anova(object, ...)

## S3 method for class 'betadisper'
scores(x, display = c("sites", "centroids"),
       choices = c(1,2), ...)

## S3 method for class 'betadisper'
eigenvals(x, ...)

## S3 method for class 'betadisper'
plot(x, axes = c(1,2), cex = 0.7,
     pch = seq_len(ng), col = NULL, lty = "solid", lwd = 1, hull = TRUE,
     ellipse = FALSE, conf,
     segments = TRUE, seg.col = "grey", seg.lty = lty, seg.lwd = lwd,
     label = TRUE, label.cex = 1,
     ylab, xlab, main, sub, ...)

## S3 method for class 'betadisper'
boxplot(x, ylab = "Distance to centroid", ...)

## S3 method for class 'betadisper'
TukeyHSD(x, which = "group", ordered = FALSE,
         conf.level = 0.95, ...)

## S3 method for class 'betadisper'
print(x, digits = max(3, getOption("digits") - 3),
                           neigen = 8, ...)

Arguments

d	a distance structure such as that returned by dist, betadiver or vegdist.
group	vector describing the group structure, usually a factor or an object that can be coerced to a factor using as.factor. Can consist of a factor with a single level (i.e., one group).
type	the type of analysis to perform. Use the spatial median or the group centroid? The spatial median is now the default.
bias.adjust	logical: adjust for small sample bias in beta diversity estimates?
sqrt.dist	Take square root of dissimilarities. This often euclidifies dissimilarities.
add	Add a constant to the non-diagonal dissimilarities such that all eigenvalues are non-negative in the underlying Principal Co-ordinates Analysis (see wcmdscale for details). Choice "lingoes" (or TRUE) use the recommended method of Legendre & Anderson (1999: “method 1”) and "cailliez" uses their “method 2”.
display	character; partial match to access scores for "sites" or "species".
object, x	an object of class "betadisper", the result of a call to betadisper.
choices, axes	the principal coordinate axes wanted.
hull	logical; should the convex hull for each group be plotted?
ellipse	logical; should the standard deviation data ellipse for each group be plotted?
conf	Expected fractions of data coverage for data ellipses, e.g. 0.95. The default is to draw a 1 standard deviation data ellipse, but if supplied, conf is multiplied with the corresponding value found from the Chi-squared distribution with 2df to provide the requested coverage (probability contour).
pch	plot symbols for the groups, a vector of length equal to the number of groups.
col	colors for the plot symbols and centroid labels for the groups, a vector of length equal to the number of groups.
lty, lwd	linetype, linewidth for convex hulls and confidence ellipses.
segments	logical; should segments joining points to their centroid be drawn?
seg.col	colour to draw segments between points and their centroid. Can be a vector, in which case one colour per group.
seg.lty, seg.lwd	linetype and line width for segments.
label	logical; should the centroids by labelled with their respective factor label?
label.cex	numeric; character expansion for centroid labels.
cex, ylab, xlab, main, sub	graphical parameters. For details, see plot.default.
which	A character vector listing terms in the fitted model for which the intervals should be calculated. Defaults to the grouping factor.
ordered	logical; see TukeyHSD.
conf.level	A numeric value between zero and one giving the family-wise confidence level to use.
digits, neigen	numeric; for the print method, sets the number of digits to use (as per print.default) and the maximum number of axes to display eigenvalues for, repsectively.
...	arguments, including graphical parameters (for plot.betadisper and boxplot.betadisper), passed to other methods.

Details

One measure of multivariate dispersion (variance) for a group of samples is to calculate the average distance of group members to the group centroid or spatial median (both referred to as 'centroid' from now on unless stated otherwise) in multivariate space. To test if the dispersions (variances) of one or more groups are different, the distances of group members to the group centroid are subject to ANOVA. This is a multivariate analogue of Levene's test for homogeneity of variances if the distances between group members and group centroids is the Euclidean distance.

However, better measures of distance than the Euclidean distance are available for ecological data. These can be accommodated by reducing the distances produced using any dissimilarity coefficient to principal coordinates, which embeds them within a Euclidean space. The analysis then proceeds by calculating the Euclidean distances between group members and the group centroid on the basis of the principal coordinate axes rather than the original distances.

Non-metric dissimilarity coefficients can produce principal coordinate axes that have negative Eigenvalues. These correspond to the imaginary, non-metric part of the distance between objects. If negative Eigenvalues are produced, we must correct for these imaginary distances.

The distance to its centroid of a point is

$z_{ij}^c = \sqrt{\Delta^2(u_{ij}^+, c_i^+) - \Delta^2(u_{ij}^-, c_i^-)},$

where $\Delta^2$ is the squared Euclidean distance between $u_{ij}$, the principal coordinate for the $j$th point in the $i$th group, and $c_i$, the coordinate of the centroid for the $i$th group. The super-scripted ‘$+$’ and ‘$-$’ indicate the real and imaginary parts respectively. This is equation (3) in Anderson (2006). If the imaginary part is greater in magnitude than the real part, then we would be taking the square root of a negative value, resulting in NaN, and these cases are changed to zero distances (with a warning). This is in line with the behaviour of Marti Anderson's PERMDISP2 programme.

To test if one or more groups is more variable than the others, ANOVA of the distances to group centroids can be performed and parametric theory used to interpret the significance of $F$. An alternative is to use a permutation test. permutest.betadisper permutes model residuals to generate a permutation distribution of $F$ under the Null hypothesis of no difference in dispersion between groups.

Pairwise comparisons of group mean dispersions can also be performed using permutest.betadisper. An alternative to the classical comparison of group dispersions, is to calculate Tukey's Honest Significant Differences between groups, via TukeyHSD.betadisper. This is a simple wrapper to TukeyHSD. The user is directed to read the help file for TukeyHSD before using this function. In particular, note the statement about using the function with unbalanced designs.

The results of the analysis can be visualised using the plot and boxplot methods. The distances of points to all centroids (group) can be found with function centredist.

One additional use of these functions is in assessing beta diversity (Anderson et al 2006). Function betadiver provides some popular dissimilarity measures for this purpose.

As noted in passing by Anderson (2006) and in a related context by O'Neill (2000), estimates of dispersion around a central location (median or centroid) that is calculated from the same data will be biased downward. This bias matters most when comparing diversity among treatments with small, unequal numbers of samples. Setting bias.adjust=TRUE when using betadisper imposes a $\sqrt{n/(n-1)}$ correction (Stier et al. 2013).

Value

The anova method returns an object of class "anova" inheriting from class "data.frame".

The scores method returns a list with one or both of the components "sites" and "centroids".

The plot function invisibly returns an object of class "ordiplot", a plotting structure which can be used by identify.ordiplot (to identify the points) or other functions in the ordiplot family.

The boxplot function invisibly returns a list whose components are documented in boxplot.

eigenvals.betadisper returns a named vector of eigenvalues.

TukeyHSD.betadisper returns a list. See TukeyHSD for further details.

betadisper returns a list of class "betadisper" with the following components:

eig	numeric; the eigenvalues of the principal coordinates analysis.
vectors	matrix; the eigenvectors of the principal coordinates analysis.
distances	numeric; the Euclidean distances in principal coordinate space between the samples and their respective group centroid or median.
group	factor; vector describing the group structure
centroids	matrix; the locations of the group centroids or medians on the principal coordinates.
group.distances	numeric; the mean distance to each group centroid or median.
call	the matched function call.

Warning

Stewart Schultz noticed that the permutation test for type="centroid" had the wrong type I error and was anti-conservative. As such, the default for type has been changed to "median", which uses the spatial median as the group centroid. Tests suggests that the permutation test for this type of analysis gives the correct error rates.

Note

If group consists of a single level or group, then the anova and permutest methods are not appropriate and if used on such data will stop with an error.

Missing values in either d or group will be removed prior to performing the analysis.

Author(s)

Gavin L. Simpson; bias correction by Adrian Stier and Ben Bolker.

References

Anderson, M.J. (2006) Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62, 245–253.

Anderson, M.J., Ellingsen, K.E. & McArdle, B.H. (2006) Multivariate dispersion as a measure of beta diversity. Ecology Letters 9, 683–693.

O'Neill, M.E. (2000) A Weighted Least Squares Approach to Levene's Test of Homogeneity of Variance. Australian & New Zealand Journal of Statistics 42, 81-–100.

Stier, A.C., Geange, S.W., Hanson, K.M., & Bolker, B.M. (2013) Predator density and timing of arrival affect reef fish community assembly. Ecology 94, 1057–1068.

See Also

permutest.betadisper, centredist, anova.lm, scores, boxplot, TukeyHSD. Further measure of beta diversity can be found in betadiver.

Examples

data(varespec)

## Bray-Curtis distances between samples
dis <- vegdist(varespec)

## First 16 sites grazed, remaining 8 sites ungrazed
groups <- factor(c(rep(1,16), rep(2,8)), labels = c("grazed","ungrazed"))

## Calculate multivariate dispersions
mod <- betadisper(dis, groups)
mod

## Perform test
anova(mod)

## Permutation test for F
permutest(mod, pairwise = TRUE, permutations = 99)

## Tukey's Honest Significant Differences
(mod.HSD <- TukeyHSD(mod))
plot(mod.HSD)

## Plot the groups and distances to centroids on the
## first two PCoA axes
plot(mod)

## with data ellipses instead of hulls
plot(mod, ellipse = TRUE, hull = FALSE) # 1 sd data ellipse
plot(mod, ellipse = TRUE, hull = FALSE, conf = 0.90) # 90% data ellipse

# plot with manual colour specification
my_cols <- c("#1b9e77", "#7570b3")
plot(mod, col = my_cols, pch = c(16,17), cex = 1.1)

## can also specify which axes to plot, ordering respected
plot(mod, axes = c(3,1), seg.col = "forestgreen", seg.lty = "dashed")

## Draw a boxplot of the distances to centroid for each group
boxplot(mod)

## `scores` and `eigenvals` also work
scrs <- scores(mod)
str(scrs)
head(scores(mod, 1:4, display = "sites"))
# group centroids/medians 
scores(mod, 1:4, display = "centroids")
# eigenvalues from the underlying principal coordinates analysis
eigenvals(mod) 

## try out bias correction; compare with mod3
(mod3B <- betadisper(dis, groups, type = "median", bias.adjust=TRUE))
anova(mod3B)
permutest(mod3B, permutations = 99)

## should always work for a single group
group <- factor(rep("grazed", NROW(varespec)))
(tmp <- betadisper(dis, group, type = "median"))
(tmp <- betadisper(dis, group, type = "centroid"))

## simulate missing values in 'd' and 'group'
## using spatial medians
groups[c(2,20)] <- NA
dis[c(2, 20)] <- NA
mod2 <- betadisper(dis, groups) ## messages
mod2
permutest(mod2, permutations = 99)
anova(mod2)
plot(mod2)
boxplot(mod2)
plot(TukeyHSD(mod2))

## Using group centroids
mod3 <- betadisper(dis, groups, type = "centroid")
mod3
permutest(mod3, permutations = 99)
anova(mod3)
plot(mod3)
boxplot(mod3)
plot(TukeyHSD(mod3))

---

Indices of beta Diversity

Description

The function estimates any of the 24 indices of beta diversity reviewed by Koleff et al. (2003). Alternatively, it finds the co-occurrence frequencies for triangular plots (Koleff et al. 2003).

Usage

betadiver(x, method = NA, order = FALSE, help = FALSE, ...)
## S3 method for class 'betadiver'
plot(x, ...)
## S3 method for class 'betadiver'
scores(x, triangular = TRUE, ...)

Arguments

x	Community data matrix, or the betadiver result for plot and scores functions.
method	The index of beta diversity as defined in Koleff et al. (2003), Table 1. You can use either the subscript of $\beta$ or the number of the index. See argument help below.
order	Order sites by increasing number of species. This will influence the configuration in the triangular plot and non-symmetric indices.
help	Show the numbers, subscript names and the defining equations of the indices and exit.
triangular	Return scores suitable for triangular plotting of proportions. If FALSE, returns a 3-column matrix of raw counts.
...	Other arguments to functions.

Details

The most commonly used index of beta diversity is $\beta_w = S/\alpha - 1$, where $S$ is the total number of species, and $\alpha$ is the average number of species per site (Whittaker 1960). A drawback of this model is that $S$ increases with sample size, but the expectation of $\alpha$ remains constant, and so the beta diversity increases with sample size. A solution to this problem is to study the beta diversity of pairs of sites (Marion et al. 2017). If we denote the number of species shared between two sites as $a$ and the numbers of unique species (not shared) as $b$ and $c$, then $S = a + b + c$ and $\alpha = (2 a + b + c)/2$ so that $\beta_w = (b+c)/(2 a + b + c)$. This is the Sørensen dissimilarity as defined in vegan function vegdist with argument binary = TRUE. Many other indices are dissimilarity indices as well.

Function betadiver finds all indices reviewed by Koleff et al. (2003). All these indices could be found with function designdist, but the current function provides a conventional shortcut. The function only finds the indices. The proper analysis must be done with functions such as betadisper, adonis2 or mantel.

The indices are directly taken from Table 1 of Koleff et al. (2003), and they can be selected either by the index number or the subscript name used by Koleff et al. The numbers, names and defining equations can be seen using betadiver(help = TRUE). In all cases where there are two alternative forms, the one with the term $-1$ is used. There are several duplicate indices, and the number of distinct alternatives is much lower than 24 formally provided. The formulations used in functions differ occasionally from those in Koleff et al. (2003), but they are still mathematically equivalent. With method = NA, no index is calculated, but instead an object of class betadiver is returned. This is a list of elements a, b and c. Function plot can be used to display the proportions of these elements in triangular plot as suggested by Koleff et al. (2003), and scores extracts the triangular coordinates or the raw scores. Function plot returns invisibly the triangular coordinates as an "ordiplot" object.

Value

With method = NA, the function returns an object of class "betadisper" with elements a, b, and c. If method is specified, the function returns a "dist" object which can be used in any function analysing dissimilarities. For beta diversity, particularly useful functions are betadisper to study the betadiversity in groups, adonis2 for any model, and mantel to compare beta diversities to other dissimilarities or distances (including geographical distances). Although betadiver returns a "dist" object, some indices are similarities and cannot be used as such in place of dissimilarities, but that is a user error. Functions 10 ("j"), 11 ("sor") and 21 ("rlb") are similarity indices. Function sets argument "maxdist" similarly as vegdist, using NA when there is no fixed upper limit, and 0 for similarities.

Warning

Some indices return similarities instead of dissimilarities.

Author(s)

Jari Oksanen

References

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

Koleff, P., Gaston, K.J. and Lennon, J.J. (2003) Measuring beta diversity for presence-absence data. Journal of Animal Ecology 72, 367–382.

Marion, Z.H., Fordyce, J.A. and Fitzpatrick, B.M. (2017) Pairwise beta diversity resolves an underappreciated source of confusion in calculating species turnover. Ecology 98, 933–939.

Whittaker, R.H. (1960) Vegetation of Siskiyou mountains, Oregon and California. Ecological Monographs 30, 279–338.

See Also

designdist can be used to implement all these functions, and also allows using notation with alpha and gamma diversities. vegdist has some canned alternatives. Functions betadisper, adonis2 and mantel can be used for analysing beta diversity objects. The returned dissimilarities can be used in any distance-based methods, such as metaMDS, capscale and dbrda. Functions nestedbetasor and nestedbetajac implement decomposition beta diversity measures (Sørensen and Jaccard) into turnover and nestedness components following Baselga (2010).

Examples

## Raw data and plotting
data(sipoo)
m <- betadiver(sipoo)
plot(m)
## The indices
betadiver(help=TRUE)
## The basic Whittaker index
d <- betadiver(sipoo, "w")
## This should be equal to Sorensen index (binary Bray-Curtis in
## vegan)
range(d - vegdist(sipoo, binary=TRUE))

---

Coefficients of Biogeographical Dispersal Direction

Description

This function computes coefficients of dispersal direction between geographically connected areas, as defined by Legendre and Legendre (1984), and also described in Legendre and Legendre (2012, section 13.3.4).

Usage

bgdispersal(mat, PAonly = FALSE, abc = FALSE)

Arguments

mat	Data frame or matrix containing a community composition data table (species presence-absence or abundance data).
PAonly	FALSE if the four types of coefficients, DD1 to DD4, are requested; TRUE if DD1 and DD2 only are sought (see Details).
abc	If TRUE, return tables a, b and c used in DD1 and DD2.

Details

The signs of the DD coefficients indicate the direction of dispersal, provided that the asymmetry is significant. A positive sign indicates dispersal from the first (row in DD tables) to the second region (column); a negative sign indicates the opposite. A McNemar test of asymmetry is computed from the presence-absence data to test the hypothesis of a significant asymmetry between the two areas under comparison.

In the input data table, the rows are sites or areas, the columns are taxa. Most often, the taxa are species, but the coefficients can be computed from genera or families as well. DD1 and DD2 only are computed for presence-absence data. The four types of coefficients are computed for quantitative data, which are converted to presence-absence for the computation of DD1 and DD2. PAonly = FALSE indicates that the four types of coefficients are requested. PAonly = TRUE if DD1 and DD2 only are sought.

Value

Function bgdispersal returns a list containing the following matrices:

DD1	$DD1_{j,k} = (a(b - c))/((a + b + c)^2)$
DD2	$DD2_{j,k} = (2 a (b - c))/((2a + b + c) (a + b + c))$ where $a$, $b$, and $c$ have the same meaning as in the computation of binary similarity coefficients.
DD3	$DD3_{j,k} = {W(A-B) / (A+B-W)^2}$
DD4	$DD4_{j,k} = 2W(A-B) / ((A+B)(A+B-W))$ where W = sum(pmin(vector1, vector2)), A = sum(vector1), B = sum(vector2)
McNemar	McNemar chi-square statistic of asymmetry (Sokal and Rohlf 1995): $2(b \log(b) + c \log(c) - (b+c) \log((b+c)/2)) / q$, where $q = 1 + 1/(2(b+c))$ (Williams correction for continuity)
prob.McNemar	probabilities associated with McNemar statistics, chi-square test. H0: no asymmetry in $(b-c)$.

Note

The function uses a more powerful alternative for the McNemar test than the classical formula. The classical formula was constructed in the spirit of Pearson's Chi-square, but the formula in this function was constructed in the spirit of Wilks Chi-square or the $G$ statistic. Function mcnemar.test uses the classical formula. The new formula was introduced in vegan version 1.10-11, and the older implementations of bgdispersal used the classical formula.

Author(s)

Pierre Legendre, Departement de Sciences Biologiques, Universite de Montreal

References

Legendre, P. and V. Legendre. 1984. Postglacial dispersal of freshwater fishes in the Québec peninsula. Can. J. Fish. Aquat. Sci. 41: 1781-1802.

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

Sokal, R. R. and F. J. Rohlf. 1995. Biometry. The principles and practice of statistics in biological research. 3rd edn. W. H. Freeman, New York.

Examples

mat <- matrix(c(32,15,14,10,70,30,100,4,10,30,25,0,18,0,40,
  0,0,20,0,0,0,0,4,0,30,20,0,0,0,0,25,74,42,1,45,89,5,16,16,20),
  4, 10, byrow=TRUE)
bgdispersal(mat)

---

Best Subset of Environmental Variables with Maximum (Rank) Correlation with Community Dissimilarities

Description

Function finds the best subset of environmental variables, so that the Euclidean distances of scaled environmental variables have the maximum (rank) correlation with community dissimilarities.

Usage

## Default S3 method:
bioenv(comm, env, method = "spearman", index = "bray",
       upto = ncol(env), trace = FALSE, partial = NULL, 
       metric = c("euclidean", "mahalanobis", "manhattan", "gower"),
       parallel = getOption("mc.cores"), ...)
## S3 method for class 'formula'
bioenv(formula, data, ...)
bioenvdist(x, which = "best")

Arguments

comm	Community data frame or a dissimilarity object or a square matrix that can be interpreted as dissimilarities.
env	Data frame of continuous environmental variables.
method	The correlation method used in cor.
index	The dissimilarity index used for community data (comm) in vegdist. This is ignored if comm are dissimilarities.
upto	Maximum number of parameters in studied subsets.
formula, data	Model formula and data.
trace	Trace the calculations
partial	Dissimilarities partialled out when inspecting variables in env.
metric	Metric used for distances of environmental distances. See Details.
parallel	Number of parallel processes or a predefined socket cluster. With parallel = 1 uses ordinary, non-parallel processing. The parallel processing is done with parallel package.
x	bioenv result object.
which	The number of the model for which the environmental distances are evaluated, or the "best" model.
...	Other arguments passed to vegdist.

Details

The function calculates a community dissimilarity matrix using vegdist. Then it selects all possible subsets of environmental variables, scales the variables, and calculates Euclidean distances for this subset using dist. The function finds the correlation between community dissimilarities and environmental distances, and for each size of subsets, saves the best result. There are $2^p-1$ subsets of $p$ variables, and an exhaustive search may take a very, very, very long time (parameter upto offers a partial relief).

The argument metric defines distances in the given set of environmental variables. With metric = "euclidean", the variables are scaled to unit variance and Euclidean distances are calculated. With metric = "mahalanobis", the Mahalanobis distances are calculated: in addition to scaling to unit variance, the matrix of the current set of environmental variables is also made orthogonal (uncorrelated). With metric = "manhattan", the variables are scaled to unit range and Manhattan distances are calculated, so that the distances are sums of differences of environmental variables. With metric = "gower", the Gower distances are calculated using function daisy. This allows also using factor variables, but with continuous variables the results are equal to metric = "manhattan".

The function can be called with a model formula where the LHS is the data matrix and RHS lists the environmental variables. The formula interface is practical in selecting or transforming environmental variables.

With argument partial you can perform “partial” analysis. The partializing item must be a dissimilarity object of class dist. The partial item can be used with any correlation method, but it is strictly correct only for Pearson.

Function bioenvdist recalculates the environmental distances used within the function. The default is to calculate distances for the best model, but the number of any model can be given.

Clarke & Ainsworth (1993) suggested this method to be used for selecting the best subset of environmental variables in interpreting results of nonmetric multidimensional scaling (NMDS). They recommended a parallel display of NMDS of community dissimilarities and NMDS of Euclidean distances from the best subset of scaled environmental variables. They warned against the use of Procrustes analysis, but to me this looks like a good way of comparing these two ordinations.

Clarke & Ainsworth wrote a computer program BIO-ENV giving the name to the current function. Presumably BIO-ENV was later incorporated in Clarke's PRIMER software (available for Windows). In addition, Clarke & Ainsworth suggested a novel method of rank correlation which is not available in the current function.

Value

The function returns an object of class bioenv with a summary method.

Note

If you want to study the ‘significance’ of bioenv results, you can use function mantel or mantel.partial which use the same definition of correlation. However, bioenv standardizes environmental variables depending on the used metric, and you must do the same in mantel for comparable results (the standardized data are returned as item x in the result object). It is safest to use bioenvdist to extract the environmental distances that really were used within bioenv. NB., bioenv selects variables to maximize the Mantel correlation, and significance tests based on a priori selection of variables are biased.

Author(s)

Jari Oksanen

References

Clarke, K. R & Ainsworth, M. 1993. A method of linking multivariate community structure to environmental variables. Marine Ecology Progress Series, 92, 205–219.

See Also

vegdist, dist, cor for underlying routines, monoMDS and metaMDS for ordination, procrustes for Procrustes analysis, protest for an alternative, and rankindex for studying alternatives to the default Bray-Curtis index.

Examples

# The method is very slow for large number of possible subsets.
# Therefore only 6 variables in this example.
data(varespec)
data(varechem)
sol <- bioenv(wisconsin(varespec) ~ log(N) + P + K + Ca + pH + Al, varechem)
sol
## IGNORE_RDIFF_BEGIN
summary(sol)
## IGNORE_RDIFF_END

---

PCA biplot

Description

Draws a PCA biplot with species scores indicated by biplot arrows

Usage

## S3 method for class 'rda'
biplot(x, choices = c(1, 2), scaling = "species",
       display = c("sites", "species"), type, xlim, ylim, col = c(1,2), 
       const, correlation = FALSE, ...)

Arguments

x	A rda result object.
choices	Axes to show.
scaling	Scaling for species and site scores. Either species (2) or site (1) scores are scaled by eigenvalues, and the other set of scores is left unscaled, or with 3 both are scaled symmetrically by square root of eigenvalues. With negative scaling values in rda, species scores are divided by standard deviation of each species and multiplied with an equalizing constant. Unscaled raw scores stored in the result can be accessed with scaling = 0. The type of scores can also be specified as one of "none", "sites", "species", or "symmetric", which correspond to the values 0, 1, 2, and 3 respectively. Argument correlation can be used in combination with these character descriptions to get the corresponding negative value.
correlation	logical; if scaling is a character description of the scaling type, correlation can be used to select correlation-like scores for PCA. See argument scaling for details.
display	Scores shown. These must some of the alternatives "species" for species scores, and/or "sites" for site scores.
type	Type of plot: partial match to text for text labels, points for points, and none for setting frames only. If omitted, text is selected for smaller data sets, and points for larger. Can be of length 2 (e.g. type = c("text", "points")), in which case the first element describes how species scores are handled, and the second how site scores are drawn.
xlim, ylim	the x and y limits (min, max) of the plot.
col	Colours used for sites and species (in this order). If only one colour is given, it is used for both.
const	General scaling constant for scores.rda.
...	Other parameters for plotting functions.

Details

Produces a plot or biplot of the results of a call to rda. It is common for the "species" scores in a PCA to be drawn as biplot arrows that point in the direction of increasing values for that variable. The biplot.rda function provides a wrapper to plot.cca to allow the easy production of such a plot.

biplot.rda is only suitable for unconstrained models. If used on an ordination object with constraints, an error is issued.

If species scores are drawn using "text", the arrows are drawn from the origin to 0.85 * species score, whilst the labels are drawn at the species score. If the type used is "points", then no labels are drawn and therefore the arrows are drawn from the origin to the actual species score.

Value

The plot function returns invisibly a plotting structure which can be used by identify.ordiplot to identify the points or other functions in the ordiplot family.

Author(s)

Gavin Simpson, based on plot.cca by Jari Oksanen.

See Also

plot.cca, rda for something to plot, ordiplot for an alternative plotting routine and more support functions, and text, points and arrows for the basic routines.

Examples

data(dune)
mod <- rda(dune, scale = TRUE)
biplot(mod, scaling = "symmetric")

## plot.cca can do the same
plot(mod, scaling = "symmetric", spe.par = list(arrows=TRUE))

## different type for species and site scores
biplot(mod, scaling = "symmetric", type = c("text", "points"))

## We can use ordiplot pipes to build similar plots with flexible
## control
plot(mod, scaling = "symmetric", type="n") |>
   points("sites", cex=0.7) |>
   text("species", arrows=TRUE, length=0.05, col=2, cex=0.7, font=3)

---

K-means partitioning using a range of values of K

Description

This function is a wrapper for the kmeans function. It creates several partitions forming a cascade from a small to a large number of groups.

Usage

cascadeKM(data, inf.gr, sup.gr, iter = 100, criterion = "calinski",
  parallel = getOption("mc.cores"))

cIndexKM(y, x, index = "all")

## S3 method for class 'cascadeKM'
plot(x, min.g, max.g, grpmts.plot = TRUE, 
     sortg = FALSE, gridcol = NA, ...)

Arguments

data	The data matrix. The objects (samples) are the rows.
inf.gr	The number of groups for the partition with the smallest number of groups of the cascade (min).
sup.gr	The number of groups for the partition with the largest number of groups of the cascade (max).
iter	The number of random starting configurations for each value of $K$.
criterion	The criterion that will be used to select the best partition. The default value is "calinski", which refers to the Calinski-Harabasz (1974) criterion. The simple structure index ("ssi") is also available. Other indices are available in package cclust. In our experience, the two indices that work best and are most likely to return their maximum value at or near the optimal number of clusters are "calinski" and "ssi".
y	Object of class "kmeans" returned by a clustering algorithm such as kmeans
x	Data matrix where columns correspond to variables and rows to observations, or the plotting object in plot
index	The available indices are: "calinski" and "ssi". Type "all" to obtain both indices. Abbreviations of these names are also accepted.
min.g, max.g	The minimum and maximum numbers of groups to be displayed.
grpmts.plot	Show the plot (TRUE or FALSE).
sortg	Sort the objects as a function of their group membership to produce a more easily interpretable graph. See Details. The original object names are kept; they are used as labels in the output table x, although not in the graph. If there were no row names, sequential row numbers are used to keep track of the original order of the objects.
gridcol	The colour of the grid lines in the plots. NA, which is the default value, removes the grid lines.
...	Other parameters to the functions (ignored).
parallel	Number of parallel processes or a predefined socket cluster. With parallel = 1 uses ordinary, non-parallel processing. The parallel processing is done with parallel package.

Details

The function creates several partitions forming a cascade from a small to a large number of groups formed by kmeans. Most of the work is performed by function cIndex which is based on the clustIndex in package cclust). Some of the criteria were removed from this version because computation errors were generated when only one object was found in a group.

The default value is "calinski", which refers to the well-known Calinski-Harabasz (1974) criterion. The other available index is the simple structure index "ssi" (Dolnicar et al. 1999). In the case of groups of equal sizes, "calinski" is generally a good criterion to indicate the correct number of groups. Users should not take its indications literally when the groups are not equal in size. Type "all" to obtain both indices. The indices are defined as:

calinski:

$(SSB/(K-1))/(SSW/(n-K))$, where $n$ is the number of data points and $K$ is the number of clusters. $SSW$ is the sum of squares within the clusters while $SSB$ is the sum of squares among the clusters. This index is simply an $F$ (ANOVA) statistic.

ssi:

the “Simple Structure Index” multiplicatively combines several elements which influence the interpretability of a partitioning solution. The best partition is indicated by the highest SSI value.

In a simulation study, Milligan and Cooper (1985) found that the Calinski-Harabasz criterion recovered the correct number of groups the most often. We recommend this criterion because, if the groups are of equal sizes, the maximum value of "calinski" usually indicates the correct number of groups. Another available index is the simple structure index "ssi". Users should not take the indications of these indices literally when the groups are not equal in size and explore the groups corresponding to other values of $K$.

Function cascadeKM has a plot method. Two plots are produced. The graph on the left has the objects in abscissa and the number of groups in ordinate. The groups are represented by colours. The graph on the right shows the values of the criterion ("calinski" or "ssi") for determining the best partition. The highest value of the criterion is marked in red. Points marked in orange, if any, indicate partitions producing an increase in the criterion value as the number of groups increases; they may represent other interesting partitions.

If sortg=TRUE, the objects are reordered by the following procedure: (1) a simple matching distance matrix is computed among the objects, based on the table of K-means assignments to groups, from $K$ = min.g to $K$ = max.g. (2) A principal coordinate analysis (PCoA, Gower 1966) is computed on the centred distance matrix. (3) The first principal coordinate is used as the new order of the objects in the graph. A simplified algorithm is used to compute the first principal coordinate only, using the iterative algorithm described in Legendre & Legendre (2012). The full distance matrix among objects is never computed; this avoids the problem of storing it when the number of objects is large. Distance values are computed as they are needed by the algorithm.

Value

Function cascadeKM returns an object of class cascadeKM with items:

partition	Table with the partitions found for different numbers of groups $K$, from $K$ = inf.gr to $K$ = sup.gr.
results	Values of the criterion to select the best partition.
criterion	The name of the criterion used.
size	The number of objects found in each group, for all partitions (columns).

Function cIndex returns a vector with the index values. The maximum value of these indices is supposed to indicate the best partition. These indices work best with groups of equal sizes. When the groups are not of equal sizes, one should not put too much faith in the maximum of these indices, and also explore the groups corresponding to other values of $K$.

Author(s)

Marie-Helene Ouellette [email protected], Sebastien Durand [email protected] and Pierre Legendre [email protected]. Parallel processing by Virgilio Gómez-Rubio. Edited for vegan by Jari Oksanen.

References

Calinski, T. and J. Harabasz. 1974. A dendrite method for cluster analysis. Commun. Stat. 3: 1–27.

Dolnicar, S., K. Grabler and J. A. Mazanec. 1999. A tale of three cities: perceptual charting for analyzing destination images. Pp. 39-62 in: Woodside, A. et al. [eds.] Consumer psychology of tourism, hospitality and leisure. CAB International, New York.

Gower, J. C. 1966. Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika 53: 325–338.

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

Milligan, G. W. & M. C. Cooper. 1985. An examination of procedures for determining the number of clusters in a data set. Psychometrika 50: 159–179.

Weingessel, A., Dimitriadou, A. and Dolnicar, S. 2002. An examination of indexes for determining the number of clusters in binary data sets. Psychometrika 67: 137–160.

See Also

kmeans.

Examples

# Partitioning a (10 x 10) data matrix of random numbers
 mat <- matrix(runif(100),10,10)
 res <- cascadeKM(mat, 2, 5, iter = 25, criterion = 'calinski') 
 toto <- plot(res)
 
 # Partitioning an autocorrelated time series
 vec <- sort(matrix(runif(30),30,1))
 res <- cascadeKM(vec, 2, 5, iter = 25, criterion = 'calinski')
 toto <- plot(res)
 
 # Partitioning a large autocorrelated time series
 # Note that we remove the grid lines
 vec <- sort(matrix(runif(1000),1000,1))
 res <- cascadeKM(vec, 2, 7, iter = 10, criterion = 'calinski')
 toto <- plot(res, gridcol=NA)

---

[Partial] [Constrained] Correspondence Analysis and Redundancy Analysis

Description

Function cca performs correspondence analysis, or optionally constrained correspondence analysis (a.k.a. canonical correspondence analysis), or optionally partial constrained correspondence analysis. Function rda performs redundancy analysis, or optionally principal components analysis. These are all very popular ordination techniques in community ecology.

Usage

## S3 method for class 'formula'
cca(formula, data, na.action = na.fail, subset = NULL,
  ...)
## S3 method for class 'formula'
rda(formula, data, scale=FALSE, na.action = na.fail,
  subset = NULL, ...)
## Default S3 method:
cca(X, Y, Z, ...)
## Default S3 method:
rda(X, Y, Z, scale=FALSE, ...)
ca(X, ...)
pca(X, scale=FALSE, ...)

Arguments

formula	Model formula, where the left hand side gives the community data matrix, right hand side gives the constraining variables, and conditioning variables can be given within a special function Condition.
data	Data frame containing the variables on the right hand side of the model formula.
X	Community data matrix.
Y	Constraining matrix, typically of environmental variables. Can be missing. If this is a data.frame, it will be expanded to a model.matrix where factors are expanded to contrasts (“dummy variables”). It is better to use formula instead of this argument, and some further analyses only work when formula was used.
Z	Conditioning matrix, the effect of which is removed (“partialled out”) before next step. Can be missing. If this is a data.frame, it is expanded similarly as constraining matrix.
scale	Scale species to unit variance (like correlations).
na.action	Handling of missing values in constraints or conditions. The default (na.fail) is to stop with missing value. Choice na.omit removes all rows with missing values. Choice na.exclude keeps all observations but gives NA for results that cannot be calculated. The WA scores of rows may be found also for missing values in constraints. Missing values are never allowed in dependent community data.
subset	Subset of data rows. This can be a logical vector which is TRUE for kept observations, or a logical expression which can contain variables in the working environment, data or species names of the community data.
...	Other arguments for print or plot functions (ignored in other functions). For pca() and ca(), arguments are passed to rda() and cca(), respectively.

Details

Since their introduction (ter Braak 1986), constrained, or canonical, correspondence analysis and its spin-off, redundancy analysis, have been the most popular ordination methods in community ecology. Functions cca and rda are similar to popular proprietary software Canoco, although the implementation is completely different. The functions are based on Legendre & Legendre's (2012) algorithm: in cca Chi-square transformed data matrix is subjected to weighted linear regression on constraining variables, and the fitted values are submitted to correspondence analysis performed via singular value decomposition (svd). Function rda is similar, but uses ordinary, unweighted linear regression and unweighted SVD. Legendre & Legendre (2012), Table 11.5 (p. 650) give a skeleton of the RDA algorithm of vegan. The algorithm of CCA is similar, but involves standardization by row and column weights.

The functions cca() and rda() can be called either with matrix-like entries for community data and constraints, or with formula interface. In general, the formula interface is preferred, because it allows a better control of the model and allows factor constraints. Some analyses of ordination results are only possible if model was fitted with formula (e.g., most cases of anova.cca, automatic model building).

In the following sections, X, Y and Z, although referred to as matrices, are more commonly data frames.

In the matrix interface, the community data matrix X must be given, but the other data matrices may be omitted, and the corresponding stage of analysis is skipped. If matrix Z is supplied, its effects are removed from the community matrix, and the residual matrix is submitted to the next stage. This is called partial correspondence or redundancy analysis. If matrix Y is supplied, it is used to constrain the ordination, resulting in constrained or canonical correspondence analysis, or redundancy analysis. Finally, the residual is submitted to ordinary correspondence analysis (or principal components analysis). If both matrices Z and Y are missing, the data matrix is analysed by ordinary correspondence analysis (or principal components analysis).

Instead of separate matrices, the model can be defined using a model formula. The left hand side must be the community data matrix (X). The right hand side defines the constraining model. The constraints can contain ordered or unordered factors, interactions among variables and functions of variables. The defined contrasts are honoured in factor variables. The constraints can also be matrices (but not data frames). The formula can include a special term Condition for conditioning variables (“covariables”) partialled out before analysis. So the following commands are equivalent: cca(X, Y, Z), cca(X ~ Y + Condition(Z)), where Y and Z refer to constraints and conditions matrices respectively.

Constrained correspondence analysis is indeed a constrained method: CCA does not try to display all variation in the data, but only the part that can be explained by the used constraints. Consequently, the results are strongly dependent on the set of constraints and their transformations or interactions among the constraints. The shotgun method is to use all environmental variables as constraints. However, such exploratory problems are better analysed with unconstrained methods such as correspondence analysis (decorana, corresp) or non-metric multidimensional scaling (metaMDS) and environmental interpretation after analysis (envfit, ordisurf). CCA is a good choice if the user has clear and strong a priori hypotheses on constraints and is not interested in the major structure in the data set.

CCA is able to correct the curve artefact commonly found in correspondence analysis by forcing the configuration into linear constraints. However, the curve artefact can be avoided only with a low number of constraints that do not have a curvilinear relation with each other. The curve can reappear even with two badly chosen constraints or a single factor. Although the formula interface makes it easy to include polynomial or interaction terms, such terms often produce curved artefacts (that are difficult to interpret), these should probably be avoided.

According to folklore, rda should be used with “short gradients” rather than cca. However, this is not based on research which finds methods based on Euclidean metric as uniformly weaker than those based on Chi-squared metric. However, standardized Euclidean distance may be an appropriate measures (see Hellinger standardization in decostand in particular).

Partial CCA (pCCA; or alternatively partial RDA) can be used to remove the effect of some conditioning or background or random variables or covariables before CCA proper. In fact, pCCA compares models cca(X ~ Z) and cca(X ~ Y + Z) and attributes their difference to the effect of Y cleansed of the effect of Z. Some people have used the method for extracting “components of variance” in CCA. However, if the effect of variables together is stronger than sum of both separately, this can increase total Chi-square after partialling out some variation, and give negative “components of variance”. In general, such components of “variance” are not to be trusted due to interactions between two sets of variables.

The unconstrained ordination methods, Principal Components Analysis (PCA) and Correspondence Analysis (CA), may be performed using pca() and ca(), which are simple wrappers around rda() and cca(), respectively. Functions pca() and ca() can only be called with matrix-like objects.

The functions have summary and plot methods which are documented separately (see plot.cca, summary.cca).

Value

Function cca returns a huge object of class cca, which is described separately in cca.object.

Function rda returns an object of class rda which inherits from class cca and is described in cca.object. The scaling used in rda scores is described in a separate vignette with this package.

Functions pca() and ca() return objects of class "vegan_pca" and "vegan_ca" respectively to avoid clashes with other packages. These classes inherit from "rda" and "cca" respectively.

Author(s)

The responsible author was Jari Oksanen, but the code borrows heavily from Dave Roberts (Montana State University, USA).

References

The original method was by ter Braak, but the current implementation follows Legendre and Legendre.

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

McCune, B. (1997) Influence of noisy environmental data on canonical correspondence analysis. Ecology 78, 2617-2623.

Palmer, M. W. (1993) Putting things in even better order: The advantages of canonical correspondence analysis. Ecology 74,2215-2230.

Ter Braak, C. J. F. (1986) Canonical Correspondence Analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology 67, 1167-1179.

See Also

This help page describes two constrained ordination functions, cca and rda and their corresponding unconstrained ordination functions, ca and pca. A related method, distance-based redundancy analysis (dbRDA) is described separately (capscale), as is dbRDA's unconstrained variant, principal coordinates analysis (PCO). All these functions return similar objects (described in cca.object). There are numerous support functions that can be used to access the result object. In the list below, functions of type cca will handle all three constrained ordination objects, and functions of rda only handle rda and capscale results.

The main plotting functions are plot.cca for all methods, and biplot.rda for RDA and dbRDA. However, generic vegan plotting functions can also handle the results. The scores can be accessed and scaled with scores.cca, and summarized with summary.cca. The eigenvalues can be accessed with eigenvals.cca and the regression coefficients for constraints with coef.cca. The eigenvalues can be plotted with screeplot.cca, and the (adjusted) $R^2$ can be found with RsquareAdj.rda. The scores can be also calculated for new data sets with predict.cca which allows adding points to ordinations. The values of constraints can be inferred from ordination and community composition with calibrate.cca.

Diagnostic statistics can be found with goodness.cca, inertcomp, spenvcor, intersetcor, tolerance.cca, and vif.cca. Function as.mlm.cca refits the result object as a multiple lm object, and this allows finding influence statistics (lm.influence, cooks.distance etc.).

Permutation based significance for the overall model, single constraining variables or axes can be found with anova.cca. Automatic model building with R step function is possible with deviance.cca, add1.cca and drop1.cca. Functions ordistep and ordiR2step (for RDA) are special functions for constrained ordination. Randomized data sets can be generated with simulate.cca.

Separate methods based on constrained ordination model are principal response curves (prc) and variance partitioning between several components (varpart).

Design decisions are explained in vignette on “Design decisions” which can be accessed with browseVignettes("vegan").

Examples

data(varespec)
data(varechem)
## Common but bad way: use all variables you happen to have in your
## environmental data matrix
vare.cca <- cca(varespec, varechem)
vare.cca
plot(vare.cca)
## Formula interface and a better model
vare.cca <- cca(varespec ~ Al + P*(K + Baresoil), data=varechem)
vare.cca
plot(vare.cca)
## Partialling out and negative components of variance
cca(varespec ~ Ca, varechem)
cca(varespec ~ Ca + Condition(pH), varechem)
## RDA
data(dune)
data(dune.env)
dune.Manure <- rda(dune ~ Manure, dune.env)
plot(dune.Manure)

---

Result Object from Constrained Ordination

Description

Ordination methods cca, rda, dbrda and capscale return similar result objects. All these methods use the same internal function ordConstrained. They differ only in (1) initial transformation of the data and in defining inertia, (2) weighting, and (3) the use of rectangular rows $\times$ columns data or symmetric rows $\times$ rows dissimilarities: rda initializes data to give variance or correlations as inertia, cca is based on double-standardized data to give Chi-square inertia and uses row and column weights, capscale maps the real part of dissimilarities to rectangular data and performs RDA, and dbrda performs an RDA-like analysis directly on symmetric dissimilarities.

Function ordConstrained returns the same result components for all these methods, and the calling function may add some more components to the final result. However, you should not access these result components directly (using $): the internal structure is not regarded as stable application interface (API), but it can change at any release. If you access the results components directly, you take a risk of breakage at any vegan release. The vegan provides a wide set of accessor functions to those components, and these functions are updated when the result object changes. This documentation gives an overview of accessor functions to the cca result object.

Usage

ordiYbar(x, model = c("CCA", "CA", "pCCA", "partial", "initial"))
## S3 method for class 'cca'
model.frame(formula, ...)
## S3 method for class 'cca'
model.matrix(object, ...)
## S3 method for class 'cca'
weights(object, display = "sites", ...)

Arguments

object, x, formula	A result object from cca, rda, dbrda, or capscale.
model	Show constrained ("CCA"), unconstrained ("CA") or conditioned “partial” ("pCCA") results. In ordiYbar the value can also be "initial" for the internal working input data, and "partial" for the internal working input data after removing the partial effects.
display	Display either "sites" or "species".
...	Other arguments passed to the the function.

Details

The internal (“working”) form of the dependent (community) data can be accessed with function ordiYbar. The form depends on the ordination method: for instance, in cca the data are weighted and Chi-square transformed, and in dbrda they are Gower-centred dissimilarities. The input data in the original (“response”) form can be accessed with fitted.cca and residuals.cca. Function predict.cca can return either working or response data, and also their lower-rank approximations.

The model matrix of independent data (“Constraints” and “Conditions”) can be extracted with model.matrix. In partial analysis, the function returns a list of design matrices called Conditions and Constraints. If either component was missing, a single matrix is returned. The redundant (aliased) terms do not appear in the model matrix. These terms can be found with alias.cca. Function model.frame tries to reconstruct the data frame from which the model matrices were derived. This is only possible if the original model was fitted with formula and data arguments, and still fails if the data are unavailable.

The number of observations can be accessed with nobs.cca, and the residual degrees of freedom with df.residual.cca. The information on observations with missing values can be accessed with na.action. The terms and formula of the fitted model can be accessed with formula and terms.

The weights used in cca can be accessed with weights. In unweighted methods (rda) all weights are equal.

The ordination results are saved in separate components for partial terms, constraints and residual unconstrained ordination. There is no guarantee that these components will have the same internal names as currently, and you should be cautious when developing scripts and functions that directly access these components.

The constrained ordination algorithm is based on QR decomposition of constraints and conditions (environmental data), and the QR component is saved separately for partial and constrained components. The QR decomposition of constraints can be accessed with qr.cca. This will also include the residual effects of partial terms (Conditions), and it should be used together with ordiYbar(x, "partial"). The environmental data are first centred in rda or weighted and centred in cca. The QR decomposition is used in many functions that access cca results, and it can be used to find many items that are not directly stored in the object. For examples, see coef.cca, coef.rda, vif.cca, permutest.cca, predict.cca, predict.rda, calibrate.cca. See qr for other possible uses of this component. For instance, the rank of the constraints can be found from the QR decomposition.

The eigenvalues of the solution can be accessed with eigenvals.cca. Eigenvalues are not evaluated for partial component, and they will only be available for constrained and residual components.

The ordination scores are internally stored as (weighted) orthonormal scores matrices. These results can be accessed with scores.cca and scores.rda functions. The ordination scores are scaled when accessed with scores functions, but internal (weighted) orthonormal scores can be accessed by setting scaling = FALSE. Unconstrained residual component has species and site scores, and constrained component has also fitted site scores or linear combination scores for sites and biplot scores and centroids for constraint variables. The biplot scores correspond to the model.matrix, and centroids are calculated for factor variables when they were used. The scores can be selected by defining the axes, and there is no direct way of accessing all scores of a certain component. The number of dimensions can be assessed from eigenvals. In addition, some other types can be derived from the results although not saved in the results. For instance, regression scores and model coefficients can be accessed with scores and coef functions. Partial component will have no scores.

Distance-based methods (dbrda, capscale) can have negative eigenvalues and associated imaginary axis scores. In addition, species scores are initially missing in dbrda and they are accessory and found after analysis in capscale (and may be misleading). Function sppscores can be used to add species scores or replace them with more meaningful ones.

Note

The latest large change of result object was made in release 2.5-1 in 2016. You can modernize ancient stray results with modernobject <- update(ancientobject).

Author(s)

Jari Oksanen

References

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

See Also

The core function is ordConstrained which is called by cca, rda, dbrda and capscale. The basic class is "cca" for all methods, and the following functions are defined for this class: RsquareAdj.cca, SSD.cca, add1.cca, alias.cca, anova.cca, as.mlm.cca, biplot.cca, bstick.cca, calibrate.cca, centredist.cca, coef.cca, cooks.distance.cca, deviance.cca, df.residual.cca, drop1.cca, eigenvals.cca, extractAIC.cca, fitted.cca, goodness.cca, hatvalues.cca, labels.cca, model.frame.cca, model.matrix.cca, nobs.cca, permutest.cca, plot.cca, points.cca, predict.cca, print.cca, qr.cca, residuals.cca, rstandard.cca, rstudent.cca, scores.cca, screeplot.cca, sigma.cca, simulate.cca, stressplot.cca, summary.cca, text.cca, tolerance.cca, vcov.cca, weights.cca. Other functions handling "cca" objects include inertcomp, intersetcor, mso, ordiresids, ordistep and ordiR2step.

---

Canonical Correlation Analysis

Description

Canonical correlation analysis, following Brian McArdle's unpublished graduate course notes, plus improvements to allow the calculations in the case of very sparse and collinear matrices, and permutation test of Pillai's trace statistic.

Usage

CCorA(Y, X, stand.Y=FALSE, stand.X=FALSE, permutations = 0, ...)

## S3 method for class 'CCorA'
biplot(x, plot.type="ov", xlabs, plot.axes = 1:2, int=0.5, 
   col.Y="red", col.X="blue", cex=c(0.7,0.9), ...)

Arguments

Y	Left matrix (object class: matrix or data.frame).
X	Right matrix (object class: matrix or data.frame).
stand.Y	Logical; should Y be standardized?
stand.X	Logical; should X be standardized?
permutations	a list of control values for the permutations as returned by the function how, or the number of permutations required, or a permutation matrix where each row gives the permuted indices.
x	CCoaR result object.
plot.type	A character string indicating which of the following plots should be produced: "objects", "variables", "ov" (separate graphs for objects and variables), or "biplots". Any unambiguous subset containing the first letters of these names can be used instead of the full names.
xlabs	Row labels. The default is to use row names, NULL uses row numbers instead, and NA suppresses plotting row names completely.
plot.axes	A vector with 2 values containing the order numbers of the canonical axes to be plotted. Default: first two axes.
int	Radius of the inner circles plotted as visual references in the plots of the variables. Default: int=0.5. With int=0, no inner circle is plotted.
col.Y	Color used for objects and variables in the first data table (Y) plots. In biplots, the objects are in black.
col.X	Color used for objects and variables in the second data table (X) plots.
cex	A vector with 2 values containing the size reduction factors for the object and variable names, respectively, in the plots. Default values: cex=c(0.7,0.9).
...	Other arguments passed to these functions. The function biplot.CCorA passes graphical arguments to biplot and biplot.default. CCorA currently ignores extra arguments.

Details

Canonical correlation analysis (Hotelling 1936) seeks linear combinations of the variables of Y that are maximally correlated to linear combinations of the variables of X. The analysis estimates the relationships and displays them in graphs. Pillai's trace statistic is computed and tested parametrically (F-test); a permutation test is also available.

Algorithmic note – The blunt approach would be to read the two matrices, compute the covariance matrices, then the matrix S12 %*% inv(S22) %*% t(S12) %*% inv(S11). Its trace is Pillai's trace statistic. This approach may fail, however, when there is heavy multicollinearity in very sparse data matrices. The safe approach is to replace all data matrices by their PCA object scores.

The function can produce different types of plots depending on the option chosen: "objects" produces two plots of the objects, one in the space of Y, the second in the space of X; "variables" produces two plots of the variables, one of the variables of Y in the space of Y, the second of the variables of X in the space of X; "ov" produces four plots, two of the objects and two of the variables; "biplots" produces two biplots, one for the first matrix (Y) and one for second matrix (X) solutions. For biplots, the function passes all arguments to biplot.default; consult its help page for configuring biplots.

Value

Function CCorA returns a list containing the following elements:

Pillai	Pillai's trace statistic = sum of the canonical eigenvalues.
Eigenvalues	Canonical eigenvalues. They are the squares of the canonical correlations.
CanCorr	Canonical correlations.
Mat.ranks	Ranks of matrices Y and X.
RDA.Rsquares	Bimultivariate redundancy coefficients (R-squares) of RDAs of Y|X and X|Y.
RDA.adj.Rsq	RDA.Rsquares adjusted for n and the number of explanatory variables.
nperm	Number of permutations.
p.Pillai	Parametric probability value associated with Pillai's trace.
p.perm	Permutational probability associated with Pillai's trace.
Cy	Object scores in Y biplot.
Cx	Object scores in X biplot.
corr.Y.Cy	Scores of Y variables in Y biplot, computed as cor(Y,Cy).
corr.X.Cx	Scores of X variables in X biplot, computed as cor(X,Cx).
corr.Y.Cx	cor(Y,Cy) available for plotting variables Y in space of X manually.
corr.X.Cy	cor(X,Cx) available for plotting variables X in space of Y manually.
control	A list of control values for the permutations as returned by the function how.
call	Call to the CCorA function.

Author(s)

Pierre Legendre, Departement de Sciences Biologiques, Universite de Montreal. Implemented in vegan with the help of Jari Oksanen.

References

Hotelling, H. 1936. Relations between two sets of variates. Biometrika 28: 321-377.

Legendre, P. 2005. Species associations: the Kendall coefficient of concordance revisited. Journal of Agricultural, Biological, and Environmental Statistics 10: 226-245.

Examples

# Example using two mite groups. The mite data are available in vegan
data(mite)
# Two mite species associations (Legendre 2005, Fig. 4)
group.1 <- c(1,2,4:8,10:15,17,19:22,24,26:30)
group.2 <- c(3,9,16,18,23,25,31:35)
# Separate Hellinger transformations of the two groups of species 
mite.hel.1 <- decostand(mite[,group.1], "hel")
mite.hel.2 <- decostand(mite[,group.2], "hel")
rownames(mite.hel.1) = paste("S",1:nrow(mite),sep="")
rownames(mite.hel.2) = paste("S",1:nrow(mite),sep="")
out <- CCorA(mite.hel.1, mite.hel.2)
out
biplot(out, "ob")                 # Two plots of objects
biplot(out, "v", cex=c(0.7,0.6))  # Two plots of variables
biplot(out, "ov", cex=c(0.7,0.6)) # Four plots (2 for objects, 2 for variables)
biplot(out, "b", cex=c(0.7,0.6))  # Two biplots
biplot(out, xlabs = NA, plot.axes = c(3,5))    # Plot axes 3, 5. No object names
biplot(out, plot.type="biplots", xlabs = NULL) # Replace object names by numbers

# Example using random numbers. No significant relationship is expected
mat1 <- matrix(rnorm(60),20,3)
mat2 <- matrix(rnorm(100),20,5)
out2 = CCorA(mat1, mat2, permutations=99)
out2
biplot(out2, "b")

---

Distances of Points to Class Centroids

Description

Function finds Euclidean or squared Mahalanobis distances of points to class centroids in betadisper or ordination results.

Usage

## Default S3 method:
centredist(x, group, distance = c("euclidean",
    "mahalanobis"), display = "sites", w, ...)
## S3 method for class 'betadisper'
centredist(x, ...)
## S3 method for class 'cca'
centredist(x, group, distance = c("euclidean", "mahalanobis"),
    display = c("sites", "species"), rank = 2, ...)
## S3 method for class 'rda'
centredist(x, group, distance = c("euclidean", "mahalanobis"),
    display = "sites", rank = 2, ...)
## S3 method for class 'dbrda'
centredist(x, group, distance = c("euclidean","mahalanobis"),
    display = "sites", rank = 2, ...)
## S3 method for class 'wcmdscale'
centredist(x, group, distance= c("euclidean","mahalanobis"),
    display = "sites", rank = 2, ...)

Arguments

x	An ordination result object, or a matrix for the default method or a betadisper result object.
group	Factor to define classes for which centroids are found.
distance	Use either simple Euclidean distance or squared Mahalanobis distances that take into account the shape of covariance ellipses.
display	Kind of scores.
w	Weights. If x has weights (such as cca result), these will be used.
rank	Number of axes used in ordination methods.
...	Other arguments to functions (ignored).

Details

Function finds either simple Euclidean distances of all points to their centroid as given by argument group or squared Mahalanobis distances of points to the centroid. In addition to the distances, it returns information of the original group and the nearest group for each point.

Currently the distances are scaled so that the sum of squared Euclidean distance will give the residual (unexplained) eigenvalue in constrained ordination with given group. The distances may not correspond to the distances in ordination plots with specific scaling and setting argument scaling gives an error. Due to restriction of this scaling, Euclidean distances for species scores are only allowed in ca and cca. These restrictions may be relaxed in the future versions of the method.

Squared Mahalanobis distances (mahalanobis) are scaled by the within group covariance matrix, and they are consistent with ordiellipse.

Mahalanobis distances can only be calculated for group that are larger than rank. For rank = 2, minimum permissible group has three points. However, even this can fail in degenerate cases, such as three points on a line. Euclidean distances can be calculated for any rank. With rank = "full" distance-based methods wcmdscale, pco, dbrda and betadisper will handle negative eigenvalues and associated eigenvectors, but with numeric rank they will only use real eigenvectors with positive eigenvalues.

There are specific methods for some ordination models, but many will be handled by the default method (including metaMDS and monoMDS).

The specific method for betadisper ignores most arguments and calculates the distances as specified in model fitting. The betadisper result with argument type = "centroid" is equal to wcmdscale, pco or dbrda with the same group and rank = "full".

Value

Function returns a list with components

centre	Centre to which a point was originally allocated.
nearest	Nearest centre to a point.
distances	Distance matrix where each observation point is a row, and column gives the distance of the point to that centre.

Note

centerdist is a synonym of centredist.

Author(s)

Jari Oksanen

See Also

betadisper, goodness.cca, mahalanobis and various ordination methods.

Examples

data(dune, dune.env)
  d <- vegdist(dune) # Bray-Curtis
  mod <- with(dune.env, betadisper(d, Management))
  centredist(mod)
  ## equal to metric scaling when using centroids instead of medians
  modc <- with(dune.env, betadisper(d, Management, type = "centroid"))
  ord <- pco(d)
  all.equal(centredist(modc), centredist(ord, dune.env$Management, rank="full"),
     check.attributes = FALSE)

---

Multinomial Species Classification Method (CLAM)

Description

The CLAM statistical approach for classifying generalists and specialists in two distinct habitats is described in Chazdon et al. (2011).

Usage

clamtest(comm, groups, coverage.limit = 10, specialization = 2/3, 
   npoints = 20, alpha = 0.05/20)
## S3 method for class 'clamtest'
summary(object, ...)
## S3 method for class 'clamtest'
plot(x, xlab, ylab, main,  pch = 21:24, col.points = 1:4, 
   col.lines = 2:4, lty = 1:3, position = "bottomright", ...)

Arguments

comm	Community matrix, consisting of counts.
groups	A vector identifying the two habitats. Must have exactly two unique values or levels. Habitat IDs in the grouping vector must match corresponding rows in the community matrix comm.
coverage.limit	Integer, the sample coverage based correction is applied to rare species with counts below this limit. Sample coverage is calculated separately for the two habitats. Sample relative abundances are used for species with higher than or equal to coverage.limit total counts per habitat.
specialization	Numeric, specialization threshold value between 0 and 1. The value of $2/3$ represents ‘supermajority’ rule, while a value of $1/2$ represents a ‘simple majority’ rule to assign shared species as habitat specialists.
npoints	Integer, number of points used to determine the boundary lines in the plots.
alpha	Numeric, nominal significance level for individual tests. The default value reduces the conventional limit of $0.05$ to account for overdispersion and multiple testing for several species simultaneously. However, the is no firm reason for exactly this limit.
x, object	Fitted model object of class "clamtest".
xlab, ylab	Labels for the plot axes.
main	Main title of the plot.
pch, col.points	Symbols and colors used in plotting species groups.
lty, col.lines	Line types and colors for boundary lines in plot to separate species groups.
position	Position of figure legend, see legend for specification details. Legend not shown if position = NULL.
...	Additional arguments passed to methods.

Details

The method uses a multinomial model based on estimated species relative abundance in two habitats (A, B). It minimizes bias due to differences in sampling intensities between two habitat types as well as bias due to insufficient sampling within each habitat. The method permits a robust statistical classification of habitat specialists and generalists, without excluding rare species a priori (Chazdon et al. 2011). Based on a user-defined specialization threshold, the model classifies species into one of four groups: (1) generalists; (2) habitat A specialists; (3) habitat B specialists; and (4) too rare to classify with confidence.

Value

A data frame (with class attribute "clamtest"), with columns:

Species:	species name (column names from comm),
Total_*A*:	total count in habitat A,
Total_*B*:	total count in habitat B,
Classes:	species classification, a factor with levels Generalist, Specialist_*A*, Specialist_*B*, and Too_rare.

*A* and *B* are placeholders for habitat names/labels found in the data.

The summary method returns descriptive statistics of the results. The plot method returns values invisibly and produces a bivariate scatterplot of species total abundances in the two habitats. Symbols and boundary lines are shown for species groups.

Note

The code was tested against standalone CLAM software provided on the website of Anne Chao (which were then at http://chao.stat.nthu.edu.tw/wordpress); minor inconsistencies were found, especially for finding the threshold for 'too rare' species. These inconsistencies are probably due to numerical differences between the two implementation. The current R implementation uses root finding for iso-lines instead of iterative search.

The original method (Chazdon et al. 2011) has two major problems:

1. It assumes that the error distribution is multinomial. This is a justified choice if individuals are freely distributed, and there is no over-dispersion or clustering of individuals. In most ecological data, the variance is much higher than multinomial assumption, and therefore test statistic are too optimistic.

2. The original authors suggest that multiple testing adjustment for multiple testing should be based on the number of points (npoints) used to draw the critical lines on the plot, whereas the adjustment should be based on the number of tests (i.e., tested species). The function uses the same numerical values as the original paper, but there is no automatic connection between npoints and alpha arguments, but you must work out the adjustment yourself.

Author(s)

Peter Solymos [email protected]

References

Chazdon, R. L., Chao, A., Colwell, R. K., Lin, S.-Y., Norden, N., Letcher, S. G., Clark, D. B., Finegan, B. and Arroyo J. P.(2011). A novel statistical method for classifying habitat generalists and specialists. Ecology 92, 1332–1343.

Examples

data(mite)
data(mite.env)
sol <- with(mite.env, clamtest(mite, Shrub=="None", alpha=0.005))
summary(sol)
head(sol)
plot(sol)

---

Create an Object for Null Model Algorithms

Description

The commsim function can be used to feed Null Model algorithms into nullmodel analysis. The make.commsim function returns various predefined algorithm types (see Details). These functions represent low level interface for community null model infrastructure in vegan with the intent of extensibility, and less emphasis on direct use by users.

Usage

commsim(method, fun, binary, isSeq, mode)
make.commsim(method)
## S3 method for class 'commsim'
print(x, ...)

Arguments

method	Character, name of the algorithm.
fun	A function. For possible formal arguments of this function see Details.
binary	Logical, if the algorithm applies to presence-absence or count matrices.
isSeq	Logical, if the algorithm is sequential (needs burnin and thinning) or not.
mode	Character, storage mode of the community matrix, either "integer" or "double".
x	An object of class commsim.
...	Additional arguments.

Details

The function fun must return an array of dim(nr, nc, n), and must take some of the following arguments:

x:

input matrix,

n:

number of permuted matrices in output,

nr:

number of rows,

nc:

number of columns,

rs:

vector of row sums,

cs:

vector of column sums,

rf:

vector of row frequencies (non-zero cells),

cf:

vector of column frequencies (non-zero cells),

s:

total sum of x,

fill:

matrix fill (non-zero cells),

thin:

thinning value for sequential algorithms,

...:

additional arguments.

You can define your own null model, but several null model algorithm are pre-defined and can be called by their name. The predefined algorithms are described in detail in the following chapters. The binary null models produce matrices of zeros (absences) and ones (presences) also when input matrix is quantitative. There are two types of quantitative data: Counts are integers with a natural unit so that individuals can be shuffled, but abundances can have real (floating point) values and do not have a natural subunit for shuffling. All quantitative models can handle counts, but only some are able to handle real values. Some of the null models are sequential so that the next matrix is derived from the current one. This makes models dependent from previous models, and usually you must thin these matrices and study the sequences for stability: see oecosimu for details and instructions.

See Examples for structural constraints imposed by each algorithm and defining your own null model.

Value

An object of class commsim with elements corresponding to the arguments (method, binary, isSeq, mode, fun).

If the input of make.comsimm is a commsim object, it is returned without further evaluation. If this is not the case, the character method argument is matched against predefined algorithm names. An error message is issued if none such is found. If the method argument is missing, the function returns names of all currently available null model algorithms as a character vector.

Binary null models

All binary null models preserve fill: number of presences or conversely the number of absences. The classic models may also preserve column (species) frequencies (c0) or row frequencies or species richness of each site (r0) and take into account commonness and rarity of species (r1, r2). Algorithms swap, tswap, curveball, quasiswap and backtracking preserve both row and column frequencies. Three first ones are sequential but the two latter are non-sequential and produce independent matrices. Basic algorithms are reviewed by Wright et al. (1998).

"r00":

non-sequential algorithm for binary matrices that only preserves the number of presences (fill).

"r0":

non-sequential algorithm for binary matrices that preserves the site (row) frequencies.

"r1":

non-sequential algorithm for binary matrices that preserves the site (row) frequencies, but uses column marginal frequencies as probabilities of selecting species.

"r2":

non-sequential algorithm for binary matrices that preserves the site (row) frequencies, and uses squared column marginal frequencies as as probabilities of selecting species.

"c0":

non-sequential algorithm for binary matrices that preserves species frequencies (Jonsson 2001).

"swap":

sequential algorithm for binary matrices that changes the matrix structure, but does not influence marginal sums (Gotelli & Entsminger 2003). This inspects $2 \times 2$ submatrices so long that a swap can be done.

"tswap":

sequential algorithm for binary matrices. Same as the "swap" algorithm, but it tries a fixed number of times and performs zero to many swaps at one step (according to the thin argument in the call). This approach was suggested by Miklós & Podani (2004) because they found that ordinary swap may lead to biased sequences, since some columns or rows are more easily swapped.

"curveball":

sequential method for binary matrices that implements the ‘Curveball’ algorithm of Strona et al. (2014). The algorithm selects two random rows and finds the set of unique species that occur only in one of these rows. The algorithm distributes the set of unique species to rows preserving the original row frequencies. Zero to several species are swapped in one step, and usually the matrix is perturbed more strongly than in other sequential methods.

"quasiswap":

non-sequential algorithm for binary matrices that implements a method where matrix is first filled honouring row and column totals, but with integers that may be larger than one. Then the method inspects random $2 \times 2$ matrices and performs a quasiswap on them. In addition to ordinary swaps, quasiswap can reduce numbers above one to ones preserving marginal totals (Miklós & Podani 2004). The method is non-sequential, but it accepts thin argument: the convergence is checked at every thin steps. This allows performing several ordinary swaps in addition to fill changing swaps which helps in reducing or removing the bias.

"greedyqswap":

A greedy variant of quasiswap. In greedy step, one element of the $2 \times 2$ matrix is taken from $> 1$ elements. The greedy steps are biased, but the method can be thinned, and only the first of thin steps is greedy. Even modest thinning (say thin = 20) removes or reduces the bias, and thin = 100 (1% greedy steps) looks completely safe and still speeds up simulation. The code is experimental and it is provided here for further scrutiny, and should be tested for bias before use.

"backtracking":

non-sequential algorithm for binary matrices that implements a filling method with constraints both for row and column frequencies (Gotelli & Entsminger 2001). The matrix is first filled randomly, but typically row and column sums are reached before all incidences are filled in. After this begins "backtracking", where some of the incidences are removed, and filling is started again, and this backtracking is done so many times that all incidences will be filled into matrix. The results may be biased and should be inspected carefully before use.

Quantitative Models for Counts with Fixed Marginal Sums

These models shuffle individuals of counts and keep marginal sums fixed, but marginal frequencies are not preserved. Algorithm r2dtable uses standard R function r2dtable also used for simulated $P$-values in chisq.test. Algorithm quasiswap_count uses the same, but preserves the original fill. Typically this means increasing numbers of zero cells and the result is zero-inflated with respect to r2dtable.

"r2dtable":

non-sequential algorithm for count matrices. This algorithm keeps matrix sum and row/column sums constant. Based on r2dtable.

"quasiswap_count":

non-sequential algorithm for count matrices. This algorithm is similar as Carsten Dormann's swap.web function in the package bipartite. First, a random matrix is generated by the r2dtable function preserving row and column sums. Then the original matrix fill is reconstructed by sequential steps to increase or decrease matrix fill in the random matrix. These steps are based on swapping $2 \times 2$ submatrices (see "swap_count" algorithm for details) to maintain row and column totals.

Quantitative Swap Models

Quantitative swap models are similar to binary swap, but they swap the largest permissible value. The models in this section all maintain the fill and perform a quantitative swap only if this can be done without changing the fill. Single step of swap often changes the matrix very little. In particular, if cell counts are variable, high values change very slowly. Checking the chain stability and independence is even more crucial than in binary swap, and very strong thinning is often needed. These models should never be used without inspecting their properties for the current data. These null models can also be defined using permatswap function.

"swap_count":

sequential algorithm for count matrices. This algorithm find $2 \times 2$ submatrices that can be swapped leaving column and row totals and fill unchanged. The algorithm finds the largest value in the submatrix that can be swapped ($d$). Swap means that the values in diagonal or antidiagonal positions are decreased by $d$, while remaining cells are increased by $d$. A swap is made only if fill does not change.

"abuswap_r":

sequential algorithm for count or nonnegative real valued matrices with fixed row frequencies (see also permatswap). The algorithm is similar to swap_count, but uses different swap value for each row of the $2 \times 2$ submatrix. Each step changes the the corresponding column sums, but honours matrix fill, row sums, and row/column frequencies (Hardy 2008; randomization scheme 2x).

"abuswap_c":

sequential algorithm for count or nonnegative real valued matrices with fixed column frequencies (see also permatswap). The algorithm is similar as the previous one, but operates on columns. Each step changes the the corresponding row sums, but honours matrix fill, column sums, and row/column frequencies (Hardy 2008; randomization scheme 3x).

Quantitative Swap and Shuffle Models

Quantitative Swap and Shuffle methods (swsh methods) preserve fill and column and row frequencies, and also either row or column sums. The methods first perform a binary quasiswap and then shuffle original quantitative data to non-zero cells. The samp methods shuffle original non-zero cell values and can be used also with non-integer data. The both methods redistribute individuals randomly among non-zero cells and can only be used with integer data. The shuffling is either free over the whole matrix, or within rows (r methods) or within columns (c methods). Shuffling within a row preserves row sums, and shuffling within a column preserves column sums. These models can also be defined with permatswap.

"swsh_samp":

non-sequential algorithm for quantitative data (either integer counts or non-integer values). Original non-zero values values are shuffled.

"swsh_both":

non-sequential algorithm for count data. Individuals are shuffled freely over non-zero cells.

"swsh_samp_r":

non-sequential algorithm for quantitative data. Non-zero values (samples) are shuffled separately for each row.

"swsh_samp_c":

non-sequential algorithm for quantitative data. Non-zero values (samples) are shuffled separately for each column.

"swsh_both_r":

non-sequential algorithm for count matrices. Individuals are shuffled freely for non-zero values within each row.

"swsh_both_c":

non-sequential algorithm for count matrices. Individuals are shuffled freely for non-zero values with each column.

Quantitative Shuffle Methods

Quantitative shuffle methods are generalizations of binary models r00, r0 and c0. The _ind methods shuffle individuals so that the grand sum, row sum or column sums are preserved. These methods are similar as r2dtable but with still slacker constraints on marginal sums. The _samp and _both methods first apply the corresponding binary model with similar restriction on marginal frequencies and then distribute quantitative values over non-zero cells. The _samp models shuffle original cell values and can therefore handle also non-count real values. The _both models shuffle individuals among non-zero values. The shuffling is over the whole matrix in r00_, and within row in r0_ and within column in c0_ in all cases.

"r00_ind":

non-sequential algorithm for count matrices. This algorithm preserves grand sum and individuals are shuffled among cells of the matrix.

"r0_ind":

non-sequential algorithm for count matrices. This algorithm preserves row sums and individuals are shuffled among cells of each row of the matrix.

"c0_ind":

non-sequential algorithm for count matrices. This algorithm preserves column sums and individuals are shuffled among cells of each column of the matrix.

"r00_samp":

non-sequential algorithm for count or nonnegative real valued (mode = "double") matrices. This algorithm preserves grand sum and cells of the matrix are shuffled.

"r0_samp":

non-sequential algorithm for count or nonnegative real valued (mode = "double") matrices. This algorithm preserves row sums and cells within each row are shuffled.

"c0_samp":

non-sequential algorithm for count or nonnegative real valued (mode = "double") matrices. This algorithm preserves column sums constant and cells within each column are shuffled.

"r00_both":

non-sequential algorithm for count matrices. This algorithm preserves grand sum and cells and individuals among cells of the matrix are shuffled.

"r0_both":

non-sequential algorithm for count matrices. This algorithm preserves grand sum and cells and individuals among cells of each row are shuffled.

"c0_both":

non-sequential algorithm for count matrices. This algorithm preserves grand sum and cells and individuals among cells of each column are shuffled.

Author(s)

Jari Oksanen and Peter Solymos

References

Gotelli, N.J. & Entsminger, N.J. (2001). Swap and fill algorithms in null model analysis: rethinking the knight's tour. Oecologia 129, 281–291.

Gotelli, N.J. & Entsminger, N.J. (2003). Swap algorithms in null model analysis. Ecology 84, 532–535.

Hardy, O. J. (2008) Testing the spatial phylogenetic structure of local communities: statistical performances of different null models and test statistics on a locally neutral community. Journal of Ecology 96, 914–926.

Jonsson, B.G. (2001) A null model for randomization tests of nestedness in species assemblages. Oecologia 127, 309–313.

Miklós, I. & Podani, J. (2004). Randomization of presence-absence matrices: comments and new algorithms. Ecology 85, 86–92.

Patefield, W. M. (1981) Algorithm AS159. An efficient method of generating r x c tables with given row and column totals. Applied Statistics 30, 91–97.

Strona, G., Nappo, D., Boccacci, F., Fattorini, S. & San-Miguel-Ayanz, J. (2014). A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals. Nature Communications 5:4114 doi:10.1038/ncomms5114.

Wright, D.H., Patterson, B.D., Mikkelson, G.M., Cutler, A. & Atmar, W. (1998). A comparative analysis of nested subset patterns of species composition. Oecologia 113, 1–20.

See Also

See permatfull, permatswap for alternative specification of quantitative null models. Function oecosimu gives a higher-level interface for applying null models in hypothesis testing and analysis of models. Function nullmodel and simulate.nullmodel are used to generate arrays of simulated null model matrices.

Examples

## write the r00 algorithm
f <- function(x, n, ...)
    array(replicate(n, sample(x)), c(dim(x), n))
(cs <- commsim("r00", fun=f, binary=TRUE,
    isSeq=FALSE, mode="integer"))

## retrieving the sequential swap algorithm
(cs <- make.commsim("swap"))

## feeding a commsim object as argument
make.commsim(cs)

## making the missing c1 model using r1 as a template
##   non-sequential algorithm for binary matrices
##   that preserves the species (column) frequencies,
##   but uses row marginal frequencies
##   as probabilities of selecting sites
f <- function (x, n, nr, nc, rs, cs, ...) {
    out <- array(0L, c(nr, nc, n))
    J <- seq_len(nc)
    storage.mode(rs) <- "double"
    for (k in seq_len(n))
        for (j in J)
            out[sample.int(nr, cs[j], prob = rs), j, k] <- 1L
    out
}
cs <- make.commsim("r1")
cs$method <- "c1"
cs$fun <- f

## structural constraints
diagfun <- function(x, y) {
    c(sum = sum(y) == sum(x),
        fill = sum(y > 0) == sum(x > 0),
        rowSums = all(rowSums(y) == rowSums(x)),
        colSums = all(colSums(y) == colSums(x)),
        rowFreq = all(rowSums(y > 0) == rowSums(x > 0)),
        colFreq = all(colSums(y > 0) == colSums(x > 0)))
}
evalfun <- function(meth, x, n) {
    m <- nullmodel(x, meth)
    y <- simulate(m, nsim=n)
    out <- rowMeans(sapply(1:dim(y)[3],
        function(i) diagfun(attr(y, "data"), y[,,i])))
    z <- as.numeric(c(attr(y, "binary"), attr(y, "isSeq"),
        attr(y, "mode") == "double"))
    names(z) <- c("binary", "isSeq", "double")
    c(z, out)
}
x <- matrix(rbinom(10*12, 1, 0.5)*rpois(10*12, 3), 12, 10)
algos <- make.commsim()
a <- t(sapply(algos, evalfun, x=x, n=10))
print(as.table(ifelse(a==1,1,0)), zero.print = ".")

---