wd---
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)
```

```{r utility, include=FALSE}

plotbar <- function(gg, plot.title=NA,x.lab=""){ 
  res <- gg + geom_bar(fill="#5aa59e")+
  stat_count(aes(group=1,label=scales::percent(..prop..,accuracy = 0.1),
                 hjust=ifelse(..prop..<0.1,-0.2,1.1),
                 color=ifelse(..prop..<0.1,"black","white")
                 ),
             geom="text",size=3)+
  scale_color_identity(guide="none")+
  coord_flip()+xlab(x.lab)+ylab("Count")+
  theme_minimal()+theme(panel.grid.major.y = element_blank(),
                        axis.title.y = element_text(angle=0))
  if( !is.na(plot.title)){
    res <- res + ggtitle(plot.title)
  }
  res
}

to_proportion <- function(f,na.rm=TRUE){
  if(is.data.frame(f)){
    f <- f$Frequency
  }
  if(na.rm){
    f <- na.omit(f)
  }
  if(is.factor(f)){
    f <- droplevels(f)
    f <- table(f)/length(f)
  }else{
      if(sum(f)>1){
        f <- f / sum(f)
      }
  }
  return(f)
}

to_proportion_na <- function(f,na.rm=TRUE){
  if(is.data.frame(f)){
    f <- f$Frequency
  }
  if(is.factor(f)){
    f <- table(f)/length(f)
  }else{
      if(sum(f)>1){
        f <- f / sum(f)
      }
  }
  return(f)
}

index_table <- function(.data,vars=colnames(.data),digits=4){
  res = NULL
  for( v in vars){
    values <- .data[[v]]
    if(is.numeric(values)){
      values <- table(values)
    }
    ix <- data.frame(
            as.list(
              round(c( "Gini"=gini(values), "Shannon"=shan(values), "Simpson"=simp(values),
                       "IIR"=iir(values), 
                       
                       "Renyi_2"=renyi(values, 2), "Renyi_4"=renyi(values, 4), "Renyi_6"=renyi(values, 6),                                                            "Renyi_8"=renyi(values, 8), "Renyi_10"=renyi(values, 10),   "Renyi_12"=renyi(values, 12),
                       "Renyi_14"=renyi(values, 14),   "Renyi_16"=renyi(values,16), "Renyi_18"=renyi(values, 18),  
                       "Renyi_20"=renyi(values, 20), "Renyi_22"=renyi(values, 22),   "Renyi_24"=renyi(values,24),
                       "Renyi_50"=renyi(values, 50), "Renyi_100"=renyi(values, 100), "Renyi_200"=renyi(values, 200), 
                       "Renyi_500"=renyi(values, 500),
                       
                       
                       "Hill_3"=hill(values, 3), "Hill_5"=hill(values, 5), "Hill_7"=hill(values, 7),
                       "Hill_9"=hill(values, 9), "Hill_11"=hill(values, 11), "Hill_13"=hill(values, 13),
                       "Hill_15"=hill(values, 15), "Hill_17"=hill(values, 17), "Hill_19"=hill(values, 19),
                       "Hill_21"=hill(values, 21), "Hill_23"=hill(values, 23), "Hill_25"=hill(values, 25),
                       "Hill_50"=hill(values, 50), "Hill_100"=hill(values, 100), "Hill_200"=hill(values, 200),
                       "Hill_500"=hill(values, 500)),
    
                    digits) #insert index here
          ))
    rownames(ix) <- names(vars)[which(vars==v)]
    res <- rbind(res,ix) 
  }
  res
}

index_table_continuous <- function(.data,vars=colnames(.data),digits=4){
  res = NULL
  for( v in vars){
    values <- .data[[v]]
    if(is.numeric(values)){
      values <- table(values)
    }
    ix <- data.frame(
            as.list(
              round(c( "Theil_T"=theil_t(values), "Theil_L"=theil_l(values), "Gini_cont"=1 - reldist::gini(values),
                       "Pietra"=pietra(values)), digits)  #insert index here
          ))
    rownames(ix) <- names(vars)[which(vars==v)]
    res <- rbind(res,ix) 
  }
  res
}
```


# 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_train <- read.csv(paste(source,"/default_uci_train.csv", sep=''), stringsAsFactors=TRUE)
default_uci_train <- default_uci_train %>% mutate(age_cat=cut(AGE, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))
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_train <- read.csv(paste(source,"/german_credit_train.csv", sep=''),stringsAsFactors=TRUE)
german_credit_train <- german_credit_train %>% mutate(age_cat=cut(Age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))
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_train <- read.csv(paste(source,"/census_income_train.csv", sep=''),stringsAsFactors=TRUE)
census_income_train <- census_income_train %>% mutate(age_cat=cut(age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))
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_train <- read.csv(paste(source,"/student_mat_train.csv", sep=''),stringsAsFactors=TRUE)
student_mat_test <- read.csv(paste(source,"/student_mat_test.csv", sep=''),stringsAsFactors=TRUE)

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

### 8. Term deposit
term_deposit_train <- read.csv(paste(source,"/term_deposit_train.csv", sep=''),stringsAsFactors=TRUE)
term_deposit_train <- term_deposit_train %>% mutate(age_cat=cut(age, breaks=c(-Inf, 25, 45, Inf), labels=c("less than 25","25-45","greater than 45")))
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_train <- read.csv(paste(source,"/titanic_train.csv", sep=''),stringsAsFactors=TRUE)
titanic_test <- read.csv(paste(source,"/titanic_test.csv", sep=''),stringsAsFactors=TRUE)

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

```

# Balance Measures (categorical attributes)

We introduce 6 balance measures in order to detect imbalance in datasets:

- Gini index
- Shannon index
- Rényi entropy (generalization of Shannon)
- Simpson index
- Hill (generalization of Simpson)
- Imbalance ratio

The closer the Imbalance Measure to 0, the more data is affected by $imbalance$ and the higher is the risk that such imbalance would cause effects in the output of an ADM system.

The closer the Imbalance Measure to 1, the higher is the heterogeneity/diversity of data (i.e., categories have similar frequencies).

### Normalized Gini-Simpson index

$$ I = \frac{m}{m-1} \cdot \left( 1-\sum_{i=1}^{m} p_i^2 \right) $$

```{r normalized Gini}
gini <- function(f){
  f <- to_proportion(f)
  
  m = length(f)
  (1 - sum(f^2))*(m/(m-1))
}
```

### Normalized Shannon index (renyi = 1)

$$ H = - \frac{1}{ \ln(m) } \cdot \sum_{j=1}^{m} p_j \cdot ln( p_j ) $$

```{r normalized shannon}
shan <- function(f){
  g <- to_proportion_na(f)
  m = length(g)
  
  f <- to_proportion(f)
  
  -sum(f*log(f))/log(m)
}  
```

### Rényi entropy 

$$ H_{\alpha} = \frac{1}{\ln(m)} \cdot \frac{1}{1 - \alpha} \cdot \ln(\sum_{i=1}^{m} p_i^{\alpha}) $$

```{r normalized Renyi}
renyi <- function(f, alpha){
  g <- to_proportion_na(f)
  m = length(g)
  
  f <- to_proportion(f)
  g <- g**alpha
  
  (log(sum(g)))/((1-alpha)*log(m))
}  
```

### Normalized Simpson index (hill = 2)

$$ D =  \frac{1}{m-1} \cdot \left( \frac{1}{ \sum_{j=1}^{m} p_j^2 } - 1\right) $$

```{r normalized simpson}
simp <- function(f){
  f <- to_proportion_na(f)
  
  m = length(f)
  ((1/sum(f^2))-1)/(m-1)
} 
```

### True diversity index (Hill numbers of order q)

$$ D_{q} = \frac{1}{m-1} \cdot [((\sum_{i=1}^{m} p_i^{q})^{\frac{1}{1-q}})-1] $$

```{r normalized Hill}
hill <- function(f, q){
  g <- to_proportion_na(f)
  m = length(g)
  exp = (1)/(1-q)
  
  g <- g**q
  g <- sum(g)
  g <- g**exp
  
  (g-1)/(m-1)
}  
```
### Inverse Imbalance Ratio

$$ IIR =  \frac{\min{p_j}}{\max{p_j}} $$

```{r inverse imbalance ratio}
iir <- function(f){
  f <- to_proportion_na(f)
  
  min(f)/max(f)
} 
```

### CONTINUOUS IMBALANCE INDICES (for the age attribute)

We introduce 4 balance measure for continuous attributes:

- Theil T
- Theil L
- Gini coefficient 
- Pietra-Ricci

### Theil T

$$ Th_T = 1 - \frac{1}{N*\log(N)} \cdot  \sum_{i=1}^{N}(\frac{x_i}{\overline{x}}\ln(\frac{x_i}{\overline{x}})) $$

```{r Theil T}
theil_t <- function(values){
  N = length(values)
  u = mean(values)
  ratio = values/u
  product = ratio*log(ratio)
  
  1 - sum(product)/(N*log(N))
} 
```

### Theil L

$$ Th_L =  1 - [- \frac{1}{N*\log(N)} \cdot  \sum_{i=1}^{N}\ln(\frac{x_i}{\overline{x}})] $$

```{r Theil L}
theil_l <- function(values){
  N = length(values)
  u = mean(values)
  ratio = values/u
  ratio <- log(ratio)
  
  1 + sum(ratio)/(N*log(N))
} 
```

### Gini coefficient

$$ G = 1 - \frac{\sum_{i=1}^{N}\sum_{j=1}^{N}|x_i - x_j|}{2N^{2}\overline{x}} $$

```{r Gini}
# 1 - reldist::gini(x)
```

### Pietra-Ricci

$$ P = 1 - \frac{\sum_{i=1}^{N}|x_i - \overline{x}|}{2N\overline{x}} $$

```{r Pietra}
pietra <- function(values){
  N = length(values)
  u = mean(values)
  diff <- values - u
  diff <- abs(diff)
  
  1 - sum(diff)/(2*N*u)
} 
```


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

```

# Imbalance detection (categorical)

```{r Imbalance detection, echo = FALSE}

clname <- c("Dataset", "Gini","Shannon","Simpson", "IIR", 
            
            "Renyi_2", "Renyi_4", "Renyi_6",
            "Renyi_8", "Renyi_10", "Renyi_12",
            "Renyi_14", "Renyi_16", "Renyi_18",
            "Renyi_20", "Renyi_22", "Renyi_24",
            "Renyi_50", "Renyi_100","Renyi_200",
            "Renyi_500",
            
            
            "Hill_3", "Hill_5", "Hill_7",
            "Hill_9", "Hill_11", "Hill_13",
            "Hill_15", "Hill_17", "Hill_19",
            "Hill_21", "Hill_23", "Hill_25",
            "Hill_50", "Hill_100", "Hill_200",
            "Hill_500"
            
            )

# 1
indexes_compas <- cbind("COMPAS",index_table(compas,vars.compas))
colnames(indexes_compas) <- clname
kable(indexes_compas)

# 2
indexes_juvenile <- cbind("Juvenile",index_table(juvenile,vars.juvenile))
colnames(indexes_juvenile) <- clname
kable(indexes_juvenile)

# 3
indexes_dccc <- cbind("Default_Uci",index_table(default_uci_train,vars.dccc))
colnames(indexes_dccc) <- clname
kable(indexes_dccc)

# 4
indexes_statlog <- cbind("Statlog",index_table(german_credit_train,vars.statlog))
colnames(indexes_statlog) <- clname
kable(indexes_statlog)

# 5
indexes_income <- cbind("Income",index_table(census_income_train,vars.income))
colnames(indexes_income) <- clname
kable(indexes_income)

# 6
indexes_student_mat <- cbind("Student_Math",index_table(student_mat_train,vars.student))
colnames(indexes_student_mat) <- clname
kable(indexes_student_mat)

#7
indexes_student_por <- cbind("Student_Port",index_table(student_por_train,vars.student))
colnames(indexes_student_por) <- clname
kable(indexes_student_por)

#8
indexes_term_depost <- cbind("Term_deposit",index_table(term_deposit_train,vars.deposit))
colnames(indexes_term_depost) <- clname
kable(indexes_term_depost)

#9
indexes_titanic <- cbind("Titanic",index_table(titanic_train,vars.titanic))
colnames(indexes_titanic) <- clname
kable(indexes_titanic)

#10
indexes_crimes <- cbind("Community_Crime",index_table(community_crime_train,vars.crimes))
colnames(indexes_crimes) <- clname
kable(indexes_crimes)
```


# Sensitive Attributes (continuos)

```{r sensitive attributes}

# PLOT AGE DISTRIBUTION

# 1. Compas
age.compas = c("Age"="age")

# 3. Credid Card Client Default dataset
age.dccc = c("Age"="AGE")

# 4. Statlog (German Credit - SCHUFA) dataset
age.statlog = c("Age"="Age")  

# 5. Census Income dataset
age.income = c("Age"="age") 

# 8. Term deposit dataset
age.deposit = c("Age"="age")

# 9. Titanic dataset
age.titanic = c("Age"="Age")

```

# Imbalance detection (continuous)

```{r Imbalance detection, echo = FALSE}

clname <- c("Dataset", "Theil_T","Theil_L","Gini_cont", "Pietra")

# 1
indexes_compas <- cbind("COMPAS",index_table_continuous(compas,age.compas))
colnames(indexes_compas) <- clname
kable(indexes_compas)

# 3
indexes_dccc <- cbind("Default_Uci",index_table_continuous(default_uci_train,age.dccc))
colnames(indexes_dccc) <- clname
kable(indexes_dccc)

# 4
indexes_statlog <- cbind("Statlog",index_table_continuous(german_credit_train,age.statlog))
colnames(indexes_statlog) <- clname
kable(indexes_statlog)

# 5
indexes_income <- cbind("Income",index_table_continuous(census_income_train,age.income))
colnames(indexes_income) <- clname
kable(indexes_income)


#8
indexes_term_depost <- cbind("Term_deposit",index_table_continuous(term_deposit_train,age.deposit))
colnames(indexes_term_depost) <- clname
kable(indexes_term_depost)

#9
indexes_titanic <- cbind("Titanic",index_table_continuous(titanic_train,age.titanic))
colnames(indexes_titanic) <- clname
kable(indexes_titanic)
```