Presence-only Habitat Suitability Model; Ecological Niche Factor Analysis

Ecological Niche Factor Analysis (ENFA) is a multivariate empirical approach to study geographic species distributions. It does not require absence data. The working principle is based on the procedures of:

- Summarizing all variables into a few uncorrelated, ecologically relevant factors, and - Computing suitability functions by comparing environmental variable values of species presence cells with respective mean values of the entire study area.

It is built on the concept of two fundamental assumptions: Marginality and specialization.

If the species distribution mean differs from the global distribution mean (m_s≠mG) (see Figure 31), this is called the marginality (M). Formally it can be shown by the mathematical equation 1.

Figure 30: Variable of distance to streambed (m) was normalized by using the Box-Cox algorithm in the BioMapper software 4. The left figure represents the distribution before the transformation and the right one the resulting histogram after the Box-Cox transformation.

Distance (m) Transformed_ Distance (m)

(eq.1) Where M = The marginality of a species

ms = The mean of the species distribution mG = The mean of the global distribution σG = The variance of the global distribution

A large value (close to one) of the marginality means that the species lives in a very particular habitat in the reference study area. Division by σ_G is needed to remove any bias introduced by the variance of the global distribution. The coefficient weighting (1.96) *σ_G assures that the marginality value lies between zero and one.

A middle value (close to 0.5) denotes species habitat which is not too different from the mean condition of the reference area. But, a larger value (close to 1) means that the species has a particular habitat preference regarding the reference area. This equation (1) mainly Figure 31: Marginality and specialization value represented for one variable. The dark area means the species distribution on that variable whereas the blue area represents the distribution for the whole set of cells. The difference in distribution means of a variable for species presence cells (m_s) and the global set of landscape cells (m_G), quantifies the species marginality. Specialization is the ratio of standard deviation of the global distribution σ_G to that of the species distribution σs (Modified from Hirzel et al., 2002).

σ

σ

m

m

Value of Ecogeographical variable

Frequency of cells

explains the principle of the method. The operational equation of marginality to be implemented using BioMapper software is as follows:

(eq. 2) Where M = The overall marginality to compare species‟ marginalities within different

study areas

mi = The marginality of the focal species on each EGV, in units of standard deviations of the global distribution

The higher the coefficient values of an EGV the further the species departs from the mean available habitat regarding the corresponding variable (Hirzel et al., 2002). Negative coefficients on the marginality factor express that smaller values of an EGV are preferred by a species whereas positive coefficients shows a species preference for higher values of the corresponding EGV.

Specialization defines how much different is the variance of the EGV values which can be found in species presence location than the global variance; it is known as the ratio of the variance of the global distribution (σ_G ) to that of the focal species (σ_s). It can be expressed by Equation 3:

(eq. 3) S = The specialization of a species

σ_G = The variance of global distribution σ_s= The variance of species distribution

The higher the specialization factor the stronger is the contribution of that EGV to species specialization. Equation (3) mainly expresses the principle of the ENFA method. The optional definition of specialization implemented in the BioMapper software is:

(eq.4)

Where S = Overall specialization (range from 1 to infinity), V = The number of EGVs

i = The ratio of the variance of the global distribution (σGi) to that of the species distribution (σsi) for any EGV condition in the model

The larger the global specialization value becomes the narrower the species niche (Bryan and Metaxas, 2007). Both global marginality and specialization values depend mainly on the reference area of the study (Derek et al., 2009).

Figure 32 (a) represents a 3 dimensional EGV space. The larger ellipsoid (yellow balloon) represents a global distribution of 3 EGVs, whereas the small violet balloon is the subset of cells of 3 EGVs at which the focal species was detected. The straight line is drawn running through the centre of the two ellipsoids and then it passes the global distribution mean and species distribution mean (μG) and species distribution mean (μs). The species marginality is the difference between global distribution mean and species distribution mean. To extract the specialization factors, two ellipsoids were projected onto a plane perpendicular to the marginality factor for changing the ellipsoids into a sphere (Hirzel et al., 2002).

Orthogonal to the marginality factor, a first specialization factor can be produced as uncorrelated factor by computing the axis that maximizes the ratio of the variance of the global distribution (yellow) to that of the species distribution (violet) (see Figure 32-b).

The other uncorrelated specialization factors were produced by extracting subsequently and restored each EGV, describing how specialized the focal species is in the available

Figure 32: Geometrical interpretation of Ecological Niche Factor Analysis (Hirzel ,2005).

(a). Extraction of marginality factor (b). Extraction of specialization factors (Modified from Hirzel et al., 2002).

(a)

(b)

condition of habitat in the study area. The successive specialization factors are ordered by decreasing coefficient value. Hence, most of the information is retained in the first few factors (Hirzel et al., 2002).

The ENFA model normally applied Idrisi raster maps which are grids and have continuous values. Each cell of a map contains the value of one variable. Before conducting the ENFA, all the EGV maps are normalized as far as possible. The species Boolean map is used to link EGVs in the analysis. To avoid model overfitting because of the large number of EGVs and to assure model reliability, Hirzel et al. (2008) suggest to categorize EGVs into groups such as land use related features, geographical features, etc. ENFA can be computed separately group by group, keeping the best EGVs from each model run. The outputs of ENFA are:

a). A score matrix (cf. Table 15) which is ranking the environmental variables based on their importance for habitat selection in a study area. It can give the information of species-environment relationship by means of marginality and specialization values. In the rows of the score matrix the EGV contributions (variable coefficients) to each factor are given.

Table15: Score matrix sorting the EGVs by decreasing coefficient values of the marginality factor. The coefficient values on the marginality and specialization factors provide the basis for the ecological interpretation of species-habitat relationships.

EGVs Factors of Marginality and Specialization Factor 1 Variable 1 Coefficient value11 Coefficient value21 --- Coefficient valuen1

Variable 2 Coefficient value₁₂ Coefficient value₂₂ --- Coefficient value_n2

--- --- --- --- ---

Variable n Coefficient value1n Coefficient value2n --- Coefficient valuenn

Global

Marginality ---- Specialization ---- Tolerance:1/S ----

In the score matrix, coefficient values of the ecological niche factors explain how marginal and specialized the species are in terms of the various relevant EGVs (Hirzel et al.,2002).

The first factor explains 100% of the marginality and it may also explain some amount of specialization. The next factors take account only for specialization. The coefficients‟ signs have meanings only for the marginality factor. These signs have no interpretation for specialization. A negative sign indicates a species‟ preferences for low value of the

respective EGV whereas a positive sign indicates a preference for a higher value. A high value of global marginality (M) means the species range is different from average conditions of all EGVs. The species‟ tolerance is measured by the inverse of the specialization factors (Sattler et al., 2007). A low value of tolerance (close to 0) indicates that the species is bound to a narrow niche whereas a high value (close to 1) means the species accepts a wide spectrum of habitat conditions (Hirzel et al., 2002). Habitat suitability of any cell for the global distribution is calculated by the first few important factors, accounting for 100% of marginality and some proportion of specialization. The best EGVs are determined by the highest coefficient values on marginality and specialization. The final ENFA model can be summarized by extracting the variables of highest scores.

b). The Habitat suitability (HS) map gives an area-wide prognosis of habitat quality/species spatial distribution. Hirzel et al. (2002) described standard robust methods to compute the suitability for the cells of the whole study area for the focal species. The detailed explanation of habitat suitability computation can be found in the published main ENFA paper of Hirzel et al. (2002). Habitat suitability maps created in the BioMapper software are based on four different habitat suitability algorithms, namely median, distance geometric mean, distance harmonic mean and minimal distance algorithms. Out of these four, the median algorithm is recommended to be used in the type of non-systematic species distribution data (Hirzel, 2004). Before BioMapper 3.0, median algorithm was the only available. It gives good results in most situations and can process quicker than the others. The other three algorithms have no assumption regarding the distribution of species points and are based on functions of the distance between the species occurrences in the environmental space. But in the case of small sample size, the Harmonic mean algorithm should be taken into account to get better results rather than the other three ones (Hirzel, 2004).

This study relies on a small number of species presence points. The Harmonic mean was suited for that small sample to create the tiger habitat suitability map (Hirzel, 2004). This algorithm is commonly used to define home ranges and activity centres from detection locations (Dixon and Chapman, 1994) in the geographic space. The function of this algorithm is:

(eq. 5)

Where H = The harmonic mean

P = Species‟ observation points

N = N-dimensional environmental space (the number of EGVs) Oi = The harmonic mean of the distances of all observation points.

The effect of this mean algorithm is to give a (too) high weight to all observations while keeping the information of observation density in the factor space. Therefore, it has a tendency to overfit the data, which might be good when in case of small sample sizes (Hirzel, 2004).

Figure 33: Computing habitat suitability by using the median algorithm; the farther the location (arrow) is from the median (dotted line), the lower its suitability (Hirzel et al., 2003). HS of any cell for the whole area is calculated from its location (arrow) relative to the species distribution (dark green) (Braunisch et al., 2008). The global suitability is derived by computing a weighted mean on these "partial suitabilities" (Modified from Hirzel et al., 2002).

Σ = (1/2) suitability y_index

Median

Ecological Niche Factor Classes

Frequency of cells

Focal class