explanations/tsv2_train_outlier_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.stats import mannwhitneyu, spearmanr, ttest_ind\n",
    "from sklearn.ensemble import IsolationForest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_agg_explanations = (\n",
    "    pd.read_csv(\"data/tsv2_train_aggregated_sg_explanations_with_dice.csv\")\n",
    "    .dropna()\n",
    "    .reset_index()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ID2LABELS = {\n",
    "    0: \"background\",\n",
    "    1: \"aorta\",\n",
    "    2: \"lung_upper_lobe_left\",\n",
    "    3: \"lung_lower_lobe_left\",\n",
    "    4: \"lung_upper_lobe_right\",\n",
    "    5: \"lung_middle_lobe_right\",\n",
    "    6: \"lung_lower_lobe_right\",\n",
    "    7: \"trachea\",\n",
    "    8: \"heart\",\n",
    "    9: \"pulmonary_vein\",\n",
    "    10: \"thyroid_gland\",\n",
    "    11: \"ribs\",\n",
    "    12: \"vertebraes\",\n",
    "    13: \"autochthon_left\",\n",
    "    14: \"autochthon_right\",\n",
    "    15: \"sternum\",\n",
    "    16: \"costal_cartilages\",\n",
    "}\n",
    "\n",
    "labels = list(ID2LABELS.values())\n",
    "\n",
    "MASS_COLUMNS_FORMAT = \"{label_1}_explanation_in_{label_2}_mass\"\n",
    "\n",
    "outliers_analysis_datasets = {}\n",
    "columns = {}\n",
    "for label in labels:\n",
    "    mass_columns = [\n",
    "        MASS_COLUMNS_FORMAT.format(label_1=label, label_2=label_2) for label_2 in labels\n",
    "    ]\n",
    "    columns[label] = mass_columns\n",
    "    get_only_first_record = False\n",
    "\n",
    "    # drop explanations with NA\n",
    "    aggregated_explanations = train_agg_explanations.dropna().reset_index(drop=True)\n",
    "\n",
    "    mass_explanations = aggregated_explanations[mass_columns].copy()\n",
    "    outliers_analysis_datasets[label] = mass_explanations\n",
    "\n",
    "label_mass_in_label_columns = [\n",
    "    f\"{label}_explanation_in_{label}_mass\" for label in labels\n",
    "]\n",
    "iou_columns = [f\"{label}_iou\" for label in labels]\n",
    "dice_columns = [f\"{label}_dice\" for label in labels]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_context = (\n",
    "    aggregated_explanations[label_mass_in_label_columns]\n",
    "    .describe()\n",
    "    .T.sort_values(by=\"mean\", ascending=False)\n",
    "    .drop(columns=[\"count\", \"std\"])\n",
    ")\n",
    "train_context.index = [\n",
    "    idx.split(\"_explanation\")[0].replace(\"_\", \" \").title()\n",
    "    for idx in train_context.index\n",
    "]\n",
    "train_context = train_context.reset_index().rename(columns={\"index\": \"Label\"})\n",
    "print(train_context.to_markdown(index=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass_ifs = {label: IsolationForest(random_state=123) for label in labels}\n",
    "for label, mass_if in mass_ifs.items():\n",
    "    mass_if.fit(outliers_analysis_datasets[label])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_agg_explanations = pd.read_csv(\n",
    "    \"data/tsv2_test_aggregated_sg_explanations_with_dice.csv\"\n",
    ")\n",
    "\n",
    "b50_agg_explanations = (\n",
    "    pd.read_csv(\"data/b50_aggregated_sg_explanations.csv\")\n",
    "    .dropna()\n",
    "    .reset_index(drop=True)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsv2_test_dataset = {label: test_agg_explanations[columns[label]] for label in labels}\n",
    "\n",
    "b50_datasets = {label: b50_agg_explanations[columns[label]] for label in labels}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How contextual are tsv2_test and b50\n",
    "\n",
    "tsv2_test_context = (\n",
    "    test_agg_explanations[label_mass_in_label_columns]\n",
    "    .describe()\n",
    "    .T.sort_values(by=\"mean\", ascending=False)\n",
    "    .drop(columns=[\"count\", \"std\"])\n",
    ")\n",
    "tsv2_test_context.index = [\n",
    "    idx.split(\"_explanation\")[0].replace(\"_\", \" \").title()\n",
    "    for idx in tsv2_test_context.index\n",
    "]\n",
    "tsv2_test_context = tsv2_test_context.reset_index().rename(columns={\"index\": \"Label\"})\n",
    "print(tsv2_test_context.to_markdown(index=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_context = (\n",
    "    b50_agg_explanations[label_mass_in_label_columns]\n",
    "    .describe()\n",
    "    .T.sort_values(by=\"mean\", ascending=False)\n",
    "    .drop(columns=[\"count\", \"std\"])\n",
    ")\n",
    "b50_context.index = [\n",
    "    idx.split(\"_explanation\")[0].replace(\"_\", \" \").title() for idx in b50_context.index\n",
    "]\n",
    "b50_context = b50_context.reset_index().rename(columns={\"index\": \"Label\"})\n",
    "print(b50_context.to_markdown(index=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_outliers_per_label = {\n",
    "    label: (\n",
    "        mass_if.predict(b50_datasets[label]),\n",
    "        mass_if.decision_function(b50_datasets[label]),\n",
    "    )\n",
    "    for label, mass_if in mass_ifs.items()\n",
    "}\n",
    "\n",
    "tsv2_test_outliers_per_label = {\n",
    "    label: (\n",
    "        mass_if.predict(tsv2_test_dataset[label]),\n",
    "        mass_if.decision_function(tsv2_test_dataset[label]),\n",
    "    )\n",
    "    for label, mass_if in mass_ifs.items()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsv2_test_outliers_per_label[\"heart\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_dice = test_agg_explanations[dice_columns].copy()\n",
    "test_iou = test_agg_explanations[iou_columns].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(\n",
    "    [outlier_scores for _, outlier_scores in tsv2_test_outliers_per_label.values()],\n",
    "    axis=0,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data = []\n",
    "\n",
    "for label, (_, outlier_scores) in tsv2_test_outliers_per_label.items():\n",
    "    test_score = spearmanr(test_dice[f\"{label}_dice\"], outlier_scores)\n",
    "    test_data.append(\n",
    "        {\n",
    "            \"label\": label,\n",
    "            \"p-value\": test_score.pvalue,\n",
    "            \"spearmanr\": test_score.correlation,\n",
    "        }\n",
    "    )\n",
    "    print(f\"p-val for {label}: {test_score.pvalue:,.3g}\")\n",
    "mean_outlier_score = np.mean(\n",
    "    [outlier_scores for _, outlier_scores in tsv2_test_outliers_per_label.values()],\n",
    "    axis=0,\n",
    ")\n",
    "mean_dice = test_dice.mean(axis=1)\n",
    "test_score = spearmanr(mean_dice, mean_outlier_score)\n",
    "print(\n",
    "    f\"p-val for mean outlier score: {test_score.pvalue:,.3g}, spearmanr: {test_score.correlation}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data = pd.DataFrame(test_data)\n",
    "print_data = test_data.copy()\n",
    "print_data.label = print_data.label.apply(lambda x: x.replace(\"_\", \" \"))\n",
    "print_data = print_data.set_index(\"label\")\n",
    "print_data[\"p-value\"] = print_data[\"p-value\"].map(\"{:,.3g}\".format)\n",
    "print_data[\"spearmanr\"] = print_data[\"spearmanr\"].map(\"{:,.3g}\".format)\n",
    "print(print_data.sort_values(\"p-value\").to_latex())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_theme()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(x=mean_outlier_score, y=mean_dice)\n",
    "plt.xlabel(\"Mean Outlier Score\")\n",
    "plt.ylabel(\"Mean Dice\")\n",
    "plt.title(\"Mean Outlier Score vs Mean Dice\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(4, 4, figsize=(12, 12))\n",
    "\n",
    "for i, (ax, label) in enumerate(zip(axs.flatten(), labels[1:]), start=1):\n",
    "    sns.scatterplot(\n",
    "        x=tsv2_test_outliers_per_label[label][1], y=test_dice[f\"{label}_dice\"], ax=ax\n",
    "    )\n",
    "    label_for_print = label.replace(\"_\", \" \").title()\n",
    "    ax.set_title(label_for_print)\n",
    "    ax.set_xlabel(\"Outlier score\")\n",
    "    ax.set_ylabel(f\"{label_for_print} Dice\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get patients selected as outliers in b50\n",
    "b50_outliers = []\n",
    "for label, (outliers, _) in b50_outliers_per_label.items():\n",
    "    b50_outliers.append(\n",
    "        np.unique(\n",
    "            b50_agg_explanations.patient.apply(lambda x: int(x.split(\"/\")[0]))[\n",
    "                outliers == -1\n",
    "            ]\n",
    "        )\n",
    "    )\n",
    "\n",
    "b50_outliers = np.unique(np.concatenate(b50_outliers))\n",
    "b50_outliers, b50_outliers.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlier_score_lable = \"lung_lower_lobe_right\"\n",
    "b50_with_outlier_score = b50_agg_explanations.copy()\n",
    "# b50_with_outlier_score[\"outlier_score\"] = np.mean(\n",
    "#     [outlier_scores for _, outlier_scores in b50_outliers_per_label.values()], axis=0\n",
    "# )\n",
    "b50_with_outlier_score[\"outlier_score\"] = b50_outliers_per_label[outlier_score_lable][1]\n",
    "b50_with_outlier_score[\"patient_id\"] = b50_with_outlier_score.patient.apply(\n",
    "    lambda x: int(x.split(\"/\")[0])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_data = pd.read_csv(\"b50/patient_data/patient_metadata_1.0/csv_1.0.csv\")\n",
    "\n",
    "with open(\"b50/patient_data/descriptions_1.0/json_1.0.json\") as f:\n",
    "    patient_descriptions = json.load(f)\n",
    "\n",
    "patient_descriptions = pd.DataFrame(patient_descriptions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_outliers_metadata = meta_data[meta_data[\"patient_id\"].isin(b50_outliers)]\n",
    "b50_inliers_metadata = meta_data[~meta_data[\"patient_id\"].isin(b50_outliers)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney U test\n",
    "u_stat, u_p_val = mannwhitneyu(\n",
    "    b50_outliers_metadata[\"age\"], b50_inliers_metadata[\"age\"]\n",
    ")\n",
    "print(f\"Mann-Whitney U test P-value: {u_p_val}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_stat, t_p_val = ttest_ind(b50_outliers_metadata[\"age\"], b50_inliers_metadata[\"age\"])\n",
    "print(f\"t-test P-value: {t_p_val}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_outliers_metadata[\"sex\"].value_counts(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_inliers_metadata[\"sex\"].value_counts(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers_descriptions = patient_descriptions[\n",
    "    patient_descriptions[\"patient_id\"].isin(b50_outliers)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers_descriptions.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_theme(style=\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_3_b50_outliers = (\n",
    "    b50_with_outlier_score.drop_duplicates(subset=[\"patient_id\"])[\n",
    "        [\"patient\", \"outlier_score\"]\n",
    "    ]\n",
    "    .sort_values(by=\"outlier_score\")\n",
    "    .iloc[:3]\n",
    ")\n",
    "top_3_b50_inliers = (\n",
    "    b50_with_outlier_score.drop_duplicates(subset=[\"patient_id\"])[\n",
    "        [\"patient\", \"outlier_score\"]\n",
    "    ]\n",
    "    .sort_values(by=\"outlier_score\")\n",
    "    .iloc[-3:]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_3_b50_outliers, top_3_b50_inliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b50_slices_outliers = [[181, 153, 85], [93, 140, 143], [174, 150, 133]]\n",
    "b50_slices_inliers = [[92, 95, 113], [100, 98, 81], [81, 89, 147]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from generation_utils import B50_FOLDER, get_transform, read_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = get_transform()\n",
    "\n",
    "\n",
    "def load_image(path):\n",
    "    return transform(read_path(path)).squeeze().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers_imgs = [\n",
    "    load_image(str(B50_FOLDER / patient_id))\n",
    "    for patient_id in top_3_b50_outliers.patient\n",
    "]\n",
    "inliers_imgs = [\n",
    "    load_image(str(B50_FOLDER / patient_id)) for patient_id in top_3_b50_inliers.patient\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 5))\n",
    "subfigs = fig.subfigures(1, 2, wspace=0.1)\n",
    "\n",
    "subfigs[0].suptitle(\"Top 3 Anomaly Scores\", y=1.01)\n",
    "subfigs[1].suptitle(\"Bottom 3 Anomaly Scores\", y=1.01)\n",
    "ax = subfigs[0].subplots(3, 3)\n",
    "plt.subplots_adjust(hspace=0.01, wspace=0.01)\n",
    "for i, (img, slices) in enumerate(zip(outliers_imgs, b50_slices_outliers)):\n",
    "    ax[i, 0].imshow(np.rot90(img[slices[0]]), cmap=\"gray\")\n",
    "    ax[i, 0].axis(\"off\")\n",
    "    ax[i, 0].set_aspect(\"auto\")\n",
    "    ax[i, 1].imshow(np.rot90(img[:, slices[1]]), cmap=\"gray\")\n",
    "    ax[i, 1].axis(\"off\")\n",
    "    ax[i, 1].set_aspect(\"auto\")\n",
    "    ax[i, 2].imshow(np.rot90(img[:, :, slices[2]]), cmap=\"gray\")\n",
    "    ax[i, 2].axis(\"off\")\n",
    "    ax[i, 2].set_aspect(\"auto\")\n",
    "\n",
    "ax = subfigs[1].subplots(3, 3)\n",
    "\n",
    "for i, (img, slices) in enumerate(zip(inliers_imgs, b50_slices_inliers)):\n",
    "    ax[i, 0].imshow(np.rot90(img[slices[0]]), cmap=\"gray\")\n",
    "    ax[i, 0].axis(\"off\")\n",
    "    ax[i, 0].set_aspect(\"auto\")\n",
    "    ax[i, 1].imshow(np.rot90(img[:, slices[1]]), cmap=\"gray\")\n",
    "    ax[i, 1].axis(\"off\")\n",
    "    ax[i, 1].set_aspect(\"auto\")\n",
    "    ax[i, 2].imshow(np.rot90(img[:, :, slices[2]]), cmap=\"gray\")\n",
    "    ax[i, 2].axis(\"off\")\n",
    "    ax[i, 2].set_aspect(\"auto\")  # Adjust the layout of the subplots\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"top3_outliers_inliers.pdf\", bbox_inches=\"tight\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "xlungs",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}