# Mind the Gap: A Causal Perspective on Bias Amplification in Prediction &
Decision-Making

Drago Plecko and Elias Bareinboim
Department of Computer Science
Columbia University
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

As society increasingly relies on AI-based tools for decision-making in
socially sensitive domains, investigating fairness and equity of such
automated systems has become a critical field of inquiry. Most of the
literature in fair machine learning focuses on defining and achieving fairness
criteria in the context of prediction, while not explicitly focusing on how
these predictions may be used later on in the pipeline. For instance, if
commonly used criteria, such as independence or sufficiency, are satisfied for
a prediction score $S$ used for binary classification, they need not be
satisfied after an application of a simple thresholding operation on $S$ (as
commonly used in practice). In this paper, we take an important step to
address this issue in numerous statistical and causal notions of fairness. We
introduce the notion of a margin complement, which measures how much a
prediction score $S$ changes due to a thresholding operation. We then
demonstrate that the marginal difference in the optimal 0/1 predictor
$\widehat{Y}$ between groups, written $P(\hat{y}\mid x_{1})-P(\hat{y}\mid
x_{0})$, can be causally decomposed into the influences of $X$ on the
$L_{2}$-optimal prediction score $S$ and the influences of $X$ on the margin
complement $M$, along different causal pathways (direct, indirect, spurious).
We then show that under suitable causal assumptions, the influences of $X$ on
the prediction score $S$ are equal to the influences of $X$ on the true
outcome $Y$. This yields a new decomposition of the disparity in the predictor
$\widehat{Y}$ that allows us to disentangle causal differences inherited from
the true outcome $Y$ that exists in the real world vs. those coming from the
optimization procedure itself. This observation highlights the need for more
regulatory oversight due to the potential for bias amplification, and to
address this issue we introduce new notions of weak and strong business
necessity, together with an algorithm for assessing whether these notions are
satisfied. We apply our method to three real-world datasets and derive new
insights on bias amplification in prediction and decision-making.

##  1 Introduction

Automated systems based on machine learning and artificial intelligence are
increasingly used for decision-making in a variety of real-world settings.
These applications include hiring decisions, university admissions, law
enforcement, credit lending and loan approvals, health care interventions, and
many other high-stakes scenarios in which the automated system may
significantly affect the well-being of individuals (Khandani et al., 2010;
Mahoney and Mohen, 2007; Brennan et al., 2009). In this context, society is
increasingly concerned about the implications and consequences of using
automated systems, compared to the currently implemented decision processes.
Prior works highlight the potential of automated systems to perpetuate or even
amplify inequities between demographic groups, with a range of examples from
decision support systems for (among others) sentencing Angwin et al. (2016),
face-detection Buolamwini and Gebru (2018), online advertising Sweeney (2013);
Datta et al. (2015), and authentication Sanburn (2015). Notably, issues of
unfairness and discrimination are also pervasive in settings in which
decisions are made by humans. Some well-studied examples include the gender
pay gap, supported by a decades-long literature (Blau and Kahn, 1992, 2017),
or the racial bias in criminal sentencing (Sweeney and Haney, 1992; Pager,
2003), just to cite a few. Therefore, AI systems designed to make decisions
may often be trained with data that contains various historical biases and
past discriminatory decisions against certain protected groups, constituting a
large part of the underlying problem. In this work, we specifically focus on
investigating when automated systems may potentially lead to an even more
discriminatory process, possibly amplifying already existing differences
between groups.

Within this context, it is useful to distinguish between different tasks
appearing in the growing literature on fair machine learning. One can
distinguish three specific and different tasks, namely (1) bias detection and
quantification for exisiting outcomes or decision policies; (2) construction
of fair predictions of an outcome; (3) construction of fair decision-making
policies that are intended to be implemented in the real-world. Interestingly,
a large portion of the literature in fair ML focuses on the second task of
fair prediction, and what is often left unaddressed is how these predictions
may be used later on in the pipeline, and what kind of consequences they may
have. For instance, consider a prediction score $S$ for a binary outcome $Y$
that satisfies well-known fairness criteria, such as independence (demographic
parity (Darlington, 1971)) or sufficiency (calibration (Chouldechova, 2017)).
After a simple thresholding operation, commonly applied in settings with a
binary outcome, the resulting predictor is no longer guaranteed to satisfy
independence or sufficiency, and the previously provided fairness guarantees
may be entirely lost. The same behavior can be observed for numerous other
measures.

These difficulties do not apply only to statistical measures of fairness.
Recently, a growing literature has explored causal approaches to fair machine
learning (Kusner et al., 2017; Kilbertus et al., 2017; Nabi and Shpitser,
2018; Zhang and Bareinboim, 2018b, a; Wu et al., 2019; Chiappa, 2019; Plečko
and Meinshausen, 2020; Plečko and Bareinboim, 2024), which have two major
benefits. First, they allow for human-understandable and interpretable
definitions and metrics of fairness, which are tied to the causal mechanisms
transmitting the change between groups. Secondly, they offer a language that
is aligned with the legal notions of discrimination, such as the disparate
impact doctrine. In particular, causal approaches allow for considerations of
business necessity – which aim to ellucidate which covariates may be
justifiably used by decision-makers even if their usage implies a disparity
between groups. However, causal approaches to fairness also suffer from the
above-discussed issues – namely, a guarantee of absence of a causal influence
from the protected attribute $X$ onto a predictor $S$ need not hold true after
the predictor is thresholded (Plečko and Bareinboim, 2024). Therefore, within
the causal approach, there is also a major need for a better understanding of
how probabilistic predictions are translated into binary predictions or
decisions.

In this work, we take an important step in the direction of addressing this
issue. We work in a setting with a binary label $Y$, and the goal is to
provide a binary prediction $\widehat{Y}$ or a binary decision $D$. Our
approach is particularly suitable for settings in which the utility of the
decision is monotonic with respect to the conditional probability of $Y$ being
positive, written $P(Y\mid\text{covariates})$, but the developed tools also
have ramifications more general utilities. Examples that fall under our scope
are numerous; for instance, the utility of admitting a student to the
university ($D$) is often monotonic in the probability that the student
successfully graduates ($Y$). In the context of criminal justice, decisions of
detention ($D$) are used to prevent recidivism, and the utility of the
decision is monotonic in the probability that the individual recidivates
($Y$). Finally, various preventive measures in healthcare ($D$, such as
vaccination, screening tests, etc.) are considered to be best applied to
individuals with the highest risk of developing a target disease or suffering
a negative outcome ($Y$). We now provide an example that illustrates one of
the key insights of our paper:

######  Example 1 (Disparities in Hiring).

Consider a company deciding to hire employees using an automated system for
the first time. From a previous hiring cycle when humans were in charge, the
company has access to data on gender $X$ ($x_{0}$ for female, $x_{1}$ for
male) and the hiring outcome $Y$ ($y_{1}$ for being hired, $y_{0}$ otherwise).
The true underlying mechanisms of the system are given by:

$\displaystyle X$ $\displaystyle\leftarrow\mathbb{1}(U_{X}<0.5)$ (1)
$\displaystyle Y$
$\displaystyle\leftarrow\begin{cases}\mathbb{1}(U_{Y}<p_{0})\text{ if
}X=x_{0}\\\ \mathbb{1}(U_{Y}<p_{1})\text{ if }X=x_{1}\end{cases}$ (2)

where $U_{X},U_{Y}\sim\text{Unif}[0,1]$. In words, an applicant is female with
a 50% probability, and the probability of being hired as a female $x_{0}$ is
$p_{0}$, whereas for males $x_{1}$ the probability is $p_{1}$. The company
finds the optimal prediction score $S$ to be $S(x)=p_{x}$. The optimal 0/1
predictor $\widehat{Y}$, which will also be the company’s decision, is given
by $\widehat{Y}(x)=\mathbbm{1}(S(x)\geq\frac{1}{2})$, meaning that the company
will threshold the predictions at $\frac{1}{2}$. Suppose that
$p_{0}=0.49,p_{1}=0.51$, and consider the visualization in Fig. 1. The
probability $p_{0}=0.49$ means that 49 out of 100 females were hired, while 51
were not. After thresholding,
$\widehat{Y}(x_{0})=\mathbbm{1}(p_{0}\geq\frac{1}{2})=0$ for each female
applicant, meaning that the thresholding operation maps the prediction of each
individual to $0$, even though 49/100 would have a positive outcome (Fig. 1
right). Similarly, for $p_{1}=0.51$, we have that 51/100 male applicants would
be hired, resulting in a thresholded predictor $\widehat{Y}(x_{1})=1$ for all
the applicants, even though 49/100 would not have a positive outcome (Fig. 1
left). Therefore, the gender disparity in hiring after introducing the
automated predictor $\widehat{Y}$ is 100%, compared to a 2% disparity in the
outcome $Y$ before introducing $\widehat{Y}$. Formally, the disparity in the
optimal 0/1 predictor $P(\widehat{y}\mid x_{1})-P(\widehat{y}\mid
x_{0})=\mathbbm{1}(p_{1}\geq\frac{1}{2})-\mathbbm{1}(p_{0}\geq\frac{1}{2})$
can be decomposed as:

$\displaystyle P(\widehat{y}\mid x_{1})-P(\widehat{y}\mid
x_{0})=\underbrace{p_{1}-p_{0}}_{\text{Term
I}}+\underbrace{(\mathbbm{1}(p_{1}\geq\frac{1}{2})-p1)-(\mathbbm{1}(p_{0}\geq\frac{1}{2})-p_{0})}_{\text{Term
II}}.$ (3)

Term I measures the disparity coming from the true outcome $Y$ (2% disparity),
which is equal to the disparity in the prediction score $S$, written $P(s\mid
x_{1})-P(s\mid x_{0})$. Term II measures the contribution coming from
thresholding the prediction score to obtain an optimal 0/1 prediction (98%
disparity).

Male ($x_{1}$)Female ($x_{0}$)Before ($Y$)After ($\widehat{Y}$)ThresholdingNot
hired ($y_{0}$)Hired ($y_{1}$)

Figure 1: Visualization of hiring disparities from Ex. 1.

The above example illustrates a canonical point in a simple setting: a small
disparity in the outcome $Y$, and consequently the prediction score $S$, may
result in a large disparity in the optimal 0/1 predictor $\widehat{Y}$, a case
we call bias amplification. Contrary to this, a large disparity in $S$ may
also result in a small disparity in $\widehat{Y}$, a case we call bias
amelioration.

In the remainder of the manuscript, our goal is to provide a decomposition of
the disparity in a thresholded predictor $\widehat{Y}$ into the disparity in
true outcome $Y$ and the disparity originating from optimization procedure,
but along each causal pathway between the protected attribute $X$ and the
predictor $\widehat{Y}$. In particular, our contributions are the following:

  1. (1)

We introduce the notion of margin complement (Def. 1), and provide a path-
specific decomposition of the disparity in the 0/1 predictor $\widehat{Y}$
into its contributions from the optimal score predictor $S$ and the margin
complement $M$ (Thm. 1),

  2. (2)

We prove that under suitable assumptions, the causal decomposition of the
optimal prediction score $S$ is equivalent with the causal decomposition of
the true outcome $Y$ (Thm. 2). This allows us to obtain a new decomposition of
the disparity in $\widehat{Y}$ into contributions from $Y$ and the margin
complement $M$ (Cor. 3),

  3. (3)

Motivated by the above decompositions, we introduce a new concept of weak and
strong business necessity (Def. 3), highlighting a new need for regulatory
instructions in the context of automated systems. We provide an algorithm for
assessing fairness under considerations of weak and strong business necessity
(Alg. 1),

  4. (4)

We provide identification, estimation, and sample influence results for all of
the quantities relevant to the above framework (Props. 4, 5). We evaluate our
approach on three real-world examples (Ex. 2-3) and provide new empirical
insights into bias amplification.

$Z$$X$$W$$Y$$\widehat{Y}$

Figure 2: Standard Fairness Model.

Our work is related to the previous literature on causal fairness and the
causal decompositions appearing in this literature (Zhang and Bareinboim,
2018b; Plečko and Bareinboim, 2024). However, our approach offers an entirely
new causal decomposition into contributions from the true outcome $Y$ and the
margin complement $M$. Our work is also related to recent results showing that
focusing purely on prediction, and ignoring decision-making aspects, may lead
to inequitable outcomes and cause harm to marginalized groups (Plečko and
Bareinboim, 2024; Nilforoshan et al., 2022; Plecko and Bareinboim, 2024a),
highlighting a need to expand focus from narrow statistical definitions of
fair predictions to a more comprehensive understanding of equity in
algorithmic decisions. Finally, our work is also related to the literature
auditing and assessing fairness of decisions made by humans (Pierson et al.,
2021; Kleinberg et al., 2018), and understanding how AI systems may help
humans overcome their biases (Imai et al., 2023).

###  1.1 Preliminaries

We use the language of structural causal models (SCMs) (Pearl, 2000). An SCM
is a tuple $\mathcal{M}:=\langle V,U,\mathcal{F},P(u)\rangle$ , where $V$, $U$
are sets of endogenous (observable) and exogenous (latent) variables,
respectively, $\mathcal{F}$ is a set of functions $f_{V_{i}}$, one for each
$V_{i}\in V$, where $V_{i}\leftarrow f_{V_{i}}(\mathrm{pa}(V_{i}),U_{V_{i}})$
for some $\mathrm{pa}(V_{i})\subseteq V$ and $U_{V_{i}}\subseteq U$. $P(u)$ is
a strictly positive probability measure over $U$. Each SCM $\mathcal{M}$ is
associated to a causal diagram $\mathcal{G}$ (Bareinboim et al., 2022) over
the node set $V$ where $V_{i}\rightarrow V_{j}$ if $V_{i}$ is an argument of
$f_{V_{j}}$, and $V_{i}\dashleftarrow\\!\\!\\!\\!\\!\\!\dashrightarrow V_{j}$
if the corresponding $U_{V_{i}},U_{V_{j}}$ are not independent. An
instantiation of the exogenous variables $U=u$ is called a unit. By $Y_{x}(u)$
we denote the potential response of $Y$ when setting $X=x$ for the unit $u$,
which is the solution for $Y(u)$ to the set of equations obtained by
evaluating the unit $u$ in the submodel $\mathcal{M}_{x}$, in which all
equations in $\mathcal{F}$ associated with $X$ are replaced by $X=x$.
Throughout the paper, we assume a specific cluster causal diagram
$\mathcal{G}_{\text{SFM}}$ known as the standard fairness model (SFM) (Plečko
and Bareinboim, 2024) over endogenous variables $\\{X,Z,W,Y,\widehat{Y}\\}$
shown in Fig. 2 (see also (Anand et al., 2023)). The SFM consists of the
following: protected attribute, labeled $X$ (e.g., gender, race, religion),
assumed to be binary; the set of confounding variables $Z$, which are not
causally influenced by the attribute $X$ (e.g., demographic information, zip
code); the set of mediator variables $W$ that are possibly causally influenced
by the attribute (e.g., educational level or other job-related information);
the outcome variable $Y$ (e.g., GPA, salary); the predictor of the outcome
$\widehat{Y}$ (e.g., predicted GPA, predicted salary). The SFM also encodes
the lack-of-confounding assumptions typically used in the causal inference
literature. The availability of the SFM and the implied assumptions are a
possible limitation of the paper, while we note that partial identification
techniques for bounding effects can be used for relaxing them (Zhang et al.,
2022).

##  2 Margin Complements

We begin by introducing a quantity that plays a key role in the results of
this paper.

######  Definition 1 (Margin Complement).

Let $U=u$ be a unit, and let $S$ denote a prediction score for a binary
outcome $Y$. Let the subscript $C$ denote a counterfactual clause, so that
$Z_{C}$ denotes a potential response. The margin complement $M$ of the score
$S$ for the unit $U=u$ and threshold $t$ is defined as:

$\displaystyle M(u)=\mathbbm{1}(S(u)\geq t)-S(u).$ (4)

A potential response of $M$, labeled $M_{C}$, is given by
$M_{C}(u)=\mathbbm{1}(S_{C}(u)\geq t)-S_{C}(u).$

In words, the margin complement for a unit $U=u$ represents the difference in
the score after thresholding vs. the score that would happen naturally.

######  Example 1 (Disparities in Hiring continued).

Consider the hiring example with the SCM in Eqs. 1-2 with
$p_{0}=0.49,p_{1}=0.51$. The unit $(U_{X},U_{Y})=(0,0)$ corresponds to a male
applicant ($X(u)=1$) who was hired ($Y(u)=1$). We have $S(u)=p_{1}=0.51$, and
$M(u)=\mathbbm{1}(S(u)\geq 0.5)-S(u)=1-0.51=0.49$. For this $u$, the margin
complement indicates that the predicted outcome
$\widehat{Y}(u)=\mathbbm{1}(S(u)\geq 0.5)$ is 49% greater than the predicted
probability $S(u)$. The same computation can be done for a female
$X(u^{\prime})=0$, in which case $M(u^{\prime})=\mathbbm{1}(p_{x_{0}}\geq
0.5)-p_{x_{0}}=-0.49$, meaning that the predicted outcome
$\widehat{Y}(u^{\prime})$ is 49% smaller than the predicted probability
$S(u^{\prime})$.

Given a prediction score $S$ and a threshold $t$, the margin complement $M$
tells us in which direction the thresholded version $\mathbbm{1}(S(u)\geq t)$
moves compared to the score $S(u)$. A positive margin complement indicates
that a thresholded predictor is larger than the probability prediction, and a
negative margin complement the opposite. A similar reasoning holds for the
potential responses of the margin complement $M_{C}$: we are interested in
what the margin complement would have been for an individual $U=u$ under
possibly different, counterfactual conditions described by $C$. As we
demonstrate shortly, margin complements (and their potential responses) play a
major role in explaining how inequities are generated between groups at the
time of decision-making. In this section, our key aim is to analyze the
optimal 0/1 predictor $\widehat{Y}$ and provide a decomposition of its total
variation measure (TV, for short), defined as
$\text{TV}_{x_{0},x_{1}}(\widehat{y})=P(\widehat{y}\mid
x_{1})-P(\widehat{y}\mid x_{0})$. When working with the causal diagram in Fig.
2, we can notice that the TV measure comprises of three types of variations
coming from $X$: the direct effect $X\to\widehat{Y}$, the mediated effect
$X\to W\to\widehat{Y}$, and the confounded effect
$X\dashleftarrow\\!\\!\\!\\!\\!\\!\dashrightarrow Z\to\widehat{Y}$. Our goal
is to construct a decomposition of the TV measure that allows us to
distinguish how much of each of the causal effects is due to a difference in
the prediction score $S$, and how much due to margin complements $M$. To
investigate this, we first introduce the known definitions of direct,
indirect, and spurious effects from the causal fairness literature:

######  Definition 2 ($x$-specific Causal Measures (Zhang and Bareinboim,
2018b; Plečko and Bareinboim, 2024)).

The $x$-specific {direct, indirect, spurious} effects of $X$ on $Y$ are
defined as $x\text{-DE}_{x_{0},x_{1}}(y\mid x)=P(y_{x_{1},W_{x_{0}}}\mid
x)-P(y_{x_{0}}\mid x)$, $x\text{-IE}_{x_{1},x_{0}}(y\mid
x)=P(y_{x_{1},W_{x_{0}}}\mid x)-P(y_{x_{1}}\mid x)$, and
$x\text{-SE}_{x_{1},x_{0}}(y)=P(y_{x_{1}}\mid x_{0})-P(y_{x_{1}}\mid x_{1})$,
respectively.

Armed with these definitions, we can prove the following result:

######  Theorem 1 (Causal Decomposition of Optimal 0/1 Predictor).

Let $\widehat{Y}$ be the optimal predictor with respect to the 0/1-loss based
on covariates $X,Z,W$. Let $S$ denote the optimal predictor with respect to
the $L_{2}$ loss. The total variation (TV, for short) measure of the predictor
$\widehat{Y}$, written as $P(\hat{y}\mid x_{1})-P(\hat{y}\mid x_{0})$, can be
decomposed into direct, indirect, and spurious effects of $X$ on the score $S$
and the margin complement $M$ as follows:

$\displaystyle\text{TV}_{x_{0},x_{1}}(\widehat{y})$
$\displaystyle={x\text{-DE}_{x_{0},x_{1}}(s\mid
x_{0})}+{x\text{-DE}_{x_{0},x_{1}}(m\mid x_{0})}$ (5)
$\displaystyle\quad-\big{(}{x\text{-IE}_{x_{1},x_{0}}(s\mid
x_{0})}+{x\text{-IE}_{x_{1},x_{0}}(m\mid x_{0})}\big{)}$ (6)
$\displaystyle\quad-\big{(}{x\text{-SE}_{x_{1},x_{0}}(s)}+{x\text{-SE}_{x_{1},x_{0}}(m)}\big{)}.$
(7)

The above theorem is the first key result of this paper. The disparity between
groups with respect to the optimal 0/1-loss predictor, measured by
TV${}_{x_{0},x_{1}}(\hat{y})$ can be decomposed into direct, indirect, and
spurious contributions coming from (i) the optimal $L_{2}$-loss predictor $S$
(e.g., term $x$-DE${}_{x_{0},x_{1}}(s\mid x_{0})$), and (ii) the margin
complement $M$ (e.g., term $x$-DE${}_{x_{0},x_{1}}(m\mid x_{0})$). This
provides a unique capability since for each causal pathway (direct, indirect,
spurious) the contribution coming from the probability prediction $S$ can be
disentangled from the contribution coming from the optimization procedure
itself (i.e., the rounding of the predictor). The former, as we will see
shortly, is simply a representation of the bias already existing in the true
outcome $Y$, whereas the latter represents a newly introduced type of bias
that is the result of using an automated system. The contribution of the
margin complement may act to both ameliorate or amplify an existing disparity,
a point we investigate later on.

######  Example 1 (Disparities in Hiring extended).

Consider the hiring example from Ex. 1 extended with a mediator $W$ indicating
whether the applicant has a PhD degree ($W=1$) or not ($W=0$). Suppose the
augmented SCM is given by:

$\displaystyle X$ $\displaystyle\leftarrow\mathbb{1}(U_{X}<0.5)$ (8)
$\displaystyle W$ $\displaystyle\leftarrow\mathbb{1}(U_{W}<0.5+\lambda X)$ (9)
$\displaystyle Y$ $\displaystyle\leftarrow\mathbb{1}(U_{Y}<0.1+\alpha X+\beta
W),$ (10)

and $P(U_{X},U_{W},U_{Y})$ is such that $U_{X},U_{W},U_{Y}$ are independent
$\text{Unif}[0,1]$ random variables. The coefficients satisfy the constraints
$\alpha,\lambda>0,0.4<\beta<\frac{0.9-\alpha}{1+\lambda}$. The optimal
probability predictor is $S(x,w)=0.1+\alpha x+\beta w$. The TV measure of the
optimal 0/1 predictor be computed as
$\text{TV}_{x_{0},x_{1}}(\hat{y})=\frac{\alpha}{\beta}+\lambda$. Using Thm. 1
we can decompose it as:

$\displaystyle\text{TV}_{x_{0},x_{1}}(\widehat{y})=\underbrace{\alpha}_{x\text{-DE}_{x_{0},x_{1}}(s\mid
x_{0})}+\underbrace{\frac{\alpha}{\beta}-\alpha}_{{x\text{-DE}_{x_{0},x_{1}}(m\mid
x_{0})}}-\underbrace{(-\lambda\beta)}_{{x\text{-IE}_{x_{1},x_{0}}(s\mid
x_{0})}}-\underbrace{(\lambda-\lambda\beta)}_{{x\text{-IE}_{x_{1},x_{0}}(m\mid
x_{0})}}.$ (11)

Along the direct effect, there is a difference of $\alpha$ between the groups
in the probability predictor $S(x,w)$ that is amplified by
$\frac{\alpha}{\beta}-\alpha$ by the margin complement, and since $\beta<1$ we
always have bias amplification. Along the indirect effect, there is an initial
difference of $\lambda\beta$, amplified by an additional
$\lambda-\lambda\beta$, again amplifying bias. The decomposition provides a
breakdown of the disparity of the optimal 0/1 predictor into its constitutive
parts.

We next consider a crucial point mentioned above. The disparity observed in
the optimal $L_{2}$ score $S$ is simply a representation of the disparity
existing in the real world, under suitable causal assumptions:

######  Theorem 2 (Relationship of $L_{2}$-score $S$ and Outcome $Y$
Decompositions).

Let $\mathcal{M}$ be an SCM compatible with the Standard Fairness Model, and
let $S$ be the optimal $L_{2}$ prediction score. Then, the causal
decompositions of the score $S$ and the true outcome $Y$ are symmetric,
meaning that ${x\text{-DE}_{x_{0},x_{1}}(s\mid
x_{0})}={x\text{-DE}_{x_{0},x_{1}}(y\mid x_{0})}$,
${x\text{-IE}_{x_{1},x_{0}}(s\mid x_{0})}={x\text{-IE}_{x_{1},x_{0}}(y\mid
x_{0})}$, and ${x\text{-SE}_{x_{1},x_{0}}(s)}={x\text{-SE}_{x_{1},x_{0}}(y)}$.

Based on this result, we can see that the causal influences from $X$ on the
true outcome $Y$ are equivalent to the causal influences from $X$ on the
optimal prediction score $S$. Combining this insight with Thm. 1 leads to the
following result:

######  Corollary 3 (Causal Decomposition of Optimal 0/1 Predictor).

Under the Standard Fairness Model, the TV measure of the optimal 0/1-loss
predictor $\widehat{Y}$ can be decomposed into contributions from the true
outcome $Y$ and the margin complement $M$:

$\displaystyle\text{TV}_{x_{0},x_{1}}(\widehat{y})$
$\displaystyle={x\text{-DE}_{x_{0},x_{1}}(y\mid
x_{0})}+{x\text{-DE}_{x_{0},x_{1}}(m\mid
x_{0})}-\big{(}{x\text{-IE}_{x_{1},x_{0}}(y\mid
x_{0})}+{x\text{-IE}_{x_{1},x_{0}}(m\mid x_{0})}\big{)}$ (12)
$\displaystyle\quad-\big{(}{x\text{-SE}_{x_{1},x_{0}}(y)}+{x\text{-SE}_{x_{1},x_{0}}(m)}\big{)}.$
(13)

This corollary provides an important insight compared to Thm. 1: the disparity
in the optimal predictor $\widehat{Y}$ can be decomposed into the contribution
inherited from the true outcome $Y$ and the contribution that arises from the
optimization procedure (rounding of the predictor).

##  3 Weak and Strong Business Necessity

In law, questions of fairness and discrimination can be interpreted based on
the disparate treatment (DT) and disparate impact (DI) doctrines of Title VII
of the Civil Rights Act of 1964 (Act, 1964). The DT doctrine, when interpreted
causally, can be seen as disallowing a direct type of effect from $X$ onto the
outcome. The DI doctrine, however, is more broad, and may prohibit any from of
discrimination (be it direct, indirect, or spurious) that results in a large
disparity between groups. The core of this doctrine is the notion of business
necessity (BN), which allows certain variables correlated with the protected
attribute to be used for prediction due to their relevance to the business
itself (Els, 1993) (or more broadly the utility of the decision-maker). Based
on the decomposition from Cor. 3, new BN considerations emerge:

######  Definition 3 (Weak and Strong Business Necessity).

Let $\mathcal{M}$ be an SCM compatible with the Standard Fairness Model. Let
CE denote a causal pathway (DE, IE, or SE). If a causal pathway does not fall
under business necessity, then we require:

$\displaystyle x\text{-CE}_{x,x^{\prime}}(s\mid
x^{\prime\prime})=x\text{-CE}_{x,x^{\prime}}(m\mid x^{\prime\prime})=0.$ (14)

A pathway is said to satisfy weak business necessity if:

$\displaystyle x\text{-CE}_{x,x^{\prime}}(s\mid
x^{\prime\prime})=x\text{-CE}_{x,x^{\prime}}(y\mid
x^{\prime\prime}),\;x\text{-CE}_{x,x^{\prime}}(m\mid x^{\prime\prime})=0.$
(15)

A pathway is said to satisfy strong business necessity if:

$\displaystyle x\text{-CE}_{x,x^{\prime}}(s\mid
x^{\prime\prime})=x\text{-CE}_{x,x^{\prime}}(y\mid
x^{\prime\prime}),\;x\text{-CE}_{x,x^{\prime}}(m\mid x^{\prime\prime})\text{
unconstrained}.$ (16)

The above definition distinguishes between three important cases, and sheds
light on a new aspect of the concept of business necessity. According to the
definition, there are three versions of BN considerations: (1) A causal
pathway is not in the BN set, and is considered discriminatory. In this case,
both the contribution of the prediction score $S$ and the margin complement
$M$ need to be equal to $0$ (i.e., no discrimination is allowed along the
pathway); (2) A causal pathway satisfies weak BN, and is not considered
discriminatory. In this case, the effect of $X$ on the prediction score $S$
needs to equal the effect of $X$ onto the true outcome $Y$ along the same
pathway (Plecko and Bareinboim, 2024b). However, the contribution of the
margin complement $M$ along the pathway needs to equal $0$; (3) A causal
pathway satisfies strong BN, and is not considered discriminatory. Similarly
as for weak necessity, the effect of $X$ on $S$ needs to equal the effect of
$X$ on $Y$, but in this case, the contribution of the margin complement $M$ is
unconstrained. The distinction between cases (2) and (3) opens the door for
new regulatory requirements and specifications. In particular, whenever a
causal effect is considered non-discriminatory, the attribute $X$ needs to
affect $S$ to the extent to which it does in the real world. However, the
system designer also needs to decide whether a difference existing in the
predicted probabilities $S$ is allowed to be amplified (or ameliorated) by
means of rounding. The latter point distinguishes between weak and strong BN,
and should be a consideration of any system designer issuing binary decisions.
In Alg. 1, we propose a formal approach for evaluating considerations of weak
and strong BN for any input of a predictor $\widehat{Y}$ and a prediction
score $S$.

Input: data $\mathcal{D}$, BN-Set
$\subseteq\\{\text{DE},\text{IE},\text{SE}\\}$, BN-Strength, predictor
$\widehat{Y}$, score $S$, SFM

Output: SUCCESS or FAIL of ensuring that disparate impact and treatment hold
under BN

1

2foreach  _$\text{CE}\in\\{\text{DE},\text{IE},\text{SE}\\}$_ do

3 Compute the effects $x\text{-CE}(y)$, $x\text{-CE}(m)$, $x\text{-CE}(s)$,
$x\text{-CE}(\widehat{y})$

4 if  _CE $\in$ Strong-BN_ then

5 Assert that $x\text{-CE}(s)=x\text{-CE}(\widehat{y})$, otherwise FAIL

6 else if  _CE $\in$ Weak-BN_ then

7 Assert that $x\text{-CE}(s)=x\text{-CE}(\widehat{y})\wedge
x\text{-CE}(m)=0$, otherwise FAIL

8 else

9 Assert that $x\text{-CE}(\widehat{s})=x\text{-CE}(\widehat{m})=0$, otherwise
FAIL

10

if not FAIL return SUCCESS

Algorithm 1 Auditing Weak & Strong Business Necessity

##  4 Identification, Estimation, and Sample Influence

In Thm. 1 and Cor. 3 the observed disparity in the TV measure of the optimal
0/1 predictor is decomposed into its constitutive components. The quantities
appearing in the decomposition are counterfactuals, and thus the question of
identification of these quantities needs to be addressed. In other words, we
need to understand whether these quantities can be uniquely computed based on
the available data and the causal assumptions. The following is a positive
answer:

######  Proposition 4 (Identification and Estimation of Causal Measures).

Let $\mathcal{M}$ be an SCM compatible with the Standard Fairness Model, and
let $P(V)$ be its observational distribution. The $x$-specific direct,
indirect, and spurious effects of $X$ on the outcome $Y$, predictor
$\widehat{Y}$, prediction score $S$, and the margin complement $M$ are
identifiable (uniquely computable) from $P(V)$ and the SFM. Denote by
$f(x,z,w)$ estimator of $\mathbbm{E}[T\mid x,z,w]$, and by $\hat{P}(x\mid
v^{\prime})$ the estimator of the probability $P(x\mid v^{\prime})$ for
different choices of $v^{\prime}$. For $T\in\\{Y,\widehat{Y},S,M\\}$, the
effects can be estimated as:

$\displaystyle x\text{-DE}^{\mathrm{est}}_{x_{0},x_{1}}(t\mid x_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}[f(x_{1},w_{i},z_{i})-f(x_{0},w_{i},z_{i})]\frac{\hat{P}(x_{0}\mid
w_{i},z_{i})}{\hat{P}(x_{0})}$ (17) $\displaystyle
x\text{-IE}^{\mathrm{est}}_{x_{1},x_{0}}(t\mid x_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}f(x_{1},w_{i},z_{i})\Big{[}\frac{\hat{P}(x_{0}\mid
w_{i},z_{i})}{\hat{P}(x_{0})}-\frac{\hat{P}(x_{1}\mid
w_{i},z_{i})}{\hat{P}(x_{1}\mid z_{i})}\frac{\hat{P}(x_{0}\mid
z_{i})}{\hat{P}(x_{0})}\Big{]}$ (18) $\displaystyle
x\text{-SE}^{\mathrm{est}}_{x_{1},x_{0}}(t)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}f(x_{1},w_{i},z_{i})\Big{[}\frac{\hat{P}(x_{1}\mid
w_{i},z_{i})}{\hat{P}(x_{1}\mid z_{i})}\frac{\hat{P}(x_{0}\mid
z_{i})}{\hat{P}(x_{0})}-\frac{\hat{P}(x_{1}\mid
w_{i},z_{i})}{\hat{P}(x_{1})}\Big{]}.$ (19)

The proof of the proposition, together with the identification expressions for
the different quantities can be found in Appendix B. We next define the sample
influences for the different estimators:

######  Definition 4 (Sample Influence).

The sample influence of the $i$-th sample on the estimator
$x\text{-CE}^{\mathrm{est}}$ of the causal effect CE is given by corresponding
term in the summations in Eqs. 17-19. For instance, the $i$-th sample
influence on $x\text{-DE}^{\mathrm{est}}_{x_{0},x_{1}}(t\mid x_{0})$ is given
by (and analogously for IE, SE terms):

$\displaystyle\text{SI-
DE}(i)=[f(x_{1},w_{i},z_{i})-f(x_{0},w_{i},z_{i})]\frac{\hat{P}(x_{0}\mid
w_{i},z_{i})}{\hat{P}(x_{0})}.$ (20)

The sample influences tell us how each of the samples contributes to the
overall estimator of the quantity. These sample-level contributions may be
interesting to investigate from the point of view of the system designer,
including identifying any subpopulations that are discriminated against. For
direct sample influences, the following proposition can be proved:

######  Proposition 5 (Direct Effect Sample Influence).

The $\text{SI-DE}(i)$ in Eq. 20 is an estimator of

$\displaystyle\mathbbm{E}[T_{x_{1},W_{x_{0}}}-T_{x_{0}}\mid
x_{0},z_{i},w_{i}]\frac{P(w_{i},z_{i}\mid x_{0})}{P(w_{i},z_{i})},$ (21)

where $\mathbbm{E}[T_{x_{1},W_{x_{0}}}-T_{x_{0}}\mid x_{0},z_{i},w_{i}]$ is
the $(x_{0},z_{i},w_{i})$-specific direct effect of $X$ on $T$.

Prop. 5 demonstrates an important point – namely that the sample influences
along the direct path are not just quantities of statistical interest, but
also causally meaningful quantities. In particular, the influence of the
$i$-th sample is proportional to the direct effect of the $x_{0}\to x_{1}$
transition for the group of units $u$ compatible with the event
$x_{0},z_{i},w_{i}$. The influence is further proportional to
$P(w_{i},z_{i}\mid x_{0})/P(w_{i},z_{i})$ that measures how much more likely
the covariates $z_{i},w_{i}$ of the $i$-th sample are in the $X=x_{0}$ group
(for which the discrimination is quantified) vs. the overall population.
Therefore, practitioners also have a causal reason for investigating these
sample influences.

Figure 3: Causal decomposition from Cor. 3 and sample influence on the MIMIC-
IV dataset.

##  5 Experiments

We analyze the MIMIC-IV (Ex. 2), COMPAS (Ex. 3), and Census (Ex. 4, Appendix
C) datasets.

######  Example 2 (Acute Care Triage on MIMIC-IV Dataset (Johnson et al.,
2023)).

Clinicians in the Beth Israel Deaconess Medical Center in Boston,
Massachusetts treat critically ill patients admitted to the intensive care
unit (ICU). For all patients, various physiological and treatment information
is collected 24 hours after admission, and the available data consists of
(grouped into the Standard Fairness model): protected attribute $X$, in this
case race ($x_{0}$ African-American, $x_{1}$ White), set of confounders $Z=$
{sex, age, chronic health status}, set of mediators $W=${lactate, SOFA score,
admission diagnosis, PaO2/FiO2 ratio, aspartate aminotransferase}.

Clinicians are interested in patients who require closer monitoring. They want
to determine the top half of the patients who are the most likely to (i) die
during their hospital stay; (ii) have an ICU stay longer than 10 days. This
combined outcome is labeled $Y$. These high-risk patients will remain in the
most acute care unit. To predict the outcome, clinicians use the electronic
health records (EHR) data of the hospital and construct score predictions $S$
and a binary predictor
$\widehat{Y}=\mathbbm{1}(S(x,z,w)>\mathrm{Quant}(0.5;S))$ that selects the top
half of the patients.

To investigate the fairness implications of the new AI-based system, they use
the decomposition described in Cor. 3 to investigate different contributions
to the resulting disparity. The decomposition is shown in Fig. 3(a), and
uncovers a number of important effects. Firstly, along the direct effect,
$x$-DE${}_{x_{0},x_{1}}(y)$ and $x$-DE${}_{x_{0},x_{1}}(m)$ are larger than
$0$, meaning that minority group individuals have a lower chance of receiving
acute care (purely based on race). Along the indirect and spurious effects,
the situation is different: $x$-IE${}_{x_{1},x_{0}}(y)$ and
$x$-SE${}_{x_{1},x_{0}}(y)$ and their respective margin complement
contributions are different than $0$ and negative – implying that minority
group individuals have a larger probability of being given acute care as a
result of confounding and mediating variables. Finally, the direct effect
sample influences (Fig. 3(b)) highlight that the margin complements are large
for a small minority of individuals, requiring further subgroup investigation
by the hospital team.

Figure 4: Application of Alg. 1 on the COMPAS dataset.

######  Example 3 (Recidivism Prevention on the COMPAS Dataset (Larson et al.,
2016)).

Courts in Broward County, Florida use machine learning algorithms, developed
by a private company called Northpointe, to predict whether individuals
released on parole are at high risk of re-offending within 2 years ($Y$). The
algorithm is based on the demographic information $Z$ ($Z_{1}$ for gender,
$Z_{2}$ for age), race $X$ ($x_{0}$ denoting White, $x_{1}$ Non-White),
juvenile offense counts $J$, prior offense count $P$, and degree of charge
$D$. The courts wish to know which individuals are highly likely to
recidivate, such that their probability of recidivism is above 50%. The
company constructs a prediction score $\hat{S}^{NP}$ and the court
subsequently uses this for deciding whether to detain individuals at high risk
of re-offending.

After a court hearing in which it was decided that the indirect and spurious
effects fall under business necessity requirements, a team from ProPublica
wishes to investigate the implications of using the automated predictions
$\hat{S}^{NP}$. They obtain the relevant data and apply Alg. 1, with the
results shown in Fig. 4. The team first compares the decompositions of the
true outcome $Y$ and the predictor $\hat{S}^{NP}$ (Fig. 4(a)). For the
spurious effect, they find that $x$-SE${}_{x_{1},x_{0}}(y)$ is not
statistically different from $x$-SE${}_{x_{1},x_{0}}(\hat{s}^{NP})$, in line
with BN requirements. For the indirect effects, they find that the indirect
effect is lower for the predictor $\hat{S}^{NP}$ compared to the true outcome
$Y$, indicating no concerning violations. However, for the direct effect,
while the $x$-DE${}_{x_{0},x_{1}}(y)$ is not statistically different from $0$,
the predictor $\hat{S}^{NP}$ has a significant direct effect of $X$, i.e.,
$x$-DE${}_{x_{0},x_{1}}(\hat{s}^{NP})\neq 0$. This indicates a violation of
the fairness requirements determined by the court.

After comparing the decompositions of $\hat{S}^{NP}$ and $Y$, the team moves
onto understanding the contributions of the margin complements (Fig. 4b). For
each effect, there is a pronounced impact of the margin complements. For the
direct effect (not under BN), the non-zero margin complement contribution
$x$-DE${}_{x_{1},x_{0}}(m)\neq 0$ represents a violation of fairness
requirements. For the indirect and spurious effects, the ProPublica team
realizes the court did not specify anything about margin complement
contributions – based on this, for the next court hearing they are preparing
an argument showing that the effects $x$-IE${}_{x_{1},x_{0}}(m)$ and
$x$-SE${}_{x_{1},x_{0}}(m)$ are significantly different from $0$, thereby
exacerbating the differences between groups. Finally, based on sample
influences (Fig. 4c), they realize that the direct effect is driven by a small
minority of individuals, and they decide to investigate this further.

##  6 Conclusion

In this paper, we developed tools for understanding the fairness impacts of
transforming a continuous prediction score $S$ into binary predictions
$\widehat{Y}$ or binary decisions $D$. In Thm. 1 and Cor. 3 we showed that the
TV measure of the optimal 0/1 predictor decomposes into direct, indirect, and
spurious contributions that are inherited from the true outcome $Y$ in the
real world, and also contributions from the margin complement $M$ (Def. 1)
arising from the automated optimization procedure. This observation motivated
new notions of weak and strong business necessity (BN) – in the former case,
differences inherited from the true outcome $Y$ are allowed to be propagated
into predictions or decisions, while any differences resulting from the
optimization procedure are disallowed. In contrast, strong BN allows both of
these differences and does not prohibit possible disparity amplification. In
Alg. 1, we developed a formal procedure for assessing weak and strong BN.
Finally, real-world examples demonstrated that the tools developed in this
paper are of genuine importance in practice, since converting continuous
predictions into binary decisions may often result in bias amplification in
practice – highlighting the need for this type of analysis, and the importance
of regulatory oversight.

## References

  * Els (1993) Elston v. Talladega County Bd. of Educ.  997 F.2d 1394 (11th Cir. 1993), 1993.  United States Court of Appeals for the Eleventh Circuit.
  * Act (1964) C. R. Act.  Civil rights act of 1964.  _Title VII, Equal Employment Opportunities_ , 1964.
  * Anand et al. (2023) T. V. Anand, A. H. Ribeiro, J. Tian, and E. Bareinboim.  Causal effect identification in cluster dags.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 37(10), pages 12172–12179, 2023.
  * Angwin et al. (2016) J. Angwin, J. Larson, S. Mattu, and L. Kirchner.  Machine bias: There’s software used across the country to predict future criminals. and it’s biased against blacks.  _ProPublica_ , 5 2016.  URL https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-sentencing.
  * Bareinboim et al. (2022) E. Bareinboim, J. D. Correa, D. Ibeling, and T. Icard.  On pearl’s hierarchy and the foundations of causal inference.  In _Probabilistic and Causal Inference: The Works of Judea Pearl_ , page 507–556. Association for Computing Machinery, New York, NY, USA, 1st edition, 2022.
  * Blau and Kahn (1992) F. D. Blau and L. M. Kahn.  The gender earnings gap: learning from international comparisons.  _The American Economic Review_ , 82(2):533–538, 1992.
  * Blau and Kahn (2017) F. D. Blau and L. M. Kahn.  The gender wage gap: Extent, trends, and explanations.  _Journal of economic literature_ , 55(3):789–865, 2017.
  * Brennan et al. (2009) T. Brennan, W. Dieterich, and B. Ehret.  Evaluating the predictive validity of the compas risk and needs assessment system.  _Criminal Justice and Behavior_ , 36(1):21–40, 2009.
  * Buolamwini and Gebru (2018) J. Buolamwini and T. Gebru.  Gender shades: Intersectional accuracy disparities in commercial gender classification.  In S. A. Friedler and C. Wilson, editors, _Proceedings of the 1st Conference on Fairness, Accountability and Transparency_ , volume 81 of _Proceedings of Machine Learning Research_ , pages 77–91, NY, USA, 2018.
  * Chiappa (2019) S. Chiappa.  Path-specific counterfactual fairness.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 7801–7808, 2019.
  * Chouldechova (2017) A. Chouldechova.  Fair prediction with disparate impact: A study of bias in recidivism prediction instruments.  Technical Report arXiv:1703.00056, arXiv.org, 2017.
  * Darlington (1971) R. B. Darlington.  Another look at “cultural fairness” 1.  _Journal of educational measurement_ , 8(2):71–82, 1971.
  * Datta et al. (2015) A. Datta, M. C. Tschantz, and A. Datta.  Automated experiments on ad privacy settings: A tale of opacity, choice, and discrimination.  _Proceedings on Privacy Enhancing Technologies_ , 2015(1):92–112, Apr. 2015.  doi: 10.1515/popets-2015-0007.
  * Imai et al. (2023) K. Imai, Z. Jiang, D. J. Greiner, R. Halen, and S. Shin.  Experimental evaluation of algorithm-assisted human decision-making: Application to pretrial public safety assessment.  _Journal of the Royal Statistical Society Series A: Statistics in Society_ , 186(2):167–189, 2023.
  * Johnson et al. (2023) A. E. Johnson, L. Bulgarelli, L. Shen, A. Gayles, A. Shammout, S. Horng, T. J. Pollard, B. Moody, B. Gow, L.-w. H. Lehman, et al.  Mimic-iv, a freely accessible electronic health record dataset.  _Scientific data_ , 10(1):1, 2023.
  * Khandani et al. (2010) A. E. Khandani, A. J. Kim, and A. W. Lo.  Consumer credit-risk models via machine-learning algorithms.  _Journal of Banking & Finance_, 34(11):2767–2787, 2010.
  * Kilbertus et al. (2017) N. Kilbertus, M. Rojas-Carulla, G. Parascandolo, M. Hardt, D. Janzing, and B. Schölkopf.  Avoiding discrimination through causal reasoning.  _arXiv preprint arXiv:1706.02744_ , 2017.
  * Kleinberg et al. (2018) J. Kleinberg, H. Lakkaraju, J. Leskovec, J. Ludwig, and S. Mullainathan.  Human decisions and machine predictions.  _The quarterly journal of economics_ , 133(1):237–293, 2018.
  * Kusner et al. (2017) M. J. Kusner, J. Loftus, C. Russell, and R. Silva.  Counterfactual fairness.  _Advances in neural information processing systems_ , 30, 2017.
  * Larson et al. (2016) J. Larson, S. Mattu, L. Kirchner, and J. Angwin.  How we analyzed the compas recidivism algorithm.  _ProPublica (5 2016)_ , 9, 2016.
  * Mahoney and Mohen (2007) J. F. Mahoney and J. M. Mohen.  Method and system for loan origination and underwriting, Oct. 23 2007.  US Patent 7,287,008.
  * Nabi and Shpitser (2018) R. Nabi and I. Shpitser.  Fair inference on outcomes.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Nilforoshan et al. (2022) H. Nilforoshan, J. D. Gaebler, R. Shroff, and S. Goel.  Causal conceptions of fairness and their consequences.  In _International Conference on Machine Learning_ , pages 16848–16887. PMLR, 2022.
  * Pager (2003) D. Pager.  The mark of a criminal record.  _American journal of sociology_ , 108(5):937–975, 2003.
  * Pearl (2000) J. Pearl.  _Causality: Models, Reasoning, and Inference_.  Cambridge University Press, New York, 2000.  2nd edition, 2009.
  * Pierson et al. (2021) E. Pierson, D. M. Cutler, J. Leskovec, S. Mullainathan, and Z. Obermeyer.  An algorithmic approach to reducing unexplained pain disparities in underserved populations.  _Nature Medicine_ , 27(1):136–140, 2021.
  * Plečko and Bareinboim (2024) D. Plečko and E. Bareinboim.  Causal fairness analysis: A causal toolkit for fair machine learning.  _Foundations and Trends® in Machine Learning_ , 17(3):304–589, 2024.
  * Plecko and Bareinboim (2024a) D. Plecko and E. Bareinboim.  Causal fairness for outcome control.  _Advances in Neural Information Processing Systems_ , 36, 2024a.
  * Plecko and Bareinboim (2024b) D. Plecko and E. Bareinboim.  Reconciling predictive and statistical parity: A causal approach.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 38 (13), pages 14625–14632, 2024b.
  * Plečko and Meinshausen (2020) D. Plečko and N. Meinshausen.  Fair data adaptation with quantile preservation.  _Journal of Machine Learning Research_ , 21:242, 2020.
  * Sanburn (2015) J. Sanburn.  Facebook thinks some native american names are inauthentic.  _Time_ , Feb. 14 2015.  URL http://time.com/3710203/facebook-native-american-names/.
  * Shpitser and Pearl (2007) I. Shpitser and J. Pearl.  What counterfactuals can be tested.  In _Proceedings of the Twenty-third Conference on Uncertainty in Artificial Intelligence_ , page 352–359, 2007.
  * Sweeney (2013) L. Sweeney.  Discrimination in online ad delivery.  Technical Report 2208240, SSRN, Jan. 28 2013.  URL http://dx.doi.org/10.2139/ssrn.2208240.
  * Sweeney and Haney (1992) L. T. Sweeney and C. Haney.  The influence of race on sentencing: A meta-analytic review of experimental studies.  _Behavioral Sciences & the Law_, 10(2):179–195, 1992.
  * Wu et al. (2019) Y. Wu, L. Zhang, X. Wu, and H. Tong.  Pc-fairness: A unified framework for measuring causality-based fairness.  _Advances in neural information processing systems_ , 32, 2019.
  * Zhang and Bareinboim (2018a) J. Zhang and E. Bareinboim.  Equality of opportunity in classification: A causal approach.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31_ , pages 3671–3681, Montreal, Canada, 2018a. Curran Associates, Inc.
  * Zhang and Bareinboim (2018b) J. Zhang and E. Bareinboim.  Fairness in decision-making—the causal explanation formula.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018b.
  * Zhang et al. (2022) J. Zhang, J. Tian, and E. Bareinboim.  Partial counterfactual identification from observational and experimental data.  In _Proceedings of the 39th International Conference on Machine Learning_ , 2022.

## Technical Appendices for Mind the Gap: A Causal Perspective on Bias
Amplification in Prediction & Decision-Making

The source code for reproducing all the experiments can be found in the
anonymized code repository. The code is also included with the supplementary
materials, in the folder source-code. All experiments were performed on a
MacBook Pro, with the M3 Pro chip and 36 GB RAM on macOS 14.1 (Sonoma). All
experiments can be run with less than 1 hour of compute on the above-described
machine or equivalent.

##  Appendix A Proofs of Key Theorems

###### Thm. 1 Proof:.

Notice that by definition we have

$\displaystyle M_{C}(u)=\mathbb{1}(S_{C}(u)\geq\frac{1}{2})-S_{C}(u)\;\forall
u\in\mathcal{U}.$ (22)

Now, let $E$ be any observed event, and let $C$ be a counterfactual clause
representing a possibly counterfactual intervention. Note that we can write:

$\displaystyle\mathbbm{E}[\widehat{Y}_{C}\mid E]\overset{\text{(def)}}{=}$
$\displaystyle\sum_{u}\widehat{Y}_{C}(u)P(u\mid
E)=\sum_{u}\mathbb{1}(S_{C}(u)\geq\frac{1}{2})P(u\mid E)$ (23)
$\displaystyle=$
$\displaystyle\sum_{u}\Big{(}\mathbb{1}(S_{C}(u)\geq\frac{1}{2})-S_{C}(u)+S_{C}(u)\Big{)}P(u\mid
E)$ (24) $\displaystyle=$
$\displaystyle\sum_{u}\Big{(}\mathbb{1}(S_{C}(u)\geq\frac{1}{2})-S_{C}(u)\Big{)}P(u\mid
E)+\sum_{u}S_{C}(u)P(u\mid E)$ (25) $\displaystyle=$
$\displaystyle\sum_{u}M_{C}(u)P(u\mid E)+\sum_{u}S_{C}(u)P(u\mid E)$ (26)
$\displaystyle=$ $\displaystyle\mathbbm{E}[M_{C}\mid E]+\mathbbm{E}[S_{C}\mid
E].$ (27)

Now, we can expand the TV measure of $\widehat{Y}$ as follows:

$\displaystyle\text{TV}_{x_{0},x_{1}}(\hat{y})$
$\displaystyle=\mathbbm{E}[\widehat{Y}\mid x_{1}]-\mathbbm{E}[\widehat{Y}\mid
x_{0}]$ (28) $\displaystyle=\mathbbm{E}[\widehat{Y}_{x_{1}}\mid
x_{1}]-\mathbbm{E}[\widehat{Y}_{x_{1}}\mid x_{0}]$ (29)
$\displaystyle\quad+\mathbbm{E}[\widehat{Y}_{x_{1}}\mid
x_{0}]-\mathbbm{E}[\widehat{Y}_{x_{1},W_{x_{0}}}\mid x_{0}]$ (30)
$\displaystyle\quad+\mathbbm{E}[\widehat{Y}_{x_{1},W_{x_{0}}}\mid
x_{0}]-\mathbbm{E}[\widehat{Y}_{x_{0}}\mid x_{0}]$ (31)
$\displaystyle=\mathbbm{E}[S_{x_{1}}\mid x_{1}]-\mathbbm{E}[S_{x_{1}}\mid
x_{0}]+\mathbbm{E}[M_{x_{1}}\mid x_{1}]-\mathbbm{E}[M_{x_{1}}\mid x_{0}]$ (32)
$\displaystyle\quad+\mathbbm{E}[S_{x_{1}}\mid
x_{0}]-\mathbbm{E}[S_{x_{1},W_{x_{0}}}\mid x_{0}]+\mathbbm{E}[M_{x_{1}}\mid
x_{0}]-\mathbbm{E}[M_{x_{1},W_{x_{0}}}\mid x_{0}]$ (33)
$\displaystyle\quad+\mathbbm{E}[S_{x_{1},W_{x_{0}}}\mid
x_{0}]-\mathbbm{E}[S_{x_{0}}\mid x_{0}]+\mathbbm{E}[M_{x_{1},W_{x_{0}}}\mid
x_{0}]-\mathbbm{E}[M_{x_{0}}\mid x_{0}]$ (34)
$\displaystyle=-{x\text{-SE}_{x_{1},x_{0}}(s)}-{x\text{-SE}_{x_{1},x_{0}}(m)}$
(35) $\displaystyle\quad-x\text{-IE}_{x_{1},x_{0}}(s\mid
x_{0})-{x\text{-IE}_{x_{1},x_{0}}(m\mid x_{0})}$ (36)
$\displaystyle\quad+{x\text{-DE}_{x_{0},x_{1}}(s\mid
x_{0})}+{x\text{-DE}_{x_{0},x_{1}}(m\mid x_{0})},$ (37)

completing the proof. ∎

###### Thm. 2 Proof:.

Let $\mathcal{M}$ be an SCM compatible with the SFM in Fig. 2 as per theorem
assumption. Let the optimal $L_{2}$ prediction score $S$ be given by
$S(x,z,w)=\mathbbm{E}[Y\mid x,z,w]$. More precisely, we can let $f_{S}$ be the
structural mechanism of $S$, taking $X,Z,W$ as inputs. The structural
mechanism of $S$ is then given by:

$\displaystyle f_{S}(x,z,w)=\mathbbm{E}\left[Y\mid x,z,w\right].$ (38)

Therefore, the score $S$ is a deterministic function of $X,Z,W$, and we can
add it to the standard fairness model as shown in Fig. 5(a). Now, note that
for a potential outcome $\mathbbm{E}[S_{x,W_{x^{\prime}}}\mid
x^{\prime\prime}]$ we can write:

$\displaystyle\mathbbm{E}[S_{x,W_{x^{\prime}}}\mid x^{\prime\prime}]$
$\displaystyle=\sum_{z}\mathbbm{E}[S_{x,W_{x^{\prime}}}\mid
x^{\prime\prime},z]P(z\mid x^{\prime\prime})$ (39)
$\displaystyle=\sum_{z,w}\mathbbm{E}[S_{x,w,z}\mathbb{1}(W_{x^{\prime}}=w)\mid
x^{\prime\prime},z]P(z\mid x^{\prime\prime})$ (40)
$\displaystyle=\sum_{z,w}\mathbbm{E}[S_{x,w,z}\mid
x^{\prime\prime},z]\mathbbm{E}[\mathbb{1}(W_{x^{\prime}}=w)\mid
x^{\prime\prime},z]P(z\mid x^{\prime\prime})$ (41)
$\displaystyle=\sum_{z,w}\mathbbm{E}[S_{x,w,z}]P(W_{x^{\prime}}=w\mid
x^{\prime\prime},z)P(z\mid x^{\prime\prime})$ (42)
$\displaystyle=\sum_{z,w}\mathbbm{E}[Y_{x,w,z}]P(W_{x^{\prime}}=w\mid
x^{\prime\prime},z)P(z\mid x^{\prime\prime})$ (43)
$\displaystyle=\sum_{z,w}\mathbbm{E}[Y_{x,w,z}\mathbb{1}(W_{x^{\prime}}=w)\mid
x^{\prime\prime},z]P(z\mid x^{\prime\prime})$ (44)
$\displaystyle=\sum_{z}\mathbbm{E}[Y_{x,W_{x^{\prime}}}\mid
x^{\prime\prime},z]P(z\mid
x^{\prime\prime})=\mathbbm{E}[Y_{x,W_{x^{\prime}}}\mid x^{\prime\prime}],$
(45)

where Eq. 41 is using the independence
$S_{x,z,w}{\perp\\!\\!\\!\perp}W_{x^{\prime}}\mid X,Z$, and Eq. 44 is using
the independence $Y_{x,z,w}{\perp\\!\\!\\!\perp}W_{x^{\prime}}\mid X,Z$, both
of which are implied by the standard fairness model extended with the node $S$
above. Based on the equality $\mathbbm{E}[S_{x,W_{x^{\prime}}}\mid
x^{\prime\prime}]=\mathbbm{E}[Y_{x,W_{x^{\prime}}}\mid x^{\prime\prime}]$ for
arbitrary values $x,x^{\prime},x^{\prime\prime}$, the claim follows from the
appropriate choices: for instance, taking
$\\{x=x_{1},x^{\prime}=x_{0},x^{\prime\prime}=x_{0}\\}$ and
$\\{x=x_{0},x^{\prime}=x_{0},x^{\prime\prime}=x_{0}\\}$, it follows that

$\displaystyle x\text{-DE}_{x_{0},x_{1}}(s)$
$\displaystyle=\mathbbm{E}[S_{x_{1},W_{x_{0}}}\mid
x_{0}]-\mathbbm{E}[S_{x_{0},W_{x_{0}}}\mid x_{0}]$ (46)
$\displaystyle=\mathbbm{E}[Y_{x_{1},W_{x_{0}}}\mid
x_{0}]-\mathbbm{E}[Y_{x_{0},W_{x_{0}}}\mid
x_{0}]=x\text{-DE}_{x_{0},x_{1}}(y).$ (47)

Similarly, analogous choices of $x,x^{\prime},x^{\prime\prime}$ can be used to
show the equality for indirect and spurious effects. ∎

##  Appendix B Proofs related to Identification & Estimation

$Z$$X$$W$$Y$$S$

(a) SFM extended with predictor $S$.

$Z$$X$$W$$T_{x_{1},w}$$W_{x_{0}}$${x_{1}}$${x_{0}}$${w}$ (b) Counterfactual
graph of the SFM.

Figure 5: Graphs used in proofs of Thm. 2 and Prop. 4.

###### Prop. 4 Proof:.

Consider the $x$-specific {direct, indirect, spurious} effects of $X$ on a
random variable $T$ given by:

$\displaystyle x\text{-DE}_{x_{0},x_{1}}(t\mid x_{0})$
$\displaystyle=P(t_{x_{1},W_{x_{0}}}\mid x_{0})-P(t_{x_{0}}\mid x_{0})$ (48)
$\displaystyle x\text{-IE}_{x_{1},x_{0}}(t\mid x_{0})$
$\displaystyle=P(t_{x_{1},W_{x_{0}}}\mid x_{0})-P(t_{x_{1}}\mid x_{0})$ (49)
$\displaystyle x\text{-SE}_{x_{1},x_{0}}(t)$ $\displaystyle=P(t_{x_{1}}\mid
x_{0})-P(t_{x_{1}}\mid x_{1}).$ (50)

We first demonstrate why each of the potential outcomes appearing is
identifiable under the Standard Fairness Model. Consider the potential outcome
$T_{x_{1},W_{x_{0}}}\mid X=x_{0}$. For proving its identifiability, we make
use of the counterfactual graph [Shpitser and Pearl, 2007] in Fig. 5(b). We
can expand $P(T_{x_{1},W_{x_{0}}}=t\mid x_{0})$ as follows:

$\displaystyle=\sum_{w}P(T_{x_{1},w}=t,W_{x_{0}}=w\mid
x_{0})\quad\text{(Counterfactual Un-nesting)}$ (51)
$\displaystyle=\sum_{z,w}P(T_{x_{1},w}=t,W_{x_{0}}=w\mid z,x_{0})P(z\mid
x_{0})\quad\text{(Law of Total Probability)}$ (52)
$\displaystyle=\sum_{z,w}P(T_{x_{1},w}=t\mid z,x_{0})P(W_{x_{0}}=w\mid
z,x_{0})P(z\mid x_{0})\quad(T_{x_{1},w}{\perp\\!\\!\\!\perp}W_{x_{0}}\mid
Z,X)$ (53) $\displaystyle=\sum_{z,w}P(T_{x_{1},w}=t\mid z,x_{0})P(W=w\mid
z,x_{0})P(z\mid x_{0})\quad\text{(Consistency)}$ (54)
$\displaystyle=\sum_{z,w}P(T_{x_{1},w}=t\mid z,w,x_{0})P(W=w\mid
z,x_{0})P(z\mid x_{0})\quad(T_{x_{1},w}{\perp\\!\\!\\!\perp}W\mid Z,X)$ (55)
$\displaystyle=\sum_{z,w}P(T_{x_{1},w}=t\mid z,w,x_{1})P(W=w\mid
z,x_{0})P(z\mid x_{0})\quad(T_{x_{1},w}{\perp\\!\\!\\!\perp}X\mid Z)$ (56)
$\displaystyle=\sum_{z,w}P(T=t\mid z,w,x_{1})P(W=w\mid z,x_{0})P(z\mid
x_{0})\quad\text{(Consistency)}$ (57) $\displaystyle=\sum_{z,w}P(t\mid
z,w,x_{1})P(w\mid z,x_{0})P(z\mid x_{0}).$ (58)

Using a similar but simplified argument, we also obtain that:

$\displaystyle P(T_{x_{1}}=t\mid x_{0})=\sum_{z}P(t\mid x_{1},z)P(z\mid
x_{0}),$ (59)

and since $P(T_{x}=t\mid x)=P(t\mid x)$, all the terms in Eqs. 48-50 can be
computed based on observational data.

We next move onto the estimation part of the proof. We again focus on the
quantity $P(T_{x_{1},W_{x_{0}}}=t\mid x_{0})$. Suppose that we have:

$\displaystyle f(x_{1},z,w)$ $\displaystyle\xrightarrow{P}\mathbbm{E}[T\mid
x_{1},z,w]\;\forall z,w$ (60) $\displaystyle\hat{P}(x_{0}\mid w,z)$
$\displaystyle\xrightarrow{P}P(x_{0}\mid w,z)\;\forall z,w$ (61)
$\displaystyle\hat{P}(w,z)$ $\displaystyle\xrightarrow{P}P(w,z)\;\forall w,z$
(62) $\displaystyle\hat{P}(x_{0})$ $\displaystyle\xrightarrow{P}P(x_{0})$ (63)

where $f,\hat{P}(x_{0}\mid w,z)$ are estimators, and
$\hat{P}(w,z),\hat{P}(x_{0})$ are the empirical distributions:

$\displaystyle\hat{P}(w,z)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}=z,W_{i}=w)$ (64)
$\displaystyle\hat{P}(x_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}(X_{i}=x_{0}).$ (65)

Now, note that we can re-write the identification expression from Eq. 58 as:

$\displaystyle\sum_{z,w}P(t\mid z,w,x_{0})P(z,w)\frac{P(x_{0}\mid
z,w)}{P(x_{0})}.$ (66)

By a coupling of quantities in Eqs. 60-63 to a joint probability space and an
application of the continuous mapping theorem we obtain that

$\displaystyle\sum_{z,w}f(x_{1},z,w)\hat{P}(z,w)\frac{\hat{P}(x_{0}\mid
z,w)}{\hat{P}(x_{0})}\xrightarrow{P}\sum_{z,w}P(t\mid
z,w,x_{1})P(z,w)\frac{P(x_{0}\mid z,w)}{P(x_{0})}.$ (67)

Finally, note that we have

$\displaystyle\sum_{z,w}f(x_{1},z,w)\hat{P}(z,w)\frac{\hat{P}(x_{0}\mid
z,w)}{\hat{P}(x_{0})}$
$\displaystyle=\sum_{z,w}f(x_{1},z,w)\Big{(}\sum_{i=1}^{n}\frac{\mathbbm{1}(Z_{i}=z,W_{i}=w)}{n}\Big{)}\frac{\hat{P}(x_{0}\mid
z,w)}{\hat{P}(x_{0})}$ (68)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}f(x_{1},z_{i},w_{i})\frac{\hat{P}(x_{0}\mid
z_{i},w_{i})}{\hat{P}(x_{0})}.$ (69)

The remaining identification expression are derived based on the same
technique, completing the proof of the proposition. ∎

###### Prop. 5 Proof:.

Consider the direct effect sample influence of the $i$-th sample given by

$\displaystyle\text{SI-
DE}(i)=[f(x_{1},w_{i},z_{i})-f(x_{0},w_{i},z_{i})]\frac{\hat{P}(x_{0}\mid
w_{i},z_{i})}{\hat{P}(x_{0})}.$ (70)

By assumption, we know that

$\displaystyle f(x,z,w)$ $\displaystyle\xrightarrow{P}\mathbbm{E}[T\mid
x_{1},z,w]\;\forall x,z,w$ (71) $\displaystyle\hat{P}(x_{0}\mid w,z)$
$\displaystyle\xrightarrow{P}P(x_{0}\mid w,z)\;\forall z,w$ (72)
$\displaystyle\hat{P}(w,z)$ $\displaystyle\xrightarrow{P}P(w,z)\;\forall w,z$
(73) $\displaystyle\hat{P}(x_{0})$ $\displaystyle\xrightarrow{P}P(x_{0})$ (74)

Therefore, by a coupling of the above quantities to a joint probability space
and an application of the continuous mapping theorem, we have that

$\displaystyle\text{SI-DE}(i)\xrightarrow{P}$
$\displaystyle\Big{[}\mathbbm{E}[T\mid x_{1},w_{i},z_{i}]-\mathbbm{E}[T\mid
x_{0},w_{i},z_{i}]\Big{]}\times\frac{P(x_{0}\mid w_{i},z_{i})}{P(x_{0})}.$
(75)

Note that

$\displaystyle\frac{P(x_{0}\mid
w_{i},z_{i})}{P(x_{0})}=\frac{P(x_{0},w_{i},z_{i})}{P(w_{i},z_{i})P(x_{0})}=\frac{P(w_{i},z_{i}\mid
x_{0})}{P(w_{i},z_{i})},$ (76)

completing the proof of the proposition. ∎

##  Appendix C Census 2018 Example

In this appendix, we analyze the Census 2018 dataset, demonstrating an
application in the context of labor and salary decisions:

Figure 6: Causal decomposition from Cor. 3 and sample influence on the Census
2018 dataset.

######  Example 4 (Salary Increase on the Census 2018 Dataset [Plečko and
Bareinboim, 2024]).

The United States Census of 2018 collected broad information about the US
Government employees, including demographic information $Z$ ($Z_{1}$ for age,
$Z_{2}$ for race, $Z_{3}$ for nationality), gender $X$ ($x_{0}$ female,
$x_{1}$ male), marital and family status $M$, education information $L$, and
work-related information $R$. The US Government wishes to use this data
prospectively to decide whether a new employee should receive a bonus starting
package upon signing the employment contract. To determine which employees
should receive such a bonus, a Government department decides to predict which
of the employees should earn a salary above the median (which is
$50,000/year). They construct a machine learning prediction score $S$ that
predicts above-median earnings ($Y$), and the department decides to allocate
the bonus to all employees who are predicted to be above-median earners with a
probability greater than 50% (i.e., $\widehat{Y}=\mathbbm{1}(S\geq 0.5)$).

A team of investigators within the Government is in charge of assessing what
the impacts of this AI system are on the gender pay gap. They collect the
required data and perform the decomposition from Cor. 3, shown in Fig. 6(a).
The decomposition indicates strong direct effects $x$-DE${}_{x_{0},x_{1}}(y)$
and $x$-DE${}_{x_{0},x_{1}}(m)$, which imply that men are more likely to
receive a starting bonus than women. The indirect effects
$x$-IE${}_{x_{1},x_{0}}(y)$, $x$-IE${}_{x_{1},x_{0}}(m)$ are both negative,
with the latter effect of the margin complement not being significant. These
effects, however, also mean that men are more likely to receive a bonus due to
mediating variables. Finally, for the spurious effects, the effects are not
significantly different from $0$. Based on these findings, the team decided to
return the predictions to the original department, with the requirement that
the direct effect of gender on the margin complement,
$x$-DE${}_{x_{0},x_{1}}(m)$, must be reduced to $0$, to avoid any possibility
of bias amplification.


11institutetext:  ⋆ King’s College London
† Siemens Healthineers

# A multi-organ point cloud registration algorithm
for abdominal CT registration

Samuel Joutard⋆,† Thomas Pheiffer† Chloe Audigier†
Patrick Wohlfahrt† Reuben Dorent⋆ Sebastien Piat†
Tom Vercauteren⋆ Marc Modat⋆ Tommaso Mansi†

###### Abstract

Registering CT images of the chest is a crucial step for several tasks such as
disease progression tracking or surgical planning. It is also a challenging
step because of the heterogeneous content of the human abdomen which implies
complex deformations. In this work, we focus on accurately registering a
subset of organs of interest. We register organ surface point clouds, as may
typically be extracted from an automatic segmentation pipeline, by expanding
the Bayesian Coherent Point Drift algorithm (BCPD). We introduce MO-BCPD, a
multi-organ version of the BCPD algorithm which explicitly models three
important aspects of this task: organ individual elastic properties, inter-
organ motion coherence and segmentation inaccuracy. This model also provides
an interpolation framework to estimate the deformation of the entire volume.
We demonstrate the efficiency of our method by registering different patients
from the LITS challenge dataset. The target registration error on anatomical
landmarks is almost twice as small for MO-BCPD compared to standard BCPD while
imposing the same constraints on individual organs deformation.

##  1 Introduction

Registering CT images of the chest is an important step for several pipelines
such as surgical planning for liver cancer resection or disease progression
tracking [2, 1, 10, 15]. This step is both crucial and challenging as the
deformations involved are large and may contain complex patterns such as
sliding motion between organs. While traditional registration methods tend to
fail on this task, learning approaches such as [12, 6, 5] obtained promising
results at the Learn2Reg 2020 challenge, task 3 [7]. Yet, traditional and
learning approaches both aims at registering the whole image content instead
of focusing on the relevant structures of interests. This introduces undesired
noise and complexity to the registration process. To tackle this issue, we
propose to exploit the recent availability of high quality automatic
segmentation pipelines such as [17, 3] and register the segmented structures.
Specifically, structures are registered using their surface point cloud
representation, allowing for exploiting meaningful geometric information of
the different organs and finely modeling their dynamic properties. We also
stress that surface point clouds are easy to derive from segmentation masks
and are a lightweight representation of the structures of interest.

The Coherent Point Drift [13] (CPD) algorithm is one of the most popular
method for deformable point cloud registration considered as state of the art
[11]. A recent work [9] extended this framework using a Bayesian formulation
and obtained more robust performances. CDP and BCPD both assume that points
move coherently as a group to preserve the structure coherence. This is mainly
because these frameworks are designed to register point clouds representing a
single object. Consequently, [13, 9] are not adapted for registering multi-
organ points clouds. In particular, the coherency assumption doesn’t stand for
organs registration as each organ-specific point cloud may move independently
to its neighbour, especially if we aim at registering inter-patient images.

In this work, we introduce a Multi-Organ Bayesian Coherent Point Drift
algorithm (MO-BCPD) that models independent coherent structures. The
contribution of this work is four-fold. Firstly, we extend the Bayesian
formulation of CPD to model more complex structures interactions such as organ
motion independence. Secondly, given that points clouds are obtained using
automated segmentations, the proposed framework models partial segmentation
errors allowing MO-BCPD to recover them. Thirdly, we model individual organ
elasticity as part of the formulation. Fourthly, extensive experiments on 104
patients (10,712 pairs of patients) from the LiTS public dataset [4]
demonstrate the effectiveness of our approach compared to BCPD. In particular,
our method achieves an average target registration error on anatomical
landmarks of 13mm compared to 22mm for the standard BCPD.

##  2 Method

$\mathbf{v}\leftarrow 0_{M,3}$, $\Sigma\leftarrow Id_{M}$, $s\leftarrow 1$,
$R\leftarrow Id_{3}$, $t\leftarrow 0_{3}$,
$<\alpha_{m}>\leftarrow\frac{1}{M}$,
$\sigma^{2}\leftarrow\frac{\gamma}{D\sum\limits_{m,n}u_{l^{y}_{m},l^{x}_{n}}}\sum\limits_{m,n}u_{l^{y}_{m},l^{x}_{n}}\left\lVert
x_{n}-y_{m}\right\rVert^{2}$,
$\theta\leftarrow(\mathbf{v},\mathbf{\alpha},\mathbf{c},\mathbf{e},\rho,\sigma^{2})$,
$P\leftarrow\frac{1}{M}\mathbf{1}_{M,N}$ $\nu^{\prime}\leftarrow 1_{N}$,
$q_{1}(.,.)\leftarrow D^{\kappa\mathbf{1}_{M}}\phi^{0,\Sigma}$,
$q_{2}(c,e)\leftarrow\prod\limits_{n=1}^{N}(1-\nu^{\prime}_{n})^{1-c_{n}}\left(\nu^{\prime}_{n}\prod\limits_{m=1}^{M}\left(\frac{p_{mn}}{\nu^{\prime}_{n}}\right)^{\delta_{n}(e_{m})}\right)^{c_{n}}$,
$q_{3}(.,.)\leftarrow\delta_{\rho}\delta_{\sigma^{2}}$

while  _$L(q_{1}q_{2}q_{3})$ increases more than $\epsilon$_ do

Update $P$ and related terms: $\forall m,n\;\phi_{m,n}\leftarrow
u_{l^{y}_{m},l^{x}_{n}}\phi^{y_{m}^{\prime},\sigma^{2}Id_{3}}(x_{n})\exp{-\frac{3s^{2}\Sigma_{m,m}}{2\sigma^{2}}}$,
$\forall m,n\;p_{m,n}\leftarrow\frac{(1-\omega)<\alpha_{m}>\phi_{m,n}}{\omega
p_{out}(x_{n})+(1-\omega)\sum\limits_{m^{\prime}}<\alpha_{m^{\prime}}>\phi_{m^{\prime},n}}$,
$\nu\leftarrow P.1_{N}$, $\nu^{\prime}\leftarrow P^{T}.1_{M}$,
$\hat{N}\leftarrow\nu^{T}.1_{M}$,
$\hat{\mathbf{x}}\leftarrow\Delta(\nu)^{-1}.P.\mathbf{x}$,

Update displacement field and related terms:
$\Sigma\leftarrow\left(G^{-1}+\frac{s^{2}}{\sigma^{2}}\Delta(\nu)\right)$,
$\forall
d\in\\{1,2,3\\}\;\mathbf{v}^{d}\leftarrow\frac{s^{2}}{\sigma^{2}}\Sigma\Delta(\nu)(\rho^{-1}(\hat{\mathbf{x}}^{d})-\mathbf{y}^{d})$,
$\mathbf{u}\leftarrow\mathbf{y}+\mathbf{v}$, $<\alpha_{m}>\leftarrow
exp\\{\psi(\kappa+\nu_{m})-\psi(\kappa M+\hat{N})\\}$

Update $\rho$ and related terms:
$\bar{x}\leftarrow\frac{1}{\hat{N}}\sum\limits_{m=1}^{M}\nu_{m}\hat{x}_{m}$,
$\bar{\sigma}^{2}\leftarrow\frac{1}{\hat{N}}\sum\limits_{m=1}^{M}\nu_{m}\sigma_{m}^{2}$,
$\bar{u}\leftarrow\frac{1}{\hat{N}}\sum\limits_{m=1}^{M}\nu_{m}u_{m}$,
$S_{xu}\leftarrow\frac{1}{\hat{N}}\sum\limits_{m=1}^{M}(\hat{x}_{m}-\bar{x})(u_{m}-\bar{u})^{T}$,
$S_{uu}\leftarrow\frac{1}{\hat{N}}\sum\limits_{m=1}^{M}(u_{m}-\bar{u})(u_{m}-\bar{u})^{T}+\bar{\sigma}^{2}Id_{3}$,
$\Phi S_{xu}^{\prime}\Psi^{T}\leftarrow svd(S_{xu})$, $R\leftarrow\Phi
d(1,\ldots,1,|\Phi\Psi|)\Psi^{T}$,
$s\leftarrow\frac{Tr(RS_{xu})}{Tr(S_{uu})}$, $t\leftarrow\bar{x}-sR\bar{u}$,
$\mathbf{y^{\prime}}\leftarrow\rho(\mathbf{y}+\mathbf{v})$
$\sigma^{2}\leftarrow\frac{1}{3\hat{N}}\sum\limits_{d=1}^{3}\left((\mathbf{x}^{d})^{T}\Delta(\nu^{\prime})\mathbf{x}^{d}-2\mathbf{x}^{d}P^{T}\mathbf{y^{\prime}}^{d}+(\mathbf{y^{\prime}}^{d})^{T}\Delta(\nu)\mathbf{y^{\prime}}^{d}\right)+s^{2}\bar{\sigma}^{2}$

Update q: $q_{1}(.,.)\leftarrow
D^{\kappa\mathbf{1}_{M}}\phi^{\mathbf{v},\Sigma}$,
$q_{2}(c,e)\leftarrow\prod\limits_{j=1}^{N}(1-\nu^{\prime}_{j})^{1-c_{j}}\left(\nu^{\prime}_{j}\prod\limits_{i=1}^{M}\left(\frac{p_{ij}}{\nu^{\prime}_{j}}\right)^{\delta_{i}(e_{j})}\right)^{c_{j}}$,
$q_{3}(.,.)\leftarrow\delta_{\rho}\delta_{\sigma^{2}}$

end while

Algorithm 1 Multi-Organ BCPD
$\;\left(\mathbf{y},\,\mathbf{x},\,\omega,\,\Lambda,\,B,\,S,\,U,\,\kappa,\,\gamma,\,\epsilon\right)$

In this section, we present our Multi-Organ Bayesian Coherent Point Drift
algorithm. Let $\mathbf{y}=[\mathbf{y}_{m}]_{m\in\\{1\ldots
M\\}}\in\mathbb{R}^{M,3}$ be the source point cloud and
$\mathbf{x}=[\mathbf{x}_{n}]_{n\in\\{1\ldots N\\}}\in\mathbb{R}^{N,3}$ be the
target point cloud where N and M are respectively the number of source and
target points. We aim at finding the transformation $\mathcal{T}$ that
realistically aligns these point clouds. In particular here, unlike in [13,
9], the considered point clouds both represent a set of organ surfaces. Hence,
each point is associated with an organ. Let
$\mathbf{l}^{y}=[l^{y}_{m}]_{m\in\\{1\ldots M\\}}\in\\{1\ldots L\\}^{M}$ be
the organ labels of the source point cloud and
$\mathbf{l}^{x}=[l^{x}_{n}]_{n\in\\{1\ldots N\\}}\in\\{1\ldots L\\}^{N}$ be
the organ labels of the target point cloud. L is the number of organs.

### Transformation model.

Similarly to the BCPD, the Multi-Organ Bayesian Coherent Point Drift (MO-
BCPD), decomposes the motion in two components: a similarity transform
$\rho:\mathbf{p}\longrightarrow s\mathbf{R}\mathbf{p}+\mathbf{t}$ and a dense
displacement field $\mathbf{v}$. Hence the deformed source point could is
$[\mathcal{T}(\mathbf{y}_{m})]_{m\in\\{1\ldots
M\\}}=[\rho(\mathbf{y}_{m}+\mathbf{v}_{m})]_{m\in\\{1\ldots M\\}}$. While this
parametrization is redundant, [9] has shown that this makes the algorithm more
robust to target rotations. Moreover, it is equivalent to performing a rigid
alignment followed by a non-rigid refinement which corresponds to the common
practice in medical image registration.

### Generative model.

As in [9], MO-BCPD assumes that all points from the target point cloud
$[\mathbf{x}_{n}]_{n\in\\{1\ldots N\\}}$ are sampled independently from a
generative model. A point $x_{n}$ from the target point cloud is either an
outlier or an inlier which is indicated by a hidden binary variable $c_{n}$.
We note the probability for a point to be an outlier $\omega$ (i.e.
$\mathcal{P}(c_{n}=0)=\omega$). If $x_{n}$ is an outlier, it is sampled from
an outlier distribution of density $p_{out}$ (typically, a uniform
distribution over a volume containing the target point cloud). If $x_{n}$ is
an inlier ($c_{n}=1$), $x_{n}$ is associated with a point $\mathcal{T}(y_{m})$
in the deformed source point cloud. Let $e_{n}$ be a multinomial variable
indicating the index of the point of the deformed source point cloud with
which $x_{n}$ is associated (i.e. $e_{n}=m$ in our example). Let $\alpha_{m}$
be the probability of selecting the point $\mathcal{T}(y_{m})$ to generate a
point of the target point cloud (i.e. $\forall
n\;\mathcal{P}(e_{n}=m|c_{n}=1)=\alpha_{m}$). $x_{n}$ is then sampled from a
Gaussian distribution with covariance-matrix $\sigma^{2}Id_{3}$ ($Id_{3}$ is
the identity matrix of $\mathbb{R}^{3}$) centered on
$\mathcal{T}(\mathbf{y}_{m})$. Finally, the organ label $l^{x}_{n}$ is sampled
according to the label transition distribution
$\mathcal{P}(l^{x}_{n}|l^{y}_{m})=u_{l^{x}_{n},l^{y}_{m}}$. The addition of
the label transition term is our contribution to the original generative model
[9]. This term encourages to map corresponding organs between the different
anatomies while allowing to recover from partial segmentation errors from the
automatic segmentation tool.

We can now write the following conditional probability density:

$p^{e}(x_{n},l^{x}_{n},c_{n},e_{n}|\mathbf{y},\mathbf{l}^{y},\mathbf{v},\alpha,\rho,\sigma^{2})\\\
=\left(\omega
p_{out}(x_{n})\right)^{1-c_{n}}\left((1-\omega)\prod\limits_{m=1}^{M}\left(\alpha_{m}u_{l^{y}_{m},l^{x}_{n}}\phi^{\mathbf{y}_{m}^{\prime},\sigma^{2}Id_{3}}(\mathbf{x}_{n})\right)^{\delta_{e_{n}=m}}\right)^{c_{n}}$
(1)

where $\phi^{\mu,\Sigma}$ is the density of a multivariate Gaussian
distribution $\mathcal{N}(\mu,\Sigma)$ and $\delta$ is the Kronecker symbol.

### Prior distributions

MO-BCPD also relies on prior distributions in order to regularize the
registration process and obtain realistic solutions. As in [9], MO-BCPD
defines two prior distributions: $p^{v}(\mathbf{v}|\mathbf{y},\mathbf{l}_{y})$
that regularizes the dense displacement field and $p^{\alpha}(\alpha)$ that
regularizes the parameters $\alpha$ of the source point cloud selection
multinomial distribution mentioned in the generative model. The prior on
$\alpha$ follows a Dirichlet distribution of parameter $\kappa\mathbf{1}_{M}$.
In practice, $\kappa$ is set to a very high value which forces
$\alpha_{m}\approx 1/M$ for all $m$. To decouple motion characteristics within
and between organs, we propose a novel formulation of the displacement field
prior $p^{v}$. Specifically, we introduce 3 parameters: a symmetric matrix
$S=[s_{l,l^{\prime}}]_{l,l^{\prime}\in\\{0\ldots L\\}}$ and two vectors
$\Lambda=[\Lambda_{l}]_{l\in\\{0\ldots L\\}}$ and $B=[B_{l}]_{l\in\\{0\ldots
L\\}}$. The matrix $S$ parametrizes the motion coherence inter-organs. The
vectors $\Lambda$ and $B$ respectively characterizes the variance of the
deformation magnitude and motion coherence bandwidth within each organ. We
define the displacement field prior for the MO-BCPD as:

$\displaystyle
p^{v}(\mathbf{v}|\mathbf{y},\mathbf{l}^{y})=\phi^{\mathbf{0},G}(\mathbf{v}^{1})\phi^{\mathbf{0},G}(\mathbf{v}^{2})\phi^{\mathbf{0},G}(\mathbf{v}^{3})$
(2) $\displaystyle
G=\left[\Lambda_{l^{y}_{i}}\Lambda_{l^{y}_{j}}S_{l^{y}_{i},l^{y}_{j}}\exp{-\frac{\left\lVert
y_{i}-y_{j}\right\rVert^{2}}{2B_{l^{y}_{i}}B_{l^{y}_{j}}}}\right]_{i,j\leq M}$
(3)

Note that $G$ must be definite-positive, leading to strictly positive values
for variance of the displacement magnitude $\Lambda_{l}$ and mild constraints
on $S$.

### Learning.

Combining equations (1) and (2), the joint probability distribution of the
variables $\mathbf{y},\mathbf{l}^{y},\mathbf{x},\mathbf{l}^{x},\theta$, where
$\theta=(\mathbf{v},\alpha,\mathbf{c},\mathbf{e},\rho,\sigma^{2})$ is defined
as:

$p(\mathbf{x},\mathbf{l}^{x},\mathbf{y},\mathbf{l}^{y},\theta)\propto
p^{v}(\mathbf{v}|\mathbf{y},\mathbf{l}_{y})p^{\alpha}(\alpha)\prod\limits_{n=1}^{N}p^{e}(x_{n},l^{x}_{n},c_{n},e_{n}|\mathbf{y},\mathbf{l}_{y},\mathbf{v},\alpha,\rho,\sigma^{2})$
(4)

As in [9], we use variational inference to approximate the posterior
distribution $p(\theta|\mathbf{x},\mathbf{y})$ with a factorized distribution
$q(\theta)=q_{1}(\mathbf{v},\alpha)q_{2}(\mathbf{c},\mathbf{e})q_{3}(\rho,\sigma^{2})$
so that
$q=\operatorname*{arg\,min}\limits_{q_{1},q_{2},q_{3}}KL(q|p(.|\mathbf{x},\mathbf{y}))$
where $KL$ is the Kullback–Leibler divergence. Similarly to [9], we derive the
MO-BCPD algorithm presented in Algorithm 1. The steps detailed in Algorithm 1
perform coordinate ascent on the evidence lower bound
$L(\theta)=\int_{\theta}q(\theta)\ln{\frac{p(\mathbf{x},\mathbf{y},\theta)}{q(\theta)}}d\theta$.
In algorithm 1, $\gamma$ is a hyper-parameter used to scale the initial
estimation of $\sigma^{2}$ and $\epsilon$ is used for stopping criteria. We
note $\Delta(\nu)$ the diagonal matrix with diagonal entries equal to $\nu$.

### Hyper-parameter setting.

The model has a large number of hyper-parameters which can impact the
performance of the algorithm. Regarding $\kappa$, the parameter of the prior
distribution $p^{\alpha}$, and $\gamma$, the scaling applied to the initial
estimation of $\sigma^{2}$, we followed the guidelines in [9]. $\omega$ is set
based on an estimate of the proportion of outliers on a representative testing
set. Regarding $B$ and $\Lambda$, respectively the vector of organ-specific
motion coherence bandwidth and expected deformation magnitude, they
characterise organs elastic properties. Concretely, a larger motion coherence
bandwidth $B_{l}$ increases the range of displacement correlation for organ
$l$ (points that are further away are encouraged to move in the same
direction). A larger expected deformation magnitude $\Lambda_{l}$ increases
the probability of larger displacements for organ $l$. These are physical
quantities expressed in $mm$ that could be set based on organs physical
properties. The inter-organ motion coherence matrix $S$ should be a symmetric
matrix containing values between 0 and 1. $S_{l,l}=1$ for all organs
$l\in\\{1\ldots L\\}$ and $S_{l,l^{\prime}}$ is closer to 0 if organs $l$ and
$l^{\prime}$ can move independently.

$u_{l,l^{\prime}}$ is the probability that a point with label $l^{\prime}$
generates a point with label $l$. As points labels are in practice obtained
from an automatic segmentation tool, we note $[g^{y}_{m}]_{m\in\\{1\ldots
M\\}}$ and $[g^{x}_{n}]_{n\in\\{1\ldots n\\}}$ respectively the unknown true
organ labels of the source and target point clouds (as opposed to the
estimated ones $[l^{y}_{m}]_{m\in\\{1\ldots M\\}}$ and
$[l^{y}_{m}]_{m\in\\{1\ldots M\\}}$). We assume that points from the deformed
source point cloud generate points with the same true labels (i.e.
$\mathcal{P}(g^{x}_{n}=g^{y}_{m}|e_{n}=m)=1$). Hence, the probability for
$y_{m}$, with estimated organ label $l^{y}_{m}$ to generate a point with label
$l^{x}_{n}$ is given by:
$u_{l^{x}_{n},l^{y}_{m}}=\sum_{k}p(g^{y}_{m}=k|l^{y}_{m})p(l^{x}_{n}|g^{x}_{n}=k)$
where $p(g^{y}_{m}=k|l^{y}_{m})$ is the probability that a point labelled
$l^{y}_{m}$ by the automatic segmentation tool has true label $k$ and
$p(l^{x}_{n}|g^{x}_{n}=k)$ is the probability that the automatic segmentation
tool predicts the organ label $l^{x}_{n}$ for a point with true label $k$.
These probabilities need to be estimated on a representative testing set. We
note that if the segmentation is error-free, the formula above gives
$U=Id_{L}$. Indeed, in that case the points organ labels correspond exactly to
the organ true labels so a point belonging to an certain organ can only
generate a point from the same organ in the target anatomy. Properly setting
the organ label transition probability matrix $U$ is crucial to recover from
potential partial segmentation errors. Fig. 1 illustrates with a toy example a
situation where the algorithm converges to an undesired state if the
segmentation error is not modeled properly.

### Interpolation.

Once the deformation on the organ point clouds is known, one might want to
interpolate the deformation back to image space in order to resample the whole
volume. As in [8] we propose to use Gaussian process regression to interpolate
the deformation obtained by the MO-BCPD algorithm. This interpolation process
can also be used to register sub-sampled point clouds to decrease computation
time as in [8].

Given a set of points
$\tilde{\mathbf{y}}=[\tilde{y}_{i}]_{i\in\\{1\ldots\tilde{M}\\}}$ with labels
$\mathbf{l}^{\tilde{y}}={[l^{\tilde{y}}_{i}]}_{i\in\\{1\ldots\tilde{M}\\}}$.
We compute the displacement for the set of points $\tilde{y}$ as:

$\displaystyle\mathbf{v}^{\tilde{\mathbf{y}}}=G^{int}(\tilde{\mathbf{y}},\mathbf{l}^{\tilde{y}},\mathbf{y},\mathbf{l}^{y},B,\Lambda,S).G^{-1}.v$
(5) $\displaystyle
G^{int}(\tilde{\mathbf{y}},\mathbf{l}^{\tilde{y}},\mathbf{y},\mathbf{l}^{y},B,\Lambda,S)_{i,j}=\Lambda_{l^{\tilde{y}}_{i}}\Lambda_{l^{y}_{i}}S_{l^{\tilde{y}}_{i},l^{y}_{i}}\exp{-\frac{\left\lVert\tilde{y}_{i}-y_{j}\right\rVert^{2}}{2\beta_{l^{\tilde{y}}_{i}}\beta_{l^{y}_{j}}}}$
(6)

Figure 1: Toy example registering a pair of organs (a blue and an orange
organ) with $\sim$10% segmentation error (corrupted input labels). Both organs
(orange and blue) of the target point cloud are shown in (a) in transparent
while the source point cloud is shown in opaque. The blue (orange) dots on the
left (right) of the figure corresponds to simulated segmentation errors. (b)
shows the registered point cloud without modeling the inter-organ segmentation
error , (c) shows the registered point cloud with segmentation error
modelization

### Acceleration.

The speed ups mentioned in [9] section 4.6 are fully transferable to the MO-
BCPD pipeline. In our experiments though, the main improvement, by far, came
from performing a low rank decomposition of $G$ at the initialization of the
algorithm. Indeed, this yielded consistent reliable $\times 10$ speed-ups with
negligible error when using $\geq 20$ eigen values. The Nystrom methods to
approximate $P$ sometime implied large error due to the stochasticity of the
method while yielding up to $\times 2$ speed-ups which is why we did not use
of it. This allows MO-BCPD to be run in a few seconds with $M,N\approx 5000$.

##  3 Experiments

Table 1: Target registration error on landmarks. Results in mm (std).  | Sim | BCPD | GMC-MO-BCPD | OMC-MO-BCPD
---|---|---|---|---
bladder | 257 (26) | 29 (15) | 30 (15) |  26 (15)
left kidney bottom | 128 (18) | 23 (11) | 22 (10) |  8 (4)
left kidney center | 107 (13) | 18 (10) | 15 (8) |  6 (3)
left kidney top | 101 (15) | 23 (12) | 21 (10) |  9 (4)
liver bottom | 114 (15) | 29 (14) | 28 (14) |  24 (13)
liver center | 65 (10) | 12 (7) | 12 (7) |  11 (7)
liver top | 123 (15) |  24 (12) | 25 (12) | 26 (14)
right kidney bottom | 98 (15) | 26 (13) | 24 (11) |  10 (6)
right kidney center | 65 (11) | 21 (12) | 17 (9) |  5 (3)
right kidney top | 64 (15) | 25 (13) | 21 (11) |  9 (4)
round ligament of liver | 95 (20) | 27 (14) | 27 (13) |  25 (13)

We evaluate the MO-BCPD algorithm by performing inter-patient registration
from the LITS challenge training dataset [4] which contains 131 chest-CT
patient images. 27 patients were removed due to different field of view,
issues with the segmentation or landmark detection. In total, 10,712
registrations were performed on all the pairs of remaining patients. The
segmentation was automatically performed using an in-house tool derived from
[17] which also provides a set of anatomical landmarks for each image which
were used for evaluation. We considered five organs of interest: the liver,
the spleen, the left and right kidneys and the bladder. We compared 4
different algorithms: registration of the point clouds with a similarity
transform (Sim), BCPD, GMC-MO-BCPD which is MO-BCPD with global motion
coherence ($S=1_{L,L}$) and OMC-MO-BCPD which is MO-BCPD with intra-organ
motion coherence only ($S=Id_{L}$). As the segmentation tool performed very
well on the considered organs, we set $U=Id_{L}$ and $\omega=0$ for both MO-
BCPD versions ($\omega=0$ for BCPD as well). We used for all organs the same
values for $\Lambda_{l}$ and $B_{l}$ respectively 10mm and 30mm as a trade off
between shape matching and preservation of individual organs appearance
($\beta=30$ and $\lambda=0.1$ for BCPD which is the equivalent configuration).
We also set, $\gamma=1$ and $\epsilon=0.1$. We compared those algorithms by
computing the registration error on the anatomical landmarks belonging to
those organs. We chose this generic, relatively simple setting (same rigidity
values for all organs, no outlier modeling, only two extreme configurations
for $S$) in order to perform large scale inter-patient registration
experiments but we would like to stress that further fine tuning of these
parameters for a specific application or even for a specific patient would
further improve the modeling and hence the registration outcome. Results are
presented in Table 1. We observe that while GMC-MO-BCPD induces some marginal
improvements with respect to BCPD, OMC-MO-BCPD allows a much more precise
registration. As illustrated in Fig. 2, the main improvement from BCPD to GMC-
MO-BCPD is that different organs no longer overlap. Indeed, as highlighted by
the green and blue ellipses, the spleen and the right kidney from the deformed
source patient overlap the liver of the target patient when using BCPD. OMC-
MO-BCPD properly aligns the different organs while preserving their shape (see
orange and purple ellipses for instance).

Figure 2: Qualitative comparison of registration output for the liver,
left/right kidneys and spleen. From left to right, BCPD, GMC-MO-BCPD, OMC-MO-
BCPD. The target organs are shown in transparency while the deformed point
cloud are shown in opaque.

##  4 Conclusion

We introduced MO-BCPD, an extension of the BCPD algorithm specifically adapted
to abdominal organ registration. We identified three limitations of the
original work [9] on this task and proposed solutions to model: the
segmentation error between neighboring organs of interest, the heterogeneous
elastic properties of the abdominal organs and the complex interaction between
various organs in terms of motion coherence. We demonstrated significant
improvements over BCPD on a large validation set (N=10,712).

Moreover, we would like to highlight that segmentation error could also be
taken into account by tuning the outlier probability distribution $p_{out}$
and the probability of being an outlier $\omega$. When the point is estimated
as a potential outlier by the algorithm, its contribution to the estimation of
the transformation $\mathcal{T}$ is lowered. Hence, the segmentation error
modelled by $p_{out}$ and $\omega$ corresponds to over/under segmentation,
i.e. when there is a confusion between an organ and another class we don’t
make use of in MO-BCPD (e.g. background). Hence, MO-BCPD introduces a finer
way of handling segmentation error by distinguishing two types of errors: mis-
labeling between classes of interest which is modelled by $U$ and over/under-
segmentation of classes of interest modeled by $\omega$ and $p_{out}$.

In this manuscript, we focused on highlighting the improvements yielded by the
MO-BCPD formulation specifically designed for multi-organ point cloud
registration. That being said, some clinical applications would require the
deformation on the whole original image volume. Hence, we present in
supplementary material preliminary results on a realistic clinical use case.
In figures 3 and 4 we see that MO-BCPD coupled with the proposed interpolation
framework obtain better alignment on the structures of interest than
traditional intensity-based baselines. It is also interesting to note that MO-
BCPD also better aligns structures that are particularly challenging to align
such as the hepatic vein while not using these structures in the MO-BCPD.

Future work will investigate how MO-BCPD could be used as a fast, accurate
initialization for image-based registration algorithms. From a modeling stand-
point, we would also like to further work on segmentation error modeling (in
particular over/under segmentation) with a more complex organ specific outlier
distribution.

## References

  * [1] Carrillo, A., Duerk, J., Lewin, J., Wilson, D.: Semiautomatic 3-d image registration as applied to interventional mri liver cancer treatment. IEEE Transactions on Medical Imaging 19(3), 175–185 (2000). https://doi.org/10.1109/42.845176
  * [2] Cash, D.M., Miga, M.I., Glasgow, S.C., Dawant, B.M., Clements, L.W., Cao, Z., Galloway, R.L., Chapman, W.C.: Concepts and preliminary data toward the realization of image-guided liver surgery. Journal of Gastrointestinal Surgery 11(7), 844–59 (07 2007), https://www.proquest.com/scholarly-journals/concepts-preliminary-data-toward-realization/docview/1112236808/se-2?accountid=11862, copyright - The Society for Surgery of the Alimentary Tract 2007; Dernière mise à jour - 2014-03-30
  * [3] Chen, X., Sun, S., Bai, N., Han, K., Liu, Q., Yao, S., Tang, H., Zhang, C., Lu, Z., Huang, Q., Zhao, G., Xu, Y., Chen, T., Xie, X., Liu, Y.: A deep learning-based auto-segmentation system for organs-at-risk on whole-body computed tomography images for radiation therapy. Radiotherapy and Oncology pp. 175–184 (Jul 2021)
  * [4] Christ, P.F., Vorontsov, E., Chlebus, G., Chen, H., Dou, Q., Fu, C., Han, X., Heng, P., Hesser, J., Kadoury, S., Konopczynski, T.K., Le, M., Li, C., Li, X., Lipková, J., Lowengrub, J.S., Meine, H., Moltz, J.H., Pal, C., Piraud, M., Qi, X., Qi, J., Rempfler, M., Roth, K., Schenk, A., Sekuboyina, A., Zhou, P., Hülsemeyer, C., Beetz, M., Ettlinger, F., Grün, F., Kaissis, G., Lohöfer, F., Braren, R., Holch, J., Hofmann, F., Sommer, W.H., Heinemann, V., Jacobs, C., Mamani, G.E.H., van Ginneken, B., Chartrand, G., Tang, A., Drozdzal, M., Ben-Cohen, A., Klang, E., Amitai, M.M., Konen, E., Greenspan, H., Moreau, J., Hostettler, A., Soler, L., Vivanti, R., Szeskin, A., Lev-Cohain, N., Sosna, J., Joskowicz, L., Menze, B.H.: The liver tumor segmentation benchmark (LiTS). CoRR abs/1901.04056 (2019), http://arxiv.org/abs/1901.04056
  * [5] Estienne, T., Vakalopoulou, M., Battistella, E., Carré, A., Henry, T., Lerousseau, M., Robert, C., Paragios, N., Deutsch, E.: Deep learning based registration using spatial gradients and noisy segmentation labels. CoRR abs/2010.10897 (2020), https://arxiv.org/abs/2010.10897
  * [6] Heinrich, M.P.: Closing the gap between deep and conventional image registration using probabilistic dense displacement networks. In: Shen, D., Liu, T., Peters, T.M., Staib, L.H., Essert, C., Zhou, S., Yap, P.T., Khan, A. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2019\. pp. 50–58. Springer International Publishing, Cham (2019)
  * [7] Hering, A., Hansen, L., Mok, T., Chung, A., Siebert, H., Häger, S., Lange, A., Kuckertz, S., Heldmann, S., Shao, W., Vesal, S., Rusu, M., Sonn, G., Estienne, T., Vakalopoulou, M., Han, L., Huang, Y., Brudfors, M., Balbastre, Y., Heinrich, M.: Learn2Reg: comprehensive multi-task medical image registration challenge, dataset and evaluation in the era of deep learning. arXiv preprint arXiv:2112.04489 (12 2021)
  * [8] Hirose, O.: Acceleration of non-rigid point set registration with downsampling and Gaussian process regression. IEEE Transactions on Pattern Analysis and Machine Intelligence 43(8), 2858–2865 (2021). https://doi.org/10.1109/TPAMI.2020.3043769
  * [9] Hirose, O.: A Bayesian formulation of coherent point drift. IEEE Transactions on Pattern Analysis and Machine Intelligence 43(7), 2269–2286 (2021). https://doi.org/10.1109/TPAMI.2020.2971687
  * [10] Lange, T., Wenckebach, T., Lamecker, H., Seebaß, M., Hünerbein, M., Eulenstein, S., Schlag, P.M.: Registration of portal and hepatic venous phase of mr/ct data for computer-assisted liver surgery planning. International Congress Series 1281, 768–772 (05 2005). https://doi.org/10.1016/j.ics.2005.03.332
  * [11] Maiseli, B., Gu, Y., Gao, H.: Recent developments and trends in point set registration methods. Journal of Visual Communication and Image Representation 46, 95–106 (2017). https://doi.org/https://doi.org/10.1016/j.jvcir.2017.03.012, https://www.sciencedirect.com/science/article/pii/S1047320317300743
  * [12] Mok, T.C.W., Chung, A.C.S.: Large deformation diffeomorphic image registration with Laplacian pyramid networks. ArXiv abs/2006.16148 (2020)
  * [13] Myronenko, A., Song, X.: Point set registration: Coherent point drift. IEEE Transactions on Pattern Analysis and Machine Intelligence 32(12), 2262–2275 (2010). https://doi.org/10.1109/TPAMI.2010.46
  * [14] Papież, B.W., Franklin, J.M., Heinrich, M.P., Gleeson, F.V., Brady, M., Schnabel, J.A.: Gifted demons: deformable image registration with local structure-preserving regularization using supervoxels for liver applications. Journal of Medical Imaging 5, 024001 (4/2018 2018). https://doi.org/10.1117/1.JMI.5.2.024001, https://doi.org/10.1117/1.JMI.5.2.024001
  * [15] Robu, M.R., Ramalhinho, J., Thompson, S.A., Gurusamy, K.S., Davidson, B.R., Hawkes, D.J., Stoyanov, D., Clarkson, M.J.: Global rigid registration of ct to video in laparoscopic liver surgery. International Journal of Computer Assisted Radiology and Surgery 13, 947 – 956 (2018)
  * [16] Thirion, J.P.: Image matching as a diffusion process: an analogy with maxwell’s demons. Medical Image Analysis 2(3), 243–260 (1998). https://doi.org/https://doi.org/10.1016/S1361-8415(98)80022-4, https://www.sciencedirect.com/science/article/pii/S1361841598800224
  * [17] Yang, D., Xu, D., Zhou, S.K., Georgescu, B., Chen, M., Grbic, S., Metaxas, D., Comaniciu, D.: Automatic liver segmentation using an adversarial image-to-image network. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2017 (07 2017)

Supplementary material:
A multi-organ point cloud registration algorithm
for abdominal CT registration

Figure 3: One protocol for hepatocellular carcinoma patients diagnosis and
follow up consist of acquiring CT scans of the patient liver at 4 different
phases of diffusion of a contrast agent: before the diffusion, during the
arterial phase, during the venous phase and 3minutes after the diffusion of
the contrast agent. This process is repeated regularly to check the disease
progression. This figure reports Dice scores on different structures when
registering an image of a certain phase with the image of the same phase
acquired during a follow-up scan of the same patient. In total, 31
registration were performed to compare 4 methods: no registration applied
(init), the demons registration algorithm [16] (Demons_Gaussian), the demons
registration algorithm with a bilateral filtering kernel [14] (Demons_bilat)
and MO-BCPD using the proposed interpolation scheme to obtain the deformation
on the whole image (MO_BCPD). We report results on 4 segmented structures used
by MO-BCPD to give a general trend of the performances of the proposed method
for whole volume registration. The 4 organs used by MO-BCPD are: the liver,
the spleen and both kidneys. We see that, as expected, MO-BCPD coupled with
the proposed interpolation framework better aligns those structures. Figure 4:
In this figure, we report Dice scores on 8 structures not used by MO-BCPD from
the same experiments as in figure 3: the body, the hepatic vein, the heart,
the aorta, the spinal cord, the skeleton, the inferior vena cava and the
portal vein. We see that, in general, MO-BCPD tends to achieve better or
competitive results compared to the other baseline studied. In particular, we
see that the interpolation of the MO-BCPD transformation to the whole volume
gives good results for vessels alignments such as the hepatic vein or the
portal vein compared to other baselines. This could be explain by the fact
that those vessels motion is correlated with the one of the organs of interest
used by MO-BCPD.