Generalized Linear Models outperform commonly used canonical analysis in estimating spatial structure of presence/absence data

Baseline datasets

We compared the two approaches to spatial variable selection using simulated community data based on three real community composition datasets with a range of properties:

1. Presence/Absence of 110 marine benthic macroalgae species from a Rapid Assessment Program for biodiversity of 42 sample sites spanning roughly 2,000 km² at Ilha Grande Bay, Rio de Janeiro, Brazil (southwest Atlantic) (Carlos-Júnior et al., 2019), permit number IBAMA/RJ:031/04);

2. Presence/Absence of 588 plant species from grassland covering 500 km² of Scotland’s coast. Data were collected from 3639 5 × 5 m quadrats from 94 sites. We used sites as our sample units, treating species as present when they occurred in at least one quadrat at a site, and absent otherwise (see Lewis, Pakeman & Marrs, 2014) for more information);

3. Presence/Absence of 47 freshwater aquatic insect species collected from 30 sample sites in five tributaries of the Guapiaçú River basin, Brazil which covers about 40 km² (Feijó-Lima in prep, permit number INEA-RJ: 019-2014).

For each of the datasets we used the geographical coordinates (maps and sampling sites in Fig. S1) to calculate spatial explanatory variables for regression (Fig. 1). We chose MEMs as our spatial variables since they are commonly used to describe spatial structure in ecological studies. Moreover, in contrast to coarser methods such as trend-surface analysis, MEMs are a flexible method, capable of describing all spatial scales provided by the sampling design (Borcard, Gillet & Legendre, 2011). They are also more flexible and powerful than the method of principal coordinates of neighbor matrices (PCNMs, a special case of distance-based MEMs) (Bauman et al., 2018a; Bauman et al., 2018b; Borcard & Legendre, 2002; Dray, Legendre & Peres-Neto, 2006). One needs two matrices to build the MEM variables for a given set of site coordinates: matrix B describing the connectivity among the geographical sampling sites and matrix A describing the weights of such connections. The Hadamard product of these two matrices generates the spatial weighting matrix (matrix W), which is then doubly centred and diagonalized, yielding eigenvectors to be used as spatial variables. For ecological studies, the processes of interest are usually those generating positive autocorrelation, and it is therefore common to use only MEMs associated with positive eigenvalues (as in this study). For studies in which negative spatial autocorrelation is also of interest (e.g., where negative interactions such as competitive exclusion, predation, etc are suspected), the eigenvectors associated with negative eigenvalues can also be separately used (Bauman et al., 2018a). We made decisions about B and A for each dataset based on our ecological knowledge of the spatial structure of these regions, since our goal was to simulate communities with ecologically sensible spatial structures. Therefore, for dataset 1 we chose the minimum spanning tree (B) with Euclidian linear distances as weights (A). Our decision was based on the shape of the bay and the fact that the main water movements make the sampling sites geographically compartmentalised in subregions where sites are likely to be minimally connected (Carlos-Júnior et al., 2019). Similarly, spatial organisation in dataset 2 could be sensibly described in terms of Delaunay triangulation (B) with Euclidian weights (A). Despite some degree of connectivity among all sites, pairs of sites could be mostly associated not to their immediate neighbours but rather as a function of their distances. This is due to cultural differences in land management. For example, northern and western islands share cultural histories, which is reflected in species composition (Lewis, Pakeman & Marrs, 2014). Directional spatial processes in ecological data, such as those observed in rivers, are well described by a special case of MEMs called asymmetric eigenvector maps (AEM, Blanchet, Legendre & Borcard, 2008b), which were used for constructing variables for dataset 3. In MEMs, larger eigenvalues are associated with broader-scale spatial structures while smaller eigenvalues represent fine-scale spatial structures. This allowed us to control the spatial scale of variation in community structure. Dataset 1 had 16 positive MEMs from 42 sites, dataset 2 had 30, and dataset 3 had 12 AEM variables with positive autocorrelation. For computation of the MEMs for the three datasets we used the packages adespatial (version 0.3–7, Dray et al., 2019) and spdep (version 0.7–4, Bivand & Piras, 2015; Bivand, Hauke & Kossowski, 2013).

Simulating communities with chosen spatial drivers

We simulated realistic communities with known spatial structure, based on the three datasets. We used spatial eigenvectors as explanatory variables. We varied the number of MEMs with non-zero coefficients and created new binary (presence/absence) communities (with the same number of sites and same expected number of species as the real ones) using two different modelling scenarios. These simulated communities reflected the effect of those MEMs with non-zero coefficients. By varying the number and ordering of the non-zero coefficients, we could therefore control the spatial structure and scale of the simulated community data (see scheme in Fig. 1 and Table 1).

In order to simulate new binary communities under the simulated presence method (SPM, in which species are always detected if present), we first estimated a coefficient matrix C of size (m variables + 1 (first) row with intercepts) × p species from each real data set. This was achieved using the manyglm function with binomial errors in R package mvabund (version 3.11.9, Wang et al., 2012), with explanatory matrix X (n sites × m positive MEMs + an initial column of 1′s). The matrix C gives the effect of each explanatory variable on the logit-transformed probabilities of presence. The mvabund package provides a GLM framework for multivariate response data.

We then created new hypothetical scenarios by generating a new coefficient matrix C^∗, of the same size as C, whose elements $c_{kj}^{\ast }$ are given by (1)${\displaystyle \left\{\begin{array}{c}{c}_{kj}^{\ast }=c_{1j},\text{if}k=1,j=1,2,\dots ,p,\left(\text{intercepts}\right)\hfill \\ c_{kj}^{\ast }\sim {\hat{F}}_b,\text{if}k-1\in K,j=1,2,...,p\hfill \\ c_{kj}^{\ast }=0,\text{otherwise},\hfill \end{array}\right.}$ where ${\hat{F}}_b$ is the empirical distribution function of c_kj (k =2, 3, …, m +1, j = 1, 2, …, p) (Evans, Hastings & Peacock, 2000), and the $b_{kj}^{\ast }$ are sampled with replacement. The set K defines to which rows of C^∗ the non-zero coefficients were allocated: we studied 14 such sets (see below and Table 1(A–C)). In other words, we used the originally-estimated intercepts in each simulation (first row of Eq. (1)), and drew those coefficients assigned to non-zero values (second row of Eq. (1)) from the empirical distribution of all the originally-estimated explanatory variable coefficients. We sampled the values of the non-zero coefficients from the empirical distribution in order to simulate plausible but not fixed spatial structures. Table 1 depicts for each dataset how the non-zero coefficients were assigned for each dataset and simulation scenario (see below).

We then calculated predicted probabilities of presence ${\hat{p}}_{ij}$ for the jth species at the ith site. Given the matrix $\hat{\mathbf{Y}}=\mathbf{X}{\mathbf{C}}^{\ast }$ (n sites × p species) of predicted logit probabilities of presence, the predicted probability of presence is (2)${\displaystyle {\hat{p}}_{ij}=\frac{exp\left({\hat{y}}_{ij}\right)}{1+exp\left({\hat{y}}_{ij}\right)}.}$

The simulated presence/absence value for species j at site i was sampled from a Bernoulli distribution with success probability ${\hat{p}}_{ij}$. The result is a community matrix with the same number of sites and the same expected number of species as the real community, and with realistic coefficients for spatial eigenvectors. As in the maximum likelihood estimation done by manyglm (Wang et al., 2012), species and sites were assumed conditionally independent when generating simulated presence/absence data, given the values of the explanatory variables. Our simulated communities correspond to the simple case in which presence/absence patterns are affected by environmental variables but not interspecific interactions. Nevertheless, interspecific interactions could well be relevant to real world systems and other models (Godsoe & Harmon, 2012; Anderson, 2017).

Since GLMs are specified correctly for presence/absence data generated this way, we would expect them to perform well. We therefore devised a second ecologically meaningful simulation method in which absences arise from the sampling protocol, called the simulated abundance method (SAM). The two simulation methods differ in whether they assume we have true absences or sampling-related absences. Note that it is not possible to simulate binary data directly using RDA, because RDA does not generate predicted probabilities of presence. Instead, we treated $\hat{\mathbf{Y}}$ as log expected abundances and exponentiated each element to get expected abundances λ. Then we calculated the probability of detecting the species under Poisson sampling (i.e., the probability of drawing a value of at least 1 from a Poisson distribution with parameter λ), which is (3)${\displaystyle {\hat{p}}_{ij}=1-e^{-\lambda }.}$

Finally, we generated a Bernoulli random variable with success probability ${\hat{p}}_{ij}$ to produce a simulated presence-absence observation. Both GLM and RDA are mis-specified for data generated in this way. Codes for both the SPM and SAM simulation frameworks and all the datasets used in our simulations are available as supplemental information (Data S1–Data S3).

We compared GLM and RDA variable selection under up to 14 different scenarios, differing in the number of non-zero coefficients (nVar) and whether these coefficients were associated with fine or broad spatial scales. We simulated up to six different choices of the number of MEM variables creating the spatial structure in the data (i.e. having non-zero coefficients): none, approximately one sixth, approximately one third, approximately half, approximately three-quarters, and all (Table 1(A–C), rows). We also simulated three different spatial scales of the patterns. As mentioned above, MEMs associated with larger eigenvalues represent broader spatial scales. We ordered the MEMs in descending order of eigenvalues and arranged the non-zero coefficients within matrix C^∗ in three different ways (Table 1 (A–C), columns): only broad-scale MEMs with non-zero coefficients (scaling 1); only fine-scale MEMs with non-zero coefficients (scaling 2); half broad-scale, half fine-scale (scaling 3). Because not every combination of number of non-zero coefficients and spatial scaling is possible (e.g., it is not possible to assign one non-zero coefficient in scaling 3), there were 14 possible combinations overall for each dataset (Table 1). The main steps of the simulation scheme are summarized in Fig. 1.

RDA and GLM

We used the default RDA function from the R package vegan (version 2.5-6, Oksanen et al., 2019), with simulated community composition as the response variable, and MEMs associated with positive eigenvalues generated from geographical coordinates of the sample sites as explanatory variables. In order to perform a transformation-based RDA (Borcard, Gillet & Legendre, 2011; Blanchet et al., 2014) we used the Ochiai coefficient, which is the Hellinger transformation analogue for binary data, as recommended by Legendre & Gallagher (2001) and Borcard, Gillet & Legendre (2011).

Binomial GLMs were fitted to the same data using the manyglm function in R package mvabund (Wang et al., 2012). We fitted our models using a logistic regression (logit link function for binomial response), with species compositional data as the multivariate response variable and MEMs as predictors. No interaction terms were included, following common practice in spatial modelling of community data.

Comparing model selection between RDA and GLM frameworks

We compared the results of model selection between the approach usually taken in the RDA and a somewhat-similar approach for GLMs. For RDA, we used the forward selection with double stopping criterion following (Bauman et al., 2018a; Bauman et al., 2018b), beginning with a global test of significance (model with all spatial predictors) and carrying on with the variable selection if the global model was significant. The forward selection itself consists of a stepwise procedure including in the model the variable contributing the most to the adjusted R². The procedure stops either when the next variable with the highest contribution is not significant (first stopping criterion) or causes the adjusted R² to be bigger than that of the global model (i.e., containing all variables; second criterion). This is implemented in the function ordiR2step in the vegan package (Oksanen et al., 2019). For GLM, we used forward selection with a stopping rule based on minimum Akaike Information Criterion (AIC) (Akaike, 1973; Wagenmakers & Farrell, 2004). The selection procedure started from a model with intercept only and added one explanatory variable at a time, until no further improvement in the sum of AIC over each of the response variables was possible. We used this approach because the usually large number of MEMs makes it difficult to compare the AIC sum over all possible GLMs.

The performance of each method on simulated data was mainly assessed by two criteria. First, we assessed how many MEMs with zero coefficients were incorrectly included in the final model. Second, we assessed how many MEMs with non-zero coefficients were incorrectly excluded from the final model. Also, we assessed overall accuracy (score) as the percentage of MEMs whose inclusion/exclusion status was correct. The goals of ecological studies are usually not directly related to the inclusion/exclusion of individual MEM variables, but instead to identify spatial pattern, represented by a linear combination of MEMs. However, since the MEMs form a basis for the space spanned by the transformed spatial weighting matrix, such a linear combination is unique (Fraleigh & Beauregard, 1995, pages 197-198). Furthermore, the MEMs are orthogonal, so that each represents a qualitatively distinct aspect of spatial pattern. Therefore, if an individual MEM is incorrectly included or excluded, the estimated spatial pattern is qualitatively wrong.

We further explored the ability of each method to capture spatial pattern using a graphical approach (Article S1). For each real dataset and each method, we haphazardly picked one simulated data set. We plotted the MEM decompositions of both the true and estimated spatial patterns. We chose the scenarios in which each method had the worst performance in terms of correctly including/excluding variables, in order to determine whether in such cases, overall spatial pattern would still be captured.

Finally, we calculated how much of the variation in response variables was explained by each method using the adjusted R² for the linear model in RDA and its analogue for GLMs, the D-value (Tjur, 2009). These two values cannot be directly compared since they are not exactly equivalent, but their results could yield interesting insights and are made available as supplemental information (see table results in Data S4).

For each of the combinations of conditions in Table 1, 1,000 simulated data sets were generated under each of SPM and SAM. For each simulated data set, spatial explanatory variables were selected using both GLM/AIC and RDA/FW.