Mahalanobis distances for ecological niche modelling and outlier detection: implications of sample size, error, and bias for selecting and parameterising a multivariate location and scatter method

Introduction

The Mahalanobis distance (Mahalanobis, 1936) is a statistical technique that can be used to measure how distant a point is from the centre of a multivariate normal distribution. Mahalanobis distances are commonly applied to problems such as classifying data into groups and determining differences between groups (Manly, 2005). Mahalanobis distances have also been used to quantify species-environment relationships through habitat and ecological niche models (Dettmers, Buehler & Bartlett, 2002; Johnson & Gillingham, 2005; Tsoar et al., 2007; Etherington et al., 2009). In this context Mahalanobis distances are classified as a presence-only technique because they do not require species absence or background environmental data (Peterson et al., 2011) and simply require a data matrix ${\displaystyle \mathbf{A}=\left[\begin{array}{ccccc}\hfill x_{1,1}\hfill & \hfill x_{1,2}\hfill & \hfill x_{1,3}\hfill & \hfill \dots \hfill & \hfill x_{1,n}\hfill \\ \hfill x_{2,1}\hfill & \hfill x_{2,2}\hfill & \hfill x_{2,3}\hfill & \hfill \dots \hfill & \hfill x_{2,n}\hfill \\ \hfill x_{3,1}\hfill & \hfill x_{3,2}\hfill & \hfill x_{3,3}\hfill & \hfill \dots \hfill & \hfill x_{3,n}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill x_{m,1}\hfill & \hfill x_{m,2}\hfill & \hfill x_{m,3}\hfill & \hfill \dots \hfill & \hfill x_{m,n}\hfill \end{array}\right]}$ that has m rows of species occurrences for which n columns of environmental variables have been obtained. The multivariate location and scatter of the data are defined by an n-dimensional vector $\hat{\mu }$ containing the sample means for each column of variables, and a sample variance-covariance matrix $\hat{\Sigma }$ of dimensions n × n that contains variances for each column along the main diagonal and pair-wise column covariances values elsewhere (Manly, 2005).

The Mahalanobis distance (1)${\displaystyle D^2\left(\mathbf{x}\right)={\left(\mathbf{x}-\hat{\mu }\right)}^{\text{T}}{\hat{\Sigma }}^{-1}\left(\mathbf{x}-\hat{\mu }\right)}$ can then be calculated for any vector x = [x₁, x₂, x₃, …, x_n] that represents a position in environmental space as defined by the n environmental variables. As D² is essentially the sum of n independent standard normal variables, the D² values from a multivariate normal population will follow a chi-squared distribution with degrees of freedom equal to the number of dimensions n (Manly, 2005). This means that a chi-squared cumulative distribution function $F_{{\chi }_n^2}\left(x\right)$ can be used to convert D² into a probability $P\left({\chi }_n^2\le D^2\right)$ that indicates if a location in environmental space has a D² that is greater than would be expected by chance (Etherington, 2019). For example, when applied to m = 20 hypothetical points in n = 2 dimensions, the calculated $P\left({\chi }_n^2\le D^2\right)$ values follow a characteristic elliptical pattern centred at the mean of each environmental variable, with $P\left({\chi }_n^2\le D^2\right)$ increasing from the centre outwards in a manner that accounts for the variability within and correlation between each environmental variable (Fig. 1).

The elliptical form of Mahalanobis distances fits well with the theoretical concept of the fundamental niche, which Hutchinson (1957) p. 416 defined as “an n-dimensional hypervolume …which corresponds to a state of the environment which would permit the species …to exist indefinitely”. Hutchinson (1957) used a rectangular model to define environment limits of the fundamental niche, but also stated that “If the variables are independent in their action on the species we may regard this area as the rectangle …but failing such independence the area will exist whatever the shape of its sides”. So any convex shape, which includes the elliptical shape of $P\left({\chi }_n^2\le D^2\right)$, would be an appropriate model of the fundamental niche. Indeed, we see the use of ellipses alongside other convex shapes in later development of the niche concept (Hutchinson, 1978). In this context $\hat{\mu }$ represents the optimal environmental conditions at the centre of the fundamental niche, and $\hat{\Sigma }$ represents both the range of and interaction between environmental conditions within the fundamental niche.

While $P\left({\chi }_n^2\le D^2\right)$ is usually inverted to $P\left({\chi }_n^2>D^2\right)$ for use in ecological niche modelling to estimate the probability of an environmental location being within a fundamental niche (Etherington, 2019), $P\left({\chi }_n^2\le D^2\right)$ is also commonly used in statistics and data science to detect outliers (Aggarwal, 2017). However, when applied to detect outliers, the use of sample means and variance-covariance to estimate D² can be problematic because these measures of multivariate location and scatter are sensitive to outliers (Fig. 1). So for outlier detection, D² can be calculated with different methods for defining the multivariate location and scatter of data, such as the minimum covariance determinant (MCD) and minimum volume ellipsoid (MVE), which are much more insensitive to outliers (Rousseeuw, 1985).

The MCD approach estimates multivariate location ${\hat{\mu }}_{\text{MCD}}$ and scatter ${\hat{\Sigma }}_{\text{MCD}}$ from a subset numbering h data points that has the smallest variance-covariance matrix determinant (Hubert & Debruyne, 2010). The MVE approach is similar to the MCD in that it works with a subset of size h data points, but the MVE estimates multivariate location ${\hat{\mu }}_{\text{MVE}}$ and scatter ${\hat{\Sigma }}_{\text{MVE}}$ from the ellipsoid of minimal volume that encapsulates the h data points (Van Aelst & Rousseeuw, 2009). D² can then be calculated using either the MCD measures of multivariate location and scatter (2)${\displaystyle D^2\left(\mathbf{x}\right)={\left(\mathbf{x}-{\hat{\mu }}_{\text{MCD}}\right)}^{\text{T}}{\hat{\Sigma }}_{\text{MCD}}^{-1}\left(\mathbf{x}-{\hat{\mu }}_{\text{MCD}}\right)}$ or the MVE measures of multivariate location and scatter (3)${\displaystyle D^2\left(\mathbf{x}\right)={\left(\mathbf{x}-{\hat{\mu }}_{\text{MVE}}\right)}^{\text{T}}{\hat{\Sigma }}_{\text{MVE}}^{-1}\left(\mathbf{x}-{\hat{\mu }}_{\text{MVE}}\right)}$ for any vector x = [x₁, x₂, x₃, …, x_n] that represents a position in n-dimensional environmental space.

The utility of the MCD and MVE methods for defining the multivariate location and scatter of data is highly relevant for ecological niche modelling based on digitally mobilised data through data-sharing networks such as the Global Biodiversity Information Facility (Edwards, Lane & Nielsen, 2000). These networks are reliant on natural history observation data that are likely to contain errors such as taxonomic misidentification or incorrect and imprecise georeferencing (Graham et al., 2004), which can result in species occurrences that are outliers in environmental space. Returning to our hypothetical example (Fig. 1), we can see that both the MCD and MVE methods ignore the apparent outlier and produce measures of multivariate location and scatter that are focussed on where the data are more concentrated. Given the robustness of MCD and MVE methods to outliers resulting from errors in natural history observation data, it is perhaps no surprise that both methods have been adopted recently for ecological niche modelling (Norris, Jackson & Betancourt, 2006; Liu, White & Newell, 2018; Soberón, Peterson & Osorio-Olvera, 2018; Yañez-Arenas et al., 2018; Qiao et al., 2019; Altamiranda-Saavedra et al., 2020; Osorio-Olvera et al., 2020; Castaño-Quintero et al., 2020).

On the other hand, natural history observation data will also contain sampling bias (Graham et al., 2004) that can be extreme and result in greater amounts of data for more charismatic species, more accessible places, more developed countries, and more recent times (Meyer, Weigelt & Kreft, 2016). If unaccounted for this could skew models results. Therefore, returning to our hypothetical example (Fig. 1), we may also have a situation in which the apparent outlier only appears to be an outlier due to sampling bias at environments with lower temperatures and rainfall. If this were the situation then the MCD and MVE methods would be providing a poorer estimate of multivariate location and scatter because they are both focussing on the data bias, whereas the method based on the sample mean and covariance provides a better estimate that relates to all the data points.

Unfortunately, the differing effects of errors and bias means there is unlikely to be a best method for estimating the multivariate location and scatter of data in all situations, but this is consistent with ecological niche modelling more generally (Qiao, Soberón & Peterson, 2015). Therefore, this paper uses a virtual ecology approach (Zurell et al., 2010) to simulate a series of ecological niche modelling experiments to understand when different methods of defining the multivariate location and scatter of the data are more appropriate as sample size, errors, and bias vary.

Materials & Methods

The virtual species used in the experiments is the Antipodean opaleye dragon, which is imagined to live in the mountain valleys of New Zealand (Scamander, 2001). I have defined the species’ fundamental niche in terms of population growth rates (Maguire, 1973) measured as the finite rate of increase λ_F of a population in a two-dimensional environmental space of temperature and rainfall. The fundamental niche is defined by three parameters: λ_max, the maximum finite rate of increase at the fundamental niche optimum; μ, a 2 × 1 column vector of means that gives the optimal temperature and rainfall condition; and Σ, a 2 × 2 variance-covariance matrix that determines the size and orientation of the fundamental niche in each dimension. With ${\displaystyle {\lambda }_{\text{max}}=2.5\text{,}\mu =\left[\begin{array}{c}\hfill 7.5\hfill \\ \hfill 1800\hfill \\ \hfill \hfill \end{array}\right]\text{, and}\Sigma =\left[\begin{array}{cc}\hfill 2\hfill & \hfill -950\hfill \\ \hfill -950\hfill & \hfill 800000\hfill \\ \hfill \hfill \end{array}\right]}$ the fundamental niche finite rate of increase λ_F(x) of a virtual species for any vector x of environmental space coordinates is then calculated as (4)${\displaystyle {\lambda }_F\left(\mathbf{x}\right)={\lambda }_{\text{max}}\times e^{-\frac{1}{2}{\left(\mathbf{x}-\mu \right)}^{\text{T}}{\Sigma }^{-1}\left(\mathbf{x}-\mu \right)}}$ that results in an elliptically shaped fundamental niche (Fig. 2A).

Generating samples of species occurrences begins with an idealised virtual sampling of the niche space. An initial set of sampling locations S = {x₁, x₂, …, x_i} consisted of a series of locations in environmental space. S was randomly generated from a two-dimensional normal distribution $S\sim {\mathcal{N}}_2\left(\mu ,\Sigma \right)$ with the same vector of means μ and variance-covariance matrix Σ as defined the fundamental niche. Then the probability that a sampled location x_i resulted in the species being both present and detected P_d(x_i) was described as a logistic function (5)${\displaystyle P_d\left({\mathbf{x}}_i\right)=\frac{1}{1+e^{-10\left({\lambda }_F\left({\mathbf{x}}_i\right)-0.5\right)}}}$ which was parameterised so that P_d(x_i) increases as λ_F(x_i) increases, but with P_d(x_i) ≈ 0 where λ_F(x_i) ≈ 0 because the virtual species population is unlikely to exist under these environmental conditions, and P_d(x_i) ≈ 1 where λ_F(x_i) ≳ 1 because above this population growth rate the population should always be present.

Using this idealised sampling, a sample of m = 100 occurrence locations can be produced that are concentrated towards the centre of the fundamental niche, and for which $P\left({\chi }_2^2>D^2\right)$ estimates for all three methods of multivariate location and scatter align with the elliptical shape of the fundamental niche (Fig. 2A). Under these idealised sampling conditions we can see that the actual λ_F values and the estimated niche probabilities $P\left({\chi }_2^2>D^2\right)$ for the occurrence samples are very highly correlated, and that all three methods of determining the multivariate location and scatter of the niche perform equally well (Fig. 2B).

Of course this idealised virtual sampling is completely unrealistic, but does serve to demonstrate that it is possible to produce useful models with enough accurate data. What is of interest is understanding under what conditions of sample size, error, and bias the predictive ability of the various multivariate location and scatter methods begin to break down. Therefore, to explore this, a series of virtual experiments were conducted.

All experiments were done using R (R Core Team, 2019) with the MASS (Venables & Ripley, 2002), virtualNicheR (Etherington & Omondiagbe, 2019), fields (Nychka et al., 2017), raster (Hijmans, 2020), and extrafont (Chang, 2014) packages.

Sample size

It is common for many species to have as few as seven unique occurrence locations within the Global Biodiversity Information Facility (Meyer, Weigelt & Kreft, 2016). However, as the MCD approach needs m ≥ n × 5 (Hubert & Debruyne, 2010), and given the experiments are two-dimensional, the smallest value of m that could be analysed is 10. Therefore, to explore the effects of sample size, the idealised sampling approach was applied, but varying the sample size m from 10 to 160 in increments of 15. As m increases, the idealised samples increasingly represent the elliptical shape of the fundamental niche (Figs. 3A–3C)

Regarding the choice of sample subset h used by the MCD and MVE methods, the standard choice is h = ⌊(m + n + 1)∕2⌋ for both the MCD and MVE methods, because this produces the most robust estimates (Hubert & Debruyne, 2010; Van Aelst & Rousseeuw, 2009). This standard choice uses just over half the sample, so for sample sizes from 10 to 160 this would mean a standard choice of h from 6 to 81 as the sample size increases. Because h varies as a function of sample size, for experimental consistency and to aid interpretation I specified the sample subset used by the MCD and MVE methods as a proportion k of the sample size, such that h = ⌊k × m⌋. The standard choice of h was represented by k = 0.55, and with k = 0.75 and k = 0.95 used to explore the effect of increasing k, and therefore h, on the performance of the MCD and MVE methods.

Each sample of size m was replicated 500 times, and was applied to each multivariate location and scatter method, with the MCD and MVE methods applied at the three different k values. The performance of each multivariate location and scatter method was measured as the correlation between the actual λ_F and estimated $P\left({\chi }_2^2>D^2\right)$ values.

Results

There were obvious consistencies and trends amongst all the experiments. First, the results for the MCD and MVE methods were very similar, making it hard to differentiate the performance of these two methods (Figs. 4–6). Second, in all cases, as k increases the MCD and MVE methods become more similar to the sample means and variance-covariance method (Figs. 4–6).

Regarding the effects of sample size, the sample means and variance-covariance method performed better, but this difference only became notable when m ≲ 50 and was less pronounced as k increased (Fig. 4). The fact that m = 100 gives good results regardless of the method is important to recognise, because as we can be confident that any performance effects in the error and bias experiments that used m = 100 will be a function of the imposed error or bias rather than the sample size.

Considering errors, in all cases the MCD and MVE methods performed better than the sample means and variance-covariance method, though this difference was only evident at lower error levels for higher k values (Fig. 5). The pattern of response to bias was more subtle, with the sample means and variance-covariance method performing better when bias was ≳50%, but with this effect being very prominent when k = 0.55, less obvious when k = 0.75, and no longer present when k = 0.95 (Fig. 6).

Discussion

Conclusions

When using D² for ecological niche modelling and outlier detection, the performance of multivariate location and scatter methods varies based on sample size, error, and bias. Comparison of the sample means and variance-covariance, MCD, and MVE multivariate location and scatter methods provides the clear conclusion that none of the methods, or individual parameterisations of any methods, can be considered universally the best. Rather, any ecological niche modeller using these techniques needs to think carefully about their data and objective to choose the method and parameterisation that are most appropriate to their individual circumstances. For those modellers who wish to explore the potential of the MCD and MVE methods, given these methods have been used widely in statistical analyses for some time, these methods should be widely available in statistical software. However, modellers using the MCD and MVE methods should carefully consider and clearly state the h or k parameter used in their analyses.

Supplemental Information