INTRODUCTION
The conversion of indigenous species habitat for human land uses is one of the leading causes of global biodiversity loss (Houghton Reference Houghton1994; Ricketts & Imhoff Reference Ricketts and Imhoff2003; Hoekstra et al. Reference Hoekstra, Boucher, Ricketts and Roberts2005; Reidsma et al. Reference Reidsma, Tekelenburg, van den Berg and Alkemade2006), and continues worldwide despite conservation efforts (Sala et al. Reference Sala, Chapin, Armesto, Berlow, Bloomfield, Dirzo, Huber-Sanwald, Huenneke, Jackson, Kinzig, Leemans, Lodge, Mooney, Oesterheld, Poff, Sykes, Walker, Walker and Wall2000; Mottet et al. Reference Mottet, Ladet, Coqué and Gibon2006; Brown et al. 2007; Kangalawe Reference Kangalawe2010). Increasing global population and greater demand for food, fodder, fibre and fuel is leading to rapid changes in land use patterns, and areas once considered impervious to human activity are increasingly coming under threat (Parks Reference Parks1995; de Koning et al. Reference de Koning, Verburg, Veldkamp and Fresco1999; Rouget et al. Reference Rouget, Cowling, Pressey and Richardson2003; Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005a).
Understanding the factors that drive habitat conversion and the spatial pattern of conversion is essential for the informed prioritization of conservation actions to minimize future biodiversity loss (Margules & Pressey Reference Margules and Pressey2000; Wilson et al. Reference Wilson, Underwood, Morrison, Klausmeyer, Murdoch, Reyers, Wardell-Johnson, Marquet, Rundel, McBride, Pressey, Bode, Hoekstra, Andelman, Looker, Rondinini, Kareiva, Shaw and Possingham2007). Habitat conversion is a major component of vulnerability (namely the likelihood or imminence of biodiversity loss to current or impending threatening processes; Pressey & Cowling Reference Pressey and Cowling2001; Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005b). Spatial statistical models provide insights into the important drivers of vulnerability and tools for spatial and temporal prediction of habitat conversion (Hall et al. Reference Hall, Tian, Qi, Pontius and Cornell1995; Pontius et al. Reference Pontius, Cornell and Hall2001). However, the distribution and rate of habitat conversion will change with time, for example through the exhaustion of formerly suitable areas or changes in global markets, technology and crops (Wilson Reference Wilson, Pressey, Newton, Burgman, Possingham and Weston2005b; Pressey et al. Reference Pressey, Cabeza, Watts, Cowling and Wilson2007). Therefore, predictions of future vulnerability based on spatial models of past habitat conversion will become less reliable over time, risking misallocation of scarce conservation resources (Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005b).
Because habitat conversion is a dynamic threat, it is important that practitioners keep abreast of conversion and regularly validate the utility of their vulnerability assumptions and models (Pressey et al. Reference Pressey, Cabeza, Watts, Cowling and Wilson2007). Validation requires testing the predictions of independent data (not those used in model parameterization) to ensure that the relationships inferred by a model are robust and the predictions reliable. Yet absence of validation is a common weakness of habitat-conversion models (Pontius et al. Reference Pontius, Huffaker and Denman2004), and the robustness of vulnerability models used in conservation planning is seldom assessed (Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005a; Pressey et al. Reference Pressey, Cabeza, Watts, Cowling and Wilson2007).
Here, we use statistical modelling techniques on observed land use conversion in New Zealand to understand drivers of conversion and how these change through time, and then test the power of the models to make spatial and temporal predictions of vulnerability. Using observed habitat conversion in indigenous grasslands in New Zealand over two recent time periods, we modelled the observed conversion against a range of potential environmental and socioeconomic explanatory variables to quantify the changing drivers and patterns of vulnerability to habitat conversion and make spatial predictions of vulnerability. We modelled and validated over three different time periods to provide temporal validation to test the ability of the models to predict future vulnerability. These results may have important implications for international conservation practitioners attempting to predict potential vulnerability to future habitat conversion.
METHODS
Study area
The study area covers c. 4.3 million hectares between latitudes 41° S and 46° S in the centre of New Zealand's South Island, east of the Southern Alps (Fig. 1). Elevation (height above sea level) ranges from sea level to 2750 m; rainfall ranges from 300 mm yr−1 at lower elevations to 7000 mm yr−1 at higher elevations. The area includes the majority (3.3 million hectares) of New Zealand's remaining indigenous ‘tussock’ grasslands (Levy Reference Levy1951; Mark Reference Mark1965, Reference Mark1993; Ashdown & Lucas Reference Ashdown and Lucas1987), including those largely induced by human fire below the natural treeline (Stevens et al. Reference Stevens, McGlone and McCulloch1988; McGlone Reference McGlone2001; Ewers et al. Reference Ewers, Kliskey, Walker, Rutledge, Harding and Didham2006) (Fig. 2).
The locations of the 4.3 million hectare study area in the interior of New Zealand's South Island. These areas encompass the largest extent of indigenous grasslands remaining in the region since 1990.
Land cover in New Zealand from pre-human (Maori) settlement (c. 800 years ago) (left) (McGlone Reference McGlone2001), and post-Maori and pre-European settlement (c. 1840) (middle) (McGlone Reference McGlone2001; Mark & McLennan Reference Mark and McLennan2005), to post-European settlement (c. 1990) (right) (Newsome Reference Newsome1987).
Within the study area, there is anecdotal evidence that conversion of indigenous grassland for pasture, forestry and urban uses is proceeding rapidly at lower elevations and on private land (Walker et al. Reference Walker, Price and Stephens2008; Mark et al. Reference Mark, Michel, Dickinson and McLennan2009), but inaccurate and outdated land use and land cover data have previously limited quantification of grassland vulnerability for conservation decision-making (Walker et al. Reference Walker, Price, Rutledge, Stephens and Lee2006). Weeks et al. (Reference Weeks, Walker, Dymond, Shepherd and Clarkson2013) showed that these data overestimated the extent of remaining grassland and underestimated recent rates of conversion. Although severely modified and invaded by exotic species (Mark & McLennan Reference Mark and McLennan2005; Mark et al. Reference Mark, Michel, Dickinson and McLennan2009) lower-elevation grasslands retain important residual indigenous biodiversity, including many threatened plants (de Lange et al. Reference de Lange, Norton, Courtney, Heenan, Barkla, Cameron, Hitchmough and Townsend2009; Walker et al. Reference Walker, Price and Stephens2008), and endemic lizards and invertebrates (Patterson Reference Patterson1992; Patrick & Dugdale Reference Patrick and Dugdale2000), and are poorly protected.
Data
Response variables
Three maps of grassland conversion were available to serve as our primary data for conversion from 1840 to 1990 (pre-1990), from 1990 to 2001, and from 2001 to 2008. Past (pre-1990) conversion was mapped by comparing the 1990 extent of remaining indigenous grasslands with the expert-based estimate of Mark and McLennan (Reference Mark and McLennan2005) of the extent of indigenous grasslands in 1840.
Recent conversion was manually digitized using three nominal dates (1990, 2001 and 2008) of Landsat TM, Landsat 7 ETM+ and SPOT-5 imagery satellite imagery (Weeks et al. Reference Weeks, Walker, Dymond, Shepherd and Clarkson2013). These maps were used to generate three different binary response variables (‘converted’ or ‘not converted’) from 1840 to 1990 (‘past’), from 1990 to 2001 (‘recent’) and from 1990 to 2008 (‘recent and current’), and one continuous response variable (‘area of conversion’ from 1990 to 2008). Grassland conversion was observed within the study area (Fig. 3). Each 25-m grid cell in a map represented either grassland that had not been converted or grassland that had been converted. We considered conversion to include exotic forestry and forestry tree weed spread, agriculture (cropland or exotic pasture) and urban development.
Observed conversion from 1840 to 2008 of indigenous grasslands to non-indigenous land uses within the study area. The insert shows detail for an expanded area in the centre of the region.
Explanatory variables
As potential predictor variables, we developed a comprehensive set of the environmental and socioeconomic variables that were potentially important drivers of conversion (Table 1). As environmental predictors, we used 11 climate, substrate and landform variables available from the Land Environments New Zealand (LENZ) database (Leathwick et al. Reference Leathwick, Morgan, Wilson, Rutledge, McLeod and Johnston2003). The 11 relatively independent environmental variables were selected from a set of 20 potential environmental variables by removal of strongly correlated potential predictors (correlation coefficient r < 0.75). Our 11 socioeconomic predictor variables represented aspects of governance, land tenure, infrastructure or productivity that had plausible explanatory links to conversion rates. Governance variables included local government administrative regions and water catchments derived from topographic maps. We also derived seven categories of land tenure from maps supplied by the New Zealand Department of Conservation (DOC). Four categories represent different designations on former Crown pastoral lease land that completed land reform by 2008 (conservation land, conservation land with a grazing licence, unencumbered privatized land, and private land encumbered by a conservation easement, colloquially a ‘covenant’). The remaining three categories were remaining Crown pastoral leases, public conservation land and private land.
Explanatory (predictor) variables used for the models, brief definitions, units, categories described in the text, and data sources. In the text, climate, substrate and landform categories are collectively referred to as environmental predictors, and governance, land tenure, infrastructure and productivity variables as socioeconomic predictors.
Six ‘infrastructure’ variables were derived from national digital topographic databases; we calculated 25-m raster layers of the distance to each infrastructure type. Because the influence of infrastructure is more important at short distances, and the majority of the landscape is at large distances, we log-transformed the distance variables to provide greater resolution of the effects at the shorter distances. A seventh variable (proximity of existing intensive agricultural activity) was calculated for the entire study area using neighbourhood statistics. For each pixel in the study area, the proximity of existing intensive agricultural activity was defined as the proportion of pixels within a 2-km radius that had been converted before 1990.
We used two different variables to represent agricultural productivity. Land use capability (LUC) was retrieved from the New Zealand Land Resource Inventory data layers (Lynn et al. Reference Lynn, Manderson, Page, Harmsworth, Eyeles, Douglas, Mackay and Newsome2009) and classifies land into eight categories according to their capability to sustain continuous production (Newsome et al. Reference Newsome, Wilde and Willoughby2000). We also used the continuous pasture productivity index created by Baisden (Reference Baisden2006), with values ranging from 41 (low productivity) to 2038 (high productivity).
The models
We created models predicting the ‘probability of conversion’ over three time periods, and a fourth model predicting the area of individual conversion events. The ‘past’ conversion model (1840–1990) and the ‘recent’ conversion model (1990–2001) estimated the probability of conversion over their respective time periods (150 years and 11 years), and were then used for temporal validation purposes by evaluating their ability to predict conversion in later time periods. The ‘past’ and ‘recent’ conversion models were also used to provide insight into the drivers of conversion in their respective time periods, and how the relative importance of these drivers changed with time. The ‘recent and current’ conversion model estimated the combined conversion over the period 1990–2008, and was used to predict the probability of future (post-2008) patterns of conversion. Unlike the ‘past’ and ‘recent’ models, the ‘recent and current’ model had no temporal boundaries. The ‘area of conversion’ model was used to investigate the drivers of the size (area) of each conversion event recorded between 1990 and 2008. Area of conversion is predicted for an 18 year time period.
Sampling design
Point data with response and predictor variables used to construct each model were generated by interrogating the GIS maps of response and predictor variables. The locations of the points were different for each model, and were determined by a stratified random sample of the study area. For each of the ‘probability of conversion’ models, 1000 points were sampled in ‘converted’ (for that time period) and 4000 points sampled in ‘not converted’ grassland. For the ‘area of conversion’ model, each of the 353 observed conversion polygons was treated as a single observation, and the area of the polygon was recorded at each random point.
Model calibration and identification of important explanatory variables
We used generalized additive models (GAMs) to model the probabilities of conversion and area in relation to the explanatory variables. As a modelling technique, GAMs provide good predictive power with relatively simple models that are easily understood. While their predictive performance is slightly less than machine learning techniques such as boosted regression trees (for example see De'ath Reference De'ath2007), the GAM models are easier to interpret and better at revealing relationships. We used generalized regression analysis and spatial prediction (GRASP), a set of functions (Lehmann et al. Reference Lehmann, Overton and Austin2002) in S-PLUS software (MathSoft 1998–1999) that facilitates the fitting, validation and cross-validation of GAM models and has been widely used for ecological modelling (see for example Guisan et al. Reference Guisan, Broennimann, Engler, Vust, Yoccoz, Lehmann and Zimmerman2006; Elith et al. Reference Elith, Leathwick and Hastie2008). We used a logistic link function for the models of the ‘probability of conversion’ (a categorical variable in which the observed outcome has the two possible values ‘converted’ and ‘not converted) and a Poisson link function for the ‘area of conversion’ model where the response variable (potential area of individual conversion events) varies continuously. In all models, a starting model including all continuous and categorical predictors was fitted first and significant predictor variables selected thereafter by backward and forwards stepwise procedure using the Bayesian information criterion (BIC) for variable selection. The chosen predictors for each final model were used to map predictions in geographic space.
To interpret the GAMs, we used plots of the partial response curves that resulted from the model, and the overall contribution of the variables to the model. Partial response curves allowed visualization of how the response variable varies as a function of the predictor variables, while the contributions allowed us to assess the relative importance of the predictor variables in explaining the variation in the response variable.
The models were exported to lookup tables and then mapped as raster surfaces in a GIS. Cell values in probability of conversion maps show likelihood that the cell was converted independent of other cells, while cell values in the area of conversion map show predicted area of conversion (in hectares) in that cell. To adjust for the bias in our sampling procedure, we rescaled the probability predictions (multiplying by 
$\frac{{\overline x oberved}}{{\overline x predicted}}$) so that the mean predicted probability was similar to the mean observed probability of conversion for the study area.
Temporal validation and spatial cross-validation
In the ‘past' and the ‘recent' conversion models, we used observed conversion from one time period to calibrate the model and observed patterns of conversion from 2001 to 2008 to temporally validate the model. We modelled conversion from 1840–1990 (‘past’) and from 1990–2001 (‘recent’) and used a temporal validation process to test the ability of these models to predict conversion from 2001–2008 (‘current’). For this temporal validation process we used receiver operating characteristic (ROC) curves (Swets Reference Swets1988; Pontius & Schneider Reference Pontius and Schneider2001; Pontius & Batchu Reference Pontius and Batchu2003) to evaluate model discrimination between ‘converted’ and ‘not converted’. The true positive rate (sensitivity) was plotted as a function of the false positives (1 – specificity) for threshold values of predicted probabilities to be considered ‘converted' or ‘not converted' ranging from 0 to 1. The area under the ROC curve measures the performance of the model. An area of 1 represents a perfect model, and a value over 0.9 is considered a very good model, while an area of 0.5 or less signifies a useless model.
Second, we assessed the accuracy of the temporal validation for each (i.e. pre-1990 and 1990–2001) model using interactive dot diagrams, that separate ‘not converted’ (0) and ‘converted (1) pixels on the horizontal axis and display probabilities of conversion from the models on the vertical axis. Thresholds identified in these dot diagrams indicate the best separation (minimal false negatives and false positives) between ‘not converted’ and ‘converted’ (Schoonjans et al. Reference Schoonjans, Zalata, Depuydt and Comhaire1995).
Our ‘recent and current' model of probability of conversion from 1990 to 2008 and our ‘area of conversion’ from 1990 to 2008 model were spatially cross-validated using k-fold cross-validation with five groups, ROC curves were again used to evaluate model performance for the probability of conversion model, and the correlation of predicted values with observed values was used for the area model (Lehmann et al. Reference Lehmann, Overton and Austin2002).
Vulnerability comparison
Finally, we also compared predicted vulnerability based on the pre-1990 and 1990–2008 conversions. Each map of vulnerability was scaled to span a range of 0 to 1 (by dividing by maximum probability), and the difference between the two maps calculated and mapped (scaled vulnerability based on 1990–2008 data minus scaled vulnerability based on pre-1990 data). A high positive difference would indicate that pre-1990 vulnerability underestimated future vulnerability, while more negative values indicate overestimates. We sampled the difference map at 5000 random points across the study area, and modelled difference in relation to our suite of predictor variables using a GAM (Table 1).
RESULTS
Model performance
ROC values of 0.913, 0.916 and 0.921 for our pre-1990, 1990–2001 and 1990–2008 conversion probability models, respectively, indicate very good predictive performance. The correlation between observed and predicted conversion was highest in the 1990–2001 model (Spearman correlation coefficient r = 0.708), followed by the 1990–2008 model (r = 0.693), and the pre-1990 model (r = 0.666). Correlation between observed and actual area of conversion was low (r = 0.467), suggesting this model explained less than a quarter of observed variation in conversion area (R 2 = 0.22).
Important explanatory variables
The ‘alone’ contributions of variables (the potential for each variable alone to explain conversion) differed among our three ‘probability of conversion’ models (Table 2). Mean annual temperature played a dominant role in the pre-1990 (‘past’) model. It was the best ‘alone’ predictor of conversion, followed by slope, rainfall, catchment group, vapour pressure deficit and annual water deficit. When dropping each predictor from the final model, mean annual temperature was the only variable whose contribution could not be compensated for by any of the other variables. Partial response curves for selected predictors in each model show that pre-1990 conversion was positively related to mean annual temperature and negatively related to slope, soil moisture deficit and vapour pressure deficit, and to rainfall < 1000 mm yr−1 (Appendix 1, Fig. S1, see supplementary material at Journals.cambridge.org/ENC).
The contribution (rounded to the nearest whole number) of selected predictors in models of ‘past’ (pre-1990), ‘recent’ (1990–2001), ‘recent and current’ (1990–2008), and area response variables. Drop contribution indicates the marginal contribution of each variable and is obtained by dropping each explanatory variable and calculating the associated change in deviance. Alone contributions reflect the potential of each variable and are calculated by creating new models with only one predictor. Difference is the scaled predicted probability of conversion in our model of ‘recent and current’ (1990–2008) conversion minus that in our model of ‘past’ (pre-1990) conversion. Note that values for the area model and the difference model cannot be compared with the values of the other models, because these models have different model families and sample sizes.
A wider selection of variables explained significant variation in conversion after 1990 (in the 1990–2001 and 1990–2008 models), and the variables were differently ranked in their ability to explain observed variation in conversion (Table 2, ‘alone' contributions). After 1990, elevation became a significant predictor of conversion in addition to temperature, slope and water balance, and significant differentiation between types of land tenure and effects of distance to roads and proximity to existing agricultural development also became apparent.
Slope was the dominant variable in both the 1990–2001 (‘recent’) and 1990–2008 (‘recent’ and ‘current’) probability of conversion models. Between 1990 and 2001, slope was the best ‘alone’ predictor of conversion, followed by elevation, distance to roads, rain, ratio of rainfall to potential evapotranspiration, land tenure, mean annual temperature, annual water deficit and proximity to existing agricultural development (Table 2). In the 1990–2008 model, the ranking was first slope, then rainfall, land tenure, distance to roads, proximity to existing agriculture, administrative region and mean annual temperature.
Conversion from 1990 to 2001 (‘recent’ conversion) was positively related to mean annual temperature and elevation, and negatively related to slope, rainfall, soil moisture deficit and distance to roads (Appendix 1, Fig. S1a, see supplementary material at Journals.cambridge.org/ENC). Probability of conversion peaked at intermediate proximity to existing agriculture (Appendix 1, Fig. S1b, see supplementary material at Journals.cambridge.org/ENC). Among land tenure categories, the probability of grassland conversion between 1990 and 2001 was highest on private land (category 6; Appendix 1, Fig. S1b, see supplementary material at Journals.cambridge.org/ENC), followed by former Crown pastoral lease land privatized by 2008 (category 3) and then by existing Crown pastoral lease land (category 7). A small portion of conversion was also predicted on former Crown pastoral lease land that was privatized and covenanted between 1992 and 2008 (category 4; Appendix 1, Fig. S1b, see supplementary material at Journals.cambridge.org/ENC), reflecting areas where the data recorded a ski field had been developed and exotic forestry tree weeds had spread. Probability of grassland conversion on recent conservation land (categories 1 and 2) and pre-existing public conservation land (category 6) was negligible.
Probability of conversion in the 18 years from 1990 to 2008 (‘recent and current’ conversion) was negatively related to slope, rainfall and distance to roads, and positively related to mean annual temperature; the model showed a similar peak at intermediate proximity to existing agricultural development as the 1990– 2001 model (Appendix 1, Fig. S1c, see supplementary material at Journals.cambridge.org/ENC). Probability of conversion was also higher within two regions (Canterbury and Otago). As with the 1990–2001 model, conversion probability was highest on private land from 1990 to 2008, followed by recently privatized former Crown pastoral leased land and then existing Crown pastoral lease land, and was negligible on public conservation land.
Proximity to existing agricultural development was the highest ranked predictor in our model of area of conversion from 1990 to 2008 (Table 2, ‘alone’ contribution), followed by slope, land tenure, October vapour pressure deficit, annual rainfall and pasture productivity. Larger areas of grassland were converted in places further from existing agricultural development (Appendix 1, Fig. S1d, see supplementary material at Journals.cambridge.org/ENC), and on flat and highly productive land (although large areas were also converted on some steep land). Area of conversion was similar to the probability of conversion in its response to land tenure, with the smallest conversion areas on public conservation land (category 6) and the largest on private land (category 7). Conversion on both conservation and private land created by recent land reform (categories 1 and 3) also tended to cover relatively large areas.
Spatial and temporal validation
Temporal validations of the ‘past' (pre-1990) and ‘recent' (1990–2000) conversion models using the observed ‘current' (2001–2008) conversion indicated that both provide good to very good ability to predict current conversion (Fig. 4). The ‘past' conversion model achieved an area under the ROC curve of 0.84 in predicting ‘current' conversion. Predictably, the ‘recent' change model (using more recent information on change) achieved a higher ROC score of 0.913 in predicting ‘current' conversion, indicating excellent ability to predict future conversion. The ‘recent and current' (1990–2008) conversion model was validated with the data from the same time period and achieved a ROC score of 0.92, indicating a very good ability to provide spatial prediction.
Temporal and spatial validations of vulnerability models. Here, ROC curves are used to assess the ability of the models to predict grassland conversion. The ability of the ‘past' (pre-1990) and ‘recent' (1990–2000) vulnerability models to predict later conversion is tested by evaluating their ability to predict the observed ‘current’ (2001–2008) conversion. The ‘recent and current' (1990–2008) model is spatially validated against data in the same time period. The diagonal line represents a model providing no discrimination.
Interactive dot diagrams were used to assess the accuracy (temporal validation) of the two models (Appendix 1, Fig. S2, see supplementary material at Journals.cambridge.org/ENC), and similarly showed that the ‘recent' conversion model had higher sensitivity (93%) and specificity (80%) than the ‘past' conversion model (86% and 71%, respectively).
Predicted vulnerability
We modelled the ‘probability of conversion’ for 150 years (Fig. 5a), 11 years (Fig. 5b), and 18 years (Fig. 5c), and area of conversion for 18 years (Fig. 5d). In all three ‘probability of conversion’ models, the predicted probability of conversion of indigenous grasslands that remain in 2008 was lowest in steep mountainous land and higher on inland basin floors and lower range slopes (Fig. 5a, b, c). However, there were differences in the amount of land predicted to be highly vulnerable to conversion, with the least land predicted to be vulnerable by the pre-1990 model (Fig. 5a), and the greatest area of land of high vulnerability to grassland conversion predicted by the 1990–2008 model (Fig. 5c).
Vulnerability of remaining indigenous grasslands to future conversion predicted by four different models. Probability of conversion was predicted based on observed patterns of conversion (a) pre-1990, (b) 1990–2001 and (c) 1990–2008. (d) Area (hectares) of conversion based on a model of observed conversion from 1990 to 2008. Areas of conversion before 1990 and areas that are not indigenous grasslands are labelled ‘non-indigenous’ (white). Legends are scaled so that spatial patterns of vulnerability can be directly compared.
The spatial distribution of sites predicted to be most vulnerable varied between the three ‘probability of conversion’ models, but most strikingly between the pre-1990 model and the two models based on conversion after 1990 (Figs. 5b, c). Comparison with the 1990–2008 model showed that predictions based on pre-1990 data underestimated indigenous grassland vulnerability on gentle slopes (< 10°) and at elevations of 500 m and 1100 m, and overestimated vulnerability at elevations above and below this range and on steeper slopes (see Table 2, difference, and Appendix 1, Fig. S1e, see supplementary material at Journals.cambridge.org/ENC). Sites predicted to be more vulnerable after 1990 also had lower mean annual temperatures, were in land-use capability classes with lower capability for cropping or pasture, were closer to roads, and were more likely to be within the Canterbury region.
The predicted area of conversion was 4–2147 ha (Fig. 5d). Large areas of conversion were predicted at lower elevations with low slopes, and small areas (4–116 ha) of conversion were predicted at higher elevations. Small-scale incremental conversion occurred near areas of existing agriculture activity, while large areas of conversion were further away from existing agriculture. In general, areas with the highest probability of conversion predicted by the 1990–2008 model of vulnerability were also the places predicted to have the largest areas of conversion. However, some places where large areas of conversion were predicted at higher elevations (Fig. 5d) had low probabilities of conversion (or low vulnerability; Fig. 5c). These places appear to represent the predicted locations of non-indigenous forestry plantations or forestry tree weed spread, rather than agricultural conversion.
DISCUSSION
Implications of validated vulnerability models
In this study, we demonstrate a rigorous empirical approach to quantifying and predicting vulnerability in space and time. We assessed the extent to which observations of past conversion can be used to predict future patterns of land conversion. This highlights the importance of using a robust validation process in estimating vulnerability (likelihood of future loss) in conservation planning.
A common weakness of land conversion modelling is the use of the same data for both calibration (making the model as consistent as possible with the data from which the parameters were estimated) and validation (assessment of the predictive power of the model) (Pontius et al. Reference Pontius, Huffaker and Denman2004). Lack of consideration of model uncertainty through rigorous validation has been shown to result in inaccurate and overconfident predictions (Pontius & Batchu Reference Pontius and Batchu2003). We calibrated our model using one set of data and compared the results to a future set, to determine how well the model used the general pattern in the calibration data to extrapolate a pattern of future conversion. Model validation is a critical aspect often overlooked in predictions of vulnerability in conservation planning (Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005b). Vulnerability assessments based on past observations assume that past patterns of conversion will remain the same in the future. However, it is well accepted that factors that are important for explaining patterns of conversion for a past period may not necessarily predict future landscape changes (Wilson et al. Reference Wilson, Newton, Echeverría, Weston and Burgman2005a). We tested the ability of models of past grassland conversion to predict conversion in a future time period using rigorous statistical validation procedures. In doing so, we demonstrated how rapidly patterns of vulnerability may change over time, using the indigenous grassland landscape of New Zealand as a test case.
By comparing patterns of conversion over time, we found that past conversion patterns provide a reasonably good guide to future conversion, yet conversion patterns also show a distinct trend. In New Zealand's indigenous grasslands, sites most vulnerable to conversion before 1990 were different from those most vulnerable between 1990 and 2001, and differed even further from those most vulnerable between 2001 and 2008. Clearly, more recent conversion data will provide the best estimates of future conversion. However, we caution that vulnerability estimates based on too narrow a time range may also provide less accurate forecasts, because they are based on an inadequate sample of conversion events. Our composite model of ‘recent and current’ conversion from 1990 to 2008 should therefore provide more accurate forecasts of grassland vulnerability to future conversion than models based solely on either the ‘recent’ (1990–2001) or ‘current’ (2001–2008) periods.
Past and present vulnerability to conversion in New Zealand's grasslands
While the distribution of more vulnerable New Zealand grasslands has changed considerably since European settlement, indigenous grasslands on steep slopes at high elevation have remained relatively safe from conversion. Historically (until 1990), grasslands on the lowest, warmest, flattest and driest land were converted most rapidly for production. Gentle slope and high mean annual temperature continue to be important predictors of grassland vulnerability, but the most vulnerable grasslands are now on higher (and therefore steeper, cooler and more marginal land for primary production) land than before 1990. Variations in infrastructure and governance have also emerged as important predictors of grassland conversion. This changing pattern of grassland vulnerability fits the ‘maximum power principle’, that people will use the most economically productive land first (Odum Reference Odum1983; Hall et al. Reference Hall, Tian, Qi, Pontius and Cornell1995. Our results indicate that the most economically viable land is becoming less available in New Zealand, while demand is rapidly increasing due to an increased demand for dairy exports (Ministry of Agriculture and Forestry 2007).
Our model predicting area of conversion was less accurate than our models of probability of conversion (vulnerability ‘exposure’). The correlation between observed and predicted areas corresponded to less than half (r2 = 0.47) of the variation in the area of conversion. Nevertheless, the area model also revealed that proximity to existing agricultural activity was an important predictor of the size of conversion events after 1990. The relatively poor predictive ability of our area model may be because environment is less important than socioeconomic drivers (such as land cost) that may be poorly represented in our models. Another complicating factor may be divergent predictors of large agricultural developments and large forestry plantations, and/or areas of forestry tree weed spread.
Ultimate and proximate threats in modelling vulnerability
Uncertainties caused by changes in ultimate threats underlie the limitations of models of vulnerability based on proximate threats (Pressey et al. Reference Pressey, Cabeza, Watts, Cowling and Wilson2007). Because land-use systems respond to a combination of proximate (biophysical) and ultimate (socioeconomic) drivers, modelling vulnerability to conversion (intensification) ideally requires a multidisciplinary approach (Veldkamp & Fresco Reference Veldkamp and Fresco1996; Lambin et al. Reference Lambin, Turner, Geist, Agbola, Angelsen, Bruce, Coomes, Dirzo, Fischer, Folke, George, Homewood, Imbernon, Leemans, Li, Moran, Mortimore, Ramakrishnan, Richards, Skånes, Steffen, Stone, Svedin, Veldkamp, Vogel and Xu2001; Rounsevell et al. Reference Rounsevell, Reginster, Araújo, Carter, Dendoncker, Ewert, House, Kankaanpää, Leemans, Metzger, Schmit, Smith and Tuck2006). Our models attempt to incorporate some dimensions of the socioeconomic environment, as well as the biophysical environment. However, incorporating socioeconomic data was challenging. Because these data are mainly held as aggregated national datasets in New Zealand, they do not represent patterns at regional and local levels. Furthermore, few are readily translated into spatial layers. We therefore used proxy variables such as administrative districts, distances to infrastructure and land tenure to represent dimensions of the socioeconomic environment such as governance, social contagion and law. These variables proved useful (and indeed necessary) for accurate predictions of conversion, especially after 1990, but did not allow us to explore causality (Veldkamp & Lambin Reference Veldkamp and Lambin2001). Therefore explanations for the correlations we observed need to be further explored, and models adopting different approaches to our biophysical models are likely to be needed. For example, models using system, actor-based and narrative approaches to incorporate endogenous variables (such as macroeconomic, land management technology, infrastructure and land use policy changes) often highlight the important role of economic opportunities and policies in driving land conversion (Lambin et al. Reference Lambin, Turner, Geist, Agbola, Angelsen, Bruce, Coomes, Dirzo, Fischer, Folke, George, Homewood, Imbernon, Leemans, Li, Moran, Mortimore, Ramakrishnan, Richards, Skånes, Steffen, Stone, Svedin, Veldkamp, Vogel and Xu2001).
Our incorporation of socioeconomic proxy variables into biophysical models did, however, provide some insights into potential higher-level drivers of grassland conversion. For example, we showed that grasslands on privately owned land, and on land that was formerly Crown pastoral lease but was privatized through land reform (‘tenure review’) between 1992 and 2008, had a high probability of conversion and is currently highly vulnerable. Current Crown pastoral lease land was also vulnerable, but less than private and recently privatized land. This is consistent with privatization through land reform increasing the vulnerability of remaining indigenous grassland habitats, as predicted by Walker et al. (Reference Walker, Price and Stephens2008). Our results also highlight that public conservation land appears to have provided complete protection against grassland conversion, except in a few cases where forestry tree weed spread was already advanced when the land was acquired for conservation through land reform. We observed an increase, after 1990, in intensive agricultural development on land previously considered suitable only for extensive grazing, which suggests a recent increase in the economic viability of irrigation on marginal land as a second potential higher-level driver. We also showed that the pattern of conversion varies regionally, which may be due to differences in the availability of developable land, historic development trends, variations in current administrative land use rules and policies, or some combination of these factors.
Although spatial regression models such as ours are useful for predicting the vulnerability of grasslands to conversion because of their robustness, we caution they do not account for temporal heterogeneity. Land-use decisions are often triggered by single events such as economic fluctuations or crises, often remote in space and time, which operate at a higher hierarchical level (Houghton Reference Houghton1994). Where temporal heterogeneity is high, process-based models or models using economic frameworks might be more appropriate and yield better representations of the decision-making process. However, for systematic conservation planning (Margules & Pressey Reference Margules and Pressey2000; Margules et al. Reference Margules, Pressey and Williams2002; Wilson et al. Reference Wilson, Underwood, Morrison, Klausmeyer, Murdoch, Reyers, Wardell-Johnson, Marquet, Rundel, McBride, Pressey, Bode, Hoekstra, Andelman, Looker, Rondinini, Kareiva, Shaw and Possingham2007) the two different types of model (process based and spatial regression) may complement each other. The strength of spatial regression models is in identifying effectively those areas vulnerable to threatening processes such as conversion, which could be targeted for increased protection, while process-based models foster understanding of the drivers of these changes, which could potentially be addressed through socioeconomic policy and instruments.
Incorporating vulnerability into conservation prioritization
Absence of robust vulnerability assessment such as this study provides may have hampered protection of New Zealand's more vulnerable indigenous grasslands. To anticipate threats, decision-makers have formerly relied on outdated land use and land cover data that overestimate the extent of remaining grassland and underestimate recent conversion rates (Weeks et al. Reference Weeks, Walker, Dymond, Shepherd and Clarkson2013). Furthermore, objective investigation of the correlates and potential drivers of grassland habitat loss for this zone have been lacking. We have also demonstrated that a vulnerability assessment such as ours may be relatively readily achieved with straightforward remote sensing and spatial regression techniques if the requisite satellite imagery and spatial predictor data are available.
However, a measure of vulnerability to conversion alone is insufficient for identifying priority areas for conservation, which must take into account not only vulnerability, but also measures of relative conservation value such as irreplaceability or significance. In applying the results of this study to conservation planning, it would also be useful to consider not only vulnerability to conversion, but also other types of vulnerability, (such as to invasive species and climate change) and the effectiveness and costs of different management approaches and activities (Carwardine et al. Reference Carwardine, Wilson, Hajkowicz, Smith, Klein, Watts and Possingham2009).
CONCLUSIONS
A variety of methods have been developed to model vulnerability to conversion. These models vary in their complexity and applicability to conservation management. Generalized additive models provide a simple robust method to explore predictors and patterns of land-use change. Our analysis demonstrates that they can also be feasibly used to predict future patterns of conversion using recent land conversion data.
Our results provide the first data-derived and statistically validated measurement of the vulnerability of New Zealand's indigenous grasslands to conversion, and show a trend to greater agricultural conversion on higher, more marginal land. They also show that models based on earlier conversion patterns performed more poorly in predicting modern conversion. Up-to-date land conversion data therefore appear crucial for accurately predicting future vulnerability to habitat conversion.
ACKNOWLEDGEMENTS
This paper was funded by The Miss E.L. Hellaby Indigenous Grasslands Research Trust (ESW) the New Zealand Foundation for Research, Science, and Technology (SW) and Ministry of Research, Science, and Technology capability funding to Landcare Research (JMO). Data were developed with support of Landcare Research. We particularly thank Robbie Price, Hamish Heke, and Janice Willoughby for their contributions to the processing and compiling of data, and Anne Austin for editing the manuscript.
