<jupyter_start><jupyter_text>---- Tutorial for Evaluating Fairness in Regression ---- Evaluating Your Own Model? Looking for a Template? **Please note: This notebook is intended as an introduction to and reference for evaluating the fairness of a set of machine learning models using FairMLHealth, and it is not designed to work with outside models.**  We recommend using our [Template-RegressionAssessment](../templates/Template-RegressionAssessment.ipynb) if you are looking for a template into which you can insert existing models.  OverviewThis notebook introduces concepts, methods, and libraries for measuring fairness in machine learning (ML) models as it relates to problems in healthcare. There are abundant other publications covering the theoretical basis for fairness metrics, with many online and academic resources covering the details of specific fairness measures (See [References (bottom)](references) and [Additional Resources (bottom)](additional_resources), or [Our Resources Page](../docs/Measures_QuickReference.md)). Many of these excellent references fall short of discussing edge cases as well as the practical and philosophical considerations when evaluating real models for real customers. Here we attempt to bridge this gap. Table of Contents[Part 1](part1) - Framing the Problem  [Part 2](part2) - Evaluating the Baseline Model[Part 3](part3) - Comparing Models   [References](references)   RequirementsTo run this notebook, please install FairMLHealth using [the instructions posted in GitHub](https://github.com/KenSciResearch/fairMLHealthinstallation_instructions). Some components of this notebook additionally require the [Fairlearn](https://github.com/fairlearn/fairlearn) package.The tutorial uses data from the MIMIC III Critical Care database, a freely accessible source of electronic health records from Beth Israel Deaconess Medical Center in Boston. To download the MIMIC III data, please use this link: [Access to MIMIC III](https://mimic.physionet.org/gettingstarted/access/) and save the data with the default directory name ("MIMIC"). No further action is required beyond remembering the download location, and you do not need to unzip any files.A basic knowledge of ML implementation in Python is assumed.  --------  Part 1 - Framing the Problem<jupyter_code># Verify that all required packages are present
from fairmlhealth.__validation import validate_tutorial_requirements
validate_tutorial_requirements()

# fairMLHealth tools
from fairmlhealth import measure, report, stat_utils

# Loaders specific to fairMLHealth tutorials
from fairmlhealth.__mimic_data import load_mimic3_example, feature_table
# Standard Libraries
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

# Load Prediction Libraries
import sklearn.metrics as sk_metric
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

# Remove limit to the number of columns and column widths displayed by pandas
pd.set_option('display.max_columns', None)
pd.set_option('display.max_colwidth', 0) # adjust column width as needed<jupyter_output><empty_output><jupyter_text>Loading the MIMIC III Data Subset As mentioned above, the MIMIC-III data download contains a folder of zipped_files. The tutorial code in the [Data Subset](datasubset) section below will automatically unzip and format all necessary data for these experiments, saving the formatted data in the same MIMIC folder. Simply enter the correct path of the MIMIC folder in the following cell to enable this feature. Your path should end with the directory "MIMIC".Example: path_to_mimic_data_folder = "~/data/MIMIC"<jupyter_code># path_to_mimic_data_folder = "[path to your downloaded data folder]"
path_to_mimic_data_folder = "~/data/MIMIC"<jupyter_output><empty_output><jupyter_text>Example models in this notebook use data from all years of the MIMIC-III dataset for patients aged 65 and older. Data are imported at the encounter level with all additional patient identification dropped. All models include an "AGE" feature, simplified to 5-year bins, as well as Boolean diagnosis and procedure features categorized through the Clinical Classifications Software system ([HCUP](https://www.hcup-us.ahrq.gov/toolssoftware/ccs/ccs.jsp)). All features other than age are one-hot encoded and prefixed with their variable type (e.g. "GENDER_", "ETHNICITY_").<jupyter_code># Load data and keep a 10K observation subset to speed processing
df = load_mimic3_example(path_to_mimic_data_folder) 
df = df.sample(n=10000, random_state=32).reset_index(drop=True)

# Subset to ages 65+
df = df.loc[df['AGE'].ge(65), :]
df.drop('GENDER_F', axis=1, inplace=True) # Redundant with GENDER_M

# Display insights
display(feature_table(df))
print("\n\n", "Below is a scrollable version of the first five rows of data:")
display(df.head())<jupyter_output>This data subset has 4704 total observations and 648 input features<jupyter_text>Generating the Model<jupyter_code># Display LOS Distributions
ax = df.groupby('LANGUAGE_ENGL')['length_of_stay'
        ].plot(kind='kde', 
               title="Distribution of Length of Stay (in Days) Relative to English Language Status")
plt.show()
# Subset Data
baseline_cols = [c for c in df.columns if (c.startswith('AGE') or c.startswith('DIAGNOSIS_') or c.startswith('PROCEDURE_') or c == "LANGUAGE_ENGL")]
X = df.loc[:, baseline_cols]
y = df['length_of_stay'].clip(0, 25)

# Split Data
splits = train_test_split(X, y, test_size=0.33, random_state=42)
X_train, X_test, y_train, y_test = splits
# Train New Model with Language Feature
baseline_model = LinearRegression()
baseline_model.fit(X_train, y_train)
y_pred = baseline_model.predict(X_test)
print("Max Y before bounding:", np.max(y_pred))


# Since the regression is unbounded, clip all values greater than 25 
class bounded_baseline(object):
    def predict(X):
        return np.clip(baseline_model.predict(X), 0, 25)

y_pred = bounded_baseline.predict(X_test)
print("Max Y after bounding:", np.max(y_pred))<jupyter_output>Max Y before bounding: 2938873593794.152
Max Y after bounding: 25.0<jupyter_text>FairMLHealth Stratified Data TableFairMLHealth includes stratified table features to aid in identifying the source of unfairness or other bias. The data table can be used to compare basic statistics specific to each feature-value, in addition to relative statistics for the target value. Since many features can be evaluated at once, it can be a useful option for identifying patterns of bias either alone or in concert with other methods (e.g., visualization).Our experiment tests sociodemographic bias as it relates to language and age. Language may also be used as a proxy for race, religion, or nationality, which are [explicitly protected attributes](https://www.eeoc.gov/employers/small-business/3-who-protected-employment-discrimination). In this example, we can see that there are similar sample sizes (within 6%) for the English speaking and Non-English speaking groups, where the mean length of stay for Non-English speakers is 10 days (vs 8.7 days for English speakers). Similarly, for age, most groups have a similar number of members with the exceptions of the 85 y.o. and 90 y.o. groups.<jupyter_code>print("FairMLHealth Stratified Data Table")
measure.data(X_test, y_test, features=['LANGUAGE_ENGL', 'AGE'])<jupyter_output>FairMLHealth Stratified Data Table<jupyter_text>FairMLHealth Stratified Performance TableThe stratified performance table compares model performance specific to each feature-value subset. In this example, we can see that the MAE for the Non-English speaking group is higher than that of the English-speaking group.<jupyter_code>print("FairMLHealth Stratified Performance Table")
measure.performance(X_test, y_test, y_pred, features=['LANGUAGE_ENGL', 'AGE'], pred_type="regression")<jupyter_output>FairMLHealth Stratified Performance Table<jupyter_text>FairMLHealth Stratified Bias TableThe stratified bias table contains prediction bias measures specific to each feature-value subset. It assumes each feature-value as the "privileged" group relative to all other possible values for the feature. For this example, let's assume conservative "fair" boundaries for our error. Since the default boundaries for MAE difference are 10% of the range of absolute error, we define a more conservative custom boundary of +/- 10% of the MAE. Similarly, for the MAE Ratio, we set is as 1 +/- 10%.<jupyter_code># The default boundaries for MAE difference are 10% of the range of absolute error. For this example,
#    we've set it to a more conservative custom boundary of +/- 10% of the MAE
mae = sk_metric.mean_absolute_error(y_test, y_pred)
max_diff = 0.1*mae
custom_fair_range = {'MAE Difference':(-max_diff, max_diff), 'MAE Ratio':(0.9, 1.1)}
print("FairMLHealth Stratified Bias Table")

# Note that flag_oor is set to False by default for this feature
measure.bias(X_test, 
             y_test, 
             y_pred, 
             features=['LANGUAGE_ENGL', 'AGE'], 
             pred_type="regression",
             flag_oor=True, 
             custom_ranges=custom_fair_range
            )<jupyter_output>FairMLHealth Stratified Bias Table<jupyter_text>We may also want to evaluate the intersectional fairness of our model relative to both attributes at once. We can see that differences in language hold across most age groups, and are particularly pronounced in the oldest two groups.<jupyter_code>measure.bias(X_test, 
             y_test, 
             y_pred, 
             features=['LANGUAGE_ENGL'], 
             pred_type="regression",
             flag_oor=True, 
             custom_ranges=custom_fair_range,
             add_overview=False,
             cohort_labels=X_test['AGE'].astype(int)
            )<jupyter_output><empty_output><jupyter_text>Significance TestingIt is generally recommended to test whether any differences in model outcomes for protected attributes are the effect of a sampling error in our test. FairMLHealth comes with a bootstrapping utility and supporting functions that can be used in statistical testing. The bootstrapping feature accepts any function that returns a p-value and will return a True or False if the p-value is greater than some alpha for a threshold number of randomly sampled trials. Three examples using the bootstrap_significance tool with built-in test functions are shown below. The first and second tests indicate that there is likely a statistically significant difference in the of length of stay between English speakers and Non-English speakers. The third and fourth test indicate no likely significant difference either in the predicted distributions or in the distributions of prediction error, even when using a conservative threshold of 50% (i.e., at least 50% of trials must have a p-value less than alpha to indicate a difference). However, these tests are not comprehensive. For the purposes of the tutorial, we'll proceed as though we're convinced of a statistical difference.<jupyter_code>speaks_engl = df['LANGUAGE_ENGL'].eq(1)

# Krusakal-Wallis (K-W) to test the likelihood of statistical difference in the mean LOS
# relative to recorded English-speaking ability in the input data
pval = stat_utils.kruskal_pval(a=df.loc[speaks_engl, 'length_of_stay'], 
                               b=df.loc[~speaks_engl, 'length_of_stay'], 
                               # If n_sample set to None, tests on full dataset rather than sample
                               n_sample=None 
                              )
print("P-Value of single K-W test:", pval)

# Bootstrapped K-W test for the same data
alpha = 0.05
reject_h0 = stat_utils.bootstrap_significance(func=stat_utils.kruskal_pval, 
                                              a=y.loc[speaks_engl], 
                                              b=y.loc[~speaks_engl],
                                              alpha=alpha,
                                              n_trials=100)
print(f"Was the p-value less than alpha in at least 1-alpha ({1-alpha}) proportion of trials?", reject_h0)
# Bootstrapped K-W Test for true values
test_speaks_engl = X_test['LANGUAGE_ENGL'].eq(1)
reject_h0 = stat_utils.bootstrap_significance(func=stat_utils.kruskal_pval, 
                                              a=y_test.loc[test_speaks_engl], 
                                              b=y_test.loc[~test_speaks_engl],
                                              alpha=0.05,
                                              n_trials=1000)
print("Can we reject the null hypothesis that targets are from the same distribution?\n", reject_h0)
# Bootstrapped K-W Test for predictions
reject_h0 = stat_utils.bootstrap_significance(func=stat_utils.kruskal_pval, 
                                              a=np.ma.masked_array(y_pred, test_speaks_engl), 
                                              b=np.ma.masked_array(y_pred, ~test_speaks_engl),
                                              alpha=0.05,
                                              threshold=0.5,
                                              n_trials=100)
print("Can we reject the null hypothesis that prediction results are from the same distribution?\n", reject_h0)
# Bootstrapped K-W Test for prediction error
reject_h0 = stat_utils.bootstrap_significance(func=stat_utils.kruskal_pval, 
                                              a=np.ma.masked_array(y_pred-y_test, test_speaks_engl), 
                                              b=np.ma.masked_array(y_pred-y_test, ~test_speaks_engl),
                                              alpha=0.05,
                                              threshold=0.5,
                                              n_trials=100)
print("Can we reject the null hypothesis that prediction errors are from the same distribution for language?\n", reject_h0)
# Bootstrapped K-W Test for prediction error
test_is90plus = X_test['AGE'].eq(90).astype(int)
reject_h0 = stat_utils.bootstrap_significance(func=stat_utils.kruskal_pval, 
                                              a=np.ma.masked_array(y_pred-y_test, test_is90plus), 
                                              b=np.ma.masked_array(y_pred-y_test, ~test_is90plus),
                                              alpha=0.05,
                                              threshold=0.5,
                                              n_trials=100)
print("Can we reject the null hypothesis that prediction errors are from the same distribution for age?\n", reject_h0)<jupyter_output>Can we reject the null hypothesis that prediction errors are from the same distribution for age?
 False<jupyter_text>-------- Part 2 - Comparing Measures In the scenario where a statistical difference was found relative to language or relative to age 90+, we may want to report our findings to our team or to our stakeholders. Using fairMLHealth we can create a custom report comparing key indicators by passing integer tags for our protected attributes.Using our previously defined boundaries, we can see that the error in our test model is outside of what we've defined as an acceptable rage.<jupyter_code>print("FairMLHealth Fairness Measure Report: English Language")
report.compare(test_data = X_test, 
                targets = y_test, 
                protected_attr = {'English Language': X_test['LANGUAGE_ENGL'], 
                                  'Age 90+': X_test['AGE'].eq(90).astype(int)}, 
                models = bounded_baseline, 
                pred_type = "regression",
                custom_ranges = custom_fair_range
               )<jupyter_output>FairMLHealth Fairness Measure Report: English Language<jupyter_text>Special CasesFor many regression problems, the effect of the error may differ depending upon the result. To use a hypothetical example, suppose that those patients for whom the length of ICU stay will be greater than 10 hours are better served by special longer-term unit. Those paients whose length of stay is predicted to be longer than 12 hours are automatically transferred to this special unit. However, cases for which the predictions are between 4 and 12 hours are reviewed on a case-by-case basis where proximity and direction to the 8-hour threshol heavily weighted factors in the decision. Since the mean length of stay in our dataset is 9 hours, the prediction error will be  critical for patients with middling values. In these cases it can be useful to evaluate the model relative to each prediction cohort. For the example plots below, data are divided into three groups as shown in the table below.| Predicted LOS (Time Group) | Action ||- |- ||0 - 6 hrs | No Action ||6 - 12 hrs | Flagged for Review; Possible Transfer to Special Unit ||12+ hrs | Automatic Transfer to Special Unit |<jupyter_code># Generate cohorts
cohorts = pd.Series(np.zeros(len(y_test)), index=y_test.index, name="Time Group")
cohorts.loc[y_test.ge(6) & y_test.lt(12)] = 1
cohorts.loc[y_test.ge(12)]  = 2
cohorts = cohorts.map({0:"0-6 hrs", 1:"6-12 hrs", 2:"12+ hrs"})

# Append cohorts to data
cohort_data = pd.concat([cohorts, y_test.rename('Length of Stay')], axis=1)
cohort_data['Error'] = y_pred-y_test<jupyter_output><empty_output><jupyter_text>As we can observe in the table and figure below, the mean error for each group leans in the direction of the population mean, which was shown above to be approximately 9 hours. The model is over-predicting for the 0-6 hour time group and under-predicting for the 12+ hour time group, but errors appear to be small enough that predictions remain within their relative groups.<jupyter_code>cohort_data.groupby('Time Group')[['Length of Stay', 'Error']].mean()
# Generate plots
g = sns.lineplot(x=cohorts, y=y_test)
sns.lineplot(x=cohorts, y=pd.Series(y_pred-y_test, name="Error")) # Adds line for the error
g.axhspan(0, 4, alpha=0.25, color='green') 
g.axhspan(4, 12, alpha=0.25, color='lightyellow')
g.axhspan(12, 25, alpha=0.25, color='blue')
plt.xticks([*range(3)], labels=["0 - 6 hrs", "6 - 12 hrs", "12+ hrs"])

g.set_ylabel("Length of Stay")
plt.legend(labels=["True LOS", "Error"])
plt.show()<jupyter_output><empty_output><jupyter_text>Using the bias table, we can verify that any between-group differences in the observed error are within the fair range for each group.<jupyter_code>measure.bias(X_test, 
             y_test, 
             y_pred, 
             features=['LANGUAGE_ENGL'], 
             pred_type="regression",
             cohort_labels=cohorts,
             add_overview=False,
             flag_oor=False
            )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0723804406233208, 'fraction_numerical': 0.012251477700161204, 'mean_word_length': 3.532635168046761, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 52, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '318888', 'license': 'MIT', 'n_tokens_mistral': 5212, 'n_tokens_neox': 4823, 'n_words': 1962, 'path': 'examples_and_tutorials/Tutorial-EvaluatingFairnessInRegression.ipynb', 'repo_name': 'masino-lab/fairMLHealth', 'size': '18610'}