Evolving protected-area impacts in Panama: impact shifts show that plans require anticipation

Protected areas (PAs) are the leading forest conservation policy, so accurate evaluation of future PA impact is critical in conservation planning. Yet by necessity impact evaluations use past data. Here we argue that forward-looking plans should blend such evaluations with anticipation of shifts in threats. Applying improved methods to evaluate past impact, we provide rigorous support for that conceptual approach by showing that PAs’ impacts on deforestation shifted with land use. We study the Republic of Panama, where species-dense tropical forest faces real pressure. Facing variation in deforestation pressure, the PAs’ impacts varied across space and time. Thus, if shifts in pressure levels and patterns could be anticipated, that could raise impact.


Introduction
Tropical forest loss is a major environmental problem (Benhin 2006, DeFries et al 2010, Rudel et al 2009. Longstanding concern about loss of species and essential ecosystem services has been joined by greater understanding and acknowledgment of links between forest loss and climate change, which has raised the urgency for forest conservation. In response, global efforts are ongoing to increase conservation efforts, for instance via incentives for generating 'REDD+' (reduced emissions from deforestation and forest degradation) (Corbera and Schroeder 2011). 5 Co-lead authors.
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Globally to date, protected areas (PAs) are the leading actions for conservation of forests. In humid tropical forests, they have expanded greatly in the past 30 years (Ervin 2003) to cover perhaps 20% of the area (Scharlemann et al 2010, Chape et al 2005, Jenkins and Joppa 2009. Looking forward, REDD will require contributions from both current and future protected areas. Scharlemann et al (2010) estimate the value of carbon sequestration provided by all current PAs within humid tropical forested regions to be about US$7 billion. To maximize contributions from additional future PAs-to REDD+ (Brandon andWells 2009, Campbell et al 2008), biodiversity or both-the allocation of conservation resources should reflect understanding of PAs' impacts.
Here we consider estimates of the reduction in deforestation due to the existence of a PA. Recently (see Joppa andPfaff 2010a, Miteva et al 2012), 'matching' methods for generating such estimates have been applied to improve similarity of comparison locations with PAs' locations, as PA locations have tended toward lands that face lower clearing pressures (Joppa and Pfaff 2009). Since PAs do not have impact without pressure to be blocked, addressing this location bias using more similar locations has tended to lower estimates of PA impact. PAs still reduce deforestation, on average, yet not all PAs have impact and those that do often have less than has been believed (see the results in, e.g., Ervin 2003; Andam et al 2008; Pfaff et al 2009;Joppa and Pfaff 2010b).
We apply 'matching' to PAs within the Republic of Panama with two leading objectives. Our first objective is focusing on Panama per se, improving estimates of PA impacts there by applying recent methods to estimate the baseline deforestation that would have occurred without protection. Following others' efforts (see further examples in Robalino 2012 or Pfaff et al 2013), the matching methods improve baseline estimates by explicitly trying to 'balance', i.e., equalize, characteristics such as slope, soil quality and distance to market that are identified as relevant for the observed rates of deforestation. For each protected location, we search unprotected locations for those with the most similar characteristics to permit the best 'apples-to-apples' comparisons.
Second, in considering impacts for multiple periods of time (1992-2000 and 2000-2008), we rigorously illustrate the need to anticipate deforestation pressure to plan future conservation. Impact evaluation uses data about the past. Yet we show that PA impacts have varied in the past within Panama, where deforestation pressure has shifted (Sloan 2008, Sloan andPelletier 2012, and  Thus, if policy makers can anticipate shifts in the level and pattern of pressures in whatever way (statistical predictions or knowledge of pending public investments), that additional perspective could be used to better locate additional future PAs for the purpose of maximizing their impacts. Below we present background on Panama, data methods, results and some additional discussion. From the 1960s to the 1980s, the area designated for legal protection rose by a factor of four, and now 27% of the country is protected.

Background on Panama
The 1980s were the peak including the establishment of Darien, La Amistad and Chagres parks-while in the 1990s, many small PAs were created (see IUCN and UNEP-WCMC 2009). Figure 1 shows PA locations as of 2006. Protection is somewhat fragmented and often somewhat remote from major roads and cities. In light of IUCN's distinctions between PAs by stringency, we note that the majority of the PAs in Panama are within IUCN category 'II' (National Parks), though some are in other categories. Without affecting our results, we combine protected areas.
Panama's PA system is administered by the National Environment Authority (ANAM), a government agency established in 1998. ANAM recently shifted from a 'command-andcontrol' to a more 'community-based' emphasis, including economic incentives (Oestreicher et al 2009). Interactions with other countries include projects on forest-carbon offsets supported by the CDM. Since 2008, funding for capacity building in forest governance has come from the World Bank's 'REDD readiness' program and UN REDD program (St-Laurent et al 2013). Studies on potential for REDD and offsets in Panama include Tschakert et al (2007), St-Laurent et al (2013), Potvin et al (2008), Coomes et al (2008), Pelletier et al (2011, Sloan and Pelletier (2012) and Asner et al (2013).

Forest change
Panama has a mix of tropical humid and dry forest in both pristine forests and secondary forests which are in succession. Deforestation measured within the earliest survey, done in 1947 was recorded mainly in Pacific coastal forests. Intensity of agricultural migration grew after the 1950s yet trends have shifted, e.g., recent rapid forest regeneration occurred due to abandonment of agricultural lands (alongside some direct afforestation efforts: Wright and Samaniego 2008). We consider the periods 1992-2000 and 2000-2008 because they are the most recent two distinct time periods for which data are readily available (not linked in any way to shifts in PA policies).
Deforestation continues in Panama. Between 1992 and 2000, the net rate of deforestation was 1.1% per year. Further, since it includes regeneration, that net rate masks some deforestation. During this period, mature forests were cleared mainly in Panama and Darién provinces-while clearing of secondary forest happened in Ngöbe-Buglé indigenous reserve (Pelletier et al 2011). Sloan and Pelletier (2012) identified as important drivers of the rate of 1992-2000 deforestation factors such as elevation, soil texture, slope, settlement density, rural income and protected areas. Others have found that within-district variation in deforestation during this period is correlated with further factors such as the history of land use and land title (Wright and Samaniego 2008).
FAO reports Panama's deforestation rate during 2000-2010 to be 0.4%, much lower than in the previous period (FAO 2011). The government of Panama has not reported official statistics of forest change after 2000 thus we use data provided by others (more below). These data show that excluding regeneration, the deforestation rates will be much higher than the average FAO reports (true in the prior time period as well). In terms of drivers, Sloan and Pelletier (2012) found some drivers of 1992-2000 deforestation diminished in relevance during this subsequent time period. Throughout these two time periods, the Panama Canal had a big role in the Panamanian economy, given its role in maritime trade, and looking forward that will continue to be the case. The watershed that supplies the canal encompasses about 4% of Panama's total land cover and is protected to maintain forest cover and minimize soil erosion within its boundaries (Oestreicher et al 2009). The canal's proximity to markets creates pressures for land conversion (Rompre et al 2008), so this protection has impact. Larger reserves are for indigenous tribes (Kuna, Emberá-Wounann and Ngöbe-Buglé) covering about 20% of national land-often in and around forests. The land use within these reserves can differ considerably from other areas (Nelson et al 2001).

Protected areas and forest change
Some prior research has considered the impacts of Panama's protected areas on its forests. As noted, Sloan and Pelletier (2012) included PAs as a factor influencing rates of deforestation. Oestreicher et al (2009), comparing nine protected areas, considered not only deforestation rates but also other characteristics of protected areas, in considering reasons for variations in impacts. They concluded that insufficient resources for management of PAs, as well as biases in resource allocation, limit protected-area capacities for effective forest conservation. Econometric analysis of Darien province by Nelson et al (2001) found that the protected areas there had less effect on land use than did indigenous reserves. We will not study indigenous reserves-which represent a complicated intervention, involving many dimensions and deserving its own distinct scrutiny-but focus on improving the estimates of average and heterogeneous impacts of protected areas.

Methodology
In measuring PA deforestation impact, we should be concerned about biased estimates if protected lands and their comparison unprotected lands are not similar (Joppa and Pfaff 2010a). If those two sets of locations differ significantly in terms of deforestation-relevant characteristics, then inferring PAs' forest impacts is not straightforward. To remove the influences of differences in relevant observable characteristics, ordinary least squares (OLS) regression is an option yet if control group characteristics differ greatly, the burden to control for factors can be considerable 6 .
Matching-a method widely used to measure intervention impact in various disciplines-addresses such differences in observable relevant characteristics by finding 'control' observations (unprotected lands in our application) with similar values for all of the relevant characteristics. What 'similar' means, however, varies across popular matching methods. In the propensity-score matching method we mainly employ (Rubin 1973, Rosenbaum andRubin 1983), regression of treatment on land characteristics predicts the likelihood or propensity an observation is 'treated' (protected in our application). That predicted probability is a summary of the characteristics that is used to identify 'similarity'. An alternate approach we employ to check robustness, covariate matching (Abadie et al 2004), uses distance in characteristics space to indicate the 'similarity'.
Importantly, the most similar points might still differ considerably in key characteristics. Generally, one must confirm whether the search for similarity actually yielded high similarity by directly comparing the distributions of the control and the matched treated observations for each of the characteristics considered. Also, regression analysis for the matched sample can address imperfections in matching. Another approach is to re-run analyses without the treated points that have poor matches. This drops data but, for the rest, it increases the balancing of characteristics.

Data
3.2.1. Deforestation. We use spatially detailed forest observations at three points in time (1992,2000,2008). To our knowledge, these are the latest years for which national spatially detailed land-cover data exist. Deforestation between 1992 and 2000 is measured using maps of land-cover inventories created by ANAM for each year. Landsat 5 and 7 imageries were used to create maps with 100 m resolution (Pelletier et al 2011). Deforestation between 2000 and 2008 is measured using a 2008 land-use map from Sloan and Pelletier (2012). That map was derived from both Landsat 7 ETM+ and ASTER satellite images which were all recorded around the year 2008 and had a resolution of 15-30 m. 7 The vegetation classification for ANAM and Sloan and Pelletier (2012) is comparable over time with the exception of a new forest-plantation category defined for 2008 (1% of forest cover). We use the forest definition in Sloan and Pelletier (2012) by including mature, secondary and intervened forests, the latter defined by forest cover altered by human activity by over 60%. Definitions of mature, secondary, and intervened correspond to the legal government definitions (la Ley Forestal, resolución J.D. 1998). Since treating intervened forest as forest has a risk of missing forest degradation, we also redid all our analyses treating that intervened category as being in 'non-forest'-not one of our central conclusions was affected by that shift.
3.2.2. Protected areas. The majority of PAs were established before 1992, so we evaluate most PAs in Panama. For each period, we drop the PAs established after 1992, i.e., we examine a consistent set of PAs over time so that any shifts in estimated PA impact are due solely to shifts in land-cover change. We also drop the indigenous reserves (Kuna-Yala, Kuna de Wargandi, Kuna de Madungandi, Embera-Woonan, and Ngöbe-Bugré), since they represent a very different intervention from the protected areas. That leaves 42 755 km 2 (57% of Panama) to study for 1992-2000 deforestation and 42 164 km 2 to study for 2000-2008 deforestation. Roughly 40% of that sample is within PAs.
3.2.3. Parcel characteristics. The factors that we consider as potential drivers of deforestation include elevation, slope, distances from roads, urban areas and rivers, provinces, lifezones/ecoregions and agricultural suitability. Elevation and slope were calculated using the Shuttle Radar Topography Mission elevation data provided by the Global Land Cover Facility. Road locations were identified based upon the VMAP Level0 dataset. Urban locations were defined based upon the Gridded Rural Urban Population (GRUMP), v1 data that has been generated by CIESIN (CIESIN et al 2011).
In order to include broad spatial controls for unobserved differences across the locations, we use indicator variables for provinces (although we run robustness checks without provinces), as well as indicators for mountainous and dry forests that are derived from the map of lifezones (Sloan and Pelletier 2012). Finally, the variable for agricultural suitability was generated by the International Institute for Applied Systems Analysis (discussion provided in Fischer et al 2002). The index provided is estimated using climate, soil type, land cover and slope. From the original integer data-running most suitable (0) to least suitable (9)-we created a dummy for range 0-4, i.e., for a relatively low value which indicates that the soil is relatively beneficial for agriculture.

Protection and deforestation statistics
4.1.1. Protected locations (versus unprotected). Table 1A shows that protection in Panama, as in many other countries around the world (Joppa and Pfaff 2009), on average is located on lands likely to face lower pressure for clearing.
At the top of the table, protected lands' elevations and slopes are considerably higher on average (note that these two variables are relatively highly correlated, so not used jointly in regressions). Protected forests are also farther from roads in particular, as well as cities, and on lands with lower agricultural suitability-all in comparison with the forests on lands that are not protected.
We also consider the fraction of land in each province. Provinces do differ in terms of the deforestation drivers, e.g., some have higher agricultural suitability and others are far from cities. Provinces also differ in ways not described by these characteristics and, to the right in table 1A, we see that they vary considerably in deforestation. Starting on the central theme of this paper, we see the deforestation pattern across provinces shifted significantly between the time periods.
Table 1B adds perspective on where protection was located. A probit regression explains what lands were protected, conditional on province (marginal effects on treatment are displayed). This confirms from table 1A that protection is on higher slopes and farther from roads but, even for those factors, differences in coefficients are smaller than the average differences in table 1A. Further, the coefficient on distance from urban area, as well as the one of agricultural suitability, shows the opposite relationship to protection's locations compared to that suggested in table 1A. The reason is that provinces differ in characteristics. Thus, the average differences in table 1A appear in table 1B in two forms: province effects; and coefficients describing where protection tends to be located within a province. Given significant province differences in where protection was located and significant  Starting from the bottom, as noted just above it is important to control for the provinces. That is confirmed by the significant coefficients for many provinces in both of these regressions (noting also that some provinces are dropped as they are very highly overlapping the demarcated indigenous lands-suggesting that it may be difficult to separately test the indigenous reserves). We also highlight significant changes in provinces' effects over time-consistent with table 1A. Three provinces' coefficients are significant in both time periods but differ in sign across periods. Other provinces' coefficients are significant in both periods but shift considerably in magnitude.
Other observed drivers of deforestation also vary over time in their estimated impacts and recall that these impacts are conditioned on provinces, i.e., refer to patterns within provinces; that means some of the average impact differences for these factors are seen in the provinces' effects. While slope discourages deforestation in both periods, and agricultural suitability's coefficient is quite stable over time, the coefficient for urban distance changes in sign across the time periods. Further, the significant discouraging effect of road distance for 1992-2000 vanishes across time.

Protection's impacts on deforestation-means and regressions
4.2.1. Simple means (by time period). An initial, simplest estimate of the impact of protection on deforestation was provided by table 1A. Comparing clearing during 1992-2000 for those locations that were forested in 1992 but unprotected (19%) with clearing in the same time period for protected forested lands (1%) suggests that almost all clearing was blocked-starting from the considerable baseline of 19%. The question remaining after this simple comparison is whether that difference in deforestation between unprotected and protected lands is due to the protection or to differences in those lands. For 2000-2008, this simplest approach to estimating protection's deforestation impact suggests even more avoided deforestation-with more clearing of protected lands but a higher baseline.

Regressions (by time period).
However, we must expect that not all of those differences are due solely to the protection. From the top of table 1A, protected locations' characteristics themselves suggest lower pressure for deforestation in those locations even had they not been protected. Further, from the bottom of table 1A, looking across all of the columns suggests that the distribution of protection across the provinces, relative to the distribution of unprotected lands, also could help explain clearing rates. Table 3 reports an OLS regression with the deforestation dummy as a dependent variable. Its results, in comparison to the simplest estimates of PA impact derived from table 1A, confirm that differences in location, or land characteristics, were responsible in part for those differences in deforestation rates between protected and unprotected in table 1A. Controlling for the drivers, covariates and provinces, the coefficients for protection in table 3 indicate lower-although still quite significant-effects in each time period, specifically 14% and 18% for these time periods (and a formal test shows these two effects differ from each other-also true without provinces).

Protection's impacts on deforestation-matching
4.3.1. Matching balances. Given the differences between protected and unprotected locations (tables 1A and 1B) in terms of factors that also influence deforestation (table 2), matching could add to the regression analysis in table 3 as another method to control for the influences of differences in the locations. However, as noted above, even the most similar unprotected locations are not always so similar 8 . Table 4 shows that matching increased the similarity of unprotected locations to compare with the protected (yet the last row shows a few protected locations did not have good matches). The final column conveys that while the matched unprotected locations are not identical in their characteristics to the protected locations, pre-matching differences become much smaller after. This may not change the impact estimates but it reduces concerns about the burden of control.

4.3.2.
Matching estimates-average impact. Table 5 presents matching estimates of average PA impacts compared to other estimates, by time period. The top two rows show the estimates derived by simply comparing mean rates of deforestation on protected and unprotected lands (from table 1A) and from regression (table 3). In comparison to the initial estimates lacking controls, not surprisingly these matching estimates are more like the regression estimates (and formal tests confirm differences across time periods). Put another way, controlling for differences in characteristics between protected and unprotected locations reduces estimated impacts by something approaching one half (robust to not using the  Table 6 starts by confirming, for Panama, a common result concerning locations' slopes: PAs on flatter land have higher impacts because there is greater clearing pressure to be blocked. This holds in both periods-as we would expect given slope's consistent effects within table 2. While a similar comparison using distance to road has much less planning guidance in this case, quite striking is the comparison for urban distance-which again follows directly from table 2 where there is a dramatic change in sign in the effect of urban distance controlling for province. Table 6 shows that for 1992-2000 deforestation, PAs farther from urban areas had higher impact. However, for the 2000-2008 land-cover change, we no longer see this spatial impacts difference and The difference in estimated impacts across time periods here is significant as pooling periods and estimating an interaction of protection and period confirms a 13.6% effect in 1992-2000 and 18.3% effect in 2000-2008, while adding that the p-value on the 4.7% difference in effect is 0.000. Also, as noted above these results are not changed by using provinces. c The difference in estimated impacts across time periods here is significant as pooling periods and estimating an interaction of protection and period confirms a 10.6% effect in 1992-2000 and 14.5% effect in 2000-2008, while adding that the p-value on the 3.9% difference in effect is 0.001. Also, removing provinces from these specifications has little influence for the first time period while raising the estimate for the second period, which emphasizes that PAs' deforestation impacts varied across periods. thus the guidance to planning future conservation depends on which past predicts the future.

Discussion
We applied matching methods to improve PA-impact estimates for Panama, while also demonstrating that changes in land use, over time, shifted the deforestation impacts of the PAs. The latter result calls into question any assumption that the future will be like any particular part of the past. That, in turn, should lead policy makers to consider any basis (empirical forecasting, public plans) upon which they might anticipate how patterns of deforestation pressure may shift.
The question 'what will be the future distribution across space of deforestation pressure, and thus PA impact?' can generate informative consideration. For instance, it raises possibilities such as 'deforestation hotspots' moving to pristine frontiers, or crops' relative prices changing. Whether modeling is formal or not, some form of attempt to anticipate future pressures seems to have value. Sufficient consensus, even if informal, could shift planning estimates of PA impacts.
It is worth emphasizing that varied and quite distinct approaches might be employed and indeed blended in anticipating future threat. One might attempt to use the limited past data over time that has spatial detail in order to fit a predictive statistical model-although that is a real challenge. Alternatively, without any statistics a policy maker might well have knowledge that there will be future public investments in roads, ports, or other infrastructures that might affect deforestation.
That question might in turn lead to another: 'what policies affect pressure's evolution?'. Almost surely shifting public investments in new development infrastructure could affect threats (e.g., Pfaff and Robalino (2012) consider beneficial forest spillovers from Brazilian Amazon PAs, consistent with dynamic models yielding multiple equilibria in migration and public investment). Naturally, which form of anticipating change is appropriate, and its planning implications, varies. Yet the results presented here for Panama, with change in both pressure and PA impact over time, motivate quite generally both past PA evaluations and anticipation of future changes in pressure.