[3]Yonatan Mintz

1]Tepper School of Business, Carnegie Mellon University, 5000 Forbes Avenue,
Pittsburgh, 15213, PA, USA

2]Maternal-Fetal Medicine University of Wisconsin – Madison, 1010 Mound
Street, Madison, 53715, WI, USA

[3]Industrial and Systems Engineering, University of Wisconsin – Madison, 1513
University Avenue, Madison, 53706, WI, USA

# Applications of 0-1 Neural Networks in Prescription and Prediction

Vrishabh Patil<EMAIL_ADDRESS>Kara Hoppe<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  *

###### Abstract

A key challenge in medical decision making is learning treatment policies for
patients with limited observational data. This challenge is particularly
evident in personalized healthcare decision-making, where models need to take
into account the intricate relationships between patient characteristics,
treatment options, and health outcomes. To address this, we introduce
prescriptive networks (PNNs), shallow 0-1 neural networks trained with mixed
integer programming that can be used with counterfactual estimation to
optimize policies in medium data settings. These models offer greater
interpretability than deep neural networks and can encode more complex
policies than common models such as decision trees. We show that PNNs can
outperform existing methods in both synthetic data experiments and in a case
study of assigning treatments for postpartum hypertension. In particular, PNNs
are shown to produce policies that could reduce peak blood pressure by 5.47 mm
Hg (p=0.02) over existing clinical practice, and by 2 mm Hg (p=0.01) over the
next best prescriptive modeling technique. Moreover PNNs were more likely than
all other models to correctly identify clinically significant features while
existing models relied on potentially dangerous features such as patient
insurance information and race that could lead to bias in treatment.

###### keywords:

mixed-integer programming, personalized healthcare, neural networks, causal
inference, interpretability.

##  1 Introduction

Advances in machine learning and data science theory with the advent of
Artificial Neural Networks (ANNs) have brought about substantial progress in
learning complex tasks in the context of classification, prediction, and
prescription. ANNs have become ubiquitous tools to address numerous
classification tasks in domain areas such as natural language processing and
image recognition. The burgeoning availability and accessibility of data in
recent years have only improved the utility of ANNs. Despite the benefits of
acquiring such data, challenges such as the high cost and time required to
collect data, the presence of biased data, or incomplete data persist [1, 2].
These challenges are particularly prevalent in areas of interest to the
operations and management science community; for example, epidemiological
research on the prescription of drugs or the affects of smoking necessitate
obtaining limited data through observational studies as it is often unethical
to obtain a large amount of patient data through randomized controlled trials
[3, 4]. In addition to limitations due to ethical concerns, studies involving
a small subset of the population, such as treating postpartum birthing persons
for hypertension, are a naturally data-scarce settings. Moreover, using data
that may be biased or incomplete can lead to inaccurate or unreliable models
[5]. The consequences of enacting policies derived from such models can be
dire. In the context of medicine and healthcare, these could include lawsuits,
increased costs of insurance, and adverse health effects to individuals
targeted by the policy. This can be particularly problematic in public health,
where bias in data for a specific demographic leads to discriminatory
policies.

In data-rich settings, deep neural networks, i.e, networks with a large number
of hidden layers, have taken center stage in learning modern complex tasks.
However, current state-of-the-art optimization tools used to train ANNs, such
as stochastic gradient descent (SGD), fall short when data availability is
limitted [6, 7, 8]. Furthermore, criticisms against the interpretability of
deep neural networks are increasingly prevalent in domains such as medicine,
public health, and social intervention where fairness and accountability are
critical [9, 10, 11]. Recent techniques emerging from operations research
literature have explored the use mixed-integer programming to address some of
the limitations of ANNs trained with SGD [12, 13, 14]. These works have
established the potential of integer programming in data-limited setting.
However, there is currently no consensus on modelling an appropriate loss
function for the mixed-integer programming (MIP) formulations.

In this paper, we explore the problem of training shallow neural networks with
limited data using MIP with a focus on personalized prescription problems.
While prediction remains a rich field for learning tasks, decision makers in
several fields are motivated to make personalized decisions tailored to
individuals. These problems arise in domains such as personalized online
advertising [15], direct mail marketing campaigns [16], and personalized
medical treatment assignments [17]. We provide empirical insight into
designing personalized prescription policies in the medical domain where
conducting randomized control trials are expensive, unethical, or otherwise
infeasible. Specifically, we consider the tasks of deciding treatments for
hypertension in postpartum birthing persons. Our methods are well-suited to
solve such problems since the availability of data is limited, and the need
for understanding which patient characteristics influence policy design is
high. In addition to these experiments, in the appendix we provide
supplementary experiments showing the use of our prescriptive networks for
designing optimal policies to prescribing warfarin for patients susceptible to
blood clots. Additionally, in the appendix, we produce empirical results
demonstrating the use of our models for prediction using the MNIST digit
recognition problem for the 2-digit case.

###  1.1 Problem Description

Machine learning is commonly used to solve prediction problems, which involve
building models that can use data to predict future events or outcomes. These
models learn from labeled data to identify patterns and relationships, which
can then be applied to new, unseen data. Successfully solving prediction
problems can provide valuable insights for decision-making in a variety of
disciplines.

Alternately, prescription problems are a growing array of tasks in
personalized decision-making. These problems arise when managers and decision-
makers need to make general policy decisions that account for the
heterogeneous responses to the policy among different individuals in the
population. In such cases, prediction alone is not sufficient as individuals
may respond to the same treatments differently. While adaptive or
reinforcement learning can be used to address prescription problems [18, 19,
20], there are cases where neither are entirely appropriate – adaptive
learning is a promising tool in online learning settings, and reinforcement
learning is particularly useful in data-rich settings. However, given limited
observational data, we aim to optimize over expected outcomes directly. In
this paper, we consider the problem of prescribing the best treatment to an
individual as a function of their characteristics such that the outcome of
such a prescription is optimized.

These problems are particularly challenging because they require us to
determine not only what the best treatments are, but also why they are the
best treatments. In other words, we need to understand the causal
relationships between treatments and outcomes. Indeed, we use techniques from
causal inference to define an objective that maximizes (or minimizes) the
outcome by estimating the counterfactual information from historical
observations. By doing so, we can gain a better understanding of the causal
relationships between actions and outcomes, allowing us to make better policy
decisions. For instance, when prescribing blood-pressure medication to
postpartum patients, we may be interested in ensuring that the model
identifies those features that attribute to a positive response to receiving
treatment as well as those that attribute to a negative response.

###  1.2 Contributions

The overarching goal of this paper is to provide a structured framework to
training ANNs using MIP techniques for prescription problems.

We make several methodological contributions to these problem domains. In the
realm of prescription problems, we design an optimization framework that uses
MIP-based ANNs, which we refer to as Prescriptive Neural Networks (PNN). PNNs
designs interpretable policies with the goal of optimizing population average
outcomes given limited observational data. This problem is non-trivial because
we need to account for the outcomes of unobserved treatments, as well as
incorporate causal effects into our framework. To address this challenge, the
PNN optimizes over-estimates of the conditional average treatment effects
(CATE) to account for unobserved counterfactual outcomes. To implement the
framework into commercial solvers, we propose a novel mixed-integer linear
objective function for the prescription problem that incorporates outcomes
estimated by the Inverse Propensity Weighting (IPW) method, Direct Method (DM)
or the Doubly Robust (DR) estimator. We also provide statistical guarantees
for the proposed formulations that address both prediction and prescription
problems. In particular, our analysis shows that these frameworks produce
statistically consistent estimators in both cases, to the best of our
knowledge this is one of the first derived statistical guarantees for ANNs
trained with MIP.

We validated our framework through a series of experiments on both simulated
and real-world data. To validate PNNs, we benchmark our methods against state-
of-the-art prescriptive methods and heterogeneous treatment effect estimation
techniques with simulated experiments common in the prescription literature
[21] along with two case studies in personalized healthcare. The first case
study entails designing prescriptive policies to treat postpartum patients
against hypertension. The second case study, which we present in the appendix,
involves incorporating prescription in personalized warfarin dosing. Our
contributions also extend to assessing the interpretability of the PNN, and
its ability to identify key features that impact in the models’ decisions. By
doing so, we derive a sense of our model’s strengths and identify where it is
most apt to use prescriptive neural networks over others. In particular, we
show that PNNs are more likely to identify relevant clinical features for
treatment decisions, unlike existing methods that may use biased surrogate
features such as patient insurance information.

While not in the main paper but in the appendix, we also present contribution
to the use of MIP ANNs in the realm of prediction. We develop a novel MIP
formulation of the negative log-likelihood loss (NLL) that incorporates a
softmax activation function applied over the output layer of an ANN formulated
by mixed-integer constraints. Existing literature that explores using MIP to
train ANNs has considered a Support Vector Machine loss, minimizing empircal
classification loss, and other linear and non-linear loss functions [12, 13,
14]. The softmax function is an appealing activation function since it
considers the output of the neural network as a probabilistic distribution
over predicted output classes [22]. By incorporating this well-studied loss
function from the literature into MIP-based ANNs, our work helps bridge the
divide between MIP-based training methods and classical training methods. To
evaluate the performance of our MIP-based ANN with the NLL objective, we
conduct numerical experiments using the MNIST digit recognition dataset. We
benchmark our performance against state-of-the-art stochastic gradient descent
(SGD) based algorithms. These experiments are conducted over a variety of ANN
architectures and hidden-layer activation functions.

##  2 Literature Review

Our work contributes to several important streams of literature in machine
learning and operations research. The MIP formulations we develop in this
paper builds upon the existing literature on using integer programming and
constraint programming to train ANNs. Our work on applying this framework for
personalized prescription problems parallel previous papers that consider the
intersection of optimization and statistical inference. The first stream of
literature pertains to training ANNs using MIPs. In the second stream, we
review works that estimate heterogeneous causal effects in the context of
learning treatments effects across subsets of the population with causal trees
and forests. Along these lines, we also discuss non-tree methods used in
statistical learning. Additionally, we present related work on learning
prescriptive trees that optimize using MIP formulations.

Using integer programming to train neural networks has gained some popularity
in recent years. [12] use Constraint Programming, MIP, and a hybrid of the two
to train Binarized Neural Networks (BNNs), a class of neural networks with the
weights and activations restricted to the set {-1, +1}. To relax the binarized
restriction on weights, [14] propose a formulation to train Integer Neural
Networks with different loss functions. [13] developed a formulation for
integer weights, and offer an iterative data splitting algorithm that uses the
k-means classification method to train on subsets of data with similar
activation patterns. For problems with continuous activation functions, [23]
introduce a mixed-integer linear programming (MILP) formulation for deep
neural networks to count the number of linear regions represented by it. While
the authors do not use the formulation to train neural networks, they utilize
a systematic one-tree approach [24] to count integer solutions of the MILP.
MIPs have also been used for already trained networks to provide adversarial
samples that can improve network stability [25, 26, 27]. Additionally, there
has been research directed towards using ANNs to solve MIPs [28].

In the context of statistical learning based approaches for causal inference,
[21] proposed causal tree estimators that use recursive partitions to estimate
conditional average treatment effects. The causal tree structure uses the
Inverse Propensity Weighting [29] to estimate counterfactual outcomes. [30]
build on the concept of causal trees by extending estimation using random
forests with a causal structure, and provide asymptotic confidence intervals
for their work. [31] further build on causal forests and introduce Pollinated
Transformed Outcome forests for the treatment effect estimation problem in
high dimensions. They additionally apply modifications to boosting and MARS
(Multivariate Adaptive Regression Splines) to use in the context of structured
causal problems. Along the lines of estimating counterfactual estimation, [32]
provides a doubly robust framework that combines inverse propensity weighting
methods with “direct methods”, i.e, directly estimating counterfactual
outcomes using learning models with a given sample. The same paper proposes
methodologies for devising optimal treatment policies using the doubly robust
framework. A study on using deep neural networks for semiparametric inference
was conducted by [33], which apply theoretical results on treatment effect and
expected welfare estimation.

We now discuss tree-based approaches for personalized prescription problems.
The works reviewed in this avenue of research build on the classification
trees [34]. A key motivation presented in literature for personalized trees is
to provide a sense of interpretability to prescription problems. [35] proposes
personalization trees with a mixed-integer optimization framework to learn to
choose the optimal treatment from a discrete set of treatments. [36] modify
this approach with by incorporating regularization for variance-reduction in
the objective function of the problem. The authors use coordinate descent with
multiple starts to solve the optimization problem numerically. A prescriptive
tree approach presented by [37] builds on a MIP-based classification tree [38]
formulated as a flow network. Both [35] and [36] provide structure for
learning optimal trees without explicitly estimating counterfactual outcomes,
while the MIP formulation detailed by [37] incorporates inverse propensity
weighting methods, direct methods, or doubly robust methods in their
optimization problem. [39] extend the work of [36] to the problem of learning
prescription policies using estimated counterfactual information, and to
problems with continuous-valued treatments. Finally, a new type of
prescription trees is presented by [40], which uses threshold-based methods to
devise decision-rules that can be used in-tandem with the optimal policy trees
[39] under budget constraints.

##  3 Mixed-Integer Programming Formulation of ANNs

In this section we consider MIP formulations for ANNs in the prescriptive
setting. However, in the appendix, we discuss how the formulation can be
adapted to the predictive setting. Consider a personalized decision-making
problem where the goal is to assign the best treatment $t\in\mathcal{T}$ for
an individual given their covariates
$X\in\mathcal{X}\subseteq\mathbb{R}^{\mathcal{F}}$ such that the individual’s
potential outcome $Y(t)\in\mathcal{Y}\subseteq\mathbb{R}$ is either minimized
or maximized over the treatment group. We assume both
$\mathcal{X},\mathcal{Y}$ are compact sets and that $\mathcal{T}$ is a finite
set. In this problem, decision-makers use observational data to design a
personalized treatment policy $s:\mathcal{X}\mapsto\mathcal{T}$ that should
then be able to prescribe an optimal treatment $s(X)$ for an individual based
on their covariates $X$. In the case where all of the potential outcomes
$\\{Y_{i}(t)\\}_{t\in\mathcal{T}}$ for an individual $i$ are known, the
problem of learning an optimal policy $\hat{s}$ is a classification problem
where the model classifies the patient to the treatment group that maximizes
or minimizes individual’s outcome. Recall that we are interested in the
observational case when prescribing a patient every possible treatment to
procure all potential outcomes can be impractical due to ethical concerns or
resource constraints. In this problem, we observe an outcome
$\\{Y_{i}(t_{i})\\}$ for individual $i\in\\{1,\ldots,n\\}$ having received
treatment $t_{i}$, but not their counterfactual outcomes
$\\{Y_{i}(t)\\}_{t\in\mathcal{T}\setminus t_{i}}$. Our goal is to learn a
policy $\hat{s}$ from $n$ independent and identically distributed observations
of participants
$\\{(X_{1},t_{1},Y_{1}(t_{1})),\ldots,(X_{n},t_{n},Y_{n}(t_{n}))\\}$.

We use the notation $[n]_{m}=\\{m,m+1,...,n\\}$ for any two integers $n>m$.
For our problem, we will think of policies $s$ being parameterized by
$\theta\in\Theta$, where $\Theta$ is a compact set, and our goal will be to
find $\theta^{*}$ such that for any $X\in\mathcal{X}$, $s(X;\theta^{*})$ is
able to closely predict the optimal treatment. Our goal will be to estimate
$\theta^{*}$ by finding a $\hat{\theta}$ that minimizes an empirical loss
$\mathcal{L}_{n}:\mathcal{T}\times\mathcal{Y}\mapsto\mathbbm{R}$ that is
$\hat{\theta}=\mathrm{argmin}_{\theta\in\Theta}\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}_{n}(s(X_{i};\theta),Y_{i})$.
We assume $s$ is an ANN that can be expressed as a functional composition of
$L$ vector valued functions that is $s(X;\theta)=h_{L}\circ
h_{L-1}\circ...\circ h_{1}\circ h_{0}(x)$. Here each function
$h_{\ell}(\cdot,\theta_{\ell})$ is referred to as a layer of the neural
network, we assume that $h_{0}:\mathcal{X}\mapsto\mathbbm{R}^{K}$,
$h_{L}:\mathbbm{R}^{K}\mapsto\mathcal{T}$, and
$h_{\ell}:\mathbbm{R}^{K}\mapsto\mathbbm{R}^{K}$, we denote each component of
$h_{\ell}$ as $h_{k,\ell}$ for $\ell\in[L-1]_{0}$ and $k\in[K]_{1}$, and each
component of $h_{L}$ as $h_{t,L}$ for $t\in\mathcal{T}$. We refer to these
components as units or neurons. We refer to the index $\ell$ as the layer
number of the ANN, and the index $k$ as the unit number. The layer
corresponding to the index $L$ is referred to as the output layer, and all
other indices are referred to as hidden layers. Each unit has the functional
form: $h_{k,0}=\sigma(\alpha_{k,0}^{\top}X+\beta_{k,0})$,
$h_{k,\ell}=\sigma(\alpha_{k,\ell}^{\top}h_{\ell-1}+\beta_{k,\ell})$ for
$\ell\in[L-1]_{1}$, $h_{t,L}=\varphi(\alpha_{t,L}^{\top}h_{L-1}+\beta_{t,L})$,
where $\sigma:\mathbbm{R}\mapsto\mathbbm{R}$ is a non-linear function applied
over the hidden layers $\ell=0,...,L-1$ called the activation function,
$\varphi:\mathbbm{R}\mapsto\mathbbm{R}$ is a potentially different non-linear
function applied over the output layer $L$, and $\alpha$ and $\beta$ are the
weights and biases matrices respectively where $(\alpha,\beta)=\theta$. In our
work, we consider $\sigma$ to be the binary activation and we assume for any
input $z\in\mathbbm{R}$, $\sigma(z)=\mathbbm{1}[z\geq 0]$, where $\mathbbm{1}$
is the indicator function. The notation $\alpha_{a,b,\ell}$ indicates the
weight being applied to unit index $a$ in layer $\ell-1$ for the evaluation of
unit index $b$ in layer $\ell$. Likewise, $\beta_{k,\ell}$ indicates the bias
associated with unit index $k$ in layer $\ell$ of the network.

We show how the above training problem and model can be formulated as a MIP.
The key to our reformulation is the introduction of decision variables
$h_{i,k,\ell}$ which correspond to the unit output of unit $k$ in layer $\ell$
when the ANN is evaluated at data point with index $i$. Having these decision
variables be data point dependant, and ensuring that
$\alpha_{k^{\prime},k,\ell},\beta_{k,\ell}$ are the same across all data
points forms the backbone of our formulation. Thus, if $X_{i,d}$ denotes the
$d^{th}$ feature value of individual $i$, the general form of the optimization
problem is:

$\displaystyle\hat{\theta}_{n}=\operatorname*{arg\,min}_{\theta\in\Theta}$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sum_{t\in
T}\mathcal{L}_{n}(h_{i,j,t},y_{i,j})$ (1a) subject to $\displaystyle
h_{i,k,0}=\sigma(\sum_{d\in\mathcal{F}}\alpha_{d,k,0}X_{i,d}+\beta_{k,0}),$
$\displaystyle\forall\ k,i\in[K]_{1}\times[n]_{1},$ (1b) $\displaystyle
h_{i,k,\ell}=\sigma(\sum_{k^{\prime}=1}^{K}\alpha_{k^{\prime},k,\ell}h_{i,k^{\prime},\ell-1}+\beta_{k,\ell}),$
$\displaystyle\forall\ \ell,k,i\in[L-1]_{1}\times[K]_{1}\times[n]_{1},$ (1c)
$\displaystyle
h_{i,t,L}=\varphi(\sum_{k^{\prime}=1}^{K}\alpha_{k^{\prime},j,L}h_{i,k^{\prime},L-1}+\beta_{t,L}),$
$\displaystyle\forall j,i\in\mathcal{T}\times[n]_{1}.$ (1d)

In this section we will discuss how the constraints of (1) can be represented
using a MIP formulation. We will also discuss how different forms of common
causal inference loss functions $\mathcal{L}_{n}$ can be formulated using MIP
constraints. We will conclude by proving that the resulting formulation
provides consistent estimates of the optimal policy. That is, that as
$n\rightarrow\infty$, $\hat{\theta}_{n}\overset{p}{\rightarrow}\theta^{*}$.
This property is crucial, since it indicates that as decision makers gather
more data our technique will generate more effective policies. Moreover, it is
not clear if many ANN architectures trained using stochastic gradient based
approaches exhibit this property in general.

###  3.1 Binary Activations

In this section, we present a MIP reformulation of the constraints of (1) for
binary activations in the hidden layers that can be solved using commercial
solvers. Our formulation can be seen as an extension of prior MIP formulations
for ANN training by using continuous instead of integer valued model weights
[12]. Specifically, we relax the bounds such that
$\alpha_{d,k,0},\alpha_{k^{\prime},k,l},\alpha_{k^{\prime},j,L}\in[\alpha^{L},\alpha^{U}]$.
We rewrite Constraints (1b),(1c),(1d) for a single unit as
$h_{i,k,0}=\mathbbm{1}[\sum_{d\in\mathcal{F}}\alpha_{d,k,0}X_{i,d}+\beta_{k,X0}\geq
0]$,
$h_{i,k,\ell}=\mathbbm{1}[\sum_{k^{\prime}=1}^{K}\alpha_{k,k^{\prime},\ell}h_{i,k^{\prime},\ell-1}+\beta_{k,\ell}\geq
0]$,
$h_{i,j,L}=\sum_{k^{\prime}=1}^{K}\alpha_{k^{\prime},j,L}h_{i,k^{\prime},L-1}+\beta_{j,\ell}$
respectively. Note that for all layers that are not the input layer, the above
constraints contain bi-linear products which make this formulation challenging
to solve. We propose the following reformulation:

######  Proposition 1.

In the binary activation case, the Constraints (1b),(1c),(1d) can be
reformulated as a set of MIP constraints. Specifically for all
$k\in[K]_{1},i\in[n]_{1}$ Constraint (1b) can be reformulated as:

$\displaystyle\sum_{d\in\mathcal{F}}(\alpha_{d,k,0}X_{i,d})+\beta_{k,0}\leq
Mh_{i,k,0},$ (2)
$\displaystyle\sum_{d\in\mathcal{F}}(\alpha_{d,k,0}X_{i,d})+\beta_{k,0}\geq\epsilon+(-M-\epsilon)(1-h_{i,k,0}),$
(3)

For all, $\ell\in[L-1]_{1},k\in[K]_{1},i\in[n]_{1}$ Constraints (1c) can be
reformulated as:

$\displaystyle\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},k,\ell})+\beta_{k,\ell}\leq
Mh_{i,k,\ell},$ (4)
$\displaystyle\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},k,\ell})+\beta_{k,\ell}\geq\epsilon+(-M-\epsilon)(1-h_{i,k,\ell}),$
(5) $\displaystyle
z_{i,k^{\prime},k,\ell}\leq\alpha_{k^{\prime},k,\ell}+M(1-h_{i,k^{\prime},\ell-1}),\
\forall\ k^{\prime}\in[K]_{1},$ (6) $\displaystyle
z_{i,k^{\prime},k,\ell}\geq\alpha_{k^{\prime},k,\ell}-M(1-h_{i,k^{\prime},\ell-1}),\
\forall\ k^{\prime}\in[K]_{1},$ (7) $\displaystyle-
Mh_{i,k^{\prime},\ell-1}\leq z_{i,k^{\prime},k,\ell}\leq
Mh_{i,k^{\prime},\ell-1},\ \forall\ k^{\prime}\in[K]_{1},$ (8)

And for all $t\in\mathcal{T},i\in[n]_{1}$, Constraints (1d) can be
reformulated as:

$\displaystyle\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},k,L})+\beta_{k,L}\leq
Mh_{i,k,L},$ (9)
$\displaystyle\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},k,L})+\beta_{k,L}\geq\epsilon+(-M-\epsilon)(1-h_{i,k,L}),$
(10) $\displaystyle
z_{i,k^{\prime},t,L}\leq\alpha_{k^{\prime},t,L}+M(1-h_{i,k^{\prime},L-1}),\
\forall\ k^{\prime}\in[K]_{1},$ (11) $\displaystyle
z_{i,k^{\prime},t,L}\geq\alpha_{k^{\prime},t,L}-M(1-h_{i,k^{\prime},L-1}),\
\forall\ k^{\prime}\in[K]_{1},$ (12) $\displaystyle-Mh_{i,k^{\prime},L-1}\leq
z_{i,k^{\prime},t,L}\leq Mh_{i,k^{\prime},L-1},\ \forall\
k^{\prime}\in[K]_{1}.$ (13)

Constraints (2) and (3) are big-M constraints that, when combined with the
integrality of $h_{i,k,0}$, impose the output of the first hidden layer to be
1 if the linear combination of the units of the input vector summed with the
bias term is greater than or equal to some small constant $\epsilon>0$, and 0
otherwise. To reformulate Constraints (1c), we introduce an auxiliary variable
$z_{i,k^{\prime},k,\ell}=\alpha_{k^{\prime},k,\ell}h_{i,k^{\prime},\ell-1}$.
Constraints (4) and (5) are then similar to Constraints (2) and (3), and
define the output of the remaining hidden layers $\ell\in[L-1]_{1}$. We
leverage the fact that the bi-linear term is a product of a continuous and
binary variable for Constraints (6), (7), and (8); $z_{i,k^{\prime},k,\ell}=0$
when $h_{i,k^{\prime},\ell-1}=0$ or
$z_{i,k^{\prime},k,\ell}=\alpha_{k^{\prime},k,\ell}$ when
$h_{i,k^{\prime},\ell-1}=1$. Similarly for the definition of Constraints
(9)-(13), which define the binary activated output of the output later, and
ensures $z_{i,k^{\prime},t,L}=\alpha_{k^{\prime},t,L}h_{i,k^{\prime},L-1}$,
the bi-linear terms associated with the output layer.

A full proof of this proposition can be found in the Appendix B. However, we
will present a sketch here. The main techniques for this proof rely on first
using a big-M formulation for disjunctive constraints [see, for example 41,
42] to model the binary activation. Then bi-linear terms are reformulated
using the techniques applied to products of binary and continuous variables.

###  3.2 Loss Functions for Prescription

In this section we design MIP losses for ANNs in the context of prescription
problems.

The key challenge of using observational data is the lack of counterfactual
information. Without the counterfactual, it would be impossible to identify
the best treatment in expectation for a given individual. Thus, we consider
this problem as a conditionally randomized experiment, by making the following
key assumptions [43].

######  Assumption 1.

For every treatment $t\in\mathcal{T}$, the outcome $Y(t)$ is independent of
observed treatment, $T$, given the covariates $X$. That is, conditional
exchangeability, $Y(t)\perp\\!\\!\\!\\!\perp T|X$, holds.

Exchangeability in this context refers to the case where for any two
treatments $t,t^{\prime}\in\mathcal{T}$, the group that receive a treatment
$t$ would observe the same outcome distribution as the group that receive a
different treatment $t^{\prime}\neq t$ had their treatments been exchanged.
Here, we specifically assume that conditional exchangeability
$Y(t)\perp\\!\\!\\!\\!\perp T|X$ holds as, a sufficient assumption to
conceptualize observational studies as conditionally randomized experiments.
The key for this assumption to hold is that the treatment groups are
conditionally exchangeable on $X$ since we consider only information from $X$
to prescribe treatments. Indeed this is a sufficiency condition. That is, the
counterfactual outcome and the observed treatments are independent conditioned
on the covariates since only the patient information is used to assign the
treatment.

######  Assumption 2.

For every $t\in\mathcal{T}$, the probability that an individual $i\in[n]_{1}$
with covariates $X_{i}$ is assigned $t$ is strictly positive. That is
positivity, $\mathbb{P}(t_{i}=t|X_{i})>0$, holds.

Both assumptions are standard in literature and are known to hold in practice.
For example, in the domain of personalized medicine, Assumption 1 holds
considering that only the patient’s covariates are available in an
observational study. Particularly, patients features are the only information
used for most medical guidelines used prescribes a treatment for their
patient. Assumption 2 holds considering that practitioners often use intuition
developed over clinical or medical training to often prescribe treatments that
deviate from treatment protocol. That is, patients need not be prescribed the
same treatment despite having the same characteristics. Moreover, the
assumption is tenable given that practitioners may disagree on treatment
plans.

We now discuss methods from the literature that estimate the average causal
effect of observational studies that satisfy the above conditions. The first
is the Inverse Propensity Weighting (IPW) method, which was introduced to
estimate the expected value of counterfactual treatment outcomes and average
treatment effect [29]. The method creates a pseudo-population where every
individual in the data is given all of the possible treatments, simulating a
randomized controlled trial. The IPW method first estimates the propensity
scores for each individual $i$, $\mathbb{\hat{P}}[T=t_{i}|X=X_{i}]$, and then
re-weights their outcome by taking its inverse. Given a policy $s$, we
leverage the conditional exchangeability assumption to estimate the
performance of $s$ as a measure of the average outcome:

$\displaystyle\hat{\pi}^{IPW}(s):=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbbm{1}(s(X_{i})=t_{i})}{\mathbb{\hat{P}}[T=t_{i}|X=X_{i}]}Y_{i}$
(14)

If the estimated propensity score converges almost surely to the true
propensity score, the IPW estimator is unbiased but suffers from high variance
for cases where $\mathbb{\hat{P}}[T=t_{i}|X=X_{i}]$ is small [43].

The second method we are interested in is the so-called Direct Method (DM) to
directly estimate the average causal effect of observational studies. DM is an
adaption of the Regress-and-Compare (RC) approach for estimating
counterfactual outcomes for optimal treatment assignment. By the RC framework,
we first partition the dataset by treatment, then for $X\in\mathcal{X}$ learn
a model $\hat{\mu}_{t}(X)$ for $\mathbb{E}[Y|X,T=t]$ using the subpopulation
that was assigned treatment $t$, and finally estimate the counterfactual
outcomes $\hat{Y}(t)$ as $\hat{\mu}_{t}(X)$. In the case where we are
maximizing over the outcomes, the RC approach would then take
$\operatorname*{arg\,max}_{t\in T}\hat{\mu}_{t}(X)$ as the prescription for an
individual [44]. The DM approach modifies this approach to estimate average
outcome of a $s$:

$\displaystyle\hat{\pi}^{DM}(s):=\frac{1}{n}\sum_{i=1}^{n}\hat{\mu}_{s(X_{i})}(X_{i})$
(15)

However, we note that the performance of RC relies heavily on the quality of
the estimate $\hat{\mu}_{t}(X),t\in T$ [45]. As is the case with IPW, if the
estimated outcome $\hat{\mu}_{t}(X)$ is unbiased, then the DM estimator is
unbiased as well.

Lastly, we discuss the Doubly Robust (DR) estimation method as proposed by
[32] and adapted by [21]. To counteract the reliance on one estimator being
unbiased, this method combines the IPW and DM approaches to reduce the errors
individually brought in by the two estimators; the DR estimator of the outcome
appropriately combines the estimates $\hat{\mu}_{t}(X)$ and
$\mathbb{\hat{P}}[T=t_{i}|X=X_{i}]$. It can then estimate an individual’s
outcome using the doubly robust framework as:

$\displaystyle\hat{\psi}_{t}(z_{i})=\hat{\mu}_{s(X_{i})}(X_{i})+\frac{\mathbbm{1}(s_{0}(X_{i})=t_{i})}{\mathbb{\hat{P}}[T=t_{i}|X=X_{i}]}(Y_{i}-\hat{\mu}_{t_{i}}(X_{i}))$
(16)

where $z_{i}=(x_{i},t_{i},Y_{i})$ is the tuple of the individual’s covariates,
treatment, and outcome from the data, and we let
$\mathcal{Z}:=\mathcal{X}\times\mathcal{T}\times\mathcal{Y}$. Under
Assumptions 1 and 2, the bias of the doubly robust estimator is small if
$\hat{\mu}(X)$ is close to $\mu(X)$ or $\hat{\mathbb{P}}[T=t_{i}|X=X_{i}]$ is
close $\mathbb{P}[T=t_{i}|X=X_{i}]$. Moreover, the bias
$\mathbb{E}[\mathbb{E}[\hat{Y}_{DR}]-\mathbb{E}[Y]]$ of the estimated outcomes
$\mathbb{E}[\hat{Y}]$ is asymptotically zero when either the IPW model or the
DM model is consistent. This doubly robust estimator benefits from second-
order bias, so estimating $\hat{\mu}(X)$ and
$\hat{\mathbb{P}}[T=t_{i}|X=X_{i}]$ with machine learning estimators can lead
to a smaller bias than standard parametric models [43].

Given the doubly robust estimator $\hat{\psi}_{t}$, [37] learn the optimal
policy by solving the optimization problem $\operatorname*{arg\,max}_{t\in
T}\hat{\psi}_{t}(z)$ in a maximization problem. Alternately, one may be
interested in optimizing over the treatment effects rather than the estimated
outcomes directly. [33] present a one such approach for evaluating a policy
$s^{\prime}$ against a baseline policy $s_{0}$. In the dichotomous treatment
case, the function $\hat{\pi}(s^{\prime},s_{0})$ can be used to evaluate
$s^{\prime}$ against the baseline $s_{0}$:

$\displaystyle\hat{\pi}(s^{\prime},s_{0})$
$\displaystyle=\mathbb{E}_{n}\big{[}[s^{\prime}(X_{i})-s_{0}(X_{i})]\hat{\psi}_{1}(z_{i})-[s^{\prime}(X_{i})-s_{0}(X_{i})]\hat{\psi}_{0}(z_{i})\big{]}$
(17)
$\displaystyle=\mathbb{E}_{n}[(s^{\prime}(X_{i})-s_{0}(X_{i}))(\hat{\psi}_{1}(z_{i})-\hat{\psi}_{0}(z_{i})]$
(18)
$\displaystyle=\mathbb{E}_{n}[(s^{\prime}(X_{i})-s_{0}(X_{i}))\hat{\tau}(X_{i})]$
(19)

where
$\hat{\tau}(X_{i})=\mathbb{E}_{n}[\hat{\psi}_{1}(z_{i})-\hat{\psi}_{0}(z_{i})]$
is the estimator for the average treatment effect. For cases where we consider
more than two treatments, we can then define the following:

$\displaystyle\hat{\pi}(s)$ $\displaystyle=\mathbb{E}_{n}[\sum_{t\in
T}\mathbbm{1}\\{s(X_{i})=t\\}\hat{\psi}_{t}(z_{i})]$ (20)
$\displaystyle\hat{\pi}(s^{\prime},s_{0})$
$\displaystyle=\mathbb{E}_{n}[\sum_{t\in
T}\mathbbm{1}\\{s^{\prime}(X_{i})=t\\}\hat{\psi}_{t}(z_{i})-\sum_{t\in
T}\mathbbm{1}\\{s_{0}(X_{i})=t\\}\hat{\psi}_{t}(z_{i})]$ (21)

By using $\hat{\pi}$ as our desired estimator, we explicitly consider the
average treatment while evaluating the performance of a policy. Such a
formulation also allows us to directly optimize against a known baseline
policy. We now propose a linear objective function derived from Equations 20
and 21 through Proposition 2:

######  Proposition 2.

Under Assumption 1 and 2, and in the case where larger outcomes are preferred,
the objective of the optimal prescription problem that optimizes over the CATE
with DR estimators can be formulated as

$\displaystyle\max\frac{1}{n}\sum_{i=1}^{n}\sum_{t\in
T}h_{i,t,L}\hat{\psi}(z_{i})_{t}$ (22)

when we do not have a baseline policy, and

$\displaystyle\max\frac{1}{n}\sum_{i=1}^{n}\sum_{t\in
T}(h_{i,t,L}-s_{0}(x_{i})_{t})\hat{\psi}(z_{i})_{t}$ (23)

when we have a baseline policy $s_{0}$.

The proof of Proposition 2 follow directly from the derivations detailed in
[33]. The proof can be found in the Appendix B for the sake of completeness.

Clearly, Proposition 2 can be adapted to formulate the objective with the IPW
and DM estimators as well. From our findings in Section 3.1 along with the
above reformulated objectives, our main result extends to the formulation of
Prescriptive Neural Networks. Matching the notation of this proposition to
that presented in (1) we can think of the loss function in this case as
$\mathcal{L}_{n}(s(X_{i}),Y_{i})=-\sum_{t\in\mathcal{T}}h_{i,t,L}\hat{\psi}_{n}(z_{i})_{t}$

We note the need to include one additional constraint to the MIP formulation
that necessitates only one treatment to be assigned for each individual:

$\displaystyle\sum_{t\in\mathcal{T}}h_{i,t,L}=1,\ \forall\ i\in[n]_{1}$ (24)

###  3.3 Consistency of Prescriptive Networks

An advantage of policies calculated using our MIP formulation is that they can
be shown to produce consistent prescriptive policies. In general, to show
consistency of estimates derived from solving an optimization problem one must
show the convergence of the minimizers of a sequence of stochastic
optimization problems to the minimizer of a target optimization problem. This
is not trivial since it requires showing that the value functions of these
optimization problems properly converge to the value function of the target
problem. For instance point-wise convergence of these functions is
insufficient to guarantee this property [46]. Thus our focus will be to show
that in the case of MIP trained prescriptive networks, we obtain a necessary
form of convergence that yields consistency. To obtain this property we will
have to make one simplifying assumption on the parameters of the input layer
of the network.

######  Assumption 3.

Parameters $\alpha_{d,k,0}$ come from a finite set $\mathcal{A}$.

Note that this assumption is only required on the parameters of the input
layer and not any of the subsequent layers. Moreover, this assumption is not
too restrictive since we are already assuming all parameters come from a
compact set and thus we can think of $\mathcal{A}$ as an arbitrarily fine grid
of small width. In our computational experiments, we did not place this
integrality restriction on our parameters and found no reduction in
performance. In addition to this assumption we need to make another standard
assumption when considering our causal inference based loss.

######  Assumption 4.

Estimators $\hat{\mu}_{n},\mathbbm{\hat{P}}_{n},\hat{\psi}_{n}$ are
consistent. That is, as $n\rightarrow\infty$,
$\hat{\mu}_{n}\overset{p}{\rightarrow}\mu,\mathbbm{\hat{P}}\overset{p}{\rightarrow}\mathbbm{P},\hat{\psi}_{n}\overset{p}{\rightarrow}\psi$,
where $\mu,P,\psi$ are the true values of the direct, indirect, and doubly
unbiased estimators as evaluated with the population parameters.

This assumption is standard when analyzing prescriptive analytics [37, 35].
Essentially we are assuming that the methods employed for counterfactual
estimation are statistically consistent, since if this is not the case then it
would be impossible to learn a consistent policy.

With these assumptions in mind we will define our target optimization problem
as
$\Theta^{*}=\operatorname*{arg\,min}_{\theta\in\Theta}\mathbb{E}[\sum_{t\in\mathcal{T}}\mathbbm{1}\\{s(X,\theta)=t\\}\psi(Z)]$,
where $Z\in\mathcal{Z}$, and the expectation is taken with respect to the true
probability law of the population, which we will call $\mathbb{P}_{0}$. Note
that, $\Theta^{*}\subset\Theta$ is a set and not a single parameter since
there could be multiple paramaeterizations for the optimal network structure.
Without loss of generality, we are considering the case where counterfactual
estimation is given by the doubly robust estimate, though in principle our
arguments will also hold for the direct and indirect estimates. Our main
result is showing the following theorem:

######  Theorem 3.

Given Assumptions 1-4 and all parameter and modeling assumptions. For all
$\epsilon>0$, $\hat{\theta}_{n}$ as computed by (1) satisfies:

$\lim_{n\rightarrow\infty}\mathbb{P}_{0}(\text{dist}(\hat{\theta}_{n},\Theta^{*})>\epsilon)=0$
(25)

Where for a point $x$ and set $S$ $\text{dist}(x,S):=\min_{y\in S}\|x-y\|$.

To show this main result we will need to prove several intermediate results on
the structure of the value functions of the target optimization problem and
the training problem. For notational brevity let
$\varphi(\theta,Z)=\sum_{t\in\mathcal{T}}\mathbbm{1}\\{s(X,\theta)=t\\}\psi(Z)$
and
$\varphi_{n}(\theta,Z)=\sum_{t\in\mathcal{T}}\mathbbm{1}\\{s(X,\theta)=t\\}\psi_{n}(Z)$.
Our first result will establish the continuity of these functions.

######  Proposition 4.

Under Assumption 3, the functions $\varphi(\theta,Z),\varphi_{n}(\theta,Z)$
are lower semicontinuous for all pairs $(\theta,Z)\in\Theta\times\mathcal{Z}$.

The full proof of this result can be found in the appendix. However, this
proposition will follow from the construction of $\varphi,\varphi_{n}$ as
value functions of parametric MIP problems with the parameters appearing as
affine terms in the constraints [47]. Lower semi continuity is a crucial
property since it ensures the closure of the epigraph of the value function,
and is one of the conditions necessary to establish proper convergence of
minima. Next we establish a point-wise condition for the convergence of the
empirical loss function to the target loss function.

######  Proposition 5.

Given Assumptions 1-4. For any $\epsilon>0$ define
$U_{\epsilon}(\theta):=\\{\theta^{\prime}\in\Theta:\|\theta^{\prime}-\theta\|<\epsilon\\}$.
Then consider the following quantities:

$\displaystyle
H_{n}^{\epsilon}(\theta,Z)=\frac{1}{n}\sum_{i=1}^{n}\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},z_{i})$ (26) $\displaystyle
H^{\epsilon}(\theta)=\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi(\theta^{\prime},Z)]$ (27)

Then for all $\theta\in\Theta,\epsilon>0$, $H^{\epsilon}(\theta)>-\infty$ and
for any $n\geq 1$ both $H_{n}^{\epsilon}(\theta,Z),H^{\epsilon}(\theta)$ exist
and are well defined. Moreover,
$H_{n}^{\epsilon}(\theta,\cdot)\overset{l-prob}{\rightarrow}H^{\epsilon}(\theta)$,
where the relation $H_{n}\overset{l-prob}{\rightarrow}H$ implies the sequence
$H_{n}$ is a lower semicontinuous approximation in probability to the target
$H$ as per Definition 4.1 in [48].

This proposition establishes a point-wise functional law of large numbers for
our MIP ANN models. Note that this result is stronger than a standard law of
large numbers since we are required to control the convergence behavior of our
value functions within an open neighborhood of the parameters. While a full
proof of the proposition will be in the appendix here we will provide a brief
sketch. Since we are assuming the data $Z$ is i.i.d, the two key arguments
needed to establish this relation are showing that the quantities in the
proposition satisfy Kolmogorov’s Law of Large Numbers, and then that a strong
enough convergence is achieved to apply Slutsky’s Theorem [49]. Note that this
is still a point-wise result since it holds for each individual $\theta$, and
it is not by itself sufficient to prove our main result. Thus we require one
additional intermediate result before proceeding to prove our main theorem.

######  Proposition 6.

Given Assumptions 1-4. Let $\Phi,\Phi_{n}:\theta\mapsto\mathbb{R}$ be defined
as
$\Phi(\bar{\theta}):=\min_{\theta\in\Theta}\big{\\{}\mathbb{E}[\varphi(\theta,Z)]:\theta=\bar{\theta}\big{\\}}$
and
$\Phi_{n}(\bar{\theta}):=\min_{\theta\in\Theta}\big{\\{}\frac{1}{n}\sum_{i=1}^{n}\varphi_{n}(\theta,z_{i}):\theta=\bar{\theta}\big{\\}}$
respectively. Then
$\Phi_{n}\underset{\Theta}{\overset{l-prob}{\rightarrow}}\Phi$, where the
operator indicates that $\Phi_{n}$ are a lower semi-continuous approximation
in probability to $\Phi$ across the entire parameter set $\Theta$.

Note that this is a stronger global result across the entire parameter set
$\Theta$ than the point-wise result in Proposition 5. Moreover, note that the
objective $\frac{1}{n}\sum_{i=1}^{n}\varphi_{n}(\theta,z_{i})$ is equivalent
to the objective function in (1a) meaning $\Phi(\theta)$ can be thought of as
the value function of the training problem and $\Phi$ as the value function of
the target problem. This proposition implies that as we collect additional
data, we have the proper form of epi-convergence between our empirical and
traget problems. The full proof of this proposition is in the appendix, and
relies on using the results of Proposition 5 in conjunction with existing
results from the convergence of stochastic optimization problems [48].
Combining this final result with properties of convergence of minima of
stochastic optimization problems yields the main result of Theorem 3.

##  4 Numerical Experiments on Prescriptive Problems

In this section, we evaluate the performance of the PNN in numerical
experiments. We use experimental data from three settings: a synthetic dataset
commonly used in literature introduced by [21], real-world data from a study
on postpartum hypertension [50], and a dataset for personalized warfarin
dosing collected by the International Warfarin Pharmacogentic Consortium [51].
We benchmark our methods against four state of the art models. These include
prescriptive tree models from [35], [36], and [37] along with the causal
forest approach from [30]. For brevity we refer to these methods as B-PT,
K-PT, J-PT, and CF respectively.

Prescriptive trees are mixed-integer optimization frameworks that learn
policies for prescription problems that can be encoded as a binary decision
trees. The notation used to describe each model follow from [37]. A limitation
shared by all prescriptive tree models is the need for binary valued features
as input parameters in the optimization framework. We provide a description of
the prescriptive tree methods below.

K-PT: This model introduced by [35] is one of the first MIP models used to
learn prescriptive trees. The formulation relies on producing trees by
recursively partitioning the covariates such that no feature falls under two
partitions. Formally, given that $\mathcal{L}$ represents the set of leaf
nodes and $\mathcal{X}_{l}$ denote the subset of covariate values for the data
points collected in leaf $l\in\mathcal{L}$, the model creates partitions
$\mathcal{X}=\cup_{l\in\mathcal{L}}\mathcal{X}_{l}$ such that
$\mathcal{X}_{l}\cap\mathcal{X}_{l^{\prime}}=\emptyset,\ l\neq l^{\prime}$. A
key assumption made to support the theoretical guarantees of the K-PT
formulation is that depth of the trees are sufficiently deep such that the
treatment assignment in the data and the covariates are independent at each
leaf node. With respect to the experiments detailed in this section, our
analysis is conducted on depth-1 and depth-2 trees, considering that the
optimality gap of the K-PT model scales poorly with size of the tree (see
Appendix D.1).

B-PT: [36] address some limitations of the K-PT model and provide a modified
MIP objective that regularizes against high variance in the outcomes of the
data points conditioned on treatment assigned in each leaf. In particular, the
authors include the term
$\mu\sum_{i=1}^{n}{\big{(}Y_{i}-\hat{Y}_{\mathcal{X}_{l(i)}:l(i)\in\mathcal{L}}(t_{i})\big{)}}^{2}$
to the K-PT objective where $\mu\geq 0$ is an appropriately tuned
hyperparameter. In our experiments, we found that the hyperparamter $\mu=0.5$
yielded strong results for the B-PT formulation, a consistent result with
experiments conducted in [36]. Again, our model is benchmarked against depth-1
and depth-2 B-PT models since these models scale poorly with deeper trres.

J-PT: To address the need for the assumptions on the depth of the perscriptive
tree, [37] directly optimize the outcomes estimated with the IPW, DM or DR
methods. The formulation proposed in their work follows as an adaptation of
the weighted classification tree formulation introduced by [38] for the
prescriptive setting. The J-PT is modelled as a flow graph with edges between
a source node and the root of a binary tree, and the leaves of the tree and a
sink node. We analyze experimental results for depth 1 and depth 2 trees to
maintain consistency with the other prescriptive trees.

We use the CF as another model to benchmark our work in our numerical
experiments. A description of the CF is given below.

CF: We benchmark against an adapted causal forest model proposed by [30]. The
original work estimates heterogeneous treatment effects using nonparametric
random forests. CF estimates the CATE by recursively splitting the feature
space such that training examples are close with respect to their treatments
and outcomes, where the measure of closeness is identified by data points
falling in the same leaf of the tree. Since CF does not explicitly prescribe
treatments, we adapt the original method to use the sign of the estimated
treatment effect for each individual to assign treatments. In experiments with
multiple treatment groups, we construct a pipeline to estimate the pairwise
individual treatment effect and assign those treatments which yield the
largest outcome for each individual. Experiments on the CF model use the
EconML package [52] with their default arguments, which builds 100 trees in
the forest.

All experiments were run on a high throughput computing cluster systems
allocated with 4 CPU Cores, 32GB of RAM, and 4GB of Disk Space, with Python
version 3.9.5 [53] and Gurobi Optimizer version 9.1.2 [54].

###  4.1 Experimental Setup

We discuss the nature of the datasets used in our experiments in this
subsection. Since J-PT and PNN optimize over doubly robust objectives, we also
provide insight into estimating the propensity scores and counterfactual
outcomes for each experiment.

####  4.1.1 Simulated Data

The experiments with the simulated data were adapted from [21] for data
simulation, and [37] for the experiment design. Given that we simulate the
counterfactual outcomes, a synthetic dataset allows us to evaluate learned
policies against the ground truth. We considered three models of varying
feature sizes $\mathcal{F}$. As given by Equation 28, the potential outcome
for each individual in the simulation
$Y_{i}:\mathbbm{R}^{\mathcal{F}}\times\\{0,1\\}\mapsto\mathbbm{R}$ is a
function of their covariates
$X_{i}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{1})\in\mathbbm{R}^{\mathcal{F}}$,
where $\boldsymbol{0}$ and $\boldsymbol{1}$ are $\mathcal{F}$-dimensional
vectors, and their assigned treatment $t_{i}\in\\{0,1\\}$:

$\displaystyle
Y_{i}(t)=\eta(X_{i})+\frac{1}{2}\cdot(2t-1)\cdot\kappa(X_{i})+\epsilon_{i}$
(28)

We define the mean effect by
$\eta:\mathbbm{R}^{\mathcal{F}_{1}\subseteq\mathcal{F}}\mapsto\mathbbm{R}$ and
the treatment effect by
$\kappa:\mathbbm{R}^{\mathcal{F}_{2}\subseteq\mathcal{F}}\mapsto\mathbbm{R}$,
and introduce a noise term $\epsilon_{i}\sim\mathcal{N}(0,0.01)$ independent
of $X_{i}$. The three model designs we consider are as follows:

  1. 1.

$\displaystyle\mathcal{F}=2;\;\eta(x)=\frac{1}{2}x_{1}+x_{2};\;\kappa(x)=\frac{1}{2}x_{1}$

  2. 2.

$\displaystyle\mathcal{F}=10;\;\eta(x)=\frac{1}{2}\sum_{f=1}^{2}x_{f}+\sum_{f=3}^{6}x_{f};\;\kappa(x)=\sum_{f=1}^{2}\mathbbm{1}\\{x_{f}>0\\}\cdot
x_{f}$

  3. 3.

$\displaystyle\mathcal{F}=20;\;\eta(x)=\frac{1}{2}\sum_{f=1}^{4}x_{f}+\sum_{f=5}^{8}x_{f};\;\kappa(x)=\sum_{f=1}^{4}\mathbbm{1}\\{x_{f}>0\\}\cdot
x_{f}$

Note that $\eta$ and $\kappa$ are functions mapped from subsets of features.
Furthermore, in designs 2 and 3, some features are not present in either
function. Such a setup should help us assess the models’ ability to train over
only those features that are relevant to the patient outcome. In designs 2 and
3, the treatment effect $\kappa(x)$ is always non-negative, and so, not
accounting for the noise, treatment $t_{i}=1$ always yields an outcome that is
equal to or better than the outcome if treatment $t_{i}=0$. The purpose of
such a design is to test the regularization capabilities of complex learning
against skewed data. Note that simpler models, such as prescriptive trees, are
incentivized to design policies where every individual receives treatment
$t_{i}=1$.

To simulate conditionally randomized treatments, we use varying parameter
$p\in\\{0.1,0.25,0.5,0.75,0.9\\}$ for the probability of correct treatment
assignment. Marginally randomized treatments are simulated at $p=0.5$. We
generate 5 random datasets with a train:test split of 100:10,000 datapoints,
and 5 with a train:test split of 500:10,000 datapoints for each $p$. We
evaluate each model by taking the average Out-of-Sample Performance (OOSP) –
the percentage of correct treatment assignments – over the 5 datasets. The
prescriptive tree models used as benchmarks necessitate binary feature for a
feature-based branching rule. Therefore, we transform the simulated real and
continuous covariates into binary-valued features by discretizing them into
deciles from a normal distribution and obtaining their one-hot encoded
vectors. The one-hot encoded vectors were further transformed by converting
$x_{i(j+1)}$ to $1$’s if $x_{ij}=1$ for all $j\in\\{1,..,10\\}$ and for all
$x_{i}\in\mathcal{X}$. We refer to these transformed covariates as adapted
binarized covariates. The adapated binarized covariates were used for all
three experimental designs for the prescriptive trees. To address the value
that adapted binarization adds to the prescriptive tree formulations, we run
one additional set of experiments on experimental design 1 with one-hot
binarized covariates.

The PNN and J-PT models necessitate counterfactual outcome estimation. We
utilize a doubly robust estimation model, where we consider linear regression
and a lasso regression model to fit the counterfactual outcomes directly, and
a decision tree and a logistic regression model to estimate the propensity
scores. The linear regression and decision tree models were found to be the
best pair of models via a model selection regime that uses 10-fold cross
validation for all simulated data and compares classification accuracies for
each pair. The parameters for the linear regression and decision tree models
were passed as arguments to the [30] causal forest model, which estimates the
propensity score and counterfactual data internally. We also used 10-fold
cross validation to tune $\ell_{0}$ and $\ell_{1}$ regularization
hyperparameters for the PNN models.

####  4.1.2 Personalized Postpartum Hypertension Treatments

We use prescriptive models to design treatment policies that minimize the
maximum postpartum systolic blood pressure (SBP) at a patient-level.
Hypertensive disorders are one of the leading causes of morbidity and
mortality worldwide [55]. Approximately one third of morbidity and mortality
occurs during the first six weeks postpartum. Partially, this is because
pregnancies complicated by hypertensive disorders will often have exacerbation
of hypertension in the postpartum period that is likely secondary to a
combination of mobilization of extravascular fluid as well as a persistence of
preeclampsia [55]. This can lead to elevations in blood pressure that are
often described in the first 3 to 7 days postpartum [56]. Thus evaluating how
different treatment policies effect SBP during this time period is a problem
of significant importance to the obstetrics community. The data comes from a
retrospective cohort study from a single United States Midwestern academic
center of all women who delivered from 2009-2018. The original data set
includes 32,645 deliveries, where 170 were readmitted in the postpartum period
due to a hypertension-related diagnosis [57]. This data set is an
observational study that includes 67 features collected from individuals
during specific milestones – baseline/pre-pregnancy through the prenatal
period (first, second, and third trimesters, pre-admission period and the day
of admission for delivery) and up to 42 days postpartum. Of the 67 features, a
subset of 18 were chosen as covariates (12 binary, 3 categorical, and 3
continuous) after eliminating correlated features by consulting with domain
experts. We considered two treatment conditions – whether or not the patient
was prescribed blood pressure medication, and which day the patient was
discharged postpartum. We take a patient’s maximum SBP postpartum to be a
function of their covariates conditioned upon the treatment that they
received. We designed three experiments with this dataset, all with the
objective to estimate policies that minimize the maximum SBP and differ by the
treatment conditions: (1) Blood pressure medication as a binary treatment
condition, (2) Discharge day postpartum as an ordinal treatment condition, (3)
Both blood pressure medication and discharge day postpartum.

Figure 1: Causal Graphs Detailing the Covariates, Treatments, and Outcomes
Along with their Causal Structures for the Three Hypertension Experiments.

Figure 1 depicts the characteristics of the three experiments in the form of
causal graphs. For experiments 1 and 2, we considered the alternate treatment
condition as a covariate (discharge day postpartum and blood pressure
medications, respectively). We conducted experiments on a subset of the
complete dataset to account for missing and incomplete data. Additionally, we
discarded a single outlier – a patient who was not prescribed any blood
pressure medication but was discharged on the 5th day postpartum – to address
the case where experiment 3 sees only one patient with this treatment regime.
The final dataset included 291 individuals from whom we had at least 38 blood
pressure measurements. The discharge day postpartum for every individual was
transformed into 4 ordinal categories: {day 1 or day 2, day 3, day 4, day 5+}.
We performed 10-fold cross validation to tune the $\ell_{0}$ regularization
hyperparameter for the PNN models. Finally, for the decision tree models, the
continuous covariates in the dataset were discretized and transformed into
adapted binarized covariates, and the categorical covariates were split into
one-hot encoded binary vectors.

We considered a logistic regression model with elastic net regression, a
random forest regressor, and the multilayer perceptron regressor to estimate
the counterfactual outcomes for these experiments and evaluated their
performance using their root-mean-squared errors. To estimate the propensity
scores, we fit a logistic regression model, decision tree classifier, and
random forest classifier and evaluated them using their classification
accuracies. We used 10-fold cross validation to conclude that the logistic
regression models produced the best estimators for experiment types 1 and 2,
with the random forest regressor yielding better estimators was better for
some cross-treatments in experiment 3. We evaluate each model using the
average performance of the learned policy $\pi(s)$ over each experiment run
with 10-fold cross-validation.

###  4.2 Results

In this section, we analyze the results of the experiments conducted on
simulated and real-world data. We build a PNN of depth 2 (one hidden layer,
and one output layer) for every experiment. PNNs of widths 3 and 10 were
benchmarked against the prescriptive tree formulations of [36] (B-PT), [35]
(K-PT), and [37] (J-PT), all with depths 1 and 2. Experiments on the Causal
Forest (CF) model use the EconML package [52] with their default arguments,
which builds 100 trees in the forest. All approaches were allotted a run time
of 1 hour.

####  4.2.1 Simulated Data

A summary of the results of the three experimental designs for the simulated
dataset can be found in Figures 2 and 3. The complete table of results can be
found in Appendix D.1. Results for the hyperparameter tuning can be found in
Appendix D.2. It is evident from these figures that the PNN models match or
outperform the OOSP of its benchmarks in every experiment. Additionally, their
performance is unperturbed by the probability of correct treatment assignment
in the data. Furthermore, the PNN achieves a strong OOSP given a 100-sample
training set as well as a 500-sample training set in the allocated time.

Figure 2: Results of the simulated dataset with experimental design 1. Panel A
visualizes the results when the data is binarized using adapted encoding.
Panel B visualizes the results when the data is binarized using one-hot
encoded binarization. The left plots are policies learned from 100 datapoints
and the right plots are policies learned from 500 datapoints.

Panel A of Figure 2 shows the experimental results of design 1 when the
covariates have the adapted binarization for the prescriptive trees. Panel B
of Figure 2 depicts the results for the same design, but with one-hot encoded
binarization. These results show that the performance of prescriptive trees
are sensitive to different featurization of the data. While the best
prescriptive tree model achieves an OOSP around 90% across all treatment
assignments when the covariates are binarized to support feature splitting
with the adapated binarization, the best prescriptive tree only yields a
policy with a 70% OOSP when discretized into deciles without the
featurization. These experiments also show that the performance of the causal
forest model is dependant on to the probability of correct treatment
assignment in the data, as the OOSP for this model tapers off at the 0.1 and
0.9 probabilities. This means the causal forest designs policies that are
stronger when the treatments in the data are marginally randomized, but are
weaker with more conditionally randomized treatments. Additionally, the B-PT
and K-PT formulations drop below a 50% OOSP when the probability of correct
treatment assignment approaches 0.9. That is, when the informed policy assigns
the correct treatment to most patients, these trees design a policy that are
akin to a random treatment assignment. Given that these models achieve a 0%
optimality gap for these datasets, this drop is performance is likely
attributed to a failure in branching on key features at the given depths.

Figure 3: Results of the simulated dataset with experimental design 2 and 3.
Panel A visualizes the results when the data is binarized using one-hot
encoding with the left learned from 100 datapoints and the right learned from
500 datapoints. Panel B visualizes the results when the data is binarized
using adapted binarization with the left learned from 100 datapoints and the
right learned from 500 datapoints.

Figure 3 shows the experimental results of designs 2 and 3 on Panels A and B
respectively, with the adapted binarization for the prescriptive trees. Recall
that for these designs, it is easy to achieve a strong OOSP from a policy that
assigns treatments liberally. In fact, the policies designed by the PNN and
the causal forest are identical in each experiment for designs 2 and 3, with
an average OOSP of around 85% and 95% respectively achieved by assigning
treatments to everyone. Unsurprisingly, the depth 1 counterpart for each
prescriptive tree model achieves a higher peak OOSP given the structure of the
datasets. We see that the J-PT model with appropriate depths achieve the same
policy as the PNNs and the causal forest when trained on a larger training
set. However, there exists a larger variance in the performance of the B-PT
and K-PT models, implying limitations in identifying key branching features
that support a treatment-forward policy. Tables 3 and 4 in Appendix D.1 show
that the policies designed by the B-PT and K-PT models have a standard
deviation in the OOSP of up to 36.55% and 16.16% respectively. The drop in
performance of the B-PT model in the 500-sample training set plot of Panel B
in Figure 3 is attributed to the fact that the policy was designed by the
heuristic solution found by the solver; an optimal solution was not found in
the allocated time limit.

####  4.2.2 Personalized Postpartum Hypertension Treatments

The experiments for the hypertension datasets involved designing policies that
reduce the maximum systolic blood pressure (SBP) of postpartum patients in the
study. Since the counterfactual outcomes for the patients are unknown, it is
difficult to evaluate and compare various policies. Hence, we computed the
expected CATE of each policy. Note that we estimate the average CATE for the
policy in the data using the doubly robust estimator to allow for an adjusted
comparison of the learned policies against the policy in data. we used 10-fold
cross validation to tune the $\ell_{0}$ regularization hyperparameters for the
PNN. Results for the hyperparameter tuning can be found in Appendix D.2.

Figure 4: Results of experiments on the hypertension dataset for out-of-sample
patients prescribed treatments based on policies learned by different models.
Each panel displays the mean out-of-sample maximum systolic blood pressure.
The left plots include the 90% confidence intervals and the right plots
include the 95% confidence intervals.

Figure 4 provides a summary of the results including the confidence intervals
for the three experiments. The complete experimental results can be found in
Appendix D.1. Recall that the experiments differ in terms of how the treatment
condition is defined: Experiment 1 considers whether or not a patient receives
blood pressure medication, Experiment 2 considers the day a patient is
discharged postpartum, and Experiment 3 considers the cross-product of the
two. We observe that the PNN model outperforms the decision tree and causal
forest benchmarks with regards to the mean out-of-sample maximum SBP in every
experiment while maintaining a low variance.

Figure 5: Out-of-sample treatment assignments for policies learned by
different models over the hypertension dataset. Each panel displays the number
of patients who were prescribed each treatment.

To further analyze the behavior of each model, we refer to Figure 5. We
observe that when a model has a smaller architecture, the policies they design
are extremely treatment averse in experiment 1 (Panel A). We posit that the
complexity of these models is not high enough to capture the causal effect in
the data. Note that the depth 2 decision trees design policies with a larger
average maximum SBP than the depth 1 trees. This anomaly could be attributed
to the fact that the K-PT and J-PT formulations do not regularize against
overfitting. Furthermore, while the B-PT formulation does incorporate a
regularization term, the depth 2 variant exhibits a large average optimality
gap at the time limit. We also note that the 10-width PNN was able to
outperform all of the other models by designing a treatment policy that is
more willing to prescribe the blood pressure medication.

Panel B of Figure 4 indicates that policies consider the discharge day as the
treatment for postpartum patients lead to a larger mean maximum SBP across all
models. The B-PT and K-PT trees perform worse than the policy in data, and
reveal a large variance, despite the B-PT formulation incorporating
regularization. This may be explained by the large optimality gap in the depth
2 tree, and the notion that the depth 1 was unable to regularize successfully
with 4 treatment groups. Panel B of Figure 5 additionally shows that the
causal forest model recommends a highly risk-averse regiment, with the
majority of the patients being discharged on or after day 4. Additionally, the
J-PT trees are more inclined to discharge patients at day 3. Considering the
performance of the J-PT trees and the PNNs, it can be interpreted that the
majority of the patients benefit from being discharged on day 3 and day 4.

Panels C of Figures 4 and 5 indicate the results of experiment 3. The general
trends of experiments 1 and 2 are seen to be present in experiment 3 as well,
with the majority of smaller architecture models designing policies that are
treatment-averse. Note that the 10-width PNN only suggests a treatment plan
where patients receive no blood pressure medication, and are discharged either
on days 1, 2 or 3. While this may seem irregular, the policy yielded the
lowest average maximum SBP. We also see that the J-PT trees maintain their
affinity for discharging patients on day 3. We hypothesize that the lack of
regularization explains the decline in performance of depth 2 trees in
comparison with their depth 1 counterparts, and that the B-PT depth 2 tree
failed to achieve a small optimality gap to design a stronger policy.

Figure 6: Paired results of experiments on the hypertension dataset for out-
of-sample patients prescribed treatments based on policies learned by
different models. Each panel displays paired evaluation of the various models.
The left heatmaps display the point differences of out-of-sample mean maximum
systolic blood pressure between pairs of models. Differences are taken as row
- column. The right plots display the hypothesis tests on the differences of
out-of-sample mean maximum systolic blood pressure between pairs of models.
Smaller p-values favor models in the rows. Cool colors denote lower blood
pressures (left) or smaller p-values (right). Axes values denote the model (PT
depth/# of trees in CF/PNN width).

Additionally, we conducted a paired hypothesis test to make statistically
significant claims on the performance of each model. The results of the
hypothesis test are presented in Figure 6, which aims to provide a
visualization that allows for a fair comparison of the models. The left
subplots show the difference in the mean out-of-sample maximum systolic blood
pressure (SBP) between each pair. The right subplots show the p-values of the
paired hypothesis tests given the null hypothesis
$N_{0}:\mu_{row}-\mu_{column}=0$, where $\mu_{row}$ is the true mean of the
model on the rows and $\mu_{columns}$ is the true mean of the model on the
columns. PNNs design policies with a mean out-of-sample maximum SBP that is
statistically significantly lower than its competitors’ policies at least one
of the chosen architecture. We note that the large variance with the B-PT and
K-PT models are translated here as relatively large p-values. Overall, our
results show that PNNs statistically outperform other causal and prescriptive
models in this setting.

###  4.3 Model Interpretation

High stakes domains such as personalized medicine in the context of assigning
treatments for postpartum patients demand policies that are fair and can be
audited for individual safety. Therefore, interpretability and accountability
are highly desirable. In this section, we assess each model’s ability to
discern important features in the personalized postpartum hypertension
treatment assignment.

Analyzing the tree-based approaches is straightforward – features that are
deemed more important are branched on higher in the tree. Given that we
experiment on trees with a depth of at most 2, these models should therefore
choose at most 3 features to branch on. Single-layer PNNs are similarly
straightforward to review – features consider more important have larger
absolute weights, and less important have weights closer to 0 along the arcs
of the network. The inclusion of the $\ell_{0}$ regularization term in the
objective further motivates sparsity in the weights, making the model easier
to interpret. However, the causal forest is much less straightforward to
interpret. In order to analyze the policies designed by all the models on a
level field, we turn to Shapley Additive Explanations, or SHAP values [58]. We
compute the global SHAP values for the features in each model to interpret the
prescriptive trees, causal forests PNNs.

To assess the policies designed by the 5 models in the personalized postpartum
hypertension treatment experiments, we ranked the features by their SHAP
values for each model. Since these experiments followed a 10-fold cross-
validation regime, we observed 10 models with different policies. To gain a
better understanding of these policies in general, the feature rankings were
then averaged across these 10 models. Furthermore, since the prescriptive
trees branch on at most 3 features, those that were not chosen to branch on
were given the lowest ranking. Recall now that the prescriptive trees optimize
on a binarized feature space. To fairly compare the PNN and CF models with the
tree models, we took the minimum rank for each binarized feature as the
reported rank. Finally, to account for varied number of features due the
binarization, the rankings were normalized between 0 and 1. We then subtract
these normalized rankings from 1 to maintain an intuitive representation of
these rankings. So, features with a larger normalized rank are considered more
important by the model. We note that rows with null values represent those
models where only treatment was assigned in each fold of the cross-validation;
due to this behavior, the results of the SHAP model are inconclusive for these
cases.

We now analyze the models based on their feature rankings. Tables 16-18 in
Appendix D.3 display these rankings for the three experiments. We see that the
features that are commonly ranked highly in the PNNs are prenatal BMI, chronic
hypertension, mode of delivery, gestational hypertension, gestational age, and
preeclampsia. Additionally, in experiments 1 and 2, discharge day postpartum
and blood pressure medication are highly ranked, respectively. These
observations align well with medical and clinical preferences [59, 50]. We
also note that social characteristics such as insurance type is often ranked
low by the PNNs. We see that the causal forest ranks similar features highly,
with the birthing person’s age often considered most important. We also
observe that the B-PT and K-PT models frequently consider features such as
mode of delivery and prenatal BMI as important. However, the J-PT model finds
it more challenging to identify medically relevant features in experiment 1
and 2. We note that ethnicity and insurance are among the highest ranked
feature by the depth 2 model in both experiments, which may be a concerning
since these attributes are generally considered protected. Note that the J-PT
model appropriately identifies important features in experiment 3. This may be
because correlations brought on by including the alternate treatments as
covariates are negated.

##  5 Managerial Implications

From our case studies, our implementation of MIP-based ANNs does not only
provide strong empirical performance, but also provides additional insights
into how our framework can be deployed in healthcare and business settings.

  1. 1.

The PNN framework yields consistently strong polices when given limited data
and optimized over small networks. The policies designed by the PNN models
used in the experiments consistently achieve the highest out-of-sample
probability of correct treatment assignment in the simulated experiments.
Moreover, when used in the real-world case study of treating patients with
hypertension during pregnancy and postpartum, the PNN models result in the
lowest average maximum systolic blood pressure postpartum when considering any
of the proposed treatment conditions. In particular, the width 3 and width 10
models were consistently the two best performing models, with the best policy
outperforming the policy in the data by 7.62 mmHG (see Figure 6) on average,
in a statistically significant manner. When compared against other
prescriptive models, the width-10 PNN designed policies that resulted in lower
point average maximum systolic blood pressure in every case. Moreover, results
for the join treatment plan in experiment 3 showed a statistically significant
improvement, with a $p$-value less than $0.05$, against all but two models:
the depth-1 J-PT ($p$-value=$0.14$) and CF with 100 trees ($p$-value=$0.11$).
Additionally, the ability to supplement the PNNs with a regularization term
leads to lower policy variance and smaller confidence intervals across
replications. The PNN models also have a better capacity to personalize
treatment as opposed to some other models which tended to provide one
treatment to a majority of patients (see Panels B and C in Figure 5).

  2. 2.

Methods that rely on supplementary models for estimating counterfactual
outcomes are sensitive to accuracy of those models. Our experiments on
warfarin dosing suggest that accuracy of the inverse propensity weighting
method and direct method for estimating counterfactual outcomes significantly
affect the performance of the policies of the models that optimize using the
doubly robust estimator. In particular, given the warfarin dataset, the
logistic regression model, which was deemed the best propensity score
estimator via cross-validation, had an average accuracy of 66.77%. Similarly,
the elastic net model, which was the best direct estimator, yielded an average
root-mean squared error of 2.056. Given that the policy in the data was
ascertained to be optimal, our results show that the models that optimize over
the estimated counterfactual outcomes designed sub-optimal policies. It is
therefore key for decision-makers to to consider such limitations when using
such prescriptive models with any dataset as they be highly sensitive to the
performance of the estimation models.

  3. 3.

The PNN is more likely to identify medically reliable features then other
prescriptive methods. The need to design interpretable policies to solve
complex problems gives rise to a tradeoff for prescriptive models. When we
have strong counterfactual estimations, a smaller neural network can provide
better interpretability on the policies that are designed. Additionally, given
that MIP models can be used solved to small optimality gaps when used to
trained networks with smaller architectures, these interpretable policies also
yield optimal outcomes. The PNN therefore achieves a strong tradeoff since it
a complex yet parsimonious training model. As seen in Tables 16-18, PNNs
identify key patient characteristics, such as prenatal BMI and chronic and
gestational hypertension as confirmed by our clinical collaborator, and
appropriately weigh these features in designing its policy. These results were
consistent across all of the experimental designs for these experiments. In
comparison, we observe that some benchmarked models consider fail to identify
medically reasonable features and instead consider social features such as
ethnicity and insurance in their treatment assignments. The interpretability
of small PNNs allows for medical professionals and policy-makers to work with
this learning models to consider the best approach for assigning treatments.

##  6 Conclusions

Our work presents the effectiveness of 0-1 neural networks in addressing the
challenges of personalized decision-making, particularly in healthcare
domains. By using mixed-integer programming and counterfactual estimation,
prescriptive networks offer a valuable tool for optimizing treatment policies
when limited patient data is available. The ability to encode complex policies
while remaining interpretable sets the PNN apart from existing approaches that
employ decision trees. The results obtained from both synthetic data
experiments and the real-world case study of postpartum hypertension treatment
assignment provide compelling evidence of the advantages offered by
prescriptive neural networks towards improving healthcare outcomes.

## References

  *   * Wiseman [1972] Wiseman, F.: Methodological bias in public opinion surveys. The Public Opinion Quarterly 36(1), 105–108 (1972)
  * Rimando et al. [2015] Rimando, M., Brace, A.M., Namageyo-Funa, A., Parr, T.L., Sealy, D.-A., Davis, T.L., Martinez, L.M., Christiana, R.W.: Data collection challenges and recommendations for early career researchers. The Qualitative Report 20(12), 2025–2036 (2015)
  * Jadad [1998] Jadad, A.R.: Randomised controlled trials: a user’s guide. Health Technology Assessment 2(13), 214 (1998)
  * Polubriaginof et al. [2017] Polubriaginof, F., Salmasian, H., Albert, D.A., Vawdrey, D.K.: Challenges with collecting smoking status in electronic health records. In: AMIA Annual Symposium Proceedings, vol. 2017, p. 1392 (2017). American Medical Informatics Association
  * Gilbert and Mintz [2019] Gilbert, T.K., Mintz, Y.: Epistemic therapy for bias in automated decision-making. In: AIES, pp. 61–67 (2019)
  * Forman and Cohen [2004] Forman, G., Cohen, I.: Learning from little: Comparison of classifiers given little training. In: PKDD, pp. 161–172 (2004). Springer
  * Aminian and Ameri [2005] Aminian, K., Ameri, S.: Application of artificial neural networks for reservoir characterization with limited data. Journal of Petroleum Science and Engineering 49(3-4), 212–222 (2005)
  * Shaikhina and Khovanova [2017] Shaikhina, T., Khovanova, N.A.: Handling limited datasets with neural networks in medical applications: A small-data approach. Artificial intelligence in medicine 75, 51–63 (2017)
  * Chakraborty et al. [2017] Chakraborty, S., Tomsett, R., Raghavendra, R., Harborne, D., Alzantot, M., Cerutti, F., Srivastava, M., Preece, A., Julier, S., Rao, R.M., et al.: Interpretability of deep learning models: A survey of results. In: IEEE smartworld/SCALCOM/UIC/ATC/CBDcom/IOP/SCI), pp. 1–6 (2017). IEEE
  * Fan et al. [2021] Fan, F.-L., Xiong, J., Li, M., Wang, G.: On interpretability of artificial neural networks: A survey. IEEE Transactions on Radiation and Plasma Medical Sciences 5(6), 741–760 (2021)
  * Zhang et al. [2021] Zhang, Y., Tiňo, P., Leonardis, A., Tang, K.: A survey on neural network interpretability. IEEE Transactions on Emerging Topics in Computational Intelligence 5(5), 726–742 (2021)
  * Toro Icarte et al. [2019] Toro Icarte, R., Illanes, L., Castro, M.P., Cire, A.A., McIlraith, S.A., Beck, J.C.: Training binarized neural networks using mip and cp. In: International Conference on Principles and Practice of Constraint Programming, pp. 401–417 (2019). Springer
  * Kurtz and Bah [2021] Kurtz, J., Bah, B.: Efficient and robust mixed-integer optimization methods for training binarized deep neural networks. arXiv preprint arXiv:2110.11382 (2021)
  * Thorbjarnarson and Yorke-Smith [2023] Thorbjarnarson, T., Yorke-Smith, N.: Optimal training of integer-valued neural networks with mixed integer programming. Plos one 18(2), 0261029 (2023)
  * Li et al. [2010] Li, L., Chu, W., Langford, J., Schapire, R.E.: A contextual-bandit approach to personalized news article recommendation. In: WWW, pp. 661–670 (2010)
  * Hitsch and Misra [2018] Hitsch, G.J., Misra, S.: Heterogeneous treatment effects and optimal targeting policy evaluation. Available at SSRN 3111957 (2018)
  * Flume et al. [2007] Flume, P.A., O’Sullivan, B.P., Robinson, K.A., Goss, C.H., Mogayzel Jr, P.J., Willey-Courand, D.B., Bujan, J., Finder, J., Lester, M., Quittell, L., et al.: Cystic fibrosis pulmonary guidelines: chronic medications for maintenance of lung health. American journal of respiratory and critical care medicine 176(10), 957–969 (2007)
  * Mintz et al. [2020] Mintz, Y., Aswani, A., Kaminsky, P., Flowers, E., Fukuoka, Y.: Nonstationary bandits with habituation and recovery dynamics. Operations Research 68(5), 1493–1516 (2020)
  * Mintz et al. [2023] Mintz, Y., Aswani, A., Kaminsky, P., Flowers, E., Fukuoka, Y.: Behavioral analytics for myopic agents. European Journal of Operational Research 310(2), 793–811 (2023)
  * He and Mintz [2023] He, Q., Mintz, Y.: Model based reinforcement learning for personalized heparin dosing. arXiv preprint arXiv:2304.10000 (2023)
  * Athey and Imbens [2016] Athey, S., Imbens, G.: Recursive partitioning for heterogeneous causal effects. Proceedings of the National Academy of Sciences 113(27), 7353–7360 (2016)
  * Bridle [1990] Bridle, J.S.: Probabilistic interpretation of feedforward classification network outputs, with relationships to statistical pattern recognition. In: Neurocomputing: Algorithms, Architectures and Applications, pp. 227–236 (1990). Springer
  * Serra et al. [2018] Serra, T., Tjandraatmadja, C., Ramalingam, S.: Bounding and counting linear regions of deep neural networks. In: ICML, pp. 4558–4566 (2018). PMLR
  * Danna et al. [2007] Danna, E., Fenelon, M., Gu, Z., Wunderling, R.: Generating multiple solutions for mixed integer programming problems. In: IPCO, pp. 280–294 (2007). Springer
  * Fischetti and Jo [2017] Fischetti, M., Jo, J.: Deep neural networks as 0-1 mixed integer linear programs: A feasibility study. arXiv preprint arXiv:1712.06174 (2017)
  * Anderson et al. [2020] Anderson, R., Huchette, J., Ma, W., Tjandraatmadja, C., Vielma, J.P.: Strong mixed-integer programming formulations for trained neural networks. MPA, 1–37 (2020)
  * Tjeng et al. [2017] Tjeng, V., Xiao, K., Tedrake, R.: Evaluating robustness of neural networks with mixed integer programming. arXiv preprint arXiv:1711.07356 (2017)
  * Nair et al. [2020] Nair, V., Bartunov, S., Gimeno, F., Glehn, I., Lichocki, P., Lobov, I., O’Donoghue, B., Sonnerat, N., Tjandraatmadja, C., Wang, P., et al.: Solving mixed integer programs using neural networks. arXiv preprint arXiv:2012.13349 (2020)
  * Horvitz and Thompson [1952] Horvitz, D.G., Thompson, D.J.: A generalization of sampling without replacement from a finite universe. JASA 47(260), 663–685 (1952)
  * Wager and Athey [2018] Wager, S., Athey, S.: Estimation and inference of heterogeneous treatment effects using random forests. JASA 113(523), 1228–1242 (2018)
  * Powers et al. [2018] Powers, S., Qian, J., Jung, K., Schuler, A., Shah, N.H., Hastie, T., Tibshirani, R.: Some methods for heterogeneous treatment effect estimation in high dimensions. Statistics in medicine 37(11), 1767–1787 (2018)
  * Dudík et al. [2011] Dudík, M., Langford, J., Li, L.: Doubly robust policy evaluation and learning. arXiv preprint arXiv:1103.4601 (2011)
  * Farrell et al. [2021] Farrell, M.H., Liang, T., Misra, S.: Deep neural networks for estimation and inference. Econometrica 89(1), 181–213 (2021)
  * Breiman et al. [1984] Breiman, L., Friedman, J., Olshen, R., Stone, C.: Cart. Classification and regression trees (1984)
  * Kallus [2017] Kallus, N.: Recursive partitioning for personalization using observational data. In: ICML, pp. 1789–1798 (2017). PMLR
  * Bertsimas et al. [2019] Bertsimas, D., Dunn, J., Mundru, N.: Optimal prescriptive trees. INFORMS Journal on Optimization 1(2), 164–183 (2019)
  * Jo et al. [2021] Jo, N., Aghaei, S., Gómez, A., Vayanos, P.: Learning optimal prescriptive trees from observational data. arXiv preprint arXiv:2108.13628 (2021)
  * Aghaei et al. [2021] Aghaei, S., Gómez, A., Vayanos, P.: Strong optimal classification trees. arXiv preprint arXiv:2103.15965 (2021)
  * Amram et al. [2022] Amram, M., Dunn, J., Zhuo, Y.D.: Optimal policy trees. ML 111(7), 2741–2768 (2022)
  * Bertsimas et al. [2022] Bertsimas, D., Klasnja, P., Murphy, S., Na, L.: Data-driven interpretable policy construction for personalized mobile health. In: 2022 IEEE ICDH, pp. 13–22 (2022). IEEE
  * Wolsey and Nemhauser [1999] Wolsey, L.A., Nemhauser, G.L.: Integer and Combinatorial Optimization vol. 55. John Wiley & Sons, ??? (1999)
  * Conforti et al. [2014] Conforti, M., Cornuéjols, G., Zambelli, G., et al.: Integer Programming vol. 271. Springer, ??? (2014)
  * Hernán and Robins [2010] Hernán, M.A., Robins, J.M.: Causal inference. CRC Boca Raton, FL (2010)
  * Qian and Murphy [2011] Qian, M., Murphy, S.A.: Performance guarantees for individualized treatment rules. Annals of statistics 39(2), 1180 (2011)
  * Beygelzimer and Langford [2009] Beygelzimer, A., Langford, J.: The offset tree for learning with partial labels. In: ACM SIGKDD, pp. 129–138 (2009)
  * Rockafellar and Wets [2009] Rockafellar, R.T., Wets, R.J.-B.: Variational Analysis vol. 317. Springer, ??? (2009)
  * Ralphs and Hassanzadeh [2014] Ralphs, T.K., Hassanzadeh, A.: On the value function of a mixed integer linear optimization problem and an algorithm for its construction. COR@ L Technical Report 14T–004 (2014)
  * Vogel and Lachout [2003] Vogel, S., Lachout, P.: On continuous convergence and epi-convergence of random functions. part i: Theory and relations. Kybernetika 39(1), 75–98 (2003)
  * Keener [2010] Keener, R.W.: Theoretical Statistics: Topics for a Core Course. Springer, ??? (2010)
  * Mukhtarova et al. [2021] Mukhtarova, N., Alagoz, O., Chen, Y.-H., Hoppe, K.: Evaluation of different blood pressure assessment strategies and cutoff values to predict postpartum hypertension-related readmissions: a retrospective cohort study. American Journal of Obstetrics & Gynecology MFM 3(1), 100252 (2021)
  * Consortium [2009] Consortium, I.W.P.: Estimation of the warfarin dose with clinical and pharmacogenetic data. NEJM 360(8), 753–764 (2009)
  * Battocchi et al. [2019] Battocchi, K., Dillon, E., Hei, M., Lewis, G., Oka, P., Oprescu, M., Syrgkanis, V.: EconML: A Python Package for ML-Based Heterogeneous Treatment Effects Estimation. https://github.com/microsoft/EconML. Version 0.13.1 (2019)
  * Van Rossum and Drake [2009] Van Rossum, G., Drake, F.L.: Python 3 Reference Manual. CreateSpace, Scotts Valley, CA (2009)
  * Gurobi Optimization, LLC [2023] Gurobi Optimization, LLC: Gurobi Optimizer Reference Manual (2023). https://www.gurobi.com
  * Garovic et al. [2022] Garovic, V.D., Dechend, R., Easterling, T., Karumanchi, S.A., McMurtry Baird, S., Magee, L.A., Rana, S., Vermunt, J.V., August, P.: Hypertension in pregnancy: diagnosis, blood pressure goals, and pharmacotherapy: a scientific statement from the american heart association. Hypertension 79(2), 21–41 (2022)
  * Hirshberg et al. [2016] Hirshberg, A., Levine, L., Srinivas, S.: Clinical factors associated with readmission for postpartum hypertension in women with pregnancy-related hypertension: a nested case control study. Journal of Perinatology 36(5), 405–409 (2016)
  * Hoppe et al. [2020] Hoppe, K.K., Thomas, N., Zernick, M., Zella, J.B., Havighurst, T., Kim, K., Williams, M., Niu, B., Lohr, A., Johnson, H.M.: Telehealth with remote blood pressure monitoring compared with standard care for postpartum hypertension. American Journal of Obstetrics & Gynecology 223(4), 585–588 (2020)
  * Lundberg and Lee [2017] Lundberg, S.M., Lee, S.-I.: A unified approach to interpreting model predictions. NeurIPS 30 (2017)
  * Hoffman et al. [2021] Hoffman, M.K., Ma, N., Roberts, A.: A machine learning algorithm for predicting maternal readmission for hypertensive disorders of pregnancy. American Journal of Obstetrics & Gynecology MFM 3(1), 100250 (2021)
  * Shalev-Shwartz and Ben-David [2014] Shalev-Shwartz, S., Ben-David, S.: Understanding Machine Learning: From Theory to Algorithms. Cambridge university press, ??? (2014)
  * Bertsimas and Tsitsiklis [1997] Bertsimas, D., Tsitsiklis, J.N.: Introduction to Linear Optimization vol. 6. Athena Scientific Belmont, MA, ??? (1997)
  * Goodfellow et al. [2017] Goodfellow, I., Bengio, Y., Courville, A.: Deep learning (adaptive computation and machine learning series). Cambridge Massachusetts, 321–359 (2017)
  * Calafiore and El Ghaoui [2014] Calafiore, G.C., El Ghaoui, L.: Optimization Models. Cambridge university press, ??? (2014)
  * Cormen et al. [2022] Cormen, T.H., Leiserson, C.E., Rivest, R.L., Stein, C.: Introduction to Algorithms. MIT press, ??? (2022)
  * Zhang and Sabuncu [2018] Zhang, Z., Sabuncu, M.: Generalized cross entropy loss for training deep neural networks with noisy labels. NeurIPS 31 (2018)
  * LeCun [1998] LeCun, Y.: The mnist database of handwritten digits. http://yann. lecun. com/exdb/mnist/ (1998)
  * Ravi and Larochelle [2017] Ravi, S., Larochelle, H.: Optimization as a model for few-shot learning. In: ICLR (2017)
  * Berge [1963] Berge, C.: Espaces Topologiques: Fonctions Multivoques. Dunod, Paris, 1959. English translation: Translated by EM Patterson. Topological Spaces. New York: MacMillan (1963)
  * Walker [1979] Walker, M.: A generalization of the maximum theorem. International Economic Review, 267–272 (1979)
  * Vogel and Lachout [2003] Vogel, S., Lachout, P.: On continuous convergence and epi-convergence of random functions. part ii: Sufficient conditions and applications. Kybernetika 39(1), 99–118 (2003)
  * Dumont et al. [2013] Dumont, Z., Mordasiewicz, M., Kosar, L., Schuster, B.: Warfarin: Its highs and lows. Canadian Family Physician 59(8), 856–860 (2013)

##  Appendix A MIP ANNs and Prediction

###  A.1 ANNs for Prediction: Negative Log Likelihood Loss

A common loss for prediction is the negative log likelihood loss (NLL), in
this section we will present a MIP formulation for this loss in the case of
classification.

Our goal is to represent Equation 1a in terms of a set of linear constraints
and linear objective function such that our ANN makes predictions on data
given disjoint classes. Note that the $L^{th}$ layer outputs are the only ones
that interact directly with the loss function, so we focus on these decision
variables and leave the discussion of the model parameters when considering
individual unit activations. The loss function proposed by [12] is motivated
by maximizing the neuron margin, akin to Support Vector Machines. However,
max-margin models are known to under-perform on noisy data with overlapping
classes [60]. Other prediction objectives that have been proposed in the
literature are $\ell_{0}$ or $\ell_{1}$ losses that can be straightforwardly
formulated through MILP linearization techniques [61, 41]. However, in
practice a common loss function used for training ANNs in prediction tasks, is
the negative log likelihood applied to a soft-max distribution [62]. Such a
loss can be considered analogous to minimizing the log likelihood of a
multinomial distribution across the classes in $J$, thus outputting a
predictive distribution over the classes. This loss is particularly appealing
as it harshly penalizes incorrect predictions and is numerically stable [22,
62]. The equation for the soft-max activation function for the output of the
final layer $h_{i,j,L}$ is
$\varphi(h_{i,j,L})=\frac{\exp(h_{i,j,L})}{\sum_{j^{\prime}\in
J}\exp(h_{i,j^{\prime},L})}$, and the resulting negative log likelihood loss
is $-\log(\prod_{i=1}^{n}(\varphi(h_{i,j,L})^{\mathbbm{1}[Y_{i}=j]}))$. We
reformulate this loss using the following result:

######  Proposition 7.

The negative log likelihood loss applied to the soft-max activation
$-\log(\prod_{i=1}^{n}(\varphi(h_{i,j,L})^{\mathbbm{1}[Y_{i}=j]}))$ is within
a constant additive error of the optimal value of the following linear
optimization problem:

$\displaystyle\min_{h_{i,j,L},\omega_{i}}\big{\\{}\sum_{i=1}^{n}\sum_{j\in
J}\mathbbm{1}[Y_{i}=j](\omega_{i}-h_{i,j,L}):\omega_{i}\geq
h_{i,j^{\prime},L},\ j^{\prime}\in[J]_{1},i\in[n]_{1}\big{\\}}$ (29)

Here we present a brief sketch of the proof, for the detailed proof of this
proposition see Appendix B. The main insight used is the fact that the
function $\log\sum_{i=1}^{n}\exp(X_{i})=\Theta(\max_{i\in[n]_{1}}X_{i})$ [63],
where $\Theta$ is the Big Theta notion that characterizes the function
$\log\sum\exp$ to be bounded asymptotically both above and below by the $max$
function [64], and the fact that the minimization of a maximum can be written
using linear constraints [41].

We note that for the output of the final layer to give a unique prediction
with the negative log likelihood soft-max loss, we need add constraints which
force the optimization problem to assign diverse values to
$h_{i,j,L},\forall\;i\in[n]_{1},j\in J$:

$\displaystyle\omega_{i}\geq h_{i,j,L}$ (30) $\displaystyle
h_{i,j,L}+h_{i,j^{\prime},L}-2h_{i,j,L}\leq-\epsilon+Mr_{i,j,j^{\prime}},\
\forall\ j,j^{\prime}\in[J]^{2}_{1},j\neq j^{\prime}$ (31) $\displaystyle
h_{i,j,L}+h_{i,j^{\prime},L}-2h_{i,j,L}\geq\epsilon-M(1-r_{i,j,j^{\prime}}),\
\forall\ j,j^{\prime}\in[J]^{2}_{1},j\neq j^{\prime}$ (32)

These constraints can be thought of as diversifying constraints, since without
them the loss would be minimized by setting all the outputs of the final
layer, $h_{i,j,L}$ to be equal. However, owing to the fact that no data point
can be labeled with multiple classes, we require unique unit outputs.
Constraints (30) forces $\omega_{i}$ to be the max of the components of the
output vector. Finally, Constraints (31) and (32) ensure that no two units in
the output layer are equal. The idea behind these two constraints relies on
the fact that the average of two variables is equal to one of the variables if
and only if the two variables are equal. So, we introduce the binary variables
$r_{i,j,j^{\prime}}\in\mathbbm{B}$ to define Constraints (31) and (32) as
big-M constraints.

Considering that the soft-max activation of the output layer is captured by
the objective function, we reformulate Constraint (1d) as:

$\displaystyle
h_{i,j,L}\leq\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},j,L})+\beta_{j,L}$ (33)
$\displaystyle
h_{i,j,L}\geq\sum_{k^{\prime}=1}^{K}(z_{i,k^{\prime},j,L})+\beta_{j,L}$ (34)

Constraints (33) and (34) guarantee that the output of the final layer is the
linear combination of the units of the activated outputs of the penultimate
layer summed with the bias term.

###  A.2 Numerical Experiments on Prediction Problems

We now evaluate the MIP formulation for training ANNs through empirical
experiments for prediction problems. We aim to show value of incorporating a
negative log likelihood loss into the MIP formulation. For the prediction
problems, we benchmark against the Stochastic Gradient Descent algorithm with
momentum and adaptive learning rates with both binary activations and ReLU in
the hidden layer. The SGD models were fit with Cross Entropy Loss which
incorporates a negative log likelihood loss over the sigmoid activation in the
output layer [65]. Experiments were done over a subset of the MNIST dataset
[66] in the few-shot learning set-up [67]. All experiments were run on a high
throughput computing cluster systems allocated with 4 CPU Cores, 32GB of RAM,
and 4GB of Disk Space, with Python version 3.9.5 [53].

####  A.2.1 Experimental Setup

We use the MNIST handwritten dataset [66] to evaluate our model’s performance
in the image recognition problem. The images in this dataset include labelled
handwritten digits from 0 through 9 in a grey-scale $28\times 28$ pixels
resolution. We design the experiments as a few-shot learning problem by taking
a subset of the labels and images for training; in these experiments, we
consider images labeled as $2$ or $7$, and sampled $100$ images for training
and $1000$ images for testing. Sampling was conducting with 5 seeds and
without replacement between the training and testing sets to evaluate the
average out-of-sample testing accuracy of each model. As a result, each model
was given $100$ reshaped images of size $28\times 28=784$ and was asked to
learn to predict dichotomous labels.

The baseline SGD with momentum models were trained on $50,000$ iterations, an
initial learning rate $=0.1$, a constant learning rate decay $=0.5$ every
$10,000$ steps, an initial momentum $=0.5$, and final momentum $=0.95$. These
hyperparameters were fine tuned through cross-validation for both binary and
ReLU activation models. Finally, we used cross-validation to fine tune the
$\ell_{0}$ regularization term in the objective of the MIP ANN model. The
results of the cross-validation experiments showed that not having a
regularization term yields the best performance for this experimental design.

We assess the performance of the MIP, SGD-Binary and SGD-ReLU ANNs over a
series of architectures, reported in the form (number of hidden layers, number
of neurons): (1,3), (1,5), (1,10), (2,3), (2,5), (2,10), (3,3), (3,5), (3,10).
The results of such a neural network architecture search are reported in this
section to better understand the settings in which a MIP-based training regime
serves as viable alternative to SGD.

####  A.2.2 Results

We analyze the results of the prediction problem experiments conducted over
the 2-Digit MNIST datasets in this section. All approaches were allotted a run
time of 1 hour.

Figure 7: Average Final training (left) and testing (right) accuracies of the
three prediction ANNs over the 2-Digit MNIST dataset.

Figure 7 shows a summary of the results averaged over the 5 seeds. The full
experimental results can be found in Appendix D.1. We observe that the MIP-
based training approach outperforms the SGD training approaches in all but one
ANN architectures when considering the average final testing accuracy. On
average, the MIP-based ANN with 3 hidden layers with 3 neurons in each hidden
layer yielded the best final testing accuracy at 90%. For the 3 hidden layers
and 10 neurons architectures, we found that the MIP model was unable to solve
to optimality in the allocated time limit, due to which we associate the
deterioration in performs. Additionally, we see that when we use the ReLU
activation function in the hidden layers for the SGD approach, the model fails
to perform better than random. In fact, our experiments showed that the
performs of the ReLU model improves when we train over deeper neural networks;
considering the few-shot learning setting, however, the binary models were
able to yield results better than random classification with much more
parsimonious neural networks. We find that these results are consistent with
that which is observed by [12, 13, 14].

##  Appendix B Proofs

###  B.1 Proof of Proposition 1

###### Proof.

There are two key challenges for the reformulation of Constraints
(1b),(1c),(1d) first is the binary activation function itself and second are
the bi-linear terms present in the hidden and output layers. Without loss of
generality, let us consider a particular unit in the neural network such that
$h_{n,k,\ell}=\mathbbm{1}[\sum_{k^{\prime}=0}^{K}\alpha_{k,k^{\prime}}h_{k^{\prime},\ell}+\beta_{k,\ell}\geq
0]$. First let us consider the case of the binary activation. This condition
can be modeled as a disjunction [41, 42], a form of constraint that can be
reformulated as follows:

$\displaystyle\sum_{k^{\prime}=0}^{K}\alpha_{k,k^{\prime}}h_{k^{\prime},\ell}+\beta_{k,\ell}\leq
Mh_{n,k,\ell}$ (35)
$\displaystyle\sum_{k^{\prime}=0}^{K}\alpha_{k,k^{\prime}}h_{k^{\prime},\ell}+\beta_{k,\ell}\geq\epsilon+(-M-\epsilon)(1-h_{n,k,\ell})$
(36)

Where $M$ is a sufficiently large constant, and $\epsilon$ is a small
constant. Next let us consider the bi-linear terms
$\alpha_{k^{\prime},k,\ell}h_{k^{\prime},\ell-1}$. Here we note that
$\alpha_{k^{\prime},k,\ell}\in[\alpha^{L},\alpha^{U}]$ is a continuous valued
variable, while $h_{k^{\prime},\ell-1}\in\mathbbm{B}$ is binary valued. This
type of product can be reformulated though a standard technique by introducing
an additional binary variable
$z_{n,k^{\prime},k,\ell}=\alpha_{k^{\prime},k,\ell}h_{k^{\prime},\ell-1}$ then
these products can be written as:

$\displaystyle
z_{n,k^{\prime},k,\ell}\leq\alpha_{k^{\prime},k,\ell}+M(1-h_{n,k^{\prime},\ell-1})$
(37) $\displaystyle
z_{n,k^{\prime},k,\ell}\geq\alpha_{k^{\prime},k,\ell}-M(1-h_{n,k^{\prime},\ell-1})$
(38) $\displaystyle-Mh_{n,k^{\prime},\ell-1}\leq z_{n,k^{\prime},k,\ell}\leq
Mh_{n,k^{\prime},\ell-1}$ (39)

where again $M$ is an appropriately picked sufficiently large constant. Using
these reformulations where appropriate and making the necessary variable
substitutions results in the formulations presented above. ∎

###  B.2 Proof of Proposition 7

###### Proof.

First, we can explicitly write the negative log likelihood as:

$\displaystyle-\log(\prod_{n=1}^{N}(\sigma(h_{n,j,L})^{\mathbbm{1}[y_{n}=j]}))=-\sum_{n=1}^{N}\mathbbm{1}[y_{n}=j]\log\Big{(}\frac{\exp(h_{n,j,L})}{\sum_{j^{\prime}\in
J}\exp(h_{n,j^{\prime},L})}\Big{)}$ (40)
$\displaystyle=-\sum_{i=n}^{N}\mathbbm{1}[y_{n}=j](h_{n,j.L}-\log\sum_{j^{\prime}\in
J}\exp(h_{n,j^{\prime},L}))$ (41)

Here we note that the expression
$\log\sum_{j^{\prime}=J}\exp(h_{n,j^{\prime},L})=\Theta(\max_{j^{\prime}\in
j}h_{n,j^{\prime},L})$ [63]. Thus using this fact we obtain:

$\displaystyle\eqref{NLL 2}$
$\displaystyle\leq-\sum_{n=1}^{N}\mathbbm{1}[y_{n}=j](h_{n,j,L}-\max_{j^{\prime}\in
J}\\{h_{n,j^{\prime},L}\\}-\log(|J|))$ (42)
$\displaystyle=\sum_{n=1}^{N}\mathbbm{1}[y_{n}=j](\max_{j^{\prime}\in
J}\\{h_{n,j^{\prime},L}\\}-h_{n,j,L}+\log(|J|))$ (43)

Using the other side of the big $\Theta$ condition yields that $\eqref{NLL
2}\geq\sum_{n=1}^{N}\mathbbm{1}[y_{n}=j](\max_{j^{\prime}\in
J}\\{h_{n,j^{\prime},L}\\}-h_{n,j,L})$. If we let
$\omega_{n}=\max_{j^{\prime}\in J}\\{h_{n,j^{\prime},L}\\}$ we can use
standard formulation techniques to reformulate it using linear constraints
[41]. In conjunction with removing the constant terms that do not depend on
$\omega_{n},\\{h_{n,j^{\prime},L}\\}$ yields the desired result. ∎

###  B.3 Proof of Proposition 2

###### Proof.

Given that we are in the dichotomous treatment case, [33] gives us the
following estimators for causal parameters for the full-population averages:

$\displaystyle\hat{\tau}$
$\displaystyle=\mathbbm{E}_{n}[\hat{\psi}_{1}(z_{i})-\hat{\psi}_{0}(z_{i})]$
(44) $\displaystyle\hat{\pi}(s)$
$\displaystyle=\mathbbm{E}_{n}[s(X_{i})\hat{\psi}_{1}(z_{i})+(1-s(X_{i}))\hat{\psi}_{0}(z_{i})]$
(45)

We can generalize Equation 45 for more than two treatments by incorporating
the indicator function such that the function returns 1 if the policy outcome
matching the treatment and 0 otherwise:

$\displaystyle\hat{\pi}(s)$ $\displaystyle=\mathbb{E}_{n}[\sum_{t\in
T}\mathbbm{1}\\{s(X_{i})=t\\}\hat{\psi}_{t}(z_{i})]$ (46)

Considering that the output of the MIP-ANN from Proposition 1 with Constraint
24 gives us one-hot encoded vector where the node associated with the assigned
treatment equals 1 and all other nodes are 0, this is precisely the indicator
function we use in Equation 46. Thus, maximizing over this function under the
parameters $z_{i}$ given in the data, we arrive at the desired objective

$\displaystyle\max\frac{1}{n}\sum_{i=1}^{n}\sum_{t\in
T}h_{i,t,L}\hat{\psi}(z_{i})_{t}$ (47)

Where $h_{i,t,L}$ is the $t^{\text{th}}$ node of the output layer for the
$i^{\text{th}}$ individual in the vector output by the ANN.

The proof of the objective given a baseline policy $s_{0}$, given by Equation
23, follows identically as above. ∎

###  B.4 Proof of Proposition 4

###### Proof.

We note that the argument for the lower semi-continuity of
$\varphi_{n}(\theta,Z)$ and $\varphi(\theta,Z)$ is essentially the same so we
will focus on showing $\varphi(\theta,Z)$ is lower semi-continuous.

Using the result of Propositions 1-2, we can see that $\varphi(\theta,Z)$ can
be expressed as the value function of a special case of the formulation with
$n=1$ and $\hat{\psi}$ replaced with $\psi$. In particular, suppose that
$\theta=(\\{\alpha^{\theta}_{d,k,0},\alpha^{\theta}_{k^{\prime},k,\ell},\alpha^{\theta}_{k^{\prime},t,L}\\}_{\forall
k,k^{\prime}\in[K]_{1},\ell\in[L-1]},\\{\beta^{\theta}_{k,\ell}\\}_{\forall
k\in[K]_{1},\ell\in[L]})$, then:

$\displaystyle\varphi(\theta,Z)=\min-\sum_{t\in\mathcal{T}}h_{t,L}\psi(Z)_{t}$
(48a) $\displaystyle\text{subject to: }\eqref{h0_define C2}-\eqref{zL_define
C3}$ (48b) $\displaystyle\sum_{t\in\mathcal{T}}h_{t,L}=1$ (48c)
$\displaystyle\alpha_{d,k,0}=\alpha^{\theta}_{d,k,0},\alpha_{k^{\prime},k,\ell}=\alpha^{\theta}_{k^{\prime},k,\ell},\alpha_{k^{\prime},t,L}=\alpha^{\theta}_{k^{\prime},t,L},$
$\displaystyle\forall k,k^{\prime}\in[K]_{1},\ell\in[L-1]$ (48d)
$\displaystyle\beta_{k,\ell}=\beta^{\theta}_{k,\ell},$ $\displaystyle\forall
k\in[K]_{1},\ell\in[L].$ (48e)

While this formulation is integer linear for fixed $Z$ or $\theta$ note that
the way the input layer is processed by Constraints (3) results in bi-linear
terms from the products $\alpha_{d,k,0}X_{d}$. However using Assumption 3 we
can reformulate these products to remove the bi-linearity. Let us index the
elements of set $\mathcal{A}$ as follows:
$\mathcal{A}=\\{\alpha_{1},\alpha_{2},...,\alpha_{m}\\}$, then let us define a
new set of binary variables $\nu_{i,k,d}\in\mathbb{B}$ and new continuous
variables $z^{0}_{d,k,i},z^{0}_{d,k}\in\mathbb{R}$ for $i\in\\{1,...,m\\}$ and
we can rewrite (3) as:

$\displaystyle\sum_{d\in\mathcal{F}}z^{0}_{d,k}+\beta_{k,0}\leq Mh_{k,0}$
$\displaystyle\forall k\in[K]_{1}$ (49)
$\displaystyle\sum_{d\in\mathcal{F}}z^{0}_{d}+\beta_{k,0}\geq\epsilon+(-M-\epsilon)(1-h_{k,0})$
$\displaystyle\forall k\in[K]_{1}$ (50)
$\displaystyle-M(1-\nu_{i,k,d})+\alpha_{i}x_{d}\leq
z^{0}_{d,k,i}\leq\alpha_{i}x_{d}+M(1-\nu_{i,k,d}),$ $\displaystyle\forall
i\in[m]_{1},d\in\mathcal{F},k\in[K]_{1}$ (51) $\displaystyle-M\nu_{i,k,d}\leq
z^{0}_{d,k,i}\leq M\nu_{i,k,d},$ $\displaystyle\forall
i\in[m]_{1},d\in\mathcal{F},k\in[K]_{1}$ (52)
$\displaystyle\sum_{i=1}^{m}\nu_{i,k,d}=1,$ $\displaystyle\forall
k\in[K]_{1},d\in\mathcal{F}$ (53)
$\displaystyle\sum_{i=1}^{m}z^{0}{d,k,i}=z^{0}_{d,k},$ $\displaystyle\forall
k\in[K]_{1},d\in\mathcal{F}.$ (54)

Finally we can reformulate the first portion of (48d) with respect to
$\nu_{i,k,d}$ as:

$\displaystyle\sum_{i=1}^{m}\alpha_{i}\nu_{i,k,d}=\alpha_{k,d}^{\theta},$
$\displaystyle\ \forall k\in[K]_{1},d\in\mathcal{F}.$ (55)

Since these reformulations imply that the corresponding optimization problem
to $\varphi(\theta,Z)$ depends on $\theta,Z$ only as affine terms in the
constraints and a continuous function in the objective, and that this
optimization problem is a mixed integer linear minimization problem applying
Berge’s Maximum Theorem [68, 69] in conjunction with the results from [47]
implies that $\varphi(\theta,Z)$ is lower semi-continuous on the domain
$\Theta\times\mathcal{Z}$. As previously mentioned a similar reformulation
argument can be used to show that $\varphi_{n}(\theta,Z)$ is also lower semi-
continuous on the domain $\Theta\times\mathcal{Z}$. ∎

###  B.5 Proof of Proposition 5

###### Proof.

To proceed with the proof, note that by construction of the training problem
in (1) (and thus by extension its reformulation as a MIP), the problem will be
feasible for any $\theta\in\Theta$, this means that the first conditions of
the proposition are all satisfied. Thus we proceed to show the convergence
result.

To prove that $H^{\epsilon}(\theta,\cdot)$ is a lower semi-continuous
approximation to $H^{\epsilon}(\theta)$ we will first need to show that for
any given sequence of data, the two quantities converge in probability. Then
we can use the result of Proposition 4 in conjunction with Proposition 5.1
from [48] to turn this initial point-wise condition into a global condition
over the set of all possible data observations.

Thus we proceed to show the point-wise convergence in probability to complete
the proof. Using the formal definition of convergence in probability our goal
is to show that
$|H_{n}^{\epsilon}(\theta,Z)-H^{\epsilon}(\theta)|\overset{p}{\rightarrow}0$.
Since the data observations are i.i.d. and $H_{n}^{\epsilon}$ is constructed
using an arithmetic mean we have the following relationship:

$\mathbb{E}[H_{n}^{\epsilon}(\theta,Z)]=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},z_{i})]=\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},Z)].$ (56)

The last equality follows from the i.i.d. assumption. We can thus upper bound
our target expression using the triangle inequality:

$|H_{n}^{\epsilon}(\theta,Z)-H^{\epsilon}(\theta)|\leq\\\
\underbrace{|H_{n}^{\epsilon}(\theta,Z)-\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},Z)]|}_{\text{(a)}}+\underbrace{|\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},Z)]-H^{\epsilon}(\theta)|}_{\text{(b)}}$
(57)

First consider expression (b), by Assumption 4 this expression must converge
in probability to zero. Thus if we can show that expression (a) also converges
in probability to zero we can apply Slutsky’s theorem to complete the proof
[49].

Observe that by assumption $\forall Z\in\mathcal{Z}$ we have that
$|\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},z)|\leq\text{diam}(\mathcal{Y})$,
where
$\text{diam}(\mathcal{Y})=\max\\{|y-y^{\prime}|:y,y^{\prime}\in\mathcal{Y}\\}$.
Therefore both $\mathbb{E}[\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},Z)]<\infty$ and
$\text{Var}(\inf_{\theta^{\prime}\in
U_{\epsilon}(\theta)}\varphi_{n}(\theta^{\prime},Z))<\infty$, which means we
can apply Kolmogorov’s law of large numbers [49] to obtain that (a) must
indeed converge in probability to zero. ∎

###  B.6 Proof of Proposition 6

###### Proof.

Since the results of Proposition 5 hold $\forall\theta\in\Theta$ and
$\forall\epsilon>0$, we have that by Theorem 3.2 from [70]
$\Phi_{n}\overset{l-prob}{\underset{\\{\theta\\}}{\rightarrow}}\Phi$, that is
it is a point-wise lower semi-continuous approximation in probability
$\forall\theta\in\Theta$. Note that as a direct consequence of Proposition 4
both $\Phi,\Phi_{n}$ are lower semi-continuous functions of $\theta$ since
semi-continuity is preserved under addition and multiplication by positive
scalars. Thus applying Proposition 5.1 from [48] proves the desired result. ∎

###  B.7 Proof of Theorem 3

###### Proof.

Because Proposition 6 implies that the empirical training problem’s value
function is lower semi-continuous approximation in probability to the target
problem, we can apply Theorem 4.3 from [48] to obtain the main result of the
proposition. ∎

##  Appendix C Prescriptive ANN Experiments on Warfarin Dosing

###  C.1 Experimental Setup

We now employ the prescriptive models to learn personalized warfarin dosing.
Warfarin is an oral anticoagulant prescribed to reduce the risk of blood
clots. While its usage is commonplace, prescribing the appropriate dosage can
be difficult because an individual’s response to the anticoagulant varies
greatly based on their genetic, clinical, and demographic characteristics
[51].

We use the publicly available dataset curated by the International Warfarin
Pharmacogenetics Consortium and published by the Pharmacogenetics and
Pharmacogenomics Knowledge Base. The dataset comprises of information for 5700
individuals including single nucleotide polymorphisms (SNPs) in genes for
cytochrome P450 2C9 (CYP2C9) and vitamin K epoxide reductase (VKORC1), genetic
information known to correlate with warfarin requirements, demographic
information such as gender, race, and age, medical information such as which
medication they were on at the time of the study and their reason for
receiving warfarin as a treatment, treatment information regarding the
patient’s target international normalised ratio (INR), the therapeutic dose of
warfarin, and whether or not their INR on receiving said dose [51]. The
therapeutic warfarin doses reported in the dataset are the optimal stable
doses for each patient, found by controlled experimentation conducted by
physicians, which allows us to compare the learned policies’ performances
against the true optimal policy given. For our analysis, we remove patients
that had missing entries in any of the covariates listed above. The final
dataset included patient information for 3398 individuals. We used a
discretization of the dosage as the treatment, as is standard with literature
that study the same dataset [51, 37, 36, 35]: $t_{i}=0\ (\leq 3\text{
mg/day}),\ 1\ (>3\text{ and}<5\text{ mg/day}),2\ (\geq 5\text{ mg/day})$. As
for the observed outcome, we consider the absolute different between INR on
Reported Therapeutic Dose of Warfarin from Target INR. Any continuous
covariate was binarized and any categorical covariate was split into one-hot
encoded binary vectors for the prescriptive tree models. Furthermore, we
performed the data imputation steps recommended by [51]. The dataset was
finally split on 5 random seeds with a 70:30 train:test ratio. We used a
5-fold cross validation scheme to fine-tune the $\ell_{0}$ regularization
hyperparameter for the PNN models.

In order to estimate the propensity scores for treatments in this dataset, we
used 10-fold cross validation to evaluate logistic regression models, decision
trees and random forests. Similarly, to estimating the counnterfactual
outcomes for patients in the dataset and compare elastic net rergessors,
random forests and multilayer perceptrons. We conclude that the logistic
regression is ideal to estimate the propensity scores, and that elastic net is
best for estimating the counterfactual outcomes of the patients.

The policies constructed by each prescriptive model are evaluated using
$\pi(s)$, where a function output closer to 0 corresponds to a policy where
most patients were prescribed a dosage that led to their INR being close to
their target INR. The results reported in the following section are the
averages of $\pi(s)$ over the 5 seeds.

###  C.2 Results

We present a summary of our results for the personalized warfarin dosing
experiments in Figures 8 and 9 and a complete table of results can be found in
Appendix D.1.

Figure 8: Results of experiments on the warfarin dataset for out-of-sample
patients prescribed treatments based on policies learned by different models.
The plots display the average outcome. The left plot include the 90%
confidence intervals and the right plot include the 95% confidence intervals.
Figure 9: Paired results of experiments on the warfarin dataset for out-of-
sample patients prescribed treatments based on policies learned by different
models. Each panel displays paired evaluation of the various models. Left
display the point differences of out-of-sample average outcome. Differences
are taken as row - column. Right displays the hypothesis tests on the
differences of out-of-sample average outcome between every pair of models.
Smaller p-values favor models in the rows. Cool colors denote lower blood
pressures (left) or smaller p-values (right). Axes values denote the model (PT
depth/# of trees in CF/PNN width).

Recall that in these experiments, the models work to find the optimal warfarin
dosing policy such that the difference between an individual’s target INR and
their INR post receiving warfarin is minimized. We see that both figures
indicate that every learning model outperforms the policy in the data
evaluated with the doubly robust estimand for the counterfactual outcome.
However, we know that the policy in the data contains the true optimal dosage
for each patient. We note that the best performing models in Figures 8 and 9,
the CF, J-PT and PNN, are those models that incorporate the doubly robust
estimator in their optimization problem. By minimizing over this estimate,
these models heavily rely on the doubly robust estimator in estimating the
counterfactual outcomes in the data. In this case, while a logistic regression
model yielded the best score for estimating the propensity scores in our
cross-validation regime, the accuracy of the model was still only 66.77%.
Similarly, the root mean squared error of the elastic net model that was used
to estimate the counterfactual outcomes with the DM method was 2.0558 on
average across the three treatments (the typical effective therapeutic INR
range for individuals taking warfarin is between 2.0–3.0 and 2.5–3.5 [71]).
This observation is maintained considering that the B-PT and K-PT models
estimate the counterfactual independent from the doubly robust estimation
model. Nevertheless, the left plot in Figure 9 indicates that the paired out-
of-sample mean differences for every pair is within $\pm 0.35$, and so, given
the typical therapeutic INR range, the policies produced by every model are
relatively similar to the optimal policy in the data.

These results shed light on the limitations of the doubly robust estimation
technique, and the models that utilize it their optimization problem. We also
observed that the 3- and 10-width PNN models terminated with an average of
74.6% and 92.78% optimality gap, respectively. This means PNNs would require
addtional computational resources to provide effectiver policies given the
high-dimensional feature vector and the large number of samples in this
setting. Considering that the performance of the PNN models are competitive
given current resources, we expect even stronger performance when additional
resources are available.

##  Appendix D Additional Experimental Results

In this section, we provide additional experimental results that supplement
our analysis.

###  D.1 Complete Experiment Results

####  D.1.1 Simulated Data

The table in this section display the average out-of-sample probability of
correct treatment assignment, average optimality gap, and average running time
in seconds of the models in the simulated data experiments. The averages are
taken over the 5 seeds in the experiments, with the standard deviation given
parenthetically. The results of experiments with 100 and 500 training points
are included.

Train Data Length | Model | Architec -ture | Probability of Correct Assignments
---|---|---|---
0.1 | 0.25 | 0.5 | 0.75 | 0.9
OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s)
100 | B-PT | 1 | 58.34 (18.42) | 0 (0) | 1.51 (0.22) | 62.37 (18.27) | 0 (0) | 1.08 (0.59) | 53.99 (8.9) | 0 (0) | 1.89 (0.38) | 50.09 (0.59) | 0 (0) | 1.27 (0.32) | 44.16 (13.49) | 0 (0) | 1.02 (0.33)
2 | 81.15 (12.77) | 0 (0) | 45.27 (3.42) | 69.44 (17.74) | 0 (0) | 64.91 (15.69) | 75.58 (14.57) | 0 (0) | 90.23 (12.92) | 53.4 (10.66) | 0 (0) | 71.5 (11.76) | 39.86 (7.86) | 0 (0) | 52.97 (12.99)
CF | 100 | 78.59 (16.92) | - | 0.65 (0.01) | 89.96 (1.48) | - | 0.63 (0.01) | 90.46 (1.1) | - | 0.63 (0.01) | 90.48 (0.57) | - | 0.62 (0.01) | 84.03 (14.95) | - | 0.64 (0.02)
J-PT | 1 | 86.99 (4.4) | 0 (0) | 0.11 (0.02) | 89.3 (1.9) | 0 (0) | 0.1 (0.02) | 89.3 (1.9) | 0 (0) | 0.09 (0.01) | 89.97 (1.86) | 0 (0) | 0.09 (0.01) | 89.97 (1.77) | 0 (0) | 0.09 (0.02)
2 | 87.33 (4.44) | 0 (0) | 1.38 (1) | 88.29 (2.02) | 0 (0) | 2.75 (1.42) | 88.76 (1.19) | 0 (0) | 1.6 (1.06) | 87.88 (1.98) | 0 (0) | 2.27 (1.57) | 88.2 (4.39) | 0 (0) | 1.7 (1.58)
K-PT | 1 | 73.33 (15.59) | 0 (0) | 0.39 (0.11) | 73.98 (16.4) | 0 (0) | 0.47 (0.08) | 76.51 (16.91) | 0 (0) | 1.06 (0.41) | 56.03 (8.95) | 0 (0) | 0.84 (0.52) | 42.07 (12.99) | 0 (0) | 0.44 (0.27)
2 | 83.58 (8.02) | 0 (0) | 16.83 (5.82) | 79.02 (12.81) | 0 (0) | 17.85 (5.36) | 83.96 (5.49) | 0 (0) | 29.77 (6.66) | 60.08 (14.19) | 0 (0) | 26.82 (3.57) | 35.37 (13.88) | 0 (0) | 14.52 (3.51)
PNN | 2-3-2 | 88.93 (3.42) | 0 (0) | 22.56 (14.52) | 89.75 (3.39) | 0 (0) | 16.02 (10.87) | 91.27 (0.15) | 0 (0) | 11.61 (2.69) | 90.08 (2.03) | 0 (0) | 14.85 (5.4) | 91.28 (0.15) | 0 (0) | 14.79 (10.02)
2-10-2 | 88.93 (3.43) | 0.02 (0.02) | 3232.61 (821.74) | 89.8 (3.32) | 0.01 (0.01) | 2885.75 (1018.32) | 91.25 (0.15) | 0.01 (0.01) | 3080.71 (1161.3) | 90.14 (1.89) | 0.02 (0.02) | 3268.12 (587.45) | 91.26 (0.15) | 0.02 (0.02) | 2193.95 (1369.05)
500 | B-PT | 1 | 49.8 (0.61) | 0 (0) | 4.77 (1.97) | 57.7 (17.1) | 0 (0) | 7.13 (3.13) | 58.53 (18.56) | 0 (0) | 5.11 (4.71) | 50.2 (0.61) | 0 (0) | 5.54 (3.55) | 50.2 (0.61) | 0 (0) | 4.04 (2.61)
2 | 74.33 (20.4) | 0 (0) | 502.1 (92.75) | 83.2 (12.61) | 0 (0) | 526.6 (146.83) | 85.52 (8.89) | 0 (0) | 570.98 (78.21) | 53.34 (2.53) | 0 (0) | 588.45 (97.26) | 49.36 (1.96) | 0 (0) | 562.64 (83.7)
CF | 100 | 90.95 (0.45) | - | 0.69 (0.02) | 91.21 (0.4) | - | 0.68 (0.01) | 91.17 (0.3) | - | 0.66 (0.02) | 90.9 (0.85) | - | 0.69 (0.02) | 89.96 (1.56) | - | 0.7 (0.01)
J-PT | 1 | 91.41 (0.24) | 0 (0) | 1.13 (0.17) | 91.41 (0.24) | 0 (0) | 0.98 (0.18) | 91.41 (0.24) | 0 (0) | 0.84 (0.03) | 90.66 (1.85) | 0 (0) | 1.03 (0.17) | 90.66 (1.85) | 0 (0) | 1.21 (0.11)
2 | 90.89 (0.44) | 0 (0) | 35.2 (8.09) | 90.95 (0.52) | 0 (0) | 32.71 (14.9) | 90.79 (0.38) | 0 (0) | 27.27 (14.23) | 90.58 (0.94) | 0 (0) | 35.11 (12.88) | 89.75 (0.72) | 0 (0) | 39.09 (3.29)
K-PT | 1 | 80.73 (17.71) | 0 (0) | 2.67 (0.68) | 88.39 (4.91) | 0 (0) | 2.4 (0.85) | 84.79 (9.71) | 0 (0) | 3.4 (0.62) | 52.17 (4.44) | 0 (0) | 3 (0.26) | 50.2 (0.61) | 0 (0) | 1.83 (1.02)
2 | 87.93 (4.11) | 0 (0) | 134.46 (47.93) | 87.71 (4.63) | 0 (0) | 148.33 (35.1) | 83.78 (9.52) | 0 (0) | 146.02 (20.47) | 58.5 (7.11) | 0 (0) | 157.8 (46.88) | 48.92 (2.04) | 0 (0) | 140.41 (34.6)
PNN | 37290 | 90.79 (0.58) | 0.05 (0.04) | 2996.29 (1350.27) | 91.21 (0.1) | 0.02 (0.03) | 2721.05 (1120.22) | 91.21 (0.08) | 0.02 (0.02) | 2764.36 (1144.47) | 91.28 (0.15) | 0.05 (0.04) | 3514.17 (192.23) | 91.21 (0.13) | 0.11 (0.04) | 3600.13 (0.1)
2-10-2 | 91 (0.34) | 0.06 (0.01) | 3600.15 (0.07) | 91.21 (0.09) | 0.06 (0.01) | 3601.7 (3.09) | 91.2 (0.08) | 0.06 (0.01) | 3600.47 (0.27) | 90.96 (0.69) | 0.06 (0.01) | 3600.93 (1.66) | 91.2 (0.12) | 0.08 (0.01) | 3600.6 (0.48)

Table 1: Complete Results Simulated Experiments with Athey Model 1 (with
Adapted Binarization).

Train Data Length | Model | Architec -ture | Probability of Correct Assignments
---|---|---|---
0.1 | 0.25 | 0.5 | 0.75 | 0.9
OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s) | OOSP (%) | Gap (%) | Runtime (s)
100 | B-PT | 1 | 51.56 (3.63) | 0 (0) | 0.93 (0.26) | 51.39 (3.99) | 0 (0) | 1.13 (0.32) | 51.83 (4.38) | 0 (0) | 0.84 (0.11) | 48.13 (4.4) | 0 (0) | 0.85 (0.06) | 43.81 (4.9) | 0 (0) | 0.85 (0.27)
2 | 58.38 (0.85) | 0 (0) | 19.41 (4.42) | 56.68 (7.8) | 0 (0) | 22.86 (3.88) | 53.85 (4.53) | 0 (0) | 24.8 (8.13) | 47.69 (4.82) | 0 (0) | 21.75 (5.7) | 35.27 (4.59) | 0 (0) | 19.8 (3.69)
CF | 100 | 78.59 (16.92) | - | 0.64 (0.01) | 89.96 (1.48) | - | 0.63 (0.01) | 90.46 (1.1) | - | 0.62 (0.02) | 90.48 (0.57) | - | 0.62 (0.01) | 84.03 (14.95) | - | 0.65 (0.01)
J-PT | 1 | 60.21 (0.59) | 0 (0) | 0.09 (0.02) | 59.81 (0.45) | 0 (0) | 0.09 (0.01) | 60.07 (0.65) | 0 (0) | 0.08 (0) | 59.67 (0.36) | 0 (0) | 0.09 (0.01) | 60.07 (0.65) | 0 (0) | 0.08 (0)
2 | 70 (0.88) | 0 (0) | 0.52 (0.19) | 69.79 (0.82) | 0 (0) | 0.43 (0.15) | 70.18 (0.69) | 0 (0) | 0.46 (0.2) | 69.83 (1.05) | 0 (0) | 0.41 (0.15) | 70.26 (0.62) | 0 (0) | 0.41 (0.14)
K-PT | 1 | 54 (5.17) | 0 (0) | 0.21 (0.2) | 56.01 (4.84) | 0 (0) | 0.4 (0.19) | 58.28 (3.39) | 0 (0) | 0.74 (0.43) | 49.46 (7.27) | 0 (0) | 0.36 (0.39) | 41.8 (2.97) | 0 (0) | 0.07 (0.01)
2 | 58.83 (7.29) | 0 (0) | 5.49 (0.85) | 61.01 (4.17) | 0 (0) | 6.02 (1.48) | 59.54 (1.87) | 0 (0) | 7.54 (1.55) | 48.7 (5.54) | 0 (0) | 7.95 (1.41) | 38.27 (2.83) | 0 (0) | 4.63 (1.66)
PNN | 2-3-2 | 88.93 (3.42) | 0 (0) | 22.56 (14.52) | 89.75 (3.39) | 0 (0) | 16.02 (10.87) | 91.27 (0.15) | 0 (0) | 11.61 (2.69) | 90.08 (2.03) | 0 (0) | 14.85 (5.4) | 91.28 (0.15) | 0 (0) | 14.79 (10.02)