tutorials/meps_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from aif360.datasets import MEPSDataset19\n",
    "from aif360.explainers import MetricTextExplainer\n",
    "from aif360.metrics import BinaryLabelDatasetMetric, ClassificationMetric\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below code downloads the required files from the UCI website if they are not already present in ai360's data directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import subprocess\n",
    "\n",
    "import aif360\n",
    "\n",
    "aif360_location = os.path.dirname(aif360.__file__)\n",
    "meps_data_dir = os.path.join(aif360_location, \"data\", \"raw\", \"meps\")\n",
    "h181_file_path = os.path.join(meps_data_dir, \"h181.csv\")\n",
    "\n",
    "if not os.path.isfile(h181_file_path):\n",
    "    r_script_path = os.path.join(meps_data_dir, \"generate_data.R\")\n",
    "    process = subprocess.Popen([\"Rscript\", r_script_path], stdin=subprocess.PIPE)\n",
    "    process.communicate(input=b\"y\\n\")\n",
    "\n",
    "    # Move the generated CSV files to meps_data_dir\n",
    "    generated_files = [\"h181.csv\", \"h192.csv\"]\n",
    "    for file_name in generated_files:\n",
    "        src_path = os.path.join(os.getcwd(), file_name)\n",
    "        dest_path = os.path.join(meps_data_dir, file_name)\n",
    "        if os.path.isfile(src_path):\n",
    "            shutil.move(src_path, dest_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocessing_w_multimorb(df):\n",
    "    \"\"\"\n",
    "    1.Create a new column, RACE that is 'White' if RACEV2X = 1 and HISPANX = 2 i.e. non Hispanic White\n",
    "      and 'non-White' otherwise\n",
    "    2. Restrict to Panel 19\n",
    "    3. RENAME all columns that are PANEL/ROUND SPECIFIC\n",
    "    4. Drop rows based on certain values of individual features that correspond to missing/unknown - generally < -1\n",
    "    5. Compute UTILIZATION, binarize it to 0 (< 10) and 1 (>= 10)\n",
    "    \"\"\"\n",
    "\n",
    "    def race(row):\n",
    "        if (row[\"HISPANX\"] == 2) and (\n",
    "            row[\"RACEV2X\"] == 1\n",
    "        ):  # non-Hispanic Whites are marked as WHITE; all others as NON-WHITE\n",
    "            return \"White\"\n",
    "        return \"Non-White\"\n",
    "\n",
    "    df[\"RACE\"] = df.apply(lambda row: race(row), axis=1)\n",
    "\n",
    "    df = df[df[\"PANEL\"] == 19]\n",
    "\n",
    "    # RENAME COLUMNS\n",
    "    df = df.rename(\n",
    "        columns={\n",
    "            \"FTSTU53X\": \"FTSTU\",\n",
    "            \"ACTDTY53\": \"ACTDTY\",\n",
    "            \"HONRDC53\": \"HONRDC\",\n",
    "            \"RTHLTH53\": \"RTHLTH\",\n",
    "            \"MNHLTH53\": \"MNHLTH\",\n",
    "            \"CHBRON53\": \"CHBRON\",\n",
    "            \"JTPAIN53\": \"JTPAIN\",\n",
    "            \"PREGNT53\": \"PREGNT\",\n",
    "            \"WLKLIM53\": \"WLKLIM\",\n",
    "            \"ACTLIM53\": \"ACTLIM\",\n",
    "            \"SOCLIM53\": \"SOCLIM\",\n",
    "            \"COGLIM53\": \"COGLIM\",\n",
    "            \"EMPST53\": \"EMPST\",\n",
    "            \"REGION53\": \"REGION\",\n",
    "            \"MARRY53X\": \"MARRY\",\n",
    "            \"AGE53X\": \"AGE\",\n",
    "            \"POVCAT15\": \"POVCAT\",\n",
    "            \"INSCOV15\": \"INSCOV\",\n",
    "        }\n",
    "    )\n",
    "\n",
    "    df = df[df[\"REGION\"] >= 0]  # remove values -1\n",
    "    df = df[df[\"AGE\"] >= 0]  # remove values -1\n",
    "\n",
    "    df = df[df[\"MARRY\"] >= 0]  # remove values -1, -7, -8, -9\n",
    "\n",
    "    df = df[df[\"ASTHDX\"] >= 0]  # remove values -1, -7, -8, -9\n",
    "\n",
    "    df = df[\n",
    "        (\n",
    "            df[\n",
    "                [\n",
    "                    \"FTSTU\",\n",
    "                    \"ACTDTY\",\n",
    "                    \"HONRDC\",\n",
    "                    \"RTHLTH\",\n",
    "                    \"MNHLTH\",\n",
    "                    \"HIBPDX\",\n",
    "                    \"CHDDX\",\n",
    "                    \"ANGIDX\",\n",
    "                    \"EDUCYR\",\n",
    "                    \"HIDEG\",\n",
    "                    \"MIDX\",\n",
    "                    \"OHRTDX\",\n",
    "                    \"STRKDX\",\n",
    "                    \"EMPHDX\",\n",
    "                    \"CHBRON\",\n",
    "                    \"CHOLDX\",\n",
    "                    \"CANCERDX\",\n",
    "                    \"DIABDX\",\n",
    "                    \"JTPAIN\",\n",
    "                    \"ARTHDX\",\n",
    "                    \"ARTHTYPE\",\n",
    "                    \"ASTHDX\",\n",
    "                    \"ADHDADDX\",\n",
    "                    \"PREGNT\",\n",
    "                    \"WLKLIM\",\n",
    "                    \"ACTLIM\",\n",
    "                    \"SOCLIM\",\n",
    "                    \"COGLIM\",\n",
    "                    \"DFHEAR42\",\n",
    "                    \"DFSEE42\",\n",
    "                    \"ADSMOK42\",\n",
    "                    \"PHQ242\",\n",
    "                    \"EMPST\",\n",
    "                    \"POVCAT\",\n",
    "                    \"INSCOV\",\n",
    "                ]\n",
    "            ]\n",
    "            >= -1\n",
    "        ).all(1)\n",
    "    ]  # for all other categorical features, remove values < -1\n",
    "\n",
    "    def utilization(row):\n",
    "        return row[\"OBTOTV15\"] + row[\"OPTOTV15\"] + row[\"ERTOT15\"] + row[\"IPNGTD15\"] + row[\"HHTOTD15\"]\n",
    "\n",
    "    df[\"TOTEXP15\"] = df.apply(lambda row: utilization(row), axis=1)\n",
    "    lessE = df[\"TOTEXP15\"] < 10.0\n",
    "    df.loc[lessE, \"TOTEXP15\"] = 0.0\n",
    "    moreE = df[\"TOTEXP15\"] >= 10.0\n",
    "    df.loc[moreE, \"TOTEXP15\"] = 1.0\n",
    "    df[\"MULTIMORBIDITY\"] = (\n",
    "        df.filter(regex=\"DX$|CHBRON$|JTPAIN$\").drop(columns=[\"ADHDADDX\"]).apply(lambda x: (x == 1).sum(), axis=1)\n",
    "    )\n",
    "\n",
    "    df = df.rename(columns={\"TOTEXP15\": \"UTILIZATION\"})\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test = MEPSDataset19(\n",
    "    custom_preprocessing=preprocessing_w_multimorb,\n",
    "    features_to_keep=[\n",
    "        \"REGION\",\n",
    "        \"AGE\",\n",
    "        \"SEX\",\n",
    "        \"RACE\",\n",
    "        \"RACEV2X\",\n",
    "        \"MARRY\",\n",
    "        \"FTSTU\",\n",
    "        \"ACTDTY\",\n",
    "        \"HONRDC\",\n",
    "        \"RTHLTH\",\n",
    "        \"MNHLTH\",\n",
    "        \"HIBPDX\",\n",
    "        \"CHDDX\",\n",
    "        \"ANGIDX\",\n",
    "        \"MIDX\",\n",
    "        \"OHRTDX\",\n",
    "        \"STRKDX\",\n",
    "        \"EMPHDX\",\n",
    "        \"CHBRON\",\n",
    "        \"CHOLDX\",\n",
    "        \"CANCERDX\",\n",
    "        \"DIABDX\",\n",
    "        \"JTPAIN\",\n",
    "        \"ARTHDX\",\n",
    "        \"ARTHTYPE\",\n",
    "        \"ASTHDX\",\n",
    "        \"ADHDADDX\",\n",
    "        \"PREGNT\",\n",
    "        \"WLKLIM\",\n",
    "        \"ACTLIM\",\n",
    "        \"SOCLIM\",\n",
    "        \"COGLIM\",\n",
    "        \"DFHEAR42\",\n",
    "        \"DFSEE42\",\n",
    "        \"ADSMOK42\",\n",
    "        \"PCS42\",\n",
    "        \"MCS42\",\n",
    "        \"K6SUM42\",\n",
    "        \"PHQ242\",\n",
    "        \"EMPST\",\n",
    "        \"POVCAT\",\n",
    "        \"INSCOV\",\n",
    "        \"UTILIZATION\",\n",
    "        \"PERWT15F\",\n",
    "        \"MULTIMORBIDITY\",\n",
    "    ],\n",
    ").split([0.5, 0.8], shuffle=True)\n",
    "\n",
    "sens_ind = 0\n",
    "sens_attr = dataset_orig_panel19_train.protected_attribute_names[sens_ind]\n",
    "\n",
    "unprivileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.unprivileged_protected_attributes[sens_ind]]\n",
    "privileged_groups = [{sens_attr: v} for v in dataset_orig_panel19_train.privileged_protected_attributes[sens_ind]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def describe(train=None, val=None, test=None):\n",
    "    if train is not None:\n",
    "        print(train.features.shape)\n",
    "    if val is not None:\n",
    "        print(val.features.shape)\n",
    "    print(test.features.shape)\n",
    "    print(test.favorable_label, test.unfavorable_label)\n",
    "    print(test.protected_attribute_names)\n",
    "    print(test.privileged_protected_attributes, test.unprivileged_protected_attributes)\n",
    "    print(test.feature_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "describe(dataset_orig_panel19_train, dataset_orig_panel19_val, dataset_orig_panel19_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metric_orig_panel19_train = BinaryLabelDatasetMetric(\n",
    "    dataset_orig_panel19_train, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups\n",
    ")\n",
    "explainer_orig_panel19_train = MetricTextExplainer(metric_orig_panel19_train)\n",
    "\n",
    "print(explainer_orig_panel19_train.disparate_impact())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = dataset_orig_panel19_train\n",
    "model = make_pipeline(StandardScaler(), LogisticRegression(solver=\"liblinear\", random_state=1))\n",
    "fit_params = {\"logisticregression__sample_weight\": dataset.instance_weights}\n",
    "\n",
    "lr_orig_panel19 = model.fit(dataset.features, dataset.labels.ravel(), **fit_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "def test(dataset, model, thresh_arr):\n",
    "    try:\n",
    "        # sklearn classifier\n",
    "        y_val_pred_prob = model.predict_proba(dataset.features)\n",
    "        pos_ind = np.where(model.classes_ == dataset.favorable_label)[0][0]\n",
    "    except AttributeError:\n",
    "        # aif360 inprocessing algorithm\n",
    "        y_val_pred_prob = model.predict(dataset).scores\n",
    "        pos_ind = 0\n",
    "\n",
    "    metric_arrs = defaultdict(list)\n",
    "    for thresh in thresh_arr:\n",
    "        y_val_pred = (y_val_pred_prob[:, pos_ind] > thresh).astype(np.float64)\n",
    "\n",
    "        dataset_pred = dataset.copy()\n",
    "        dataset_pred.labels = y_val_pred\n",
    "        metric = ClassificationMetric(\n",
    "            dataset, dataset_pred, unprivileged_groups=unprivileged_groups, privileged_groups=privileged_groups\n",
    "        )\n",
    "\n",
    "        metric_arrs[\"bal_acc\"].append((metric.true_positive_rate() + metric.true_negative_rate()) / 2)\n",
    "        metric_arrs[\"avg_odds_diff\"].append(metric.average_odds_difference())\n",
    "        metric_arrs[\"disp_imp\"].append(metric.disparate_impact())\n",
    "        metric_arrs[\"stat_par_diff\"].append(metric.statistical_parity_difference())\n",
    "        metric_arrs[\"eq_opp_diff\"].append(metric.equal_opportunity_difference())\n",
    "        metric_arrs[\"theil_ind\"].append(metric.theil_index())\n",
    "\n",
    "    return metric_arrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh_arr = np.linspace(0.01, 0.5, 50)\n",
    "val_metrics = test(dataset=dataset_orig_panel19_val, model=lr_orig_panel19, thresh_arr=thresh_arr)\n",
    "lr_orig_best_ind = np.argmax(val_metrics[\"bal_acc\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def describe_metrics(metrics, thresh_arr):\n",
    "    best_ind = np.argmax(metrics[\"bal_acc\"])\n",
    "    print(\"Threshold corresponding to Best balanced accuracy: {:6.4f}\".format(thresh_arr[best_ind]))\n",
    "    print(\"Best balanced accuracy: {:6.4f}\".format(metrics[\"bal_acc\"][best_ind]))\n",
    "    disp_imp_at_best_ind = 1 - min(metrics[\"disp_imp\"][best_ind], 1 / metrics[\"disp_imp\"][best_ind])\n",
    "    print(\"Corresponding 1-min(DI, 1/DI) value: {:6.4f}\".format(disp_imp_at_best_ind))\n",
    "    print(\"Corresponding average odds difference value: {:6.4f}\".format(metrics[\"avg_odds_diff\"][best_ind]))\n",
    "    print(\"Corresponding statistical parity difference value: {:6.4f}\".format(metrics[\"stat_par_diff\"][best_ind]))\n",
    "    print(\"Corresponding equal opportunity difference value: {:6.4f}\".format(metrics[\"eq_opp_diff\"][best_ind]))\n",
    "    print(\"Corresponding Theil index value: {:6.4f}\".format(metrics[\"theil_ind\"][best_ind]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "describe_metrics(val_metrics, thresh_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_orig_metrics = test(\n",
    "    dataset=dataset_orig_panel19_test, model=lr_orig_panel19, thresh_arr=[thresh_arr[lr_orig_best_ind]]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "describe_metrics(lr_orig_metrics, [thresh_arr[lr_orig_best_ind]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dataset.convert_to_dataframe()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"MULTIMORBIDITY\"].hist()\n",
    "df_add = df.copy()\n",
    "df_add[\"RTHLTH\"] = df_add.filter(regex=\"RTHLTH\").idxmax(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_add.boxplot(column=\"MULTIMORBIDITY\", by=\"RTHLTH\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from daindex import DAIndex, DeteriorationFeature, Group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first specify our deterioration feature and groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "det_feature = DeteriorationFeature(col=\"MULTIMORBIDITY\", threshold=1, is_discrete=True)\n",
    "groups = [\n",
    "    Group(\"White\", 1, \"RACEV2X\"),\n",
    "    Group(\n",
    "        \"Non-White\", [12, 5, 10, 2, 3, 6, 4], \"RACEV2X\", det_threshold=1.2\n",
    "    ),  # Modify deterioration threshold only for this group\n",
    "    Group(\"Black\", 2, \"RACEV2X\"),\n",
    "    Group(\"Asian\", [4, 5, 6, 10], \"RACEV2X\"),\n",
    "]\n",
    "print(det_feature)\n",
    "print(groups[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can call a group on the cohort to get the group representation in the cohort:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(groups[0](df).RACEV2X.value_counts())\n",
    "groups[0](df).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_list = dataset.feature_names\n",
    "index = DAIndex(\n",
    "    cohort=df,\n",
    "    groups=groups,\n",
    "    det_feature=det_feature,\n",
    ")\n",
    "index.evaluate_all_groups_by_models(\n",
    "    models=lr_orig_panel19, feature_list=feature_list, reference_group=\"White\", n_jobs=6\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index.present_all_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index.get_all_ratios()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "daindex",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}