{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 🩺 EDA 📊","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In this Notebook we get to know the data we are working with in this competition. The main steps are:\n\nI. [Clean tabular data](#clean_data)\n\nII. [Information about patients](#patients)\n\nIII. [Investigate the progress of some patients](#progress)\n\nIV. [Data from DICOM files](#dicom)","execution_count":null},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install missingpy","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nfrom IPython.display import Image as show_gif\nfrom sklearn.cluster import KMeans\nfrom skimage.transform import resize\nimport scipy.ndimage as ndimage\nfrom plotly.tools import FigureFactory as FF\nfrom skimage import measure, morphology, segmentation\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom tqdm.notebook import tqdm\nfrom sklearn.metrics import roc_auc_score, mean_squared_error\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.linear_model import LinearRegression\nfrom missingpy import MissForest\nimport lightgbm as lgb\nfrom scipy import stats\nimport copy\nimport pydicom\nimport glob\nimport re\nimport os\nimport scipy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets have a first look at the data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# I. Clean tabular data <a id='clean_data'></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## I. Duplicates","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":" train[train.duplicated(subset=['Patient','Weeks'], keep=False)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks like the dataset contains duplicate rows with different FVC and Percent values. We remove those rows by calculating the average for FVC and Percent","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train = train.groupby(['Patient', 'Weeks']).agg({ \n    'FVC': 'mean', \n    'Percent': 'mean', \n    'Age': 'first',\n    'Sex': 'first',\n    'SmokingStatus': 'first'\n}).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Number of duplicates: %s'%len(train[train.duplicated(subset=['Patient','Weeks'], keep=False)]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## II. Missing data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train.isna().any()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So none of the columns contain any missing values! ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# II. Information about patients <a id='patients'></a>","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Number of unique patients: %s\"%(train.Patient.nunique()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_obs_per_patient = train.groupby('Patient').count()['Weeks'].sort_values()\n\nfig, axs = plt.subplots(2, 1, figsize=(15, 10))\ndense = sns.kdeplot(num_obs_per_patient, bw=.5, ax=axs[0])\ndense.get_legend().remove()\ndense.set_xlabel(\"Number of Observations\")\nbar = sns.barplot(x=list(range(len(num_obs_per_patient))), y=num_obs_per_patient, ax=axs[1])\nbar.axes.get_xaxis().set_visible(False);\nbar.set_ylabel(\"Oberservations\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most patients have 9 oberservations in the dataset providing FVC values over time. ","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(2, 2, figsize=(15, 10))\n\naxs[0,0].set_title('Weeks', fontsize=18)\nsns.kdeplot(train.Weeks, shade=True, ax=axs[0,0])\naxs[0,0].get_legend().remove()\n\naxs[0,1].set_title('FVC', fontsize=18)\nsns.kdeplot(train.FVC, shade=True, ax=axs[0,1])\naxs[0,1].get_legend().remove()\n\naxs[1,0].set_title('Percent', fontsize=18)\nsns.kdeplot(train.Percent, shade=True, ax=axs[1,0])\naxs[1,0].get_legend().remove()\n\naxs[1,1].set_title('Age', fontsize=18)\nsns.kdeplot(train.Age, shade=True, ax=axs[1,1])\naxs[1,1].get_legend().remove()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def corrfunc(x, y, **kws):\n    r, _ = stats.pearsonr(x, y)\n    ax = plt.gca()\n    ax.annotate(\"r = {:.2f}\".format(r),\n                xy=(.1, .9), xycoords=ax.transAxes)\n\ng = sns.PairGrid(train, palette=[\"red\"])\ng.map_upper(plt.scatter, s=10)\ng.map_diag(sns.distplot, kde=False)\ng.map_lower(sns.kdeplot, cmap=\"Blues_d\")\ng.map_lower(corrfunc)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As seen in the correlation plot a certain correlation between FVC and Percent can be observed (0.67).","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# III. Investigate the progress of some patients <a id='progress'></a>","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"id_patients_most_weeks = train.groupby('Patient').Weeks.count().sort_values(ascending=False).iloc[:5].index\ntrain_patients = train[train.Patient.isin(id_patients_most_weeks)].sort_values('Weeks')\n\nfig, ax = plt.subplots(figsize=(16,6))\n\nfor name, group in train_patients.groupby('Patient'):\n    color = next(ax._get_lines.prop_cycler)['color']\n    group.sort_values('Weeks').plot(x='Weeks', y='FVC', ax=ax, label=name, color=color)\n    reg = LinearRegression().fit(np.array(group.Weeks).reshape(-1, 1), np.array(group.FVC))\n    ax.plot(group.Weeks,reg.predict(np.array(group.Weeks).reshape(-1, 1)),'--', color=color)\n    \nax.set_title('Progress of patients with the most FVC measurements', fontsize=18)\nax.set_xlabel('Weeks', fontsize=12)\nax.set_ylabel('FVC', fontsize=12);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see a steady decline in all patient lungs over time, as expected.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# IV. Data from DICOM files <a id='dicom'></a>","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"example_patient_file_path = \"../input/osic-pulmonary-fibrosis-progression/train/ID00014637202177757139317/\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## I. General <a id='dicom_general'></a>","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"num_scans = {}\nfor patient_id in train.Patient.unique():\n    files = glob.glob(\"../input/osic-pulmonary-fibrosis-progression/train/%s/*.dcm\"%patient_id)\n    num_scans[patient_id] = len(files)\ndf_scans = pd.DataFrame.from_dict(num_scans, orient='index', columns=['num_scans'])\ndf_scans.index = df_scans.index.rename('Patient_id')\ndf_scans = df_scans.sort_values('num_scans')\n\nfig, axs = plt.subplots(2, 1, figsize=(15, 10))\nfig.suptitle('Number of CT-Scans for every Patient id', fontsize=16)\nbar = sns.barplot(x=df_scans.index, y=df_scans.num_scans, ax=axs[0])\nbar.axes.get_xaxis().set_visible(False)\nbar.set_ylabel(\"Number of CT-Scan Files\")\n\ndense = sns.kdeplot(df_scans.num_scans, bw=50, ax=axs[1])\ndense.get_legend().remove()\ndense.set_xlabel('Numober of CT-Scan Files');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_scans.num_scans.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## II. Metadata","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### I. Extract ","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#https://github.com/pydicom/pydicom/issues/319\ndef dictify(ds):\n    \"\"\"Turn a pydicom Dataset into a dict with keys derived from the Element tags.\n\n    Parameters\n    ----------\n    ds : pydicom.dataset.Dataset\n        The Dataset to dictify\n\n    Returns\n    -------\n    output : dict\n    \"\"\"\n    output = dict()\n    for elem in ds:\n        # skip the image data\n        if elem.name=='Pixel Data':\n            continue\n        if elem.VR != 'SQ':\n            output[elem.name] = elem.value\n        else:\n            output[elem.name] = [dictify(item) for item in elem]\n    return output","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ds1 = dictify(pydicom.filereader.dcmread(example_patient_file_path+\"1.dcm\"))\nds2 = dictify(pydicom.filereader.dcmread(example_patient_file_path+\"2.dcm\"))\nds1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In each DICOM-File a lot of metadata can be found. The basic idea of the next steps is to determine attributes that don't change over DICOM-Files of the same patient and extract additional features that can be merged with the train.csv","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_diff(ds1, ds2):\n    diff_keys = []\n    for key in ds1.keys():\n        if ds1[key]!=ds2[key]:\n            diff_keys.append(key)\n    return diff_keys\n\nchanging_keys = get_diff(ds1, ds2)\nstable_keys = ds1.keys()-changing_keys","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stable_keys","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are **49** attributes that (propably) don't change over image files from the same patient. Now we extract all of those attributes by reading the first file of each patient. The attributes are stored in a data frame which will later get merged with the provided train.csv file","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"additional_metadata = []\nfor patient_id in tqdm(train.Patient.unique()):\n    file = glob.glob(\"../input/osic-pulmonary-fibrosis-progression/train/%s/*.dcm\"%patient_id)[0]\n    ds_dict = dictify(pydicom.filereader.dcmread(file))\n    metadata = {k:ds_dict[k] for k in stable_keys if k in ds_dict}\n    additional_metadata.append(metadata)\n    \ndf_add_meta = pd.DataFrame(additional_metadata)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for column in df_add_meta.columns:\n    if df_add_meta[column].dtypes!='float64':\n        try:\n            print(column, \":\", df_add_meta[column].unique()[:5], \"-\", df_add_meta[column].dtypes)\n        except:\n            print(column, \":\", df_add_meta[column].iloc[:10].values, \"-\", df_add_meta[column].dtypes)\n        print(\"-\"*100)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After a closer look of non numeric column some of them seem to have strange values. Take \"Convolution Kernel\" as an example:\n\nPandas reports that no unique values can be extracted because the values are of type 'MultiValue'. But if we show all values we can see that there is only one faulty entry","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_add_meta['Convolution Kernel'].values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One entry contains an array with an additional value. We can fix this by only using the first element of the array as value.","execution_count":null},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"column = 'Convolution Kernel'\n# strange way to get the single faulty value but it works \ndf_add_meta[column][(pd.isna(df_add_meta[column].str.contains('I50f'))).values] ='I50f'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\ncolumns_to_keep = [\"Patient Position\", \"Patient ID\", \"Manufacturer\", \"Convolution Kernel\"]\ncolumns_to_keep = list(df_add_meta.select_dtypes(include=numerics).columns.values) + columns_to_keep\ndf_add_meta = df_add_meta[columns_to_keep]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Check for missing values","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"perc_missing_cols = (df_add_meta.isna().sum()/len(df_add_meta)).sort_values(ascending=False)\nprint(\"Percentage of missing values in each column:\")\nprint(\"-\"*50)\nprint(perc_missing_cols[perc_missing_cols!=0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Many colums seem to have too many missing values as they would contain any relevant information. Therefore all columns are dropped that contain more than 10% missing values.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"perc_missing_cols = perc_missing_cols[perc_missing_cols<0.1]\ndf_add_meta = df_add_meta[df_add_meta.columns.intersection(list(perc_missing_cols.index))]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have a \"cleaned\" dataset from the additional metadata we can join them directly over the patient id.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train_join = pd.merge(train, df_add_meta, left_on='Patient', right_on=\"Patient ID\", how='inner')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next we need to check if those columns contain new missing values. This can happen when there is no value(s) for a patient. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"perc_missing_cols = (df_train_join.isna().sum()/len(df_train_join)).sort_values(ascending=False)\nprint(\"Percentage of missing values in each column:\")\nprint(\"-\"*50)\nprint(perc_missing_cols[perc_missing_cols!=0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Impute missing values using the MissForest-Algorithm","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"imputer = MissForest()\nnumerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n# create a new dataframe containing only numeric values which can be imputed\ndf_train_join_numeric = df_train_join.select_dtypes(include=numerics)\nimputed_matrix = imputer.fit_transform(df_train_join_numeric)\n# replace imputed columns in the df_add_meta dataframe\ndf_train_join_numeric = pd.DataFrame(imputed_matrix, index=df_train_join_numeric.index, columns=df_train_join_numeric.columns)\nfor column in df_train_join_numeric.columns:\n    df_train_join[column] = df_train_join_numeric[column]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The last step is now to reduce the memory size of the dataset by changing columns that are of type float to int if the column contains only integer values","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"for column in df_train_join.select_dtypes(include=numerics).columns:\n    if(df_train_join[column].apply(float.is_integer).all()):\n        df_train_join[column] = df_train_join[column].astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def evaluate_features(X_train, y_train, X_test, y_test, params, metrices):\n    \"\"\"\n    Trains a simple gradient boosting model and evaluates its feature importances (if multiple columns provided).\n    Furthermore the trained model is evaluated with the provided metric(es).\n    :param X_train:\n    :param y_train:\n    :param X_test:\n    :param y_test:\n    :param params:\n    :param metrices:\n    :return:\n    \"\"\"\n\n    for col in X_train.select_dtypes(include='object').columns:\n        le = LabelEncoder()\n        le.fit(list(X_train[col].astype(str).values) + list(X_test[col].astype(str).values))\n        X_train[col] = le.transform(list(X_train[col].astype(str).values))\n        X_test[col] = le.transform(list(X_test[col].astype(str).values))\n\n\n    clf = lgb.LGBMRegressor(**params)\n    clf.fit(X_train.values, y_train)\n\n    importances = clf.feature_importances_\n    indices = np.argsort(importances)[::-1]\n\n    features_to_show = len(X_train.columns)\n\n    plt.figure(figsize=(15,10))\n    plt.title(\"Feature importances\")\n    plt.bar(range(features_to_show), importances[indices][:features_to_show],\n            color=\"r\", align=\"center\")\n    feature_names = [X_train.columns[indices[f]] for f in range(features_to_show)]\n    plt.xticks(range(features_to_show), feature_names, rotation='vertical')\n    plt.xlim([-1, features_to_show])\n    plt.show()\n\n    scores = get_model_scores(clf, X_train, y_train, X_test, y_test, metrices, True)\n\n    df_feature_importance = pd.DataFrame({'column':X_train.columns[indices], 'importance':importances[indices]})\n    return (df_feature_importance, scores)\n\ndef get_model_scores(model, x_train, y_train, x_test, y_test, metrices, print_values=True):\n    scores = {}\n    for metric in metrices:\n        try:\n            score_train = metric(model.predict(x_train), y_train)\n            score_test = metric(model.predict(x_test), y_test)\n            if print_values:\n                print(metric.__name__, \"(train):\", score_train)\n                print(metric.__name__, \"(test):\", score_test)\n                print(\"------------------------------------------------------------\")\n            scores[metric.__name__] = [score_train, score_test]\n        except:\n            print(\"Could not calculate score\", metric.__name__)\n            print(\"------------------------------------------------------------\")\n            scores[metric.__name__] = [None, None]\n    return scores","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(df_train_join.drop(columns=['FVC']), df_train_join.FVC)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"evaluate_features(X_train, y_train, X_test, y_test, {}, [mean_squared_error])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Based on the feature importance of the Regression Tree many additional metadata attributes seem to hold a certain level of information to predict the <font size=\"4\">**current**</font> FVC value. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train_join.to_csv('train_merged_and_cleaned.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## III. Image data","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Next we look at the provided image data from the DICOM-Files.","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def center_crop(img, new_width=512, new_height=512):        \n\n    width = img.shape[1]\n    height = img.shape[0]\n\n    if new_width is None:\n        new_width = min(width, height)\n\n    if new_height is None:\n        new_height = min(width, height)\n\n    left = int(np.ceil((width - new_width) / 2))\n    right = width - int(np.floor((width - new_width) / 2))\n\n    top = int(np.ceil((height - new_height) / 2))\n    bottom = height - int(np.floor((height - new_height) / 2))\n\n    if len(img.shape) == 2:\n        center_cropped_img = img[top:bottom, left:right]\n    else:\n        center_cropped_img = img[top:bottom, left:right, ...]\n\n    return center_cropped_img\n\ndef slices_as_gif(slices, fps=10):\n    fig = plt.figure()\n    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)\n    plt.axis('off')\n    slices = [[plt.imshow(img, cmap='gray')] for img in slices]\n    ani = animation.ArtistAnimation(fig, slices, interval=200, repeat_delay=0)\n    ani.save('test_anim.gif', writer='imagemagick', fps=fps)\n    plt.close()\n    return 'test_anim.gif'\n    \ndef read_dicoms(file, with_mask=False):\n    img = pydicom.filereader.dcmread(file).pixel_array\n    img = center_crop(img)\n    if with_mask:\n        img = make_lungmask(img)\n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"files = glob.glob(example_patient_file_path+\"*.dcm\")\nfiles = sorted(files, key=lambda x:float(re.findall(\"(\\d+)\",x)[-1]))\nfig=plt.figure(figsize=(18, 6))\ncolumns = 10\nrows = 3\nfor i in range(30):\n    ds = pydicom.filereader.dcmread(files[i])\n    fig.add_subplot(rows, columns, i+1)\n    plt.imshow(center_crop(ds.pixel_array), cmap='gray')\n    plt.title(os.path.basename(files[i]))\n    plt.axis('off')\n    plt.grid(b=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slices = [read_dicoms(file) for file in files]\nshow_gif(filename=slices_as_gif(slices), format='png', width=512, height=512)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now that we have an overview of the data we can proceed with preprocessing the data and the images!**\n\n**Click [here](https://www.kaggle.com/foodaholic/image-preprocessing-tfrecords-with-3d-scan-2-3/edit) for my following notebook**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}