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.
51,047 characters · 12 sections · 55 citation commands
A sliced Wasserstein and diffusion approach to random coefficient models
{5pt} {5pt} {5pt} {5pt} \hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}
{\bf Keywords:} sliced Wasserstein distance, minimum-distance estimation, nearest neighbor methods, diffusion process.
Consider a linear random coefficient model (RCM) involving a set of random variables $(X, Y, \beta)$ with the underlying structure
Here the covariates $X$ and the random coefficients $\beta$ are both supported on $\mathbb{R}^d$. Throughout this paper, we assume that $X$ and $\beta$ are independent. This paper focuses on estimating the distribution of $\beta$ based solely on the observations $\{(X_i, Y_i); i \in [n]\}$, where $[n] := \{1, 2, \ldots, n\}$ and $\{(X_i, Y_i, \beta_i); i \in [n]\}$ are $n$ independent realizations of $(X, Y, \beta)$.
Linear RCMs of the form (ref) play a fundamental role in many quantitative studies; see, for instance, lewbel2017unobserved and Bonhomme2024 for some recent reviews on the applications of RCMs in economics. Traditionally, research on RCMs has primarily focused on estimating the mean and variance of $\beta$ rubin1950note, hildreth1968some, swamy1970efficient. While these two moments summarize average marginal effects and variances, they often fail to capture all information of interest to researchers.
To address this limitation, Beran and Hall beran92 pioneered the study of nonparametrically identifying and estimating the distribution of $\beta$. Since then, various approaches have been proposed, including minimum-distance estimation Beran76,beran1993semiparametric, series expansion methods for (conditional) density functions beran1996nonparametric, Hohmann2016, Dunker19, Dunker25, gaillac2021nonparametric, Gaillac22, kernel methods Hoderlein2010, Holzmann20, nearest neighbor (NN) methods holzmann2024multivariate, and discrete-grid-based constrained least squares methods fox2011, fox2016simple, heiss2022nonparametric.
In this paper, we revisit the minimum-distance estimators proposed by Beran and Millar Beran76 and enhance their framework by incorporating nearest neighbor (NN) methods lin2023estimation,holzmann2024multivariate and the sliced Wasserstein (SW) metric, the latter of which has recently gained prominence in generative artificial intelligence. Compared to generative modeling approaches based on the original Wasserstein metric arjovsky2017wasserstein, tolstikhin2017wasserstein, bousquet2017optimal, the SW metric has been observed to provide similar empirical performance while being computationally more attractive and less sensitive to the dimension of the data kolouri2018sliced, liutkus2019sliced,fan2024minimum. Computationally, optimization algorithms towards minimizing the SW metric have also been introduced and validated in bonneel2015sliced and tanguy2024properties, among many others.
Building on these insights, this paper proposes a novel SW-based minimum-distance estimator for the distribution of $\beta$, where weights are adaptively chosen over an NN graph to reduce bias. To the best of our knowledge, this formulation is new. Computationally, following the framework of tanguy2024properties, we develop block gradient descent algorithms that achieve efficient optimization with a time complexity that scales polynomially in $d$, whereas existing methods typically require an exponentially large number of grid points in $d$. We further establish theoretical guarantees for the proposed estimator and derive its rates of convergence. Lastly, drawing inspiration from parallel research on SW flows bonnotte2013unidimensional, liutkus2019sliced and causal random coefficient models heckman1997making, we bridge our formulation with diffusion processes and causal inference. For the former, we introduce diffusion process-based methods to approximate the distribution of $\beta$. For the latter, we incorporate entropic regularization into the model of heckman1997making, enhancing the identification and estimation of the treatment effect distribution.
{\bf Paper organization.} Section (ref) presents the proposed minimum-distance estimator and Section (ref) illustrates the (approximate) block coordinate decent algorithms for implementation. Section (ref) presents the main theory concerning estimation and computation. Sections (ref) and (ref) connect the studied problem to diffusion processes and casual inference, respectively. The empirical performance of the algorithms is reported in Section (ref). Proofs of all theoretical results are relegated to the Appendix.
{\bf Notation.} Let $(\mathbb{R}^{d}, \|\cdot\|_2)$ represent the $d$-dimensional real space coupled with the Euclidean metric. Write $\mathbb{S}^{d-1} = \{V \in \mathbb{R}^{d}\,\vert\, \Vert V \Vert_{2}=1\}$, $B_{r}(a) = \{ U\in \mathbb{R}^{d} \,\vert \,\Vert U - a\Vert_{2}< r \}$, and $\widebar{B}_{r}(a)$ to represent the sphere, the radius-$r$ open ball, and the corresponding closed ball in $(\mathbb{R},\|\cdot\|_2)$, respectively. Let $\{e_{1},\dots, e_{d}\}$ be the standard bases in $\mathbb{R}^{d}$ with regard to the Euclidean inner product $\langle\cdot,\cdot\rangle$, and let $\sigma$ be the Haar measure in $\mathbb{S}^{d-1}$. For any vectors $v_1,\ldots,v_m\in\mathbb{R}^d$, write $(v_1,\ldots,v_m)=(v_1^\top,\ldots,v_m^\top)^\top$. For vectors $U, \tilde{U}\in \mathbb{R}^{d}$, we also use $P^{U}(\tilde{U}) = \langle U, \tilde{U}\rangle$ to represent their inner product. Let $\mathcal{P}(K)$ denote the set of all probability measures on $K$ and $\delta_{x}\in \mathcal{P}(K)$ be the Dirac measure on $x \in K$. Lastly, for arbitrary measure $\mu \in \mathcal{P}(\mathbb{R}^{d})$, we introduce its projected measure with respect to the direction $V\in\mathbb{S}^{d-1}$ as $ \mu^{V} = \mu \circ (P^{V})^{-1}$.
Consider $\{(X_i, Y_i, \beta_i); i \in [n]\}$ to be $n$ realizations of (ref). The objective of interest is to estimate the distribution of $\beta$, denoted as $\mu_\beta$, based solely on the observations $(X_i, Y_i)$'s.
This paper is concerned with minimum-distance estimation Wolfowitz53, Wolfowitz57, which involves an objective functional and a domain of minimization. In our formulation, the objective functional is based on the Wasserstein $W_2$ distance, defined as
where $\mu$ and $\tilde{\mu}$ are two probability measures on $\mathbb{R}^d$, and $\Pi(\mu, \tilde{\mu})$ denotes the set of all couplings of $\mu$ and $\tilde{\mu}$, i.e., probability measures on $\mathbb{R}^d \times \mathbb{R}^d$ with marginal distributions $\mu$ and $\tilde{\mu}$, respectively.
For each direction $V \in \mathbb{S}^{d-1}$, we introduce a $k$-NN empirical measure $\widebar{\mu}^V_k$ for the objective functional. Specifically, for each $i \in [n]$, first let \[ \tilde{X}_{i} = \frac{X_{i}}{\Vert X_{i} \Vert_{2}}, \quad \tilde{Y}_{i} = \frac{Y_{i}}{\Vert X_{i} \Vert_{2}} \] be the {\it normalized} counterparts of $X_i$ and $Y_i$, respectively. Let ${\mathds 1}(\cdot)$ represent the indicator function and \[ \overline{S}_X(V, k) := \Big\{ j \in [n]: \sum_{i=1}^n {\mathds 1}\big(\| \tilde{X}_i - V \|_2 < \| \tilde{X}_j - V \|_2 \big) < k \Big\} \] denote the set of indices corresponding to the $k$-NNs of $V$ in $\{\tilde{X}_i; i \in [n]\}$. Since $X_i$ may take discrete values, ties can occur, causing the cardinality of $\overline{S}_X(V, k)$ to exceed $k$. In such cases, consistent with lin2024failure, lin2024consistency, we select an {\it arbitrary} subset \[ S_X(V, k) \subset \overline{S}_X(V, k) \quad \text{such that its cardinality}~ \big| S_X(V, k) \big| = k. \] When $\tilde{X}_i$ is continuous, as discussed in Section (ref), we have $S_X(V,k) = \overline{S}_X(V,k)$.
Using this subset and noting that, for any direction $V \in \mathbb{S}^{d-1}$ and $i \in [n]$,
we define the NN-induced random measure as \[ \widebar{\mu}^V_k = \frac{1}{k} \sum_{i \in S_X(V, k)} \delta_{P^V(\tilde{Y}_i \tilde{X}_i)} = \frac{1}{k} \sum_{i \in S_X(V, k)} \delta_{\tilde{Y}_i \tilde{X}_i^\top V}. \] In light of (ref), for arbitrary $V$ and all sufficiently small $k$, the laws of $\widebar{\mu}^V_k$ and $\beta^\top V$ should be close to each other.
Encouraged by the above insight, we introduce the following objective functional
which integrates the squared $W_2$ distance between the NN random measure $\widebar{\mu}^V_k$ and the projected measure $\mu^V$ over all directions $V \in \mathbb{S}^{d-1}$.
Lastly, we propose the following {\it constrained} $k$-NN-based minimum-distance estimator
Here $R$ is a pre-specified radius to enforce a compact space, and the domain $\mathcal{P}_N(\widebar{B}_R(0))$ consists of all size-$N$ discrete measures in $\mathcal{P}(\widebar{B}_R(0))$, defined as
The output is then a discrete approximation to any $\mu_\beta\in \mathcal{P}(\widebar{B}_R(0))$.
The computation of $\hat{\mu}_{\beta}$ in (ref) reduces to finding a size-$N$ set of elements in $\widebar{B}_R(0)$:
Optimizing (ref) can thus also be interpreted as minimizing $\mathcal{F}_{k} \circ \eta$ over $\widebar{B}_{R}(0)^{N}$, which is identifiable up to a permutation. However, solving this minimization problem requires continuous integration over the Haar measure, which is computationally infeasible.
To address this issue, we adopt a common trick and propose a discrete Monte Carlo approximation to the original $\mathcal{F}_k$:
where $V_{1}, \dots, V_{m}$ are sampled independently from the Haar measure $\sigma$, independent of the system.
We accordingly introduce two algorithms for optimizing (ref) using the Monte Carlo approximation (ref), namely, Algorithms (ref) and (ref). Without loss of generality, in these two algorithms we set $N = k$ since otherwise, we could choose the larger of the two. In addition, in these algorithms, we use functions for {\it sorting, argsorting}, and {\it projecting} and define them as follows. For a real sequence $a = (a_{q})_{q\in[k]} = (a_{1}, \dots, a_{k})$, we define \[ \text{sort}(a) = (a_{\nu(q)})_{q\in[k]} = (a_{\nu(1)}, \dots, a_{\nu(k)}) \] and \[ \text{argsort}(a) = (\nu(q))_{q\in[k]} = (\nu(1), \dots, \nu(k)), \] where $(\nu(1), \dots, \nu(k))$ is an arbitrary permutation of $[k]$ such that $a_{\nu(1)} \leq \cdots \leq a_{\nu(k)}$. Additionally, the projection of a vector $\psi \in \mathbb{R}^{kd}$ onto the set $\widebar{B}_{R}(0)^{k}$ in Euclidean space is defined as \[ \text{Proj}_{\widebar{B}_{R}(0)^{k}}(\psi) = \argmin_{w \in \widebar{B}_{R}(0)^{k}} \Vert \psi - w \Vert_{2}. \]
Algorithm (ref) is a block coordinate descent method for computing $\hat{w}$. This algorithm is a revised version of tanguy2024properties, incorporating the $k$-NN graph into the computation. In each iteration, the local optimal transport law is updated using Lemma (ref) in the Appendix, leveraging the vector $\psi_{\ell}$ and the $k$-NN matrix $D$. The vector $\psi_{\ell+1}$ is obtained as the minimizer of the local optimal transport law within the domain $\widebar{B}_{R}(0)^{k}$.
The computation of $\psi_{\ell+1}$, however, involves convex optimization over a constrained domain, which can be computationally expensive. Algorithm (ref) addresses this computational challenge by introducing a projected gradient descent variant. In this algorithm, the optimal transport law from Algorithm (ref) is approximated using an $L^{2}$-distance formulation. This enables the minimization problem to be solved using linear algebra, significantly reducing computational complexity. In addition, it is possible to check that the computational complexity is $O(mdn\log n+tmk \log k)$; see Appendix Section (ref) ahead for details. This is substantially faster than the existing grid-based algorithms that has computational complexity exponentially growing with the dimension.
This section is divided into three parts. The first part discusses the assumptions, while the other two address the rates of convergence of the estimator and the convergence of the proposed algorithms.
We focus on the case $d \geq 2$, as the case $d=1$ is trivial. Define the normalized versions of $X$ and $Y$ as $\tilde{X} = X / \Vert X \Vert_{2}$ and $\tilde{Y} = Y / \Vert X \Vert_{2}$, respectively, such that $\tilde{X} \in \mathbb{S}^{d-1}$. This normalization preserves the inner product structure, as $\tilde{Y} = \langle \beta, \tilde{X} \rangle$.
The first assumption specifies the data-generating structure, which implies that $\{(\tilde{X}_i, \tilde{Y}_i, \beta_i); i \in [n]\}$ are independent and identically distributed as $(\tilde{X}, \tilde{Y}, \beta)$.
The next assumption imposes constraints on the distribution of $\beta$, ensuring that it has a compact support. This setting is similar to those made in Beran76, Hoderlein2010, and Holzmann20, and also common in other nonparametric mixture model analyses Miao24,Han23,lim2024smoothed.
The following two assumptions provide restrictions on the distribution of $\tilde{X}$.
Assumption (ref) is equivalent to Assumption 4 in Hoderlein2010, which naturally holds when $X$ has a full-dimensional support. In particular, this assumption is valid when (ref) does not include an intercept term, i.e., $X_{i,1} = 1$ for all $i \in [n]$. Such scenarios are plausible in many cases, including those discussed in Hoderlein2010 and lewbel2017unobserved. However, it no longer holds if the model includes an intercept term. To address such cases, we introduce the following general assumption.
Assumption (ref) is designed to be weak, although at the expense of some technical complexity. To aid understanding, we present the following result, which simplifies Assumption (ref) to a tail condition on $X$.
Our first main theoretical result addresses the landscape of the objective function and the existence of minimizers to (ref).
The next theorem establishes the rates of convergence of $\hat\mu_\beta$ to $\mu_\beta$ under the $SW_2$ distance. The first part pertains to the case under Assumption (ref), and the second part addresses the case under Assumption (ref).
It is worth noting that Theorem (ref) guarantees the existence of global minimizer(s) for (ref), but it does not provide a general guarantee of uniqueness unless $N = k = n$. Theorem (ref) further demonstrates that any global minimizer of (ref) achieves a polynomial rate of convergence to the true parameter $\mu_\beta$. For related discussions, see Miao24, which observe similar phenomena in nonparametric Poisson mixtures where multiple global minimizers (in that case, nonparametric maximum likelihood estimators) may exist but all converge to the estimand.
At this stage, it remains unclear whether the derived rate is minimax optimal or not. Moreover, as different metrics are employed, the obtained results are not directly comparable to those in the literature. Nevertheless, it is worth highlighting that, unlike prior works including Hoderlein2010, Holzmann20, and Gaillac22, our analysis avoids imposing {\it any} smoothness conditions on the densities of $X$ or $\tilde{X}$, which underscores an advantage of the proposed distance-based approach.
Section (ref) discusses the statistical properties of the estimator $\hat\mu_\beta$ in (ref). However, from a practical standpoint, it is also crucial to understand the procedure from a computational perspective, which is the focus of this section. In this section, the algorithms of interest for approximating $\hat\mu_{\beta}$ are Algorithms (ref) and (ref) presented in Section (ref).
Our first theoretical result in this section examines the landscape of the approximate objective function (ref), paralleling Theorem (ref).
Next, we aim to quantify the algorithmic convergence. To this end, we introduce the following proposition, which analyzes the approximation of $\widebar{\mathcal{F}}_{k}\circ\eta$ to $\mathcal{F}_{k}\circ\eta$.
The main theorem of this section establishes the validity of Algorithms (ref) and (ref).
In this section, building on a recent breakthrough connecting optimal transport-based generative modeling to diffusion processes liutkus2019sliced, we propose a diffusion process-based approximation to the optimization problem (ref). While rigorous theoretical support for this approximation is currently unavailable, its empirical effectiveness has been studied and will be demonstrated in Section (ref).
Specifically, consider the following minimization problem involving a regularized functional $\tilde{\mathcal{F}}_{k}^{\lambda}(\mu)$ over the space of probability measures $(\mathcal{P}(\widebar{B}_{R}(0)), W_{2})$:
where $\mathcal{F}_{k}(\cdot)$ is introduced in (ref), $\lambda > 0$ is a tuning parameter, and $\mathcal{H}(\mu)$ is the regularization functional defined as:
The following proposition provides the theoretical background for formulating the associated continuity equation in this setting. For functions \(f\colon \mathbb{R}^{d} \to \mathbb{R}^{d}\) and \(\tilde{f}\colon \mathbb{R}^{d} \to \mathbb{R}\), define: \[ \textnormal{div} f \coloneqq \sum_{i=1}^{d} \frac{\partial f_{i}}{\partial x_{i}}, \quad \Delta \tilde{f} \coloneqq \sum_{i=1}^{d} \frac{\partial^{2} \tilde{f}}{\partial x_{i}^{2}}. \]
Equation (ref) resembles the Fokker-Planck equation pavliotis2014stochastic for the stochastic differential equation
where \((W_{t})_{t}\) is a standard Brownian motion. We then propose Algorithm (ref) to approximate (ref) using particle systems liutkus2019sliced and Euler-Maruyama discretization pavliotis2014stochastic. For \(L\) particles at step \(\ell\), the iteration is given by:
where \(Z_{\ell h}^{i}\) is standard Gaussian noise. The drift \(v(Q_{\ell h}^{i}, \tilde{\mu}_{\ell h}^{L})\) is computed as:
where \(F_{\tilde{\varrho}_{\ell h}^{V}}\) is the empirical CDF of \((\tilde{\mu}_{\ell h}^{L})^{V}\), and \(F_{\widebar{\mu}^{V}_{k}}^{-1}\) is the empirical quantile function for \(\widebar{\mu}^{V}_{k}\). The algorithm outputs \(Q_{(t+1)h}^{i}\) as samples of the unobserved coefficient. If \(m\) scales polynomially with \(d\), the algorithm has polynomial complexity in \(d\).
In this section, we connect the model (ref) to the causal random coefficient model introduced in heckman1997making:
Here $Z \in \mathbb{R}^d$ represents pretreatment covariates, $w \in \mathcal{W} \subset \mathbb{R}$ is the (continuous or discrete) treatment, $U$ is an unspecified noise term, and $\{Y(w); w \in \mathcal{W} \}$ denotes the potential outcomes corresponding to different treatment levels. It follows that \[ R = Y(w+1) - Y(w), \quad \text{for any } w \in \mathcal{W}. \] Thus, $R$ denotes the unit-level treatment effect. In the special case where $\mathcal{W} = \{0,1\}$, we obtain \[ R = Y(1) - Y(0), \] which represents the unit-level treatment effect of receiving the treatment compared to control.
Notably speaking, we have made two key assumptions in Model (ref). First, we have assumed that the untreated potential outcome satisfies the random coefficient structure: $ Y(0) = \langle Z, \beta_Z \rangle + U$. Second, we have assumed that the unit-level treatment effect $R$ is independent of the observed covariates $Z$, though, it may be correlated with the unobserved noise $U$.
In practice, however, not all potential outcomes $\{Y(w); w \in \mathcal{W} \}$ are observable. In the context of a {\it randomized controlled trial}, we introduce a random treatment variable $W \in \mathcal{W}$ and only observe the outcome $Y : = Y(W)$, which corresponds to the potential outcome under the assigned treatment. The observed outcome follows the model:
Clearly, Model (ref) is a special case of (ref), with
This facilitates the estimation of $R$'s distribution in view of the results established in previous sections. However, to make the arguments formal, we still need to verify Assumption (ref), which is the focus of the following results.
The following proposition addresses the case of a continuous treatment $W$.
This indicates that the Assumption (ref) holds for the case of continuous and unbounded treatment, which aligns with the support assumption regarding the identification in heckman1997making.
However, assuming a continuous and unbounded treatment may be restrictive in practice. To address this, we introduce a regularized working model, analogous to ridge regression in linear models and entropic regularization in optimal transport:
where \[ W_{\epsilon} = W + \tilde{W}_{\epsilon}, \quad \tilde{W}_{\epsilon} \sim \text{Cauchy}(0, \epsilon) \text{ and independent of the system}. \]
As $\epsilon \to 0$, the working model (ref) converges to the original causal random coefficient model (ref), ensuring that the estimated distribution of $R$ remains approximately unchanged. Meanwhile, the additional noise term $\tilde{W}_{\epsilon}$ enhances model flexibility.
The following proposition formalizes this idea.
Proofs are provided in Section (ref) and empirical results for this causal approach are put in Section (ref).
This section evaluates the performance of the algorithms using the folllowing numerical simulations. In the following, covariates are drawn from the von Mises-Fisher distribution with a concentration parameter $0.1$ and a mean direction $(d^{-1/2},\dots,d^{-1/2})$, satisfying Assumption (ref). We consider three types of underlying uniform distributions for the coefficient samples: \[ \mu_{\text{sph}},~~ \mu_{\text{deg}},~~ {\rm and}~~ \mu_{\text{dis}}, \] with supports of different dimensions and $R = 10$ as
The simulations are performed on an 11th Gen Intel Core i7-11700 @ 2.50GHz CPU, and all reproducible codes are available at \url{https://github.com/keunwoolim/rcm_sw}.
The performance of the algorithms is measured by computation time (in seconds) and the empirical distance (in $SW_{2}$ distance) between the original samples and the output, averaged over $100$ experiments. The empirical $SW_{2}$ distance using the Monte Carlo method is computed with the POT package flamary2021pot with $100$ unit vector samples.
For implementing Algorithms (ref) and (ref), we perform $20$ iterations, setting $m = 50$ and $m = 1000$ respectively. We set $k = \lceil{n^{d/(2d-1)}} \rceil$, which corresponds to the general $d\geq 6$ case in Theorem (ref) and is empirically observed to be more robust than the suggested diverging rates in Theorem (ref). For Algorithm (ref), we set $L = k$, $m = 50$, $t = 20$, $h = 1$, and $\lambda = 0.01$. The convex optimization problem for Algorithm (ref) is solved using the CVXPY package diamond2016cvxpy.
Table (ref) reports the results of Algorithms (ref)-(ref) for dimensions $2 \leq d \leq 5$ with $n = 500$ samples from $\mu_{\text{sph}}$. This demonstrates the efficiency of Algorithm (ref), as it is averagely $2.1$ times faster than Algorithm (ref) while maintaining similar empirical accuracy. The approximation accuracy of Algorithm (ref) is lower but not substantially worse than the first two algorithms. This justifies its validity and also aligns with the purpose of generative modeling.
Table (ref) reports the performance of Algorithm (ref) for $n = 500, 1000, 1500, 2000$ samples from $\mu_{\text{sph}}$, $\mu_{\text{deg}}$, and $\mu_{\text{dis}}$, with $d = 2$. For all types of probability measures, the average empirical distance decreases as the sample size increases, while maintaining empirically robust computation time.
Lastly, Figure (ref) plots the coefficient samples and outputs (i.e., $\hat w$'s introduced in (ref)) of Algorithms (ref)-(ref) for $100$ experiments. For Algorithms (ref) and (ref), we performed $20$ iterations with $n = 2000$ and $k = 159$, setting $m = 50$ and $m = 1000$ respectively. For Algorithm (ref), the parameters are $L = 159$, $m = 10$, $t = 20$, $h = 1$, and $\lambda = 0.01$.
Figure (ref) plots the outputs of $100$ experiments of Algorithm (ref) depending on $L$ and $\lambda$, with fixed parameters $n = 2000$, $m = 10$, $t = 20$, and $h = 1$. In Figure (ref), $\lambda$ is fixed as $0.01$, and $L$ varies as $10$, $20$, $40$, and $80$. Meanwhile, in Figure (ref), $L$ is fixed as $20$, and $\lambda$ varies as $0.01$, $0.02$, $0.04$, and $0.08$. Intuitively, the particle systems approximation becomes more precise as $L$ increases, and the generated samples become more diverse as $\lambda$ increases.
Based on the regularized working model (ref) for causal inference, we analyze the AIDS Clinical Trials Group Study 175 (ACTG 175) data provided in the R speff2trial package speff2trial to evaluate the treatment effects of therapy for human immunodeficiency virus type I. Specifically, using the monotherapy of zidovudine (Therapy 1) as the control group, we examine the effects of combination therapies—zidovudine and didanosine (Therapy 2) and zidovudine and zalcitabine (Therapy 3)—as well as monotherapy with didanosine (Therapy 4), as considered in hammer1996trial.
We analyze each investigational treatment (Therapies 2-4) separately against the control. In each analysis, the random treatment variable $W$ is assigned as $W = 0$ for the control group and $W = 1$ for the treatment group. The pretreatment covariates $Z$ include the baseline CD4 T cell count and the number of days of prior antiretroviral therapy, analogous to the setting in 10.1093/biomet/asad045. After standard normalization, this design produces three sets of covariate samples, each with $d = 4$ and respective sample sizes of $n = 1054, 1056, 1093$ for Therapies 2–4. The outcome variable $Y$ is the CD4 T cell count measured at $20 \pm 5$ weeks, which is also normalized.
We implemented a modified version of Algorithm (ref) and conducted 100 experiments for each of the three cases. The parameter $k = \lceil{n^{d/(2d-1)}} \rceil$ was computed as $k = 54, 54, 55$, respectively, and additional parameters were set as $\epsilon = 0.005$, $R = 10$, and $m = 1000$, with 20 iterations.
Figure (ref) presents the box plots for the estimated treatment coefficients corresponding to Therapies 2–4. The estimation results also show the estimated percentiles of the unit treatment effect. For example, the estimated median treatment effects for Therapies 2-4 are $0.113$, $0.047$, and $0.045$, respectively. These results suggest that all three therapies are effective compared to zidovudine monotherapy, with the combination of zidovudine and didanosine showing particularly notable effectiveness. These findings are generally consistent with the results of hammer1996trial.