{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Why do we care about finding data issues in the VinBigData Chest X-ray competition?\n\nIn the [VinBigData Chest X-ray Abnormalities Detection competition](https://www.kaggle.com/c/vinbigdata-chest-xray-abnormalities-detection/overview), we are given bounding boxes and their assignments to classes assigned independently by 3 radiologists as our training data. On the test data, there is an additional step, where 2 more experienced radiologists review what the 3 initial assessments were. These two more experienced ones then decide what to make of the disagreements (and can discuss with each other while doing so). The process is described on the [paper linked on the competition overview page](https://arxiv.org/pdf/2012.15029.pdf):\n> For the test set, 5 radiologists involved into a two-stage labeling process. During the first stage, each image was independently annotated by 3 radiologists. In the second stage, 2 other\nradiologists, who have a higher level of experience, reviewed the annotations of the 3 previous annotators and communicated with each other in order to decide the final labels. The disagreements among initial annotators were carefully discussed and resolved by the 2 reviewers. Finally, the consensus of their opinions will serve as reference ground-truth.\n\n(see also some further discussion on this in [this discussion thread](https://www.kaggle.com/c/vinbigdata-chest-xray-abnormalities-detection/discussion/211035), which also ends up discussing some data/labeling issues). \n\n**Bottom line:** our training data is not generated from the same process as the test data and a lot of the differences will be mistakes in the training data that would likely have been fixed, if the image were in the test data. Thus, it makes sense to try and fix the issues we can fix."},{"metadata":{},"cell_type":"markdown","source":"# What mistakes in the data do we already know?\nSome of the [forum discussions](https://www.kaggle.com/c/vinbigdata-chest-xray-abnormalities-detection/discussion/211035) and some [EDA notebooks](https://www.kaggle.com/bjoernholzhauer/eda-dicom-reading-vinbigdata-chest-x-ray#5.-How-big-do-bounding-boxes-tend-to-be-for-different-classes?-How-many-are-there?) already suggested that at least some 'aortic enlargement' and 'cardiomegaly' labels for bounding boxes are wrong. To me this looks like these were maybe less mis-assessments by the radiologists, but rather software usage issues, where they created a new bounding box, but accidentally gave it the same class as for the previous bounding box they had labelled. For these two classes this is relatively easy to see (see images in [data issues Section](https://www.kaggle.com/bjoernholzhauer/eda-dicom-reading-vinbigdata-chest-x-ray#11.Possible-data-issues) of [this EDA notebook](https://www.kaggle.com/bjoernholzhauer/eda-dicom-reading-vinbigdata-chest-x-ray#5.-How-big-do-bounding-boxes-tend-to-be-for-different-classes?-How-many-are-there?). However, it will also be interesting to see whether other approaches also identify these cases"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import os\nimport re\nimport pandas as pd\n#from fastai.medical.imaging import *\n#from fastai.vision.all import *\nimport numpy as np\nfrom pathlib import Path\nimport matplotlib\nimport matplotlib.patches as ptc\nfrom tqdm import tqdm # for getting a progress bar on loops\nfrom sklearn.model_selection import KFold\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport shap\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Rectangle\nimport shelve","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"dicom_meta = pd.read_csv('../input/eda-dicom-reading-vinbigdata-chest-x-ray/train_dicom_properties.csv.bz2').rename(columns={'file':'image_id'})\ntrain = pd.read_csv('../input/vinbigdata-chest-xray-abnormalities-detection/train.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"bygroup = train[['image_id', 'class_name', 'rad_id','class_id']].copy().groupby(['image_id', 'class_name', 'rad_id']).count().reset_index()\nstrange_cases = bygroup.loc[ [ (cn[0] in ['Aortic enlargement', 'Cardiomegaly']) & cn[1]  for cn in zip(bygroup['class_name'].values, list(bygroup['class_id']==2)) ], ['image_id', 'rad_id']]\nstrange_cases = pd.merge(train, strange_cases, on=['image_id', 'rad_id'], how='inner').sort_values(['image_id', 'rad_id', 'class_name']).reset_index(drop=True)\nstrange_cases","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have also [been told](https://www.kaggle.com/c/vinbigdata-chest-xray-abnormalities-detection/discussion/208111) that `PatientAge` in the `.dicom` meta data may be wrong. When I [did some exploratory data analysis](https://www.kaggle.com/bjoernholzhauer/eda-dicom-reading-vinbigdata-chest-x-ray#6.-What-is-in-the-.dicom-meta-data?), I saw 839 .dicom files with PatientAge given as 000Y (plus 000D = 0 days?), which could in theory be yet another placeholder for missing or might really mean 0 years of age (i.e. < 1 year-old). There's the issue of there being empty strings `''`, `000` as well as `NaN` string. Additionally, some of the ages are below < 18 years, e.g. there's 107 images with 001Y to 017Y. This is obbviously interesting, because obviously some diagnoses are just much more likely at certain ages. From the paper describing the methods to create the dataset, it appears that all data should be from adults\n> \"The collected raw data was mostly of adult PA-view CXRs, but also included a significant amount of outliers such as images of body parts other than chest (due to mismatched DICOM tags), pediatric scans, low-quality images, or lateral CXRs. […] All outliers were automatically excluded from the dataset […].\"\n\nThat prompted me to ask for clarification on the forum and the response suggested that these are data issues, which makes you wonder how many other ages above 18 may also be wrong.\n\nCertainly, several of the ages substantially above 100 years may also be doubtful."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def string_to_float_process(string):\n    if (string=='') | (string=='000D') | (string=='000'):\n        return float(\"NaN\")\n    else:\n        return float(string)\n\nPatientAge = [string_to_float_process(string) for string in dicom_meta['PatientAge'].replace(re.compile('Y$'), '') ]\nsns.distplot(PatientAge);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Bounding boxes with more than one label\nA lot of these are probably not really data issues. Quite often the same bounding boxes are given a lot of labels. As per the competition host this is legitimate. As pointed out in a [discussion thread](https://www.kaggle.com/c/vinbigdata-chest-xray-abnormalities-detection/discussion/212859#1161967), they have stated:\n> Some dense lesion area contains many labels, our radiologists created a box then assign many labels to that box.\n\nThe only question is whether some of these multiple labels are accidental or erroneous, after all. Let's have a look at a list of how often each combination of labels occurs:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"tmp1 = train.copy()\ntmp_cols = ['image_id', 'rad_id', 'x_min', 'x_max', 'y_min', 'y_max']\ntmp1 = tmp1[tmp_cols+[ 'class_id']]\\\n    .groupby(tmp_cols)\\\n    .count()\\\n    .reset_index()\\\n    .rename(columns={'class_id': 'Number of records'})\n#tmp1 = tmp1[tmp1['Number of records']>1]\ntmp1 = pd.merge(train.copy(), tmp1, \n                on=tmp_cols, \n                how='inner')\\\n    .sort_values(tmp_cols+['class_id'])\\\n    .reset_index(drop=True)\ntmp1['class_id_name'] = [str(row['class_id']) + ': ' + row['class_name'] for idx, row in tmp1.iterrows()]\n\n\ntmp1[\"grouprank\"] = tmp1.groupby(tmp_cols)['class_id']\\\n    .rank(method=\"first\", ascending=True)\\\n    .astype(np.int8)\n\ntmp1 = tmp1.pivot(index=tmp_cols+['Number of records'],\n           columns='grouprank',values='class_id_name')\\\n    .add_prefix('class_')\\\n    .reset_index()\n   \ndef isNaN(string):\n    return string != string\n\ntmp1['class_combinations'] = [ ' + '.join([  row['class_'+str(i)]  for i in range(1,6) if (not( isNaN(row['class_'+str(i)])))])  for idx, row in tmp1.iterrows() ]\npd.set_option('display.max_rows', 1000)\ntmp1['class_combinations'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see the most common thing a single bounding box is labelled as when there's more than one label is **10: Pleural effusion + 11: Pleural thickening** (1285 cases) followed by **6: Infiltration + 13: Pulmonary fibrosis** (318 cases), and **6: Infiltration + 7: Lung Opacity** (216 cases). Does this fall under \"dense lesion areas\"?\n\nSome of the combinations are clearly mistakes such as **0: Aortic enlargement + 3: Cardiomegaly** and do not seem very plausible such as **3: Cardiomegaly + 8: Nodule/Mass**."},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# What's the idea in this notebook for finding the issues and mislabeled bounding boxes?\n\nWhat are we trying to address? Two possible questions are:\n* Are the labels wrong for an existing bounding box (metrics: Do we predict this label given the characteristics of the bounding box? This is an obvious question is whether we can more systematically do what humans did there more or less manually: identify implausible bounding box labels. In that case we make the assumption that the position of bounding boxes are right, we are just not sure whether they are mislabeled.\n* Another question could be whether certain bounding boxes are in the wrong place. We could look at that by looking at whether other radiologists give this label for overlapping boxes or even have overlapping boxes, at all. A metric could be IOU.\n\nIn this notebook I use LightGBM to try and find misclassified bounding boxes. For classifying the `class_id` for a bounding box, categorical cross entropy loss is one obvious choice, but multi-class focal loss could be an alternative. Obviously, when we optimize hyperparameters and evaluate performance we need to do cross-validation by `image_id` to avoid issues with leakage. On the other hand, we want to **stay away** from `rad_id` as a feature - it may well be a great feature, but we don't want to train a model to repeat a radiologist's mistakes.\n\n# Other possible approaches I did not try, yet\n\nAlternatively, we could try neural networks, where we could even have an `image_id` embedding or even use the actual image as an input. One possibility would simply be a classification model and to look at it as a multi-label problem - and then look for the labels the model considers the most \"surprising\" (aka most confused predictions).\n\nWe could also try to predict IOU or some other metric on average vs. the two other radiologists that labelled the same image - possibly with the actually labelled class for the bounding box as an input. \n\n# Other possible uses of these approaches\nAnother thought that occured to me is that we could dobule-check the predictions we get from models for plausibility using a very different approach. Sure, you'd have to think quite carefully how skeptical these models would have to be for you to override a model prediction from a vision model using the actual image, but with some careful tuning this might very well be useful."},{"metadata":{},"cell_type":"markdown","source":"# Feature engineering"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create a list of classes and a dictionary\nclasses = train[['class_id', 'class_name', 'rad_id']]\\\n    .groupby(['class_id', 'class_name'])\\\n    .count()\\\n    .rename(columns={'rad_id': 'Number of records'})\\\n    .reset_index()\n\nfor index, row in classes.iterrows():\n    if index==0:\n        label_dict = {row['class_id']: row['class_name']}\n    else:\n        label_dict.update({row['class_id']: row['class_name']})\n\ntrain = pd.merge(train, dicom_meta, on='image_id', how='left')\ntrain['x_max'] = train['x_max']/train['Columns']\ntrain['x_min'] = train['x_min']/train['Columns']\ntrain['y_max'] = train['y_max']/train['Rows']\ntrain['y_min'] = train['y_min']/train['Rows']\ntrain['width'] = (train['x_max']-train['x_min'])\ntrain['height'] = (train['y_max']-train['y_min'])\ntrain['area'] = train['height']*train['width']\ntrain['x_center'] = (train['x_max']+train['x_min'])/2\ntrain['y_center'] = (train['y_max']+train['y_min'])/2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# In this cell, we create a feature for how many other classes a radiologist has already assigned \n# for the same image.\n\ntmp1 = pd.merge(train, train, on=['image_id', 'rad_id'], how='left').fillna(0)\n#tmp1 = tmp1[ (tmp1['class_id_x']!=tmp1['class_id_y']) | (tmp1['x_min_x']!=tmp1['x_min_y']) | (tmp1['x_max_x']!=tmp1['x_max_y']) | (tmp1['y_min_x']!=tmp1['y_min_y']) | (tmp1['y_max_x']!=tmp1['y_max_y'])]\n\ntmp_cols = ['image_id', 'rad_id', 'class_id_x', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x']\ntmp1 = tmp1[ tmp_cols + ['class_id_y', 'x_max_y']]\\\n    .groupby(tmp_cols+['class_id_y'])\\\n    .count()\\\n    .reset_index()\\\n    .pivot(index=tmp_cols,\n           columns='class_id_y', values='x_max_y')\\\n    .add_prefix('other_')\\\n    .reset_index()\\\n    .rename(columns={'class_id_x':'class_id',\n                     'class_id_x': 'class_id',\n                     'x_min_x': 'x_min',\n                     'y_min_x': 'y_min',\n                     'x_max_x': 'x_max',\n                     'y_max_x': 'y_max'})\n\ntrain = pd.merge(train, tmp1, \n                 on=['image_id', 'rad_id', 'class_id', 'x_min', 'y_min', 'x_max', 'y_max'],\n                 how='left')\ntrain[['other_'+str(i) for i in range(15)]] = train[['other_'+str(i) for i in range(15)]]\\\n    .fillna(0)\\\n    .astype(np.int)\n\n# Finally, we subtract the extra count of +1 for each label itself (when we want to predict it, we\n# do not want to have a feature that leaks the label, which it otherwise would).\nfor idx, row in train.iterrows():\n    if row['class_id']<14:        \n        train['other_' + str(row['class_id'])].values[idx] += -1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Features based on where bounding boxes for a class tend to be"},{"metadata":{"trusted":true},"cell_type":"code","source":"locations = np.zeros((14, 1000, 1000))\nfor index, row in tqdm(train.iterrows(), total=train.shape[0]):\n    if row['class_id']<14:\n        locations[row['class_id'], \n                  ((np.round(row['y_min'],3)*1000).astype(np.int)):((np.round(row['y_max'],3)*1000).astype(np.int)), \n                  ((np.round(row['x_min'],3)*1000).astype(np.int)):((np.round(row['x_max'],3)*1000).astype(np.int))] += 1\n        \nclasscounts = train[['image_id', 'rad_id', 'class_id','class_name']]\\\n    .groupby(['image_id', 'rad_id', 'class_id'])\\\n    .count()\\\n    .reset_index()\\\n    .pivot(index=['image_id', 'rad_id'], columns='class_id', values='class_name')\\\n    .rename(columns={i:'n_class'+str(i) for i in range(15)})\\\n    .fillna(0)\n\nclassareas = train[['image_id', 'rad_id', 'class_id','area']]\\\n    .groupby(['image_id', 'rad_id', 'class_id'])\\\n    .sum()\\\n    .reset_index()\\\n    .pivot(index=['image_id', 'rad_id'], columns='class_id', values='area')\\\n    .rename(columns={i:'area_class'+str(i) for i in range(15)})\\\n    .fillna(0)\n\ntrain = pd.merge( pd.merge( train, classcounts, on=['image_id', 'rad_id'], how='left'), \n                  classareas, on=['image_id', 'rad_id'], how='left')\ntrain = train[train['class_id']!=14]\n\nclasses = train[['class_id', 'class_name', 'rad_id']]\\\n    .groupby(['class_id', 'class_name'])\\\n    .count()\\\n    .rename(columns={'rad_id': 'Number of records'})\\\n    .reset_index()\n\nfor index, row in classes.iterrows():\n    if index==0:\n        label_dict = {row['class_id']: row['class_name']}\n    else:\n        label_dict.update({row['class_id']: row['class_name']})\n        \nf, axs = plt.subplots(5, 3, sharey=True, sharex=True, figsize=(16,28));\n\nfor class_id in range(14):\n    axs[class_id // 3, class_id - 3*(class_id // 3)].imshow(locations[class_id], cmap='inferno', interpolation='nearest');\n    axs[class_id // 3, class_id - 3*(class_id // 3)].set_title(str(class_id) + ': ' + label_dict[class_id])\n    \nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's zero out really low counts (<=3) - especialyl to avoid target leakage\nfor class_id in range(14):\n    for x in range(1000):\n        for y in range(1000):\n            if locations[class_id,x,y]<=3:\n                locations[class_id,x,y] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for class_id in range(14):\n    train['overlap'+str(class_id)] = 0\n    \nfor index, row in tqdm(train.iterrows(), total=train.shape[0]):    \n    for class_id in range(14):\n        train.loc[index, 'overlap'+str(class_id)] = np.mean( locations[class_id, int(np.floor(1000*row['x_min'])):int(np.ceil(1000*row['x_max'])), int(np.floor(1000*row['y_min'])):int(np.ceil(1000*row['y_max']))] )\n        \n#overlaps = ['overlap'+str(class_id) for class_id in range(14)]\n#[ (class_id, train[overlaps[class_id]].values[idx]) for idx, class_id in enumerate(train['class_id'].values) if (class_id!=14) & (train[overlaps[class_id]].values[idx]>0)]        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Creating our features\n\nNow we get to features based what the other two radiologists said that looked at the same image:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef bbox_overlap(bbox1, bbox2):\n    area1 = (bbox1['x_max']-bbox1['x_min'])*(bbox1['y_max']-bbox1['y_min'])\n    area2 = (bbox2['x_max']-bbox2['x_min'])*(bbox2['y_max']-bbox2['y_min'])\n    intersection = max(0, min(bbox1['x_max'], bbox2['x_max']) - max(bbox1['x_min'], bbox2['x_min'])) * max(0, min(bbox1['y_max'], bbox2['y_max']) - max(bbox1['y_min'], bbox2['y_min']))\n    return intersection / ( area1 + area2 - intersection)\n\ntmp1 = pd.merge(train, train, on='image_id', how='inner')\n\n# Feature counting other classes assigne by the other two radiologists reviewing the image\nclass_counts = tmp1.loc[ (tmp1['rad_id_x']!=tmp1['rad_id_y']), ['image_id', 'rad_id_x', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x', 'rad_id_y']]\\\n    .groupby(['image_id', 'rad_id_x', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x'])\\\n    .count()\\\n    .reset_index()\\\n    .rename(columns={'rad_id_x': 'rad_id', 'class_id_x': 'class_id', 'x_min_x': 'x_min', 'x_max_x':'x_max', 'y_min_x': 'y_min', 'y_max_x': 'y_max', 'rad_id_y': 'class_count'})\\\n    .pivot(index=['image_id', 'rad_id', 'class_id', 'x_min', 'y_min', 'x_max', 'y_max'], columns='class_id_y', values='class_count')\\\n    .fillna(0)\\\n    .reset_index()\n\ntmp1 = tmp1[ (tmp1['rad_id_x']!=tmp1['rad_id_y']) & (tmp1['class_id_x']!=14)  & (tmp1['class_id_y']!=14)]\ntmp1 = tmp1[['image_id', 'rad_id_x', 'rad_id_y', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x', 'x_min_y', 'y_min_y', 'x_max_y', 'y_max_y']]\n\n# Derive intersection-over-union vs. other bounding boxes by the other radiologists\ntmp1['iou'] = [bbox_overlap({'x_min': row['x_min_x'], 'x_max': row['x_max_x'], 'y_min':row['y_min_x'], 'y_max':row['y_max_x']}, {'x_min': row['x_min_y'], 'x_max': row['x_max_y'], 'y_min':row['y_min_y'], 'y_max':row['y_max_y']}) for idx, row in tmp1.iterrows()]\n\n# What is the maximum IOU?\nmaxiou = tmp1[['image_id', 'rad_id_x', 'rad_id_y', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x', 'iou']]\\\n    .groupby(['image_id', 'rad_id_x', 'rad_id_y', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x'])\\\n    .max()\\\n    .reset_index()\\\n    .rename(columns={'rad_id_x': 'rad_id', 'class_id_x': 'class_id', 'x_min_x': 'x_min', 'x_max_x':'x_max', 'y_min_x': 'y_min', 'y_max_x': 'y_max'})\n\n# What is the sum of IOU?\nsumiou = tmp1[['image_id', 'rad_id_x', 'rad_id_y', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x', 'iou']]\\\n    .groupby(['image_id', 'rad_id_x', 'rad_id_y', 'class_id_x', 'class_id_y', 'x_min_x', 'y_min_x', 'x_max_x', 'y_max_x'])\\\n    .sum()\\\n    .reset_index()\\\n    .rename(columns={'rad_id_x': 'rad_id', 'class_id_x': 'class_id', 'x_min_x': 'x_min', 'x_max_x':'x_max', 'y_min_x': 'y_min', 'y_max_x': 'y_max'})\n\ntmp_cols = ['image_id', 'rad_id', 'class_id', 'class_id_y', 'x_min', 'y_min', 'x_max', 'y_max']\n# Take maximum or min or mean over the two other radiologists for these two metrics.\nmaxsumiou = sumiou[tmp_cols+['iou']].groupby(tmp_cols).max().reset_index()\nminsumiou = sumiou[tmp_cols+['iou']].groupby(tmp_cols).min().reset_index()\nmeansumiou = sumiou[tmp_cols+['iou']].groupby(tmp_cols).mean().reset_index()\nmaxmaxiou = maxiou[tmp_cols+['iou']].groupby(tmp_cols).max().reset_index()\nminmaxiou = maxiou[tmp_cols+['iou']].groupby(tmp_cols).min().reset_index()\nmeanmaxiou = maxiou[tmp_cols+['iou']].groupby(tmp_cols).mean().reset_index()\n\ntmp_cols = ['image_id', 'rad_id', 'class_id', 'x_min', 'y_min', 'x_max', 'y_max']\nmaxsumiou = maxsumiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'maxsumiou'+str(i) for i in range(14)})\nminsumiou = minsumiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'minsumiou'+str(i) for i in range(14)})\nmeansumiou = meansumiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'meansumiou'+str(i) for i in range(14)})\nmaxmaxiou = maxmaxiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'maxmaxiou'+str(i) for i in range(14)})\nminmaxiou = minmaxiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'minmaxiou'+str(i) for i in range(14)})\nmeanmaxiou = meanmaxiou\\\n    .pivot(index=tmp_cols, columns='class_id_y', values='iou')\\\n    .fillna(0)\\\n    .reset_index()\\\n    .rename(columns={i:'meanmaxiou'+str(i) for i in range(14)})\n\ntrain = pd.merge(pd.merge(pd.merge(pd.merge(pd.merge(pd.merge(pd.merge(train,\n    class_counts, on=tmp_cols, how='left'),\n    maxsumiou, on=tmp_cols, how='left'),\n    minsumiou, on=tmp_cols, how='left'),\n    meansumiou,on=tmp_cols, how='left'),\n    maxmaxiou, on=tmp_cols, how='left'),\n    minmaxiou, on=tmp_cols, how='left'),\n    meanmaxiou,on=tmp_cols, how='left')\ntrain = train[train['class_id']!=14].fillna(0)\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#PatientAge, PatientSex, PatientSize, PatientWeight\n#list_of_ages = [string_to_float_process(string) for string in train_files['PatientAge'].replace(re.compile('Y$'), '') ]\n#int(train['PatientAge'])\n\ndef string_to_float_process(string):\n    if (string=='') | (string=='000D') | (string=='000'):\n        return float(\"NaN\")\n    else:\n        return float(string)\n\ntrain['Age'] = [string_to_float_process(string) for string in train['PatientAge'].replace(re.compile('Y$'), '') ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['Sex'] = train['PatientSex'].fillna('Q').map({'M':0, 'F':3, 'O':1, 'Q':2}).astype('category')\ntrain = train.reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Set-up for modelling including CV based on image_id"},{"metadata":{"trusted":true},"cell_type":"code","source":"kf = KFold(n_splits=5)\n\nimage_ids = np.unique(train['image_id'].values)\nfold_assignments = { image_ids[val_idx]: fold for fold, idxs in enumerate(kf.split(image_ids)) for val_idx in idxs[1]} #train_index, test_index\ntrain['fold'] = train['image_id'].map(fold_assignments)\n\n# LightGBM CV takes folds (generator or iterator of (train_idx, test_idx) tuples\nfold_splits = [ [train[train['fold']!=fff].index.tolist(), train[train['fold']==fff].index.tolist()] for fff in range(5) ]\n\ncont_features = ['Age', 'x_min', 'y_min', 'x_max', 'y_max', 'width', 'height', 'area', 'x_center', 'y_center']\\\n    + [f'minsumiou{i}' for i in range(14)]\\\n    + [f'maxsumiou{i}' for i in range(14)]\\\n    + [f'meansumiou{i}' for i in range(14)]\\\n    + [f'minmaxiou{i}' for i in range(14)]\\\n    + [f'maxmaxiou{i}' for i in range(14)]\\\n    + [f'meanmaxiou{i}' for i in range(14)]\\\n    + [f'overlap{i}' for i in range(14)]\\\n    + [f'other_{i}' for i in range(14)]\n#cat_features = ['Sex']\n\nX = np.array( train[cont_features] )#+ cat_features\ny = np.array( train['class_id'] ).flatten()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# LightGBM modeling\n\nI did an initial hyperparameter search using the `optuna.LightGBMTunerCV` function that implements a sensible tuning strategy for LightGBM. I'm now commenting this out, because it's a bit time-consuming.\n\nIf you want to find out more about tuning LightGBM with `optuna`, I've got [a notebook on this](https://www.kaggle.com/bjoernholzhauer/lightgbm-tuning-with-optuna) using the Titanic dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"#import optuna.integration.lightgbm as lgb\n#import optuna\n\n#dtrain = lgb.Dataset(X, label=y)\n\n#params = {\n#    'objective': 'multiclass',\n#    'boosting_type': 'gbdt',    \n#    'num_class':14,\n#    'metric': 'multi_logloss',#['multi_logloss', 'multi_error'],\n#    'learning_rate': 0.05, #0.005,\n#}\n\n#study_tuner = optuna.create_study(direction='minimize')\n\n#tuner = lgb.LightGBMTunerCV(params, dtrain,                             \n#                            study=study_tuner,\n#                            verbose_eval=10, #True, #False,\n#                            #categorical_feature=cat_features,\n#                            early_stopping_rounds=25,\n#                            time_budget=19800, # Time budget of 5 hours, we will not really need it\n#                            seed = 42,\n#                            folds=fold_splits,\n#                            num_boost_round=800, #10000,\n#                            #callbacks=[lgb.reset_parameter(learning_rate = [0.005]*200 + [0.001]*9800) ] #[0.1]*5 + [0.05]*15 + [0.01]*45 + \n#                           )\n#tuner.run()\n#print(tuner.best_params)\n#print(tuner.best_score)\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The results I received were\n> {'objective': 'multiclass', \n>  'boosting_type': 'gbdt', \n>  'num_class': 14, \n>  'metric': 'multi_logloss', \n>  'learning_rate': 0.05, \n>  'feature_pre_filter': False, \n>  'lambda_l1': 4.053841983827334, \n>  'lambda_l2': 8.075199666030122e-06, \n>  'num_leaves': 25, \n>  'feature_fraction': 0.5, \n>  'bagging_fraction': 0.875336404331162, \n>  'bagging_freq': 6, \n>  'min_child_samples': 100}\n\nAnd the best log-loss was:\n> Best Score: 0.7220672051541743"},{"metadata":{},"cell_type":"markdown","source":"Now let's round these parameter values a bit and fit the model with a lower learning rate:"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\n\nparams = {'objective': 'multiclass', \n 'boosting_type': 'gbdt', \n 'num_class': 14, 'metric': ['multi_logloss', 'multi_error'],\n 'learning_rate': 0.005, \n 'feature_pre_filter': False, \n 'lambda_l1': 4, \n 'lambda_l2': 8, \n 'num_leaves': 25, \n 'feature_fraction': 0.5, \n 'bagging_fraction': 0.875, \n 'bagging_freq': 6, \n 'min_child_samples': 100,\n 'num_threads': 4}\n\ndtrain = lgb.Dataset(X, label=y)\n\n#lgbcv = lgb.cv(params, dtrain, num_boost_round=5000, folds=fold_splits,verbose_eval=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Doing cross-validation to see a good number of iterations suggests around 3950:\n> [3948]\tcv_agg's multi_logloss: 0.719552 + 0.019352\tcv_agg's multi_error: 0.250182 + 0.00657187\n\n> [3949]\tcv_agg's multi_logloss: 0.719551 + 0.0193541\tcv_agg's multi_error: 0.250182 + 0.00657187\n\n> [3950]\tcv_agg's multi_logloss: 0.719552 + 0.0193557\tcv_agg's multi_error: 0.250209 + 0.00659723\n\n> [3951]\tcv_agg's multi_logloss: 0.719553 + 0.019359\tcv_agg's multi_error: 0.250181 + 0.00660947\n\nSo we'll do that a refit without CV. We'll show the training error (since we no longer have a validation dataset)."},{"metadata":{"trusted":true},"cell_type":"code","source":"lgbfit = lgb.train(params, dtrain, num_boost_round=3950, verbose_eval=250,valid_sets=[dtrain])\n# We only add the valid_sets to get some updates as the model is trained.\n# It only shows performance on training set.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What labels does the model find weird?"},{"metadata":{},"cell_type":"markdown","source":"As we can see in the plots below, some bounding boxes seem to be very distinctive (i.e. the model manages to assign high probabilities to them ones with this `class_id`). \n\nThese are primarily **Aortic enlargment**, **Cardiomegaly** and **Pneumothorax** and for these classes there are then some observations with very low probabilities that we should regard as under suspicion for being mislabelled. Of course, this could also be a split into easy and difficult to label cases/bounding boxes. I'm especially unsure for pneumothorax, because there's quite a lot of low-probability records.\n\nThen there's somewhat bimodal distributions of probabilities e.g. for **Pleural effusion** and **Pleural thickening**. Such a phenomenon may also be occuring for **Other lesion**, **Calcification** and **Consolidation**.\n\nThe other classes seem to have more of a smooth distribution."},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = lgbfit.predict(X)\ntrain['prob'] = [ preds[y] for idx, (y, preds) in enumerate(zip(y, preds))]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['first_pred'] = np.argsort(preds, axis=1)[:,-1:].flatten()\ntrain['second_pred'] = np.argsort(preds, axis=1)[:,-2:-1].flatten()\ntrain['alternative'] = (train['first_pred']==train['class_id']) * train['second_pred']\\\n    + (train['first_pred']!=train['class_id'])*train['first_pred']\ntrain['first_prob'] = np.sort(preds, axis=1)[:,-1:].flatten()\ntrain['second_prob'] =np.sort(preds, axis=1)[:,-2:-1].flatten()\ntrain['alternative_prob'] = (train['first_pred']==train['class_id']) * train['second_prob']\\\n    + (train['first_pred']!=train['class_id'])*train['first_prob']\n\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = sns.distplot(train[\"prob\"], kde=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#for class_id in range(7):\n#    ax = sns.kdeplot(train.rename(columns={'prob':f'Class_{class_id}'}).loc[train['class_id']==class_id, f\"Class_{class_id}\"]);\n#plt.ylim(0, 5);\n#plt.xlim(0, 1);\n\ntrain['class_id_name'] = [f'{class_id}: {class_name}' for class_id, class_name in zip(train['class_id'], train['class_name'])]\ng = sns.FacetGrid(train.sort_values('class_id'), col='class_id_name', \n                  sharex=False, sharey=False,\n                  col_wrap=3, height=5);\ng.map(plt.hist, 'prob', bins=30);\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['logit'] = np.log(train['prob']) - np.log1p(-train['prob'])\ng = sns.FacetGrid(train.sort_values('class_id'), col='class_id_name', \n                  sharex=False, sharey=False,\n                  col_wrap=3, height=5);\ng.map(plt.hist, 'logit', bins=50);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['flagged'] = ((train['class_id'].isin([0,3])) & (train['logit']< -2))\\\n    | ((train['class_id'].isin([1,2,4,5,7,9,11,13])) & (train['logit']< -6))\\\n    | ((train['class_id'].isin([6,8,10,12])) & (train['logit']< -5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.to_csv('train_with_flags.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[train['flagged']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Getting SHAP explanations for cases the model considers implausible\n\nTo find out more about SHAP, have a look at [this notebook](https://www.kaggle.com/bjoernholzhauer/will-rose-or-jack-survive-lightgbm-and-shap) explaining LightGBM and SHAP.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"explainer = shap.TreeExplainer(lgbfit)\nshap_values = explainer.shap_values(train[cont_features].sample(n=10000, replace=False, random_state=42))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Color-scheme from http://www.randalolson.com/2014/06/28/how-to-make-beautiful-data-visualizations-in-python-with-matplotlib/\ntableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),  \n             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),  \n             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),  \n             (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),  \n             (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]\nfor i in range(len(tableau20)):  \n    r, g, b = tableau20[i]  \n    tableau20[i] = (r / 255., g / 255., b / 255.)\n    \nshap.summary_plot(shap_values, \n                  train[cont_features], \n                  color=matplotlib.colors.ListedColormap(tableau20))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"flagged_idxs = np.where( train['flagged'] == True)[0]\nflagged_idxs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_image(image_id, rad_id):\n    with shelve.open('../input/eda-dicom-reading-vinbigdata-chest-x-ray/training_data.db', \n                     flag='r', writeback=False) as myshelf:\n        tmpdict = myshelf[image_id]                        \n    which_indices = [idx for idx, val in enumerate(tmpdict['rad_id']) if val==rad_id]\n    image = np.stack([tmpdict['image']]*3).transpose(1,2,0)        \n    bboxes = tmpdict['bboxes'][which_indices]\n    class_labels = tmpdict['class_labels'][which_indices]        \n    return {'image': image, 'bboxes': bboxes, 'class_labels': class_labels}\n\nfor index, row in train.loc[train['flagged'], \n                            ['image_id', 'rad_id', 'class_id', \n                             'class_name', 'prob', 'alternative', 'alternative_prob']]\\\n        .drop_duplicates()\\\n        .reset_index(drop=True)\\\n        .iterrows():\n    image_id = row['image_id']    \n    problem_class_id = row['class_id']\n    problem_class = str(row['class_id']) + \": \" + row['class_name']\n    rad_id = int(re.findall(r'\\d+', row['rad_id'])[0])\n    im = get_image(image_id, rad_id)\n    alternative = str(row['alternative']) + \": \" + label_dict[row['alternative']]\n    alt_prob = np.round(row['alternative_prob'],8)\n    prob = np.round(row['prob'], 8)\n\n    plt.figure(figsize=(20,10));    \n    plt.imshow(im['image']);\n    plt.suptitle(f'Image {image_id}, radiologist {rad_id}\\n red: {problem_class} (prob={prob})\\nAlternative: {alternative} (prob={alt_prob})', fontsize=16)    \n    ax = plt.gca();\n\n    for bbox, class_id in zip(im['bboxes'], im['class_labels']):\n        #print(f'{image_id}, {rad_id}, {bbox}, {class_id}')        \n        if class_id==problem_class_id: \n            ecol='r'\n            linestyle = 'dashed'\n            offset=12\n        else:\n            ecol='b'\n            linestyle = 'dotted'\n            offset=24\n\n        # Create a Rectangle patch\n        rect = Rectangle((bbox[0],bbox[1]),bbox[2]-bbox[0],bbox[3]-bbox[1],\n                         linewidth=1,linestyle=linestyle,\n                         edgecolor=ecol,facecolor='none');\n        # Add the patch to the Axes\n        ax.add_patch(rect);\n        ax.annotate(str(class_id) + ': ' + label_dict[class_id], \n                    xy=((bbox[2]+bbox[0])/2, bbox[1]+offset), \n                    xycoords='data', color=ecol);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Individual explanations for predictions\nLet's also look at 20 specific SHAP explanations for model predictions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"shap.initjs()\ndef shap_explain_flgidx(flgidx):\n    print(train.iloc[flgidx])\n    class_id = train.iloc[flgidx]['class_id']\n    return shap.force_plot(explainer.expected_value[class_id],\n                shap_values[class_id][flgidx],\n                train[cont_features].iloc[flgidx])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[4])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[6])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[7])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[8])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[9])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[10])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[11])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[12])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[13])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[14])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[15])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[16])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[17])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[18])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shap_explain_flgidx(flagged_idxs[19])","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}