Model-based unconstrained ordination

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

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 

`geom_smooth()` using formula = 'y ~ x'

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 manages only 0.54. The reason is in the second pair of numbers: the NMDS axis correlates 0.94 with log total abundance, so its leading dimension is very nearly the sampling-intensity gradient in disguise. The model-based axis, with the site term already holding total abundance, keeps that correlation down at 0.31.

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, comfortably in the range usually called a good ordination. It faithfully preserved the dissimilarity ranks. The trouble is that those ranks encode sampling intensity, so a technically good ordination points at the wrong thing. 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. What this example does show is concrete: 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. Model-based approaches to unconstrained ordination. 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. So many variables: joint modeling in community ecology. Trends in Ecology and Evolution 30:766-779 (10.1016/j.tree.2015.09.007).

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

Minchin PR 1987. An evaluation of the relative robustness of techniques for ecological ordination. Vegetatio 69:89-107 (10.1007/BF00038690).

Kruskal JB 1964. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. 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