---
title: "Reproducibility Package for 'A Data Quality Approach to the Identification of Discrimination Risk in Automated Decision Making Systems' _ RQ2 investigation"
author: "Antonio Vetrò, Marco Torchiano, Mariachiara Mecati"
date: 2021-04-12  
output:
 html_document: default
 output_dir: '../results/'
---

This documents contains the reproducibility scripts for the paper

**A Data Quality Approach to the Identification of Discrimination Risk in Automated Decision Making Systems**

by Mariachiara Mecati, Antonio Vetrò, Marco Torchiano

submitted/accepted for publication at Government Information Quarterly.

The goal of our study is to anticipate possible discrimination risks on the basis of the balance features of the training data employed by ADM systems. To this end we formulated two research questions that will drive our investigation:

RQ 1. _How are existing measures able to detect imbalance among the classes of a given attribute in a dataset?_

RQ 2. _Are existing measures able to reveal a discrimination risk when an ADM system is trained with such data?_

In this script we present our code to start investigating the second Research Question. 
The final results are summarized in unfair-balance.xlsx .

```{r setup, warning=FALSE, error=FALSE, message=FALSE, include=FALSE}
library(ggplot2)
library(tidyverse)
library(knitr)
library(dplyr)
```

# Data loading

We read data from "../Tesi/my code/data/DATASET"

Note that we do NOT exclude category `NA` from the data and we focus only on Categorical attributes.

```{r data, echo = FALSE}

source = 'run5'

### 1. COMPAS Recidivism racial bias dataset
compas <- read.csv(paste(source,"/d_compas.csv", sep=''), stringsAsFactors=TRUE)

### 2. Recidivism in Juvenile Justice dataset
juvenile <- read.csv(paste(source,"/d_juvenile.csv", sep=''),stringsAsFactors=TRUE)

### 3. Credit card default dataset
default_uci_test <- read.csv(paste(source,"/default_uci_test.csv", sep=''), stringsAsFactors=TRUE)
default_uci_test <- default_uci_test %>% mutate(age_cat=cut(AGE, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))

### 4. Statlog (German Credit - SCHUFA) dataset
german_credit_test <- read.csv(paste(source,"/german_credit_test.csv", sep=''),stringsAsFactors=TRUE)
german_credit_test <- german_credit_test %>% mutate(age_cat=cut(Age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))

### 5. Census Income dataset
census_income_test <- read.csv(paste(source,"/census_income_test.csv", sep=''),stringsAsFactors=TRUE)
census_income_test <- census_income_test %>% mutate(age_cat=cut(age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))

### 6. Student Mat
student_mat_test <- read.csv(paste(source,"/student_mat_test.csv", sep=''),stringsAsFactors=TRUE)

### 7. Stundent Por
student_por_test <- read.csv(paste(source,"/student_port_test.csv", sep=''),stringsAsFactors=TRUE)

### 8. Term deposit
term_deposit_test <- read.csv(paste(source,"/term_deposit_test.csv", sep=''),stringsAsFactors=TRUE)
term_deposit_test <- term_deposit_test %>% mutate(age_cat=cut(age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))

### 9. Titanic
titanic_test <- read.csv(paste(source,"/titanic_test.csv", sep=''),stringsAsFactors=TRUE)
titanic_test <- titanic_test %>% mutate(age_cat=cut(Age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))

### 10. Community and Crimes 
community_crime_test <- read.csv(paste(source,"/community_crime_test.csv", sep=''),stringsAsFactors=TRUE)

```


# Fairness criteria

We aim to check whether the ADM system, which assigned a predicted class, behaved fairly w.r.t. the different values of a sensitive attribute.

We assessed the unfairness of automated classifications relying on two criteria formalized in https://fairmlbook.org/ : 
we consider a sensitive categorical attribute $A$ that can assume values $(a, b, ...)$, a target variable $Y$, and a predicted class $R$.

### Independence criterion

The Independence criterion requires the acceptance rate to be the same in all groups, where acceptance correspond to the event R=1. In term of probability, this condition correspond to the following constraint:

P( R = 1 \mid A = a ) = P( R = 1 \mid A = b ) = ... 

if $A$ is binary, we can compute the Independence \emph{unfairness} measure as:

$$ \mathfrak{U}_I(a_1,a_2) = | P( R = 1 \mid A = a_1 ) - P( R = 1 \mid A = a_2 ) | $$

### Separation criterion

The Separation criterion, when R is binary (i.e., R=0 or R=1 and thus Y=0 or Y=1), requires the equivalence of true positive rate and false positive rate for each level of the protected attribute under analysis:

TPR
* P( R = 1 \mid Y = 1, A = a_1 ) = P( R = 1 \mid Y = 1, A = a_2 ) = ...

FPR
* P( R = 1 \mid Y = 0, A = a_1 ) = P( R = 1 \mid Y = 0, A = a_2 ) = ...

If A is binary (that is, $A$ = $a_1$ or $a_2$), the we can compute two Separation \emph{unfairness} measures ($\mathfrak{U}$) as:

equality of opportunity with Y=1 (TPR)
$$ \mathfrak{U}_{S\_TPR}(a_1,a_2) = | P( R = 1 \mid Y = 1, A = a_1 ) - P( R = 1 \mid Y = 1, A = a_2 ) | $$
equality of opportunity with Y=0 (FPR)
$$ \mathfrak{U}_{S\_FPR}(a_1,a_2) = |P( R = 1 \mid Y = 0, A = a_1 ) - P( R = 1 \mid Y = 0, A = a_2 )| $$

### Sufficiency criterion
The suffiency criterion, where R is binary (i.e., R=0 or R=1 and thus Y=0 or Y=1), requires the equality of PPV and NPV for each level of the protected attribute under analysis:

PPV
* P( Y = 1 \mid R = 1, A = a_1 ) = P( Y = 1 \mid R = 1, A = a_2 ) = ...

NPV
* P( Y = 1 \mid R = 0, A = a_1 ) = P( Y = 1 \mid R = 0, A = a_2 ) = ...

if $A$ is binary, we can compute the Independence \emph{unfairness} measure as:

PPV
$$ \mathfrak{U}_{S\_PPV}(a_1,a_2) = | P( Y = 1 \mid R = 1, A = a_1 ) - P( Y = 1 \mid R = 1, A = a_2 ) | $$
NPV
$$ \mathfrak{U}_{S\_NPV}(a_1,a_2) = |P( Y = 1 \mid R = 0, A = a_1 ) - P( Y = 1 \mid R = 0, A = a_2 )| $$

### Overall accuracy equality
A classifier satisfies this definition if the subject in the two groups (binary case) have equal prediction accuracy, that is, the probability of a subject from one class to be assigned to it.
* P( R = Y \mid A = a_1 ) = P( R = Y \mid A = a_2 ) = ...

if $A$ is binary, we can compute the Independence \emph{unfairness} measure as:

Group Accuracy
$$ \mathfrak{U}_{S\_GroupAccuracy}(a_1,a_2) = | P( R = Y \mid A = a_1 ) - P( R = Y \mid A = a_2 ) | $$

### Treatment equality
A classifier satisfies this definition if the subjects in two groups (binary case) have an equal ratio of FN and FP.
* \frac{FN_{A=a1}}{FP_{A=a1}} = \frac{FN_{A=a2}}{FP_{A=a2}} = ...

if $A$ is binary, we can compute the Independence \emph{unfairness} measure as:

Group ratio
$$ \mathfrak{U}_{S\_GroupRatio}(a_1,a_2) = | \frac{FN_{A=a1}}{FP_{A=a1}} - \frac{FN_{A=a2}}{FP_{A=a2}}) | $$

### To summarize

The definitions above can be easily extended to the case of non-binary attributes -- i.e. $m>2$ -- by taking the mean of indexes computed considering all the possible pairs of levels in $A$:

$$ \mathfrak{U}(a_1, ..., a_m) = \frac{2}{m(m-1)}\sum_{i=1}^{m-1} \sum_{j=i+1}^m \mathfrak{U}(a_i,a_j) $$

The unfairness measures defined above range in the interval $[0,1]$. They assume values equal to zero for a perfectly fair classification and higher values for unfair behavior.

```{r utility - Fairness, include=FALSE}

unfairness <- function(.data,dataset_name, vars=colnames(.data),
                     target=colnames(.data),score=colnames(.data),digits=4){
  
  target  <- .data[[target]]
  score  <- .data[[score]]

  res_final <- NULL
  
  for( v in vars){  # v = sensitive attribute
    
    values <- .data[[v]]
    #if(is.numeric(values)){
    #  values <- table(values)
    #}
    
    l_values <- NULL
    l_values <- levels(as.factor(values))
    len_values <- length(l_values)
    #print(len_values)
    #print('LEVELS',str(l_values))
    #print('VALUE',str(values))
    
    res <- NULL
    for(i in 1:(len_values-1)){
      
    for(j in (i+1):(len_values)){
    
    ix <- data.frame(
      
              ####### INDEPENDENCE criterion:
              
              "Ind.a" =  ind.a <-
                ifelse((filter(.data, values==l_values[i]) %>% nrow()) !=0,
                filter(.data, score==1 & values==l_values[i]) %>% nrow() /
                filter(.data, values==l_values[i]) %>% nrow(), 0),

              "Ind.b" = ind.b <-
                ifelse((filter(.data, values==l_values[j]) %>% nrow()) !=0,
                filter(.data, score==1 & values==l_values[j]) %>% nrow() /
                filter(.data, values==l_values[j]) %>% nrow(), 0),
          
              "Diff.Ind" = diff.ind <- 
                round(abs(ind.a-ind.b),digits),

              ####### TRUE POSITIVE RATE (Separation criterion):
                
              "TP.a" =  tp.a <- 
                ifelse((filter(.data, target==1 & values==l_values[i]) %>% nrow()) !=0,
                filter(.data, score==1 &  target==1 &
                   values==l_values[i]) %>% nrow() /
                filter(.data, target==1 & values==l_values[i]) %>% nrow(), 0),

              "TP.b" = tp.b <-
                ifelse((filter(.data, target==1 & values==l_values[j]) %>% nrow()) != 0,
                filter(.data, score==1 &  target==1 &
                   values==l_values[j]) %>% nrow() /
                filter(.data, target==1 & values==l_values[j]) %>% nrow(), 0),

              "Diff.TP" = diff.tp <-
                round(abs(tp.a-tp.b),digits),
              
              ######## FALSE POSITIVE RATE (Separation criterion):
              
              "FP.a" =  fp.a <-
                ifelse((filter(.data, target==0 & values==l_values[i]) %>% nrow()) != 0,
                filter(.data, score==1 &  target==0 &
                   values==l_values[i]) %>% nrow() /
                filter(.data, target==0 & values==l_values[i]) %>% nrow(), 0),
              
              "FP.b" = fp.b <-
                ifelse((filter(.data, target==0 & values==l_values[j]) %>% nrow()) !=0,
                filter(.data, score==1 &  target==0 &
                   values==l_values[j]) %>% nrow() /
                filter(.data, target==0 & values==l_values[j]) %>% nrow(), 0),
              
              "Diff.FP" = diff.fp <- 
                round(abs(fp.a-fp.b),digits),
              
              
              ######### POSITIVE PREDICTIVE VALUES (Sufficiency criteria):
                
              "PPV.a" =  ppv.a <- 
                ifelse((filter(.data, score==1 & values==l_values[i]) %>% nrow()) !=0,
                filter(.data, score==1 &  target==1 &
                   values==l_values[i]) %>% nrow() /
                filter(.data, score==1 & values==l_values[i]) %>% nrow(), 0),

              "PPV.b" = ppv.b <-
                ifelse((filter(.data, score==1 & values==l_values[j]) %>% nrow()) != 0,
                filter(.data, score==1 &  target==1 &
                   values==l_values[j]) %>% nrow() /
                filter(.data, score==1 & values==l_values[j]) %>% nrow(), 0),

              "Diff.PPV" = diff.ppv <-
                round(abs(ppv.a-ppv.b),digits),
              
              
              ######## NEGATIVE PREDICTIVE VALUES (Sufficiency criteria):
                
              "NPV.a" =  npv.a <- 
                ifelse((filter(.data, score==0 & values==l_values[i]) %>% nrow()) !=0,
                filter(.data, score==0 &  target==1 &
                   values==l_values[i]) %>% nrow() /
                filter(.data, score==0 & values==l_values[i]) %>% nrow(), 0),

              "NPV.b" = npv.b <-
                ifelse((filter(.data, score==0 & values==l_values[j]) %>% nrow()) != 0,
                filter(.data, score==0 &  target==1 &
                   values==l_values[j]) %>% nrow() /
                filter(.data, score==0 & values==l_values[j]) %>% nrow(), 0),

              "Diff.NPV" = diff.npv <-
                round(abs(npv.a-npv.b),digits),
              
              
              ######## OVERALL ACCURACY EQUALITY:
              
              "OAE.a" =  oae.a <-
                ifelse((filter(.data, values==l_values[i]) %>% nrow()) !=0,
                filter(.data, score==target & values==l_values[i]) %>% nrow() /
                filter(.data, values==l_values[i]) %>% nrow(), 0),

              "OAE.b" = oae.b <-
                ifelse((filter(.data, values==l_values[j]) %>% nrow()) !=0,
                filter(.data, score==target & values==l_values[j]) %>% nrow() /
                filter(.data, values==l_values[j]) %>% nrow(), 0),
          
              "Diff.OAE" = diff.oae <- 
                round(abs(oae.a-oae.b),digits)
            

              ######## TREATMENT EQUALITY: FN/FP
              #
              #"FN.a" =  fn.a <-
              #  ifelse((filter(.data, target==1 & values==l_values[i]) %>% nrow()) != 0,
              #  filter(.data, score==0 &  target==1 &
              #     values==l_values[i]) %>% nrow(), 0),
#
              #"FP.a" =  fp.a <-
              #  ifelse((filter(.data, target==0 & values==l_values[i]) %>% nrow()) != 0,
              #  filter(.data, score==1 &  target==0 &
              #     values==l_values[i]) %>% nrow(), 0),
#
              #"FN.b" =  fn.b <-
              #  ifelse((filter(.data, target==1 & values==l_values[j]) %>% nrow()) != 0,
              #  filter(.data, score==0 &  target==1 &
              #     values==l_values[j]) %>% nrow(), 0),
#
              #"FP.b" =  fp.b <-
              #  ifelse((filter(.data, target==0 & values==l_values[j]) %>% nrow()) != 0,
              #  filter(.data, score==1 &  target==0 &
              #     values==l_values[j]) %>% nrow(), 0),
              #
              #"Diff.RATIO" = diff.ratio <- 
              #  round(abs(fn.a/fp.a - fn.b/fp.b),digits)
            
                )

    res <- rbind(res,ix)
    
  }
  } 
  
  res_final <- rbind( res_final, 
                      data.frame(dataset_name, names(vars)[which(vars==v)],
                                round(mean(res$Diff.Ind),digits),
                                round(mean(res$Diff.TP),digits),
                                round(mean(res$Diff.FP),digits),
                                round(mean(res$Diff.PPV),digits),
                                round(mean(res$Diff.NPV),digits),
                                round(mean(res$Diff.OAE),digits)
                                #,100*round(mean(res$Diff.RATIO),digits)
                                ))
                                    
  }
  
  colnames(res_final) <- c("Dataset","Attribute", "Diff Ind" ,"Diff TPR","Diff FPR",
                           "Diff PPV","Diff NPV", "Diff OAE")#, "Diff Ratio")
  res_final
}
```
Note that, in order to apply the fairness function, for each dataset we set the target variable equal to 0 or 1 and, analogously, we set the predictions in such a way to obtain a binary classification, with score equal to 0 or 1.

# Sensitive Attributes (categorical)

```{r sensitive attributes}

# 1. COMPAS Recidivism racial bias dataset
vars.compas = c("Ethnicity"="race","Sex"="sex","Age category"="age_cat")

# 2. Recidivism in Juvenile Justice dataset
vars.juvenile = c("Sex"="V1_sexe","Stranger"="V2_estranger","Country of origin"="V3_nacionalitat","Area of origin"="V4_nacionalitat_agrupat","Age category"="V5_edat_fet_agrupat","Age"="V8_edat_fet")

# 3. Credid Card Client Default dataset
vars.dccc = c("Sex"="SEX","Education"="EDUCATION", "Marriage"="MARRIAGE", "Age_cat"="age_cat")

# 4. Statlog (German Credit - SCHUFA) dataset
vars.statlog = c("Status"="Status","Sex"="Sex","Foreign worker"="Foreign_Worker", "Age_cat"="age_cat")  

# 5. Census Income dataset
vars.income = c("Education"="education","Race"="race","Sex"="sex","Native country"="native.country", "Age_cat"="age_cat") 

# 6/7. Student Performance dataset

### All Sensitive Attributes (divided in 2 groups):
#1. Categorical Attributes: 
#   "Sex"="sex","Mother's job"="Mjob","Father's job"="Fjob"
#2. Continuous Attributes converted into Categorical : 
#   "Age"="age_f","Mother's education"="Medu_f","Father's education"="Fedu_f"

vars.student = c("Sex"="sex","Age"="age_f",
        "Mother's job"="Mjob","Father's job"="Fjob",
        "Mother's education"="Medu_f","Father's education"="Fedu_f")

# 8. Term deposit dataset
vars.deposit = c("Education"="education", "Marriage"="marital", "Age_cat"="age_cat")

# 9. Titanic dataset
vars.titanic = c("Sex"="Sex")

# 10. Community and Crime dataset
vars.crimes = c("Race"="Afro.americans","Age"="X12.29.people")

```

# Fairness assessment

```{r Fairness assessment, echo = FALSE}
unfair = NULL

model_predictions = 'rf_pred_SMOTE'   # Choose model among lr, svm, knn, rf (model_pred)

# 1. COMPAS Recidivism racial bias dataset

print("Target Compas:")
print(levels(as.factor(compas$two_year_recid)))

print("Score Compas:")
print(levels((as.factor(compas$riskScoreBin))))

unfair <- rbind( unfair,               
               unfairness(compas,"COMPAS -",vars.compas,
               "two_year_recid","riskScoreBin"))

# 2. Recidivism in Juvenile Justice dataset

print("Target Juvenile:")
print(levels(as.factor(juvenile$recidivism_juv)))

print("Score Juvenile:")
print(levels((as.factor(juvenile$recidivism_juv_score))))

unfair <- rbind( unfair,               
               unfairness(juvenile,"Juvenile -",vars.juvenile, 
               "recidivism_juv","recidivism_juv_score"))

# 3. Credid Card Client Default dataset

print("Target dccc:")
print(levels(as.factor(default_uci_test$default.payment.next.month)))

print("Score dccc:")
print(levels((as.factor(default_uci_test$lr_pred))))

unfair <- rbind( unfair,      
               unfairness(default_uci_test,"dccc -",vars.dccc,
                "default.payment.next.month",model_predictions))

# 4. Statlog (German Credit - SCHUFA) dataset

print("Target Statlog:")
print(levels(as.factor(german_credit_test$costMatrix)))

print("Score Statlog:")
print(levels((as.factor(german_credit_test$lr_pred))))

unfair <- rbind( unfair,          
               unfairness(german_credit_test,"Statlog -",vars.statlog,
                "costMatrix",model_predictions))

# 5. Census Income dataset

print("Target Income:")
print(levels(as.factor(census_income_test$test_income)))

print("Score Income:")
print(levels((as.factor(census_income_test$lr_pred))))

unfair <- rbind( unfair,              
               unfairness(census_income_test,"Income -",vars.income,
               "test_income",model_predictions))

# 6/7. Student Performance dataset

# Students - Mathematics

print("Target Students Math:")
print(levels(as.factor(student_mat_test$G3_target)))

print("Score Students Math:")
print(levels((as.factor(student_mat_test$lr_pred))))

unfair <- rbind( unfair,           
               unfairness(student_mat_test,"Student_Math -", 
               vars.student,"G3_target",model_predictions))

# Students - Portuguese

print("Target Students Port:")
print(levels(as.factor(student_por_test$G3_target)))

print("Score Students Port:")
print(levels((as.factor(student_por_test$lr_pred))))

unfair <- rbind( unfair,              
               unfairness(student_por_test,"Student_Port -",
               vars.student,"G3_target",model_predictions))

# 8. Term deposit dataset

print("Target Term deposit:")
print(levels(as.factor(term_deposit_test$y)))

print("Score Term deposit:")
print(levels((as.factor(term_deposit_test$lr_pred))))

unfair <- rbind( unfair,              
               unfairness(term_deposit_test,"Term deposit -",vars.deposit,
               "y",model_predictions))

# 9. Titanic dataset

print("Target Titanic:")
print(levels(as.factor(titanic_test$Survived)))

print("Score Titanic:")
print(levels((as.factor(titanic_test$lr_pred))))

unfair <- rbind( unfair,              
               unfairness(titanic_test,"Titanic -",vars.titanic,
               "Survived",model_predictions))

## 10. Community and Crime dataset
#
#print("Target Crimes:")
#print(levels(as.factor(community_crime_test$y)))
#
#print("Score Crimes:")
#print(levels((as.factor(community_crime_test$lr_pred))))
#
#unfair <- rbind( unfair,              
#               unfairness(community_crime_test,"Community crime -",vars.crimes,
#               "y",model_predictions))

kable(unfair)
```