Model-based unconstrained ordination

An NMDS axis can track sampling intensity rather than ecology when the mean-variance link is strong. A model-based ordination on residuals recovers it.

Non-metric multidimensional scaling is how most community ecologists look at their data: turn the site-by-species table into Bray-Curtis dissimilarities, squeeze them into two dimensions, and read the picture. It is an algorithm for preserving the rank order of dissimilarities, and it works hard at that. What it does not do is account for the statistical properties of counts, and that omission has a cost. Because Bray-Curtis is dominated by the most abundant species, and because abundance and variance travel together, the first axis often ends up tracking how many individuals were caught rather than which species were there.

A model-based ordination starts from the other end. Fit a generalised linear model to the table, with a site effect that absorbs total abundance and a species effect that absorbs how common each species is, then ordinate what is left over: the Dunn-Smyth residuals. Those residuals live on a common scale for every species, so a rare informative species counts as much as a common uninformative one. This post recovers a known gradient with both methods and shows the difference.

A gradient hidden behind sampling intensity

We simulate 60 sites and 25 species. A single latent gradient z drives composition, but the species that respond most strongly to it are the rarer ones; the abundant species that dominate a dissimilarity barely track it. On top of that, sites differ in sampling intensity, a multiplier on every species that has nothing to do with ecology. This is the situation where the two approaches part company.

set.seed(4223)
n <- 60; S <- 25
z <- runif(n, -2, 2)                          # true latent gradient (unknown to the methods)
b0     <- sort(rnorm(S, mean = log(c(rep(16,6), rep(5,7), rep(1.6,7), rep(0.7,5))), sd = 0.2),
               decreasing = TRUE)
lambda <- c(rnorm(6, 0, 0.05), rnorm(7, 0, 0.5), rnorm(7, 0, 1.1), rnorm(5, 0, 1.3))
theta_true <- 1.5
samp <- rnorm(n, 0, 0.6)                       # sampling intensity, uncorrelated with z

eta <- outer(rep(1,n), b0) + outer(z, lambda) + outer(samp, rep(1,S))
Y   <- matrix(rnbinom(n * S, mu = as.vector(exp(eta)), size = theta_true), n, S)
colnames(Y) <- paste0("sp", 1:S)
c(total = sum(Y), min_rowsum = min(rowSums(Y)), max_rowsum = max(rowSums(Y)))

     total min_rowsum max_rowsum 
     13371         24       1851 

Site totals range from 24 to 1851 individuals, a spread driven by sampling intensity rather than by the gradient we care about.

Distance-based ordination

We build Bray-Curtis by hand and run NMDS from a classical scaling start, iterating a monotone regression of the configuration distances against the dissimilarity ranks (the Kruskal idea, solved with a Guttman transform).

bray <- function(M) { n <- nrow(M); D <- matrix(0, n, n)
  for (i in 1:(n-1)) for (j in (i+1):n) {
    d <- sum(abs(M[i,] - M[j,])) / sum(M[i,] + M[j,]); D[i,j] <- D[j,i] <- d }; D }
Dbc <- bray(Y)

nmds <- function(D, k = 2, maxit = 300, seed = 1) {
  set.seed(seed); n <- nrow(D)
  X <- scale(cmdscale(D, k), scale = FALSE)
  low <- lower.tri(D); Dl <- D[low]; o <- order(Dl); Sprev <- Inf; Sval <- NA
  for (it in 1:maxit) {
    dx <- as.matrix(dist(X)); dl <- dx[low]
    iso <- isoreg(dl[o])$yf; dhat <- numeric(length(dl)); dhat[o] <- iso
    Sval <- sqrt(sum((dl - dhat)^2) / sum(dl^2))
    Dh <- matrix(0, n, n); Dh[low] <- dhat; Dh <- Dh + t(Dh)
    B <- matrix(0, n, n); nz <- dx > 1e-9; B[nz] <- -Dh[nz] / dx[nz]; diag(B) <- -rowSums(B)
    X <- scale((1/n) * B %*% X, scale = FALSE)
    if (abs(Sprev - Sval) < 1e-7) break; Sprev <- Sval
  }
  list(points = X, stress = Sval)
}
mds <- nmds(Dbc, k = 2, seed = 7)
db1 <- mds$points[,1]; db2 <- mds$points[,2]
mds_stress <- mds$stress
round(mds_stress, 3)

[1] 0.157

Model-based ordination

Now the model. We stack the table into long form, fit a negative binomial GLM with a site term and a species term, and take the Dunn-Smyth residuals: the fitted cumulative probability at each count, jittered across its probability mass and pushed through the normal quantile function. Reshaped back into a matrix, those residuals carry the site-by-species structure that the main effects could not explain, and a singular value decomposition reads off the axes.

long <- data.frame(count = as.vector(Y),
                   site  = factor(rep(1:n, times = S)),
                   sp    = factor(rep(1:S, each = n)))
fit    <- glm.nb(count ~ site + sp, data = long)   # site total + species mean, removed
mu_hat <- fitted(fit); th_hat <- fit$theta

ds_resid <- function(y, mu, theta, u) {
  a <- pnbinom(y - 1, mu = mu, size = theta)
  b <- pnbinom(y,     mu = mu, size = theta)
  qnorm(a + u * (b - a))
}
set.seed(99)
R  <- matrix(ds_resid(long$count, mu_hat, th_hat, runif(nrow(long))), n, S)
sv <- svd(scale(R, center = TRUE, scale = FALSE))
mb1 <- sv$u[,1] * sv$d[1]; mb2 <- sv$u[,2] * sv$d[2]
round(th_hat, 2)

[1] 0.86

Which axis found the gradient

The test is simple: correlate each ordination’s first axis with the true z, and also with log total abundance, the nuisance we hoped to be rid of.

cor_mb <- abs(cor(mb1, z)); cor_db <- abs(cor(db1, z))
lrs <- log(rowSums(Y))
cor_mb_tot <- abs(cor(mb1, lrs)); cor_db_tot <- abs(cor(db1, lrs))
c(mb_vs_z = round(cor_mb,3), nmds_vs_z = round(cor_db,3),
  mb_vs_total = round(cor_mb_tot,3), nmds_vs_total = round(cor_db_tot,3))

      mb_vs_z     nmds_vs_z   mb_vs_total nmds_vs_total 
        0.939         0.537         0.309         0.941 

Figure 1: Each point is a site. A tighter line means the ordination axis found the real gradient.

The model-based axis correlates 0.94 with the true gradient; the NMDS axis reaches 0.54. The NMDS axis also correlates 0.94 with log total abundance. The model-based axis, with the site term already absorbing total abundance, correlates 0.31 with it.

Figure 2: Sites coloured by the true gradient. A clean sweep means the gradient sits on an axis.

Note that the NMDS fit is not broken in its own terms: its stress is 0.157, in the band conventionally read as usable rather than good, where the broad layout can be believed but the fine detail cannot. It faithfully preserved the dissimilarity ranks. Low stress is a statement about the fit to the dissimilarities, not about ecology.

Honest limits

The model-based ordination is not free of assumptions, it simply makes them explicit. It assumes a negative binomial mean-variance relationship; if that is wrong the residuals are wrong and the axes drift, which is exactly what the residual diagnostics in the companion post are for. The latent axes are identified only up to rotation and sign, so the numbers on them mean nothing absolute and two runs can look mirrored. Whether to soak up total abundance with a site term is itself a modelling choice: if differences in total abundance are the ecology, removing them is the wrong move. And an ordination of any kind is a description of pattern, not a test of cause. An NMDS axis can be a sampling-intensity axis wearing an ecological label, and putting species on a common residual scale is what pulls the real gradient back into view.

References

Hui FKC, Taskinen S, Pledger S, Foster SD, Warton DI 2015. Methods in Ecology and Evolution 6:399-411 (10.1111/2041-210X.12236).

Warton DI, Blanchet FG, O’Hara RB, Ovaskainen O, Taskinen S, Walker SC, Hui FKC 2015. Trends in Ecology and Evolution 30:766-779 (10.1016/j.tree.2015.09.007).

Dunn PK, Smyth GK 1996. Journal of Computational and Graphical Statistics 5:236-244 (10.1080/10618600.1996.10474708).

Minchin PR 1987. Vegetatio 69:89-107 (10.1007/BF00038690).

Kruskal JB 1964. Psychometrika 29:1-27 (10.1007/BF02289565).

Legendre P, Legendre L 2012. Numerical Ecology, 3rd edition. Elsevier. ISBN 978-0-444-53868-0.

Related tutorials

• Model-based multivariate abundance
• Fourth-corner trait-environment analysis
• Checking a model-based multivariate model
• Common PERMANOVA mistakes