From 91cb7476456e1d407ef66d4461988a1bb26a886a Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:06:25 +0100 Subject: [PATCH 1/6] Update docs to make it clear the decision boundary is an argument --- daindex/main.py | 2 +- daindex/util.py | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/daindex/main.py b/daindex/main.py index 12822cd..3bd1510 100644 --- a/daindex/main.py +++ b/daindex/main.py @@ -160,7 +160,7 @@ class DAIndex(object): weight_sum_steps: The number of bins for the weighted sum of k-step cutoffs. n_jobs: The number of jobs to run in parallel. model_name: The name of the model to use in the plots. - decision_boundary: The decision boundary for the DA curve. + decision_boundary: The decision boundary for the DA curve. Defaults to 0.5. Methods: setup_groups: Set up the groups for the DAI calculation. diff --git a/daindex/util.py b/daindex/util.py index 9e3ce79..4d7045b 100644 --- a/daindex/util.py +++ b/daindex/util.py @@ -86,8 +86,13 @@ def compare_two_groups( def area_under_curve(w_data: np.ndarray, decision_boundary: float = 0.5) -> tuple[float, float]: - """ - calculate the area under curve - do NOT do interpolation + """Calculate the area under curve - do NOT do interpolation + + Args: + w_data: The weighted data, should be a 2D array with shape (n, 3) where + the first column is the x values, the second column is the y values, + and the third column is the weights. + decision_boundary: The decision boundary for calculating the area under curve. """ prev = None area = 0.0 From ddedb960440267e03832b1727f6c659eac92e3f0 Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:19:10 +0100 Subject: [PATCH 2/6] Small improvement to docs --- daindex/main.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/daindex/main.py b/daindex/main.py index 3bd1510..12e21d7 100644 --- a/daindex/main.py +++ b/daindex/main.py @@ -141,8 +141,8 @@ class DAIndex(object): `evaluate_all_groups_by_models` methods can be called to calculate the DAI for all groups relative to a single reference group. - The results can be accessed using the `get_group_ratios` and `get_group_figures` methods, or printed using - the `present_results` and `present_all_results` methods. + The results can be accessed using the `get_group_ratio` and `get_group_figure` methods (with a reference and + other group), or printed using the `present_results` and `present_all_results` methods. Args: cohort: The DataFrame containing the data for the cohort. @@ -174,8 +174,8 @@ class DAIndex(object): model objects. present_results: Print the DAI results for a pair of groups. present_all_results: Print the DAI results for all group pairs. - get_group_ratios: Get the DAI ratios for a pair of groups. - get_group_figures: Get the DAI figures for a pair of groups. + get_group_ratio: Get the DAI ratios for a pair of groups. + get_group_figure: Get the DAI figure for a pair of groups. get_all_ratios: Get the DAI ratios for all group pairs. get_all_figures: Get the DAI figures for all group pairs. """ @@ -745,10 +745,10 @@ def present_all_results(self) -> None: for group_pair in self.group_figures.keys(): self.present_results(*group_pair) - def get_group_ratios(self, reference_group: str, other_group: str) -> dict[str, float | str]: + def get_group_ratio(self, reference_group: str, other_group: str) -> dict[str, float | str]: return self.group_ratios[(reference_group, other_group)] - def get_group_figures(self, reference_group: str, other_group: str) -> plt.Figure: + def get_group_figure(self, reference_group: str, other_group: str) -> plt.Figure: return self.group_figures[(reference_group, other_group)] def get_all_ratios(self) -> pd.DataFrame: From 1b0db7ff695f4a44c4016ba545031af25ba779c9 Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:46:45 +0100 Subject: [PATCH 3/6] Add methods to return all, failed and sub-optimal steps per group --- daindex/main.py | 98 +++++++++ pyproject.toml | 2 +- tutorials/meps_tutorial.ipynb | 31 ++- tutorials/meps_tutorial.py | 389 ++++++++++++++++++++++++++++++++++ uv.lock | 2 +- 5 files changed, 517 insertions(+), 5 deletions(-) create mode 100644 tutorials/meps_tutorial.py diff --git a/daindex/main.py b/daindex/main.py index 12e21d7..c252890 100644 --- a/daindex/main.py +++ b/daindex/main.py @@ -178,6 +178,12 @@ class DAIndex(object): get_group_figure: Get the DAI figure for a pair of groups. get_all_ratios: Get the DAI ratios for all group pairs. get_all_figures: Get the DAI figures for all group pairs. + get_group_step_samples: Get the number of samples used for each step for a specific group. + get_group_sub_optimal_steps: Get information about sub-optimal steps for a specific group. + get_group_failed_steps: Get information about failed steps for a specific group. + get_all_groups_step_samples: Get step sample information for all evaluated groups. + get_all_groups_sub_optimal_steps: Get sub-optimal step information for all evaluated groups. + get_all_groups_failed_steps: Get failed step information for all evaluated groups. """ def __init__( @@ -221,6 +227,10 @@ def __init__( self.group_ratios = {} self.group_figures = {} self.issues = [] + # New attributes to track step information + self.group_sub_optimal_steps = {} # Stores sub-optimal steps for each group + self.group_failed_steps = {} # Stores failed steps for each group + self.group_step_samples = {} # Stores sample counts for each step for each group def setup_groups(self, groups: Group | list[Group], group_col: str = None) -> None: if not isinstance(groups, list): @@ -485,6 +495,16 @@ def process_step(s: int) -> tuple[float, int, float, bool, bool]: delayed(process_step)(s) for s in tqdm(self.step_scores, desc=f"Calculating k-steps for '{group}' group", position=1, leave=False) ) + + # Store step information internally + sub_opt_steps = {s[0]: s[1] for s in ret if s[3]} + failed_steps = {s[0]: s[1] for s in ret if s[4]} + step_samples = {s[0]: s[1] for s in ret} + + self.group_sub_optimal_steps[group] = sub_opt_steps + self.group_failed_steps[group] = failed_steps + self.group_step_samples[group] = step_samples + sub_opt_list = ", ".join([f"{s[0]}: {s[1]}" for s in ret if s[3]]) failed_list = ", ".join([f"{s[0]}: {s[1]}" for s in ret if s[4]]) if sub_opt_list or failed_list: @@ -758,3 +778,81 @@ def get_all_ratios(self) -> pd.DataFrame: def get_all_figures(self) -> dict[tuple[str, str], plt.Figure]: return self.group_figures + + def get_group_step_samples(self, group: str) -> dict[float, int]: + """ + Returns the number of samples used for each step for the specified group. + + Args: + group: The name of the group to get step sample information for. + + Returns: + A dictionary mapping step scores to the number of samples used. + """ + if group not in self.group_step_samples: + raise ValueError(f"No step information available for group '{group}'. Run evaluation first.") + return self.group_step_samples[group] + + def get_group_sub_optimal_steps(self, group: str) -> dict[float, int]: + """ + Returns information about sub-optimal steps for the specified group. + These are steps where the number of samples was below the preferred threshold + but still acceptable. + + Args: + group: The name of the group to get sub-optimal step information for. + + Returns: + A dictionary mapping step scores to the number of samples used for + sub-optimal steps. + """ + if group not in self.group_sub_optimal_steps: + raise ValueError(f"No step information available for group '{group}'. Run evaluation first.") + return self.group_sub_optimal_steps[group] + + def get_group_failed_steps(self, group: str) -> dict[float, int]: + """ + Returns information about failed steps for the specified group. + These are steps where the number of samples was below the minimum threshold + and the step was excluded from the analysis. + + Args: + group: The name of the group to get failed step information for. + + Returns: + A dictionary mapping step scores to the number of samples used for + failed steps. + """ + if group not in self.group_failed_steps: + raise ValueError(f"No step information available for group '{group}'. Run evaluation first.") + return self.group_failed_steps[group] + + def get_all_groups_step_samples(self) -> dict[str, dict[float, int]]: + """ + Returns the number of samples used for each step for all evaluated groups. + + Returns: + A dictionary mapping group names to dictionaries of step scores and + their corresponding sample counts. + """ + return self.group_step_samples.copy() + + def get_all_groups_sub_optimal_steps(self) -> dict[str, dict[float, int]]: + """ + Returns information about sub-optimal steps for all evaluated groups. + + Returns: + A dictionary mapping group names to dictionaries of step scores and + their corresponding sample counts for sub-optimal steps. + """ + return self.group_sub_optimal_steps.copy() + + def get_all_groups_failed_steps(self) -> dict[str, dict[float, int]]: + """ + Returns information about failed steps for all evaluated groups. + + Returns: + A dictionary mapping group names to dictionaries of step scores and + their corresponding sample counts for failed steps. + """ + return self.group_failed_steps.copy() diff --git a/pyproject.toml b/pyproject.toml index 312f7b9..fd4e09c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,7 @@ authors = [ { name = "Harrison Wilde", email = "h.wilde@ucl.ac.uk" }, { name = "Honghan Wu", email = "honghan.wu@gmail.com" }, ] -version = "0.7.3" +version = "0.8.0" description = "Deterioration Allocation Index Framework" license = "MIT" readme = "README.md" diff --git a/tutorials/meps_tutorial.ipynb b/tutorials/meps_tutorial.ipynb index eb37326..27b8998 100644 --- a/tutorials/meps_tutorial.ipynb +++ b/tutorials/meps_tutorial.ipynb @@ -494,9 +494,13 @@ " cohort=df,\n", " groups=groups,\n", " det_feature=det_feature,\n", + " decision_boundary=0.75, # Set decision boundary for DA curve\n", ")\n", "index.evaluate_all_groups_by_models(\n", - " models=lr_orig_panel19, feature_list=feature_list, reference_group=\"White\", n_jobs=6\n", + " models=lr_orig_panel19,\n", + " feature_list=feature_list,\n", + " reference_group=\"White\",\n", + " n_jobs=1, # Up n_jobs to parallelize if needed\n", ")" ] }, @@ -518,6 +522,24 @@ "index.get_all_ratios()" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "index.get_group_ratios(\"White\", \"Asian\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "index.get_all_groups_failed_steps()" + ] + }, { "cell_type": "code", "execution_count": null, @@ -527,8 +549,11 @@ } ], "metadata": { + "jupytext": { + "formats": "ipynb,py:percent" + }, "kernelspec": { - "display_name": "daindex", + "display_name": ".venv", "language": "python", "name": "python3" }, @@ -542,7 +567,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.5" + "version": "3.10.18" } }, "nbformat": 4, diff --git a/tutorials/meps_tutorial.py b/tutorials/meps_tutorial.py new file mode 100644 index 0000000..82249f8 --- /dev/null +++ b/tutorials/meps_tutorial.py @@ -0,0 +1,389 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.17.2 +# kernelspec: +# display_name: .venv +# language: python +# name: python3 +# --- + +# %% [markdown] +# We will use the [`aif360`](https://aif360.readthedocs.io/en/latest/Getting%20Started.html) package to load the UCI adult dataset, fit a simple model and then analyse the fairness of the model using the DA-AUC. + +# %% +import numpy as np +from aif360.datasets import MEPSDataset19 +from aif360.explainers import MetricTextExplainer +from aif360.metrics import BinaryLabelDatasetMetric, ClassificationMetric +from sklearn.linear_model import LogisticRegression +from sklearn.pipeline import make_pipeline +from sklearn.preprocessing import StandardScaler + +np.random.seed(1) + +# %% [markdown] +# The below code downloads the required files from the UCI website if they are not already present in ai360's data directory. + +# %% +import os +import shutil +import subprocess + +import aif360 + +aif360_location = os.path.dirname(aif360.__file__) +meps_data_dir = os.path.join(aif360_location, "data", "raw", "meps") +h181_file_path = os.path.join(meps_data_dir, "h181.csv") + +if not os.path.isfile(h181_file_path): + r_script_path = os.path.join(meps_data_dir, "generate_data.R") + process = subprocess.Popen(["Rscript", r_script_path], stdin=subprocess.PIPE) + process.communicate(input=b"y\n") + + # Move the generated CSV files to meps_data_dir + generated_files = ["h181.csv", "h192.csv"] + for file_name in generated_files: + src_path = os.path.join(os.getcwd(), file_name) + dest_path = os.path.join(meps_data_dir, file_name) + if os.path.isfile(src_path): + shutil.move(src_path, dest_path) + + +# %% +def preprocessing_w_multimorb(df): + """ + 1.Create a new column, RACE that is 'White' if RACEV2X = 1 and HISPANX = 2 i.e. non Hispanic White + and 'non-White' otherwise + 2. Restrict to Panel 19 + 3. RENAME all columns that are PANEL/ROUND SPECIFIC + 4. Drop rows based on certain values of individual features that correspond to missing/unknown - generally < -1 + 5. Compute UTILIZATION, binarize it to 0 (< 10) and 1 (>= 10) + """ + + def race(row): + if (row["HISPANX"] == 2) and ( + row["RACEV2X"] == 1 + ): # non-Hispanic Whites are marked as WHITE; all others as NON-WHITE + return "White" + return "Non-White" + + df["RACE"] = df.apply(lambda row: race(row), axis=1) + + df = df[df["PANEL"] == 19] + + # RENAME COLUMNS + df = df.rename( + columns={ + "FTSTU53X": "FTSTU", + "ACTDTY53": "ACTDTY", + "HONRDC53": "HONRDC", + "RTHLTH53": "RTHLTH", + "MNHLTH53": "MNHLTH", + "CHBRON53": "CHBRON", + "JTPAIN53": "JTPAIN", + "PREGNT53": "PREGNT", + "WLKLIM53": "WLKLIM", + "ACTLIM53": "ACTLIM", + "SOCLIM53": "SOCLIM", + "COGLIM53": "COGLIM", + "EMPST53": "EMPST", + "REGION53": "REGION", + "MARRY53X": "MARRY", + "AGE53X": "AGE", + "POVCAT15": "POVCAT", + "INSCOV15": "INSCOV", + } + ) + + df = df[df["REGION"] >= 0] # remove values -1 + df = df[df["AGE"] >= 0] # remove values -1 + + df = df[df["MARRY"] >= 0] # remove values -1, -7, -8, -9 + + df = df[df["ASTHDX"] >= 0] # remove values -1, -7, -8, -9 + + df = df[ + ( + df[ + [ + "FTSTU", + "ACTDTY", + "HONRDC", + "RTHLTH", + "MNHLTH", + "HIBPDX", + "CHDDX", + "ANGIDX", + "EDUCYR", + "HIDEG", + "MIDX", + "OHRTDX", + "STRKDX", + "EMPHDX", + "CHBRON", + "CHOLDX", + "CANCERDX", + "DIABDX", + "JTPAIN", + "ARTHDX", + "ARTHTYPE", + "ASTHDX", + "ADHDADDX", + "PREGNT", + "WLKLIM", + "ACTLIM", + "SOCLIM", + "COGLIM", + "DFHEAR42", + "DFSEE42", + "ADSMOK42", + "PHQ242", + "EMPST", + "POVCAT", + "INSCOV", + ] + ] + >= -1 + ).all(1) + ] # for all other categorical features, remove values < -1 + + def utilization(row): + return row["OBTOTV15"] + row["OPTOTV15"] + row["ERTOT15"] + row["IPNGTD15"] + row["HHTOTD15"] + + df["TOTEXP15"] = df.apply(lambda row: utilization(row), axis=1) + lessE = df["TOTEXP15"] < 10.0 + df.loc[lessE, "TOTEXP15"] = 0.0 + moreE = df["TOTEXP15"] >= 10.0 + df.loc[moreE, "TOTEXP15"] = 1.0 + df["MULTIMORBIDITY"] = ( + df.filter(regex="DX$|CHBRON$|JTPAIN$").drop(columns=["ADHDADDX"]).apply(lambda x: (x == 1).sum(), axis=1) + ) + + df = df.rename(columns={"TOTEXP15": "UTILIZATION"}) + return df + + +# %% +dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test = MEPSDataset19( + custom_preprocessing=preprocessing_w_multimorb, + features_to_keep=[ + "REGION", + "AGE", + "SEX", + "RACE", + "RACEV2X", + "MARRY", + "FTSTU", + "ACTDTY", + "HONRDC", + "RTHLTH", + "MNHLTH", + "HIBPDX", + "CHDDX", + "ANGIDX", + "MIDX", + "OHRTDX", + "STRKDX", + "EMPHDX", + "CHBRON", + "CHOLDX", + "CANCERDX", + "DIABDX", + "JTPAIN", + "ARTHDX", + "ARTHTYPE", + "ASTHDX", + "ADHDADDX", + "PREGNT", + "WLKLIM", + "ACTLIM", + "SOCLIM", + "COGLIM", + "DFHEAR42", + "DFSEE42", + "ADSMOK42", + "PCS42", + "MCS42", + "K6SUM42", + "PHQ242", + "EMPST", + "POVCAT", + "INSCOV", + "UTILIZATION", + "PERWT15F", + "MULTIMORBIDITY", + ], +).split([0.5, 0.8], shuffle=True) + +sens_ind = 0 +sens_attr = dataset_orig_panel19_train.protected_attribute_names[sens_ind] + +unprivileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.unprivileged_protected_attributes[sens_ind]] +privileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.privileged_protected_attributes[sens_ind]] + + +# %% +def describe(train=None, val=None, test=None): + if train is not None: + print(train.features.shape) + if val is not None: + print(val.features.shape) + print(test.features.shape) + print(test.favorable_label, test.unfavorable_label) + print(test.protected_attribute_names) + print(test.privileged_protected_attributes, test.unprivileged_protected_attributes) + print(test.feature_names) + + +# %% +describe(dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test) + +# %% +metric_orig_panel19_train = BinaryLabelDatasetMetric( + dataset_orig_panel19_train, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups +) +explainer_orig_panel19_train = MetricTextExplainer(metric_orig_panel19_train) + +print(explainer_orig_panel19_train.disparate_impact()) + +# %% +dataset = dataset_orig_panel19_train +model = make_pipeline(StandardScaler(), LogisticRegression(solver="liblinear", random_state=1)) +fit_params = {"logisticregression__sample_weight": dataset.instance_weights} + +lr_orig_panel19 = model.fit(dataset.features, dataset.labels.ravel(), **fit_params) + +# %% +from collections import defaultdict + + +def test(dataset, model, thresh_arr): + try: + # sklearn classifier + y_val_pred_prob = model.predict_proba(dataset.features) + pos_ind = np.where(model.classes_ == dataset.favorable_label)[0][0] + except AttributeError: + # aif360 inprocessing algorithm + y_val_pred_prob = model.predict(dataset).scores + pos_ind = 0 + + metric_arrs = defaultdict(list) + for thresh in thresh_arr: + y_val_pred = (y_val_pred_prob[:, pos_ind] > thresh).astype(np.float64) + + dataset_pred = dataset.copy() + dataset_pred.labels = y_val_pred + metric = ClassificationMetric( + dataset, dataset_pred, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups + ) + + metric_arrs["bal_acc"].append((metric.true_positive_rate() + metric.true_negative_rate()) / 2) + metric_arrs["avg_odds_diff"].append(metric.average_odds_difference()) + metric_arrs["disp_imp"].append(metric.disparate_impact()) + metric_arrs["stat_par_diff"].append(metric.statistical_parity_difference()) + metric_arrs["eq_opp_diff"].append(metric.equal_opportunity_difference()) + metric_arrs["theil_ind"].append(metric.theil_index()) + + return metric_arrs + + +# %% +thresh_arr = np.linspace(0.01, 0.5, 50) +val_metrics = test(dataset=dataset_orig_panel19_val, model=lr_orig_panel19, thresh_arr=thresh_arr) +lr_orig_best_ind = np.argmax(val_metrics["bal_acc"]) + + +# %% +def describe_metrics(metrics, thresh_arr): + best_ind = np.argmax(metrics["bal_acc"]) + print("Threshold corresponding to Best balanced accuracy: {:6.4f}".format(thresh_arr[best_ind])) + print("Best balanced accuracy: {:6.4f}".format(metrics["bal_acc"][best_ind])) + disp_imp_at_best_ind = 1 - min(metrics["disp_imp"][best_ind], 1 / metrics["disp_imp"][best_ind]) + print("Corresponding 1-min(DI, 1/DI) value: {:6.4f}".format(disp_imp_at_best_ind)) + print("Corresponding average odds difference value: {:6.4f}".format(metrics["avg_odds_diff"][best_ind])) + print("Corresponding statistical parity difference value: {:6.4f}".format(metrics["stat_par_diff"][best_ind])) + print("Corresponding equal opportunity difference value: {:6.4f}".format(metrics["eq_opp_diff"][best_ind])) + print("Corresponding Theil index value: {:6.4f}".format(metrics["theil_ind"][best_ind])) + + +# %% +describe_metrics(val_metrics, thresh_arr) + +# %% +lr_orig_metrics = test( + dataset=dataset_orig_panel19_test, model=lr_orig_panel19, thresh_arr=[thresh_arr[lr_orig_best_ind]] +) + +# %% +describe_metrics(lr_orig_metrics, [thresh_arr[lr_orig_best_ind]]) + +# %% +df = dataset.convert_to_dataframe()[0] + +# %% +df["MULTIMORBIDITY"].hist() +df_add = df.copy() +df_add["RTHLTH"] = df_add.filter(regex="RTHLTH").idxmax(axis=1) + +# %% +df_add.boxplot(column="MULTIMORBIDITY", by="RTHLTH") + +# %% +from daindex import DAIndex, DeteriorationFeature, Group + +# %% [markdown] +# Let's first specify our deterioration feature and groups. + +# %% +det_feature = DeteriorationFeature(col="MULTIMORBIDITY", threshold=1, is_discrete=True) +groups = [ + Group("White", 1, "RACEV2X"), + Group( + "Non-White", [12, 5, 10, 2, 3, 6, 4], "RACEV2X", det_threshold=1.2 + ), # Modify deterioration threshold only for this group + Group("Black", 2, "RACEV2X"), + Group("Asian", [4, 5, 6, 10], "RACEV2X"), +] +print(det_feature) +print(groups[1]) + +# %% [markdown] +# We can call a group on the cohort to get the group representation in the cohort: + +# %% +print(groups[0](df).RACEV2X.value_counts()) +groups[0](df).head() + +# %% [markdown] +# We can now instantiate the `DAIndex` with the components above and evaluate all groups via the model we have made. Note this could also be done with a list of arbitrary models. The package also supports evaluation via an existing predictions column in the cohort, see `evaluate_all_groups_from_predictions`. We could also just compare two specific groups with `evaluate_group_pair_by_models` and `evaluate_group_pair_from_predictions`. + +# %% +feature_list = dataset.feature_names +index = DAIndex( + cohort=df, + groups=groups, + det_feature=det_feature, + decision_boundary=0.75, # Set decision boundary for DA curve +) +index.evaluate_all_groups_by_models( + models=lr_orig_panel19, + feature_list=feature_list, + reference_group="White", + n_jobs=1, # Up n_jobs to parallelize if needed +) + +# %% +index.present_all_results() + +# %% +index.get_all_ratios() + +# %% +index.get_group_ratios("White", "Asian") + +# %% diff --git a/uv.lock b/uv.lock index 0d1f821..eec866a 100644 --- a/uv.lock +++ b/uv.lock @@ -505,7 +505,7 @@ wheels = [ [[package]] name = "daindex" -version = "0.7.0" +version = "0.7.3" source = { editable = "." } dependencies = [ { name = "joblib" }, From 8f57c069f793c5010ec8e68b2df313d1d37ba181 Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:51:58 +0100 Subject: [PATCH 4/6] Clean up --- .pre-commit-config.yaml | 10 +++++----- tutorials/meps_tutorial.ipynb | 2 +- tutorials/meps_tutorial.py | 3 +++ uv.lock | 2 +- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6a171dd..1e3478a 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -11,14 +11,14 @@ repos: - id: check-yaml - id: check-toml - id: check-added-large-files -- repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.12.5 - hooks: - - id: ruff-check - - id: ruff-format - repo: https://github.com/nestauk/pre-commit-hooks rev: v1.2.0 hooks: - id: nbstripout-preserve-timestamp - id: jupytext-enforce-pairing - id: jupytext-smart-sync +- repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.12.5 + hooks: + - id: ruff-check + - id: ruff-format diff --git a/tutorials/meps_tutorial.ipynb b/tutorials/meps_tutorial.ipynb index 27b8998..8e61197 100644 --- a/tutorials/meps_tutorial.ipynb +++ b/tutorials/meps_tutorial.ipynb @@ -528,7 +528,7 @@ "metadata": {}, "outputs": [], "source": [ - "index.get_group_ratios(\"White\", \"Asian\")" + "index.get_group_ratio(\"White\", \"Asian\")" ] }, { diff --git a/tutorials/meps_tutorial.py b/tutorials/meps_tutorial.py index 82249f8..1f653f1 100644 --- a/tutorials/meps_tutorial.py +++ b/tutorials/meps_tutorial.py @@ -387,3 +387,6 @@ def describe_metrics(metrics, thresh_arr): index.get_group_ratios("White", "Asian") # %% +index.get_group_ratios("White", "Non-White") + +# %% diff --git a/uv.lock b/uv.lock index eec866a..652b2fa 100644 --- a/uv.lock +++ b/uv.lock @@ -505,7 +505,7 @@ wheels = [ [[package]] name = "daindex" -version = "0.7.3" +version = "0.8.0" source = { editable = "." } dependencies = [ { name = "joblib" }, From d16573b66bb09c96f0be2f3d1b54e7d8769bbf23 Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:52:55 +0100 Subject: [PATCH 5/6] Clean up --- tutorials/meps_tutorial.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/meps_tutorial.py b/tutorials/meps_tutorial.py index 1f653f1..430fbf5 100644 --- a/tutorials/meps_tutorial.py +++ b/tutorials/meps_tutorial.py @@ -387,6 +387,6 @@ def describe_metrics(metrics, thresh_arr): index.get_group_ratios("White", "Asian") # %% -index.get_group_ratios("White", "Non-White") +index.get_all_groups_failed_steps() # %% From c33400bf728f734994c51fd0448cde19fa4a2211 Mon Sep 17 00:00:00 2001 From: Harry Wilde Date: Fri, 25 Jul 2025 23:53:10 +0100 Subject: [PATCH 6/6] Clean up --- tutorials/meps_tutorial.py | 392 ------------------------------------- 1 file changed, 392 deletions(-) delete mode 100644 tutorials/meps_tutorial.py diff --git a/tutorials/meps_tutorial.py b/tutorials/meps_tutorial.py deleted file mode 100644 index 430fbf5..0000000 --- a/tutorials/meps_tutorial.py +++ /dev/null @@ -1,392 +0,0 @@ -# --- -# jupyter: -# jupytext: -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.17.2 -# kernelspec: -# display_name: .venv -# language: python -# name: python3 -# --- - -# %% [markdown] -# We will use the [`aif360`](https://aif360.readthedocs.io/en/latest/Getting%20Started.html) package to load the UCI adult dataset, fit a simple model and then analyse the fairness of the model using the DA-AUC. - -# %% -import numpy as np -from aif360.datasets import MEPSDataset19 -from aif360.explainers import MetricTextExplainer -from aif360.metrics import BinaryLabelDatasetMetric, ClassificationMetric -from sklearn.linear_model import LogisticRegression -from sklearn.pipeline import make_pipeline -from sklearn.preprocessing import StandardScaler - -np.random.seed(1) - -# %% [markdown] -# The below code downloads the required files from the UCI website if they are not already present in ai360's data directory. - -# %% -import os -import shutil -import subprocess - -import aif360 - -aif360_location = os.path.dirname(aif360.__file__) -meps_data_dir = os.path.join(aif360_location, "data", "raw", "meps") -h181_file_path = os.path.join(meps_data_dir, "h181.csv") - -if not os.path.isfile(h181_file_path): - r_script_path = os.path.join(meps_data_dir, "generate_data.R") - process = subprocess.Popen(["Rscript", r_script_path], stdin=subprocess.PIPE) - process.communicate(input=b"y\n") - - # Move the generated CSV files to meps_data_dir - generated_files = ["h181.csv", "h192.csv"] - for file_name in generated_files: - src_path = os.path.join(os.getcwd(), file_name) - dest_path = os.path.join(meps_data_dir, file_name) - if os.path.isfile(src_path): - shutil.move(src_path, dest_path) - - -# %% -def preprocessing_w_multimorb(df): - """ - 1.Create a new column, RACE that is 'White' if RACEV2X = 1 and HISPANX = 2 i.e. non Hispanic White - and 'non-White' otherwise - 2. Restrict to Panel 19 - 3. RENAME all columns that are PANEL/ROUND SPECIFIC - 4. Drop rows based on certain values of individual features that correspond to missing/unknown - generally < -1 - 5. Compute UTILIZATION, binarize it to 0 (< 10) and 1 (>= 10) - """ - - def race(row): - if (row["HISPANX"] == 2) and ( - row["RACEV2X"] == 1 - ): # non-Hispanic Whites are marked as WHITE; all others as NON-WHITE - return "White" - return "Non-White" - - df["RACE"] = df.apply(lambda row: race(row), axis=1) - - df = df[df["PANEL"] == 19] - - # RENAME COLUMNS - df = df.rename( - columns={ - "FTSTU53X": "FTSTU", - "ACTDTY53": "ACTDTY", - "HONRDC53": "HONRDC", - "RTHLTH53": "RTHLTH", - "MNHLTH53": "MNHLTH", - "CHBRON53": "CHBRON", - "JTPAIN53": "JTPAIN", - "PREGNT53": "PREGNT", - "WLKLIM53": "WLKLIM", - "ACTLIM53": "ACTLIM", - "SOCLIM53": "SOCLIM", - "COGLIM53": "COGLIM", - "EMPST53": "EMPST", - "REGION53": "REGION", - "MARRY53X": "MARRY", - "AGE53X": "AGE", - "POVCAT15": "POVCAT", - "INSCOV15": "INSCOV", - } - ) - - df = df[df["REGION"] >= 0] # remove values -1 - df = df[df["AGE"] >= 0] # remove values -1 - - df = df[df["MARRY"] >= 0] # remove values -1, -7, -8, -9 - - df = df[df["ASTHDX"] >= 0] # remove values -1, -7, -8, -9 - - df = df[ - ( - df[ - [ - "FTSTU", - "ACTDTY", - "HONRDC", - "RTHLTH", - "MNHLTH", - "HIBPDX", - "CHDDX", - "ANGIDX", - "EDUCYR", - "HIDEG", - "MIDX", - "OHRTDX", - "STRKDX", - "EMPHDX", - "CHBRON", - "CHOLDX", - "CANCERDX", - "DIABDX", - "JTPAIN", - "ARTHDX", - "ARTHTYPE", - "ASTHDX", - "ADHDADDX", - "PREGNT", - "WLKLIM", - "ACTLIM", - "SOCLIM", - "COGLIM", - "DFHEAR42", - "DFSEE42", - "ADSMOK42", - "PHQ242", - "EMPST", - "POVCAT", - "INSCOV", - ] - ] - >= -1 - ).all(1) - ] # for all other categorical features, remove values < -1 - - def utilization(row): - return row["OBTOTV15"] + row["OPTOTV15"] + row["ERTOT15"] + row["IPNGTD15"] + row["HHTOTD15"] - - df["TOTEXP15"] = df.apply(lambda row: utilization(row), axis=1) - lessE = df["TOTEXP15"] < 10.0 - df.loc[lessE, "TOTEXP15"] = 0.0 - moreE = df["TOTEXP15"] >= 10.0 - df.loc[moreE, "TOTEXP15"] = 1.0 - df["MULTIMORBIDITY"] = ( - df.filter(regex="DX$|CHBRON$|JTPAIN$").drop(columns=["ADHDADDX"]).apply(lambda x: (x == 1).sum(), axis=1) - ) - - df = df.rename(columns={"TOTEXP15": "UTILIZATION"}) - return df - - -# %% -dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test = MEPSDataset19( - custom_preprocessing=preprocessing_w_multimorb, - features_to_keep=[ - "REGION", - "AGE", - "SEX", - "RACE", - "RACEV2X", - "MARRY", - "FTSTU", - "ACTDTY", - "HONRDC", - "RTHLTH", - "MNHLTH", - "HIBPDX", - "CHDDX", - "ANGIDX", - "MIDX", - "OHRTDX", - "STRKDX", - "EMPHDX", - "CHBRON", - "CHOLDX", - "CANCERDX", - "DIABDX", - "JTPAIN", - "ARTHDX", - "ARTHTYPE", - "ASTHDX", - "ADHDADDX", - "PREGNT", - "WLKLIM", - "ACTLIM", - "SOCLIM", - "COGLIM", - "DFHEAR42", - "DFSEE42", - "ADSMOK42", - "PCS42", - "MCS42", - "K6SUM42", - "PHQ242", - "EMPST", - "POVCAT", - "INSCOV", - "UTILIZATION", - "PERWT15F", - "MULTIMORBIDITY", - ], -).split([0.5, 0.8], shuffle=True) - -sens_ind = 0 -sens_attr = dataset_orig_panel19_train.protected_attribute_names[sens_ind] - -unprivileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.unprivileged_protected_attributes[sens_ind]] -privileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.privileged_protected_attributes[sens_ind]] - - -# %% -def describe(train=None, val=None, test=None): - if train is not None: - print(train.features.shape) - if val is not None: - print(val.features.shape) - print(test.features.shape) - print(test.favorable_label, test.unfavorable_label) - print(test.protected_attribute_names) - print(test.privileged_protected_attributes, test.unprivileged_protected_attributes) - print(test.feature_names) - - -# %% -describe(dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test) - -# %% -metric_orig_panel19_train = BinaryLabelDatasetMetric( - dataset_orig_panel19_train, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups -) -explainer_orig_panel19_train = MetricTextExplainer(metric_orig_panel19_train) - -print(explainer_orig_panel19_train.disparate_impact()) - -# %% -dataset = dataset_orig_panel19_train -model = make_pipeline(StandardScaler(), LogisticRegression(solver="liblinear", random_state=1)) -fit_params = {"logisticregression__sample_weight": dataset.instance_weights} - -lr_orig_panel19 = model.fit(dataset.features, dataset.labels.ravel(), **fit_params) - -# %% -from collections import defaultdict - - -def test(dataset, model, thresh_arr): - try: - # sklearn classifier - y_val_pred_prob = model.predict_proba(dataset.features) - pos_ind = np.where(model.classes_ == dataset.favorable_label)[0][0] - except AttributeError: - # aif360 inprocessing algorithm - y_val_pred_prob = model.predict(dataset).scores - pos_ind = 0 - - metric_arrs = defaultdict(list) - for thresh in thresh_arr: - y_val_pred = (y_val_pred_prob[:, pos_ind] > thresh).astype(np.float64) - - dataset_pred = dataset.copy() - dataset_pred.labels = y_val_pred - metric = ClassificationMetric( - dataset, dataset_pred, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups - ) - - metric_arrs["bal_acc"].append((metric.true_positive_rate() + metric.true_negative_rate()) / 2) - metric_arrs["avg_odds_diff"].append(metric.average_odds_difference()) - metric_arrs["disp_imp"].append(metric.disparate_impact()) - metric_arrs["stat_par_diff"].append(metric.statistical_parity_difference()) - metric_arrs["eq_opp_diff"].append(metric.equal_opportunity_difference()) - metric_arrs["theil_ind"].append(metric.theil_index()) - - return metric_arrs - - -# %% -thresh_arr = np.linspace(0.01, 0.5, 50) -val_metrics = test(dataset=dataset_orig_panel19_val, model=lr_orig_panel19, thresh_arr=thresh_arr) -lr_orig_best_ind = np.argmax(val_metrics["bal_acc"]) - - -# %% -def describe_metrics(metrics, thresh_arr): - best_ind = np.argmax(metrics["bal_acc"]) - print("Threshold corresponding to Best balanced accuracy: {:6.4f}".format(thresh_arr[best_ind])) - print("Best balanced accuracy: {:6.4f}".format(metrics["bal_acc"][best_ind])) - disp_imp_at_best_ind = 1 - min(metrics["disp_imp"][best_ind], 1 / metrics["disp_imp"][best_ind]) - print("Corresponding 1-min(DI, 1/DI) value: {:6.4f}".format(disp_imp_at_best_ind)) - print("Corresponding average odds difference value: {:6.4f}".format(metrics["avg_odds_diff"][best_ind])) - print("Corresponding statistical parity difference value: {:6.4f}".format(metrics["stat_par_diff"][best_ind])) - print("Corresponding equal opportunity difference value: {:6.4f}".format(metrics["eq_opp_diff"][best_ind])) - print("Corresponding Theil index value: {:6.4f}".format(metrics["theil_ind"][best_ind])) - - -# %% -describe_metrics(val_metrics, thresh_arr) - -# %% -lr_orig_metrics = test( - dataset=dataset_orig_panel19_test, model=lr_orig_panel19, thresh_arr=[thresh_arr[lr_orig_best_ind]] -) - -# %% -describe_metrics(lr_orig_metrics, [thresh_arr[lr_orig_best_ind]]) - -# %% -df = dataset.convert_to_dataframe()[0] - -# %% -df["MULTIMORBIDITY"].hist() -df_add = df.copy() -df_add["RTHLTH"] = df_add.filter(regex="RTHLTH").idxmax(axis=1) - -# %% -df_add.boxplot(column="MULTIMORBIDITY", by="RTHLTH") - -# %% -from daindex import DAIndex, DeteriorationFeature, Group - -# %% [markdown] -# Let's first specify our deterioration feature and groups. - -# %% -det_feature = DeteriorationFeature(col="MULTIMORBIDITY", threshold=1, is_discrete=True) -groups = [ - Group("White", 1, "RACEV2X"), - Group( - "Non-White", [12, 5, 10, 2, 3, 6, 4], "RACEV2X", det_threshold=1.2 - ), # Modify deterioration threshold only for this group - Group("Black", 2, "RACEV2X"), - Group("Asian", [4, 5, 6, 10], "RACEV2X"), -] -print(det_feature) -print(groups[1]) - -# %% [markdown] -# We can call a group on the cohort to get the group representation in the cohort: - -# %% -print(groups[0](df).RACEV2X.value_counts()) -groups[0](df).head() - -# %% [markdown] -# We can now instantiate the `DAIndex` with the components above and evaluate all groups via the model we have made. Note this could also be done with a list of arbitrary models. The package also supports evaluation via an existing predictions column in the cohort, see `evaluate_all_groups_from_predictions`. We could also just compare two specific groups with `evaluate_group_pair_by_models` and `evaluate_group_pair_from_predictions`. - -# %% -feature_list = dataset.feature_names -index = DAIndex( - cohort=df, - groups=groups, - det_feature=det_feature, - decision_boundary=0.75, # Set decision boundary for DA curve -) -index.evaluate_all_groups_by_models( - models=lr_orig_panel19, - feature_list=feature_list, - reference_group="White", - n_jobs=1, # Up n_jobs to parallelize if needed -) - -# %% -index.present_all_results() - -# %% -index.get_all_ratios() - -# %% -index.get_group_ratios("White", "Asian") - -# %% -index.get_all_groups_failed_steps() - -# %%