Study site and context
The BV program was created in 2011. Its objectives were to increase and diversify the incomes of poor households in eligible conservation sites, while incentivizing conservation and sustainable management of ecosystems 53–56. Eligible sites consisted of i) inhabited federal protected allowing the sustainable use of natural resources, ii) settlements allocated by agrarian reforms to traditional populations involved in agro-extractive activities and sustainable forest management, and iii) recognized traditional communities living on Brazilian public lands (Fig. S1).5 The three different tenure systems are respectively regulated by distinct Brazilian federal agencies, and property rights are subjected to legal restrictions and an approved management plan on the use of natural resources.6 Within eligible sites, eligible households are those with livelihoods reliant on natural resources, and a monthly income per capita under 85 Brazilian reais (~ 24 USD in 2016). Households interested in enrolling on the program also had to be included in the governmental list of households eligible to receive governmental social subsidies (Bolsa Familia program). At the beginning of the program (October 2011), BV initially enrolled households from eligible sites in the Legal Amazon, corresponding to Amazonia and some parts of the Cerrado biome. BV then expanded, in 2012, to other biomes in Brazil (Caatinga, Cerrado, Marinho Costeiro, and the Mata Atlântica) while continuing to impose similar social criteria in terms of eligibility.
Between 2011 and 2018, almost 80,000 households enrolled, though some for less than the full duration of the program (Fig. S2). They had to comply with a list of activities that foster sustainable management and conservation, in exchange for quarterly payments of 300 Brazilian reais (about USD 86 in 2016). Complementary funding supported training in diversifying sources of income and complying with local environmental regulations, which was monitored by officers 37–39. The number of enrolled households started to fall after 2016, notably following the non-renewal of many contracts because the implemented agency argued that many households were no longer meeting the extreme poverty criteria. BV ended in early 2018. Many sites experienced some deforestation during the implementation of the program (Fig. S3) and particularly after 2015, given a gradual weakening of environmental law enforcement in all Brazil 40–42.
Reports based on data from BV sites note positive socioeconomic impacts (via cash transfers, training, income diversification) as well as environmental impacts (via contracts’ incentives, clarification of rules governing the management of natural resources, and better reporting of illegal activities) 37,38. One rigorous impact evaluation concluded that BV reduced forest loss by about half of the counterfactual for Settlements and PAs in Amazonia during 2011 and 2015 39, while causally attributing most of the impact not to the enrolled households’ behaviors but instead better monitoring of and increased fines for illegal deforestation by outsiders. An evaluation using a socio-economic survey in 2015, for all biomes, also concluded BV did not impact enrolled households’ logging practices but, rather, raised their familiarity with environmental regulations and, for supporting local livelihoods, facilitated their acquisition of new productive equipment 37. Our analysis complements these studies by considering BV’s entire duration and distinguishing its impact by biome, while also separating the sample of settlements from protected areas. We extend all prior studies by considering the possibility of post-program impacts on forests.
Evaluating BV Impacts
We evaluate impact at the site level, given BV’s individual and site-level elements (e.g., both payment transfers and training to increase the compliance with site management plans). Examining the polygons for inhabited federal protected, and those for settlements allocated to traditional populations, then comparing them with respective control sites, is consistent with BV having impact via these conjectured causal mechanisms: a) reduction of deforestation in enrolled households’ properties resulting from a compensation for opportunity costs; b) diversification of income that reduces extensive land-use change; c) better awareness and compliance with local site management plans; and d) better reporting and sanctioning of illegal deforestation activities that are carried out by trespassers. As BV did not include all the households in each site, impact may involve raising compliance with site regulations regarding resources management − especially through increased reporting of illegal deforestation by non-BV actors.
Our key assumption is that, without BV, BV’s enrolled sites would have had on average the same loss as “the most similar” unenrolled sites. The validity of this assumption depends upon how similar the controls are, i.e., whether they represent the counterfactual situation of what would have happened in the absence of BV. Unenrolled sites are selected from a sample of sites of the same tenure system and located in the same biome that each sample of enrolled sited. We further select sites with the most similar values of a summary statistic for relevant characteristics, the estimated probability (or ‘propensity’) of being enrolled in BV, knowing that the program targeted poor households whose livelihoods depend on natural resources. Propensity scores are one way to process the sample selection for “apples-to-apples” comparisons of enrolled with unenrolled sites, such that post-processing differences in group outcomes are reasonable estimates of BV impact.
BV has been targeting distinct tenure systems and has been implemented in several biomes of Brazil (see S4 for covariate distribution across tenure systems and biomes). Rigorous studies of the impacts of conservation interventions have shown that the impacts of a single intervention are heterogeneous across regions, notably because the magnitude of deforestation threats is not the same across regions and sites9,44,58. Measuring the average impact of the whole program, as often measured in impact evaluation of small-scale conservation programs, would not be informative in our case given this regional heterogeneity.7 Thus, in our analysis we distinguish Amazonian sites from sites in other biomes (Caatinga, Cerrado, and Mata Atlântica) not only because of the distinctive drivers of deforestation and regulations in the former but also due to the largest number of enrolled households and area included in the program, as well as earlier time of initial enrolment in the program. We also distinguish the impact of the program on protected areas from the impacts of the program on settlements, as these two distinct tenure systems are likely to differ in terms of magnitude of the drivers of deforestation and regulation (Figs. S3 and S4). We therefore evaluate BV’s impacts by biome and tenure system: Amazonia PAs; non-Amazonian PAs, Amazonia Settlements; and non-Amazonian Settlements.
In each region, we evaluate the cumulative BV impact (i.e., the sum of annually avoided tree-cover loss during program implementation 2011-2017) by computing a generalized Difference-in-Differences, i.e., comparing changes over time in tree-cover loss (expressed in percentage of site forest cover in 2000) of in enrolled sites with those in similar unenrolled sites. We take into account that enrolled sites started in different years and that BV’s impacts could vary across years, e.g., due to shifts in local deforestation pressures. We also evaluate the post-program (2018-2020) differences in annual forest cover loss.
Unit of analysis
BV targeted households below the poverty line defined by the program and located in eligible sites. Individual recipients’ plots were not publicly available, but the Ministry of Environment’s publicly available database provided the name of the sites containing enrolled households. We identify 70 PAs and 902 settlements containing beneficiary households. We match this database with the shapefile of all federal PA recognized by the Instituto Chico Mendes de Conservação da Biodiversidade, ICMBio, which is the Brazilian Ministry of the Environment’s administrative agency (version of September 2018) as well as the shapefile of all the settlements created by the Instituto Nacional de Colonização e Reforma Agrária, INCRA, which is the Brazilian agency administering land reforms (version of September 2018). We remove all sites that were only legally recognized by the government after 2011 as well as State PA and federal "strict" PA because these two types were not eligible to the program. We also filter data on Settlements to only keep the management categories eligible to the program (“environmental distinctive settlements). Our working sample contains 113 PAs and 6048 settlements, respectively 60 enrolled PA and 53 non-enrolled PA as well as 843 enrolled settlements and 5205 non-enrolled settlements. In Amazonia, our working sample of enrolled sites includes 51 PA and 299 settlements.
We use a sample of never-enrolled sites that correspond to sites belonging to same types of sustainable use federal PAs and environmentally distinctive settlements eligible to the program, but where no households have been enrolled in the program. Explanation about why no household has been enrolled in the program in these sites includes the fact that no household met the poverty eligibility criteria, that the sites were not properly registered in the social program database, or that households were not interested in joining the program. Never-enrolled sites faced in principle the same obligation to comply with sites’ management plans than enrolled sites, even if compliance and/or enforcement may vary across time and between sites, independently of the presence or absence of BV.
Our analysis at site level allows us to consider eventual internal leakages, i.e., displacement of deforestation from one plot enrolled in the program to another plot within the same site. However, our unit of analysis may in some cases underestimate impact if only a small proportion of each site has been affected by the program.
Treatment status
The pre-matched sample of non-enrolled sites corresponds to sites where no households have been enrolled in the program but belonging to categories eligible to the program. Within each enrolled site, we compute the number of enrolled households each year. We do not distinguish the month of enrollment within each year but we “group” sites according to the year of first enrollment in the program. We also do not distinguish impacts as a function of the evolving number of enrolled households within each site.
Outcome variables
We used MapBiomas version 6.0 as a data source for computing tree cover loss. MapBiomas provides high resolution (30*30m pixels) annual raster built from Landsat satellite images. MapBiomas covers the six biomes of continental Brazil (Amazônia, Cerrado, Caatinga, Mata Atlântica, Pantanal and Pampa). We extracted pixels corresponding to natural forests (Forest, Savanna, Mangrove, and “restinga” coastal forests) and excluded forest plantations and other non-forest land covers. We excluded reforested pixels (i.e., any pixel that is classified as natural forest in a year while it was not classified as natural forest the previous years) to focus only on tree cover loss. We removed site that contained less than 10 ha of forests in 2010. We only kept annual raster corresponding to the years between 2000 and 2020 and calculate the outcome variable as the annual tree cover loss divided by the forest cover in 2000 in each site expressed in percentage for all the (enrolled and not enrolled) PA and settlements in our database. The program legally ended at the beginning of 2018, when the last enrolled households received their last quarterly payments. Therefore, our main econometric specifications are performed using the year 2017 as endline, while the years 2018 to 2020 corresponds to post-intervention years.
Confounding variables
We calculate the travel time to nearest city of 50 000 or more people in 2000 (Joint Research Centre of the European Commission‘s global map of Accessibility). We also measure the linear distance between each conservation site‘s centroid and the nearest road (using a database of road infrastructure in 2005, produced by the Brazilian Institute for Geography and Statistics (IBGE), as well as to the nearest river World Rivers from Harvard’s WorldMap. We add a discrete map of topographical regions with 5 classes: flat (less than 3% slope), undulated (3 to 20 % slope), very undulated (20 to 45%), mountainous (45 to 75%), and rugged (more than 75%) derived from IBGE’s Digital Elevation Model. We calculate the total forested area in 2010 in hectare for each site. To reflect the temporal variations of deforestation, we also calculate pre-program tree cover loss split into several periods in our base models: 2000-2004, 2004-2007, and 2007-2010, and for each of these periods we divided the total loss by the forest cover in 2000, expressed in percentage. We also used the 2010 IBGE census to obtain the Human Development Index, although this metric is at the municipal rather than site level. In the case of settlements, we also use household density (number of registered households on site divided by forest cover in 2010). The number of registered households per site is not available systematically for PA, so we do not use this confounding variable for PA systems. Finally, due to the small sample size for PAs outside Amazonia, the matching variables are limited to the distance to roads, and the 3 periods of past tree cover loss.
Comparing annual tree cover loss in enrolled and never enrolled sites
To adjust for factors affecting the probability of deforestation, we compute characteristics for each conservation site polygon in our database, and these variables are used to estimate a Propensity Score, that is the probability that a site would contain households enrolled in the program given the set of observed characteristics. Similar but never-enrolled enrolled sites were identified by performing an Inverse Probability Weighting based on the estimated Propensity Score in respective samples. Tables S1 to S4 provide covariate balance tables before and after weighting in each of the region. In Fig. 2, we plot the annual average tree cover loss in enrolled sites and the annual weighted average in similar but never-enrolled sites to comparatively visualize tree cover loss before, during and after the intervention. The plots are simply illustrating annual tree cover loss of the enrolled and matched controls, so eventual post-2011 differences between the curves cannot be directly interpreted as an indication of a causal impact because the observed difference is not adjusted for the staggered enrollment into the program, or the eventual differences in tree cover loss at baseline, as well as do not display the statistical uncertainty around the differences between enrolled and never-enrolled sites.
Econometric identification strategy
Our main econometric estimates are based on a generalized Difference-in-Differences estimation (Doubly Robust Difference-in-Differences with multiple time periods) 59,60. Doubly Robust estimators combines a Difference-in-Differences design (to model the evolution of the outcome) with an Inverse Probability Weighting based on the estimated Propensity Score. This estimator is less susceptible of bias than canonical Difference-in-Differences or matching procedures because the Doubly Robust property ensure that the models are correctly specified even if either (but not both) the Propensity Score and the Difference-in-Differences outcome regression models are miss-specified.
Our generalized DID first dis-aggregates treatment effect into group-time Average Treatment Effects (ATT), i.e., the ATT of each group at each time period, where a group is defined by the year when each site is first enrolled in the program (Fig. S4) 61,62. Group-time ATT are obtained by performing a simple DID for each group (defined by year of first enrollment) and each time periods, and these group-time ATT can later be aggregated into a more precisely estimated cumulative ATT, corresponding to the sum of calendar year effects during program duration. The cumulative ATT of the intervention across all groups until the end of the program is the sum of the annual impacts of the program on tree cover loss (Fig. S7). Under the Parallel Trend Assumption that both enrolled and non-enrolled sites would have followed the same tree cover loss trend, conditional on observed confounders, we obtain an unbiased estimation of the impact of the program on the average accumulated avoided tree cover loss.
All our estimations were performed with R version 4.1.2 (https://cran.r-project.org/). We implemented our econometric estimates using the R package did, version 2.1.1 (https://bcallaway11.github.io/did/index.html), and R package DRDID version 1.0.3 (https://psantanna.com/DRDID/index.html). Robust analytical standard errors are clustered at the individual level in each specification in the main text.
We provide disaggregated figures of Average Treatment Effects in SI for our key region of interest (Settlements in Amazonia), namely a) Fig.S6: the complete group-time ATT (annual ATT for each pre and post initial enrollment effect, with respect to different year of first enrollment), b) Fig. S7 event-study graph (the annual impact of the program expressed for different duration before and after the first year of enrollment into the program, c) Fig. S8: Group ATT ( the aggregation of ATT with respect to groups defined as sites having the same year of initial enrollment, and d) Fig. S10: Calendar ATT, (the aggregation of ATT with respect to calendar year) (See also Supplementary Text). InS11 to S14, we also provide specification charts testing the sensitivity of our main estimates in all regions to variations in econometric specifications and subgroups obtained from the spec_chart R package (https://github.com/ArielOrtizBobea/spec_chart) 63,64. In Table S1 to Table S4, we also provide covariate balance tables (obtained from cobalt 4.3.2 R package https://ngreifer.github.io/cobalt) to visualize the performance of our initial Inverse Probability Weighting.
In Supplementary Text, Figs. S6 to S11, and table S1, we show complementary tests of the robustness of the estimated cumulative Average Treatment Effect. We conclude that our econometric specifications are robust to various set of covariates, but that variables characterizing pre-program forest cover loss play a preeminent role in our main estimate. We also observe important heterogeneity across years of initial enrollment into the program, as the average effect seems to be mostly driven by the cohort that started enrollment in 2012 (Figures S6, S8, and S9). We also assess that the program is not necessarily more effective in areas closer to roads or in sites where a higher proportion of households has been enrolled (Figure S11), although rigorously isolating the moderating effect of each variable on the effectiveness of BV would require further analysis.
Converting Average Treatment Effect into Standardized Forest cover effect sizes
Our econometric estimates correspond to the accumulated impact of the program on tree cover loss, expressed in percentage points of tree cover in 2000. To enable comparisons with evaluation performed in other contexts, we want to transform our statistical estimate into an interpretable Standardized Forest cover effect size (12). To do so, we multiply the cumulative ATT by the sum of tree cover in 2000 to obtain an estimate of the total avoided tree cover loss in hectare attributable to the program, under the assumption of homogeneous effect of the program across years and sites. Adapting the formula from 12, we note Fcounterfactual, 2017 = Fobserved, 2017 + (ATT * F observed, 2000) the counterfactual total forest cover at the end of the program in enrolled sites, and Fobserved, 2017 and F observed, 2000 respectively the total observed forest cover at the end of the program and in 2000 in these sites. ATT is negative when the program did avoid deforestation compared to the counterfactual scenario. The standardized forest cover effect size (i.e. the annualized effect of the program on forest cover change) corresponds to the quantity ra-rc = robserved -rcounterfactual = ln (Fobserved, 2017 / Fcounterfactual, 2017) / (2017-2010), as we use 2010 as pre-intervention baseline year, with approximate Standard Error = SE(ATT) / [Fcounterfactual, 2017 /(2017-2010)]. ra= robserved= ln (Fobserved, 2017 / Fobserved , 2010) / (2017-2010) and rc=rcounterfactual = ln (Fcounterfactual, 2017 / Fobserved , 2010) / (2017-2010) are respectively the annual observed deforestation and the counterfactual deforestation rate in enrolled sites, and Fobserved , 2010 is the total observed forest cover in 2010 in this sites.
Benefit-cost analysis
Our estimate of the avoided CO2 emissions is based on BV’s average deforestation impact. It does not take into account any avoided degradation, avoided emissions in non-forest woody vegetation, shifts in either belowground or soil carbon biomass, or new CO2 storage in tree regrowth. It assumes that when a tree is cut, all the CO2 eventually will be released to the atmosphere. We assume that CO2 emissions due to the program (induced tree cover loss or deforestation spatial spillovers or from gasoline for vehicles during program implementation) are negligible.
For cost data, we are using data on the annual expenditure of BV (not just the payment to enrolled households but the total program expenditures, including operational costs but not including non-budgeted costs eventually supported by enrollees and other stakeholders) in the Non-contributory Social Protection Programmes Database of the Economic Commission for Latin America and the Caribbean (https://dds.cepal.org/bpsnc/programme?id=60).
The estimation of the climate mitigation benefits of the program is obtained using the following formula from 16: Benefits per tons of avoided CO2 emissions = SCC* [1/ (1 + rS)] * [1 - [1 / ( (1 +r) D) ]]. Our estimate of climate-mitigation benefits is essentially a linear function of the number of years of delay for deforestation, and the slope is determined by various parameters linking forest cover to CO2 sequestration. Extrapolating from the Technical Support Document released by the US Interagency Working Group on Social Cost of Greenhouse Gases (February 2021), the Social Cost of Carbon (SCC) in 2016 with a 3% discount rate corresponds to USD 40.64 (SE=5.65) in constant US dollars 2016 (SCC= 61.88 with SE=7.77 when discount rate is 2.5%). We choose a discount rate of 3%, consistent with the suggestions of the United States Office of Management and Budget circular A-4 (2003), but also provide estimates for a discount rate of 2.5% (Figs. S15 and S16). We extrapolate the annual growth rate of the SCC from the evolution of this SCC between 2020 and 2050 and estimate the annual rate of growth to be about 2.19%. The effective discount rate r is then equal to (1.03/1.022)-1=0.079. Finally, we follow 16 in assuming that there is a D=10 years gap between tree-clearing and the complete release of CO2 into the atmosphere (45% of the biomass released immediately, 45% decomposing on average over 15 years, and 10% used as lumber allowing carbon to be stored for 30 years.). We have not tried to monetize other environmental (e.g., carbon stored in ground and belowground biomass but also other biodiversity and hydrological benefits) and social benefits (e.g., poverty alleviation) of the program, which could have further increased the cost-effectiveness of BV.
We use 2016 as the reference year in our calculation, i.e., all monetary values are expressed in constant 2016 USD dollars while the post-2016 annual costs and benefits of the program are discounted (using a 3% discount rate).
5 We did not include recognized traditional communities living on Brazilian public lands in our analysis because we were not able to find georeferenced information about these sites. Enrolled households living in these sites represented about 7% of the total number of households enrolled in the program.
6 In each tenure system, management plans notably state the rules, decision-making, and benefit-sharing processes related to collective and individual natural resources management. The compliance with the management plan and the degree of collective action can vary from one site to another 57.
7 The first row of Fig. S5 presents (statistically non-significant) findings at program scale, respectively in Settlements and Protected Areas, without taking into account systematic differences between Amazonian and non-Amazonian biomes.