Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,671 characters · 18 sections · 82 citation commands
Heterogeneous Responses to Continuous Treatments: A Cluster-Based Causal Framework
\def\spacingset#1{ {#1}} \spacingset{1}
\if1\anon \fi
\if0\anon {
} \fi
\abstract{When treatments are non-randomly assigned, continuous, and yield heterogeneous effects at the same intensity, causal identification becomes particularly challenging. In such contexts, existing approaches often fail to provide policy-relevant estimates of the relationship between treatment intensity and outcomes, especially in the presence of limited common support. To fill this gap, we introduce the Clustered Dose-Response Function (Cl-DRF), a novel estimator designed to uncover the continuous causal relationship between treatment intensity and the dependent variable across distinct subgroups. Our approach leverages both theoretical and data-driven sources of heterogeneity, relying on relaxed versions of the conditional-independence and positivity assumptions that are plausible across various observational settings. We apply the Cl-DRF estimator to estimate subgroup-specific dose–response relationships between European Cohesion Funds and economic growth. In contrast to much of the literature, higher funding increases growth in more developed regions without diminishing returns, while limited absorptive capacity prevents other regions from fully benefiting.}
Continuous treatment, Heterogeneous treatment effects, Clustering, Potential outcome, European Union funds
Many treatments are non-binary and continuously distributed, making the identification of causal effects a nontrivial problem athey2017state,fong2018covariate. Moreover, the complexity of applied empirical contexts typically implies that treatments are allocated based on factors correlated with potential outcomes, violating random assignment and further complicating identification of the causal relationship between treatment intensity and outcomes.
The literature has introduced several estimators specifically designed to address the challenges posed by continuous treatments, significantly broadening the methodological toolkit available to researchers hirano2004propensity, imai2004causal, kennedy2017non, fong2018covariate, zhao2020propensity, wu2022matching, huling2023independence. These estimators are designed to recover an unbiased estimate of the causal Average Dose–Response Function (ADRF), a population-based estimand describing how the outcome of the average unit in the population varies with the treatment dosage tubbicke2023ebct. Identification relies on two fairly stringent assumptions: weak unconfoundedness and positivity\footnote{Weak unconfoundedness is also referred to as ignorability, no unmeasured confounding, or the conditional independence assumption; the positivity assumption is also referred to as overlap or the common support assumption.}. The first assumption states that, once we account for the observed pre-treatment characteristics, the potential outcomes are independent of the treatment assignment at each value of the treatment intensity, while positivity requires that, for every combination of these characteristics, each unit could, in principle, be observed at any treatment level. Under these assumptions, existing estimators for continuous treatments identify the ADRF, while accommodating treatment effect heterogeneity that is fully captured by observed covariates. However, these estimators are not suitable for empirical contexts characterized by treatment effect heterogeneity across units and limited common support across units. Indeed, we show that under these circumstances they fail to consistently recover the true ADRF and, moreover, that the ADRF itself is not the most policy-relevant estimand.
Treatment effect heterogeneity may arise when the causal effect differs systematically across subpopulations, even when units receive the same treatment intensity. These variations may reflect differences in baseline characteristics, contextual factors, or differential responses to the treatment. Although the presence of heterogeneous effects is well established in the binary treatment literature athey2016recursive, wager2018estimation, de2020two, callaway2021difference, it remains largely overlooked in the continuous treatment literature. Revealing systematic differences in treatment responses necessitates subgroup-specific estimates that are informative about heterogeneous effects and relevant for policy design.
To bridge this methodological gap, we introduce the Clustered Dose Response Function (Cl-DRF) approach. The Cl-DRF characterizes the continuous causal relationship between treatment intensity and outcomes by allowing for multiple dose response functions (DRFs) indexed by clusters, rather than maintaining a single population average ADRF. We first define a novel causal estimand, namely cluster-specific DRFs, and then propose an estimator to recover them. The Cl-DRF framework allows the causal relationship between treatment intensity and outcomes to vary across distinct subgroups or clusters\footnote{Throughout the paper, we use the terms subgroups and clusters interchangeably.}, each governed by a potentially different underlying mechanism. Following hirano2004propensity, we adopt a parametric strategy to model treatment assignment and outcomes. Although less flexible than fully nonparametric alternatives kennedy2017non, colangelo2025double, this choice provides the tractability needed to jointly estimate the clustering structure and the cluster-specific DRFs. A key feature of the Cl-DRF estimator is its reliance on assumptions that are plausible across various observational settings. Indeed, while the Cl-DRF estimator adheres to the weak unconfoundedness and positivity framework, it relies on a less stringent version of these assumptions, requiring their validity exclusively within the identified subgroups. For instance, while it is highly unlikely for every unit to maintain a non-zero conditional density across all treatment levels branson2023causal, meeting the positivity assumption becomes more plausible within a subgroup characterized by similar covariate values, contextual factors, responses of the outcome to the treatment and, possibly, a smaller range of treatment intensities. By mitigating concerns about assumptions' validity, Cl-DRF can be applied to a broader range of empirical contexts, particularly those characterized by likely heterogeneous treatment effects.
Since in most empirical applications treatment effect heterogeneity cannot be known a priori, the Cl-DRF approach accommodates both theory-driven subgroup definitions, grounded in observed covariates, and data-driven procedures for subgroup identification. When implementing the data-driven version of the approach via the Cl-DRF estimator, we do not assume the number of clusters to be known in advance. Instead, the number of clusters is treated as an unknown parameter and jointly estimated within the Cl-DRF approach. Moreover, we do not assume a priori that the likely heterogeneity depends solely on one or more covariates. Instead, we consider that it could stem from the varying response of the outcome to the treatment, which is a function of the covariates. Our focus is on settings where not only the magnitude but also the shape of the dose–response relationship varies across subgroups, in ways that cannot be captured by simple interactions between covariates and treatment intensity. To determine the clusters and estimate the corresponding DRFs, we develop an iterative procedure that alternates between updating cluster membership and estimating cluster-specific treatment assignment models, summarized by the generalized propensity score (GPS), and cluster-specific outcome regressions. This structure mirrors standard clusterwise regression approaches spath1979algorithmus,demidenko2018next,sugasawa2021grouped, allowing heterogeneity to emerge from differences in the treatment--outcome relationship across clusters.
To illustrate the properties of the Cl-DRF estimator, we conduct a series of simulations designed to assess its performance under various scenarios. We next apply the Cl-DRF estimator to quantify the impact of the European Union (EU) Cohesion Policy on economic growth. This setting is particularly suitable for our framework, as EU funds are continuously distributed, non-random, and likely to generate heterogeneous returns depending on regions’ pre-treatment characteristics and absorptive capacity. This empirical application not only showcases the estimator’s ability to handle complex, real-world data but, more importantly, deepens our understanding of the economic effects of EU place-based interventions. The resulting causal estimates enable policymakers to fine-tune the allocation of funds by identifying clusters of regions where marginal returns are highest and those that are unable to benefit from additional support. These insights, in turn, help direct resources toward the most promising areas while mitigating the risk of inefficiently subsidizing regions with limited growth potential.
The paper structure is as follows. Section (ref) illustrates a motivating example with simulated data, while Section (ref) discusses previous studies dealing with causal inference for continuous treatment and treatment effects heterogeneity. Section (ref) develops the proposed methodology by formalizing the causal framework and the target estimand, introducing the Cl-DRF estimator under its identifying assumptions, and addressing the testing of the cluster structure hypothesis and the selection of the number of clusters. Section (ref) applies the proposed estimator to derive policy-relevant estimates guiding the allocation of European Union (EU) funds assignment to stimulate regional economic growth, while Section (ref) concludes with final remarks.
Observational studies in the social sciences often involve continuous treatments, heterogeneous treatment effects, and non-random treatment assignment that induces systematic associations between treatment intensity and observed covariates. These associations can generate severe limitations in common support: for some covariate profiles, it is difficult, both conceptually and empirically, to observe units receiving widely different treatment intensities.\footnote{In other disciplines, this limitation is generally less significant, as it is often easier to identify similar units subjected to widely varying treatment intensities. For example, in medicine, patients hospitalized for a heart attack may share similar characteristics yet exhibit markedly different levels of creatinine, representing the quantitative exposure austin2019assessing. Similarly, in environmental health, comparable areas may be affected by pollution at significantly different dosages wu2022matching.} Consider, for instance, public subsidies to private firms. Subsidy amounts are continuously distributed and depend on regional and sectoral characteristics, as well as on managerial quality. As a result, the relationship between treatment intensity and covariates may differ across clusters. Moreover, the effect of subsidies on firm outcomes, such as productivity, employment, or innovation, is likely to vary across these clusters, reflecting differences in market conditions and in the quality of investment projects. In such settings, estimating a single ADRF over a wide range of treatment intensities may not be credible. By contrast, it may be feasible to construct theory-based or data-driven clusters, each characterized by a more limited range of treatment intensities. Under these circumstances, cluster-specific DRFs can be estimated. The motivating example below mimics this setting.
Suppose there are four groups of units that differ both in treatment assignment and in the way outcomes respond to treatment. Within each group, treatment intensities lie in a relatively narrow range, common support is limited across groups, and treatment effects may differ in sign and exhibit nonlinearities. If this underlying structure is ignored and a single ADRF is estimated, these heterogeneous response patterns are averaged out. Moreover, the resulting ADRF is likely to be biased due to violations of the positivity assumption. Let $i= 1, \dots, n$ index the observational units in the sample. Let $c=1,\dots,C$ index the clusters, and assume that each unit $i$ belongs to exactly one cluster. Let $\mathbf{y}$ $(y_i, i=1,\dots,n)$ denote the $n$-dimensional vector of observed outcomes. Moreover, let $\mathbf{t}$ $(t_i, i=1,\dots,n)$ be the $n$-dimensional vector of observed treatments, where $t_i \in \mathcal{T} = [t_{\min}, t_{\max}] \subset \mathbb{R}$, and let $\mathbf{X}$ be the $n\times p$ matrix of pre-treatment covariates, where $\mathbf{x}_i$ is the $p$-dimensional vector of covariates associated with each $i$-th unit. We consider a sample size of $n=800$ units that are clustered into $C=4$ equally sized groups and $p=2$ vectors of pre-treatment covariates. For notational convenience, we augment the covariate vector with an intercept, $\tilde{\mathbf{x}}_i = (1, x_{1i}, x_{2i})'$, so that the treatment model
has cluster-specific parameter vectors $\boldsymbol{\beta}_c\in\mathbb{R}^3$. We set $\boldsymbol{\beta}_1 = [3,1.5,1.2]$, $\boldsymbol{\beta}_2 = [3,1.8,1.5]$, $\boldsymbol{\beta}_3 = [3,2,1.8]$ and $\boldsymbol{\beta}_4 = [3,2.2,2]$. We simulate $x_{1i} \sim U [0, 0.4]$ and $x_{2i} \sim U [0, 0.5]$ for Cluster 1, $x_{1i} \sim U [0.2, 0.6]$ and $x_{2i} \sim U [0.3, 0.6]$ for Cluster 2, $x_{1i} \sim U [0.5, 0.8]$ and $x_{2i} \sim U [0.5, 0.9]$ for Cluster 3, and $x_{1i} \sim U [0.7, 1]$ and $x_{2i} \sim U [0.7, 1]$ for Cluster 4, where $U[a,b]$ denotes the continuous uniform distribution with support $[a,b]$. Then, we simulate the outcome variable as a function of the treatment while allowing for confounding through the covariates,\footnote{In Appendix (ref) we model the outcome as a function of the treatment only. Since the treatment depends on covariates, the outcome implicitly depends on covariates as well. This reflects a strategy in which covariate–outcome relationships arise indirectly through the treatment assignment mechanism.} and assume the following within–cluster relationships:
where $e_i \sim N(0,1)$ for all $i$, and $\left(\gamma_{1c}, \gamma_{2c}\right)$ are cluster–specific moderation parameters determining how the covariates $(x_{1i}, x_{2i})$ modify the marginal effect of the treatment within cluster $c$. We set $(\gamma_{1c})_{c=1}^4 = (0.5,\,0.8,\,1.1,\,1.4)$ and $(\gamma_{2c})_{c=1}^4 = (0.3,\,0.6,\,0.9,\,1.2)$.
We remark that in this setting we assume two sources of heterogeneity: one in the covariate–treatment relationship and another in the treatment–outcome relationship. The first arises because treatment assignment depends on covariates in a cluster-specific (i.e., non-uniform) way. The second reflects the fact that units in different clusters exhibit heterogeneous responses to the same level of treatment. Both sources of heterogeneity are likely to coexist in practical settings. For example, subsidy amounts may depend on regional and sectoral characteristics as well as managerial quality; consequently, the relationship between treatment intensity and covariates may vary across clusters. Moreover, the effect of subsidies on firm outcomes, such as productivity, employment, or innovation, is likely to differ across these clusters, reflecting variation in market conditions and in the quality of investment projects.
Figure (ref) illustrates the association between treatment and outcome for each cluster. Notice that although the support for treatment varies among clusters, we display the true relationship between outcome and treatment also in areas outside the common support to take into account that theoretically each unit could have received any possible treatment intensity level, i.e., we assume that the positivity condition holds across the entire range of possible treatment intensities. However, due to the cluster-specific assignment mechanism, some treatment values occur with extremely low probability for particular groups, leading to very limited practical overlap in the observed data. The dashed black line represents the averaged relationship between outcome and treatment for the average unit (without considering the clustering information), i.e., the population ADRF.
In this scenario, we argue that estimating a single ADRF for the relationship between treatment intensity and the outcome obscures heterogeneity in the causal link and thus lacks policy relevance. Indeed, this averaging process overlooks the existence of clusters where a strongly positive (or negative) relationship between treatment and outcome is evident. After generating the data, we first estimate the ADRF by using both the parametric GPS approach of hirano2004propensity and the weighted nonparametric approach of huling2023independence (see Section (ref) for further details) to the whole set of 800 units, i.e., without taking into account the presence of clusters. The estimated cluster-based DRFs are displayed in Figure (ref).
The ADRF estimated using the hirano2004propensity estimator exhibits a convex shape and deviates substantially from the true simulated population ADRF, depicted by the dashed line. Similarly, applying the estimator proposed by huling2023independence yields comparable shortcomings. In this case, the estimated ADRF suggests a mildly nonlinear relationship between the outcome and treatment intensity, whereas the true relationship is linear and flat. From a policy perspective, this pooled ADRF would be interpreted as indicating modest returns at low treatment intensities and negligible or even negative effects at higher intensities. However, this pattern is purely an artifact of averaging across groups with sharply heterogeneous responses, with some groups benefiting from higher doses and others not. Therefore, in the setting under analysis, the ADRF is not only an estimand of limited policy relevance but also one that cannot be consistently estimated using existing methods.
Let us assume knowledge of the partition of the units into four groups as in Figure (ref), that is, we know exactly what units share similar treatment-outcome relationships. Under this assumption, we can apply the hirano2004propensity approach within each of these known groups to estimate cluster-specific DRFs. As shown in Figure (ref), the estimated DRFs (dashed black lines) perfectly resemble the true relationships (colored lines) within each cluster. The same result holds when using the approach of huling2023independence (see Figure (ref) in Appendix (ref)). This exercise shows that, conditional on knowing the true cluster structure, standard estimators are perfectly capable of recovering the correct DRF within each group. The challenge, therefore, is not the lack of flexible estimators per se, but the absence of a framework that jointly discovers the relevant clusters and estimates their corresponding DRFs. However, policymakers often lack precise knowledge about which groups may exhibit distinct treatment–outcome relationships, meaning that the assumption of known cluster structure is not realistic in practice. This raises the need for a novel modeling framework that contemporaneously identifies the grouping structure of the units in terms of the treatment-outcome relationship and estimates a DRF within each of these groups. The Cl-DRF approach aims to fill this need. In Section (ref), we apply the Cl-DRF estimator to the same simulated setting without assuming prior knowledge of the cluster structure. We show that the estimator correctly identifies the true number of clusters, assigns units to clusters with high accuracy, and closely recovers the underlying simulated dose–response relationships.
In the context of policy evaluation, the simplification of treatment to binary terms and the assumption of homogeneous treatment effects across units have long represented the prevailing evaluation framework, aiding in the tractability and interpretation of empirical findings. However, in most circumstances, these simplifications fail to capture the complex nature of policy interventions, which can vary in intensity (e.g., firms receiving public subsidies are generally financed with different amounts), and produce heterogeneous effects even at a fixed treatment intensity, i.e., heterogeneous responses to a homogeneous treatment (e.g., training programs for the unemployed might assist some individuals in finding a job, while potentially disadvantaging others who might have secured a similar job more quickly without the training).
Addressing these limitations, two parallel strands of literature have emerged. The first line of research—reviewed in Section (ref)—models the relationship between continuous treatments and outcomes hirano2004propensity, fong2018covariate, wu2022matching without accounting for the potential heterogeneity of treatment effects. Another strand of literature examines the varied responses to the same treatment across different units athey2016recursive, wager2018estimation, de2020two, callaway2021difference within the binary treatment framework and it is reviewed in Section (ref).
While dichotomizing a continuous treatment variable is known to lead to information loss and potentially compromise data analysis quality fong2018covariate, managing the inherent complexity of a continuous treatment framework poses significant challenges, especially in observational studies where the continuous treatment is not randomly assigned. A common approach to reducing confounding by observed variables is to use the propensity score generalized for the setting of continuous treatments hirano2004propensity. This method aims to provide an unbiased estimate of the ADRF, which relates changes in the outcome variable to changes in treatment for the average unit. Each point of the ADRF represents the expected value of the potential outcome variable at a given treatment level of interest. The estimation process can be split in two steps:
The difficulty of correctly specifying a conditional distribution is exacerbated by the increased dimension of the pre-treatment covariates to be controlled for, that is the confounders. To tackle this issue, wu2022matching adapt nonparametric matching techniques to the context of continuous treatment. One advantage of this approach is that it clearly separates the design phase from the analysis phase, thereby preserving objectivity because no post-treatment information is used during design. Additionally, as misspecification of a conditional density model can be difficult to diagnose and assess, several studies have focused on directly estimating weights to reduce the correlation between marginal moments of (pre-treatment) covariates and the treatment. Approaches along this line of work include, among others, the generalized covariate balancing propensity score approach fong2018covariate, covariate association eliminating weights yiu2018covariate, entropy balancing weights tubbicke2022entropy and independence weights huling2023independence. These approaches prioritize the direct estimation of weights over the explicit estimation and subsequent inversion of a conditional density, and they generally demonstrate superior empirical effectiveness compared to directly modeling the GPS huling2023independence, particularly in small samples cork2025methods.
All these estimators rely on the conditional independence and the positivity assumptions (see Section (ref) for more details). Given the stringency of these assumptions, the use of such estimators to recover the population ADRF should be restricted to settings in which treatment intensity can be regarded as quasi-randomly assigned, conditional on a set of pre-treatment observable characteristics, and in which each unit can plausibly be exposed to a wide range of treatment intensities.
It is well-established that units tend to respond differently to binary treatments athey2016recursive. This highlights the increasing need for estimators that can assess heterogeneous treatment effects (HTE), focusing not just on average treatment effects but on capturing the diverse responses of different units to a binary treatment. Traditionally, HTE estimation involves repeating the analysis for different subgroups expected to exhibit heterogeneous effects based on prevailing theories (e.g., dividing between young and old individuals or between developed and underdeveloped geographical areas). However, this approach could lead to cherry-picking and multiple testing issues athey2016recursive. In addition, it might miss the identification of the most relevant subgroups of units wager2018estimation. Recently, the emphasis has shifted towards data-driven estimation of HTE, with an approach that adopts an algorithm to determine relevant features for estimating treatment effects through the application of machine learning techniques. The goal of this approach is to identify how the effect of a binary treatment varies across a population. These estimators adapt the objective function of well-known machine learning methods, such as decision trees or random forests, with the goal of exploring the heterogeneity of treatment effects relative to a certain set of covariates.
For example, the causal trees introduced by athey2016recursive are decision trees tailored to uncover the heterogeneity of treatment effects. In a causal tree, the aim is to achieve the best prediction of the treatment effect and to do this, the algorithm divides the data to minimize the heterogeneity of treatment effects within the leaves (i.e., the differences in potential outcomes), rather than minimizing the heterogeneity of observed outcomes within the leaves. A more complex approach is represented by the Generalized Random Forest (GRF) developed by athey2019grf. The GRF iterates over random subsets of the data, constructing a causal tree for each subset. By combining numerous trees, the GRF produces robust and consistent treatment effect estimates under the conditional independence assumption.\footnote{Although GRF can be adapted to continuous‐treatment settings, it estimates average partial effects rather than DRFs and requires untreated units for the analysis.}
In summary, both strands of the literature have produced effective estimators within their respective domains. However, they fall short in observational settings with continuous treatments and heterogeneous unit effects, even at the same treatment intensity. Applying these estimators in such contexts could lead to oversimplified policy recommendations, obscuring the complex impacts of policy measures. In the following section, we introduce the Cl-DRF estimator, a method expressly developed for this evaluation context.
The potential outcomes model for observational studies rubin1974estimating is typically employed in binary contexts. Let $t_i\in \{0,1\}$ be the binary treatment variable for unit $i$, and $y_i$ be the outcome variable for unit $i$. Here, $y_i(1)$ represents the potential outcome under treatment, while $y_i(0)$ denotes the potential outcome in the absence of treatment. Consequently, $y_i(1) - y_i(0)$ can be defined as the unit-level causal effect of a binary treatment. In a more general scenario with continuous treatment, $t_i$ assumes values within a real interval $\mathcal{T} = [t_{\min}, t_{\max}]$. Therefore, $y_i(t)$, where $t \in \mathcal{T}$, signifies the outcomes for unit $i$ when exposed to treatment level $t$, varying within the set $\mathcal{T}$. This is often referred to as the unit-dose response tubbicke2022entropy. So, in a continuous setting, we are interested in investigating the trajectory illustrating how $y_i(t)$ changes across relevant values of $t$. Usually, we are interested in population causal effects, so we consider the expected potential outcomes,
also called the average dose-response function (ADRF). Figure (ref) shows the representation of the potential outcomes path for each unit with grey lines in a simulated example in which we know the actual dose-response function of each unit.\footnote{In practice, we do not observe the potential outcomes path for each unit but only the single value $y_i(t_i)$ corresponding to the observed outcome $y_i$ at the realized treatment intensity $t_i$.} The vertical lines at a generic level $t$ allow us to approximate the ADRF for each level $t$. Thus, the population ADRF is the solid black line.
Under random assignment of treatment intensity, differences in average outcomes across treatment levels identify the population ADRF. In contrast, most empirical analyses in the social sciences rely on observational data, where treatment is not randomly assigned. In such studies, treatment intensity typically depends on a set of observed covariates, so careful adjustment for confounding is essential. Therefore, we define the vector of observed pre-treatment covariates for each unit $i$ as $X_i$. With observational studies, we need to rely on some strong assumptions to obtain consistent estimates of the ADRF. Most available estimators for the ADRF in observational settings are based on the same set of assumptions. We refer to A1--A3 as the classical assumptions for ADRF identification.
Assumption A1 (Consistency). For each unit $i$ and treatment level $t$, if $T_i = t$ then $Y_i = Y_i(t)$.
According to the consistency assumption the potential outcome of each unit depends solely on the treatment level it receives. In particular, this assumption implies that potential outcomes do not depend on the treatment assigned to other units, that is, there is no interference cox1958planning, nor on how the treatment is delivered, provided the treatment level is the same, that is, there are no multiple versions of the treatment. This set of conditions is commonly referred to as the Stable Unit Treatment Value Assumption (SUTVA, rubin1980).
Assumption A2 (Weak unconfoundedness). The weak unconfoundedness assumption\footnote{Unlike the unconfoundedness assumption, weak unconfoundedness does not require the joint independence of all potential outcomes; instead, it requires conditional independence to hold at given treatment levels becker2012too. However, in the context under analysis, the difference between the weak unconfoundedness and the unconfoundedness assumptions is more formal than empirically relevant.} implies that $Y_i(t) \perp T_i \mid X_i$, i.e., potential outcomes are independent of treatment conditional on covariates.
The weak unconfoundedness assumption states that once we control for observed covariates, how much treatment a unit receives provides no additional information about what its outcome would be at any specific dose. As argued by bia2012assessing, such assumption is not testable, difficult to satisfy, and not generally applicable.
Assumption A3 (Positivity). For every value of $X_i$, the conditional density of receiving any treatment level $t \in \mathcal{T}$ is bounded away from zero, \[ f_{T \mid X}(t \mid X_i) \geq v > 0 . \]
Also the positivity assumption is crucial for the estimation of causal effects because it guarantees that there is sufficient overlap in the covariate distributions across different levels of the treatment, allowing for meaningful comparisons and generalization of causal effects across the population. In other words, the positivity assumption requires that, for every combination of covariates observed in the data, the conditional density of receiving any treatment level of interest is strictly positive. Such an assumption is arguably less tenable when the treatment is continuous because it is unlikely that every unit has a non-zero conditional density at every treatment value branson2023causal. Violations of the positivity assumption arise when certain treatment levels are rarely or never observed within specific regions of the covariate space. This lack of overlap can lead to extrapolation and a biased estimate of the ADRF. In observational studies, researchers must be particularly vigilant for areas where overlap might be insufficient. branson2023causal introduce a framework for trimming and smoothing techniques that allows for robust estimation even in the presence of limited overlap.\footnote{However, even this approach cannot fully compensate for fundamental violations of the positivity assumption. Furthermore, even if positivity holds, propensity scores close to zero adversely affect the bias and variance of most estimators li2018balancing. See also zhang2025doubly, who show that propensity score weighting and doubly robust estimators can be inconsistent in the absence of positivity, even when the average dose–response function remains identifiable under strong structural assumptions on the potential outcomes. } Moreover, recent works have explored alternative estimands and novel identification strategies that circumvent the need for strict positivity by trading off either the scope of the causal estimand or the strength of the underlying identifying assumptions.\footnote{In particular, schindl2024incremental introduce incremental effects, an estimand defined under stochastic policy interventions that marginally tilt the observed treatment distribution while preserving its support. While a dose–response function characterizes how outcomes vary across treatment intensities, incremental effects collapse this information into a single average effect of a support-preserving policy shift, abstracting from dose-specific heterogeneity in the identified estimand. Complementarily, zhang2024nonparametric retain the ADRF as the target estimand and show that it can be identified without global positivity by leveraging local derivative information, at the cost of imposing stronger structural assumptions on the confounding process.}
As shown above, conventional approaches to continuous treatments adopt a population-level identification strategy and thus obscure heterogeneity in treatment responses. We propose an alternative strategy based on local identification within causal regimes, yielding empirically credible and policy-relevant dose–response relationships.
Specifically, if the sample is partitioned into $C \geq 2$ clusters, either based on policy considerations, such as firm size categories (small, medium and large), or using a data-driven procedure, the target estimands become $C$ cluster-specific DRFs. These estimands remain policy-relevant and identifiable under weaker conditions, as weak unconfoundedness and positivity are only required to hold within clusters rather than in the full population (see the remainder of this section). More precisely, we target a cluster-specific DRF, which constitutes the primary causal estimand of the proposed framework.
This definition makes explicit that the estimand of interest is a subpopulation-specific causal response, rather than a population-average quantity. The Cl-DRF therefore allows treatment effects to vary systematically across clusters, while retaining a clear causal interpretation within each cluster.
Identification and estimation are conducted at the cluster level under the assumptions A1 and B1–B3.
\noindentAssumption B1 (Existence of causal regimes). The population can be partitioned into $C\geq 2$ clusters indexed by $C_i\in\{1,\dots,C\}$ such that a cluster-specific dose--response function $\mu_c(t)=\mathbb E[Y_i(t)\mid C_i=c]$ exists and is well defined for all $t\in\mathcal T_c$. Here, $\mathcal T_c$ denotes the support of the treatment distribution within cluster $c$.
Assumption B2 (Weak unconfoundedness within clusters). For each cluster $c \in \{1,\dots,C\}$ and for all $t \in \mathcal{T}_c$, \[ Y_i(t) \perp T_i \mid X_i,\, C_i = c . \] This condition allows treatment assignment to be confounded at the population level, while requiring conditional independence to hold locally within each causal regime.
Assumption B3 (Positivity within clusters). This assumption rules out extrapolation within clusters and is substantially weaker than requiring overlap over the entire population. For every value of $X_i$ and for each cluster $c$, the conditional density of receiving any treatment level $t \in \mathcal T_c$ is uniformly bounded away from zero, \[ f_{T \mid X, C}(t \mid X_i, C_i = c) \geq v > 0 . \]
We remark that Assumptions B2--B3 can be seen as the cluster-specific counterpart of Assumptions A2--A3. Under the assumptions A1 and B1--B3, the cluster-specific dose--response function is identified. In particular, Assumption B3 ensures that each cluster-specific dose--response function $\mu_c(t)$ is identified only over the treatment support actually observed within cluster $c$, thereby ruling out extrapolation across heterogeneous regimes. The following lemma establishes identification.
Proof of Lemma (ref) is shown in Appendix (ref). While Lemma (ref) provides a valid identification result, direct conditioning on the full covariate vector $X_i$ is not convenient for causal adjustment in continuous-treatment settings. In particular, identification of DRFs requires comparing units at the same treatment level while ensuring covariate balance. To this end, hirano2004propensity introduce the generalized propensity score (GPS) as a balancing score that summarizes all confounding information relevant for treatment assignment at each treatment level.
The cluster-specific GPS summarizes all confounding information contained in the covariates into a scalar function of the treatment level, conditional on cluster membership. This is the natural restriction of the GPS to cluster $c$. The usefulness of the cluster-specific GPS is based on its balancing property within cluster. Indeed, it directly follows from hirano2004propensity that, within each cluster, conditioning on the GPS suffices to control for confounding in the treatment assignment.
Proof of Lemma (ref) is shown in Appendix (ref). This representation provides the basis for estimation strategies that condition on the treatment level and the cluster-specific GPS, rather than on the full set of covariates.
\color{black}
Estimation of the cluster-specific DRFs requires specifying how cluster membership is handled in practice. Two conceptually distinct approaches are in principle available:
We focus on the latter case by employing a cluster-wise extension of the GPS approach of hirano2004propensity\footnote{We recall that the idea behind the GPS approach by hirano2004propensity is to predict missing potential outcomes at specific values of $t$ by fitting a parametric model of the outcome that includes both the treatment and the GPS as covariates, where the GPS is the conditional density of the treatment given a set of covariates.}, building on ideas from grouped and clustered regression models demidenko2018next,sugasawa2021grouped. This procedure allows us to simultaneously estimate cluster membership and cluster-specific DRFs, and will be referred to as the Cl-DRF estimator. In the Cl-DRF framework, units assigned to the same cluster share common parameters governing the treatment--outcome relationship, while clusters are recovered through an iterative procedure.
Following hirano2004propensity, we posit a parametric working model for the treatment assignment mechanism. However, cluster membership is unknown and the treatment assignment model is allowed to vary across clusters. Specifically, we assume a cluster-specific Gaussian specification,
which induces a cluster-specific GPS. Evaluated at the observed treatment level, the GPS is given by
In the second stage, we model the conditional expectation of $Y_i$ given $T_i$ and the GPS. Specifically, we posit a parametric working model for
where the functional form $f(\cdot)$ is common across clusters, while the parameter vector $\boldsymbol{\alpha}_c$ is cluster-specific. Without loss of generality, we can consider a linear model,
though higher-order terms and interactions can be incorporated without altering the structure of the estimator. Because the GPS depends on cluster membership, the regressors entering the outcome regression are cluster-specific.
Motivated by the identification result in Lemma (ref), we model the conditional expectation of the outcome given the treatment level and the cluster-specific GPS. Within each cluster, conditioning on $(T_i,R_{i,c})$ is sufficient for causal adjustment, and the resulting regression can therefore be used to estimate the cluster-specific DRFs. In particular, we remark that, within each cluster, the GPS plays the same role as in the standard framework of hirano2004propensity, thereby justifying this regression-based estimation strategy. However, the parameters $\boldsymbol{\alpha}_c$ could be estimated by ordinary least squares conditional on cluster membership and on the GPS implied by the treatment assignment model. Indeed, in the present setting cluster membership is unknown and the outcome regression is allowed to vary across clusters. As a result, estimation does not reduce to a single linear regression, but instead takes the form of a cluster-wise linear regression problem spath1979algorithmus,hathaway1993switching,park2017algorithms.
To formalize this structure, let $\kappa_{i,c}$ denote cluster membership indicators, with $\kappa_{i,c}=1$ if unit $i$ belongs to cluster $c$ and $\kappa_{i,c}=0$ otherwise. Given a partition of the sample, estimation of the outcome regression parameters reduces to least-squares fitting within each cluster. When cluster membership is unknown, however, the indicators $\{\kappa_{i,c}\}$ constitute additional quantities that must be inferred from the data. Because the outcome regression parameters and the cluster indicators are not jointly identified in closed form under this specification, estimation proceeds by minimizing the following objective function
subject to $\sum_{c=1}^C \kappa_{i,c} = 1$ for all $i$. This formulation naturally leads to an iterative estimation procedure that alternates between updating the outcome regression parameters and updating cluster membership. The resulting block-wise updates admit closed-form expressions, which are summarized in the following proposition.
Taken together, Proposition (ref) characterises a block-coordinate descent procedure in which outcome regression parameters, cluster membership, and the GPS are updated sequentially until convergence. The proof of Proposition (ref) is shown in Appendix (ref).
Taken together, Proposition (ref) and Remark (ref) define an iterative estimation procedure that belongs to the broad class of $k$-means--type and grouped regression algorithms. Theoretical properties of such procedures, including consistency of the estimated partition and of the cluster-specific regression parameters, have been extensively studied in the literature under suitable regularity conditions pollard1981strong,pollard1982central,bonhomme2015grouped,sugasawa2021grouped. We anticipate that the simulation study reported in detail in Appendix (ref) provides strong evidence that the proposed algorithm accurately recovers the cluster structure and the corresponding cluster-specific dose--response functions.
\color{black}
The Cl-DRF estimator builds on the $k$-means regression methodology and therefore it requires specifying the number of clusters $C$ in advance. To this end, we adopt a penalized empirical risk criterion inspired by information criteria commonly used in grouped regression models sugasawa2021grouped. For a given number of clusters $C$, we define
where $c_i$ denotes the cluster assignment of unit $i$ and $\operatorname{dim}(\widehat{\boldsymbol{\alpha}}_c)$ denotes the number of parameters in the cluster-specific regression. Throughout, we set $\zeta_n = \log n$. We notice that the criterion (ref) corresponds to a penalized empirical risk based on the within-cluster sum of squared residuals of the outcome model. The penalty term controls model complexity through the number of clusters and the dimensionality of the cluster-specific regression functions. This construction is analogous to information criteria used in grouped and $k$-means regression settings sugasawa2021grouped, and is therefore referred to as BIC-like.
The criterion (ref) is evaluated for $C = 1,\dots,C_{\max}$, where $C=1$ corresponds to a baseline model without clustering. Since the objective function typically decreases monotonically with $C$, we select the number of clusters using an elbow-type rule. Specifically, we identify the value of $C$ that maximizes the vertical distance between $\operatorname{IC}(C)$ and the straight line connecting $\operatorname{IC}(1)$ and $\operatorname{IC}(C_{\max})$. This procedure captures the point at which further increases in $C$ yield only marginal improvements in fit relative to the increase in model complexity. We emphasize that this selection rule is heuristic in nature and is not claimed to be consistent for the true number of clusters. Nevertheless, the simulation results in Appendix (ref) show that it performs very well in finite samples and accurately recovers the underlying cluster structure across a wide range of scenarios.
Let us consider the motivating example discussed in Section (ref). We now apply the Cl-DRF estimator and the results are shown in Figure (ref). To apply the Cl-DRF estimator we need to choose the number of clusters $C$. Considering different values in the set $C = \{1,2,3,4,5,6,7\}$, we compute the values of the BIC-like criterion as shown in (ref). As shown in Figure (ref), an elbow is detected at $C=4$, which is the true number of clusters. The Rand Index\footnote{The Rand Index measures the similarity between two data clusterings by assessing the proportion of pairs of elements that are consistently grouped together or separately in both partitions rand1971objective. Once a true partition is available, the Rand Index can be used to assess the similarity between the obtained partition with the true one. The index ranges from 0 to 1, where 0 indicates that the clusterings do not agree on any pair of points and 1 indicates that the two partitions under comparison are identical. The higher the index, the better the partition.}, which provides a measure of agreement between the obtained partition with the proposed Cl-DRF estimator and the true one, is equal to 0.98, which suggests an almost perfect identification of the simulated clusters. As shown in Figure (ref), the DRFs reproduced with the Cl-DRF estimator (dashed black line) are very close to the true DRFs (colored lined), demonstrating the efficacy of the Cl-DRF estimator in this example.
In line with the motivating example discussed in Section (ref), we conduct an extensive simulation study to assess the finite-sample performance of the proposed Cl-DRF estimator. The simulations are designed to evaluate two key aspects of the methodology. First, we assess the ability of the information-criterion (ref) to correctly select the number of clusters $C$. Because cluster membership is unobserved, accurate recovery of the underlying partition is a crucial prerequisite for meaningful estimation of cluster-specific dose--response functions. To this end, we compute the Rand Index between the estimated and true partitions across repeated simulations, and we examine the distribution of the selected number of clusters obtained via the IC-based Elbow criterion discussed in Section (ref). Second, conditional on the selected number of clusters, we evaluate the ability of the Cl-DRF estimator to recover the true cluster-specific dose--response functions.
To this aim, we consider two complementary data-generating scenarios. In the first scenario, shown in Appendix (ref), treatment assignment depends on observed covariates through a cluster-specific assignment mechanism, inducing confounding that must be addressed via the GPS. In the second scenario, shown in Appendix (ref), treatment is randomly assigned and therefore independent of covariates. Comparing these settings allows us to isolate the role of the treatment assignment model and to verify that the Cl-DRF estimator behaves appropriately both in the presence and in the absence of confounding. More precisely, we design a simulation study that retains the key structural ingredients of the motivating example—namely, the presence of unobserved cluster-specific treatment-response relationships and heterogeneity in covariate distributions across clusters—while simplifying some components to better isolate the core features of the method. Specifically, the simulation design preserves the violation of the positivity assumption: although treatment is generated conditionally on covariates within each cluster, the supports of the covariates differ across clusters, making some treatment values uncommon or even absent in certain subpopulations. While the positivity violation is present in both settings, in the simulation setting the treatment ranges exhibit more overlap, making cluster boundaries less distinct. This design choice increases the difficulty of the estimation problem, providing a more conservative test of the Cl-DRF estimator’s ability to detect heterogeneity under less favorable, but more realistic, conditions. Another key distinction involves the functional form of the outcome model. In contrast to the motivating example in Section (ref), here we follow the simplified specification introduced in Appendix (ref), where the outcome depends only on the treatment. Covariates therefore induce confounding solely through the assignment mechanism, but do not enter the outcome equation directly. Moreover, we simplify matters by assuming linear outcome models within each cluster.
Overall, the simulation results show that the proposed IC-based procedure selects the correct number of clusters in the vast majority of replications, and that the Cl-DRF estimator recovers the cluster structure with a high Rand Index. Moreover, the estimated dose--response functions closely track the true cluster-specific relationships, supporting the use of the Cl-DRF approach for policy-relevant analysis in settings with clustered heterogeneity.
\color{black}
In this section, we apply the Cl-DRF estimator to assess the causal impact of the European Union’s Cohesion Policy—the world’s largest regional development program in terms of geographic scope, budget, and longevity.
While the bulk of the Cohesion Policy aims at supporting the growth processes of the poorest EU regions, it also finances projects in richer regions. This results in a situation where all regions receive some treatment, but the intensity of treatment is vastly heterogeneous in the amount of funds. For instance, from 1994 to 2006, the North-Holland region received an annual average per capita transfer close to €9, whereas the Região Autónoma dos Açores (PT) received €773, almost 85 times more cerqua2018we. The Cohesion Policy assignment process suggests that regions receiving large amounts of funds are fundamentally different from those receiving few funds. Estimating a single ADRF under these circumstances could result in comparisons between areas with significant disparities in pre-treatment economic growth, per capita GDP, investment rates, and other relevant factors.
Previous literature has already analyzed and demonstrated the presence of heterogeneity of the Cohesion Policy effects with respect to the intensity of treatment becker2012too, cerqua2018we. These studies show that transfers enable faster economic growth in the recipient regions, but the transfer intensity exceeds the aggregate efficiency maximizing level. At the same time, other studies becker2013absorptive, rodriguez2015quality have found that the effect of public transfers on economic growth depends on the absorptive capacity of recipient regions (proxied by human capital endowments and quality of institutions). Considering these findings together, they suggest that one should ideally account for both the heterogeneity of treatment in terms of transfer intensity and the heterogeneous responses to the same amount of funds.\footnote{rodriguez2015quality make an attempt in this direction, but rather than estimating a DRF for different subgroups of regions, they use an approach not grounded in the potential outcomes framework. Indeed, they estimate a simple parametric method in which they model regional economic growth as a function of the amount of funds received, the regional quality of institutions and the interaction among these terms.}
In this analysis we will focus on the relationship between EU funds and economic growth for the period from 1994 to 2015. The main data source is the Annual Regional Database of the European Commission’s Directorate General for Regional and Urban Policy (ARDECO) dataset, which provides comprehensive socio-economic and demographic data for regions across the European Union.
In particular, we will focus on the NUTS-2 regional level and use the GDP per capita adjusted at Purchasing Power Standard (PPS) compound annual growth rate from 1994 to 2015 as the dependent variable.
We will control for the 1994 values of the following variables: i) per capita GDP in PPS; ii) population density; iii) share of gross value added (GVA) in the primary sector; iv) average yearly hours worked per employee; v) average yearly compensation per employee; vi) workplace employment rate. Another key variable for our analysis is the yearly estimate of the amount of funds per capita received over the programming period cerqua2023will that will be used as the treatment intensity variable.
In this section, we estimate the causal relationship between EU funding and economic growth using three approaches for continuous treatments: the GPS of hirano2004propensity, the independence weights of huling2023independence, and the Cl-DRF estimator. The corresponding DRFs are reported in Figure (ref) for the first two methods and Figure (ref) for the Cl-DRF estimator.
Figure (ref) shows a slightly positive relationship between EU funds and economic growth, which flattens out for regions receiving more than 100 euros per capita per year and turns even negative for the most subsidized regions. Not surprisingly, this finding is in line with the findings of becker2012too who have used the GPS estimator to study the same relationship over the two programming periods 1994-1999 and 2000-2006.
On the other hand, looking at Figure (ref) and Table (ref), we see that the Cl-DRF algorithm selects three clusters: Cluster 2 exhibits a relationship between treatment intensity and economic growth similar to that reported in Figure (ref). In contrast, Cluster 1 shows a steep positive relationship between funds and growth, while regions in Cluster 3 display a slight negative relationship. Notably, the common support of the first cluster is much smaller and does not include any of the most subsidized areas. Among the three clusters, regions in Cluster 1 tend to have higher population densities, higher employment rates, greater wealth, and lower subsidy levels. Importantly, the difference between the DRFs reported in Figure (ref) and those in Figure (ref) suggests that the lack of consideration for the likely presence of HTE in the case of continuous treatments might lead to biased estimates and, consequently, to incorrect policy recommendations. For instance, our analysis reveals that, at least for some of the wealthiest EU regions in Cluster 1, there is no evidence supporting the hypothesis that a maximum funding threshold exists beyond which additional resources fail to stimulate further economic growth becker2012too, cerqua2018we. This finding suggests that policymakers can minimize deadweight losses by tailoring subsidy schemes to the heterogeneous marginal returns identified across clusters—avoiding uniform funding caps and aligning allocations with region-specific growth potentials. At the same time, recognizing that some regions cannot fully benefit from additional funding should prompt policymakers to invest in developing those regions’ absorptive capacity.
Looking at the composition of the clusters displayed in Figure (ref), we see that Cluster 1 includes regions from Ireland, Southern Germany, and Southern England, while Cluster 2 predominantly consists of regions in Spain, Northern Germany, the Netherlands, Sweden, and Finland. Lastly, Cluster 3 is composed of many French regions as well as nearly all regions in Italy and Greece. These latter countries were among the hardest hit during the Great Recession, which may explain the counterintuitive finding that higher EU funding led to less growth in these regions. Indeed, particularly during the years 2008-2015, regional inequalities in Greece and Italy widened fingleton2015shocking, despite the larger influx of EU funds to their least developed regions.
Continuous treatments are common in applied work, yet existing methods typically target a single population ADRF. This approach is fully appropriate only when treatment responses are homogeneous across units and overlap is sufficient across treatment intensities. In many observational settings, however, treatment effects and treatment support vary systematically across subpopulations. In such cases, pooling units and estimating a single ADRF generally fails to recover the true relationship, potentially leading to misleading policy conclusions. We address these limitations by proposing the Clustered Dose--Response Function estimator, which replaces population-level identification with cluster-specific DRFs capturing distinct causal response patterns.
By integrating clustering methodologies into the analysis of continuous treatments, this paper expands the methodological toolkit for causal inference in settings characterized by treatment effect heterogeneity and limited overlap. The Cl-DRF estimator provides more detailed and policy-relevant information to decision makers by uncovering distinct causal response patterns that are obscured by population-level analyses. This integrative approach reflects the increasing complexity of policy interventions and responds to the growing demand for more sophisticated tools in policy evaluation, as emphasized in recent contributions to the literature wu2022matching.
The Cl-DRF offers several advantages. First, it relies on weaker versions of the weak unconfoundedness and positivity assumptions by requiring them to hold only within clusters. Second, it enhances interpretability by delivering cluster-specific dose--response relationships. Third, it reduces the scope for specification search and p-hacking by making heterogeneity explicit. Finally, it naturally accommodates data-driven clustering procedures.
We apply the Cl-DRF estimator to evaluate the effects of EU Cohesion Policy in EU-15 NUTS-2 regions. When the clustering structure is ignored, standard estimators suggest a modest positive effect of EU funds on economic growth that plateaus beyond €100 per capita and becomes negative for the most highly funded regions. In contrast, the Cl-DRF uncovers three distinct response patterns: one group exhibits a flat or mildly positive effect, a second displays a strong positive response, and a third shows a slightly negative effect. Ignoring this heterogeneity risks producing biased estimates and misleading policy conclusions. By contrast, our analysis demonstrates that more targeted funding strategies, combined with investments aimed at strengthening regional absorptive capacity, can reduce deadweight losses and enhance policy effectiveness.
This framework opens up several promising directions for future research. In the current implementation, the outcome model is assumed to be parametric and to share the same functional form across clusters. This modeling choice reflects the main focus of the paper, namely the identification of DRFs while allowing for heterogeneity across clusters through cluster-specific estimation of model parameters. While the parameters of the DRF are permitted to vary across clusters, its functional form is specified ex-ante and held fixed. An important avenue for future research is therefore to relax this restriction, either by allowing the functional form itself to differ across clusters or by adopting more flexible, potentially nonparametric, estimation strategies. In addition, the Cl-DRF approach is developed under the working assumption that units can be partitioned into a finite number of clusters characterized by distinct dose--response relationships. While such an assumption is common in the literature on grouped heterogeneity bonhomme2015grouped,sugasawa2021grouped, developing formal statistical tests to assess its empirical validity represents another important direction for future research.
The authors confirm that the data supporting the findings of this study are available within the article or its supplementary materials.