stan_modelling_empirical_SEfixed_useprediction.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) group_amount : How many groups if Jung-inspired model is used.
# (3) stan_code_file_name : Name of file containing the stan model code.
# (4) sigma_tau : Values of prior variance for the Jung-inspired model.
# (5) data_path : File of compas data.

"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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]

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

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

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

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

# Prefix for the figures and log files

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

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)

def standardError(p, n):
    denominator = p * (1 - p)
    
    return np.sqrt(denominator / n)

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

# ## 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], D_q.shape[0]


# ### Evaluators

def trueEvaluationEvaluator(df, resultY_col, r):

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

    to_release = int(round(df.shape[0] * r))
    
    failed = df[resultY_col][0:to_release] == 0

    return np.sum(failed) / df.shape[0], df.shape[0]


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

    df_observed = df.loc[df[decisionT_col] == 1, :]

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

    to_release = int(round(df_observed.shape[0] * r))
    
    failed = df_observed[resultY_col][0:to_release] == 0

    if adjusted:
        return np.mean(failed), df.shape[0]

    return np.sum(failed) / df.shape[0], df.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, 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)]

    failed = released[resultY_col] == 0
    
    # Get their failure rate, aka ratio of reoffenders to number of people judged in total
    return np.sum(failed) / correct_leniency_list.shape[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))
error_Ns = np.zeros((8, 6))

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

# Assign same observations to unlabeled data
test_unlabeled = compas_unlabeled.iloc[test_labeled.index.values]


# 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_labeled[feature_cols], 0)

test_labeled = test_labeled.assign(B_prob_0_model=predictions)
test_unlabeled = test_unlabeled.assign(B_prob_0_model=predictions)

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

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

# Find observed values
observed = test_labeled['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_labeled.loc[observed, 'judge_ID']+1,
           jj_cens=test_labeled.loc[~observed, 'judge_ID']+1,
           kk_obs=kk_array[observed]+1,
           kk_cens=kk_array[~observed]+1,
           dec_obs=test_labeled.loc[observed, 'decision_T'],
           dec_cens=test_labeled.loc[~observed, 'decision_T'],
           X_obs=test_labeled.loc[observed, 'B_prob_0_model'].values.reshape(-1, 1),
           X_cens=test_labeled.loc[~observed, 'B_prob_0_model'].values.reshape(-1, 1),
           y_obs=test_labeled.loc[observed, 'result_Y'].astype(int))

sm = pystan.StanModel(file=stan_code_file_name)

fit = sm.sampling(data=dat, chains=5, iter=6000, 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_labeled.shape[0]))

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

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

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

to_release_list = np.round(test_labeled.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)):
    
    failed = np.sum(y_imp[:, 0:to_release_list[i]], axis=1)
    
    est_failure_rates =  failed / test_labeled.shape[0]
            
    failure_rates[i, 5] = np.mean(est_failure_rates)
    
    error_Ns[i, 5] = test_labeled.shape[0]

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

    print(".", end="")

    # True evaluation

    FR, N = trueEvaluationEvaluator(test_unlabeled, 'result_Y', r / 10)
    
    failure_rates[r - 1, 0] = FR
    error_Ns[r - 1, 0] = N

    # Labeled outcomes only

    FR, N = labeledOutcomesEvaluator(test_labeled, 'decision_T', 'result_Y', r / 10)

    failure_rates[r - 1, 1] = FR
    error_Ns[r - 1, 1] = N
    
    # Adjusted labeled outcomes

    FR, N = labeledOutcomesEvaluator(test_labeled, 'decision_T', 'result_Y', r / 10,
                                     adjusted=True)

    failure_rates[r - 1, 2] = FR
    error_Ns[r - 1, 2] = N

    # Human evaluation

    FR, N = humanEvaluationEvaluator(test_labeled, 'judge_ID', 'decision_T', 
                                     'result_Y',  'judge_leniency', 
                                     r / 10)

    failure_rates[r - 1, 3] = FR
    error_Ns[r - 1, 3] = N

    # Contraction

    FR, N = contraction(test_labeled, 'judge_ID', 'decision_T', 'result_Y', 
                        'B_prob_0_model', 'judge_leniency', r / 10)

    failure_rates[r - 1, 4] = FR
    error_Ns[r - 1, 4] = N

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

failure_SEs = standardError(failure_rates, error_Ns)

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_SEs[:, 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))