<figure><img src="image_1.png" /></figure>

HAL
open science

# A Generalized Propensity Score Estimation Methodology for Discrete and Continuous Treatments

Felipe L. Angelim Vieira, Alessandro Leite

<figure><img src="image_2.png" /></figure>

## To cite this version:

Felipe L. Angelim Vieira, Alessandro Leite. A Generalized Propensity Score Estimation Methodology for Discrete and Continuous Treatments. 56es Journées de Statistique de la SFdS, Jun 2025, Marseille, France. hal-05492049

HAL Id: hal-05492049

https://hal.science/hal-05492049v1

Submitted on 3 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A GENERALIZED PROPENSITY SCORE ESTIMATION
METHODOLOGY FOR DISCRETE AND CONTINUOUS
TREATMENTS

Felipe L. Angelim Vieira¹ & Alessandro Leite²

¹Mercado Libre, Brazil, felipeangelim@gmail.com

²INSA Rouen Normandie, Normandie Univ, LITIS UR 4108, 76000 Rouen, France

**Résumé.** Les scores de propension (SP) constituent une stratégie efficace pour contrôler les biais lors de l'estimation des effets causaux à partir de données observationnelles. Néanmoins, les méthodes traditionnelles rencontrent souvent des difficultés avec des structures de traitements complexes, allant au-delà des traitements binaires ou discrets. Dans ce contexte, ce travail propose une méthodologie d'estimation des scores de propension qui prend en compte des structures de traitements généralisées, sans imposer d'hypothèses paramétriques restrictives. Nous adoptons un modèle de classification probabiliste pour estimer les poids stabilisés et les scores de propension via une estimation du rapport de densité, rendant notre méthode adaptable à diverses formes de traitements. Les résultats expérimentaux montrent des améliorations notables en termes de précision et de stabilité sur différents ensembles de données synthétiques et semi-synthétiques, même lorsque le nombre de covariables augmente.

**Mots-clés.** Score de propension, Traitements continus, Estimation dose-réponse, Effet moyen du traitement.

**Abstract.** Propensity scores (PS) comprise an effective strategy to control for biases when estimating causal effects based on observational data. Nevertheless, traditional methods often struggle with complex treatment structures beyond binary or discrete treatments. In this context, this work introduces a methodology for propensity score estimation that accommodates generalized treatment structures without imposing restrictive parametric assumptions. We follow a probabilistic classification model to estimate the stabilized weights and the propensity scores via density ratio estimation, making our method adaptable to various treatment forms. Experimental results show notable performance improvements in accuracy and stability across different synthetic and semi-synthetic datasets, even as the number of covariates increases. This suggests that our approach is flexible enough to be considered for challenging real-world applications when estimating causal effects, such as healthcare, policy analysis, and personalized marketing.

**Keywords.** Propensity Score, Continuous treatments, dose-response estimation, average treatment effect
# 1 Introduction

Causal inference is fundamental in both industrial and research applications that rely on observational data to answer counterfactual questions [1, 2]. In the absence of randomized control trials, observational data can provide valuable causal insights, provided that appropriate causal methods are applied [3]. However, treatment assignment in observational data is often influenced by covariates, and failing to account for this bias may lead to misleading conclusions. A widely used approach to mitigate this issue is the propensity score (PS), which represents the probability of assigning a treatment value given the covariates. Propensity scores can be employed to re-weight dataset samples or replace covariates as regressors when estimating response functions. While probabilistic classification models are effective for estimating propensity scores in discrete treatment settings, continuous and mixed treatment regimes typically require additional parametric assumptions.

This work introduces a novel methodology to approximate the PS independent of the treatment structure. The proposed approach can be applied with any binary classification model or integrated as a loss function in neural network architectures, such as VCNet [4]. Through experimental analysis, we demonstrate the efficacy of our method in handling both discrete and continuous treatments, achieving promising results. Furthermore, the generalization capacity of the methodology enhances its applicability for practitioners aiming to build scalable solutions for obtaining causal insights.

# 2 Problem statement

We consider a dataset $\mathcal{D} = \{(x_i, t_i, y_i)\}_{i=1}^n$, consisting of i.i.d. samples from an observational study where treatment assignment depends on covariates. Each sample $(X, T, Y) \in \mathcal{X} \times \mathcal{T} \times \mathcal{Y}$ represents covariates, treatment level, and observed outcome, respectively, drawn from the joint distribution $P_s(Y | X, T)P_s(T | X)P_s(X)$. Using Rubin's potential outcome framework [5], we define $Y^*(t)$ as the hypothetical outcome if treatment $t$ were applied, irrespective of the actual treatment received. However, only $y_i = Y^*(t_i)$ is observed for each sample, making it impossible to directly estimate $Y^*(t')$ for $t' \neq t_i$. Additionally, observational studies face bias due to treatment assignments influenced by $x_i$, necessitating methods to address these biases. The goal is to estimate the average dose-response function (ADRF), $\mu(t)$, which captures how the average outcome varies with treatment levels $t \in \mathcal{T}$. It is given by:

$$ \mu(t) = \mathbb{E}_{Y' \sim P_s(Y^*(t))}[Y'] = \mathbb{E}_{X \sim P_s(X)} \mathbb{E}_{Y' \sim P_s(Y^*(t)|X)}[Y'] $$

under the positivity and ignorability standard causal inference assumptions.

# 3 Algorithm for estimating propensity score

The core idea of the proposed methodology is to build a synthetic binary classification dataset and to obtain estimates for the density ratio, which can be equal to the generalized propensity score (GPS) or stabilized weights depending on the construction of the dataset. It is deeply
connected to covariate shift methods in machine learning that use classifiers to approximate importance sampling weights in which a model is trained to discriminate between the source and the target distributions $P_s$ and $P_t$ and use the estimated probabilities to approximate the density ratio.

**Lemma 1** Let $n_t, n_s \in \mathbb{N}$. Consider the joint distribution of features $Z \in \mathcal{Z}$ and labels $Q(Z, C) = Q(Z|C)Q(C)$, where $C \sim \text{Bern}(\frac{n_t}{n_t+n_s})$. In addition, let class 0 represent the source distribution $Q(Z|C=0) = P_s(Z)$ and class 1 represent the target $Q(Z|C=1) = P_t(Z)$. The density ratio $\rho(Z) = \frac{P_t(Z)}{P_s(Z)}$ equals the ratio of conditional class probabilities, adjusted by the prior ratio:

$$ \rho(Z) = \frac{Q(C = 1 | Z) n_s}{Q(C = 0 | Z) n_t} \qquad (1) $$

In the case of causal inference, the same methodology can be applied by considering the distribution of the observational data set as a source and as a target, a distribution where the assignment of treatment is independent of the covariates. Letting $\mathcal{Z} = \mathcal{X} \times \mathcal{T}$, the following corollary of Lemma 1 shows how forcing $P_t$ to have independence of treatment and covariates leads to a useful ratio of probabilities of receiving a given treatment.

**Corollary 2** Let $X \in \mathcal{X}$ be covariates, $T \in \mathcal{T}$ a continuous treatment, $P_s(Z) = P_s(X, T) = P_s(T|X)P_s(X)$ the data generating process of an observational experiment and $P_t(Z) = P_s(X)P_{cf}(T)$. Then the density ratio $\rho(X, T)$ is:

$$ \rho(X, T) = \frac{P_{cf}(T)}{P_s(T|X)} $$

Furthermore, creating a binary classification dataset according to Lemma 1 leads to:

$$ \rho(X, T) = \frac{P_{cf}(T)}{P_s(T|X)} = \frac{Q(C = 1|X, T) n_s}{Q(C = 0|X, T) n_t} \qquad (2) $$

In practice, the treatment probability $P_{cf}(T)$ can be set to any probability we can sample from. The Theorem 3 states how specific choices of $P_{cf}$ lead us to the GPS and stabilized weights.

**Theorem 3** If $P_{cf}(T) = \mathcal{U}(\mathcal{T})$ and we account for the volume of treatment space, then the density ratio estimates the inverse of the propensity score:

$$ \rho(X, T) = \frac{Q(C = 1|X, T) n_s}{Q(C = 0|X, T) n_t} = \frac{\mathcal{U}(\mathcal{T})}{P_s(T|X)} = \frac{1}{|\mathcal{T}|} \frac{1}{P_s(T|X)} = \frac{w(X, T)}{|\mathcal{T}|} \qquad (3) $$

In addition, if we let $P_{cf}(T) = P_s(T)$, we obtain the stabilized weights:

$$ \rho(X, T) = \frac{Q(C = 1|X, T) n_s}{Q(C = 0|X, T) n_t} = \frac{P_s(T)}{P_s(T|X)} = w_s(X, T) \qquad (4) $$

While the source distribution corresponds to our observational dataset $\mathcal{D}_s = \{(x_i, t_i)\}$, the target distribution can be constructed synthetically $\mathcal{D}_t = \{(x_i, t'_i) : x_i \sim P_s(X), t'_i \sim P_{cf}(T)\}_{i=1}^{n_r}$. When estimating the propensity score, we assign treatment values sampled from a uniform distribution to samples $x_i$ and, to estimate the stabilized weight, we randomly sample the treatment from the treatment values in $\mathcal{D}_s$. A classifier can then be trained to discriminate whether a sample comes from $\mathcal{D}_s$ or $\mathcal{D}_t$.
### 3.1 Losses for propensity score estimation

One concern when estimating classification probability to obtain density ratio estimation may be the convergence guarantees we have. Fortunately, as demonstrated by Menon and Ong [6], minimizing any Class Probability Estimation (CPE) loss is equivalent to minimizing a Bregman divergence to the true density ratio, with different composite proper losses focusing on different ranges of the density ratio. In other words, the choice of CPE loss dictates a weight function that describes the cost of density ratio error. The least-squares importance fitting (LSIF) loss has the desired property of uniformly weighting the error of density ratio estimation. This loss, when applied to the proposed classification problem in Theorem 3, is defined as:

$$LSIF(\hat{\rho}) = \mathbb{E}_{P_s(X,T)}(\rho(X, T) - \hat{\rho}(X, T))^2 \quad (5)$$

$$= \mathbb{E}_{P_s(X,T)}\hat{\rho}^2(X, T) - 2\mathbb{E}_{P_s(X)P_{cf}(T)}\hat{\rho}(X, T) + C_{\rho} \quad (6)$$

$$= \ell_{lsif}(\hat{\rho}) + C_{\rho} \quad (7)$$

where $\rho$ is the true density ratio and $C_{\rho}$ is a constant that can be ignored when optimizing for $\hat{\rho}$, since:

$$\rho = \arg \min_{\hat{\rho}} LSIF(\hat{\rho}) = \arg \min_{\hat{\rho}} \ell_{lsif}(\hat{\rho}) \quad (8)$$

The term $\ell_{lsif}(\hat{\rho})$ can be approximated using the observational dataset $\mathcal{D}_s$ by:

$$\ell_{lsif}(\hat{\rho}) \approx \frac{1}{|\mathcal{D}_s|} \sum_{(x_i, t_i) \in \mathcal{D}_s} \hat{\rho}^2(x_i, t_i) - 2 \frac{1}{|\mathcal{D}_s|} \sum_{\substack{x_i \in \mathcal{D}_s \\ t'_i \sim P_{cf}(T)}} \hat{\rho}(x_i, t'_i) \quad (9)$$

This loss can be easily incorporated into the training loop of neural network models that estimate balancing weights, such as VCNet [4]. Similarly, the logistic loss can be used to approximate the weights:

$$\ell_{logistic}(\hat{\rho}) \approx -\frac{1}{|\mathcal{D}_s|} \sum_{(x_i, t_i) \in \mathcal{D}_s} \log \hat{\rho}(x_i, t_i) - \frac{1}{|\mathcal{D}_s|} \sum_{\substack{x_i \in \mathcal{D}_s \\ t'_i \sim P_{cf}(T)}} \log (1 - \hat{\rho}(x_i, t_i)) \quad (10)$$

where, $\hat{\rho}(x_i, t_i) := (\hat{\rho}(x_i, t_i) + 1)^{-1}$. When $P_{cf}(T)$ is set to $P_s(T)$ we call stabilized LSIF (LSIF-S) or logistic loss, and uniform when $P_{cf}(T)$ is set to $\mathcal{U}(T)$.

## 4 Experimental Results

We evaluate the effectiveness of the proposed propensity score estimation approach in both continuous and binary treatment settings. To quantify performance, we use the expected mean squared error (EMSE) to compare the estimated function $\hat{\mu}(T)$ with the true function $\mu(T)$. This metric is computed over either the observational treatment distribution $P_s(T)$ or a uniform distribution $\mathcal{U}(T)$, offering insights into model performance across different treatment distributions. For continuous treatments, we employ the GPS method [7] and the VCNet architecture with targeted regularization [4] to estimate the average dose-response function. For binary treatments, we use GPS, standard Inverse Propensity Score Weighting (IPW), and Doubly Robust (DR) estimation. Each method is evaluated using both the proposed approach and baseline propensity score methodologies. In the continuous treatment setting,
we adopt the piecewise linear interpolation method from [4] as a baseline for propensity score estimation, alongside Gradient Boosting and a feed-forward neural network classifier for binary treatment datasets. Our implementation employs feed-forward neural networks trained with the losses defined in Equations (9) and (10), in both uniform and stabilized forms. We pre-process treatment variables for continuous treatments using a spline basis of degree 2 with 10 knots uniformly placed between 0.1 and 0.9 instead of using the raw treatment values.<fnref n="1" /> To evaluate the versatility of our methodology, we also use tree-based models, such as the Gradient Boosting Classifier, on synthetic datasets generated according to Equations (3) and (4). A challenge with tree models is their tendency, depending on hyperparameters, to ignore the treatment variable, making predictions solely based on covariates and returning constant class probabilities. If this occurs, we resample treatments from $P_{cf}(T)$ and refit the classifier, repeating up to five iterations until the model produces non-constant probabilities.

**Continuous Treatment:** Table 1 summarizes the results for the continuous treatment settings. Focusing first on neural network PS estimators, our GPS methodology consistently outperforms the baseline in all cases except for LRS in the Synthetic B dataset. For Synthetic A, all neural network-based PS methods achieved at least a 37% reduction in the EMSE compared to the GPS baseline and a 33% reduction relative to the uniform distribution (EMSE Uniform). Regarding the VCNet + TR estimator, the impact of density ratio estimation varies across datasets. For example, in the Synthetic A dataset with the LSIF loss, it worsens performance, while in other cases, it improves results. Notably, LSIFS consistently surpassed the baseline across all scenarios, demonstrating its robustness. For tree-based classifiers, while the standard Gradient Boosting Classifier (GBC) showed suboptimal performance, its stabilized variant improved errors over the baseline in all datasets, highlighting its effectiveness in this context.

**Discrete Treatment:** Table 2 describes the results for the proposed propensity score approach in the discrete treatment setting. Except for LRS and GBC, we observe performance improvements over both the neural network and tree-based baselines. To investigate the high error observed with IPW-BR, we found that the neural network struggles to improve accuracy on the validation set, quickly reaching the patience threshold and returning an underfitted model. In contrast, for the other methods, the loss improves for a longer duration before plateauing, though they eventually reach a similar validation loss level.

<figure><img src="image_4.png" /><figcaption>Figure 1: Results for Synthetic A data with varying feature dimensions and a fixed sample size of 1000 using the GPS estimation procedure.</figcaption></figure>

**Sensitivity to feature dimension:** Figure 1 highlights that although the EMSE increases

<footnote n="1">¹Spline pre-processing allows for greater flexibility in capturing non-linear relationships.</footnote>
<table><caption>Table 1: Comparing methods' performance in the continuous settings</caption><thead><tr><th rowspan="2">PS<br>Method</th><th colspan="2">GPS</th><th colspan="2">VCNet TR</th></tr><tr><th>√EMSE</th><th>√EMSE (Uniform)</th><th>√EMSE</th><th>√EMSE (Uniform)</th></tr></thead><tbody><tr><td colspan="5"><strong>Synthetic A (1000 samples, 10 features)</strong></td></tr><tr><td>BR</td><td>0.2845 ± 0.0469</td><td>0.2725 ± 0.0373</td><td>0.1189 ± 0.0519</td><td><strong>0.1148 ± 0.0556</strong></td></tr><tr><td>LR</td><td>0.1293 ± 0.0277</td><td>0.1282 ± 0.0261</td><td>0.1266 ± 0.0751</td><td><strong>0.1241 ± 0.0674</strong></td></tr><tr><td>LRS</td><td>0.1208 ± 0.0355</td><td><strong>0.1181 ± 0.0319</strong></td><td>0.1458 ± 0.0652</td><td>0.1458 ± 0.0604</td></tr><tr><td>LSIF</td><td><strong>0.1768 ± 0.0543</strong></td><td>0.1807 ± 0.0542</td><td>0.4315 ± 0.1628</td><td>0.4353 ± 0.1626</td></tr><tr><td>LSIF-S</td><td>0.1495 ± 0.0394</td><td>0.1563 ± 0.0419</td><td>0.1061 ± 0.0387</td><td><strong>0.1055 ± 0.0379</strong></td></tr><tr><td>GBC</td><td>0.4384 ± 0.0410</td><td>0.4520 ± 0.0340</td><td>-</td><td>-</td></tr><tr><td>GBC-S</td><td>0.2016 ± 0.0381</td><td>0.2017 ± 0.0353</td><td>-</td><td>-</td></tr><tr><td colspan="5"><strong>Synthetic B [8]</strong></td></tr><tr><td>BR</td><td>0.0955 ± 0.0137</td><td>0.1220 ± 0.0196</td><td>0.1144 ± 0.0405</td><td><strong>0.0781 ± 0.0164</strong></td></tr><tr><td>LR</td><td><strong>0.0577 ± 0.0146</strong></td><td>0.0711 ± 0.0214</td><td>0.0740 ± 0.0184</td><td>0.0698 ± 0.0154</td></tr><tr><td>LRS</td><td>0.1571 ± 0.0445</td><td>0.2423 ± 0.0719</td><td><strong>0.0707 ± 0.0170</strong></td><td>0.0811 ± 0.0265</td></tr><tr><td>LSIF</td><td><strong>0.0598 ± 0.0155</strong></td><td>0.0772 ± 0.0235</td><td>0.0736 ± 0.0185</td><td>0.0784 ± 0.0206</td></tr><tr><td>LSIF-S</td><td>0.0891 ± 0.0267</td><td>0.1181 ± 0.0491</td><td>0.1077 ± 0.0286</td><td><strong>0.0724 ± 0.0183</strong></td></tr><tr><td>GBC</td><td><strong>0.1649 ± 0.0451</strong></td><td>0.2410 ± 0.0697</td><td>-</td><td>-</td></tr><tr><td>GBC-S</td><td><strong>0.0707 ± 0.0195</strong></td><td>0.0912 ± 0.0366</td><td>-</td><td>-</td></tr><tr><td colspan="5"><strong>News (semi-synthetic)</strong></td></tr><tr><td>BR</td><td>0.1762 ± 0.0741</td><td>0.1556 ± 0.0545</td><td>0.1539 ± 0.0614</td><td><strong>0.1346 ± 0.0426</strong></td></tr><tr><td>LR</td><td>0.1320 ± 0.0672</td><td><strong>0.1045 ± 0.0507</strong></td><td>0.1410 ± 0.0583</td><td>0.1285 ± 0.0494</td></tr><tr><td>LRS</td><td>0.1299 ± 0.0703</td><td><strong>0.1026 ± 0.0511</strong></td><td>0.1421 ± 0.0531</td><td>0.1314 ± 0.0407</td></tr><tr><td>LSIF</td><td>0.1220 ± 0.0592</td><td><strong>0.1022 ± 0.0432</strong></td><td>0.1448 ± 0.0601</td><td>0.1392 ± 0.0534</td></tr><tr><td>LSIF-S</td><td>0.1293 ± 0.0714</td><td><strong>0.1016 ± 0.0513</strong></td><td>0.1334 ± 0.0653</td><td>0.1249 ± 0.0649</td></tr><tr><td>GBC</td><td>0.1606 ± 0.0538</td><td><strong>0.1487 ± 0.0450</strong></td><td>-</td><td>-</td></tr><tr><td>GBC-S</td><td>0.1510 ± 0.0689</td><td><strong>0.1182 ± 0.0530</strong></td><td>-</td><td>-</td></tr></tbody></table>

BR: baseline, LR: logistic regression, LRS: logistic regression stabilized, LSIF (our): least-squares importance fitting

LSIF-S (our): LSIF stabilized, GBC-S: Gradient Boosting Classifier (stabilized), GBC (our): Gradient Boosting Classifier

with the dimensionality, all neural network methods consistently outperform the baseline, except for the LSIF loss at 100 features. Logistic losses, in particular, show significantly lower errors compared to the baseline across all feature dimensions. For gradient-boosting trees, the results indicate that while GBC-S performs well in low-dimensional settings, its performance deteriorates as the number of features increases. Additionally, the performance of GBC suggests that it is not well-suited for density ratio estimation in high-dimensional problems.

# 5 Conclusion

This study introduces a novel methodology for estimating propensity scores that supports diverse treatment structures, including continuous and binary regimes. Experimental results demonstrate that the proposed approach consistently outperforms state-of-the-art methods, even in high-dimensional covariate settings. These findings have significant implications for both scientific and industrial applications, such as advancing causal analysis in healthcare studies of drug dosages or improving causal insights from complex observational data in industry. Despite these promising results, several challenges remain. For example, tree-based classifiers may require careful pre-processing of the joint covariate-treatment space to
<table><caption>Table 2: Comparing the different methods' performance on the News 1D-Discrete dataset</caption><thead><tr><th rowspan="2">PS<br>Method</th><th colspan="2">GPS</th><th colspan="2">IPW</th><th colspan="2">DR</th></tr><tr><th>√EMSE</th><th>√EMSE<br>(Uniform)</th><th>√EMSE</th><th>√EMSE<br>(Uniform)</th><th>√EMSE</th><th>√EMSE<br>(Uniform)</th></tr></thead><tbody><tr><td>BR</td><td>0.0975 ± 0.0426</td><td>0.1330 ± 0.0613</td><td>0.5087 ± 0.0260</td><td>0.7156 ± 0.0081</td><td>0.0671 ± 0.0545</td><td><b>0.0441 ± 0.0242</b></td></tr><tr><td>BTC</td><td>0.1366 ± 0.0632</td><td>0.1824 ± 0.0894</td><td>0.3955 ± 0.0528</td><td>0.5559 ± 0.0732</td><td><b>0.0284 ± 0.0110</b></td><td>0.0356 ± 0.0182</td></tr><tr><td>LR</td><td>0.0833 ± 0.0454</td><td>0.1125 ± 0.0646</td><td><b>0.0241 ± 0.0098</b></td><td>0.0260 ± 0.0171</td><td>0.0348 ± 0.0129</td><td>0.0388 ± 0.0218</td></tr><tr><td>LRS</td><td>0.1215 ± 0.0366</td><td>0.1691 ± 0.0515</td><td>0.0459 ± 0.0181</td><td>0.0612 ± 0.0269</td><td><b>0.0373 ± 0.0172</b></td><td>0.0388 ± 0.0229</td></tr><tr><td>LSIF</td><td>0.0865 ± 0.0451</td><td>0.1146 ± 0.0659</td><td><b>0.0269 ± 0.0123</b></td><td>0.0308 ± 0.0207</td><td>0.0349 ± 0.0135</td><td>0.0393 ± 0.0227</td></tr><tr><td>LSIF-S</td><td>0.0870 ± 0.0393</td><td>0.1199 ± 0.0551</td><td><b>0.0279 ± 0.0057</b></td><td>0.0333 ± 0.0101</td><td>0.0340 ± 0.0134</td><td>0.0379 ± 0.0225</td></tr><tr><td>GBC</td><td>0.0663 ± 0.0493</td><td>0.0876 ± 0.0710</td><td>0.1007 ± 0.0609</td><td>0.1371 ± 0.0851</td><td>0.0551 ± 0.0352</td><td><b>0.0458 ± 0.0283</b></td></tr><tr><td>GBC-S</td><td>0.0595 ± 0.0410</td><td>0.0758 ± 0.0611</td><td>0.0434 ± 0.0183</td><td>0.0575 ± 0.0266</td><td>0.0410 ± 0.0233</td><td><b>0.0391 ± 0.0236</b></td></tr></tbody></table>

BR: baseline, LR: logistic regression, LRS: logistic regression stabilized, LSIF (our): least-squares importance fitting, LSIF-S (our): LSIF stabilized, BTC: Baseline Tree classifier, GBC-S: Gradient Boosting Classifier (stabilized), GBC (our): Gradient Boosting Classifier

balance essential interaction modeling and overfitting prevention. Furthermore, the current methodology does not extend to structured data inputs (e.g., images or text) or sequential and multi-dimensional treatments. Future research will focus on addressing these limitations and exploring the applicability of the approach in more complex settings and challenging applications.

## Acknowledgments

This research was partially funded by the French National Research Agency (ANR) (Grant No. ANR-23-CPJ1-0099-01) and the France 2030 program under the reference 23-PEIA-004.

## References

[1] J. Pearl, *Causality: Models, Reasoning and Inference*, 2nd ed. Cambridge University Press, 2009.

[2] J. Pearl and D. Mackenzie, *The book of Why: the new Science of Cause and Effect*. Basic books, 2018.

[3] D. Jacob, “CATE Meets ML: conditional average treatment effect and machine learning,” *Digital Finance*, vol. 3, no. 2, pp. 99–148, 2021.

[4] L. Nie, M. Ye, Q. Liu, and D. Nicolae, “VCNet and functional targeted regularization for learning causal effects of continuous treatments,” in *International Conference on Learning Representations*, 2021.

[5] D. B. Rubin, “Estimating causal effects of treatments in randomized and nonrandomized studies.” *Journal of Educational Psychology*, vol. 66, no. 5, p. 688, 1974.

[6] A. Menon and C. S. Ong, “Linking losses for density ratio and class-probability estimation,” in *International Conference on Machine Learning*, 2016, pp. 304–313.

[7] K. Hirano and G. W. Imbens, “The propensity score with continuous treatments,” *Applied Bayesian modeling and causal inference from incomplete-data perspectives*, vol. 226164, pp. 73–84, 2004.

[8] X. Wang, S. Lyu, X. Wu, T. Wu, and H. Chen, “Generalization bounds for estimating causal effects of continuous treatments,” *Advances in Neural Information Processing Systems*, vol. 35, pp. 8605–8617, 2022.