Selection biases and spillovers from collective conservation incentives in the Peruvian Amazon

In 2012, 51% of total greenhouse gas emissions in Peru originated from deforestation in the Amazon [24], primarily driven by shifting agriculture [25], gold mining [26], and cash-crop plantations such as oil palm [27] and coca (Erythroxylum spp) [28]. Estimates of deforestation suggested an increasing trend [18], with an average of 160 000 ha per year between 2011 and 2016 [21]. As a contribution to climate change mitigation, the Government communicated a zero-deforestation target to the United Nations Framework Convention on Climate Change by 2021 [29]. In 2010, the Peruvian Ministry of Environment created the NFCP 'to contribute to the conservation of tropical forests and the generation of income for the most vulnerable, poor and marginalized peoples' [30] (author's translation). The NFCP seeks to: (i) map forestlands, (ii) promote sustainable production systems, and (iii) strengthen forest conservation capacities [30]. Given the government's lack of experience in paying cash to landholders for not deforesting, conditional 'projects' had to be implemented to provide local compensatory benefits, while also striving to 'green' local livelihoods. This collective PES-cum-ICDP intervention intended to align conservation with poverty alleviation goals, piloted in selected Amazon indigenous communities [31]—some of the poorest population groups in Peru [32]. From the approximately 1300 titled native communities [32] controlling roughly 12 million hectares of forests (figure 1), 50 communities were enrolled between 2011 and 2013 for the pilot phase (table 1). These communities were selected non-randomly, using criteria ranging from forest conditions to accessibility indices [33], and subsequently applied at two spatial-administrative levels: first, at the province level (second highest sub-national political unit in Peru), and second, at the community level. The logic behind this approach was to first prioritize the provinces with the highest threats of deforestation and then to select communities within them. However, eventually those criteria were not implemented consistently and transparently (SM-Targeting available online at stacks.iop.org/ERL/14/045004/mmedia), leaving room for discretional targeting decisions. Together with the fact that communities voluntarily decide to participate (SM-Engagement), institutional selection created de facto a source of adverse selection bias [34]: as we show under Results, communities with historically higher deforestation were underrepresented in the NFCP.

Zoom In Zoom Out Reset image size

The NFCP provides collective payments of 10.00 Peruvian Soles (1.00 Peruvian Sol ∼ USD 0.29 in 2015) per year and hectare of forest enrolled under five year contracts, supplemented by technical assistance. The payment is publicly funded and is conditional upon (i) its spending on a collectively agreed investment plan to finance forest-friendly production (e.g. agroforestry, aquaculture, and small animal husbandry), community forest patrolling, and public services or infrastructure, and (ii) the maintenance of forest cover in 'conservation forest zones' (CFZ) that communities define themselves. This community self-selection of land constitutes a second source of adverse selection bias [34]. The remaining land, i.e. 'other use zones' (OUZ), is not subject to land use restrictions and typically contains homesteads, agricultural fields, and secondary as well as primary forests remnants (figure 2). Our empirical strategy seeks to measure the NFCP's impact during its five initial years (2011–2015) on deforestation in the community lands, as a whole, and in both the CFZ (primary effect) and OUZ (spillovers). We use the term 'spillover' to denote that the NFCP's intervention to avoid deforestation, targeted at the enrolled CFZ, may also have indirect yet measurable impacts on the unenrolled subareas of the treated communities (OUZ).

Zoom In Zoom Out Reset image size

We created two datasets, one for community polygons and one for cells. The first is comprised of 992 communities (50 treated and 942 non-treated), for which we could gather geographical, biophysical, and socioeconomic attributes (table S2 in the SM). The second only includes cells that partially or fully overlap with polygons representing community lands, including those in CFZ and OUZ (figure 2), and that present a deforestation risk larger than 1% (SM-Forest risk model). We developed a deforestation risk model to cope with the fact that the distribution of the observed annual deforestation is skewed towards zero and thus to focus only on cells with a minimum risk of 1% of having been deforested between 2001 and 2010. We fit a logistic regression model, where the outcome is measured as the fraction of a cell deforested between 2001 and 2010, and covariates include: population, number of houses, area of coca plantations, number of population centers, community area, slope, precipitation, distance to protected areas, forest loss density in 2010, distance to population centers within communities, internal distance to the community's boundaries, share area of forest in 2010, spatially lagged biomass, and distance to deforestation outside communities in 2010. Using the fitted values of the model, we discard all cells with a fitted value smaller than or equal to 0.01. This effectively reduced the number of total cells and treated cells with no deforestation in any year between 2001–2015 to only 11%. Only the trimmed cell-dataset is used for impact evaluation.

We define the outcome variable as the total area (ha) of annual gross deforestation within community lands for both units of analysis. We used a deforestation dataset, covering 15 years (2001–2015), within the Peruvian Amazon [18–21], and define deforestation as the complete removal of forest cover from a Landsat pixel [18].

We use a quasi-experimental approach that combines matching and double difference regression models to estimate the average treatment effect on the treated (ATT). As a pre-processing approach, matching reduces the selection bias due to the non-random selection of enrolled communities and CFZ, an approach that is increasingly used to measure conservation outcomes [36, 40–45], including in the context of the Peruvian Amazon [32, 45, 46]. We measure effects at different spatial scales in both enrolled and non-enrolled community lands (SM-Treated units). We start by comparing outcomes at the scale of decision units using spatial polygons and cells (225 ha each) of treated and untreated-communities (SM-Cells and table S1). We chose an aggregated area of 225 ha because the minimum sizes of (i) treated community polygons and (ii) CFZ-polygons are 2700 ha and 1500 ha, respectively. Thus, there are at least approximately 12 and 6 units covering each of these polygons. We believe that aggregating the outcome variable from 0.09 to 225 ha improves the representation of the land use decision unit and reduces potential bias arising from spatial autocorrelation [15, 47].

This approach, however, aggregates over intra-community effects, given that payments restrict deforestation only in CFZ. Avoided CFZ deforestation could be partially outweighed by added OUZ deforestation (negative spillover, i.e. leakage). Alternatively, participants could shift resources to fulfill the NFCP's rules and thereby reduce deforestation in OUZ (positive spillover). Thus we also estimate the ATT in CFZ and OUZ cells, separately. For this, we model at the cell level which areas of non-participating communities would have most likely been selected as CFZ and OUZ from which to withdraw a control group (SM-Modeling untreated CFZ and OUZ).

We use one-to-one nearest neighbor matching with replacement to find a control group from the pool of untreated communities and cells (SM-Matching), using the Genetic algorithm [48] of the R Match function and the Mahalanobis distance for polygons and cells, respectively, and the exact distance for the Department identifier in both. We use a set of covariates from a pre-treatment data set of geo-biophysical, land-use and land-cover, infrastructure, and socioeconomic variables (table S2 and Covariates in the SM). These covariates include¹ : (i) biophysical: elevation^† (m), slope^† (°), above ground live woody biomass^†† (Mg/ha), temperature^† (°C), precipitation^† (mm), distance to rivers^† (m); (ii) infrastructure: distance to roads^† (m), accessibility^†† (index), distance to district's capitals^† (m), distance to population centers^† (m); (iii) land use/land cover: forest cover area in 2010^†† (%), deforestation density in 2010^† (ha km^–2), community´s total area^††† (ha), distance to deforestation outside communities^††† (m), distance to protected areas^†† (m), internal distance to community´s boundary^† (m), deforestation risk^†; (iv) spatial lags^††† for: deforestation (ha), forest cover area in 2010 (%), slope (°), above ground live woody biomass (Mg/ha), elevation (m); (v) socioeconomic: density of coca plantations^†† (ha km^–2), years passed since communal land titled^††† (years), population^†† (person), number of houses^††† (house), access to drinking water^††† (%), access to electricity^† (%), population centers within a community^†† (center), per capita income^†† (PEN), human development index^† (index), total poverty^†† (%), and extreme poverty^††† (%). These covariates are likely to affect both the selection of participating communities and deforestation, but are not affected by the treatment [53]. Treated and control units were drawn from a total number of 992 communities (polygons) and 18 319 community-cells (cells). The number of unique treated and control units for each level of analyses is: communities (polygons, 50 treated and 36 untreated); communities (cells, 986 treated and 495 untreated); CFZ (cells, 523 treated and 304 untreated); and OUZ (cells, 655 treated and 305 untreated). The treatment variable is continuous from 0 to 1 in both polygons and cells as units of analysis. In the first case it denotes the fraction of the year in which a community polygon has been treated. In the second case it denotes the share area of polygons representing treated communities/CFZ/OUZ within a cell and the fraction of the year in which the cell has been treated. The control units always have a treatment value of zero.

We use the matched dataset and apply a fixed-effect regression model to a panel dataset of 15 years (2001–2015). The original model at the community level is:

$$\begin{eqnarray}&&\begin{array}{l}{Y}_{ctd}=\beta {D}_{ct}+t{X^{\prime} }_{c}\delta +\gamma Out\_pre{s}_{ct}\\ \,\,+\,{\varphi }_{t}+{\alpha }_{c}+t{\omega }_{d}+{u}_{cdt}.\end{array}\end{eqnarray} \tag{ 1 }$$

And at the cell level:

$$\begin{eqnarray}&&\begin{array}{l}{Y}_{ictd}=\beta {D}_{ict}+t{X^{\prime} }_{i}\delta +\gamma Out\_pre{s}_{ict}\\ \,\,+\,tW^{\prime} {Z^{\prime} }_{i}\lambda +{\varphi }_{t}+{\alpha }_{i}+t{\omega }_{d}+{u}_{icdt},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${Y}_{ctd}$ is the annual deforestation in each community c, in year t, located in department d (highest sub-national political-administrative unit in Peru). ${D}_{ct}\,$ is the treatment indicator at the community level, turning positive in year t when a community is enrolled (equation (1)). Enrolled communities become a treatment indicator larger than 0 and up to 1 in the year of enrollment, depending on the number of treated months during that year. Thereafter, the indicator is 1 if the community is still treated. Similarly, ${Y}_{ictd}$ is the annual deforestation in each cell i, located in community c, in year t, in Department d. ${D}_{ict}$ indicates treatment at the cell level using the treated share of cell i in year t (equation (2)). Cells within communities have a share of one, and cells at the margin of community borders have a share between 0 and 1. X' is the vector of time-invariant pre-treatment characteristics, such as slope. These variables can influence the deforestation trend, and are therefore interacted with time indicator t. $Out\_pre{s}_{it}\,$is the average distance from community i in year t, to deforestation patches (>1 ha) located outside communities, and represents external deforestation pressure. In equation (2), $WZ^{\prime}$ represents a vector of spatially lagged covariates weighted by a standardized queen contiguity matrix (W) [22]. δ, γ, and λ are coefficients to be estimated. Year fixed effects, denoted by ${\varphi }_{t},$ control for yearly factors influencing all units of analysis equally, such as, policy changes to the Peruvian forestry rules [32]. Individual fixed effects (${\alpha }_{c}$ or ${\alpha }_{i}$) represent the individual unobserved time-invariant heterogeneity (e.g. soil quality). ${\omega }_{d}$ is introduced to capture department-specific forest conservation efforts. Finally, ${u}_{cdt}$ or ${u}_{icdt}\,$denote the idiosyncratic errors [54].

Taking first differences (FD) of equations (1) and (2) eliminates all time-invariant unobserved heterogeneity (${\alpha }_{c}$ and ${\alpha }_{i}$), which could have biased our estimates (SM-Specification Tests). The coefficient β represents the ATT on annual deforestation changes for all years after treatment. For communities, we estimate the FD of equation (1):

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{Y}_{cdt}=\beta {\rm{\Delta }}{D}_{ct}+{X^{\prime} }_{c}\delta +\gamma {\rm{\Delta }}Out\_pre{s}_{ct}\\ \,\,\,+\,{\rm{\Delta }}{\varphi }_{t}+{\omega }_{d}+{\rm{\Delta }}{u}_{cdt}.\end{array}\end{eqnarray} \tag{ 3 }$$

For cells we estimate the FD of equation (2):

$$\begin{eqnarray}&&\begin{array}{l}{Y}_{ictd}=\beta {\rm{\Delta }}{D}_{ict}+{X^{\prime} }_{i}\delta +\gamma {\rm{\Delta }}Out\_pre{s}_{ict}+W^{\prime} {X^{\prime} }_{i}\lambda \\ \,\,\,\,\,+\,{\rm{\Delta }}{\varphi }_{t}+{\omega }_{d}+{\rm{\Delta }}{u}_{icdt}.\end{array}\end{eqnarray} \tag{ 4 }$$

By comparing fixed and random effects models using the Hausman test [55], we rejected the null hypothesis (p < 0.01) in all cases, thus indicating that a fixed effects model was more appropriate. Given that treatment assignment occurs at a higher level than the cell, namely, at the community level, we cluster standard errors of $\beta$ at the level of the community polygon to avoid inconsistent variance-covariance matrices due to heteroscedasticity and autocorrelation in the error terms [56, 57] (see SM-Specification tests). In addition to the overall effect across 2011–2015, we also estimate the effect of the NFCP over the years after enrollment by replacing the treatment variable in equations (3) and (4) with five new treatment variables, each one representing year zero through year four after enrollment. The estimated coefficients for each of these variables are interpreted as the average effect of the NFCP after t years of enrollment [58]. Used models are presented in SM-Effect over time.

After matching, balance between treated and non-treated groups improved in almost all covariates for all levels of analysis (tables S3–S6 in SM). We use the normalized difference, the mean difference in the standardized empirical-QQ plot (as proposed by [59]) and the mean difference between treated and control groups to assess covariate balance after matching. For a few time-invariant covariates normalized differences remained above the rule of thumb of 35% even after matching [60]. This potential source of bias is addressed by using the fixed effects model to estimate treatment effects [61].

Annual averages of deforestation between 2001 and 2015 of all treated communities, non-treated communities, and the matched control group increase over time (figure 3). The NFCP predominantly selected communities with relatively lower deforestation threats. Figure S2 in the SM confirms this by comparing three alternative measures of pre-treatment deforestation levels among the 50 treated communities and top-50 non-treated communities.

Zoom In Zoom Out Reset image size

Trends were similar at the cell level, where treated cells, CFZ-cells and OUZ-cells, and their corresponding matched control groups, exhibited increased deforestation (figure 3). Non-zero, but low levels of deforestation in the CFZ of treated communities (figures S3 in the SM) also suggest mild levels of non-compliance in most cases. In eight communities, however, the percentage of the deforested area relative to the area of the CFZ exceeds the threshold of 0.3% [62], above which a community is allegedly to be evicted from the NFCP. Nonetheless, from these eight communities, only two were expelled in 2014, the rest remained enrolled. Another eight communities were also evicted (table 1), but due to causes not related to deforestation in their CFZ (e.g. non-compliance with their investment plan).

In addition, we observe that there are substantial differences between the CFZ and the OUZ of participating communities regarding characteristics that could affect both the selection of an area as a CFZ and deforestation outcomes. These include: slope, elevation, distance to rivers, distance to population centers within the community, deforestation previous to the start of the NFCP (2010) and deforestation risk (table 2).

The ATT of the NFCP at community scale is statistically significant and negative (table 3, column 1) at polygon scale (equation (3)), but insignificant at cell scale (table 3, column 2, and equation (4)). However, our results in columns 1 and 2 might have been affected by the independent effects occurred within CFZ and OUZ leading to a relatively less precise ATT [15]. Therefore in column 3 and 4 we present the independent effects of the NFCP within CFZ and OUZ, respectively. Assessing intra-community effects (equation (4)), we only find statistically significant and negative effects in OUZ, not in CFZ (table 3, column 5). This implies that the NFCP may have avoided deforestation within OUZ-cells by an average of 0.4 ± 0.2 ha y⁻¹ (Mean±SE), in every subsequent year after treatment. This estimate represents a total of 557 ha (considering SE: 59-1, 056 ha) of avoided deforestation between 2011–2015, corresponding to a 5.8% reduction (0.61%–11.1%).

In addition to the overall effect for each scale of analysis, we also explore how the effect evolves over time (figure 4). When analyzing community effects, we find statistically significant and negative ATTs only in the first year after enrollment at the polygon (−8.5 ± 3.5 ha) and cell scales (−0.27 ± 0.13 ha). When analyzing intra-community effects at the cell scale, we only find significant and negative ATTs in the second (−0.21 ± 0.1 ha) and first years after enrollment (−0.45 ± 0.17 ha) for CFZ and OUZ cells, respectively.

Zoom In Zoom Out Reset image size

These results do not change our previous conclusions regarding overall effects [58], but provide additional clues to understand the potential impact channels. The fact that significant effects throughout are only present in the initial years, dissipating thereafter, suggests that the NFCP might have induced a behavioral change, but probably only for a short period.