stan_modelling_empirical.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
# Author: Riku Laine
# Date: 26JUL2019
# Project name: Potential outcomes in model evaluation
# Description: This script creates the figures and results used 
#              in empirical data experiments.
#
# Parameters:
# -----------
# (1) figure_path : file name for saving the created figures.
# (2) nIter : Number of train test splits to perform on the data.
# (3) group_amount : How many groups if Jung-inspired model is used.
# (4) stan_code_file_name : Name of file containing the stan model code.
# (5) sigma_tau : Values of prior variance for the Jung-inspired model.
# (6) data_path : File of compas data.

"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as scs
import scipy.special as ssp
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
import pystan

plt.switch_backend('agg')

import sys

# figure storage name
figure_path = sys.argv[1]

# Iterations
nIter = int(sys.argv[2])

# How many groups if jung model is used
group_amount = int(sys.argv[3])

# Name of stan model code file
stan_code_file_name = sys.argv[4]

# Variance prior
sigma_tau = float(sys.argv[5])

# figure storage name
data_path = sys.argv[6]

# Prefix for the figures and log files

print("These results have been obtained with the following settings:")

print("Number of train test splits :", nIter)

print("Number of groups:", group_amount)

print("Prior for the variances:", sigma_tau)

save_name = "sl_compas"

def inv_logit(x):
    return 1.0 / (1.0 + np.exp(-1.0 * x))


def logit(x):
    return np.log(x) - np.log(1.0 - x)


def inverse_cumulative(x, mu, sigma):
    '''Compute the inverse of the cumulative distribution of logit-normal
    distribution at x with parameters mu and sigma (mean and st.dev.).'''

    return inv_logit(ssp.erfinv(2 * x - 1) * np.sqrt(2 * sigma**2) - mu)

##########################

# ## Evaluator modules

# ### Convenience functions

def fitPredictiveModel(x_train, y_train, x_test, class_value, model_type=None):
    '''
    Fit a predictive model (default logistic regression) with given training 
    instances and return probabilities for test instances to obtain a given 
    class label.
    
    Arguments:
    ----------
    
    x_train -- x values of training instances
    y_train -- y values of training instances
    x_test -- x values of test instances
    class_value -- class label for which the probabilities are counted for.
    model_type -- type of model to be fitted.
    
    Returns:
    --------
    (1) Trained predictive model
    (2) Probabilities for given test inputs for given class.
    '''

    if model_type is None or model_type in ["logistic_regression", "lr"]:
        # Instantiate the model (using the default parameters)
        logreg = LogisticRegression(solver='lbfgs')

        # Check shape and fit the model.
        if x_train.ndim == 1:
            logreg = logreg.fit(x_train.values.reshape(-1, 1), y_train)
        else:
            logreg = logreg.fit(x_train, y_train)

        label_probs_logreg = getProbabilityForClass(x_test, logreg,
                                                    class_value)

        return logreg, label_probs_logreg

    elif model_type in ["random_forest", "rf"]:
        # Instantiate the model
        forest = RandomForestClassifier(n_estimators=100, max_depth=3)

        # Check shape and fit the model.
        if x_train.ndim == 1:
            forest = forest.fit(x_train.values.reshape(-1, 1), y_train)
        else:
            forest = forest.fit(x_train, y_train)

        label_probs_forest = getProbabilityForClass(x_test, forest,
                                                    class_value)

        return forest, label_probs_forest

    elif model_type == "fully_random":

        label_probs = np.ones_like(x_test) / 2

        model_object = lambda x: 0.5

        return model_object, label_probs
    else:
        raise ValueError("Invalid model_type!", model_type)


def getProbabilityForClass(x, model, class_value):
    '''
    Function (wrapper) for obtaining the probability of a class given x and a 
    predictive model.

    Arguments:
    -----------
    x -- individual features, an array of shape (observations, features)
    model -- a trained sklearn model. Predicts probabilities for given x. 
        Should accept input of shape (observations, features)
    class_value -- the resulting class to predict (usually 0 or 1).

    Returns:
    --------
    (1) The probabilities of given class label for each x.
    '''
    if x.ndim == 1:
        # if x is vector, transform to column matrix.
        f_values = model.predict_proba(np.array(x).reshape(-1, 1))
    else:
        f_values = model.predict_proba(x)

    # Get correct column of predicted class, remove extra dimensions and return.
    return f_values[:, model.classes_ == class_value].flatten()


# ### Contraction algorithm
# 
# Below is an implementation of Lakkaraju's team's algorithm presented in 
# [their paper](https://helka.finna.fi/PrimoRecord/pci.acm3098066). Relevant
# parameters to be passed to the function are presented in the description.

def contraction(df, judgeIDJ_col, decisionT_col, resultY_col, modelProbS_col,
                accRateR_col, r):
    '''
    This is an implementation of the algorithm presented by Lakkaraju
    et al. in their paper "The Selective Labels Problem: Evaluating 
    Algorithmic Predictions in the Presence of Unobservables" (2017).

    Arguments:
    ----------
    df -- The (Pandas) data frame containing the data, judge decisions,
        judge IDs, results and probability scores.
    judgeIDJ_col -- String, the name of the column containing the judges' IDs
        in df.
    decisionT_col -- String, the name of the column containing the judges' decisions
    resultY_col -- String, the name of the column containing the realization
    modelProbS_col -- String, the name of the column containing the probability
        scores from the black-box model B.
    accRateR_col -- String, the name of the column containing the judges' 
        acceptance rates
    r -- Float between 0 and 1, the given acceptance rate.

    Returns:
    --------
    (1) The estimated failure rate at acceptance rate r.
    '''
    # Get ID of the most lenient judge.
    most_lenient_ID_q = df[judgeIDJ_col].loc[df[accRateR_col].idxmax()]

    # Subset. "D_q is the set of all observations judged by q."
    D_q = df[df[judgeIDJ_col] == most_lenient_ID_q].copy()

    # All observations of R_q have observed outcome labels.
    # "R_q is the set of observations in D_q with observed outcome labels."
    R_q = D_q[D_q[decisionT_col] == 1].copy()

    # Sort observations in R_q in descending order of confidence scores S and
    # assign to R_sort_q.
    # "Observations deemed as high risk by B are at the top of this list"
    R_sort_q = R_q.sort_values(by=modelProbS_col, ascending=False)

    number_to_remove = int(
        round((1.0 - r) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])))

    # "R_B is the list of observations assigned to t = 1 by B"
    R_B = R_sort_q[number_to_remove:R_sort_q.shape[0]]

    return np.sum(R_B[resultY_col] == 0) / D_q.shape[0]


# ### Evaluators

def contractionEvaluator(df, featureX_col, judgeIDJ_col, decisionT_col,
                         resultY_col, accRateR_col, r):

    train, test = train_test_split(df, test_size=0.5)

    B_model, predictions = fitPredictiveModel(
        train.loc[train[decisionT_col] == 1, featureX_col],
        train.loc[train[decisionT_col] == 1, resultY_col], test[featureX_col],
        0)

    test = test.assign(B_prob_0_model=predictions)

    # Invoke the original contraction.
    FR = contraction(test,
                     judgeIDJ_col=judgeIDJ_col,
                     decisionT_col=decisionT_col,
                     resultY_col=resultY_col,
                     modelProbS_col="B_prob_0_model",
                     accRateR_col=accRateR_col,
                     r=r)

    return FR


def trueEvaluationEvaluator(df, featureX_col, decisionT_col, resultY_col, r,
                            fit_model=True):

    train, test = train_test_split(df, test_size=0.5)

    if fit_model:
        B_model, predictions = fitPredictiveModel(train[featureX_col],
                                                  train[resultY_col],
                                                  test[featureX_col], 0)
    
        test = test.assign(B_prob_0_model=predictions)

    test.sort_values(by='B_prob_0_model', inplace=True, ascending=True)

    to_release = int(round(test.shape[0] * r))

    return np.sum(test[resultY_col][0:to_release] == 0) / test.shape[0]


def labeledOutcomesEvaluator(df,
                             featureX_col,
                             decisionT_col,
                             resultY_col,
                             r,
                             adjusted=False,
                             fit_model=True):

    train, test = train_test_split(df, test_size=0.5)

    if fit_model:
        B_model, predictions = fitPredictiveModel(
            train.loc[train[decisionT_col] == 1, featureX_col],
            train.loc[train[decisionT_col] == 1, resultY_col], test[featureX_col],
            0)

        test = test.assign(B_prob_0_model=predictions)

    test_observed = test.loc[test[decisionT_col] == 1, :]

    test_observed = test_observed.sort_values(by='B_prob_0_model',
                                              inplace=False,
                                              ascending=True)

    to_release = int(round(test_observed.shape[0] * r))

    if adjusted:
        return np.mean(test_observed[resultY_col][0:to_release] == 0)

    return np.sum(
        test_observed[resultY_col][0:to_release] == 0) / test.shape[0]


def humanEvaluationEvaluator(df, judgeIDJ_col, decisionT_col, resultY_col,
                             accRateR_col, r):

    # Get judges with correct leniency as list
    is_correct_leniency = df[accRateR_col].round(1) == r

    # No judges with correct leniency
    if np.sum(is_correct_leniency) == 0:
        return np.nan

    correct_leniency_list = df.loc[is_correct_leniency, judgeIDJ_col]

    # Released are the people they judged and released, T = 1
    released = df[df[judgeIDJ_col].isin(correct_leniency_list)
                  & (df[decisionT_col] == 1)]

    # Get their failure rate, aka ratio of reoffenders to number of people judged in total
    return np.sum(released[resultY_col] == 0) / correct_leniency_list.shape[0]


###################
    
# Read in the data
compas_raw = pd.read_csv(data_path)

# Select columns
compas = compas_raw[[
    'age', 'c_charge_degree', 'race', 'age_cat', 'score_text', 'sex',
    'priors_count', 'days_b_screening_arrest', 'decile_score', 'is_recid',
    'two_year_recid', 'c_jail_in', 'c_jail_out'
]]

# Subset values, see reasons in ProPublica methodology.
compas = compas.query('days_b_screening_arrest <= 30 and \
                      days_b_screening_arrest >= -30 and \
                      is_recid != -1 and \
                      c_charge_degree != "O"')

# Drop row if score_text is na
compas = compas[compas.score_text.notnull()]

# Recode recidivism values to fit earlier notation
# So here result_Y = 1 - is_recid (inverted binary coding).
compas['result_Y'] = np.where(compas['is_recid'] == 1, 0, 1)

# Convert string values to dummies, drop first so full rank
compas_dummy = pd.get_dummies(
    compas,
    columns=['c_charge_degree', 'race', 'age_cat', 'sex'],
    drop_first=True)

compas_dummy.drop(columns=[
    'age', 'days_b_screening_arrest', 'c_jail_in', 'c_jail_out',
    'two_year_recid', 'score_text', 'is_recid'
],
    inplace=True)

# Shuffle rows for random judge assignment
compas_shuffled = compas_dummy.sample(frac=1)

nJudges_M = 9

# Assign judges as evenly as possible
judge_ID = pd.qcut(np.arange(len(compas_shuffled)), nJudges_M, labels=False)

# Assign fixed leniencies from 0.1 to 0.9
judge_leniency = np.arange(1, 10) / 10

judge_leniency = judge_leniency[judge_ID]

compas_shuffled = compas_shuffled.assign(judge_ID=judge_ID,
                                         judge_leniency=judge_leniency)

# Sort by judges then probabilities in decreasing order
# Most dangerous for each judge are at the top.
compas_shuffled.sort_values(by=["judge_ID", "decile_score"],
                            ascending=False,
                            inplace=True)

# Iterate over the data. Subject will be given a negative decision
# if they are in the top (1-r)*100% of the individuals the judge will judge.
# I.e. if their within-judge-index is under 1 - acceptance threshold times
# the number of subjects assigned to each judge they will receive a
# negative decision.
compas_shuffled.reset_index(drop=True, inplace=True)

subjects_allocated = compas_shuffled.judge_ID.value_counts()

compas_shuffled['judge_index'] = compas_shuffled.groupby('judge_ID').cumcount()

compas_shuffled['decision_T'] = np.where(
    compas_shuffled['judge_index'] < (1 - compas_shuffled['judge_leniency']) *
    subjects_allocated[compas_shuffled['judge_ID']].values, 0, 1)

compas_labeled = compas_shuffled.copy()
compas_unlabeled = compas_shuffled.copy()

# Hide unobserved
compas_labeled.loc[compas_labeled.decision_T == 0, 'result_Y'] = np.nan

# Choose feature_columns
feature_cols = ~compas_labeled.columns.isin(
    ['result_Y', 'decile_score', 'judge_ID', 'judge_leniency', 'judge_index', 'decision_T'])

feature_cols = compas_labeled.columns[feature_cols]


####### Draw figures ###########

failure_rates = np.zeros((8, 6))
failure_sems = np.zeros((8, 6))

f_rate_true = np.zeros((nIter, 8))
f_rate_label = np.zeros((nIter, 8))
f_rate_label_adj = np.zeros((nIter, 8))
f_rate_human = np.zeros((nIter, 8))
f_rate_cont = np.zeros((nIter, 8))
f_rate_caus = np.zeros((nIter, 8))


# Split data
train, test = train_test_split(compas_labeled, test_size=0.5)

# Train a logistic regression model
B_model, predictions = fitPredictiveModel(
    train.loc[train['decision_T'] == 1, feature_cols],
    train.loc[train['decision_T'] == 1, 'result_Y'], test[feature_cols], 0)

test = test.assign(B_prob_0_model=predictions)

test.sort_values(by='B_prob_0_model', inplace=True, ascending=True)

kk_array = pd.qcut(test['B_prob_0_model'], group_amount, labels=False)

# Find observed values
observed = test['decision_T'] == 1

# Assign data to the model
dat = dict(D=1,
           N_obs=np.sum(observed),
           N_cens=np.sum(~observed),
           K=group_amount,
           sigma_tau=sigma_tau,
           M=len(set(compas_labeled['judge_ID'])),
           jj_obs=test.loc[observed, 'judge_ID']+1,
           jj_cens=test.loc[~observed, 'judge_ID']+1,
           kk_obs=kk_array[observed]+1,
           kk_cens=kk_array[~observed]+1,
           dec_obs=test.loc[observed, 'decision_T'],
           dec_cens=test.loc[~observed, 'decision_T'],
           X_obs=test.loc[observed, 'B_prob_0_model'].values.reshape(-1,1),
           X_cens=test.loc[~observed, 'B_prob_0_model'].values.reshape(-1,1),
           y_obs=test.loc[observed, 'result_Y'].astype(int))

sm = pystan.StanModel(file=stan_code_file_name)

fit = sm.sampling(data=dat, chains=5, iter=4000, control = dict(adapt_delta=0.9))

pars = fit.extract()

print(fit,  file=open(save_name + '_stan_fit_diagnostics.txt', 'w'))

plt.figure(figsize=(15,30))

fit.plot();

plt.savefig(save_name + '_stan_diagnostic_plot')

plt.show()
plt.close('all')

# Bayes

# Alusta matriisi, rivillä yksi otos posteriorista
# sarakkeet havaintoja
y_imp = np.ones((pars['y_est'].shape[0], test.shape[0]))

# Täydennetään havaitsemattomat estimoiduilla
y_imp[:, ~observed] = 1-pars['y_est']

# Täydennetään havaitut havaituilla
y_imp[:, observed] = 1-test.loc[observed, 'result_Y']

Rs = np.arange(.1, .9, .1)

to_release_list = np.round(test.shape[0] * Rs).astype(int)

f_rate_bayes = np.full((pars['y_est'].shape[0], 8), np.nan)

for i in range(len(to_release_list)):
    est_failure_rates = np.sum(y_imp[:, 0:to_release_list[i]], axis=1) / test.shape[0]
    
    f_rate_bayes[:, i] = est_failure_rates
    
    failure_rates[i, 5] = np.mean(est_failure_rates)    
        
for i in range(nIter):
    
    print(" [", i, "] ", sep='', end="")

    for r in np.arange(1, 9):

        print(".", end="")

        # True evaluation

        f_rate_true[i, r - 1] = trueEvaluationEvaluator(
            compas_unlabeled, feature_cols, 'decision_T', 'result_Y', r / 10)

        # Labeled outcomes only

        f_rate_label[i, r - 1] = labeledOutcomesEvaluator(
            compas_labeled, feature_cols, 'decision_T', 'result_Y', r / 10)

        # Adjusted labeled outcomes

        f_rate_label_adj[i, r - 1] = labeledOutcomesEvaluator(
            compas_labeled,
            feature_cols,
            'decision_T',
            'result_Y',
            r / 10,
            adjusted=True)

        # Human evaluation

        f_rate_human[i, r - 1] = humanEvaluationEvaluator(
            compas_labeled, 'judge_ID', 'decision_T', 'result_Y',
            'judge_leniency', r / 10)

        # Contraction

        f_rate_cont[i, r - 1] = contractionEvaluator(
            compas_labeled, feature_cols, 'judge_ID', 'decision_T', 'result_Y',
            'judge_leniency', r / 10)

failure_rates[:, 0] = np.mean(f_rate_true, axis=0)
failure_rates[:, 1] = np.mean(f_rate_label, axis=0)
failure_rates[:, 2] = np.mean(f_rate_label_adj, axis=0)
failure_rates[:, 3] = np.mean(f_rate_human, axis=0)
failure_rates[:, 4] = np.mean(f_rate_cont, axis=0)

failure_sems[:, 0] = scs.sem(f_rate_true, axis=0)
failure_sems[:, 1] = scs.sem(f_rate_label, axis=0)
failure_sems[:, 2] = scs.sem(f_rate_label_adj, axis=0)
failure_sems[:, 3] = scs.sem(f_rate_human, axis=0)
failure_sems[:, 4] = scs.sem(f_rate_cont, axis=0)
failure_sems[:, 5] = scs.sem(f_rate_bayes, axis=0, nan_policy='omit')

x_ax = np.arange(0.1, 0.9, 0.1)

labels = [
    'True Evaluation', 'Labeled outcomes', 'Labeled outcomes, adj.',
    'Human evaluation', 'Contraction', 'Potential outcomes'
]
colours = ['g', 'magenta', 'darkviolet', 'r', 'b', 'c']

for i in range(failure_rates.shape[1]):
    plt.errorbar(x_ax,
                 failure_rates[:, i],
                 label=labels[i],
                 c=colours[i],
                 yerr=failure_sems[:, i])

plt.title('Failure rate vs. Acceptance rate')
plt.xlabel('Acceptance rate')
plt.ylabel('Failure rate')
plt.legend()
plt.grid()

plt.savefig(save_name + '_all')

plt.show()

print("\nFailure rates:")
print(np.array2string(failure_rates, formatter={'float_kind':lambda x: "%.5f" % x}))

print("\nMean absolute errors:")
for i in range(1, failure_rates.shape[1]):
    print(
        labels[i].ljust(len(max(labels, key=len))),
        np.round(
            np.mean(np.abs(failure_rates[:, 0] - failure_rates[:, i])), 5))
 
    
# Draw diagnostic figures
f_rates= [f_rate_true, f_rate_label, f_rate_label_adj, f_rate_human, 
          f_rate_cont, f_rate_bayes]

cols = 2
rows = np.ceil(len(f_rates) / cols)

plt.figure(figsize=(16, 4.5*rows+1))

ax = plt.subplot(rows, cols, 1)
x_ax = np.arange(1, 9, 1) / 10

plt.boxplot(f_rates[0], labels=x_ax)

plt.title(labels[0])
plt.xlabel('Acceptance rate')
plt.ylabel('Failure rate')
plt.grid()

for i in range(len(f_rates)):
    plt.subplot(rows, cols, i + 1, sharey=ax)

    plt.boxplot(f_rates[i], labels=x_ax)

    plt.title(labels[i])
    plt.xlabel('Acceptance rate')
    plt.ylabel('Failure rate')
    plt.grid()

plt.tight_layout()
plt.subplots_adjust(top=0.89)

title = "COMPAS data set\nFluctuation of failure rate estimates per method"

plt.suptitle(title, y=0.96, weight='bold')

plt.savefig(save_name + '_diagnostic_plot')

plt.show()