Contents:

• Introduction
• ReDundancy Analysis (RDA)
• Darlingtonia Example (RDA)
• Distance-based ReDundancy Analysis (dbRDA)
• Oak Example (dbRDA)
• Interpreting RDA and dbRDA
• Conclusions
• References

Introduction

Other than CCA, all of the ordination techniques that we’ve discussed thus far are unconstrained or indirect gradient analyses – the only information that is incorporated into the procedure is the matrix of response variables.  Constrained ordinations are direct gradient analyses.  These are analogous to regression, where the objective is to correlate the response matrix and one or more explanatory variables.

An example might help here.  Imagine the sample units as a 3-D point cloud.  We could view this 3-D point cloud from different 2-D vantage points.  For example, we could view it from the perspective from which it is as wide and as tall as possible – this is what an unconstrained ordination does by showing as much variation as possible.  Now imagine that the points are also colored to reflect different groups.  Depending on where we view the point cloud from, these groups might appear more or less discrete.  A constrained ordination shows us the view from the spot where the groups are as distinct as possible.  These different perspectives are illustrated in the chapter comparing ordination techniques.

In summary, different types of analyses provide different ways to view the point cloud:

• Unconstrained ordinations identify the perspective that shows as much of the variation in the point cloud as possible.
• Constrained ordinations view or interpret the point cloud from the perspective of the chosen explanatory variable(s).

We will consider two similar types of constrained ordinations, ReDundancy Analysis (RDA) and distance-based RDA (dbRDA).  RDA is the standard approach, and dbRDA involves an additional preliminary step to ensure that the distance matrix can be handled by RDA.  The steps are illustrated in the below figure.

ReDundancy Analysis (RDA)

Approach

Legendre & Legendre (2012; Section 11.1) describe ReDundancy Analysis (RDA) as “the direct extension of multiple regression to the modelling of multivariate response data” (p. 629).  It is based on Euclidean distances. Transformations and standardizations can be conducted at the user’s discretion before the analysis.

RDA uses two common techniques in sequence: “RDA is a multivariate (meaning multiresponse) multiple linear regression followed by a PCA of the matrix of fitted values” (Borcard et al. 2018, p. 205).

1. Regression: Each variable in the response matrix is fit to the variables in the explanatory matrix.  The regression model assumes that the explanatory and response variables are linearly related to one another.  This may require transformations of variables.  The fitted values of each response variable with each explanatory variable are calculated.  The residuals (i.e., variation not accounted for by the explanatory variables) are also calculated – these are simply the original values minus the fitted values.

2. PCA: Two separate PCAs are conducted:
   1. A PCA on the matrix of fitted values. By conducting the PCA on the fitted values rather than the original values, the resulting eigenvalues and eigenvectors are constrained to be linear combinations of the explanatory variables.
   2. A PCA on the matrix of residuals.  These residuals are conditional: they are the variation that was not explained by the explanatory variables.

Notes:

• If no explanatory variables are provided, then the first step is meaningless: there are no fitted values, and the second step results in a regular PCA.
• The explanatory variables should not be highly correlated with one another. If two such variables are highly correlated, consider omitting one or using PCA first to reduce them to a single principal component.

The result of this process is a set of principal components conditioned (or weighted) by the explanatory variables.  “In RDA, one can truly say that the axes explain or model (in the statistical sense) the variation of the dependent [response] matrix” (Borcard et al. 2018, p. 205).

The degrees of freedom (df) required for the explanatory variables indicate how many axes are constrained by those variables.  As usual, a single continuous variable requires 1 df, and a factor with m levels requires m-1 df.  The axes that are constrained by these variables are referred to as canonical axes or constrained axes.  As with a regular PCA, each of these axes is orthogonal to the others, and can be tested and evaluated individually.   The eigenvalues associated with these axes measure the amount of variance explained by the RDA model (Borcard et al. 2018).

Once the constrained axes have accounted for all of the variation that they can, the remaining variation is arrayed along unconstrained axes.  These are interpreted in the same way as in a PCA, except that they are explaining the variation that remains after having accounted for as much variation as possible through the constrained axes.   The eigenvalues associated with these axes measure the amount of variance represented by these axes, but not explained by any included explanatory variable (Borcard et al. 2018).

We have noted several times that PCs are always a declining function: the first explains as much variation as it can, the second as much as it can of the remaining, etc.  Constrained ordinations are a partial exception to this rule.  The constrained axes are a declining function and the unconstrained axes are also a declining function, but the transition from the last constrained axis to the first unconstrained axis may be associated with a (sometimes dramatic) increase in the magnitude of the eigenvalue.

A partial RDA allows the user to relate a set of responses to one or more explanatory variables while accounting for the effects of one or more other explanatory variables.  For example, you might want to account for the effect of environmental variables while accounting for spatial coordinates.  Borcard et al. (2018, ch. 6.3.2.5) illustrates this approach.

RDA in R (vegan::rda())

In R, RDA can be conducted using function rda() in vegan.  Its usage is:

rda(formula,
data,
scale = FALSE,
na.action = na.fail,
subset = NULL,
...
)

The arguments are:

• formula – model formula, with the community data matrix on the left hand side and the explanatory (constraining) variables on the right hand side.  Factors are converted to a series of linear contrasts (one per df).  Additional explanatory variables can be included that ‘condition’ the response – their effect is removed before testing the constraining variables.  See the help files for details of how to include conditioning variables.  The explanatory and response variables can also be specified as separate arguments, but the formula approach is preferred as it allows you to easily control which variables are being used.
• data – the data frame containing the explanatory variables.
• scale – whether to scale species (i.e., columns in the response matrix) to a variance of 1 (like a correlation). Default is to not do so (FALSE).  While this often makes sense with other types of responses, Borcard et al. (2018) state that this should not be used with community data.
• na.action – what to do if there are missing values in the explanatory variables. Several options:
  • na.fail – stop analysis. The default.
  • na.omit – drop rows with missing values and then proceed with analysis.
  • na.exclude – keep all observations but return NA for values that cannot be calculated.
• subset – whether to subset the data, and how to do so.

The resulting object contains lots of different elements.  An object of class ‘rda’ can be subject to pre-defined methods for plot(), print(), and summary().  Other types of available functions include:

• Extracting components: eigenvals(), fitted(),  model.frame(), model.matrix(), nobs(), residuals(), RsquareAdj(), scores()
• Model diagnostics: hatvalues(), rstandard(), rstudent()
• Model fitting: anova(), permutest()
• Model selection: add1(), drop1(), ordistep(), ordiR2step(), vif.cca()
• Predictions: predict(), calibrate()
• Graphing: plot(), points(), text()

Extracting Coordinates (vegan::scores())

We can also extract the coordinates to create our own graphics, etc.  The scores() function that we’ve used before also applies to objects of class ‘rda’.  Its usage for this class of object is:

scores(
x,
choices = c(1,2),
display = c("sp","wa","bp","cn"),
scaling = "species",
const,
correlation = FALSE,
tidy = FALSE,
...
)

The key arguments are:

• x – the object produced by rda()
• choices – which axes to return.  Default is to return just the first two axes.
• display – which coordinates to return.  Options:
  • species – coordinates of each response variable
  • wa or sites – coordinates of each sample unit.  Borcard et al. (2018) describe these as the weighted sums of species (i.e., response variables) scores and “coordinates of the sites [i.e., sample units] as expressed in the space of the response variables” (p. 211).
  • lc – linear constraints – the fitted value for each sample unit based on its explanatory values.  Borcard et al. (2018) describe these linear combinations of the constraining variables and as the “coordinates of the sites [i.e., sample units] in the space of the explanatory variables” (p. 211).  Note that, when fitted to a factor, this predicts the same value for all sample units that belong to the same level.
  • bp – coordinates that can be used to locate the head of biplot arrows (vectors) associated with explanatory variables
  • cn – for an explanatory variable that is a factor, the location of each level of the factor on each axis
  • reg – regression coefficients to find LC scores from constraints
  • all – return all sets of coordinates.  Appears to be the only option when tidy = TRUE.
• tidy – whether to compile all coordinates in a data frame.  Default is to not do so (FALSE).  My testing indicates that the display argument is ignored and treated as display = "all" if this is TRUE.

The tidy = TRUE argument provides various sets of coordinates for each axis specified.  The resulting object also includes a score column that identifies the set of coordinates and a label column that uniquely identifies each value (sample unit, species, etc.).  The score column contains more interpretable names than are used in the display argument above:

• species – same as display = "species".  One row per response variable.
• sites – same as display = "sites" or display = "wa".  One row per sample unit.
• constraints – same as display = "lc".  One row per sample unit.
• biplot – same as display = "bp", but only for continuously distributed explanatory variables.  Coordinates that can be used to add biplot arrows to a graphic. One row per df.  Only present for constrained axes (0 for unconstrained axes).
• regression – same as display = "reg".  One row per df of the explanatory variables.  Only present for constrained axes (0 for unconstrained axes).
• centroids – same as display = "cn".  One row per level of explanatory variables that are factors.  Only present for constrained axes (0 for unconstrained axes).
• factorbiplot – same as display = "bp", but only for categorical explanatory variables.  Coordinates that can be used to add biplot arrows to a graphic.  One row per df.  Only present for constrained axes (0 for unconstrained axes).

Darlingtonia Example (RDA)

We used a dataset about Darlingtonia (pitcher plants) to illustrate PCA by analyzing plant traits.  The dataset also includes information about which of four sites each plant was from.  Let’s use RDA to illustrate how plants from these sites differ in plant traits.

Import Data

We begin by importing the data.  For clarity later, we’ll assign a rowname to each line of data.

library(tidyverse, vegan)

darl <- read.csv("data/Darlingtonia_GE_Table12.1.csv")

rownames(darl) <- with(darl, paste(site, plant, sep = "_"))

Test for Differences Among Sites

Now we create an object that contains the responses.  Since the plant traits are measured in very different units, we will scale (i.e., normalize) them as we did during the PCA.

darl.resp <- darl |>
select(-c(site, plant)) |>
scale()

Do the set of plant traits differ among plants from different sites?

adonis2(
darl.resp ~ site,
data = darl,
method = "euclidean")

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = darl.resp ~ site, data = darl, method = "euclidean")
         Df SumOfSqs      R2      F Pr(>F)    
site      3   181.87 0.21148 7.4202  0.001 ***
Residual 83   678.13 0.78852                  
Total    86   860.00 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Yes, there are significant differences between sites.

Conduct RDA

Now, let’s conduct a RDA to help us visualize those differences.

darl_rda <- rda(
darl.resp ~ site,
data = darl,
scale = TRUE)

View the results:

print(darl_rda)

Call: rda(formula = darl.resp ~ site, data = darl, scale = TRUE)

              Inertia Proportion Rank
Total         10.0000     1.0000     
Constrained    2.1148     0.2115    3
Unconstrained  7.8852     0.7885   10
Inertia is correlations 

Eigenvalues for constrained axes:
  RDA1   RDA2   RDA3 
1.5075 0.5074 0.1000 

Eigenvalues for unconstrained axes:
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10 
3.441 1.606 0.910 0.663 0.359 0.328 0.222 0.199 0.112 0.046 

Inertia is another name for variance – these are the same sums of squares that are reported in a PERMANOVA of these data (note correspondence with above analysis!).  In this case, 21.15% of the variation in plant traits is accounted for by site identity.  The other 78.85% of variation in plant traits is unconstrained. 

Since our explanatory variable (site) had 3 df, three axes are constrained to maximize the correlation between the response matrix and the explanatory variable.  The variation explained by the explanatory variable was subject to PCA, resulting in axes (RDA1 through RDA3) that are in declining order of importance.  These axes are, “in successive order, a series of linear combination of the explanatory variables that best explain the variation of the response matrix” (Borcard et al. 2018, p.155).  In other words, each axis represents a multiple regression model based on all explanatory variables, that yield fitted values for the response.  There is an eigenvalue associated with each axis.

The variation that could not be explained by site is unconstrained.  This variation was also subject to PCA, resulting in axes (PC1 through PC10) that are in declining order of importance.  There is an eigenvalue associated with each axis.

Recall that the eigenvalues measure the variance accounted for by each axis – the proportion of variance explained by each eigenvalue is its value divided by the sum of all of the eigenvalues.  The full set of eigenvalues can be returned using eigenvals().  Applying summary() to this object also returns the proportion of variance explained by each term:

eigenvals(darl_rda) |>
summary()

Importance of components:
                        RDA1    RDA2     RDA3    PC1    PC2     PC3     PC4     PC5     PC6     PC7     PC8     PC9     PC10
Eigenvalue            1.5075 0.50737 0.099965 3.4409 1.6058 0.90969 0.66269 0.35878 0.32794 0.22207 0.19885 0.11210 0.046337
Proportion Explained  0.1507 0.05074 0.009997 0.3441 0.1606 0.09097 0.06627 0.03588 0.03279 0.02221 0.01989 0.01121 0.004634
Cumulative Proportion 0.1507 0.20148 0.211480 0.5556 0.7162 0.80712 0.87339 0.90927 0.94206 0.96427 0.98416 0.99537 1.000000 

The cumulative proportion explained sums to 1 after all axes have been considered.  The RDA* axes were constrained and the PC* axes were unconstrained.  Each set of axes is arranged as a declining function, but note the jump – in magnitude of eigenvalue and thus proportion explained – between RDA3 and PC1.

Create Visualization from RDA

To graph these data, we’ll use the gg_ordiplot() function from the ggordiplots package.  See the ‘Visualizing and Interpreting Ordinations’ chapter for more information about this package and other ways to customize ordination graphics.

library(ggordiplots)

gg_ordiplot(
darl_rda,
groups = darl$site,
hull = TRUE,
ellipse = FALSE)

Plants from the HD site appear to be quite different than those from other sites, and there is very little difference between plants from the DG and TJH sites.

In this case, there are 3 constrained axes associated with differences among sites.  The eigenvalues of these axes (shown above) sum to 0.2115.  The eigenvalues of the two axes that we’ve graphed here sum to 0.2015 (there are slight differences due to rounding in the values reported in the graphic).  These two axes therefore account for 0.2015 / 0.2115 = 0.953 or 95.3% of the variation that can be explained by site.

Extract RDA Coordinates

Let’s extract the coordinates for the first six axes (i.e., RDA1 through RDA3, plus PC1 through PC3):

darl.points <- scores(
darl_rda,
choices = c(1:6),
tidy = TRUE)

This object contains multiple types of coordinates:

darl.points |>
count(score) 

         score  n
1    centroids  4
2  constraints 87
3 factorbiplot  3
4   regression  3
5        sites 87
6      species 10

Whenever we’re going to use these points, we’ll have to filter to the desired set of coordinates.

We can use these points in multiple ways.  For example, we can filter to the site scores and then merge them with the explanatory variables:

darl.points.sites <- darl.points |>
filter(score == "sites") |>
merge(y = darl[ , c("site", "plant")],
  by = "row.names")

We now have everything we need to create our own graphics.

The ‘constraints’ set of coordinates are the fitted values for each sample unit based on its explanatory values.  Since the last step in the algorithm was a PCA, what is the expectation regarding correlations between these axes?

darl.points |>
filter(score == "constraints") |>
select( ! any_of(c("score", "label"))) |>
cor() |>
round(digits = 3)

     RDA1 RDA2 RDA3 PC1 PC2 PC3
RDA1    1    0    0   0   0   0
RDA2    0    1    0   0   0   0
RDA3    0    0    1   0   0   0
PC1     0    0    0   1   0   0
PC2     0    0    0   0   1   0
PC3     0    0    0   0   0   1

The constrained axes and unconstrained axes are uncorrelated with one another, as we saw previously with PCA and PCoA.

The ‘regression’ coordinates:

darl.points |>
filter(score == "regression")

              RDA1       RDA2      RDA3 PC1 PC2 PC3      score   label
siteHD  -1.0733949 -0.3199807 0.1459375   0   0   0 regression  siteHD
siteLEH -0.1468306 -0.8103160 0.8643588   0   0   0 regression siteLEH
siteTJH -0.4601139  0.3262829 1.0521986   0   0   0 regression siteTJH

Note that these coordinates are only present (non-zero) for the constrained axes and are absent from the unconstrained axes.  With a categorical explanatory variable as in this case, the levels are represented by dummy variables representing a series of contrasts – in this case, each contrast distinguishes one site from the first site in the dataset, DG.  These coefficients are interpreted like the loadings from a PCA – they are the eigenvectors associated with each eigenvalue.  The loadings can be matrix multiplied by either the original response matrix or the matrix of fitted values to locate each sample in the two ordination spaces mentioned above:

• Site scores – result of multiplying loadings by the original variables. Stored as darl.points |> filter(score == "sites").
• Site constraints – result of multiplying loadings by the fitted values (constrained axes) and residuals (unconstrained axes). Stored as darl.points |> filter(score == "constraints").

distance-based ReDundancy Analysis (dbRDA)

Approach

ReDundancy Analysis (RDA) is a constrained ordination that uses Euclidean distances.  Distance-based ReDundancy Analysis (dbRDA or db-RDA; Legendre & Anderson 1999) is a more general way to do the same type of constrained ordination; it can be applied to any distance or dissimilarity matrix.

The connection between RDA and dbRDA is simply the step of expressing a non-Euclidean dissimilarity matrix in a Euclidean space.  Sound familiar?  That’s what a PCoA does.

As a result, the first step of a dbRDA is to apply a PCoA to the distance matrix and to keep all principal coordinates.  As discussed in that chapter, this can be done regardless of the distance measure used.  By keeping all of the principal coordinates we have rotated the distance matrix and expressed it in Euclidean distances without reducing its dimensionality.

The second step of a dbRDA is to conduct a RDA on the principal coordinates.  In other words, the principal coordinates are subject to a multivariate multiple linear regression and then the fitted values and residuals are subject to separate PCAs.  This procedure maximizes the fit between the distance matrix containing the PCoA coordinates and the explanatory variables or constraints.

Notes:

1. If applied to a Euclidean distance matrix, the PCoA has no effect and the dbRDA is identical to RDA.
2. If no explanatory variables are provided, a dbRDA is identical to a PCoA – the first step results in a regular PCoA and the second step is meaningless.
3. Canonical Analysis of Principal coordinates (CAP) is another technique that has been proposed (McArdle & Anderson 2001; Anderson & Willis 2003). According to Borcard et al. (2018), it analyzes the dissimilarity matrix directly and does not require the initial PCoA. However, the help file for vegan::capscale() indicates that it is similar to dbRDA, and it appears that it is not maintained as a separate function.

db-RDA in R (vegan::dbrda() and vegan::capscale())

The help file describes these functions as ‘alternative implementations of dbRDA’.  These functions differ somewhat in their approach – including how they deal with and report on negative eigenvalues – but generally should give very similar results; see the help file for details.  Both are wrapper functions that call on other functions for various steps in the analysis.  I focus here on dbrda().

The usage of dbrda() is:

dbrda(formula,
data,
distance = "euclidean",
sqrt.dist = FALSE,
add = FALSE,
dfun = vegdist,
metaMDSdist = FALSE,
na.action = na.fail,
subset = NULL,
...
)

The key arguments are:

• formula – model formula, with a community data frame or distance matrix on the left hand side and the explanatory (constraining) variables on the right hand side. Additional explanatory variables can be included that condition the response – their effect is removed before testing the constraining variables.  The explanatory and response variables can also be specified as separate arguments, but the formula approach is preferred as it allows you to easily control which variables are being used.
• data – the data frame containing the explanatory variables.
• distance – type of distance matrix to be calculated if the response is a data frame. Default is to use Euclidean distance.
• sqrt.dist – whether to transform the distance matrix by taking the square root of all dissimilarities. This can be helpful to ‘euclidify’ some distance measures.  Default is to not do so (FALSE).
• add – when conducting the PCoA, whether to add a constant to the dissimiliarities so that none of them are non-negative. Same options as for wcmdscale():  FALSE or TRUE or "lingoes" or "cailliez", with the default to not add a constant (FALSE).
• dfun – name of function used to calculate distance matrix. Default is vegdist (i.e., vegdist()).
• metaMDSdist – whether to use the same initial steps as in metaMDS() (auto-transformation, stepacross for extended dissimilarities). Default is to not do so (FALSE).
• na.action – what to do if there are missing values in the constraints (explanatory variables). Same options as for wcmdscale().  Default is to stop (na.fail); other options are to drop sample units with missing values (na.omit or na.exclude).
• subset – whether to subset the data, and how to do so.

The resulting object contains lots of different elements.  Objects of class ‘dbrda’ or ‘capscale’ can be subject to pre-defined methods for plot(), print(), and summary().  There are not as many other functions available for objects of these classes as for objects of class ‘rda’.

By default, only the first 6 axes are reported – since the objective is to view patterns related to one or more explanatory variables, the user is unlikely to need or use other axes.

The scores() function will extract the sets of coordinates produced by dbrda().

Oak Example (dbRDA)

Let’s use dbRDA to visualize the relationship between oak plant community composition and geographic variables.

source("scripts/load.oak.data.R")

geog.stand <- Oak %>%
select(c(LatAppx, LongAppx, Elev.m)) %>%
decostand("range")

Do these explanatory variables relate to the composition of the stands?  We want to consider the set of explanatory variables en masse so we set by = NULL:

adonis2(
Oak1 ~ Elev.m + LatAppx + LongAppx,
data = geog.stand,
by = NULL)

Some variation is associated with these variables, but the effect is marginally significant.

Although the patterns are marginally significant, for illustration purposes we’ll proceed with the dbRDA.

Oak1_dbrda <- dbrda(Oak1 ~ Elev.m + LatAppx + LongAppx,
data = geog.stand,
distance = "bray")

print(Oak1_dbrda)

Call: dbrda(formula = Oak1 ~ Elev.m + LatAppx + LongAppx, data =
geog.stand, distance = "bray")

               Inertia Proportion Rank RealDims
Total         11.59490    1.00000              
Constrained    0.97790    0.08434    3        3
Unconstrained 10.61700    0.91566   43       33
Inertia is squared Bray distance 

Eigenvalues for constrained axes:
dbRDA1 dbRDA2 dbRDA3 
0.4568 0.3043 0.2168 

Eigenvalues for unconstrained axes:
  MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8 
1.7901 1.1519 1.0933 0.7744 0.7069 0.6514 0.5212 0.4344 
(Showing 8 of 43 unconstrained eigenvalues)

In this case, 8.43% of the variance is accounted for by the constraining variables, the set of geographic variables.  The other 91.57% of the variance is unconstrained – it cannot be explained by the spatial location of the sample units.

The axes are interpreted as for a RDA:

• The constrained axes (dbRDA1 through dbRDA3) are in declining importance
• The unconstrained axes (MDS1 and on) are in declining importance.  Since the dbRDA was applied to a distance matrix that does not satisfy the triangle inequality, some imaginary axes were created during the PCoA, and these imaginary axes have negative eigenvalues.  The ‘RealDims’ element in the results indicates how many of the axes have positive eigenvalues; the remaining 43 – 33 = 10 axes had negative eigenvalues and are not shown.
• The magnitude of the eigenvalue can jump – sometimes dramatically! – between the last constrained axis and the first unconstrained axis

To view the resulting ordination:

plot(Oak1_dbrda)

Notes:

• The first two axes are shown by default. If two or more axes are constrained, these are what will be shown.  However, if the dbRDA was based on a single explanatory variable with one df (i.e., only one constrained axis), the constrained axis and the first unconstrained axis will be shown.
• Since the explanatory variables in this analysis are continuously distributed, they are shown as vectors on the ordination.  Each vector shows the direction in which that explanatory variable is most strongly associated with the axes.
• The location of each sample unit is shown by its name.

The scores() function can extract the scores as with rda().   for use in ggplot2.  However, the resulting object is a bit more complex than we’ve seen with other types of ordinations.

Oak1.dbrda.points <- scores(
Oak1_dbrda,
tidy = TRUE)

Oak1.dbrda.points |>
count(score)

        score  n
1      biplot  3
2 constraints 47
3  regression  3
4       sites 47

This object does not contain all of the same elements as we saw above for rda().  Some of the missing elements reflect the fact that this analysis had continuous rather than categorical explanatory variables.  The ‘species’ scores are missing because we converted the data matrix to a distance matrix before analysis.

Interpreting RDA and dbRDA

RDA uses Euclidean distances, so it is most appropriate when this distance measure is appropriate for the response matrix.  If no constraining variables are included, the result is a regular PCA.

Although RDA and dbRDA begin with a multivariate linear regression model, Borcard et al. (2018, ch. 6.3.2.10) note that it is possible to also test second-degree (i.e., quadratic) variables.  This can be useful if responses are unimodal, for example.

A major difficulty with all constrained ordinations is choosing which explanatory variables to constrain the ordination with.  Model-testing and model-building functions can be applied to identify the optimal set of explanatory variables to include.  See Borcard et al. (2018, ch. 6.3.2.6) for details.  As more constraining variables are included, the model ends up increasingly behaving like an unconstrained ordination because there is more flexibility within the constraints.

If an explanatory variable with one df is tested (e.g., a continuous variable or a factor with two levels), there will only be one constrained axis.  By default, the resulting graphics will therefore include the constrained axis and the first unconstrained axis.

Finally, note that the set of explanatory variables is considered en masse; the loadings demonstrate how the score on each constrained axis is a linear combination of those explanatory variables.  These techniques are therefore not designed to evaluate which explanatory variables would be important if considered independently.  Instead, I think they are more appropriately used to visualize the variables identified as important through separate statistical analyses.

Conclusions

Several of these techniques are closely related.  For example, a RDA without constraining explanatory variables is a PCA, and a dbRDA without constraining explanatory variables is a PCoA.

RDA and dbRDA blur the distinction between statistical testing and visualization.  They can be used to create visuals that support other statistical tests, but they can also be used to do statistical testing, including model selection to identify parsimonious sets of explanatory variables.

Like Mantel tests and RRPP, RDA and dbRDA have the ability to control for one set of explanatory variables while examining the relationship of another set to the response variables.

RDA does not have CCA’s issues related to the use of the chi-square distance measure, though it does assume Euclidean distances and therefore is appropriate for variables where that distance measure is appropriate.

dbRDA can be applied to a distance matrix obtained using any distance measure.

References

Anderson, M.J., and T.J. Willis. 2003. Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology. Ecology 84:511-525.

Borcard, D., F. Gillet, and P. Legendre. 2018. Numerical Ecology with R. 2nd edition. Springer, New York, NY.

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

Legendre, P., and L. Legendre. 2012. Numerical Ecology. Third English edition. Elsevier, Amsterdam, The Netherlands.

Manly, B.F.J., and J.A. Navarro Alberto. 2017. Multivariate Statistical Methods: A Primer. Fourth edition. CRC Press, Boca Raton, FL.

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.