---
title: "capscale vs dbrda: same axes, different test"
description: "vegan has two functions for distance-based constrained ordination. On Bray-Curtis they return the same axes but a different pseudo-F: the negative eigenvalues."
date: "2026-05-03 13:00"
categories: [R, ordination, vegan, multivariate, ecology tutorial]
image: thumbnail.png
image-alt: "A bar chart of residual eigenvalues with a long positive tail and a shorter negative tail at the end."
date-modified: "2026-10-07"
---
*Corrected 7 October 2026: the post presented the lower adjusted R-squared of capscale, and its proportions taken against the real total, as general behaviour; they belong to newer vegan releases, and the post now says that in older ones, including the release in its References, the two functions give almost the same adjusted R-squared.*
*Corrected 27 September 2026: the recommendation said the two functions agree whenever the dissimilarity is metric, and called Jaccard semimetric; the condition is that the dissimilarity is Euclidean, and Jaccard is metric but not Euclidean, so it still belongs with Bray-Curtis on the dbrda side.*
The [post on constrained ordination](../constrained-ordination-dbrda/) used `dbrda` to put environmental constraints on a Bray-Curtis dissimilarity. vegan has a second function for the same job, `capscale`, and most of the time the two are interchangeable. On a semimetric dissimilarity like Bray-Curtis they are not quite. They return the same constrained ordination, but a different pseudo-F, and the split happens at one specific place: the negative eigenvalues.
The two functions take different routes. `capscale` (Legendre and Anderson 1999) runs principal coordinates analysis on the dissimilarities, then ordinary RDA on the resulting coordinates. `dbrda` (McArdle and Anderson 2001) works directly on the double-centred distance matrix, without the intermediate coordinates. For a Euclidean distance the routes are equivalent. For a non-Euclidean one they diverge, because the dissimilarities cannot be embedded in real Euclidean space without an imaginary part, and that part shows up as negative eigenvalues.
## A community and its dissimilarity
Two predictors drive composition; counts are Poisson around Gaussian niches. Bray-Curtis on this kind of data is semimetric, so it carries negative eigenvalues.
```{r}
#| label: setup
#| code-fold: true
#| code-summary: "Setup: packages, helper functions, colours and plot theme (run this first)"
#| results: hide
#| message: false
#| warning: false
library(ggplot2)
pal <- list(paper="#f5f4ee", ink="#16241d", body="#2c3a31", forest="#275139",
line="#dad9ca", sage="#93a87f", gold="#cda23f", terra="#b5534e",
faint="#5d6b61", low="#c9b458", high="#1d5b4e")
th <- function(base=12) theme_minimal(base_size=base) +
theme(panel.grid.minor=element_blank(),
panel.grid.major=element_line(color=pal$line, linewidth=.3),
plot.background=element_rect(fill=pal$paper, color=NA),
panel.background=element_rect(fill=pal$paper, color=NA),
axis.title=element_text(color=pal$body), axis.text=element_text(color=pal$faint),
legend.title=element_text(color=pal$body), legend.text=element_text(color=pal$faint),
plot.title=element_text(color=pal$ink, face="bold"), legend.position="top")
```
```{r}
#| label: data
#| warning: false
#| message: false
library(vegan)
set.seed(37)
S <- 24; n <- 60
x1 <- runif(n, -1, 1); x2 <- runif(n, -1, 1)
opt1 <- seq(-1.2, 1.2, length.out = S)
opt2 <- sample(seq(-1.2, 1.2, length.out = S))
lambda <- sapply(1:S, function(j)
18 * exp(-0.5*((x1 - opt1[j])/0.6)^2 - 0.5*((x2 - opt2[j])/0.7)^2))
comm <- matrix(rpois(n * S, lambda), n, S)
comm[rowSums(comm) == 0, 1] <- 1
env <- data.frame(x1, x2)
D <- vegdist(comm, "bray")
```
How non-Euclidean is this dissimilarity? Principal coordinates returns one eigenvalue per axis, and the negative ones are the imaginary part.
```{r}
#| label: eigen
#| warning: false
#| message: false
ev <- eigen(-0.5 * (diag(n) - 1/n) %*% as.matrix(D^2) %*% (diag(n) - 1/n),
symmetric = TRUE, only.values = TRUE)$values
c(positive = sum(ev > 1e-8), negative = sum(ev < -1e-8),
pct_negative = round(sum(abs(ev[ev < 0])) / sum(ev[ev > 0]) * 100, 1))
```
Many of the axes carry negative eigenvalues, and together they come to `r sprintf("%.1f", sum(abs(ev[ev < 0])) / sum(ev[ev > 0]) * 100)` per cent of the positive inertia. Small, but it is exactly what the two functions handle differently.
## The same axes
```{r}
#| label: fit
#| warning: false
#| message: false
cap <- capscale(D ~ x1 + x2, data = env, comm = comm) # comm only feeds species scores
db <- dbrda(D ~ x1 + x2, data = env)
cap
```
```{r}
#| label: print-db
#| warning: false
#| message: false
db
```
Read the two partition tables together. Total inertia is the same. The constrained inertia is the same to three figures, a little over seventy per cent of the `Total` line in both; the proportion column of the `capscale` table above is taken against its real total instead, as explained below. The constrained eigenvalues match to two decimal places: `CAP1` `r sprintf("%.3f", unname(cap$CCA$eig[1]))` against `dbRDA1` `r sprintf("%.3f", unname(db$CCA$eig[1]))`, `CAP2` `r sprintf("%.3f", unname(cap$CCA$eig[2]))` against `dbRDA2` `r sprintf("%.3f", unname(db$CCA$eig[2]))`. As a constrained ordination the two results are the same picture.
The difference is in the bookkeeping below the constrained part. `capscale` reports an `Imaginary` line: it has set the negative-eigenvalue inertia aside, outside both the constrained and the real unconstrained totals. `dbrda` has no imaginary line; its unconstrained total already nets the negative eigenvalues against the positive ones, which is why its unconstrained inertia is smaller.
```{r}
#| label: compare
#| warning: false
#| message: false
rbind(
capscale = c(adjR2 = RsquareAdj(cap)$adj.r.squared,
resid = cap$CA$tot.chi),
dbrda = c(adjR2 = RsquareAdj(db)$adj.r.squared,
resid = db$CA$tot.chi))
```
The adjusted R-squared values are not the same: 0.65 for `capscale` against 0.72 for `dbrda`. The constrained inertia behind them agrees to three figures; what differs is the total each function divides it by. `capscale` divides by the real total, which leaves the imaginary inertia outside and so is the larger denominator, while `dbrda` divides by a total that has already netted the negative eigenvalues away. The residual inertia splits for the same reason. The adjusted R-squared of `capscale` depends on the vegan release: newer releases, which produced the output above, take its printed proportions and its adjusted R-squared from the real components only, while older ones, including the release listed in the References, divide by the same total as `dbrda`, so that the two adjusted values come out almost equal.
## A different denominator, a different F
The permutation test divides constrained variation by residual variation. Same numerator, different residual, different pseudo-F.
```{r}
#| label: anova
#| warning: false
#| message: false
set.seed(1); ac <- anova(cap, permutations = 999)
set.seed(1); ad <- anova(db, permutations = 999)
c(capscale_F = round(ac$F[1], 1), dbrda_F = round(ad$F[1], 1),
capscale_p = ac$`Pr(>F)`[1], dbrda_p = ad$`Pr(>F)`[1])
```
The dbRDA pseudo-F comes out larger, near seventy-seven against fifty-six. The mechanism is direct: `dbrda` subtracts the negative-eigenvalue inertia from the residual, so its denominator is smaller, so its F is larger. Here both p values land at the permutation floor and the conclusion is the same, but in a borderline case the choice of function could move the p value across a threshold. When the test is the point, that matters.
This is the residual structure that drives the gap:
```{r}
#| label: fig-eig
#| fig-cap: "Residual (unconstrained) eigenvalues. The positive tail is shared. The negative tail at the end is the non-Euclidean part: capscale sets it aside as imaginary, dbrda keeps it in the residual and so divides by a smaller number."
#| fig-alt: "Bar chart of residual eigenvalues: a long row of positive green bars followed by a long run of short red negative bars that deepen towards the right end."
#| fig-width: 6.6
#| fig-height: 3.4
#| warning: false
#| message: false
ca <- db$CA$eig
dfe <- data.frame(idx = seq_along(ca), eig = as.numeric(ca),
sign = ifelse(ca < 0, "negative", "positive"))
ggplot(dfe, aes(idx, eig, fill = sign)) +
geom_col(width = .85) +
geom_hline(yintercept = 0, color = pal$body, linewidth = .3) +
scale_fill_manual(values = c(positive = pal$forest, negative = pal$terra), name = NULL) +
labs(x = "Residual (unconstrained) axis", y = "Eigenvalue") + th()
```
And the constrained ordination itself, which both functions agree on:
```{r}
#| label: fig-ord
#| fig-cap: "The constrained ordination (from dbrda; capscale gives the same axes). Sites are coloured by the first predictor, which lines up along the first axis."
#| fig-alt: "Distance-based RDA ordination scatter with sites coloured along a gradient by the first predictor and two constraint arrows, the first predictor's running along the first axis and the second predictor's along the second."
#| fig-width: 6.2
#| fig-height: 4.6
#| warning: false
#| message: false
si <- as.data.frame(scores(db, display = "sites", choices = 1:2)); si$x1 <- x1
bp <- as.data.frame(scores(db, display = "bp", choices = 1:2))
mul <- 0.9 * max(abs(unlist(si[, 1:2]))) / max(sqrt(rowSums(bp^2)))
bp <- bp * mul
eig <- round(100 * db$CCA$eig / db$tot.chi, 1)
ggplot(si, aes(dbRDA1, dbRDA2, color = x1)) +
geom_point(size = 2.4) +
scale_color_gradient(low = pal$low, high = pal$high, name = "first predictor") +
geom_segment(data = bp, aes(0, 0, xend = dbRDA1, yend = dbRDA2), inherit.aes = FALSE,
arrow = arrow(length = unit(.18, "cm")), color = pal$terra, linewidth = .6) +
geom_text(data = bp, aes(dbRDA1, dbRDA2,
label = c(x1 = "first predictor", x2 = "second predictor")[rownames(bp)]),
inherit.aes = FALSE, color = pal$terra, vjust = -.5, size = 3.5) +
labs(x = paste0("dbRDA1 (", eig[1], "%)"), y = paste0("dbRDA2 (", eig[2], "%)")) + th()
```
## With a Euclidean distance they agree exactly
Hellinger distance is Euclidean, so it produces no negative eigenvalues and no imaginary part. With nothing to treat differently, the two functions return the same constrained inertia to numerical precision.
```{r}
#| label: euclid
#| warning: false
#| message: false
DH <- dist(decostand(comm, "hellinger"))
capH <- capscale(DH ~ x1 + x2, data = env)
dbH <- dbrda(DH ~ x1 + x2, data = env)
c(capscale = capH$CCA$tot.chi, dbrda = dbH$CCA$tot.chi,
difference = abs(capH$CCA$tot.chi - dbH$CCA$tot.chi))
```
The difference is at the rounding limit. The whole divergence on Bray-Curtis was the negative eigenvalues, and a Euclidean distance has none.
## One practical difference, and what to use
`capscale` can give you species scores, which is convenient for a biplot, but not from a dissimilarity object alone. When the left-hand side of the formula is a `dist`, the species table has to be handed over separately as `comm`, which is why the fit above was called with `comm = comm`. `dbrda`, the function of the [dbRDA post](../constrained-ordination-dbrda/), has no such argument at all: the analysis runs on the double-centred distance matrix, not on a species table, so species coordinates only appear if you add them afterwards with `sppscores`.
```{r}
#| label: spp
#| warning: false
#| message: false
cap_bare <- capscale(D ~ x1 + x2, data = env) # same fit, no species table handed over
grab <- function(m) tryCatch(scores(m, display = "species", choices = 1:2),
error = function(e) NULL)
count <- function(s) if (is.null(s) || all(is.na(s))) "none" else nrow(s)
list(capscale_without_comm = count(grab(cap_bare)),
capscale_with_comm = count(grab(cap)),
dbrda_species_scores = count(grab(db)))
```
The recommendation follows the mechanism. When the dissimilarity is Euclidean (Hellinger, chord, Euclidean), the two functions agree, so use either; `capscale` with `comm` supplied is convenient if you want species scores for the biplot. When the dissimilarity is non-Euclidean (Bray-Curtis, Jaccard) and the test is what you report, prefer `dbrda`: it partitions the variation directly from the distance matrix, including the non-Euclidean part, which is the partition McArdle and Anderson (2001) argued for. If you would rather stay with `capscale` on a non-Euclidean dissimilarity, the negative eigenvalues can be removed first, by taking square roots of the dissimilarities or by adding a constant (`add = TRUE`); the vegan help documents both.
A note on the name. `capscale` takes its label from CAP, the canonical analysis of principal coordinates of Anderson and Willis (2003). That method is a discriminant or correlation analysis on principal coordinates, aimed at an a priori hypothesis, and it is not identical to what `capscale` computes, which is RDA on principal coordinates. The family is the same; the exact engine is not.
## References
Legendre P, Anderson MJ 1999 Ecological Monographs 69(1):1-24 (10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2)
McArdle BH, Anderson MJ 2001 Ecology 82(1):290-297 (10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2)
Anderson MJ, Willis TJ 2003 Ecology 84(2):511-525 (10.1890/0012-9658(2003)084[0511:CAOPCA]2.0.CO;2)
Oksanen J, Simpson GL, Blanchet FG, et al. 2022 vegan: Community Ecology Package, R package version 2.6-4 (https://CRAN.R-project.org/package=vegan)
## Related tutorials
- [Constrained ordination with dbRDA](../constrained-ordination-dbrda/)
- [RDA vs CCA and gradient length](../rda-vs-cca-gradient-length/)
- [PCA on environmental data](../pca-environmental-data/)
- [Variance partitioning](../variance-partitioning-varpart/)
