Analysis of relative abundances with zeros on environmental gradients: a multinomial regression model

Data aggregation

Due to the rarity of barnacles and Stomphia coccinea (three and one individuals respectively), these two taxa were excluded from the analysis. Points where these taxa were sampled were not redrawn, leaving one still with 91 points, three with 99, and the remainder with 100 points. The remaining taxa were combined into eight categories, consisting of organisms that were ecologically similar and/or could not be reliably distinguished: algae (red and green), Aurelia aurita polyps, Bugula spp., colonial ascidians (Botryllus schlosseri, Botrylloides leachii and Botrylloides violaceus), Diadumene cincta, solitary ascidians (Ciona intestinalis and Styela clava), sponges (Halichondria spp. and others), Mytilus edulis. We also included the “bare wall” category (for the absence of macroscopic organisms, although usually there was a biofilm of microscopic algae and bacteria, or a layer of detritus).

Statistical model

Let the counts in the ith observation (still image) be y_i = (y_i,1, y_i,2, …, y_i,9)^T, where y_i,j is the observed count of the jth taxon in the ith observation, and let $n_i={∑}_{j=1}^9{y}_{ij}$ be the total number of points counted for the ith observation (usually 100 in our data). Our model is ${\displaystyle {\mathbf{y}}_i∼\text{multinomial}\left(n_i,{\mathrm{ρ}}_i\right),}$ ${\displaystyle {\mathrm{ρ}}_i={\text{ilr}}^{−1}{\mathbf{x}}_i,}$ (1)${\displaystyle {\mathbf{x}}_i={\mathrm{β}}_0+{\mathrm{β}}_1{z}_i+{\mathrm{β}}_2{z}_i^2+{\mathrm{É›}}_i,}$ ${\displaystyle {\mathrm{É›}}_i∼N\left(\mathbf{0},\mathit{Î\pounds }\right).}$ In a non-Bayesian context, this model would be referred to as a multivariate generalized linear mixed model (Agresti, 2002, p.492), with a multinomial response distribution, an isometric logratio (ilr: Egozcue et al., 2003) link function, linear predictor x_i and random effects É›_i. The vector ρ_i is the expected relative abundance of each taxon. The multinomial observation model arises from the assumption that individual points within a still are drawn independently from a categorical distribution with probabilities ρ_i (Johnson, Kotz & Balakrishnan, 1997, p. 33). The ilr link function transforms the 8-simplex into an unconstrained 8-dimensional real space, with an ilr coordinate system described below. The linear predictor x_i is an 8-dimensional vector in ilr coordinates, and depends on β₀, β₁ and β₂, the unknown 8-dimensional intercept and linear and quadratic depth coefficient vectors respectively, and on z_i, the centred and scaled depth for the ith observation. The observation-specific intercepts É›_i are drawn from an 8-dimensional multivariate normal distribution in ilr coordinates, with mean vector 0 and covariance matrix Σ. These intercepts deal with extra-multinomial variation (overdispersion) arising from factors such as clustering due to the spatial extension of organisms and unmeasured covariates (McCullagh & Nelder, 1989, pp. 124–125, 174). In particular, in our data, variation in the distance of the ROV from the wall is likely to lead to varying amounts of overdispersion among stills. This treatment of overdispersion leads to a normal distribution of expected values on the simplex, in the sense of Pawlowsky-Glahn, Egozcue & Tolosana-Delgado (2015, p. 114). This distribution is much more flexible than, for example, a Dirichlet distribution, although there are many other reasonable choices.

It is important that observations y_i including zero counts are in the support of the multinomial distribution, and that fitting the model involves back-transforming the linear predictor (which is always in the domain of ilr^−1), not an ilr transformation of y_i. Thus, no special treatment of zeros (such as pseudocounts) is necessary. We fitted this model using Bayesian estimation via the NUTS algorithm (Hoffman & Gelman, 2014). NUTS is derived from Hamiltonian Monte Carlo, in which the problem of sampling from the posterior distribution of interest is formulated in terms of simulating the dynamics of a physical system with position, potential energy and momentum (Neal, 2011). This can explore the state space much more rapidly than random-walk methods such as the Metropolis–Hastings algorithm. NUTS improves on Hamiltonian Monte Carlo by requiring much less fine-tuning, and is implemented in the Stan programming language (Carpenter et al., 2017). We give more details in Section S3. We checked the model’s performance using a simulation study (Section S4). We used a Bayesian approach, despite the additional computation it involves, because it leads almost automatically to estimates of uncertainty in the compositional analyses described below. We compared the performance of this model against models with only a linear depth effect and with a cubic depth effect, using leave-one-out cross-validation to estimate the expected log predictive density for a new data set (Section S5).

The vector ρ_i consists of non-negative elements with a fixed sum of 1, and is therefore a composition. The sum constraint, and associated constraints on the covariance structure of compositions, make it difficult and inconvenient to specify sufficiently flexible parametric models for untransformed compositions (Aitchison, 1986, chapter 3). The most popular modern approach to analysis of compositional data is to transform an s-part composition into an unconstrained real space with s − 1 dimensions. We chose an isometric logratio transformation (Egozcue et al., 2003), which is an isomorphism (so that perturbation and powering in the simplex correspond to ordinary vector addition and scalar multiplication in the real space) and an isometry (so that distances under an appropriate norm in the simplex correspond to Euclidean distances in the real space).

The coordinates in an ilr coordinate system represent logcontrasts between groups of taxa (loglinear combinations of relative abundances whose coefficients sum to zero: Aitchison, 1986, p. 84). The ilr transformation is defined by a basis matrix, constructed from a set of s − 1 orthogonal logcontrasts. In principle, such logcontrasts can be very informative biologically. For example, in our study we would expect the logcontrast between algae and animals to decrease with depth, because algae were the only photosynthetic organisms included. We would expect the logcontrast between predatory and nonpredatory animals to increase with depth, because predatory animals do not rely on photosynthetic food, and we would expect the logcontrast between the two predators, A. aurita and D. cincta, to increase with depth because A. aurita polyps have a strong preference for dark locations (Ishii & Shioi, 2003).

In order to fit the model, we used the isometric logratio transformation with the default basis matrix in the R package ‘compositions’, version 1.40-1 (Van den Boogaart & Tolosana-Delgado, 2008). Our results do not depend on this choice of basis, but if it is important to be able to interpret logratio coordinates, an appropriate basis can be chosen by sequential binary partition (Egozcue & Pawlowsky-Glahn, 2005). We describe such a basis in the Section S6. Meaningful bases can also be constructed from hierarchical clustering of environmental preferences (Morton et al., 2017) or from a phylogeny (Silverman et al., 2017). Advantages and disadvantages of the ilr transformation, compared to other transformations, are discussed in Bacon-Shone (2011, section 1.5).

Because the isometric logratio transformation is an isomorphism between the simplex with Aitchison geometry (Pawlowsky-Glahn & Egozcue, 2001) and the ordinary real space, we can back-transform the deterministic part of Eq. (1) to obtain an expression in terms of perturbation and powering in the simplex: ${\displaystyle M\left({\mathrm{ρ}}_i\right)={\text{ilr}}^{−1}\left({\mathrm{β}}_0+{\mathrm{β}}_1{z}_i+{\mathrm{β}}_2{z}_i^2\right)={\mathrm{γ}}_0⊕\left(z_i⊙{\mathrm{γ}}_1\right)⊕\left(z_i^2⊙{\mathrm{γ}}_2\right),}$ where γ_j = ilr^−1(β_j),  j = 0 , 1, 2. The composition M(ρ_i) is the metric centre (Pawlowsky-Glahn & Egozcue, 2001) of the distribution of ρ_i, an appropriate measure of location for compositions (Aitchison, 1989).

To make the behaviour of the predictions for rare taxa more obvious, we also examined the predictions on a centred logratio (clr) scale, in which the value on the y-axis is the log of the ratio of the corresponding component to the geometric mean of all components (Aitchison, 1986, p. 79). A constant slope on the clr scale corresponds to constant proportional change in the relative abundance of a given taxon. This is also true of the ilr scale, but not of the original proportions. We use the clr scale here because, unlike the ilr scale, it has one coordinate associated with each taxon. For the same reason, clr coordinates are usually chosen as rays in a compositional biplot (Aitchison & Greenacre, 2002). However, it is important to remember that slopes on the clr scale are dependent on the set of taxa analyzed. In addition, although there are s clr coordinates, points in the clr space are constrained to lie in an (s − 1)-dimensional hyperplane in which the sum of the coordinates is zero. This means, that, for example, covariance matrices in the clr scale are singular (Aitchison, 1986, pp. 78–81).

Comparison with non-metric multidimensional scaling

We contrasted our approach with what is likely to be the most popular alternative in marine ecology, a non-metric multidimensional scaling of the raw counts. We used the metaMDS() function in R package ‘vegan’, with default options (square root transformation, Wisconsin double standardization, Bray–Curtisrtis dissimilarity). For comparison, we plotted the first two principal components of the posterior mean still-specific predictions in ilr coordinates.

Alternative models

We also considered multinomial regression fitted by penalized likelihood using ‘glmnet’ (Friedman, Hastie & Tibshirani, 2010), and two naive models that are easy to fit: overdispersed Poisson regression using HMSC (Ovaskainen et al., 2017), which does not respect the multinomial sums, and multivariate linear regression on ilr-transformed counts with the addition of three different kinds of pseudocount (Martín-Fernandez, Palarea-Albaladejo & Olea, 2011). For details, see Section S7.

Community dissimilarity

As described above, most of the common measures of dissimilarity between communities are not perturbation invariant. In the Aitchison geometry, the obvious perturbation invariant measure of difference between two s-part compositions is the Aitchison distance (the Aitchison norm of the compositional difference), defined by ${\displaystyle d_a\left({\mathrm{ρ}}_1,{\mathrm{ρ}}_2\right)=âˆ\yen {\mathrm{ρ}}_1⊖{\mathrm{ρ}}_2{âˆ\yen }_a={\left[\underset{i=1}{\overset{s}{∑}}{\left(log\frac{{\mathrm{ρ}}_{1,i}}{g\left({\mathrm{ρ}}_1\right)}−log\frac{{\mathrm{ρ}}_{2,i}}{g\left({\mathrm{ρ}}_2\right)}\right)}^2\right]}^{1/2}={\left[\underset{j=1}{\overset{s−1}{∑}}{\left(x_{1,j}−x_{2,j}\right)}^2\right]}^{1/2}}$ (Aitchison, 1992; Egozcue et al., 2003), where g(ρ_k) denotes the geometric mean of the parts of a composition, and x_k,j denotes the jth ilr coordinate of x_k = ilr(ρ_k), k = 1, 2. The last line gives the Aitchison distance as the Euclidean distance in ilr coordinates (Egozcue et al., 2003). It is immediately obvious that the Aitchison distance is perturbation invariant, because (a⊕ρ₁)⊖(a⊕ρ₂) = ρ₁⊖ρ₂, by the associative, commutative and identity properties of the vector space. Under this approach, the dissimilarity between the expected compositions ρ₁, ρ₂ is given by (2)${\displaystyle âˆ\yen {{\mathrm{ρ}}_1⊖{\mathrm{ρ}}_2âˆ\yen }_a=âˆ\yen {\left[{\mathrm{γ}}_0⊕\left(z_1⊙{\mathrm{γ}}_1\right)⊕\left(z_1^2⊙{\mathrm{γ}}_2\right)\right]⊖\left[{\mathrm{γ}}_0⊕\left(z_2⊙{\mathrm{γ}}_1\right)⊕\left(z_2^2⊙{\mathrm{γ}}_2\right)\right]âˆ\yen }_a=|z_1−z_2|âˆ\yen {{\mathrm{γ}}_1⊕\left[\left(z_1+z_2\right)⊙{\mathrm{γ}}_2\right]âˆ\yen }_a,}$ using the identity, commutative, associative and distributive properties of the vector space to simplify.

The Aitchison distance has a biological meaning in terms of population growth. In temporal comparisons, the Aitchison distance between two sets of relative abundances is proportional to the among-taxon standard deviation of proportional population growth rates (Spencer, 2015). In spatial comparisons, we can therefore think of the Aitchison distance as measuring the among-taxon variability in proportional population growth rates that is needed to transform one set of relative abundances into another, over a given time interval. This property is important because in a closed system, population growth is the only way to transform one set of relative abundances into another. No other measure of community dissimilarity has this interpretation.

The simplex with Aitchison geometry is a normed vector space (Egozcue et al., 2003). Thus |ρ₁⊖ρ₂|_a = 0 if and only if ρ₁⊖ρ₂ = 0, where 0 is the identity element in the simplex (e.g. Horn & Johnson, 1985, p. 259). From Eq. (2), assuming that γ₁ ≠ 0 and γ₂ ≠ 0, this happens when either z₁ = z₂ (the two compositions are at the same depth) or ${\mathrm{γ}}_2=\left(−\frac{1}{z_1+z_2}\right)⊙{\mathrm{γ}}_1$ (the coefficient of squared depth is a powering of the coefficient of depth). Thus, if we plot dissimilarity on a grid of depths, there will always be zeros on the main diagonal, because communities at the same depth have the same expected composition. There may also be communities at different depths with the same expected composition, along a counter-diagonal where centred and scaled depth has a constant sum, but only in the special case where γ₂ is a powering of γ₁ (or equivalently, where β₂ is a scalar multiple of β₁ in ilr coordinates).

We calculated posterior distributions of dissimilarities among 100 equally-spaced expected compositions between the minimum and maximum depths, both including and excluding bare wall. We plotted the posterior mean dissimilarity matrix, and the widths of the 95% highest posterior density intervals. We only report the results including bare wall here, because those excluding bare wall were very similar. Note that it is valid to exclude some parts of the composition if necessary, because the subcompositional coherence property means that such exclusion will not affect relationships among the remaining parts (Aitchison, 1994).

Rate of change of community composition with depth

The community is changing rapidly with respect to depth if a small increase in depth leads to a large difference in composition. In order to correctly evaluate this change, we need an appropriate definition of difference in composition. Given the geometry of the simplex, the difference in composition between depths z and z + h is naturally expressed as f(z + h)⊖f(z). Then letting h go to zero leads to the obvious definition of the derivative D^⊕f of a simplex-valued function f, ${\displaystyle D^{⊕}\mathbf{f}\left(z\right)=\underset{h→0}{lim}\left(\frac{1}{h}⊙\left(\mathbf{f}\left(z+h\right)⊖\mathbf{f}\left(z\right)\right)\right),}$ provided this limit exists (Egozcue, Jarauta-Bragulat & Díaz-Barrero, 2011, section 12.2.2). Using the rules for differentiation of simplex-valued functions (Egozcue, Jarauta-Bragulat & Díaz-Barrero, 2011, section 12.2.2), in our model, the derivative of expected community composition M with respect to depth, at a depth of z, is ${\displaystyle D^{⊕}M\left(z\right)={\mathrm{γ}}_1⊕\left(2z⊙{\mathrm{γ}}_2\right).}$ This is itself a composition. If we want a scalar measure of rate of change, the obvious choice is the norm of this derivative. It is intuitively obvious that the usual Euclidean norm is not appropriate, because the zero element for compositions (with all parts equal, corresponding to no change in composition with respect to depth) does not have zero Euclidean norm. Instead, we use the Aitchison norm |D^⊕M(z)|_a (Egozcue et al., 2003), which is zero in the situation where there is no change in composition with respect to depth, and is used in the definition of a limit in the simplex (Egozcue, Jarauta-Bragulat & Díaz-Barrero, 2011, Definition 12.2.1). The easiest way to think of this norm is that it is equal to the Euclidean norm of the derivative in isometric logratio coordinates. We evaluated the posterior distribution of this scalar measure of rate of change at 100 equally-spaced depths over the observed depth range.

It is important to remember that we are measuring proportional change: doubling of relative abundance means the same thing whether the initial relative abundance is low or high. This is an essential property, because relative abundances have meaning only in relative terms. In addition, an increase in relative abundance of a taxon may occur in several different ways. For example, the absolute abundance of a taxon may increase while absolute abundances of other taxa remain constant, or the absolute abundance of a taxon may decrease while absolute abundances of other taxa decrease more. In compositional data analysis (and in ecological situations where the focus is on relative abundances), these situations are equivalent.

In order to show how the compositional approach leads to different results from widely-used approaches in ecology, we plotted Bray–Curtis dissimilarities between adjacent predicted compositions (on a grid of 100 equally-spaced depths) against depth (Section S8). This gives a rough estimate of the relationship between rate of change in community composition and depth, because the depth intervals are small. In order to show that this is a potentially general result, we performed a similar analysis for the mite data set of Borcard, Legendre & Drapeau (1992). We fitted a compositional regression model with linear effects of substrate density and water content, with the same multinomial observation model as for the marine community data, and plotted Bray–Curtis dissimilarities between adjacent predicted compositions at equally-spaced values of each explanatory variable, with the other variable held constant (Section S9).

Depth-integrated relative abundances

Over a vertical slice from surface to bottom, a taxon that has high relative abundance over a small range of depths may be unimportant compared to a taxon that has moderate relative abundance at all depths. We therefore want some measure of the “mean” relative abundances over a vertical slice. The arithmetic mean is not appropriate for compositional data. For example, with a banana-shaped distribution, the arithmetic mean may lie completely outside the cloud of observations. The metric centre is a more appropriate measure of the centre of a compositional distribution which avoids these problems (Aitchison, 1989). However, taking a sample estimate of the metric centre over all depths is problematic when there are zero counts. Zeros are difficult to deal with in compositional data analysis (Martín-Fernandez, Palarea-Albaladejo & Olea, 2011), and in this context, will lead to the estimate of the centre being undefined. In addition, if the depth distribution of samples is not uniform, the sample estimate of the centre will be biased. Thus, integrating the model-estimated composition over the full range of depths may be a better way to summarize the structure of the community.

The mean of a real function f of one variable over the interval [a, b] is ${\displaystyle \frac{1}{b−a}{∫}_a^{b}f\left(x\right)dx,}$ which can be thought of as the value of the constant function whose integral over [a, b] is the same as that of f over the same interval (Riley, Hobson & Bence, 2002, pp. 73–74). If we treat community composition as a simplex-valued function of depth, then the analogous mean of this function over the full range of depths gives the composition representing the relative abundance of each part over a vertical slice from top to bottom of the dock wall. Let [S, D] be the depth range, from shallow to deep. Using the rules for integration of simplex-valued functions (Egozcue, Jarauta-Bragulat & Díaz-Barrero, 2011, section 12.3.2), the required mean value is ${\displaystyle \frac{1}{D−S}⊙{\left[\left(z⊙{\mathrm{γ}}_0\right)⊕\left(\frac{z^2}{2}⊙{\mathrm{γ}}_1\right)⊕\left(\frac{z^3}{3}⊙{\mathrm{γ}}_2\right)\right]}_S^D.}$ We evaluated the posterior distribution of this mean value.