How Differently Do Farms Respond to Agri-environmental Policies? A Probabilistic Machine-Learning Approach

Observational Dataset

We use information from the Italian Farm Accountancy Data Network (FADN), which represents the only source of microeconomic agricultural data that is harmonized at EU level and collects physical, structural, economic, and financial data on farms in all EU member states (European Council 2009). The survey is representative of the farms that can be considered professional and market oriented, due to their economic size (that is equal or more than €8,000 of standard output). In Italy these correspond to 95% of utilized agricultural area, 97% of the value of standard production, 92% of labor units, and 91% of livestock units. The representativeness of the dataset is ensured on three dimensions, namely region, economic size, and farm typology. For these reasons, the FADN is the most (and only) widely used farm-level dataset for, among others, CAP evaluations and specifically for the assessments of the AEP impacts (among others, Arata and Sckokai 2016; Bartolini et al. 2021; Stetter, Mennig, and Sauer 2022).

Our research focuses on the 2014–2020 programming period of the CAP.12 However, unlike Stetter, Mennig, and Sauer (2022), we exclude the initial year (2014) for two reasons: first, payments of one of the policies under consideration (the GP) only started in 2015; second, many of the farms observed in 2014 may still benefit from measures of the previous programming period. We thus focus on the 2015–2020 period, although we only have detailed and validated information until 2018. Therefore, our initial sample consists of a representative collection of Italian commercial farms that produces an unbalanced panel consisting of 9,580, 10,135, 10,792, and 10,386 observations in 2015, 2016, 2017, and 2018, respectively. Because our analysis does not address regime-switching dynamics, we only consider farms for which the treatment status did not change over the period analyzed; that is, T_it,k = T_i,k for all i ∈ (1,…, N). For this reason, we first extract a balanced panel consisting of 5,836 units observed over 2015–2018 and then drop all entries satisfying T_it,k ≠ T_is,k for any s, t ∈ {2015, …, 2018} and s ≠ t.13 The resulting dataset consists of 4,001 farms repeated over four years, for a total of 16,004 observations. Compared with other related works (Bertoni et al. 2020; Stetter, Mennig, and Sauer 2022), our study provides wide coverage of the agricultural sector by focusing on the entire national area instead of a single region. Furthermore, since the treatments presented in Section 4 are likely to affect the agri-environment over several years, our outcome variable uses information from the last two years in the series to account for potential accumulation effects (see Section 4 for details).

Definition of Treatments

As mentioned in Section 2, the 2015–2020 CAP AEP design is based primarily on two main policy instruments that belong to CAP’s pillar 1, pillar 2, or both. On the one hand, we observe pillar 1 subsidies that are conditional on a set of compulsory requirements (i.e., CC and the GP) with which farmers must comply to preserve the DP. On the other hand, we have voluntary measures aimed at compensating farmers for income losses or increased costs resulting from the voluntary adoption of more sustainable farming practices (i.e., the AEM of pillar 2). Consequently, farms are subscribed to—in fact, they voluntarily choose—one of three possible policy alternatives, which effectively reflect the interplay between the two pillars of the CAP: (1) farms failing to meet all the CC and GP requirements, that is, farms receiving neither pillar 1 nor pillar 2 payments; (2) farmers receiving both pillar 1 (DP and GP) and pillar 2 (AEM) payments; and (3) farms complying with the CC and GP requirements but not adopting any AEM.

Table 1 indicates how the farms in our sample are distributed across the three policy categories. The third cohort is the largest group, which includes approximately 71% of the observed farms (2,841 units). Using the terminology introduced in Section 3, we consider the corresponding policy option as the baseline treatment, T_i,0, associated with the netput vector y_it,0. Next, all farms choosing not to benefit from pillar 1 and pillar 2 payments (i.e., the first cohort, corresponding to approximately 13% of the sample) take up the first treatment, T_i,k=1, which implies giving up both pillar 1 and pillar 2 resources. We assume that this decision follows the behavioral model stylized in Section 3, according to which, conditional on X_it, T_i,1 produces higher profits than T_i,0. Similarly, farms applying for pillar 2 AEM supports (i.e., the second cohort, corresponding to approximately 16% of the sample) choose treatment T_i,k=2 through the same profit-maximizing mechanism. In this respect, our work extends the analysis in Stetter, Mennig, and Sauer (2022) by distinguishing between the two different AEPs (i.e., the AEMs and the pillar 1 environmental requirements).

We postulate that treatments T_i,1 and T_i,2 belong to two nonoverlapping choice sets; in other words, we rule out a multiple treatment setup by positing treatment T_i,1 as infeasible for farms choosing T_i,2 and vice versa. Although this assumption is quite strong, it is necessary to identify the treatment effects of interest. Given that T_i,1 and T_i,2 represent two ends of a rather wide spectrum of policy options, it is plausible that both treatments may appeal to (i.e., are feasible for) farms with very distinctive characteristics. Conversely, our setup implies that both T_i,1 and T_i,2 are feasible alternatives to the baseline treatment T_i,0. This presupposes that farms in the control group are characterized by features X_it that overlap with the characteristic of the units in T_i,1 or T_i,2. That is, we can always find comparable farms in either of the two groups in different strata of X_it, that is, 0 < Pr(T_i,k = 1 | X_it = x_it) < 1. This restriction is also commonly known as common support (or positivity), and as we discuss in Section 5 and Appendix E, it limits extrapolation issues, thus preventing unreliable treatment effects.

One caveat in our setup is that unlike T_i,1 farms choosing T_i,2 may in fact opt for one among four treatment versions. As outlined in Section 2, T_i,2 aggregates measure 10 and 11 which in turn can be decomposed in two submeasures: agri-environment-climate commitments (10.1); conservation and sustainable use and development of genetic resources in agriculture (10.2); payment to convert to organic farming practices and methods (11.1); and payment to maintain organic farming practices and methods (11.2). While measure 10.2 only concerns a small share of farms (roughly 3% of our sample) and can be thus excluded or safely merged into measure 10.1 (our current choice), submeasures 11.1 and 11.2 are substantially equivalent in terms of farmer behavior, the only difference being the amount of support granted. For this reason, we de facto consider submeasure 11.1 and 11.2 as a unique measure (i.e., measure 11). As put forward in Sections 2 and 3, disregarding such distinctions may greatly affect the interpretation of the HTEs via treatment heterogeneity.

It is also worth mentioning that in principle, the submeasures could be further disaggregated into specific actions (using the RDP jargon). Unfortunately, the Italian FADN data do not provide enough information on AEM actions. In fact, to our knowledge, there are no high-quality representative datasets that can provide more detail on AEMs (e.g., Stetter, Mennig, and Sauer 2022, who use the German version of our dataset). Had this level of disaggregation been observable, it would imply a very large number of actions (i.e., treatment versions), as evidenced by the 21 RDPs implemented in Italy.14 Clearly, expanding the treatment options well beyond the four submeasures would greatly affect the sample size of each subgroup and challenge the estimation of any HTE under the standard conditions discussed in Section 5 (Heiler and Knaus 2022). Finally, focusing on more specific measures does not necessarily imply a more refined outcome variable (see Section 4 for further discussion).15

Since organic farming (measure 11) is homogeneous across the RDPs and involves a reasonable number of farms (271), we repeat our analysis by redefining treatment T₂ as a two-versions treatment T₂ = (T_2o, T_2n), where o = organic and n = nonorganic. Given our initial definition of the treatments, T_2n coincides with measure 10 which, unlike measure 11, is not entirely homogeneous across RDPs and could be exposed to further treatment heterogeneity. We therefore estimate the HTE of T₂ under two different setups: (1) we analyze the HTE of participating to AEMs as in Stetter, Mennig, and Sauer (2022); (2) we break down the treatment in setup 1 into T_2o and T_2n and obtain the corresponding HTE; and (3) we compare the results from setups 1 and 2 and discuss their implications for the interpretation of the HTE of interest (see Section 6).

Outcome Variable

The theoretical framework presented in Section 3 expresses the farm response to the treatment as Δy_it,k, that is, a vector whose nonzero elements represent all of the farmer’s production choices associated with the treatment in terms of both input and output.16 These elements may consist of a long list of the farmer’s specific production decisions, ranging from crop and livestock management practices to water and nutrient use (Burton and Schwarz 2013; Guerrero 2021, 11). One way to reduce the dimensionality of Δy_it,k consists of identifying and extracting the elementary indicators expressing the change in farming practices toward extensification or environmentally friendly practices. However, as discussed in Sections 2 and 3, focusing on elementary indicators might cause ambiguity when interpreting treatment effect heterogeneity because of the OTH problem. Given the potential correlation among the components of Δy_it,k, one way to retain all the information in the netput vector while avoiding multiple marginal evaluations is to perform dimension reduction (Chipman and Gu 2005) to obtain composite dimensional indices (Bartolini et al. 2021). This strategy not only provides an insulation against OTH but also resonates the need for a comprehensive evaluation of complex policy instruments such as the AEM discussed in Sections 2 and 4. As also argued by Stetter, Mennig, and Sauer (2022, 727), despite the articulation of AEMs in specific submeasures, the goal of the AEPs remains more general, aiming to improve the overall environmental performance of the agricultural sector. Although many studies have tried to evaluate the effectiveness of distinct AEPs with respect to specific policy targets (e.g., the impact on biodiversity), the integrated assessment of multifaceted goals involving, for example, soil and water protection and the curbing of greenhouse gas (GHG) emissions have received relatively little attention until recently (Hudec et al. 2007; Zhen et al. 2022). However, the literature has long suggested that the intricate and ecosystemic nature of the agri-environment requires that any assessment should be based on a comprehensive integration of indicators across many environmental dimensions (Wascher 2003; Purvis et al. 2009).

In this respect, Purvis et al. (2009) propose an interesting, harmonized approach to evaluating AEMs: the so-called agri-environmental footprint index (AFI). The AFI expresses a multidimensional assessment as a univariate index that can be flexibly adapted to diverse contexts. We use the AFI framework as adapted by Westbury et al. (2011) with the FADN data. We refer to this methodology as FADN-AFI, as the resulting index uses elementary information included in the FADN dataset. We extend the FADN-AFI to evaluate whether and to what extent the implementation of the CC requirements, GPs, and AEMs meet the CAP 2015–2020 environmental objectives.17

Table 2 presents the elementary components of our FADN-AFI (see Appendix Table B2). The land use diversity indicator (the Shannon index) is detailed in Appendix A. Appendix B discusses the definition of a farm-level GHG emissions indicator using farm-level information. This measure should provide a reliable proxy of the contribution of a farm’s practices to climate change mitigation (Dabkiene, Balezentis, and Streimikiene 2021). The FADN-AFI’s elementary components are then standardized to obtain dimensionless z-scores that we eventually aggregate using the weights indicated in the last column of Table 2 (i.e., giving a positive or negative sign for positive or negative environmental externalities, respectively).18 The resulting FADN-AFI is monotonic in farms’ environmental performance in that higher FADN-AFI scores correspond to “better” environmental performance. Since the range of the FADN-AFI is not bounded, the index might be difficult to interpret per se. However, since HTEs are defined through pairwise differences, these can easily be understood comparatively. Finally, we average the FADN-AFI in 2017–2018 to provide more stable values for the outcome variable.19

Confounding Variables

As discussed in Section 3, the choice of covariates entering the X_it vector becomes crucial for identifying the HTEs of interest. These should encompass farm heterogeneity as extensively as possible, thereby allowing fair comparisons between treated and untreated units. Selecting all the relevant confounders such that the assumptions outlined in Section 5 are satisfied may follow multiple routes. On the one hand, one may construct a very large collection of internal farm characteristics and external socioeconomic indicators that might explain the individual decision of adopting one of the treatments. In this case, we would let ML algorithm choose which feature contributes the most to predict farmer behavior through a regularization mechanism. However, as recently outlined by Hünermund, Louw, and Caspi (2023), this strategy may lead to severely biased treatment effects if the covariate set includes potentially endogenous confounding variables. Ultimately, the authors advocate that when the goal is conducting CI, researchers need to justify the controls they want to include and, more importantly, make sure that these are exogenous (i.e., pretreatment).

For these reasons, we begin by defining the confounders in S_i and V_it through an extensive literature review covering several empirical studies addressing farmer participation in AEPs and the impact of AEPs on farms’ economic and environmental performance. The results of this survey are displayed in Table 3, where the list of covariates resulting from this desk research is classified using the taxonomy elaborated by Brown et al. (2021) and discussed in Section 3. We invite the reader to refer to the individual studies for a throughout explanation of how these regressors are relevant for the research questions. The abundance of controls compiled in this long list might suggest some form of preliminary selection to avoid redundancy and achieve a more parsimonious set of variables. Nevertheless, unlike most parametric econometric tools, forest-based ML algorithms can easily accommodate multiple overlapping information sources and use them to either create intermediate features or discard redundant ones through regularization. Therefore, our empirical analysis makes use of all the covariates in Table 3.20

To satisfy the identifying conditions anticipated in Section 3, the time-varying controls, V_it, must be exogenous with respect to the treatment (i.e., predetermined). In theory, this would preclude the use of certain direct measures of farm physical and economic size, such as utilized arable land, profit, revenue, costs, and total workforce. To circumvent this issue, some authors suggest using covariates measured before the introduction of the treatment (for studies assessing AEMs, see, e.g., Bertoni et al. 2020; Uehleke, Petrick, and Hüttel 2022; Stetter, Mennig, and Sauer 2022). However, this strategy is sometimes infeasible, as such measurements may not be available if the policies under investigation were introduced several years before the outcome is measured. When this happens, going back in time may imply a major loss of observations. This concern is particularly relevant for our application, as the rotating structure of the Italian FADN panel shows that 582 farms (approximately 15% of the sample) included in the 2015–2018 dataset are not present in the 2014 data. Therefore, our choice is to follow the strategy of Arata and Sckokai (2016) and Pufhal and Weiss (2009), which consists of using the first year since the introduction of the policy as the pretreatment period (2015, in this case).21 Notice that since our outcome variable is calculated using the years 2017 and 2018, V_it contains lagged (by two years) elementary components of the FADN-AFI. Moreover, since farms usually sign up for participating in certain AEMs over several years (Bertoni et al. 2020; Uehleke, Petrick, and Hüttel 2022), we also include information on previous participation to such programs in V_it (Chabé-Ferret and Subervie 2013). Appendix Tables C1 and C2 report descriptive statistics for the outcome variable and all the control variables discussed above.

Unobservable Characteristics

The theoretical derivation in Section 3 provides the behavioral foundation of the farmer’s treatment choice and response to the treatment. This behavior depends on some observable characteristics but also on unobservable farm characteristics, u_i. The conditional independence between any of the treatments and the corresponding potential outcomes also hinges on the last component of the conditioning vector Z_it, namely, the unobservable farm characteristics, u_i. If these latent features influence the choice between T_i,1 and T_i,2 and the corresponding potential outcomes, the identification of the HTE becomes challenging because of the violation of unconfoundedness. Even though X_it can be extended to collect as many observable farm characteristics as possible, this strategy may be insufficient to insulate against selection-on-unobservable. Policy conclusions drawn from the HTE estimation could be problematic and even erroneous if the relevance of these unobservables and their possible association with the observable characteristics are not properly investigated and understood.

In these situations, ML methods (including BCFs) can help in identifying automatically creating nonlinearities and complex interactions among the variables in X_it, generating artificial strata that allow more precise comparisons between treated/untreated units and their counterfactuals. These “synthetic traits” not only greatly expand the initial set of confounders but also correlate with the unobservable characteristics, thereby making the unconfoundedness assumption more credible. This argument is also put forward by Stetter, Mennig, and Sauer (2022, 738–39, 744), who provide a nice example of how this property of ML techniques may help to control for farmer attitudes toward environmental issues.22 Since this is not directly testable, we check the robustness of the above propositions through several sensitivity analysis tests. As illustrated in Appendix H, we probe the stability of our results in the presence of omitted variable bias from unobserved endogenous heterogeneity by introducing synthetically generated u_i into the covariate set. See Section 6 for more details and caveats of this approach.

Estimating Treatment Effects via BCF

The estimation of HTEs using the BCF algorithm requires the usual assumptions of unconfoundedness and SUTVA, which can be expressed as follows: 1 where Y_i represents the FADN-AFI defined in Section 4, X_i indicates the vector of confounders defined in Section 4, while Y_i(1) and Y_i(0) indicate potential outcomes for individuals in a treatment group (T_i,k = 1) or control group (T_i,k = 0), respectively (Imbens and Rubin 2015, ch. 1). Notice that SUTVA implies no hidden variations of the treatment. As discussed in Sections 2 and 3, binarized multiple-versions treatments can lead to violations of this assumption unless one imposes stringent restrictions on the treatment assignment mechanism. For example, in case any individual i with characteristics X_i can only choose one of the hidden treatments, SUTVA is still a credible assumption (VanderWeele and Hernán 2013; Lopez and Gutman 2017). As previously discussed, we make this assumption for the treatments defined in Section 4, except for the distinction between organic and nonorganic farming. We therefore set k to k ∈ {1,2} such that T_i,k = 1 indicates either T_i,1 = 1 or T_i,2 = 1, while T_i,k = 0 always refers to farms in the control group. We discuss the implication for disaggregating T_i,2 into T_i,2o and T_i,2n in Section 6. For notational convenience, we drop the subscript k. Of these elements, we only observe the potential outcome that corresponds to the realized T_i, namely, Y_i = T_iY_i(1) + (1 – T_i)Y_i(0). Equation [1] postulates independence between the potential outcomes and the treatment, conditional on the set of exogenous variables, X_i.

Combining unconfoundedness, SUTVA, and overlap (as discussed in Section 4) allows the estimation of causal effects via strong ignorability; that is, E[Y_i(t) | X_i = x_i] = E[Y_i | T_i = t_i, X_i = x_i], with t_i ∈ {0,1}. The latter implies that the estimand of interest is simply the difference between two conditional expectation functions: 2 where τ(x_i) is typically referred to as a conditional average treatment effect (CATE). Since one can use μ_T(x_i) to impute conditional treatment effects at the individual level, equation [2] is sometimes referred to as individualized average treatment effect (IATE) (Lechner 2018; Knaus, Lechner, and Strittmatter 2021, 2022). This estimand represents the most disaggregated form of HTE.

Often researchers may be interested in subgroups or intermediate aggregation levels of the exogenous covariates, leading to the definition of group average treatment effects (GATEs): 3 where ϕ(.) represents a generic probability density of mass function, G_i denotes the collection of possible groups, and g_i denotes one such group. GATEs have recently gained considerable attention in the applied literature as treatment effect heterogeneity is often better understood for subsets of the population (Lechner 2018; Lee, Bargagli-Stoffi, and Dominici, 2020). ATEs can also be obtained by averaging the IATEs over the full distribution of X_i: 4 To estimate the IATEs (and then the GATEs and ATEs), we assume that the data-generating process for follows a stochastic process defined as follows: 5 where f indicates an arbitrarily complex function25 and ε_i represents an additive idiosyncratic error term ε_i ~ N (0,σ²), independently distributed.

In this context, E[Y_i | T_i = t_i, X_i = x_i] = f (x_i, t_i) therefore, at least in principle, τ(x_i) can be estimated by the simple difference f (x_i, t_i = 1) – f (x_i, t_i = 0) = μ₁(x_i) – μ₀(x_i), as illustrated above. However, as discussed by Künzel et al. (2019) and Nie and Wager (2021), training two separate conditional mean functions and taking their difference may produce highly unstable estimates. For this reason, Hahn, Murray, and Carvalho (2020) proposed a slightly different approach, wherein the expected value of the outcome of interest has two components: a prognostic function, 𝓂(x_i) plus an additive heterogeneous treatment effect, τ(x_i): 6 where both 𝓂(.) and τ(.) represent stochastic functions with BART priors, namely, 𝓂 ~ BART(θ | PS(x_i), x_i) and τ ~ BART(ϑ | x_i), and PS(x_i) indicates the estimated PS. The two vectors θ and ϑ collect the hyperparameters regulating the number of trees in the BART ensembles, their depth, and the splitting rule associated with each single tree (see Appendix F for details). As previously mentioned, the specification in equation [6] allows regularizing τ(x_i) directly and independently, thereby reducing the noisiness of the IATEs with respect to the same estimates obtained from simple differences in conditional mean functions. Furthermore, the additive nature of equation [6] ensures that the prior on f (x_i, t_i) is also a BART (Chipman, George, and McCulloch 2010; Hill, Linero, and Murphy 2020). Finally, notice that the model presented in equation [6] also appears in Nie and Wager (2021), who propose a frequentist approach to estimating τ(x_i). In contrast to the setup discussed above, however, the authors propose a residuals-on-residuals reparameterization of equation [6] which is then used to obtain (regularized) consistent estimates of τ(x_i) via a two-stage optimization procedure.

The full Bayesian model requires the definition of a likelihood function for the outcome variable (Gelman et al. 2013; McElreath 2020). Consistent with equation [5] and Chipman, George, and McCullogh (2010), Hill (2011), and Hahn, Murray, and Carvalho (2020), we employ a normal model for Y_i, along with a semiconjugate inverse chi square prior for its variance: 7 where ω is set following Chipman, George, and McCullogh (2010) (see Appendix F for further details). Samples from the posterior distribution of τ(x_i) are obtained via Markov chain Monte Carlo sampling, as implemented in the R package bcf. We indicate posterior draws from ϕ(τ(x_i) | x_i, t_i, y_i,…, y_N) as {τ ^s(x_i)}^S_s=1, where S indicates the number of Markov chain Monte Carlo simulations.

Subgroup Search via Shallow Regression Trees

The approximated posterior {τ ^s(x_i)}^S_s=1 is a multivariate probability distribution over a complex P-dimensional function, and as such, it might be difficult to interpret directly. One way to compress such information consists of obtaining marginal distributions of the IATEs for one covariate of interest and plotting them against the full range of that variable. A similar approach was adopted by Stetter, Mennig, and Sauer (2022), who used Shapley values (Shapley 1953) to identify the marginal contributions of several treatment effect drivers and used these indicators to construct partial dependence plots. Another sensible approach to investigating IATE heterogeneity consists of comparing farm subgroups obtained by projecting the full posterior distribution onto a lower-dimensional covariate space. In this respect, we follow the work of Yeager et al. (2019), Hahn, Murray, and Carvalho (2020), Woody, Carvalho, and Murray (2021) and (and partially Lee, Bargagli-Stoffi, and Dominici 2020), who suggest eliciting the relvant subgroups by partitioning the IATE maximum a posteriori (MAP) estimates, τ_i = S ^–1∑^S_s=1 τ ^s_q(x_i), using shallow regression trees (CART) (Breiman 1984). Specifically, the authors propose to split τ_i along w_i, where w_i ⊆ x_i indicates a vector of policy-relevant variables and setting w_i ⊂ x_i implies using domain knowledge to enforce an initial regularization of the resulting tree. We restrict our attention to a subset of simple and understandable characteristics that policy makers might find helpful to improve the targeting of AEMs (see Section 6). Once farm subgroups have been identified, GATEs can be obtained as weighted averages of the IATEs that fall into each cluster. This approach to calculating GATEs is also consistent with Lechner (2018) in that group-level effects are obtained as convex combinations of the IATEs. In our application, however, weighting is automatically performed when fitting a tree to τ_i.

Finally, for some potential effect moderator x_p ∈ x, the comparison between pointwise estimates (or intervals) computed at different levels of x_p ignores any potential correlation between IATEs along other variables x_l, for all p, l ∈ {1,…, P}. In other words, the marginal distribution of τ(x_p,i) disregards the information encoded in the correlation between τ(x_l,i) and τ(x_l,–i) when x_l,i and x_l,–i are close. This might lead to misleading comparisons along x_p and, consequently, unreliable policy implications. Therefore, once the relevant subgroups have been identified, one can obtain the full posterior distribution of each pairwise difference as: ϕ_g1,g2 = ϕ(τ_i|i∈g1 – τ_i|i∈g2), where g₁ and g₂ indicate any two subsets of τ_i.

IATEs

The two graphs numbered “1” in Figure 1 display the MAP; that is, the average over the S samples from the posterior distribution of τ(x_i) estimates and corresponding 95% confidence intervals (CrI) of the IATEs over the two treatment comparisons. These are ordered across the respective samples from the lowest to the highest individual value. We start our discussion by presenting the results for T₂, the treatment that is more frequently addressed by the literature. First, it is worth noting that overall, the modal direction of the responses to the treatment (T₂) is fully consistent with theoretical expectations: adding the AEM to the environmental standards implied by the CC and the GP (Figure 1a [1]) induces an improvement in the FADN-AFI, that is, in the farm-level environmental performance. The opposite response is observed when the environmental standards implied by the CC and GP are dropped (i.e., treatment T₁) (Figure 1b [1]). Whereas in the first case, most estimated IATEs exhibit CrI not including zero (black dots), the converse applies to the second comparison group, for which a large proportion of farms have inconclusive individual-level TEs (light gray dots). The first graph in Figure 1b also indicates that some farms might even exhibit opposite responses, although the corresponding IATEs appear quite noisy. This evidence is presented in greater detail in Table 4, which provides descriptive summaries of our main results.

The two graphs numbered “2” in Figure 1 show the IATE’s MAP frequency distribution for the two cases. These plots highlight the variability of the responses, with few cases showing a treatment effect direction that conflicts with the expected direction (despite exhibiting CrI including positive and negative values). Apart from these rare extreme cases, however, our MAP estimates range between roughly 0.1 and 1.0 for treatment T₂ and between approximately −3 and 1.5 for treatment T₁. The nature and determinants of these different patterns can be further investigated by estimating GATEs, as addressed in the next section.

The irregularity of farms’ responses to the treatments is a clear sign of heterogeneity, one that would be lost by the mere inspection of ATEs (see the two graphs numbered “3” in Figure 1). Whereas these latter aggregated estimands provide clear indications of policy effectiveness (as both show an effect in the expected direction), the inspection of the IATEs tells a different and more subtle story. This is especially true for the treatment T₁, whereas the responses seem more homogeneous when studying the treatment effect of implementing CC and GP requirements together with AEMs (treatment T₂).

Finally, for each individualized treatment effect, we calculate the posterior probability that the corresponding IATE is either greater than zero or lower than zero for the T₂ and T₁, respectively. Our results show that when comparing farms implementing CC and GP requirements plus AEMs with the control group, most of the IATEs’ posterior distributions lie above zero. For example, the proportions of IATEs with at least 60%, 75%, and 90% positive posterior are 100%, 88.5%, and 5%, respectively. Conversely, when comparing the control group to farms with no adherence to CC or GP, the posterior distributions of their IATEs are largely negative. In this case, the proportions of IATEs with at least 60%, 75%, and 90% negative posterior are 83%, 15%, and 0%, respectively.

Notice that all the results discussed thus far are based on observations satisfying the common support as defined by rule I in Appendix E. Under such a restriction to the range of X_i, however, our dataset does not suffer drops. The sensitivity of these figures to different exclusion rules is discussed in Section 6 (robustness check), in which the selection method we used based on the estimated PS is also discussed.

Figure 1

Estimated Individualized Average Treatment Effects: (a) Treatment Group T₂ (Farmers Implementing Agri-environmental Measures); (b) Treatment Group T₁ (Farmers Not Fulfilling Conditionality Restraints or Implementing Agri-environmental Measures)

GATEs

We partition the posterior distribution of τ(x_i) using a set of policy-relevant measures w_i covering the most relevant dimensions of heterogeneity, as evidenced by the measures of feature importance produced by the BCF. Our characterization of w_i involves (1) examining the variable importance metrics generated as a by-product of the fitting model [7],26 and (2) choosing the 10 most predictive dimensions that policy makers might target to improve the effectiveness of AEMs. We fit a CART algorithm to τ_i using the attributes selected using the procedure illustrated above: latitude, longitude, altitude (geographical location); total arable land, share of rented land, revenue (physical or economic size); farm specialization (relative importance of the first and second crop, farms specialized in livestock, crop and livestock farms, farms specialized in annual crops, and farm specialized in perennial crops). The results for the two treatments are shown in Figures 4 and 5, wherein, for the sake of interpretability, we do not allow the trees to split more than three times.

When we consider the adoption an AEM in addition to CC and GP requirements (treatment T₂) (Figure 2b [1]), we find that TE heterogeneity is mostly associated with five variables: latitude, physical farm size, altitude, crop specialization (share of the second crop in the crop mix), and livestock intensity. These covariates trace out eight subgroups with different levels of treatment effects. For example, subgroup g₈ exhibits the lowest treatment effect and consists of farms in southern Italy with less than 85 ha arable land. On the opposite end of the spectrum, we find subgroup g₁₅, which comprises crop-specialized farms in northern Italy with low livestock intensity. One can then obtain the full posterior distribution of g₁₅ – g₈ with 95% CrI between −0.27 and 0.49 (Figure 2b [2]), which indicates that the difference between the two subgroups is in fact small, if not zero. Interestingly, if we repeat this exercise across all the leaves defined by the tree in Figure 2b [1], no group differences emerge (see Appendix G). These results are consistent with our discussion in Section 6 (i.e., our preliminary findings suggested limited treatment effect heterogeneity for treatment T₂).

In the case of treatment T₁ (Figure 2a [1]), we see that the shallow tree picks up four moderating variables: specialization in perennial crops, latitude, altitude, and livestock intensity. In this case, the subgroup with the strongest TE is g₈, which consists of farms specialized in perennial crops in Italy’s southmost regions. Subgroup g₁₅ includes observations from farms in the Po Valley that are not specialized in perennial crops. The difference in TE between these subgroups lies approximately between −2.2 and −0.41 (95% CrI; Figure 2a [2]), indicating the presence of treatment effect heterogeneity. Repeating this exercise across all the terminal nodes, we find that unlike treatment T₂, when the treatment consists of dropping both CC and GP requirements, many groups exhibit diversified responses. These further details are provided in Appendix G, where we also provide a deeper tree to gain further insights into these HTEs and a graphic representation of the geographical distribution of the IATEs.

It is finally worth stressing that although our main goal is to explore which observable farm characteristics exhibit a greater heterogeneity of response, some of these features might not be easily addressed by AEPs due to cost constraints or infeasibility or because they could potentially lead to discriminatory outcomes. From a policy perspective, it would be more useful to evaluate the level of heterogeneity associated with covariates that can be targeted more easily and effectively through policy measures. Most of the geographical features considered in our study, along with variables indicating long-term farm production specialization, appear particularly suitable for this purpose. In this respect, our results confirm that most of these geographical features significantly contribute to the observe heterogeneity of response. Similarly, the presence of perennial crops, crop specialization, and livestock density, all of which relate to distinct and consistent farming practices, pinpoint to patterns of strong heterogeneity. This suggests that AEPs could significantly enhance their effectiveness by specifically targeting these features. For a more detailed discussion on this matter, please refer to Appendix G.

Figure 2

Shallow Regression Tree Fitted to the Maximum a Posteriori (MAP) Individualized Average Treatment Effects: (a) Treatment Group T₂; (b) Treatment Group T₁

Robustness Checks

We check the consistency of our results to the assumptions formulated in Sections 4 and 5. Our first robustness check concerns the common support condition. As anticipated in Section 5 and further detailed in Appendix E, we use both the posterior distribution of the BART algorithm and a PS-based algorithm to investigate common support. Our tests show that the results presented in Section 6 are robust to these different methods to achieve overlap (see Appendix Tables H1 and H2).

We perform a battery of tests that largely encompass those discussed by Stetter, Mennig, and Sauer (2022) in that we reestimate our BCF multiple times, each time manipulating different model features. We begin by probing unconfoundedness through a recursive procedure in which we fit model [7] after dropping: (1) the most important feature in terms of relative frequency in the forest, (2) the three most important features, and (3) the five most important features. As detailed in Appendix Figures H1–H3), this exercise yields the first indication that the BCF in equation [7] is fundamentally resilient against unobserved heterogeneity as long as this is associated with the set of observed confounders. Put differently, the complex interactions and nonlinearities generated by the tree ensemble seem to work as additional synthetic controls associated with the left-out covariates, thus compensating for their absence in the model. However, this line of reasoning hinges on the (strong) assumption that the most predictive features are also associated with both Y and T_k. In case this assumption fails, the procedure discussed above cannot be interpreted as a robustness check for unconfoundedness. For this reason, we build on these preliminary results and devise an additional test targeting endogenous unobserved heterogeneity directly. Our strategy consists of generating a random variable correlated with both Y and T_k, forming the vector Z_it as described in Section 3, and rerunning the model. As shown in Appendix Figure H4, our results do not change substantially, even under a strong imposed association between the unobserved variable and (Y, T_k). This stability could result from the properties of the BART ensemble in that when the forest is dense, the marginal contribution of each covariate becomes increasingly small (Chipman, George, and McCulloch 2010). Alternatively, it could be that the correlation between the nonlinear interactions generated by the BCF and the new confounder is strong enough to prevent distortions in the IATEs. In either case, it is worth warning that treatment effect estimates might deteriorate quickly when unobserved heterogeneity is more abundant and complex. This test is in fact only restricted to a single unobserved factor, which we model as linearly associated with Y and T_k (i.e., through correlations, which do not necessarily imply a direct effect of the synthetic u_i on the outcome or the treatment). We thus expect that in presence of multiple endogenous latent confounders, possibly related to the treatment and the outcome (or other elements of X_i in a nonlinear fashion, our estimates might tun out sensibly different. Although the literature offers other methods to perform sensitivity analysis with respect to omitted confounders (Dorie et al. 2016; VanderWeele and Ding 2018), we believe that they either do not overcome the limitations discussed above or they remain difficult to implement in HTE estimation. Therefore, despite the promising results presented so far, we stress that these only hold if several important restrictions are met.

The following robustness check consists of creating both a placebo treatment and a placebo outcome, replacing their observed counterparts in equation [7], and fitting the model two more times. If the model is correctly specified, the IATEs resulting from these “fake” variables should be uncorrelated with τ(x_i). As Appendix Figures H5 and H6 show, the new results obtained through placebo treatments and outcomes not only have no correlation with our estimated IATEs but also produce zero ATE with minimal treatment effect heterogeneity.

Finally, we assess the robustness of the estimated IATEs with respect to the OTH problem discussed in previous sections. We proceed by replacing the FADN-AFI with its elementary components and reestimating model [7] as discussed. Appendix Figure H7 suggests that focusing on marginal indicators produces TEs whose individual directions are essentially in line with those presented in Section 6. For example, implementing AEMs seems to yield lower GHGs, higher crop diversity, lower fertilizer expenditure, and more woodland areas. Nonetheless, a noteworthy difference emerges in terms of treatment effect heterogeneity. Whereas adopting the FADN-AFI points to a limited diversity across farms, using marginal measurements would suggest that treatment T₂ is environmentally beneficial only when the treatment effect is large. For this reason, our results invite to caution when it comes to choosing the dependent variable of model [7]. Although addressing individual indicators may appear more attractive and interpretable, it is worth stressing that missing out on the potential correlation or interdependence among them can affect the TE estimates in a nontrivial way.

Figure 3

Individualized Average Treatment Effects Differentials for the Two Versions of T₂

Role of Heterogeneous Treatments

As discussed in Sections 3 and 4, one potential limitation of our results (as well as other works investigating HTE of aggregated treatments) is that part of the estimated treatment effect heterogeneity in T₂ might be a statistical artifact. This would result from the fact that T₂ is a multiple-versions treatment as it aggregates two distinct measures which admit, in turn, several submeasures (see Section 4). As introduced in Section 5, the presence of treatment heterogeneity may affect our results by violating SUTVA (Heiler and Knaus 2022). Since in this case, the resulting interpretation of τ(x) would be misleading, we reestimate model [7] replacing T₂ with the two respective measures (measure 11 and measure 10) and approach the problem from a multiple-versions treatment perspective as discussed in Lopez and Gutman (2017).

To assess the possible bias in HTE estimation due to treatment heterogeneity, we compare the posterior distribution of the IATEs presented in Section 6 with the posterior density of the IATEs estimated using either T_2o, τ_2o(x_i), or T_2n, τ_2n(x_i). Figure 3 shows the 95% CrI for the differences τ_2o(x_i) – τ₂(x_i) and τ_2n(x_i) – τ₂(x_i), respectively, where τ₂(x_i) indicates the IATE for individual i under treatment T₂. As we can see from these plots, the difference between our initial estimates and those obtained by substituting T₂ with T_2o are minimal. Indeed, although τ_2o(x_i) is on average (black line in the left graph in Figure 3) slightly smaller than τ₂(x_i) for all i ∈ N_o, where N_o indicates the number of units choosing T_2o, all the CrI include both positive and negative values. At the same time, when focusing on T_2n, we see that τ_2n(x_i) – τ₂(x_i) are on average higher than zero for all i ∈ N_n, where N_n indicates the units choosing T_2n. However, the CrI once again includes zero for all such comparisons, although they are all moderately skewed toward positive values. Moreover, as mentioned in Section 4, T_2n could still entail some degree of treatment heterogeneity, which recommends caution when interpreting the corresponding estimates. Overall, examining the two measures separately highlights that the posterior distribution of the IATEs does not seem to change markedly when the aggregated (T₂) or the disaggregated (T_2o, T_2n) treatment is considered. This would suggest a limited impact of treatment heterogeneity on our interpretation of the HTEs discussed above. Nonetheless, further research effort remains desirable to better clarify the possible role of multiple versions in the correct identification and estimation of the HTE.