{"cells":[{"metadata":{},"cell_type":"markdown","source":"The goal of this notebook is to present an initial EDA regarding the data of **OSIC Pulmonary Fibrosis Progression** Kaggle challenge.\n\nFurthermore, in the end of the notebook, a **XGBoost** model is trained with the goal of predicting the FVC.This simple approach only takes in consideration the Patient's data (and not the CT scans), therefore, it is just an introductory and baseline approach, that shoud be improved in the future.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nfrom collections import Counter\nimport pydicom\nimport os\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.filters import threshold_otsu, median\nfrom scipy.ndimage import binary_fill_holes\nfrom skimage.segmentation import clear_border\nimport xgboost\nfrom sklearn.model_selection import train_test_split, GridSearchCV\nfrom sklearn.preprocessing import LabelEncoder, OneHotEncoder\nimport random\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        if 'ID00419637202311204720264' in dirname:\n#            print(os.path.join(dirname, filename))\n            \n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv')\ntest_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Exploratory Data Analysis","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Variables Distribution","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Patient\n\npatients_dist = Counter(train_df['Patient'].values).most_common()\n\nplt.figure(figsize=(10, 6))\nplt.bar(range(len(patients_dist)), list(dict(patients_dist).values()))\n\nplt.xlabel('Ids Patient')\nplt.ylabel('Observations per Patient')\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"info_df = pd.DataFrame(columns=['Patient','Age','Sex', 'SmokingStatus'])\n\nfor ind, row in train_df.groupby('Patient'):\n    new_row = {'Patient': ind, 'Age': row.iloc[0]['Age'], 'Sex': row.iloc[0]['Sex'], 'SmokingStatus': row.iloc[0]['SmokingStatus']}\n    info_df.loc[len(info_df)] = new_row\n    \n    \ninfo_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Age\n\nplt.figure(figsize = (10, 7))\nplt.hist(info_df['Age'], 10)\nplt.xlabel('Ages')\nplt.ylabel('Frequency')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#SmokingStatus\n\nplt.figure(figsize = (10, 7))\ninfo_df['SmokingStatus'].value_counts().plot(kind='bar');\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Gender\n\nplt.figure(figsize = (10, 7))\ninfo_df['Sex'].value_counts().plot(kind='bar');\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Weeks","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12, 8))\n\nfor ind, patient_data in train_df.groupby('Patient'):\n    plt.plot(patient_data['Weeks'], patient_data['FVC'], '.', label=ind)\n\nplt.xlabel('Weeks')\nplt.ylabel('FVC')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As observed, there is not a clear pattern on the FVC values over the number of weeks. This is understandable, as for each patient, the week = 0 correspond to a different stage of the decline in the lung function. This means that the weeks information cannot be crossed between different patients.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC evolution over time","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12, 8))\n\nc = ['blue', 'red', 'green']\ni = 0\n\nfor ind_smoke, status_data in train_df.groupby('SmokingStatus'):\n    counter = 0\n    print(ind_smoke,':', c[i], 'lines')\n    \n    for ind, patient_data in status_data.groupby('Patient'):\n        plt.plot(range(len(patient_data)), patient_data['FVC'], 'o-', label=ind_smoke, color=c[i])\n\n        counter += 1\n\n        if counter == 10:\n            break\n    i += 1\n\nplt.xlabel('Time')\nplt.ylabel('FVC')\n#plt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Overral, there is a decreasing trend on the FVC value over time","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Initial value of FVC","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"c = ['blue', 'red', 'green']\ni = 0\n\nfirst_FVC = {}\nfor ind_smoke, status_data in train_df.groupby('SmokingStatus'):\n    \n    first_FVC[ind_smoke] = []\n    for ind, patient_data in status_data.groupby('Patient'):\n        first_FVC[ind_smoke].append(patient_data.iloc[0]['FVC'])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12, 8))\n\nplt.hist(first_FVC['Ex-smoker'], bins = 15, label='Ex-smoker')\nplt.hist(first_FVC['Never smoked'], bins = 15, label='Never smoked')\nplt.hist(first_FVC['Currently smokes'], bins = 5, label='Currently smokes')\n\n\nplt.legend()\nplt.xlabel('FVC')\nplt.ylabel('Frequency')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"By analysing the first value of FVC of the different types of person (regarding the smoking habits), we can see that, generically, the FVC values are higher in the *Ex-smoker* and *Currently smokes* persons","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Percent","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\nplt.figure(figsize=(12, 8))\n\nfor ind, patient_data in train_df.groupby('Patient'):\n    plt.plot(patient_data['FVC'], patient_data['Percent'], '.', label=ind)\n\nplt.xlabel('FVC')\nplt.ylabel('Percent')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, there is a linear relationship between the FVC and Percent (which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics). This relationship may help to model the FVC","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## CT scans Analysis","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Vizualize the CTs","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ct_patient = 'ID00007637202177411956430'\n#ct_patient = 'ID00015637202177877247924'\n\nimg_paths = []\nfor dirname, _, filenames in os.walk('/kaggle/input/osic-pulmonary-fibrosis-progression/train'):\n    if ct_patient in dirname:\n        filenames.sort(key = lambda x: int(x.split('.')[0]))\n        for filename in filenames:\n            img_paths.append(os.path.join(dirname, filename))\n\n#img_paths","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#First image\n\nds = pydicom.dcmread(img_paths[0])\n\nplt.figure(figsize = (8,8))\nplt.imshow(ds.pixel_array, cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### CT Animation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Sequence of images\nfrom IPython.display import clear_output\n\n\nfor i in img_paths:\n    plt.figure(figsize = (8,8))\n    ds = pydicom.dcmread(i)\n    \n    plt.imshow(ds.pixel_array, cmap=plt.cm.bone)\n    plt.pause(0.1)\n    clear_output(wait=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Image Segmentation\n\nThe goal is to isolate the lungs region","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#from https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n\ndef load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 1\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n\npatient = load_scan('/kaggle/input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/')\nimgs = get_pixels_hu(patient)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"index_example = 20\n\nplt.imshow(imgs[index_example], cmap='gray')\nplt.title('HU Image')\nplt.show()\n\nds = pydicom.dcmread(img_paths[index_example])\nplt.imshow(ds.pixel_array, cmap='gray')\nplt.title('Original Image')\nplt.show()\n\nnew_img = imgs[index_example] < -500\n\nplt.imshow(new_img, cmap='gray')\nplt.title('Applying HU threshold')\nplt.show()\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see in the last image, the outside region is being considered as lung section, due to the HU pixels values. Therefore, a mask will be applied to filter outside boundary section.\n\nThen the lung volume will be computed (as suggested in https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165727)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"lungs = median(clear_border(new_img))\nlungs = morphology.binary_closing(lungs, selem=morphology.disk(7))\nmask = binary_fill_holes(lungs)\n\n\ndef compute_lung_volume(ds, mask):\n    return np.sum(mask) * float(ds.SliceThickness) * ds.PixelSpacing[0] * ds.PixelSpacing[1]\n\n\nprint(\"Lung Volume: \", compute_lung_volume(ds, mask))\n\nplt.imshow(mask, cmap='gray')\nplt.title('Applying HU threshold and Masking the outside region')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model to predict FVC\n\nA XGBoost model will be trained to try predicting the FVC, only using the patients data (that is, without analysing the CT images).\n\nSome preprocessing steps will be applied, like oversamling, features categorization and data normalization.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#oversampling, assuming that FVC follows a Linear Regression according to the Week value\n\nfor i, d in train_df.groupby('Patient'):\n    \n    week_data = d['Weeks']\n    \n    \n    xx = [c for c in range(min(week_data), max(week_data)) if c not in week_data]\n    \n    reg = LinearRegression().fit(week_data.values.reshape(-1,1), d['FVC'].values)\n    res = reg.predict(np.array(xx).reshape(-1,1))\n    \n    for j in range(len(xx)):\n        train_df.loc[len(train_df)] = [d.iloc[0]['Patient'], xx[j], res[j], d.iloc[0]['Percent'], d.iloc[0]['Age'], d.iloc[0]['Sex'], d.iloc[0]['SmokingStatus']]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Transforming Training Set\nX_train = pd.concat([train_df, test_df], axis = 0, ignore_index = True)\n\ny_train = X_train['FVC']\nX_train = X_train.drop(['FVC', 'Percent'], axis = 1)\n\nX_train.reset_index(inplace=True, drop=True)\n\n\n#Transforming Test Set\nX_test = pd.DataFrame(columns = ['Patient', 'Weeks', 'Age', 'Sex', 'SmokingStatus'])\n\nfor ind, row in test_df.iterrows():\n    \n    for i in range(-12, 133+1):\n        new_row = [row.Patient, i, row.Age, row.Sex, row.SmokingStatus]\n        X_test.loc[len(X_test)] = new_row","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ncategorical_vars = ['Patient', 'Sex', 'SmokingStatus']\nint_vars = ['Weeks', 'Age']\n\n\n#Categorize features\n\nencoder = OneHotEncoder(categories = 'auto', handle_unknown = 'ignore', sparse = False)\nencoder.fit(X_train[categorical_vars])\n\nX_train = pd.concat([X_train.drop(categorical_vars, axis = 1), pd.DataFrame(encoder.transform(X_train[categorical_vars]), \n                                           columns = encoder.get_feature_names())], axis=1, sort=False)\nX_test = pd.concat([X_test.drop(categorical_vars, axis = 1), pd.DataFrame(encoder.transform(X_test[categorical_vars]),\n                                         columns = encoder.get_feature_names())], axis=1, sort=False)\n\n\nX_test.Weeks = pd.to_numeric(X_test.Weeks)\nX_test.Age = pd.to_numeric(X_test.Age)\n\n\n\n#Normalization\n\nsc = StandardScaler()\nsc.fit(X_train[int_vars])\nX_train[int_vars] = sc.transform(X_train[int_vars])\nX_test[int_vars] = sc.transform(X_test[int_vars])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(X_train.head())\ndisplay(X_test.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Train the model\n\nparameters = {\n    'max_depth': [5, 10],\n    'n_estimators': [500, 1000],\n    'learning_rate': [0.01, 0.005]\n}\n\n#Cross Validation\ngrid_search = GridSearchCV(\n    estimator=xgboost.XGBRegressor(random_state = i),\n    param_grid=parameters,\n    n_jobs = 5,\n    cv = 5\n)\n\ngrid_search.fit(X_train, y_train)\nprint(\"Best Parms\", grid_search.best_params_)\nprint(\"Best Score: \", grid_search.best_score_)\nbest_model = grid_search.best_estimator_    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next Steps/Future Work:\n- Compute the **Confidence** in each prediction\n- Create **submission.csv** file\n- Possibly, incorporate the CT data in the model\n    ","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}