---
title:
  plain: "fairadapt: Causal Reasoning for Fair Data Pre-processing"
  formatted: "\\pkg{fairadapt}: Causal Reasoning for Fair Data Pre-processing"
  short: "\\pkg{fairadapt}: Fair Data Adaptation"
author:
  - name: Drago Plečko
    affiliation: ETH Zürich
    address: >
      Seminar for Statistics
      Rämistrasse 101
      CH-8092 Zurich
    email: \email{drago.plecko@stat.math.ethz.ch}
  - name: Nicolas Bennett
    affiliation: ETH Zürich
    address: >
      Seminar for Statistics
      Rämistrasse 101
      CH-8092 Zurich
    email: \email{nicolas.bennett@stat.math.ethz.ch}
  - name: Nicolai Meinshausen
    affiliation: ETH Zürich
    address: >
      Seminar for Statistics
      Rämistrasse 101
      CH-8092 Zurich
    email: \email{meinshausen@stat.math.ethz.ch}
abstract: >
  Machine learning algorithms are useful for various predictions tasks, but they can also learn how to discriminate, based on gender, race or other sensitive attributes. This realization gave rise to the field of fair machine learning, which aims to measure and mitigate such algorithmic bias. This manuscript describes the \proglang{R}-package \pkg{fairadapt}, which implements a causal inference pre-processing method. By making use of a causal graphical model and the observed data, the method can be used to address hypothetical questions of the form "What would my salary have been, had I been of a different gender/race?". Such individual level counterfactual reasoning can help eliminate discrimination and help justify fair decisions. We also discuss appropriate relaxations which assume certain causal pathways from the sensitive attribute to the outcome are not discriminatory.
keywords:
  formatted: [algorithmic fairness, causal inference, machine learning]
  plain: [algorithmic fairness, causal inference, machine learning]
preamble: >
  \usepackage{amsmath}
  \usepackage[ruled]{algorithm2e}
  \usepackage{bbm}
  \usepackage{array}
  \usepackage{enumerate}
  \usepackage{booktabs}
  \usepackage{makecell}
  \usepackage{threeparttablex}
  \newtheorem{definition}{Definition}
  \newcommand{\pa}{\mathrm{pa}}
  \newcommand{\Pa}{\mathrm{Pa}}
  \newcommand{\de}{\mathrm{de}}
  \newcommand{\ch}{\mathrm{ch}}
  \newcommand{\an}{\mathrm{an}}
  \newcommand{\pr}{\mathbbm{P}}
  \newcommand{\ex}{\mathbbm{E}}
  \renewcommand{\tilde}[1]{ {#1}^{(fp)}}
  \def\ci{{\perp\!\!\!\perp}}
vignette: >
  %\VignetteIndexEntry{Fair Data Adaptation (Plecko et al., JSS 2021)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output: >
  if (packageVersion("rticles") < 0.5 || rmarkdown::pandoc_version() >= 2)
    rticles::jss_article else rmarkdown::html_vignette
documentclass: jss
classoption: nojss
bibliography: jss.bib
pkgdown:
  as_is: true
  extension: pdf
---

```{r, setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(fairadapt)
library(data.table)
library(ggplot2)
library(ggraph)
options(tinytex.verbose = TRUE)

options(
  prompt = "R> ",
  continue = "+ ",
  width = 70,
  useFancyQuotes = FALSE
)
```

```{tikz, setup-graph, eval = FALSE, echo = FALSE}
\tikzset{
  >=stealth,
  rv/.style = {circle, draw, thick, minimum size = 6mm},
  rvc/.style = {triangle, draw, thick, minimum size = 10mm},
  node distance = 18mm
}
```

# Introduction

Machine learning algorithms have become prevalent tools for decision-making in socially sensitive situations, such as determining credit-score ratings or predicting recidivism during parole. It has been recognized that algorithms are capable of learning societal biases, for example with respect to race \citep{larson2016recidivism} or gender \citep{lambrecht2019algorithmic, blau2003pay}, and this realization seeded an important debate in the machine learning community about fairness of algorithms and their impact on decision-making.

In order to define and measure discrimination, existing intuitive notions have been statistically formalized, thereby providing fairness metrics. For example, *demographic parity* \citep{darlington1971fairness} requires the protected attribute $A$ (gender/race/religion etc.) to be independent of a constructed classifier or regressor $\widehat{Y}$, written as $\widehat{Y} \ci A$. Another notion, termed *equality of odds* \citep{hardt2016eosl}, requires equal false positive and false negative rates of classifier $\widehat{Y}$ between different groups (females and males for example), written as $\widehat{Y} \ci A \mid Y$. To this day, various different notions of fairness exist, which are sometimes incompatible \citep{corbett2018measure}, meaning not of all of them can be achieved for a predictor $\widehat{Y}$ simultaneously. There is still no consensus on which notion of fairness is the correct one.

The discussion on algorithmic fairness is, however, not restricted to the machine learning domain. There are many legal and philosophical aspects that have arisen. For example, the legal distinction between disparate impact and disparate treatment \citep{mcginley2011ricci} is important for assessing fairness from a judicial point of view. This in turn emphasizes the importance of the interpretation behind the decision-making process, which is often not the case with black-box machine learning algorithms. For this reason, research in fairness through a causal inference lens has gained attention.

A possible approach to fairness is the use of counterfactual reasoning \citep{galles1998axiomatic}, which allows for arguing what might have happened under different circumstances that never actually materialized, thereby providing a tool for understanding and quantifying discrimination. For example, one might ask how a change in sex would affect the probability of a specific candidate being accepted for a given job opening. This approach has motivated another notion of fairness, termed *counterfactual fairness* \citep{kusner2017counterfactual}, which states that the decision made, should remain fixed, even if, hypothetically, some parameters such as race or gender were to be changed (this can be written succinctly as $\widehat{Y}(a) = \widehat{Y}(a')$ in the potential outcomes notation). Causal inference can also be used for decomposition of the parity gap measure \citep{zhang2018fairness}, $\pr(\widehat{Y} = 1 \mid A = a) - \pr(\widehat{Y} = 1 \mid A = a')$, into the direct, indirect and spurious components (yielding further insights into the demographic parity as a criterion), as well as the introduction of so-called resolving variables \cite{kilbertus2017avoiding}, in order to relax the possibly prohibitively strong notion of demographic parity.

The following sections describe an implementation of the fair data adaptation method outlined in \cite{plecko2020fair}, which combines the notions of counterfactual fairness and resolving variables, and explicitly computes counterfactual values for individuals. The implementation is available as \proglang{R}-package \pkg{fairadapt} from CRAN. Currently there are only few packages distributed via CRAN that relate to fair machine learning. These include \pkg{fairml} \citep{scutari2021fairml}, which implements the non-convex method of \cite{komiyama2018nonconvex}, as well as packages \pkg{fairness} \citep{kozodoi2021fairness} and \pkg{fairmodels} \citep{wisniewski2021fairmodels}, which serve as diagnostic tools for measuring algorithmic bias and provide several pre- and post-processing methods for bias mitigation. The only causal method, however, is presented by \pkg{fairadapt}. Even though theory in fair machine learning is being expanded at an accelerating pace, good quality implementations of the developed methods are often not available.

The rest of the manuscript is organized as follows: In Section \ref{methodology} we describe the methodology behind \pkg{fairadapt}, together with quickly reviewing some important concepts of causal inference. In Section \ref{implementation} we discuss implementation details and provide some general user guidance, followed by Section \ref{illustration}, which illustrates the usage of \pkg{fairadapt} through a large, real-world dataset and a hypothetical fairness application. Finally, in Section \ref{extensions} we elaborate on some extensions, such as Semi-Markovian models and resolving variables.

# Methodology

First, the intuition behind \pkg{fairadapt} is described using an example, followed by a more rigorous mathematical formulation, using Markovian structural causal models (SCMs). Some relevant extensions, such as the Semi-Markovian case and the introduction of so called *resolving variables*, are discussed in Section \ref{extensions}.

## University Admission Example

Consider the example of university admission based on previous educational achievement and an admissions test. Variable $A$ is the protected attribute, describing candidate gender, with $A = a$ corresponding to females and $A = a'$ to males. Furthermore, let $E$ be educational achievement (measured for example by grades achieved in school) and $T$ the result of an admissions test for further education. Finally, let $Y$ be the outcome of interest (final score) upon which admission to further education is decided. Edges in the graph in Figure \ref{fig:uni-adm} indicate how variables affect one another.

```{tikz, uni-adm, fig.cap = "University admission based on previous educational achievement $E$ combined with and an admissions test score $T$. The protected attribute $A$, encoding gender, has an unwanted causal effect on  $E$, $T$, as well as $Y$, which represents the final score used for the admission decision.", fig.ext = "png", cache = TRUE, echo = FALSE, out.width = "50%"}

<<setup-graph>>

\begin{tikzpicture}

  \pgfsetarrows{latex-latex};

  \begin{scope}
    \node[rv] (1) at (-2,0) {$A$};
    \node[rv] (2) at (0,0) {$E$};
    \node[rv] (3) at (2,0) {$T$};
    \node[rv] (4) at (4,0) {$Y$};
    \draw[->] (1) -- (2);
    \draw[->] (1) edge[bend left = 20] (3);
    \draw[->] (2) -- (3);
    \draw[->] (2) -- (3);
    \draw[->] (3) -- (4);
    \draw[->] (2) edge[bend right = 25] (4);
  \end{scope}

\end{tikzpicture}
```

Attribute $A$, gender, has a causal effect on variables $E$, $T$, as well as $Y$, and we wish to eliminate this effect. For each individual with observed values $(a, e, t, y)$ we want to find a mapping

$$(a, e, t, y) \longrightarrow  (\tilde{a}, \tilde{e}, \tilde{t}, \tilde{y}),$$

which represents the value the person would have obtained in an alternative world where everyone was female. Explicitly, to a male person with education value $e$, we assign the transformed value $\tilde{e}$ chosen such that

$$\pr(E \geq e \mid A = a') = \pr(E \geq \tilde{e} \mid A = a).$$

The key idea is that the *relative educational achievement within the subgroup* remains constant if the protected attribute gender is modified. If, for example, a male has a higher educational achievement value than 70\% of males in the dataset, we assume that he would also be better than 70\% of females had he been female^[This assumption of course is not empirically testable, as it is impossible to observe both a female and a male version of the same individual.]. After computing transformed educational achievement values corresponding to the *female* world ($\tilde{E}$), the transformed test score values $\tilde{T}$ can be calculated in a similar fashion, but conditioned on educational achievement. That is, a male with values $(E, T) = (e, t)$ is assigned a test score $\tilde{t}$ such that

$$\pr(T \geq t \mid E = e, A = a') = \pr(T \geq \tilde{t} \mid E = \tilde{e}, A = a),$$

where the value $\tilde{e}$ was obtained in the previous step. This step can be visualized as shown in Figure \ref{fig:rel-edu}.

```{tikz, rel-edu, fig.cap = "A graphical visualization of the quantile matching procedure. Given a male with a test score corresponding to the 70\\% quantile, we would hypothesize, that if the gender was changed, the individual would have achieved a test score corresponding to the 70\\% quantile of the female distribution.", fig.ext = "png", cache = TRUE, echo = FALSE, engine.opts=list(extra.preamble = c("\\usepackage{bbm}", "\\usepackage{pgfplots}", "\\pgfplotsset{compat=newest}"))}

\newcommand{\pr}{\mathbbm{P}}
\renewcommand{\tilde}[1]{ {#1}^{(fp)}}

\usetikzlibrary{arrows.meta}

\pgfmathdeclarefunction{gauss}{2}{\pgfmathparse{1/(#2*sqrt(2*pi))*exp(-((x-#1)^2)/(2*#2^2))}}

\begin{tikzpicture}

  \begin{axis}[
    no markers, domain=0:10, samples=100,
    axis lines = left,
    xlabel = $t$,
    ylabel = $\pr(t)$,
    every axis y label/.style={at=(current axis.above origin),anchor=south},
    every axis x label/.style={at=(current axis.right of origin),anchor=west},
    height=5cm, width=12cm,
    xtick=\empty, ytick=\empty,
    enlargelimits=false, clip=false, axis on top,
    grid = major
  ]

    \addplot [very thick,green!50!black] {gauss(4,1)};
    \addplot [very thick,blue!50!black] {gauss(6.5,0.8)};

    \draw[-{Latex[length=3mm,width=2mm]}, dashed] (axis cs:2.718,0.175) to[bend right = 30] (axis cs:5.474,0.219);
    \draw[-{Latex[length=3mm,width=2mm]}, dashed] (axis cs:4.524, 0.348) to[bend right = 30] (axis cs:6.920, 0.435);

    \node at (axis cs:2.718,0.175) [above, left]
      {\shortstack{$10\%$-quantile \\ (Male)}};
    \node at (axis cs:5.474,0.219) [below = 0.3cm, right = 0cm]
      {\shortstack{$10\%$-quantile \\ (Female)}};
    \node at (axis cs:4.524, 0.348) [above = 0.3cm]
      {\shortstack{$70\%$-quantile \\ (Male)}};
    \node at (axis cs:6.920, 0.435) [right]
      {\shortstack{$70\%$-quantile \\ (Female)}};
    \node at (axis cs:4,0.5) [below = 0.65cm, left = 0.5cm, green!50!black]
      {$T \mid E = e, A = a'$};
    \node at (axis cs:6.5,0.5) [below = 1.3cm, right = 1cm, blue!50!black]
      {$T \mid E = \tilde{e}, A = a$
    };

  \end{axis}

\end{tikzpicture}
```

As a final step, the outcome variable $Y$ remains to be adjusted. The adaptation is based on the same principle as above, using transformed values of both education and the test score. The resulting value $\tilde{y}$ of $Y = y$ satisfies

\begin{equation} \label{eq:labeltransform}
	\pr(Y \geq y \mid E = e, T = t, A = a') = \pr(Y \geq \tilde{y} \mid E = \tilde{e}, T = \tilde{t}, A = a).
\end{equation}

This form of counterfactual correction is known as *recursive substitution* \citep[Chapter~7]{pearl2009causality} and is described more formally in the following sections. The reader who is satisfied with the intuitive notion provided by the above example is encouraged to go straight to Section \ref{implementation}.

## Structural Causal Models

In order to describe the causal mechanisms of a system, a *structural causal model* (SCM) can be hypothesized, which fully encodes the assumed data-generating process. An SCM is represented by a 4-tuple $\langle V, U, \mathcal{F}, \pr(u) \rangle$, where

* $V = \lbrace V_1, \ldots, V_n \rbrace$ is the set of observed (endogenous) variables.
* $U = \lbrace U_1, \ldots, U_n \rbrace$ are latent (exogenous) variables.
* $\mathcal{F} = \lbrace f_1, \ldots, f_n \rbrace$ is the set of functions determining $V$, $v_i \gets f_i(\pa(v_i), u_i)$, where $\pa(V_i) \subset V, U_i \subset U$ are the functional arguments of $f_i$ and $\pa(V_i)$ denotes the parent vertices of $V_i$.
* $\pr(u)$ is a distribution over the exogenous variables $U$.

Any particular SCM is accompanied by a graphical model $\mathcal{G}$ (a directed acyclic graph), which summarizes which functional arguments are necessary for computing the values of each $V_i$ and therefore, how variables affect one another. We assume throughout, without loss of generality, that

(i) $f_i(\pa(v_i), u_i)$ is increasing in $u_i$ for every fixed $\pa(v_i)$.
(ii) Exogenous variables $U_i$ are uniformly distributed on $[0, 1]$.

In the following section, we discuss the Markovian case in which all exogenous variables $U_i$ are mutually independent. The Semi-Markovian case, where variables $U_i$ are allowed to have a mutual dependency structure, alongside the extension introducing *resolving variables*, are discussed in Section \ref{extensions}.

## Markovian SCM Formulation

Let $Y$ take values in $\mathbbm{R}$ and represent an outcome of interest and $A$ be the protected attribute taking two values $a, a'$. The goal is to describe a pre-processing method which transforms the entire data $V$ into its fair version $\tilde{V}$. This can be achieved by computing the counterfactual values $V(A = a)$, which would have been observed if the protected attribute was fixed to a baseline value $A = a$ for the entire sample.

More formally, going back to the *university admission* example above, we want to align the distributions

$$V_i \mid \pa(V_i), A = a \text{ and } V_i \mid \pa(V_i), A = a',$$

meaning that the distribution of $V_i$ should be indistinguishable for both female and male applicants, for every variable $V_i$. Since each function $f_i$ of the original SCM is reparametrized so that $f_i(\pa(v_i), u_i)$ is increasing in $u_i$ for every fixed $\pa(v_i)$, and also due to variables $U_i$ being uniformly distributed on $[0, 1]$, variables $U_i$ can be seen as the latent *quantiles*.

The algorithm proposed for data adaption proceeds by fixing $A = a$, followed by iterating over descendants of the protected attribute $A$, sorted in topological order. For each $V_i$, the assignment function $f_i$ and the corresponding quantiles $U_i$ are inferred. Finally, transformed values $\tilde{V_i}$ are obtained by evaluating $f_i$, using quantiles $U_i$ and the transformed parents $\tilde{\pa(V_i)}$ (see Algorithm \ref{algo:fairadapt}).

\begin{algorithm}
	\DontPrintSemicolon
	\KwIn{$V$, causal graph $\mathcal{G}$}
	set $A \gets a$ for everyone\\
	\For{$V_i \in \de(A)$ in topological order}{
	  learn function $V_i \gets f_i(\pa(V_i), U_i)$ \;
		infer quantiles $U_i$ associated with the variable $V_i$\;
		transform values as $\tilde{V_i} \gets f_i (\tilde{\pa(V_i)}, U_i)$
		 \;
  }
  \Return{$\tilde{V}$}
	\caption{Fair Data Adaptation}
	\label{algo:fairadapt}
\end{algorithm}

The assignment functions $f_i$ of the SCM are of course unknown, but are non-parametrically inferred at each step. Algorithm \ref{algo:fairadapt} obtains the counterfactual values $V(A = a)$ under the $do(A = a)$ intervention for each individual, while keeping the latent quantiles $U$ fixed. In the case of continuous variables, the latent quantiles $U$ can be determined exactly, while for the discrete case, this is more subtle and described in \citet[Section~5]{plecko2020fair}.

# Implementation

In order to perform fair data adaption using the \pkg{fairadapt} package, the function `fairadapt()` is exported, which returns an object of class `fairadapt`. Implementations of the base \proglang{R} S3 generics `print()`, `plot()` and `predict()`, as well as the generic `autoplot()`, exported from \pkg{ggplot2} \citep{wickham2016ggplot2}, are provided for `fairadapt` objects, alongside `fairadapt`-specific implementations of S3 generics `visualizeGraph()`, `adaptedData()` and `fairTwins()`. Finally, an extension mechanism is available via the S3 generic function `computeQuants()`, which is used for performing the quantile learning step.

The following sections show how the listed methods relate to one another alongside their intended use, beginning with constructing a call to `fairadapt()`. The most important arguments of `fairadapt()` include:

* `formula`: Argument of type `formula`, specifying the dependent and explanatory variables.
* `adj.mat`: Argument of type `matrix`, encoding the adjacency matrix.
* `train.data` and `test.data`: Both of type `data.frame`, representing the respective datasets.
* `prot.attr`: Scalar-valued argument of type `character` identifying the protected attribute.

As a quick demonstration of fair data adaption using, we load the `uni_admission` dataset provided by \pkg{fairadapt}, consisting of synthetic university admission data of `r nrow(uni_admission)` students. We subset this data, using the first `n_samp` rows as training data (`uni_trn`) and the following `n_samp` rows as testing data (`uni_tst`). Furthermore, we construct an adjacency matrix `uni_adj` with edges $\texttt{gender} \to \texttt{edu}$, $\texttt{gender} \to \texttt{test}$, $\texttt{edu} \to \texttt{test}$, $\texttt{edu} \to \texttt{score}$, and $\texttt{test} \to \texttt{score}$. As the protected attribute, we choose `gender`.

```{r, basic}
n_samp <- 200

uni_dat <- data("uni_admission", package = "fairadapt")
uni_dat <- uni_admission[seq_len(2 * n_samp), ]

head(uni_dat)

uni_trn <- head(uni_dat, n = n_samp)
uni_tst <- tail(uni_dat, n = n_samp)

uni_dim <- c(       "gender", "edu", "test", "score")
uni_adj <- matrix(c(       0,     0,      0,       0,
                            1,     0,      0,       0,
                            1,     1,      0,       0,
                            0,     1,      1,       0),
                  ncol = length(uni_dim),
                  dimnames = rep(list(uni_dim), 2))

basic <- fairadapt(score ~ ., train.data = uni_trn,
                    test.data = uni_tst, adj.mat = uni_adj,
                    prot.attr = "gender")

basic
```

The implicitly called `print()` method in the previous code block displays some information about how `fairadapt()` was called, such as number of variables, the protected attribute and also the total variation before and after adaptation, defined as

$$\ex [Y \mid A = a] - \ex [Y \mid A = a'] \text{ and } \ex [\tilde{Y} \mid A = a] - \ex [\tilde{Y} \mid A = a'],$$

respectively, where $Y$ denotes the outcome variable. Total variation in the case of a binary outcome $Y$, corresponds to the parity gap.

## Specifying the Graphical Model

As the algorithm used for fair data adaption in `fairadapt()` requires access to the underlying graphical causal model $\mathcal{G}$ (see Algorithm \ref{algo:fairadapt}), a corresponding adjacency matrix can be passed as `adj.mat` argument. The convenience function `graphModel()` turns a graph specified as an adjacency matrix into an annotated graph using the \pkg{igraph} package \citep{csardi2006igraph}. While exported for the user to invoke manually, this function is called as part of the `fairadapt()` routine and the resulting `igraph` object can be visualized by calling the S3 generic `visualizeGraph()`, exported from `fairadapt` on an object of class `fairadapt`.

```{r, graph-model}
uni_graph <- graphModel(uni_adj)
```

```{r, graph-plot, echo = FALSE, fig.width = 6, fig.height = 3, fig.cap = "The underlying graphical model corresponding to the university admission example (also shown in Figure \\ref{fig:uni-adm}).", out.width = "60%"}

set.seed(11)

ggraph(graphModel(uni_adj), "igraph", algorithm = "sugiyama") +
  geom_edge_link(arrow = arrow(length = unit(4, "mm"), angle = 15,
                               type = "closed"),
                 end_cap = circle(8, "mm"),
                 color = "grey20") +
  geom_node_point(color = "grey80", size = 21) +
  geom_node_text(aes(x = x, y = y, label = name), size = 5) +
  theme_bw() +
  theme(plot.margin = margin(30, 30, 30, 30), panel.border = element_blank(),
        panel.background = element_blank(), axis.title = element_blank(),
        axis.text = element_blank(), axis.line = element_blank(),
        axis.ticks = element_blank(), panel.grid = element_blank()) +
  coord_cartesian(clip = "off")
```

A visualization of the `igraph` object returned by `graphModel()` is available from Figure \ref{fig:graph-plot}. The graph shown is equivalent to that of Figure \ref{fig:uni-adm} as they both represent the same causal model.

## Quantile Learning Step

The training step in `fairadapt()` can be carried out in two slightly distinct ways: Either by specifying training and testing data simultaneously, or by only passing training data (and at a later stage calling `predict()` on the returned `fairadapt` object in order to perform data adaption on new test data). The advantage of the former option is that the quantile regression step is performed on a larger sample size, which can lead to more precise inference in practice.

The two data frames passed as `train.data` and `test.data` are required to have column names which also appear in the row and column names of the adjacency matrix, alongside the protected attribute $A$, passed as scalar-valued character vector `prot.attr`.

The quantile learning step of Algorithm \ref{algo:fairadapt} can in principle be carried out by several methods, three of which are implemented in \pkg{fairadapt}:

* Quantile Regression Forests \citep{meinshausen2006qrf, wright2015ranger}.
* Non-crossing quantile neural networks \citep{cannon2018non, cannon2015package}.
* Linear Quantile Regression \citep{koenker2001qr, koenker2018package}.

Using linear quantile regression is the most efficient option in terms of runtime, while for non-parametric models and mixed data, the random forest approach is well-suited, at the expense of a slight increase in runtime. The neural network approach is, substantially slower when compared to linear and random forest estimators and consequently does not scale well to large sample sizes. As default, the random forest based approach is implemented, due to its non-parametric nature and computational speed. However, for smaller sample sizes, the neural network approach can also demonstrate competitive performance. A quick summary outlining some differences between the three natively supported methods is available from Table \ref{tab:qmethods}.

\begin{table}
\centering
\begin{threeparttable}
\begin{tabular}[t]{llll}
\toprule
  & Random Forests & Neural Networks & Linear Regression\\
\midrule
    \proglang{R}-package & \pkg{ranger} & \pkg{qrnn} & \pkg{quantreg} \\
    \addlinespace[0.3em]
    \texttt{quant.method} & \code{rangerQuants} & \code{mcqrnnQuants} & \code{linearQuants} \\
    \addlinespace[0.3em]
    complexity & $O(np\log n)$ & $O(npn_{\text{epochs}})$ & $O(p^2n)$ \\
    \addlinespace[0.3em]
    \makecell[l]{default\\parameters} & \makecell[l]{$ntrees = 500$\\$mtry = \sqrt{p}$} & \makecell[l]{2-layer fully\\connected\\feed-forward\\network} & \makecell[l]{\code{"br"} method of\\Barrodale and\\Roberts used for\\fitting} \\
    \addlinespace[0.3em]
    $T(200, 4)$ & $0.4$ sec & $89$ sec & $0.3$ sec \\
    \addlinespace[0.3em]
    $T(500, 4)$ & $0.9$ sec & $202$ sec & $0.5$ sec \\
\bottomrule
\end{tabular}
\caption{Summary table of different quantile regression methods. $n$ is the number of samples, $p$ number of covariates, $n_{\text{epochs}}$ number of training epochs for the neural network. $T(n, 4)$ denotes the runtime of different methods on the university admission dataset, with $n$ training and testing samples.}
  \label{tab:qmethods}
\end{threeparttable}
\end{table}

The above set of methods is not exhaustive. Further options are conceivable and therefore \pkg{fairadapt} provides an extension mechanism to account for this. The `fairadapt()` argument `quant.method` expects a function to be passed, a call to which will be constructed with three unnamed arguments:

1. A `data.frame` containing data to be used for quantile regression. This will either be the `data.frame` passed as `train.data`, or depending on whether `test.data` was specified, a row-bound version of train and test datasets.
1. A logical flag, indicating whether the protected attribute is the root node of the causal graph. If the attribute $A$ is a root node, we know that $\pr(X \mid \text{do}(A = a)) = \pr(X \mid A = a)$. Therefore, the interventional and conditional distributions are in this case the same, and we can leverage this knowledge in the quantile learning procedure, by splitting the data into $A = 0$ and $A = 1$ groups.
1. A `logical` vector of length `nrow(data)`, indicating which rows in the `data.frame` passed as `data` correspond to samples with baseline values of the protected attribute.

Arguments passed as `...` to `fairadapt()` will be forwarded to the function specified as `quant.method` and passed after the first three fixed arguments listed above. The return value of the function passed as `quant.method` is expected to be an S3-classed object. This object should represent the conditional distribution $V_i \mid \pa(V_i)$ (see function `rangerQuants()` for an example). Additionally, the object should have an implementation of the S3 generic function `computeQuants()` available. For each row $(v_i, \pa(v_i))$ of the `data` argument, the `computeQuants()` function uses the S3 object to (i) infer the quantile of $v_i \mid \pa(v_i)$; (ii) compute the counterfactual value $\tilde{v}_i$ under the change of protected attribute, using the counterfactual values of parents $\pa(\tilde{v}_i)$ computed in previous steps (values $\pa(\tilde{v}_i)$ are contained in the `newdata` argument). For an example, see `computeQuants.ranger()` method. 

## Fair-Twin Inspection

The university admission example presented in Section \ref{methodology} demonstrates how to compute counterfactual values for an individual while preserving their relative educational achievement. Setting candidate gender as the protected attribute and gender level *female* as baseline value, for a *male* student with values $(a, e, t, y)$, his *fair-twin* values $(\tilde{a}, \tilde{e}, \tilde{t}, \tilde{y})$, i.e., the values the student would have obtained, had he been *female*, are computed. These values can be retrieved from a `fairadapt` object by calling the S3-generic function `fairTwins()` as:

```{r, fairtwin-uni}
ft_basic <- fairTwins(basic, train.id = seq_len(n_samp))
head(ft_basic, n = 3)
```

In this example, we compute the values in a *female* world. Therefore, for *female* applicants, the values remain fixed, while for *male* applicants the values are adapted, as can be seen from the output.

# Illustration

As a hypothetical real-world use of \pkg{fairadapt}, suppose that after a legislative change the US government has decided to adjust the salary of all of its female employees in order to remove both disparate treatment and disparate impact effects. To this end, the government wants to compute the counterfactual salary values of all female employees, that is the salaries that female employees would obtain, had they been male.

To do this, the government is using data from the 2018 American Community Survey by the US Census Bureau, available in pre-processed form as a package dataset from \pkg{fairadapt}. Columns are grouped into demographic (`dem`), familial (`fam`), educational (`edu`) and occupational (`occ`) categories and finally, salary is selected as response (`res`) and sex as the protected attribute (`prt`):

```{r, load-census}
gov_dat <- data("gov_census", package = "fairadapt")
gov_dat <- get(gov_dat)

head(gov_dat)

dem <- c("age", "race", "hispanic_origin", "citizenship",
         "nativity", "economic_region")
fam <- c("marital", "family_size", "children")
edu <- c("education_level", "english_level")
occ <- c("hours_worked", "weeks_worked", "occupation",
         "industry")

prt <- "sex"
res <- "salary"
```

The hypothesized causal graph for the dataset is given in Figure \ref{fig:census-tikz}. According to this, the causal graph can be specified as an adjacency matrix `gov_adj` and as confounding matrix `gov_cfd`:

```{tikz, census-tikz, fig.cap = "The causal graph for the government-census dataset. $D$ are demographic features, $A$ is gender, $F$ represents marital and family information, $E$ education, $W$ work-related information and $Y$ the salary, which is also the outcome of interest.", fig.ext = "png", cache = TRUE, echo = FALSE, out.width = "50%"}

<<setup-graph>>

\begin{tikzpicture}

  \pgfsetarrows{latex-latex};

  \begin{scope}
    \node[rv] (c) at (2,2) {$D$};
    \node[rv] (a) at (-2,2) {$A$};
    \node[rv] (m) at (-3,0) {$F$};
    \node[rv] (l) at (-1,0) {$E$};
    \node[rv] (r) at (1,0) {$W$};
    \node[rv] (y) at (3,0) {$Y$};
    \draw[->] (c) -- (m);
    \draw[->] (c) -- (l);
    \draw[->] (c) -- (r);
    \draw[->] (c) -- (y);
    \draw[->] (a) -- (m);
    \draw[->] (m) -- (l);
    \draw[->] (l) -- (r);
    \draw[->] (r) -- (y);
    \path[->] (a) edge[bend left = 0] (l);
    \path[->] (a) edge[bend left = 0] (r);
    \path[->] (a) edge[bend left = 0] (y);
    \path[->] (m) edge[bend right = 20] (r);
    \path[->] (m) edge[bend right = 30] (y);
    \path[->] (r) edge[bend right = 20] (y);
    \path[->, dashed] (a) edge[bend left = 10] (c);
  \end{scope}

\end{tikzpicture}
```

```{r, census-adj}
cols <- c(dem, fam, edu, occ, prt, res)

gov_adj <- matrix(0, nrow = length(cols), ncol = length(cols),
                  dimnames = rep(list(cols), 2))
gov_cfd <- gov_adj

gov_adj[dem, c(fam, edu, occ, res)] <- 1
gov_adj[fam, c(     edu, occ, res)] <- 1
gov_adj[edu, c(          occ, res)] <- 1
gov_adj[occ,                  res ] <- 1

gov_adj[prt, c(fam, edu, occ, res)] <- 1

gov_cfd[prt, dem] <- 1
gov_cfd[dem, prt] <- 1

gov_grph <- graphModel(gov_adj, gov_cfd)
```

```{r, census-graph, echo = FALSE, fig.height = 9, fig.width = 9, fig.cap = "Full causal graph for the government census dataset, expanding the grouped view presented in Figure \\ref{fig:census-tikz}. \\textit{Demographic} features include age (\\textbf{ag}), race (\\textbf{ra}), whether an employee is of Hispanic origin (\\textbf{hi}), is US citizen (\\textbf{ci}), whether the citizenship is native (\\textbf{na}), alongside the corresponding economic region (\\textbf{ec}). \\textit{Familial} features are marital status (\\textbf{ma}), family size (\\textbf{fa}) and number of children (\\textbf{ch}), \\textit{educational} features include education (\\textbf{ed}) and English language levels (\\textbf{en}), and \\textit{occupational} features, weekly working hours (\\textbf{ho}), yearly working weeks (\\textbf{we}), job (\\textbf{oc}) and industry identifiers (\\textbf{in}). Finally, the yearly salary (\\textbf{sa}) is used as the \\textit{response} variable and employee sex (\\textbf{se}) as the \\textit{protected} attribute variable.", out.width = "90%"}
set.seed(11)
gov_tmp <- graphModel(
  `dimnames<-`(gov_adj, lapply(dimnames(gov_adj), substr, 1L, 2L)),
  `dimnames<-`(gov_cfd, lapply(dimnames(gov_cfd), substr, 1L, 2L))
)

grp <- list(dem, fam, edu, occ, prt, res)
grp <- `names<-`(
  rep(c("demographic", "familial", "educational", "occupational",
        "protected", "response"), lengths(grp)),
  substr(unlist(grp), 1L, 2L)
)

igraph::V(gov_tmp)$color <- grp[names(igraph::V(gov_tmp))]

gov_tmp <- igraph::delete_edges(
  gov_tmp, which(igraph::E(gov_tmp)$lty == "blank")
)

lty_selector <- function(lty) {
  function(layout) {
    get_all <- get_edges()
    edges <- get_all(layout)
    res <- edges[edges$lty == lty, ]
    res
  }
}

ggraph(gov_tmp, "igraph", algorithm = "fr") +
  geom_edge_arc(data = lty_selector("dashed"),
                arrow = arrow(length = unit(4, "mm"), angle = 15,
                              ends = "both", type = "closed"),
                start_cap = circle(6.5, "mm"),
                end_cap = circle(6.5, "mm"),
                strength = 0.25,
                color = "grey20",
                linetype = "dashed") +
  geom_edge_link(data = lty_selector("solid"),
                 arrow = arrow(length = unit(4, "mm"), angle = 15,
                               type = "closed"),
                 end_cap = circle(6.5, "mm"),
                 color = "grey20",
                 linetype = "solid") +
  geom_node_point(aes(color = color), size = 15) +
  geom_node_text(aes(x = x, y = y, label = name), size = 5) +
  scale_edge_linetype_manual(values = c(dashed = "dashed", solid = "solid")) +
  scale_color_discrete(name = "Grouping") +
  theme_bw() +
  theme(plot.margin = margin(30, 30, 30, 30), legend.position = "bottom",
        panel.background = element_blank(), axis.title = element_blank(),
        axis.text = element_blank(), axis.line = element_blank(),
        axis.ticks = element_blank(), panel.grid = element_blank(),
        panel.border = element_blank()) +
  coord_cartesian(clip = "off")
```

Before applying `fairadapt()`, we first log-transform the salaries and look at respective densities by sex group. We subset the data by using `n_samp` samples for training and `n_pred` samples for predicting and plot the data before performing the adaption.

```{r, log-sub}
gov_dat$salary <- log(gov_dat$salary)

n_samp <- 30000
n_pred <- 5

gov_trn <- head(gov_dat, n = n_samp)
gov_prd <- tail(gov_dat, n = n_pred)
```

```{r, before-adapt, echo = FALSE, fig.height = 3, fig.cap = "Visualization of salary densities grouped by employee sex, indicating a shift to higher values for male employees. This uses the US government-census dataset and shows the situation before applying fair data adaption, while Figure \\ref{fig:vis-adapt} presents transformed salary data."}
ggplot(gov_trn, aes(x = salary, fill = sex)) +
  geom_density(alpha = 0.4)  +
  theme_bw() +
  ggtitle("Salary density by gender")
```

There is a clear shift between the two sexes, indicating that `male` employees are currently better compensated when compared to `female` employees. However, this differences in `salary` could, in principle, be attributed to factors apart form gender inequality, such as the economic region in which an employee works. This needs to be accounted for as well, i.e., we do not wish to remove differences in salary between economic regions.

```{r, census-adapt}
gov_ada <- fairadapt(salary ~ ., train.data = gov_trn,
                     adj.mat = gov_adj, prot.attr = prt)
```

After performing the adaptation, we can investigate whether the salary gap has shrunk. The densities after adaptation can be visualized using the \pkg{ggplot2}-exported S3 generic function `autoplot()`:

```{r, vis-adapt, fig.height = 3, fig.cap = "The salary gap between male and female employees of the US government according to the government-census dataset is clearly reduced when comparing raw data (see Figure \\ref{fig:before-adapt}) to transformed salary data as yielded by applying fair data adaption using employee sex as the protected attribute and assuming a causal graph as shown in Figure \\ref{fig:census-graph}."}
autoplot(gov_ada, when = "after") +
  theme_bw() +
  ggtitle("Adapted salary density by gender")
```

If we are provided with additional testing data, and wish to adapt this as well, we can use the base \proglang{R} S3 generic function `predict()`:

```{r, census-predict}
predict(gov_ada, newdata = gov_prd)
```

Finally, we can do fair-twin inspection using the `fairTwins()` function of \pkg{fairadapt}, to retrieve counterfactual feature values:

```{r, census-twins}
fairTwins(gov_ada, train.id = 1:5,
          cols = c("sex", "age", "education_level", "salary"))
```

Values are unchanged for female individuals (as *female* is used as baseline level), as is the case for `age`, which is not a descendant of the protected attribute `sex` (see Figure \ref{fig:census-graph}). However, variables `education_level` and `salary` do change for males, as they are descendants of the protected attribute `sex`.

The variable `hours_worked` is also a descendant of $A$, and one might argue that this variable should *not* be adapted in the procedure, i.e., it should remain the same, irrespective of employee sex. This is the idea behind *resolving variables*, which we discuss next, in Section \ref{adding-resolving-variables}. It is worth emphasizing that we are not trying to answer the question of which choice of resolving variables is the correct one in the above example - this choice is left to social scientists deeply familiar with context and specifics of the above described dataset.

# Extensions

Several extensions to the basic Markovian SCM formulation introduced in Section \ref{markovian-scm-formulation} exist, some of which are available for use in `fairadapt()` and are outlined in the following sections.

## Adding Resolving Variables

\cite{kilbertus2017avoiding} discuss that in some situations the protected attribute $A$ can affect variables in a non-discriminatory way. For instance, in the Berkeley admissions dataset \citep{bickel1975sex} we observe that females often apply for departments with lower admission rates and consequently have a lower admission probability. However, we perhaps would not wish to account for this difference in the adaptation procedure, if we were to argue that applying to a certain department is a choice everybody is free to make. Such examples motivated the idea of *resolving variables* by \citet{kilbertus2017avoiding}. A variable $R$ is called resolving if

(i) $R \in \de(A)$, where $\de(A)$ are the descendants of $A$ in the causal graph $\mathcal{G}$.
(i) The causal effect of $A$ on $R$ is considered to be non-discriminatory.

In presence of resolving variables, computation of the counterfactual is carried out under the more involved intervention do$(A = a, R = R(a'))$. The potential outcome value $V(A = a, R = R(a'))$ is obtained by setting $A = a$ and computing the counterfactual while keeping the values of resolving variables to those they *attained naturally*. This is a nested counterfactual and the difference in Algorithm \ref{algo:fairadapt} is simply that resolving variables $R$ are skipped in the for-loop. In order to perform fair data adaptation with the variable `test` being resolving in the `uni_admission` dataset used in Section \ref{implementation}, the string `"test"` can be passed as `res.vars` to `fairadapt()`.

```{r, res-uni}
fairadapt(score ~ ., train.data = uni_trn, test.data = uni_tst,
          adj.mat = uni_adj, prot.attr = "gender", res.vars = "test")
```

As can be seen from the respective model summary outputs, the total variation after adaptation, in this case, is larger than in the `basic` example from Section \ref{implementation}, with no resolving variables. The intuitive reasoning here is that resolving variables allow for some discrimination, so we expect to see a larger total variation between the groups.

```{r, res-assign}
uni_res <- graphModel(uni_adj, res.vars = "test")
```

```{r, res-graph, echo = FALSE, fig.width = 6, fig.height = 3, fig.cap = "Visualization of the causal graph corresponding to the university admissions example introduced in Section \\ref{introduction} with the variable \\texttt{test} chosen as a \\textit{resolving variable} and therefore highlighted in red.", out.width = "60%"}

set.seed(11)

ggraph(graphModel(uni_adj, res.vars = "test"), "igraph",
       algorithm = "sugiyama") +
  geom_edge_link(arrow = arrow(length = unit(4, "mm"), angle = 15,
                               type = "closed"),
                 end_cap = circle(8, "mm"),
                 color = "grey20") +
  geom_node_point(aes(color = color), size = 21, show.legend = FALSE) +
  geom_node_text(aes(x = x, y = y, label = name), size = 5) +
  theme_bw() +
  theme(plot.margin = margin(30, 30, 30, 30), panel.border = element_blank(),
        panel.background = element_blank(), axis.title = element_blank(),
        axis.text = element_blank(), axis.line = element_blank(),
        axis.ticks = element_blank(), panel.grid = element_blank()) +
  coord_cartesian(clip = "off")
```

A visualization of the corresponding graph is available from Figure \ref{fig:res-graph}, which highlights the resolving variable `test` in red. Apart from color, the graphical model remains unchanged from what is shown in Figure \ref{fig:graph-plot}.

## Semi-Markovian and Topological Ordering Variant

Section \ref{methodology} focuses on the Markovian case, which assumes that all exogenous variables $U_i$ are mutually independent. However, in practice this requirement is often not satisfied. If a mutual dependency structure between variables $U_i$ exists, this is called a Semi-Markovian model. In the university admission example, we could, for example, have $U_{\text{test}} \not\!\perp\!\!\!\perp U_{\text{score}}$, i.e., latent variables corresponding to variables test and final score being correlated. Such dependencies between latent variables can be represented by dashed, bidirected arrows in the causal diagram, as shown in Figures \ref{fig:semi-markov} and \ref{fig:semi-graph}.

```{tikz, semi-markov, fig.cap = "Causal graphical model corresponding to a Semi-Markovian variant of the university admissions example, introduced in Section \\ref{implementation}.  and visualized in its basic form in Figures \\ref{fig:uni-adm} and \\ref{fig:graph-plot}. Here, we allow for the possibility of a mutual dependency between the latent variables corresponding to variables test and final score.", fig.ext = "png", cache = TRUE, echo = FALSE, out.width = "50%"}

<<setup-graph>>

\begin{tikzpicture}

  \pgfsetarrows{latex-latex};

  \begin{scope}
    \node[rv] (a) at (-3,0) {$A$};
    \node[rv] (v1) at (-1,0) {$E$};
    \node[rv] (v2) at (1,0) {$T$};
    \node[rv] (y) at (3,0) {$Y$};
    \draw[->] (a) -- (v1);
    \draw[->] (a) edge[bend left = 30] (v2);
    \draw[->] (v1) -- (v2);
    \draw[->] (v1) edge[bend left = 30] (y);
    \draw[->] (v2) -- (y);
    \path[<->, dashed] (v2) edge[bend right = 20] (y);
  \end{scope}

\end{tikzpicture}
```

There is an important difference in the adaptation procedure for Semi-Markovian case: when inferring the latent quantiles $U_i$ of variable $V_i$, in the Markovian case, only the direct parents $\pa(V_i)$ are needed. In the Semi-Markovian case, due to correlation of latent variables, using only the $\pa(V_i)$ can lead to biased estimates of the $U_i$. Instead, the set of direct parents needs to be extended, as described in more detail by \citet{tian2002general}. A brief sketch of the argument goes as follows: Let the *C-components* be a partition of the set $V$, such that each *C-component* contains a set of variables which are mutually connected by bidirectional edges. Let $C(V_i)$ denote the entire *C-component* of variable $V_i$. We then define the set of extended parents as

$$\Pa(V_i) := (C(V_i) \cup pa(C(V_i))) \cap \an(V_i),$$

where $\an(V_i)$ is the set of ancestors of $V_i$. The adaptation procedure in the Semi-Markovian case in principle remains the same as outlined in Algorithm \ref{algo:fairadapt}, with the difference that the set of direct parents $\pa(V_i)$ is replaced by $\Pa(V_i)$ at each step.

To include the bidirectional confounding edges in the adaptation, we can pass a `matrix` as `cfd.mat` argument to `fairadapt()` such that:

* `cfd.mat` has the same dimension, column and row names as `adj.mat`.
* `cfd.mat` is symmetric.
* As is the case with the adjacency matrix passed as `adj.mat`, an entry `cfd.mat[i, j] == 1` indicates that there is a bidirectional edge between variables `i` and `j`.

The following code performs fair data adaptation of the Semi-Markovian university admission variant with a mutual dependency between the variables representing test and final scores. For this, we create a matrix `uni_cfd` with the same attributes as the adjacency matrix `uni_adj` and set the entries representing the bidirected edge between vertices `test` and `score` to $1$. Finally, we can pass this confounding matrix as `cfd.mat` to `fairadapt()`. A visualization of the resulting causal graph is available from Figure \ref{fig:semi-graph}.

```{r, semi-markov-uni}
uni_cfd <- matrix(0, nrow = nrow(uni_adj), ncol = ncol(uni_adj),
                  dimnames = dimnames(uni_adj))

uni_cfd["test", "score"] <- 1
uni_cfd["score", "test"] <- 1

semi <- fairadapt(score ~ ., train.data = uni_trn,
                  test.data = uni_tst, adj.mat = uni_adj,
                  cfd.mat = uni_cfd, prot.attr = "gender")
```

```{r, semi-graph, echo = FALSE, fig.width = 6, fig.height = 3, fig.cap = "Visualization of the causal graphical model also shown in Figure \\ref{fig:semi-markov}, obtained when passing a confounding matrix indicating a bidirectional edge between vertices \\texttt{test} and \\texttt{score} to \\texttt{fairadapt()}. The resulting Semi-Markovian setting is also handled by \\texttt{fairadapt()}, extending the basic Markovian formulation introduced in Section \\ref{markovian-scm-formulation}.", out.width = "60%"}

set.seed(17)

ggraph(semi$graph, "igraph", algorithm = "fr") +
  geom_edge_arc(data = lty_selector("dashed"),
                arrow = arrow(length = unit(4, "mm"), angle = 15,
                              ends = "both", type = "closed"),
                start_cap = circle(8, "mm"),
                end_cap = circle(8, "mm"),
                strength = 0.25,
                color = "grey20",
                linetype = "dashed") +
  geom_edge_link(data = lty_selector("solid"),
                 arrow = arrow(length = unit(4, "mm"), angle = 15,
                               type = "closed"),
                 end_cap = circle(8, "mm"),
                 color = "grey20",
                 linetype = "solid") +
  geom_node_point(color = "grey80", size = 21) +
  geom_node_text(aes(x = x, y = y, label = name), size = 5) +
  theme_bw() +
  theme(plot.margin = margin(30, 30, 30, 30), panel.border = element_blank(),
        panel.background = element_blank(), axis.title = element_blank(),
        axis.text = element_blank(), axis.line = element_blank(),
        axis.ticks = element_blank(), panel.grid = element_blank()) +
  coord_cartesian(clip = "off")
```

Alternatively, instead of using the extended parent set $\Pa(V_i)$, we could also use the entire set of ancestors $\an(V_i)$. This approach is implemented as well, and available by specifying a topological ordering. This is achieved by passing a `character` vector, containing the correct ordering of the names appearing in `names(train.data)` as `top.ord` argument to `fairadapt()`. The benefit of using this option is that the specific edges of the causal model $\mathcal{G}$ need not be specified. However, in the linear case, specifying the edges of the graph, so that the quantiles are inferred using only the set of parents, will in principle have better performance.

## Questions of Identifiability

So far we did not discuss whether it is always possible to carry out the counterfactual inference described in Section \ref{methodology}. In the causal literature, an intervention is termed *identifiable* if it can be computed uniquely using the data and the assumptions encoded in the graphical model $\mathcal{G}$. An important result by \cite{tian2002general} states that an intervention do$(X = x)$ on a singleton variable $X$ is identifiable if there is no bidirected path between $X$ and $\ch(X)$. Therefore, our intervention of interest is identifiable if one of the two following conditions are met:

* The model is Markovian.
* The model is Semi-Markovian and, 
  (i) there is no bidirected path between $A$ and $\ch(A)$ and,
  (i) there is no bidirected path between $R_i$ and $\ch(R_i)$ for any resolving variable $R_i$.

Based on this, the `fairadapt()` function may return an error, if the specified intervention is not possible to compute. An additional limitation is that \pkg{fairadapt} currently does not support *front-door identification* \citep[Chapter~3]{pearl2009causality}, meaning that certain special cases which are in principle identifiable are not currently handled. We hope to include this case in a future version.

```{r session-info, include = FALSE}
sessionInfo()
```

