# A note on coupled constraint Nash games

Orestes Bueno Universidad del Pacífico. Lima, Perú. Email: {o.buenotangoa,
<EMAIL_ADDRESS>Carlos Calderon Instituto de Matemática y Ciencias
Afines. Email<EMAIL_ADDRESS>John Cotrina11footnotemark: 1

###### Abstract

In this note we are interested in a relevant generalized Nash equilibrium
problem, which was proposed by Rosen in 1965. An existence result is
established in the general setting of quasiconvexity, which is independent
from the one given by Aussel and Dutta in 2008.

Keywords:  Generalized Nash Equilibrium; Generalized Convexity; Variational
Inequalities

MSC (2010): 49J53, 91A10, 91B50

##  1 Introduction

Let $N$ be a finite and non-empty set, representing a set of players. Let us
assume that each player $\nu\in N$ chooses a strategy $x^{\nu}$ in a strategy
set $K_{\nu}$ of $\mathbb{R}^{n_{\nu}}$.Consider the sets

$K=\prod_{\nu\in N}K_{\nu},\quad\mathbb{R}^{n}=\prod_{\nu\in
N}\mathbb{R}^{n_{\nu}},\quad K_{-\nu}=\prod_{\mu\in
N\setminus\\{\nu\\}}K_{\mu},\quad\mathbb{R}^{-n_{\nu}}=\prod_{\mu\in
N\setminus\\{\nu\\}}\mathbb{R}^{n_{\mu}}$

Also, we can write $x=(x^{\nu},x^{-\nu})\in K$ in order to emphasize the
strategy of player $\nu$, $x^{\nu}\in K_{\nu}$, and the strategy of the rival
players $x^{-\nu}\in K_{-\nu}$.

For each player $\nu$, given the rivals’ strategy $x^{-\nu}$, the player $\nu$
chooses a strategy $x^{\nu}$ such that it solves the following optimization
problem

$\textstyle{\operatorname*{NEP}_{\nu}}(x^{-\nu})\>:\>\displaystyle{\min_{K_{\nu}}\theta_{\nu}(\cdot,x^{-\nu})},$
(1)

where $\theta_{\nu}:\mathbb{R}^{n}\to\mathbb{R}$ and
$\theta_{\nu}(x^{\nu},x^{-\nu})$ denotes the _loss_ that player $\nu$ suffers
when his own strategy is $x^{\nu}$ and the rivals’ strategies are $x^{-\nu}$.
Thus, a _Nash equilibrium_ [26], associated to the loss functions
$(\theta_{\nu})_{\nu}$ and the constraint sets $(K_{\nu})_{\nu}$, is a vector
$\hat{x}\in K$ such that $\hat{x}^{\nu}$ solves
$\operatorname*{NEP}_{\nu}(\hat{x}^{-\nu})$, for any $\nu$.

It is well known [18, Theorem 1] in the literature that a Nash equilibrium
exists when, for each $\nu\in N$,

  * •

the set $K_{\nu}$ is compact convex and non-empty,

  * •

the function $\theta_{\nu}$ is continuous and is quasiconvex on its player’s
variable.

Since then, many authors extended this result in different ways, see for
instance [18, 32, 30, 28, 9], where continuity assumption is relaxed, or [29,
22], where the quasiconvexity is weakened.

On the other hand, in a Generalized Nash Equilibrium Problem (GNEP), each
player $\nu$ posseses a _constraint map_
$\mathfrak{X}_{\nu}:K_{-\nu}\rightrightarrows K_{\nu}$, so their strategy must
belong to the set $\mathfrak{X}_{\nu}(x^{-\nu})\subset K_{\nu}$, where the
sets $\mathfrak{X}_{\nu}(x^{-\nu})$ clearly depend on the rival players’
strategies. The aim of player $\nu$, given the rival players’ strategies
$x^{-\nu}$, is to choose a strategy $x^{\nu}$ that solves the minimization
problem

$\textstyle{\operatorname*{GNEP}_{\nu}}(x^{-\nu})\>:\>\displaystyle{\min_{\mathfrak{X}_{\nu}(x^{-\nu})}\theta_{\nu}(\cdot,x^{-\nu})}.$
(2)

A _Generalized Nash Equilibrium_ is a vector $\hat{x}\in K$ such that
$\hat{x}^{\nu}$ solves $\operatorname*{GNEP}_{\nu}(\hat{x}^{-\nu})$, for any
$\nu$. Arrow and Debreu [3, Lemma 2.5] proved the existence of Generalized
Nash Equilibria when for all $\nu\in N$,

  * •

the set $K_{\nu}$ is compact convex and non-empty,

  * •

the function $\theta_{\nu}$ is continuous and quasiconvex on its player’s
variable,

  * •

the map $\mathfrak{X}_{\nu}$ is lower semicontinuous with closed graph and
convex and non-empty values.

Clearly this existence result generalizes the first one, by considering
$\mathfrak{X}_{\nu}$ as the constant map
$\mathfrak{X}_{\nu}(x^{-\nu})=K_{\nu}$. Since then, many authors dealt with
this kind of problems, see for instance [21, 33, 32, 17, 8, 6, 16, 7, 15, 14,
13].

An important instance of GNEPs was presented by Rosen in [31], which is known
as the _jointly convex GNEP_. A jointly convex GNEP is a GNEP in which the
graph of each map constraint $\mathfrak{X}_{\nu}$ is a certain convex, compact
and non-empty subset of $\mathbb{R}^{n}$, denoted by $\mathfrak{X}$. In other
words

$\mathfrak{X}_{\nu}(x^{-\nu})=\\{x^{\nu}\in\mathbb{R}^{n_{\nu}}\>:\>(x^{\nu},x^{-\nu})\in\mathfrak{X}\\},$

When this is the case, $K_{\nu}$ is the projection of $\mathfrak{X}$ onto
$\mathbb{R}^{n_{\nu}}$. This kind of GNEP have recently begun to gain more and
more attention as it models problems from economics such as electricity
markets, environmental games and bilateral exchange of bads, see [11, 23, 12,
24, 34].

On the other hand, jointly convex GNEPs can be reformulated as variational
inequalities [19] when assuming pseudoconvexity and differentiability. An
extension of this reformulation was presented by Aussel and Dutta in [5],
using normal cones, under quasiconvexity and continuity assumptions. However,
their main result concerning the existence of equilibria was established under
semi strict quasiconvexity.

In this note, our main result consists to assume only quasiconvexity in order
to show the existence of equilibria for the jointly convex GNEP. In Section 2
we recall some definitions of generalized convexity, continuity for functions
and continuity for set-valued maps. Additionally some preliminary results are
established concerning closedness and lower semicontinuity for special maps.
Finally, in Section 3 we reformulate jointly convex GNEPs as variational
inequality problems and show our main result.

##  2 Definitions, Notations and Preliminary Results

Let $U,V$ be non-empty sets. A _multivalued_ or _set-valued_ map
$T:U\rightrightarrows V$ is an application $T:U\to\mathcal{P}(V)$, that is,
for $u\in U$, $T(u)\subset V$. The _graph_ of $T$ is defined as

$\displaystyle\operatorname{gra}(T)=\big{\\{}(u,v)\in U\times V\>:\>v\in
T(u)\big{\\}}.$

Let $X,Y$ be two real Banach spaces. Recall that a set-valued map
$T:X\rightrightarrows Y$ is

  * •

_lower semicontinuous_ (lsc) at $x_{0}\in X$ when, for any open set $V$ such
$T(x_{0})\cap V\neq\emptyset$, there exists a neighborhood $U$ of $x_{0}$ such
that $T(x)\cap V\neq\emptyset$, for every $x\in U$;

  * •

_upper semicontinuous_ (usc) at $x_{0}\in X$ when, for any neighborhood $V$ of
$T(x_{0})$, there exists a neighborhood $U$ of $x_{0}$ such that $T(x)\subset
V$, for all $x\in U$;

  * •

_lower_ (respectively, _upper_) _semicontinuous_ when it is lower (resp.
upper) semicontinous at every $x_{0}\in X$;

  * •

_closed_ when its graph is a closed subset of $X\times Y$.

Moreover, $T:X\rightrightarrows Y$ is lower semicontinuous at $x_{0}$ if, and
only if, $T$ is _sequentially lower semicontinuous_ at $x_{0}$, that is, for
any sequence $(x_{n})_{n\in\mathbb{N}}$ of $X$ converging to $x_{0}$ and any
$y_{0}\in T(x_{0})$, there exists a sequence $(y_{n})_{n\in\mathbb{N}}$ of $Y$
converging to $y_{0}$ such that $y_{n}\in T(x_{n})$, for any $n$. The
interested reader can find a proof of this fact in [20, Proposition 2.5.6].

Let $X$ be a Banach space, $X^{*}$ its topological dual and
$\langle\cdot,\cdot\rangle$ the duality pairing. Given $A$ a subset of $X$,
the _closure_ , the _interior_ and the _convex hull_ of $A$ are denoted by
$\bar{A}$, $\operatorname*{int}(A)$ and $\operatorname*{conv}(A)$,
respectively. In addition, $-A$ denotes the set $-A=\\{x\in X\>:\>-x\in A\\}$.
The _polar cone_ of $A$ is the set

$A^{-}:=\\{x^{*}\in X^{*}\>:\>\langle x^{*},x\rangle\leq 0,\,\forall\,x\in
A\\},$

which is a convex cone on $X^{*}$. The _normal cone_ $N_{A}$ of $A$ at $x\in
X$ is defined as

$N_{A}(x):=\\{x^{*}\in X^{*}\>:\>\langle x^{*},y-x\rangle\leq
0,\,\forall\,y\in A\\}.$

We will consider $N_{A}(x)=X^{*}$, whenever $A=\emptyset$. It is common in
this definition to impose $A$ convex and $x\in A$, however, we won’t consider
such conditions in this work. See [10] for another instance where such
conditions are not considered. Note that, for all $x\in X$,

$N_{A}(x)=(A-x)^{-}=N_{\bar{A}}(x)=N_{\operatorname*{conv}(A)}(x).$

The following lemma will be used in the proof of Proposition 2.4.

######  Lemma 2.1.

Let $A$ be a subset of a real Banach space $X$. If
$\operatorname*{int}(A)\neq\emptyset$, then $A^{-}$ is a pointed cone, i.e.
$A^{-}\cap-A^{-}=\\{0\\}$. Moreover, if $x\in\operatorname*{int}(A)$ and
$x^{*}\in A^{-}$, $x^{*}\neq 0$, then $\langle x^{*},x\rangle<0$.

###### Proof.

Note that $A^{-}\cap-A^{-}=A^{\perp}=\\{x^{*}\in X^{*}\>:\>\langle
x,x^{*}\rangle=0,\,\forall\,x\in A\\}$. Since
$\operatorname*{int}(A)\neq\emptyset$, $\operatorname{span}(A)=X$, so

$A^{-}\cap-A^{-}=A^{\perp}=\operatorname{span}(A)^{\perp}=X^{\perp}=\\{0\\}.$

To conclude the proof, let $x\in\operatorname*{int}(A)$, $x^{*}\in A^{-}$,
$x^{*}\neq 0$ and assume that $\langle x,x^{*}\rangle\geq 0$, hence $\langle
x,x^{*}\rangle=0$. Since $x\in\operatorname*{int}(A)$, for some
$\varepsilon>0$ and all $v\in B(0,1)$, $x+\varepsilon v\in A$, which implies
$\langle x+\varepsilon v,x^{*}\rangle\leq 0$. Therefore $\langle
v,x^{*}\rangle\leq 0$, for all $v\in B(0,1)$, which in turn implies that
$\langle v,x^{*}\rangle=0$, for all $v\in B(0,1)$. This is a contradiction,
since $x^{*}\neq 0$. ∎

Let us now recall some classical definitions of generalized convexity. Given a
Banach space $X$, a real-valued function $h:X\to\mathbb{R}$ is

  * •

_quasiconvex_ if, for any $x,y\in X$ and $t\in[0,1]$,

$h(tx+(1-t)y)\leq\max\\{h(x),h(y)\\};$

  * •

_semistrictly quasiconvex_ if it is quasiconvex and, for any $x,y\in X$ and
$t\in]0,1[$,

$h(x)\neq h(y)\quad\Rightarrow\quad h(tx+(1-t)y)<\max\\{h(x),h(y)\\}.$

Our definition of semistrictly quasiconvex function was taken from [2, Chapter
5], and it is equivalent to the definition of pseudoconvex function given in
[27, Definition 3.1]. Some authors may not include the quasiconvex condition
in this definition, see for instance [2, Chapter 4].

It is well known that a function $h$ is quasiconvex if, and only if, the
_sublevel sets_ $S_{h,\lambda}=\\{y\in X\>:\>h(y)\leq\lambda\\}$ are convex,
for all $\lambda\in\mathbb{R}$, and also if and only if the _strict sublevel
sets_ $S^{<}_{h,\lambda}=\\{y\in X\>:\>h(y)<\lambda\\}$ are convex, for all
$\lambda\in\mathbb{R}$.

Let $X$ be a Banach space. Given $x\in X$, $\mathscr{V}_{x}$ will denote an
open neighbourhood of $x$ in $X$. A function $h:X\to\mathbb{R}$ is called:

  * •

_lower semicontinuous_ , if for all $x\in X$ and $\lambda\in\mathbb{R}$ with
$\lambda<h(x)$, there exists $\mathscr{V}_{x}$ such that
$\lambda<h(x^{\prime})$, for all $x^{\prime}\in\mathscr{V}_{x}$;

  * •

_upper semicontinuous_ , if for all $x\in X$ and $\lambda\in\mathbb{R}$ with
$\lambda>h(x)$, there exists $\mathscr{V}_{x}$ such that
$\lambda>h(x^{\prime})$, for all $x^{\prime}\in\mathscr{V}_{x}$;

  * •

_lower pseudocontinuous_ [25], if for all $x,y\in X$ with $h(y)<h(x)$, there
exists $\mathscr{V}_{x}$ such that $h(y)<h(x^{\prime})$, for all
$x^{\prime}\in V_{x}$;

  * •

_upper pseudocontinuous_ , if for all $x,y\in X$ with $h(y)>h(x)$, there
exists $\mathscr{V}_{x}$ such that $h(y)>h(x^{\prime})$, for all
$x^{\prime}\in V_{x}$.

It is not difficult to prove that if $h$ is lower semicontinuous then it is
lower pseudocontinuous. On the other hand, it is well known that $h$ is lower
semicontinuous if, and only if, $S_{h,\lambda}$ is closed, for all
$\lambda\in\mathbb{R}$. In the same way, considering for $x\in X$ the sets
$S_{h}(x):=S_{h,h(x)}$, $S^{<}_{h}(x):=S^{<}_{h,h(x)}$, $h$ is lower
pseudocontinuous if, and only if, $S_{h}(x)$ is closed, for all $x\in X$.
Also, $h$ is upper pseudocontinuous if, and only if, $S^{<}_{h}(x)$ is open,
for all $x\in X$.

Let $X$, $Y$ be two real Banach spaces and $f:X\times Y\to\mathbb{R}$ be a
function. We define the following set-valued maps $L^{<}_{f},L_{f}:X\times
Y\rightrightarrows X$ as

$L^{<}_{f}(x,y):=S^{<}_{f(\cdot,y)}(x)\quad\text{and}\quad
L_{f}(x,y):=S_{f(\cdot,y)}(x).$ (3)

and the set-valued maps $\mathscr{N}_{f}^{<},\mathscr{N}_{f}:X\times Y\to
X^{*}$

$\mathscr{N}_{f}^{<}(x,y):=N_{L^{<}_{f}(x,y)}(x)\quad\text{and}\quad\mathscr{N}_{f}(x,y):=N_{L_{f}(x,y)}(x).$
(4)

It is important to note that $L^{<}_{f}(x,y)$ is not merely the projection
onto $X$ of $S^{<}_{f}(x,y)$. Consider for instance the function
$f:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$, $f(x,y)=x^{2}+y^{2}$, and take
$(0,1)$. It is not difficult to see that $L_{f}^{<}(0,1)=\emptyset$, but the
projection of $S^{<}_{f}(0,1)$ onto $\mathbb{R}$ is the interval $]-1,1[$.

In general, we have the following relation between $L_{f}^{<}(x,y)$ and
$S_{f}^{<}(x,y)$:

$L^{<}_{f}(x,y)\times\\{y\\}=S^{<}_{f(x,y)}\cap(X\times\\{y\\}).$

The following result establishes the lower semicontinuity of $L^{<}_{f}$,
which is inspired by [1].

######  Lemma 2.2.

Let $X,Y$ be two real Banach spaces and $f:X\times Y\to\mathbb{R}$ be a
function. If $f$ is lower pseudocontinuous on its first argument and
continuous on its second one, then the map $L^{<}_{f}$ is lower
semicontinuous.

###### Proof.

Assume, on the contrary, that $L^{<}_{f}$ is not lower semicontinuous at some
$(x_{0},y_{0})\in X\times Y$. Then, there exists an open set $V\subset X$ with
$L^{<}_{f}(x_{0},y_{0})\cap V\neq\emptyset$ and, for every $n,m\in\mathbb{N}$,
there exist $x_{n}\in B(x_{0},1/n)$, $y_{m}\in B(y_{0},1/m)$, such that

$L^{<}_{f}(x_{n},y_{m})\cap V=\emptyset.$

This means that $f(z,y_{m})\geq f(x_{n},y_{m})$, for all $z\in V$ and
$n,m\in\mathbb{N}$. Since $f$ is continuous on its second argument, taking
limit when $m\to\infty$, we obtain $f(z,y_{0})\geq f(x_{n},y_{0})$, that is,
$x_{n}\in L_{f}(z,y_{0})$, for all $n\in\mathbb{N}$. We now take limit when
$n\to\infty$ to conclude $x_{0}\in L_{f}(z,y_{0})$, since $L_{f}(z,y_{0})$ is
closed, as $f$ is lower pseudocontinuous on its first argument. Therefore
$f(x_{0},y_{0})\leq f(z,y_{0})$, for all $z\in V$. The lemma follows by
observing that, $f(z_{0},y_{0})<f(x_{0},y_{0})$, for some $z_{0}\in V$. ∎

######  Lemma 2.3.

Let $X,Y$ be two real Banach spaces and $T:X\times Y\rightrightarrows X$ be a
set-valued map. If $T$ is lower semicontinuous, then the set-valued map
$\mathcal{N}_{T}:X\times Y\rightrightarrows X^{*}$ defined as

$\mathcal{N}_{T}(x,y):=N_{T(x,y)}(x)$

is closed.

###### Proof.

Let $(x_{n},y_{n},x_{n}^{*})_{n\in\mathbb{N}}$ be a sequence in the graph of
$\mathcal{N}_{T}$ converging to $(x_{0},y_{0},x_{0}^{*})$. We aim to show that
$x_{0}^{*}\in\mathcal{N}_{T(x_{0},y_{0})}(y_{0})$. If
$T(x_{0},y_{0})=\emptyset$, there is nothing to prove. Now, assume that
$T(x_{0},y_{0})\neq\emptyset$ and take any $x\in T(x_{0},y_{0})$. By lower
semicontinuity of $T$ there exists a sequence $(z_{n})_{n\in\mathbb{N}}$
converging to $x$ such that $z_{n}\in T(x_{n},y_{n})$, for all
$n\in\mathbb{N}$. Therefore

$\langle x_{n}^{*},z_{n}-x_{n}\rangle\leq 0.$

The lemma follows by letting $n$ tend to $\infty$, to obtain $\langle
x_{0}^{*},x-x_{0}\rangle\leq 0$. ∎

######  Proposition 2.4.

Let $X,Y$ be two real Banach spaces and $f:X\times Y\to\mathbb{R}$ be a
function.

  1. 1.

If $f$ is lower pseudocontinuous on its first argument and continuous on its
second one, then the map $\mathscr{N}^{<}_{f}$ is closed.

  2. 2.

If $f$ is upper pseudocontinuous on its first argument then $L^{<}_{f}(x,y)$
is open, for all $(x,y)\in X\times Y$. In particular, if
$w^{*}\in\mathscr{N}^{<}_{f}(x,y)$, $w^{*}\neq 0$, then $\langle
w^{*},z-x\rangle<0$, for all $z\in L^{<}_{f}(x,y)$.

  3. 3.

If $f$ is upper pseudocontinuous and quasiconvex, both on its first argument,
then $\mathscr{N}^{<}_{f}(x,y)\neq\\{0\\}$, for all $(x,y)\in X\times Y$.

###### Proof.

  1. 1.

Closedness of $\mathscr{N}^{<}_{f}$ follows from Lemmas 2.2 and 2.3, and the
fact that $\mathscr{N}^{<}_{f}(x,y)=\mathcal{N}_{L^{<}_{f}}(x,y)$, for any
$(x,y)\in X\times Y$.

  2. 2.

The set $L^{<}_{f}(x,y)$ is open, since $L^{<}_{f}(x,y)=S^{<}_{f(\cdot,y)}(x)$
and $f(\cdot,y)$ is upper pseudocontinuous. The second claim follows from the
fact that $\mathscr{N}^{<}_{f}(x,y)=(L^{<}_{f}(x,y)-x)^{-}$ and Lemma 2.1.

  3. 3.

Take $(x,y)\in X\times Y$. If $x\in\operatorname*{arg\,min}_{X}f(\cdot,y)$
then $L^{<}_{f}(x,y)=\emptyset$ and the claim trivially follows. Now assume
that $x\notin\operatorname*{arg\,min}_{X}f(\cdot,y)$. The quasiconvexity of
$f(\cdot,y)$ implies $S^{<}_{f(\cdot,y)}(x)=L^{<}_{f}(x,y)$ is non-empty and
convex. Moreover $x\notin L^{<}_{f}(x,y)$ so, using the Hahn-Banach separation
theorem, there exists $x^{*}\in X^{*}$, $x^{*}\neq 0$ such that

$\langle x^{*},z\rangle\leq\langle x^{*},x\rangle,\quad\forall z\in
L^{<}_{f}(x,y).$

This implies that $x^{*}\in\mathscr{N}^{<}_{f}(x,y)$.∎

######  Remark 2.5.

The pseudocontinuity condition on $f$ in item 3 can be dropped if $X$ is
finite-dimensional.

The following example shows that continuity and convexity only in the first
argument of $f$ is not enough to guarantee the closedness of
$\mathscr{N}^{<}_{f}$, neither the lower semicontinuity of $L^{<}_{f}$.

######  Example 2.6.

Define $\Theta:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ as

$\Theta(x,y)=\begin{cases}x,&\text{when }y\neq 1,\\\ -x,&\text{when
}y=1.\end{cases}$

It is not difficult to see that $\Theta$ is continuous and convex with respect
to its first argument. Moreover, for each $(x,y)\in\mathbb{R}^{2}$,

$L^{<}_{\Theta}(x,y)=\begin{cases}]-\infty,x[,&\text{when }y\neq 1,\\\
]x,+\infty[,&\text{when }y=1.\end{cases}$

Clearly the map $L^{<}_{\Theta}$ is not lower semicontinuous at $(0,1)$. In
addition, the strict normal operator associated to $\Theta$ is given by

$\mathscr{N}^{<}_{\Theta}(x,y)=\begin{cases}[0,+\infty[,&\text{when }y\neq
1\\\ ]-\infty,0],&\text{when }y=1.\end{cases}$

which is not closed.

##  3 GNEPs and Variational Inequality Problems

Throughout this section, we will consider a jointly convex GNEP, given by a
convex set $\mathfrak{X}\subset\mathbb{R}^{n}$. We will also consider the loss
functions $\theta_{\nu}$ as bifunctions from
$\mathbb{R}^{n_{\nu}}\times\mathbb{R}^{-n_{\nu}}$ to $\mathbb{R}$.

For each $\nu\in N$ and $x\in\mathbb{R}^{n}$, consider the set [5, 27]

$D_{\nu}(x)=\operatorname*{conv}\displaystyle\left(\mathscr{N}^{<}_{\theta_{\nu}}(x^{\nu},x^{-\nu})\cap
S_{\nu}[0,1]\right),$

where $S_{\nu}[0,1]$ is the unit sphere in $\mathbb{R}^{n_{\nu}}$. This allows
us to define the set-valued map
$T:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ as

$T(x):=\prod_{\nu\in N}D_{\nu}(x).$ (5)

Clearly, for all $x\in\mathbb{R}^{n}$, $T(x)$ is compact and convex, as each
$D_{\nu}(x)$ is compact and convex aswell. However, the operator $T$ may not
be a closed operator, as shown by the following example.

######  Example 3.1.

Consider a two player NEP where each player $\nu\in\\{1,2\\}$ has a loss
function defined as in Example 2.6, namely

$\theta_{\nu}(x^{\nu},x^{-\nu})=\Theta(x^{\nu},x^{-\nu})=\begin{cases}x^{\nu},&\text{when
}x^{-\nu}\neq 1,\\\ -x^{\nu},&\text{when }x^{-\nu}=1.\end{cases}$

Therefore

$\mathscr{N}^{<}_{\theta_{\nu}}(x)=\begin{cases}[0,+\infty[,&\text{when
}x^{-\nu}\neq 1,\\\ ]-\infty,0],&\text{when }x^{-\nu}=1,\end{cases}$

and

$D_{\nu}(x)=\begin{cases}\\{1\\},&\text{when }x^{-\nu}\neq 1,\\\
\\{-1\\},&\text{when }x^{-\nu}=1,\end{cases}$

which is not closed, so neither is $T$.

In view of the previous example, we need additional conditions to guarantee
the closedness of $T$.

######  Proposition 3.2.

Let $T:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ be defined as in (5).

  1. 1.

If each loss function $\theta_{\nu}$ is lower pseudocontinuous on its own
player’s variable and continuous with respect to its rivals’ variables, then
$T$ is closed.

  2. 2.

If each loss function $\theta_{\nu}$ is quasiconvex with respect to its own
player’s variable then $T$ is non-empty valued.

###### Proof.

  1. 1.

From Proposition 2.4, item 1, we have that each map
$\mathscr{N}^{<}_{\theta_{\nu}}$ is closed with convex values. This implies
that the map $D_{\nu}$ is closed with convex values. The result follows from
the fact that the Cartesian product of closed maps is a closed map.

  2. 2.

In view of Remark 2.5, we can use Proposition 2.4, item 3, to conclude that
the set $\mathscr{N}^{<}_{\theta_{\nu}}(x^{\nu},x^{-\nu})\setminus\\{0\\}$ is
non-empty, for all $(x^{\nu},x^{-\nu})\in\mathbb{R}^{n}$. Thus, the map
$D_{\nu}$ is non-empty valued and so is $T$.

∎

Let $S(T,\mathfrak{X})$ be the solution set of the Variational Inequality
Problem associated to the map $T$ (given as in (5)) and the set
$\mathfrak{X}$, namely

$S(T,\mathfrak{X})=\\{x\in\mathfrak{X}\>:\>\exists\,w_{0}\in T(x),\,\langle
w_{0},y-x\rangle\geq 0,\,\forall\,y\in\mathfrak{X}\\}$

The following result establishes a link between the jointly convex GNEP and
its associated variational inequality problem.

######  Proposition 3.3.

Assume that $\mathfrak{X}$ is non-empty. If every loss function $\theta_{\nu}$
is upper pseudocontinuous on its own player’s variable, then every point in
$S(T,\mathfrak{X})$ is a generalized Nash equilibrium.

###### Proof.

Take $x\in S(T,\mathfrak{X})$, so there exists $w_{0}\in T(x)$ such that

$\displaystyle\langle w_{0},y-x\rangle\geq 0,\quad\forall\,y\in\mathfrak{X}.$
(6)

Assume that $x$ is not a solution of the jointly convex GNEP. Then there exist
a player $\nu$, and $z^{\nu}\in\mathfrak{X}_{\nu}(x^{-\nu})$, such that

$\theta_{\nu}(z^{\nu},x^{-\nu})<\theta_{\nu}(x^{\nu},x^{-\nu}),$

that is $z^{\nu}\in L^{<}_{\theta_{\nu}}(x^{\nu},x^{-\nu})$ and
$y=(z^{\nu},x^{-\nu})\in\mathfrak{X}$. In view of (6),

$\langle w_{0}^{\nu},z^{\nu}-x^{\nu}\rangle=\langle w_{0},y-x\rangle\geq 0.$

This in turn implies that $w_{0}^{\nu}=0$, by Proposition 2.4, item 2, and the
fact that $w_{0}^{\nu}\in\mathscr{N}^{<}_{\theta_{\nu}}(x^{\nu},x^{-\nu})$.

On the other hand, since $w_{0}^{\nu}\in D_{\nu}(x)$, there exist
$w^{\nu}_{1},\ldots,w^{\nu}_{p}\in\mathscr{N}^{<}_{\theta_{\nu}}(x^{\nu},x^{-\nu})\cap
S_{\nu}[0,1]$, $t_{1},\ldots,t_{p}\geq 0$,
$\displaystyle\sum_{i=1}^{p}t_{i}=1$, such that
$0=w_{0}^{\nu}=\displaystyle\sum_{i=1}^{p}t_{i}w_{i}^{\nu}$. Take $i_{0}$ such
that $t_{i_{0}}>0$, so we have

$0=\sum_{i\neq i_{0}}t_{i}w_{i}^{\nu}+t_{i_{0}}w^{\nu}_{i_{0}}$

and this implies

$-w_{i_{0}}^{\nu}=\sum_{i\neq i_{0}}\frac{t_{i}}{t_{i_{0}}}w_{i}^{\nu}.$

As $\mathscr{N}^{<}_{\theta_{\nu}}(x)$ is a convex cone,
$-w_{i_{0}}^{\nu}\in\mathscr{N}^{<}_{\theta_{\nu}}(x)$, hence

$w_{i_{0}}^{\nu}\in\mathscr{N}^{<}_{\theta_{\nu}}(x)\cap-\mathscr{N}^{<}_{\theta_{\nu}}(x).$

However, since $L^{<}_{\theta_{\nu}}(x)$ is open, Lemma 2.1 implies that
$w_{i_{0}}=0$, a contradiction. The proposition follows. ∎

Proposition 3.3 is strongly related to [5, Theorem 3.1]. However, the authors
considered continuity instead of upper pseudocontinuity.

######  Proposition 3.4.

Assume that $\mathfrak{X}$ is convex, compact and non-empty and let $T$ be
defined as in (5). If for all $\nu\in N$ the following hold

  1. 1.

the function $\theta_{\nu}$ is lower pseudocontinuous on its own player’s
variable and continuous with respect to its rivals’ variables,

  2. 2.

the function $\theta_{\nu}$ is quasiconvex with respect to its own player’s
variable;

then $S(T,\mathfrak{X})$ is non-empty.

###### Proof.

Thanks to Proposition 3.2, the map $T$ is closed with convex, compact and
nonempty values. Thus, $T$ is upper semicontinuous. Finally, the result
follows from [4, Theorem 9.9]. ∎

Finally, we are ready for our main result, which establishes the existence of
solution for jointly convex GNEPs.

######  Theorem 3.5.

Assume that $\mathfrak{X}$ is convex, compact and non-empty. If for all
$\nu\in N$ the following hold

  1. 1.

the function $\theta_{\nu}$ is pseudocontinuous on its own player’s variable
and continuous on its rivals’ variables,

  2. 2.

the function $\theta_{\nu}$ is quasiconvex on its own player’s variable,

then there exists a generalized Nash equilibrium.

###### Proof.

It is a consequence of Propositions 3.3 and 3.4. ∎

######  Remark 3.6.

Theorem 3.5 is not a consequence of [3, Theorem 2]. In fact, Theorem 3.5
improves [5, Theorem 4.2], [19, Theorem 2.1], [31, Theorem 1] and [18, Theorem
1]

The following example shows that we cannot drop the continuity of each loss
function with respect to its rivals’ variables.

######  Example 3.7.

Given the functions $\theta_{1},\theta_{2}:\mathbb{R}^{2}\to\mathbb{R}$
defined as

$\theta_{1}(x^{1},x^{2}):=\left\\{\begin{matrix}\left(x^{1}-\frac{\sqrt{2}}{2}\right)^{2},&x^{2}\in\mathbb{Q}\\\
(2x^{1}-x^{2})^{2},&x^{2}\notin\mathbb{Q}\end{matrix}\right.\mbox{ and
}\theta_{2}(x^{1},x^{2}):=\left\\{\begin{matrix}\left(x^{2}-x^{1}\right)^{2},&x^{1}\notin\mathbb{Q}\\\
\left(2x^{2}-x^{1}\right)^{2},&x^{1}\in\mathbb{Q}\end{matrix}\right.$

Clearly, each loss function is convex and continuous with respect to its own
variable. For $X=[0,1]^{2}$, the GNEP reduces to the classic Nash equilibrium
problem. Furthermore, it is no difficult to see that there is not solution to
this GNEP.

We must note that the previous example directly contradicts [27, Corollary
4.3].

## References

  * [1] S. Al-Homidan, N. Hadjisavvas, and L. Shaalan.  Transformation of quasiconvex functions to eliminate local minima.  J. Optim. Theory and Appl., 177(1):93–105, 2018.
  * [2] S. A. R. Al-Mezel, F. R. M. Al-Solamy, and Q. H. Ansari, editors.  Fixed point theory, variational analysis, and optimization.  CRC Press, Boca Raton, FL, 2014.
  * [3] K. J. Arrow and G. Debreu.  Existence of an equilibrium for a competitive economy.  Econometrica, 22(3):265–290, 1954.
  * [4] J. P. Aubin.  Optima and equilibria : an introduction to nonlinear analysis / Jean-Pierre Aubin ; translated from the French by Stephen Wilson. Graduate texts in mathematics, 140. Springer, Berlin ;, 2nd ed. edition, 1998.
  * [5] D. Aussel and J. Dutta.  Generalized Nash Equilibrium Problem, Variational Inequality and Quasiconvexity.  Oper. Res. Lett., 36(4):461–464, 2008.
  * [6] D. Aussel, R. Gupta, and A. Mehra.  Evolutionary Variational Inequality Formulation of the Generalized Nash Equilibrium Problem.  J. Optim. Theory Appl., 169(1):74–90, 2016.
  * [7] D. Aussel and A. Sultana.  Quasi-variational inequality problems with non-compact valued constraint maps.  J. Math Anal. Appl., 456(2):1482–1494, 2017.
  * [8] D. Aussel, A. Sultana, and V. Vetrivel.  On the existence of projected solutions of quasi-variational inequalities and generalized Nash equilibrium problems.  J. Optim. Theory Appl., 2016.
  * [9] J. Banks and J. Duggan.  Existence of Nash Equilibria on Convex Sets.  Technical report, 2004.
  * [10] O. Bueno and J. Cotrina.  On maximality of quasimonotone operators.  Set-Valued and Variational Analysis, 27(1):87–101, 2019.
  * [11] J. Contreras, M. Klusch, and J. B. Krawczyk.  Numerical solutions to nash-cournot equilibria in coupled constraint electricity markets.  IEEE Transactions on Power Systems, 19(1):195–206, 2004.
  * [12] J. Contreras, J. B. Krawczyk, and J. Zuccollo.  Economics of collective monitoring: a study of environmentally constrained electricity generators.  Computational Management Science, 13(3):349–369, 2016.
  * [13] J. Cotrina, A. Hantoute, and A. Svensson.  Existence of quasi-equilibria on unbounded constraint sets.  Optimization, pages 1–18, 2020.
  * [14] J. Cotrina and A. Svensson.  The finite intersection property for equilibrium problems.  J Global Optim, 2020.
  * [15] J. Cotrina and J. Zúñiga.  Time-dependent generalized Nash equilibrium problem.  J. Optim. Theory Appl., 179:1054–1064, 2018.
  * [16] J. Cotrina and J. Zúñiga.  Quasi-equilibrium problems with non-self constraint map.  J. Glob Optim, 75:177–197, 2019.
  * [17] P. Cubiotti.  Existence of nash equilibria for generalized games without upper semicontinuity.  International Journal of Game Theory, 26(2):267–273, 1997.
  * [18] P. Dasgupta and E. Maskin.  The existence of equilibrium in discontinuous economic games, part i (theory).  Review of Economic Studies, 53(1):1–26, 1986.  Reprinted in K. Binmore and P. Dasgupta (eds.), Economic Organizations as Games, Oxford: Basil Blackwell, 1986, pp. 48-82.
  * [19] F. Facchinei, A. Fischer, and V. Piccialli.  On generalized nash games and variational inequalities.  Oper. Res. Let, 35(2):159 – 164, 2007.
  * [20] A. Göpfert, H. Riahi, C. Tammer, and C. Zălinescu.  Variational Methods in Partially Ordered Spaces.  CMS Books in Mathematics. Springer, New York, NY, 2003.
  * [21] P. T. Harker.  Generalized Nash games and quasi-variational inequalities.  European Journal of Operational Research, 54(1):81 – 94, 1991.
  * [22] P. Q. Khanh and N. H. Quan.  Versions of the Weierstrass theorem for bifunctions and solution existence in optimization.  SIAM J. Optim, 29(2):1502–1523, 2019.
  * [23] J. B. Krawczyk.  Coupled constraint nash equilibria in environmental games.  Resource and Energy Economics, 27(2):157 – 181, 2005.
  * [24] J. B. Krawczyk and M. Tidball.  Economic problems with constraints: How efficiency relates to equilibrium.  International Game Theory Review, 18(04):1650011, 2016.
  * [25] J. Morgan and V. Scalzo.  Pseudocontinuous functions and existence of nash equilibria.  Journal of Mathematical Economics, 43(2):174–183, 2007.
  * [26] J. Nash.  Non-cooperative games.  Annals of Mathematics, 54(2):286–295, 1951.
  * [27] M. Nasri and W. Sosa.  Equilibrium problems and generalized Nash games.  Optimization, 60:1161–1170, 2011.
  * [28] R. Nessah and G. Tian.  Existence of Solution of Minimax Inequalities, Equilibria in Games and Fixed Points Without Convexity and Compactness Assumptions.  J. Optim. Theory Appl., 157:756–95, 2013.
  * [29] K. Nishimura and J. Friedman.  Existence of nash equilibrium in n person games without quasi-concavity.  International Economic Review, 22(3):637–648, 1981.
  * [30] P. J. Reny.  On the Existence of Pure and Mixed Strategy Nash Equilibria in Discontinuous Games.  Econometrica, 67(5):1029–1056, September 1999.
  * [31] J. B. Rosen.  Existence and uniqueness of equilibrium points for concave n-person games.  Econometrica, 33(3):520–534, 1965.
  * [32] G. Tian and J. Zhou.  Transfer continuities, generalizations of the Weierstrass and maximum theorems: a full characterization.  J. Math. Econ., 24:281–303, 1995.
  * [33] G. Q. Tian.  Generalization of the KKM Theorem and the Ky Fan Minimax Inequality, with Applications to Maximal elements, Price equilibrium, and Complementary.  J. Math. Anal. Appl., 170:457–471, 1992.
  * [34] B. Vardar and G. Zaccour.  Strategic bilateral exchange of a bad.  Operations Research Letters, 47(4):235 – 240, 2019.


:

# Instability in clinical risk stratification models
using deep learning

Daniel Lopez-Martinez<EMAIL_ADDRESS>
Alex Yakubovich<EMAIL_ADDRESS>
Martin Seneviratne<EMAIL_ADDRESS>
Adam D. Lelkes<EMAIL_ADDRESS>
Akshit Tyagi<EMAIL_ADDRESS>
Jonas Kemp<EMAIL_ADDRESS>
Google Research Corresponding author. Mountain View  CA  USA  Ethan Steinberg
<EMAIL_ADDRESS>
N. Lance Downing<EMAIL_ADDRESS>
Ron C. Li<EMAIL_ADDRESS>
Keith E. Morse<EMAIL_ADDRESS>
Nigam H. Shah<EMAIL_ADDRESS>
Department of Medicine  Stanford University  Stanford  CA  USA  Ming-Jun Chen
<EMAIL_ADDRESS>
Google Research  Mountain View  CA  USA

###### Abstract

While it has been well known in the ML community that deep learning models
suffer from instability, the consequences for healthcare deployments are under
characterised. We study the stability of different model architectures trained
on electronic health records, using a set of outpatient prediction tasks as a
case study.

We show that repeated training runs of the same deep learning model on the
same training data can result in significantly different outcomes at a patient
level even though global performance metrics remain stable.

We propose two stability metrics for measuring the effect of randomness of
model training, as well as mitigation strategies for improving model
stability.

###### keywords:

Model Stability, Nondeterminism, Neural Networks, Model Ensembles, Electronic
Health Records

##  1 Introduction

Machine learning methods have shown great promise in a wide range of
healthcare applications, from diagnosing conditions from imaging, to
personalizing care, to improving administrative workflows (see e.g. Davenport
and Kalakota (2019)). One particular use case that has attracted a lot of
research interest is _risk stratification_ ; i.e., assigning each patient a
risk status using a machine learning model trained to predict a patient’s
individual risk of some adverse event. Previously studied risk prediction
tasks cover a range of inpatient and outpatient risks, from acute inpatient
predictions such as acute kidney injury (Mohamadlou et al., 2018) to longer
term outpatient predictions such as readmissions (Huang et al., 2021).

These risk stratification use cases are typically motivated by the promise of
improving overall patient care by helping clinicians prioritize the highest
risk patients. Therefore, they can be used to inform the allocation of scarce
and valuable resources, such as a care manager’s time, among a group of
patients. Whenever algorithmic tools are used for resource allocation,
especially in such an important and sensitive area as healthcare, ensuring the
robustness and fairness of the algorithmic decision-making process is of
utmost importance.

In this paper, we want to draw the community’s attention to a particular
source of concern in this respect: model instability. As in many other areas
of machine learning applications, deep learning methods often demonstrate
state-of-the-art performance for risk prediction tasks. It is well known in
the machine learning community that deep learning models are more sensitive to
random phenomena during model training, and their outputs are often less
stable than those of traditional, simpler models. However, the practical
implications of model instability for risk stratification haven’t been well
studied.

As a case study for this phenomenon, we train feed-forward neural network
model on EHR data to predict all-cause emergency admissions, along with a set
of disease-specific admission events. These models are trained an data from
the EHR of an academic medical center. Importantly, we evaluate the impact of
instability in a pragmatic fashion, by simulating the prospective deployment
setting of the model—namely, as a method to select high-risk patients for a
care manager to contact. We compare deep learning (DL) models against simpler
logistic regression (LR) models on both performance and stability metrics. We
conclude by proposing an ensembling approach to mitigate against DL
instability and suggest directions for future methodological work.

While we used unplanned hospital admission predictions as our case study, the
concept of model stability and our methods to evaluate it are generalizable to
any other models where the ordinal nature of the risk scores is
important—essentially any clinical workflow where resources (e.g. capacity to
intervene or follow up) are restricted.

##  2 Related work

Risk stratification models have been built for a number of prediction tasks,
such as forecasting mortality (Luo et al., 2022), patient deterioration (Mann
et al., 2021), unplanned hospital admissions (King et al., 2022), or the need
for medical procedures (Lopez-Martinez et al., 2022). However, the performance
of these models is usually evaluated using traditional metrics of predictive
performance, such as the area under the curves (AUC) of the receiver operating
characteristic (ROC-AUC) and the precision-recall curve (PR-AUC). Although
having merit, recent work has highlighted the need for additional metrics to
evaluate models.

For example, Mišić et al. (2021) has noted that traditional metrics lack
interpretability and provide an incomplete picture of the overall benefit of
the model as they do not take into account how providers may use them. For
example, if used to allocate scarce medical resources, resource constraints
should inform model evaluation. King et al. (2022) has also noted that
building a model for real-world implementation involves several additional
challenges to those encountered when simply optimising a model’s technical
performance, such as evaluating the drift in model performance over time in
deployment simulations. Model fairness (Rajkomar et al., 2018a; Mhasawade et
al., 2021) should also be key aspect of model evaluation, to avoid disparities
in resource allocation. Therefore, if risk-stratification models are to be
used operationally, model evaluation should go beyond the traditional metrics.

Reproducibility is also a key aspect of the assessment of a model, and a
minimal prerequisite for any machine learning model applied to human health
(McKinney et al. (2020); Beam et al. (2020); Liu et al. (2022)). Lack of
reproducibility undermines the credibility of a prediction, and diminishes
user trust in the machine learning model. Moreover, regulatory agencies are
implementing rules governing medical use of machine learning (Vollmer et al.
(2020)), and may reasonably expect being able to duplicate the results of
model’s evaluation using the same elements as were used in the prior study.

Unfortunately, even models trained with the same training algorithm and
dataset can have different model parameters and yield different predictions.
This is due to nondeterminism during training (Nagarajan and Warnell, ;
Summers and Dinneen, 2021; Madhyastha and Jain, 2019), which includes:

  * •

Parameter initialization Prior to training, the weights of a NN are
initialized randomly (often from a Gaussian or uniform distribution centered
at zero) (Glorot and Bengio, 2010).

  * •

Stochastic optimization NNs are trained with stochastic gradient descent (SGD)
or other similar stochastic optimization techniques, which typically sample
small batches of examples from the training set at random in each iteration.

  * •

Stochastic regularization layers Regularization techniques such as dropout or
noisy activation functions (Gulcehre et al., 2016).

In addition to randomness due to stochastic training algorithms, the hardware
and software tooling used for model training can also result in differences
across models trained with the same data and algorithm. These often trade
determinisim and precision for speed, introducing nondeterminism e.g. by the
random completion order of parallel tasks (Zhuang et al., 2022).

In the EHR literature, model stability issues have already been identified.
D’Amour et al. (2020) and Dusenberry et al. (2020) have noted that changing
random seeds in training led to substantially different test set predictions.
Therefore, Dusenberry et al. (2020) highlighted that model uncertainty is a
key factor one should pay attention to while adopting a DL approach to predict
events from EHR data. Despite these insights, the instability of deep learning
models with EHR data, and specifically its impact in the real-world use of
risk stratification models, remains poorly understood. For example, DL
training instability has been shown to lead to variability in model fairness
(Qian et al., 2021), but this has not been studied in clinical models.

##  3 Methods

###  3.1 Dataset

We used a limited dataset provided by the STAnford medicine Research data
Repository (STARR-OMOP) (Datta et al., 2020), which consists of Stanford’s EHR
data in the Observational Medical Outcomes Partnership (OMOP) Common Data
Model (CDM). Data from Lucile Packard Children’s Hospital (LPCH), which is a
pediatric institution, was excluded from this study. The dataset was utilized
in accordance with HIPAA and Stanford Medicine Research IT processes,
including IRB review.

To leverage existing FHIR-based modeling infrastructure, we manually mapped
the dataset from the OMOP CDM to the HL7 Fast Healthcare Interoperability
Resources (FHIR) standard (version R4) (Bender and Sartipi (2013)).

###  3.2 Model development

To explore the stability of risk stratification methods, we implemented 4
outpatient deterioration tasks as a case study.

####  3.2.1 Task definition

We predicted the 30-day risk of all-cause emergency hospitalizations and three
condition-specific potentially preventable unplanned hospitalizations: chronic
obstructive pulmonary disease (COPD), diabetes and heart failure (HF).

Emergency admissions were identified by examining the class and admit source
of the medical encounter. The 3 disease-specific unplanned hospitalizations
were defined using the inclusion and exclusion criteria of the Agency for
Healthcare Research and Quality’s Patient Quality Indicators (PQIs) for
ambulatory care-sensitive conditions (Agency for Healthcare Research and
Quality, 2016, 2018), and identified by examining the diagnosis and procedure
codes corresponding to the PQI definition (Quan et al., 2005).

####  3.2.2 Label extraction

Following Rajkomar et al. (2018b), we represented each patient’s history as an
ordered timeline of FHIR resources. As in Lopez-Martinez et al. (2022), for
each patient we generated a prediction on the first day of every month within
the study period, if and only if the patient had at least one FHIR resource
within the preceding year. Each example in our dataset corresponded to one
such monthly prediction: we included all resources occurring prior to the
prediction time as features, and we computed a binary label denoting whether
the patient had an emergency admission within the first 30 days following the
prediction time.

####  3.2.3 Feature extraction

For each patient, we extracted information about the demographics, conditions,
procedures, medications, observations and encounters.

The demographics features were gender and age. The condition-based features
contained all International Classification of Diseases (ICD) codes associated
with each encounter. The procedure-based features included all Current
Procedural Terminology (CPT) codes, as well as ICD procedure codes. The
medication-based features contained RxNorm codes (Liu et al., 2005) for
medication requests. The observation-based features included indicators for
measurements of patients’ vital signs and lab results and were constructed by
concatenating the observation code (LOINC (McDonald et al., 2003) and SNOMED
(Cornet and de Keizer, 2008) codes) and the unit of measurement. The
encounter-based feature was the class of encounter.

Following Lopez-Martinez et al. (2022), we bucketed the feature values into
one of four disjoint time intervals indicating how long before the prediction
time that feature value occurred: less than 30 days before, 30-90 days before,
90 days-1 year before, or 1-10 years before. The demographic features (age and
gender) were not bucketed. Each feature was represented as a multi-hot vector.

####  3.2.4 Model architecture

Using the aforementioned features, we designed independent logistic regression
(LR) and feedforward neural network (FNN) machine learning models for each of
the 4 tasks. The FNN consisted of three hidden layers of 128, 64, and 32 units
respectively, with ReLU activations and dropout after each hidden layer,
followed by a final prediction layer. Weights for the FNN model were randomly
initialized using Glorot uniform initialization (Glorot and Bengio, 2010) on
each training run; weights for the LR model were simply initialized to zero.
All models were trained using Adam (Kingma and Ba, 2014), a stochastic
gradient-based optimizer method, with batch size 512. Predictions were
rendered on the 1st day of every month on those patients who had at least one
health interaction within the preceding year and were aged 18 years or older.

All models were implemented and trained with TensorFlow in Python Abadi et al.
(2016). Model hyperparameters, including L1/L2 regularization coefficients,
dropout rates, and the parameters of an exponential learning rate decay, were
tuned using Google Vizier Golovin et al. (2017) to minimize loss on the
validation set.

####  3.2.5 Model calibration

While it is desirable for our models to output well-calibrated risk
probabilities, our LR and FNN models, like most supervised learning methods,
produce scores that can be used to rank the predictions from highest to lowest
risk, but are not accurate probability estimates.

To address this limitation, we used isotonic regression (Zadrozny and Elkan,
2002; Niculescu-Mizil and Caruana, 2005) to calibrate these models and
transform ranking scores into probability estimates. Specifically, we fit a
non-parametric isotonic (monotonically increasing) step-wise function $m$ to
the original model scores $f_{i}$ such that $y_{i}=m(f_{i})+\epsilon_{i}$,
where $y_{i}$ denotes the true targets. Isotonic regression works by finding
the isotonic function $\hat{m}$ such that:

$\hat{m}=\arg\min_{z}\sum(y_{i}-z(f_{i}))^{2}$ (1)

subject to $z(f_{i})\geq z(f_{j})$ whenever $f_{i}\geq f_{j}$.

All risk predictions reported in this papers are calibrated probabilities
obtained from after applying isotonic regression.

###  3.3 Model evaluation

Model evaluation was done on an fixed independent held out test set comprising
a 10% random selection of the available patients. No patients in the test set
were used for model training or validation.

Performance was evaluated using two approaches. In the first one, we
aggregated all test set predictions to compute the ROC-AUC and PR-AUC. Note
that given the low prevalence of triggers with positive labels, the PR-AUC
metric is more informative than ROC-AUC (Saito and Rehmsmeier, 2015). In the
second approach, we simulated the real-word deployment of the model by
evaluating the test set predictions of each month independently (from January
2009 to December 2021), and then averaging these monthly metrics. This last
evaluation approach allows us to introduce capacity constraints by examining
the top K highest risk patients each month only.

###  3.4 Stability analysis

To quantify model stability, in addition to the standard deviation of the ROC-
AUC and PR-AUC metrics, we implemented the following metrics, which are
computed using $N=10$ model runs with exactly the same training data,
configurations and model hyperparameters, but different random seeds.

####  3.4.1 Top K Jaccard index

In risk stratification methods, users may care about the $K$ highest-risk
patients, as they are often limited in the amount of patients they can act
upon. Therefore, it is reasonable to expect that the list of top $K$ patients
should be the same across different models.

To assess the similarity of the $N$ patient lists, we use the pairwise Jaccard
index which calculates the average Jaccard index across all pairwise
comparisons. The Jaccard index is given by

$j({X^{K}_{i},X^{K}_{j}})=\frac{|X^{K}_{i}\cap X^{K}_{j}|}{|X^{K}_{i}\cup
X^{K}_{j}|}$ (2)

where $X^{K}_{i}$ and $X^{K}_{j}$ denote the top $K$ predictions in models $i$
and $j$ respectively.

For $N$ models, there exist ${N}\choose{2}$ pairwise comparisons. Therefore,
the top $K$ pairwise Jaccard index is defined as follows:

$J=\frac{1}{{{N}\choose{2}}}\sum_{i<j}j({X^{K}_{i},X^{K}_{j}})$ (3)

The choice of $K$ depends on the specific use case.

####  3.4.2 Rank correlation coefficient

Risk scores are often used to prioritize patients, so that those with higher
risk scores are acted upon earlier. In these scenarios, the ordinal nature of
the patient lists is important. Users would expect models to produce similarly
ranked patient lists.

Kendall’s $\tau$ is commonly used to compare the ordinal correlation between
two ranked lists (Kendall, 1938, 1945).

It is based on computing the number of concordant ($|C|$) and discordant
($|D|$) pairs:

$\tau=\frac{|C|-|D|}{{{N}\choose{2}}}$ (4)

where a pair $((x_{i},y_{i}),(x_{j},y_{j}))$ is defined to be concordant if
the sort order of the $x$’s is the same as the sort order of the $y$’s,
otherwise it is discordant.

We also consider a weighted version of Kendall’s tau with hyperbolic
weighting, in which correlations between top ranked patients have more
influence on the metric (Vigna (2015)).

In this formulation, a datapoint of rank $r$ is assigned a weight of
${(r+1)}^{-1}$. The weight of exchanging two datapoints with ranks $r$ and $s$
is the sum of their weights, i.e. ${(r+1)}^{-1}+{(s+1)}^{-1}$.

Weighting captures the assumption that the relative ordering of very low-risk
patients is less important since interventions are typically aimed at high-
risk patients.

###  3.5 Ensembling

Ensembling is a popular approach for improving prediction performance across
many different tasks, and it also reduces run-to-run variablity. Averaging the
outputs of multiple models reduces the random noise which might be causing the
unnecessary variance in the predictive distribution (Huang et al., 2017; Xie
et al., 2013). It has also been shown to give better-calibrated estimates for
uncertainty in the predictive distributions (Lakshminarayanan et al., 2016).
In our experiments we use a simple arithmetic mean of model outputs as a
baseline; more complex averaging schemes could potentially yield further
improvement.

##  4 Results

###  4.1 Dataset

Task | # predictions | # positives | ROC-AUC | PR-AUC
---|---|---|---|---
LR | FNN | LR | FNN
Emergency | 60,259,528 | 104,123 (0.17%) | 0.9052 | 0.9116 | 0.0734 | 0.0858
COPD | 5,759,224 | 1,228 (0.02%) | 0.8894 | 0.9234 | 0.0047 | 0.0069
Diabetes | 5,034,065 | 701 (0.01%) | 0.8273 | 0.8919 | 0.0177 | 0.0057
Heart failure | 2,161,112 | 2,369 (0.11%) | 0.8717 | 0.8789 | 0.0160 | 0.0170
Table 1: Number of predictions and positive label prevalence in the
development (train and validation) set, together with model performance on the
aggregated test set predictions for single ($N=1$) LR and FNN models. For each
metric, the best performing model is highlighted in bold. ROC: receiver
operating characteristic; PR: precision-recall curve; AUC: area under the
curve.

A total of 1,576,215 individuals were present in the dataset. Of these, 53.66%
were females. White was the most predominant race (44%), followed by asian
(18.13%), black (3.85%), hispanic (2.02%) and native american (1.12%).

From these patients, we extracted a total of 66,977,919, 6,412,742, 5,588,138,
and 2,400,145 predictions for the emergency, COPD, diabetes and heart failure
tasks respectively. These predictions spanned a time period from January 1st
2009 to December 1st 2021. The differences in the number of predictions for
each task can be explained by the different cohort inclusion criteria for
prediction rendering.

Using this data, we generated training, validation and test sets by randomly
selecting 80%, 10% and 10% of patients respectively. Positive label prevalence
in the development set (training and validation sets) for the 4 tasks is shown
in Table 1.

###  4.2 Model performance

Table 1 shows the performance of the LR and FNN models for all 4 tasks in
terms of ROC-AUC and PR-AUC. Note that given the imbalance of the data
(positive labels were overwhelmingly underrepresented), PR-AUC is more
informative (Saito and Rehmsmeier (2015)).

The best performance was obtained by the FNN model architecture, with the
exception of the diabetes task, where it was outperformed by the linear model.
In both cases, the emergency task achieved the best performance, and did
significantly better than the 3 remaining tasks. This may be explained by the
significantly higher label prevalence of the emergency task.

###  4.3 Stability analysis

In order to conduct the model stability analysis described in Sec.3.4, for
each task and model architecture we trained $N=10$ models using the same model
hyperparameters and training data. Differences in model parameters and
predictions across model runs may only be explained by nondeterminism in the
training algorithm and hardware. Note that the choice of $N$ was based on
visual inspection of the convergence of the stability metrics as we increased
$N$. We determined that increasing it beyond $N=10$ would not have provided
any additional value.

First, we computed the standard deviation of the ROC-AUC and PR-AUC metrics
across the $N$ model runs. This is depicted in Fig.3, which shows higher
variability for the FNN architecture, with the exception of the PR-AUC metric
of the diabetes model. Also, note that all PR-AUC standard deviation metrics
are smaller than $10^{-2}$, and these are difficult to contextualize and
interpret. Hence, the standard deviations of the ROC-AUC and PR-AUC metrics
may be inadequate for the assessment of model stability.

Figure 1: Standard deviation of the ROC-AUC and PR-AUC metrics.

Figure 2: Average pairwise Jaccard index $J(K)$ for all models.

Figure 3: Comparison of calibrated risks produced by two FNN models.

A better assessment of model stability may be provided by the pairwise Jaccard
index $J(K)$, which is shown in Fig.3 for different values of $K$. By
illustrating average differences in top $K$ patient lists, it provides a more
interpretable stability assessment. Of note is that the linear architecture
achieves the best stability when compared to the FNN one, for all tasks and
values of $K$. In addition to this, note that $J(K)$ decreases inversely with
$K$. This suggests that stability is less of an issue for large values of $K$,
but when resources are constrained and $K$ represents a small subset of the
patient population, patient lists may be impacted by instability.

Note that given that the test set is comprised of a random sample of 10% of
the patients available in the STARR-OMOP dataset, a top K value of e.g.
$K=100$ in the test set would correspond to 1,000 monthly patients in the
entire dataset.

A different insight into model stability is provided by Kendall’s $\tau$ rank
correlation coefficient, which, unlike the Jaccard index, takes into account
the relative ordering of the patient risk scores. We consider both the
standard Kendall’s $\tau$ and one with hyberbolic weighting, where low-risk
patients are downweighted (Vigna (2015)). We observed that the FNN model
resulted in lower values of both weighted and unweighted $\tau$ with the
exception of the unweighted $\tau$ of the COPD task, which confirms the
previous insights that the LR model is more stable. In addition to this, we
noted that using hyperbolic weighting resulted in larger values of $\tau$,
suggesting that ordinal relationships are more preserved for the highest
ranked patients.

Task | Model |  Unweighted $\tau$ |  Weighted $\tau$
---|---|---|---
Emergency | LR | 0.9647 (0.0027) | 0.9912 (0.0029)
FNN | 0.9314 (0.0026) | 0.9817 (0.0032)
EFNN | 0.9796 (0.0017) | 0.9947 (0.0008)
COPD | LR | 0.9267 (0.0214) | 0.9864 (0.0049)
FNN | 0.9540 (0.0013) | 0.9831 (0.0020)
Diabetes | LR | 0.9060 (0.0188) | 0.9797 (0.0037)
FNN | 0.5883 (0.0700) | 0.8312 (0.0704)
HF | LR | 0.9870 (0.0010) | 0.9950 (0.0007)
FNN | 0.8847 (0.0197) | 0.9704 (0.0038)
Table 2: Average pairwise Kendall’s rank correlation coefficient ($\tau$), and
standard deviations, with uniform and hyperbolic weighting.

Finally, we explore the differences in calibrated risk scores. Fig.3 compares
two randomly selected models for each task. Of note is that an individual’s
calibrated risk score can vary significantly due to training variability. For
the emergency model, across all pairwise comparisons we observed that the
highest ranked patients experienced the largest variability (see Fig.3, top
left, for an example). This behavior was occasionally present in the disease-
specific tasks, but the risks of the mid-ranked patients varied more more as
illustrated in the comparison of the COPD and diabetes tasks.

###  4.4 Ensemble

We trained $N=10$ averaging ensembles, each one comprised of $M=10$ models
whose final outputs are averaged. Due to the cost of training $N\times M$
models, we focused our analysis on the FNN emergency model only.

Ensembling did not significantly impact the average ROC-AUC, which was 0.9111
for both the FNN and ensembled FNN (E-FNN), but did improve the average PR-AUC
from 0.0843 to 0.0848, and the variance of these AUC metrics, as shown in
Fig.3. Kendall’s $\tau$ (Table 2) and the Jaccard index ( Fig.5) also improved
significantly, surpassing the linear architecture in stability.

Figure 4: Average pairwise Jaccard index $J(K)$ for the LR, FNN and E-FNN
emergency models. The shadow represents one standard deviation computed over
all prediction times.

Figure 5: Range in gender (left) and race (right) composition of the top K
patients for the $N$ runs of the FNN (solid lines) and E-FNN (dotted lines)
emergency models.

###  4.5 Fairness analysis

Finally, we investigate whether the stability issues discussed in Sec.4.3
impact the fairness of the model, specifically, whether it differentially
impacted the gender and race representation in the monthly lists of top K
highest-risk patients.

For each patient and prediction time, we computed the range in representation
across $N=10$ models. This is illustrated in Fig.5 for the FNN emergency model
and different values of $K$. It shows that in addition to leading to different
patient selections and individual risk scores, instability can also lead to
gender and race differences in patient selection. Note that while smaller
differences were observed for the Asian and Black race categories, that is
explained by their reduced prevalence in the dataset (see Sec.4.1).

Finally, we also noted that the E-FNN reduced both gender and race
variability. For example, for $K=30$, the range in gender across the $N$ FNN
model runs was 10.04%, whereas the range across the $N$ E-FNN model runs
(where each E-FNN model was derived from $M=10$ independent FNN models) was
4.94%. Similarly, the range in race composition decreased with the ensemble
from 9.92% to 4.88%, 6.02% to 2.44%, and 4.42% to 1.67% for white, asian and
black races respectively.

##  5 Discussion

This work examined the instability of risk stratification models resulting
from nondeterminism in the training process. Using a set of unplanned hospital
admission prediction tasks as a case study, we confirmed previous observations
that repeated training runs can lead to significantly different predictions at
a patient level even if traditional model accuracy metrics remain unchanged.

We have proposed two similarity metrics to quantify this run-to-run
variability among risk stratification models in a way that accounts for real-
world usage under resource constraints. The first one, the Jaccard index,
evaluates the consistency in patients selected for intervention when the model
is used to allocate limited resources. The second one, Kendall’s rank
correlation coefficient, evaluates whether the ordinal relationships are
preserved; therefore, it is broadly useful to evaluate models that prioritize
patients for intervention. These stability metrics can also be used to select
the most stable models.

We have also shown that the observed instability can lead to variance in the
representation of subgroups among the top K predictions, which may have
implications for equitable allocation of scare resource and could in turn
diminish trust in the models.

Finally, we have evaluated a simple strategy (model ensembling) to mitigate
model instability. There are other potential strategies, such as implementing
techniques for addressing data imbalance (e.g. synthetic minority oversampling
or weighted loss functions), decaying the learning rate, replacing the ReLU
activation function with SmeLU (Shamir et al., 2020) etc., that can also
improve the model stability. However, producing a throughout review of
strategies to improve stability is not in the scope of this work; instead, we
want to provide an example to show how stability issue can be mitigated. Our
averaging ensemble improved all stability metrics considered while preserving
overall model performance. However, this benefit of creating an ensemble comes
with an extra cost of multiple model runs in memory and disk space. Future
work should consider more efficient ensemble approaches, such as batch
ensembling (Wen et al., 2020).

The acceptable level of instability will vary among different application use
cases depending on a complex set of factors such as capacity to intervene,
subgroup representation, and frequency of training runs; we do not make any
claim as to whether any specific models we have examined are stable enough for
particular applications. However, we believe that model instability has been
an overlooked factor that should receive more consideration when deciding
whether to deploy a model to a healthcare setting. In order for potential
users to be able to make informed decisions, researchers need to transparently
report and discuss the stability of their proposed methods. We hope that our
work will lead to more transparency and that it will start a discussion about
the importance of model stability in the research community.

We are grateful to Lan Huong Nguyen, Andreas Schagerer, Eric Loreaux and
Bakhtiar Rehman for their contributions to make this work possible.

## References

  * Abadi et al. (2016) Martin Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, Manjunath Kudlur, Josh Levenberg, Rajat Monga, Sherry Moore, Derek G Murray, Benoit Steiner, Paul Tucker, Vijay Vasudevan, Pete Warden, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: A system for large-scale machine learning.  In _12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16)_ , pages 265–283, 2016.
  * Agency for Healthcare Research and Quality (2016) Agency for Healthcare Research and Quality.  Prevention quality indicators technical specifications - version 6.0, september 2016.  https://www.qualityindicators.ahrq.gov/modules/pqi_techspec.aspx, 2016.
  * Agency for Healthcare Research and Quality (2018) Agency for Healthcare Research and Quality.  Prevention quality indicators technical specifications updates - version v2018 and v2018.0.1 (ICD 10-CM/PCS).  https://www.qualityindicators.ahrq.gov/Modules/PQI_TechSpec_ICD10_v2018.aspx, 2018\.
  * Beam et al. (2020) Andrew L Beam, Arjun K Manrai, and Marzyeh Ghassemi.  Challenges to the reproducibility of machine learning models in health care.  _JAMA_ , 323(4):305–306, January 2020.
  * Bender and Sartipi (2013) Duane Bender and Kamran Sartipi.  HL7 FHIR: An agile and RESTful approach to healthcare information exchange, 2013.
  * Cornet and de Keizer (2008) Ronald Cornet and Nicolette de Keizer.  Forty years of SNOMED: a literature review.  _BMC Med. Inform. Decis. Mak._ , 8 Suppl 1:S2, October 2008\.
  * D’Amour et al. (2020) Alexander D’Amour, Katherine Heller, Dan Moldovan, Ben Adlam, Babak Alipanahi, Alex Beutel, Christina Chen, Jonathan Deaton, Jacob Eisenstein, Matthew D Hoffman, Farhad Hormozdiari, Neil Houlsby, Shaobo Hou, Ghassen Jerfel, Alan Karthikesalingam, Mario Lucic, Yian Ma, Cory McLean, Diana Mincu, Akinori Mitani, Andrea Montanari, Zachary Nado, Vivek Natarajan, Christopher Nielson, Thomas F Osborne, Rajiv Raman, Kim Ramasamy, Rory Sayres, Jessica Schrouff, Martin Seneviratne, Shannon Sequeira, Harini Suresh, Victor Veitch, Max Vladymyrov, Xuezhi Wang, Kellie Webster, Steve Yadlowsky, Taedong Yun, Xiaohua Zhai, and D Sculley.  Underspecification presents challenges for credibility in modern machine learning.  _J. Mach. Learn. Res._ , 2020.
  * Datta et al. (2020) Somalee Datta, Jose Posada, Garrick Olson, Wencheng Li, Ciaran O’Reilly, Deepa Balraj, Joseph Mesterhazy, Joseph Pallas, Priyamvada Desai, and Nigam Shah.  A new paradigm for accelerating clinical data science at stanford medicine.  March 2020.
  * Davenport and Kalakota (2019) Thomas Davenport and Ravi Kalakota.  The potential for artificial intelligence in healthcare.  _Future Healthc J_ , 6(2):94–98, June 2019.
  * Dusenberry et al. (2020) Michael W Dusenberry, Dustin Tran, Edward Choi, Jonas Kemp, Jeremy Nixon, Ghassen Jerfel, Katherine Heller, and Andrew M Dai.  Analyzing the role of model uncertainty for electronic health records.  In _Proceedings of the ACM Conference on Health, Inference, and Learning_ , CHIL ’20, pages 204–213, New York, NY, USA, April 2020. Association for Computing Machinery.
  * Glorot and Bengio (2010) Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, PMLR_ , 9:249–256, 2010.
  * Golovin et al. (2017) Daniel Golovin, Benjamin Solnik, Subhodeep Moitra, Greg Kochanski, John Karro, and D Sculley.  Google vizier: A service for Black-Box optimization.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’17, pages 1487–1495, New York, NY, USA, August 2017. Association for Computing Machinery.
  * Gulcehre et al. (2016) Caglar Gulcehre, Marcin Moczulski, Misha Denil, and Yoshua Bengio.  Noisy activation functions.  March 2016.
  * Huang et al. (2017) Gao Huang, Yixuan Li, Geoff Pleiss, Zhuang Liu, John E Hopcroft, and Kilian Q Weinberger.  Snapshot ensembles: Train 1, get M for free.  April 2017.
  * Huang et al. (2021) Yinan Huang, Ashna Talwar, Satabdi Chatterjee, and Rajender R Aparasu.  Application of machine learning in predicting hospital readmissions: a scoping review of the literature.  _BMC Med. Res. Methodol._ , 21(1):96, May 2021\.
  * Kendall (1938) M G Kendall.  A new measure of rank correlation.  _Biometrika_ , 30(1/2):81–93, 1938.
  * Kendall (1945) M G Kendall.  The treatment of ties in ranking problems.  _Biometrika_ , 33(3):239, November 1945.
  * King et al. (2022) Zella King, Joseph Farrington, Martin Utley, Enoch Kung, Samer Elkhodair, Steve Harris, Richard Sekula, Jonathan Gillham, Kezhi Li, and Sonya Crowe.  Machine learning for real-time aggregated prediction of hospital admission for emergency patients.  _NPJ Digit Med_ , 5(1):104, July 2022.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  December 2014.
  * Lakshminarayanan et al. (2016) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  December 2016.
  * Liu et al. (2005) S Liu, Wei Ma, R Moore, V Ganesan, and S Nelson.  RxNorm: prescription for electronic drug information exchange.  _IT Prof._ , 7(5):17–23, September 2005.
  * Liu et al. (2022) Zhichao Liu, Ting Li, Skylar Connor, Shraddha Thakkar, Ruth Roberts, and Weida Tong.  Best practice and reproducible science are required to advance artificial intelligence in real-world applications.  _Brief. Bioinform._ , 23(4), July 2022.
  * Lopez-Martinez et al. (2022) Daniel Lopez-Martinez, Christina Chen, and Ming-Jun Chen.  Machine learning for dynamically predicting the onset of renal replacement therapy in chronic kidney disease patients using claims data.  In _Applications of Medical Artificial Intelligence_ , pages 18–28. Springer Nature Switzerland, 2022.
  * Luo et al. (2022) Cida Luo, Yi Zhu, Zhou Zhu, Ranxi Li, Guoqin Chen, and Zhang Wang.  A machine learning-based risk stratification tool for in-hospital mortality of intensive care unit patients with heart failure.  _J. Transl. Med._ , 20(1):136, March 2022.
  * Madhyastha and Jain (2019) Pranava Madhyastha and Rishabh Jain.  On model stability as a function of random seed.  In _The SIGNLL Conference on Computational Natural Language Learning_ , November 2019.
  * Mann et al. (2021) Kay D Mann, Norm M Good, Farhad Fatehi, Sankalp Khanna, Victoria Campbell, Roger Conway, Clair Sullivan, Andrew Staib, Christopher Joyce, and David Cook.  Predicting patient deterioration: A review of tools in the digital hospital setting.  _J. Med. Internet Res._ , 23(9):e28209, September 2021.
  * McDonald et al. (2003) Clement J McDonald, Stanley M Huff, Jeffrey G Suico, Gilbert Hill, Dennis Leavelle, Raymond Aller, Arden Forrey, Kathy Mercer, Georges DeMoor, John Hook, Warren Williams, James Case, and Pat Maloney.  LOINC, a universal standard for identifying laboratory observations: A 5-year update.  _Clin. Chem._ , 49(4):624–633, April 2003.
  * McKinney et al. (2020) Scott Mayer McKinney, Alan Karthikesalingam, Daniel Tse, Christopher J Kelly, Yun Liu, Greg S Corrado, and Shravya Shetty.  Reply to: Transparency and reproducibility in artificial intelligence, 2020.
  * Mhasawade et al. (2021) Vishwali Mhasawade, Yuan Zhao, and Rumi Chunara.  Machine learning and algorithmic fairness in public and population health.  _Nature Machine Intelligence_ , 3(8):659–666, July 2021.
  * Mišić et al. (2021) Velibor V Mišić, Kumar Rajaram, and Eilon Gabel.  A simulation-based evaluation of machine learning models for clinical decision support: application and analysis using hospital readmission.  _NPJ Digit Med_ , 4(1):98, June 2021.
  * Mohamadlou et al. (2018) Hamid Mohamadlou, Anna Lynn-Palevsky, Christopher Barton, Uli Chettipally, Lisa Shieh, Jacob Calvert, Nicholas R Saber, and Ritankar Das.  Prediction of acute kidney injury with a machine learning algorithm using electronic health record data.  _Can J Kidney Health Dis_ , 5:2054358118776326, June 2018\.
  * (32) Prabhat Nagarajan and Garrett Warnell.  The impact of nondeterminism on reproducibility in deep reinforcement learning.
  * Niculescu-Mizil and Caruana (2005) Alexandru Niculescu-Mizil and Rich Caruana.  Predicting good probabilities with supervised learning, 2005.
  * Qian et al. (2021) Shangshu Qian, Viet Hung Pham, Thibaud Lutellier, Zeou Hu, Jungwon Kim, Lin Tan, Yaoliang Yu, Jiahao Chen, and Sameena Shah.  Are my deep learning systems fair? an empirical study of Fixed-Seed training.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 30211–30227, 2021.
  * Quan et al. (2005) Hude Quan, Vijaya Sundararajan, Patricia Halfon, Andrew Fong, Bernard Burnand, Jean-Christophe Luthi, L Duncan Saunders, Cynthia A Beck, Thomas E Feasby, and William A Ghali.  Coding algorithms for defining comorbidities in ICD-9-CM and ICD-10 administrative data.  _Med. Care_ , 43(11):1130–1139, November 2005\.
  * Rajkomar et al. (2018a) Alvin Rajkomar, Michaela Hardt, Michael D Howell, Greg Corrado, and Marshall H Chin.  Ensuring fairness in machine learning to advance health equity.  _Ann. Intern. Med._ , 169(12):866–872, December 2018a.
  * Rajkomar et al. (2018b) Alvin Rajkomar, Eyal Oren, Kai Chen, Andrew M Dai, Nissan Hajaj, Michaela Hardt, Peter J Liu, Xiaobing Liu, Jake Marcus, Mimi Sun, Patrik Sundberg, Hector Yee, Kun Zhang, Yi Zhang, Gerardo Flores, Gavin E Duggan, Jamie Irvine, Quoc Le, Kurt Litsch, Alexander Mossin, Justin Tansuwan, De Wang, James Wexler, Jimbo Wilson, Dana Ludwig, Samuel L Volchenboum, Katherine Chou, Michael Pearson, Srinivasan Madabushi, Nigam H Shah, Atul J Butte, Michael D Howell, Claire Cui, Greg S Corrado, and Jeffrey Dean.  Scalable and accurate deep learning with electronic health records.  _NPJ Digit Med_ , 1:18, May 2018b.
  * Saito and Rehmsmeier (2015) Takaya Saito and Marc Rehmsmeier.  The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets.  _PLoS One_ , 10(3):e0118432, March 2015.
  * Shamir et al. (2020) Gil I. Shamir, Dong Lin, and Lorenzo Coviello.  Smooth activations and reproducibility in deep networks.  2020\.  10.48550/ARXIV.2010.09931.  URL https://arxiv.org/abs/2010.09931.
  * Summers and Dinneen (2021) Cecilia Summers and Michael J Dinneen.  Nondeterminism and instability in neural network optimization.  In _International Conference on Machine Learning (ICML)_ , July 2021\.
  * Vigna (2015) Sebastiano Vigna.  A weighted correlation index for rankings with ties.  In _Proceedings of the 24th International Conference on World Wide Web_ , WWW ’15, pages 1166–1176, Republic and Canton of Geneva, CHE, May 2015\. International World Wide Web Conferences Steering Committee.
  * Vollmer et al. (2020) Sebastian Vollmer, Bilal A Mateen, Gergo Bohner, Franz J Király, Rayid Ghani, Pall Jonsson, Sarah Cumbers, Adrian Jonas, Katherine S L McAllister, Puja Myles, David Granger, Mark Birse, Richard Branson, Karel G M Moons, Gary S Collins, John P A Ioannidis, Chris Holmes, and Harry Hemingway.  Machine learning and artificial intelligence research for patient benefit: 20 critical questions on transparency, replicability, ethics, and effectiveness.  _BMJ_ , 368:l6927, March 2020.
  * Wen et al. (2020) Yeming Wen, Dustin Tran, and Jimmy Ba.  BatchEnsemble: An alternative approach to efficient ensemble and lifelong learning.  February 2020.
  * Xie et al. (2013) Jingjing Xie, Bing Xu, and Zhang Chuang.  Horizontal and vertical ensemble with deep representation for classification.  June 2013.
  * Zadrozny and Elkan (2002) Bianca Zadrozny and Charles Elkan.  Transforming classifier scores into accurate multiclass probability estimates, 2002.
  * Zhuang et al. (2022) Donglin Zhuang, Xingyao Zhang, Shuaiwen Leon Song, and Sara Hooker.  Randomness in neural network training: Characterizing the impact of tooling.  In _Proceedings of the 5th MLSys Conference_ , 2022.