Points, cells, or polygons? On the choice of spatial units in forest conservation policy impact evaluation

A fast-growing literature uses remotely sensed land-cover data along with quasi-experimental statistical methods to assess the efficacy of forest conservation interventions. A critical modeling choice is the spatial unit of analysis—points, grid cells, and polygons are all commonly used. Yet little is known about the implications of this choice for treatment effect estimates and for their interpretation. We demonstrate that point-level data can generate treatment effect estimates substantially different from those based on polygon-level data when (i) a disproportionate share of sample points is drawn from relatively large, treated polygons as a result of random or quasi-random spatial sampling, and (ii) the intervention analyzed has heterogeneous effects that depend on treatment polygon size. Our paper has four parts. First, using real-world data (on the award of timber extraction permits to forest management units in Mexico) that meet the two aforementioned criteria, we demonstrate that point- and polygon-level data generate qualitatively different results, and we propose a simple method for weighting the point-level data to recover the polygon-level results. Second, we conduct a Monte Carlo simulation to clarify the mechanism that causes this phenomenon and to provide reassurance that it is not driven by unobserved confounding factors. Third, we present new evidence (on Mesoamerican and Dominican protected areas) suggesting this phenomenon is not uncommon. Finally, we discuss the implications of our findings for the design and interpretation of spatial evaluations of forest conservation interventions. Although our analysis focuses on point- versus polygon-level data, the mechanism we describe also applies to grid cell- versus polygon-level data.


Introduction
In the two decades since the Millennium Ecosystem Assessment found that 'few well-designed empirical analyses assess even the most common biodiversity conservation measures' (MEA 2005) and social and natural scientists called for research to fill that gap (Sutherland et al 2004, Ferraro andPattanayak 2006), the evidence base on the efficacy of conservation policies has grown rapidly (Börner et al 2020, Busch andFerretti-Gallon 2023).A mainstay of this literature is so-called spatial evaluation, which entails using remote sensing data on land-cover change along with quasi-experimental statistical methods to assess the efficacy of interventions such as protected areas and community forestry that 'treat' some areas but not others (Blackman 2013).Impacts are measured by comparing tree cover changes in treated and untreated (control) areas, controlling for confounding factors.
A critical modeling choice in such analyses is the spatial unit of analysis.Points (defined by latitudelongitude coordinates; e.g.Baragwanath andBayi 2020, Vélez et al 2020), polygons (administrative units such as communal properties; e.g.Sims andAlix-Garcia 2017, Blackman andVillalobos 2021), and cells (from a user-defined grid; e.g.Somanathan et al 2009, Ben Yishay et al 2017) are all commonly used.The choice is often driven by data availability: points and cells are often used because polygon-level data are not defined for untreated (control) units.But does this choice have implications for treatment effect estimates and for the interpretation of these estimates?
To our knowledge, the only analysis focused squarely on these questions is Avelino et al (2016) 3 .Using both real-world and simulated data (i.e. a Monte Carlos analysis), the authors evaluate the effect of a forest conservation intervention using data aggregated into grid cells of varying sizes.They find that using cells that are either too large or too small relative to the spatial unit at which land use and land cover change decisions are made can result in inefficient and/or biased treatment effect estimates.On one hand, using cells that are too small can introduce unnecessary noise into the analysis, thereby biasing treatment effect estimates toward zero-an errorsin-variables problem.An extreme example would be using cells so small as to pick up meter-by-meter variation in slope that has minimal influence on actual land use and land cover decisions.On the other hand, using cells that are too large relative to a binary treatment can also generate an errors-in-variables problem.An example would be an analysis of the effect of protected areas that used cells so large as to encompass both areas that were legally protected (treated) and areas that were not protected.Using cells that are too large can also result in a loss of efficiency by restricting sample size.
Here, we explore an additional way in which the choice of spatial unit could influence treatment effect estimates for conservation interventions.We show that point-level data can generate estimates substantially different from those based on polygon-level data-that is, estimates that differ qualitatively as well as quantitatively-when the intervention analyzed has heterogeneous effects that depend on polygon size.The reason is that in such cases, statistical models used to generate treatment effect estimates using point-level data will give greater weight to the effects of the intervention in relatively large polygons simply because any random or quasi-random selection of sample points will draw a disproportionate number of points from such polygons.By contrast, treatment effects generated using polygon-level data will weight effects on all sample polygons equally.Because this phenomenon can occur whenever the representation of units of analysis in the study sample depends on the size of the treated polygons (i.e., whenever the sample includes a disproportionate number of units from large polygons), it can also arise when cells are used as spatial units4 .We propose a simple weighting method that can be used to control for this type of selection and render a point-or cell-level spatial analysis more comparable to a polygon-level analysis.We believe these findings have important implications for the design of spatial evaluations and for the interpretation of their results.
Our paper has four parts.First, we use real-world data to compare quasi-experimental treatment effect estimates generated using polygon-level versus pointlevel data.The treatment we analyze with both data sets is the award of timber extraction permits to communal forest management units (FMUs) in Mexico, an intervention that its advocates say helps stem forest loss by strengthening local agents' conservation incentives.Using polygon (FMU)-level data, we we are not able to discern a statistically significant effect on tree cover loss.But using point-level data, we find that they have a statistically and economically significant negative effect.In addition, we find that the size of the polygons (FMUs) from which points are drawn moderates the effect of permits: permits do not have a discernible effect in small polygons but reduce tree cover loss in large ones.
Second, we conduct a Monte Carlo simulation to clarify the hypothesized mechanism for the phenomenon described above and to provide reassurance that our case study results are not instead driven by unobserved confounding factors.We demonstrate that it is possible to replicate the qualitative results from our case study using simulated data in which confounding factors are not present.
An implication of the first two parts of the paper is that the phenomenon we describe here is not confined to our case study.In principle, it should arise in either point-or grid cell-level spatial analyses whenever the two factors noted abovedisproportionate representation of spatial units from large treatment polygons and effect heterogeneity across treatment polygons of different sizes-are present.The third part of the paper explores evidence on whether this phenomenon is, in fact, common.Relying on a recently published systematic review of spatial evaluations of conservation interventions (Busch and Ferretti-Gallon 2023), we show that unfortunately, we have virtually no evidence one way or the other, simply because researchers have not conducted the requisite heterogeneity tests.To help fill that gap, we test for and find treatment effect heterogeneity across polygons of different sizes using pointlevel data on national protected areas in Mesoamerica and the Dominican Republic.
The fourth and final part of the paper explores these implications of our findings for spatial analysis of conservation and development interventions.Replication data for the first three parts of the paper are available online (Blackman et al 2024).

Case study
We test the increasingly common Coasian argument that awarding timber extraction permits to local FMUs can help stem tree cover loss by giving managers of these FMUs economic incentives to protect forests so as to ensure a long-run stream of future profits from timber sales (Griscom et al 2017, Runting et al 2019) 5 .Our case study focuses on the award of permits to communal FMUs in Mexico (called ejidos and comunidades) that manage an estimated 55%-80% of the country's forestland (Bray andMerino-Pérez 2002, Madrid et al 2010).Blackman and Villalobos (2021) present an FMUlevel analysis of this intervention.Here, we present a complementary point-level analysis and compare the two sets of results.We only briefly sketch the background for the analysis and the broad empirical approach since these topics are covered in detail in Blackman and Villalobos (2021).The brief summaries in the next three sections are drawn from that article.For the remainder of this paper, we refer to Mexican communal FMUs as polygons in order to link the analysis to the more general issue of the choice of spatial unit.

Empirical approach
To identify the effect of timber extraction permits on tree cover loss, both our polygon-level and pointlevel models use 2001-2012 panel data along with a two-stage empirical strategy.In the first stage, we use propensity score matching to preprocess our data, a step that helps control for misspecification and omitted variables bias.This procedure entails dropping from our sample all never-permitted polygons that are not similar to ever-permitted polygons in terms of time-invariant observable characteristics that affect tree cover loss.In the second stage, we use the preprocessed data to estimate a two-way fixed effects model, specified as where u indexes the type of spatial unit (either points, p, or polygons, f ), i indexes spatial units, t indexes years, Y is either the percentage of the polygon with tree cover loss (u = f ) or a binary indicator of such loss (u = p), γ are spatial unit fixed effects, δ are year fixed effects, D is a binary dummy variable indicating that a permit was in force in year t, X is a vector of time-varying control variables, βs are parameters or vectors of parameters to be estimated, and ε is an error term.The parameter β 1 measures permits' net effect on tree cover loss-formally, the average treatment effect on the treated.We estimate equation (1) using ordinary least squares and we cluster standard errors at the spatial unit level.The principal identifying assumption is that absent permits, trends in tree cover loss would be the same for both the treatment and the control groups.The supplementary material presents tests of this assumption (table S5, figure S1).

Sample
For our polygon-level data set, the universe of all possible observations from which our regression sample is drawn is the set of all polygons in national cadastral data (RAN 2010).For our point-level data set, the universe is a set of points selected by overlaying on a map of Mexico a grid with horizontal and vertical gridlines 350 m apart and selecting points where gridlines cross.In both our polygon-level and pointlevel data sets, we limit the sample to spatial units in the 16 Mexican states that have significant tree cover (figure 1)6 .In addition, we limit these samples to polygons and points with communal tenure.Finally, we limit our sample of points to those with tree cover at baseline (i.e. in year 2000).
We use a permit database (INECC 2013) to identify treated communal polygons in our 16-state study area-that is, polygons that had timber extraction permits during our 2001-2012 study period.See Blackman and Villalobos (2021) for a discussion of the criteria used to define our samples of polygons and points.
Our polygon-level regression sample comprises 9527 polygons, of which 650 are ever-permitted and 8877 are never-permitted (figure 1).From these cross-sectional data, we create a 12-year unbalanced panel spanning 2001-2012.The panel covers 113 805 polygon-years, 7301 of which are treated (under permit) and 106 504 are untreated.
Our point-level regression sample comprises 1585 104 points, of which 340 413 are ever-permitted and 1246 965 are never-permitted.The pointlevel panel data set contains 18 678 729 pointyears, of which 3715 257 are treated and 14 963 472 are untreated.Figure 2 illustrates the relationship between our point-level and polygon-level data sets.
Note that the distribution of the size of the 9527 polygons represented in both our polygon-level and our point-level data sets is decidedly skewed.The majority of polygons are smaller than 5000 ha but a significant number are much larger (figure 3, panel (B)).

Variables
Table 1 lists the variables in our empirical analysis, including their names, units, definitions, sources, spatial scales, years, and for each data set, means and standard deviations.The time-varying variables are included directly in our two-way fixed effects regressions, whereas the time-invariant variables (with the exception of large and small) are used to generate propensity scores, which in turn are used to match ever-treated and never-treated spatial units7 .Typically, the number of external agents who have some type of de facto land rights.
e Scale = 1:50 000-1:1000 000.Variable means in the polygon-level data set and point-level data set differ because a disproportionate share of observations in the point-level data set comes from large polygons.For example, the mean of size in the polygon-level data is 3370 ha versus 27 710 ha in the point-level data set, an eight-fold difference.

Results
Because the focus of our analysis is comparing treatment effects generated using our polygon-level and point-level data sets, and not interpreting the effects themselves, we relegate to the supplementary material most of the intermediate statistics that underpin these estimates, including propensity score results, balance checks, and parallel trends tests.
Consistent with the broad qualitative result in Blackman and Villalobos (2021), the treatment effect estimated using our polygon-level data is not statistically significant (table 2, Model 1).The implication is that the award of a permit does not have a discernible effect on tree cover loss.
By contrast, the treatment effect estimated using our point-level data is negative and significant at the one percent level (table 2, Model 2).The estimated coefficient (−0.0004) implies that the award of a permit reduces the average annual probability that a sample point is cleared by four one-hundredths of a percentage point, a 17% reduction compared with the counterfactual probability of clearing (0.0024).
As noted above, we hypothesize that the treatment effect estimates generated by the two data sets differ because (i) the point-level data set contains a disproportionate share of observations drawn from large polygons, and (ii) large and small polygons have differential effects on tree cover loss.Two additional sets of econometric results support that hypothesis.
First, if we weight all treated observations in the point-level data so that collectively, all points from any given ever-treated polygon have the same weight as all points from any other ever-treated polygon, then we can recover the null result from the polygonlevel regression (table 2, Model 3).Specifically, we obtain this result if we assign each ever-treated observation a weight equal to the total number of points in all polygons divided by the number of points in the polygon with which that observation is associated.That is, we recover the polygon-level result by estimating where and w is the sampling weight, N is the total number of points in all polygons, and n is the number of points in the polygon where the point is located.Second, tests indicate that polygon size moderates the effect of permits on tree cover loss.To test for moderating effects, we include in our model two interaction terms: one that interacts permit with large, an indicator of whether an ever-treated point is drawn from a polygon larger than median size, and a second that interacts permit with small, an indicator of whether the point is drawn from a polygon smaller than median size.We find that the first interaction term is negative and significant and the second is not significant (table 2, Model 4).These results suggest that permits reduce tree cover loss in large polygons but do not have a discernible effect in small ones.In addition, expressed as a function of polygon size, the marginal effect of a permit on tree cover loss is negative and significant for polygons smaller than about 70 000 hectares, a range that includes more than 99% of all polygons and 89% of all points in our sample (figure 3, panel (A)).Moreover, the marginal effect is decreasing in polygon size.These results are consistent with our hypothesis that the negative effect of the treatment on the probability of clearing in Model 2 is driven by the observations from relatively large polygons 8 . 8Blackman and Villalobos (2021) include a theory of change that describes mechanisms having to do with (i) FMU governance capacity; (ii) the demand for cleared land; and (iii) external monitoring and enforcement of land use and land cover regulations.The theory of change hypothesizes that the award of permits increases tree cover loss in FMUs where governance and/or external monitoring are weak and/or where demand for cleared land is strong; conversely, permits decrease tree cover loss where the opposite is true.Any or all of these three mechanisms could explain the differential effects of permits in large and small FMUs if they were correlated with FMU size-that is, if governance and/or external monitoring were weak and demand for cleared land was strong in small FMUs, and vice versa.Although Blackman and Villalobos (2021) do not test for these correlations, subsequent analysis of the data used for that study confirms the first two correlations: in small FMUs, governance (proxied by a measure of socioeconomic level) is relatively

Monte Carlo simulation
This section presents a Monte Carlo simulation that clarifies the hypothesized mechanism for the phenomenon we have described.In addition, by demonstrating that it is possible to replicate the qualitative results from our case study using simulated data in which unobserved confounding factors are not present, it provides reassurance that these results are not driven by such factors.In what follows, we present an overview of the Monte Carlo simulation and its main results.In the interest of concise exposition, details are relegated to the supplementary material.
The Monte Carlo simulations has two stages.In the first stage, we define two data-generating weak, and demand for cleared land (proxied by a measure of the average return to agriculture) is strong.The data do not provide support one way or the other for a correlation between external monitoring and FMU size simply because they do not include a measure of (or proxy for) external monitoring.processes (DGPs) and use them to create thousands of simulated data sets.The DGPs and the simulated data sets are constructed so as to reflect the 'true' relationships among tree cover loss (our dependent variable), permits (our independent variable of interest), and polygon size (a moderator), which together constitute the mechanism we believe explains our real-world results.More specifically, we first define a DGP to create simulated polygon-level data in which permits have heterogeneous effects  depending on polygon size.Next, we define a second DGP to generate simulated point-level data so as to mimic the process of using a quasi-random sampling grid to select points within treatment and control polygons.Both DGPs are constructed to generate simulated data with the same basic structure as our real-world data.For example, these data sets feature the same number of time steps (12 years), approximately the same number of polygons (1000), and roughly the same percentage of treated versus control polygons (55%).We use the two DGPs to generate 1000 polygon-level data sets and 1000 point-level data sets.
In the second stage of the Monte Carlo simulation, for each of our simulated data sets, we estimate a treatment effect using a two-way fixed effects estimator at either the polygon level (n = 1000) or the point level (n = 1000), similar to equation (1).
Histograms and summary statistics provide a first indication that mean values of these estimated treatment effects comport with the broad qualitative results from our real-world analysis (figures 4 and 5).For the polygon-level data, the distribution of treatment effects is centered around zero and the mean is not statistically different from zero.For the point-level data, the distribution is centered around −0.004 and the mean is significantly different from zero.And finally, for the weighted point-level data, the distribution is once again centered around zero and the mean is not statistically different from zero.Tests for differences in means confirm that mean treatment effect estimates for the polygon-level and point-level data are different but estimates from the polygon-level data and weighted point-level data are not (table 3).

How common is this phenomenon?
The phenomenon described above arises in the context of a spatial evaluation of a conservation intervention when the following two conditions are met: i. the spatial unit of analysis is either grid cells or randomly or quasi-randomly selected points such that the sample contains a disproportionate share of units drawn from relatively large, treated polygons, and ii. the intervention analyzed has heterogeneous effects that depend on the size of the polygons that map affected areas (e.g.protected areas, community forests, payments for environmental services parcels).
How often are these criteria met?Unfortunately, we have very little evidence one way or the other simply because researchers have not conducted the heterogeneity tests needed to establish criterion (ii).That conclusion is based on a review of the compiled for Busch and Ferretti-Gallon (2023), a systematic meta-analysis of peer-reviewed spatial evaluations of conservation and development interventions9 .Of the 320 articles included in that review, 155 use grid cells or points as the spatial unit of analysis (criterion i).But of these 155 articles, none use interaction terms to test whether the treatment effect is moderated by the size of treatment polygons (criterion ii).
Given the scarcity of evidence about the prevalence of the problem described above, we searched for new evidence using point-level data on national protected areas (NPAs) in Mesoamerica and the Dominican Republic10 .By definition, these data meet criteria (i) listed above.Econometric analysis indicates that they also meet criteria (ii): NPAs have heterogeneous effects that depend on their size (for details, please see the supplementary material).Note that here, however, the moderating effect of treatment polygon size is the opposite from our Mexican case study: the effect of the treatment is stronger in small polygons, not large ones (table S9).Nevertheless, the general phenomenon we have described is the same: average treatment effect estimates using point-level data disproportionately reflect differential effects that occur in large treatment polygons.Overall, this analysis suggests that this phenomenon is not unique to the case study described in section 2.

Discussion
What are the implications of our findings for spatial analysis of conservation interventions?We offer five brief observations.First, for the reasons discussed in section 4, we believe our findings are applicable beyond our Mexico case study.
Second, the fact that point-and cell-level data can generate qualitatively different treatment effect estimates does not necessarily imply that one estimate or the other is 'correct.'Rather, it simply implies that they need to be interpreted differently.A treatment effect estimate from a point-or grid cell-level analysis should be interpreted as the per-unit-of-land effect of the treatment, whereas a treatment effect estimate from a polygon-level analysis should be interpreted as the per-polygon effect.Either one may be appropriate, depending on the goal of the analysis.If the goal is to evaluate the average effect of the treatment on all units of treated land-perhaps to gauge the total ecological effect of the program-then a point-or celllevel analysis would be called for.But if, on the other hand, the goal is to shed light on the average effect of an intervention on a polygon-perhaps to assess the expected marginal effect of expanding or continuing the intervention-then a polygon-level analysis would be better.To make this distinction more concrete, we reference our case study.If our goal was to evaluate ex post the average effect of awarding timber harvest permits to communal polygons on all of the five million hectares of treated land represented in our study sample, then using (unweighted) pointlevel data would be appropriate.But if, on the other hand, our goal was to understand the average effect on a polygon, then using polygon-level data would be appropriate.
Third, by using point-or cell-level analysis along with weighting, it is possible derive both types of estimates.Unweighted data can provide a per-unitof-land treatment effect estimate, and weighted data can provide a per-polygon estimate.
Fourth, whatever spatial unit of analysis is chosen, it is important both to make explicit the appropriate interpretation and, in the case of point-level analyses, to explore and test for possible heterogeneous treatment effects for large and small treatment polygons.For example, referring once again to our case study, in an analysis that relied on point-level data, it would be important to clarify that treatment effect estimates are per-unit-of-land; that the award of permits has different effects on large and small polygons; and that, as a result, the per-polygon effect might well be substantially different.
Finally, we recognize that a variety of considerations influence the choice of spatial unit of analysis.These include, for example, the availability of cadastral data for both treated and untreated polygons (which is a requisite of polygon-level analysis); the spatial scale at which the intervention being analyzed is expected to affect outcomes according to the theory of change that underpins the evaluation; and the challenge of developing polygon-level dependent and independent variables by aggregating data at a finer spatial scale (Avelino et al 2016. Blackman 2013).We believe that the additional factors elucidated here should be added to this list of considerations.

Figure 1 .
Figure 1.Study area and regression sample.

Figure 2 .
Figure 2. Representative points and polygons in Oaxaca state.

Figure 3 .
Figure 3.Estimated marginal effect of permits on tree cover loss using point-level data as a function of polygon size (panel (A)); and distribution of points by size of polygons from which they are drawn (panel (B)).

Figure 4 .
Figure 4. Effect of timber extraction permits on tree cover loss in Mexico; distribution of treatment effects estimated using simulated polygon-level data (panel (A)); simulated point-level data (panel (B)); and weighted simulated point-level data (panel (C)).

Figure
FigureEffect of timber extraction permits on tree cover loss in Mexico; mean estimated treatment effects using simulated polygon-level data, simulated point-level data, and simulated weighted point-level data.Dots are point estimates, whiskers are 95% confidence intervals.

Table 1 .
Effect of timber extraction permits on tree cover loss in Mexico; variables and sample means.
a Original scale = 0-6, rescaled by adding 2.2 to all values.b Ejido or comunidad members with formal land rights.c Adults who are formally recognized as residents of the ejido/comunidad but who do not have formal land rights.d

Table 2 .
Effect of timber extraction permits on tree cover loss in Mexico; treatment effect estimates from two-way fixed effects regressions; dependent variable is tree cover loss (s.e.).covariates are temperature and rainfall.Time-invariant covariates used to select matched control variables are listed in table 1.Standard errors are clustered at the level of the spatial unit of analysis.a Predicted outcome with all treatment variables is set equal to zero, computed using delta-method.

Table 3 .
Effect of timber extraction permits on tree cover loss in Mexico; tests for difference in mean estimated treatment effect using simulated data.