Using multivariate statistics to explore trade& ... - Wiley Online Library

Journal of Applied Ecology 2014, 51, 1504–1514

doi: 10.1111/1365-2664.12345

Using multivariate statistics to explore trade-offs among spatial planning scenarios Linda R. Harris1*, Matthew E. Watts2, Ronel Nel1, David S. Schoeman1,3 and Hugh P. Possingham2,4 1

Coastal and Marine Research Unit, Department of Zoology, Nelson Mandela Metropolitan University, PO Box 77000, Port Elizabeth 6031, South Africa; 2ARC Centre of Excellence for Environmental Decisions, School of Biological Sciences, The University of Queensland, Brisbane, Qld 4072, Australia; 3School of Science & Engineering, University of the Sunshine Coast, Maroochydore, Qld 4558, Australia; and 4Department of Life Sciences, Imperial College London, Silwood Park, Ascot SL5 7PY Berkshire, UK

Summary 1. Scenario planning can be useful to guide decision-making under uncertainty. While sys-

tematic conservation planning can create protected-area networks for multiple and complex reserve–design scenarios, planners rarely compare different reserve networks explicitly, or quantify trade-offs among scenarios. 2. We demonstrate the use of multivariate statistics traditionally applied in community ecology to compare reserves designed under different scenarios, using conservation planning for beaches in South Africa as an example. Twelve reserve–design scenarios were run in Marxan in a hierarchical experimental design with three levels: including/excluding the probability of site destruction; two different cost types; and three different configurations of existing terrestrial and marine reserves. 3. Multivariate statistics proved to be useful tools in the conservation planning context. In our case study, they showed that the trade-off associated with including the probability of site destruction during coastal reserve design depended on the cost type: if the cost is related to the site-destruction probability then reserves are significantly larger; if not, then reserves are significantly more costly. In both cases, the configuration of existing reserves locked a priori into the solutions was more important and resulted in significantly larger and more costly reserves. 4. Synthesis and applications. This study demonstrates a novel application of multivariate statistical tools to robustly quantify potential trade-offs among diverse sets of reserve–design scenarios. These statistics can be applied: to support negotiations with stakeholders and decision-makers regarding reserve configurations in the face of uncertainty; in reserve–design sensitivity analyses; and in priority setting for future research and data collection to improve conservation plans. Key-words: complete hierarchical cluster analysis, Marxan, non-metric multidimensional scaling, protected areas, reserve design, sandy beaches, scenario planning, site destruction, spatial prioritization, systematic conservation planning

Introduction Between cumulative impacts from human activities and stressors from global change (Halpern et al. 2008b; IPCC 2013), the intensity of environmental impacts is escalating at unprecedented rates (Butchart et al. 2010). Consequently, networks of protected areas and complementary ecosystem-based management programmes beyond reserves *Correspondence author. E-mail: [email protected]

are becoming increasingly important tools for conservation (e.g. Halpern et al. 2008a; Stokes et al. 2014). The placement and design of these biodiversity priority areas thus need to be carefully planned. We focus on reserve design for the remainder of the paper, but the principles and tools apply equally to other spatially explicit ecosystem-based management programmes. Constraints in reserve design are that space and resources are limited and that multiple stakeholders compete for exclusive rights to portions of land and

© 2014 The Authors. Journal of Applied Ecology © 2014 British Ecological Society

Multivariate statistics in spatial planning 1505 sea. Therefore, reserves must be efficient to reduce conflict with competing user groups, thereby increasing the likelihood of implementation. Systematic conservation planning (SCP; Margules & Pressey 2000; Moilanen, Wilson & Possingham 2009) is a well-established method for efficient reserve design that can meet the needs of multiple users (e.g. Watts et al. 2009). In this process, targets for representing key features are met, while the size and cost of reserves are minimized. One advantage of SCP is that it offers a variety of alternative configurations of protected-area networks rather than a single answer (e.g. Linke et al. 2011), allowing for flexibility and negotiation during decision-making. Apart from visual inspection, however, conservation planners have few tools to quantify trade-offs among different reserve– network designs (but see Weeks et al. 2010) and to be transparent and explicit about why some sites are selected more frequently for inclusion in reserve networks than others. Multivariate statistics are widely used to determine the patterns in complex data sets and to identify the corresponding explanatory variables. Such analyses have been of recent interest to conservation planners and are starting to be applied in reserve design and negotiation. For example, outputs of conservation planning exercises include a user-defined number of alternative protected-area networks (solutions) to a reserve–design problem (scenario). Studies have sought to distinguish distinct groups of similar solutions to a scenario, thereby providing alternative reserve configurations to decision-makers (Airame et al. 2003; Linke et al. 2011). However, because decisions are made under uncertainty, a scenario-planning approach is a prudent course of action (Peterson, Cumming & Carpenter 2003; Varum & Melo 2010). Consequently, ideally in consultation with key stakeholders, planners can construct multiple scenarios, for example with different feature targets, data inputs or trajectories of global change, that might address various conservation or management strategies for authorities to consider (e.g. Weeks et al. 2010; Levy & Ban 2013). Methods to compare solutions both within and among scenarios would therefore be useful, particularly if the underlying reasons for the different reserve configurations and their consequences/implications can be revealed simultaneously. To explore the use of multivariate statistics in comparing reserve networks under different scenarios, we use conservation planning for sandy shores in South Africa as a case study. This tripartite ecosystem comprises dunes, the intertidal and surf zone as a single geomorphic unit and requires land–sea reserves to be protected (Harris et al. 2014). Beaches support diverse assemblages of unique, largely endemic species (Harris et al. 2014) and provide key goods and services (Barbier et al. 2011). However, they are exposed to many ever-intensifying threats from human activities and global change (Defeo et al. 2009). Therefore, addressing their currently poor representation in coastal reserves is both critical and

urgent (Harris et al. 2014). Furthermore, coastal reserves for beach conservation must take resilience into account because sea-level rise poses a real threat to beaches where the upper shore is permanently constrained by development. In such cases, beaches lie trapped in a coastal squeeze, and the habitat, associated diversity and services will eventually be inundated and entirely lost (Dugan et al. 2008; Fish et al. 2008). Consequently, some beaches face a high probability of complete destruction in the near future; including these sites in coastal reserves would be inappropriate (see also Game et al. 2008). Reserve design for beaches thus requires scenario planning, serving here as a useful test case for our analyses. Overall, we aim to investigate the utility of multivariate statistics in quantifying potential trade-offs among different reserve–design scenarios. Note that although we use Marxan (Ball, Possingham & Watts 2009) to design the reserve networks, the methods should apply equally well to outputs derived from other conservation planning software, for example C-Plan (Pressey et al. 2009) or Zonation (Moilanen, Kujala & Leathwick 2009). In the case study, we aim to quantify the trade-offs of planning for resilience, that is including the probability of site destruction in reserve design. Specifically, we test whether including the probability of site destruction increases the size and/or cost of the solutions. These effects are considered in the context of using different proxies for cost and different conservation-status scenarios (i.e. locking in different combinations of existing protected areas). Given that including the probability of site destruction is a more constrained problem formulation than those excluding it, we therefore hypothesize that more resilient reserve networks (those that include the probability of site destruction) will need to comprise more planning units to meet the conservation targets and consequently will be larger and more costly. This study makes a contribution to the theory and practice of SCP by providing a set of statistical tools to compare reserve–design scenarios in a robust and defensible way.

Materials and methods CONSERVATION PLANNING IN MARXAN

Marxan solves reserve–design problems formulated as an algorithm. Briefly (for details see Ball, Possingham & Watts 2009), the study region is divided into planning units. Features (e.g. species distributions) and a cost metric (e.g. land-acquisition prices) are coded to these planning units and targets for each feature are set (e.g. proportion of a species’ distribution required to be included in protected areas). Marxan then generates a userdefined number of alternative reserve networks by selecting collections of planning units that cost effectively meet the features’ targets. Marxan may not meet the targets for certain features if they occur consistently in planning units with high cost values. To force those features’ targets to be met, a species penalty factor

© 2014 The Authors. Journal of Applied Ecology © 2014 British Ecological Society, Journal of Applied Ecology, 51, 1504–1514

Multivariate statistics in spatial planning 1507 The most severe threats are those that are both unstoppable and irreversible: sites exposed to such stressors have a very high probability of being destroyed, so their selection should be avoided during reserve design (see Game et al. 2008). In our example, the probability of site destruction should be informed by existing coastal development (stoppable, but irreversible) and the associated threat of coastal squeeze (unstoppable, irreversible); the remaining threats are stoppable and reversible (see Defeo et al. 2009; Harris 2012). The probability of site destruction was thus based on the intensity of coastal development and coastal squeeze per planning unit and an expert-based quantification of functional impact, scaled 0–1 (see Harris 2012). These were combined to give an overall probability of site destruction per planning unit (PUT), calculated as: PUT ¼ 1 ð1 PDEV Þð1 PCSQ Þ

eqn 1

where: PDEV = coastal development threat probability and PCSQ = coastal squeeze threat probability. The threat-avoidance target (ptarget1d) was set to 09 in all cases, that is, including sites in reserves where confidence that they will persist into the future is at least 90 %.

Level 2: cost Cost, in the context of SCP, relates to the relative demand on each site by competing sectors: the more interest or investment in a site by a competing sector(s), the greater the cost to secure that site for conservation. However, there is no standard cost metric; instead, the appropriate proxy depends on the context of the conservation plan (Naidoo et al. 2006; Ban & Klein 2009). Some studies used area to represent cost (e.g. Airame et al. 2003), although this has been criticized because not all areas are equally costly (Naidoo et al. 2006). Other studies explicitly incorporated economic values in planning, for example land-acquisition prices or management costs (e.g. Carwardine et al. 2010). But the cost of acquiring land or sea for conservation can extend beyond financial considerations by, for example, excluding certain activities from particular areas. In such cases, metrics for socioeconomic or opportunity costs may be used, for example, fisheries’ catch per unit effort (e.g. Lombard et al. 2007), or forgone revenues (e.g. Richardson et al. 2006). While these cost metrics can be expressed in financial terms, alternative metrics are required where such data are unavailable, for example the number of activities/people affected or displaced by formally conserving a site (e.g. Yates & Schoeman 2013). We lacked economic data, so considered two alternative cost metrics at Level 2: size of the coastal population (CCP) and stoppable/reversible threats (CRT). For coastal population, we assumed that the greater the number of people living adjacent to a site, the more the site is used and thus the greater its relative cost. The number of people living adjacent to the beach per planning unit (South African Census 2004 data; Statistics South Africa) was scaled 0–100 to set up cost metrics of directly comparable magnitude. All zero values were replaced with a negligible value (1 9 108) to avoid offering Marxan sites with zero cost. The stoppable/reversible threats to sandy beaches are all anthropogenic and are often associated (directly or indirectly) with recreational activities (see Defeo et al. 2009 for a review; and Harris 2012). Therefore, an integrated score representing all of these threats per planning unit is a good proxy for cost. To generate these cost data, the threatimpact scores (Ti; as calculated in Harris 2012) for all appropriate

threats were summed per planning unit and scaled 0–100 (see Appendix S2, Supporting information for the 15 threats included).

Level 3: conservation status To test if locking in different configurations of existing protected areas affected reserve design, three conservation-status scenarios were considered at Level 3: locking in all existing protected areas (both marine and terrestrial; CSAll); locking in only existing terrestrial protected areas (CSTerr); and not locking any protected areas into the solutions (CSNone). The data used to assign conservation status to the planning units included the following: formal and informal terrestrial reserves; priority areas in the National Protected Areas Expansion Strategy (Government of South Africa 2010); and marine protected areas (available at: http:// bgis.sanbi.org).

INPUT PARAMETERS AND EXECUTION OF SYSTEMATIC CONSERVATION PLANNING

All reserve networks were constructed in Marxan with Zones (version 2.01 64-bit Windows; Ball, Possingham & Watts 2009; Watts et al. 2009), executed in Zonae Cogito (Segan et al. 2011). Input data sets were compiled following standard procedures (Watts et al. 2008) and feature targets were taken from Harris et al. (2013). Prior to running each scenario, the following input parameters required calibration to ensure the scenarios met their targets without disproportionately increasing the cost: boundary length modifier (BLM); species penalty factor (spf) and threat probability weighting (PROBABILITYWEIGHTING). Given the complexity of the interacting, multi-parameter calibration, we simplified the problem by setting BLM to zero in all cases. This is because it would be difficult to justify that the level of clustering of planning units is equal across scenarios after adjusting the other two parameters, particularly for such a fragmented habitat. Virtually all targets (>90–95%) were met with spf = 1. Of those targets missed, only two were of key importance: African penguins and leatherback turtles. These two features had very high targets because of their respective endangered and critically endangered threat status (Harris et al. 2013). Meeting these targets was therefore deemed imperative, so the spf for these features was increased until the targets were met (spf = 100). The few other features that missed their targets were vegetation types. We were not concerned by this because vegetation extends well inland of the planning domain, and we were satisfied that these habitat types were only partially represented in the reserve network (generally, >70 % of the target was met). Finally, PROBABILITYWEIGHTING was set to 45; at this point, the threat-avoidance target was met within 1–6% (i.e. targets achieved at 94–99 %).

SCENARIO COMPARISON USING MULTIVARIATE STATISTICS

The statistical methods applied in this study are conventionally used to analyse patterns in community ecology, but are explored here as potential tools to analyse Marxan outputs. Effectively, each solution is treated analogously to a sample from a biological community, where ‘species’ (planning units) are either present (selected in the reserve network) or absent (not selected). Under this premise, each scenario represents a sampling site, and scenario–design inputs and solution results are used as explanatory


Multivariate statistics in spatial planning 1507 The most severe threats are those that are both unstoppable and irreversible: sites exposed to such stressors have a very high probability of being destroyed, so their selection should be avoided during reserve design (see Game et al. 2008). In our example, the probability of site destruction should be informed by existing coastal development (stoppable, but irreversible) and the associated threat of coastal squeeze (unstoppable, irreversible); the remaining threats are stoppable and reversible (see Defeo et al. 2009; Harris 2012). The probability of site destruction was thus based on the intensity of coastal development and coastal squeeze per planning unit and an expert-based quantification of functional impact, scaled 0–1 (see Harris 2012). These were combined to give an overall probability of site destruction per planning unit (PUT), calculated as: PUT ¼ 1 ð1 PDEV Þð1 PCSQ Þ

eqn 1

where: PDEV = coastal development threat probability and PCSQ = coastal squeeze threat probability. The threat-avoidance target (ptarget1d) was set to 09 in all cases, that is, including sites in reserves where confidence that they will persist into the future is at least 90 %.

Level 2: cost Cost, in the context of SCP, relates to the relative demand on each site by competing sectors: the more interest or investment in a site by a competing sector(s), the greater the cost to secure that site for conservation. However, there is no standard cost metric; instead, the appropriate proxy depends on the context of the conservation plan (Naidoo et al. 2006; Ban & Klein 2009). Some studies used area to represent cost (e.g. Airame et al. 2003), although this has been criticized because not all areas are equally costly (Naidoo et al. 2006). Other studies explicitly incorporated economic values in planning, for example land-acquisition prices or management costs (e.g. Carwardine et al. 2010). But the cost of acquiring land or sea for conservation can extend beyond financial considerations by, for example, excluding certain activities from particular areas. In such cases, metrics for socioeconomic or opportunity costs may be used, for example, fisheries’ catch per unit effort (e.g. Lombard et al. 2007), or forgone revenues (e.g. Richardson et al. 2006). While these cost metrics can be expressed in financial terms, alternative metrics are required where such data are unavailable, for example the number of activities/people affected or displaced by formally conserving a site (e.g. Yates & Schoeman 2013). We lacked economic data, so considered two alternative cost metrics at Level 2: size of the coastal population (CCP) and stoppable/reversible threats (CRT). For coastal population, we assumed that the greater the number of people living adjacent to a site, the more the site is used and thus the greater its relative cost. The number of people living adjacent to the beach per planning unit (South African Census 2004 data; Statistics South Africa) was scaled 0–100 to set up cost metrics of directly comparable magnitude. All zero values were replaced with a negligible value (1 9 108) to avoid offering Marxan sites with zero cost. The stoppable/reversible threats to sandy beaches are all anthropogenic and are often associated (directly or indirectly) with recreational activities (see Defeo et al. 2009 for a review; and Harris 2012). Therefore, an integrated score representing all of these threats per planning unit is a good proxy for cost. To generate these cost data, the threatimpact scores (Ti; as calculated in Harris 2012) for all appropriate

threats were summed per planning unit and scaled 0–100 (see Appendix S2, Supporting information for the 15 threats included).

Level 3: conservation status To test if locking in different configurations of existing protected areas affected reserve design, three conservation-status scenarios were considered at Level 3: locking in all existing protected areas (both marine and terrestrial; CSAll); locking in only existing terrestrial protected areas (CSTerr); and not locking any protected areas into the solutions (CSNone). The data used to assign conservation status to the planning units included the following: formal and informal terrestrial reserves; priority areas in the National Protected Areas Expansion Strategy (Government of South Africa 2010); and marine protected areas (available at: http:// bgis.sanbi.org).

INPUT PARAMETERS AND EXECUTION OF SYSTEMATIC CONSERVATION PLANNING

All reserve networks were constructed in Marxan with Zones (version 2.01 64-bit Windows; Ball, Possingham & Watts 2009; Watts et al. 2009), executed in Zonae Cogito (Segan et al. 2011). Input data sets were compiled following standard procedures (Watts et al. 2008) and feature targets were taken from Harris et al. (2013). Prior to running each scenario, the following input parameters required calibration to ensure the scenarios met their targets without disproportionately increasing the cost: boundary length modifier (BLM); species penalty factor (spf) and threat probability weighting (PROBABILITYWEIGHTING). Given the complexity of the interacting, multi-parameter calibration, we simplified the problem by setting BLM to zero in all cases. This is because it would be difficult to justify that the level of clustering of planning units is equal across scenarios after adjusting the other two parameters, particularly for such a fragmented habitat. Virtually all targets (>90–95%) were met with spf = 1. Of those targets missed, only two were of key importance: African penguins and leatherback turtles. These two features had very high targets because of their respective endangered and critically endangered threat status (Harris et al. 2013). Meeting these targets was therefore deemed imperative, so the spf for these features was increased until the targets were met (spf = 100). The few other features that missed their targets were vegetation types. We were not concerned by this because vegetation extends well inland of the planning domain, and we were satisfied that these habitat types were only partially represented in the reserve network (generally, >70 % of the target was met). Finally, PROBABILITYWEIGHTING was set to 45; at this point, the threat-avoidance target was met within 1–6% (i.e. targets achieved at 94–99 %).

SCENARIO COMPARISON USING MULTIVARIATE STATISTICS

The statistical methods applied in this study are conventionally used to analyse patterns in community ecology, but are explored here as potential tools to analyse Marxan outputs. Effectively, each solution is treated analogously to a sample from a biological community, where ‘species’ (planning units) are either present (selected in the reserve network) or absent (not selected). Under this premise, each scenario represents a sampling site, and scenario–design inputs and solution results are used as explanatory


1508 L. R. Harris et al. environmental variables. This allows computation of multivariate statistics based on resemblance matrices, which in turn facilitate determining which reserve–design parameters contribute most to the variability and differences among solutions. A complete hierarchical cluster analysis was performed in R version 2.15.0 (R Development Core Team 2012), using the hclust function (stats package), to compare the resemblance within and among scenarios, using all solutions as samples of each scenario. This hierarchical cluster analysis was based on a Jaccard resemblance matrix constructed using the vegdist function (vegan package; Oksanen et al. 2012). Linke et al. (2011) suggest that this resemblance matrix should be constructed using the Bray–Curtis method on the premise that this method excludes joint absences. However, the Jaccard method should be preferred, because it similarly ignores joint absences but is constructed explicitly for binary (presence/absence) data, which is the output format of Marxan solutions. A dendrogram of the complete hierarchical cluster analysis was constructed using the ColorDendrogram function (sparcl package; Witten & Tibshirani 2010, 2012). To further visualize the data, a non-metric multi-dimensional scaling (nMDS) ordination of the solutions per scenario was rendered using the metaMDS function (vegan package; Oksanen et al. 2012), based on the Jaccard resemblance matrix. An envfit analysis (vegan package; Oksanen et al. 2012) was run on a data set of metrics associated with the solutions to each scenario to determine which of the explanatory variables are best correlated across the nMDS ordination. These variables included the discrete factors: cost; conservation status; and threatprobability scenario – that is the scenario–design categories (Table 1); and the continuous vectors: number of planning units selected in the solution; total boundary length of the reserve network; total cost of the solution; and Marxan penalty value. The variables that were significantly correlated across the ordination surface (at a = 005) were plotted as centroids for factors and GAM spline isopleths for vectors on the nMDS surface. Colinearity among explanatory vectors was assessed using a pairs plot, with relationships illustrated in the form of LOESS smoothers and corresponding correlation coefficients. A three-factor analysis of variance (ANOVA; stats package) was used to determine whether there were significant differences among scenarios in terms of mean reserve size (number of planning units and reserve boundary length) and cost. Attempts were made to eliminate higher order interactions using likelihood ratio tests and the Akaike information criterion, but both main effects and all interactions were significant and models could not be simplified. Therefore, a univariate ANOVA was used to evaluate variables by a concatenated factor (i.e. comparing data among the 12 scenarios) and with planned contrasts fitted from the gmodels package (Warnes 2012) because analysing all pairwise comparisons and selecting only those of interest would increase the rate of Type I errors (see Ruxton & Beauchamp 2008). All our R code is available in Appendix S3 (Supporting information).

Results EXPLORATORY PLOTS OF MARXAN SOLUTIONS

Clustering of solutions and scenarios The complete hierarchical clustering (represented as a dendrogram) and ordination (nMDS biplot) of the

Marxan output solutions per scenario proved to be very useful tools for data exploration and visualization (Fig. 1). Scenarios 3 (pTICRTCSNone) and 9 (pTECRTCSNone), and 6 (pTICCPCSNone) and 12 (pTECCPCSNone) were most dissimilar to the other scenarios (Fig. 1). These scenarios do not have any existing reserves locked into their solutions (Table 1). The second major split divides solutions by their cost type. These child nodes are split further depending on their conservation status, that is whether all reserves or only terrestrial reserves were locked into the solutions. The final level of solution clustering depends on the cost type. Reserve networks designed with stoppable, reversible threats as the cost show no clear differences in configuration between scenarios that include or exclude site-destruction probabilities (no distinct clustering of the solutions into single bands of colour). In contrast, those designed with coastal population as the cost cluster neatly into two dichotomous groups; those that include

Fig. 1. Relationships among solutions for each of the 12 scenarios. Data are presented as (a) a dendrogram from a complete hierarchical cluster analysis and (b) an nMDS biplot based on a Jaccard resemblance matrix. S1 – S12 refer to the 12 scenarios (Table 1). Solutions are coloured to represent pairs of scenarios, with those that include site-destruction probabilities (S1 – S6) coloured in a darker shade of the same colour than their counterpart solutions that exclude site-destruction probabilities (S7 – S12) – refer to the article online for the coloured version of this figure.


Multivariate statistics in spatial planning 1509 site-destruction probabilities and those that exclude them. This is most distinct for Scenarios 6 and 12 (Fig. 1a). Note also that as the conservation status was relaxed (fewer reserves locked into the solutions), the variability among solutions increased, the range of dissimilarity values increases per scenario (Fig. 1a), and the spread of data points in the ordination is greater (Fig. 1b). This suggests that conservation status and cost type appear to drive greater differences in reserve configurations than does avoiding sites with high probability of site destruction. Significant explanatory factors and vectors: causes and effects of different reserve configurations The greater contribution of conservation status and cost to the dissimilarity of solutions among scenarios was confirmed by the envfit analysis (Fig. 2a). Of the factors, conservation status had the strongest correlation with the nMDS ordination structure (r2 = 0663; P = 0001), distantly followed by cost type (r2 = 0120; P = 0001); threat probability was not significantly correlated with the nMDS surface (r2 < 0001; P = 0996). All vectors from the solution summaries were significantly correlated with the nMDS ordination structure (P = 0001), three of which had high correlation coefficients (Fig. 2b–e): Marxan penalty score (r2 = 0927), number of planning units selected (r2 = 0833) and reserve boundary length (r2 = 0815). In contrast, cost (Fig. 2f) was only weakly correlated with the nMDS surface (r2 = 0183). However, there are strong correlations among some of these vectors (Fig. 3), suggesting that it is very likely that they reflect slight variations of the same effect on the nMDS surface, rather than independent trends. Specifically, the more planning units selected, the greater the reserve boundary length and the lower the Marxan penalty. The factors are used to design reserve networks, and thus they drive the trends across the ordination surface. The vectors arise as a result of reserve design and therefore represent the effects that the factors have on reserve networks. For example, conservation status drives a significant separation of data points along the first ordination axis (NMDS1), whereas cost types drive significant separation along the second ordination axis (NMDS2; Fig. 2a). The increase in the cost vector (Fig. 2b) and isopleths (Fig. 2f) diagonally towards the bottom left of the ordination suggests that solutions with reserves locked in and with reversible threats as the cost type have the most costly reserve solutions. Performance of Marxan with site-destruction probabilities All variables tested were significantly different among scenarios (one-way ANOVA: P < 0001 in all cases), with mostly significant differences among scenario pairs in the planned contrasts test (Fig. 4). In confirmation of results above, the more the conservation status is relaxed, the

smaller and less costly reserve networks become. However, also as indicated above, the relative performance of scenarios including site-destruction probabilities differs depending on the cost type. When reversible threats are the cost type, reserve networks that include site-destruction probabilities are not necessarily larger than their counterpart scenarios that exclude site-destruction probabilities, but are significantly more costly. In other words, when site-destruction probabilities are included, Marxan selects a similar number of more costly sites to achieve the targets. In contrast, when coastal population is the cost type, reserve networks that include site-destruction probabilities during their design are significantly larger, but not more costly than their counterpart scenarios that exclude site-destruction probabilities. This means that when site-destruction probabilities are included, Marxan selects a greater number of less costly sites to achieve the targets.

Discussion MULTIVARIATE STATISTICS: TOOLS TO COMPARE SOLUTIONS AND QUANTIFY TRADE-OFFS

Meta models are increasingly being used to aid the understanding of complex outputs from ever more sophisticated modelling approaches (e.g. Coutts & Yokomizo 2014). Here, we demonstrated that multivariate statistics are easily adapted from their traditional use in community ecology to provide novel insights into scenario evaluation in conservation planning. Specifically, their application enabled us to determine which scenarios had the most flexible solutions (greatest variability among solutions), which scenarios resulted in statistically different reserve networks from others and why; and what the trade-offs among scenarios might be. These capabilities have important applications in reserve design and implementation. The simplest of these is in scenario planning. Planners can present alternatives to decision-makers and highlight what the likely trade-offs are depending on which solutions/scenarios they implement (see also Airame et al. 2003; Linke et al. 2011). This is especially useful in contemporary conservation planning where reserve networks are designed under uncertainty, can include complex data sets and design constraints, and account for a variety of features, processes and socio-economic information. Despite the complexity, it is still possible to provide simple and defensible explanations for why some sites are selected over others and what the implications of particular scenarios are. Another key application is in sensitivity analyses. By comparing scenarios where some of the design parameters are altered each time, for example using different targets or surrogates for biodiversity/costs, planners can identify which data sets or input parameters drive the configuration of reserve networks, which might highlight where key uncertainties are in the plans (see also Peterson, Cumming


Connectivity No_PUs

CostReversible

−0·4

−0·2

0·0

0·2

0·4

0·6

0·1

NMDS2

ConsStatusNone

ConsStatusAll ConsStatusTerrestrial

−0·3 −0·2 −0·1 0·0

0·1

CostPopulation

–0·3 −0·2 −0·1 0·0

NMDS2

(b)

0·2

0·3

(a)

0·2

0·3

1510 L. R. Harris et al.

Cost Score

Penalty

−0·4

−0·2

0·0

0·3

(c)

0·6

0·2

0·4

0·6

0·2

0·4

0·6

NMDS2 −0·2

0·0

0·2

0·4

0·6

0·1

0·2

(d)

−0·4

−0·2

0·0

NMDS1

NMDS1

0·0

0·2

0·4

0·6

NMDS1

0·1

NMDS2 −0·2

−0·3 −0·2 −0·1 0·0

0·1 −0·3 −0·2 −0·1 0·0 −0·4

(f)

0·2

0·3

(e)

0·2

0·3

0·4

−0·3 −0·2 −0·1 0·0

0·1

NMDS2

−0·3 −0·2 −0·1 0·0 −0·4

NMDS2

0·2

NMDS1

0·2

0·3

NMDS1

−0·4

−0·2

0·0

NMDS1

Fig. 2. A non-metric multi-dimensional scaling (nMDS) plot of the solutions based on a Jaccard resemblance matrix. (a) Factors and (b) vectors that are significantly correlated with the nMDS surface (at a = 005) are plotted. Factors are represented as text labels on their respective centroids. Arrows indicate the direction in which the vector increases, with the length of the arrow representing the relative importance of the vector (i.e. longer arrows have higher correlation coefficients). Isopleths of the vectors that are significantly correlated to the nMDS surface are plotted: (c) Marxan penalty score (r2 = 0927, P = 0001); (d) number of planning units selected (r2 = 0833, P = 0001); (e) reserve boundary length (km; r2 = 0815, P = 0001); and (f) cost (r2 = 0183, P = 0001).

& Carpenter 2003; Coutts & Yokomizo 2014). If the outputs are driven by the data sets of greatest uncertainty or lowest quality, our approach could challenge how these are weighted in the planning process and/or inform where

resources should be directed to obtain better data. Note that the answer to this kind of question may also depend on the conservation goal (Perhans et al. 2008). Similarly, our analytical approach could be useful in negotiations if


Multivariate statistics in spatial planning 1511 800

0 100 200 300 400

0·19

0·35

0·19

No_PUs

0·96

0·91

Connectivity

0·95

300

5

7

9

11

400 600 800

Cost

600

0 5000 15 000

400

0 100

Penalty

0 5000

15 000

5 6 7 8 9 10

Fig. 3. Correlations among Marxan output values. Data plotted as scatterplots with a LOESS smoother below the diagonal and corresponding Pearson’s correlation coefficients given above the diagonal. Cost = cost; No_PUs = number of planning units selected in the reserve network; Connectivity = reserve boundary length; Penalty = Marxan penalty score.

planners can show decision-makers and stakeholders that the most important biodiversity features, for example critically endangered vertebrates or rare threatened habitats, drive reserve configurations. Statistical comparison of scenarios also has theoretical applications that might guide future research priorities. For example, studies might determine a minimum resolution for input data required to design suitable reserve networks for a particular ecosystem. This might help to answer the question of whether designing reserve networks with habitats as a surrogate for biodiversity are significantly different from those designed with actual species’ distributions. The resolution of the cost data has similarly been shown to be important for designing reserves efficiently (Richardson et al. 2006) and some surrogates may fail to accurately represent socio-economic costs (Weeks et al. 2010). Alternatively, it might be worth exploring whether including data describing fine-scale biodiversity patterns, fine-scale costs, dynamic processes or global-change phenomena have a bigger influence on reserve configurations. This could have important implications for how limited research funds are spent in order to achieve the most robust and efficient conservation plans (see also Perhans et al. 2008). In short, regardless of the context within which the multivariate statistics are applied, they provide potentially useful new tools for conservation planners. Furthermore, assuming most practitioners are ecologists by training, uptake and use of these ecology-familiar statistics should come more easily than use of the optimization software itself.

RESILIENT RESERVES: TRADE-OFFS OF INCLUDING SITE-DESTRUCTION PROBABILITIES

Including site-destruction probabilities in the Marxan algorithm had the weakest effect on reserve–network solutions among the scenario–design parameters explored; conservation status had the greatest effect. This suggests that, in our case study, planning for resilience in the future has much less of an effect on reserve design than does the legacy of ad hoc proclamation of protected areas, where the latter are inefficient and more costly in the long term than if they had been designed a priori (e.g. Pressey 1994; but see Hansen et al. 2011). The explanation for the different trade-offs of including/excluding site-destruction probabilities between scenarios that use coastal population (larger reserves) and stoppable, reversible threats (more costly reserves) as the cost, relates to the relationship between the cost measure and probability of site destruction. Where coastal population is the cost, sites with high site-destruction probabilities (greater coastal development) tend to have higher costs (higher coastal population densities). Therefore, avoiding sites with a high probability of destruction also means avoiding sites of high cost. Consequently, reserves that account for the probability of site destruction (with coastal population as the cost metric) are larger but not more costly. This is not true when stoppable, reversible threats are the cost type, because sites that have low sitedestruction probabilities (undeveloped/remote sites) might still carry a high cost (because they support a number of


Coastal Population

t = 7·50, P < 0·001

Cost (arbitrary units)

10 000

t = 12·56, P < 0·001

400

500

t = 2·75, P = 0·006

t = 0·09, P = 0·928

S1

S7

S2

S8

S3

S9

S5

S11

S6

S12

10·5

(d)

10·0

t = 76·75, P < 0·001

9·5

t = 9·49, P < 0·001

t = 79·17, P < 0·001

S9

t = 3·00, P = 0·003

S10

S4

S10

S5

S11

S6

S12

(f)

850

S3

t = 1·63, P = 0·104

S4

750

S8

(e)

600

700

S2

t = 1·50, P = 0·135

9·0

6

t = −0·74, P = 0·457

S7

t = 1·48, P = 0·138

S9

Reserve boundary length (km x 103)

S3

Number of PUs in reserve network

S8

t = 3·37, P < 0·001 t = 1·76, P = 0·079

7

8

S2

(c)

S1

Number of PUs in reserve network

S7

5

Reserve boundary length (km x 103)

S1

(b)

t = 83·60, P < 0·001

t = 18·84, P < 0·001

t = 79·72, P < 0·001

650

14 000

t = 11·17, P < 0·001

100 150 200 250 300 350

Reversible Threats (a)

6000

Cost (arbitrary units)

18 000

1512 L. R. Harris et al.

S4

S10

S5

S11

S6

S12

Fig. 4. Paired boxplots showing the differences in solution (a, b) cost, (c, d) reserve boundary length and (e, f) number of planning units (PUs) selected in the reserve network, between scenarios that include site-destruction probabilities and those that exclude them. Scenarios are arranged in scenario pairs on the x-axis (see Table 1) and per cost type (reversible threats: a, c, e, g and coastal population: b, d, f, h) to facilitate visualization because of the differences in scale (y-axis). Data are presented as the median (thick black line), interquartile range (box; grey-filled for scenarios including site-destruction probabilities; white-filled for scenarios excluding site-destruction probabilities), 25th and 75th percentiles (whiskers) and outliers (dots). Results from the planned contrasts test plotted on the figure; significant differences indicated with an asterisk.

different activities, e.g. fishing, beach driving and subsistence harvesting). In this latter case, selecting sites with a low probability of destruction does not necessarily mean simultaneously avoiding sites of high cost, and hence the reserve networks including the probability of site destruction are more costly. The conclusion therefore seems to be this: if the cost is positively related to the site-destruction probability, even if not directly, then reserves will need to be larger but not more costly; if, conversely, cost is unrelated to the site-destruction probability, then

reserves will need to be more costly, but not necessarily larger. Whether this holds as a generalization or not will require further testing. IMPLICATIONS FOR POLICY AND MANAGEMENT

The specific implications for policy and management will relate largely to the local context and the details of the conservation plan. Generally, however, because scenario planning can guide formulation of conservation policies


Multivariate statistics in spatial planning 1513 that are more resilient in the face of uncertainty (Peterson, Cumming & Carpenter 2003), being able to compare scenarios statistically will add more rigour to the decisionmaking process. This is particularly important in the context of global change, where it is unclear which trajectory of socio-economic development and climate change is more likely (see IPCC 2013) and ecological systems can sometimes respond to change in drastic and unexpected ways (e.g. Christensen et al. 2006). In terms of the example of conservation planning for sandy shores, when human use of the shore (stoppable, reversible threats) is the cost, including the probability of site destruction results in reserve networks that are more costly. In other words, if we are to conserve beaches, the cost to society of developing the shoreline inappropriately (thus trapping beaches in a coastal squeeze) is that fewer activities will be permissible on sandy shores. The converse opportunity, however, is that greater use of the shore generally is permissible if coastal development proceeds in an ecologically sensible way. In coastal reserves for beaches specifically, management actions should be aimed at eliminating, or at the very least ameliorating the impacts of reversible threats, particularly where overlapping activities have negative and synergistic impacts. Perhaps most importantly, though, managers should strive to preclude the presence of other stoppable threats in the reserves, particularly coastal development, because it has irreversible impacts. CONCLUSION

Multivariate statistics can be used to compare different reserve–design scenarios and quantify the potential tradeoffs among them. These analyses can be adapted to compare sets of solutions for a variety of scenarios produced by any conservation planning software that gives multiple solutions. The potential application of these multivariate statistics are vast and include the following: supporting negotiations during decision-making and reserve implementation in the face of uncertainty; analysing reserve– design sensitivity; and guiding research prioritization and future collection of data.

Acknowledgements We are grateful to three anonymous reviewers who helped to improve this paper substantially. Financial support for this research (for LRH) was provided by the National Research Foundation (NRF) and the Nelson Mandela Metropolitan University (Research Capacity Development and Science Faculty). Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF. The authors declare no conflict of interest. The authors also acknowledged funding from the Australian Research Council Centre of Excellence for Environmental Decisions, National Environmental Research Program Environmental Decisions Hub (M.E.W., H.P.P.) and Australian Government through the Collaborative Research Network (D.S.S.).

References Airame, S., Dugan, J.E., Lafferty, K.D., Leslie, H., McArdle, D.A. & Warner, R.R. (2003) Applying ecological criteria to marine reserve

design: a case study from the California Channel Islands. Ecological Applications, 13, S170–S184. Ball, I., Possingham, H. & Watts, M. (2009) Marxan and relatives: software for spatial conservation prioritization. Spatial Conservation Prioritization (eds A. Moilanen, K. Wilson & H. Possingham), pp. 185–195. Oxford University Press, Oxford. Ban, N.C. & Klein, C.J. (2009) Spatial socioeconomic data as a cost in systematic marine conservation planning. Conservation Letters, 2, 206–215. Barbier, E.B., Hacker, S.D., Kennedy, C., Koch, E.W., Stier, A.C. & Silliman, B.R. (2011) The value of estuarine and coastal ecosystem services. Ecological Monographs, 81, 169–193. Butchart, S.H.M., Walpole, M., Collen, B., van Strien, A., Scharlemann, J.P.W., Almond, R.E.A. et al. (2010) Global biodiversity: indicators of recent declines. Science, 328, 1164–1168. Campbell, E.E. & Bate, G.C. (1997) Coastal features associated with diatom disoloration of surf-zones. Botanica Marina, 40, 179–185. Carwardine, J., Wilson, K.A., Hajkowicz, S.A., Smith, R.J., Klein, C.J., Watts, M. & Possingham, H.P. (2010) Conservation planning when costs are uncertain. Conservation Biology, 24, 1529–1537. Christensen, M.R., Graham, M.D., Vinebrooke, R.D., Findlay, D.L., Paterson, M.J. & Turner, M.A. (2006) Multiple anthropogenic stressors cause ecological surprises in boreal lakes. Global Change Biology, 12, 2316–2322. Coutts, S.R. & Yokomizo, H. (2014) Meta-models as a straightforward approach to the sensitivity analysis of complex models. Population Ecology, 56, 7–19. Defeo, O., McLachlan, A., Schoeman, D.S., Schlacher, T.A., Dugan, J., Jones, A., Lastra, M. & Scapini, F. (2009) Threats to sandy beach ecosystems: a review. Estuarine, Coastal and Shelf Science, 81, 1–12. Dugan, J., Hubbard, D.M., Rodil, I.F., Revell, D. & Schroeter, S. (2008) Ecological effects of coastal armouring on sandy beaches. Marine Ecology, 29, 160–170. Fish, M.R., Cote, I.M., Horrocks, J.A., Mulligan, B., Watkinson, A.R. & Jones, A.P. (2008) Construction setback regulations and sea-level rise: mitigating sea turtle nesting beach loss. Ocean & Coastal Management, 51, 330–341. Game, E.T., Watts, M.E., Wooldridge, S. & Possingham, H.P. (2008) Planning for persistence in marine reserves: a question of catastrophic importance. Ecological Applications, 18, 670–680. Government of South Africa (2010) National Protected Areas Expansion Strategy for South Africa 2008: Priorities for Expanding the Protected Area Network for Ecological Sustainability and Climate Change Adaptation. Government of South Africa, Pretoria. Halpern, B.S., McLeod, K.L., Rosenberg, A.A. & Crowder, L.B. (2008a) Managing for cumulative impacts in ecosystem-based management through ocean zoning. Ocean & Coastal Management, 51, 203– 211. Halpern, B.S., Walbridge, S., Selkoe, K.A., Kappel, C.V., Micheli, F., D’Agrosa, C. et al. (2008b) A global map of human impact on marine ecosystems. Science, 319, 948–952. Hansen, G.J.A., Ban, N.C., Jones, M.L., Kaufman, L., Panes, H.M., Yasue, M. & Vincent, A.C.J. (2011) Hindsight in marine protected area selection: a comparison of ecological representation arising from opportunistic and systematic approaches. Biological Conservation, 144, 1866–1875. Harris, L.R. (2012) An ecosystem-based spatial conservation plan for the South African sandy beaches. PhD thesis, Nelson Mandela Metropolitan University, South Africa. Harris, L.R., Nel, R. & Schoeman, D.S. (2011) Mapping beach morphodynamics remotely: a novel application tested on South African sandy shores. Estuarine, Coastal and Shelf Science, 92, 78–89. Harris, L., Nel, R., Holness, S., Sink, K. & Schoeman, D. (2013) Setting conservation targets for sandy beach ecosystems. Estuarine, Coastal and Shelf Science, doi: 10.1016/j.ecss.2013.05.016 Harris, L.R., Campbell, E.E., Nel, R. & Schoeman, D.S. (2014) Rich diversity, strong endemism, but poor protection: addressing the neglect of sandy beach ecosystems in coastal conservation planning. Diversity and Distributions, 20, 1120–1135. IPCC (2013) Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. (eds T.F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S.K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex & P.M. Midgley), Cambridge, UK and New York City, New York, USA.


Journal of Applied Ecology 2014, 51, 1504–1514

doi: 10.1111/1365-2664.12345

Using multivariate statistics to explore trade-offs among spatial planning scenarios Linda R. Harris1*, Matthew E. Watts2, Ronel Nel1, David S. Schoeman1,3 and Hugh P. Possingham2,4 1

Coastal and Marine Research Unit, Department of Zoology, Nelson Mandela Metropolitan University, PO Box 77000, Port Elizabeth 6031, South Africa; 2ARC Centre of Excellence for Environmental Decisions, School of Biological Sciences, The University of Queensland, Brisbane, Qld 4072, Australia; 3School of Science & Engineering, University of the Sunshine Coast, Maroochydore, Qld 4558, Australia; and 4Department of Life Sciences, Imperial College London, Silwood Park, Ascot SL5 7PY Berkshire, UK

Summary 1. Scenario planning can be useful to guide decision-making under uncertainty. While sys-

tematic conservation planning can create protected-area networks for multiple and complex reserve–design scenarios, planners rarely compare different reserve networks explicitly, or quantify trade-offs among scenarios. 2. We demonstrate the use of multivariate statistics traditionally applied in community ecology to compare reserves designed under different scenarios, using conservation planning for beaches in South Africa as an example. Twelve reserve–design scenarios were run in Marxan in a hierarchical experimental design with three levels: including/excluding the probability of site destruction; two different cost types; and three different configurations of existing terrestrial and marine reserves. 3. Multivariate statistics proved to be useful tools in the conservation planning context. In our case study, they showed that the trade-off associated with including the probability of site destruction during coastal reserve design depended on the cost type: if the cost is related to the site-destruction probability then reserves are significantly larger; if not, then reserves are significantly more costly. In both cases, the configuration of existing reserves locked a priori into the solutions was more important and resulted in significantly larger and more costly reserves. 4. Synthesis and applications. This study demonstrates a novel application of multivariate statistical tools to robustly quantify potential trade-offs among diverse sets of reserve–design scenarios. These statistics can be applied: to support negotiations with stakeholders and decision-makers regarding reserve configurations in the face of uncertainty; in reserve–design sensitivity analyses; and in priority setting for future research and data collection to improve conservation plans. Key-words: complete hierarchical cluster analysis, Marxan, non-metric multidimensional scaling, protected areas, reserve design, sandy beaches, scenario planning, site destruction, spatial prioritization, systematic conservation planning

Introduction Between cumulative impacts from human activities and stressors from global change (Halpern et al. 2008b; IPCC 2013), the intensity of environmental impacts is escalating at unprecedented rates (Butchart et al. 2010). Consequently, networks of protected areas and complementary ecosystem-based management programmes beyond reserves *Correspondence author. E-mail: [email protected]

are becoming increasingly important tools for conservation (e.g. Halpern et al. 2008a; Stokes et al. 2014). The placement and design of these biodiversity priority areas thus need to be carefully planned. We focus on reserve design for the remainder of the paper, but the principles and tools apply equally to other spatially explicit ecosystem-based management programmes. Constraints in reserve design are that space and resources are limited and that multiple stakeholders compete for exclusive rights to portions of land and

© 2014 The Authors. Journal of Applied Ecology © 2014 British Ecological Society