Definitions of fairness

The first step towards understanding algorithmic fairness is about providing a formal definition of what fairness (or discrimination) means. In light of this, existing intuitive notions have been mathematically formalized, thereby also providing fairness metrics that can be used to quantify discrimination. However, various different notions of fairness exist, and these are sometimes mutually incompatible , meaning they cannot be satisfied at the same time for a given predictor Ŷ. In fact, there is currently no consensus on which notion of fairness is the correct one. Among the many proposal discussed in the literature, the most commonly considered ones include:

1. Demographic parity , which requires the protected attribute A (gender, race, religion etc.) to be independent of a constructed classifier or regressor Ŷ, written as Ŷ⊥   ⊥A.
2. Equality of odds , which requires equal false positive and false negative rates of classifier Ŷ between different groups (females and males for example) written as Ŷ⊥   ⊥A ∣ Y.
3. Calibration , which requires the protected attribute to be independent of the actual outcome given the prediction Y⊥   ⊥A ∣ Ŷ. Intuitively, this means that the protected attribute A does not offer additional information about the outcome Y once we know what the prediction Ŷ is.

Fairness tasks

Apart from choosing a notion of fairness most appropriate to the setting that is analyzed, it is instructive to distinguish between two somewhat different tasks in fairness analysis. The first, usually simpler task is that of bias quantification, or bias measurement. In this case, we are interested in computing metrics from our dataset, in order to verify whether a definition is satisfied or not. Concretely, suppose we have a dataset in which Ŷ is the predicted salary of an employee and A is sex. If we are interested in demographic parity, Ŷ⊥   ⊥A, we can compute the average difference in salary between the male/female groups, that is $$\mathbbm{E}[\widehat{Y} \mid A = \text{male}] - \mathbbm{E}[\widehat{Y} \mid A = \text{female}].$$ Based on this quantity (which is known as total variation, or TV for short), we might decide that either there is a disparity between sexes, or perhaps not. For other definitions of fairness, different metrics would be appropriate.

Sometimes performing the first step of bias measurement is not the end goal. Suppose that we found a large salary gap between sexes in the example above, raising issues about possible discrimination. We then might be interested in correcting this bias, by computing new, more fair predictions (in which, say, the salary gap is lower). This second task is that of bias removal, or bias mitigation, in which we want to remove the undesired bias from our predictions.

Different software tools can be used to perform the two fairness tasks described above. In Figure we provide a graphical overview of some of the available software in and , which are the languages most commonly used for fairness analysis. For a given task and outcome/definition of fairness, we show the existing software packages that the reader might want to use. Since calibration is often satisfied by fitting an unconstrained model, our focus is on demographic parity and equality of odds.

Software options for fair data analysis in and . In the left column, nodes correspond to -packages and in the right column to modules. DP stands for demographic parity, EO equality of odds.

The software landscape for fair ML in is more mature. Currently three well developed packages exist, capable of computing various fairness metrics. Also, these tools support training predictors that satisfy either equality of odds or demographic parity. These repositories are maintained by IBM, maintained by Microsoft and . A further package, , is narrower in scope but suitable for computing fairness metrics on very large datasets.

For , i.e., distributed via CRAN, there are fewer packages that relate to fair machine learning (see Figure ). Available packages include , which implements the non-convex method of , as well as and , which serve as diagnostic tools for measuring algorithmic bias and provide several pre- and post-processing methods for bias mitigation. The package described in this manuscript is a bias removal method which is able to interpolate between demographic parity and calibration notions, and is applicable to both regression and classification settings. In particular, is the only software in Figure which is causally aware. This means that the bias removal performed by can be explained by and related to the causal mechanisms that generated the discrimination in the first place.

A causal approach

The discussion on algorithmic fairness is, however, not restricted to the machine learning domain. There are many legal and philosophical aspects that are paramount. For example, the legal distinction between the disparate impact and disparate treatment doctrines is important for assessing fairness from a legal standpoint. 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 gained attention.

A possible approach to fairness is the use of counterfactual reasoning , 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 , which states that the decision made, should remain fixed, even if, hypothetically, the protected attribute such as race or gender were to be changed (this can be written succinctly as Ŷ_i(a) = Ŷ_i(a′) in the potential outcomes notation). Causal inference can also be used for decomposing the total variation measure into its direct, mediated, and confounded contributions , yielding further insights into demographic parity as a criterion. Furthermore, by introducing the notion of so-called resolving variables, described relaxations of demographic parity, which can possibly be a prohibitively strong notion.

The following sections describe an implementation of the fair data adaptation method outlined in , which combines the notions of counterfactual fairness and resolving variables, and explicitly computes counterfactual instances for individuals. The implementation is available as the -package from CRAN.

Novelty in the package

A first version of was published with the original manuscript . The software has since been developed further and novelty in the package presented in this manuscript includes the following:

• The methodology has been extended from the Markovian to the Semi-Markovian case (allowing noise variables to be correlated), which generalizes the scope of applications.
• Backdoor paths into the protected attribute A are now allowed, meaning that the attribute A does not need to be a root node of the causal graph.
• The user is provided with functionality for uncertainty quantification of the estimates.
• Introduction of S3 classes fairadapt and fairadaptBoot, alongside associated methods, provides a more formalized implementation.
• More flexibility is allowed in the quantile learning step, including different algorithms for quantile regression (linear, forest based, neural networks). Additionally, there is also a possibility of specifying a custom quantile learning function (utilizing S3 dispatch).
• The user is provided with functionality for assessing the quality of the quantile regression fit, which allows for tuning the hyperparameters of the more flexible algorithms.

The rest of the manuscript is organized as follows. In Section we describe the methodology behind , together with reviewing some important concepts of causal inference. In Section we discuss implementation details and provide some general user guidance, followed by Section in which we discuss how to perform uncertainty quantification with . Section illustrates the use of through a large, real-world dataset and a hypothetical fairness application. Finally, in Section we elaborate on some extensions, such as Semi-Markovian models and resolving variables.

University admission example

Consider for example a dataset about students applying for university admission. Let variable A, gender, be the protected attribute (A = a corresponding to females and A = a′ to males). Let E denote 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. A visual representation of how the variables affect each other is given in Figure .

A visual representation of the university admission example. A denotes gender, E previous educational achievement, and T an admissions test score. The outcome Y represents the final score used for the admission decision. The protected attribute A has a discriminatory causal effect on variables E, T, and Y, which we wish to remove.

Attribute A, gender, has a causal effect on variables E, T, and Y. We consider this effect discriminatory and wish to eliminate it, in the following way. For each individual with observed values (a, e, t, y) we want to find a mapping

(a, e, t, y) → (a^(fp), e^(fp), t^(fp), y^(fp)),

which represents the value the person would have obtained in an alternative world where everyone was female. To construct such a mapping, we adapt (transform) the variables in the dataset in order. Explicitly, to a male person with education value e, we assign the transformed value e^(fp), chosen such that

$$\mathbbm{P}(E \geq e \mid A = a') = \mathbbm{P}(E \geq {e}^{(fp)} \mid A = a).$$

The key idea is that the relative educational achievement within the subgroup remains preserved if the protected attribute gender is changed. If, for example, a male person 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¹. After computing transformed educational achievement values corresponding to the female world (E^(fp)), the transformed test score values T^(fp) can be calculated in a similar fashion, but conditional on educational achievement. That is, a male with values (E, T) = (e, t) is assigned a test score t^(fp) such that

$$\mathbbm{P}(T \geq t \mid E = e, A = a') = \mathbbm{P}(T \geq {t}^{(fp)} \mid E = {e}^{(fp)}, A = a),$$

where the value e^(fp) was obtained in the previous step. This step can be visualized as shown in Figure .

In the final step, the outcome variable Y needs to be adjusted. The adaptation is based on the same principle as above, using transformed values of both education and the test score. The transformed value y^(fp) of Y = y is chosen to satisfy

$$\mathbbm{P}(Y \geq y \mid E = e, T = t, A = a') = \mathbbm{P}(Y \geq {y}^{(fp)} \mid E = {e}^{(fp)}, T = {t}^{(fp)}, A = a).$$

The form of counterfactual correction described above is known as recursive substitution . We formalize this approach 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 .

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}, \mathbbm{P}(u) \rangle$, where

• V = {V₁, …, V_n} is the set of observed (endogenous) variables.
• U = {U₁, …, U_n} are latent (exogenous) variables.
• ℱ = {f₁, …, f_n} is the set of functions determining V, v_i ← f_i(pa(v_i), u_i), where pa(V_i) ⊂ V, U_i ⊂ U are the functional arguments of f_i and pa(V_i) denotes the parent vertices of V_i.
• $\mathbbm{P}(u)$ is a distribution over the exogenous variables U.

Any particular SCM is accompanied by a graphical model 𝒢 (a directed acyclic graph, DAG). The set of variables that are inputs of the mechanism f_Vi are the parents of V_i denoted by pa(V_i). Therefore, the graph encodes how variables affect one another. Furthermore, we also write ch(V_i), de(V_i), an(V_i) for the children, descendants, and ancestors of V_i in the graph 𝒢. We assume throughout, without loss of generality, that

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

Equipped with the notion of an SCM and the assumptions (i)-(ii), we can describe the adaptation procedure in the Markovian case, in which all exogenous variables U_i are mutually independent (for the Semi-Markovian case, where variables U_i are allowed to share information, see Section ).

Markovian SCM formulation

Let Y be the outcome of interest, taking values in $\mathbbm{R}$. Let A be the binary protected attribute taking two values a, a′. Denote by X the remaining covariates, and let V = (A, X, Y) denote the observed variables. Our goal is to describe a pre-processing method which transforms the observed variables V into their fair version V^(fp). This is 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 ∣ pa(V_i), A = a and V_i ∣ pa(V_i), A = a′, meaning that the distribution of V_i conditional on pa(V_i) should be indistinguishable between female and male groups (and this should hold 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 since variables U_i are uniformly distributed on [0, 1], the U_i values can be interpreted as the latent quantiles associated with V_i. These latent quantiles are assumed to be preserved when performing the adaptation procedure.

The fair data adaption algorithms starts by fixing A = a for all individuals. After this, the algorithm iterates over descendants of the protected attribute A, in any valid topological order (this topological order is inferred from the causal graph 𝒢, which is also an input of the algorithm). For each V_i, the assignment function f_i and the corresponding quantiles U_i are inferred. Finally, transformed values V_i^(fp) are obtained by evaluating f_i, using quantiles U_i and the transformed parents pa(V_i)^(fp) (see Algorithm ).

The assignment functions f_i of the SCM are always assumed to be unknown, but are inferred non-parametrically at each step. Algorithm 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, the situation is more subtle. A detailed discussion can be found in .

Specifying the graphical model

The algorithm used for fair data adaption in fairadapt() is based on graphical causal model 𝒢 (see Algorithm ). To specify the causal graph, we pass the corresponding adjacency matrix as the adj.mat argument. The convenience function graphModel() turns a graph specified as an adjacency matrix into an annotated graph using the package . Alternatively, by calling the S3 generic visualizeGraph() on a fairadapt object, the user can also inspect the graphical model that was used for the data adaptation.

uni_graph <- graphModel(uni_adj)

The underlying graphical model corresponding to the university admission example (also shown in Figure ).

A visualization of the igraph object returned by graphModel() is shown in Figure . The graph is the same as that in Figure . However, specifying the causal graph is not the only option to perform the data adaptation. A possible alternative is to specify a valid topological ordering over the observable variables V, and specify it as a character vector, using the top.ord argument.

Quantile learning step

The most common and recommended usage of fairadapt() follows a typical machine learning framework. We start by calling the fairadapt() function that performs the quantile learning step and counterfactual correction, based on the train.data argument. Following this, we can call the predict() function on the returned fairadapt S3 object, in order to perform data adaption on new test data. In such a workflow, the adaptation of the training and testing data is done separately. In specific situations, it might be desirable to input both train.data and test.data arguments directly to fairadapt(), which then transforms both the training and testing data jointly. This one-step procedure might be considered when the proportion of test samples compared to train samples is large, and when the train.data has a relatively small sample size. The benefit of this approach is that, even though the outcome Y is not available, other attributes X of test.data can be used for quantile learning step.

The 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. The protected attribute A, passed as scalar-valued character vector prot.attr, should also appear in the column names of train.data and test.data. The test.data argument defaults to NULL, with the intention that test.data is specified as an input to the predict() function at a later stage.

set.seed(22)
fit_qual <- fairadapt(score ~ ., train.data = uni_trn,
                      adj.mat = uni_adj, prot.attr = "gender",
                      eval.qfit = 3L)

quantFit(fit_qual)

##       edu      test     score 
## 0.3446551 0.2796783 0.3428743

Fair-twin inspection

We now turn to a useful property of , which allows the user to explore counterfactual instances for different individuals in the dataset. The university admission example presented in Section 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 (a^(fp), e^(fp), t^(fp), y^(fp)), 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:

ft_basic <- fairTwins(basic, train.id = seq_len(n_samp))
head(ft_basic, n = 3)

##   gender      score score_adapted        edu edu_adapted       test
## 1      1  1.9501728      0.427458  1.3499572   0.8167695  1.6177397
## 2      0 -2.3502495     -2.350250 -1.9779234  -1.9779234 -3.1217962
## 3      1  0.6285619      1.098413  0.6263626   0.1834470  0.5300347
##   test_adapted
## 1    0.8557834
## 2   -3.1217962
## 3    0.1141050

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. Having access to explicit counterfactual instances as above may help justify fair decisions in practice or help guide the choice of the assumed causal model and resolving variables (see Section for resolving variables).

Analyzing the uncertainty

In order to analyze the different sets of predictions Ŷ_test^fair, two slightly different perspectives can be taken, and we elaborate on both in the following sections.

Decision-maker analysis

The first way to analyze the uncertainty of the predictions is from the viewpoint of the decision-maker. By decision-maker, in this context we refer to the individual performing the analysis and obtaining a set of fair predictions. For a decision-maker, it is important to understand how sensitive the outcome of the classification is to uncertainties induced by both finite sample size and the inherent uncertainty induced by discrete variables. Let p_A be the vector of predicted probabilities on the test.data, with length nrow(test.data). Denote nrow(test.data) with n_test. Let p_B be another vector of predicted probabilities (under a different bootstrap repetition). We can consider the following four metrics of uncertainty:

jac_frm <- function(x, modes = "single") {

  jac <- function(a, b) {
    intersection <- length(intersect(a, b))
    union <- length(a) + length(b) - intersection
    intersection / union
  }

  res <- lapply(
    seq(quantile(x[, 1L], 0.05), quantile(x[, 1L], 0.95), 0.01),
    function(tsh) {

      ret <- replicate(100L, {
        col <- sample(ncol(x), 2L)
        jac(which(x[, col[1L]] > tsh), which(x[, col[2L]] > tsh))
      })

      data.frame(tsh = tsh, y = mean(ret), sd = sd(ret),
                 mode = modes)
    }
  )

  do.call(rbind, res)
}

jac_df <- rbind(jac_frm(pred_fin, "Finite Sample"),
                jac_frm(pred_quant, "Quantiles"))

ord_ind <- function(x, modes = "single") {

  res <- replicate(5000L, {
    row <- sample(nrow(x), 2)
    ord <- mean(x[row[1], ] > x[row[2], ])
    max(ord, 1 - ord)
  })

  data.frame(res = res, mode = modes)
}

ord_df <- rbind(ord_ind(pred_fin, "Finite Sample"),
                ord_ind(pred_quant, "Quantiles"))

inv_frm <- function(x, modes = "single") {

  gt <- function(x) x[1L] > x[2L]

  res <- replicate(100L, {
    col <- sample(ncol(x), 2L)
    prm <- order(x[, col[2L]][order(x[, col[1L]])])
    sum(combn(prm, 2L, gt)) / choose(length(prm), 2L)
  })

  data.frame(res = res, mode = modes)
}

inv_df <- rbind(inv_frm(pred_fin, "Finite Sample"),
                inv_frm(pred_quant, "Quantiles"))

prb_frm <- function(x, modes = "single") {
  qnt <- apply(x, 1L, quantile, probs = c(0.05, 0.95))
  data.frame(width = qnt[2L, ] - qnt[1L, ], mode = modes)
}

prb_df <- rbind(prb_frm(pred_fin, "Finite Sample"),
                prb_frm(pred_quant, "Quantiles"))

The results of these four metrics applied to different bootstrap repetitions on the COMPAS dataset are shown in Figure .

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in
## ggplot2 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()`
##   instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning
## was generated.

## Warning: Removed 2 rows containing missing values or values outside the scale
## range (`geom_point()`).

## Warning: Removed 2 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 137 rows containing non-finite outside the scale range
## (`stat_density()`).

Analyzing uncertainty of predictions in the COMPAS dataset from decision-maker’s point of view. Panel A shows how the Jaccard similarity of two repetitions varies depending on the decision threshold. Panel B shows the cumulative distribution of the random variable that indicates whether two randomly selected individuals preserve order (in terms of predicted probabilities) in bootstrap repetitions. Panel C shows the density of the normalized inversion number of between predicted probabilities in bootstrap repetitions. Panel D shows the cumulative distribution function of the 95% confidence interval (CI) width for the predicted probability of different individuals.

Individual analysis

The other point of view we can take in quantifying uncertainty is that of the individual experiencing fair decisions. In this case, we are not so much interested in how much the overall decision set changes (as was the case above), but rather how much variation there is in the estimate for the specific individual. To this end, one might investigate the spread of the predicted values for a specific individual. The spread of the predictions for the first three individuals can be obtained as follows:

ind_prb <- data.frame(
  prob = as.vector(t(pred_quant[seq_len(3), ])),
  individual = rep(c(1, 2, 3), each = fa_boot_quant$n.boot)
)

This spread is visualized in Figure . While it might be tempting to take the mean of the different individual predictions, to have more stable results, this should not be done, as such an approach does not guarantee the fairness constraint fairadapt() aims to achieve. Hence, in order to achieve the desired fairness criteria overall, we have to accept some level of randomization at the level of individual predictions. It would be interesting for future work to quantify and optimize explicitly the trade-off between the desired fairness criteria and the necessary level of randomization.

Analyzing the spread of individual predictions in the COMPAS dataset, resulting from different bootstrap repetitions.

Adding resolving variables

As we mentioned earlier, in some situations the protected attribute A might affect other variables in a non-discriminatory way. For instance, in the Berkeley admissions dataset 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 . A variable R is called resolving if

1. R ∈ de(A), where de(A) are the descendants of A in the causal graph 𝒢.
2. 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 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 , the string "test" can be passed as res.vars to fairadapt().

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

## 
## Call:
## fairadapt(formula = score ~ ., prot.attr = "gender", adj.mat = uni_adj, 
##     train.data = uni_trn, test.data = uni_tst, res.vars = "test")
## 
## Protected attribute:                 gender
## Protected attribute levels:          0, 1
## Adapted variables:                   edu, score
## Resolving variables:                 test
## 
## Number of training samples:          500
## Number of test samples:              500
## Quantile method:                     rangerQuants
## 
## Total variation (before adaptation): -0.7045
## Total variation (after adaptation):  -0.3027

As can be seen from the respective model summary outputs, the total variation after adaptation, in this case, is larger than in the example from Section , 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.

uni_res <- graphModel(uni_adj, res.vars = "test")

Visualization of the causal graph corresponding to the university admissions example introduced in Section with the variable chosen as a and therefore highlighted in red.

A visualization of the corresponding graph is available from Figure , which highlights the resolving variable test in red, but the underlying graphical model remains the same.

Semi-Markovian and topological ordering variant

In Section we focused 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, we are speaking about Semi-Markovian models. In the university admission example, we might have that $U_{\text{test}} \not\!\perp\!\!\!\perp U_{\text{score}}$. That is, 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 and .

Causal graphical model corresponding to a Semi-Markovian variant of the university admissions example, introduced in Section . Here, we allow for the possibility of a mutual dependency between the latent variables corresponding to variables test and score.

There is an important difference in the adaptation procedure for the 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 . 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 bidirected 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) ∪ pa(C(V_i))) ∩ 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 , with the difference that the set of direct parents pa(V_i) is replaced by Pa(V_i) at each step.

To include the bidirected 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 with the adjacency matrix adj.mat, an entry cfd.mat[i, j] == 1 indicates that there is a bidirected 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 .

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")

Visualization of the causal graphical model also shown in Figure , obtained when passing a confounding matrix indicating a bidirected edge between vertices and to . The resulting Semi-Markovian model can also be handled by , extending the basic Markovian formulation introduced in Section .

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 𝒢 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. The topological variant can be invoked as follows:

set.seed(2022)
top_ord <- fairadapt(score ~ ., train.data = uni_trn, test.data = uni_tst, 
                      top.ord = c("gender", "edu", "test", "score"),
                      prot.attr = "gender")

summary(top_ord)

## 
## Call:
## fairadapt(formula = score ~ ., prot.attr = "gender", train.data = uni_trn, 
##     test.data = uni_tst, top.ord = c("gender", "edu", "test", 
##         "score"))
## 
## Protected attribute:                 gender
## Protected attribute levels:          0, 1
## Adapted variables:                   edu, test, score
## 
## Number of training samples:          500
## Number of test samples:              500
## Quantile method:                     rangerQuants
## 
## Total variation (before adaptation): -0.7045
## Total variation (after adaptation):  0.02508

Note that the topological variant for the university admissions dataset is the same as supplying the adjacency matrix, since the causal graph has no missing edges.

Questions of identifiability

So far we did not discuss whether it is always possible to carry out the counterfactual inference described in Section . 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 𝒢. An important result by states that an intervention do(X = x) on a 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,
  1. there is no bidirected path between A and ch(A) and,
  2. there is no bidirected path between R_i and ch(R_i) for any resolving variable R_i.

Based on this, the fairadapt() function might return an error, if the specified intervention is not possible to compute. An additional limitation is that currently does not support front-door identification , meaning that certain special cases, which are in principle identifiable, are currently not handled.

Future avenues to be explored

We conclude with a brief look at the possible extensions of the package, which we hope to consider in future work:

1. Spurious pathways: the package allows for correcting discrimination along causal pathways from the protected attribute A to outcome Y. However, it is also possible that A and Y are associated through a confounding mechanism, and correcting for such spurious association poses an interesting methodological and practical challenge.

2. General identification: as discussed above, there are certain cases of causal graphs in which our do(A = a, R = R(a′)) intervention is identifiable, but the fairadapt() function currently does not support doing so. One of such examples is front-door identification, mentioned above. In a future version of , we hope to cover all scenarios in which identification is possible.

3. Path-specific effects: when using resolving variables (Section ), the user decides to label these variables as “non-discriminatory”, that is, the algorithm is free to distinguish between groups based on these variables. In full generality, a user might be interested in considering all path-specific effects . Such an approach would offer even more flexibility in modeling, since for every attribute-outcome path A → ... → Y, the user could decide whether it is fair or not.

4. Selection bias: a commonly considered problem in causal inference is that of selection bias , when inclusion of individuals into the dataset depends on the observed variables in the dataset. In fairness applications, the presence of selection bias could invalidate our conclusions about discrimination and make our fair predictions biased. Recovering from selection bias algorithmically would therefore be a desirable feature in the package.

5. Non-binary attribute A: one assumption used currently is that the protected attribute A is binary. When considering possible protected attributes as socioeconomic status or education, this assumption may need to be relaxed. We hope to address this in future work.