MOSIC: Model-Agnostic Optimal Subgroup Identification with Multi-Constraint for Improved Reliability

## Abstract

Abstract Current subgroup identification methods typically follows a two-step approach: first estimate conditional average treatment effects (CATEs) and then
apply thresholding or rule-based procedures to define subgroups.
While intuitive,
this decoupled approach fails to incorporate key constraints essential for real-world clinical decision-making—such as subgroup size and propensity overlap. These constraints operate on fundamentally different axes than CATE estimation and are not naturally accommodated within existing frameworks, thereby limiting the practical applicability of these methods.
We propose a unified optimization framework that directly solves the primal constrained optimization problem to identify optimal subgroups.
Our key innovation is a reformulation of the constrained primal problem as an unconstrained differentiable min-max objective , solved via a
gradient descent-ascent algorithm.
We theoretically establish that our solution converges to a feasible and locally optimal solution.
Unlike threshold-based CATE methods that apply constraints as post-hoc filters, our approach enforces them directly during optimization.
The framework is model-agnostic , compatible with a wide range of CATE estimators, and extensible to additional constraints like cost limits or fairness criteria.
Extensive experiments on synthetic and real-world datasets demonstrate its effectiveness in identifying high-benefit subgroups while maintaining better satisfaction of constraints.

## 1 Introduction

In precision medicine, a fundamental challenge in observational studies is identifying patient subgroups that benefit most from specific treatments, where heterogeneous effects must be estimated from non-randomized data  (Kosorok and Laber 2019; Kravitz, Duan, and Braslow 2004).
Most existing methods adopt a two-step paradigm: they first estimate CATEs using machine learning methods (Athey and Imbens 2016; Shalit, Johansson, and Sontag 2017; Shi, Blei, and Veitch 2019), then deriving subgroups through either thresholding (VanderWeele et al. 2019; Lipkovich et al. 2024) or simplified rule-based models (Cai et al. 2022).

However, this two-stage approach falls short when applied to real-world settings, where subgroup identification must account for diverse and interacting constraints.
Clinical deployment often requires satisfying statistical conditions like minimum subgroup size and overlap (VanderWeele et al. 2019; Crump et al. 2009), as well as operational and ethical considerations such as budget limits, safety thresholds, and fairness criteria.
Existing methods typically treat these constraints as post-hoc filters rather than integrating them into the optimization process. As a result, they *struggle to jointly satisfy multiple constraints*, leading to instability and poor performance.
These challenges highlight a deeper disconnect between the continuous nature of CATE estimates and the discrete, constraint-driven structure of clinically actionable subgroup identification. This motivates the need for new frameworks that can incorporate and optimize over multiple real-world constraints in a unified and principled way.

We propose MOSIC (Model-agnostic Optimal Subgroup Identification with multi-Constraints), a *novel optimization framework* that identifies subgroups with maximal CATE while satisfying group size and overlap constraints—with flexibility to incorporate additional constraints. Our approach addresses the key challenge of nonconvex/nonconcave optimization, and our contributions are threefold:

- •
Problem Reformulation for Stable Solutions:
We develop a stable optimization procedure that: (1) formulates the task as a constrained problem
(Section [3.1](https://arxiv.org/html/2504.20908v2#S3.SS1)),
(2) absorbs constraints into the objective via a reformulation
(Section [3.2](https://arxiv.org/html/2504.20908v2#S3.SS2)), and
(3) solves it using a gradient descent-ascent (Section [4.1](https://arxiv.org/html/2504.20908v2#S4.SS1)). To improve stability, we modify the objective (Section [4.2](https://arxiv.org/html/2504.20908v2#S4.SS2)) and
establish that the resulting solution is locally optimal and feasible.
- •
Flexibility:
MOSIC maintains unique flexibility across three dimensions: (1) supporting multiple subgroup model architectures (e.g., multilayer perceptrons, decision trees) for different interpretability-performance tradeoffs, (2) compatibility with various ATE estimators, and (3) extensibility to diverse clinical constraints beyond our focus on size and overlap requirements.
- •
Comprehensive Evaluation: We extensively evaluate our framework on both synthetic and real-world data, demonstrating its effectiveness in optimal subgroup identification under multiple constraints (Section [5](https://arxiv.org/html/2504.20908v2#S5)). Our implementation is publicly available at
https://anonymous.4open.science/r/MOSIC3-8F13.

## 2 Related Work

We review treatment effect estimation, overlap handling, subgroup identification, and constrained optimization methods, highlighting that multi-constraints subgroup identification remains under-explored.

##### Treatment Effect Estimation

MOSIC accommodates various average treatment effect (ATE) estimators. Traditional methods like IPTW and meta-learners (e.g., S-, T-, X-learners, Künzel et al. 2019) rely on either the treatment or outcome model, making them sensitive to misspecification. In contrast, doubly robust estimators such as AIPTW (Robins, Rotnitzky, and Zhao 1995), targeted maximum likelihood estimation (TMLE) (Van Der Laan and Rubin 2006), and DML (Chernozhukov et al. 2018) require only one model (treatment or outcome) to be correctly specified. In this paper, we adopt AIPTW for illustration.

While ATE captures population-level effects, CATE estimation enables subgroup-specific recommendations.
Modern methods include 1) tree-based methods like Causal Tree (CT) (Athey and Imbens 2016) and Causal Forest (CF) (Wager and Athey 2018), and 2) neural network-based approaches, such as TARNet and CFR (Shalit, Johansson, and Sontag 2017), Dragonnet (Shi, Blei, and Veitch 2019), and BCAUSS (Tesei et al. 2023). Notably, Dragonnet includes a TMLE-inspired regularization term for enhanced finite-sample stability.
In this work, to construct our AIPTW estimator, we use Dragonnet to estimate the potential outcome models and use LR to estimate the propensity score model.

##### Dealing with Limited Overlap

Limited sample overlap can bias treatment effect estimates or inflate variance
(Crump et al. 2009).
Common solutions include truncating propensity scores (Gruber et al. 2022; Cole and Hernán 2008) and excluding
low-overlap units
(Crump et al. 2009; Schweisthal et al. 2024; Kallus 2020).
To ensure reliable ATE estimation
on the identified subgroup,
we incorporate a set of constraints to avoid regions with limited overlap.

##### Optimal Subgroup Identification

Existing approaches fall into three concise categories: (i) *baseline methods* without interpretability or constraints, (ii) methods emphasizing *interpretability*, and (iii) methods designed to incorporate domain-specific *constraints*.

Common
subgroup identification methods either: (1)
rank patients by estimated CATE values (Cai et al. 2011; VanderWeele et al. 2019) (benchmarked this approach employing CT, CF, and Dragonnet as
CATE estimators in
Section [5](https://arxiv.org/html/2504.20908v2#S5)), or (2)
optimize
individual treatment
rules using methods like
Outcome-Weighted Learning (OWL) (Zhao et al. 2012; Liu et al. 2018),
which frames treatment as a weighted classification problem.
While providing useful baselines, these methods do not readily accommodate additional constraints, a limitation our approach addresses.

For *interpretability*,
*Decision Tree* (DT)-based methods are widely used (Foster, Taylor, and Ruberg 2011; Lipkovich et al. 2011; Dusseldorp, Doove, and Mechelen 2016; Huang et al. 2017; Athey and Wager 2021).
A representative example is
Virtual Twins (VT) (Foster, Taylor, and Ruberg 2011). It first estimates CATE and then applies a DT for interpretable subgroup identification.
Beyond trees, rule learning approaches have been adopted (Wang and Rudin 2021; Zhou et al. 2024).
However, these methods rely on combinatorial searches and do not scale. Our method can instead leverage DTs as the backbone model, achieving the same level of interpretability while remaining scalable. We compare its performance against other DT-based methods in Section [5.4](https://arxiv.org/html/2504.20908v2#S5.SS4.SSSx1).

While single *constraints* like minimum group size can be applied post hoc to CATE estimates (VanderWeele et al. 2019),
multi-constraint problems often require
directly optimization
to yield better outcomes.
CAPITAL (Cai et al. 2022) is the most closely related approach to ours:
it
maximize subgroup size
under an average CATE threshold
and allow extension to a single constraint using
Lagrangian relaxation. However, it
struggles with multiple constraints
due to instability and lack of feasibility guarantees.
In Section [5](https://arxiv.org/html/2504.20908v2#S5), we compare our approach with VT, OWL, and CAPITAL.

##### Constrained Optimization

Constrained optimization algorithms
vary by problem convexity. For *convex objective and constraints*,
classical methods such as Projected Gradient Descent (PGD), Frank-Wolfe (FW), Interior Point Methods (IPM), and Lagrangian-based methods (e.g., Alternating Direction Method of Multipliers, ADMM, Boyd et al. 2011) are effective.
For *non-convex objectives with convex feasible regions*, global optimality is NP-hard and convergence guarantees weaken
(Lacoste-Julien 2016; Wang, Yin, and Zeng 2019).
Our setting—*non-convex objective with non-convex constraints*—poses greater challenges:
PGD struggles with complex prediction, IPM scales poorly with constraint count, FW assumes convexity, and Lagrangian methods often lack stability and feasibility guarantees. ADMM fails here due to the non-separable objective. While ADMM variants exist for structured problems (Gao, Goldfarb, and Curtis 2020), none apply to our setting.
In contrast, our method guarantees both constraint feasibility and local optimality, critical for real-world deployment.

## 3 Problem Settings

Let $\boldsymbol{X}\in\mathcal{X}\subseteq\mathds{R}^{d}$ denote baseline covariates, $A\in\mathcal{A}=\{0,1\}$ the treatment (0: control, 1: treatment), and $Y\in\mathcal{Y}$ the outcome, which may be binary ($\mathcal{Y}={0,1}$) or continuous ($\mathcal{Y}\subseteq\mathds{R}$). The observational dataset consists of $n$ samples ${(\boldsymbol{x}i,a_{i},y_{i})}_{i=1}^{n}$. Let $Y(0)$ and $Y(1)$ denote the potential outcomes under control and treatment. The propensity score is denoted as $e(\boldsymbol{x})=P(A=1|\boldsymbol{X}=\boldsymbol{x})$, and potential outcomes as $\mu_{a}(\boldsymbol{x})=\mathds{E}[Y|\boldsymbol{X}=\boldsymbol{x},A=a]$.

We adopt standard causal inference assumptions: 1) Stable Unit Treatment Value Assumption (SUTVA): $Y=A\cdot Y(1)+(1-A)\cdot Y(0)$; 2) Unconfoundedness: $\{Y(0),Y(1)\}\perp A|\boldsymbol{X}$; 3) Overlap: $0<e(\boldsymbol{x})<1,\forall\boldsymbol{x}\in\mathcal{X}$.

### 3.1 Reliable Subgroup Identification Framework

In this framework, our goal is to identify a patient subgroup with the largest ATE while ensuring reliability. We achieve this by learning a subgroup identification model $\tilde{S}:\mathds{R}^{d}\mapsto\{0,1\}$, which assigns a patient with covariates $\boldsymbol{X}=\boldsymbol{x}$ to the subgroup ($\tilde{S}(\boldsymbol{x})=1$) or excludes them ($\tilde{S}(\boldsymbol{x})=0$). We define the subgroup’s ATE as:

$$ $\displaystyle\mathds{E}[Y(1)-Y(0)|\tilde{S}(\boldsymbol{X})=1]=$ $\displaystyle\mathds{E}\left[\mathds{E}[Y|A=1,\boldsymbol{X}]-\mathds{E}[Y|A=0,\boldsymbol{X}]|\tilde{S}(\boldsymbol{X})=1\right],$ (1) $$

following from the assumptions of SUTVA and unconfoundedness, and the inner expectations in Equation ([3.1](https://arxiv.org/html/2504.20908v2#S3.Ex1)) represent the pseudo-outcomes on each sample.

Let $\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})$ denote the estimated pseudo-outcomes for each sample, a general estimator for the subgroup’s ATE, can then be written as:

$$ $\frac{\sum_{i=1}^{n}\mathds{1}(\tilde{S}(\boldsymbol{X})=1)\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})}{\sum_{i=1}^{n}\mathds{1}(\tilde{S}(\boldsymbol{X})=1)},$ $$

where $\mathds{1}(\cdot)$ is the indicator function.

In addition to maximizing the treatment effect on the identified subgroup, we introduce the following constraints:

- •
Minimum subgroup size requirement. A sufficiently large subgroup is essential for ensuring reliable treatment effect estimation and practical applicability. For instance, maintaining adequate statistical power is crucial for robust analysis, and in drug repurposing, having enough responders makes a treatment economically viable.
- •
Each sample in the identified subgroup adheres to the overlap assumption. Similar to excluding samples with extreme propensity scores (Crump et al. 2009), we impose that *each sample* in the selected subgroup has a propensity score bounded away from 0 and 1.

The above task can be formally stated as follows:

$$ $\displaystyle\underset{\tilde{S}}{\min}\;$ $\displaystyle-\frac{\sum_{i=1}^{n}\mathds{1}(\tilde{S}(\boldsymbol{X})=1)\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})}{\sum_{i=1}^{n}\mathds{1}(\tilde{S}(\boldsymbol{X})=1)}$ ( Problem I ) s.t. $\displaystyle\quad\mathds{E}\left[\tilde{S}(\boldsymbol{X})\right]\geq c$ $\displaystyle\quad\alpha\leq e(\boldsymbol{x})\leq 1-\alpha,\forall\boldsymbol{x}:\tilde{S}(\boldsymbol{x})=1,$ (2) $$

where $c\in(0,1)$ is the desired size of the subgroup, and $\alpha\in[0,0.5)$ is the threshold controlling the overlap constraint. While we focus on the subgroup size and overlap constraints here, our method naturally accommodates more general linear and ratio-form constraints, which we formally introduce in Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2) and Remark [2](https://arxiv.org/html/2504.20908v2#Thmremark2) (Section [4.2](https://arxiv.org/html/2504.20908v2#S4.SS2)).

### 3.2 Relaxing the Combinatorial Formulation

Since the [Problem I](https://arxiv.org/html/2504.20908v2#S3.Ex3) is
combinatorial and difficult to optimize, we relax the subgroup identification to a probabilistic assignment. This relaxation is implemented using a parametric surrogate model $S:\mathds{R}^{d}\mapsto(0,1)$ with parameters $\boldsymbol{\theta}$.
The ATE on the identified subgroup can be expressed as

$$ $f(\boldsymbol{\theta})=\frac{\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\boldsymbol{\theta})\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})}{\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\boldsymbol{\theta})}.$ (3) $$

Let $\hat{e}(\boldsymbol{x}_{i})$, $\hat{\mu}_{a}(\boldsymbol{x}_{i})$ denote the estimated $e(\boldsymbol{x}_{i})$ and $\mu_{a}(\boldsymbol{x}_{i})$. We then adopt the AIPTW estimator,
which can be expressed as functions of
$\hat{e}(\boldsymbol{x}_{i})$ and $\hat{\mu}_{a}(\boldsymbol{x}_{i})$:

$$ $\displaystyle\hat{\phi}_{\mathrm{aiptw}}(\boldsymbol{x}_{i},a_{i},y_{i})=$ $\displaystyle\hat{\mu}_{1}(\boldsymbol{x}_{i})-\hat{\mu}_{0}(\boldsymbol{x}_{i})+\frac{a_{i}}{\hat{e}(\boldsymbol{x}_{i})}(y_{i}-\hat{\mu}_{1}(\boldsymbol{x}_{i}))$ $\displaystyle-\frac{1-a_{i}}{1-\hat{e}(\boldsymbol{x}_{i})}(y_{i}-\hat{\mu}_{0}(\boldsymbol{x}_{i})).$ (4) $$

While this study focus on AIPTW, $\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})$ can be derived using other ATE estimators, making this formulation flexible. More details are illustrated in Appendix [C.1](https://arxiv.org/html/2504.20908v2#A3.SS1).

Due to the relaxation of subgroup identification into a probabilistic assignment, Constraint ([2](https://arxiv.org/html/2504.20908v2#S3.E2)), originally designed for discrete subgroup selection, must be adapted.
To achieve this, we introduce $h(\boldsymbol{x}_{i},\alpha)$, a surrogate function that reformulates the overlap constraint from [Problem I](https://arxiv.org/html/2504.20908v2#S3.Ex3):

$$ $h(\boldsymbol{x}_{i},\alpha)=1-\frac{\hat{e}(\boldsymbol{x}_{i})(1-\hat{e}(\boldsymbol{x}_{i}))}{\alpha(1-\alpha)}.$ (5) $$

The following result, Lemma [1](https://arxiv.org/html/2504.20908v2#Thmlemma1) (proof in Appendix [B.1](https://arxiv.org/html/2504.20908v2#A2.SS1)), establishes that Constraint ([2](https://arxiv.org/html/2504.20908v2#S3.E2)) in [Problem I](https://arxiv.org/html/2504.20908v2#S3.Ex3) can be replaced by a constraint on $h(\boldsymbol{x}_{i},\alpha)$:

###### Lemma 1 .

$S(\boldsymbol{x}_{i};\boldsymbol{\theta})h(\boldsymbol{x}_{i},\alpha)\leq 0$ if and only if $\alpha\leq\hat{e}(\boldsymbol{x}_{i})\leq 1-\alpha$.

With Lemma [1](https://arxiv.org/html/2504.20908v2#Thmlemma1), our optimization can be reformulated as:

$$ $\displaystyle\underset{\boldsymbol{\theta}}{\min}\quad$ $\displaystyle-f(\boldsymbol{\theta})$ ( Problem II ) s.t. $\displaystyle\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\boldsymbol{\theta})\geq c$ (6) $\displaystyle S(\boldsymbol{x}_{i};\boldsymbol{\theta})h(\boldsymbol{x}_{i},\alpha)\leq 0,\forall i.$ (7) $$

Solving [Problem II](https://arxiv.org/html/2504.20908v2#S3.Ex6) remains a significant challenge.
As shown in Equation ([3](https://arxiv.org/html/2504.20908v2#S3.E3)), the parametric
model $S(\boldsymbol{x_{i}};\boldsymbol{\theta})$ appears in both the numerator and denominator, making the objective function $f(\boldsymbol{\theta})$ neither convex nor concave, even for simple models like logistic regression. This nonconvexity also extends to the feasible set, rendering standard convex optimization methods, such as ADMM and FW unsuitable. Although it could, in principle, be solved using Lagrangian relaxation, doing so would require tuning a separate multiplier for each constraint, making the hyperparameter tuning process impractical at scale.
To address this, we reformulate [Problem II](https://arxiv.org/html/2504.20908v2#S3.Ex6) and present the final framework in Section [4](https://arxiv.org/html/2504.20908v2#S4).

###### Lemma 1 .

## 4 Optimization Methods

Section [4.1](https://arxiv.org/html/2504.20908v2#S4.SS1) reformulates the task
as a min-max optimization
and adopts the *$\gamma$-Gradient Descent Ascent* ($\gamma$-GDA) algorithm (Schweisthal et al. 2024).
To improve numerical stability, we modified the formulation
and establish feasibility guarantees in Section [4.2](https://arxiv.org/html/2504.20908v2#S4.SS2), showing that MOSIC can identify the optimal subgroup while satisfying multiple constraints.

### 4.1 Min-Max Formulation and GDA

Since neither the objective nor the feasible region is convex,
we first rewrite [Problem II](https://arxiv.org/html/2504.20908v2#S3.Ex6) using the saddle point formulation in [Problem III](https://arxiv.org/html/2504.20908v2#S4.Ex7), converting it to a min-max problem. We then solve it by applying $\gamma$-GDA, as described in Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1).
Our goal is to minimize
the resulting
*Lagrangian*, $L(\boldsymbol{\theta},\boldsymbol{\lambda})$, associated with [Problem II](https://arxiv.org/html/2504.20908v2#S3.Ex6):

$$ $\displaystyle L(\boldsymbol{\theta},$ $\displaystyle\boldsymbol{\lambda}):=-f(\boldsymbol{\theta})+\boldsymbol{\lambda}^{T}\boldsymbol{g}(\boldsymbol{\theta}),$ (8) $$

where

$$ $\boldsymbol{\lambda}=\begin{bmatrix}\lambda_{1}\\ \vdots\\ \lambda_{n}\\ \lambda_{n+1}\end{bmatrix},\quad\boldsymbol{g}(\boldsymbol{\theta})=\begin{bmatrix}S(x_{1};\boldsymbol{\theta})h(x_{1};\alpha)\\ \vdots\\ S(x_{n};\boldsymbol{\theta})h(x_{n};\alpha)\\ c-\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\end{bmatrix}.$ (9) $$

For notational convenience, we write $\boldsymbol{g}(\boldsymbol{\theta})$ without explicitly indicating its dependence on fixed constraint values $c$ and $\alpha$ (or those for any additional constraint, if present).
We solve the primal constrained problem directly through a *min-max Lagrangian formulation*, unlike Lagrangian relaxation which operates on the dual problem. This *saddle-point approach* preserves feasibility guarantees while avoiding the instability of dual methods under non-convex constraints. The solution to [Problem II](https://arxiv.org/html/2504.20908v2#S3.Ex6) is equivalent to solving (Boyd and Vandenberghe 2004):

$$ $\underset{\boldsymbol{\theta}}{\min}\;\underset{\boldsymbol{\lambda}\geq 0}{\max}\;L(\boldsymbol{\theta},\boldsymbol{\lambda}).$ ( Problem III ) $$

We adopt $\gamma$-GDA algorithm (Jin, Netrapalli, and Jordan 2020) to solve this problem.
As shown in Alrorithm [1](https://arxiv.org/html/2504.20908v2#alg1), each iteration performs a gradient descent step on $\boldsymbol{\theta}$ followed by a gradient ascent step on $\boldsymbol{\lambda}$. The descent step’s learning rate decays over time, as required by Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1), which underpins Proposition [1](https://arxiv.org/html/2504.20908v2#Thmproposition1) and Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2).

Figure: Algorithm 1 $\gamma$-Gradient Descent Ascent ($\gamma$-GDA)

While appearing effective, Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1) exhibits notable numerical instability. To illustrate, consider a scenario with a single group size constraint (Constraint ([6](https://arxiv.org/html/2504.20908v2#S3.E6))) and its Lagrange multiplier $\lambda_{n+1}$.
Initially, the constraint is violated because the group size is below the threshold,
leading to a positive gradient
that increments $\lambda_{n+1}$ at each iteration.
Once the constraint is met, the gradient becomes negative, driving $\lambda_{n+1}$ toward 0.
However, due to the small learning rate, $\lambda_{n+1}$ remains nonzero in subsequent iterations, continuing to penalize the objective and lead to artificially shrunk feasible regions.
This instability extends to multiple constraints, causing the algorithm to oscillate between a strictly feasible but suboptimal $\boldsymbol{\theta}$ and an optimal but infeasible $\boldsymbol{\theta}$, undermining convergence. To address this, we introduce a modified objective in Section [4.2](https://arxiv.org/html/2504.20908v2#S4.SS2) and provide theoretical guarantees.

### 4.2 Modified Min-Max problem and Feasibility Guarantees

We make two additional modifications to [Problem III](https://arxiv.org/html/2504.20908v2#S4.Ex7):

- a)
Deactivating the corresponding penalty term once a constraint is satisfied, and
- b)
Incorporating additional penalty terms to ensure convergence to a feasible, locally optimal solution.

For [a)](https://arxiv.org/html/2504.20908v2#S4.I1.i1), we wrap the constraint terms with ReLU function:

$$ $L^{\prime}(\boldsymbol{\theta},\boldsymbol{\lambda})=-f(\boldsymbol{\theta})+\boldsymbol{\lambda}^{T}ReLU(\boldsymbol{g}(\boldsymbol{\theta})).$ (10) $$

When $g_{i}(\boldsymbol{\theta})\leq 0$ (i.e., the constraint is satisfied), the penalty term vanishes (i.e., $ReLU(g_{i}(\boldsymbol{\theta}))=0$), preventing further penalization of the objective. This approach excludes satisfied constraints from optimization, improving stability.

For [b)](https://arxiv.org/html/2504.20908v2#S4.I1.i2), we begin by defining local optimality in the context of a nonconvex-nonconcave min-max problem, introducing the concept of *local minmax point* (Definition [1](https://arxiv.org/html/2504.20908v2#Thmdefinition1)). We then modify the objective
to ensure that
Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1) converges to a local minmax point (Proposition [1](https://arxiv.org/html/2504.20908v2#Thmproposition1)).
Finally,
we establish the
feasibility of the returned solution in
Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2).

A local minmax point (Definition [1](https://arxiv.org/html/2504.20908v2#Thmdefinition1)), informally, is a fixed point where the objective remains stable under small perturbations in $\boldsymbol{\theta}$ (the parameters over which we minimize) and worst-case perturbations in $\boldsymbol{\lambda}$ (the parameters over which we maximize).
Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1) (Jin, Netrapalli, and Jordan 2020) states that if the min-max objective function is twice differentiable, then the $\gamma$-GDA algorithm, upon convergence, reaches either a local minmax point or a stationary point with a degenerate (i.e., not invertible) Hessian.

However, the Hessian of our objective in Equation ([10](https://arxiv.org/html/2504.20908v2#S4.E10)) is inherently degenerate:
$\nabla^{2}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}L(\boldsymbol{\theta},\boldsymbol{\lambda})=\boldsymbol{0}$.
This degeneracy prevents us from proving that the algorithm converges to local minmax points,
and further, hinders our ability to guarantee the feasibility of the resulting solution.
Unlike Nandwani et al. (2019), who did not exclude degenerate points–potentially invalidating their results–we introduce a
modification, as shown below, to
eliminates degenerate points, leading to the final optimization objective of MOSIC:

$$ $\tilde{L}(\boldsymbol{\theta},\boldsymbol{\lambda})=-f(\boldsymbol{\theta})+\boldsymbol{\lambda}^{T}ReLU(\boldsymbol{g}(\boldsymbol{\theta}))-\frac{\beta}{2}\boldsymbol{\lambda}^{2}.$ (11) $$

Now our problem reduces to solving [Problem IV](https://arxiv.org/html/2504.20908v2#S4.Ex8):

$$ $\underset{\boldsymbol{\theta}}{\min}\;\underset{\boldsymbol{\lambda}\geq 0}{\max}\;\tilde{L}(\boldsymbol{\theta}_{,}\boldsymbol{\lambda}).$ ( Problem IV ) $$

With this update, we establish that, if Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1) converges, it converges to a strict local minmax point:

###### Proposition 1 (Local Optimality) .

Let $(\boldsymbol{\theta}^{\prime},\boldsymbol{\lambda}^{\prime})$ be a *linearly stable point* (formally defined in Definition [2](https://arxiv.org/html/2504.20908v2#Thmdefinition2)) of Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1). Then, $(\boldsymbol{\theta}^{\prime},\boldsymbol{\lambda}^{\prime})$ must be a *strict local minmax point*.

The proof of Proposition [1](https://arxiv.org/html/2504.20908v2#Thmproposition1) (Appendix [B.2](https://arxiv.org/html/2504.20908v2#A2.SS2)) relaxes the twice-differentiability requirement in Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1), and shows that the algorithm cannot converge to degenerate points. Building on this,
the following result, Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2), shows that upon convergence, all constraints are approximately satisfied within a small tolerance. This guarantee applies not only to the group size constraint but also to arbitrary linear constraints, which include the overlap constraint.

###### Lemma 2 .

Let $(\boldsymbol{\theta}^{*},\boldsymbol{\lambda}^{*})$ be a strict local min-max point obtained by Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1). Suppose the constraint set includes:

- •
a group size constraint:
$g^{\text{size}}(\boldsymbol{\theta})=c-\frac{1}{n}\sum_{i=1}^{n}S(x_{i};\boldsymbol{\theta}),$
- •
and $K$ ($K\geq 0$) additional linear constraints of the form:
$g^{k}(\boldsymbol{\theta})=a^{k}+\sum_{i=1}^{n}b^{k}_{i}S(x_{i};\boldsymbol{\theta}),$

which together define the constraint vector:

$$ $\boldsymbol{g}(\boldsymbol{\theta}):=\begin{bmatrix}a^{1}+\sum_{i=1}^{n}b^{1}_{i}S(x_{i};\boldsymbol{\theta})\\ \vdots\\ a^{K}+\sum_{i=1}^{n}b^{K}_{i}S(x_{i};\boldsymbol{\theta})\\ c-\frac{1}{n}\sum_{i=1}^{n}S(x_{i};\boldsymbol{\theta})\end{bmatrix},$ $$

where $a^{k},b^{k}_{i}\in\mathbb{R}$, and, $\forall k,|\sum_{i=1}^{n}b^{k}_{i}|>0$.

Suppose $\beta$ satisfies:

$$ $\beta<\frac{\xi(c-\xi)|\mu_{\Delta}|}{2\,\phi_{\max}L},$ $$

where $\xi>0$ is a small tolerance error; $\phi_{\max}=\max_{i}\hat{\phi}(x_{i},a_{i},y_{i})$; $L=\sup|\frac{\partial S(\cdot;\boldsymbol{\theta})}{\partial\theta_{j}}|$, the coordinate-wise Lipschitz constant of $S(\cdot;\boldsymbol{\theta})$, $\mu_{\Delta}=\mathbb{E}\frac{\partial S(x_{i};\boldsymbol{\theta})}{\partial\theta_{j}}$, and $j=\underset{j}{\arg\max}\frac{|\mu_{\Delta}|}{L}$.

Then one of the following holds:

- 1.
Model collapse: $\frac{1}{n}\sum_{i=1}^{n}S(x_{i};\boldsymbol{\theta}^{*})<\xi.$
- 2.
Each constraint is approximately satisfied. For the group size constraint,
$g^{\text{size}}(\boldsymbol{\theta}^{*})\leq\xi.$
For the $k$-th linear constraint, with probability at least $1-\delta,$
$g^{k}(\boldsymbol{\theta}^{*})\leq\frac{\xi}{|\sum_{i=1}^{n}b_{i}^{k}|(1+\frac{L}{|\mu_{\Delta}|\sqrt{n}}\sqrt{\log\frac{2}{\delta}})}.$

The proof of Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2) (Appendix [B.3](https://arxiv.org/html/2504.20908v2#A2.SS3)) analyzes each constraint at strict local minmax points and establishes that, with $\beta$ properly chosen,(^2^22 In practice, relaxing $\beta$ does not substantially increase constraint violation. We recommend setting $\beta$ between $10^{-5}$ and 0.01.)
the constraints are either fully satisfied or violated up to a small tolerance error, if the model does not collapse. The proof further shows that, as long as the feasible region is non-negligible, the model is unlikely to collapse. In practice, if the model converges to a near-zero subgroup size, the run is restarted using a different random seed. Persistent collapse suggests that the feasible region defined by the multiple constraints should be carefully re-evaluated.

###### Remark 1 (Implication on the overlap constraint) .

When $n\rightarrow\infty$,
the bound on constraint violation is governed by $|\sum_{i=1}^{n}b_{i}^{k}|$. For the overlap constraint on sample $j$, we set $b^{k}_{i}=\mathds{1}(i=j)h(x_{j};\alpha)$, which gives $\sum_{i=1}^{n}b_{i}^{k}=h(\boldsymbol{x_{j}};\alpha)$. When $h(\boldsymbol{x_{j}};\alpha)$ is small, this denominator inflates the bound, suggesting potentially high violation. However, a small $h(\boldsymbol{x_{j}};\alpha)$ indicates that the corresponding violation of the overlap condition is itself negligible. Thus, such constraints can be safely ignored in practice.

###### Remark 2 (Extension to ratio constraints) .

Because the model maximizes the subgroup ATE, it naturally favors smaller subgroups by excluding samples with relatively low estimated CATE, while the group size constraint enforces a group size lower bound $c$. As a result, the group size typically converges to the size threshold $c$, i.e., $\sum_{i=1}^{n}S(x_{i};\boldsymbol{\theta})\approx c$ at termination. This observation allows us to extend the linear constraint to ratio constraints

$$ $g^{k}(\boldsymbol{\theta})=a^{k}+\frac{\sum_{i=1}^{n}b_{i}^{k}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}{\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}.$ $$

To retain compatibility with Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2), we block gradient flow through the denominator during backpropagation, effectively treating it as a constant.

The linear and ratio families capture many practical constraints relevant to healthcare applications. Linear constraints include the overlap constraint and budget constraints (each patient incurs a treatment cost under limited resources). Ratio constraints cover safety constraints (e.g., risk levels associated with patients) and fairness constraints. Section [5.4](https://arxiv.org/html/2504.20908v2#S5.SS4.SSSx2) assesses MOSIC’s extendability to these constraints. Due to the nonconvexity of both the objective and the feasible region, extending the guarantee to more complex constraints remains an open direction for future work.

Finally, we summarize our overall framework, MOSIC, in Algorithm [2](https://arxiv.org/html/2504.20908v2#alg2). It first estimates the nuisances and computes the pseudo-outcomes for each sample, which are then fed into the objective function ([Problem IV](https://arxiv.org/html/2504.20908v2#S4.Ex8)). The optimization is subsequently performed using the $\gamma$-GDA algorithm.

Figure: Algorithm 2 MOSIC

###### Proposition 1 (Local Optimality) .

###### Lemma 2 .

###### Remark 1 (Implication on the overlap constraint) .

###### Remark 2 (Extension to ratio constraints) .

## 5 Experiments

We evaluate MOSIC on both synthetic and real-world data. Section [5.1](https://arxiv.org/html/2504.20908v2#S5.SS1) describes the datasets, and Section [5.2](https://arxiv.org/html/2504.20908v2#S5.SS2) outlines the experimental setup.
Results in Section [5.3](https://arxiv.org/html/2504.20908v2#S5.SS3) demonstrate that MOSIC achieves high ATE while maintaining comparable covariate balance, or vice versa.
Finally, Section [5.4](https://arxiv.org/html/2504.20908v2#S5.SS4) presents ablation studies, highlighting that MOSIC outperforms baselines under the same interpretability requirement and can be extended to multiple additional constraints.

### 5.1 Datasets

##### Synthetic Data

We generate synthetic data following a procedure adapted from Assaad et al. (2021) (details in Appendix [C.2](https://arxiv.org/html/2504.20908v2#A3.SS2)).
We introduce an *imbalance parameter* $\tilde{\omega}\geq 0$ to determine the strength of confounding bias.
In particular, we generated two datasets with (1) no confounding bias ($\tilde{\omega}=0$); and (2) high confounding bias ($\tilde{\omega}=5$).

##### Real-World Data

We use two datasets of de-identified health records from intensive care units (ICU): eICU (Pollard et al. 2018) and MIMIC-IV (Johnson et al. 2023). Our study focuses on evaluating the effectiveness of Glucocorticoids in treating sepsis patients in the ICU. The treated group ($A=1$) consists of patients who received an initial Glucocorticoids dose of 160mg
within 10 hours before to 24 hours after ICU admission. The outcome $Y$ represents 7-day survival, with $Y=1$ indicating survival and $Y=0$ otherwise. The input $\boldsymbol{X}$ contains baseline covariates, including lab test results, vital signs, and sequential organ failure assessment (SOFA) scores (Vincent et al. 1996).

**Table 1: Summary of data sets. n: sample size; d: data dimension; Y: outcome.**
| Dataset | $n$ | $d$ | $\mathrm{Y}$ |
| --- | --- | --- | --- |
| Synthetic | 5,000 | 10 | Continuous |
| eICU | 13,361 | 23 | Binary, Survival |
| MIMIC-IV | 6,516 | 23 | Binary, Survival |

### 5.2 Setups

#### Baselines

We compare MOSIC with two categories from Section [2](https://arxiv.org/html/2504.20908v2#S2): Those designed for subgroup identification and those adapted from CATE estimation algorithms.

- •
Dedicated subgroup identification methods: We evaluate (1) CAPITAL, (2) OWL, and (3) VT, modifying them to incorporate a group size constraint via thresholding
(Details in Appendix [C.3](https://arxiv.org/html/2504.20908v2#A3.SS3)).
These methods prioritize interpretability, making it unclear whether performance limits stem from this trade-off or poor CATE estimation. To address this, we consider the next category of benchmarks.
- •
CATE estimation algorithms: We evaluate three methods: (1) CT, (2) CF, (3) Dragonnet.
Patients are ranked by the estimated CATE values and the top subgroup of the desired size is selected (VanderWeele et al. 2019).
Although these methods are not designed for subgroup identification, they serve as a natural baseline. If CATE estimation were accurate, this approach would identify the optimal subgroup, allowing us to separate the impact of estimation reliability from the interpretability trade-off.

#### Evaluation Metrics

We assess performance using two metrics: 1) Subgroup ATE, which measures whether identified subgroups achieve high ATE at the desired size (using ground‑truth ATE for synthetic data and the difference between subgroup and overall AIPTW estimates for real‑world data). 2) ATE Estimation Reliability measures the AIPTW estimation error on synthetic data and the number of unbalanced features on real-world data because the ground truth ATE is not observed. Although MOSIC can enforce the overlap constraint on the training set by design (Figure LABEL:fig:syn-error), directly evaluating constraint violations on the test set is less meaningful under distribution shifts. Instead, we assess ATE estimation reliability, which captures the intended effect of the overlap constraint.

A feature is considered unbalanced (Cohen 2013) if its standardized mean difference (SMD) $>0.2$ after IPTW-reweighting (Austin 2009; Zang et al. 2023). A higher number of unbalanced features indicates larger estimation uncertainty and lower reliability of the subgroup.
To validate feature imbalance as a proxy for estimation error, we further report the feature imbalance in synthetic data (Figure LABEL:fig:syn-error).

#### Implementation Details

For all datasets, we perform 100 random splits of the training and test sets. For each split, we first conduct a 5-fold cross-validation on the training set to determine the optimal hyperparameters (see Appendix [C.5](https://arxiv.org/html/2504.20908v2#A3.SS5) for details). The model is then retrained on the entire training set using the selected hyperparameters and evaluated on the corresponding test set. The final results are reported using the mean and standard deviations across all 100 evaluations.

We implement the subgroup identification model ($S$) using MLP and DT. In the next sections, ‘MOSIC’ refers to MOSIC-MLP, unless stated otherwise. Following (Marton et al. 2024), we represent the DT as a neural network. These models are trained using Algorithm [1](https://arxiv.org/html/2504.20908v2#alg1), with L1 regularization applied in the loss function. For nuisance function estimation, we use LR for the propensity score model $\hat{e}(\cdot)$, and Dragonnet for the outcome models $\hat{\mu}_{a}(\cdot)$ (Appendix [C.4](https://arxiv.org/html/2504.20908v2#A3.SS4)).
(^3^33Empirically, estimating propensity scores with Dragonnet results in a large number of unbalanced features; we therefore adopt LR, demonstrating the flexibility of MOSIC.)

### 5.3 Results

#### Synthetic Data

To assess the impact of overlap constraints, we compare MOSIC with ($\alpha=0.02$) and without ($\alpha=0$) them.
Figure [1](https://arxiv.org/html/2504.20908v2#S5.F1) and Figure [3](https://arxiv.org/html/2504.20908v2#S5.F3) present results from the confounded and unconfounded synthetic data.
Compared with baselines, the left column indicates that MOSIC consistently identifies subgroups with the highest ATE across all group sizes, and the right column shows that it achieves the lowest ATE estimation errors.

Further, we numerically verify that MOSIC can indeed satisfy the overlap constraint and investigate the correlation between the estimation error and the number of unbalanced features (Figure LABEL:fig:syn-error). We also investigate the statistical properties of the proposed procedure in Appendix [F](https://arxiv.org/html/2504.20908v2#A6) and  [G](https://arxiv.org/html/2504.20908v2#A7).

Figure: Figure 1: True ATE and ATE estimation error across different group sizes on confounded synthetic data($\omega=5$).
Refer to caption: Figures/syn-gamma-5.png

#### Real-World Data

Figures LABEL:fig:real-world demonstrates that MOSIC consistently outperforms other methods. We highlight the importance of jointly evaluating subgroup ATE and covariate balance when assessing performance. A large ATE alone is insufficient if covariate imbalance persists, as it may indicate unreliable estimates. As shown in Figures LABEL:fig:real-world (left), MOSIC achieves higher subgroup ATEs at comparable group sizes in most cases. In scenarios where the ATE advantage is not statistically significant (e.g., MOSIC ($\alpha=0.01$) vs. CF at $c=0.6$ on MIMIC, p = 0.68), MOSIC still exhibits improved covariate balance (e.g., MOSIC ($\alpha=0.01$) vs. CF at $c=0.6$ on MIMIC, p = 0.000052). Full statistical test results are provided in Appendix [D](https://arxiv.org/html/2504.20908v2#A4). Additionally, the overlap constraint in MOSIC significantly improves covariate balance on eICU (Figures LABEL:subfig:real-world-eicu, right).

### 5.4 Ablation Studies

#### Decision Tree as Base Model

While MOSIC-MLP achieves stronger performance than the baselines, its black-box nature limits interpretability. In contrast, decision trees are favored in clinical settings for their transparency, allowing stakeholders to understand and trust subgroup recommendations (Cai et al. 2022).
We compare MOSIC with DT implementation to other DT-based baselines (CT, CF) on eICU. VT is excluded due to its high ATE estimation error on synthetic data (Figure [1](https://arxiv.org/html/2504.20908v2#S5.F1)). To impose equal interpretability requirement, we fix the tree depth to 5 for both MOSIC (denoted as MOSIC-DT) and CT, and use ensembles of 3 trees of depth 5 for MOSIC (denoted as MOSIC-Forest) and CF. As shown in Figure [E.2.1](https://arxiv.org/html/2504.20908v2#A5.SS2.F1), MOSIC consistently outperforms CT and CF under the same interpretability requirement. Notably, despite the limited model capacity, MOSIC effectively enforces both group size and overlap constraints, leading to improved covariate balance.

#### Extension to Additional Constraints

While this study focuses on the group size and overlap constraints, MOSIC can be readily extended to multiple constraints. To illustrate, we experimented on synthetic data to progressively enforce safety, budget, and fairness constraints (Appendix [E.4](https://arxiv.org/html/2504.20908v2#A5.SS4)). We also experimented on eICU to discourage the selection of patients with injured central nervous system (CNS), as studies have shown that glucocorticoids may exacerbate neural damage (Hill and Spencer-Segal 2021). Specifically, we required that the proportion of patients with a Glasgow Coma Scale (GCS) score $<6$ remain below 0.05. A GCS$<6$ corresponds to the most severe CNS dysfunction in the SOFA measurement, a widely used ICU metric. This constraint was applied in addition to the group size constraint ($c=0.4$) and overlap constraint ($\alpha=0.01$). Table [2](https://arxiv.org/html/2504.20908v2#S5.T2) and Appendix [E.5](https://arxiv.org/html/2504.20908v2#A5.SS5) demonstrate that MOSIC can additionally satisfy this safety constraint. The findings demonstrate that MOSIC can effectively accommodate multiple constraints simultaneously.

**Table 2: eICU with GCS score constraint**
| Metric | before | after |
| --- | --- | --- |
| Group Size | $0.46\pm 0.10$ | $0.44\pm 0.10$ |
| # Unbalanced Features | $2.0\pm 1.72$ | $2.2\pm 2.28$ |
| ATE Improvement | $0.02\pm 0.02$ | $0.02\pm 0.02$ |
| Proportion of GCS$<6$ | $0.15\pm 0.05$ | $0.03\pm 0.03$ |

## 6 Conclusion

We propose MOSIC, a model-agnostic framework for optimal subgroup identification that handles multiple constraints with feasibility guarantees. It demonstrates strong empirical performance under group size and overlap constraints, flexibly extends to additional clinical constraints, and supports diverse models to deliver interpretable solutions.

## References

- Assaad et al. (2021)
Assaad, S.; Zeng, S.; Tao, C.; Datta, S.; Mehta, N.; Henao, R.; Li, F.; and Carin, L. 2021.
Counterfactual representation learning with balancing weights.
In *International Conference on Artificial Intelligence and Statistics*, 1972–1980. PMLR.
- Athey and Imbens (2016)
Athey, S.; and Imbens, G. 2016.
Recursive partitioning for heterogeneous causal effects.
*Proceedings of the National Academy of Sciences*, 113(27): 7353–7360.
- Athey and Wager (2021)
Athey, S.; and Wager, S. 2021.
Policy learning with observational data.
*Econometrica*, 89(1): 133–161.
- Austin (2009)
Austin, P. C. 2009.
Using the standardized difference to compare the prevalence of a binary variable between two groups in observational research.
*Communications in statistics-simulation and computation*, 38(6): 1228–1234.
- Boyd et al. (2011)
Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; Eckstein, J.; et al. 2011.
Distributed optimization and statistical learning via the alternating direction method of multipliers.
*Foundations and Trends® in Machine learning*, 3(1): 1–122.
- Boyd and Vandenberghe (2004)
Boyd, S.; and Vandenberghe, L. 2004.
*Convex optimization*.
Cambridge university press.
- Cai et al. (2022)
Cai, H.; Lu, W.; Marceau West, R.; Mehrotra, D. V.; and Huang, L. 2022.
CAPITAL: Optimal subgroup identification via constrained policy tree search.
*Statistics in medicine*, 41(21): 4227–4244.
- Cai et al. (2011)
Cai, T.; Tian, L.; Wong, P. H.; and Wei, L.-J. 2011.
Analysis of randomized comparative clinical trial data for personalized treatment selections.
*Biostatistics*, 12(2): 270–282.
- Chernozhukov et al. (2018)
Chernozhukov, V.; Chetverikov, D.; Demirer, M.; Duflo, E.; Hansen, C.; Newey, W.; and Robins, J. 2018.
Double/debiased machine learning for treatment and structural parameters.
- Cohen (2013)
Cohen, J. 2013.
*Statistical power analysis for the behavioral sciences*.
routledge.
- Cole and Hernán (2008)
Cole, S. R.; and Hernán, M. A. 2008.
Constructing inverse probability weights for marginal structural models.
*American journal of epidemiology*, 168(6): 656–664.
- Crump et al. (2009)
Crump, R. K.; Hotz, V. J.; Imbens, G. W.; and Mitnik, O. A. 2009.
Dealing with limited overlap in estimation of average treatment effects.
*Biometrika*, 96(1): 187–199.
- Doubleday et al. (2022)
Doubleday, K.; Zhou, J.; Zhou, H.; and Fu, H. 2022.
Risk controlled decision trees and random forests for precision Medicine.
*Statistics in medicine*, 41(4): 719–735.
- Dusseldorp, Doove, and Mechelen (2016)
Dusseldorp, E.; Doove, L.; and Mechelen, I. v. 2016.
Quint: An R package for the identification of subgroups of clients who differ in which treatment alternative is best for them.
*Behavior research methods*, 48: 650–663.
- Foster, Taylor, and Ruberg (2011)
Foster, J. C.; Taylor, J. M.; and Ruberg, S. J. 2011.
Subgroup identification from randomized clinical trial data.
*Statistics in medicine*, 30(24): 2867–2880.
- Gao, Goldfarb, and Curtis (2020)
Gao, W.; Goldfarb, D.; and Curtis, F. E. 2020.
ADMM for multiaffine constrained optimization.
*Optimization Methods and Software*, 35(2): 257–303.
- Gruber et al. (2022)
Gruber, S.; Phillips, R. V.; Lee, H.; and van der Laan, M. J. 2022.
Data-Adaptive Selection of the Propensity Score Truncation Level for Inverse-Probability–Weighted and Targeted Maximum Likelihood Estimators of Marginal Point Treatment Effects.
*American Journal of Epidemiology*, 191(9): 1640–1651.
- Hill and Spencer-Segal (2021)
Hill, A. R.; and Spencer-Segal, J. L. 2021.
Glucocorticoids and the brain after critical illness.
*Endocrinology*, 162(3): bqaa242.
- Huang et al. (2017)
Huang, X.; Sun, Y.; Trow, P.; Chatterjee, S.; Chakravartty, A.; Tian, L.; and Devanarayan, V. 2017.
Patient subgroup identification for clinical drug development.
*Statistics in medicine*, 36(9): 1414–1428.
- Jin, Netrapalli, and Jordan (2020)
Jin, C.; Netrapalli, P.; and Jordan, M. 2020.
What is Local Optimality in Nonconvex-Nonconcave Minimax Optimization?
In III, H. D.; and Singh, A., eds., *Proceedings of the 37th International Conference on Machine Learning*, volume 119 of *Proceedings of Machine Learning Research*, 4880–4889. PMLR.
- Johnson et al. (2023)
Johnson, A. E.; Bulgarelli, L.; Shen, L.; Gayles, A.; Shammout, A.; Horng, S.; Pollard, T. J.; Hao, S.; Moody, B.; Gow, B.; et al. 2023.
MIMIC-IV, a freely accessible electronic health record dataset.
*Scientific data*, 10(1): 1.
- Kallus (2020)
Kallus, N. 2020.
More Efficient Policy Learning via Optimal Retargeting.
arXiv:1906.08611.
- Kiriakidou and Diou (2022)
Kiriakidou, N.; and Diou, C. 2022.
An evaluation framework for comparing causal inference models.
In *Proceedings of the 12th Hellenic Conference on Artificial Intelligence*, 1–9.
- Kosorok and Laber (2019)
Kosorok, M. R.; and Laber, E. B. 2019.
Precision medicine.
*Annual review of statistics and its application*, 6(1): 263–286.
- Kravitz, Duan, and Braslow (2004)
Kravitz, R. L.; Duan, N.; and Braslow, J. 2004.
Evidence-based medicine, heterogeneity of treatment effects, and the trouble with averages.
*The Milbank Quarterly*, 82(4): 661–687.
- Künzel et al. (2019)
Künzel, S. R.; Sekhon, J. S.; Bickel, P. J.; and Yu, B. 2019.
Metalearners for estimating heterogeneous treatment effects using machine learning.
*Proceedings of the national academy of sciences*, 116(10): 4156–4165.
- Lacoste-Julien (2016)
Lacoste-Julien, S. 2016.
Convergence rate of frank-wolfe for non-convex objectives.
*arXiv preprint arXiv:1607.00345*.
- Lipkovich et al. (2011)
Lipkovich, I.; Dmitrienko, A.; Denne, J.; and Enas, G. 2011.
Subgroup identification based on differential effect search—a recursive partitioning method for establishing response to treatment in patient subpopulations.
*Statistics in medicine*, 30(21): 2601–2621.
- Lipkovich et al. (2024)
Lipkovich, I.; Svensson, D.; Ratitch, B.; and Dmitrienko, A. 2024.
Modern approaches for evaluating treatment effect heterogeneity from clinical trials and observational data.
*Statistics in Medicine*, 43(22): 4388–4436.
- Liu et al. (2018)
Liu, Y.; Wang, Y.; Kosorok, M. R.; Zhao, Y.; and Zeng, D. 2018.
Augmented outcome-weighted learning for estimating optimal dynamic treatment regimens.
*Statistics in medicine*, 37(26): 3776–3788.
- Marton et al. (2024)
Marton, S.; Lüdtke, S.; Bartelt, C.; and Stuckenschmidt, H. 2024.
GradTree: Learning axis-aligned decision trees with gradient descent.
In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 38, 14323–14331.
- Mehrabi et al. (2021)
Mehrabi, N.; Morstatter, F.; Saxena, N.; Lerman, K.; and Galstyan, A. 2021.
A survey on bias and fairness in machine learning.
*ACM computing surveys (CSUR)*, 54(6): 1–35.
- Nandwani et al. (2019)
Nandwani, Y.; Pathak, A.; Mausam; and Singla, P. 2019.
A primal-dual formulation for deep learning with constraints.
In *Proceedings of the 33rd International Conference on Neural Information Processing Systems*, 12179–12190.
- Pollard et al. (2018)
Pollard, T. J.; Johnson, A. E.; Raffa, J. D.; Celi, L. A.; Mark, R. G.; and Badawi, O. 2018.
The eICU Collaborative Research Database, a freely available multi-center database for critical care research.
*Scientific data*, 5(1): 1–13.
- Qiu, Carone, and Luedtke (2022)
Qiu, H.; Carone, M.; and Luedtke, A. 2022.
Individualized treatment rules under stochastic treatment cost constraints.
*Journal of causal inference*, 10(1): 480–493.
- Robins, Rotnitzky, and Zhao (1995)
Robins, J. M.; Rotnitzky, A.; and Zhao, L. P. 1995.
Analysis of semiparametric regression models for repeated outcomes in the presence of missing data.
*Journal of the american statistical association*, 90(429): 106–121.
- Schweisthal et al. (2024)
Schweisthal, J.; Frauen, D.; Melnychuk, V.; and Feuerriegel, S. 2024.
Reliable off-policy learning for dosage combinations.
In *Proceedings of the 37th International Conference on Neural Information Processing Systems*, NIPS ’23. Red Hook, NY, USA: Curran Associates Inc.
- Shalit, Johansson, and Sontag (2017)
Shalit, U.; Johansson, F. D.; and Sontag, D. 2017.
Estimating individual treatment effect: generalization bounds and algorithms.
In *International conference on machine learning*, 3076–3085. PMLR.
- Shi, Blei, and Veitch (2019)
Shi, C.; Blei, D.; and Veitch, V. 2019.
Adapting neural networks for the estimation of treatment effects.
*Advances in neural information processing systems*, 32.
- Tesei et al. (2023)
Tesei, G.; Giampanis, S.; Shi, J.; and Norgeot, B. 2023.
Learning end-to-end patient representations through self-supervised covariate balancing for causal treatment effect estimation.
*Journal of Biomedical Informatics*, 140: 104339.
- Van Der Laan and Rubin (2006)
Van Der Laan, M. J.; and Rubin, D. 2006.
Targeted maximum likelihood learning.
*The international journal of biostatistics*, 2(1).
- VanderWeele et al. (2019)
VanderWeele, T. J.; Luedtke, A. R.; van der Laan, M. J.; and Kessler, R. C. 2019.
Selecting optimal subgroups for treatment using many covariates.
*Epidemiology*, 30(3): 334–341.
- Vincent et al. (1996)
Vincent, J. L.; Moreno, R.; Takala, J.; Willatts, S.; De Mendonça, A.; Bruining, H.; Reinhart, C.; Suter, P.; and Thijs, L. G. 1996.
The SOFA (Sepsis-related Organ Failure Assessment) score to describe organ dysfunction/failure: On behalf of the Working Group on Sepsis-Related Problems of the European Society of Intensive Care Medicine (see contributors to the project in the appendix).
- Wager and Athey (2018)
Wager, S.; and Athey, S. 2018.
Estimation and inference of heterogeneous treatment effects using random forests.
*Journal of the American Statistical Association*, 113(523): 1228–1242.
- Wang and Rudin (2021)
Wang, T.; and Rudin, C. 2021.
Causal Rule Sets for Identifying Subgroups with Enhanced Treatment Effect.
arXiv:1710.05426.
- Wang, Yin, and Zeng (2019)
Wang, Y.; Yin, W.; and Zeng, J. 2019.
Global convergence of ADMM in nonconvex nonsmooth optimization.
*Journal of Scientific Computing*, 78: 29–63.
- Zang et al. (2023)
Zang, C.; Zhang, H.; Xu, J.; Zhang, H.; Fouladvand, S.; Havaldar, S.; Cheng, F.; Chen, K.; Chen, Y.; Glicksberg, B. S.; et al. 2023.
High-throughput target trial emulation for Alzheimer’s disease drug repurposing with real-world data.
*Nature communications*, 14(1): 8180.
- Zhao et al. (2012)
Zhao, Y.; Zeng, D.; Rush, A. J.; and Kosorok, M. R. 2012.
Estimating individualized treatment rules using outcome weighted learning.
*Journal of the American Statistical Association*, 107(499): 1106–1118.
- Zhou et al. (2024)
Zhou, J.; Yang, L.; Liu, X.; Gu, X.; Sun, L.; and Chen, W. 2024.
CURLS: Causal Rule Learning for Subgroups with Significant Treatment Effect.
In *Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, KDD ’24, 4619–4630. New York, NY, USA: Association for Computing Machinery.
ISBN 9798400704901.

## Appendix A Preliminary Definitions and Theorems

### A.1 Local minmax point

###### Definition 1 (Local minmax point) .

A point $(\mathbf{\theta}^{\star},\mathbf{\lambda}^{\star})$ is said to be a local minmax point of $L$, if there exists $\delta_{0}>0$ and a continuous function $h$ satisfying $h(\delta)\to 0$ as $\delta\to 0$, such that for any $\delta\leq\delta_{0}$, and any $(\mathbf{\theta},\mathbf{\lambda})$ satisfying

$$ $\|\mathbf{\theta}-\mathbf{\theta}^{\star}\|\leq\delta\quad\text{and}\quad\|\mathbf{\lambda}-\mathbf{\lambda}^{\star}\|\leq h(\delta),$ $$

we have

$$ $L(\mathbf{\theta}^{\star},\mathbf{\lambda})\leq L(\mathbf{\theta}^{\star},\mathbf{\lambda}^{\star})\leq\max_{\mathbf{\lambda}^{\prime}:\|\mathbf{\lambda}^{\prime}-\mathbf{\lambda}^{\star}\|\leq h(\delta)}L(\mathbf{\theta},\mathbf{\lambda}^{\prime}).$ $$

If a point $(\theta^{\star},\lambda^{\star})$ satisfy

$$ $[\nabla^{2}_{\theta\theta}L-\nabla^{2}_{\theta\lambda}L(\nabla^{2}_{\lambda\lambda}L)^{-1}\nabla^{2}_{\lambda\theta}L]\succ 0,$ $$

we call it a strict local minmax point (Jin, Netrapalli, and Jordan 2020).

###### Definition 1 (Local minmax point) .

### A.2 Linearly Stable Point

###### Definition 2 (Linearly Stable Point) .

For a differentiable dynamical system $\mathbf{w}$, a fixed point $\mathbf{z}^{\star}$ is a linearly stable point of $\mathbf{w}$ if its Jacobian matrix $\mathbf{J}(\mathbf{z}^{\star}):=\left(\frac{\partial\mathbf{w}}{\partial\mathbf{z}}\right)(\mathbf{z}^{\star})$ has spectral radius $\rho(\mathbf{J}(\mathbf{z}^{\star}))\leq 1$.

###### Definition 2 (Linearly Stable Point) .

### A.3 Convergence of γ \gamma italic_γ -GDA algorithm

###### Theorem 1 .

(Jin et al (Jin, Netrapalli, and Jordan 2020), Theorem 26): Given an objective: $\underset{x\ in\mathcal{X}}{\min}\underset{y\in\mathcal{Y}}{\max}f(x,y)$. For any twice-differentiable f, the strict linearly stable limit points of the $\gamma$-GDA flow are {strict local minmax points} $\cup$ {($\theta$,$\lambda$) — ($\theta$,$\lambda$) is stationary and $\nabla^{2}_{\lambda\lambda}f(\theta,\lambda)$ is degenerate} as $\gamma\rightarrow\infty$.

###### Theorem 1 .

## Appendix B Proofs

### B.1 Proof of Lemma 1

See [1](https://arxiv.org/html/2504.20908v2#Thmlemma1)

###### Proof.

We proceed by proving both directions separately.

Forward Direction: Suppose $S(x_{i};\theta)h(x_{i})\leq 0$. We aim to show that this implies $\alpha\leq e(x_{i})\leq 1-\alpha$.

$$ $\displaystyle S(x_{i};\theta)h(x_{i})\leq 0$ $\displaystyle\implies$ $\displaystyle h(x_{i})\leq 0\quad\quad\hfill{\text{(Since }S(x_{i})>0\text{)}}$ $\displaystyle\implies$ $\displaystyle 1-\frac{e(x_{i})(1-e(x_{i}))}{\alpha(1-\alpha)}\leq 0\quad\quad\hfill{\text{(Substituting }h(x_{i})\text{)}}$ $\displaystyle\implies$ $\displaystyle\alpha(1-\alpha)-e(x_{i})(1-e(x_{i}))\leq 0$ $\displaystyle\implies$ $\displaystyle e(x_{i})-\alpha-\left(e(x_{i})^{2}-\alpha^{2}\right)\geq 0$ $\displaystyle\implies$ $\displaystyle(e(x_{i})-\alpha)\left(1-(e(x_{i})+\alpha)\right)\geq 0$ $\displaystyle\implies$ $\displaystyle(e(x_{i})-\alpha)\left(1-\alpha-e(x_{i})\right)\geq 0$ $\displaystyle\implies$ $\displaystyle\alpha\leq e(x_{i})\leq 1-\alpha.$ $$

Backward Direction: Suppose $\alpha\leq e(x_{i})\leq 1-\alpha$. We need to show that this implies $S(x_{i};\theta)h(x_{i})\leq 0$.

Define the auxiliary function $q(z)=z(1-z)$, whose derivative is given by:

$$ $q^{\prime}(z)=1-2z.$ $$

Thus, $q(z)$ attains its maximum at $z=0.5$.

When $0<\alpha\leq e(x_{i})\leq 0.5$, we have $q^{\prime}(e(x_{i}))=1-2e(x_{i})\geq 0$, which implies that $q(e(x_{i}))$ is non-decreasing on $[0,0.5]$. Consequently, from the assumption $\alpha\leq e(x_{i})$, we obtain:

$$ $\displaystyle q(e(x_{i}))\geq q(\alpha)\implies e(x_{i})(1-e(x_{i}))\geq\alpha(1-\alpha).$ $$

- Similarly, for $0.5<e(x_{i})\leq 1-\alpha<1$, we have $e(x_{i})(1-e(x_{i}))\geq\alpha(1-\alpha)$.

By combining both cases, we conclude that:

$$ $e(x_{i})(1-e(x_{i}))\geq\alpha(1-\alpha).$ $$

Dividing both sides by $\alpha(1-\alpha)$ yields:

$$ $\displaystyle 1-\frac{e(x_{i})(1-e(x_{i}))}{\alpha(1-\alpha)}=h(x_{i})\leq 0.$ $$

Since $S(x_{i})>0$, it follows that:

$$ $S(x_{i})h(x_{i})\leq 0.$ $$

This completes the proof.
∎

###### Proof.

### B.2 Proof of Proposition 1

See [1](https://arxiv.org/html/2504.20908v2#Thmproposition1)

###### Proof.

The derivation of Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1) relies on the Jacobian matrix and requires twice differentiability. However, our objective function incorporates ReLU activations, introducing non-differentiability at the origin. This prevents us from directly applying Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1) in its standard form. Nonetheless, since the probability of encountering exact zero inputs to ReLU is negligible, our objective function remains effectively differentiable in practice when using gradient-based optimization.

Thus, the key arguments of Theorem [1](https://arxiv.org/html/2504.20908v2#Thmtheorem1) extend to our objective ([Problem IV](https://arxiv.org/html/2504.20908v2#S4.Ex8)), implying that its stable limit points must either be local minmax points or points where the second derivative is degenerate.

Further, we compute the first and second derivatives of $L$ with respect to $\lambda$:

$$ $\frac{\partial L}{\partial\lambda_{i}}=\operatorname{ReLU}(g_{i}(\theta))-\beta\lambda_{i},$ (12) $$

$$ $\frac{\partial^{2}L}{\partial\lambda_{i}\lambda_{j}}=\begin{cases}0&\text{if }i\neq j,\\ -\beta&\text{if }i=j.\end{cases}$ (13) $$

Since $\beta>0$, the Hessian matrix $\frac{\partial^{2}L}{\partial\lambda_{i}\lambda_{j}}$ is never degenerate. Applying Theorem [2](https://arxiv.org/html/2504.20908v2#Thmlemma2), we conclude that
the strict linearly stable limit points of the $\gamma$-GDA flow are precisely the set of local minmax points.
∎

###### Proof.

### B.3 Proof of Lemma 2

See [2](https://arxiv.org/html/2504.20908v2#Thmlemma2)

###### Proof.

At convergence, the following condition holds:

$$ $\frac{\partial L}{\partial\boldsymbol{\lambda}}\Big{|}_{\boldsymbol{\lambda}=\boldsymbol{\lambda^{*}}}=ReLU(\boldsymbol{g}(\boldsymbol{\theta^{*}}))-\beta\boldsymbol{\lambda^{*}}=\boldsymbol{0},$ (14) $$

$$ $\frac{\partial L}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}=-\frac{\partial f}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}+\boldsymbol{\lambda^{*}}\frac{\partial ReLU}{\partial\boldsymbol{g}}\cdot\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}=\boldsymbol{0}$ (15) $$

Since Eq. ([14](https://arxiv.org/html/2504.20908v2#A2.E14)) applised component-wise to $\boldsymbol{g}(\boldsymbol{\theta^{*}})$, we analyze each individual component $g(\boldsymbol{\theta^{*}})$ separately. Due to ReLU, we distinguish between two cases: $g(\boldsymbol{\theta^{*}})\leq 0$ and $g(\boldsymbol{\theta^{*}})>0$. When $g(\boldsymbol{\theta^{*}})\leq 0$, Lemma [2](https://arxiv.org/html/2504.20908v2#Thmlemma2) holds trivially. Therefore, we focus on the remaining case where $g(\boldsymbol{\theta^{*}})>0$, which means the constraint is violated and Eq. ([14](https://arxiv.org/html/2504.20908v2#A2.E14)) and Eq.([15](https://arxiv.org/html/2504.20908v2#A2.E15)) simplifies to

$$ $g(\boldsymbol{\theta^{*}})=\beta\lambda^{*},$ (16) $$

$$ $-\frac{\partial f}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}+\lambda^{*}\frac{\partial g}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}=\boldsymbol{0}$ (17) $$

Substituting Eq.([16](https://arxiv.org/html/2504.20908v2#A2.E16)) into Eq.([17](https://arxiv.org/html/2504.20908v2#A2.E17)), we obtain

$$ $-\frac{\partial f}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}+\frac{g(\boldsymbol{\theta^{*}})}{\beta}\frac{\partial g}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}=0.$ (18) $$

We note that
$g(\boldsymbol{\theta^{*}})$ is a scalar, meaning that
each elements of $\frac{\partial f}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}$ is proportional to the corresponding elements in $\frac{\partial g}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}$ by the same constant $\frac{g(\boldsymbol{\theta^{*}})}{\beta}>0$. This allows us to pick any element of $\boldsymbol{\theta}$ to analyze $g(\boldsymbol{\theta^{*}})$. (Note that when $\theta\in\mathbb{R}^{d}$ is high-dimensional, Eq.([18](https://arxiv.org/html/2504.20908v2#A2.E18)) is challenging to achieve, because we need one scalar to satisfy all $d$ equations. In practice, we observe small oscillations when the algorithm converges, which implies the algorithm finds it hard to exactly satisfy all equations.)

Denote $\phi_{i}=\phi(\boldsymbol{x_{i}},a_{i},y_{i})$, $\Delta_{i}=\frac{\partial S(\boldsymbol{x_{i}};\boldsymbol{\theta})}{\partial\theta_{j}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}\leq L$, $\mu_{\Delta}=\mathbb{E}\Delta_{i}$, where $j=\underset{j}{\arg\max}\frac{\mu_{\Delta}}{L}$. The definition of function $f$ (Eq.([3](https://arxiv.org/html/2504.20908v2#S3.E3))) results in

$$ $\displaystyle\Big{|}\frac{\partial f}{\partial\theta_{j}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}\Big{|}$ $\displaystyle=\Big{|}\frac{\sum_{i=1}^{n}\phi_{i}\Delta_{i}\sum_{s=1}^{n}S(\boldsymbol{x_{s}};\boldsymbol{\theta})}{\left(\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\right)^{2}}-$ $\displaystyle\frac{\sum_{i=1}^{n}\Delta_{i}\sum_{s=1}^{n}\phi_{s}S(\boldsymbol{x_{s}};\boldsymbol{\theta})}{\left(\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\right)^{2}}\Big{|}$ $\displaystyle=\Big{|}\frac{\sum_{i=1}^{n}\sum_{s=1}^{n}\Delta_{i}S(\boldsymbol{x_{s}};\boldsymbol{\theta})(\phi_{i}-\phi_{s})}{\left(\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\right)^{2}}\Big{|}$ $$

With $\phi_{max}=\max_{i}|\hat{\phi}(x_{i},a_{i},y_{i})|$, we have,

$$ $\displaystyle\left|\frac{\partial f}{\partial\theta_{j}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}\right|$ $\displaystyle\leq 2\phi_{max}\Big{|}\frac{\sum_{i=1}^{n}\sum_{s=1}^{n}\Delta_{i}S(\boldsymbol{x_{s}};\boldsymbol{\theta})}{\left(\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\right)^{2}}\Big{|}$ $\displaystyle=2\phi_{max}\Big{|}\frac{\sum_{i=1}^{n}\Delta_{i}\sum_{s=1}^{n}S(\boldsymbol{x_{s}};\boldsymbol{\theta})}{\left(\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})\right)^{2}}\Big{|}$ $\displaystyle=2\phi_{max}\frac{\Big{|}\sum_{i=1}^{n}\Delta_{i}\Big{|}}{\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ (19) $$

Substituting Eq. ([19](https://arxiv.org/html/2504.20908v2#A2.E19)) into Eq. ([18](https://arxiv.org/html/2504.20908v2#A2.E18)) and take assoluate value, we obtain

$$ $g(\boldsymbol{\theta^{*}})\leq\beta\left(\left|\frac{\partial g}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}\right|\right)^{-1}\frac{2\phi_{max}\left|\sum_{i=1}^{n}\Delta_{i}\right|}{\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ (20) $$

Next, we will prove the upper bound of the group size constraint violation and the general linear constraints violation by analyzing $\frac{\partial g}{\partial\boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}$.

Case 1: $g(\boldsymbol{\theta^{*}})=c-\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta^{*}})>0$: the group size constraint is violated.

$$ $\displaystyle\frac{\partial g}{\partial\theta_{j}}|_{\boldsymbol{\theta}=\boldsymbol{\theta^{*}}}=-\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}$ Substituting this into Eq. ( 20 ): $\displaystyle g(\boldsymbol{\theta^{*}})\leq\beta\left(\left|-\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\right|\right)^{-1}\frac{2\phi_{max}\left|\sum_{i=1}^{n}\Delta_{i}\right|}{\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ $\displaystyle=\beta\frac{2\phi_{max}}{\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ (21) $$

Substitute $\beta\leq\frac{\xi(c-\xi)\mu_{\Delta}}{2\phi_{max}L}$ into Eq. ([21](https://arxiv.org/html/2504.20908v2#A2.E21)), we obtain

$$ $g(\boldsymbol{\theta^{*}})\leq\frac{\xi(c-\xi)\mu_{\Delta}}{\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})L}\leq\frac{\xi(c-\xi)}{\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}=\frac{\xi(c-\xi)}{c-g(\boldsymbol{\theta^{*}})}$ $$

Solving for $g(\boldsymbol{\theta^{*}})$, we obtain $g(\boldsymbol{\theta^{*}})\leq\xi$ (group size constraint satisfied) or $g(\boldsymbol{\theta^{*}})\geq c-\xi$ (the model collapses, group size is near zero).

Next, we show that the model is improbable to collapse if the feasible region is non-negligible.

Define the collapsed region as

$$ $\boldsymbol{\Theta}_{\text{collapse}}:=\left\{\boldsymbol{\theta}\in\Theta:\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\boldsymbol{\theta})<\xi\right\},$ $$

and the feasible region as

$$ $\displaystyle\boldsymbol{\Theta}_{\text{feasible}}:=\Big{\{}\boldsymbol{\theta}\in\Theta:\;$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\boldsymbol{\theta})\geq c$ $\displaystyle\text{and other constraints are met}\Big{\}}$ $$

for some thresholds $0<\xi\ll c<1$. Suppose that the volume of the feasible region dominates that of the collapsed region, i.e.,

$$ $|\boldsymbol{\Theta}_{\text{collapse}}|\ll|\boldsymbol{\Theta}_{\text{feasible}}|.$ $$

Given $\forall i,0\leq S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})\leq 1$, using Hoeffding inequality, we obtain

$$ $P\left(\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})-\mathbb{E}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})\leq-t\right)\leq\exp(-\frac{2t^{2}}{n}).$ $$

Let $t=\sqrt{\frac{n\log(1/\delta)}{2}}$, with probability at least $1-\delta$, we have

$$ $\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})-\mathbb{E}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})\geq-\sqrt{\frac{n\log(1/\delta)}{2}}$ $$

Assume $\boldsymbol{\theta}$ is sampled uniformly from $\Theta_{\text{collapase}}\cup\Theta_{\text{feasible}}$, $|\Theta_{collapase}|\ll|\Theta_{feasible}|\Rightarrow P(\boldsymbol{\theta}\in\Theta_{\text{collapse}})\ll P(\boldsymbol{\theta}\in\Theta_{\text{feasible}})$. Then,

$$ $\displaystyle\mathbb{E}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})$ $\displaystyle=P(\boldsymbol{\theta}\in\Theta_{collapse})\cdot\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})$ $\displaystyle+P(\boldsymbol{\theta}\in\Theta_{feasible})\cdot\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})$ $$

Given that the contribution from $\Theta_{\text{collapse}}$ is negligible, we approximate:

$$ $\mathbb{E}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})\approx\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})|\{\boldsymbol{\theta}\in\Theta_{feasible}\}\geq nc$ $$

Therefore, with probability at least $1-\delta$, we have

$$ $\displaystyle\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})$ $\displaystyle\geq\mathbb{E}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})-\sqrt{\frac{n\log(1/\delta)}{2}}$ $\displaystyle\geq nc-\sqrt{\frac{n\log(1/\delta)}{2}}$ $\displaystyle\implies\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{\boldsymbol{x_{i}};\boldsymbol{\theta}})$ $\displaystyle\geq c-\sqrt{\frac{\log(1/\delta)}{2n}},$ $$

i.e. the model is improbable to collapse if the feasible region is non-negligible.

Case 2 : $g(\boldsymbol{\theta^{*}})=a^{k}+\boldsymbol{b}^{k\top}\boldsymbol{S}(\boldsymbol{x};\boldsymbol{\theta})>0$:

Similarly, we have

$$ $g(\boldsymbol{\theta^{*}})\leq\beta\left(\left|\sum_{i=1}^{n}b_{i}^{k}\Delta_{i}\right|\right)^{-1}\frac{2\phi_{max}\left|\sum_{i=1}^{n}\Delta_{i}\right|}{\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}.$ (22) $$

Let $Z=\sum_{i=1}^{n}b_{i}^{k}\Delta_{i}$, then $\mathbb{E}Z=\sum_{i=1}^{n}b_{i}^{k}\mathbb{E}\Delta_{i}=\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}.$

By the triangle inequality, we have $\Big{|}|Z|-|\mathbb{E}Z|\Big{|}\leq|Z-\mathbb{E}Z|$. So, by the Hoeffding inequality, we obtain

$$ $\displaystyle P(\Big{|}|Z|-|\mathbb{E}Z|\Big{|}\geq t)$ $\displaystyle\leq P(|Z-\mathbb{E}Z|\geq t)$ $\displaystyle\leq 2\exp(-\frac{2t^{2}}{4L^{2}\sum_{i=1}^{n}(b_{i}^{k})^{2}})$ $$

Let $t=L\sqrt{2\sum_{i=1}^{n}(b_{i}^{k})^{2}\log\frac{2}{\delta}}$. Then with probability at least $1-\delta$,

$$ $|Z|\geq|\mathbb{E}Z|+t=|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t.$ $$

Substituting this into Eq. ([22](https://arxiv.org/html/2504.20908v2#A2.E22)), we obtain

$$ $\displaystyle g(\boldsymbol{\theta^{*}})$ $\displaystyle\leq\beta\frac{2\phi_{max}\left|\sum_{i=1}^{n}\Delta_{i}\right|}{\left(|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t\right)\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ $\displaystyle\leq\beta\frac{2\phi_{max}\sum_{i=1}^{n}\left|\Delta_{i}\right|}{\left(|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t\right)\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ $\displaystyle\leq\beta\frac{2\phi_{max}L}{\left(|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t\right)\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})}$ $$

As analyzed in Case 1, when the feasible region is non-negligible, $\frac{1}{n}\sum_{i=1}^{n}S(\boldsymbol{x_{i}};\boldsymbol{\theta})>c-\xi$ with high probability. Therefore,

$$ $\displaystyle g(\boldsymbol{\theta^{*}})$ $\displaystyle\leq\beta\frac{2\phi_{max}L}{\left(|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t\right)(c-\xi)}$ $\displaystyle\text{Plug in }\beta\leq\frac{\xi(c-\xi)|\mu_{\Delta}|}{2\phi_{max}L},$ $\displaystyle\leq\frac{\xi|\mu_{\Delta}|}{|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+t}$ $\displaystyle\leq\frac{\xi|\mu_{\Delta}|}{|\mu_{\Delta}\sum_{i=1}^{n}b_{i}^{k}|+L\sqrt{2\sum_{i=1}^{n}(b_{i}^{k})^{2}\log\frac{2}{\delta}}}$ $\displaystyle\leq\frac{\xi}{|\sum_{i=1}^{n}b_{i}^{k}|+\frac{L}{|\mu_{\Delta}|}\sqrt{2\sum_{i=1}^{n}(b_{i}^{k})^{2}\log\frac{2}{\delta}}}$ $$

By Cauchy-Schawarz inequality, we have $\sqrt{\sum_{i=1}^{n}(b_{i}^{k})^{2}}\geq\frac{|\sum_{i=1}^{n}b_{i}^{k}|}{\sqrt{n}}$, therefore,

$$ $\displaystyle g(\boldsymbol{\theta^{*}})$ $\displaystyle\leq\frac{\xi}{|\sum_{i=1}^{n}b_{i}^{k}|+\frac{L}{|\mu_{\Delta}|}\frac{|\sum_{i=1}^{n}b_{i}^{k}|}{\sqrt{n}}\sqrt{2\log\frac{2}{\delta}}}$ $\displaystyle=\frac{\xi}{|\sum_{i=1}^{n}b_{i}^{k}|\left(1+\frac{L}{|\mu_{\Delta}|\sqrt{n}}\sqrt{\log\frac{2}{\delta}}\right)}$ $$

Put together Case 1 and Case 2, we finish the proof.

∎

###### Proof.

## Appendix C Implementation Details

### C.1 Selection of ϕ ^ \hat{\phi} over^ start_ARG italic_ϕ end_ARG

$\frac{\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\theta)\hat{\phi}(\boldsymbol{x}_{i},a_{i},y_{i})}{\sum_{i=1}^{n}S(\boldsymbol{x}_{i};\theta)}$ can be used to represent the ATE estimated by different methods. And the following $\hat{\phi}_{\mathrm{iptw}}$ and $\hat{\phi}_{\mathrm{aiptw}}$,
correspond to the ATE estimated using IPTW and AIPTW, respectively.
*inverse probability of treatment weighting* (IPTW) and *augmented inverse probability of treatment weighting* (AIPTW), respectively.
To clarify, $\phi(\boldsymbol{x}_{i},a_{i},y_{i})$ does not depend on the parameter $\boldsymbol{\theta}$, these estimates are separated from the parametric surrogate model $S$.

$$ $\displaystyle\hat{\phi}_{\mathrm{iptw}}(\boldsymbol{x}_{i},a_{i},y_{i})=$ $\displaystyle\frac{a_{i}}{\hat{e}(\boldsymbol{x}_{i})}y_{i}-\frac{1-a_{i}}{1-\hat{e}(\boldsymbol{x}_{i})}y_{i},$ $\displaystyle\hat{\phi}_{\mathrm{aiptw}}(\boldsymbol{x}_{i},a_{i},y_{i})=$ $\displaystyle\hat{\mu}_{1}(\boldsymbol{x}_{i})-\hat{\mu}_{0}(\boldsymbol{x}_{i})+\frac{a_{i}}{\hat{e}(\boldsymbol{x}_{i})}(y_{i}-\hat{\mu}_{1}(\boldsymbol{x}_{i}))$ $\displaystyle-\frac{1-a_{i}}{1-\hat{e}(\boldsymbol{x}_{i})}(y_{i}-\hat{\mu}_{0}(\boldsymbol{x}_{i})).$ $$

### C.2 Synthetic Data Generation

Let
$\sigma_{X},\sigma_{Y},\rho\in\mathbb{R}$ be fixed constants,
and let $\beta_{1},\beta_{\tau},\omega\in\mathbb{R}^{p}$.
Draw $\boldsymbol{X}$ according to a multi-variate normal distribution, and $A,Y(0),Y(1)$ as follows:

$$ $\displaystyle\boldsymbol{X}\sim\mathcal{MVN}(0,\sigma_{X}^{2}[(1-\rho)I_{p}+\rho 1_{p}1_{p}^{T}]),$ (23) $\displaystyle A|\boldsymbol{X}\sim\textup{Bernoulli}(\sigma(\boldsymbol{X}^{T}\omega)),$ (24) $\displaystyle\epsilon\sim\mathcal{N}(0,\sigma_{Y}^{2}),\hskip 7.22743pt$ (25) $\displaystyle Y(0)=(\sin(10*\boldsymbol{X})+5*\boldsymbol{X}^{2})^{T}\beta_{1}+\epsilon,$ (26) $\displaystyle Y(1)=(\sin(10*\boldsymbol{X})+5*\boldsymbol{X}^{2})^{T}\beta_{1}+\boldsymbol{X}^{T}\beta_{\tau}+\epsilon.$ (27) $$

We set the parameters as follows:

$$ $\displaystyle p$ $\displaystyle=10,\quad\sigma_{X}=\sigma_{Y}=0.1,\quad\rho=0.3,$ $\displaystyle\beta_{1}$ $\displaystyle=[0,0,0,0,2,0,0,0,0,0],$ $\displaystyle\beta_{\tau}$ $\displaystyle=[0.5,0.5,0.5,0.5,0,0,0,0,0,0],$ $\displaystyle\omega$ $\displaystyle=[0,-1\cdot\tilde{\omega},-1\cdot\tilde{\omega},1\cdot\tilde{\omega},1\cdot\tilde{\omega},-2\cdot\tilde{\omega},0,0,0,0].$ $$

The *imbalance parameter*, $\tilde{\omega}\geq 0$, scales the magnitude of the treatment assignment weight vector $\omega$. In particular, we generated two synthetic datasets with (1) no confounding bias ($\tilde{\omega}=0$); and (2) high confounding bias ($\tilde{\omega}=5$). As there is no limitation to generate synthetic data, we set up the total sample size for synthetic data as 5,000 to study the model performance under finite samples and use a balance 50/50 train-test split ratio. As for real-world datasets, we use a 70/30 train-test split.

### C.3 Baseline Implementation

CAPITAL identifies a subgroup by maximizing its size while ensuring the CATE exceeds a predefined threshold. Since subgroup size is not directly controlled in this setting, we vary the CATE threshold and construct the group size vs. ATE curve for comparison.

OWL is originally designed for individual treatment rule estimation but can be viewed as a subgroup identification method without interpretability constraints. It assigns scores between 0 and 1, which we threshold to obtain subgroups of a desired size. We implement OWL using the DTRlearn2 R package.

VT The original VT supports both binary and continuous outcomes, but the commonly used R package aVirtualTwins handles only binary outcomes. Since our experiments require both, we implemented our own VT following (Foster, Taylor, and Ruberg 2011): first estimating treatment effects with a random forest, and then fitting a regression tree to assign subgroup scores. Subgroups of different sizes are obtained by thresholding these scores.

Dragonnet & CT & CF We adopts the implementations from the Python package causalml.

**Table C.3: Number of unbalance features after reweighting by propensity scores obtained by propensity models**
| Dataset | Dragonnet | LR |
| --- | --- | --- |
| Synthetic | 0.10 | 0.10 |
| MIMIC-IV | 2.32 | 1.57 |
| eICU | 1.08 | 1.10 |

### C.4 Model Selection for Nuisance Function Estimation

We consider LR and Dragonnet for estimating the propensity score model. We following the literature (Zang et al. 2023) to use the number of unbalanced features after inverse propensity score weighing as the metric (the lower the better) to select the propensity score model. The results are shown in Table [C.3](https://arxiv.org/html/2504.20908v2#A3.T3), which shows that LR performs better or equivalently than Dragonnet. Thus we select LR as the propensity score model for all datasets.

For the potential outcome model, we consider Dragonnet, CF, and CT, as they are the CATE estimators we are comparing against in subgroup identification. Theoretically, CF improves upon CT by reducing bias and producing smoother decision boundaries through aggregating CTs. As for Dragonnet and CF, prior work (Kiriakidou and Diou 2022) has shown that Dragonnet outperforms CF.

While many other methods exist for nuisance function estimation, our goal is to demonstrate MOSIC’s flexibility rather than prescribe a specific estimator. If more effective models are developed or if a particular method performs better in a given setting, we encourage adapting those for nuisance function estimation.

### C.5 Hyperparameter Tuning

CAPITAL optimizes the policy tree depth with $max_depth\in\{2,3\}$. VT selects decision tree depth for subgroup identification from $max_depth\in\{3,5,7,10\}$. Similarly, CT uses $max_depth\in\{3,5,7,10\}$, while CF considers both tree depth $max_depth\in\{3,5,7,10\}$ and the number of trees $num_tree\in\{5,10,20,50,100\}$. Moreover, Dragonnet tunes the hidden layer size with $hidden_size\in\{50,100,200\}$. Finally, the hyperparameters for our method include $\beta\in\{10^{-2},10^{-3},10^{-4},10^{-5}\}$ and those depending on the implementation of the subgroup identification model. For MOSIC-MLP, we tune the hidden layer size with $hidden_size\in\{50,100,200\}$.

## Appendix D Significance Tests on Read-world data

Tables [D.4](https://arxiv.org/html/2504.20908v2#A4.T4) and [D.5](https://arxiv.org/html/2504.20908v2#A4.T5) present the p-values comparing MOSIC ($\alpha=0.01$) against baseline methods, corresponding to the hypothesis test that MOSIC is different than the baselines. These results directly support the performance trends shown in Figure LABEL:fig:real-world. For instance, in eICU experiments with c=0.5:

- •
MOSIC achieves significantly better ATE than Dragonnet (p=0.0057)
- •
MOSIC maintains significantly fewer unbalanced features than Dragonnet (p=9.2E-04)

**Table D.4: eICU: p-values for ATE and feature balance comparisons between MOSIC and baseline methods.CA*: CAPITAL; DR*: Dragonnet.**
|  |  | $c=0.4$ | $c=0.5$ | $c=0.6$ | $c=0.7$ | $c=0.8$ |
| --- | --- | --- | --- | --- | --- | --- |
| ATE | CT | 0.71 | 0.065 | 0.0011 | 0.0046 | 0.025 |
| CF | 0.19 | 0.0048 | 0.00019 | 0.0012 | 0.0057 |  |
| CA* | - | - | - | 3.5E-05 | - |  |
| DR* | 0.30 | 0.0057 | 5.5E-04 | 0.011 | 0.11 |  |
| OWL | 0.27 | 0.0019 | 0.0010 | 0.084 | 0.16 |  |
| VT | 0.38 | 0.07 | 0.039 | 0.31 | 0.46 |  |
| Balance | CT | 0.16 | 0.58 | 0.56 | 0.65 | 0.04 |
| CF | 0.35 | 0.12 | 0.48 | 0.47 | 0.16 |  |
| CA* | - | - | - | 0.59 | - |  |
| DR* | 0.022 | 9.2E-04 | 0.14 | 0.33 | 0.066 |  |
| OWL | 0.0088 | 0.0077 | 0.13 | 0.22 | 0.046 |  |
| VT | 0.0011 | 0.0060 | 0.036 | 0.23 | 0.39 |  |

**Table D.5: MIMIC: p-values for ATE and feature balance comparisons between MOSIC and baseline methods. CA*: CAPITAL; DR*: Dragonnet.**
|  |  | $c=0.5$ | $c=0.6$ | $c=0.7$ | $c=0.8$ |
| --- | --- | --- | --- | --- | --- |
| ATE | CT | - | - | - | 0.02 |
| CF | 0.42 | 0.68 | 0.35 | 0.16 |  |
| CA* | 0.021 | - | 0.30 | 0.72 |  |
| DR* | 0.039 | 0.10 | 0.21 | 0.30 |  |
| VT | 0.10 | 0.16 | 0.15 | 0.15 |  |
| Balance | CT | - | - | - | 0.63 |
| CF | 0.19 | 5.1E-04 | 1.3E-04 | 1.1E-04 |  |
| CA* | 6.3E-07 | - | 0.17 | 3.3E-05 |  |
| DR* | 0.062 | 0.80 | 0.28 | 0.48 |  |
| VT | 0.033 | 0.0020 | 0.018 | 0.20 |  |

## Appendix E Additional Experiment Results

### E.1 Subgruop Identification on Unconfounded Synthetic data( ω ~ = 0 \tilde{\omega}=0 over~ start_ARG italic_ω end_ARG = 0 )

Figure: Figure E.1.0: True ATE and ATE estimation error across different group sizes on confounded synthetic data($\omega=0$).
Refer to caption: Figures/syn-gamma-0.png

### E.2 MOSIC with interpretability restriction

Figure: Figure E.2.1: Results on eICU. MOSIC-DT and CT use a single tree with depth = 5; MOSIC-Forest and CF use 3 trees with depth = 5. Mean $\pm$ standard error across 100 random splits are reported.
Refer to caption: Figures/eicu_dt_all_methods.png

### E.3 Overlap constraint on synthetic data

We numerically verify that MOSIC can indeed accommodate the overlap constraint (Figure LABEL:fig:syn-error). To assess its effect on estimation error, we conduct experiments on the confounded synthetic data with varying $\alpha$ values. Figure LABEL:subfig:syn-constraint is generated by a random instance of the 100 random train-test splits presented in Figure LABEL:subfig:syn-quan. In Figure LABEL:subfig:syn-constraint, each dot represents an individual, with the x-axis denoting the estimated propensity score, the y-axis representing the true individual treatment effect (ITE), and the vertical lines indicating the desired overlap threshold. Without overlap constraints, MOSIC selects patients with the highest ITEs.

With overlap constraints (Figure LABEL:subfig:syn-constraint, left), MOSIC continues to select patients with large ITEs but systematically excludes those who violate the overlap constraint ($\hat{e}(x)<0.05$ or $\hat{e}(x)>0.95$).
Figure LABEL:subfig:syn-quan further quantifies the effect of excluding patients with limited overlap by reporting the estimation error and the number of unbalanced features.

As the desired overlap threshold $\alpha$ increases, the estimation error of the subgroup ATE decreases, suggesting that enforcing stronger overlap improves estimation reliability.
Meanwhile, stronger overlap also correlates with a lower number of unbalanced features.

### E.4 Extension to Safety, Budget, and Fairness Constraint on Synthetic Data

To demonstrate MOSIC’s flexibility to handle additional constraints, we extend the synthetic data generation process described in Appendix [C.2](https://arxiv.org/html/2504.20908v2#A3.SS2) by introducing a safety, budget, and fairness constraint.
In particular:

- •
Safety constraint: Following Doubleday et al. (2022), each sample is assigned a risk score $r_{i}=1/(1+\exp(10*x_{i}[10]+1))$. We require the average risk of the selected subgroup to be no greater than 0.05:
$\frac{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)r_{i}}{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)}\leq 0.05.$
- •
Budget constraint: Following Qiu, Carone, and Luedtke (2022), each sample is assigned a treatmentd cost value $cost_{i}=(x_{i}[3]+5)/5$. The total cost of the selected subgroup must not exceed half the cost of treating the entire population (assuming a unit cost per sample):
$\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)cost_{i}\leq 0.5*n.$
The test contains 2500 samples, so the total cost limit $0.5*n=1250$ in this case.
- •
Fairness constraint: Let a binary sensitive attribute: $sens_{i}=\mathds{1}(x[3]>0.5)$. We adopt the conditional statistical parity metric (Mehrabi et al. 2021). In our setting, this corresponds to maintaining a sensitive‑group proportion of 0.5 in the selected subgroup.:
$\Big{|}\frac{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)sens_{i}}{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)}-0.5\Big{|}\leq 0.01,$
where we allow a violation of 0.01. We can then convert this to two ratio-form constraints:
$-0.01\leq\frac{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)sens_{i}}{\sum_{i=1}^{n}\mathds{1}(S(x_{i})>0.5)}-0.5\leq 0.01.$

In addition to the group size constraint ($c=0.5$) and overlap constraint($\alpha=0.02$), we progressively add the following constraints: 1) Plus safety constraint; 2) Plus safety and budget constraint; 3) Plus safety, budget, and fairness constraint. Results in Table [E.4.1](https://arxiv.org/html/2504.20908v2#A5.SS4.T1) show that MOSIC can effectively enforce all of these constraints.

**Table E.4.1: Performance on synthetic data ($\tilde{\omega}=5$) under multiple additional constraints.**
| Metric | Group Size and Overlap | +Safety | +Safety+Budget | +Safety+Budget+Fairness |
| --- | --- | --- | --- | --- |
| Group Size | 0.50 ± 0.01 | 0.50 ± 0.05 | 0.48 ± 0.05 | 0.49 ± 0.02 |
| # Unbalance | 0.14 ± 0.40 | 0.20 ± 0.45 | 0.23 ± 0.57 | 0.23 ± 0.45 |
| True CATE | 0.10 ± 0.01 | 0.10 ± 0.01 | 0.09 ± 0.01 | 0.09 ± 0.01 |
| ATE Error | 0.01 ± 0.01 | 0.01 ± 0.01 | 0.01 ± 0.01 | 0.01 ± 0.01 |
| Average Risk | 0.08 ± 0.01 | 0.05 ± 0.01 | 0.05 ± 0.01 | 0.05 ± 0.01 |
| Total Cost | 1377.41 ± 41.57 | 1365.72 ± 144.95 | 1252.78 ± 132.59 | 1268.40 ± 42.44 |
| Sensitive Group Ratio | 0.45 ± 0.03 | 0.45 ± 0.03 | 0.46 ± 0.03 | 0.49 ± 0.02 |

### E.5 Extension to Safety Constraint on eICU

The Glasgow Coma Scale (GCS), a core component of the Sequential Organ Failure Assessment (SOFA) score, measures level of consciousness, with lower scores indicating more severe dysfunction. We targeted patients with GCS $<6$ because this threshold represents the most severe central nervous system dysfunction in the SOFA score, the most commonly-used ICU metric. Without a constraint on the GCS score, we observed that the selected subgroup contains a nontrivial number of patients with GCS$<6$. As observational studies have shown that glucocorticoids may exacerbate neural system damage, such a subgroup raises a safety concern. Therefore, we required that the proportion of patients with GCS $<6$ remain below 0.05. We did not enforce a strict exclusion (i.e., a threshold of 0) because it is possible that some patients with severe neural dysfunction could still derive meaningful benefit from the treatment.

To impose the interpretability requirement, we implement MOSIC-DT for this setting. Similarly, we run experiments on 100 random train-test splits and use 5-fold cross-validation to select the tree depth among {3,5,7}. Results in Table [2](https://arxiv.org/html/2504.20908v2#S5.T2) show MOSIC can effectively exclude patients with GCS$<6$. Additionally, an example (Figure LABEL:fig:eicu_dt_cns_constraint_example) shows that MOSIC indeed learned to exclude such high-risk patients.

## Appendix F Evaluation on Binary Subgroup Setting

Unlike approaches that assume the existence of two or a finite number of subgroups with distinct ATEs, our study focuses on identifying a subset of the population with maximal ATE under real-world constraints, without making structural assumptions about the underlying heterogeneity. This design enables our method to generalize to continuous or complex heterogeneity. Nevertheless, we recognize that such structural assumptions, such as that binary subgroups with distinct ATE exist, are plausible in certain real-world scenarios. In such cases, the real-world constraint will naturally introduce a trade-off for identification performance.

To illustrate, we modify the DGP of synthetic data by replacing $Y(1)$ in Appendix [C.2](https://arxiv.org/html/2504.20908v2#A3.SS2) to:

$$ $Y(1)=(\sin(10*\boldsymbol{X})+5*\boldsymbol{X}^{2})^{T}\beta_{1}+\mathds{1}(\boldsymbol{X}^{T}>0.05)\beta_{\tau}+\epsilon.$ $$

That says, only patients with covariates $>0.05$ at positive $\beta_{\tau}$ indices receive a positive effect; others receive none. This yields a positive subgroup comprising  68% of samples. Using this DGP, we test MOSIC with beta = 1e-5, alpha = 0, and c in {0.6, 0.7, 0.8}. The precision and recall of subgroup identification are evaluated. Results in Table [F.0.1](https://arxiv.org/html/2504.20908v2#A6.SS0.T1) are reported as mean $\pm$ standard deviation over 100 runs. Performance aligns with theory: precision/recall degrade when $c$ (the group size constraint) exceeds/falls below the true subgroup size, reflecting constraint-driven trade-offs.

**Table F.0.1: Performance of binary subgroup identification**
| c | ATE | Group Size | Precision | Recall |
| --- | --- | --- | --- | --- |
| 0.6 | 0.92$\pm$0.07 | 0.60$\pm$0.01 | 0.93$\pm$0.05 | 0.83$\pm$0.05 |
| 0.7 | 0.84$\pm$0.05 | 0.70$\pm$0.01 | 0.88$\pm$0.04 | 0.92$\pm$0.04 |
| 0.8 | 0.75$\pm$0.04 | 0.80$\pm$0.01 | 0.80$\pm$0.03 | 0.96$\pm$0.04 |

## Appendix G Evaluation of Type I error

Although our primary focus is on multi-constraint subgroup identification rather than statistical inference, we also evaluate Type I error of the identified subgroup. We use a data-splitting approach to test whether identified subgroups arise spuriously under the null hypothesis, where all individuals have zero treatment effect. Using data splitting, Type I error of the selected subgroup can be evaluated after we build the subgroup assignment model. Specifically:

- 1.
Split the dataset (e.g., 50-50) into training (subgroup selection) and holdout (inference);
- 2.
Train MOSIC on the training set to learn the subgroup model;
- 3.
Apply the model to the holdout set to identify the subgroup, then compute its subgroup ATE, denoted as $ATE_{hold_out}$
- 4.
Test the null hypothesis on the holdout set:
(a)
Construct the distribution of subgroup ATE under the null (here by directly sampling from the test distribution, with bootstrap subsample size equaling the target subgroup size specified by the parameter $c$)
(b)
Compare $ATE_{hold_out}$ to this distribution, and determine whether the null hypothesis should be rejected.
- 5.
Repeat steps 1-4 on additional synthetic data instances, aggregate results, and estimate type-I error.

We generated 100 synthetic datasets using the same DGP in Appendix [C.2](https://arxiv.org/html/2504.20908v2#A3.SS2) except that $\beta_{\tau}$ is set to $\boldsymbol{0}$. This aligns with the null hypothesis: all individuals have zero treatment effect. For each instance, we set the bootstrap iterations to 10000. Type I error rate is computed as the proportion of instances in which the null hypothesis is rejected by 5% significance level.

As shown in Table [G.0.1](https://arxiv.org/html/2504.20908v2#A7.SS0.T1), the Type I error rate increases when the parameter $c$ is small, consistent with theoretical expectations. Smaller subgroups lead to higher variance in ATE estimates, highlighting a fundamental trade-off between real-world constraints and statistical reliability.

**Table G.0.1: Type I error under different constraints for group size.**
| c | ATE | Group Size | Type I error |
| --- | --- | --- | --- |
| 0.4 | -0.0011 $\pm$ 0.0401 | 0.3991 $\pm$ 0.0141 | 0.12 |
| 0.6 | -0.0018 $\pm$ 0.0298 | 0.6020 $\pm$ 0.0141 | 0.00 |
| 0.8 | -0.0000 $\pm$ 0.0212 | 0.8038 $\pm$ 0.0121 | 0.00 |

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. How do you ensure unbiased AIPTW estimates if subgroup selection and nuisance estimation share data? Was any cross-fitting or sample splitting applied?
2. Can you evaluate overlap on test data, e.g., showing the distribution of estimated propensities or proportion outside [0.05, 0.95] within selected subgroups?
3. Why not compare against overlap-weighted baselines or modern CATE learners like DR-/R-learner or BART?
4. Could you report uncertainty for subgroup ATEs (e.g., influence-function or bootstrap CIs) and include a sensitivity analysis for unmeasured confounding in the ICU data?
soundness: 2
strengths: The paper proposes a clear optimization framework for subgroup identification under multiple constraints. The formulation is well defined and the paper is clearly written. The experiments are consistent across datasets and show some improvement in subgroup ATE and balance, though mainly within controlled and limited settings
summary: The paper proposes MOSIC, a model-agnostic framework for identifying treatment-benefit subgroups under multiple constraints such as minimum size and overlap. It formulates a constrained optimization with ReLU-gated Lagrangian penalties and provides local min–max guarantees. Experiments on synthetic data and two ICU datasets show that MOSIC achieves higher subgroup treatment effects with fewer imbalance violations than prior subgroup discovery methods.
weaknesses: 1. The largest concern is the limited baselines: the paper includes CAPITAL, causal trees, causal forests, and outcome-weighted learning, but exclude stronger HTE estimators such as DR-learner, R-learner, BART, and locally centered causal forests. These are standard in healthcare and could change the ranking. Overlap-weighted estimators are also missing and align directly with the overlap constraint. It would also be helpful to discuss recent constrained policy methods from off-policy evaluation, safe RL, and constrained contextual bandits that optimize under budgets or risk

2. Experimental setup lacks robustness and diagnostic depth.
The evaluation relies on AIPTW without clear cross-fitting, which can bias subgroup ATEs when selection and estimation share nuisance models. The reliability metric is only the count of covariates with standardized mean difference>0.2 after IPTW, which is coarse and does not measure magnitude or joint imbalance. Overlap enforcement is only evaluated during training, not validated on test sets. 

3. Causal validity in ICU cohorts is not well supported.
The treatment window (10 hours before to 24 hours after ICU admission) risks immortal time bias if early deaths are excluded. Confounding from indication severity is possible, especially for steroids, yet no sensitivity analysis or negative control outcomes are reported. The added safety constraint on GCS is useful, but without external validation or clinician input, it’s unclear whether the selected subgroups are clinically meaningful or plausible.

---

confidence: 3
contribution: 2
initial_rating: 4
presentation: 3
questions: - The imbalance criterion should be SMD > 0.1 rather than 0.2 according to Austin (2009).  
- The abbreviation “LR” in line 079 likely refers to Logistic Regression, which should be clearly defined when it first appears.
- What is the computational complexity and scalability of $\gamma$-GDA when applied to high-dimensional data or a large number of constraints?
- Can the authors discuss how subgroup interpretability is maintained when using black-box backbones such as MLPs and whether any post-hoc explainability tools were explored?
soundness: 3
strengths: - Introduces an optimization-based framework that integrates multiple constraints directly into subgroup identification, addressing the limitations of post-hoc filtering methods.  
- Provides theoretical guarantees for convergence, feasibility, and approximate satisfaction of nonconvex constraints.  
- Empirical evaluations show improvements in subgroup ATE, size control, and covariate balance on both synthetic and real-world datasets.  
- Offers an open-source implementation with transparent algorithmic details, enhancing reproducibility.
summary: This paper introduces MOSIC, a unified optimization framework for identifying optimal treatment subgroups under multiple practical constraints such as subgroup size and propensity overlap. Unlike conventional two-step CATE-based approaches, MOSIC formulates the task as a constrained min–max optimization and solves it using a modified Gradient Descent–Ascent algorithm. The authors establish theoretical guarantees for feasibility and local optimality and validate the approach on both synthetic and ICU datasets (eICU, MIMIC-IV). The framework is model-agnostic, compatible with various CATE estimators and model architectures, and can incorporate additional constraints such as fairness and safety.
weaknesses: - The method’s scalability in high-dimensional settings is not fully addressed, especially given the use of $\gamma$-GDA with many constraints.  
- Limited exploration is provided on uncertainty quantification or confidence intervals for subgroup ATEs, which is critical for decision-making in clinical applications.
- Numerical stability concerns of $\gamma$-GDA are only discussed theoretically; empirical results showing training dynamics (e.g., convergence plots, constraint violation over iterations) would substantiate the stability claim.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: 1. Line 156: The paper claims that the method naturally accommodates more general linear and ratio-form constraints. Could the authors clarify what kinds of constraints can be naturally supported by MOSIC, and what kinds cannot? Additionally, how does the computational or optimization complexity scale as more constraints are added, does it grow linearly, or is there a risk of exponential blow-up? I also wonder would there be a finite-sample guarantees with constraints? 

2. I understand that strict local min–max is a standard target for nonconvex optimization. However, since we are dealing with patient data, robustness is a critical concern. You mention running experiments over 100 random train/test splits, but do you also evaluate stability across different random initialization seeds of the optimizer (not just different data splits)? In other words, do different seeds lead to meaningfully different subgroups being selected, even if the average ATE is similar? Is the identified subgroup S(x) stable across initializations?

3. As mentioned in the weaknesses section, providing runtime and interpretability analysis would significantly strengthen the paper, particularly for clinical decision. Even an approximate runtime comparison against two-stage baselines, and at least one example of a human-readable decision-tree subgroup rule, would be very helpful.
soundness: 3
strengths: **Good Motivation**

The paper transfers the subgroup identification problem from the standard two-stage CATE estimation and filtering into an optimization problem with constraints. The idea is elegant, and this kind of shift is decision-focused and real-world aware, aligning with the emerging trend in causal machine learning.

**Theoretical Proof**

I like that the paper provides several theoretical lemmas to support their algorithms, e.g., the strict local min–max point in Lemma 2. I didn't check the full details of the proof, but they seem correct (although I could be mistaken). The authors also provide experiments on both synthetic and real-world data to show the proposed method’s performance compared with baselines.
summary: This paper proposes MOSIC, a model-agnostic optimization framework for treatment subgroup identification under multiple real world constraints, including minimum subgroup size, propensity score overlap, and optionally safety or fairness constraints. Unlike standard two-stage approaches that first estimate CATE and then apply post-hoc thresholding or filtering, MOSIC directly formulates subgroup identification as a constrained optimization problem. Experiments on both synthetic datasets and real data show that MOSIC consistently achieves higher subgroup treatment effects while satisfying the specified constraints.
weaknesses: My main concern is about the practical utility of the proposed methods. The reason why the standard two-stage CATE methods are widely accepted is because they are easy to understand, easy to use, and easy to interpret. While I appreciate the theoretical and empirical results that MOSIC presents in the paper, the key question is: is the gain worth the additional complexity?

**Easy to understand?**

The idea is straightforward: shifting from estimation to optimization. However, the transition from Problem I to Problem IV sometimes feels reactive rather than intentional. Section 4 in particular reads like a narrated debugging process (“we tried X but it was unstable, so we added Y...”), rather than a clean forward design narrative. This hurts clarity, especially for theoretical readers.

**Easy to use?**

There is no discussion of computational cost or efficiency. The two real-world datasets used have sample sizes of 13,361 and 6,516, but there is no reporting of runtime or resource overhead relative to simpler two-stage baselines. This omission limits the reader’s ability to judge whether the method is practical for real-time or large-scale deployment. And for many clinicians who may need to use this method, will it still be practical if they do not have sufficient computational resources?

**Easy to interpret?**

While the method supports decision-tree backbones in Section 5.3, the paper never shows a concrete learned subgroup rule. For example, a clinician might prefer an interpretable patient subgroup like: “patients aged > 60 should not take the drug.” But this kind of explicit rule is missing from the paper. For deployment and clinical trust, this is a significant gap.


Decision: Reject