Main

The UN Sustainable Development Goals underline the importance of economic growth and infrastructure development in alleviating poverty, while at the same time emphasizing that halting biodiversity loss is vital for global prosperity1,2. Policies aimed at delivering no net loss (NNL) of biodiversity, in theory, allow development to proceed while avoiding environmental damage3,4. NNL depends on implementation of the mitigation hierarchy: damage to biodiversity resulting from development must first be avoided, minimized and restored5, and any residual biodiversity loss offset through equivalent gains elsewhere6. At present, 101 countries either mandate some form of biodiversity compensation or support voluntary measures7. In countries with less-established environmental governance, lender requirements, such as the International Finance Corporation performance standards, are important drivers of NNL commitments8,9. Over 12,000 biodiversity offsets exist worldwide10, yet evaluations of their effectiveness are rare and most do not use robust methods11.

Offsets generate gains in biodiversity by creating or restoring habitat or protecting existing habitat which would have otherwise been lost (so called ‘averted loss’ offsets12). Offsets are controversial due to questions of permanence3, equivalence6 and equity13,14 and for generating gain against a background rate of biodiversity decline4,15. However, where high-quality habitat remains but is threatened by unregulated sectors, averted loss offsets may result in the best-possible biodiversity outcomes16. Biodiversity is an inherently complex concept so proxy measures are used to calculate losses and gains6. In forested ecosystems where most species are forest-dependent, forest loss can be a useful measure.

Quantifying the biodiversity gains from averted loss offsets requires estimation of the counterfactual scenario—the loss that would have occurred without protection15. While the counterfactual is inherently unknowable, statistical approaches exist to approximate it and consequently evaluate the impact of interventions on outcomes such as deforestation17,18,19. Statistical matching is commonly used to estimate the counterfactual on the basis of outcomes in matched control units, yet can be contingent on arbitrary modelling choices20. Recent advances, which test the robustness of estimates to a range of valid, alternative matching model specifications20 and different regression models18,21, can improve the quality of inference.

The Ambatovy nickel and cobalt mine (Fig. 1) is one of the largest lateritic nickel mines in the world. It is located within the biodiverse eastern rainforests of Madagascar which are highly threatened by deforestation, driven principally by shifting agriculture22,23. From the outset, Ambatovy promoted itself as a world-leader in sustainable mining and committed to ensure NNL, and preferably net gain, of biodiversity24,25. Its offset strategy was a pilot for the influential Business and Biodiversity Offset Programme24 which shaped guidelines widely used in mitigating biodiversity loss from development16,25. We use statistical matching and regression models to estimate the avoided deforestation achieved by Ambatovy’s four biodiversity offsets and check the robustness of our results to 116 alternative matching model specifications (Fig. 2). We provide encouraging evidence that this high-profile project, in one of the world’s hottest biodiversity hotspots, is on track to achieve NNL of forest and critically reflect on this finding in the broader context of NNL.

Fig. 1: Study area in eastern Madagascar showing the location of Ambatovy’s biodiversity offsets and our study design.
figure 1

a, The study area is the former province of Toamasina. Control pixels were sampled from pixels that were forested at baseline in 2000 (grey), excluding those within 10 km of a biodiversity offset or within established protected areas (grey hatched). The Corridor Ankeniheny-Zahamena (CAZ) new protected area was included in sampling (Methods). Inset area expanded in b. b, Ambatovy’s four biodiversity offsets: the Conservation Zone (yellow) which is within the mine concession area, the Corridor Forestier Analamay-Mantadia (CFAM; green), Torotorofotsy (blue) and Ankerana (orange). The 10 km buffer zone (which excludes established protected areas) around each offset is shown in lighter shades and was used to explore deforestation leakage. c, Our grid-based sampling strategy (Methods). The top layer illustrates the selection of our subsample of pixels. Data layers labelled X represent the outcome variable and covariates; all data used in this study are publicly available (Supplementary Table 4).

Fig. 2: Flowchart of methods.
figure 2

Statistical matching was used to match sampled pixels from each offset to control pixels sampled from the wider forested landscape with similar exposure to drivers of deforestation (Supplementary Table 4). Difference-in-differences regressions were run for each matched offset–control sample to estimate the effect of protection within each offset (termed site-based difference-in-differences). Pooled data were used in a fixed-effects panel regression to estimate the impact of protection across the whole offset portfolio. Resulting estimates were converted into hectares of avoided deforestation. a,b, To test the robustness of results to arbitrary modelling choices, the matching and outcome regressions were repeated using 116 alternative matching model specifications (a) to produce a range of estimates (b). The statistical distance measure used in matching (for example, Mahalanobis), caliper size, ratio of matched control to treated units and matching with or without replacement (shades of blue/purple) were varied in all 54 possible combinations. Holding these choices constant, we constructed 31 models comprising all possible combinations of five additional covariates (shown in shades of red/orange) with a core set of five essential covariates (green). Finally, we explore the robustness of the results to 31 randomly selected combinations of distance measure, model parameters and additional covariates.

Results

Ambatovy’s offset strategy is based on averted loss. It aims to generate biodiversity gains to offset the losses incurred at the mine site by preventing an equivalent amount of biodiversity loss within four biodiversity offset sites (which face a high rate of deforestation from shifting agriculture)24. To this end, the company and its non-governmental organization (NGO) partners implemented conservation activities aimed at slowing forest clearance within the four offsets. These included ecological monitoring, establishing community forest management associations and supporting them with the monitoring and enforcement of resource-use restrictions, environmental education programmes and promoting alternative income-generating activities in surrounding communities26,27. Occasionally the local police are brought in to assist with enforcement27.

According to our site-based difference-in-differences regressions (Methods), of the four biodiversity offsets associated with the Ambatovy mine, two significantly reduced deforestation relative to the counterfactual (Ankerana and the Conservation Zone ; P < 0.01). Protection reduced deforestation by an average of 96% (95% confidence interval (CI) 89–98%; P < 0.001, n = 38) per year in Ankerana and 66% (27–84%; P < 0.01, n = 38) per year in the Conservation Zone (Fig. 3 and Supplementary Table 9). One offset showed no significant effect (Torotorofotsy; −41 to +510%; P = 0.28, n = 38), while the remaining offset (Corridor Forestier Analamay-Mantadia (CFAM)) could not be assessed due to the lack of parallel trends in outcomes between the offset and matched control sample in the pre-intervention period—a critical assumption in difference-in-differences analyses. In CFAM, there was a significant declining trend in deforestation before protection while the matched control sample showed a significant increasing trend (Supplementary Fig. 5). Therefore, CFAM could not be used in the difference-in-differences analysis.

Fig. 3: The estimated percentage reduction in annual deforestation within each offset (from the site-based difference-in-differences regressions) and overall, across the entire offset portfolio (from the fixed-effects panel regression).
figure 3

The treatment effect is expressed as the average percentage difference in annual deforestation between the offset(s) and the estimated counterfactual following protection. Error bars represent 95% CIs (the upper bound for TTF extends to +510%). The width of the bar is proportional to the area of forest within each offset at the year of protection (Supplementary Table 2). ANK, Ankerana (orange); CZ, the Conservation Zone (yellow); TTF, Torotorofotsy (blue). CFAM (green) could not be included in the site-based difference-in-differences analysis due to lack of parallel trends in the pre-intervention period (Supplementary Fig. 5). n = 38 for Ankerana, the Conservation Zone and Torotorofotsy, and n = 152 for the overall result.

Including all four offsets in a single analysis using a fixed-effects panel regression (Methods), we estimate that protection reduced deforestation by an average of 58% per year (95% CI 37–73%, n = 152) across all four biodiversity offsets, relative to the estimated counterfactual (Fig. 3). We also tested the effect of excluding CFAM and estimate a greater reduction in deforestation of 72% per year (54–83%, n = 114; Supplementary Table 12 and Supplementary Fig. 8). Given that the two estimates are not significantly different (Z-test, P > 0.2), we present the more conservative estimate, which incorporates the effect of all four offsets, as our main result. Our results are also robust to the alternative specification of site and year as random effects (−53%, −27% to −69%; Supplementary Table 12).

Results robust to alternative model specifications

Arbitrary modelling choices, particularly associated with the decisions made in a matching analysis, are inevitable yet can exert a strong influence on estimated impacts28. Following Desbureaux20 we show that our results are robust to 116 alternative matching model specifications, all of which are a priori valid (Fig. 4). The vast majority of models for both Ankerana and the Conservation Zone confirm the results from the main model specification (Methods), presented in Fig. 3, of significant avoided deforestation. Where some models show an insignificant result (for example, for the Conservation Zone), in most cases these models are not a posteriori valid. By this we mean that >90% of treated units were unmatched (that is, a match within the caliper of the statistical distance measure could not be found), mean covariate balance exceeded the accepted threshold or parallel trends were not achieved. Exploring alternative model specifications also did not substantially change our results for Torotorofotsy; 78 of the 79 a posteriori valid models showed no significant impact of protection on deforestation, one suggested protection was associated with an increase in deforestation. For CFAM, the vast majority of alternative specifications, like our main model, were not a posteriori valid as they failed the parallel trends test. Of the seven a posteriori valid models, six showed no significant effect while one showed that protection was associated with a significant increase in deforestation relative to the counterfactual. Our result of a significant overall reduction in deforestation across all four offsets from the fixed-effects panel regression was robust for 106/116 alternative model specifications and none showed a significant increase in deforestation. Therefore, the evidence of avoided deforestation presented in Fig. 3 is robust.

Fig. 4: Raw estimates of treatment effect (points) and corresponding 95% CIs (bars) derived from 116 alternative matching model specifications.
figure 4

The alternative specifications included 54 possible combinations of matching distance measure and model parameters, 31 possible combinations of the 5 additional covariates with the core set of essential covariates, and 31 randomly selected combinations of distance measure, model parameters and additional covariates (Methods). Results from our main model specification, presented in Fig. 3, are shown in black. The asterisk indicates that the main model was not a posteriori valid. All alternative specifications are a priori valid but models that are not a posteriori valid (that is, >90% of treated units were unmatched, acceptable covariate balance or parallel trends were not achieved) are shown in lighter shades. See Supplementary Figs. 11 and 12 for full details of parameters and covariates associated with each result. Values are reported untransformed and represent the effect of treatment on the log(y + 1)-transformed count of annual deforestation.

We explored which modelling choices had the greatest influence on estimated impacts and found that the choice of statistical distance measure and model parameters had the most consistent, significant effect while the effect of including additional covariates is mixed (Supplementary Table 13).

NNL of forest nearly achieved by the offsets

The mine has destroyed or substantially degraded 2,064 ha of natural forest at the footprint and upper reaches of the slurry pipeline (henceforth mine site)24. The offsets have been in operation for between 7 and 12 years. Using site-based difference-in-differences regressions we estimate that between the year of protection and January 2020, 1,922 ha (95% CI 669–5,260 ha) of deforestation has been avoided within Ankerana and 26 ha (5–71 ha) has been avoided within the Conservation Zone (Fig. 5; Supplementary Methods). This equates to 1,948 ha of total avoided deforestation (over 94% of the forest loss caused by the mine), with most achieved in Ankerana. Using the fixed-effects panel regression incorporating all four offsets, we estimate an overall reduction in deforestation of 1,644 ha (674–3,122 ha) between 2009, when the first offset was protected, and January 2020 (Fig. 5). This represents over 79% (33–151%) of the forest loss caused by the mine. From 2014, when all the offsets became protected, an average of 265 ha of deforestation was avoided each year until 2020. If this rate continued, by the end of 2021, 2,174 ha of deforestation will have been avoided, fully offsetting forest loss at the mine site. Using the upper and lower bounds of estimated avoided deforestation (674 and 3,122 ha) suggests that NNL will be achieved between 2018 and 2033. In 2014, the company estimated that they would achieve NNL between 2022 and 203524. Our data therefore suggest that Ambatovy is on track to achieve NNL of forest earlier than anticipated.

Fig. 5: The total observed, counterfactual and the resulting estimate of avoided deforestation within each offset (estimated using site-based difference-in-differences regressions) and overall (using the fixed-effects panel regression) between the year of protection and January 2020.
figure 5

The counterfactual is an estimate of the deforestation that would have occurred in the absence of protection and was calculated using the estimated treatment effect (n = 38; Supplementary Methods). Avoided deforestation is the difference between the observed and counterfactual deforestation; negative values indicate that the offset resulted in a reduction in deforestation. The error bars show the 95% CI of the estimates of counterfactual deforestation (derived from the upper and lower CIs of the treatment effect) and the resulting estimates of avoided deforestation. The green dashed line indicates the 2,064 ha of forest loss caused by the mine itself. The number of years following protection is 9 yr for Ankerana, 11 yr for the Conservation Zone, 6 yr for Torotorofotsy and 11 yr overall (deforestation within later protected offsets is only counted from the year of protection).

Our estimate of the reduction in deforestation achieved within the Conservation Zone (26 ha, 1.6% of the total reduction in deforestation achieved within the offsets) is probably attributable to a combination of conservation management and the site’s location within the mining concession. The company and its predecessor (Phelps Dodge Madagascar) have been present in the concession area since the early 1990s, albeit with a hiatus from 1998 to 2003 (Supplementary Fig. 1). Therefore, for most of the 19-yr study period, access to the concession area, including the Conservation Zone, has been restricted27. This de facto protection reduced deforestation within the Conservation Zone to low levels before it was officially designated as an offset (Fig. 6).

Fig. 6: Comparison of the annual deforestation rate within the sample of pixels from each offset and the matched controls over the whole study period.
figure 6

The offset sample is shown in colour while the matched control sample is shown in grey. The dashed line indicates the year of protection. The offset and matched control samples contain an equal number of pixels (2,862 for Ankerana, 2,626 for CFAM, 1,340 for the Conservation Zone and 1,170 for Torotorofotsy) as the ratio of treated to control units in the matching was set to 1:1. For each offset, n = 38.

A number of studies have documented leakage effects from conservation interventions whereby impacts within the project area are simply displaced outside the boundaries, negating the effect of the intervention at the landscape-scale29. These leakage effects are not observed in our analysis of Ambatovy’s offsets (Supplementary Results) as we found that protection of the biodiversity offsets had no significant effect on deforestation within a 10 km radius (Supplementary Table 16; P = 0.15).

Putting these results in a broader context

Despite two-thirds of the 12,000+ biodiversity offsets that have been implemented worldwide occurring within forested ecosystems10, by 2019, <0.05% of these had been evaluated to assess the effectiveness of forest offsets at achieving NNL, and none of these evaluations used robust methods11 (although there have been several robust evaluations of wider offset policies12,30). This makes our estimation of the effectiveness of Ambatovy’s biodiversity offsets at avoiding deforestation valuable. Our results suggest that, by January 2020, the mine had offset 79% (33–151%) of the forest loss incurred at the mine site and is on track to achieve NNL by the end of 2021.

In recent years, there has been an explosion of studies using robust counterfactual methods to evaluate the effectiveness of other conservation interventions aimed at slowing tropical deforestation. Börner et al.19 synthesize these findings, using Cohen’s d normalized effect sizes to compare the effectiveness of 136 conservation interventions at reducing deforestation. Converting our estimate of the total avoided deforestation achieved by Ambatovy’s biodiversity offset policy (1,644 ha according to the fixed-effects model) to a Cohen’s d effect size yielded an estimate of −0.51 (classed as a ‘medium effect’31; Supplementary Results). This increases to −1.03 for the individual effect of Ankerana and −0.63 for the Conservation Zone (classed as ‘large effects’31). Comparison to the normalized effect sizes of the 136 other conservation interventions compiled by Börner et al.19 shows that overall Ambatovy’s biodiversity offsets were more effective at reducing deforestation than 97% of the other interventions and all bar one of the protected area interventions (Supplementary Fig. 10).

Discussion

We lack the empirical evidence to explain why Ambatovy’s offsets, as a whole, were so successful at reducing deforestation compared to other forest conservation interventions. We speculate that this may stem from the fact that offsetting is inherently centred on achieving measurable impact (NNL). All activities are designed specifically to meet this goal and progress can be regularly evaluated. Furthermore, large companies may possess sufficient funds to ensure, when they are committed, that they deliver this outcome. In contrast, public protected areas tend to be more focussed on measures such as coverage and investment and less explicitly impact-oriented32. Another important question is why conservation efforts were so successful in Ankerana but not in Torotorofotsy. It may be that enforcement of conservation restrictions was particularly effective within Ankerana, supported by evidence that local communities lost access to resources after the site was protected27 (discussed in more detail below).

Methodological caveats

An important caveat to our positive central result relates to the uncertainty inherent in impact evaluation using observational data33. The validity of causal inference rests on our ability to accurately model the counterfactual deforestation in the offset sites (what would have happened in the absence of the intervention) using data from matched pixels in the wider landscape which were not protected as offsets. In difference-in-differences analyses this assumes that all important factors influencing selection to treatment and the outcome of interest have been controlled for (or proxied) in the matching, so that the matched offset and control samples have similar trends in deforestation before the intervention33. Omitted variables may leave outstanding differences between the two samples, which can bias results33. Our choice of matching covariates is based on a good understanding of the local drivers of deforestation and selection to the treatment22,23 (Supplementary Methods) and our robustness checks demonstrate that our results are robust to alternative specifications (Fig. 4).

Our small sample size (n = 38 for the difference-in-differences regressions), limited by the length of the time series of the deforestation data34, reduces the precision of our estimates. In addition, methods for impact evaluation using observational data are constantly evolving with recent research highlighting the challenges of evaluating projects with staggered implementation dates35. Despite these caveats, which are the result of inherent challenges from such a real-world evaluation, our methodology represents a substantial advance in impact evaluation applied to biodiversity offsets. While our results seem relatively robust to alternative modelling specifications, this is only one case study. We hope that this work will stimulate further impact evaluations of biodiversity offsetting and emphasize the importance for future researchers to take considerable care over data selection and modelling choices (particularly the matching covariates, distance measure and model parameters) to ensure analyses are context-specific, appropriate and robust.

Wider concerns with offsetting

Biodiversity offsets in general, and averted loss offsets in particular, are controversial15,16,36. General criticisms include whether a concept as complex as biodiversity can be meaningfully reduced to proxies, questions of permanence3,37 and the potential social and equity issues of trading biodiversity (including access to ecosystem services) in one place for that in another13,14. Specific criticisms of averted loss offsets focus on the accuracy of counterfactual scenarios of loss against which gains are measured4,15 and the mismatch between stakeholder expectations and how much averted loss offsets can actually deliver16,36. We explore each of these criticisms in turn. In all cases they present clear and important caveats to our positive central result.

The aim of Ambatovy’s offset policy is to achieve NNL of biodiversity, whereas our study uses forest cover as an imperfect proxy. Rarely are the appropriate biodiversity data at the required spatial and temporal scale available to facilitate independent evaluation of NNL commitments. In forested ecosystems where most species are forest-dependent38, forest loss is a transparent and crucially measurable34 proxy for biodiversity loss. Furthermore, offsetting development-induced deforestation to achieve NNL of forest is a desirable outcome in itself, given its implications for biodiversity, ecosystem services and carbon storage. However, our measure of deforestation34 does not capture damage to forest biodiversity occurring at smaller scales, from activities such as selective logging, artisanal mining and harvesting of forest products for food, fuel and building materials39. Moreover, our method does not capture outcomes for species. In a context of high microendemism with many threatened species there is a real risk that large developments such as Ambatovy could lead to species extinction. To mitigate this risk, the company surveyed areas scheduled for clearance to identify, catch and relocate priority species to conservation areas outside the mine footprint (Supplementary Methods gives other mitigation measures) and conducted follow-up monitoring of certain species24. Whether the impacts of the mine on biodiversity are truly offset will depend on species responses to the changing pressures as well as the presence and efficacy of protection of these species within the offsets, which we were unable to capture in our analysis.

While we present strong evidence that Ambatovy has effectively conserved forest within its biodiversity offsets, questions remain regarding the permanence of this achievement. Although Ankerana and Torotorofotsy have been incorporated into the national protected area network and CFAM has been proposed as a new protected area26, continued effective management after the mine’s involvement ceases remains in doubt, given chronic underinvestment in Madagascar’s protected areas40. If the offsets become de facto unprotected after the company pulls out (expected between 2040 and 205024), deforestation is likely to resume and forest within the previously protected offsets may be lost. Offsets are intended to persist for as long as the impacts of the development remain3. Although Ambatovy have committed to restoring the impact site and have taken steps to prepare, tropical forest restoration is notoriously difficult41. If restoration fails and the offsets are no longer protected, a future acceleration in biodiversity loss will jeopardize Ambatovy’s claims to NNL.

Communities around Madagascar’s forests depend on forests for land to practice shifting agriculture and to provide wild products for food, fuel and building materials22,27. The mine and its associated biodiversity offsets have removed or reduced access to these provisioning ecosystem services. To compensate for this loss of access, Ambatovy invested in promoting alternative income-generating activities (including training and the provision of materials) in communities around the mine site and offsets26,27. However, research conducted within four affected communities (two near the Conservation Zone and two near Ankerana) found that local people did not consider these benefits to outweigh the substantial opportunity costs of the conservation restrictions27. The compensatory activities failed to reach those most affected by the restrictions and there was a temporal mismatch between the immediate loss of access to resources following establishment of the offsets and the time required for the alternatives to yield benefits27. This indicates that poor, rural communities living around the biodiversity offsets are bearing the cost of achieving NNL. For infrastructure developments such as Ambatovy to truly contribute towards sustainable development, SDG 1 (No Poverty) cannot be traded-off for SDG 15 (Life on Land). Instead, project proponents should strive to achieve NNL for both people and planet14.

An important criticism of averted loss offsets focuses on the accuracy of estimation of the counterfactual scenario; the baseline against which biodiversity losses and gains are measured4. Many offset policies use historical background rates of deforestation to define the counterfactual but previous studies have shown that this can overestimate the deforestation that would have occurred and consequently overstate the impact of the intervention17,37. We found that the baseline deforestation rates used by Ambatovy in their loss–gain calculations (based on the highest and lowest background deforestation rates at the district level between 1990 and 201024) are actually lower than the counterfactual rates we estimate here using robust methods for impact evaluation, meaning their estimates were conservative (Supplementary Table 1). However, there is an important caveat to this: the mine resulted in in-migration to the region26,27 which may have indirectly increased pressures on forest resources within the wider landscape, as observed with Rio Tinto’s QMM ilmenite mine in Southern Madagascar37. If any mine-related pressures were captured within the period used to define the ‘background’ rate of deforestation this would no longer represent baseline conditions in the absence of the mine and inflate the counterfactual (and the resulting estimates of biodiversity gains). Ambatovy employs ~9,000 people26, many of whom moved to the area from other regions of Madagascar26,27. The influx of migrant workers probably increased local demand for food, charcoal and fuelwood, which may have increased forest clearance and bushmeat hunting26,42. Such indirect impacts associated with industrial development are notoriously difficult to quantify and therefore offset43. Neither our approach, nor Ambatovy’s loss–gain calculations, could account for the indirect impacts of the mine on regional deforestation.

Another criticism of averted loss offsets is that they are premised on a background rate of biodiversity decline which can be slowed to generate the required biodiversity gains4,16. Therefore, even if NNL as defined by best-practice guidelines8 is achieved, loss of biodiversity has still occurred15,36. This is not what many stakeholders would understand as NNL of biodiversity44. However, given Madagascar’s high rates of deforestation45 and poor outcomes from tropical forest restoration41, averted loss is likely to be the better offsetting option16. Yet Madagascar has little remaining forest left to lose. Given the importance of the country’s biodiversity and the multitude of threats facing it40, future developments could aim to go beyond NNL and contribute towards the overall conservation of Madagascar’s remaining biodiversity16.

Hope for mitigating the environmental impacts of mines

There are over 6,000 industrial mines operating worldwide, covering an estimated 57,000 km2 (ref. 46) and impacting around 10% of global forested lands47. Low-income countries, like Madagascar, desperately need economic development. Mining, if well-regulated, can be part of the solution. From the start, Ambatovy promoted itself as a world-leader in sustainable mining and has some of the strongest commitments to conservation among 29 large-scale mines operating within forests47. Given this, and the resulting substantial investment the company made in NNL, failure would have been worrying for the concept of mitigating biodiversity loss from development. However, the achievements are notable, especially considering the challenging institutional and political context48 in which Ambatovy operates. Our results provide encouraging evidence that Ambatovy’s economic contributions to Madagascar49 (tens of millions of US dollars a year), were made while minimizing trade-offs with the island’s precious remaining forest habitat. There are many important caveats to this finding, as to any claim of NNL achieved through offsetting, however, the result certainly demonstrates the value of high aspirations combined with substantial investment in mitigating the biodiversity impacts of mining.

Methods

Study site and context

Ambatovy is a very large nickel, cobalt and ammonium sulphate mine in central-eastern Madagascar owned by a consortium of international mining companies50. It represents the largest ever foreign investment in the country24 (US$8 billion by 201650) and a substantial source of fiscal income49. In 2018, the company contributed ~US$50 million in taxes, tariffs, royalties and other payments49 and employed >9,000 people (93% of whom were Malagasy)51. Commercial production began in January 201424 (Supplementary Fig. 1). As key components in batteries, supplies of nickel and cobalt are critical to the green energy transition and demand for these metals is predicted to increase notably in future52.

The mining concession covers an area of 7,700 ha located in the eastern rainforests of Madagascar (Fig. 1) which have very high levels of biodiversity and endemism53,54. After avoidance and minimization measures were applied (Supplementary Methods) the mine was predicted to clear or substantially degrade 2,064 ha of high-quality natural forest at the mine footprint and upper pipeline24. Any impacts on plantations or secondary habitat are not included in this estimate. Losses at the impact site were not discounted in relation to a background rate of decline, meaning that the company took responsibility for the full area of forest lost25. Independent verification by our team (by measuring the size of the mine footprint on Google Earth) confirms the extent of forest loss at the mine footprint (Supplementary Fig. 2). Clearance of the footprint accounts for most of the forest loss associated with the mine as losses associated with the pipeline are small54.

Ambatovy aims to generate biodiversity gains to offset the mine-induced losses by slowing deforestation driven by shifting agriculture elsewhere26. To this end the company designated four sites, totalling 28,740 ha, to be protected as biodiversity offsets; Ankerana, Corridor Forestier Analamay-Mantadia (CFAM), the Conservation Zone and Torotorofotsy54 (Fig. 1). The offsets are considered like-for-like30 and were selected on the basis of similarity to the impact site in terms of forest structure and type, geology, climate and altitude24. The large combined area of the offsets relative to the impacted area was designed to allow flexibility, account for uncertainty and incorporate as many of the affected biodiversity components as possible24. Ankerana is the flagship offset, selected on the basis of its size, connectivity to the Corridor Ankeniheny-Zahamena (CAZ) forest corridor and the presence of ultramafic outcrops thought to support the same rare type of azonal forest lost at the mine site54. Extensive surveys conducted within Ankerana to establish biological similarity concluded the offset to be of higher conservation significance than the forests of the mine site due to the presence of rare lowland tropical forest24.

The Conservation Zone is directly managed by the company, given its location within the concession area, while the other offsets are managed in partnership with local and international NGOs24,25. Ambatovy funds the management of Ankerana by Conservation International and local NGO partners (although before 2015 Ankerana was directly managed by Ambatovy via a Memorandum of Understanding with Conservation International24), supports BirdLife partner Asity with the management of Torotorofotsy and a number of local NGOs including Voary Voakajy25 are involved in CFAM26. The company is also working to secure formal, legal protection for CFAM26 as part of a proposed Torotorofotsy–CFAM complex new protected area (although progress on this has stalled).

Overview of methods

To estimate the impact of the offsets on deforestation and determine whether this has prevented enough deforestation to offset forest loss at the mine site, we combined several complementary methods for robust impact evaluation. First, we used statistical matching to match a sample of pixels from each biodiversity offset to pixels from the wider forested landscape with similar exposure to drivers of deforestation. Then we used a site-based difference-in-differences regression for each matched offset–control sample and a fixed-effects panel regression on the pooled data, to estimate the effect of protection. We systematically explored how arbitrary modelling choices (including the statistical distance measure used in matching, caliper size, ratio of control to treated units, matching with or without replacement and which, if any, additional covariates were included) affected our inference, exploring the robustness of our results to 116 alternative model specifications.

Matching

The former province of Toamasina was selected as the geographic area from which control pixels were sampled as it encompasses forests of the same type as the concession area with varying degrees of intactness and accessibility. The four biodiversity offsets are located within this province (Fig. 1).

The unit of analysis is a 30 × 30 m2 pixel that was forested in the baseline year 200045,55. It is important that the scale of analysis aligns with the scale at which the drivers of deforestation (in this case, small-scale shifting agriculture) operate56. The median agricultural plot size (from 564 measured plots) in the study region is ~36 × 36 m2 (ref. 57). We took a subsample of pixels to reduce computational effort while maintaining the capacity for robust statistical inference58,59. We used a grid-based sampling strategy ensuring a minimum distance between sample units to reduce spatial autocorrelation60 and equal coverage of the study area58. A 150 × 150 m2 resolution grid, aligned to the other 30-m resolution data layers (Fig. 1c), was overlaid on the province and the 30 × 30 m2 pixel at the centre of each grid square was extracted to produce a subsample of pixels that are 120 m away from their nearest neighbour. The 120 m is larger than the minimum distance between units used in another matching study in Madagascar (68 m; ref. 59) but smaller than that used in other studies (200 m; ref. 61) and so strikes an appropriate balance between the avoidance of spatial autocorrelation and maximizing the possible sample cells.

Protected areas in the study area managed by Madagascar National Parks were excluded from our control sample as they are actively managed and therefore do not represent counterfactual outcomes for the biodiversity offsets in the absence of protection (Fig. 1). However, control pixels were sampled from within the CAZ new protected area as legal protection was only granted in 2015 and resources for management are limited and thinly spread62. Additionally, Ankerana and parts of CFAM overlap with the CAZ and would have experienced the same management, and likely trajectory, as the rest of the CAZ, had they not been designated biodiversity offsets. Areas within 10 km of an offset boundary were excluded from the control sample to reduce the chance of leakage (where pressures are displaced rather than avoided) biasing results17,29. The 10 km was selected as it is a commonly used buffer zone within the literature17,58.

To test for leakage effects, we used Veronoi polygons to partition the buffer area for CFAM, the Conservation Zone and Torotorofotsy (which overlap) into three individual buffer areas according to the nearest offset centroid and took a subsample of pixels from each (Fig. 1). Areas that overlapped with the established protected areas of Mantadia National Park and Analamazotra Special Reserve were excluded from the buffer zones.

The outcome variable is the annual deforestation rate sourced from the Global Forest Change (GFC) dataset34. Following Vieilledent et al.45 these data were restricted to only include pixels classed as forest in a forest cover map of Madagascar for the year 200045,55, reducing the probability of false positives (whereby tree loss is identified in pixels that were not forested). The resulting tree loss raster was snapped to the forest cover 2000 layer to align cells, resulting in a maximum spatial error of 15 m. The GFC product34 has been shown to perform reasonably well at detecting deforestation in humid tropical forests63. In the north-eastern rainforests of Madagascar, Burivalova et al.39 found that GFC data performed comparably to a local classification of very high resolution satellite imagery at detecting forest clearance for shifting agriculture (although it was not effective at detecting forest degradation from selective logging). As clearance for shifting agriculture is considered the principal agent of deforestation in the study area22 and the forests of the study area are tropical humid (>75% canopy cover), the GFC data are an appropriate tool for quantifying forest loss. Although recent evidence suggests that GFC data may have temporal biases64, this phenomenon affects our control and treated samples equally and so is unlikely to impact our results.

The choice of covariates is extremely important in matching analyses. They must include, or proxy, all important factors influencing selection to treatment and the outcome of interest so that the matched control sample is sufficiently similar to the treated sample in these characteristics to constitute a plausible counterfactual, otherwise the resulting estimates may not be valid33. On the basis of the literature and a local theory of change we selected five covariates that we believe capture or proxy for the aspects of accessibility, demand and agricultural suitability that drive deforestation in the study area22,59,65,66. These are slope, elevation, distance to main road, distance to forest edge and distance to deforestation (Supplementary Methods). These five essential covariates comprise the main matching specification and form the core set used in all alternative specifications that we tested in the robustness checks. We also defined five additional variables (annual precipitation, distance to river, distance to cart track, distance to settlement and population density) and tested the effect of including these in the robustness checks. The additional covariates were so defined because they were of poorer data quality (population density and distance to settlement), correlated with an essential variable (annual precipitation and population density) or simply considered less influential (distance to river and distance to cart track; Supplementary Methods).

Statistical matching was conducted in R statistics using the MatchIt package v.4.1 (ref. 67). To improve efficiency and produce closer matches we cleaned the data before matching to remove control units with values outside the calipers of the treated sample in any of the essential covariates (see Supplementary Methods for caliper definition). Following the recommendations of Schleicher et al.68 we tested several matching specifications and selected the one that maximized the trade-off between the number of treated units matched and the closeness of matches as the main specification (Supplementary Table 7). This was 1:1 nearest-neighbour matching without replacement, using Mahalanobis distance and a caliper of 1 s.d. This specification produced acceptable matches (within 1 s.d. of the Mahalanobis distance) for all treated units within all offsets. The maximum postmatching standardized difference in mean covariate values between treated and control samples was 0.05, well below the threshold of 0.25 considered to constitute an acceptable match69. This indicates that, on average, treated and control units were very well matched across all covariates.

Matching was run separately for each offset. The resulting matched datasets were aggregated by treated status (offset or control) and year to produce a matrix of the count of pixels that were deforested each year (2001–2019) in the offset and the matched control sample. Converting the outcome variable to a continuous measure of deforestation avoids the problem of attrition associated with binary measures of deforestation and is better suited to the framework of the subsequent regressions70.

Robustness checks

Statistical matching requires various choices to be made68, many of which are essentially arbitrary. There therefore exists a range of possible alternative specifications that are all a priori valid (although some may be better suited to the data and study objectives69) but which could influence the results20,28. We tested the robustness of our results to 116 different matching model specifications (Fig. 4). First, we tested the robustness of the estimates to the use of three alternative matching distance measures (Mahalanobis, standard propensity score matching using generalized linear model regressions with a logit distribution and propensity score matching using RandomForest), three different calipers (0.25, 0.5 and 1 s.d.), different ratios of control to treated units (one, five and ten nearest neighbours) and matching with/without replacement. Holding the choice of covariates constant (using only the essential covariates), the combination of these led to the estimation of 54 different models. Second, we tested the robustness of results to the inclusion of the five additional covariates. Holding the choice of distance measure and model parameters constant, we constructed 31 models comprising all possible combinations of additional covariates with the core set of essential covariates. Finally, we explore the robustness of results for 31 randomly selected combinations of distance measure, model parameters and additional covariates. All 116 specifications are a priori valid, assuming that the covariates capture or proxy for all important factors influencing outcomes, but may fail to satisfy the parallel trends condition or produce matches for insufficient numbers of treated observations (<10%), rendering them a posteriori invalid. It remains important to test the assumptions of the alternative models as failure to do so may lead to erroneous conclusions about effect size and direction being drawn from invalid models. Results are presented through specification graphs based on codes developed in Ortiz-Bobea et al.71.

Additionally, we tested the robustness of our results from the site-based difference-in-differences regressions to alternative temporal specifications using an equal number of years before and after the intervention (8 yr for Ankerana and the Conservation Zone, 6 yr for CFAM and 5 yr for Torotorofotsy) and dropping individual years from the analysis. This did not change the significance or magnitude of our results (Supplementary Table 10 and Supplementary Figs. 6 and 7).

Outcome regressions

Deriving estimates of causal effect from statistical comparisons of outcomes between treated and control samples relies on the assumption that the latter is a robust counterfactual for the former. In a difference-in-differences analysis this assumes that, in the absence of the intervention, the treated sample would have experienced the same average change in outcomes over the before–after period as the control sample72. Parallel trends in outcomes between treated and control before the intervention are an essential prerequisite for this assumption. We tested this for each matched offset–control dataset using the following formula:

$$\begin{array}{l}{{{\mathrm{log}}}}({\mathrm{Count}}\;{\mathrm{of}}\;{\mathrm{deforestation}} + 1)_{i,t} =\\ \quad \beta _0 + \beta _1{\mathrm{Year}}_t + \beta _2{\mathrm{CI}}_i + \beta _3({\mathrm{Year}} \times{\mathrm{CI}}_{it})+ \in _{i,t}\end{array}$$
(1)

where the outcome is the log(y + 1)-transformed count of deforestation within sample i at year t and CI is a binary variable indicating whether the observation is from the offset (1) or control (0) sample.

Parallel trends in deforestation between offset and matched control samples in the years before the intervention were present for all offsets except for CFAM (Supplementary Fig. 5). Consequently, CFAM could not be used in the site-based difference-in-differences analysis. However, its effect is still captured in the results from the fixed-effects panel regression as this is not based on an identifying assumption of parallel trends between groups in the pretreatment period72.

To estimate the impact of protection within each individual offset, we ran an ordinary-least-squares difference-in-differences regression for each matched offset–control dataset using the following formula:

$$\begin{array}{l}{{{\mathrm{log}}}}({\mathrm{Count}}\;{\mathrm{of}}\;{\mathrm{deforestation}} + 1)_{i,t} =\\ \beta _0 + \beta _1{\mathrm{BA}}_t + \beta _2{\mathrm{CI}}_i + \beta _3({\mathrm{BA}} \times{\mathrm{CI}})_{i,t}+ \in _{i,t}\end{array}$$
(2)

where BA and CI are binary variables indicating whether the observation occurred before (0) or after (1) the intervention, in the offset (1) or control sample (0). Given the non-normal properties of count data and the presence of zero values, a log(y + 1) transformation was applied to the outcome variable70,73. The coefficient of BA × CI and the corresponding confidence intervals were back-transformed (Supplementary Table 9) to obtain an estimate of the percentage difference in average annual deforestation between the offset and the matched control sample after protection, controlling for pre-intervention differences between samples (that is, the estimated counterfactual).

To estimate the overall impact of Ambatovy’s biodiversity offset policy at reducing deforestation we pooled the data for all four offsets and their corresponding matched control samples and ran a fixed-effects panel regression. The pooled data (n = 152) comprise an observation for each site (i = 8, 4 offset and 4 control) for each year (t = 19). The fixed-effects panel regression quantifies the effect of protection on the log-transformed count of deforestation controlling for site and year fixed effects, according to the following formula :

$${{{\mathrm{log}}}}( {\mathrm{Count}}\;{\mathrm{of}}\;{{\mathrm{deforestation}} + 1} )_{i,t} = \beta _0 + \beta _1{\mathrm{Tr}}_{i,t} + \propto _i + \gamma _t + {\it{\epsilon }}_{it}$$
(3)

where Tr is a binary measure indicating the treated status of sample i in year t (Tr = 1 for observations from offset sites in the years following protection and 0 for all other observations), i and γt represent site and year fixed effects, respectively, and it represents the composite error. The coefficient of interest (β1) and the associated confidence intervals were back-transformed to obtain the percentage difference in average annual deforestation across all four biodiversity offsets following protection (the treatment effect).

Evaluating deforestation leakage

To determine whether protection of the four biodiversity offsets simply displaced deforestation into the surrounding forested landscape, we repeated the matching and outcome regressions with the subsample of units from each buffer zone assigned as the treated group17,58 (Supplementary Results).

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.