Predictive modelling and ground validation of the spatial distribution of the New Zealand long-tailed bat (Chalinolobus tuberculatus)

Predictive modelling and ground validation of the spatial distribution of the New Zealand long-tailed bat (Chalinolobus tuberculatus)

B I O L O G I CA L C O N S E RVAT I O N 1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1 available at www.sciencedirect.com journal homepage: www.elsevier.com/locate...

560KB Sizes 0 Downloads 25 Views

B I O L O G I CA L C O N S E RVAT I O N

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/biocon

Predictive modelling and ground validation of the spatial distribution of the New Zealand long-tailed bat (Chalinolobus tuberculatus) Glen J. Greavesa, Renaud Mathieub,*, Philip J. Seddona a

Department of Zoology, University of Otago, P.O. Box 56, Dunedin, New Zealand School of Surveying, University of Otago, P.O. Box 56, Dunedin, New Zealand

b

A R T I C L E I N F O

A B S T R A C T

Article history:

The use of predictive models is continually increasing, but few models are subsequently

Received 26 March 2005

field-checked and evaluated. This study evaluates the statistical strength and usefulness

Received in revised form

for conservation purposes of a predictive habitat use model developed for Chalinolobus

27 March 2006

tuberculatus, a threatened microchiropteran bat species found in the temperate rainforests

Accepted 2 April 2006

of New Zealand. The relationship between various environmental variables and the pres-

Available online 12 June 2006

ence/absence of the species was investigated using generalised linear modelling. The model developed was coupled with GIS data to develop maps of predicted occurrence

Keywords:

within the West Coast region of New Zealand’s South Island. It was found that distance

Bat conservation

to forest boundary, slope, presence of Nothofagus, general land cover, variability in mean

Habitat selection

annual solar radiation, and mean ambient winter minimum temperature were significantly

Chalinolobus tuberculatus

associated with the occurrence of the species. Evaluation of the statistical strength of the

GIS modelling

distribution model with independent data of species’ occurrence collected at 152 sites

New Zealand

found that the C. tuberculatus model showed a moderate ability to predict both species presence and absence (s(b) coefficient = 0.37). The field detection rate (0.45) using this model was significantly higher than that of historical surveys (0.12). The value of the species habitat model and the need to evaluate its utility in the development of conservation strategies is discussed.  2006 Elsevier Ltd. All rights reserved.

1.

Introduction

Recognition of the factors that influence zoogeographical patterns is undoubtedly critical for the success of conservation management (Wang et al., 2003). Hence, mapping of species–habitat relationships has long been a key provider of graphical information for wildlife conservation projects. Yet, the applicability of traditional map products is limited by their two-dimensional, ‘static’ nature. The introduction of digital maps and Geographic Information Systems (GIS) tech-

nology dramatically changed the capacity of spatial data representation, as such enhancing their utility in predicting the distribution of species (Greenberg et al., 2002). Its development allowed wildlife biologists to analyse variables at larger spatial scales, with greater precision and accuracy than previous methods of spatial analysis. GIS-derived modelling of species distribution has been successfully applied to many animal and plant groups (Martı´nez-Salvador et al., 2005; Powell et al., 2005), including insectivorous bats (Jaberg and Guisan, 2001; Wang et al., 2003).

* Corresponding author: Tel.: +64 3 479 7698; fax: +64 3 479 7586. E-mail addresses: [email protected] (G.J. Greaves), [email protected] (R. Mathieu), philip.seddon@ stonebow.otago.ac.nz (P.J. Seddon). 0006-3207/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.biocon.2006.04.016

212

B I O L O G I C A L C O N S E RVAT I O N

However, a conspicuous feature of many of these earlier models is their lack of independent assessment of model performance. Instead, researchers have used alternative techniques based on resampling of input data (Seoane et al., 2004) (e.g. bootstrapping, resubstitution, randomisation, or partitioning), or often undertaken no substantiative testing at all (Fielding and Bell, 1997; Johnson et al., 2004). Although computer-intensive validation procedures have some merit in simulating species occurrence, they fail to provide the same degree of confidence as using an independent dataset (Pearce and Ferrier, 2000). In our view, the best means of objectively assessing model performance is to use an independent set of locality records and quantitative accuracy measures. Many recent studies have achieved this by withholding a sub-sample of the original dataset, to which validation tests are applied (Boyce et al., 2002). However, when the original dataset is limited in size and sufficient resources are available, validation of the model through comprehensive field survey is considered the best method. New Zealand’s unique endemic flora and fauna has suffered greatly following the arrival of humans and their associated exotic species introductions (Atkinson, 2001; Dowding and Murphy, 2001; Fitzgerald and Gibb, 2001). This has resulted in a highly disproportionate number of threatened species, extending widely across plant and animal groups (Atkinson, 2001; Dowding and Murphy, 2001; Nugent et al., 2001; Towns et al., 2001). It is only the relatively recent date of this colonisation and the success of intensive management, particularly the relocation of remnant populations to predator-free offshore islands, that has saved the majority of faunal species from extinction (Armstrong and Ewen, 2001; Towns and Ferreira, 2001). The country’s only endemic terrestrial mammals, three species of microchiropteran bat, have not escaped the onslaught, but until recently did not receive the attention afforded to the suite of threatened bird, reptile, and amphibian species. For the Greater Short-tailed bat (Mystacina robusta), this attention has likely come too late (the last recorded sighting was 1967). Despite being fully protected under the Wildlife Act (1953) in response to their significant post-settlement reduction in abundance and distribution (Daniel, 1990; O’Donnell and Sedgeley, 1994; Molloy, 1995; O’Donnell, 2000), New Zealand’s two remaining bat species – Long-tailed bat (Chalinolobus tuberculatus) and Lesser Short-tailed bat (Mystacina tuberculata) – continue to decline. The cause of this decline is uncertain, but is likely to be a combination of habitat loss, predation by introduced mammals, and competition for food and roosting sites by introduced mammals, birds, and wasps (O’Donnell, 2000, 2001a). Both C. tuberculatus and M. tuberculata are under serious threat of following M. robusta into extinction in the medium term (O’Donnell, 2000). Before the future management of New Zealand’s bat species can be effectively planned, a thorough evaluation of their status is required. Improved knowledge of the distribution of bats will enable sites to be selected for ongoing management, monitoring, and advocacy (Molloy, 1995). However, the existing methods of surveying for New Zealand bats have a limited efficiency, with a detection rate of only 12% for C. tuberculatus surveys within the current study area (West Coast Bat Database, Department of Conservation, DOC). The primary pur-

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

pose of this research was to improve this efficiency by modelling the likely spatial distribution of bats, thus allowing concentration of surveys within areas where, on the basis of habitat suitability, bats are more likely to occur. Historical surveys have also focused on sites believed to be suitable, but the definition of suitability has been anecdotal. There is a need for more rigorous assessment of habitat suitability, using an objective measure of resource use. The implicit assumption with animal distribution surveys is that suitable habitat is saturated with individuals, i.e. the habitat is at carrying capacity. Therefore, the first reason why we may not detect a bat in suitable habitat relates to the non-detection of bats at a site despite their presence (Type II error), giving a probability of detection that is less than one. Without specific modelling the actual probability of detection is unknown, however we can assume that it is constant for both historical and validation surveys. The second reason may be attributed to an imperfect knowledge of what constitutes suitable habitat for a bat, so surveying for bats would occur in places where they have chosen not to be. This is the key error as it dictates the improvement of an objective spatial model over a subjective historical method, and as such, presents the two key questions to be answered by this research: (1) is C. tuberculatus more likely to be detected in sites considered suitable by the predictive distribution model, compared with sites the model deems unsuitable? (2) is the predictive distribution model a better basis for identifying suitable bat habitat for surveys than the historical approach? This second question is particularly critical as it provides the justification for the construction of, and is therefore a direct measure of the success of, the predictive distribution model.

2.

Materials and methods

2.1.

Study area

This research took place within the New Zealand Department of Conservation’s West Coast Tai Poutini Conservancy (Fig. 1). Several minor and localized extensions to the eastern boundary were added to incorporate a few significant bat surveys that bestride the conservancy boundary. Located on the west coast of the South Island of New Zealand, the West Coast Conservancy manages more than 1.9 million hectares of land. The conservancy contains within its boundaries two national parks, parts of three others and a World Heritage Area. Spanning over 600 km in length, the West Coast region can be regarded as an ‘ecological island’, confined by the Southern Alps to the east, and the Tasman Sea to the west. At its widest, the distance from the sea to the mountain tops (which regularly peak over 3000 m) is 90 km, and is often much less. As such, there is marked diversity in habitat type from coastal temperate rainforest through to a true alpine zone, and short transition distances between these distinct biomes. Indigenous forest covers about three quarters of the region, characterised by a mixture of Podocarp/hardwood at low altitudes, gradually giving way to Nothofagus dominant forest from

B I O L O G I C A L C O N S E RVAT I O N

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

213

Fig. 1 – Study area, West Coast Conservancy, South Island, New Zealand. The general focal areas for the predictive model validation are displayed (Karamea, Paparoa, Maruia, Shenandoah).

about 500 m above sea level (McGlone et al., 1996; Leathwick, 1998).

2.2. 2.2.1.

Development of the predictive distribution map Bat field data

The database used for this research was sourced from the Department of Conservation and comprised bat presence/absence records obtained between 1994 and 2003, derived from a combination of stationary and transect surveying. Where more than one survey was done at a particular site, the surplus surveys were removed from the analysis to avoid pseudo replication. To reduce the influence of false negative (Type II) error on the bat presence analysis, all ‘absence’ sites within 100 m of a ‘presence’ site were removed. Likewise, presence sites within 100 m of another presence site were removed. Due to the mobility of C. tuberculatus, it was considered highly likely that bats would also be located in habitat within the 100 m radius around a presence site, and may be detected given further surveying effort. Moreover, measuring habitat selection at smaller scales than this is beyond the scope of this study. In addition to the removal of these replicate sites, a number of data were excluded due to recording inaccuracies found following assessment of geographical co-ordinates. As

a result, of the over 2500 records available, a total of 1033 records were deemed suitable for input into the model.

2.2.2.

Environmental data

Environmental descriptors were derived from three databases – Land Environments New Zealand (Ministry for the Environment), New Zealand Forestry Service Maps (Series 6) (New Zealand Forestry Service), and Land Cover Database Version 1 (Ministry for the Environment). While the Land Environments New Zealand database was provided in raster format (100 · 100 m grid), the Land Cover Database and New Zealand Forestry Service database were provided in vector format and therefore were converted to raster for ease of analysis, with a discrete code allocated to each variable. A total of 14 potential environmental predictors were selected following analysis of current knowledge of C. tuberculatus ecology (O’Donnell, 1999, 2001a) and habitat data availability. Table 1 provides a summary of the variables investigated. Twelve of the 14 variables were extracted directly from the databases, however, ‘distance to forest boundary’ and ‘presence of beech (Nothofagus)’ descriptors were created within ArcGIS (Environmental Science Research Institute Inc.) using data from Land Cover Database Version 1 and New Zealand Forestry Service Maps (Series 6) respectively. The environmental value allocated to each bat

214

B I O L O G I C A L C O N S E RVAT I O N

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

Table 1 – Source, scale, and estimated errors of habitat variables chosen for inclusion into the regression analysis Habitat variable Vegetation Variables

Underlying Variables

Source

Scale

Errora

I II III IV V VI VII

Land cover Forest type Presence of beech Area of indigenous forest patch Area of forest patch Perimeter of forest patch Distance to forest boundary

LCDB1 NZFS6 NZFS6 LCDB1 NZFS6 LCDB1 LCDB1

1/50,000 1/50,000 1/50,000 1/50,000 1/50,000 1/50,000 1/50,000

6.1% Insignificant Insignificant 6.1% Insignificant 6.1% 6.1%

Climate

VIII IX X XI XII XIII

Mean annual temperature Mean winter minimum temp Mean annual solar radiation Mean winter solar radiation October vapour pressure deficit Monthly water balance ratio

LENZ LENZ LENZ LENZ LENZ LENZ

14 stations 14 stations 7 stations 7 stations 14 stations 65 stations

±0.3–0.4 C ±0.4–0.6 C ±0.21–0.27 MJ/m2/day ±0.15–0.55 MJ/m2/day ±0.04–0.6 kPa ±0.512

Physical

XIV

Slope

LENZ

1/50,000

NA

LCDB1 = Land-Cover Database Version 1; NZFS6 = New Zealand Forestry Service Maps (Series 6); LENZ = Land Environments New Zealand. a Estimated error (LCDB1, NZFS6) or standard error (LENZ).

sample was calculated as the average of values contained within a 100 m radius buffer where the measured variable was continuous, and the value which contributed the largest area within the same buffer where the variable was discrete.

2.2.3.

Spatial accuracy of data

Spatial accuracy of the bat location data, both in the source dataset and the field validation phase of this research, may be compromised by error in Global Positioning System (GPS) plotting, or in data transcription. GPS error is associated with factors such as signal degradation, signal multipath, poor signal reception, or atmospheric delays. Prior to 2000 GPS signals were artificially degraded through United States Defence applied Selective Availability (SA), resulting in an error of up to 100 m. As the majority of the data used in the model was sourced between 1996 and 2000, prior to the removal of SA, a 100 m buffer to all data sites was applied. Transcription error was effectively removed by cross-validation of all coordinates with in-depth site location description. The two potential errors associated with categorical environmental datasets such as those used in this research may be related to either: (1) positional accuracy, or (2) semantic accuracy. For Land Cover Database Version 1, the satellite imagery was orthorectified to produce a positional accuracy of ±25 m (Terralink International Ltd.). The semantic accuracy of this dataset was stated as 90% by Terralink International Ltd. In an accuracy assessment undertaken by New Zealand Forest Research in 2000, overall map accuracy was estimated at 93.9% using a simple accuracy percentage statistic. To minimise errors in the New Zealand Forestry Service digitisation process, unused and unfolded source maps were utilised, and any significant linear differences and incorrect polygon tagging were noted and corrected before the final database was released. The map producers (New Zealand Forestry Service), and therefore the authors of this paper, consider the error value to be insignificant. The Land Environments New Zealand database was the source for the seven underlying climate and physical habitat characteristics used in this research. As this database predicts the variation in continuous environmental variables across

space (interpolated from a number of measuring stations), each of which is displayed as discrete layers, the associated errors occur within each independent variable and not over the database as a whole. For the majority of layers, the potential standard error, for the whole of New Zealand, is stated by the database source (Ministry for the Environment). The two temperature variables (Table 1) use a total of 300 meteorological stations nationwide from which the layers are derived, of which only 14 are located within the study area. This low number for the West Coast study area is reflected in the standard error, which, ranging from 0.30 to 0.40 C (mean annual temperature; mean = 10.1 C, SD = 2.03 C) and 0.40–0.60 C (mean winter minimum temperature; mean = 0.86 C, SD = 3.0 C), is slightly higher than the majority of the remainder of New Zealand (Leathwick et al., 2002). The potential errors associated with the two Land Environments New Zealand variables estimating solar radiation are proportionally smaller than that for the temperature variables discussed above (mean annual solar radiation estimated error = 0.210–0.270 MJ/m2/day, mean = 14.2 MJ/m2/day, SD = 0.95 MJ/m2/day; mean winter solar radiation estimated error = 0.35–0.55 MJ/m2/day, mean = 5.6 MJ/m2/ day, SD = 0.92 MJ/m2/day) (Leathwick et al., 2002), representing 0.017 and 0.080 for annual solar radiation and winter solar radiation measures, respectively. Standard errors associated with the predictions of October vapour pressure deficit are predominantly below 0.05 kPa within the study area (mean = 0.33 kPa, SD = 0.12 kPa), and significantly lower in the mid to low altitudes occupied by C. tuberculatus. This corresponds to a potential error of 15% around the mean, which is high compared with the other variables included for analysis. No assessment has been made of the errors associated with the Land Environments New Zealand slope layer (Leathwick et al., 2002). However, mean elevation error for the Digital Elevation Model (DEM) was found to be 0.41 m (standard error: 6.13 m) using over 2500 independent GPS data points taken from throughout the South Island (Landcare Research, unpublished data).

2.2.4.

Statistical analysis and spatial prediction

This research used a resource selection function (RSF) model to characterise the selection of resources by C. tuberculatus in the

B I O L O G I C A L C O N S E RVAT I O N

West Coast region of New Zealand. RSF models are defined as any function that is proportional to the probability of use of a resource unit by an organism (Manly et al., 1993), and are often fitted using generalised linear models (GLM) (Boyce et al., 2002). As the bat location data was in the form of presence/absence, logistic regression was considered the most appropriate method of generating the RSF model (Guisan and Theurillat, 2000; Rushton et al., 2004). Broadly applied in ecological research, modern regression approaches have been found to be particularly useful in modelling the spatial distribution of species (Osborne and Tigar, 1992; Guisan et al., 2002), including bats (Jaberg and Guisan, 2001). Logistic regression uses a linear combination of independent variables to explain the variance of a dependant binary. It applies maximum likelihood estimation after transforming the dependant into a logit variable. In this way, logistic regression estimates the probability of a certain event occurring (Eq. (1)). This research used hypothesised habitat relations (14 independent variables), based on the known ecology of C. tuberculatus, to describe its known presence/absence (dependant binary), thereby providing a prediction of the likelihood of bat occurrence at any 100 · 100 m cell within the study area.



eðaþb1 x1 þb2 x2 þþbi xi Þ 1 þ eðaþb1 x1 þb2 x2 þþbi xi Þ

ð1Þ

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

215

where H is the likelihood of bat occurrence, a is the constant of the equation, x is the predictor variable, b is the logit coefficient of the predictor variable. The potential error increasing effects of multicollinearity between the 14 variables was screened for by computing the squared multiple correlation of each variable, and calculating the Variance Inflation Factor (VIF). In all cases the r-value calculated was less than 0.80, and VIF did not exceed 5, so removal or modification of data was not considered necessary (Montgomery and Peck, 1982). To determine which combination of variables best explains the data in a RSF model some form of information criteria is used, such as Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) (Anderson and Burnham, 1998; Boyce et al., 2002). Based on the principle of parsimony, these criteria help to determine which model provides the best compromise of high explanatory factor and minimal descriptor inclusion. AIC was considered the better of the two criteria because of the small sample size (Hastie et al., 2001). The best GLM, as selected by AIC, enabled the creation of the predictive distribution map (Fig. 2). Following the acquisition of predictive distribution variables from the forward-stepwise algorithm, ArcGIS, ARCInfo (Environmental Science Research Institute Inc.), and MapInfo (MapInfo Corp.) software was used to process and map the predictive

Fig. 2 – Predictive distribution map for C. tuberculatus within the West Coast Conservancy.

216

B I O L O G I C A L C O N S E RVAT I O N

bat occurrence within the study area. The application of the logistic regression equation (Eq. (1)) to the raster calculator function of ArcGIS (using the logit coefficient output from the forward-stepwise algorithm), yielded a map that displayed the predicted likelihood of bat occurrence for each 100 · 100 m cell of the GIS grid (Fig. 2).

2.3.

Model validation

2.3.1.

Field surveys

The predictive distribution map was evaluated using standard New Zealand bat surveying techniques (O’Donnell and Sedgeley, 2001). Surveying took place during the months of October and November, the period of high bat activity corresponding to the highest ambient temperatures (O’Donnell and Sedgeley, 2001). Surveys did not take place during periods of extreme weather, such as heavy rain or strong winds to minimise the bias due to reduced detectability of bats under such conditions. Surveying was restricted to four selected sectors of the study area – Karamea, Shenandoah, Maruia, and Paparoa (Fig. 1) – due to accessibility and high number of areas predicted as being suitable for bats. Effort was further focused to testing only those areas attributed with extreme suitability values of 0–0.1 (predicted as unsuitable), and 0.9–1.0 (predicted as suitable). The presence of C. tuberculatus at previously unsurveyed sites was tested using a combination of stationary (42%) and mobile (transect walking) (58%) Batbox III detectors (Stag Electronics, Sussex, United Kingdom), the standard techniques for obtaining a simple presence/absence measure (O’Donnell and Sedgeley, 1994, 2001). Stationary detectors have the ability to continuously survey one site for many hours, but cannot detect bats outside their maximum range (a Batbox III at full volume can detect C. tuberculatus on average 43.5 ± 9.8 m away when recorded along the forest margin (O’Donnell, unpublished data)), and are therefore spatially constrained. By comparison, transect walking allows the coverage of large areas, but relatively little time spent at each site. Therefore a combination of the two methods in one study area allows the undertaking of thorough surveys over large tracts of land. Temperature and weather measures were used to avoid those survey periods when bat activity was highly unlikely, such as dusk temperature below 5 C, or heavy rain. In all cases, the Batbox III detector was set at 40 kHz, the frequency accepted as the optimal for the surveying of C. tuberculatus (O’Donnell and Sedgeley, 2001). Stationary boxes were left on site for a minimum of three suitable nights, as dictated by the maximum length of battery life. To reduce the likelihood of false detection, a minimum of two passes was required before a site was confirmed as being positive for bat presence. Also, to reduce the risk of the same individual being detected and recorded on more than one occasion, and to standardise the data collected with that used to produce the model, a 100 m survey exclusion zone was applied around each site found to have bats present. Transects of 1 km were walked in about 20 min, and the distance between adjacent transects was determined using computer-derived random numbers. A 100 m buffer was applied where a bat was detected during a particular transect, followed by the use of a random number to determine the distance from the edge of

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

Table 2 – A confusion matrix (modified from Fielding and Bell, 1997) Actual (field) Predicted (model) + 

+



a c

b d

the buffer to the start of the next transect. Each transect was walked four times, twice each night over two nights.

2.3.2. Validation statistical analysis 2.3.2.1. Model accuracy. The overall accuracy of a probabilistic model derived from presence/suspected absence data can be assessed using a confusion matrix (Fielding and Bell, 1997) (Table 2). The accuracy measures utilised in this research (correct classification rate, positive and negative predictive powers, and s(b) coefficient (Eqs. (2)–(5))) were chosen for their ability to evaluate the predictive accuracy of the models. The unbalanced nature of the dataset (1:6.5 presence versus suspected absence ratio in the West Coast Bat Database), presence of zeros in the matrix, and small sample size precluded the use of some commonly used measures including odds ratio, and j coefficient. Ma and Redmond (1995) found the s coefficient, which measures the improvement of a classification over a random assignment of pixels to groups, to better adjust percentage agreement than j, and easier to calculate and interpret. Correct classification rate ¼ ða þ dÞ=N

ð2Þ

Positive predictive power ¼ a=ða þ bÞ

ð3Þ

Negative predictive power ¼ d=ðc þ dÞ sðbÞ coefficient T ¼ p0  pr =1  pr

where pr ¼ 1=N

ð4Þ M X 2

ni xi ð5Þ

i¼1

where p0 is the correct classification rate; M is the number of categories of the model classification (presence/absence); i is the ith category; N is the total number of sites; ni is the row total for category i; xi is the diagonal value for category i (i.e. number of correct assignments for habitat i); a, b, c, d refer to Table 2.

2.3.2.2. Comparative success. v2 analysis was used to establish whether the distribution models provide a better basis for identifying suitable bat habitat for surveys than the historical approach. In this example, the null hypothesis is that there is no significant difference between the predictive distribution model and historical techniques in the ability to determine a site’s suitability as a bat foraging zone. The expected frequency was calculated as the number of sites surveyed divided by the historical success rate, with a predetermined cut-off probability of 0.05 applied. SPSS (Version 13.0) (SPSS Inc.) was the statistical software package used throughout this research.

3.

Results

3.1.

Model construction

Six habitat variables were found to improve the accuracy of the best model significantly (p < 0.10) in predicting bat pres-

B I O L O G I C A L C O N S E RVAT I O N

217

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

Table 3 – Effects of each habitat variable on the likelihood of C. tuberculatus detection, as assessed by fitting a logistic regression model to known bat presence/suspected absence bat location data Step

1 2 3 4 5

6

Variable sequentially added

Mean winter minimum temp Mean annual solar radiation Presence of Nothofagus Slope Land cover 1. Urban 2. Pastoral 3. Indigenous 4. Scrub 5. Inland water 6. Inland wetland 7. Mines/Dumps 8. Planted forest Distance to forest boundary

b

Wald’s test p-value

0.067 0.265 1.766 0.096 3.044 1.076 1.214 4.118 0.641 1.802 1.186 11.007 0.115

AIC Value

D AIC

Prop

0.000 0.000 0.000 0.000 0.061

750.511 710.644 674.425 649.690 637.948

– 39.867 36.219 24.735 11.742

– 0.056 0.051 0.037 0.019

0.002

628.502

9.446

0.015

With b = logit coefficient of the predictor variables, Walds p-value = probability value associated with the Wald statistic, DAIC = change in AIC, found by difference between previous and present AIC values, Prop = Change in AIC value expressed as a proportion of the AIC value.

ence/suspected absence (Table 3). Applying this 6-step model to the source bat presence/suspected absence data, it correctly predicted the occurrence of bats for 76% of the validation sites. The key environmental predictors of C. tuberculatus distribution, as extracted by the GLM analysis, indicate a preference for feeding zones very near the forest boundary (presence  v ¼ 34 m outside of the forest boundary, rv ¼ 2:8 m, n = 137; suspected absence v ¼ 146 m inside of the forest boundary, rv ¼ 9:0 m, n = 896), of low slope (presence  v ¼ 5:98 , rv ¼ 0:45 , n = 137; suspected absence  v ¼ 7:12 , rv ¼ 0:22 , n = 896), moderate mean minimum winter temperature (presence v ¼ 0:25  C, rv ¼ 0:15  C, n = 137; suspected absence v ¼ 1:00  C, rv ¼ 0:05  C, n = 896), relatively high mean solar radiation (presence v ¼ 137 MJ=m2 =day, rv ¼ 0:43 MJ=m2 =day, n = 137; suspected absence v ¼ 133 MJ= m2 =day, rv ¼ 0:15 MJ=m2 =day, n = 896), and presence of Nothofagus as a major forest constituent (69% of bat ‘presence’ sites).

3.2.

Model validation

Over the 10 week survey period, a total of 961 h of surveying was completed, at 99 different sites located within the West Coast Conservancy study area. A further 24 sites (192 h) surveyed were removed from the following analysis due to system failure, adverse weather conditions, or external interference. Of the 80 sites predicted as having C. tuberculatus presence, 36 were confirmed as such over the survey period, giving an overall success rate of 0.45 (Table 4). No bats were detected at sites where they were not expected (19 sites, 132 h of surveying).

3.2.1.

Assessing model accuracy

Applying the unlikely assumptions that (1): C. tuberculatus utilises the entire range of its suitable habitat (i.e. the habitat is at carrying capacity), and (2): the chance of detection at a utilised site equals one; the survey results suggest that

Table 4 – Confusion matrix comparing expected bat presence/absence (based on predicted chance of bat occurrence from distribution model) and observed bat presence/absence (from field surveys) for C. tuberculatus Actual (field)

Predicted (model)

Presence

Absence

Total

Presence Absence

36 0

44 19

80 19

Total

36

63

99

within the total of 99 survey sites, 44 (44%) of the sites were wrongly classified for bat occurrence (for all sites bats were not confirmed to be present where expected). This gives a correct classification rate value of 0.56. Evaluating the models’ ability to predict the occurrence of bats at a site provides predictive power values of 0.45 and 1.0 for presence and suspected absence respectively. This suggests a greater ability to predict the sites of bat absence than presence (although having a smaller sample). Kendall’s s(b) coefficient for this data indicates that 37% more pixels were classified correctly than would be expected by chance alone. This suggests that the model performs to a reasonable standard, but much improvement could be expected through future refinement.

3.2.2.

Comparative success of historical and present surveys

The ability of the distribution model developed in this research to determine suitable feeding habitat (0.45) compares well against the success rate of traditional means of locating such habitat (0.12) calculated, using the data of the DOC bat database, as the ratio number of sites with bat presence/total number of sites. However, as this validation data is drawn from a total of only 80 sites, and surveying effort was focused toward more accessible areas, further testing would be required to establish a more accurate figure.

218

B I O L O G I C A L C O N S E RVAT I O N

4.

Discussion

4.1.

Key environmental predictors

The regression analysis, and subsequent field surveying, suggests that C. tuberculatus feeding habitat is most closely linked to indigenous forest, particularly Nothofagus forest that occurs in areas of high solar radiation, low winter minimum temperatures, and gently sloping topography. However, it must be stressed that the actual differences between the presence and suspected absence datasets of these variables were small. Yet, given that all source data surveying was restricted to areas of indigenous forest, and the clear change in dominant Nothofagus species with subtle environmental variation (Leathwick, 1998), this differentiation may be considered ecologically significant. Nothofagus forest generally out-competes other forest types below a temperature threshold recognised as the lowest for New Zealand forest classes, so it is not surprising that moderately low winter temperature and occurrence of beech are both cited as important bat distribution predictors. Selection of sites with higher solar radiation levels and lower topography may reflect their influence via species differentiation of Nothofagus forest or effects on primary prey species. The West Coast region has five of the total seven Nothofagus species present in New Zealand, the differentiation of which may be determined fairly accurately by a site’s environmental characteristics (Leathwick, 1998; Leathwick and Austin, 2001). An undisturbed West Coast Nothofagus site with a relatively moderate temperature regime, low profile, and high solar radiation levels, will invariably be dominated by N. fusca (McGlone et al., 1996). The inclusion of the distance to forest boundary element would confirm a non-uniform spatial distribution within the forest, tending towards the forest edge zone (±200 m). This zone is characterised by a general absence of medium to tall plant specimens, thereby providing a relatively uncluttered flight-space. Flying macro-invertebrate density also tends to be higher in these areas (Major et al., 2003). For these reasons bats may be able to fulfil their feeding requirements at a lower energetic cost than foraging in the forest interior would allow. As would be expected due to the site selection process, all validation sites had visible environmental characteristics that differentiated them from those sites that were ‘unsuitable’, and reflected the characteristics discussed above. All were typified by being at, or close to, a major forest boundary. The forest, without exception, consisted of pure or co-dominant Nothofagus, and was predominantly in a pristine state with a high density of mature trees. The sites displayed all the attributes of a climax state Nothofagus forest including a relatively open sub-canopy, low species diversity, and thick humus layer. Old-age red beech was present, or in close proximity (within 1 km) at all but two of the sites where C. tuberculatus was confirmed as being present. At both of these non-red beech sites, mature Dacrydium cypressinum were numerous within adjacent forested area, suggesting an alternative roost tree species may be used in the absence of red beech. Although predominantly found roosting in Nothofagus specimens, C. tuberculatus has been frequently found feeding, and roosting in areas where red beech, and the entire Nothofagus genus is absent (O’Donnell, 2001b). A prime example is

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

the population that inhabits the forested shores of Lake Kaniere, Hokitika (West Coast Bat Database, Department of Conservation). Just north of this population marks the northern extremity of the Beech Gap, a narrow zone of about 150 km between mountains and sea where Nothofagus is entirely absent as a constituent of the forest diversity, despite dominating the environment immediately north and south of this zone (McGlone et al., 1996). Analysis has shown the area to have once supported beech in large numbers, and found the current environment to be highly beech compatible (McGlone et al., 1996). It is thought that during the last ice age (20,000– 11,000 B.P.), glacial advances stripped this narrow zone of vegetation, and beech has since failed to reinvade due mainly to its slow propagational properties (McGlone et al., 1996). In its place now stand the faster spreading podocarp, hardwood and softwood species such as D. cypressinum, Dacrycarpus dacrydioides, and Metrosideros umbellata. Outside the forested area the habitat type varied, with riparian, pastoral, coastal, and road-edge environments well represented. As such there was an array of vegetation species present, both native and introduced, as well as non-vegetated land cover (concrete, river, alluvial gravel and boulders, and lake). This would suggest that C. tuberculatus shows little preference towards feeding habitat type, provided it is an open environment in close proximity to a stand of old-growth red beech forest. By comparison, those sites surveyed within areas predicted as unsuitable for bat habitation were predominantly situated greater than 2 km from the nearest forest stand (76%), and generally on pastoral or exotic plantation land (81%). Others were located deep into the forest interior (>2 km), further emphasising the importance of the distance to boundary relationship. These observations strongly support those of past studies on C. tuberculatus habitat use and morphometric analysis. O’Donnell (2001a) and others have stressed the apparent importance of mature N. fusca forest and forest edge habitats. A bat’s flight mode may be deduced from morphological characteristics such as wing shape and echolocation call structure (Norberg, 1987; Webb et al., 1998; O’Donnell, 1999), with these characteristics suggesting C. tuberculatus to be a forest periphery or gap forager.

4.2.

Predictive power of the model

Assessing the accuracy of the C. tuberculatus model in predicting the occurrence of bats within the study area suggests a moderate overall performance. The combination of high negative predictive power (1.0), and modest positive predictive power (0.45), to give a total correct classification rate of sites of 0.55, suggests that the model is significantly better at discerning bat absence than presence. However, these two predictive power values are likely to be moderated by the high ‘false negative’ error rate associated with the data. Survey error (particularly through the detection probability being less than one), and perhaps most importantly, the non-utilisation of suitable habitat by bats due to external factors such as predation, are likely to underestimate significantly the potential presence of bats and therefore undermine the model’s performance. Yet, this result cannot be entirely attributed to the presence of high Type II error; there is likely to be a combina-

B I O L O G I C A L C O N S E RVAT I O N

tion of errors associated with both model construction, and field validation phases of this research. The future addition of key environmental variables absent from this study (primarily through the unavailability of quality data) may significantly improve the accuracy of the predictive distribution model. A measure of forest canopy maturity, direct measure of dominant canopy species, and forest density are likely to be critical factors in determining a sites’ suitability for bat foraging and general habitation. Also likely to provide further improvement in the predictive ability of spatial modelling is the scale at which habitat characteristics are measured, and the number of stations from which each layer is extrapolated. The future refinement of this model as such data becomes available is strongly urged, as is the inclusion of future survey outcomes. Further optimisation of the model validation surveying methods may be achieved by sampling later in the summer season than the present study, when warmer and more stable weather conditions prevail. Also, any increase in the survey effort, both numerically and spatially, would better evaluate the models’ predictive ability. The results from the evaluation of the predictive distribution model suggest that the C. tuberculatus predictive distribution model is a significant improvement over historical techniques, and therefore has immediate application in future surveying for this species. Over 829 h of surveying at 80 sites predicted as ‘suitable’ C. tuberculatus feeding zones, 36 of these sites were found to be frequented by bats (0.45). By comparison, past surveys within the study area, which predominantly used measurable vegetation variables as predictors, exhibited a success rate of 12 bat detections per 100 sites surveyed (West Coast Bat Database, Department of Conservation). Despite this apparent improvement, the rate of detection (0.45) during evaluation surveys must still be considered as being low, and therefore again suggests scope for further refinement of the model. Although extensively used in making spatial predictions from point data, GLM is only one of a number of methods available for use (Yee and Mitchell, 1991; Lehmann et al., 2002). There is considerable debate over the ability of GLM to approximate the true regression surface adequately, particularly when modelling flora (Yee and Mitchell, 1991; Austin and Meyers, 1996). Other popular alternative methods of analysis include the method of least squared (LS), generalised regression analysis and spatial prediction (GRASP), and the generalised additive model (GAM) (Hastie and Tibshirani, 1990), which is an extension of the GLM (Leathwick, 1998; Lehmann et al., 2002). The classical LS regression approach is theoretically valid only when the response variable is normally distributed and the variance does not change as a function of the mean (Guisan and Zimmermann, 2000), so was not applicable to this research. However, some research has found the use of GAM techniques to allow greater flexibility in the fitting of statistical models than that provided by GLM, although these too have their limitations (Yee and Mitchell, 1991; Lehmann et al., 2002). GAMs estimate response curves with a non-parametric smoothing function, instead of the parametric terms required by GLM (Yee and Mitchell, 1991; Leathwick, 1998; Lehmann et al., 2002). This enhanced flexibility has particular importance when the relationships cannot easily be

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

219

explained using parametric models, as commonly occurs with nature-derived data (Leathwick, 1995; Lehmann et al., 2002). In light of this debate the authors suggest testing the appropriateness of a statistical model by running a series of tests (Marshall et al., 1995) or by graphical methods (e.g. diagnostic plots) (Guisan and Zimmermann, 2000).

4.3.

Conservation implications

Threatened species management generally has two broad objectives. In the short term the aim is to minimize the risk of extinction, while the longer term aspiration is to promote conditions in which a species’ viability is apparent without the need for intensive management. A viable population is one that has a high probability of persisting over a long time period (Schaffer, 1981). The analysis of the status of a species population via population viability analysis has become a core tool in conservation practice, applying principles of population ecology to assess their chances of survival (Morris and Doak, 2002). It does this by allowing assessment of vulnerability, and by targeting the research towards factors that may have the greatest influence on extinction probabilities. Of particular relevance to this analysis is population density, which while not directly provided by presence/absence measures, may be estimated through methods such as multiple regression, simple polynomial analyses (Hone, 1999), or power analyses (O’Donnell and Langton, 2003), and species/population geographical range which can only be discerned following intensive and extensive survey. Through the analysis of current population status, and the continued monitoring over time, an indication of trends is gained (Buckland and Elston, 1993), thereby allowing assessment of the long term vulnerability of the species. However, the requirement of quality census data presents a serious challenge for population viability analysis, particularly for endangered species where the data are often plagued by limiting factors such as sampling errors, density dependence, and non-stable age structure, which may mask the true population status (Holmes, 2004). As such, there is a high degree of uncertainty regarding the accuracy of any assessment of a populations’ vulnerability. To minimize the risk of such errors, and therefore maximize the benefit of population viability analysis to conservation management, continued improvement in the accuracy in measuring population characteristics is required. The introduction of GIS-based habitat use analysis into the mapping of species distribution is proving to be of great benefit in determining the spatial extent of target populations, and thereby providing an improving basis for obtaining the dynamic population characteristics required by population viability analysis. This success comes despite limited quality and quantity of input data, an observation further verified by the relative success of the C. tuberculatus model in this research. A more focused observation of C. tuberculatus population trends, permitted through the implementation of the predictive distribution model, will help in the assessment of the immediate and long-term vulnerability, and therefore the extent and immediacy of the management required.

220

B I O L O G I C A L C O N S E RVAT I O N

Acknowledgements The study benefited from the advice and support provided by Dr. Colin O’Donnell (Department of Conservation); Claire Cameron and Brian Niven (Department of Mathematics, University of Otago); the Department of Conservation, Area Managers and support staff, particularly John Lyall (Hokitika Area Office); private landowners; and the many volunteers. This work was supported by Bat Conservation International through a research grant to Glen Greaves and by logistical and financial support provided by the Department of Zoology, University of Otago. Katherine Dickinson and Stuart Parsons commented on an earlier version of this manuscript.

R E F E R E N C E S

Anderson, K.P., Burnham, D.R., 1998. Model Selection and Inference: A Practical Information – Theoretic Approach. Springer, New York. Armstrong, D.P., Ewen, J.G., 2001. Assessing the value of follow-up translocations: a case study using New Zealand robins. Biological Conservation 101, 239–247. Atkinson, I.A.E., 2001. Introduced mammals and models for restoration. Biological Conservation 99, 81–96. Austin, M.P., Meyers, J.A., 1996. Current approaches to modeling the environmental niche of eucalypts: implications for management of forest biodiversity. Forest Ecology and Management 85, 95–106. Boyce, M.S., Vernier, P.R., Nielsen, S.E., Schmiegelow, F.K.A., 2002. Evaluating resource selection functions. Ecological Modelling 157, 281–300. Buckland, S.T., Elston, D.A., 1993. Empirical models for the spatial distribution of wildlife. Journal of Applied Ecology 30, 478–495. Daniel, M.J., 1990. Bats: Order Chiroptera. In: King, C.M. (Ed.), The Handbook of New Zealand Mammals. Auckland University Press, Auckland, pp. 114–137. Dowding, J.E., Murphy, E.C., 2001. The impact of predation by introduced mammals on endemic shorebirds in New Zealand: a conservation perspective. Biological Conservation 99, 47–64. Fielding, A.J., Bell, J.F., 1997. A review of methods for the assessment of prediction errors in conservation presence/ absence models. Environmental Conservation 24, 38–49. Fitzgerald, B.M., Gibb, J.A., 2001. Introduced mammals in a New Zealand forest: long-term research in the Orongorongo Valley. Biological Conservation 99, 97–108. Greenberg, J.D., Logsdon, M.G., Franklin, J.F., 2002. Introduction to Geographic Information Systems (GIS). In: Learning Landscape Ecology, A Practical Guide to Concepts and Techniques. Springer Verlag, New York (Chapter 3). Guisan, A., Theurillat, J.P., 2000. Equilibrium modelling of alpine plant distribution and climate change: how far can we go? Phytocoenologia 30, 353–384. Guisan, A., Zimmermann, N.E., 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135, 147–186. Guisan, A., Edwards Jr, T.C., Hastie, T., 2002. Generalised linear and generalised additive models in studies of species distributions: setting the scene. Journal of Ecological Modelling 157, 89–100. Hastie, T., Tibshirani, T., 1990. Generalised Additive Models. Chapman and Hall, London. Hastie, T., Tibshirani, T., Friedman, J., 2001. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, New York.

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

Holmes, E.E., 2004. Beyond theory to application and evaluation: Diffusion approximations for population viability analysis. Ecological Applications 14, 1272–1293. Hone, J., 1999. Accuracy of multiple regression method for estimating population density in strip transects. Australian Wildlife Research 13, 121–126. Jaberg, C., Guisan, A., 2001. Modelling the distribution of bats in relation to landscape structure in a temperate mountain environment. Journal of Applied Ecology 38, 1169–1181. Johnson, C.J., Seip, D.R., Boyce, M.S., 2004. A quantitative approach to conservation planning: using resource selection functions to map the distribution of mountain caribou at multiple spatial scales. Journal of Applied Ecology 41, 238–251. Leathwick, J.R., 1995. Climatic relationships of some New Zealand forest tree species. Journal of Vegetation Science 3, 601–616. Leathwick, J.R., 1998. Are New Zealand’s Nothofagus species in equilibrium with their environment? Journal of Vegetation Science 9, 719–732. Leathwick, J.R., Austin, M.P., 2001. Competitive interactions between tree species in New Zealand’s old growth indigenous forests. Ecology 82, 2560–2573. Leathwick, J.R., Morgan, F., Wilson, G., Rutledge, D., McLeod, M., Johnson, K., 2002. Land Environments of New Zealand (Ngia Taiao o Aotearoa): A technical guide. Ministry for the Environment, Wellington. Lehmann, A., Overton, J.M., Leathwick, J.R., 2002. GRASP: generalized regression analysis and spatial prediction. Ecological Modelling 157, 189–207. Ma, Z.K., Redmond, R.L., 1995. Tau-coefficients for accuracy assessment of classification of remote-sensed data. Photogrammetric Engineering and Remote Sensing 61, 435–439. Major, R.E., Christie, F.J., Gowing, G., Cassus, G., Reid, C.A.M., 2003. The effect of habitat configuration on arboreal insects in fragmented woodlands of south-eastern Australia. Biological Conservation 113, 35–48. Manly, B., McDonald, L., Thomas, D., 1993. Resource selection by animals: statistical design and analysis for field studies. Chapman and Hall, London. Marshall, P., Szikszai, T., LeMay, V., Kozak, A., 1995. Testing the distributional assumptions of least squares linear regression. Forest Chronicle 71, 213–218. Martı´nez-Salvador, M., Valdez-Cepeda, R., Arias, H.R., Beltra´n-Morales, L.F., Murillo-Amador, B., Troyo-Die´guez, E., Ortega-Rubio, A., 2005. Distribution and density of maguey plants in the arid Zacatecas Plateau, Mexico. Journal of Arid Environments 61, 525–534. McGlone, M.S., Mildenhall, D.C., Pole, M.S., 1996. History and paleoecology of New Zealand Nothofagus forests. In: Veblen, T.T., Hill, R.S., Read, J. (Eds.), The Ecology and Biogeography of Nothofagus Forests. Yale University Press, New Haven, pp. 83–130. Molloy, J., 1995. Bat (Peka Peka) Recovery Plan. Threatened Species Recovery Plan, series number 15 Threatened Species Unit. Department of Conservation, Wellington. Montgomery, D.C., Peck, E.A., 1982. Introduction to Linear Regression Analysis. Wiley, London. Morris, W.F., Doak, D.F., 2002. Quantitative conservation biology: Theory and practice of population viability analysis. Sinauer Associates, Sunderland. Norberg, U.M., 1987. Wing form and flight mode in bats. In: Fenton, M.B., Racey, P., Rayner, J.M.V. (Eds.), Recent Advances in the Study of Bats. Cambridge University Press, Cambridge, pp. 43–56. Nugent, G., Fraser, W., Sweetapple, P., 2001. Top down or bottom up? Comparing the impacts of introduced arboreal possums and ‘terrestrial’ ruminants on native forests in New Zealand. Biological Conservation 99, 65–79.

B I O L O G I C A L C O N S E RVAT I O N

O’Donnell C.F.J., 1999. The ecology, spatial organisation and conservation of long-tailed bats Chalinolobus tuberculatus. Ph.D. thesis, University of Otago, Dunedin, New Zealand. O’Donnell, C.F.J., 2000. Conservation status and causes of decline of the threatened New Zealand long-tailed bat Chalinolobus tuberculatus (Chiroptera: Vespertilionidae). Mammal Review 30, 89–106. O’Donnell, C.F.J., 2001a. Advances in New Zealand mammalogy 1990–2000: Long-tailed bat. Journal of the Royal Society of New Zealand 31, 43–57. O’Donnell, C.F.J., 2001b. Home range and use of space by Chalinolobus tuberculatus, a temperate rainforest bat from New Zealand. Journal of Zoology, London 253, 253–264. O’Donnell, C.F.J., Langton, S. 2003. Power to detect trends in abundance of long-tailed bats (Chalinolobus tuberculatus) using counts on transects. Science for Conservation 224. Department of Conservation, Wellington. O’Donnell, C.F.J., Sedgeley, J.A., 1994. An automatic monitoring system for recording bat activity. Department of Conservation Technical Series 5. Department of Conservation, Wellington. O’Donnell, C.F.J., Sedgeley, J.A., 2001. Guidelines for surveying and monitoring long-tailed bat populations using line transects. Department of Conservation Science Internal Series 12. Department of Conservation, Wellington. Osborne, P.E., Tigar, B.J., 1992. Interpreting bird atlas data using logistic models: an example from Lesotho, Southern Africa. Journal of Applied Ecology 29, 55–62. Pearce, J., Ferrier, S., 2000. Evaluating the predictive performance of habitat models developed using logistic regression. Ecological Modelling 133, 225–245.

1 3 2 ( 2 0 0 6 ) 2 1 1 –2 2 1

221

Powell, M., Accad, A., Shapcott, A., 2005. Geographic Information System (GIS) predictions of past, present habitat distribution and areas for re-introduction of the endangered subtropical rainforest shrub Triunia robusta (Proteaceae) from south-east Queensland, Australia. Biological Conservation 123, 165–175. Rushton, S.P., Ormerod, S.J., Kerby, G., 2004. New paradigms for modelling species distributions? Journal of Applied Ecology 41, 193–200. Schaffer, W.M., 1981. On reproductive value and fitness. Ecology 62, 1683–1685. Seoane, J., Bustamante, J., Diaz-Delgado, J., 2004. Are existing vegetation maps adequate to predict bird distributions? Ecological Modelling 175, 137–149. Towns, D.R., Ferreira, S.M., 2001. Conservation of New Zealand lizards (Lacertilia: Scincidae) by translocation of small populations. Biological Conservation 98, 211–222. Towns, D.R., Daugherty, C.H., Cree, A., 2001. Raising the prospects for a forgotten fauna: a review of 10 years of conservation effort for New Zealand reptiles. Biological Conservation 99, 3–16. Wang, G.H., Owen, R.D., Sanches-Hernandez, C., RomeroAlmaraz, M.D.L., 2003. Ecological characterization of bat species distributions in Michoacan, Mexico, using a geographic information system. Global Ecology & Biogeography 12, 65–85. Webb, P.I., Sedgeley, J.A., O’Donnell, C.F.J., 1998. Wing shape in New Zealand lesser short-tailed bats (Mystacina tuberculata). The Journal of Zoology, London 246, 462–465. Yee, T.W., Mitchell, N.D., 1991. Generalized additive models in plant ecology. Journal of Vegetation Science 2, 587–602.